Steady-State Process Optimization with Guaranteed Robust Stability

Ind. Eng. Chem. Res. 2005, 44, 2737-2753

2737

Steady-State Process Optimization with Guaranteed Robust Stability and Flexibility: Application to HDA Reaction Section M. Mo1 nnigmann and W. Marquardt* Process Systems Engineering, RWTH Aachen University, Templergraben 55, D-52056 Aachen, Germany

This paper extends a recently suggested approach for the steady-state process optimization with guaranteed robust stability and flexibility under parametric uncertainty. The approach is based on measuring the distance of a candidate point of operation to so-called critical manifolds. Critical manifolds locally separate regions of the space of the process and controller design parameters with desired process traits from those regions with undesired traits. Pairs of desired and undesired traits may, for example, be feasibility and infeasibility of operation, stability and instability, or more generally, desired dynamical behavior and undesired dynamical behavior. While in previous applications knowledge on the existence and location of critical manifolds were assumed to be available before the optimization was attempted, the present paper presents an algorithm in which critical manifolds are automatically detected as the optimization proceeds. This algorithm, which is the conceptual contribution of the paper, allows the application of the critical manifold-based approach to processes for which no a priori information on the existence and location of critical manifolds exists. As a proof of concept the algorithm is applied to the reaction section of the HDA process. An analysis of the critical manifolds of this process model is not available. Since 12 uncertain parameters exists, analyzing the critical manifolds would be tedious. While an analysis of the 12 uncertain parameters is not practical, the critical manifoldbased optimization approach can be applied to models with this number of parameters and beyond. Introduction A new approach to taking constraints on the process dynamics into account in steady-state process optimization has recently been suggested.1-3 The key features of this new approach are that (i) it can be applied to a broad class of dynamical features, among them stability and disturbance rejection; (ii) it guarantees the desired dynamical behavior despite parametric uncertainty in the process model; and (iii) it can be applied to ensuring the desired dynamics and process feasibility simultaneously. While the new approach is more general, a typical application is to guarantee stability and feasibility of the optimal point of operation despite parametric uncertainty of the process model when optimizing the process with respect to an economic objective function.1 The two central concepts behind the new approach are the notion of a critical manifold and the normal space to that manifold. Readers not familiar with the notion of a manifold may think of the nonlinear generalization of a linear subspace of a vector space. Just as a plane defined by a single linear equation is a linear subspace of the linear vector space R3, nonlinear equations or sets of equations can define sets of points in Rn. These so-called manifolds can locally be approximated by linear subspaces. Globally, however, nontrivial nonlinear manifolds may bend and fold to form nonlinear geometrical objects. The stability boundary is a typical example of a critical manifold. In a simple case, the stability boundary separates the part of the process design parameter space in which stable steady states exist from parts with unstable steady states. More * To whom correspondence should be addressed. Tel.: +49 241 804668. Fax: +49 241 802326. E-mail: marquardt@ lpt.rwth-aachen.de.

generally, critical manifolds locally separate parts of the parameter space for which the system has desired traits from parts with undesired traits. The normal space of the critical manifold is of great importance, because the closest connection between a candidate point of operation and a critical manifold occurs along a direction that is normal to the critical manifold. The normal direction can therefore be used to quantify the closest distance between a candidate point and the critical manifold. The concepts of the critical manifold and the normal vector are detailed in the next section. In applications of the critical manifold-based approach reported so far, a priori knowledge on the existence of critical manifolds had to be available for the given process model. Moreover, the location of the critical manifolds in the parameter space had to be disclosed before the approach could be applied. Methods for the analysis of critical manifolds by numerical continuation are established,4,5 and mature implementations are available.6-8 In fact these tools are well-accepted in the community and continue to be used in different contexts.9-14 Such an analysis is, however, often tedious if not impossible for processes with more than two or three uncertain parameters.15 The present paper suggests an algorithm for the steady-state process optimization that can be applied to process models with uncertain parameters for which no a priori information on the existence and location of the critical manifolds is available. To our knowledge such an algorithm has not been presented or even suggested before. The method avoids an analysis by continuation and therefore is not restricted to process models with a few uncertain parameters. As a proof of concept, the algorithm is applied to the reaction section of the HDA process.16 The process model consists of

10.1021/ie0495776 CCC: $30.25 © 2005 American Chemical Society Published on Web 02/04/2005

2738

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

about 100 ordinary differential equations and 370 algebraic equations. Twelve uncertain parameters exist. The paper is organized as follows. The next section states the system and problem class which the algorithm can be applied to. It is stressed that these subsections are a summary of more detailed discussions of the system and problem class given elsewhere.1,2 The material is repeated here to the extent necessary to understand the algorithm. The third and fourth section explain the building blocks of the algorithm and the algorithm itself, respectively. The results obtained by an application of the suggested algorithm to the reaction section of the HDA process are then reported in a separate section. Finally, a summary and future perspectives are given in the last section. Background System Class. This section summarizes the requirements that must hold for an application of the approach presented here. Assume a model of the form

x˘ d(t) ) f d(xd(t), xa(t), r(t), u(t), d (t), ϑ), xd(0) ) xd0 0 ) f a(xd(t), xa(t), r(t), u(t), d (t), ϑ) d

(1)

a

y(t) ) h(x (t), x (t), r(t), u(t), d (t), ϑ) is given. In eq 1, x ) (xdT, xaT)T ∈ Rnx, r ∈Rnr, u ∈ Rnu, d ∈Rnd, ϑ ∈ Rnϑ, and y ∈Rny are state variables, reference signals, inputs, disturbances, parameters, and outputs of the system. Furthermore, t denotes time, f ) (f dT, f aT)T, and h are smooth functions which map from some subset U ⊂ Rnx × Rnr × Rnu × Rnd × Rnϑ × R into Rnx and Rny, respectively. The index of this differentialalgebraic equation system is assumed to be one, and the initial conditions are understood to be consistent. The upper labels d and a distinguish dynamic state variables xd and differential equations f d, on one hand, from algebraic variables xa and algebraic equations f a on the other hand. According to Mo¨nnigmann and Marquardt,1 both open-loop and closed-loop systems can be modeled in the form of eq 1. In the closed-loop case, the ODEs in eq 1 comprise both ODEs of the controller such as an integrator of a PI-controller and the ODEs of the process model. Likewise, the algebraic equations in eq 1 comprise algebraic equations for the controller and the process model. Note that some inputs u may exist in the closed-loop system that are not occupied by any controller. The requirements for an application of the critical manifold-based approach can now be introduced: 1. The linearization of the algebraic equations f a in eq 1 is assumed to have full rank with respect to xa at all steady states of interest. 2. Inputs u(t) can be partitioned into a constant mean value and bounded time-dependent variations:

j i + ηui(t), ηulower e ηui(t) e ηuupper ∀t ui(t) ) u i i

(2)

where Re(λj) denotes the real parts of the eigenvalues λj of the Jacobian of f with respect to x. Condition 3 must hold at all steady states of interest. 3. The assumptions on u(t) must also hold for d(t) and r(t). Note that the these assumptions involve linearizations of the model (eq 1). It is stressed, however, that the method presented in this work captures the dynamics of the nonlinear system. The first assumption implies that the algebraic variables can locally be expressed as functions of the remaining variables according to the implicit function theorem (IFT): IFT

0 ) f a(xd, xa, r, u, d, ϑ, t) 98 xa ) xa(xd, r, u, d, ϑ, t) (4) More precisely, the implicit function theorem states that, for each steady state of the system in eq 1, an open neighborhood exists on which the algebraic variables can be expressed as a function of the remaining variables. On any such neighborhood the DAE system can therefore be reduced to an ODE system by substituting eq 4 into the nonlinear system (eq 1): In eq 5, the label

x˘ ) f(x, r, u, d, ϑ, t), x(0) ) x0 y ) h(x, r, u, d, ϑ, t)

d for the dynamic variables is omitted for simplicity. Furthermore, the time dependencies are omitted in x(t) to simplify the notation, and by a slight abuse of notation, we use the same symbols f and h in eqs 1 and 5. We stress that we are interested in reducing the DAE system (eq 1) to an ODE system from a conceptual point of view only, because it allows us to make use of results from the ODE bifurcation theory in the next section. From an implementational point of view on the other hand, it is not advisable to substitute explicit equations for algebraic variables in to the differential equations. This is briefly discussed below when the algorithm is introduced. The second and third assumption allow for the further simplification of the notation. According to eqs 2 and 3, inputs u, reference signals r, and disturbances d vary only quasistatically compared to the system dynamics. They can therefore be described by a mean value and a parametric uncertainty. Parameters ϑ are not timedependent and can be modeled in the same manner. Some of the parameters are, however, typically not uncertain, or the uncertainty can be neglected in the given model. Inputs, references, disturbances, and parameters are therefore partitioned into those that are uncertain and those for which uncertainty does not exist or can be neglected. Denoting uncertain inputs u, reference signals r, and disturbances d by R and denoting all remaining parameters by p, the differential equations of the system in eq 5 can be rewritten as

x˘ ) f (x, R, p), x(0) ) x0 y ) h(x, R, p)

Furthermore, the variations are slow compared to the time scales of the system (1), i.e.

1 dηui | | , |Re(λj)|, ∀i, j ηui dt

(3)

(5)

(6)

Since the inputs, references, disturbances, and uncertain parameters denoted by R are all assumed to be bounded, they can be described by intervals (0) Ri ∈ [R(0) i - ∆Ri, Ri + ∆Ri], i ) 1, ..., nR

(7)

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2739

manifolds due to dynamical traits. Boundaries in the space of the uncertain process parameters R that are critical for dynamical traits of a system can in general not be described by inequalities of the form given in eq 8. In contrast, defining equations for critical manifolds for dynamical traits can be stated. Bifurcation theory provides so-called augmented systems that can be used to define critical manifolds of the form given in eq 9 for stability boundaries and other types of dynamical boundaries.3,17 Fortunately, the concept of a critical manifold therefore allows us to treat both feasibility boundaries and dynamical constraints in a unified manner. The symbol Rc defined by

Rc ) ∪iRc,i

Figure 1. Abstract parametric flexibility problem.

It is stressed that the critical manifold-based approach can deal with more general descriptions of the uncertain parameters than the hypercube (eq 7) with a centered nominal value R(0)2,1. Since the application to the HDA reaction section uses the simple type of uncertainty (eq 7), the discussion is restricted to this case for simplicity. The remainder of the paper deals with systems of the form in eq 6 and uncertain parameters of the form in eq 7. Problem Class. The approach presented here is based on the idea of representing both feasibility constraints and constraints on the dynamical behavior of the system (eq 6) by so-called critical manifolds. The idea of a critical manifold is first illustrated with a feasibility constraint. A large class of feasibility constraints can be stated in the form

0 e g(x, R, p)

(8)

where g is defined on an appropriate subset of Rnx × RnR × Rnp, it maps into R, and it is smooth with respect to x, R, and p. The points at which the constraint (eq 8) is active define the critical manifold:

Mc(p(0)) ) {(xT, RT)T ∈ Rnx × RnR|0 ) f (x, R, p(0)), 0 ) g(x, R, p(0))} (9) where p(0) denotes the nominal values of the parameters that are not subject to uncertainty. The linearization of the defining equations with respect to x and R is assumed to have full rank.3 The critical manifold (eq 9) separates the part of the process parameters R in which the constraint (eq 8) holds from the part in which the constraint (eq 8) is violated. Usually more than one critical manifold exist. Figure 1 sketches a case with two such manifolds: Mc,1 and Mc,2. The regions in which the respective constraint is violated are denoted by Rc,1 and Rc,2. In Figure 1 and throughout the paper, the symbol R(0) denotes the nominal values of the uncertain parameters, and Rr refers to the values the uncertain parameters can attain. The hypercube (eq 7) of the uncertain parameters is overestimated by the ellipsoid

{

|∑( nR

Ri

i)1

∆Ri

M r ) R ∈ R nR

) }

(11)

is introduced to denote the union of all regions of critical parameter values bounded by feasibility and dynamical critical manifolds. Based on Rr and Rc, the process optimization problem addressed in the present paper can be stated:

φ(x(0), R(0), p(0))

(12a)

0 ) f (x(0), R(0), p(0))

(12b)

min

x(0), R(0), p(0)

s.t.

L ) Rr(p(0)) ∩ Rc(p(0))

(12c)

x(0) ∈ X, R(0) ∈ A, p(0) ∈ P

(12d)

The symbols X ⊂ Rnx, A ⊂ RnR, P ⊂ Rnp denote box constraints on the respective variables and parameters. The upper index zero denotes the nominal point of operation. In particular R(0) are the nominal values of the uncertain parameters. Equation 12b guarantees that the optimal point of operation is a steady state of eq 6. Equation 12c ensures that the regions bounded by Mc and Mr have an empty intersection; therefore, Mc is not crossed for any value R ∈ Rr. In this sense, the optimal point of operation resulting from eq 12 is robust despite the uncertainties (eq 7). It is noted that more general robustness regions than the hyperellipsoids (eq 10) can be treated. In the most general case, Rr is allowed to vary in shape and size when the nominal point R(0) is varied and both Rr and Rc may depend on p.1,2 Here we restrict ourselves to the case of a robustness ellipsoid, since this case is used in the application to the HDA reaction section. Building Blocks of the Algorithm Measuring Distance in the Space of the Uncertain Parameters. To state the constraint (eq 12c), repeated here for convenience,

L ) Rr(p(0)) ∩ Rc(p(0))

2

) nR

(10)

in order to obtain a smoothly bounded region Rr. Having illustrated the concept of a critical manifold with a feasibility constraint, we now turn to critical

in a form that can be implemented, the distance of Rr and Rc in the space of the uncertain parameters R must be quantified. The parameter space distance can be measured with the aid of normal vectors to the critical manifolds. This idea is discussed in the present subsection.

2740


Figure 2. Normal vector to critical manifolds.

First the uncertain parameters R are scaled according to

Ri Ri f , i ) 1,..., nR ∆Ri

R(0) ) R(1) + l (13)

It is stressed that this scaling is not introduced for numerical but for physical reasons. In general, the uncertain process model parameters will not all have the same physical units. In the HDA reaction section model, for example, both a feed rate (measured in kmol min-1) and a feed temperature (measured in K) are uncertain. From a physical point of view it is not meaningful to measure distance in such a space of parameters of heterogeneous units. Scaling the parameters according to eq 13 renders them dimensionless, however, and distance can be measured with the 2-norm in the dimensionless parameter space. After scaling, a unit step along a coordinate axis of the rescaled parameter number i corresponds to step of ∆Ri in the unscaled parameter space. It is stressed that the choice of the uncertainties ∆Ri has a strong impact on the result of the optimization; therefore, uncertainties should be assigned with care. From a physical point of view, changing the uncertainties ∆Ri corresponds to requesting a different robustness and therefore will generally result in a different point of robust process operation. From a mathematical point of view, changing ∆Ri will result in a different scaling of the parameters (eq 13), which in turn gives rise to different critical manifolds and normal vectors. The observation that normal vectors and critical manifolds are not invariant under scaling has been discussed in detail in the context of hydraulic servo-systems.18 For simplicity, it is assumed in the sequel that the parameters have been scaled according to eq 13. Note that in the scaled parameters, the hyperellipsoid (eq 10) that bounds the robustness region becomes the hyperball nR

Mr ) {R ∈ RnR|

∑i Ri2 ) nR}

the nR-dimensional ball in Figure 2a. The sketch in Figure 2b illustrates that the critical point and the candidate point of operation have the same values of the parameters not affected by uncertainty p. While obvious intuitively, the idea of measuring the closest distance to a critical manifold along the normal vector can be substantiated theoretically by comparing the normal vector to the Lagrange multipliers of a optimization problem for the closest distance between Mr and Mc.1,2 Assume a normal vector n is known that spans the direction along which the closest connection exists between Mc and the center R(0) of the hyperball (eq 14). Denoting by R(1) the point at which this normal direction intersects Mc, the constraint

(14)

The locally closest connections between the hyperball (eq 14) and the critical manifolds Mc,i occur along directions that are normal to the critical manifolds. See Figure 2a for a sketch. The parameter values of the uncertain parameters R at the candidate point of operation and the closest critical point are denoted with an upper index of zero and one, respectively. The normal vector is given by the short bold line in Figure 2 that is normal to the critical manifold Mc,1. The distance between the two points is equal to the radius xnR of

n ||n||2

l g xnR, l ∈ R

(15)

is a necessary condition for L ) Rr(p(0)) ∩ Rc(p(0)) to hold. The inequality constraint in eq 15 ensures that the distance between R(0) and R(1) is larger than or equal to the radius xnR of the hyperball (eq 14). Parametric robustness is ensured by eq 15 in the sense that the critical manifold will not be crossed regardless of the actual value of R ∈ Mr where Mr is defined by eq 14. This is illustrated in Figure 2a. While the sketch in this figure is two-dimensional, the concept of imposing a lower bound on the distance to critical manifolds carries over to spaces of the uncertain parameters of arbitrary dimension nR. It is stressed that the number of equations in the constraint (eq 15) grows only linearly in the number of uncertain parameters nR. Critical Manifolds for Generalized Constraints on Eigenvalues. According to the previous section feasibility constraints of the form in eq 8 give rise to critical manifolds of the type in eq 9, repeated here for convenience

Mc(p(0)) ) {(xT, RT)T ∈ Rnx × RnR| 0 ) f (x, R, p(0)), 0 ) g(x, R, p(0))} While a simple constraint of the form in eq 8

0 e g(x, R, p)

(16)

does not exist for the critical dynamical points such as stability boundaries, a meaningful definition of a critical manifold does exist. The present section gives details on a particular critical manifold that allows us to impose a bound on the rate of decay to steady state after a disturbance. As detailed below, this type of manifold is related to the location of the eigenvalues of the linearized system. Since this type of manifold furthermore can be looked upon as a generalization of the stability boundary,2,19 it is referred to as the critical manifolds for generalized constraints on eigenvalues in the sequel. The focus is on this kind of critical manifold, since it is this kind of critical manifold that is used in the application to the HDA reaction section. It is noted, however, that the critical manifold-based approach can be applied to a more general class of critical manifolds than the particular type discussed here. References to more general discussions are given at the end of this section.


symbol ˜f is introduced in eq 18 for ease of reference, and the arguments of f and ˜f have been omitted for simplicity. The first two parts of ˜f in eq 18 result from splitting the eigensystem

fxv ) λv

Figure 3. According to proposition 1, the decay rate can be bounded above to linear order. This amounts to shifting the upper bound on the real parts to the left (a). As a simple extension both real and imaginary part can be bounded (b).

The critical manifolds for generalized constraints on eigenvalues are based on the following result from nonlinear systems theory, which is cited without proof: Proposition 1. Assume x(t) ) x(ss)(R(ss), p(ss)) for all t > t0 is a steady state of eq 6 for fixed values of the parameters R ) R(ss) and p ) p(ss). By A(x, R, p) denote the Jacobian of f in eq 6 with respect to the state variables x. If the real parts of the eigenvalues of A(x(ss), R(ss), p(ss)) are bounded above by a real number σ0 < 0, then for all > 0 there exists a δ > 0 such that for all initial conditions x(t0) ) x(0) with

||x(0) - x(ss)|| < δ the solution x(t) satisfies

for all t g t0. For a proof, the reader is referred to the book by Perko.20 Proposition 1 suggests to split the complex plane at z ) -|σ| i, for some σ ∈ R, σ > 0. See Figure 3a for an illustration. Similarly to bounding the real part above, both imaginary and real parts of the eigenvalues may be bounded by restricting their location to, for example, a trapezoidal sector of the complex plane or by a hyperbolic approximation of such a trapezoidal sector (cf. Figure 3b). The boundaries shown in Figure 3a,b can be described by a function

(17)

where m is smooth and maps from some subset of R2 into R, and σ and ω are the real and imaginary part of z ) σ + iω. Without loss of generality the function m can be assumed to be symmetric about the real axis, since the eigenvalues of A(x, R, p) are real or occur in complex conjugate pairs. Steady states of the process model (eq 6), for which an eigenvalue on the boundary defined by eq 17 exists, fulfill the system of equations

(

f)0

fxw(1) + ωw(2) - σw(1) fxw(2) - ωw(1) - σw(2)

)

˜f ) w(1)Tw(1) + w(2)Tw(2) - 1 ) 0 w(1)Tw(2) m(σ,ω)

m ) σ - σ0 The example of the hyperbolic sector sketched in Figure 3b can be described by

m)

x2 y2 -1 c12 c22

where c1 and c2 are real numbers that control the shape and location of the hyperbola. The augmented system (eq 18) is used to define the critical manifold of those steady states of the process model (eq 6) for which one or more eigenvalues lie on the boundary in the complex plane defined by eq 17. This critical manifold

M(p(0)) ) {(xT, x˜ T, RT)T ∈ U|f (x, x˜ , R, p(0)) ) 0, ˜f (x, x˜ , R, p(0)) ) 0} (19)

||x(t) - x(ss)|| e exp(-|σ0|(t - t0))

0 ) m(σ,ω)

into its real and imaginary part for the eigenvalue λ ) σ + i ω, σ, ω ∈ R and the corresponding eigenvector v ) w(1) + i w(2), where w(1), w(2) ∈ Rnx. The remaining equations in eq 18 normalize the eigenvector. The case of an upper bound σ < σ0 shown in Figure 3a is obtained for the choice

(18)

where f refers to the model equations in eq 6. The

has the same form as the critical manifolds for feasibility (eq 9). In eq 19, the auxiliary variable x˜ is defined (2) (1) (2) f denote the to be x˜ ) (w(1) 1 , ..., wnx , w1 , ..., wnx ), f and ˜ functions given in eq 18, and the linearization of f and ˜f with respect to x, x˜ , and R is assumed to have full rank. It is stressed that descriptions of the form of eq 18 and the corresponding critical manifolds of the form of eq 19 can be derived for several other types of critical points. The most prominent cases of critical points are saddle-node and Hopf bifurcations. These bifurcation points are of great practical importance because they constitute the stability boundary of a given process model.1,3 In fact, the analysis of saddle-node and Hopf bifurcations by numerical parameter continuation and applied bifurcation theory4 gave rise to the basic idea behind the approach presented here.3 The notion of a manifold of critical points is, however, more general in that it can be applied to critical points of higher codimension.3 Furthermore the concept of a critical manifold and the distance to such a manifold in the process parameter space can presumably be applied to critical points of the transient behavior as briefly discussed in the discussion of the HDA reaction section results. In the present paper, however, only the critical manifold (eq 19) is discussed in detail since this type of critical manifold is used in the application to the HDA reaction section. Details on other types of critical points along with applications in process optimization can be found elsewhere.2,3,21 Determination of Normal Vectors. Having identified a system of defining equations for the critical manifold, a system of equations for the calculation of the normal vector to this critical manifold must be

2742


derived. Mo¨nnigmann and Marquardt3 present a scheme for the derivation of normal vector systems for a general class of critical manifolds. This general class comprises the critical manifolds due to feasibility constraints (eq 9) as well as the critical manifold introduced in the previous section. When applied to the critical manifold for generalized constraints on eigenvalues (eq 19), the resulting normal vector system is

(

G(eig)(x, xj(eig), R, p, n) ) eq (18) fTx v(1) - ωv(2) - σv(1) + γ1w(1) - γ2w(2) fTx v(2) + ωv(1) - σv(2) + γ1w(2) + γ2w(1) v(1)Tw(1) + v(2)Tw(2) - mσ v

(1)T

fTx q

w

(2)

-v

(1)T

+v

n-

(2)T

(fTR q

fxxw

w

(1)

(1)T

+v

(1)

+ mω (2)T

+v

fxRw

(1)

fxxw

(2)

(2)T

+v

(2)

fxRw )

)

) 0 (20) Figure 4. More than one locally closest connection may exist due to the nonconvexity of critical manifolds or because more than one critical manifold exists.

where n is the desired normal vector and xj(eig) is short for the auxiliary variables (1) (2) (1) (2) (1) xj(eig) ) (w(1) 1 , ..., wnx , w1 , ..., wnx , ω, σ, v1 , ..., vnx , (2) T v(2) 1 , ..., vnx , γ1, γ2, q1, ..., qnx)

In eq 20, w ) w(1) + iw(2), w(1) ∈ Rnx, w(2) ∈ Rnx denotes the eigenvector of fx that corresponds to the eigenvalue σ + iω. Similarly, v ) v(1) + iv(2), v(1) ∈ Rnx, v(2) ∈ Rnx denotes the eigenvector of fTx that corresponds to the eigenvalue σ - iω. The symbols mσ and mω are short for the partial derivatives of the function m(σ, ω) with respect to σ and ω, respectively. Finally, q ∈ Rnx, γ1 ∈ R, and γ2 ∈ R are auxiliary variables that are necessary to state a square regular system.3 For brevity, the arguments (x, R, p) of the derivatives of f are omitted in (20). The notation G(eig)(x, xj(eig), R, p, n) ) 0 in eq 20 anticipates that normal vectors n can in general be calculated from systems of equations of the form

G(τ)(x, xj(τ), R, p, r) ) 0

(21)

where G is defined on some appropriate subset of, and maps into, Rnx × Rnxj × RnR × Rnp.3 The index τ denotes the type of critical manifold. The symbol xj(τ) denotes auxiliary variables that are necessary to calculate the normal vector n but are not used otherwise. Since these auxiliary variables depend on the type of the critical manifold, they are labeled with an upper index τ. Normal vector system of the form in eq 21 have been derived for various types of critical manifolds, among them feasibility constraints of the form in eq 8 and different types of bifurcation points of process models of the form eq 6. Normal vector systems can also be derived, however, for bifurcation points and singularities of higher order such as isolas and cusp bifurcations.3 These systems have been used in optimization problems similar to the one solved in the present paper. The reader is referred to previous works2,19,22,23 for applications to stability and feasibility boundaries, and for applications to higher codimension bifurcation points.3,21 Optimization with Normal Vector Constraints. In general multiple closest connections between Mr and critical manifolds Mc,i exist. Multiple closest connections

may arise due to the nonconvexity of a critical manifold, or because more than one critical manifold exists. These situations are sketched in Figure 4. Multiple closest connections are taken into account by stating a constraint of type (eq 15) for each locally closest connection. Let the upper index i ) 1, ..., imax enumerate the critical points that must be considered. The nonlinear program for the optimization of the process model (eq 6) with respect to the cost function φ then reads

φ(x(0), R(0), p(0))

(22)

s.t. 0 ) f (x(0), R(0), p(0))

(23)

min

x(0), R(0), p(0)

0 ) G(τi,i)(x(i), xj(τi,i), R(i), p(0), r(i)), ∀i ∈ I (24a) r(i) 0 ) R(0) - R(i) + l (i) , ∀i ∈ I ||r ||2

(24b)

0 e l(i) - xnR, ∀i ∈ I

(24c)

x ∈ X, R ∈ A, p ∈ P

(25)

where X, A, and P denote appropriate box constraints for the respective variables, and I abbreviates I ) {1, ..., imax}. Equation 23 ensures that the optimal point of operation is a steady state of the process model (eq 6). The set of constraints (eq 24) ensures the distance to the critical points. Note that n(i), x(i), xj(τi,i), R(i), and l(i) are degrees of freedom in the optimization in addition to x(0), R(0), and p(0). These variables are not listed as optimization variables in eq 22 in order to avoid a too tedious notation. The next subsection addresses the determination of the index set I ) {1, ..., imax}. Identification of Critical Points. In the optimization problem (eqs 22-25) the index set I is assumed to be known. This may be the case if an analysis of critical boundaries has already been conducted for the process of interest. In general, however, such an analysis is not available. In fact the time and effort spent on analyzing a process model becomes prohibitive for process models with many parameters, since the number of diagrams to be computed grows exponentially in the number of process parameters to be analyzed. Ideally a constructive method can therefore be applied to a process model for which no information on the


Figure 5. Detection of critical boundaries as the optimization proceeds. In panel a, the critical point (marked with a filled dot) is detected along the path from the upper right to the lower left corner. In panel b, the normal vector constraint has been added after the detection. Note that the detected critical point can in general not be connected to the center of the robustness ellipse along a normal direction. This is discussed in the text.

location and form of critical manifolds is available from an a priori analysis. In this case the index set I is unknown. The present section outlines a procedure for the optimization of process models for which an analysis is not available. Roughly speaking, the approach outlined below is a procedure for building up I ) {1, ..., imax} as the optimization of the process model proceeds. This procedure is illustrated with the aid of Figures 5 and 6. Detection of Critical Points along Paths. Assume that a process model is to be optimized for which critical boundaries exist that are, however, not known. For ease of discussion assume that these critical boundaries are stability boundaries. Since the critical boundaries are not known a priori, the index set I is empty to begin with. The optimization problem (eqs 22-25) therefore reduces to the steadystate optimization given by eqs 22, 23, and 25, without any additional constraints for parametric robustness (eq 24). Assume that this nonlinear program (NLP) is solved after setting the process variables x and parameters R, p to the values of a known, stable point of operation in order to provide the NLP solver with a feasible solution to start with. Figure 5a sketches a situation in which the optimizer starts at a stable point of operation and ends at an optimal but unstable point of operation. Along a path of feasible points, as sketched by the bold dashed line in Figure 5a, the stability boundary is crossed. The crossing of the stability boundary can be detected by monitoring the eigenvalues of the linearized process model. Alternatively, test functions used to

detect critical points in parameter continuation4,5 can be employed for several types of critical manifolds. Once a critical manifold has been found by detecting a critical point along the path in Figure 5a, the normal vector constraints for the corresponding critical manifold can be added to the NLP. This is sketched in Figure 5b. Note that the corresponding new normal vector constraint need not necessarily be active, i.e., the robustness ellipse need not touch the critical manifold Mc,1 (cf. Figure 5b). The normal vector constraints may be initialized for any point along the path for which the robustness ellipse entirely lies on the desired side of Mc,1. The index set I now comprises one element, and a new NLP can be solved which takes the location of the detected critical boundary into account (cf. Figure 6a). When this second NLP is solved, other critical points may be passed. The sequence of (i) optimization with best currently known approximation to the index set I, (ii) detection of critical boundaries crossed, and (iii) switching on of new normal vector constraints for the critical boundaries crossed, must be repeated until an optimal point is found which is at a safe distance to all critical manifolds (cf. Figure 6c). The detection of critical points along paths has been investigated in the context of numerical bifurcation analysis by continuation.4,5 The present section does not digress to discuss the test functions in more detail but only asserts that appropriate tests for critical points such as saddle-node and Hopf bifurcations and critical points of higher order such as cusps or isolas exist and can be implemented.4,5,24 In the application used here, eigenvalues are monitored explicitely to detect the crossing of this type of critical manifold. Before the algorithm that has just been sketched with the aid of Figures 5 and 6 is stated, several technical details deserve attention. Initializing Normal Vector Constraints. The NLP solver used to solve the optimization problem (eqs 2225) needs to be initialized with a solution that is feasible with respect to the constraints. In particular this implies that the normal vectors n(i) which connect the point of operation to the critical points must be known. The detection of critical points as described in the previous subsection, however, only provides a close critical point. This close critical point is in general not the point on the critical manifold which can be connected to the point of operation by a normal direction (cf. Figure 5b). For a given fixed candidate point of operation and a detected critical point, the locally closest critical point can be obtained from an optimization problem.

Figure 6. Independent or nonconvex critical manifolds may result in more than one normal vector constraint. A second critical boundary is detected in panel a. The corresponding normal vector constraint is switched on subsequently (b), until finally all critical boundaries are taken into account (c).

2744


crossed for some of the parameter values within the uncertainty ellipsoid because the critical manifold enters the robustness ellipse, but the nominal point does not cross the critical manifold. This problem can be solved with global solvers for nonlinear sets of equations. The current implementation uses interval Newton algorithms.25 Unlike algorithms that use floating point arithmetics, interval Newton solvers allow us to find all solutions for a set of nonlinear equations by a numerical calculation or to prove that no solution exists in a given domain.26 For simplicity, assume solutions to a set of nξ nonlinear equations Figure 7. Test functions fail to signal the crossing of a critical manifold if the critical manifold enters the robustness ellipse without passing the center.

In the situation sketched in Figure 5b, the closest point on the critical manifold to the candidate point of operation can be found by minimizing the distance ||R - R(0)||2 over all R ∈ Mc for fixed R(0), where Mc is the critical manifold. Minimizing ||R - R(0)||2 amounts to solving the NLP

min

x, x˜ , R

x(R - R(0))T(R - R(0))

s.t. 0 ) f (x, x˜ , R, p(0)) 0 ) ˜f (x, x˜ , R, p(0))

(26)

where R(0) and p(0) refer to the values at the fixed candidate point of operation. The constraints 0 ) f and 0 ) ˜f in eq 26 guarantee that R corresponds to the uncertain parameters of a steady state on the critical manifold (eq 19). It can be shown that the Lagrange multipliers needed to state the first-order KKT conditions of eq 26 can be used to initialize the normal vector n in constraints of the form in eq 15.2 For ease of reference, close critical points, in the vicinity of which locally closest critical points must be identified, are collected in a set denoted by J. Solving the optimization (eq 26) amounts to processing points in J to the set I. Switching Off Normal Vector Constraints. While methods for the detection of critical points and the subsequent turning on of normal vector constraints have been discussed, nothing has been said about switching off normal vector constraints. Situations may arise in which a stability boundary first imposes a constraint on the optimization, but later becomes inactive again in the sense that the distance to the corresponding critical manifold becomes larger than the required minimal distance for parametric robustness. A threshold ˆl is introduced, and normal vector constraints number i is switched off if

l(i) > ˆl,

xnR < ˆl < ∞

The threshold value ˆl is a tunable parameter of the algorithm. Its value depends on the particular problem at hand. Rigorous Search for Critical Points. A complication to the detection of critical points arises because the test functions are only evaluated at the nominal point of operation. Figure 7 sketches a situation in which the test function fails to signal that the critical manifold is

0 ) ψ(ξ), ψ: X f Rnξ are sought in a hypercube

ξ ∈ X ) X1 × X2 × ... × Xnξ ⊂ Rnξ where

X1 ) [ξlower , ξupper ], ..., Xnξ ) [ξnlower , ξnupper ] 1 1 ξ ξ are intervals in R. Interval arithmetics provide an inclusion for the functions on a given domain. Instead of calculating the value of a function ψi at a point ξ ∈ X by floating point arithmetics, interval arithmetics result in bounds on the function values on a given hypercube

(Iψ)(X) ) [ψlower, ψupper]

(27)

where (Iψ) denotes the interval extension of ψ. The interval extension of a function is created by replacing the elementary arithmetics operations (+, -, ×, /) and elementary functions by appropriate rules for the calculation with intervals. The extension of the elementary addition of two real numbers is replaced by the addition of two intervals,

, ξupper ] + [ξ˜ lower , ξ˜ upper ]) [ξlower i i i i + ξ˜ lower , ξupper + ξ˜ upper ] [ξlower i i i i , ξupper , ξ˜ lower , and ξ˜ upper are for example, where ξlower i i i i real numbers. Just as a floating point instance of a nonlinear function can be programmed in terms of the elementary floating point operations and intrinsic functions of a programming language, interval extensions can be built up from the interval extensions of the elementary arithmetics operations and intrinsic functions. Several implementations are available.26 In the search for solutions to the nonlinear systems 0 ) ψ(ξ) on X, the hypercube X is subdivided into smaller hypercubes X(k) ⊂ X until the existence or nonexistence of a solution can be shown for an X(k). If the inclusion of ψ on X(k) yields

0 ∉ X(k) i for any component i of the interval, then the hypercube X(k) cannot contain any solution. Otherwise X(k) needs to be subdivided into smaller parts and the inclusion test must be repeated for these smaller parts. For sufficiently small parts X(k) ⊂ X, for which the inclusion does not imply that no solution exists, the existence of a solution in X(k) can be verified with automated fixed point theorems.26 Since interval arithmetics can be used to show that no solution to a set of nonlinear equations exists on a


given hypercube, they can be used to prove that no critical points exist in the robustness hypercube (eq 7). On the other hand, if critical points exists in the hypercube of interest, the interval arithmetics solver can be used to identify one or more of these critical points as a solution. These critical points can then be put into the set J and processed to I. Since interval Newton solvers are computationally demanding, their use becomes prohibitive for large search domains and large numbers of equations. It is stressed here, however, that these solvers are only used to search for critical points in the robustness region. This robustness region typically is a very small part of the process parameter space only. Furthermore, the rigorous search with interval arithmetics is only conducted at a late stage of the algorithm given below. Nevertheless, the rigorous search with interval Newton solvers currently constitutes the bottleneck for an application of the algorithm. Since the normal vector constraints can be applied to much larger problems than the rigorous search, alternatives to a check for critical points based on interval arithmetics must be sought. While no approaches exist that provide the same mathematical rigor as interval arithmetics, workarounds are available. These workarounds are discussed at the end of the next section. Algorithm Having stated the necessary building blocks, the algorithm can now be outlined as follows. In a narrow sense the term critical manifold refers to the type of manifold (eq 19). More generally, the optimization algorithm can be applied to process models for which any combination of feasibility boundaries, stability boundaries, boundaries due to generalized constraints on eigenvalues, or critical manifolds due to bifurcation points of higher order such as cusps or nontransversal Hopf points exists.2 The term critical point refers to a point on any such manifold. 1. Initialization: Choose an equilibrium point that is feasible and has the desired nonlinear dynamics properties. If close critical points are known a priori, put them into J. Choose a value for ˆl. 2. Update of I: For any point in J solve eq 26 to find a locally closest point and the corresponding normal vector that connects the candidate point of operation and the critical manifold. Add the result to the index set I and remove the point from J. Remove those normal vector constraints from the index set I for which l(i) > ˆl. 3. Optimization: Solve problem eqs 22-25 using the current index set I. 4. Analysis: Analyze the path between the starting and the end point of step 3. If a point on the path exists at which a critical manifold is crossed, put it into J and return to step 2. Otherwise proceed to the rigorous search. 5. Rigorous Search: Check for critical points in the robustness region by running an interval arithmeticsbased rigorous check on the defining system of equations for any type of critical manifold to be taken into account for the process at hand. If critical points exist in the robustness region that have not been detected in step 4, put them into J and return to step 2. If no critical points are found, stop. An optimal steady state has been identified which is parametrically robust with respect to the feasibility constraints and with

respect to the desired nonlinear dynamics properties. The algorithm is therefore terminated. Several details of the implementation deserve discussion. These details are treated in the following subsections. Implementational Issues Linear Path versus Feasible Path Detection. The test function for critical points are only defined at steady-state solutions of the model (eq 6).5 Similarly, it is only meaningful to calculate eigenvalues for the Jacobian evaluated at a steady state of the model. The detection of critical points along a path can therefore only be carried out if the optimizer moves along a path of feasible points. While feasible path optimizers for nonlinear programs exist,27,28 this is in fact an unnecessarily strong restriction. In fact all problems which the algorithm has been applied to so far1-3,19,21-23 as well as the large-scale example presented here have successfully been optimized without a feasible path optimizer. Instead of forcing the optimizer to proceed along such a path, the starting point of the optimization and the candidate optimal point to which the optimizer proceeds are connected by a linear path in the process parameter space. The test functions are then evaluated along this linear path. Critical points detected along this linear path are put into J and processed to I by solving the optimization problem discussed in the section on normal vector constraint initialization. Clearly, this is only a heuristic shortcut to feasible path optimization. Since this test is cheaper, however, it can be used first. If it fails, one must retreat to feasible path optimization. Practical Test for Overlooked Critical Manifolds. It has been pointed out that the rigorous search becomes prohibitive for large search domains and a large number of nonlinear equations. Step 5 therefore turns out to be the bottleneck for the application of the algorithm to large-scale examples. Since the normal vector constraints are anticipated to work for much larger examples than tractable with interval arithmetics, the search for critical points in the robustness region in step 5 must be conducted with local methods. Recall that the rigorous test for critical points is necessary because the test functions are only evaluated at the nominal point of the robustness ellipsoid. If a critical manifold enters the robustness region without passing the nominal point, the test function evaluation fails to signal the existence of critical points in the robustness ellipsoid. The obvious idea is to evaluate the test functions for a sample of points in the robustness ellipsoid. A natural choice for sample points are the corner points of the robustness hypercube (eq 7) or the centers of the faces of the hyperplanes that bound the hypercube. For the large-scale example presented below, the latter set of sample points was used. Generation of Higher Order Derivatives. In the particular normal vector system used in the present paper (eq 20) as well as in the normal vector systems in general,3 higher order derivatives of the model equations with respect to the state variables and the uncertain parameters R occur. Derivatives are currently calculated by a combination of symbolic and automatic differentiation with MAPLE29 and ADIFOR.30 First, the derivatives of the model equations that occur in the normal vector system (eq 20) are calculated with a combination of symbolic and automatic differentiation by MAPLE. The current implementation automatically

2746


Figure 8. Flow sheet of the reaction section of the HDA process. The splitter and the purge are modeled in terms of algebraic equations only. Dynamical models are used for all other units. The compressor is not part of the dynamical process model but only enters the cost function. The cooler between the hot side of the heat exchanger and the flash is part of the flash model. It is drawn separately for ease of presentation only.

unravels common subexpressions in the functions and their derivatives. Intermediate variables are then introduced to make sure that common subexpressions are only evaluated once at runtime. Similarly, algebraic equations that can be solved for algebraic variables globally are solved symbolically. The resulting explicit equations for algebraic variables are treated as common subexpressions that can be calculated once and used repeatedly in subsequent algebraic or dynamic equations. This model structure is superior to a model that results after substituting explicit algebraic equations into the ODEs for two reasons. For one, common subexpressions are only evaluated once. More importantly, however, this structure is exploited when generating derivatives by automatic differentiation as the sparsity structure of matrices of first order and tensors of second-order derivatives can be taken into account. Fortran code is generated for NLPs of the form in eqs 22-25 using MAPLES’s codegen device. Finally, the derivatives needed by the NLP solver are generated by applying ADIFOR to this fortran code. The NLPs are currently solved with NPSOL.31 In the current implementation, the entire optimization can be controlled from within a MAPLE sheet. Code generated by MAPLE’s codegen device as well as third party numerical code such as NPSOL or PITCON for parameter continuation is interfaced to the MAPLE sheet with MAPLE’s define_external command.

Application to Reaction Section of HDA Process The application presented in this section deals with the reaction section of a process for the production of benzene. Since benzene is produced by the hydrodealkylation of toluene, the process is usually referred to as the HDA process in the literature. Process Model. In the HDA process, benzene is produced from the raw materials hydrogen and toluene. Methane and diphenyl arise as byproducts. The two gasphase reactions are r1

toluene + H2 98 benzene + methane r2

z diphenyl + H2 2benzene y\ r 3

(28)

Figure 8 shows a flowsheet of the reaction section of the HDA process. The hydrogen stream is mixed with the recycle stream in the mixer, and the resulting stream is preheated in the feed effluent heat exchanger. The fresh toluene stream is fed to the cold side of the heat exchanger. A furnace is used to further heat the cold side exit stream of the heat exchanger to the desired reactor inlet temperature. The reaction takes place in a plug flow reactor that is modeled as a series of continuous stirred tank reactors. The effluent of the reactor is then cooled to prevent coking in the hot side of the heat exchanger. Most of the hydrogen and methane is separated from the aromatics contained in


the reactor effluent in a flash. The liquid stream leaving the flash is further processed in a separation section that is not modeled in detail in the present work. Here it is assumed that the desired benzene can be separated from the liquid stream leaving the flash in a splitter and that the unreacted toluene can be separated and fed back to the reaction section. The hydrogen present in the flash exit vapor stream is recycled to the reaction section as well. The methane leaving the flash in the vapor stream must be considered an undesired impurity. A purge stream is needed to avoid accumulation of methane. Douglas16 lists process conditions to be maintained in the HDA process. From these conditions constraints for process optimization can be inferred as follows. The reactions take place in the temperature range between about 894 K (1150 °F) and 978 K (1300 °F). Below the given temperature range, the reaction rate is too low, and above the upper bound, hydrocracking becomes significant.16 These temperature ranges impose the bounds

894 K e TReactor,i e 975 K ∀i

(29)

where TReactor,i denotes the exit temperature of lump number i. The lower bounds are in fact insignificant since the reactor feed stream is heated to the temperature 894 K, the reactor is operated adiabatically, and the main reaction is exothermic. The requirement of an inlet temperature of 894 K to the reactor establishes the setpoint of the furnace temperature controller:

TFurnace,setpoint ) 894 K

(30)

According to Douglas,16 the reaction is carried out at a pressure of 30 × 105Pa ) 30 bar. This requirement implies the value

pMixer,setpoint ) 30 bar

(31)

for the setpoint of the mixer pressure controller. A 5:1 ratio of hydrogen to aromatics at the reactor inlet is needed to prevent coking.16 This imposes the constraint

yH,0 g5 yB,0 + yT,0 + yD,0

(32)

where yi,0 denotes the molar fractions of the components at the reactor inlet. The indices H, B, T, and D are short for hydrogen, benzene, toluene, and diphenyl, respectively. The index M denotes methane in the sequel. The fresh H2 feed is assumed to contain 5% methane as an impurity.16 To prevent coking in the heat exchanger, the reactor effluent must be cooled to 894 K (1150 °F).16 This determines the setpoint of the temperature controller in the cooler to be

TCooler,setpoint ) 894 K

(33)

Finally, the production rate of benzene is required to be 2 kmol/min according to Douglas.16 This imposes the constraint

xFlash,B LFlash ) 2

kmol min

(34)

where xFlash,B and LFlash denote the concentration of benzene in the liquid stream leaving the flash and the flow rate of this stream, respectively. In constraint eq 34, it is again assumed that practically all benzene can be separated from the hydrogen, methane, toluene, and diphenyl. The dynamical models of the units are discussed in the Appendix. While simulations of the HDA process have been carried out before, a dynamical model is not available in the literature. Luyben and Tyreus32 use a model implemented in an in-house simulator at DuPont. The authors give a rough description of the model only but do not present the model equations. Cao et al.33 use a model implemented in SPEEDUP. According to the authors the model is not available in a form that allows it to be used with currently available simulation software.34 The model used for the optimization presented here therefore had to be implemented from scratch.2 Some important details of the model are summarized here: • Dynamical models with material and energy balances are used for the mixer, heat exchanger, furnace, reactor, cooler, and flash. The purge and the splitter are modeled in terms of algebraic equations that capture the conservation of mass. • Rigorous mass and energy balances for the reactor and the heat exchanger would result in partial differential equations. Here the reactor is modeled as a series of five continuous stirred tank reactors. Similarly, the heat exchanger is modeled as five continuous stirred tank reactors each on the hot and cold side. The numbers of lumps were chosen to obtain reasonable computation times and memory requirements in the optimization runs. • The separation in the flash is modeled with constant relative volatilities. All other units contain gas only. • The separation section of the HDA process is currently not modeled. The algebraic splitter model assumes that the toluene to be recycled can be separated completely from the liquid stream leaving the flash. Likewise, it is assumed that the benzene contained in the liquid stream leaving the flash can be separated completely from the remaining diphenyl and methane. • The process model comprises five PI controllers as indicated in the process flowsheet in Figure 8. The optimization problem solved in the next section integrates process and controller design in the sense that the normal vector-based approach allows to optimize the process and tune the controllers simultaneously. • The algebraic process model equations can be solved for the algebraic variables explicitly in this particular example. This implies that a solution of the type in eq 4 exists not only locally but globally. The process is optimized with respect to the total annualized cost

φTAC ) (capital cost mixer + capital cost heat exchanger + capital cost furnace + capital cost reactor + capital cost cooler + capital cost heat exchanger + capital cost flash)/Tpayback + operating cost furnace + operating cost compressor + cost toluene + cost hydrogen (35) where Tpayback refers to a payback period, which was

2748


Table 1. Uncertain Parameters r of the HDA Reaction Section uncertain parameter Ri reaction rate (main reaction) heat of reaction (main reaction) fresh toluene feed rate fresh toluene feed temperature fresh hydrogen feed temperature heat capacity of FEHE metal FEHE heat transfer coefficient relative volatility hydrogen relative volatility methane relative volatility benzene relative volatility toluene relative volatility diphenyl

nominal value

Table 2. Parameters p Not Affected by Uncertaintya

uncertainty (5% (1% (0.01 kmol min-1 (3 K (3 K (10% (10% (10% (10% (10% (10% (10%

bar3/2

Table 3. Flow Rates, Temperatures, and Molar Fractions at the Optimum for Some Unitsa

reactor volume heat exchanger area purge stream flow rate mixer output flow rate kTC,Cooler, τTC,Cooler kTC,Furnace, τTC,Furnace kTC,Flash, τTC,Flash kPC,Mixer, τPC,Mixer kLC,Flash, τTC,Flash a The indices TC, PC, and LC are short for temperature, pressure, and level control, respectively. The parameters ki and τi denote the proportional and integral time constants of the respective controller.

chosen to be 3 yr. The contributions to the cost function are detailed in the Appendix. Optimization. In this section, the algorithm for the optimization with normal vector constraints is applied to the reaction section of the HDA process. The generalized eigenvalue normal vector system (eq 20) is used to find an optimal point of operation for which a userspecified decay rate to linear order can be guaranteed despite parametric uncertainty. This amounts to setting the function m in the generalized eigenvalue system (eq 20) to

m ) σ - σ0

min-1

3.261 × kmol 49976 kJ kmol-1 subject to optimization 303.15 K 303.15 K 0.46 kJ kmol-1 K-1 17.21 kJ m-2 min-1 K-1 1.0 0.13575 1.1665 × 10-4 4.003 × 10-5 2.0488 × 10-5 109

(36)

The bound on the decay rate is chosen to be σ0 ) -1/30 min-1. Twelve parameters are considered to be uncertain in the optimization. Table 1 lists these uncertain parameters R. Note that two different types of uncertain parameters exist. All parameters but the fresh toluene feed rate are uncertain but have fixed nominal values. The nominal value of the fresh toluene feed rate, in contrast, is subject to optimization. It must be stressed that the uncertainties given in Table 1 are not motivated by a real process, but they are merely reasonable values that are used to demonstrate the applicability of the developed approach to a large-scale model. The parameters p not subject to uncertainty are given in Table 2. Note that the mixer output flow rate is considered a degree of freedom in the optimization. Increasing or decreasing this flow rate causes the pressure in the mixer to change. The pressure changes in turn cause the mixer pressure controller to adjust the fresh hydrogen flow rate to larger or smaller values accordingly. Freeing the mixer output flow rate therefore amounts to indirectly considering the fresh hydrogen flow rate to be a degree of freedom. The fresh hydrogen stream cannot be considered to be a free variable directly since it is a manipulated variable. Since an analysis of the dynamical behavior of the HDA reaction section model is not available, the location

flow rate temperature H2 methane benzene toluene diphenyl




mixer gas recycle

mixer fresh H2

mixer out

25.324 318.15 0.43186 0.56672 0.12684e-2 0.14677e-3 0.78119e-7

3.5108 303.15 0.95000 0.5000e-1 0.0 0.0 0.0

28.835 316.52 0.49495 0.50381 0.11140e-2 0.12890e-3 0.687e-7

FEHE cold in

FEHE recycled toluene

FEHE fresh toluene

28.835 316.52 0.49495 0.50381 0.11140e--2 0.12890e--3 0.68700e--7

0.67441 318.15 0.0 0.0 0.0 1. 0.0

2.145 303.15 0.0 0.0 0.0 1.0 0.0

FEHE cold out

FEHE hot in

FEHE hot out

31.654 717.01 0.45088 0.45895 0.10148e-2 0.89160e-1 0.625e-7

31.654 893.61 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

31.654 514.47 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

furnace in

furnace out

31.654 717.01 0.45088 0.45895 0.10148e-2 0.89160e-1 0.625e-7

31.654 893.61 0.45088 0.45895 0.10148e-2 0.89160e-1 0.625e-7

a See Table 4 for information on the other units. Flow rates and temperatures are given in kmol min- 1 and K, respectively. The abbreviation FEHE is short for feed effluent heat exchanger.

of the critical boundaries due to eq 36 is unknown. The sets I and J are therefore empty in step 1 of the algorithm. The starting point chosen in step 1 has a dominant decay rate of -1/26.23 min-1 and therefore complies with the constraint in eq 36. The value of ˆl is chosen to be ˆl ) 2xnR. In fact this value turns out be insignificant since no normal vector constraints ever have to be switched off. Since I and J are empty, step 2 can be skipped. The optimization is therefore carried out with an empty index set I in step 3 of the algorithm, and the crossing of critical boundaries is detected as in step 4. The

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2749 Table 4. Flow Rates, Temperatures, and Molar Fractions of Some Units at the Optimuma flow rate temperature H2 methane benzene toluene diphenyl


flow rate temperature H2 methane benzene toluene diphenyl furnace power FEHE power cooler power reactor volume heat exchanger volume heat exchanger area

reactor in

reactor out

cooler in

cooler out

31.654 893.61 0.45088 0.45895 0.10148e-2 0.89160e-1 0.625e-7

31.654 943.74 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

31.654 943.74 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

31.654 893.61 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

flash in

flash liquid out

flash vapor out

31.654 514.47 0.38537 0.52667 0.64308e-1 0.21436e-1 0.22156e-2

3.5918 318.15 0.22114e-1 0.21378 0.55682 0.18776 0.19525e-1

28.062 318.15 0.43186 0.56672 0.12684e-2 0.14677e-3 0.78119e-7

purge in

purge purged

purge out

28.062 318.15 0.43186 0.56672 0.12684e-2 0.14677e-3 0.78119e-7 5.938363 MW 1.610605 MW -1.781581 MW 127.203774 m3 4.87819 m3 (same on hot and cold side) 957.430727 m2

2.7375 318.15 0.43186 0.56672 0.12684e-2 0.14677e-3 0.78119e-7

25.324 318.15 0.43186 0.56672 0.12684e-2 0.14677e-3 0.78119e-7

a Some additional data is given on unit sizes and energy consumption. For the units not listed here, see Table 3. Flow rates and temperatures are given in kmol min-1 and K, respectively. The abbreviation FEHE is short for feed effluent heat exchanger.

optimization with an empty index set I started at the feasible point chosen in step 1 results in an optimal point of operation with a dominant decay rate of -1/48.4 min-1 and a value of $18.9 million/year for the cost function. Clearly, the critical boundary induced by eq 36 for σ j ) -1/30 min-1 must have been crossed. The evaluation of the leading eigenvalue along the linear connection between the starting and the end point of the optimization allows us to identify a critical point that is put into the set J, and the algorithm returns to step 2. The next pass through the algorithm starts by processing the point in J to I in step 2. The optimization of the process model is then repeated with the generalized eigenvalue constraints for this new point in I in step 3. The result is a steady state with a dominant decay rate of -1/27.4 min-1 and a value of the cost function (eq 35) of $20.1 million/year. In step 4, the evaluation of the leading eigenvalue along the path between the starting point and end point of the optimization shows that no critical boundary has been crossed. The algorithm therefore proceeds to step 5. As anticipated in the section on the algorithm, the HDA process model is too large for a search of critical points in the robustness ellipse with interval arithmetics. The existence of critical points in the robustness ellipse is therefore tested for by continuation. Starting at the candidate optimal point, a continuation to both the upper bound (Ri + ∆Ri) and the lower bound (Ri - ∆Ri) is conducted for each of the uncertain parameters listed in Table 1. Along the obtained curves, the leading eigenvalue is calculated. For all curves, the leading eigenvalue does not attain nor exceed the user-specified value of -1/30 min-1. From this result the robustness of the obtained optimum is inferred, and the algorithm

is terminated. The obtained optimum is documented in Tables 3 and 4. It is further analyzed with the aid of step responses in the next subsection. Analysis of the Result. In step 5 of the proposed algorithm, the eigenvalues of the candidate optimal points and of points in the vicinity are evaluated. These checks reveal whether the user-specified dynamic behavior (eq 36) can be guaranteed despite the presence of the parametric uncertainties given in Table 1. While the generalized eigenvalue constraints for the userspecified dynamic behavior (eq 36) can be used to guarantee the desired decay rate to linear order, the actual nonlinear response of the process (i.e., the trajectory to steady state after a disturbance) in general deviates considerably from that of the linearized process. Apart from the checks of the eigenvalues conducted in the previous section, the optimal point of operation is therefore further analyzed by simulation after disturbances. If losses due to the purge stream and the side reaction of benzene to diphenyl are neglected for the moment, the benzene production rate is set by the fresh toluene feed rate (cf. eq 28). The toluene feed rate therefore is the natural input variable to control the benzene production rate. Since changes in the benzene production rate are of interest from a practical point of view, step responses to changes in the fresh toluene feed rate are investigated. Figure 9 shows the step response to an increase of the toluene feed rate by 10%. These results were obtained by setting the initial conditions to the values at the optimum as summarized in Tables 3 and 4. The step change is employed at time t ) 0. Since the decay rate is chosen to be -1/30 min-1, the amplitude of the disturbance due to the step change must have decayed

2750


Figure 9. Step response to an increase of the fresh toluene feed rate by 10%. The symbols Qi denote the power transferred to (Qi > 0) or from (Qi < 0) the respective units. The symbols LLC,Flash and FPC,Mixer denote the flow rates of the liquid output stream leaving the flash and the fresh H2 inlet stream to the mixer, respectively. The dashed boxes indicate the bounds of a decay to 1/e. See the text for a discussion.

to within exp(-1) ) 36.8% after 30 min. The dashed boxes shown in Figure 9 indicate the corresponding upper and lower bounds on the respective variable which must not be exceeded. When calculating these bounds, the new steady-state values are approximated by the values reached at t ) 180 min. From the results shown in Figure 9, it can be inferred that the user-specified bound on the performance is met by the process. Note that the benzene production increases by slightly less than 10% as anticipated in the brief discussion of the benzene losses above. Figure 10 shows the step response to a decrease of the toluene feed rate by 10%. These results were obtained in the same way as the ones shown in Figure 10. Again it can be inferred that the user-specified bound on the performance is met by the process. In this case, the benzene production rate drops by approximately 10% as expected. Summary and Future Perspectives In this paper, a new approach to the optimizationbased design of steady-state processes subject to parametric uncertainty was presented and applied to a sample process model. The new approach is based on the recently reported idea1-3,17 of measuring the dis-

tance between a candidate point of operation and the critical manifolds in its vicinity with the normal vector to the critical manifolds. By imposing a lower bound on the distances to the critical manifolds, parametric robustness can be ensured in process optimization. Based on these ideas, an algorithm for the optimization under parametric uncertainty has been developed that is presented and applied in the present paper for the first time. The notion of a critical manifold and the notion of the distance between a candidate point of operation to a critical manifold can be applied to stability boundaries, other boundaries on the dynamical properties, and feasibility boundaries alike. Due to the generality of the concept of a critical manifold, parametric robustness can be guaranteed with respect to a variety of system properties. While examples for the versatility of the new approach have been presented elsewhere,1-3,19,21,22 the aim of the present paper is to demonstrate that the new approach can be applied to dynamical systems with (i) several hundred equations for which (ii) no a priori knowledge on the shape and location of the existing critical manifolds is available. In fact an analysis of the treated process model, the reaction section of the HDA process,16 would have been tedious since 12 uncertain parameters would have had to be analyzed. While in


Figure 10. Step response to a decrease of the fresh toluene feed rate by 10%. The same variables as in Figure 9 are shown. The dashed boxes indicate the bounds of a decay to 1/e. See the text for a discussion.

principle the information on the location of the critical manifolds can be disclosed by an analysis of the model by parameter continuation and bifurcation theory, such an analysis becomes cumbersome for more than four or five parameters. The example therefore demonstrates that the critical manifold-based approach can be applied to a number of parameters that can hardly be addressed by continuation analysis. The tractable system size is currently limited by high runtime memory requirements that are caused by treating the emerging matrices of derivatives as dense matrices. The use of sparse derivatives is anticipated to reduce the memory requirements of the critical manifold-based approach considerably. This will allow the application of the critical manifold-based approach developed in this thesis to models even larger than a few hundred equations. In particular, the complete HDA process including a detailed description of the separation section can be expected to be tractable. The presented approach can currently only be applied to steady-state process optimization. Moreover, disturbances, inputs and reference signals to the process systems must be presumed to vary only quasistatically with respect to the system dynamics. This clearly is a strong assumption and a drawback of the method. Future work will deal with relaxing these prerequisites. In fact, the concepts of a critical manifold and the distance to the critical manifold are likely to be ap-

plicable to critical points for the transient behavior.2 This idea is currently being investigated. Appendix Cost Function. The HDA process is optimized with respect to the total annualized cost for a fixed production rate. The total annualized cost is the sum of the annual capital cost of the mixer, heat exchanger, furnace, reactor, cooler, flash, and compressor. In addition, the operating cost due to the fuel needed in the furnace, and the operating cost of the compressor is accounted for. The annual capital cost is assumed to be the capital costs of the units divided by a payback period of 3 years. The cost functions for the unit operations are taken from Reyes and Luyben35 and Douglas.16 Marshal and Swift indices for chemical engineering units as of the year 2002 are used to update the cost correlations.36 The index values are IMS ) 697.8 for the compressor and IMS ) 353.8 for all other units. Capital Cost of Reactor, Mixer, and Flash. According to Douglas,16 the installed cost of the vessels of these units can be approximated by

5485 d1.066 h0.802(Fm + Fp) [$]

(37)

where d and h denote the diameter and height in meters, respectively. Note that the diameter and the

2752


height of the units could be considered independent optimization variables. Since the cost function provide rough estimates of the actual cost according to Douglas16 only, however, the diameter-to-height ratio is set to 1. This amounts to considering the volume to be a free variable only. The coefficients Fm and Fp account for the material the vessel is build from and for the pressure. They are set to Fm ) 1 for steel and Fp ) 1.45 for pressures of about 30 bar. Capital Cost of the Feed-Effluent Heat Exchanger and the Cooler. The cost of the heat exchangers is a function of the heat transfer area. According to Douglas,16 the cost can be estimated by

c ) 2285 A0.65 [$]

(38)

where A is the heat exchanger area in square meters. While the feed-effluent heat exchanger does not have an operating cost, the cooler introduces an operating cost due to the necessary cooling water. A steady-state calculation assuming a heat transfer coefficient value k ) 100 W m-2 K-1,37 a cooling water temperature of T ) 25°K, and a cooling water price of $0.01/m3 shows that the cost of cooling water would be a few hundred dollars per year and thus can be neglected. Operating and Capital Cost of the Furnace. Following Reyes and Luyben,35 the cost of fuel is assumed to be $5/Btu ) 0.4742 × 10-5 $/kJ. This amounts to an operating cost of

c1 ) 2.499 Q $ min-1 where Q is measured in kJ min-1. According to Douglas,16 the installed cost of the furnace is a function of the heat to be transferred

c2 ) 3.904 Q0.85 [$]

(39)

Operating and Capital Cost of the Compressor. According to Douglas,16 the operating and installed cost of the compressor can be stated as functions of compressor power. They read

c1 ) 13.66P [$/year] and

c2 ) 177.7P0.82 [$]

(40)

respectively, where P denotes the compressor power. The compressor power is calculated from the equation for reversible adiabatic compression of an ideal gas35 according to

P)

(( )

γRT pout γ - 1 pin

(γ-1)/γ

)

-1

(41)

where γ ) cp/cv. Literature Cited (1) Mo¨nnigmann, M.; Marquardt, W. Steady-state process optimization with guaranteed robust stability and feasibility. AIChE J. 2003, 49 (12), 3110-3126. (2) Mo¨nnigmann, M. Constructive Nonlinear Dynamics Methods for the Design of Chemical Engineering Processes. Ph.D. Thesis, Aachen University, 2003.

(3) Mo¨nnigmann, M.; Marquardt, W. Normal vectors on manifolds of critical points for parametric robustness of equilibrium solutions of ODE systems. J. Nonlinear Sci. 2002, 12, 85-112. (4) Kuznetsov, Y. A. Elements of Applied Bifurcation Theory, 2nd ed.; Springer-Verlag: New York, 1999. (5) Beyn, W.-J.; Champneys, A.; Doedel, E.; Govaerts, W.; Kuznetsov, Y. A.; Sanstede, B. Numerical continuation, and computation of normal forms. In Handbook of Dynamical Systems, Vol. 2; Fiedler, B., Ed.; North Holland/Elsevier: Amsterdam, 2002; pp 149-219. (6) Mangold, M.; Kienle, A.; Gilles, E. D.; Mohl, K. D. Nonlinear computation in DIVAsmethods and applications. Chem. Eng. Sci. 2000, 55 (2), 441-454. (7) Kuznetsov, Y. A. CONTENTsIntegrated environment for analysis of dynamical systems. Tutorial; Dynamical Systems Laboratory, Centrum voor Wiskunde en Informatica, Amsterdam, 1998; ftp.cwi.nl/pub/CONTENT. (8) Doedel, E.; Champneys, A. R.; Fairgrieve, T. F.; Kuznetsov, Y. A.; Sandstede, B.; Wang, X.-J. AUTO97: Continuation and bifurcation software for ordinary differential equations (with HomCont), 1997; Computer Science, Concordia University, Montreal, Canada; ftp.cs.concordia.ca/pub/doedel/auto. (9) Bildea, C. S.; Dimian, A. C.; Iedema, P. D. Multiplicity and stability approach to the design of heat-integrated multibed plug flow reactor. Comput. Chem. Eng. 2001, 25, 41-48. (10) Bildea, C. S.; Dimian, A. C.; Iedema, P. D. Nonlinear behavior of reactor-separator-recycle systems. Comput. Chem. Eng. 2000, 24, 209-215. (11) Bildea, C. S.; Dimian, A. C. Singularity theory approach to ideal binary distillation. AIChE J. 1999, 45 (12), 2662-2666. (12) Bildea, C. S.; Dimian, A. C. Stability and multiplicity approach to the design of heat-integrated PFR. AIChE J. 1998, 44, 2703-2712. (13) Khinast, J.; Luss, D.; Harold, M.; Ostermaier, J.; McGill, R. Continuously stirred decanting reactor: Operability and stability considerations. AIChE J. 1998, 44 (2), 372-387. (14) Harold, M. P.; Ostermaier, J. J.; Drew, D. W.; Lerou, J. J.; Daniel, L., Jr. The continuously-stirred decanting reactor: Steady state and dynamic features. Chem. Eng. Sci. 1996, 51, 1777-1786. (15) Mo¨nnigmann, M.; Hahn, J.; Marquardt, W. Towards constructive nonlinear dynamicsscase studies in chemical process design. In Nonlinear Dynamics of Production Systems; Radons, G., Neugebauer, R., Eds.; Wiley-VCH: Weinheim, 2003. (16) Douglas, J. M. Conceptual Design of Chemical Processes, 1st ed.; McGraw-Hill Chemical Engineering Series; McGraw-Hill: New York, 1988. (17) Dobson, I. Computing a closest bifurcation instability in multidimensional parameter space. J. Nonlinear Sci. 1993, 3, 307-327. (18) Kremer, G. G.; Thompson, D. F. A bifurcation-based procedure for designing and analysing robustly stable non-linear hydraulic servo systems. Proc. Inst. Mech. Eng. Part I 1998, 212 (15), 383-394. (19) Mo¨nnigmann, M.; Marquardt, W. Parametrically Robust Control-Integrated Design of Nonlinear Systems. In Proceedings of American Control Conference, Anchorage, Alaska, USA, Vol. 6; 2002; pp 4321-4326. (20) Perko, L. Differential Equations and Dynamical Systems, 2nd ed.; Springer-Verlag: New York, 1996. (21) Gerhard, J.; Mo¨nnigmann, M.; Marquardt, W. Robust stable nolinear control and design of a CSTR in a large operating range. Accepted for presentation at DYCOPS 7, Cambridge, MA, 2004. (22) Grosch, R.; Mo¨nnigmann, M.; Briesen, H.; Marquardt, W. Optimal design of a continuous crystallizer with guaranteed parametric robust stability. In European Symposium on Computer Aided Process Engineering-12; Elsevier: Amsterdam, 2002; pp 205-210. (23) Grosch, R.; Mo¨nnigmann, M.; Marquardt, W. Integrated robust optimal design and control of an MSMPR crystallization process. Presented at the AIChE Annual Meeting, San Franscisco, CA, November 16-21, 2003. (24) Golubitsky, M.; Schaeffer, D. Singularities and Groups in Bifurcation Theory, Volume I; Springer-Verlag: New York, 1985. (25) Bischof, C. H.; Lang, B.; Marquardt, W.; Mo¨nnigmann, M. Verified determination of singularitites in chemical processes. Proceedings of SCAN 2000, 9th GAMM-IMACS Interational

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2753 Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics, September 18-22, 2000, Karlsruhe, Germany. (26) Kearfott, R. B. Rigorous Global Search: Continuous Problems; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (27) Rawlings, J. B. Private communication to W. Marquardt, 2003. (28) Wright, S. J.; Tenny, M. J. A feasible trust-region sequential quadratic programming algorithm. Technical report, Computer Science Dept., University of Wisconsin, 2002. (29) Monagan, M. B.; Geddes, K. O.; Heal, K. M.; Labahn, G.; Vorkoetter, S. M.; McCarron, J. Maple6 Programming Guide; Waterloo Maple Inc.: Waterloo, ON, Canada, 2000. (30) Bischof, C.; Carle, A.; Hovland, P.; Khademi, P.; Mauer, A. ADIFOR 2.0 User’s Guide (Revision D). Technical Report 192/ CRPC-95516-S, Mathematics and Computer Science Division Technical Memorandum and Center for Research on Parallel Computation, 1998. (31) Gill, P.; Murray, W.; Saunders, M.; Wright, M. User’s Guide for NPSOL, Version 4.0; Systems Optimization Laboratory, Stanford University: Stanford, CA, 1986.

(32) Luyben, M. L.; Tyreus, B. D. An industrial design/control study for the vinyl acetate monomer process. Comput. Chem. Eng. 1998, 22 (7-8), 867-877. (33) Cao, Y.; Rossiter, D.; Edwards, D. W.; Knechtel, J.; Owens, D. Modelling issues for control structure selection in a chemical process. Comput. Chem. Eng. 1998, 22, 411-418. (34) Cao, Y. Private communication, 2002. (35) Reyes, F.; Luyben, W. L. Steady-state and dynamic effects of design alternatives in heat-exchanger/furnace/reactor processes. Ind. Eng. Chem. Res. 2000, 39, 3335-3346. (36) Chemical Engineering Magazine Online. 2002; www.che.com. (37) Verein Deutscher Ingenieure. VDI-Wa¨ rmeatlas: Berechnungsbla¨ tter fu¨ r den Wa¨ rmeu¨ bergang, 6th ed.; VDI-Verlag: Dordrecht, The Netherlands, 1991.

Received for review May 18, 2004 Revised manuscript received October 16, 2004 Accepted October 19, 2004 IE0495776

Steady-State Process Optimization with Guaranteed Robust Stability

Recommend Documents