Min-Max Predictive Control of a Pilot Plant using a QP Approach

Min-Max Predictive Control of a Pilot Plant using a QP Approach J.K. Gruber, D.R. Ramirez, T. Alamo, C. Bordons, E.F. Camacho Abstract—The practical implementation of Min-Max MPC (MMMPC) controllers is limited by the computational burden required to compute the control law. This problem can be circumvented by using approximate solutions or upper bounds of the worst possible case of the performance index. In a previous work, the authors presented a computationally efficient MMMPC control strategy in which a close approximation of the solution of the min-max problem is computed using a quadratic programming problem. In this paper, this approach is validated through its application to a pilot plant in which the temperature of a reactor is controlled. The pilot plant is operated with a Simatic IT SCADA system. The controller has been implemented in Matlab and connected to the SCADA using OPC. Realistic values of the parameters of the MMMPC controller have been used. The behavior of the system and the controller is illustrated by means of experimental results. I. INTRODUCTION The idea behind min-max model predictive controllers (MMMPC) is not new ([4]). In these controllers, the control signal is computed for the worst case of a cost function that considers the effect of process model uncertainties and disturbances in the controller performance. The main drawback of this approach is the computational burden that takes to compute the control signal. This usually involves the solution of an NP-hard min-max problem ([10], [16]). As a result, the number of applications of these control strategies is very small, even when there is evidence that they work better than standard predictive controllers in processes with uncertain dynamics [5], [8]. Multi-parametric programming has been applied to show that the MMMPC control law is piecewise affine when a quadratic ([14]) or 1-norm based criterion ([2], [7]) is usedas cost function. Thus, explicit forms of the control law can be built. Such explicit forms can be evaluated very fast provided that the complexity of the state space partition is moderate, which is the case for many applications. However, if the process model or the controller tuning parameters change, the computation of the controller has to be redone. Acommonsolutiontothecomputationalburdenissueis to use an upper bound of the worst case cost instead of computing it explicitly. This upper bound can be computed by using linear matrix inequalities (LMI) techniques such as in [9], [11]. However, the LMI problems have a computational burden that cannot be neglected in certain applications. In [1] adifferentapproachbasedonacomputationallycheapupper bound of the worst case cost is presented. In that work, the The authors are with the Dept. de Ingenier´ıa de Sistemas y Autom´atica, Escuela Superior de Ingenieros, University ofSeville,Spain {jgruber, danirr, alamo}@cartuja.us.es, {bordons, eduardo}@esi.us.es min-max problem is replaced by a quadratic programming (QP) problem that provides a close approximation to the solution of the original min-max problem. The computational burden is much lower than that of the min-max problem and is comparable to that of a standard constrained MPC based on a quadratic cost function. Thus, it can be easily implemented in almost any platform capable to run a constrained MPC. Also, stability of the proposed approach is guaranteed. In this work, the approach presented in [1] has been validated by means of its application to a pilot plant. The pilot plant is used to simulate an exothermic chemical reaction with nonlinear dynamics. This process has been used in previous works, thus the experimental results presented can be compared with other strategies such as nonlinear and linear predictive control [6]. In the experiments, restrictions in the control action and the output have been considered. The results obtained prove the validity of the control strategy. The low computational burden of the control strategy applied to the pilot plant allows realistic values for the control and prediction horizons (i.e., the parameters on which the computational burden depends). It is also noteworthy that the computer on which the predictive control algorithm is implemented and simultaneously executes the SCADA Simatic-IT, does not have sufficient calculation power to implement a conventional min-max predictive control strategy. Therefore acomputationallyefficientstrategyasthatusedinthiswork is a good choice if the use of this type of control is desired. The paper is organized as follows: section II presents the MMMPC strategy. Section III presents the proposed implementation strategy. In section IV a detailed description of the used pilot plant is given. The strategy is illustrated by means of experimental results of the pilot plant in section V. Finally, section VI presents some conclusions. II. MIN-MAX MPC WITH BOUNDED ADDITIVE UNCERTAINTIES Consider the following state space model with bounded additive uncertainties ([3]): x(t+1)=Ax(t)+Bu(t)+D θ (t+1)(1) with x(t)∈Rdimx the state vector, u(t)∈Rdimu the input vector and θ (t)∈{ θ ∈Rdim θ :∥ θ ∥∞≤ ε }the uncertainty, that is supposed to be bounded. The system is subject to p state and input time invariant constraints Fuu(t)+Fxx(t)≤g where Fu∈Rp×dimu and Fx∈Rp×dimx.Itisassumedasemifeedback approach in which the control input is given by u(t)=−Kx(t)+v(t),(2) Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. where the feedback matrix Kis chosen to achieve some desired property such as nominal stability or LQR optimality without constraints. The MMMPC controller will compute the optimal sequence of correction control inputs v(t).The state equation of system (1) can be rewritten as x(t+1)=ACLx(t)+Bv(t)+D θ (t+1),ACL =(A−BK).(3) The proposed strategy also works without semi-feedback approach (i.e., u(t)=v(t)). All the computational advantages of the strategy remain the same. Furthermore if the process is open-loop stable (as in the case of pilot plant used in this work) the stabilizing conditions, which will be discussed later, can be used without problems. The cost function is a quadratic performance index: V(x,v, θ θ θ )= N−1 ∑ j=0 x(t+j|t)TQx(t+j|t)+ (4) N−1 ∑ j=0 u(t+j|t)TRu(t+j|t) +x(t+N|t)TPx(t+N|t) where x(t|t)=x,x(t+j|t)is the prediction of the state for t+jmade at tand u(t+j|t)=−Kx(t+j|t)+v(t+ j|t).Notethatthesevaluesdependonthefuturevalues of the uncertainty. The sequence of future values of θ (t)over a prediction horizon Nis denoted by θ θ θ = ! θ (t+1)T,···, θ (t+N)T"T,andΘ Θ Θ={ θ θ θ ∈RN·dim θ :∥ θ θ θ ∥∞≤ ε }is the set of possible uncertainty trajectories. On the other hand, v=!v(t|t)T,···,v(t+N−1|t)T"Tis the control correction sequence. Matrices Q,P∈Rdimx×dimx and R∈ Rdimu×dimu are symmetric positive definite matrices used as weighting parameters. Min-Max MPC ([4]) minimizes the cost function for the worst possible case of the predicted future evolution of the process state or output signal. This is accomplished through the solution of a min-max problem: v∗(x)= argmin vmax θ θ θ ∈ΘV(x,v, θ θ θ ) s.t.Fuu(t+j|t)+Fxx(t+j|t)≤g, j=0,...,N,∀ θ θ θ ∈Θ, x(t+N|t)∈Ω,∀ θ θ θ ∈Θ, (5) Aterminalregionconstraintx(t+N|t)∈Ω,whereΩis a polyhedron, is included to assure stability of the control law ([12]). The predictions x(t+j|t)and u(t+j|t)depend linearly on x,vand θ θ θ .Thismeansthatitispossibletofindavector d∈Rpand matrices Gx,Gvand G θ ,suchthatalltherobust linear constraints of problem (5) can be rewritten as: Gi xx+Gi vv+Gi θθ θ θ ≤di,i=1...,p,∀ θ θ θ ∈Θ, where Gi x,Gi v,Gi θ denote the i-th rows of Gx,Gvand G θ respectively and diis the i-th component of d∈Rp.Denote now ∥Gi θ ∥1the sum of the absolute values of row Gi θ . Taking into account that max θ θ θ ∈ΘGi θθ θ θ =max∥ θ θ θ ∥∞≤ ε Gi θθ θ θ = ε ∥Gi θ ∥1,therobustfulfillmentoftheconstraintsissatisfied if and only if Gi xx+Gi vv+ ε ∥Gi θ ∥1≤di,i=1,...,p. Therefore, to guarantee robust constraint satisfaction, the following set of linear constraints must be satisfied: Gxx+Gvv≤d ε . where the i-th component of d ε is equal to di− ε ∥Gi θ ∥1.Note that this is a necessary and sufficient condition. Taking into account (3),(2) and (4), the cost function can be evaluated as a quadratic function: V(x,v, θ θ θ )=vTMvvv+ θ θ θ TM θθθ θ θ +2 θ θ θ TM θ vv +2xTMT vfv+2xTMT θ f θ θ θ +xTMffx(6) where the matrices can be obtained from the system and the control parameters ([3]). Due to the convexity properties of V(x,v, θ θ θ ),problem(5)isequivalentto([3]) v∗(x)=arg min vmax θ θ θ ∈vert(Θ Θ Θ)V(x,v, θ θ θ ) s.t.Gxx+Gvv≤d ε (7) where vert(Θ Θ Θ)is the set of vertices of Θ Θ Θ. The terminal region Ωis assumed to satisfy the following conditions: •C1:Ifx∈Ωthen ACLx+D θ ∈Ω,forevery θ ∈{ θ ∈ Rdim θ :∥ θ ∥∞≤ ε }. •C2:Ifx∈Ωthen u(x)=−Kx ∈U,whereU!{u: Fuu+Fxx≤g}. Moreover, matrix Pthat characterizes the terminal cost is assumed to satisfy •C3:P−AT CLPACL >Q+KTRK. The stability of ACL guarantees the existence of a positive definite matrix Psatisfying C3. The maximum cost for a given xand vis denoted as V∗(x,v)= max θ θ θ ∈vert(Θ Θ Θ)V(x,v, θ θ θ )=V(x,v,0)(8) +max θ θ θ ∈vert(Θ) θ θ θ TH θ θ θ +2 θ θ θ Tq(x,v) where H=M θθ ,q(x,v)=M θ vv+M θ fxand V(x,v,0)= vTMvvv+2xTMT vfv+xTMffxis the part of the cost that does not depend on the uncertainty. With this definition, problem (7) can be rewritten as v∗(x)= argmin vV∗(x,v) s.t.Gxx+Gvv≤d ε ,(9) and the system is controlled by KMPC(x(t)) = −Kx(t)+ v∗(t|t),wherev∗(x(t)) = !v∗(t|t)T,···,v∗(t+N−1|t)T"T. III. A QP APPROACH TO MIN-MAX MPC In this section the main results of [1] are presented briefly. In that work, it is shown how the min-max problem (9) can be replaced by a tractable QP problem which provides a close approximation of the solution of the original problem. This can be accomplished with the following steps: 1) Obtain an initial guess of the solution of (9), denoted ˜v ∗.Asseenlater,thiscanbeachievedbysolvingaQP problem. 2) Using ˜v ∗,obtainaquadraticfunctionofvthat bounds the worst case cost. Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. 3) Compute the control law. This involves the solution of aQPproblem. All these steps will be detailed in the following. A. Computing ˜v ∗ Given Hdefined as in equation (8), denote Ti= ∑N·dim θ j=1|Hij|,whereHij denotes the (i,j)-th component of matrix H.Then,definethediagonalmatrixTas T=diag(T1,···,Tn)(10) Because of how matrix Tis defined, T−His a symmetric diagonally dominant real matrix with nonnegative diagonal entries, thus T−H≥0whichimpliesthatT≥H.Let ˜ V(x,v, θ θ θ )be: ˜ V(x,v, θ θ θ )=V(x,v,0)+ θ θ θ TT θ θ θ +2qT(x,v) θ θ θ (11) From the inequality T≥Hit is inferred that ˜ V(x,v, θ θ θ )≥ V(x,v, θ θ θ ).Themaximumof ˜ V(x,v, θ θ θ )can be computed as ˜ V∗(x,v)=max θ θ θ ∈Θ Θ Θ ˜ V(x,v, θ θ θ ) =V(x,v,0)+trace(T) ε 2+2 ε ∥q(x,v)∥1 =V(x,v,0)+∥H∥s ε 2+2 ε ∥q(x,v)∥1(12) where ∥H∥sdenotes the sum of the absolute values of the elements of H.Thenaninitialguessofthesolutionof(9) can be obtained as ˜v ∗(x)=argmin ˜v ˜ V∗(x,˜v ) s.t.Gxx+Gv˜v ≤d ε ,(13) This problem can be casted as a QP problem by making use of slack variables to deal with the 1-norm term in ˜ V∗(x,v). B. Obtaining an upper bound of the worst case cost The upper-bound of the maximum will be obtained in two steps. In the first one we compute a set of parameters from ˜v ∗that allows us later, in the second step, to compute the bound as a quadratic function of v. 1) Computing the parameter vector α (v):Note that: V∗(x,v)= max θ θ θ ∈vert(Θ Θ Θ)# θ θ θ 1$T#Hq(x,v) qT(x,v)V(x,v,0)$# θ θ θ 1$ =max ∥z∥∞≤1 zTM(v)z(14) with z=% θ θ θ T ε 1 &T ,M(v)= # ε 2H ε q(x,v) ε qT(x,v)V(x,v,0)$∈Rn×n, where n=N·dim θ +1. The following procedure provides an upper bound of the worst case cost for a given v.Itcomputes α (v)= [ α 1(v),..., α n−1(v)]Tand a diagonal matrix Γ(v)≥M(v) such that its trace is an upper bound of the worst case cost for v(see property 1 of [1]). Procedure 1: Computation of α (v)= [ α 1(v),..., α n−1(v)]Tand Γ(v). 1) Let S(0=M(v)∈Rn×n. 2) For k=1ton−1 3) Let M(k−1 sub =[S(k−1 ij ]for i,j=k···n. 4) Obtain the partition M(k−1 sub =#ab T bM r$,where a∈R,b∈Rn−kand Mr∈R(n−k)×(n−k). 5) Make α k(v)='∥b∥1. 6) If α k(v)=0thenS(k=S(k−1,elseS(k=S(k−1+ (0T k−1,1 α k(v)−bT α k(v))T(0T k−1,1 α k(v)−bT α k(v)). 7) end for 8) Make Γ(v)=S(n−1. Note that in the previous procedure, 0m,ndenotes a (m×n) matrix of zeros. Property 1 of [1] shows that the trace of Γ(v)constitutes an improved upper bound of V∗(x,v).That is, V∗(x,v)≤trace(Γ(v)) ≤˜ V∗(x,v). 2) Obtaining the bound as a quadratic function on v: The diagonalization process shown in procedure 1 can be used to obtain a matrix denoted by ˆ Γ(v),whichallowsone to obtain a bound of the maximum that can be computed as aquadraticfunctionofv.Thisisachievedbymeansofthe following procedure: Procedure 2: Obtaining the matrix ˆ Γ(v). 1) Obtain ˜v ∗from the QP problem defined in eq. (13). 2) Compute α (˜v ∗)by procedure 1. 3) Let ˆ S(0(v)=M(v)∈Rn×n. 4) For k=1ton−1 5) Let ˆ Msub(v)=[ˆ S(k−1 ij (v)] for i,j=k···n. 6) Obtain the partition ˆ Msub(v)=#a(v)bT(v) b(v)Mr(v)$, where a(v)∈R. 7) If α k(˜v ∗)=0then ˆ S(k(v)= ˆ S(k−1(v),else ˆ S(k(v)= ˆ S(k−1(v)+ (0T k−1,1 α k(˜v ∗)−b(v)T α k(˜v ∗))T(0T k−1,1 α k(˜v ∗)−b(v)T α k(˜v ∗)). 8) end for 9) Make ˆ Γ(v)= ˆ S(n−1(v). Denote that ˆ V∗(x,v)=trace (ˆ Γ(v)).Theorem1of[1] shows that ˆ V∗(x,v)is a quadratic function on vand also an upper bound of the original worst case cost V∗(x,v). C. Computing the control law The value of the control signal is obtained by solving the following QP optimization problem ˆv ∗(x)=argmin ˆv ˆ V∗(x,ˆv ) s.t.Gxx+Gvˆv ≤d ε ,(15) and the system is controlled by ˆ KMPC(x(t)) = −Kx(t)+ ˆv∗(t|t),whereˆv∗(t|t)is the first element of ˆ v∗(x). The computational burden of the proposed strategy is much lower than that of the exact MMMPC. This computational burden is mostly due to the QP problems that must be solved to obtain the initial guess and the proposed solution itself. Note that in both cases, the complexity of each problem is the same as that of a standard constrained MPC using a quadratic cost function. In contrast to evaluate the maximum cost V∗(x,v)it is necessary to evaluate the function for all the 2N∗dim θ vertices of Θ.Notethatthisisa well known NP-hard problem. Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. Fig. 1. Pilot plant used to apply the MMMPC. Heat exchanger Water tank v8 Fj Fj Tj,out Tj,in Ff Ff TT2 T Fig. 2. Diagram of the pilot plant with its four main elements:reactor, heat exchanger, cooling jacket and valve. IV. PROCESS DESCRIPTION Arealprocessrepresentedbyapilotplanthasbeenchosen for the application of the proposed algorithm. The process has been studied previously by several authors [6], [17]. A. Laboratory process The pilot plant (see Fig. 1) is used to simulate exothermic chemical reactions based on temperature changes. It has been used as a benchmark for control purposes by several researchers [13], [6]. The main elements of the pilot plant are the reactor, the heat exchanger, the cooling jacket and the valve to manipulate the flow rate through the cooling jacket (see Fig. 2). Acoolingjacketisusedtoreducethereactortemperature. The heat dissipation can be regulated by the valve v8which manipulates the flow rate Fjthrough the cooling jacket. The cooling fluid, water, enters the cooling jacket with a constant temperature. The reactive is supplied to the reactor by the feed Ff,in to keep the chemical reaction active. Before entering the reactor, the feed passes through a heat exchanger in order to adopt the temperature of the reactor content. The outflow Ff,out is used to keep the volume of the reactor content constant. To simulate exothermic reactions, the reactor possesses an electrical resistance in order to supply caloric energy. The energy to be supplied by the 14.4kW electrical resistance is calculated by means of a mathematical model of the simulated reaction. The use of a resistance means that no chemical reaction takes place in the reactor, instead the reaction is emulated on basis of temperature changes, as done in [15]. B. Mathematical model Although it is not necessary to have a mathematical model for the design of the min-max predictive controller, this section shows the process model to emphasize its nonlinear character. The mathematical model also justifies the way to emulate the heat generated by the chemical reaction with the aid of the resistance. The emulated chemical reaction, representing a refinement process, was used previously in [6]. With Ff=Ff,in =Ff,out and a constant volume, the model of the chemical reaction can be defined as: dT dt=−Fj V(Tj,in −Tj,out) +(−∆H)·V MCp k0e−E/(RT)C2 A(16) dCA dt=Ff V(CA,in−CA)−k0e−E/(RT)C2 A(17) denoting Fj,Tj,in and Tj,out the flow rate through the jacket and the temperature of the water entering and leaving the cooling jacket, respectively. CAand CA,in represent the reactive concentration in the reactor and in the feed, respectively. The feed passes through the heat exchanger and enters the reactor nearly with the temperature of the reactor content. Thus it is assumed that no heat is neither removed or supplied due to the feed. The heat exchange in the cooling jacket is given by the following empirical model: Fj·(Tj,out −Tj,in)=T− α β (1−e− γ Fj)(18) with α =292.19K, β =14.94s/land γ =13.18s/l. The chemical reaction is nonlinear in the dynamics of the temperature and the concentration due to the quadratic terms of the concentration in the model equations (16) and (17). For further details on the model parameters see [6]. V. EXPERIMENTAL RESULTS The strategy described in section III has been applied to the refinement process. In this section the experimental results will be exposed and discussed. CARIMA type prediction models with bounded additive uncertainties were used in the experiments. This type of model extends the concept of noise in traditional CARIMA models so that an uncertainty is considered: A(z−1)y(t)=z−dB(z−1)u(t−1)+C(z−1) θ (t) ∆(19) Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. with ∆=1−z−1, θ (t)∈{ θ ∈Rdimy :∥ θ ∥∞≤ ε },anddimy the dimension of y(t).Theuseofthistypeofprediction models results in a control law without error in steady state. There are few differences between implementing the algorithm of section III for a state space model and a CARIMA model with bounded additive uncertainties. The main difference is the method used to find the matrices of the prediction equation [3]. The cost function is the same as in (6). Appendix A presents practical aspects related to the use of a CARIMA model (19) for prediction. In the following sections the control system in the pilot plant will be described and the necessary steps to obtain apredictionmodelwillbepresented.Finally,experimental results will be exposed. A. Description of the control system The sensors and actuators in the plant are connected to a PMC-10 control unit. The PMC-10 is connected by ARCnet to a personal computer that runs the control and monitoring system Simatic-IT.Thecontrolalgorithmhas been implemented directly in Matlab and the communication with Simatic-IT is done using the OPC protocol (OLE for Process Control). Both Simatic-IT and the controller run on the same personal computer, based on a Pentium II processor at 300 Mhz. This computer does not have enough computational power to solve exactly the min-max problem of a typical MMMPC, but can compute the control action using the proposed strategy. B. Identification of the prediction model APRMSS (Pseudo-Random Multilevel Step Sequence)has been applied to the recirculation valve with the objective of collecting data for the parameter identification of the prediction model. The periods of the PRMSS have been chosen sufficiently long to observe the reaction of the pilot plant to changes in the input (see Fig. 3). It can be seen that the temperature of the tank reaches steady state in each step in something more than two hours, although the variations in steady state are of several degrees. The reagent concentration also suffers variations in steady state. It can be observed that the input–output gain is negative and clearly variable (greater gain for low openings of v8). A first order transfer function model with delay is proposed as prediction model. This low order model cannot correctly describe the dynamics of the plant, but it is a good approach to check the robustness of the controller in presence of uncertainties and disturbances. Using the data of Fig. 3 the following model has been identified: G(s)= −0.975 950s+1e−31.25s(20) This model was discretized with a sampling time of Ts=60s. The delay was rounded to 1 sampling time in order to avoid approximations of the time delay, e.g. Pad´e approximation. Thereby, the following CARIMA model was obtained: y(t+1)=0.939y(t)−0.0597u(t−1)+ θ (t) ∆(21) with the noise polynomial C(z−1)=1. v8[%] CA[mol/l] T[oC] t[min] 0 0 0 100 100 100 200 200 200 300 300 300 400 400 400 500 500 500 600 600 600 30 40 40 50 60 60 70 80 80 0 0.1 0.2 0.3 0.4 Fig. 3. Experiment for the prediction model identification. From top to bottom: Tank temperature (T), valve opening (v8)yreagentconcentration (CA). error t[min] 100 200 300 400 500 600 −0.5 −0.25 0 0 0.25 0.5 Fig. 4. One step ahead prediction error during the experimentforthe model identification. C. Experimental results of the controller The proposed control strategy was applied to the pilot plant described in section IV-A using (21) as a prediction model. For the prediction and control horizons values of N= 15 and Nu=12 have been used1.Thereforetheprediction horizons includes approximately one time constant of the process, a common value for this parameter in predictive control. The weighting factor for the control effort has been chosen equal to Rj=2. Based upon the one step ahead prediction error (see Fig. 4) the parameter ε has been chosen to ε =0.25. As a result, in 97% of the samples the one step ahead prediction error is bounded by the chosen value. Finally, in order to restrict the system input and output in the experiments, the following constraints have been used: 30 ≤ˆy(t+j|t)≤70,j=2,...,16,∀ θ θ θ ∈vert(Θ Θ Θ) 5≤u(t+j|t)≤100,j=0,...,11 −20 ≤∆u(t+j|t)≤20,j=0,...,11 Note that in the output restrictions the effect of the uncertainty has to be considered. In order to analyse the system behaviour, several experiments with reference changes and disturbance rejection have been made using the proposed control strategy. Fig. 5showstheresultsofthetrackingexperimentwithref- erences different enough to result in control actions in a large interval. After the first reference change no overshoot 1The model delay implies that the prediction horizon starts int+2and ends in t+16. Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. v8[%]T[oC] mol/ t[min] 0 0 150 150 35 40 40 45 50 50 50 55 60 60 65 20 80 100 100 100 TT2[oC]CA[mol/l] t[min] 0 0 150 150 50 50 100 100 0.1 0.15 0.2 0.25 0.3 17 17.5 18 18.5 Fig. 5. Reference tracking experiment. From top to bottom: Tank temperature (T), valve opening (v8), reagent concentration (CA)andcold water temperature (TT2). appears in spite of a quite fast controller reaction. After the second reference change a small overshoot (of about −0.5oC), justified by the nonlinear process behaviour, can be observed. In steady state the controller shows small changes in the control action necessary to stabilize the output on the reference in presence of variations in the generated heat and the cold water temperature. In second place a disturbance rejection experiment was carried out. The results of a disturbance in the system input, the opening of the valve v8,arepresentedinFig. 6. As can be seen, after approximately 70min a constant disturbance in the input of ∆v8=15% was applied. The controller reacts rapidly and rejects the perturbation in about 20 minutes. After the disappearance of the perturbation in t= 101min the controlled system shows the same behaviour and reaches steady state in approximately 20 minutes. Neither the temperature nor the control action show oscillations after the perturbation. The third experiment, using an additive disturbance in the feeding Ff,isshowninFig.7.Int=70min a change in the feeding flow of ∆Ff=0.0125l/s, which corresponds to an error of 25%, has been applied. With an increasing error, the controller reduces the opening of the valve and reaches acompensationofthedivergenceafter15minutes.Inthis experiment an overshoot of −0.50oCcanbeobserved.The oscillation in the temperature and the control action is quite small and seems acceptable due to the strong disturbance. The reference tracking experiment was repeated with a linear constrained predictive controller (GPC) to allow the comparison between the proposed strategy and a standard MPC method. The GPC is based on the linear model (21) v8[%]T[oC] t[min] 0 0 120 120 40 40 40 45 50 50 55 60 60 60 60 65 70 20 20 80 80 100 100 30 u[%] TT2[oC]CA[mol/l] t[min] 0 0 0 120 120 120 40 40 40 40 50 60 60 60 60 70 20 20 20 80 80 80 100 100 100 0.05 0.1 0.15 0.2 17 17.5 18 18.5 30 Fig. 6. Experiment with input disturbance rejection. From top to bottom: Tank temperature (T), valve opening(v8), controller output (u), reagent concentration (CA)andcoldwatertemperature(TT2). v8[%]T[oC] t[min] 0 0 40 40 40 40 45 50 50 55 60 60 60 60 65 70 20 20 80 80 80 100 100 30 TT2[oC]CA[mol/l] t[min] 0 0 40 40 60 60 20 20 80 80 100 100 0.05 0.1 0.15 0.2 0.25 17 17.5 18 18.5 Fig. 7. Experiment with disturbance rejection in the feed flow. From top to bottom: Tank temperature (T), valve opening(v8), reagent concentration (CA)andcoldwatertemperature(TT2). Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. MMMPC MPC ref MMMPC MPC v8[%]T[oC] t[min] t[min] 20 80 100 35 40 40 45 50 55 60 60 65 0 0 25 25 50 50 75 75 100 100 125 125 150 150 Fig. 8. Reference tracking results of the MMMPC and the GPC. From top to bottom: Tank temperature (T), valve opening (v8). and was used with the same parameters as the MMMPC. It can be observed in the results (see Fig. 8) that the process controlled by the GPC exhibits significant oscillations in the temperature and the control action after the reference changes. The comparison of the results shows that the MMMPC stabilizes the temperature more efficiently and with less oscillations in the opening of the valve. Finally, it is important to mention that the calculation of the control signal took place without problems within the chosen sampling time (60 seconds). During the experiments the average computation time was 5.64 seconds, with a maximum of 9.90 seconds and a minimum of 1.863 seconds. VI. CONCLUSIONS In this paper an MMMPC based on an tractable QP problem was applied to a pilot plant. The results showed a good system behaviour and the stabilisation of the plant temperature around the operation point. After reference changes the controller quickly compensates deviations. Furthermore, the MMMPC showed its capacity to compensate errors caused by the disturbances. The application to a process shown in this work joins the small number of MMMPC applications reported in specialised literature. The low computational requirements of the proposed control strategy allowed the use of appropriate sampling times and realistic prediction and control horizons. Thereby it is shown that the use of proposed strategy allows the application of this kind of controllers to a larger number of processes. REFERENCES [1] T. Alamo, D.R. Ram´ırez, D. Mu˜noz de la Pe˜na, and E.F. Camacho. Min-max MPC using a tractable QP problem. Automatica,43:693– 700, 2007. [2] A. Bemporad, F. Borrelli, and M. Morari. Min-max Control of Constrained Uncertain Discrete-Time Linear Systems. IEEE Transactions on Automatic Control,48(9):1600–1606,2003. [3] E.F. Camacho and C. Bord´ons. Model Predictive Control.Springer- Verlag, second edition, 2004. [4] P.J. Campo and M. Morari. Robust Model Predictive Control. In Proc. American Control Conference,pages1021–1026,June10-121987. [5] D. Mu˜noz de la Pe˜na, D.R. Ram´ırez, E.F. Camacho, and T. ´ Alamo. Application of an explicit min-max mpc to a scaled laboratoryprocess. Control Engineering Practice,13:1463–1471,2005. [6] J.K. Gruber and C. Bordons. Control Predictivo no Lineal Basado en Modelos de Volterra. Aplicaci´on a una Planta Piloto. Revista Iberoamericana de Autom´atica e Inform´atica Industrial (in spanish), 4(3):34–45, 2007. [7] E.C. Kerrigan and J.M. Maciejowski. Feedback min-max Model predictive Control Using a Single Linear Program: Robust Stability and the Explicit Solution. International Journal of Robust Nonlinear Control,14:395–413,2004. [8] Y.H. Kim and W.H. Kwon. An Application of Min-Max Generalized Predictive Control to Sintering Processes. Control Engineering Practice,6:999–1007,1998. [9] M.V. Kothare, V. Balakrishnan, and M. Morari. Robust constrained model predictive control using linear model inequalities. Automatica, 32(10):1361–1379, 1996. [10] J.H. Lee and Zhenghong Yu. Worst-case formulations of model predictive control for systems with bounded parameters. Automatica, 33(5):763–781, 1997. [11] Y. Lu and Y. Arkun. Quasi-Min-Max MPC Algorithms for LPV systems. Automatica,36(4):527–540,2000. [12] D.Q. Mayne, J.B. Rawlings, C.V. Rao, and P.O.M. Scokaert. Constrained model predictive control: Stability and optimality. Automatica, 36:789–814, 2000. [13] D.R. Ram´ırez, M.R. Arahal, and E.F. Camacho. Min-Max Predictive Control of a Heat Exchanger using a Neural Network Solver. IEEE Trans. on Control Systems Technology,12(5):776–786,2004. [14] D.R. Ramirez and E.F. Camacho. Piecewise Affinity of Min-Max MPC with bounded additive uncertainties and a quadratic criterion. Automatica,42:295–302,2006. [15] L.O. Santos, P.A.F.N.A. Afonso, J.A.A.M. Castro, N.M.C. Oliveira, and L.T. Biegler. On-line implementation of nonlinear MPC: an experimental case study. Control Engineering Practice,9(8):847–857, 2001. [16] P.O.M. Scokaert and D.Q. Mayne. Min-max feedback model predictive control for constrained linear systems. IEEE Transactions on Automatic Control,43(8):1136–1142,1998. [17] F. Szeifert, T. Chovan, and L. Nagy. Process dynamics andtemperature control of fed-batch reactors. Computers & Chemical Engineering, 19(1):447–452, 1995. APPENDIX A. PRACTICAL ASPECTS RELATED TO THE INPUT/OUTPUT DESCRIPTION This appendix presents the main differences between the implementation of the proposed strategy with a state space model (1) and an input/output model like (19). In this case, the evolution of the output along the prediction horizon can be described in condensed form as [3]: y=Guu+G θθ θ θ +Fxx(22) The matrices Gu,G θ and Fxcan be obtained from the original model using several different methods, e.g. the Diophantine equation [3]. A more intuitive method is the use of the step response coefficients [3] to form matrix Gu.Thefree response vector Fxxcan easily be computed by iterating with the process model with the assumption of a constant control signal along the prediction horizon. Matrix G θ is computed by treating the uncertainty as an additional system input weighted by a unit polynomial in the model. Thus, using the method of the Diophantine equation, G θ can be computed as Gu,butassumingthatB(z−1)=1. Also, G θ can be formed from the coefficients of the system response to a step in the uncertainty θ (t).Thevectorucontains the future control increments if a CARIMA model is used [3]. In this case state vector xcontains the present and past output values as Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. well as the past control increments: x=[y(t),...,y(t−na),∆u(t−1),...,∆u(t−nb)]T denoting naand nbthe polynomial orders of (1−z−1)A(z−1) and B(z−1),respectively. The different model implies changes in the cost function, but, as shown in the following, it can be rewritten as in (6). The usual form of the cost function with a prediction model like (19) is: V(x,u, θ θ θ )= N ∑ j=1 e(t+j|t)TQe(t+j|t)+ Nu−1 ∑ j=0 ∆u(t+j)TRm∆u(t+j)(23) being e(t+j|t)the predicted reference tracking error at time t+j: e(t+j|t)=(y(t+j|t)−w(t+j)) where w(t+j)is the reference at time t+j.Withthe prediction equation (22) and the assumption Q=1thecost function can be written as: V(x,u, θ θ θ )=(Guu+G θθ θ θ +Fxx−w)T ·(Guu+G θθ θ θ +Fxx−w) +uTRu(24) with Radiagonalmatrixoftheform: R=diag(Rm,...,Rm)(25) The vector wcontains the future values of the reference trajectory and is defined as w=[w(t+1),...,w(t+N)]T. If a zero reference is used, i.e. w(t+j)=0for j=1,...,N, the cost function can be expressed as in (6) being the matrices Muu =GT uGu+R,M θθ =GT θ G θ ,M θ u=GT θ Gu, Muf =GT uFx,M θ f=GT θ Fxand Mff =FT xFx.Inthecaseof anonzeroreference,thecomputationofthecostfunction can be carried out in an analogous way and does not change the optimisation algorithm. From this point, the proposed strategy can be applied in the same way as with state space models. Preprint submitted to 47th IEEE Conference on Decision and Control. Received March 10, 2008. View publication statsView publication stats