Dynamic Real-time Optimization of Batch Processes using Pontryagin's Minimum Principle and Set-membership Adaptation

TL;DR

This paper introduces a real-time optimization method for batch processes that combines Pontryagin's Minimum Principle with set-membership adaptation to handle model mismatch efficiently.

Contribution

It proposes a novel approach that reduces computational load by integrating reachable sets and online parameter estimation in a predictive control framework.

Findings

Effective handling of model uncertainty using reachable sets.

Reduced computational burden compared to traditional re-optimization.

Successful application to batch membrane separation processes.

Abstract

This paper studies a dynamic real-time optimization in the context of model-based time-optimal operation of batch processes under parametric model mismatch. In order to tackle the model-mismatch issue, a receding-horizon policy is usually followed with frequent re-optimization. The main problem addressed in this study is the high computational burden that is usually required by such schemes. We propose an approach that uses parameterized conditions of optimality in the adaptive predictive-control fashion. The uncertainty in the model predictions is treated explicitly using reachable sets that are projected into the optimality conditions. Adaptation of model parameters is performed online using set-membership estimation. A class of batch membrane separation processes is in the scope of the presented applications, where the benefits of the presented approach are outlined.

Tables

$t$	time [h]
$𝒙$	vector of state variables
$𝒑$	vector of model parameters
$u$	manipulated variable
$𝒥$	objective functional
$F_0$	constant-in-control part of the Lagrange term
$F_u$	multiplier of linear-in-control part of the Lagrange term
${𝒇}_0$	constant-in-control term of the model
${𝒇}_u$	multiplier of linear-in-control term of the model
$𝑷$	interval box
$H$	Hamiltonian (function)
$H_0$	constant-in-control term of the Hamiltonian
$H_u$	multiplier of linear-in-control term of the Hamiltonian
$\mu$	Lagrange multiplier of bound on manipulated variable
$𝝂$	Lagrange multiplier of final conditions
$𝝀(\cdot )$	vector of adjoint variables
$S$	state-dependent switching function
$𝝅$	vector parameterizing the optimal control strategy
$𝒈$	measurement function
$𝒚$	vector of predicted plant outputs
${𝒚}_p\left(t\right)$	vector of true outputs of the plant
${𝒚}_m\left(t\right)$	vector of measured outputs
$𝝈$	magnitude of measurement noise
mid	mid-point of the interval box
$\epsilon$	user-defined tolerance
$A$	membrane area [m2]
$V$	volume of the processed solution
$c_1$	macro-solute concentration
$c_2$	micro-solute concentration
$q$	permeate flux, flow rate through the membrane
${\gamma }_1$	mass-transfer coefficient
${\gamma }_2$	limiting concentration of the macro-solute
${\gamma }_3$	dimensionless non-ideality factor
$𝜸$	parameters of the original permeate-flux model
$p_1$, $p_2$, $p_3$	parameters of the reparametrized permeate-flux model
$𝒰(a,b)$	uniform distribution bounded by $a$ and $b$
Subscripts
0	initial, control-independent
$u$	linear-in-control part
$k$	sampling instant of the plant
f	final
$s$	singular
Superscripts
$L$	lower bound
$U$	upper bound
opt	optimal
nom	nominal
adapt	adaptive
rob	robust

Equations

$$\underset{u(t),t_{\text{f}}}{\min }$$

$$\displaystyle\dot{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})+\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})u(t),$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(0)=\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{0},\quad\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t_{\text{f}},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{\text{f}},$$

$$u(t)\in [u^L,u^U],$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}(t)=\mathchoice{\mbox{\boldmath$\displaystyle g$}}{\mbox{\boldmath$\textstyle g$}}{\mbox{\boldmath$\scriptstyle g$}}{\mbox{\boldmath$\scriptscriptstyle g$}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}}),$$

$$\displaystyle H:=\mu^{L}(u^{L}-u)+\mu^{U}(u-u^{U})+\underbrace{F_{0}+\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}^{T}\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{0}}_{H_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t,\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$}})}\!+\underbrace{\left(F_{u}+\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}^{T}\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{u}\right)}_{H_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t,\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$}})}u,$$

$$\displaystyle\dot{\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=-\frac{\partial H}{\partial\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\quad\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t_{\text{f}},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle\nu$}}{\mbox{\boldmath$\textstyle\nu$}}{\mbox{\boldmath$\scriptstyle\nu$}}{\mbox{\boldmath$\scriptscriptstyle\nu$}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),$$

$$\displaystyle\frac{\partial H}{\partial u}:=H_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t,\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$}})-\mu^{L}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})+\mu^{U}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})$$

$$\displaystyle H(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t,\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$}},u(t),\mu^{L}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mu^{U}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}))$$

$$\displaystyle H_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\mathchoice{\mbox{\boldmath$\displaystyle\lambda$}}{\mbox{\boldmath$\textstyle\lambda$}}{\mbox{\boldmath$\scriptstyle\lambda$}}{\mbox{\boldmath$\scriptscriptstyle\lambda$}}(t,\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$}})=0,\quad\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t_{\text{f}},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})-\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{\text{f}}$$

$$\displaystyle\frac{\text{d}S}{\text{d}t}=\frac{\partial S}{\partial\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}^{T}}\frac{\text{d}\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}{\text{d}t}=$$

$$\displaystyle\quad\Rightarrow\quad u_{\text{s}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})=$$

$$\displaystyle u^{\ast}(t,\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}}):=\begin{cases}u^{L},&t\in[0,t_{1}),\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})>0,\\ u^{U},&t\in[0,t_{1}),\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})<0,\\ u_{\text{s}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}}),&t\in[t_{1},t_{2}),\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})=0,\\ u^{L},&t\in[t_{2},t_{\text{f}}],\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{\text{f}},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})<0,\\ u^{U},&t\in[t_{2},t_{\text{f}}],\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{\text{f}},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})>0,\end{cases}$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}_{\text{f}}=\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t_{2},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})+\int_{t_{2}}^{t_{\text{f}}}\!\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})+\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})u^{\ast}(t,\mathchoice{\mbox{\boldmath$\displaystyle\pi$}}{\mbox{\boldmath$\textstyle\pi$}}{\mbox{\boldmath$\scriptstyle\pi$}}{\mbox{\boldmath$\scriptscriptstyle\pi$}})\,\text{d}t,$$

$$\displaystyle|\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(t)-\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{p}(t)|\leq\mathchoice{\mbox{\boldmath$\displaystyle\sigma$}}{\mbox{\boldmath$\textstyle\sigma$}}{\mbox{\boldmath$\scriptstyle\sigma$}}{\mbox{\boldmath$\scriptscriptstyle\sigma$}},$$

$$\displaystyle|\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(t)-\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}(t)|\leq\mathchoice{\mbox{\boldmath$\displaystyle\sigma$}}{\mbox{\boldmath$\textstyle\sigma$}}{\mbox{\boldmath$\scriptstyle\sigma$}}{\mbox{\boldmath$\scriptscriptstyle\sigma$}}.$$

$$\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}}_{k}\subseteq\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}}_{k-1}\subseteq\dots\subseteq\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}}_{1}\subseteq\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}}_{0},$$

$$p_{k,j}^L/$$

$$\displaystyle\dot{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})+\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})u(t),\ \forall t\in[0,t_{k}],$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(0,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle h$}}{\mbox{\boldmath$\textstyle h$}}{\mbox{\boldmath$\scriptstyle h$}}{\mbox{\boldmath$\scriptscriptstyle h$}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),$$

$$\displaystyle\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}(t_{i},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle g$}}{\mbox{\boldmath$\textstyle g$}}{\mbox{\boldmath$\scriptstyle g$}}{\mbox{\boldmath$\scriptscriptstyle g$}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t_{i},\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$}}),\ \forall i\in\{1,\dots,k\},$$

$$\displaystyle-\mathchoice{\mbox{\boldmath$\displaystyle\sigma$}}{\mbox{\boldmath$\textstyle\sigma$}}{\mbox{\boldmath$\scriptstyle\sigma$}}{\mbox{\boldmath$\scriptscriptstyle\sigma$}}\leq\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}(t_{i},\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})-\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}_{m}(t_{i})\leq\mathchoice{\mbox{\boldmath$\displaystyle\sigma$}}{\mbox{\boldmath$\textstyle\sigma$}}{\mbox{\boldmath$\scriptstyle\sigma$}}{\mbox{\boldmath$\scriptscriptstyle\sigma$}},\ \forall i\in\{1,\dots,k\},$$

$$\displaystyle t_{l}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})$$

$$\displaystyle u_{s}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})$$

$$T_{1}:=\left\{t_{1}\,\left|\,\begin{array}[]{c}\exists\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}\in\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}},\forall t\in[0,t_{1}]:\\ \mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(0,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle h$}}{\mbox{\boldmath$\textstyle h$}}{\mbox{\boldmath$\scriptstyle h$}}{\mbox{\boldmath$\scriptscriptstyle h$}}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),\ S(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t_{1},\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$}})=0,\\ \dot{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})=\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{0}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})+\mathchoice{\mbox{\boldmath$\displaystyle f$}}{\mbox{\boldmath$\textstyle f$}}{\mbox{\boldmath$\scriptstyle f$}}{\mbox{\boldmath$\scriptscriptstyle f$}}_{u}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}(t,\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$}})u^{L}\end{array}\right.\right\}.$$

$$\displaystyle\min_{\stackrel{{\scriptstyle u_{s}(t,\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})\in U^{\text{opt}}(t),\forall t\in[t_{1}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}),t_{2}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})]}}{{t_{l}\in T_{l},\ \forall l\in\{1,2,\text{f}\}}}}\max_{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}\in\mathchoice{\mbox{\boldmath$\displaystyle P$}}{\mbox{\boldmath$\textstyle P$}}{\mbox{\boldmath$\scriptstyle P$}}{\mbox{\boldmath$\scriptscriptstyle P$}}_{0}}\ \|\mathcal{J}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})-\mathcal{J}(\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}^{\text{nom}}_{0})\|_{2}^{2}\qquad\text{s.t.}\ \eqref{eq:model},\eqref{eq:singular_ctrl_switch},\eqref{eq:fin_cond},$$

$$\displaystyle q(\underbrace{c_{1}(t),c_{2}(t)}_{\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t)},\underbrace{\gamma_{1},\gamma_{2},\gamma_{3}}_{\mathchoice{\mbox{\boldmath$\displaystyle\gamma$}}{\mbox{\boldmath$\textstyle\gamma$}}{\mbox{\boldmath$\scriptstyle\gamma$}}{\mbox{\boldmath$\scriptscriptstyle\gamma$}}})$$

$$=A{\gamma }_1\left[\ln ({\gamma }_2)-\ln (c_1(t))-{\gamma }_3\ln (c_2(t))\right],$$

$$\displaystyle q(\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t),\underbrace{p_{1},p_{2},p_{3}}_{\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}}})$$

$$\underset{t_{\text{f}},u(t)}{\min }{\int }_0^{t_{\text{f}}}1\text{d}t,$$

$$\text{s.t. }{\dot{c}}_1(t)$$

$${\dot{c}}_2(t)$$

$$\displaystyle q(\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t),\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})$$

$$u(t)$$

$$\displaystyle u_{s}(\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t,\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$}})$$

$$\displaystyle S(\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t,\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$}})$$

$$\displaystyle:=\underbrace{p_{1}-p_{2}\ln(c_{1})-p_{3}\,\ln(c_{2})}_{q(\mathchoice{\mbox{\boldmath$\displaystyle c$}}{\mbox{\boldmath$\textstyle c$}}{\mbox{\boldmath$\scriptstyle c$}}{\mbox{\boldmath$\scriptscriptstyle c$}}(t),\mathchoice{\mbox{\boldmath$\displaystyle p$}}{\mbox{\boldmath$\textstyle p$}}{\mbox{\boldmath$\scriptstyle p$}}{\mbox{\boldmath$\scriptscriptstyle p$}})}-p_{2}-p_{3}.$$

Full text

††journal: Computers & Chemical Engineering

1 Introduction

Optimization of operations of batch processes is a rich field of research. One of the main goals is to reduce the variability among the produced batches despite the present uncertainties and disturbances. This problem struck the attention of many research groups [21, 34, 1, 10, 19, 20, 15].

In this paper, we consider a real-time implementation of a control policy under parametric plant-model mismatch that optimizes a batch process by assigning dynamic degrees of freedom such that a certain performance index is optimized. Similar problems were studied in many previous works using on-line or batch-to-batch adaptation of the optimality conditions [11, 10], by mid-course correction [38, 14] or by design of robust controller for tracking the conditions of optimality [21]. Another set of approaches to the problem uses advanced robust strategies in the framework of model predictive control [19]. This paper proposes an adaptation of these approaches to the problem of dynamic real-time optimization of batch processes. This task is not straightforward because if one uses a receding-horizon control strategy, the prediction horizons used need to be quite long, because of the usual presence of terminal constraints, which might compromise the real-time feasibility of the scheme.

We base the presented methodology on the parameterization of the optimal operation using the optimality conditions given by Pontryagin’s minimum principle. As the cost is usually insensitive w.r.t. a precise singular control trajectory [34], the parameterization of the optimal policy makes the real-time decision problem to mainly boil down to identification of switching times of the optimal control policy. Such approach reduces computational burden while allowing for the use of sufficiently long prediction horizons when projecting the parametric uncertainty in controller performance and feasibility, particularly w.r.t. terminal time conditions. Robustness w.r.t. parametric uncertainty is addressed by taking into account the imprecision of parameter estimates, which is projected into the uncertainty of the switching times. In order to improve performance of such a controller, i.e., to reduce conservatism introduced by uncertain switching times, we use on-line parameter estimation. While having the optimal control policy explicitly parameterized in the uncertain parameters, one can tailor the real-time implementation of the optimal operation, e.g., in a way that minimizes the number of on-line calculations.

The novelty of this paper lies foremost in an effective combination of Pontryagin’s minimum principle and set-based techniques (set-membership estimation and reachability analysis). This gives rise to a methodology capable of projecting the propagation of the uncertainty in model parameters into uncertainty in the optimal operation of a plant. Using this methodology, efficient and effective real-time optimization of a plant can be established.

The outline of the paper is as follows. Section 3 present preliminary theoretical knowledge on Pontryagin’s minimum principle [26] and on set-membership estimation [31, 9]. The former is used to parameterize the optimality conditions of the dynamic optimization problem, while the latter technique is used for adaptation of the model parameters based on the measured data along the process run. Next we propose the implementation of the real-time optimization using parameter adaptation. Finally, we present a case study from chemical engineering domain and discuss various aspects of the obtained results.

2 Problem definition

In this paper, we consider a real-time implementation of a control policy that optimizes a process by assigning dynamic degrees of freedom such that a certain performance index is optimized:

[TABLE]

where $t$ is time with $t\in[0,t_{\text{f}}]$, \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}(\cdot) is an $n$-dimensional vector of state variables, $\textstyle p$ is an $m$-dimensional vector of model parameters, $u(t)$ is a (scalar) manipulated variable, $F_{0}(\cdot)$, $F_{u}(\cdot)$, \mathchoice{\mbox{\boldmath\displaystyle f}}{\mbox{\boldmath\textstyle f}}{\mbox{\boldmath\scriptstyle f}}{\mbox{\boldmath\scriptscriptstyle f}}_{0}(\cdot), and \mathchoice{\mbox{\boldmath\displaystyle f}}{\mbox{\boldmath\textstyle f}}{\mbox{\boldmath\scriptstyle f}}{\mbox{\boldmath\scriptscriptstyle f}}_{u}(\cdot) are continuously differentiable functions, \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{0} represents a vector of initial conditions, and \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{\text{f}} are specified final conditions. We note here that an inclusion of multi-input and/or state-constrained cases is a straightforward extension but it is not considered in this study for the sake of simplicity of the presentation. We also note that the specific class of input-affine systems is a suitable representation for a large variety of the controlled systems [13]. For a general nonlinear model, one may use simple manipulations to rearrange the model into input-affine structure [33], which might though increase the number states of the problem. In the domain of chemical engineering, it is, however, very common to encounter input-affine problems [2] (e.g., when the optimized variable is a reactor feed) or to reformulate the model and arrive at the input-affine structure [18].

We will assume that the plant behavior is known qualitatively and that the only source of uncertainty is present in the unknown values of model parameters. Only a prior knowledge is assumed about the parameters, i.e., the true values of the parameters lie in the a priori known interval box \mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}:=[\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}_{0}^{L},\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}_{0}^{U}], where superscripts $L$ and $U$ denote the lower and upper bounds of $\textstyle p$. The nominal realization of the uncertainty will be assumed as \mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}^{\text{nom}}:=\text{mid}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}), where $\text{mid}(\cdot)$ indicates a mid-point of the interval box.

We will also assume that certain measurements are available from the plant. Their corresponding model-based predictions are

[TABLE]

where \mathchoice{\mbox{\boldmath\displaystyle g}}{\mbox{\boldmath\textstyle g}}{\mbox{\boldmath\scriptstyle g}}{\mbox{\boldmath\scriptscriptstyle g}}(\cdot) is a continuously differentiable vector function.

3 Preliminaries

3.1 Conditions for Optimality

Pontryagin’s minimum principle can be used [17, 35, 23, 24] to identify the optimal solution to (1) via enforcing the necessary conditions for minimization of a Hamiltonian

[TABLE]

where \mathchoice{\mbox{\boldmath\displaystyle\lambda}}{\mbox{\boldmath\textstyle\lambda}}{\mbox{\boldmath\scriptstyle\lambda}}{\mbox{\boldmath\scriptscriptstyle\lambda}}(\cdot) is a vector of adjoint variables, which are defined through

[TABLE]

and \mu^{L}(t,\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}), \mu^{U}(t,\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}), and \mathchoice{\mbox{\boldmath\displaystyle\nu}}{\mbox{\boldmath\textstyle\nu}}{\mbox{\boldmath\scriptstyle\nu}}{\mbox{\boldmath\scriptscriptstyle\nu}}(\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}) are the corresponding Lagrange multipliers. The optimality conditions of (1) can then be stated as [35]: $\forall t\in[0,t_{\text{f}}]$,

[TABLE]

The condition $H=0$ arises from the transversality, since the final time is free [26], and from the fact that the optimal Hamiltonian is constant over the whole time horizon, as it is not an explicit function of time. The condition $H_{0}=0$ is the consequence of the former two conditions. Since the Hamiltonian is affine in input (see (3)), the optimal trajectory of control variable is either determined by active input constraints or it evolves inside the feasible region.

Assume that for some point $t$ we have $H_{u}=0$ and $u^{L}<u(t)<u^{U}$. It follows from (5) that the optimal control maintains $H_{u}(\cdot)=0$. Such control is traditionally denoted as singular. Further properties of the singular arc, such as switching conditions or state-feedback control trajectory can be obtained by differentiation of $H_{u}$ with respect to time (sufficiently many times) and by requiring the derivatives to be zero. The time derivatives of $H$ and $H_{0}$ must be equal to zero as well. Earlier results on derivation of optimal control for input-affine systems [35] suggest that it is possible to eliminate \mathchoice{\mbox{\boldmath\displaystyle\lambda}}{\mbox{\boldmath\textstyle\lambda}}{\mbox{\boldmath\scriptstyle\lambda}}{\mbox{\boldmath\scriptscriptstyle\lambda}}(\cdot) from the optimality conditions and thus to arrive at analytical characterization of switching conditions between singular and saturated-control arcs.

As the optimality conditions obtained by the differentiation w.r.t. time are linear in the adjoint variables, the differentiation of $H_{u}$ (or $H_{0}$) can be carried out until it is possible to transform the obtained conditions to a pure state-dependent switching function S(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}(t),\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}). It is usually convenient to use a determinant of the coefficient matrix of the equation system \mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}\mathchoice{\mbox{\boldmath\displaystyle\lambda}}{\mbox{\boldmath\textstyle\lambda}}{\mbox{\boldmath\scriptstyle\lambda}}{\mbox{\boldmath\scriptscriptstyle\lambda}}=\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}} for this. The singular control u_{\text{s}}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}(t),\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}) can be found from differentiation of switching function w.r.t. time as

[TABLE]

The resulting optimal-control policy is then given as a step-wise strategy [24] by

[TABLE]

where \mathchoice{\mbox{\boldmath\displaystyle\pi}}{\mbox{\boldmath\textstyle\pi}}{\mbox{\boldmath\scriptstyle\pi}}{\mbox{\boldmath\scriptscriptstyle\pi}}:=(\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}^{T},t_{1},t_{2},t_{\text{f}})^{T} is the vector that parameterizes the optimal control strategy. Note that the presented optimal-control strategy determines implicitly the switching times $t_{1}$, $t_{2}$ and the terminal time $t_{\text{f}}$ as functions of model parameters $\textstyle p$.

In case that the use of the minimum principle turns out to be too complex (e.g., many derivatives are needed to characterize the solution), a numerical identification of the control arcs [30, 29] or a recently presented parsimonious input parameterization [3, 28] can be used.

3.2 Set-membership estimation

In order to estimate the model parameters, we will make use of plant outputs (measurements), whose predictions are expressed as in (2). We will assume that the true output of the plant \mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}_{p}(t) is corrupted with a (sensor) noise that is bounded with a known magnitude $\textstyle\sigma$. Thus, the measured output \mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}_{m}(t) is such that

[TABLE]

where the absolute value is understood component-wise. In turn, the set-membership constraints for predicted output \mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}(t) apply in the form:

[TABLE]

We are interested in the determination of parametric bounds such that

[TABLE]

where $k$ is the ordinal number of a measurement taken. The parametric bounds can be determined through solution of a series of optimization problems as [12, 37]:

[TABLE]

for given $u(t)$, where $j\in\{1,\dots,n_{p}\}$ indicates the $j^{\text{th}}$ element of a vector.

4 Dynamic real-time optimization

As the optimal control structure is a function of uncertain parameters, the uncertainty should be taken into account when devising a real-time implementation of the optimal control of the plant.

4.1 Projection of parametric uncertainty into solution strategy

Given the structure of the optimal-control policy (9), one can project the parametric uncertainty into uncertainty of the switching times and singular control as (\forall\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}\in\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}})

[TABLE]

This can be achieved either by using some set-theoretic techniques for calculating reachable sets [5] or by sampling approaches. Figure 1 provides an illustration, where the reachable sets are shown over time for the switching function $S(\cdot)$.

Formally, the problem of determination of (15) can be cast as a set-inversion problem [16]. As an example, let us consider that S(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}(0,\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}})>0,\ \forall\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}\in\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}. The interval $T_{1}$ can then be defined as:

[TABLE]

The rest of the uncertain intervals and controls can be defined and determined analogously. Efficient set-inversion techniques exist [25] and can be used herein. The problem might also be reformulated to a bound-determining optimization problem, similarly to the estimation problem in (14). Using such a reformulation, it is also possible to merge the problems of reachability analysis and set-membership estimation and to formulate the reachability-analysis problem directly over the collected data. This way several possible deficiencies (such as those arising from outbounding the parameter set by a box) can be eliminated.

Figure 2 illustrates the parameterization (15) for a simple case, where the singular control is constant. Note that the parameterization reveals time intervals (i.e., [0,t_{1}^{L}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}})], [t_{1}^{U}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}),t_{2}^{L}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}})], and [t_{2}^{U}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}),t_{\text{f}}^{L}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}})]), which are parts of the optimal solution for any realization of uncertain parameters and are thus invariant to the presence of uncertainty.

A particular technical advantage can be exploited for determination of the switching intervals i.e., that the integration in (10) can be done backwards in time from the final condition. As the batch processes exhibit inherently unstable dynamics, their backward integration is stable [4]. Such a feature can readily be exploited by modern reachability analysis approaches for parametric ordinary differential-algebraic equations [36].

4.2 Robust approach to real-time optimization

The result (15), in practice, establishes a parametric solution to the real-time optimization problem. Its implementation can be performed in a robust fashion to determine the parameters of the optimal-control structure that lead to the best performance in the worst case. We can then solve

[TABLE]

for a given \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}(0)=\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}_{0} and \mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}. Here we propose to minimize the variance of the objective w.r.t. nominal scenario under the worst-case realization of \mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}\in\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}, which can also be modified to \|\mathcal{J}(\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}})-\min_{p^{\text{opt}}\in\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}}\mathcal{J}(\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}^{\text{opt}})\|_{2}^{2}. Note that this goal goes in line with the efforts of practical batch process control, where the reduction of the batch-to-batch variability is one of the main targets of decision making.

4.3 Robust adaptive approach to real-time optimization

In order to reduce conservatism of a robust scheme, parameter estimation can be used for exploitation of data gathered along the process run. The employed parameter estimation scheme should take into account the presence of noise in the measurements. Here we propose to use set-membership strategy outlined in Section 3.2.

The problem (17) can then be resolved with the initial state conditions \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}}_{k} and with updated parameter bounds \mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{k} in a shrinking-horizon fashion. Computational efficiency of this real-time optimization scheme can be achieved by exploiting the fact that the re-optimization does not need to be done at each sampling time of the plant (i.e., when new measurements become available) but can be scheduled before a consecutive switching event must be realized. As an example, consider Fig. 2, where one can start the operation on the lower-bound of the input variable and estimation of parameter bounds and re-optimization can be scheduled in the sampling instant of the plant just before minimal value of the time $t_{1}$, t_{1}^{L}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}). The re-optimization with updated bounds on parameters would then update (possibly increase) the value of t_{1}^{L}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}). Further re-optimizations can then follow based on this scheme.

Once the optimal value of the objective function of (17) reaches \|\mathcal{J}(\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}})-\mathcal{J}(\mathchoice{\mbox{\boldmath\displaystyle p}}{\mbox{\boldmath\textstyle p}}{\mbox{\boldmath\scriptstyle p}}{\mbox{\boldmath\scriptscriptstyle p}}_{k}^{\text{nom}})\|_{2}^{2}<\varepsilon, where $\varepsilon>0$ represents user-defined tolerance for the worst-case cost variation, the calculated control actions can be implemented until the end of the batch, e.g., with a feedback scheme [10], until the terminal conditions are met.

Note that because of the switching nature of the optimal control strategy, the proposed problem might show discontinuity when the set of active constraints changes. This can be remedied by the adaptation of continuous-formulation technique presented in [7].

A pseudo-algorithm can be devised at this point to summarize the proposed approach:

Given the problem setup (Eqs. (1) and (2)), identify the solution structure using Pontryagin’s minimum principle (as shown in Section 3.1). 2. 2.

Given \mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}, use reachability analysis to project the uncertainty in the parameters to uncertainty about the solution structure (as shown in Section 4.1). 3. 3.

Apply the optimal policy until the next uncertain switching time and collect the measurements along. 4. 4.

Solve problem (14) to determine new interval box $\textstyle P$ and re-calculate the uncertain solution structure. 5. 5.

If significant reduction in the uncertainty of the switching time is achieved, go to Step 3. Else solve problem (17) to determine the switching times. Apply the first switching in the control, collect the data along and go to Step 4.

5 Case study

We demonstrate the findings of this study on an example of time-optimal control of a batch diafiltration process [6]. This is a membrane-based separation process designed for a simultaneous concentration of valuable products in the liquid solutions (referred to as macro-solutes) and a wash-out of the impurities (referred to as micro-solute).

A simplified scheme of the plant is shown in Fig. 3. After the separated solution with initial volume ($V_{0}$) comprising a macro-solute (high molecular weight component) and a micro-solute (low molecular weight component) with initial concentrations $c_{1,0}$ and $c_{2,0}$, respectively, is transferred to the feed tank, the operation of the process is launched. Solution containing diluant (solvent), micro-solute and macro-solute is taken from the feed tank to the membrane module. The installed membrane is designed in a way to allow passage of micro-solute and to retain macro-solute. Permeate stream then leaves the system with the flow rate $q$ which is specific for a given membrane, operating conditions and is often a function of actual concentrations of separated species. Retentate stream is then introduced back into the feed tank. Once the final conditions, which are the prescribed final concentrations of the species $c_{1,\text{f}}$ and $c_{2,\text{f}}$, are met, the process is terminated and the solution is withdrawn from the system. During the operation, the transmembrane pressure is controlled at a constant value. The temperature of the solution is maintained around a constant value using a heat exchanger (not shown in Fig. 3 for the sake of simplicity). The manipulated variable $u(t)$ is the ratio between fresh water inflow into the tank and the permeate outflow $q$.

In our study, the outflow $q$ is measured at intervals of one minute and its model is given by

[TABLE]

Here the parameters $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ can be related to phenomenological constants; $\gamma_{1}$ stands for the mass-transfer coefficient, $\gamma_{2}$ is the limiting concentration of the macro-solute, and $\gamma_{3}$ is a dimensionless non-ideality factor. This model is proposed in [27] as a generalization of a limiting-flux model, which originates from film-theory of mass transfer [8] and where $\gamma_{3}=0$. We study situations where the permeate flux at the plant obeys either one of these models.

As Eq. (5) suggests, an equivalent re-parameterization of the model is possible, which gives the model linear in parameters $p_{1}$, $p_{2}$, and $p_{3}$. This is convenient for parameter estimation. In this work, we will assume that the concentrations $c_{1}(t)$ and $c_{2}(t)$ can be measured perfectly (i.e., their measurement sensors are noise-free). As the values of state variables are known exactly, the dynamic equations can be eliminated from the problem (14). Hence, the problem of estimating bounds of the model parameters (14) boils down to a problem of linear programming, since dynamics can be excluded and since the re-parameterization of the model yields linear-in-parameters structure of the model. The measurement noise associated with $q$ is assumed to be bounded $\sigma=0.1\mathrm{L/h}$ and the measurements are available each second. In the simulation studies below, the realization of noise will be taken from a uniform distribution $\mathcal{U}(-\sigma,\sigma)$.

The objective is to find $u(t)$, which guarantees the transition from the given initial ($c_{1,0}$ and $c_{2,0}$) to final ($c_{1,\text{f}}$ and $c_{2,\text{f}}$) concentrations in minimum time. This problem can be formulated as:

[TABLE]

The parameters of the problem are $c_{1,0}=50\,\mathrm{g/L}$, $c_{1,\text{f}}=150\,\mathrm{g/L}$, $c_{2,0}=50\,\mathrm{g/L}$, $c_{2,\text{f}}=0.05\,\mathrm{g/L}$, $V_{0}=20\,\mathrm{L}$, and $A=1\,\mathrm{m^{2}}$. Note that the extremal values of $u(t)$ in (19e) stand for a mode with no water addition, when $u(t)=0$ and pure dilution, i.e., a certain amount of water is added at a single time instant, $u(t)=\infty$.

The parameterized optimal control of this process can be identified using Pontryagin’s minimum principle [26] as (9) where the singular control and the respective switching function can be found explicitly [23] as

[TABLE]

The structure of the optimal-control policy clearly reveals that the singular arc condition gives a constant value for the permeate flux (equal to $p_{2}+p_{3}$) and that the singular control is a constant that depends on the value of $\gamma_{3}=p_{3}/p_{2}$. This shows that if one devises a feedback-based real-time optimization scheme, precise estimation of parameters $p_{2}$ and $p_{3}$ is of paramount interest.

For the simulation-based studies on the implementation of the outlined optimal-control policy, we will assume that the nominal values of the parameters are $\gamma_{1}=3\times 10^{-2}\,\mathrm{L/h}$, $\gamma_{2}=1000\,\mathrm{g/L}$, and $\gamma_{3}=0.1$. Similar values of the parameters were observed to validate the model against experimental data in [32]. The uncertainty at the initial point in time (\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}_{0}) will be assumed as $\pm 10\%$ of the nominal values. The true realization of the parameter values will be taken randomly from a uniform distribution $\mathcal{U}(p_{0}^{L},p_{0}^{U})$. As we deal with a time-minimization problem and the sampling time of the plant is 1 second, we naturally select the $\varepsilon=1\,\mathrm{s^{2}}$. For this simple example, the reachability analysis can be performed explicitly using the expressions for switching times provided in [22].

5.1 Plant under limiting-flux conditions

We first study the case when the plant is under limiting-flux conditions, i.e., the flux obeys Eq. (5) with $\gamma_{3}=p_{3}=0$. Based on the values of initial conditions and range of uncertainty in the parameters, the optimal-control policy boils down to three arcs:

Use $u(t)=0$ until $t_{1}$, when $c_{1}(t_{1})=\gamma_{2}/\text{e}$. 2. 2.

Use $u(t)=u_{s}(t)=1$ until $c_{1}(t_{2})/c_{2}(t_{2})=c_{1,\text{f}}/c_{2,\text{f}}$. 3. 3.

Dilute the solution (use $u(t)=\infty$ instantaneously) to arrive at the final concentrations.

Taking into account that we measure both concentrations precisely, the only uncertainty in this case lies in the switching times $t_{1}$, where $t_{1}$ depends on the value of $\gamma_{2}$.

The implementation of the scheme, where one is aware of the true value of $\gamma_{2}$, results in $t_{1}^{\text{opt}}=2.515\,\mathrm{h}$ and $t_{\text{f}}^{\text{opt}}=8.284\,\mathrm{h}$. The worst-case minimization of the batch variability (methodology described in Section 4.2) coincides in this case with the nominal strategy, where one takes $\gamma_{2}=\gamma_{2}^{\text{nom}}$. When applied to the plant, this strategy results in $t_{1}^{\text{rob}}=t_{1}^{\text{nom}}=2.625\,\mathrm{h}$ and $t_{\text{f}}^{\text{rob}}=t_{\text{f}}^{\text{nom}}=8.327\,\mathrm{h}$. The adaptive real-time dynamic optimization (described in Section 4.3) results in $t_{1}^{\text{adapt}}=2.533\,\mathrm{h}$ and $t_{\text{f}}^{\text{adapt}}=8.301\,\mathrm{h}$, which is only a slight improvement compared to the robust (and nominal) strategy.

Figure 4 presents performance of the estimation (in terms of estimated parameter bounds) throughout the run of the batch. It is clear that the bounds on both parameters are dramatically reduced around the time point of 2 h, which precedes the time point $t_{1}^{\text{opt}}$, when the switch in the control input should be executed. The bottom plot of Fig. 4 also shows the evolution of the uncertainty in $t_{1}$, which is projected using interval-based calculations (as discussed in Section 4.1). It should be noted here that the adaptive approach is successful mainly since the applied control input in the first arc coincides with an input that would result from a dynamic optimal-experiment design study. Here, $u(t)=0$ ensures the fastest possible increase of concentration $c_{1}(t)$, which reveals the most informative measurements about $\gamma_{2}$.

Finally, we evaluate a statistical performance of the presented dynamic real-time optimization schemes. This is realized by running 1,000 simulated batches with different true values of parameters. The resulting statistics is shown in Fig. 5. One can clearly notice here that the robust and nominal strategies perform very well on average and even the standard deviations of their performance are not significantly increased compared to the optimal performance. On the other hand, distributions of the final batch times of the nominal and robust strategies have long tails, which points towards the existence of rare cases where the batch time obtained by application of robust dynamic real-time optimization increases significantly compared to the truly optimal solution. This behavior corresponds to the situations, where the first step (with $u(t)=0$) over-concentrates the solution too much so that the subsequent diafiltration step (with $u(t)=1$) requires long time to reach the desired condition ($c_{1}(t_{2})/c_{2}(t_{2})=c_{1,\text{f}}/c_{2,\text{f}}$). This study reveals that the strategy, which uses estimation of parameter bounds, has a clear merit as it does not show this type of inconsistency in the performance (long tails) and it results overall in the batch times very close to the optimal ones.

Drawing a comparison in the computational time, robust and nominal strategies only require a single optimization before the batch starts. Adaptive strategy uses only a single re-optimization (scheduled just before $t_{1}^{\text{nom}}$) in this case, which shows a significant reduction in computational burden w.r.t. receding-horizon strategies. Due to explicit nature of the optimal control strategy, the computation of this re-optimization mostly lies in the estimation step. As the estimation problem can be boiled down to an LP, its solution is available in order of milliseconds using MATLAB’s linprog routine. Due to only a single re-optimization, one can also interpret this scheme as a mid-course correction [39] with optimally timed adaptation.

5.2 Plant under generalized limiting-flux conditions

Based on the values of initial conditions and range of uncertainty in the parameters, the optimal-control policy again boils down to three arcs:

Use $u(t)=0$ until $t_{1}$, when (21) is zero. 2. 2.

Use $u(t)=u_{s}(t)=1/(1+\gamma_{3})$ until $c_{1}(t_{2})/c_{2}(t_{2})=c_{1,\text{f}}/c_{2,\text{f}}$. 3. 3.

Dilute the solution (use $u(t)=\infty$ instantaneously) to arrive at the final concentrations.

The uncertainty is in this case extended even on the value of singular control input and it is clear that one needs good estimates of both the values of $\gamma_{2}$ (which mostly influences the switching time $t_{1}$ as in the previous case) and $\gamma_{3}$ (which influences the quality of the singular control) to achieve a good performance. Knowledge of the a precise value of the parameter $\gamma_{1}$ is of minor importance as this parameter can be factored out of the optimality conditions.

We use the same values of the uncertain parameters $\gamma_{1}$ and $\gamma_{2}$ as in the previous case. The performance of the studied schemes is as follows:

Optimal strategy: $t_{1}^{\text{opt}}=2.501\,\mathrm{h}$, $t_{\text{f}}^{\text{opt}}=9.254\,\mathrm{h}$

• 2.

Adaptive strategy: $t_{1}^{\text{adapt}}=2.510\,\mathrm{h}$ and $t_{\text{f}}^{\text{adapt}}=9.271\,\mathrm{h}$

• 3.

Nominal strategy: $t_{1}^{\text{nom}}=2.561\,\mathrm{h}$ and $t_{\text{f}}^{\text{nom}}=9.277\,\mathrm{h}$

• 4.

Robust strategy: $t_{1}^{\text{rob}}=2.417\,\mathrm{h}$ and $t_{\text{f}}^{\text{rob}}=9.269\,\mathrm{h}$

We can observe a similar differences between the strategies as in the previous case. Robust strategy performs on an acceptable level and even marginally outperforms the nominal and the adaptive strategy, which results from the fact that the plant parameters coincide with the worst-case parameters.

Figure 6 presents performance of the estimation (in terms of estimated parameter bounds) throughout the run of the batch. Similarly to the previous case, the estimation performed in the first control arc helps in determination of the value of the first switching time before the optimal switching instant occurs. Here the determining parameter is $p_{2}$, whose estimation performance was discussed in the previous case and the same conclusions hold here.

When the controller applies $u(t)=0$, the parameter $\gamma_{3}$ (or $p_{3}$) is unidentifiable as the concentration $c_{2}(t)$ remains constant. This can be seen in Fig. 6 as the bounds on $p_{3}$ remain constant from the beginning of the operation until the time when control input is switched to singular. It is also shown that the uncertainty in $p_{3}$ results in a relatively small uncertainty in the value of the singular control, so a precise knowledge of $p_{3}$ is not paramount for the application of the optimal control policy.

When the statistical performance is evaluated, we can first conclude that all the strategies perform almost identically on average. The biggest differences arise when one evaluates the outliers of the distribution of the achieved batch times. It is evident that the nominal strategy achieves the worst performance and that the robust strategy reduces the batch variability to a good extent (given the wide range of the uncertainty). The adaptive strategy is clearly superior here as it reduces the batch variability much further, compared to the robust scheme, and the achieved performance is practically indistinguishable from the truly optimal one.

6 Conclusion

We have presented a methodology for dynamic real-time optimization of batch processes via parameterization of the optimal controller using Pontryagin’s minimum principle. The employed parameterization greatly reduces the computational burden to guarantee feasibility of the operation compared to receding-horizon strategies. In order to address parametric plant-model mismatch issue, we have suggested a robust approach, which consisted in projection of the plant uncertainty into optimality conditions using reachability analysis. This again greatly reduces on-line computational burden as one can exploit the uncertainty in the switching to schedule the on-line re-optimization. As the uncertainty in parameters can greatly affect the optimality of the batch, we have proposed an adaptive scheme that makes use of parameter estimation and, as shown in the case study, can greatly assist in reducing variability in the batch performance subject to parametric uncertainty. The adaptive scheme turned out to be key in reduction of batch-to-batch variability. The future work will consider implementation of the proposed strategy on a laboratory plant.

Acknowledgments

The authors gratefully acknowledge the contribution of the Scientific Grant Agency of the Slovak Republic under the grant 1/0004/17, of the Slovak Research and Development Agency under the project APVV 15-0007 and of the European Commission under the grant 790017 (GuEst). This publication is also a partial result 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. This work was also supported by the funding from Slovak Ministry of Education, Science, Research and Sport under the project STU as the Leader of Digital Coalition 002STU-2-1/2018.

Notation