Global Optimization of Highly Nonlinear Dynamic ... - ACS Publications

Ind. Eng. Chem. Res. 2008, 47, 2643-2655

2643

Global Optimization of Highly Nonlinear Dynamic Systems Antonio Flores-Tlacuahuac* and Sebastian Terrazas Moreno Departamento de Ingenierı´a y Ciencias Quı´micas, UniVersidad Iberoamericana, Prolongacio´ n Paseo de la Reforma 880, Me´ xico D.F. 01210, Me´ xico

Lorenz T. Biegler Department of Chemical Engineering, Carnegie-Mellon UniVersity, 5000 Forbes AVenue, Pittsburgh, PennsylVania 15213

We address the determination of globally optimal solutions for a class of dynamic systems with highly nonlinear behavior. Specifically, to determine optimal transition trajectories between steady-state operating points, we consider the impact of output multiplicities in obtaining multiple optima around regions featuring severe nonlinearities. For dynamic optimization, we found that current global optimization solvers tend to demand much more computational effort when compared to local optimization solvers. Therefore, with appropriate variable bounds, good initialization strategies, and, if necessary, multiple restart strategies, local optimization solvers are often competitive, in terms of CPU time, when faced with the decision of finding global optimal solutions. Moreover, while the presence of output multiplicities does not seem to be strictly necessary for multiple local optima, we find that nonunique solutions may emerge for transitions around nonlinear regions. In order to quantify the presence of multiple local optima, we propose an empirical indicator that seems able to predict the potential for the emergence of unique optimal solutions. 1. Introduction Dynamic optimization tools have been used to improve the operation of continuous and batch chemical processes. In this framework, time-dependent values are computed for a set of manipulated variables that allow a given objective function to be maximized/minimized. Commonly, the objective function is a function of the controlled variables or is related to them in some way. To address dynamic optimization problems featuring smallscale systems, it has been customary to apply the Pontryagin maximum principle.1 However, this approach becomes difficult to handle for larger systems or with path constraints. To address these problems, methodologies based on sequential2 and simultaneous3 approaches have been widely used. Here, we focus on the simultaneous approach, where the set of index 1 differential algebraic equations (DAEs), comprising the mathematical model which is part of the dynamic optimization problem, is discretized and transformed into a set of algebraic equations. Hence, the dynamic optimization problem is transformed into a nonlinear program (NLP), which can be solved by standard techniques.4 Because dynamic optimization problems may lead to nonconvex NLPs that may be difficult to solve, it is likely that an optimal solution determined with a local solver suffices for many practical purposes. However, in some cases, it is important to ensure the global dynamic solution. In this work, we consider global optimization strategies for an important class of dynamic optimization problems, the determination of optimal transitions from one operating point to another. Global optimization strategies for dynamic systems have been developed recently with sequential optimization strategies. The RBB algorithm was described in refs 5 and 6 and applied in refs 7-11. In particular, in the paper by Esposito * To whom correspondence should be addressed. E-mail: [email protected]. Phone/Fax: +52(55)59504074. http:// 200.13.98.241/∼antonio.

and Floudas,8 the sequential approach was compared with a simultaneous collocation approach on several small-parameter estimation problems. They found that, while local solutions could be determined much faster with the simultaneous approach, finding global solutions with RBB was much more expensive with this approach. Global underestimators and algorithms were also developed for linear dynamic systems in ref 12 and generalized to mixed-integer dynamic optimization problems in ref 13. However, using either simultaneous or sequential strategies, optimal transition problems can lead to large NLPs with many degrees of freedom, which can exceed the capabilities of current global solvers. Also, highly nonlinear steady-state and dynamic behavior with output multiplicities, oscillatory behavior, and high parameter sensitivity may be present and operating regions could exhibit both stable and unstable open-loop steady states. These features lead to well-known difficulties for the sequential approach. As seen in the phase-plane diagram shown in Figure 1, the problem class considered in this paper has significant regions of instability at which the sequential approach can fail. As a result, we need to apply the simultaneous approach, as this is the only method guaranteed to find optimal solutions at open-loop unstable points.3 Recently, a heuristic combination of stochastic/deterministic methods for global dynamic optimization has been proposed.14 However, no convergence theoretical proof (toward global solution) is provided. Moreover, the algorithm seems to be problem oriented in the sense that it must be tuned for each example and, because a sequential method is used, it tends to demand large CPU time. One of our aims in this study is to evaluate the advantages and disadvantages of state-of-the-art global optimization solvers (BARON and OQNLP) compared to well-known local NLP solvers (CONOPT and MINOS) when faced with the dynamic optimization of nonlinear systems. This comparison was

10.1021/ie070379z CCC: $40.75 © 2008 American Chemical Society Published on Web 03/15/2008

2644

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 1. Bifurcation diagrams for the addressed cases; only plots for the controlled variable are shown: (a) consecutive reaction system, (b) biochemical reactor, (c) Hicks CSTR with output multiplicities, and (d) two CSTRs with recycle. The continuous lines denote open-loop stable steady states, while dashed lines denote open-loop unstable steady states. The A, B, C, etc. letters stand for the desired products in each case.

performed within the GAMS15 modeling system using default options for all of the solvers. The branch-and-reduce optimization navigator (BARON)16 finds the global solution of nonlinear (NLP) and mixed-integer nonlinear programs (MINLP). It is based on a spatial branchand-bound strategy with a variety of constraint propagation and duality techniques that reduce variable ranges over the course of the solution process. In GAMS, BARON uses the MINOS local solver. Also, BARON is the only solver that provides both upper and lower bounds on the objective function value for the global optimum. OQNLP17 is a multistart heuristic algorithm designed to find the global optimum of constrained NLPs. By “multistart”, we mean that the algorithm calls an NLP solver from multiple starting points, keeps track of all feasible solutions found by that solver, and reports back the best one as its final solution. The starting points are computed by a scatter-search implementation called OptQuest (see www.opttek.com). In this approach, random starting points are determined and refined so that they remain reasonably far apart from previous solutions and distribute evenly in the search space. This is accomplished through the use of distance and merit filters, which ensure that only a few of the randomly generated starting points are

Table 1. Nominal Steady States for the First Case Study: Consecutive Reaction Systema

x1 x2 x†3 x4 qc

A

B

C

0.0016 0.1387 8.5188 2.1729 2

0.5917 0.4249 1.7306 0.0013 1.7

0.8495 0.1503 0.2871 -0.6323 2.5

a The dagger (†) stands for the controlled variable, and q is the c manipulated variable.

acceptable. As a result, fewer NLP problems are solved, but with a higher probability of eventually finding the global solution. In this study, the CONOPT local solver was used in OQNLP. While it does not guarantee global solutions nor determine lower bounds, it is claimed to perform well in finding the global solution. A more general case study comparison of global methods on a large number of test problems was recently published in ref 18. In our study, a set of five case studies featuring different types of continuous stirred tank reactor models are used to demonstrate the difficulty in computing solutions using global optimization solvers. Here it should be noted that good local solvers and a carefully stated optimization formulation (i.e.,

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2645 Table 2. Summary of Results for the First Case Study: Consecutive Reaction CSTR, with ((Objective Function Value)/(CPU s))a transition

CONOPT

MINOS

OQNLP

BARON

ηij

AfB BfA AfC CfA BfC CfB

6.3702/0.04 4.4278/5.358 9.5958/2.454 11.5673/6.109 0.2657/0.17 0.4594/0.681

I I I I 0.2657/0.081 0.4594/0.04

6.3702/1079.76* 4.4278/1430.8* 9.5958/1028.56* 11.5677/1035.05* 0.2657/1302.6* 0.4594/1191.72*

N/1002.28 4.5055/689.58 N/421.07 11.7849/829.09 0.2657/1001.57 0.4594/1001.99

5.3 22.23 3.84 143.36 1.76 15.71

a Solved on a 1.6 GHz Pentium M laptop. In all the cases, 20 finite elements and 3 internal collocation points were used, leading to 741 variables and 681 constraints. I means an infeasible problem, N indicates that the solver could not find a solution, and the asterisk (*) means that the solution reported was obtained when the elapsed time exceeded 1000 s.

using a good initialization strategy, adequate variable bounds, etc.) are necessary even if they serve as components for global solvers. 2. Case Studies To compare the performance of local and global optimization solvers, we consider five case studies that feature different degrees of nonlinearity. In all cases, the dynamic optimization problem can be stated as follows:

∫Θ0 |z(t) - zˆ|2 dt

(1)

dz(t) ) F(z(t), u(t), t)z(0) ) z0 dt

(2)

min z,u

s.t.

zl E z E z u

(3)

u l E u E uu

(4)

where F is the vector of right-hand sides in the ordinary differential equation (ODE) model, z is the differential state vector with initial conditions z0, zˆ is the new desired operating state, u is the control profile vector, and Θ is the transition horizon. The superscripts l and u denote lower and upper bounds. The dynamic optimization problem is converted into an NLP by approximating state and control profiles by a family of polynomials on finite elements (0 < t1 < ... < tne ) Θ). Here, we use the following monomial basis representation for the differential profiles,

( )

ncol

z(t) ) zi-1 + hi

∑

Ωq

q)1

ti - ti-1 dz hi

dti,q

(5)

where zi-1 is the value of the differential variable at the beginning of element i, hi is the length of element i, and dz/ dti,q is the value of its first derivative in element i at the collocation point q. Also, Ωq is the polynomial of order ncol, satisfying Ωq(0) ) 0 and Ω′q(Fr) ) δq,r where q, r ) 1, ..., ncol and Fr is the location of the rth collocation point within each element. Continuity of the differential profiles is enforced by

where coefficients ui,q represent the values of the control variables, respectively, in element i at collocation point q. ψq is the Lagrange polynomial of order ncol satisfying ψq(Fr) ) δq,r for q, r ) 1, ..., ncol. Note from eq 5 that differential variables are required to be continuous throughout the time horizon, while control variables are allowed to have discontinuities at the element boundaries. Polynomial approximations 5 and 7 are used to convert the dynamic optimization problems (eqs 1-4) to an NLP. Specifically, eq 1 is replaced by the appropriate quadrature formula and eq 2 is enforced only at collocation points in each finite element. The continuity equation for the differential state profiles (eq 6) is also added. Therefore, the resulting NLP has a convex objective function and linear interpolation constraints, and the only source of nonconvexity is from the nonlinear dynamic model. The dynamic models for each case study are given next. As seen in the phase-plane diagram shown in Figure 1, all of these models have significant regions of instability at which the sequential approach can fail. As a result, we need to apply the simultaneous approach. In all cases, we verified that 20 equally spaced finite elements with three internal collocation points (ncol ) 3) yielded an adequate approximation for these dynamic problems. Results for these cases are discussed in the next section. 2.1. Case I: Consecutive Reaction System, X f Y f Z. The first example refers to a nonisothermal continuous stirred tank reactor (CSTR) where the consecutive reactions X f Y f Z take place.19 The dynamic behavior of the cooling jacket has also been included. The dimensionless dynamic model is given as follows.

dx1 ) q(x1f - x1) - x1κ1(x3)φ dt

(8)

dx2 ) q(x2f - x2) - x2φSκ2(x3) + x1φκ1(x3) dt

(9)

dx3 ) q(x3f - x3) + δ(x4 - x3) + βφ[x1κ1(x3) + dt Rx2κ2(x3)S] (10) dx4 ) δ1(qc(x4f - x4) + δδ2(x3 - x4)) dt

(11)

(6)

κ1(x3) ) e(x3)/(1+(x3/γ))

(12)

Radau collocation points are used because they allow constraints to be set easily at the end of each element. The control profiles are approximated using a similar representation of the form

κ2(x3) ) e(ψx3)/(1+(x3/γ))

(13)

ncol

zi ) zi-1 + hi

∑ Ωq

q)1

ncol

u(t) )

( )

ti - ti-1 dz hi

( )

∑ ψq q)1

dti,q

ti - ti-1 hi

ui,q

(7)

where x1 is the dimensionless concentration of reactant X, x2 is the dimensionless concentration of reactant Y, x3 is the dimensionless reactor temperature, and x4 is the dimensionless cooling jacket temperature. Additional dimensionless parameter values are given as follows: q ) 1, x1f ) 1, x2f ) x3f ) 0, x4f

2646


Figure 2. Optimal solutions for the first case study: consecutive reaction system, X f Y f Z; transitions, A f B, A f C, B f C, and C f B.

Figure 3. Local and global solutions for the first case study: consecutive reaction system, X f Y f Z; transition, B f A.

) -1, β ) 8, φ ) 0.133, δ ) 1, R ) 1, S ) 0.01, ψ ) 1, δ1 ) 10, δ2 ) 1, and γ ) 1000. The consecutive reaction system features highly nonlinear behavior, as can be seen from Figure 1a, where the bifurcation diagram of this system for the x3 state is presented using the cooling water flow rate (qc) as the continuation parameter. Three steady states, A, B, and C, are defined and shown in Table 1 around the nonlinear operating region. In Table 2, a summary of both the local and global optimization calculations is shown. In all the result tables, the last column refers to an empirical indicator (ηij) that we define and use later in Section 3. From Table 2, the transitions A f B, A f C, B f C, and C f B appear to exhibit a single optimal solution, and because the solution is the same using any solver, we believe the solutions

Table 3. Nominal Steady States for the Second Case Study: Bioreactor Systema

x1 x†2 Da

A

B

C

D

0.1683 0.8926 1

0.2685 0.6343 0.56

0.2335 0.3363 0.65

0.1354 0.1581 0.8

a The dagger (†) Stands for the controlled variable, and Da is the manipulated variable.

in Figure 2 are unique. In contrast, the transitions B f A and C f A feature more than one optimal solution; they are plotted in separate graphs (see Figures 3 and 4, respectively) to better appreciate different optimal solutions. 2.2. Case II: Biochemical Reactor. The second example deals with an isothermal biochemical reactor whose bifurcation

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2647

Figure 4. Local and global solutions for the first case study: consecutive reaction system, X f Y f Z; transition, C f A. Table 4. Summary of Results for the Second Case Study: Bioreactor System, with ((Objective Function Value)/(CPU s))a transition

CONOPT

MINOS

OQNLP

BARON

ηij

AfB BfA AfC CfA AfD DfA BfC CfB BfD DfB CfD DfC

0.0074/0.625 0.0029/0.191 0.0204/0.563 0.0173/0.125 0.0505/0.488 0.0612/0.414 0.0097/0.375 0.0124/0.488 0.0423/0.688 0.0701/0.375 0.0156/0.34 0.0206/0.219

0.0074/0.531 0.0029/0.527 0.0204/0.406 0.0173/0.438 0.0505/0.52 0.0934/0.547 0.0097/0.293 0.0124/0.238 0.0423/0.238 0.0701/0.582 0.0156/0.203 0.0206/0.531

0.0074/251.422 0.0029/193.219 0.0204/152.984 0.0173/162.031 0.0505/150.062 0.0612/170.656 0.0097/252.344 0.0124/211.906 0.0423/280.234 0.0701/176.703 0.0156/179.781 0.0206/324.078

0.0074/1002.48* 0.0029/1000.39* 0.0204/1000.81* 0.0173/1000.45* 0.0505/1001.12* 0.0612/1001.11* 0.0097/1001.36* 0.0124/1001.27* 0.0423/1000.91* 0.0701/1000.3* 0.0156/1000.16* 0.0206/1000.48*

0.658 0.52 1.78 3.07 4.11 18.58 2.92 6.40 1.75 10.04 2.3 6.01

a Solved on a 1.6 GHz Pentium M laptop. In all the cases, 20 finite elements and 3 internal collocation points were used, leading to 341 variables and 281 constraints. The asterisk (*) means that the solution reported was obtained when the elapsed time exceeded 1000 s.

analysis was undertaken by Agrawal et al.20 The dimensionless model is given as follows,

dx1 ) -x1 + DaM(x2)x1 dt

(14)

dx2 ) -x2 + DaΣ(x2)x1 dt

(15)

where x1 is the normalized cell concentration, x2 is the substrate conversion, Da is the Damkohler number, M(x2) ) (1 - x2) e-x2/γ is the normalized specific growth rate, and Σ(x2) ) M(x2)(1 + β)/(1 + β - x2) is the normalized substrate consumption rate. Using β ) 0.1 and γ ) 0.4, we obtain the bifurcation diagram depicted in Figure 1b for the x2 state, while the nominal steady states are shown in Table 3. 2.3. Case III: CSTR with Output Multiplicities. The third example, and its corresponding embedded transitions, was selected in order to cover a highly nonlinear operating region where global solutions, according to the empirical indicator proposed in Section 3, would be expected. In order to compute transition trajectories, the CSTR model as proposed by Hicks and Ray21 was used. Because the original parameter set used by these authors did not lead to multiple

Table 5. Nominal Steady States for the Third Case Study: Hicks CSTR, with Output Multiplicitiesa y†1 y2 u

A

B

C

D

0.0944 0.7766 340

0.1367 0.7293 390

0.1926 0.6881 430

0.2632 0.6519 455

a The dagger (†) stands for the controlled variable, and u is the manipulated variable.

steady states, some of the values were modified in order to end up with a multiplicity map. In dimensionless form, the model is given by

dy1 1 - y1 ) - k10 e-N/y2y1 dt θ

(16)

dy2 yf - y2 ) + k10 e-N/y2y1 - Ru(y2 - yc) dt θ

(17)

where y1 stands for dimensionless concentration (c/cf), y2 is the dimensionless temperature (T/Jcf), the dimensionless coolant temperature (Tc/Jcf), yf is the dimensionless feed temperature (Tf/Jcf), and u is the cooling flowrate. The parameter values are as follows: θ ) 20 min, J ) 100 K‚L/mol, cf ) 7.6 mol/L,

2648


Table 6. Summary of Results for the Third Case Study: CSTR with Output Multiplicities, with ((Objective Function Value)/(CPU s)a transition

CONOPT

MINOS

OQNLP

BARON

ηij

AfB BfA AfC CfA AfD DfA BfC CfB BfD DfB CfD DfC

302.7021/0.16 292.411/0.129 2245.9108/0.168 2013.5326/0.109 8952.4359/0.148 7068.1908/0.18 593.3461/0.117 543.5424/0.098 4382.7592/0.23 3597.5496/0.191 1087.532/0.18 940.4929/0.121

304.0221/0.082 293.283/0.102 2255.8604/0.07 2019.2966/0.07 8697.2679/0.07 7097.5597/0.059 597.7491/0.098 545.9115/0.07 4410.9352/0.07 3616.2631/0.078 1095.753/0.07 946.3959/0.113

302.7019/39.807 292.4109/91.601 2245.9104/104.46 2013.5325/94.596 8952.4357/93.184 7068.1905/101.065 593.3463/35.32 543.5424/115.225 4328.7586/108.746 3597.5497/42.631 1087.5316/82.448 940.4928/37.03

366.8274/171.78 346.1666/45.61 2389.9938/115.99 2050.4504/41.62 9277.5929/63.55 7952.0554/46.85 627.6879/35.68 693.5091/81.03 4328.8433/88.73 3771.7455/50.88 1156.9256/87.64 1096.7578/73.16

3.04 2.41 3.93 2.43 5.29 2.54 3.98 3.12 5.55 3.36 6.3 4.88

a Solved on a 1.6 GHz Pentium M laptop. In all the cases, 20 finite elements and 3 internal collocation points were used, leading to 401 variables and 341 constraints.

Table 7. Nominal Steady States for the Fourth Case Study: Sequence of Two CSTRs with Recyclea product

x1

θ1

x2

θ†2

Da

A B1 B2 C1 C2 D1 D2 E1 E2 F

0.3629 0.0979 0.3566 0.0985 0.3799 0.1048 0.3907 0.8846 0.3533 0.9722

2.3480 0.4049 2.2594 0.3596 2.3774 0.3553 2.4049 5.8964 2.0872 6.4840

0.5125 0.6001 0.6008 0.7008 0.7004 0.8002 0.8001 0.9002 0.9005 0.9809

1.8795 3.8178 2.5435 4.5371 3.1421 5.2180 3.8042 2.0798 4.7090 2.2257

0.047 0.028 0.04837 0.022 0.04664 0.01937 0.04629 0.0196 0.0507 0.05

a The dagger (†) stands for the controlled variable, while Da is the manipulated variable.

R ) 1.95 × 10-4 1/L, Tf ) 300 K, k10 ) 300 1/min, Tc ) 290 K, and N ) 5. This set of parameter values leads to the multiplicity region shown in Figure 1c. Note that products A and B are manufactured around open-loop stable steady states. Product C is located at a stability transition, which is also a Hopf bifurcation point. Finally, product D is manufactured within an unstable open-loop operating region. Process operating conditions leading to each product are shown in Table 5. Local and global optimization results are shown in Table 6. 2.4. Case IV: A Sequence of Two CSTRs with Mass and Energy Recycle. The next example deals with a system of two cascaded CSTRs where mass and energy is recycled to achieve larger reactant conversion at lower reactor temperatures.22 The dimensionless dynamic mathematical model is given as follows:

(

dx1 θ1 ) (1 - λ)x2 - x1 + Da1(1 - x1) exp dt 1 + θ1/γ

)

(18)

dθ1 ) (1 - λ)θ2 - θ1 + Da1B(1 - x1) × dt θ1 exp - β1(θ1 - θc1) (19) 1 + θ1/γ

(

)

(

θ2 dx2 ) x1 - x2 + Da2(1 - x2) exp dt 1 + θ2/γ

)

(20)

dθ2 ) θ1 - θ2 + Da2B(1 - x2) × dt θ2 - β2(θ2 - θc2) (21) exp 1 + θ2/γ

(

)

where x and θ are the dimensionless concentration and temperature, respectively; the subindex 1 refers to the first reactor, whereas the subindex 2 refers to the second reactor. Da

stands for the Dahmkoler number, θc is the dimensionless cooling temperature, λ is the recycle ratio, B is the dimension-less heat of reaction, γ is the dimensionless activation energy, and β is the dimensionless heat transfer coefficient. The bifurcation diagram of this two-reactor system is displayed in Figure 1d. Additional parameter values are given as follows: λ ) 0.9, γ ) 1000, B ) 22, θc1 ) θc2 ) 0, β1 ) 2, β2 ) 2. Using the bifurcation diagram for product design, we carry out a series of product transitions around the highly nonlinear operating region as displayed in Figure 1d, with the scaled temperature of the second reactor as the controlled variable. The operating conditions for manufacturing the set of products are shown in Table 7. The statistics for local and global optimal solutions are shown in Table 8. 2.5. Case V: A High-Impact Polystyrene CSTR. In this example, we use a polymerization CSTR displaying highly nonlinear behavior. We have previously used this example23 to demonstrate a simultaneous dynamic optimization approach for grade transitions involving open-loop unstable steady states. The process consists of a CSTR where the high-impact styrene polymerization takes place. Moreover, industrial validation of the model has been assessed.24 The problem is interesting since it represents a typical industrial-scale situation, and although we have simplified the kinetic mechanism, the model features strong nonlinearities to make the computation of global dynamic transitions a challenging task. Because of wide variations in the numerical magnitude of the states and manipulated variables, we have taken the HIPS dynamic model as described in ref 23 and applied logarithmic scaling, i.e., xj ) ln(x), where xj stands for the scaled variable. Using this type of scaling, all the decision variables have similar magnitudes. The resulting dimensionless dynamic model is given as follows.

dC h i eFh i+Ch oi -Ch i - eFh T h ) - Ad e-Ed/Re dt V

(22)

dC h m eFh +Ch om-Ch m - eFh o o T h ) - Ap e-Ep/Re (eµj r + eµj b) dt V

(23)

dT h eFh +Th o-Th - eFh ) dt V o o T h ∆HAp e-Ep/Re eCh m(eµj r + eµj b) FsCps eTh

-

UA(1 - eTh j-Th ) (24) FsCpsV

dT h j eFh j+Th oj -Th j - eFh UA(eTh j-Th -1) ) + dt V FjCpjVj

(25)


dC h b eFh +Ch ob-Ch b - eFh o T h T h ) - (Ai2 e-Ei2/Re eCh r + Afs e-Efs/Re eµh r + dt V o T h Afb e-Efb/Re eµj b) (26) dC hr T h T h ) 2f*Ad e-Ed/Re eCh i-Ch r - (Ai1 e-Ei1/Re eCh m + dt T h Ai2 e-Ei2/Re eCh b) (27) dC h br T h ) eCh b-Ch br[Ai2 e-Ei2/Re eCh r + dt o o T h T h Afb e-Efb/Re (eµj r + eµj b)] - [Ai3 e-Ei3/Re eCh m + T h

At e-Et/Re (eµj r + eµj b + eCh br)] (28) o

o

dµ j or o o T h T h C h ) 2Aio e-Eio/Re e3Ch m-µj r + Ai1 e-Ei1/Re eCh r+e m-µj r + dt o o o T h T h Afs e-Efs/Re eCh m-µj r (eµj r + eµj b) - [Ap e-Ep/Re eCh m + T h

T h

At e-Et/Re (eµj r + eµj b +eCh br) + Afs e-Efs/Re eCh m + o

o

T h

T h

Afb e-Efb/Re eCh b] + Ap e-Ep/Re eCh m (29) dµ j ob o T h T h ) Ai3 e-Ei3/Re eCh br+Ch m-µj b - [Ap e-Ep/Re eCh m + dt o o T h T h At e-Et/Re (eµj r + eµj b + eCh br) + Afs e-Efs/Re eCh m + T h

T h

Afb e-Efb/Re eCh b] + Ap e-Ep/Re eCh m (30) In terms of the logarithmic scaling, C h i is the initiator concentration, C h m is the monomer concentration, C h b is the butadiene concentration, C h r is the radicals concentration, C h br is the branched radicals concentration, T h is the reactor temperature, T h j is the jacket temperature, µ j or is the radical h i is zeroth moment, µ j ob is the butadiene zeroth moment, F the initiator feed stream volumetric flow rate, and F h is the monomer feed stream volumetric flow rate. Moreover, Table 9 contains parameter values of the HIPS CSTR in addition to its

Figure 5. Optimal solutions for the second case study: biochemical reactor.

physical measurement units. Optimal dynamic transitions are displayed in Figure 12. 3. Observations on Multiple Solutions In the above examples, we found that the application of global solvers was expensive even on these small examples (as measured in terms of the number of states and manipulated variables and not in terms of the size of the discretized system). We note again that the significant instability of these examples precluded the use of a sequential approach. For larger dynamic examples, it is often prohibitive. On the other hand, we observed that the appearance of local solutions appears to follow a pattern, which may be useful as a predictive tool for assessing whether global optimization is necessary for transitions. It is interesting that some transition cases do not exhibit multiple optimal solutions. Although presently the relationship between nonlinear behavior and multiple optimal solutions is not fully clear, the emergence of nonlinear behavior might be a prerequisite for a system to exhibit multiple optima; if a given system displays multiple optimal dynamic trajectories, these trajectories must be embedded around operating regions that exhibit highly nonlinear behavior. However, if we operate the system around these regions, this is not a guarantee that such a system would necessarily feature multiple optimal solutions. For transition optimization, there is a second requirement that has to do with the magnitude of the relative change of the controlled variable (∆x) related to the relative change in the manipulated variable (∆u). Here we propose an index ηij to consider the transition from steady state i to steady state j (corresponding in eq 1 from z(0) to zˆ, respectively),

ηij )

∆x ∆u

(31)

where ∆x ) (|xi-xj|)/(|xi|) and ∆u ) (|ui-uj|)/(|ui|). The index ηij is a gain-like variable that measures the size of

2650


Figure 6. Optimal results for the third case study: biochemical reactor, (a) A f B and (b) B f A transitions.

Figure 7. Optimal results for the third case study: biochemical reactor, (a) A f C and (b) C f A transitions.

the system response for a unit change in the system input. Note that evaluation of this index requires only steadystatea priori information. Such an indicator provides a quantitative measure of the system sensitivity to changes in the manipulated variable. In the tables containing the operating conditions for manufacturing each product, we indicate the “controlled” variable (x) with respect to which the ηij index was computed. The manipulated variable (u) is also indicated in the same table. From the cases below, we observe that transitions with low sensitivities generally have unique

transitions. In other words, low values of ηij correspond to cases that do not have multiple optima and do not require global solvers. We next discuss the five cases from the previous section. In assessing the presence of multiple optima, we apply two local solvers, two global solvers, and multiple starting points. Except for the BARON results, all of the reported solutions in these cases satisfy the KKT conditions, with each element of the reduced gradients being no greater than 10-6. While this is not an exhaustive test for multiple optima, it nevertheless provides


Figure 8. Optimal results for the third case study: biochemical reactor, (a) A f D and (b) D f A transitions.

Figure 9. Optimal results for the third case study: biochemical reactor, (a) B f C and (b) C f B transitions.

a reasonable guideline of which optimal solutions can be found with a local solver. 3.1. Case I: Consecutive Reaction System, X f Y f Z. For the CSTR with consecutive reactions, the transitions (A f B and A f C cases) do not exhibit multiple optima, although they do in the opposite direction (B f A and C f A cases). Moreover, there are also two transitions (B S C) that exhibit unique optimal solutions in both directions. From the information shown in Table 1, we see that the transition that features the smaller η value (B f C) corresponds to the one that does

not exhibit multiple optimal dynamic solutions. On the other hand, large ηij values seem to indicate the potential existence of multiple optima (as was indeed the case for the B f A and C f A transitions). Although it is difficult to postulate the ηij value for which multiple optima emerge, an approximate value 2 seems to be a good heuristic indicator. Of course, no proof of this statement can be offered at this time. 3.2. Case II: Biochemical Reactor. In Table 4, a summary of dynamic optimization results is shown using both local and global solvers. From this table, we notice that local optimization

2652


Table 8. Summary of Results for the Fourth Case Study: Sequence of Two CSTRs with Recycle, ((Objective Function Value)/(CPU s)a transition

CONOPT

MINOS

OQNLP

BARON

ηij

A f B1 B1 f A A f B2 B2 f A A f C1 C1 f A A f C2 C2 f A A f D1 D1 f A A f E1 E1 f A AfF FfA B1 f C1 C1 f B1 B1 f D1 D1 f B1 B1 f E1 E1 f B1 B1 f F F f B1 B2 f C2 C2 f B2 B2 f D2 D2 f B2 B2 f E2 E2 f B2 B2 f F F f B2 C1 f D1 D1 f C1 C1 f E1 E1 f C1 C1 f F F f C1 D1 f E1 E1 f D1 C2 f D2 D2 f C2 C2 f E2 E2 f C2 C2 f F F f C2 D1 f F F f D1 D2 f E2 E2 f D2 D2 f F F f D2 E1 f F F f E1 E2 f F F f E2

1.0746/1.914 0.7109/1.762 0.0173/0.570 0.0286/0.699 1.2684/1.973 0.3886/1.145 0.1246/1.688 0.0352/1.172 0.4168/0.922 0.4848/1.145 0.0277/1.818 0.2419/1.066 1.5554/5.516 0.3271/0.828 0.0004/0.313 0.0005/0.262 0.0015/0.422 0.0019/0.516 4.2672/1.055 3.4098/0.609 0.822/2.078 I 0.0589/0.328 0.0089/0.41 0.1175/1.719 0.1138/0.781 0.2377/1.66 0.0936/1.672 I 0.9443/0.379 0.0004/0.301 0.0004/0.258 4.2494/1.75 2.315/4.141 0.7979/1.402 3.4788/2.152 4.2047/1.598 3.5406/5.449 0.1123/1.02 0.0355/0.441 0.1759/0.699 0.0458/1.047 0.0472/1.082 1.9203/0.836 0.2965/1.422 4.1353/2.75 0.0215/0.891 0.0524/2.215 1.1521/0.738 1.7898/1.871 0.0017/0.367 0.0013/0.25 0.1427/0.563 3.9152/1.961

0.9532/1.270 2.4493/1.160 0.0173/0.891 0.0166/0.520 I 2.336/1.191 0.1246/2.531 I I 2.2143/1.625 2.4256/4 0.2419/3.738 I 0.3271/2.406 0.0004/0.793 0.0005/0.633 0.0015/3.156 0.0019/0.922 I I I 4.8395/24.484 0.0589/0.582 0.1213/9.203 0.1175/17.547 I I 0.0936/6.973 I 0.9085/13.902 0.0004/0.754 0.0004/0.613 I I I I I 7.3585/9.219 0.1123/1.152 0.023/1.254 I 0.0458/3.914 I 1.6294/3.492 I 7.2594/4.824 0.0215/1.609 0.0524/2.027 I 2.8175/3.594 0.0017/0.867 0.0013/0.813 I 4.7854/20.273

0.1612/875.118 0.3507/1179.486 0.0173/604.649 0.0334/574.826 0.1264/978.877 0.3881/1466.549 0.1246/1057.313 0.0352/617.053 0.1032/821.811 0.4307/1109.775 0.0277/1012.063 0.2419/1034.687 0.0421/1428.594 0.3271/1237.422 0.0004/557.071 0.0005/784.758 0.0015/832.578 0.0019/1110.391 0.1425/1225.171 1.6865/957.817 0.1271/2246.234 3.7549/1912.296 0.038/491.777 0.0089/624.377 0.1175/610.167 0.1138/552.514 0.3604/734.756 0.0936/661.16 0.0385/1159.357 0.793/802.564 0.0004/555.248 0.0004/571.041 0.2423/1742.458 2.0147/6144.262 0.2027/1239.612 5.2591/1283.171 0.3552/1816.968 3.0585/1476 0.1123/393.175 0.023/474.932 0.4761/971.036 0.0458/574.165 0.0472/950.546 1.6294/674.539 0.2965/2604.437 7.1293/1222.593 0.0215/1052.39 0.0524/972.078 0.075/958.553 1.5033/921.354 0.0017/655.512 0.0013/582.297 0.1427/925.01 4.7202/982.763

0.1612*/10003.2 2.3507*/10001.88 0.0173*/10000 0.0166*/10001.37 0.308*/10003.38 0.4447*/10000.12 0.1246*/10004.98 0.0352*/10000.47 0.4168*/10002.12 0.4306*/10000.49 0.0277/202.45 0.2419*/10004.12 0.0421/1189.41 0.3271/10004.8 0.0004/407.72 0.0005/1342.62 0.0015/2965.16 0.0019/3288.88 0.1425/3255.49 2.7231*/10003.7 0.1271/1956 2.045*/10000.34 0.0589*/10000.46 0.0089*/10000.37 0.1175*/10001.74 0.0649*/10002.21 0.1964*/10002.46 0.0936/10002.8 0.0385/652.96 0.8984*/10002.82 0.0004/1656.62 0.0004/1493.45 0.2423/3272.12 2.315/10005.03 0.2027/1521.95 4.1214*/10001.22 0.3553/3496.84 3.0585*/10006.2 0.1123*/10000.81 0.023*/10000.18 0.1019*/10001.08 0.0458*/10002.83 0.0472/1336.35 1.1077*/10000.72 0.2965/2671.06 3.1976*/10003.88 0.0215*/10005.08 0.0524*/10003.89 0.075/1166.97 2.3522*/10000.85 0.0017/2285.31 0.0013/583.06 0.1427/385.6 4.2047*/10000.08

2.55 0.75 12.12 9.22 2.6583 0.50 87.7039 52.06 3.0215 0.45 0.1828 0.07 2.8858 2.59 0.8792 0.45 1.1899 0.60 1.5175 1.95 0.5308 0.63 6.5801 5.14 11.5264 7.38 17.6745 10.01 3.7078 4.38 1.2554 1.22 4.9647 8.94 0.4003 1.71 50.6499 128.58 28.0798 23.02 5.7287 4.16 4.0484 6.13 0.3626 2.19 2.4965 2.21 5.1772 9.56 0.0452 0.11 38.1953 79.70

a Solved on a 1.6 GHz Pentium M laptop. In all the cases, 20 finite elements and 3 internal collocation points were used, leading to 621 variables and 561 constraints. The asterisk (*) means that the solution reported was obtained when the elapsed time exceeded 10 000 s, while I means infeasible problem.

Table 9. Nominal Parameter and Scaling Values for the Fifth Case Study: HIPS CSTR parameter

value

F Cfi Fi Cfm V f* Fe Fw A R Ei0 Ed Ei1 Ei2 Ei3 Ep Et Efs Efb

1.1411813934 0.9814814815 0.0015 8.6310347459 9450 0.57 0.915 1 19.5 1.9858 27 340 29 508 7 067 7 067 7 067 7 067 843 14 400 18 000

units L/s mol/L L/s mol/L L kg/L kg/L m2 cal/(mol-K) cal/mol cal/mol cal/mol cal/mol cal/mol cal/mol cal/mol cal/mol cal/mol

parameter

value

units

Tcf Tf Fcw Cfb Vc ∆h Cpe Cpw U Ai0 Ad Ai1 Ai2 Ai3 Ap At Afs Afb

294 333 1 1.0547874055 2000 -69919.56 1647.265 4045.7048 80 1.1 × 105 9.1 × 1013 1 × 107 2 × 106 1 × 107 1 × 107 1.7 × 109 6.6 × 107 2.3 × 109

K K L/s mol/L L J/mol J/(kg‚K) J/(kg‚K) J/(s‚K‚m2) L2/(mol2‚s) 1/s L/(mol‚s) L/(mol‚s) L/(mol‚s) L/(mol‚s) L/(mol‚s) L/(mol‚s) L/(mol‚s)

Table 10. Nominal Steady States for the Fifth Case Study: HIPS CSTRa product

Ci (mol/L)

Cm (mol/L)

T†r (K)

Tj (K)

Qw (L/s)

N2 A1 A2 A4 A5

10-5

6.0722 6.6231 7.7705 8.5994 0.55115

389.31 381.39 360.5 320.76 567.67

320.52 332.05 346.8 305.65 413.16

1 0.5 0.1 0.5 0.5

6× 1.3 × 10-4 6.6 × 10-4 1.28 × 10-3 0

a The dagger (†) stands for the controlled variable, while Q is the w manipulated variable.

solvers were always able to find the same global dynamic solution as reported by the global ones. Moreover, the local solvers’ CPU time is at least 3 orders of magnitude less than the corresponding CPU time of the global solvers. Because both local and global solvers lead to the same optimal solution (except for a small difference in the D f A transition), no distinction is made among them in Figure 5, which also displays optimal

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2653 Table 11. Summary of Results for Fifth Case Study: HIPS CSTR, with ((Objective Function Value)/(CPU s))a transition

CONOPT

MINOS

OQNLP

BARON

ηij

N2 f A1 A1 f N2 N2 f A2 A2 f N2 N2 f A4 A4 f N2 N2 f A5

29.1958/3.477 42.8528/2.789 169.1607/7.281 440.6729/16.273 12390.1334/16.586 15016.0385/52.734 18083.2048/537.867

I I I I I I I

29.1958/869.74 42.8525/951.237 169.1606/3380.16 440.6729/2261.762 12389.6439/2000.636 15016.0385/2239.97 18083.2048/3748.279

I I I I I I I

0.0407 0.0207 0.0822 0.0088 0.3522 0.2137 0.9163

a Solved on a 1.6 GHz Pentium M laptop. In all the cases, 20 finite elements and 3 internal collocation points were used, leading to 341 variables and 281 constraints. The asterisk (*) means that the solution reported was obtained when the elapsed time exceeded 10 000 s, while the dagger () means function evaluation error limit exceeded. I means infeasible problem.

Figure 10. Optimal results for the third case study: biochemical reactor, (a) B f D and (b) D f B transitions.

dynamic transitions between nonlinear operating regions. Even though MINOS detects a local solution for the D f A transition, we do not include it because it is similar to the ones found by the other solvers. The global solution was correctly determined by both BARON and OQNLP. Interestingly, the CONOPT local solver is also able to detect the global solution but at a quite modest computational cost when compared against the BARON CPU time. Regarding the prediction of global solutions, as seen from Table 4, the D f A transition is the only one that exhibits multiple solutions. Therefore, the D f A transition has a high η value and may feature a multiple optimal dynamic trajectory, as is indeed the case. High η values also occur for the C f B, D f B, and D f C transitions; these may also possess multiple solutions. However, as can be seen in this table, all of the optima appear to be the same. 3.3. Case III: CSTR with Output Multiplicities. From Table 6, we see that the CONOPT local solver is always able to find an optimal solution as good as OQNLP and better than the best upper bound found by the BARON multiple optimal solver. We also note that a different local solver in BARON should be able to produce better upper bounds. In addition, the CPU times for MINOS and CONOPT were about 2-3 orders of magnitude lower for the corresponding OQNLP and BARON solvers.

As can be noticed from the results shown in Table 6, the ηij index correctly predicts the presence of multiple solutions for all of these product transitions. This case study shows when multiple optimal solutions are likely to occur when product transitions are embedded within highly nonlinear regions. All the optimal dynamic trajectories are displayed in Figures 6-11. 3.4. Case IV: Sequence of Two CSTRs with Mass and Energy Recycle. Because of the highly nonlinear behavior embedded around the operating region (see Figure 1d), we perform a large number of product transitions. This allows more general conclusions to be made on the difficulty of finding the global solution and the effectiveness of the heuristic index. As shown in Table 8, for this case study, a total of 54 product transitions were computed. In most of the transitions, the ηij index correctly predicted the presence of unique solutions. However, there were a small set of product transitions (A f E1, B1 f E1, B1 f F, C1 f F, and D1 f F) for which the empirical index failed to detect unique solutions. In this case, we stress that local solutions were not far from the global ones; these were computed using a modest fraction of global CPU time. Because of the large number of product transitions, no dynamic transition trajectory plots are presented. 3.5. Case V: A High-Impact Polystyrene CSTR. This case study was chosen to stress that systems operating around nonlinear regions do not necessarily exhibit multiple optima.

2654


Figure 11. Optimal results for the third case study: biochemical reactor, (a) C f D and (b) D f C transitions.

Figure 12. Optimal dynamic transition trajectories for the fourth case study: HIPS reactor.

As can be seen from Figure 2a in ref 23 and product definition in Table 10, the addressed transitions feature low-sensitivity operating regions. Therefore, multiple optima should not necessarily be expected even for the largest N2 f A5 product transition, which still has a small ηij value. This situation is

confirmed by the ηij values in Table 11. The empirical indicator correctly predicts the lack of multiple optima. Hence, highly sensitive nonlinear regions seem to be a requisite to detect multiple solutions. Dynamic optimal transitions are displayed in Figure 12.


4. Conclusions Recent work has extended global optimization algorithms to dynamic systems. In particular, these approaches have been applied to dynamic optimization problems with few degrees of freedom and embedded DAE solvers to evaluate the dynamic model. Here, we consider instead optimal control problems for highly nonlinear systems. Because of the instabilities in these problems (shown in Figure 1), only simultaneous dynamic optimization strategies can be applied reliably. As shown on five case studies, we found that current global dynamic solvers seem to demand large computational loads; therefore, further advances are needed if we expect to use them on a routine basis. For transition optimization, however, it appears that local optimal solvers, with good bounds, initialization strategies, and possibly multiple restarts, are often a very inexpensive and reasonable alternative to global solvers. While the relationship between nonlinear behavior and multiple optima is not fully clear, we observe that steady-state transitions around operating regions with low sensitivity tend to display unique optimal solutions. Moreover, we believe that the intended steady-state transition should feature a large enough transition (measured by an empirical index ηij) in order to exhibit multiple optima. The case studies discussed in this paper support this statement, although we recognize the difficulty of stating general conclusions from a small set of examples. Although we did not consider cases where multiple optima occur far from nonlinear operating regions (because our aim was to explore the link between nonlinear behavior and multiple optima), it is clear that such a possibility might happen in systems exhibiting fast dynamic behavior (without giving rise to multiplicity regions). On the other hand, based on computations in this work, multiple optima would probably emerge for large ηij values. The reason why multiple optima are more likely around nonlinear operating regions is because a large ηij indicator signifies large changes in output values even for small input changes. In summary, even though mathematical models of our case studies lead to a nonconvex optimization problems, we do not necessarily find multiple optima. These appear to also depend upon the values of the operating parameters and not just on the nonconvex structure itself. Literature Cited (1) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Taylor and Francis: London, 1981. (2) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994, 33 (9), 2111-2122. (3) Kameswaran, S.; Biegler, L. T. Simultaneous Dynamic Optimization Strategies: Recent Advances and Challenges. Comput. Chem. Eng. 2006, 30 (10-12), 1560-1575. (4) Cervantes, A.; Biegler, L. T. Large-Scale DAE Optimization using a simultaneous NLP formulation. AIChE J. 1998, 44, 1038-1041.

(5) Adjiman, C. S.; Dallwig, S.; Floudas, C. A.; Neumaier, A. A Global Optimization Method, RBB, for General Twice-Differentiable Constrained NLPs. I. Theoretical Advances. Comput. Chem. Eng. 1998, 22, 1137-1158. (6) Adjiman, C. S.; Androulakis, I. P.; Floudas, C. A. A Global Optimization Method, RBB, for General Twice-Differentiable Constrained NLPs-II. Implementation and Computational Results. Comput. Chem. Eng. 1998, 22, 1159-1179. (7) Esposito, W. R.; Floudas, C. A. Global Optimization in Parameter Estimation of Nonlinear Algebraic Models via the Error-Invariables Approach. Ind. Eng. Chem. Res. 1998, 37, 1841-1858. (8) Esposito, W. R.; Floudas, C. A. Global Optimization for the Parameter Estimation of Differential-Algebraic Systems. Ind. Eng. Chem. Res. 2000, 39 (5), 1291-1310. (9) Esposito, W. R.; Floudas, C. A. Deterministic Global Optimization in Nonlinear Optimal Control Problems. J. Global Optimization 2000, 17, 97-126. (10) Esposito, W. R.; Floudas, C. A. Deterministic Global Optimization in Isothermal Reactor Network Synthesis. J. Global Optimization 2002, 22, 59-95. (11) Papamichail, C. S.; Adjiman, I. Global Optimization of Dynamic Systems. Comput. Chem. Eng. 2004, 28 (3), 403-415. (12) Singer, A. B.; Barton, P. I. Global Optimization of Dynamic Systems. J. Optimization Theory Appl. 2004, 121 (3), 613. (13) Chachuat, B.; Singer, A. B.; Barton, P. I. Global Mixed-Integer Dynamic Optimization. AIChE J. 2005, 51 (8), 2235-2253. (14) Balsa-Canto, E.; Vassiliadis, V. S.; Banga, J. R. Dynamic Optimization of Single- and Multi-Stage Systems using a Hybrid StochasticDeterministic Method. Ind. Eng. Chem. Res. 2005, 44, 1514-1523. (15) Brooke, A.; Kendrick, D.; Meeraus; Raman, R. A. GAMS: A User’s Guide; GAMS Development Corporation: Washington, DC, 1998. (16) Tawarmalani, M.; Sahinidis, N. V. ConVexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming: Theory, Algorithms, Software, and Applications; Kluwer Academic Publishers: Norwell, MA, 2002. (17) Lasdon, L.; Plummer, J.; Ugray, Z.; Bussieck, M. Working paper. MSIS MSIS Department, McCombs College of Business, May 2004. Math. Program. 2004, submitted for publication. (18) Neumaier, A.; Shcherbina, O.; Huyer, W.; Vinko, T. A Comparison of Complete Global Optimization Solvers. Math. Program. 2005, 103, 353356. (19) Gamboa-Torres, A. E.; Flores-Tlacuahuac, A. Effect of Process Modeling on the Nonlinear Behaviour of a CSTR Reacions A f B f C. Chem. Eng. J. 2000, 77, 153-164. (20) Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical Investigations of Dynamic Behavior of Isothermal Continuous Stirred Tank Biological Reactors. Chem. Eng. Sci. 1982, 37 (3), 453-462. (21) Hicks, G. A.; Ray, W. H. Approximation Methods for Optimal Control Synthesis. Can. J. Chem. Eng. 1971, 40, 522-529. (22) Kubicek, M.; Hofmann, H.; Hlavacek, V. Multiplicity and Stability in a Sequence of Two Nonadiabatic Nonisothermal CSTR. Chem. Eng. Sci. 1980, 35, 987-996. (23) Flores-Tlacuahuac, A. L.; Biegler, T.; Saldı´var-Guerra, E. Dynamic optimization of hips open-loop unstable polymerization reactors. Ind. Eng. Chem. Res. 2005, 44 (8), 2659-2674. (24) Verazaluce-Garcı´a, J.; Flores-Tlacuahuac, C. A.; Saldı´var-Guerra, E. Steady-State Nonlinear Bifurcation Analysis of a High-Impact Polystyrene Continuous Stirred Tank Reactor. Ind. Eng. Chem. Res. 2000, 39, 1972.

ReceiVed for reView March 13, 2007 ReVised manuscript receiVed September 25, 2007 Accepted February 4, 2008 IE070379Z

Global Optimization of Highly Nonlinear Dynamic ... - ACS Publications

Recommend Documents