Model Predictive Control of Hammerstein Systems with Multivariable

168

Ind. Eng. Chem. Res. 2007, 46, 168-180

Model Predictive Control of Hammerstein Systems with Multivariable Nonlinearities Kwong Ho Chan and Jie Bao* School of Chemical Sciences and Engineering, The UniVersity of New South Wales (UNSW), Sydney, NSW 2052, Australia

In this paper, a model predictive control (MPC) design method is proposed for Hammerstein systems with multivariable nonlinearities. Iterative inversion is introduced such that classical linear MPC design can be used. The condition for which the numerical inversion is guaranteed to converge is derived. From this condition, the input space at which the nonlinear function is invertible can be obtained, which poses extra input constraints for MPC design. The overall constraints are transformed and approximated by a set of linear constraints by means of convex hull construction. Illustrative examples are presented to demonstrate the effectiveness of the proposed control method. 1. Introduction In many process control applications, the key nonlinearities and dynamics of processes can be sufficiently approximated using linear dynamic models with static nonlinearities.1 This type of process model is called a Hammerstein model, which has a block structure, as shown in Figure 1. It consists of the input nonlinearity σ:Rm f Rm and the linear dynamics H:Rm f Rm. Because many highly nonlinear processes can be sufficiently modeled using Hammerstein models (see, for example, the work of Tian et al.2), control of Hammerstein processes is an area of interest for control engineers. There are quite a few control methods for single-input-single-output (SISO) and multi-input-multi-output (MIMO) Hammerstein systems reported in the literature (see, for example, refs 3-9). Most of these methods use direct inversion to remove the static nonlinearities and convert the nonlinear control problem to a linear one. This approach is often effective for SISO processes, because the analytical inverse of static nonlinearity can often be obtained. However, for most MIMO processes with multivariable static nonlinear functions, the analytical inverse of the nonlinearity usually does not exist. Even if it does exist, it is often very difficult to find. Therefore, many of the inversion approaches can only be applied to process systems with decoupled or very simple nonlinearity.3-6,10,11 Among the existing methods, model predictive control (MPC)3-6,9 is most widely used in control practice,12,13 because it can handle process constraints. Patwardhan et al.3 and Fruzzetti et al.4 decoupled the nonlinearity so that the inverse and transformed constraints can be obtained easily. Abonyi et al.5 approximated the nonlinearity functions using fuzzy models, whose inverse can be obtained in a straightforward analytical procedure. Bloemen et al.6 transformed the nonlinearities to polytopic descriptions such that robust linear MPC techniques can be used to control these systems. The main disadvantage of this method is the rapid increase of the number of linear matrix inequalities (LMIs) for general nonlinearities. To address the static nonlinearity, Haddad and Chellaboina7 introduced a nonlinear compensator, to ensure closed-loop stability. However, this approach only applies to a small family of input nonlinearity σ(u) ) [σ1(u), ..., σi(u), ..., σm(u)]T (where * To whom correspondence should be addressed. Tel: +61 (2) 93856755. Fax: +61 (2) 9385-5966. E-mail: [email protected].

Figure 1. Hammerstein model.

u ) [u1, ..., ui, ..., um]T) such that σi(u) ) 0 if ui ) 0. Recently, Chan et al.8 proposed an extended internal model control (IMC) framework where the linear and nonlinear subsystems are addressed by multiple IMC controllers. The concept of passive systems was used to construct the IMC controllers, which approximate the inverses of the linear and nonlinear subsystems to achieve dynamic control performance. This approach applies to multivariable coupled nonlinearities. However, it cannot handle constraints, which are present in virtually all real-life applications. For MIMO Hammerstein processes, handling the constraints is not straightforward. The physical input constraints must be transformed to the nonlinear constraints for the linear MPC. The existing MPC methods mentioned previously often assume simple and/or decoupled nonlinearity, so that constraints can usually be transformed to a new set of linear constraints merely by inspection. However, when the nonlinearity is more complex (e.g., being multivariable), the transformed constraints are nonlinear and cannot be used in linear MPC design directly. The entire problem becomes more complicated if an analytical inverse does not exist. In this paper, an MPC approach is developed for Hammerstein processes with multivariable nonlinearities, which are often encountered in chemical process control. The proposed approach developed in this paper uses iterative numerical inverse, instead of relying on analytical solutions. The condition for which the numerical inversion is guaranteed to converge is derived. Based on these conditions, the input variable regime in which the numerical inverse can be obtained is identified. This constitutes new constraints on the input space, in addition to the physical input constraints. The total input constraints are transformed to linear input constraints using polytopic descriptions and incorporated in the MPC design. This paper is organized as follows: First, the framework is briefly explained in section 2. The constraint transformation is then described in section 3. Section 4 describes the control synthesis procedure. Examples are presented in section 5, to illustrate the implementation and

10.1021/ie0609113 CCC: $37.00 © 2007 American Chemical Society Published on Web 12/07/2006

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 169

Figure 2. Model predictive control (MPC) with direct inverse.

R, with V(0) ) 0, called the storage function, such that for all u(k) ∈ Rm and all k ) {0, 1, 2, ...},

V(x(k + 1)) - V(x(k)) e w(u(k),y(k))

(4)

In particular, if the system is dissipative, with respect to the following supply rate,

w(u(k),y(k)) ) y(k)Tu(k) Figure 3. One-degree-of-freedom MPC architecture.

effectiveness of the proposed method, followed by the conclusion in section 6.

Consider the Hammerstein model in Figure 1 to be represented by

x(k + 1) ) Ax(k) + Bσ(u(k))

(1a)

y(k) ) Cx(k)

(1b)

where x ∈ u∈ y∈ and the parameters A, B, and C are stabilizable and detectable. Assume the nonlinear mapping σ:Rm f Rm is such that σ(0) ) 0, and none of the eigenvalues of A are equal to 1. If an inverse σ-1 exists and eq 1 represents the true process, then, by implementing a direct inverse, a linear model predictive controller can be used (as shown in Figure 2). In this case, the configurations in Figures 2 and 3 are equivalent. In general, analytical inverse σ-1 does not often exist. Therefore, a numerical inverse approach is necessary. 2.1. Nonlinear Subsystem. 2.1.1. Numerical Inverse. In this section, the issue of finding the inverse by numerical iteration is addressed. The conditions of convergency of the numerical iteration are established on the basis of the passivity theory.14-19 The following preliminaries are useful in our discussion. Definition 1 (Zero-State Detectability and Observability).14 Consider a nonlinear discrete-time system Ω of the form Rn,

then it is said to be passive. The passivity of linear passive systems can be determined by the following condition. Theorem 1. Necessary and Sufficient Condition for Discrete PositiVe Realness.16 Consider a linear discrete-time system H

H:

2. Framework

Rm,

Ω:

R m,

{

x(k + 1) ) f(x(k), u(k)) y(k) ) h(x(k), u(k))

(2)

with zero input (u(k) ) 0) and let F ⊂ Rn be its largest positively invariant set contained in

x(k + 1) ∈ Rn|y(k) ) h(x(k + 1),0) ) 0

(3)

Ω is zero-state detectable (ZSD) if x(k) ) 0 is asymptotically stable conditionally to F. If F ) {0}, then Ω is zero-state observable (ZSO). Definition 2 (Passivity).15 The nonlinear system described in eq 2 is said to be dissipative with a supply rate w(u(k),y(k)) (which is a function of the input and output) if there exists a positive semidefinite C1 function of state variables V(x): Rn f

(5)

{

x(k + 1) ) Ax(k) + Bu(k) y(k) ) Cx(k) + Du(k)

(6)

with a transfer function that is giVen by

H(z) ) C(zI - A)-1B + D

(7)

where the parameters A, B, C, and D represent the stabilizable and detectable minimal realization, and where no poles of H(z) are in |z| > 1 (simple poles only on |z| ) 1). The necessary and sufficient conditions for system H to be passiVe (or, equiValently, H(z) to be discrete positiVe real (DPR)) are that there exists a real symmetric positiVe definite matrix P and real matrices L and W, such that

P - ATPA ) LTL

(8)

C - A PB ) L W

(9)

D + D - B PB ) W W

(10)

T

T

T

T

T

T

Remark 1. When H(z) have all poles on |z| ) 1, the conditions for H(z) to be DPR are given in Theorem 1 with L ) 0.17 For static nonlinearity, the state space is void, i.e., V(x(k)) ) 0 for all k ) {0, 1, 2, ...}. Therefore, definition 2 is reduced to the following condition. Definition 3 (Passive Static Nonlinearity).18 A static nonlinearity σ: Rm f Rm is said to be passive if

σ(u)Tu g 0

∀u ∈ Rm

(11)

Passivity is an input/output property. A process can have excessive passivity or a shortage of passivity, which can be defined as follows. Definition 4 (Excess/Shortage of Passivity).14 A system Ω of the form in eq 2 with storage function V is said to be input feedforward passive (IFP) if it is dissipative, with respect to the supply rate,

w(u(k),y(k)) ) y(k)Tu(k) - νu(k)Tu(k)p for some ν ∈ R, or, equivalently, when system Ω with constant

170

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

V2(x1(k + 1)) - V2(x1(k)) ) 0 e u2(k)Ty2(k) ν*u2(k)Tu2(k) - Fy2(k)Ty2(k) (17) Define ∆V ) V1(x1(k + 1)) - V1(x1(k)). To simplify the equations, the counting index k is dropped in the following derivation. Adding eqs 16 and 17 gives

Figure 4. Numerical inverse.

negative feedforward νI is passive; or it is said to be output feedback passive (OFP) if it is dissipative, with respect to the supply rate T

for some F ∈ R, or, equivalently, when system Ω with constant positive feedback FI is passive. If ν > 0 (or F > 0), then the system is in excess of input feedforward (or output feedback) passivity. In this case, the system is said to be strictly input (or output) passive. If ν < 0 (or F < 0), then the system is in shortage of input feedforward (or output feedback) passivity. A system can also be characterized by IFP and OFP simultaneously, e.g., be dissipative, with respect to the supply rate of w(u(k), y(k)) ) y(k)Tu(k) νu(k)Tu(k) - Fy(k)Ty(k). The condition on the convergence of the iterative numerical inverse can be derived from the aforementioned passivity conditions. Consider the following iterative procedure:

u(k + 1) ) u(k) + B(q - σ(u(k)))

(12)

where B is nonsingular (e.g., B ) ζIm×m). If the iteration does converge, as k f ∞, u(k + 1) ) u(k) ) u(∞). Hence,

B(q - σ(u(∞))) ) 0

(13)

u(∞) ) σ-1(q)

(14)

i.e.,

The aforementioned iteration calculates u ) σ-1(q). It is equivalent to the feedback system depicted in Figure 4, where K is a discrete-time integrator (which can be described by eq 6 with A ) C ) I and D ) 0). Clearly, the convergence of the iteration is equivalent to the closed-loop stability of the system in Figure 4. It can be shown that the closed-loop system is stable if the following passivity-based condition is satisfied. Proposition 1. Consider a linear system K:u1 f y1:

x1(k + 1) ) Ax1(k) + Bu1(k) y1(k) ) Cx1(k) + Du1(k)

(18)

e uT1 u2 - νuT1 u1 + uT2 (-u1) - ν*yT1 y1 - F(-u1)T(-u1) (19)

w(u(k),y(k)) ) y(k) u(k) - Fy(k) y(k) T

∆V e uT1 y1 - νuT1 u1 + uT2 y2 - ν*uT2 u2 - FyT2 y2

(15)

which is ZSD, and a static nonlinearity σ:u2 f y2. If K has IFP(ν) and σ is dissipative, with respect to the supply rate, uT2 y2 - ν*uT2 u2 - FyT2 y2, then the equilibrium of the closed-loop configuration in Figure 4, x1 ) 0, is asymptotically stable (when q(k) ) 0) if ν* > 0 and ν + F > 0. Proof. Consider Figure 4, where q(k) ) 0 for all k ) {0, 1, 2, ...}. Because K has IFP(ν), we have

V1(x1(k + 1)) - V1(x1(k)) e u1(k)Ty1(k) - νu1(k)Tu1(k) (16) For the nonlinear static subsystem, because its state space in eq 2 is void, the only choice for its storage function is V2(x(k)) ) 0 for all k ) {0, 1, 2, ...}. Because σ is dissipative, with respect to the supply rate uT2 y2 - ν*uT2 u2 - FyT2 y2, then

e -(ν + F) uT1 u1 - ν*yT1 y1

(20)

If ν* > 0 and ν + F > 0, then

∆V < 0

∀ y1 * 0

(21)

The bounded solution of x1 is confined in {x1|y1 ) 0}. Because K is ZSD, which implies x1 f 0 when y1 f 0, the equilibrium of the closed-loop system, x1 ) 0, is asymptotically stable (when q(k) ) 0). 9 Consider the diagonal discrete-time integrator K shown in Figure 4 with the following transfer function:

K ) ζ[(z - 1)-1I]

(22)

where ζ > 0 is the effective integrator gain. This is equivalent to the iteration in eq 12 with B ) ζI. Denote K* as a system with a positive feedforward of 1/2ζI, (i.e., K* ) K + 1/2ζI); K* then can be described by eq 6 with A ) C ) I, B ) ζI, and D ) 1/2ζI. According to Theorem 1, K* is DPR with L ) W ) 0 and the positive definite matrix P ) 1/ζI (ζ > 0). Therefore, the integrator K has IFP(-ζ/2) by Definition 4. From Proposition 1, if the static nonlinearity σ is dissipative, with respect to the supply rate uTy - ν*uTu - FyTy (ν* > 0, F > 0), then the integrator gain ζ must be small enough such that it satisfies the condition ζ < 2F to guarantee the asymptotical stability of the closed-loop system in Figure 4. Proposition 1 shows the condition for which an accurate numerical inverse, with respect to q(k) ) 0, can be obtained. With the iterative inverse, an MPC control structure given in Figure 2 can be implemented. However, uOPT in Figure 2 is calculated at every time step k. Hence, the setpoint q(k) ) uOPT in Figure 4 is different at every kth step. Therefore, the stability condition given in Proposition 1 must be satisfied at different setpoints q(k) ) q* for the closed-loop system. Define ∆x1 ) x1 - x1*, ∆u1 ) u1 - u1*, ∆y1 ) y1 - y1*, ∆u2 ) u2 - u2*, and ∆y2 ) y2 - y2*, where x1*, u1*, u2*, y1* and y2* are steadystate values of x1, u1, u2, y1, and y2 that correspond to q(k) ) q* ) σ(u2*). When q* * 0, if ∆x1 ) 0 is asymptotically stable, then ∆u2 f 0, which implies that the iteration converges to u2* ) σ-1(q*). The linear system K will not be affected by the aforementioned coordinate transformation; however, the nonlinear function σ will be. For the nonlinear function, the passivity of the mapping ∆σ:∆u2 f ∆y2 must be considered. This requires the introduction the following definition. Definition 5 (Incremental Passivity).19 A static nonlinearity σ:u f σ(u) is said to be strictly incrementally passive if

(u - u*)T(σ(u) - σ(u*)) > 0

∀ u * u* and u, u* ∈ Rm (23)


Lemma 1. Consider a static nonlinear function σ:u f y and the function deriVed from coordinate transformation, where ∆u ) u - u* and ∆y ) σ(u) - σ(u*). If σ is strictly incrementally passiVe and bounded, i.e.,

σ(u)Tσ(u) < γ < ∞

∀u ∈ Rm

(24)

where γ is a constant; then, the function is dissipatiVe, with respect to the supply rate, ∆uT∆y - ν*∆uT∆u - F∆yT∆y, for all setpoints (u*, σ(u*)) ∈ Rm, such that ν*, F > 0. Proof. Because σ is strictly incrementally passive (i.e., satisfies eq 23),

∆uT∆y )

(for some > 0)

(25)

Because both u and σ(u) are bounded, there always exist constants η and γ1 such that 0 < ∆uT∆u < η < ∞ and 0 < ∆yT∆y < γ1 < ∞. Therefore, there always exist 1 and 2, such that

∆uT∆y ) (∆u ∆u)(∆yT∆y) (∆uT∆u)(∆yT∆y) T

(26)

1 2 ∆uT∆y > T + T T (∆u ∆u)(∆y ∆y) ∆y ∆y ∆u ∆u (for some 1, 2 > 0) (27) T

Hence, ν* and F can be chosen to be 1 and 2 (∀u, u* ∈ respectively, leading to

∆uT∆y - ν*∆uT∆u - F∆yT∆y > 0

Rm),

(28)

i.e., ∆σ is dissipative, with respect to the supply rate given in eq 28 for all setpoints (u*, σ(u*)), such that ν*, F > 0. 9 Proposition 2. If a nonlinear system σ is strictly incrementally passiVe and bounded, then the iteratiVe numerical inVerse σ-1 exists, using eq 12. Proof. Because σ is strictly incrementally passive and bounded, then, from Lemma 1, ∆σ is dissipative, with respect to the supply rate ∆uT∆y - ν*∆uT∆u - F∆yT∆y for all setpoints (u*, σ(u*)) such that ν*,F > 0. Hence, a discretetime integrator of the form in eq 22 with 0 < ζ < 2F can be used to obtain the numerical inverse σ-1 iteratively, according to Proposition 1. 9 For most multivariable nonlinearities defined on an input space U ⊂ Rm, there usually exists a region Up ⊂ U where the nonlinearity is strictly incrementally passive. Therefore, the first step is to identify the region Up by numerically testing the condition given in eq 23. If the nonlinear function is bounded, the input space Up guarantees the convergency of the inverse estimate. As shown in eq 25, if the strictly incrementally passive condition is satisfied, we have ∆uT∆y > 0. Therefore, there always exists a positive ) ∆uT∆y. Consequently, there exist positive 1 and 2 such that eq 27 is satisfied. Rather than examine both variables simultaneously, 1 (and, hence, ν*) is chosen to have a very small and positive value, which leaves only the condition F + ν > 0 to be considered. The next step is to find the smallest F for all u, u* ∈ Up (u * u*). For a given ν*, this can be performed by solving the following problem numerically.

Problem 1.

∆uT∆y - ν*∆uT∆u u,u* ∈ Up ∆yT∆y

F ) min

(29)

The integrator in eq 22 then can be designed. Clearly, the larger the value of ζ, the faster the convergence of the closed loop. However, notice that the aforementioned integrator has IFP (-ζ/ 2). From Proposition 1, the gain ζ must be chosen to be small enough such that it satisfies the condition ζ/2 < F for the closed loop in Figure 4 to be asymptotically stable. Hence, it limits how fast the estimate will converge. The iterative numerical inversion must be completed in each sampling period of the MPC control. Therefore, a large number of iterations may be required within each sampling period to obtain an accurate inverse approximation. This would not be a serious practical problem, because MPC is usually used as a second-level controller, which often has sampling periods on the order of minutes. Thousands of iterations would be possible within each sampling period, to give an accurate estimate of the nonlinear inverse. 2.1.2. Input/Output Scaling. Because passivity is a phaserelated property, the input and/or output of a given process can be scaled to “improve’’ its passivity.20 For example, consider a very simple SISO Hammerstein system,

H (z) )

-z z-a

F(u) ) -u

(30) (31)

where |a| < 1. It is obvious that both subsystems are not passive and, thus, Proposition 1 cannot be applied. However, simply multiplying both subsystems by -1 on the input side gives two passive subsystems, to which the proposed numerical inverse can be applied. For MIMO processes, scaling matrix operators are introduced (as shown in Figure 5). A static nonsingular real matrix DJ can be chosen such that [σ(DJu)]Tu > 0 at the desired setpoint. The numerical inverse can be designed based on σD(u) ) σ(DJu). There are, of course, many choices for DJ. Define

Ja )

|

∂σ ∂u u)uj

(32)

to be the nonsingular Jacobian matrices of σ around the desired steady state. One obvious choice is DJ ) Ja-1, because the new steady-state gain matrix of σD at the desired setpoint become the identity matrix Im×m, which satisfies the strict passivity condition at steady state. This input/output scaling treatment can significantly reduce the conservativeness of the proposed control design approach. 2.2. Linear Subsystem. Consider the classic MPC design in Figure 3, where QP is the optimizer using quadratic programming (see Appendix A), and the observer will calculate the estimated current states and disturbances. Notice that the input to the observer is e(k) ) y(k) + (d(k) - ys) from Figure 3 (where ys denotes the setpoint). Assuming the model is perfect and the inverse σ-1 exists, the same expression for e(k) holds from the configuration in Figure 5. Hence, a linear MPC will stabilize the closed loop and drive the process output to the desired setpoint, provided that the model is perfect and inverse is exact. However, such an MPC controller cannot be designed easily mainly due to the existence of constraints. 3. Constraints Using Polytopic Descriptions 3.1. Input Constraints. Variables in all physical systems are inherently bounded. It is well-understood that if the constraints

172


Figure 5. Modified MPC structure with numerical inverse.

Figure 6. Constraints transformation.

Figure 7. Graphical depictions of the two approaches to be used: (left) Approach 1 and (right) Approach 2.

of the variables are not considered in the design of controllers, the performance of the controllers can be very poor, particularly when the system is operated near the constraints. Addressing input/output constraints is one of the main motivations behind the development of MPC. With the presence of the input nonlinearity, the physical input constraints will have to be converted to the limits on the linear MPC output. The output constraints can be addressed by the linear MPC itself. The existing MPC designs for Hammerstein models focus on obtaining an inverse of the nonlinearity. However, very little has been reported on how to transform the actual physical constraints to that of the controller output constraints of the MPC. Let KMPC denote the operator of the MPC controller. From Figure 5, u ) DJσD-1(q). Therefore, q ) σD(DJ-1u), where DJ is the scaling matrix defined in eq 32. Given the physical constraints u ∈ U ⊂ Rm, the corresponding permissible region of the linear MPC output q ∈ Q ⊂ Rm could be determined based on σD (see Figure 6). This region should be approximated with linear constraints in the form of Lq e b (L ∈ Rd×m and b

∈ Rd, for some positive integer d) so that the constraints can be immediately used in linear MPC design. To simplify the controller design, optimization of a convex problem is preferred. Hence, the constraints q ∈ Q should also be linear and convex (Q ) {q|Lq e b}). Generally, the constraints q ∈ Q are arbitrary. To examine this problem, two approaches were considered (see Figure 7): Approach 1, which involves finding the largest set of convex constraints such that it is within the original set of constraints, and Approach 2, which involves finding the smallest set of convex constraints such it contains the original set of constraints. Approach 1 guarantees that the controller output will satisfy the input constraints. However, it is a conservative approach, because it does not fully utilize all the operating conditions within the constraints. In addition, obtaining such a set of constraints is often very difficult, because of the nonconvex nature of the problem. Approach 2 is much more amendable, through the use of linear programming (LP). Furthermore, the controller action will not be conservative, because it covers the desired constraints. However, sometimes it will actually violate


the desired constraints. In this case, a saturation function will be used. The implication is a loss of optimal operating conditions. 3.2. Passivity Contraints. As discussed in Section 2.1.1, σD must be strictly incrementally passive, so that the iterative numerical inversion can converge. Unless the Jacobian matrix in eq 32 is singular, there always exists a region Up in the neighborhood of a desired steady state where eq 23 is satisfied. A search program can easily be written to determine the region Up and the values 1 and 2 by simply checking the condition in eq 27 for all u, u* ∈ Up. Therefore, an extra set of constraints u ∈ Up is introduced to guarantee the convergency of the inverse estimate. Following Approach 2 in Section 3.1, the transformed constraints q ∈ Qp can be determined where Qp ) {q|Lpq e bp}. Hence, the overall set of constraints on q becomes Qs ) {q|Lsq e bs}, where Ls ) [LT,LTp ]T and bs ) [bT,bTp ]T. 3.3. Convex Hull. This section shows how the set of constraints Qs ) {q|Lsq e bs} is obtained. Usually, the transformed region Qs is bounded by curves rather than “flat’’ facets. To approximate curves accurately, Us ) U∩Up (where U and Up represent the process input spaces confined by the physical constraints and the passivity condition, respectively) is divided into a mesh grid U* and transformed to Q*, which results in many piecewise-continuous “flat’’ facets on the boundaries. This will result in a large matrix Ls and vector bs. Finding such Ls and bs from Q* is a classical problem in computational geometry. For a subset S of Rm, the convex hull, denoted as conv(S), is defined as the smallest convex set in Rm that contains S. The convex hull computation means the determination of conv(S) for a given finite set of R points S ) {p1, p2, ..., pR} in Rm. One way to determine conv(S) is to represent it as the intersection of half-spaces, or, more precisely, as a set of solutions to a minimal system of linear inequalities, which is called an H-representation, or, equivalently, conv(S) ) {q|Lsq ebs}. When conv(S) is full-dimensional, each (nonredundant) inequality corresponds to a facet of conv(S). The convex hull computation can also be defined as the determination of extreme points of conv(S) (called a V-representation), or, equivalently, the determination of redundant points in S to determine conv(S). This can be done by solving O(R) (order of R time complexity) linear programs and, thus, is polynomial solvable. A standard method to check whether a given point p* is in conv(S) uses the linear programming (LP) technique. An LP problem is formulated for the question as follows. Problem 2.21 Let S ) {p1,p2,...,pR}. For a given p*, find scalars λi such that, for all i ) 1, ..., R, R

p* )

λipi ∑ i)1

(33)

R

λi ) 1 ∑ i)1

(34)

λi g 0

(35)

This is a linear feasibility problem, which is polynomially equivalent to a general LP. The feasibility problem has a solution if and only if the following has no solution. Problem 3.21 Find z0 ∈ R and z ∈ Rm such that for all i ) 1, ..., R,

zTpi e z0

(36)

zTp* > z0

(37)

If there is a solution (z0,z), then the set H ) {x ∈ Rm: zTx ) z0} is a hyperplane in Rm separating the polytope conv(S) from the inquiry point p*. Therefore, the existence of the separation means that p* is not redundant. To actually solve the problem, the following LP is set up. Problem 4.21 Find f* ) maxz,z0 zTp* - z0 such that for all i ) 1, ..., R,

z T pi - z 0 e 0

(38)

zTp* - z0 e 1

(39)

Equation 39 is artificially added so that the LP has a bounded solution. It is easy to see that the point p* is nonredundant if and only if the optimal value f* of the LP is strictly positive. Therefore, conv(S) can be determined using the following procedure: (1) Create a simplex of m + 1 points (convex hull of m + 1 affinely independent points). (2) Redundant points can be identified by solving Problem 4. (3) Create a larger convex hull with nonredundant points. (4) Remove redundant points from the new convex hull. (5) Identify non-redundant points and create a new convex hull. (6) Repeat Steps (4) and (5) until there is no redundant points. One of the many ways to construct a convex hull is to use the Quickhull algorithm,22 which is based on the following theorem. Theorem 2 (Simplified Beneath-Beyond).23 Let H be a conVex hull in Rm, and let p* ∈ Rm be a point outside H. F then is a facet of conV({p*} ∪ H) if and only if (1) F is a facet of H, and p* is below F; or (2) F is not a facet of H, and its Vertices are p* and the Vertices of a ridge of H with one incident facet below p* and the other incident facet aboVe p*. An outline of an algorithm for constructing a convex hull is presented in Appendix B. The algorithm will give a convex hull in the V-representation. Ls and bs can be obtained easily by finding the hyperplane equation through the m vertices of each facet. Typical industrial processes may have hundreds or even thousands of manipulated variables. Therefore, the approximation of the constraints is very important to reduce the number of constraints handled by the controller. One approach is to identify and use the important constraints only. Ideally, the most important subset of constraints should be determined such that the exclusion of other constraints will lead to the least deterioration of the control performance. However, this is highly impractical, because this analysis itself is more computationally complex than implementation of the full set of constraints and certainly cannot be performed online. A more practical way is to construct a new convex hull that has the smallest “volume’’ from a subset of the original constraints. The steps are as follows: (1) Consider that the convex hull H is represented in H-representation with Ls and bs. A new convex hull Hi ) {q|L h iq e bhi} is formed by omitting the ith constraint (ls,i1q1 + ‚‚‚ + ls,ikqk ebs,i) for i ∈ 1, 2, ..., ds. (2) The “volume’’ of the new convex hull (vol(Hi)) is calculated using the Multi-Parametric Toolbox24 in MATLAB for every i ∈ 1, 2, ..., ds. The larger vol(Hi), the more important

174


the missing ith constraint. Let H* be the convex hull such that vol(Hi) is a minimum for i ∈ 1, 2, ..., ds.

Example 1: Input Constraints. Consider the following Hammerstein system in the form of eq 1, presented as follows:

[

(3) Assume that dm is the maximum number of constraints that the controller can handle. Let dn ) min [dm, ds] and H ) H*. Repeat Steps (1) and (2) until ds ) dn.

0.3780 0.1331 -0.1651 -0.4715 0.2182 0.2275 0.5139 -0.1475 A) 0.4626 -0.0954 -0.0255 0.2939 -0.0753 -0.5538 0.1346 -0.1260

4. Control Design Procedure The control design procedure can be summarized as follows: (1) Identify process models, the static nonlinearity σ, and linear dynamics H.

C)

(2) Calculate the scaling matrix according to eq 32. (3) Find the region Up by testing σD(DJ-1u) from eq 23. (4) For the given input constraints U, calculate the F value required for the region Us ) Up ∩ U from Problem 1. (5) From the F value obtained from Step (4), design the integrator K. The effective gain ζ must be chosen such that F + ζ/2 > 0. (6) Map the mesh grid U* that represents the overall constraints Us into the corresponding region Q* because q ) σD(DJ-1u). (7) Calculate the approximated convex set of constraints that contains the transformed constraints, as described in Section 3.3. (8) The linear MPC controller KMPC is designed based on H using the linear MPC design in Appendix A, together with the approximated constraints in MPC design in Step 7. (9) As shown in Figure 5, the overall controller includes the linear MPC, the iterative inverter and the static rescaling matrix. Saturation functions should be implemented immediately after the rescaling matrix to ensure that the physical constraints are not violated. Remark 2. The proposed approach works for Wiener systems in principle. In this case, it is the output (rather than the input) constraints that must be transformed and approximated. However, when dealing with such nonlinear systems with output disturbance d, the error becomes e(k) ) σD-1(y(k) + d(k)) σD-1(ys). The relationship between the error and disturbance is no longer linear, which is problematic in linear MPC design, where the unknown disturbance must be estimated using a linear observer (as shown in Figure 3). Therefore, the proposed approach cannot be directly applied to Wiener systems unless no output disturbance exists.

[

-10.6197 9.9973 11.5842 -10.4204 B) -0.5028 0.6781 -7.9532 7.3762

[

]

0.4723 0.8878 0.1904 0.5784 0.5252 0.4876 0.6927 0.0413

]

]

and

σ(u) )

[ ] σ1(u) σ2(u)

σ1(u) ) u1 + 0.2u1u2 - 0.1u1u22 + u2

(40)

σ2(u) ) -u1 - 0.2u1u2 + 0.5u12u2 + 2u2

(41)

with input constraints, i.e., input space U,

-3 e u1 e 3

(42)

-3 e u2 e 3

(43)

Following the procedures in section 4, it is determined that no scaling matrix is required, i.e., σD(u) ) σ(u). The strictly incrementally passive region Up is then identified numerically by testing the condition given in eq 23. The overall region Us ) U ∩ Up is determined and illustrated in Figure 8. Further-

Remark 3. The extra constraints to ensure invertibility can lead to a certain level of conservatism, because the incremental passivity condition is a sufficient (but not a necessary) condition for invertibility. The degree of conservatism is dependent on the curvature of the nonlinear functions. Quantification of the conservatism of the incremental passivity condition requires further investigation.

5. Illustrative Examples To illustrate the proposed method, two case studies are performed: Example 1, which involves only input constraints, and Example 2, which involves both input and output constraints.

Figure 8. Graphic depiction of Example 1, showing the overall constraints Us ) Up ∩ U: U, boundary of the graph; Up, dotted region.


Figure 9. Graphic depiction of Example 1, showing the transformed region: σD(DJ-1u), where u ∈ Us (shown in Figure 8). Note: Circles are vertices of the convex hull of the transformed region.

Figure 10. Approximated region of Example 1.

more, F < 0.0048 is obtained by solving Problem 1. Hence, it is required that ν g -0.0048. Therefore, if the integrator K in Figure 4 is chosen to be

[

0.0096 0 K(z) ) z - 1 0.0096 0 z-1

]

(44)

for which ν ) -ζ/2 ) -0.0048, the closed loop with σ(u) ) σD(u) in Figure 4 will be asymptotically stable for all u ∈ Us. The tolerance of the iterative inversion is set to be 0.01, which, in this case, requires ∼2000 iterations per sampling period of the MPC. Existing MPC methods cannot handle the given constraints, because of the nonlinearity transformation. Hence, the proposed approximation method, which provides constraints in linear form (Lsq e bs), is used. The transformed region bounded by the constraints is depicted approximately in Figure 9. Following the approximation method in section 3.1, the approximated constraints (see Figure 10) are obtained using the Multi-Parametric Toolbox24 in MATLAB:

L s q e bs where

[ ] []

-0.9890 0.1477 -0.9989 0.0478 0.9461 -0.3239 Ls ) 0.9806 -0.1961 -0.8721 -0.4893 0.7198 0.6942 3.7661 3.1401 3.9462 bs ) 3.6713 5.3590 8.8224

Figure 11. Controller output of Example 1 (using the proposed method).

The linear MPC design described in Appendix A is applied with an error weighting matrix of Ψ ) 1000I, an input weighting matrix of Φ ) I, a control horizon of M ) 20, and a prediction horizon of N ) 30 in eq 56. In addition, the terminal state weighting matrix Ψf in eq 56 is chosen to be the solution to the following discrete-time Riccati equation:

ATΨfA - Ψf - (ATΨfB)(Φ + BTΨfB)(BTΨfA) + Ψ ) 0 (45) Let the setpoint be [1, -1]T. The responses are shown in Figures 11 and 12. To illustrate the effectiveness of the proposed constraints approximation approach, it is compared with a rough rectangular

176


Figure 14. Output response of Example 1 (using the minimum-maximum constraints).

Figure 12. Output response of Example 1 (using the proposed method).

An MPC control is designed on the basis of the aforementioned input constraints. The responses are shown in Figures 13 and 14. The proposed MPC design performs much better than the MPC controller with minimum-maximum constraints. The responses for the MPC design with minimum-maximum constraints give an integral square error (ISE) of 9.01, whereas the ISE for the proposed MPC design is 3.79 over 60 sampling periods. Example 2: Input and Output Constraints. This example shows how the proposed approach handles both the input and output constraints. As mentioned in the paper, the output constraints can be addressed by the linear MPC itself. The presence of both constraints will lead to more-conservative control action. Consider the following Hammerstein system, which has the same linear part as the Hammerstein system in Example 1:

σ1(u) ) 4u1 + 0.5u1u2 - 0.1u1u22 + u2

(48)

σ2(u) ) -u1 - 0.2u1u2 + 0.1u1 u2 + 2u2

(49)

2

with input constraints Figure 13. Controller output of Example 1 (using the minimum-maximum constraints).

approximation of the transformed region. Consider the transformed region in Figure 9 to be of the form

-6 e u1 e 6

(50)

-6 e u2 e 6

(51)

and output constraints

q1,min e q1 e q1,max

(46)

0 e y1 e 1.1

(52)

q2,min e q2 e q2,max

(47)

-1.1 e y2 e 0

(53)

The minimum and maximum values of the transformed input variables are determined to be q1,min ) -3.13, q2,min ) -11.26, q1,max ) 4.91, and q2,max ) 14.36, i.e., Lmq e bm, where

[ ] [ ]

1 0 0 1 Lm ) -1 0 0 -1 4.91 14.36 bm ) -3.13 -11.26

with the setpoint [1, -1]T, i.e., a maximum 10% overshoot. Following the procedures in section 4, it is observed that no scaling matrix is required. The region Up is then identified numerically. The overall region Us is determined and illustrated in Figure 15. Furthermore, F < 0.025 is obtained by solving Problem 1. Hence, it is required that ν g -0.025. Therefore, the integrator K in Figure 4 is chosen to be

[ ]

0.05 0 z K(z) ) - 1 0.05 0 z-1

(54)

The tolerance of the iterative inversion is set to be 0.01, which,


Figure 15. Graphic depiction of Example 2, showing the overall constraints Us ) Up ∩ U: U, boundary of the graph; Up, dotted region.

Figure 16. Graphic depiction of Example 2, showing the transformed region: σD(DJ-1u), where u ∈ Us (shown in Figure 15). Note: Circles are vertices of the convex hull of the transformed region.

in this case, requires ∼1000 iterations per sampling period of the MPC. The transformed region bounded by the constraints is depicted approximately in Figure 16. Following the approximation method in Section 3.1, the approximated constraints (see Figure 17) are obtained using the Multi-Parametric Toolbox24 in MATLAB:

L s q e bs where

[ ] []

0.9670 0.2549 -0.9721 -0.2345 -0.8553 0.5182 Ls ) 0.8122 -0.5833 -0.3693 -0.9293 0.5433 0.8396 32.3992 21.9710 36.9777 bs ) 23.2275 15.0450 31.4691

Figure 17. Approximated region of Example 2.

The linear MPC is designed according to Appendix A, with Ψ ) 1000I, Φ ) I, and M ) N ) 5 in eq 56. To illustrate how the proposed approach handles the output constraints, the control and prediction horizons are chosen to be small so that a large overshot in the output will appear (and, thus, violate the output constraints) if no output constraints are considered in the MPC design. In addition, Ψf in eq 56 is chosen to be the solution to the discrete-time Riccati equation in eq 45. The responses are shown in Figures 18 and 19. The responses give an ISE of 3.21 over 60 sampling periods, with all the constraints satisfied. 6. Conclusion A new model predictive control (MPC) for Hammerstein systems, based on the numerical inverse, was developed for process models with multivariable nonlinearities. The most attractive feature of this method is the ability to handle input constraints with the presence of coupled input nonlinearities. An iterative numerical inversion is proposed to address the

Figure 18. Controller output of Example 2 (using the proposed method).

nonlinearities, together with a condition that guarantees the convergence of the iteration. The dynamical control performance is achieved using classical linear MPC, after the numerical inverse cancels the effect of the input nonlinearity of the Hammerstein process. Approximated linear constraints are

178


where

[] []

x(1) x(2) X) · ·· x(N) xs x Xs ) ·s ·· xs

[ ]

u(0) u(1) U) ·· · u(M - 1)

[]

A A2 Λ) l AN

Figure 19. Output response of Example 2 (using the proposed method).

[

obtained using polytopic descriptions, which keeps the entire problem linear. The proposed control design procedure seems complex to a certain extent, especially the transformation of constraints. However, the complexity can be hidden from the users, because the solution procedure can be handled automatically by a computer. The performance of the proposed control design is exemplified by two 2 × 2 Hammerstein processes, which illustrate how both the input and output constraints are handled.

··· 0 0 ··· 0 0 ·· ··· ·· · · M-1 M-2 B Γ ) A B A B ··· AB AMB AM-1B ··· A2B AB ·· ·· ·· ··· ·· · · · · N-1 N-2 N-M A B A B ··· ··· A B

Appendix A. Linear MPC Algorithm25 Referring to Figure 3, assume that the linear subsystem H is described by eq 6 with D ) 0 (most physical systems are strictly proper). The error,

e(k) ) y(k) + d(k) - ys

B AB ·· ·

0 B ·· ·

Xs )

(55)

is to be regulated to zero. Given the current on-line estimated state measurement, x(0), minimize the finite horizon performance index:

∑

AN-Mxs

Then, let

N-1

J0 ) [x(N) - xs]TΨf[x(N) - xs] +

[] xs xs ·· · xs Axs ·· ·

]

eT(k)Ψe(k) +

k)0 M-1

∑ [u(k) - us]TΦ[u(k) - us]

(56)

k)0

Ψ h ) diag [CTΨC, ..., CTΨC, Ψf]

(60)

Φ h ) diag [Φ, ..., Φ]

(61)

Us ) [us

T

where Ψ g0, Φ g0, Ψf g0, M (the control horizon) eN (the prediction horizon). Assume d(k) is the output disturbance and de(k) ) d(k) - ys. Let dhe, us, and xs denote the corresponding steady-state values of de, u, and x, respectively. Then,

us ) -[C(I - A)-1 B]-1dhe

(57)

xs ) (I - A)-1Bus

(58)

T

us

‚‚‚ us ]

T T

(62)

D ˜ ) [de(0)T -dheT de(1)T -dheT ‚‚‚ de(N - 1)T -dheT 0 ]T (63) Z ) diag [CTΨ, CTΨ, ..., CTΨ]

(64)

Hence, eq 56 can be expressed as

J0 ) eT(0)Ψe(0) + (X - Xs)TΨ h (X - Xs) +

Assume that

h (U - Us) + 2(X - Xs)T ZD ˜ + (U - Us)TΦ

u(k) ) us ∀ k gM

D ˜ T diag [Ψ, ..., Ψ] (65) ) Jh0 + UTWU + 2UTV

Therefore, the states in eq 6 can be written as

X - Xs ) ΓU + Λx(0) - X hs

(59)

where Jh0 is independent of U, and

(66)


W ) Γ TΨ hΓ+Φ h V ) Γ TΨ h [Λx(0) - X h s] - Φ h Us + ΓTZD ˜ For the given constraints in the form Lu(k) eb for k ) 0, 1, ..., M - 1 and ymin e y(k) e ymax for all components and k ) 1, 2, ..., N - 1, the optimal solution can be numerically computed via the following static-optimization problem:

UOPT ) arg min UTWU + 2UTV

(67)

U L˜ Ueb˜

where

[ ] [ ]

Du L˜ ) Dy -Dy

(68a)

bu b˜ ) by,max by,min

(68b)

Du is the following dnM × mM matrix:

Du ) diag[L,...,L]

(69)

[

bu ) [bT ... bT ]T

Dy is the following m(N - 1) × mM:

[

0 ···

···

0 ·· · ··· 0 · · · CAB CB ··· ··· CAB ·· · N-M N-M-1 · · · CA B CA B

by,max and by,min are given by the following vectors:

by,max )

by,min )

[

]

ymax - CAx0 - d1 l N-M-2

ymax - CA

N-1

x0 - dN-1 -

∑ i)0

CAiBus

-ymin + CAx0 + d1 l N-M-2

-ymin + CAN-1x0 + dN-1 +

• Create a simplex of m + 1 points (convex hull of m + 1 affinely independent points) • For each facet F: For each unassigned point p*, If p* is above F Assign p* to F’s outside set • For each facet F with a non-empty outside set: Select the furthest point p* of F’s outside set Initialize the visible set V to F For all unvisited neighbors N of facets in V: If p* is above N, Add N to V The boundary of V is the set of horizon ridges H For each ridge R in H: Create a new facet from R and p* Link the new facet to its neighbors For each new facet F′: For each unassigned point p* in an outside set of a facet in V: If p* is above F′, Assign p* to F′’s outside set Delete the facets in V

Acknowledgment

bu is the following mM × 1 matrix:

CB CAB ·· · M-1 B CA Dy ) M CA B ·· · N-2 CA B

for each point. If a point is above multiple new facets, one of the new facets is selected. If it is below all of the new facets, the point is inside the convex hull and can be discarded. Partitioning also records the furthest point of each outside set. An outline of the algorithm is given below:

∑ i)0

CAiBus

(70)

The first author of this paper (K.H.C.) would like to acknowledge the financial support from the Australian Postgraduate Awards. He would also like to thank Dr. Osvaldo J. Rojas for clarifying a few issues regarding MPC design. Literature Cited

(71)

] ]

(72a)

(72b)

This optimization problem can be solved by standard quadratic programming (QP) algorithms in MATLAB. B. Quickhull Algorithm.22 Quickhull selects a nondegenerate set of points (there are no m + 1 given points that are on a common hyperplane) for the initial simplex. If possible, it selects points with either a maximum or minimum coordinate. After initialization, Quickhull assigns each unprocessed point to an outside set. By definition, the corresponding facet is visible from the point. When Quickhull creates a cone of new facets, it builds new outside sets from the outside sets of the visible facets. This process, which is called partitioning, locates a visible new facet

(1) Greblicki, W. Recursive identification of continuous-time Hammerstein systems. Int. J. Syst. Sci. 2002, 33, 969-977. (2) Tian, Y.-C.; Zhao, F.; Bisowarno, B. H.; Tade, M. O. Pattern-based predictive control for ETBE reactive distillation. J. Process Control 2003, 13 (1), 57-67. (3) Patwardhan, R. S.; Lakshminarayanan, S.; Shah, S. L. Constrained nonlinear MPC using Hammerstein and Wiener models: PLS framework. AIChE J. 1998, 44, 1611-1622. (4) Fruzzetti, K. P.; Palazolu, A.; McDonald, K. A. Nonlinear model predictive control using Hammerstein models. J. Process Control 1997, 7, 31-41. (5) Abonyi, J.; Babusˇka, R.; Botto, M. A.; Szeifert, F.; Nagy, L. Identification and control of nonlinear systems using fuzzy Hammerstein models. Ind. Eng. Chem. Res. 2000, 39, 4302-4314. (6) Bloemen, H. H. J.; Van den Boom, T. J. J.; Verbruggen, H. B. Modelbased predictive control for Hammerstein-Wiener systems. Int. J. Control 2001, 74, 482-495. (7) Haddad, W.; Chellaboina, V. Nonlinear control of Hammerstein systems with passive nonlinear dynamics. IEEE Trans. Autom. Control 2001, 46, 1630-1634. (8) Chan, K. H.; Bao, J.; Whiten, W. J. A new approach to control of MIMO processes with static nonlinearities using an extended IMC framework. Comput. Chem. Eng. 2005, 30, 329-342. (9) Msahli, F.; Abdennour, R. B.; Ksouri, M. Nonlinear model-based predictive control using a generalised Hammerstein model and its application to a semi-batch reactor. Int. J. AdV. Manuf. Technol. 2002, 20, 844-852. (10) Ni, X.; Verhaegen, M.; Krijgsman, A. J.; Verbruggen, H. B. A new method for identification and control of nonlinear dynamic systems. Eng. Appl. Artif. Intell. 1996, 9, 231-243. (11) Cervantes, A. L.; Agamennoni, O. E.; Figueroa, J. L. A nonlinear model predictive control system based on Wiener piecewise linear models. J. Process Control 2003, 13, 655-666. (12) Huang, H.; Riggs, J. B. Comparison of PI and MPC for control of a gas recovery unit. J. Process Control 2002, 12 (1), 163-173. (13) Jia, C.; Rohani, S.; Jutan, A. FCC unit modeling, identification and model predictive control, a simulation study. Chem. Eng. Process. 2003, 42 (4), 311-325.

180


(14) Sepulchre, R.; Jankoviaˆc, M.; Kokotoviaˆc, P. V. Constructive nonlinear control. In Communications and Control Engineering; Springer: London, New York, 1997. (15) Lin, W.; Byrnes, C. I. Passivity and absolute stabilization of a class of discrete-time nonlinear systems. Automatica 1995, 31 (2), 263-267. (16) Hitz, L.; Anderson, B. D. O. Discrete positive-real functions and their application to system stability. Proc. IEE 1969, 116, 153-155. (17) Xiao, C.; Hill, D. J. Generalization and new proof of the discretetime positive real lemma and bounded real lemma. IEEE Trans. Circuits Syst. I: Fundam. Theory Appl. 1999, 46, 740-743. (18) Khalil, H. K. Nonlinear Systems, 3rd Edition; Prentice Hall: Upper Saddle River, NJ, 2002. (19) Desoer, C. A.; Vidyasagar, M. Feedback Systems: Input-Output Properties; Academic Press: New York, 1975. (20) Su, S. W.; Bao, J.; Lee, P. L. Analysis of decentralized integral controllability for nonlinear systems. Comput. Chem. Eng. 2004, 28, 17811787. (21) Schrijver, A. Theory of Linear and Integer Programming; Wiley: New York, 1986.

(22) Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. T. The Quickhull algorithm for convex hulls. ACM Trans. Mathemat. Software 1996, 22, 469483. (23) Gru¨nbaum, B. Measure of symmetry for convex sets. In Proceedings of the 7th Symposium in Pure Mathematics of the American Mathematical Society, Symposium on ConVexity, (Seattle, WA); American Mathematical Society (AMS): Providence, RI, 1961; pp 233-270. (24) Kvasnica, M.; Grieder, P.; Baoti&cacute, M. “Multi-Parametric Toolbox (MPT),’’ 2004. (25) Goodwin, G. C.; Graebe, S. F.; Salgado, M. E. Control System Design, International Edition; Prentice Hall: Upper Saddle River, NJ, 2001.

ReceiVed for reView July 14, 2006 ReVised manuscript receiVed October 9, 2006 Accepted October 11, 2006 IE0609113

Model Predictive Control of Hammerstein Systems with Multivariable

Recommend Documents