Hybrid NMPC of supermarket display cases

ARTICLE IN PRESS Control Engineering Practice 17 (2009) 428–441

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Hybrid NMPC of supermarket display cases Daniel Sarabia a,, Flavio Capraro b, Lars F.S. Larsen c, Ce´sar de Prada a a b c

Department of Systems Engineering and Automatic Control, University of Valladolid, c/Real de Burgos s/n, 47011 Valladolid, Spain ´tica (INAUT), Universidad Nacional de San Juan, Avda. San Martı´n 1109 (oeste), 5400 San Juan, Argentina Instituto de Automa Central R&D—RA—DP, Danfoss A/S, Nordborg, Denmark

a r t i c l e in f o

a b s t r a c t

Article history: Received 3 November 2006 Accepted 4 September 2008 Available online 6 November 2008

This paper presents an NMPC of a supermarket refrigeration system. This is a hybrid process involving switched nonlinear dynamics and discrete events, on/off manipulated variables (valves and compressors), continuous controlled variables (goods temperatures) and, finally, several operation constraints. The hybrid controller is based on a parameterization of the on/off control signals in terms of time of occurrence of events instead of using directly binary values; in this way, the optimization problem can be re-formulated as an NLP problem. A rigorous model of a real supermarket refrigeration system is presented, as well as results of the hybrid controller operating on it. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Predictive control Optimization problems Nonlinear programming Variable-structure systems Discrete-event systems Variable valve timing control

1. Introduction A supermarket refrigeration system consists of a central compressor bank that, besides other process units, maintains the required flow of refrigerant to several refrigerated display cases located in the supermarket sales area. Each display case operates with an inlet on/off refrigerant valve to keep the air temperature in the display case within a range, in order to ensure a high quality of the goods. For many years, the control of supermarket refrigeration systems has been based on standard control systems. In particular, each display case is equipped with an independent hysteresis controller that regulates the air temperature in it by manipulating the inlet valve. The major drawback, however, is that the control loops are vulnerable to self-inflicted disturbances caused by the interaction between the distributed control loops. In particular, practice and simulations show that the distributed hysteresis controllers have the tendency to synchronize (Larsen, 2007; Larsen, Geyer, & Morari, 2005), meaning that the opening and closing actions of the valves coincide. Consequently, the compressor periodically has to work harder to keep up the required flow of refrigerant, which results in low efficiency, inferior control performance and a higher wear.

 Corresponding author. Tel.: +34 983 423162; fax: +34 983 423161.

E-mail addresses: [email protected] (D. Sarabia), [email protected] (F. Capraro), [email protected] (L.F.S. Larsen), [email protected] (C. de Prada). 0967-0661/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.conengprac.2008.09.003

The use of other control techniques, such as standard Model Predictive Control (MPC), is problematical because many of the control inputs are restricted to integer values, such as the opening/closing of the inlet valves and the turning on/off of the compressors. Furthermore, the system features switched nonlinear dynamics, including also slow and fast dynamics and several constraints on controlled and manipulated variables. As a result, classical control does not fit very well with the overall operation of this kind of process and other approaches must be considered. These processes that involve continuous as well as on/off variables and certain logic of operation can be formulated as hybrid systems where continuous variables are represented by real variables and the logic and on/off variables can be converted into binary variables in a set of inequalities, using for example the Mixed Logical Dynamical (MLD) framework (Bemporad & Morari, 1999). This leads, in the context of MPC, to a mixed-integer optimization problem that must be solved every sampling time, see, for instance, Prada, Sarabia, and Cristea (2005) and Larsen et al. (2005). This kind of formulation in the discrete time domain using linear models allows for a systematic controller design, but in practice, due to the high dimensionality of the problem, except in a limited number of cases, it is not possible to solve the associated MILP or MIQP optimization problem on-line. Additionally, linear models may not provide an adequate representation of the nonlinear dynamics that takes place in refrigeration systems that convert the optimization problem in an MINLP one, not suitable for real-time control. This paper presents an alternative approach of the hybrid controller where certain on/off actions are formulated in terms of

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

time instants of occurrence of events, instead of as binary values on each sampling time. This approach has been tested by the authors in other applications (Sarabia, Prada, Cristea, Mazaeda, & Colmenares, 2007), and has the advantage that all decision variables are of real type, and its number is reduced significantly, which means that the associated optimization problem is NLP and can be solved more efficiently. The controller uses the framework of nonlinear MPC, but uses a nonlinear continuous-time model, instead to a discrete one, to compute the output predictions. The control actions are calculated by means of a dynamic optimization problem solved on-line every sample time. The paper is organized as follows. In Section 2 the process is described while Section 3 presents the nonlinear model obtained from first principles. Section 4 gives an overview of the traditional control of supermarket refrigeration systems and Sections 5 and 6 describe the control objectives and the hybrid MPC, respectively. Results of some test for evaluation of the proposed approach are given in Section 7, showing the process response to setpoint changes, disturbances, etc. Finally, the paper ends with some conclusions and references.

429

In Fig. 2 a cross-section of an open display case is depicted. The refrigerant is led into the evaporator located at the bottom of the display case, where the refrigerant evaporates while absorbing heat from the wall and, ultimately, from the surrounding air circulating around the evaporator. The resulting airflow creates an air-curtain at the front of the display case. The fact that the air-curtain is colder than the goods (when the inlet valve is open) leads to a heat transfer Qgoodsair from the goods to the air-curtain. Also, there is another heat flow Qairload from the surrounding space to the air-curtain. Inside the display cases, a temperature sensor is located that measures the air temperature Tair close to the goods. This temperature measurement is used in the control loop as an indirect measure of the goods’ temperature Tgoods,

Goods

Qairload

Qgoods−air

Tgoods

2. Supermarket refrigeration systems

Air-curtain In a typical supermarket, some of the goods need to be refrigerated to ensure preservation for consumption. These goods are normally placed in open refrigerated display cases that are located in the sales area for self-service. A simplified supermarket refrigeration circuit with two display cases is shown in Fig. 1. The heart of the system is the compressor rack, which consists of a number of compressors connected in parallel. They supply the flow of refrigerant to the system by compressing the low-pressure refrigerant from the suction manifold, which receives the vapours returning from the display cases. The compressors keep a certain constant pressure in the suction manifold, thus ensuring the desired evaporation temperature. From the compressors, the refrigerant flows to the condenser and further on to the liquid manifold. The evaporators inside the display cases are fed in parallel from the liquid manifold through an expansion valve. The outlets of the evaporators lead to the suction manifold and back to the compressors, thus closing the circuit.

Tair Evaporator Twall Te Qe

Refrigerant in

Qair−wall

Refrigerant out Suction line

Inlet valve Fig. 2. Cross-section of a refrigerated display case.

hysteresis controller

Liquid manifold

Tair,1 Display case

suction line Suction manifold

hysteresis controller Tair,2 Display case

Psuc

suction line

PI-controller

Condenser unit

Compressor rack

Fig. 1. A simplified layout of a supermarket refrigeration system.

ARTICLE IN PRESS 430

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

which is not measured directly. Furthermore, an on/off inlet valve is located at the refrigerant inlet of the evaporator, which is used very often to control the temperature in the display case.

3. Nonlinear continuous-time model In this section, the main features of the nonlinear continuoustime model of the refrigeration system are presented. More details can be found in Larsen (2007) and Larsen, Izadi-Zamanabadi, Wisniewski, and Christian Sonntag (2007). The model of the supermarket refrigeration system is composed of individual models for the display cases, the suction manifold, the compressor rack and the condensing unit. The main emphasis in the modelling is laid on the suction manifold and the display cases that contain the most relevant dynamics. 3.1. The display cases The dynamic in each display case i can be described by four states, namely the goods temperature (Tgoods,i), the temperature of the evaporator wall (Twall,i), the air temperature (Tair,i) and the mass of liquefied refrigerant in the evaporator (Mref,i). The model of the display cases encompasses three parts, namely the goods, the evaporator and the air-curtain in between (see Fig. 2). By setting up the energy balance for the goods, the air-curtain and the evaporator, the following three state equations can be derived, assuming a lumped temperature model: Q goodsair;i dT goods;i ¼ dt Mgoods;i Cpgoods;i

(1)

dT wall;i Q airwall;i  Q e;i ¼ dt M wall Cpwall

(2)

dT air;i Q goodsair;i þ Q airload  Q airwall;i ¼ dt M air Cpair

(3)

where Qairload is a given external heat load on the air-curtain that is assumed equal for all display cases, and, Q goodsair;i ¼ UAgoodsair ðT goods;i  T air;i Þ

(4)

Q airwall;i ¼ UAairwall ðT air;i  T wall;i Þ

(5)

Q e;i ¼ UAwallref ;i ðM ref ;i ÞðT wall;i  T e;i Þ

(6)

where UA is the heat transfer coefficient times surface with the subscript denoting the media between which the heat is transferred, M denotes the mass and Cp the heat capacity, where the subscript also denotes the media (wall, air and refrigerant ‘‘ref’’). Moreover, Te,i is the evaporation temperature and is a refrigerant-dependent function of the evaporation pressure Pe,i. Assuming that there are no pressure drops in the suction line, the suction pressure Psuc is equal to the evaporation pressure (Pe,i ¼ Psuc); so, Te,i is the same for all evaporators (in the following Te) and depends on suction pressure. Qe,i is the cooling capacity, and, as indicated in (6), the heat transfer coefficient between the refrigerant and the evaporator wall (UAwallref,i) is a function of the mass (Mref,i) of the liquefied refrigerant in the evaporator i. It is assumed here that this relation can be described by the following linear function: UAwallref ;i ðM ref ;i Þ ¼ UAwallrefrig;max

M ref ;i M ref ;max

(7)

where Mref,i ¼ Mref,max when the evaporator is completely filled with liquid refrigerant. The accumulation of refrigerant in the

evaporator has been modelled by 8M ref ;max M ref ;i > if udisp;i ¼ 1 > tfill > dM ref ;i < Q e;i ¼  Dh if udisp;i ¼ 0 and M ref ;i 40 lg > dt > > :0 if udisp;i ¼ 0 and M ref ;i ¼ 0

(8)

The amount of liquid refrigerant in the evaporator Mref,i follows a switched nonlinear dynamic governed by (8). tfill is a parameter describing the filling time of the evaporator when the valve is opened (udisp,i ¼ 1). Furthermore, when the valve is closed (udisp,i ¼ 0) and all of the enclosed refrigerant has evaporated (Mref,i ¼ 0), then dMref,i/dt ¼ 0. So, the on/off manipulated variables are udisp,i (i ¼ 1,2,y,nd), where nd is the total number of display cases. Dhlg is the specific latent heat of the remaining liquefied refrigerant in the evaporator, which is a nonlinear function of the suction pressure. When the refrigerant leaves the evaporator it flows on to the suction manifold. Assuming that all the heat Qe,i is transferred in the two-phase region of the evaporator, the mass flow fi out of the evaporator can be computed as fi ¼

Q e;i

Dhlg

(9)

3.2. The suction manifold The dynamic in the suction manifold is modelled by one state, namely the suction pressure Psuc. Setting up the mass balance, the corresponding state equation is obtained dP suc f insuc þ f ref ;const  F comp rsuc ¼ dt V suc drsuc =dPsuc

(10)

where Vsuc is the total volume of the suction manifold, Fcomp is the volume flow out of the suction manifold determined by the P compressors, f insuc ¼ nd i¼1 f i is the total mass flow from the display cases to the suction manifold, nd is the total number of display cases and fi is the mass flow of the refrigerant out of the ith display case given by (9). fref,const is a constant mass flow into the manifold originated from other unmodelled refrigerated entities, such as cold storage rooms. rsuc is the density of the refrigerant in the suction manifold, which is a nonlinear refrigerant-dependent function of the suction pressure and the superheat in the manifold Tsh. 3.3. The compressors In most refrigeration systems, the compressor capacity is discrete-valued, as compressors can be switched only either on or off. Let nc denote the total number of compressors. The compressor bank can be modelled using a constant volumetric efficiency Zvol and the maximal displacement volume Vsl. Thus, the volume flow Fcomp,i out of the suction manifold created by the ith compressor can be determined as follows: 1 F comp;i ¼ ucomp;i C comp;i 100 Zvol V sl ;

i ¼ 1; 2; . . . ; nc (11) P where Ccomp,i is the ith compressor capacity ( nc C ¼ 100), i¼1 comp;i P and F comp ¼ nc i¼1 F comp;i is the total volume flow and, finally, ucomp,i is the on/off manipulated variable associated with each compressor i (i ¼ 1,2,y,nc). The compressors consume the major part of the energy in a refrigeration system; the power consumption is given by Powcomp ¼ F comp rsuc ðhoc  hic Þ

(12)

where Fcomp is the total volume flow mentioned before, and hic and hoc are the enthalpies of the refrigerant in and out of the

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

compressor. To get a true measure on how close the compression process is to the theoretical-most efficient, that is, to a reversible isentropic process, the isentropic efficiency his is introduced. The isentropic efficiency Zis is defined, for a process where work is added (Sonntag, Borgnakke, & Van Wylen, 1998), as

Zis ¼

W is ðh  hic Þ ¼ is W real ðhoc  hic Þ

(13)

where Wis is the required work for performing an isentropic compression process and Wreal is the real work added. Finally, the power consumption can be computed using (12) and (13) Powcomp ¼ F comp rsuc

ðhis  hic Þ

(14)

Zis

431

Psuc (bar) 2.2 2.1 2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0

2000

4000

6000

8000

10000

12000

14000

Tair,1 and Tair,2 (°C)

3.4. Refrigerant properties In the model described before, a number of refrigerantdependent properties are used, for example, Te, Dhlg, his, hie drsuc/dPsuc and rsuc, these properties can be computed by using the freeware software package ‘‘RefEqns’’ (Skovrup) available for languages C, C++, Matlab and others. The corresponding libraries in C, RefrigC.lib and RefrigC.dll, incorporating them in the simulation of the process, have been used. For a more detailed description of these thermodynamics properties, and the explanation of the vapour compression cycle in an h-log(p) diagram, see Larsen (2007).

0

2000

4000

6000

8000

10000

12000

14000

TIME (sec.) Fig. 3. Effects of synchronization on the suction pressure (Psuc) and on air temperatures (Tair,1 in black and Tair,2 in grey) using the traditional control.

4. Traditional control The control systems used in today’s supermarket refrigeration systems are, as previously mentioned, decentralized. Each of the display cases is equipped with an air temperature controller and the compressor rack is equipped with a suction pressure controller, see Fig. 1. The air temperature Tair in each display case is controlled by a hysteresis controller that opens and closes the inlet valve. This means that when Tair reaches a certain upper temperature bound T max air , the valve is opened and Tair decreases until the lower temperature bound T min air is reached and the valve is closed again: 8 0 if T air;i 4T max > air;i > < min 1 if T oT (15) udisp;i ðkÞ ¼ air;i air;i > > :u ðk  1Þ if T min oT oT max disp;i

6 5 4 3 2 1 0 -1 -2 -3 -4

air;i

air;i

air;i

When the valve is closed, the air temperature continues to decrease below the lower bound, creating an ‘‘undershoot’’. This is caused by two things: first, the remaining refrigerant contained in the evaporator evaporates and, secondly, the thermal capacity of the evaporator wall adds an extra order to the system dynamics, min   see bottom side of Fig. 3, where T max air ¼ 5 C and T air ¼ 2 C. This figure shows a simulation of a refrigeration system with two display cases and two compressors, equipped with the traditional control setup. The hysteresis controller is typically operated at a sample time of 1 s. In the supermarket, many of the display cases are alike in design and they work under uniform conditions. As a result, the inlet valves of the display cases are switched with very similar switching frequencies. Therefore, the valves have a tendency to synchronize, leading to periodic high and low flow of evaporated refrigerant into the suction manifold, thus creating large variations in the suction pressure. These effects are displayed in Fig. 3. Recall that the compressors control the suction pressure by turning off and on compressors in the compressor rack. The

suction pressure controller is normally a Proportional-Integrator controller (PI controller) with a dead band (7DB) around the reference of the suction pressure Pref suc : Z eðtÞ uPI ðtÞ ¼ K p eðtÞ þ dt (16) Ki where ( eðtÞ ¼

P suc  P ref suc

  if eðtÞ4DB

0

otherwise

(17)

e(t) is the control error and uPI is output from the PI controller. If the compressor capacity is discrete valued, then the applied control signal ucomp,i (i ¼ 1,2,y,nc) is obtained as follows: ðucomp;1 ; ucomp;2 ; . . . ; ucomp;nc Þ  8 ð0; 0; . . . ; 0Þ if uPI oC comp;1 =2 > > > > ( > > uPI XC comp;1 =2 > > > > ð1; 0; . . . ; 0Þ if > > > uPI pC comp;1 þ C comp;2 =2 < ¼ .. .. > > > . . > > > ( > > nc1 > P > > C comp;i þ C comp;nc =2 > > ð1; 1; . . . ; 1Þ if uPI X :

(18)

i¼1

The suction pressure controller typically is operated at a sample time of 60 s. The dead band ensures that the integration is stopped as long as the suction pressure is within the dead band. When the pressure exceeds the upper bound, the integration is started and one or more additional compressors are turned on, eventually, to reduce the pressure, and vice versa when the pressure falls below the lower bound. Hence, moderate changes in the suction pressure and small biases from the reference do not initiate a compressor switching. Nevertheless, pronounced synchronization effects lead

ARTICLE IN PRESS 432

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

to frequent compressor switching, causing large fluctuations in the suction pressure and a high wear on the compressors, see upper graph of Fig. 3, where P ref suc is equal to 1.4 bar and the dead band DB is fixed in 0.2. All of these problems would however be solved if a continuous control of the air temperature was applied, e.g., a cascaded control with a superheat control in the inner loop and an air temperature control in the outer loop. The inner loop would control the filling level of the evaporator (the superheat is an indirect measure of the filling level) and the outer loop would adjust the filling continuously to maintain the required air temperature. Nevertheless, there are several reasons, besides costs, that explain why hysteresis controllers are used in the display cases, despite the severe side-effects. If a mechanical superheat controller (a thermostatic expansion valve) is installed in the display cases, it is of course not easy to apply a cascaded control structure. Furthermore, if the evaporator in the display cases is operated with a low degree of filling, then the temperature distribution across the display case will be nonhomogenous, thus deteriorating the goods quality. Finally, humid air will cause a frost built up on the coldest areas of the evaporator, thus deteriorating the heat transfer and blocking the airflow through the evaporator. The hysteresis control has the positive side-effect, that in the off periods the frost will melt (Fahle´n, 1996); hence, defrosting of the evaporators can be carried out less often using hysteresis control.

the number of compressor switches and to minimize the power consumption. According to each control objective, the performance of the process is evaluated on 3 measures, gcon, gswitch and gpow. These performance measured have been initially proposed in Larsen et al. (2007).

 Constraints: It is assumed that the upper level of the suction pressure P max suc should not be surpassed elsewhere in the system, and can be treated as a soft constraint using a penalty term technique. The goods in the display case should be kept within the specified temperature bounds in order to preserve them for consumption. The air temperature is used as an indirect measure of the goods temperature. As long as the temperature max of the air is within the limits (Tair,i A[T min air ,T air ]), the temperature is considered satisfactory; however, if the bounds are crossed, the quality of the goods is reduced and this is therefore associated with a cost. This means the bounds on the

udisp,1 and udisp,2 (on/off) 1

5. Control problem The control challenge of this supermarket refrigeration system is to develop an efficient and optimal control strategy that presents better performance than the traditional control presented in Section 4 using similar instrumentation.

0 0

2000

4000

 The controlled variables are respectively the suction pressure





 



Psuc and the air temperature Tair,i in each display case (i ¼ 1,y,nd). The pressure and the temperature measurement are sampled at a sample time of 60 and 1 s, respectively. The sensors are indicated in Fig. 1. The on/off manipulated variables are the start/stop of compressors (ucomp,i 8i ¼ 1,y,nc) and the opening/closing of the inlet valves (udisp,i 8i ¼ 1,y,nd). The valves and the compressors are manipulated locally by a device with a communication system that has a band limit of 60 s. Hence, the compressors can only be manipulated every 60 s and, even if the valves can be manipulated every second, the control sequence to the valves can only be updated every 60 s. The measured disturbances are the heat load on the display cases Qairload and the additional flow of refrigerant into the suction manifold mref,const and the nonmeasured are the mass of goods Mgoods,i in each display case. Other measured state is the temperature of the wall in each display case Twall,i. The goods temperature Tgoods,i and the mass of the refrigerant Mref,i in each display case are nonmeasured states. The temperature of goods is controlled indirectly using the air temperature. Several constraints apply on suction pressure Psuc and on air temperatures Tair,i.

6000

8000

10000

12000

14000

12000

14000

ucomp,1 and ucomp,2 (on/off) 1

0 0

2000

4000

6000

8000

10000

TIME (sec.) Fig. 4. Manipulated variables in traditional control. The open/close of valves in each display case (udisp,1 in black and udisp,2 in grey), and the start/stop of the compressors (ucomp,1 and ucomp,2).

u1(t)

pulse 1

pulse 2

pulse k

1 off

on

0 5.1. Control objectives and control performance The overall control objectives are to respect the constraints on the display case temperatures and suction pressure, to minimize

T1 off1 T1 on1 T 2 off1

T 2 on1

T k off1

T k on1

Fig. 5. Parameterization of integer decision variable u1(t). Definition of the pulse and the new real decision variables Ton1k and Toff1k along time.

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

consumption gpow is an obvious performance measure:

air temperature can also be considered soft constrains. The following performance measure gcon can be interpreted as the averaged square violation of constrains per second: ! , Z tT nd 1 X 2 2 gcon ¼ Psuc ðtÞ þ T ðtÞ dt ðtT  t0 Þ (19) nd i¼1 air;i t0 where t0 and tT are the begin and end time of the given test sequence (the test sequences will be specified later under the test scenarios in Section 7); furthermore ( P suc  Pmax if Psuc 4P max suc suc (20) Psuc ðtÞ ¼ 0 otherwise and 8 max > < T air;i  T air T air;i ðtÞ ¼ T min air  T air;i > : 0

if T air;i 4T max air if T air;i oT min air

(21)

otherwise

 Switches: As mentioned earlier, one of the main objectives was to reduce the number of compressor switches. A general rule of thumb says that a compressor approximately lasts 100,000 start/stops before it wears out, i.e., by reducing the number of compressor switches the lifetime of the compressors can be prolonged. It is less critical to open/close an inlet valve than to start a compressor because on/off valves are cheap and are designed to operate in this way. It is assumed that the cost of opening a valve is 100 times less expensive than starting a compressor. The following performance measure is introduced:  Z tT   gswitch ¼ N comp ðtÞ þ Ndisp ðtÞ=100 dt ðt T  t 0 Þ (22) t0



where Ncomp(t) denotes the number of compressors started and stopped at time t and Ndisp(t) denotes the number of valves that are opened and closed at time t. The gswitch is the number of compressor switches per second, where 100 valve opening counts as one compressor switch. Power consumption: The compressors consume most of the energy in a refrigeration system. It is therefore important to reduce its power consumption. The average power

u1(t)

pulse 1

pulse 2

433

gpow ¼

Z

tT

Powcomp ðtÞdt

 ðt T  t 0 Þ

where Powcomp can be computed using (14). 6. Hybrid MPC A natural approach to many decision-making problems is MPC. An internal model of the process is used to predict its future behaviour as a function of the present and future control actions, which are selected in order to minimize some performance index or cost function J. The optimal control signals corresponding to the present time are applied to the process and the whole procedure is repeated every sampling period. These concepts require certain adaptation in order to deal with hybrid systems. Next, a description of the components of the proposed Hybrid MPC is presented: the internal model, the predictions and control horizons, the cost function and the architecture of the controller. 6.1. The internal model The internal model used in the hybrid MPC is of continuous nature and is based on the first principles one presented in Section 3. However, it is formulated in a different way. There, all manipulated variables were integers, corresponding to the opening/closing of the valves in each display case udisp,i(t) and the start/ stop of each compressor ucomp,i(t) at every sampling time along the prediction horizon. Hence, the optimization of a certain cost function J in terms of these variables implies solving a MINLP optimization problem every sampling time, which is very difficult to perform in real time. Nevertheless, it is possible to avoid the use of these integer variables associated with each sampling instant in the prediction horizon by means of a new parameterization converting the integer decision variables into real decision ones that correspond to the time instants when the binary actuator changes its state. This approach allows, besides reducing the number of decision variables, using conventional nonlinear optimization techniques (in terms only on real variables) instead

pulse 3

pulse 4

pulse 5

1 0 T1off1

T1on1

T2off1

T2on1

T 5off1

T 5on1

Nb1=1 u2(t)

pulse 1

pulse 2

pulse 3

pulse 4

1 0 T1on2

T1off 2

T 2on2

T 2off 2 Nb2=2

Start prediction

(23)

t0

T 4on2

T 4off 2 Np=4 End prediction

Fig. 6. Prediction (Np) and control horizons (Nb1, Nb2) for integer variables u1 and u2 in hybrid MPC.

ARTICLE IN PRESS 434

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

of mixed-integer nonlinear programming, decreasing the complexity of the optimization procedure and saving computation time. In the following, and for simplicity, the manipulated variable ui combines the variables associated with valves and compressors, with an index i ranging from 1 to nd and from nd+1 to nd+nc, where the first nd variables are associated with the display case valves, named before udisp,i, and the second nc variables are associated with the compressors, that is ucomp,i. The proposed parameterization defines two new real variables Ton,ik(t) and Toff,ik(t) (kAZ) for each integer manipulated variable ui (i ¼ 1,y,nd). Ton,ik and Toff,ik denote the duration of the integer variable ui when this variable is one and when it is zero, respectively, for each pulse k, see Fig. 5. The times Ton,ik and Toff,ik for all future predicted pulses along the prediction horizon are the decision variables in the optimization problem. This parameterization is based on the assumption that ‘‘the normal operation’’ of the manipulated variables, once the transient due to a disturbance or setpoint change has passed, that is at the end of the prediction horizon, will follow a periodical behaviour along the time. Notice that a ‘‘stable’’ pattern of operation cannot be reached with a constant value of the actuators, see for instance Fig. 4. Although there is no need to formulate the problem in this way, a further simplification has been considered. The hybrid MPC controller will manipulate the valves through the independent real variables Ton,ik(t) and Toff,ik(t) (i ¼ 1,y,nd) and the start/stop of the compressors (ui, i ¼ nd+1,y,nd+nc) will be manipulated with the PI controller described before. So, the internal model consists not only of the dynamical equations of the process (1)–(14) but also of the equations of the PI controller on the suction pressure (16)–(18) and the parameterization mentioned above. In principle, and from the optimization point of view, there are less degrees of freedom, but taking into account that there is also a strong relationship between the operation of the valves and the dynamics of the suction pressure, it is possible to influence its behaviour through the valves. Notice that, this approach maintains the hybrid model and computes predictions using both values ON and OFF of every position of the integer manipulated variables (on/off valves). This means that real sub-model switching takes place along the prediction horizon every time an event is activated by the valves or other changes, even if the decision variables, the commutation times, are continuous.

before, the on/off manipulated variables of the refrigeration problem will never settle in a stationary point, but will oscillate, eventually reaching a stable operational pattern. So, again, instead of specifying the interval where the predictions will be made, it is more convenient to talk in terms of Np, number of pulses of the slower changing integer manipulated variable ui, assumed to be required for obtaining a stable operational pattern. This means that, in order to compute a certain cost function J, the internal model will be integrated until the slowest train of pulses completes Np pulses. Notice that each integer manipulated

6.2. Predictions and control horizons The concept of control horizon Nu in continuous MPC corresponds to the time interval where present and future control actions are computed. After Nu sampling times, the controller maintains the last control signal computed at Nu until the end of the prediction horizon N2 in order to compute the output predictions. In the hybrid environment, and referring to binary manipulated variables, it is more convenient to formulate the control horizon Nb,i as the number of pulses of each control variable ui for which the controller computes the decision variables Tkon,i and Toff,ik (on and off periods of every pulse (k ¼ 1,2,3,y)). From the pulse number Nb,i until the end of the prediction horizon Tp ¼ N2  Ts, the values of the decision variables Ton,ik and Toff,ik will be taken equal to the ones of the pulse Nb,i; hence, the schedule of ui will be forced to follow the pattern of the Nb,i pulse, assuming that the system should reach a final stationary schedule. The prediction horizon N2 in MPC corresponds to the number of future sampling periods, of duration Ts, used to compute predictions with the internal model. In standard MPC, it is chosen longer than the process settling time. Notice that, as mentioned

Fig. 7. (a) How to calculate the term of the cost function associated with air temperatures when these variables have a periodical oscillation; (b) using a classical cost function with a penalization of the error between setpoint and air temperature; (c) the same kind of cost function but with a value lower than the previous one JcoJb. With this kind of cost function, the optimization algorithm tries to obtain this second solution Jc although in terms of periodicity of the air temperature it is worse than Jb.

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

variable can perform a different number of pulses in the same period of time Tp. Fig. 6 displays an example with two integer variables u1 and u2 and Np ¼ 4 where the prediction horizon is the time needed by u2 to implement four full pulses. Notice that u1 makes more pulses than u2 (5 pulses). In this example, the variable u1 has Nb,1 ¼ 1; 1 1 hence, only the times Ton,1 and Toff,1 that correspond to the first pulse will be computed and the same pattern will be applied in k 1 k 1 the following ones; that is to say, Ton,1 ¼ Ton,1 and Toff,1 ¼ Toff,1 8k4Nb,1 ¼ 1. For u2 Nb,2 ¼ 2 has been assigned, so the controller 1 1 2 2 computes Ton,2 , Toff,2 , Ton,2 and Toff,2 and from the second pulse on, it 2 2 follows the pattern defined by the last Ton,2 and Toff,2 computed k 2 k 2 times. Hence, Ton,2 ¼ Ton,2 and Toff,2 ¼ Toff,2 8k4Nb,2 ¼ 2.

minimum amplitude ¼

Tp

0





max þ a1 min P suc þ P suc

min

Þ2 þðT^ air;i  T refmin air;i



nd X



^ max 1þi ðT air;i

a

 T refmax Þ2 air;i (24)

(

max Psuc ðtÞ ¼

2 ðP suc  P min suc Þ

if P suc oP min suc

0

otherwise

2 ðP suc  P max suc Þ

if P suc 4P max suc

0

otherwise

(desired

Fig. 7(a) shows how to evaluate

 Constraint satisfaction: gcon is reduced with the minimization of

 (

min T ref air;i

are parameters of the hybrid controller and they are the equivalent form in the periodical environment of the setpoints, to be reached, in a classical continuous process. In general terms, the cost function (24) carries out the performance measures (19), (22) and (23) described in Section 5.1:

with

min Psuc ðtÞ ¼

and

max min and T ref before in the control performance Section 5.1. T ref air;i air;i

i¼1

dt



min T ref ). air;i

max T ref air;i

max min min T ref and T ref are different from T max air;i and T air;i described air;i air;i

The performance index J to be minimized (according to the control objectives described in Section 5) is the function Z

max T ref air;i

allowed

this term and Figs. 7(b) and (c) display the differences between this kind of penalization and the typical one, in terms of the periodicity of the controlled variables, in this case, the air temperatures, assuming that the ‘‘natural’’ behaviour of these variables is a periodical oscillation. Notice that the variables

6.3. The cost function

1 J¼ TP

values

435

(25)



the first and second terms in (24) and taking into account max ref min ref max T ref oT max 4T min air;i and T air;i air;i . If these ‘‘setpoints’’ (T air;i air;i min and T ref ) are followed, then air temperatures are always air;i min within the extended range T max air;i and T air;i used to calculate the performance index in (21). Number of switches: Decreasing the amplitude (parameter defined by the user) of air temperature oscillation in the second term of (24) reduces the variation of suction pressure with respect to its setpoint and compressors need to start/stop less times, so gswitch is reduced. Obviously, this effect is also obtained minimizing the first term in (24). Power consumption: It is not possible to reduce gpow directly, because the hybrid controller does not manipulate the

max where P min suc and P suc are the minimum and maximum values max min allowed for the suction pressure, T ref and T ref are the air;i air;i

Table 1 Parameters values for the nonlinear model

maximum and minimum values allowed for the air temperatures of each display case i, ai (i ¼ 1,y,1+nd) are appropriate weight factors and Tp is the total time of the prediction (equivalent to Np pulses, see Section 6.2). Although the hybrid controller does not manipulate the compressors, path constraints on suction pressure have been included in the first term of the cost function, where max min Psuc and Psuc in (25) penalize the values of the suction pressure below the minimum and over the maximum, respectively. On the other hand, the minimization of the second term in (24) tries to obtain a stable periodical pattern of the air temperatures in the

Parameter

Value

Parameter

Value

Display case

UAgoodsair UAairwall UAwallref,max Cp,wall Cp,air Cp,goods

300 W/K 500 W/K 900 W/K 385 J/kg K 1000 J/kg K 1000 J/kg K

Mwall Mref,max Tsh Mair

260 kg 0.6 kg 40 s 10 K 50 kg

Vsuc

5 m3

Z

0.9

Vsl Thermodynamical properties of R134a

0.08 m3/s

Zvol

0.81

Suction manifold Compressor Refrigerant

max min display cases: T^ air;i and T^ air;i are the maximum and minimum values reached by Tair,i along the prediction, and they are calculated during each prediction. The idea is to obtain waves with a certain constant amplitude fixed by the maximum and

HMPC CONTROLLER

Component

They are equal for all display cases.

Setpoints

Simulator Internal model: DAE equations

T kon T koff Constraints

Objective Function to minimize

J

Non-linear Optimizer min J

Tk

on

Tk

off

Converter real to on/off signals

ui(t)

Fig. 8. Nonlinear MPC controller implementation.

Process + PI Controller

y(t)

tfil

ARTICLE IN PRESS 436

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

compressors (they are not decision variables). But, when the compressors are on, the power has a linear dependency on suction pressure, so if suction pressure does not increase so much a small reducing of the power consumption is possible. All terms in (24) are normalized between 0 and 100 (0 if there is no error and 100 when the error reaches the maximum permitted) and notice that the influence of the free integration time Tp, favouring shorter switching times, has been compensated normalizing (24) by division by Tp.

6.4. Adaptation strategy In order to reduce the problems associated with modelling errors and unknown states and disturbances, an additional term E for each controlled variable (suction pressure and air temperatures) has been incorporated into the output equations of the internal model that are used to make the predictions: dx ¼ f ðxðtÞ; uðtÞÞ dt y ¼ gðxðtÞ; uðtÞÞ þ E

(26)

Psuc (bar) 2.2 2.1 2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9

ref

Psuc

max

Psuc

0

2000

4000

6000

8000

10000

12000

14000

Tair,1 and Tair,2 (°C) max

6

Tair

ref max

Tair

5 4 3 2 1

min

T air

ref min

Tair

0 -1 -2 -3 -4 0

2000

4000

6000

8000

10000

12000

14000

10000

12000

14000

Tgoods,1 and Tgoods,2 (°C) 6 5 4 3 2 1 0 -1 -2 -3 -4 0

2000

4000

6000

8000 TIME (sec.)

Fig. 9. Controlled variables for the first simulation: suction pressure and air temperatures (Tair,1 in black and Tair,2 in grey). The last figure shows the goods temperature (Tgoods,1 in black and Tgoods,2 in grey).

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

This term is estimated every sampling time using past and present information from the process, using a well-known strategy: as the past control is known, the model is integrated from sampling period t1 up to t and then E is obtained as the difference between the measured output ym(t) and the one predicted y^ ðtÞ by the model at time t: E ¼ ym ðtÞ  y^ ðtÞ ¼ ym ðtÞ  gðxðtÞ; uðtÞÞ

(27)

In this way, modelling errors or disturbances can be compensated in the future predictions and integral action is added to the controller. 6.5. The NMPC controller The aim of Hybrid NMPC is to minimize (24) under the dynamics (1)–(14) and PI controller on suction pressure (16)–(18) with relation to the decision variables Ton,ik and Toff,ik (k ¼ 1,y,Nb,i and i ¼ 1,y,nd) and satisfying static constraints (28) on the decision variables: max T min on;i pT on;i pT on;i max T min off ;i pT off ;i pT off ;i

8i ¼ 1; . . . ; nd

437

the internal model equations along the prediction horizon Np taking as initial conditions the current process state and evaluating the formulated cost function J at the end of the integration. The internal model was simulated using the EcosimPro simulation language (EA Int, 1999), which allows combining DAE equations with events, while performing a correct integration in spite of the model discontinuities. It belongs to the category of the so-called ‘‘modelling languages’’ that are object oriented and allows an easy re-use of the components. It generates automatically C++ code corresponding to the simulation. This code of the internal model is embedded in the controller, which has been programmed in C++. From a technical point of view, it is necessary to add a converter between the controller and the process. Hybrid MPC calculates the exact time instants Ton,ik and Toff,ik of changes for each on/off manipulated variable, and the converter transforms these times in a discrete manipulated signal ui, when the events will occur, see Fig. 8. Notice that this feature, along with the continuous nature of the internal model, allows decoupling the application of the on/off signals from the sampling time.

(28)

Pnd The optimization problem involves i¼1 2  N b;i number of decision variables. The resulting optimization problem is solved periodically, with a given sampling period using an SQP algorithm (Lawrence & Tits, 2001) (implemented in a commercial library NAG for C) using a schematic like the one of Fig. 8, which corresponds to a sequential approach where the cost function J is computed by integration of the dynamical internal model. The simulation package integrates

udisp,1 and udisp,2 (on/off) 1

7. Simulations and control results The process model described in Section 3 has been simulated and tested using also the simulation environment EcosimPro. To Mgoods,1 and Mgoods,2 (Kg) 220 180 160 140 120 100 80 60 40 20 0 0

2000

4000

6000

8000

10000

12000

14000

TIME (sec.) Fig. 11. Evolution of the nonmeasured disturbances during the second simulation, the mass of goods (Mgoods,1 in black and Mgoods,1 in grey) in each display case.

0 0

2000

4000

6000

8000

1000

12000

14000

ucomp,1 and ucomp,2 (on/off)

Mref,1 and its estimation (Kg) 0.7

1

valve open valve close

0.6 0.5 0.4 0.3 0.2 0.1

0 0

2000

4000

6000

8000

1000

12000

14000

TIME (sec.) Fig. 10. Manipulated variables for the first simulation. The open/close of valves in each display case (udisp,1 in black and udisp,2 in grey) are calculated by the hybrid controller, and the start/stop of the compressors (ucomp,1 and ucomp,2) is calculated by the PI controller. Only one compressor needs to be used (ucomp,1), the second compressor is always off (ucomp,2 ¼ 0).

0 0

50

100

150

200 250 TIME (sec.)

300

350

400

Fig. 12. Calculation of the nonmeasured state Mref,1 (mass of refrigerant) in each sampling time, using a fixed triangular pulse (dashed line). In continuous line, the exact evolution of the mass of refrigerant. The same kind of estimation has been used for Mref,2.

ARTICLE IN PRESS 438

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

illustrate the behaviour and performance of the hybrid MPC, several experiments have been made with a supermarket refrigeration system corresponding to a day/night operation. The supermarket is composed of two compressors (with the same capacity Ccomp,1 ¼ Ccomp,2 ¼ 50%) and two equally sized display cases, so that the dimension of the internal model is 58 DAEs, 15 of them are differential equations and 43 are algebraic equations and there are 1082 variables. In real supermarkets and during the nights, the display cases are covered by ‘‘night-covers’’, hence effectively decreasing the load around 40% (Qairload ¼ 3000–1800

J/s) and reducing the constant mass flow into the suction manifold (fref,const ¼ 0.2–0.0 kg/s). As a result, because of the reduced load, the maximal suction pressure is allowed to increase (P max suc ¼ 1:721:9 bar) and the setpoint and the lower limit of suction pressure are also changed in the same proportion min (P ref suc ¼ 1:421:6 bar and P suc ¼ 1:121:3 bar). The upper and lower bound for the air temperature (only to calculate the corresponding performance index) and for all display cases remain unchanged between night and day time operation, being set at min   T max air ¼ 5 C and T air ¼ 2 C. The changes between day and night

Psuc (bar) 2.2 2.1 2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0

2000

4000

6000

8000

10000

12000

14000

10000

12000

14000

10000

12000

14000

Tair,1 and Tair,2 (°C) 6 5 4 3 2 1 0 -1 -2 -3 -4 0

2000

4000

6000

8000

Tgoods,1 and Tgoods,2 (°C) 6 5 4 3 2 1 0 -1 -2 -3 -4 0

2000

4000

6000

8000 TIME (sec.)

Fig. 13. Controlled variables for the second simulation: suction pressure and air temperatures (Tair,1 in black and Tair,2 in grey). The lower figure shows the goods temperature (Tgoods,1 in black and Tgoods,2 in grey).

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

of Qairload and fref,const are well-known and do not need to be estimated. The complete test sequence covers 4 h (14,400 s.), 2 h of day and 2 h of night operation. The tuning parameters of the suction pressure PI controller were Kp ¼ 75.0, tI ¼ 50.0 s., DB ¼ 0.2 bar and a sampling time of 60 s. The initial values of the states for all simulations (including the simulation with traditional control) are Twall,1 ¼ Twall,1 ¼ 0.0 1C, Tair,1 ¼ 5.1 1C, Tair,2 ¼ 0.0 1C, Tgoods,1 ¼ Tgoods,2 ¼ 4.0 1C, Mref,1 ¼ 0.0 kg, Mref,2 ¼ 1.0 kg and Psuc ¼ 1.4 bar. Other parameter values and constants used in the nonlinear model are shown in Table 1. Regarding the hybrid MPC controller, the sampling time was set at 60 s, the prediction horizon was fixed in Np ¼ 6, which corresponds to predictions of around half an hour (1800 s), and the control horizons for every manipulated integer variable were selected as Nb,1 ¼ Nb,2 ¼ 2. So the optimization problem has 8 decision variables. Finally, the weights in cost function (24) are a1 ¼ 1.5, a2 ¼ 1 and a3 ¼ 1. Next two simulation cases are presented, the first one with all states measured (Tgoods,i and Mref,i), nor disturbances either changes in the maximum an minimum references on air max min temperatures (T ref and T ref ) and no modelling errors, air air that is, the internal model used by NMPC and the model of the process are exactly the same. However, the second simulation is used to test the robustness of the proposed approach against disturbances and lack of process information. In this way, there are nonmeasured states, Tgoods,i and Mref,i being estimated, several nonmeasured disturbances (Mgoods,i) and differences in the parameters of the internal model and the process model. All the rest values and parameters are equal in both simulations.

439

(see Fig. 4), because only the first compressor is used and the second one is always stopped (udisp,2 ¼ 0). 7.2. Nonmeasured states and nonmeasured disturbances max min and T ref are fixed in 4.4 and 2.6 1C, respectively. Now, T ref air air Moreover, the mass of goods in both display cases follows different patterns along the day (decreasing from 200 to 100 kg in display case one and from 200 to 40 kg in display case two), see Fig. 11. Notice that the internal model uses always the nominal and constant values for these nonmeasured disturbances (Mgoods,1 ¼ Mgoods,2 ¼ 200 kg) to make the predictions. On the other hand, at the beginning of each prediction, the internal max min model fix in 3.5 1C (the average between T ref and T ref ), the air air initial value for the temperature of goods Tgoods,1 and Tgoods,2 (nonmeasured states). For the other nonmeasured states Mref,i have been calculated assuming a triangular fix pulse, see Fig. 12, with a positive constant slope when the valve is opened and a negative one, when the valve is closed. Finally, a certain parameter of the model has been changed in the internal model with regard to the process model to take into account modelling errors. The value of heat transfer coefficient UAgoodsair in Eq. (4) from goods to air in both display cases was fixed in 366 W/K in the internal model and in 300 W/K in the model of the process, which corresponds to a difference of 20%. Once again, the controlled variables (see Fig. 13) present a good behaviour in presence of the disturbances on the mass of goods. In this second simulation sometimes the air temperatures decrease min below their minimum reference T ref ; this is due to the fact that air the mass refrigerant Mref,i is estimated each sampling time, and the controller does not use its exact value, the value of this state

7.1. Measured states

udisp,1 and udisp,2 (on/off) The purpose of this first simulation, where all states are measured and there is no disturbance, is to test the hybrid MPC during the day and night operation and undergo changes on the maximum and minimum references of the air temperatures. The evolution of the controlled variables, that is, suction pressure and air temperatures and their limits, is shown in Fig. 9. It is possible to observe the good behaviour and small range of these variables, being always within the imposed bounds. On the other hand and during the day, both air temperatures follow the max min exact periodical oscillation fixed by T ref and T ref in the cost air air function (24). However, during the night, small violations of these limits take place. This is due to the great amplitude in the oscillation of the suction pressure imposed by the dynamical behaviour of the process during the night. This makes difficult the task of optimization because the actions to start/stop the compressors are not decision variables in the optimization procedure. The temperature of goods (the nonmeasured states), which are actually the real variables to be controlled, is around max min the average between T ref and T ref , see the lower graph air air in Fig. 9. It is possible to compare these results with the traditional control ones in the same conditions (see Fig. 3), and view the drastic decrease of amplitude of the oscillations of all variables and what is more important, the synchronization effect on air temperature never occurs (see Fig. 9). The upper graph of Fig. 10 shows the manipulated on/off variables: opening/closing the valves and the perfect periodicity in all signals. This is a direct effect produced by the policy used to force repeated pulses at the end of the predictions (see Section 6.2). On the other hand, the bottom graph of Fig. 10 shows the start/stop of the compressors. In this approach, it is necessary to make less starts/stops of the compressors than the traditional case

1

0 0

2000

4000

6000

8000

1000

12000

14000

12000

14000

ucomp,1 and ucomp,2 (on/off) 1

0 0

2000

4000

6000

8000

1000

TIME (sec.) Fig. 14. Manipulated variables for the second simulation. The open/close movements of valves in each display case (udisp,1 in black and udisp,2 in grey) are calculated by the hybrid controller, and the start/stop of the compressors (ucomp,1 and ucomp,2) are calculated by the PI controller. Only one compressor has to be used (ucomp,1); the second compressor is always off (ucomp,2 ¼ 0).

ARTICLE IN PRESS 440

D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

(ucomp,i and udisp,i) are shown in Fig. 14 with a similar behaviour to the previous simulation. It is interesting to compare these results with the traditional control in the same conditions (see Fig. 3), and view the drastic decrease of amplitude of the oscillations of all variables and, of course, the fact that synchronization effect on air temperature never occurs (see Fig. 14). Finally Table 2 shows the corresponding performance measured of the process (proposed before in Section 5.1) during the day and night and for traditional control and hybrid MPC. All indexes obtained with the hybrid NMPC have been reduced with respect to the traditional control, the decreasing of gcon (satisfaction of the constraints) about 99% being more significant both during the day and during the night. The switching times gswitch have been decreased during the day in 89%, and the saved power consumption, gpow, in 11% during the night. In the rest of the situations, the indexes have the same valour than the traditional control. The advantages of the proposed controller are clear in relation to traditional control in terms of avoiding synchronization, decreasing variance of the air temperatures, better constraint satisfaction and decreased compressor switching. On the other hand, the use of a continuous nonlinear internal model, based on

and the modelling error added in the internal model affect strongly the predictions. The temperature of goods (the nonmeasured states), which are actually the real variables to be max controlled, is around 3.5 1C (the average between T ref and air ref min T air ), see the bottom side of Fig. 14. Manipulated variables

Table 2 Comparison of the performance measured when a traditional control is used and when the hybrid MPC is used in two different scenarios

gcon (–)

gswitch (–)

gpow (J/s)

Traditional control Day Night

1.6135 2.0375

0.0381 0.0212

13,008.5096 1003.8621

Hybrid MPC control, Simulation 1 Day Night

0.0031 0.0011

0.0042 0.0219

12,985.6424 887.7792

Hybrid MPC control, Simulation 2 Day Night

0.0031 0.0030

0.0042 0.0219

12,989.3648 890.7601

MPC Suction Pressure [bar]

4.6 4.4 4.2 4 3.8 0

20

40

60

80

100

120

140

160

180

time [min] MPC Air temperature [°C]

5 Tair1 Tair2

0

-5 0

20

40

60

80

100

120

140

160

180

120

140

160

180

time [min]

Compressor Capacity [%]

MPC 100 80 60 40 20 0 0

20

40

60

80

100

time [min] Fig. 15. Controlled variables (suction pressure and air temperatures) for MPC with MLD internal model (Larsen’s previous results) and manipulated variables (start/stop of the compressors) in the same simulation.

ARTICLE IN PRESS D. Sarabia et al. / Control Engineering Practice 17 (2009) 428–441

first principles with variable structure, allows the proposed hybrid controller to operate in a wider range and provide better performance if compared with a linear MPC. In order to compare the results obtained with the hybrid controller and those obtained with MLD plus MILP, it is important to remark that the published results with the last one (Larsen et al., 2005; Larsen et al., 2007) do not incorporate disturbances or changes day/night and all states are measured and known; even so, the oscillations in the air temperature are reduced by 64%, the oscillations in suction pressure by 75% and the number of compressor switches from 8 to 0, but in both approaches it is only necessary to use only one compressor. Compare the first 120 min in Fig. 15 corresponding to Larsen’s previous results with the day operation in Figs. 13 and 14. However, from the computational point of view, the approach presented in this paper is much more time consuming, around 50 s to find the minimum instead of 10 s of MPC with MLD internal model and MILP optimization.

8. Conclusions The control approach followed in this paper (NMPC with continuous first principles internal model, re-parameterization of integer manipulated variables and forcing the integer variables to follow a periodical pattern at the end of the control horizon) presents good results on a realistic simulation of a supermarket refrigeration system consisting of two display cases and two compressors, opening the door to economic optimization by changing the operating point while respecting constraints. In this approach, to add more devices (display cases or compressors) implies an increase in the decision variables but not in an exponential form (like in MPC with MLD models and integer decision variables, see Larsen et al., 2005), making possible the implementation of this controller to larger refrigeration systems. Another advantage is the decoupling of the manipulated variables from the sampling time; this leads to the possibility of using larger sampling times, but the control signals are always applied in it corresponding to exact instant calculated by the hybrid controller. Notice that, even though only half of the manipulated variables have been included in the scheme of the hybrid

441

controller (the open/close of the valves), the overall behaviour of the process has been advanced. A field of further research is to improve the SQP method used to find the minimum, trying to exploit the structure of the problem and to incorporate in an efficient way the manipulation of the compressors like the decision variables in the optimization procedure.

Acknowledgement This work was supported in part by the Spanish Ministry of Education and Science (former MYCT) through project DPI20030013, as well as the European Commission through VI FP NoE HYCON, Contract no. 511368. References Bemporad, A., & Morari, M. (1999). Control of systems integrating logic, dynamic, and constraints. Automatica, 35, 407–427. EA Int. (1999). EcosimPro User Manual. Available /http://www.ecosimpro.comS. Fahle´n, P. (1996). Frosting and defrosting of air-coils. Ph.D. thesis, Go¨teborg: Department of Building Services Engineering, Chalmers University of Technology. Larsen, L. F. S. (2007). Model based control of refrigeration systems. Ph.D. thesis, Aalborg, Denmark: Department of Control Engineering, Aalborg University. Larsen, L. F. S., Geyer, T., & Morari, M. (2005). Hybrid MPC in supermarket refrigeration systems. In: 16th IFAC World Congress, Prague, Czech Republic. Larsen, L. F. S., Izadi-Zamanabadi, R., Wisniewski, R., & Christian Sonntag, C. (2007). Supermarket refrigeration systems—A benchmark for the optimal control of hybrid systems control. HYCON WP2 Benchmark in /http://astwww.bci.tu-dortmund. de/hycon4b/S. Lawrence, C. T., & Tits, A. L. (2001). A computationally efficient feasible sequential quadratic programming algorithm. SIAM Journal on Optimization, 11(4), 1092–1118. Prada, C., Sarabia, D., & Cristea, S. (2005). Modelling and control of four tanks Mixed Logic Dynamical (MLD) system. Dynamics of Continuous, Discrete, and Impulsive Systems B: Applications and Algorithms, 12(2), 243–264. Sarabia, D., Prada, C., Cristea, S., Mazaeda, R., & Colmenares, W. (2007). Hybrid NMPC control of a sugar house. In: Rolf Findeisen, Frank Allgo¨wer, & Lorenz Biegler (Eds.), Assessment and future directions of nonlinear model predictive control, series: Lecture Notes in Control and Information Sciences (Vol. 358, pp. 495–502). Berlin: Springer. Skovrup, M. J. Thermodynamic and thermophysical properties of refrigerants. Software package in Borland Delphi, Version 3.00. Sonntag, R. E., Borgnakke, C., & Van Wylen, G. J. (1998). Fundamentals of thermodynamics (5th ed.). New York: Wiley.