Ind. Eng. Chem. Res. 2009, 48, 5767–5778
5767
A Computationally Efficient Scheme for Model Predictive Control of Nonlinear Hybrid Systems Using Generalized Outer Approximation† Naresh N. Nandola and Sharad Bhartiya* Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai 400 076, India
This paper presents an efficient optimization algorithm suitable for online solution of mixed integer nonlinear programs resulting from the model predictive control (MPC) of nonlinear hybrid systems. The system model is based on a recently proposed multiple partially linear (MPL) modeling scheme. The algorithm, based on generalized outer approximation (GOA), uses structural information of the canonical MPL framework as well as analytical expressions for the objective function and constraints of a relatively simple primal problem as well as the master problem. Specifically, the primal problem of GOA reduces to a quadratic program when MPL models are used in MPC. Computational efficiency of the algorithm over the branch and bound strategy is demonstrated using a simulated benchmark three-spherical tank system and a hydraulic process plant. 1. Introduction Hybrid systems refer to processes that involve continuous dynamics in addition to discrete (logical) decisions1,2 and have found applications in manufacturing systems, automobile control, and computer disk drive control among others. Although a framework for modeling and control of chemical processes as hybrid systems has emerged only recently, large continuous plants have always used logic controllers to implement safety features such as the triggering of a coolant pump and the various safety interlocks. However, current trends in the chemical process industry emphasize the need for flexible processing, which invariably necessitates a greater degree of logical decision-making along with the continuous control laws even during the course of normal operation. A number of modeling formalisms that represent nonlinear hybrid dynamical systems (NHDS) have been proposed in literature.1,3,4 These formalisms can be broadly classified into the following three categories:5 (i) a discrete formalism, such as finite automata, that can be extended with continuous variables resulting in a hybrid framework such as timed automata and hybrid Petri Nets, (ii) a continuous formalism that can accommodate discrete variables or logical conditions by appropriately switching between system dynamics, and (iii) an approach that directly combines the continuous subsystem with its discrete counterpart through an interface. These models play a key role in various applications of hybrid systems such as simulation, verification, identification, optimization, and control. Among these, optimal feedback control of NHDS, such as model predictive control, is particularly challenging since it requires an online solution of a mixed integer nonlinear/quadratic program (MINLP/MIQP) within a small fraction of the sampling period. This impediment to online control of NHDS can be addressed along three paths: (i) efficient representation of the NHDS, (ii) efficient algorithms for solution of MINLP/MIQP, and (iii) enhanced computer speed. Our earlier work addressed the problem along the first path by modeling the NHDS using a multiple, partially linear (MPL) scheme,6 which is related to the well-known mixed logical dynamical (MLD) model.2 In the MPL scheme, each linear † An abbreviated version of this paper has been presented in the 17th IFAC World Congress 2008, Seoul, South Korea. * To whom correspondence should be addressed. E-mail: bhartiya@ che.iitb.ac.in. Tel.: +91 22 2576 7225. Fax: +91 22 2572 6895.
model is a local representation of all locations of the hybrid system. These models are then combined using Bayes’ theorem to describe the nonlinear hybrid system. The multiple models, which consist of continuous as well as discrete variables, are used for synthesis of a model predictive control (MPC) law. In the absence of discrete states or inputs, the MPC formulation reduces to that used for discrete-time control of a continuous variable system. Although implementation of the control law for hybrid systems using the MPL model requires an online solution of an MINLP, the optimization problem has a canonical structure with certain computational advantages. These advantages accrue because of the use of the MPL model rather than the MLD model and were demonstrated for control of a benchmark three tank system and a PLC driven hydraulic process plant, where all integer programs were solved using a branch and bound (BB) strategy.6 The current work overcomes the second impediment to control of NHDS by exploiting the canonical structure of the MINLP resulting from the MPL model based predictive control law. We show that this MINLP is particularly suited for efficient solution by reduced space decomposition approaches such as generalized outer approximation (GOA) and generalized bender decomposition (GBD). Specifically, we show that in solving the MINLP resulting from the MPL based control law, the primal problem at each iteration of the GOA algorithm reduces to a quadratic program (QP) while the master problem retains its mixed integer linear program (MILP) form. Thus the GOA method obtains the solution of the original MINLP by solving a series of QPs and MILPs. Both of these subproblems are comparatively less expensive than the NLP solutions needed when using the BB strategy. Moreover, the canonical structure of the MINLP enables us to derive analytical expressions for the gradients of the nonlinear objective function and constraints, which can be used in the master problem to further speed up the solution. The benefits of this algorithm relative to the BB strategy are demonstrated through simulation examples. The results confirm significant computational advantage relative to our previous results in ref 6 that used a BB strategy. The current paper is organized as follows: section 2 reviews modeling and control of nonlinear hybrid dynamic systems based on MPL models; section 3 discusses the GOA algorithm in perspective of our MPL model based predictive control; section 4 demonstrates the superior performance of GOA algorithm over the
10.1021/ie801384y CCC: $40.75 2009 American Chemical Society Published on Web 05/18/2009
5768
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
BB strategy for MPL model based control of NHDS on two applications, namely, a three-spherical-tank system and a hydraulic process plant. Finally, the work is concluded in section 5. 2. Modeling and Control of Hybrid System Using MPL Models: An Overview Hybrid systems may involve both continuous and discrete states as well as continuous and discrete inputs. Here, the flowfield describing the evolution of continuous states is dependent on discrete phenomena characterized by discrete events such as state events and control events. A formal description of hybrid systems is represented by a hybrid automaton.4,7 The choice of the specific model structure of the hybrid automaton, however, is dependent on the application at hand. For example, a simpler model, usually an overapproximation of the original model, is employed in formal verification. Here, timed automata8 or rectangular automata5 are common choices. On the other hand, control as well as fault detection and diagnosis of hybrid systems have led to development of various linear modeling paradigms such as piecewise affine9 (PWA) and mixed logical dynamic2 (MLD) models. The equivalence between MLD and a variety of other linear hybrid models has also been proven.10 In the area of nonlinear hybrid systems, Buss et al.3 model the nonlinear hybrid system by introducing discrete states (xd) as well as discrete control inputs in addition to the continuous states (xc) and control inputs (uc). Their model, known as the hybrid state model (HSM), forms the starting point for the development of multiple partially linear (MPL) model,6 which we briefly review in the remainder of this section. Here, the continuous states of the hybrid system evolve based on the flow-field fl, which is dependent on the location, l, of the system. Upon occurrence of an event, the system transits to a new location l′ and evolves based on a new flow-field fl′. To identify the different locations of the systems and the transitions between them, suitable event-generating functions are defined. When one or more of these functions take on a value of zero, a discrete event is said to occur, which includes state as well as control events. To obtain a control relevant model that represents all locations of the NHDS by a global flow-field, the HSM may be modified by assigning binary indicator variables δj ∈ {0,1}, j ) 1, ..., ns to each event generating function. The resulting logical expressions are then converted into linear inequalities using equivalence with propositional logic expressions (such as Big-M constraints; see, for example, Williams11 and Raman and Grossmann12). Thus, the modified HSM is written as follows: x˙c ) fg(xc, uc, δ)
(1)
xd(t) ) bd(δ(t))
(2)
E1uc(t) + E2δ(t) + E3x(t) e E4
(3)
xc(t+) ) xc(t-)
(4)
xd(t+) ) bd(δ(t+))
(5)
where x ) [xc xd]T, fg is a global flow-field that subsumes all location-dependent flow-fields fl and is parametrized by the indicator vector δ, whose elements are determined by inequality 3 that captures the events of the hybrid system. Note that the discrete state xd is related to the continuous variable through the indicator variable δ. Superscripts + and - indicate values
of states just after and before the occurrence of an event, respectively. Equation 4 captures the fact that we are modeling switched systems. Matrices Ei (i ) 1, 2, 3, 4) are constant coefficient matrices and vector bd is a function of binary variables δ. Note that the vector δ, consisting of ns binary elements, can describe 2ns locations. A change in the status of one or more elements of δ corresponds to an event that may be triggered by discontinuity in states (that is, a state event (SE)) and/or discontinuity in inputs (that is, a control event (CE)). Remark 1: Elements of δ describe discrete control inputs in addition to discrete states (see Equation 2) of the NHDS. Thus, occurrence of control eVents and state eVents can be identified from the status of binary indicator Variables δj ∈ {0,1},j ) 1, ..., ns. Model predictive control of NHDS described by model 1-5 requires an online solution of an MINLP, which is generally difficult. To overcome this problem, the model was simplified by linearization of the flow-field at different locations and operating conditions followed by their reaggregation, which results in the multiple partially linearized (MPL) model scheme.6 2.1. Elements of the MPL Model. The MPL model development begins with performing a Taylor series expansion of 1 at an operating point characterized by the continuous variables (xc, uc) and retaining the binary variables δ as a parameter. Thus, one obtains a linear model whose system matrices are a function of the binary variables δ. The continuous-time linear model is discretized, for all possible instantiations of the binary variables δ, thereby obtaining discrete-time models corresponding to all locations of the hybrid systems. Since ns binary variables result in 2ns possible locations, one can obtain 2ns discrete time linear models from the continuous time model represented by 1 and 2. These 2ns models are then combined using a corresponding scalar logical multiplier li, i )1, ..., 2ns to produce a composite, discrete-time representation of 1-3 as follows,6,13 2ns
xk+1 ) (
∑ l (δ )Φ )x i
i)1
k
i
k
2ns
+(
∑
2ns
li(δk)Γl)uck
+(
i)1
∑ l (δ )f i
k di)
i)1
(6) E1uck + E2δk + E3xk e E4
(7)
The logical multiplier li is defined using the binary variables δk, and is designed to take on a value 1 if and only if the ith mode of the hybrid system represented by the ith combination of the binary variables is encountered and zero otherwise. Remark 2: The nondeViation form of Variables has been used as this allows representation of nonequilibrium operation, a common feature in hybrid systems applications, resulting in the affine representation. Note that the RHS of eq 6 consists of multiplicative terms between binary variable δk, state variables xck, and inputs uck, and thus the model is nonlinear. These multiplicative terms can be masked by introducing auxiliary binary and auxiliary continuous variables and their corresponding constraints, thereby reformulating eqs 6 and 7 into the well-known MLD model.2 However, the increased size of the MLD model imposes a large computational burden in its use for optimal control. On the other hand, retaining it in the nonlinear form is computationally efficient as it requires fewer number of variables and constraints.6 The model in eqs 6 and 7 can be expressed in a compact form using vector notation as follows: ¯ )xk + (L¯k(δk)Γ¯ )uck + L¯k(δk)f¯ xk+1 ) (L¯k(δk)Φ
(8)
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
E1uck
+ E2δk + E3xk e E4
(9)
j k, Φ j, Γ j , and jf are constituted where the augmented matrices L j k may be interpreted as from li,k, Φi, Γi, and fdi, respectively. L an operator that retrieves the system matrices corresponding to the current location as identified by δk from the composite j,Γ j , and jf. The outputs of the linearized model system matrices Φ may be written as follows: yk ) Cxk
(10)
Note that this model accounts for all locations of the nonlinear hybrid system in the vicinity of a single operating point. Similar linear discrete-time models may be obtained at different operating points characterized by the continuous variables (xc, uc).14 These models are then combined using a multiple model approach such as Bayesian weighting6,15 to reconstitute the original nonlinear model. Thus, the multiple-model approach may be represented in a canonical form as follows, ¯ avg)xk + (L¯k(δk)Γ¯ avg)uck + L¯k(δk)f¯avg xk+1 ) (L¯k(δk)Φ (11) E1uck + E2δk + E3xk e E4
(12)
yk ) Cxk
(13)
j avg, jfavg are the blended system matrices that depend j avg, Γ where Φ on weighting of the different models. Equations 11-13 are referred to as the multiple partially linearized model, and it approximates the nonlinear operating range as well as all locations of the hybrid system. The canonical structure of the MPL model can potentially describe any arbitrary NHDS. In addition, this model becomes linear at a fixed location. Nandola and Bhartiya6 formulated the MPC problem for the NHDS described by MPL models as follows:
(
¯ 1k(δ¯ k)xk + H ¯ 2k(δ¯ k)µck + H ¯ 3k(δ¯ k) - ψref)T (H c ¯ 1k(δ¯ k)xk + H ¯ 2k(δ¯ k)µk + H ¯ 3k(δ¯ k) - ψref) + min J ) Wy(H µck,δk
c c (Rµck - R0uk-1 )TWu(Rµck - R0uk-1 )
)
(14)
such that, g1: E¯1µck + E¯2δ¯ k + E¯3(H1k-1(δ¯ k-1)xk + H2k-1(δ¯ k-1)µck + H3k-1(δ¯ k-1)) e E¯4 (15) ¯ 1k(δ¯ k)xk + H ¯ 2k(δ¯ k)µck + H ¯ 3k(δ¯ k) e ψmax g2: ψmin e H (16) µcmin e µck ∈ Rnc e µcmax,
δ¯ k ∈ {0, 1}nb
(17)
j 2k(δj k), H j 3k(δj k), H1k-1(δj k-1), H2k-1(δj k-1), j 1k(δj k), H where H j k-1) are variable coefficient matrices which depend on H3k-1(δ j k. Expressions for these matrices the discrete decision variables δ are presented in ref 6 and are summarized in appendix I for j k consist of continuous inputs and completion. Vectors µck and δ binary variables based on a p step ahead prediction for m control moves using the model 11-13 and constitute the decision variables of the MINLP. Vector ψref represents the set point j 3, and E j 4 are constant coefficient j 2, E j 1, E trajectory. Matrices E matrices made up of matrices E1, E2, E3, and E4, respectively. Note that the nonlinear dependence of the objective function 14 and constraints, g1 and g2, on the binary variables, makes the above optimization problem an MINLP. While we have not
5769
considered stabilizing constraints in the above, the MPC formulation can be modified to include the terminal cost and constraint set as discussed in ref 16. A popular approach for solving the MINLP, such as in eqs 14-17, uses branch and bound (BB) strategy. The BB strategy performs a tree search in the space of integer variables and solves a relaxed NLP at each node where a subset of binary variables is fixed. The relaxed NLP provides a lower bound to the original MINLP and this information is used to fathom the BB nodes. This method is computationally expensive for high dimensional binary variable search space and can be of use only for few binary variables or when the relaxed NLP is easy to solve.17 On the other hand, decomposition based algorithms such as generalized bender decomposition (GBD)18,19 and generalized outer approximation (GOA)20 have proved computationally efficient for a certain class of MINLP such as mixed-integer dynamic optimization (MIDO)21 and generalized disjunctive programming.22 In the next section, we discuss the GOA method in light of the MINLP resulting from formulation of the MPLmodel-based control law. The chief advantage of the MPL model is that although the control law continues to be an MINLP, it takes on a canonical form for which the GOA algorithm is particularly suitable. 3. Generalized Outer Approximation (GOA) for MPL-Model-Based Control The key idea of the GOA algorithm20 is generation of a nonincreasing upper bound and a non-decreasing lower bound by solving a series of primal and master problems. The primal problem, obtained by fixing the binary variables, is an NLP whose solution represents an upper bound to the MINLP. The master problem, based on support functions of the nonlinear objective functions as well as an outer approximation of the nonlinear constraints that are active at the solution of the primal, results in an MILP whose solution is a lower bound to the original MINLP if the MINLP is convex. The solution of the master problem also produces the value of the binary variables for the next iteration. When the difference between the lower and upper bounds lies within a user-defined tolerance, the algorithm terminates and the current solution of the primal problem and corresponding binary variables are considered optimal solution for the MINLP.19 Viswanathan and Grossmann23 consider a rise in the value of the upper bound as a termination criterion for nonconvex MINLP. In such a case the optimal solution corresponds to the solution of primal problem at previous iteration. We have used both the above criteria for termination in the simulation examples reported below. For nonconvex problems, however, the linear approximations needed in the master problem may not result in an outer approximation of the original MINLP and hence may sometimes result in a suboptimal solution. In such problems, convexification methods for MINLP24 or augmented penalty based methods23 may be used prior to application of this algorithm. In the current work, however, remedial measures for nonconvexity have not been pursued. 3.1. Primal Problem. The objective function 14 of this j k and quadratic in MINLP is nonlinear in binary variables δ c continuous variables µk. This feature arises because the MPL model takes on an affine form at a fixed location of hybrid system. The constraints represented by 15 and 16 are nonlinear j k but linear with respect to continuous in binary variables δ c variables µk. Therefore, upon fixing the binary variables in 14-16, the primal problem reduces to the following QP:
5770
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
min Jp,i ) µck
( 21 (µ ) Q (µ ) + F µ ) c T k
i
c k
T c i k
(18)
such that, Aiµck e bi
(19)
µcmin e µck e µcmax
(20)
respectively. While these gradients may be calculated numerically, the canonical structure of the MINLP enables one to write analytical gradients, which can be used for enhanced accuracy as well as efficient computation. In the simulation examples presented later, we use the analytical gradients which are summarized in the remainder of this section. 3.2.1. Analytical Gradient of Objective Function. The gradient of the objective function 14 can be written as
[]
∂J ∂µck ∇J ) ∂J ∂δ¯
where Qi, Ai, Fi, and bi are constant matrices and can be easily j k,i in eq 14-16 derived by fixing the value of binary variables δ as follows: ¯ T2k,iWyH ¯ 2k,i + RTWuR) Qi ) 2(H
[
j1 + E j 3H2k-1,i E ¯ Ai ) H2k,i ¯ 2k,i -H FTi
(
(21)
]
(22)
[
j4 - E j 2δ¯ k,i - E j 3H1k-1,ixk - H3k-1,i) (E ¯ ¯ 3k,i - dk) bi ) (ψmax - H1k,ixk - H ¯ 1k,ixk + H ¯ 3k,i + dk) (-ψmin + H
]
)
(23)
(24)
In the above QP, suffix i denotes the GOA iteration and suffix k indicates the time instant. Note that the solution of this problem c and the yields the current value of the continuous variables µk,i upper bound JUB,i ) Jp,i of the MINLP. 3.2. Master Problem. For the case of a feasible primal problem, the master problem represents the outer approximation of the nonlinear objective function 14 as well as the active nonlinear constraints 15 and 16 obtained by linearization at (µck,i, j k,i) and the constraints on the bound eq 17. On the other hand, δ if the primal problem is infeasible, then the master problem is obtained by linearization of the active nonlinear constraints alone.19 This is achieved by solving a feasibility problem which minimizes some norm of the constraint violation. Thus, the master problem reduces to the following mixed integer linear program (MILP) min Jm,i ) RGOA
[ ] [ ] [ ] }
e RGOA ¯δ - δ¯ k k,i µck - µck,i ga(µck,i, δ¯ k,i) + ∇gTa e0 δ¯ - δ¯
J(µck,i, δ¯ k,i) + ∇J
c T µk
µck,i
k
ga(µck,i, δ¯ k,i) + ∇gTa
k,i
∂J ¯ 1k(δ¯ k)xk + H ¯ 2k(δ¯ k)µck + H ¯ 3k(δ¯ k) ) 2(H ∂µck c ¯ 2k(δ¯ k) + 2(Rµck - R0uk-1 )TWuR ψref)TWyH ∂J ¯ 1k(δ¯ k)xk + H ¯ 2k(δ¯ k)µck + H ¯ 3k(δ¯ k) ) 2(H ¯ ∂δk ¯ 1k(δ¯ k) ¯ 2k(δ¯ k) ¯ 3k(δ¯ k) ∂H ∂H ∂H xk + µck + ψref)TWy ∂δ¯ k ∂δ¯ k ∂δ¯ k
(
(30)
)
(31)
j k))/(∂δ j k), (∂H j 2k(δ j k))/ j 1k(δ Analytical expressions for derivatives (∂H j j j j (∂δk), (∂H3k(δk))/(∂δk) are provided in appendix II. 3.2.2. Gradient of Nonlinear Constraints g1 and g2. The gradient of the nonlinear constraints g1 in eq 15 can be obtained j k, we as, taking partial derivative of 15 with respect to µck and δ obtain
[] ∂g1
∇g1 )
[
∂µkc ∂g1 ∂δ¯
k
j1 + E j 3H2k-1(δ¯ k-1) E ∂H (δ¯ ) ∂H (δ¯ ) ∂H (δ¯ ) ) j2 + E j 3 1k-1 k-1 xk + E j 3 2k-1 k-1 µkc + E j 3 3k-1 k-1 E ¯ ¯ ¯ ∂δk ∂δk ∂δk
}
]
(32)
(25)
Similarly, gradient for constraints g2 in eq 16 can be obtained as follows: ∀i ∈ FP
µck - µck,i e 0 ∀i ∈ NFP δ¯ k - δ¯ k,i
µcmin e µck ∈ Rnc e µcmax,
k
where
¯ T1k,iWyH ¯ 2k,i + H ¯ T3k,iWyH ¯ 2k,i + xTk H )2 T T c ¯ ¯ )TRT0 WuR dk,iWyH2k,i - ψrefWyH2k,i - (uk-1
µck,δk,RGOA
(29)
δ¯ k ∈ {0, 1}nb
(26) ∇g2 ) (27)
[] ∂g2
[[ [
(28)
where variable RGOA represents an upper bound of the support of the objective function. FP and NFP stand for feasible and nonfeasible primal problem, respectively.As previously, suffix i denotes the GOA iteration and suffix k indicates the time instant. Vector ga represents all nonlinear constraints that are active at the solution of the primal while ∇J, ∇ga are gradients of objective function 14 and active nonlinear constraints, ga,
∂µck ∂g2 ∂δ¯
)
k
¯ 2k(δ¯ k) H ¯ 2k(δ¯ k) -H ¯ ∂H1k(δ¯ k) ∂δ¯
]
k
¯ 1k(δ¯ k) ∂H ∂δ¯ k
] [ ] [ ]] xk +
¯ 2k(δ¯ k) ∂H ∂δ¯ k
¯ 2k(δ¯ k) ∂H ∂δ¯ k
µck +
¯ 3k(δ¯ k) ∂H ∂δ¯ k
¯ 3k(δ¯ k) ∂H ∂δ¯ k
(33)
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
πh2(hmax - h2)
5771
dh2 ) (Q2 - V23Q23V23 - V2Q23V2) dt (36)
πh3(hmax - h3)
Figure 1. Schematic of the 3-spherical tank problem.
Combining eqs 32 and 33, we get gradient of nonlinear constraints,∇g as follows, ∇g )
[∇gT1
∇gT2 ]T
where V1, V2, V13, V23, VL1, and VN3 represents binary indicator variables for corresponding valves, and hmax represents height of the spherical tank (0.6 m). Variables Qi represent flow rates through valves Vi and may be evaluated using the following constitutive equations: Qi3Vi3 ) azSi3 sign(hi - h3)√ |2g(hi - h3)|,
i ) 1, 2 (38)
(34)
j k-1))/∂δ j k, Expressions for evaluation of matrices (∂H1k-1(δ j j j j (∂H2k-1(δk-1))/∂δk, and (∂H3k-1(δk-1))/∂δk are documented in appendix II. In summary, the primal problem is a QP given by eqs 18-20, while the master problem is an MILP, presented in eqs 25-28. The linearization of the objective functions and active nonlinear constraints required in the master are represented by eqs 29-34. 4. Application We present two simulation examples to demonstrate the computational advantages of the GOA algorithm for online solution of the MINLP that arises during the MPL-based model predictive control. Specifically, we compare the average computation time needed to solve the identical control problem with the BB strategy. A solver quadprog of MATLAB 6.5 (Mathworks Inc., Natick, MA) based on the active set strategy was used to solve the primal QP while the MILP resulting from the master problem was solved using Tomlab/Cplex. The BB implementation was based on a Matlab based solver BNB20 (Matlab Central). BNB20 uses a BB shell integrated with Matlab’s NLP solver fmincon, which is an implementation of the npsol algorithm. We found that BNB20 with fmincon was computationally more efficient than the Tomlab/minlpbb solver. However, Tomlab/Cplex produced better computational performance of MILP than BNB20 with linprog of Matlab which is probably due to different pivoting strategies and was therefore used to solve the master problem. All NLP and QP algorithms were based on the active set method. All simulations have been performed on a 3.0 GHz P-IV machine with 1 GB RAM. For comparison between the GOA and BB strategies, the computation time was obtained using the subroutine cputime in matlab. 4.1. Spherical Three-tank System. As shown in Figure 1, the system consists of two independent pumps that deliver the liquid flow rates Q1 and Q2 to tank 1 and tank 2, respectively through the two control valves. Six independent solenoid (on/ off) valves (V1, V2, V13, V23, VL1, and VN3) can be manipulated to interrupt the flows into or out of the three tanks. Tank 1 and tank 3 as well as tank 2 and tank 3 are connected through upper and lower pipes. To enhance the nonlinear behavior, we replaced the cylindrical tanks as reported in the benchmark problem25 with spherical tanks.6 The first principles model of the setup is briefly described below: πh1(hmax - h1)
dh3 ) (V13Q13V13 + V23Q23V23 + V1Q13V1 + dt V2Q23V2 - VN3QN) (37)
dh1 ) (Q1 - V13Q13V13 - V1Q13V1 - VL1QL1) dt (35)
QL1 ) azSL1√2gh1
(39)
QN ) azSN√2gh3
(40)
Qi3Vi ) (azSi sign(max{hi, hV} max{h3, hV}) × √2g|(max{hi, hV} - max{h3, hV})|),
i ) 1, 2 (41)
where Si, Si3, SL1, and SN are cross sectional areas of valves and assumed identical for all valves (0.95 cm2), hV (see Figure 1) is height of upper pipe from bottom (0.3 m), and az is the discharge coefficient, which is assumed to be unity. Next, we consider three partially linearized models for the three-spherical-tank system. The points of linearization are listed below: h1 ) h2 ) 0.15 m,
(i) Model I: h3 ) 0.14 m, Q1 ) Q2 ) 0
(ii) Model II: h3 ) 0.24 m,
Q 1 ) Q2 ) 0
(iii) Model III: h1 ) h2 ) 0.35 m, h3 ) 0.34 m,
Q 1 ) Q2 ) 0
h1 ) h2 ) 0.25 m,
Model I and model II correspond to levels below the upper pipe connections while model III corresponds to a level above the upper pipe. The three points correspond to low, medium, and high levels in the three tanks. The points of linearization were chosen such that the continuity of the max function is maintained. Alternatively, smooth approximation of the max function may be used. We use the MPL model 11-13 based MPC for a set-point tracking control problem, which involves the filling of empty tanks to desired levels, followed by multiple set-point changes. To study the computational efficiency of the GOA algorithm, we compare the computation time with the results obtained using BB. Both cases use a sampling time ts ) 3 s, prediction horizon, p ) 5 (that is, 15 s) and control horizon, m) 2 (that is, 6 s). We consider the first principles nonlinear model 35-41 as the plant model and the MPL framework 11-13 as the controller model. Figure 2 documents the results of the MPL-model-based control of the levels in the three tanks at different target levels using both the BB algorithm (dotted line) as well as the GOA algorithm (solid line). Note that the state trajectories for both cases are found to be similar and hence almost indistinguishable in the figure. Set point changes were deliberately chosen so that
5772
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
Figure 2. MPL-model-based predictive control of levels h1, h2, and h3 in the 3-spherical tank system using GOA (solid line) and using BB (dotted line). The upper pipes connecting the three tanks are located at 0.3 m. The set-points (dashed-dotted line) for the three levels [h1, h2, h3] are as follows: at 0 s [0.45 0.3 0.2], at 153 s [0.3 0.4 0.3], at 363 s [0.45 0.25 0.2], at 903 s [0.25 0.45 0.3].
Figure 3. Control moves for the level control problem in the 3-spherical tank system using GOA (solid line) and using BB (dotted line).
the discrete state events, characterized by the crossing of the upper pipes at 0.3 m by the water levels, occur. Figure 3 shows the corresponding control moves for both cases. The noisy behavior of the flow rate to tank 1 between 400 and 750 s, Q1, corresponds to the chattering in the solenoid valve VL1. During this period, the valves V1 and V13 connecting tank 1 to tank 3 remain closed and the discharge of water through the solenoid valve VL1 is matched by the filling through the continuous valve Q1. Both the algorithms produce similar results when used for the solution of the MINLP resulting from MPL-model-based control. However, the main advantage of GOA algorithm lies in its computational efficiency. The BB algorithm solves a relaxed NLP at each node while at each iteration of the GOA algorithm, a QP and an MILP are solved. Thus, the GOA algorithm solves comparatively easier optimization problems than BB but it requires additional computation time for linearization of the nonlinear objective function as well as constraints to obtain the master problem. However, availability of analytical gradients reduces the computation effort. Figure 4 documents the computation time required to solve the control problem at each sampling time using BB and GOA algorithms. The average computation time and standard deviation for
Figure 4. Comparison of computation time needed for online solution of the MINLP using GOA (+) and BB (o) strategies. The horizontal dash-dotted line indicates the sampling period of 3 s. The vertical dotted lines indicate the samples where set-point changes were implemented.
computation of the control move determined through the BB strategy are 1.99 and 2.8 s, respectively. Further, about 16% of the optimization problems using the BB algorithm required computation time greater than one sampling period. This clearly shows that the BB algorithm is impractical for online implementation. On the other hand, the average computation time and standard deviation for the GOA algorithm are noted as 0.12 and 0.083 s, respectively. Moreover, the GOA algorithm takes less than one sampling time for all instants. Ideally, values of the manipulated variables should be injected immediately upon availability of the measurement. This requires that the MINLP must be solved within a very small fraction of the sampling time. Since the GOA algorithm takes only about 4% of the sampling period on an average, it is likely to show better practical behavior relative to the BB algorithm. To test whether the benefits of the GOA algorithm are observed consistently, we varied the control and prediction horizons in order to vary the complexity of the MINLP. An increase in the control horizon results in an increase in the number of binary variables. Since in this example, the six binary variables correspond to the manipulated variables, the number of binary variables is 6m, where m is the control horizon. On the other hand, an increase in the prediction horizon results in an increase in number of nonlinear constraints. The nonlinear constraints correspond to the bounds on the levels of the three tanks and equal 6p (3 upper bounds and 3 lower bounds), where p is the prediction horizon. Table 1 documents the average computation time required per MINLP solution using the BB and GOA algorithms, for a variable control horizon between 2 and 8 and a fixed prediction horizon of 8. This corresponds to an increase in the binary variables from 12 to 48 with fixed number of constraints, 48. It can be seen that the computational time for both GOA and BB typically increases with increasing number of binary variables. However, the GOA algorithm scales well with the increase in binary variables and outperforms the BB algorithm in all the cases described in Table 1. Table 1 also documents the average computation time per MINLP solution for a variable prediction horizon between 3 and 8 while keeping the number of binary variables constant at 12 in all cases. This corresponds to a variation in nonlinear constraints between 18 and 48. This results in an increase in the computational effort for both these
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
5773
Table 1. Speedup Factor and Average Computation Time Using BB (tBB) and GOA (tGOA) for Various Cases p m nonlinear constraints binary variables tGOA (s) tBB (s) 8 2 3 4 5 6 7 8
48
12 18 24 30 36 42 48
0.17 0.31 0.32 0.27 0.37 0.35 0.40
3 2 4 5 6 7 8
18 24 30 36 42 48
12
0.07 0.10 0.12 0.13 0.16 0.17
S
5.44 31.46 16.65 53.35 44.44 137.59 106.52 388.75 236.65 637.86 369.93 1053.92 565.88 1428.99 0.72 1.46 1.99 4.50 4.16 5.44
9.75 14.73 16.47 36.04 26.34 31.46
techniques. However, the increase in the computational time for the BB algorithm is more significant as compared to the GOA algorithm. This can be explained by the fact that in the BB algorithm, the increase in the number of nonlinear constraints increases the computational effort for solving the relaxed NLP. On the other hand in the GOA algorithm, only active nonlinear constraints are linearized leading to a smaller MILP subproblem. Moreover, the nonlinear constraints become linear upon fixing the binary variables in the primal problem, and the resulting QP solution is relatively cheap. Further, the use of analytical gradients for nonlinear constraints34 and the objective function 29 make the GOA computations efficient. We use a speed up factor to quantify the computational benefits of GOA defined as follows: S)
Figure 5. Schematic diagram of the hydraulic process plant for temperature and level control in tank 1 and tank 3 (adapted from refs 6 and 13).
average cpu time required for BB algorithm average cpu time required for GOA algorithm (42)
Thus, S determines the increase in the computational efficiency and has been computed for various cases in Table 1. The smallest speedup value indicates an order of magnitude increase in computational efficiency of the GOA algorithm over the BB algorithm while the largest speedup value indicates a 3 orders of magnitude increase in the efficiency of GOA. It is clear that the GOA algorithm is especially suitable for real time implementation of MPL-model-based predictive control of NHDS. It is of interest to investigate the performance and computational attributes of the MLD-based MPC algorithm for the problem discussed above. The advantage of using an MLD model lies in the fact that the resulting control problem is an MIQP with a guaranteed globally optimal solution. However, the disadvantage is the larger problem size as discussed previously and therefore increased computational time. The MIQP resulting from the MPC problem using the MLD model with prediction and control horizons of 5 and 2, respectively, consists of 112 binary variables, 1189 continuous variables, and 5919 linear inequalities. The equivalent MINLP resulting from the MPC problem using the MPL model consists of 12 binary variables, 4 continuous variables, and 62 inequalities (30 nonlinear and 32 linear). It was observed that for several samples, the time needed for computing the control moves using the MIQP formulation solved with CPLEX, was in excess of an hour and thus impractical. It has also been shown in literature that the computational complexity of the BB algorithm implemented in CPLEX is exponential with increase in the prediction horizon and is therefore not scalable.26 Thomas et al.27 also report large computation times for the predictive control of a cylindrical 3-tank system modeled as an MLD system despite
Figure 6. Control of temperature (T3) and level (h3) in tank 3 and level (h1) in tank 1 of the hydraulic process plant for various set-point changes. The cyclic behavior is due to the semicontinuous nature of the process.
making simplifying assumptions. They report a mean computation time of 36 s for solution of a MIQP with 18 continuous variables, 14 binary variables and 104 constraints using the BB strategy. On the other hand, the mean computation time for the MINLP using the GOA approach as described above was only 0.12 s. An implementation of the GOA algorithm on an experimental setup has been described in ref 28. Despite these computational advantages, the MINLP/GOA approach may be susceptible to local optima and must be therefore used with caution. Ultimately, the trade-off is between an implementable but possibly suboptimal solution on the one hand and an impractical but global solution on the other. 4.2. A PLC-Based Composite Discrete/Continuous Control of a Hydraulic Process. A schematic of the hydraulic process plant, whose objective is supply of warm water, is
5774
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
shown in Figure 5. Tank 1 receives cold water with a flow rate q0 via an on/off pump and feeds water to tank 21 and tank 22 at rates q11 and q12, respectively, using on/off pumps. Water is heated in tank 21 and tank 22 by two separate on/off electric heaters namely, QH1 and QH2, respectively. The warm water is then fed from these two tanks to tank 3 by means of two separate on/off pumps with flow rates q21 and q22. Finally, the warm water from tank 3 discharges at a constant rate (q32 ) 1) as the final product. One control valve is used for continuous manipulation of a recycle flow (q31) from tank 3 to tank 1 and another for fresh cold water to tank 3 (q23). A first principles model has been provided by Colmenares et al.13 and is reproduced below. dh1 ) 0.1154(q31 + q0 - q11 - q12) dt
(43)
dh3 ) 0.1154(q21 + q22 + q23 - q31 - q32) dt
(44)
h1 h3
dT1 ) 0.1154(q31(T3 - T1) + q0(10 - T1)) dt
dT3 ) 0.1154(q21(T21 - T3) + q22(T22 - T3) + dt q23(10 - T3)) dh21 ) 0.1346(q11 - q21) dt h21
dT21 ) 0.1346(7.102QH1 + q11(T1 - T21)) dt dh22 ) 0.1346(q12 - q22) dt
h22
dT22 ) 0.1346(7.102QH2 + q12(T1 - T22)) dt
(45)
(46) (47)
(48)
(49)
(50)
where q0, q11, q12, q21, q22, and q32 represent flow rates that have been normalized to take a value of 1 when they are active and 0 otherwise (that is, they are discrete inputs). Flow rates q23 and q31 represent continuous normalized values and the rate of heat addition by the heaters namely QH1 and QH2 take discrete values of 0 or 1. The control objective of this process is to obtain water at a desired temperature while maintaining levels in tank 1 and tank 3 within prescribed upper and lower limits. Each of the two tanks must fill up to a height of 15 cm and be heated to a temperature of 25 °C starting from 10 °C before they can discharge to tank 3. Simultaneous filling and discharging of a particular tank is not permitted. Further, during the discharge of a tank, the corresponding heater is turned off. Also, tank 21 and tank 22 cannot be simultaneously discharged but simultaneous filling is allowed. However, once a particular tank begins discharging it cannot be interrupted until the water level drops to 5 cm. This process is an example of composite discrete/ continuous processes, whose control is achieved by integration of a PLC with MPC. The PLC determines the filling (q11, q12) and heating (QH1, QH2) in tank 21 and tank 22 using discrete decisions. The MPC controller has the task of controlling levels in tank 1 (H1) and tank 3 (H3) as well as temperature in tank 3 (T3). It accomplishes these tasks by continuous manipulation of flow rates q31 and q23 and making discrete decisions on the status of q0, q21, and q22. Depending on the set points that need
Figure 7. Levels (h21 and h22) and temperature (T21 and T22) in tank 21 and tank 22 of the hydraulic process plant for control of temperature and level in tank 3 and level in tank 1. (For clarity, simulation results corresponding to 0-2.5 h are shown from the total simulation duration of 5 h.)
to be achieved, the MPC controller decides whether tank 21 and tank 22 should discharge once they have been readied for discharge by the PLC or if they should merely wait to discharge. The hydraulic process has been modeled using the MPL model as presented in our previous work.6 It consists of eight continuous states (H1, H3, T1, T3, h21, T21, h22, T22), two continuous inputs(q23 and q31), and seven discrete inputs(q0, q11, q12, q21, q22, QH1, QH2) with five logical (binary) variables (δ1, δ2, δ3, δ4, δ5). Of these logical variables, variables δ1 and δ3 indicate the maximum levels in tank 21 and tank 22, respectively, while δ2 and δ4 represent the corresponding minimum levels. Variable δ5 corresponds to the permission to begin offloading either tank 21 or tank 22. As this permission applies to either of the two tanks, offloading of tank 21 is assigned a higher priority in the case of a tie. Two additional discrete states whose status determines if tank 21 and tank 22 are ready for discharge are defined. These discrete states depend on indictor variables δ1 to δ5 as well as discrete manipulated variables QH1 and QH2. The above logic results in 52 linear inequalities and these form part of the MINLP constraints. The corresponding set of expressions representing the logic as discussed above is provided in ref 6. The original nonlinear dynamic eqs 43-50 are considered as the “plant” model and the corresponding partially linear model as the controller model. Since the PLC decision-making scheme is dependent on the states of the system, it necessitates immediate changes in the status of the discrete inputs as soon as the states reach their threshold. This in turn compels selection of prediction horizon (p) equal to control horizon (m) in the MPC problem. We use a sampling time, ts ) 30 s, and prediction and control horizons of unity, which corresponds to an online optimization problem with 12 binary variables (5 state and 7 input variables), 2 continuous input variables and 16 nonlinear constraints (bounds on the 8 state variables), 52 linear constraints and 28 bound constraints (bounds on the 14 decision variables). Model predictions are corrected using constant output disturbance feedback strategy. We compare the computational efficiencies of the GOA algorithm with the BB algorithm for the composite control of the hydraulic process. The set-points in the temperature of the exit flow are implemented as follows: 17 °C during first hour, 20 °C during second hour, 16 °C for the next 2 hours, 18 °C for the fifth hour. The target levels in tank 3 and tank 1 over the entire duration is 15 cm. The simulation results using GOA
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
are shown in Figure 6. The oscillatory nature of all variables is due to the semicontinuous nature of the process characterized by the filling and emptying of tank 21 and tank 22. The water levels in tank 1 and tank 3 also oscillate about their set-point. The simulation results for the manipulated inputs consisting of the two continuous and seven binary variables are not shown for brevity. Temperatures and levels in tank 21 and tank 22 are shown in Figure 7. The control inputs calculated using the BB algorithm are found to be similar to GOA results and have been omitted for brevity. The mean and standard deviation for the computation time using GOA are 1.18 s and 1.29, respectively, and for BB they are 8.16 and 5.02 s, respectively. Thus, the computational benefits of the GOA algorithm noted in the previous example are also seen here. In summary, the GOA approach is computationally superior and scales well with the size of the MINLP. However, as discussed previously, the outer approximation using linearization may be invalid in presence of nonconvexities. In such cases, the GOA algorithm may yield a suboptimal solution or result in integer infeasibilities and needs to be used with caution.
5775
online feedback predictive control of hybrid systems. As evident from the current work, the issue of computational complexity depends on the model formulation, the optimization algorithm as well as the solver employed. In this work, we have exploited the fixed structure of MPL framework and demonstrated that the GOA algorithm is well-suited for the online solution of the MINLP resulting from the optimal feedback control of nonlinear hybrid systems. In particular, the primal problem of the GOA algorithm reduces to a comparatively simpler QP owing to the structure of the MPL model. Moreover, use of analytical gradients for the objective function and nonlinear constraints needed in the master problem of GOA also aids in enhanced efficiency. The benefits for online control of the nonlinear hybrid dynamic system are demonstrated using the three-spherical tank benchmark system and a PLC-based hydraulic process plant. A significant component of the online computational burden in the GOA algorithm is due to the repeated solution of the MILP subproblem. Multiparametric MILP techniques precompute solutions corresponding to the parameter space in the form of a look-up table.29,30 Here, the online burden consists of locating the region in the parameter space and deploying the relevant piecewise affine control law. An issue that can be further explored is integration of these multiparametric approaches with the master problem of the GOA. The reduced space approach developed in the optimization literature is relevant to online implementation of predictive control of hybrid systems and can be extended to predictive control formulations that incorporate stability and robustness.31,32
5. Conclusions Applications of hybrid systems are becoming increasingly common in the process industry. The main hurdle in the optimal control of hybrid systems is the requirement of an online solution to an MINLP/MIQP. This is probably one of the main reasons for the scarcity of literature on practical implementation of
Appendix I. Expressions for Matrices Used in Equations 14-16
[
¯ avg) (L¯kΦ ¯ avg)(L¯kΦ ¯ avg) (L¯k+1Φ
H2k )
[
¯ avg)(L¯k+1Φ ¯ avg)...(L¯kΦ ¯ avg) H1k ) (L¯k+2Φ : ¯ avg)(L¯k+p-2Φ ¯ avg)...(L¯kΦ ¯ avg) (L¯k+p-1Φ
]
(L¯kΓ¯ avg) ¯ avg)(L¯kΓ¯ avg) (L¯k+1Φ
[0] (L¯k+1Γ¯ avg)
[0]
[0]
[0]
[0]
¯ avg)(L¯k+1Φ ¯ avg)(L¯kΓ¯ avg) (L¯k+2Φ : :
¯ avg)(L¯k+1Γ¯ avg) (L¯k+2Φ : :
:
:
: : (L¯k+m-2Γ¯ avg) ¯ avg)(L¯k+m-2Γ¯ avg) (L¯k+m-1Φ
[0] : : (L¯k+m-1Γ¯ avg)
: :
: :
: :
¯ avg)(L¯k+m-1Γ¯ avg) + (L¯k+mΓ¯ avg) (L¯k+mΦ : ¯ avg)(L¯k+p-2Φ ¯ avg)...(L¯k+m-1Γ¯ avg) + (L¯k+p-1Φ
¯ avg)(L¯k+p-2Φ ¯ avg)...(L¯kΓ¯ avg) (L¯k+p-1Φ ¯ avg)(L¯k+p-2Φ ¯ avg)...(L¯k+1Γ¯ avg) (L¯k+p-1Φ ¯ avg)(L¯k+p-2Φ ¯ avg)...(L¯k+m-2Γ¯ avg) (L¯ ¯ ¯ ¯ ¯ ¯ (L¯k+p-1Φ k+p-1Φavg)(Lk+p-2Φavg)...(Lk+mΓavg) + ... + ¯ avg)(L¯k+p-2Γ¯ avg) + (L¯k+p-1Γ¯ avg) (L¯k+p-1Φ
[
(L¯kf¯avg) ¯ avg)(L¯kf¯avg)] + [(L¯k+1f¯avg)] [(L¯k+1Φ ¯ avg)(L¯k+1f¯avg)] + [(L¯k+2f¯avg)] ¯ avg)(L¯k+1Φ ¯ avg)(L¯kf¯avg)] + [(L¯k+2Φ H3k ) [(L¯k+2Φ : ¯ avg)...(L¯k+1f¯avg)] + ... + [(L¯k+p-1f¯avg)] ¯ avg)...(L¯kf¯avg)] + [(L¯k+p-1Φ ¯ avg)(L¯k+p-2Φ [(L¯k+p-1Φ
]
]
j 1k, H j 2k, and H j 3k by matrix C, j 2k, and H j 3k may be obtained by premultiplying each element block of H j 1k, H Coefficient matrices H j j j j j j respectively. Matrices, H1k-1, H2k-1, and H3k-1 can also be constructed using H1k, H2k, and H3k, respectively, as follows:
5776
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
H1k-1 )
[ ]
Inx , H1k,r
H2k-1 )
[
]
[0]nx×nuc , H2k,r
H3k-1 )
[ ] [0]nx×1 H3k,r
where H1k,r, H2k,r, and H3k,r represents matrices made from first (p - 1)nx rows of the matrices H1k, H2k, and H3k, respectively. Variables nx andnuc represent the number of states and continuous inputs, respectively and p represents prediction horizon. II. Expressions for Gradient of Matrices Needed in Equations 29-34. j 1k. 1. Gradient of H
[] ¯ 1k ∂H ∂δk ¯ 1k ∂H
L¯k+i )
[
(1 - δ1,k+i)(1 - δ2,k+i)Inx
(δ1,k+i)(1 - δ2,k+i)Inx
∂L¯k+i ) ∂δk+i -(1 - δ2,k+i)Inx -(δ2,k+i)Inx (1 - δ2,k+i)Inx (δ2,k+i)Inx -(1 - δ1,k+i)Inx (1 - δ1,k+i)Inx -(δ1,k+i)Inx (δ1,k+i)Inx
[
[ ] (
¯ avg) C(L¯k+1Φ
∂L¯k ¯ Φ ∂δk avg
)
¯ 1k ∂H ) ∂L¯k ∂δk ¯ avg)(L¯k+1Φ ¯ avg) ¯ C(L¯k+2Φ Φ ∂δk avg : ∂L¯k ¯ avg)(L¯k+p-2Φ ¯ avg)... ¯ Φ C(L¯k+p-1Φ ∂δk avg
(
[]
¯ 2k ∂H ) ∂δk+1 ∂δ¯ k : ¯ 2k ∂H
∂δk+p-1
)
j 2k with respect to δk+i, i ) 0, 1, ..., The partial derivative of H p - 1 may be written as follows:
[ [ ] (
( )
)
¯ 2k ∂H ∂δk+i
[0] : [0](ithrow) ∂L¯k+i ¯ ¯ )...(L¯kΦ ¯ avg) C Φ (L¯ Φ ∂δk+i avg k+i-1 avg ∂L¯k+i ¯ avg) ¯ ) C(L¯k+i+1Φ Φ ∂δk+i avg ¯ avg)...(L¯kΦ ¯ avg) (L¯k+i-1Φ : ∂L¯k+i ¯ avg)(L¯k+p-2Φ ¯ avg)... ¯ Φ C(L¯k+p-1Φ ∂δk+i avg ¯ avg)...(L¯kΦ ¯ avg) (L¯k+i-1Φ
(
)
(
)
(
,
)
Vector, δk+i, i ) 1, 2, ..., p - 1 is a vector of binary variables at (k + i)th time step which can be written as follows: δk+l ) [δ1,k+i δ2,k+i ... δnb,k+i ]
[0] : [0](ithrow) ∂L¯k+i ¯ ¯ )...(L¯k+j-1Γ¯ avg) C Φ (L¯ Φ ∂δk+i avg k+i-1 avg ∂L¯k+i ¯ avg) ¯ ¯ )...(L¯k+j-1Γ¯ avg) Φ (L¯ Φ C(L¯k+i+1Φ ∂δk+i avg k+i-1 avg : ∂L¯k+i ¯ avg)...(L¯k+i+1Φ ¯ avg) ¯ Φ C(L¯k+p-1Φ ∂δk+i avg
(
)
(
)
(
)
[ ] ¯ avg)...(L¯k+j-1Γ¯ avg) (L¯k+i-1Φ
where i ) 0, 1, ..., p - 1; j ) 1, 2, ..., i, and j e m - 1.
i ) 1...p - 1
T
)
jthcol
j k+i, i ) 0, 2, ..., p - 1 is function of only δk+i. Note that L j Therefore, ∂Lk+i/∂δk+j ) 0,∀i * j. j 1k with respect to δk+i, i ) 1,2, ..., Thus, partial derivative of H p - 1 can be written as follows:
¯ 1k ∂H ∂δk+i
]
¯ 2k ∂H ∂δk ¯ 2k ∂H
j 1k with respect to δk results: Partial derivative of H
)
]
where Inx is an identity matrix of the size of number of states (that is nx). Thus
∂δk+p-1
(
(δ1,k+i)(δ2,k+i)Inx
The above derivative follows the same structure for any number of binary variables. Therefore, a generalized computer program can be developed, using the symbolic toolbox of MATLAB that produces the above derivatives. j 2k. 2. Gradient of H
¯ 1k ∂H ) ∂δk+1 ∂δ¯ k : ¯ 1k ∂H
∂L¯k ¯ C Φ ∂δk avg
(1 - δ1,k+i)(δ2,k+i)Inx
(61)
j k+i is a function of above binary variables. Construction and L j of Lk+i for a two binary variables systems, δk+i ) [δ1,k+i δ2,k+i]T is as follows:
( ) ¯ 2k ∂H ∂δk+i
) jthcol
[0] : [0](ithrow) ∂L¯k+i C Γ¯ ∂δk+i avg
(
)
(
¯ avg) C(L¯k+i+1Φ :
∂L¯k+i Γ¯ ∂δk+i avg
)
(
¯ avg)...(L¯k+i+1Φ ¯ avg) C(L¯k+p-1Φ
∂L¯k+i Γ¯ ∂δk+i avg
where i ) 0, 1, ..., p - 1; j ) i + 1, and j e m - 1.
)
]
[]
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
( ) ¯ 2k ∂H ∂δk+i
¯ 3k ∂H ∂δk ¯ 3k ∂H ¯ 3k ∂H ) ∂δk+1 ∂δ¯ k : ¯ 3k ∂H
) [0] jthcol
where i ) 0, 1, ..., p - 1; j > i + 1, and j e m - 1.
( ) ¯ 2k ∂H ∂δk+i
[ ]
) [0]
∂δk+p-1
mthcol
( )
where i ) 0, 1, ..., m - 2.
(
)
¯ 2k ∂H ∂δk+m-1
[
C
) mthcol
[0] : [0](m-1throw) ∂L¯k+m-1 C Γ¯ ∂δk+m-1 avg ∂L¯k+m-1 ¯ avg) C(L¯k+mΦ Γ¯ ∂δk+m-1 avg :
[
( ) ¯ 2k ∂H ∂δk+m
)
)
(
(
)
(
)
∂L¯k+m-1 Γ¯ ∂δk+m-1 avg
( ) ( )
( )
¯ avg)(L¯k+p-2Φ ¯ avg)... C(L¯k+p-1Φ
∂L¯k+m Γ¯ ∂δk+m avg
Similarly we can write
[
)
mthcol
(
)
)
∂L¯k+p-1 ¯ )...(L¯k+mΓ¯ avg) + ¯ (L¯ Φ Φ ∂δk+p-1 avg k+p-2 avg
... + C
¯ 3k ∂H ) ∂δk+1
( )
( )
( )
( ) ( ) ( )
∂L¯k+1 f¯ ∂δk+1 avg ∂L¯k+1 ¯ avg)(L¯k+p-2Φ ¯ avg)... C(L¯k+p-1Φ ∂δk+1 ¯ avg) C(L¯k+2Φ
∂L¯k+1 ¯Φavg(L¯k+1f¯avg) + C(L¯k+p-1Φ ¯ avg)... f¯ ∂δk+1 avg
and
¯ 3k ∂H ) ∂δk+p-1
[0] [0] (p - 1)throw ∂L¯k+p-1 ¯ (L¯ ¯ )...(L¯kf¯avg) + C Φ Φ ∂δk+p-11 avg k+p-2 avg
(
C
) ( )
]
∂L¯k+p-1 ¯ (L¯ ¯ )...(L¯k+1f¯avg) + Φ Φ ∂δk+p-1 avg k+p-2 avg
(
... + C
)
∂L¯k+p-1 f¯ ∂δk+p-1 avg
]
Literature Cited
[0] : [0](p-1thraw) ∂L¯k+p-1 ¯ )...(L¯k+m-1Γ¯ avg) + ¯ (L¯ Φ C Φ ∂δk+p-1 avg k+p-2 avg C
( ) ( )
Similarly derivatives ∂H1k-1/∂δk, ∂G2k-1/∂δk, ∂H3k-1/∂δk can be written.
)
(
∂L¯k f ∂δk avg
∂L¯k+1 ∂L¯k+1 ¯ avg(L¯k+1f¯avg) + C Φ f¯ ∂δk+1 ∂δk+1 avg ∂L¯k+1 ¯ avg) ¯ (L¯ f¯ ) + C(L¯k+2Φ Φ ∂δk+1 avg k+1 avg
C
mthcol
¯ 2k ∂H ∂δk+p-1
( )
¯ avg) C(L¯k+1Φ
[0]
[0] : [0](mthrow) ∂L¯k+m ∂L¯k+m ¯ avg(L¯k+m-1Γ¯ avg) + C C Φ Γ¯ ∂δk+m ∂δk+m avg : ∂L¯k+m ¯ avg)(L¯k+p-2Φ ¯ avg)... ¯ (L¯ Φ Γ¯ ) + C(L¯k+p-1Φ ∂δk+m avg k+m-1 avg
(
] [ ] [ ]
∂L¯k f ∂δk avg
¯ 3k ∂H ) ∂L¯k ∂δk ¯ avg)(L¯k+1Φ ¯ avg) C(L¯k+2Φ f ∂δk avg : ∂L¯k ¯ avg)(L¯k+p-2Φ ¯ avg)... C(L¯k+p-1Φ f ∂δk avg
¯ avg)(L¯k+p-2Φ ¯ avg)... C(L¯k+p-1Φ
( )
5777
(
)
(
)
∂L¯k+p-1 ∂L¯k+p-1 ¯ avg(L¯k+p-2Φ ¯ avg)... Φ Γ¯ ∂δk+p-1 ∂δk+p-1 avg
j 3k. 2. Gradient of H
(1) Branicky, M. S.; Borkar, V. S.; Mitter, S. K. A unified framework for hybrid control: model and optimal control theory. IEEE Transa. Automat. Control 1998, 43 (1), 31. (2) Bemporad, A.; Morari, M. Control of systems integrating logic, dynamics, and constraints. Automatica 1999, 35 (3), 407. (3) Buss, M.; Glocker, M.; Hardt, M.; von Stryk, O.; Bulirsch, R.; Schmidt, G. Nonlinear hybrid dynamical systems: modeling, optimal control, and applications. In Modelling, Analysis, and Design of Hybrid Systems; Sebastian, E.; Goran, F.; Schnieder, E., Eds.; Springer-Verlag: New York, 2002; p 311. (4) van der Schaft, A.; Schumacher, H. An Introduction to Hybrid Dynamical Systems; Springer-Verlag: London, 2000. (5) Kowalewski, S., Introduction to the analysis and verification of hybrid systems. In Modelling, Analysis, and Design of Hybrid Systems, Sebastian, E.; Goran, F.; Schnieder, E., Eds.; Springer-Verlag: New York, 2002; 153.
5778
Ind. Eng. Chem. Res., Vol. 48, No. 12, 2009
(6) Nandola, N. N.; Bhartiya, S. A multiple model approach for predictive control of nonlinear hybrid systems. J. Process Control 2008, 18 (2), 131. (7) Alur, R.; Courcoubetis, C.; Halbwachs, N.; Henzinger, T. A.; Ho, P.-H.; Nicollin, X.; Olivero, A.; Sifakis, J.; Yovine, S. The algorithmic analysis of hybrid systems. Theor. Comput. Sci. 1995, 138 (1), 3. (8) Alur, R.; Dill, D. L. A theory of timed automata. Theor. Comput. Sci. 1994, 126 (2), 183. (9) Sontag, E. D. Nonlinear regulation: the piecewise linear approach. IEEE Trans. Automat. Control 1981, AC-26 (2), 346. (10) Heemels, W. P. M. H.; De Schutter, B.; Bemporad, A. Equivalence of hybrid dynamical models. Automatica 2001, 37 (7), 1085. (11) Williams, H. P. Model Building in Mathematical Programming, 3rd ed.; John Wiley & Sons: New York, 1993. (12) Raman, R.; Grossmann, I. E. Relation between MILP modelling and logical inference for chemical process synthesis. Comput. Chem. Eng. 1991, 15 (2), 73. (13) Colmenares, W.; Cristea, S.; De Prada, C.; Villegas, T. MLD systems: modeling and control. Experience with a pilot process, Proceedings of the 2001 IEEE International Conference on Control Applications (CCA’01), 5-7 Sept. 2001; IEEE: Mexico City, Mexico, 2001, 618. (14) Ozkan, L.; Kothare, M. V.; Georgakis, C. Model predictive control of nonlinear systems using piecewise linear models. Comput. Chem. Eng. 2000, 24 (2), 793. (15) Schott, K. D.; Bequette, B. W., Multiple model adaptive control. In Multiple Model Approaches to Modelling and Control; Taylor & Francis Inc.: Oxford, U.K., 1997; p 269. (16) Lazar, M.; Heemels, W. P. M. H.; Weiland, S.; Bemporad, A. Stabilizing Model Predictive Control of Hybrid Systems. IEEE Trans. Automat. Control 2006, 51 (11), 1813. (17) Grossmann, I. E.; Kravanja, Z. Mixed-integer nonlinear programming techniques for process systems engineering. Comput. Chem. Eng. 1995, 19 (1), 189. (18) Geoffrion, A. M. Generalized benders decomposition. J. Optim. Theor. Appl. 1972, 10 (4), 237. (19) Floudas, C. A. Nonlinear and Mixed-Integer Optimization: Fundamentals and Applications; Oxford University Press: New York, 1995. (20) Fletcher, R.; Leyffer, S. Solving Mixed Integer Nonlinear Programs by Outer Approximation. Math. Program. 1994, 66, 327. (21) Bansal, V.; Sakizlis, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. New algorithms for mixed-integer dynamic optimization. Comput. Chem. Eng. 2003, 27 (5), 647.
(22) Jackson, J. R.; Grossmann, I. E. A disjunctive programming approach for the optimal design of reactive distillation columns. Comput. Chem. Eng. 2001, 25 (11-12), 1661. (23) Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14 (7), 769. (24) Porn, R.; Harjunkoski, I.; Westerlund, T. Convexification of different classes of non-convex MINLP problems. Comput. Chem. Eng. 1999, 23 (3), 439. (25) Villa, J. L.; Duque, M.; Gauthier, A.; Rakoto-Ravalontsalama, N., A New Algorithm for Translating MLD Systems into PWA Systems, Proceedings of the 2004 American Control Conference (ACC), Jun 30-Jul 2 2004, Institute of Electrical and Electronics Engineers Inc.: Piscataway, NJ, 2004; p 1208. (26) Till, J.; Engell, S.; Panek, S.; Stursberg, O. Applied hybrid system optimization: An empirical investigation of complexity. Control Eng. Pract. 2004, 12 (10), 1291. (27) Thomas, J.; Dumur, D.; Buisson, J. PredictiVe Control of Hybrid Systems under a Multi-MLD Formalism with State Space Polyhedral Partition, Proceedings of the 2004 American Control Conference, 30 June-2 July 2004, IEEE: Boston, MA, 2004; p 2516. (28) Nandola, N. N.; Bhartiya, S. Hybrid system identification using a structural approach and its model based control: An experimental validation. Nonlinear Anal., Hybrid Syst., 2009, 3, 87. (29) Sakizlis, V.; Dua, V.; Perkins, J. D.; Pistikopoulos, E. N. The Explicit Control Law for Hybrid Systems Via Parametric Programing, American Control Conference, Anchorage AK, 2002, 674. (30) Dua, V.; Nikolaos, A. B.; Pistikopoulos, E. N. A multiparametric programming approach for mixed integer quadratic engineering problems. Comput. Chem. Eng. 2002, 26, 715. (31) Mhaskar, P.; El-Farra, N. H.; Christofides, P. D. Robust predictive control of switched systems: Satisfying uncertain schedules subject to state and control constraints. Int. J. Adapt. Control Signal Process. 2008, 22 (2), 161. (32) Mhaskar, P. Robust model predictive control design for fault-tolerant control of process systems. Ind. Eng. Chem. Res. 2006, 45 (25), 8565.
ReceiVed for reView September 13, 2008 ReVised manuscript receiVed April 23, 2009 Accepted April 28, 2009 IE801384Y
