mintOC - User contributions [en]

Supermarket refrigeration system

2009-08-11T14:35:21Z

PaulineBoizet: /* Optimica */

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

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

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

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

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

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

== Mathematical formulation ==

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


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



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

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

\end{array}
</math>


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

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

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

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


<math>
\begin{array}{rcl}

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


== Parameters ==

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

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

== Reference Solutions ==

For the relaxed problem (we only demand <math>u_i(t) \in [0,1]</math> instead of <math>u_i(t) \in \{0,1\}</math> the optimal solution is 12072.45.
If we restrict ourselves to a solution with only integer controls we obtain the optimum objective value 12252.81.

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
// number of compressors
#define NVALVES 2

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

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

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

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

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

double m_in_suc = 0.0;

int i;

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

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

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

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

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

// temperatures:

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

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

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

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

</source>
=== Optimica ===
<source lang="optimica">
package Supermarket_pack

optimization Supermarket_Opt (objective = cost(tf),
startTime = 0,
finalTime=1)

// Parameters
parameter Real Q_airload = 3000.00;
parameter Real m_ref = 0.20;
parameter Real M_goods=200.00;
parameter Real Cp_goods=1000.00;
parameter Real UA_goods=300.00;
parameter Real M_wall=260.00;
parameter Real Cp_wall=385.00;
parameter Real UA_air=500.00;
parameter Real M_air=50.00;
parameter Real Cp_air=1000.00;
parameter Real UA_wall=4000.00;
parameter Real t_fill=40.00;
parameter Real T_sh=10.00;
parameter Real M_ref=1.00;
parameter Real V_suc=5.00;
parameter Real V_sl=0.08;
parameter Real eta_vol=0.81;
parameter Real tf(free=true,min=650,max=750)=650;

// The states
Real x0;
Real x1;
Real x2;
Real x3;
Real x4;
Real x5;
Real x6;
Real x7;
Real x8;
Real cost(start=0);

// The control signal
input Real u0(free=true,min=0,max=1);
input Real u1(free=true,min=0,max=1);
input Real u2(free=true,min=0,max=1);

//Te=(-4.3544*x0*x0+29.224*x0-51.2005)
//delta_h=((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)
//rho_suc=(4.6073*x0+0.3798)
//drho=(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)
//f=((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)

equation
der(x0) = ((1.0/(V_suc*(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)))*((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*(x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005))+x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))+M_ref-eta_vol*V_sl*0.5*(u2)*(4.6073*x0+0.3798)))*tf;
der(x1) =(-UA_goods*(x1-x3)/M_goods*Cp_goods)*tf;
der(x2)=((UA_air*(x3-x2)-(UA_wall/M_ref)*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x3)=((UA_goods*(x1-x3)+Q_airload-UA_air*(x3-x2))/(M_air*Cp_air))*tf;
der(x4)=(((M_ref-x4)/t_fill)*u0-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u0))*tf;
der(x5)=(-UA_goods*(x5-x7)/(M_goods*Cp_goods))*tf;
der(x6)=((UA_air*(x7-x6)-(UA_wall/M_ref)*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x7)=((UA_goods*(x5-x7)+Q_airload-UA_air*(x7-x6))/(M_air*Cp_air))*tf;
der(x8)=(((M_ref-x8)/t_fill)*u1-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u1))*tf;
der(cost) = ((-1.0/((1/*time*/)))*cost+(1.0/((1/*time*/)))*((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;
// der(cost) = 1.0 / 650 * (((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;

constraint
x3>=2.0;
x3<=5.0;
x7>=2.0;
x7<=5.0;
x0<=1.7;
x0(tf)=x0(startTime);
x1(tf)=x1(startTime);
x2(tf)=x2(startTime);
x3(tf)=x3(startTime);
x4(tf)=x4(startTime);
x5(tf)=x5(startTime);
x6(tf)=x6(startTime);
x7(tf)=x7(startTime);
x8(tf)=x8(startTime);

end Supermarket_Opt;

end Supermarket_pack;
</source>

== Variants ==

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

In the paper <bibref>Larsen2007</bibref> mentioned above, the problem was stated slightly different:

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

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

\end{cases}
</math>

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

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

== References ==
<bibreferences/>


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

Supermarket refrigeration system

2009-08-11T14:28:53Z

PaulineBoizet: /* Source Code */

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

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

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

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

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

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

== Mathematical formulation ==

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


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



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

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

\end{array}
</math>


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

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

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

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


<math>
\begin{array}{rcl}

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


== Parameters ==

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

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

== Reference Solutions ==

For the relaxed problem (we only demand <math>u_i(t) \in [0,1]</math> instead of <math>u_i(t) \in \{0,1\}</math> the optimal solution is 12072.45.
If we restrict ourselves to a solution with only integer controls we obtain the optimum objective value 12252.81.

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
// number of compressors
#define NVALVES 2

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

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

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

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

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

double m_in_suc = 0.0;

int i;

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

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

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

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

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

// temperatures:

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

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

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

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

</source>
=== Optimica ===
<source lang="optimica">
package Supermarket_pack

optimization Supermarket_Opt (objective = cost(tf),
startTime = 0,
finalTime=1)

// Parameters
parameter Real Q_airload = 3000.00;
parameter Real m_ref = 0.20;
parameter Real M_goods=200.00;
parameter Real Cp_goods=1000.00;
parameter Real UA_goods=300.00;
parameter Real M_wall=260.00;
parameter Real Cp_wall=385.00;
parameter Real UA_air=500.00;
parameter Real M_air=50.00;
parameter Real Cp_air=1000.00;
parameter Real UA_wall=4000.00;
parameter Real t_fill=40.00;
parameter Real T_sh=10.00;
parameter Real M_ref=1.00;
parameter Real V_suc=5.00;
parameter Real V_sl=0.08;
parameter Real eta_vol=0.81;
parameter Real tf(free=true,min=650,max=750)=650;

// The states
Real x0;
Real x1;
Real x2;
Real x3;
Real x4;
Real x5;
Real x6;
Real x7;
Real x8;
Real cost(start=0);

// The control signal
input Real u0(free=true,min=0,max=1);
input Real u1(free=true,min=0,max=1);
input Real u2(free=true,min=0,max=1);

//Te=(-4.3544*x0*x0+29.224*x0-51.2005)
//delta_h=((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)
//rho_suc=(4.6073*x0+0.3798)
//drho=(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)
//f=((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)

equation
der(x0) = ((1.0/(V_suc*(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)))*((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*(x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005))+x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))+M_ref-eta_vol*V_sl*0.5*(u2)*(4.6073*x0+0.3798)))*tf;
der(x1) =(-UA_goods*(x1-x3)/M_goods*Cp_goods)*tf;
der(x2)=((UA_air*(x3-x2)-(UA_wall/M_ref)*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x3)=((UA_goods*(x1-x3)+Q_airload-UA_air*(x3-x2))/(M_air*Cp_air))*tf;
der(x4)=(((M_ref-x4)/t_fill)*u0-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u0))*tf;
der(x5)=(-UA_goods*(x5-x7)/(M_goods*Cp_goods))*tf;
der(x6)=((UA_air*(x7-x6)-(UA_wall/M_ref)*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x7)=((UA_goods*(x5-x7)+Q_airload-UA_air*(x7-x6))/(M_air*Cp_air))*tf;
der(x8)=(((M_ref-x8)/t_fill)*u1-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u1))*tf;
// der(cost) = ((-1.0/((1.0/*time*/)))*cost+(1.0/((1.0/*time*/)))*((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;
der(cost) = 1.0 / 650 * (((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;

constraint
x3>=2.0;
x3<=5.0;
x7>=2.0;
x7<=5.0;
x0<=1.7;
x0(tf)=x0(startTime);
x1(tf)=x1(startTime);
x2(tf)=x2(startTime);
x3(tf)=x3(startTime);
x4(tf)=x4(startTime);
x5(tf)=x5(startTime);
x6(tf)=x6(startTime);
x7(tf)=x7(startTime);
x8(tf)=x8(startTime);

end Supermarket_Opt;

end Supermarket_pack;
</source>

== Variants ==

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

In the paper <bibref>Larsen2007</bibref> mentioned above, the problem was stated slightly different:

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

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

\end{cases}
</math>

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

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

== References ==
<bibreferences/>


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

Supermarket refrigeration system

2009-08-11T14:25:44Z

PaulineBoizet: /* Source Code */

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

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

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

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

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

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

== Mathematical formulation ==

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


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



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

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

\end{array}
</math>


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

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

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

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


<math>
\begin{array}{rcl}

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


== Parameters ==

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

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

== Reference Solutions ==

For the relaxed problem (we only demand <math>u_i(t) \in [0,1]</math> instead of <math>u_i(t) \in \{0,1\}</math> the optimal solution is 12072.45.
If we restrict ourselves to a solution with only integer controls we obtain the optimum objective value 12252.81.

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
// number of compressors
#define NVALVES 2

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

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

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

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

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

double m_in_suc = 0.0;

int i;

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

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

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

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

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

// temperatures:

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

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

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

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

=== Optimica ===
package Supermarket_pack

optimization Supermarket_Opt (objective = cost(tf),
startTime = 0,
finalTime=1)

// Parameters
parameter Real Q_airload = 3000.00;
parameter Real m_ref = 0.20;
parameter Real M_goods=200.00;
parameter Real Cp_goods=1000.00;
parameter Real UA_goods=300.00;
parameter Real M_wall=260.00;
parameter Real Cp_wall=385.00;
parameter Real UA_air=500.00;
parameter Real M_air=50.00;
parameter Real Cp_air=1000.00;
parameter Real UA_wall=4000.00;
parameter Real t_fill=40.00;
parameter Real T_sh=10.00;
parameter Real M_ref=1.00;
parameter Real V_suc=5.00;
parameter Real V_sl=0.08;
parameter Real eta_vol=0.81;
parameter Real tf(free=true,min=650,max=750)=650;

// The states
Real x0;
Real x1;
Real x2;
Real x3;
Real x4;
Real x5;
Real x6;
Real x7;
Real x8;
Real cost(start=0);

// The control signal
input Real u0(free=true,min=0,max=1);
input Real u1(free=true,min=0,max=1);
input Real u2(free=true,min=0,max=1);

//Te=(-4.3544*x0*x0+29.224*x0-51.2005)
//delta_h=((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)
//rho_suc=(4.6073*x0+0.3798)
//drho=(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)
//f=((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)

equation
der(x0) = ((1.0/(V_suc*(-0.0329*x0*x0*x0+0.2161*x0*x0-0.4742*x0+5.4817)))*((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*(x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005))+x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))+M_ref-eta_vol*V_sl*0.5*(u2)*(4.6073*x0+0.3798)))*tf;
der(x1) =(-UA_goods*(x1-x3)/M_goods*Cp_goods)*tf;
der(x2)=((UA_air*(x3-x2)-(UA_wall/M_ref)*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x3)=((UA_goods*(x1-x3)+Q_airload-UA_air*(x3-x2))/(M_air*Cp_air))*tf;
der(x4)=(((M_ref-x4)/t_fill)*u0-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x4*(x2-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u0))*tf;
der(x5)=(-UA_goods*(x5-x7)/(M_goods*Cp_goods))*tf;
der(x6)=((UA_air*(x7-x6)-(UA_wall/M_ref)*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))/(M_wall*Cp_wall))*tf;
der(x7)=((UA_goods*(x5-x7)+Q_airload-UA_air*(x7-x6))/(M_air*Cp_air))*tf;
der(x8)=(((M_ref-x8)/t_fill)*u1-((UA_wall/(M_ref*((0.0217*x0*x0-0.1704*x0+2.2988)*10^5)))*x8*(x6-(-4.3544*x0*x0+29.224*x0-51.2005)))*(1-u1))*tf;
// der(cost) = ((-1.0/((1.0/*time*/)))*cost+(1.0/((1.0/*time*/)))*((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;
der(cost) = 1.0 / 650 * (((u2)*0.5*eta_vol*V_sl*((0.0265*x0*x0*x0-0.4346*x0*x0+2.4923*x0+1.2189)*10^5)))*tf;

constraint
x3>=2.0;
x3<=5.0;
x7>=2.0;
x7<=5.0;
x0<=1.7;
x0(tf)=x0(startTime);
x1(tf)=x1(startTime);
x2(tf)=x2(startTime);
x3(tf)=x3(startTime);
x4(tf)=x4(startTime);
x5(tf)=x5(startTime);
x6(tf)=x6(startTime);
x7(tf)=x7(startTime);
x8(tf)=x8(startTime);

end Supermarket_Opt;

end Supermarket_pack;
</source>

== Variants ==

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

In the paper <bibref>Larsen2007</bibref> mentioned above, the problem was stated slightly different:

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

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

\end{cases}
</math>

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

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

== References ==
<bibreferences/>


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

Annihilation of calcium oscillations

2009-08-11T14:18:30Z

PaulineBoizet: /* Optimica */

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

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

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

== Biological motivation ==

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

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

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

== Mathematical formulation ==

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

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

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

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

== Parameters ==

These fixed values are used within the model.

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

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

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

== Optimization ==

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

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

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
rhs[0] = p[0]
+ p[1]*xd[0]
- p[2]*xd[0]*xd[1]/(xd[0]+p[3])
- p[4]*xd[0]*xd[2]/(xd[0]+p[5]);

rhs[1] = p[6]*xd[0]
- p[7]*xd[1]/(xd[1]+p[8]);

rhs[2] = p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])
+ p[18]*p[11]*xd[1]
+ p[12]*xd[0]
- p[13]*xd[2]/(u[0]*xd[2]+p[14])
- p[15]*xd[2]/(xd[2]+p[16])
+ 0.1*xd[3];

rhs[3] = -p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])
+ p[15]*xd[2]/(xd[2]+p[16])
- 0.1*xd[3];
</source>

=== Optimica ===

The optimica code of the problem:
<source lang="optimica">

package Calcium_pack

optimization Calcium_Opt (objective = cost(finalTime),
startTime = 0,
finalTime = 22)

// Parameters
parameter Real k1 = 0.09; // Parameter k1
parameter Real k2 = 2.30066; // Parameter k2
parameter Real k3 = 0.64; // Parameter k3
parameter Real K4 = 0.19; // Parameter K4
parameter Real k5 = 4.88; // Parameter k5
parameter Real K6 = 1.18; // Parameter K6
parameter Real k7 = 2.08; // Parameter k7
parameter Real k8 = 32.24; // Parameter k8
parameter Real K9 = 29.09; // Parameter K9
parameter Real k10 = 5.0; // Parameter k10
parameter Real K11 = 2.67; // Parameter K11
parameter Real k12 = 0.7; // Parameter k12
parameter Real k13 = 13.58; // Parameter k13
parameter Real k14 = 153.0; // Parameter k14
parameter Real K15 = 0.16; // Parameter K15
parameter Real k16 = 4.85; // Parameter k16
parameter Real K17 = 0.05; // Parameter K17
parameter Real p1 = 100; // Parameter p1
parameter Real x0tilde = 6.78677; // Parameter x0tilde
parameter Real x1tilde = 22.65836; // Parameter x1tilde
parameter Real x2tilde = 0.384306; // Parameter x2tilde
parameter Real x3tilde = 0.28977; // Parameter x3tilde

// The states
Real x0(start=0.03966);
Real x1(start=1.09799);
Real x2(start=0.00142);
Real x3(start=1.65431);

// The control signal
input Real u(free=true);
//parameter Real umax(free=true,min=1.1,max=1.3);
//parameter Real umax=1.3;

Real cost(start=0);

equation
der(x0) = k1+k2*x0-((k3*x0*x1)/(8*x0+K4))-((k5*x0*x2)/(x0+K6));
der(x1) = k7*x0 - ((k8*x1)/(x1+K9));
der(x2)=(k10*x1*x2*x3)/(x3+K11)+k12*x1+k13*x0-(k14*x2)/(u*x2+K15)-(k16*x2)/(x2+K17)+x3/10;
der(x3)=-(k10*x1*x2*x3)/(x3+K11)+(k16*x2)/(x2+K17)-x3/10;
der(cost)=sqrt((x0-x0tilde)*(x0-x0tilde)+(x1-x1tilde)*(x1-x1tilde)+
(x2-x2tilde)*(x2-x2tilde)+(x3-x3tilde)*(x3-x3tilde))+p1*u;
constraint
x0>=0;
x1>=0;
x2>=0;
x3>=0;
// u={0,1.1,1.3};
u>=1;
u<=1.3;
end Calcium_Opt;

end Calcium_pack;

</source>

== Variants ==

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

== References ==
<bibreferences/>


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

Supermarket refrigeration system

2009-08-11T13:59:34Z

PaulineBoizet: Removed control u_3 from model formulation

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

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

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

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

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

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

== Mathematical formulation ==

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


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



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

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

\end{array}
</math>


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

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

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

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


<math>
\begin{array}{rcl}

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


== Parameters ==

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

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

== Reference Solutions ==

For the relaxed problem (we only demand <math>u_i(t) \in [0,1]</math> instead of <math>u_i(t) \in \{0,1\}</math> the optimal solution is 12072.45.
If we restrict ourselves to a solution with only integer controls we obtain the optimum objective value 12252.81.

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
// number of compressors
#define NVALVES 2

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

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

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

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

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

double m_in_suc = 0.0;

int i;

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

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

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

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

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

// temperatures:

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

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

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

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

== Variants ==

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

In the paper <bibref>Larsen2007</bibref> mentioned above, the problem was stated slightly different:

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

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

\end{cases}
</math>

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

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

== References ==
<bibreferences/>


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

Annihilation of calcium oscillations

2009-08-11T12:58:04Z

PaulineBoizet: /* Source Code */

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

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

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

== Biological motivation ==

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

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

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

== Mathematical formulation ==

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

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

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

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

== Parameters ==

These fixed values are used within the model.

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

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

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

== Optimization ==

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

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

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

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
rhs[0] = p[0]
+ p[1]*xd[0]
- p[2]*xd[0]*xd[1]/(xd[0]+p[3])
- p[4]*xd[0]*xd[2]/(xd[0]+p[5]);

rhs[1] = p[6]*xd[0]
- p[7]*xd[1]/(xd[1]+p[8]);

rhs[2] = p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])
+ p[18]*p[11]*xd[1]
+ p[12]*xd[0]
- p[13]*xd[2]/(u[0]*xd[2]+p[14])
- p[15]*xd[2]/(xd[2]+p[16])
+ 0.1*xd[3];

rhs[3] = -p[9]*xd[1]*xd[2]*xd[3]/(xd[3]+p[10])
+ p[15]*xd[2]/(xd[2]+p[16])
- 0.1*xd[3];
</source>

=== Optimica ===

The optimica code of the problem:
<source lang="optimica">

package Calcium_pack

optimization Calcium_Opt (objective = cost(finalTime),
startTime = 0,
finalTime = 22)

// Parameters
parameter Real k1 = 0.09; // Parameter k1
parameter Real k2 = 2.30066; // Parameter k2
parameter Real k3 = 0.64; // Parameter k3
parameter Real K4 = 0.19; // Parameter K4
parameter Real k5 = 4.88; // Parameter k5
parameter Real K6 = 1.18; // Parameter K6
parameter Real k7 = 2.08; // Parameter k7
parameter Real k8 = 32.24; // Parameter k8
parameter Real K9 = 29.09; // Parameter K9
parameter Real k10 = 5.0; // Parameter k10
parameter Real K11 = 2.67; // Parameter K11
parameter Real k12 = 0.7; // Parameter k12
parameter Real k13 = 13.58; // Parameter k13
parameter Real k14 = 153.0; // Parameter k14
parameter Real K15 = 0.16; // Parameter K15
parameter Real k16 = 4.85; // Parameter k16
parameter Real K17 = 0.05; // Parameter K17
parameter Real p1 = 100; // Parameter p1
parameter Real x0tilde = 6.78677; // Parameter x0tilde
parameter Real x1tilde = 22.65836; // Parameter
parameter Real x2tilde = 0.384306; // Parameter x2tilde
parameter Real x3tilde = 0.28977; // Parameter x3tilde

// The states
Real x0(start=0.03966);
Real x1(start=1.09799);
Real x2(start=0.00142);
Real x3(start=1.65431);

// The control signal
input Real u(free=true);
//parameter Real umax(free=true,min=1.1,max=1.3);
//parameter Real umax=1.3;

Real cost(start=0);

equation
der(x0) = k1+k2*x0-((k3*x0*x1)/(8*x0+K4))-((k5*x0*x2)/(x0+K6));
der(x1) = k7*x0 - ((k8*x1)/(x1+K9));
der(x2)=(k10*x1*x2*x3)/(x3+K11)+k12*x1+k13*x0-(k14*x2)/(u*x2+K15)-(k16*x2)/(x2+K17)+x3/10;
der(x3)=-(k10*x1*x2*x3)/(x3+K11)+(k16*x2)/(x2+K17)-x3/10;
der(cost)=sqrt((x0-x0tilde)*(x0-x0tilde)+(x1-x1tilde)*(x1-x1tilde)+
(x2-x2tilde)*(x2-x2tilde)+(x3-x3tilde)*(x3-x3tilde))+p1*u;
constraint
x0>=0;
x1>=0;
x2>=0;
x3>=0;
// u={0,1.1,1.3};
u>=1;
u<=1.3;
end Calcium_Opt;

end Calcium_pack;

</source>

== Variants ==

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

== References ==
<bibreferences/>


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

Fuller's problem

2009-08-11T11:39:33Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==

===Solutions obtained with optimica===

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|200px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-11T11:39:06Z

PaulineBoizet: /* Reference Solutions */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==

===Solutions obtained with optimica===

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|200px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|200px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Lotka Volterra fishing problem

2009-08-11T09:57:02Z

PaulineBoizet: /* Mathematical formulation */

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

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

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

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

== Mathematical formulation ==

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


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


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

== Parameters ==

These fixed values are used within the model.

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

== Reference Solutions ==

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

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

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

== Source Code ==

Model descriptions are available in

* [[:Category:C | C code]] at [[Lotka Volterra fishing problem (C)]]
* [[:Category:AMPL | AMPL]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:optimica | optimica]] at [[Lotka Volterra fishing problem (optimica)]]

== Variants ==

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

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

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



== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Population dynamics]]

F-8 aircraft

2009-08-10T08:44:16Z

PaulineBoizet: /* Optimica */

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

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.

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.

The optimal, relaxed control function shows [[:Category:Bang bang|bang bang]] behavior. The problem is furthermore interesting as it should be reformulated equivalently.

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - 0.215 \; w + 0.28 \; x_0^2 \; w + 0.47 \; x_0 \; w^2 + 0.63 \; w^3 \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - 20.967 \; w + 6.265 \; x_0^2 \; w + 46 \; x_0 \; w^2 + 61.4 \; w^3 \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{-0.05236,0.05236\}.
\end{array}
</math>

<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<bibref>Garrard1977</bibref>.

In the control problem, both initial and terminal values of the differential states are fixed.

== Reformulation ==

The control w(t) is restricted to take values from a finite set only. Hence, the control problem can be reformulated equivalently to

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \right) \; w \\
&&& - \left( - 0.215 \; \xi + 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 + 0.63 \; \xi^3 \right) \; (1 - w) \\
& &=& -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 \\
&&& + 0.215 \; \xi - 0.28 \; x_0^2 \; \xi + 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.63 \; \xi^3 \right) \; 2 w \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \right) \; w \\
&&& - \left( - 20.967 \; \xi + 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 + 61.4 \; \xi^3 \right) \; (1 - w) \\
& &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& + 20.967 \; \xi - 6.265 \; x_0^2 \; \xi + 46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi - 61.4 \; \xi^3 \right) \; 2 w \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{0,1\},
\end{array}
</math>

with <math>\xi = 0.05236</math>. Note that there is a bijection between optimal solutions of the two problems.

== Reference solutions ==

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>.

{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Arc !! w(t) !! Lee et al.<bibref>Lee1997a</bibref> !! Kaya<bibref>Kaya2003</bibref> !! Sager<bibref>Sager2005</bibref> !! [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] !! Sager
|-
| align=center | 1 || align=center | 1 || 0.00000 || 0.10292 || 0.10235 || 0.0 || 1.13492
|-
| align=center | 2 || align=center | 0 || 2.18800 || 1.92793 || 1.92812 || 0.608750 || 0.34703
|-
| align=center | 3 || align=center | 1 || 0.16400 || 0.16687 || 0.16645 || 3.136514 || 1.60721
|-
| align=center | 4 || align=center | 0 || 2.88100 || 2.74338 || 2.73071 || 0.654550 || 0.69169
|-
| align=center | 5 || align=center | 1 || 0.33000 || 0.32992 || 0.32994 || 0.0 || 0.0
|-
| align=center | 6 || align=center | 0 || 0.47200 || 0.47116 || 0.47107 || 0.0 || 0.0
|-
| Claim: Infeasibility || align=center | - || 1.00E-10 || 7.30E-06 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|-
| Claim: Objective || align=center | - || 6.03500 || 5.74217 || 5.72864 || 4.39981 || 3.78086
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Infeasibility || align=center | - || 1.75E-03 || 1.64E-03 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Objective || align=center | - || 6.03500 || 5.74218 || 5.72864 || 4.39981 || 3.78086
|}

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>.

<gallery caption="Reference solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraftRelaxedCoarse.png| Locally optimal relaxed control on a coarse control discretization grid and corresponding differential states.
Image:f8aircraftRelaxedFine.png| Locally optimal relaxed control on a fine, adaptively chosen control discretization grid and corresponding differential states.
Image:f8aircraftInteger.png| Integer control and corresponding differential states from Sager<bibref>Sager2005</bibref>.

Image:f8aircraftSchlueter.png| Integer control and corresponding differential states from [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] solution.
Image:f8aircraftSager2009.png| Best known integer control and corresponding differential states from Sager solution.
</gallery>

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
double x1 = xd[0];
double x2 = xd[1];
double x3 = xd[2];

#ifdef CONVEXIFIED
double xi = 0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2 -x1*x1*x3 + 3.846*x1*x1*x1
+ 0.215*xi - 0.28*x1*x1*xi + 0.47*x1*xi*xi - 0.63*xi*xi*xi
- 2*u[0] * ( 0.215*xi - 0.28*x1*x1*xi - 0.63*xi*xi*xi );
rhs[1] = x3;
rhs[2] = - 4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
+ 20.967*xi - 6.265*x1*x1*xi + 46*x1*xi*xi - 61.4*xi*xi*xi
- 2*u[0] * ( 20.967*xi - 6.265*x1*x1*xi - 61.4*xi*xi*xi );
#else
double u_ = -0.05236 + u[0]*2.0*0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2
-x1*x1*x3 + 3.846*x1*x1*x1 - 0.215*u_ + 0.28*x1*x1*u_ + 0.47*x1*u_*u_ + 0.63*u_*u_*u_;
rhs[1] = x3;
rhs[2] = -4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
- 20.967*u_ + 6.265*x1*x1*u_ + 46*x1*u_*u_ + 61.4*u_*u_*u_;
#endif

</source>

=== Optimica ===
Model description in Optimica.Solver=Ipopt
Objective : 5.12799232
infeasibility : 6.2235588037251599e-10

<source lang="optimica">
package F8_aircraft_pack

optimization F8_aircraft_Opt (objective = tf,
startTime = 0,
finalTime = 1)

// The time is scaled so that the new
// time variable, call it \tau,
// is related to the original time,
// call it t as \tau=tf*t where tf
// is a parameter to be minimized
// subject to the dynamics
// \dot x = tf*f(x,w)

// Parameters
parameter Real tf(free=true,min=0.1);
parameter Real ksi=0.05236;

// The states
Real x0(start=0.4655); //angle of attack
Real x1(start=0); //pitch angle
Real x2(start=0); //pitch rate

// The control signal
input Real w(min=0,max=1);

equation

der(x0)= tf*(-0.877*x0 + x2 - 0.088*x0*x2 + 0.47*x0^2 - 0.019*x1^2 -x0^2*x2 + 3.846*x0^3
-(0.215*ksi -0.28*x0^2*ksi -0.47*x0*ksi^2 - 0.63*ksi^3)*w
-(-0.215*ksi + 0.28*x0^2 -0.47*x0*ksi^2 + 0.63*ksi^3)*(1-w));
der(x1) = tf*x2;
der(x2) = tf*(-4.208*x0 - 0.396*x2 - 0.47*x0^2 - 3.564*x0^3
+20.967*ksi - 6.265*x0^2*ksi + 46*x0*ksi^2 - 61.4*ksi^3
-(20.967*ksi - 6.265*x0^2*ksi - 61.4*ksi^3)*2*w);
constraint
x0(finalTime)=0;
x1(finalTime)=0;
x2(finalTime)=0;

end F8_aircraft_Opt;

end F8_aircraft_pack;
</source>
<gallery caption="Obtained solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraft1.png|.
Image:f8aircraft2.png|.

</gallery>

== Variants ==

* a prescribed time grid for the control function, see [[F-8 aircraft (AMPL)]],

== Miscellaneous and Further Reading ==
See <bibref>Kaya2003</bibref> and <bibref>Sager2005</bibref> for details.

== References ==
<bibreferences/>


[[Category:Bang bang]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Aeronautics]]

F-8 aircraft

2009-08-10T08:40:53Z

PaulineBoizet: /* Source Code */

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

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.

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.

The optimal, relaxed control function shows [[:Category:Bang bang|bang bang]] behavior. The problem is furthermore interesting as it should be reformulated equivalently.

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - 0.215 \; w + 0.28 \; x_0^2 \; w + 0.47 \; x_0 \; w^2 + 0.63 \; w^3 \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - 20.967 \; w + 6.265 \; x_0^2 \; w + 46 \; x_0 \; w^2 + 61.4 \; w^3 \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{-0.05236,0.05236\}.
\end{array}
</math>

<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<bibref>Garrard1977</bibref>.

In the control problem, both initial and terminal values of the differential states are fixed.

== Reformulation ==

The control w(t) is restricted to take values from a finite set only. Hence, the control problem can be reformulated equivalently to

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \right) \; w \\
&&& - \left( - 0.215 \; \xi + 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 + 0.63 \; \xi^3 \right) \; (1 - w) \\
& &=& -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 \\
&&& + 0.215 \; \xi - 0.28 \; x_0^2 \; \xi + 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.63 \; \xi^3 \right) \; 2 w \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \right) \; w \\
&&& - \left( - 20.967 \; \xi + 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 + 61.4 \; \xi^3 \right) \; (1 - w) \\
& &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& + 20.967 \; \xi - 6.265 \; x_0^2 \; \xi + 46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi - 61.4 \; \xi^3 \right) \; 2 w \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{0,1\},
\end{array}
</math>

with <math>\xi = 0.05236</math>. Note that there is a bijection between optimal solutions of the two problems.

== Reference solutions ==

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>.

{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Arc !! w(t) !! Lee et al.<bibref>Lee1997a</bibref> !! Kaya<bibref>Kaya2003</bibref> !! Sager<bibref>Sager2005</bibref> !! [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] !! Sager
|-
| align=center | 1 || align=center | 1 || 0.00000 || 0.10292 || 0.10235 || 0.0 || 1.13492
|-
| align=center | 2 || align=center | 0 || 2.18800 || 1.92793 || 1.92812 || 0.608750 || 0.34703
|-
| align=center | 3 || align=center | 1 || 0.16400 || 0.16687 || 0.16645 || 3.136514 || 1.60721
|-
| align=center | 4 || align=center | 0 || 2.88100 || 2.74338 || 2.73071 || 0.654550 || 0.69169
|-
| align=center | 5 || align=center | 1 || 0.33000 || 0.32992 || 0.32994 || 0.0 || 0.0
|-
| align=center | 6 || align=center | 0 || 0.47200 || 0.47116 || 0.47107 || 0.0 || 0.0
|-
| Claim: Infeasibility || align=center | - || 1.00E-10 || 7.30E-06 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|-
| Claim: Objective || align=center | - || 6.03500 || 5.74217 || 5.72864 || 4.39981 || 3.78086
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Infeasibility || align=center | - || 1.75E-03 || 1.64E-03 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Objective || align=center | - || 6.03500 || 5.74218 || 5.72864 || 4.39981 || 3.78086
|}

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>.

<gallery caption="Reference solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraftRelaxedCoarse.png| Locally optimal relaxed control on a coarse control discretization grid and corresponding differential states.
Image:f8aircraftRelaxedFine.png| Locally optimal relaxed control on a fine, adaptively chosen control discretization grid and corresponding differential states.
Image:f8aircraftInteger.png| Integer control and corresponding differential states from Sager<bibref>Sager2005</bibref>.

Image:f8aircraftSchlueter.png| Integer control and corresponding differential states from [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] solution.
Image:f8aircraftSager2009.png| Best known integer control and corresponding differential states from Sager solution.
</gallery>

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
double x1 = xd[0];
double x2 = xd[1];
double x3 = xd[2];

#ifdef CONVEXIFIED
double xi = 0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2 -x1*x1*x3 + 3.846*x1*x1*x1
+ 0.215*xi - 0.28*x1*x1*xi + 0.47*x1*xi*xi - 0.63*xi*xi*xi
- 2*u[0] * ( 0.215*xi - 0.28*x1*x1*xi - 0.63*xi*xi*xi );
rhs[1] = x3;
rhs[2] = - 4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
+ 20.967*xi - 6.265*x1*x1*xi + 46*x1*xi*xi - 61.4*xi*xi*xi
- 2*u[0] * ( 20.967*xi - 6.265*x1*x1*xi - 61.4*xi*xi*xi );
#else
double u_ = -0.05236 + u[0]*2.0*0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2
-x1*x1*x3 + 3.846*x1*x1*x1 - 0.215*u_ + 0.28*x1*x1*u_ + 0.47*x1*u_*u_ + 0.63*u_*u_*u_;
rhs[1] = x3;
rhs[2] = -4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
- 20.967*u_ + 6.265*x1*x1*u_ + 46*x1*u_*u_ + 61.4*u_*u_*u_;
#endif

</source>

=== Optimica ===
<source lang="optimica">
package F8_aircraft_pack

optimization F8_aircraft_Opt (objective = tf,
startTime = 0,
finalTime = 1)

// The time is scaled so that the new
// time variable, call it \tau,
// is related to the original time,
// call it t as \tau=tf*t where tf
// is a parameter to be minimized
// subject to the dynamics
// \dot x = tf*f(x,w)

// Parameters
parameter Real tf(free=true,min=0.1);
parameter Real ksi=0.05236;

// The states
Real x0(start=0.4655); //angle of attack
Real x1(start=0); //pitch angle
Real x2(start=0); //pitch rate

// The control signal
input Real w(min=0,max=1);

equation

der(x0)= tf*(-0.877*x0 + x2 - 0.088*x0*x2 + 0.47*x0^2 - 0.019*x1^2 -x0^2*x2 + 3.846*x0^3
-(0.215*ksi -0.28*x0^2*ksi -0.47*x0*ksi^2 - 0.63*ksi^3)*w
-(-0.215*ksi + 0.28*x0^2 -0.47*x0*ksi^2 + 0.63*ksi^3)*(1-w));
der(x1) = tf*x2;
der(x2) = tf*(-4.208*x0 - 0.396*x2 - 0.47*x0^2 - 3.564*x0^3
+20.967*ksi - 6.265*x0^2*ksi + 46*x0*ksi^2 - 61.4*ksi^3
-(20.967*ksi - 6.265*x0^2*ksi - 61.4*ksi^3)*2*w);
constraint
x0(finalTime)=0;
x1(finalTime)=0;
x2(finalTime)=0;

end F8_aircraft_Opt;

end F8_aircraft_pack;
</source>
<gallery caption="Obtained solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraft1.png|.
Image:f8aircraft2.png|.

</gallery>

== Variants ==

* a prescribed time grid for the control function, see [[F-8 aircraft (AMPL)]],

== Miscellaneous and Further Reading ==
See <bibref>Kaya2003</bibref> and <bibref>Sager2005</bibref> for details.

== References ==
<bibreferences/>


[[Category:Bang bang]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Aeronautics]]

F-8 aircraft

2009-08-10T08:39:40Z

PaulineBoizet: /* Optimica */

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

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.

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.

The optimal, relaxed control function shows [[:Category:Bang bang|bang bang]] behavior. The problem is furthermore interesting as it should be reformulated equivalently.

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - 0.215 \; w + 0.28 \; x_0^2 \; w + 0.47 \; x_0 \; w^2 + 0.63 \; w^3 \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - 20.967 \; w + 6.265 \; x_0^2 \; w + 46 \; x_0 \; w^2 + 61.4 \; w^3 \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{-0.05236,0.05236\}.
\end{array}
</math>

<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<bibref>Garrard1977</bibref>.

In the control problem, both initial and terminal values of the differential states are fixed.

== Reformulation ==

The control w(t) is restricted to take values from a finite set only. Hence, the control problem can be reformulated equivalently to

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \right) \; w \\
&&& - \left( - 0.215 \; \xi + 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 + 0.63 \; \xi^3 \right) \; (1 - w) \\
& &=& -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 \\
&&& + 0.215 \; \xi - 0.28 \; x_0^2 \; \xi + 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.63 \; \xi^3 \right) \; 2 w \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \right) \; w \\
&&& - \left( - 20.967 \; \xi + 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 + 61.4 \; \xi^3 \right) \; (1 - w) \\
& &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& + 20.967 \; \xi - 6.265 \; x_0^2 \; \xi + 46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi - 61.4 \; \xi^3 \right) \; 2 w \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{0,1\},
\end{array}
</math>

with <math>\xi = 0.05236</math>. Note that there is a bijection between optimal solutions of the two problems.

== Reference solutions ==

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>.

{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Arc !! w(t) !! Lee et al.<bibref>Lee1997a</bibref> !! Kaya<bibref>Kaya2003</bibref> !! Sager<bibref>Sager2005</bibref> !! [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] !! Sager
|-
| align=center | 1 || align=center | 1 || 0.00000 || 0.10292 || 0.10235 || 0.0 || 1.13492
|-
| align=center | 2 || align=center | 0 || 2.18800 || 1.92793 || 1.92812 || 0.608750 || 0.34703
|-
| align=center | 3 || align=center | 1 || 0.16400 || 0.16687 || 0.16645 || 3.136514 || 1.60721
|-
| align=center | 4 || align=center | 0 || 2.88100 || 2.74338 || 2.73071 || 0.654550 || 0.69169
|-
| align=center | 5 || align=center | 1 || 0.33000 || 0.32992 || 0.32994 || 0.0 || 0.0
|-
| align=center | 6 || align=center | 0 || 0.47200 || 0.47116 || 0.47107 || 0.0 || 0.0
|-
| Claim: Infeasibility || align=center | - || 1.00E-10 || 7.30E-06 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|-
| Claim: Objective || align=center | - || 6.03500 || 5.74217 || 5.72864 || 4.39981 || 3.78086
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Infeasibility || align=center | - || 1.75E-03 || 1.64E-03 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Objective || align=center | - || 6.03500 || 5.74218 || 5.72864 || 4.39981 || 3.78086
|}

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>.

<gallery caption="Reference solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraftRelaxedCoarse.png| Locally optimal relaxed control on a coarse control discretization grid and corresponding differential states.
Image:f8aircraftRelaxedFine.png| Locally optimal relaxed control on a fine, adaptively chosen control discretization grid and corresponding differential states.
Image:f8aircraftInteger.png| Integer control and corresponding differential states from Sager<bibref>Sager2005</bibref>.

Image:f8aircraftSchlueter.png| Integer control and corresponding differential states from [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] solution.
Image:f8aircraftSager2009.png| Best known integer control and corresponding differential states from Sager solution.
</gallery>

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
double x1 = xd[0];
double x2 = xd[1];
double x3 = xd[2];

#ifdef CONVEXIFIED
double xi = 0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2 -x1*x1*x3 + 3.846*x1*x1*x1
+ 0.215*xi - 0.28*x1*x1*xi + 0.47*x1*xi*xi - 0.63*xi*xi*xi
- 2*u[0] * ( 0.215*xi - 0.28*x1*x1*xi - 0.63*xi*xi*xi );
rhs[1] = x3;
rhs[2] = - 4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
+ 20.967*xi - 6.265*x1*x1*xi + 46*x1*xi*xi - 61.4*xi*xi*xi
- 2*u[0] * ( 20.967*xi - 6.265*x1*x1*xi - 61.4*xi*xi*xi );
#else
double u_ = -0.05236 + u[0]*2.0*0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2
-x1*x1*x3 + 3.846*x1*x1*x1 - 0.215*u_ + 0.28*x1*x1*u_ + 0.47*x1*u_*u_ + 0.63*u_*u_*u_;
rhs[1] = x3;
rhs[2] = -4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
- 20.967*u_ + 6.265*x1*x1*u_ + 46*x1*u_*u_ + 61.4*u_*u_*u_;
#endif

</source>

=== Optimica ===
package F8_aircraft_pack

optimization F8_aircraft_Opt (objective = tf,
startTime = 0,
finalTime = 1)

// The time is scaled so that the new
// time variable, call it \tau,
// is related to the original time,
// call it t as \tau=tf*t where tf
// is a parameter to be minimized
// subject to the dynamics
// \dot x = tf*f(x,w)

// Parameters
parameter Real tf(free=true,min=0.1);
parameter Real ksi=0.05236;

// The states
Real x0(start=0.4655); //angle of attack
Real x1(start=0); //pitch angle
Real x2(start=0); //pitch rate

// The control signal
input Real w(min=0,max=1);

equation

der(x0)= tf*(-0.877*x0 + x2 - 0.088*x0*x2 + 0.47*x0^2 - 0.019*x1^2 -x0^2*x2 + 3.846*x0^3
-(0.215*ksi -0.28*x0^2*ksi -0.47*x0*ksi^2 - 0.63*ksi^3)*w
-(-0.215*ksi + 0.28*x0^2 -0.47*x0*ksi^2 + 0.63*ksi^3)*(1-w));
der(x1) = tf*x2;
der(x2) = tf*(-4.208*x0 - 0.396*x2 - 0.47*x0^2 - 3.564*x0^3
+20.967*ksi - 6.265*x0^2*ksi + 46*x0*ksi^2 - 61.4*ksi^3
-(20.967*ksi - 6.265*x0^2*ksi - 61.4*ksi^3)*2*w);
constraint
x0(finalTime)=0;
x1(finalTime)=0;
x2(finalTime)=0;

end F8_aircraft_Opt;

end F8_aircraft_pack;

<gallery caption="Obtained solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraft1.png|.
Image:f8aircraft2.png|.

</gallery>

== Variants ==

* a prescribed time grid for the control function, see [[F-8 aircraft (AMPL)]],

== Miscellaneous and Further Reading ==
See <bibref>Kaya2003</bibref> and <bibref>Sager2005</bibref> for details.

== References ==
<bibreferences/>


[[Category:Bang bang]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Aeronautics]]

F-8 aircraft

2009-08-10T08:38:44Z

PaulineBoizet: /* C */

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

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.

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.

The optimal, relaxed control function shows [[:Category:Bang bang|bang bang]] behavior. The problem is furthermore interesting as it should be reformulated equivalently.

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - 0.215 \; w + 0.28 \; x_0^2 \; w + 0.47 \; x_0 \; w^2 + 0.63 \; w^3 \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - 20.967 \; w + 6.265 \; x_0^2 \; w + 46 \; x_0 \; w^2 + 61.4 \; w^3 \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{-0.05236,0.05236\}.
\end{array}
</math>

<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<bibref>Garrard1977</bibref>.

In the control problem, both initial and terminal values of the differential states are fixed.

== Reformulation ==

The control w(t) is restricted to take values from a finite set only. Hence, the control problem can be reformulated equivalently to

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w, T} & T \\[1.5ex]
\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 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \right) \; w \\
&&& - \left( - 0.215 \; \xi + 0.28 \; x_0^2 \; \xi - 0.47 \; x_0 \; \xi^2 + 0.63 \; \xi^3 \right) \; (1 - w) \\
& &=& -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 \\
&&& + 0.215 \; \xi - 0.28 \; x_0^2 \; \xi + 0.47 \; x_0 \; \xi^2 - 0.63 \; \xi^3 \\
&&& - \left( 0.215 \; \xi - 0.28 \; x_0^2 \; \xi - 0.63 \; \xi^3 \right) \; 2 w \\
& \dot{x}_1 &=& x_2 \\
& \dot{x}_2 &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \right) \; w \\
&&& - \left( - 20.967 \; \xi + 6.265 \; x_0^2 \; \xi -46 \; x_0 \; \xi^2 + 61.4 \; \xi^3 \right) \; (1 - w) \\
& &=& -4.208 \; x_0 - 0.396 \; x_2 - 0.47 \; x_0^2 - 3.564 \; x_0^3 \\
&&& + 20.967 \; \xi - 6.265 \; x_0^2 \; \xi + 46 \; x_0 \; \xi^2 - 61.4 \; \xi^3 \\
&&& - \left( 20.967 \; \xi - 6.265 \; x_0^2 \; \xi - 61.4 \; \xi^3 \right) \; 2 w \\
& x(0) &=& (0.4655,0,0)^T, \\
& x(T) &=& (0,0,0)^T, \\
& w(t) &\in& \{0,1\},
\end{array}
</math>

with <math>\xi = 0.05236</math>. Note that there is a bijection between optimal solutions of the two problems.

== Reference solutions ==

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>.

{| border="1" align="center" cellpadding="5" cellspacing="0"
|- bgcolor=#c7c7c7
! Arc !! w(t) !! Lee et al.<bibref>Lee1997a</bibref> !! Kaya<bibref>Kaya2003</bibref> !! Sager<bibref>Sager2005</bibref> !! [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] !! Sager
|-
| align=center | 1 || align=center | 1 || 0.00000 || 0.10292 || 0.10235 || 0.0 || 1.13492
|-
| align=center | 2 || align=center | 0 || 2.18800 || 1.92793 || 1.92812 || 0.608750 || 0.34703
|-
| align=center | 3 || align=center | 1 || 0.16400 || 0.16687 || 0.16645 || 3.136514 || 1.60721
|-
| align=center | 4 || align=center | 0 || 2.88100 || 2.74338 || 2.73071 || 0.654550 || 0.69169
|-
| align=center | 5 || align=center | 1 || 0.33000 || 0.32992 || 0.32994 || 0.0 || 0.0
|-
| align=center | 6 || align=center | 0 || 0.47200 || 0.47116 || 0.47107 || 0.0 || 0.0
|-
| Claim: Infeasibility || align=center | - || 1.00E-10 || 7.30E-06 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|-
| Claim: Objective || align=center | - || 6.03500 || 5.74217 || 5.72864 || 4.39981 || 3.78086
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Infeasibility || align=center | - || 1.75E-03 || 1.64E-03 || 5.90E-06 || 3.29169e-06 || 2.21723e-07
|- style="font-style:italic;background-color:#ffffcc;"
! Simulation: Objective || align=center | - || 6.03500 || 5.74218 || 5.72864 || 4.39981 || 3.78086
|}

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>.

<gallery caption="Reference solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraftRelaxedCoarse.png| Locally optimal relaxed control on a coarse control discretization grid and corresponding differential states.
Image:f8aircraftRelaxedFine.png| Locally optimal relaxed control on a fine, adaptively chosen control discretization grid and corresponding differential states.
Image:f8aircraftInteger.png| Integer control and corresponding differential states from Sager<bibref>Sager2005</bibref>.

Image:f8aircraftSchlueter.png| Integer control and corresponding differential states from [[User:MartinSchlueter | Schlueter]]/[[User:Matthias.gerdts | Gerdts]] solution.
Image:f8aircraftSager2009.png| Best known integer control and corresponding differential states from Sager solution.
</gallery>

== Source Code ==

=== C ===

The differential equations in C code:
<source lang="cpp">
double x1 = xd[0];
double x2 = xd[1];
double x3 = xd[2];

#ifdef CONVEXIFIED
double xi = 0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2 -x1*x1*x3 + 3.846*x1*x1*x1
+ 0.215*xi - 0.28*x1*x1*xi + 0.47*x1*xi*xi - 0.63*xi*xi*xi
- 2*u[0] * ( 0.215*xi - 0.28*x1*x1*xi - 0.63*xi*xi*xi );
rhs[1] = x3;
rhs[2] = - 4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
+ 20.967*xi - 6.265*x1*x1*xi + 46*x1*xi*xi - 61.4*xi*xi*xi
- 2*u[0] * ( 20.967*xi - 6.265*x1*x1*xi - 61.4*xi*xi*xi );
#else
double u_ = -0.05236 + u[0]*2.0*0.05236;

rhs[0] = -0.877*x1 + x3 - 0.088*x1*x3 + 0.47*x1*x1 - 0.019*x2*x2
-x1*x1*x3 + 3.846*x1*x1*x1 - 0.215*u_ + 0.28*x1*x1*u_ + 0.47*x1*u_*u_ + 0.63*u_*u_*u_;
rhs[1] = x3;
rhs[2] = -4.208*x1 - 0.396*x3 - 0.47*x1*x1 - 3.564*x1*x1*x1
- 20.967*u_ + 6.265*x1*x1*u_ + 46*x1*u_*u_ + 61.4*u_*u_*u_;
#endif

</source>

=== Optimica ===
ackage F8_aircraft_pack

optimization F8_aircraft_Opt (objective = tf,
startTime = 0,
finalTime = 1)

// The time is scaled so that the new
// time variable, call it \tau,
// is related to the original time,
// call it t as \tau=tf*t where tf
// is a parameter to be minimized
// subject to the dynamics
// \dot x = tf*f(x,w)

// Parameters
parameter Real tf(free=true,min=0.1);
parameter Real ksi=0.05236;

// The states
Real x0(start=0.4655); //angle of attack
Real x1(start=0); //pitch angle
Real x2(start=0); //pitch rate

// The control signal
input Real w(min=0,max=1);

equation

der(x0)= tf*(-0.877*x0 + x2 - 0.088*x0*x2 + 0.47*x0^2 - 0.019*x1^2 -x0^2*x2 + 3.846*x0^3
-(0.215*ksi -0.28*x0^2*ksi -0.47*x0*ksi^2 - 0.63*ksi^3)*w
-(-0.215*ksi + 0.28*x0^2 -0.47*x0*ksi^2 + 0.63*ksi^3)*(1-w));
der(x1) = tf*x2;
der(x2) = tf*(-4.208*x0 - 0.396*x2 - 0.47*x0^2 - 3.564*x0^3
+20.967*ksi - 6.265*x0^2*ksi + 46*x0*ksi^2 - 61.4*ksi^3
-(20.967*ksi - 6.265*x0^2*ksi - 61.4*ksi^3)*2*w);
constraint
x0(finalTime)=0;
x1(finalTime)=0;
x2(finalTime)=0;

end F8_aircraft_Opt;

end F8_aircraft_pack;

<gallery caption="Obtained solution plots" widths="240px" heights="167px" perrow="3">
Image:f8aircraft1.png|.
Image:f8aircraft2.png|.

</gallery>

== Variants ==

* a prescribed time grid for the control function, see [[F-8 aircraft (AMPL)]],

== Miscellaneous and Further Reading ==
See <bibref>Kaya2003</bibref> and <bibref>Sager2005</bibref> for details.

== References ==
<bibreferences/>


[[Category:Bang bang]]
[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Aeronautics]]

File:F8aircraft2.png

2009-08-10T08:38:08Z

PaulineBoizet:

File:F8aircraft1.png

2009-08-10T08:35:34Z

PaulineBoizet:

Fuller's problem

2009-08-07T10:13:34Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|200px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:12:56Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|500px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:11:53Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|1500px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:08:01Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|1500px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:06:56Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|500px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:06:34Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|left|270px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:04:57Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|270px|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T10:04:06Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[File:Fullerspng.png|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:46:04Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[Fichier:Fullerspng.png|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:45:23Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.
[[Fichier:Fullerspng|thumb|alt=a graph with the optimal solution of the Fuller's Problem with Optimica and Ipopt|Solution of the Fuller's Problem with Optimica and Ipopt]]
 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

File:Fullerspng.png

2009-08-07T09:38:21Z

PaulineBoizet: Solution of Fuller's Problem with Optimica and Ipopt

Solution of Fuller's Problem with Optimica and Ipopt

Fuller's problem

2009-08-07T09:32:34Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.

 

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:30:34Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:29:58Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The solution found for the relaxed Fuller's problem with optimica using the solver Ipopt (with the linear solver MA27) is obtained with 12 iterations and the objective is 1.5296058259296967e-05.

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:20:34Z

PaulineBoizet: /* C Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T09:03:38Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

* [[:Category:optimica | optimica]] at [[Fuller's Problem (optimica)]]

=== C Code ===

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T08:58:32Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

=== '''Optimica Code''' ===
----

=== C Code ===

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T08:56:58Z

PaulineBoizet: /* Optimica Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

=== Optimica Code ===
----

=== C Code ===

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Fuller's problem

2009-08-07T08:55:51Z

PaulineBoizet: /* Source Code */

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

The first control problem with an optimal [[:Category:Chattering|chattering]] solution was given by <bibref>Fuller1963</bibref>. An optimal trajectory does exist for all initial and terminal values in a vicinity of the origin. As Fuller showed, this optimal trajectory contains a bang-bang control function that switches infinitely often.

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

== Mathematical formulation ==

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

<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} & \int_{0}^{1} x_0^2 \; \mathrm{d} t \\[1.5ex]
\mbox{s.t.} & \dot{x}_0(t) & = & x_1(t), \\
& \dot{x}_1(t) & = & 1 - 2 \; w(t), \\[1.5ex]
& x(0) &=& x_S, \\
& x(t_f) &=& x_T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>

== Parameters ==

We use <math>x_S = x_T = (0.01, 0)^T</math>.

== Reference Solutions ==



== Source Code ==

=== Optimica Code ===

=== C Code ===

The differential equations in C code:
<source lang="cpp">

</source>

== Miscellaneous and further reading ==

An extensive analytical investigation of this problem and a discussion of the ubiquity of Fuller's problem can be found in <bibref>Zelikin1994</bibref>, a recent investigation of chattering controls in relay feedback systems in <bibref>Johansson2002</bibref>.

== References ==
<bibreferences/>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Bang bang]]

Template:Current News

2009-07-31T11:55:35Z

PaulineBoizet:

{| style="background:#EEEEFF;"

{{News|2009/07/31|New category [[:Category:Optimica | optimica]] introduced}}
{{News|2009/07/17|New results for [[Lotka Volterra fishing problem (AMPL)]] with Knitro}}
{{News|2009/07/09|[[F-8 aircraft (AMPL)]] page}}
{{News|2009/07/07|[[:Category:AMPL]] revised page with discretized MINLPs in AMPL format}}
{{News|2009/04/08|[[F-8 aircraft]] new solutions found}}

|}

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

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

|}
</source>

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

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

[[Category:News]]

</noinclude>

Main Page

2009-07-31T11:40:21Z

PaulineBoizet: Added optimica and C code characterization

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


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

{{:Current News}}

</td>
</tr>
</table>

<table width="100%">


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

via ''[[:Category:Model characterization|mathematical model]]'' - via ''[[:Category:Solution characterization|optimal solution]]'' - via ''[[:Category:Application|application area]]''

''[[:Category:ODE model|ODE model]]'' - ''[[:Category:Bang bang|Bang bang]]'' - ''[[:Category:Chattering|Chattering]]'' - ''[[:Category:AMPL|AMPL]]'' - ''[[:Category:C|C code]]'' - ''[[:Category:Optimica|optimica]]'' - ''[[:Category:Periodic|Periodicity]]''
</td>


<td width="50%" >

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

''[[Lotka Volterra fishing problem]]'' - [[F-8 aircraft]] - [[Annihilation of calcium oscillations | Calcium]] - [[Annihilation of calcium oscillations with PLC activation inhibition| Calcium 2]]

''[[Supermarket refrigeration system]] - ''[[Car testdrive]]'' - ''[[Fuller's problem]]''
</td>
</tr>

<tr valign="top">


<td width="50%" >

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

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


<td width="50%" >

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

''Have a look at the contributors to this benchmark library. Feel encouraged to participate!''

</td>
</tr>

</table>

[[Category:Main]]