Novel Method for Considering Process Flexibility and Stability

Article pubs.acs.org/IECR

Novel Method for Considering Process Flexibility and Stability Simultaneously Hao Jiang, Bingzhen Chen,* Hangzhou Wang, Tong Qiu, and Jinsong Zhao Department of Chemical Engineering, Tsinghua University, Beijing 100084, China S Supporting Information *

ABSTRACT: Flexibility and stability are key components of chemical process operability and are important in practical applications. In previous research, flexibility and stability have been studied separately, resulting in flexible regions containing unstable parts that should be avoided in terms of operability. This paper considers flexibility and stability simultaneously, with the aim of obtaining a stable flexible region for the process. On the basis of existing algorithms for flexibility analysis, the vertex algorithm and active constraint algorithm are presented for the simultaneous analysis of flexibility and stability. Certain limitations of these algorithms are explored and discussed. Then, using the eigenvalue optimization method, we propose a new formulation of the simultaneous analysis of process flexibility and stability. Case studies of the methyl methacrylate polymerization and 1,3-propanediol fermentation are provided. The obtained results illustrate the need to consider flexibility and stability simultaneously and demonstrate the effectiveness of the novel method presented in this paper.

1. INTRODUCTION Many uncertainties exist in chemical process design, e.g., thermodynamic and kinetic parameters, equipment performance, market forecasts, and raw material and product prices.1 The aim of chemical process design is to obtain an optimal design with respect to cost and performance while taking into account these variables. The concept of flexibility was introduced to describe the ability to operate over a range of uncertain conditions. Any chemical process can be described by the following equations and inequalities:2,3 ⎧ ⎪ h(d , z , x , θ ) = 0 ⎨ ⎪ ⎩ g (d , z , x , θ ) ≤ 0

F = max δ s. t. max min max f j (d , z , θ ) ≤ 0 θ ∈ T (δ)

z

j∈J

T (δ) = {θ|(θ N − δ Δθ −) ≤ θ ≤ (θ N + δ Δθ +)}

(4)

Figure 1 shows a graphical representation of the flexibility index. The irregular region represents the feasible region, which

(1)

Here, d is the vector of design variables that define equipment sizes and is determined at the design stage, z is the vector of control variables that represent the degrees of freedom, x is the vector of state variables having the same dimension as h, and θ represents the uncertain parameters. When the form of h is simple, for a fixed design, x can be expressed as a function of θ and z using x = x(d, z, θ), so that the chemical process can be described by the following inequalities: g (d , z , x(d , z , θ ), θ ) = f (d , z , θ ) ≤ 0

Figure 1. Schematic diagram of flexibility index.

(2)

Each of the uncertain parameters has a nominal operating point denoted by θN. Their maximum positive and negative deviations are determined at the design stage, denoted by Δθ+ and, Δθ−, respectively. A non-negative scalar δ is used to represent the degree of deviation from the nominal operating point, and the range of uncertain parameters is described by a hyperrectangle T(δ): T(δ) = {θ|(θ N − δ Δθ −) ≤ θ ≤ (θ N + δ Δθ +)}

is the allowable variation range of θ, and the flexible region corresponds to the largest inscribed rectangle T. For calculation of the flexibility index, Swaney and Grossmann proposed the heuristic search algorithm and implicit enumeration algorithm,4 Grossmann and Floudas proposed the active constraint algorithm,5 Ostrovsky et al. proposed the ULB (upper

(3)

Received: Revised: Accepted: Published:

To describe flexibility quantitatively, Swaney and Grossmann introduced the flexibility index which is the maximum value of δ defined as follows:2 © 2014 American Chemical Society

14765

April 18, 2014 August 20, 2014 August 27, 2014 September 15, 2014 dx.doi.org/10.1021/ie501602d | Ind. Eng. Chem. Res. 2014, 53, 14765−14775

Industrial & Engineering Chemistry Research

Article

and lower bound) method,6 Pistikopoulos and Ierapetritou introduced a novel method for optimal process design under uncertainty, 7 Ostrovsky et al. proposed the BB-Active algorithm,8 Raspanti, Bandoni, and Biegler proposed a method based on the KS function and smoothing function,9 and Floudas, Gumus, and Ierapetritou proposed a flexibility index calculation algorithm based on a global optimization algorithm.10 Ostrovsky et al. developed a general approach to solving the three subproblems based on the split and bound strategy,11 and extended the feasibility test and the two-stage optimization problem.12 Chemical processes may be nonlinear systems,13,14 and chemical production systems may experience a large number of disturbances. Stability is used to represent the tolerance for small perturbations in a production system and emphasizes the dynamic characteristics of a system. Besides nonlinearity, systems may have multiple steady states.15−18 A set of differential equations is used for the dynamic description of chemical processes. dx = F (x ) dt

Figure 2. Stable and unstable part of flexible region.

2. NECESSITY FOR THE SIMULTANEOUS ANALYSIS OF FLEXIBILITY AND STABILITY To illustrate the necessity for this study, a case study of the MMA polymerization process is provided where flexibility only and simultaneous flexibility and stability are studied. 2.1. Methyl Methacrylate Polymerization Reaction. The polymerization reaction takes place in a CSTR with jacketed cooling by water. Monomer and initiator are added to the top of the reactor and the polymer is removed from the bottom of the reactor. A mathematical model of the process includes six differential equations obtained from literature:20,21

(5)

When F(x) = 0, we obtain the steady-state solution of eq 5. When tiny perturbations Δx are added to x*, eq 5 becomes ⎧ dx = F (x ) ⎪ ⎨ dt ⎪ ⎩ x(0) = x * + Δx

F(Cmin − Cm) dCm = −(k p + k fm)CmP0 + V dt FC − FC I dC I = −kIC I + I Iin V dt

(6)

If the solution x(t) of eq 6 satisfies the condition limt→∞x(t) = x*, then the solution x* is stable; otherwise, x* is unstable. Stability can be determined from the eigenvalues of the Jacobian matrix J of F(x): (1) if the real parts of all eigenvalues of J are negative, the steady-state solution is stable; (2) if the eigenvalues of J have a positive real part, the steady-state solution is unstable; (3) if all the real parts of the eigenvalues of J are nonpositive and have zero real part eigenvalues, the stability of the steady-state solution cannot be judged by J. Various indicators of operability have been studied separately, and many achievements have been made in the calculation of flexibility index and stability. Speed and scope of application of the algorithms have improved significantly. In recent years, researchers have begun to combine various operability indicators, including the combination of flexibility and controllability.19 However, previous studies have not considered stability when analyzing flexibility. This may result in flexible regions containing unstable parts that should be avoided in terms of operability, as shown in Figure 2. Although the flexible region meets production requirements, control during unstable times is difficult. We therefore want to find a region that not only adapts to changes in uncertain parameters, but also ensures stable operation. For designers, this can simplify the design of control system, and for operators, this can promote the safety of process and cut down the cost of operating. This can be achieved by considering flexibility and stability simultaneously, such as by recalculating the flexible region while considering stability as a constraint. Section 2 in this paper highlights the necessity for this research work. This is followed by a study of algorithms for analyzing flexibility and stability simultaneously. In section 4, a novel eigenvalue optimization algorithm is investigated. In section 5, two cases of methyl methacrylate (MMA) polymerization and 1,3-propanediol fermentation are studied. In section 6, the scopes of application and limitations of the three algorithms are discussed. The conclusions are presented in section 7.

−ΔHk pCm F(Tin − T ) UA dT = P0 − (T − Tj) + V ρCp ρCpV dt FD0 dD0 = (0.5k tc + k td)P02 + k fmCmP0 − V dt FD1 dD1 = M m(k p + k fm)CmP0 − V dt dTj dt

=

Fcw(Tw0 − Tj) V0

+

UA(T − Tj) ρw CpwV0 (7)

Here, P0 = ((2f *CIkI)/(ktd + ktc))1/2, kr = Ar e−Er/RT, r = p, fm, I, td, tc. The values of the parameters are available in the Supporting Information. According to eq 7, the state variables are Cm, CI, T, D0, D1, and Tj. To simplify the system without loss of generality, one operating variable (cooling water flow Fcw) and two uncertain parameters (monomer feed concentration Cmin and monomer feed flow F) are considered. Given the nominal operating point, CNmin = 6.4678 kmol m−3 and FN = 1 m3 h−1, the maximum positive and negative deviations are Δθ+1 = 0.3, Δθ−1 = 0.2, Δθ+2 = 0.2, and Δθ−2 = 0.2. The range of process parameters is given in Table 1. 2.2. Feasible and Flexible Regions of MMA Polymerization. To illustrate the necessity for the simultaneous analysis of flexibility and stability, the feasible and flexible regions of MMA polymerization were studied. The feasible region R is a range within which exist some operating variables z which ensure all process indicators meet requirements when the uncertain 14766

dx.doi.org/10.1021/ie501602d | Ind. Eng. Chem. Res. 2014, 53, 14765−14775


Article

determined by the range of parameters in Table 1. φ is discretized from 0 to 2π, the optimization problem expressed by eq 10 is solved for every φ and every θ on the diagram is described with a point. The real part of the eigenvalues of the Jacobian matrix are checked for every point. If they are smaller than 0, the point is stable, otherwise it is unstable. The feasible and flexible regions are shown in Figure 4a. If the stability constraints are embedded in eq 10, we obtain the stable feasible region. The formulation is as follows:

Table 1. Range of Process Parameters process params

range

Cm/kmol m−3 CI/kmol m−3 T/K D0/kmol m−3 D1/kg m−3 Tj/K Fcw/m3 h−1

5.4141−6.0171 0−0.03 340−360 0−0.003 40−60 293.2−350 0.1−0.3

δ *(θ )̃ = max δ

parameters are changed arbitrarily. The mathematical formulation is as follows: ∀ θ ∈ R { ∃ z|f (d , z , θ ) ≤ 0}

s. t. h(d , x , z , θ ) = 0 g (d , x , z , θ ) ≤ 0

(8)

Re(eig(J )) < 0

A more convenient description of the feasible region is provided by Halemane and Grossmann who defined function ψ(d, θ).22 The set of θ that meets the criterion ψ(d, θ) ≤ 0 is the feasible region, with the definition of ψ(d,θ) as follows:

θ = θ N + δθ ̃ θ ̃ = [sin φ , cos φ]T φ ∈ [0, 2π ]

ψ (d , θ ) = min u u,z

s. t. f j (d , z , θ ) ≤ u , j ∈ J

Re(eig(J)) < 0 means that the real parts of the eigenvalues of the Jacobian matrix J of equality constraints h are all smaller than 0. Thus, the solution to eq 11 can be guaranteed to be stable. The feasible and flexible regions are shown in Figure 4b. As shown in Figure 4, the asterisk in the center represents the nominal operating point (CNmin,FN), and the rectangle represents the flexible region. The boxes and the solid dots on the boundary represent the unstable and stable points, respectively. 2.3. Comparison of Two Cases. To compare the difference between the above-mentioned two cases, the two images are superimposed, and the difference between the two flexible regions is shaded as shown in Figure 5. The larger outside region and larger rectangle correspond to the case when stability is not considered, and the smaller region and smaller rectangle correspond to the case when stability constraints are considered. The difference between the two cases is significant: when stability constraints are added, the boundary of feasible region becomes stable, and the flexible region shrinks accordingly. Unstable parts therefore exist in the traditional flexible region, represented by the shaded area in Figure 5, which demonstrates the necessity for analyzing process flexibility and stability simultaneously. Two points in Figure 5 are chosen to illustrate the difference between stable and unstable cases. Point 1 follows: Cmin = 6.42 kmol m−3 and F = 1 m3 h−1, which corresponds to the stable case. Point 2 follows: Cmin = 6.6 kmol m−3 and F = 1 m3 h−1, which corresponds to the unstable case. The bifurcation diagrams versus Fcw are shown in Figure 6a,b, respectively. The shaded area represents the feasible region, within which exists the stable part in Figure 6a. In Figure 6b, however, the line in shaded region is unstable, meaning that there does not exist an Fcw that satisfies performance specifications and stability conditions simultaneously.

(9)

The feasible region is bounded by ψ(d, θ) = 0, but this form is difficult to solve, so another formulation2 is used: the offset direction from the nominal point is fixed at first, then the maximum offset distance is found and marked with points. If the maximum offset distance is found for all directions, the feasible region is bounded by all the points, as shown in Figure 3.

Figure 3. Alternate formulation of feasible region.

Here, two cases are discussed: (1) not adding stability constraints but checking the stability of the boundary; (2) adding stability constraints to obtain a stable feasible region. When stability is not considered, for a fixed direction, the maximum offset distance is determined from δ *(θ )̃ = max δ s. t. h(d , x , z , θ ) = 0 g (d , x , z , θ ) ≤ 0 θ = θ N + δθ ̃

3. ALGORITHMS FOR ANALYZING FLEXIBILITY AND STABILITY SIMULTANEOUSLY 3.1. Vertex Algorithm. The basic idea of algorithms for the simultaneous analysis of flexibility and stability is to embed stability constraints into the optimization model for calculating the flexibility index. To simplify the problem, the assumption that all values of x should be stable instead of at least one x being stable was made. Flexibility index algorithms have been provided

θ N = [Cmin , F ]T θ ̃ = [sin φ , cos φ]T φ ∈ [0, 2π ]

(11)

(10)

Here, h represents the equality constraints in eq 7 with the left term being equal to 0, and g represents the inequality constraints 14767



Article

Figure 4. Feasible and flexible regions: (a) without considering stability and (b) considering stability.

There is no special treatment for the stability constraints; all real parts of the eigenvalues of the Jacobian matrix J of equality constraints h must be smaller than 0. Thus, the problem is defined as follows: k

δ *(θ ̃ ) = max δ s. t. h(d , x , z , θ ) = 0 g (d , x , z , θ ) ≤ 0 Re(eig(J )) < 0 k

θ = θ N + δθ ̃

A schematic diagram of the vertex algorithm is shown in Figure 7. As shown in Figure 4, for the MMA polymerization process, the critical point lies at the hyper-rectangle vertex, so the vertex algorithm is suitable. This process has two uncertain parameters, so four vertices exist. Cmin and F are denoted as θ1 and θ2, and the maximum positive and negative deviations of the two uncertain parameters are denoted as Δθ+1 , Δθ−1 , Δθ+2 , and Δθ−2 with the specific values given in section 2.1. The four vertices are denoted as θ̃1(Δθ+1 , Δθ+2 ), θ̃2(Δθ−1 , Δθ+2 ), θ̃3(Δθ+1 , Δθ−2 ), θ̃4(Δθ−1 , Δθ−2 ). For comparison, the maximum offset distances of the four vertices considering stability and without considering stability are all listed in Table 2. The flexibility index F corresponds to the minimum of δ*(θ̃k) and is equal to 0.274 when considering stability. However, the

Figure 5. Comparison of two cases.

in the Introduction. Here, the vertex algorithm and active constraint algorithm are studied. In the vertex algorithm, if the critical point corresponds to a hyper-rectangle vertex, then the θ̃ of eq 11 corresponds to the vertex direction. For the case of p uncertain parameters, there are 2p vertices in total, so θ̃ ∈ {θ̃k, k = 1,2, ..., 2p}, and the element in the kth direction is Δθ+j or −Δθ−j . As the number of offset directions is limited, so the maximum offset distance δ*(θ̃k) can be obtained by enumerating all offset directions, and the flexibility index corresponds to the minimum of δ*(θ̃k). k

F = min δ*(θ ̃ ), k = 1, ..., 2p k

(13)

(12)

Figure 6. Bifurcation diagram versus Fcw: (a) stable case and (b) unstable case. Note: LP = limit point, H = Hopf point. 14768



Article

obtained for the matrices larger than 4 × 4.23 The simultaneous calculation algorithm of flexibility and stability is described in Figure 8.

Figure 7. Vertex algorithm.

Table 2. Four Vertices of MMA Polymerization Reaction vertex θ̃k ̃θ1(Δθ+1 , Δθ+2 ) θ̃2(Δθ−1 , Δθ+2 ) θ̃3(Δθ+1 , Δθ−2 ) θ̃4(Δθ−1 , Δθ−2 )

δ*(θ̃k) with stability

δ*(θ̃k) without stability

0.274 1.729 0.325 0.846

0.495 1.729 0.495 0.858

Figure 8. Active constraint algorithm.

The traditional flexibility index is solved directly without consideration of stability. The stability of the obtained flexible region must be checked. If all real parts of the eigenvalues of the Jacobian matrix J of equality constraints h are smaller than 0, then the region is stable. If not, it is unstable, and stability constraints should be embedded in the problem. If the stability constraint is added, then there exists an active constraint where the real part of the eigenvalue of the Jacobian matrix J of equality constraints h is equal to 0. Here we discuss only the case where one of the real parts of the eigenvalues is equal to 0. Since one cannot specify an eigenvalue to be equal to 0, the active constraint is max Re(eig(J)) = ε. ε is a small negative value; therefore, the real parts of the other eigenvalues will be smaller than 0. nz ordinary active constraints need to be determined, and there is no uniform method to find these nz active constraints besides using the enumeration approach. δ should be calculated for each group of active constraints, and the flexibility index is equal to the minimum of δ. We first discuss the case without considering stability and define the partial derivatives matrix of f j with respect to z as follows:

flexibility index without considering stability is 0.495 and corresponds to Figure 5. This further illustrates that it is necessary to consider stability when calculating flexibility. The vertex algorithm is simple but is only suitable for cases where the critical point lies at the vertex. The scale of the problem increases exponentially along with increase in the number of uncertain parameters. Therefore, we need to develop an algorithm that does not rely on the assumption that the critical point lies at the vertex of the hyper-rectangle. 3.2. Active Constraint Algorithm. The active constraint algorithm was presented by Grossmann and Floudas.5 Using this algorithm, the simultaneous calculation of process flexibility and stability are studied. The general form of the simultaneous calculation of flexibility and stability is as follows: F = min δ s. t. ψ (d , θ ) = 0 ψ (d , θ ) = min u

⎛ ∂f1 ∂f2 ∂f ⎞ , ..., m ⎟ ⎜ , ∂z ⎠ ⎝ ∂z ∂z

z

s. t. hi(d , z , x , θ ) = 0 g j (d , z , x , θ ) ≤ u

(15)

If all the nz × nz submatrices are full rank, then the number of active constraints ( f j(d,z,θ) = u) is equal to nz + 1, where nz represents the number of control variables. For a fixed θ, the number of unknowns in the nz + 1 active constraints is equal to nz + 1; therefore, z and u can be solved. For the remaining inactive constraints, ( f j(d, z, θ) < u), so ψ(d, θ) = u. For the flexibility index problem, the optimal solution θc must lie at the boundary, so ψ(d, θc) = u and u is equal to 0 in the optimal solution.

Re(eig(J )) < 0 T (δ) = {θ|(θ N − δ Δθ −) ≤ θ ≤ (θ N + δ Δθ +)}

m ≥ nz + 1

(14)

This is a two-level optimization problem. Generally, the inner problem can be transformed into a single-level optimization problem using the Karush−Kuhn−Tucker(KKT) condition; however, explicit expressions for the eigenvalues cannot be 14769



Article

Equation 9 can be expressed in terms of the KKT condition as follows: (a)

F=

s. t. hi(d , x , z , θ) = 0

∑ γj = 1

g j (d , z , x , θ ) ≤ 0

j∈J

(b)

∑ γj j∈J

∂f j ∂z

∑ γj = 1

=0

j∈J

∑ γj

(c) sj = u − f j (d , z , θ )

j∈J

(d) γjsj = 0

∑ γj

(e) γj ≥ 0, sj ≥ 0

j∈J

(16)

δ≥0

T/K

Tj/K

δ

353.5

0.0029

60

333.8

0.495

γj ≥ 0

sj ≥ 0

(18)

Table 6. Active Constraints and Resultsa

Table 3. Results of MMA Polymerization Process 1

0.025

yj = {0, 1}

For the MMA polymerization process, when not considering stability, the SBB solver of GAMS is used to solve this MINLP problem. Here only partial variable values are given, as shown in Tables 3 and 4. As shown in Table 3, the flexibility index without considering stability is 0.495, which corresponds to the previous result. Since there is only one control variable, there should be two active constraints. As shown in Table 4, the active constraints g2 and g10 are Cm = 6.0171 kmol m−3 and D1 = 60 kg m−3 which corresponds with Table 3. To test the stability of the result, the eigenvalues of the Jacobian matrix J of equality constraints h should be calculated. The real parts of its eigenvalues are listed in Table 5. The real part of eigenvalue larger than 0 in Table 5 indicates that the results are unstable. Therefore, the stability constraint should be added and the problem resolved. For the MMA polymerization process, there is only one control variable, so two constraints are active. Let max Re(eig(J)) = −0.001. The other active constraints should be determined. Since the total number of inequalities is 14, there are 14 possible active constraints. All possible active constraints and results are listed in Table 6.

In eq 17, if yj = 1, then γj ≥ 0 and sj = 0, which indicates that f j is active. If yj = 0, then γj = 0 and sj ≥ 0, which indicates that f j is inactive. This also illustrates that eqs 17 and (16d), (16e) are equivalent. The original flexibility index problem can be transformed into a mixed integer nonlinear programming (MINLP) problem by eq 16 and eq 17. In general, the state variables x are difficult to express explicitly by equality constraints; therefore, the full form is as follows:

6.0171

i∈I

∂hi =0 ∂x

θ − δ Δθ − ≤ θ ≤ θ N + δ Δθ +

(17)

D1/ kg m−3

∑ μi

j

i∈J

D0/ kg m−3

+

∂hi =0 ∂z

N

∑ yj = nz + 1

CI / kmol m−3

∑ μi i∈I

∂gj ∂x

+

∑ yj ≤ nz + 1

sj − U (1 − yj ) ≤ 0

Cm/ kmol m−3

∂z

sj − U (1 − yj ) ≤ 0

γj − yj ≤ 0

γj , sj ≥ 0

∂gj

γj − yj ≤ 0

sj is the slack variables, and γj is the KKT multipliers for constraint j. Equation 16, parts a and b, show the stationary Lagrange conditions with respect to u and z, respectively. Equation 16, part c, defines the slack variables, and parts d and e represent the complementary conditions.5 When nz + 1 constraints are active, then eq 16 is the sufficient and necessary condition for a local minimum in eq 9, and this applies to both convex and nonconvex constraints. When the constraints are quasiconvex with respect to z, then eq 16 defines the global minimum. To illustrate the choice of active constraints explicitly, a set of 0−1 variables yj will be defined. yj = 0 indicates that f j is not an active constraint, and yj = 1 indicates that f j is active. Thus, eq 16, parts d and e, can be replaced by

yj = 0, 1

δ

min

θ , x , z , δ , sj , μj , γj , yj

active constraints flexibility index active constraints flexibility index a

Table 4. Results of MMA Polymerization Process 2

Cm LB 2.456 D0 UB NS

Cm UB 0.274 D1 LB NS

CI LB NS D1 UB NS

CI UB 0.838 Tj LB NS

T LB NS Tj UB NS

T UB NS Fcw LB 1.595

D0 LB NS Fcw UB NS

LB = lower bound, UB = upper bound, and NS = no solution.

y1

y2

y3

y4

y5

y6

y7

y8

y9

y10

y11

y12

y13

y14

0

1

0

0

0

0

0

0

0

1

0

0

0

0

According to Table 6, the flexibility index is 0.274 when adding stability constraints. This corresponds with the results from the vertex algorithm in section 3.1.

4. EIGENVALUE OPTIMIZATION ALGORITHM Eigenvalues play an important role in structural design, quantum mechanics, and other fields. When eigenvalues become constraints or objective functions, this kind of problem is termed eigenvalue optimization. In chemical engineering, eigenvalue

Table 5. Real Parts of Optimal Solution Eigenvalues real parts

1

2

3

4

5

6

−10.99

−10.99

−29.73

1.28

−11.11

−11.11 14770



Article

det(Pi) represents the successive principal minors of P. This is a twolevel optimization problem, similar to the active constraint algorithm. Equation 20 can be expressed in terms of the KKT condition.

optimization techniques are used in dynamic system design, analysis, and control. As already mentioned, the eigenvalue expression for a matrix larger than 4 × 4 cannot be obtained. Furthermore, even if an analytical expression could be obtained, they are usually complex and nonconvex.24 With reference to the eigenvalue optimization method,25−29 this paper proposes a new formulation for the simultaneous analysis of process flexibility and stability. Through Lyapunov’s matrix equality shown in eq 19, a matrix A (here, the Jacobian matrix J of equality constraints h) is associated with a real symmetric matrix P. For any positive definite matrix Q (usually Q = I, I is the identity matrix), if eq 19 has a unique solution, then the positive definiteness of P is equivalent to the real parts of the eigenvalues of A which are all smaller than 0. The eigenvalues of A are denoted as {λ1, λ2, ..., λn}. According to the theorem,28 eq 19 has a unique solution if and only if λi + λj ≠ 0, ∀i,j. Therefore, when the sum of any two eigenvalues of A is not equal to zero, the condition that the real part of the eigenvalues is smaller than zero can be transformed into the positive definiteness of a real symmetric matrix. T

A P + PA + Q = 0

F = min δ s. t. hi(d , z , x , θ ) = 0 g j (d , z , x , θ ) ≤ 0 J T P + PJ + I = 0 di < 0

∑ γj = 1 j

∑ γj j

∑ γj

(19)

j

The definiteness of P can be judged by the successive principal minors. If the Jacobian matrix J is termed A and Q = I, then a new form of the simultaneous analysis of flexibility and stability can be obtained.

∂g ′j ∂z ∂g ′j ∂x′

+

∑ μi i

+

∑ μi i

∂hi′ =0 ∂z ∂hi′∂h =0 ∂x′

γjg ′j = 0 θ N − δ Δθ − ≤ θ ≤ θ N + δ Δθ + δ≥0

F = min δ

γj ≥ 0

(21)

θ ∈ T (δ)

g j (d , z , x , θ ) ≤ u

di denotes −det(Pi), and gj′ denotes the inequalities, including the product quality requirements gj and successive principal minors of P larger than zero. hj′ denotes the equalities, including the original six equalities hi and the Lyapunov matrix equality. x′ denotes the six state variables and variables in the Lyapunov matrix equality. J denotes the Jacobian matrix of equalities.

J T P + PJ + I = 0

5. CASE STUDIES

s. t. ψ (d , θ ) = 0 ψ (d , θ ) = min u z

s. t. hi(d , z , x , θ ) = 0

det(Pi) > 0

(20)

5.1. MMA Polymerization. There are six state variables, Cm, CI, T, D0, D1, Tj, one control variable Fcw, and two uncertain parameters Cmin and F in the MMA polymerization process. The jacobian matrix is as follows:

Since P is a real symmetric matrix, for the n-dimensional Jacobian matrix, there are (n2 + n)/2 variables denoted as p1, p2,..., p21 in the Lyapunov matrix equality. Matrix P is defined in eq 22. A Lagrange multiplier exists for each equality and inequality constraint in the inner optimization. As there are 27 equality

constraints (including six original differential equations and 21 Lyapunov matrix equalities), the Lagrange multipliers for the equalities are denoted as μ1, μ2,..., μ27. There are 20 inequality constraints (including 14 inequalities for the range of process parameters in Table 1 and six for the successive principal minors

N

−

N

+

T (δ) = {θ|θ − δ Δθ ≤ θ ≤ θ + δ Δθ } δ≥0

14771



Article

The first three terms of the successive principal minors are

larger than zero). The Lagrange multipliers for the inequalities are denoted γ1, γ2,..., γ20. Besides the constraints mentioned, there are also 20 complementary conditions. ⎡ p1 p2 p3 p4 p5 p6 ⎤ ⎢ ⎥ ⎢ p2 p7 p8 p9 p10 p11 ⎥ ⎢p p p p p p ⎥ ⎢ 3 8 12 13 14 15 ⎥ P=⎢ p p p p p p ⎥ ⎢ 4 9 13 16 17 18 ⎥ ⎢ p5 p10 p14 p17 p19 p20 ⎥ ⎢ ⎥ ⎣ p6 p11 p15 p18 p20 p21 ⎦ (22)

shown in eq 23, and the first two terms of the Lyapunov matrix equalities are shown in eq 24. ⎧ d1 = −det(P1) = −p 1 ⎪ ⎪ d = −det(P ) = p 2 − p p 2 2 1 7 ⎪ 2 ⎨ ⎪ d3 = −det(P3) ⎪ = p 2 p − 2p p p + p2 p + p p2 − p p p 2 12 2 3 8 3 7 1 8 1 7 12 ⎪ ⎩ ...

⎧ ⎡F 2P0ΔHk p ⎤ ⎪−2⎢ + P0(k fm + k p)⎥p1 − p3 + 2P0k fmp4 + 2M mP0(k fm + k p)p5 + 1 = 0 ⎣ ⎦ Cpρ V ⎪ ⎪ Cm(k fm + k p)P0 P0ΔHk p ⎪ ⎡F ⎤ ⎪− p1 ⎨ ⎢⎣ V + P0(k fm + k p)⎥⎦p2 − C ρ p8 + P0k fmp9 + M mP0(k fm + k p)p10 − 2C I p ⎪ ⎪ CmM m(k fm + k p)P0 C ΔHk pP0 [C k + 1e − 6P0(k tc + 2k td)]P0 ⎪−⎜⎛k + F ⎟⎞p − m p5 = 0 p4 + p3 + m fm I 2 ⎪ ⎝ 2CpρC I 2C I 2C I V⎠ ⎪ ⎩ ...

Equation 21 is a complex nonlinear programming problem, and it is difficult to obtain a global optimal solution. The strategy is to use MultiStart algorithm in MATLAB combined with a SQP local solver. This is more likely to yield the global optimal solution. There are 50 random initial points given, and the smallest local optimum is chosen as the optimal solution. The final result is 0.274, which is consistent with the vertex algorithm and active constraint algorithm. It takes about 315 s to solve this problem on a 3.4 GHz quad-core CPU. 5.2. 1,3-Propanediol Fermentation. The mathematical model of 1,3-propandiol fermentation is as follows, and the values of the parameters are available in the Supporting Information.30

⎛ dy y ⎞ = d(y0 − y) − xφ1⎜β1 + γ1u + ⎟ dτ y + α1 ⎠ ⎝

6. DISCUSSION To compare the scope of application of the novel algorithm with the two other algorithms, the following discussion is given. The applicable condition for the vertex algorithm is that the critical point should correspond to the vertex of a hyper-rectangle.

⎛ y ⎞ dz = xφ2⎜β2 + γ2u + ⎟ − zd dτ y + α2 ⎠ ⎝ y (1 − y)(1 − z) y + α3

(24)

is 0.345, taking about 24 s on a 3.4 GHz quad-core CPU. For comparison, the maximum offset distances of the four vertices considering stability and without considering stability are all listed in Table 8. As shown in Table 8, the flexibility index considering stability is 0.345, which is consistent with the eigenvalue optimization algorithm. However, the flexibility index without considering stability is 0.453. The boundaries of feasible region and flexible region are depicted in Figure 9. The left boundary of feasible region without considering stability is the unstable boundary in Figure 9, and the right boundary is outside the figure. Both of the boundaries are unstable. When considering stability, the right boundary moves to the stable boundary in Figure 9, but the left boundary is almost the same. The larger rectangle corresponds to the case when stability is not considered, and the smaller rectangle corresponds to the case when stability is considered.

dx = x(u − d) dτ

u=

(23)

(25)

Table 7. Ranges of Process Parameters

The state variables are x, y, and z, which mean dimensionless biomass concentration, substrate concentration, and product concentration, respectively. The operating variable is dimensionless dilution rate d, and the two uncertain parameters are dimensionless feed substrate concentration y0 and dimensionless parameter α3. Given the nominal operating point, yN0 = 0.228, αN3 = 1.4 × 10−4, the maximum positive and negative deviations are Δθ+1 = 0.02, Δθ−1 = 0.02, Δθ+2 = 2 × 10−5, and Δθ−2 = 2 × 10−5. The range of process parameters is given in Table 7. Using the eigenvalue optimization algorithm proposed in section 4, the flexibility index considering stability is obtained. Similarly, 50 random initial points are given, and the smallest local optimum is chosen as the optimal solution. The final result

parameters

ranges

x z d

0.29−0.6 ≥0.58 0.3−0.4

Table 8. Four Vertices of 1,3-Propanediol Fermentation vertex θ̃k ̃θ1(Δθ+1 , Δθ+2 ) θ̃2(Δθ−1 , Δθ+2 ) θ̃3(Δθ+1 , Δθ−2 ) θ̃4(Δθ−1 , Δθ−2 ) 14772

δ*(θ̃k) with stability

δ*(θ̃k) without stability

0.345 0.346 0.456 0.453

3.861 3.945 0.456 0.453



Article

1 048 576.4 For large-scale problems, the vertex algorithm is unrealistic. The key of the active constraint algorithm is to find all possible active constraints, and the requirement that the critical point lies at the vertex no longer exists. The active constraint algorithm is accomplished in two steps, the first of which is solving the MINLP problem defined in eq 18. Conditions for the global optimal solution are given in the literature:5 f j(d, z, θ) are jointly quasiconcave in z and θ and strictly quasiconvex in z for fixed θ. For the second step of the active constraint algorithm, one of the active constraints will be that the eigenvalue is equal to zero. The situation where only one eigenvalue is equal to zero is considered, and the active constraint is max Re(eig(J)) = ε. There is no general method to find the remaining nz active constraints besides trying all possible combinations. For cases where the number of inequality constraints and control variables is small, this algorithm is suitable, but when many inequality constraints and control variables exist, this algorithm can lead to significant computations. The eigenvalue optimization algorithm is a general method, in which the necessary condition of stability where the real part of the eigenvalues is smaller than zero is transformed into the positive definiteness of a real symmetric matrix so that the stability condition is expressed explicitly. When the inner optimization is convex, the KKT condition is necessary and sufficient to obtain the extremum. However, if the inner optimization is nonconvex, the KKT condition is merely a necessary condition, and the result from this algorithm is the upper bound of the flexibility index. Referring to the global optimization algorithm α-BB,10,31,32 the true flexibility index can be obtained. As shown in Figure 10, first, the original nonlinear programming problem is solved, and the upper bound of the flexibility index is obtained. After convexification, the lower bound can be obtained by the same method, then branch on θ and z until the gap between the upper bound and lower bound is smaller than a given value. The convexification of Lyapunov equation and the positive definiteness of the successive principal minors of matrix P is our future work. Finally, the scopes of application and limitations of the three algorithms are summarized in Table 9.

Figure 9. Flexible regions with and without stability.

Figure 10. Flowchart of eigenvalue optimization method to obtain global optimal solution.

Convexity of the feasible region is not required, but the feasible region must be one-dimensional quasiconvex (1-DQC),2 meaning that the boundary function ψ(d, θ) is convex along the coordinate axis direction. Generally, an explicit expression of ψ(d, θ) is difficult to obtain, so we must judge according to the constraint function f j(d, z, θ). If f j(d, z, θ) are jointly quasiconvex in z and one-dimensional quasiconvex in θ, then the function ψ(d, θ) is one-dimensional quasiconvex in θ. The critical point lies at the vertex, and the vertex algorithm is suitable. Besides the limitations of constraints, since the vertex algorithm is an enumeration method, the scale of the problem should not be too large. The number of vertices which must be evaluated is given by 2p. When p = 10, 1024 subproblems defined by eq 13 are required. When p = 20, the number of subproblems grows to

7. CONCLUSIONS Because previous studies have not considered stability when calculating flexibility, this paper first illustrates the necessity for the simultaneous analysis of flexibility and stability. Then, two algorithms developed by embedding stability constraints in the optimization model for calculating the flexibility index are studied and a new algorithm is proposed. Finally, the applicable scope of the three algorithms is compared. In summary, the main results are as follows: (1) We have illustrated the necessity for the simultaneous analysis of process flexibility and stability. For the MMA polymerization process, the cases with and without consideration of stability were analyzed from which we verified the necessity for this research work. (2)

Table 9. Scopes of Application and Limitations of the Three Algorithms algorithms vertex algorithm active constraint algorithm eigenvalue optimization algorithm

scopes of application

limitations

The critical point should correspond to the vertex of the hyper-rectangle. The number of vertices should not be too large. The condition that the critical point corresponds to the vertex is not required, but There is no unified method to find the active constraints that there should be only one active constraint that contains eigenvalue. contain eigenvalues. It is a programming method without the assumption that the critical point lies at The inner optimization should be convex; otherwise, the upper the vertex and can be used to solve the problem effectively. bound is obtained. But this can be solved by the algorithm in Figure 10. 14773



Article

We studied the vertex algorithm and active constraint algorithm for considering flexibility and stability simultaneously, and we explored the inadequacies of the two algorithms. (3) For the eigenvalue optimization method, the necessary condition of stability that the real part of the eigenvalues is smaller than zero is transformed into the positive definiteness of a real symmetric matrix so that the stability condition is expressed explicitly. (4) The scope of application is discussed for the three algorithms mentioned above. The vertex algorithm is an enumeration method for the small-scale problem under the condition that the critical point lies at the vertex of the hyper-rectangle. The active constraint algorithm does not require the condition that the critical point corresponds to the vertex, but there is no unified method to find the active constraints that contain eigenvalues. The proposed eigenvalue optimization algorithm is a programming method without the assumption that the critical point lies at the vertex, and the problem of the simultaneous analysis of process flexibility and stability can be solved effectively.

■

■

REFERENCES

(1) Grossmann, I. E.; Morari, M. Operability, resiliency, and flexibility: Process design objectives for a changing world.Proceedings of the 2nd International Conference on Foundations of Computer Aided Process Design; Westerberg, A., Chien, H. H., Eds.; American Institute of Chemical Engineers: Ann Arbor, MI, 1983. (2) Swaney, R. E.; Grossmann, I. E. An index for operational flexibility in chemical process design. Part I. AIChE J. 1985, 31 (4), 621−630. (3) Grossmann, I. E.; Straub, D. A. Recent developments in the evaluation and optimization of flexible chemical processes. Batch Processing Syst. Eng. 1996, 143, 495−516. (4) Swaney, R. E.; Grossmann, I. E. An index for operational flexibility in chemical process design. Part II. AIChE J. 1985, 31 (4), 631−641. (5) Grossmann, I. E.; Floudas, C. A. Active constraint strategy for flexibility analysis in chemical process. Comput. Chem. Eng. 1987, 11 (6), 675−693. (6) Ostrovsky, G.; Volin, Y. M.; Barit, E.; Senyavin, M. Flexibility analysis and optimization of chemical plants with uncertain parameters. Comput. Chem. Eng. 1994, 18 (8), 755−767. (7) Pistikopoulos, E.; Ierapetritou, M. Novel approach for optimal process design under uncertainty. Comput. Chem. Eng. 1995, 19 (10), 1089−1110. (8) Ostrovsky, G.; Achenie, L.; Wang, Y.; Volin, Y. A new algorithm for computing process flexibility. Ind. Eng. Chem. Res. 2000, 39 (7), 2368− 2377. (9) Raspanti, C.; Bandoni, J.; Biegler, L. New strategies for flexibility analysis and design under uncertainty. Comput. Chem. Eng. 2000, 24 (9), 2193−2209. (10) Floudas, C. A.; Gumus, Z. H.; Ierapetritou, M. G. Global optimization in design under uncertainty: Feasibility test and flexibility index problems. Ind. Eng. Chem. Res. 2001, 40 (20), 4267−4282. (11) Ostrovsky, G.; Achenie, L.; Wang, Y.; Volin, Y. A unique approach for solving sub-problems in flexibility analysis. Chem. Eng. Commun. 2002, 189 (1), 125−149. (12) Ostrovsky, G.; Datskov, I.; Achenie, L.; Volin, Y. M. Process uncertainty: case of insufficient process data at the operation stage. AIChE J. 2003, 49 (5), 1216−1232. (13) Wang, H. Z.; Chen, B. Z.; He, X. R.; Zhao, J. S.; Qiu, T. Numerical analysis tool for obtaining steady-state solutions and analyzing their stability characteristics for nonlinear dynamic systems. J. Chem. Eng. Jpn. 2010, 43 (4), 394−400. (14) Wang, H.; Yuan, Z.; Chen, B.; He, X.; Zhao, J.; Qiu, T. Analysis of the stability and controllability of chemical processes. Comput. Chem. Eng. 2011, 35 (6), 1101−1109. (15) Gritton, K. S.; Seader, J.; Lin, W. J. Global homotopy continuation procedures for seeking all roots of a nonlinear equation. Comput. Chem. Eng. 2001, 25 (7), 1003−1019. (16) Elnashaie, S. S.; Grace, J. R. Complexity, bifurcation and chaos in natural and man-made lumped and distributed systems. Comput. Chem. Eng. 2007, 62 (13), 3295−3325.

ASSOCIATED CONTENT

S Supporting Information *

Tables containing additional parameter values. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

f * = initiator efficiency F = monomer feed stream volumetric flow rate (m3/h) FI = initiator feed stream volumetric flow rate (m3/h) Fcw = cooling fluid feed stream volumetric flow rate (m3/h) Mm = monomer molecular weight (kg/kmol) R = universal ideal gas law constant (kJ/kmol K) T = reactor temperature (K) Tin = monomer and initiator feed stream temperature (K) Tj = cooling fluid temperature (K) Tw0 = cooling fluid feed stream temperature (K) t = time (h) U = global heat-transfer coefficient (kJ/h K m2) V = reactor volume (m3) V0 = cooling fluid volume (m3) ΔH = propagation reaction heat of reaction (kJ/kmol)

AUTHOR INFORMATION

Corresponding Author

*Tel.: +86 10 62781499. Fax: +86 10 62770304. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge support from the National Basic Research Program (No. 2012CB720500) and the National Natural Science Foundation of China (No. 21306100).

■

NOMENCLATURE A = heat transfer area (m2) Afm = monomer transfer reaction Arrhenius pre-exponential factor (m3/kmol h) AI = initiation reaction Arrhenius pre-exponential factor (h−1) Ap = propagation reaction Arrhenius pre-exponential factor (m3/kmol h) Atc = coupling termination reaction Arrhenius pre-exponential factor (m3/kmol h) Atd = disproportionation termination reaction Arrhenius preexponential factor (m3/kmol h) CI = initiator concentration (kmol/m3) CIin = initiator feed stream concentration (kmol/m3) Cm = monomer concentration (kmol/m3) Cmin = monomer feed stream concentration (kmol/m3) Cp = reaction mixture heat capacity(kJ/kg K) Cpw = cooling fluid heat capacity (kJ/kg K) D0 = molar concentration of polymer death chains (kmol/m3) D1 = mass concentration of polymer death chains (kg/m3) Efm = monomer transfer reaction activation energy (kJ/kmol) Ep = propagation reaction activation energy (kJ/kmol) EI = initiation reaction activation energy (kJ/kmol) Etc = coupling termination reaction activation energy (kJ/kmol) Etd = disproportionation termination reaction activation energy (kJ/kmol) 14774



Article

(17) Jalali, F.; Seader, J.; Khaleghi, S. Global solution approaches in equilibrium and stability analysis using homotopy continuation in the complex domain. Comput. Chem. Eng. 2008, 32 (10), 2333−2345. (18) Wang, H. Z.; Zhang, N.; Qiu, T.; Zhao, J. S.; He, X. R.; Chen, B. Z. A process design framework for considering the stability of steady state operating points and Hopf singularity points in chemical processes. Chem. Eng. Sci. 2013, 99, 252−264. (19) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Integrated flexibility and controllability analysis in design of chemical processes. AIChE J. 1997, 43 (4), 997−1015. (20) Wang, H. Z. Inherently Safer Design Oriented Analysis of SteadyState Multiplicity and Stability of Chemical Processes. Ph.D. Thesis, Tsinghua University, Beijing, 2010. (21) Silva-Beard, A.; Flores-Tlacuahuac, A. Effect of process design/ operation on the steady-state operability of a methyl methacrylate polymerization reactor. Ind. Eng. Chem. Res. 1999, 38 (12), 4790−4804. (22) Halemane, K. P.; Grossmann, I. E. Optimal process design under uncertainty. AIChE J. 1983, 29 (3), 425−433. (23) Blanco, A. M.; Bandoni, J. A. Design for operability: A singularvalue optimization approach within a multiple-objective framework. Ind. Eng. Chem. Res. 2003, 42 (19), 4340−4347. (24) Blanco, A. M.; Bandoni, J. A. Eigenvalue optimization-based formulations for nonlinear dynamics and control problems. Chem. Eng. Process: Process Intensification 2007, 46 (11), 1192−1199. (25) Lewis, A. S.; Overton, M. L. Eigenvalue optimization. Acta Numer. 1996, 5 (1), 149−190. (26) Ringertz, U. T. Eigenvalues in optimum structural design. Largescale Optimization with Applications; Springer: New York, 1997; pp 135− 149. (27) Kanno, Y.; Ohsaki, M. Necessary and sufficient conditions for global optimality of eigenvalue optimization problems. Struct. Multidiscip. Optim. 2001, 22 (3), 248−252. (28) Vidyasagar, M. Nonlinear Systems Analysis; Society for Industrial and Applied Mathematics: Philadelphia, 2002. (29) Blanco, A. b. M.; Bandoni, J. A. Interaction between process design and process operability of chemical processes: An eigenvalue optimization approach. Comput. Chem. Eng. 2003, 27 (8), 1291−1301. (30) Xiu, Z. L.; Zeng, A. P.; Deckwer, W. D. Multiplicity and stability analysis of microorganisms in continuous culture: Effects of metabolic overflow and growth inhibition. Biotechnol. Bioeng. 1998, 57 (3), 251− 261. (31) Adjiman, C. S.; Dallwig, S.; Floudas, C. A.; Neumaier, A. A global optimization method αBB for general twice-differentiable constrained NLPs-I. Theoretical advances. Comput. Chem. Eng. 1998, 22 (9), 1137− 1158. (32) Adjiman, C.; Androulakis, I.; Floudas, C. A global optimization method αBB for general twice-differentiable constrained NLPs-II. Implementation and computational results. Comput. Chem. Eng. 1998, 22 (9), 1159−1179.

14775


Novel Method for Considering Process Flexibility and Stability

Recommend Documents