Supermarket refrigeration system

Supermarket refrigeration system
State dimension:	1
Differential states:	9
Discrete control functions:	4
Interior point inequalities:	5

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.

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 ODE model. The initial values of the differential states are not fixed but periodicity of the whole process is required.

The optimal integer control functions shows chattering behavior, making the supermarket refrigeration system problem a candidate for benchmarking of algorithms.

Mathematical formulation

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

$\min_{x,u} \frac {1}{t_f - t_0}\int_{t_0}^{t_f} \left( (u_2 + u_3)\cdot 0.5 \cdot \eta_{vol} \cdot V_{sl} \cdot f \right) dt$

$\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+u_3\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 3\}, \\ & t_f &\in& [ 650, 750 ]. \end{array}$

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

$(x_5, x_6, x_7, x_8)$ describe the corresponding states in the second display case.

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

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

$\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}$

Parameters

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

$\begin{array}{lcrr} t_0 &=& 0.00, & s\\ \dot{Q}_{airload} &=& 3000.00, & \frac{J}{s} \\ \dot{m}_{ref,const} &=& 0.20, & \frac{kg}{s} \\ M_{goods} &=& 200.00, & kg\\ C_{p,goods} &=& 1000.00, & \frac{J}{kg \cdot K} \\ UA_{goods-air} &=& 300.00, & \frac{J}{s \cdot K} \\ M_{wall} &=& 260.00, & kg\\ C_{p,wall} &=& 385.00, & \frac{J}{kg \cdot K} \\ UA_{air-wall} &=& 500.00, & \frac{J}{s \cdot K} \\ M_{air} &=& 50.00, & kg\\ C_{p,air} &=& 1000.00, & \frac{J}{kg \cdot K} \\ UA_{wall-ref,max} &=& 4000.00, & \frac{J}{s \cdot K} \\ \tau_{fill} &=& 40.00, & s\\ T_{SH} &=& 10.00, & K\\ M_{ref,max} &=& 1.00, & kg\\ V_{suc} &=& 5.00, & m^3\\ V_{sl} &=& 0.08, & \frac{m^3}{s} \\ \eta_{vol} &=& 0.81, & -\\ \end{array}$

Reference Solutions

For the relaxed problem (we only demand $u_i(t) \in [0,1]$ instead of $u_i(t) \in \{0,1\}$ 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.

Source Code

C

The differential equations in C code:

// 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]));
}

Variants

Since the compressors are parallel connected one can introduce a single control $u_2 \in \{0,1,2\}$ instead of two equivalent controls. The same holds for scenarions with $n$ 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.
• A night scenario with two different parameters was given. At night the following parameters change their value:

$\begin{array}{lcrr} \dot{Q}_{airload} &=& 1800.00 & \frac{J}{s}, \\ \dot{m}_{ref,const} &=& 0.00 & \frac{kg}{s}, \\ \end{array}$

Additionally the constraint on the suction pressure $x_0(t)$ is softened to $x_0(t) \leq 1.9$.

• No periodicity was required but the solution on a fixed time horizon 4 hours - 2 in day scenario and 2 in night scenario - with $t_f = 14400$ 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 $V_{sl} = 0.095 \frac{m^3}{s}$. Unfortunately this constant is only given for these two cases although Larsen proposed scenarios with more compressors and display cases.

References

<bibreferences/>