Interior Point Solution of Multilevel Quadratic Programming Problems

Ind. Eng. Chem. Res. 2008, 47, 81-91

81

Interior Point Solution of Multilevel Quadratic Programming Problems in Constrained Model Predictive Control Applications R. Baker and C. L. E. Swartz* Department of Chemical Engineering, McMaster UniVersity, 1280 Main Street West, Hamilton, Ontario, Canada L8S 4L7

This paper examines the use of an interior point strategy to solve multilevel optimization problems that arise from the inclusion of the closed-loop response of constrained, linear model predictive control (MPC) within a primary quadratic or linear programming problem. We motivate the formulation through its application to optimizing control problems, although the strategy is applicable to several problem types. The problem is cast as a dynamic optimization problem in which an optimal steady-state operating point is sought, subject to constraints on the closed-loop response of the system under constrained predictive control. Because a quadratic programming (QP) problem must be solved at every sampling period, the resulting problem is multilevel in nature. The formulation approach used in this paper is to include the Karush-Kuhn-Tucker (KKT) conditions that correspond to the MPC quadratic programming subproblems as constraints within a single-level optimization problem. The resulting complementarity constrained optimization problem is shown to be reliably and efficiently solved using an interior point approach. The method is applied to two case studies, and its performance is compared to an alternative mixed-integer programming solution approach. 1. Introduction A key objective in chemical process operation and design is to determine an operating point and design parameters that correspond to the economic optimum, subject to the prevailing operational, safety, and environmental constraints. However, all processes are subject to disturbances, and determination of operating conditions based on steady-state models alone may consequently lead to constraint violations in practice. This has led to the incorporation of dynamics in process operation and design calculations. Narraway and co-workers1,2 recognized that the economic optimum generally lies at the intersection of process constraints and, for linear dynamics, determine the required amount of backoff from constraints for feasible operation for assumed disturbance characteristics. From this observation, an approximation of the economic loss or giveaway is determined. This back-off concept is illustrated in Figure 1. The approach is used to assess candidate process designs1 and control structures2 on an economic basis that takes into account the dynamic performance of the system. The authors use a frequency response approach and assume that controlled variables are perfectly controlled, with the disturbance variation propagated to other variables. An extension to nonlinear dynamic models is described in Narraway et al.,3 where a formulation for the synthesis of a multiloop proportional-integral (PI) control system is developed. The closed-loop dynamics of the system are accounted for through inclusion of the controller equations within a dynamic optimization framework. Figueroa et al.4 considered constraint backoff within a dynamic optimization framework, where a steady-state economic objective is optimized, subject to constraints on the closed-loop response. They computed what they defined as the maximum percentage recovery, which reflects the degree to which the nominal steady-state optimum can be approached for a given control system, subject to constraint feasibility for disturbances * To whom correspondence should be addressed. Tel.: (905) 525 9140. Fax: (905) 521 1350. E-mail: [email protected].

Figure 1. Backoff of the steady-state operating point from constraints.

within a specified set. A sequential solution approach is applied, and the method used to assess the performance of alternative control schemes. Bahri et al.5 apply a similar strategy using, instead, a simultaneous strategy for solution of the dynamic optimization problem. Mohideen et al.6 consider integrated plant and control system design within a dynamic optimization framework, where the decision space includes equipment design parameters, the steady-state operating point, the control structure, and PI controller tuning parameters. Several studies following related approaches have appeared in the literature, accounts of which may be found in the work of van Schijndel and Pistikopoulos.7 Strategies for solving the mixed-integer dynamic optimization problems that arise in such formulations are advanced in the work of Schweiger and Floudas8 and that of Pistikopoulos and Sakizlis.9 Most prior work on the direct inclusion of dynamic effects in process design and/or steady-state operating point calculations considered perfect control,1,2 no control,10 or linear control.4-6,8 Young et al.11 considered actuator saturation effects in optimizing control and showed that failure to account for model discontinuities induced by actuator saturation could result in

10.1021/ie070270r CCC: $40.75 © 2008 American Chemical Society Published on Web 11/29/2007

82 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008

conservative steady-state operating points. Baker and Swartz12 incorporated the actuator saturation formulation within an optimization-based integrated plant and control system design framework and, in another work,13 they proposed an interior point solution approach for such problems. The attractive features of model predictive control (MPC) and its wide acceptance for advanced control in the chemical process industry14 make it a natural choice for inclusion within dynamic optimization problems of the type previously discussed. Loeblein and Perkins15 considered unconstrained MPC in the structural design of real-time optimization systems. Brengel and Seider16 considered nonlinear predictive control in an optimization-based integrated design and control formulation. The MPC calculation subproblems are replaced by their Karush-KuhnTucker (KKT) optimality conditions and a homotopy-continuation method was applied to the solution of the overall problem. Sakizlis et al.17 included constrained MPC within a simultaneous design and control framework using a parametric formulation of MPC.18 Steep hyperbolic tangent functions are used to approximate the Boolean algebraic conditions that define the critical regions in the state space where the piecewise affine control functions hold. In this paper, we consider calculation of an optimal steadystate operating point, subject to constraints on the dynamic response of the process under constrained MPC. Because a quadratic programming (QP) problem must be solved at every controller execution, this results in a multilevel optimization problem. We replace the convex MPC QP subproblems by the necessary and sufficient KKT conditions, giving a single-level mathematical program with complementarity constraints (MPCC). We solve this problem using an interior point technique that has been tailored for MPCCs and compare results obtained to those obtained using an alternative mixed-integer quadratic (linear) programming formulation that is applicable if the objective function of the primary problem is quadratic (linear). The paper is organized as follows. Section 2 presents the problem formulation, where we describe the MPC formulation used, as well as the overall multilevel optimization problem. Section 3 focuses on the solution approach, where the KKT equations that correspond to the MPC subproblems are described and the interior point solution approach is discussed. An alternative mixed-integer programming solution strategy, which we utilize in performance comparisons, is also described. Section 4 presents two case studies: the first involves a singleinput/single-output (SISO) system and, the second involves a multi-input/multi-output (MIMO) system that is based on the Shell Standard Control Problem. This is followed by conclusions.

x(k + 1) ) Ax(k) + Bu(k) y(k) ) Cx(k) where x ∈ Rnx is a state vector, u ∈ Rnu is an input vector, and y ∈ Rny is a vector of measured outputs. A ∈ Rnx×nx, B ∈ Rnx×nu, and C ∈ Rny×nx are linear(ized), discrete-time, state-space matrices. The optimization problem to be solved at each sampling period takes the following form: P

min φ )

|yˆ (k + i|k) - r(k + i|k)|Q + ∑ i)1 i

M-1

|∆uˆ (k + i|k)|R ∑ i)0

i

(1)

subject to

xˆ (k + i|k) ) Axˆ (k + i - 1|k) + Buˆ (k + i - 1|k) (for i ) 1, ..., M) (2) xˆ (k + i|k) ) Axˆ (k + i - 1|k) + Buˆ (k + M - 1|k) (for i ) M + 1, ..., P) (3) yˆ (k + i|k) ) Cxˆ (k + i|k) + dˆ (k + i|k) (for i ) 1, ..., P) (4) ∆uˆ (k|k) ) uˆ (k|k) - u(k - 1)

(5)

∆uˆ (k + i|k) ) uˆ (k + i|k) - uˆ (k + i - 1|k) (for i ) 1, ..., M - 1) (6) umin e uˆ (k + i|k) e umax

(for i ) 0, ..., M - 1) (7)

The norms in eq 1 are defined as

|x|Q ) xTQx and yˆ (k + i|k) represents the predicted value of the outputs at time step k + i, based on information available at time step k. Remark 1. We consider here the dynamic matrix control (DMC) disturbance estimation scheme,19,20 where the disturbance estimate is assumed constant over the prediction horizon and computed as the difference between the measured output and that predicted based on information available at the previous time step. The corresponding equations are

dˆ (k|k) ) y(k) - Cxˆ (k|k - 1) 2. Problem Formulation 2.1. Model Predictive Control (MPC). Model predictive control (MPC) is a control algorithm in which an optimal trajectory of future manipulated inputs is computed, based on the predicted response generated by an internal process model. The inputs corresponding to the first sampling period are implemented and the calculation process repeated at each sampling period. The difference between measured and predicted process outputs is used to adjust the model prediction and provides a feedback mechanism. Features of MPC include the capability for direct inclusion of input and output constraints within the control algorithm, and implicit recognition of process interactions through a multivariable model. In this study, we use a discrete-time state-space model of the form

dˆ (k + i|k) ) dˆ (k|k)

(for i ) 1, ..., P)

The predicted states at k may be obtained from

xˆ (k|k) ) xˆ (k|k - 1) ) Axˆ (k - 1|k - 1) + Bu(k - 1) initialized with xˆ (0|0) ) 0 for the system initially at steady state. This is consistent with the DMC formulation as originally proposed19 (with the use of a state-space model, instead of a step response model). However, more sophisticated state estimation schemes may be desirable;20 these are readily accommodated within the proposed framework. Remark 2. We may express the aforementioned optimization problem in terms of the variables ∆uˆ k alone, where

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 83

∆uˆ k ) [∆uˆ (k|k)T, ..., ∆uˆ (k + M - 1|k)T]T

(8)

as shown in Maciejowski.20 Although this results in a quadratic programming problem of smaller dimension, it also involves the computation of Ai, which for large prediction horizons could lead to numerical problems.20 Therefore, we use the direct form of the model equations as stated previously. We do not include hard output constraints within the MPC optimization problem, but we do include them instead in the outer optimization problem, as described in the next section. Hard output constraints within the MPC calculation could lead to closed-loop instability.21 Remark 3. For clarity of exposition, we have used a basic formulation of MPC. However, various enhancements such as soft output constraints, constraints on input changes, a terminal penalty, and terminal constraints, are readily included within the framework presented. 2.2. Multilevel Dynamic Optimization Problem. Our goal is to determine an optimal steady-state operating point, subject to constraints on the steady-state inputs and outputs, as well as on the closed-loop input and output trajectories with the process controlled via constrained MPC. The formulation below considers a single, unmeasured disturbance, but it is readily extended to accommodate disturbance ranges through, for example, a multiperiod formulation;11 this is applied in the second case study. The controller variables are considered as deviations around the steady-state operating point, which is itself included in the optimization decision space. Therefore, we distinguish between bounds that are fixed and those that apply to the deviation variables through use of a tilde (∼) above the variable symbol for the former (y˜ max, y˜ min, etc.). The optimization problem is stated as follows.

MPC execution at time step k, P

min φ )

|yˆ (k + i|k) - r(k + i|k)|Q + ∑ i)1 i

M-1

|∆uˆ (k + i|k)|R ∑ i)0

i

subject to

xˆ (k + i|k) ) Axˆ (k + i - 1|k) + Buˆ (k + i - 1|k) (for i ) 1, ..., M) xˆ (k + i|k) ) Axˆ (k + i - 1|k) + Buˆ (k + M - 1|k) (for i ) M + 1, ..., P) yˆ (k + i|k) ) Cxˆ (k + i|k) + dˆ (k + i|k)

(for i ) 1, ... P)

∆uˆ (k|k) ) uˆ (k|k) - u(k - 1) ∆uˆ (k + i|k) ) uˆ (k + i|k) - uˆ (k + i - 1|k) (for i ) 1, ... M - 1) umin e uˆ (k + i|k) e umax

(for i ) 0, ... M - 1)

where

xˆ (k|k) )

{

0, k)0 Axˆ (k - 1|k - 1) + Bu(k - 1),

k>0

dˆ (k|k) ) y(k) - Cxˆ (k|k) dˆ (k + i|k) ) dˆ (k|k)

(for i ) 1, ..., P)

min Φ(uss, yss) Plant input:

subject to: Steady-state model:

f(uss, yss, dss) ) 0 Bounds on steady-state operating point:

y˜ min e yss e y˜ max u˜ min e uss e u˜ max

u(k) ) uˆ (k|k) Plant output:

x(k + 1) ) Ax(k) + Bu(k) y(k + 1) ) Cx(k + 1) + d(k + 1) (for k ) 0, ..., Ns - 1)

Bounds on closed-loop response:

y˜ min e y(k) + yss e y˜ max Controller bounds as a function of the steady-state operating point:

umin ) u˜ min - uss umax ) u˜ max - uss Initialization of plant state:

x(0) ) 0 Initial plant output:

y(0) ) Cx(0) + d(0)

In the aforementioned equations, Ns is the horizon over which the closed-loop response is computed. We note that the disturbance sequence, {d(k)}, is specified. The disturbance estimates, dˆ (k + i|k), are determined as discussed in Section 2.1. Also, the model used for the evolution of the plant outputs does not need to be the same as that used within the MPC controller; thus, different coefficient matrices are used. We observe that the MPC subproblems required to generate the closed-loop response result in a multilevel optimization problem. The MPC QP subproblems are dependent on the inputs and model outputs generated at the previous step, as well as the steady-state inputs, which define the constraint limits in deviation form, in relation to the fixed, physical limits. The primary (outer) problem constraints involve the plant outputs at the end of each sampling period. Input constraints are not


included at the outer level, because they are enforced by the MPC controller. Multilevel optimization problems share several characteristics of the more commonly treated bilevel optimization problems, which involve only two levels of interaction. Interactions between the inner and outer optimization problems can lead to issues such as multiple optima of the inner problem, nonconvexity, nondifferentiability, and degeneracy.22 The general form of a bilevel problem may be stated as follows:

min F(w,z) w

subject to

H(w,z) ) 0 G(w,z) e 0 min f(w,z) z

subject to

h(w,z) ) 0 g(w,z) e 0 where the uppercase letters refer to the outer problem and the lowercase letters refer to the inner problem. Several solution approaches involve replacement of the inner optimization problem by its KKT conditions.23,24 Consider a QP inner problem of the following general form:

min z

(21)z Hz + g z T

T

subject to

Az ) b zg0

(9)

where z ∈ Rn. The KKT conditions for this QP problem can be written as

Hz + g - ATλ - ν ) 0 Az - b ) 0 ziνi ) 0

(for i ) 1, ..., n) (z,ν) g 0

(10)

where λ and ν are the associated Lagrange multipliers. For convex QPs, the KKT conditions are both necessary and sufficient for global optimality of the inner problem. The transformation of the QP optimization problem to its KKT firstorder optimality conditions introduces complementarity constraints, which may cause computational difficulties. Solution approaches for bilevel programming problems have used different strategies to address these constraints. Fortuny-Amat and McCarl23 considered problems that are quadratic in both the outer and inner problems, and transform the complementarity constraints of the inner problem KKT conditions into mixed-integer constraints through the

introduction of binary variables. This approach is able to yield the global optimum; however, the combinatorial nature of the solution procedure means that the worst-case solution time grows exponentially with problem size. Bard and Moore24 proposed a similar strategy for linear/quadratic bilevel problems, but with a specialized branch-and-bound scheme. Although significantly superior solution times are observed, relative to the previous approach, a combinatorial element to the solution strategy remains. Clark and Westerberg22 proposed two algorithms for finding local solutions of nonlinear bilevel programming problems. One involves relaxation of the complementarity constraints to be less than or equal to a positive constant, with the nonlinear problem solved for a series of decreasing values of the relaxation parameter, whereas the other is an active set strategy that avoids addressing the complementarity constraints directly. Gu¨mu¨s¸ and Floudas25 presented an approach for the global solution of general nonlinear bilevel problems that involve twice-differentiable functions. They applied relaxation of the feasible region by convex underestimation within a branch-and-bound framework. Several alternative solution strategies for bilevel programming problems are described in the work of Gu¨mu¨s¸ and Floudas25 and Bard.26 The inclusion within the main optimization problem of the first-order optimality conditions of the inner optimization problem, with associated complementarity constraints, leads to a mathematical program with complementarity constraints (MPCC). Because of previous numerical performance, solving MPCCs using nonlinear programming solvers had been considered to be numerically unsafe;27,28 however, recent work has demonstrated that it is possible to solve a large class of MPCCs reliably, using nonlinear programming (NLP) solvers with reformulation of the complementarity constraints27 or modifications to the algorithm.29 In particular, SQP algorithms were observed to perform well, while interior point codes fared less favorably.27 In this paper, we apply an interior point solution strategy, but one in which the complementarity constraints are explicitly considered. Soliman30 considered a steady-state constraint back-off problem similar to that treated in the present paper, where an optimal steady-state operating point is sought such that feasibility of the closed-loop response under constrained MPC is maintained. The MPC QP subproblems are replaced by their KKT conditions, and the complementarity constraints that result from the KKT conditions are formulated as mixed-integer linear constraints. The rapid growth in solution time of the resulting mixed-integer quadratic programming problem (MIQP) with an increasing number of integer variables limits the dimension of problems to which the solution technique can be effectively applied. This motivates the alternative interior point strategy described in the next section.

3. Solution Strategy 3.1. Transformation to a Single-Level Complementarity Constrained Problem. In this section, we will derive the equations for the first-order optimality conditions of a linear, discrete-time, state-space constrained model predictive controller. These may be used to replace the inner MPC subproblems in the multilevel problem described in Section 2.2. The Lagrangian of the constrained MPC optimization problem defined by eqs 1-7 is


uˆ (k|k) - umin - νL(k|k) ) 0

P

L)

(yˆ (k + i|k) - r(k + i|k))TQi(yˆ (k + i|k) ∑ i)1

l

M-1

r(k + i|k)) +

(20)

∆uˆ (k + i|k)TRi∆uˆ (k + i|k) ∑ i)0

uˆ (k + M - 1|k) - umin - νL(k + M - 1|k) ) 0

(21)

umax - uˆ (k|k) - νU(k|k) ) 0

(22)

λ1(k + 1|k)T(xˆ (k + 1|k) - Axˆ (k|k) - Buˆ (k|k)) -

l

... - λ1(k + P|k)T(xˆ (k + P|k) - Axˆ (k + P - 1|k) Buˆ (k + M - 1|k)) - λ2(k + 1|k) (yˆ (k + 1|k) Cxˆ (k + 1|k) - dˆ (k + 1|k)) - ... λ2(k + P|k)T(yˆ (k + P|k) - Cxˆ (k + P|k) -

umax - uˆ (k + M - 1|k) - ν (k + M - 1|k) ) 0 (23)

T

U

Complementarity constraints:

NL(k|k)λ4(k|k) ) 0

dˆ (k + P|k)) - λ3(k|k)T(uˆ (k|k) - u(k - 1) ∆uˆ (k|k)) - ... - λ3(k + M - 1|k) (uˆ (k + M - 1|k) uˆ (k + M - 2|k) - ∆uˆ (k + M - 1|k)) λ4(k|k)T(uˆ (k|k) - umin) - ... T

(24)

l N (k + M - 1|k)λ4(k + M - 1|k) ) 0

(25)

NU(k|k)λ5(k|k) ) 0

(26)

L

λ4(k + M - 1|k)T(uˆ (k + M - 1|k) - umin) -

l

λ5(k|k)T(umax - uˆ (k|k)) - ... λ5(k + M - 1|k) (umax - uˆ (k + M - 1|k))

N (k + M - 1|k)λ5(k + M - 1|k) ) 0

(27)

Including slack variables νL and νU for the lower and upper inequality input constraints, respectively, the KKT optimality conditions for the problem are as follows. Gradient of L, with respect to ∆uˆ :

νL(k|k), ..., νL(k + M - 1|k) g 0

(28)

ν (k|k), ..., ν (k + M - 1|k) g 0

(29)

λ4(k|k), ..., λ4(k + M - 1|k) g 0

(30)

2R0∆uˆ (k|k) + λ3(k|k) ) 0 (11)

λ5(k|k), ..., λ5(k + M - 1|k) g 0

(31)

T

U

U

U

l 2RM-1∆uˆ (k + M - 1|k) + λ3(k + M - 1|k) ) 0 (12) Gradient of L, with respect to yˆ :

2Q1yˆ (k + 1|k) - 2Q1r(k + 1|k) - λ2(k + 1|k) ) 0

(13)

l 2QPyˆ (k + P|k) - 2QPr(k + P|k) - λ2(k + P|k) ) 0

(14)

Gradient of L, with respect to xˆ :

-λ1(k + 1|k) + ATλ1(k + 2|k) + CTλ2(k + 1|k) ) 0 (15) l -λ1(k + P - 1|k) + ATλ1(k + P|k) + CTλ2(k + P - 1|k) ) 0 (16) -λ1(k + P|k) + CTλ2(k + P|k) ) 0 (17) Gradient of L, with respect to uˆ :

BTλ1(k + 1|k) - λ3(k|k) + λ3(k + 1|k) - λ4(k|k) + λ5(k|k) ) 0 (18)

where NL ) diag(νL) and NU ) diag(νU). Equations 2-6 are also included in the KKT conditions. By embedding these KKT conditions within the steady-state optimization problem in place of the inner QP subproblems, we obtain a single-level optimization problem. We note, in addition, that, because the Hessian of an MPC QP subproblem is positive semidefinite, the KKT optimality conditions are necessary and sufficient, and the transformed problem is equivalent to the original, multilevel problem. However, the KKT conditions introduce complementarity constraints, which may lead to solution difficulties if not addressed appropriately. We describe next an effective solution strategy for problems of this type. 3.2. Interior Point Solution. Interior point methods have generated considerable interest and activity in mathematical programming theory and practice over the past two decades or so. Initial developments have focused on applications to linear programming, where polynomial worst-case growth in solution time with problem size was established, as opposed to worst-case exponential growth of the simplex method (although in practice, the performance of the simplex method is generally considerably superior to the theoretical worst case). Very efficient interior point codes have also been developed. Interior point methods have subsequently been applied to several other problem types, including nonconvex NLP problems.31-33 We outline briefly the application of a primal-dual interior point approach to the solution of the NLP problem

l B λ1(k + P|k) + ... + B λ1(k + M|k) - λ3(k + M - 1|k) λ4(k + M - 1|k) + λ5(k + M - 1|k) ) 0 (19) T

T

Inequality constraints rewritten with slacks:

min φ(w)

(32a)

h(w) ) 0

(32b)

wg0

(32c)

subject to


where w ∈ Rn and h(w): Rn f Rm. Problems with general inequality constraints can be transformed to the structure of eq 32 through the introduction of slack variables. The KKT conditions that correspond to eq 32 are

∇wL ) ∇wφ(w) - HTw(w)λ - θ ) 0

(33a)

h(w) ) 0

(33b)

θiwi ) 0

(for i ) 1, ..., n) (w,θ) g 0

(33c) (33d)

where λ and θ are the Lagrange multipliers for the equality and non-negativity constraints, respectively. Hw(w) is the Jacobian of the function h(w), with respect to the vector of variables w. The primal-dual interior point approach involves the application of Newton iterates to the KKT equality constraints with perturbed (relaxed) complementarity conditions,

θiwi ) µ

(for i ) 1, ..., n)

where z ∈ Rnc and v ∈ Rnc are complementarity pairs, w ∈ Rns are variables not associated with complementarity constraints, and h(w,z,v): R2nc+ns f Rnh. As mentioned previously, general inequality constraints can be handled through the introduction of slack variables. The complementarity constraints are a subset of the general equality constraints but are defined separately for the purpose of the discussion that follows. The KKT conditions that correspond to eq 36 result from two sources: the non-negativity constraints (eq 36d) and complementarity constraints present in the original problem (eq 36c). By applying an interior point constraint relaxation to both sets of complementarity constraints, we hope to avoid problems with degeneracy that may occur when one or more of the variables in eq 36c approaches zero too soon. Baker and Swartz37 proposed a strategy of this type in which the same relaxation as that in eq 34 is applied to both sets of complementarity constraints. Raghunathan and Biegler,38 on the other hand, used a complementarity constraint relaxation of the form

ziVi e δµ

(34)

with the parameter µ reduced after each iteration. Strict feasibility of w and θ is maintained, and, under certain assumptions, in the limit as µ f 0, a local solution of the NLP problem is obtained. Key ingredients of the algorithm include a mechanism for reducing the parameter µ and a line search, which, in addition to maintaining the feasibility of w and θ, may be used to impose a sufficient decrease in the objective function and constraint violations. The interior point approach may be viewed as a barrier method in which the sequence of problems

for eq 36c, where µ > 0 is the barrier parameter and δ > 0 is a fixed constant. This was implemented in the software code, IPOPT-C, based on the successful IPOPT interior point code33 for NLPs. Theoretical details are described in the work by Raghunathan and Biegler.39 We use IPOPT-C for the case studies in Section 4. 3.3. Mixed-Integer Programming Formulation and Solution. The inequality and complementarity constraints in the single-level formulation of Section 3.1 may be reformulated as mixed-integer constraints by replacing eqs 20-31 with the following:

n

min φ(w) - µ w∈Rn

ln(wi) ∑ i)1

(35a)

0 e uˆ (k + M - 1|k) - umin e (e - zL(k + M - 1|k))β (38) (35b)

is solved for decreasing values of the barrier parameter µ. In the aforementioned discussion, wi are the scalar components of w. The KKT conditions for eq 35 can be readily shown to correspond to eqs 33a, 33b, and 34. The interior point approach offers a potentially attractive alternative to active set strategies for problems with large numbers of inequality constraints and is implemented in several robust NLP codes.31-33 The references also describe several implementation considerations that have not been taken into account in the previous discussion. Thorough discussions of the interior point method may be found, inter alia, in the work of Roos et al.34 and Wright.35 We consider now the application of an interior point approach to optimization problems with complementarity constraints, into which category our transformed multilevel problem falls. Such problems also occur in game theory, transportation planning, economics, and engineering design.36 We consider the problem

min φ(w,z,v)

(36a)

subject to

h(w,z,v) ) 0 ci(zi,Vi) ) ziVi ) 0, i ) 1,2,...,nc (w,z,v) g 0

(37)

l

subject to

h(w) ) 0

0 e uˆ (k|k) - umin e (e - zL(k|k))β

(36b) (36c) (36d)

0 e umax - uˆ (k|k) e (e - zU(k|k))β

(39)

l 0 e umax - uˆ (k + M - 1|k) e (e - zU(k + M - 1|k))β

(40)

0 e λ4(k|k) e zL(k|k)β

(41)

l 0 e λ4(k + M - 1|k) e zL(k + M - 1|k)β

(42)

0 e λ5(k|k) e zU(k|k)β

(43)

l 0 e λ5(k + M - 1|k) e zU(k + M - 1|k)β

(44)

zL(k|k) ... zL(k + M - 1|k) ∈ {0, 1}nu

(45)

z (k|k) ... z (k + M - 1|k) ∈ {0, 1}

(46)

U

U

nu

where e ) [1, ..., 1]T and β is a suitably large constant. If the remaining constraints are linear, and the objective function is quadratic or linear, we obtain a mixed-integer quadratic program (MIQP) or mixed-integer linear program (MILP) respectively, for which commercial solvers are available (such as CPLEX). We observe further that if the quadratic objective function is convex, the resulting solution in both cases will be a global optimum. Thus, the mixed-integer programming solutions form


Figure 2. Case study 1: mixing tanks in series.

Figure 3. Optimal trajectories for a disturbance change of 1.58%.

a useful benchmark for the complementarity constrained approach; both methods are applied to the case studies presented in the next section. 4. Case Studies 4.1. Case Study 1: Single-Input, Single-Output (SISO) System. 4.1.1. Description. The first case study is a singleinput, single-output (SISO) system that is comprised of three mixing tanks in series, based on an example in the work by Marlin.40 A graphical representation of the system is depicted in Figure 2. The objective is to find an operating point that minimizes the gap between the steady-state composition of component A exiting the final tank, yss () xA3,ss), and its ideal target steady-state value in the absence of disturbances. The manipulated variable (u) is the percentage valve opening for the flow of stream F2, whereas the disturbance (d) is a change in the concentration of component A in stream F1. The transfer function relating the manipulated input to the output is

Gp(s) )

0.039 (5.0s + 1)3

and that which describes the effect of the disturbance is given by

Figure 4. Case study 2: Shell Standard Control Problem.

Table 1. Parameters for Mixing Tank Model Predictive Control (MPC) parameter

description

value

M P Q R ∆t

manipulated variable moves prediction horizon controlled variable weight manipulated variable weight sampling period

2 10 100 0.1 3.0

1.0 Gd(s) ) (5.0s + 1)3 The steady-state relationship between the inputs and the output is

yss ) 0.039uss + dss where the subscript “ss” indicates a steady-state value. The target setpoint is 2% of component A. The following path constraints are imposed:

0% e y e 2.5% 0% e u e 100% The parameters used for the MPC controller are given in Table 1. The reference signal and the controller tuning weights (Q and R) were assumed to remain constant over the prediction horizon. 4.1.2. Results and Discussion. The time horizon for the problem is 33 time steps, and the disturbance is a 1.58% step

Table 2. Solution Time Determined Using the CPLEX and IPOPT-C Programs CPU Time (s) M

complementarities/ integer variables

1 2 3 4 5

68 136 204 272 340

CPLEX

IPOPT-C

1.748 0.9679 3.318

0.5249 0.8209 0.9349 1.3838 1.1638

increase in the composition of component A in stream F1. The case study was solved using both the MIQP and MPCC formulations. Both were coded in AMPL, with the CPLEX solver being used for the former and IPOPT-C for the latter. All variables were initialized to zero.


4.2. Case Study 2: Mutli-Input, Multi-Output (MIMO) System. 4.2.1. Description. The next case study is a subset of the Shell Standard Control Problem41 (SSCP), which is illustrated in Figure 4; this is a benchmark problem that is comprised of a heavy oil fractionator with three product draws and three side circulating loops. The state-space model used was derived from the transfer function model. The control objectives are to maintain the top and side draw product at a steady state of zero (in deviation variables), to minimize the steady-state bottoms reflux duty (u3) and reject unmeasured disturbances in the intermediate reflux (d) that result from changes to the heat duty requirements of other columns. The bottoms reflux temperature (y3) must remain above a specified lower limit. The bottoms reflux duty is directly related to the operating cost; therefore, its minimization corresponds to optimization of an economic objective. We consider here a 2 × 2 controller, in which the top draw end point (y1) and side product end point (y2) are controlled by manipulating the top draw (u1) and the side draw (u2). We address the objectives of bottoms reflux duty minimization and constraint satisfaction of the bottoms reflux temperature in the primary level of the multilevel optimization problem:

Min u3ss subject to

plant and controller equations defining closed-loop dynamics path constraints on the closed-loop response of y1 and y2 path constraints on y3

Figure 5. Optimal steady-state trajectories for a disturbance step of 0.2 from -0.1.

Figure 3 shows that the optimal, closed-loop manipulated variable trajectory saturates at the lower bound for a period of time. Also, to remain feasible, the optimal steady-state setpoint has backed off from the desired target of 2% to 1.60%, so that the upper bound on the concentration of component A is not exceeded. Using IPOPT-C to solve the QP problem with complementarity constraints returned the same objective function value and trajectory as using CPLEX 9 to solve the MIQP problem when CPLEX was able to return a solution. Moreover, because the MIQP formulation is able to guarantee global optimality for this problem, we observe that IPOPT-C is determining the global optimum. We next explore the effect on solution time of increasing problem size. The manipulated variable horizon (M) was increased, resulting in an increase in the number of complementarity pairs and integer variables required. The solution times for both formulations are reported in Table 2. From this observation, we see a general trend: as the manipulated variable horizon M is increased, the solution time increases. However, the solution time of the MIQP grows significantly faster than that of the complementarity problem. For M ) 4 and M ) 5, CPLEX did not return a solution after 12 h. The rapid growth of the solution time of the MIQP is expected, because, in the worst case, mixed-integer solution times grow exponentially. This is in contrast to the modest growth in solution time with the interior point method. Note that no attempt was made to optimize the formulation or tune the algorithm for either the mixed-integer or interior-point approaches.

steady-state model relating inputs and outputs constraints on steady-state inputs and outputs dynamic model relating y3 to inputs Input constraints are handled by the MPC controller. We wish to find an economically optimal operating point, subject to constraints that include bounds on the closed-loop response for a specified disturbance (or disturbance set). The transfer function model of the open-loop multi-input, multi-output (MIMO) process is

[ ]

4.05e-27s 50s + 1 y1(s) -18s 5.39e y2(s) ) 50s + 1 y3(s) 4.38e-20s 33s + 1

[ ]

1.77e-28s 60s + 1 5.72e-14s u1(s) + 60s + 1 u2(s) 4.42e-22s 44s + 1 1.20 -27s e 45s + 1 1.52 -15s e d(s) (47) 25s + 1 1.14 27s + 1

[ ]

[ ]

The steady-state model of the MIMO system is given by

[][

][ ] [ ]

yss1 4.05 1.77 5.88 uss1 1.20 yss2 ) 5.39 5.72 6.90 uss2 + 1.52 dss yss3 4.38 4.42 7.20 uss3 1.14

(48)


Figure 6. Optimal steady-state trajectories for a disturbance step change of -0.2 from 0.1 (scenario 1).

Figure 7. Optimal steady-state trajectories for a disturbance step change of 0.2 from -0.1 (scenario 2).

Table 3. MPC Parameters for the Shell Standard Control Problem

This case study was run using both the MILP and MPCC formulations, with the MPCC formulation being solved using IPOPT-C, and the MILP problem being solved using CPLEX. Both IPOPT-C and CPLEX were initialized with all zeros. The MPCC and MILP formulations contained 9060 complementarity constraints or integer variables, respectively. A solution for the MPCC formulation was determined in 248 s; however, the MILP did not return an integer feasible solution after 12 h. The optimal bottoms reflux duty determined by the MPCC formulation was -0.293. The optimal controlled and manipulated variable trajectories are given in Figure 5. The top draw is located at its upper bound for four sampling periods, beginning at t ) 25 min. The lower bound on y3 is also active at t ) 105 min. Because the MILP did not return an integer feasible solution, it is not possible to verify that the MPCC formulation found the global optimum. No attempt was made to optimize the formulation or tune the algorithm for either the mixed-integer or interior-point approaches. We now consider two distinct disturbance scenarios running in parallel, with the same controller parameters for both scenarios. The first disturbance is a step change in the intermediate reflux of -0.2 from an initial value of 0.1, whereas the second disturbance is a step change to the intermediate reflux of 0.2 from an initial value of -0.1. An optimal bottoms reflux duty is sought, for which feasibility can be maintained for both disturbance scenarios. This case study was run for both the MILP and MPCC formulations, using the CPLEX and IPOPT-C solvers, respec-

parameter

description

value

M P Qy1 Qy2 Ru1 Ru2 ∆t

manipulated variable moves prediction horizon top draw end product weight side draw end product weight top draw weight side draw weight sampling period (min)

15 20 10 10 1 1 5.0

The constraints on the system are

-0.5 e y1 e 0.5

(49a)

-0.5 e y2 e 0.5

(49b)

-0.5 e y3

(49c)

-0.5 e u1 e 0.5

(49d)

-0.5 e u2 e 0.5

(49e)

-0.5 e u3 e 0.5

(49f)

The tuning parameters for the MPC controller can be found in Table 3. The reference signal and controller tuning weights (Q and R) were assumed to remain constant over the prediction horizon. 4.2.2. Results and Discussion. In the first test case, the time horizon for the problem was 150 time steps and the disturbance was a step change to the intermediate reflux of 0.2 from an initial value of -0.1. The discrete time step size was 5 min.


tively. The MPCC and MILP formulations contained 18120 complementarity constraints or integer variables, respectively, and both were initialized with all zeros. A solution for the MPCC formulation was obtained in 205 s; the MILP did not return an integer feasible solution after 12 h. The optimal steady-state bottoms reflux duty determined by the MPCC formulation was -0.205, which is less favorable than the optimal objective function value obtained for the singledisturbance case. The optimal controlled and manipulated variable trajectories for the two disturbance scenarios are given in Figures 6 and 7. We observe that the bottoms reflux temperature constraint is active for the first disturbance scenario, but not the second; the latter corresponds to the disturbance used in the prior single-disturbance formulation. Thus, the inclusion of the more restrictive disturbance scenario results in a decrease in steady-state performance. Inclusion of both of these scenarios within the same optimization problem ensures feasibility for a common bottoms reflux duty. 5. Conclusion A framework for including a constrained model predictive controller within a dynamic optimization problem was presented. The optimization problem is multilevel in nature, because the controller output at each sampling period is determined through the solution of a quadratic programming (QP) problem. Here, we have considered a simultaneous solution approach, in which the convex QP subproblems are replaced by their equivalent Karush-Kuhn-Tucker (KKT) conditions. An interior-point optimization strategy, as implemented in the software code IPOPT-C, was proposed for the solution of the resulting singlelevel complementarity constrained problem. An alternative mixed-integer formulation of the complementarity constraints was presented. With linear models, and a linear or quadratic objective function at the primary level, the complementarity constrained problem may be reformulated as a mixed-integer linear programming (MILP) problem or a mixed-integer quadratic programming (MIQP) problem, respectively. For the problems considered here, the MILP and MIQP solutions are globally optimal, which provides a useful check for the interior-point approach. The interior-point strategy was determined to be highly efficient and to converge to the global optimum for all cases for which this could be confirmed. Although the interior-point approach was observed to significantly outperform the mixed-integer formulation, we recognize that the solution time for mixed-integer problems is often strongly dependent on the manner in which the problem is formulated. A performance comparison of different mixedinteger programming formulations for problems of this type would be interesting. Although the interior-point approach was observed to converge consistently to the global optimum for the class of problems considered, we do not have any guarantee of this yet. A theoretical analysis of this result would be a useful avenue for further study. Although the applications considered in this paper have been confined to linear dynamic models, the method is readily extended to nonlinear systems, discretized using an approach such as orthogonal collocation on finite elements. Application to nonlinear systems in the context of integrated design and control will be the subject of a future communication. Literature Cited (1) Narraway, L. T.; Perkins, J. D.; Barton, G. W. Interaction between process design and process control: Economic analysis of process dynamics. J. Process Control 1991, 1, 243-250.

(2) Narraway, L. T.; Perkins, J. D. Selection of process control structure based on linear dynamic economics. Ind. Eng. Chem. Res. 1993, 32, 26812692. (3) Narraway, L.; Perkins, J. Selection of process control structure based on economics. Comput. Chem. Eng. 1994, 18 (Suppl.), S511-S515. (4) Figueroa, J. L.; Bahri, P.; Bandoni, J. A.; Romagnoli, J. A. Economic impact of disturbances and uncertain parameters in chemical processessA dynamic back-off analysis. Comput. Chem. Eng. 1996, 20 (4), 453-461. (5) Bahri, P. A.; Bandoni, J. A.; Barton, G. W.; Romagnoli, J. A. Backoff calculations in optimising control: A dynamic approach. Comput. Chem. Eng. 1995, 19 (Suppl.), S699-S708. (6) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal design of dynamic systems under uncertainty. AIChE J. 1996, 42 (8), 2251-2272. (7) van Schijndel, J. M. G.; Pistikopoulos, E. N. Towards the integration of process design, process control & process operabilitysCurrent status & future trends. In Foundations of Computer-Aided Process Design; Malone, M. F., Trainham, J. A.; Carnahan, B., Eds.; American Institute of Chemical Engineers: New York, 2000; pp 99-112. (8) Schweiger, C. A.; Floudas, C. A. Interactions of design and control: Optimization with dynamic models. In Optimal Control: Theory, Algorithms, and Applications; Hager, W. W., Pardalos, P., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1998; pp 388-435. (9) Pistikopoulos, E. N.; Sakizlis, V. Simultaneous design and control optimization under uncertainty in reaction/separation systems. In Chemical Process Control-VI: Assessment and New Directions for Research; Rawlings, J. B.; Ogunnaike, B. A., Eaton, J. W., Eds.; AIChE Symposium Series No. 326, Volume 98; CACHE and American Institute of Chemical Engineers: New York, 2002; pp 223-238. (10) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Effect of disturbances in optimizing control: Steady-state open-loop backoff problem. AIChE J. 1996, 42 (4), 983-994. (11) Young, J. C. C.; Baker, R.; Swartz, C. L. E. Input saturation effects in optimizing controlsInclusion within a simultaneous optimization framework. Comput. Chem. Eng. 2004, 28, 1347-1360. (12) Baker, R.; Swartz, C. L. E. Rigorous handling of input saturation in the design of dynamically operable plants. Ind. Eng. Chem. Res. 2004, 43, 5880-5887. (13) Baker, R.; Swartz, C. L. E. Simultaneous solution strategies for inclusion of input saturation in the optimal design of dynamically operable plants. Optim. Eng. 2004, 5, 5-24. (14) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control Eng. Pract. 2003, 11, 733-764. (15) Loeblein, C.; Perkins, J. D. Structural design for on-line process optimization: I. Dynamic economics of MPC. AIChE J. 1999, 45 (5), 10181029. (16) Brengel, D. D.; Seider, W. D. Coordinated design and control optimization of nonlinear processes. Comput. Chem. Eng. 1992, 16 (9), 861-886. (17) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. N. Parametric controllers in simultaneous process and control design optimization. Ind. Eng. Chem. Res. 2003, 42, 4545-4563. (18) Pistikopoulos, E. N.; Dua, V.; Bozinis, N. A.; Bemporad, A.; Morari, M. On-line optimization via off-line parametric optimization tools. Comput. Chem. Eng. 2002, 26, 175-185. (19) Cutler, C. R.; Ramaker, B. L. Dynamic Matrix ControlsA computer control algorithm. Presented at the AIChE National Meeting, Houston, TX, 1979. (20) Maciejowski, J. M. PredictiVe Control with Constraints; Prentice Hall: Harlow, U.K., 2002. (21) Zafiriou, E. Robust model predictive control of processes with hard constraints. Comput. Chem. Eng. 1990, 14, 359-371. (22) Clark, P. A.; Westerberg, A. W. Bilevel programming for steadystate chemical process designsI. Fundamentals and algorithms. Comput. Chem. Eng. 1990, 14 (1), 87-97. (23) Fortuny-Amat, J.; McCarl, B. A representation and economic interpretation of a two-level programming problem. J. Oper. Res. Soc. 1981, 32, 783-792. (24) Bard, J. F.; Moore, J. T., A branch and bound algorithm for the bilevel programming problem. SIAM J. Sci. Stat. Comput. 1990, 11, 281292. (25) Gu¨mu¨s¸ , Z. H.; Floudas, C. A. Global optimization of nonlinear bilevel programming problems. J. Global Optim. 2001, 20, 1-31. (26) Bard, J. F. Practical BileVel Optimization: Algorithms and Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1998. (27) Fletcher, R.; Leyffer, S. Solving mathematical programs with complementarity constraints as nonlinear programs. Optim. Methods Software 2004, 19 (1), 15-40. (28) Leyffer, S. Complementarity constraints as nonlinear equations: Theory and numerical experience. In Optimization with MultiValued

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 91 Mappings: Theory, Applications, and Algorithms; Dempe, S., Kalashnikov, V., Eds.; Springer Optimization and Its Applications, Vol. 2; Springer: New York, 2006; pp 169-208. (29) Anitescu, M. On using the elastic mode in nonlinear programming approaches to mathematical programs with complementarity constraints. SIAM J. Optim. 2005, 15 (4), 1203-1236. (30) Soliman, M. A Strategy for Inclusion of Closed Loop Dynamics in Real Time Optimization with Application to an Oxygen Delignification Reactor, Master’s Thesis, Department of Chemical Engineering, McMaster University, Hamilton, Ontario, Canada, 2004. (31) Vanderbei, R. J.; Shanno, D. F. An interior-point algorithm for nonconvex nonlinear programming. Comput. Optim. Appl. 1999, 13, 231252. (32) Byrd, R. H.; Hribar, M. E.; Nocedal, J. An interior point algorithm for large-scale nonlinear programming. SIAM J. Optim. 1999, 9 (4), 877900. (33) Wa¨chter, A.; Biegler, L. T. On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming. Math. Progr. 2006, 106 (1), 25-57. (34) Roos, C.; Terlaky, T.; Vial, J.-P. Theory and Algorithms for Linear Optimization: An Interior Point Approach; Springer Science: New York, 2005. (35) Wright, S. J. Primal-Dual Interior Point Methods; SIAM: Philadelphia, PA, 1997.

(36) Grossman, I. E.; Biegler, L. T. Part II: Future perspective on optimization. Comput. Chem. Eng. 2004, 28, 1193-1218. (37) Baker, R.; Swartz, C. L. E., Simultaneous solution strategies for inclusion of input saturation in the optimal design of dynamically operable plants. Optim. Eng. 2004, 5 (1), 5-24. (38) Raghunathan, A. U.; Biegler, L. T. Mathematical programs with equlibrium constraints (MPECs) in process engineering. Comput. Chem. Eng. 2003, 27, 1381-1392. (39) Raghunathan, A. U.; Biegler, L. T. An interior point method for mathematical programs with complementarity constraints (MPCCs). SIAM J. Optim. 2005, 15 (3), 720-750. (40) Marlin, T. E. Process Control: Designing Processes and Control Systems for Dynamic Performance, 2nd Ed.; McGraw-Hill: New York, 2000. (41) Prett, D. M.; Garcia, C. E.; Ramaker, B. L. The Second Shell Process Control Workshop: Solutions to the Shell Standard Control Problem; Butterworths: Boston, 2000.

ReceiVed for reView February 20, 2007 ReVised manuscript receiVed September 14, 2007 Accepted September 20, 2007 IE070270R

Interior Point Solution of Multilevel Quadratic Programming Problems

Recommend Documents