Model predictive control with soft constraints with application to lime kiln control

Computers and Chemical Engineering 23 (1999) 791 – 806

Model predictive control with soft constraints with application to lime kiln control Rolando Zanovello, Hector Budman* Department of Chemical Engineering, Uni6ersity of Waterloo, Waterloo, Ontario, Canada Received 23 May 1996; received in revised form 24 November 1997; accepted 27 January 1999

Abstract A model predictive control algorithm with soft constraints is presented. The technique, finite number of weights model predictive control (FNWMPC), is based on the selection of a proper combination of weights when constraints are violated. Robust stability of the FNWMPC technique is studied using m-structured singular value techniques. The algorithm is tested through simulations of a lime kiln control problem and is compared to other existing model predictive control techniques. © 1999 Elsevier Science Ltd. All rights reserved. Keywords: Lime kiln; Model predictive control; Soft constraints

1. Introduction This paper deals with the robustness and constraint properties of model predictive control (MPC), which is a widely accepted control algorithm for multivariable systems with constraints. Although the techniques and comparisons presented in this work are general and applicable to any multivariable process, a lime kiln was selected as a specific case study to test our methodology. The kiln is an ideal example to test our proposed technique since the issues of constraint handling and robustness are essential for the success of the control strategy. MPC was previously applied to a lime kiln process by Charos, Taylor and Arkun (1991). However no specific problems regarding constraints handling or robustness were reported in that work. In our lime kiln application, infeasibilities in the presence of constraints and issues of sensitivity to model mismatch occurred during normal operation of the kiln and had a clear impact on the performance of the control system. This motivated us to search for modifications of existing MPC techniques which can deal more effectively with the aforementioned problems. * Corresponding author.

Different alternatives exist to deal with infeasibilities of the resulting constrained optimization problem. For example, it is common practice to select a constraint window (Morari, 1989) which is a specific number of steps in the horizon for which the corresponding constraints are being removed. This is especially effective when large time delays are associated to specific output variables. However, this is an ad hoc method which does not satisfy any optimal criterion. A second option to deal with infeasibilities is to soften the constraints. This approach was adopted in the present work since it can be done systematically in contrast with the ad hoc constraint window approach. Some work on this topic was recently reported (Zafiriou, 1990; De Oliveira & Biegler, 1994; Feher & Erickson, 1993). The idea of constraint softening is to allow the controlled variables to violate their constraints by small amounts, whereas no constraint violations are allows in the hard constrained case. Constraints may be imposed on manipulated, controlled or associated variables. Associated variables are those which have to be maintained in between limits but are not required to be controlled at a specific target value. Violation of constraints by small amounts will generally be accepted in the chemical process industry, especially when dealing with constraints on the output or associated variables which will often temporarily deviate from their setpoint

0098-1354/99/$ - see front matter © 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 9 9 ) 0 0 0 0 8 - 3

792

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

or prespecified limits. For example, the oxygen level in the exit gases from the kiln should be maintained above a concentration limit of 1% O2 due to environmental restrictions. However, deviations of 10–15% from this limit (i.e. 0.85 – 0.9% O2) are still considered acceptable for safe process operation. When dealing with softening of constraints, a basic question to answer is to what degree of violation could be allowed for each one of the controlled or constrained variables during process operation. This is especially critical for nonsquare processes such as the lime kiln where the number of outputs and constrained variables is larger than the number of manipulated variables. For these type of processes, situations may arise where offsets will occur in all of the outputs and constrained variables due to the insufficient number of degrees of freedom in the process. However, in practice, some of the controlled variables or constrained variables have a larger effect than others on the final quality of the product. For instance, in the lime kiln, the hot end temperature has a larger impact on the lime conversion than the cold end temperature. Therefore different variables should be assigned different constraint compliance priorities in the control strategy. One way to solve this problem is by adding slack variables into the cost function to be minimized (Ricker, Subrahmanian, & Sim, 1989). In this technique a scalar error term is added which represents the maximum weighted output constraint violation to the cost function. However, no action is being taken to correct for other output variables which may have significant violation since the error term only accounts for the largest constraint violation. In addition, in this method, it is not easy to analyse the system for stability and robustness. Another natural way to assign this relative priority in MPC is by assigning different weight values to the controlled or associated variables in the cost function to be minimized should violation occur. This will be the basis for the soft constrained technique presented in this work. Two different approaches were previously proposed to assign weights for constraint softening. In the first approach proposed by Feher and Erickson (1993), to be referred to as the FE algorithm, weight values are assigned to each variable and are maintained constant in time. In the second approach proposed by Kouvaritakis, Rositer and Chang (1992), the weights are allowed to change with time. The key disadvantage in the first approach is that the weight values are fixed. Thus, the fact that different combinations of disturbances and setpoints occurring during process operation which may affect the outputs differently, cannot be taken into account for the selection of the weight values. This problem is directly tackled by Rossiter where the weights are continuously adjusted with

time. Although, Rossiter developed his approach for the hard constrained case, it results into a soft constrained algorithm whenever infeasibilities are encountered. Rossiter’s algorithm, CSGPC, has guaranteed stability properties; however, in the case of input constraints stability can only be guaranteed under certain feasibility assumptions (Rossiter & Kouvaritakis, 1993). In the case of nonsquare systems, Rossiter’s technique will minimize a maximal constraint violation condition by using a special optimization technique. This optimization has to be implemented in conjunction with an algorithm, stable generalized predictive control (SGPC) which forms a nominally stable closed loop. Due to these, the CSGPC algorithm cannot be easily integrated into conventional unconstrained model predictive controllers commercially available. This was an important consideration for the industrial application considered in this study since it facilitated eventual implementation of the strategy. In this work we will show a technique which is an extension of the FE technique and which blends together some positive characteristics from the methods mentioned above. Finite number of weights model predictive controller, to be referred as FNWMPC can be easily integrated into a conventional unconstrained model predictive controller based on a least squares solution of the cost function. In the proposed method the weights are permitted to change with time following different set point changes or disturbances occurring in the process. Each weight in the cost function can assume only a finite number of values, assigned off-line. Consequently, in this work we expand the FE algorithm and we also address the stability and robustness of the new algorithm which were not previously studied by Feher. The nominal stability of the system, is tested off line for each possible combination of weights. Robustness is an additional important topic to be considered for the successful application of the control strategies. It will be shown in this paper that our technique can be cast into a form suitable for analysis with m-structured singular value theory. Finally, our technique is compared to some existing soft and hard constraint algorithms for a lime kiln case study and conclusions will be drawn. This paper will be organized as follows: In Section 2 we will present the FNWMPC technique. In Section 3 we will discuss nominal and robust stability of the algorithm. Subsequently we will present the application of the methods presented in the previous sections to the lime kiln problem in Section 4. Comparisons between our technique to quadratic dynamic matrix control (QDMC) and the FE (Feher & Erickson, 1993) technique will be shown in Section 5.

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

2. Soft constrained model predictive control algorithm In this section we will present a new technique for model predictive control with soft constraints. Consider a linear multivariable system with nu inputs and ny outputs to be controlled with a model predictive controller using a prediction horizon n and a control horizon p. For simplicity, we consider throughout this work that one future input move is computed at each time interval. The typical cost function to be minimized by the model predictive controller is:

*

n*

−

V = G Y (k+1 k)−R(k +1 k)

+ LDu(k) 22

(1)

Where G and L are diagonal weight matrices, Y, Du and R are the p-step ahead prediction vector, the manipulated variable move and the desired reference trajectory respectively as defined in Appendix A. In most available MPC algorithms the weights in this function are selected off line and are kept constant during the operation of the controller. L is called the suppression move factor and its purpose is to penalize large moves in the manipulated variables while G is known as the output weights vector and it is used to prioritize one controlled variable over another. An additional possible utilization of weights in the cost function is for the handling of soft constraints (Feher and Erickson, 1993). In this approach, whenever a variable violates its constraint limit, the set point for this variable is reset to this limit. If the constraint violation occurs in an associated variable, this variable is automatically treated as a controlled variable with set point equal to the constraint limit. Following this reasoning, whenever a constraint violation occurs, the cost function is changed accordingly to:

*

−

−

V% = V+ pj Y (k+ 1/k) − Y cj

n*

2

(2) 2

where V is defined in Eq. (1), Ycj is the constraint limit and pj is the weight which penalizes the jth-output value in the horizon. The last term in Eq. (2) was added to penalize the constraint violation occurring in controlled or associated variables. The value of the weights pj will determine the amount by which the manipulated, controlled or associated variables violate their respective constraint limits. For instance, if one of the weights pj is set to infinity, then a hard constraint is imposed on the corresponding variable. On the other hand, setting one of these weights to zero implies that the constraint is being removed from the problem. The soft constraint situation corresponds to cases where the weights are set to values between zero and infinity. This idea is the starting point for the MPC algorithm proposed in this work. The algorithm will be referred to as FNWMPC.

793

2.1. Finite number of weights model predicti6e controller This algorithm is an extension of the FE technique (Feher & Erickson, 1993). The algorithm consists of minimizing the following objective function: min V[u(k)]+P[u(k)] u(k)

(3)

V is the cost function in Eq. (1), P is a penalty function defined below and u is the manipulated variable value. The principal difference between the proposed FNWMPC technique and the FE method is in the penalty function P which is now supplemented with additional penalty terms as follows: nu

P[u(k)]= % [pILi (u)(e Ti u− aLi )2 + pIUi (u)(e Ti u−aUi )2] i=1

nl

no ·p

+ % %



t=1 j=1

 

−

pOLjt e Tj Y (k+1/k)− bLjt

−

+ pOUjt e Tj Y (k+1/k)− bUjt

n



2

2

(4)

In the FE algorithm, the penalty function P(u) does not include the double summation term in Eq. (4). The vector elements ei and pIj are given as: e T1 = [10···0] e T2 = [01···0]

···etc.

(5)

and: pIb = pIb if aLjt or aUjt are violated, pIb = 0

no violation

b=Li, Ui (6)

The double summations represent the fact that different levels or ranges of variables deviating from their constraint limits are penalized with different weight values, pOb. These weights are now defined as follows: pOb = pOb if bLjt or bUjt are violated and 15t5 nl, b= Ljt, Ujt pOb = 0

no violation

(7)

There are now n1 pairs of upper/lower limits per variable as compared to one single pair considered in the FE technique. The advantage obtained with the addition of the terms in the double summation in Eq. (9) is that now different weights can be assigned to different variable ranges according to their relative importance in the process. The values of the weights pOb and pIb have to be specified by the user and they are used both to scale the variables and to prioritize some constraint over others. For example, consider a controlled variable x with three possible levels of constraints, i.e. n1 =3, and where these levels are defined as follows: x 1 5 x5x2 x 2 5 x5x3 x 3 5 x5x4

p Ob = w1 p Ob = w2 p Ob = w3

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

794

The weights can be equal or have different values. However, in general, since it is desired to minimize constraint deviations, larger weight values should be applied to the levels furthest away from the lowest constraint limit x1, i.e. we select w3 \w2 \w1. For example, assume that the current value xj in the horizon of x lies in between x1 and x2, then according to the definitions in Eqs. (6) and (7) above: pOUj 1 =w1,

pOUj 2 =0

and pOUj 3 =0

Thus, the first level of constraints (t = l) is active and a corresponding term will be added in the cost function while the second level and third level (t =2, 3) are inactive and will not contribute any terms to the cost function. Note that the FNWMPC methodology improves over the FE technique since the number of weights is larger than for the FE technique allowing greater flexibility in dealing with constraints. Each possible combination of weights in the FNWMPC will be referred to hereafter as component Ji. It should be noticed that each component may correspond to constraint violation in either only one or more variables. Although the optimization problem is nonlinear, the cost function has a unique optimum. The proof readily follows from the fact that the function V(u) is convex and P(u) is convex, hence, the sum V(u) + P(u) is also convex, and the solution to the minimization problem will be unique (Fletcher, 1987). The key difficulty to compute the global minimum of the cost function V(u) +P(u) is that the weights are a function of the solution which is not known a priori. A method originally proposed by Feher was also adopted in this work to solve this problem. It consists of a systematic search of the minimum through all the components starting from the unconstrained component. This search may be summarized by the following steps: Step 1. Compute the unconstrained solution. If implementation of this solution produces prediction which do not violate constraints, the search is stopped and the solution is accepted. Otherwise, if violation is predicted, proceed to step 2. Step 2. Determine in which component the prediction computed in step 1, which was found to violate constraints, resides. Verify that this component was not implemented before. If this is true calculate the weights corresponding to this component and proceed to step 3. However, if this component was previously selected, ignore it and assume a new component which has not been previously used. Compute the weights for this component and proceed to step 3. Step 3. Compute the least squares solution with the weights calculated in step 2. If the calculated solution resides in the last assumed component then the iteration converged and the solution has been found. If not, the iteration is repeated from step 2.

Since this search is conducted until a solution is found, checking all components if necessary, convergence to the global minimum is guaranteed.

3. Nominal and robust stability The technique presented above, produce, at each time interval, a linear system (controller+plant). However, the implemented closed loop system is overall nonlinear since it is composed of a series of different linear systems. Zafiriou (1990) developed sufficient conditions for robust stability of the QDMC algorithm which solves a quadratic objective function with hard constraints. In that work it was shown that QDMC also results in a linear system at each time interval Therefore, the stability conditions proposed by Zafiriou for piecewise linear systems can be also applied to the FNWMPC method presented in this work. Following Zafiriou’s notation we define an operator which maps the vector of states of the system at time k to the vector of states at time k+ 1 as follows: X(k+ 1)= FJi,nom[X(k)], i= 1, . . . , nc

(8)

where Ji is the current active component, Fji ,nom is the linear operator describing the system at the current time and its derivation is shown in Appendix A—the subscript nom indicates that the operator is computed using the nominal impulse response coefficients of the process—nc total number of possible components Ji. For the case that the manipulated variable u(k) is selected as the state x(k), Zafiriou show that a sufficient condition for asymptotic stability is given by:

F%Ji B 1,

Öi,

i= 1, . . . , nc

(9)

A more easily computable condition is given by: (10)

DF%Ji, nomD − 1 B 1, Öi, i= 1, . . . , nc ’ where is a consistent matrix norm and D is a scaling matrix. If this condition is satisfied, Fji ,nom is said to be a contraction. The difficulty in satisfying this equation resides in the fact that the same norm with the same scaling matrix D has to be applied to all components and has to be smaller than one. The numerical example in Section 5 will clearly show that in a practical problem this condition is difficult to satisfy and may lead to very conservative results. To avoid conservativeness, a less stringent test of stability is desired in practical applications. Since the spectral radius r is a lower bound to the matrix norm in Eq. (9) a sufficient condition for instability of the overall nonlinear system is: r(F%Ji,nom)\ 1,

for any i,

i =1, . . . , nc

(11)

When this condition occurs a local linear system (controller+ plant) is unstable, and, consequently the overall nonlinear system is also unstable.

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

795

Fig. 1. Schematic block diagram of the MPC controller and the connection matrix M.

For the FNWMPC algorithm this condition is tested off line for all possible components Ji, since each component and consequently each F operator, are known a priori. The block diagram in Fig. 1 represents the controller and the process subdivided as a nominal model with uncertainty. This diagram can also be represented by the feedback connection in Fig. 2 where the function M is indicated in Fig. 1 and calculated in Appendix B. Using the definition of the FJi operator in Appendix A and the definition of M, the input vector U(k) selected as the system’s state vector evolve via the following linear difference equation:

has to be satisfied for all the family of plants which define the uncertainty set, H. Clearly, this condition is more conservative than the nominal stability test which is conducted using the nominal impulse response model A less conservative yet weaker condition for robustness is to verify a necessary condition for robust stability based on the spectral radius, which is a lower bound of the norm in Eq. (14):

U(k +1)= FJi,true(M, Dk ) U(k)

m(MJi )B 1,

(12)

FJi ,true is the operator calculated with the true impulse response model coefficients of the process, Dk is the uncertainty in the model for each time k. The following assumptions are made:

r(F%Ji,true)B 1, sup H

Öi,

i= 1, . . . , nc

From the results of the main loop theorem (Doyle and Packard, 1988) the following equation can be shown to be equivalent to the inequality in Eq. (15): Öi,

i= 1, . . . , nc

(13)

Then from Zafiriou’s work, a sufficient condition for robust stability is:

F%Ji,true B1, sup H

Öi,

i= 1, . . . , nc

(14)

where H is the family of impulse response models which represents the true plant and are identified from experiments (Webb, Budman, & Morari, 1995). Condition (14) is the nominal robust stability condition which

(16)

where M is given in Appendix B and m is the structured singular value of M computed with respect to the following uncertainty structure:

i s¯ (Dk )5 1 ii Dk is fixed, does not vary with time

(15)

Fig. 2. M– D feedback connection structure.

796

D1 =

!

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

d 1In s 0m*n s

n

"

0n *m s : d1  C, D  C m*m  C (ns + m)*(ns + m) D (17)

(i.e. number of states) and di 51. where ns = n *n u This theorem implies that the m test has to be performed for every component which may become active during system operation. The m calculation is performed using MUTOOLS (Matlab toolbox). The uncertainty description used in this work and given in Appendix B in conjunction with Fig. 1 is a diagonal matrix of the form: D =[d1, d2, . . . , dn*nu* ny ]. The corresponding M matrix is a square matrix with dimensions (ns + n*n*u ny ) by (ns +n*n*u ny ) and can become very large. For example, for the case study investigated in this work (n = 40, nu =2 and ny =3) the M matrix is typically 320 by 320. It was found that the m computation for such matrices was slow and potentially inaccurate. For this reason we use the frequency domain counterpart of theorem (21) given by (Doyle and Packard, 1988):

m[G(v)] B1

Öv

(18)

where m is computed with respect to D and where G(v) =M22 + M21(e jvTI −M11) − 1M12 The conversion of the M matrix from the time/state space representation to the frequency domain representation is conducted using the function FRSP in the MUTOOLS toolbox mentioned above. This transformation enables m computations for equivalent matrices with the advantage of working with lower dimensions. For example, the 320 by 320 matrices mentioned above are now of dimensions 240 by 240. If correlations in the uncertainty elements are assumed, as will be explained below, the dimension of these matrices are further reduced to 6 by 6. The disadvantage with this approach is that since the matrix elements change as a function of frequency, the calculation has to be repeated at each frequency in the relevant range. The m tests presented in this section have to be conducted for each one of the components Ji that become active during operation. All possible components are known a priori for the FNWMPC, enabling the m computation to be performed off line for all the possible components.

4. Case study—lime kiln control The lime kiln is an endothermic reactor where calcium carbonate is being continuously converted to calcium oxide by means of the heat produced by a large flame at the front end of the kiln bed. The kiln is one of the most energy intensive processes in a pulp mill. Therefore, tight control of this process would result not only in quality improvement, but also in significant

operating cost savings. The quality of the product is determined by the amount of conversion of CaCO3 to CaO which is a direct function of the exit (or front end) temperature. It, therefore, becomes imperative to maintain the front temperature at such a level that complete, or almost complete, conversion will be ensured. On the other hand, overly high temperature may result in damage to the refractories lining the inner surface of the reactor vessel. Control of the cold end temperature is also important since low values of this variable may result in agglomeration of the feed material in large balls. This will impede uniform exposure of the material to the external heat and consequently lower conversion may result. High inlet temperatures may result in damage to the system of chains located at the inlet, which are installed to increase the overall heat transfer between the hot induced air and the feed. A third variable to be controlled is the oxygen concentration in the exit gases. Excessive oxygen levels indicate that fan power is being wasted. On the other hand, low oxygen level indicate the presence of total reduced sulphur (TRS) especially H2S which is considered to be a toxic waste. Therefore specific levels of TRS, and concurrently of oxygen, cannot be exceeded. Based on the process description it is desired to control the cold and front end temperatures at specific values while allowing the oxygen to vary between two prespecified limits. Given the definitions described in the introduction, it follows that the end temperatures are the controlled variables while the oxygen is an associated variable. The two manipulated variables are the flow of air through the kiln induced by a pump located at the cold end and the flowrate of natural gas used for combustion. The principal disturbance to the process is the change of the lime feed located at the cold end section of the kiln. The FNWMPC method presented in Section 2 was used to design different controllers for the system. The transfer functions used for model prediction were identified from a pulp mill located in Prince George, BC, Canada. The transfer matrix for the kiln is given by:

Æ 80 e − 2s Ç 1.1 e − 15s à à T 152.5s+1 (40s+ 1)(50s+ 1) Æ hot Ç Ã Ã u ID ÃT coldà = à − 13.65 e − 2s 0.3 e − 2s à u NG È O É Ã 6.5s+1 à 100s+1 2 Ã Ã É È −0.72 −0.0075

 n

Æ 0.65 e − 149s Ç Ã Ã 25s+ 1 à − 0.785 e − 25sà +à à [dMUD] à 82.5s+ 1 à − 25s à − 0.011 e Ã È É 33s+ 1

(19)

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

where the variables are denoted as: Thot, Tcold for hot and cold end temperatures, O2 for oxygen, uID, uNG for air and natural gas flowrates and dMUD for the mud (feed) disturbance. Typical operating conditions, are Thot =1000°C, Tcold =200°C, 0.6BO2 B5%, uNG =1100 m3 h − 1, uID =3 mm Hg pressure (induced vacuum), MUD (feed) =800 ton day − 1, percent CaO recovery= 90%. The air flowrate is regulated through the vacuum induced by the fan; uID is given in terms of this pressure. A 3 mmHg pressure corresponds approximately to 2500 ton day − 1 of air. In order to test the stability and robustness properties of FNWMPC, a first controller, to be referred in the discussion as FNWMPC-A, was designed. This controller uses a control horizon p =15, a prediction horizon n= 40 and a sampling interval of 5 min. The constraint limits and the corresponding chosen weights were as follows: −0.4% \ O2 \ − 0.6% O2 B − 0.6 Thot B20°C

pOb =50

pOb =200 pOb = 3

(20)

All constraint values are given in deviations from the initial steady state. The local nominal stability of the FNWMPC-A controller was tested using the spectral radius test (Eq. (11)) off-line for all possible components nc. Regarding nc, in principle, a penalty weight value could be assigned to the error of a variable at each interval in the horizon for which it is found to violate either its upper or lower constraint. Consider, for example, a control horizon p of 15, and a given variable has both an upper and lower limit. The resulting number of components or possible combinations along the horizon corresponding to upper limit active, lower limit active or constraint inactive is 315. Consequently, the iterative procedure required to find the solution for the FNWMPC controller may become, at a given time interval, very lengthy. Moreover, the nominal and robust stability tests also become time consuming since they must be conducted for each one of the possible components. To simplify the simulations for the purpose of this paper, only a smaller subset of four weights in the horizon was considered and all other weights were assumed to be zero. The location of the nonzero weights is important since it will determine which errors in the horizon are considered in the optimization. The response in oxygen was not found to be sensitive to the location of the weights in the horizon due to the fact that oxygen responds very fast to changes in either manipulated variables or disturbances. On the other hand, it was found in additional simulations that the location of the weights significantly affected the response of the system when constraints

797

in Thot or Tcold were considered. Since the minimization of steady state offsets was of primary importance for the kiln, it was found that the best location of the nonzero weights is at the last section of the prediction horizon of Thot and Tcold, (i.e. close to the steady state situation). Considering that only 4 nonzero weights were used for each constraint level throughout the horizon, the total number of components amounted for our case study to 1296. As shown in Fig. 3, the spectral radii for all these components are smaller than one, indicating that each local control system is stable. However, since the spectral radius condition is only necessary, no definite conclusions can be made for the stability of the overall control system. An interesting result arising from this plot is that the spectral radius values are clustered around a value of 0.93 indicating that it may not be necessary to verify all the components for stability assessment. A similar result was recently observed by De Oliveira and Biegler (1994), but at the moment no systematic method exists to find the combination which will produce the upper bound on r. To test the conservativeness of the global nominal stability condition, an additional controller, FNWMPC-B was designed to satisfy the contraction mapping condition (Eq. (10)). The scaling matrix D is assumed to be the identity matrix. To satisfy Eq. (15), the control horizon was increased from p=15 to p=40, and the first 33 elements of Thot had to be set to zero, resulting in a very conservative weighing of the errors in Thot. The weighing matrix G was selected as follows: G= [0, . . . , 0 1, . . . , 1 1, . . . , 1 0, . . . , 0] ¿ËÀ ¿ËÀ ¿ËÀ ¿ËÀ 33 7 40 40

(21)

where the first 40 values in this matrix weigh the errors in Thot, the next 40 values weigh the errors in Tcold, and the remaining 40 weigh the errors in O2. Based on this selection of G, the first 33 errors in the horizon of Thot are not considered in the least squares solution and this will obviously result in inefficient control of the hot end temperature. The controllers designed based on the local and global stability criteria, FNWMPC-A and FNWMPC-B, respectively, are compared in Fig. 4 for a set point change of 10°C in Tcold. Clearly the controller designed based on the contraction mapping condition is more conservative resulting in poorer performance in Thot. It should be recognized that the selection of 0’s and 1’s in the weights composing the matrix G (Eq. (21)) may not be optimal. The search for an optimal set of weights which provides good performance and satisfies the contraction mapping condition is a difficult nonlinear optimization problem which was beyond the scope of this paper.

798

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

Fig. 3. Spectral radius plot for each component of the FNWMPC-A controller.

The robust stability of the FNWMPC-A controller is tested by computing m for each one of the 1296 components. The model uncertainty was obtained from several step tests performed on the system around different operating conditions. From the different step responses, uncertainties in the different transfer function parameters were identified. The uncertainty in the time constant was neglected since it was found to be insignificant (less than 10% variation in the whole window of operation). The uncertainties in the time delays and steady state gains are summarized in Table 1. From the parameter variations a nominal step or impulse model and the corresponding variation of the step or impulse model along the window of operation may be computed. The following observation was made from the identification of the process uncertainty: the entire step response curve corresponding to each transfer function shifts either up or down from the nominal response curve. This indicated that the individual variations in the step coefficients are changing in an approximately correlated fashion. Such correlation in the step coefficients variation is demonstrated in Fig. 5 where repeated step tests in the ID fan speed result in Thot response curve which lie entirely either above or below the nominal curve. Thus, the actual step response coefficients may be represented as follows: Si = S nominal +dS*d i i

d 51

(22)

The advantage of having correlated uncertainties in the step or the impulse response is that it is no longer necessary to consider individual perturbations for each

step or impulse response coefficient for the robust analysis calculations, which would result in a very conservative design. One single perturbation element may now be used for each one of the step response models, thus, avoiding overly conservative designs. For example, for the uncertainty description considered in Appendix B, if correlations are ignored the uncertainty block is as stated above D= [d1, d2, . . . , dn*nu* ny ]. When correlations are taken into account, the uncertainty matrix becomes D= [d *I 1 n, d *I 2 n, . . . , dnu* nyIn ]. Thus, for the lime kiln with two inputs and three output, only six scalar perturbations have to be considered when the observed correlations are accounted for. The m computation as a function of frequency for the FNWMPC-A controller is shown in Fig. 6. The maximum value in the curve was 0.357, indicating that robust stability of the local linear systems is ensured. Given that the m values are much smaller than one, it may be possible to design a less conservative controller.

5. Comparison of the finite number of weights model predictive controller to other model predictive control techniques To demonstrate the need for a soft constraint predictive control strategy, feasibility problems arising from the application of a standard QDMC, are illustrated in Fig. 7. These plots show the response of the variables to a step change of 80 ton day − 1 in the feed occurring at t= 0. The upper and lower constraint limits, in devia-

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

tion values from the initial steady state, were − 0.4 and 4% for O2 and − 20 and 20°C for Thot. The parameters for the QDMC controller are n = 40, p =15, Dt = 5, G=I and L= I. After 230 min, when the low constraint in Thot and O2 become active at the same time, the optimization problem to be solved by QDMC becomes unfeasible. For comparison purposes, a con-

799

troller based on the FE algorithm with one set of lower and upper constraints per variable was designed and included in Fig. 8. The tuning parameters n, p, Dt, G and L are the same as those selected for the QDMC controller and the weight which penalizes the constraint violation in oxygen is pOb = 50. In this case the feasibility problems are avoided but at the expense of a

Fig. 4. Comparison of the FNWMPC-A controller (solid) to the FNWMPC-B controller (dash) for a set point change of 10°C in Tcold.

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

800 Table 1 Uncertainty in the parameters y/u

Gain

Delay

Thot/ID Thot/NG Tcold/ID Tcold/NG O2/ID O2/NG

755K585 15K51.2 −125K5−15 0.25K50.4 −0.75K5−0.74 −0.0075K5−0.008

15tD53 135tD517 0.15tD50.3 15tD53

significant violation of the low limit on O2. The steady state value of O2. using the FE technique is −0.8%. The response of the FNWMPC-A controller is compared to the FE controller in Fig. 8. By using two levels of constraints in O2, the offset in this variable is considerably reduced as compared to the FE control. The final offset is − 0.52% in O2. The reduction of the oxygen constraint violation is achieved at the expense of a larger offset in the hot end temperature of approximately −20°C. However, this offset is still acceptable in view of the fact that the constraints in Thot are − 20°C. The constraint violation of the oxygen lower limit could be decreased for the FE controller by selecting a larger weight, e.g. pOb =200 which corresponds to the weight imposed in the second constraint level in FNWMPC (see Eq. (20)). However, this weight value is not optimal for disturbances smaller than 80 ton day − 1. For example, if the FNWMPC and FE controllers are simulated for a disturbance of 70 ton day − 1 the following steady state offsets are found: DO2 = −0.40% and DThot = −14.1°C for the FE controller as compared to DO2 = −0.45% and DThot = −6.4°C for the FNWMPC controller. The practically insignificant and almost unmeasurable improvement in DO2 obtained by the FE technique (0.05%) is achieved at the

Fig. 5. Example of correlated uncertainty for the step response between uID and Thot.

expense of a significant difference in the temperature offsets. If the oxygen weight will be further increased, the FE controller will respond for the smallest constraint violation, as a hard constraint technique such as QDMC. The advantages of constraint softening will then be lost. The FNWMPC-A controller was again compared to a standard QDMC controller for the same disturbance in the feedrate and the results are shown in Fig. 9. In this case, the constraints in Thot were removed from the QDMC formulation to avoid the feasibility problems mentioned above. The weight matrix G for the unconstrained situation is selected for both controllers equal to the identity matrix. It can be seen from the figure that although strict compliance with the oxygen constraint limit was achieved by using the QDMC, the response of Thot is significantly worst using QDMC as compared to the response obtained with FNWMPC-A. Thus, by softening the constraints, the FNWMPC controller provides more flexibility to balance the errors between the different variables according to priorities set by the operation. This is especially important for nonsquare system where due to insufficient degrees of freedom, steady state offsets cannot be completely removed. This last feature is one of the main advantages of the FNWMPC technique over conventional QDMC.

6. Conclusions A predictive control algorithm with soft constraints is proposed. The technique is based on the computation of proper weights in the objective function of the predictive control problem. The required computation time is found to be short and convergence to a global optimum is guaranteed. Moreover, since the number of weight combinations is finite and known a priori, the issues of nominal and robust stability of the control scheme can be assessed off-line. The proposed algorithm which is an extension of a technique previously proposed in the literature is studied for nominal stability and robustness. These issues were not studied before by the authors of the original algorithm on which our technique is based. Nominal stability is determined using either a sufficient condition based on contraction mapping concepts or by testing a sufficient instability condition based on a spectral radius test. Robust stability analysis is conducted using a m-structured singular value test. The FNWMPC controllers was designed and simulated for a lime kiln control problem. The nominal stability and robust stability tests were applied off-line for the FNWMPC controller. The sufficient nominal stability test was shown to be conservative and hence, it was concluded that in general applications, the spectral radius condition should be used. The existence of corre-

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

801

Fig. 6. m-Structured singular value plot for each component of the FNWMPC-A controller.

lation in the uncertainty for the kiln was taken into account to reduce the conservativeness of the design. Finally, FNWMPC was compared to the FE (Feher & Erickson, 1993) and QDMC algorithms. In both comparisons, FNWMPC proved to be more efficient in distributing the amount of error between the different output variables.

Appendix A

ÆS u1Ç ÃS u2à S u = à · Ã; ֈ ÈS unÉ ÆS u1,1,i S u1,2,i . . . . . . . . . . . S u1,nu,i Ç Ã Ã . S ui = à Ã, . à u à ÈS ny, 1, i . . . . . . . . . . . . . . . S uny,nu,iÉ

Development of the operator F Define the model update vector Y( (k) as follows: Y( (k) =MIY( (k− 1)+ S Du(k −1) +S Dd(k −1) u

d

(a.1) Y( (k) is a vector of length n. The matrix shift MI is given by:

Æ0 I ny 0 . . . . . . . . 0 Ç Ã0 0 I ny . . . . . . 0 à MI = à ...0 à à . I ny à È0 0. . . 0 I ny É u

and the step response matrix S is given by:

i= 1, . . . , n

and S d is defined as S u, but the step response coefficients of S d correspond to the effect of the measured disturbances on the outputs. Then define the p-step ahead prediction vector Y( (k+ 1/k): Y( (k+1/k)= MpY( (k)+S upDu(k)+ S dpDd(k) + W(k+ 1/k)

(a.2)

(a.3)

(a.4)

The matrix Mp is given by: Mp = [Ip·ny × p·ny 0]MI S up and S dp are the first p submatrices of S u and S d, respectively as defined in Eq. (a.3).

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

802

The vector W(k + 1/k) represents the unmeasured disturbances and error due to plant-model mismatch. For step like disturbances W(k + 1/k) is given by: W(k +1/k) = N2[y(k) −y¯ (k)]

(a.5)

where y(k) is the current ny output measured values, y¯ (k) are the first ny elements in the vector Y( (k) and: N2 =[Iny Iny. . . . .Iny ]T

(a.6)

y(k) and y¯ (k) are both of length ny and the matrix N2 is of dimensions ny by n *p. For the nominal stability y analysis we use W(k + 1/k) =0. Defining o(k + 1/k) as the feedback corrected vector of future output deviations from the set point trajectory, assuming all future moves are zero, as: o(k + 1/k)=R(k + 1) − MpY( (k) − S dDd(k) −W(k +1/k)

(a.7)

The prediction error is filtered by a diagonal filter Li which is included for the purpose of additional tuning. This type of filter which is a standard feature in some commercial predictive control algorithms (Richalet, Rault, Testud, & Papon, 1978) permits to shape the closed loop response of the system. The diagonal elements of Li are given by: ljj =

e − tdj s tpj s +1

j= 1, . . . , ny

(a.8)

The time constants of the filter were selected in our work as 100, 100 and 3 min for Thot, Tcold and O2, respectively. The filter also contains time delays which were selected approximately equal to the smallest delay associated to each one of the controlled variables. The delays are added in this way to avoid aggressive control actions. The delays used in our study were 5 min for Thot and Tcold and 0 for O2. Then, the least squares solution of the minimization problem with cost function defined by Eq. (1) is given by: Du(k)= KMPCLio(k+ 1/k)

(a.9)

where KMPC is given by: T T − 1 uT T u KMPC = (S uT Sp G G p G GS p + L L)

(a.10)

The objective of the development that follows is to obtain an expression for the operator F which relates U(k) to U(k−l), where U(k− 1) is a vector of length defined by the RHS of the following equation: n *n u x 1(k) D u(k− 1) x 2(k) D u(k− 2) · · · x n (k) D u(k− n)

(a.11)

Fig. 7. Feasibility difficulties encountered with a QDMC controller (solid) as compared to an unconstrained MPC controller (dash) for a step disturbance of 80 ton day − 1 in mud feedrate occurring at t =0. Constraints in both Thot and O2 are considered.

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

803

Fig. 8. Comparison of the FNWMPC-A controller (solid) to the FE controller (dash) for a step disturbance of 80 ton day − 1 in the mud feedrate.

Hence, U(k −1)= [u(k −1), u(k −2) . . . u(k− n)]T. Each u(k − i ) element in a.11 is a vector of length nu. In the absence of disturbances and assuming a perfect model and R=0 without loss of generality from Eqs. (a.7), (a.9) and (a.10), we obtain: U(k) = T2U(k −1)+T1KMPCLi [− MpY( (k)] T2 and T1 are given by:

(a.12)

ÆI nu ÃI nu T 2 = Ã0 Ã È0

0 0 0 0 I nu 0

. . . . I nu

. . . 0Ç . . . 0Ã . . . 0Ã . 0Ã 0É

ÆI nuÇ Ã0à T1 = à · à ÷à È0É

(a.13) T2 is of dimensions n *n by n *n and T1 is of dimensions u u

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

804

n *n by nu. Using the relation between step and impulse u response coefficients: j

Sj = % hi

j=1, . . . n

(a.14)

i=1

and assuming the system reaches steady state in n steps, the model update vector y(k) is given by:

Æ y(k) Ç Ã y(x+ 1) Ã Ã Ã · Y( (k)= Ã Ã · Ã Ã · Ã Ã · Èy(k+ n+ 1)É

Fig. 9. Comparison of the FNWMPC-A controller (solid) to a QDMC controller (dash) for a step disturbance of 80 ton day − 1 in mud feedrate. Only constraints in O2 are considered.

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

Æ u(k − 1) Ç Æs 1 h 2 . . . h nÇ Ã ÃÃ u(k − 2) Ã s h . . . hn 0 Ã2 3 ÃÃ · Ã . =Ã ÃÃ · ·Ã . Ã ÃÃ · Ã . Ã Ã Ã Ã · Ès n 0 . . . 0 É Èu(k −n) É = HU(k− 1)

where N1 has dimensions ny by n*ny. The uncertainty weights matrices W1 and W2, W1 = {dH1,1dH1,2 . . . dH1,nu dH2,1dH2,1 . . . dH2,nu . . . dHny,1dHny,2 . . . dHny,nu }T dH1, 1 = diag[dH1,1,i, . . . dH1,1,n ] (a.15)

dH1,2,i = [0 (a.16)

(a.17)

X(k) = F[X(k − 1)]

(a.18)

or equivalently from definition (a. 11). U(k) =f[U(k − 1)]

(a.19)

In order to apply the contraction mapping criteria we have to found F%, i.e. the derivative of F which is a matrix made of elements f %ij given by: f %ij =

(Fij (xi

(a.20)

In our case from Eqs. (a.17) and (a.19) F% is given by: F% = T2 −T1KMPCLi MpH

(a.21)

The notation F%Ji ,nom is used for the nominal case component Ji. For the nominal case, the nominal hi coefficients are used to compute the matrix H. The notation F%Ji ,true is a function of the actual hi ’s and is used for robustness analysis.

0

0

...

dh1,2,i

0

0

0] ...

(b.4) 0]

Æ B O ··· O Ç Ã Ã O B ··· O W2 = à à ·· — · à à ÈO ··· B É

1

···

1]

N1 = [Iny

0 0 . . . 0]

(b.1)

(b.7)

U(k)= M11U(k− 1)+M12w(k)

(b.8)

z(k)= M21U(k− 1)+M22w(k)

(b.9)

These state space equations are schematically represented by the M-delta interconnection structure given in Fig. 2. This structure is known in the robust control literature as linear fractional transformation or LFT. By comparison of Fig. 2 and Fig. 1 the elements of the matrix M are given as follows: (b.10)

M11 = F%Ji and F%Ji is given in Appendix A. M12 = − T1KMPCLi N2W2

(b.11)

M21 = W1

(b.12)

M22 = [0]

(b.13)

Appendix C

The closed loop system with the model predictive controller used in this work may be represented by the block structure shown in Fig. 1. The blocks Fi, KMPC, MI, N2. T1, T2, H and Mp are defined in Appendix A. The matrix N1,

(b.6)

dHi,j,k is 1 by nu, W2 is ny by (ny · nu · n) and B is 1 by (nu · n). The block diagram in Fig. 1 may be represented by the following state space equations:

Appendix B Development of the M matrix

(b.5)

where: B= [1

The operator F is defined in our work as:

(b.3)

and dhi,j,k is the uncertainty or maximal variation of the impulse response coefficient hi,j,k over the window of operating conditions and for example, dH1,1,i = [dh1,1,i

Then using Eqs. (a.12) and (a.15) we obtain U(k) = (T2 −T1KMPCLi MpH)U(k − 1)

(b.2)

where,

where,

Æ h 1,1,i . . . . h1,nu,i Ç Ã Ã · Ã Ã · hi = Ã Ã · Ã Ã · Ã Ã Èh ny,1,i . . . hny,nu,i É

805

Notation a D e f %ij

deviation from constraint scaling matrix (Eq. (15)) unit vectors (Eq. (10)) coefficient in the operator F (Eq. (a.15))

806

F G hi H Ji KMPC Li ljj M MI Mp n nu ny no ns N1 N2 p pij P R Si Su Sd T T1 T2 u U V W W1 W2 x X y

R. Zano6ello, H. Budman / Computers and Chemical Engineering 23 (1999) 791–806

operator (Appendix A) transfer function impulse response coefficient impulse response matrix (Eq. (a.15)) component (specific combination of p weights) model predictive controller gain matrix diagonal filter operating on the prediction errors elements of the filter Li matrix defined in Fig. 1 and Eqs. (b.10) and (b.11) shift matrix truncated shift matrix prediction horizon number of inputs number of outputs n *n y n *n (total number of states) u Eq. (b.1) Eq. (a.6) control horizon inputs or outputs weights in the soft constrained problem penalty function for the soft constrained problem set point vector step response coefficients step response model (manipulated variables to outputs) step response model (disturbances to outputs) temperature Eq. (a.13) Eq. (a.13) manipulated variable vector of manipulated variables objective function of the unconstrained problem unmeasured disturbances vector weighting matrix (Eq. (b.2)) for uncertainty block weighting matrix (Eq. (b.6)) for uncertainty block state vector of x states output

.

Y Yc Greek letters a b s m d D o r L G

output prediction vector constraints vector constraints on inputs (Eq. (4)) constraints on outputs (Eq. (4)) singular value structured singular value scalar uncertainty element uncertainty block (Eq. (17)) feedback corrected error vector spectral radius input weights matrix for the unconstrained problem output weights matrix for the unconstrained problem

References Charos, G. N., Taylor, R. A., & Arkun, Y. (1991). Model predictive control of an industrial lime kiln. Tappi Journal, 74 (2), 203. De Oliveira, N. M., & Biegler, L. T. (1994). Constraint handling and stability properties of model predictive control. American Institute of Chemical Engineers Journal, 40 (7), 1138. Doyle, J., & Packard, A. (1988). Robust control of multi6ariable and large scale systems. Washington D.C.: USAF-Office of Scientific Research. Feher, J. D., & Erickson, K. T. (1993). Solving the model predictive control problem with soft constraints, Proceedings of the American Control Conference (pp. 377 – 378), San Francisco, CA. Fletcher, R. (1987). Practical methods of optimization (2nd ed.). New York: Wiley. Kouvaritakis, B., Rossiter, J. A., & Chang, A. O. T (1992). Stable generalised predictive control: an algorithm with guaranteed stability. IEE Proceedings D Control Theory and Applications, 139 (4), 349. Morari, M. (1989). Industrial course notes on model predicti6e control. Pasadena, CA: California Institute of Technology. Rossiter, J., & Kouvaritakis, B. (1993). Constrained stable generalised predictive control. IEE Proceedings D Control Theory and Applications, 140 (4), 243. Richalet, J., Rault, A., Testud, J. L., & Papon, J. (1978). Model Predictive Heuristic Control: Application to Industrial Processes. Automatica, 14, 413. Ricker, N., Subrahmanian, T., & Sim, T. (1989). Case studies of model predictive control in pulp and paper production. In Model based control (pp. 13 – 22). Oxford: Pergamon Press. Webb, C. J., Budman, H. M., & Morari, M. (1995). Identification of uncertainty bounds for robust control with applications to a fixed bed reactor. Industrial Engineering and Chemical Research, 34, 1743. Zafiriou, E. (1990). Robust model predictive control of processes with hard constraints. Computers and Chemical Engineering, 14 (5), 359.