# Sequential optimal experimental design for vapor-liquid equilibrium modeling

Martin Bubel1 Jochen Schmid1 Volodymyr Kozachynskyi2
Erik Esche2 Michael Bortz1
1Fraunhofer Institute for Industrial Mathematics (ITWM) 67663 Kaiserslautern
Germany
2Technische Universität Berlin Process Dynamics and Operations Group  10623
Berlin  Germany
<EMAIL_ADDRESS>
###### Abstract

We propose a general methodology of sequential locally optimal design of
experiments for explicit or implicit nonlinear models, as they abound in
chemical engineering and, in particular, in vapor-liquid equilibrium modeling.
As a sequential design method, our method iteratively alternates between
performing experiments, updating parameter estimates, and computing new
experiments. Specifically, our sequential design method computes a whole batch
of new experiments in each iteration and this batch of new experiments is
designed in a two-stage locally optimal manner. In essence, this means that in
every iteration the combined information content of the newly proposed
experiments and of the already performed experiments is maximized. In order to
solve these two-stage locally optimal design problems, a recent and efficient
adaptive discretization algorithm is used. We demonstrate the benefits of the
proposed methodology on the example of of the parameter estimation for the
non-random two-liquid model for narrow azeotropic vapor-liquid equilibria. As
it turns out, our sequential optimal design method requires substantially
fewer experiments than traditional factorial design to achieve the same model
precision and prediction quality. Consequently, our method can contribute to a
substantially reduced experimental effort in vapor-liquid equilibrium modeling
and beyond.

##  1 Introduction

In order to reliably simulate and optimize engineering processes, one needs
sufficiently precise mathematical models describing the functional
relationships between relevant input and output quantities of the considered
process. In many cases, these models come in the form of parametric models,
that is, their predictions $f(x,\theta)$ for the output quantity $y$ of
interest not only depend on the input quantity values $x$ but also on
additional model parameters $\theta$. A priori, the true value of these model
parameters is unknown and therefore has to be estimated from appropriate
experimental data. Specifically, at suitably selected input points
$x_{1},\dots,x_{n}$ one measures the values $y_{1},\dots,y_{n}$ of the output
quantity and then adjusts the model parameters $\theta$ in such a way that the
discrepancy between the model’s predictions $f(x_{i},\theta)$ and the
measurements $y_{i}$ becomes minimal in the least-squares sense.

Intuitively, it is clear that in general the quality of the thus obtained
parameter estimate strongly depends on the selected experimental points
$x_{1},\dots,x_{n}$ or, in other words, on where one measures and on how often
one measures. It is also clear that in practice one can often perform only a
fairly limited number $n$ of experiments, just because experiments in a lab
are typically time- and cost-intensive. Consequently, it is of crucial
importance to carefully plan and select the actually performed experiments. In
particular, it is important that (i) the parameter estimate and the resulting
predictions based on the selected experiments $x_{1},\dots,x_{n}$ are
sufficiently precise and that (ii) the number $n$ of experiments is affordably
small. Coming up with experimental plans that meet these criteria is the core
task of experimental design.

In the extensive literature on experimental design, two main approaches can be
distinguished, namely factorial experimental design on the one hand and
optimal experimental design on the other hand. In factorial design, one
spreads the experiments in a uniform and grid-like manner in the entire input
space and, as a consequence, ends up with a relatively large number of
experiments already for moderate input space dimensions. In optimal design, by
constrast, one tries to make experiments only at places in the input space
where enough information can be gained about the model parameters to be
estimated. Specifically, one tries to distribute the experimental points such
that the information content – and hence the precision – of the corresponding
parameter estimate becomes optimal. As a consequence, with optimal design one
can expect to obtain the same parameter and prediction quality as with
factorial design, but with substantially fewer experiments. And this effect
can be expected to be all the more pronounced, the more nonlinear the
considered model $f$ is, as a function of the inputs $x$ or of the model
parameters $\theta$. In light of this, optimal experimental design is a
promising approach especially for vapor-liquid equilibrium modeling with its
intrinsically nonlinear models for activity coefficients, for instance.

In the present paper, we propose a general methodology of sequential locally
optimal experimental design [10, 39, 15, 30] for general nonlinear models, and
these nonlinear models can be defined either by explicit or by implicit
algebraic equations. As a sequential design method, our method iteratively
alternates between carrying out experiments, estimating the parameters, and
designing new experiments. Specifically, the sequential design method proposed
here is characterized by the following key features, which also distinguish it
from the sequential design method recently proposed in [12].

A first essential feature of our method is that, in every iteration, it
computes a whole batch of new experiments, the maximal batch size being
specified by the user. In contrast, the design methodology of [12] computes
new experiments one by one (except for the initial iteration). We decided for
batches of new experiments because this will generally lead to fewer
alternations between the experimental and the computational phase of the
procedure and thus to a smoother workflow in practice. In particular, this is
beneficial when each experimental phase requires extensive set-up or ramp-up
times or when the experiments and the computations are performed by different
people at different institutions.

A second important feature of our method is that, in every iteration, the
batch of new experiments is computed in a locally optimal manner, namely by
solving a two-stage locally optimal experimental design problem. Specifically,
the new experiments are designed such that the combined information content of
the already performed experiments (stage one) and of the new experiments to be
performed next (stage two) becomes maximal. As usual in two-stage locally
optimal design, the combined information content is measured in terms of a
suitable design criterion $\Psi$ applied to the information matrix of a
suitably linearized version of the model $f$ at the combined experimental
design. Also, the linearization of $f$, in every iteration, is performed about
the current best parameter estimate, namely, about the least-squares parameter
estimate based on all experimental data collected until the current iteration.
As the design criterion $\Psi$, in turn, one can choose each of the usual
standard criteria provided only it is differentiable. In particular, this
includes the A- and the D-criterion, while it excludes essentially only the
E-criterion.

A third distinctive feature of our method is that it solves the occurring two-
stage locally optimal experimental design problems by means of the recent
adaptive discretization algorithm from [48], which is specifically tailored to
such design problems. And, in fact, this algorithm exhibits very good
convergence properties, both from a theoretical and from a practical point of
view. We terminate our iterative alternation between experiments, parameter
estimation, and locally optimal experimental design, if the total number of
already performed and of newly proposed experiments exceeds a user-specified
upper bound or if the newly proposed experiments do not differ substantially
anymore from already performed ones.

In order to demonstrate the benefits of the proposed sequential design
methodology for practical use, we apply it to the estimation of non-random
two-liquid parameters. As is well-known, these parameters describe the binary
interactions in the liquid phase of a multi-component mixture and therefore
they play an essential role in vapor-liquid equilibrium modeling. It is also
well-known that the estimation of these parameters is generally a challenging
task, especially for azeotropic systems. In fact, the corresponding sum of
squared errors often features flat local minima, multiple local minima, or
even multiple global minima [41, 22, 46].

In the present paper, we specifically consider the binary system consisting of
propanol and propyl acetate and we measure the vapor composition and the
temperature of this system at vapor-liquid equilibrium for different liquid
compositions and pressures. It is known that this system exhibits narrow
azeotropic behavior and therefore the estimation of its non-random two-liquid
parameters is already a rather challenging task. Applying our sequential
design method to this system, we end up with a locally optimal experimental
design consisting of $15$ experiments. And, as it turns out, these $15$
optimally designed experiments lead to the same prediction quality as $27$
experiments from a traditional factorial experimental design. So, compared to
traditional factorial design, our sequential locally optimal design method
substantially reduces the experimental effort – and thus helps the
experimenter to save substantial amounts of time and money – without impairing
the model prediction quality. Additionally, our method also helps the
experimenter in deciding what to do in case there still is some experimental
budget left, after all the experiments proposed by a first run of our
sequential design method have been performed. Indeed, by its very definition,
our method always indicates where making a batch of further experiments is
most promising and where it is not.

Apart from the papers [11, 12], we are not aware of any other works that apply
optimal design techniques to the estimation of thermodynamic property
parameters. It should be noticed that in the mentioned papers [11, 12] no real
lab experiments are performed, but simulated data are used instead, based on
literature values for the non-random two-liquid parameters to be estimated. In
sharp contrast to that, we perform actual lab experiments and only rely on the
corresponding real experimental data. In particular, we do not rely on
literature data or on literature values for the considered parameters at all.

##  2 Setting and preliminaries

In this section, we specify the setting and the basic terminology used
throughout the paper. We also collect those facts from parameter estimation
and experimental design that are relevant for the sequential locally optimal
experimental design method proposed in this paper.

###  2.1 Setting and basic terminology

In the entire paper, we assume that the considered functional relationship
between input quantities $x$ and output quantities $y$ can be described by the
predictions $f(x,\theta)$ of a parametric model

$\displaystyle f:\mathcal{X}\times\Theta\to\mathcal{Y},$ (2.1)

where $\mathcal{X}\subset\mathbb{R}^{d_{\mathrm{x}}}$ and
$\Theta\subset\mathbb{R}^{d_{\mathrm{\theta}}}$ are the sets of admissible
input values and of admissible model parameter values respectively, while
$\mathcal{Y}:=\mathbb{R}^{d_{\mathrm{y}}}$ is a superset of the set of
possibly observable output values.

Such a parametric model $f$ can be linear or nonlinear and, moreover, it can
be explicit or implicit. As usual, $f$ is called _linear_ iff it is affine-
linear w.r.t. $\theta$, that is,

$\displaystyle
f(x,\theta)=c(x)+J(x)\theta\qquad((x,\theta)\in\mathcal{X}\times\mathbb{R}^{d_{\mathrm{\theta}}})$
(2.2)

for some arbitrary functions $c:\mathcal{X}\to\mathbb{R}^{d_{\mathrm{y}}}$ and
$J:\mathcal{X}\to\mathbb{R}^{d_{\mathrm{y}}\times d_{\mathrm{\theta}}}$. In
particular, a linear model $f$ can be arbitrarily nonlinear w.r.t. the inputs
$x\in\mathcal{X}$. Simple examples of linear models are furnished by the
polynomial functions in $x$ with the monomial coefficients making up the model
parameter vector $\theta$. Additionally, $f$ is called _explicit_ iff it is
explicitly defined by a closed-form expression for $f(x,\theta)$ for every
$(x,\theta)\in\mathcal{X}\times\Theta$. And correspondingly, $f$ is called
_implicit_ iff one has no closed-form expression for $f(x,\theta)$ but instead
$f(x,\theta)$ for every $(x,\theta)\in\mathcal{X}\times\Theta$ is implicitly
defined as the solution $y=h(s(x,\theta),x,\theta)$ of a system of algebraic
equations of the form

$\displaystyle g(s,x,\theta)=0\qquad\text{and}\qquad y=h(s,x,\theta)$ (2.3)

with explicitly defined functions $g$ and $h$. In this system of equations,
the first equation (2.3.a) determines the state $s=s(x,\theta)$ as a function
of $(x,\theta)$, while the second equation (2.3.b) determines the output
quantity $y=h(s(x,\theta),x,\theta)$ as a function of the state.

Simple examples of explicit nonlinear models from chemical engineering are
given by the Antoine model for the saturation vapor pressure of a pure
substance or by the non-random two-liquid model for the activity coefficients
of a mixture of substances. Similarly, a simple chemical engineering example
of an implicit nonlinear model arises in the prediction of the bubble-point
temperature of a mixture of substances. With this particular example of an
implicit nonlinear model, we will illustrate the general methodology proposed
in this paper. See Section 4 below.

We point out, though, that the proposed general methodology is by no means
restricted to thermodynamic models but rather can be applied in a wide range
of circumstances, to explicit and implicit nonlinear models alike. All that
needs to be satisfied for the methodology to apply are a few rather mild
regularity and representability assumptions on the considered model $f$.
Specifically, the model $f$ needs to be sufficiently regular (Assumption 1)
and sufficiently representative (Assumption 2) in the sense that, up to
additive random measurement errors, all possible output measurement values are
explained by the model’s predictions with one specific model parameter value
$\theta^{*}$. Spelled out in detail, these assumptions take the following
form.

######  Assumption 1.

  * (i)

$\mathcal{X}$ is a compact subset of $\mathbb{R}^{d_{\mathrm{x}}}$ and
$\underline{\mathcal{X}}$ is a compact subset of $\mathcal{X}$

  * (ii)

$\Theta$ is a compact subset of $\mathbb{R}^{d_{\mathrm{\theta}}}$ that is
equal to the closure of its interior, in short,
$\overline{\Theta^{\circ}}=\Theta$

  * (iii)

$\Theta\ni\theta\mapsto f(x,\theta)$ is continuously differentiable for every
$x\in\mathcal{X}$ and $\mathcal{X}\ni x\mapsto
D_{\theta}f(x,\overline{\theta})$ is continuous for every
$\overline{\theta}\in\Theta$. As usual, $D_{\theta}f(x,\overline{\theta})$
denotes the derivative (Jacobi matrix) of $f(x,\cdot)$ at $\overline{\theta}$.

Conventionally, the subset $\underline{\mathcal{X}}$ of the input space
$\mathcal{X}$ from the assumption above is called the _experimental design
space_. It consists of all potential experiments, that is, all input points
$x\in\mathcal{X}$ at which the experimenter can measure the output quantity.
We point out that the experimental design space $\underline{\mathcal{X}}$ can
be discrete (finite) or continuous (infinite). In particular,
$\underline{\mathcal{X}}$ can be equal to $\mathcal{X}$, but in most
experimental design algorithms $\underline{\mathcal{X}}$ is chosen to be
finite for the sake of compuational feasibility [15, 30, 48, 12, 12]. We also
point out that the assumed regularity conditions are very mild and usually
easy to verify. Indeed, for an explicit model one has only to investigate the
explicit model parametrization $f$, and for an implicit model one has only to
investigate the implicit model-defining functions $g$ and $h$ and invoke the
implicit function theorem.

######  Assumption 2.

A parameter value $\theta^{*}\in\Theta$ exists such that for every
$n\in\mathbb{N}$ and every set of input values
$x_{1},\dots,x_{n}\in\underline{\mathcal{X}}$, the corresponding measured
values $y_{1},\dots,y_{n}\in\mathcal{Y}$ of the output quantity are given by
the predictions $f(x_{1},\theta^{*}),\dots,f(x_{n},\theta^{*})$ of the model
$f(\cdot,\theta^{*})$ up to independent normally distributed measurement
errors $\epsilon_{1},\dots,\epsilon_{n}$, that is,

$\displaystyle
y_{i}=f(x_{i},\theta^{*})+\epsilon_{i}\qquad(i\in\\{1,\dots,n\\})$ (2.4)

where $\epsilon_{1},\dots,\epsilon_{n}$ are realizations of independent and
normally distributed measurement errors
$\varepsilon_{1},\dots,\varepsilon_{n}$ having mean $0$ and a known positive
definite covariance matrix $\varsigma\in\mathbb{R}^{d_{\mathrm{y}}\times
d_{\mathrm{y}}}_{\mathrm{pd}}$.

Any value $\theta^{*}$ as in the above assumption is called (the) _true value
of the model parameter_ and, accordingly, the corresponding model
$f(\cdot,\theta^{*})$ is called (the) _true model_. It is important to notice
that the assumption does not require such a true parameter value to be known
to the modeler. In fact, quite to the contrary, it is always unknown in
practice and therefore has to be estimated. All the assumption requires is
that such a true (but unknown) parameter value exists and that the additive
random measurement errors are independent and normally distributed. See [6]
(Section 1.3.1) for a critical discussion and justification of these standard
assumptions in parameter estimation [6, 36] and in experimental design [15,
26, 30].

With the basic assumptions set out, we now formally introduce the notion of
(unweighted and weighted) experimental designs along with the pertinent
notation that will be continually used throughout the paper. An _experimental
design_ on $\underline{\mathcal{X}}$ simply is a tuple

$\displaystyle\tilde{x}=(x_{1},\dots,x_{n})\in\underline{\mathcal{X}}^{n}$
(2.5)

consisting of a natural number $n$ of input points
$x_{1},\dots,x_{n}\in\underline{\mathcal{X}}$ which are called the
_experimental points of $\tilde{x}$_. Also, $n$ is called the _size of the
experimental design $\tilde{x}$_. It should be noticed that the experimental
points of a given experimental design $\tilde{x}$ are not required to be
pairwise distinct. If they are, then $\tilde{x}$ is called an _experimental
design without replications (repetitions)_ , and otherwise $\tilde{x}$ is
called an _experimental design with replications (repetitions)_. Since
performing actual experiments in a lab is typically expensive, it is often
important to take into account experiments
$\tilde{x}^{-}=(x^{-}_{1},\dots,x^{-}_{n^{-}})$ that have already be performed
at an earlier stage and combine them with newly performed experiments
$\tilde{x}^{+}=(x^{+}_{1},\dots,x^{+}_{n^{+}})$. We will denote the resulting
_combined (or refined) experimental design_ by

$\displaystyle\tilde{x}^{-}\,\&\,\tilde{x}^{+}:=(x^{-}_{1},\dots,x^{-}_{n^{-}},x^{+}_{1},\dots,x^{+}_{n^{+}})\in\underline{\mathcal{X}}^{n^{-}+n^{+}}.$
(2.6)

In order to distinguish experimental designs in the sense of (2.5) from
weighted experimental designs defined next, we will sometimes refer to the
former as _unweighted experimental designs_.

If $\tilde{x}=(x_{1},\dots,x_{n})$ is a given unweighted experimental design,
it is often useful to aggregate the repeated experimental points, to indicate
the number of replications, and then to normalize by the size $n$ of the
design. In this manner, one can associate to every unweighted experimental
design $\tilde{x}=(x_{1},\dots,x_{n})$ the corresponding weighted experimental
design

$\displaystyle\xi_{\tilde{x}}:=\big{(}(w_{x},x))_{x\in\underline{\mathcal{X}}}\qquad\text{with}\qquad
w_{x}:=r_{x}/n,$ (2.7)

where $r_{x}\in\\{0,\dots,n\\}$ for every potential experiment
$x\in\underline{\mathcal{X}}$ indicates how often it occurs in the design
$\tilde{x}$. Sometimes, it is further useful to allow for general real weights
$w_{x}\in[0,1]$ instead of the specific rational weights in (2.7).
Specifically, for a finite design space $\underline{\mathcal{X}}$, it is
sometimes useful to consider tuples
$((w_{x},x))_{x\in\underline{\mathcal{X}}}$ of general weight-experiment
pairs. Such tuples are called _weighted experimental designs_ on
$\underline{\mathcal{X}}$ and are usually denoted by the letter $\xi$ in the
literature. Accordingly, the set of all weighted experimental designs on a
finite design space $\underline{\mathcal{X}}$ is denoted by

$\displaystyle\Xi(\underline{\mathcal{X}}):=\bigg{\\{}\big{(}(w_{x},x))_{x\in\underline{\mathcal{X}}}:w_{x}\in[0,1]\text{
for every }x\in\underline{\mathcal{X}}\text{ and
}\sum_{x\in\underline{\mathcal{X}}}w_{x}=1\bigg{\\}}.$ (2.8)

In parameter estimation and in experimental design alike, information matrices
play a fundamental role. As usual, the _one-point information matrix of $f$ at
a design point $x\in\underline{\mathcal{X}}$_ for the reference parameter
$\overline{\theta}$ is defined as

$\displaystyle
m_{f}(x,\overline{\theta}):=D_{\theta}f(x,\overline{\theta})^{\top}\varsigma^{-1}D_{\theta}f(x,\overline{\theta})\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}},$ (2.9)

where $\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$ denotes the set of all positive semidefinite
$d_{\mathrm{\theta}}\times d_{\mathrm{\theta}}$ matrices. Analogously, the
_information matrix of an experimental design_ – either an unweighted
experimental design $\tilde{x}=(x_{1},\dots,x_{n})$ or a weighted experimental
design
$\xi=((w_{x},x))_{x\in\underline{\mathcal{X}}}\in\Xi(\underline{\mathcal{X}})$
– is defined as the sum or, respectively, as the weighted sum of the
corresponding one-point information matrices. In formulas,

$\displaystyle
M_{f}(\tilde{x},\overline{\theta}):=\sum_{i=1}^{n}D_{\theta}f(x_{i},\overline{\theta})^{\top}\varsigma^{-1}D_{\theta}f(x_{i},\overline{\theta})\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}},$ (2.10) $\displaystyle
M_{f}(\xi,\overline{\theta}):=\sum_{x\in\underline{\mathcal{X}}}w_{x}D_{\theta}f(x,\overline{\theta})^{\top}\varsigma^{-1}D_{\theta}f(x,\overline{\theta})\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}.$ (2.11)

An extension of the tilde notation (2.5) for experimental designs to vector-
or matrix-valued functions on $\mathcal{X}$ or $\mathcal{X}\times\Theta$ is
very convenient both in parameter estimation and in experimental design.
Specifically, if $a:\mathcal{X}\to\mathcal{V}$ and
$b:\mathcal{X}\times\Theta\to\mathcal{V}$ are given functions from
$\mathcal{X}$ or $\mathcal{X}\times\Theta$ to some space $\mathcal{V}$ of real
column vectors or of real matrices and if $\tilde{x}=(x_{1},\dots,x_{n})$ is a
given experimental design of size $n$, then we will repeatedly use the
shorthand notations

$\displaystyle\tilde{a}(\tilde{x}):=\begin{pmatrix}a(x_{1})\\\ \vdots\\\
a(x_{n})\end{pmatrix}\in\mathcal{V}^{n}\qquad\text{and}\qquad\tilde{b}(\tilde{x},\theta):=\begin{pmatrix}b(x_{1},\theta)\\\
\vdots\\\ b(x_{n},\theta)\end{pmatrix}\in\mathcal{V}^{n}$ (2.12)

to denote the outcome of applying $a$ or $b$ simultaneously to all
experimental points of $\tilde{x}$. With this notation, Assumption 2 can be
concisely reformulated as follows: for every experimental design
$\tilde{x}=(x_{1},\dots,x_{n})$ on $\underline{\mathcal{X}}$, the
corresponding experimental results
$\tilde{y}=(y_{1}^{\top},\dots,y_{n}^{\top})^{\top}$ are realizations of the
$nd_{\mathrm{y}}$-dimensional random variable

$\displaystyle\tilde{y}(\tilde{x}):=\tilde{f}(\tilde{x},\theta^{*})+\tilde{\varepsilon}(\tilde{x}),$
(2.13)

where
$\tilde{\varepsilon}(\tilde{x})=(\varepsilon_{1}^{\top},\dots,\varepsilon_{n}^{\top})^{\top}$
is an $nd_{\mathrm{y}}$-dimensional normally distributed random variable with
expectation value $0$ and block-diagonal covariance matrix

$\displaystyle\tilde{\varsigma}:=\operatorname{diag}(\varsigma,\dots,\varsigma)=\begin{pmatrix}\varsigma&\dotsb&0\\\
\vdots&\ddots&\vdots\\\
0&\dotsb&\varsigma\end{pmatrix}\in\mathbb{R}_{\mathrm{pd}}^{nd_{\mathrm{y}}\times
nd_{\mathrm{y}}}.$ (2.14)

Similarly, applying the above tilde notation to the function
$J:\mathcal{X}\times\Theta\to\mathbb{R}^{d_{\mathrm{y}}\times
d_{\mathrm{\theta}}}$ defined by
$J(x,\overline{\theta}):=D_{\theta}f(x,\overline{\theta})$, the information
matrix of an experimental design $\tilde{x}=(x_{1},\dots,x_{n})$ can be
rewritten in the form

$\displaystyle
M_{f}(\tilde{x},\overline{\theta})=\tilde{J}(\tilde{x},\overline{\theta})^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x},\overline{\theta}).$
(2.15)

It is clear from this representation that the information matrix
$M_{f}(\tilde{x},\overline{\theta})$ is positive definite if and only if the
matrix
$\tilde{J}(\tilde{x},\overline{\theta})\in\mathbb{R}^{nd_{\mathrm{y}}\times
d_{\mathrm{\theta}}}$ has a trivial nullspace. In particular, for the
information matrix $M_{f}(\tilde{x},\overline{\theta})$ of
$\tilde{x}=(x_{1},\dots,x_{n})$ to be invertible, it is necessary that

$\displaystyle ld_{\mathrm{y}}\geq
d_{\mathrm{\theta}}\qquad(l:=\\#\\{x_{1},\dots,x_{n}\\})$ (2.16)

or, in other words, that the number $l$ of different experimental points in
$\tilde{x}$ is at least as large as $d_{\mathrm{\theta}}/d_{\mathrm{y}}$, the
number of model parameters divided by the number of outputs. As is well-known,
the invertibility of $M_{f}(\tilde{x},\overline{\theta})$ plays an essential
role in the estimability of the model parameters from the experimental design
$\tilde{x}$, see (2.18) and (2.20) below and [30] (Section 7.2).

###  2.2 Central facts from parameter estimation

Since the true value $\theta^{*}$ of the model parameter (Assumption 2) is
unknown, it needs to be estimated. A standard way of doing so is least-squares
estimation, that is, one measures the values $y_{1},\dots,y_{n}\in\mathcal{Y}$
of the output quantity at certain experimental points
$x_{1},\dots,x_{n}\in\underline{\mathcal{X}}$ and then chooses the model
parameters $\theta$ such that the model’s predictions
$f(x_{1},\theta),\dots,f(x_{n},\theta)$ deviate as little as possible from the
observed measurements $y_{1},\dots,y_{n}$. Specifically, one chooses the model
parameter $\theta\in\Theta$ as a solution to the least-squares estimation
problem

$\displaystyle\min_{\theta\in\Theta}\sum_{i=1}^{n}\big{(}f(x_{i},\theta)-y_{i}\big{)}^{\top}\varsigma^{-1}\big{(}f(x_{i},\theta)-y_{i}\big{)},$
(2.17)

where $\varsigma\in\mathbb{R}_{\mathrm{pd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$ is the covariance matrix of the measurement errors
(Assumption 2). Any solution to (2.17) is called a _least-squares estimate for
$f$ based on the experimental design $\tilde{x}=(x_{1},\dots,x_{n})$ with
corresponding experimental results
$\tilde{y}=(y_{1}^{\top},\dots,y_{n}^{\top})^{\top}$_. Since $\Theta$ is
assumed to be compact and $f(x_{i},\cdot)$ is continuous (Assumption 1), the
existence of a least-squares estimate is guaranteed for arbitrary experimental
data $\tilde{x}$ and $\tilde{y}$. If the model is further assumed to be
globally least-squares estimable for a given experimental design $\tilde{x}$,
then uniqueness of the least-squares estimate for $\tilde{x}$ can also be
established (Section 7.2 of [30]). We do not need such global estimability
assumptions, however, and will simply write
$\widehat{\theta}_{f}(\tilde{x},\tilde{y})$ for any of the potentially
multiple least-squares estimates for $f$ based on $\tilde{x}$ and $\tilde{y}$.

####  2.2.1 Computing least-squares estimates

In the special case of linear models, one has a closed-form expression for the
least-squares estimates. Specifically, for $f$ as in (2.2), the least-squares
estimate is unique and given by the formula

$\displaystyle\widehat{\theta}_{f}(\tilde{x},\tilde{y})=\big{(}\tilde{J}(\tilde{x})^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x})\big{)}^{-1}\tilde{J}(\tilde{x})^{\top}\tilde{\varsigma}^{-1}\big{(}\tilde{y}-\tilde{c}(\tilde{x})\big{)},$
(2.18)

provided that the occurring inverse of
$\tilde{J}(\tilde{x})^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x})$
exists, see (2.16). If it does not, then one can replace it in (2.18) by the
Moore-Penrose pseudoinverse to get the unique minimal-norm least-squares
estimate of $f$ for $\tilde{x}$ and $\tilde{y}$. As is well-known, the
pseudoinverse can be computed by means of the singular value decomposition or
by various regularizations. See [29] (Section 10.2) and [18] (Chapter 5) for
an overview of concrete methods to compute least-squares estimates for linear
models.

In the general case of nonlinear models, by contrast, one usually has no
closed-form expression for the least-squares estimates anymore and therefore
one has to resort to suitable approximations. A standard approximation is
based on linearizing $f$ around a suitable reference parameter value
$\overline{\theta}\in\Theta$. Specifically, in the least-squares problem
(2.17) one replaces the general nonlinear model $f$ by the linearized model
$f_{\overline{\theta}}^{\mathrm{lin}}$ defined by

$\displaystyle
f_{\overline{\theta}}^{\mathrm{lin}}(x,\theta):=f(x,\overline{\theta})+J(x,\overline{\theta})(\theta-\overline{\theta})\qquad((x,\theta)\in\mathcal{X}\times\mathbb{R}^{d_{\mathrm{\theta}}})$
(2.19)

with $J(x,\overline{\theta}):=D_{\theta}f(x,\overline{\theta})$. Applying then
(2.18), one obtains the approximate formula

$\displaystyle\widehat{\theta}_{f}(\tilde{x},\tilde{y})\approx\widehat{\theta}_{f_{\overline{\theta}}^{\mathrm{lin}}}(\tilde{x},\tilde{y})=\overline{\theta}+\big{(}\tilde{J}(\tilde{x},\overline{\theta})^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x},\overline{\theta})\big{)}^{-1}\tilde{J}(\tilde{x},\overline{\theta})^{\top}\tilde{\varsigma}^{-1}\big{(}\tilde{y}-\tilde{f}(\tilde{x},\overline{\theta})\big{)}$
(2.20)

provided that the occurring inverse of
$\tilde{J}(\tilde{x},\overline{\theta})^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x},\overline{\theta})=M_{f}(\tilde{x},\overline{\theta})$
exists, see (2.15) and (2.16) above. If one iterates this approximation (2.20)
starting from an initial $\overline{\theta}^{0}$ and if in each iteration one
chooses an appropriate step size, then one arrives at the Gauß-Newton method.
As is well-known, this iterative method converges to a critical point of the
sum-of-squares objective function in (2.17), provided that the matrices
$\tilde{J}(\tilde{x},\theta)^{\top}\tilde{\varsigma}^{-1}\tilde{J}(\tilde{x},\theta)$
are invertible for all $\theta$ in the region of interest (Theorem 10.1 of
[29]). If this uniform invertibility is not satisfied, one can use the
Levenberg-Marquardt method instead (Theorem 10.3 of [29]). We refer to [29]
(Section 10.3) and to [36] (Chapter 14) for a an extensive overview of methods
to compute least-squares estimates for nonlinear models.

####  2.2.2 Covariance of least-squares estimates

Since the experimental results $\tilde{y}$ for any given experimental design
$\tilde{x}$ are subject to random measurement errors, the same is true for the
corresponding least-squares parameter estimates
$\widehat{\theta}_{f}(\tilde{x},\tilde{y})$. In other words, the least-squares
estimates for any given experimental design $\tilde{x}$ are uncertain. A
standard measure to quantify this uncertainty is the covariance matrix

$\displaystyle\operatorname{Cov}\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$ (2.21)

of the least-squares estimate for $\tilde{x}$ and the random observations
$\tilde{y}(\tilde{x})$ from (2.13) (Assumption 2). In the case of linear
models, by (2.18) and (2.15) one has the exact identity

$\displaystyle\operatorname{Cov}\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))=M_{f}(\tilde{x},\overline{\theta})^{-1}$
(2.22)

for arbitrarily chosen reference parameter values
$\overline{\theta}\in\Theta$. Indeed, the information matrix
$M_{f}(\tilde{x},\overline{\theta})$ actually does not at all depend on
$\overline{\theta}$ for linear models (2.2). In the case of nonlinear models,
by (2.20) and (2.15) one has at least the approximate identity

$\displaystyle\operatorname{Cov}\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))\approx
M_{f}(\tilde{x},\overline{\theta})^{-1}$ (2.23)

for suitably chosen reference parameter values $\overline{\theta}\in\Theta$.
As opposed to linear models, the information matrix
$M_{f}(\tilde{x},\overline{\theta})$ of nonlinear models generally strongly
depends on the reference parameter value $\overline{\theta}$. In general, the
approximation (2.23) will be the better, the closer $\overline{\theta}$ is to
true parameter value $\theta^{*}$. A rigorous asymptotic justification of this
is provided by the classical asymptotic normality results for least-squares
estimation (Section 12.2.3 of [36] or Section 3.1.3 of [30]).

###  2.3 Central facts from experimental design

Since by (2.22) and (2.23) the uncertainty of the least-squares estimate
$\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))$ depends on the
underlying experimental design $\tilde{x}$, the prediction uncertainty of the
estimated model
$f(\cdot,\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x})))$ depends on the
experimental design as well. And so, in order to get a sufficiently precise
estimated model that is sufficiently certain about its predictions, one has to
carefully choose the experimental design $\tilde{x}$ on which the estimation
is based. In experimental design, one tries to come up with experimental
designs $\tilde{x}=(x_{1},\dots,x_{n})$ that on the one hand have a fesible
size $n$ and on the other hand lead to an estimated model that is sufficiently
certain about its predictions. Commonly, this uncertainty is quantified in
terms of a suitable scalar function

$\displaystyle\Psi(\operatorname{Cov}\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))^{-1})$
(2.24)

of the inverse of the estimate’s covariance. Such functions $\Psi$ are called
_design criteria_ and for the proposed methodology to work, we have to impose
only the following assumptions on $\Psi$ which are standard in convex design
theory [24, 39, 32, 2, 15, 30]. A reader not familiar with that theory can
simply always think of the particularly important special case where $\Psi$ is
the D- or the A-criterion, see (2.25) below.

######  Assumption 3.

$\Psi:\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}\to\mathbb{R}\cup\\{\infty\\}$ is convex, antitonic, lower
semicontinuos function whose restriction to
$\mathbb{R}_{\mathrm{pd}}^{d_{\mathrm{\theta}}\times d_{\mathrm{\theta}}}$
takes only finite values and is continuously differentiable.

All standard design criteria $\Psi_{p}$ except for the E-criterion
$\Psi_{\infty}$ satisfy the above assumptions and thus the methodology
proposed here can be applied with all standard design criteria except for the
E-criterion. In particular, the assumptions are satisfied for the D-criterion
$\Psi_{0}$ and for the A-criterion $\Psi_{1}$, which are defined as follows:

$\displaystyle\Psi_{0}(M):=\ln\det(M^{-1})\qquad\text{and}\qquad\Psi_{1}(M):=\operatorname{tr}(M^{-1})$
(2.25)

for $M\in\mathbb{R}_{\mathrm{pd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$ and $\Psi_{p}(M):=\infty$ for
$M\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}\setminus\mathbb{R}_{\mathrm{pd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$. See [35] for precise definitions of the notions used in
Assumption 3 and of the standard design criteria $\Psi_{p}$ as well as for
detailed verifications.

In the extensive literature on experimental design, two main approaches can be
distinguished, namely factorial design [16, 8, 26] on the one hand and optimal
design [23, 14, 24, 39, 32, 2, 15, 30] on the other hand. In essence, these
approaches differ in their primary goal and in their generality and scope of
applicability.

####  2.3.1 Factorial experimental design

In factorial experimental design, the primary goal is to find a factorial
experimental design $\tilde{x}^{*}$, that is, a geometrically regular design
whose experimental points constitute either a full grid in
$\underline{\mathcal{X}}$ or a regularly formed fraction of a full grid in
$\underline{\mathcal{X}}$. Accordingly, one speaks of full factorial designs
in the first case and of fractional factorial designs in the second case. In
contrast to optimal experimental design, the optimality of $\tilde{x}^{*}$
plays only a secondary role here. As it turns out, for certain fairly special
models $f$ and for cuboidal design spaces

$\displaystyle\underline{\mathcal{X}}=[a_{1},b_{1}]\times\dotsb\times[a_{d_{\mathrm{x}}},b_{d_{\mathrm{x}}}],$
(2.26)

suitably chosen factorial designs are also optimal (w.r.t. suitable design
criteria $\Psi$). Specifically, for linear models $f$ that are also linear or
at least multilinear in the inputs $x$, the standard full factorial designs
consisting of the corner points of (2.26) are A-, D-, and E-optimal (Section
3.3.4 of [4]). Similarly, for linear models $f$ that are quadratic in the
inputs $x$, fractional factorial designs can be found [25] that are D-optimal.
It should be noted, however, that these designs are weighted designs and their
support is fairly large even for small input dimensions (Section 2.2.2 of
[4]). It should also be noted that, unlike the specialized designs from [25],
the standard factorial designs for linear models that are quadratic in the
inputs are not D-optimal ([28] or Section 3.4.7 of [4]). As soon as the model
$f$ becomes more general (higher-order in $x$ or nonlinear in $\theta$),
factorial designs will no longer be optimal, however. In fact, this can
already be seen from third-order polynomial regression in one input dimension
(Section 2.2.1 of [4] or Section 2.6 of [15]) or in two input dimensions
(Section 2.2.2 of [4]). So, for applications described by higher-order linear
models or by nonlinear models, factorial design is not the method of choice
and optimal experimental design is recommended instead (Section 11.4.4 of
[26]). Similarly, optimal experimental design is recommended for models on
non-cuboidal input spaces $\mathcal{X}$ (Section 11.4.4 of [26]).

####  2.3.2 Optimal experimental design

In optimal experimental design, the primary goal is to find an optimal
experimental design $\tilde{x}^{*}$ w.r.t. the chosen design criterion $\Psi$,
that is, a solution to the optimization problem

$\displaystyle\min_{\tilde{x}\in\underline{\mathcal{X}}^{n}}\Psi(\operatorname{Cov}\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x}))^{-1}),$
(2.27)

where the design size $n$ can be fixed or allowed to vary in some given range.
In contrast to factorial design, the geometric regularity of $\tilde{x}^{*}$
plays no role here. And in fact, by the remarks in the previous paragraph,
optimal experimental designs are generally not factorial. As has been pointed
out above, except for linear models $f$, one usually has no exact expressions
for the estimate’s covariance and hence for the objective function in (2.27).
So, for nonlinear models $f$, one has to work with appropriate approximations
of the covariance. In many cases, the linearization-based approximation (2.23)
is used, which turns the original optimal experimental design problem (2.27)
into the substitute problem

$\displaystyle\min_{\tilde{x}\in\underline{\mathcal{X}}^{n}}\Psi(M_{f}(\tilde{x},\overline{\theta}))$
(2.28)

with some appropriate reference parameter value $\overline{\theta}\in\Theta$
and with a fixed or a variable design size $n$. It should be noticed that
(2.28) is a nonlinear optimization problem which is discrete, continuous, or
discrete-continuous depending on whether the design space
$\underline{\mathcal{X}}$ is discrete or continuous and on whether $n$ is
fixed or not. Conventionally, problems like (2.28) are called locally optimal
experimental design problems [9] and, accordingly, any approach based on
substituting (2.27) by (2.28) is referred to as locally optimal design. It is
important to notice that the linearization-based approximation (2.23) is by no
means restricted to locally optimal design. In fact, quite to the contrary,
the vast majority of approaches to optimal experimental design for nonlinear
models relies, in one way or another, on the approximation (2.23). In
particular, this is true for average optimal design and minimax optimal
design, which are defined by the problems

$\displaystyle\min_{\tilde{x}\in\underline{\mathcal{X}}^{n}}\int_{\Theta}\Psi(M_{f}(\tilde{x},\theta))\,\mathrm{d}\pi(\theta)\qquad\text{and}\qquad\min_{\tilde{x}\in\underline{\mathcal{X}}^{n}}\max_{\theta\in\Theta}\Psi(M_{f}(\tilde{x},\theta))$
(2.29)

with some prior probability distribution $\pi$ on $\Theta$. Additionally, this
is true for the class of sequential design methods [10, 39, 15, 30], to which
the method proposed here belongs. We refer to [15, 30] for comprehensive
expositions of these linearization-based approaches to optimal experimental
design for nonlinear models and to [17] for a concise overview of locally
optimal experimental design specifically tailored to dynamic nonlinear models.
We also refer to [19, 43, 27] for approaches that – unlike the first-order
approaches mentioned so far – propose higher-order approximations of the
covariance or of the related confidence regions of the least-squares estimate.

##  3 Sequential locally optimal experimental design

In this section, we formally introduce our sequential locally optimal
experimental design method. In rough terms, the method repeatedly solves
locally optimal experimental design problems of the kind

$\displaystyle\min_{\tilde{x}\in\underline{\mathcal{X}}^{n}}\Psi(M_{f}(\tilde{x},\overline{\theta}))$
(3.1)

with successively improved reference parameter values $\overline{\theta}$.
Specifically, in every iteration, $\overline{\theta}$ is chosen to be the
current best parameter estimate, based on all experiments performed until
then, and the locally optimal experimental design problem with that value of
$\overline{\theta}$ is solved by means of the recent adaptive discretization
algorithm from [48]. We continue this iteration until the total number of
performed experiments reaches some upper bound or the newly proposed
experiments do not differ substantially anymore from already performed
experiments. As is usual [15, 30, 48] for concrete experimental design
algorithms, we choose the design space $\underline{\mathcal{X}}$ to be a
finite set and work with weighted designs for the sake of computational
feasibility.

######  Assumption 4.

$\underline{\mathcal{X}}$ is a finite set and for every given reference
parameter $\overline{\theta}\in\Theta$ all the one-point information matrices
$m_{f}(x,\overline{\theta})=D_{\theta}f(x,\overline{\theta})^{\top}\varsigma^{-1}D_{\theta}f(x,\overline{\theta})$
for $x\in\underline{\mathcal{X}}$ can be computed sufficiently fast.

In the case of an explicit model $f$, the computability assumption above poses
no restriction. In the case of an implicit model $f$, in turn, it boils down
to requiring a sufficiently fast numerical solvability of the equation (2.3.a)
for $s=s(x,\overline{\theta})$ along with a sufficiently fast numerical
invertibility of the Jacobi matrices
$D_{s}g(s(x,\overline{\theta}),x,\overline{\theta})$ for all
$x\in\underline{\mathcal{X}}$.

###  3.1 Computing locally optimal experimental designs

At the heart of our sequential experimental design method is the solution of
suitable variants of the locally optimal experimental design problem (3.1). In
essence, these variants are obtained by means of two modifications.

As a first modification to (3.1), we extend the search space from the
unweighted experimental designs of a given size $n$ to all weighted
experimental designs on $\underline{\mathcal{X}}$. In this manner, we arrive
at the one-stage locally optimal experimental design problem

$\displaystyle\min_{\xi\in\Xi(\underline{\mathcal{X}})}\Psi\big{(}M_{f}(\xi,\overline{\theta})\big{)}.$
(3.2)

Since the weighted designs $\xi\in\Xi(\underline{\mathcal{X}})$ are uniquely
determined by their weights $w_{x}$ for $x\in\underline{\mathcal{X}}$, the
design problem (3.2) is actually a weight optimization problem. Since,
moreover, the design criterion $\Psi$ is convex (Assumption 3), the weight
optimization problem (3.2) is a convex nonlinear optimization problem and
therefore has a much simpler structure than the original general nonlinear
design problem (3.1). It is precisely this substantial structural
simplification that makes the introduction of weighted experimental designs
standard in optimal experimental design. In the literature, problems like
(3.2) are referred to as one-stage experimental design problems because their
design criterion does not explicitly take into account experiments that have
already been performed at an earlier stage. Such preliminary experiments,
however, are quite often available in experimental practice.

As a further modification to (3.1), we therefore explicitly incorporate such
previous experiments into the design criterion. Specifically, we consider the
two-stage locally optimal experimental design problem

$\displaystyle\min_{\xi\in\Xi(\underline{\mathcal{X}})}\Psi\big{(}\alpha
M_{f}(\xi_{\tilde{x}^{-}},\overline{\theta})+(1-\alpha)M_{f}(\xi,\overline{\theta})\big{)},$
(3.3)

where $\tilde{x}^{-}=(x_{1}^{-},\dots,x_{n^{-}}^{-})$ are experiments that
have already been performed at a previous stage, $\xi_{\tilde{x}^{-}}$ is the
corresponding weighted design (2.7), and $\alpha\in[0,1]$ is an importance
factor for those previously performed experiments. In the limiting case
$\alpha=0$, the information contained in the already performed experiments is
completely ignored in the design of the new experiments, and thus (3.3)
reduces to (3.2) then. In the limiting case $\alpha=1$, by contrast, the
information contained in the new experiments is completely ignored which makes
the optimization (3.3) pointless, of course. If one chooses $\alpha=0.5$, then
the uncertainty induced by the previous experiments and the uncertainty
induced by the prospective experiments are both taken into account in equal
measure. Just like (3.2), the two-stage experimental design problem (3.3) is
actually a weight optimization problem and, by the convexity of $\Psi$
(Assumption 3), it is a convex continuous optimization problem.

In order to solve one- or two-stage locally optimal design problems of the
kind (3.2) or (3.3), one can resort to a large variety of tailor-made
iterative algorithms from the literature like the algorithms from [14, 47, 3,
38, 7, 49] or [20], for instance. See also [2, 15, 30] for a good overview of
the pertinent algorithms known until 2013. In this paper, we use the more
recent adaptive discretization algorithm from [48] because of its very good
convergence properties, both from a theoretical and from a practical point of
view [48, 35, 37, 42]. Instead of solving (3.3) in one shot, this algorithm
repeatedly solves discretized versions of (3.3) with successively refined
discretizations. And since these discretizations are refined adaptively, the
discretized problems to be solved in the course of the algorithm remain fairly
small, typically. Specifically, we use the following version of the algorithm
[35].

######  Algorithm 1.

Input: an unweighted design $\tilde{x}^{-}$ on $\underline{\mathcal{X}}$ with
an associated importance factor $\alpha\in[0,1)$, a reference parameter value
$\overline{\theta}\in\Theta$, the one-point information matrices
$m_{f}(x,\overline{\theta})$ for all $x\in\underline{\mathcal{X}}$, a small
finite subset $\underline{\mathcal{X}}^{0}$ of $\underline{\mathcal{X}}$ such
that $\alpha
M_{f}(\xi_{\tilde{x}^{-}},\overline{\theta})+(1-\alpha)M_{f}(\xi,\overline{\theta})$
is invertible for some design $\xi$ on $\underline{\mathcal{X}}^{0}$, and an
optimality tolerance $\varepsilon\in(0,\infty)$. Initialize $k=0$. With these
inputs, proceed in the following steps.

  * 1.

Compute a solution $\xi^{k}\in\Xi(\underline{\mathcal{X}}^{k})$ of the
discretized version

$\displaystyle\min_{\xi\in\Xi(\underline{\mathcal{X}}^{k})}\Psi^{(\alpha)}(M_{f}(\xi,\overline{\theta}))$
(3.4)

of the two-stage locally optimal design problem (3.3), where $\Psi^{(\alpha)}$
is the two-stage design criterion defined by $\Psi^{(\alpha)}(M):=\Psi(\alpha
M_{f}(\xi_{\tilde{x}^{-}},\overline{\theta})+(1-\alpha)M)$ for all
$M\in\mathbb{R}_{\mathrm{psd}}^{d_{\mathrm{\theta}}\times
d_{\mathrm{\theta}}}$.

  * 2.

Check how far the computed solution $\xi^{k}$ of the discretized problem (3.4)
is from being optimal for the original problem (3.3). In order to do so,
compute a solution $x^{k}\in\underline{\mathcal{X}}$ of the strongest
optimality violator problem

$\displaystyle\min_{x\in\underline{\mathcal{X}}}\psi^{(\alpha)}(M_{f}(\xi^{k},\overline{\theta}),x),$
(3.5)

where $\psi^{(\alpha)}$ is the sensitivity function of $\Psi^{(\alpha)}$
w.r.t. $M_{f}(\Xi(\underline{\mathcal{X}}),\overline{\theta})$ (which is a
certain directional derivative of $\Psi^{(\alpha)}$).

    * –

In case
$\psi^{(\alpha)}(M_{f}(\xi^{k},\overline{\theta}),x^{k})<-\varepsilon$, refine
the discretization by setting
$\underline{\mathcal{X}}^{k+1}:=\underline{\mathcal{X}}^{k}\cup\\{x^{k}\\}$
and return to Step 1 with $k$ replaced by $k+1$.

    * –

In the opposite case, return $\xi^{k}$ and terminate.

As is shown in [35], the algorithm above terminates after finitely many
iterations with an approximately optimal design
$\xi^{+}\in\Xi(\underline{\mathcal{X}})$ in the sense that

$\displaystyle\Psi\big{(}\alpha
M_{f}(\xi_{\tilde{x}^{-}},\overline{\theta})+(1-\alpha)M_{f}(\xi^{+},\overline{\theta})\big{)}\leq\min_{\xi\in\Xi(\underline{\mathcal{X}})}\Psi\big{(}\alpha
M_{f}(\xi_{\tilde{x}^{-}},\overline{\theta})+(1-\alpha)M_{f}(\xi,\overline{\theta})\big{)}+\varepsilon.$
(3.6)

In words, this relation says that the design criterion value for $\xi^{+}$ is
guaranteed to exceed the optimal design criterion value by at most
$\varepsilon$. It is important to notice that this approximation error bound
$\varepsilon$ can be freely set by the user. In the following, we will refer
to a design $\xi^{+}\in\Xi(\underline{\mathcal{X}})$ that satisfies (3.6) as a
_locally $\varepsilon$-optimal weighted design for $f$ on
$\underline{\mathcal{X}}$ around $\overline{\theta}$ with $\tilde{x}^{-}$ as
the previous-stage design of importance $\alpha$_. Algorithm 1 is an adaptive
discretization algorithm because the discretization
$\underline{\mathcal{X}}^{k+1}$ is obtained from the previous discretization
$\underline{\mathcal{X}}^{k}$ by an adaptive refinement, namely by adding a
strongest violator $x^{k}$ of the optimality condition for the solution
$\xi^{k}$ of the previous discretized design problem (3.4). It is because of
this adaptive refinement that Algorithm 1 typically does not need many
iterations until termination. In particular, the discretized design problems
(3.4) to be solved in the course of the algorithm are typically fairly small
and hence easy to solve. After all, these discretized design problems (3.4)
are convex weight optimization problems with the number of weight variables
being equal to the size of the discretization $\underline{\mathcal{X}}^{k}$.
Also, by the finiteness of $\underline{\mathcal{X}}$ (Assumption 4), the
strongest violator problems (3.5) can simply be solved by means of
enumeration, using the explicit formulas for the sensitivity function
$\psi^{(\alpha)}$ of two-stage design criteria $\Psi^{(\alpha)}$ provided in
[35], for instance.

As is indicated already by their name, the locally $\varepsilon$-optimal
weighted designs computed by Algorithm 1 are weighted experimental designs
$\xi^{+}=((w_{x}^{+},x))_{x\in\underline{\mathcal{X}}}\in\Xi(\underline{\mathcal{X}})$.
And in most cases, the weights $w_{x}^{+}$ will not be of the form
$w_{x}^{+}=r_{x}^{+}/n^{+}$ with a practically realizable number $n^{+}$ of
new experiments and replication numbers $r_{x}^{+}\in\\{0,\dots,n^{+}\\}$
indicating how often the experimenter should measure at the points
$x\in\underline{\mathcal{X}}$. In other words, the weighted experimental
designs $\xi^{+}$ computed by Algorithm 1 usually do not directly correspond
to an unweighted experimental design
$\tilde{x}^{+}=(x_{1}^{+},\dots,x_{n^{+}}^{+})$ of a practically realizable
size $n^{+}$. Instead, they first have to be converted to such a practically
realizable unweighted design. Standard methods of doing so are based on
suitable rounding procedures for the (pseudo)quotas $n^{+}w_{x}^{+}$ [15]
(Section 3.3), $(n^{+}-l^{+})w_{x}^{+}$ [14] (Section 3.1), or
$(n^{+}-l^{+}/2)w_{x}^{+}$ [33] (Section 6), where $l^{+}$ is the number of
experimental points in $\xi^{+}$. In the case where the affordable number
$n^{+}$ of new experiments is large, these rounding procedures produce
unweighted designs $\tilde{x}^{+}$ that are efficient in the sense that the
design criterion value of $\xi_{\tilde{x}^{+}}$ exceeds that of $\xi^{+}$ only
by a little bit. Indeed, for large experimental budgets $n^{+}$, the
efficiency bounds from [14, 15, 33] become tight. In the practically relevant
case of small experimental budgets $n^{+}$, however, the aforementioned
efficiency bounds do not say much about the actually achievable efficiency.
Additionally, the mentioned efficiency bounds are valid only for homogeneous
design criteria [15], but two-stage design criteria – as the criteria
$\Psi^{(\alpha)}$ used here – are not homogeneous. We therefore propose a
different conversion procedure in this paper, namely one that is not based on
rounding but instead directly aims at producing efficient unweighted designs
$\tilde{x}^{+}$ from $\xi^{+}$. Specifically, among the experimental points of
$\xi^{+}$ with non-negligible weights, we simply select $n^{+}$ points such
that the design criterion value becomes as small as possible.

Summarizing, we arrive at the following two-stage locally optimal experimental
design algorithm. Starting from an unweighted design $\tilde{x}^{-}$
containing experiments from a previous stage, it computes an unweighted
experimental design $\tilde{x}^{+}$ containing the experiments for the next
stage.

######  Algorithm 2.

Input: an unweighted experimental design $\tilde{x}^{-}$ on
$\underline{\mathcal{X}}$ with an associated importance $\alpha\in[0,1)$, the
one-point information matrices $m_{f}(x,\overline{\theta})$ for all
$x\in\underline{\mathcal{X}}$, a reference parameter value
$\overline{\theta}\in\Theta$, an optimality tolerance
$\varepsilon\in(0,\infty)$, a minimal total weight
$\underline{w}^{+}\in(0,1]$, and a maximal number $\overline{n}^{+}$ of new
experiments. With these inputs, proceed in the following steps.

  * 1.

Compute a locally $\varepsilon$-optimal weighted design

$\displaystyle\xi^{+}=((w_{x}^{+},x))_{x\in\underline{\mathcal{X}}}$ (3.7)

for $f$ on $\underline{\mathcal{X}}$ around $\overline{\theta}$ with
$\tilde{x}^{-}$ as the previous-stage design of importance $\alpha$. In order
to do so, use Algorithm 1.

  * 2.

Convert the computed weighted design $\xi^{+}$ to an unweighted design
$\tilde{x}^{+}$ consisting of at most $\overline{n}^{+}$ experiments. In order
to do so, successively sieve out the least important experiment from $\xi^{+}$
until the total weight of the remaining experiments only just remains above
the prescribed threshold $\underline{w}^{+}$. Write

$\displaystyle\tilde{\underline{x}}^{+}=(\underline{x}_{1}^{+},\dots,\underline{x}_{\underline{n}^{+}}^{+})$
(3.8)

for the remaining unweighted design and then compare the number
$\underline{n}^{+}$ of experiments in $\tilde{\underline{x}}^{+}$ with the
prescribed upper bound $\overline{n}^{+}$.

    * –

In case $\underline{n}^{+}>\overline{n}^{+}$, choose an appropriate subdesign
$\tilde{x}^{+}$ of the intermediate design $\tilde{\underline{x}}^{+}$,
consisting of exactly $\overline{n}^{+}$ experiments. Specifically, among all
subdesigns $\tilde{z}^{+}$ of $\tilde{\underline{x}}^{+}$ consisting of
exactly $\overline{n}^{+}$ experiments, choose one with minimal design
criterion value

$\displaystyle\Psi\big{(}\alpha
M_{f}(\xi^{-},\overline{\theta})+(1-\alpha)M_{f}(\xi_{\tilde{z}^{+}},\overline{\theta})\big{)}$
(3.9)

    * –

In the opposite case, simply choose $\tilde{x}^{+}$ to be equal to the
intermediate design $\tilde{\underline{x}}^{+}$.

Clearly, the output $\tilde{x}^{+}$ of the above algorithm is an unweighted
experimental design without replications consisting of at most
$\overline{n}^{+}$ experiments in $\underline{\mathcal{X}}$. In the following,
we will refer to it as a _locally $\varepsilon$-optimal design for $f$ on
$\underline{\mathcal{X}}$ around $\overline{\theta}$ consisting of at most
$\overline{n}^{+}$ experiments with total weight at least $\underline{w}^{+}$
and with $\tilde{x}^{-}$ as the previous-stage design of importance $\alpha$_.

###  3.2 Iterating experiments, parameter estimation, and locally optimal
experimental design

As has been pointed out above, the quality of the approximation (2.23)
generally strongly depends on the specific choice of the reference parameter
$\overline{\theta}$ about which the model is linearized. So, as the locally
optimal experimental design problems (3.2) and (3.3) are built precisely on
this approximation, it is essential to come up with good choices for the
reference parameter $\overline{\theta}$. A very natural and established
strategy here is sequential design [10, 39, 15, 30] where $\overline{\theta}$
is simply always chosen to be the current best parameter estimate
$\widehat{\theta}_{f}(\tilde{x}^{-},\tilde{y}^{-})$, based on all experimental
data $\tilde{x}^{-}$ and $\tilde{y}^{-}$ that are available at the current
stage. With this estimate $\overline{\theta}$, we then apply our locally
optimal experimental design algorithm (Algorithm 2) to get the most
informative batch $\tilde{x}^{+}$ of new experiments. We repeat this three-
step procedure of performing experiments, estimating the parameters, and
computing locally optimal experimental designs until the total number of
experiments exceeds a user-specified upper bound $\overline{n}$ or each of the
newly proposed experiments differs from an already performed experiment by
less than a user-specified progress tolerance $\delta$. We measure this
distance between two experimental points
$x=(x_{j})_{j\in\\{1,\dots,d_{\mathrm{x}}\\}}$ and
$x^{\prime}=(x^{\prime}_{j})_{j\in\\{1,\dots,d_{\mathrm{x}}\\}}$ in terms of a
suitably scaled maximum norm defined by

$\displaystyle\left\|x-x^{\prime}\right\|_{\infty}:=\max_{j\in\\{1,\dots,d_{\mathrm{x}}\\}}|x_{j}-x^{\prime}_{j}|/\lambda_{j},$
(3.10)

where $\lambda_{1},\dots,\lambda_{d_{\mathrm{x}}}$ are the side lengths of the
smallest cuboid (axis-aligned box) containing the design space
$\underline{\mathcal{X}}$. (In other words,
$\lambda_{1},\dots,\lambda_{d_{\mathrm{x}}}$ describe the extension of the
design space in the individual input directions.) In detail, our sequential
design procedure is described in Algorithm 3. A schematic illustration of the
procedure is given in Figure 1.

Figure 1: Schematic of the proposed sequential locally optimal experimental
design method.

######  Algorithm 3.

Input: an initial experimental design $\tilde{x}^{0}$ on
$\underline{\mathcal{X}}$, an importance factor $\alpha\in[0,1)$ for existing
experiments, an optimality tolerance $\varepsilon\in(0,\infty)$, a minimal
total weight $\underline{w}^{+}\in(0,1]$, a maximal number $\overline{n}^{+}$
of new experiments per iteration, a maximal total number $\overline{n}$ of
experiments, and a progress tolerance $\delta\in(0,\infty)$. Initialize $k=0$.
With these inputs, proceed in the following steps.

  * 1.

Carry out the experiments that have been newly proposed in the previous
iteration, that is, measure the output quantity at those experimental points
of $\tilde{x}^{k}=(x_{1}^{k},\dots,x_{n_{k}}^{k})$ that have been newly
proposed. Write

$\displaystyle\tilde{y}^{k}:=((y_{1}^{k})^{\top},\dots,(y_{n_{k}}^{k})^{\top})^{\top}.$
(3.11)

for all experimental results collected so far (the newly and the previously
collected ones).

  * 2.

Compute a least-squares estimate
$\widehat{\theta}^{k}:=\widehat{\theta}_{f}(\tilde{x}^{k},\tilde{y}^{k})$ for
$f$ based on all experimental data $(\tilde{x}^{k},\tilde{y}^{k})$ collected
so far, and then compute the one-point information matrices
$m_{f}(x,\widehat{\theta}^{k})$ for all candidate design points
$x\in\underline{\mathcal{X}}$.

  * 3.

Compute a locally $\varepsilon$-optimal design $\tilde{x}^{k,+}$ for $f$ on
$\underline{\mathcal{X}}$ around $\widehat{\theta}^{k}$ consisting of at most
$\overline{n}^{+}$ experiments with total weight at least $\underline{w}^{+}$
and with the current design $\tilde{x}^{k}$ as the previous-stage design of
importance $\alpha$. In order to do so, use Algorithm 2. Write

$\displaystyle\tilde{x}^{k+1}:=\tilde{x}^{k}\,\&\,\tilde{x}^{k,+}$ (3.12)

for the combined design consisting of the current experiments and the newly
computed experiments, $n_{k+1}:=n_{k}+n_{k}^{+}$ for the size of
$\tilde{x}^{k+1}$, and
$\operatorname{dist}(x,\tilde{x}^{k}):=\min_{i\in\\{1,\dots,n_{k}\\}}\left\|x-x_{i}\right\|_{\infty}$
for the distance of any given $x$ from the experimental points of
$\tilde{x}^{k}$.

    * –

In case $n_{k+1}>\overline{n}$ or
$\operatorname{dist}(x,\tilde{x}^{k})<\delta$ for every experimental point $x$
of $\tilde{x}^{k,+}$, terminate.

    * –

In the opposite case, return to Step 1 with $k$ replaced by $k+1$.

Clearly, if the progress tolerance $\delta$ is chosen small enough, then the
progress-based termination criterion of the above algorithm reduces to the
condition that each of the newly computed experimental points of
$\tilde{x}^{k,+}$ is actually equal to one of the experimental points of the
previous design $\tilde{x}^{k}$. Specifically, this is the case if $\delta$ is
chosen to be smaller than the mesh size (granularity)

$\displaystyle\min\big{\\{}\left\|x-x^{\prime}\right\|_{\infty}:x,x^{\prime}\in\underline{\mathcal{X}}\text{
with }x\neq x^{\prime}\big{\\}}$ (3.13)

of the chosen finite design space $\underline{\mathcal{X}}$ (Assumption 4).

###  3.3 Assessing the quality of experimental designs

In this section, we introduce the criteria that will be used to assess the
quality of a given experimental design. We use three different criteria, and
all of them express the quality of an experimental design
$\tilde{x}=(x_{1},\dots,x_{n})$ in terms of the prediction quality of the
corresponding estimated model
$f(\cdot,\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x})))$, that is, the
model that has been trained based on $\tilde{x}$.

As a first quality measure of a given experimental design $\tilde{x}$ (with
associated experimental results $\tilde{y}$), we take the sum of squared
prediction errors that the trained model
$f(\cdot,\widehat{\theta}_{f}(\tilde{x},\tilde{y}))$ makes on all available
experimental data $(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$.
Specifically, we take the root mean squared errors

$\displaystyle\rho_{j}(\tilde{x}):=\rho_{j}(\tilde{x},\tilde{y}):=\bigg{(}\frac{1}{n^{\mathrm{tot}}}\sum_{i=1}^{n^{\mathrm{tot}}}\big{(}f_{j}(x^{\mathrm{tot}}_{i},\widehat{\theta}_{f}(\tilde{x},\tilde{y}))-y^{\mathrm{tot}}_{ij}\big{)}^{2}\bigg{)}^{1/2}\qquad(j\in\\{1,\dots,d_{\mathrm{y}}\\})$
(3.14)

of the trained model’s components on the totality of available data. As our
second and third quality measures for a given experimental design $\tilde{x}$,
we take two different approximations of the prediction variance

$\displaystyle
V_{j}(x,\tilde{x}):=\operatorname{Var}f_{j}(x,\widehat{\theta}_{f}(\tilde{x},\tilde{y}(\tilde{x})))$
(3.15)

of the estimated model’s components based on $\tilde{x}$, namely a
linearization-based approximation $V_{j}^{\mathrm{lin}}$ on the one hand and a
sampling-based approximation $V_{j}^{\mathrm{sam}}$ on the other hand.
Specifically, the linearization-based prediction variance
$V_{j}^{\mathrm{lin}}(x,\tilde{x})$ is defined as the variance of the
estimated linearized model at $x$, that is,

$\displaystyle
V_{j}^{\mathrm{lin}}(x,\tilde{x}):=\operatorname{Var}f^{\mathrm{lin}}_{\overline{\theta},j}(x,\widehat{\theta}_{f^{\mathrm{lin}}_{\overline{\theta}}}(\tilde{x},\tilde{y}(\tilde{x})))\qquad(j\in\\{1,\dots,d_{\mathrm{y}}\\}),$
(3.16)

where $x\in\mathcal{X}$ is an arbitrary input point,
$\tilde{y}(\tilde{x}):=\tilde{f}(\tilde{x},\theta^{*})+\tilde{\varepsilon}(\tilde{x})$
is the random observation model (2.13), and the model is linearized around
$\overline{\theta}:=\widehat{\theta}_{f}(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$,
which is the best parameter estimate based on the totality of available data.
In view of (2.19) and (2.22),

$\displaystyle
V_{j}^{\mathrm{lin}}(x,\tilde{x})=\nabla_{\theta}f_{j}(x,\overline{\theta})^{\top}M_{f}(\tilde{x},\overline{\theta})^{-1}\nabla_{\theta}f_{j}(x,\overline{\theta}).$
(3.17)

Additionally, the sampling-based prediction variance
$V_{j}^{\mathrm{sam}}(x,\tilde{x})$ is defined as the sample variance of the
estimated nonlinear model at $x$, that is,

$\displaystyle
V_{j}^{\mathrm{sam}}(x,\tilde{x}):=\frac{1}{n^{\mathrm{sam}}}\sum_{s=1}^{n^{\mathrm{sam}}}\big{(}f_{j}(x,\widehat{\theta}_{f}(\tilde{x},\tilde{y}^{(s)}))-\mu_{j}^{\mathrm{sam}}(x,\tilde{x})\big{)}^{2}\qquad(j\in\\{1,\dots,d_{\mathrm{y}}\\}),$
(3.18)
$\displaystyle\mu_{j}^{\mathrm{sam}}(x,\tilde{x}):=\frac{1}{n^{\mathrm{sam}}}\sum_{s=1}^{n^{\mathrm{sam}}}f_{j}(x,\widehat{\theta}_{f}(\tilde{x},\tilde{y}^{(s)})),$
(3.19)

where $x\in\mathcal{X}$ is an arbitrary input point, $\tilde{y}^{(s)}$ are
$n^{\mathrm{sam}}$ samples independently drawn from the random observation
model
$\tilde{y}(\tilde{x},\overline{\theta}):=\tilde{f}(\tilde{x},\overline{\theta})+\tilde{\varepsilon}(\tilde{x})$
with
$\overline{\theta}:=\widehat{\theta}_{f}(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$.
Since this $\overline{\theta}$ is the best parameter estimate based on the
totality of available data, the employed observation model
$\tilde{y}(\tilde{x},\overline{\theta})$ is the best available substitute for
the true observation model (2.13). Similarly to (3.14), we will use the square
roots

$\displaystyle\sigma_{j}^{\mathrm{lin}}(x,\tilde{x}):=\big{(}V_{j}^{\mathrm{lin}}(x,\tilde{x})\big{)}^{1/2}\qquad\text{and}\qquad\sigma_{j}^{\mathrm{sam}}(x,\tilde{x}):=\big{(}V_{j}^{\mathrm{sam}}(x,\tilde{x})\big{)}^{1/2}$
(3.20)

of the prediction variances (3.16) and (3.18) at given input points
$x\in\mathcal{X}$. In fact, in in our comparisons of different designs, we
will actually use the $x$-independent worst-case uncertainties

$\displaystyle\sigma_{j}^{\mathrm{lin}}(\tilde{x}):=\max_{x\in\mathcal{X}}\sigma_{j}^{\mathrm{lin}}(x,\tilde{x})\qquad\text{and}\qquad\sigma_{j}^{\mathrm{sam}}(\tilde{x}):=\max_{x\in\mathcal{X}}\sigma_{j}^{\mathrm{sam}}(x,\tilde{x}).$
(3.21)

We will refer to (3.21) as the worst-case linearization- and the worst-case
sampling-based prediction uncertainties, respectively.

It is clear that for substantially nonlinear models $f$, the sampling-based
prediction uncertainty $\sigma_{j}^{\mathrm{sam}}(x,\tilde{x})$ is a better
approximation to the exact prediction uncertainty than the linearization-based
prediction uncertainty $\sigma_{j}^{\mathrm{lin}}(x,\tilde{x})$, just because
the former does not rely on any linearization. It is also clear by the very
definitions (3.14) and (3.20) of our three quality measures for experimental
designs that they all depend not only on the experimental design $\tilde{x}$
alone but also on the total set of gathered experimental data
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$. It should be noticed
that this dependence on $(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$
is strongest for the root mean squared errors (3.14), while it is only
indirect for the linearization-based prediction uncertainty (3.20.a) and even
more so for the sampling-based prediction unctertainty (3.20.b). In order to
keep the notation simple, we suppress this additional dependence on
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$ in the following.

##  4 Application to vapor-liquid equilibrium modeling

In this section, we apply the proposed general sequential design method to the
estimation of non-random two-liquid parameters. Specifically, we consider the
binary system consisting of propanol and propyl acetate. It is known that this
system exhibits narrow azeotropic behavior and therefore the estimation of its
non-random two-liquid parameters is already a challenging task. With our
sequential locally optimal design method, we end up with a locally optimal
design consisting of $15$ experiments, four of which are repeated experiments.
As we will see, traditional factorial design, by contrast, already requires
$27$ experiments to achieve roughly the same prediction quality.

###  4.1 Considered system

As explained above, the system considered here is the narrow azeotropic binary
mixture consisting of propanol and propyl acetate. In the following, propanol
will always be considered as component $1$ and propyl acetate as component
$2$. In order to estimate the binary interaction parameters of this two-
component mixture, we observe it at different liquid compositions $(l,1-l)$
and at different pressures $P$ and, for all these conditions, we measure the
corresponding vapor composition $(v,1-v)$ and the temperature $T$ at vapor-
liquid equilibrium. In our general terminology, this means that the input
quantities $x$ are the liquid mole fraction $l$ of propanol along with the
system pressure $P$, while the output quantities $y$ are the vapor mole
fraction $v$ of propanol along with the system temperature $T$ at vapor-liquid
equilibrium. In short,

$\displaystyle
x:=(l,P)\in\mathcal{X}:=[0,1]\times[\underline{P},\overline{P}]\qquad\text{and}\qquad
y:=(v,T)\in\mathcal{Y}:=[0,1]\times[\underline{T},\overline{T}],$ (4.1)

where $\underline{P},\overline{P}$ and $\underline{T},\overline{T}$ determine
the pressure and temperature ranges of interest. In our case, pressure and
temperature are consistently measured, respectively, in units of Pascal (Pa)
and Kelvin (K), and mole fractions are measured in the unit mol/mol, of
course. Also, $\underline{P}:=1\cdot 10^{5}$ and $\overline{P}:=3\cdot 10^{5}$
in our case.

In order to model the functional relationship between the input and output
quantities (4.1), we assume Raoult’s law in its extended version without
Poynting correction. In particular, we assume non-ideal mixture behavior in
the liquid phase but an ideal vapor phase (Section 13.3 of [40]). Spelled out,
this means that for every given input condition $x=(l,P)\in\mathcal{X}$, the
equilibrium temperature $T$ is implicitly defined by the equation

$\displaystyle\frac{P_{1}(T)}{P}\cdot\gamma_{1}(l,T)l+\frac{P_{2}(T)}{P}\cdot\gamma_{2}(l,T)(1-l)=1,$
(4.2)

and the equilibrium vapor mole fraction $v$ of propanol, in turn, is defined
by the equation

$\displaystyle v=\frac{P_{1}(T)}{P}\cdot\gamma_{1}(l,T)l.$ (4.3)

In the above equations, $P_{i}(T)$ denotes the saturation vapor pressure of
pure component $i$ at temperature $T$ and $\gamma_{i}(l,T)$ denotes the
activity coefficient of component $i$ in the mixture with liquid composition
$(l,1-l)$ at temperature $T$. Concerning the saturation vapor pressure, we
assume it to be given by the Antoine model (Section 6.5 of [40]), that is,

$\displaystyle P_{i}(T):=10^{5}\cdot
10^{A_{i}-B_{i}/(T+C_{i})}\qquad(i\in\\{1,2\\}),$ (4.4)

and for the Antoine parameters $A_{i},B_{i},C_{i}$ we take the estimated
values from Table 1. Concerning the activity coefficients, we assume them to
be given by the non-random two-liquid model of Renon and Prausnitz (Section
13.4 of [40] or Section 6.15 of [31]), that is,

$\displaystyle\gamma_{i}(l,T):=\exp\bigg{(}\frac{\sum_{k=1}^{2}l_{k}\tau_{ki}(T)G_{ki}(T)}{\sum_{k=1}^{2}l_{k}G_{ki}(T)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\sum_{j=1}^{2}\frac{l_{j}G_{ij}(T)}{\sum_{k=1}^{2}l_{k}G_{kj}(T)}\bigg{(}\tau_{ij}(T)-\frac{\sum_{k=1}^{2}l_{k}\tau_{kj}(T)G_{kj}(T)}{\sum_{k=1}^{2}l_{k}G_{kj}(T)}\bigg{)}\bigg{)},$
(4.5)

where $(l_{1},l_{2}):=(l,1-l)$ and where $G_{ij}$ and $\tau_{ij}$ are defined
by

$\displaystyle
G_{ij}(T):=\exp(-\alpha_{ij}(T)\tau_{ij}(T))\qquad(i,j\in\\{1,2\\}),$ (4.6)
$\displaystyle\tau_{ij}(T):=a_{ij}+\frac{b_{ij}}{T}\qquad(i\neq
j)\qquad\text{and}\qquad\tau_{ii}(T):=0\qquad(i\in\\{1,2\\})$ (4.7)
$\displaystyle\alpha_{ij}(T):=c_{ij}+d_{ij}T\qquad(i\neq
j)\qquad\text{and}\qquad\alpha_{ii}(T):=0\qquad(i\in\\{1,2\\}).$ (4.8)

As in the original model proposed in [34], we further assume the non-
randomness terms $\alpha_{ij}$ to be symmetric and temperature-independent,
that is,

$\displaystyle c_{ij}=c_{ji}\qquad\text{and}\qquad
d_{ij}=0\qquad(i,j\in\\{1,2\\}).$ (4.9)

Since the Antoine parameters have already been estimated by a prior parameter
estimation based on pure-component measurements, the parameters that remain to
be estimated for our two-component mixture are the five non-random two-liquid
parameters

$\displaystyle\theta:=(a_{12},a_{21},b_{12},b_{21},c_{12})\in\Theta\subset\mathbb{R}^{5}.$
(4.10)

In contrast to [34] and [12], we explicitly incorporate the non-randomness
parameter $c_{12}=c_{21}$ into the least-squares parameter estimation. In
[34], instead, this parameter is recommended to be manually selected as one of
the four values from the set $\\{0.2,0.3,0.4,0.47\\}$, based on the
characteristics of the mixture’s vapor-liquid equilibrium diagram.

Table 1: Antoine parameters of the two components, estimated from pure-component measurements in a prior parameter estimation step. Component | $A_{i}$ | $B_{i}$ | $C_{i}$
---|---|---|---
propanol (1) | 4.65413 | 1292.869 | -91.992
propyl acetate (2) | 3.84871 | 1088.392 | -90.571

Summarizing, the parametric model $f:\mathcal{X}\times\Theta\to\mathcal{Y}$
for the functional relationship between our input and output quantities (4.1)
is defined by the equations (4.2) and (4.3) in conjunction with (4.4)-(4.9)
and, moreover, the unknown parameters $\theta$ of $f$ that are to be estimated
are given by (4.10). Since (4.2) and (4.3) are algebraic equations of the form
(2.3), the model $f$ is an implicit parametric model in the sense of our
general framework (Section 2.1). Also, in view of the strong nonlinearity of
the non-random two-liquid model (4.1) w.r.t. $\theta$, it is clear that the
resulting implicit model $f$ is strongly nonlinear as well. In fact, $f$ is
intrinsically nonlinear to a substantial degree in the sense of [5, 36].

###  4.2 Apparatus and experimental setup

In this section, we describe the measurement apparatus and process along with
the initial experimental design.

####  4.2.1 Apparatus and procedure

All vapor-liquid equilibrium experiments for the considered two-component
system were carried out in an ILUDEST® measurement device, namely the FISCHER®
LABODEST® VLE 602 apparatus. Also, the two components propanol (CAS 71-23-8)
and propyl acetate (CAS 109-60-4) were taken from Geyer Chemsolute and from
Sigma-Aldrich, respectively, and both had a purity of at least
$0.995\,\mathrm{mol}/\mathrm{mol}$. At the beginning of each experiment, the
apparatus is cleaned, checked for leaks, and inerted. A mixture of the two
components is then poured into a heating zone and boiled at a constant heat
rate and pressure. A Cottrell pump maintains direct contact between liquid and
vapor to ensure equilibrium in a phase separation chamber. Vapor exiting the
phase separation chamber is condensed, mixed with liquid, and returned to the
heating zone. We assume vapor-liquid equilibrium to be reached when the
pressure and the temperature of the system have not changed for at least $30$
minutes. The system pressure is measured by a Wika P30 sensor with an average
absolute error of

$\displaystyle\sigma_{P}:=50\,\mathrm{Pa}$ (4.11)

and is controlled by valves connected to a nitrogen bottle and a vacuum pump.
And the system temperature, in turn, is measured by a Wika Pt100 temperature
sensor, located in the phase separation chamber, with an accuracy of

$\displaystyle\sigma_{T}:=0.03\,\mathrm{K}.$ (4.12)

At equilibrium, the liquid and vapor phases are sampled after the phase
separation chamber, and then their compositions are determined with an average
absolute error of

$\displaystyle\sigma_{l},\sigma_{v}:=0.0015\,\mathrm{mol/mol},$ (4.13)

according to the procedure described next. So, for each experiment, there are
measurement data for liquid composition $(l,1-l)$ and pressure $P$ on the one
hand and for the vapor composition $(v,1-v)$ and the equilibrium temperature
$T$ on the other hand.

####  4.2.2 Sample analytics

In order to determine the compositions of the liquid and vapor phases at
equilibrium, the samples taken from the phase separation chamber are first
diluted with ethanol (CAS 64-17-5) at a $1:4$ sample-to-ethanol mass ratio and
then measured on an Agilent gas chromatograph 6850 equipped with a thermal
conductivity detector and a gas chromatography column, namely an Agilent HP-
PLOT Q column of length $30\,\mathrm{m}$, inner diameter $0.32\,\mathrm{mm}$,
and film thickness $0.02\,\mathrm{mm}$. The diluting ethanol was obtained from
VWR Chemicals and its purity was at least $0.999\,\mathrm{mol}/\mathrm{mol}$.
The column inlet was heated to $230\,^{\circ}\mathrm{C}$, helium was used as a
carrier gas, and the split ratio was set to $40:1$. The oven stayed heated at
$190\,^{\circ}\mathrm{C}$ for $10$ minutes. And finally, two calibration
curves for lower and higher concentrations of both components were used to
provide concentration measurements with an average absolute error as specified
in (4.13) above.

####  4.2.3 Initial design and measurement results

We performed the first batch of experiments according to the fractional
factorial design init which consists of the four corner points and the center
point of the box

$\displaystyle[\underline{l},\overline{l}]\times[\underline{P},\overline{P}]:=[0.05,0.95]\times[1\cdot
10^{5},3\cdot 10^{5}]\subset\mathcal{X}$ (4.14)

and of one additional point next to the center point (Figure 2). Specifically,
the initial design init consists of the six $(l,P)$ points listed in Table 2.
With its $6>5=d_{\mathrm{\theta}}$ experimental points, the initial design of
course satisfies the necessary condition (2.16) for the invertibility of the
information matrices of $f$ at init. It is important to notice that the
experimental points of init – and, for that matter, of any other experimental
design – cannot be realized exactly in the lab. And so there is always a
mismatch between planned and actually performed experiments (which is often
referred to as errors in variables in the context of parameter estimation
[36]). In particular, this is true for the planned and the actually realized
liquid compositions. Apart from the planned experimental points $(l,P)$, we
therefore also list the actual measurement locations $(l^{\prime},P^{\prime})$
in Table 2. And, along with them, we list the measured values $(v,T)$ of our
output quantities.

Figure 2: Illustration of the experimental points from the initial design init.  Table 2: Initial experiments conducted in the study. | $l$ | $l^{\prime}$ | $P$ | $P^{\prime}$ | $v$ | $T$ | $\sigma_{v}$ | $\sigma_{T}$
---|---|---|---|---|---|---|---|---
init | 0.050000 | 0.0456 | 100000.0 | 99990.0 | 0.0813 | 372.21 | 0.0015 | 0.03
init | 0.050000 | 0.6961 | 300000.0 | 299970.0 | 0.7243 | 401.50 | 0.0015 | 0.03
init | 0.500000 | 0.4466 | 200000.0 | 199950.0 | 0.5257 | 389.12 | 0.0015 | 0.03
init | 0.950000 | 0.9728 | 100000.0 | 99990.0 | 0.9611 | 369.16 | 0.0015 | 0.03
init | 0.950000 | 0.9642 | 300000.0 | 299981.6 | 0.9546 | 401.46 | 0.0015 | 0.03
init | 0.612500 | 0.6426 | 200000.0 | 199950.0 | 0.6730 | 388.14 | 0.0015 | 0.03

###  4.3 Application-specific settings for parameter estimation and optimal
design

As a general modeling and optimization framework supporting automatic
differentiation, we choose the software package Casadi [1]. Additionally, in
the present application example, we use the following specific settings for
the occurring parameter estimation and optimal design problems. In keeping
with the accuracies (4.12) and (4.13) of the employed measurement devices, we
take the covariance matrix of the measurement errors to be

$\displaystyle\varsigma:=\begin{pmatrix}\sigma_{1}^{2}&0\\\
0&\sigma_{2}^{2}\end{pmatrix}=\begin{pmatrix}\sigma_{v}^{2}&0\\\
0&\sigma_{T}^{2}\end{pmatrix}=\begin{pmatrix}2.25\cdot 10^{-6}&0\\\ 0&9\cdot
10^{-4}\end{pmatrix}\in\mathbb{R}_{\mathrm{pd}}^{d_{\mathrm{y}}\times
d_{\mathrm{y}}}.$ (4.15)

In particular, by the diagonality of this matrix, we implicitly assume the
measurement errors for vapor composition and the measurement errors for
temperature to be independent of each other. It should be recalled from (2.10)
and (2.17) that the covariance of the measurement errors enters all least-
squares estimation problems and – via the information matrices – also all
locally optimal experimental design problems. And so the above covariance
matrix (4.15) features in all optimization problems of our application
example.

####  4.3.1 Settings for parameter estimation

Since the covariance matrix (4.15) of the measurement errors is diagonal, all
least-squares estimation problems in our case study take the special form

$\displaystyle\min_{\theta\in\Theta}\sum_{i=1}^{n}\sum_{j=1}^{2}(f_{j}(x_{i},\theta)-y_{ij})^{2}/\sigma_{j}^{2}$
(4.16)

of weighted least-squares problems. In order to solve these problems, we use
the interior-point optimizer Ipopt [45] (version 3.10.4) – installed by Casadi
– together with the linear solver MA57 [13, 21]. Specifically, the optimizer
settings were chosen according to the default settings of Casadi (version
3.5.5). All model derivatives required for the solution of (4.16) are provided
by Casadi’s automatic differentiation functionality. In order to avoid ending
up with suboptimal local minima for (4.16), we use a simple multistart
procedure which randomly selects candidate start values $\theta^{0}$ and then
filters them according to their objective function value. Specifically, if the
objective function value at a candidate start value $\theta^{0}$ is
sufficiently small, it is accepted as a start value for Ipopt; otherwise it is
rejected.

####  4.3.2 Settings for optimal design

In all locally optimal experimental design problems of our case study, we use
the D-criterion as the underlying design criterion $\Psi$, that is,

$\displaystyle\Psi:=\Psi_{0}$ (4.17)

as defined in (2.25.a). As the importance factor for already performed
experiments, we consistently use $\alpha:=0.5$ throughout the case study. In
this manner, we ensure that previous and prospective experiments are given
equal importance throughout the entire sequential design procedure. As the
maximal tolerable optimality error for the locally $\varepsilon$-optimal
weighted designs computed by Algorithm 1, we choose $\varepsilon:=5\cdot
10^{-5}$. In order to solve the discretized versions (3.4) of the weight
optimization problem (3.3), we use the solver SLSQP from scipy.optimize [44]
with its default settings (version 12.1.0). (It should be noted that other
solvers could have been used just as well. After all, as has been explained
after Algorithm 1, the discretized problems (3.4) are just convex weight
optimization problems with only as many variables as the size of the
discretization $\underline{X}^{k}$.) All model derivatives required for the
computation of the information matrices in (3.4) and (3.5) are provided by
Casadi’s automatic differentiation functionality. As the minimal total weight
of the candidate experiments that remain after our sorting-out procedure
(Algorithm 2), we choose $\underline{w}^{+}:=0.95$. And as the maximal number
of new experiments per iteration of our sequential design algorithm, we choose
$\overline{n}^{+}:=3$. We make this choice because, from a comparison
perspective, it fits well with our reference factorial designs of the sizes
$3\cdot 3$, $5\cdot 3$, and $9\cdot 3$, respectively. In order to initialize
our sequential design procedure (Algorithm 3), we choose the $6$-point
fractional factorial design

$\displaystyle\tilde{x}^{0}:=\texttt{init}$ (4.18)

introduced above (Figure 2 and Table 2). And finally, the maximal number of
iterations and the progress tolerance of our sequential design algorithm are
chosen to be $\overline{n}:=27$ and $\delta:=0.1$ in our case study.
Summarizing, we apply Algorithm 3 with the following inputs:

$\displaystyle\alpha:=0.5,\qquad\varepsilon:=0.5\cdot
10^{-4},\qquad\underline{w}^{+}:=0.95,$ (4.19)
$\displaystyle\overline{n}^{+}:=3,\qquad\overline{n}:=27,\qquad\delta:=0.1.$
(4.20)

###  4.4 Considered designs, measurement results, and paremeter estimates

After performing the six experiments from the initial fractional factorial
design init (Table 2), we continue to measure the output quantities (4.1) at a
number of additional input points, namely first at a sequence

$\displaystyle\texttt{fed}^{1}=\texttt{fed}^{0}\,\&\,\texttt{fed}^{0,+},\qquad\texttt{fed}^{2}=\texttt{fed}^{1}\,\&\,\texttt{fed}^{1,+},\qquad\texttt{fed}^{3}=\texttt{fed}^{2}\,\&\,\texttt{fed}^{2,+}$
(4.21)

of factorial designs $\texttt{fed}^{1},\texttt{fed}^{2},\texttt{fed}^{3}$, and
then at a sequence

$\displaystyle\texttt{oed}^{1}=\texttt{oed}^{0}\,\&\,\texttt{oed}^{0,+},\qquad\texttt{oed}^{2}=\texttt{oed}^{1}\,\&\,\texttt{oed}^{1,+},\qquad\texttt{oed}^{3}=\texttt{oed}^{2}\,\&\,\texttt{oed}^{2,+}$
(4.22)

of locally optimal designs
$\texttt{oed}^{1},\texttt{oed}^{2},\texttt{oed}^{3}$. We used the notation
(2.6) for refined designs here. In particular, all designs from (4.21) arise
by a successive refinement of the $5$-point fractional factorial design
$\texttt{fed}^{0}$ consisting of the first five experiments from Table 2 (all
points in Figure 2 except for the off-center point). Similarly, all designs
from (4.22) arise by a successive refinement of the $6$-point fractional
factorial design $\texttt{oed}^{0}=\texttt{init}$ consisting of the six
experiments from Table 2.

As the designs $\texttt{fed}^{1},\texttt{fed}^{2},\texttt{fed}^{3}$ in (4.21),
we choose full factorial designs on the $9\cdot 3$ equidistant grid

$\displaystyle\underline{\mathcal{X}}^{\mathrm{fed}}:=\big{\\{}(\underline{l}+i(\overline{l}-\underline{l})/8,\underline{P}+j(\overline{P}-\underline{P})/2):i\in\\{0,1,\dots,8\\}\text{
and }j\in\\{0,1,2\\}\big{\\}}\subset\mathcal{X},$ (4.23)

where $\underline{l}:=0.05$ and $\overline{l}:=0.95$ are the mole fraction
bounds from (4.14) and $\underline{P},\overline{P}$ are the pressure bounds
from (4.1). Specifically, $\texttt{fed}^{1},\texttt{fed}^{2},\texttt{fed}^{3}$
are the $3\cdot 3$, $5\cdot 3$, and $9\cdot 3$ full factorial designs depicted
in Figure 3 and defined by the $l$\- and $P$-columns in Table 3.

Figure 3: Illustration of the experimental points from the factorial designs $\texttt{fed}^{1}$ (blue triangles), $\texttt{fed}^{2}=\texttt{fed}^{1}\,\&\,\texttt{fed}^{1,+}$, and $\texttt{fed}^{3}=\texttt{fed}^{2}\,\&\,\texttt{fed}^{2,+}$. The experimental points of the designs $\texttt{fed}^{1,+}$ and $\texttt{fed}^{2,+}$ are the red pentagons and the green stars, respectively.  Table 3: Sequence of factorial experiments (4.21) conducted in the study (following the five initial experiments of $\texttt{fed}^{0}$, that is, the first five experiments from init). | $l$ | $l^{\prime}$ | $P$ | $P^{\prime}$ | $v$ | $T$ | $\sigma_{v}$ | $\sigma_{T}$
---|---|---|---|---|---|---|---|---
$\texttt{fed}^{0,+}$ | 0.050000 | 0.0444 | 200000.0 | 199950.0 | 0.0815 | 396.06 | 0.0015 | 0.03
$\texttt{fed}^{0,+}$ | 0.500000 | 0.4469 | 100000.0 | 99990.0 | 0.4979 | 367.32 | 0.0015 | 0.03
$\texttt{fed}^{0,+}$ | 0.500000 | 0.3538 | 300000.0 | 299981.6 | 0.4560 | 404.32 | 0.0015 | 0.03
$\texttt{fed}^{0,+}$ | 0.950000 | 0.9685 | 200000.0 | 199950.0 | 0.9525 | 388.73 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.275000 | 0.2440 | 100000.0 | 99990.0 | 0.3341 | 368.70 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.275000 | 0.2429 | 200000.0 | 199950.0 | 0.3431 | 391.58 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.275000 | 0.1168 | 300000.0 | 299981.6 | 0.1933 | 409.82 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.725000 | 0.6974 | 100000.0 | 99990.0 | 0.6750 | 367.12 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.725000 | 0.6958 | 200000.0 | 199950.0 | 0.7099 | 387.82 | 0.0015 | 0.03
$\texttt{fed}^{1,+}$ | 0.725000 | 0.5809 | 300000.0 | 299981.6 | 0.6499 | 402.05 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.162500 | 0.1503 | 100000.0 | 99990.0 | 0.2385 | 369.95 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.162500 | 0.1486 | 200000.0 | 199950.0 | 0.2432 | 393.22 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.162500 | 0.0454 | 300000.0 | 299981.6 | 0.0815 | 411.82 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.387500 | 0.3707 | 100000.0 | 99990.0 | 0.4378 | 367.46 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.387500 | 0.3697 | 200000.0 | 199950.0 | 0.4430 | 389.46 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.387500 | 0.2346 | 300000.0 | 299981.6 | 0.3424 | 406.86 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.612500 | 0.6037 | 100000.0 | 99990.0 | 0.6026 | 367.05 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.612500 | 0.6426 | 200000.0 | 199950.0 | 0.6730 | 388.14 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.612500 | 0.4356 | 300000.0 | 299981.6 | 0.5353 | 403.59 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.837500 | 0.9168 | 100000.0 | 99990.0 | 0.8869 | 368.41 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.837500 | 0.9037 | 200000.0 | 199950.0 | 0.8887 | 388.39 | 0.0015 | 0.03
$\texttt{fed}^{2,+}$ | 0.837500 | 0.9067 | 300000.0 | 299981.6 | 0.8903 | 401.36 | 0.0015 | 0.03

As the designs $\texttt{oed}^{1},\texttt{oed}^{2},\texttt{oed}^{3}$ in (4.22),
we choose the locally optimal designs computed in the course of our sequential
design algorithm (Algorithm 3) with design criterion (4.17), inputs
(4.18)-(4.20), and with the $10\cdot 10$ equidistant grid

$\displaystyle\underline{\mathcal{X}}^{\mathrm{oed}}:=\big{\\{}(i/9,\underline{P}+j(\overline{P}-\underline{P})/9):i,j\in\\{0,1,\dots,9\\}\big{\\}}\subset\mathcal{X},$
(4.24)

as the underlying design space, where again $\underline{P},\overline{P}$ are
the pressure bounds from (4.1). With these settings, the sequential design
algorithm terminates after just three iterations because all experimental
points of the design $\texttt{oed}^{2,+}$ proposed in the third iteration
already belong to the previous design $\texttt{oed}^{2}$. And therefore we end
up with a sequence of just three locally optimal designs
$\texttt{oed}^{1},\texttt{oed}^{2},\texttt{oed}^{3}$. See Figure 4 and the
$l$\- and $P$-columns in Table 4 for the experimental points making up these
designs. In particular, the designs
$\texttt{oed}^{1},\texttt{oed}^{2},\texttt{oed}^{3}$ have the sizes $9$, $12$
and $15$. Also, while $\texttt{oed}^{1}$ is a design without replications,
$\texttt{oed}^{2}$ and $\texttt{oed}^{3}$ contain replicated – or more
precisely, duplicated – experiments. Specifically, $\texttt{oed}^{1}$ contains
one duplicated experiment and $\texttt{oed}^{2}$ contains four duplicated
experiments. See the overlaps in the figure.

Figure 4: Illustration of the experimental points from the optimal designs $\texttt{oed}^{0}=\texttt{init}$ (gray circles), and $\texttt{oed}^{1}=\texttt{oed}^{0}\,\&\,\texttt{oed}^{0,+}$, $\texttt{oed}^{2}=\texttt{oed}^{1}\,\&\,\texttt{oed}^{1,+}$, and $\texttt{oed}^{3}=\texttt{oed}^{2}\,\&\,\texttt{oed}^{2,+}$. The experimental points of the designs $\texttt{oed}^{0,+}$, $\texttt{oed}^{1,+}$ and $\texttt{oed}^{2,+}$ are the blue diamonds, the green pluses and the red triangles, respectively.  Table 4: Sequence of optimal experiments (4.22) conducted in the study (following the six initial experiments from $\texttt{oed}^{0}=\texttt{init}$). | $l$ | $l^{\prime}$ | $P$ | $P^{\prime}$ | $v$ | $T$ | $\sigma_{v}$ | $\sigma_{T}$
---|---|---|---|---|---|---|---|---
$\texttt{oed}^{0,+}$ | 0.222222 | 0.1731 | 300000.0 | 299981.6 | 0.2748 | 407.73 | 0.0015 | 0.03
$\texttt{oed}^{0,+}$ | 0.444444 | 0.3878 | 100000.0 | 99990.0 | 0.4561 | 367.49 | 0.0015 | 0.03
$\texttt{oed}^{0,+}$ | 0.666667 | 0.6520 | 300000.0 | 299981.6 | 0.6952 | 401.62 | 0.0015 | 0.03
$\texttt{oed}^{1,+}$ | 0.111111 | 0.1018 | 100000.0 | 99990.0 | 0.1721 | 370.81 | 0.0015 | 0.03
$\texttt{oed}^{1,+}$ | 0.444444 | 0.4518 | 100000.0 | 99990.0 | 0.4992 | 367.23 | 0.0015 | 0.03
$\texttt{oed}^{1,+}$ | 0.888889 | 0.8845 | 100000.0 | 99990.0 | 0.8459 | 368.03 | 0.0015 | 0.03
$\texttt{oed}^{2,+}$ | 0.111111 | 0.1094 | 100000.0 | 99990.0 | 0.1803 | 370.68 | 0.0015 | 0.03
$\texttt{oed}^{2,+}$ | 0.222222 | 0.2385 | 300000.0 | 300000.0 | 0.3466 | 406.26 | 0.0015 | 0.03
$\texttt{oed}^{2,+}$ | 0.666667 | 0.7372 | 300000.0 | 300000.0 | 0.7586 | 401.28 | 0.0015 | 0.03

Just like for the initial design, Tables 3 and 4 not only specify the planned
experimental designs $\texttt{fed}^{k}$ and $\texttt{oed}^{k}$, but also
detail the input conditions $(l^{\prime},P^{\prime})$ that were actually
realized in the lab experiments, along with the corresponding measurement
results $(v,T)$ for the outputs. In particular, we see that the planned
replications in $\texttt{oed}^{2}$ and in $\texttt{oed}^{3}$ represented by
overlaps in Figure 4 are fanned out a bit, in reality.

In the following, we write
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$ for the totality of the
collected experimental data, that is, for the combination of all the $(l,P)$\-
and $(v,T)$-values from Tables 2, 3 and 4. We also write
$\texttt{tot}:=\tilde{x}^{\mathrm{tot}}$ for the total combined design. In
formulas,

$\displaystyle\tilde{x}^{\mathrm{tot}}:=\texttt{tot}:=\texttt{fed}^{0}\,\&\,\texttt{fed}^{0,+}\,\&\,\texttt{fed}^{1,+}\,\&\,\texttt{fed}^{2,+}\,\&\,\texttt{oed}^{0,+}\,\&\,\texttt{oed}^{1,+}\,\&\,\texttt{oed}^{2,+},$
(4.25)

so that the total combined design tot has the size $n^{\mathrm{tot}}=36$.
Table 5 records our final least-squares estimate for the non-random two liquid
parameters of propanol and propyl acetate, that is, the least-squares estimate
$\widehat{\theta}^{\texttt{tot}}:=\widehat{\theta}_{f}(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$
based on all collected data
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$. Figure 5 shows the
corresponding predicted bubble- and dew-point temperature curves for the three
pressures considered in the experiments, along with the measured temperature
values. As can be seen, the predictions are in perfect agreement with the
measurements. Additionally, the figure nicely illustrates the narrow
azeotropic behavior of our mixture.

Table 5: The final least-squares estimate $\widehat{\theta}^{\texttt{tot}}:=\widehat{\theta}_{f}(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$ for the non-random two liquid parameters of propanol (1) and propyl acetate (2), based on all collected data $(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$. $a_{1,2}$ | $a_{2,1}$ | $b_{1,2}$ | $b_{2,1}$ | $c_{1,2}=c_{2,1}$
---|---|---|---|---
9.396525 | -10.305843 | -786.446701 | 1510.352034 | 0.010000
Figure 5: The predicted and measured bubble- and dew-point temperatures for
propanol and propyl acetate at the pressures $1\,\mathrm{bar}$ (bottom-most
curves), $2\,\mathrm{bar}$ (middle curves), and $3\,\mathrm{bar}$ (top-most
curves). The predictions are based on the final least-squares estimate
$\widehat{\theta}^{\texttt{tot}}$ from Table 5.

###  4.5 Comparison of the considered designs and of the resulting models

In this section, we finally compare the prediction quality of the models
trained on the factorial designs $\texttt{fed}^{k}$ with the prediction
quality of the models trained on the locally optimal designs
$\texttt{oed}^{k}$ computed by our sequential design algorithm. In order to do
so, we use the three prediction quality measures (3.14), (3.21.a) and (3.21.b)
from Section 3.3. As has been pointed out there, these prediction quality
measures all depend on the totality
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$ of the collected
experimental data as reference data, see (4.25) above.

We begin by considering the root mean squared prediction errors
$\rho_{j}(\tilde{x})$ for the models trained on the different designs, see
(3.14). Since our model has two outputs (4.1), there are two such prediction
errors for each design $\tilde{x}$ and we denote them by
$\rho_{v}(\tilde{x}):=\rho_{1}(\tilde{x})$ and
$\rho_{T}(\tilde{x}):=\rho_{2}(\tilde{x})$, respectively. Table 6 lists these
values for all considered designs $\tilde{x}$, and Figure 6 plots them as
functions of the design size.

Table 6: Root mean squared prediction errors for the models trained on the different designs | init | $\texttt{oed}^{1}$ | $\texttt{fed}^{1}$ | $\texttt{oed}^{2}$ | $\texttt{fed}^{2}$ | $\texttt{oed}^{3}$ | $\texttt{fed}^{3}$ | tot
---|---|---|---|---|---|---|---|---
$\operatorname{size}(\tilde{x})$ | 6 | 9 | 9 | 12 | 15 | 15 | 27 | 36
$\rho_{v}(\tilde{x})/(10^{-4})$ | 72.07 | 59.96 | 65.74 | 59.59 | 63.21 | 59.61 | 59.86 | 58.95
$\rho_{T}(\tilde{x})/(10^{-2})$ | 24.80 | 15.83 | 19.50 | 15.10 | 18.26 | 16.18 | 15.75 | 14.63
Figure 6: Root mean squared prediction errors for the models trained on the
different designs. The open squares and open circles represent respectively
the designs $\texttt{fed}^{k}$ and $\texttt{oed}^{k}$, while the closed circle
and the closed square represent the design init and tot, respectively. As can
be seen, $\texttt{oed}^{1}$ with its $9$ experiments already achieves the same
prediction error as $\texttt{fed}^{3}$ with its $27$ experiments.

As can be seen from this figure, the prediction error
$\rho_{j}(\texttt{oed}^{3})$ based on the size-$15$ locally optimal design
$\texttt{oed}^{3}$ is practically the same as the prediction error
$\rho_{j}(\texttt{fed}^{3})$ based on the size-$27$ full factorial design
$\texttt{fed}^{3}$. In short,

$\displaystyle\rho_{j}(\texttt{oed}^{3})\approx\rho_{j}(\texttt{fed}^{3})\qquad(j\in\\{1,2\\})$
(4.26)

even though $\texttt{oed}^{3}$ is only about half the size of
$\texttt{fed}^{3}$. Spelled out, (4.26) means that on the totality of
collected data, the model trained just on the $15$ optimally designed
experiments from $\texttt{oed}^{3}$ makes essentially the same prediction
error as the model trained on the $27$ factorially designed experiments from
$\texttt{fed}^{3}$.

Additionally, the figure and the table show that the prediction errors
$\rho_{j}(\texttt{oed}^{3})$ and $\rho_{j}(\texttt{fed}^{3})$ from (4.26) are
not only essentially equal to each other, but also almost equal to
$\rho_{j}(\texttt{tot})$ which, in turn, is the smallest possible prediction
error value on our total data set. After all,

$\displaystyle\sum_{j=1}^{2}\rho_{j}(\texttt{tot})^{2}/\sigma_{j}^{2}=\min_{\theta\in\Theta}\frac{1}{n^{\mathrm{tot}}}\sum_{i=1}^{n^{\mathrm{tot}}}\sum_{j=1}^{2}(f_{j}(x^{\mathrm{tot}}_{i},\theta)-y^{\mathrm{tot}}_{ij})^{2}/\sigma_{j}^{2}$
(4.27)

by the very definition of the least-squares estimate
$\widehat{\theta}_{f}(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$ for
$(\tilde{x}^{\mathrm{tot}},\tilde{y}^{\mathrm{tot}})$.

And finally, the figure and the table reveal that the prediction error is
essentially constant along the designs $\texttt{oed}^{k}$ computed by our
algorithm, while it is strictly decreasing along our factorial designs
$\texttt{fed}^{k}$. In particular, the terminal prediction error
$\rho_{j}(\texttt{oed}^{3})\approx\rho_{j}(\texttt{fed}^{3})$ is reached
already with $\texttt{oed}^{1}$ and $\texttt{oed}^{2}$, while it is missed
with $\texttt{fed}^{1}$ or $\texttt{fed}^{2}$. So, from a pure prediction-
error perspective, it would have made sense to terminate our algorithm even
earlier. Yet, from a prediction-uncertainty perspective, a premature
termination turns out to be not advisable, as we will see in a moment.

We continue by considering the worst-case prediction uncertainties
$\sigma_{j}^{\mathrm{lin}}(\tilde{x})$ and
$\sigma_{j}^{\mathrm{sam}}(\tilde{x})$ for the models trained on the different
designs, see (3.21). In analogy to the root mean squared prediction errors, we
denote them by
$\sigma^{\mathrm{lin}}_{v}(\tilde{x}):=\sigma_{1}^{\mathrm{lin}}(\tilde{x})$,
$\sigma^{\mathrm{lin}}_{T}(\tilde{x}):=\sigma_{2}^{\mathrm{lin}}(\tilde{x})$
and by
$\sigma^{\mathrm{sam}}_{v}(\tilde{x}):=\sigma_{1}^{\mathrm{sam}}(\tilde{x})$,
$\sigma^{\mathrm{sam}}_{T}(\tilde{x}):=\sigma_{2}^{\mathrm{sam}}(\tilde{x})$,
respectively, and for the sampling-based uncertainties we used samples of the
size

$\displaystyle n^{\mathrm{sam}}:=1000.$ (4.28)

Tables 7 and 8 record these values for all designs $\tilde{x}$ considered
here.

Table 7: Worst-case prediction uncertainties (linearization-based) for the models trained on the different designs | init | $\texttt{oed}^{1}$ | $\texttt{fed}^{1}$ | $\texttt{oed}^{2}$ | $\texttt{fed}^{2}$ | $\texttt{oed}^{3}$ | $\texttt{fed}^{3}$ | tot
---|---|---|---|---|---|---|---|---
$\operatorname{size}(\tilde{x})$ | 6 | 9 | 9 | 12 | 15 | 15 | 27 | 36
$\sigma^{\mathrm{lin}}_{v}(\tilde{x})/(10^{-4})$ | 67.90 | 33.08 | 32.08 | 28.11 | 27.46 | 25.47 | 24.92 | 23.07
$\sigma^{\mathrm{lin}}_{T}(\tilde{x})/(10^{-2})$ | 34.14 | 10.79 | 10.87 | 10.40 | 8.96 | 8.26 | 8.37 | 7.85
Table 8: Worst-case prediction uncertainties (sampling-based) for the models trained on the different designs | init | $\texttt{oed}^{1}$ | $\texttt{fed}^{1}$ | $\texttt{oed}^{2}$ | $\texttt{fed}^{2}$ | $\texttt{oed}^{3}$ | $\texttt{fed}^{3}$ | tot
---|---|---|---|---|---|---|---|---
$\operatorname{size}(\tilde{x})$ | 6 | 9 | 9 | 12 | 15 | 15 | 27 | 36
$\sigma^{\mathrm{sam}}_{v}(\tilde{x})/(10^{-4})$ | 17.3 | 10.78 | 9.84 | 6.48 | 6.98 | 5.60 | 4.49 | 3.71
$\sigma^{\mathrm{sam}}_{T}(\tilde{x})/(10^{-2})$ | 7.56 | 3.08 | 3.75 | 2.53 | 2.39 | 1.94 | 1.36 | 1.13

As in the prediction-error case considered above, it can be seen from the
prediction-uncertainty tables that the worst-case prediction uncertainties
$\sigma_{j}^{\mathrm{lin}}(\texttt{oed}^{3})$ and
$\sigma_{j}^{\mathrm{sam}}(\texttt{oed}^{3})$ based on the size-$15$ locally
optimal design $\texttt{oed}^{3}$ are practically the same as the worst-case
prediction uncertainties $\sigma_{j}^{\mathrm{lin}}(\texttt{fed}^{3})$ and
$\sigma_{j}^{\mathrm{sam}}(\texttt{fed}^{3})$ based on the size-$27$ full
factorial design $\texttt{fed}^{3}$. In short,

$\displaystyle\sigma_{j}^{\mathrm{lin}}(\texttt{oed}^{3})\approx\sigma_{j}^{\mathrm{lin}}(\texttt{fed}^{3})\qquad\text{and}\qquad\sigma_{j}^{\mathrm{sam}}(\texttt{oed}^{3})\approx\sigma_{j}^{\mathrm{sam}}(\texttt{fed}^{3})$
(4.29)

despite the substantially different sizes of $\texttt{oed}^{3}$ and
$\texttt{fed}^{3}$. In words, (4.29) means that the predictions of the model
trained just on the $15$ optimally designed experiments from
$\texttt{oed}^{3}$ can be trusted basically as much as the predictions of the
model trained on the $27$ factorially designed experiments from
$\texttt{fed}^{3}$.

As opposed to the prediction-error case, the terminal prediction uncertainties
from (4.29) are no longer reached already for $\texttt{oed}^{1}$ or
$\texttt{oed}^{2}$. Instead, the prediction uncertainties are strictly
improving not only along our factorial designs $\texttt{fed}^{k}$, but also
along the optimal designs $\texttt{oed}^{k}$ computed by our algorithm. And
therefore, from a prediction-uncertainty perspective, a premature termination
of our algorithm would be disadvantageous.

Figure 7: Worst-pressure-case prediction uncertainties (linearization-based)
for the models trained on the different designs. The prediction uncertaintiy
curves (4.30) for init, $\texttt{oed}^{1}$, $\texttt{fed}^{1}$,
$\texttt{oed}^{2}$, $\texttt{fed}^{2}$, $\texttt{oed}^{3}$, $\texttt{fed}^{3}$
and tot are represented, respectively, by the gray circles, the blue diamonds,
the green pluses, the blue circles, the red triangles, the green stars, the
gray triangles, and the red squares. The orange lines indicate the measurment
uncertainties $\sigma_{v}$ and $\sigma_{T}$. Comparing the maxima of the
green-stars and the gray-triangles curves, we see that the worst-case
prediction uncertainty achieved with the size-$15$ design $\texttt{oed}^{3}$
is practically the same as the worst-case prediction uncertainty achieved with
the size-$27$ design $\texttt{fed}^{3}$. Figure 8: Worst-pressure-case
prediction uncertainties (sampling-based) for the models trained on the
different designs. The prediction uncertainty curves (4.30) for init,
$\texttt{oed}^{1}$, $\texttt{fed}^{1}$, $\texttt{oed}^{2}$,
$\texttt{fed}^{2}$, $\texttt{oed}^{3}$, $\texttt{fed}^{3}$ and tot are
represented, respectively, by the gray circles, the blue diamonds, the green
pluses, the blue circles, the red triangles, the green stars, the gray
triangles, and the red squares. The orange lines indicate the measurment
uncertainties $\sigma_{v}$ and $\sigma_{T}$. Comparing the maxima of the
green-stars and the gray-triangles curves, we see that the worst-case
prediction uncertainty achieved with the size-$15$ design $\texttt{oed}^{3}$
is practically the same as the worst-case prediction uncertainty achieved with
the size-$27$ design $\texttt{fed}^{3}$.

In order to highlight the variability of the prediction uncertainties (3.20),
we plot them as functions of one of our inputs $x=(l,P)$, namely the liquid
mole fraction $l$ of propanol (Figures 7 and 8). Specifically, the curves in
Figures 7 and 8 depict the worst-pressure-case prediction uncertainties

$\displaystyle
l\mapsto\max_{P\in[\underline{P},\overline{P}]}\sigma_{j}^{\mathrm{lin}}((l,P),\tilde{x})\qquad\text{and}\qquad
l\mapsto\max_{P\in[\underline{P},\overline{P}]}\sigma_{j}^{\mathrm{sam}}((l,P),\tilde{x})$
(4.30)

as functions of $l$ for all the considered designs $\tilde{x}$. Clearly, the
maxima of these worst-pressure-case curves are nothing but the worst-case
prediction uncertainties $\sigma_{j}^{\mathrm{lin}}(\tilde{x})$ and
$\sigma_{j}^{\mathrm{sam}}(\tilde{x})$ that we considered above (Tables 7 and
8).

We immediately observe from the figures that the worst-pressure-case
prediction uncertainty curves (4.30) are strongly varying with $l$, in
general. And so, the prediction uncertainties (3.20) are strongly dependent on
the inputs $x$ as well, in general. We also observe from the figures that the
variability of the depicted curves (4.30) becomes less and less pronounced,
the larger the size of the design $\tilde{x}$. In other words, the curves
(4.30) get closer and closer to being constants equal to the worst-case
prediction uncertainties $\sigma_{j}^{\mathrm{lin}}(\tilde{x})$ and
$\sigma_{j}^{\mathrm{sam}}(\tilde{x})$. It should be noticed, however, that
the curves (4.30) will alway retain some variability – they can never become
exactly constant just because at the boundaries $l=0$ (pure propyl acetate)
and $l=1$ (pure propanol) they take exactly the value $0$. In short,

$\displaystyle\sigma_{j}^{\mathrm{lin}}((0,P),\tilde{x})=0=\sigma_{j}^{\mathrm{lin}}((1,P),\tilde{x})\qquad\text{and}\qquad\sigma_{j}^{\mathrm{sam}}((0,P),\tilde{x})=0=\sigma_{j}^{\mathrm{sam}}((1,P),\tilde{x})$
(4.31)

for every arbitrary pressure value $P$ and for every arbitrary experimental
design $\tilde{x}$. Indeed, the simple reason behind (4.31) is that for the
pure-component cases $l=0$ and $l=1$, the non-random two-liquid model (4.1)
does not depend on the parameters $\theta$ at all. And therefore, the same is
true for the solutions $T(l,P,\theta)$ and $v(l,P,\theta)$ of (4.2) and (4.3).
So, for the boundary values $l=0$ and $l=1$, the model predictions

$\displaystyle f(l,P,\theta):=(v(l,P,\theta),T(l,P,\theta))^{\top}$ (4.32)

are actually completely independent of $\theta$ and this, in conjunction with
the definitions (3.20), immediately implies the relations claimed in (4.31)
above and evidenced by the figures.

##  5 Conclusions and outlook

In this paper, we have proposed a general method of sequential locally optimal
experimental design for nonlinear models. It can be applied to estimate the
parameters of nonlinear models as they abound in vapor-liquid equilibrium
modeling, for instance. And these nonlinear models can be defined either by
explicit or by implicit algebraic equations. In every iteration of our method,
a whole batch of new experiments is computed in a two-stage locally optimal
manner. In other words, a two-stage locally optimal experimental design
problem is solved in every iteration of our method, based on the previously
collected experimental data and on the corresponding least-squares estimate.
As the underlying design criterion, any of the standard differentiable design
criteria can be chosen, in particular, the A- and the D-criterion. Also, the
number of new experiments computed in each iteration can be freely specified
by the user in such a way that it fits the concrete experimental circumstances
best.

We have demonstrated the benefits of the general method by applying it to the
estimation of non-random two-liquid parameters of a narrow azeotropic system.
In particular, we have found that the proposed sequential locally optimal
experimental design method requires substantially fewer experiments than
traditional factorial experimental design to achieve the same model precision
and prediction quality. Specifically, compared to factorial design, our
sequential design method required only about half the experiments and thus
helps to considerably reduce the experimental effort for the non-random two-
liquid parameter estimation.

In future work, we plan to demonstrate the benefits of the proposed sequential
design methodolgy also on other application examples, especially, on systems
with liquid-liquid equilibria and on systems with chemical reactions. We also
plan to extend the proposed sequential design methodolgy to experimental
design problems with constraints on the designs, like cost or location
constraints. In practice, such constraints are highly relevant, but are rarely
taken into account systematically. And finally, we will investigate how
experimental design can further be improved with the help of other design
criteria than the linearization-based design criteria used in this work. In
particular, it will be interesting to see how and to what extent experimental
design for nonlinear models improves when the designs are computed using
suitable sampling-based approximations to the parameter uncertainties or to
the prediction uncertainties, as opposed to the successive linearization-based
approximations used in this work.

## Acknowledgments

We gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – project number 466397921 – within the
Priority Programme “SPP 2331: Machine Learning in Chemical Engineering”.

## Author contributions

M. Bubel: conceptualization, software, validation, data curation,
visualization, writing - review and editing. J. Schmid: conceptualization,
methodology, formal analysis, validation, writing - original draft, writing -
review and editing. V. Kozachynskyi: conceptualization, software,
investigation, validation, resources, data curation, visualization, writing -
review and editing. E. Esche: conceptualization, investigation, resources,
writing - review and editing, supervision, project administration, funding
acquisition. M. Bortz: conceptualization, writing - review and editing,
supervision, project administration, funding acquisition.

## References

  * [1] J.A.E. Andersson, J. Gillis, G. Horn, J.B. Rawlings, M. Diehl: _CasADi – A software framework for nonlinear optimization and optimal control._ Math. Progr. Comp. 11 (2019), 1-36
  * [2] A.C. Atkinson, A.N. Donev, R.D. Tobias: Optimum experimental design, with SAS. Oxford University Press (2007)
  * [3] C. Atwood: _Sequences converging to D-optimal designs of experiments._ Ann. Statist. 1 (1973), 342-352
  * [4] H. Bandemer, A. Bellmann, W. Jung, K. Richter: Optimale Versuchsplanung. 2nd edition, Akademie-Verlag (1976)
  * [5] D.M. Bates, D.G. Watts: _Relative curvature measures of nonlinearity._ J. R. Statist. Soc. B 42 (1980), 1-25
  * [6] D.M. Bates, D.G. Watts: Nonlinear regression analysis and its applications. Wiley (1988)
  * [7] D. Böhning: _A vertex-exchange method in D-optimal design theory._ Metrika 33 (1986), 337-347
  * [8] G.E.P. Box, J.S. Hunter, W.G. Hunter: Statistics for experimenters. Design, innovation, and discovery. 2nd edition, Wiley (2005)
  * [9] H. Chernoff: _Locally optimal designs for estimating parameters._ Ann. Math. Stat. 24 (1953), 586-602
  * [10] H. Chernoff: _Sequential design of experiments._ Ann. Math. Stat. 30 (1959), 755-770
  * [11] B.P.M. Duarte, A.C. Atkinson, J.F.O. Granjo, N.M.C. Oliveira: _Optimal design of experiments for implicit models._ J. Ann. Statist. Assoc. 117 (2020), 1-38
  * [12] B.P.M. Duarte, A.C. Atkinson, J.F.O. Granjo, N.M.C. Oliveira: _A model-based framework assisting the design of vapor-liquid equilibrium experimental plans._ Comp. Chem. Eng. 145 (2021), 107168
  * [13] I.S. Duff: _MA57 – a code for the solution of sparse symmetric definite and indefinite systems._ ACM Transactions on Math. Software 30 (2004), 118-144
  * [14] V.V. Fedorov: Theory of optimal experiments. Academic Press (1972)
  * [15] V.V. Fedorov, S.L. Leonov: Optimal design for nonlinear response models. Chapman & Hall (2014)
  * [16] R.A. Fisher: The design of experiments. 8th edition, Hafner (1966)
  * [17] G. Franceschini, S. Machietto: _Model-based design of experiments for parameter precision: state of the art._ Chem. Eng. Sci. 63 (2008), 4846-4872
  * [18] G.H. Golub, C.F. Van Loan: Matrix Computations. 3rd edition, Johns Hopkins University Press (1996)
  * [19] D.C. Hamilton, D.G. Watts: _A quadratic design criterion for precise estimation in nonlinear regression models._ Technometrics 27 (1985), 241-250
  * [20] R. Harman, L. Filová and P. Richtárik: _A randomized exchange algorithm for computing optimal approximate designs of experiments._ J. Am. Stat. Assoc. 115 (2020), 348-361
  * [21] HSL, a collection of Fortran codes for large-scale scientific computation. See https://www.hsl.rl.ac.uk
  * [22] J. Jacq, L. Asselineau: _Binary liquid-liquid equilibria, multiple solutions for the NRTL equation._ Fluid Phase Equilib. 14 (1983), 185-192
  * [23] J. Kiefer, J. Wolfowitz: _Optimum designs in regression problems._ Ann. Math. Stat. 30 (1959), 271-294
  * [24] J. Kiefer: _General equivalence theory for optimum designs (approximate theory)._ Ann. Stat. 2 (1974), 849-879
  * [25] K. Kôno: _Optimum design for quadratic regression on $k$-cube._ Mem. Fac. Sci., Kyushu Univ., Ser. A 16 (1962), 114-122
  * [26] D.C. Montgomery: Design and analysis of experiments. 8th edition, Wiley (2009)
  * [27] A.R.G. Mukkula, R. Paulen: _Optimal experimental design in nonlinear parameter estimation with exact confidence regions._ J. Proc. Control 83 (2019), 187-195
  * [28] V.V. Nalimov, T.I. Golikova, N.G. Mikeshina: _On practical use of the concept of D-optimality._ Technometrics 12 (1970), 799-812
  * [29] J. Nocedal, S.J. Wright: Numerical optimization. 2nd edition, Springer (2006)
  * [30] L. Pronzato, A. Pázman: Design of experiments in nonlinear models. Asymptotic normality, optimality criteria and small-sample properties. Springer (2013)
  * [31] J.M. Prausnitz, R.M. Lichtenthaler, E.G. de Azevedo: Molecular thermodynamics of fluid-phase equilibria. 3rd edition, Prentice Hall (1999)
  * [32] F. Pukelsheim: Optimal design of experiments. SIAM Classics of Applied Mathematics (2006)
  * [33] F. Pukelsheim, S. Rieder: _Efficient rounding of approximate designs._ Biometrika 79 (1992), 763-770
  * [34] H. Renon, J.M. Prausnitz: _Local compositions in thermodynamic excess functions for liquid mixtures._ AIChE J. 14 (1968), 135-144
  * [35] J. Schmid, P. Seufert, M. Bortz: _Adaptive discretization algorithms for locally optimal experimental design._ In preparation
  * [36] G.A.F. Seber C.J. Wild: Nonlinear regression. Wiley (1989)
  * [37] P. Seufert, J. Schwientek, M. Bortz: _Model-based design of experiments for high-dimensional inputs supported by machine-learning methods._ Processes 9 (2021), 508
  * [38] S.D. Silvey, D.H. Titterington, B. Torsney: _An algorithm for optimal designs on a design space._ Comm. Statist. Theory Methods 14 (1978), 1379-1389
  * [39] S.D. Silvey: Optimal design. An introduction to the theory for parameter estimation. Chapman & Hall (1980)
  * [40] J.M. Smith, H.C. Van Ness, M.M. Abbott, M.T. Swihart: Introduction to chemical engineering thermodynamics. 8th edition, McGraw-Hill (2018)
  * [41] D. Tassios: _The number of roots in the NRTL and LEMF equations and the effect on their performance._ Ind. Eng. Chem. Process Des. Dev. 18 (1979), 182-186
  * [42] C. Vanaret, P. Seufert, J. Schwientek, G. Karpov, G. Ryzhakov, I. Oseledets, N. Asprion, M. Bortz: _Two-phase approaches to optimal model-based design of experiments: how many experiments and which ones?_ Comp. Chem. Eng. 146 (2021), 107218
  * [43] J.-P. Vila, J.-P. Gauchi: _Optimal designs based on exact confidence regions for parameter estimation of a nonlinear regression model._ J. Stat. Plann. Inference 137 (2007), 2935–2953.
  * [44] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, İ. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, A.E. Quintero, C.R. Harris, A.M. Archibald, A. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors: _SciPy 1.0: Fundamental algorithms for scientific computing in Python._ Nature Methods 17 (2020), 261-272
  * [45] A. Wächter, L.T. Biegler: _On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming._ Math. Progr. 106 (2005), 25-57
  * [46] J. Werner, T. Seidel, R. Jafar, R. Heese, H. Hasse, M. Bortz: _Multiplicities in thermodynamic activity coefficients._ AIChE J. (2023), 18251
  * [47] H.P. Wynn: _The sequential generation of D-Optimal experimental designs._ Ann. Math. Statist. 41 (1970), 1655-1664
  * [48] M. Yang, S. Biedermann, E. Tang: _On optimal designs for nonlinear models: a general and efficient algorithm._ J. Am. Stat. Assoc. 108 (2013), 1411-1420
  * [49] Y. Yu: _D-optimal designs via a cocktail algorithm._ Stat. Comp. 21 (2011), 475-481