Optimal Experiment Design in Nonlinear Parameter Estimation with Exact Confidence Regions

TL;DR

This paper develops a new optimal experiment design method for nonlinear systems that explicitly considers exact confidence regions, using novel approximation techniques to improve computational efficiency and accuracy.

Contribution

It introduces a framework for OED based on exact CRs with ellipsoidal approximations, advancing beyond linearized approaches.

Findings

Exact CR-based OED improves measurement informativeness.

Ellipsoidal approximation reduces computational complexity.

Comparison shows advantages over linearized methods.

Abstract

A model-based optimal experiment design (OED) of nonlinear systems is studied. OED represents a methodology for optimizing the geometry of the parametric joint-confidence regions (CRs), which are obtained in an a posteriori analysis of the least-squares parameter estimates. The optimal design is achieved by using the available (experimental) degrees of freedom such that more informative measurements are obtained. Unlike the commonly used approaches, which base the OED procedure upon the linearized CRs, we explore a path where we explicitly consider the exact CRs in the OED framework. We use a methodology for a finite parametrization of the exact CRs within the OED problem and we introduce a novel approximation technique of the exact CRs using inner- and outer-approximating ellipsoids as a computationally less demanding alternative. The employed techniques give the OED problem as a finite-dimensional mathematical program of bilevel nature. We use two small-scale illustrative case studies to study various OED criteria and compare the resulting optimal designs with the commonly used linearization-based approach. We also assess the performance of two simple heuristic numerical schemes for bilevel optimization within the studied problems.

Tables

Table 1. Table 1 : Optimal designs ( 𝒰 ∗ superscript 𝒰 \mathcal{U}^{*} ) as identified for classical, ellipsoidal and the exact OED using N = { 4 , 5 } 𝑁 4 5 N=\{4,5\} and the values of objective function of exact designs ( ϕ ​ ( ⋅ ) italic-ϕ ⋅ \phi(\cdot) ) evaluated at the identified optimal designs ( 𝒰 ∗ superscript 𝒰 \mathcal{U}^{*} ). In the case of D design, ϕ ​ ( 𝒰 ∗ ) := ϕ ^ D ​ ( 𝒰 ∗ ) assign italic-ϕ superscript 𝒰 subscript ^ italic-ϕ 𝐷 superscript 𝒰 \phi(\mathcal{U}^{*}):=\hat{\phi}_{D}(\mathcal{U}^{*}) .

	Design	$N$	Solution (${𝒰}^{\ast }$)	$\varphi \left({𝒰}^{\ast }\right)$
Classical OED	A	4	$\{1.69,1.69,20,20\}$	1.610
	A	5	$\{1.77,1.77,20,20,20\}$	0.940
	D	4	$\{2,2,20,20\}$	0.425
	D	5	$\{2,2,20,20,20\}$	0.155
	E	4	$\{1.61,20,20,20\}$	1.016
	E	5	$\{1.75,20,20,20,20\}$	0.365
Ellips.	D	4	$\{1.42,1.42,20,20\}$	0.414
	D	5	$\{1.69,1.69,19.99,19.99,20\}$	0.154
Exact OED	A	4	$\{1.37,1.37,20,20\}$	1.585
	A	5	$\{1.60,1.60,20,20,20\}$	0.938
	D	4	$\{1.62,1.62,20,20\}$	0.409
	D	5	$\{1.81,1.82,1.83,19.99,19.99\}$	0.154
	E	4	$\{1.04,1.04,20,20\}$	0.974
	E	5	$\{1.22,1.23,20,20,20\}$	0.322

Table 2. Table 2 : Optimal designs ( 𝒰 ∗ superscript 𝒰 \mathcal{U}^{*} ) as identified for classical, ellipsoidal and the exact OED using N = { 2 , 3 , 4 } 𝑁 2 3 4 N=\{2,3,4\} and the values of objective function of exact designs ( ϕ ​ ( ⋅ ) italic-ϕ ⋅ \phi(\cdot) ) evaluated at the identified optimal designs ( 𝒰 ∗ superscript 𝒰 \mathcal{U}^{*} ). In the case of D design, ϕ ​ ( 𝒰 ∗ ) := ϕ ^ D ​ ( 𝒰 ∗ ) assign italic-ϕ superscript 𝒰 subscript ^ italic-ϕ 𝐷 superscript 𝒰 \phi(\mathcal{U}^{*}):=\hat{\phi}_{D}(\mathcal{U}^{*}) .

	Design	$N$	Solution (${𝒰}^{\ast }$)	$\varphi \left({𝒰}^{\ast }\right)$
Classical OED	A	2	$\{1.91,10\}$	1.666
	A	3	$\{1.86,1.86,10\}$	1.151
	A	4	$\{1.81,1.81,1.81,10\}$	0.974
	D	2	$\{2,10\}$	0.386
	D	3	$\{2,2,10\}$	0.231
	D	4	$\{2,2,10,10\}$	0.148
	E	2	$\{1.90,10\}$	1.225
	E	3	$\{1.82,1.82,10\}$	0.520
	E	4	$\{1.74,1.74,1.74,10\}$	0.341
Ellipsoid	D	2	$\{1.70,10\}$	0.363
	D	3	$\{1.73,1.73,10\}$	0.219
	D	4	$\{1.82,1.82,10,10\}$	0.144
Exact OED	A	2	$\{1.63,10\}$	1.584
	A	3	$\{1.67,1.67,10\}$	1.132
	A	4	$\{1.66,1.66,1.67,10\}$	0.966
	D	2	$\{1.61,10\}$	0.344
	D	3	$\{1.65,1.66,10\}$	0.218
	D	4	$\{1.74,1.77,10,10\}$	0.144
	E	2	$\{1.62,10\}$	1.094
	E	3	$\{1.63,1.63,10\}$	0.497
	E	4	$\{1.59,1.59,1.59,10\}$	0.331

Equations

$$\quad\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{l}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau})=\mathbf{Q}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau})\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},$$

$$\displaystyle\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}}$$

$$\displaystyle J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})$$

$$(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}-\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})^{T}\mathbb{FIM}(\mathcal{U})(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}-\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}},$$

$$\mathbb{FIM}(\mathcal{U})=\sum_{\tau=\tau_{1}}^{\tau_{N}}\mathbf{Q}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau})^{T}\text{diag}^{-1}(\sigma_{1}^{2},\ldots,\sigma_{n_{y}}^{2})\mathbf{Q}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle s^{2}:=\frac{J(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})}{N-n_{p}},$$

$$\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}}\in\underset{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}}{\operatorname{arg}\!\operatorname{min}}\;J(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\sum_{i=1}^{n_{y}}\sum_{\tau=\tau_{1}}^{\tau_{N}}(y_{i,m}(\tau)-\hat{{y}}_{i}({\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}},\tau))^{2}.$$

$$\displaystyle(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}-\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})^{T}\left(\sum_{\tau=\tau_{1}}^{\tau_{N}}\mathbf{Q}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau})^{T}\left(s^{2}\mathbb{I}\right)^{-1}\mathbf{Q}(\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau})\right)(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}-\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})$$

$$\le n_p{\mathcal{F}}_{n_p,N-n_p,\alpha },$$

$$\displaystyle\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{nl}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})-J_{w}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}},$$

$$\displaystyle J(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})-J(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq n_{p}s^{2}\mathcal{F}_{n_{p},N-n_{p},\alpha}.$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}\!:=\!\left\{\!\!\left(\!\!\begin{array}[]{c}p_{1}^{L}\\ p_{2}^{1,L}\\ \vdots\\ p_{n_{p}}^{1,L}\end{array}\!\!\right)\!\!,\!\left(\!\!\begin{array}[]{c}p_{1}^{U}\\ p_{2}^{1,U}\\ \vdots\\ p_{n_{p}}^{1,U}\end{array}\!\!\right)\!\!,\left(\!\!\begin{array}[]{c}p_{1}^{2,L}\\ p_{2}^{L}\\ \vdots\\ p_{n_{p}}^{2,L}\end{array}\!\!\right)\!,\cdots\!\>,\!\left(\!\!\begin{array}[]{c}p_{1}^{n_{p},U}\\ p_{2}^{n_{p},U}\\ \vdots\\ p_{n_{p}}^{U}\end{array}\!\!\right)\!\!\right\},$$

$$\displaystyle\phi_{A}(\mathcal{U}):=\max_{\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}}$$

$$\forall j\in \{1,\dots ,2n_p\},\forall \tau \in \{{\tau }_1,\dots ,{\tau }_N\}:$$

$$\displaystyle\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{nl}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(\tau)=\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}},\tau),$$

$$\displaystyle J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j})-J_{w}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}}.$$

$$\min_{\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}],\,\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}\phi_{A}(\mathcal{U}),$$

$$\displaystyle\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\tilde{\phi}_{\text{A}}(\mathbb{FIM}):=\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\text{trace}(\mathbb{FIM}^{-1})$$

$$\displaystyle\qquad\text{s.t. }\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{l}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),\ \forall\tau\in\{\tau_{1},\dots,\tau_{N}\}.$$

$$\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\phi_{D}(\mathcal{U}):=\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\idotsint\limits_{P_{w}(\mathcal{U})}\text{d}p_{1}\dots\text{d}p_{n_{p}}.$$

$$\displaystyle\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\hat{\phi}_{D}(\mathcal{U}):=\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\sum_{\forall i\in\mathcal{I}_{\Pi}}\delta_{i}$$

$$\text{s.t. }{\delta }_i=\{\begin{array}{ll}1,& \text{if }{\Pi }_i\in P_w\\ 0,& \text{otherwise}\end{array}$$

$$\Pi =\{p_1^L,p_1^L+ϵ,p_1^L+2ϵ,\dots ,p_1^U\}\times$$

$$\{p_2^L,p_2^L+ϵ,p_2^L+2ϵ,\dots ,p_2^U\}\times \dots$$

$$\times \{p_{n_p}^L,p_{n_p}^L+ϵ,p_{n_p}^L+2ϵ,\dots ,p_{n_p}^U\},$$

$$\displaystyle\qquad\max_{\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}}\sum_{j=1}^{n_{p}}p_{j}^{U}-p_{j}^{L},$$

$$\text{s.t. }\forall j\in \{1,\dots ,2n_p\},\forall \tau \in \{{\tau }_1,\dots ,{\tau }_N\}:$$

$$\displaystyle\qquad\quad\quad\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{nl}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle\qquad\quad\quad\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(\tau)=\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}},\tau),$$

$$\displaystyle\qquad\quad\quad J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}_{j})-J_{w}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}},$$

$$\displaystyle\min_{\stackrel{{\scriptstyle\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}\in[\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{L},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}^{U}]}}{{\forall\tau\in\{\tau_{1},\dots,\tau_{N}\}}}}\frac{\text{det}(k_{out}\mathbb{FIM}^{-1})}{k_{out}}+\frac{\text{det}(k_{in}\mathbb{FIM}^{-1})}{k_{in}}$$

$$\displaystyle\text{s.t. }\max_{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{out},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{in},k_{out},k_{in}}k_{out}-k_{in}$$

$$\text{s.t. }\forall \tau \in \{{\tau }_1,\dots ,{\tau }_N\}:$$

$$\displaystyle\qquad\quad\quad\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{out},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{nl}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{out},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle\qquad\quad\quad\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{in},\tau)=\mathchoice{\mbox{\boldmath$\displaystyle F$}}{\mbox{\boldmath$\textstyle F$}}{\mbox{\boldmath$\scriptstyle F$}}{\mbox{\boldmath$\scriptscriptstyle F$}}_{nl}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{in},\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}_{\tau}),$$

$$\displaystyle\qquad\quad\quad\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(\tau)=\hat{\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}},\tau),$$

$$\displaystyle\qquad\quad\quad J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{out})-J_{w}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}},$$

$$\displaystyle\qquad\quad\quad J_{w}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}_{in})-J_{w}(\hat{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})\leq\chi^{2}_{\alpha,n_{p}},$$

Full text

††journal: Journal of Process Control

1 Introduction

At present, advanced industrial engineering and management strive for resource- and energy-efficient design and operation of systems, plants, and processes. Here a use of the model-based techniques is a leading paradigm. The employed models, whether mechanistic or data-based, include a finite number of parameters, whose values are related to the particular natural and system-wide phenomena and are thus commonly only known to belong to some interval or unknown completely. Therefore, as-precise-as-possible determination of the unknown (uncertain) model parameters is crucial for the deployment of the effective model-based solutions.

In the world, where the sensing technology becomes virtually present everywhere, the simplest way of obtaining the parameter values is to conduct a series of observations, experiments, during which the measurements of some quantities (output variables) are gathered. Subsequently, an estimation procedure is employed to find parameter values such that the model outputs match the observed data. As the measurements from a real plant are corrupted with some (random) noise, the uncertain model parameters cannot be identified precisely. This a posteriori uncertainty can be represented by a joint-confidence region (CR) that encompasses all the likely parameter estimates, given the probability distribution of the measurement noise.

The parametric uncertainty can be reduced by performing an experiment that appropriately sets the plant into an operating region, where more informative measurements are obtained. A way how to identify such experiment consists in performing an optimal experiment design (OED) [1, 2]. OED identifies the best experiment in terms of initial conditions, control inputs, sampling times and/or locations of measurement devices. The model-based OED problem is usually formulated as a mathematical program, where a certain measure of the CR, such as volume, is minimized.

Some well-established methods exist for the OED for linear systems [3], wherein CR is an ellipsoid [4], so the task of its optimal shaping is greatly simplified. For nonlinear systems, the most common, yet approximate, approach is to resort to a system linearization and a use of the linear techniques. Beale [5] presented methodology for assessment of system nonlinearity in this respect. Other approaches, covered in [6] for convex problems, are of more or less approximate character and mostly rely on the convexity properties of the OED problem. Bates and Watts [7] presented a framework based on local nonlinear curvature properties, which is also an approximate technique. The use of a finite parametrization of the exact CRs [8, 9, 10] is a relatively recent subject of study.

In this contribution, we study the framework for OED of nonlinear systems that is based on explicit consideration of the exact CRs. We also present a computationally simpler variant that is based on simultaneous inner- and outer-approximations of the CR by ellipsoids. The preliminary findings of this work were presented in the conference contribution [11]. We study various common design criteria and compare the performance of the presented techniques with the linearized OED.

We organized the paper as follows. The concepts of linear and nonlinear parameter estimation and construction of CRs are discussed first. Next, the formulation is presented of the experiment design criteria that directly use exact CRs (exact designs). We also present two simple heuristic numerical approaches taken from literature to solve the arising bilevel optimization problems. Finally results of two illustrative case studies are presented and discussed.

2 Parameter Estimation Problem

2.1 Mathematical Model

In this paper, a mathematical model of a system is represented by

[TABLE]

with \hat{\mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}} as $n_{y}$ measured variables, \mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}_{\tau} as $n_{u}$ degrees of freedom and \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}} as $n_{p}$ uncertain parameters. Here $\tau$ represents an ordinal number of the data point taken in one or more experiments. Function \mathchoice{\mbox{\boldmath\displaystyle F}}{\mbox{\boldmath\textstyle F}}{\mbox{\boldmath\scriptstyle F}}{\mbox{\boldmath\scriptscriptstyle F}}:\mathbb{R}^{n_{p}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{y}} is twice continuously differentiable mapping. Throughout the paper, we resort to a representation (1), which considers the system model as static and explicit. However, the presented findings can straightforwardly be extended to dynamic and implicit models.

We will assume throughout the paper that the model is not over-parametrized and that the parameters are identifiable. We consider that, upon the realization of an experiment or several experiments, $N$ instances are gathered of $n_{y}$-dimensional vector of plant measurements \mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}_{m} and are subsequently used for the estimation of unknown parameters. Throughout the paper, we will assume the Gaussian noise to be corrupting the measurements. In the following subsections, existing frameworks are presented for the identification of the unknown parameters and the corresponding exact CRs for both the linear and the nonlinear parameter estimation problems.

2.2 Linear Parameter Estimation

Parameter estimation with a linear model involved is a well-studied topic in the literature [4]. Assume a mathematical model with a mapping function \mathchoice{\mbox{\boldmath\displaystyle F}}{\mbox{\boldmath\textstyle F}}{\mbox{\boldmath\scriptstyle F}}{\mbox{\boldmath\scriptscriptstyle F}}_{l} of the form

[TABLE]

where \mathbf{Q}(\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}_{\tau}) is a so-called regressors matrix with appropriate dimensions.

Under the assumption of uncorrelated and normally distributed measurement noise with known standard deviation vector $\textstyle\sigma$, the maximum-likelihood estimate is found via the least-squares estimation as

[TABLE]

with

[TABLE]

The CR of parameter estimates is then given by an ellipsoid [4]

[TABLE]

where $\mathbb{FIM}$ is the so-called Fisher information matrix,

[TABLE]

\mathcal{U}:=(\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}_{\tau_{1}}^{T},\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}_{\tau_{2}}^{T},\dots,\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}_{\tau_{N}}^{T})^{T}, and $\chi^{2}_{\alpha,n_{p}}$ represents the upper $\alpha$ quantile of the chi-squared statistical distribution with $n_{p}$ degrees of freedom.

If the variance of the measurement noise is unknown, the covariance matrix $\text{diag}(\sigma_{1}^{2},\ldots,\sigma_{n_{y}}^{2})$ is normally approximated by $s^{2}\mathbb{I}$ with

[TABLE]

where the expected (most-likely) value of parameters \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}} is identified by solving

[TABLE]

The joint-confidence ellipsoid is then given by

[TABLE]

where $\mathbb{FIM}$ is replaced by its corresponding approximation and $\mathcal{F}$ represents the upper $\alpha$ quantile of the Fisher distribution with $n_{p}$ and $N-n_{p}$ degrees of freedom in the numerator and in the denominator, respectively.

2.3 Nonlinear Parameter Estimation

Given a static nonlinear mathematical model

[TABLE]

one can identify the (exact) CR dependent upon the availability of information about the variance of the measurement noise. If the variance is known, the exact CR is given by [4]

[TABLE]

while, when the variance of the measurement noise is unknown, the exact CR is given by [4]

[TABLE]

At this point we can define the sets P_{w}:=\{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}\,|\,\text{Eq.}~{}\eqref{eq:NLPE_with_variance}\} and P:=\{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}\,|\,\text{Eq.}~{}\eqref{eq:NLPE_without_variance}\}. Unlike in the linear parameter estimation, the CR does not generally have a shape of an ellipsoid due to nonlinearity.

3 Optimal Experiment Design

We present a methodology for OED for both linear and nonlinear parameter estimation problems. We will assume that the CR is connected. For disjoint exact CRs, which result from non-identifiability issues, it is normally suggested to perform a re-parameterization of the model [7]. We will also assume that an estimate \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}} is available. The final assumptions which are inherent to the standard experiment design techniques is that there exists no structural plant-model mismatch and that the expected realization of the measurement noise is 0. In turn this results in \mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}_{m}(\tau)=\hat{\mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}}(\hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}},\tau),\,\forall\tau.

Several design criteria are proposed in the literature [3] such as A, D, E, Modified E, V, Q, M and so on. Each of these designs aims to tune a particular property of the confidence region. We will focus our study on the most used criteria, i.e., A, D, and E, but other design criteria might be considered as well using the ideas presented herein.

3.1 A-optimal design

The idea behind the A design criterion is to minimize the perimeter of the box that encloses the exact CR [3], i.e., to minimize the sum of projections of the CR on the parameter axes. This idea is sketched in Fig. 1, where the shaded set represents a CR. The enclosing box is given by the four anchor points, marked as squares ( ).

For a general CR, one can identify $2n_{p}$ anchor points

[TABLE]

where each point gives a lower or an upper limit of the value of a particular parameter in the exact CR. The anchor points can be identified by solving the following optimization problem

[TABLE]

Note that the problem of identifying the orthotopic enclosure of the CR is formulated for the case when the measurement-noise variance is known using the expression for the CR (11). This would simply be exchanged with the expression (12) in case when the variance is unknown.

Note that the problem (30) is separable and highly structured. On the other hand it is non-convex in general and the number of its optimization variables ($2n_{p}^{2}$) grows quadratically with the number of parameters. This means that identification of an orthotope might get challenging for the state-of-the-art solvers and for the high-dimensional problems.

The A-optimal experiment design (A design) can be identified by solving

[TABLE]

which is a special case of a bilevel program. The bounds \mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}^{L} and \mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}^{U} represent the lower and upper limits of the experimental degrees of freedom.

In the linear case, the A design is identified by

[TABLE]

In case of unknown variance of measurement noise, the approximate $\mathbb{FIM}$ from (9) will be used. We denote this linearization-based approach as classical in this study.

3.2 D-optimal design

The D design aims to find the experimental conditions such that the exact CR would have minimum volume. The D-optimal design problem then reads as

[TABLE]

As there is no finite-dimensional parameterization of the set $P_{w}(\mathcal{U})$ available in general, it is very hard in general to evaluate the volume integral. We propose to use a gridding-based approach, where the grid is evaluated inside the aforementioned orthotopic enclosure of the CR. The proposed optimization problem may be written as

[TABLE]

where $\epsilon>0$ is the tuning parameter that determines the accuracy of the approximation and $\mathcal{I}_{\Pi}$ is the index set of $\Pi$. This approach for approximating the volume of the CR is illustrated in Fig. 1 as a grid in the shaded set.

In principle, the identification of the orthotopic enclosure of the CR can be removed and the problem (34) can be modified to a single-level mathematical program. Nonetheless the optimization problem (34) is non-smooth due to the presence of indicator function (34b) and thus it can get very challenging and computationally highly expensive, especially in higher dimensions. An alternative approach can be exploited by approximation of the volume using the semi-algebraic sets [12]. Approaches to inner approximation of the CR based on an orthotope and on successive SDP approximations are presented by Streif et al. [13]. We propose a simpler approximation here, which uses an idea similar to the Löwner-John ellipsoids [14]. We construct the inner- and outer-approximating ellipsoids, which are the scaled counterparts of the linearized CRs. The proposed approximate D design is found by

[TABLE]

where \mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}_{out} and \mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}_{in} are intersection points between outer-/inner-approximating ellipsoids and the exact CR. The scaling factors $k_{out}$ and $k_{in}$ express the magnitude of deviation of the outer- and, respectively, inner-approximating ellipsoid from the linearized CR. The weighting in the cost function is then introduced to penalize the contribution of the most deviating ellipsoid. This prevents the design procedure from concentrating on shaping the ellipsoid that is potentially a very loose approximation of the CR and in practice avoids numerical and irregularity problems. Hence that the proposed problem also scales well, i.e., linearly, w.r.t. the number of the parameters as the lower-level problem optimizes $2n_{p}+2$ variables. The proposed approximate D-optimal design is therefore computationally a less expensive problem when compared to the exact D-optimal design using the gridding-based approach. We will denote the proposed approximate D-optimal design approach as the ellipsoidal D design. The idea behind this method, slightly modified, could also be used for an approximate A design but we do not explore this path in the present study explicitly.

Note also that, if the CR can be expressed exactly as (5), the proposed optimization problem boils down to the classical D design where one solves

[TABLE]

3.3 E-optimal design

Objective of the E design is to minimize the Euclidean distance (\|\mathchoice{\mbox{\boldmath\displaystyle\varphi}}{\mbox{\boldmath\textstyle\varphi}}{\mbox{\boldmath\scriptstyle\varphi}}{\mbox{\boldmath\scriptscriptstyle\varphi}}_{1}-\mathchoice{\mbox{\boldmath\displaystyle\varphi}}{\mbox{\boldmath\textstyle\varphi}}{\mbox{\boldmath\scriptstyle\varphi}}{\mbox{\boldmath\scriptscriptstyle\varphi}}_{2}\|_{2}) between the two points ( in Fig. 1) that belong to the CR and their Euclidean distance is maximal. This can be expressed as

[TABLE]

The E design is also known as a decorrelating design as it aims at finding the experimental conditions such that the CR is as spherical as possible. This criterion is illustrated in Fig. 1. It is noteworthy that the lower-level problem of (37) scales linearly w.r.t. the number of parameters as it optimizes $2n_{p}$ variables.

In the linear case the (classical) E design is identified by

[TABLE]

where $\lambda_{i}(\cdot)$ is the $i$-th eigenvalue of a matrix.

4 Numerical Implementation

In this section we discuss the possible ways to solve the proposed optimization problems. We will exploit BARON [15] in this work in order to guarantee global optimality of the classical OED problems (32), (36), and (38). Special attention is devoted to the presented bilevel programs as the classical OED problems are single-level optimization problems and can straightforwardly be solved using a nonlinear program solver.

We present two simple heuristic approaches taken from literature that can be used to solve the presented bilevel problems, which can be generalized in the form

[TABLE]

A special care has to be taken w.r.t. the non-convexity of lower-level problem. Its global optimum has to be identified in order to guarantee feasibility of the upper-level problem [16].

4.1 Nested approach

The following nested approach is inspired by the formulation proposed in [17, 18] for solving dynamic optimization problems and in [19] for solving a coordination control algorithm using a price-driven coordination technique. The nested approach splits the bilevel optimization problem into a lower-level optimization problem

[TABLE]

that is solved for a given \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1}^{*} using a global solver (e.g. BARON) and a upper-level optimization problem

[TABLE]

that can then be solved for a given \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2}^{*} using a local solver. We use IPOPT [20] as a local solver in this work.

The individual optimization problems are interconnected by the copy variables \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1}^{*} and \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2}^{*} that are exchanged between the problems. The problems are solved repeatedly and the convergence of the nested approach is claimed once the consecutive values of the copy variables satisfy \|\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}^{*}_{1,k+1}-\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}^{*}_{1,k}\|\approx 0 and \|\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}^{*}_{2,k+1}-\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}^{*}_{2,k}\|\approx 0, where $k$ is an iteration counter.

If a gradient-based solver is used to determine the local optimum of the problem (41), the objective and constraint gradients are to be supplied. An approach from [21] can be used in this respect with \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{k}:=(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1,k}^{T},\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2,k}^{T})^{T}

[TABLE]

where $L$ represents the Lagrangian of the lower-level problem and \mathchoice{\mbox{\boldmath\displaystyle\nu}}{\mbox{\boldmath\textstyle\nu}}{\mbox{\boldmath\scriptstyle\nu}}{\mbox{\boldmath\scriptscriptstyle\nu}}:=(\mathchoice{\mbox{\boldmath\displaystyle\nu}}{\mbox{\boldmath\textstyle\nu}}{\mbox{\boldmath\scriptstyle\nu}}{\mbox{\boldmath\scriptscriptstyle\nu}}_{E}^{T},\mathchoice{\mbox{\boldmath\displaystyle\nu}}{\mbox{\boldmath\textstyle\nu}}{\mbox{\boldmath\scriptstyle\nu}}{\mbox{\boldmath\scriptscriptstyle\nu}}_{I,i}^{T})^{T},\forall i\in\mathcal{I}_{A} is the vector of multipliers corresponding to the equality and active inequality constraints \mathchoice{\mbox{\boldmath\displaystyle h}}{\mbox{\boldmath\textstyle h}}{\mbox{\boldmath\scriptstyle h}}{\mbox{\boldmath\scriptscriptstyle h}}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1}^{*},\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2}):=(\mathchoice{\mbox{\boldmath\displaystyle h}}{\mbox{\boldmath\textstyle h}}{\mbox{\boldmath\scriptstyle h}}{\mbox{\boldmath\scriptscriptstyle h}}_{E}^{T}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1}^{*},\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2}),\mathchoice{\mbox{\boldmath\displaystyle h}}{\mbox{\boldmath\textstyle h}}{\mbox{\boldmath\scriptstyle h}}{\mbox{\boldmath\scriptscriptstyle h}}_{I,i}^{T}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{1}^{*},\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{2}))^{T},\forall i\in\mathcal{I}_{A} of the lower-level problem and $\mathcal{I}_{A}$ is an index set of the active inequality constraints. It is not guaranteed that the nested approach always converge to a local minimum [18]. Instead it may sometimes converge to a local maximum or a saddle point. The obtained solution can be verified by evaluating the necessary and sufficient conditions for optimality.

4.2 KKT-based approach

Another heuristic approach for solving a bilevel optimization problem is to reformulate the lower-level problem using the KKT conditions [22, 23, 24]. The reformulated problem reads as

[TABLE]

where $\mathcal{I}_{I}$ is the index set of the inequality constraints of (39). As discussed above, the lower-level problem has to be solved to global optimality, which is not guaranteed by satisfaction of the KKT conditions. Reaching of global optimum of the lower-level problem has to be assured upon convergence in order to guarantee feasibility of the lower-level problem and thus a local optimum of the bilevel program. The feasibility test can be performed by solving (40) globally or by gridding or by set inversion [25] techniques with a subsequent comparison of obtained values for variables of lower-level problem.

We note that there are other approaches that can be employed to solve the problem (39). The solution methods proposed by Dutta et al. in [26] and by Dempe et al. in [27] assume a convex inner level optimization problem. Bard et al., Dempe et al. and Mitsos et al. in [28], [24], and [29] proposed solution methods considering a nonconvex inner level optimization problem. It is generally well known that there is a close connection between bilevel problem and semi-infinite programming (SIP) as discussed in [30]. Cutting-plane SIP algorithm is proposed by Blankenship and Falk [31], branch and bound algorithm by Bard and Moore [32] or double penalty function method by Ishizuka and Aiyoshi [33]. Recently, Walz et al. [34] presented a global SIP algorithm proposed in [35] that could be used in the context of optimal experiment design. Reference therein give a complete picture about the use of SIP algorithms for the solution of bilevel programs.

Also various stochastic methods, such as genetic algorithms, simulated annealing, etc., could be used in principle, where these might be especially interesting for the D design problem (34) because of its non-smoothness. We only exploit the described nested and KKT-based approaches in this study.

5 Case Studies

We apply the presented methodologies for finding OED for two small-scale illustrative case studies. The employed models are in the form of explicit step responses of linear time-invariant dynamic systems and the optimal experiment design should reveal the best sampling instants. We denote the designs that are based on the exact CRs (problems (31), (34), and (37)) as exact designs. The OED problems are solved for $2\sigma$-confidence level ($\alpha=0.9545$) using the following approaches:

Classical OED problems are solved globally. 2. 2.

The designs based on exact CRs and ellipsoidal D design are solved using nested approach where the lower-level problem is solved globally and the upper-level problem is solved using a local solver. 3. 3.

In order to study numerical efficacy of different algorithms, the A-optimal design is solved by KKT-based approach globally.

5.1 Case Study 1

The mathematical model for biological oxygen demand (BOD) [7] is considered. The cumulative BOD of microorganisms at incubation time $u_{\tau}$ is given by

[TABLE]

which can also be interpreted as a step response of the first-order linear time-invariant system with static gain $p_{1}$ and time constant $1/p_{2}$. At this point it can be observed that $p_{1}$ enters the output equation linearly while $p_{2}$ enters nonlinearly.

The true values for the parameters $p_{1}$ and $p_{2}$ are, respectively, 2.5 and 0.5. These are also considered as expected least-squares estimates \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}} for all the studied OED problems. The measurements $y_{m}(\tau)$ are assumed to be corrupted by a zero-mean Gaussian noise with the standard deviation 0.1. We will assume here a case where variance or the standard deviation of the measurement noise is unknown. The exact CRs are then defined by (12). Additionally, we consider $J(\hat{p})=0$. The tolerance $\epsilon$ for the exact D design is set to $5\times 10^{-3}$.

The optimal sampling times $\mathcal{U}^{*}$ with $N=\{4,5\}$ for the classical and exact A-, D- and E-optimal designs and the ellipsoidal D design are reported in Tab. 1. The values of the objective function for exact designs $\phi(\cdot)$ are evaluated at the identified optimal solution $\mathcal{U}^{*}$ for each design with $N=\{4,5\}$. For all the designs, the exact OED has a superior performance when compared to the linearized design. The $\mathcal{U}^{*}$ as identified by the classical and the exact OED contain multiple common repetitive measurements at $u_{\tau}^{*}=20$. This agreement between the designs can be attributed to the linear entry of $p_{1}$ into the model $\hat{y}(\tau)$. It can also further be reasoned by an obvious fact; that one can infer the steady-state gain irrespective of the value of time constant closer to the steady state. The classical and the exact A OED had same number of repetitions of $u_{\tau}^{*}=20$ using $N=\{4,5\}$, however this is not the case for D-optimal design where only using $N=4$ the number of repetitions match. In case of the E-optimal design, the classical OED identified $\mathcal{U}^{*}$ in which $u_{\tau}^{*}=20$ is repeated once more when compared to the solution identified by the exact OED, which again points to the linear decorrelation between $p_{1}$ and $p_{2}$ that classical E design tries to achieve. The performance of the proposed D-optimal design based on the inner- and the outer-approximation ellipsoids is better when compared to the classical D design for all values of $N$. It achieves a relatively small loss in performance when compared to the exact D design. This suggests that the proposed ellipsoidal D design is an interesting framework to perform approximate D design as compared to linearization-based alternative.

The exact CRs for the classical ( ) and the exact ( ) A OED for four measurements are compared in Fig. 2. The orthotopes enclosing the exact CRs are plotted using the anchor points $\textstyle\pi$ represented by . It is clear that the orthotope that encloses the exact CR identified by the classical A design ( ) has a larger perimeter when one compares it with the orthotope identified by the exact A design. The reason for this can be found when looking at the linearized CRs for both designs (dashed ellipsoids). The linearized ellipsoid clearly does not approximate the exact CRs well, where, as it can be expected, the approximation is looser w.r.t. $p_{2}$ that enters nonlinearly in output equation (44). It is an interesting observation that the presented linearized CRs are very similar to each other, despite the fact that they give significantly different exact CRs.

Figure 3 shows the exact CRs for the classical ( ), ellipsoidal ( ) and the exact ( ) D designs using $N=4$. The linearized CRs (dashed ellipsoids) are again very similar to each other. Among them the ellipsoid from the classical design ( ) has the smallest volume, as might be expected. However, the exact CR for the classical design is the largest one (see Tab. 1), which again comes from disregarding of nonlinearity of the output equation in $p_{2}$ by the classical design.

We also present the inner- and the outer-approximation ellipsoids for all the three OED approaches by dash-dotted and dotted ellipsoids, respectively. The corresponding intersection points are marked by and respectively in Fig. 3. Here we can observe the benefits of weighting introduced in the objective function of the problem (35). While looking at the sizes of outer-approximating ellipsoids (especially the one constructed for exact D design), it might appear reasonable to minimize the volume of the over-approximating ellipsoid as a good approximation of the exact D design. This would correspond to setting $1/k_{in}\to 0$ while solving the problem (35). We have explored this path in our earlier study [11], but the obtained design results were unsatisfactory, since the over-approximation by an ellipsoid might become very loose. The proposed ellipsoidal D design therefore balances out the concentration on the size of the ellipsoid and the appropriateness of the over-approximation by an ellipsoid. It is clearly visible that the exact CR for the $\mathcal{U}^{*}$ identified by the proposed ellipsoidal OED has much smaller volume than the exact CR identified by the classical approach and, at the same time, it is very close to the optimal exact OED.

The exact CRs for the classical ( ) and the exact ( ) E designs are compared in Fig. 4 using $N=4$ measurements. in Fig. 4 mark \mathchoice{\mbox{\boldmath\displaystyle\varphi}}{\mbox{\boldmath\textstyle\varphi}}{\mbox{\boldmath\scriptstyle\varphi}}{\mbox{\boldmath\scriptscriptstyle\varphi}}_{1} and \mathchoice{\mbox{\boldmath\displaystyle\varphi}}{\mbox{\boldmath\textstyle\varphi}}{\mbox{\boldmath\scriptstyle\varphi}}{\mbox{\boldmath\scriptscriptstyle\varphi}}_{2} (see (37)) obtained for the classical and the exact E-optimal designs. In this case, unlike for the previous designs, we observe a major discrepancy in the orientation between the linearized ellipsoids obtained for the classical and exact designs, despite the largest semi-axes, which correspond to the largest eigenvalues (see (38)) are very similar. Again this is attributed to the nonlinearity in $p_{2}$. The resulting difference of distances between the most distant points that belong to $P$ is significant, however, among the two designs. Similar behavior can be seen for the case with $N=5$ (see Tab. 1).

Next, we study the performance of the obtained designs against a number of simulated experiments. The aim here is to evaluate robustness of the designs against random realization of noise that would be present in the real experiment. We exclude the dependence on the least-squares estimate here, i.e., we will use the nominal values for \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}}. Such dependence is the subject of study for robust OED, which is not in the scope here. We simulated 1,000 experiments with each studied design using the obtained optimal incubation (measurement) times $\mathcal{U}^{*}$ and we corrupted the measurements y_{m}(\tau):=\hat{y}({\hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}}},\mathchoice{\mbox{\boldmath\displaystyle\tau}}{\mbox{\boldmath\textstyle\tau}}{\mbox{\boldmath\scriptstyle\tau}}{\mbox{\boldmath\scriptscriptstyle\tau}})+e with a Gaussian noise $e$ of standard deviation 0.1.

Figure 5 shows the mean and the variance of the objective value of the exact A, D and E designs. In this plot, we also include the nominal values (when noise-free measurements are gathered) of the different designs using dashed lines. With respect to the mean values of obtained $\phi_{A}$, $\hat{\phi}_{D}$ and $\phi_{E}$, it is confirmed that exact OED is the best option on average, despite the mean values do not match the expected nominal values of the design criteria. Regarding the obtained variances, it is noteworthy that the classical design exhibits the strongest sensitivity to noise as it can be concluded from the magnitude of the variances and thus appears to be the worst option. This again underpins the importance of consideration of the nonlinearity in the OED and it can be documented by comparing the worst-case value of the exact E design w.r.t. nominal mean obtained for the linearized design (right plot). The last interesting observation is that the variance of $\hat{\phi}_{D}$ obtained with ellipsoidal design is slightly smaller than the variance of the exact OED. This might point to the approximation error introduced (tolerance $\epsilon$ in (34)) in calculating $\hat{\phi}_{D}$, which is however not severe as the mean and worst-case values follow the expected trend.

5.2 Case study 2

Here we consider a problem where the system output can be modeled as

[TABLE]

which can also be interpreted as a step response of the second-order linear time-invariant system with a double pole at $-p_{2}$ and a zero at $p_{1}$. The corresponding transfer function can be given by

[TABLE]

Clearly the constant $b_{0}$ is a parameter that influences the steady-state gain of the system and we will assume it, for simplicity, to be known $b_{0}=-4$.

Notice that both $p_{1}$ and $p_{2}$ enter the output equation nonlinearly, so this problem can be considered as more challenging and even greater discrepancy might be expected between linearized and exact designs. The true values of the parameters $p_{1}$ and $p_{2}$ are $0.5$ and $1.0$, respectively, and are equal to \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}}. The measurements $y_{m}(\tau)$ are assumed to be corrupted with a zero-mean Gaussian noise with known standard deviation of $0.4$. For a fair comparison of the proposed framework, we assume J_{w}(\hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}})=0. The tolerance ($\epsilon$) for exact D design is $7.5\times 10^{-2}$.

The classical and the proposed OED frameworks are applied to identify the optimal sampling times $\mathcal{U}^{*}$, where $u_{\tau}^{*}\in[0,10]$. The previously discussed OED problems are solved for $N=\{2,3,4\}$ and $\alpha=0.9545$ ($2\sigma$-CR) with the same numerical techniques as in case study 1. The results are presented in Tab. 2. The trends regarding the performance of the different designs are the same as described in case study 1. We observe in this case a more significant benefit of employing an ellipsoidal D design, which, for $N=\{3,4\}$ almost reaches the performance of exact OED and is superior to classical OED.

In Fig. 6, we compare the exact and linearized CRs obtained for the classical ( ) and the exact ( ) A design using $N=2$. We can see that the linearized CR captures well the orientation of the exact CR but due to nonlinearity the approximation is relatively poor. Similarly to the first case study we can observe that, despite very similar orientation of the linearized confidence regions, significant benefits of the exact design over the linearized one are obtained. Unlike in case study 1, we obtain reduction in the range on both parametric axes, which is caused by the nonlinearity of the output equation w.r.t. both parameters.

Figure 7 shows the resulting CRs for D design criterion. This shows the reason for the good performance of the proposed ellipsoidal technique, which is able to tackle the nonlinearity of the CR far better than the linearization-based design.

The E-optimal designs for the classical and the exact OED with $N=2$ are compared in Fig. 8. Despite the fact that the linearized CRs show great similarity and they capture the orientation of the exact CR, the exact design tackles the nonlinearity far better and shows clear benefits w.r.t. the linearization-based counterpart.

Robustness of the obtained designs was tested against the random realization of the measurement noise for $N=4$ in the same way as in the previous case study. It is clear again that exact designs outperform the classical OED, which reaches the largest variances and the inferior means. The performance of the ellipsoidal D design is practically the same the performance of exact D design. In comparison with with the case study 1, we observe larger variance of the exact D design, which we attribute to the higher nonlinearity.

5.3 Discussion

We studied cases where the CR is found for 2$\sigma$ confidence. For a greater confidence, the CR increases in size and nonlinearity affects it stronger. That is why even bigger differences can be expected between classical and exact OED and even greater benefits can be obtained by using exact design (see [11]).

Regarding the computational efficacy of the different studied problems, it must be clearly pointed out that there exists a high sensitivity of model outputs w.r.t. the model parameters in both case studies due to the presence of highly nonlinear exponential terms. The presence of highly nonlinear terms and the need for inversion of the Fisher information matrix to formulate the objective function in classical OED prohibit BARON from closing the optimality gap and thus it returns locally optimal and sub-optimal solutions unless the optimization problem is properly initialized. However, the classical OED can be solved very efficiently, in few minutes on standard desktop workstation if initialized properly. The solution times for exact designs followed the expectations that result from the aforementioned complexity analysis (see Section 3). On average the solution procedure for exact A, D, and E design using nested approach took less than 10 min, $\approx$6 h, and 15–30 min, respectively. This shows that the optimal exact A and E designs can be obtained with practically the same effort as in the case of the classical design for the small-scale problems. The exact A design procedure scales quadratically in $n_{p}$ so it can get much more time-consuming in higher dimensions. We note that the reduction of CPU time obtained using ellipsoidal D design w.r.t. to exact D design was two-fold ($\approx$3 h), which is, on one hand, a considerable time saving but, on the other hand, it puts the potential user in question, whether the benefits prevail over the costs. We note, for completeness, that the KKT-based approach applied to problem of exact A design required the solution time of approximately one hour, which makes this approach clearly inferior.

In problems with small number of samples, it might be problematic to identify approximate (experimental) variance or to satisfy the asymptotic properties under which the CRs are defined. In this case, one might think of different approaches to experiment design. One such approach might be to relax the assumption of the presence of a white Gaussian noise in the measurements. This might in turn lead to set-membership estimation approach, also commonly known as guaranteed parameter estimation. A step in the direction of experimental design in set-membership context was taken in [36] and in the recent studies [34, 37, 38].

6 Conclusions

In this paper, exact and linearization-based methods were presented for the optimal experiment design of a nonlinear parameter estimation problem. The ellipsoidal method is proposed as a computationally less demanding counterpart to the exact D design, which is a computationally demanding problem since it requires a good approximation for the volume of a set. Two simple heuristic numerical methods are used here to solve the corresponding optimization problems, which are of bilevel nature. The OED framework is tested upon two illustrative small-scale nonlinear case studies, where the benefits of the exact design are showcased. The proposed ellipsoidal technique is shown to perform very well. Despite this study treated the case when the system model describes a static system in an explicit form, the methodology is straightforwardly applicable to dynamic systems and implicit model equations. An interesting direction for the future work lies, on one hand, in increasing the efficiency of the solution of the bilevel programs and, on the other hand, in the study of robust OED that relaxes the assumption of known (expected) least-squares estimates \hat{\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}}, which might be relevant in practical tasks.

Acknowledgments

The research leading to these results has received funding from the European Commission under grant agreement number 291458 (ERC Advanced Investigator Grant MOBOCON). RP acknowledges the contribution of the European Commission under the project GuEst (grant agreement No. 790017), of the Slovak Research and Development Agency under the project APVV 15-0007, of the Scientific GrantAgency of the Slovak Republic under the grant 1/0004/17, and of the Research & Development Operational Programme for the project University Scientific Park STU in Bratislava, ITMS 26240220084, supported by the Research 7 Development Operational Programme funded by the ERDF.