Ind. Eng. Chem. Res. 2004, 43, 7807-7814
7807
On the Tuning of Predictive Controllers: Inverse Optimality and the Minimum Variance Covariance Constrained Control Problem Donald J. Chmielewski* and Amit M. Manthanwar Department of Chemical and Environmental Engineering, Illinois Institute of Technology, Chicago, Illinois 60616
This paper presents a systematic tuning approach for linear model predictive controllers based on the computationally attractive minimum variance covariance constrained control (MVC3) problem. Unfortunately, the linear feedback policy generated by the MVC3 problem is incompatible with the algorithmic framework of predictive control, in which the primary tuning vehicle is the selection of objective function weights. The main result of this paper is to show that all linear feedbacks generated by the MVC3 problem exhibit the property of inverse optimality with respect to an appropriately defined linear quadratic regulator (LQR) problem. Thus, the proposed tuning scheme is a two-step procedure: application of the MVC3 problem to achieve tuning objectives, followed by application of inverse optimality to determine the predictive control weights from the MVC3-generated linear feedback. 1. Introduction The advent of model predictive control (MPC) has revolutionized the chemical process control community. Its ability to incorporate meaningful limits on manipulated as well as control variables has allowed the industry to move away from traditional regulation-type issues and focus on the economics of operating point selection. Traditional designs were typically guided by a fear of constraints, resulting in the “bird house” effect (i.e., the operating point and its dynamic region being placed at the center of a constraint polytope; see Figure 1). As constrained MPC has reduced this fear of constraints, the industry has moved toward the more profitable regions along the edges of the constraint polytope, the most profitable being termed the economically optimal operating point. Although operation at this optimum is desired, it typically cannot be achieved with simultaneous satisfaction of all constraints, because of the influence of external disturbances. Loeblein and Perkins1 recognized this issue and advocated the selection of a backed-off operating point. The amount of backoff (which is proportional to the loss of profit) is governed by the size and shape of the dynamic operating region. This interplay between the dynamic operating region and the set of process constraints forms the basis of our automated tuning scheme. Although the literature on controller tuning is vast, only a handful of efforts have addressed the tuning of predictive controllers. The early works focused on sensitivity studies of the closed-loop system with respect to the standard dynamic matrix control tuning parameters of horizon size and move suppression (see, for example, Marchetti et al.2). Tuning strategies based on the methods of principle component analysis were proposed by both Maurath et al.3 and Challaghan and Lee.4 A number of works have focused exclusively on stability and robustness (for stability, see Clark and Scattolini5 or Rawlings and Muske;6 for robustness, see Ohshima et al.7 or Lee and Yu8). In an effort to simplify the tuning issue, Shridhar and Cooper9 presented a * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: 312-567-3537. Fax: 312-567-8874.
Figure 1. Dynamic economics of process control.
tuning strategy based on first-order plus dead time approximations of the open-loop system. Al-Ghazzawi et al.10 proposed an on-line adaptation strategy of MPC parameters for step response models using the timedomain response envelopes as the performance specification. Each of the above methods focuses on obtaining a predictive controller with “well-behaved trajectories” in the face of set-point changes and step-like disturbance inputs. Although these inputs are important, the typical time profile of a chemical process is such that shortduration step changes are followed by long periods in which the only excitation is due to relatively small stochastic inputs. On the basis of this view of chemical processes, Loeblein and Perkins1 introduced the dynamic operating region and operating point back-off concepts. Their controller design method was based on an iterative search scheme in which the dynamic region was calculated manually at each step. van Hessem et al.11 also approached the tuning problem from this perspective. Their design formulation focused on minimizing the standard deviation distance between the backed-off and economically optimal operating points. Their approach was to identify near-optimal linear feedbacks and then apply an unspecified inverse-optimal synthesis scheme to determine the objective function weights. Chmielewski12 proposed a tuning method that did not require an inverse-optimal synthesis scheme.
10.1021/ie030686e CCC: $27.50 © 2004 American Chemical Society Published on Web 10/16/2004
7808
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004
Unfortunately, this method did not result in a tight tuning algorithm (i.e., the computational conditions were sufficient but not necessary). Regarding inverse optimality, Kalman13 introduced the subject by posing the question: For a given feedback gain, do there exist corresponding tuning values of state and control weight functions? Other efforts pertaining to the LQR inverse optimality question include those of Anderson and Moore,14 Yokoyama and Kinnen,15 and Park and Lee.16 Of particular note is the result found in Boyd et al.,17 which employs linear matrix inequality (LMI) methods to answer the above question. In the next section, we describe the tuning framework and provide motivation for the results to follow. In section 3, we define the MVC3 problem and highlight its computational features. In section 4, we present the main results of the paper. In particular, it is shown that all MVC3 problem solutions exhibit the property of being LQR inverse optimal. In sections 5 and 6, we discuss inverse-optimal synthesis as well as extensions to the discrete-time and partial-state-information cases. Finally, section 7 illustrates the proposed tuning scheme with a couple of examples. 2. Tuning Framework Consider a linear feedback system with outputs z(t) and manipulated variables u(t). Additionally assume that the system will be excited by a class of disturbance inputs w(t). Given some characterization of these disturbances, one can easily calculate the expected responses of z(t) and u(t). This is frequently done in PID tuning, where step input disturbances of a particular magnitude are expected and tuning parameters are selected so as to maintain both z(t) and u(t) within some bounds. Now consider a linear time-invariant statespace system, x˘ (t) ) Ax(t) + Bu(t) + Gw(t), along with a set of pointwise-in-time constraints, |zi(t)| < zji for all t g 0, where zi is the ith element of the vector z(t) ) Cxx(t) + Cuu(t). If we then apply a quadratic infinitetime objective function, ∫∞0 (xTQx + 2xTMu + uTRu) dt, the resulting tuning question is as follows: Determine a triple (Q, R, M) such that the unconstrained optimal controller will satisfy the imposed signal bounds in the face of expected disturbances. Clearly, this is a rather general statement of the tuning activity and will strongly depend on one’s interpretation of the phrases “satisfy signal bounds” and “expected disturbances”. In the current article, we assume the class of disturbances w(t) to be stationary, Gaussian, white-noise processes with zero mean and covariance Σw. Additionally, “satisfaction of the signal bounds” is defined as the variance of the signal zi(t) satisfying the constraints σzi2 < zji2. The manual approach to our tuning charge is given as follows:1 Assume that we have selected a triple Q, R, and M. From this triple, calculate the unconstrained linear feedback u(t) ) Lx(t), where L ) -R-1(MT + BTP), where P is the positive semidefinite solution to the algebraic Riccati equation ATP + PA + Q - (PB + M)R-1(MT + BTP) ) 0. Next calculate Σx g 0 from the Lyapunov equation (A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0. Finally, the steady-state covariance matrix of the signal z(t) is calculated as Σz ) (Cx + CuL)Σx(Cx + CuL)T. From the diagonal elements of Σz, we can now pick off the individual variance values of the vector z(t), σzi2 ) φiΣzφiT, where φi is the ith row of an appropriately dimensioned identity matrix. Thus, the final step in the manual tuning algorithm is to check whether the
Figure 2. Covariance ellipsoids for tuning.
constraints σzi2 < zji2 are all satisfied. If any of the constraints are violated, then we must reselect the Q, R, and M weights and recalculate the σzi2 terms. A more graphical interpretation of the relationship between the covariance Σz and the parameters zji is given as follows: Define the dynamic operating region as the elliptical set Z ) {ζ ∈ Rp|ζTΣz-1ζ e 1}, which represents the region in which the signal z(t) will be found 63% of the time. If we now contrast this set with the bounds zji, defining a box in the Rp space, then we find a situation similar to that of Figure 2. Clearly, a Z region similar to that in Figure 2a represents a poor tuning, whereas Figure 2b represents acceptable tuning. Thus, the graphical interpretation of our tuning challenge is to select a triple Q, R, and M such that the resulting dynamic operating region Z will be contained in the box set defined by the bounds zji. In the context of constrained infinite-time linear quadratic optimal control (CITLQOC) as defined in refs 18-20, we find the above notion of tuning to be of relevance. Insofar as the constrained controller will be similar to its corresponding unconstrained counterpart, it will be advantageous to have tuning parameters that are cognizant of the constraints to be imposed. That is, given the expected disturbances, we desire an unconstrained closed-loop system predisposed to the observation of constraints. If we then impose hard constraints on this already amiable controller, then we will find it less likely to encounter infeasible conditions (as discussed in ref 21-23). 3. Minimum Variance Covariance Constrained Control The MVC3 problem is defined as p
min {
∑diζi}
(1)
Σxg0,L,ζi i)1
s.t.
(A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0 φi(Cx + CuL)Σx(Cx + CuL)TφiT ) ζi, i ) 1, ..., p ζi < zji2, i ) 1, ..., p
where ζi represents the variances σzi2 of signal zi. The goal of the MVC3 problem is to minimize the weighted, di, sum of the output variances ζi. The constraints ζi < zji2 ensure that the dynamic operating region is within the constraint box. It is important to note that the central optimization variable of this formulation is the linear feedback gain, L, which generates the covariance matrix, Σz, defining the dynamic operating region.
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004 7809
The MVC3 problem was originally introduced as a linear feedback existence problem (see Skelton et al.24 for details). There, the question posed was: Does there exist a linear feedback L such that the steady-state covariance matrix Σz satisfies a set of bounds. A key result of ref 24 was to show that this nonlinear existence question has a convex form owing to the technology of LMIs. Here, we present a slight generalization of the results found in ref 25 (see theorem 6.15), the proof of which can be found in the Appendix. Theorem 3.1. ∃ stabilizing L and Σx g 0 s.t. ζi, i ) 1, ..., p
(A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0 T
T
φi(Cx + CuL)Σx(Cx + CuL) φi ) ζi, i ) 1, ..., p ζi < zji2, i ) 1, ..., p
[
ζ′i T
T
(CxX + CuY) φi
Consider a dynamic system x˘ ) Ax + Bu along with an output signal z ) Cxx + Cuu. The classic LQR-optimal control problem is stated as
∫0∞zTDz dt}
s.t.
x˘ ) Ax + Bu(x), x(0) ) x0 z ) Cxx + Cuu(x)
(8)
∫TE[zTDz] dt}
(10)
x˘ ) Ax + Bu(x) + Gw
Again, the optimal feedback function is given by u(x) ) -R-1(MT + BTP)x(t), where P, Q, R, and M are defined as above. Our reformulation of the LQR problem will start by exploiting the seminal certainty equivalence results of Kalman and Koepcke,29 so as to limit the search space to that of linear feedbacks (i.e., u ) Lx). The next step is to observe that the objective function of problem 10 can be converted as follows
lim
4. Reformulation of the LQR Problem
1
Tf∞T 0
z ) Cxx + Cuu(x)
(7)
9 Application of theorem 3.1 to problem 1 will result in an exact conversion to a convex form (as LMIs are widely known to be convex constraints; see ref 26 or 27). The desired linear feedback is obtained as L* ) Y*X*-1, where Y* and X* denote the optimal solution matrices of the new LMI-constrained minimization problem. It is interesting to note that conditions 5-7 are exactly those required to determine the feasibility of problem 1. If the problem turns out to be infeasible, then the only recourse is to enlarge the bounding region defined by the zji terms. Developing economic-based methods of reselecting bounds is an important aspect of future research. Although the MVC3 problem appears to meet the tuning challenge, there is one very important disconnect. That is, MPC tuning requires the specification of objective function weights (i.e., Q, R, and M) because of its on-line optimization-based implementation. Thus, the disconnect is that we have no guarantee that there will exist LQR weights corresponding to the MVC3 linear feedback. This brings us to the primary objective of the paper. In the next section, we show that the linear feedback solution generated by an MVC3 problem exhibits the property of being inverse optimal.
u(‚)
s.t.
(4)
(5)
{
min lim
(3)
]
min{
Now, consider a stochastic process x˘ ) Ax + Bu + Gw along with an output signal z ) Cxx + Cuu. The classic linear quadratic Gaussian (LQG) optimal control problem can be stated as
u(‚)
φi(CxX + CuY) > 0, i ) 1, ..., p X (6)
ζ′i < zji2, i ) 1, ..., p
Q ) CxTDCx, R ) CuTDCu, and M ) CxTDCu (9)
(2)
if and only if ∃ Y, X > 0 and ζ′i, i ) 1, ..., p s.t.
(AX + BY) + (AX + BY)T + GΣwGT < 0
p diφiTφi, di g 0. The optimal feedback where D ) ∑i)1 28 function is given by u(x) ) -R-1(MT + BTP)x(t), where P is the positive definite solution to the algebraic Riccati equation ATP + PA + Q - (PB + M)R-1(MT + BTP) ) 0 and
Tf∞
1
∞ E[zTDz] dt ) limE[zTDz] ) ∫ 0 tf∞ T p
Tr{DΣz} )
diφiΣzφiT ∑ i)1
(11)
Now consider the following minimum variance control (MVC) problem p
min {
∑diζi}
(12)
Σxg0,L,ζi i)1
s.t.
(A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0 φi(Cx + CuL)Σx(Cx + CuL)TφiT ) ζi, i ) 1, ..., p
along with the following hypotheses: Assumption 1. R > 0 and Q - MR-1MT g 0. Assumption 2. The pairs (A - BR-1MT, B) and (A BR-1MT, Q) are stabilizable and detectable, respectively. Assumption 3. The pair (A, GΣwGT) is controllable. Theorem 4.1. Let hypotheses 1-3 hold; then the solution to MVC problem 12 is coincident with the solution to LQR problem 8. Proof. Hypotheses 1 and 2 indicate that the solutions to problems 8 and 10 are unique and stabilizing. Hypotheses 1-3 ensure that the solution to problem 12 is unique and stabilizing. Next we observe that, from our construction of problem 12, it is equivalent to problem 10), in the sense that the linear feedback generated by 12 is the solution to problem 10. Finally, certainty equivalence between problems 10 and 8 completes the proof. 9 Now let us return to the MVC3 problem. The difficulty with this problem is that its equivalence with the LQR problem appears to be disrupted by the constraints ζi < zji2. However, the following theorem indicates that such a disruption will not preclude the existence of
7810
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004
corresponding LQR weights, but will only serve to modify them. Theorem 4.2. Let hypotheses 1-3 hold; then the solution, L*, to problem 1 is coincident with the solution to an appropriately defined LQR problem. Proof. Building on the convexity results of theorem 3.1 and the duality theorem of ref 30, problem 1 can be exactly restated (i.e., no duality gap exists) as
max λig0
{
p
min L,Σxg0,ζi
s.t.
{
∑
p
d iζ i +
i)1
∑λ (ζ - zj i
i
2
i
)}
i)1
(A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0 φi(Cx + CuL)Σx(Cx + CuL)TφiT ) ζi
}
(13)
If we rewrite the minimization objective function as p
p
(di + λi)ζi} - ∑λizji2 ∑ i)1 i)1
min {
L,Σxg0,ζi
then it is clear that the assumptions of theorem 4.1 are satisfied for all values of λi g 0. This indicates that all λi-dependent solutions, L*(λi), coincide with the solution to some LQR problem (d h i ) di + λi). Thus, L* ) L*(λ/i ) / (where λi is the solution to the max problem) is the solution to some LQR problem. 9 5. Inverse-Optimal Control Theorem 4.2 guarantees that the MVC3 problem will generate a linear feedback such that there exists an LQR problem for which it is a solution. Unfortunately, the exact form of this LQR problem is unclear, unless the MVC3 solution procedure provides us with the optimal Lagrangians λ/i . In the absence of such parameters, we are left with an LQR-optimal feedback, L*, and charged with the task of determining the objective function weights (Q, R, and M) such that the Riccati equation will regenerate L*. Toward the development of a computationally attractive synthesis procedure, we present the following theorem. Theorem 5.1. If ∃ P > 0 and R > 0 s.t.
{
T LTRL - ATP - PA -(L R + PB) T R -(RL + B P)
}
>0
then Q ) ˆ LTRL - ATP - PA and M ) ˆ - (LTR + PB) will be s.t.
[
]
Q M >0 MT R
and P and L will satisfy
0 ) ATP + PA + Q - (M + BP)R-1(MT + PBT) L ) -R-1(MT + PBT) Proof. The result is easily verified by simple substitutions. 9 It is noted that this theorem is not tight in the sense that, if L is the result of an LQR problem in which Q is only positive semidefinite but (A, Q) is detectable, then the above procedure might not find acceptable weights.
However, in the present context, we can avoid this possibility by ensuring that CxTDCx > 0 in the MVC3 problem. At this point, it is instructive to review the proposed method of implementation. The envisioned scheme is to allow the economic LP (or RTO) to select the steadystate operating point with sufficient back-off from the constraints. This selection of operating point is usually infrequent, but in some implementations, it can occur at every time step. From this operating point and its relative closeness to the constraints, one has enough information to run the MVC3 algorithm and subsequent inverse-optimality scheme. The result is the generation of the triple (Q, R, M) along with the solution, P, to the appropriate Riccati equation. These four matrices can now be used in the objective function of a numerically solved, constrained infinite-horizon controller,19,20 with the implicit assumption that the MVC3-generated feedback, L*, will be used as the controller in the terminal region. As outlined above, these weights will be such that the unconstrained controller and its constraints are in harmony. If the steady-state operating point and the constraints do not change between sample periods, then there will be no motivation to recalculate the weights. However, if the operating point does move with respect to the constraints, then a recalculation of the weights will be desired. In fact, when such movement occurs, it will be desirable to obtain new weights as soon as possible, preferably within one sample period. Given the convexity of our formulation, it is reasonable to assume that such objectives can be met. An additional note concerning assumption 3 is also in order. In the absence of assumption 3, a subspace of the state will not be excited by the disturbance w in the open loop. If problem 1 finds no advantage to activating this subspace through the feedback u, then elements of L might have no impact on the minimized objective function of 1. Unfortunately, this nonuniqueness with respect to L is not an option for problem 8, which is clearly unaffected by GΣwGT. Thus, without assumption 3 it is possible for problem 1 to generate a non-LQR feedback, which will further cause the inverse-optimal synthesis routine to fail. If this is the case, then it will manifest through the singularity of Σx. If singularity is found, then it is suggested that a state-space reduction be performed that will remove the unexcited states of the closed-loop system. An alternative to state-space reduction is to artificially add a small amount of disturbance to every state. If sufficiently small, this additional excitation will have a negligible impact on the variance magnitudes but will guarantee the uniqueness of L. An inspection of the proof of theorem 3.1 indicates that the LMI version of problem 1 implicitly adds an amount of disturbance to all states. Thus, an additional advantage of the LMI version of problem 1 is the ability to remove assumption 3 from theorems 4.1 and 4.2. 6. Extension to the Discrete-Time and PSI Cases In most MPC implementations, the process model is in a discrete-time form, and the system state is usually generated via a Kalman filter augmentation. In this section, we indicate how the above results can be extended to these scenarios. Conceptually, the discretetime full-state-information case is identical to the above, the only difference being the specific form of the equations involved. As an aid to the reader, we present the following discrete-time versions of theorems 3.1 and 5.1.
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004 7811
Theorem 6.1. ∃ stabilizing L and Σx g 0 s.t. ζi, i ) 1, ..., p
p
min { Σxg0,L,ζi
(A + BL)Σx(A + BL)T + GΣwGT ) Σx
(14)
φi(Cx + CuL)Σx(Cx + CuL)TφiT ) ζi; i ) 1, ..., p (15) ζi < zji2; i ) 1, ..., p
(16)
diζi} ∑ i)1
(20)
s.t. (A + BL)Σxˆ + Σxˆ (A + BL)T + ΣeCTΣv-1CΣe ) 0 φi(Cx + CuL)Σxˆ (Cx + CuL)TφiT + φiCxΣeCxT φiT ) ζi, i ) 1, ..., p ζi < zji2, i ) 1, ..., p
if and only if ∃ Y, X > 0 and ζ′i, i ) 1, ..., p s.t.
[
]
X - GΣwGT (AX + BY) >0 (AX + BY)T X
[
ζ′i T
T
(CxX + CuY) φi
]
(17)
φi(CxX + CuY) > 0, i ) 1, ..., p X (18)
ζ′i < zji2, i ) 1, ..., p
(19)
[
7. Design Examples Example 1. Consider an underdamped mass spring damper system described by
9
Theorem 6.2. If ∃ P > 0 and R > 0 s.t. P - ATPA + LT(R + BTPB)L -LT(R + BTPB) - ATPB R -(R + BTPB)L - BTPA
where Σv and Σe are known matrix constants (from the steady-state Kalman filter solution). Once the linear feedback, L*, has been generated, the inverse-optimality procedure is the same as before. The discretetime version of the PSI results are given in the Appendix.
A)
]
>0
]
Q M >0 MT R
Cx )
L ) -(R + BTPB)-1(M + ATPB)T 9 Next, we consider the partial-state-information (PSI) case. Consider the system x˘ ) Ax + Bu + Gw, along with the measurement equation y ) Cx + v, where v(t) is a zero-mean, Gaussian white-noise process with covariance matrix Σv. In the PSI case, the generation of a feedback signal is done with respect to the state estimate u ) Lxˆ , because the actual state, x, is unavailable. If the state estimate is generated by the Kalman ˘ where K ) ΣeCTΣv-1 and filter xˆ˘ ) xˆ˘ + Bu + K(y - Cxˆ), T Σe > 0 satisfies AΣe + ΣeA + GΣwGT + ΣeCTΣv-1CΣe ) 0, then the innovations signal, defined as η ) y - Cxˆ , is a zero-mean, Gaussian white-noise process with covariance Ση ) Σv.31 Additionally, the orthogonality between the Kalman state estimate and its error signal (i.e., E[(x - xˆ ) xˆ T)] ) 0) leads to the identity Σx ) Σxˆ + Σe. As a result, controller design can be performed with respect to the estimator system xˆ˘ ) (A + BL)xˆ + ΣeCTΣv-1η, while still keeping track of the true state covariance Σx. Thus, the partial information MVC3 problem is stated as
]
B)
[] 0 1
G)
[] 0 1
(21)
[
-3 -2 0 0
]
Cu )
[] 1 1
If we further define d1 ) 1 and d2 ) 0.01, then the resulting LQR weights are calculated via eq 8 as
and P and L will satisfy
P ) ATPA + Q (M + ATPB)(R + BTPB)-1(M + ATPB)T
0 1 -3 -2
The disturbance force is assumed to be of magnitude Σw ) 0.12. Next, define the output, z(t), by the matrices
then Q ) ˆ P - ATPA + LT(R + BTPB)L and M ) ˆ - LT(R + BTPB) - ATPB will be s.t.
[
[
Q)
[ ]
9 6 , 6 4
R ) 1.01,
MT )
[ ] -3 -2
which results in the generation of the unconstrained linear feedback L* ) [2.7015 1.2402]. The resulting output closed-loop covariance matrix is
Σz )
[
0.0058 -0.0239 -0.0239 0.1677
]
which is translated into a dynamic operating region indicated by the dotted line of Figure 3 (denoted by zj2 ) ∞). Next we consider the impact of physical limits on the actuator device. These will manifest as a variance constraint on the second output signal (i.e., ζ2 < zj22). We can interpret this constraint as a reflection of the quality of the purchased actuator. For the cases of zj2 ) 0.2 and zj2 ) 0.05, the resulting dynamic operating regions are depicted in Figure 3. The final step in our tuning scheme is to convert the linear feedback gains of each case into the appropriate objective function weights. These are presented in Table 1. It is highlighted that the objective function weights of the unconstrained case (i.e., zj2 ) ∞) differ from those above. This illustrates the nonuniqueness of the inverseoptimal synthesis scheme.
7812
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004
Figure 3. Covariance ellipsoids for example 1.
Figure 5. Reactor temperature vs furnace temperature.
Figure 4. Preheating furnace reactor system. Table 1. Tuning Parameters for Example 1 zj2
L*
∞
L ) [2.7015 1.2402 ]
0.2
L ) [2.3263 0.8833 ]
0.05 L ) [1.1779 0.2899 ]
( ( (
Q
1.0600 0.0591 0.0591 1.0641 1.1118 0.0549 0.0549 1.0933 2.0183 0.1016 0.1016 1.8847
) ) )
R 0.1079 0.1519 0.8193
( ( (
M 0.2454 -0.1139 0.3050 -0.1337 0.8182 -0.2899
) ) )
Figure 6. Reactor temperature vs O2 concentration.
Example 2. Consider the preheating furnace reactor system of Figure 4, which is described by the following matrices
[
]
-8 0 0 0 2 -1.5 0 0 A) × 103 0 0 -5 0 0 0 0 -5
[
-75 75000 0 -25 0 0 B) 0 -8500 8.5 × 105 0 0 -5 × 107 10000 0 G) 0 0
[ ]
]
(22)
Figure 7. Reactor temperature vs CO concentration.
The disturbance size is assumed to be Σw ) 0.27952. The open-loop dynamic operating region is given in Figure 5 (dotted line). Next define Cx and Cu such that [Cx l Cu] equals the identity matrix. Then application of eq 9 will result in Q ) I4×4, R ) I3×3, and M ) 0. The LQR dynamic operating region is depicted in Figures 5-10
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004 7813
[
Q) 494.5008 508.3485 1184.7015 28.1097 508.3485 7959.7980 -1634.2245 15.5088 1184.7015 -1634.2245 27084.2987 472.0720 28.1097 15.5088 472.0720 8.7327
[
495.7674 -270.4322 37.9161 R ) -270.4322 950.6331 150.6843 37.9161 150.6843 6444.2294
[
]
-307.2749 -100.0023 -77.6836 -1663.6668 1413.6481 -16.1815 M) -574.5418 -690.2072 -304.6972 -22.5859 -4.5071 -4.0311
]
]
8. Conclusion
Figure 8. Reactor temperature vs reactant feed rate.
In this work, we present an automated tuning scheme for linear model predictive controllers excited by stochastic disturbances. It is proposed that such a scheme will serve to create an alignment between the controller weights (i.e., tuning parameters) and the system constraints. It should be noted that the stochastic framework of the method will require an additional level of modeling before implementation. However, this type of disturbance modeling is beginning to gain support within the chemical process control community and might become common practice in the coming years. Computationally, the scheme is extremely efficient and robust (owing to the convexity of the formulations) and will likely be a very small fraction of the system sampling interval. Additionally, the proposed scheme is easily extended to the discrete-time and partial-stateinformation cases. Acknowledgment
Figure 9. Reactor temperature vs fuel feed rate.
This work was supported by the Department of Chemical and Environmental Engineering and the Armour College of Engineering and Science, both at the Illinois Institute of Technology. Notation
Figure 10. Reactor temperature vs vent position.
(solid lines). Application of the MVC3 problem in which the box constraints of Figures 5-10 are imposed results in the MVC3 dynamic operating region of Figures 5-10 (dashed lines). Finally, application of the inverseoptimality synthesis problem results in the weighting matrices
x, xi ) state variables u, ui ) manipulated variables w, wi ) disturbance input z, zi ) output variables v, vi ) measurement noise A, B, C, G ) system matrices Cx, Cu ) output matrices Q, R, M ) weighting matrices P, L ) Riccati solution and resulting gain Σe, K ) Kalman filter solution and resulting gain Σx, Σu, Σv, Σw ) covariance matrices σxi2, σui2, σzi2 ) variances of signals xi, ui, and zi, respectively zji ) bounds on the ith output Z ) elliptical dynamic operating region ζ, ζi ) dummy variables for the variances of z and zi, respectively di, D ) diag(di) ) cost associated with ith variable ) positive constant X, Y ) LMI variables φi ) ith row of an identity matrix λi ) Lagrange multiplier
Appendix Proof of Theorem 3.1. “Only If” Part. By assumption, (A + BL) is stable. Thus
7814
X)
Ind. Eng. Chem. Res., Vol. 43, No. 24, 2004
∫0∞e(A+BL)t(GΣwGT + 2I)e(A+BL) t dt ) T
Σx + 2E > 0 (A-1) will exist and satisfy (A + BL)X + X(A + BL)T + GΣwGT + 2I ) 0. If we define Y ) LX, then (AX + BY) + (AX + BY)T + GΣwGT ) -2I < 0 will result. From X ) Σx + 2E, one easily finds φi(CxX + CuY)X-1(CxX + CuY)TφiT < zji2 + φi(Cx + CuL)2E(Cx + CuL)TφiT. Then, for small enough 2, one finds φi(CxX + CuY)X-1(CxX + CuY)TφiT < zji2, because E > 0. Finally, introduction of a variable ζ′i between the strict inequality will yield eqs 6 and 7. “If” Part. By assumption, ∃ X > 0, Y s.t. (AX + BY) + (AX + BY)T + GΣwGT < 0. Let L ) YX-1, then (A + BL)X + X(A + BL)T < 0 w (A + BL) is stable. Define
Σx )
∫0∞e(A+BL)tGΣwGTe(A+BL)Tt dt
Then, Σx g 0 satisfies (A + BL)Σx + Σx(A + BL)T + GΣwGT ) 0. Subtracting this from the above inequality, we find (A + BL)(X - Σx) + (X - Σx)(A + BL)T < 0. This along with (A + BL) being stable indicates that X - Σx > 0 w X > Σx. From this point, the remaining inequalities are trivially observed. 9 Discrete-Time PSI Version of Theorem 6.1. Theorem A.1. ∃ stabilizing L and Σx g 0 s.t. ζi, i ) 1, ..., p
(A + BL)Σxˆ (A + BL)T + AΣeCT(CΣeCT + Σv)-1CΣeAT ) Σxˆ φi(Cx + CuL)Σxˆ (Cx + CuL)TφiT + φiCxΣeCxT φiT ) ζi, i ) 1, ..., p ζi < zji2, i ) 1, ..., p if and only if ∃ Y, X > 0 and ζ′i, i ) 1, ..., p s.t.
(
(
)
X - AΣeCT(CΣeCT + Σv)-1CΣeAT (AX + BY) >0 X (AX + BY)T
ζ′i - φiCxΣeCxT φiT φi(CxX + CuY) (CxX + CuY)TφiT
X
)
> 0, i ) 1, ..., p
ζ′i < zji2, i ) 1, ..., p 9 Literature Cited (1) Loeblein, C.; Perkins, J. D. Stuctural design for on-line process optimization: I. Dynamic economics of MPC. AIChE J. 1999, 45, 1018. (2) Marchetti, J. L.; Mellichamp, D. A.; Seborg, D. E. Predictive control based on discrete convolution models. Ind. Eng. Chem. Process Des. Dev. 1983, 22, 488. (3) Maurath, P. R.; Laub, A. J.; Mellichamp, D. A.; Seborg, D. E. Predictive controller design by principle component analysis. Ind. Eng. Chem. Res. 1988, 27, 1204. (4) Callaghan, P. J.; Lee, P. L. An experimental investigation of predictive controller design by principle component analysis. Chem. Eng. Res. Des. 1988, 66, 345. (5) Clark, D. W.; Scattolini, R. Constrained receding horizon predictive control. IEE Proc. 1991, 138, 347. (6) Rawlings, J.; Muske, K. Stability of constrained receding horizon control. IEEE Trans. Autom. Control 1993, 38, 1512.
(7) Ohshima, M.; Hashimoto, I.; Takamatsu, T.; Ohno, H. Robust stability of model predictive control (MPC). Int. Chem. Eng. 1991, 31, 119. (8) Lee, J. H.; Yu, Z. H. Tuning of model predictive controllers for robust performance. Comput. Chem. Eng. 1994, 18, 15. (9) Shridhar, R.; Cooper, D. J. A tuning strategy for multivariable model predictive control. Ind. Eng. Chem. Res. 1998, 37, 4003. (10) Al-Ghazzawi, A.; Ali, E.; Nouh, A.; Zafiriou, E. On-line tuning strategy for model predictive controllers. J. Process Control 2000, 11, 265. (11) van Hessem, D. H.; Scherer, C. W.; Bosgra, O. H. LMIbased closed-loop economic optimization of stochastic process operation under state and input constraints. In Proceedings of the 40th IEEE Conference on Decision and Control; IEEE Press: Piscataway, NJ, 2001; p 4228. (12) Chmielewski, D. J. Tuning and steady-state target selection in predictive control: A covariance assignment approach. Presented at the AIChE Annual Meeting, Los Angeles, CA, Nov 12-17, 2000; Paper 257g. (13) Kalman, R. E. When is a linear control system optimal? Trans. ASME J. Basic Eng. 1964, 51. (14) Anderson, B. D. O.; Moore, J. B. Optimal Control: Linear Quadratic Methods; Prentice Hall Inc.: Englewood Cliffs, NJ, 1990. (15) Yokoyama, R.; Kinnen, E. The inverse problem of the optimal regulator. IEEE Trans. Autom. Control 1972, 17, 497. (16) Park, J. G.; Lee, K. Y. An inverse optimal control problem and its application to the choice of performance index for economic stabilization policy. IEEE Trans. Syst., Man Cybernet. 1975, 5, 64. (17) Boyd, S.; Balakrishnan, V.; Feron, E.; El Ghaoui, L. Control system analysis and synthesis via linear matrix inequalities. In Proceedings of the American Control Conference; IEEE Press: Piscataway, NJ, 1993; 2147. (18) Sznaier, M.; Damborg, M. J. Suboptimal control of linear systems with state and control inequality constraints. In Proceedings of the 26th IEEE Conference on Decision and Control; IEEE Press: Piscataway, NJ, 1987; p 761. (19) Chmielewski, D. J.; Manousiouthakis, V. I. On constrained infinite-time quadratic optimal control. Syst. Control Lett. 1996, 29, 121. (20) Scokaert, P. O. M.; Rawlings, J. B. Constrained linear quadratic regulation. IEEE Trans. Autom. Control 1998, 43, 1163. (21) Zheng, A. Z.; Morari, M. Stability of model predictive control with mixed constraints. IEEE Trans. Autom. Control 1995, 40, 1818. (22) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. M. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789. (23) Scokaert, P. O. M.; Mayne, D. Q.; Rawlings, J. B. Suboptimal model predictive control (feasibility implies stability). IEEE Trans. Autom. Control 1999, 44, 648. (24) Skelton, R. E.; Iwasaki, T.; Grigoriadis, K. M. A Unified Algebraic Approach to Linear Control Design; Taylor & Francis: New York, 1999. (25) Dullerud, G. E.; Paganini, F. A Course in Robust Control Theory: A Convex Approach; Springer-Verlag: New York, 1999. (26) Boyd, S. P.; El Ghaoui, L.; Feron, E.; Balakrishnan, V. Linear Matrix Inequalities in System and Control Theory; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, 1994. (27) VanAntwerp, J. G.; Braatz, R. B. A tutorial on linear and bilinear matrix inequalities. J. Process Control 2000, 10, 363. (28) Russell, D. L. Mathematics of Finite-Dimensional Control Systems; Marcel Dekker: New York, 1979. (29) Kalman, R. E.; Koepcke, R. W. Optimal synthesis of linear sampling control systems using generalized performance indexes. Trans. ASME 1958, 80, 1820. (30) Luenberger, D. G. Optimization by Vector Space Methods; John Wiley & Sons: New York, 1969. (31) Burl, J. B. Linear Optimal Control: H2 and H∞ Methods; Addison-Wesley: New York, 1999.
Received for review August 22, 2003 Revised manuscript received June 25, 2004 Accepted August 30, 2004 IE030686E
