Ind. Eng. Chem. Res. 1996, 35, 4281-4290
4281
GENERAL RESEARCH Improved Techniques for the Development of Quadratic Perturbation Models Sachin C. Patwardhan†,‡ and K. P. Madhavan*,§ Systems and Control Engineering and Department of Chemical Engineering, Indian Institute of Technology, Powai, Bombay 400 076, India
In this paper, the quadratic model developed by Patwardhan and Madhavan (1993) using the regular perturbation approach has been further refined to produce a better estimate of secondorder effects. As an alternate to the regular perturbation model, a structurally better form of the quadratic model has also been derived that does not require any of the simplifying assumptions that were used in our previous work (Patwardhan and Madhavan, 1993). The coefficients of this model are computed using the solutions of the first- and second-order sensitivity equations. A recursive form of quadratic perturbation models is also developed for the systems with a nonlinear output map. The computational advantages of the proposed quadratic models have been presented. The improved prediction capabilities of the quadratic models obtained using the proposed modifications have been demonstrated using a simulation example. The analysis of the model structure and the simulation results indicate that the quadratic model developed using the sensitivity equation approach is a better choice for approximating the local behavior of a highly nonlinear system and, consequently, for the development of an nonlinear model predictive control scheme. 1. Introduction The nonlinear model predictive control (NMPC) schemes proposed in the literature essentially differ in their approach to the handling of nonlinear ODE constraints describing the plant behavior. Three basic approaches of handling ODE constraints can be identified (Bequette, 1991): (A) direct use of ODE solvers; (B) use of orthogonal collocation; and (C) successive local linearization of the nonlinear ODE. The NMPC strategies that employ explicit ODE solvers require extensive computational effort because of the iterative nature of control law calculations. Orthogonal collocation techniques have the advantage of expressing the dynamic model in an algebraic equation framework so that optimal control laws can be calculated by efficient optimization algorithms. However, the nonlinearity problem remains so that iterative calculations are required for control law development. Hence, direct use of ODE solvers or collocation techniques can pose computational time constraints if realtime implementation is contemplated. On the other hand, the successive linearization based approach requires considerably less computations and appears to be a more attractive option if the real-time implementation is contemplated. However, a locally linearized model may not give good predictions over a long prediction horizon for a strongly nonlinear system. Consequently, an NMPC strategy that uses such prediction models may lead to a suboptimal closed-loop perfor* Author to whom correspondence should be addressed. Email:
[email protected]. † Current address: Chemical Engineering, Indian Institute of Technology, Madras 600 036, India. ‡ Systems and Control Engineering. § Department of Chemical Engineering.
S0888-5885(95)00206-5 CCC: $12.00
mance unless the linear model coefficients are recomputed along the projected trajectory as required by the multistep pseudo-Newton strategy proposed by Li and Biegler (1989). Also, a linear approximation based NMPC strategy may run into difficulties if the system being controlled exhibits a change in the sign of the steady-state gain in the neighborhood of the desired operating point (Patwardhan and Madhavan, 1993). Thus, with the motivation of evolving a prediction model that captures the local nonlinear behavior reasonably accurately and calls for less computational effort, a quadratic perturbation model was developed by Patwardhan and Madhavan (1993). The proposed quadratic model could account for nonlinear behavior along the prediction horizon without requiring recomputation of model coefficients. Also, the quadratic model could be used recursively over the prediction horizon for future trajectory prediction which simplified the optimization calculations at each sampling instant. The application of the Morse lemma revealed that the proposed quadratic approximation has a distinct advantage over a linear approximation in describing the linear local dynamic behavior of a nonlinear system characterized by extremum (maximum or minimum) in their output(s) and attendant change in the sign of the steady-state gain(s). The effectiveness of the NMPC algorithm that uses the quadratic prediction model was demonstrated by simulating two benchmark control problems. In both cases, the objective was to control the system (CSTR/fermenter) at its optimum operating point where a change in the sign of the steady-state gain(s) is observed. The quadratic perturbation models developed using the approach proposed by Patwardhan and Madhavan (1993) had the following drawbacks: 1. Model Structure. The use of regular perturba© 1996 American Chemical Society
4282 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
tion approach for arriving at an approximate form of the quadratic model resulted in a model structure in which the first- and second-order perturbation effects were considered separately. The second-order effects were approximated using the first-order perturbation terms. 2. Average Value Approximation. An average value of the first-order perturbation terms over a sampling interval was used for computing the model coefficients of the second-order terms. This, in general, gave rise to high approximation errors. 3. Output Map. It was assumed that the output map is linear. Further investigations carried out with the motivation of reducing the approximation errors arising from these simplifying assumptions have resulted in the following new developments (Patwardhan, 1994): 1. Improved Regular Perturbation Approach. An improved method for computing second-order coefficients that does not require the average value approximation to be used while developing a quadratic model using the regular perturbation approach. 2. Model Structure. A better recursive form of the quadratic model that adopts the more rigorous method of computing the coefficients of the quadratic model using the solutions of sensitivity equations. 3. A recursive form of the quadratic prediction model has also been evolved for the systems that are characterized by a nonlinear output map. The aim of the present paper is to present the theoretical basis behind the development of quadratic models of greater rigor and demonstrate the improved prediction capabilities of these models on some strongly nonlinear reactor systems. 2. Mathematical Preliminaries The systems under consideration are general “multipleinput multiple-output” (MIMO) lumped parameter processes modeled by a vector differential equation represented as
dX/dt ) F(X,U(t)), X(t0) ) X0
(1)
where X(t) ∈ Rn represent the system states, U(t) ∈ Rm for t ∈ [t0, ∞) represent the manipulated inputs, and F represents an n-dimensional function vector. The system outputs Y ∈ Rr are given by the relation
Y(t) ) G(X(t))
and
Xk+1 ) χ(tk+T;Xk,Uk) ) χk ) χ(T;Xk,Uk) ) χ(Xk,Uk) (3) 3. Quadratic Perturbation Models An explicit functional relationship of the form (3) can be obtained only when the system of ODEs representing the plant can be integrated analytically, and this is seldom possible for a general nonlinear system. Thus, it is necessary to develop a model that approximates the nonlinear operator χ at least locally. In the present study, it is desired to develop quadratic perturbation models to approximate the local behavior of the operator χ in the neighborhood of a given point, say (X h, U h ). Regarding the choice of the basis point (X h, U h ), the following two possibilities are considered: 1. Successive Taylor Expansion. By this approach, at the k sampling instant, the basis point (X h, U h ) is chosen as X h ) Xk, U h ) Uk-1, i.e., the current operating point. Effectively, the model coefficient matrices and bilinear matrices are recomputed at every sampling instant. 2. Fixed Coefficient Model Obtained at a Steady State. By this approach, the basis point (X h, U h ) is chosen as a fixed steady-state point of eq 3, say X h ) Xs and U h ) Us, preferably the steady state corresponding to the setpoint. The model coefficient matrices and bilinear matrices are computed only once, and the resulting perturbation model is used recursively. In the development that follows, the second-order derivatives of operator F or χ appear frequently in the equations. These are three-dimensional arrays and are highlighted by inclusion in curly brackets. The notations used for representing these second-order operator derivatives (also referred to as bilinear matrices in the rest of the text) and their various operations with matrices and vectors are adopted from Economou (1985) (also see Appendix I). These notations provide an alternate and more convenient way of representing the quadratic terms than the notations used in our previous work (Patwardhan and Madhavan, 1993). 3.1. Modified Regular Perturbation Approach. The nonlinear function vector F in eq 1 can be expanded as a multivariable Taylor series in the neighborhood of point (X h, U h ) as
(2)
where G represents an r-dimensional function vector. In the present study, it is assumed that unique solutions of the system of ODEs (2) exist for every U(t) ∈ Rm, t ∈ [t0, ∞), and these solutions are continuously differentiable with respect to all its arguments. Also, it is assumed that the manipulated inputs are piecewise constant functions. Adapting notation from Economou (1985), let the state transition function χ(t;Xk,Uk) denote the solution of eq 1 over a sampling interval tk e t < tk + T with the initial condition
χ(tk;Xk,Uk) ) Xk
dropped from the parameter list, resulting in the following convention:
U(t) ) Uk for tk e t < tk + T
where T denotes the sampling time. Here, letter k is used as a superscript to mark discrete time intervals. Since eq 1 does not involve time explicitly, the term t is
d(X h + x(t)) ) hf + J h xx(t) + J h Uu(t) + dt 1 1 {H h }(x(t),x(t)) + {H h }(x(t),u(t)) + 2 XX 2 UX 1 1 {H h }(u(t),x(t)) + {H h }(u(t),u(t)) + ... 2 XU 2 UU where
x(t) ) X(t) - X h
and
u(t) ) U(t) - U h k k for t ∈ [t , t + T) (4)
h (initial condition) x(tk) ) Xk - X h ,U h )), hf ) F h (X h ,U h ), J h X ) (F(1) X (X J h U ) (F(1) h ,U h )), U (X
{H h XX} ) {F(2) h ,U h )} XX(X
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4283
{H h XU} ) {F(2) h ,U h )}, XU(X
{H h UX} ) {F(2) h ,U h )}, UX(X (2) {H h UU} ) {FUU (X h ,U h )}
Here, the superscript indicates the order of derivative of operator F and the bar indicates that the derivatives have been computed at the point (X h, U h ). As the inputs are assumed to be piecewise constant, we have
u(t) ) uk ) Uk - U h
for t ∈ [tk, tk + T)
k
x1(t ) )
{Φ h 0} )
∫0TΦ(T - τ)•{Hh XX}oΨ(τ)•Ψ(τ) dτ
(10)
∃ ) (1/2){Φ h 1} is (b) Φ h 1}(fh) is a (n × n) matrix where {Φ solution of the integral equation
{Φ h 1} )
The use of regular perturbation technique results in the following set of coupled differential equations representing first- and second-order perturbation effects (see Patwardhan and Madhavan, 1993)
dx1 ) hf + J h Xx1 + J h Uuk, dt
where (a) xj02 ) (1/2){Φ h 0}(fh,fh) is a (n × 1) vector where {Φ h 0} is the solution of the integral equation
∫0TΦ(T - τ)•{Hh XX}oΨ(τ)•Φ(τ) dτ
(11)
(c) Γ∃ ) (1/2){Γ h 1}(fh) is a (n × m) matrix where {Γ h 1} is the solution of the integral equation
∫0TΦ(T - τ)•{{Hh XX}oΨ(τ)•Γ(τ) +
{Γ h 1} ) 2
xk1
(5)
{H h UX}oΨ(τ)} dτ (12)
dx2 1 1 )J h Xx2 + {H h }(x (t),x1(t)) + {H h }(x (t),uk) + dt 2 XX 1 2 UX 1 1 1 {H h }(uk,x1(t)) + {H h }(uk,uk) x2(tk) ) xk2 (6) 2 XU 2 UU
(d) Bilinear matrix {Φ h X} is the solution of the integral equation
Equation 5 can be integrated over a sampling interval to obtain a discrete linear perturbation model as follows:
(e) Bilinear matrix {Φ h U} is the solution of the integral equation
) xj01 + Φ h xk1 + Γ h uk xk+1 1
(7)
{Φ h X} )
{Φ h U} )
∫0TΦ(T - τ)•{Hh XX}oΦ(τ)•Φ(τ) dτ
(13)
∫0TΦ(T - τ)•{{Hh XX}oΓ(τ)•Φ(τ) +
{H h XU}•Φ(τ)} dτ (14)
where
Φ h ) exp(J h XT), Ψ h )
Γ h)
∫0Texp(Jh Xτ) Jh U dτ,
∫0 exp(Jh Xτ) dτ, T
and
xj01 ) Ψ h hf
Note that eq 6 is linear in x2(t) but has a nonhomogeneous term depending on the previous term x1(t) as well as ({H h UU}(uk,uk)). Integrating eq 6 over a sampling interval, we obtain
xk+1 )Φ h xk2 + 2
∫0TΦ(T - τ)[21{Hh XX}(x1(τ),x1(τ)) +
1 1 {H h }(x (τ),uk) + {H h }(uk,x1(τ)) + 2 UX 1 2 XU 1 {H h }(uk,uk) dτ (8) 2 UU
]
where x1(τ) for τ ) t - tk is defined as
x1(τ) ) x01(τ) + Φ(τ) Φ(τ) ) exp(J h Xτ),
xk1
k
+ Γ(τ) u
Γ(τ) )
(f) Bilinear matrix {Γ h X} is the solution of the integral equation
{Γ h X} )
∫0TΦ(T - τ)•{{Hh XX}oΦ(τ)•Γ(τ) +
{H h UX}oΦ(τ)} dτ
(g) Bilinear matrix {Γ h U} is the solution of the integral equation
{Γ h U} )
∫0TΦ(T - τ)•{{Hh XX}oΓ(τ)•Γ(τ) +
{H h XU}•Γ(τ) {H h UX}oΓ(τ) + {H h UU}} dτ (15)
Thus, defining
xk ) xk1 + xk2 (0 e τ < T)
∫0 exp(Jh Xq) Jh u dq τ
In order to simplify eq 8 further, it is necessary to introduce certain algebraic operations (such as right dot product, left dot product, and circle product) involving matrices and bilinear matrices. The definitions of these operations along with some algebraic identities required for achieving the desired simplification are given in Appendix I. Thus, using the bilinearity property of second-order derivatives and Results I.1-I.2 proved in Appendix I, eq 8 can be rearranged as follows:
1 ∃xk + Φ xk+1 h }(xk,xk) + ) xj02 + Φ h xk2 + Γ∃uk + {Φ 2 1 2 X 1 1 1 1 1 {Φ h }(uk,xk1) + {Γ h }(xk,uk) + {Γ h }(uk,uk) (9) 2 U 2 X 1 2 U
and
xj0 ) xj01 + xj02
(16)
the regular perturbation quadratic model can be represented as
xk1 ) xj01 + Φ h xk1 + Γ h uk ∃xk + 1{Φ h }(xk,xk) + xk+1 ) xj0 + Φ h xk + (Γ∃ + Γ h )uk + Φ 1 2 X 1 1 1 h }(uk,uk) (17) {Φ h U}(uk,xk1) + {Γ 2 U When employed in an NMPC formulation, these equations can be used recursively over the prediction horizon for the future trajectory prediction. 3.2. Sensitivity Equations Approach. The multivariable Taylor series expansion of the discrete state transition model given by eq 3 in the neighborhood of a point (X h, U h ), maintaining only up to second-order terms, can be written as
4284 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
1 h }(xk,xk) + Xk+1 ) χj0 + Φ h xk + Γ h uk + {Φ 2 X 1 1 1 {Φ h }(uk,xk) + {Γ h }(xk,uk) + {Γ h }(uk,uk) (18) 2 U 2 X 2 U where k
h x )X -X def
def
{Φ h X} )
def
{Γ h U} )
k
and def
χj0 ) χ(T,X h ,U h ),
{ } { }
∂2χk , ∂(Xk)2
Φ h )
k-1 J h U(t) ) (F(1) )) U (χN(t),U
k-1 {H h XX(t)} ) {F(2) )}, X (χN(t),U
∂χk , ∂Xk
def
{Γ h X} )
def
Γ h )
k-1 {H h UU(t)} ) {F(2) )} U (χN(t),U
∂χk ∂Uk
{ } { } { } { }
∂ ∂χk ∂2χk ) ∂Uk ∂Xk ∂Uk∂Xk
def
∂χ , ∂(Uk)2
k-1 {H h UX(t)} ) {F(2) )}, UX(χN(t),U
u )U -U h
{Φ h U} )
2 k
k
k
∂ ∂χ ∂χ ) k k ∂X ∂U ∂Xk∂Uk
xj0 ) χj0 - X h
(19)
and using the identity
the following quadratic perturbation model is obtained:
1 h }(xk,xk) + ) xj0 + Φ hx + Γ h u + {Φ 2 X 1 h }(uk,uk) (20) {Φ h U}(uk,xk) + {Γ 2 U
Φ(tk) ) I
(23)
Input sensitivity matrix Γ h is the solution at t ) tk+1 of the initial value problem
∂Γ(t) ) JX(t) Γ(t) + JU(t), ∂t
Γ(tk) ) [0 h]
(24)
where [0 h ] represents a (n × m) null matrix. def
The bilinear matrix {Φ h X} ) {∂Φ/∂Xk} is the solution at t ) tk+1 of the initial value problem
{ΦX(tk)} ) {0 h}
(25)
k
Note that, in general, the vector xj0 * 0h except at the steady-state conditions and depends on the choice of (X h, U h ). When the successive Taylor expansion approach is used, the term χj0 can be interpreted as the state variable vector obtained by choosing uk ) 0h , i.e., holding Uk-1 constant for tk e t < tk + T. Thus, χj0 can be obtained by integrating eq 1 as
χj0 ) Xk +
∂Φ(t) ) JX(t) Φ(t), ∂t
∂{ΦX(t)} ) JX(t)•{ΦX(t)} + {HXX(t)}oΦ(t)•Φ(t) ∂t
h X}(xk,uk) {Φ h U}(uk,xk) ) {Γ
k
the state transition matrix Φ h is the solution at t ) tk+1 of the initial value problem
2 k
represent the first- and second-order derivatives of the operator χk evaluated at (X h, U h ). Subtracting Xk from both sides of eq 18 and defining xk as
x
k-1 J h X(t) ) (F(1) )), X (χN(t),U
k-1 {H h XU(t)} ) {F(2) )}, XU(χN(t),U k
k+1
Defining
∫tt
K+!
K
F((X(t),Uk-1) dt
(21)
and xj0 is recomputed at each sampling instant. When the perturbation model is developed in the neighborhood of a steady state, since at a steady state F(Xs,Us) ) 0 h, we have χj0 ) Xs and xj0 ) 0h . The first-order derivatives of operator χk, i.e., Φ h and Γ h , can be computed by integrating the classical (firstorder) sensitivity equations over a sampling interval. The differential equations that can be used for evaluation of second-order derivatives of operator χk have been derived by Economou (1985). These differential equations, referred to as second-order sensitivity equations in the present text, can be integrated over a sampling interval to obtain the bilinear coefficient matrices of the quadratic model (20). The original equations are adopted here with the following modification: The first- and second-order derivatives of function vector F appearing in these sensitivity equations are computed along a nominal trajectory χN(t) over interval tk e t < tk + T defined as
χjN ) χ(t;X h ,U h ) ) X(tk) +
∫ttF((X(t),Uh ) dt k
(22)
def
The bilinear matrix {Φ h U} ) {∂Φ/∂Uk} is the solution at t ) tk+1 of the initial value problem
∂{ΦU(t)} ) JX(t)•{ΦU(t)} + {HXX(t)}oΓ(t)•Φ(t) + ∂t {HXU(t)}•Φ(t) {ΦU(tk)} ) {0 h}
(26)
def
The bilinear matrix {ΓU} ) {∂Φ/∂Uk} is the solution at t ) tk+1 of the initial value problem
∂{ΓU(t)} ) JX(t)•{ΓU(t)} + {HXX(t)}oΓ(t)•Γ(t) + ∂t {HXU(t)}•Γ(t) + {HUX(t)}oΓ(t) + {HUU(t)} (27) {ΓU(tk)} ) {0 h} When the repeated Taylor expansion approach is adopted to develop a quadratic model, the coupled differential equations (23)-(27) have to be solved simultaneously with initial value problem
dX ) F(X(t),Uk-1), dt
X(tk) ) Xk
(28)
over a sampling interval tk e t < tk + T. This can be achieved using any standard nonlinear integration technique. The solution procedure is considerably simplified when the quadratic perturbation model is developed at a steady state, say (Xs, Us). In this case the nominal trajectory reduces to
χjN(t) ) χ(t;Xs,Us) ) Xs
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4285
and, consequently, the derivatives of operator F appearing in eqs 23-27 have to be computed only at one point, i.e., (Xs, Us). 3.3. Output Prediction. When the state-output map G in eq 2 is nonlinear, the quadratic perturbation model has to be modified if it is desired to predict the future output trajectory by recursive use of the prediction model. Consider the single step output prediction model given as
Yk+1 ) G(Xk+1) ) G(χ(Xk,Uk))
(29)
Using Results I.3 and I.4 stated in Appendix I, the Taylor series expansion of the output vector Yk+1 in the neighborhood of a point (X h, U h ), maintaining only up to second-order terms, can be written as
1 yk+1 = yk0 + C h X}oΦ hΦ h xk + C hΓ h xk + {{C h •Φ h + 2 1 h X}oΓ h •Φ h+ C h •{Γ h X}}(xk,uk) + C h •{Φ h X}}(xk,xk) + {{C 2 1 1 {{C h X}oΦ h X}oΓ h •Γ h+C h •{Φ h U}}(uk,xk) + {{C h •Γ h+ 2 2 C h •{Γ h U}}(uk,uk) (30) where
h; yk+1 ) Yk+1 - Y def
|
∂G C h ) k ; ∂χ χK)χj0
yk0 ) G(χj0) - Y h;
Y h ) G(X h)
{ | }
∂2 G {C h X} ) ∂(χk)2 χK)χj0 def
Though eq 30 gives the most general form of the quadratic expansion, in practice, the computational complexity can be reduced if the system under consideration belongs to either of the following two categories. Case A: System with Nonlinear State Dynamics and Linear Output Map.
Yk+1 ) CXk+1
and
{C h X} ) {0 h}
where C is a (n × r) output matrix. Subsequently, eq 30 reduces to
1 h •{Φ h X}}(xk,xk) + h xk + CΓ h xk + {C yk+1 = yk0 + CΦ 2 1 1 {C h •{Φ h U}}(uk,xk) + {C h •{Γ h X}}(xk,uk) + 2 2 1 {C h •{Γ h U}}(uk,uk) (31) 2
assumptions that the models have been developed at a steady state X h ) Xs and U h ) Us of eq 3. In this analysis, the notations used for representing the quadratic model operator are as follows: a. Sensitivity Equation Approach.
xk+1 ) ξ(xk,uk) 1 ξ(xk,uk) ) Φsxk + Γsuk + {ΦsX}(xk,xk) + 2 1 {ΦsU}(uk,xk) + {ΓsU}(uk,uk) (33) 2 b. Regular Perturbation Approach.
) Φsxk1 + Γsuk xk+1 1 xk+1 ) ξ(xk1,xk,uk)
1 ξ(xk1,xk,uk) ) Φsxk + Γsuk + {ΦsX}(xk1,xk1) + 2 1 {ΦsU}(uk,xk1) + {ΓsU}(uk,uk) 2 4.1. Theoretical Basis for Studying Open-Loop Stability. Consider the open-loop operator O: Rn f Rn defined for a constant input u* o
xk 98 xk+1 ) ξ(xk,u*)
1 h X}oΦ h xk + CΓ h xk + {{C h •Φ h }(xk,xk) + yk+1 = yk0 + CΦ 2 1 1 {{C h X}oΓ h X}oΦ h •Φ h }(xk,uk) + {{C h •Γ h }(uk,xk) + 2 2 1 {{C h X}oΓ h •Γ h }(uk,uk) (32) 2 4. Open-Loop Stability of Quadratic Perturbation Models In this section, the open-loop stability analysis of the quadratic perturbation models is presented under the
(35)
and an equilibrium state x*. Theorem 1 (Economou, 1985). Consider the openloop operator given by eq 35 and an equilibrium state x*. If
|
|
∂ξ(x,u*) eR1 ∂x
for every x ∈ Ω
where F denotes the spectral radius, then the region Ω does not contain any equilibrium state of operator ξ. 4.2. Stability Analysis of Sensitivity Equation Based Quadratic Model. Using Liapunov’s first (indirect) method, it can be shown that the equilibrium state (x*, u*) of the open-loop quadratic operator in eq 33 is locally asymptotically stable if the following condition holds
F(Φs + {ΦsU}u* + {ΦsX}x*) < 1 where F indicates the spectral radius. Theorem 1 can be used for the characterization of the region of attraction around x*. Thus, if there exists a nonzero neighborhood of ϑ j (x*,r) of x* such that
4286 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
| |
∂xk+1 eR 1
for every x ∈ Ω
where F is the spectral radius of the operator in the brackets. 4.3. Stability Analysis of Regular Perturbation Based Quadratic Model. Using Liapunov’s first (indirect) method, it can be shown that the equilibrium state (x*1, x*, u*) of an open-loop quadratic operator in eq 34 is asymptotically stable if the following condition holds:
(
∂(xk1,xk)
| )
[
]
k+1 ∂(xk+1 ) 1 ,x
F
