https://mintoc.de/api.php?action=feedcontributions&user=SebastianSager&feedformat=atommintOC - User contributions [en]2024-03-28T14:36:21ZUser contributionsMediaWiki 1.25.2https://mintoc.de/index.php?title=LinearMetabolic&diff=2366LinearMetabolic2023-10-20T11:35:17Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nu = 3<br />
|nc = 1<br />
|nri = 0<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
The gIOC is given by<br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|<br />
</math><br />
subject to<br />
<math><br />
(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 <br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& \displaystyle 0 \\<br />
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\<br />
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\<br />
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\<br />
x(0) &=& (1, 0, 0, 0) \\<br />
x_4(T) &=& 0.9 \\<br />
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\ <br />
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\<br />
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\<br />
w_1, w_2 \in [0,1] \\<br />
w_1 + w_2 = 1<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Data ==<br />
<br />
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.<br />
<br />
{| class="wikitable"<br />
|+Data<br />
|-<br />
|weights <math>w_1, w_2</math><br />
|parameters <math>p_1, p_2, p_3</math><br />
|noise<br />
|data <math>\eta</math><br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
| 0<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
|0.2<br />
|<br />
[[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],<br />
[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],<br />
[ 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],<br />
[ 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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0<br />
|<br />
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0.1<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|}<br />
<br />
== Examples ==<br />
The plots show noise free data for different weights from forward optimal control problems.<br />
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3"><br />
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].<br />
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].<br />
</gallery><br />
<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2365LinearMetabolic2023-10-20T11:28:15Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 3<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
The gIOC is given by<br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|<br />
</math><br />
subject to<br />
<math><br />
(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 <br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& \displaystyle 0 \\<br />
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\<br />
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\<br />
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\<br />
x(0) &=& (1, 0, 0, 0) \\<br />
x_4(T) &=& 0.9 \\<br />
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\ <br />
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\<br />
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\<br />
w_1, w_2 \in [0,1] \\<br />
w_1 + w_2 = 1<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Data ==<br />
<br />
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.<br />
<br />
{| class="wikitable"<br />
|+Data<br />
|-<br />
|weights <math>w_1, w_2</math><br />
|parameters <math>p_1, p_2, p_3</math><br />
|noise<br />
|data <math>\eta</math><br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
| 0<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
|0.2<br />
|<br />
[[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],<br />
[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],<br />
[ 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],<br />
[ 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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0<br />
|<br />
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0.1<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|}<br />
<br />
== Examples ==<br />
The plots show noise free data for different weights from forward optimal control problems.<br />
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3"><br />
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].<br />
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].<br />
</gallery><br />
<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2364LinearMetabolic2023-10-20T11:27:04Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 3<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
The gIOC is given by<br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|<br />
</math><br />
subject to<br />
<math><br />
(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 <br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& \displaystyle 0 \\<br />
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\<br />
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\<br />
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\<br />
x(0) &=& (1, 0, 0, 0) \\<br />
x_4(T) &=& 0.9 \\<br />
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\ <br />
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\<br />
1 &\ge& u_1(t) + u_2(t) + u_3(t) \\<br />
w_1, w_2 \in [0,1] \\<br />
w_1 + w_2 = 1<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Data ==<br />
<br />
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.<br />
<br />
{| class="wikitable"<br />
|+Data<br />
|-<br />
|weights <math>w_1, w_2</math><br />
|parameters <math>p_1, p_2, p_3</math><br />
|noise<br />
|data <math>\eta</math><br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
| 0<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
|0.2<br />
|<br />
[[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],<br />
[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],<br />
[ 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],<br />
[ 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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0<br />
|<br />
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0.1<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|}<br />
<br />
== Examples ==<br />
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.<br />
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3"><br />
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].<br />
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].<br />
</gallery><br />
<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2363LinearMetabolic2023-10-20T11:26:21Z<p>SebastianSager: Edited the data table</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 3<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
The gIOC is given by<br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|<br />
</math><br />
subject to<br />
<math><br />
(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 <br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& \displaystyle 0 \\<br />
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\<br />
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\<br />
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\<br />
x(0) &=& (1, 0, 0, 0) \\<br />
x_4(T) &=& 0.9 \\<br />
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\ <br />
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\<br />
1 &\ge& u_1 + u_2 + u_3 \\<br />
w_1, w_2 \in [0,1] \\<br />
w_1 + w_2 = 1<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Data ==<br />
<br />
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.<br />
<br />
{| class="wikitable"<br />
|+Data<br />
|-<br />
|weights <math>w_1, w_2</math><br />
|parameters <math>p_1, p_2, p_3</math><br />
|noise<br />
|data <math>\eta</math><br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
| 0<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
|0.2<br />
|<br />
[[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],<br />
[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],<br />
[ 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],<br />
[ 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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0<br />
|<br />
[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0.1<br />
|<br />
[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
|-<br />
|}<br />
<br />
== Examples ==<br />
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.<br />
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3"><br />
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].<br />
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].<br />
</gallery><br />
<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2361LinearMetabolic2023-10-20T06:46:23Z<p>SebastianSager: Removed old code and p</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 3<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
Under construction...<br />
<br />
The gIOC is given by<br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*, T^*)} \; \| x^* - \eta \|<br />
</math><br />
subject to<br />
<math><br />
(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 <br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& \displaystyle 0 \\<br />
\dot{x_2}(t) &=& \displaystyle p_1 \cdot x_1 \cdot u_1 - p_2 \cdot x_2 \cdot u_2 \\<br />
\dot{x_3}(t) &=& \displaystyle p_2 \cdot x_2 \cdot u_2 - p_3 \cdot x_3 \cdot u_3 \\<br />
\dot{x_4}(t) &=& \displaystyle p_3 \cdot x_3 \cdot u_3 \\<br />
x(0) &=& (1, 0, 0, 0) \\<br />
x_4(T) &=& 0.9 \\<br />
x_i(t) &\ge& 0 \; \forall \; i \in \{1,2,3,4\} \\ <br />
u_i(t) &\in& [0,1] \; \forall \; i \in \{1,2,3\} \\<br />
w_1, w_2 \in [0,1] \\<br />
w_1 + w_2 = 1<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Data ==<br />
{| class="wikitable"<br />
|+Data<br />
|-<br />
|weights<br />
|parameters<br />
|noise<br />
|data<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
| 0<br />
|[[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],<br />
[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],<br />
[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],<br />
[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]]<br />
<br />
|-<br />
|[0.232, 0.768]<br />
|[1,1,1]<br />
|0.2<br />
|[[1. ,0.9590812 ,0.82266497,0.95877308,0.60614596,1.20781224,<br />
0.92205471,1.0215116 ,0.89080449,1.05747735,1.04019676,1.06665347,<br />
0.78545416,1.17084962,1.08658091,1.01467309,0.91098074,1.10349004,<br />
0.83411596,0.73187308,0.70399766],[0. ,0.37493181,0.10514573,0.79448823,0.44666767,1.19676033,<br />
1.14532119,1.64453389,1.93572019,1.37435635,0.92790967,1.11928728,<br />
0.52061388,0.43837628,0.31518386,0.50300858,0.52891948,0.74255261,<br />
0.07589526,0.64051491,0.47184952],[ 0. , 0.23680072, 0.05047157,-0.03770305, 0.34603556,-0.45259701,<br />
0.25340369, 0.25177705, 0.16210031, 0.34367395, 0.96197961, 1.10224845,<br />
1.21656021, 1.38442255, 1.02149302, 0.68818982, 0.7921558 , 0.32238544,<br />
0.30229594, 0.43876605, 0.22073986],[ 0. ,-0.01591557,-0.12432336,-0.17669569,-0.21325853,-0.2489822 ,<br />
-0.22989873,-0.0564741 ,-0.19961696,-0.16036058,-0.04386598, 0.25004781,<br />
-0.07279584, 0.09573081,-0.12329308,-0.03451486, 0.0173815 , 0.79895333,<br />
0.8178261 , 0.53498205, 1.1021996 ]]<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0<br />
|[[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],<br />
[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],<br />
[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],<br />
[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]]<br />
<br />
|-<br />
|[0,1]<br />
|[1,1,1]<br />
|0.1<br />
|[[1. ,0.98392517,1.02825092,1.0127091 ,1.21814838,1.03987996,<br />
1.17367563,1.0721497 ,1.11254303,0.91560159,0.92664064,0.93069442,<br />
1.01285111,1.10017347,0.9952073 ,1.05616016,1.12564585,1.02743444,<br />
0.93306189,1.04808763,0.87659083],[0. ,0.19112834,0.6016258 ,0.70894548,0.86077117,0.95316772,<br />
1.12244906,1.51535841,1.48548905,1.87282215,2.13181979,1.48523121,<br />
1.14624017,1.06839771,0.8838448 ,0.72466348,0.61379013,0.70096009,<br />
0.73549769,0.62452413,0.61907767],[ 0. , 0.04514029,-0.06510522,-0.00304812,-0.01620013, 0.006101 ,<br />
0.04992798, 0.29317191, 0.11305888, 0.2285637 , 0.26142926, 0.50457228,<br />
0.87829835, 0.90900075, 1.32279588, 1.42225272, 1.18086941, 0.67720058,<br />
0.61562944, 0.6339817 , 0.36613472],[ 0. ,-0.13882966, 0.18200648,-0.06629627,-0.01496019, 0.21349797,<br />
-0.10124824,-0.05505019, 0.23235014, 0.09448396, 0.06059911,-0.12137716,<br />
-0.00781841, 0.04355511, 0.1369212 ,-0.13011339, 0.2027646 , 0.51512985,<br />
0.65331917, 0.78324169, 0.81304223]]<br />
<br />
|-<br />
|}<br />
<br />
== Examples ==<br />
Here are two plots with noise free data for weights [0.232, 0.768] to compare with the solution of different optimal control problems.<br />
<gallery caption="Solution plots" widths="280px" heights="240px" perrow="3"><br />
Image:LPN3B_min_sum_data.png| OCP with weights [0.232, 0.768].<br />
Image:LPN3B_minT_x_data.png| OCP with weights [0,1].<br />
</gallery><br />
<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Main_Page&diff=2352Main Page2023-10-19T13:21:22Z<p>SebastianSager: </p>
<hr />
<div>__NOTOC__<br />
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.<br />
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.<br />
<br />
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.<br />
<br />
<!-- TOP TABLE WITH NEWS --><br />
==[[:Category:News|News]] <span style="font-variant:small-caps" style="font-size:10px">[[Current News|(add)]]</span>==<br />
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px"><br />
<tr valign="top"><br />
<td width="80%" bgcolor="#EEEEFF"><br />
<!-- The actual news are being included from the page Current News with the following line. --><br />
{{Current News}}<br />
</td><br />
</tr><br />
</table><br />
<br />
<table width="100%"><br />
<br />
<!-- CATEGORY FIELD--><br />
<tr valign="top"><br />
<td width="50%"><br />
==[[:Category:Problem characterization|Problem characterization]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a problem characterization|(add)]]</span>==<br />
<br />
*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]]<br />
<br />
*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]]<br />
<br />
*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]<br />
<br />
*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]<br />
</td><br />
<br />
<!-- PROBLEM FIELD--><br />
<td width="50%" ><br />
<br />
==Direct links to [[:Category:MIOCP|Problems]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a Problem|(add)]]</span>==<br />
<br />
[[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]]''<br />
</td><br />
</tr><br />
<br />
<tr valign="top"><br />
<br />
<!-- HELP FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Help|Help]] <span style="font-variant:small-caps" style="font-size:10px">[[User:SebastianSager|(contact)]]</span>==<br />
<br />
''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]<br />
</td><br />
<br />
<!-- COMMUNITY FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Community|Community]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a User|(add)]]</span>==<br />
<br />
[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''<br />
<br />
</td><br />
</tr><br />
<br />
</table><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Main]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Main_Page&diff=2351Main Page2023-10-19T13:19:38Z<p>SebastianSager: added GIOC direct link</p>
<hr />
<div>__NOTOC__<br />
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.<br />
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.<br />
<br />
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.<br />
<br />
<!-- TOP TABLE WITH NEWS --><br />
==[[:Category:News|News]] <span style="font-variant:small-caps" style="font-size:10px">[[Current News|(add)]]</span>==<br />
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px"><br />
<tr valign="top"><br />
<td width="80%" bgcolor="#EEEEFF"><br />
<!-- The actual news are being included from the page Current News with the following line. --><br />
{{Current News}}<br />
</td><br />
</tr><br />
</table><br />
<br />
<table width="100%"><br />
<br />
<!-- CATEGORY FIELD--><br />
<tr valign="top"><br />
<td width="50%"><br />
==[[:Category:Problem characterization|Problem characterization]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a problem characterization|(add)]]</span>==<br />
<br />
*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]<br />
<br />
*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]]<br />
<br />
*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]<br />
<br />
*via ''[[:Category: Implementation | implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]<br />
</td><br />
<br />
<!-- PROBLEM FIELD--><br />
<td width="50%" ><br />
<br />
==Direct links to [[:Category:MIOCP|Problems]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a Problem|(add)]]</span>==<br />
<br />
[[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]]''<br />
</td><br />
</tr><br />
<br />
<tr valign="top"><br />
<br />
<!-- HELP FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Help|Help]] <span style="font-variant:small-caps" style="font-size:10px">[[User:SebastianSager|(contact)]]</span>==<br />
<br />
''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]<br />
</td><br />
<br />
<!-- COMMUNITY FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Community|Community]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a User|(add)]]</span>==<br />
<br />
[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''<br />
<br />
</td><br />
</tr><br />
<br />
</table><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Main]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Template:Current_News&diff=2350Template:Current News2023-10-19T13:18:17Z<p>SebastianSager: </p>
<hr />
<div>{| style="background:#EEEEFF;"<br />
<!-- Add news below this comment --><br />
{{News|2023/10/19| Added new category [[:Category:GIOC | generalized inverse optimal control]]}}<br />
{{News|2023/10/19| Added [[LinearMetabolic]] problem}}<br />
<!-- Add news above this comment --><br />
|}<br />
<!-- Text within the following tag is not shown when this page is included in a different page (e.g. the main page, "News" section --><br />
<noinclude> <br />
== Older news ==<br />
<br />
{| style="background:#eaecd0;"<br />
<!-- Copy old news below this comment --><br />
{{News|2019/3/14| Added multiple control problems with [[:Category: Gekko | Python GEKKO]] implementation}}<br />
{{News|2018/9/12| Added [[Control of Transmission Lines]] problem}}<br />
{{News|2016/7/26 | Added [[Industrial robot]] problem}}<br />
{{News|2016/7/24 | Added [[Continuously Stirred Tank Reactor problem | CSTR problem]]}}<br />
{{News|2016/6/30 | Added [[Electric Car]] problem}}<br />
{{News|2016/5/5| Added multiple control problems with [[:Category: AMPL/TACO | AMPL with TACO]] implementation}}<br />
{{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}}<br />
{{News|2016/1/10|Added [[Van der Pol Oscillator]], [[Batch reactor]], [[Cushioned Oscillation]] and [[Goddart's rocket problem]] control problems}}<br />
{{News|2015/11/10|Moved mintoc.de to a new server}}<br />
{{News|2012/09/01|Added the [[Lotka Experimental Design]] and a [[:Category:Optimum Experimental Design | category for experimental design problems]]}}<br />
{{News|2011/09/29|Added the [[:Category:AMPL/TACO | first set of AMPL optimal control problems]] using the TACO toolkit}}<br />
{{News|2010/11/21|Added New York [[Subway ride]] control problem}}<br />
{{News|2010/11/18|Extended description of [[:Category:Problem characterization | problem characterization]]}}<br />
{{News|2010/11/18|Description of benchmark library as [http://mathopt.de/PUBLICATIONS/Sager2011b.pdf pdf file] preprint}}<br />
{{News|2010/08/16|Added [[Bang-bang approximation of a traveling wave]] 1D PDE example}}<br />
{{News|2010/02/11|Launch of EU project [http://embocon.org embocon.org]}}<br />
{{News|2009/11/20|Added [[External Links]] page}}<br />
{{News|2009/10/26|Revision and Correction of several [[:Category:Optimica | optimica]] models}}<br />
{{News|2009/08/26|Knitro solution for [[F-8 aircraft (AMPL)]]}}<br />
{{News|2009/08/11|[[F-8 aircraft]] new local minimum found with [[:Category:Optimica | optimica]]/ipopt}}<br />
{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}<br />
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}<br />
<!-- Add news above this comment --><br />
|}<br />
<br />
== Adding items ==<br />
If you want to add events you have to edit this page. At the top of the page you will see something like this:<br />
<source lang="text"><br />
{| style="background:#EEEEFF;"<br />
<!-- Add items below this comment --><br />
{{News|2008/08/14, 15:15|Presentation of the restructured group wiki in the group meeting, room 432.}}<br />
<!-- Add items above this comment --><br />
|}<br />
</source><br />
<br />
To add a new item you have to enter a line of the form:<br />
<source lang="text"><br />
{{News|date|content}}<br />
</source><br />
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.<br />
<br />
''Please keep the items in a chronological order and remove older ones to keep the list short.''<br />
<br />
[[Category:News]]<br />
<br />
</noinclude></div>SebastianSagerhttps://mintoc.de/index.php?title=Category:GIOC&diff=2349Category:GIOC2023-10-19T12:53:11Z<p>SebastianSager: integer controls</p>
<hr />
<div>This category includes all generalized inverse optimal control problems.<br />
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> <br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*) \in \Omega_1} \; \| h(x^*, u^*) - \eta \| + R(p,w)<br />
</math><br />
subject to<br />
<math><br />
(x^*,u^*) \in \displaystyle \arg \min_{x, u} \sum_{i \in M_1} w_i \; \phi_i[x,u,p]<br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x}(t) &=& \displaystyle \sum_{i \in M_2} w_i f_i(x(t), u(t), p) \\<br />
0 &\le& \displaystyle w_i \; g_i(x(t), u(t), p) \quad \forall \; i \in M_3\\<br />
(x,u,p,w) &\in& \Omega_2<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
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.<br />
<br />
== Extensions ==<br />
* The functions <math>u \in \mathcal{U}</math> may also include integer controls.<br />
* For some problems the functions may as well depend explicitely on the time <math>t</math>.<br />
* The differential equations might depend on [[:Category:State dependent switches | state-dependent switches]].<br />
* The variables may include [[:Category:Boolean variables | boolean variables]].<br />
* The underlying process might be a [[:Category:Multistage process | multistage process]].<br />
* The dynamics might be [[:Category:Unstable | unstable]].<br />
* There might be an underlying [[:Category:Network topology | network topology]].<br />
* The integer control functions might have been (re)formulated by means of an [[:Category:Outer convexification|outer convexification]].<br />
<br />
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.<br />
<br />
[[Category:Model characterization]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2348LinearMetabolic2023-10-19T12:48:26Z<p>SebastianSager: Introductory text</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 1<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
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.<br />
<br />
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.<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
Under construction...<br />
<br />
The gIOC is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\<br />
& \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)), \\<br />
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\<br />
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]<br />
& x(t_0) &=& (0,0,0,0)^T, \\<br />
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\<br />
& x_0(t) & \in & [-i_{max}, i_{max}], \\<br />
& u(t) &\in& \{-1, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
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.<br />
<br />
== Parameters ==<br />
These fixed values are used within the model.<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|Coefficient of reduction<br />
|<math>K_r</math><br />
|10<br />
|[-]<br />
|-<br />
|Air density<br />
|<math>\rho</math><br />
|1.293<br />
|<math> kg/m^3 </math><br />
|-<br />
|Aerodynamic coefficient<br />
|<math>C_x</math><br />
|0.4<br />
|[-]<br />
|-<br />
|Area in the front of the vehicle<br />
|<math>S</math><br />
|2<br />
|<math> m^2 </math><br />
|-<br />
|Radius of the wheel<br />
|<math>r</math><br />
|0.33<br />
|<math> m </math><br />
|-<br />
|Constant representing the friction of the wheels on the road<br />
|<math>K_f </math><br />
|0.03<br />
|[-]<br />
|-<br />
|Coefficient of the motor torque<br />
|<math>K_m</math><br />
|0.27<br />
|[-]<br />
|-<br />
|Inductor resistance<br />
|<math>R_m</math><br />
|0.03<br />
|<math> Ohms </math><br />
|-<br />
|inductance of the rotor<br />
|<math>L_m</math><br />
|0.05<br />
|[-]<br />
|-<br />
|Mass<br />
|<math>M</math><br />
|250<br />
|<math> kg </math><br />
|-<br />
|Gravity constant<br />
|<math>g</math><br />
|9,81<br />
|[-]<br />
|-<br />
|Battery voltage<br />
|<math>V_{alim}</math><br />
|150<br />
|<math> V </math><br />
|-<br />
|Resistance of the battery<br />
|<math>R_{bat}</math><br />
|0.05<br />
|<math> Ohms </math><br />
|-<br />
|}<br />
<br />
== Reference Solutions ==<br />
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.<br />
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.<br />
<br />
Applying the sum up rounding strategy results in an integer-feasible chattering solution.<br />
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.<br />
<br />
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>.<br />
<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3"><br />
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>.<br />
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.<br />
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>. <br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]<br />
<br />
== Variants ==<br />
<br />
There are several alternative formulations and variants of the above problem, in particular<br />
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math><br />
* 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><br />
* 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>.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2347LinearMetabolic2023-10-19T12:41:07Z<p>SebastianSager: Tsiantis reference</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 1<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
Under construction...<br />
<br />
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 ...<br />
<br />
The problem is discussed in detail in <bib id="Tsiantis2018" />.<br />
<br />
== Mathematical formulation ==<br />
<br />
The mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\<br />
& \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)), \\<br />
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\<br />
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]<br />
& x(t_0) &=& (0,0,0,0)^T, \\<br />
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\<br />
& x_0(t) & \in & [-i_{max}, i_{max}], \\<br />
& u(t) &\in& \{-1, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
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.<br />
<br />
== Parameters ==<br />
These fixed values are used within the model.<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|Coefficient of reduction<br />
|<math>K_r</math><br />
|10<br />
|[-]<br />
|-<br />
|Air density<br />
|<math>\rho</math><br />
|1.293<br />
|<math> kg/m^3 </math><br />
|-<br />
|Aerodynamic coefficient<br />
|<math>C_x</math><br />
|0.4<br />
|[-]<br />
|-<br />
|Area in the front of the vehicle<br />
|<math>S</math><br />
|2<br />
|<math> m^2 </math><br />
|-<br />
|Radius of the wheel<br />
|<math>r</math><br />
|0.33<br />
|<math> m </math><br />
|-<br />
|Constant representing the friction of the wheels on the road<br />
|<math>K_f </math><br />
|0.03<br />
|[-]<br />
|-<br />
|Coefficient of the motor torque<br />
|<math>K_m</math><br />
|0.27<br />
|[-]<br />
|-<br />
|Inductor resistance<br />
|<math>R_m</math><br />
|0.03<br />
|<math> Ohms </math><br />
|-<br />
|inductance of the rotor<br />
|<math>L_m</math><br />
|0.05<br />
|[-]<br />
|-<br />
|Mass<br />
|<math>M</math><br />
|250<br />
|<math> kg </math><br />
|-<br />
|Gravity constant<br />
|<math>g</math><br />
|9,81<br />
|[-]<br />
|-<br />
|Battery voltage<br />
|<math>V_{alim}</math><br />
|150<br />
|<math> V </math><br />
|-<br />
|Resistance of the battery<br />
|<math>R_{bat}</math><br />
|0.05<br />
|<math> Ohms </math><br />
|-<br />
|}<br />
<br />
== Reference Solutions ==<br />
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.<br />
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.<br />
<br />
Applying the sum up rounding strategy results in an integer-feasible chattering solution.<br />
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.<br />
<br />
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>.<br />
<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3"><br />
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>.<br />
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.<br />
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>. <br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]<br />
<br />
== Variants ==<br />
<br />
There are several alternative formulations and variants of the above problem, in particular<br />
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math><br />
* 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><br />
* 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>.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Template:Current_News&diff=2346Template:Current News2023-10-19T12:27:14Z<p>SebastianSager: </p>
<hr />
<div>{| style="background:#EEEEFF;"<br />
<!-- Add news below this comment --><br />
{{News|2023/10/19| Added new category [[:Category:GIOC | generalized inverse optimal control]]}}<br />
{{News|2023/10/19| Added [[LinearMetabolic]] problem}}<br />
{{News|2019/3/14| Added multiple control problems with [[:Category: Gekko | Python GEKKO]] implementation}}<br />
{{News|2018/9/12| Added [[Control of Transmission Lines]] problem}}<br />
<!-- Add news above this comment --><br />
|}<br />
<!-- Text within the following tag is not shown when this page is included in a different page (e.g. the main page, "News" section --><br />
<noinclude> <br />
== Older news ==<br />
<br />
{| style="background:#eaecd0;"<br />
<!-- Copy old news below this comment --><br />
{{News|2016/7/26 | Added [[Industrial robot]] problem}}<br />
{{News|2016/7/24 | Added [[Continuously Stirred Tank Reactor problem | CSTR problem]]}}<br />
{{News|2016/6/30 | Added [[Electric Car]] problem}}<br />
{{News|2016/5/5| Added multiple control problems with [[:Category: AMPL/TACO | AMPL with TACO]] implementation}}<br />
{{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}}<br />
{{News|2016/1/10|Added [[Van der Pol Oscillator]], [[Batch reactor]], [[Cushioned Oscillation]] and [[Goddart's rocket problem]] control problems}}<br />
{{News|2015/11/10|Moved mintoc.de to a new server}}<br />
{{News|2012/09/01|Added the [[Lotka Experimental Design]] and a [[:Category:Optimum Experimental Design | category for experimental design problems]]}}<br />
{{News|2011/09/29|Added the [[:Category:AMPL/TACO | first set of AMPL optimal control problems]] using the TACO toolkit}}<br />
{{News|2010/11/21|Added New York [[Subway ride]] control problem}}<br />
{{News|2010/11/18|Extended description of [[:Category:Problem characterization | problem characterization]]}}<br />
{{News|2010/11/18|Description of benchmark library as [http://mathopt.de/PUBLICATIONS/Sager2011b.pdf pdf file] preprint}}<br />
{{News|2010/08/16|Added [[Bang-bang approximation of a traveling wave]] 1D PDE example}}<br />
{{News|2010/02/11|Launch of EU project [http://embocon.org embocon.org]}}<br />
{{News|2009/11/20|Added [[External Links]] page}}<br />
{{News|2009/10/26|Revision and Correction of several [[:Category:Optimica | optimica]] models}}<br />
{{News|2009/08/26|Knitro solution for [[F-8 aircraft (AMPL)]]}}<br />
{{News|2009/08/11|[[F-8 aircraft]] new local minimum found with [[:Category:Optimica | optimica]]/ipopt}}<br />
{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}<br />
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}<br />
<!-- Add news above this comment --><br />
|}<br />
<br />
== Adding items ==<br />
If you want to add events you have to edit this page. At the top of the page you will see something like this:<br />
<source lang="text"><br />
{| style="background:#EEEEFF;"<br />
<!-- Add items below this comment --><br />
{{News|2008/08/14, 15:15|Presentation of the restructured group wiki in the group meeting, room 432.}}<br />
<!-- Add items above this comment --><br />
|}<br />
</source><br />
<br />
To add a new item you have to enter a line of the form:<br />
<source lang="text"><br />
{{News|date|content}}<br />
</source><br />
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.<br />
<br />
''Please keep the items in a chronological order and remove older ones to keep the list short.''<br />
<br />
[[Category:News]]<br />
<br />
</noinclude></div>SebastianSagerhttps://mintoc.de/index.php?title=Category:GIOC&diff=2345Category:GIOC2023-10-19T12:24:08Z<p>SebastianSager: Initial setup of gIOC class</p>
<hr />
<div>This category includes all generalized inverse optimal control problems.<br />
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> <br />
<br />
<math><br />
\displaystyle \min_{(p, w, x^*, u^*) \in \Omega_1} \; \| h(x^*, u^*) - \eta \| + R(p,w)<br />
</math><br />
subject to<br />
<math><br />
(x^*,u^*) \in \displaystyle \arg \min_{x, u} \sum_{i \in M_1} w_i \; \phi_i[x,u,p]<br />
</math><br />
subject to<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x}(t) &=& \displaystyle \sum_{i \in M_2} w_i f_i(x(t), u(t), p) \\<br />
0 &\le& \displaystyle w_i \; g_i(x(t), u(t), p) \quad \forall \; i \in M_3\\<br />
(x,u,p,w) &\in& \Omega_2<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
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.<br />
<br />
== Extensions ==<br />
* For some problems the functions may as well depend explicitely on the time <math>t</math>.<br />
* The differential equations might depend on [[:Category:State dependent switches | state-dependent switches]].<br />
* The variables may include [[:Category:Boolean variables | boolean variables]].<br />
* The underlying process might be a [[:Category:Multistage process | multistage process]].<br />
* The dynamics might be [[:Category:Unstable | unstable]].<br />
* There might be an underlying [[:Category:Network topology | network topology]].<br />
* The integer control functions might have been (re)formulated by means of an [[:Category:Outer convexification|outer convexification]].<br />
<br />
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.<br />
<br />
[[Category:Model characterization]]</div>SebastianSagerhttps://mintoc.de/index.php?title=LinearMetabolic&diff=2344LinearMetabolic2023-10-19T11:33:00Z<p>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..."</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 1<br />
|nc = 2<br />
|nri = 2<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
--><br />
<br />
Under construction...<br />
<br />
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 ...<br />
<br />
The problem is discussed in detail in <bib id="Sager2015" />.<br />
<br />
== Mathematical formulation ==<br />
<br />
The mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{x, u} & x_3(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_0 & = & (V_{alim} u - R_m x_0 - K_m x_1) / L_m, \\<br />
& \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)), \\<br />
& \dot{x}_2 & = & \frac{r}{K_r} x_1, \\<br />
& \dot{x}_3 & = & V_{alim} u x_0 + R_{bat} x_0^2, \\[1.5ex]<br />
& x(t_0) &=& (0,0,0,0)^T, \\<br />
& x(t_f) & \in & \mathcal{T} \subseteq \mathbb{R}^4,\\<br />
& x_0(t) & \in & [-i_{max}, i_{max}], \\<br />
& u(t) &\in& \{-1, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
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.<br />
<br />
== Parameters ==<br />
These fixed values are used within the model.<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|Coefficient of reduction<br />
|<math>K_r</math><br />
|10<br />
|[-]<br />
|-<br />
|Air density<br />
|<math>\rho</math><br />
|1.293<br />
|<math> kg/m^3 </math><br />
|-<br />
|Aerodynamic coefficient<br />
|<math>C_x</math><br />
|0.4<br />
|[-]<br />
|-<br />
|Area in the front of the vehicle<br />
|<math>S</math><br />
|2<br />
|<math> m^2 </math><br />
|-<br />
|Radius of the wheel<br />
|<math>r</math><br />
|0.33<br />
|<math> m </math><br />
|-<br />
|Constant representing the friction of the wheels on the road<br />
|<math>K_f </math><br />
|0.03<br />
|[-]<br />
|-<br />
|Coefficient of the motor torque<br />
|<math>K_m</math><br />
|0.27<br />
|[-]<br />
|-<br />
|Inductor resistance<br />
|<math>R_m</math><br />
|0.03<br />
|<math> Ohms </math><br />
|-<br />
|inductance of the rotor<br />
|<math>L_m</math><br />
|0.05<br />
|[-]<br />
|-<br />
|Mass<br />
|<math>M</math><br />
|250<br />
|<math> kg </math><br />
|-<br />
|Gravity constant<br />
|<math>g</math><br />
|9,81<br />
|[-]<br />
|-<br />
|Battery voltage<br />
|<math>V_{alim}</math><br />
|150<br />
|<math> V </math><br />
|-<br />
|Resistance of the battery<br />
|<math>R_{bat}</math><br />
|0.05<br />
|<math> Ohms </math><br />
|-<br />
|}<br />
<br />
== Reference Solutions ==<br />
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.<br />
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.<br />
<br />
Applying the sum up rounding strategy results in an integer-feasible chattering solution.<br />
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.<br />
<br />
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>.<br />
<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="3"><br />
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>.<br />
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.<br />
Image:ElectricCar PrimalStates SumUpRounding.PNG| Primal states with Sum up rounding on grid with <math> \Delta t = 10^{-2} </math>. <br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL code]] at [[Electric car problem (AMPL)]]<br />
<br />
== Variants ==<br />
<br />
There are several alternative formulations and variants of the above problem, in particular<br />
* fixed final velocity, <math> \mathcal{T} = \mathbb{R} \times \{ 50 \frac{K_r}{3.6r} \} \times \{100\} \times \mathbb{R}, t_f = 10s </math><br />
* 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><br />
* 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>.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:gIOC]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=External_Links&diff=2342External Links2022-03-02T14:38:48Z<p>SebastianSager: /* Work Groups */</p>
<hr />
<div>= Work Groups =<br />
<br />
An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.<br />
<br />
* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager<br />
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl<br />
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts<br />
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin<br />
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich<br />
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell<br />
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt<br />
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler<br />
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering<br />
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann<br />
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson<br />
* [https://www.tu-braunschweig.de/mo/team/kirches Mathematical Optimization], led by Christian Kirches<br />
* [http://www.mathematik.tu-clausthal.de/arbeitsgruppen/kontinuierliche-optimierung/ Continuous Optimization], led by Andreas Potschka<br />
* [https://www.mathematik.hu-berlin.de/de/forschung/forschungsgebiete/mathematische-optimierung/Falk-hante/falk-hante-fp Applied Mathematics], led by Falk Hante<br />
<br />
Feel free to update this list!<br />
<br />
= Benchmark Libraries =<br />
<br />
=== Embedded Optimization and Control ===<br />
<br />
* European FP 7 project [http://embocon.org EMBOCON]<br />
<br />
=== Hybrid Systems ===<br />
<br />
* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.<br />
<br />
=== Optimal Control ===<br />
<br />
* 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. <br />
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems. <br />
* 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.<br />
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems<br />
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization<br />
* 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<br />
<br />
=== MINLP ===<br />
<br />
* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]<br />
* MINLPlib <bib id="Bussieck2003a" /> <br />
* [http://www.gamsworld.org/performance GAMSWORLD]<br />
<br />
=== Generic Optimization ===<br />
<br />
* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.<br />
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)<br />
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)<br />
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)<br />
<br />
== References ==<br />
<biblist /></div>SebastianSagerhttps://mintoc.de/index.php?title=Oil_Shale_Pyrolysis&diff=1663Oil Shale Pyrolysis2016-02-04T19:40:52Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nw = 0<br />
|nre = 2<br />
}}<br />
<br />
The following problem is an example from the global optimal control literature and was introduced in <bib id="Wen1977" />. <br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
<math><br />
\begin{array}{ll}<br />
\displaystyle \min_{u} & \displaystyle -x_1(t_N)^2 \\[1.5ex]<br />
\mbox{s.t.} & \displaystyle \dot{x}_0(t) = -k_0x_0(t)-(k_2+k_3+k_4)x_0(t)x_1(t)\\<br />
& \displaystyle \dot{x}_1(t) = k_0x_0(t)-k_1x_1(t) + k_2x_0(t)x_1(t)\\<br />
& \displaystyle k_i = a_i e^{-u(t)\frac{b_i}{R}},\quad \forall i\in \{1,\dots,5\} \\ [1.5ex]<br />
& \displaystyle t \in \left[t_0,t_N\right] \\<br />
& \displaystyle u(t) \in \left[698.15/748.15,1\right]\\<br />
& \displaystyle x(t_0) = (1,0)^T\\<br />
\end{array} <br />
</math><br />
<br />
where this is the normalized form with <br />
<br />
<math> u(t)= \frac{1}{u_{temp}} </math>, with <br />
<br />
<math> u_{temp} \in \left[698.15,748.15\right] </math> <br />
<br />
== Parameters ==<br />
<br />
<br />
<br />
{| class="wikitable"<br />
|+State variables<br />
|-<br />
|Symbol<br />
|Initial value (<math>t_0</math>)<br />
|-<br />
|<math>x_0(t)</math><br />
|<math>1</math><br />
|-<br />
|<math>x_1(t)</math><br />
|<math>0 </math><br />
|}<br />
<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Symbol<br />
|Value<br />
|-<br />
|<math>a_1</math><br />
|<math>8.86</math><br />
|-<br />
|<math>a_2</math><br />
|<math>24.25</math><br />
|-<br />
|<math>a_3</math><br />
|<math>23.67</math><br />
|-<br />
|<math>a_4</math><br />
|<math>18.75</math><br />
|-<br />
|<math>a_5</math><br />
|<math>20.7</math><br />
|-<br />
|<math>b_1</math><br />
|<math>20.3</math><br />
|-<br />
|<math>b_2</math><br />
|<math>37.4</math><br />
|-<br />
|<math>b_3</math><br />
|<math>33.8</math><br />
|-<br />
|<math>b_4</math><br />
|<math>28.2</math><br />
|-<br />
|<math>b_5</math><br />
|<math>31.0</math><br />
|}<br />
<br />
{| class="wikitable"<br />
|+Control variable<br />
|-<br />
|Symbol<br />
|Interval<br />
|-<br />
|<math>u(t)</math><br />
|[698.15/748.15,1]<br />
|}<br />
<br />
'''Measurement grid'''<br />
<br />
== Reference solution ==<br />
<br />
Coming soon.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are not yet available.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, or presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem&diff=1662Lotka Volterra fishing problem2016-02-04T19:38:59Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nw = 1<br />
|nre = 3<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->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. <br />
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.<br />
<br />
<br />
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.<br />
<br />
The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0(t) & = & x_0(t) - x_0(t) x_1(t) - \; c_0 x_0(t) \; w(t), \\<br />
& \dot{x}_1(t) & = & - x_1(t) + x_0(t) x_1(t) - \; c_1 x_1(t) \; w(t), \\<br />
& \dot{x}_2(t) & = & (x_0(t) - 1)^2 + (x_1(t) - 1)^2, \\[1.5ex]<br />
& x(0) &=& (0.5, 0.7, 0)^T, \\<br />
& w(t) &\in& \{0, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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>.<br />
<br />
== Parameters ==<br />
<br />
These fixed values are used within the model.<br />
<br />
<math><br />
\begin{array}{rcl}<br />
[t_0, t_f] &=& [0, 12],\\<br />
(c_0, c_1) &=& (0.4, 0.2).<br />
\end{array}<br />
</math><br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
<br />
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>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.<br />
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.<br />
Image:lotka2Switches.png| Control and differential states with only two switches.<br />
Image:lotka56Switches.png| Control and differential states with 56 switches.<br />
</gallery><br />
<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]<br />
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]<br />
* [[:Category:optimica | optimica code]] at [[Lotka Volterra fishing problem (optimica)]]<br />
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]<br />
* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]<br />
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]<br />
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]<br />
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]<br />
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]<br />
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]<br />
<br />
== Variants ==<br />
<br />
There are several alternative formulations and variants of the above problem, in particular<br />
<br />
* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],<br />
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,<br />
* 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>,<br />
* different fishing control functions for the two species,<br />
* different parameters and start values.<br />
<br />
== Miscellaneous and Further Reading ==<br />
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.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Tracking objective]]<br />
[[Category:Chattering]]<br />
[[Category:Sensitivity-seeking arcs]]<br />
[[Category:Population dynamics]]<br />
<br />
<br />
<!--Testing Graphviz<br />
<graphviz border='frame' format='svg'><br />
digraph G {Hello->World!}<br />
</graphviz><br />
--></div>SebastianSagerhttps://mintoc.de/index.php?title=External_Links&diff=1591External Links2016-01-28T10:09:13Z<p>SebastianSager: </p>
<hr />
<div>= Work Groups =<br />
<br />
An incomplete list of work groups interested in aspects of mixed-integer dynamic optimization.<br />
<br />
* [http://mathopt.de Mathematical and Computational Optimization], Uni Magdeburg, led by Sebastian Sager<br />
* [https://www.imtek.de/professuren/systemtheorie Systems Control and Optimization Laboratory], Uni Freiburg, led by Moritz Diehl<br />
* [http://www.unibw.de/lrt1/gerdts Optimal Control], Universität der Bundeswehr München, led by Matthias Gerdts<br />
* [http://www.am.uni-erlangen.de/wima Discrete Optimization], Uni Erlangen, led by Alexander Martin<br />
* [http://www3.mathematik.tu-darmstadt.de/ags/optimization/research/nonlinear-optimization.html Nonlinear Optimization], Uni Darmstadt, led by Stefan Ulbrich<br />
* [http://www.dyn.bci.tu-dortmund.de/ Process Dynamics and Operations], Uni Dortmund, led by Sebastian Engell<br />
* [http://www.aices.rwth-aachen.de/organization/institutes/procsyseng Process Systems Engineering], RWTH Aachen, led by Wolfgang Marquardt<br />
* [http://numero.cheme.cmu.edu/ Biegler Research Group], Carnegie Mellon, led by Larry Biegler<br />
* [http://www.am.uni-erlangen.de/home/leugering/ Applied Mathematics II], Uni Erlangen, led by Günter Leugering<br />
* [http://egon.cheme.cmu.edu/ Grossmann Research Group], Carnegie Mellon, led by Ignacio Grossmann<br />
* [http://www.control.lth.se/user/johan.akesson/ Automatic Control], Lund University, led by Johan Akesson<br />
<br />
Feel free to update this list!<br />
<br />
= Benchmark Libraries =<br />
<br />
=== Embedded Optimization and Control ===<br />
<br />
* European FP 7 project [http://embocon.org EMBOCON]<br />
<br />
=== Hybrid Systems ===<br />
<br />
* [http://www.ist-hycon.org/ HYCON] network: two benchmark problems have been defined.<br />
<br />
=== Optimal Control ===<br />
<br />
* 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. <br />
* With the software package [http://abs-5.me.washington.edu/noc/dsoa.html dsoa] come currently 77 test problems. <br />
* 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.<br />
* The [http://www.mcs.anl.gov/~more/cops/ COPS library] contains around 10 dynamic control and parameter estimation problems<br />
* The [http://www.optpde.net OPTPDE collection] contains problems in PDE-constrained optimization<br />
* 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<br />
<br />
=== MINLP ===<br />
<br />
* CMU-IBM Cyber-Infrastructure for MINLP [http://minlp.org/ collaborative site]<br />
* MINLPlib <bib id="Bussieck2003a" /> <br />
* [http://www.gamsworld.org/performance GAMSWORLD]<br />
<br />
=== Generic Optimization ===<br />
<br />
* [http://cuter.rl.ac.uk/cuter-www/ CUTEr] is a versatile testing environment for optimization and linear algebra solvers.<br />
* [http://www.netlib.org/lp/ NETLIB] for linear programming (LP)<br />
* [http://miplib.zib.de/ MIPLIB] for mixed-integer linear programming (MILP)<br />
* Schittkowski Library <bib id="Schittkowski2002" /> for nonlinear programming (LP)<br />
<br />
== References ==<br />
<biblist /></div>SebastianSagerhttps://mintoc.de/index.php?title=Main_Page&diff=1590Main Page2016-01-28T10:01:22Z<p>SebastianSager: </p>
<hr />
<div>__NOTOC__<br />
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.<br />
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.<br />
<br />
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.<br />
<br />
<!-- TOP TABLE WITH NEWS --><br />
==[[:Category:News|News]] <span style="font-variant:small-caps" style="font-size:10px">[[Current News|(add)]]</span>==<br />
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px"><br />
<tr valign="top"><br />
<td width="80%" bgcolor="#EEEEFF"><br />
<!-- The actual news are being included from the page Current News with the following line. --><br />
{{Current News}}<br />
</td><br />
</tr><br />
</table><br />
<br />
<table width="100%"><br />
<br />
<!-- CATEGORY FIELD--><br />
<tr valign="top"><br />
<td width="50%"><br />
==[[:Category:Problem characterization|Problem characterization]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a problem characterization|(add)]]</span>==<br />
<br />
*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]<br />
<br />
*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]]<br />
<br />
*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]<br />
<br />
*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category:Casadi | Casadi]] - [[:Category: Switch | Switch]]<br />
</td><br />
<br />
<!-- PROBLEM FIELD--><br />
<td width="50%" ><br />
<br />
==Direct links to [[:Category:MIOCP|Problems]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a Problem|(add)]]</span>==<br />
<br />
[[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]]''<br />
</td><br />
</tr><br />
<br />
<tr valign="top"><br />
<br />
<!-- HELP FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Help|Help]] <span style="font-variant:small-caps" style="font-size:10px">[[User:SebastianSager|(contact)]]</span>==<br />
<br />
''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]<br />
</td><br />
<br />
<!-- COMMUNITY FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Community|Community]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a User|(add)]]</span>==<br />
<br />
[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''<br />
<br />
</td><br />
</tr><br />
<br />
</table><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Main]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Main_Page&diff=1589Main Page2016-01-28T09:53:28Z<p>SebastianSager: </p>
<hr />
<div>__NOTOC__<br />
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.<br />
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.<br />
<br />
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.<br />
<br />
<!-- TOP TABLE WITH NEWS --><br />
==[[:Category:News|News]] <span style="font-variant:small-caps" style="font-size:10px">[[Current News|(add)]]</span>==<br />
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px"><br />
<tr valign="top"><br />
<td width="80%" bgcolor="#EEEEFF"><br />
<!-- The actual news are being included from the page Current News with the following line. --><br />
{{Current News}}<br />
</td><br />
</tr><br />
</table><br />
<br />
<table width="100%"><br />
<br />
<!-- CATEGORY FIELD--><br />
<tr valign="top"><br />
<td width="50%"><br />
==[[:Category:Problem characterization|Problem characterization]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a problem characterization|(add)]]</span>==<br />
<br />
*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Chemical engineering | Chemical engineering]] - [[:Category:Medicine | Medicine]] - [[:Category:Transport | Transport]]<br />
<br />
*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]]<br />
<br />
*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]<br />
<br />
*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category: Muscod | Muscod]] - [[:Category: Switch | Switch]]<br />
</td><br />
<br />
<!-- PROBLEM FIELD--><br />
<td width="50%" ><br />
<br />
==Direct links to [[:Category:MIOCP|Problems]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a Problem|(add)]]</span>==<br />
<br />
[[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]]''<br />
</td><br />
</tr><br />
<br />
<tr valign="top"><br />
<br />
<!-- HELP FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Help|Help]] <span style="font-variant:small-caps" style="font-size:10px">[[User:SebastianSager|(contact)]]</span>==<br />
<br />
''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]<br />
</td><br />
<br />
<!-- COMMUNITY FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Community|Community]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a User|(add)]]</span>==<br />
<br />
[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''<br />
<br />
</td><br />
</tr><br />
<br />
</table><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Main]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Annihilation_of_calcium_oscillations_with_PLC_activation_inhibition&diff=1588Annihilation of calcium oscillations with PLC activation inhibition2016-01-28T09:51:48Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 2<br />
|nre = 4<br />
}}<br />
<br />
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.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<math><br />
\begin{array}{llcl}<br />
\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]<br />
\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} \\<br />
& \dot{x}_1 & = & (1 - w_2) k_7 x_0 - \frac{k_8 x_1}{x_1 + K_9} \\<br />
& \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} \\<br />
& \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]<br />
& x(0) &=& (0.03966, 1.09799, 0.00142, 1.65431)^T, \\<br />
& x(t) & \ge & 0.0, \\<br />
& w_1(t) &\in& \{1, w^{\mathrm{max}}\}, \\<br />
& w_2(t) &\in& \{0, 1\}, \\<br />
& w^{\mathrm{max}} & \ge & 1.1, \\<br />
& w^{\mathrm{max}} & \le & 1.3.<br />
\end{array} <br />
</math><br />
<br />
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>. <br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:calcium2Controls.png| Optimal stimuli determined by mixed-integer optimal control.<br />
Image:calcium2States.png| Corresponding differential states.<br />
</gallery><br />
<br />
== Variants ==<br />
<br />
[[Annihilation of calcium oscillations| Only one control function]].<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, or presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:Medicine]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Unstable]]<br />
[[Category:Tracking objective]]<br />
[[Category:Bang bang]]<br />
[[Category:Systems biology]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Annihilation_of_calcium_oscillations&diff=1587Annihilation of calcium oscillations2016-01-28T09:51:36Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 1<br />
|nre = 4<br />
}}<br />
<br />
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.<br />
<br />
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]].<br />
<br />
== Biological motivation ==<br />
<br />
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.<br />
<br />
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.<br />
<br />
As a source of external control a temporally varying concentration of an uncompetitive inhibitor of the PMCA ion pump is considered.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\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]<br />
\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} \\<br />
& \dot{x}_1 & = & k_7 x_0 - \frac{k_8 x_1}{x_1 + K_9} \\<br />
& \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} \\<br />
& \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]<br />
& x(0) &=& (0.03966, 1.09799, 0.00142, 1.65431)^T, \\<br />
& x(t) & \ge & 0.0, \\<br />
& w(t) &\in& \{1, w^{\mathrm{max}}\}, \\<br />
& w^{\mathrm{max}} & \ge & 1.1, \\<br />
& w^{\mathrm{max}} & \le & 1.3.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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. <br />
<br />
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.<br />
<br />
== Parameters ==<br />
<br />
These fixed values are used within the model.<br />
<br />
<p><br />
<math><br />
\begin{array}{rcl}<br />
[t_0, t_f] &=& [0, 22],\\<br />
k_1 &=& 0.09, \\<br />
k_2 &=& 2.30066, \\<br />
k_3 &=& 0.64, \\<br />
K_4 &=& 0.19, \\<br />
k_5 &=& 4.88, \\<br />
K_6 &=& 1.18, \\<br />
k_7 &=& 2.08, \\<br />
k_8 &=& 32.24, \\<br />
K_9 &=& 29.09, \\<br />
k_{10} &=& 5.0, \\<br />
K_{11} &=& 2.67, \\<br />
k_{12} &=& 0.7, \\<br />
k_{13} &=& 13.58, \\<br />
k_{14} &=& 153.0, \\<br />
K_{15} &=& 0.16, \\<br />
k_{16} &=& 4.85, \\<br />
K_{17} &=& 0.05, \\<br />
p_1 &=& 100, \\<br />
\tilde{x}_0 &=& 6.78677, \\<br />
\tilde{x}_1 &=& 22.65836, \\<br />
\tilde{x}_2 &=& 0.384306, \\<br />
\tilde{x}_3 &=& 0.28977.<br />
\end{array}<br />
</math><br />
</p><br />
<br />
== Reference Solutions ==<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:calciumControl.png| Optimal stimuli determined by mixed-integer optimal control.<br />
Image:calciumStates.png| Corresponding differential states.<br />
Image:calciumStatesPerturbed.png| Slightly perturbed control: stimulus 0.001 too early.<br />
Image:calciumStatesLongTime.png| Long time behavior of optimal solution: numerical rounding errors lead to transition back from unstable steady-state to stable limit-cycle.<br />
</gallery><br />
<br />
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.<br />
<br />
== Optimization ==<br />
<br />
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.<br />
<br />
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.<br />
<br />
<gallery caption="Objective landscape simulation plots" widths="180px" heights="140px" perrow="2"><br />
Image:calciumObjectiveSimulation1.png| Brute-force simulation of objective landscape for different values of stimulus beginning and length.<br />
Image:calciumObjectiveSimulation2.png| Zoom in.<br />
Image:calciumObjectiveSimulation3.png| Zoom in.<br />
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.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[Annihilation of calcium oscillations (Muscod)]]<br />
* [[:Category:optimica | optimica]] at [[Annihilation of calcium oscillations (optimica)]]<br />
<br />
== Variants ==<br />
<br />
<ul><br />
<li> [[Annihilation of calcium oscillations with PLC activation inhibition| Use of two control functions]].<br />
<li> Scaled deviation in objective function.<br />
</ul><br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, or presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:Medicine]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Unstable]]<br />
[[Category:Bang bang]]<br />
[[Category:Systems biology]]<br />
[[Category:Tracking objective]]</div>SebastianSagerhttps://mintoc.de/index.php?title=De_Pillis_chemotherapy_model&diff=1586De Pillis chemotherapy model2016-01-28T09:51:12Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 6<br />
|nu = 3<br />
}}<br />
<br />
The model by de Pillis combines chemotherapy with immunotherapy.<br />
<br />
== Mathematical formulation ==<br />
<br />
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<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\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]<br />
\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 ,\\<br />
& \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, \\<br />
& \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 \\<br />
& & & - 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, \\<br />
& \dot{x}_3(t) & = & \alpha - \beta x_3 - K_C (1- \mbox{e}^{- x_4}) x_3,\\<br />
& \dot{x}_4(t) & = & - \gamma x_4 + u_0,\\<br />
& \dot{x}_5(t) & = & - \mu_I x_5 + u_1.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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.<br />
<br />
== Parameters ==<br />
<br />
This set of parameters can be found as “patient 9” in [].<br />
<br />
<math><br />
\begin{array}{lll}<br />
a = 4.31 \times 10^{-1} & b = 1.02 \times 10^{-9} & c = 6.41 \times 10^{-11}\\<br />
d = 2.34 & e = 2.08 \times 10^{-7} & f = 4.12 \times 10^{-2}\\<br />
g = 1.25 \times 10^{-2} & h = 2.02 \times 10^{7} & j = 2.49 \times 10^{-2}\\<br />
k = 3.66 \times 10^{7} & l = 2.09 & m = 2.04 \times 10^{-1}\\<br />
q = 1.42 \times 10^{-6} & p = 3.42 \times 10^{-6} & s = 8.39 \times 10^{-2}\\<br />
r_1 = 1.01 \times 10^{-7} & r_2 = 6.50 \times 10^{-11} & u = 3.00 \times 10^{-10}\\<br />
\alpha = 7.50 \times 10^{8} & \beta = 1.20 \times 10^{-2} & \gamma = 9.00 \times 10^{-1}\\<br />
p_I = 1.25 \times 10^{-1} & g_I = 2.00 \times 10^{7} & \mu_I = 1.00 \times 10^{1}\\<br />
K_T = 9.00 \times 10^{-1} & K_N = K_L = K_C = 6 \times 10^{-1}<br />
\end{array}<br />
</math><br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
<br />
Objective function 1<br />
<math><br />
\begin{array}{rclr}<br />
p_0 &=& 1,&\\<br />
p_2 &=& 1,&\\<br />
p_i &=& 0,&otherwise.\\<br />
\end{array}<br />
</math><br />
<br />
Objective function 2<br />
<math><br />
\begin{array}{rclr}<br />
p_0 &=& -1,&\\<br />
p_2 &=& 1,&\\<br />
p_i &=& 0,&otherwise.\\<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
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.<br />
Image:DePillis_2.png| Optimal controls and states with objective function 2.<br />
</gallery><br />
<br />
==Source Code==<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[De Pillis chemotherapy model (Muscod)]]<br />
<br />
== References ==<br />
<br />
[[Category:Medicine]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Category:Medicine&diff=1585Category:Medicine2016-01-28T09:50:15Z<p>SebastianSager: Initial setup</p>
<hr />
<div>This category includes all control problems that stem from a medical application. <br />
<br />
[[Category:Application]]</div>SebastianSagerhttps://mintoc.de/index.php?title=D%27Onofrio_chemotherapy_model&diff=1584D'Onofrio chemotherapy model2016-01-28T09:49:51Z<p>SebastianSager: added medicine</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nu = 2<br />
}}<br />
<br />
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]].<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> the optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} & x_0(t_f) &+& \alpha \int_{t_0}^{t_f} u_0(t)^2 \text{d}t \\[1.5ex]<br />
\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), \\<br />
& \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), \\<br />
& \dot{x}_2(t) & = & u_0(t), \\<br />
& \dot{x}_3(t) & = & u_1(t), \\ [1.5ex]<br />
& 0 & \leq & u_0(t) \leq u_0^{max}, \\<br />
& 0 & \leq & u_1(t) \leq u_1^{max}, \\<br />
& x_2(t) & \leq & x_2^{max}, \\<br />
& x_3(t) & \leq & x_3^{max}.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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.<br />
<br />
== Parameters ==<br />
<br />
In the model these parameters are fixed.<br />
<br />
<math><br />
\begin{array}{rcl}<br />
t_0 &=& 0,\\<br />
(\zeta, b, \mu, d, G) &=& (0.192, 5.85, 0.0, 0.00873, 0.15),\\<br />
(x_2(0), x_3(0), u_0^{max}, x_2^{max}) &=& (0,0,75,300).<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Reference Solutions ==<br />
<br />
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 <br />
<br />
Parameter set 1<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\<br />
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 2<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\<br />
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 3<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\<br />
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 4<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\<br />
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\<br />
\end{array}<br />
</math><br />
<br />
Furthermore in the objective function <math>\alpha =0</math> is chosen.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
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.<br />
Image:d'Onofrio 2.png| Optimal controls and states with parameter set 2.<br />
Image:d'Onofrio 3.png| Optimal controls and states with parameter set 3.<br />
Image:d'Onofrio 4.png| Optimal controls and states with parameter set 4.<br />
</gallery><br />
<br />
==Source Code==<br />
<br />
[[:Category:Muscod | Muscod code]] <br />
<source lang="cpp"><br />
/* volume of tumor */<br />
rhs[0] = - zeta*x[0]*log(x[0]/x[1])-F*x[0]*u1_;<br />
/* volume of neighboring blood vessels */<br />
rhs[1] = b*x[0] - (mu + d*pow(x[0],2.0/3.0))*x[1] - G*u0_*x[1] - eta*x[1]*u1_;<br />
/* amount of anti-angiogetic drug */<br />
rhs[2] = u0_;<br />
/* amount of cytostatic drug */<br />
rhs[3] = u1_;<br />
</source><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Medicine]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Cushioned_Oscillation&diff=1583Cushioned Oscillation2016-01-28T09:41:42Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nre = 4<br />
}}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.<br />
<br />
== Model formulation == <br />
<br />
An object with mass <math> m </math> is attached to a spring with stiffness constant <math> c </math>.<br />
<br />
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:<br />
<br />
<br />
<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>)<br />
<br />
<br />
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.<br />
<br />
Through external force, the object has been put into an initial state :<br />
<br />
<math>(x(0),v(0)) = (x_0,v_0)</math><br />
<br />
The goal is to reset position and velocity of the object as fast as possible, meaning:<br />
<br />
<math>(x(t_f),v(t_f)) = (0,0)</math>,<br />
<br />
with the objective function:<br />
<br />
<math>\min\limits_{t_f} t_f</math><br />
<br />
== Optimal Control Problem (OCP) Formulation ==<br />
<br />
The above results in the following OCP <br />
<br />
<math> \begin{array}{llr}<br />
\min\limits_{x,v,u,t_f} & t_f\\ <br />
<br />
s.t. & \dot x (t) = v(t), & \forall t \in [0,t_f]\\<br />
<br />
& \dot v (t)= \frac{1}{m}(u(t) - cx(t)), & \forall t \in [0,t_f]\\<br />
\\ <br />
& x(0)=x_0,\\<br />
& v(0)=v_0,\\<br />
& x(t_f)=0,\\<br />
& v(t_f)=0,\\<br />
\\<br />
& -u_{mm} \le u(t) \le u_{mm}, & \forall t \in [0,t_f]\\ <br />
\end{array}</math><br />
<br />
== Parameters and Reference Solution ==<br />
<br />
The following parameters were used, to create the reference solution below, with an almost optimal final time <math> t_f = 8.98 s</math>:<br />
<br />
<math> m=5, </math><br />
<math> c=10, </math><br />
<math> x_0=2, </math><br />
<math> v_0=5, </math><br />
<math> u_{mm}=5.</math><br />
<br />
== Reference Solution ==<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="1"><br />
Image:Ref_sol_plot_cushioned_oscillation_m5.png| States and Controls<br />
</gallery><br />
<br />
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).<br />
<br />
== Source Code ==<br />
<br />
* A MATLAB script using [[:Category:TomDyn/PROPT | PROPT]] can be found in: [[Cushioned Oscillation (PROPT)]]<br />
<br />
== References ==<br />
<br />
[[Category:MIOCP]] <br />
[[Category: Minimum time]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Cushioned_Oscillation&diff=1582Cushioned Oscillation2016-01-28T09:40:21Z<p>SebastianSager: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nre = 4<br />
}}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.<br />
<br />
== Model formulation == <br />
<br />
An object with mass <math> m </math> is attached to a spring with stiffness constant <math> c </math>.<br />
<br />
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:<br />
<br />
<br />
<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>)<br />
<br />
<br />
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.<br />
<br />
Through external force, the object has been put into an initial state :<br />
<br />
<math>(x(0),v(0)) = (x_0,v_0)</math><br />
<br />
The goal is to reset position and velocity of the object as fast as possible, meaning:<br />
<br />
<math>(x(t_f),v(t_f)) = (0,0)</math>,<br />
<br />
with the objective function:<br />
<br />
<math>\min\limits_{t_f} t_f</math><br />
<br />
== Optimal Control Problem (OCP) Formulation ==<br />
<br />
The above results in the following OCP <br />
<br />
<math> \begin{array}{llr}<br />
\min\limits_{x,v,u,t_f} & t_f\\ <br />
<br />
s.t. & \dot x (t) = v(t), & \forall t \in [0,t_f]\\<br />
<br />
& \dot v (t)= \frac{1}{m}(u(t) - cx(t)), & \forall t \in [0,t_f]\\<br />
\\ <br />
& x(0)=x_0,\\<br />
& v(0)=v_0,\\<br />
& x(t_f)=0,\\<br />
& v(t_f)=0,\\<br />
\\<br />
& -u_{mm} \le u(t) \le u_{mm}, & \forall t \in [0,t_f]\\ <br />
\end{array}</math><br />
<br />
== Parameters and Reference Solution ==<br />
<br />
The following parameters were used, to create the reference solution below, with an almost optimal final time <math> t_f = 8.98 s</math>:<br />
<br />
<math> m=5, </math><br />
<math> c=10, </math><br />
<math> x_0=2, </math><br />
<math> v_0=5, </math><br />
<math> u_{mm}=5.</math><br />
<br />
== Reference Solution ==<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="1"><br />
Image:Ref_sol_plot_cushioned_oscillation_m5.png| States and Controls<br />
</gallery><br />
<br />
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).<br />
<br />
== Source Code ==<br />
<br />
* A MATLAB script using [[Category:TomDyn/PROPT]] can be found in: [[Cushioned Oscillation (PROPT)]]<br />
<br />
== References ==<br />
<br />
[[Category:MIOCP]] <br />
[[Category: Minimum time]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Main_Page&diff=1581Main Page2016-01-28T09:33:52Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>__NOTOC__<br />
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.<br />
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.<br />
<br />
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.<br />
<br />
<!-- TOP TABLE WITH NEWS --><br />
==[[:Category:News|News]] <span style="font-variant:small-caps" style="font-size:10px">[[Current News|(add)]]</span>==<br />
<table bgcolor="#EEEEFF" width="100%" cellspacing="10px"><br />
<tr valign="top"><br />
<td width="80%" bgcolor="#EEEEFF"><br />
<!-- The actual news are being included from the page Current News with the following line. --><br />
{{Current News}}<br />
</td><br />
</tr><br />
</table><br />
<br />
<table width="100%"><br />
<br />
<!-- CATEGORY FIELD--><br />
<tr valign="top"><br />
<td width="50%"><br />
==[[:Category:Problem characterization|Problem characterization]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a problem characterization|(add)]]</span>==<br />
<br />
*via ''[[:Category:Application|application area]]'', e.g., [[:Category: Aeronautics | Aeronautics]] - [[:Category: Population dynamics | Population dynamics]] - [[:Category: Chemical engineering | Chemical engineering]]<br />
<br />
*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]]<br />
<br />
*via ''[[:Category:Solution characterization|optimal solution]]'', e.g., [[:Category: Bang bang | Bang bang]] - [[:Category: Chattering | Chattering]] - [[:Category: Sensitivity-seeking arcs | Sensitivity-seeking arcs]]<br />
<br />
*via ''[[:Category: Implementation | Implementation]], e.g., [[:Category: ACADO | ACADO]] - [[:Category: AMPL | AMPL]] - [[:Category: AMPL/TACO | AMPL with TACO]] - [[:Category: Muscod | Muscod]] - [[:Category: Switch | Switch]]<br />
</td><br />
<br />
<!-- PROBLEM FIELD--><br />
<td width="50%" ><br />
<br />
==Direct links to [[:Category:MIOCP|Problems]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a Problem|(add)]]</span>==<br />
<br />
[[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]]''<br />
</td><br />
</tr><br />
<br />
<tr valign="top"><br />
<br />
<!-- HELP FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Help|Help]] <span style="font-variant:small-caps" style="font-size:10px">[[User:SebastianSager|(contact)]]</span>==<br />
<br />
''[[Help:How to contribute|How to contribute]]'' - ''[[Help:How to cite|How to cite]]'' - ''[[Help:Using LaTeX|LaTeX]]'' - ''[[Help:Description guidelines | Guidelines]]<br />
</td><br />
<br />
<!-- COMMUNITY FIELD--><br />
<td width="50%" ><br />
<br />
==[[:Category:Community|Community]] <span style="font-variant:small-caps" style="font-size:10px">[[Help:Adding a User|(add)]]</span>==<br />
<br />
[[:Category:Community|Contributors]] - ''[[External Links]]'' - ''Feel encouraged to participate!''<br />
<br />
</td><br />
</tr><br />
<br />
</table><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:Main]]</div>SebastianSagerhttps://mintoc.de/index.php?title=De_Pillis_chemotherapy_model&diff=1572De Pillis chemotherapy model2016-01-28T09:29:49Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 6<br />
|nu = 3<br />
}}<br />
<br />
The model by de Pillis combines chemotherapy with immunotherapy.<br />
<br />
== Mathematical formulation ==<br />
<br />
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<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\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]<br />
\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 ,\\<br />
& \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, \\<br />
& \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 \\<br />
& & & - 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, \\<br />
& \dot{x}_3(t) & = & \alpha - \beta x_3 - K_C (1- \mbox{e}^{- x_4}) x_3,\\<br />
& \dot{x}_4(t) & = & - \gamma x_4 + u_0,\\<br />
& \dot{x}_5(t) & = & - \mu_I x_5 + u_1.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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.<br />
<br />
== Parameters ==<br />
<br />
This set of parameters can be found as “patient 9” in [].<br />
<br />
<math><br />
\begin{array}{lll}<br />
a = 4.31 \times 10^{-1} & b = 1.02 \times 10^{-9} & c = 6.41 \times 10^{-11}\\<br />
d = 2.34 & e = 2.08 \times 10^{-7} & f = 4.12 \times 10^{-2}\\<br />
g = 1.25 \times 10^{-2} & h = 2.02 \times 10^{7} & j = 2.49 \times 10^{-2}\\<br />
k = 3.66 \times 10^{7} & l = 2.09 & m = 2.04 \times 10^{-1}\\<br />
q = 1.42 \times 10^{-6} & p = 3.42 \times 10^{-6} & s = 8.39 \times 10^{-2}\\<br />
r_1 = 1.01 \times 10^{-7} & r_2 = 6.50 \times 10^{-11} & u = 3.00 \times 10^{-10}\\<br />
\alpha = 7.50 \times 10^{8} & \beta = 1.20 \times 10^{-2} & \gamma = 9.00 \times 10^{-1}\\<br />
p_I = 1.25 \times 10^{-1} & g_I = 2.00 \times 10^{7} & \mu_I = 1.00 \times 10^{1}\\<br />
K_T = 9.00 \times 10^{-1} & K_N = K_L = K_C = 6 \times 10^{-1}<br />
\end{array}<br />
</math><br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
<br />
Objective function 1<br />
<math><br />
\begin{array}{rclr}<br />
p_0 &=& 1,&\\<br />
p_2 &=& 1,&\\<br />
p_i &=& 0,&otherwise.\\<br />
\end{array}<br />
</math><br />
<br />
Objective function 2<br />
<math><br />
\begin{array}{rclr}<br />
p_0 &=& -1,&\\<br />
p_2 &=& 1,&\\<br />
p_i &=& 0,&otherwise.\\<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
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.<br />
Image:DePillis_2.png| Optimal controls and states with objective function 2.<br />
</gallery><br />
<br />
==Source Code==<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[De Pillis chemotherapy model (Muscod)]]<br />
<br />
== References ==<br />
<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]</div>SebastianSagerhttps://mintoc.de/index.php?title=De_Pillis_chemotherapy_model_(Muscod)&diff=1570De Pillis chemotherapy model (Muscod)2016-01-28T09:29:12Z<p>SebastianSager: SebastianSager moved page De Pillis chemotherapy model (C) to De Pillis chemotherapy model (Muscod): Moved C to Muscod</p>
<hr />
<div>This is source code written in [[:Category:Muscod | Muscod]] for the [[De Pillis chemotherapy model]].<br />
<br />
<source lang="cpp"><br />
/* volume of tumor in absolute cell count*/<br />
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];<br />
<br />
/* volume of unspecific immune cells */<br />
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];<br />
<br />
/* volume of tumor-specific cytotoxic T-cells */<br />
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] ...<br />
- 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];<br />
<br />
/* volume of circulating lymphocytes */<br />
rhs[3] = alpha - beta *x[3] - K_C *(1- exp(- x[4])) *x[3];<br />
<br />
/* amount of chemotherapeutic drug */<br />
rhs[4] = u[0] - gamma*x[4];<br />
<br />
/* amount of immunotherapeutic drug */<br />
rhs[5] = u[1] - mu_I*x[5];<br />
<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=De_Pillis_chemotherapy_model_(C)&diff=1571De Pillis chemotherapy model (C)2016-01-28T09:29:12Z<p>SebastianSager: SebastianSager moved page De Pillis chemotherapy model (C) to De Pillis chemotherapy model (Muscod): Moved C to Muscod</p>
<hr />
<div>#REDIRECT [[De Pillis chemotherapy model (Muscod)]]</div>SebastianSagerhttps://mintoc.de/index.php?title=De_Pillis_chemotherapy_model_(Muscod)&diff=1569De Pillis chemotherapy model (Muscod)2016-01-28T09:29:01Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>This is source code written in [[:Category:Muscod | Muscod]] for the [[De Pillis chemotherapy model]].<br />
<br />
<source lang="cpp"><br />
/* volume of tumor in absolute cell count*/<br />
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];<br />
<br />
/* volume of unspecific immune cells */<br />
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];<br />
<br />
/* volume of tumor-specific cytotoxic T-cells */<br />
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] ...<br />
- 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];<br />
<br />
/* volume of circulating lymphocytes */<br />
rhs[3] = alpha - beta *x[3] - K_C *(1- exp(- x[4])) *x[3];<br />
<br />
/* amount of chemotherapeutic drug */<br />
rhs[4] = u[0] - gamma*x[4];<br />
<br />
/* amount of immunotherapeutic drug */<br />
rhs[5] = u[1] - mu_I*x[5];<br />
<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=D%27Onofrio_chemotherapy_model&diff=1568D'Onofrio chemotherapy model2016-01-28T09:27:13Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nu = 2<br />
}}<br />
<br />
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]].<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> the optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} & x_0(t_f) &+& \alpha \int_{t_0}^{t_f} u_0(t)^2 \text{d}t \\[1.5ex]<br />
\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), \\<br />
& \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), \\<br />
& \dot{x}_2(t) & = & u_0(t), \\<br />
& \dot{x}_3(t) & = & u_1(t), \\ [1.5ex]<br />
& 0 & \leq & u_0(t) \leq u_0^{max}, \\<br />
& 0 & \leq & u_1(t) \leq u_1^{max}, \\<br />
& x_2(t) & \leq & x_2^{max}, \\<br />
& x_3(t) & \leq & x_3^{max}.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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.<br />
<br />
== Parameters ==<br />
<br />
In the model these parameters are fixed.<br />
<br />
<math><br />
\begin{array}{rcl}<br />
t_0 &=& 0,\\<br />
(\zeta, b, \mu, d, G) &=& (0.192, 5.85, 0.0, 0.00873, 0.15),\\<br />
(x_2(0), x_3(0), u_0^{max}, x_2^{max}) &=& (0,0,75,300).<br />
\end{array}<br />
</math><br />
<br />
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.<br />
<br />
== Reference Solutions ==<br />
<br />
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 <br />
<br />
Parameter set 1<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\<br />
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 2<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 12000,& x_1(0) &=& 15000,\\<br />
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 3<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\<br />
u_1^{max} &=& 1,& x_3^{max} &=& 2.\\<br />
\end{array}<br />
</math><br />
<br />
Parameter set 4<br />
<br />
<math><br />
\begin{array}{rclrcl}<br />
x_0(0) &=& 14000,& x_1(0) &=& 5000,\\<br />
u_1^{max} &=& 2,& x_3^{max} &=& 10.\\<br />
\end{array}<br />
</math><br />
<br />
Furthermore in the objective function <math>\alpha =0</math> is chosen.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
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.<br />
Image:d'Onofrio 2.png| Optimal controls and states with parameter set 2.<br />
Image:d'Onofrio 3.png| Optimal controls and states with parameter set 3.<br />
Image:d'Onofrio 4.png| Optimal controls and states with parameter set 4.<br />
</gallery><br />
<br />
==Source Code==<br />
<br />
[[:Category:Muscod | Muscod code]] <br />
<source lang="cpp"><br />
/* volume of tumor */<br />
rhs[0] = - zeta*x[0]*log(x[0]/x[1])-F*x[0]*u1_;<br />
/* volume of neighboring blood vessels */<br />
rhs[1] = b*x[0] - (mu + d*pow(x[0],2.0/3.0))*x[1] - G*u0_*x[1] - eta*x[1]*u1_;<br />
/* amount of anti-angiogetic drug */<br />
rhs[2] = u0_;<br />
/* amount of cytostatic drug */<br />
rhs[3] = u1_;<br />
</source><br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category: MIOCP]]<br />
[[Category:ODE model]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem_(TACO)&diff=1567Lotka Volterra fishing problem (TACO)2016-01-28T09:25:35Z<p>SebastianSager: </p>
<hr />
<div>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" />.<br />
Note that you will need to include a generic [[support AMPL files|AMPL/TACO support file]], OptimalControl.mod.<br />
To solve this model, you require an optimal control or MINLP code that uses the TACO toolkit to support the AMPL optimal control extensions.<br />
<br />
=== AMPL ===<br />
<br />
This is the source file lotka_taco.mod<br />
<source lang="AMPL"><br />
# ----------------------------------------------------------------<br />
# Lotka Volterra fishing problem using AMPL and TACO<br />
# (c) Christian Kirches, Sven Leyffer<br />
# ----------------------------------------------------------------<br />
include OptimalControl.mod;<br />
<br />
var t;<br />
var xd0 := 0.5, >= 0, <= 20;<br />
var xd1 := 0.7, >= 0, <= 20;<br />
var dev := 0.0, >= 0, <= 20;<br />
<br />
var u := 1, >= 0, <= 1 integer suffix type "u0";<br />
<br />
param p0 := 0.4;<br />
param p1 := 0.2;<br />
param ref0 := 1.0;<br />
param ref1 := 1.0;<br />
<br />
# Minimize accumulated deviation from reference after 12 time units<br />
minimize Mayer: eval(dev,12);<br />
<br />
subject to<br />
<br />
# ODE system<br />
ODE_0: diff(xd0,t) = xd0 - xd0*xd1 - p0*u*xd0; # prey<br />
ODE_1: diff(xd1,t) = -xd1 + xd0*xd1 - p1*u*xd1; # predator<br />
ODE_2: diff(dev,t) = (xd0-ref0)^2 + (xd1-ref1)^2; # deviation from reference<br />
<br />
# initial value constraints<br />
IVC_0: eval(xd0,0) = 0.5;<br />
IVC_1: eval(xd1,0) = 0.7;<br />
IVC_2: eval(dev,0) = 0.0;<br />
<br />
option solver ...;<br />
<br />
solve;<br />
<br />
</source><br />
<br />
== References ==<br />
<biblist /> <br />
<br />
[[Category:AMPL/TACO]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Supermarket_refrigeration_system&diff=1566Supermarket refrigeration system2016-01-28T09:24:31Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 9<br />
|nw = 3<br />
|nre = 9<br />
}}<br />
<br />
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.<br />
<br />
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.<br />
<br />
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. <br />
<br />
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.<br />
<br />
The optimal integer control function shows [[:Category:Chattering|chattering]] behavior, making the supermarket refrigeration system problem a candidate for benchmarking of algorithms.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\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 <br />
</math><br />
</p><br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle <br />
\mbox{s.t.} & <br />
\dot{x_0}(t) &=& \dfrac{1}{V_{suc} \cdot \frac{d\rho_{suc}}{dP_{suc}}(x_0)} \cdot \bigg[ <br />
\left(\dfrac{UA_{wall-ref, max}}{M_{ref, max} \cdot<br />
\Delta h_{lg}(x_0)}\right) \Big( x_4 \big( x_2 - T_e(x_0) \big)\\ <br />
& && + \, 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) <br />
\bigg] \\<br />
&\dot{x_1}(t) &=& - \dfrac{<br />
UA_{goods-air} \left( x_1 - x_3 \right)<br />
}{<br />
M_{goods} \cdot C_{p,goods} <br />
} \\ <br />
&\dot{x_2}(t) &=& \dfrac{<br />
UA_{air-wall} \left( x_3-x_2 \right)<br />
- \dfrac{UA_{wall-ref,max}}{M_{ref,max}}<br />
\, x_4 \big( x_2 - T_e(x_0) \big)<br />
}{<br />
M_{wall} \cdot C_{p,wall}<br />
} \\ [2.5ex]<br />
&\dot{x_3}(t) &=& \dfrac{<br />
UA_{goods-air} \left( x_1-x_3 \right) + \dot{Q}_{airload}<br />
- UA_{air-wall} \, (x_3-x_2)<br />
}{<br />
M_{air} \cdot C_{p,air}<br />
} \\ [2.5ex]<br />
&\dot{x_4}(t) &=& \left(\dfrac{M_{ref,max} - x_4}{\tau_{fill}} \right) u_0<br />
- \left( \dfrac{UA_{wall-ref,max}}{M_{ref,max} \cdot \Delta h_{lg}(x_0)} \,<br />
x_4 \big( x_2 - T_e(x_0) \big) \right) (1-u_0)<br />
\\ \\<br />
&\dot{x_5}(t) &=& - \dfrac{<br />
UA_{goods-air} \left( x_5 - x_7 \right)<br />
}{<br />
M_{goods} \cdot C_{p,goods} <br />
} \\ <br />
&\dot{x_6}(t) &=& \dfrac{<br />
UA_{air-wall} \left( x_7-x_6 \right)<br />
- \dfrac{UA_{wall-ref,max}}{M_{ref,max}}<br />
\, x_8 \big( x_6 - T_e(x_0) \big)<br />
}{<br />
M_{wall} \cdot C_{p,wall}<br />
} \\ [2.5ex]<br />
&\dot{x_7}(t) &=& \dfrac{<br />
UA_{goods-air} \left( x_5-x_7 \right) + \dot{Q}_{airload}<br />
- UA_{air-wall} \, (x_7-x_6)<br />
}{<br />
M_{air} \cdot C_{p,air}<br />
} \\ [2.5ex]<br />
&\dot{x_8}(t) &=& \left(\dfrac{M_{ref,max} - x_8}{\tau_{fill}} \right) u_1<br />
- \left( \dfrac{UA_{wall-ref,max}}{M_{ref,max} \cdot \Delta h_{lg}(x_0)} \,<br />
x_8 \big( x_6 - T_e(x_0) \big) \right) (1-u_1)<br />
\\ [4ex]<br />
<br />
& x_3(t) &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_3(t) &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_7(t) &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_7(t) &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_0(t) &\leq& 1.7 \quad \forall t \in [t_0, t_f], \\<br />
& x_i(t_0) &=& free \quad \forall i \in \{0,\dots, 8\}, \\<br />
& x_i(t_f) &=& x_i(t_0) \quad \forall i \in \{0,\dots, 8\}, \\<br />
& u_i(t) &\in& \{0, 1\} \quad \forall i \in \{0,\dots, 2\}, \\<br />
& t_f &\in& [ 650, 750 ]. <br />
<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
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).<br />
<br />
<math>(x_5, x_6, x_7, x_8)</math> describe the corresponding states in the second display case.<br />
<br />
<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.<br />
<math>u_2, u_3</math> denote the activity of a single compressor.<br />
<br />
The following polynomial functions are used in the model description originating from interpolations:<br />
<br />
<p><br />
<math><br />
\begin{array}{rcl}<br />
<br />
T_e(x_0) &=& -4.3544 x_0^2 + 29.224 x_0 - 51.2005,\\<br />
\Delta h_{lg}(x_0) &=& (0.0217 x_0^2 - 0.1704 x_0 + 2.2988)\cdot 10^5, \\<br />
\rho_{suc}(x_0) &=& 4.6073 x_0 + 0.3798, \\<br />
\frac{d\rho_{suc}}{dP_{suc}}(x_0) &=& -0.0329 {x_0}^3 + 0.2161 {x_0}^2 - 0.4742 x_0 + 5.4817,\\<br />
f(x_0) &=& (0.0265 x_0^3 - 0.4346 x_0^2 + 2.4923 x_0 + 1.2189) \cdot 10^5.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
These fixed values are used within the model for the day scenario. A night scenario is also available, see [[#Variants|Variants]].<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Symbol !! Value !! Unit !! Description<br />
|- <br />
| align=center | <math>\dot{Q}_{airload}</math> || align=right | 3000.00 || <math>\frac{J}{s}</math> || Disturbance, heat transfer from outside the display case<br />
|-<br />
| align=center | <math>\dot{m}_{ref,const}</math> || align=right | 0.20 || <math>\frac{kg}{s}</math> || Disturbance, constant mass flow of refrigerant<br />
from unmodeled entities <br />
|-<br />
| align=center | <math>M_{goods}</math> || align=right | 200.00 || <math>kg</math> || Mass of goods <br />
|-<br />
| align=center | <math>C_{p,goods}</math> || align=right | 1000.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of goods<br />
|-<br />
| align=center | <math><br />
UA_{goods-air} </math> || align=right | 300.00 || <math>\frac{J}{s \cdot K}</math> || Heat transfer coefficient between goods<br />
and air <br />
|-<br />
| align=center | <math>M_{wall} </math> || align=right | 260.00 || <math>kg</math> || Mass of evaporator wall <br />
|-<br />
| align=center | <math>C_{p,wall} </math> || align=right | 385.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of evaporator wall<br />
|-<br />
| align=center | <math>UA_{air-wall}</math> || align=right | 500.00 || <math>\frac{J}{s \cdot K}</math> || Heat transfer coefficient between air and<br />
wall<br />
|-<br />
| align=center | <math>M_{air}</math> || align=right | 50.00 || <math>kg</math> || Mass of air in display case <br />
|-<br />
| align=center | <math>C_{p,air}</math> || align=right | 1000.00 || <math>\frac{J}{kg \cdot K}</math> || Heat capacity of air<br />
|-<br />
| align=center | <math>UA_{wall-ref,max}</math> || align=right | 4000.00 || <math>\frac{J}{s \cdot K}</math> || Maximum heat transfer coefficient between<br />
refrigerant and evaporator wall<br />
|-<br />
| align=center | <math>\tau_{fill}</math> || align=right | 40.00 || <math>s</math> || Parameter describing the filling time of the<br />
evaporator under opened valve <br />
|-<br />
| align=center | <math>T_{SH}</math> || align=right | 10.00 || <math>K</math> || Superheat in the suction manifold<br />
|-<br />
| align=center | <math>M_{ref,max}</math> || align=right | 1.00 || <math>kg</math> || Maximum mass of refrigerant in evaporator<br />
|-<br />
| align=center | <math>V_{suc}</math> || align=right | 5.00 || <math>m^3</math> || Total volume of suction manifold<br />
|-<br />
| align=center | <math>V_{sl}</math> || align=right | 0.08 || <math>\frac{m^3}{s}</math> || Total displacement volume<br />
|-<br />
| align=center | <math>\eta_{vol}</math> || align=right | 0.81 || <math>-</math> || Volumetric efficiency<br />
|}<br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
The illustrated solution with integer controls has a (suboptimal) objective function value of 12252.81.<br />
<br />
<gallery caption="Reference solution plots" widths="350px" heights="300px" perrow="2"><br />
Image:FridgeControlsRelaxed.png| Optimal relaxed control.<br />
Image:FridgeControls.png| Optimal integer control.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[Supermarket refrigeration system (Muscod)]]<br />
* [[:Category:optimica | optimica]] at [[Supermarket refrigeration system (optimica)]]<br />
<br />
== Variants ==<br />
<br />
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.<br />
<br />
In the paper <bib id="Larsen2007" /> mentioned above, the problem was stated slightly different:<br />
<br />
* 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.<br />
* 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:<br />
<math><br />
\dot{x_4} = \begin{cases}<br />
<br />
\dfrac{M_{ref,max} - x_4}{\tau_{fill}} & \qquad \text{if} \quad u_0 = 1 \\ \\<br />
- \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 \\ \\ <br />
0 & \qquad \text{if} \quad u_0 = 0 \quad \text{and} \quad x_4 = 0<br />
<br />
\end{cases}<br />
</math><br />
<br />
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>.<br />
* A night scenario with two different parameters was given. At night the following parameters change their value:<br />
<math><br />
\begin{array}{lcrr}<br />
\dot{Q}_{airload} &=& 1800.00 & \frac{J}{s}, \\<br />
\dot{m}_{ref,const} &=& 0.00 & \frac{kg}{s}, \\<br />
\end{array}<br />
</math><br />
<br />
Additionally the constraint on the suction pressure <math>x_0(t)</math> is softened to <math>x_0(t) \leq 1.9</math>.<br />
* 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.<br />
* 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.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Chattering]]<br />
[[Category:Periodic]]<br />
[[Category:Minimum energy]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Supermarket_refrigeration_system_(Muscod)&diff=1565Supermarket refrigeration system (Muscod)2016-01-28T09:24:03Z<p>SebastianSager: SebastianSager moved page Supermarket refrigeration system (C) to Supermarket refrigeration system (Muscod) without leaving a redirect: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Supermarket refrigeration system]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
// number of compressors<br />
#define NVALVES 2<br />
<br />
// constants<br />
#define M_GOODS 200.0<br />
#define C_P_GOODS 1000.0<br />
#define UA_GOODS_AIR 300.0<br />
#define M_WALL 260.0<br />
#define C_P_WALL 385.0<br />
#define UA_AIR_WALL 500.0<br />
#define M_AIR 50.0<br />
#define C_P_AIR 1000.0<br />
#define UA_WALL_REF_MAX 4000.0<br />
#define M_REF_MAX 1.0<br />
#define TAU_FILL 40.0<br />
#define T_SH 10.0<br />
#define V_SUC 5.0<br />
#define V_SL 0.08 // 2 display cases - 2 compressors<br />
// #define V_SL 0.095 // 3 display cases - 3 compressors<br />
#define ETA_VOL 0.81<br />
<br />
// disturbances - day scenario<br />
#define Q_AIRLOAD 3000.0<br />
#define M_REF_CONST 0.2<br />
<br />
// disturbances - night scenario<br />
// #define Q_AIRLOAD 1800.0<br />
// #define M_REF_CONST 0.0<br />
<br />
double delta_h = (0.0217*xd[0]*xd[0] - 0.1704*xd[0] + 2.2988)*1e5;<br />
double T_e = -4.3544*xd[0]*xd[0] + 29.224*xd[0] - 51.2005;<br />
double rho_suc = 4.6073*xd[0] + 0.3798;<br />
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;<br />
double f = (0.0265*xd[0]*xd[0]*xd[0] -0.4346*xd[0]*xd[0] + 2.4923*xd[0] + 1.2189)*1e5;<br />
<br />
double Q_goods_air[NVALVES];<br />
double Q_air_wall[NVALVES];<br />
double UA_wall_ref[NVALVES];<br />
double Q_e[NVALVES];<br />
double m[NVALVES];<br />
<br />
double m_in_suc = 0.0;<br />
<br />
int i;<br />
<br />
for (i=0; i<NVALVES; i++){<br />
Q_goods_air[i] = UA_GOODS_AIR*(xd[1 + i*4] - xd[3 + i*4]);<br />
Q_air_wall[i] = UA_AIR_WALL*(xd[3 + i*4] - xd[2 + i*4]);<br />
UA_wall_ref[i] = UA_WALL_REF_MAX * xd[4 + 4*i]/M_REF_MAX;<br />
Q_e[i] = UA_wall_ref[i]*(xd[2 + 4*i] - T_e);<br />
<br />
m[i] = Q_e[i]/delta_h;<br />
m_in_suc += m[i];<br />
}<br />
<br />
double V_comp = 0.0;<br />
double comp_scale = (double) 1.0/NCOMPS;<br />
V_comp = comp_scale*u[NVALVES]*ETA_VOL*V_SL;<br />
<br />
// suction pressure<br />
rhs[0] = (m_in_suc + M_REF_CONST - V_comp*rho_suc) / (V_SUC*rho_suc__P_suc);<br />
<br />
// for each display/valve<br />
for (i=0; i<NVALVES; i++){<br />
<br />
// temperatures:<br />
<br />
// goods<br />
rhs[1 + i*4] = - Q_goods_air[i]/(M_GOODS*C_P_GOODS);<br />
<br />
// wall<br />
rhs[2 + i*4] = (Q_air_wall[i] - Q_e[i])/(M_WALL*C_P_WALL);<br />
<br />
// air<br />
rhs[3 + i*4] = ((Q_goods_air[i] + Q_AIRLOAD - Q_air_wall[i]) /(M_AIR*C_P_AIR));<br />
<br />
// mass of liquefied refrigerant:<br />
rhs[4 + i*4] = ((M_REF_MAX - xd[4 + 4*i])/TAU_FILL * u[i] - m[i] * (1 - u[i]));<br />
}<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Supermarket_refrigeration_system_(Muscod)&diff=1564Supermarket refrigeration system (Muscod)2016-01-28T09:23:53Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Supermarket refrigeration system]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
// number of compressors<br />
#define NVALVES 2<br />
<br />
// constants<br />
#define M_GOODS 200.0<br />
#define C_P_GOODS 1000.0<br />
#define UA_GOODS_AIR 300.0<br />
#define M_WALL 260.0<br />
#define C_P_WALL 385.0<br />
#define UA_AIR_WALL 500.0<br />
#define M_AIR 50.0<br />
#define C_P_AIR 1000.0<br />
#define UA_WALL_REF_MAX 4000.0<br />
#define M_REF_MAX 1.0<br />
#define TAU_FILL 40.0<br />
#define T_SH 10.0<br />
#define V_SUC 5.0<br />
#define V_SL 0.08 // 2 display cases - 2 compressors<br />
// #define V_SL 0.095 // 3 display cases - 3 compressors<br />
#define ETA_VOL 0.81<br />
<br />
// disturbances - day scenario<br />
#define Q_AIRLOAD 3000.0<br />
#define M_REF_CONST 0.2<br />
<br />
// disturbances - night scenario<br />
// #define Q_AIRLOAD 1800.0<br />
// #define M_REF_CONST 0.0<br />
<br />
double delta_h = (0.0217*xd[0]*xd[0] - 0.1704*xd[0] + 2.2988)*1e5;<br />
double T_e = -4.3544*xd[0]*xd[0] + 29.224*xd[0] - 51.2005;<br />
double rho_suc = 4.6073*xd[0] + 0.3798;<br />
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;<br />
double f = (0.0265*xd[0]*xd[0]*xd[0] -0.4346*xd[0]*xd[0] + 2.4923*xd[0] + 1.2189)*1e5;<br />
<br />
double Q_goods_air[NVALVES];<br />
double Q_air_wall[NVALVES];<br />
double UA_wall_ref[NVALVES];<br />
double Q_e[NVALVES];<br />
double m[NVALVES];<br />
<br />
double m_in_suc = 0.0;<br />
<br />
int i;<br />
<br />
for (i=0; i<NVALVES; i++){<br />
Q_goods_air[i] = UA_GOODS_AIR*(xd[1 + i*4] - xd[3 + i*4]);<br />
Q_air_wall[i] = UA_AIR_WALL*(xd[3 + i*4] - xd[2 + i*4]);<br />
UA_wall_ref[i] = UA_WALL_REF_MAX * xd[4 + 4*i]/M_REF_MAX;<br />
Q_e[i] = UA_wall_ref[i]*(xd[2 + 4*i] - T_e);<br />
<br />
m[i] = Q_e[i]/delta_h;<br />
m_in_suc += m[i];<br />
}<br />
<br />
double V_comp = 0.0;<br />
double comp_scale = (double) 1.0/NCOMPS;<br />
V_comp = comp_scale*u[NVALVES]*ETA_VOL*V_SL;<br />
<br />
// suction pressure<br />
rhs[0] = (m_in_suc + M_REF_CONST - V_comp*rho_suc) / (V_SUC*rho_suc__P_suc);<br />
<br />
// for each display/valve<br />
for (i=0; i<NVALVES; i++){<br />
<br />
// temperatures:<br />
<br />
// goods<br />
rhs[1 + i*4] = - Q_goods_air[i]/(M_GOODS*C_P_GOODS);<br />
<br />
// wall<br />
rhs[2 + i*4] = (Q_air_wall[i] - Q_e[i])/(M_WALL*C_P_WALL);<br />
<br />
// air<br />
rhs[3 + i*4] = ((Q_goods_air[i] + Q_AIRLOAD - Q_air_wall[i]) /(M_AIR*C_P_AIR));<br />
<br />
// mass of liquefied refrigerant:<br />
rhs[4 + i*4] = ((M_REF_MAX - xd[4 + 4*i])/TAU_FILL * u[i] - m[i] * (1 - u[i]));<br />
}<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Subway_ride&diff=1563Subway ride2016-01-28T09:23:18Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nw = 1<br />
|nre = 4<br />
}}<br />
<br />
The ''subway ride'' optimal control problem goes back to work of <bib id="Bock1982" /> for the city of New York. In an extension, also velocity limits that lead to path-constrained arcs appear.<br />
The aim is to minimize the energy used for a subway ride from one station to another, taking into account boundary conditions and a restriction on the time. <br />
<br />
== Mathematical formulation ==<br />
<br />
The MIOCP reads as<br />
<br />
<center><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w} & & & \int_{0}^{t_f} L(x, w) \; \mathrm{d} t \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\<br />
& \dot{x}_1(t) & = & f_1(x,w), \\<br />
& x(0) &=& (0, 0)^T, \\<br />
& x(t_f) &=& (2112, 0)^T, \\<br />
& w(t) &\in& \{1, 2, 3, 4\}, \quad t \in [0, t_f].<br />
\end{array} <br />
</math><br />
</center><br />
<br />
The terminal time <math>t_f=65</math> denotes the time of arrival of a subway train in the next station. The differential states <math>x_0(\cdot)</math> and <math>x_1(\cdot)</math> describe position and velocity of the train, respectively. The train can be operated in one of four different modes, <math>w(\cdot) = 1</math> series, <math>w(\cdot) = 2</math> parallel, <math>w(\cdot) = 3</math> coasting, or <math>w(\cdot) = 4</math> braking that accelerate or decelerate the train and have different energy consumption. <br />
<br />
Acceleration and energy comsumption are velocity-dependent. Hence, we will need switching functions <math>\sigma_i( x_1 ) = v_i - x_1</math> for given velocities <math>v_i, i=1..3</math>.<br />
<br />
The Lagrange term reads as<br />
<center><math><br />
\begin{array}{rcl}<br />
L(x, 1) &=& \left\{ \begin{array}{cl} e \; p_1 & \mbox{if } \sigma_1 \ge 0 \\ e \; p_2 & \mbox{else if } \sigma_2 \ge 0 \\ e \; \sum_{i=0}^{5} c_i (1) \left( \frac{1}{10} \gamma\ x_1 \right)^{-i} \quad & \mbox{else} \end{array} \right. \\<br />
L(x, 2) &=& \left\{ \begin{array}{cl} \infty & \mbox{if } \sigma_2 \ge 0\\ e \; p_3 & \mbox{else if } \sigma_3 \ge 0 \\ <br />
e \; \sum_{i=0}^{5} c_i (2) \left( \frac{1}{10} \gamma\ x_1 - 1 \right)^{-i} & \mbox{else} \end{array} \right. \\<br />
L(x, 3) &=& L(x, 4) = 0.<br />
\end{array}<br />
</math></center><br />
<br />
The right hand side function <math>f_1(x,w)</math> reads as<br />
<center><math><br />
\begin{array}{rcl}<br />
f_1(x, 1) &=& \left\{ \begin{array}{cl} f_1^{1A} := \frac{g \; e \; a_1}{W_{\mathrm{eff}}} & \mbox{if } \sigma_1 \ge 0 \\ f_1^{1B} :=\frac{g \; e \; a_2}{W_{\mathrm{eff}}} & \mbox{else if } \sigma_2 \ge 0 \\ f_1^{1C} := \frac{g \; (e \; T(x_1, 1) - R(x_1)}{W_{\mathrm{eff}}} & \mbox{else} \end{array} \right. \\<br />
f_1(x, 2) &=& \left\{ \begin{array}{cl} 0 & \mbox{if } \sigma_2 \ge 0\\ f_1^{2B} := \frac{g \; e \; a_3}{W_{\mathrm{eff}}} & \mbox{else if } \sigma_3 \ge 0 \\ <br />
f_1^{2C} := \frac{g \; (e \; T(x_1, 2) - R(x_1)}{W_{\mathrm{eff}}} & \mbox{else} \end{array} \right. \\<br />
f_1(x, 3) &=& - \frac{g \; R(x_1)}{W_{\mathrm{eff}}} - C, \\<br />
f_1(x, 4) &=& - u = - u_{\mathrm{max}}.<br />
\end{array}<br />
</math></center><br />
<br />
The braking deceleration <math>u(\cdot)</math> can be varied between <math>0</math> and a given <math>u_{\mathrm{max}}</math>. It can be shown that only maximal braking can be optimal, hence we fixed <math>u(\cdot)</math> to <math>u_{\mathrm{max}}</math> without loss of generality. <br />
<br />
Occurring forces are<br />
<center><math><br />
\begin{array}{rcl}<br />
R(x_1) &=& c a \; {\gamma}^2 {x_1}^2 + b W {\gamma} x_1 + \frac{1.3}{2000}W + 116, \\<br />
T(x_1,1) &=& \sum_{i=0}^{5}b_i(1) \left( \frac{1}{10} {\gamma} x_1 - 0.3 \right) ^{-i}, \\<br />
T(x_1,2) &=& \sum_{i=0}^{5}b_i(2) \left( \frac{1}{10} {\gamma} x_1 - 1 \right) ^{-i}.<br />
\end{array}<br />
</math></center><br />
<br />
Details about the derivation of this model and the assumptions made can be found in <bib id="Bock1982" /> or in <bib id="Kraemer-Eis1985" />.<br />
<br />
== Parameters ==<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Symbol !! Value !! Unit !! Symbol !! Value !! Unit<br />
|- <br />
!<math>W</math> !! 78000 !! lbs !! <math>v_1</math> !! 0.979474 !! mph<br />
|- <br />
!<math>W_{\mathrm{eff}}</math> !! 85200 !! lbs !! <math>v_2</math> !! 6.73211 !! mph<br />
|- <br />
!<math>S</math> !! 2112 !! ft !! <math>v_3</math> !! 14.2658 !! mph<br />
|- <br />
!<math>S_4</math> !! 700 !! ft !! <math>v_4</math> !! 22.0 !! mph <br />
|- <br />
!<math>S_5</math> !! 1200 !! ft !! <math>v_5</math> !! 24.0 !! mph <br />
|- <br />
!<math>\gamma</math> !! <math>\frac{3600}{5280}</math> !! <math>\frac{\mbox{sec}}{\mbox{h}} / \frac{\mbox{ft}}{\mbox{mile}}</math> !! <math>a_1</math> !! 6017.611205 !! lbs <br />
|- <br />
!<math>a</math> !! 100 !! ft<math>^2</math> !! <math>a_2</math> !! 12348.34865 !! lbs <br />
|- <br />
!<math>n_{\mathrm{wag}}</math> !! 10 !! - !! <math>a_3</math> !! 11124.63729 !! lbs <br />
|- <br />
!<math>b</math> !! 0.045 !! - !! <math>u_{\mathrm{max}}</math> !! 4.4 !! ft / sec<math>^2</math> <br />
|- <br />
!<math>C</math> !! 0.367 !! - !! <math>p_1</math> !! 106.1951102 !! - <br />
|- <br />
!<math>g</math> !! 32.2 !! <math>\frac{\mbox{ft}}{\mbox{sec}^2}</math> !! <math>p_2</math> !! 180.9758408 !! - <br />
|- <br />
!<math>e</math> !! 1.0 !! - !! <math>p_3</math> !! 354.136479 !! - <br />
|}<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Parameter !! Value !! Parameter !! Value<br />
|- <br />
! b_0(1) !! -0.1983670410E02 !! c_0(1) !! 0.3629738340E02<br />
|- <br />
! b_1(1) !! 0.1952738055E03 !! c_1(1) !! -0.2115281047E03<br />
|- <br />
! b_2(1) !! 0.2061789974E04 !! c_2(1) !! 0.7488955419E03<br />
|- <br />
! b_3(1) !! -0.7684409308E03 !! c_3(1) !! -0.9511076467E03<br />
|- <br />
! b_4(1) !! 0.2677869201E03 !! c_4(1) !! 0.5710015123E03<br />
|- <br />
! b_5(1) !! -0.3159629687E02 !! c_5(1) !! -0.1221306465E03<br />
|- <br />
! b_0(2) !! -0.1577169936E03 !! c_0(2) !! 0.4120568887E02 <br />
|- <br />
! b_1(2) !! 0.3389010339E04 !! c_1(2) !! 0.3408049202E03 <br />
|- <br />
! b_2(2) !! 0.6202054610E04 !! c_1(2) !! 0.3408049202E03 <br />
|- <br />
! b_3(2) !! -0.4608734450E04 !! c_3(2) !! 0.8108316584E02 <br />
|- <br />
! b_4(2) !! 0.2207757061E04 !! c_4(2) !! -0.5689703073E01 <br />
|- <br />
! b_5(2) !! -0.3673344160E03 !! c_5(2) !! -0.2191905731E01<br />
|}<br />
<br />
== Reference Solutions ==<br />
<br />
The optimal trajectory for this problem has been calculated by means of an indirect approach in <bib id="Bock1982" /><bib id="Kraemer-Eis1985" />, and based on the direct multiple shooting method in <bib id="Sager2009" />. <br />
<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Time <math>t</math> !! <math>w(\cdot)</math> !! <math>f_1 = </math> !! <math>x_0</math> [ft] !! <math>x_1</math> [mph] !! <math>x_1</math> [ftps] !! Energy<br />
|- <br />
! <math>0.00000</math> !! 1 !! <math>f_1^{1A}</math> !! 0.0 !! 0.0 !! 0.0 !! 0.0<br />
|- <br />
! <math>0.63166</math> !! 1 !! <math>f_1^{1B}</math> !! 0.453711 !! 0.979474 !! 1.43656 !! 0.0186331<br />
|- <br />
! <math>2.43955</math> !! 1 !! <math>f_1^{1C}</math> !! 10.6776 !! 6.73211 !! 9.87375 !! 0.109518<br />
|- <br />
! <math>3.64338</math> !! 2 !! <math>f_1^{2B}</math> !! 24.4836 !! 8.65723 !! 12.6973 !! 0.147387<br />
|- <br />
! <math>5.59988</math> !! 2 !! <math>f_1^{2C}</math> !! 57.3729 !! 14.2658 !! 20.9232 !! 0.339851<br />
|- <br />
! <math>12.6070</math> !! 1 !! <math>f_1^{1C}</math> !! 277.711 !! 25.6452 !! 37.6129 !! 0.93519<br />
|- <br />
! <math>45.7827</math> !! 3 !! <math>f_1(3)</math> !! 1556.5 !! 26.8579 !! 39.3915 !! 1.14569<br />
|- <br />
! <math>46.8938</math> !! 3 !! <math>f_1(3)</math> !! 1600 !! 26.5306 !! 38.9115 !! 1.14569<br />
|- <br />
! <math>57.1600</math> !! 4 !! <math>f_1(4)</math> !! 1976.78 !! 23.5201 !! 34.4961 !! 1.14569<br />
|- <br />
! <math>65.0000</math> !! - !! <math>-</math> !! 2112 !! 0.0 !! 0.0 !! 1.14569<br />
|}<br />
<br />
<br />
== Variants ==<br />
<br />
The given parameters have to be modified to match different parts of the track, subway train types, or amount of passengers. A minimization of travel time might also be considered.<br />
<br />
The problem becomes more challenging, when additional point or path constraints are considered. <br />
<br />
=== Point constraint ===<br />
We consider the point constraint<br />
<center><math><br />
x_1 \le v_4 \mbox{ if } x_0 = S_4<br />
</math></center><br />
for a given distance <math>0 < S_4 < S</math> and velocity <math>v_4 > v_3</math>. Note that the state <math>x_0(\cdot)</math> is strictly monotonically increasing with time, as <math>\dot{x}_0 = x_1 > 0</math> for all <math>t \in (0, T)</math>.<br />
<br />
The optimal order of gears for <math>S_4 = 1200</math> and <math>v_4 = 22 / \gamma</math> with the additional interior point constraints (\ref{FASOPOINTCON}) is<br />
<math>1, 2, 1, 3, 4, 2, 1, 3, 4</math>. The stage lengths between switches are 2.86362, 10.722, 15.3108, 5.81821, 1.18383, 2.72451, 12.917, 5.47402, and 7.98594<br />
with <math>\Phi = 1.3978</math>. For different parameters <math>S_4 = 700</math> and <math>v_4 = 22 / \gamma</math> we obtain the gear choice 1, 2, 1, 3, 2, 1, 3, 4 and stage lengths<br />
2.98084, 6.28428, 11.0714, 4.77575, 6.0483, 18.6081, 6.4893, and 8.74202<br />
with <math>\Phi = 1.32518</math>.<br />
<br />
=== Path constraint ===<br />
A more practical restriction are path constraints on subsets of the track. We will consider a problem with additional path constraints<br />
<center><math><br />
x_1 \le v_5 \; \mbox{ if } \; x_0 \ge S_5.<br />
</math></center><br />
<br />
The additional path constraint changes the qualitative behavior of the relaxed solution. While all solutions considered this far were bang-bang and the main work consisted in finding the switching points, we now have a path-constrained arc. The optimal solutions for refined grids yield a series of monotonically decreasing objective function values, where the limit is the best value that can be approximated by an integer feasible solution. In our case we obtain 1.33108, 1.31070, 1.31058, 1.31058, ...<br />
<br />
The plot shows two possible integer realizations, with a trade-off between energy consumption and number of switches. Note that the solutions approximate the optimal driving behavior (a convex combination of two operation modes) by switching between the two and causing a touching of the velocity constraint from below as many times as we switch.<br />
<br />
<gallery caption="Reference solution plots" widths="271px" heights="379px" perrow="2"><br />
Image:subway78000PathConstraint1600-24Opt1.png| Optimal solution with 1 touch point.<br />
Image:subway78000PathConstraint1600-24Opt2.png| Optimal solution with 3 touch points.<br />
</gallery><br />
<br />
The differential state ''velocity'' of a subway train over time. The dotted vertical line indicates the beginning of the path constraint, the horizontal line the maximum velocity. Left: one switch leading to one touch point. Right: optimal solution for three switches. The energy-optimal solution needs to stay as close as possible to the maximum velocity on this time interval to avoid even higher energy-intensive accelerations in the start-up phase to match the terminal time constraint <math>t_f \le 65</math> to reach the next station.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[Subway ride (Muscod)]]<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:State dependent switches]]<br />
[[Category:Minimum energy]]<br />
[[Category:Chattering]]<br />
[[Category:Path-constrained arcs]]<br />
[[Category:Bang bang]]<br />
[[Category:Transport]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Subway_ride_(Muscod)&diff=1562Subway ride (Muscod)2016-01-28T09:22:49Z<p>SebastianSager: SebastianSager moved page Subway ride (C) to Subway ride (Muscod) without leaving a redirect: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Subway ride]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
<br />
#define V4 21.0/GAMMA<br />
#define S4 1000.0<br />
#define V5 22.0/GAMMA<br />
<br />
// LOCAL<br />
#define S 2112.0<br />
#define TEND 65.0<br />
//#define TEND 70.0<br />
// EXPRESS<br />
//#define S 6336.0<br />
//#define TEND 151.0<br />
/* User changes end */<br />
<br />
// Weight<br />
#define W 78000.0<br />
#define WEFF (W + 72000.0 * 0.1)<br />
/* 3600.0 / 5280.0 */<br />
#define GAMMA (3600.0 / 5280.0)<br />
#define A 100.0<br />
#define N 10.0<br />
#define B 0.045<br />
/* [ 0.24 + 0.034 * (N-1) ] / (100*N) */<br />
#define C 0.000546<br />
#define CC (0.367 * GAMMA)<br />
#define G (32.2 * GAMMA)<br />
#define E 1.0<br />
#define V1 (1.03 - (W - 72000) / ( 110000 - 72000) * (1.03 - 0.71) ) / GAMMA<br />
#define V2 (6.86 - (W - 72000) / ( 110000 - 72000) * (6.86 - 6.05) ) / GAMMA<br />
#define V3 (14.49 - (W - 72000) / ( 110000 - 72000) * (14.49 - 13.07) ) / GAMMA<br />
#define A1 (5998.6162 + (W - 72000) / ( 110000 - 72000) * (6118.9179 - 5998.6162) )<br />
#define A2 (11440.7968 + (W - 72000) / ( 110000 - 72000) * (17188.6252 - 11440.7968) )<br />
#define A3 (10280.0514 + (W - 72000) / ( 110000 - 72000) * (15629.0954 - 10280.0514) )<br />
#define P1 (105.880645 + (W - 72000) / ( 110000 - 72000) * (107.872258 - 105.880645) )<br />
#define P2 (168.931957 + (W - 72000) / ( 110000 - 72000) * (245.209888 - 168.931957) )<br />
#define P3 (334.626716 + (W - 72000) / ( 110000 - 72000) * (458.188550 - 334.626716) )<br />
<br />
#define B01 -0.1983670410E02<br />
#define B11 0.1952738055E03<br />
#define B21 0.2061789974E04<br />
#define B31 -0.7684409308E03<br />
#define B41 0.2677869201E03<br />
#define B51 -0.3159629687E02<br />
#define B02 -0.1577169936E03<br />
#define B12 0.3389010339E04<br />
#define B22 0.6202054610E04<br />
#define B32 -0.4608734450E04<br />
#define B42 0.2207757061E04<br />
#define B52 -0.3673344160E03<br />
<br />
#define C01 0.3629738340E02<br />
#define C11 -0.2115281047E03<br />
#define C21 0.7488955419E03<br />
#define C31 -0.9511076467E03<br />
#define C41 0.5710015123E03<br />
#define C51 -0.1221306465E03<br />
#define C02 0.4120568887E02<br />
#define C12 0.3408049202E03<br />
#define C22 -0.1436283271E03<br />
#define C32 0.8108316584E02<br />
#define C42 -0.5689703073E01<br />
#define C52 -0.2191905731E01<br />
<br />
/* -------------------------------------------- */<br />
<br />
static void mfcn(double *ts, double *xd, double *xa, double *p, double *pr, <br />
double *mval, long *dpnd, long *info)<br />
{<br />
if (*dpnd) { *dpnd = MFCN_DPND(0, 0, *xd, 0, 0); return; }<br />
*mval = xd[2];<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Seriell */<br />
static void ffcn1a( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A1 / WEFF / GAMMA;<br />
rhs[2] = P1;<br />
}<br />
static void ffcn1b( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A2 / WEFF / GAMMA;<br />
rhs[2] = P2;<br />
}<br />
static void ffcn1c( double *xd, double *rhs )<br />
{<br />
double T, R, temp, temp2, tempb, tempb2;<br />
<br />
rhs[0] = xd[1];<br />
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 0.3);<br />
temp2 = temp * temp;<br />
T = B01<br />
+ B11 * temp <br />
+ B21 * temp2<br />
+ B31 * temp2 * temp<br />
+ B41 * temp2 * temp2<br />
+ B51 * temp2 * temp2 * temp<br />
;<br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116; <br />
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;<br />
<br />
tempb = 1.0 / ( 0.1 * GAMMA * xd[1]);<br />
tempb2 = tempb * tempb;<br />
rhs[2] = ( C01 <br />
+ C11 * tempb <br />
+ C21 * tempb2<br />
+ C31 * tempb2 * tempb<br />
+ C41 * tempb2 * tempb2<br />
+ C51 * tempb2 * tempb2 * tempb )<br />
; <br />
if (xd[1] < V2 + 1e-2) {<br />
rhs[1] = G * E * A2 / WEFF / GAMMA;<br />
rhs[2] = P2;<br />
}<br />
<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Parallel */<br />
static void ffcn2b( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A3 / WEFF / GAMMA;<br />
rhs[2] = P3;<br />
}<br />
<br />
static void ffcn2c( double *xd, double *rhs )<br />
{<br />
double T, R, temp, temp2, tempb, tempb2;<br />
<br />
rhs[0] = xd[1];<br />
<br />
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 1.0);<br />
temp2 = temp * temp;<br />
T = B02<br />
+ B12 * temp <br />
+ B22 * temp2<br />
+ B32 * temp2 * temp<br />
+ B42 * temp2 * temp2<br />
+ B52 * temp2 * temp2 * temp<br />
; <br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;<br />
<br />
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;<br />
<br />
tempb = 1.0 / ( 0.1 * GAMMA * xd[1] - 1);<br />
tempb2 = tempb * tempb;<br />
rhs[2] = ( C02 <br />
+ C12 * tempb <br />
+ C22 * tempb2<br />
+ C32 * tempb2 * tempb<br />
+ C42 * tempb2 * tempb2<br />
+ C52 * tempb2 * tempb2 * tempb )<br />
; <br />
<br />
<br />
// Achtung<br />
if (xd[1] < V3 + 1e-2) {<br />
rhs[1] = G * E * A3 / WEFF / GAMMA;<br />
rhs[2] = P3;<br />
}<br />
<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Rollen */<br />
static void ffcn3( double *xd, double *rhs )<br />
{<br />
double R;<br />
<br />
rhs[0] = xd[1];<br />
<br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;<br />
rhs[1] = ( - G * R / WEFF - CC / GAMMA ) ;<br />
<br />
rhs[2] = 0;<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Bremsen */<br />
static void ffcn4( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = - 3.0 / GAMMA;<br />
rhs[2] = 0;<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for v < V1*/<br />
static void ffcn_A(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1a( xd, rhs0 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V1 < v < V2*/<br />
static void ffcn_B(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1b( xd, rhs0 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V2 < v < V3*/<br />
static void ffcn_C(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1c( xd, rhs0 );<br />
ffcn2b( xd, rhs1 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V3 < v */<br />
static void ffcn_D(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1c( xd, rhs0 );<br />
ffcn2c( xd, rhs1 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Pure coasting resp. braking on last stage to avoid nondifferentiabilities*/<br />
static void ffcn_E(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Subway_ride_(Muscod)&diff=1561Subway ride (Muscod)2016-01-28T09:22:34Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Subway ride]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
<br />
#define V4 21.0/GAMMA<br />
#define S4 1000.0<br />
#define V5 22.0/GAMMA<br />
<br />
// LOCAL<br />
#define S 2112.0<br />
#define TEND 65.0<br />
//#define TEND 70.0<br />
// EXPRESS<br />
//#define S 6336.0<br />
//#define TEND 151.0<br />
/* User changes end */<br />
<br />
// Weight<br />
#define W 78000.0<br />
#define WEFF (W + 72000.0 * 0.1)<br />
/* 3600.0 / 5280.0 */<br />
#define GAMMA (3600.0 / 5280.0)<br />
#define A 100.0<br />
#define N 10.0<br />
#define B 0.045<br />
/* [ 0.24 + 0.034 * (N-1) ] / (100*N) */<br />
#define C 0.000546<br />
#define CC (0.367 * GAMMA)<br />
#define G (32.2 * GAMMA)<br />
#define E 1.0<br />
#define V1 (1.03 - (W - 72000) / ( 110000 - 72000) * (1.03 - 0.71) ) / GAMMA<br />
#define V2 (6.86 - (W - 72000) / ( 110000 - 72000) * (6.86 - 6.05) ) / GAMMA<br />
#define V3 (14.49 - (W - 72000) / ( 110000 - 72000) * (14.49 - 13.07) ) / GAMMA<br />
#define A1 (5998.6162 + (W - 72000) / ( 110000 - 72000) * (6118.9179 - 5998.6162) )<br />
#define A2 (11440.7968 + (W - 72000) / ( 110000 - 72000) * (17188.6252 - 11440.7968) )<br />
#define A3 (10280.0514 + (W - 72000) / ( 110000 - 72000) * (15629.0954 - 10280.0514) )<br />
#define P1 (105.880645 + (W - 72000) / ( 110000 - 72000) * (107.872258 - 105.880645) )<br />
#define P2 (168.931957 + (W - 72000) / ( 110000 - 72000) * (245.209888 - 168.931957) )<br />
#define P3 (334.626716 + (W - 72000) / ( 110000 - 72000) * (458.188550 - 334.626716) )<br />
<br />
#define B01 -0.1983670410E02<br />
#define B11 0.1952738055E03<br />
#define B21 0.2061789974E04<br />
#define B31 -0.7684409308E03<br />
#define B41 0.2677869201E03<br />
#define B51 -0.3159629687E02<br />
#define B02 -0.1577169936E03<br />
#define B12 0.3389010339E04<br />
#define B22 0.6202054610E04<br />
#define B32 -0.4608734450E04<br />
#define B42 0.2207757061E04<br />
#define B52 -0.3673344160E03<br />
<br />
#define C01 0.3629738340E02<br />
#define C11 -0.2115281047E03<br />
#define C21 0.7488955419E03<br />
#define C31 -0.9511076467E03<br />
#define C41 0.5710015123E03<br />
#define C51 -0.1221306465E03<br />
#define C02 0.4120568887E02<br />
#define C12 0.3408049202E03<br />
#define C22 -0.1436283271E03<br />
#define C32 0.8108316584E02<br />
#define C42 -0.5689703073E01<br />
#define C52 -0.2191905731E01<br />
<br />
/* -------------------------------------------- */<br />
<br />
static void mfcn(double *ts, double *xd, double *xa, double *p, double *pr, <br />
double *mval, long *dpnd, long *info)<br />
{<br />
if (*dpnd) { *dpnd = MFCN_DPND(0, 0, *xd, 0, 0); return; }<br />
*mval = xd[2];<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Seriell */<br />
static void ffcn1a( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A1 / WEFF / GAMMA;<br />
rhs[2] = P1;<br />
}<br />
static void ffcn1b( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A2 / WEFF / GAMMA;<br />
rhs[2] = P2;<br />
}<br />
static void ffcn1c( double *xd, double *rhs )<br />
{<br />
double T, R, temp, temp2, tempb, tempb2;<br />
<br />
rhs[0] = xd[1];<br />
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 0.3);<br />
temp2 = temp * temp;<br />
T = B01<br />
+ B11 * temp <br />
+ B21 * temp2<br />
+ B31 * temp2 * temp<br />
+ B41 * temp2 * temp2<br />
+ B51 * temp2 * temp2 * temp<br />
;<br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116; <br />
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;<br />
<br />
tempb = 1.0 / ( 0.1 * GAMMA * xd[1]);<br />
tempb2 = tempb * tempb;<br />
rhs[2] = ( C01 <br />
+ C11 * tempb <br />
+ C21 * tempb2<br />
+ C31 * tempb2 * tempb<br />
+ C41 * tempb2 * tempb2<br />
+ C51 * tempb2 * tempb2 * tempb )<br />
; <br />
if (xd[1] < V2 + 1e-2) {<br />
rhs[1] = G * E * A2 / WEFF / GAMMA;<br />
rhs[2] = P2;<br />
}<br />
<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Parallel */<br />
static void ffcn2b( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = G * E * A3 / WEFF / GAMMA;<br />
rhs[2] = P3;<br />
}<br />
<br />
static void ffcn2c( double *xd, double *rhs )<br />
{<br />
double T, R, temp, temp2, tempb, tempb2;<br />
<br />
rhs[0] = xd[1];<br />
<br />
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 1.0);<br />
temp2 = temp * temp;<br />
T = B02<br />
+ B12 * temp <br />
+ B22 * temp2<br />
+ B32 * temp2 * temp<br />
+ B42 * temp2 * temp2<br />
+ B52 * temp2 * temp2 * temp<br />
; <br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;<br />
<br />
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;<br />
<br />
tempb = 1.0 / ( 0.1 * GAMMA * xd[1] - 1);<br />
tempb2 = tempb * tempb;<br />
rhs[2] = ( C02 <br />
+ C12 * tempb <br />
+ C22 * tempb2<br />
+ C32 * tempb2 * tempb<br />
+ C42 * tempb2 * tempb2<br />
+ C52 * tempb2 * tempb2 * tempb )<br />
; <br />
<br />
<br />
// Achtung<br />
if (xd[1] < V3 + 1e-2) {<br />
rhs[1] = G * E * A3 / WEFF / GAMMA;<br />
rhs[2] = P3;<br />
}<br />
<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Rollen */<br />
static void ffcn3( double *xd, double *rhs )<br />
{<br />
double R;<br />
<br />
rhs[0] = xd[1];<br />
<br />
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;<br />
rhs[1] = ( - G * R / WEFF - CC / GAMMA ) ;<br />
<br />
rhs[2] = 0;<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Bremsen */<br />
static void ffcn4( double *xd, double *rhs )<br />
{<br />
rhs[0] = xd[1];<br />
rhs[1] = - 3.0 / GAMMA;<br />
rhs[2] = 0;<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for v < V1*/<br />
static void ffcn_A(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1a( xd, rhs0 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V1 < v < V2*/<br />
static void ffcn_B(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1b( xd, rhs0 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V2 < v < V3*/<br />
static void ffcn_C(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1c( xd, rhs0 );<br />
ffcn2b( xd, rhs1 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
/* -------------------------------------------- */<br />
/* Overall right hand side for V3 < v */<br />
static void ffcn_D(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn1c( xd, rhs0 );<br />
ffcn2c( xd, rhs1 );<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
<br />
/* -------------------------------------------- */<br />
/* Pure coasting resp. braking on last stage to avoid nondifferentiabilities*/<br />
static void ffcn_E(double *t, double *xd, double *xa, double *u,<br />
double *p, double *rhs, double *rwh, long *iwh, long *info)<br />
{<br />
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];<br />
long i;<br />
<br />
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }<br />
<br />
ffcn3( xd, rhs2 );<br />
ffcn4( xd, rhs3 );<br />
<br />
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];<br />
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];<br />
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];<br />
}<br />
<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem&diff=1560Lotka Volterra fishing problem2016-01-28T09:21:04Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nw = 1<br />
|nre = 3<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->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. <br />
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.<br />
<br />
<br />
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.<br />
<br />
The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0(t) & = & x_0(t) - x_0(t) x_1(t) - \; c_0 x_0(t) \; w(t), \\<br />
& \dot{x}_1(t) & = & - x_1(t) + x_0(t) x_1(t) - \; c_1 x_1(t) \; w(t), \\<br />
& \dot{x}_2(t) & = & (x_0(t) - 1)^2 + (x_1(t) - 1)^2, \\[1.5ex]<br />
& x(0) &=& (0.5, 0.7, 0)^T, \\<br />
& w(t) &\in& \{0, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
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>.<br />
<br />
== Parameters ==<br />
<br />
These fixed values are used within the model.<br />
<br />
<math><br />
\begin{array}{rcl}<br />
[t_0, t_f] &=& [0, 12],\\<br />
(c_0, c_1) &=& (0.4, 0.2).<br />
\end{array}<br />
</math><br />
<br />
== Reference Solutions ==<br />
<br />
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.<br />
<br />
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>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.<br />
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.<br />
Image:lotka2Switches.png| Control and differential states with only two switches.<br />
Image:lotka56Switches.png| Control and differential states with 56 switches.<br />
</gallery><br />
<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]<br />
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]<br />
* [[:Category:optimica | optimica code]] at [[Lotka Volterra fishing problem (optimica)]]<br />
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]<br />
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]<br />
* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]<br />
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]<br />
* [[:Category:GloOptCon | GloOptCon code]] at [[Lotka Volterra fishing problem (GloOptCon)]]<br />
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]<br />
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]<br />
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]<br />
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]<br />
<br />
== Variants ==<br />
<br />
There are several alternative formulations and variants of the above problem, in particular<br />
<br />
* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],<br />
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,<br />
* 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>,<br />
* different fishing control functions for the two species,<br />
* different parameters and start values.<br />
<br />
== Miscellaneous and Further Reading ==<br />
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.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Tracking objective]]<br />
[[Category:Chattering]]<br />
[[Category:Sensitivity-seeking arcs]]<br />
[[Category:Population dynamics]]<br />
<br />
<br />
<!--Testing Graphviz<br />
<graphviz border='frame' format='svg'><br />
digraph G {Hello->World!}<br />
</graphviz><br />
--></div>SebastianSagerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem_(Muscod)&diff=1559Lotka Volterra fishing problem (Muscod)2016-01-28T09:20:34Z<p>SebastianSager: SebastianSager moved page Lotka Volterra fishing problem (C) to Lotka Volterra fishing problem (Muscod) without leaving a redirect: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Lotka Volterra fishing problem]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
/* steady state with u == 0 */<br />
double ref0 = 1, ref1 = 1;<br />
<br />
/* Biomass of prey */<br />
rhs[0] = xd[0] - xd[0]*xd[1] - p[0]*u[0]*xd[0];<br />
/* Biomass of predator */<br />
rhs[1] = - xd[1] + xd[0]*xd[1] - p[1]*u[0]*xd[1];<br />
/* Deviation from reference trajectory */<br />
rhs[2] = (xd[0]-ref0)*(xd[0]-ref0) + (xd[1]-ref1)*(xd[1]-ref1);<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=F-8_aircraft_(Muscod)&diff=1558F-8 aircraft (Muscod)2016-01-28T09:20:05Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[F-8 aircraft]] problem in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
double x1 = xd[0];<br />
double x2 = xd[1];<br />
double x3 = xd[2];<br />
<br />
#ifdef CONVEXIFIED<br />
double xi = 0.05236;<br />
<br />
rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2 -x1*x1*x3 + 3.846*x1*x1*x1<br />
+ 0.215*xi - 0.28*x1*x1*xi + 0.47*x1*xi*xi - 0.63*xi*xi*xi<br />
- 2*u[0] * ( 0.215*xi - 0.28*x1*x1*xi - 0.63*xi*xi*xi );<br />
rhs[1] = x3;<br />
rhs[2] = - 4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1<br />
+ 20.967*xi - 6.265*x1*x1*xi + 46*x1*xi*xi - 61.4*xi*xi*xi<br />
- 2*u[0] * ( 20.967*xi - 6.265*x1*x1*xi - 61.4*xi*xi*xi );<br />
#else<br />
double u_ = -0.05236 + u[0]*2.0*0.05236;<br />
<br />
rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2<br />
-x1*x1*x3 + 3.846*x1*x1*x1 - 0.215*u_ + 0.28*x1*x1*u_ + 0.47*x1*u_*u_ + 0.63*u_*u_*u_;<br />
rhs[1] = x3;<br />
rhs[2] = -4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1<br />
- 20.967*u_ + 6.265*x1*x1*u_ + 46*x1*u_*u_ + 61.4*u_*u_*u_;<br />
#endif<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Annihilation_of_calcium_oscillations_(Muscod)&diff=1557Annihilation of calcium oscillations (Muscod)2016-01-28T09:19:48Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for [[Annihilation of calcium oscillations]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
rhs[0] = p[0]<br />
+ p[1]*xd[0]<br />
- p[2]*xd[0]*xd[1]/(xd[0]+p[3])<br />
- p[4]*xd[0]*xd[2]/(xd[0]+p[5]);<br />
<br />
rhs[1] = p[6]*xd[0]<br />
- p[7]*xd[1]/(xd[1]+p[8]);<br />
<br />
rhs[2] = p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])<br />
+ p[18]*p[11]*xd[1]<br />
+ p[12]*xd[0]<br />
- p[13]*xd[2]/(u[0]*xd[2]+p[14])<br />
- p[15]*xd[2]/(xd[2]+p[16])<br />
+ 0.1*xd[3];<br />
<br />
rhs[3] = -p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])<br />
+ p[15]*xd[2]/(xd[2]+p[16])<br />
- 0.1*xd[3];<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Fuller%27s_Problem_(Muscod)&diff=1556Fuller's Problem (Muscod)2016-01-28T09:19:32Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for [[Fuller's problem]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
double myu = -1 + 2*u[0];<br />
<br />
rhs[0] = xd[1];<br />
rhs[1] = myu;<br />
rhs[2] = xd[0]*xd[0];<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem_(Muscod)&diff=1555Lotka Volterra fishing problem (Muscod)2016-01-28T09:18:56Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for the [[Lotka Volterra fishing problem]] in [[:Category:Muscod | Muscod code]] read as follows<br />
<br />
<source lang="cpp"><br />
/* steady state with u == 0 */<br />
double ref0 = 1, ref1 = 1;<br />
<br />
/* Biomass of prey */<br />
rhs[0] = xd[0] - xd[0]*xd[1] - p[0]*u[0]*xd[0];<br />
/* Biomass of predator */<br />
rhs[1] = - xd[1] + xd[0]*xd[1] - p[1]*u[0]*xd[1];<br />
/* Deviation from reference trajectory */<br />
rhs[2] = (xd[0]-ref0)*(xd[0]-ref0) + (xd[1]-ref1)*(xd[1]-ref1);<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Fuller%27s_problem&diff=1554Fuller's problem2016-01-28T09:18:30Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nw = 1<br />
|nre = 4<br />
}}<br />
<br />
The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bib id="Fuller1963" />. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.<br />
<br />
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix initial and terminal values of the differential states.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [t_0, t_f]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\<br />
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]<br />
& x(0) &=& x_S, \\<br />
& x(t_f) &=& x_T, \\<br />
& w(t) &\in& \{0, 1\}.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
We use <math>x_S = x_T = (0.01, 0)^T</math>.<br />
<br />
== Reference Solutions ==<br />
<br />
===Solutions obtained with optimica===<br />
<br />
The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.<br />
[[File:Fullerspng.png|left|200px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]<br />
<br /><br />
<br style="clear:both;"> <br />
<br />
<!--<br />
The optimal trajectory for the relaxed control problem on an equidistant grid $\mathcal{G}^0$ with $n_\mathrm{ms} = 19$ is shown in the top row of figure \ref{GCSEX3}. Note that this solution is not bang--bang due to the discretization of the control space. Even if this discretization is made very fine, a trajectory with $w(t) = 0.5$ on an interval in the middle of $[0,1]$ will be found as a minimum. <br />
<br />
The switching time optimization method yields an objective value of $1.89 \; 10^{-4}$. As the objective function value of the relaxed problem is smaller, $\Phi^{\mathrm{RL}} = 1.53 \; 10^{-5}$, one might want to reduce the function value further, e.g. closer than $\varepsilon = 10^{-6}$ to $\Phi^{\mathrm{RL}}$. If we apply algorithm \ref{AlgMIOC}, we obtain the trajectory shown in the bottommost row of figure \ref{GCSEX3} that yields an objective function value of $1.52 \; 10^{-5}$ and switches $35$ times.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.<br />
</gallery><br />
--><br />
<br />
== Source Code ==<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[Fuller's Problem (Muscod)]]<br />
* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]<br />
<br />
== Miscellaneous and further reading ==<br />
<br />
An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bib id="Zelikin1994" />, a recent investigation of chattering controls in relay feedback systems in <bib id="Johansson2002" />.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Tracking objective]]<br />
[[Category:Chattering]]<br />
[[Category:Bang bang]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Fuller%27s_Problem_(Muscod)&diff=1553Fuller's Problem (Muscod)2016-01-28T09:18:05Z<p>SebastianSager: SebastianSager moved page Fuller's Problem (C) to Fuller's Problem (Muscod) without leaving a redirect: Moved C to Muscod</p>
<hr />
<div>The differential equations for [[Fuller's problem]] in [[:Category:C | C code]] read as follows<br />
<br />
<source lang="cpp"><br />
double myu = -1 + 2*u[0];<br />
<br />
rhs[0] = xd[1];<br />
rhs[1] = myu;<br />
rhs[2] = xd[0]*xd[0];<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=Fuller%27s_Problem_(Muscod)&diff=1552Fuller's Problem (Muscod)2016-01-28T09:17:53Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>The differential equations for [[Fuller's problem]] in [[:Category:C | C code]] read as follows<br />
<br />
<source lang="cpp"><br />
double myu = -1 + 2*u[0];<br />
<br />
rhs[0] = xd[1];<br />
rhs[1] = myu;<br />
rhs[2] = xd[0]*xd[0];<br />
</source><br />
<br />
[[Category:Muscod]]</div>SebastianSagerhttps://mintoc.de/index.php?title=F-8_aircraft&diff=1551F-8 aircraft2016-01-28T09:17:15Z<p>SebastianSager: Moved C to Muscod</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nw = 1<br />
|nre = 6<br />
}}<br />
<br />
The F-8 aircraft control problem is based on a very simple aircraft model. The control problem was introduced by Kaya and Noakes and aims at controlling an aircraft in a time-optimal way from an initial state to a terminal state.<br />
<br />
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix both initial and terminal values of the differential states.<br />
<br />
The optimal, relaxed control function shows [[:Category:Bang bang|bang bang]] behavior. The problem is furthermore interesting as it should be reformulated equivalently.<br />
<br />
== Mathematical formulation ==<br />
<br />
For <math>t \in [0, T]</math> almost everywhere the mixed-integer optimal control problem is given by<br />
<br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w, T} & T \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0 &=& -0.877 \; x_0 + x_2 - 0.088 \; x_0 \; x_2 + 0.47 \; x_0^2 - 0.019 \; x_1^2 - x_0^2 \; x_2 + 3.846 \; x_0^3 \\<br />
&&& - 0.215 \; w + 0.28 \; x_0^2 \; w + 0.47 \; x_0 \; w^2 + 0.63 \; w^3 \\ <br />
& \dot{x}_1 &=& x_2 \\<br />
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\<br />
&&& - 20.967 \; w + 6.265 \; x_0^2 \; w + 46 \; x_0 \; w^2 + 61.4 \; w^3 \\ <br />
& x(0) &=& (0.4655,0,0)^T, \\<br />
& x(T) &=& (0,0,0)^T, \\<br />
& w(t) &\in& \{-0.05236,0.05236\}.<br />
\end{array} <br />
</math><br />
<br />
<math>x_0</math> is the angle of attack in radians, <math>x_1</math> is the pitch angle, <math>x_2</math> is the pitch rate in rad/s, and the control function <math>w = w(t)</math> is the tail deflection angle in radians. This model goes back to Garrard<bib id="Garrard1977" />.<br />
<br />
In the control problem, both initial and terminal values of the differential states are fixed.<br />
<br />
== Reformulation ==<br />
<br />
The control w(t) is restricted to take values from a finite set only. Hence, the control problem can be reformulated equivalently to<br />
<br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, w, T} & T \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_0 &=& -0.877 \; x_0 + x_2 - 0.088 \; x_0 \; x_2 + 0.47 \; x_0^2 - 0.019 \; x_1^2 - x_0^2 \; x_2 + 3.846 \; x_0^3 \\<br />
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \right) \; w \\ <br />
&&& - \left( - 0.215 \; \xi + 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 + 0.63 \; \xi^3 \right) \; (1 - w) \\ <br />
& &=& -0.877 \; x_0 + x_2 - 0.088 \; x_0 \; x_2 + 0.47 \; x_0^2 - 0.019 \; x_1^2 - x_0^2 \; x_2 + 3.846 \; x_0^3 \\<br />
&&& + 0.215 \; \xi - 0.28 \; x_0^2 \; \xi + 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \\ <br />
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.63 \; \xi^3 \right) \; 2 w \\ <br />
& \dot{x}_1 &=& x_2 \\<br />
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\<br />
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \right) \; w \\ <br />
&&& - \left( - 20.967 \; \xi + 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 + 61.4 \; \xi^3 \right) \; (1 - w) \\ <br />
& &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\<br />
&&& + 20.967 \; \xi - 6.265 \; x_0^2 \; \xi + 46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \\ <br />
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi - 61.4 \; \xi^3 \right) \; 2 w \\ <br />
& x(0) &=& (0.4655,0,0)^T, \\<br />
& x(T) &=& (0,0,0)^T, \\<br />
& w(t) &\in& \{0,1\},<br />
\end{array} <br />
</math><br />
<br />
with <math>\xi = 0.05236</math>. Note that there is a bijection between optimal solutions of the two problems.<br />
<br />
== Reference solutions ==<br />
<br />
We provide here a comparison of different solutions reported in the literature. The numbers show the respective lengths <math>t_i - t_{i-1}</math> of the switching arcs with the value of <math>w(t)</math> on the upper or lower bound (given in the second column). ''Claim'' denotes what is stated in the respective publication, ''Simulation'' shows values obtained by a simulation with a Runge-Kutta-Fehlberg method of 4th/5th order and an integration tolerance of <math>10^{-8}</math>.<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Arc !! w(t) !! Lee et al.<bib id="Lee1997a" /> !! Kaya<bib id="Kaya2003" /> !! Sager<bib id="Sager2005" /> !! [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] !! Sager<br />
|- <br />
| align=center | 1 || align=center | 1 || 0.00000 || 0.10292 || 0.10235 || 0.0 || 1.13492<br />
|-<br />
| align=center | 2 || align=center | 0 || 2.18800 || 1.92793 || 1.92812 || 0.608750 || 0.34703<br />
|-<br />
| align=center | 3 || align=center | 1 || 0.16400 || 0.16687 || 0.16645 || 3.136514 || 1.60721<br />
|-<br />
| align=center | 4 || align=center | 0 || 2.88100 || 2.74338 || 2.73071 || 0.654550 || 0.69169<br />
|-<br />
| align=center | 5 || align=center | 1 || 0.33000 || 0.32992 || 0.32994 || 0.0 || 0.0<br />
|-<br />
| align=center | 6 || align=center | 0 || 0.47200 || 0.47116 || 0.47107 || 0.0 || 0.0<br />
|-<br />
| Claim: Infeasibility || align=center | - || 1.00E-10 || 7.30E-06 || 5.90E-06 || 3.29169e-06 || 2.21723e-07<br />
|-<br />
| Claim: Objective || align=center | - || 6.03500 || 5.74217 || 5.72864 || 4.39981 || 3.78086<br />
|- style="font-style:italic;background-color:#ffffcc;"<br />
! Simulation: Infeasibility || align=center | - || 1.75E-03 || 1.64E-03 || 5.90E-06 || 3.29169e-06 || 2.21723e-07<br />
|- style="font-style:italic;background-color:#ffffcc;"<br />
! Simulation: Objective || align=center | - || 6.03500 || 5.74218 || 5.72864 || 4.39981 || 3.78086<br />
|}<br />
<br />
The best known optimal objective value of this problem given is given by <math>T = 3.78086</math>. The corresponding solution is shown in the rightmost plot. The solution of bang-bang type switches three times, starting with <math>w(t) = 1</math>.<br />
<br />
<gallery caption="Reference solution plots" widths="240px" heights="167px" perrow="3"><br />
Image:f8aircraftRelaxedCoarse.png| Locally optimal relaxed control on a coarse control discretization grid and corresponding differential states.<br />
Image:f8aircraftRelaxedFine.png| Locally optimal relaxed control on a fine, adaptively chosen control discretization grid and corresponding differential states.<br />
Image:f8aircraftInteger.png| Integer control and corresponding differential states from Sager<bib id="Sager2005" />.<br />
<br />
Image:f8aircraftSchlueter.png| Integer control and corresponding differential states from [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] solution.<br />
Image:f8aircraftSager2009.png| Best known integer control and corresponding differential states from Sager solution.<br />
</gallery><br />
<br />
=== Optimica ===<br />
<br />
Objective : 5.12799232<br />
infeasibility : 6.2235588037251599e-10<br />
<br />
<gallery caption="Obtained solution plots" widths="240px" heights="167px" perrow="3"><br />
Image:f8aircraft1.png| (Probably sub-)Optimal control.<br />
Image:f8aircraft2.png| Corresponding differential states.<br />
<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL]] at [[F-8 aircraft (AMPL)]]<br />
* [[:Category:Muscod | Muscod code]] at [[F-8 aircraft (Muscod)]]<br />
* [[:Category:optimica | optimica]] at [[F-8 aircraft (optimica)]]<br />
<br />
== Variants ==<br />
<br />
* a prescribed time grid for the control function, see [[F-8 aircraft (AMPL)]],<br />
<br />
== Miscellaneous and Further Reading ==<br />
See <bib id="Kaya2003" /> and <bib id="Sager2005" /> for details.<br />
<br />
== References ==<br />
<biblist /><br />
<br />
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here --><br />
[[Category:Bang bang]]<br />
[[Category:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Aeronautics]]<br />
[[Category:Outer convexification]]<br />
[[Category: Minimum time]]</div>SebastianSager