Copyright 6th IFAC Symposium on Cost Oriented Automation Berlin, Germany, 2001
ROBUST PREDICTIVE CONTROL FOR A GREENHOUSE USING INPUT/OUTPUT LINEARIZA TION AND LINEAR MATRIX INEQUALITIES
Sandra Piiion, MigueJ Peiia, Eduardo F. Ca macho and Benjamln Kuchen
instituto de Automatiea. Fae. ingenieria. Univ. Nae. de San Juan. Av. Libertador 1109 oeste. (4500) Argentina. email:
[email protected]
Abstract: The present wor.k p:oposes a scheme for temperature control of a greenhouse. An appr~ach based .on a com~l~atlOn of two different control schemes: Input/Output Linearization an L.mear Matnx Inequahtles based pr~dicti~e :ontroller design is proposed. The algorithm consl~ts . of three steps, first a feedback hneanzatlOn, secondly the derivation of an uncertainty descnptlOn and finally the LMI-based optimization. Several simulations shows the effect of the combination of its techniques in the performance of the inside temperature for a greenhouse. Keywords: non.li.near systems, model predictive control, feedback linearization, constrained control Linear MatrIx Inequahttes. '
I . INTRODUCTION
The use the feedback linearization (FL) in Model Predictive Control (MPC) is motivated by the fact that the non linear systems description is restated as a linear one. The computational burden is reduced considerably, in spite of the linear input constraints becoming non linear (Nevistic and Morari, 1995). The main disadvantage of input/output linearization is the poor robustness that may occur in the case of model mismatch. A robust controller design method using Linear Matrix Inequalities (MPC+LMI) was described by Kothare et af. (1996). An explicit model uncertainty description is used and the predictive control problem is solved using convex optimization methods based on Linear Matirx Inequalities (LMI's). A disadvantage of the method is that the nonlinearly is handled as uncertainty on a linear model, which may lead to very conservative solutions. In this paper we propose a design method for robust non linear predictive control applying the LMI-based techniques on a input/output linearizing greenhouse system and that way combine the advantages of both methods. (As was proposed in Kothare et al., 1995). Other reason for introducing this combination approach is to increase the computing efficiency, by linearizing the plant and by reformulating the MPC problem on the new linearized coordinates. This is not
easy to solve, unless the non-linear constraints and the objective function are convex (Nevistic and Morari, 1995). Thus, the originally non-linear model used for the MPC prediction become linear, which leads to an easier implementation of the MPC algorithm and a significant reduction of the computations involved in the non-linear optimization problem. To be able to solve the robust Predictive Control Problem (PCP) with LMI-based algorithms, an uncertainty descriptions for the input/output linearized system and the feedback law have to be formulated. The algorithm suggested in this work (MPC+FL+LMI), follows the next three steps: I. Apply input/output linearization, using a feedback law that is computed based on the nominal model. 2. Find an uncertainty description for the linearized system and the nonlinear feedback law. 3. Solve the PCP at each time instant using LMI techniques. The paper is organized as follows. In section 2 the non linear greenhouse system is described together 2 as with the assumptions on signal-sets, uncertainty and boundedness. The input/output linearization of the systems is described in section 3. Section 4 gives a specific realization for the uncertainty description of the feedback linearized system in the form as described Kothare et af. (1996). The solutions of the MPC problem using LMI's is described in section 5. Finally, in section 6 shows the performance of the MPC+FL+LMI to greenhouse system.
lZ
2.
NONLINEAR MODEL FOR THE GREENHOUSE
In the paper it will be used a greenhouse model with a single layer cover. Takakura (1993) formerly discussed the greenhouse nonlinear model selected. The model takes into account the absorbed solar radiation, radiation exchange between the sky and the cover and between the cover and the interior space of the greenhouse. Heat exchanges due to convection are regarded as well: between the inside air and the cover; between the cover and the outside air, and finally, the latent heat liberated by the water vapor condensed. The energy balance of the inside air includes the convective exchange with the cover, with the seedlings and with the floor, as well as the heat exchange with the outside air. For simplicity, the air mass inside the greenhouse is considered with a homogenous temperature distribution. The system has been established as a SISO system discretized by Euler method each 3.6 sec long (tlT= 0.001 hours) and several real disturbances are considered. The state vector x(k)=(te t; tf tp ] is former by: greenhouse cover temperature te. temperature of the inside air t;, floor 1;- and plant temperature tp. The complete dynamics of the temperature inside the greenhouse get into account some perturbations such as: direct solar radiation rad, diffuse solar radiation rads> external temperature to, and the temperature under a 10 cm-depth soil layer tbl' These five variables conform the perturbations vector d(x) of (I). The temperatures te> t;, tf are measured variables. Due to the seedling temperature can not be directly measured, it is estimated by a Extended Kalman Filter. The details of its implementation are showed in Pinon et al. (1998b). Tables I and 2 showing the parameters and symbols used in the model. The discrete greenhouse model can be expressed as x(k + I) = lo(x(k» + go(k)u(k) + Po(k)d(k) (I) where XEfIt is the state vector, u E \R is the input, dE!It is the disturbances vector and y E \R is the output. /0: !It ~!ltxl, go: !It ~!ltxl, Po: !It ~!ltX4 are smooth vector fields, and ho:!It~!R is a smooth function. So x(k+I)=k(k+l) t;(k+l) tf(k+l) tp (k+l)f
and fo(x(k» = Ui(x(k»
fi(x(k)) h(x(k»
h(x(k))Y
r
alf
radiation [0.05 ND] Stefan-Boltzman cons!. [5.67xlO' s Wlm 2/C] Transmissivity of cover for long wave radiation [ND] Lower boundary of soil temperature. [20°C] Absorptivity of the inner cover for solar radiation [0.1 ND] Absorptivity of the floor for solar radiation
alp
Absorptivity of the plant for solar radiation
S;g 11,
IbI ale
[0.8ND] [O.SND]
10 Wo
Outside air temperature [0C] (Measured) Humidity ratio in the greenhouse [kglkg Dry Air] (Measured) lrun Transmissivity of the covering material [ND] Zo- ZI Depth of soil layer [0.86 m] Table 2: Parameters used in the model -h,-ho Cl = s;g.epsc-epsf C5 = s;g-epsc-epsu C? = a/c +a/c( l-alf}lds C9 = -2h;-ca.qh-2. h;.ar1o .la; CII = Cu·qh CIl = -2k/(zo+zl )-h; CI5 = d..ce Cl
=
el7 = afftran. at101 CI9 = -hlg.km.gs;
is described by II (x(k» =
Table i: List of svrnbols for model Symbol Definition abso Absorptivity of the covering material [ND] ah Average height of a greenhouse [2 m] a"o Area ratio of total plant leaf area to floor area. a"ol=i-a"o [ND] Ca Volum. heat capacity of air [1.164 KJlml/°C] Cc Heat capacity of cover [50 KJlmlo9 cp Heat capacity of a plant [41S0KJlm 0c] Cs. Heat capacity of floor [2000 KJlml0C] dex Thickness of covers [0.00 I m] Water vapor resistance [0. 132hlm ] gs; Dryness factor of the surface [0-1 ND] h; Convective heat transfer in the lOner cover [14.4 KJlmlhaC] (Nominal value) ho Convective heat transfer in the outer cover [25.2 KJlm1h°C] hlg Latent heat for evaporation [2501 KJIKg] 2 km Mass transfer coeff. [3.6h;lL.caP Kglm h] Le Lewis number [3.6 h;lkm.calp] ks Soil thermal conductivity [5.5 KJlhm°C] la; Leaf area index [0-1 ND] qh Air flow rate due to ventilation [m l /m 2h] rad Direct solar radiation [Wlm2] ~Measured) rads Diffused solar radiation lW/m ] (Calculated) p Density of air [1.164 kglml] r me Reflectivity of the outer cover for long wave
4
4
(Cite (k) + C2te (k) + C3t f (k) + cst; (k)) I CI5
12 (x(k» = (cste (k) + c 9 t; (k) + cst f (k) + ClOt p (k» I c 25 4
13 (x(k» = (Cl2tc \k) + cst;(k) + C\3t f(k) + c I4 t f (k) + C I9 WSI! (tf (k)))1 C30
CZI = Sig.epsp.epsc
C2l = - 2/a ;· h; C25 = caah C2? = alp-Iran C29 = -2p1y.hlg.la; Cl!= a,Jcqvp)
C2
= -(eps~ I)
ho C6=abso+(I-alf}lrun.ale Cs = h; CIO = 2h;.arlO .la; CI2 = S;g-eps{epse CI4 = -S;g.epsA I-rmJ CI6 = s;g-epsfepsa·ll, CIS = alflds.a"ol C20 = 2k/(zo+z/} Cn = 2/a ;·h; C24 = -s;g-epsp ( I-rme} C4 =
C26 = Sjg-epsp·epsa.l,v C2S
=
ClO =
alp·lds CsZo
4
4
f4 (x(k» = (CZlte (k) + czzt;(k) + CZ3t p(k) + c Z4 t p (k) + CZ9 WSlp (t p (k»)a rlo / C31
We assume thatfo y t1fto be affme in x(k), so fo (x(k»
and the term for latent heat from solid surfaces is, w,,/k)
O.62e
0 62
P(tt( k))
w (k) =
1.10' _eP(lt(k))
",
_
1.10' _eP(t,(k))
f3 is a factor empirical which depends on the non linear shape of the temperature profile. The nominal discrete time nonlinear system give by the affine state-space description can be formulated by The column vector g(x(k» can be formulated as go(x(k»=[O
t1T1c25
0
oy
(2)
and the control variable is the heat flow given off by the heater/cooler (3) u(k) = Qcu/{k) The perturbations vector is formulated as d{k) = [rad(k) rads{k) to{k) tblY
For the present case, the inside temperature has been chosen as the controllable variable of the system ho(x(k)) = I,(k) and the heat flow Qcol as the control variable and the matrix po(x) is given by
Po(x(k»
=
c6KI'
c711T
Cl5
Cl5
0
0
(e4 + c,.t~)IlT Cl5 cllllT
cI6·t~.IlT
cl.·KI' C30
C30 c27 11T
e2K· a r,o· AT
C30 C26'(~ ,arto ·IlT
C31
CJI
C31
0 0 0
Considering uncertainties in the transference coefficients hi and qh, parameters, the model uncertainty is described by x(k + I) = f(x(k» + g(x(k»u(k) + Pd (k)d(k) where f{x{k» go{k)=go{k)
= f{x{k»
and
+ t-.J{x(k»,
p{x{k» = p{x{k» + t-.p{x{k»
If{x{k)) = [lfl (x) tif2{x) lf3{x) lf4(X)Y
and
0 0 0 0
0 &11 t1 T
t-.p(x(k» =
a denotes the largest singular value.
i~
;,'~2
4
6
8
W
U
W
~
~
~
o
4
6
8
W
U
W
~
~
~
2
1012W1618~
"""(a.,s) 1-21 JJy 1!HJ
This work uses the measurements collected by the commercial Davis Weather Station (© David Instruments Corporation) over 21 days of July-August. This station is situated in San Juan, Argentina (31 °32 lat., 68°31 long.) Fig. I shows the external climate recorder. Temperature [DC), outside humidity ratio [kg/kg Dry Air] and wind speed [Km/h] are some of the collected data that have been introduced to the system in each control time. 3.
INPUT/OUTPUT LINEARIZACION OF A NOMINAL MODEL OF A GREENHOUSE
In this work, the input/output linearization of a greenhouse is fulfilled. The linearization adopted here may be refereed for details in Pinon et al. (1998). In such a work, a relative degree of the greenhouse system is r= I was obtained, a diffeomorphism [';(k)T,I](k)T] =
T[x{k)] where ';I(k)=tj{k),
I](k)=~e{k)
the disturbances were decoupled. Therefore, the transformed system can be expressed in normal form as ';1(k+I)=A';(k)+Bv{k) (4a)
0
0 0
0
0
and
I'lfi (x{k» = (&1 te (k) + &8 t;(k» / Cl5 !'if2(x(k) = (&8tc(k) + &9t;Ck) + &8IJ(k) + &loIp(k») I C25
and
where
a(po(x(k» < 00
0
CZ5 0
!:Jf4(x(k») =(&2Mk) + &2ip(k))1C31
_0_
a(go(x(k)) <00
0
0 0
4fj(x(k)) = (~ti(k)+Mrh(k))1 C30
Further we assume that Fo, Fa, go YPo are bounded for allx6X, so a{Fo (x(k» < 00 a(Fa (x{k» < 00
Fig. I: Climate of the weather station 0
C25 cI7·11T
Fo (x(k»x(x)
t-.J(x{k» = Fa (x{k»x{k)
P(l, (k))
._ e---;;:-;-;-::-
=
tl(k) tp{k)/ was defined as well and
I](k + I) = q(';{k), I](k» + k(';(k), I](k»d(k) (4b) y(k) =';1 (k)
(4c)
where A = Ao +M and B = Bo +M1. A feedback law was defined u(k) = a{x{k» +b(x{k»v(k) + y(x(k»d(k) (5)
where a(x(k» = -12 (x(k»C 25 b(x(k» = C25 r(x(k» = -kll
by means of the nominal linearized model is transformed to ~r(k+l)
= Ao~(k)+Bov(k)
rank k m;n $ k modeled as
k m= ' The uncertainty parameter
$
o;Ck)
where Ao= I, Bo=!:.T12 and Co= I. Have been introduced the smooth functions q and k. The details on how these functions depend on f. g. P and h can be seen in Henson (1997). In this case q=[rl(x)
1 3(x)
14(X)Y
k=[PI( x )
P3(x)
P4(X)Y ·
and Thus, the control law exactly linearizes the map between the transformed input v(k) and the output y(k). Consequently, a linear controller can be designed to satisfY control objectives such as setpoint tracking. Is important to note that the " zero dynamics remain nonlinear although in (Piii6n et al. 1998a) have been showed that it is stable. 4. UNCERTAINTY DESCRIPTION
= k;-k;nom k;m{d
k
where
= k imax -k imin
2
i mid
,,(k + I) = q(~(k), ,,(k» + k(~(k), ,,(k»d(k)
y(k) = Co~(k)
IS
k,
k
'
= hi and
k imax +k i min 2
inom
= qh
k2
Then, let be the system in LFR for A= I, B=O.I h, Cv> = I Cu=-144, Dllv=2 .328, and let Cq =[2.011 210f .
=[35.13
Du.
= [0.619
D
0
qw
35.13 43.92 0.619 0
0.773 0
40.22] 0] 2.108
Considering these values, a specific realization of the feedback uncertainty description to a 20 % of uncertainty in hi and 50% in qh, is given and substituting in (6) together with the zero dynamics leads to (I) and (5). 5. SOLUTION TO PREDICTIVE CONTROL PROBLEM USING LMI The MPC+LMI control problem is to find, in each time-instant k, a state feedback v(k)= Kx(k) such that the criterion 00
The starting point of the MPC+LMI (as described Kothare et al. 1996), is a specific uncertainty description. In this section the existence of this uncertainty description for the greenhouse linearized system is given by construction of a realization using Linear Fractional Representation (LFR). Then, since of the section 2, together with the input/output linearization law (5). A feedback uncertainty description for the perturbed linearized system exist, and can be given in the form
x(k + 1) = Ax(k) + Bv(k) + Bpp(k) + Bww(k) y(k) = Cyx(k) u(k) = Cux(k) + Duvv(k) + Duww(k) (6) q(k) = Cqx(k) + Dqww(k) p(k) = Aq(k) where x(k) is the state, v(k) the command input, z(k) the controlled output and w(k) the disturbances. !:. is a block diagonal operator with 11~ ; (k)112 = a(~ ; (k)) $1 ,
i = 1,2 •...• ni ' k ~ 0
and ni is the parameter number with uncertainties in the linearized system. Further
min max J(k) = Le;'o, (k + i I k)Q,e"o, (k + ilk) lItk)
d
;= 1
+vT(k+i-I)Rv(k+i-l)
with error (k) = y(k) -
Yref
q;(k),
is optimized subject to
a output constrain maxlly(k+i)II, ::; Y mo, , 6
'v'i > 0
(8)
The problem MPC+LMI can be reduced to a convex optimization problem using linear matrix inequalities (Kothare et al. 1996) and turn out to be a linear objective minimization problem. Theorem: Given a system with an uncertainty description showed in the previous section and a control problem of minimizing (7) subject to constraints (8) and (9). The state feedback F= YQ-I , where Q and Yare obtained from the solution (if it exists) to the following linear objective minimization problem, minimizing the worst-case J(k), can be found by solving:
r,
min
r
(9)
r .Q,Y.T.T,-
subject to the following constraints Q_QT=O
Pi (k) = 0; (k)
(k),
(7)
Q >O
(10)
i = I, ... , n
We will assume for this problem that exact measurement of the state of the system, that is k(k) li(k) IJ(k) IpCk)]r, is available and that also assume to uncertainties parameters is inside of a
(LMI for stability:)
l e ;'orQ (k Ik)] [ error (k Ik)
~0
(11)
(LMI for robust perfonnance:) Q yTRI /2 QQ I2 (CqQ+DqvY/ (AQ+BX/ 0 0 0 * j1
* * *
* * *
j1
0
* *
A
0
0
Q-BpAB;
0
Figure 2b shows the inside temperature error. Notice that it is constrained to lIerror 11 $ 3°C during the day and Ilerror
if
with
r
CqQ+~'y
(AQ+-IP)T CzT1
~
0
*
I-CyB/'yff" C,
~O
(13)
t/r, T= y
25
f?/
,,-
20
--;;.,."
r-
--/
~
.'
/~.---- Outside T empcr.tlurc
" ,'"
~15
(LMI for output constraint:)
* *
15°C during the night.
~o
(12)
Y!.Q
11 $
- - ' .. _ _ _ . .. /
.j
=~::"ar.:,1
10'----'--'---L.--'---'---'-----'--~--'---'
o
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
~~-'C:~:::t:1 o
0.2
D..
0.6
0.8
1
1.2
1.4
1.6
1.6
500~';__ __
rl==~
·5000L-~O.~2-0~.4~~O~.6-~O.~8-~,--,~.2-~,.~4-~,.6~-,~.8~~ Tome (Days) 1 _ 2 July 1999
Figure 2: Adjust of the controller to different parameters.
r
Proof: If P = iQ -I and using Schur-Complements after some simplification, this establishes LMI's (10-13). The derivation is analogous to the derivation of the LMI's in Kothare et al. (1996). 6.
SIMULA nONS RESULTS
Several simulation tests were applied to get a better knowledge on the perfonnance properties and efficiency of the Robust Model Predictive Control via VO Linearization (RMPCFL) or MPC+FL+LMI. All cases were implemented using Mathlab 5.3 software in a Pentium 550 MHz, 128 MByte RAM.
Qd = I Rd = 10 Qn = 0 Rn = 0 Qd= 1 Rd=O.OOI Qn=O Rn=O Qd= I Rd= 0.001 Qn = I Rn = 10
Figure 2c exhibit the energy consumption of the heater/cooler necessary to achieve the control objective. It is clear that the positive part of the control action belongs to the consumption of the heater and the negative part to the cooler. Both have been measured in W/m 2 . Keeping in mind the control objective the parameter with better perfonnance are Qd= I, Rd= 0.001, Qn = I and Rn = \0.
Based on the experimental results from the INT A (National Institute for Agricultural and Livestock Technology and Research), San Juan, Argentina, and on a previous work (Fullana et al. 1999) to tomatoes seedling, is necessary that the inside temperature remain around 25 oe. However, for saving energy during the night the inside temperature may be until
lOoe.
In this application, the error and the control are weight in different fonn during the day and the night to get a great energy save. Qd, Qn are defined factors that weigh the error during the day and the night, respectively. Whereas Rd , Rn denote the control's weigh . Figure 2 show some experiments carry out to adjust the controller parameters. In order to obtaining lower energy consumption during the night hours and smaller error during day hours (control objective), the control was weighed in different way. The Figure 2a shows the inside temperature, desired inside temperature (set-point) and the outside temperature.
o
~
M
M
M
U
W
ffl
W
Tirre (Days) 1-2July19ll
Figure 3: Influence of the uncertainties in hi when qh take the nominal value. hi uncertainties are of: 0 % (solid fine line); 10% (solid strong line); 20% (dashed line) and 30% (dash-doted line).
Figure 3 shows the influence of the uncertainty in hi (convective heat transfer coefficient at the inner cover). In Figure 3a can be appreciated the evolution of the inside temperature. The hi uncertainties has a sinusoidal form with peak amplitudes of 10, 20 and 30 % respectively (see Figure 3d). In this simulation, the qh parameter (airflow rate due to ventilation) is nominal constant value. When hi uncertainty is equal to 30 %, the inside temperature error is larger than the bounded error. This is due to a considerable influence of this parameter in the system behavior. It is should be noted that in every case the nominal response is shown. Rl=QCDI. Q1=1. Rl=1Q 01=1. SruiUmtirty
remarkable the importance of pondering in different ways the control action during the day and the night. 6. CONCLUSIONS A nonlinear model of a greenhouse subjected to constraints and real disturbances was analyzed. By decreasing the computation burden and making the performance analysis easier, the discussed hybrid control structure, MPC+FL+LMI, enables a general approach to the solutions of the nonlinear control problems. The proposed approach produces non linear state dependent constraints. It means that an optimization problem for a nonlinear greenhouse system subject to constraints in the output has been transformed to a new optimization problem for a Iinearized greenhouse. Due to its relative computational efficiency, the MPC+FL+LMI strategy seems to be attractive for a class of feedback Iinearizable systems. This study is intended to be continued in a future practical implementation, where the physical limitations of the heater/cooler will be considered as well as the necessary restrictions for an optimal performance.
~ ~~tj'l :' o
~
M
o
~
M
;1
M
CJ M
M
1
,11. M
1
)tJ!'
1
,~
I
U
U
l!rne (Days) 1-2July
W
W
ffl
U
ffl
U
REFERENCES
Figure 4: Influence over qh of a 20% uncertainty in hi' Figure 4 shows the influence of the uncertainties in both parameters (hi and qh). Figure 4a shows the response of the inside temperature. In this case the hi uncertainty has also a sinusoidal form with 20% of peak amplitude. The qh uncertainty, when the windows are opening and closing is assumed as 50%. Figure 4 b shows the inside temperature error. The solid line is the theoretical output error and the dashed line is the real output error. The real error does not fulfill the constraint strictly, because of the zeros dynamic is not consider in the controller. In Fig. 4c, the heater/cooler consumption of energy is shown. It should be pointed out that the energy consumption is lower than the shown in Fig. 3, although the uncertainty is higher. This is due to the qh parameter, in Fig. 3, has a nominal value (the windows are middle-opened so qh 3 2 = 5.4 m /m h, and in Fig. 4 the windows are used as a part of the control action. Figure 4d shows the control gain in the experimental period. Notice that during the hours without solar radiation, the controller has a small gain (Fcontrol "" 0). This is one of the reasons for achieving a lower energetic consumption. Is
Fullana, R., and C. Schugurensky, (1999). "Optimization of constrained systems by iterative dynamical programming" (In Spanish). VIII RPIC, Argentina, 23-25 Setiembre. 3: 11- JO - 17-10. Henson M., and D. Seborg, Nonlinear process control, Prentice Hall, New Jersey (1997). Kothare, M. V., V. Balakrishnan and M. Morari, (1996). "Robust Constrained Predictive Control for non linear systems". Automatica. Nevistic, V. and M. Morari, (1995). "Constrained control of Feedback-Linearized Systems,". ECC. Piii6n S., C. Soria, C. Garcia and B. Kuchen (1998a). "Optimal Control of a greenhouse by feedback linearization". IEEE 98 INDUSCON. 261-265. Sao Paulo. Piii6n S., C. Soria, and B. Kuchen (I 998b). "Nonlinear Greenhouse's state variables optimal estimation ". XVI Argentine Congress of Automatic Control, 338343. Takakura T. (1993) Climate under Cover, Digital Dynamic Simulation in Plant Bio-Engineering, Kluwer Academic Publishers, The Netherlands.