Parametric Controllers in Simultaneous Process and Control Design

Ind. Eng. Chem. Res. 2003, 42, 4545-4563

4545

Parametric Controllers in Simultaneous Process and Control Design Optimization Vassilis Sakizlis, John D. Perkins, and Efstratios N. Pistikopoulos* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, U.K.

While design optimization has been employed effectively in the past few years to address the integration of design and control, the implementation of advanced model-based control techniques in such a framework has largely been ignored despite their advantages. The reason is that the advanced control schemes involve solving an optimization problem on-line, leading to a drastic increase in the complexity of the design framework. In this work a novel method is presented for integrating advanced model-based predictive controllers in a simultaneous process and control design methodology. Parametric programming is used for the controller derivation, giving rise to an explicit controller structure and removing the need for solving an optimization problem online. The resulting parametric controller is readily incorporated in the design optimization framework, bringing about significant economic and operability benefits. The key features and advantages of our approach are highlighted via binary distillation and evaporator motivating examples. 1. Introduction During the past decades, there has been growing awareness both in academia and in industry that operability issues need to be considered explicitly at the early stages of process design. As a result, a number of methodologies have been developed for addressing the interactions between process design and process control (for a review see, work by Van Schijndel and Pistikopoulos1). A large proportion of the work in this field has concentrated on the application of controllability metrics based, for instance, on the open- and closed-loop stability analysis,2,3 open-loop output/input achievable performance (output controllability index),4 or linear quadratic Guassian criteria.5,6 These measures provide some assessment of a system’s operability but may not relate directly and unambiguously to real performance requirements, while they do not address directly the selection of the best process and control design. The most systematic approaches in that field determine simultaneously the optimal process and control design using advanced optimization algorithms that accommodate both continuous and discrete decisions about the plant design, control, and operation (for a review see, work by Bansal et al.7). However, most of these approaches utilize conventional proportional-integralderivative controllers that usually feature a single-input single-output structure and do not explicitly handle process specifications, i.e., the environmental, safety, and operational constraints. The shortcomings of conventional control schemes can be overcome by pursuing an advanced model-based predictive control (MPC) scheme. MPC exhibits performance benefits for a wide class of complex multivariable processes because it handles efficiently the presence of operating constraints. The work of other groups on the incorporation of online optimization-based controllers in process design8-10 and control structure selection11,12 is indicative of the potential of MPC schemes to improve the process * To whom correspondence should be addressed. Tel.: (44) (0) 20 7594 6620. Fax: (44) (0) 20 7594 6606. E-mail: [email protected].

economics and operability. However, most of those approaches solve the resulting multilevel process design optimization problems by making significant simplifications in the MPC controller structure, thereby removing some of its advantageous features. The development of a new type of advanced controllers, the so-called modelbased parametric controllers, can potentially address these issues. This control methodology is based on newly proposed parametric programming algorithms, developed at Imperial College,13 which derive the explicit mapping of the optimal control actions in the space of the state measurements. Thus, a simple explicit state feedback controller is derived that moves off-line the embedded online control optimization and preserves all of the beneficial features of MPC. The aim of this work is to incorporate model-based parametric controllers into a simultaneous process and control design framework. This is performed in an iterative manner. In the first step, a controller is derived based on a linearized version of the process model for fixed tunings and discrete decision variables. Then the explicit form of the controller is included in the rigorous dynamic model of the process, and the resulting closedloop system is optimized in terms of economic criteria. A control master follows, providing new controller tuning variables. When no tunings can be found that improve the system performance, a structural master problem is solved, providing new values for the discrete design decision variables. The rest of the paper is organized as follows: In the next section, an overview of the newly developed methods for (i) process and control design under uncertainty, (ii) solution of mixed-integer dynamic optimization (MIDO) problems, and (iii) derivation of the explicit model-based optimizing control law is presented, with the perspective of incorporating them in the proposed simultaneous process and control design framework of this paper. In the following paragraph, the problem formulation is stated. In section 4, the simultaneous process and control design algorithm featuring parametric controllers is described via a motivating example. Results for the examples are presented in sections 5 and

10.1021/ie0209273 CCC: $25.00 © 2003 American Chemical Society Published on Web 05/23/2003

4546 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003

6, demonstrating the significant economic and operability benefits of our approach. In the last section, conclusions are discussed. 2. Preliminaries 2.1. Framework for Simultaneous Design and Control under Uncertainty. Here, the simultaneous process and control design framework of Mohideen et al.14 and Bansal et al.7 is reviewed as one of the most suitable approaches in this field. The problem of the integration of process design, process control, and process operability can be conceptually posed as follows:

minimize

expected total annualized cost

subject to

differential-algebraic process model, inequality path constraints control scheme equations, process design equations feasibility of operation (over time) process variability constraints

P

A general, algorithmic framework for solving (P) is schematically shown in Figure 1 and can be summarized as follows: Step 1. Choose an initial set of scenarios for the uncertain parameters. Step 2. For the current set of scenarios, determine the optimal process and control design by solving the (multiperiod) MIDO problem: ns

min d,δ,vv1,vv2,...,vvns i

wiφ[x˘ i(tf),xi(tf),xai(tf),vi(tf),viv,θi(tf),d,δ,tf] ∑ i)1

(1)

i

s.t. fd[x˘ (t),x

(t),xia(t),vi(t),viv,θi(t),d,δ]

)0

fa[xi(t),xia(t),vi(t),viv,θi(t),d,δ] ) 0 g[x˘ i(t),xi(t),xia(t),vi(t),viv,θi(t),d,δ] e 0 f0[x˘ i(t0),xi(t0),xia(t0),vi(t0),viv,θi(t0),d,δ] ) 0 t 0 e t e tf

i ) 1, ..., ns

where x denotes the differential state variables of the system pertaining, e.g., to molar holdup and internal energies, xa includes the algebraic variables, e.g., geometric characteristics, thermodynamic variables, and internal flow rates, δ ∈ Y ≡ {0, 1}Nδ comprises the binary variables for the process and the control structure (corresponding to, e.g., the number of trays in a distillation column or whether a manipulated variable is paired with a particular controlled variable or not), v denotes the vector of time-varying manipulated variables, such as the utility flow rates, vv is the set of timeinvariant operating variables (e.g., set point of controllers, utility flows), d are the design variables that cannot be readjusted during operation (e.g., equipment size), i is the index set for the scenarios of the uncertain parameters θ(t) that can be time-varying or timeinvariant, ns is the number of scenarios, wi, i ) 1, ..., ns, are discrete probabilities for the selected scenarios ns (∑i)1 wi ) 1), φ is usually an economic index of performance that may include weights on the operability or the environmental impact on the plant, fd represents the differential equations corresponding, e.g., to mass and energy balances, fa denotes the algebraic equations pertaining, e.g., to thermodynamic and hydraulic rela-

Figure 1. Decomposition algorithm of Mohideen et al.14 and Bansal et al.7

tions, f0 are the initial conditions of the dynamic system, e.g., the system is initially at steady state x˘ ) 0, g e 0 represents the set of constraints (end, point, and path) that must be satisfied for feasible operation (e.g., purity specifications, environmental restrictions). Step 3. Test the process and control design from step 2 for feasibility over the whole range of uncertain parameters by solving the dynamic feasibility test problem:

χ ) max min vv

θ

max gl(‚)

l∈q, t∈[0,tf]

(2)

s.t. fa(‚), fd(‚), f0(‚) ) 0 If χ(d,δ) e 0, feasible operation can be ensured dynamically for all values of θ within the given ranges. In this case, the algorithm terminates; otherwise, the solution of (2) defines a critical scenario that is added to the current set of scenarios before returning to step 2. Remarks. 1. The formulation (P) is an exact closedloop, dynamic analogue of the steady-state problem of optimal design with a fixed degree of flexibility.15 The solution strategy shown in Figure 1 and described above is a closed-loop dynamic analogue of the flexible design algorithm of Grossmann et al.16 2. The integrated design and control problem requires the solution of MIDO problems in steps 2 and 3. Until recently, there were no reliable methods for dealing with such problems. In the next paragraph, a newly developed MIDO algorithm17,18 is briefly outlined. 2.2. MIDO. Consider a general MIDO formulation similar to (1) (single period for brevity):

min φ[x˘ (tf),x(tf),xa(tf),vv,d,δ,tf] vv,δ,d

s.t. 0 ) fd[x(t),x(t),xa(t),vv,d,δ,t] 0 ) fa[x(t),xa(t),vv,d,δ,t] 0 ) f0[x(t0),x˘ (t0),xa(t0),vv,d,δ,t0] 0 g g[x˘ (tf),x(tf),xa(tf),vv,d,δ,tf] t 0 e t e tf

(3)

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4547

The binary variables δ participate only in a linear form in the objective function, the differential system, and the constraints. Here, the path constraints (e.g., minimum allowable temperature difference in the heat exchangers over the complete horizon) are equivalently transformed to end-point constraints,19 as is clearly shown through the example in section 5. The MIDO approach reviewed here is based on generalized Benders decomposition (GBD), where the original MIDO problem is decomposed to a primal and a master subproblem. The primal problem is constructed by fixing the binaries to a specific value δ ) δk. Then problem (3) becomes an optimal control problem. In GBD-based approaches, the master problem is constructed using the dual information of the primal problem at the optimum solution. The dual information is embedded in the Lagrange multipliers µ of the constraints g and the adjoint time-dependent variables λ(t) and p(t) that are associated with the differential system of equations, i.e., fd and fa. The evaluation of the adjoint variables requires an extra integration of the so-called adjoint differential algebraic equation (DAE) system. After the adjoint functions are calculated, the master problem is constructed and has the following form:

min η δ,η

s.t.

[]

[] ∫[ ] [ ]

f η g φ + (µκ)T‚g + (ωfκ)T‚ fd a (Fκ)T‚f0 + κ∈K

f + (ω0κ)T‚ fd a f

0

T f λκ(t) ‚ fd κ a p (t)

t

tf

t0

+ dt (4)

where F, ωf, and ω0 are multipliers that are evaluated from the first-order optimality conditions of the optimal control primal problem.20 The demanding adjoint DAE system solution can be eliminated by introducing an extra set of continuous optimization variables δd, in the primal problem, that are fixed according to the equality constraint δd - δκ ) 0. This gives rise to the following primal optimal control problem:

min φ[x˘ i(tf),x(tf),xa(tf),vv,d,δd,tf]

(5)

vv,d,δd

The master problem is clearly a mixed-integer linear program (MILP) because the integer δ and the continuous variables η appear linearly in the problem structure. The solution of the master problem, apart from being lower bound to the MIDO problem, also provides a new integer realization. If the lower bound evaluated at the master problem and the upper bound calculated in the primal cross, then the solution is found and is equal to the primal problem, whereas if they do not cross, the new integer set is added to the primal problem and the algorithm recommences. Summary of the MIDO Algorithm. The steps of the algorithm are briefly summarized below: (1) Fix the values of the binary variables, δ ) δκ, and solve a standard dynamic optimization (DO) problem (eq 3, κth primal problem). An upper bound, UB, is thus obtained. (2) Re-solve the primal problem at the optimal solution (eq 5) with additional constraints of the form δd δκ ) 0, where δd is a set of continuous search variables and δκ is the set of (complicating) binary variables. Convergence is achieved in one iteration. Obtain the Lagrange multipliers, Ωκ, corresponding to the new constraints. (3) Construct the κth relaxed master problem from the κth primal solution, φκ, and the Lagrange multipliers, Ωκ (eq 6). This corresponds to a MILP. The solution of the master, ηκ, gives a lower bound, LB, on the MIDO solution. If UB - LB is less than a specified tolerance or the master problem is infeasible, the algorithm terminates and the solution to the MIDO problem is given by UB. Otherwise, set κ ) κ + 1, set δκ+1 equal to the integer solution of the master, and return to step 1. 2.3. Explicit Model-Based Parametric Controller for Discrete-Time Dynamic Systems. Consider a linear discrete-time state-space description of a dynamic plant:

xt+1 ) A1xt + A2vt yt ) B1xt + B2vt

(7)

yt ∈ Rm is the vector of the output variables that we aim to drive to their set points. The following receding horizon open-loop optimal control problem is formulated in order to derive the explicit model-based control law for such a system:21

s.t. 0 ) fd[x˘ i(t),x(t),xa(t),vv,d,δd,t]

N-1

0 ) fa(•,δd,•), 0 ) f0(•,δd,•)

φˆ (xt|t) ) min xt+N|tTPxt+N|t + vN∈VN

0 g g(•,δd,•)

[yt+k|tTQyt+k|tT + ∑ k)0 T Rvt+k] (8) vt+k

0 ) δd - δ κ

s.t. xt+k+1|t ) A1xt+k|t + A2vt+k yt+k|t ) B1xt+k|t + B2vt+k

t 0 e t e tf The corresponding master problem is then simplified to the following equation:

0 g g(yt+k|t,xt+k|t,vt+k) ) C0yt+k|t + C1xt+k|t + C2vt+k + C3

min η

k ) 0, 1, 2, ..., N - 1

(6)

δ,η

s.t. η g

φ[x˘ κd(tf),xκ(tf),xκa(tf),vκv,dκ,δκd,tf]

+Ω

κ

(δκd

- δ)

κ∈K In the modified equivalent master problem (6), all of the terms are calculated at the solution of the primal problem and no adjoint calculations are required.

0 g ψe(xt+N|t) ) D1xt+N|t + D2 vt|k ) Kxt+k|t

Nek

xt|t ) x* where t is the time when a measurement is taken; xt+k|t denotes the future prediction of the state vector at time


t + k when starting from the state xt|t ) x* and applying T the control sequence vt, vt+1, ..., vt+k-1. vN ) [vTt , vt+1 , ..., T T vt+N-1] denotes the sequence of the control vector over the receding horizon. The constraints g: XV f Rq, ψe: X f RQe and bounds on x, y, and v define the feasible operating region. Q g 0 and R > 0 are symmetric matrices and P g 0 is the terminal cost that satisfies either the Riccati22 (K ) KLQ) or the Lyapunov equation23 (K ) 0) to guarantee stability. The optimal tuning of the input horizon N is usually based on stability and feasibility criteria.24,25 We assume that the pair (A1, A2) is stabilizable and the pair (A1, B1) detectable. By k-1 substitution of xt+k|t ) Ak1x* + ∑j)0 (Aj1A2vt+k-1-j) for the states and treatment of the current states x* ∈ X as parameters, problem (8) is recast as a multiparametric quadratic program (mp-QP) of the form

{

Figure 2. Schematic description of the explicit parametric controller.

1 φˆ (x*) ) min L1 + LT2 v + LT3 x* + vTL4v + v 2 1 (x*)TL5v + (x*)TL6x* (9) 2

}

s.t. G1v e G2 + G3x* The solution of (9)26 consists of a set of affine control functions in terms of the states and a set of regions where these functions are valid. This mapping of the manipulated inputs in the state space constitutes a feedback parametric control law for the system. The mathematical form of the parametric controller is as follows:

vˆ t(x*) ) Acx* + Bc; if CR1c x* + CR2c e 0 for c ) 1, ..., Nc

(10)

where Ac, CR1c and Bc, CR2c are constant vectors and matrices, respectively, the index c signifies that different control expressions apply in different regions in the state space, and Nc denotes the number of such regions. The vector vˆ t is the first element of the optimal control sequence, whereas similar expressions are derived for the rest of the control elements. Note that only the first element of the open-loop control trajectory is implemented on the plant. The control action at the next time interval corresponds to the first element of the control sequence pertaining to the newly updated state realization. In realistic process applications, the state of the system is estimated from a reduced set of measurements. Some of the choices for a state observer include a Kalman filter, an extended Kalman filter,27 and a moving horizon estimator.28 In this work we assume that we have a perfect estimation of the values of the current states. The performance of this model-based parametric controller is optimal in terms of the given performance criteria, the plant model, and the imposed constraints. The implementation of the parametric controller is based merely on simple function evaluations, rather than solving an optimization problem online, which makes it attractive for a wide range of systems. A summary of the principle of parametric controllers is shown in Figure 2.

Figure 3. Binary distillation system for the process and control design example.

3. Problem Formulation To illustrate the simultaneous process and control design paradigm, consider the binary distillation column in Figure 3, which has been primarily studied in the literature.7,29 The model of the system is described in detail in appendix A. The system is subject to a highfrequency sinusoid disturbance in the feed composition with an uncertain mean value. The goal is to obtain the economically optimum process and control design for this system that ensures feasibility, while accounting directly for both continuous and discrete design decisions. The conceptual design framework for this example is outlined in Table 1. The process model is adapted from Schweiger and Floudas29 featuring (i) dynamic material balances for the trays, the reboiler, the condenser, and the reflux drum., (ii) material holdup dynamics and tray liquid hydraulics, (iii) ideal thermodynamics with constant relative volatility, (iv) fast temperature/energy dynamics, and (v) perfect inventory and pressure control. The design optimization problem features a closed-loop dynamic system that incorporates a receding horizon model-based controller. The math-

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4549 Table 1. Formulation of Process and Control Design Optimization for the Binary Distillation Example minimize expected total annualized cost s.t. differential-algebraic process model Inequality Constraints xbenz,D(t) g 0.98, xbenz,B(t) e 0.02 minimum diameter due to flooding Dc - Dcmin(t) g 0 Tout,R(t) - TB(t) g 5 (reboiler), TD(t) - Tout,w(t) g 5 (condenser) Tin,R(t) - Tout,R(t) g 7

purity of the product flooding restrictions minimum temperature iff. subcooling conditions: reboiler

Degrees of Freedom model-based controller tray diameter, Dc; heat-transfer areas, AR (reboiler), AC (condenser); set set points, xset benz,D, xbenz,B Q, R number of trays: reflux tray location ) δrk, feed tray location ) δfk feed composition zbenz,f(t) t0 e t e tf where t0 ) 0, tf ) 720 min

control scheme equations process design variables d control tunings q j discrete decisions δ disturbance/uncertainty time horizon

ematical representation of this problem is

1 min Eθ∈Θ (Ccolumn + Creb + Ccond) + OpCost 3

d,δ,q j ,d h

[

]

s.t. 0 ) fd[x˘ (t),x(t),xa(t),v(t),θ(t),d,δ,t] 0 ) fa[x(t),xa(t),v(t),θ(t),d,δ,t] y(t) ) fy[x(t),xa(t),v(t),θ(t),d,δ,t] j ) 1, ..., Nf, t0 e t e tf (11)

0 g g[x(tj),xa(tj),v(tj),θ(tj),d,δ,tj]

N-1

∑ [(yt+k|t -

h ,δ) ) min xt+N|tTPxt+N|t + φˆ (xt|t,θt,yset,d vN

k)0

yset)TQ(q j ) (yt+k|t - yset) + vt+kTR(q j ) vt+k] s.t.

h ,δ) xt+k|t + A2(d h ,δ) vt+k + W1(d h ,δ) θt xt+k+1|t ) A1(d

yt+k|t ) B1(d h ,δ) xt+k|t + B2(d h ,δ) vt+k + W2(d h ,δ) θt 0 g C0(d h ,δ) yt+k|t + C1(d h ,δ) xt+k|t + C2(d h ,δ) vt+k + h ,δ) C3(d k ) 0, 1, 2, ..., N - 1 0 g D1(d h ,δ) xt+N|t + D2(d h ,δ)

(12)

h ,δ) ) v(t) - vl, xt|t ) x(t) - xl, vˆ t(xt|t,θt,yset,d θt ) θ(t) - θs d h ) d, xl ) xl(d h ,δ), vl ) vl(d h ,δ)

(13)

where

( )[ ( )[

( ) ] ( ) ]

∂fd ∂x˘

-1

∂fd ∂fd ∂fa ∂x ∂xa ∂xa

-1

Ac1 ) -

∂fa ∂x

∂fd ∂x˘

-1

∂fd ∂fd ∂fa ∂v ∂xa ∂xa

-1

Ac2 ) -

∂fa ∂v

∫0∆teA

A1 ) eA 1‚∆t, A2 ) Ac2( c

1‚τ

c

(14)

dτ)

where ∆t is the sampling time that is used to convert the linear continuous-time dynamic system to a discrete-

time representation. W•, B•, C•, W c•, B c•, and C c• are defined accordingly; the superscript c denotes the matrices of the continuous-time linear dynamic system while the matrices A, B, C, and W that do not have this superscript correspond to their discrete-time counterparts.30 d is the vector of continuous process design variables, q j is the vector of control tuning parameters, and δ comprises the binary variables as shown in Table 1. For our example x ) {xB, xD, M1, ..., MNtrays, xbenz,1, ..., xbenz,Ntrays}, where M denotes the tray holdup, xB and xD denote the benzene mole fraction in the top and bottom products, and xbenz pertains to the benzene mole fraction on each tray; y is the vector of output controlled variables, i.e., the top and bottom mole fractions; v represents here the vector of the manipulated variables, i.e., reflux flow rate Refl and boilup flow rate V that are fully determined from the controller equations (12). Note that in this example there are no operating variables vv that can be readjusted during operation. Note also that here we assume that the controller has complete information about the values of the current system states and the values of the disturbances. This information can be obtained from direct measurements or implicitly, via inferential measurements and the use of a Kalman filter estimator. The process constraints g e 0 are shown in Table 1. (xl, vl, θs) is the linearization point of the system, and yset is the vector of output set points. The input disturbance in the feed composition is modeled as

2π t) + θ (t - θ ) (240

zbenz,f(t) ) θ(t) ) θ h + θa sin

i

ω

0.45 e θ h e 0.5, θω ) 248, θa ) 0.08 θi: pulse of magnitude 0.05

(15)

The presence of the uncertainty θ h and the binary variables δr,f makes the design formulation (11)-(13) a stochastic MIDO problem. Note also that this problem involves a bilevel optimization problem, where the process design optimization (11) is the leader and the control optimization problem (12) is the follower (see work by Gumus and Floudas31 for a review of bilevel programming). Brengel and Seider8 formulated a similar problem and relaxed the Karush Kuhn Tucker optimality conditions of the inner optimization problem to enable its solution. This method may not always account properly for switching in the active constraints, while it does not formally address the determination of the control design variables q j . While there are methods


Figure 4. Simultaneous process and control design algorithm featuring model-based predictive controllers. Table 2. Initial Linearization Points for Structure δr26 ) 1, δf12 ) 1 lin. Rl lin. βl

zf

Refll

Vl

xbenz,D,l

xbenz,B,l

AC,l

AR,l

Dc,l

0.45 0.5

3.218 3.199

5.459 5.701

0.98 0.98

0.191 0.191

110 110

280 280

1.65 1.65

for the treatment of uncertainty14 and for solving the underlying MIDO problem7,32 as shown in section 2, there is no efficient method for integrating the receding horizon controller into the simultaneous process and control design framework. In the next section we present an approach for that purpose that is based on parametric programming for deriving the controller structure (see section 2.3) and on an outer approximation technique to determine its design. 4. Solution Procedure A summary of the simultaneous process and control design algorithm is shown in Figure 4. 4.1. Step 0: Initialization. Discretize the uncertainty space θ h ∈ [0.45, 0.5] to two uncertainty scenarios θ h ) {0.45, 0.5}. Choose an initial process structure, e.g., number of trays ) 26 (δr26 ) 1) and feed location ) 12th tray (δf12 ) 1). Select as an initial guess two linearization points Rl and βl shown in Table 2 based on the uncertainty scenarios. Parametrize the diagonal penalty

matrices in the control design (12) to Q(1,1) ) q j 1, Q(2,2) ) 1, and R(1,1) ) R(2,2) ) q j 2, and select as an initial guess for their values q j 1 ) 1.2 and q j 2 ) 10-4. Set iteration counters to κ ) 1 and l ) 1. 4.2. Step I: Structural Primal. 4.2.1. Step I1: Control Design Primal. Linearize the open-loop dynamic model of the process at the two fixed points in Table 2. Then perform model reduction33,34 to derive two reduced linear four-state models. From the Jacobian of the full-state model, balanced truncation was used to reduce the model order. The error introduced by truncating the 2Ntrays + 2-state model with transfer function Hf(jω) to an n ) 4 state representation Hr(jω) is computed from 2Ntrays+2

error )

||Hf - Hr||h∞ ||Hf||h∞

2 )

∑

σH i

i)n+1

||Hf||h∞

(16)

where σH i is the ith Hankel singular value of the original model, while the H-infinity norm of a transfer function is defined as ||H ˆ ||h∞ ) maxω∈R σmax[H ˆ (jω)], where σmax is the maximum singular value of matrix H ˆ . Compute the discrete model matrices, keep the control designs fixed at q j)q j l, and formulate an openloop receding horizon problem (12). For this particular process design, the appropriate values of sampling time and time horizon according to heuristics (Seborg et al.,25


Chapter 27) are ∆t ) 0.3 min and N ) 6, respectively. Note that these values are allowed to change during the solution procedure according to the current design and linearization point. The terminal cost P is the solution of the Lyapunov discrete algebraic equation, while only the purity constraints are considered in the controller design. Instead of solving this problem online in the traditional MPC fashion, parametric programming13 is used to derive the explicit state feedback control law for the plant (see section 2.3). The current states xt|t, the set points yset, and the disturbances θt are treated as parameters and the control inputs as optimization variables; therein, problem (12) is recast as a mp-QP. The solution of this program results in a parametric controller (Parco) for the distillation process that comprises a set of piecewise affine control functions and the critical regions in the state space where these functions hold:

{

[V, Refl]T ) n)4

D

R R set (Ai,c xri) + CRc zbenz,f + ∑ (Di,c xbenz,i) + BRc ∑ i)1 i)B set set if CRRc (xr,zbenz,f,xbenz,D ,xbenz,B ) e 0, R c ) 1, ..., Nc , zbenz,f ∈ [0.37, 0.475]

n)4

∑ i)1

D

β (Ai,c xri)

+

β Ci,c zbenz,f

+

∑

β set (Di,c xbenz,i)

i)B set set if CRβc (xr,zbenz,f,xbenz,D ,xbenz,B ) e 0, β c ) 1, ..., Nc , zbenz,f ∈ [0.475, 0.6]

+

Bβc

|

-4.94101 e 10xr1 +

+

s.t. 0 ) fd(x˘ i(t),xi(t),xia(t),[V(t),Refl(t)],d,θ h i) 0 ) fa(•), 0 ) f0(•), 0 g g(•), i ) 1, 2 where [V(t)i, Refl(t)i] is given in (17)

xir ) (TiLT)T‚xi

(17)

where xr are the reduced states. For instance, in the region described by the following inequalities set 0.0076xbenz,D

Note that (17) replaces exactly (12). The next step is to substitute (17) into (11) and (13) and treat only d as an optimization variable in the solution of the resulting multiperiod primal problem:

where TLT is a matrix derived from the model reduction that represents the mapping of the real states to the reduced states. The solution of the DO problem (18) provides an upper bound UPc inasmuch as the fixed controller tunings and the linearization point values are not necessarily optimal. 4.2.2. Step I2: Control Design Master. Compute the gradients of the objective function and the constraints with respect to the (i) linearization point d h , (ii) control designs q j , and (iii) process designs d to construct and solve the following control master problem:

min η q j ,d h ,d

set 0.0105xbenz,B e

+5.54101

s.t. η g φl +

set set -15.5127 e +15xr2 - 0.0418xbenz,D + 0.0200xbenz,B e

-14.9127 + 20xr3 -

set 0.0575xbenz,D

+

e +13.5137

set set - 0.5958xbenz,B e -204.557 + 10xr1 + 0.9655xbenz,D 763.158xr2 - 735.539xr3 - 613.728xr4 set set - 0.5202xbenz,B e370.639zbenz,f + 1.7432xbenz,D 2064.91

T

l

l

T

l

l

T

set 0.0375xbenz,B

set set 27.9961 e +30xr4 - 0.1219xbenz,D + 0.0923xbenz,B e +29.1961 + 20zbenz,f e +12 + 53.1505xr1 + 267.148xr2 - 48.4428xr3 + 30xr4 + 60.3675zbenz,f +

dφ dφ (d h-d h)+( (q j-q j)+ (dd h ) dq j ) dφ (dd ) (d - d ) (19)

0 g gl +

l

l

dg dg (d h-d h)+( (q j-q j)+ (dd h ) dq j ) dg (dd ) (d - d ) (20) T

l

l

T

l

l

T

l

e g |dLk - d h |, g |vLv k - vj v|

l

(21)

l ) 1, ..., Lκ The control expression has the form

V ) +6.10246xr1 + 30.6725xr2 - 5.56195xr3 + set 3.44444xr4 + 6.93108zbenz,f + 110.69797xbenz,D set + 29.7862 68.2136xbenz,B

Refl ) -37.144xr1 + 52.495xr2 + 3.36925xr3 + set set - 44.0635xbenz,B + 6.32022xr4 + 190.3907xbenz,D 8.03056zbenz,f + 63.1623

where Lκ is the number of feasible control primal problems that have been solved so far. Inequalities (21) ensure that the new process linearization point in terms of d h is close enough to the process design point that is derived by the previous control primal problem. These inequalities are, in fact, a relaxation of the equalities (13). Note that if the primal problem is infeasible, a constraint minimization problem is solved instead. Accordingly, in the master problem the inequality pertinent to the objective is omitted. An example of a


control design master problem for the structure number of trays ) 25 and feed location ) 12, is shown here:

min

set set set set q j 1,q j 2,A h C,A h R,D h c,xjbenz,D ,xjbenz,B ,AC,AR,Dc,xbenz,D ,xbenz,B

η

set - 0.98038) - 2.5e η g 6.2283 + 13.159(xjbenz,D set - 0.0156) + 1e - 4(A h C - 115) + 1e 2(xjbenz,B 3(A h R - 268) + 1e - 5(D h c - 1.64) + 4.04e - 2(q j1 1.2) + 6.05e - 3(q j 2 - 1e - 4) + 3.03e - 3(AC 131.92) + 1.21e - 2(AR - 277.09) + 1.21e - 4(Dc set set - 0.9815) + (-4.6)(xbenz,B 1.66) + 10(xbenz,D 0.01782)

set - 0.0156) - 1e - 4(A h C - 115) + 1e 3.1152(xjbenz,B h c - 1.64) - 1.48e - 3(q j1 4(A h R - 268) - 1e - 5(D 1.2) + 4.29(q j 2 - 1e - 4) + 1e - 7(AC - 131.92) 1e - 8(AR - 277.09) + 2.96e - 1(Dc - 1.66) set set - 0.9815) - 3.77e - (xbenz,B - 0.01782) 90.6(xbenz,D

l set 0 g 1.4866e - 3 - 1.7912(xjbenz,D - 0.98038) + set - 0.0156) + 1e - 6(A h C - 115) + 7.635e - 2(xjbenz,B h c - 1.64) - 1.3e 1e - 7(A h R - 268) + 1e - 6(D 5(q j 1 - 1.2) + 5.13(q j 2 - 1e - 4) + 3.7e - 7(AC 131.92) - 2.82e - 7(AR - 277.09) + 1.2e - (Dc set set - 0.9815) - 4.92e - (xbenz,B 1.66) + 3.75(xbenz,D 0.01782)

Note that problem (19) is a linear program (LP) and its solution provides a lower bound to the structural primal and a new realization for d h and q j . If UPcl e LOcl , then stop; the solution is equal to the upper bound. s , Update the structural upper bound UPsκ ) min{UPκ-1 c c c j l+1 ) q UPl } and go to step II. If UPl > LOl , then set q j and d h l+1 ) d h and update the counter l ) l + 1. Then go to step I1 and update the values of the control tunings and the linearization point. 4.3. Step II: Resolve Session of Structural Primal. Introduce a set of continuous variables δd that replace the binaries while adding the equality constraints ψeδ ) δd - δ to the formulation (18). Resolve the structural primal in one iteration to get the multipliers Ωκ of the extra constraint.7 4.4. Step III: Structural Master. Use the multipliers Ωκ to formulate the master problem:

(22)

min η

Ntrays

η g TotalCost + κ

r f [Ωrl (δl,d - δrl ) + Ωfl(δl,d - δfl)] ∑ l)1 k

k

k

Ntrays

1)

δfk ∑ k)1

Ntrays

1) 0 g δfk -

∑ Ntraysδk′r

k′)k

χ ) max {χl|χl ) l)1,...,q

set 0 g 4.24199e - 4 + 3.2745(xjbenz,D - 0.98038) +

ηδkr ,δkf

where Kfeas is the number of feasible primal solutions. The solution of the MILP master problem (22) is a lower bound LOs for the design MIDO problem and provides a new integer realization. If LOs g UPs, then stop; the solution of the MIDO problem corresponds to the upper bound; else, go to step I; set κ ) κ + 1 and update the integer values. 4.5. Step IV: Feasibility Test. Check if there are any constraint violations for all θ h ∈ [0.45, 0.5]. Note that there are no operating variables in this process because the output set points are treated as design variables. Hence, the feasibility problem (2) reduces to

δrk ∑ k)1

k ) 1, ..., Ntrays, κ ∈ Kfeas}

max

θ h ∈[0.45,0.5]

gl(‚)}

(23)

If χ < 0, the optimal design remains operable for all possible uncertainty realizations. Otherwise, augment the critical uncertain values that bottleneck the system feasibility to the set of uncertainty scenarios and return to step I. 4.6. Remarks. (a) The key feature in the solution of the structural primal problem is the derivation of the parametric controller that enables the elimination of the online control optimization. The method employed for the tuning of this controller as described here is based on an outer approximation technique, where the control tuning parameters and the linearization point are viewed as the complicating variables. This approach generates recursively a set of supporting functions of the form of (19)-(21) for the structural primal (18). The control master problem is thus a lower bound to the structural primal under certain assumptions that are discussed in appendix B. (b) The design obtained with our approach is optimal with respect to the economic performance index and also guarantees operability in the presence of time-varying uncertainties. (c) The parametric controller (17) has an explicit form, but it involves “if-then” statements. In the optimization problem (18), these Boolean algebra operations are replaced by steep hyperbolic tangent functions as shown in appendix C. (d) The control structure selection, i.e., the optimal pairing between the output and the manipulated inputs, is accounted for in this work via the controller design. For instance, if at the end of step I2 the control tuning q j 1 is approaching zero, it implies that the output xbenz,D should not be participating in the control structure. In such a case, in the next control design primal problem q j 1 would be fixed to zero. The next two sections elucidate in detail the application of the algorithm, described in this section, to two illustrative examples. 5. Evaporation Process Example This example is concerned with deriving the optimal process and control design of an evaporation process from the literature.35,36 The schematic representation of the process is shown in Figure 5. The simultaneous process and control design optimization problem for this example is posed in Table 3. 5.1. Objective Function. The total annualized cost is partitioned into the following term:

Ctot ) Ccap + Cop

(24)

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4553 Table 3. Formulation of the Process and Control Design Optimization for the Evaporation Process minimize total annualized cost s.t. differential-algebraic process model36,50 purity of product, pressure of product, utility consumption bounds: 25% e C2(t) 40 kPa e P2(t) e 80 kPa 0 kg/min e F200 e 400 kg/min (85) 0 kPa e P100 e 400 kPa control scheme equations: parametric controller design variables: condenser area 200 m2 g Ac g 10 m2, heater area 200 m2 g Ae g 18 m2, controller tunings Q g 0, R g 0 operating variables: controller set points set 25% e Cset 2 , 40 kPa e P2 e 80 kPa feed conditions: F1(t), C2(t) time horizon: t0 e t e tf where t0 ) 0, tf ) 1440 min

inequality constraints

degrees of freedom

disturbances

idea in this approach is to define a variable Χ h that tracks the maximum value of the path constraint Χ h (t) e 0 over time and then to enforce an end-point constraint on this new variable. Hence:

Χ h (t) ) max χ(th) t∈[t0,t]

χj(tf) e 0

(27)

The optimization problem (27) is approximated by

dΧ h dΧ dΧ ) 0.25[1 + tanh(106Float)] 1 + tanh 106 dt dt dt Float ) Χ - Χ h (28) Float(t0) ) 0

[

(

)]

If χ corresponds to a differential state, then its derivative is explicitly available for use in (28). Otherwise, dΧ/ dt can be accurately approximated by dΧ ˆ /dt, where

Figure 5. Flowsheet of the evaporation process.

The total annualized capital cost is37

Ccap ) Ce + Cc Cc ) 0.6(1/3)(M&S/280) ×

(

)

0.65 Ac × 3.29 × 1.35 101.3 104 144 × 2.542 Ce ) 0.6(1/3)(M&S/280) × 0.65 Ae × 3.29 × 1.00 (25) 101.3 104 144 × 2.542

(

)

where Ce is the cost of the heater and Cc is the cost of the condenser. The annualized operating cost is

Cop )

t (0.6F200 + 1.009(F2 + F3) + ∫t)0 f

8150 600F100) dt (26) tf × 1000 where F100 is the flow of the steam, F200 the flow rate of cooling water, and F2 + F3 the recycle flow rate. The Marshall and Swift (M&S) index for year 2001 when this example was solved is 1100.38 5.2. Inequality Constraints. The inequality path constraints (85) are modified to end-point constraints via the method developed in Bansal et al.7 The basic

Χ)Χ ˆ + 10-8

|

dΧ ˆ dΧ ˆ , dt dt

)0

(29)

t)t0

5.3. Disturbances. The disturbances that trigger the dynamic behavior have a sinusoid variation around their nominal point (F1l ) 10 kg/min, C1l ) 5%) described by the following equations:

( ) ( )

F1 ) 10 + 0.8 sin

2π t tf/4

C1 ) 5 + 0.5 sin

2π t tf/3

(30)

No bounded unknown uncertainty is considered in this example because the disturbances have a fixed amplitude, mean value, and frequency. 5.4. Decision Variables. In this problem no discrete process design decisions are considered. The value of the concentration set point varies periodically between 25% and 27% as a result of operating decisions. This variation is allowed to occur between time t1 ) 360 min and t2 ) 740 min on a daily basis. This transition is modeled as a piecewise linear variation, and its duration is a free optimization decision. To derive the parametric controller, a discrete-time model-based predictive control problem is formulated based on the linearized version of the open-loop model. The manipulating variables v ) {P100, F200} are the control optimization variables, whereas the outputs y ) {C2, P2} correspond

4554 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 4. Results of the Design Decomposition Approach for the Evaporation Example variables F200-BIAS (kg/min) P100-BIAS (kPa) Ae (m2) Ac (m2) cap. cost (104$/yr) oper. cost (104$/yr) obj. (104$/yr) UP

initial guess

iteration 1

Control Primal 207.515 189.023 193.372 211.471 19.2 18 13.68 13.90 1.7631 4.62960 6.39278 +∞ 6.39278

iteration 2 189.061 213.261 18 13.81 1.7603 4.63187 6.39217 6.39217

iteration 3 182.851 217.555 18 13.76 1.7587 4.63347 6.39209 6.39209

Control Master control designs q j1 q j2 linearization point F200 (kg/min) P100 (kPa) Ae (m2) Ac (m2) obj (104$/yr) LO is UP-LO < 0?

1 10-3

100 0.24

55.32 0.2519

207.515 193.372 19.2 13.68

198.676 185.406 19 14 6.3849 6.3849 no

212.838 180.602 18 14 6.39215 6.39215 no

-∞

yes

to the system states. The objective function is a quadratic index of performance as in (12) where

[

]

[ ]

q j /(C - 30) 0 1 0 Q ) Q ) 1 2l , R)q j 2‚ 1/(P2l - 80) 0 1 0 (31) where (P2l, C2l) is the linearization point. In this study we consider q j 1 and q j 2 the control design variables to be determined by our optimization procedure. The discrete sampling time is taken as ∆t ) 1 min, and the prediction horizon is considered to be N ) 3. The terminal cost P is computed from the solution of the Lyapunov equation. The constraints enforced on the prediction horizon are shown in Table 3, eq 85. The model-based control optimization problem is recast as a mp-QP by treating the current states, the set points, and the disturbances as a set of parameters. The resulting state feedback parametric controller has a piecewise affine structure. 5.5. Results. The simultaneous process and control design problem in Table 3 is solved via the outer approximation decomposition algorithm described in section 4. Convergence was achieved in three iterations, and the results are shown in Table 4. The initial guess corresponds to the steady-state optimal point. The profiles for the set points in each iteration are shown in Tables 5 and 6. The optimal controller for this process is portrayed in Figure 6, in the two-dimensional state space. Each colored region corresponds to a different function for the control variables in terms of the states. The evolution of the system over the time horizon of interest is also shown in Figure 6. The time trajectories of the outputs and the inputs are shown in Figure 7. 5.6. Remarks. (a) A large value for q j 1 implies tight control on the product concentration and favors constraint satisfaction. However, there is a tradeoff because it has an adverse effect on the economic objective. (b) A high q j 2 implies tight control on the energy input, thus favoring economic performance. The end result contradicts in this particular example the common MPC heuristic39 of imposing a small weight on the input deviation. Clearly, the reason is that economics and not tight output control is the performance requirement in this problem.

Table 5. Optimal Concentration Set-Point Variability in the Design Decomposition Approach for the Evaporation Example concn of Cset 2 (%) at the beginning of the interval

concn of Cset 2 (%) at the end of the interval

interval

duration of the interval

1 2 3 4 5

360.000 18.13 341.87 20 700.000

Iteration 1 25.030 25.030 27 27 25.053

25.030 27 27 25.053 25.012

1 2 3 4 5

360.000 14.84 345.16 20 700.000

Iteration 2 25.004 25.004 27 27 25.016

25.004 27 27 25.189 25

1 2 3 4 5

360.000 14.84 345.16 20 700.000

Iteration 3 25.005 25.005 27 27 25.074

25.005 27 27 25 25.008

Table 6. Optimal Value of the Pressure Set Point in the Design Decomposition Approach for the Evaporation Example

pressure

Pset 2

(kPa)

iteration 1

iteration 2

iteration 3

51.260

51.744

52.840

(c) The solution, for the given capital and operating cost objective, settles to a point where a tight, almost perfect, control is enforced on the product composition, whereas the other output, i.e., product pressure, is loosely controlled. One advantage of this technique is that it provides two more feasible solutions to the operator apart from the optimal. Thus, the operator has the ability to select between an extremely tight control on the concentration (primal 2 and 3) and an equally distributed control performance on the two outputs (primal 1). 6. Distillation Process Example The distillation process example (Figure 3) that served as an illustration of the theoretical developments in this paper is revisited here. 6.1. Results. Two uncertainty periods were selected as discussed in section 4, and the multiperiod design MIDO problem is solved first. The algorithm converged in three iterations between the structural master and primal solutions, as shown in Table 7. Each structural primal problem requires two iterations between the control primal and master solutions, as shown in Table 8 for the optimal structure. In all three MIDO iterations, the control master determines for every structure that matrix Q is positive-definite, i.e., that both outputs participate in the optimal control structure. In the process structure, there are 30 possible feed tray locations, and for the feed located on tray k, there are 31 k alternatives for the reflux tray locations. Hence, the 30 (31 - k) ) total number of discrete alternatives is ∑k)1 465. Despite this large number of alternative discrete decisions, the algorithm converged in three iterations between the structural master and primal problems. A reduced model of four states was used in all three process structures for control design purposes. The error from the reduction was error ) 0.17% for the optimal structure, while the frequency responses of the singular


Figure 6. Schematic description of the parametric controller for the evaporation example.

Figure 7. Time trajectories of the inputs and outputs in the evaporation example.

values of the reduced versus full model are shown in Figure 8. Note that for frequencies above 100 rad/min the responses of the two models deviate. However, in this study we are interested in disturbance frequencies of less than ω e 1-10 rad/min, where the reduced model portrays the dynamic behavior satisfactorily. The profiles of the control inputs and outputs are depicted in Figure 9, whereas the time trajectory for the minimum allowable column diameter is shown in Figure 10. After the solution of the simultaneous process and control design MIDO problem, the feasibility problem (23) is solved. Its solution is shown in Table 9. The maximum constraint value χ ) -8.6 × 10-3 corresponds to the bottoms benzene mole fraction and is less than

zero. Thus, the optimal process and control design obtained in Table 7 is indeed feasible for all possible bounded uncertainty scenarios. 6.2. Remarks. (a) Table 10 compares different designs to demonstrate the benefits of pursuing a simultaneous approach in process and control design rather than the traditional sequential approach. The steadystate design in column 1 (Table 10) is feasible at steady state but inoperable at transient conditions because it cannot satisfy the specifications. In an attempt to make the steady-state design operable, the equipment size was increased by 10%, resulting in a still inoperable process as a result of violations of composition and thermodynamic driving force constraints. The equip-

4556 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 7. Progress of the Iterations for the Multiperiod MIDO Design Problem in the Distillation Example iteration 1 no. of trays feed location control scheme process design Dc (m) AR (m2) AC (m2) controller tunings q j1 q j2 set point xset benz, D set xbenz, B costs capital cost ($100K yr-1) operating cost ($100K yr-1) expected cost ($100K yr-1) UP no. of trays feed location control scheme LO a

iteration 2

iteration 3

Primal Solutions 26 25 24 12 12 12 full structurea full structure full structure 1.654 275.15 131.61

1.657 276.16 131.95

1.663 278.74 133.63

16 10-4

15 10-4

11 10-4

0.9815 0.01782

0.98146 0.01770

0.98149 0.01716

1.9953

1.9720

1.9543

4.2327

4.2507

4.2914

6.2280

6.2227

6.2457

6.2280

6.2227

6.2227

Master Solutions 25 24 12 12 full structure full structure 6.2118 6.2205

25 11 full structure 6.2339

Full structure: xbenz,D, xbenz,B - Refl, V.

Table 8. Progress of the Iterations in the Control Design Procedure for Structure Ntrays ) 25, Feed Location ) 12 (Distillation Example) iteration 1 xbenz,D xbenz,B AC AR Dc

Linearization Point 0.98038 0.0156 115 268 1.64

iteration 2 0.98025 0.01916 110 280 1.65

zbenz,f Vl Refll

Lin R 0.45 5.530 3.280

0.45 5.480 3.238

zbenz,f Vl Refll

Lin β 0.5 5.725 3.224

0.5 5.238 3.223

Control Primal 0.9815 xset benz,D set 0.01782 xbenz,B 131.92 AC AR 277.09 1.66 Dc q j1 1.2 q j2 10-4 cap. cost ($100K yr-1) 1.9738 oper. cost ($100K yr-1) 4.2545 exp. cost ($100K yr-1) 6.2283 UP 6.2283 q j1 q j2 LO is LO > UP?

Control Master 15 10-4 6.2218 no

0.9814 0.0177 131.95 276.16 1.65 15 10-4 1.9720 4.2507 6.2227 6.2227 5 10-4 6.2360 yes

ment size was then increased by 20%, leading to an operable but expensive design presented in the second column of Table 10. As opposed to this ad hoc sequential overdesign procedure, the systematic optimization methods of columns 3 and 4 lead to a selective increase in

Figure 8. Singular values of a full 52-state model (Ntrays ) 25 and feed location ) 12) vs a reduced 4-state model for the distillation example. Table 9. Feasibility Test Results in the Distillation Example constraint

χl

critical value of uncertainty θ h

0.98 - xbenz,D e 0 -(0.02 - xbenz,B) e 0 7 - Tin,R(t) + Tout,R(t) e 0 5 - TD(t) + Tout,w(t) e 0 5 - Tout,R(t) + TB(t) e 0 -Dc + Dcmin e 0

-1.2 × 10-3 -8.6 × 10-3 -1.1 × 10-2 -7.6 × 10-1 -2.2 × 10-1 -2.7 × 10-3

0.45 0.45 0.5 0.45 0.45 0.5

the equipment size of the reboiler and the column diameter, whereas the condenser size remains almost the same. As a result, economic advantages on the order of 3% (column 3 vs column 2) to 5-6% (column 4 vs column 2) are achieved with feasible dynamic performance. Thereby, it is guaranteed that the resulting design meets all of the production specifications and operational restrictions despite the presence of rapidly varying disturbances. (b) The economically optimum design point usually lies at constraint intersections.40 The parametric controller is particularly effective in dealing with constraints, while it contains a feedforward element to compensate for the disturbance effect. These features allow the plant to operate closer to the constraint limit, as opposed to the operating point derived when using PI control. This property, in combination with the reduced equipment size, explains why the parametric controller leads to total economic benefits of 2-3% compared to the PI-based design as shown in columns 3 and 4 of Table 10. Another significant benefit of the parametric controller is shown in Figure 11, where we examine the scenario of an increased disturbance amplitude and impulse (θa ) 0.095, θi ) 0.07, and θω ) 228 min in (15)). The design with the parametric controller exhibits half the size of the overshoot compared to the PI controller, while it avoids underdamped oscillatory behavior. This demonstrates that this novel control law respects the process constraints, thus enhancing operational plant performance. (c) The employment of advanced controllers in the design framework enables the direct accommodation of stability performance criteria in the controller synthesis phase. This is done via computing the terminal cost and the time horizon length in (17) according to established literature criteria.23-25,41 This is a clear advantage over


Figure 9. Profiles of the manipulating inputs and controlled outputs for the distillation example, for θ h ) 0.45. Table 10. Comparison of Different Designs in the Distillation Example

design variables no. of trays feed location Dc (m) AR (m2) AC (m2) nominal utility consumption Fw (kg/min) Fst (kg/min) set points xset benz,D set xbenz,B controller tunings KcxD-Refl KcxB-V τcxD-Refl τcxB-V q j1 q j2 costs capital cost ($100K yr-1) operating cost ($100K yr-1) expected cost ($100K yr-1)

steady-state inoperable

sequential approach SISO-PI

25 12 1.6061 255.34 133.66

25 12 1.93 306.43 156.41

970 83

922 85

0.98 0.02

0.985 0.01009 50 -50 80.03 3.55

simultaneous approach SISO-PI

simultaneous approach Parco

26 13 1.68 289.06 132.55

25 12 1.65 276.16 131.95

1029 84 0.985 0.0099

1004 83 0.9814 0.0177

50 -50 80 5.59 15 10-4

1.9076 4.2304 6.1380

the PI case, where stability either is not accounted for7,14 or is guaranteed via overconservative complex constraints.42 7. Conclusions In this paper, an optimization strategy is considered for simultaneous process and control design using advanced parametric controllers. The explicit control functions that comprise the controller allow its direct incorporation in the design framework. Our approach is based on a parametric programming technique for deriving the controller structure, combined with an outer approximation decomposition method for simultaneously identifying the optimal process and control design decisions. The presence of discrete decisions

2.274 4.3024 6.5764

2.034 4.3353 6.3691

1.9720 4.2507 6.2227

about the system design and operation is accommodated by formulating and solving a MIDO problem, while in the face of parametric uncertainty the two-stage decomposition algorithm of Mohideen et al.14 and Bansal et al.7 is adopted. The advantage of the algorithm is that it introduces formally for the first time advanced optimizing modelbased controllers in a complete design framework. The benefits from this approach include (i) improved economic performance, (ii) enhancement of the system’s dynamic performance, (iii) guaranteed operability in the face of uncertainties, and (iv) improved system stability characteristics. Our method indicates that simultaneous process and advanced control design can lead to economically attractive designs with guaranteed feasible operation.

4558 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 Table 11. Major Variables of the Distillation Example variable

Figure 10. Profiles of the minimum allowable column diameters for θ h ) 0.45 and 0.5 for the distillation example.

Figure 11. Mole fraction vs time for an aggressive disturbance realization. (distillation example).

feed flow rate in tray k total reflux flow rate reflux flow rate in tray k vapor flow rate leaving tray k liquid flow rate leaving tray k vapor boilup distillate flow rate bottoms flow rate vapor mole fraction of component i in tray k liquid mole fraction of component i in tray k distillate mole fraction bottoms mole fraction feed model fraction of component i binary indicating the existence of tray k binary determining whether tray k is the feed tray molar holdup of component i in tray k total molar holdup in tray k molar holdup of component i at the bottom of the column molar holdup of component i in the reflux drum volumetric hold up in tray k molar density of the liquid mixture in tray k minimum flooding diameter tray diameter level of liquid in tray k length of weir tray area heat-transfer area of the reboiler heat-transfer area of the condenser vapor pressure of component i in bottoms/distillate mean relative volatility temperature of utility stream leaving the reboiler temperature of utility cooling water leaving the condenser bottoms temperature distillate temperature heat input into the condenser heat input into the reboiler column installation cost reboiler installation cost condenser installation cost operating cost total cost

units Fk Refl Rk Vk Lk VB D B yi,k

kmol/min . . . . . . . dimensionless

xi,k

.

xi,D xi,B zi,f δrk

. . . .

δfk

.

Mi,k Mk Mi,B

kmol . .

Mi,D

.

Volk FLmix,k

m3 kmol/m3

Dcmin Dc Levelk Lweir Atray AR AC P°i,B/D

m . . m m2 . m2 kPa

a Tout,R

dimensionless K

Tout,w

.

TB TD QD QB Ccolumn Creb Ccond OpCost TotalCost

. . kJ/min . $ $ $ $/yr $/yr

help in software implementation issues. The financial support from the Department of Energy Transport and the Regions (DETR-ETSU), Shell Chemicals, and Air Products and Chemicals Inc. (APCI) is gratefully acknowledged. Appendix A: Model Description for the Binary Distillation Column Example

Figure 12. Distillation process superstructure.

Acknowledgment The authors thank Dr. Vivek Dua, Dr. Myrian Schenk, and Dr. Nikos Bozinis for useful discussions and their

The process superstructure is shown in Figure 12. The integer modeling to represent the discrete process alternatives is performed along the lines of Bansal et al.,7 whereas the DAE dynamic model is formulated according to Schweiger and Floudas29 and Kookos and Perkins.43 The components of the binary mixture are toluene and benzene. The key features involved in this model are as follows: (i) dynamic material balances for trays, reboiler, condenser, and reflux drum, (ii) dynamic liquid material holdups, (iii) liquid hydraulics for each tray, and (iv) flooding correlations. The assumptions of this model are the following: (i) perfect control on reflux


Equations Involving Binary Variables.

Table 12. Major Variables of the Distillation Example parameter

value (units)

feed flow maximum no. of trays mass density of benzene mass density of toluene molecular weight of benzene molecular weight of toluene height of weir reboiler and reflux drum holdup latent heat of the bottoms mixture latent heat of the distillate mixture latent heat of steam height of the distillation column distance between two trays inlet temperature of the steam inlet temperature of cooling water heat capacity of condensed steam heat capacity of cooling water heat-transfer coefficient at the condenser heat-transfer coefficient at the reboiler

F N Fbenz Ftol MWbenz MWtol Hweir MB, MD

5 kmol/min 30 885 kg/m3 867kg/m3 78.108 kg/kmol 92.134 kg/kmol 5 cm 1.75 kmol

∆HB

359‚92.134 kJ/kmol

397.4‚78.108 kJ/kmol ∆Hst 2182 kJ/kg Heightstack Space(∑(kδrk) - 1) (ft) Space 2 ft Tin,R 400 K

Component molar balances for each tray: Ntrays

(

k′)1

298 K

CpwB

4.24 kJ/kg‚K

Ntrays

r δk′ ) ∑ k′)1

(

UB

) Lk+1xi,k+1 + Vk-1yi,k-1 + Fkzi,f + dt Rkxi,d - Lkxi,k - Vkyi,k, i ) benz, k ) 2, ..., Ntrays - 1 (41)

dMi,1

) L2xi,2 + VByi,B + F1zi,f + R1xi,d -

dt

L1xi,1 - V1yi,1, i ) benz (42)

Last tray from bottom Ntrays

∑

( Cpw UC

dMi,k

First tray from bottom

∆HD

Tin,w

∑

r δk′ )

4.18 kJ/kg‚K 0.7 × 60 kJ/ K‚m3‚min 0.9 × 60 kJ/ K‚m3‚min

r δk′ )

dMi,Ntrays

) VNtrays-1yi,Ntrays-1 + FNtrayszi,f + dt RNtraysxi,d - LNtraysxi,Ntrays - VNtraysyi,Ntrays,

k′)Ntrays

i ) benz (43)

drum and reboiler via product flow rates, (ii) perfect pressure control via coolant flow, (iii) constant relative volatility, (iv) fast temperature/energy dynamics. The major variables in this process are shown in Tables 11 and 12. Integer Decisions.

Total molar balance (this is a difference from the model in Schweiger and Floudas29 where constant molar holdup is assumed): Ntrays

(

∑

r δk′ )

k′)k

dMk dt

) Lk+1 + Vk-1 + Fk + Rk - Lk - Vk, k ) 2, ..., Ntrays - 1 (44)

Feed location integers: Fk ) Fδfk, k ) 1, ..., Ntrays

(32) First tray from bottom

Ntrays

δfk ) 1 ∑ k)1

Reflux location integers: Rk ) Refl ×

δrk,

Ntrays

(33)

(

dM1

r δk′ ) ∑ k′)1

dt

) L2 + VB + F1 + R1 - L1 - V1

(45)

Last tray from bottom k ) 1, ..., Ntrays

(34)

Ntrays

(

∑

k′)Ntrays

Ntrays

δrk ) 1 ∑ k)1

r δk′ )

dMNtrays dt

) VNtrays-1 + FNtrays + RNtrays LNtrays - VNtrays (46)

(35) Other Equations.

Feed location only below reflux: Ntrays

δrk -

r δk′ e 0, ∑ k′)k

k ) 1, ..., Ntrays

(36)

Vk ) Vk-1 ) VB, k ) 2, ..., Ntrays

Control structure selection integers are as follows:

δc(ji) ) 1 if both vj and yi belong to the control structure (37) c

else δ (ji) ) 0, j ) 1, 2, i ) 1, 2 where v1 ) V - Vl, v2 ) Refl - Refll

Constant vapor molar flow

(38) (39)

y1 ) xbenz,D - xbenz,D,l, y2 ) xbenz,B - xbenz,B,l (40)

(47)

Holdup equation Mk,i ) Mkxk,i

(48)

Volk ) Mk/FLmix,k

(49)

Volk ) (Atray/3.32)Levelk

(50)

i ) benz, k ) 1, ..., Ntrays


Reflux drum molar balances

Francis’ weir formula Levelk )

Lk2/3 + Hweir 1.84FLmix,kLweir

(51)

dMi,D ) VNtrays(yi,Ntrays - xi,D) dt

(65)

Ntrays


0 ) VD -

∑ Rk - D

(66)

k′)1

MD ) Mi,D/xi,D

Density of liquid mixture 1 FLmix,k

)

xi,k 1 - xi,k + Fi Ftol MWbenz MWtol

i ) benz (52) Reboiler energy balance (subcooling in steam; details of the heat exchangers models are not incorporated in the Scheiger and Floudas29 and Kookos and Perkins43 models; the equations here are taken from Bansal et al.7)

i ) benz, k ) 1, ..., Ntrays Minimum tray diameter due to flooding Dcmin ) 0.6719V0.5

QB ) V∆HB

(53) Fst )

Tray area and weir length 2

Atray ) 0.8πDc /4

(54)

Lweir ) 0.77Dc

(55)

QB ∆Hst + CpwB(Tin,R - Tout,R)

QB ) UBAR

Reboiler vapor-liquid equilibrium (bubble point) 1)

P°benz,Bxi,B + P°tol,B(1 - xi,B) P

(56)

Condenser vapor-liquid equilibrium (dew point)

(

)

1 - xi,D xi,D + P°benz,D P°tol,D

P°i,D ) P°i,D(TD)

(58)

P°i,B ) P°i,B(TB)

(59)

i ) benz, tol Constant relative volatility R)

x

yi,k )

P°benz,DP°benz,B P°tol,DP°tol,B Rxi,k

1 + xi,k(R - 1)

(Tin,R - TB) - (Tout,R - TB) Tin,R - TB ln Tout,R - TB

(70)

QD ) FwCpw(Tout,w - Tin,w)

(71)

QD ) UcAc

Vapor pressure (Antoine equation)44

(TD - Tin,w) - (TD - Tout,w) TD - Tin,w ln TD - Tout,w

First-order delay in bottoms and distillate composition measurements and in the estimation of the rest of the states dx´ i,B xi,B - x´ i,B ) dt τB

(72)

i ) benz (60) (61)


dx´ i,D xi,D - x´ i,D ) dt τD

(73)

i ) benz dx´ i,k xi,k - x´ i,k ) dt τk

Reboiler molar balances dMi,B ) L1xi,1 - Bxi,B - VByi,B dt 0 ) L1 - B - VB

(62)


(63)

MB ) Mi,B/xi,B

(64)

dM Ä k Mk - M Äk ) dt τk

i ) benz

(69)

QD ) V∆HD

(57)

i ) benz

(68)

Condenser energy balance (assuming no subcooling in the process stream)

i ) benz

1)P

(67)

k ) 1, ..., Ntrays

(74)

(75)


1. (Notation). Let the vectors yˆ and vˆ be defined as

Initial conditions (steady state)

|

d(.) dt Cost

)0

yˆ ) [d h, q j ]T, vˆ ) [d]T

(76)

t)0

After the appropriate substitutions, let problem (11)(13) for fixed binaries be reformulated as follows:

Correlations7,37.

{

min φ h (yˇ ,vˇ )

Cost of the column Ccolumn ) 1.4 × 3.18 × M&S 100Dc 1.066 (Heightstack + 15)0.802 + 101.9 280 12 × 2.54 100Dc 1.55 M&S Heightstack (77) 1.4 × 4.7 280 12 × 2.54

(

)

(

Cost of the reboiler

)

(

104AR M&S 280 144 × 2.542

Creb ) 0.6 × 101.3

Cost of the condenser

(

4

)

10 Ac

0.65

× 3.29 × 1.35

)

M&S 280 144 × 2.542

Ccond ) 0.6 × 101.3

|

h ) in which the BNLP(yˆ j) is defined as the problem (P complicating variables yˆ are fixed yˆ j. We aim here to find a representation of P h as an intersection of an infinite collection of sets. This is performed by means of supporting hyperplanes providing a polyhedral representation of the continuous feasible space of program (P h ). Such a representation implies linearity in the optimization variables enabling the replacement of the difficult (11)-(13) program with a LP:

{

min η

0.65

× 3.29

Annual operating cost OpCost )

yˇ ,vˇ

(P h ) s.t. 0 e gj (yˇ ,vˇ ) yˇ ∈ yˇ , vˇ ∈ vˇ

∫t)tt)t [Fst × 0.01 × 8150 × 60 + Fw × 4 × f

0

10-5 × 8150 × 60] dt/tf Objective: total annual cost 1 TotalCost ) OpCost + (Ccolumn + Creb + Ccond) 3 Appendix B: Outer Approximation Approach for the Tuning of the Explicit Model-Based Parametric Controller It is interesting to note that by treating the set of vectors (d h, q j , x*) in problem (12) as parameters we can recast problem (12) as a parametric program. The solution of this problem will provide explicit expressions for the control action in terms of the state x* and the tuning parameters (d h and q j ). However, the mapping d h, q j f v is, in general, nonconvex and nonlinear.45 Thus, its parametric representation may lead to a large number of complex piecewise affine functions that will complicate and hamper the solution of the design problem. Even if the linearization point remains fixed, the determination of the controller design (tuning) parameters q j still remains a challenging task. Most of the literature methods46-49 for computing q j focus on controllability rather than economic criteria, rendering these techniques unsuitable for a process design purpose. The methods that use economic criteria35 decompose the problem to a steady state and a dynamic subproblem and determine the tuning parameters based on the steady-state solution. Here, the optimal tuning is performed by adapting an outer approximation based approach50 to our case. Variables vjv, d h , and q j are treated as the complicating variables, thus giving rise to the following sequence of transformations:

yˇ ,vˇ ,η

h /∂yjj)T(yˇ - yˇ j) + (∂φ h /∂vˇ j)T(vˇ - vˇ j) 0eφ h (yˇ j,vˇ j) + (∂φ

s.t.

(M h ) 0 e gj (yˇ ,vˇ j) + (∂gj /∂yˇ j)T(yˇ - yˇ j) + (∂gj /∂vˇ j)T(vˇ - vˇ j),

∀j∈J

0 e gj (yˇ j,vˇ j) + (∂gj /∂yˇ j)T(yˇ - yˇ j) + (∂gj /∂vˇ j)T(vˇ - vˇ j), yˇ ∈ Y h , vˇ ∈ V h

∀j∈I

j

|

where J and I is the collection of all of the feasible and infeasible solutions of problem BNLP(yˆ ), ∀ yˇ ∈ Y h , vˇ ∈ V h. Problem (M h ) is equivalent to (P h ) only if the following assumptions hold: A1. Y h and V h are nonempty compact polyhedral sets. The functions φ h and gj are convex. A2. φ h and gj are continuously differentiable. A3. The normal vectors of the active constraints are linearly independent. Problem (M h ) is infinitely dimensional because it requires all BNLP(yˇ j/yˇ i) to be solved first. The key idea to simplify it is to decompose it into a control design primal problem, where the complicating variables are fixed, and a control design master problem, where the complicating variables are free optimization decisions. This decomposition algorithm is presented in section 4.2 in this paper and guarantees local optimality of the control primal problem for fixed values of the complicating variables. Appendix C: Modeling Aspects of the Parametric Controller The implementation of the parametric controller in (17) involves the selection of the optimal control function according to the operating region where the current system resides. This implies that a logical decision has to be made online on what function will be used based on the current values of the system states. Here, we propose the modeling of this logical rule via steep exponential functions. Thus, problem (17) is written as follows:

φ ) min φ(•,v,•) vv,d

s.t. 0 ) fd(•,v,•), fa(•,v,•), f0(•,v,•) y(t) ) fy(•,v,•), 0 g g(•,v,•)


where Nc

v ) [V, Refl]T ) vB +

∑vc[xr,zbenz,f,yset]gˇ c

(78)

c)1

vc ) Acxr + Cczbenz,f + Dcyset + Bc

(79)

hic ) CRc(xr,zbenz,f,yset) ) CR1icxr(t) + CR3iczbenz,f(t) + CR4icyset + CR2ic (80) Nineqc 1 - tanh(hic × 103) aic, aic ) , gˇ c ) 2 i)1 c ) 1, ..., Nc (81)

∏

vj B ) 0, err ) yset - y, err|t)t0 ) 0

(82)

xr(t) ) (TLT)Tx(t)

(83)

i ) 1, ..., Nineqc

(84)

set T where yset ) [zset benz,D, zbenz,B] . Note that when xr(t) ∈CRc)cˆ, then all hicˆ, where i ) 1, ..., Nineqcˆ, are negative; thus, via (81) all aicˆ are unity; henceforth, gˇ cˆ is also unity. Otherwise, if xr(t) ∉CRcˆ , then at least one hicˆ > 0, so at least one aicˆ ) 0, which yields gˇ cˆ ) 0. Hence, once all gˇ c, where c ) 1, ..., Nc, are substituted back into (78), only the control function pertaining to the region where xr(t) resides contributes with a nonzero coefficient to the value of the control variable. vB is the control bias while equalities (82) are used to establish that the closed-loop system initializes at steady state.

Literature Cited (1) Van Schijndel, J. M. G.; Pistikopoulos, E. N. Towards the Integration of Process Design, Process Control & Process OperabilitysCurrent Status & Future Trends. In Foundations of Computer-Aided Process Design; Malone, M. F., Trainham, J. A., Carnahan, B., Eds.; American Institute of Chemical Engineers: Breckenridge, CO, 2000; Vol. 96, p 99. (2) Luyben, M. L.; Tyreus, B. D.; Luyben, W. L. Analysis of Control Structure for Reactions/Separation/Recycle Processes with Second-Order Reactions. Ind. Eng. Chem. Res. 1996, 35, 758. (3) Yi, C. K.; Luyben, W. L. Design and Control of Coupled Reactor/Column SystemssPart 1. A Binary coupled Reactor/ Rectifier System. Comput. Chem. Eng. 1997, 21, 24. (4) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. J. Process Control 2000, 10, 185. (5) Tousain, R. L.; Meeuse, F. M. Closed loop controllability analysis of process designs: Application to distillation column design. In ESCAPE-11; Gani, R., Jorgensen, S. B., Eds.; 2001; p 799. (6) Seferlis, P.; Grievink, J. Process design and control screening based on economic and static controllability criteria. Comput. Chem. Eng. 2001, 25, 177. (7) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N. A Case Study in Simultaneous Design and Control using Rigorous, MixedInteger Dynamic Optimization Models. Ind. Eng. Chem. Res. 2002, 41, 760. (8) Brengel, D. D.; Seider, W. D. Co-ordinated Design and Control Optimization of Non Linear Processes. Comput. Chem. Eng. 1992, 16, 861. (9) Loeblein, C.; Perkins, J. D. Structural Design for On-Line Process Optimization: I. Dynamic Economics of MPC. AIChE J. 1999, 45, 1018.

(10) Swartz, C. L. E.; Perkins, J. D.; Pistikopoulos, E. N. Incorporation of Controllability in Process Design Through Controller Parameterization. Process Control Instrum. 2000, 49-54. (11) Zhu, G. Y.; Henson, M. A. Model Predictive Control of Interconnected Linear and Nonlinear Processes. Ind. Eng. Chem. Res. 2002, 41, 801. (12) Ricker, N. L. Decentralized control of the Tennessee Eastman Challenge Process. J. Process Control 1996, 6, 205. (13) Pistikopoulos, E. N.; Bozinis, N. A.; Dua, V. POP, a MATLAB implementation of parametric programming algorithms; Technical report; Centre for Process Systems Engineering, Imperial College: London, 1999-2002. (14) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal design of Dynamic Systems under Uncertainty. AIChE J. 1996, 42, 2251. (15) Pistikopoulos, E. N.; Grossmann, I. E. Optimal Retrofit Design for Improving Process Flexibility in Linear Systems. Comput. Chem. Eng. 1988, 12, 719. (16) Grossmann, I. E.; Halemane, K. P.; Swaney, R. E. Optimization Strategies for Flexible Chemical Processes. Comput. Chem. Eng. 1983, 7, 439. (17) Bansal, V. Analysis, Design and Control Optimization of Process Systems under Uncertainty. Ph.D. Dissertation, Imperial College of Science, Technology and Medicine, London, 2000. (18) Bansal, V.; Sakizlis, R.; Ross, V.; Perkins, J. D.; Pistikopoulos, E. N. New Algorithms for Mixed Integer Dynamic Optimization. Comput. Chem. Eng. 2003, 27, 647. (19) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints. Ind. Eng. Chem. Res. 1994, 33, 2123. (20) Vassiliadis, V. Computational Solution of Dynamic Optimization Problems with General Differential-Algebraic Constraints. Ph.D. Dissertation, Imperial College of Science, Technology and Medicine, London, 1993. (21) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. M. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789. (22) Scokaert, P. O. M.; Rawlings, J. B. Constrained Linear Quadratic Regulation. IEEE Trans. Autom. Control 1998, 43, 1163. (23) Rawlings, J. B.; Muske, K. R. The Stability of Constrained Receding Horizon Control. IEEE Trans. Autom. Control 1993, 38, 1512. (24) Chmielewski, D.; Manousiouthakis, V. On constrained infinite-time linear quadratic optimal control. Syst. Control Lett. 1996, 29, 121. (25) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control; Wiley and Sons: New York, 1989. (26) Pistikopoulos, E. N.; Dua, V.; Bozinis, N. A.; Bemporad, A.; Morari, M. On-line optimization via off-line parametric optimization tools. Comput. Chem. Eng. 2002, 26, 175. (27) Muske, K. R.; Edgar, T. F. Nonlinear State Estimation. Nonlinear Process Control; Prentice Hall: Englewood Cliffs, NJ, 1997; p 233. (28) Rao, C. V.; Rawlings, J. B. Constrained process monitoring: Moving-horizon approach. AIChE J. 2001, 48, 97. (29) Schweiger, C. A.; Floudas, C. A. Interaction of Design and Control: optimization with Dynamic models. In Optimal Control: Theory, Algorithms and Applications; Hager, W. W., Pardalos, P. M., Eds.; 1997; p 388. (30) Kwakernaak, H.; Sivan, R. Linear Optimal Control Systems; John Wiley and Sons Inc.: New York, 1972. (31) Gumus, Z. H.; Floudas, C. A. Global Optimization of Nonlinear Bilevel Programming Problems. J. Global Optim. 2001, 20, 1. (32) Pistikopoulos, E. N.; Sakizlis, V. Simultaneous Design and Control Optimization under Uncertainty in Reaction/Separation Systems. In CPC-VI Proceedings; Rawlings, J., Ogunnaike, T., Eaton, J., Eds.; AIChE Symposium Series; AIChE: New York, 2002; Vol. 98, p 223. (33) Moore, B. C. Principal Component Analysis in Linear Systems, Controllability, Observability, and Model Reduction. IEEE Trans. Autom. Control 1981, 26, 17. (34) Jaimoukha, I. M.; Kasenally, E. M.; Limebeer, D. J. N. Numerical solution of large scale Lyapunov equations using Krylov subspace methods. 31st IEEE Conf. Decision Control 1992, 2, 1927. (35) Kookos, I. K.; Perkins, J. D. An algorithmic method for the selection of multivariable process control structures. J. Process Control 2002, 12, 85.

Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4563 (36) Newell, R. B.; Lee, P. L. Applied Process ControlsA Case Study; Prentice-Hall: Sydney, Australia, 1989. (37) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw-Hill International Edition: New York, 1988. (38) Chemical Engineering; McGraw-Hill: New York, 2001. (39) Meadows, E. S.; Rawlings, J. B. Model Predictive Control. Nonlinear Process Control; Prentice Hall: Englewood Cliffs, NJ, 1997; p 233. (40) Narraway, L.; Perkins, J. D. Selection of Process control Structure based on Linear Dynamic Economics. Ind. Eng. Chem. Res. 1993, 32, 2681. (41) Keerthi, S. S.; Gilbert, E. G. Optimal infinite - horizon feedback law for a general class of constrained discrete-time systems: Stability and moving horizon approximations. J. Optim. Theory Appl. 1988, 54, 265. (42) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Robust stability considerations in optimal design of dynamic systems under uncertainty. J. Process Control 1997, 7, 371. (43) Kookos, I. K.; Perkins, J. D. An algorithm for simultaneous process design and control. Ind. Eng. Chem. Res. 2001, 40, 4079. (44) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The properties of gases and liquid; McGraw-Hill: New York, 1987.

(45) Dua, V.; Papalexandri, K. P.; Pistikopoulos, E. N. Global Optimization Issues in Multiparametric Continuous and MixedInteger Optimization Problems. J. Global Optim. 2002, in press. (46) Ge, M.; Chiu, M.; Wang, Q. Robust PID controller design via LMI approach. J. Process Control 2002, 12, 3. (47) Tan, K. K.; Lee, T. H.; Leu, F. M. Optimal Smith-Predictor Design Based on a GPC Approach. Ind. Eng. Chem. Res. 2002, 41, 1242. (48) Al-Ghazzawi, A.; Ali, E.; Nouh, A.; Zafiriou, E. On-line tuning strategy for model predictive controllers. J. Process Control 2001, 11, 265. (49) Shridhar, R.; Cooper, D. J. A Tuning Strategy for Unconstrained Multivariable Model Predictive Control. Ind. Eng. Chem. Res. 1998, 37, 4003. (50) Fletcher, R.; Leyffer, S. Solving mixed integer nonlinear programs by outer approximation. Math. Program. 1994, 66, 327.

Received for review November 18, 2002 Revised manuscript received March 4, 2003 Accepted March 6, 2003 IE0209273

Parametric Controllers in Simultaneous Process and Control Design

Recommend Documents