Simultaneous Design and Control under Uncertainty Using Model

Article pubs.acs.org/IECR

Simultaneous Design and Control under Uncertainty Using Model Predictive Control Kelvyn B. Sanchez-Sanchez and Luis A. Ricardez-Sandoval* Department of Chemical Engineering, University of Waterloo, N2L 3G1, Waterloo Canada S Supporting Information *

ABSTRACT: This work presents a new methodology for the simultaneous process flowsheet and control design of dynamic systems under uncertainty using a model predictive control (MPC) strategy. Although several methodologies that include structural decisions in the analysis have been reported, only a few have considered a model-based control strategy such as MPC. An iterative decomposition framework that includes a dynamic flexibility analysis, a robust dynamic feasibility test, a nominal stability analysis, and a robust asymptotic stability test is presented for optimal process flowsheet selection. While previous methodologies have formulated the dynamic feasibility and stability analysis as a mixed-integer nonlinear programming (MINLP) problem, the present method formulates these analyses as convex problems for which efficient numerical algorithms exist. The simultaneous process flowsheet and MPC design method was tested on a system of Continuous Stirred Tank Reactors (CSTRs) with multiple inlet streams. The results show that the optimal design obtained by the present method remained feasible and asymptotically stable in the presence of the critical realizations in the disturbances. Comparisons between the designs obtained by the present MPC-based method and those obtained with other design approaches, i.e., optimal steady-state design and simultaneous design and control using a multiloop proportional and integral (PI) control scheme, are presented and discussed in this work.

1. INTRODUCTION It has been widely recognized that efficient and low-cost integrated chemical plants with tight process constraints and stringent environmental and safety regulations are needed to satisfy the current market demands. The integration of chemical processes is not an easy task to perform since the transient operation of a system is usually affected by factors such as disturbances with high and low frequency content, e.g., sudden or seasonal variability in the upstream processes, and uncertainty in the physical parameters, e.g., uncertainty in the density of a fluid. These undesirable though realistic operating conditions may lead to product variability, process infeasibility, and process instability. Previous studies on process design have shown that a controllability analysis is essential to design economically attractive chemical plants that can accommodate critical realizations in the process disturbances and uncertainty in the physical parameters.1−5 Integration of design and control has been proposed as an attractive alternative to overcome the issues posed by the sequential process design approach. The fundamental idea in simultaneous design and control is to search for both the process flowsheet topology and the control structure that complies with the process operational constraints while meeting the product specifications at minimum cost. Several methodologies have been proposed for simultaneous design and control. An in-depth review on the features and limitations of those methodologies can be found elsewhere.6−8 While simultaneous design and control methodologies that incorporate advanced control strategies within the analysis are currently available, e.g., model predictive control (MPC),9−13 model-based parametric controllers,6,14,15 embedded control optimization,16,17 they have not been widely studied as those based on conventional feedback controllers.18−24 Also, most of © 2013 American Chemical Society

the MPC-based approaches reported in the literature for optimal design and control have only considered the optimal sizing of the process units,9−13 i.e., the process flowsheet remained fixed during the calculations. To the author’s knowledge, the approach proposed by Sakizlis and coworkers6,14 is the only method that has embedded structural decisions and an MPC strategy within a simultaneous design and control framework. This paper presents a new integration of process flowsheet and control design methodology that incorporates a model predictive control (MPC) strategy in the analysis. An iterative decomposition strategy comprising of a dynamic flexibility analysis, a robust dynamic feasibility analysis, a nominal stability analysis, and a robust asymptotic stability analysis is presented here for optimal process flowsheet design using a multivariable MPC scheme. Both the robust dynamic feasibility analysis and the robust asymptotic stability analysis are formulated as convex problems. Thus, the proposed method is computationally attractive since recent optimization-based methodologies formulate the dynamic feasibility and stability analyses as a mixed-integer nonlinear programming (MINLP) problem.6,21,22 The robust dynamic feasibility test implements a norm-bounded analysis that computes the critical realizations in the disturbances that produce the largest (worst-case) variability in the outputs. Likewise, process asymptotic stability is enforced by adding a formal asymptotic stability test. These analyses enable the specification of an optimal design that Received: Revised: Accepted: Published: 4815

August 17, 2012 January 28, 2013 February 16, 2013 February 17, 2013 dx.doi.org/10.1021/ie302215c | Ind. Eng. Chem. Res. 2013, 52, 4815−4833

Industrial & Engineering Chemistry Research

Article

Figure 1. Iterative decomposition algorithm for simultaneous process synthesis and control design.

the most promising process flowsheet candidates using heuristic methods or formal mathematical formulations, which are based on steady-state information. Instead, the present method aims to perform an advanced (detailed) process synthesis analysis on those promising process flowsheets, usually obtained from a preliminary (steady-state) process synthesis analysis, such that, when combined with a model-based control strategy such as MPC, an optimal process flowsheet that can accommodate the most critical time-varying realizations in the disturbances can be identified. This paper is organized as follows: section 2 presents the iterative decomposition algorithm proposed in this work for the optimal process flowsheet and control design under uncertainty

remains feasible and asymptotically stable despite critical realizations in the disturbances. A few simultaneous design and control approaches have been previously developed by one of the authors.25−28 In those works, the process flowsheet and the control structure was assumed to be fixed during the analysis. The approach presented here expands upon those previous methodologies by considering structural decisions combined with a multivariable MPC strategy for optimal process design under uncertainty. While the present methodology incorporates structural decisions in the analysis for the selection of an optimal process flowsheet, the approach proposed here is not intended to replace an early process synthesis analysis that aims to identify 4816

dx.doi.org/10.1021/ie302215c | Ind. Eng. Chem. Res. 2013, 52, 4815−4833


Article

interval is available. Likewise, the vector Δûn+q|n represents the future moves in the nu manipulated variables at the (n+q)th time interval that need to be implemented in the actual plant to maintain the output (controlled) variables y close to the reference signal r. The scalars P and Q are the prediction and control horizon, respectively, whereas r̂ ∈ 9 1xnr represents the reference signals in deviation form and are assumed to remain constant during the solution of the MPC formulation. The reference signals (r ∈ 9 1xnr) represent optimization variables in the present simultaneous design and control formulation and will be obtained by the dynamic flexibility analysis. As shown in (1), the prediction (P time steps ahead) in the output (controlled) variables will be obtained from a linear discrete state space model. The linear state space model represents the internal model used by the MPC algorithm to provide estimates for the future control moves, i.e., Δû. This model captures the key process dynamics characteristics around a nominal operating condition specified by d, the vector of real process design variables that cannot be adjusted during the operation (d ∈ 9 1xnd); θ, the vector of nθ integer variables that includes the different process synthesis alternatives (θ ∈ ), e.g., the number of trays in a distillation column; unom, the nominal (optimum) values in the manipulated variables (unom ∈ 9 1xnu); and the process set points (r). The states of the linear discrete state space model at the nth sampling interval are given by x̂n. A linear discrete state observer based on the linear internal process model used by the MPC algorithm is employed in (1) to estimate the system’s states at any sampling interval (see section 2 in the work of Maciejowski29). The key inputs to the linear state observer are the output measurements from the system at the nth (current) time step. Following the MPC algorithm (1), ûn and V̂ n are vectors that represent the nu manipulated variables to be used by the MPC algorithm and the nv disturbances affecting the process at the nth time interval, respectively. While the manipulated variables are assumed to be changing up until the Qth time step predictions ahead, the process disturbances in the MPC algorithm are assumed to remain constant for the Q time intervals ahead and equal to its value at the nth (current) sampling interval. As it will be described in section 2.2, the realizations of the process disturbances at the nth time step will be obtained from the block matrix Vk, which is an input to the dynamic flexibility analysis (see Figure 1). The internal MPC model can be obtained from systems identification techniques such the leastsquares method30 or from Taylor series expansions.31 While ŷ represents the prediction obtained for the ny outputs from the linear state space model at different time intervals, y(t) ∈ 9 ny represents the actual ny output measurements obtained from the complete mechanistic process model, i.e., ŷn represents an approximation to the actual output measurements at the nth time step. Note that ûn, ŷn, Δûn, and V̂ n are in deviation variables and defined with respect to a nominal operating point defined by d, θ, unom, and r, respectively. Also, linear inequality constraints on ûn, ŷn, and Δûn are included in (1) to ensure that the control moves (Δû) specified by the MPC scheme satisfy the operational limits on these variables. The elements of the positive semidefinite matrices Ψ and Λ represent the weights assigned to the controlled and the manipulated variables, respectively. These weights represent the tuning parameters in the MPC scheme and are considered as optimization variables in the dynamic flexibility analysis. The matrices Ψ and Λ are typically assumed to be diagonal matrices and are defined as follows:

using MPC. Section 3 introduces a system of Continuous Stirred Tank Reactors (CSTRs) in series, which is used to illustrate the application of the present simultaneous design and control methodology. The results obtained for the present case study as well as comparisons made with other approaches are presented and discussed at the end of that section. Concluding remarks are stated at the end. Notation. 9 denotes the set of real numbers whereas  represent the set of integer numbers. Variables in bold, e.g., d or Δ, denote a vector or a matrix whereas variables denoted as x(t) indicate the realizations of a single variable x as a function of time t, i.e., x(t) ∈ 9 . For example, x(t) can represent the time variations of a state variable x due to time-varying changes in the disturbances. Likewise, x(t) is a time-dependent vector that includes the realizations of all the variables x (denoted as nx) as a function of time t, i.e., x(t) ∈ 9 nx. For example, x(t) can represent the time response of all the nx state variables of a system due to time-varying fluctuations in the disturbances.

2. ITERATIVE DECOMPOSITION FRAMEWORK This section presents the approach proposed in this work to search for the optimal process flowsheet design and the model predictive control (MPC) strategy that can accommodate the critical realizations in the process disturbances at minimum cost. As shown in Figure 1, the iterative decomposition strategy proposed in this work consists of four analyses: a dynamic flexibility analysis, a robust dynamic feasibility analysis, a nominal stability analysis, and a robust asymptotic stability analysis. The assumptions and the initialization step are presented next followed by the mathematical formulations proposed for each analysis included in the algorithmic framework. 2.1. Assumptions and Initialization. The present methodology assumes that the dynamic behavior of the process flowsheet alternatives can be obtained from simulations of their corresponding process models. Also, the disturbances are considered as random time-varying perturbations bounded by upper and lower limits. The control scheme included in the present method is assumed to be a multivariable linear constrained model predictive controller. This control strategy is mathematically formulated as follows: P

min

∑ (yn̂ + p | n − r)̂ T Ψ(yn̂ + p | n − r)̂

Δû p=1

Q −1

+

∑ Δû n+ q| n TΛΔû n+ q| n q=0

s.t. x̂ n + 1 = A(d, θ , u nom , r)x̂ n + B(d, θ , u nom , r)ĥ n yn̂ = C(d, θ , u nom , r)x̂ n ĥ n = [û n, Vn̂ ]T û min ≤ û n ≤ û max Δû min ≤ Δû n ≤ Δû max ymin ̂ ≤ yn̂ ≤ ymax ̂

(1)

where ŷn+p|n is a vector that represents the prediction for the ny output (controlled) variables at the (n+p)th time interval assuming that the actual output measurement at the nth time 4817


Industrial & Engineering Chemistry Research λ = [diag(Ψ), diag(Λ)]T

Article

analysis but that failed the robust dynamic feasibility test. Likewise, the nominal or asymptotically unstable solutions included in the block matrix Φ are identified from the nominal or the robust stability tests as it will be shown in section 2.4. The following constraints are included in the dynamic flexibility analysis to avoid the reselection of solutions that did not comply with the stability tests:

(2)

where λ is a vector that includes the weights on the manipulated and controlled variables (λ ∈ 9 1xnλ). The linear constrained MPC formulation (1) can be cast as a constrained least-squares optimization problem that can be reformulated as a quadratic optimization (convex) problem29 that can be efficiently solved using standard quadratic programming (QP) optimization algorithms, e.g., CPLEX, MINQ, KNITRO.32 Also, numerical subroutines that implement linear constrained MPC algorithms are currently available on commercial software packages, e.g., the MPC toolbox in MATLAB.33 As shown in Figure 1, previous solutions considered as optimal that did not comply with the nominal or the robust asymptotic stability tests cannot be considered in the next dynamic flexibility analyses. Thus, the present algorithm introduces a block matrix (Φ) that includes those process flowsheet alternatives and MPC schemes that did not satisfy the stability tests. The block matrix Φ is defined as follows: Φ = [Φ(1), ..., Φ(ii), ..., Φ(ι)]

∀ ii = 1, ..., ι

⎛ λ − λ (ii) ⎞ Φ − ελ⎟⎟ ≥ 0 (1 − Ζii)⎜⎜ ⎝ max(λ , λ Φ(ii)) ⎠

∀ ii = 1, ..., ι

⎛ u ⎞ nom − unomΦ(ii) − εu⎟⎟ ≥ 0 (1 − Ζii)⎜⎜ ⎝ max(unom , unomΦ(ii)) ⎠ ∀ ii = 1, ..., ι ⎛ r − r (ii) ⎞ Φ − εr ⎟⎟ ≥ 0 (1 − Ζii)⎜⎜ ⎝ max(r , rΦ(ii)) ⎠

ii = 1, ... ι

Φ(ii) = [d Φ(ii), θΦ(ii), λΦ(ii), u nomΦ(ii), rΦ(ii)]T

⎛ d − d (ii) ⎞ Φ − εd⎟⎟ ≥ 0 (1 − Ζii)⎜⎜ ⎝ max(d , dΦ(ii)) ⎠

(3)

where Φ(ii) is a process synthesis and MPC design configuration obtained from a previous dynamic flexibility analysis that did not satisfy the stability tests included in the proposed method. As shown in (3) and in Figure 1, in the event that the current optimal solution obtained from the dynamic flexibility analysis at the kth iteration (ηk) is identified as unstable, then the index set ι is increased by one and the matrix Φ is updated (augmented) with a new column that includes the current solution, i.e., Φ(ι) = ηk. The index set ι shown in (3) is used to keep track of the solutions that failed the stability tests and is different from the overall iteration index k used in the iterative decomposition framework (see Figure 1). Based on the above, the present algorithm is initialized as follows: • Set k = 0, ii = 0, ι = 0, Φ = []. • Specify the upper and lower bounds for the disturbances (vlow and vup) and initial time-discrete realizations for the nv disturbances considered in the analysis, which is denoted as V0. • Define the manipulated and controlled variables included in the MPC algorithm, i.e., u and y. • Specify an initial process synthesis configuration and weights for the MPC tuning parameters, i.e., d0, θ0, r0, unom0, λ0. For simplicity, the initial guesses for these variables are referred to as η0. 2.2. Dynamic Flexibility Analysis. The dynamic flexibility analysis aims to search for the process flowsheet alternative and the MPC tuning parameters that comply with the process constraints at minimum cost. As shown in Figure 1, the two key inputs to perform this analysis are the critical time-discrete realizations in the nv disturbances that have been collected up to the kth iteration (denoted as Vk), and the block matrix Φ, which includes the previous process flowsheet and MPC design alternatives that failed the stability tests. As shown in section 2.3, the critical realizations included in Vk are identified from the robust dynamic feasibility analysis and include the previous and the current critical realizations in the disturbances that produce violations to the process feasibility constraints. This specification for Vk avoids the reselection of previous solutions that were identified as optimal by the dynamic flexibility

∀ ii = 1, ..., ι

Ζ = [Ζ1, ..., Ζii, ..., Ζι] Ζii = tanh{M(∑ (θ − θΦ(ii))2 )} θ

(4)

The symbol Zii is a smooth approximation function that activates the inequality constraints on the process design parameters (d), the MPC tuning parameters (λ), and the nominal operating conditions (r, unom) when a process synthesis alternative (θ) is the same to that included in the block matrix Φ. The scalar parameters εd, εr, εu, and ελ represent tolerance criterions used to determine how close the dynamic flexibility test can select values for d, r, and λ with a process synthesis alternative (θ) that is included in Φ. The values assigned to εd, εr, εu, and ελ depend on the problem under consideration, the computational resources available, and the desired accuracy in the optimal process synthesis and MPC design. Setting small values for these tolerance parameters will enable the dynamic flexibility analysis to select values for d, r, unom, and λ closer to previous unstable solutions. However, this also increases the chances of specifying a solution that may be unstable and thus increasing the computational costs for this problem since more iterations will be needed to achieve an optimal feasible and stable solution. On the other hand, setting relatively large tolerance parameter values will reduce the computational demands at the expense of specifying a feasible and stable suboptimal solution. Accordingly, there is a tradeoff in the specification of the values for εd, εr, εu, and ελ. On the basis of the case studies considered with the present algorithm, setting a value between 0.05 and 0.1 for εd, εr, εu, and ελ seems to be appropriate. For example, if the dynamic flexibility analysis selects a process flowsheet configuration θ that is included in Φ, e.g., θΦ(ι), then the inequality constraints shown in (4) ensure that the dynamic flexibility analysis can only assign values to d, λ, unom, and r that are at least 5−10% different from dΦ(ι), λΦ(ι), unomΦ(ι), and rΦ(ι), respectively. The parameter M included in Zii is a sufficiently large number. Enforcing constraints (4) ensure that previous solutions identified as unstable cannot be reselected in the next flexibility analysis. 4818



Article

On the basis of the above, the dynamic flexibility analysis formulation proposed in this method is as follows: min

d, r, θ , λ , u nom

critical realizations of a specific disturbance that have been identified up until the (k − 1)th robust dynamic feasibility tests. Likewise, each column in Vk represents the time-discrete realizations in the disturbances that produce violations to the feasibility constraints at a given iteration step. As shown in section 2.3, the worst-case scenario in the process feasibility constraints is estimated under the assumption that the critical realizations in the disturbances are acting on the system at different times. Thus, the cost function (Θ) and the feasibility constraints (gκ) in (5) are evaluated in the presence of the critical realizations of each individual disturbance. The procedure to evaluate Θ and gκ in (5) is presented in the Supporting Information. To reduce (5)’s dimensionality, the process model equations f l’s can be approximated to a finite number of discrete intervals using a piecewise (typically constant) function (control vector parametrization) or using orthogonal collocation on finite elements (complete discretization of the dynamic system). This approximation reduces (5) into a multilevel mixed-integer nonlinear optimization problem, which can be solved using efficient optimization algorithms, e.g., BARON, outer approximation, DICOPT.32,34 Nevertheless, the main shortcoming of the most efficient optimization methods available for multilevel MINLP is that they may return local minima. Although multistart techniques or global MINLP optimization algorithms are available,32 the main shortcoming of those methods is their computational costs. Thus, the quality of the solution obtained from the dynamic flexibility analysis depends on the problem under consideration, the computational resources available and the algorithm used to solve the multilevel MINLP problem posed in (5). Simultaneous design and control methods that reduce the complexity in the calculations in the dynamic flexibility analysis have recently emerged in the literature16,17 and are the current subject of research in this area. The dynamic flexibility analysis formulation (5) returns a process flowsheet configuration and an MPC strategy that comply with the process constraints (gκ’s) at minimum cost (Θ) for a given (fixed) critical time-discrete profile in the disturbances, Vk. Hence, the optimal design obtained from the dynamic flexibility analysis at the kth iteration (ηk) needs to be further tested to ensure that other critical realizations in the disturbances will not violate the process feasibility constraints. 2.3. Robust Dynamic Feasibility Analysis. A key challenge in simultaneous design and control is the specification of the critical time-dependent trajectories in the disturbances that generates the largest (worst-case) deviation in the time-varying process variables, i.e., the worst-case scenario. The specification of the critical realizations in the disturbances that produce the worst-case scenario in the process variables is not a simple task to perform since it involves the search of the critical realizations of the disturbances at each time interval, i.e., a critical time-trajectory, that, when simulated in closed-loop, produce the largest deviation in the process variables with respect to a nominal operating (reference) point. The search of such critical time-dependent disturbance profile is challenging even for simple systems.25,26 Previous simultaneous design and control methodologies have addressed this issue by assuming a particular disturbance profile in the analysis, i.e., the disturbance dynamics are assumed to know a priori, e.g., a series of steps15,23 or sinusoidal function with known amplitude and frequency.9 Other methods have used a predefined disturbance time-dependent function with unknown (critical) model parameters, e.g., a sinusoidal function with unknown amplitude

Θ(ẋ(t ), d, u nom , y(t ), θ , Vk , x(t ), r, u(t ))

subject to fl (ẋ(t ), d, r, θ , u nom , y(t ), x(t ), Vk , u(t )) = 0, l∈L gκ (ẋ(t ), d, r, θ , u nom , y(t ), x(t ), Vk , u(t )) ≤ 0, κ∈Κ MPC scheme (1) constraints (4) d = {d|dl ≤ d ≤ du} r = {r |rl ≤ r ≤ ru} λ = {λ|λl ≤ λ ≤ λ u} u nom = {unom|unom,l ≤ unom ≤ unom,u} θ = {0, 1}1xnθ t = [0, tf ]

(5)

where the function f l denotes the nonlinear differential and algebraic equations that describe the behavior of the process design alternatives, e.g., mass, energy, and mass component balance equations; gκ represents the process feasibility constraints, e.g., the temperature inside a reactor must not exceed the maximum allowed working temperature at any time t to avoid a thermal runaway; u(t) ∈ 9 nu represents the nu manipulated variables that can be adjusted for control; x(t) ∈ 9 nx represents the nx state variables of the system, e.g., x(t) can be a molar hold-up. As shown in (5), the integer process design variables θ are reformulated as binary variables, i.e., θ ∈ {0, 1}. The cost function Θ is an economic index that evaluates different process design alternatives. This function is problemspecific and usually defined in terms of the process steady-state economics (capital and operating costs) and the performance costs of the system in closed loop (variability costs).26,27 The block-matrix Vk, which is an input to (5), includes the critical time-discrete realizations in the disturbances and is defined as follows: ⎡ v1,1 v2,1 ⎢ ⋮ ⎢ ⋮ ⎢ v v 2, i Vk = 1, i ⎢ ⋮ ⎢ ⋮ ⎢v v ⎣ 1, nv 1, nv

... vk − 2,1 vk − 1,1 ⎤ ⎥ ⋱ ⋮ ⋮ ⎥ ... vk − 2, i vk − 1, i ⎥ ⎥ ⋱ ⋮ ⋮ ⎥ ... vk − 2, nv vk − 1, nv ⎥⎦

(6)

where vk−1,i is a row vector of length nc(Nk−1 + 1) whose elements include the time-discrete realizations in the ith disturbance that produce violations to nc process feasibility constraints and that are identified from the robust dynamic feasibility analysis performed at the (k − 1)th iteration step. Nk−1 represents the final sampling instant that corresponds to the closed-loop process settling time of the process design and MPC scheme configuration obtained from the dynamic flexibility analysis at the (k − 1)th iteration step, i.e., ηk−1 (see Figure 1). As shown in (6), each row in Vk denotes the 4819



Article

and frequency.6,21,24,35 In this work, the critical time-dependent trajectories in the disturbances are obtained from the solution of a structured singular value (SSV) analysis formulation, which is posed as a convex problem that can be solved using computationally efficient numerical algorithms. To implement the robust dynamic feasibility analysis, the process disturbances are defined as follows: V = [v1, ..., vi, ..., vnv]T ;

viup + vilow ; 2

sj , i = sj̅ , i ± δsj , i

δvi =

n = 0, 1, 2, ..., N viup − vilow 2

vi̅ = {vi̅ (n)| vi̅ (n) = vi(n) − vi ,nom} (7)

where vi(n) is a scalar that represents the ith process disturbance at the nth time interval; vi,nom and δvi are scalars that represent the nominal (steady-state) value and the maximal deviation expected for the ith disturbance, respectively. N represents the final sampling instant, i.e., tf = NΔt; whereas vi̅ represents the ith external perturbation in deviation form. As shown in (7), the disturbances are not limited to follow a particular time trajectory as it has been proposed in previous methodologies.6,9,10,21−24,35 Note that V in (7) represents a general time-discrete description of the disturbances, i.e., the disturbance values at any sampling instant are not specified, whereas the disturbance description Vk presented in (6) and used to solve (5) represents a particular time-discrete realization of the disturbances at a specific iteration step, i.e., the disturbance values at each sampling instant are known. The goal in the present robust dynamic feasibility analysis is to obtain the critical realizations in the nv disturbances at the kth (current) iteration step, i.e., vk,1, ..., vk,i, ..., vk,nv, that produces the largest (worst-case) deviation in the process variables with respect to a nominal operating point, i.e., the worst-case scenario. Consider a single time-dependent output (controlled) variable y(t) ∈ y(t); the critical realizations in the disturbances that produces the worst-case scenario for y(t) is denoted as vy* ∈ 9 (N+1)xnv. To calculate vy*, the present method requires that the process dynamics that describes the closedloop transient behavior of the system specified by ηk can be represented by a linear impulse response process model complemented with uncertain model parameters, i.e., an uncertain finite impulse process model. This uncertain discrete-time model is as follows:

max || y(t )||∞

v(t ) ∈ V(t )

s.t. fl (ẋ(t ), ηk , y(t ), x(t ), V(t ), u(t )) = 0,

∑ sj ,ivi̅(n − j);

l∈L

MPC scheme (1) (10)

where ∥y(t)∥∞ represents the infinity norm over all the timedependent realizations in an output (controlled) variable y due to critical time-dependent realizations in the disturbances at the kth iteration step around a process flowsheet and MPC scheme configuration specified by ηk. The direct solution of (10) can become challenging for large-scale complex systems. The main difficulty lies in the fact that the dynamic behavior of the process synthesis alternatives are usually represent by nonlinear process models. However, the uncertain impulse response model shown in (8) has a linear model structure and represents the process model equations and the MPC strategy included in (10). Hence, (10) can be reformulated as follows:

N

y imp( n) = i

(9)

where sj̅ ,i is scalar that represents the nominal value of the impulse response coefficient at the jth time step due to changes in the ith disturbance. Likewise, δsj,i is the corresponding uncertainty associated to that impulse response coefficient. The uncertain terms δsj,i’s included in the impulse response model account for the process nonlinearities that cannot be accurately described by the nominal impulse response coefficients, sj̅ ,i. According to (9), sj,i can take any value that satisfies that model parameter description; this implies that the uncertain dynamic model shown in (8) includes an infinite set of linear impulse response models, each evaluated at specific values of sj,i that satisfies (9). Similarly, the use of an uncertain dynamic model to represent the actual closed-loop dynamic behavior of the process and control design alternative specified by ηk adds conservatism due to the addition of the uncertain terms δsj,i into the impulse response model. However, the conservatism in the solution is expected to be relatively small because the control actions implemented by the MPC scheme aim to maintain the system closed to its nominal operating point thus reducing the system’s variability and consequently the magnitude of the uncertain model parameters,δsj,i’s. The output y is an output (controlled) variable that appears in a process feasibility constraint, i.e., gκ’s in (5). Note that the only inputs to the uncertain dynamic models included in yimp are the disturbances affecting the process. Thus, the uncertain dynamic model (8) is only valid around a nominal operating point defined by the current process flowsheet and MPC design alternative (ηk). The procedure to identify uncertain models such as (8) is presented in the Supporting Information. The critical realizations in the disturbances that produce the worst case scenario in y around an operating condition ηk can be obtained as follows:

i = 1, ... nv

vi = {vi(n)|vilow ≤ vi(n) ≤ viup}; vi ,nom =

disturbance. The impulse response model (8) is uncertain because the impulse response coefficients sj,i are uncertain model parameters, i.e.,

∀ n = 0, 1, 2, ..., N

j=0

(8)

where yimpi(n) is a scalar that denotes the impulse response of output y (in deviation form) at the nth sampling instant due to changes in the ith disturbance whereas vi̅ (n − j) represents the ith disturbance value at the (n − j)th time interval in deviation form. Therefore, yimpi represents the uncertain finite impulse response model that describes the system’s closed-behavior between a controlled variable y and the ith disturbance vi̅ . Accordingly, yimp: 9 nv → 9 is a vector function that represents a family of uncertain impulse response models; each uncertain model included in yimp represents the closed-loop dynamic behavior between an output y and a specific process

max || yimp ||∞ vi ∈ V

(11)

where V and vi are the time-discrete realizations of the disturbances defined in (7). (11) cannot be directly solved since each of uncertain models included in yimp poses an 4820



Article

As shown in (12), the structure of the interconnection matrix H includes the scalar parameter χ, which represents the only unknown in the SSV formulation. Moreover, the μΔ(H) problem in (12) is linear since the elements included in H and Δ are linear.38 Thus, the SSV problem (12) can be cast as a convex roots problem that returns a single (global) solution. This is a key feature in the present approach since efficient numerical methods, which require a minimum amount of computational effort, are currently available for solving convex roots problems, e.g., Newton−Raphson method, Secant method. The value in χ that satisfies (12) represents a bound on the worst-case variability expected for a time-varying process output (controlled) variable around a nominal operating point. The bound, denoted as χy for an output variable y, can be used to measure the worst-case variability for the output y due to the critical time-discrete realizations in the disturbances identified for the current process design and MPC scheme (ηk). Nonetheless, the elements in Δ that produce the worst-case scenario can be extracted from the evaluation of μΔ(H) at χ = χy; hence, the norm-bounded elements in Δv can be extracted and used to specify the actual critical realizations in the disturbances that produces the worst-case variability in y, i.e.,

infinite set of linear impulse response models. However, analytical bounds for ∥yimp∥∞ can be estimated if (11) is reformulated as a structured singular value (SSV) analysis problem,36,37 i.e., max || yimp ||∞ ≤ χ ⇔ μΔ (H) − χ = 0 vi ∈ V

(12)

where the real scalar variable χ represents a bound on the maximum deviation expected for y with respect to a process design and MPC scheme specified by ηk. The SSV analysis, i.e., μΔ(H) in (12), is a linear technique38 that has been widely used in control theory to design robust controllers.39−41 In the present analysis, the application of this technique is extended to compute bounds on the largest variability of a process variable, which can then be used to perform a robust dynamic feasibility analysis, i.e., evaluate the compliance of the process feasibility constraints in the presence of the worst-case scenario. As shown in (12), the SSV formulation is a function of the interconnection matrix H, whose structure is shown in (13). Π is a vector of ones of length N + 1; δvi, the maximal deviation in the ith disturbance, is calculated from (7) whereas sj̅ ,i and δsj,i are obtained from systems identification. In addition, I is an identity matrix of dimensions nv(N + 1) × nv(N + 1).

v *y = χ y H13·Δ*v

⎡ 0 0 χ H13 ⎤ ⎢ ⎥ 0 ⎥ H = ⎢ χI 0 ⎢ ⎥ 0 ⎦ ⎣ H31 H32

where Δ*v represents the solution obtained from μΔ(H) when χ = χy. The dot symbol in (15) denotes an element by element multiplication. The column vector vy* of length nv(N + 1) includes the actual critical realizations in the nv disturbances that produce the worst-case variability for the output y using the process design and MPC configuration specified by ηk. For example, the first N + 1 elements in vy* correspond to the critical time-discrete realizations in the first disturbance, the elements from N + 2 up until 2(N + 1) represent the critical realizations of the second disturbance, and so on. Thus, vy* can be partitioned as follows:

H13 = [δv1Π , ..., δvi Π , ...δvnv Π]T H31 = [Σ1, ..., Σi , ..., Σnv] H32 = [Ο1, ..., Οi , ..., Οnv ] Οi = [δs0, i , ..., δsj , i , ..., δsN , i] Σi = [ s0,̅ i , ..., sj̅ , i , ..., sN̅ , i]

v *y = [v *y ,1, ..., v *y , i, ... v *y , nv]T

(13)

The structure of the perturbation block Δ in the SSV formulation presented in (12) is as follows:

Δv ∈ 9 nv(N + 1)

Δs = {Δs ||Δs | ≤ 1};

Δs ∈ 9 nv(N + 1)

|δp| ≤ 1;

δp ∈ C

(16)

where vy,i * represents the critical time-discrete realizations in the ith disturbance for the output y. Accordingly, the actual largest (worst-case) variability in y (ywc) due to the critical realizations in the disturbances (vy*) can be obtained from simulations of the actual closed-loop nonlinear dynamic process model, i.e., the f l’s and the MPC strategy shown in (10). The scalar ywc is obtained using a family of uncertain impulse response models (yimp); each of the uncertain models, which describe the closed-loop dynamics between an output y and a single disturbance, were obtained from a single-input−single-output (SISO) identification test. When multiple disturbances are considered in the analysis, a multiple-input−single-output (MISO) identification involving the nv disturbances and a single output y need be performed to obtain a single uncertain model (yimp) that captures the simultaneous effect of the nv disturbances acting on y. Thus, the computation of the critical disturbance ywc is performed through the use of a MISO uncertain model. While that approach will explicitly account for the simultaneous occurrence of most critical time-varying realizations in the disturbances, the implementation of that approach may return conservative designs. Also, the computational demands for the present method will increase since the identification of MISO uncertain models, and the subsequent computation of ywc is more intensive and, thus, requires an additional computational effort, than that needed for the

⎡ Δv 0 0 ⎤ ⎢ ⎥ Δ = ⎢ 0 Δs 0 ⎥ ⎢ ⎥ ⎢⎣ 0 0 δp ⎥⎦ Δv = {Δv ||Δv | ≤ 1};

(15)

(14)

where Δv and Δs are norm-bounded real scalars that are elements of the vectors Δv and Δs, respectively. The vector Δv represents the uncertainty associated with the disturbances in deviation form (vi̅ ) that generates the largest variability in y whereas Δs represents the uncertainty with each of the impulse response coefficients (δsj,i’s). Moreover, δp is a norm-bounded complex scalar that measures the system’s dynamic performance. The elements in the perturbation block Δ represent the unknowns in the structured singular value problem.38−41 Standard commercial packages such as MATLAB42 and CVX43 include efficient numerical subroutines that compute μΔ(H). 4821



Article

(controlled) variable, i.e., ywc, due to the most critical realizations in one of the disturbances included in v*y . Thus, vy* in (18) also represents the time-discrete disturbance vector that generates the worst-case deviation for that process constraint, i.e., vg*κ . Note that the constraint function gwc κ is evaluated around the process flowsheet and MPC scheme configuration specified by ηk. (17) may represent a maximum temperature constraint on a process unit (see (33) in section 3). Dynamic path constraints involving multiple process timevarying variables can be posed as follows:

identification of SISO uncertain models and the corresponding calculation of ywc. The approach proposed here relaxes that condition and assumes that the most critical time-varying realizations in the disturbances are calculated for each disturbance individually. Hence, ywc is estimated under the assumption that only one disturbance, e.g., the ith disturbance, with critical realizations given by vy,i * , is affecting the system while the rest of the disturbances remain at their nominal values. Note that the critical realizations in the disturbances obtained by this approach will not necessarily be the same than those obtained by considering the simultaneous effect of the disturbances. On the basis of the above, the procedure proposed here to calculate ywc is as follows: 1. For each process disturbance i (i = 1, ..., nv): (i) Assign the time-discrete realizations vy,i * to the ith disturbance, i.e., vi in (7), while maintaining the rest of the disturbances at their corresponding nominal values. (ii) Compute the response of output y by simulating the closed-loop dynamic behavior of the system using the critical realizations in the ith disturbance (v*y,i) assigned to the time-discrete vector vi as an input. The closed-loop model must be simulated using the process flowsheet alternative and the MPC scheme obtained from the dynamic flexibility analysis at the kth iteration step ηk (see section 2.2). 2. The above procedure is repeated up until the response of output y to the critical realizations in the nv disturbances included in v*y have been collected from simulations. 3. The output data is then used to search for the largest deviation in the output y, i.e., ywc. Estimates for ywc are used to evaluate process feasibility constraints that are a direct function of y, e.g., gκ (ηk , y(t ), V(t )) ≤ 0

gκ (ηk , y(t ), w(t ), z(t ), V(t )) ≤ 0

where y(t) ∈ y(t), w(t) ∈ y(t), and z(t) ∈ y(t). To evaluate (19), the actual nonlinear dynamic constraint function that involves the multiple process outputs (gκ) needs to considered as an output function in the analysis. That is, uncertain finite impulse response models such as (8), which describe the closed-loop dynamic behavior between each disturbance and the function gκ in (19), needs to be obtained following the systems identification procedure described in the Supporting Information. As shown in (13) and (14), sj̅ ,i and δsj,i are the key entries to compute an analytical bound on the worst-case deviation expected gκ by solving the SSV problem shown in (12). The bound obtained for the constraint function gκ (χgκ) and the norm-bounded elements in Δ*v , calculated from μΔ(H) when χ = χgκ, are the key inputs used in (15) to compute the critical realizations on the disturbances that produce the worstcase scenario for gκ, i.e., vg*κ . The procedure to compute the actual worst-case deviation expected in the constraint gκ due to wc and has been described vg*κ , i.e., gwc κ , is the same used for y above. On the basis of the above, the proposed robust feasibility test evaluates the compliance of each of the feasibility constraints by identifying, for each constraint gκ, the critical realizations in the disturbances (v*gκ ) that generate the worst-case deviation on each constraint gκ, i.e., gwc κ . This test is formulated as follows:

(17)

where gκ represents a dynamic path feasibility constraint that only depends on one time-varying process variable, e.g., a state variable x(t) a manipulated variable u(t) or an output (controlled) variable y(t). To evaluate (17), the worst-case deviation expected for that constraint will be a direct function of the worst-case variability expected for the output y, i.e., ywc. Hence, a family of uncertain finite impulse models (yimp), each describing the closed-loop dynamic behavior between each disturbance and the output y, need to be identified following the systems identification procedure described above. The estimates of the uncertain model coefficients, i.e., sj̅ ,i and δsj,i, are used to build the interconnection matrix H (13) and the perturbation block Δ (14), which are the key matrices needed to solve the SSV problem shown in (12). The solution of the SSV problem returns an analytical bound on the output y (χy) and the normalized realizations in the disturbances that produces the worst-case scenario Δν*, which are the key entries in (15) to compute the actual critical realizations in the disturbances (v*y ). The actual worst-case variability in the output y due to v*y (ywc) is obtained from the procedure described above. Hence, (17) can be reformulated as follows: gκwc(ηk , y wc , v *y ) ≤ 0

(19)

ϑ* = max ϑ ϑ∈ϑ

s.t.

ϑκ = gκwc

κ∈Κ

(20)

The main difference between gwc κ in (20) and the feasibility constraints gκ shown in (5) is that the first is evaluated using the actual critical time-discrete trajectories in the disturbances identified for the κth constraint using the current process and MPC scheme designs specified at the kth iteration step (ηk), i.e., vg*κ . On the other hand, gκ in (5) is evaluated using the block-matrix Vk, which includes the time-discrete disturbance realizations that were identified for previous process and MPC scheme designs, e.g., η1, ..., ηk−1. Hence, gwc κ and gκ are evaluated using vg*κ and Vk, respectively. The real scalar ϑ* is an index used to determine the compliance of the process feasibility constraints. As shown in Figure 1, a positive value in ϑ* means that, for the solution obtained at the kth iteration step (ηk), there are critical disturbance realizations that generate a violation in one or a few of the process constraints. Thus, the feasibility test (20) is rejected, and a new process flowsheet and MPC design alternative must be sought from the dynamic flexibility analysis shown in (5) using the previous and the current critical realizations in the disturbances. The critical time-discrete disturbance trajectories that produce constraint

(18)

gwc κ

where represents a dynamic path constraint function that depends on single time-varying output y and that is evaluated using the worst-case variability expected for that output 4822



Article

violations at the current kth iteration step are identified as follows:

A ss = A̅ ss ± δ A ss A̅ ss = {aα̅ , β }

v* = [v 1*, ..., v *i , ..., v *nv]T v *i = {v *g , i|ϑκ > 0} κ

δ A ss = {δaα , β } ∀ α = 1, ..., J ;

(21)

(22)

where Vk+1 represents the block matrix that includes the previous and the current time-discrete realizations in the nv disturbances identified for the first k process flowsheet and MPC scheme configurations obtained from the dynamic flexibility formulation shown in (5). Note that Vk+1 is the same as Vk but augmented with a column vector that includes the critical time-discrete disturbance realizations identified from the present robust feasibility analysis at the kth iteration step. The disturbance description shown in (22) ensures that the dynamic flexibility analysis can account for the previous and the most current critical scenarios identified for the disturbances. Also, this disturbance description avoids the reselection of previous solutions that were found as optimal by previous dynamic flexibility analyses. Following (20), a nonpositive value in ϑ* indicates that the feasibility constraints are satisfied in the presence of the critical realizations in the disturbances (v*), i.e., the robust dynamic feasibility test is accepted. The present robust dynamic feasibility analysis ensures that ηk satisfies the feasibility constraints in the presence of the critical realizations in the disturbances. 2.4. Stability Analyses. The present iterative decomposition algorithm implements two stability tests that aim to evaluate the nominal and asymptotic stability of those process flowsheet and MPC designs that satisfy the robust dynamic feasibility test. The stability tests included in the present analysis requires the identification of an uncertain state space model, which captures the transient characteristics of the actual system in closed-loop. The uncertain state space model is as follows:

ψ* < 0 ψ * = max(Re(ψ )) ψ∈ ψ

(25)

where the symbol Re(ψ) denotes the real parts of the eigenvalues of matrix A̅ ss, which are included in ψ and are calculated as follows:31 ψ = {ψ |det(ψ Iss − A̅ ss) = 0}

(26)

where the symbol det denotes the evaluation of the determinant of the matrix (ψIss − A̅ ss); Iss is an identity matrix of the same dimensions as matrix A̅ ss. The nominal stability test shown in (25) states that process flowsheet and MPC design alternative specified by ηk is nominally stable if and only if all the eigenvalues of A̅ ss have negative real parts. If the nominal stability condition (25) is satisfied, then a robust stability test can be performed on ηk to evaluate the system’s asymptotic stability. The linear model structure posed by the uncertain model (23) enables the application of robust stability tests that have been widely used in robust control theory and that are based on the direct method developed by Lyapunov.44 Lyapunov’s direct method requires the specification of an energy-like scalar function (Ω) such that it approaches zero along every nonzero trajectory of the states of the system ζ. Quadratic Lyapunov functions have been widely used for asymptotic stability

ζ (̇ ηk) = A ssζ (ηk) + BssV ycl̅ = Cssζ (ηk)

(24)

where the matrix A̅ ss contains the nominal model parameter estimates for the uncertain Ass state space matrix, i.e., aα̅ ,β represents the nominal value in the (α, β) element in Ass. Similarly, the matrix δAss includes the uncertain parameters associated with each of the elements in Ass, i.e., δaα,β denotes the uncertainty associated with the (α, β) element in Ass. The matrices Bss and Css pose the same structure to that shown in (24) for Ass and are not shown here for brevity. As in the robust feasibility analysis, the uncertainty in the state space model parameters account for the process nonlinearities due to changes in the disturbances that cannot be accurately described with the linear model parameters, e.g., aα̅ ,β’s. Estimates for the nominal model parameters, e.g., aα̅ ,β, and the uncertainty associated with each model parameters, e.g., δaα,β, can be obtained following the systems identification procedure described in the Supporting Information. Note that state space model shown in (23) is different from that used in the MPC algorithm shown in (1) because the first represents the response of output y in closed-loop (yc̅ l) whereas the latter captures the dynamic behavior of the system’s outputs in openloop (ŷ). Also, the only inputs considered in the uncertain state space model (23) are the magnitude-bounded disturbances. The output variable in the state space model (yc̅ l) can be a process output variable (y) or a manipulated variable (u). Furthermore, note that the uncertain state space model (23) has a linear model structure. The nominal stability test embedded in the present method requires the evaluation of the eigenvalues of the nominal model parameters estimates of the state space matrix Ass, i.e., A̅ ss in (24). The nominal stability condition is as follows:31

where vi* is a column vector that includes all the critical timediscrete realizations for the ith disturbance that produce a constraint violation for the system specified by ηk. v*i is of length nc(N + 1) where nc represents the number of process constraints that do not meet the constraint feasibility specifications for the current solution ηk; hence, v* is a column vector of length nvnc(N + 1). Accordingly, the time-discrete disturbance realizations that need to be used as inputs for the (k + 1)th dynamic flexibility analysis is as follows: Vk + 1 = [Vk , v*]

∀ β = 1, ..., J

(23)

where ζ(ηk) is a vector that represents the differential states of the uncertain state space model obtained around a nominal point specified by ηk. Likewise, yc̅ l represents the closed-loop response of an output (controlled) variable (in deviation form) evaluated around ηk. The matrices Ass, Bss, and Css are the corresponding state space matrices of the closed-loop system identified around ηk. Each of the elements included in the state space matrices is assumed to be given by a nominal value complemented with model parameter uncertainty. For example, the matrix Ass can be represented as follows: 4823



Article

particular form of the Lyapunov function, i.e., a quadratic energy-like function Ω as shown in (27). Thus, (31) is not a necessary condition for the robust stability of the system since stability conditions based on other forms of the Lyapunov function may return a different outcome. The stability condition (31) represents a convex set of LMI’s that returns a single (global) solution to this problem: whether the uncertain state space model (23) is asymptotically stable or not. The set of LMI’s posed in (31) can be efficiently solved using software such as MATLAB or CVX. Although the robust asymptotic stability test shown in (31) also ensures nominal stability, the evaluation of condition (31) may be intensive for systems that include a relatively large number of states whereas the evaluation of nominal stability condition (25) can be rapidly performed for large-scale systems. The addition of the nominal stability analysis in the present algorithm enables a quick initial screening on ηk to determine if it is locally stable or not. Thus, the present algorithm avoids the implementation of a computationally intensive robust asymptotic stability test to process flowsheet and MPC design configurations (ηk) that are locally unstable. In brief, the nominal stability has been included to reduce the computational costs in the present algorithm. In the event that ηk does not satisfy the nominal stability test (25) or the robust asymptotic stability test (31), then the index ι is updated, i.e., ι = ι + 1, and the block matrix Φ described in (3) is updated as follows:

analysis due to its simplicity and the convex characteristics of that function.45−47 Thus, the energy-like (quadratic) Lyapunov function used in the present robust stability analysis is as follows: Ω = ζ TΣζ

(27)

where Σ, the Lyapunov matrix, is considered to be independent of the states ζ included in the uncertain model (23). Accordingly, the robust asymptotic stability criterion can be posed as follows:44−47 Ω̇ < 0 Σ > 0;

Σ = ΣT

(28)

The asymptotic stability condition indicates that the uncertain model (23) is asymptotically stable if and only if there exists a symmetric positive definitive matrix Σ that is independent of ζ and that makes the energy-like function Ω to approach zero for every nonzero trajectory in the system’s states ζ. (28) ensures that all the possible trajectories in ζ are bounded and converge to zero as time approaches infinity. (28) can be reformulated using the definition (27) as follows:47 A ss TΣ + ΣA ss < 0 Σ > 0;

Σ = ΣT

(29)

The first inequality in (27) represents an infinite set of linear matrix inequalities (LMIs) since each of the elements in Ass can take an infinite number of values that satisfy the following description:

Φ = {Φ(1), ..., Φ(ii), ..., Φ(ι)} Φ(ι) = ηk

A ss = {aα , β |aα , β = [aα̅ , β − δaα , β , aα̅ , β + δaα , β ]} ∀ α = 1, ..., J ;

∀ β = 1, ..., J

(30)

A ss(ρjj )T Σ + ΣA ss(ρjj ) < 0 JJ = 2J

ρjj ∈ {aα̅ , β − δaα , β , aα̅ , β + δaα , β } ∀ β = 1, ..., J Σ > 0;

(32)

This procedure ensures that solutions, previously obtained from the dynamic flexibility formulation (5) as optimal, cannot be reselected in the next iterations since they did not satisfy the stability tests included in the present method. If the robust stability condition (31) is satisfied, then the algorithm is terminated. This means that an optimal process flowsheet and MPC design configuration that remains feasible and asymptotically stable in the presence of critical magnitude-bounded timevarying realizations in the disturbances has been identified, i.e., η* = ηk. 2.5. Summary. The methodology presented in this section can be summarized as follows (see Figure 1): 1. Set k = 0, ι = 0, Φ = []. Define V0, u, y, and η0. At each iteration k: 2. Given Vk and Φ, solve (5), the dynamic flexibility analysis. The solution to this problem returns an optimal process synthesis and MPC design alternative, ηk. 3. Implement the robust dynamic feasibility analysis on ηk. If (20) is satisfied, then proceed to step 4; otherwise, compute Vk+1 as shown in (22), set k = k+1, and go back to step 2. 4. Evaluate the nominal stability condition (25) on ηk. If nominally stability is guaranteed, then go to step 5; otherwise, update ι and Φ as shown in (32), compute Vk+1 as shown in (22), set k = k + 1, and go back to step 2. 5. Evaluate the asymptotic stability of ηk using the robust asymptotic stability condition shown in (31). If robust asymptotic stability is not guaranteed, then update ι and Φ as shown in (32), update Vk+1 as shown in (22), set k = k + 1, and go back to step 2. Otherwise, STOP; a

That is, each of the elements aα,β in Ass is bounded by extreme values. Due to the convex properties of the quadratic Lyapunov function (Ω), only the evaluation of those LMIs that involve the different combinations between the extreme values of each aα,β included in Ass is needed to ensure robust asymptotic stability.45−47 The set of combinations formed between the extreme values of each aα,β in Ass specifies a hyper-rectangle (ρ); each corner or vertex of ρ specifies a particular arrangement between extreme values of each of the model parameters aα,β included in Ass. This simplification reduces the robust asymptotic stability test presented in (29) into a finite set of LMIs, i.e.,

ρ = {ρ1 , ..., ρjj , ..., ρJJ };

∀ ii = 1, ... ι

2

∀ α = 1, ..., J ;

Σ = ΣT (31)

The formulation shown in (31) represents the robust asymptotic stability test included in the present methodology to evaluate the asymptotic stability of the system specified by ηk due to magnitude-bounded fluctuations in the disturbances. Note that (31) is only a sufficient condition to ensure the asymptotic stability of the system since it is based on a 4824



Article

Figure 2. Process flowsheet: systems of CSTRs in series with multiple inlet streams.

process flowsheet and MPC design configuration that can accommodate the critical realizations in the process disturbances at minimum costs has been found, i.e., η* = ηk. As mentioned in the Introduction, the only methodology that has considered a simultaneous process synthesis analysis and an MPC-based control strategy is that proposed by Sakizlis and co-workers.6,14 The key difference between the present method and that described in refs 6 and 14 is that the first enables the evaluation of the dynamic feasibility and asymptotic stability of the system using convex problems whereas the latter method requires the solution of an MINLP problem, which can become computationally challenging for large-scale complex systems. While the present method defines disturbances as random time-varying realizations bounded by lower and upper limits, see definition (7), Sakizlis et al.’s method assumes that disturbances are given by a user-defined time-dependent function with uncertain (critical) model parameters, e.g., a sinusoidal function with unknown frequency. Thus, Sakizlis et al.’s method6,14 is limited since there is no guarantee that other realizations in the disturbances that do not follow the dynamics imposed by the uncertain disturbance function may produce violations to the process feasibility or stability constraints. Accordingly, the method proposed in this work is more robust than Sakizlis et al.’s method because it enables the specification of an optimal process and MPC design that can accommodate critical time-discrete realizations in the disturbances. The implementation of the methodology proposed in this work is illustrated next.

3. CASE STUDY: SYSTEM OF CSTRS IN SERIES WITH MULTIPLE INLET STREAMS The iterative decomposition algorithm shown in Figure 1 has been applied to assess the optimal process flowsheet and MPC design of a system of continuously stirred tank reactors (CSTR’s) connected in series with multiple inlet streams. As shown in Figure 2, the process flowsheet superstructure under analysis involves multiple reactors and inlet feedstreams, respectively. The present analysis assumes that an exothermic irreversible reaction that transforms reactant A into product B takes place on each CSTR, i.e., A → B.The outlet flowrate of each reactor and the cooling heat flow supplied to each unit represent the process manipulated variables u(t) whereas the volume hold-up, the temperature, and the concentration of reactant A inside each reactor are variables that can be measured online and represent the controlled (output) variables y(t). The temperature and flowrate of the inlet feedstreams are subject to random time-varying fluctuations due to sudden (or sustained) variations in their corresponding upstream processes. Therefore, these input variables represent the process disturbances V(t). The concentration of each inlet stream was assumed to be constant and equal to their corresponding nominal values. Table 1 presents the process synthesis superstructure model specified for this case study. Table 2 shows the case study’s process variables, their corresponding match with those used in the methodology section, and the inlet feedstreams specifications. As shown in Table 1, the binary variable ωin,rn determines the existence of the inth inlet feed stream entering the rnth reactor. Likewise, the binary variable Rrn determines the existence of the rnth 4825



Article

the inlet streams and the reactor number, respectively. To ensure process flowsheet feasibility, the existence of the first reactor has been specified a priori by setting R1 = 1. Both ωin,rn and Rrn determine the process flowsheet structure; hence, θsys includes the IN × RN binary variables specified for ωin,rn and the 1 × RN − 1 binary variables specified for Rrn, respectively. The first logical constraint shown in Table 1 ensures that the inth inlet feed stream only enters one reactor whereas the last logical constraint ensures that the rnth reactor exists only if the (rn − 1)th reactor exists. The present case study assumes that the concentration of reactant A needs to be below 0.03 mol/L for the product’s stream at any time t whereas the reactors’ temperature cannot go beyond or below the maximum and the minimum allowed working temperatures at any time t, i.e., 500 and 400 K, respectively. These process feasibility constraints are formulated as follows:

Table 1. Process Design Superstructure for the System of CSTRs in Series first reactor (rn = 1): IN

̇ rn = Vol

∑ ωin,rnqF ,in − qrn in = 1

IN

Trṅ =

∑

ωin,rnqF ,in(TF ,in − Trn) Vol rn

in = 1 IN

̇ = CA,rn

∑

+

ΔHk 0CA,rn exp(−E / RT1)

ωin,rnqF ,in(C AF ,in − CA,rn) Vol rn

in = 1

ϑCp

− Q c ,rn

− k 0CA,rn exp(−E / RTrn)

last reactor (rn = RN): IN

̇ rn = R rn{ ∑ ωin,rnq + q Vol − qrn} F ,in rn − 1 in = 1 IN ω ⎧ q (Trn − 1 − Trn) ⎪ in,rnqF ,in(TF ,in − Trn) Trṅ = R rn ⎨ ∑ + rn − 1 ⎪ Vol Vol rn rn ⎩ in = 1

+

̇ CA,rn

ΔHk 0CA,rn exp(−E / RTrn) ϑCp

g≤0

⎫ ⎪ − Q c ,rn ⎬ ⎪ ⎭

g = {g Tinc, g Tdec, gC } A

g Tinc = {g Tinc,rn |g Tinc,rn = Trninc − 500};

IN ω q ⎧ q (CA,rn − 1 − CA,rn) ⎪ in F ,in(C AF ,in − CA,rn) ∑ = R rn ⎨ + rn − 1 ⎪ Vol rn Vol rn ⎩ in = 1

∀ rn = 1, ..., RN g Tdec = {g Tdec,rn |g Tdec,rn = 400 − Trndec};

⎫ ⎬ ⎪ ⎭

(−E / RTrn)⎪

− k 0CA,rn exp

∀ rn = 1, ..., RN gC = C AincPS − 0.03

RN

∑ ωin,rnR rn = 1

g∈g

∀ in = 1, ..., IN

(33)

A

rn = 1

R rn ≤ R rn − 1

dec where Tinc rn and Trn represent the largest variability of the temperature inside the rnth reactor in the positive (worst-case increase) and in the negative (worst-case decrease) direction, respectively. Similarly, Cinc APS represents the maximum variability of concentration of reactant A in the product stream. These constraints are included in the case study’s dynamic flexibility analysis and are estimated from the simultaneous solution of the process model equations and the MPC algorithm using as inputs the critical realizations imposed on the temperature and the flowrate of the inlet feedstreams at the kth iteration (Vk). The annualized cost function considered for the present case study (Θsys) is represented as the addition of the capital costs (CC) and the variability costs of the system in closed-loop (VC). The annualized capital cost function (34) is estimated from the capacities of each reactor included in the process flowsheet configuration. The capacity of the rnth reactor is obtained from the maximal variability in the volume hold-up observed on that reactor (Volmax rn ) due to the critical realizations in the disturbances included in the block matrix Vk. Each reactor was assumed to be an atmospheric vertical cylindrical vessel made of carbon steel with design parameters p1 (1.066) and p2 (0.802), respectively.48 Also, the cost of this process was updated to year 2011 using the Marshall and Swift cost index (MS = 1597.7).49 In addition, an annual return on investment (roi) of 20% was assumed for this process.

∀ rn = 2, ..., RN

rn = 2, ..., RN;

ωin,rn = {0, 1}IN × RN ;

R rn = {0, 1}1 × RN − 1

Table 2. Variable Categories and Inlet Stream Specifications: System of CSTRs in Series variable type process flowsheet configuration (θ): states of the system (x): manipulated variables (u): controlled variables (y): process disturbances (V): MPC tuning parameters (λ): constant model parameters:

inlet feedstreams specifications inlet stream 1 (in = 1)

inlet stream 2 (in = 2)

process variables ωin,rn, Rrn Volrn, Trn, CA,rn qrn, Qc,rn Volrn, Trn, CA,rn TF,in, qF,in Ψy, Λu E (activation energy) = 83145 J/mol k0 (pre-exponential factor) = 7.2 ×10−10 R (gas constant) = 8.3144 J/mol·K Cp (liquid’s heat capacity) = 0.239 J/g·K ΔH0 (heat of reaction) = 4.78 × 104 J/ mol ϑ (fluid’s density) = 1 × 103 kg/m3 process variables vi,nom δvi TF,1 (K) qF,1 (L/min) CA,F,1 (mol/min) TF,2 (K) qF,2 (L/min) CA,F,2 (mol/min)

350 200 1 320 250 2

10 10 10 10

⎛ MS ⎞ max ⎟(Vol + CC = 60.22roi⎜ 1 ⎝ 280 ⎠

RN

∑

(p1 + p2 )/3 R rn Vol max rn )

rn = 2

(34)

reactor in the process flowsheet. For example, setting ω2,3 = 1 means that the second inlet stream is fed to the third reactor whereas R3 = 1 indicates that at least 3 CSTRs are included in the process flowsheet structure. The indexes in and rn denotes

The process needs to continuously produce a stream that maintains a molar flowrate of product B at its target value (wtar B ). This specification can be used to measure the system’s 4826



Article

dynamic performance. Thus, the annualized variability costs (VC) for this process are calculated as follows: var VC = 5.256 × 105ϕBwB,PS

Table 3. Performance of the Optimal Process and Control Design Algorithm for the System of CSTRs

(35)

iterations

wvar B,PS

where represents the maximal variability in the molar flowrate of component B at the outlet stream of the last reactor with respect to its nominal (target) specification, i.e., var inc dec wB,PS = max{wB,PS − wBtar , wBtar − wB,PS }

winc B,PS

1

(36)

2

wdec B,PS

where and represent the largest variability in the positive direction (worst-case increase) and in the negative direction (worst-case decrease) of component B at the product dec stream (PS), respectively. Estimates for winc B,PS and wB,PS can be obtained from the simultaneous solution of the process model equations (see Table 1) and the MPC algorithm for specific realizations in the disturbances. The product stream (PS) is as follows:

3

a

RN

PS = 1 +

∑ rn = 2

R rn

(37)

process flowsheeta 2 CSTRs: IS(1) enters CSTR 2 IS(2) enters CSTR 1 3 CSTRs: IS(1) enters CSTR 2 IS(2) enters CSTR 1 3 CSTRs: IS(1) enters CSTR 2 IS(2) enters CSTR 1

Θsys ($/y)

robust feasibility test

stability tests

6710

failed

not performed

7176

failed

not performed

7534

satisfied

satisfied

IS(i): ith inlet stream.

the presence of critical realizations in the disturbances. Note that the stability tests were satisfied for the first process design and MPC scheme configuration that satisfied the robust feasibility analysis, i.e., the matrix Φ shown in (3) and in Figure 1 remained empty during the calculations. For the present case study, V0 was initialized with a series of step changes using the disturbances’ upper and lower limits. Thus, a larger number of iterations to that reported in Table 3 would have been needed if V0 would have been initialized using the disturbances’ nominal (steady-state) values. Table 4 shows the optimal process and MPC design obtained with the present method (optimal design, MPC-based). A process flowsheet with three CSTRs with the first and the second inlet streams entering the second CSTR and the first CSTR, respectively, was specified for this case study. The results show that the optimal set-point temperature for the second reactor has been specified close to its maximum allowed working temperature whereas the temperature of the other two CSTRs were specified farther from the maximum temperature operating constraint. When simulating the optimal process design in open-loop, it was found that the nominal operating point of this system of CSTR’s is a stable steady-state operating point. However, a thermal runaway was observed in the open-loop system when step changes in the temperature of inlet feed stream entering the first reactor is simulated with the present design. This result shows the need to incorporate a control algorithm such as an MPC controller to avoid this undesirable condition in the system. The results also show that the variability costs (VC) are dominant for this process since it represents almost 70% of the total annualized costs. Figure 3 shows the realizations in the inlet feedstream conditions (the disturbances) obtained at solution point, i.e., Vk. As shown in that figure, the critical realizations for the disturbances include large deviations with high and low frequency content. The results obtained by the present analysis were validated by simulating the optimal process flowsheet and MPC design using the critical realizations in the disturbances shown in Figure 3. Figure 4 shows that the temperature in the first and in the second CSTR remained stable and at their corresponding maximum allowed operating limit (500 K) in the presence of the worst-case increase in temperature due to the most critical realizations in the temperature and flowrate of the inlet feedstreams. The rest of the process operational constraints

As shown in (35), the present variability cost function (VC) assigns a cost (φB = 0.001$/mol of B) to the largest deviations (worst-case variability) of component B’s molar flowrate at the product stream with respect to its target specification (wtar B ). That is, the variability cost function penalizes the largest variability in any direction in B’s molar flowrate due to critical realizations in the disturbances. While a worst-case decrease in B’s flowrate at PS denotes an economic loss since it does not comply with the process goals, a large variability in B’s flowrate in the positive direction represents a loss in profit since the process is producing a product with a quality that is higher than that required. Hence, the variability cost function VC aims to weight (in economic terms) the dynamic variability of product B in closed-loop due to changes in the disturbances. The goal of this case study is to identify the process flowsheet and the MPC scheme that maintains the process asymptotically stable and close to its design, controllability, and operability specifications in the presence of critical (timevarying) realizations in the disturbances at minimum cost. The simultaneous design and control of this case study is challenging because: (i) It involves multiple structural decisions since it needs to specify the number of CSTRs and the corresponding feedstream(s) that enter each reactor, respectively. (ii) The process is subject to multiple process disturbances since the temperature and flowrate of each inlet stream are considered as random (time-varying) independent disturbances. (iii) A linear constrained model predictive controller with several manipulated and controlled variables needs to be properly tuned to maintain the feasible and stable operation of the process despite critical (magnitudebounded) fluctuations in the disturbances. The methodology presented in the previous section was applied to this case study; the step-by-step implementation of the proposed methodology is presented in the Supporting Information. Table 3 summarizes the evolution of the optimal process flowsheet and MPC control design methodology when applied to the present case study. As shown in that table, three iterations were required to achieve a process flowsheet and MPC scheme that remains feasible and asymptotically stable in 4827



Article

that reported by the present method (optimal design, MPCbased). Also, the process flowsheet specified by the steady-state design is slightly different to that obtained by the present method since it specifies that the two inlet streams enter the first CSTR. However, the solution returned by the steady-state design is not dynamically feasible since constraint violations were observed when that design is simulated using time-varying changes in the disturbances, e.g., a step change in the disturbances. This result shows the need to account for a controllability analysis at the earlier stages in process synthesis. In order to compare the benefits between the implementation of a multivariable (model-based) control scheme such as MPC and a decentralized control scheme, the method presented in this work was applied for the simultaneous process synthesis and control design of the system of CSTRs using a multiloop PI control superstructure. To evaluate this scenario, the MPC algorithm shown in (1) was replaced by a general multiloop PI control superstructure.21 Consequently, the tuning parameters of the PI controllers included in the decentralized control superstructure used for the present case study were considered as optimization variables in that simultaneous design and control formulation. Although the implementation of a multiloop control scheme avoids the need to solve a convex optimization problem at each time step as in the case of a linear MPC, the decentralized control strategy adds new binary decision variables in the dynamic flexibility analysis formulation for the specification of the pairings between the set of controlled variables y and manipulated variables u, which increases the dimensionality of the MINLP formulation in the dynamic flexibility analysis. As shown in Table 4 (optimal design, PI-based), the implementation of a decentralized control structure within the simultaneous design and control method returns a process flowsheet that is slightly different to that obtained when the control scheme is based on a linear constrained MPC, i.e., the first and the second inlet streams enter the first and the second CSTRs, respectively. The solution obtained by the multiloop PI control structure remains feasible and stable in the presence of any realization in the disturbances since it was obtained using the iterative algorithm presented in this work. However, the total annual cost (Θsys) obtained for the decentralized control scheme is 8.2% higher than that obtained by the MPC algorithm (optimal design, MPC-based). As mentioned above, the variability costs dominate the process economics for this case study. The results shown in Table 4 indicate that the variability costs in the decentralized control scheme are approximately 18% higher than that obtained by the optimal process flowsheet and MPC design. This result indicates that a multivariable control scheme based on a linear MPC returns an economically attractive process design for the present case study. While the MPC control algorithm is a multivariable control structure than can implement control actions while simultaneously considering the variability in all the controlled variables, the decentralized control scheme is based on local PI controllers, which actions are independent of each other. Therefore, it is expected that the MPC control scheme will provide efficient control moves to maintain the system at its near optimal operating conditions in the presence of critical realizations in the disturbances. Accordingly, the variability of the process variables in the presence of disturbances is expected to be reduced when the control scheme is a multivariable controller such as MPC. However, a multiloop PI control scheme is expected to provide similar results to MPC when the

Table 4. Optimal Designs Using Different Approaches: System of CSTRs in Series optimal design, MPCbaseda process flowsheet 3 CSTRs IS(1) enters CSTR 2

optimal SS

optimal design, PIbased

3 CSTRs IS(1) and (2) enter CSTR 1

3 CSTRs IS(1) enters CSTR 1 IS(2) enters CSTR 2

not specified

PI control loops:

IS(2) enters CSTR 1 control structure Λq1 = 94.80 Λq2 = 99.58 Λq3 = 56.34 ΛQc,1 = 70.11 ΛQc,2 = 2 × 10−4 ΛQc,3 = 99.85

operating conditions 486.80

T1* (K)

ΨVol,1 = 94.96 ΨT1 = 18.31 ΨCA,1 = 79.76 ΨVol,2 = 81.15

q1−Vol1: Kc = 80.12, τI = 70.03 q2−Vol2: Kc = 76.71, τI = 35.01 q3−Vol3: Kc = 84.31, τI = 67.33

ΨT2 = 46.47

Qc,1−T1: Kc = 80.30, τI = 70.03

ΨCA,2 = 99.69 ΨVol,3 = 78.84 ΨT3 = 99.86 ΨCA,3 = 84.02

Qc,2−T2: Kc = 88.28, τI = 80.44 Qc,3−T3: Kc = 68.97, τI = 84.40

500

486.56

* (mol/L) CA,1

0.038

0.02

0.041

Vol* 1 (L)

147.09

20.55

54.26

T2* (K)

496.24

500

492.72

* (mol/L) CA,2

0.012

0.0437

0.034

Vol*2 (L)

127.60

10.84

130.47

T3* (K)

485.68

500

493.67

* (mol/L) CA,3

0.01

0.01

0.01

Vol*3 (L)

1.38

10.21

9.48

CC ($/y)

2277.86

700.52

1831.62

VC ($/y)

5256.76

0.0

6376.05

Θsys ($/y)

7534.63

700.52

8207.67

annualized costs

a

IS(i): ith inlet stream.

were validated in the same fashion and are not shown here for brevity. These results show that the implementation of the proposed methodology returns an optimal process flowsheet and MPC design alternative that remains stable and feasible in the presence of the most critical realizations in the disturbances. To illustrate the benefits of the proposed approach, the present case study was also solved using two different design approaches: an optimal steady-state process design (optimal SS) and a simultaneous process and control design using a multiloop PI control structure (optimal design, PI-based). In the steady-state design approach, the process dynamics were neglected and the disturbances affecting the process were assumed to be fixed and equal to their nominal values. As shown in Table 4 (optimal SS), the annualized cost for the steady-state design approach is orders of magnitude lower than 4828



Article

Figure 3. Critical realizations for the temperature and flowrate of the inlet streams.

Figure 4. Worst-case variability of the temperature inside CSTR 1 and CSTR 2, respectively.

case study demonstrates the potential capabilities of a multivariable model-based control strategy such as MPC for simultaneous process flowsheet and control design.

level of interactions between the manipulated and controlled variables is relative small or insignificant. Also, the weights assigned in the cost function to each of the aspects included in the process economics are key in the selection of a control strategy that returns an optimal process design. On the basis of the above, the qualitative conclusion obtained for the present case study, i.e., a multivariable MPC is more economically attractive than that obtained by a decentralized PI control scheme, is problem-specific and cannot be generalized to other case studies. Nevertheless, the result obtained for the present

■

CONCLUSIONS A new iterative decomposition framework for the simultaneous process flowsheet and MPC design was presented. While previous methodologies evaluate the dynamic feasibility and asymptotic stability using an MINLP formulation, the key benefit of the present method is that both, the robust dynamic 4829



■

feasibility analysis and the robust asymptotic stability analysis, are formulated as convex problems that can be efficiently solved using numerical subroutines included in commercial software packages. The robust dynamic feasibility analysis estimates the largest (worst-case) variability expected for the process variables included in the process feasibility constraints under the assumption that the disturbances are random time-varying (magnitude-bounded) variables. This is a key feature introduced by the present method since previous optimization-based methodologies that included structural decisions in the analysis require the specification of a time-dependent disturbance function with unknown (critical) model parameters, which is used to describe the transient behavior of this variable. Thus, the optimal process flowsheet and MPC designs specified by the algorithm presented in this work is more general since it is not restricted to a particular (user-defined) disturbance function and can accommodate any low and high frequency magnitude-bounded realizations in the process disturbances. The implementation of the present method requires the identification of uncertain models around a nominal operating condition. While the identification of those models enables the formulation of the robust feasibility and robust stability tests as convex optimization problems, those model descriptions represent approximations to the full nonlinear closed-loop process model. Thus, the critical time-discrete disturbance realizations identified by the present method may be deviate from those obtained with the full nonlinear closed-loop process model. Hence, the implementation of the critical disturbance trajectories computed from the nonlinear dynamic model on the process design specified by the present method may result in constraint violations. Although the computation of the critical disturbance trajectories using uncertain impulse response models adds a computational cost to the method, which may become significant for complex systems with a large number of constraints, the direct identification of the critical disturbance trajectories using the full nonlinear process model is a computationally intensive task even for simple processes. Thus, the implementation of the latter approach to compute the critical disturbance realizations may become prohibitive for systems with a large number of constraints. Moreover, the present method may return conservative process design configurations due to the use of the uncertain process models. However, the uncertainty in the models is captured in closedloop; thus, it is expected that the present methodology may not specify overly conservative designs though this particular aspect of the method depends on the level of complexity and nonlinearity of the dynamic system under analysis. The results obtained from the case study suggest that the proposed MPC-based methodology is a practical tool that has potential to address the optimal process synthesis and control design of industrial processes.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 519-888-4567, x38667. Fax: 519-888-4347. Notes

The authors declare no competing financial interest.

■

ASSOCIATED CONTENT

S Supporting Information *

Information that supports the proposed methodology and the step-by-step implementation of the proposed method into the case study are presented as supporting information. This material is available free of charge via the Internet at http:// pubs.acs.org. 4830

NOMENCLATURE A = state space matrix of the internal process model in the MPC algorithm Ass = state space matrix of the closed-loop process model around ηk aα̅ ,β = nominal values in the Ass state space matrix B = state space matrix of the internal process model in the MPC algorithm Bss = state space matrix of the closed-loop process model around ηk C = state space matrix of the internal process model in the MPC algorithm Cinc A,PS = worst-case increase of of reactant A in the product stream, system of CSTRs in series CA,rn = concentration of the reactant A in the rnth reactor, system of CSTRs in series CA,F,in = concentration of reactant A in the inth inlet feed stream, system of CSTRs in series Css = state space matrix of the closed-loop process model around ηk CC = capital cost, system of CSTRs in series d = vector of process design variables f l = process differential equations f MPC = representation of the MPC algorithm gκ = process feasibility constraint functions gwc κ = feasibility constraint functions evaluated with the critical realizations in the disturbances H = interconnection matrix in the SSV analysis Iu = identity matrix for the manipulated variables, system of CSTRs in series Iy = identity matrix for the controlled variables, system of CSTRs in series ii = auxiliary iteration index used to store unstable process and control design configurations in = inlet streams index k = iteration index in the simultaneous design and control algorithm l = index set of process differential equations n = sampling instant N = final sampling instant M = large constant number MS = Marshall and Swift cost index P = prediction horizon in the MPC scheme p1 = CSTR design parameter, system of CSTRs in series p2 = CSTRs design parameter, system of CSTRs in series PS = product stream, system of CSTRs in series Q = control horizon in the MPC scheme Qc,rn = heat cooling rate supplied to the rnth reactor, system of CSTRs in series qF,in = mass flowrate of the inth inlet feed stream, system of CSTRs in series qrn = outlet mass flowrate of the rnth reactor, system of CSTRs in series r = vector for set points dx.doi.org/10.1021/ie302215c | Ind. Eng. Chem. Res. 2013, 52, 4815−4833


Article

r̂ = vector of set points in deviation form used in the MPC formulation Rrn = binary variable that determines the existence of the rnth reactor, system of CSTRs in series rn = reactor number index, system of CSTRs in series roi = annual return on investment, system of CSTRs in series sj,i = uncertain impulse response model coefficient sj̅ ,i = nominal value of the j,ith impulse response model coefficient t = time tf = process settling time Trn = temperature of the rnth reactor, system of CSTRs in series TF,in = temperature of the inth inlet feed stream, system of CSTRs in series Tdec rn = worst-case decrease of the temperature inside the rnth reactor, system of CSTRs in series Tinc rn = worst-case increase of the temperature inside the rnth reactor, system of CSTRs in series Tssrn = vector that includes the steady-state temperatures of each reactor, system of CSTRs in series u = vector of manipulated variables ûn = vector of manipulated variables’ values in deviation form at the nth sampling instant unom = vector that includes the nominal (steady-state) values in the manipulated variables V(t) = vector of time-dependent (continuous) process disturbances V = vector of time-discrete (random) process disturbances v* = critical time-discrete realizations in the disturbances at the kth iteration step vi = time-discrete realizations in the ith disturbance vi,nom = nominal values of the ith disturbance Vk = block matrix that includes the critical realizations in the disturbances up until the (k − 1)th iteration vk−1,i = critical time-discrete realizations of the ith disturbance at the (k − 1)th iteration step V̂ n = time-discrete disturbances realizations (in deviation form) to be used by the MPC formulation Vy* = critical time-discrete realizations of the disturbances for a controlled variable y VC = variability cost, system of CSTRs in series Volrn = volume hold-up of the rnth reactor, system of CSTRs in series Volmax m = maximal variability in the volume hold-up in the rnth reactor, system of CSTRs in series Volssrn = volume hold-up of all the reactors at steady-state, system of CSTRs in series wtar B = molar flowrate of product B at its target value, system of CSTRs in series wvar B,PS = largest variability in component B’s molar flowrate at the product stream, system of CSTRs in series wdec B,PS = worst-case decrease of component B at the product stream, system of CSTRs in series winc B,PS = worst-case increase of component B at the product stream, system of CSTRs in series x = vector of time differential state variables x̂n = vector of state variables’ values of the internal MPC process model at the nth sampling instant y = vector of time-dependent process output (controlled) variables yc̅ l = closed-loop response of output y around ηk

ŷn = vector of controlled variables’ values of the internal MPC process model at the nth sampling instant yimp = vector function that represents a family of uncertain impulse response models ywc = worst-case scenario for output variable y Z = vector of smooth approximation functions Greek Letters

Δûn = vector of future moves in the manipulated variables at the nth sampling instant, MPC algorithm Δt = sampling period Δ = perturbation block in the SSV analysis Δ*v = norm-bounded vector that includes the critical timediscrete realizations in the disturbances δaα,β = uncertainty in the elements of the A state space matrix (Ass) δsj,i = uncertainty associated with the impulse response model coefficient si,j δvi = maximal deviation of the ith disturbance ε = tolerance criteria for optimization variables of type real ζ = states of the system around ηk η0 = initial guess for the set of decision variables η = set of decision variables in the simultaneous process design and control scheme ηk = solution from the dynamic flexibility analysis at the kth iteration step η* = optimal process and MPC scheme configuration Θ = objective function in the simultaneous process and control design formulation Θsys = annualized cost function, system of CSTRs in series θ = vector of integer variables associated with the topology of the process flowsheet θsys = vector of binary process design variables, system of CSTRs in series ι = iteration index that keeps track of unstable process and control design configurations κ = index set of process inequality constraints Λ = matrix that includes the weights of the manipulated variables in the MPC scheme Λsys = weights of the manipulated variables in the MPC scheme, system of CSTRs in series λ = vector of the MPC tuning parameters λsys = MPC tuning parameters, system of CSTRs in series ρ = hyper-rectangle that includes all the combinations between the uncertain elements in Ass φB = cost assigned to the worst-case variability of B at the product stream, system of CSTRs in series Φ = block matrix that stores unstable process and control design configurations Φ = subscript assigned to those optimization variables that specify an unstable process and control design configuration χy = bound on the worst-case scenario for an output variable y ϑ* = index used to determine the compliance of the robust feasibility test ψ = eigenvalues of the state space matrix Ass Ψ = matrix that includes the weights of the controlled variables in the MPC scheme Ψsys = weights of the controlled variables in the MPC scheme, system of CSTRs in series Ω = energy-like Lyapunov function ωin,rn = binary variable that determines the existence of the inth inlet feed stream entering the rnth reactor Σ = Lyapunov matrix

4831



■

Article

(22) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N. A case study in simultaneous design and control using rigorous, mixed-integer dynamic optimization models. Ind. Eng. Chem. Res. 2002, 41, 760−778. (23) 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. (24) Kookos, I. K.; Perkins, J. D. An algorithm for simultaneous process design and control. Ind. Eng. Chem. Res. 2001, 40, 4079−4088. (25) 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. (26) 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, 801−813. (27) Ricardez-Sandoval, L. A.; Budman, H. M.; Douglas, P. L. Simultaneous design and control: A new approach and comparison with existing methodologies. Ind. Eng. Chem. Res. 2010, 49, 2822− 2833. (28) Ricardez-Sandoval, L. A.; M., H.; Douglas, H. M.; Budman, H. M. A methodology for the simultaneous design and control of largescale systems under process parameter uncertainty. Comput. Chem. Eng. 2011, 35, 307−318. (29) Maciejowski, J. M. Predictive Control with Constraints; Prentice Hall: London, 2002. (30) Ljung, L. System Identification: Theory for the User, second ed.; Prentice-Hall: Upper Saddle River, NJ, 1999. (31) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A.; Doyle, F. J. Process dynamics and control; John Wiley & Sons: New Jersey, 2011. (32) Floudas, C. A.; Pardalos, P. M. Encyclopedia of optimization; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2001. (33) Bemporad, A.; Morari, M.; Lawrence, R. L. N. Model predictive control toolbox. User’s guide; The MathWorks, Inc.: Natick, MA, 2012. (34) Floudas, C. A. Nonlinear and mixed-integer optimization; Oxford University Press: New York, 1995. (35) 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. (36) Braatz, R.; Young, P. M.; Doyle, J. C.; Morari, M. Computational-complexity of mu-calculation. IEEE Trans. Autm. Cont. 1994, 39, 1000−1002. (37) 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. (38) Doyle, J. Analysis of feedback systems with structured uncertainties. IEE Proc., Part D: Control Theory Appl. 1982, 129, 242−250. (39) Morari, M.; Zafiriou, E. Robust process control; Prentice Hall: Upper Saddle River, NJ, 1989. (40) Skogestad, S.; Postlethwaite, I. Multivariable feedback control: Analysis and design, second ed.; John Wiley & Sons, Ltd.: New Jersey, 2005. (41) Zhou, K.; Doyle, J. C.; Glover, K. Robust and optimal control; Prentice Hall: Upper Saddle River, NJ, 1996. (42) Balas, G.; Chiang, R.; Packard, A.; Safonov, M. Robust Control Toolbox User’s Guide; The MathWorks Inc.: Natick, MA, 2012. (43) Grant, M. C.; Boyd, S. P. Graph implementations for non smooth convex programs. In Recent advances in learning and control; Blondel, V., Boyd, S., Kimura, H., Eds.; Lecture notes in control and information sciences; Springer: New York, 2008. (44) Lyapunov, A. M. The General Problem of motion stability; Princeton University Press: Princeton, NJ, 1949. (45) Slotine, J. J. E.; Li, W. Applied nonlinear control; Prentice-Hall: Upper Saddle River, NJ, 1991. (46) Khalil, H. K. Nonlinear systems; Prentice-Hall: Upper Saddle River, NJ, 2002. (47) Boyd, S.; Yang, Q. P. Structured and simultaneous Lyapunov functions for system stability problems. Int. J. Control 1989, 49, 2215− 2240.

REFERENCES

(1) Buckley, P. S. Techniques of process control; Wiley: New York. 1964. (2) Lenhoff, A. M.; Morari, M. Design of resilient processing plants 0.1. Process design under consideration of dynamic aspects. Chem. Eng. Sci. 1982, 37, 245−258. (3) Luyben, W. L. The need for simultaneous design education. In The integration of process design and control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam, The Netherlands, 2004; pp 10−41. (4) Seider, W. D.; Seader, J. D. Lewin, D. R.; Widagdo, S. Producto and process design principles: synthesis, analysis and design, third ed.; John Wiley & Sons: Hoboken, NJ, 2009. (5) Shinskey, F. G. Uncontrollable processes and what to do about them. Hydrocarbon Process. 1983, 62, 179−182. (6) 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. (7) Seferlis, P., Georgiadis, M. C. The integration of process design and control; Elsevier: Amsterdam, The Netherlands, 2004. (8) 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. (9) 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−46. (10) Brengel, D. D.; Seider, W. D. Coordinated design and control optimization of nonlinear processes. Comput. Chem. Eng. 1992, 16, 861−886. (11) Loeblein, C.; Perkins, J. D. Structural design for on-line process optimization: 1. Dynamic economics of MPC. AIChE J. 1999, 45, 1018−1029. (12) Soliman, M.; Swartz, C. L.E.; Baker, R. A mixed-integer formulation for back-off under constrained predictive control. Comput. Chem. Eng. 2008, 32, 2409−2419. (13) Francisco, M.; Vega, P.; Alvarez, H. Robust integrated design of process with terminal penalty model predictive controllers. Ind. Eng. Chem. Res. 2011, 89, 1011−1024. (14) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. N. Simultaneous process and control design using mixed-integer dynamic optimization and parametric programming. In The integration of process design and control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam, The Netherlands, 2004; pp 187−215. (15) Swartz, C. L. E. The use of controller parameterization in the integration of design and control. In The integration of process design and control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam, The Netherlands, 2004; pp 239−263. (16) 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. (17) Moon, J.; Kim, S.; Linninger, A. A. Integrated design and control under uncertainty: Embedded control optimization for plantwide processes. Comput. Chem. Eng. 2011, 35, 1718−1724. (18) Alhammadi, H.; Romagnoli, J. 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, The Netherlands, 2004. (19) 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. (20) Ricardez-Sandoval, L. Optimal design and control of dynamic systems under uncertainty: A probabilistic approach. Comput. Chem. Eng. 2012, 43, 91−107. (21) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal design of dynamic systems under uncertainty. AIChE J. 1996, 42, 2251−2272. 4832



Article

(48) Douglas, J. M. Conceptual design of chemical processes; Mc-Graw Hill: Columbus, OH, 1988. (49) Economic indicators. Chem. Eng. 2012, April.

4833


Simultaneous Design and Control under Uncertainty Using Model

Recommend Documents