# Randomly weighted CNNs for (music) audio classification

###### Abstract

The computer vision literature shows that randomly weighted neural networks
perform reasonably as feature extractors. Following this idea, we study how
non-trained (randomly weighted) convolutional neural networks perform as
feature extractors for (music) audio classification tasks. We use features
extracted from the embeddings of deep architectures as input to a classifier –
with the goal to compare classification accuracies when using different
randomly weighted architectures. By following this methodology, we run a
comprehensive evaluation of the current deep architectures for audio
classification, and provide evidence that the architectures alone are an
important piece for resolving (music) audio problems using deep neural
networks.

##  1 MOTIVATION – FROM PREVIOUS WORKS

Some intriguing properties of deep neural networks are periodically showing up
in the scientific literature. Examples of those are: (i) perceptually non-
relevant signal perturbations that dramatically affect the predictions of an
image classifier [12, 49]; (ii) although there is no guarantee of converging
to a global minima that might generalize, image classification models perform
well with unseen data [25, 14]; or (iii) non-trained deep neural networks are
able to perform reasonably well as image feature extractors [43, 51, 41]. In
this work, we exploit one of the above listed properties (iii) to evaluate how
discriminative deep audio architectures are before training.

Figure 1: CNN front-ends for audio classification tasks – with examples of
possible configurations for every paradigm.

Previous works already explored the idea of empirically studying the qualities
of non-trained (randomly weighted) networks, but mainly in the computer vision
field:

$\mathbf{\cdot}$ Saxe et al. [43] studied how discriminative are the
architectures themselves by evaluating the classification performance of SVMs
fed with features extracted from non-trained (random) CNNs.111CNNs stands for
Convolutional Neural Networks.They showed that a surprising fraction of the
performance in deep image classifiers can be attributed to the architecture
alone. Therefore, the key to good performance lies not only on improving the
learning algorithms but also in searching for the most suitable architectures.
Further, they showed that the (classification) performance delivered by random
CNN features is correlated with the results of their end-to-end trained
counterparts – this result means, in practice, that one can bypass the time-
consuming process of learning for evaluating a given architecture. We build on
top of this result to evaluate current CNN architectures for audio
classification.

$\mathbf{\cdot}$ Rosenfeld and Tsotsos [41] fixed most of the model’s weights
to be random, and only allowed a small portion of them to be learned. By
following this methodology, they showed a small decrease in image
classification performance when these models were compared to their fully
trained counterparts. Further, the performance of their non fully-trained
models can be summarized as follows:

DenseNet[17] $\gg$ ResNet[15] $>$ VGG[48] $\gg$ AlexNet[25]

What matches previous works reporting how these (fully trained) models
perform[17, 15, 48], confirming the performance correlation between randomly
weighted models and their trained counterparts found by Saxe et al. [43]

$\mathbf{\cdot}$ Adebayo et al. [1] empirically assessed the local
explanations of deep image classifiers to find that randomly weighted models
produce explanations similar to those produced by models with learned weights.
They conclude that the architectures introduce a strong prior which affects
the learned (and not learned) representations.

$\mathbf{\cdot}$ Ulyanov et al. [51] also showed that the structure of a
network (the non-trained architecture) is sufficient to capture useful
features for the tasks of image denoising, super-resolution and inpainting.
They think of any designed architecture as a hand-crafted model where prior
information is embedded in the structure of the network itself. This way of
thinking resonates with the rationale behind the family of audio models
designed considering domain knowledge (see section 2) – what denotes that in
both audio and image fields it exists the interest of bringing together the
end-to-end learning literature and previous research.

Few related works exist in the audio field – and every randomly weighted
neural network we found in the audio literature was a mere baseline [7, 24,
2]. Inspired by previous computer vision works, we study which audio
architectures work the best via evaluating how non-trained CNNs perform as
feature extractors. To this end, we use the CNNs’ embeddings to construct
feature vectors for a classifier – with the goal to compare classification
performances when different randomly weighted architectures are used to
extract features. To the best of our knowledge, this is the first
comprehensive evaluation of randomly weighted CNNs for (music) audio
classification.

Extreme learning machines (ELMs) [47, 32, 18] and echo state networks (ESNs)
[19] are also closely related to our work. In short, ELMs are
classification/regression models222Support Vector Machines are also
classification/regression models.that are based on a single-layer feed-forward
neural network with random weights. They work as follows: first, ELMs randomly
project the input into a latent space; and then, learn how to predict the
output via a least-square fit. More formally, we aim to predict:

$\hat{Y}=W_{2}\hskip 2.84526pt\sigma(W_{1}X),$

where $W_{1}$ is the (randomly weighted) matrix of input-to-hidden-layer
weights, $\sigma$ is the non-linearity, $W_{2}$ is the matrix of hidden-to-
output-layer weights, and X represents the input. The training algorithm is as
follows: 1) set $W_{1}$ with random values; 2) estimate $W_{2}$ via a least-
squares fit:

$W_{2}=\sigma(W_{1}X)^{+}Y$

where + denotes the Moore-Penrose inverse. Since no iterative process is
required for learning the weights, training is faster than stochastic gradient
descent [18]. Provided that we process audio signals with randomly weighted
CNNs, ELM-based classifiers are a natural choice for our study – so that all
the pipeline (except the last layer) is based on random projections that are
only constrained by the structure of the neural network. Although ELMs are not
widely used by the audio community, they have been used for speech emotion
recognition[13, 21], or for music audio classification[44, 29, 23]. ESNs
differ from ELMs in that their random projections use recurrent connections.
Given that the audio models we aim to study are not recurrent, we leave for
future work using ESNs – see [45, 16] for audio applications of ESNs.

##  2 ARCHITECTURES

In this work we evaluate the most used deep learning architectures for (music)
audio classification. In order to facilitate the discussion around these
architectures, we divide the deep learning pipeline into two parts: front-end
and back-end, see Figure 2. The front-end is the part that interacts with the
input signal in order to map it into a latent-space, and the back-end predicts
the output given the representation obtained by the front-end. Note that one
can interpret the front-end as a “feature extractor” and the back-end as a
“classifier”. Given that we compare how several non-trained (random) CNNs
perform as feature extractors, and we will use out-of-the-box classifiers to
predict the classes: this literature review focuses in introducing the main
deep learning front-ends for audio classification.

Figure 2: The deep learning pipeline.

Front-ends — These are generally conformed by CNNs [9, 6, 53, 38, 39], since
these can encode efficient representations by sharing weights333Which
constitute the (eventually learnt) feature representations.along the signal.
Figure 1 depicts six different CNN front-end paradigms, which can be divided
into two groups depending on the used input signal: waveforms [9, 53, 27] or
spectrograms [6, 38, 39]. Further, the design of the filters can be either
based on domain knowledge or not. For example, one leverages domain knowledge
when the frame-level single-shape444Italicized names correspond to the front-
end types in Figure 1. front-end for waveforms is designed so that the length
of the filter is set to be the same as the window length in a STFT [9]. Or for
a spectrogram front-end, it is used vertical filters to learn spectral
representations [26] or horizontal filters to learn longer temporal cues [46].
Generally, a single filter shape is used in the first CNN layer [9, 6, 26,
46], but some recent work reported performance gains when using several filter
shapes in the first layer [53, 38, 39, 34, 36, 5]. Using many filters promotes
a more rich feature extraction in the first layer, and facilitates leveraging
domain knowledge for designing the filters’ shape. For example: a frame-level
many-shapes front-end for waveforms can be motivated from a multi-resolution
time-frequency transform555The Constant-Q Transform [3] is an example of such
transform.perspective – with window sizes varying inversely with frequency
[53]; or since it is known that some patterns in spectrograms are occurring at
different time-frequency scales, one can intuitively incorporate many vertical
and/or horizontal filters to efficiently capture those in a spectrogram front-
end [38, 39, 36, 34]. As seen, using domain knowledge when designing the
models allows to naturally connect the deep learning literature with previous
relevant signal processing work. On the other hand, when domain knowledge is
not used, it is common to employ a deep stack of small filters, e.g.:
3$\times$1 in the sample-level front-end used for waveforms [27, 52, 40], or
3$\times$3 in the small-rectangular filters front-end used for spectrograms
[6]. These VGG-like666VGG: a computer vision model based on a deep stack of
3$\times$3 filters.models make minimal assumptions over the local
stationarities of the signal, so that any structure can be learnt via
hierarchically combining small-context representations.

##  3 METHOD

Our goal is to study which CNN front-ends work best via evaluating how non-
trained models perform as feature extractors. Our evaluation is based on the
traditional pipeline of features extraction \+ classifier. We use the
embeddings of non-trained (random) CNNs as features: for every audio clip, we
compute the average of each feature map (in every layer) and concatenate these
values to construct a feature vector[7]. The baseline feature vector is
constructed from 20 MFCCs, their $\Delta$s and $\Delta\Delta$s. We compute
their mean and standard deviation through time, and the resulting feature
vector is of size 120. We set the widely used MFCCs + SVM setup as baseline.
To allow a fair comparison with the baseline, CNN models have $\approx$ 120
feature maps – so that the resulting feature vectors have a similar size as
the MFCCs vector. Further, we evaluate an alternative configuration with more
feature maps ($\approx$3500) to show the potential of this approach. Model’s
description omit the number of filters per layer for simplicity – full
implementation details are accessible online.88footnotemark: 8

###  3.1 Features: randomly weighted CNNs

Except for the VGG model that uses ELUs as non-linearities [6, 8], the rest
use ReLUs [10] – and we do not use batch normalization, see discussion in
section 3.4. We use waveforms and spectrograms as input to our CNNs:

Waveform inputs — are of $\approx$ 29sec (350,000 samples at 12kHz) and the
following architectures are under study:

$\cdot$ Sample-level: is based on a stack of 7 blocks that are composed by a
1D-CNN layer (filter length: 3, stride: 1) and a max-pool layer (size: 3,
stride: 3) – with the exception of the input block which has no max-pooling
and its 1D-CNN layer has a stride of 3 [27]. Averages to construct the feature
vector are computed after every pooling layer, except for the first layer that
are computed after the CNN.

$\cdot$ Frame-level many-shapes: consists of a 1D-CNN layer with five filter
lengths: 512, 256, 128, 64, 32 [53]. Every filter’s stride is of 32 and we use
same padding – to easily concatenate feature maps of the same size. Note that
out of this single 1D-CNN layer, five feature maps (resulting of the different
filter length convolutions) are concatenated. For that reason, every feature
map needs to have the same (temporal) size. Since this model is single-layered
and it might be in clear disadvantage when compared to the sample-level CNN,
we increase its depth via adding three more 1D-CNN layers (length: 7, stride:
1) – where the last two layers have residual connections, and the penultimate
layer’s feature map is down-sampled by two (MP x2), see Figure 3. Averages to
construct the feature vector are calculated for each feature map after every
1D-CNN layer.

$\cdot$ Frame-level: consists of a 1D-CNN layer with a filter of length 512
[9]. Stride is set to be 32 to allow a fair comparison with the frame-level
many-shapes architecture. As in frame-level many-shapes, we increase the depth
of the model via adding three more 1D-CNN layers – as in Figure 3. Averages to
construct the feature vector are calculated for each feature map after every
1D-CNN layer.

Figure 3: Additional layers for the frame-level & frame-level many-shapes
architectures, similar as in [9, 37] – where MP stands for max pooling.

Spectrogram inputs — are set to be log-mel spectrograms(spectrograms size:
1376$\times$96777STFT parameters: window_size = 512, hop_size=256, and
fs=12kHz., being $\approx$ 29sec of signal).Differently from waveform models,
spectrogram architectures use no additional layers to deepen single-layered
CNNs because these already deliver a reasonable classification performance.
Unless we state the contrary, every CNN layer used for processing spectrograms
is set to have stride 1. As for the frame-level many-shapes model, we use same
padding when many filter shapes are used in the same layer. The following
spectrogram models are studied:

$\cdot$ 7$\times$96: consists of a single 1D-CNN layer with filters of length
7 that convolve through the time axis [9]. As a result: CNN filters are
vertical and of shape 7$\times$96\. Therefore, these filters can encode
spectral (timbral) representations. Averages to construct the feature vector
are calculated for each feature map after the 1D-CNN layer.

$\cdot$ 7$\times$86: consists of a single 2D-CNN layer with vertical filters
of shape 7$\times$86 [39, 36]. Since its vertical shape is smaller than the
input (86$<$96), filters can also convolve through the frequency axis – what
can be seen as “pitch shifting” the filter. Consequently, 7$\times$86 filters
can encode pitch-invariant timbral representations [39, 36]. Further, since
the resulting activations can carry pitch-related information, we max-pool the
frequency axis to get pitch-invariant features (max-pool shape: 1$\times$11).
Averages to construct the feature vector are calculated for each feature map
after the max-pool layer.

$\cdot$ Timbral: consists of a single 2D-CNN layer with many vertical filters
of shapes: 7$\times$86, 3$\times$86, 1$\times$86, 7$\times$38, 3$\times$38,
1$\times$38, see Figure 4 (top)[35, 11, 39]. These filters can also convolve
through the frequency axis and therefore, these can encode pitch-invariant
representations. Several filter shapes are used to efficiently capture
different timbrically relevant time-frequency patterns, e.g.: kick-drums (can
be captured with 7$\times$38 filters representing sub-band information for a
short period of time), or string ensemble instruments (can be captured with
1$\times$86 filters representing spectral patterns spread in the frequency
axis). Further, since the resulting activations can carry pitch-related
information, we max-pool the frequency axis to get pitch-invariant features
(max-pool shapes: 1$\times$11 or 1$\times$59). Averages to construct the
feature vector are calculated for each feature map after the max-pool layer.

$\cdot$ Temporal: several 1D-CNN filters (of lengths: 165, 128, 64, 32)
operate over an energy envelope obtained via mean-pooling the frequency-axis
of the spectrogram, see Figure 4 (bottom). By computing the energy envelope in
that way, we are considering high and low frequencies together while
minimizing the computations of the model. Observe that this single-layered
1D-CNN is not operating over a 2D spectrogram, but over a 1D energy envelope –
therefore no vertical convolutions are performed, only 1D (temporal)
convolutions are computed. As seen, domain knowledge can also provide guidance
to (effectively) minimize the computations of the model. Averages to construct
the feature vector are calculated for each feature map after the CNN layer.

$\cdot$ Timbral+temporal: combines both timbral and temporal CNNs in a single
(but wide) layer, see Figure 4[37]. Averages to construct the feature vector
are calculated in the same way as for timbral and temporal architectures.

$\cdot$ VGG: is a computer vision model based on a stack of 5 blocks combining
2D-CNN layers (with small rectangular filters of 3$\times$3) and max-pooling
(of shapes: 4$\times$2, 4$\times$3, 5$\times$2, 4$\times$2, 4$\times$4,
respectively)[6]. Averages to construct the feature vector are computed after
every pooling layer.

Figure 4: Timbral+temporal architecture. MP: max-pool.

As seen, studied architectures are representative of the audio classification
state-of-the-art – introduced in section 2. For further details about the
models under study: the code is accessible
online888https://github.com/jordipons/elmarc, and a graphical
conceptualization of the models is available in Figures 1, 3 and 4.

###  3.2 Classifiers: SVM and ELM

We study how several feature vectors (computed considering different CNNs)
perform for a given set of classifiers: SVMs and ELMs. We discarded the use of
other classifiers since their performance was not competitive when compared to
those. SVMs and ELMs are hyper-parameter sensitive, for that reason we perform
a grid search:

$\cdot$ SVM hyper-parameters: we consider both linear and rbf kernels. For the
rbf kernel, we set $\gamma$ to: $2^{-3}$, $2^{-5}$, $2^{-7}$, $2^{-9}$,
$2^{-11}$, $2^{-13}$, #features-1; and for every kernel configuration, we try
several C’s (penalty parameter): 0.1, 2, 8, 32. We use scikit-learn’s SVM
implementation [33].

$\cdot$ ELM’s main hyper-parameter is the number of hidden units: 100, 250,
500, 1200, 1800, 2500. We use ReLUs as non-linearity, and we use a public ELM
implementation.999https://github.com/zygmuntz/Python-ELM

###  3.3 Datasets: music and acoustic events

$\cdot$ GTZAN fault-filtered version [50, 22]. Training songs: 443, validation
songs: 197, and test songs: 290 – divided in 10 classes. We use this dataset
to study how randomly weighted CNNs perform for music genre classification.

$\cdot$ Extended Ballroom [30, 4] – 4,180 songs divided in 13 classes; 10
stratified folds are randomly generated for cross-validation. We use this
dataset to study how randomly weighted CNNs classify rhythm/tempo classes.

$\cdot$ Urban Sound 8K [42] – 8732 acoustic events divided in 10 classes; 10
folds are already defined by the dataset authors for cross-validation. Since
urban sounds are shorter than 4 seconds and our models accepts $\approx$ 29sec
inputs, the signal is repeated to create inputs of the same length. We use
this dataset to study how randomly weighted CNNs perform to classify natural
(non-music) sounds.

###  3.4 Reproducing former results to discuss our method

Choi et al. [7] used random CNN features as baseline for their work, and found
that (in most cases) these random CNN features perform better than MFCCs.
Motivated by this result, we pursue this idea for studying how different deep
architectures perform when resolving audio problems. To start, we first
reproduce one of their experiments using random CNNs – under the same
conditions101010https://github.com/keunwoochoi/transfer_learning_music/ (more
information is available in issue #2): the GTZAN dataset is split in 10
stratified folds used for cross-validation111111Our work does not use this
split, we use the fault-filtered version., a VGG architecture with batch
normalization is employed, and the classifier is an SVM. We found that results
can vary depending on the batch size if, when computing the feature vectors,
layers are normalized with the statistics of current batch (batch
normalization). For example: if audio-features of the same genre are batch-
normalized by the same factor, one can create an artificial genre cue that
might affect the results. One can observe this phenomena in Figure 5, where
the best results are achieved when all songs of the same genre fill a full
batch (batch size of 100).121212The GTZAN has 10 genres with 100 audios each,
one can fill batches of 100 with audios of the same genre via sorting the data
by genres.We also observe that small batch sizes are harming the model’s
performance – see in Figure 5 when batch sizes are set to 1 and 10. And
finally, when batch normalization is not used, no matter what batch size we
use that the results remain the same – ANOVA test with $H_{0}$ being that
averages are equal(p-value=0.491). Since it is not desirable that performance
depends on the batch size, and that the feature vector of an audio depends on
other audios (that are present in the batch): we decided not to use batch
normalization.

Figure 5: Random CNN features behavior when using (or not) batch
normalization. Dataset: GTZAN, 10-fold cross-validation. Performance metric
(%): average accuracies (and standard deviations) across 3 runs. Classifier:
SVM.

##  4 Results

Figures show average accuracies across 3 runs for every feature type (listed
on the right with the length of the feature vector in parenthesis). We use a
t-test to reveal which models are performing the best – $H_{0}$: averages are
equal.

###  4.1 GTZAN: music genre recognition

Figure 6: Accuracy (%) results for the GTZAN dataset with random CNN feature
vectors of length $\approx$ 120. Figure 7: Accuracy (%) results for the GTZAN
dataset with random CNN feature vectors of length $\approx$ 3500.

The sample-level waveform model always performs better than frame-level many-
shapes (t-test: p-value$\ll$0.05). The two best performing spectrogram-based
models are: timbral+temporal and VGG – with a remarkable performance of the
timbral model alone. The timbral+temporal CNN performs better than VGG when
using the ELM ($\approx$3500) classifier (t-test: p-value=0.017); but in other
cases, both models perform equivalently (t-test: p-value$>$0.05). Moreover,
the 7x86 model performs better than 7x96 when using SVMs (t-test:
p-value$<$0.05), but when using ELMs: 7x96 and 7x86 perform equivalently
(t-test: p-value$\gg$0.05). The best VGG and timbral+temporal models achieve
the following (average) accuracies: 59.65% and 56.89%, respectively – both
with an SVM ($\approx$3500) classifier. Both models outperform the MFCCs
baseline: 53.44% (t-test: p-value$<$0.05), but these random CNNs perform worse
than a CNN pre-trained with the Million Song Dataset: 82.1% [28]. Finally,
note that although timbral and timbral+temporal models are single-layered,
these are able to achieve remarkable performances – showing that single-
layered spectrogram front-ends (attending to musically relevant contexts) can
do a reasonable job without paying the cost of going deep [39, 36]. Thus
meaning, e.g., that the saved capacity can now be employed by the back-end to
learn (some more) representations.

###  4.2 Extended Ballroom: rhythm/tempo classification

Figure 8: Accuracy (%) results for the Extended Ballroom dataset with random
CNN feature vectors of length $\approx$ 120. Figure 9: Accuracy (%) results
for the Extended Ballroom dataset with random CNN feature vectors of length
$\approx$ 3500

The sample-level waveform model always performs better than frame-level many-
shapes (t-test: p-value$\ll$0.05). The two best performing spectrogram-based
models are: temporal and timbral+temporal, but the temporal model performs
better than timbral+temporal in all cases (t-test: p-value$\ll$0.05) –
denoting that spectral cues can be a confounding factor for this dataset.
Moreover, the 7x86 model performs better than 7x96 in all cases (t-test:
p-value$<$0.05). The best (average) accuracy score is obtained using temporal
models and SVMs ($\approx$3500): 89.82%. Note that the temporal model clearly
outperforms the MFCCs baseline: 63.25%(t-test: p-value$\ll$0.05) and,
interestingly, it performs slightly worse than a trained CNN: 93.7% [20]. This
result confirms that the architectures (alone) introduce a strong prior which
can significantly affect the performance of an audio model. Thus meaning that,
besides learning, designing effective architectures might be key for resolving
(music) audio tasks with deep learning. In line with that, note that the
temporal architecture is designed considering musical domain knowledge – in
this case: how tempo & rhythm are expressed in spectrograms. Hence, its good
performance also validates the design strategy of using musically motivated
architectures as a way to intuitively navigate through the network parameters
space [36, 39, 38].

###  4.3 Urban Sounds 8K: acoustic event detection

For these experiments we do not use the temporal model (with 1D-CNNs of length
165, 128, 64, 32). Instead, we study the temporal+time model – where time
follows the same design as temporal but with filters of length: 64, 32, 16, 8.
This change is motivated by the fact that temporal cues in (natural) sounds
are shorter and less important than temporal cues in music (i.e., rhythm or
tempo).

Figure 10: Accuracy (%) results for the Urban Sounds 8k dataset with random
CNN feature vectors of length $\approx$ 120. Figure 11: Accuracy (%) results
for the Urban Sounds 8k dataset with random CNN feature vectors of length
$\approx$ 3500

The sample-level waveform model always performs better than frame-level many-
shapes (t-test: p-value$\ll$0.05). The two best performing spectrogram-based
models are: VGG and timbral+time – but VGG performs better than timbral+time
in all cases (t-test: p-value$\ll$0.05). Also, the 7x86 model performs better
than 7x96 in all cases (t-test: p-value$<$0.075). The best (average) accuracy
score is obtained using VGG and SVMs ($\approx$3500): 70.74% – outperforming
the MFCCs baseline: 65.49%(t-test: p-value$<$0.05), and performing slightly
worse than a trained CNN: 73%131313The same CNN achieves 79% when trained with
data augmentation.[28]. Finally, note that VGGs achieved remarkable results
when recognizing genres and detecting acoustic events – tasks where timbre is
an important cue. As a result: one could argue that VGGs are good at
representing spectral features. Hence, these might be of utility for tasks
where spectral cues are relevant.

##  5 CONCLUSIONS

This study builds on top of prior works showing that the (classification)
performance delivered by random CNN features is correlated with the results of
their end-to-end trained counterparts [43, 41]. We use this property to run a
comprehensive evaluation of current deep architectures for (music) audio. Our
method is as follows: first, we extract a feature vector from the embeddings
of a randomly weighted CNN; and then, we input these features to a classifier
– which can be an SVM or an ELM. Our goal is to compare the obtained
classification accuracies when using different CNN architectures. The results
we obtain are far from random, since: (i) randomly weighted CNNs are (in some
cases) close to match the accuracies obtained by trained CNNs; and (ii) these
are able to outperform MFCCs. This result denotes that the architectures alone
are an important piece of the deep learning solution and therefore, searching
for efficient architectures capable to encode the specificities of (music)
audio signals might help advancing the state of our field. In line with that,
we have shown that (musical) priors embedded in the structure of the model can
facilitate capturing useful (temporal) cues for classifying rhythm/tempo
classes. Besides, we show that for waveform front-ends: sample-level $\gg$
frame-level many-shapes $>$ frame-level – as noted in the (trained) literature
[27, 53, 52]. The differential aspect of the sample-level front-end is that
its representational power is constructed via hierarchically combining small-
context representations, not by exploiting prior knowledge about waveforms.
Further, we show that for spectrogram front-ends:7x96$<$7x86 – as shown in
prior (trained) works [36, 31]. By allowing the filters to convolve through
the frequency axis, the architecture itself facilitates capturing pitch-
invariant timbral representations. Finally: timbral (+temporal/time) and VGG
spectrogram front-ends achieve remarkable results for tasks where timbre is
important – as previously noted in the (trained) literature [39]. Their
respective advantages being that: (i) timbral (+temporal/time) architectures
are single-layered front-ends which explicitly capture acoustically relevant
receptive fields – which can be known via exploiting prior knowledge about the
task; and (ii) VGG front-ends require no prior domain knowledge about the task
for its design. Although our main conclusions are backed by additional results
in the (trained) literature, we leave for future work consolidating those via
doing a similar study considering trained models.

##  6 ACKNOWLEDGMENTS

This work was partially supported by the Maria de Maeztu Units of Excellence
Programme (MDM-2015-0502) – and we are grateful for the GPUs donated by
NVidia.

## References

  * [1] Julius Adebayo, Justin Gilmer, Ian Goodfellow, and Been Kim.  Local explanation methods for deep neural networks lack sensitivity to parameter values.  In International Conference on Learning Representations Workshop (ICLR), 2018.
  * [2] Relja Arandjelovic and Andrew Zisserman.  Look, listen and learn.  In IEEE International Conference on Computer Vision (ICCV), pages 609–617. IEEE, 2017.
  * [3] Judith C Brown.  Calculation of a constant q spectral transform.  The Journal of the Acoustical Society of America, 89(1):425–434, 1991.
  * [4] Pedro Cano, Emilia Gómez, Fabien Gouyon, Perfecto Herrera, Markus Koppenberger, Beesuan Ong, Xavier Serra, Sebastian Streich, and Nicolas Wack.  ISMIR 2004 audio description contest.  Music Technology Group of the Universitat Pompeu Fabra, Technical Report, 2006.
  * [5] Ning Chen and Shijun Wang.  High-level music descriptor extraction algorithm based on combination of multi-channel cnns and lstm.  In International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * [6] Keunwoo Choi, George Fazekas, and Mark Sandler.  Automatic tagging using deep convolutional neural networks.  International Society for Music Information Retrieval Conference (ISMIR), 2016.
  * [7] Keunwoo Choi, György Fazekas, Mark Sandler, and Kyunghyun Cho.  Transfer learning for music classification and regression tasks.  International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * [8] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289, 2015.
  * [9] Sander Dieleman and Benjamin Schrauwen.  End-to-end learning for music audio.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014.
  * [10] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In International Conference on Artificial Intelligence and Statistics, pages 315–323, 2011.
  * [11] Rong Gong, Jordi Pons, and Xavier Serra.  Audio to score matching by combining phonetic and duration information.  International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * [12] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572, 2014.
  * [13] Kun Han, Dong Yu, and Ivan Tashev.  Speech emotion recognition using deep neural network and extreme learning machine.  In Annual Conference of the International Speech Communication Association, 2014.
  * [14] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In IEEE International Conference on Computer Vision (ICCV), pages 1026–1034, 2015.
  * [15] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In IEEE conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [16] Georg Holzmann.  Reservoir computing: A powerful black-box framework for nonlinear audio processing.  In International Conference on Digital Audio Effect (DAFx), 2009\.
  * [17] Gao Huang, Zhuang Liu, Kilian Q Weinberger, and Laurens van der Maaten.  Densely connected convolutional networks.  In IEEE conference on Computer Vision and Pattern Recognition (CVPR), volume 1, page 3, 2017.
  * [18] Guang-Bin Huang, Qin-Yu Zhu, and Chee-Kheong Siew.  Extreme learning machine: theory and applications.  Neurocomputing, 70(1-3):489–501, 2006.
  * [19] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  German National Research Center for Information Technology, Technical Report, 148(34):13, 2001.
  * [20] Yeonwoo Jeong, Keunwoo Choi, and Hosan Jeong.  Dlr: Toward a deep learned rhythmic representation for music content analysis.  arXiv preprint arXiv:1712.05119, 2017.
  * [21] Heysem Kaya and Albert Ali Salah.  Combining modality-specific extreme learning machines for emotion recognition in the wild.  Journal on Multimodal User Interfaces, 10(2):139–149, 2016.
  * [22] Corey Kereliuk, Bob L Sturm, and Jan Larsen.  Deep learning and music adversaries.  IEEE Transactions on Multimedia, 17(11):2059–2071, 2015.
  * [23] Suisin Khoo, Zhihong Man, and Zhenwei Cao.  Automatic han chinese folk song classification using extreme learning machines.  In Australasian Joint Conference on Artificial Intelligence, pages 49–60. Springer, 2012.
  * [24] Jaehun Kim, Julián Urbano, Cynthia Liem, and Alan Hanjalic.  One deep music representation to rule them all?: A comparative analysis of different representation learning strategies.  arXiv preprint arXiv:1802.04051, 2018.
  * [25] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in Neural Information Processing Systems (NIPS), pages 1097–1105, 2012.
  * [26] Honglak Lee, Peter Pham, Yan Largman, and Andrew Y Ng.  Unsupervised feature learning for audio classification using convolutional deep belief networks.  In Advances in Neural Information Processing Systems (NIPS), 2009\.
  * [27] Jongpil Lee, Jiyoung Park, Keunhyoung Luke Kim, and Juhan Nam.  Sample-level deep convolutional neural networks for music auto-tagging using raw waveforms.  Sound and Music Computing Conference (SMC), 2017.
  * [28] Jongpil Lee, Jiyoung Park, Keunhyoung Luke Kim, and Juhan Nam.  Samplecnn: End-to-end deep convolutional neural networks using very small filters for music classification.  Applied Sciences, 8(1):150, 2018.
  * [29] Qi-Jun Benedict Loh and Sabu Emmanuel.  Elm for the classification of music genres.  In International Conference on Control, Automation, Robotics and Vision, pages 1–6. IEEE, 2006.
  * [30] Ugo Marchand and Geoffroy Peeters.  Scale and shift invariant time/frequency representation using auditory statistics: Application to rhythm description.  In International Workshop on Machine Learning for Signal Processing, pages 1–6. IEEE, 2016.
  * [31] Sergio Oramas, Oriol Nieto, Francesco Barbieri, and Xavier Serra.  Multi-label music genre classification from audio, text, and images using deep features.  International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * [32] Yoh-Han Pao, Gwang-Hoon Park, and Dejan J Sobajic.  Learning and generalization characteristics of the random vector functional-link net.  Neurocomputing, 6(2):163–180, 1994.
  * [33] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [34] Huy Phan, Lars Hertel, Marco Maass, and Alfred Mertins.  Robust audio event recognition with 1-max pooling convolutional neural networks.  arXiv preprint arXiv:1604.06338, 2016.
  * [35] Jordi Pons, Rong Gong, and Xavier Serra.  Score-informed syllable segmentation for a cappella singing voice with convolutional neural networks.  International Society for Music Information Retrieval Conference (ISMIR), 2017.
  * [36] Jordi Pons, Thomas Lidy, and Xavier Serra.  Experimenting with musically motivated convolutional neural networks.  In International Workshop on Content-Based Multimedia Indexing (CBMI), pages 1–6. IEEE, 2016.
  * [37] Jordi Pons, Oriol Nieto, Matthew Prockup, Erik M Schmidt, Andreas F Ehmann, and Xavier Serra.  End-to-end learning for music audio tagging at scale.  Workshop on Machine Learning for Audio Signal Processing (ML4Audio) at NIPS, 2017.
  * [38] Jordi Pons and Xavier Serra.  Designing efficient architectures for modeling temporal features with convolutional neural networks.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017.
  * [39] Jordi Pons, Olga Slizovskaia, Rong Gong, Emilia Gómez, and Xavier Serra.  Timbre analysis of music audio signals with convolutional neural networks.  European Signal Processing Conference (EUSIPCO), 2017.
  * [40] Dario Rethage, Jordi Pons, and Xavier Serra.  A Wavenet for speech denoising.  IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2018.
  * [41] Amir Rosenfeld and John K Tsotsos.  Intriguing properties of randomly weighted networks: Generalizing while learning next to nothing.  arXiv preprint arXiv:1802.00844, 2018.
  * [42] J. Salamon, C. Jacoby, and J. P. Bello.  A dataset and taxonomy for urban sound research.  In ACM International Conference on Multimedia, Orlando, FL, USA, Nov. 2014.
  * [43] Andrew M Saxe, Pang Wei Koh, Zhenghao Chen, Maneesh Bhand, Bipin Suresh, and Andrew Y Ng.  On random weights and unsupervised feature learning.  In International Conference on Machine Learning (ICML), pages 1089–1096, 2011.
  * [44] Simone Scardapane, Danilo Comminiello, Michele Scarpiniti, and Aurelio Uncini.  Music classification using extreme learning machines.  In International Symposium on Image and Signal Processing and Analysis (ISPA), pages 377–381. IEEE, 2013.
  * [45] Simone Scardapane and Aurelio Uncini.  Semi-supervised echo state networks for audio classification.  Cognitive Computation, 9(1):125–135, 2017.
  * [46] Jan Schluter and Sebastian Bock.  Improved musical onset detection with convolutional neural networks.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014.
  * [47] Wouter F Schmidt, Martin A Kraaijveld, and Robert PW Duin.  Feedforward neural networks with random weights.  In International Conference on Pattern Recognition, pages 1–4. IEEE, 1992.
  * [48] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [49] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  arXiv preprint arXiv:1312.6199, 2013.
  * [50] George Tzanetakis and Perry Cook.  Musical genre classification of audio signals.  IEEE Transactions on speech and audio processing, 10(5):293–302, 2002.
  * [51] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep image prior.  arXiv preprint arXiv:1711.10925, 2017.
  * [52] Aaron Van Den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  Wavenet: A generative model for raw audio.  arXiv preprint arXiv:1609.03499, 2016.
  * [53] Zhenyao Zhu, Jesse H Engel, and Awni Hannun.  Learning multiscale features directly from waveforms.  arXiv preprint arXiv:1603.09509, 2016.


# Quality Control in Particle Precipitation via Robust Optimization

Martina Kuchlbauer  Jana Dienstbier  Adeel Muneer  Hanna Hedges  Michael
Stingl  Frauke Liers  Lukas Pflug  Friedrich-Alexander-Universität Erlangen-
Nürnberg (FAU), Optimization under Uncertainty and Data Analysis, Departement
of Data Science, Cauerstraße 11, 91058 Erlangen, Germany  Friedrich-Alexander-
Universität Erlangen-Nürnberg (FAU), Chair of Applied Mathematics (Continuous
Optimization), Cauerstraße 11, 91058 Erlangen, Germany  Friedrich-Alexander-
Universität Erlangen-Nürnberg (FAU), Central Institute for Scientific
Computing, Martensstr. 5a, 91058 Erlangen, Germany

###### Abstract

We propose a robust optimization approach to mitigate the impact of
uncertainties in particle precipitation. Our model of particle synthesis
incorporates, as partial differential equations, nonlinear and nonlocal
population balance equations. The optimization goal is to design products with
desired size distributions. Recognizing the impact of uncertainties, we extend
the model to robustly hedge against them to ensure tailored particle sizes.
For the robust problem, we enhance an adaptive bundle framework for nonlinear
robust optimization integrating the exact method of moments approach for the
population balance equations. Computational experiments focus on uncertainties
in precursor inflow rates, which greatly influence the resulting product’s
quality. Using realistic parameter values for quantum dot synthesis, we
demonstrate the algorithm’s efficiency and find that unprotected processes
fail to achieve desired particle sizes, even for small uncertainties,
highlighting the need for robust processes. The latter outperforms the
unprotected process concerning the product’s quality in perturbed scenarios.

###### keywords:

particle design , nonlinear robust optimization , process optimization ,
population balance equation , exact method of moments , adaptive bundle
framework

##  1 Introduction

Designing (nano)particles with optimal properties involves numerous
challenging tasks, such as establishing suitable experimental setups and
processes to produce the desired products. Indeed, the optimization of
particle size and shape has a wide variety of applications, ranging from
development of cosmetics and pharmaceutical usage for oral delivery of peptide
drugs to the performance improvement of lithium batteries [21].
Achieving optimal results in particle synthesis is particularly challenging
due to the presence of uncertainties in, e.g., material parameters,
temperature, particle nucleation, and growth, to name just a few, that can
have a detrimental impact on the quality of the final product. Since even
small disturbances can lead to a deterioration of product quality [19], it is
essential to ensure that possible uncertainties or fluctuations affect the
product as little as possible. Neglecting uncertainties in particle design can
have severe consequences, as even slight variations in their dispersity can
significantly impact the performance of the final product.

As a consequence, robust mathematical optimization offers an approach to
overcome this issue by enabling the development of optimized processes that
are best possible, even under uncertain conditions. Developing optimization
frameworks for particle synthesis and employing suitable methods is of utmost
importance in achieving efficient and practical manufacturing processes.
Robust design is particularly crucial in industrial contexts to guarantee
desired characteristics such as absorption, emission properties, quantum
yield, and maximum mass while optimizing production costs. Uncertainties,
encompassing factors like solubility, flow rates, and temperature, profoundly
influence nanoparticle synthesis even at the laboratory level. These
uncertainties, coupled with challenges in characterization, significantly
affect the overall quality of the product. Consequently, effectively managing
and mitigating uncertainties and characterization issues becomes imperative to
attain nanoparticles of high quality. For hedging against uncertainty, the
currently very active research area of robust optimization is very well
suited. However, robust optimization adds additional challenges to the
optimization of process designs. In this work, we contribute to this area by
introducing a practically efficient robust solution algorithm that is able to
determine optimal processes for particle precipitation. We show its
effectiveness for robust processes.

### Our Contribution

We present a quality control system for particle precipitation modeled by
population balance equations (PBEs). For this nonlinear, nonlocal class of
partial differential equations (PDEs), no analytical closed-form solutions are
available. However, semi-analytical solutions can be used by the exact method
of moment (eMoM) approach. For hedging against uncertainties, we integrate
eMoM into a novel adaptive bundle framework for robust (nonlinear)
optimization. We use it in particular for optimizing particle size
distributions (PSDs) under uncertain precursor inflow rates as the latter has
a strong impact on the result of the synthesis process. The computational
results convincingly show the need for robust protection. In addition, it
turns out that a robust process yields products of desired sizes while keeping
the cost of robust protection very low.

### Related Literature

The impact of dispersity on product performance is particularly critical for
small nanoparticles and new materials, where slight variations can
significantly affect properties such as band gap and quantum yield [1, 27, 20,
5, 30, 19]. However, the mathematical optimization of robust nanoparticle
design is a yet largely unexplored area, and we aim at contributing to filling
this gap. Mathematical optimization provides a unique opportunity to establish
optimized processes that ensure product quality under uncertainty which is, in
particular, relevant for up-scaling the processes. The PBEs that are used to
model particle precipitation are nonlinear PDEs for which typical solution
approaches are finite-volume type methods [24, 10], methods based on
characteristics [8, 29], and methods based on moments of the solution such as
MoM [26], DQ-MoM [18], and eMoM [22] which is based on theoretical insights
from [12, 13]. In this work, we set up an optimization problem for particle
precipitation that uses the eMoM representation. We also go one step further
and hedge the optimized process against uncertainties.

Indeed, in many applications, uncertainties largely impact the quality of the
obtained solutions and can easily render an obtained result useless. Typical
reasons are that uncertainties can lead to high costs, to a low-quality
solution, or the solutions obtained may even be infeasible as they do not
satisfy the side constraints that are affected by the uncertainty. Several
solution approaches have been presented in the literature. Stochastic
optimization is a well-known paradigm used to address optimization problems
with uncertain parameters. This approach involves the use of random variables,
such as random objective functions, random constraints, and stochastic
quantities like expected values or higher moments. In various literature, such
as [23], [4], and [25], these concepts are discussed in detail. It should be
noted that stochastic optimization requires knowledge of the distributions of
uncertain parameters, which is often not available. In addition, protection
against uncertainties is only guaranteed in a probabilistic sense and with a
certain probability. In contrast, in the case of nanoparticle design, it is
essential to guarantee quality requirements with certainty. Moreover, the
underlying distributions are themselves uncertain. This requires a different
and robust approach.

Robust optimization does not require knowledge of the underlying probability
distributions. It addresses these needs by predefining so-called uncertainty
sets against which protection is sought. The resulting robust problem ensures
feasibility for all realizations of the uncertain parameters within the
uncertainty set and optimizes the guaranteed cost. For continuous uncertainty
sets, robust optimization models are typically represented as semi-infinite
problems (SIPs) with finitely many variables and infinitely many constraints.
In order to obtain algorithmically tractable robust counterparts of these
SIPs, elegant theoretical as well as algorithmic concepts have been developed.
For a focus on robust linear and convex optimization where - among others -
reformulations to tractable counterparts have been presented, we refer to the
books [2], [3]. Robust optimization approaches have been developed and applied
in the context of process optimization, e.g., in [17, 31, 14]. In particular,
for a restricted synthesis process assuming seeded growth and linearized PBEs,
a convex robust optimization task could be established, together with an
algorithmically tractable reformulated robust counterpart in [7].

More generally, when including nucleation and growth and using the general
nonlinear PBE, we face nonlinear and non-convex robust optimization problems.
For this challenging class of optimization problems, reformulation approaches
are typically not available, and approximation schemes, such as [6, 11], or
decomposition methods are a potential choice. Indeed, only a few general
approaches are known for this setting. In this work, we enhance a specific
decomposition algorithm. The latter consists in a novel adaptive bundle method
that has been developed as a general framework for nonlinear robust
optimization in [15]. The method has been integrated into an outer
approximation scheme for discrete-continuous decisions in [16]. In this work,
we make it concrete for optimal particle precipitation. To this end, we
present the adaptive bundle method, which integrates the eMoM approach and
uses it for a practically efficient solution of the underlying PBE.

### Outline

We introduce the setting in Section 2 and summarize the mathematical model for
particle synthesis via PBEs including nucleation and growth. We revisit the
solution approach via eMoM [22] and also introduce the nominal optimization
problem for particle precipitation. Section 3 extends the model by considering
mass uncertainties and presents the robust optimization problem that hedges
the precipitation process against uncertainties. Section 4 summarizes the
adaptive bundle method [15] and presents how to integrate the eMoM that is
used for the solution of the PBE. Computational results for the unprotected as
well as for the robust synthesis process are presented and discussed in 6. We
conclude with a discussion and an outlook for future research.

##  2 Optimized Synthesis modeled by Population Balance Equations

The evolution of the distribution of particle sizes is based on nucleation and
growth, which are driven by the supersaturation of the reactant
concentrations. Particle precipitation can be described by a PBE. With the
diameter of crystals with spherical shape denoted by $x$, process time by $t$
and the particle size distribution denoted by $q$ is a solution of this PBE,
see [22] and [9].

An important goal in particle design is to establish the best possible
synthesis processes with respect to an objective functional, such as
maximizing the amount of product that satisfies specific size or, more
generally, quality criteria. We will later use the time-dependent inflow rate
of precursor as the control variable because it has a significant effect on
the final particle size distribution.

###  2.1 Population Balance Equation

The precursor concentration $c_{pre}$ and the total concentration of reduced
educt species $c_{tot}$ are governed by the control $v$, i.e. the precursor
inflow rate, and prescribed by the following set of ordinary differential
equations:

$\displaystyle\partial_{t}c_{pre}(t)$ $\displaystyle=v(t)-k_{r}c_{pre}(t),$
$\displaystyle c_{pre}(0)$ $\displaystyle=0,$ (1)
$\displaystyle\partial_{t}c_{tot}(t)$ $\displaystyle=k_{r}c_{pre}(t),$
$\displaystyle c_{tot}(0)$ $\displaystyle=0,$

where $k_{r}\in\mathbbm{R}_{>0}$ is the reduction rate of the precursor
species.

Let $\mathcal{G}_{0}$ be the concentration-dependent and $\mathcal{G}_{1}$ the
size-dependent part of the growth rate, $q_{0}$ the initial number density
distribution, $\mathcal{N}$ the concentration-dependent nucleation rate, and
$x_{n}\in\mathbb{R}_{>0}$ the nucleation size of particles. The solution of
the following system of equations then gives the number density distribution
$q$:

$\displaystyle\partial_{t}q(t,x)+\partial_{x}\big{(}\mathcal{G}_{0}(c(t))\mathcal{G}_{1}(x)q(t,x)\big{)}$
$\displaystyle=0,$ (2) $\displaystyle q(0,x)$ $\displaystyle=q_{0}(x),$
$\displaystyle\mathcal{G}_{0}\left(c(t)\right)\mathcal{G}_{1}(x_{n})q(t,x_{n})$
$\displaystyle=\mathcal{N}\left(c(t)\right).$

We consider only nucleation and no preexisting particles in this contribution,
thus $q_{0}\equiv 0$. The total concentration $c_{tot}$, the solution in the
liquid phase $c$, and the PSD $q$ are connected by the mass balance as
follows:

$\displaystyle
c(t)=c_{tot}(t)-\tfrac{\rho\pi}{6V}\int_{x_{n}}^{\infty}x^{3}q(t,x)\,\mathrm{d}x.$
(3)

where $V$ denotes the total volume of the system, $x_{n}$ represents the
nucleation size, while $\rho>0$ signifies the nanoparticle density. As the
system of equations governing the evolution of the precursor concentration and
the total educt concentration can be solved analytically, we obtain:

$\displaystyle
c(t)=\int_{0}^{t}\\!\\!\big{(}1-e^{-k_{r}(t-\tau)}\big{)}v(\tau)\,\mathrm{d}\tau-\tfrac{\rho\pi}{6V}\\!\\!\int_{x_{n}}^{\infty}\\!\\!\\!\\!q(t,x)x^{3}\,\mathrm{d}x.$
(4)

In addition, it is assumed that the size-dependent growth rate
$\mathcal{G}_{1}$ lies in the class of growth kinetics of the form
$\mathcal{G}_{1}=x^{\beta}$, for $\beta\in\mathbbm{R}\setminus\\{1\\}$. The
examples presented in this paper focus on diffusion-limited growth, i.e.
$\beta=-1$. However, we would like to point out that the generality of the
presented approach also covers other relevant growth kinetics in nanoparticle
precipitation, such as those shown in [28].

###  2.2 Exact Method of Moments

The exact method of moments, as derived in [22], approximates the evolution of
the concentration in the liquid phase $c$, instead of solving the PBE. For
$t\in\left[0,t_{n}\right]$, we can define:

$\displaystyle c(t)$
$\displaystyle=\int_{0}^{t}\big{(}1-e^{-k_{r}(t-\tau)}\big{)}v(\tau)\,\mathrm{d}\tau$
(5)
$\displaystyle\qquad-\tfrac{\rho\pi}{6V}\int_{0}^{t}\mathcal{N}\left(c(\tau)\right)\psi\left[\tau,x_{n}\right](t)^{3}\,\mathrm{d}\tau,$

where

$\displaystyle\psi\left[\tau,x\right](t)$
$\displaystyle:=\left(x^{1-\beta}+(1-\beta)\int_{\tau}^{t}\mathcal{G}_{0}(c(s))\,\mathrm{d}s\right)^{\frac{1}{1-\beta}}.$
(6)

Here, $\psi\left[\tau,x\right](t)$ describes the size of a particle at the
time $t$, with the condition that this particle was of size $x$ at the time
$\tau$. As shown in [22], given a solution $c$ of (5), we obtain at a solution
formula for the PBE as follows.

For $(t,x)\in\mathbbm{R}_{>0}\times\mathbb{R}_{\geq x_{n}}$, the solution of
the PBE in eq. 2 is given by:

$\displaystyle q(t,x)$
$\displaystyle=x^{-\beta}\begin{cases}\frac{\mathcal{N}\left(c(\tau_{t,x})\right)}{\mathcal{G}_{0}\left(c(\tau_{t,x})\right)}&\text{
for }x\leq\psi\left[0,x_{n}\right](t)\\\ 0&\text{otherwise}\end{cases},$ (7)

with $c$ satisfying equation (5) and

$\tau_{t,x}:=\max\left\\{\nu\in\left[0,t\right]:\psi\left[\nu,x_{n}\right](t)=x\right\\},$
(8)

for all $t\in\left[0,t_{n}\right]$ and $x\leq\psi_{m}\left[0,x_{n}\right](t)$.

In order to formulate desired properties of the final PSD in terms of their
moments, e.g. the mean and the variance of the final PSD, the moments of the
PSD can also be given directly from the concentration profile. This means that
the p-th moment $m_{p}$ of the number density distribution $q$ at time
$t\in\mathbbm{R}_{\geq 0}$ is given by:

$\displaystyle m_{p}(t):$
$\displaystyle=\int_{x_{n}}^{\infty}q(t,x)x^{p}\,\mathrm{d}x,$ (9)
$\displaystyle=\int_{0}^{t}\mathcal{N}\left(c(\tau)\right)\psi\left[\tau,x_{n}\right](t)^{p}\,\mathrm{d}\tau,$

where $\psi$ is again given by (6).

###  2.3 An Optimization Model for Particle Precipitation

In order to establish an optimization problem for synthesis processes leading
to tailored PSDs, we aim here to optimize the number-weighted mean and
variance of $q$ at the end of the synthesis process. More specifically, the
goal is to control the precursor feed to the system so that the deviation of
the mean of the distribution from a desired mean and the deviation of the
standard deviation of the distribution from a desired standard deviation are
minimized.

Let $\bar{m}_{0}$ be the desired number-weighted mean and assume that we want
to control the standard deviation to zero. Then, the objective functional $J$
is given by:

$J:=\omega_{1}\left(\tfrac{m_{1}(T)}{m_{0}(T)}-\bar{m}_{0}\right)^{2}\\!\\!+\omega_{2}\left(\tfrac{m_{2}(T)}{m_{0}(T)}-\left(\tfrac{m_{1}(T)}{m_{0}(T)}\right)^{2}\right).$
(10)

Note that the squared standard deviation, i.e. the variance, of a random
variable $X$ is given by $E[X^{2}]-E[X]^{2}$ while the weights
$\omega_{1},\omega_{2}\in\mathbbm{R}_{\geq 0}$ scalarize the two criteria.

####  2.3.1 Discretization

In order to numerically approximate the solution of the PBE and its moments,
we use a discretization analoguously to [22] for (5) and (9). For this, we
discretize the time horizon into $N_{t}\in\mathbbm{N}_{>0}$ time-points
$0=\mathbf{t}_{0}<\mathbf{t}_{1}<\ldots<\mathbf{t}_{N_{t}}=T$ and define the
length of each time-interval by
$\boldsymbol{\delta}_{k}:=\mathbf{t}_{k+1}-\mathbf{t}_{k}$. Based on this, we
approximate the control $v$ by its values at these time-points, i.e.
$\mathbf{v}_{i}=v(\mathbf{t}_{i})\ \forall i$. In this way, we can approximate
the total mass in the system $c_{tot}$ as follows:

$\displaystyle\mathbf{c}_{tot,k}$
$\displaystyle=\sum_{\ell=1}^{k-1}\mathbf{v}_{\ell}\Big{(}\boldsymbol{\delta}_{\ell}+\tfrac{1-e^{k_{r}\boldsymbol{\delta}_{\ell}}}{k_{r}e^{k_{r}(\mathbf{t}_{k}-\mathbf{t}_{\ell})}}\Big{)},$
(11)

and based on this the evolution of the approximated concentrations
$\mathbf{c}_{k}\approx c(\mathbf{t}_{k})$ as follows:

$\displaystyle\mathbf{c}_{k+1}$
$\displaystyle=\mathbf{c}_{k}+\gamma_{1}\big{(}\mathbf{c}_{tot,k+1}-\mathbf{c}_{tot,k}\big{)}-\gamma_{2}\mathcal{N}(\mathbf{c}_{k})\boldsymbol{\delta}_{k}\mathcal{G}_{k,k}^{3}$
(12)
$\displaystyle\qquad-\gamma_{2}\sum_{\ell=1}^{k-1}\mathcal{N}(\mathbf{c}_{\ell})\boldsymbol{\delta}_{\ell}\big{(}\mathcal{G}_{\ell,k}^{3}-\mathcal{G}_{\ell,k-1}^{3}\big{)},$
$\displaystyle\mathbf{c}_{1}$ $\displaystyle=0,$

with the approximated characteristics:

$\displaystyle\mathcal{G}_{\ell,k}:=\Big{(}x_{n}^{1-\beta}+(1-\beta)\sum_{m=\ell}^{k}\mathcal{G}_{0}(\mathbf{c}_{m})\delta_{m}\Big{)}^{\frac{1}{1-\beta}}.$
(13)

All of these discretizations are based on dividing the integral involved into
the time intervals and then approximating the integrand by its value at the
lower bound of each interval.

Combining the derived governing equations of the discretized concentration
$\mathbf{c}$ based on the discretized control $\mathbf{v}$, we obtain:

$\displaystyle F(\mathbf{v},\mathbf{c})=0,$ (14)

with
$F:\mathbbm{R}^{N_{t}}\times\mathbbm{R}^{N_{t}}\mapsto\mathbbm{R}^{N_{t}}$.
Let $C:\mathbbm{R}^{N_{t}}\mapsto\mathbbm{R}^{N_{t}}$ be the mapping from the
control $\mathbf{v}$ to the respective concentration $\mathbf{c}$, i.e.
$F(\mathbf{v},C(\mathbf{v}))=0$.

####  2.3.2 Nominal Optimization Problem

Let

$\displaystyle\mathcal{V}:=\left\\{\mathbf{v}\in[l,u]^{N_{t}}:\textstyle\sum\limits_{k=1}^{N_{t}}\boldsymbol{\delta}_{k}\mathbf{v}_{k}=V_{tot}\right\\},$
(15)

be the set of admissible controls, with lower and upper bounds $0\leq l<u$ and
total added precursor $V_{tot}$. The optimization problem can be formulated in
the following way:

$\min_{\mathbf{v}\in\mathcal{V}}J(\mathbf{v},C(\mathbf{v})),$ (16)

where, as stated above, $C(\mathbf{v})\in\mathbb{R}^{n}$ satisfies
$F(\mathbf{v},C(\mathbf{v}))=0$.

To solve this optimization problem by first-order methods, we apply the
implicit function theorem to $F$ to obtain the derivative:

$\nabla
C(\mathbf{v})=-\nabla_{2}F(\mathbf{v},C(\mathbf{v}))^{-1}\nabla_{1}F(\mathbf{v},C(\mathbf{v})),$

and consequently

$\displaystyle\nabla_{\mathbf{v}}J(\mathbf{v},C(\mathbf{v}))=$
$\displaystyle\qquad\nabla_{1}J(\mathbf{v},C(\mathbf{v}))+\nabla_{2}J(\mathbf{v},C(\mathbf{v}))\nabla
C(\mathbf{v})).$ (17)

##  3 Particle Precipitation under Uncertainties

We acknowledge that in a synthesis process set up in a laboratory, it is not
possible to exactly control precursor inflow rates. This means that the
resulting uncertainties can have a significant impact on the outcome of the
process. Therefore, our focus will be primarily on adjusting this parameter.
As a consequence, we will formulate a robust optimization problem. To solve
it, we will extend an adaptive bundle framework from [15] and combine it with
the exact methods of moments. To model uncertainties in the mass flow, we use
a multiplicative uncertainty set that captures relative parameter variations.
The latter allows us to relate the parameter variations to their magnitude,
taking into account that a smaller parameter may have smaller variations
compared to a larger parameter. We define our model according to Section 2 in
a time-discretized setting. The aim is to find a solution that is robustly
protected against relevant deviations from a nominal, i.e. typical, scenario.

Thus, the uncertainty set naturally includes the nominal scenario (without
disturbance), which corresponds to a factor of $1$ for each time point. A
commonly used uncertainty set is given by a $\|\cdot\|_{\infty}$-bounded
function of the time interval, resulting in a box-shaped uncertainty set.
Assume that the worst cases scenario if of ”bang-bang” type, i.e., the worst
cases scenarios can be found in a vertex of the uncertainty set. If we further
assume that pushing the inflow rate in the nucleation phase at the beginning
of the nanoparticle precipitation and the growth phase in the second part of
the synthesis to opposite boundaries causes the strongest changes in the
objective functional we end up with a monotone uncertainty with values of
either the upper or the lower bound. To also encounter for the nominal case,
i.e. uncertainty $\mathbf{u}=1$, we consider all monotone perturbations with
function values $u_{l},1,u_{u}$ and a single jump within the process time.

Figure 1: Conceptual visualization of the used uncertainty set (18)

In more detail, we define the uncertainty set as follows.

Let $n$ be the number of discretization steps and $u_{l},u_{u}\in\mathbb{R}$
such that $0<u_{l}<1<u_{u}<\infty$. The uncertainty set is given by:

$\displaystyle\mathcal{U}:=\Big{\\{}$
$\displaystyle\mathbf{u}\in\mathbb{R}^{n}|u_{1},u_{2}\in\\{u_{l},1,u_{u}\\},j\in\left\\{1,...,n\right\\},$
(18) $\displaystyle\mathbf{u}_{i}=u_{1}\text{ for }i\leq
j,\mathbf{u}_{i}=u_{2}\text{ for }i>j\Big{\\}}.$ (19)

This uncertainty set then represents the basic building block for more general
uncertainty sets where more than one change in uncertainty is manifested over
the full-time period. The uncertainty considered here can easily be extended
in the future to account for multiple changes in the factor over time.

Next, we obtain a robust optimization model as follows. Let $\mathcal{V}$ be
the compact set of admissible discretized controls as defined in eq. 15 and
the cost function $J$ be chosen as in eq. 10 depending on the discretized
state operator $F$ and the discretized concentration function $C$ introduced
in Section 2.3.1 above. Then the robust optimization problem is:

$\displaystyle\min_{\mathbf{v}\in\mathcal{V}}\max_{\mathbf{u}\in\mathcal{U}}$
$\displaystyle J(\mathbf{u}\bullet\mathbf{v},C(\mathbf{u}\bullet\mathbf{v})),$
(20)

where $\mathbf{u}\bullet\mathbf{v}$ denotes the entrywise product of
$\mathbf{u}$ and $\mathbf{v}$, i.e.,

$(\mathbf{u}\bullet\mathbf{v})_{i}=\mathbf{u}_{i}\mathbf{v}_{i}\qquad\forall
i=1,\ldots,n.$

Moreover, as before, the concentration $C(\mathbf{u}\bullet\mathbf{v})$ is
implicitly defined as the unique solution of the discretized state equation:

$\displaystyle
F(\mathbf{v}\bullet\mathbf{u},C(\mathbf{u}\bullet\mathbf{v}))=0.$ (21)

Since this problem is nonlinear and non-convex, there is no standard procedure
for a practical or theoretical solution. In the next section we will focus on
its solution via an adaptive bundle framework by [15].

##  4 Robust Optimization via Bundle Method for Nonconvex Optimization

We tackle the robust PSD problem (20) by the adaptive bundle method from [15].
This bundle method is applicable to robust optimization problems of the form:

$\min_{\mathbf{v}\in\mathcal{V}}\max_{\mathbf{u}\in\mathcal{U}}h(\mathbf{v},\mathbf{u}),$
(22)

with a compact uncertainty set $\mathcal{U}$ and a cost function
$h:\mathbb{R}^{n+m}\rightarrow\mathbb{R}$ that is locally Lipschitz continuous
and lower $C^{1}$. Choosing
$h(\mathbf{v},\mathbf{u}):=J(\mathbf{u}\bullet\mathbf{v},C(\mathbf{u}\bullet\mathbf{v}))$,
it is readily seen that the robust PSD (20) fits into this setting.

More generally, the adaptive bundle method solves the problem (22) by treating
it as a nonsmooth minimization problem with the optimal value function:

$\hat{h}(\mathbf{v}):=\max_{\mathbf{u}\in\mathcal{U}}h(\mathbf{v},\mathbf{u}),$

as an objective function. An evaluation of the objective function, therefore,
requires solving the inbuilt maximization problem. This maximization problem
aims to find the worst-case realization of the uncertain parameter. We thus
call it the adversarial problem. For the PSD problem, the adversarial problem
is to find the value of the uncertain parameter $\mathbf{u}$ that maximizes
the deviation from the desired volume-weighed mean and standard deviation as
defined in eq. 10.

The bundle method requires access to approximate solutions to the adversarial
maximization problem in the following sense.

######  Assumption 1.

We assume that for every $\mathbf{v}\in\mathcal{V}$ and every choice of
$\epsilon_{f}>0$, we have access to a $\mathbf{u}_{\mathbf{v}}\in\mathcal{U}$
such that
$h(\mathbf{v},\mathbf{u}_{\mathbf{v}})\geq\max_{\mathbf{u}\in\mathcal{U}}h(\mathbf{v},\mathbf{u})-\epsilon_{f}$,
to the corresponding function value $h(\mathbf{v},\mathbf{u}_{\mathbf{v}})$
and to a subgradient
$g_{\mathbf{v}}\in\partial_{\mathbf{v}}h(\mathbf{v},\mathbf{u}_{\mathbf{v}})$.

Therefore, it is necessary to compute an approximate worst-case for every
decision and tolerance. The bundle method does not require any assumptions on
the structure of the underlying optimization models, which allows it to be
used for nonlinear robust optimization. In addition, it is also not restricted
to finite-dimensional optimization problems but can be extended more generally
to robust optimization of models containing PDEs, as we do here for robust
particle precipitation, including PBE. According to the assumption, the method
requires an approximate worst-case realization for a given tolerance
$\epsilon_{f}$:

$\mathbf{u}_{\mathbf{v}}\in\left\\{\mathbf{u}\in\mathcal{U}\ |\
h(\mathbf{u},\mathbf{u}_{\mathbf{v}})\geq\max_{\mathbf{u}\in\mathcal{U}}h(\mathbf{u},\mathbf{u}_{\mathbf{v}})-\epsilon_{f}\right\\}.$
(23)

For the first-order information required by the bundle method, we can use an
approximate subgradient:

$g_{\mathbf{v}}\in\partial_{\mathbf{v}}h(\mathbf{v},\mathbf{u}_{\mathbf{v}}).$
(24)

Note that in this paper the error $\epsilon_{f}$ can be chosen to be $0$ since
we only allow jumps of the uncertainty function at a finite number of points,
which are given by the time discretization. As a consequence of this, the
required subgradient is computed directly by evaluating section 2.3.2 in the
perturbed control $(\mathbf{v},\mathbf{u}_{\mathbf{v}})$.

We mention, however, that the adaptivity in the bundle approach can be
exploited in future work to speed up computations for very fine time
discretizations, or, even more interestingly, if a continuous uncertainty set
were used. An example of the latter would be obtained, if we allowed our
monotonous uncertainty functions to jump at arbitrary points. As will be shown
in the computational results, the approach used for our purposes already leads
to interesting results, so these algorithmic improvements are postponed to
future research.

We will now explain the algorithmic concept of the bundle method. The
iterations are divided into inner and outer loops. The outer loops generate
serious iterates while the inner loops generate trial iterates around the
current serious iterate. These are then potential candidates for the next
serious iterate. To distinguish between the inner and outer iterations in the
algorithm and in the following statements, we denote the inner loop counter by
$k$ and the respective trial iterations by $\mathbf{v}^{k}$, as well as the
outer loop counter by $j$ and the respective serious iterations by
$\mathbf{v}_{j}$. The algorithm also contains the proximity control parameter
$\tau_{k}$, which depends on the current quality of the working model. Details
can be found in [15]. We initiate the algorithm by running the nominal
optimization. The solution is our first serious iteration $\mathbf{v}_{1}$.

New trial iterates are generated by optimizing a convex working model. The
working model around a serious iterate $\mathbf{v}$ is constructed as a
maximum of cutting planes, i.e. as a piecewise linear model. To do this, the
method requires approximate function values and subgradients. For each trial
iteration, we therefore solve the adversarial maximization problem up to a
given error. The required error bound is an adaptive one and depends on the
current distance between iterations.

Algorithm 1 Adaptive bundle method [15]

1:Initial $\mathbf{v}_{0}$: Run nominal optimization. Set $j=1$

2:while dist$(0,\tilde{\partial}_{a}\hat{h}(\mathbf{v}_{j})>\tilde{\epsilon}$
do

3: Start inner loop with serious iterate $\mathbf{v}_{j}$; set $k=1$.

4: Find new trial iterate $\mathbf{v}^{k}$ by solving
$\min_{\mathbf{v}^{k}\in\mathcal{V}}\Phi_{k}(\mathbf{v}^{k},\mathbf{v}_{j})+\frac{\tau_{k}}{2}||\mathbf{v}^{k}-\mathbf{v}_{j}||^{2}$

5: Approximate $\hat{h}(\mathbf{v}^{k})$

6: if accept then

7: $\mathbf{v}_{j+1}\leftarrow\mathbf{v}^{k}$

8: else

9: Generate cutting plane with
$g_{k}\in\tilde{\partial}_{a}\hat{h}(\mathbf{v}^{k})$

10: Build a new working model

11: Increase $k$ and go to Step 5.

12: Increase j.

##  5 Implementation

The adaptive bundle method was implemented in MATLAB. The experiments were
performed on a machine with an Intel Core i7-8550U (4 cores) and 16GB of RAM.
For the computation of our robust solutions, we encountered running times
between $5$ and $20$ minutes.

The initial value for the bundle method was chosen as the result of the
nominal optimization. The algorithm used for this is an adjusted version of
[9].
In step 5 of the bundle algorithm, we solve the adversarial maximization
problem over the (discretized) uncertainty set by a straight-forward
enumeration of all possible discretized values.

##  6 Computational Results

In this section, we present the numerical results for optimizing the particle
synthesis process.

Table 1 shows the parameter values that were used in the optimization model.
To choose realistic values, we assume that the quantum dots (QDs) should have
a mean value of 4 nm, and that the process should end after 10 minutes. The
other values can be found in the table. The optimization goal is to minimize
the quadratic difference of the mean size and the variance from a given value,
where their effects are weighted by two parameters $\omega_{1},\omega_{2}$,
see (10).

In the computational evaluations presented here, we weight both contributions
to the objective equally, see table. However, it turns out that the results
are qualitatively similar even for different choices, as can be seen in
Appendix 8.1, where we present results for a different choice.

name | symbol | value | unit
---|---|---|---
desired NP size | $\bar{y}$ | 4 | nm
final control time | $t_{n}$ | 10 | min
nuclei size | $x_{\star}$ | 1 | nm
# control points | n | 100 | -
growth law exponent | $\beta$ | -1 | -
nanoparticle density | $\rho$ | 1 | -
obj. func. factor size | $\omega_{1}$ | 1 | -
obj. func. factor var. | $\omega_{2}$ | 1 | -
Table 1: Optimization parameters as used in the computational results.

In particular, we focus on two optimized solutions: the nominal optimum
synthesis process (16) from Section 2.3.2 without protection against
uncertainties, and the robust process (20) which is designed to hedge against
uncertainties. We are interested in evaluating these processes when used in
the nominal, i.e., unperturbed, scenario which represents ideal laboratory
conditions where all $\mathbf{u}_{i}$ values are set to 1. In addition, we
study perturbed scenarios where uncertainties come into play. To study extreme
cases, we analyze the worst scenario for both the nominal and robust synthesis
processes. The nominal worst case corresponds to the uncertainty realization
that results in the largest objective function value for the nominal process.
Similarly, the robust worst case represents the uncertainty where the robust
process achieves the worst objective function value.

Figure 2:  Objective function values for the uncertainty set with
$[u_{l},u_{u}]=[0.9,1.1]$. Lines are guides for the eyes only.

Figure 10 illustrates the objective function values for the nominal and robust
processes, where the robust process is obtained with a maximum uncertainty of
10% ($[u_{l},u_{u}]=[0.9,1.1]$). The results are presented for three
scenarios: the nominal, the nominal worst case, and the robust worst case.

In the nominal scenario, as expected, the nominal process, specifically
optimized for this case, outperforms the robust one. However, the situation
changes significantly when considering perturbed scenarios, particularly when
the uncertainty manifests within the 10% uncertainty set to its nominal worst
case. Assuming that this scenario manifests itself, the nominal objective
value then deteriorates by more than a factor of two, when compared to the
nominal scenario. This indicates that the nominal process is highly sensitive
to uncertainty. When compared to the robust process, it is still approximately
$\frac{1}{3}$ worse than the robust one. Conversely, the robust process
demonstrates stable objective function values that exhibit only slight changes
from the nominal to the robust worst-case scenario. Also in the robust worst-
case scenario, the robust process performs slightly better than the nominal
one.

Figure 3:  Objective function values for both the nominal and the robust
process, are evaluated in different scenarios, as a function of the size of
the uncertainty set. The size of the uncertainty is given in percent deviation
from the nominal scenario.

Figure 3 displays the objectives for the nominal process, evaluated both for
the nominal and the nominal worst-case scenario, and for the robust process,
evaluated both for the nominal and the robust worst-case scenario. It is clear
that the nominal optimum process (blue line) stays at a constant objective
value in the nominal, i.e., unperturbed, scenario, irrespective of the size of
the uncertainty set. However, when the nominal process is evaluated in the
worst nominal scenarios, the objective value of the latter significantly
increases to roughly a factor of 5 for the largest uncertainty of 20%. This
means that if this uncertainty scenario is realised, the quality of the
synthesised product obtained will be significantly reduced, as the mean and
variance of the products obtained will be significantly different from the
given values. Comparing these results with the robust process (red lines), the
obtained robust objectives are significantly better than those of the nominal
processes. In fact, the robust process basically coincides with the nominal
process in the nominal scenario for an uncertainty size of up to 5%. Only at
about 10% uncertainty the objective reaches larger values than the nominal
process in the nominal scenario. Even for the largest uncertainty of 20%, the
robust objective for the worst robust scenario is at most about a factor of
two worse than that of the nominal process in the nominal setting. This shows
that the cost of robustness, i.e., the loss of protection against uncertainty,
is low.

Figure 4: The mass flow for the nominal process (blue), for the nominal
scenario and the nominal worst-case scenario, and for the robust synthesis
process (red), for the nominal scenario and the nominal worst-case scenario.
The nominal worst-case scenario is a change of the factor from $0.9$ to $1.1$
at time $2.2$ and the robust worst-case scenario is a change of the factor
from $1.1$ to $0.9$ at time $3.9$.

Having considered the values of the objective function, we now take a closer
look at the processes over time. Figure 4 shows the mass flows as a function
of time in the optimized nominal and robust processes in their respective
nominal and worst-case scenarios. It is clear that the nominal (and robust)
solutions show similar patterns in both the nominal (robust) and worst-case
scenarios. However, there are notable differences when comparing the nominal
and robust processes. This is most prominent during process times between two
and four minutes, where the worst-case scenarios also manifest themselves
(nominal worst-case at time 2.2, robust worst case at time 3.9). In
particular, the robust process shows two additional peaks around these times
that are not present in the nominal process. Consequently, the robust process
displays distinct characteristics that set it apart from the nominal,
unprotected process.

Figure 5: The total mass for the nominal process (blue), for the nominal
scenario and the nominal worst-case scenario, and for the robust process
(red), for the nominal scenario and the nominal worst-case scenario.

Figure 5 illustrates the increase in total mass over time during the synthesis
process. It can be observed that the total mass increases continuously until
it reaches a saturation point at about 10 minutes, after which it remains
relatively constant.

Figure 6: The concentration for the nominal process (blue), for the nominal
scenario and the nominal worst-case scenario, and for the robust process
(red), for the nominal scenario and the nominal worst-case scenario.

In Figure 6 we present the concentration dynamics of the product over time
during the process. Initially, the concentration rises rapidly and reaches a
peak at approximately three minutes. It then decreases and stabilises at a
plateau between five and ten minutes. After ten minutes, the concentration
continues to fall and asymptotically reach the equilibrium concentration of 1.

Figure 7: The nucleation for the nominal process (blue), for the nominal
scenario and the nominal worst-case scenario, and for the robust process
(red), for the nominal scenario and the nominal worst-case scenario.

To study nucleation and growth, the nucleation function is shown in Figure 7.
Whereas it takes some time before nucleation takes place, nucleation is
present in the range of 2 to 5 minutes, exhibits a sharp peak around 3
minutes, and remains approximately at zero after about 5 minutes. The robust
solution starts nucleating a bit later and less pronounced, when compared to
the nominal process.

Figure 8: The growth rate for the nominal process (blue), for the nominal
scenario and the nominal worst-case scenario, and for the robust process
(red), for the nominal scenario and the nominal worst-case scenario.

To study the growth of the particle after nucleation has taken place, Figure 8
shows the growth rate over time. It initially increases rapidly and reaches a
peak at about 3 minutes. It then declines until it reaches roughly a plateau
at about 5 minutes. Thereafter, the particles continue to grow at a roughly
constant rate until about 10 minutes of process time. The growth rates differ
slightly between the nominal and the robust process between 5 and 10 minutes.
As the nominal growth rate slowly decreases slightly, the robust growth rate
increases slightly. Subsequently, all growth rates decrease rapidly to zero,
resulting in the end product having reached its final size distribution.

Figure 9: The particle size distribution for the nominal process (blue), for
the nominal scenario and the nominal worst-case scenario, and for the robust
process (red), for the nominal scenario and the nominal worst-case scenario.

Finally, we examine the particle size distributions in Figure 9, which shows
the particle number density as a function of particle size. The objective of
our optimization problem is to achieve a mean particle size of 4 nm, which is
successfully achieved by both the nominal and the robust process in the
nominal scenario. However, this changes drastically when evaluating the
processes in disturbed settings. The robust process then produces
nanoparticles with an average size of 3.85 nm, even in the worst case of
uncertainty. This mean value is therefore still close to the desired size of 4
nm, even under uncertainty. In comparison, applying the nominal process in the
nominal worst case results in significantly larger particles (4.58 nm). This
yields an end product that is larger than the desired mean, making it
potentially useless. In contrast, the robust process achieves the desired mean
size. This clearly demonstrates the superiority of the robust process, which
protects against uncertainties, over the unprotected nominal process.

Finally, we study the price that we need to pay when implementing the robust
process, which is typically referred to as the price of robustness. The latter
here means to compare the nominal process with the robust process when however
applied in the unperturbed nominal scenario, i.e., in the ideal case of
perfect laboratory conditions with no uncertainties. The mean values of the
PSDs almost coincide and are both about 4.15 nm for the nominal and robust
process. This means that with respect to the means of the PSDs, there is no
cost of robustness. However, when looking at their variances, it can be seen
that the robust process leads to a slightly larger variance than that of the
nominal process when the latter is used in the nominal scenario. In summary,
this demonstrates the low case of robust protection and its obvious advantages
when using the robust solution.

##  7 Conclusions and Outlook

In this work, we have introduced a robust optimization framework for particle
synthesis, including nucleation and growth, modeled with population balance
equations. We have formulated an optimization problem that focuses on the
time-dependent precursor inflow rate. It minimizes the quadratic deviation
from a given mean particle size and its variance. The resulting robust minimax
problem assumes that there are uncertainties in the inflow rate and hedges
against them. The numerical optimization was enabled by an advanced
integration of a recently introduced adaptive bundle framework, which can
solve nonconvex and nonlinear robust optimization problems, together with the
exact method of moments. The algorithm was applied to particle synthesis
processes with realistic parameter values taken from quantum dot synthesis.
The computational results showed that the bundle method delivers good results
in a short time while keeping the price of robustness low.

We have investigated the nominal and robust processes, specifically evaluating
them in their worst-case scenarios. It was found that the nominal solution is
highly sensitive to disturbances and is unable to synthesize a suitable
product under uncertainties. In contrast, the robust process successfully
produces particles of the desired size at a relatively low cost of robustness.
This is particularly evident in Figure 3, where the cost of the robust
solution is only marginally worse than that of the nominal solution, whereas
the cost of the nominal solution is very large in the nominal worst case. This
is even more striking in Figure 9, where it is clearly visible that the robust
process leads to better PSDs, while simultaneusly keeping the price of robust
protection low. In fact, even in the nominal scenario, which represents
unperturbed, ideal laboratory conditions, the robust procedure leads to
similar mean PSDs compared to the nominal solution, while the robust variance
is only slightly larger than that of the nominal one. The computational
results show that robust optimization of particle synthesis under uncertainty
is crucial for the production of high-quality particles. We emphasize that the
algorithm presented here is able to compute robust optimal processes in a
short time.

In the future, the adaptive bundle framework will be further extended by
adaptive discretization widths with the aim of allowing very fine
discretizations or continuous uncertainty sets. In conclusion, our
computational results show that the adaptive bundle method, together with the
exact methods of moments, is an efficient tool for dealing with uncertainties
in the nanoparticle synthesis process, which can be computed quickly and at a
low cost in terms of robustness. Finally, hedging mathematical optimization
problems against uncertainty, especially when partial differential equations
(PDEs) are involved, is also very relevant also in other contexts.
Applications are numerous, ranging from particle problems to energy problems
to problems in the natural sciences. The bundle method integrated by PDE-
related solution approaches, as presented here, is not limited to particle
synthesis. Thus, the approach to particle precipitation presented here will be
developed for other such challenging problems in the future.

CRediT author statement
Martina Kuchlbauer: Methodology, Software, Validation, Investigation, Writing,
Visualization
Jana Dienstbier: Methodology, Writing - Review & Editing
Adeel Muneer: Software, Writing - Review & Editing
Hanna Hedges: Software, Writing - First Draft
Michael Stingl: Conceptualization, Methodology, Investigation, Validation,
Writing - Review & Editing
Frauke Liers: Conceptualization, Methodology, Validation, Investigation,
Writing, Supervision
Lukas Pflug: Conceptualization, Software, Methodology, Validation, Formal
Analysis, Investigation, Writing, Visualization

Funding
This paper has been funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - Project-ID 416229255 - SFB 1411 and by Project B06 in
CRC TRR 154.

## References

  * [1] A. P. Alivisatos.  Semiconductor clusters, nanocrystals, and quantum dots.  science, 271(5251):933–937, 1996.
  * [2] A. Ben-Tal, L. El Ghaoui, and A. Nemirovski.  Robust Optimization.  Princeton Series in Applied Mathematics. Princeton University Press, October 2009.
  * [3] D. Bertsimas and D. Den Hertog.  Robust and Adaptive Optimization.  Dynamic Ideas LLC, Belmont, MA, 2022.
  * [4] J. Birge and F. Louveaux.  Introduction to Stochastic Programming.  Springer Series in Operations Research and Financial Engineering. Springer New York, 2006.
  * [5] C. Burda, X. Chen, R. Narayanan, and M. A. El-Sayed.  Chemistry and properties of nanocrystals of different shapes.  Chemical reviews, 105(4):1025–1102, 2005.
  * [6] M. Diehl, H. G. Bock, and E. Kostina.  An approximation technique for robust nonlinear optimization.  Math. Program., 107(1-2, Ser. B):213–230, 2006.
  * [7] J. Dienstbier, K.-M. Aigner, J. Rolfes, W. Peukert, D. Segets, L. Pflug, and F. Liers.  Robust optimization in nanoparticle technology: A proof of principle by quantum dot growth in a residence time reactor.  Computers & Chemical Engineering, 157:107618, 2022.
  * [8] F. Févotte and G. Févotte.  A method of characteristics for solving population balance equations (pbe) describing the adsorption of impurities during crystallization processes.  Chemical Engineering Science, 65(10):3191–3198, 2010.
  * [9] A. Gilch.  Optimization of nanoparticle synthesis using exact method of moments.  Master’s thesis, Friedrich-Alexander Universität Erlangen-Nürnberg, 2021\.
  * [10] R. Gunawan, I. Fusman, and R. D. Braatz.  High resolution algorithms for multidimensional population balance equations.  AIChE Journal, 50(11):2738–2749, 2004.
  * [11] B. Houska and M. Diehl.  Nonlinear robust optimization via sequential convex bilevel programming.  Math. Program., 142(1-2, Ser. A):539–577, 2013.
  * [12] A. Keimer and L. Pflug.  Existence, uniqueness and regularity results on nonlocal balance laws.  Journal of Differential Equations, 263(7):4023–4069, 2017.
  * [13] A. Keimer, L. Pflug, and M. Spinola.  Nonlocal scalar conservation laws on bounded domains and applications in traffic flow.  SIAM Journal on Mathematical Analysis, 50(6):6271–6306, 2018.
  * [14] R. W. Koller, L. A. Ricardez-Sandoval, and L. T. Biegler.  Stochastic back-off algorithm for simultaneous design, control, and scheduling of multiproduct systems under uncertainty.  AIChE Journal, 64(7):2379–2389, 2018.
  * [15] M. Kuchlbauer, F. Liers, and M. Stingl.  Adaptive bundle methods for nonlinear robust optimization.  INFORMS Journal on Computing, 34(4):2106–2124, 2022.
  * [16] M. Kuchlbauer, F. Liers, and M. Stingl.  Outer approximation for mixed-integer nonlinear robust optimization.  Journal of Optimization Theory and Applications, 195, 2022.
  * [17] Z. Li and M. G. Ierapetritou.  Robust optimization for process scheduling under uncertainty.  Industrial & Engineering Chemistry Research, 47(12):4148–4157, 2008.
  * [18] D. L. Marchisio and R. O. Fox.  Solution of population balance equations using the direct quadrature method of moments.  Journal of Aerosol Science, 36(1):43–73, 2005.
  * [19] M. L. Mastronardi, F. Maier-Flaig, D. Faulkner, E. J. Henderson, C. Kübel, U. Lemmer, and G. A. Ozin.  Size-dependent absolute quantum yields for size-separated colloidally-stable silicon nanocrystals.  Nano letters, 12(1):337–342, 2012.
  * [20] O. I. Micic, C. J. Curtis, K. M. Jones, J. R. Sprague, and A. J. Nozik.  Synthesis and characterization of inp quantum dots.  The Journal of Physical Chemistry, 98(19):4966–4969, 1994.
  * [21] M. Naito, T. Yokoyama, K. Hosokawa, and K. Nogi.  Nanoparticle technology handbook.  Elsevier, 2018.
  * [22] L. Pflug, T. Schikarski, A. Keimer, W. Peukert, and M. Stingl.  emom: Exact method of moments—nucleation and size dependent growth of nanoparticles.  Computers & Chemical Engineering, 136:106775, 2020.
  * [23] A. Prekopa, B. Vizvari, and T. Badics.  Programming under probabilistic constraint with discrete random variable.  In New trends in mathematical programming, pages 235–255. Springer, 1998.
  * [24] S. Qamar, G. Warnecke, and M. P. Elsner.  On the solution of population balances for nucleation, growth, aggregation and breakage processes.  Chemical Engineering Science, 64(9):2088–2095, 2009.
  * [25] A. Ruszczyński and A. Shapiro.  Stochastic Programming.  Handbooks in Operations Research and Management Science. Elsevier, 2003\.
  * [26] H.-C. Schwarzer, F. Schwertfirm, M. Manhart, H.-J. Schmid, and W. Peukert.  Predictive simulation of nanoparticle precipitation based on the population balance equation.  Chemical Engineering Science, 61(1):167–181, 2006.
  * [27] D. V. Talapin, A. L. Rogach, A. Kornowski, M. Haase, and H. Weller.  Highly luminescent monodisperse cdse and cdse/zns nanocrystals synthesized in a hexadecylamine- trioctylphosphine oxide- trioctylphospine mixture.  Nano letters, 1(4):207–211, 2001.
  * [28] N. T. Thanh, N. Maclean, and S. Mahiddine.  Mechanisms of nucleation and growth of nanoparticles in solution.  Chemical reviews, 114(15):7610–7630, 2014.
  * [29] S. M. ur Rehman and S. Qamar.  Application of the method of characteristics to population balance models considering growth and nucleation phenomena.  Applied Mathematics, 2014, 2014.
  * [30] R. Viswanatha, S. Sapra, B. Satpati, P. Satyam, B. Dev, and D. Sarma.  Understanding the quantum size effects in zno nanocrystals.  Journal of Materials Chemistry, 14(4):661–668, 2004.
  * [31] J. Wiebe, I. Cecílio, and R. Misener.  Robust optimization for the pooling problem.  Industrial & Engineering Chemistry Research, 58(28):12712–12722, 2019.

##  8 Appendix

###  8.1 Additional Computational Results

Figure 10: For the case of $\omega_{1}=10$, $\omega_{2}=1$: The mass flow for
the nominal solution (blue), for the nominal scenario and the nominal worst-
case scenario, and for the robust solution (red), for the nominal scenario and
the nominal worst-case scenario.  Figure 11: For the case of $\omega_{1}=10$,
$\omega_{2}=1$: The particle size distribution for the nominal solution
(blue), for the nominal scenario and the nominal worst-case scenario, and for
the robust solution (red), for the nominal scenario and the nominal worst-case
scenario.