Ind. Eng. Chem. Res. 2001, 40, 4079-4088
4079
PROCESS DESIGN AND CONTROL An Algorithm for Simultaneous Process Design and Control Ioannis K. Kookos and John D. Perkins* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College, London SW7 2BY, U.K.
The need to develop algorithms for the simultaneous solution of process design and control design problems has been identified by many researchers during the past two decades. However, the multifacetted nature of the problem and the computational burden of the resulting mathematical formulations have restricted progress in this field. This paper introduces a new algorithm for the solution of the simultaneous process and control design problem. The basic idea is to analyze a sequence of combined configurations making use of a bounding scheme to successively reduce the size of the search space. The bounding scheme proposed is based on a steady-state relaxation of the combined problem that generates valid lower bounds and a restriction problem that generates valid upper bounds. Application of the proposed algorithm is illustrated with a number of case studies. 1. Introduction Process design involves the selection of economically optimal process equipment and operating conditions. As in most synthesis problems, the design of chemical plants includes stochastic elements that stem from uncertainties in the operating environment and the knowledge of the process itself. The uncertainties that can affect the economic performance of a chemical plant during its lifetime have an extremely wide range and variety of sources. Changes in energy prices, legislation, environmental regulations, environmental conditions, feedstock quality, and equipment failure and wear, as well as imprecise knowledge of the physical and chemical phenomena involved, can dramatically affect the economic performance of a plant. Thus, the challenge of process design is to develop economically optimal systems that can be operated safely and reliably with acceptable environmental impact under a wide range of conditions and uncertain parameters. The need to incorporate elements of uncertainty into the design and operation of chemical plants was identified as early as the 1950s. The first works considered stochastic uncertainties and aimed to propose formal methodologies for the calculation of safety factors commonly used to overdesign process equipment to account for uncertainty.1-3 A number of interesting works appeared in the 1970s, including the works by Takamatsu et al.4 and Dittmar and Hartmann5 where the distinction between design and control variables is proposed for the first time. The methodology presented by Grossmann and Sargent6 and extended later by Grossmann and co-workers7 * To whom all correspondence should be addressed. Tel.: +44(0)207 5945556. Fax: +44(0)207 5945639. E-mail:
[email protected].
has dominated research in this field. This methodology identifies the need to guarantee the existence of a feasible region of operation for the specified range of uncertain parameters and, in addition, includes the first mathematical formulation of the distinction between static degrees of freedom (design variables) and dynamic degrees of freedom (control variables). The distinction between the dynamic and static degrees of freedom is of primary importance and has contributed significantly to a better understanding of the problem of process design under uncertainty. The main implication that stems from this distinction is that, during the operation of a plant, the process control system plays a key role in guaranteeing feasible and optimal operation of the plant. Thus, design decisions should be directed toward assisting the control system in achieving this task. It is not surprising that the process control community identified at the same time the need to consider control aspects during process design in order to improve controllability of the resulting chemical plants.8,9 Walsh and Perkins10-12 and Narraway and Perkins13 were among the first to consider general mathematical programming techniques for the simultaneous design and control problem using dynamic process models. Dimitriadis and Pistikopoulos14 applied the ideas of Halemane and Grossmann15 to systems described by sets of differential and algebraic equations. Mohideen et al.,16 on the basis of the aforementioned methodologies, proposed a unified process design framework for obtaining integrated process and control system designs under uncertainty. Bahri et al.17,18 and Schweiger and Floudas19 proposed alternative formulations and solution procedures. The main drawback of these methodologies stems from the fact that the combined design and control problem, a mixed integer dynamic optimization problem (MIDO), is solved as a mixed integer nonlinear programming problem (MINLP) by trans-
10.1021/ie000622t CCC: $20.00 © 2001 American Chemical Society Published on Web 08/24/2001
4080
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001
forming the system of differential and algebraic equations into algebraic equations by using either full discretization or integration. In the former case, the size of the resulting MINLP is explosive, and in the latter case, extensive dual information is needed to formulate a (mixed integer linear) master problem that is everincreasing in size. Recently, Bansal20 proposed a promising algorithm that alleviates most of these drawbacks. This paper presents a decomposition algorithm for solving the combined process and control design problem. The main idea is to systematically generate progressively tighter lower and upper bounds on the optimal economic performance of the plant. Lower bounds are generated by solving an improved steadystate flexibility problem (relaxation), where the topology of the plant and the controller are optimized, while the upper bounds are generated by solving a dynamic optimization problem (restriction), where the topology is fixed and the time-invariant design parameters are optimized. The structure of this paper is as follows. In section 2, the conceptual mathematical formulation of the combined process and control design problem is presented, and the proposed decomposition is outlined. In section 3, the proposed algorithm is presented, followed by case studies in section 4. Finally, this work is concluded in section 5.
of potential measured variables, f is the vector function of the differential equations of the process, h is the vector function of the algebraic equations of the process, g is the vector function of the inequality constraints that define the feasible operation, φ is the vector function of the differential equations of the controller, η is the vector function of the algebraic equations of the controller, and µ is the functional relationships of the measurements and the differential and algebraic equations. J is the objective function whose expected value (E) is to be minimized over the time period of interest. The above problem aims to choose the optimal process topology (X) and the optimal process design (d), as well as the controller structure (XC) and parameters (dC) that minimize the expected value of the objective function (E(J(.)) and at the same time satisfy the feasibility constraints over the given time horizon. Therefore, it is an infinite-dimensional, stochastic, mixed integer dynamic optimization problem (SMIDO). The requirement that feasibility of the solution be ensured for the whole set of the uncertain parameters under dynamic conditions, i.e.14
2. Conceptual Mathematical Formulation
complicates the solution of the process and control design problem. To impose this feasibility constraint, a multiperiod formulation is used in this study as proposed by Grossman and Sargent.6 In this approach, a finite set of uncertain parameters is selected (scenarios), corresponding, in most cases, to the vertexes of the space of the uncertain parameters, and a multiperiod problem is formulated where each realization of the uncertainty vector corresponds to a discrete period. Then, the solution to the design problem is obtained. In the final step, a constraint maximization problem is solved with fixed design, and the set of scenarios is updated to include the realization(s) of the uncertain vector that cause the most severe violation(s) of the constraints. The algorithm terminates when feasibility is ensured for all possible realizations of the uncertain parameters. First, we assume that the vector of uncertainty can be expressed as a function of a finite number of timeinvariant but uncertain parameters. It should be noted that θ(t) consists of disturbances (such as feed flow rate, feed composition, and temperature, etc.) and timevarying parameters (such as preexponential constants or activation energies in reaction rate expressions, heat transfer coefficients, etc.). A suitable description is the following13
The mathematical formulation of the process synthesis problem for continuous plants can be stated as follows16,21
min E[J(x(t),z(t),u(t),d,X,XC,θ(t))] d,X,dC,XC s.t. f(x˘ (t),x(t),z(t),u(t),d,X,θ(t)) ) 0 h(x(t),z(t),u(t),d,X,θ(t)) ) 0 g(x(t),z(t),u(t),d,X,θ(t)) e 0 φ(χ˘ (t),χ(t),ζ(t),y(t),u(t),dC,XC) ) 0 η(χ(t),ζ(t),y(t),u(t),dC,XC) ) 0 µ(x(t),z(t),y(t)) ) 0
(P1)
x(t) ∈ X, z(t) ∈ Z u(t) ∈ U, θ(t) ∈ Θ d ∈ D, dC ∈ DC XC ∈ {0,1}nXC, X ∈ {0,1}nX t ∈ [0,τf] where x(t) is the vector of differential variables, z(t) is the vector of algebraic variables, u(t) is the vector of control variables (variables whose values can be changed during operation), d is the vector of process design variables, X is the vector of integer variables associated with the decisions that define the topology of the process, θ(t) is the vector of uncertain parameters, dC is the vector of controller design variables, XC is the vector of integer variables associated with the decisions that define the controller structure, χ(t) is the vector of differential variables of the controller, ζ(t) is the vector of algebraic variables of the controller, y(t) is the vector
max min max {gk(x(t),z(t),u(t),θ(t))} e 0 θ(t) u(t) k t∈[0,τf]
θl(t) ) θN l + ∆θl sin(ωlt + λl2π) L ωl eωl eωU l 0 eλl e1
}
∀l
(1)
(2)
where θN l is the nominal value of the lth uncertain parameter, ∆θl is the maximum expected deviation from the nominal value, ωl is the frequency of variation of the uncertain parameter, and λl2π its phase. This is a particularly rich description as it can be used to describe fast and slow variations of the disturbances and has the desirable characteristic that it is a continuous function with continuous derivatives.
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001 4081
We further assume that the process economics are determined at nominal conditions and that the stochastic problem P1 can be written in the form
lem (MINLP). Because formulation P2SS is obtained from formulation P2 by restricting the time horizon, it follows that
min JDYN ) J(xN,zN,uN,d,X,XC,θN) d,X,dC,XC
J/SS e J/DYN
s.t.
where * denotes the minimum value. In other words, because the dynamics of the plant can only further restrict the feasible space (when compared to the steady state), the dynamic economics cannot be better than the steady-state economics. It is noteworthy that any feasible solution of problem P2 is also a feasible solution for problem P2SS. This means that any structure that is feasible under dynamic conditions (t ∈ [0,τf]) is also feasible at steady state (t ) 0) because the latter is a subset of the former. However, the converse is not true. That is, a structure that is feasible at steady state can be dynamically infeasible. Those structures (as well as parameters) that correspond to steady-state feasible but dynamically infeasible solutions have to be excluded from the feasible space of problem P2SS for the feasible space of P2 to be obtained. On the basis of these arguments, we conclude that the feasible space of problem P2 is a subset of the feasible space of problem P2SS. In fact, P2SS is a relaxation of P2 that is a sufficient condition for eq 3 to hold. The steady-state formulation (P2SS) includes equations describing the controllers. Thus, to be able to solve this problem, a simplified representation of the specifications imposed by the controller used at steady state has to be devised. If we assume that the controller used incorporates integral action, then its steady-state description is greatly simplified and, in addition, becomes independent of the specific form of the controller used (linear or nonlinear, centralized or decentralized, etc.). Because the steady-state description of any controller that incorporates integral action is equivalent to the steady-state description of the perfect controller, the following mixed integer linear formulation of the perfect controller can be used13
p
p
p
p
p
f(x˘ (t),x (t),z (t),u (t),d,X,θ (t)) ) 0 h(xp(t),zp(t),up(t),d,X,θp(t)) ) 0 g(xp(t),zp(t),up(t),d,X,θp(t)) e 0 φ(χ˘ p(t),χp(t),ζp(t),yp(t),up(t),dC,XC) ) 0 η(χp(t),ζp(t),yp(t),up(t),dC,XC) ) 0 µ(xp(t),zp(t),yp(t)) ) 0 p p θpl (t) - θN l - ∆θl sin(ωl t + λl 2π) ) 0, ∀ l p
p
p
x (t) ∈ X,z (t) ∈ Z,u (t) ∈ U
}
∀p
(P2)
θp(t) ∈ Θ, d ∈ D, dC ∈ DC p ωLl eωpl eωU l , 0 e λl e1, ∀ l, p
XC ∈ {0,1}nXC, X ∈ {0,1}nX t ∈ [0,τf] It should be noted that the objective function is based on the steady-state economic performance of the process, while the multiperiod constraints (indexed by p) ensure feasibility of operation for all periods corresponding to different realizations of the uncertain parameters. Formulation P2 is a mixed integer dynamic optimization problem that is difficult to solve. However, this problem is simplified significantly if consideration is restricted to the steady state. In this case, problem P2 can be written as
min JSS ) J(xN,zN,uN,d,X,XC,θN) d,X,dC,XC s.t. f(xp,zp,up,d,X,θp) ) 0 h(xp,zp,up,d,X,θp) ) 0 g(xp,zp,up,d,X,θp) e 0 φ(χp,ζp,yp,up,dC,XC) ) 0 η(χp,ζp,yp,up,dC,XC) ) 0 µ(xp,zp,yp) ) 0 p θpl - (θN l + ∆θl sin(λl 2π)) ) 0, ∀ l
}
p SP U yLj (1 - δj) + ySP j e yj e yj + (1 - δj)yj U yLj δj e ySP j e yj δj
e upi e uOPT + πiuU uLi πi + uOPT i i i L OPT U (1 - πi)ui e ui e ui (1 - πi) ∀p
(3)
(P2SS)
xp ∈ X, zp ∈ Z, up ∈ U θp ∈ Θ, d ∈ D, dC ∈ DC 0 e λpl e1, ∀ l, p XC ∈ {0,1}nXC, X ∈ {0,1}nX which is a mixed integer nonlinear programming prob-
}
}
∀ j, p (4)
∀ i, p
(5)
where δ and π are binary variables denoting the selection of potential measurement j (δj ) 1) or the selection of the potential manipulated variable i (πi ) 1). If we restrict the class of controllers to be square, then
∑j δj - ∑i πi ) 0
(6)
or, in the most general case where it is required that more controls are available than measurements (a necessary condition for the perfect control to be feasible)
∑j δj - ∑i πi e 0 Using eq 6, formulation P2SS becomes
(7)
4082
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001
min JSS ) J(xN,zN,uN,d,X,XC,θN) OPT ,u ,d,X,π,δ ySP j i s.t. p
p
p
f(x ,z ,u ,d,X,θp) ) 0 h(xp,zp,up,d,X,θp) ) 0 g(xp,zp,up,d,X,θp) e 0 µ(xp,zp,yp) ) 0 p θpl - (θN l + ∆θl sin(λl 2π)) ) 0, ∀ l
∑j yLj (1
-
δj) + ySP j yLj δj
e e
δj -
∑i
πi ) 0
ypj e ySP j + ySP e yU j j δj
(1 -
}
∀p
[
min E min J(x,z,u,d,X,θ) d θ∈Θ u
δj)yU j
uLi πi + uOPT e upi e uOPT + πiuU i i i U (1 - πi)uLi e uOPT e u (1 π ) i i i
}
(P3)
}
∀ j, p
πi, δj ∈ {0,1}nXC, X ∈ {0,1}nX The solution of problem P3 gives an optimum topology of the plant (XSS), an optimum value for the design vector (d), and an optimum topology for the control system (structure) as defined by the integer vectors δSS, πSS, and J/SS, which is a valid lower bound to the actual objective J/DYN. By using the topologies of the plant and the controller suggested by the solution of problem P3, we can obtain the actual value of the objective function by solving the following problem
s.t. f(x˘ (t),x (t),z (t),u (t),d,XSS,θp(t)) ) 0 h(xp(t),zp(t),up(t),d,XSS,θp(t)) ) 0 g(xp(t),zp(t),up(t),d,XSS,θp(t)) e 0 p p φ(χ˘ (t),χ (t),ζp(t),yp(t),up(t),dC,δSS,πSS) ) 0 η(χp(t),ζp(t),yp(t),up(t),dC,δSS,πSS) ) 0 p
N
θl (t) - (θl
p
p
p
µ(xp(t),zp(t),yp(t)) ) 0 + ∆θl sin(ωlpt + λlp2π)) ) 0, ∀l p
x (t) ∈ X, z (t) ∈ Z, t ∈ [0,τf]
}
h(x,z,u,d,X,θ) ) 0 max min max{gk(x,z,u,d,X,θ)} e 0, ∀ k (P5) θ∈Θ u k
0 eλlp e1, ∀ l, p
min JDYN ) J(xN(0),zN(0),uN(0),d,θN) d,dC
f(x,z,u,d,X,θ) ) 0 g(x,z,u,d,X,θ) e 0
∀ i, p
θp ∈ Θ, d ∈ D, dC ∈ DC
p
]
s.t.
xp ∈ X, zp ∈ Z, up ∈ U
p
more importantly, a valid upper bound on the optimum of problem P2 (as P4 is derived from P2 by fixing the binary variables). It is noteworthy that the formulation of the steadystate specifications imposed by the regulatory control system, given by eqs 4 and 5, can be used to improve the mathematical formulation of the classical flexibility problem. Grossmann et al.7 proposed the following twostage programming formulation for dealing with uncertainty in design
However, this formulation is optimistic because an ideal adaptation of all available manipulated variables to partially known or even unmeasurable disturbances is assumed. Grossmann et al.7 identified this shortcoming and stated that “including more detailed information on the control scheme would make the problem virtually unmanageable”. As shown in this work, substantial information regarding regulatory control can be incorporated into this formulation in the form of linear constraints. In addition, the solution of this feasibility problem provides invaluable information concerning the topology (structure) of the regulatory control system that is able to reject all major disturbances with minimum cost. Thus, formulation P3 can be seen as an improved flexibility test problem that can be used to obtain tighter lower bounds on the achievable economic performance of a chemical plant in the presence of uncertainty. 3. Algorithm for Simultaneous Design and Control
∀p
(P4)
up(t) ∈ U, d ∈ D, dC ∈ DC Formulation P4 is derived from formulation P2 by fixing all integer variables. Thus, in the dynamic optimization problem P4, only the continuous and time-invariant variables are optimized. It should be noted that the solution of problem P4 gives an optimum value for the design variables of the process and the controller and,
If the number of scenarios in the problem is fixed, a general algorithm for the simultaneous solution of the design and control problems based on formulations P3 and P4 is as follows: Step 0. Set I ) 1, choose a set of periods and values for the uncertain parameters, and set J/DYN ) +∞. Step 1. Set I ) I + 1. Step 2. Solve problem P3 with the constraint JISS e / JDYN. If no feasible solution can be found, then the algorithm terminates, and the optimum is given by J/DYN. If a feasible solution is found, proceed to step 3. Step 3. Solve problem P4. If the problem is infeasible, go to step 5. If a feasible solution is found, then proceed to step 4. Step 4. If JIDYN e J/DYN, then set J/DYN ) JIDYN. Proceed to step 5. Step 5. Add an integer cut to exclude the current structure, and go to step 1. The proposed algorithm terminates at step 2 when no steady-state solution can be found with economics better than the current best dynamic economics.
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001 4083
It is more commonly the case that the uncertainty is defined parametrically with ranges specified for the uncertain parameters (see for example eq 2). In that case, step 3 of the algorithm has to be modified to impose the feasibility restriction for the whole set of uncertainty space. First, problem P4 is solved, and the design parameters of the plant and the controller are fixed. Then, the following problem is solved
ψIk ) max gk(x(t),z(t),u(t),d,X,θ(t)) ω,λ s.t. f(x˘ (t),x(t),z(t),u(t),d,XSS,θ(t)) ) 0 h(x(t),z(t),u(t),d,XSS,θ(t)) ) 0 φ(χ˘ (t),χ(t),ζ(t),y(t),u(t),dC,δSS,πSS) ) 0
Figure 1. Evaporator system. Table 1. Major Variables of the Evaporation Process Model
η(χ(t),ζ(t),y(t),u(t),dC,δSS,πSS) ) 0 µ(x(t),z(t),y(t)) ) 0
F1 F2 F3 F4 F5 F100 F200 T1 T2 T3 T100 T200 T201 C1 C2 P2 P100 L2
- ∆θl sin(ωlt + λl2π) ) 0, ∀ l θl(t) - θNOM l ωLl e ωl e ωU l , 0 e λl e 1
(P6)
If maxk ψIk e 0, the system is feasible, and the algorithm does not change. If maxk ψIk > 0, however, then the closed-loop system is clearly infeasible for at least one realization of the uncertainty vector. In this case, the corresponding realization that causes the maximum violation is added to the set of scenarios, and the algorithm iterates at step 3 until maxk ψIk e 0 is satisfied. Because, in all case studies presented in section 4 of this paper, the disturbances are defined parametrically, the enhanced version of the algorithm will be applied. It should be noted that the proposed algorithm is a general methodology for solving SMIDO problems with static objective functions (such as those that include only installation costs). As a result, a number of process design problems (e.g., design of batch plants) with objective functions of dynamic nature cannot be handled using this framework. It is interesting, before closing this section, to discuss the class of controllers that can be incorporated in the proposed formulation. It should be noted that the formulation is independent of the type of controller used provided that the controller can be represented using a finite number of time-invariant parameters. This general class encompasses all linear time-invariant controllers, as well as a large number of nonlinear controllers. All of the constraints included in, for example, a MPC formulation are incorporated in the problem definition, and as a result, the controller used is a constrained controller. A MPC controller cannot be incorporated in the proposed formulation because of its time-varying nature (the same holds true for any time-varying controller). 4. Case Studies Formulation P4 requires the form of the controller to be specified. In the case studies that follow, the controller is assumed to be multivariable. Each element of this
feed flow rate (kg/min) product flow rate (kg/min) circulation flow rate (kg/min) vapor flow rate (kg/min) condensate flow rate (kg/min) steam flow rate (kg/min) cooling water flow rate (kg/min) feed temperature (°C) product temperature (°C) vapor temperature (°C) steam temperature (°C) cooling water inlet temperature (°C) cooling water outlet temperature (°C) feed composition (%) product composition (%) operating pressure (kPa) steam pressure (kPa) separator level (m)
}
multivariable controller is a proportional-integral controller having the following mathematical description21
ui(t) ) uiOPT + ∆uij(t) )
{
Kij ζj(t) +
0,
d
(
∑j ∆uij(t) 1 τij
)
χj(t) , (i,j) ∈ IA × JA (i,j) ∉ IA × JA
(uiOPT) ) 0
dt
ζj(t) ) ySP j - yj(t) d (χ (t)) ) ζj(t) dt j
}
∀i
(8)
∀j
ζj(0) ) 0, ∀ j ∈ JA, χj(0) ) 0, ∀ j where IA){i: πi ) 1} and JA){j: δj ) 1}. The dif)/dt ) 0 is introduced to allow ferential equation d(uOPT i consistent initialization of the resulting system of differential algebraic equations.20 In all case studies, the MINLP problems were solved using GAMS/DICOPT22 based on the combined penalty function and outer approximation method developed by Viswanathan and Grossmann.23 Dynamic optimization problems were solved using gPROMS/gOPT.24 4.1. Control Structure Selection of an Evaporator Process. The evaporation process examined in this
4084
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001
case study is shown in Figure 1.25,26 This is a process that removes a volatile liquid from a nonvolatile solute, thus concentrating the solution. It consists of a heatexchange vessel with a recirculating pump. The overhead vapor is condensed by the use of a process heat exchanger. The major variables of interest are summarized in Table 1. The dynamic model of the evaporator is given by (a slight variation of the equations proposed by Newell and Lee25)
Differential equations M
dC2 ) F1C1 - F2C2 dt dP2 ) F4 - F5 C dt
(9) (10)
Algebraic equations F1 - F4 - F2 ) 0
(11)
F1CpT1 - F4(λ + CpT4) - F2CpT2 + Q100 ) 0
(12)
0.5616P2 + 0.3126C2 + 48.43 - T2 ) 0
(13)
0.5070P2 + 55 - T4 ) 0
(14)
0.1538P100 + 90 - T100 ) 0
(15)
Q100 - UA1(T100 - T2) ) 0
(16)
Q100 - F100λS ) 0
(17)
Q200 - F200Cp(T201 - T200) ) 0
(18)
Q200 - UA2(T4 - (T201 + T200)/2) ) 0
(19)
F5λ - Q200 ) 0
(20)
Inequality constraints g1 ) 25 - C2 e 0
(21)
g2 ) 40 - P2 e 0
(22)
g3 ) P2 - 80 e 0
(23)
g4 ) P100 - 400 e 0
(24)
g5 ) F200 - 600 e 0
(25)
control system for the evaporator. The available measurements and manipulated variables are the following
[] [ ]
[ ] [ ]
y C u P y ) y1 ) P 2 , u ) u1 ) F100 2 2 2 200
(29)
where C2 is the product composition, P2 is the operating pressure, P100 is the steam pressure, and F200 the cooling water flow rate. The steady-state optimal operating point is found first. Both disturbances are assumed to obtain their nominal values. At the optimal steady-state operating point, the input and output variables have the values
[ ] [ ]
[ ] [
C P 25.00 193.45 , u ) F100 ) y ) P2 ) 50.57 207.32 2 200
]
(30)
and the value of the objective function is 5651.16 × 10-3 $/h. Then, the steady-state multiperiod problem is solved. Five discrete periods are considered. The first period corresponds to the case where the disturbances have their nominal values. In the remaining four periods, the disturbances have values that correspond to the four vertexes of the uncertainty space
θ1 )
[ ]
[ ]
[ ]
[ ]
11 11 9 9 , θ2 ) , θ3 ) , θ4 ) 5.5 4.5 4.5 5.5
(31)
where θ ) [F1 C1]T. The solution of the multiperiod problem gives the u1,u2-y1,y2 control structure as the best structure and a slightly increased objective function (JSS ) 5651.17 × 10-3 $/h). Next, the multiperiod dynamic optimization problem is solved, where the parameters of the multivariable PI controller, corresponding to the control structure u1,u2y1,y2, are to be determined. An initial set of two periods (plus the nominal period) is chosen to obtain an initial set of controller parameters. In the first period, the disturbances have the values
(π2)
F11(t) ) 10 + sin
(32)
[ () ] 3π t + τf 2
C11(t) ) 5 + 0.5 sin 2π
(33)
while in the second period, they have the values
-3
Objective function (operating cost in 10 $/h) J ) 0.6F200 + 1.009F2 + 600F100
where T1 ) 40, T200 ) 25, M ) 20, C ) 4, Cp ) 0.07, λ ) 38.5, λS ) 36.6, UA1 ) 9.6, and UA2 ) 6.84. The feed flow rate and composition are uncertain parameters having the following description
F1 ) 10 + sin(2πλF1), 0 eλF1 e1
[
()
C1(t) ) 5 + 0.5 sin 2πωC1t
(π2)
F21(t) ) 10 + sin
(26)
(27)
]
t + 2πλC1 , τf
100 e ωC1 e 102, 0 e λC1 e 1 (28) and τf ) 1440 min was chosen. It should be noted that a high-frequency (noise-like) variation is assumed for the feed composition disturbance, and a low-frequency (step-like) variation is assumed for the feed flow rate. The aim of this case study is to design a regulatory
[
C21(t) ) 5 + 0.5 sin 2π100
(34)
() ] 3π t + τf 2
(35)
The solution corresponds to a dynamic economics of JDYN ) 5652.55 × 10-3 $/h. The controller parameters are
K)
[
131.85 57.08
]
[
0.00 1.0 1.0 , τ) -189.08 9.9 9.6
]
(36)
while the setpoints are
C2SP ) 25.0243 P2SP ) 50.6
(37)
Then, the constraint maximization problem is solved
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001 4085 Table 2. Summary of the Solution of the Evaporator Case Study solution of problem P3 iteration 0 1 2
manipulated controlled variables variables u 1, u 2 u1
y 1, y 2 y1
solution of problem P4 JSS
JDYN
/ JDYN
5.65116 5.65318
5.65255
+∞ 5.65255
with fixed controller structure and parameters, where each constraint gk, k ) 1, 2, ..., 5, is maximized with respect to λF1, ωC1, and λC1. All ψk satisfy ψk < 0, and as a result, the corresponding closed-loop system is feasible for the whole range of uncertainty assumed. At this point, iteration 1 terminates. An integer cut is added to exclude structure u1,u2-y1,y2 from further consideration, and iteration 2 starts. The steady-state multiperiod problem is solved again for the same set of periods. The best control structure obtained is u1-y1 (P100-C2), and the value of the objective function is 5653.78 × 10-3 $/h. Because the steady-state objective is worse than the current best dynamic economics (JDYN ) 5652.55 × 10-3 $/h) obtained in iteration 1, the algorithm terminates. The optimum solution corresponds to the 2-by-2 structure P100,F200-C2,P2. Table 2 summarizes the results of this case study. 4.2. Simultaneous Design and Control of a Binary Distillation Column. The distillation column considered in this case study is shown in Figure 2. Details of the mathematical model can be found in Luyben and Floudas27 and Schweiger and Floudas.19 The aim of this case study is to design the distillation column (number of trays, feed tray location, and column diameter), as well as the structure and parameters of the control system, to account for 20% uncertainty in the feed flow rate and 10% uncertainty in feed composition. The constraints selected in this case study are the following
0.98 - xD e 0
(38)
xB - 0.02 e 0
(39)
0.6719xV - DC e 0
(40)
where xD is the top product composition, xB is the bottom product composition, V is the vapor flow rate, and DC the column diameter. The potential measured and manipulated variables are
[ ]
[]
x R y ) xD , u ) V B
(41)
The liquid level in the reflux drum (reboiler) is controlled using the distillate (bottom product) flow rate. The description of the uncertainty is as follows
F ) 1 + 0.2 sin(2πλF), 0 e λF e 1
[ ()
z(t) ) 0.45 + 0.05 sin 2πωz
(42)
]
t + 2πλz , τf
10-2 e ωz e 102, 0 e λz e 1 (43) where F is the feed flow rate and z(t) the feed composition.
Figure 2. Binary distillation column.
Figure 3. Reactor-separator with recycle process.
The results and the steps of the proposed algorithm are summarized in Table 3. First, the nominal case (no disturbances) is solved, and the minimum cost obtained is 34 182 $/year. Using this as the initial point, formulation P3 is solved, and the structure obtained corresponds to a column that is 0.858 m in diameter and has 17 trays, with tray number 9 as the feed tray. The 2-by-2 control structure R,V-xD,xB is obtained at this step. Then, formulation P4 is solved with fixed topology, and the column diameters and controller parameters are optimized. Two periods are considered corresponding to
[][ ][][ ] λ1F 0.25 ω1z ) 10-2 , 0.25 λ1z
λ2F 0.25 ω2z ) 102 0.25 λ2z
(44)
The controller parameters are first obtained, and the dynamic feasibility problem is then solved. A third period, where
4086
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001
Table 3. Summary of the Solution of the Binary Distillation Column Case Study iteration nominal case 0.7535 16 8 34182
1
2
column diameter number of trays feed tray location total cost control structure
0.8583 17 9 36813 V, R - xB, xD
formulation P3 0.8404 18 9 36881 V, R - xB, xD
column diameter xB set point xD set point total cost
0.8625 0.9801 0.0194 36972
formulation P4 0.8456 0.9801 0.0194 37082
3
4
5
0.8610 17 8 36921 V, R - xB, xD
0.8810 16 8 36929 V, R - xB, xD
0.8867 16 9 37305 V, R - xB, xD
0.8669 0.9801 0.0194 37150
0.8856 0.9801 0.0194 37104
Table 4. Summary of the Solution of the Reactor-Separator with Recycle Case Study iteration nominal case
1
1.015 23 23 4.310 2009 487 774
2
3
column diameter number of trays feed tray location reactor diameter reactor volume total cost control structure
formulation P3 1.168 20 20 4.374 2105 494 271 V, F - xB, VR
1.166 21 21 4.361 2087 494 327 V, F - xB, VR
1.171 19 19 4.391 2130 494 570 V, F - xB, VR
column diameter reactor diameter reactor volume total cost
formulation P4 1.168 4.374 2107 494 453
1.166 4.363 2090 494 533
[][ ] λ3F 0.25 ω3z ) 20.72 0.25 λ3z
(45)
(that results in violation of the constraints) is added, and the optimal controller parameters and column diameter are calculated again. Then, the dynamic feasibility problem is addressed. The new closed-loop system is found to be feasible for the whole set of uncertainty considered. The design variables that correspond to the final design are given in Table 3, and the controller parameters are
K)
[
50.0 -18.0
]
[
-17.8 10.0 , τ) -50.0 143.2
180.0 1.0
]
(46)
It should be noted that the controller gains are restricted to lie in the interval [-50, 50], while the reset times must lie in the interval [1, 180]. This design corresponds to a total cost of 36 972 $/year. This is a valid upper bound for the objective function in the subsequent iterations. As can be seen from Table 3, the algorithm converges in the fifth iteration, because the solution obtained for the steady-state multiperiod problem has steady-state economics (37 305 $/year) that is worse than the current best dynamic economics (36 972 $/year). It is interesting to note that, in all intermediate solutions, the control structure chosen is the same and corresponds to the twopoint composition control structure. Finally, it should be noted that the final economic penalty is 2790 $/year (8.2% of the nominal minimum cost). 4.3. Simultaneous Design and Control of a Reactor-Separator-Recycle System. In this case study, the simultaneous process and control design of the reactor-separator-recycle system shown in Figure 3
is considered. Details of the mathematical model can be found in Luyben and Floudas28 and Schweiger and Floudas.19 Fresh feed, 90% in A, is fed into the reactor, where the first-order irreversible reaction A f B takes place. Reactor product is fed into the distillation column. Unreacted A is recycled to the reactor, while product B (containing less than 0.0105 A) is obtained as the bottom product. The aim of this case study is to design the reactor and the distillation column (reactor diameter, number of trays, feed tray location, and column diameter), as well as the structure and parameters of the control system, to account for 20% uncertainty in the feed flow rate and 5% uncertainty in the feed composition
F0 ) 1.811 67 + 0.3623 sin(2πλF0),0 e λF0 e 1 (47)
[ ()
z0(t) ) 0.90 + 0.05 sin 2πωz0 10
-2
]
t + λz02π , τf 1
e ωz0 e 10 , 0 e λz0 e 1 (48)
where F0 is the feed flow rate and z0(t) is the mole fraction of A in fresh feed. The constraints considered are related to the maximum composition of A in the product (0.0105), the column flooding (eq 40), and an upper bound on the reactor level. The potential measurements selected are (1) the bottoms product composition (xB), (2) the top product composition (xD), (3) the reactor composition (xR), and (4) the reactor level (VR). The potential manipulated variables are the steam flow rate in the reboiler (V), the recycle flow rate (FR), and the reactor product flow rate (F). The liquid level in the bottom of the column is assumed to be controlled using the bottoms product flow rate, and the liquid level in the reflux drum is controlled using the reflux flow rate. The optimum process design
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001 4087
algorithm converges on the third iteration when the steady-state economics (494 570 $/year) becomes greater than the best dynamic economics (494 453 $/year) obtained in iteration 1. The best solution is also shown in Figure 4. Finally, it should be noted that the economic penalty to compensate for the effect of disturbances is 6679 $/year (1.4% of the nominal minimum cost). 5. Conclusions In this paper, a new algorithm for the simultaneous design of process and control systems is presented. This algorithm is based on the systematic generation of lower and upper bounds on the best achievable dynamic economics of the combined plant to reduce effectively the size of the search space. The lower bounds are generated by solving a relaxed problem that is based on the best steady-state economic performance of the plant. The upper bounds are generated by solving a restricted problem in which the structure of the process and the controller (topology) are fixed and the parameters are optimized. The proposed algorithm can effectively be used to answer important questions arising in simultaneous design and control, as was demonstrated using a number of case studies. Nomenclature Figure 4. Optimum solution for the reactor-separator with recycle process.
solution (shown in Table 4) is obtained first by assuming that the disturbances have their nominal values. Using the optimum solution of the nominal case as the initial point, formulation P3 is solved. The solution obtained is shown in Table 4. It is interesting to note that the optimum structure obtained at this step corresponds to a stripper column that is 1.168 m in diameter and has 20 trays and a reactor that is 4.374 m in diameter and has a 2105-kmol holdup. The composition of the bottoms product and the reactor volume are selected as the controlled variables, and the reactor effluent and steam flow rate are selected as the manipulated variables (structure V,F-xB,VR). The total cost is 494 271 $/year. Using the structural information obtained from the solution of the steady-state multiperiod problem P3, the dynamic optimization problem P4 is solved, and the cost obtained and the column and reactor diameters are given in Table 4. The controller parameters and the uncertainty scenarios that correspond to the final system are
K)
[
]
[
-10.0 0 1.0 360 , τ) -10.0 -9.8 360 1.0
]
(49)
[ ][ ][ ][ ][ ][ ] λF01 0.75 1 ωz0 ) 100 , 0.75 λz01
λF02 0.75 2 ωz0 ) 102 , 0.75 λz02
λF03 0.25 ωz03 ) 4.72 0.75 λz03 (50)
The total cost is 494 453 $/year. This cost is an upper limit for the objective function of the steady-state multiperiod problem in the subsequent iterations. The solutions of the steady-state multiperiod problem and the dynamic optimization problems obtained in the subsequent iterations are given in Table 4. The reactor volume and bottoms composition are selected as controlled variables in all intermediate solutions. The
d ) vector of process design variables dC ) vector of controller design variables E ) expected value f ) vector function of the differential equations of the process g ) vector function of the inequality constraints that define the feasible operation h ) vector function of the algebraic equations of the process J ) cost function p ) index of periods t ) time u ) vector of control variables (variables that can change during operation) x ) vector of differential variables X ) vector of integer variables associated with the decisions that define the topology of the process XC ) vector of integer variables associated with the decisions that define the structure of the process control system y ) vector of potential measured variables z ) vector of algebraic variables Greek Letters ∆ ) deviation from the nominal value δ ) integer vector defined in eq 4 ζ ) vector of algebraic variables of the controller η ) vector function of the algebraic equations of the controller θ ) vector of uncertain parameters λ ) parameter defined in eq 2 µ ) functional relationships of the measurements and the differential and algebraic equations π ) integer vector defined in eq 5 φ ) vector function of the differential equations of the controller χ ) vector of differential variables of the controller ω ) frequency Subscripts SS ) steady state DYN ) dynamic
4088
Ind. Eng. Chem. Res., Vol. 40, No. 19, 2001
Superscripts OPT ) optimum value SP ) set point
Literature Cited (1) Buckley, P. S. Sizing Process Equipment by Statistical Methods. Chem. Eng. 1950, 112. (2) Saletan, D. I.; Caselli, A. V. Optimum design capacity of new plants. Chem. Eng. Prog. 1963, 59, 69. (3) Kittrell, J. R.; Watson, C. C. Don’t overdesign process equipment. Chem. Eng. Prog. 1966, 62, 79. (4) Takamatsu, T.; Hashimoto, I.; Shioya, S. 1973. On design margin for process system with parameter uncertainty. J. Chem. Eng. Jpn. 1973, 6, 453. (5) Dittmar, R.; Hartmann, K. Calculation of optimal design margins for compensation of parameter uncertainty. Chem. Eng. Sci. 1976, 31, 563. (6) Grossmann, I. E.; Sargent, R. W. Optimum Design of Chemical Plants with Uncertain Parameters. AIChE J. 1978, 24, 1021. (7) Grossmann, I. E.; Halemane, K. P.; Swaney, R. E. Optimization strategies for flexible chemical processes. Comput. Chem. Eng. 1983, 7, 439. (8) Morari, M. Flexibility and resiliency of process systems. Comput. Chem. Eng. 1983, 7, 423. (9) Perkins, J. D.; Wong, M. P. F. Assessing controllability of chemical plants. Chem. Eng. Res. Des. 1985, 63, 358. (10) Walsh, S.; Perkins, J. D. Integrated design of effluent treatment systems. In Interactions between Process Design and Process Control; Perkins, J. D., Ed.; Pergamon Press: London, U.K., 1992. (11) Walsh, S.; Perkins, J. D. Application of integrated process and control system design to waste water neutralisation. Comput. Chem. Eng. 1994, 18, S183. (12) Walsh, S.; Perkins, J. D. Operability and Control in Process Synthesis and Design. In Advances in Chemical Engineering; Academic Press: New York, 1996; Vol. 23, p 172. (13) Narraway, L. T.; Perkins, J. D. Selection of Process Control Structure Based on Linear Dynamic Economics. Ind. Eng. Chem. Res. 1993, 32, 2681. (14) Dimitriadis, V. D.; Pistikopoulos, E. N. Flexibility Analysis of Dynamic Systems. Ind. Eng. Chem. Res. 1995, 34, 4451. (15) Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1983, 29, 425.
(16) Mohideen, J. M.; Perkins, J. D.; Pistikopoulos, E. N. Optimal Design of Dynamic Systems under Uncertainty. AIChE J. 1996, 42, 2251. (17) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Effect of Disturbances in Optimizing Control: Steady-State Open Loop Back-Off Problem AIChE J. 1996, 42, 983. (18) Bahri, P.; Bandoni, J.; Romagnoli, J. Integrated Flexibility and Controllability Analysis in Design of Chemical Processes. AIChE J. 1997, 43, 997. (19) Schweiger, C. A.; Floudas, C. A. Interaction of Design and Control: Optimisation with Dynamic Models. In Optimal Control: Theory, Algorithms and Applications; Hager, W. W., Pardalos, P. M., Eds.; Kluwer Academic Publishers: Norwell, MA, 1997. (20) Bansal, V. Analysis, Design and Control Optimization of Process Systems under Uncertainty. Ph.D. Dissertation, University of London, London, U.K., 2000. (21) Narraway, L. T.; Perkins, J. D. Selection of process control structure based on economics. Comput. Chem. Eng. 1994, 18, S511. (22) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS Release 2.25: A User’s Guide; The Scientific Press: San Francisco, CA, 1992. (23) Viswanathan, J.; Grossmann, I. E. A Combined penalty function and outer-approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14, 769. (24) gPROMS Advanced User’s Guide; Process Systems Enterprise Ltd.: London, 1999. (25) Newell, R. B.; Lee, P. L. Applied Process ControlsA Case Study; Prentice Hall: New York, 1989. (26) Wang, F. Y.; Cameron, I. T. Control studies on a model evaporation processsConstrained state driving with conventional and higher relative degree systems. J. Process Control 1994, 4, 59. (27) Luyben, M. L.; Floudas C. A. Analyzing the Interaction of Design and Controls1. A Multiobjective Framework and Application to Binary Distillation Synthesis. Comput. Chem. Eng. 1994, 18, 933. (28) Luyben, M. L.; Floudas C. A. Analyzing the Interaction of Design and Controls2. Reactor-Separator-Recycle System. Comput. Chem. Eng. 1994, 18, 971.
Received for review June 29, 2000 Revised manuscript received July 3, 2001 Accepted July 3, 2001 IE000622T
