https://mintoc.de/api.php?action=feedcontributions&user=FelixMueller&feedformat=atommintOC - User contributions [en]2024-03-28T09:37:25ZUser contributionsMediaWiki 1.25.2https://mintoc.de/index.php?title=Time_optimal_car_problem&diff=2148Time optimal car problem2016-09-29T19:22:57Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 6<br />
|nre = 4<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Time optimal car problem "consists of starting and stopping a car in minimum for a fixed distance (300 units)" and can be found e.g. in <bib id="Cuthrell1987" /> and <bib id="Logsdon1992" />.<br />
<br />
== Mathematical formulation ==<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{z, u, t_f} & t_f \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{z}_1 & = & z_2, \\<br />
& \dot{z}_2 & = & u, \\[1.5ex]<br />
& z(0) &=& (0,0)^T, \\<br />
& z(t_f) &=& (300,0)^T, \\[1.5ex]<br />
& z_1(t) & \in & [0,33],\\<br />
& z_2(t) & \in & [0,330],\\<br />
& u(t) & \in & [-2,1].\\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
where <math> z = (z_1, z_2)^T </math>.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Time optimal car problem (TACO)]]<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Minimum time]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Isomerization_of_Alpha-Pinene_problem&diff=2147Isomerization of Alpha-Pinene problem2016-09-29T19:22:33Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 5<br />
|np = 5<br />
|nc = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Isomerization of Alpha-Pinene problem tries to determine "reaction coefficients in the thermal isometrization of <math> \alpha </math>-Pinene." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{8} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_2) y_1, \\<br />
& \dot{y}_2 & = & \theta_1 y_1, \\<br />
& \dot{y}_3 & = & \theta_2 y_1 - (\theta_3 + \theta_4) y_3 + \theta_5 y_5, \\<br />
& \dot{y}_4 & = & \theta_3 y_3, \\<br />
& \dot{y}_5 & = & \theta_4 y_3 - \theta_5 y_5, \\<br />
& \theta_i & \geq & 0 \quad i = 1,...,5.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_8 </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Isomerization of Alpha-Pinene problem (TACO)]]<br />
<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:DAE model]]<br />
[[Category:Chemical engineering]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Catalytic_cracking_problem&diff=2146Catalytic cracking problem2016-09-29T19:21:15Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 2<br />
|np = 3<br />
|nc = 3<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalytic cracking problem tries to determine "reaction coefficients for the catalytic cracking of gas oil into gas and other byproducts." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{21} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_3) y_1^2, \\<br />
& \dot{y}_2 & = & \theta_1 y_1^2 - \theta_2 y_2, \\<br />
& \theta_i & \geq & 0 \quad i = 1,...,3.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_{21} </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalytic cracking problem (TACO)]]<br />
<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:DAE model]]<br />
[[Category:Chemical engineering]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Catalyst_mixing_problem&diff=2145Catalyst mixing problem2016-09-29T19:21:02Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 2<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalyst mixing problem seeks an optimal policy for mixing two catalysts "along the length of a tubular plug <br />
ow reactor involving several reactions". (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} &-1 + x_1(t_f) + x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_1 & = & u ( 10 x_2 - x_1), \\<br />
& \dot{x}_2 & = & u ( x_1 - 10 x_2) - (1 - u \, x_2) , \\<br />
& x(t_0) &=& (1, 0)^T, \\<br />
& u(t) &\in& [0,1].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
In this model the parameters used are <math> t_0 = 0, \, \, t_f = 1 </math>.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalyst mixing problem (TACO)]]<br />
<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:Chemical engineering]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Category:Julia/JuMP&diff=2144Category:Julia/JuMP2016-08-11T09:57:52Z<p>FelixMueller: </p>
<hr />
<div>This category contains the list of problems for which JuMP code has been added.<br />
<br />
JuMP stands for 'Julia for Mathematical Programming' and is a modelling language similar to [[:Category:AMPL | AMPL]]. It is written in the programming language Julia. <br />
<br />
Julia/JuMP is an open source project started in 2012 at MIT, Cambridge. Its aim is to be highly intuitive and as quick as commercial modelling languages while still maintaining the open source aspect of programming. As JuMP (like AMPL) is only able to solve Mixed-Integer Nonlinear Optimization Problems (MINLPs) you cannot simply pass an MIOCP to JuMP but have to somehow reduce it to a MINLP e.g. by discretizing by hand, which can then be solved by JuMP.<br />
<br />
For more information on Julia see the [http://julialang.org/ Julia homepage], the [http://www.juliaopt.org/ homepage of the Julia optimization group] or the [http://julia.readthedocs.org/en/latest/ documentation page of Julia].<br />
<br />
More info regarding JuMP can be found on [https://jump.readthedocs.org/en/latest/ the documentation page of JuMP].<br />
<br />
[[Artificial Pancreas Experimental Design]]<br />
<br />
<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: Implementation]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Time_optimal_car_problem&diff=2141Time optimal car problem2016-07-31T07:24:56Z<p>FelixMueller: Created page with "{{Dimensions |nd = 1 |nx = 2 |nu = 1 |nc = 6 |nre = 4 }}<!-- Do not insert line break here or Dimensions Box moves up in the layou..."</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 6<br />
|nre = 4<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Time optimal car problem "consists of starting and stopping a car in minimum for a fixed distance (300 units)" and can be found e.g. in <bib id="Cuthrell1987" /> and <bib id="Logsdon1992" />.<br />
<br />
== Mathematical formulation ==<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{z, u, t_f} & t_f \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{z}_1 & = & z_2, \\<br />
& \dot{z}_2 & = & u, \\[1.5ex]<br />
& z(0) &=& (0,0)^T, \\<br />
& z(t_f) &=& (300,0)^T, \\[1.5ex]<br />
& z_1(t) & \in & [0,33],\\<br />
& z_2(t) & \in & [0,330],\\<br />
& u(t) & \in & [-2,1].\\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
where <math> z = (z_1, z_2)^T </math>.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Time optimal car problem (TACO)]]<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Minimum time]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Robot_arm_problem&diff=2140Robot arm problem2016-07-31T07:15:48Z<p>FelixMueller: Created page with "{{Dimensions |nd = 1 |nx = 3 |nu = 3 |nc = 12 |nre = 12 }}<!-- Do not insert line break here or Dimensions Box moves up in the lay..."</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nu = 3<br />
|nc = 12<br />
|nre = 12<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The robot arm problem focuses on minimizing the time used by a robot arm to move from an origin to a destination.<br />
The arm is a bar of length <math> L </math> and sticks out distance <math> \rho </math> from its moving axis, while sticking out distance <math> L - \rho </math> in the other direction. The problem can be found in <bib id="Moessner1995" /> or in the [http://www.mcs.anl.gov/~more/cops/ COPS library].<br />
<br />
== Model formulation ==<br />
<br />
The problem is set up using the length <math> \rho </math>, "the vertical angles <math> (\theta, \Phi) </math> from the horizontal plane, the controls <math> u=(u_{\rho},u_{\theta},u_{\Phi}) </math> and the final time <math> t_f </math>".<br />
<br />
The moving robot is modelled with the following equations:<br />
<br />
<math> \ddot{\rho} = \frac{u_{\rho}}{L}, \qquad \ddot{\theta} = \frac{u_{\theta}}{I_{\theta}}, \qquad \ddot{\Phi} = \frac{u_{\Phi}}{I_{\Phi}}</math><br />
<br />
where <math> I </math> characterizes the moment of inertia, i.e.<br />
<br />
<math> <br />
\begin{array}{ccl}<br />
I_{\theta} & = & \frac{((L-\rho)^3 + \rho^3)}{3} \cdot \sin(\Phi)^2, \\<br />
I_{\Phi} & = & \frac{((L-\rho)^3 + \rho^3)}{3}.<br />
\end{array}<br />
</math><br />
<br />
The path constraints on the states <math> x= (\rho, \theta, \Phi) </math> and on the controls <math> u = (u_{\rho},u_{\theta},u_{\Phi}) </math> as well as the boundary conditions can be seen in the optimization problem further down.<br />
<br />
== Optimization problem ==<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{x, u, t_f} & t_f \\[1.5ex]<br />
\mbox{s.t.} <br />
& \ddot{\rho} & = & \frac{u_{\rho}}{L}, \\<br />
& \ddot{\theta} & = & \frac{u_{\theta}}{I_{\theta}}, \\<br />
& \ddot{\Phi} & = & \frac{u_{\Phi}}{I_{\Phi}}, \\[1.5ex]<br />
& x(0) &=& (4.5, 0, \frac{\pi}{4})^T, \\<br />
& x(t_f) &=& (4.5, \frac{2\pi}{3}, \frac{\pi}{4})^T, \\<br />
& \dot{x}(0) &=& (0,0,0)^T, \\<br />
& \dot{x}(t_f) &=& (0,0,0)^T, \\[1.5ex]<br />
& \rho(t) & \in & [0,L],\\<br />
& \theta(t) & \in & [-\pi, \pi],\\<br />
& \Phi(t) & \in & [0, \pi],\\<br />
& u_{\rho} & \leq & 1,\\<br />
& u_{\theta} & \leq & 1,\\<br />
& u_{\Phi} & \leq & 1.\\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
where <math> I </math> is the moment of inertia as above.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Robot arm problem (TACO)]]<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Minimum time]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Lotka_Volterra_fishing_problem&diff=2139Lotka Volterra fishing problem2016-07-31T06:57:06Z<p>FelixMueller: /* Mathematical formulation */</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 />
The mixed-integer optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\<br />
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\<br />
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 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: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>FelixMuellerhttps://mintoc.de/index.php?title=Category:Julia/JuMP&diff=2138Category:Julia/JuMP2016-07-31T06:45:37Z<p>FelixMueller: </p>
<hr />
<div>This category contains the list of problems for which JuMP code has been added.<br />
<br />
JuMP stands for 'Julia for Mathematical Programming' and is a modelling language similar to [[:Category:AMPL | AMPL]]. It is written in the programming language Julia. <br />
<br />
Julia/JuMP is an open source project started in 2012 at MIT, Cambridge. Its aim is to be highly intuitive and as quick as commercial modelling languages while still maintaining the open source aspect of programming. As JuMP (like AMPL) is only able to solve Mixed-Integer Nonlinear Optimization Problems (MINLPs) you cannot simply pass an MIOCP to JuMP but have to somehow reduce it to a MINLP e.g. by discretizing by hand, which can then be solved by JuMP.<br />
<br />
For more information on Julia see the [http://julialang.org/ Julia homepage], the [http://www.juliaopt.org/ homepage of the Julia optimization group] or the [http://julia.readthedocs.org/en/latest/ documentation page of Julia].<br />
<br />
More info regarding JuMP can be found on [https://jump.readthedocs.org/en/latest/ the documentation page of JuMP].<br />
<br />
<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: Implementation]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Category:Julia/JuMP&diff=2137Category:Julia/JuMP2016-07-29T15:58:15Z<p>FelixMueller: </p>
<hr />
<div>This category contains the list of problems for which JuMP code has been added.<br />
<br />
JuMP stands for 'Julia for Mathematical Programming' and is a modelling language similar to [[:Category:AMPL | AMPL]]. It is written in the programming language Julia. <br />
<br />
Julia/JuMP is an open source project started in 2012 at MIT, Cambridge. Its aim is to be highly intuitive and as quick as commercial modelling languages while still maintaining the open source aspect of programming. As JuMP (like AMPL) is only able to solve Mixed-Integer Nonlinear Optimization Problems (MINLPs) you cannot simply pass an MIOCP to JuMP but have to somehow reduce it to a MINLP e.g. by discretizing by hand, which can then be solved by JuMP.<br />
<br />
For more information on Julia see the [http://julialang.org/ Julia homepage], the [http://www.juliaopt.org/ homepage of the Julia optimization group] or the [http://julia.readthedocs.org/en/latest/ documentation page of Julia].<br />
<br />
More info regarding JuMP can be found on [https://jump.readthedocs.org/en/latest/ the documentation page of JuMP].<br />
<br />
Testing: [[Artificial pancreas Experimental Design]]<br />
<br />
<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: Implementation]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Diels-Alder_Reaction_Experimental_Design&diff=2136Diels-Alder Reaction Experimental Design2016-07-29T14:17:49Z<p>FelixMueller: /* Source Code */</p>
<hr />
<div>The '''Diels-Alder Reaction''' is an organic chemical reaction. <br />
A conjugated diene and a substituted alkene react and form a substituted cyclohexene system.<br />
Stefan Körkel used this model in his PhD thesis to compute optimum experimental designs.<br />
<br />
<br />
== Model Formulation ==<br />
<br />
The reactionkinetics can be modelled by the following differential equation system:<br />
<br />
<math><br />
\begin{array}{rcl}<br />
\dot{n_1}(t) &=& -k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}}, \\<br />
& & \\<br />
\dot{n_2}(t) &=& -k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}}, \\<br />
& & \\<br />
\dot{n_3}(t) &=& \ \ k \cdot \frac{n_1(t) \ \cdot \ n_2(t)}{m_{tot}} \\<br />
& & \\<br />
\dot{n_4}(t) &=& 0<br />
\end{array} <br />
</math><br />
<br />
The reaction velocity constant <math>k</math> consists of two parts. One part reflects the non-catalytic and the other the catalytic reaction. The velocity law follows the Arrhenius relation<br />
<br />
<math><br />
k = k_1 \ \cdot \ exp(- \frac{E_1}{R} \ \cdot \ (\frac{1}{T(t)} \ - \ \frac{1}{T_{ref}}) ) \ + \ k_{cat} \ \cdot \ c_{cat} \ \cdot \ exp(-\lambda \ \cdot \ t) \ \cdot \ exp( - \frac{E_{cat}}{R} \ \cdot \ (\ \frac{1}{T(t)} \ - \ \frac{1}{T_{ref}}) )<br />
</math><br />
<br />
Total mass: <br />
<br />
<math><br />
m_{tot} = n_1 \ \cdot \ M_1 \ + \ n_2 \ \cdot \ M_2 \ + \ n_3 \ \cdot \ M_3 \ + \ n_4 \ \cdot \ M_4 <br />
</math><br />
<br />
Temperature in Kelvin:<br />
<br />
<math><br />
T(t) = \vartheta (t) + 273<br />
</math><br />
<br />
The ODE system is summarized to:<br />
<br />
<math><br />
\begin{array}{rcl}<br />
\dot{x}(t) &=& f(x(t), u(t), p) <br />
\end{array} <br />
</math><br />
<br />
== Constraints ==<br />
<br />
The control variables are constrained with respect to the mass of sample weights (initial mass):<br />
<br />
<br />
<math><br />
\begin{array}{cll}<br />
0.1 & \le & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \le 10 <br />
\end{array} <br />
</math><br />
<br />
and to the mass of active ingredient content (fraction of active substances):<br />
<br />
<p><br />
<math><br />
\begin{array}{cll}<br />
0.1 & \le & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \le 0.7<br />
\end{array} <br />
</math><br />
<br />
== Optimum Experimental Design Problem ==<br />
<br />
The aim is to compute an optimal experimental design <math>\xi = (q,w)</math> which minimizes the uncertainties of the parameters <math>k_1, k_{cat}, E_1, E_{cat}, \lambda</math>. So, we have to solve the following optimum experimental design problem:<br />
<br />
<br />
<math><br />
\begin{array}{cll}<br />
\displaystyle \min_{x^i,\ G^i,\ F^i,\ Tc^i,\ n_{a1}^i,\ n_{a2}^i,\ n_{a4}^i,\ c_{kat}^i,\ \vartheta(t)^i} && trace(F^{-1} (t_{end})) \\[1.5ex]<br />
\mbox{s.t.} \\<br />
\dot{x}^i(t) & = & f(x^i(t), u^i(t),p), \\<br />
\\<br />
h(t) & = & \frac{n_3(t) \ \cdot \ M_3}{m_{tot}} \ \cdot \ 100 \\<br />
\\<br />
\dot{G}^i(t) & = & f_x(x^i(t),u^i(t),p)G^i(t) \ + \ f_p(x^i(t),u^i(t),p) \\<br />
\\<br />
\dot{F}(t) & = & \sum\limits_{i=1}^{4} w^i(t) (h^i_x(x^i(t),u^i(t),p)G^i(t))^T (h^i_x(x^i(t),u(t),p)G^i(t)) \\<br />
\\<br />
0.1 & \le & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \\<br />
\\<br />
10 & \ge & n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 \\<br />
\\<br />
0.1 & \le & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \\<br />
\\<br />
0.7 & \ge & \frac{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 }{ n_{a1} \ \cdot \ M_1 \ + \ n_{a2} \ \cdot \ M_2 \ + \ n_{a4} \ \cdot \ M_4 } \\<br />
\\<br />
\vartheta(t) & = & \left\{ \begin{array}{cl} \vartheta_{lo} + 273 & t \in [t_0,2] \\ <br />
\vartheta_{lo} + \frac{t-2}{6} ( \vartheta_{up} - \vartheta_{lo} ) + 273 & t \in [2,8] \\<br />
\vartheta_{up} + 273 & t \in [8,t_{end}]<br />
\end{array} \right. \\<br />
& & x \in \mathcal{X},\,u \in \mathcal{U},\, p \in P \\<br />
\dot{z}^i(t) & = & w^i(t) \\<br />
z(0) & = & 0 \\<br />
w^i(t) &\in& [0,1] \\<br />
0 & \le & 4 - z^i(t_f). \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
<br />
{| class="wikitable"<br />
|+State variables<br />
|-<br />
|Name<br />
|Symbol<br />
|Initial value (<math>t_0</math>)<br />
|-<br />
|Molar number 1<br />
|<math>n_1(t)</math><br />
|<math>n_1(t_0) = n_{a1} </math><br />
|-<br />
|Molar number 2<br />
|<math>n_2(t)</math><br />
|<math>n_2(t_0) = n_{a2} </math><br />
|-<br />
|Molar number 3<br />
|<math>n_3(t)</math><br />
|<math>n_3(t_0) = 0 </math><br />
|-<br />
|Solvent<br />
|<math>n_4(t)</math><br />
|<math>n_4(t_0) = n_{a4} </math><br />
|}<br />
<br />
<br />
{| class="wikitable"<br />
|+Constants<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|-<br />
|Molar Mass<br />
|<math>M_1</math><br />
|0.1362<br />
|-<br />
|Molar Mass<br />
|<math>M_2</math><br />
|0.09806<br />
|-<br />
|Molar Mass<br />
|<math>M_3</math><br />
|0.23426<br />
|-<br />
|Molar Mass<br />
|<math>M_4</math><br />
|0.236<br />
|-<br />
|Universal gas constant<br />
|<math>R</math><br />
|8.314<br />
|-<br />
|Reference temperature<br />
|<math>T_{ref}</math><br />
|293<br />
|-<br />
|St.dev of measurement error<br />
|<math>\sigma</math><br />
|1<br />
|}<br />
<br />
Remember, in an optimum experimental design problem the parameters of the model are fixed. But, we minimize the parameter's uncertainties by optimizing over the control variables and functions.<br />
<br />
{| class="wikitable"<br />
|+Fixed parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|-<br />
|Steric factor<br />
|<math>k_1</math><br />
|<math>p_1 \cdot 0.01</math><br />
|-<br />
|Steric factor<br />
|<math>k_{kat}</math><br />
|<math>p_2 \cdot 0.10</math><br />
|-<br />
|Activation energie<br />
|<math>E_1</math><br />
|<math>p_3 \cdot 60000</math><br />
|-<br />
|Activation energie<br />
|<math>E_{kat}</math><br />
|<math>p_4 \cdot 40000</math><br />
|-<br />
|Catalyst deactivation coefficient<br />
|<math>\lambda</math><br />
|<math>p_5 \cdot 0.25</math><br />
|}<br />
with <math>p_j = 1, \ j =1, \dots, 5</math><br />
<br />
{| class="wikitable"<br />
|+Optimization/control variables<br />
|-<br />
|Name<br />
|Symbol<br />
|Interval<br />
|Initial value Exp 1<br />
|Initial value Exp 2<br />
|Initial value Exp 3<br />
|Initial value Exp 4<br />
|-<br />
|Initial molar number 1<br />
|<math>n_{a1}</math><br />
|[0,10.0]<br />
|1.0<br />
|1.0<br />
|1.0<br />
|1.0<br />
|-<br />
|Initial molar number 2<br />
|<math>n_{a2}</math><br />
|[0,10.0]<br />
|1.0<br />
|1.0<br />
|1.0<br />
|1.0<br />
|-<br />
|Initial molar number 4<br />
|<math>n_{a4}</math><br />
|[0.4,9.0]<br />
|2.0<br />
|2.0<br />
|2.0<br />
|2.0<br />
|-<br />
|Concentration of the catalyst<br />
|<math>c_{kat}</math><br />
|[0,10.0]<br />
|0.0<br />
|1.0<br />
|2.0<br />
|3.0<br />
|}<br />
<br />
{| class="wikitable"<br />
|+Control function<br />
|-<br />
|Name<br />
|Symbol<br />
|Time interval<br />
|Value interval<br />
|Initial value Exp 1<br />
|Initial value Exp 2<br />
|Initial value Exp 3<br />
|Initial value Exp 4<br />
|-<br />
|Initial molar number 1<br />
|<math>\vartheta(t)</math><br />
|<math>[t_0,2]</math><br />
|[20.0,100.0]<br />
|20.0<br />
|60.0<br />
|40.0<br />
|20.0<br />
|-<br />
|Initial molar number 1<br />
|<math>\vartheta(t)</math><br />
|<math>[2,8]</math><br />
|[20.0,100.0]<br />
|20.0<br />
|60.0<br />
|40.0<br />
|20.0<br />
|-<br />
|Initial molar number 1<br />
|<math>\vartheta(t)</math><br />
|<math>[8,t_{end}]</math><br />
|[20.0,100.0]<br />
|20.0<br />
|60.0<br />
|40.0<br />
|20.0<br />
|}<br />
<br />
'''Measurement grid'''<br />
<br />
<math><br />
\begin{array}{llll}<br />
t_0 = 0 & & & \\<br />
t_{end} = 20 & & & \\<br />
t_j = j/3, & j = 1,\dots, 15, & t_j = j - 10, & j = 16, \dots, 20.<br />
\end{array} <br />
</math><br />
<br />
== Source Code ==<br />
<br />
* The VPLAN code using [[:Category: VPLAN | VPLAN]] can be found in: [[Diels-Alder Reaction Experimental Design (VPLAN)]]<br />
<br />
== References ==<br />
<br />
R. T. Morrison and R.N. Boyd. Organic Chemistry. Allyn and Bacon, Inc., 4th edition, 1983<br />
S. Körkel. Numerische Methoden für Optimale Versuchsplanungsprobleme bei nichtlinearen DAE-Modellen.PhD thesis, Universität Heidelberg, Heidelber,2002<br />
<br />
<br />
[[Category:Optimum Experimental Design]]<br />
[[Category:ODE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2135Industrial robot2016-07-29T13:34:07Z<p>FelixMueller: /* Virtual system */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable, the virtual input, <math> v </math> is integrated twice to produce the virtual state <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
The virtual state and input are constrained by the sets <math> \mathcal{Z} </math> and <math> \mathcal{V} </math> respectively.<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2134Industrial robot2016-07-29T13:33:24Z<p>FelixMueller: /* Virtual system */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable <math> v </math> is integrated twice to produce <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
The virtual state and input are constrained by the sets <math> \mathcal{Z} </math> and <math> \mathcal{V} </math> respectively.<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2133Industrial robot2016-07-29T13:33:05Z<p>FelixMueller: /* Virtual system */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable <math> v </math> is integrated twice to produce <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
The virtual state and input are constrained by the sets <math> \mathcal{Z} \text{ and } \mathcal{V} </math> respectively.<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2132Industrial robot2016-07-29T13:32:52Z<p>FelixMueller: /* Virtual system */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable <math> v </math> is integrated twice to produce <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
The virtual state and input are constrained by the sets <math> \mathcal{Z} \text and } \mathcal{V} </math> respectively.<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2131Industrial robot2016-07-29T13:32:29Z<p>FelixMueller: /* Virtual system */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable <math> v </math> is integrated twice to produce <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
The virtual state and input are constrained by the sets \mathcal{Z} and \mathcal{V} respectively.<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Industrial_robot&diff=2130Industrial robot2016-07-29T13:30:27Z<p>FelixMueller: Added parameter sets</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 2<br />
|nc = 4<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Industrial robot problem focuses on the KUKA LWR IV drawing a pre-defined curve on a white board, which results in a path-following problem. This means that a system should be steered along a given path while there are no constraints regarding the position at a specific time <math> t > t_0 </math>. This concept arises often where autonomous vehicles or steering of robots are involved. The source for this problem can be found online on [https://arxiv.org/abs/1506.09084 arXiv "Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot"].<br />
<br />
<br />
== Path-following problems ==<br />
<br />
As the robot should follow a geometric curve but is not constrained as to when it has to be where, a path-following problem arises. The path is given by <br />
<math> \mathcal{P} = \{ y \in \mathcal{R}^n: \theta \in [ \theta_0, \theta_1 ] \mapsto y = p(\theta)\} </math><br />
where “<math> \theta </math> is called path parameter and <math> p(\theta) </math> is a parametrization of <math> \mathcal{P} </math>”. The path parameter does depend on the time however the exact timing is not known beforehand but rather integrated into the optimization problem by making use of a virtual system. <br />
<br />
<br />
== Virtual system ==<br />
<br />
A virtual system is introduced to take care of the timing of the robot. This is an additional degree of freedom which is used here as an input into an integrator. This additional optimization variable <math> v </math> is integrated twice to produce <math> z_1 </math> which is used as the time parameter for the path (i.e. <math> p(z_1) </math>).<br />
<br />
<math> \left(<br />
\begin{array}{c}<br />
\dot{z_1}\\<br />
\dot{z_2}\\<br />
\end{array}<br />
\right)<br />
=<br />
\dot{z}<br />
=<br />
l(z, v)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
z_2\\<br />
v\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
<br />
== Model ==<br />
The system for the motions of the robot are given by the function <math> f </math>:<br />
<br />
<math> f(x) = E^{-1} <br />
\left(<br />
\begin{array}{c}<br />
x_2\\<br />
u - C(x_1,x_2) x_2 - g(x_1) - \tau_F(x_2)\\<br />
\end{array}<br />
\right)<br />
</math><br />
<br />
with<br />
<br />
<math> E = diag(I, B(x_1)) </math><br />
<br />
where <math> B, C, E, g, \tau_F </math> are given in the Parameters section.<br />
<br />
The objective function <math> J </math> will be the integral over the so-called cost function <math> F: \mathbb{R}^{n_u} \times \mathbb{R} \times \mathcal{U} \times \mathcal{V} \mapsto \mathbb{R}_0^{+} </math>, i.e.<br />
<math> J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) = \displaystyle \int_{t_k}^{t_k+T} F(\bar{e}, \bar{z}, \bar{u}, \bar{v}) d \tau </math><br />
where <math> T \in (\delta, \infty) </math> is the prediction horizon and the function <math> F </math> is given by <math> F(e,z,u,v) = || (e,z_1-\theta_1, z_2 - \dot{\theta}_{ref})^T||_Q^2 + ||(u,v)^T||_R^2 </math> with positive semi definite matrix <math> Q </math> and positive definite matrix <math> R </math>, which can be found in the original paper cited above.<br />
<br />
<br />
The path-following is introduced into the problem by optimizing the error <math> e </math> measuring the deviation of the robot from the path: <br />
<br />
<math> e(x, z) = h(x) - p(z_1) </math><br />
<br />
where <math> p </math> is the parametrization of <math> \mathcal{P} </math> as above and <math> h </math> is the position of the robot given by:<br />
<br />
<math><br />
h(x_1)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
-\cos(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
-\sin(q_1)(d_3 \sin(q_2) + \bar{d}_5 + \sin(q_2-q_3))\\<br />
d_1 + d_3 \cos(q_2) + \bar{d}_5 \cos(q_2-q_3)\\ <br />
\end{array}<br />
\right)<br />
</math><br />
<br />
where <math> q = x_1 </math> and <math> \bar{d}_5 = d_5 + d_{tool} </math> with <math> d_i, i \in \{ 1,3,5 \} </math> being the lengths of the links of the robot and <math> d_{tool} </math> is the length of the pen and tool. These values can be found in the Parameters section as well.<br />
<br />
== OCP ==<br />
The optimization problem which is solved repetitively is given by<br />
<p><br />
<math><br />
\begin{array}{llclr}<br />
\displaystyle \min_{\bar{u}, \bar{v}} & J &=& J(x(t_k),\bar{e}, \bar{u}, \bar{z}, \bar{v}) \\[1.5ex]<br />
\mbox{s.t.}<br />
& \dot{\bar{x}} & = & f(\bar{x},\bar{u}), \\<br />
& \dot{\bar{z}} & = & l(\bar{z},\bar{v}), \\<br />
& \bar{e} & = & h(\bar{x}) - p(\bar{z}_1), \\[1ex]<br />
& \bar{x}(t_k) & = & x(t_k), \\<br />
& \bar{z}(t_k) & = & z(t_k), \\[1ex]<br />
& \bar{x} & \in & \mathcal{X}, \\<br />
& \bar{u} & \in & \mathcal{U}, \\<br />
& \bar{z} & \in & \mathcal{Z}, \\<br />
& \bar{v} & \in & \mathcal{V}. \\<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
== Parameters ==<br />
<br />
The matlab script which implements the implicit formulation for <math> f </math>, i.e. <math> F(x) = E \dot{x} - f(x) </math> can be found on [[Parameters for Industrial Robot]].<br />
<br />
{| class="wikitable"<br />
|+Robot-specific constants<br />
|-<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|<math> d_1 </math><br />
|<math> 0.31</math><br />
| <math> m </math><br />
|-<br />
|<math> d_3 </math><br />
|<math> 0.40</math><br />
| <math> m </math><br />
|-<br />
|<math> d_5 </math><br />
|<math> 0.39</math><br />
| <math> m </math><br />
|-<br />
|<math> d_{tool} </math><br />
|<math> 0.15</math><br />
| <math> m </math><br />
|}<br />
<br />
The Implementation Parameters for the MPFC scheme are as follows.<br />
<br />
* <math> \mathcal{X} = [-\bar{x}, \bar{x}], \text{ where } \bar{x} = (\bar{q},\bar{q},\bar{q},\bar{\dot{q}},\bar{\dot{q}},\bar{\dot{q}})^T \text{ with } \bar{q} = \infty \text{ rad/s}, \quad \bar{\dot{q}} = 0.5 \text{ rad/s} </math><br />
* <math> \mathcal{U} = [-\bar{u}, \bar{u}], \text{ where } \bar{u} = - (\bar{\tau},\bar{\tau},\bar{\tau})^T \text{ with } \bar{\tau} = 60 \text{ Nm} </math><br />
* <math> \mathcal{Z} = [\underline{z}, \bar{z}], \text{ where } \underline{z} = (0,0)^T, \bar{z} = (1750, \infty)^T </math><br />
* <math> \mathcal{V} = [-10^4, 8 \cdot 10^3] </math><br />
<br />
== Reference Solutions ==<br />
coming soon<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Batch_reactor&diff=2123Batch reactor2016-07-27T08:42:58Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 2<br />
|nre = 2<br />
}}<br />
<br />
This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.<br />
<br />
The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.<br />
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.<br />
<br />
== Mathematical formulation ==<br />
<br />
The optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \max_{x, u} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\<br />
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\<br />
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\<br />
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]<br />
& x(0) &=& (1, 0)^T, \\<br />
& T(t) &\in& [298, 398].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.<br />
<br />
== Parameters ==<br />
<br />
The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.<br />
<br />
== Reference Solutions ==<br />
<br />
This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].<br />
<br />
The optimal objective value of the problem is <math> x_2(t_f) = -0.611715 </math>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.<br />
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]<br />
<br />
== References ==<br />
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].<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:Chemical engineering]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Parameters_for_Industrial_Robot&diff=2122Parameters for Industrial Robot2016-07-27T08:37:10Z<p>FelixMueller: </p>
<hr />
<div>This page contains a matlab script for the implicit formulation of the parameters for the [[Industrial robot]] problem.<br />
<br />
<source lang="Matlab"><br />
% Robot dynamics with 3 DOF, this function should be equal to zero at all<br />
% times. e.g. f(t,x=[q,dq],dx=[dq,ddq])=0! This model is only containing<br />
% the inerta Matrix and the torques of the joints. Friction, gravity and coriolis<br />
% effects are incorporated.<br />
<br />
function diff=f_full_effects(~,y,dy,tau)<br />
<br />
dq=zeros(3,1);<br />
dotq=zeros(3,1);<br />
ddq=zeros(3,1);<br />
<br />
q=y(1:3);<br />
dq(1:3,1)=y(4:6);<br />
dotq(1:3,1)=dy(1:3);<br />
ddq(1:3,1)=dy(4:6);<br />
<br />
% intermediate states<br />
<br />
cosq=zeros(3,1);<br />
<br />
cosq(1,1) = cos(q(1));<br />
cosq(2,1) = cos(q(2));<br />
cosq(3,1) = cos(q(3));<br />
<br />
sinq=zeros(3,1);<br />
<br />
sinq(1,1) = sin(q(1));<br />
sinq(2,1) = sin(q(2));<br />
sinq(3,1) = sin(q(3));<br />
<br />
cosq2q3 = cos(q(2) - q(3));<br />
<br />
sinq1q2= sin(q(1) - q(2));<br />
sinq2q3= sin(q(2) - q(3));<br />
<br />
sin2q2= sin(2*q(2));<br />
<br />
sin2q22q3= sin(2*q(2) - 2*q(3));<br />
sin2q2q3= sin(2*q(2) - q(3));<br />
<br />
cos2q2= cos(2*q(2));<br />
<br />
% Modell parameters<br />
<br />
b_1 = 2.6491674286195;<br />
b_2 = 0.0168799610481578;<br />
b_3 = 1.43009916890323;<br />
b_4 = 0.53524835732352;<br />
b_5 = 0.791116620901176;<br />
b_6 = 0.626240578068518;<br />
b_7 = 1.57199673626875;<br />
b_8 = 0.0175195488386143;<br />
b_9 = 0.0204978721411787;<br />
b_10= 1.07049671464704;<br />
b_11= 1.81771647041588;<br />
b_12= 0.422954084243519;<br />
<br />
b=zeros(3,3); % Inerta Matrix of Robot Modell (3DOF)<br />
<br />
b(1,1) = b_1*(sinq(2)^2) - b_2*sin2q2 + b_3*(cosq(2)^2) + b_4*cosq(3) + b_5*(cosq2q3^2) - b_6*cos(2*q(2) - 2*q(3)) - b_4*cos(2*q(2) - q(3)) - b_7;<br />
b(2,1) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(3,1) = -b_9*cosq2q3;<br />
b(1,2) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(2,2) = b_10*cosq(3) + b_11;<br />
b(3,2) = -b_4*cosq(3) - b_12;<br />
b(1,3) = -b_9*cosq2q3;<br />
b(2,3) = -b_4*cosq(3) - b_12;<br />
b(3,3) = b_12;<br />
<br />
g=zeros(3,1); % gravity torque vector<br />
<br />
g_1 = 32.4854155934467;<br />
g_2 = 13.1269659633593;<br />
g_3 = 0.41398104470607;<br />
<br />
%g(1) = 0;<br />
%g(2) =-g_1*sinq(2) - g_2*sinq2q3 + g_3*cosq(2);<br />
%g(3) = g_2*sinq2q3;<br />
<br />
f=zeros(3,1); % friction torque vector<br />
<br />
f(1) = 0.186461043875677*atan(100*dq(1));<br />
f(2) = 0.184489618216264*atan(100*dq(2));<br />
f(3) = 0.156066608371045*atan(100*dq(3));<br />
<br />
c=zeros(3,3); % coriolis matrix<br />
<br />
c_1 = 0.609534129858137;<br />
c_2 = 0.23068226761793;<br />
c_3 = b_4;<br />
c_4 = b_2;<br />
c_5 = 0.26762417866176;<br />
c_6 = 152.383532464534;<br />
c_7 = 57.6705669044826;<br />
c_8 = 133.81208933088;<br />
c_9 = 4.21999026203945;<br />
c_10= 66.90604466544;<br />
c_11= b_8;<br />
c_12= b_9;<br />
c_13= 1.3381208933088;<br />
<br />
c(1,1) = c_1*dq(2)*sin2q2 + c_2*dq(2)*sin2q22q3 + c_3*dq(2)*sin2q2q3 - c_4*dq(2)*cos2q2 - c_5*dq(3)*sinq(3) - c_2*dq(3)*sin2q22q3 - c_5*dq(3)*sin2q2q3;<br />
c(2,1) = -(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250;<br />
c(3,1) = -(-1)*(-1)*dq(1)*(-c_10*sinq(3) - c_7*sin2q22q3 - c_10*sin2q2q3)/250;<br />
c(1,2) = -(-1)*(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250 + dq(2)*(-c_11*sinq(2) - c_12*sinq2q3) + c_12*dq(3)*sinq2q3;<br />
c(2,2) = -c_3*dq(3)*sinq(3);<br />
c(3,2) = c_3*dq(2)*sinq(3);<br />
c(1,3) = -c_5*dq(1)*sinq(3) - c_2*dq(1)*sin2q22q3 - c_5*dq(1)*sin2q2q3 + c_12*dq(2)*sinq2q3 - c_12*dq(3)*sinq2q3;<br />
c(2,3) = c_13*(-(-1)*(-2)*dq(2)/5 - (-1)*2*dq(3)/5)*sinq(3);<br />
c(3,3) = 0; <br />
<br />
<br />
diff=[dq-dotq;...<br />
b*ddq+g+c*dq+f-tau'];<br />
<br />
</source></div>FelixMuellerhttps://mintoc.de/index.php?title=Van_der_Pol_Oscillator&diff=2121Van der Pol Oscillator2016-07-27T08:34:16Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 1<br />
|nre = 2<br />
}}<br />
<br />
The Van der Pol Oscillator is an oscillating system with non-linear damping and self regulation. The System was first introduced by the Dutch physician Balthasar Van der Pol in 1927. The aim is to control the oscillation such that the system stays in a mean position. <br />
<br />
<br />
== Model formulation ==<br />
<br />
The Van der Pol Oscillator evolves over time according to the second order differential equation:<br />
<br />
<math>{d^2x \over dt^2}-u(1-x^2){dx \over dt}+x= 0</math><br />
<br />
where <math>x</math> is the position coordinate, which is a function of the time <math>t</math>, and <math>u</math> is a scalar parameter indicating the non-linearity and the strength of the damping.<br />
<br />
For <math>u>>1</math> the oscillator is being damped, whereas for <math>u<<1</math> energy is added to the system.<br />
<br />
Based on the transformation <math>y = \dot x</math> the problem can be reformulated:<br />
<br />
<math> \dot x = y</math><br />
<br />
<math>\dot y = u(1-x^2) y-x</math><br />
<br />
<br />
The Optimal Control Problem arises by adding the objective function:<br />
<br />
<math>\min\limits_{u}\int\limits_{t_0}^{t_f}(x(t)^2+y(t)^2+u(t)^2) dt</math><br />
<br />
<br />
<br />
== Optimal Control Problem ==<br />
The Optimal Control Problem with the aim to minimize the deflection can be formulated as follows:<br />
<br />
<math> <br />
\begin{array}{lll}<br />
\min\limits_{u} & \int\limits_{t_0}^{t_f} & (x(t)^2+y(t)^2+u(t)^2) dt\\<br />
s.t. & \dot x & = y,\\<br />
& \dot y & = u(1-x^2) y-x,\\<br />
& x(0) & =1,\\<br />
& y(0) & =0,\\<br />
& u(t) & \le 0.75.\\<br />
\end{array}</math><br />
<br />
== Parameters ==<br />
These fixed values are used within the model:<br />
<br />
<math>[t_0,t_f]=[0,20]</math><br />
<br />
<br />
<br />
== Reference Solution ==<br />
<br />
The following reference solution was generated using JModelica with the automatic differentiation tool CasADI and the solver IPOPT. The Optimica code used to solve the problem can be found under [[Van der Pol Oscillator (JModelica)]] The optimal value of the objective function is 3.1762. <br />
<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:VDP_Plot_control.png| Control u over time (damping).<br />
Image:VDP_Plot_states.png| Position coordinate x and it's derivative y.<br />
Image:VDP_Plot_derivatives.png| Derivative <math>\dot x</math> and second derivative <math>\dot y</math>.<br />
</gallery><br />
<br />
== Source Code ==<br />
Model descriptions are available in:<br />
<br />
* [[:Category: JModelica | JModelica code]] at [[Van der Pol Oscillator (JModelica)]]<br />
* [[:Category: Julia/JuMP | JuMP code]] (of a slightly modified Van der Pol oscillator problem) at [[Van der Pol Oscillator (Jump)]]<br />
<br />
== References ==<br />
The Problem can be found under the following [https://en.wikipedia.org/wiki/Van_der_Pol_oscillator link] or in the [http://www.jmodelica.org/api-docs/usersguide/1.4.0/ch08s02.html JModelica Users Guide].<br />
<br />
[[Category:MIOCP]]<br />
[[Category:Bang bang]]<br />
[[Category:Path-constrained arcs]]<br />
[[Category:ODE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Supermarket_refrigeration_system&diff=2120Supermarket refrigeration system2016-07-27T08:33:51Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 9<br />
|nw = 3<br />
|nc = 7<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} &=& \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} &=& - \dfrac{<br />
UA_{goods-air} \left( x_1 - x_3 \right)<br />
}{<br />
M_{goods} \cdot C_{p,goods} <br />
} \\ <br />
&\dot{x_2} &=& \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} &=& \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} &=& \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} &=& - \dfrac{<br />
UA_{goods-air} \left( x_5 - x_7 \right)<br />
}{<br />
M_{goods} \cdot C_{p,goods} <br />
} \\ <br />
&\dot{x_6} &=& \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} &=& \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} &=& \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 &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_3 &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_7 &\geq& 2.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_7 &\leq& 5.0 \quad \forall t \in [t_0, t_f],\\<br />
& x_0 &\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:JModelica | JModelica code]] at [[Supermarket refrigeration system (JModelica)]]<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>FelixMuellerhttps://mintoc.de/index.php?title=Particle_steering_problem&diff=2119Particle steering problem2016-07-27T08:33:24Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nw = 1<br />
|nc = 2<br />
|nre = 7<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Particle steering problem minimizes "the time taken for a particle, acted upon by a thrust of constant magnitude, to achieve a given altitude and terminal velocity." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u, t_f} & t_f \\[1.5ex]<br />
\mbox{s.t.} <br />
& \ddot{x}_1 & = & a \cos (u), \\<br />
& \ddot{x}_2 & = & a \sin (u), \\<br />
& x(0) &=& (0, 0)^T, \\<br />
& \dot{x}(0) &=& (0, 0)^T, \\<br />
& x_2 (t_f) &=& 5, \\<br />
& \dot{x}(t_f) &=& (45, 0)^T, \\<br />
& u(t) &\in& [-\frac{\pi}{2},\frac{\pi}{2}].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<br />
where <math> (x_1, x_2) </math> is the position of the particle, <math> u </math> is the control angle and <math> a </math> is the constant magnitude of thrust.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Particle steering problem (TACO)]]<br />
<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:Minimum time]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Oil_Shale_Pyrolysis&diff=2118Oil Shale Pyrolysis2016-07-27T08:33:15Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nw = 0<br />
|nc = 4<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}{lll}<br />
\displaystyle \min_{u} & \displaystyle &-x_1(t_N)^2 \\[1.5ex]<br />
\mbox{s.t.} & \displaystyle \dot{x}_0 &= -k_0x_0-(k_2+k_3+k_4)x_0x_1\\<br />
& \displaystyle \dot{x}_1 &= k_0x_0-k_1x_1 + k_2x_0x_1\\<br />
& \displaystyle k_i &= a_i e^{-u\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.15748.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>FelixMuellerhttps://mintoc.de/index.php?title=Catalytic_cracking_problem&diff=2117Catalytic cracking problem2016-07-27T08:32:44Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 2<br />
|np = 3<br />
|nc = 3<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalytic cracking problem tries to determine "reaction coefficients for the catalytic cracking of gas oil into gas and other byproducts." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{21} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_3) y_1^2, \\<br />
& \dot{y}_2 & = & \theta_1 y_1^2 - \theta_2 y_2, \\<br />
& \theta_i & \geq & 0 \quad i = 1,...,3.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_{21} </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalytic cracking problem (TACO)]]<br />
<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:DAE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Isomerization_of_Alpha-Pinene_problem&diff=2116Isomerization of Alpha-Pinene problem2016-07-27T08:32:36Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 5<br />
|np = 5<br />
|nc = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Isomerization of Alpha-Pinene problem tries to determine "reaction coefficients in the thermal isometrization of <math> \alpha </math>-Pinene." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{8} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_2) y_1, \\<br />
& \dot{y}_2 & = & \theta_1 y_1, \\<br />
& \dot{y}_3 & = & \theta_2 y_1 - (\theta_3 + \theta_4) y_3 + \theta_5 y_5, \\<br />
& \dot{y}_4 & = & \theta_3 y_3, \\<br />
& \dot{y}_5 & = & \theta_4 y_3 - \theta_5 y_5, \\<br />
& \theta_i & \geq & 0 \quad i = 1,...,5.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_8 </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Isomerization of Alpha-Pinene problem (TACO)]]<br />
<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:DAE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Methanol_to_Hydrocarbons_problem&diff=2115Methanol to Hydrocarbons problem2016-07-27T08:32:00Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 3<br />
|np = 5<br />
|nc = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Methanol to Hydrocarbons problem tries to determine "reaction coefficients for the conversion of methanol into various hydrocarbons." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{16} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -( 2 \theta_2 - \frac{\theta_1 y_2}{(\theta_2 + \theta_5) y_1 + y_2} + \theta_3 + \theta_4) y_1, \\[0.3cm]<br />
& \dot{y}_2 & = & \frac{\theta_1 y_1 (\theta_2 y_1 - y_2)}{(\theta_2 + \theta_5) y_1 + y_2} + \theta_3 y_1, \\[0.3cm]<br />
& \dot{y}_3 & = & \frac{\theta_1 y_1 (y_2 + \theta_5 y_1}{(\theta_2 + \theta_5) y_1 + y_2} + \theta_4 y_1, \\[0.3cm]<br />
& \theta_i & \geq & 0 \quad i = 1,...,5.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_{16} </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Methanol to Hydrocarbons problem (TACO)]]<br />
<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:Chemical engineering]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Marine_population_dynamics_problem&diff=2114Marine population dynamics problem2016-07-27T08:31:49Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = <math> n_s </math><br />
|np = <math> 2 n_s </math><br />
|nc = <math> 4 n_s </math><br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Marine population dynamics problem estimates growth and mortality rates of a marine species at each stage (for example ages or development stage) given the population as a function of time.( Problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{g, m} & \sum\limits_{j=1}^{n_s} &&||y(\tau_j; g, m) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_j & = & g_{j-1} y_{j-1} - (m_j + g_j) y_j \qquad \forall j \in 1, ..., n_s,\\<br />
& g_j, m_j &\in& [0,1].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
where <math> g_j </math> and <math> m_j </math> are the growth and mortality rates at stage <math> j </math> respectively and the initial conditions are unknown.<br />
The error between computed and observed data is minimized.<br />
<br />
== Parameters ==<br />
There are <math> n_s </math> stages and <math> n_m </math> timepoints at which the error is minimized.<br />
<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Marine population dynamics problem (TACO)]]<br />
<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:Population dynamics]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Lotka_Experimental_Design&diff=2113Lotka Experimental Design2016-07-27T08:31:11Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 11<br />
|nw = 2<br />
|nc = 4<br />
|nre = 11<br />
}}<br />
<br />
The '''Lotka Experimental Design problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to control the biomasses of both predator as prey fish as described in [[Lotka Volterra fishing problem]]. The goal here, however, is to minimize the uncertainty of a follow-up parameter estimation problem. In addition to the fishing decision the question when measurements of the two state variables are performed becomes a degree of freedom.<br />
<br />
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The ODE from [[Lotka Volterra fishing problem]] is extended such that it also includes state sensitivities, the Fisher information matrix entries and integrated sampling states.<br />
<br />
The optimal integer control functions shows [[:Category:Chattering|bang bang]] behavior.<br />
<br />
== Mathematical formulation ==<br />
<br />
We are interested in estimating the parameters <math>p_2</math> and <math>p_4</math> of the Lotka-Volterra type predator-prey fish initial value problem<br />
<p><br />
<math><br />
\begin{array}{rcl}<br />
\dot{x_1}(t) &=& p_1 \; x_1(t) - p_2 x_1(t) x_2(t) - p_5 u(t) x_1(t), \; t \in [0,t_f], \quad x_1(0) = 0.5, \\<br />
\dot{x_2}(t) &=& - p_3 \; x_2(t) + p_4 x_1(t) x_2(t) - p_6 u(t) x_2(t), \; t \in [0,t_f], \quad x_2(0) = 0.7,<br />
\end{array} <br />
</math><br />
</p><br />
where <math>u(\cdot)</math> is a fishing control that may or may not be fixed. The other parameters, the initial values and <math>t_f = 12</math> are fixed. We are interested in how to fish and when to measure, with an upper bound <math>M</math> on the measuring time. We can measure the states directly, <math>h^1(x(t)) = x_1(t)</math> and <math>h^2(x(t)) = x_2(t)</math>. We use two different sampling functions, <math>w^1(\cdot)</math> and <math>w^2(\cdot)</math> in the same experimental setting. This can be seen either as a two-dimensional measurement function <math>h(x(t))</math>, or as a special case of a multiple experiment, in which <math>u(\cdot), x(\cdot)</math>, and <math>G(\cdot)</math> are identical. The experimental design problem then reads<br />
<br />
<p><br />
<math><br />
\begin{array}{lll}<br />
\displaystyle \min_{x,G,F,z^1,z^2,u,w^1,w^2} && \text{trace} \; \left( F^{-1}(t_f) \right) \\<br />
\text{subject to} \\<br />
\quad \dot{x_1}(t) & = & p_1 \; x_1(t) - p_2 x_1(t) x_2(t) - p_5 u(t) x_1(t),\\<br />
\quad \dot{x_2}(t) & = & - p_3 \; x_2(t) + p_4 x_1(t) x_2(t) - p_6 u(t) x_2(t),\\<br />
\quad \dot{G_{11}}(t) & = & f_{x11}(\cdot) \; G_{11}(t) + f_{x12}(\cdot) \; G_{21}(t) + f_{p12}(\cdot), \\<br />
\quad \dot{G_{12}}(t) & = & f_{x11}(\cdot) \; G_{12}(t) + f_{x12}(\cdot) \; G_{22}(t), \\<br />
\quad \dot{G_{21}}(t) & = & f_{x21}(\cdot) \; G_{11}(t) + f_{x22}(\cdot) \; G_{21}(t), \\<br />
\quad \dot{G_{22}}(t) & = & f_{x21}(\cdot) \; G_{12}(t) + f_{x22}(\cdot) \; G_{22}(t) + f_{p24}(\cdot), \\<br />
\quad \dot{F_{11}}(t) & = & w^1(t) G_{11}(t)^2 + w^2(t) G_{21}(t)^2, \\<br />
\quad \dot{F_{12}}(t) & = & w^1(t) G_{11}(t) G_{12}(t) + w^2(t) G_{21}(t) G_{22}(t), \\<br />
\quad \dot{F_{22}}(t) & = & w^1(t) G_{12}(t)^2 + w^2(t) G_{22}(t)^2, \\<br />
\quad \dot{z^1}(t) & = & w^1(t), \\<br />
\quad \dot{z^2}(t) & = & w^2(t), \\[1.5ex]<br />
\quad x(0) &=& (0.5, 0.7), \\<br />
\quad G(0) &=& F(0) = 0, \\<br />
\quad z^1(0) &=& z^2(0) = 0, \\ [1.5ex]<br />
\quad u(t) & \in & \mathcal{U}, \; w^1(t) \in \mathcal{W}, \; w^2(t) \in \mathcal{W}, \\<br />
\quad 0 & \le & M - z(t_f)<br />
\end{array}<br />
</math><br />
</p><br />
<br />
with <br />
<math>f_{x11}(\cdot) = \partial f_1(\cdot) / \partial x_1 = p_1 - p_2 x_2(t) - p_5 u(t)</math>, <br />
<math>f_{x12}(\cdot) = - p_2 x_1(t)</math>, <br />
<math>f_{x21}(\cdot) = p_4 x_2(t)</math>, <br />
<math>f_{x22}(\cdot) = -p_3 + p_4 x_1(t) - p_6 u(t)</math>, and <br />
<math>f_{p12}(\cdot) = \partial f_1(\cdot) / \partial p_2 = -x_1(t) x_2(t)</math>, <br />
<math>f_{p24}(\cdot) = \partial f_2(\cdot) / \partial p_4 = x_1(t) x_2(t)</math>. <br />
<br />
Note that the state <math>F_{21}(\cdot) = F_{12}(\cdot)</math> has been left out for reasons of symmetry.<br />
<br />
== Parameters ==<br />
<br />
We use <math>t_f=12</math>, <math>p_1 = p_2 = p_3 = p_4 = 1</math>, and <math>p_5 = 0.4</math>, <math>p_6 = 0.2</math>. The upper bound on the measurement time intervals is chosen as <math>M=4</math>.<br />
<br />
== Reference Solutions ==<br />
<br />
Here is one local solution to the above control problem.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:vplotkaUStates.png| States and fishing control.<br />
Image:vplotkaUSensitivities.png| Sensitivities G().<br />
Image:vplotkaUMeas1.png| Sampling function for first state.<br />
Image:vplotkaUMeas2.png| Sampling function for second state.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL | AMPL]] at [[Lotka Experimental Design (AMPL)]]<br />
* [[:Category: VPLAN | VPLAN code]] at [[Lotka Experimental Design (VPLAN)]]<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 Experimental Design (AMPL)]],<br />
* no fishing, i.e., <math>u \equiv 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. The experimental design problem has been described in the habilitation thesis of Sager, <bib id="Sager2011d" />.<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:Optimum Experimental Design]]<br />
[[Category:ODE model]]<br />
[[Category:Bang bang]]<br />
[[Category:Population dynamics]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Isomerization_of_Alpha-Pinene_problem&diff=2112Isomerization of Alpha-Pinene problem2016-07-27T08:30:15Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 5<br />
|np = 5<br />
|nc = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Isomerization of Alpha-Pinene problem tries to determine "reaction coefficients in the thermal isometrization of <math> \alpha </math>-Pinene." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{8} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_2) y_1, \\<br />
& \dot{y}_2 & = & \theta_1 y_1, \\<br />
& \dot{y}_3 & = & \theta_2 y_1 - (\theta_3 + \theta_4) y_3 + \theta_5 y_5, \\<br />
& \dot{y}_4 & = & \theta_3 y_3, \\<br />
& \dot{y}_5 & = & \theta_4 y_3 - \theta_5 y_5, \\<br />
& \theta_i & \geq & 0.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_8 </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Isomerization of Alpha-Pinene problem (TACO)]]<br />
<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:DAE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Hanging_chain_problem&diff=2111Hanging chain problem2016-07-27T08:29:52Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nu = 1<br />
|nc = 4<br />
|nre = 5<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Hanging chain problem is concerned with finding a chain (of uniform density) of length <math> L </math> suspendend between two points <math> a, b </math> with minimal potential energy. (Problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_1 & = & u, \\<br />
& \dot{x}_2 & = & x_1 (1+u^2)^{1/2}, \\<br />
& \dot{x}_3 & = & (1+u^2)^{1/2}, \\<br />
& x(t_0) &=& (a,0,0)^T, \\<br />
& x_1(t_f) &=& b, \\<br />
& x_3(t_f) &=& Lp, \\<br />
& x(t) &\in& [0,10], \\<br />
& u(t) &\in& [-10,20].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
In this model the parameters used are<br />
<br />
<math><br />
\begin{array}{rcl}<br />
[t_0, t_f] &=& [0, 1],\\<br />
(a,b) &=& (1,3),\\<br />
Lp &=& 4.<br />
\end{array}<br />
</math><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Hanging chain problem (TACO)]]<br />
<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:Minimum energy]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Gravity_Turn_Maneuver&diff=2110Gravity Turn Maneuver2016-07-27T08:29:42Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 5<br />
|nu = 1<br />
|nc = 11<br />
|nri = 2<br />
|nre = 7<br />
}}<br />
<br />
The '''gravity turn''' or '''zero lift turn''' is a common maneuver used to launch spacecraft into orbit from bodies that have non-negligible atmospheres. The goal of the maneuver is to minimize atmospheric drag by always orienting the vehicle along the velocity vector. In this maneuver, the vehicle's pitch is determined solely by the change of the velocity vector through gravitational acceleration and thrust. The goal is to find a launch configuration and thrust control strategy that achieves a specific orbit with minimal fuel consumption.<br />
<br />
== Physical description and model derivation ==<br />
<br />
For the purposes of this model, we start with the following ODE system proposed by Culler et. al. in <bib id="Culler1957" />:<br />
<br />
<dd><br />
<math><br />
\begin{array}{rcl}<br />
\dot{v} &=& \frac{F}{m} - g \cdot \cos \beta \\<br />
v \dot{\beta} &=& g \cdot \sin{\beta}<br />
\end{array}<br />
</math><br />
</dd><br />
<br />
where <math>v</math> is the speed of the vehicle, <math>g</math> is the gravitational acceleration at the vehicle's current altitude, <math>F</math> is the accelerating force and <math>\beta</math> is the angle between the vertical and the vehicle's velocity vector. In the original version of the model, the authors neglect aspects of the problem:<br />
* Variation of <math>g</math> over altitude<br />
* Decrease of vehicle mass due to fuel consumption<br />
* Curvature of the surface<br />
* Atmospheric drag<br />
<br />
=== Changes in gravitational acceleration ===<br />
<br />
To account for changes in <math>g</math>, we make the following substitution:<br />
<dd><br />
<math><br />
g = g_0 \cdot \left(\frac{r_0}{r_0 + h}\right)^2<br />
</math><br />
</dd><br />
where <math>g_0</math> is the gravitational acceleration at altitude zero and <math>r_0</math> is the distance of altitude zero from the center of the reference body.<br />
<br />
=== Decrease in vehicle mass ===<br />
<br />
To account for changes in vehicle mass, we consider <math>m</math> a differential state with the following derivative:<br />
<dd><br />
<math><br />
\dot{m} = -\frac{F}{I_{sp} \cdot g_0}<br />
</math><br />
</dd><br />
where <math>I_{sp}</math> is the [https://en.wikipedia.org/wiki/Specific_impulse specific impulse] of the vehicle's engine. Specific impulse is a measure of engine efficiency. For rocket engines, it directly correlates with the engine's exhaust velocity and may vary with atmospheric pressure, velocity, engine temperature and combustion dynamics. For the purposes of this model, we will assume it to be constant.<br />
<br />
The vehicle's fuel reserve is modelled by two parameters: <math>m_0</math> denotes the launch mass (with fuel) and <math>m_1</math> denotes the dry mass (without fuel).<br />
<br />
=== Curvature of the reference body's surface ===<br />
<br />
To accomodate the reference body's curvature, we introduce an additional differential state <math>\theta</math> which represents the change in the vehicle's polar angle with respect to the launch site. The derivative is given by<br />
<br />
<dd><br />
<math><br />
\dot{\theta} = \frac{v \cdot \sin \beta}{r_0 + h}.<br />
</math><br />
</dd><br />
<br />
Note that the vertical changes as the vehicle moves around the reference body meaning that the derivative of <math>\beta</math> must be changed as well:<br />
<br />
<dd><br />
<math><br />
\dot{\beta} = \frac{g \cdot \sin{\beta}}{v} - \frac{v \cdot \sin \beta}{r_0 + h}.<br />
</math><br />
</dd><br />
<br />
=== Atmospheric drag ===<br />
<br />
To model [https://en.wikipedia.org/wiki/Drag_%28physics%29 atmospheric drag], we assume that the vehicles draf coefficient <math>c_d</math> is constant. The drag force is given by<br />
<br />
<dd><br />
<math><br />
F_{drag} = \frac{1}{2} \rho A c_d v^2<br />
</math><br />
</dd><br />
<br />
where <math>\rho</math> is the density of the atmosphere and <math>A</math> is the vehicle's reference area. We assume that atmospheric density decays exponentially with altitude:<br />
<br />
<dd><br />
<math><br />
\rho = \rho_0 \cdot e^{-\frac{h}{H}}<br />
</math><br />
</dd><br />
<br />
where <math>\rho_0</math> is the atmospheric density at altitude zero and <math>H</math> is the scale height of the atmosphere. The [drag force] is introduced into the acceleration term:<br />
<br />
<dd><br />
<math><br />
\dot{v} = \frac{F - F_{drag}}{m} - g \cdot \cos \beta.<br />
</math><br />
</dd><br />
<br />
Note that if the vehicle is axially symmetric and oriented in such a way that its symmetry axis is parallel to the velocity vector, it does not experience any lift forces. This model is simplified. It does not account for changes in temperature and atmospheric composition with altitude. Also, <math>c_d</math> varies with fluid viscosity and vehicle velocity. Specifically, drastic changes in <math>c_d</math> occur as the vehicle breaks the sound barrier. This is not accounted for in this model.<br />
<br />
== Mathematical formulation ==<br />
<br />
The resulting optimal control problem is given by<br />
<br />
<dd><br />
<math><br />
\begin{array}{llcll}<br />
\displaystyle \min_{T, m, v, \beta, h, \theta, u} & m_0 - m(T) \\[1.5ex]<br />
\mbox{s.t.} & \dot{m} & = & -\frac{F_{max}}{I_{sp} \cdot g_0} \cdot u ,\\<br />
& \dot{v} & = & (F - \frac{1}{2} \rho_0 \cdot e^{-\frac{h}{H}} A c_d v^2) \cdot \frac{1}{m} - g_0 \cdot \left(\frac{r_0}{r_0 + h}\right)^2 \cos \beta ,\\<br />
& \dot{\beta} & = & g_0 \cdot \left(\frac{r_0}{r_0 + h}\right)^2 \cdot \frac{\sin{\beta}}{v} - \frac{v \cdot \sin \beta}{r_0 + h} ,\\<br />
& \dot{h} & = & v \cdot \cos \beta ,\\<br />
& \dot{\theta} & = & \frac{v \cdot \sin \beta}{r_0 + h},\\[1.5ex]<br />
& m(0) &=& m_0 , \\<br />
& v(0) &=& \varepsilon , \\<br />
& \beta(0) &\in& \left[ 0,\frac{\pi}{2} \right], \\<br />
& h(0) &=& 0, \\<br />
& \theta(0) &=& 0, \\[1.5ex]<br />
& \beta(T) &=& \hat{\beta} , \\<br />
& h(T) &=& \hat{h} , \\<br />
& v(T) &=& \hat{v} , \\[1.5ex]<br />
& m(t) &\in& [m_1, m_0] \qquad & \forall t \in [0,T] , \\<br />
& v(t) &\geq& \varepsilon \qquad & \forall t \in [0,T] , \\<br />
& \beta(t) &\in& [0, \pi] \qquad & \forall t \in [0,T] , \\<br />
& h(t) &\geq& 0 \qquad & \forall t \in [0,T] , \\<br />
& \theta(t) &\geq& 0 \qquad & \forall t \in [0,T] , \\[1.5ex]<br />
& u(t) &\in& [0,1] \qquad & \forall t \in [0,T] , \\<br />
& T &\in& [T_{min}, T_{max}].<br />
\end{array} <br />
</math><br />
</dd><br />
<br />
where <math>F_{max}</math> is the maximum thrust of the vehicle's engine and <math>\varepsilon</math> is a small number that is strictly greater than zero. <br />
The differential states of the problem are <math> m, v, \beta, h, \theta </math> while <math> u </math> is the control function<br />
<br />
== Parameters ==<br />
<br />
For testing purposes, the following parameters were chosen:<br />
<br />
<br />
<dd><br />
<math><br />
\begin{array}{rcl}<br />
m_0 &=& 11.3 \; t \\<br />
m_1 &=& 1.3 \; t \\<br />
I_{sp} &=& 300 \; s \\<br />
F_{max} &=& 0.6 \; MN \\<br />
c_d &=& 0.021 \\<br />
A &=& 1 \; m^2 \\[1.5ex]<br />
g_0 &=& 9.81 \cdot 10^{-3} \, \frac{km}{s^2} \\<br />
r_0 &=& 600.0 \; km \\<br />
H &=& 5.6 \; km \\<br />
\rho_0 &=& 1.2230948554874 \; \frac{kg}{m^3} \\[1.5ex]<br />
\hat{\beta} &=& \; \frac{\pi}{2} \\<br />
\hat{v} &=& 2.287 \; \frac{km}{s} \\<br />
\hat{h} &=& 75 \; km \\[1.5ex]<br />
T_{min} &=& 120 \; s \\<br />
T_{max} &=& 600 \; s \\<br />
\varepsilon &=& 10^{-6} \; \frac{km}{s}<br />
\end{array}<br />
</math><br />
</dd><br />
<br />
== Reference solution ==<br />
<br />
The reference solution was generated using a direct multiple shooting approach with 300 shooting intervals placed equidistantly along the variable-length time horizon. The resulting NLP was implemented using CasADi 2.4.2 using CVodes as an integrator and IPOPT as a solver. The exact code used to solve the problem alongside detailed solution data can be found under [[Gravity Turn Maneuver (Casadi)]]. The solution achieves the desired orbit in <math>T = 338.042 \; s</math> spending <math>m_0 - m(T) = 9.6508 \; t</math> of fuel. The launch angle is approximately <math>\beta(0) = 2.7786081986275378 \cdot 10^{-6}</math>. The downrange motion throughout the entire launch amounts to a change in polar angle of approximately <math>\theta(T) = 0.3383474346340064 \approx 19.39^{\circ}</math>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
Image:gravity-turn-altitude.png| Vehicle altitude over time.<br />
Image:gravity-turn-beta.png| Angle <math>\beta</math> over time.<br />
Image:gravity-turn-mass.png| Vehicle mass over time.<br />
Image:gravity-turn-theta.png| Angle <math>\theta</math> over time.<br />
Image:gravity-turn-velocity.png| Vehicle velocity over time.<br />
Image:gravity-turn-control.png| Discretized engine control over time.<br />
Image:gravity-turn-trajectory.png| Resulting vehicle trajectory relative to reference body surface.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category: Casadi | Casadi code]] at [[Gravity Turn Maneuver (Casadi)]]<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:Aeronautics]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=F-8_aircraft&diff=2109F-8 aircraft2016-07-27T08:29:12Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nw = 1<br />
|nc = 2<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 />
=== jModelica ===<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 code]] at [[F-8 aircraft (AMPL)]]<br />
* [[:Category:Muscod | Muscod code]] at [[F-8 aircraft (Muscod)]]<br />
* [[:Category:JModelica | JModelica code]] at [[F-8 aircraft (JModelica)]]<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>FelixMuellerhttps://mintoc.de/index.php?title=Electric_Car&diff=2108Electric Car2016-07-27T08:29:02Z<p>FelixMueller: </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 />
-->The '''Electric car problem''' tries to find an optimal driving policy for an electric car. The goal is to use minimal energy to finish a given distance. As the car can be driven in two discrete modes, which either cause acceleration (and thereby consumption of energy) or the recharging of the battery, the control variable <math> u(t) </math> is supposed to be integer. Additionally the model for the electric car itself contains nonlinearities.<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:MIOCP]]<br />
[[Category:ODE model]]<br />
[[Category:Chattering]]<br />
[[Category: Bang bang]]<br />
[[Category: Path-constrained arcs]]<br />
[[Category:Sensitivity-seeking arcs]]<br />
[[Category: Transport]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Catalyst_mixing_problem&diff=2106Catalyst mixing problem2016-07-27T08:28:28Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 2<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalyst mixing problem seeks an optimal policy for mixing two catalysts "along the length of a tubular plug <br />
ow reactor involving several reactions". (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} &-1 + x_1(t_f) + x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_1 & = & u ( 10 x_2 - x_1), \\<br />
& \dot{x}_2 & = & u ( x_1 - 10 x_2) - (1 - u \, x_2) , \\<br />
& x(t_0) &=& (1, 0)^T, \\<br />
& u(t) &\in& [0,1].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
In this model the parameters used are <math> t_0 = 0, \, \, t_f = 1 </math>.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalyst mixing problem (TACO)]]<br />
<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]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Batch_reactor&diff=2105Batch reactor2016-07-27T08:28:09Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 2<br />
|nre = 2<br />
}}<br />
<br />
This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.<br />
<br />
The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.<br />
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.<br />
<br />
== Mathematical formulation ==<br />
<br />
The optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \max_{x, u} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\<br />
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\<br />
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\<br />
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]<br />
& x(0) &=& (1, 0)^T, \\<br />
& T(t) &\in& [298, 398].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.<br />
<br />
== Parameters ==<br />
<br />
The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.<br />
<br />
== Reference Solutions ==<br />
<br />
This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].<br />
<br />
The optimal objective value of the problem is <math> x_2(t_f) = -0.611715 </math>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.<br />
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]<br />
<br />
== References ==<br />
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].<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:Chemical engineering]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Bioreactor&diff=2104Bioreactor2016-07-27T08:27:50Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nu = 1<br />
|nc = 2<br />
|nre = 3<br />
}}<br />
<br />
The bioreactor problem describes an substrate that is converted to a product by the biomass in the reactor. It has three states and a control that is describing the feed concentration of the substrate. The problem is taken from the examples folder of the ACADO toolkit described in:<br />
<bib id="Houska2011a" /><br />
<br />
Houska, Boris, Hans Joachim Ferreau, and Moritz Diehl. <br />
"ACADO toolkit—An open‐source framework for automatic control and dynamic optimization."<br />
Optimal Control Applications and Methods 32.3 (2011): 298-312.<br />
<br />
Originally the problem seems to be motivated by:<br />
<bib id="Versyck1999" /><br />
<br />
VERSYCK, KARINA J., and JAN F. VAN IMPE. <br />
"Feed rate optimization for fed-batch bioreactors: From optimal process performance to optimal parameter estimation."<br />
Chemical Engineering Communications 172.1 (1999): 107-124.<br />
<br />
== Model Formulation ==<br />
<br />
The dynamic model is an [[:Category:ODE model|ODE model]]:<br />
<br />
<p><br />
<math><br />
\begin{array}{rcl}<br />
\dot{X}&=&-DX+\mu X \\<br />
\dot{S}&=& D(S_{f}-S)-\mu /Y_{xs} X \\<br />
\dot{P}&=&-DP+ (\alpha \mu +\beta) X.<br />
\end{array}<br />
</math><br />
</p><br />
<br />
The right-hand side of these equations will be summed up in <math> f(x, S_f) </math>.<br />
<br />
The three states describe the concentration of the biomass (<math>X</math>), the substrate (<math>S</math>), and the product (<math>P</math>) in the reactor. In steady state the feed and outlet are equal and dilute all three concentrations with a ratio <math>D</math>. The biomass grows with a rate<br />
<math>\mu</math>, while it eats up the substrate with the rate <math>\mu/Y_{xs}</math> and produces product at a rate <math>(\alpha \mu +\beta)</math>. The rate <math>\mu</math> is given by:<br />
<br />
<math>\mu = \mu_{m}*(1-P/P_{m})*S/(K_m+S+S^2/K_i)</math><br />
<br />
The fixed parameters (constants) of the model are as follows.<br />
<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|Dilution<br />
|<math>D</math><br />
|0.15<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>K_i</math><br />
|22<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>K_m</math><br />
|1.2<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>P_m</math><br />
|50<br />
|[-]<br />
|-<br />
|Substrate to Biomass rate<br />
|<math>Y_{xs}</math><br />
|0.4<br />
|[-]<br />
|-<br />
|Linear slope<br />
|<math>\alpha</math><br />
|2.2<br />
|[-]<br />
|-<br />
|Linear intercept<br />
|<math>\beta</math><br />
|0.2<br />
|[-]<br />
|-<br />
|Maximal growth rate<br />
|<math>\mu_m</math><br />
|0.48<br />
|[-]<br />
|-<br />
|}<br />
<br />
== Mathematical formulation ==<br />
<br />
Writing shortly for the states in vector notation <math>x=(X,S,P)^T</math> the OCP reads:<br />
<br />
<p><br />
<math><br />
\begin{array}{clcl}<br />
\displaystyle \min_{x,S_f} & J(x,S_f)\\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x} & = & f(x,S_f)\\<br />
& x(0) & = & (6.5,12,22)^T \\<br />
& S_f & \in &[28.7,40].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
=== Objective ===<br />
<p><br />
<math><br />
J(x,S_f)=\int_0^{48}D(S_f-P)^2dt<br />
</math><br />
</p><br />
<br />
== Reference Solution ==<br />
<br />
Here we present the reference solution of the reimplemented example in the ACADO code generation with matlab. The source code is given in the next section.<br />
<br />
<gallery caption="Reference solution" widths="551px" heights="390px" perrow="1"><br />
Image:ACADO_bioreactor.png| Optimal solution for the ACADO example.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:ACADO | ACADO code]] at [[Bioreactor (ACADO)]]<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]] [[Category: ODE model]] [[Category:Chemical engineering]]<br />
[[Category:Bang bang]]</div>FelixMuellerhttps://mintoc.de/index.php?title=D%27Onofrio_chemotherapy_model&diff=2103D'Onofrio chemotherapy model2016-07-27T08:26:54Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nu = 2<br />
|nc = 4<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 & = & - \zeta x_0 \text{ln} \left( \frac{x_0}{x_1} \right) - F \; x_0 u_1, \\<br />
& \dot{x}_1 & = & b x_0 - \mu x_1 - d x_0^{\frac{2}{3}}x_1 - G u_0 x_1 - \eta x_1 u_1, \\<br />
& \dot{x}_2 & = & u_0, \\<br />
& \dot{x}_3 & = & u_1, \\ [1.5ex]<br />
& u_0 & \in & [0,u_0^{max}],\\<br />
& u_1 & \in & [0,u_1^{max}],\\<br />
& x_2 & \leq & x_2^{max}, \\<br />
& x_3 & \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]] at [[D'Onofrio chemotherapy model (Muscod)]]<br />
<br />
== References ==<br />
<biblist /><br />
<br />
[[Category:MIOCP]]<br />
[[Category:Medicine]]<br />
[[Category:ODE model]]<br />
[[Category:Bang bang]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Catalytic_cracking_problem&diff=2102Catalytic cracking problem2016-07-27T08:26:13Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nz = 2<br />
|np = 3<br />
|nc = 3<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalytic cracking problem tries to determine "reaction coefficients for the catalytic cracking of gas oil into gas and other byproducts." (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{\theta} &\sum\limits_{j=1}^{21} &&||y(\tau_j; \theta) - z_j||^2 \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{y}_1 & = & -(\theta_1 + \theta_3) y_1^2, \\<br />
& \dot{y}_2 & = & \theta_1 y_1^2 - \theta_2 y_2, \\<br />
& \theta_i & \geq & 0.<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
The values <math> z_j </math> are measurements for the concentration for <math> y </math> at time points <math> \tau_1, ..., \tau_{21} </math> and initial conditions are known.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalytic cracking problem (TACO)]]<br />
<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:DAE model]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Catalyst_mixing_problem&diff=2101Catalyst mixing problem2016-07-27T08:25:46Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 1<br />
|nre = 2<br />
}}<!-- Do not insert line break here or Dimensions Box moves up in the layout...<br />
<br />
-->The Catalyst mixing problem seeks an optimal policy for mixing two catalysts "along the length of a tubular plug <br />
ow reactor involving several reactions". (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])<br />
<br />
<br />
<br />
== Mathematical formulation ==<br />
<br />
The problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \min_{x, u} &-1 + x_1(t_f) + x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x}_1 & = & u ( 10 x_2 - x_1), \\<br />
& \dot{x}_2 & = & u ( x_1 - 10 x_2) - (1 - u \, x_2) , \\<br />
& x(t_0) &=& (1, 0)^T, \\<br />
& u(t) &\in& [0,1].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
== Parameters ==<br />
<br />
In this model the parameters used are <math> t_0 = 0, \, \, t_f = 1 </math>.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalyst mixing problem (TACO)]]<br />
<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]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Car_testdrive_(lane_change_manoeuvre)&diff=2100Car testdrive (lane change manoeuvre)2016-07-27T08:25:36Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 7<br />
|nu = 3<br />
|nw = 1<br />
|nc = 1<br />
|nri = 7<br />
}}<br />
<br />
The testdrive control problem is a time optimal double lane change maneouvre with gear shift. It has been introduced as a benchmark problem for mixed-integer optimal control by <bib id="Gerdts2005" />. <br />
<br />
== Mathematical formulation ==<br />
<br />
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. <br />
<br />
The vehicle dynamics are based on a single-track model, derived under the simplifying assumption that rolling and pitching of the car body can be neglected. Consequentially, only a single front and rear wheel is modeled, located in the virtual center of the original two wheels. Motion of the car body is considered on the horizontal plane only.<br />
<br />
Four controls represent the driver's choice on steering and velocity. We denote with <math>w_\delta</math> the steering wheel's angular velocity. The force <math>F_\text{B}</math> controls the total braking force, while the accelerator pedal position <math>\phi</math> is translated into an accelerating force. Finally, the selected gear <math>\mu</math> influences the effective engine torque's transmission.<br />
<br />
<br />
<br />
== Resulting MIOCP ==<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(\cdot), u(\cdot), \mu(\cdot)} & t_\text{f} \\[1.5ex]<br />
\mbox{s.t.} & \dot{x} & = & f(t, x, u, \mu), \\<br />
& x(t_0) &=& x_0, \\<br />
& r(t,x,u) &\geq& 0, \\<br />
& \mu(t) &\in& \{1, 2, 3, 4, 5\}.<br />
\end{array} <br />
</math><br />
<br />
== Parameters ==<br />
<br />
These fixed values are used within the model.<br />
<br />
{| border="1" align="center" cellpadding="5" cellspacing="0"<br />
|- bgcolor=#c7c7c7<br />
! Symbol !! Value !! Unit !! Description<br />
|- <br />
| align=center | <math>m</math> || align=right | 1.239e+3 || kg || Mass of the car<br />
|-<br />
| align=center | <math>g</math> || align=right | 9.81 || m/s^2 || Gravity constant<br />
|-<br />
| align=center | <math>l_\text{f}</math> || align=right | 1.19016 || m || Front wheel distance to center of gravity<br />
|-<br />
| align=center | <math>l_\text{r}</math> || align=right | 1.37484 || m || Rear wheel distance to center of gravity<br />
|-<br />
| align=center | <math>e_\text{SP}</math> || align=right | 0.5 || m || Drag mount point distance to center of gravity<br />
|-<br />
| align=center | <math>R</math> || align=right | 0.302 || m || Wheel radius<br />
|-<br />
| align=center | <math>I_\text{zz}</math> || align=right | 1.752e+3 || kg m^2 || Moment of inertia<br />
|-<br />
| align=center | <math>c_\text{w}</math> || align=right | 0.3 || - || Air drag coefficient<br />
|-<br />
| align=center | <math>\varrho</math> || align=right | 1.249512 || kg/m^3 || Air density<br />
|-<br />
| align=center | <math>A</math> || align=right | 1.4378946874 || m^2 || Effective flow surface<br />
|-<br />
| align=center | <math>i_\text{g}</math> || align=right | 3.09, 2.002, 1.33, 1.0, 0.805 || - || Transmission ratios for the five gears<br />
|-<br />
| align=center | <math>i_\text{t}</math> || align=right | 3.91 || - || Engine transmission ratio<br />
|-<br />
| align=center | <math>B_\text{f}</math> || align=right | 1.096e+1 || - || Pacejka coefficients (stiffness)<br />
|-<br />
| align=center | <math>B_\text{r}</math> || align=right | 1.267e+1 || - || <br />
|-<br />
| align=center | <math>C_\text{f}</math> || align=right | 1.3 || - || Pacejka coefficients (shape)<br />
|-<br />
| align=center | <math>C_\text{r}</math> || align=right | 1.3 || - || <br />
|-<br />
| align=center | <math>D_\text{f}</math> || align=right | 4.5604e+3 || - || Pacejka coefficients (peak)<br />
|-<br />
| align=center | <math>D_\text{r}</math> || align=right | 3.94781e+3 || - || <br />
|-<br />
| align=center | <math>E_\text{f}</math> || align=right | -0.5 || - || Pacejka coefficients (curvature)<br />
|-<br />
| align=center | <math>E_\text{r}</math> || align=right | -0.5 || - || <br />
|}<br />
<br />
== Test course ==<br />
<br />
The double-lane change manoeuvre presented in <bib id="Gerdts2005" /> is realized by constraining the car's position onto a prescribed track at any time <math>t\in[t_0,t_\text{f}]</math>. Starting in the left position with an initial prescribed velocity, the driver is asked to manage a change of lanes modeled by an offset of 3.5 meters in the track. Afterwards he is asked to return to the starting lane. This manoeuvre can be regarded as an overtaking move or as an evasive action taken to avoid hitting an obstacle suddenly appearing on the starting lane.<br />
<br />
From a mathematical point of view, the test track is described by setting up piecewise cubic spline functions <math>P_\text{l}(x)</math> and <math>P_\text{r}(x)</math> modeling the top and bottom track boundary, given a horizontal position <math>x</math>.<br />
<br />
<math><br />
\begin{align}<br />
P_\text{l}(x) &:=& \left\{ <br />
\begin{array}{llrcl}<br />
0 & \text{if } & & x & \leq 44, \\<br />
4\; h_2\; (x-44)^3 & \text{if } & 44 < & x & \leq 44.5, \\<br />
4\; h_2\; (x-45)^3 + h_2 & \text{if } & 44.5 < & x & \leq 45, \\<br />
h_2 & \text{if } & 45 < & x & \leq 70, \\<br />
4\; h_2\; (70-x)^3 + h_2 & \text{if } & 70 < & x & \leq 70.5, \\<br />
4\; h_2\; (71-x)^3 & \text{if } & 70.5 < & x & \leq 71, \\<br />
0 & \text{if } & 71 < & x. & \\<br />
\end{array} \right. \\<br />
P_\text{u}(x) &:=& \left\{ <br />
\begin{array}{llrcl}<br />
h_1 & \text{if } & & x & \leq 15, \\<br />
4\; (h_3-h_1)\; (x-15)^3 + h_1 & \text{if } & 15 < & x & \leq 15.5, \\<br />
4\; (h_3-h_1)\; (x-16)^3 + h_3 & \text{if } & 15.5 < & x & \leq 16, \\<br />
h_3 & \text{if } & 16 < & x & \leq 94, \\<br />
4\; (h_3-h_4)\; (94-x)^3 + h_3 & \text{if } & 94 < & x & \leq 94.5, \\<br />
4\; (h_3-h_4)\; (95-x)^3 + h_4 & \text{if } & 94.5 < & x & \leq 95, \\<br />
h_4 & \text{if } & 95 < & x. & \\<br />
\end{array} \right. <br />
\end{align}<br />
</math><br />
<br />
where <math>B=1.5\;\text{m}</math> is the car's width and<br />
<br />
<math><br />
h_1 := 1.1\; B + 0.25, \quad<br />
h_2 := 3.5, \quad<br />
h_3 := 1.2\; B + 3.75,\quad<br />
h_4 := 1.3\; B + 0.25.<br />
</math><br />
<br />
[[Image:test-course.png|Test course for the double lane change manoeuvre]]<br />
<br />
== Reference Solutions ==<br />
<br />
Reference solutions for the case of a fixed time-grid are given in <bib id="Gerdts2005" />. Solutions for a non-fixed time grid are given in <bib id="Gerdts2006" />.<br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:Muscod | Muscod code]] at [[Car testdrive (lane change manoeuvre) (Muscod)]]<br />
<br />
== Variants ==<br />
<br />
See testdrive [[Car testdrive | overview page]].<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:Minimum time]]<br />
[[Category:Bang bang]]<br />
[[Category:Transport]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Car_testdrive_(elliptic_track)&diff=2099Car testdrive (elliptic track)2016-07-27T08:25:25Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 7<br />
|nu = 3<br />
|nw = 1<br />
|nc = 1<br />
|nri = 7<br />
}}<br />
<br />
The elliptic track testdrive problem is a time optimal periodic control problem with gear shift, first introduced in <bib id="Sager2009a" />. <br />
<br />
== Mathematical formulation ==<br />
<br />
The mathematical equations form a small-scale [[:Category:ODE model|ODE model]] as presented for the [[Car testdrive (lane change manoeuvre) | lane change manoeuvre]].<br />
<br />
The vehicle dynamics are based on a single-track model, derived under the simplifying assumption that rolling and pitching of the car body can be neglected. Consequentially, only a single front and rear wheel is modeled, located in the virtual center of the original two wheels. Motion of the car body is considered on the horizontal plane only.<br />
<br />
Four controls represent the driver's choice on steering and velocity. We denote with <math>w_\delta</math> the steering wheel's angular velocity. The force <math>F_\text{B}</math> controls the total braking force, while the accelerator pedal position <math>\phi</math> is translated into an accelerating force. Finally, the selected gear <math>\mu</math> influences the effective engine torque's transmission.<br />
<br />
== Resulting MIOCP ==<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(\cdot), u(\cdot), \mu(\cdot)} & t_\text{f} \\[1.5ex]<br />
\mbox{s.t.} & \dot{x} & = & f(t, x, u, \mu), \\<br />
& c_\text{x}(t_0) &=& c_\text{x}(t_f), \\<br />
& c_\text{y}(t_0) &=& c_\text{y}(t_f), \\<br />
& v(t_0) &=& v(t_f), \\<br />
& \beta(t_0) &=& \beta(t_f) - 2\pi, \\<br />
& \psi(t_0) &=& \psi(t_f), \\<br />
& \delta(t_0) &=& \delta(t_f), \\<br />
& r(t,x,u) &\geq& 0, \\<br />
& \mu(t) &\in& \{1, 2, 3, 4, 5\}.<br />
\end{array} <br />
</math><br />
<br />
== Variants ==<br />
<br />
See testdrive [[Car testdrive | overview page]].<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:Minimum time]]<br />
[[Category:Bang bang]]<br />
[[Category:Transport]]<br />
[[Category:Periodic]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Batch_reactor&diff=2098Batch reactor2016-07-27T08:24:57Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 2<br />
|nu = 1<br />
|nc = 1<br />
|nre = 2<br />
}}<br />
<br />
This batch reactor problem describes the consecutive reaction of some substance A via substance B into a desired product C.<br />
<br />
The system is interacted with via the control function <math> T(t) </math> which stands for the temperature.<br />
The goal is to produce as much of substance B (which can then be converted into product C) as possible within the time limit.<br />
<br />
== Mathematical formulation ==<br />
<br />
The optimal control problem is given by<br />
<br />
<p><br />
<math><br />
\begin{array}{llcl}<br />
\displaystyle \max_{x, u} & x_2(t_f) \\[1.5ex]<br />
\mbox{s.t.} & \dot{x}_1 & = & -k_1 x_1^2.\\<br />
& \dot{x}_2 & = & k_1 x_1^2 - k_2 x_2,\\<br />
& k_1 & = & 4000 \; e^{(-2500/T(t))}, \\<br />
& k_2 & = & 620000 \; e^{(-5000/T(t))}, \\[1.5ex]<br />
& x(0) &=& (1, 0)^T, \\<br />
& T(t) &\in& [298, 398].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
<math> x_1(t) </math> and <math> x_2(t) </math> represent the concentrations of A and B at timepoint <math> t </math> respectively. The control function <math> T(t) </math> represents the temperature.<br />
<br />
== Parameters ==<br />
<br />
The starting time and end time are given by <math> [t_0, t_f] = [0, 1] </math>.<br />
<br />
== Reference Solutions ==<br />
<br />
This solution was computed using JuMP with a collocation method and 300 discretization points. The differential equations were solved using the explicit Euler Method. The source code can be found at [[Batch reactor (JuMP)]].<br />
<br />
The optimal objective value of the problem is <math> x_2(t_f) = -0.611715 </math>.<br />
<br />
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2"><br />
File:batch_reactor_states_plot.png| Optimal states <math> x_1(t)</math> and <math>x_2(t)</math>.<br />
File:batch_reactor_control_plot.png| Optimal control <math>u(t)</math>.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
* [[:Category: Julia/JuMP | JuMP code]] at [[Batch reactor (JuMP)]]<br />
<br />
== References ==<br />
The problem can be found in [http://www.tomopt.com/docs/TOMLAB_PROPT.pdf the Tomlab PROPT guide] or in [http://www.kirp.chtf.stuba.sk/moodle/mod/resource/view.php?id=5464 the Dynopt guide].<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:Chemical engineering]]<br />
[[Category:Path-constrained arcs]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Bioreactor&diff=2097Bioreactor2016-07-27T08:24:33Z<p>FelixMueller: /* Mathematical formulation */</p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 3<br />
|nu = 1<br />
|nc = 1<br />
|nre = 3<br />
}}<br />
<br />
The bioreactor problem describes an substrate that is converted to a product by the biomass in the reactor. It has three states and a control that is describing the feed concentration of the substrate. The problem is taken from the examples folder of the ACADO toolkit described in:<br />
<bib id="Houska2011a" /><br />
<br />
Houska, Boris, Hans Joachim Ferreau, and Moritz Diehl. <br />
"ACADO toolkit—An open‐source framework for automatic control and dynamic optimization."<br />
Optimal Control Applications and Methods 32.3 (2011): 298-312.<br />
<br />
Originally the problem seems to be motivated by:<br />
<bib id="Versyck1999" /><br />
<br />
VERSYCK, KARINA J., and JAN F. VAN IMPE. <br />
"Feed rate optimization for fed-batch bioreactors: From optimal process performance to optimal parameter estimation."<br />
Chemical Engineering Communications 172.1 (1999): 107-124.<br />
<br />
== Model Formulation ==<br />
<br />
The dynamic model is an [[:Category:ODE model|ODE model]]:<br />
<br />
<p><br />
<math><br />
\begin{array}{rcl}<br />
\dot{X}&=&-DX+\mu X \\<br />
\dot{S}&=& D(S_{f}-S)-\mu /Y_{xs} X \\<br />
\dot{P}&=&-DP+ (\alpha \mu +\beta) X.<br />
\end{array}<br />
</math><br />
</p><br />
<br />
The right-hand side of these equations will be summed up in <math> f(x, S_f) </math>.<br />
<br />
The three states describe the concentration of the biomass (<math>X</math>), the substrate (<math>S</math>), and the product (<math>P</math>) in the reactor. In steady state the feed and outlet are equal and dilute all three concentrations with a ratio <math>D</math>. The biomass grows with a rate<br />
<math>\mu</math>, while it eats up the substrate with the rate <math>\mu/Y_{xs}</math> and produces product at a rate <math>(\alpha \mu +\beta)</math>. The rate <math>\mu</math> is given by:<br />
<br />
<math>\mu = \mu_{m}*(1-P/P_{m})*S/(K_m+S+S^2/K_i)</math><br />
<br />
The fixed parameters (constants) of the model are as follows.<br />
<br />
{| class="wikitable"<br />
|+Parameters<br />
|-<br />
|Name<br />
|Symbol<br />
|Value<br />
|Unit<br />
|-<br />
|Dilution<br />
|<math>D</math><br />
|0.15<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>K_i</math><br />
|22<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>K_m</math><br />
|1.2<br />
|[-]<br />
|-<br />
|Rate coefficient<br />
|<math>P_m</math><br />
|50<br />
|[-]<br />
|-<br />
|Substrate to Biomass rate<br />
|<math>Y_{xs}</math><br />
|0.4<br />
|[-]<br />
|-<br />
|Linear slope<br />
|<math>\alpha</math><br />
|2.2<br />
|[-]<br />
|-<br />
|Linear intercept<br />
|<math>\beta</math><br />
|0.2<br />
|[-]<br />
|-<br />
|Maximal growth rate<br />
|<math>\mu_m</math><br />
|0.48<br />
|[-]<br />
|-<br />
|}<br />
<br />
== Mathematical formulation ==<br />
<br />
Writing shortly for the states in vector notation <math>x=(X,S,P)^T</math> the OCP reads:<br />
<br />
<p><br />
<math><br />
\begin{array}{clcl}<br />
\displaystyle \min_{x,S_f} & J(x,S_f)\\[1.5ex]<br />
\mbox{s.t.} <br />
& \dot{x} & = & f(x,S_f)\\<br />
& x(0) & = & (6.5,12,22)^T \\<br />
& S_f & \in &[28.7,40].<br />
\end{array} <br />
</math><br />
</p><br />
<br />
=== Objective ===<br />
<p><br />
<math><br />
J(x,S_f)=\int_0^{48}D(S_f-P)^2dt<br />
</math><br />
</p><br />
<br />
== Reference Solution ==<br />
<br />
Here we present the reference solution of the reimplemented example in the ACADO code generation with matlab. The source code is given in the next section.<br />
<br />
<gallery caption="Reference solution" widths="551px" heights="390px" perrow="1"><br />
Image:ACADO_bioreactor.png| Optimal solution for the ACADO example.<br />
</gallery><br />
<br />
== Source Code ==<br />
<br />
Model descriptions are available in<br />
<br />
* [[:Category:ACADO | ACADO code]] at [[Bioreactor (ACADO)]]<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]] [[Category: ODE model]] [[Category:Chemical engineering]]<br />
[[Category:Bang bang]]</div>FelixMuellerhttps://mintoc.de/index.php?title=Annihilation_of_calcium_oscillations_with_PLC_activation_inhibition&diff=2096Annihilation of calcium oscillations with PLC activation inhibition2016-07-27T08:22:48Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 2<br />
|nc = 1<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>FelixMuellerhttps://mintoc.de/index.php?title=Annihilation_of_calcium_oscillations&diff=2095Annihilation of calcium oscillations2016-07-27T08:22:39Z<p>FelixMueller: </p>
<hr />
<div>{{Dimensions<br />
|nd = 1<br />
|nx = 4<br />
|nw = 1<br />
|nc = 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:JModelica | JModelica code]] at [[Annihilation of calcium oscillations (jModelica)]]<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>FelixMuellerhttps://mintoc.de/index.php?title=Main_Page&diff=2094Main Page2016-07-27T08:04:18Z<p>FelixMueller: /* Direct links to Problems (add) */</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]] - [[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>FelixMuellerhttps://mintoc.de/index.php?title=Template:Current_News&diff=2093Template:Current News2016-07-27T08:03:40Z<p>FelixMueller: </p>
<hr />
<div>{| style="background:#EEEEFF;"<br />
<!-- Add 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 />
<!-- 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/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>FelixMuellerhttps://mintoc.de/index.php?title=Parameters_for_Industrial_Robot&diff=2092Parameters for Industrial Robot2016-07-27T08:02:51Z<p>FelixMueller: </p>
<hr />
<div>This page contains a matlab script for the implicit formulation of the parameters for the [[Industrial robot]] problem.<br />
<br />
<source lang="Matlab"><br />
%initialization<br />
<br />
dq=zeros(3,1);<br />
dotq=zeros(3,1);<br />
ddq=zeros(3,1);<br />
<br />
q=y(1:3);<br />
dq(1:3,1)=y(4:6);<br />
dotq(1:3,1)=dy(1:3);<br />
ddq(1:3,1)=dy(4:6);<br />
<br />
% intermediate states<br />
<br />
cosq=zeros(3,1);<br />
<br />
cosq(1,1) = cos(q(1));<br />
cosq(2,1) = cos(q(2));<br />
cosq(3,1) = cos(q(3));<br />
<br />
sinq=zeros(3,1);<br />
<br />
sinq(1,1) = sin(q(1));<br />
sinq(2,1) = sin(q(2));<br />
sinq(3,1) = sin(q(3));<br />
<br />
cosq2q3 = cos(q(2) - q(3));<br />
<br />
sinq1q2= sin(q(1) - q(2));<br />
sinq2q3= sin(q(2) - q(3));<br />
<br />
sin2q2= sin(2*q(2));<br />
<br />
sin2q22q3= sin(2*q(2) - 2*q(3));<br />
sin2q2q3= sin(2*q(2) - q(3));<br />
<br />
cos2q2= cos(2*q(2));<br />
<br />
% Modell parameters<br />
<br />
b_1 = 2.6491674286195;<br />
b_2 = 0.0168799610481578;<br />
b_3 = 1.43009916890323;<br />
b_4 = 0.53524835732352;<br />
b_5 = 0.791116620901176;<br />
b_6 = 0.626240578068518;<br />
b_7 = 1.57199673626875;<br />
b_8 = 0.0175195488386143;<br />
b_9 = 0.0204978721411787;<br />
b_10= 1.07049671464704;<br />
b_11= 1.81771647041588;<br />
b_12= 0.422954084243519;<br />
<br />
b=zeros(3,3); % Inerta Matrix of Robot Modell (3DOF)<br />
<br />
b(1,1) = b_1*(sinq(2)^2) - b_2*sin2q2 + b_3*(cosq(2)^2) + b_4*cosq(3) + b_5*(cosq2q3^2) - b_6*cos(2*q(2) - 2*q(3)) - b_4*cos(2*q(2) - q(3)) - b_7;<br />
b(2,1) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(3,1) = -b_9*cosq2q3;<br />
b(1,2) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(2,2) = b_10*cosq(3) + b_11;<br />
b(3,2) = -b_4*cosq(3) - b_12;<br />
b(1,3) = -b_9*cosq2q3;<br />
b(2,3) = -b_4*cosq(3) - b_12;<br />
b(3,3) = b_12;<br />
<br />
g=zeros(3,1); % gravity torque vector<br />
<br />
g_1 = 32.4854155934467;<br />
g_2 = 13.1269659633593;<br />
g_3 = 0.41398104470607;<br />
<br />
g(1) = 0;<br />
g(2) =-g_1*sinq(2) - g_2*sinq2q3 + g_3*cosq(2);<br />
g(3) = g_2*sinq2q3;<br />
<br />
f=zeros(3,1); % friction torque vector<br />
<br />
f(1) = 0.186461043875677*atan(100*dq(1));<br />
f(2) = 0.184489618216264*atan(100*dq(2));<br />
f(3) = 0.156066608371045*atan(100*dq(3));<br />
<br />
c=zeros(3,3); % coriolis matrix<br />
<br />
c_1 = 0.609534129858137;<br />
c_2 = 0.23068226761793;<br />
c_3 = b_4;<br />
c_4 = b_2;<br />
c_5 = 0.26762417866176;<br />
c_6 = 152.383532464534;<br />
c_7 = 57.6705669044826;<br />
c_8 = 133.81208933088;<br />
c_9 = 4.21999026203945;<br />
c_10= 66.90604466544;<br />
c_11= b_8;<br />
c_12= b_9;<br />
c_13= 1.3381208933088;<br />
<br />
c(1,1) = c_1*dq(2)*sin2q2 + c_2*dq(2)*sin2q22q3 + c_3*dq(2)*sin2q2q3 - c_4*dq(2)*cos2q2 - c_5*dq(3)*sinq(3) - c_2*dq(3)*sin2q22q3 - c_5*dq(3)*sin2q2q3;<br />
c(2,1) = -(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250;<br />
c(3,1) = -(-1)*(-1)*dq(1)*(-c_10*sinq(3) - c_7*sin2q22q3 - c_10*sin2q2q3)/250;<br />
c(1,2) = -(-1)*(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250 + dq(2)*(-c_11*sinq(2) - c_12*sinq2q3) + c_12*dq(3)*sinq2q3;<br />
c(2,2) = -c_3*dq(3)*sinq(3);<br />
c(3,2) = c_3*dq(2)*sinq(3);<br />
c(1,3) = -c_5*dq(1)*sinq(3) - c_2*dq(1)*sin2q22q3 - c_5*dq(1)*sin2q2q3 + c_12*dq(2)*sinq2q3 - c_12*dq(3)*sinq2q3;<br />
c(2,3) = c_13*(-(-1)*(-2)*dq(2)/5 - (-1)*2*dq(3)/5)*sinq(3);<br />
c(3,3) = 0; <br />
<br />
</source></div>FelixMuellerhttps://mintoc.de/index.php?title=Parameters_for_Industrial_Robot&diff=2091Parameters for Industrial Robot2016-07-26T15:54:09Z<p>FelixMueller: </p>
<hr />
<div>This page contains a matlab script for the implicit formulation of the parameters for the [[Industrial robot]] problem.<br />
<br />
<source lang="Matlab"><br />
%initialization<br />
<br />
dq=zeros(3,1);<br />
dotq=zeros(3,1);<br />
ddq=zeros(3,1);<br />
<br />
q=y(1:3);<br />
dq(1:3,1)=y(4:6);<br />
dotq(1:3,1)=dy(1:3);<br />
ddq(1:3,1)=dy(4:6);<br />
<br />
% intermediate states<br />
<br />
cosq=zeros(3,1);<br />
<br />
cosq(1,1) = cos(q(1));<br />
cosq(2,1) = cos(q(2));<br />
cosq(3,1) = cos(q(3));<br />
<br />
sinq=zeros(3,1);<br />
<br />
sinq(1,1) = sin(q(1));<br />
sinq(2,1) = sin(q(2));<br />
sinq(3,1) = sin(q(3));<br />
<br />
cosq2q3 = cos(q(2) - q(3));<br />
<br />
sinq1q2= sin(q(1) - q(2));<br />
sinq2q3= sin(q(2) - q(3));<br />
<br />
sin2q2= sin(2*q(2));<br />
<br />
sin2q22q3= sin(2*q(2) - 2*q(3));<br />
sin2q2q3= sin(2*q(2) - q(3));<br />
<br />
cos2q2= cos(2*q(2));<br />
<br />
% Modell parameters<br />
<br />
b_1 = 2.6491674286195;<br />
b_2 = 0.0168799610481578;<br />
b_3 = 1.43009916890323;<br />
b_4 = 0.53524835732352;<br />
b_5 = 0.791116620901176;<br />
b_6 = 0.626240578068518;<br />
b_7 = 1.57199673626875;<br />
b_8 = 0.0175195488386143;<br />
b_9 = 0.0204978721411787;<br />
b_10= 1.07049671464704;<br />
b_11= 1.81771647041588;<br />
b_12= 0.422954084243519;<br />
<br />
b=zeros(3,3); % Inerta Matrix of Robot Modell (3DOF)<br />
<br />
b(1,1) = b_1*(sinq(2)^2) - b_2*sin2q2 + b_3*(cosq(2)^2) + b_4*cosq(3) + b_5*(cosq2q3^2) - b_6*cos(2*q(2) - 2*q(3)) - b_4*cos(2*q(2) - q(3)) - b_7;<br />
b(2,1) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(3,1) = -b_9*cosq2q3;<br />
b(1,2) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(2,2) = b_10*cosq(3) + b_11;<br />
b(3,2) = -b_4*cosq(3) - b_12;<br />
b(1,3) = -b_9*cosq2q3;<br />
b(2,3) = -b_4*cosq(3) - b_12;<br />
b(3,3) = b_12;<br />
<br />
g=zeros(3,1); % gravity torque vector<br />
<br />
g_1 = 32.4854155934467;<br />
g_2 = 13.1269659633593;<br />
g_3 = 0.41398104470607;<br />
<br />
%g(1) = 0;<br />
%g(2) =-g_1*sinq(2) - g_2*sinq2q3 + g_3*cosq(2);<br />
%g(3) = g_2*sinq2q3;<br />
<br />
f=zeros(3,1); % friction torque vector<br />
<br />
f(1) = 0.186461043875677*atan(100*dq(1));<br />
f(2) = 0.184489618216264*atan(100*dq(2));<br />
f(3) = 0.156066608371045*atan(100*dq(3));<br />
<br />
c=zeros(3,3); % coriolis matrix<br />
<br />
c_1 = 0.609534129858137;<br />
c_2 = 0.23068226761793;<br />
c_3 = b_4;<br />
c_4 = b_2;<br />
c_5 = 0.26762417866176;<br />
c_6 = 152.383532464534;<br />
c_7 = 57.6705669044826;<br />
c_8 = 133.81208933088;<br />
c_9 = 4.21999026203945;<br />
c_10= 66.90604466544;<br />
c_11= b_8;<br />
c_12= b_9;<br />
c_13= 1.3381208933088;<br />
<br />
c(1,1) = c_1*dq(2)*sin2q2 + c_2*dq(2)*sin2q22q3 + c_3*dq(2)*sin2q2q3 - c_4*dq(2)*cos2q2 - c_5*dq(3)*sinq(3) - c_2*dq(3)*sin2q22q3 - c_5*dq(3)*sin2q2q3;<br />
c(2,1) = -(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250;<br />
c(3,1) = -(-1)*(-1)*dq(1)*(-c_10*sinq(3) - c_7*sin2q22q3 - c_10*sin2q2q3)/250;<br />
c(1,2) = -(-1)*(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250 + dq(2)*(-c_11*sinq(2) - c_12*sinq2q3) + c_12*dq(3)*sinq2q3;<br />
c(2,2) = -c_3*dq(3)*sinq(3);<br />
c(3,2) = c_3*dq(2)*sinq(3);<br />
c(1,3) = -c_5*dq(1)*sinq(3) - c_2*dq(1)*sin2q22q3 - c_5*dq(1)*sin2q2q3 + c_12*dq(2)*sinq2q3 - c_12*dq(3)*sinq2q3;<br />
c(2,3) = c_13*(-(-1)*(-2)*dq(2)/5 - (-1)*2*dq(3)/5)*sinq(3);<br />
c(3,3) = 0; <br />
<br />
</source></div>FelixMuellerhttps://mintoc.de/index.php?title=Parameters_for_Industrial_Robot&diff=2090Parameters for Industrial Robot2016-07-26T15:54:02Z<p>FelixMueller: Created page with "This page contains a matlab script for the implicit formulation of the parameters for the Industrial robot problem. <source lang="Matlab"> %initialization dq=zeros..."</p>
<hr />
<div>This page contains a matlab script for the implicit formulation of the parameters for the [[Industrial robot]] problem.<br />
<br />
<source lang="Matlab"><br />
%initialization<br />
<br />
dq=zeros(3,1);<br />
dotq=zeros(3,1);<br />
ddq=zeros(3,1);<br />
<br />
q=y(1:3);<br />
dq(1:3,1)=y(4:6);<br />
dotq(1:3,1)=dy(1:3);<br />
ddq(1:3,1)=dy(4:6);<br />
<br />
% intermediate states<br />
<br />
cosq=zeros(3,1);<br />
<br />
cosq(1,1) = cos(q(1));<br />
cosq(2,1) = cos(q(2));<br />
cosq(3,1) = cos(q(3));<br />
<br />
sinq=zeros(3,1);<br />
<br />
sinq(1,1) = sin(q(1));<br />
sinq(2,1) = sin(q(2));<br />
sinq(3,1) = sin(q(3));<br />
<br />
cosq2q3 = cos(q(2) - q(3));<br />
<br />
sinq1q2= sin(q(1) - q(2));<br />
sinq2q3= sin(q(2) - q(3));<br />
<br />
sin2q2= sin(2*q(2));<br />
<br />
sin2q22q3= sin(2*q(2) - 2*q(3));<br />
sin2q2q3= sin(2*q(2) - q(3));<br />
<br />
cos2q2= cos(2*q(2));<br />
<br />
% Modell parameters<br />
<br />
b_1 = 2.6491674286195;<br />
b_2 = 0.0168799610481578;<br />
b_3 = 1.43009916890323;<br />
b_4 = 0.53524835732352;<br />
b_5 = 0.791116620901176;<br />
b_6 = 0.626240578068518;<br />
b_7 = 1.57199673626875;<br />
b_8 = 0.0175195488386143;<br />
b_9 = 0.0204978721411787;<br />
b_10= 1.07049671464704;<br />
b_11= 1.81771647041588;<br />
b_12= 0.422954084243519;<br />
<br />
b=zeros(3,3); % Inerta Matrix of Robot Modell (3DOF)<br />
<br />
b(1,1) = b_1*(sinq(2)^2) - b_2*sin2q2 + b_3*(cosq(2)^2) + b_4*cosq(3) + b_5*(cosq2q3^2) - b_6*cos(2*q(2) - 2*q(3)) - b_4*cos(2*q(2) - q(3)) - b_7;<br />
b(2,1) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(3,1) = -b_9*cosq2q3;<br />
b(1,2) = b_8*cosq(2) + b_9*cosq2q3;<br />
b(2,2) = b_10*cosq(3) + b_11;<br />
b(3,2) = -b_4*cosq(3) - b_12;<br />
b(1,3) = -b_9*cosq2q3;<br />
b(2,3) = -b_4*cosq(3) - b_12;<br />
b(3,3) = b_12;<br />
<br />
g=zeros(3,1); % gravity torque vector<br />
<br />
g_1 = 32.4854155934467;<br />
g_2 = 13.1269659633593;<br />
g_3 = 0.41398104470607;<br />
<br />
%g(1) = 0;<br />
%g(2) =-g_1*sinq(2) - g_2*sinq2q3 + g_3*cosq(2);<br />
%g(3) = g_2*sinq2q3;<br />
<br />
f=zeros(3,1); % friction torque vector<br />
<br />
f(1) = 0.186461043875677*atan(100*dq(1));<br />
f(2) = 0.184489618216264*atan(100*dq(2));<br />
f(3) = 0.156066608371045*atan(100*dq(3));<br />
<br />
c=zeros(3,3); % coriolis matrix<br />
<br />
c_1 = 0.609534129858137;<br />
c_2 = 0.23068226761793;<br />
c_3 = b_4;<br />
c_4 = b_2;<br />
c_5 = 0.26762417866176;<br />
c_6 = 152.383532464534;<br />
c_7 = 57.6705669044826;<br />
c_8 = 133.81208933088;<br />
c_9 = 4.21999026203945;<br />
c_10= 66.90604466544;<br />
c_11= b_8;<br />
c_12= b_9;<br />
c_13= 1.3381208933088;<br />
<br />
c(1,1) = c_1*dq(2)*sin2q2 + c_2*dq(2)*sin2q22q3 + c_3*dq(2)*sin2q2q3 - c_4*dq(2)*cos2q2 - c_5*dq(3)*sinq(3) - c_2*dq(3)*sin2q22q3 - c_5*dq(3)*sin2q2q3;<br />
c(2,1) = -(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250;<br />
c(3,1) = -(-1)*(-1)*dq(1)*(-c_10*sinq(3) - c_7*sin2q22q3 - c_10*sin2q2q3)/250;<br />
c(1,2) = -(-1)*(-1)*dq(1)*(-c_6*sin2q2 - c_7*sin2q22q3 - c_8*sin2q2q3 + c_9*cos2q2)/250 + dq(2)*(-c_11*sinq(2) - c_12*sinq2q3) + c_12*dq(3)*sinq2q3;<br />
c(2,2) = -c_3*dq(3)*sinq(3);<br />
c(3,2) = c_3*dq(2)*sinq(3);<br />
c(1,3) = -c_5*dq(1)*sinq(3) - c_2*dq(1)*sin2q22q3 - c_5*dq(1)*sin2q2q3 + c_12*dq(2)*sinq2q3 - c_12*dq(3)*sinq2q3;<br />
c(2,3) = c_13*(-(-1)*(-2)*dq(2)/5 - (-1)*2*dq(3)/5)*sinq(3);<br />
c(3,3) = 0; <br />
<br />
</source></div>FelixMueller