Simultaneous Design and Control: A New Approach and Comparisons

2822

Ind. Eng. Chem. Res. 2010, 49, 2822–2833

Simultaneous Design and Control: A New Approach and Comparisons with Existing Methodologies Luis A. Ricardez-Sandoval,* Hector M. Budman, and Peter L. Douglas Department of Chemical Engineering, UniVersity of Waterloo, Waterloo, Ontario, N2L 3G1, Canada

In this work, a new methodology to conduct simultaneous design and control of large-scale systems is presented. The proposed method uses a structured singular value norm calculation (µ) to estimate the worst-case disturbance profile. This profile is then used to simulate the closed-loop nonlinear dynamic process model for obtaining the worst-case output variability and to test the process feasibility constraints. Thus, the proposed method referred to as the hybrid worst-case approach (HWA) methodology combines the analytical µ calculation of the worst-case disturbance and dynamic simulations using the mechanistic closed-loop process model to calculate variability. To test the proposed HWA method the integration of design and control of the reactor section of the Tennessee Eastman process was analyzed. This case study was also used to compare the proposed HWA method to other previously reported methodologies. Although the results obtained by the present HWA method for the case study are slightly conservative, the computational demand required by the present method is found to be 1 order of magnitude smaller than that required by a dynamic optimization-based methodology. 1. Introduction Integration of design and control refers to the design of a chemical process by considering both the steady state and the dynamic closed-loop operation of the process under the effect of external disturbances and process uncertainties. This concept requires the minimization of an objective function related to the process economics subject to process feasibility constraints. Due to the mathematical complexity that arises when both the steady-state and dynamic decisions are simultaneously considered for optimization, a general unified framework that completely addresses this issue is not currently available. Consequently, several methodologies have been proposed to address different partial aspects of this problem. The first methodologies developed in this field measured the variability in the process using controllability metrics such as relative gain array, condition number, disturbance condition number, or integral squared error (ISE).1-5 Some of the drawbacks in the use of these metrics are that they are applicable to steady-state or linear dynamic process models and that these metrics cannot be assigned an explicit economic value. To circumvent the issues associated with the controllability indices’-based methods, a dynamic optimization-based approach was proposed to tackle this problem.6-15 In this approach, the closed-loop nonlinear dynamic process model is used to perform a dynamic flexibility analysis. This analysis produces the critical realizations in the disturbance and process parameter uncertainties that result in the worst variability in variables that are most relevant to the objective function and to the constraints. The disturbance and process parameter uncertainty profiles thus obtained are then used to select the best design and controller parameters that minimize an expected value-based cost function. Although these methods estimate the critical disturbance profile in the analysis, the computational demands required by these techniques are significant even when a small number of process units are considered in the design. Thus, the application of these methods to large-scale systems is constrained by the computational resources available. To alleviate the computational demands associated with the dynamic optimization-based methodologies, robust approach-based methodologies have been recently proposed.16-20 The key idea in this approach is to represent * To whom correspondence should be addressed. Tel.: (519) 888 4567 (ext 38667). Fax: (519) 746 4979. E-mail: [email protected].

the process nonlinear dynamic model as uncertain models that can be used to calculate bounds on the variables that are involved in the objective function and the constraints of the problem under consideration. Although the use of this robust approach may be conservative and result in suboptimal designs, it is attractive from the computational point of view. Also, apart from Ricardez-Sandoval et al.,19,20 the rest of the methods applied to the integration of the control and design problem require a priori specification of a time-varying function for the disturbance variable. A comprehensive review of methodologies that were used for integration of design and control can be found elsewhere.12,21,22 Due to the high computational demand of previously reported methods, most of the examples involved a small number of process units. The methods that have been applied to largescale systems2,23 required the a priori specification of a disturbance function. Thus, these methods did not search for the worst-case scenario, and therefore, the resulting final design and control scheme obtained by these methods cannot guarantee feasibility in the face of other realizations in the disturbance variables. The main limitation usually faced when dealing with largescale systems is the computational burden involved in the search for a worst-case scenario, i.e., the search of the particular timedependent disturbance profile that drives the system’s output to its largest variability around a nominal operating state. This paper presents a simultaneous design and control methodology that applies robust control tools to estimate the critical realizations in the disturbance variable that produce the largest output variable, i.e., the worst-case scenario. The worstcase disturbance is then simulated with the process nonlinear dynamic model to evaluate the process variability in the closed loop and test the process constraints. Also, a stability test24 is included in the analysis to guarantee that the final design obtained by the proposed method is stable. In a previous study recently conducted by the authors,24 the bounds on the output variability and constraints-related variables were calculated by using vector norms. The key disadvantage with that earlier approach is the conservatism resulting from the use of analytical bounds. Instead, in the present work, only the critical realizations in the disturbance variables are found by an analytical technique whereas the calculation of the

10.1021/ie9010707  2010 American Chemical Society Published on Web 02/16/2010

Ind. Eng. Chem. Res., Vol. 49, No. 6, 2010

maximal variability in outputs and constraint related variables’ is done by simulating the full nonlinear dynamic model with the calculated critical disturbance realizations. Accordingly, since the current approach is comprised of an analytical calculation of the worst disturbance combined with a simulationbased calculation of the responses to that worst disturbance, it will be referred heretofore as a hybrid worst-case approach (HWA). In comparison with the previous bound-based worstcase approach (BWA) that uses analytical bounds of the output variables, the proposed hybrid method is expected to reduce the conservatism in the final design at the expense of additional computational time in the calculations. To verify this point, the results obtained by the HWA and those obtained by the BWA for a case study are compared. Moreover, the current proposed method is also compared to a dynamic optimization-based strategy and an optimal tuning/set-point selection problem. These comparisons are made on the basis of both optimal solutions and the computational demand to obtain these solutions. The case study involved the simultaneous design and control of the reactor section of the Tennessee Eastman (TE) process.25 A closed-loop nonlinear dynamic model of this process is available for simulations, and the process is representative of large-scale chemical systems since it involves five major process units with a recycle stream. Although studies have shown that the current capacities of the process units involved in this plant limit the closed-loop performance of this system,26,27 only one study, recently presented by the authors,24 has formally addressed the optimal integration of design and control of this process. This paper is organized as follows: Section 2 presents the mathematical description of the hybrid worst-case methodology proposed in this work. Section 3 briefly describes the Tennessee Eastman process and reviews the control structure selected for this work. Section 4 presents the application of the methodology presented in section 2 to the reactor section of the TE process. The simultaneous design and control of the TE process using a bound-based worst-case methodology and a dynamic optimization-based method are presented in this section. Comparisons in terms of the resulting optimal solutions as well as the computational demands required by each method are also presented in this section. Section 5 presents a second scenario proposed to solve the integration of design and control of the reactor section of the TE process. This second scenario was solved using the three methods used in section 4. Also, for comparison purposes, an additional problem was solved where only the tuning parameters and set points were optimized whereas the size of the reactor was fixed to its current value. This last problem will be referred heretofore as the optimal tuning/set-point selection problem. A discussion between the optimal tuning/set-point problem and the design and control approach is also presented in this section. Concluding remarks are presented in section 6. 2. Simultaneous Design and Control Methodology This section introduces first the calculation method used to estimate the critical realization in the disturbance variables. This is followed by the presentation of the proposed hybrid simultaneous design and control methodology. 2.1. Disturbance Description. The set of disturbances d affecting the process is defined in eq 1. Each element in the set d, e.g., dq, represents a particular disturbance in the process and is defined as a vector of length N + 1, where N is the process settling time corresponding to the slowest disturbance. As shown in eq 1, the vector dq contains the disturbance values of the disturbance qth at different time intervals, j.

2823

d ) {d1, ..., dq, ... dm} dq ) {dq(j)|dlower (j) e dq(j) e dlower (j)}, for j ) 0, 1, 2, ..., N q q (1) The present work assumes that each disturbance vector is magnitude bounded by lower and upper limits at each time interval j, i.e., dqlower(j) and dqupper(j), respectively. Thus, the set of the maximal deviation in the disturbance variables (δd) is defined as follows δd ) {δd1, ..., δdq, ... δdm}

{

|}

|

dlower (j) - dupper (j) q q δd (j)|δd (j) e , q q δdq ) 2 for j ) 0, 1, 2, ..., N

(2)

where the vector δdq contains the maximal deviation of the qth disturbances at any time interval, i.e., δdq is of length N + 1. 2.2. Uncertain FIR Model. The present methodology assumes that a mechanistic nonlinear dynamic model is available for simulations. This mechanistic nonlinear model is generally based on first principles. Then, the key idea of the proposed hybrid method is that for the purpose of estimating the worstcase disturbance, the closed-loop process behavior is approximated by an uncertain finite impulse response (FIR) model where each of the impulse response coefficients in the model is bounded between lower and upper limits. Since the FIR models used for estimating the worst-case disturbance describes the closed-loop behavior of the system, the inputs into this model are the disturbances affecting the process whereas the outputs are process variables that can be used to estimate the worstcase scenario or variables that need to be kept within constraints during closed-loop operation. Thus, both manipulated and controlled variables are considered as outputs of closed-loop FIR models used in the formulation of the optimization problem. The uncertainty in the model parameters is used to capture the differences between the actual process nonlinear model and the FIR model due to nonlinearity and truncation errors. The FIR models are identified around a nominal operating point, specified by the design and controller parameters, which are also optimization variables (η) in the integration of design and control methodology. Thus, the uncertain models must be identified for each new set of decision variables tested by the optimization algorithm. The uncertain FIR model that describes the dynamics between the qth disturbance and a particular process variable yq at a time interval, j, is defined as follows j

yq(j) )

∑ (h

kq

+ δhkq)δdq(j - k);

for

j ) 0, 1, 2, ..., N

k)0

(3) where yq(j) represents the output response at time interval j to a unit impulse in the qth disturbance. The term hkq represents the nominal value in the impulse response coefficients due to changes in the qth disturbance at time interval k. Likewise, δhkq is the uncertainty associated with the impulse response coefficients of the qth disturbance at kth time interval. The scalars yq(j), hkq, and δhkq are lumped in vectors yq, hq, and δhq of lengths N + 1, respectively. Since more than one disturbance may be affecting the process at any time interval, the output response to m disturbances can be estimated as follows: m

y(j) )

j

∑ ∑ (h

kq

+ δhkq)δdq(j - k);

for

j ) 0, 1, 2, ..., N

q)1 k)0

(4)

2824


where y(j) represents the output response of the system to a unit pulse in the m disturbances at time interval j. Therefore, the process output response to changes in the disturbance set is modeled by a set of uncertain FIR models. The linear impulse response coefficients hkq are obtained from identification of the closed-loop process model using the leastsquares method.28 This identification is based on closed-loop simulations of the original dynamic nonlinear model under closed-loop control. The uncertainty associated with each of these linear model parameters δhkq is assumed to be the variance of the linear model parameters’ estimates. The detailed procedure used to obtain the uncertain FIR model shown in eq 3 has been reported in a previous study.24 2.3. Worst-Case Disturbance Analysis. In the current work, a robust variability measure was used to search for the worstcase disturbance entering the process. The challenge is to calculate this worst-case disturbance and its effect on the variables that have a direct impact on the objective function and the process constraints. The robust measure is based on a mixed-structured singular value (µ) analysis, which is commonly used in the design of robust controllers.29,30 This robust control analysis tool has been used in previous integration of design and control studies as a screening tool to select the most promissory control structures in multiloop control in the presence of model uncertainty.31,32 On the basis of this analytical tool, the search for the largest variability in a process variable y at any time interval j with respect to a nominal operating condition, η, can be formulated as follows:

(5)

max |y| d

Although the solution of this problem with the full dynamic model is highly computationally demanding, an effective calculation can be done instead with the uncertain FIR model defined above that approximates the behavior of the actual system as given by the mechanistic nonlinear model. On the basis of the FIR model, a bound for y can be reformulated in terms of a mixed structured singular value (µ) analysis as follows33,34

the worst output variability. Although the use of simulations with the full dynamic model increases the computational time as compared to the use of the bounds-based approach (γ), the reduction in conservatism was found to be very significant as shown later in the case studies. The resulting worst-case profile in the disturbance set (d*) can be estimated from the optimal value of γ (γ*) as follows d* ) {d* 1, ..., d* q, ...,d* m} *(j)|d *(j) ) {d ) γ*δdq(j)R*j,q} for j ) 0, 1, 2, ..., N d* q q q

(7) where the term R*j,q represents a normalized value of the qth disturbance set at time interval j such that when combined with the rest of the normalized values for the qth disturbance returns the worst-case normalized temporal profile for the qth disturbance. The complete set of R’s are lumped in a vector ∆r1, which represents the first m(N + 1) diagonal elements in the perturbation matrix ∆ (see eq A5 and A6 in Appendix A). Each of the elements in ∆ are estimated from the µ-analysis calculation29 for which the optimal γ (γ*) and the worst-case normalized values for the disturbances (∆r1*) are obtained. Each of the vectors considered in the set d* is of length N + 1. Once the worst-case profile for each disturbance considered in the analysis has been calculated, then the system response to this worst-case disturbance is obtained by simulating the process nonlinear dynamic model with the set d* as an input. The results of this simulation are used to obtain the actual largest variability in a process variable y and other variables that have to be maintained within constraints in the system to be designed. 2.4. Hybrid Worst-Case Approach (HWA) Formulation. The hybrid formulation proposed in the current work to perform the simultaneous design and control of large-scale systems is as follows min CC(ξ, g¯, Γ) + OP(ξ, g¯, Γ) + VC(ξ, g¯, λ, Γ)

η)[ξ,g¯,λ]

s.t. Γ ) Γ¯ + γ*ξ g¯ ( γ*g e G

(6)

ηl e η e ηu γ*ξ ) [γ*ξ,1, γ*ξ,2, ..., γ*ξ,L]T

where the parameter γ denotes an infinite-time horizon bound on the worst-case output variability problem specified in eq 5. The matrices M and ∆ represent the corresponding interconnection and perturbation matrices of the mixed µ analysis. The structures of these matrices are presented in Appendix A. As shown in this Appendix, the elements of matrix M are a function of the model parameters of the uncertain FIR model, hkq and δhkq, the maximal deviation in the disturbance set, δdq(j), and the parameter γ, whose value is calculated at each iteration step by the optimization algorithm used to solve the problem shown in eq 6. The problem posed in eq 6 is also known in the literature as the skewed-µ calculation30 and can be solved with off-theshelf software. The solution of this problem returns the following. (i) The value in the parameter γ that denotes an analytical bound on the largest process output variability for the variable y around a nominal operating condition. (ii) The critical realizations in the disturbance set d that generates the largest variability for the variable y around a nominal operating condition, η. In previous work done by the authors24 the parameter γ has been used as an estimate of the worst output variability within the optimization problem. This approach was found in the previous work to lead to conservative solutions. Instead, in the current work, the critical realization in d is simulated with the full nonlinear dynamic model of the process to calculate

γ*g ) [γ*g,1, γ*g,2, ..., γ*g,Q]T

max |y| g γ S max γ d

µ∆(M)gγ

(8)

Although the uncertain FIR and nonlinear mechanistic closedloop models do not explicitly appear in eq 8, they are implicitly used at each calculation step since the calculation of the γ’s, eq 6, and the d*’s, eq 7, is done with the uncertain FIR model parameters, which are obtained from simulations of the nonlinear mechanistic closed-loop model. The cost function in this problem is defined as the addition of the capital (CC), operating (OP), and variability costs (VC) of the process to be designed. These costs are functions of the design parameters (ξ), the controller parameters (λ), and the nominal values in those process variables that must remain within bounds at any time (gj). These variables, lumped in a vector η, represent the optimization variables in problem 8. The capital and operating costs represent the steady-state costs of the process, whereas the process variability cost (VC) depends on the closed-loop dynamic performance. The variability cost is estimated with a variability function, Γ, which is process specific and must be defined in terms of the goals to attain by the design. As shown j and in eq 8, the process variability function is a function of Γ j γξ*. Each of the elements in Γ represents a steady-state value of a process variable considered in the calculation of Γ, and this value is estimated using the steady-state process model and the current values in η. Similarly, each of the parameters considered


2825

Figure 1. Decentralized control configuration.

in γ*ξ represent the maximum expected variability for a process variable considered in the formulation of Γ. Each γ*ξ is calculated from the µ-analysis formulation presented in eq 6. For example, if it is desired to design the volume of a tank in the presence of feed-flow disturbances, the process output variability function Γ will be equal to the steady-state volume plus the maximum variability of the volume with respect to the steady-state value j , and γξ* is due to flow disturbances. The length of vectors Γ, Γ L, where L denotes the number of process variables considered in the formulation of Γ. Thus, L is problem specific. Once Γ has been calculated, a variability cost function (VC in eq 8) that accounts for loss profit due to closed-loop performance can be estimated by assigning a cost to each of the elements in Γ. As shown in problem 8, the worst-case disturbance analysis is also used in the proposed HWA method to test the compliance with the process constraints around a nominal steady-state operating point. The term gj represents the nominal value in a particular process variable that is required to be within prespecified input limits given by G. The worst-case variability for g around gj is estimated from the worst-case disturbance analysis. denote the largest variability for a particular The parameters γ*’s g process variable that must be kept within constrains. Each γ*g is calculated using the µ-analysis formulation presented in eq 6. The length of vectors gj, G, and γg* is Q, where Q denotes to the number of process constraints considered in the optimization formulation. The two different subscripts ξ and g that appear in worst-case scenario formulation in eq 8 are used to denote that different critical realizations in d can be obtained for the worst-case process variability for different variables associated to either the objective function or to the constraints. To guarantee that the resulting design is stable in the face of the worst-case scenario, a stability test, based on the uncertain FIR model and the parameter γ obtained by the solution of problem 6, is included as a constraint within the HWA optimization.24 The HWA method assumes that the process flow sheet and control structure have been defined a priori. Although integer decisions can be added into the optimization problem to obtain an optimal flow sheet and controller structure, the computational

demands for this problem may be significant. The addition of these topics within the proposed simultaneous design and control methodology is under current analysis by the authors. Since the optimization problem presented in eq 8 requires at each iteration step dynamic simulations to estimate uncertain FIR models and the worst-case process variability, the problem shown in eq 8 can be considered as a dynamic optimization problem. However, the proposed formulation (HWA) estimates rapidly the worstcase disturbance by an analytical µ-based technique rather than finding this worst disturbance from a computationally expensive search with dynamic optimization. This issue is further discussed in section 4.2, where a comparison of the computational times required by the different methodologies is presented. 3. Tennessee Eastman Process The Tennessee Eastman (TE) process25 is used in this work to test the proposed methodology and for its comparison to other approaches. This process produces two products, G and H, and a byproduct, F, using four reactants, A, C, D, and E, and an inert component, B. Figure 1 shows the TE process flow sheet. Reactants A, C, D, and E are fed to the reactor unit together with a mixture of gaseous components that comes from a recycle stream. The components react inside the reactor to form the liquid products G and H. The products and unreacted components leave the reactor as a vapor. The reactor’s outlet stream is passed through a partial condenser unit. This two-phase stream is separated within a flash unit. The noncondensed components are recycled back to the reactor’s feed stream through a centrifugal compressor. Also, a purge stream is used to remove the excess of the inert component (B) and the byproduct (F) from the vapor stream of the flash unit. The liquid collected in the separator is pumped to a stripper column for further refinement. In this unit, a mixture of reactants A and C and inert B, which are fed to the process through the stripper base, is used as a solvent to separate the products from the reactants. The vapor stream from the stripper is mixed with the recycle stream. The modes of operation of this chemical plant are

2826


Table 1. Control Configuration Characteristics35 controllers’ parameters loop

controlled variable

manipulated variable

set point (base-case value)

Kc

τI (min)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18a 19a

A feed rate (stream 1) D feed rate (stream 2) E feed rate (stream 3) C feed rate (stream 4) Purge rate (stream 9) Sep. Liq. rate (str. 10) Strip. Liq. rate (str. 11) production rate stripper liquid level separator liquid level reactor liquid level reactor pressure mol % G in stream 11 yA yAC reactor temperature separator temperature max react. press. reactor level override

A feed flow (stream 1) D feed flow (stream2) E feed flow (stream 3) C feed flow (stream 4) Purge valve (stream 9) Sep. pot Liq. flow (stream 10) Strip. Liq. flow (stream 11) production index (Fp) ratio in loop 7 (r7) ratio in loop 6 (r6) set point of loop 17 (l17) ratio in loop 5 (r5) ratio in loop 2 (r2) and 3 (r3) ratio in loop 1 (r1) ratio in loop 4 (r4) reactor coolant valve condenser coolant valve production index (Fp) recycle valve

r1 × Fp r2 × Fp r3 × Fp r4 × Fp r5 × Fp r6 × Fp r7 × Fp production set point (22.89 m3/h) strip. level set point (50%) sep. level set point (50%) reac. level set point (65%) reactor pressure (2800 kPa) mol % G set point (53.83%) yA set point (63.1373%) yAC set point (51%) reac. temp. set point (122.9 °C) l17 production set point reac. level set point

0.01 1.6 × 10-6 1.8 × 10-6 0.003 0.01 4 × 10-4 4 × 10-4 3.2 -2 × 10-4 -1 × 10-3 0.8 -1 × 10-4 -0.4 2 × 10-4 3 × 10-4 -8.0 -4 2.0 × 10-6 1.0 × 10-6

0.001 0.001 0.001 0.001 0.001 0.001 0.001 120 200 200 60 20 100 60 120 7.5 15 0.001 1.0 × 10-5

a

Override loops not shown in Figure 1.

specified according to the desired products’ mix and mass flow rates in the product’s stream. In this work, the TE’s base-case mode of operation, given by a product mass ratio of 50/50 (G/ H) with a production rate of 7038 kg/h for each product, was the only operating mode considered in the analysis. Also, for safety reasons, six process operational constraints have been specified for the reactor, the flash, and the stripper units. Furthermore, variability constraints in the products stream, e.g., (10% in the flow rate and (5% in the mol G’s composition, have been also specified for this process.25 A nonlinear mechanistic model coded in Fortran that mimics the open-loop behavior of the real chemical plant is available. The process model equations in the TE code are formulated in standard nonlinear state space form as follows x˙(t) ) f(x, g, d, t) z(t) ) l(x, t)

(9)

where x represents the process states variables, z is the measured outputs vector, and u is the available manipulated variables vector. The vectors x, z, and u are of lengths 50, 41, and 12, respectively. Twenty possible process disturbances (d) have been proposed for the TE process model, and a function describing the operating costs has been also reported for the basic operating mode.25 The TE process is an open-loop unstable process and exhibits a complex dynamic behavior with a high degree of interaction between the process variables. Also, the plant includes a recycle stream that introduces an inherent positive feedback leading to further control challenges. To apply the proposed design and control method, it is necessary to specify a suitable control structure. In the present work, a previously defined decentralized control strategy35 was used to control this process consisting of 17 PI controllers. Figure 1 shows the TE process flow sheet and decentralized control configuration. This control configuration was designed based on process heuristics. This particular control structure was chosen because it can regulate the process during transitions between its different operating modes and can reject almost all the disturbances specified for this process. Table 1 shows the control configuration characteristics. As shown, there are nine variables whose corresponding set points can be potentially defined as design variables (ξ) and whose values must remain within input limits (gj). To account for changes in the design of the TE process, the capacities of the process units, e.g., the reactor, are considered to be part of the design variables for this problem. Similarly, the controllers’ tuning parameters (λ) can be found from optimization.

Therefore, a total of 46 decision variables (η ) [ξ, g¯, λ]) can be potentially specified for the integration of design and control of the TE process. The objective in the present problem is to find the values in the design and control parameters of the process (η) that minimize a cost function similar to that specified in problem 8 such that the process operational constraints are met in the face of the worst-case disturbance (d*). The following sections present the application of the HWA methodology to the TE process. 4. Scenario A: Reactor Section of the TE ProcesssRandom Variation in A and C Feed Composition The first scenario considers the simultaneous design and control of the reactor section of the TE process. Thus, the only process design variable that was optimized is the reactor’s capacity, whereas the other process units’ capacities remained fixed at its current values.25 This could correspond to a hypothetical situation where it is desired to retrofit the process by upgrading the reactor section without altering the rest of the plant. Also, the PI controllers’ tuning parameters that are involved in the operation of this unit were optimized. The controllers are listed in Table 1 as loop 5, 11, 12, and 16. Accordingly, the set points for loops 8, 11, 12, and 16 were considered only in this analysis (see Table 1 and Figure 1). Thus, the set of decision variables for variables (η) is composed of a total of 13 process variables. Although Downs and Vogel25 considered a set of 20 disturbances in their original work, for simplicity, in this case study a random change in A and C feed composition in stream 4 (see Figure 1) was selected as the disturbance for this problem. Previous studies have acknowledged that this perturbation is one of the most challenging disturbances that can be considered for the system.35,36 The set of disturbances specified for this case study are as follows d ) {dA, dC} dA ) {dA, 0.475 e dA e 0.495} dC ) {dC, 0.480 e dC e 0.540}

(10)

where dA and dC represent a change in stream 4’s A and C feed composition, respectively. Each composition is magnitude bounded by the lower and upper limit values specified in eq 10. Although the maximal change in each of these components seems to be small, simulations have shown that the effect of


these changes on the plant’s dynamic behavior is significant. For this scenario, the composition of reactant B in this stream was kept constant and equal to 0.005. The cost function used for this analysis is as follows: CFTE ) CC + 8760(OP + VC) CC ) 1359.1D2.1(19.82 - 12.55 ln P + 2.36(ln P)2) OP ) (purge costs)(purge rate) + (product stream costs) × (product rate) + (compressor costs) × (compressor work) + (steam costs)(steam rate) VC ) cpGΓG + cpHΓH (11) The annualized capital cost function (CC) assumes that the TE plant was built in 1992 with a 20% rate of return. Also, the reactor was assumed to be a pressured vessel made of carbon steel with a length/diameter ratio of 4. The cost of the reactor is obtained from a correlation reported in the literature37 that is a function of the vessel’s diameter (D) and working pressure (P), respectively. The OP cost function used in this study is the same as that reported by Downs and Vogel.25 Since the production goal is to produce products G and H with a desired mix and at a given mass flow rate, the present study measured the process variability as the difference between the largest deviation in the products’ mass flow rates, estimated from simulations of the worst-case disturbance set (d*), and the products’ mass flow rate production targets (G* and H*). Thus, the process worst-case variability functions for each product are defined as follows j + γ*G - G* ΓG ) G j ΓH ) H + γ*H - H*

(12)

¯ and H ¯ represent the set points for each variable and where G ¯ and γ*G and γ*H are the maximum expected variability around G ¯ , respectively. To solve for γ*G and γ*, H H uncertain FIR models like (3) must be identified from the disturbances dA and dC to the products’ mass flow rates in streams 11, G, and H. The solution of these maximization problems gives the analytical bound γ* and the elements in vector ∆r1* (r*’s), used in equation (7) to estimate the worst-case disturbance profiles in * and dC* are then used dA and dC, i.e. dA* and d*. c The profiles dA to quantify the maximum variability in G and H (γ*G and γ*) H from simulations of the mechanistic nonlinear model of the TE process under closed-loop control. Downs and Vogel25 quantify the production costs at 0.22 $/kg of product. Accordingly, the costs of producing G and H, defined as cpG and cpH in eq 11, are assumed to be equal to 0.22 $/kg multiplied by the calculated maximum variability in the production rates of G and H, respectively. The constant term in the cost function CFTE in eq 11 is used to convert OP and VC to annualized costs. According to the TE process description, the process must operate within operational constraints. To account for these restrictions in the current design and control formulation, constraints like the ones shown in eq 8 were considered in the analysis. These constraints are required for safe operation of the reactor. Also, process variability constraints in the products’ mass flow rates in the product stream and the composition of product G in the products’ stream (stream 11 in Figure 1) were also considered in the present study. Moreover, the present HWA formulation added steady-state constraints in the design and control formulation to guarantee that the production target goals specified for the base-case mode of operation were satisfied at the steady state. On the basis of the above, the formulation proposed to perform the simultaneous design and control of the reactor section of the TE process is shown in eq 13. The terms xmG

2827

and xmH are the mass fractions of products G and H in the products stream. The values for xmG and xmH are calculated from the actual values in the decision variables set (η) and the steady-state solution of the mass and energy balances equations specified for the TE process model. The term jxG represents the set point in the G’s mole fraction in the products stream. jxG was not considered for optimization, and therefore, its value was fixed to the value shown in Table 1. The constraint in xG was added in the analysis to ensure that the variability in j r, this variable does not exceed (5 mol %.25 The variables P j r represent the reactor pressure, liquid level, and j r, and T L temperature set points, respectively. The inequalities specified in eq 13 for each of these process variables guarantee that the reactor will operate within the feasible region specified for this process unit in the presence of disturbances dA and dC. Also, j r, and T j r are also included in the vector η, and consequently, j r, L P they are obtained from the solution of problem eq 13. minCFTE η

s.t. j + γ*G - G* ΓG ) G j ΓH ) H + γ*H - H* xmG - 1.02xmH e 0 0.98xmH - xmG e 0 0.98 - xmG - xmH e 0 j + γ*G e 0 G* - G j H* - H + γ*H e 0 - 0.05xjG + γ*xG e 0 j Pr + γ*P - 2895 e 0

(13)

r

j r + γ*L - 100 e 0 L r j r + γ*L e 0 30 - L r

j r + γ*T - 150 e 0 T r η l e η e ηu For each new set of decision variables η selected by the optimization algorithm, uncertain FIR models must be identified from disturbances considered in the analysis, dA and dC, to process variables G, H, xG, Pr, Lr, and Tr. The disturbance specifications and model parameters of each uncertain model are used in the µ analysis (eq 6) to calculate the optimal analytical bound γ* and the elements in the vector ∆r1* (R*’s) for each process constraint and at each function evaluation. These variables are used in eq 7 to calculate the worst-case profile in the A and C’s feed composition, d*A and d*. C Therefore, there will be a critical realization dA* and dC* for each constraint considered in the formulation. The profiles dA* and dC* obtained for each process variable are then used to simulate the closedloop behavior of the TE process using the actual nonlinear mechanistic model. In that way the actual maximum variability in that process variable for the worst dA* and dC* were obtained. The worst-case variability for each process variable is used to test the inequality constraints in eq 13. Moreover, a robust stability test24 was performed at each function evaluation to guarantee that the final solution of this problem is asymptotically stable. The formulation presented in eq 13 was implemented in MATLAB and solved using the fmincon built-in function as the optimization algorithm.38 The initial guesses for the set points and controller parameters considered in the decision variable set η were taken from Table 1. The initial guess for the reactor’s capacity was set to the original value,25 i.e., 1300 ft3. The lower and upper bounds for each of the decision variables were found from preliminary simulations of the closedloop TE process model. Table 2 summarizes the decision

2828


Table 2. Design Parameters: Scenario A design and control parameters (η)

HWA

BWA

dynamic optimization approach

reactor’s design capacity (ft3) reactor’s pressure set point (KPa) reactor’s liquid level set point (%) reactor’s temperature set point (°C) production set point (m3/h) Kc, purge valve (loop 5) τI, purge valve (loop 5) Kc, reactor liquid level(loop 11) τI, reactor liquid level (loop 11) Kc, reactor pressure (loop 12) τI, reactor pressure (loop 12) Kc, reactor temperature (loop 16) τI, reactor temperature (loop 16)

3000 2888.0 38.0 117.00 23.034 0.02 6.8e-4 1.28 60.34 -9.9e-4 15.9 -20.0 14.71

2972.40 2600.0 63.2 117.00 23.11 0.10 1.0e-2 1.50 70.62 -5.0e-5 60.0 -20.0 11.58

2986.20 2870.0 38.0 117 23.031 0.07 0.01 1.50 73.97 -5.0 × 10-5 60.0 -10.5 6.29

costs breakdowns (MM$/year) annualized capital cost operating cost variability cost plant’s annualized cost

0.1069 0.5296 0.8032 1.4397

0.1066 0.9554 0.6193 1.6813

0.1067 0.5338 0.7786 1.4192

variable values and the plant’s costs obtained from the optimization of Scenario A’s formulation (HWA). As shown in this table, the resulting reactor’s design capacity is 2.3 times larger than the current capacity reported for this unit.25 The use of a larger reactor will allow for the specification of a higher pressure set point and a lower liquid level and temperature set point for the reactor. The combination of all these factors will contribute to reduce the amount of raw materials purged from the plant since a larger reactor has been specified. Therefore, an increase in the annualized capitals costs will result in a significant reduction of the operating costs. In the decentralized control strategy study,35 the annualized capital and operating costs for this plant were estimated to be 0.087 and 0.999 MM$/year, respectively. Although the present scenario specified a plant that has a higher annualized capital cost, the plant’s steady-state costs, i.e., capital are operating costs, were reduced by 40% when compared to the costs specified in that study. Therefore, the design and control scheme in this case study presents a more economically attractive design than that specified by a decentralized control strategy study.35 To validate the design in this case study, the optimal decision variable values η* obtained from the solution of problem 13 were used to simulate the closed-loop TE process model. The input to this simulation was the critical time-dependent profile in the A and C feed composition in stream 4, d*A and d*, C that produced the worstcase output variability in the product G at the optimal solution η* (Figure 2). Although the critical profiles in dA* and dC* obtained for other process variables could have been used to validate the design, the profiles shown in Figure 2 for dA* and dC* correspond to the situation where the minimum mass flow rate inequality constraint for product G is active at the optimal solution. Thus, it is expected that changes in this process variable may produce constraint violations. As shown in Figure 3, the actual maximum variability deviation in product’s G mass flow rate is equivalent to the minimum allowed product’s flow rate. The other set of critical realizations in dA* and dC* for the rest of the inequalities posed in problem 13 were also tested in the same fashion, but no constraint violations were observed. The results obtained by the present worstcase scenario-based solution method have corroborated that the proposed approach generates a design and control scheme that is robustly stable and feasible in the presence of the worst-case realizations in the A and C feed compositions of stream 4. 4.1. Bound-Based Worst-Case Approach (BWA). In the previous problem, the solution to the µ analysis shown in eq 6 was used to generate the worst-case realization in the disturbance set dA and dC, i.e., dA* and dC*. Then, these profiles were simulated to estimate the maximum deviation in a process variable around a nominal operating condition. Although the

Figure 2. Worst-case profiles in the disturbance set. Scenario A: worstcase scenario approach.

Figure 3. Product H mass flow rates. Scenario A: worst-case scenario approach.

hybrid worst-case methodology (HWA) reduces the conservatism in the solution, it increases the computational demands for this problem since the complete TE process nonlinear dynamic model must be simulated using the profiles in dA* and dC* at each function evaluation. A less computationally expensive alternative is to use the analytical bound estimated from the µ calculation (γ in problem 6) as the index to measure the variability of a process variable referred to in this manuscript as the BWA methodology.24 The BWA calculation method is expected to reduce the computational demands in the simultaneous design and control methodology by using analytical bounds on the variability instead of simulated values obtained with the mechanistic model but at the expense of added conservatism due to the use of a bound. To compare the worstcase scenario calculation method to the analytical bound-based approach, the problem formulated in eq 13 was solved using the analytical bounds (γ’s) as the performance index for each process variable within the cost function and constraints. The initial guesses and decision variables’ lower and upper limits were the same as those used in the HWA calculation. Table 2 shows a summary of the decision variables and the economic analysis obtained by using the BWA instead of the HWA presented in the previous section. As expected, this calculation


2829

Figure 4. Worst-case profiles in the disturbance set. Scenario A: analytical bound approach.

method specified a plant that is 14% more expensive than that defined by the HWA. This is because the BWA methodology uses analytical bounds (γ’s in eq 13) to calculate the worstcase process output variability and the process feasibility and variability inequality constraints, whereas the HWA method estimates these values from simulations of the actual closedloop model with the worst-case disturbance profiles, obtained from the analytical µ calculation. Thus, the application of the BWA strategy resulted in a more conservative design as compared to the HWA method. When compared to the HWA, the BWA calculation specified a plant with a slightly smaller reactor, a very low reactor’s pressure set point, and relatively high reactor’s liquid level set point. The combinations of these factors resulted in an operating cost that is 44% larger than that specified by the HWA calculation. To validate the design obtained by the BWA calculation method, the optimal decision variables values were used to simulate the closed-loop TE process model. The inputs to this simulation were the worst-case realizations in the A and C feed compositions in stream 4 (dA* and dC*), calculated from eq 7. The time-dependent profile used to describe the disturbance behavior is shown in Figure 4 and corresponds to the critical profile that produces the maximum variability in product H. The critical profile in dA* and dC* were obtained from the evaluation of the worst-case process variability for product H at the optimal solution, which corresponds to the case where the minimum mass flow rate constraint for product H specified in the problem’s formulation is active at the optimal solution. As shown in Figure 5, the maximum variability in H is far from the minimum constraint value specified for this variable. The rest of the process constraints considered in problem 11 were validated following the same procedure. 4.2. Dynamic Optimization-Based Approach. To illustrate the computational advantages of the proposed hybrid worstcase methodology (HWA), problem 13 was solved for comparison by a dynamic optimization problem. More specifically, the search for the worst disturbances and the corresponding variability in the variables involved in the cost function and the constraints was done by a dynamic optimization-type search instead of the µ-analysis-based search used in the HWA approach. Accordingly, at each iteration of the optimization, a numerical search was conducted for the critical values in the vector of disturbances dA and dC, e.g., d*A and d*, C that produces

Figure 5. Product G mass flow rates. Scenario A: analytical bound approach.

the maximum deviation of a given output for each new set of values of decision variables η. This is in clear contrast with the proposed HWA method that finds the critical disturbances based on analytical calculations. In the dynamic optimization approach both the critical realizations in the disturbances and the maximal variability resulting from these realizations are found by a numerical search. On the basis of the above descriptions, the optimization scenario was reformulated in eq 14 as a dynamic optimization problem. As shown, the decision variables for each dynamic optimization problem are the values at each sampling period i of the disturbance variables selected for this problem, i.e., dA and dC. It should be noted that the outer optimization problem is a nonlinear optimization problem since the vector of decision variables η is time independent. Therefore, the problem shown in eq 14 can be solved in MATLAB using the built-in function fmincon as the optimization solver as follows: (1) Given a set of values for the vector η and a set of values for the disturbance vectors dA and dC simulate the closed-loop dynamic behavior of the TE process using the process model equations and the controller algorithm equations specified for this system. (2) The largest variability in one particular output is sought from the simulation results. Since this value represents the problem’s objective function to be maximized, it is used by the optimiza-

2830


tion solver, fmincon, to choose a new set of critical values for the disturbance vectors. Go back to step 1. (3) This procedure is repeated until a convergence criterion is met. The solution to this optimization problem returns the worst-case realizations in the disturbance variables, d*A and d*C and the maximum deviation in a particular process variable that was produced by dA* and dC*. min CFTE η

s.t. TE Process model equations (9) 17 PI controller equations (Table 1) j + max G(t) - G* ΓG ) G dA∈dA

dC∈dC

j + max H(t) - H* ΓH ) H dA∈dA

Figure 6. Worst-case profiles in the disturbance set. Scenario A: dynamic optimization approach.

dC∈dC

xmG - 1.02xmH e 0 0.98xmH - xmG e 0 0.98 - xmG - xmH e 0 j + max G(t) e 0 G* - G dA∈dA

dC∈dC

j + max H(t) e 0 H* - H dA∈dA

dC∈dC

-0.05xjG + max xG(t) e 0

(14)

dA∈dA

dC∈dC

j r + max Pr(t) - 2895 e 0 P dA∈dA

dC∈dC

j r + max Lr(t) - 100 e 0 L dA∈dA

dC∈dC

j r + max Lr(t) e 0 30 - L dA∈dA

dC∈dC

j r + max Tr(t) - 150 e 0 T dA∈dA

dC∈dC

dA ) {dA(t), 0.475 e dA(t) e 0.495} dC ) {dC(t), 0.480 e dC(t) e 0.540} ηl e η e η u t ∈ [0, tf] It should be noticed that in this approach, for computational reasons, the stability of the control loops could not be explicitly assessed. Since the models are nonlinear, such stability checks will require the formulation of multiple Lyapunov-stability calculations that will result in prohibitive computational times. Problem 14 was solved using a sampling period (i) of 0.5 h and a finite time horizon (tf) of 20 h. Thus, each disturbance vector (dA and dC) is composed of 40 elements that are the decision variables for each of the dynamic optimization problems posed in eq 14. The optimization problem was initialized using the starting point used in the previous methods. Table 2 shows the results obtained by this method (dynamic optimization approach). As shown, the cost of the plant specified by this method is 15.6% less expensive than that specified by the bound-based worst-case approach (BWA) but only 1.4% less expensive than that proposed by the hybrid approach (HWA) proposed in the current work. Thus, the conservatism associated with either the BWA or the HWA approach is not significant for the present case study. Moreover, the results obtained by the dynamic optimization approach are consistent with those obtained with the previous methodologies. For instance, the dynamic optimization method specified a larger

Figure 7. Product G mass flow rates. Scenario A: dynamic optimization approach.

reactor which allows for the specification of a higher reactor pressure set point and a lower liquid level set point. Also, the production set point, which is directly related to the plant’s economics, was slightly reduced when compared to the specifications obtained by the other methodologies. The results obtained by the dynamic optimization-based methodology were validated by extensive simulations following the same procedure used for the previous methods studied in this work. Figure 6 shows the worst-case disturbance profile for A and C’s feed composition. These profiles generate the largest variability in the H’s mass flow rate at the solution point. This process constraint is active at the solution. This can be verified in Figure 7, where the profiles in Figure 6 drive H’s mass flow rate to the minimum allowed mass flow rate for this variable shown as red dashed lines in Figure 7. Simulations of the TE closedloop process model to the critical disturbances found by the dynamic optimization approach showed that the design and control parameters obtained by this method are asymptotically stable and feasible for the particular worst disturbance realizations. However, it cannot be guaranteed that other profiles in the disturbance set may generate an instability condition in the plant. This is because a formal stability check could not be included in the dynamic optimization formulation. 4.3. Analysis of the Computational Burden. The key advantage of the proposed HWA method as compared to the dynamic optimization approach is related to the computational burden. To illustrate the computational loads required by the solution methods studied in this work, the CPU time required to evaluate the cost function and the constraints imposed on each of

Ind. Eng. Chem. Res., Vol. 49, No. 6, 2010 Table 4. Design and Control Scheme: Scenario B

Table 3. Computational Burden: Scenario A

CPU time first function evaluation (s) averaged CPU time (s) averaged CPU time ratio total CPU time (h)

2831

HWA

BWA

dynamic optimization approach

232.34

64.25

12116.38

193.53 4.01 19.19

48.26 1.00 7.54

7198.32 149.15 206.8

the proposed optimization problems was recorded. The calculation methods used to solve the problem proposed in this section were solved in a Core 2 Duo CPU at 3.17 GHz and 2 GB of RAM. Table 3 shows a comparison of the different CPU times obtained by the solution methods studied in this work. The averaged CPU time was obtained from the CPU time needed to compute all the function evaluations divided by the total number of function evaluations required by the method to reach the optimal solution. As shown in Table 3, the averaged CPU time was less than the CPU time for first function evaluation in all cases. This is because the first function evaluation was started using an arbitrary guess for the disturbance vector. The subsequent function evaluations on each calculation method used as the starting point the critical realizations obtained in the previous function evaluations. Therefore, a shorter computational time is expected to decrease as the optimization progresses since the initial guess used for the inner optimization problems in eqs 13 and 14 is expected to be near the solution. As shown in Table 3, the BWA is as expected the most computationally efficient method. However, the solution obtained by this calculation method was shown to be conservative when compared to the other two methods. This suggests that the BWA could be used to quickly estimate a suboptimal design that can be further used as an initial guess for more computationally demanding methodologies, e.g., HWA or dynamic optimization-based strategy. On the other hand, the dynamic optimization method is by far the most computationally demanding method. Although the solution obtained by dynamic optimization is marginally less conservative than the HWA design by 1.42%, the average CPU time for the dynamic-based method is 37 times larger than the HWA method. Therefore, it can be concluded that the HWA approach estimated a plant that is economically attractive and did not require a significant amount of computational time. Table 3 also shows that the BWA requires 4 times less execution time than that required by the HWA but at the cost of conservatism in the solution. Furthermore, Table 3 shows that while the methodologies that use the robust models, i.e., BWA and HWA, required a CPU time in the order of hours to reach the optimum, the dynamic optimization method needed days to accomplish this task. Therefore, the application of dynamic optimization-based methodologies to larger systems may be challenging. The dynamic optimization method used in this work to compare the robust approach-based methods is based on the numerical integration of a nonlinear system for a finite period of time, i.e., sequential simulation and optimization. Although this may not be the most computationally efficient approach to solve the problem, it can be readily implemented. Thus, it is expected that the solution of the same problem using other dynamic optimization techniques39 may improve the results presented in this work. Similarly, the use of other dynamic optimization solvers, e.g., MINOPT, may influence the outcome of the comparisons.

design and control parameters (η)

HWA

BWA

dynamic optimization

reactor’s design capacity (ft3) reactor’s pressure set point (KPa) reactor’s liquid level set point (%) reactor’s temperature set point (°C) production set point (m3/h) Kc, purge valve (loop 5) τI, purge valve (loop 5) Kc, reactor liquid level(loop 11) τI, reactor liquid level (loop 11) Kc, reactor pressure (loop 12) τI, reactor pressure (loop 12) Kc, reactor temperature (loop 16) τI, reactor temperature (loop 16)

2972.8 2856.2 70.0 124.92 23.4 6.86 × 10-2 5.00 × 10-4 1.13 78.60 -8.21 × 10-4 16.34 -19.38 4.15

3000.0 2711.0 47.0 117.2 23.2 3.45e-2 9.75e-4 1.44 77.33 -2.22e-4 54.03 -1.50 15.0

1850.3 2888.5 57.5 118.1 23.1 0.005 2.7e-4 1.37 80 -1e-3 60.00 -1.00 8.10

cost breakdowns (MM$/year) annualized capital cost operating cost variability cost plant’s annualized cost

0.1066 0.6966 0.8945 1.6977

0.1069 0.5767 1.2271 1.9107

0.0936 0.6895 0.8901 1.6731

5. Scenario B: Reactor Section of the TE ProcesssChanges in A, B, and C Feed Composition The previous scenario assumed that the mole fraction of component B in stream 4 remained constant during the analysis, i.e., dB ) 0.005. To study the effect of additional disturbances on the simultaneous design and control problem, this second scenario considers the simultaneous design and control of the TE process reactor by considering a random magnitude bounded variation in reactants A and B and the inert component C in stream 4’s feed composition. The nominal mole fraction values and the upper and lower mole fraction limits specified for this scenario are as follows: d ) {dA, dB, dC} dA ) {dA, 0.475 e dA e 0.495} dB ) {dB, 0.002 e dB e 0.008} dC ) {dC, 0.480 e dC e 0.540}

(15)

Thus, at each function evaluation, robust FIR models have to be identified from the disturbance set dA, dB, and dC to the products’ mass flow rate, G and H, the reactor’s pressure temperature and liquid level, Pr, Lr, and Tr, and the mole fraction of G in the product’s stream, xG. Clearly, the addition of disturbances in the analysis increases the optimization problem’s curse of dimensionality. As in Scenario A, the present scenario was solved using MATLAB’s built-in function fmincon as the optimization solver. This problem was solved first using the BWA. The initial guesses for this problem were taken from the solution obtained from Scenario A for this calculation method. The optimal solution obtained by this method was then used as the starting point to solve Scenario B’s optimization problem using the HWA. Similarly, the solution obtained by the latter calculation method was used as the initial guess to solve the dynamic optimization method proposed in the previous section. Table 4 shows a summary of the results obtained by each solution method. As shown, the results show that the plant’s costs increased in the three cases as compared to the results in Scenario A. This result was expected since the addition of disturbances within the analysis impact the process economics. In this particular case study, it is necessary to increase the plant’s annualized cost by an average of 12% to account for the B composition’s variability in stream 4. The BWA estimated the most expensive plant, while the dynamic optimization strategy defined a plant that is 12.43% more economically attractive than that defined by the BWA method. However, the computational resources required to achieve the more economical solution were significant as in the previous

2832


Table 5. Optimal Tuning/Set-Point Problem: Scenario B decision variables (η)

HWA

reactor’s pressure set point (KPa) reactor’s liquid level set point (%) reactor’s temperature set point (°C) production set point (m3/h) Kc, purge valve (loop 5) τI, purge valve (loop 5) Kc, reactor liquid level(loop 11) τI, reactor liquid level (loop 11) Kc, reactor pressure (loop 12) τI, reactor pressure (loop 12) Kc, reactor temperature (loop 16) τI, reactor temperature (loop 16)

2893.83 38.16 125.0 23.32 1.76 × 10-2 7.74 × 10-4 1.5 75.63 -1 × 10-3 59.52 -8.42 15

cost breakdowns (MM$/year) annualized capital cost operating cost variability cost plant’s annualized cost

0.0869 0.8742 0.9578 1.9190

scenario. Moreover, the solution obtained by the HWA method defined a plant that is only 1.45% more expensive than the one defined by the dynamic optimization method. Likewise, the computational load required to perform one function evaluation using the HWA is on the order of minutes, while the dynamic optimization required hours to perform the same evaluation. Thus, the proposed hybrid worst-case methodology is found to be an attractive method that can be applied to solve the simultaneous design and control of large-scale systems. The results obtained by each method were validated by extensive simulations but are not shown here for brevity. 5.1. Scenario C: Optimal Tuning/Set-Point Selection ProblemsChanges in A, B, and C Feed Composition. To test the impact of redesigning the reactor capacity to reduce costs, the proposed hybrid worst-case method was used to solve a problem where only the controller tuning parameters and set points are optimized whereas the reactor capacity was kept constant at the original reported value equal to 1300 ft3. As in the previous scenarios, this problem was solved using the MATLAB tools. In addition, this optimization problem was initialized with the values given in Table 1. Table 5 presents the solution obtained from the optimal tuning/set-point selection problem with the reactor’s capacity fixed to 1300 ft3. The process operability and controllability parameters presented in Table 5 defined a plant with an annual cost of 1.919 that is approximately 11.54% higher than the annual cost of the plant obtained for Scenario B shown in Table 4 (HWA) of 1.6977. Thus, a more economically attractive design was obtained when both process design and process control were performed simultaneously as in Scenario B. This result highlights the fact that the interactions between design and control must be taken into account at the earlier stages of the design since they can significantly affect the plant’s economics. The major savings resulting from the simultaneous design and control approach are related to the operating costs, i.e., $178 000/year as compared to total savings of $221 300/year. As in the previous case studies, compliance with constraints using the optimal tuning/set-point parameters was verified through extensive simulations. 6. Conclusions This paper presented a new simultaneous design and control methodology (HWA) for large-scale systems based on the simulation of the worst-case scenario. To apply the method, the closedloop process model is represented here as an uncertain model. The uncertain model parameters and maximal variations expected in the disturbance set are used in a structured singular value (µ) analysis to calculate the worst-case realizations in the disturbance

set. Then, simulations of the closed-loop process model were performed using the worst-case disturbances profiles to obtain the actual maximum variability in the process variables that constrains the operation of the process. This procedure is used to evaluate the worst-case process output variability and test the process feasibility constraints in the simultaneous design and control approach. The proposed solution method was applied to the simultaneous design and control of the reactor section of the Tennessee Eastman process. The proposed worst-case approach (HWA) scenario method was compared to an analytical boundbased worst-case method (BWA) and a dynamic optimizationbased technique. Although the design specified by the HWA method is slightly conservative when compared to the solution obtained by the dynamic optimization-based method, the computational time required by the dynamic optimization technique is significantly higher than that required by the proposed HWA method. Therefore, the worst-case approach methodology presented in this work is a practical and computationally efficient tool for performing the simultaneous design and control of large-scale systems. Finally, the solution of an optimal controller tuning/setpoint problem illustrated that the integration of design and control may result in significant savings. Acknowledgment The authors would like to acknowledge the Mexican National Council for Science and Technology (CONACYT) and the Natural Science and Engineering Research Council of Canada (NSERC) for their financial support. Note Added after ASAP Publication: The version of this paper that was published online February 16, 2010 contained a production error in the tuning parameter discussion in section 3. The corrected version was reposted February 18, 2010. Appendix A. Structure of the Interconnection and the Perturbation Matrices M and ∆. The worst-case disturbance analysis test presented in Section 2.3 requires the specification of the interconnection and perturbation matrices, M and ∆. The structure of the interconnection matrix M is as follows M)

[

γ[Θ1, ..., Θq, ..., Θm]T 0 0 γI 0 0 [Ψ1, ..., Ψq, ..., Ψm] [Ω1, ..., Ωq, ..., Ωm] 0

]

(A1)

where the elements in this matrix are defined as follows Θq ) [δdq(0) δdq(1) ... δdq(j) ... δdq(N) ]

(A2)

Ωq ) [δh0,q δh1,q ... δhj,q ... δhN,q ]

(A3)

Ψq ) [h0,q h1,q ... hj,q ... hN,q ]

(A4)

where the matrix I is an (m(N + 1) × m(N + 1)) identity matrix. Note that M is an explicit function of the parameter γ, which value is obtained by the optimization algorithm used to solve the skewed µ-calculation problem shown at the RHS of eq 6. The perturbation block ∆ used in the µ-analysis has the following structure ∆ ) diag(∆r1, ∆r2, δc)

(A5)

where ∆r1 and ∆r2 are independent real scalar column vectors defined as follows


∆r1 ) [R1,1, ..., RN+1,1, ..., R1,q, ..., RN+1,q, ..., R1,m, ..., RN+1,m] ∆r2 ) [β1,1, ..., βN+1,1, ..., β1,q, ..., βN+1,q, ..., β1,m, ..., βN+1,m] (A6) Similarly, the term δc is a complex scalar. Both ∆r1 and ∆r2 are of lengths m(N + 1). The dimensions of the matrices M and ∆ are (2mN + 3) × (2mN + 3), respectively. The term βj,q in A6 represents a normalized value of the impulse response coefficient to the disturbance q at time interval j. That is, the vector ∆r2 represents the uncertainty associated with the impulse response coefficients, δhq. Literature Cited (1) Brengel, D. D.; Seider, W. D. Coordinated design and control optimization of nonlinear processes. Comput. Chem. Eng. 1992, 16, 861– 886. (2) Alhammadi, H. Y.; Romagnoli, J. A. Process design and operation Incorporating environmental, profitability, heat integration and controllability considerations. In The Integration of Process Design and Control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam, 2004; pp 264-305. (3) Luyben, M. L.; Floudas, C. A. Analyzing the interaction of design and control-1. A multiobjective framework and application to binary distillation synthesis. Comput. Chem. Eng. 1994, 18, 933–969. (4) Palazoglu, A.; Arkun, Y. A multiobjective approach to design chemical plants with robust dynamic operability characteristics. Comput. Chem. Eng. 1986, 10, 567–75. (5) Lenhoff, A. M.; Morari, M. Design of Resilient processing plants1. Process design under consideration of dynamic aspects. Chem. Eng. Sci. 1982, 37, 245–258. (6) Perkins, J. D.; Walsh, S. P. K. Optimization as a tool for design/ control integration. Comput. Chem. Eng. 1996, 20, 315–323. (7) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal design of dynamic systems under uncertainty. AIChE J. 1996, 42, 2251–2272. (8) Kookos, I. K.; Perkins, J. D. An algorithm for simultaneous process design and control. Ind. Eng. Chem. Res. 2001, 40, 4079–4088. (9) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Integrated flexibility and controllability analysis in design of chemical processes. AIChE J. 1997, 43, 997–1015. (10) Swartz, C. L. E. The use of controller parametrization in the integration of design and control. In The Integration of Process Design and Control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam 2004; pp 239-263. (11) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Ross, R.; van Schijndel, J. M. G. Simultaneous design and control optimization under uncertainty. Comput. Chem. Eng. 2000, 24, 261–266. (12) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. N. Recent advances in optimization-based simultaneous process and control design. Comput. Chem. Eng. 2004, 28, 2069–2086. (13) Schweiger, C. A.; Floudas, C. A. Interaction of design and control: optimization with dynamic models. In Optimal Control; Hager, W. H., Padarlos, P. M., Eds.; Kluwer Academic Publishers: The Netherlands, 1998; pp 388-435. (14) Fuente, D. L.; Flores-Tlacuahuac, A. Recent advances in optimization-based simultaneous process and control design. Ind. Eng. Chem. Res. 2009, 48, 1933–1943. (15) Malcolm, A.; Polan, J.; Zhang, L.; Ogunnaike, B. A.; Linninger, A. A. Integrating systems design and control using dynamic flexibility analysis. AIChE J. 2007, 53, 2048–2061. (16) Chawankul, N.; Ricardez Sandoval, L. A.; Budman, H.; Douglas, P. L. Integration of design and control: A robust control approach using MPC. Can. J. Chem. Eng. 2007, 85, 433–446.

2833

(17) Monnigmann, M.; Marquardt, W. Steady-state process optimization with guaranteed robust stability and feasibility. AIChE J. 2003, 49, 3110– 3126. (18) Grosch, R.; Monnigmann, M.; Marquardt, W. Integrated design and control for robust performance: Application to an MSMPR crystallizer. J. Process Control 2008, 18, 173–188. (19) Ricardez Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Simultaneous design and control of processes under uncertainty: A robust modelling approach. J. Process Control 2008, 18, 735–752. (20) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Application of robust control tools to the simultaneous design and control of dynamic systems. Ind. Eng. Chem. Res. 2009, 48, 813–813. (21) Seferlis, P.; Georgiadis, M. C. In The Integration of Process Design and Control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier, Ltd.: Amsterdam, The Netherlands, 2004. (22) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Integration of design and control for chemical processes: A review of the literature and some recent results. Annu. ReV. Control 2009, 33, 158–171. (23) Monnigmann, M.; Marquardt, W. Steady-state process optimization with guaranteed robust stability and flexibility: Application to HDA reaction section. Ind. Eng. Chem. Res. 2005, 44, 2737–2753. (24) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Simultaneous design and control of chemical processes with application to the Tennessee Eastman process. J. Process Control 2009, 19, 1377–1391. (25) Downs, J. J.; Vogel, E. F. Plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245–255. (26) Ramirez, E.; Gani, R. Methodology for the design and analysis of reaction separation systems with recycle. 2.- Design and control integration. Ind. Eng. Chem. Res. 2007, 46, 8084–8100. (27) Wu, K.; Yu, C. Operability for processes with recycles: interaction between design and operation with application to the Tennessee Eastman challenge process. Ind. Eng. Chem. Res. 1997, 36, 2239–2251. (28) Ljung, L. System Identification: Theory for the User, 1st ed.; Prentice Hall: New York, 1987. (29) Doyle, J. Analysis of feedback systems with structured uncertainties. IEEE Proc., Part D: Control Theory Appl. 1982, 129, 242–250. (30) Fan, M. K. H.; Tits, A. L.; Doyle, J. C. Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics. IEEE Trans. Autom. Control 1991, 36, 25–38. (31) Braatz, R. D.; Lee, J. H.; Morari, M. Screening plant designs and control structures for uncertain systems. Comput. Chem. Eng. 1996, 20, 463–468. (32) Lee, J. H.; Braatz, R. D.; Morari, M.; Packard, A. Screening tools for robust control structure selection. Automatica 1995, 31, 229–235. (33) Ma, D. L.; Chung, S. H.; Braatz, R. D. Worst-case performance analysis of optimal batch control trajectories. AIChE J. 1999, 45, 1469– 1476. (34) Nagy, Z. K.; Braatz, R. D. Worst-case and distributional robustness analysis of finite-time control trajectories for nonlinear distributed parameter systems. IEEE Trans. Control Syst. Technol. 2003, 11, 694–704. (35) Ricker, N. L. Decentralized control of the Tennessee Eastman challenge process. J. Process Control 1996, 6, 205–221. (36) Larsson, T.; Hestetun, K.; Hovland, E.; Skogestad, S. Selfoptimizing control of a large-scale plant: The tennessee eastman process. Ind. Eng. Chem. Res. 2001, 40, 4889–4901. (37) Seider, W. D.; Seader, J. D.; Lewin, D. R. Product and process design principles, 2nd ed.; John Wiley and Sons, Inc.: New York, 2004. (38) Optimization Toolbox, For Use with MATLAB; The MathWorks, Inc.: Natick, MA, 2006. (39) Bertsekas, D. P. Dynamic Programming and Optimal Control; Athena Scientific: Nashua, 2007; Vols. I and II.

ReceiVed for reView July 2, 2009 ReVised manuscript receiVed October 19, 2009 Accepted January 22, 2010 IE9010707

Simultaneous Design and Control: A New Approach and Comparisons

Recommend Documents