Ind. Eng. Chem. Res. 2009, 48, 1933–1943
1933
PROCESS DESIGN AND CONTROL Integrated Design and Control Using a Simultaneous Mixed-Integer Dynamic Optimization Approach Rodrigo Lo´pez-Negrete de la Fuente and Antonio Flores-Tlacuahuac* Departamento de Ingenierı´a y Ciencias Quı´micas, UniVersidad Iberoamericana, Prolongacio´n Paseo de la Reforma 880, Me´xico DF, 01210, Me´xico
In this work, we address the simultaneous design and control of a binary distillation column. The problem is first formulated as a mixed-integer dynamic optimization problem that is then transformed into a mixedinteger nonlinear programming problem using the simultaneous dynamic optimization approach (i.e., full discretization). This formulation is capable of designing the optimal feed tray location, tray sizing, optimal operating steady states, the optimal open-loop trajectory between them, and also the best controller paring and parameters that does the best tracking of the open-loop trajectory. For solving the complex mixedinteger dynamic optimization problem, an optimization decomposition strategy is proposed. The solution strategy is based on solving relaxed versions of the optimization problem and using the results to initialize complex problem versions. The full space problem was solved with the Bonmin solver. Two cases were analyzed, and in both the solver and the proposed decomposition strategy were capable of solving the problems successfully. 1. Introduction Traditionally process design and control problems have been addressed sequentially.1 However, the need to approach process design and control (PDC) problems simultaneously has long been recognized.2-5 The natural separation between design and control problems happens because design problems are associated with steady states, whereas control problems are treated as dynamic processes. In sequential methods, usually first the steady-state design is accomplished, which in turn is then highly optimized (reducing degrees of freedom for process control), and then a control scheme is designed to make sure the process does not deviate from the previously designed operating point. In general, because processes have a nonlinear nature, the sequential approach for design and control tends to embed operability problems in the resulting process. For example, the process could end up being operated in a highly sensitive region (e.g., Lo´pez-Negrete et al.6) or regions where other multiplicities might arise such as input/output multiplicities, sustained oscillations, chaotic behavior, and so forth.7-9 Thus, in the end, the process might be operated in a suboptimal region because modifications were needed to achieve good operability properties. The simultaneous approach to solve PDC problems has a great advantage over the sequential method because operability constraints may be added to avoid certain undesired operation regions.10 In the present work, the proposed simultaneous approach requires that the problem be stated as an optimization problem in which continuous and binary variables occur (a mixed-integer nonlinear problem or MINLP). The continuous variables are associated with process state profiles and manipulated variables, whereas binary variables are used to model logical decisions such as whether to use a certain piece of equipment or to choose between different possible process * To whom correspondence should be addressed. E-mail:
[email protected]. Phone/Fax: (+52) 55 59504074. http:// 200.13.98.241/∼antonio.
structures (feed tray locations in distillation columns, possible controlled-manipulated variables pairings, etc.). In the literature, some papers have been published that discuss the simultaneous design and control of distillation columns.4,11,12 In their work, they successfully accomplish the design and control of the addressed columns using simple SISO PI or parametric controllers. They apply the sequential approach for dynamic optimization in which they discretize the control variables and integrate the DAE model at every iteration of the optimization process. Unfortunately, this methodology might be at a disadvantage because it cannot evaluate transitions where one (or all) of the steady states considered is (are) unstable.13 Some other recent attempts to approach the design and control problem using a dynamic flexibility analysis have been recently reported.14 In this work, an optimization formulation is proposed to simultaneously address design and control problems. The problem is cast in terms of a mixed-integer dynamic optimization (MIDO) problem. The contributions of the present work are: (1) The MIDO problem is solved using a solution strategy that up to our best knowledge has not been reported for the purposes of the work and (2) a way to decompose and solve MIDO problems that works efficiently is proposed. This strategy is based on solving relaxed versions of the optimization problem and using the results to initialize complex problem versions. For testing the optimization formulation and the MIDO solution strategy,thedesignandcontrol(SDC)problemofamethanol-water distillation column is addressed. The MIDO solution methodology requires that both the control and state variables are discretized, and orthogonal collocation methods are used to approximate the differentials at each collocation point. This method has great advantages over the sequential method, for example unstable steady states can be used as initial or final points in the transition.15 It is an extension on our previous work related to optimal operation distillation columns (in that case with chemical reaction).16 In this approach, purity specifications
10.1021/ie801353c CCC: $40.75 2009 American Chemical Society Published on Web 01/27/2009
1934 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
are specified for the product streams, as well as other numerical parameters (detailed below), and using a detailed tray-by-tray model the SDC is done. The point is to calculate the optimal steady-state design that meets the purity specifications (at start and end points) and has the optimal feed tray location, the optimal transition between them, and the PI controller parameters that does the best tracking of the optimal control trajectory. Because our formulation involves continuous and binary variables, the underlying design and control problem during steady-state transitions is a MIDO problem. Solving this type of problems using the simultaneous approach13 demands good bounds and initialization schemes. Also, a robust MINLP solver for large-scale problems is required. Recently, the Bonmin17 (basic open-source nonlinear mixed-integer programming) solver, that uses Ipopt18 (an interior point optimizer) as a nonlinear problem (NLP) solver, has been released. In this work, both solvers are used to obtain the final solutions of the SDC problems and for the initialization schemes also described here. 2. Problem Statement Given are: (a) the desired composition specifications of the distillate and bottoms flow rates at the start and the end of the transition, (b) the total number of trays, (c) the total feed streamflow rate, (d) its composition, and (e) upper and lower bounds for the minimum molar hold-up in each tray. The problem is to compute the optimal steady states at the start and end points, that is, the compositions, temperatures, molar hold-ups, liquid and vapor flow rates for all stages, as well as the best design composed by the minimum hold up per tray, the reboiler heat duty, optimal feed tray location, and the condenser split fraction (CSF also defined as variable R) at each point. Also, the optimal open-loop transition between these steady states will be simultaneously computed, as well as the best control pairing and controller parameters that can carry out the best tracking of the optimal trajectory. The manipulated variables are CSF and the reboiler thermal duty QR, whereas the controlled variables are the methanol mol fraction in the distillate stream and the water mol fraction in the bottoms’ stream. Therefore, the control system is a two-point control structure where simultaneous closed-loop composition control of both product streams is imposed. The problem is then formulated as a mixed-integer nonlinear program (MINLP) and solved using the Bonmin and Ipopt packages. Although there are simple design methods (i.e., McCabeThiele method) for feed tray location in binary distillation, the interplay between all of the addressed decision variables could render different tray locations that those suggested using shortcut methods. We cannot state in advance, before solving the design problem, whether or not fixing the feed tray location will render a distillation column with better operability characteristics. This is the reason why we decided to consider the feed tray location as a decision variable, and also because we suspect that proper feed tray location could improve column operability. We would like to highlight that our purpose is to exploit the relationship between design and control to improve process operability, and that probably this objective could not be achieved if we decide to fix some parts related to column design such as feed tray location. Moreover, even when the methanol-water distillation column is a well-known example of a binary separation that has been widely used for both design and control purposes, our aim in the article is to propose a simultaneous design and control methodology able to deal with different separation systems (i.e., to exploit the natural interaction between design and control for improving process oper-
ability) and not necessarily to provide new results beyond those reported regarding the addressed separation system. 3. Dynamic Mathematical Model In this section, the mathematical model that drives the dynamic behavior of a distillation column is described. This model consists in the tray-by-tray application of the MESH equations which are mass balances (M), thermodynamic equilibrium (E), sum composition (S), and enthalpy balances (H). The modeling assumptions made are: negligible vapor-phase dynamics, nonconstant molar holdups, liquid flow modeled by the Francis Weir equation, thermodynamic equilibrium between phases with a nonideal liquid phase modeled by the NRTL equation, and a total condenser. Because trays are counted from the bottom to the top, then the model is represented by eqs 1-8, where S and C are the set of stages (including reboiler and condenser) and the set of components, respectively. All symbols are defined in the nomenclature section. • Total and component mass balances dMk ) Fk + Lk+1 + Vk-1 - Lk - Vk, k ∈ S dt Mk
(1)
dxi,k ) Fk(zi,k - xi,k) + Lk+1(xi,k+1 - xi,k) + dt Vk-1(xi,k-1ki,k-1 - xi,k) + Vkxi,k(1 - ki,k), k ∈ S, i ∈ C (2)
• Energy balance Mk
dHLk L V ) Fk(HFk - HLk ) + Lk+1(Hk+1 - HLk ) + Vk-1(Hk-1 dt HLk ) + Vk(HLk - HVk ) + Qk, k ∈ S (3)
• Thermodynamic equilibrium and NRTL equation
[(
ln γi,k ) x2j,k τji,k
γi,kPsat i,k , k ∈ S, i ∈ C Ki,k ) P Gji,k xi,k + xj,kGji,k
)
2
+
Gij,kτji,k
(4)
]
, (xj,k + xi,kGij,k)2 k ∈ S, i, j ∈ C, and j * i (5) • Sum equations (because yi,k ) xi,kKi,k) NC
NC
∑x
i,k -
i
∑x
i,kKi,k ) 0, k ∈ S, i ∈ C
(6)
i
• Francis Weir equation (liquid flow) 1.5 Lk)τk(Mk - M min k ) ,k∈S
(7)
• Condenser liquid split fraction CSF )
D L0 + D
(8)
where, NC
HLk )
∑x
l i,khi,k(Tk), k ∈ S, i ∈ C
(9)
i
NC
HVk )
∑x
v i,kKi,khi,k(Tk), k ∈ S, i ∈ C
(10)
i
The condenser split fraction (eq 8) is used instead of the reflux ratio because it has finite lower and upper bounds, that is 0
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1935
e CSF e 1. This property makes it better suited for the optimization formulation because it has tighter bounds. Eq 7 representing liquid flow dynamics was taken from Aditya and Daoutidis.19 The dynamic model of the distillation column described above comprises a set of differential and algebraic equations (DAE) system featuring a high-index problem. Before this model can be used for modeling, optimization, and control purposes, the index ought to be reduced down to an index one DAE system. The DAE high-index reduction procedure is described in the appendix section. 4. Optimization Formulation The optimal open-loop transition between steady states known a priori for a reactive distillation column has been previously addressed.16 Because in that work column design optimal steady states were known, these were not considered as decision variables during the optimization formulation. In the present work, we extend the results reported in Lo´pez-Negrete and Flores-Tlacuahuac.16 We assume that the only known information is product specification at the start and end transition points, and, in consequence, design and steady-state equation constraints must be added to the problem. Moreover, controller design and closed-loop tracking of the controlled variables will be addressed simultaneously. Therefore, a closed-loop dynamic model must be added to the constraints set. Because optimal solutions from MIDO problems are known to be difficult to achieve, we propose a decomposition strategy based on solving complex versions of the MIDO problem by progressively relaxing some assumptions. We have found that this MIDO solution strategy works reasonably well for the addressed optimization problems. Of course, more formal optimization decomposition strategies20,21 can also be used for the purpose of solving complex MIDO problems. 4.1. Objective Function. The objective function used (eq 11) expresses the aim to carry out the transition as fast as possible, therefore minimizing the production of off-spec materials invariably formed during product transitions. The first two terms take into account the differential variables in the purely open-loop optimal control dynamic model, and the following terms are the manipulated variables that are added to the objective function to accomplish smooth transitions. The decision variables are the time behavior of the system states (differential and algebraic), the manipulated variables, design variables, the feed tray location, and control structure. The terms for the molar hold-ups and for the reboiler heat duty are scaled so that weight parameters are used to highlight the importance of them. Hence, the objective function reads as follows:
min
xoc,u,xss,xcl,Y 0 NT
∑ k
{
NC NT
∫ ∑ ∑ R [x tf
[
RM k
x i,k
i
OC ss 2 i,k (t) - xi,k,end +
k
ss Moc k (t) - Mk,end
MScal k
]
]
2 2 + Rcfs[CSFOC(t) - CSFSS end] +
R
QR
[
ss Qoc R (t) - QRend
Qscal R
]} 2
dt (11)
where, as mentioned previously, the decision variables are all the open-loop system states xoc ) [Mkocxi, kocTkoc], manipulated variables u ) [CSF,QR], steady-state design variables xss ) [Mkssxi, kssTkss], closed-loop system states xcl ) [Mkclxi, kclTkcl], and
binary variables Y for the selection of the best control structure and feed tray location. The super index oc stands for variables associated to the optimal control part, ss is for steady-state variables, cl denotes variables related to closed-loop behavior, scal stands for a scaling factor, whereas the subindex end denotes values at the end of a transition. Ri,kx, RkM, and RQR are the weights of the liquid composition profile, holdups, and reboiler duty, respectively. QR is the reboiler thermal duty, Nc is the number of quadrature points, NT is the number of trays, and tf is the transition time to be implicitly minimized. Because of the presence of dynamic continuous and discrete variables, the present optimization problems turns out to be a mixed-integer dynamic optimization (MIDO) problem. 4.2. Optimal Control Formulation. When no binary decision variables are involved, the simultaneous approach for dynamic optimization transforms the dynamic optimization problem into an NLP by approximating the state and control variable profiles with a family of polynomials on finite elements.13,15 Because transition profiles can feature wide variations, in the simultaneous approach for handling dynamic optimization problems, the time horizon is divided into a series of finite elements. Within each finite element the dynamic profiles are discretized around collocation points. Continuity between adjacent finite elements is also enforced. Each one of these parts of the dynamic optimization formulation is described next. • Model discretization Ncp
o,oc + tf Hi xoc i,j ) xi
∑ k)1
Ak,Ncp
dxoc i,k , i ∈ {1..Nfe}; j, k ∈ {1..Ncp} dt (12)
Inside each finite element, the states profile xi,joc is approximated in terms of their value at the initial point of the finite element xio,oc and in terms of a collocation matrix Ak,Ncp and their derivative values Ncp. The collocation points for model approximation are normally the roots of a Lagrange polynomial of proper order. Ncp is the total number of internal collocation points, Nfe is the total number of finite elements, whereas Hi is the length of the ith finite element. • Continuity between finite elements Ncp
o,oc ) xi-1 + tf Hi-1 xo,oc i
∑ k)1
Ak,Ncp
oc dxi-1,k , i ∈ {2..Nfe}; j, k ∈ {1..Ncp} dt
(13) To enforce continuity of the states between adjacent finite elements, the value of the states at the beginning of each finite element xio,oc ought to be equal to its corresponding value in the past finite element xi -1o,oc plus the states contribution that gives rise to the states value at the end of the previous finite element. In other words, the states value at the first collocation point of the ith finite element ought to be equal to its value at the last collocation point of the previous finite element. • Approximation of the dynamic model at each collocation point dxoc i,k oc ) F(xoc i,j , d, ui,j ), i ∈ {1..Nfe}; j, k ∈ {1..Ncp} dt
(14)
The value of the derivatives of the system states (dxi,koc)/(dt) is calculated at each collocation point within each finite element
1936 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
as function of the states values xi,joc, design parameters d, and manipulated variables values ui,joc. • Initial conditions ss xo,oc 1 ) xi,k,start
(15)
The value of the states at the beginning of the first finite ss element x1o,oc is equal to the steady-state values xi,k,start from which a desired transition trajectory takes place. Notice that in ss are also decision variables. the present formulation xi,k,start Eqs 12-15 represent the discretization of the differential equations in the dynamic model (eqs 1 and 2); although the algebraic equations and variables are also discretized. The next set of restrictions were added to make sure that the values of the changes in the manipulated variables were smooth enough, meaning that no aggressive control actions would be allowed. Therefore, only a certain percent of change is allowed between collocation points. If Radau collocation points are used and only two internal collocation points are considered, then: oc oc CSFoc i,1 - CSFi,2 e ∈manCSFi,1, ∀ i ∈ {1..Nfe} oc oc CSFi,1 - CSFi,2 g -∈manCSFoc i,1, ∀ i ∈ {1..Nfe} oc oc CSFoc i,2 - CSFi,3 e ∈manCSFi,2, ∀ i ∈ {1..Nfe} oc oc CSFoc i,3 - CSFi,3 g -∈manCSFi,2, ∀ i ∈ {1..Nfe} oc oc CSFoc i,3 - CSFi+1,1 e ∈manCSFi,3, ∀ i ∈ {1..Nfe - 1} oc oc CSFoc i,3 - CSFi+1,1 g -∈manCSFi,3, ∀ i ∈ {1..Nfe - 1} (16)
QRoci,1 - QRoci,2 e ∈manQRoci,1, ∀ i ∈ {1..Nfe} QRoci,1 - QRoci,2 g -∈manQRoci,1, QRoci,2 - QRoci,3 e ∈manQRoci,2,
ss des,dis - xCH e ∈ss xCH 4O,i 4O,i
∀ i ∈ {1..Nfe}
ss des,dis xCH - xCH g -∈ss, i ∈ {start, end} (19) 4O,i 4O,i
∀ i ∈ {1..Nfe}
xHss2O,i - xHdes,bot e ∈ss 2O,i
QRoci,3 - QRoci,3 g -∈manQRoci,2, ∀ i ∈ {1..Nfe}
xHss2O,i - xHdes,bot g -∈ss, i ∈ {start, end} (20) 2O,i
QRoci,3 - QRoci+1,1 e ∈manQRoci,3, ∀ i ∈ {1..Nfe - 1} QRoci,3 - QRoci+1,1 g -∈manQRoci,3, ∀ i ∈ {1..Nfe - 1}
(17)
The numerical value of the εman parameter defined in eqs 16 and 17 was found iteratively. First, we set an initial value and see if an optimal solution can be determined. If so, then we use our process experience and decide if the manipulated variables profile looks smooth enough. If so, then we have probably found a good εman value. Otherwise, we need to guess a new εman value and repeat this procedure. As noticed, the set of constraints given by eqs 16 and 17 dictates smooth control actions between adjacent collocation points within each finite element. 4.3. Steady-State Design Formulation. The constraints of this part allow us to compute the steady-state optimal design conditions. For doing so, we only need to specify the required product purity, whereas the optimization formulation will compute the steady-state operating conditions. These constraints have the following form: ss G(xss i , d, ui ) ) 0, i ∈ {start, end}
Figure 1. Optimal feed tray location superstructure.
(18)
where the decision variables of this part are the optimal steadystate values of the states xss i , the optimal values of the design variables d, and the values of the manipulated variables state values of the states uss i . The set {start, end} refers to the start point and the end point each of which is a steady state. The constraints symbolized by G are the equations of the dynamic model (i.e., 1, 2, 4, 5, 6, 7, 8, and 41) with all (dX)/(dt) ) 0. In this section, the desired composition specifications of methanol at the condenser and water at the reboiler are also included. These constraints read as follows:
Eqs 19 and 20 state that the optimization formulation should compute both steady-state product compositions (as denoted by xCH4, 0ss and xH2O, 0ss) at the condenser and reboiler and that they should be computed with a precision given by ss. The specified values of the methanol mole fraction at the distillate and the water mole fraction at the bottoms are given by xCH4O,ides,dis and xH2O,ides,bot, respectively. Up to this point, all decision variables are of continuous nature. This means that the underlying optimization problem can be cast as a dynamic NLP. However, when integer decision variables are introduced the NLP will become a MIDO problem because binary decision variables will now be part of the results that the optimization formulation ought to determine. In this work, integer decision variables will be used to find the optimal location of the feed tray and to address the selection of the best control structure. The feed tray location is defined by the following set of constraints, where FT is the set of possible feed tray locations (not including the condenser and reboiler). These were taken from Viswanathan and Grossman22 and modeled after the superstructure shown in Figure 1. All of the intermediate feed streams Fk must add up to Fmax, the total amount of feed available: FT
∑F )F k
max
(21)
k
and there must only be one feed tray: FT
∑Y )1 k
k
(22)
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1937
Figure 2. Control loop selection superstructure.
where Yk stands for the binary variable associated to feed tray location. 4.4. Controller Tracking and Design Formulation. In this subsection, the control aspects of the design problem are described. The controllers to be designed are simple proportionalintegral (PI) controllers but there are two possible configurations. The first control configuration consists of controlling the composition of methanol in the distillate using the CSF, and the composition of water in the bottoms flow rate with the reboiler heat duty (Figure 2). The second possible control configuration consists of controlling the composition of methanol in the distillate flow rate with the reboiler heat duty, and the composition of water in the bottoms’ flow rate with the CSF. Although the second control configuration does not look familiar for the two-point composition control of distillation columns, we have allowed such a control structure partially to check the integrity of the optimization formulation results, meaning that what we normally would expect as optimal control configuration is the first control structure. Should we obtain the second control structure as the optimal one, it would mean that either our optimization formulation was wrongly formulated, or that, in fact, a new way to control the column is possible. Of course, this last result would demand a full and clear explanation of the observed result before accepting unfamiliar forms to approach the two-point control of the addressed distillation column. Although some other common control structures for two-point control could have been examined, we decided not to do so to maintain both the MIDO formulation and the CPU computation times as small as possible. In our experience, the proposed MIDO formulation would be able to compute optimal solutions with a larger number of binary variables at expenses of demanding larger CPU time. Unfortunately, MIDO problems
feature exponential complexity meaning that the CPU solution time tends to exhibit exponential growth with respect to the number of binary variables. In this regard, the use of decomposition optimization approaches20,21 represents a promising avenue for addressing the solution of large scale MIDO problems. The control structure decisions were modeled by the use of Big-M constraints. It is important to notice that the sign of the controller gains were known a priori because they were determined by fitting dynamic simulations to first-order transfer functions. The column dynamic model is used here again but introducing new variables defined as a closed-loop model and denoted by the super index cl. Then the discretized closed-loop system of equations is defined as follows. Because of the model discretization, continuity between finite elements, approximation of the dynamic model at each collocation point, and initial conditions’ constraints are similar to those used previously for the openloop optimal control part, we will omit the explanation of those constraints. In this case, the differential equations added by the control loops are also discretized and defined below. • Model discretization Ncp
xcli,j ) xo,cl i + tf Hi
dxcli,k , i ∈ {1..Nfe}; j, k ∈ {1..Ncp} k,Ncp dt
∑A k)1
(23) • Continuity between finite elements Ncp
o,cl xo,cl i ) xi-1 + tf
Hi-1
∑A k)1
k,Ncp
cl dxi-1,k , i ∈ {2..Nfe}; j, k ∈ {1..Ncp} dt
(24) • Approximation of the dynamic model at each collocation point dxcli,k cl ) F(xcli,j, d, ui,j ), i ∈ {1..Nfe}; j ∈ {1..Ncp} dt • Initial conditions ss xo,cl 1 ) xi,k,start
• Condenser liquid split fraction control structure
Figure 3. (a) Internal methanol and water compositions at both steady states, (b) minimum liquid hold-ups per tray for 1 f 2 transition.
(25)
(26)
1938 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 Table 1. Product Purity Specifications for Both Transitions Transition
1f2
variable/SS
Table 5. Feed Flow Rates and Associated Binary Variables for 1 f 2 Transition
2f1
start
end
start
k
end
1
2
3
4
5
6
7
ft
xCH4O,ides,dis xH2O,ides,dis
0.90 0.90
0.99 0.99
0.99 0.99
0.90 0.90
Table 2. Objective Function Weight Values for Both Transitions weight
Ri,kx
Ri,2x
Ri,8x
RkM
R2M
R8M
RQR
RCFS
value
1
106
106
1
102
102
1
1
Table 3. Problem Sizes and Solution Times
Yk 0 1 0 0 0 0 0 flow rate 0.0 4.5 × 104 0.0 0.0 0.0 4.97 × 10-11 5.69 × 10-9 [gmol/min] Table 6. Controller Parameters and Associated Binary Variables for 1 f 2 Transition i
Yicl
KC,iCSF
KI,iCSF
KC,iQR
KI,iQR
1 2
1 0
5.2561 × 101 0.0
1.5123 × 101 0.0
2.6760 × 1010 0.0
9.6883 × 1010 0.0
transition continuous variables binary variables constraints CPU (h) 1f2 2f1
33 485 33 485
9 9
34 924 34 924
Table 4. Manipulated Variables at Both Steady States for 1 f 2 Transition steady state
CSF
Qr [KJ/min]
1 2
6.5977 × 10-1 1.9405 × 10-1
1.2865 × 106 4.4120 × 106
ss dis CSF dis CSFcl(t) ) CSFstart - K CSF C,1 E CH4O(t) - K I,1 I CH4O(t) + bot CSF bot K CSF C,2 E H2O(t) + K I ,2 IH2O(t) (27)
where CSFcl(t) stands for the closed value of the condenser liquid split fraction, CSFstartss is its corresponding bias value, KC,1CSF dis is the controller gain of the CSF f xCH4Odist loop, ECH is the 4O control error signal between the set-point and measured values of xCH4Odist, KI,1CSF is the integral part of the CSF f xCH4Odist loop, ICH4Odis(t) is the integral of the ECH4Odis(t) error signal. Similarly, KC,2CSF is the controller gain of the CSF f xH2Obot loop, EHdis2O(t) is the control error signal between the set-point and measured values of xH2Obot, KI,2CSF is the integral part of the CSF f xH2Obot loop and IH2Obot(t) is the integral of the EH2Obot(t) error signal. • Reboiler thermal duty control structure QR bot dis ss + KQC,1R EHbot2O(t) + KI,1 IH2O(t) - KCQ2RECH (t) QclR (t) ) QR,start 4O QR dis KI,2 ICH4O(t)
(28)
where QRcl(t) stands for the closed-value of the reboiler thermal duty, QR,startss is its corresponding bias value, KC,1QR is the controller gain of the QR f xH2Obot loop, EH2Obot(t) is the control error signal between the set-point and measured values of xH2Obot, KI,1QR is the integral part of the QR f xH2Obot loop, IH2Obot(t) is the integral of the EH2Obot(t) error signal. Similarly, KC2QR is the controller gain of the QR f xCH4Odist loop, ECH4Odis(t) is the control error signal between the set-point and measured values of xCH4Odist, KI,2QR is the integral part of the QR f xCH4Odist loop, and ICH4Odis(t) is the integral of the ECH4Odis(t) error signal. The control error signals read as follows, dis oc cl (t) ) xCH (t) - xCH (t) ECH 4O 4O,dis 4O,dis
[ ] EHbotO(t) ) [xHoc O,bot(t) - xHcl O,bot(t)] 2
2
2
dis dICH 4O
20.10 33.47
dt
dis ) ECH (t) 4O
(31)
dIHbot2O
) EHbot2O(t) (32) dt • Maximum product purity deviation. To guarantee that good product purity is obtained at the end of the tracking transition, the following constraints are added. The maximum deviation between open-loop and closed-loop product purities is set by the ∈ trterm: oc cl oc [xCH O,dis(t) - xCH O,dis(t)] e ∈tr[xCH O,dis(t)] oc cl oc [xCH O,dis(t) - xCH O,dis(t)] g -∈tr[xCH O,dis(t)] (33) [xHoc O,bot(t) - xHcl O,bot(t)] e ∈tr[xHoc O,bot(t)] [xHoc O,bot(t) - xHcl O,bot(t)] g -∈tr[xHoc O,bot(t)] (34) 4
4
4
4
2
4
2
4
2
2
2
2
• Big-M constraints. Big-M constraints were used to enforce only one control structure among the two postulated control structures. Hence, if P stands for the set of control structures (i.e., P ) {1,2}), then the set of Big-M constraints are: CSF e Ycli MKCSF KC,i C CSF CSF CSF KI,i e Ycli MKCSF , KC,i , KI,i g 0; i ∈ P C
(35)
QR KC,i e Ycli MKQCR QR QR QR KI,i e Ycli MKQCR, KC,i , KI,i g 0; i ∈ P
(36)
where M stands for an upper bound normally found on a trial and error basis, Yicl is the binary variable used to denote the selection of a given control loop. Finally, because only one of the pairings should be chosen then, NP
∑Y
cl i )1
(37)
i
where Np is the number of control loops.
(29)
5. Solution Methodology
(30)
The simultaneous approach for solving dynamic optimization problems requires good initialization strategies. In a first solution approach, the addressed full MIDO problems were run in a single step using reasonably good guesses of the decision variables. However, after several trials no optimal solution was found. This is the main reason why the solution of the addressed full MIDO problems were decomposed into a series of individual optimization problems trying to find good initialization points to solve more complex versions of the same optimization
where, as mentioned previously, the super index oc stands for open-loop dynamic optimal tracking profiles, whereas the super index cl represents the corresponding values under closed-loop conditions. Therefore, the error signals simply represent the difference between desired and actually implemented control tracking signals. The integral part of the PI controllers is cast as system states in terms of the control error signals as follows:
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1939
Figure 4. (a) Distillate methanol composition, (b) water-methanol composition, (c) reboiler heat duty, and (d) condenser split fraction for 1 f 2 transition. Table 7. Manipulated Variables at Both Steady States for 2 f 1 Transition steady state
CSF
Qr [KJ/min]
1 2
1.9405 × 10-1 6.5977 × 10-1
4.4120 × 106 1.2865 × 106
problem. Unfortunately, most of the complex optimization problems tend to demand similar solution strategies because they show a strong dependence on the initial guess of the decision variables. Fortunately, when solving optimization problems, good guesses of such variables can always be provided relying in the physical knowledge of the addressed problem. We would like to highlight that in general the solution of MIDO problems is complex, and it may demand specific solution procedures as the one described in this section. From this point of view, the decomposition algorithm proposed in this part is a major contribution of the present work. The method followed to get an optimal solution of the proposed problem was to first calculate two steady states that meet the purity constraints in eqs 19 and 20. Once these are known, the pure optimal control problem (POC) is solved separately using the methods and algorithms described elsewhere (e.g., Kameswaran and Biegler,13 Flores-Tlacuahuac et al.,15
Biegler et al.23). To initialize the POC problem, a dynamic simulation is done first, and the initial guess is obtained by interpolating at each of the collocation points. The NLPs formed during these stages are solved using Ipopt. Afterward, the steadystate design constraints are added to the POC problem to simultaneously calculate the design and optimal control values. The initial guess is obtained in a similar fashion as in the POC problem for the dynamic variables, and the calculated steady states in the first step are used for the new added constraints. In parallel, the solution of the controller design problem or closed-loop problem (CLP) was addressed. This problem is similar to the POC, except that controller equations are added, and the controller parameters are used as degrees of freedom. At this stage, intuition is used to propose the variable paring that seems the best, and after fixing it the problem is solved. Again, a dynamic simulation is used to initialize the variables. Once both problems are solved separately (i.e., the simultaneous optimal control and design problem and the closed-loop problem), one can proceed to join them. The previous optimal solutions of both problems are used as initial guesses for this new problem: the simultaneous design and control problem (SDCP). So far, all of the solved problems were NLPs, and no decisions (binary variables) exist. The addition of these has been
1940 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
Figure 5. (a) Internal methanol and water compositions at both steady states, (b) minimum liquid hold-ups per tray for 2 f 1 transition. Table 8. Feed Flow Rates and Associated Binary Variables for 2 f 1 Transition k
1
2
3
4
5
6
7
Ykft flow rate [gmol/min]
0 0.0
1 4.5 × 104
0 0.0
0 0.0
0 0.0
0 0.0
0 3.46 × 10-8
Table 9. Controller Parameters and Associated Binary Variables for 2 f 1 Transition i
Yicl
KC,iCSF
KI,iCSF
KC,iQR
KI,iQR
1 2
1 0
6.3800 × 102 5.92 × 10-7
3.8541 × 103 0.0
7.0627 × 1010 0.0
1.0000 × 1011 0.0
delayed as much as possible because the solution of MINLPs tends to be more complicated. Now, we add the binary variables and logical decisions proposed in subsections 4.1 and 4.3. To keep the problems as simple as possible, the controller paring decision is added first (only two binary variables) using the previous optimal solution as an initial guess. At this point, the problems are transformed into MINLPs that are solved using Bonmin. Finally, after the controller decision problem is solved, the feed tray location decision is added and solved. The algorithm is summarized as follows. 1) Calculate two steady states 2) Calculate the dynamic simulation between them 3) Solve (maybe in parallel): • Open-loop POC problem (NLP) • Closed-loop problem with a tentative control structure (NLP) 4) Add steady-state design constraints to POC and solve (NLP) 5) Join the open-loop and closed-loop problems and solve (NLP) 6) Add controller paring decision and solve (MINLP) 7) Add feed tray location decision and solve (MINLP). 6. Results and Discussion The transitions were defined by the purity specifications given in Table 1 that were used in eqs 19 and 20. The values of the objective function weights in Table 2 were obtained during step 3 of the algorithm outlined in section 5. These values are used to highlight the importance of certain parts of the objective
function and are tuned manually to obtain smooth transitions between steady states in a pure optimal control problem. These values remained constant for all NLP and MINLP problems solved to obtain the results detailed below. Finally, a summary of problem sizes and solution times is shown in Table 3. Bonmin has several algorithms implemented to solve MINLP problems (details in Bonami et al.17). The problems in this work only contain very few binary variables (only 9, see Table 3); therefore, the Branch and Bound algorithm was chosen to solve the MINLP design and control problems. All problems were solved correctly by the solver, and it is clear that good starting points are required to obtain solutions in as little time as possible. In this case, it took from 20 to 33 CPU hours to obtain such solutions, and this is due both to the highly nonlinear and largescale nature of the problem at hand. Additional information pertaining to problem solution is given as follows: the allowed change between collocation points for the open-loop manipulated variables (i.e., eqs 16 and 17) used is ( 20%, the tolerance used in eqs 19 and 20 is ( 1 × 10-5, and finally, the allowed tracking error in eqs 33 and 34 is ( 1%. In both cases, only 20 finite elements were used. 6.1. 1 f 2 Transition. In this case, the column should carry out a transition between a steady that has low purity in both product streams to one in which they both should have at least 99% of methanol in the distillate stream and 99% of water in the bottoms’ stream. The steady states are defined by the values of the manipulated variables shown in Table 4. Figure 3 displays the tray by tray composition profiles of both components (part a of Figure 3) for the start and end steady states, and the minimum liquid molar hold up per tray for the designed column. All trays have been designed to feature a minimum liquid holdup that is the lower bound allowed for that variable, except that of tray 4. The fact that the values lie on the lower bound can be explained by realizing that with smaller holdups the column internal flow rates increase, and therefore, faster transitions can be achieved. Tray 3 has a larger minimum holdup in this case because the feed stream has been located on the tray below it (Table 5). In this way, the internal flow rate between tray 1 and 2 is large, and the molar holdup in tray 2 is maintained by the feed stream and smaller internal flow rate coming from tray 3. In Table 5, the values of the binary variables associated with the feed stream location are shown as
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1941
Figure 6. (a) Distillate methanol composition, (b) bottoms water composition, (c) reboiler heat duty, and (d) condenser split fraction for 2 f 1 transition.
well as the value of the selected feed stream. As mentioned before, the selected optimal feed stream location is at tray 2; hence, the value of the binary variable Y2ft is 1, and the rest of those binary variables are 0 (only one feed stream is allowed). In Table 6, the values of the binary variables associated with the closed-loop parings and the values of the designed controller parameters are shown. Because the value of Y1cl is 1, then the chosen SISO control loops are to control methanol concentration in the distillate stream by manipulating CSF, and the second one will control the water concentration in the bottoms flow rate by manipulating the reboiler heat duty. The dynamic transition and closed-loop tracking results are shown in Figure 4. Parts a and b of Figure 4 show the dynamic transition of methanol composition in the distillate stream and the water composition of the bottoms stream respectively, whereas parts c and d of Figure 4 show the transition of the QR and the CSF, respectively. It is shown that the optimal open-loop transition between the calculated steady states is done in 10 min for the distillate stream, whereas the bottoms stream completes it in 15 min. In contrast, the closed-loop tracking takes at most 20 min to complete the transition for one of the control loops. Also, the closed loops track the optimal trajectory very closely. Although it is almost imperceptible, one loop is tighter than the other, that is the methanol/CSF loop is tighter than the other
(part a of Figure 4). Having one loop tighter than the other when controlling compositions is a common feature of distillation control that has been reported in the open literature, and it is due to the negative interaction between the loops.11 6.2. 2 f 1 Transition. In this case, the designed column must carry out a transition from a high-purity to a low-purity steady state. In Table 7, the manipulated variable values that define such operating points are shown. These points correspond exactly to the ones calculated in the above example, and this is so because they are defined in the same way. In Figure 5, the tray by tray composition profiles and the minimum liquid molar holdups are shown for the new design. In this case, all values of the holdups end up in the lower bound. Again, in Table 8 it is shown that the chosen feed tray location is tray 2, and Table 9 shows the values of the designed controller parameter values. The control structure for the present addressed transition is the same as the one found for the previous transition. The dynamic transition is shown in Figure 6. In this case, it is more clear how the compositions are controlled using one tight and one loose loop. Again, the tighter loop is the one controlling the composition of methanol in the distillate stream. The manipulated variables are also shown (parts c and d of Figure 6), and even though the optimal control trajectory shows
1942 Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009
an erratic behavior the closed-loop profile is smooth and accomplishes the transition in 20 min. 7. Conclusions In this work, a strategy for the simultaneous solution of design and control problems was proposed. A binary distillation column carrying out the separation of the methanol-water system was used to illustrate the main points in the application of the proposed solution strategy. The optimization formulation consisted of a MIDO problem that was transformed into a MINLP using the simultaneous dynamic optimization method (i.e., by discretizing the controlled and manipulated variables). Also, a decomposition algorithm for solving the MIDO problem based on solving relaxed versions of the optimization problem and using the results to initialize complex problem versions was described and successfully used. Using the proposed MIDO formulation, we were able to simultaneously design both the column and the controller that did the best job in tracking the also computed optimal openloop trajectory. The solution time required for the full space problem was very large. This could be possibly fixed to some extent by using some form of decomposition strategy.20,21 This point will be addressed in future work. Furthermore, in the present work it was very clear that the use of proper initialization strategies, such as the one used here, are crucial to solve large scale MINLPs. Appendix: High Index Reduction The mathematical model described in section 3 is a high index system of differential and algebraic equations (DAE); actually it is an index 2 DAE system. This can be noticed quickly by the fact that the tray vapor flow rates are algebraic variables that do not appear in algebraic equations. Because the DAE solver (DASKR) and the simultaneous dynamic optimization formulation require an index 1 system, then the index reduction scheme used by was applied here. The first step in the index reduction procedure is to write the left-hand side (LHS) of the energy balance as follows, taking into consideration that HkL ) f(xi,k,Tk):25
[
Mk
∂HLk dTk + ∂Tk dt
Nc
]
∂HLk dxi,k ) RHSEB i,k dt
∑ ∂x i
(38)
The partial derivatives of the liquid enthalpy are evaluated from eq 9, whereas the temperature derivative is obtained from the sum equations (eq 6) as follows. It is important to remember that in this case Ki,k ) f(xi,k, Tk). d dt
(∑ NC
NC
Ki,kxi,k -
i
∑x
i,k
i
) ∑[ NC
)
NC
xi,k
i
∂Ki,k dTk ∂Ki,k dxj,k + xi,k + dt ∂Tk dt j,k
∑ ∂x
]
j
dxi,k Ki,k dt
NC
∑ i
dxi.k d ) (0) ) 0 (39) dt dt
Because the thermodynamic constant is a function of tray composition, then its derivative with respect to time is calculated with the Jacobian matrix of partial derivatives of each constant with respect to the compositions. Finally, solving for (dTk)/(dt) dTk )dt
∑
NC i
[
xi,k
∑
NC j
dxi,k ∂Ki.k dxj,k + (Ki,k - 1) ∂xj,k dt dt ∂K NC i,k xi,k i ∂Tk
∑
]
(40)
Combining the RHS of the energy balance (eq 3) with the modified LHS in eq 38, and with the temperature differential in eq 40 combined with component mass balances (eq 2), a modified energy balance is obtained (eq 41). dxi,k NC NC ∂Ki,k dxj,k + (Ki,k - 1) xi,k j i ∂HLk ∂xj,k dt dt + Mk ∂Tk ∂Ki,k NC xi,k i ∂Tk
{
∑
[
]
∑
∑
NC
∑ i
}
∂HLk dxi,k ) RHSEB (41) ∂xi,k dt
where, L RHSEB ) Fk(HFk - HLk ) + Lk+1(Hk+1 - HLk ) +
V Vk-1(Hk-1 - HLk ) + Vk(HLk - HVk ) + Qk
This new equation is algebraic and it contains the vapor flow rate, thus reducing the index of the DAE system. Finally, the new, index reduced system is defined by eqs 1, 2, 4, 5, 6, 7, 8, and 41. Nomenclature Variables Tk ) temperature of stage k Mk ) total liquid molar holdup at tray k Fk ) feed stream at tray k Lk ) liquid flow rate at tray k Vk ) vapor flow rate at tray k CSF ) condenser split fraction xi,k ) molar fraction of component i at stage k zi,k ) molar fraction of component i at feed stage k Mkmin ) minimum liquid holdup for stage k Qk ) energy added or removed from stage k Ki,k ) thermodynamic equilibrium constant for component i at stage k γi,k ) activity coefficient of component i at stage k Pi,ksat ) vapor pressure of component i at stage k P ) operating pressure τji, k ) NRTL equation parameter of component i with respect to component j at stage k Gij, k, Gji, k ) NRTL equation parameter of one component with respect to the other HkL ) liquid enthalpy at stage k HkF ) liquid enthalpy of feed to tray k HkV ) vapor enthalpy at stage k Hi,kl ) liquid enthalpy of component i at stage k HkV ) vapor enthalpy of component i at stage k τk ) time constant of the Francis Weir equation L0 ) refluxed liquid stream D ) distillate stream NC ) number of components NT ) number of stages including reboiler and condenser Ncp ) number of collocation points Nfe ) number of finite elements Superindexes oc ) optimal control ss ) steady state cl ) closed-loop
Literature Cited (1) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Desing; Prentice Hall International Series in
Ind. Eng. Chem. Res., Vol. 48, No. 4, 2009 1943 the Physical and Chemical Engineering Series; Prentice Hall: United States of America, 1997. (2) Ziegler, J. G.; Nichols, N. B. Optimum Settings for Automatic Controllers. Transactions of ASME 1942, 64, 759–768. (3) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Integrated Flexibility and Controllability Analysis in Design of Chemical Processes. Comput. Chem. Eng. 2004, 43 (4), 997. (4) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. N. Recent Advandces in Optimization-Based Simultaneous Process and Control Design. Comput. Chem. Eng. 2004, 28, 2069. (5) Schweiger, C. A.; Floudas, C. A. Interaction of Design and Control: Optimization with Dynamic Models. Optimal Control: Theory Algorithms, and Applications 1998, 388. (6) Lo´pez-Negrete, R.; Lo´pez-Rubio, J.; Saldı´var-Guerra, E.; FloresTlacuahuac, A. Steady-State Multiplicity Behavior Analysis of a HighImpact Polystyrene Continuous Stirred Tank Reactor Using a Bifunctional Initiator. Ind. Eng. Chem. Res. 2006, 45 (5), 1689. (7) Lewin, D.; Bogle, D. Controllability Analysis of an Industrial Polymerization Reactor. Comput. Chem. Eng. 1996, 20, S871. (8) Ray, W. H.; Villa, C. Nonlinear Dynamics Found in Polymerization Processes s a review. Chem. Eng. Sci. 2000, 55, 275. (9) Ray, W. H. Polymerization Reactor Control 1986, 3–8. (10) Oldenburg, J.; Marquardt, W.; Heinz, D.; Leineweber, D. B. MixedLogic Dynamic Optimization Applied to Batch Distilation Process Design. AIChE J. 2003, 49 (11), 2900. (11) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N. A Case Study in Simultaneous Design and Control Using Rigurous, Mixed-Integer Dynamic Optimization Models. Ind. Eng. Chem. Res. 2002, 41, 760. (12) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal Design of Dynamic Systems Under Uncertainty. AIChE J. 1996, 42 (8), 2251. (13) Kameswaran, S.; Biegler, L. T. Simultaneous Dynamic Optimization Strategies: Recent Advances and Challenges. Comput. Chem. Eng. 2006, 30 (10-12), 1560. (14) Malcolm, A.; Polan, J.; Zhang, L.; Ogunnaike, B. A.; Linninger, A. A. Integrating Systems Design and Control using Dynamic Flexibility Analysis. AIChE J. 2007, 53 (8), 2048.
(15) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldı´var-Guerra, E. Dynamic Optimization of HIPS Open-Loop Unstable Polymerization Reactors. Ind. Eng. Chem. Res. 2005, 44 (8), 2659. (16) Lo´pez-Negrete, R.; Flores-Tlacuahuac, A. Optimal Operating Policies of a Reactive Distillation Column. Ind. Eng. Chem. Res. 2007, 46 (7), 2092. (17) Bonami, P.; Biegler, L. T.; Conn, A. R.; Cornue´jols, G.; Grossmann, I. E.; Laird, C. D.; Lee, J.; Lodi, A.; Margot, F.; Sawaya, N.; Wa¨chter, A. An Algorithmic Framework for Convex Mixed Integer Nonlinear Programs rep. IBM Research Report RC23771 2005. (18) Wa¨chter, A.; Biegler, L. T. On the Implementation of a PrimalDual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming. Mathematical Programming 2006, 106 (1), 25. (19) Aditya, K.; Daoutidis, P. Modeling, Analysis and Control of Ethylene Glycol Reactive Distillation Column. AIChE J. 1999, 45 (1), 51. (20) Fisher, M. An Application Oriented Guide to Lagrangian Relaxation. Interfaces 1985, 15 (2), 10. (21) Guignard, M.; Kim, S. Lagrangean Decomposition: A model yielding Stronger Lagrangean Bounds. Mathematical Programming 1987, 39, 215. (22) Viswanathan, J.; Grossmann, I. E. Optimal Feed Locations and Number of Trays for Distillation Columns with Multiple Feeds. Ind. Eng. Chem. Res. 1993, 32, 2942. (23) Biegler, L. T.; Cervantes, A. M.; Wa¨chter, A. Advances in Simultaneous Strategies for Dynamic Process Optimization. Chem. Eng. Sci. 2002, 57 (4), 575. (24) Cervantes, A. M.; Biegler, L. T. Large-Scale DAE Optimization Using a Simultaneous NLP Formulation. AIChE J. 1998, 44 (5), 1038. (25) Gear, C. W. Differential-Algebraic Equation Index Transformation. SIAM J. Sci. Comput. 1988, 9, 39.
ReceiVed for reView September 8, 2008 ReVised manuscript receiVed November 8, 2008 Accepted December 9, 2008 IE801353C
