mintOC - User contributions [en]

Batch reactor

2026-04-16T12:10:55Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nc = 2
|nre = 2
}}

This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.

The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.

== Mathematical formulation ==

The optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{x, T} & x_2(t_f) \\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]
& x(0) &=& (1, 0)^T, \\
& T(t) &\in& [298, 398].
\end{array}
</math>


<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.

== Parameters ==

The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.

== Reference Solutions ==

This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].

The optimal objective value of the problem is <math>x_2(t_f) = 0.611715</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.
</gallery>

== Source Code ==

Model descriptions are available in
* [[:Category:Gekko | GEKKO Python code]] at [[Batch reactor (GEKKO)]]
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]

== References ==
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Oil Shale Pyrolysis

2025-07-24T13:18:33Z

SebastianSager: /* Mathematical formulation */

{{Dimensions
|nd = 4
|nx = 4
|nu = 1
|nw = 0
|nc = 4
|nre = 0
}}

The following problem is an example from the global optimal control literature and was introduced in <bib id="Wen1977" />. The process starts with kerogen and is decomposed into pyrolytic bitumen, oil and gas, and residual carbon. The objective is to maximize the fraction of pyrolytic bitumen. There are 5 reactions including:

<math>A_1 \xrightarrow{k_1} A_2</math>

<math>A_2 \xrightarrow{k_2} A_3</math>

<math>A_1 + A_2 \xrightarrow{k_3} A_2 + A_2</math>

<math>A_1 + A_2 \xrightarrow{k_4} A_3 + A_2</math>

<math>A_1 + A_2 \xrightarrow{k_5} A_4 + A_2</math>

Each reaction is governed by a rate described by:

<math>k_i = k_{i0} \exp{\left(-E_i/RT\right)}, (i=1,2,3,4,5)</math>

== Mathematical formulation ==

<math>
\begin{array}{llll}
\displaystyle \max_{T} & \displaystyle x_2(t_N) \\[1.5ex]
\mbox{s.t.} & \displaystyle \dot{x}_1 &=& -k_1x_1-(k_3+k_4+k_5)x_1x_2\\
& \displaystyle \dot{x}_2 &=& k_1x_1-k_2x_2 + k_3x_1x_2\\
& \displaystyle \dot{x}_3 &=& k_2x_2 + k_4x_1x_2\\
& \displaystyle \dot{x}_4 &=& k_5x_1x_2\\
& \displaystyle k_i &=& a_i e^{-\frac{b_i}{RT}},\quad \forall i\in \{1,\dots,5\} \\[1.5ex]
& \displaystyle t &\in& \left[t_0,t_N\right] \\
& \displaystyle T(t) &\in& \left[698.15K,748.15K\right]\\
& \displaystyle x(t_0) &=& (1,0,0,0)^T\\
\end{array}
</math>

== Parameters ==

{| class="wikitable"
|+State variables
|-
|Symbol
|Initial value (<math>t_0</math>)
|-
|<math>x_1(t_0)</math>
|<math>1</math>
|-
|<math>x_2(t_0)</math>
|<math>0 </math>
|-
|<math>x_3(t_0)</math>
|<math>0 </math>
|-
|<math>x_4(t_0)</math>
|<math>0 </math>
|}

{| class="wikitable"
|+Parameters
|-
|Symbol
|Value
|-
|<math>a_1</math>
|<math>20.3</math>
|-
|<math>a_2</math>
|<math>37.4</math>
|-
|<math>a_3</math>
|<math>33.8</math>
|-
|<math>a_4</math>
|<math>28.2</math>
|-
|<math>a_5</math>
|<math>31.0</math>
|-
|<math>b_1</math>
|<math>\exp{(8.86)}</math>
|-
|<math>b_2</math>
|<math>\exp{(24.25)}</math>
|-
|<math>b_3</math>
|<math>\exp{(23.67)}</math>
|-
|<math>b_4</math>
|<math>\exp{(18.75)}</math>
|-
|<math>b_5</math>
|<math>\exp{(20.7)}</math>
|}

{| class="wikitable"
|+Control variable
|-
|Symbol
|Interval
|-
|<math>T(t)</math>
|[698.15,748.15]
|}

== Solutions ==

* [[:Category:Gekko | GEKKO Python code]] at [[Oil shale pyrolysis (GEKKO)]]

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]

Diels-Alder Reaction Experimental Design

2025-07-24T12:49:19Z

SebastianSager:

The '''Diels-Alder Reaction''' is an organic chemical reaction.
A conjugated diene and a substituted alkene react and form a substituted cyclohexene system.
Stefan Körkel used this model in his PhD thesis to compute optimum experimental designs.

== Model Formulation ==

The reactionkinetics can be modelled by the following differential equation system:

<math>
\begin{array}{rcl}
\dot{n_1}(t) &=& -k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}}, \\
& & \\
\dot{n_2}(t) &=& -k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}}, \\
& & \\
\dot{n_3}(t) &=& \ \ k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}} \\
& & \\
\dot{n_4}(t) &=& 0
\end{array}
</math>

The reaction velocity constant <math>k</math> consists of two parts. One part reflects the non-catalytic and the other the catalytic reaction. The velocity law follows the Arrhenius relation

<math>
k = k_1 \ \cdot \ exp(- \frac{E_1}{R} \ \cdot \ (\frac{1}{T(t)} \ - \ \frac{1}{T_{ref}}) ) \ + \ k_{cat} \ \cdot \ c_{cat} \ \cdot \ exp(-\lambda \ \cdot \ t) \ \cdot \ exp( - \frac{E_{cat}}{R} \ \cdot \ (\ \frac{1}{T(t)} \ - \ \frac{1}{T_{ref}}) )
</math>

Total mass:

<math>
m_{tot} = n_1 \ \cdot \ M_1 \ + \ n_2 \ \cdot \ M_2 \ + \ n_3 \ \cdot \ M_3 \ + \ n_4 \ \cdot \ M_4
</math>

Temperature in Kelvin:

<math>
T(t) = \vartheta (t) + 273
</math>

The ODE system is summarized to:

<math>
\begin{array}{rcl}
\dot{x}(t) &=& f(x(t), u(t), p)
\end{array}
</math>

== Constraints ==

The control variables are constrained with respect to the mass of sample weights (initial mass):

<math>
\begin{array}{cll}
0.1 & \le & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \le 10
\end{array}
</math>

and to the mass of active ingredient content (fraction of active substances):


<math>
\begin{array}{cll}
0.1 & \le & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \le 0.7
\end{array}
</math>

== Optimum Experimental Design Problem ==

The aim is to compute an optimal experimental design <math>\xi = (q,w)</math> which minimizes the uncertainties of the parameters <math>k_1, k_{cat}, E_1, E_{cat}, \lambda</math>. So, we have to solve the following optimum experimental design problem:

<math>
\begin{array}{cll}
\displaystyle \min_{x^i,\ G^i,\ F^i,\ Tc^i,\ n_{a1}^i,\ n_{a2}^i,\ n_{a4}^i,\ c_{kat}^i,\ \vartheta(t)^i} && trace(F^{-1} (t_{end})) \\[1.5ex]
\mbox{s.t.} \\
\dot{x}^i(t) & = & f(x^i(t), u^i(t),p), \\
\\
h(t) & = & \frac{n_3(t) \ \cdot \ M_3}{m_{tot}} \ \cdot \ 100 \\
\\
\dot{G}^i(t) & = & f_x(x^i(t),u^i(t),p)G^i(t) \ + \ f_p(x^i(t),u^i(t),p) \\
\\
\dot{F}(t) & = & \sum\limits_{i=1}^{4} w^i(t) (h^i_x(x^i(t),u^i(t),p)G^i(t))^T (h^i_x(x^i(t),u(t),p)G^i(t)) \\
\\
0.1 & \le & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \\
\\
10 & \ge & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \\
\\
0.1 & \le & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \\
\\
0.7 & \ge & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \\
\\
\vartheta(t) & = & \left\{ \begin{array}{cl} \vartheta_{lo} + 273 & t \in [t_0,2] \\
\vartheta_{lo} + \frac{t-2}{6} ( \vartheta_{up} - \vartheta_{lo} ) + 273 & t \in [2,8] \\
\vartheta_{up} + 273 & t \in [8,t_{end}]
\end{array} \right. \\
& & x \in \mathcal{X},\,u \in \mathcal{U},\, p \in P \\
\dot{z}^i(t) & = & w^i(t) \\
z(0) & = & 0 \\
w^i(t) &\in& [0,1] \\
0 & \le & 4 - z^i(t_f). \\
\end{array}
</math>


{| class="wikitable"
|+State variables
|-
|Name
|Symbol
|Initial value (<math>t_0</math>)
|-
|Molar number 1
|<math>n_1(t)</math>
|<math>n_1(t_0) = n_{a1} </math>
|-
|Molar number 2
|<math>n_2(t)</math>
|<math>n_2(t_0) = n_{a2} </math>
|-
|Molar number 3
|<math>n_3(t)</math>
|<math>n_3(t_0) = 0 </math>
|-
|Solvent
|<math>n_4(t)</math>
|<math>n_4(t_0) = n_{a4} </math>
|}

{| class="wikitable"
|+Constants
|-
|Name
|Symbol
|Value
|-
|Molar Mass
|<math>M_1</math>
|0.1362
|-
|Molar Mass
|<math>M_2</math>
|0.09806
|-
|Molar Mass
|<math>M_3</math>
|0.23426
|-
|Molar Mass
|<math>M_4</math>
|0.236
|-
|Universal gas constant
|<math>R</math>
|8.314
|-
|Reference temperature
|<math>T_{ref}</math>
|293
|-
|St.dev of measurement error
|<math>\sigma</math>
|1
|}

Remember, in an optimum experimental design problem the parameters of the model are fixed. But, we minimize the parameter's uncertainties by optimizing over the control variables and functions.

{| class="wikitable"
|+Fixed parameters
|-
|Name
|Symbol
|Value
|-
|Steric factor
|<math>k_1</math>
|<math>p_1 \cdot 0.01</math>
|-
|Steric factor
|<math>k_{kat}</math>
|<math>p_2 \cdot 0.10</math>
|-
|Activation energie
|<math>E_1</math>
|<math>p_3 \cdot 60000</math>
|-
|Activation energie
|<math>E_{kat}</math>
|<math>p_4 \cdot 40000</math>
|-
|Catalyst deactivation coefficient
|<math>\lambda</math>
|<math>p_5 \cdot 0.25</math>
|}
with <math>p_j = 1, \ j =1, \dots, 5</math>

{| class="wikitable"
|+Optimization/control variables
|-
|Name
|Symbol
|Interval
|Initial value Exp 1
|Initial value Exp 2
|Initial value Exp 3
|Initial value Exp 4
|-
|Initial molar number 1
|<math>n_{a1}</math>
|[0,10.0]
|1.0
|1.0
|1.0
|1.0
|-
|Initial molar number 2
|<math>n_{a2}</math>
|[0,10.0]
|1.0
|1.0
|1.0
|1.0
|-
|Initial molar number 4
|<math>n_{a4}</math>
|[0.4,9.0]
|2.0
|2.0
|2.0
|2.0
|-
|Concentration of the catalyst
|<math>c_{kat}</math>
|[0,10.0]
|0.0
|1.0
|2.0
|3.0
|}

{| class="wikitable"
|+Control function
|-
|Name
|Symbol
|Time interval
|Value interval
|Initial value Exp 1
|Initial value Exp 2
|Initial value Exp 3
|Initial value Exp 4
|-
|Initial molar number 1
|<math>\vartheta(t)</math>
|<math>[t_0,2]</math>
|[20.0,100.0]
|20.0
|60.0
|40.0
|20.0
|-
|Initial molar number 1
|<math>\vartheta(t)</math>
|<math>[2,8]</math>
|[20.0,100.0]
|20.0
|60.0
|40.0
|20.0
|-
|Initial molar number 1
|<math>\vartheta(t)</math>
|<math>[8,t_{end}]</math>
|[20.0,100.0]
|20.0
|60.0
|40.0
|20.0
|}

'''Measurement grid'''

<math>
\begin{array}{llll}
t_0 = 0 & & & \\
t_{end} = 20 & & & \\
t_j = j/3, & j = 1,\dots, 15, & t_j = j - 10, & j = 16, \dots, 20.
\end{array}
</math>

== Source Code ==

* The VPLAN code using [[:Category: VPLAN | VPLAN]] can be found in: [[Diels-Alder Reaction Experimental Design (VPLAN)]]

== References ==

R. T. Morrison and R.N. Boyd. Organic Chemistry. Allyn and Bacon, Inc., 4th edition, 1983
S. Körkel. Numerische Methoden für Optimale Versuchsplanungsprobleme bei nichtlinearen DAE-Modellen.PhD thesis, Universität Heidelberg, Heidelber,2002

[[Category:MIOCP]]
[[Category:Optimum Experimental Design]]
[[Category:ODE model]]

User:SebastianSager

2025-06-05T12:21:49Z

SebastianSager: /* Contact */

__NOTOC__


[[Image:Sager.jpg|thumb]]

== Personal ==
{|
|-
| '''Name''': || Prof. Dr. Sebastian Sager
|-
| '''Affiliation 1''': || [https://mathopt.de Mathematical Algorithmic Optimization], [https://ovgu.de University of Magdeburg]
|-
| '''Affiliation 2''': || [https://www.mpi-magdeburg.mpg.de/oml Optimization and Machine Learning], [https://www.mpi-magdeburg.mpg.de/ Max-Planck Institute Magdeburg]
|-
| '''Background''': || Mathematics
|}

== Contact ==
{|
|-
| '''Address''': || Universitätsplatz 2, 39108 Magdeburg, Germany
|-
| '''Web''': || work group [https://mathopt.de https://mathopt.de]
|-
| '''Email''': || sager -- ovgu.de
|}

== Miscellaneous ==
I started these pages and am the owner of the domain. Hence, if you have any questions or comments concerning mintOC.de, feel free to write me an email.


[[Category:Community]]

External Links

2025-04-25T06:49:28Z

SebastianSager: /* Benchmark Libraries */

= Work Groups =

An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.

* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson
* [https://www.tu-braunschweig.de/mo/team/kirches Mathematical Optimization], led by Christian Kirches
* [http://www.mathematik.tu-clausthal.de/arbeitsgruppen/kontinuierliche-optimierung/ Continuous Optimization], led by Andreas Potschka
* [https://www.mathematik.hu-berlin.de/de/forschung/forschungsgebiete/mathematische-optimierung/Falk-hante/falk-hante-fp Applied Mathematics], led by Falk Hante

Feel free to update this list!

= Benchmark Libraries =

=== Hybrid Systems ===

* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.

=== Optimal Control ===

* The [http://tomdyn.com/ PROPT homepage] (a matlab toolkit for dynamic optimization using collocation) states over 100 test cases from different applications with their results and computation time.
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems.
* The ESA provides a [https://www.esa.int/gsp/ACT/projects/gtop/ test set of global optimization spacecraft trajectory problems] and their best putative solutions.
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization
* The [http://plato.asu.edu/pdecon.html PDECON collection] contains discretized PDE problems in AMPL as part of a [http://plato.asu.edu/bench.html larger benchmark effort] by Hans Mittelmann

=== MINLP ===

* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]
* MINLPlib <bib id="Bussieck2003a" />
* [http://www.gamsworld.org/performance GAMSWORLD]

=== Generic Optimization ===

* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)

== References ==
<biblist />

DOW Experimental Design

2024-12-09T08:37:32Z

SebastianSager: gmol to mol

{{Dimensions
|nd = 1
|nx = 11
|nw = 2
|nc = 4
|nre = 11
}}

The '''DOW Experimental Design problem''' models the OED problem for the parameter estimation problem formulated by the DOW Chemical Co. in 1981. The following formulation is taken from [[#Biegler1986|[1]]].

== Chemical background ==

The chemical species are disguised for proprietary reasons and the desired reaction is given by <math>HA + 2BM \rightarrow AB + HBMH</math>, where <math>AB</math> is the desired product. The reactions are described as follows:

Slow Kinetic Reactions:

<math>
\begin{array}{c}
M^- + BM \underset{k_{-1}}{\leftarrow} \overset{k_1}{\rightarrow} MBM^- \\
A^- + BM \overset{k_2}{\rightarrow} ABM^- \\
M^- + AB \underset{k_{-3}}{\leftarrow} \overset{k_3}{\rightarrow} ABM^-
\end{array}
</math>


Acid-Base Reactions:

<math>
\begin{array}{c}
MBMH \leftarrow K_1 \rightarrow MBM^- + H^+ \\
HA \leftarrow K_2 \rightarrow A^- + H^+ \\
HABM \leftarrow K_3 \rightarrow ABM^- + H^+
\end{array}
</math>

In order to devise a model to account for these reactions, it is
first necessary to distinguish between the overall concentration of
a species and the concentration of its neutral form. Overall concentrations are defined for three components based on neutral
and ionic species

<math>
\begin{align}
\left[HBMH\right] &= \left[ (MBMH)_N \right] + \left[ MBM^- \right] \\
\left[HA\right] &= \left[ (HA)_N \right] + \left[ A^- \right] \\
\left[HABM\right] &= \left[ (HABM)_N \right] + \left[ ABM^- \right]
\end{align}
</math>

Here <math>[ \ ]</math> denotes the concentration of the species in <math>\operatorname{mol}/\operatorname{kg}</math>.
By assuming the rapid acid-base reactions are at equilibrium, the equilibrium constants <math>K_1,K_2,K_3</math> can be defined as

<math>
\begin{align}
K_1 &= \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\
K_2 &= \frac{[A^-][H^+]}{[(HA)_N]} \\
K_3 &= \frac{[ABM^-][H^+]}{[(HABM)_N]}
\end{align}
</math>


The anionic species may then be represented by

<math>
\begin{align}
\left[MBM^-\right] &= \frac{K_1[MBMH]}{K_1 + [H^+]} && \quad (a) \\
\left[A^-\right] &= \frac{K_2[HA]}{K_2 + [H^+]} && \quad (b)\\
\left[ABM^-\right] &= \frac{K_3[HABM]}{K_3 + [H^+]} && \quad (c)
\end{align}
</math>

Material balance equations for the three reactants in the slow
kinetic reactions yield:

<math>
\begin{align}
\frac{d[M^-]}{dt} &= -k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] - k_3 \left[ M^- \right]\left[ AB \right] + k_{-1} \left[ ABM^- \right] && \quad (d)\\
\frac{d[BM]}{dt} &= -k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] - k_2 \left[ A^- \right]\left[ BM \right] && \quad (e)\\
\frac{d[AB]}{dt} &= -k_3 \left[ M^- \right] \left[ AB \right] + k_{-3} \left[ ABM^- \right] && \quad (f)
\end{align}
</math>

From stoichiometry, rate expressions can also be written for
the total species:

<math>
\begin{align}
\frac{d[MBMH]}{dt} &= k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] && \quad (g)\\
\frac{d[HA]}{dt} &= k_2 \left[ A^- \right] \left[ BM \right] && \quad (h)\\
\frac{d[HABM]}{dt} &= k_2 \left[ A^- \right] \left[ BM \right] + k_3 \left[ M^- \right] \left[ AB \right] - k_{-3} \left[ ABM^- \right] && \quad (i)
\end{align}
</math>

An electroneutrality constraint gives the hydrogen ion concentration <math>\left[ H^+ \right]</math> as

<math>
\left[ H^+ \right] + \left[ Q^+ \right] = \left[ M^- \right] + \left[ MBM^- \right] + \left[ A^- \right] + \left[ ABM^- \right] \quad \quad (j)
</math>


Based on similarities of reacting species, we assume

<math> k_3 = k_1, \quad k_{-3} = \frac{1}{2} k_{-1} </math>

With these assumptions, the three rate constants <math>k_1,k_2</math> and <math>k_{-1}</math> must be estimated. Each of these can be fitted with two adjustable model parameters, assuming an Arrhenius temperature dependence. That is

<math>
\begin{array}{c}
k_i = \alpha_i \exp(-E_i /(RT)), \quad i \in \{-1,1,2\}
\end{array}
</math>

Here <math>R \approx 1.98720425864083 \ cal/(K \cdot mol) </math> is the gas constant and <math>T</math> is the reaction temperature in Kelvins. The parameters <math>\alpha_i</math>, given in <math>mol/( kg \cdot h)</math>, represent the pre-exponential factors and the <math>E_i</math>, with unit <math>cal/mol</math>, are the activation energies.

== Mathematical formulation ==

The chemical processes <math>(a)-(j)</math> can be expressed mathematically as six differential equations and four algebraic equations:

<math>
\begin{align}
\dot{y}_1 &= -k_2 y_8 y_2 && \quad (1),(h) \\
\dot{y}_2 &= -k_1 y_6 y_2 + k_{-1} y_{10} - k_2 y_8 y_2 && \quad (2),(e) \\
\dot{y}_3 &= k_2 y_8 y_2 + k_1 y_6 y_4 - \frac{1}{2} k_{-1} y_9 && \quad (3),(i) \\
\dot{y}_4 &= -k_1 y_6 y_4 + \frac{1}{2} k_{-1} y_9 && \quad (4),(f) \\
\dot{y}_5 &= k_1 y_6 y_2 - k_{-1} y_{10} && \quad (5),(g) \\
\dot{y}_6 &= -k_1 (y_6 y_2 + y_6 y_4) + k_{-1} (y_{10} + \frac{1}{2} y_9) && \quad (6),(d) \\
y_7 &= -\left[ Q^+ \right] + y_6 + y_8 + y_9 + y_{10} && \quad (7),(j)\\
y_8 &= \frac{\theta_8 y_1}{\theta_8 + y_7} && \quad (8),(b)\\
y_9 &= \frac{\theta_9 y_3}{\theta_9 + y_7} && \quad (9),(c)\\
y_{10} &= \frac{\theta_7 y_5}{\theta_7 + y_7} && \quad (10),(a)
\end{align}
</math>

Here the letters in parentheses stand for the corresponding chemical process and the quantity <math>\left[Q^+\right]=0.00131</math> is a constant during the reaction.
The nine parameters form the vector

<math>
\theta = (\alpha_1,E_1,\alpha_2,E_2,\alpha_{-1},E_{-1},K_1,K_2,K_3)
</math>

The predicted concentrations form the vector

<math>
y = (HA,BM,HABM,AB,MBMH,M^-,H^+,A^-,ABM^-,MBM^-)
</math>


Let <math>f_k(\cdot)</math> denote the right hand side of equation <math>(k)</math> for <math>k=1,\ldots,6</math>.
We reformulate the last four algebraic equations as differential ones:

<math>
\begin{align}
\dot{y}_7 &= f_7 = f_6 + f_8 + f_9 + f_{10} && \quad (7') \\
\dot{y}_8 &= f_8 = \frac{\theta_8 f_1 \cdot (\theta_8 + y_7) - \theta_8 y_1 f_7}{(\theta_8 + y_7)^2} && \quad (8') \\
\dot{y}_9 &= f_9 = \frac{\theta_9 f_3 \cdot (\theta_9 + y_7) - \theta_9 y_3 f_7}{(\theta_9 + y_7)^2} && \quad (9')\\
\dot{y}_{10} &= f_{10} = \frac{\theta_7 f_5 \cdot (\theta_7 + y_7) - \theta_7 y_5 f_7}{(\theta_7 + y_7)^2} && \quad (10')
\end{align}
</math>

The right hand sides of <math>(1)-(6)</math> and <math>(7')-(9')</math> are summarized as the vector-valued function <math>f</math>.
Moreover, let

<math>
\begin{align}
f_{y,m,n}(\cdot) &= \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,10\} \\
f_{\theta,m,n}(\cdot) &= \frac{\partial f_m(\cdot)}{\partial \theta_n}, \quad m \in \{1,\ldots,10\}; \ n\in \{1,\ldots,9\}
\end{align}
</math>


== Parameters ==

The initial parameter estimates are:
<table style="border-collapse: collapse;" border="1.">
<tr>
<td style="text-align: center; padding:5pt"> <math> \alpha_1 </math> </td>
<td style="text-align: center; padding:5pt"> <math>2.0 \cdot 10^{13} \ mol / (kg \cdot h)</math> </td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>\alpha_2</math></td>
<td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{13} \ mol / (kg \cdot h)</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>\alpha_{-1}</math></td>
<td style="text-align: center; padding:5pt"><math>4.3 \cdot 10^{15} \ mol / (kg \cdot h)</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>E_1</math></td>
<td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>E_2</math></td>
<td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol </math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>E_{-1}</math></td>
<td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>K_1</math></td>
<td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-17} \ mol / kg</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>K_2</math></td>
<td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-11} \ mol / kg</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt"><math>K_3</math></td>
<td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-17} \ mol / kg</math></td>
</tr>
</table>

There are three datasets for different temperatures <math>T</math>, with corresponding starting values
<table style="border-collapse: collapse;" border="1.">
<tr style="border-bottom: 2pt solid black;">
<th style="border-right: 2pt solid black;"></th>
<th style="text-align: center; padding:5pt"> <math>40^\circ C </math></th>
<th style="text-align: center; padding:5pt"> <math>67^\circ C </math></th>
<th style="text-align: center; padding:5pt"> <math>100^\circ C </math></th>
</tr>
<tr>
<td style="text-align: center; padding:5pt; border-right: 2pt solid black;"> <math> y_1(0) </math> </td>
<td style="text-align: center; padding:5pt"> <math> 1.7066 </math> </td>
<td style="text-align: center; padding:5pt"> <math> 1.6749 </math> </td>
<td style="text-align: center; padding:5pt"> <math> 1.5608 </math> </td>
</tr>
<tr>
<td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_2(0)</math></td>
<td style="text-align: center; padding:5pt"><math>8.32</math></td>
<td style="text-align: center; padding:5pt"><math>8.2262</math></td>
<td style="text-align: center; padding:5pt"><math>8.3546</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_3(0)</math></td>
<td style="text-align: center; padding:5pt"><math>0.01</math></td>
<td style="text-align: center; padding:5pt"><math>0.0104</math></td>
<td style="text-align: center; padding:5pt"><math>0.0082</math></td>
</tr>
<tr>
<td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_4(0)</math></td>
<td style="text-align: center; padding:5pt"><math>0</math></td>
<td style="text-align: center; padding:5pt"><math>0.0017</math></td>
<td style="text-align: center; padding:5pt"><math>0.0086</math></td>
</tr>
</table>

The initial model conditions in addition to those given in the
data sets are:

<math>
\begin{array}{l}
y_5 = 0 \\
y_6 = 0.0131 \\
y_7 = \frac{1}{2} \cdot \left( -K_2 + \sqrt{K_2^2 + 4K_2 y_1(0)} \right) \\
y_8 = y_7 \\
y_9 = 0 \\
y_{10} = 0
\end{array}
</math>


== Optimal Experimental Design Problem ==
To be specified.


We are interested in when to measure (with an upper bound <math>M_i</math> on the measuring time for each observable).
We define

<math>
\begin{align}
f_y(\cdot) &\in \mathbb{R}^{10 \times 10} \quad &&\text{ with } (f_y)_{i,j} = f_{y,i,j}, \\
f_\theta(\cdot) &\in \mathbb{R}^{10 \times 9} \quad &&\text{ with } (f_\theta)_{i,j} = f_{\theta,i,j}
\end{align}
</math>

In this approach, we add the so-called sensitivities <math>G=dy/d\theta</math>. For the differential equations this means

<math>
\dot{G}(t) = f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta), \quad G(0) = \frac{\partial y(0)}{\partial \theta}
</math>


Now we formulate the OED problem as described in [[#OEDUDE | [2]]].

<math>
\begin{array}{lll}
\displaystyle \min_{y,G,F,z,w} && \text{trace} \; \left( F^{-1}(t_f) \right) \\
\text{subject to} \\
\quad \dot{y}(t) & = & f(y(t),\theta) \\
\quad \dot{G}(t) & = & f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta) \\
\quad \dot{F}(t) & = & \sum_{i=1}^{n_o} w_i(t)(h^i_y(y(t))G(t))^T(h^i_y(y(t))G(t)) \\
\quad \dot{z}(t) & = & w(t), \\
\quad y(0) & = & y_0 \\
\quad G(0) & = & \frac{\partial y(0)}{\partial \theta} \\
\quad F(0) & = & 0, \\
\quad z(0) & = & 0 \\
\quad w(t) & \in & \mathcal{W} \\
\quad z_i(t_f) & \leq & M_i
\end{array}
</math>


Here <math>h:\mathbb{R}^{10} \rightarrow \mathbb{R}^{n_o}</math> is the observed function. The
evolution of the symmetric matrix <math>F: \left[0,t_f \right] \rightarrow \mathbb{R}^{9 \times 9}</math> is given by the weighted sum of observability Gramians
<math>h^i_y (y(t)) G(t), \ i = 1,\ldots,n_o</math> for each observed function of states. The weights <math>w_i (t) \in \{0, 1\}, \ i = 1,\ldots,n_o</math>
are the (binary) sampling decisions, where <math> w_i (t) = 1 </math> denotes the decision to perform a measurement at
time <math>t</math>.

== Miscellaneous and Further Reading ==
To be specified.

== References ==
[1] "Nonlinear Parameter Estimation: a Case Study Comparison" by L. T. Biegler and J. J. Damiano 
[2] "Optimal Experimental Design for Universal Differential Equations" by C. Plate, C.J. Martensen and S. Sager 
[3] "Parameter estimation in nonlinear systems" by W.J.H. Stortelder 


[[Category:MIOCP]]
[[Category:Optimum Experimental Design]]
[[Category:ODE model]]

DOW Experimental Design

2024-09-24T10:36:43Z

SebastianSager: Some readjustments

{{Dimensions
|nd = 1
|nx = 11
|nw = 2
|nc = 4
|nre = 11
}}

The '''DOW Experimental Design problem''' models the OED problem for the parameter estimation problem formulated by the DOW Chemical Co. in 1981. The following formulation is taken from "Nonlinear Parameter Estimation:
a Case Study Comparison" by L. T. Biegler and J. J. Damiano add quote.

== Chemical background ==

The chemical species are disguised for proprietary reasons and the desired reaction is given by <math>HA + 2BM \rightarrow AB + HBMH</math>, where <math>AB</math> is the desired product. The reactions are described as follows:

Slow Kinetic Reactions:

<math>
\begin{array}{c}
M^- + BM \underset{k_{-1}}{\leftarrow} \overset{k_1}{\rightarrow} MBM^- \\
A^- + BM \overset{k_2}{\rightarrow} ABM^- \\
M^- + AB \underset{k_{-3}}{\leftarrow} \overset{k_3}{\rightarrow} ABM^-
\end{array}
</math>


Acid-Base Reactions:

<math>
\begin{array}{c}
MBMH \leftarrow K_1 \rightarrow MBM^- + H^+ \\
HA \leftarrow K_2 \rightarrow A^- + H^+ \\
HABM \leftarrow K_2 \rightarrow ABM^- + H^+
\end{array}
</math>

In order to devise a model to account for these reactions, it is
first necessary to distinguish between the overall concentration of
a species and the concentration of its neutral form. Overall concentrations are defined for three components based on neutral
and ionic species

<math>
\begin{array}{c}
\left[HBMH\right] = \left[ (MBMH)_N \right] + \left[ MBM^- \right] \\
\left[HA\right] = \left[ (HA)_N \right] + \left[ A^- \right] \\
\left[HABM\right] = \left[ (HABM)_N \right] + \left[ ABM^- \right]
\end{array}
</math>

Here <math>[ \ ]</math> denotes the concentration of the species in <math>\operatorname{gmol}/\operatorname{kg}</math>.
By assuming the rapid acid-base reactions are at equilibrium, the equilibrium constants <math>K_1,K_2,K_3</math> can be defined as

<math>
\begin{array}{c}
K_1 = \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\
K_1 = \frac{[A^-][H^+]}{[(HA)_N]} \\
K_1 = \frac{[ABM^-][H^+]}{[(HABM)_N]}
\end{array}
</math>


The anionic species may then be represented by

<math>
\begin{array}{cr}
\left[MBM^-\right] = \frac{K_1[MBMH]}{K_1 + [H^+]} & \quad (a) \\
\left[A^-\right] = \frac{K_2[HA]}{K_2 + [H^+]} & \quad (b)\\
\left[ABM^-\right] = \frac{K_3[HABM]}{K_3 + [H^+]}& \quad (c)
\end{array}
</math>

Material balance equations for the three reactants in the slow
kinetic reactions yield:

<math>
\begin{array}{cr}
\frac{d[M^-]}{dt} = -k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] - k_3 \left[ M^- \right]\left[ AB \right] + k_{-1} \left[ ABM^- \right] & \quad (d)\\
\frac{d[BM]}{dt} = -k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] - k_2 \left[ A^- \right]\left[ BM \right] & \quad (e)\\
\frac{d[AB]}{dt} = -k_3 \left[ M^- \right] \left[ AB \right] + k_{-3} \left[ ABM^- \right] & \quad (f)
\end{array}
</math>

From stoichiometry, rate expressions can also be written for
the total species:

<math>
\begin{array}{cr}
\frac{d[MBMH]}{dt} = k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right] & \quad (g)\\
\frac{d[HA]}{dt} = k_2 \left[ A^- \right] \left[ BM \right] & \quad (h)\\
\frac{d[HABM]}{dt} = k_2 \left[ A^- \right] \left[ BM \right] + k_3 \left[ M^- \right] \left[ AB \right] - k_{-3} \left[ ABM^- \right]& \quad (i)
\end{array}
</math>

An electroneutrality constraint gives the hydrogen ion con-
centration <math>\left[ H^+ \right]</math> as

<math>
\left[ H^+ \right] + \left[ Q^+ \right] = \left[ M^- \right] + \left[ MBM^- \right] + \left[ A^- \right] + \left[ ABM^- \right] \quad \quad (j)
</math>


Based on similarities of reacting species, we assume

<math> k_3 = k_1, \quad k_{-3} = \frac{1}{2} k_{-1} </math>

With these assumptions, the three rate constants <math>k_1,k_2</math> and <math>k_3</math> must be estimated. Each of these can be fitted with two adjustable model parameters, assuming an Arrhenius temperature dependence. That is

<math>
\begin{array}{c}
k_i = \alpha_i \exp(-E_i /(RT)), \quad i \in \{-1,1,2\}
\end{array}
</math>

Here <math>R</math> is the gas constant and <math>T</math> is reaction temperature in Kelvins. The parameter <math>\alpha_i</math>, given in <math>\operatorname{mol}/( \operatorname{kg} \cdot \operatorname{h})</math>, represent the pre-exponential factor and <math>E_i</math>, with unit <math>\operatorname{cal}/{\operatorname{mol}}</math>, is the activation energy.

== Mathematical formulation ==

The chemical processes <math>(a)-(j)</math> can be expressed mathematically as six differential equations and four algebraic equations:

<math>
\begin{array}{lr}
\dot{y}_1 = -k_2 y_8 y_2 & \quad (1),(h) \\
\dot{y}_2 = -k_1 y_6 y_2 + k_{-1} y_{10} - k_2 y_8 y_2 & \quad (2),(e) \\
\dot{y}_3 = -k_2 y_8 y_2 + k_1 y_6 y_4 - \frac{1}{2} k_{-1} y_9 & \quad (3),(i) \\
\dot{y}_4 = -k_1 y_6 y_4 + \frac{1}{2} k_{-1} y_9 & \quad (4),(f) \\
\dot{y}_5 = -k_1 y_6 y_2 + k_{-1} y_{10} & \quad (5),(g) \\
\dot{y}_6 = -k_1 (y_6 y_2 + y_6 y_4) + k_{-1} (y_{10} + \frac{1}{2} y_9) & \quad (6),(d) \\
y_7 = -\left[ Q^+ \right] + y_6 + y_8 + y_9 + y_{10} & \quad (7),(j)\\
y_8 = \frac{\theta_8 y_1}{\theta_8 + y_7} & \quad (8),(b)\\
y_9 = \frac{\theta_9 y_3}{\theta_9 + y_7} & \quad (9),(c)\\
y_{10} = \frac{\theta_7 y_5}{\theta_7 + y_7} & \quad (10),(a)
\end{array}
</math>

Here the letters in parentheses stand for the corresponding chemical process and the quantity <math>\left[Q^+\right]</math> is a constant during the reaction.
The nine parameters form the vector

<math>
\theta = (\alpha_1,E_1,\alpha_2,E_2,\alpha_{-1},E_{-1},K_1,K_2,K_3)
</math>

The predicted concentrations form the vector

<math>
y = (HA,BM,HABM,AB,MBMH,M^-,H^+,A^-,ABM^-,MBM^-)
</math>


Let <math>f_k(\cdot)</math> denote the right hand side of equation <math>(k)</math> for <math>k=1,\ldots,6</math>.
We reformulate the last four algebraic equations as differential ones:

<math>
\begin{array}{l r}
\dot{y}_7 = f_7 = f_6 + f_8 + f_9 + f_{10} & \quad (7') \\
\dot{y}_8 = f_8 = \frac{\theta_8 f_1 \cdot (\theta_8 + f_7) - \theta_8 y_1 f_7}{(\theta_8 + y_7)^2} & \quad (8') \\
\dot{y}_9 = f_9 = \frac{\theta_9 f_3 \cdot (\theta_9 + f_7) - \theta_9 y_3 f_7}{(\theta_9 + y_7)^2} & \quad (9')\\
\dot{y}_{10} = f_{10} = \frac{\theta_7 f_5 \cdot (\theta_7 + f_7) - \theta_7 y_5 f_7}{(\theta_7 + y_7)^2} & \quad (10')
\end{array}
</math>

The right hand sides of <math>(1)-(6)</math> and <math>(7')-(9')</math> are summarized as the vector-valued function <math>f</math>.
Moreover, let

<math>
\begin{array}{l}
f_{y,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,10\} \\
f_{\theta,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial \theta_n}, \quad m \in \{1,\ldots,10\}; \ n\in \{1,\ldots,9\}
\end{array}
</math>


== Parameters ==

The initial parameter estimates are:

<math>
\begin{array}{c}
\alpha_1 = 2.0 \times 10^{13} \\
E_1 = 2.0 \times 10^{4} \\
\alpha_2 = 2.0 \times 10^{13} \\
E_2 = 2.0 \times 10^{4} \\
\alpha_{-1} = 4.3 \times 10^{15} \\
E_{-1} = 2.0 \times 10^{4} \\
K_1 = 1.0 \times 10^{-17} \\
K_2 = 1.0 \times 10^{-11} \\
K_3 = 1.0 \times 10^{-17} \\
\end{array}
</math>

There are four datasets for different temperatures <math>T</math>, with corresponding starting values

<math>
\begin{array}{lccc}
& 40^\circ C & 67^\circ C & 100^\circ C \\
y_1(0) = & 1.7066 & 1.6749 & 1.5608 \\
y_2(0) = & 8.32 & 8.2262 & 8.3546 \\
y_3(0) = & 0.01 & 0.0104 & 0.0082 \\
y_4 (0) = & 0 & 0.0017 & 0.0086
\end{array}
</math>

The initial model conditions in addition to those given in the
data sets are:

<math>
\begin{array}{l}
y_5 = 0 \\
y_6 = 0.0131 \\
y_7 = \frac{1}{2} \cdot \left( -K_2 + \sqrt{K_2^2 + 4K_2 y_1(0)} \right) \\
y_8 = y_7 \\
y_9 = 0 \\
y_{10} = 0
\end{array}
</math>


== Optimal Experimental Design Problem ==
To be specified.


We are interested in when to measure (with an upper bound <math>M_i</math> on the measuring time for each observable).
We define

<math>
\begin{array}{l}
f_y(\cdot) \in \mathbb{R}^{10 \times 10} \quad \text{ with } (f_y)_{i,j} = f_{y,i,j}, \\
f_\theta(\cdot) \in \mathbb{R}^{10 \times 9} \quad \text{ with } (f_\theta)_{i,j} = \frac{\partial f_i}{\partial \theta_j}
\end{array}
</math>

In this approach, we add the so-called sensitivities <math>G=dy/d\theta</math>. For the differential equations this means

<math>
\dot{G}(t) = f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta), \quad G(0) = \frac{\partial y(0)}{\partial \theta}
</math>


Now we formulate the OED problem as described in (Optimal Experimental Design for Universal Differential Equations add quote)

<math>
\begin{array}{lll}
\displaystyle \min_{y,G,F,z,w} && \text{trace} \; \left( F^{-1}(t_f) \right) \\
\text{subject to} \\
\quad \dot{y}(t) & = & f(y(t),\theta) \\
\quad \dot{G}(t) & = & f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta) \\
\quad \dot{F}(t) & = & \sum_{i=1}^{n_o} w_i(t)(h^i_y(y(t))G(t))^T(h^i_y(y(t))G(t)) \\
\quad \dot{z}(t) & = & w(t), \\
\quad y(0) & = & y_0 \\
\quad G(0) & = & \frac{\partial y(0)}{\partial \theta} \\
\quad F(0) & = & 0, \\
\quad z(0) & = & 0 \\
\quad w(t) & \in & \mathcal{W} \\
\quad z_i(t_f) & \leq & M_i
\end{array}
</math>


Here <math>h:\mathbb{R}^{10} \rightarrow \mathbb{R}^{n_o}</math> is the observed function. The
evolution of the symmetric matrix <math>F: \left[0,t_f \right] \rightarrow \mathbb{R}^{9 \times 9}</math> is given by the weighted sum of observability Gramians
<math>h^i_y (y(t)) G(t), \ i = 1,\ldots,n_o</math> for each observed function of states. The weights <math>w_i (t) \in \{0, 1\}, \ i = 1,\ldots,n_o</math>
are the (binary) sampling decisions, where <math> w_i (t) = 1 </math> denotes the decision to perform a measurement at
time <math>t</math>.

== Miscellaneous and Further Reading ==
To be specified.

== References ==
To be added.



[[Category:MIOCP]]
[[Category:Optimum Experimental Design]]
[[Category:ODE model]]

User:SebastianSager

2024-09-24T09:10:06Z

SebastianSager: /* Personal */

__NOTOC__


[[Image:Sager.jpg|thumb]]

== Personal ==
{|
|-
| '''Name''': || Prof. Dr. Sebastian Sager
|-
| '''Affiliation 1''': || [https://mathopt.de Mathematical Algorithmic Optimization], [https://ovgu.de University of Magdeburg]
|-
| '''Affiliation 2''': || [https://www.mpi-magdeburg.mpg.de/oml Optimization and Machine Learning], [https://www.mpi-magdeburg.mpg.de/ Max-Planck Institute Magdeburg]
|-
| '''Background''': || Mathematics
|}

== Contact ==
{|
|-
| '''Address''': || Universitätsplatz 2, 39108 Magdeburg, Germany
|-
| '''Web''': || work group [https://mathopt.de http://mathopt.de]
|-
| '''Email''': || sager -- ovgu.de
|}

== Miscellaneous ==
I started these pages and am the owner of the domain. Hence, if you have any questions or comments concerning mintOC.de, feel free to write me an email.


[[Category:Community]]

User:SebastianSager

2024-09-24T09:07:23Z

SebastianSager: /* Contact */

__NOTOC__


[[Image:Sager.jpg|thumb]]

== Personal ==
{|
|-
| '''Name''': || Prof. Dr. Sebastian Sager
|-
| '''Affiliation''': || [http://mathopt.de Mathematical Algorithmic Optimization], [http://ovgu University of Magdeburg]
|-
| '''Background''': || Mathematics
|}

== Contact ==
{|
|-
| '''Address''': || Universitätsplatz 2, 39108 Magdeburg, Germany
|-
| '''Web''': || work group [https://mathopt.de http://mathopt.de]
|-
| '''Email''': || sager -- ovgu.de
|}

== Miscellaneous ==
I started these pages and am the owner of the domain. Hence, if you have any questions or comments concerning mintOC.de, feel free to write me an email.


[[Category:Community]]

Main Page

2024-09-24T09:06:27Z

SebastianSager: /* Direct links to Problems (add) */

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Systems biology | Systems biology]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category:GIOC | GIOC]] - [[:Category:Optimum_Experimental_Design | OED]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Control of Heat Equation with Actuator Placement | Actuator Placement]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[DOW Experimental Design | DOW OED]] - [[Electric Car]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]] - [[Industrial robot]] - [[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

Template:Current News

2024-09-24T09:04:20Z

SebastianSager:

{| style="background:#EEEEFF;"

{{News|2024/09/24| Added [[DOW_Experimental_Design | DOW]] optimal experimental design problem}}
{{News|2023/10/19| Added new category [[:Category:GIOC | generalized inverse optimal control]]}}
{{News|2023/10/19| Added [[LinearMetabolic]] problem}}

|}

<noinclude>
== Older news ==

{| style="background:#eaecd0;"

{{News|2019/3/14| Added multiple control problems with [[:Category: Gekko | Python GEKKO]] implementation}}
{{News|2018/9/12| Added [[Control of Transmission Lines]] problem}}
{{News|2016/7/26 | Added [[Industrial robot]] problem}}
{{News|2016/7/24 | Added [[Continuously Stirred Tank Reactor problem | CSTR problem]]}}
{{News|2016/6/30 | Added [[Electric Car]] problem}}
{{News|2016/5/5| Added multiple control problems with [[:Category: AMPL/TACO | AMPL with TACO]] implementation}}
{{News|2016/2/23|Added [[Control of Heat Equation with Actuator Placement | Actuator Placement]] control problem and [[:Category:Elliptic | Elliptic]], [[:Category:Parabolic | Parabolic]] and [[:Category:Hyperbolic | Hyperbolic]] categories}}
{{News|2016/1/10|Added [[Van der Pol Oscillator]], [[Batch reactor]], [[Cushioned Oscillation]] and [[Goddart's rocket problem]] control problems}}
{{News|2015/11/10|Moved mintoc.de to a new server}}
{{News|2012/09/01|Added the [[Lotka Experimental Design]] and a [[:Category:Optimum Experimental Design | category for experimental design problems]]}}
{{News|2011/09/29|Added the [[:Category:AMPL/TACO | first set of AMPL optimal control problems]] using the TACO toolkit}}
{{News|2010/11/21|Added New York [[Subway ride]] control problem}}
{{News|2010/11/18|Extended description of [[:Category:Problem characterization | problem characterization]]}}
{{News|2010/11/18|Description of benchmark library as [http://mathopt.de/PUBLICATIONS/Sager2011b.pdf pdf file] preprint}}
{{News|2010/08/16|Added [[Bang-bang approximation of a traveling wave]] 1D PDE example}}
{{News|2010/02/11|Launch of EU project [http://embocon.org embocon.org]}}
{{News|2009/11/20|Added [[External Links]] page}}
{{News|2009/10/26|Revision and Correction of several [[:Category:Optimica | optimica]] models}}
{{News|2009/08/26|Knitro solution for [[F-8 aircraft (AMPL)]]}}
{{News|2009/08/11|[[F-8 aircraft]] new local minimum found with [[:Category:Optimica | optimica]]/ipopt}}
{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}

|}

== Adding items ==
If you want to add events you have to edit this page. At the top of the page you will see something like this:
<source lang="text">
{| style="background:#EEEEFF;"

{{News|2008/08/14, 15:15|Presentation of the restructured group wiki in the group meeting, room 432.}}

|}
</source>

To add a new item you have to enter a line of the form:
<source lang="text">
{{News|date|content}}
</source>
where <code>date</code> is the date and <code>content</code> is the subject and further information. The field <code>News</code> is mandatory. Please write the <code>date</code> in the format '''YYYY/MM/DD, HH:MM''' to keep it nice and clean. The field <code>content</code> can contain any information you like (but please keep it short) and even wiki syntax is possible.

''Please keep the items in a chronological order and remove older ones to keep the list short.''

[[Category:News]]

</noinclude>

Main Page

2024-09-24T09:02:23Z

SebastianSager: /* Problem characterization (add) */

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Systems biology | Systems biology]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category:GIOC | GIOC]] - [[:Category:Optimum_Experimental_Design | OED]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Control of Heat Equation with Actuator Placement | Actuator Placement]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[Electric Car]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]] - [[Industrial robot]] - [[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

External Links

2024-06-11T12:06:48Z

SebastianSager:

= Work Groups =

An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.

* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson
* [https://www.tu-braunschweig.de/mo/team/kirches Mathematical Optimization], led by Christian Kirches
* [http://www.mathematik.tu-clausthal.de/arbeitsgruppen/kontinuierliche-optimierung/ Continuous Optimization], led by Andreas Potschka
* [https://www.mathematik.hu-berlin.de/de/forschung/forschungsgebiete/mathematische-optimierung/Falk-hante/falk-hante-fp Applied Mathematics], led by Falk Hante

Feel free to update this list!

= Benchmark Libraries =

=== Embedded Optimization and Control ===

* European FP 7 project [http://embocon.org EMBOCON]

=== Hybrid Systems ===

* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.

=== Optimal Control ===

* The [http://tomdyn.com/ PROPT homepage] (a matlab toolkit for dynamic optimization using collocation) states over 100 test cases from different applications with their results and computation time.
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems.
* The ESA provides a [https://www.esa.int/gsp/ACT/projects/gtop/ test set of global optimization spacecraft trajectory problems] and their best putative solutions.
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization
* The [http://plato.asu.edu/pdecon.html PDECON collection] contains discretized PDE problems in AMPL as part of a [http://plato.asu.edu/bench.html larger benchmark effort] by Hans Mittelmann

=== MINLP ===

* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]
* MINLPlib <bib id="Bussieck2003a" />
* [http://www.gamsworld.org/performance GAMSWORLD]

=== Generic Optimization ===

* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)

== References ==
<biblist />

LinearMetabolic

2023-10-20T11:35:17Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 4
|nu = 3
|nc = 1
|nri = 0
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

The gIOC is given by

<math>
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|
</math>
subject to
<math>
(x^*,u^*, T^*) \in \displaystyle \arg \min_{x, u, T} w_1 \cdot \int_{0}^{T} (x_2 + x_3) \, dt + w_2 \cdot \int_{0}^{T} 1 \,dt
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x_1}(t) &=& \displaystyle 0 \\
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\
x(0) &=& (1, 0, 0, 0) \\
x_4(T) &=& 0.9 \\
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\
w_1, w_2 \in [0,1] \\
w_1 + w_2 = 1
\end{array}
</math>

Here the four differential states and three control functions stand for the metabolite concentrations (<math> x_1, x_2, x_3, x_4 </math>) and the enzyme concentrations (<math> u_1, u_2, u_3 </math>). The two objective function candidates are the sum of intermediate metabolite concentrations and the transition time. The three model parameters <math>p_1, p_2, p_3</math> are kinetic parameters in mass action expressions.

== Data ==

The following table gives different values for <math>\eta</math>, componentwise for all four entries of <math>x</math> as time series, for different combinations of weights and noise levels that were used to produce the synthetic data via solution of forward optimal control problems.

{| class="wikitable"
|+Data
|-
|weights <math>w_1, w_2</math>
|parameters <math>p_1, p_2, p_3</math>
|noise
|data <math>\eta</math>
|-
|[0.232, 0.768]
|[1,1,1]
| 0
|
[[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
[0, 0.215, 0.4298, 0.6446, 0.8593, 1.074, 1.2885, 1.503, 1.6591, 1.3381, 1.0792, 0.8704, 0.7021, 0.5662, 0.4849, 0.4846, 0.4844, 0.4842, 0.484, 0.4838, 0.4836],
[0, 0.0, 0.0002, 0.0004, 0.0007, 0.001, 0.0015, 0.002, 0.0391, 0.36, 0.6186, 0.827, 0.9948, 1.1301, 1.1405, 0.92, 0.7423, 0.5989, 0.4832, 0.3899, 0.3147],
[0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0008, 0.0013, 0.0019, 0.0729, 0.2935, 0.4715, 0.6152, 0.731, 0.8245, 0.9]]
|-
|[0.232, 0.768]
|[1,1,1]
|0.2
|
[[1., 0.9590812,0.82266497,0.95877308,0.60614596,1.20781224, 0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347, 0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004, 0.83411596,0.73187308,0.70399766],
[0.,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728, 0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261, 0.07589526,0.64051491,0.47184952],
[ 0., 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845, 1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544, 0.30229594, 0.43876605, 0.22073986],
[ 0.,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 , -0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781, -0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333, 0.8178261 , 0.53498205, 1.1021996]]
|-
|[0,1]
|[1,1,1]
|0
|
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0.211317, 0.422635, 0.633952, 0.845269, 1.05659, 1.2679, 1.47922, 1.69054, 1.90186, 1.95601, 1.58342, 1.2818, 1.03764, 0.839981, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458],
[0, 3.40252e-10, 2.02262e-09, 6.04129e-09, 1.32631e-08, 2.47434e-08, 4.23248e-08, 6.98021e-08, 1.16723e-07, 2.20847e-07, 0.103509, 0.476095, 0.777716, 1.02188, 1.21954, 1.36367, 1.10391, 0.893635, 0.72341, 0.58561, 0.474061],
[0, 2.79389e-19, 2.05361e-18, 7.26553e-18, 1.82666e-17, 3.77412e-17, 6.87107e-17, 1.14894e-16, 1.82155e-16, 2.86358e-16, 3.73778e-11, 7.31532e-10, 5.75723e-09, 2.28945e-08, 7.3467e-08, 0.0103927, 0.270148, 0.480425, 0.65065, 0.78845, 0.9]]
|-
|[0,1]
|[1,1,1]
|0.1
|
[[1.,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996, 1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442, 1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444, 0.93306189,1.04808763,0.87659083],
[0.,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772, 1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121, 1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009, 0.73549769,0.62452413,0.61907767],
[0., 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101, 0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228, 0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058, 0.61562944, 0.6339817 , 0.36613472],
[0.,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716, -0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985, 0.65331917, 0.78324169, 0.81304223]]
|-
|}

== Examples ==
The plots show noise free data for different weights from forward optimal control problems.
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3">
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].
</gallery>

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

LinearMetabolic

2023-10-20T11:28:15Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 4
|nw = 3
|nc = 2
|nri = 2
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

The gIOC is given by

<math>
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|
</math>
subject to
<math>
(x^*,u^*, T^*) \in \displaystyle \arg \min_{x, u, T} w_1 \cdot \int_{0}^{T} (x_2 + x_3) \, dt + w_2 \cdot \int_{0}^{T} 1 \,dt
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x_1}(t) &=& \displaystyle 0 \\
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\
x(0) &=& (1, 0, 0, 0) \\
x_4(T) &=& 0.9 \\
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\
w_1, w_2 \in [0,1] \\
w_1 + w_2 = 1
\end{array}
</math>

Here the four differential states and three control functions stand for the metabolite concentrations (<math> x_1, x_2, x_3, x_4 </math>) and the enzyme concentrations (<math> u_1, u_2, u_3 </math>). The objective function candidates are the sum of intermediate metabolite concentrations and the transition time. The parameters (<math> k_1, k_2, k_3 </math>) are kinetic parameters in mass action expressions.

== Data ==

The following table gives different values for <math>\eta</math>, componentwise for all four entries of <math>x</math> as time series, for different combinations of weights and noise levels that were used to produce the synthetic data via solution of forward optimal control problems.

{| class="wikitable"
|+Data
|-
|weights <math>w_1, w_2</math>
|parameters <math>p_1, p_2, p_3</math>
|noise
|data <math>\eta</math>
|-
|[0.232, 0.768]
|[1,1,1]
| 0
|
[[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
[0, 0.215, 0.4298, 0.6446, 0.8593, 1.074, 1.2885, 1.503, 1.6591, 1.3381, 1.0792, 0.8704, 0.7021, 0.5662, 0.4849, 0.4846, 0.4844, 0.4842, 0.484, 0.4838, 0.4836],
[0, 0.0, 0.0002, 0.0004, 0.0007, 0.001, 0.0015, 0.002, 0.0391, 0.36, 0.6186, 0.827, 0.9948, 1.1301, 1.1405, 0.92, 0.7423, 0.5989, 0.4832, 0.3899, 0.3147],
[0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0008, 0.0013, 0.0019, 0.0729, 0.2935, 0.4715, 0.6152, 0.731, 0.8245, 0.9]]
|-
|[0.232, 0.768]
|[1,1,1]
|0.2
|
[[1., 0.9590812,0.82266497,0.95877308,0.60614596,1.20781224, 0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347, 0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004, 0.83411596,0.73187308,0.70399766],
[0.,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728, 0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261, 0.07589526,0.64051491,0.47184952],
[ 0., 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845, 1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544, 0.30229594, 0.43876605, 0.22073986],
[ 0.,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 , -0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781, -0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333, 0.8178261 , 0.53498205, 1.1021996]]
|-
|[0,1]
|[1,1,1]
|0
|
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0.211317, 0.422635, 0.633952, 0.845269, 1.05659, 1.2679, 1.47922, 1.69054, 1.90186, 1.95601, 1.58342, 1.2818, 1.03764, 0.839981, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458],
[0, 3.40252e-10, 2.02262e-09, 6.04129e-09, 1.32631e-08, 2.47434e-08, 4.23248e-08, 6.98021e-08, 1.16723e-07, 2.20847e-07, 0.103509, 0.476095, 0.777716, 1.02188, 1.21954, 1.36367, 1.10391, 0.893635, 0.72341, 0.58561, 0.474061],
[0, 2.79389e-19, 2.05361e-18, 7.26553e-18, 1.82666e-17, 3.77412e-17, 6.87107e-17, 1.14894e-16, 1.82155e-16, 2.86358e-16, 3.73778e-11, 7.31532e-10, 5.75723e-09, 2.28945e-08, 7.3467e-08, 0.0103927, 0.270148, 0.480425, 0.65065, 0.78845, 0.9]]
|-
|[0,1]
|[1,1,1]
|0.1
|
[[1.,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996, 1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442, 1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444, 0.93306189,1.04808763,0.87659083],
[0.,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772, 1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121, 1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009, 0.73549769,0.62452413,0.61907767],
[0., 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101, 0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228, 0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058, 0.61562944, 0.6339817 , 0.36613472],
[0.,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716, -0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985, 0.65331917, 0.78324169, 0.81304223]]
|-
|}

== Examples ==
The plots show noise free data for different weights from forward optimal control problems.
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3">
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].
</gallery>

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

LinearMetabolic

2023-10-20T11:27:04Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 4
|nw = 3
|nc = 2
|nri = 2
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

The gIOC is given by

<math>
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|
</math>
subject to
<math>
(x^*,u^*, T^*) \in \displaystyle \arg \min_{x, u, T} w_1 \cdot \int_{0}^{T} (x_2 + x_3) \, dt + w_2 \cdot \int_{0}^{T} 1 \,dt
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x_1}(t) &=& \displaystyle 0 \\
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\
x(0) &=& (1, 0, 0, 0) \\
x_4(T) &=& 0.9 \\
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\
w_1, w_2 \in [0,1] \\
w_1 + w_2 = 1
\end{array}
</math>

Here the four differential states and three control functions stand for the metabolite concentrations (<math> x_1, x_2, x_3, x_4 </math>) and the enzyme concentrations (<math> u_1, u_2, u_3 </math>). The objective function candidates are the sum of intermediate metabolite concentrations and the transition time. The parameters (<math> k_1, k_2, k_3 </math>) are kinetic parameters in mass action expressions.

== Data ==

The following table gives different values for <math>\eta</math>, componentwise for all four entries of <math>x</math> as time series, for different combinations of weights and noise levels that were used to produce the synthetic data via solution of forward optimal control problems.

{| class="wikitable"
|+Data
|-
|weights <math>w_1, w_2</math>
|parameters <math>p_1, p_2, p_3</math>
|noise
|data <math>\eta</math>
|-
|[0.232, 0.768]
|[1,1,1]
| 0
|
[[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
[0, 0.215, 0.4298, 0.6446, 0.8593, 1.074, 1.2885, 1.503, 1.6591, 1.3381, 1.0792, 0.8704, 0.7021, 0.5662, 0.4849, 0.4846, 0.4844, 0.4842, 0.484, 0.4838, 0.4836],
[0, 0.0, 0.0002, 0.0004, 0.0007, 0.001, 0.0015, 0.002, 0.0391, 0.36, 0.6186, 0.827, 0.9948, 1.1301, 1.1405, 0.92, 0.7423, 0.5989, 0.4832, 0.3899, 0.3147],
[0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0008, 0.0013, 0.0019, 0.0729, 0.2935, 0.4715, 0.6152, 0.731, 0.8245, 0.9]]
|-
|[0.232, 0.768]
|[1,1,1]
|0.2
|
[[1., 0.9590812,0.82266497,0.95877308,0.60614596,1.20781224, 0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347, 0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004, 0.83411596,0.73187308,0.70399766],
[0.,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728, 0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261, 0.07589526,0.64051491,0.47184952],
[ 0., 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845, 1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544, 0.30229594, 0.43876605, 0.22073986],
[ 0.,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 , -0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781, -0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333, 0.8178261 , 0.53498205, 1.1021996]]
|-
|[0,1]
|[1,1,1]
|0
|
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0.211317, 0.422635, 0.633952, 0.845269, 1.05659, 1.2679, 1.47922, 1.69054, 1.90186, 1.95601, 1.58342, 1.2818, 1.03764, 0.839981, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458],
[0, 3.40252e-10, 2.02262e-09, 6.04129e-09, 1.32631e-08, 2.47434e-08, 4.23248e-08, 6.98021e-08, 1.16723e-07, 2.20847e-07, 0.103509, 0.476095, 0.777716, 1.02188, 1.21954, 1.36367, 1.10391, 0.893635, 0.72341, 0.58561, 0.474061],
[0, 2.79389e-19, 2.05361e-18, 7.26553e-18, 1.82666e-17, 3.77412e-17, 6.87107e-17, 1.14894e-16, 1.82155e-16, 2.86358e-16, 3.73778e-11, 7.31532e-10, 5.75723e-09, 2.28945e-08, 7.3467e-08, 0.0103927, 0.270148, 0.480425, 0.65065, 0.78845, 0.9]]
|-
|[0,1]
|[1,1,1]
|0.1
|
[[1.,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996, 1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442, 1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444, 0.93306189,1.04808763,0.87659083],
[0.,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772, 1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121, 1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009, 0.73549769,0.62452413,0.61907767],
[0., 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101, 0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228, 0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058, 0.61562944, 0.6339817 , 0.36613472],
[0.,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716, -0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985, 0.65331917, 0.78324169, 0.81304223]]
|-
|}

== Examples ==
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3">
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].
</gallery>

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

LinearMetabolic

2023-10-20T11:26:21Z

SebastianSager: Edited the data table

{{Dimensions
|nd = 1
|nx = 4
|nw = 3
|nc = 2
|nri = 2
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

The gIOC is given by

<math>
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|
</math>
subject to
<math>
(x^*,u^*, T^*) \in \displaystyle \arg \min_{x, u, T} w_1 \cdot \int_{0}^{T} (x_2 + x_3) \, dt + w_2 \cdot \int_{0}^{T} 1 \,dt
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x_1}(t) &=& \displaystyle 0 \\
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\
x(0) &=& (1, 0, 0, 0) \\
x_4(T) &=& 0.9 \\
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\
1 &\ge& u_1 + u_2 + u_3 \\
w_1, w_2 \in [0,1] \\
w_1 + w_2 = 1
\end{array}
</math>

Here the four differential states and three control functions stand for the metabolite concentrations (<math> x_1, x_2, x_3, x_4 </math>) and the enzyme concentrations (<math> u_1, u_2, u_3 </math>). The objective function candidates are the sum of intermediate metabolite concentrations and the transition time. The parameters (<math> k_1, k_2, k_3 </math>) are kinetic parameters in mass action expressions.

== Data ==

The following table gives different values for <math>\eta</math>, componentwise for all four entries of <math>x</math> as time series, for different combinations of weights and noise levels that were used to produce the synthetic data via solution of forward optimal control problems.

{| class="wikitable"
|+Data
|-
|weights <math>w_1, w_2</math>
|parameters <math>p_1, p_2, p_3</math>
|noise
|data <math>\eta</math>
|-
|[0.232, 0.768]
|[1,1,1]
| 0
|
[[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
[0, 0.215, 0.4298, 0.6446, 0.8593, 1.074, 1.2885, 1.503, 1.6591, 1.3381, 1.0792, 0.8704, 0.7021, 0.5662, 0.4849, 0.4846, 0.4844, 0.4842, 0.484, 0.4838, 0.4836],
[0, 0.0, 0.0002, 0.0004, 0.0007, 0.001, 0.0015, 0.002, 0.0391, 0.36, 0.6186, 0.827, 0.9948, 1.1301, 1.1405, 0.92, 0.7423, 0.5989, 0.4832, 0.3899, 0.3147],
[0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0008, 0.0013, 0.0019, 0.0729, 0.2935, 0.4715, 0.6152, 0.731, 0.8245, 0.9]]
|-
|[0.232, 0.768]
|[1,1,1]
|0.2
|
[[1., 0.9590812,0.82266497,0.95877308,0.60614596,1.20781224, 0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347, 0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004, 0.83411596,0.73187308,0.70399766],
[0.,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728, 0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261, 0.07589526,0.64051491,0.47184952],
[ 0., 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845, 1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544, 0.30229594, 0.43876605, 0.22073986],
[ 0.,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 , -0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781, -0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333, 0.8178261 , 0.53498205, 1.1021996]]
|-
|[0,1]
|[1,1,1]
|0
|
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0.211317, 0.422635, 0.633952, 0.845269, 1.05659, 1.2679, 1.47922, 1.69054, 1.90186, 1.95601, 1.58342, 1.2818, 1.03764, 0.839981, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458],
[0, 3.40252e-10, 2.02262e-09, 6.04129e-09, 1.32631e-08, 2.47434e-08, 4.23248e-08, 6.98021e-08, 1.16723e-07, 2.20847e-07, 0.103509, 0.476095, 0.777716, 1.02188, 1.21954, 1.36367, 1.10391, 0.893635, 0.72341, 0.58561, 0.474061],
[0, 2.79389e-19, 2.05361e-18, 7.26553e-18, 1.82666e-17, 3.77412e-17, 6.87107e-17, 1.14894e-16, 1.82155e-16, 2.86358e-16, 3.73778e-11, 7.31532e-10, 5.75723e-09, 2.28945e-08, 7.3467e-08, 0.0103927, 0.270148, 0.480425, 0.65065, 0.78845, 0.9]]
|-
|[0,1]
|[1,1,1]
|0.1
|
[[1.,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996, 1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442, 1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444, 0.93306189,1.04808763,0.87659083],
[0.,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772, 1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121, 1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009, 0.73549769,0.62452413,0.61907767],
[0., 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101, 0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228, 0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058, 0.61562944, 0.6339817 , 0.36613472],
[0.,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716, -0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985, 0.65331917, 0.78324169, 0.81304223]]
|-
|}

== Examples ==
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3">
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].
</gallery>

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

LinearMetabolic

2023-10-20T06:46:23Z

SebastianSager: Removed old code and p

{{Dimensions
|nd = 1
|nx = 4
|nw = 3
|nc = 2
|nri = 2
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

Under construction...

The gIOC is given by

<math>
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|
</math>
subject to
<math>
(x^*,u^*, T^*) \in \displaystyle \arg \min_{x, u, T} w_1 \cdot \int_{0}^{T} (x_2 + x_3) \, dt + w_2 \cdot \int_{0}^{T} 1 \,dt
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x_1}(t) &=& \displaystyle 0 \\
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\
x(0) &=& (1, 0, 0, 0) \\
x_4(T) &=& 0.9 \\
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\
w_1, w_2 \in [0,1] \\
w_1 + w_2 = 1
\end{array}
</math>

Here the four differential states and three control functions stand for the metabolite concentrations (<math> x_1, x_2, x_3, x_4 </math>) and the enzyme concentrations (<math> u_1, u_2, u_3 </math>). The objective function candidates are the sum of intermediate metabolite concentrations and the transition time. The parameters (<math> k_1, k_2, k_3 </math>) are kinetic parameters in mass action expressions.

== Data ==
{| class="wikitable"
|+Data
|-
|weights
|parameters
|noise
|data
|-
|[0.232, 0.768]
|[1,1,1]
| 0
|[[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
[0, 0.215, 0.4298, 0.6446, 0.8593, 1.074, 1.2885, 1.503, 1.6591, 1.3381, 1.0792, 0.8704, 0.7021, 0.5662, 0.4849, 0.4846, 0.4844, 0.4842, 0.484, 0.4838, 0.4836],
[0, 0.0, 0.0002, 0.0004, 0.0007, 0.001, 0.0015, 0.002, 0.0391, 0.36, 0.6186, 0.827, 0.9948, 1.1301, 1.1405, 0.92, 0.7423, 0.5989, 0.4832, 0.3899, 0.3147],
[0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0008, 0.0013, 0.0019, 0.0729, 0.2935, 0.4715, 0.6152, 0.731, 0.8245, 0.9]]

|-
|[0.232, 0.768]
|[1,1,1]
|0.2
|[[1. ,0.9590812 ,0.82266497,0.95877308,0.60614596,1.20781224,
0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347,
0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004,
0.83411596,0.73187308,0.70399766],[0. ,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,
1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728,
0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261,
0.07589526,0.64051491,0.47184952],[ 0. , 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,
0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845,
1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544,
0.30229594, 0.43876605, 0.22073986],[ 0. ,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 ,
-0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781,
-0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333,
0.8178261 , 0.53498205, 1.1021996 ]]
|-
|[0,1]
|[1,1,1]
|0
|[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0.211317, 0.422635, 0.633952, 0.845269, 1.05659, 1.2679, 1.47922, 1.69054, 1.90186, 1.95601, 1.58342, 1.2818, 1.03764, 0.839981, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458, 0.685458],
[0, 3.40252e-10, 2.02262e-09, 6.04129e-09, 1.32631e-08, 2.47434e-08, 4.23248e-08, 6.98021e-08, 1.16723e-07, 2.20847e-07, 0.103509, 0.476095, 0.777716, 1.02188, 1.21954, 1.36367, 1.10391, 0.893635, 0.72341, 0.58561, 0.474061],
[0, 2.79389e-19, 2.05361e-18, 7.26553e-18, 1.82666e-17, 3.77412e-17, 6.87107e-17, 1.14894e-16, 1.82155e-16, 2.86358e-16, 3.73778e-11, 7.31532e-10, 5.75723e-09, 2.28945e-08, 7.3467e-08, 0.0103927, 0.270148, 0.480425, 0.65065, 0.78845, 0.9]]

|-
|[0,1]
|[1,1,1]
|0.1
|[[1. ,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996,
1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442,
1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444,
0.93306189,1.04808763,0.87659083],[0. ,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772,
1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121,
1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009,
0.73549769,0.62452413,0.61907767],[ 0. , 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101 ,
0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228,
0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058,
0.61562944, 0.6339817 , 0.36613472],[ 0. ,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,
-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716,
-0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985,
0.65331917, 0.78324169, 0.81304223]]

|-
|}

== Examples ==
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3">
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].
</gallery>

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

Main Page

2023-10-19T13:21:22Z

SebastianSager:

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Systems biology | Systems biology]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category:GIOC | GIOC]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Control of Heat Equation with Actuator Placement | Actuator Placement]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[Electric Car]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]] - [[Industrial robot]] - [[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

Main Page

2023-10-19T13:19:38Z

SebastianSager: added GIOC direct link

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category:GIOC | GIOC]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Control of Heat Equation with Actuator Placement | Actuator Placement]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[Electric Car]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]] - [[Industrial robot]] - [[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

Template:Current News

2023-10-19T13:18:17Z

SebastianSager:

{| style="background:#EEEEFF;"

{{News|2023/10/19| Added new category [[:Category:GIOC | generalized inverse optimal control]]}}
{{News|2023/10/19| Added [[LinearMetabolic]] problem}}

|}

<noinclude>
== Older news ==

{| style="background:#eaecd0;"

{{News|2019/3/14| Added multiple control problems with [[:Category: Gekko | Python GEKKO]] implementation}}
{{News|2018/9/12| Added [[Control of Transmission Lines]] problem}}
{{News|2016/7/26 | Added [[Industrial robot]] problem}}
{{News|2016/7/24 | Added [[Continuously Stirred Tank Reactor problem | CSTR problem]]}}
{{News|2016/6/30 | Added [[Electric Car]] problem}}
{{News|2016/5/5| Added multiple control problems with [[:Category: AMPL/TACO | AMPL with TACO]] implementation}}
{{News|2016/2/23|Added [[Control of Heat Equation with Actuator Placement | Actuator Placement]] control problem and [[:Category:Elliptic | Elliptic]], [[:Category:Parabolic | Parabolic]] and [[:Category:Hyperbolic | Hyperbolic]] categories}}
{{News|2016/1/10|Added [[Van der Pol Oscillator]], [[Batch reactor]], [[Cushioned Oscillation]] and [[Goddart's rocket problem]] control problems}}
{{News|2015/11/10|Moved mintoc.de to a new server}}
{{News|2012/09/01|Added the [[Lotka Experimental Design]] and a [[:Category:Optimum Experimental Design | category for experimental design problems]]}}
{{News|2011/09/29|Added the [[:Category:AMPL/TACO | first set of AMPL optimal control problems]] using the TACO toolkit}}
{{News|2010/11/21|Added New York [[Subway ride]] control problem}}
{{News|2010/11/18|Extended description of [[:Category:Problem characterization | problem characterization]]}}
{{News|2010/11/18|Description of benchmark library as [http://mathopt.de/PUBLICATIONS/Sager2011b.pdf pdf file] preprint}}
{{News|2010/08/16|Added [[Bang-bang approximation of a traveling wave]] 1D PDE example}}
{{News|2010/02/11|Launch of EU project [http://embocon.org embocon.org]}}
{{News|2009/11/20|Added [[External Links]] page}}
{{News|2009/10/26|Revision and Correction of several [[:Category:Optimica | optimica]] models}}
{{News|2009/08/26|Knitro solution for [[F-8 aircraft (AMPL)]]}}
{{News|2009/08/11|[[F-8 aircraft]] new local minimum found with [[:Category:Optimica | optimica]]/ipopt}}
{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}

|}

== Adding items ==
If you want to add events you have to edit this page. At the top of the page you will see something like this:
<source lang="text">
{| style="background:#EEEEFF;"

{{News|2008/08/14, 15:15|Presentation of the restructured group wiki in the group meeting, room 432.}}

|}
</source>

To add a new item you have to enter a line of the form:
<source lang="text">
{{News|date|content}}
</source>
where <code>date</code> is the date and <code>content</code> is the subject and further information. The field <code>News</code> is mandatory. Please write the <code>date</code> in the format '''YYYY/MM/DD, HH:MM''' to keep it nice and clean. The field <code>content</code> can contain any information you like (but please keep it short) and even wiki syntax is possible.

''Please keep the items in a chronological order and remove older ones to keep the list short.''

[[Category:News]]

</noinclude>

Category:GIOC

2023-10-19T12:53:11Z

SebastianSager: integer controls

This category includes all generalized inverse optimal control problems.
To formalize this problem class, we define the following bilevel problem with differential states <math>x \in \mathcal{X}</math>, controls <math>u \in \mathcal{U}</math>, model parameters <math>p \in \mathbb{R}^{n_p}</math>, and convex multipliers <math>w \in \mathcal{W}</math>

<math>
\displaystyle \min_{(p, w, x^*, u^*) \in \Omega_1} \; \| h(x^*, u^*) - \eta \| + R(p,w)
</math>
subject to
<math>
(x^*,u^*) \in \displaystyle \arg \min_{x, u} \sum_{i \in M_1} w_i \; \phi_i[x,u,p]
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x}(t) &=& \displaystyle \sum_{i \in M_2} w_i f_i(x(t), u(t), p) \\
0 &\le& \displaystyle w_i \; g_i(x(t), u(t), p) \quad \forall \; i \in M_3\\
(x,u,p,w) &\in& \Omega_2
\end{array}
</math>

as a generalized inverse optimal control problem. Here <math>\mathcal{X}</math> and <math>\mathcal{U}</math> are properly defined function spaces. The variables <math>w</math> indicate which objective functions, right hand side functions, and constraints are relevant in the inner problem. The variables are normalized for given index sets <math>M_1, M_2, M_3</math> that partition the indices from <math>1</math> to <math>n_w</math>. For normalization, we define the feasible set <math>\mathcal{W} := \{ w \in [0,1]^{n_w}: \; \textstyle \sum^{i \in M_j} w_i = 1 \text{ for } j \in \{1,2,3\}\}</math>. On the outer level, the feasible set is <math>\Omega_1 := \mathbb{R}^{n_p} \times \mathcal{W} \times \mathcal{X} \times \mathcal{U}</math>, while on the inner level <math>\Omega_2</math> contains bounds, boundary conditions, mixed path and control constraints, and more involved constraints such as dwell time constraints. We have observational data <math>\eta \in \mathbb{R}^{n_\eta}</math>, a measurement function <math>h: \mathcal{X} \times \mathcal{U} \mapsto \mathbb{R}^{n_\eta}</math>, a regularization function with a priori knowledge on parameters and weights <math>R: \mathbb{R}^{n_p} \times \mathcal{W} \mapsto \mathbb{R}</math> and candidate functionals <math>\phi_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}</math> and functions <math>f_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}^{n_x}</math> and <math>g_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}^{n_g}</math> for the unknown objective function, dynamics, and constraints, respectively. Two cases are of practical interest: first, the manual, often cumbersome and trial-and-error based a priori definition of all candidates <math>\phi_i, f_i, g_i</math> by experts and second, a systematic, but often challenging automatic symbolic regression of these unknown functions.

On the outer level, a norm <math>\| \cdot \|</math> and the regularization term <math>R</math> define a data fit (regression) problem and relate to prior knowledge and statistical assumptions. On the inner level, the above bilevel problem is constrained by a possibly nonconvex optimal control problem. The unknown parts of this inner level optimal control problem are modeled as convex combinations of a finite set of candidates (and a multiplication of constraints <math>g_i</math> with <math>w_i</math> that can be either zero or strictly positive). On the one hand the problem formulation is restrictive in the interest of a clearer presentation and might be further generalized, e.g., to multi-stage formulations involving differential-algebraic or partial differential equations. On the other hand, the problem class is quite generic and allows, e.g., the consideration of switched systems, periodic processes, different underlying function spaces <math>\mathcal{X}</math> and <math>\mathcal{U}</math>, and the usage of universal approximators such as neural networks as candidate functions.

== Extensions ==
* The functions <math>u \in \mathcal{U}</math> may also include integer controls.
* For some problems the functions may as well depend explicitely on the time <math>t</math>.
* The differential equations might depend on [[:Category:State dependent switches | state-dependent switches]].
* The variables may include [[:Category:Boolean variables | boolean variables]].
* The underlying process might be a [[:Category:Multistage process | multistage process]].
* The dynamics might be [[:Category:Unstable | unstable]].
* There might be an underlying [[:Category:Network topology | network topology]].
* The integer control functions might have been (re)formulated by means of an [[:Category:Outer convexification|outer convexification]].

Note that a Lagrange term <math>\int_{t_0}^{t_f} L( x(t), u(t), v(t), q, \rho) \; \mathrm{d} t</math> can be transformed into a Mayer-type objective functional.

[[Category:Model characterization]]

LinearMetabolic

2023-10-19T12:48:26Z

SebastianSager: Introductory text

{{Dimensions
|nd = 1
|nx = 4
|nw = 1
|nc = 2
|nri = 2
|nre = 5
}}

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem formulated and investigated in <bib id="Tsiantis2018" />. It tries to identify an a priori unknown objective function from data.

The problem is a generalization of the one studied by de Hijas-Liste et al. (2014), where it was considered as a standard optimal control problem (OCP). Here, we take the solution reference of the inner problem as the multi-objective OCP described in de Hijas-Liste et al. (2014), selecting a specific point of the resulting Pareto front. This case study is interesting, because it includes path constraints on the states and the inputs.

A 3-step linear metabolic pathway with mass action kinetics is considered. The differential states <math>x</math> are metabolite concentrations, the time-dependent control functions <math>u</math> are enzyme concentrations, and the model parameters <math>p</math> are kinetic parameters in mass action expressions. Candidate objective functionals <math>\phi</math> are the final time and the time-integral of the intermediate metabolite concentrations. The measurement function is a map to the values of the differential states and comprises measurements of metabolite concentrations. The differential equations are assumed to be fully known as standard mass action kinetics. An inequality path constraint is present (but potentially unknown in such settings) on the inner level and critical from a biological point of view: limitations due to molecular crowding impose an upper bound on the maximum total concentration of enzymes (controls) at any given time. Boundary conditions are fixed initial and terminal values of the metabolite concentrations.

== Mathematical formulation ==

Under construction...

The gIOC is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\
& \dot{x}_1 & = & \frac{K_r^2}{Mr^2} (K_m x_0 - \frac{r}{K_r} ( M g K_f + \frac{1}{2} \rho S C_x \frac{r^2}{K_r^2} x_1^2)), \\
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]
& x(t_0) &=& (0,0,0,0)^T, \\
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\
& x_0(t) & \in & [-i_{max}, i_{max}], \\
& u(t) &\in& \{-1, 1\}.
\end{array}
</math>

Here the four differential states stand for the electrical current (<math> x_0 </math>), the angular velocity (<math> x_1 </math>), the position of the car (<math> x_2 </math>), and the consumed energy (<math> x_3 </math>). The objective function <math> x_3(t_f) </math>is just a reformulation of the Lagrange-type objective function tracking the used energy over time.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Coefficient of reduction
|<math>K_r</math>
|10
|[-]
|-
|Air density
|<math>\rho</math>
|1.293
|<math> kg/m^3 </math>
|-
|Aerodynamic coefficient
|<math>C_x</math>
|0.4
|[-]
|-
|Area in the front of the vehicle
|<math>S</math>
|2
|<math> m^2 </math>
|-
|Radius of the wheel
|<math>r</math>
|0.33
|<math> m </math>
|-
|Constant representing the friction of the wheels on the road
|<math>K_f </math>
|0.03
|[-]
|-
|Coefficient of the motor torque
|<math>K_m</math>
|0.27
|[-]
|-
|Inductor resistance
|<math>R_m</math>
|0.03
|<math> Ohms </math>
|-
|inductance of the rotor
|<math>L_m</math>
|0.05
|[-]
|-
|Mass
|<math>M</math>
|250
|<math> kg </math>
|-
|Gravity constant
|<math>g</math>
|9,81
|[-]
|-
|Battery voltage
|<math>V_{alim}</math>
|150
|<math> V </math>
|-
|Resistance of the battery
|<math>R_{bat}</math>
|0.05
|<math> Ohms </math>
|-
|}

== Reference Solutions ==
We look at the particular instance of our problem with <math> t_f = 10s </math>and target set <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R} </math>, in which the car needs to cover 100m in 10s.
Figure 1 shows a plot of the differential states of the optimal trajectory of the relaxed problem (i.e. <math> u \in [-1,1] </math>instead of <math> u \in \{-1,1\} </math>) for <math> N = 1000, N </math>being the number of time discretization points. The current <math> x_0 </math>increases to its maximal value of 150A, stays there for a certain time, decreases on its minimal value of -150A, stays on this value and eventually increases slightly. This behavior corresponds to the different arcs bang, path-constrained, singular, path-constrained, bang and can be observed also in Figure 2. It shows the corresponding switching function and the optimal control. Note that the plots show data from the solution with the indirect approach.

Applying the sum up rounding strategy results in an integer-feasible chattering solution.
The resulting primal states are shown in Figure 3. One observes the high-frequency zig-zagging of the current <math> x_0 </math>that results from the fast switches in the control.

The direct and indirect approaches are local optimization techniques and only provide upper bounds for the relaxed problem and hence for the original problem. Here the indirect solution of the relaxed problem gives us a bound of <math> x_3(t_f) = 22777.2 </math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3">
Image:ElectricCar PrimalStates Plot.PNG| Primal states of an optimal trajectory for the relaxed problem on a control discretization grid with <math> N = 1,000 </math>.
Image:ElectricCar Control PathConstraints SwitchingFunction Plot.PNG| The optimal control and switching function. The dotted vertical lines show the switching times <math> \tau_i </math> where transitions between different kinds of arcs occur.
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math>
* bounded velocity, <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R}, t_f = 10s, x_2(t) \leq 45 \frac{K_r}{3.6r} \forall t </math>
* fixed final velocity, bounded velocity, longer time horizon, <math> x_2(t) \leq 30 \frac{K_r}{3.6r} \forall t, \mathcal{T} = \mathbb{R} \times \{ 30 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 15s </math>.

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

LinearMetabolic

2023-10-19T12:41:07Z

SebastianSager: Tsiantis reference

{{Dimensions
|nd = 1
|nx = 4
|nw = 1
|nc = 2
|nri = 2
|nre = 5
}}

Under construction...

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem. It tries to identify an a priori unknown objective function from data. The control unction <math> u(t) </math> is ...

The problem is discussed in detail in <bib id="Tsiantis2018" />.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\
& \dot{x}_1 & = & \frac{K_r^2}{Mr^2} (K_m x_0 - \frac{r}{K_r} ( M g K_f + \frac{1}{2} \rho S C_x \frac{r^2}{K_r^2} x_1^2)), \\
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]
& x(t_0) &=& (0,0,0,0)^T, \\
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\
& x_0(t) & \in & [-i_{max}, i_{max}], \\
& u(t) &\in& \{-1, 1\}.
\end{array}
</math>

Here the four differential states stand for the electrical current (<math> x_0 </math>), the angular velocity (<math> x_1 </math>), the position of the car (<math> x_2 </math>), and the consumed energy (<math> x_3 </math>). The objective function <math> x_3(t_f) </math>is just a reformulation of the Lagrange-type objective function tracking the used energy over time.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Coefficient of reduction
|<math>K_r</math>
|10
|[-]
|-
|Air density
|<math>\rho</math>
|1.293
|<math> kg/m^3 </math>
|-
|Aerodynamic coefficient
|<math>C_x</math>
|0.4
|[-]
|-
|Area in the front of the vehicle
|<math>S</math>
|2
|<math> m^2 </math>
|-
|Radius of the wheel
|<math>r</math>
|0.33
|<math> m </math>
|-
|Constant representing the friction of the wheels on the road
|<math>K_f </math>
|0.03
|[-]
|-
|Coefficient of the motor torque
|<math>K_m</math>
|0.27
|[-]
|-
|Inductor resistance
|<math>R_m</math>
|0.03
|<math> Ohms </math>
|-
|inductance of the rotor
|<math>L_m</math>
|0.05
|[-]
|-
|Mass
|<math>M</math>
|250
|<math> kg </math>
|-
|Gravity constant
|<math>g</math>
|9,81
|[-]
|-
|Battery voltage
|<math>V_{alim}</math>
|150
|<math> V </math>
|-
|Resistance of the battery
|<math>R_{bat}</math>
|0.05
|<math> Ohms </math>
|-
|}

== Reference Solutions ==
We look at the particular instance of our problem with <math> t_f = 10s </math>and target set <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R} </math>, in which the car needs to cover 100m in 10s.
Figure 1 shows a plot of the differential states of the optimal trajectory of the relaxed problem (i.e. <math> u \in [-1,1] </math>instead of <math> u \in \{-1,1\} </math>) for <math> N = 1000, N </math>being the number of time discretization points. The current <math> x_0 </math>increases to its maximal value of 150A, stays there for a certain time, decreases on its minimal value of -150A, stays on this value and eventually increases slightly. This behavior corresponds to the different arcs bang, path-constrained, singular, path-constrained, bang and can be observed also in Figure 2. It shows the corresponding switching function and the optimal control. Note that the plots show data from the solution with the indirect approach.

Applying the sum up rounding strategy results in an integer-feasible chattering solution.
The resulting primal states are shown in Figure 3. One observes the high-frequency zig-zagging of the current <math> x_0 </math>that results from the fast switches in the control.

The direct and indirect approaches are local optimization techniques and only provide upper bounds for the relaxed problem and hence for the original problem. Here the indirect solution of the relaxed problem gives us a bound of <math> x_3(t_f) = 22777.2 </math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3">
Image:ElectricCar PrimalStates Plot.PNG| Primal states of an optimal trajectory for the relaxed problem on a control discretization grid with <math> N = 1,000 </math>.
Image:ElectricCar Control PathConstraints SwitchingFunction Plot.PNG| The optimal control and switching function. The dotted vertical lines show the switching times <math> \tau_i </math> where transitions between different kinds of arcs occur.
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math>
* bounded velocity, <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R}, t_f = 10s, x_2(t) \leq 45 \frac{K_r}{3.6r} \forall t </math>
* fixed final velocity, bounded velocity, longer time horizon, <math> x_2(t) \leq 30 \frac{K_r}{3.6r} \forall t, \mathcal{T} = \mathbb{R} \times \{ 30 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 15s </math>.

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

Template:Current News

2023-10-19T12:27:14Z

SebastianSager:

{| style="background:#EEEEFF;"

{{News|2023/10/19| Added new category [[:Category:GIOC | generalized inverse optimal control]]}}
{{News|2023/10/19| Added [[LinearMetabolic]] problem}}
{{News|2019/3/14| Added multiple control problems with [[:Category: Gekko | Python GEKKO]] implementation}}
{{News|2018/9/12| Added [[Control of Transmission Lines]] problem}}

|}

<noinclude>
== Older news ==

{| style="background:#eaecd0;"

{{News|2016/7/26 | Added [[Industrial robot]] problem}}
{{News|2016/7/24 | Added [[Continuously Stirred Tank Reactor problem | CSTR problem]]}}
{{News|2016/6/30 | Added [[Electric Car]] problem}}
{{News|2016/5/5| Added multiple control problems with [[:Category: AMPL/TACO | AMPL with TACO]] implementation}}
{{News|2016/2/23|Added [[Control of Heat Equation with Actuator Placement | Actuator Placement]] control problem and [[:Category:Elliptic | Elliptic]], [[:Category:Parabolic | Parabolic]] and [[:Category:Hyperbolic | Hyperbolic]] categories}}
{{News|2016/1/10|Added [[Van der Pol Oscillator]], [[Batch reactor]], [[Cushioned Oscillation]] and [[Goddart's rocket problem]] control problems}}
{{News|2015/11/10|Moved mintoc.de to a new server}}
{{News|2012/09/01|Added the [[Lotka Experimental Design]] and a [[:Category:Optimum Experimental Design | category for experimental design problems]]}}
{{News|2011/09/29|Added the [[:Category:AMPL/TACO | first set of AMPL optimal control problems]] using the TACO toolkit}}
{{News|2010/11/21|Added New York [[Subway ride]] control problem}}
{{News|2010/11/18|Extended description of [[:Category:Problem characterization | problem characterization]]}}
{{News|2010/11/18|Description of benchmark library as [http://mathopt.de/PUBLICATIONS/Sager2011b.pdf pdf file] preprint}}
{{News|2010/08/16|Added [[Bang-bang approximation of a traveling wave]] 1D PDE example}}
{{News|2010/02/11|Launch of EU project [http://embocon.org embocon.org]}}
{{News|2009/11/20|Added [[External Links]] page}}
{{News|2009/10/26|Revision and Correction of several [[:Category:Optimica | optimica]] models}}
{{News|2009/08/26|Knitro solution for [[F-8 aircraft (AMPL)]]}}
{{News|2009/08/11|[[F-8 aircraft]] new local minimum found with [[:Category:Optimica | optimica]]/ipopt}}
{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}

|}

== Adding items ==
If you want to add events you have to edit this page. At the top of the page you will see something like this:
<source lang="text">
{| style="background:#EEEEFF;"

{{News|2008/08/14, 15:15|Presentation of the restructured group wiki in the group meeting, room 432.}}

|}
</source>

To add a new item you have to enter a line of the form:
<source lang="text">
{{News|date|content}}
</source>
where <code>date</code> is the date and <code>content</code> is the subject and further information. The field <code>News</code> is mandatory. Please write the <code>date</code> in the format '''YYYY/MM/DD, HH:MM''' to keep it nice and clean. The field <code>content</code> can contain any information you like (but please keep it short) and even wiki syntax is possible.

''Please keep the items in a chronological order and remove older ones to keep the list short.''

[[Category:News]]

</noinclude>

Category:GIOC

2023-10-19T12:24:08Z

SebastianSager: Initial setup of gIOC class

This category includes all generalized inverse optimal control problems.
To formalize this problem class, we define the following bilevel problem with differential states <math>x \in \mathcal{X}</math>, controls <math>u \in \mathcal{U}</math>, model parameters <math>p \in \mathbb{R}^{n_p}</math>, and convex multipliers <math>w \in \mathcal{W}</math>

<math>
\displaystyle \min_{(p, w, x^*, u^*) \in \Omega_1} \; \| h(x^*, u^*) - \eta \| + R(p,w)
</math>
subject to
<math>
(x^*,u^*) \in \displaystyle \arg \min_{x, u} \sum_{i \in M_1} w_i \; \phi_i[x,u,p]
</math>
subject to
<math>
\begin{array}{rcl}
\dot{x}(t) &=& \displaystyle \sum_{i \in M_2} w_i f_i(x(t), u(t), p) \\
0 &\le& \displaystyle w_i \; g_i(x(t), u(t), p) \quad \forall \; i \in M_3\\
(x,u,p,w) &\in& \Omega_2
\end{array}
</math>

as a generalized inverse optimal control problem. Here <math>\mathcal{X}</math> and <math>\mathcal{U}</math> are properly defined function spaces. The variables <math>w</math> indicate which objective functions, right hand side functions, and constraints are relevant in the inner problem. The variables are normalized for given index sets <math>M_1, M_2, M_3</math> that partition the indices from <math>1</math> to <math>n_w</math>. For normalization, we define the feasible set <math>\mathcal{W} := \{ w \in [0,1]^{n_w}: \; \textstyle \sum^{i \in M_j} w_i = 1 \text{ for } j \in \{1,2,3\}\}</math>. On the outer level, the feasible set is <math>\Omega_1 := \mathbb{R}^{n_p} \times \mathcal{W} \times \mathcal{X} \times \mathcal{U}</math>, while on the inner level <math>\Omega_2</math> contains bounds, boundary conditions, mixed path and control constraints, and more involved constraints such as dwell time constraints. We have observational data <math>\eta \in \mathbb{R}^{n_\eta}</math>, a measurement function <math>h: \mathcal{X} \times \mathcal{U} \mapsto \mathbb{R}^{n_\eta}</math>, a regularization function with a priori knowledge on parameters and weights <math>R: \mathbb{R}^{n_p} \times \mathcal{W} \mapsto \mathbb{R}</math> and candidate functionals <math>\phi_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}</math> and functions <math>f_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}^{n_x}</math> and <math>g_i: \mathcal{X} \times \mathcal {U} \times \mathbb{R}^{n_p} \mapsto \mathbb{R}^{n_g}</math> for the unknown objective function, dynamics, and constraints, respectively. Two cases are of practical interest: first, the manual, often cumbersome and trial-and-error based a priori definition of all candidates <math>\phi_i, f_i, g_i</math> by experts and second, a systematic, but often challenging automatic symbolic regression of these unknown functions.

On the outer level, a norm <math>\| \cdot \|</math> and the regularization term <math>R</math> define a data fit (regression) problem and relate to prior knowledge and statistical assumptions. On the inner level, the above bilevel problem is constrained by a possibly nonconvex optimal control problem. The unknown parts of this inner level optimal control problem are modeled as convex combinations of a finite set of candidates (and a multiplication of constraints <math>g_i</math> with <math>w_i</math> that can be either zero or strictly positive). On the one hand the problem formulation is restrictive in the interest of a clearer presentation and might be further generalized, e.g., to multi-stage formulations involving differential-algebraic or partial differential equations. On the other hand, the problem class is quite generic and allows, e.g., the consideration of switched systems, periodic processes, different underlying function spaces <math>\mathcal{X}</math> and <math>\mathcal{U}</math>, and the usage of universal approximators such as neural networks as candidate functions.

== Extensions ==
* For some problems the functions may as well depend explicitely on the time <math>t</math>.
* The differential equations might depend on [[:Category:State dependent switches | state-dependent switches]].
* The variables may include [[:Category:Boolean variables | boolean variables]].
* The underlying process might be a [[:Category:Multistage process | multistage process]].
* The dynamics might be [[:Category:Unstable | unstable]].
* There might be an underlying [[:Category:Network topology | network topology]].
* The integer control functions might have been (re)formulated by means of an [[:Category:Outer convexification|outer convexification]].

Note that a Lagrange term <math>\int_{t_0}^{t_f} L( x(t), u(t), v(t), q, \rho) \; \mathrm{d} t</math> can be transformed into a Mayer-type objective functional.

[[Category:Model characterization]]

LinearMetabolic

2023-10-19T11:33:00Z

SebastianSager: Created page with "{{Dimensions |nd = 1 |nx = 4 |nw = 1 |nc = 2 |nri = 2 |nre = 5 }}<!-- Do not insert line break here or Dimensions Box moves up in the..."

{{Dimensions
|nd = 1
|nx = 4
|nw = 1
|nc = 2
|nri = 2
|nre = 5
}}

Under construction...

The '''Linear Metabolic problem''' is a generalized inverse optimal control problem. It tries to identify an a priori unknown objective function from data. The control unction <math> u(t) </math> is ...

The problem is discussed in detail in <bib id="Sager2015" />.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\
& \dot{x}_1 & = & \frac{K_r^2}{Mr^2} (K_m x_0 - \frac{r}{K_r} ( M g K_f + \frac{1}{2} \rho S C_x \frac{r^2}{K_r^2} x_1^2)), \\
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]
& x(t_0) &=& (0,0,0,0)^T, \\
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\
& x_0(t) & \in & [-i_{max}, i_{max}], \\
& u(t) &\in& \{-1, 1\}.
\end{array}
</math>

Here the four differential states stand for the electrical current (<math> x_0 </math>), the angular velocity (<math> x_1 </math>), the position of the car (<math> x_2 </math>), and the consumed energy (<math> x_3 </math>). The objective function <math> x_3(t_f) </math>is just a reformulation of the Lagrange-type objective function tracking the used energy over time.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Coefficient of reduction
|<math>K_r</math>
|10
|[-]
|-
|Air density
|<math>\rho</math>
|1.293
|<math> kg/m^3 </math>
|-
|Aerodynamic coefficient
|<math>C_x</math>
|0.4
|[-]
|-
|Area in the front of the vehicle
|<math>S</math>
|2
|<math> m^2 </math>
|-
|Radius of the wheel
|<math>r</math>
|0.33
|<math> m </math>
|-
|Constant representing the friction of the wheels on the road
|<math>K_f </math>
|0.03
|[-]
|-
|Coefficient of the motor torque
|<math>K_m</math>
|0.27
|[-]
|-
|Inductor resistance
|<math>R_m</math>
|0.03
|<math> Ohms </math>
|-
|inductance of the rotor
|<math>L_m</math>
|0.05
|[-]
|-
|Mass
|<math>M</math>
|250
|<math> kg </math>
|-
|Gravity constant
|<math>g</math>
|9,81
|[-]
|-
|Battery voltage
|<math>V_{alim}</math>
|150
|<math> V </math>
|-
|Resistance of the battery
|<math>R_{bat}</math>
|0.05
|<math> Ohms </math>
|-
|}

== Reference Solutions ==
We look at the particular instance of our problem with <math> t_f = 10s </math>and target set <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R} </math>, in which the car needs to cover 100m in 10s.
Figure 1 shows a plot of the differential states of the optimal trajectory of the relaxed problem (i.e. <math> u \in [-1,1] </math>instead of <math> u \in \{-1,1\} </math>) for <math> N = 1000, N </math>being the number of time discretization points. The current <math> x_0 </math>increases to its maximal value of 150A, stays there for a certain time, decreases on its minimal value of -150A, stays on this value and eventually increases slightly. This behavior corresponds to the different arcs bang, path-constrained, singular, path-constrained, bang and can be observed also in Figure 2. It shows the corresponding switching function and the optimal control. Note that the plots show data from the solution with the indirect approach.

Applying the sum up rounding strategy results in an integer-feasible chattering solution.
The resulting primal states are shown in Figure 3. One observes the high-frequency zig-zagging of the current <math> x_0 </math>that results from the fast switches in the control.

The direct and indirect approaches are local optimization techniques and only provide upper bounds for the relaxed problem and hence for the original problem. Here the indirect solution of the relaxed problem gives us a bound of <math> x_3(t_f) = 22777.2 </math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3">
Image:ElectricCar PrimalStates Plot.PNG| Primal states of an optimal trajectory for the relaxed problem on a control discretization grid with <math> N = 1,000 </math>.
Image:ElectricCar Control PathConstraints SwitchingFunction Plot.PNG| The optimal control and switching function. The dotted vertical lines show the switching times <math> \tau_i </math> where transitions between different kinds of arcs occur.
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math>
* bounded velocity, <math> \mathcal{T} = \mathbb{R} \times \mathbb{R} \times \{100\} \times \mathbb{R}, t_f = 10s, x_2(t) \leq 45 \frac{K_r}{3.6r} \forall t </math>
* fixed final velocity, bounded velocity, longer time horizon, <math> x_2(t) \leq 30 \frac{K_r}{3.6r} \forall t, \mathcal{T} = \mathbb{R} \times \{ 30 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 15s </math>.

== References ==
<biblist />


[[Category:gIOC]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Bang bang]]
[[Category: Path-constrained arcs]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

External Links

2022-03-02T14:38:48Z

SebastianSager: /* Work Groups */

= Work Groups =

An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.

* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson
* [https://www.tu-braunschweig.de/mo/team/kirches Mathematical Optimization], led by Christian Kirches
* [http://www.mathematik.tu-clausthal.de/arbeitsgruppen/kontinuierliche-optimierung/ Continuous Optimization], led by Andreas Potschka
* [https://www.mathematik.hu-berlin.de/de/forschung/forschungsgebiete/mathematische-optimierung/Falk-hante/falk-hante-fp Applied Mathematics], led by Falk Hante

Feel free to update this list!

= Benchmark Libraries =

=== Embedded Optimization and Control ===

* European FP 7 project [http://embocon.org EMBOCON]

=== Hybrid Systems ===

* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.

=== Optimal Control ===

* The [http://tomdyn.com/ PROPT homepage] (a matlab toolkit for dynamic optimization using collocation) states over 100 test cases from different applications with their results and computation time.
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems.
* The ESA provides a [http://www.esa.int/gsp/ACT/inf/op/globopt.htm test set of global optimization spacecraft trajectory problems] and their best putative solutions.
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization
* The [http://plato.asu.edu/pdecon.html PDECON collection] contains discretized PDE problems in AMPL as part of a [http://plato.asu.edu/bench.html larger benchmark effort] by Hans Mittelmann

=== MINLP ===

* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]
* MINLPlib <bib id="Bussieck2003a" />
* [http://www.gamsworld.org/performance GAMSWORLD]

=== Generic Optimization ===

* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)

== References ==
<biblist />

Oil Shale Pyrolysis

2016-02-04T19:40:52Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nw = 0
|nre = 2
}}

The following problem is an example from the global optimal control literature and was introduced in <bib id="Wen1977" />.

== Mathematical formulation ==

<math>
\begin{array}{ll}
\displaystyle \min_{u} & \displaystyle -x_1(t_N)^2 \\[1.5ex]
\mbox{s.t.} & \displaystyle \dot{x}_0(t) = -k_0x_0(t)-(k_2+k_3+k_4)x_0(t)x_1(t)\\
& \displaystyle \dot{x}_1(t) = k_0x_0(t)-k_1x_1(t) + k_2x_0(t)x_1(t)\\
& \displaystyle k_i = a_i e^{-u(t)\frac{b_i}{R}},\quad \forall i\in \{1,\dots,5\} \\ [1.5ex]
& \displaystyle t \in \left[t_0,t_N\right] \\
& \displaystyle u(t) \in \left[698.15/748.15,1\right]\\
& \displaystyle x(t_0) = (1,0)^T\\
\end{array}
</math>

where this is the normalized form with

<math> u(t)= \frac{1}{u_{temp}} </math>, with

<math> u_{temp} \in \left[698.15,748.15\right] </math>

== Parameters ==

{| class="wikitable"
|+State variables
|-
|Symbol
|Initial value (<math>t_0</math>)
|-
|<math>x_0(t)</math>
|<math>1</math>
|-
|<math>x_1(t)</math>
|<math>0 </math>
|}

{| class="wikitable"
|+Parameters
|-
|Symbol
|Value
|-
|<math>a_1</math>
|<math>8.86</math>
|-
|<math>a_2</math>
|<math>24.25</math>
|-
|<math>a_3</math>
|<math>23.67</math>
|-
|<math>a_4</math>
|<math>18.75</math>
|-
|<math>a_5</math>
|<math>20.7</math>
|-
|<math>b_1</math>
|<math>20.3</math>
|-
|<math>b_2</math>
|<math>37.4</math>
|-
|<math>b_3</math>
|<math>33.8</math>
|-
|<math>b_4</math>
|<math>28.2</math>
|-
|<math>b_5</math>
|<math>31.0</math>
|}

{| class="wikitable"
|+Control variable
|-
|Symbol
|Interval
|-
|<math>u(t)</math>
|[698.15/748.15,1]
|}

'''Measurement grid'''

== Reference solution ==

Coming soon.

== Source Code ==

Model descriptions are not yet available.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]

Lotka Volterra fishing problem

2016-02-04T19:38:59Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_0(t) - x_0(t) x_1(t) - \; c_0 x_0(t) \; w(t), \\
& \dot{x}_1(t) & = & - x_1(t) + x_0(t) x_1(t) - \; c_1 x_1(t) \; w(t), \\
& \dot{x}_2(t) & = & (x_0(t) - 1)^2 + (x_1(t) - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:optimica | optimica code]] at [[Lotka Volterra fishing problem (optimica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>,
* different fishing control functions for the two species,
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

External Links

2016-01-28T10:09:13Z

SebastianSager:

= Work Groups =

An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.

* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson

Feel free to update this list!

= Benchmark Libraries =

=== Embedded Optimization and Control ===

* European FP 7 project [http://embocon.org EMBOCON]

=== Hybrid Systems ===

* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.

=== Optimal Control ===

* The [http://tomdyn.com/ PROPT homepage] (a matlab toolkit for dynamic optimization using collocation) states over 100 test cases from different applications with their results and computation time.
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems.
* The ESA provides a [http://www.esa.int/gsp/ACT/inf/op/globopt.htm test set of global optimization spacecraft trajectory problems] and their best putative solutions.
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization
* The [http://plato.asu.edu/pdecon.html PDECON collection] contains discretized PDE problems in AMPL as part of a [http://plato.asu.edu/bench.html larger benchmark effort] by Hans Mittelmann

=== MINLP ===

* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]
* MINLPlib <bib id="Bussieck2003a" />
* [http://www.gamsworld.org/performance GAMSWORLD]

=== Generic Optimization ===

* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)

== References ==
<biblist />

Main Page

2016-01-28T10:01:22Z

SebastianSager:

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]]'' - ''[[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

Main Page

2016-01-28T09:53:28Z

SebastianSager:

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category: Muscod | Muscod]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]]'' - ''[[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

Annihilation of calcium oscillations with PLC activation inhibition

2016-01-28T09:51:48Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 4
|nw = 2
|nre = 4
}}

This control problem is closely related to [[Annihilation of calcium oscillations]]. The only difference is an additional control function, the inhibition of PLC activation. We state only the differences in this article.

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, w^{\mathrm{max}}} & & & {\int_{t_0}^{t_f} || x(\tau) - \tilde{x} ||_2^2 + p_1 w_1(\tau) + p_2 w_2(\tau) \; \mathrm{d}\tau} \\[1.5ex]
\mbox{s.t.} & \dot{x}_0 & = & k_1 + k_2 x_0 - \frac{k_3 x_0 x_1}{x_0 + K_4} - \frac{k_5 x_0 x_2}{x_0 + K_6} \\
& \dot{x}_1 & = & (1 - w_2) k_7 x_0 - \frac{k_8 x_1}{x_1 + K_9} \\
& \dot{x}_2 & = & \frac{k_{10} x_1 x_2 x_3}{x_3 + K_{11}} + k_{12} x_1 + k_{13} x_0 - \frac{k_{14} x_2}{w_1 \cdot x_2 + K_{15}} - \frac{k_{16} x_2}{x_2 + K_{17}} + \frac{x_3}{10} \\
& \dot{x}_3 & = & - \frac{k_{10} x_1 x_2 x_3}{x_3 + K_{11}} + \frac{k_{16} x_2}{x_2 + K_{17}} - \frac{x_3}{10} \\[1.5ex]
& x(0) &=& (0.03966, 1.09799, 0.00142, 1.65431)^T, \\
& x(t) & \ge & 0.0, \\
& w_1(t) &\in& \{1, w^{\mathrm{max}}\}, \\
& w_2(t) &\in& \{0, 1\}, \\
& w^{\mathrm{max}} & \ge & 1.1, \\
& w^{\mathrm{max}} & \le & 1.3.
\end{array}
</math>

Note that we write <math>w_1(t)</math> instead of <math>w(t)</math> and have an additional control function <math>w_2(t)</math>. The regularization parameters are set to <math>p_1 = p_2 = 100</math>.

== Reference Solutions ==

A solution for this problem is described in <bib id="Lebiedz2005" />. A local minimum that is actually slightly worse than the solution provided for only one control, is shown in the next plots.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:calcium2Controls.png| Optimal stimuli determined by mixed-integer optimal control.
Image:calcium2States.png| Corresponding differential states.
</gallery>

== Variants ==

[[Annihilation of calcium oscillations| Only one control function]].

== References ==
<biblist />


[[Category:Medicine]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Unstable]]
[[Category:Tracking objective]]
[[Category:Bang bang]]
[[Category:Systems biology]]

Annihilation of calcium oscillations

2016-01-28T09:51:36Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 4
|nw = 1
|nre = 4
}}

The aim of the control problem is to identify strength and timing of inhibitor stimuli that lead to a phase singularity which annihilates intracellular calcium oscillations. This is formulated as an objective function that aims at minimizing the state deviation from a desired unstable steady state, integrated over time.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states. The problem is, despite of its low dimension, very hard to solve, as the target state is [[:Category:Unstable | unstable]].

== Biological motivation ==

Biological rhythms as impressing manifestations of self-organized dynamics associated with the phenomenon life have been of particular interest since quite a long time. Even before the mechanistic basis of certain biochemical oscillators was elucidated by molecular biology techniques, their investigation and the issue of perturbation by external stimuli has attracted much attention.

A calcium oscillator model describing intracellular calcium spiking in hepatocytes induced by an extracellular increase in adenosine triphosphate (ATP) concentration is described. The calcium signaling pathway is initiated via a receptor activated G-protein inducing the intracellular release of inositol triphoshate (IP3) by phospholipase C. The IP3 triggers the opening of endoplasmic reticulum and plasma membrane calcium channels and a subsequent inflow of calcium ions from intracellular and extracellular stores leading to transient calcium spikes.

As a source of external control a temporally varying concentration of an uncompetitive inhibitor of the PMCA ion pump is considered.

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, w^{\mathrm{max}}} & & & {\int_{t_0}^{t_f} || x(\tau) - \tilde{x} ||_2^2 + p_1 w(\tau) \; \mathrm{d}\tau} \\[1.5ex]
\mbox{s.t.} & \dot{x}_0 & = & k_1 + k_2 x_0 - \frac{k_3 x_0 x_1}{x_0 + K_4} - \frac{k_5 x_0 x_2}{x_0 + K_6} \\
& \dot{x}_1 & = & k_7 x_0 - \frac{k_8 x_1}{x_1 + K_9} \\
& \dot{x}_2 & = & \frac{k_{10} x_1 x_2 x_3}{x_3 + K_{11}} + k_{12} x_1 + k_{13} x_0 - \frac{k_{14} x_2}{w \cdot x_2 + K_{15}} - \frac{k_{16} x_2}{x_2 + K_{17}} + \frac{x_3}{10} \\
& \dot{x}_3 & = & - \frac{k_{10} x_1 x_2 x_3}{x_3 + K_{11}} + \frac{k_{16} x_2}{x_2 + K_{17}} - \frac{x_3}{10} \\[1.5ex]
& x(0) &=& (0.03966, 1.09799, 0.00142, 1.65431)^T, \\
& x(t) & \ge & 0.0, \\
& w(t) &\in& \{1, w^{\mathrm{max}}\}, \\
& w^{\mathrm{max}} & \ge & 1.1, \\
& w^{\mathrm{max}} & \le & 1.3.
\end{array}
</math>


Here the differential states <math>(x_0, x_1, x_2, x_3)</math> describe concentrations of activated G-proteins, active phospholipase C, intracellular calcium, and intra-ER calcium, respectively.

Modeling details including a comprehensive discussion of parameter values and the dynamical behavior observed in simulations with a comparison to experimental observations can be found in <bib id="Kummer2000" />. In the given equations that stem from <bib id="Lebiedz2005" />, the model is identical to the one derived there, except for an additional first-order leakage flow of calcium from the ER back to the cytoplasm, which is modeled by <math>\frac{x_3}{10}</math> in equations 3 and 4. It reproduces well experimental observations of cytoplasmic calcium oscillations as well as bursting behavior and in particular the frequency encoding of the triggering stimulus strength, which is a well known mechanism for signal processing in cell biology.

== Parameters ==

These fixed values are used within the model.


<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 22],\\
k_1 &=& 0.09, \\
k_2 &=& 2.30066, \\
k_3 &=& 0.64, \\
K_4 &=& 0.19, \\
k_5 &=& 4.88, \\
K_6 &=& 1.18, \\
k_7 &=& 2.08, \\
k_8 &=& 32.24, \\
K_9 &=& 29.09, \\
k_{10} &=& 5.0, \\
K_{11} &=& 2.67, \\
k_{12} &=& 0.7, \\
k_{13} &=& 13.58, \\
k_{14} &=& 153.0, \\
K_{15} &=& 0.16, \\
k_{16} &=& 4.85, \\
K_{17} &=& 0.05, \\
p_1 &=& 100, \\
\tilde{x}_0 &=& 6.78677, \\
\tilde{x}_1 &=& 22.65836, \\
\tilde{x}_2 &=& 0.384306, \\
\tilde{x}_3 &=& 0.28977.
\end{array}
</math>


== Reference Solutions ==
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:calciumControl.png| Optimal stimuli determined by mixed-integer optimal control.
Image:calciumStates.png| Corresponding differential states.
Image:calciumStatesPerturbed.png| Slightly perturbed control: stimulus 0.001 too early.
Image:calciumStatesLongTime.png| Long time behavior of optimal solution: numerical rounding errors lead to transition back from unstable steady-state to stable limit-cycle.
</gallery>

The depicted optimal solution consists of a stimulus of <math>w^{\mathrm{max}}=1.3</math> and a timing given by the stage lengths 4.6947115, 0.1491038, and 17.1561845. The optimal objective function value is 1610.654. As can be seen from the additional plots, this solution is extremely unstable. A small perturbation in the control, or simply rounding errors on a longer time horizon lead to a transition back to the stable limit-cycle oscillations.

== Optimization ==

The determination of the stimulus by means of optimization is quite hard for two reasons. First, the unstable target steady-state. Only a stable all-at-once algorithm such as multiple shooting or collocation can be applied successfully.

Second, the objective landscape of the problem in switching time formulation (this is, for a fixed stimulus strength and modifying only beginning and length of the stimulus) is quite nasty, as the following visualizations obtained by simulation show.

<gallery caption="Objective landscape simulation plots" widths="180px" heights="140px" perrow="2">
Image:calciumObjectiveSimulation1.png| Brute-force simulation of objective landscape for different values of stimulus beginning and length.
Image:calciumObjectiveSimulation2.png| Zoom in.
Image:calciumObjectiveSimulation3.png| Zoom in.
Image:calciumObjectiveSimulation4.png| Note that the border line of the small channel is a local maximum, hence an initialization slightly away will not lead to the correct global minimum.
</gallery>

== Source Code ==

* [[:Category:Muscod | Muscod code]] at [[Annihilation of calcium oscillations (Muscod)]]
* [[:Category:optimica | optimica]] at [[Annihilation of calcium oscillations (optimica)]]

== Variants ==

<ul>
<li> [[Annihilation of calcium oscillations with PLC activation inhibition| Use of two control functions]].
<li> Scaled deviation in objective function.
</ul>

== References ==
<biblist />


[[Category:Medicine]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Unstable]]
[[Category:Bang bang]]
[[Category:Systems biology]]
[[Category:Tracking objective]]

De Pillis chemotherapy model

2016-01-28T09:51:12Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 6
|nu = 3
}}

The model by de Pillis combines chemotherapy with immunotherapy.

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> and <math>D = d \frac{(x_2/x_0)^l}{s+(x_2/x_0)^l}</math> the optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \min_{x, u} & p_0 x_0(t_f) & + & \int_{t_0}^{t_f} p_1 x_0(t)^2 \mbox{d}t + \sum_{i=0}^{3} \int_{t_0}^{t_f} p_{i+2} u_i(t)^2 \mbox{d}t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & a x_0 (1-b x_0) -c x_1 x_0 - D x_0 - K_T (1- \mbox{e}^{- x_4}) x_0 ,\\
& \dot{x}_1(t) & = & e x_3 - f x_1 + g \frac{x_0^2}{h+x_0^2}-p x_1 x_0 - K_N (1- \mbox{e}^{- x_4}) x_1, \\
& \dot{x}_2(t) & = & -m x_2 + j \frac{D^2 x_0^2}{k+ D^2 x_0^2} x_2 - q x_1 x_2 + (r_1 x_1 + r_2 x_3) x_0 \\
& & & - v x_1 x_2^2 - K_L (1- \mbox{e}^{- x_4}) x_2 + \frac{p_I x_2 x_5}{g_I + x_5} + u_2, \\
& \dot{x}_3(t) & = & \alpha - \beta x_3 - K_C (1- \mbox{e}^{- x_4}) x_3,\\
& \dot{x}_4(t) & = & - \gamma x_4 + u_0,\\
& \dot{x}_5(t) & = & - \mu_I x_5 + u_1.
\end{array}
</math>


The states <math>x_0</math> to <math>x_3</math> are measured in absolute cell counts, where <math>x_0</math> describes the number of tumor cells, <math>x_1</math> of unspecific immune cells, <math>x_2</math> of tumor-specific cytotoxic T-cells (CD<math>8^+</math> T) and <math>x_3</math> of circulating lymphocytes. The chemotherapeutic drug concentration is given by <math>x_4</math> and the immunotherapeutic by <math>x_5</math> (Interleukin-2) respectively.

== Parameters ==

This set of parameters can be found as “patient 9” in [].

<math>
\begin{array}{lll}
a = 4.31 \times 10^{-1} & b = 1.02 \times 10^{-9} & c = 6.41 \times 10^{-11}\\
d = 2.34 & e = 2.08 \times 10^{-7} & f = 4.12 \times 10^{-2}\\
g = 1.25 \times 10^{-2} & h = 2.02 \times 10^{7} & j = 2.49 \times 10^{-2}\\
k = 3.66 \times 10^{7} & l = 2.09 & m = 2.04 \times 10^{-1}\\
q = 1.42 \times 10^{-6} & p = 3.42 \times 10^{-6} & s = 8.39 \times 10^{-2}\\
r_1 = 1.01 \times 10^{-7} & r_2 = 6.50 \times 10^{-11} & u = 3.00 \times 10^{-10}\\
\alpha = 7.50 \times 10^{8} & \beta = 1.20 \times 10^{-2} & \gamma = 9.00 \times 10^{-1}\\
p_I = 1.25 \times 10^{-1} & g_I = 2.00 \times 10^{7} & \mu_I = 1.00 \times 10^{1}\\
K_T = 9.00 \times 10^{-1} & K_N = K_L = K_C = 6 \times 10^{-1}
\end{array}
</math>

== Reference Solutions ==

The problem can be solved with the [multiple shooting method]. For the following solutions the control functions and states are discretized on the same grid, with 100 nodes. For the objective function parameters have been chosen from the following sets to grand the results shown in the graphics below.

Objective function 1
<math>
\begin{array}{rclr}
p_0 &=& 1,&\\
p_2 &=& 1,&\\
p_i &=& 0,&otherwise.\\
\end{array}
</math>

Objective function 2
<math>
\begin{array}{rclr}
p_0 &=& -1,&\\
p_2 &=& 1,&\\
p_i &=& 0,&otherwise.\\
\end{array}
</math>

In both objective function the amount of chemotherapeutic drugs is penalized. The objective function 2 describes the worst case scenario of the tumor growth at end time.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:DePillis_1.png| Optimal controls and states computed with a multiple shooting approach on the same discretization grid with 100 nodes and objective function 1.
Image:DePillis_2.png| Optimal controls and states with objective function 2.
</gallery>

==Source Code==

* [[:Category:Muscod | Muscod code]] at [[De Pillis chemotherapy model (Muscod)]]

== References ==

[[Category:Medicine]]
[[Category:MIOCP]]
[[Category:ODE model]]

Category:Medicine

2016-01-28T09:50:15Z

SebastianSager: Initial setup

This category includes all control problems that stem from a medical application.

[[Category:Application]]

D'Onofrio chemotherapy model

2016-01-28T09:49:51Z

SebastianSager: added medicine

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
}}

This cancer chemotherapy model is based on the work of d'Onofrio. The corresponding dynamic describes the effect of two different drugs administered to the patient. An anti-angiogetic drug is used to suppress the formation of blood vessels from existing vessels and thereby starving the tumors supply of proliferating vessels. In addition a cytostatic drug effects the proliferation of the tumor cells directly. The dynamic of the problem is given by an [[:Category:ODE model|ODE model]].

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> the optimal control problem is given by


<math>
\begin{array}{llcl}
\displaystyle \min_{x, u} & x_0(t_f) &+& \alpha \int_{t_0}^{t_f} u_0(t)^2 \text{d}t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & - \zeta x_0(t) \text{ln} \left( \frac{x_0(t)}{x_1(t)} \right) - F \; x_0(t) u_1(t), \\
& \dot{x}_1(t) & = & b x_0(t) - \mu x_1(t) - d x_0(t)^{\frac{2}{3}}x_1(t) - G u_0(t) x_1(t) - \eta x_1(t) u_1(t), \\
& \dot{x}_2(t) & = & u_0(t), \\
& \dot{x}_3(t) & = & u_1(t), \\ [1.5ex]
& 0 & \leq & u_0(t) \leq u_0^{max}, \\
& 0 & \leq & u_1(t) \leq u_1^{max}, \\
& x_2(t) & \leq & x_2^{max}, \\
& x_3(t) & \leq & x_3^{max}.
\end{array}
</math>


where the control <math>u_0</math> denotes the administered amount of anti-angiogetic drugs and <math>u_1</math> the amount of cytostatic drugs. The state <math>x_0</math> describes the volume of tumor and <math>x_1</math> the volume of neighboring blood vessels. The remaining states <math>x_2</math> and <math>x_3</math> are used to constraint the maximum amount of drugs over the duration of the therapy.

== Parameters ==

In the model these parameters are fixed.

<math>
\begin{array}{rcl}
t_0 &=& 0,\\
(\zeta, b, \mu, d, G) &=& (0.192, 5.85, 0.0, 0.00873, 0.15),\\
(x_2(0), x_3(0), u_0^{max}, x_2^{max}) &=& (0,0,75,300).
\end{array}
</math>

The parameters <math>(x_0(0), x_1(0), u_1^{max}, x_3^{max})</math> can be taken from the parameter sets shown in the following section. To the remaining parameters <math>(F, \eta)</math> exists no experimental data.

== Reference Solutions ==

The problem can be solved with the [multiple shooting method]. For the following solutions the control functions and states are discretized on the same grid, with 100 nodes. The unknown parameters are chosen from the following parameter sets

Parameter set 1

<math>
\begin{array}{rclrcl}
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\
\end{array}
</math>

Parameter set 2

<math>
\begin{array}{rclrcl}
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\
\end{array}
</math>

Parameter set 3

<math>
\begin{array}{rclrcl}
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\
\end{array}
</math>

Parameter set 4

<math>
\begin{array}{rclrcl}
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\
\end{array}
</math>

Furthermore in the objective function <math>\alpha =0</math> is chosen.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:d'Onofrio 1.png| Optimal controls and states computed with a multiple shooting approach on the same discretization grid with 100 nodes and parameter set 1.
Image:d'Onofrio 2.png| Optimal controls and states with parameter set 2.
Image:d'Onofrio 3.png| Optimal controls and states with parameter set 3.
Image:d'Onofrio 4.png| Optimal controls and states with parameter set 4.
</gallery>

==Source Code==

[[:Category:Muscod | Muscod code]]
<source lang="cpp">
/* volume of tumor */
rhs[0] = - zeta*x[0]*log(x[0]/x[1])-F*x[0]*u1_;
/* volume of neighboring blood vessels */
rhs[1] = b*x[0] - (mu + d*pow(x[0],2.0/3.0))*x[1] - G*u0_*x[1] - eta*x[1]*u1_;
/* amount of anti-angiogetic drug */
rhs[2] = u0_;
/* amount of cytostatic drug */
rhs[3] = u1_;
</source>

== References ==
<biblist />

[[Category:Medicine]]
[[Category:MIOCP]]
[[Category:ODE model]]

Cushioned Oscillation

2016-01-28T09:41:42Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nre = 4
}}The Cushioned Oscillation is a simplified model of time optimal "stopping" of an oscillating object attached to a spring by applying a control and moving it back into the relaxed position and zero velocity.

== Model formulation ==

An object with mass <math> m </math> is attached to a spring with stiffness constant <math> c </math>.

If the resetting spring force is proportional to the deviation <math>x=x(t)</math>, an oscillation, induced by an external force <math>u(t)</math>, satisfies:

<math> m\dot v (t) + cx(t) = u(t)</math> (which is equivalent to <math>\dot v (t) = \frac{1}{m}(u(t) - cx(t))</math>)

where <math>x(t)</math> denotes the deviation to the relaxed position and <math> v(t)=\dot x (t) </math> the velocity of the oscillating object.

Through external force, the object has been put into an initial state :

<math>(x(0),v(0)) = (x_0,v_0)</math>

The goal is to reset position and velocity of the object as fast as possible, meaning:

<math>(x(t_f),v(t_f)) = (0,0)</math>,

with the objective function:

<math>\min\limits_{t_f} t_f</math>

== Optimal Control Problem (OCP) Formulation ==

The above results in the following OCP

<math> \begin{array}{llr}
\min\limits_{x,v,u,t_f} & t_f\\

s.t. & \dot x (t) = v(t), & \forall t \in [0,t_f]\\

& \dot v (t)= \frac{1}{m}(u(t) - cx(t)), & \forall t \in [0,t_f]\\
\\
& x(0)=x_0,\\
& v(0)=v_0,\\
& x(t_f)=0,\\
& v(t_f)=0,\\
\\
& -u_{mm} \le u(t) \le u_{mm}, & \forall t \in [0,t_f]\\
\end{array}</math>

== Parameters and Reference Solution ==

The following parameters were used, to create the reference solution below, with an almost optimal final time <math> t_f = 8.98 s</math>:

<math> m=5, </math>
<math> c=10, </math>
<math> x_0=2, </math>
<math> v_0=5, </math>
<math> u_{mm}=5.</math>

== Reference Solution ==

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="1">
Image:Ref_sol_plot_cushioned_oscillation_m5.png| States and Controls
</gallery>

The OCP was solved within MATLAB R2015b, using the TOMLAB Optimization Package. PROPT reformulates such problems with the direct collocation approach (n=80 collocation points) and automatically finds a suiting solver included in the TOMLAB Optimization Package (in this case, SNOPT was used).

== Source Code ==

* A MATLAB script using [[:Category:TomDyn/PROPT | PROPT]] can be found in: [[Cushioned Oscillation (PROPT)]]

== References ==

[[Category:MIOCP]]
[[Category: Minimum time]]

Cushioned Oscillation

2016-01-28T09:40:21Z

SebastianSager:

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nre = 4
}}The Cushioned Oscillation is a simplified model of time optimal "stopping" of an oscillating object attached to a spring by applying a control and moving it back into the relaxed position and zero velocity.

== Model formulation ==

An object with mass <math> m </math> is attached to a spring with stiffness constant <math> c </math>.

If the resetting spring force is proportional to the deviation <math>x=x(t)</math>, an oscillation, induced by an external force <math>u(t)</math>, satisfies:

<math> m\dot v (t) + cx(t) = u(t)</math> (which is equivalent to <math>\dot v (t) = \frac{1}{m}(u(t) - cx(t))</math>)

where <math>x(t)</math> denotes the deviation to the relaxed position and <math> v(t)=\dot x (t) </math> the velocity of the oscillating object.

Through external force, the object has been put into an initial state :

<math>(x(0),v(0)) = (x_0,v_0)</math>

The goal is to reset position and velocity of the object as fast as possible, meaning:

<math>(x(t_f),v(t_f)) = (0,0)</math>,

with the objective function:

<math>\min\limits_{t_f} t_f</math>

== Optimal Control Problem (OCP) Formulation ==

The above results in the following OCP

<math> \begin{array}{llr}
\min\limits_{x,v,u,t_f} & t_f\\

s.t. & \dot x (t) = v(t), & \forall t \in [0,t_f]\\

& \dot v (t)= \frac{1}{m}(u(t) - cx(t)), & \forall t \in [0,t_f]\\
\\
& x(0)=x_0,\\
& v(0)=v_0,\\
& x(t_f)=0,\\
& v(t_f)=0,\\
\\
& -u_{mm} \le u(t) \le u_{mm}, & \forall t \in [0,t_f]\\
\end{array}</math>

== Parameters and Reference Solution ==

The following parameters were used, to create the reference solution below, with an almost optimal final time <math> t_f = 8.98 s</math>:

<math> m=5, </math>
<math> c=10, </math>
<math> x_0=2, </math>
<math> v_0=5, </math>
<math> u_{mm}=5.</math>

== Reference Solution ==

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="1">
Image:Ref_sol_plot_cushioned_oscillation_m5.png| States and Controls
</gallery>

The OCP was solved within MATLAB R2015b, using the TOMLAB Optimization Package. PROPT reformulates such problems with the direct collocation approach (n=80 collocation points) and automatically finds a suiting solver included in the TOMLAB Optimization Package (in this case, SNOPT was used).

== Source Code ==

* A MATLAB script using [[Category:TomDyn/PROPT]] can be found in: [[Cushioned Oscillation (PROPT)]]

== References ==

[[Category:MIOCP]]
[[Category: Minimum time]]

Main Page

2016-01-28T09:33:52Z

SebastianSager: Moved C to Muscod

__NOTOC__
This wiki contains a '''benchmark library''' of '''mixed-integer optimal control problems'''. The main intention is to provide algorithm developers with a set of challenging problems to evaluate their numerical optimization methods. An important focus is given on '''reproducibility''' of optimal solutions.
As, in contrast to say linear programming, there are no standard formats for the formulation of such problems, and they often show completely different characteristics, these pages dedicate some space for a thorough description of problem and solutions.

A more detailed description of the underlying concepts of this library can be found in the article <bib id="Sager2012b" />, of which a [http://mathopt.de/PUBLICATIONS/Sager2012b.pdf preprint pdf] is available.


==[[:Category:News|News]] [[Current News|(add)]]==
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px">
<tr valign="top">
<td width="80%" bgcolor="#EEEEFF">

{{Current News}}
</td>
</tr>
</table>

<table width="100%">


<tr valign="top">
<td width="50%">
==[[:Category:Problem characterization|Problem characterization]] [[Help:Adding a problem characterization|(add)]]==

*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Population dynamics | Population dynamics]] - [[:Category: Chemical engineering | Chemical engineering]]

*via ''[[:Category:Model characterization|mathematical model]]'', e.g., [[:Category: ODE model | ODE]], [[:Category: PDE model | PDE]] or [[:Category: DAE model | DAE model]] - [[:Category: Multistage process | Multistage process]] - [[:Category: State dependent switches | State dependent switches]]

*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]

*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category: Muscod | Muscod]] - [[:Category: Switch | Switch]]
</td>


<td width="50%" >

==Direct links to [[:Category:MIOCP|Problems]] [[Help:Adding a Problem|(add)]]==

[[Bang-bang approximation of a traveling wave | 1D wave]] - [[Bioreactor]] - [[Batch reactor]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]] - [[Car testdrive]] - [[Cushioned Oscillation]] - [[Diels-Alder Reaction Experimental Design | Diels-Alder OED]] - [[Double Tank]] - [[F-8 aircraft]] - [[Fuller's problem]] - [[Gravity Turn Maneuver | Gravity Turn]] - [[Controlled Heating | Heating]]'' - ''[[Lotka Volterra fishing problem | Lotka]] - [[Lotka Experimental Design | Lotka OED]] - [[Van der Pol Oscillator | Oscillator]] - [[Oil Shale Pyrolysis | Pyrolysis]] - [[Goddart's rocket problem | Rocket]] - [[Source Inversion]] - [[Subway ride | Subway]] - ''[[Supermarket refrigeration system | Supermarket]] - [[Truck cruise control | Truck]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

==[[:Category:Help|Help]] [[User:SebastianSager|(contact)]]==

''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]
</td>


<td width="50%" >

==[[:Category:Community|Community]] [[Help:Adding a User|(add)]]==

[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''

</td>
</tr>

</table>

== References ==
<biblist />

[[Category:Main]]

De Pillis chemotherapy model

2016-01-28T09:29:49Z

SebastianSager: Moved C to Muscod

De Pillis chemotherapy model (C)

2016-01-28T09:29:12Z

SebastianSager: SebastianSager moved page De Pillis chemotherapy model (C) to De Pillis chemotherapy model (Muscod): Moved C to Muscod

#REDIRECT [[De Pillis chemotherapy model (Muscod)]]

De Pillis chemotherapy model (Muscod)

2016-01-28T09:29:12Z

SebastianSager: SebastianSager moved page De Pillis chemotherapy model (C) to De Pillis chemotherapy model (Muscod): Moved C to Muscod

This is source code written in [[:Category:Muscod | Muscod]] for the [[De Pillis chemotherapy model]].

<source lang="cpp">
/* volume of tumor in absolute cell count*/
rhs[0] = - a*x[0] (1-b*x[0]) -c*x[1]*x[0] - D*x[0] - K_T*(1- exp(- x[4]))*x[0];

/* volume of unspecific immune cells */
rhs[1] = e*x[3] - f*x[1] + g*x[0]*x[0]/(h+x[0]*x[0])-p*x[1]*x[0] - K_N*(1- exp(- x[4]))*x[1];

/* volume of tumor-specific cytotoxic T-cells */
rhs[2] = -m*x[2] + j*D*D*x[0]*x[0]/(k+ D*D*x[0]*x[0])*x[2] - q*x[1]*x[2] + (r_1*x[1] + r_2*x[3])*x[0] ...
- v*x[1]* x[2]*x[2] - K_L* (1- exp(- x[4]))*x[2] + p_I*x[2]*x[5]/(g_I + x[5]) + u[2];

/* volume of circulating lymphocytes */
rhs[3] = alpha - beta *x[3] - K_C *(1- exp(- x[4])) *x[3];

/* amount of chemotherapeutic drug */
rhs[4] = u[0] - gamma*x[4];

/* amount of immunotherapeutic drug */
rhs[5] = u[1] - mu_I*x[5];

</source>

[[Category:Muscod]]

De Pillis chemotherapy model (Muscod)

2016-01-28T09:29:01Z

SebastianSager: Moved C to Muscod

D'Onofrio chemotherapy model

2016-01-28T09:27:13Z

SebastianSager: Moved C to Muscod

Lotka Volterra fishing problem (TACO)

2016-01-28T09:25:35Z

SebastianSager:

This page contains a model of the MIOCP [[Lotka Volterra fishing problem]] in [http://www.ampl.org AMPL] format, making use of the TACO toolkit for AMPL control optimization extensions. This problem is due to <bib id="Sager2005" />.
Note that you will need to include a generic [[support AMPL files|AMPL/TACO support file]], OptimalControl.mod.
To solve this model, you require an optimal control or MINLP code that uses the TACO toolkit to support the AMPL optimal control extensions.

=== AMPL ===

This is the source file lotka_taco.mod
<source lang="AMPL">
# ----------------------------------------------------------------
# Lotka Volterra fishing problem using AMPL and TACO
# (c) Christian Kirches, Sven Leyffer
# ----------------------------------------------------------------
include OptimalControl.mod;

var t;
var xd0 := 0.5, >= 0, <= 20;
var xd1 := 0.7, >= 0, <= 20;
var dev := 0.0, >= 0, <= 20;

var u := 1, >= 0, <= 1 integer suffix type "u0";

param p0 := 0.4;
param p1 := 0.2;
param ref0 := 1.0;
param ref1 := 1.0;

# Minimize accumulated deviation from reference after 12 time units
minimize Mayer: eval(dev,12);

subject to

# ODE system
ODE_0: diff(xd0,t) = xd0 - xd0*xd1 - p0*u*xd0; # prey
ODE_1: diff(xd1,t) = -xd1 + xd0*xd1 - p1*u*xd1; # predator
ODE_2: diff(dev,t) = (xd0-ref0)^2 + (xd1-ref1)^2; # deviation from reference

# initial value constraints
IVC_0: eval(xd0,0) = 0.5;
IVC_1: eval(xd1,0) = 0.7;
IVC_2: eval(dev,0) = 0.0;

option solver ...;

solve;

</source>

== References ==
<biblist />

[[Category:AMPL/TACO]]

Supermarket refrigeration system

2016-01-28T09:24:31Z

SebastianSager: Moved C to Muscod

{{Dimensions
|nd = 1
|nx = 9
|nw = 3
|nre = 9
}}

The '''supermarket refrigeration system problem''' is based on a model describing a refrigeration system with 2 parallel connected compressors (called a compressor rack) which only can be controlled stepwise (each single compressor can be turned on or off) and 2 open refrigerated display cases containing goods needed to be refrigerated. Each display case is connected to the refrigeration circuit through an expansion valve which also can only be closed or opened. This valve controls the flow of refrigerant into the evaporator, where it absorbs heat from the surrounding air. The refrigerated air then creates the well-known air-curtain at the front of the display case.

The air temperatures surrounding the goods in each display case are modeled by one differential state each. These states have to be bounded, so that the goods are properly refrigerated.

The model was published by Larsen et. al. in 2007 <bib id="Larsen2007" />. The main goal is to control the refirgeration system energy-optimal. The problem was set up as a benchmark problem for MIOCPs.

The mathematical equations form an [[:Category:ODE model|ODE model]]. The initial values of the differential states are not fixed but periodicity of the whole process is required.

The optimal integer control function shows [[:Category:Chattering|chattering]] behavior, making the supermarket refrigeration system problem a candidate for benchmarking of algorithms.

== Mathematical formulation ==

For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by


<math>
\min_{x,u} \frac {1}{t_f - t_0}\int_{t_0}^{t_f} \left( u_2 \cdot 0.5 \cdot \eta_{vol} \cdot V_{sl} \cdot f \right) dt
</math>



<math>
\begin{array}{llcl}
\displaystyle
\mbox{s.t.} &
\dot{x_0}(t) &=& \dfrac{1}{V_{suc} \cdot \frac{d\rho_{suc}}{dP_{suc}}(x_0)} \cdot \bigg[
\left(\dfrac{UA_{wall-ref, max}}{M_{ref, max} \cdot
\Delta h_{lg}(x_0)}\right) \Big( x_4 \big( x_2 - T_e(x_0) \big)\\
& && + \, x_8 \big( x_6 - T_e(x_0) \big) \Big) + \, M_{ref,const} - \eta_{vol} \cdot V_{sl} \cdot 0.5 \, u_2 \rho_{suc}(x_0)
\bigg] \\
&\dot{x_1}(t) &=& - \dfrac{
UA_{goods-air} \left( x_1 - x_3 \right)
}{
M_{goods} \cdot C_{p,goods}
} \\
&\dot{x_2}(t) &=& \dfrac{
UA_{air-wall} \left( x_3-x_2 \right)
- \dfrac{UA_{wall-ref,max}}{M_{ref,max}}
\, x_4 \big( x_2 - T_e(x_0) \big)
}{
M_{wall} \cdot C_{p,wall}
} \\ [2.5ex]
&\dot{x_3}(t) &=& \dfrac{
UA_{goods-air} \left( x_1-x_3 \right) + \dot{Q}_{airload}
- UA_{air-wall} \, (x_3-x_2)
}{
M_{air} \cdot C_{p,air}
} \\ [2.5ex]
&\dot{x_4}(t) &=& \left(\dfrac{M_{ref,max} - x_4}{\tau_{fill}} \right) u_0
- \left( \dfrac{UA_{wall-ref,max}}{M_{ref,max} \cdot \Delta h_{lg}(x_0)} \,
x_4 \big( x_2 - T_e(x_0) \big) \right) (1-u_0)
\\ \\
&\dot{x_5}(t) &=& - \dfrac{
UA_{goods-air} \left( x_5 - x_7 \right)
}{
M_{goods} \cdot C_{p,goods}
} \\
&\dot{x_6}(t) &=& \dfrac{
UA_{air-wall} \left( x_7-x_6 \right)
- \dfrac{UA_{wall-ref,max}}{M_{ref,max}}
\, x_8 \big( x_6 - T_e(x_0) \big)
}{
M_{wall} \cdot C_{p,wall}
} \\ [2.5ex]
&\dot{x_7}(t) &=& \dfrac{
UA_{goods-air} \left( x_5-x_7 \right) + \dot{Q}_{airload}
- UA_{air-wall} \, (x_7-x_6)
}{
M_{air} \cdot C_{p,air}
} \\ [2.5ex]
&\dot{x_8}(t) &=& \left(\dfrac{M_{ref,max} - x_8}{\tau_{fill}} \right) u_1
- \left( \dfrac{UA_{wall-ref,max}}{M_{ref,max} \cdot \Delta h_{lg}(x_0)} \,
x_8 \big( x_6 - T_e(x_0) \big) \right) (1-u_1)
\\ [4ex]

& x_3(t) &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\
& x_3(t) &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\
& x_7(t) &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\
& x_7(t) &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\
& x_0(t) &\leq& 1.7 \quad \forall t \in [t_0, t_f], \\
& x_i(t_0) &=& free \quad \forall i \in \{0,\dots, 8\}, \\
& x_i(t_f) &=& x_i(t_0) \quad \forall i \in \{0,\dots, 8\}, \\
& u_i(t) &\in& \{0, 1\} \quad \forall i \in \{0,\dots, 2\}, \\
& t_f &\in& [ 650, 750 ].

\end{array}
</math>


Here the differential state <math>x_0</math> describes the suction pressure in the suction manifold (in bar). The next three states model temperatures in the first display case (in °C). <math>x_1</math> is the goods' temperature, <math>x_2</math> the one of the evaporator wall and <math>x_3</math> the air temperature surrounding the goods. <math>x_4</math> then models the mass of the liquefied refrigerant in the evaporator (in kg).

<math>(x_5, x_6, x_7, x_8)</math> describe the corresponding states in the second display case.

<math>u_0</math> describes the inlet valve of the first display case, <math>u_1</math> respectively the valve of the second display case.
<math>u_2, u_3</math> denote the activity of a single compressor.

The following polynomial functions are used in the model description originating from interpolations:


<math>
\begin{array}{rcl}

T_e(x_0) &=& -4.3544 x_0^2 + 29.224 x_0 - 51.2005,\\
\Delta h_{lg}(x_0) &=& (0.0217 x_0^2 - 0.1704 x_0 + 2.2988)\cdot 10^5, \\
\rho_{suc}(x_0) &=& 4.6073 x_0 + 0.3798, \\
\frac{d\rho_{suc}}{dP_{suc}}(x_0) &=& -0.0329 {x_0}^3 + 0.2161 {x_0}^2 - 0.4742 x_0 + 5.4817,\\
f(x_0) &=& (0.0265 x_0^3 - 0.4346 x_0^2 + 2.4923 x_0 + 1.2189) \cdot 10^5.
\end{array}
</math>


== Parameters ==

These fixed values are used within the model for the day scenario. A night scenario is also available, see [[#Variants|Variants]].

{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Symbol !! Value !! Unit !! Description
|-
| align=center | <math>\dot{Q}_{airload}</math> || align=right | 3000.00 || <math>\frac{J}{s}</math> || Disturbance, heat transfer from outside the display case
|-
| align=center | <math>\dot{m}_{ref,const}</math> || align=right | 0.20 || <math>\frac{kg}{s}</math> || Disturbance, constant mass flow of refrigerant
from unmodeled entities
|-
| align=center | <math>M_{goods}</math> || align=right | 200.00 || <math>kg</math> || Mass of goods
|-
| align=center | <math>C_{p,goods}</math> || align=right | 1000.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of goods
|-
| align=center | <math>
UA_{goods-air} </math> || align=right | 300.00 || <math>\frac{J}{s \cdot K}</math> || Heat transfer coefficient between goods
and air
|-
| align=center | <math>M_{wall} </math> || align=right | 260.00 || <math>kg</math> || Mass of evaporator wall
|-
| align=center | <math>C_{p,wall} </math> || align=right | 385.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of evaporator wall
|-
| align=center | <math>UA_{air-wall}</math> || align=right | 500.00 || <math>\frac{J}{s \cdot K}</math> || Heat transfer coefficient between air and
wall
|-
| align=center | <math>M_{air}</math> || align=right | 50.00 || <math>kg</math> || Mass of air in display case
|-
| align=center | <math>C_{p,air}</math> || align=right | 1000.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of air
|-
| align=center | <math>UA_{wall-ref,max}</math> || align=right | 4000.00 || <math>\frac{J}{s \cdot K}</math> || Maximum heat transfer coefficient between
refrigerant and evaporator wall
|-
| align=center | <math>\tau_{fill}</math> || align=right | 40.00 || <math>s</math> || Parameter describing the filling time of the
evaporator under opened valve
|-
| align=center | <math>T_{SH}</math> || align=right | 10.00 || <math>K</math> || Superheat in the suction manifold
|-
| align=center | <math>M_{ref,max}</math> || align=right | 1.00 || <math>kg</math> || Maximum mass of refrigerant in evaporator
|-
| align=center | <math>V_{suc}</math> || align=right | 5.00 || <math>m^3</math> || Total volume of suction manifold
|-
| align=center | <math>V_{sl}</math> || align=right | 0.08 || <math>\frac{m^3}{s}</math> || Total displacement volume
|-
| align=center | <math>\eta_{vol}</math> || align=right | 0.81 || <math>-</math> || Volumetric efficiency
|}

== Reference Solutions ==

For the relaxed problem (we only demand <math>u_i(t) \in [0,1]</math> instead of <math>u_i(t) \in \{0,1\}</math> the optimal solution is 12072.45.
The illustrated solution with integer controls has a (suboptimal) objective function value of 12252.81.

<gallery caption="Reference solution plots" widths="350px" heights="300px" perrow="2">
Image:FridgeControlsRelaxed.png| Optimal relaxed control.
Image:FridgeControls.png| Optimal integer control.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:Muscod | Muscod code]] at [[Supermarket refrigeration system (Muscod)]]
* [[:Category:optimica | optimica]] at [[Supermarket refrigeration system (optimica)]]

== Variants ==

Since the compressors are parallel connected one can introduce a single control <math> u_2 \in \{0,1,2\}</math> instead of two equivalent controls. The same holds for scenarions with <math> n </math> parallel connected compressors.

In the paper <bib id="Larsen2007" /> mentioned above, the problem was stated slightly different:

* The temperature constraints weren't hard bounds but there was a penalization term added to the objective function to minimize the violation of these constraints.
* The differential equation for the mass of the refrigerant had another switch, if the valve (e.g. <math>u_0</math>) is closed. It was formulated this way:
<math>
\dot{x_4} = \begin{cases}

\dfrac{M_{ref,max} - x_4}{\tau_{fill}} & \qquad \text{if} \quad u_0 = 1 \\ \\
- \dfrac{UA_{wall-ref,max}}{M_{ref,max} \cdot \Delta h_{lg}(x_0)} x_4 \big( x_2 - T_e(x_0) \big) & \qquad \text{if} \quad u_0 = 0 \quad \text{and}\quad x_4 > 0 \\ \\
0 & \qquad \text{if} \quad u_0 = 0 \quad \text{and} \quad x_4 = 0

\end{cases}
</math>

This additional switch is redundant because the mass itself is a factor on the right hand side and so the complete right hand side is 0 if <math>x_4 = 0</math>.
* A night scenario with two different parameters was given. At night the following parameters change their value:
<math>
\begin{array}{lcrr}
\dot{Q}_{airload} &=& 1800.00 & \frac{J}{s}, \\
\dot{m}_{ref,const} &=& 0.00 & \frac{kg}{s}, \\
\end{array}
</math>

Additionally the constraint on the suction pressure <math>x_0(t)</math> is softened to <math>x_0(t) \leq 1.9</math>.
* No periodicity was required but the solution on a fixed time horizon 4 hours - 2 in day scenario and 2 in night scenario - with <math>t_f = 14400</math> was asked.
* The number of compressors and display cases is not fixed. Larsen also proposed the problem with 3 compressors and 3 display cases. This leads to a change in the compressor rack's preformance to <math>V_{sl} = 0.095 \frac{m^3}{s}</math>. Unfortunately this constant is only given for these two cases although Larsen proposed scenarios with more compressors and display cases.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chattering]]
[[Category:Periodic]]
[[Category:Minimum energy]]

Supermarket refrigeration system (Muscod)

2016-01-28T09:24:03Z

SebastianSager: SebastianSager moved page Supermarket refrigeration system (C) to Supermarket refrigeration system (Muscod) without leaving a redirect: Moved C to Muscod

The differential equations for the [[Supermarket refrigeration system]] in [[:Category:Muscod | Muscod code]] read as follows

<source lang="cpp">
// number of compressors
#define NVALVES 2

// constants
#define M_GOODS 200.0
#define C_P_GOODS 1000.0
#define UA_GOODS_AIR 300.0
#define M_WALL 260.0
#define C_P_WALL 385.0
#define UA_AIR_WALL 500.0
#define M_AIR 50.0
#define C_P_AIR 1000.0
#define UA_WALL_REF_MAX 4000.0
#define M_REF_MAX 1.0
#define TAU_FILL 40.0
#define T_SH 10.0
#define V_SUC 5.0
#define V_SL 0.08 // 2 display cases - 2 compressors
// #define V_SL 0.095 // 3 display cases - 3 compressors
#define ETA_VOL 0.81

// disturbances - day scenario
#define Q_AIRLOAD 3000.0
#define M_REF_CONST 0.2

// disturbances - night scenario
// #define Q_AIRLOAD 1800.0
// #define M_REF_CONST 0.0

double delta_h = (0.0217*xd[0]*xd[0] - 0.1704*xd[0] + 2.2988)*1e5;
double T_e = -4.3544*xd[0]*xd[0] + 29.224*xd[0] - 51.2005;
double rho_suc = 4.6073*xd[0] + 0.3798;
double rho_suc__P_suc = -0.0329*xd[0]*xd[0]*xd[0] + 0.2161*xd[0]*xd[0] - 0.4742*xd[0] + 5.4817;
double f = (0.0265*xd[0]*xd[0]*xd[0] -0.4346*xd[0]*xd[0] + 2.4923*xd[0] + 1.2189)*1e5;

double Q_goods_air[NVALVES];
double Q_air_wall[NVALVES];
double UA_wall_ref[NVALVES];
double Q_e[NVALVES];
double m[NVALVES];

double m_in_suc = 0.0;

int i;

for (i=0; i<NVALVES; i++){
Q_goods_air[i] = UA_GOODS_AIR*(xd[1 + i*4] - xd[3 + i*4]);
Q_air_wall[i] = UA_AIR_WALL*(xd[3 + i*4] - xd[2 + i*4]);
UA_wall_ref[i] = UA_WALL_REF_MAX * xd[4 + 4*i]/M_REF_MAX;
Q_e[i] = UA_wall_ref[i]*(xd[2 + 4*i] - T_e);

m[i] = Q_e[i]/delta_h;
m_in_suc += m[i];
}

double V_comp = 0.0;
double comp_scale = (double) 1.0/NCOMPS;
V_comp = comp_scale*u[NVALVES]*ETA_VOL*V_SL;

// suction pressure
rhs[0] = (m_in_suc + M_REF_CONST - V_comp*rho_suc) / (V_SUC*rho_suc__P_suc);

// for each display/valve
for (i=0; i<NVALVES; i++){

// temperatures:

// goods
rhs[1 + i*4] = - Q_goods_air[i]/(M_GOODS*C_P_GOODS);

// wall
rhs[2 + i*4] = (Q_air_wall[i] - Q_e[i])/(M_WALL*C_P_WALL);

// air
rhs[3 + i*4] = ((Q_goods_air[i] + Q_AIRLOAD - Q_air_wall[i]) /(M_AIR*C_P_AIR));

// mass of liquefied refrigerant:
rhs[4 + i*4] = ((M_REF_MAX - xd[4 + 4*i])/TAU_FILL * u[i] - m[i] * (1 - u[i]));
}
</source>

[[Category:Muscod]]

Supermarket refrigeration system (Muscod)

2016-01-28T09:23:53Z

SebastianSager: Moved C to Muscod