A Coupled Strategy for the Solution of NLP and MINLP Optimization

Jan 30, 2009 - This problem was partially circumvented through the application of a more ...... Because SIMOP allows this type of customization, the m...
0 downloads 0 Views 590KB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 9611–9621

9611

A Coupled Strategy for the Solution of NLP and MINLP Optimization Problems: Benefits and Pitfalls Helder G. Silva† and Romualdo L. R. Salcedo*,†,‡ Laboratory for Process, EnVironmental and Energy Engineering, Faculdade de Engenharia da UniVersidade do Porto (LEPAE, FEUP) and Rua Dr. Roberto Frias s/n Porto, Departament of Chemical Engineering, UniVersity of Porto, 4200 465 Porto, Portugal

This paper presents a strategy for the solution of nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) problems, based on the coupling of the equation-oriented simulator ASCEND IV with the stochastic optimizers MSGA and MSIMPSA. Both NLP and MINLP formulations of a reactive distillation example were explored. The results show that the connection between the two software blocks was successfully established. Despite the highly nonlinear behavior of the proposed example, the results presented generally agree with those previously obtained for the same case study using other approaches, namely, SIMOP, which is a FORTRAN 77-based NLP simulation tool developed for optimization problems. Increasing the column dimension or switching the vapor-liquid equilibrium (VLE) description to a nonideal situation, which produces nonlinear models with a much larger number of equations that must be solved simultaneously, show that the results deteriorate. This is attributed to the lack of user control over the subsystems’ structures (viz, tearing variables). This problem was partially circumvented through the application of a more elaborate scheme initialization for the column variables. Introduction Nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) problems are two types of challenges within the wide family of problems related to the optimization field. They can be easily found in chemical engineering, predominantly in situations that are related to the design and control of processes. Typically, these examples may not be easy to solve and demand efficient tools able to cope with numerical difficulties such as nonlinearities, ill-conditioning, and singularitiessissues that are often associated with large-scale systems of nonlinear equations and multiple possible discrete configuration sets. Nevertheless, a suitable solution method is dependent quite often on several characteristics of the problem to be solved, and these should be taken into account to ensure convergence. Stochastic algorithms have been previously applied for the solution of NLP and MINLP problems.1-6 A common drawback of this class of algorithms (viz, the large number of objective function evaluations that generally are required) is progressively balanced by several characteristics of modern computing facilities, making these methods more attractive for process applications. The MSGA (adaptive random search) and MSIMPSA (simulated annealing based) algorithms are examples of stochastic methods that have been successfully used to manage the solution of NLP and MINLP case studies.1,3,7 These algorithms have proven to be less sensitive to ill-conditioning and nonconvexities in both the objective function and constraints, when compared to gradient-based approaches. An important requirement of stochastic optimizers is the robust evaluation of the process model, which is needed to compute the corresponding value of the objective function. These optimizers frequently rely on the capabilities of numerical solvers for constrained modeling, which are called after eliminating the degrees of freedom in the problem. Failures in * To whom correspondence should be addressed. Tel.: +351 22 508 1644. Fax: +351 22 508 1632. E-mail address: [email protected]. † LEPAE, FEUP. ‡ Rua Dr. Roberto Frias s/n 4200.

obtaining numerical solutions for given sets of decision variables can significantly impair the quality of the results, because, from the perspective of the optimizer, the corresponding points become indistinguishable from infeasible regions. Suitable problem representations using process simulators are strongly desired, because of their convergence capabilities. In a stochasticbased optimization framework, the coupling of both optimization and simulation tools therefore becomes a vital objective for a wider application of this methodology. The objective of this paper is to present a coupled strategy that combines the stochastic optimizers MSGA and MSIMPSA with the ASCEND IV8-10 equation-oriented simulation environment. To solve the NLP and MINLP formulations of the example problem, a reactive-distillation case study that was proposed in the literature11 has been studied. Results obtained by other simulation approaches are also presented for comparison purposes, with particular relevance to the NLP conclusions obtained by the application of the SIMOP6 methodology. The ASCEND Simulation Environment The acronym ASCEND stands for Advanced System for Computations in ENgineering Design. This project is set forth as follows: “ASCEND is both a large scale object oriented mathematical modelling environment and a strongly typed mathematical language”.12 The ASCEND environment is an equation-oriented simulator, where model writing and manipulation tasks are based on an object-oriented language. The application has also a user-friendly graphical interface, to manage the creation, manipulation, and navigation within the models. Furthermore, because of its open-source nature, it is also possible to extract virtually all of the simulation’s data from the problems being solved, which further simplifies their analysis. After the usual learning stage, it is quite straightforward to code models under the ASCEND environment. A real advantage is that the use of this tool implies compliance with good programming practices. That is because ASCEND’s simulation

10.1021/ie801613e CCC: $40.75  2009 American Chemical Society Published on Web 01/30/2009

9612

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

models must accomplish physical dimensional consistency for all included equations. The construction of dimensionally consistent models is quite useful for debugging and reusing purposes, especially in fields such as chemical engineering, where many variable types can coexist in the same problem. In addition to the information presented above, readers are referred elsewhere,8-10,12 where a deeper and detailed information about the ASCEND IV environment is presented. The Coupled Strategy for Simulation and Optimization Within the class of optimizers used in this work, the preferred mode of operation is the following: (1) An analysis of the number of degrees of freedom is first performed on the problem. (2) Based on the previous analysis, the user postulates a set of decision variables. (3) The optimizer will produce an initial point that reduces the process model to a square incidence matrix configuration, viz, a state where the number of unknown variables matches the number of equations. (4) The optimizer requires the numerical solution of the previous square system and computes the corresponding value of the objective function. (5) Based on the value of the objective calculated in the past model evaluations, and before convergence on the optimizer side is obtained, the algorithm updates the solution vector and goes back to step 3 above. Using this procedure, the optimization problem that is considered has only inequality constraints, which facilitates the implementation of alternative adaptation strategies. Also, to determine the feasible points in highly constrained problems, dynamic penalty schemes for constraint handling can also be activated.3 As a simulator, the ASCEND Solver QRSlv13 is prepared to receive square problem representations generated in step 3 of the previous approach. The ASCEND’s solver uses its own partitioning and precedence algorithm to divide the system of equations into several blocks. Subsequently, the simulator goes through the blocks sequentially and, depending on the block’s nature, uses the following strategy to solve the equation(s) within each one: (4.1) First, an attempt is made for the symbolic solution of equations in 1 × 1 matrix blocks, i.e., expliciting variables from single equations; (4.2) A numerical root finding (Brent method14) is performed on the single equations whenever explicit decomposition is not appropriate; (4.3) For blocks with more than one equation, a bounded feasible-path modified Newton method is used for the simultaneous solution of the corresponding algebraic equations. After convergence, the strategy moves to the next block. Because the environment is able to perform all the system structural decomposition by itself, the user only needs to consider the solver’s numerical parameters (e.g., maximum number of iterations, relative and absolute tolerances of the convergence criteria, etc.). The MSGA and MSIMPSA algorithms are written in FORTRAN 77 programming language. The ASCEND simulator has a hybrid nature that is based in Tcl/Tk procedures used to manage its numerical routines written in the C-language. Coupling both software packages is a typical mixed-programming task where, as is quite common nowadays, scripting languages are used to bridge software components written at different levels.15 In this particular case, it seems appropriate

Figure 1. Stochastic optimizers and ASCEND: the scheme of the implemented communication strategy.

to take advantage of the possibilities offered by the embedded ASCEND Tcl/Tk interpreter. This feature allows a scriptingbased connection to exchange information between the external application (optimizer) and the simulator. The communication scheme that follows these principles is represented in Figure 1, which displays the two main blocks, as well as the data exchanged between them. In this case, the optimization is performed via successive objective function value requests, where each one is associated with a set of design variable values suggested by the optimizers. MSGA and MSIMPSA then use the returned information to perform the optimization step. The connection also provides general model support, which easily allows model switching,

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

when necessary. This feature is handled by the optprocedures.a4s script, which is an implementation of instructions devoted to the connection itself, data reading, and model compilation. Apart from the general information, it is also necessary to set a priori particular information which is problem-dependent, namely: • The problem’s model (equations, constraints, and variable attributes as bounds and nominal values); • The RUNSIM procedure, which defines the problem dependent statements; • Data for the configuration of the solution to be conducted (optimizer, model file, tolerances established for the solver, etc.). In brief, the connection has two different main aspects. First, tasks related to the creation, compilation, and solution of the systems of equations associated with the model are performed by the ASCEND simulation environment. Second, one obtains a final solution through the definition of a feasible optimization path, based on the search performed by the optimizers MSGA and MSIMPSA, because the stochastic optimizers decide by how much each variable is changed. This combined strategy was tested and validated previously,16 using simple problems, both steady-state and dynamic. In the examples considered before, which included problems with several local optima, the global optimum could always be detected. Preliminary tests on the problem presented in this paper were published and the results obtained suggest a satisfactory behavior for the coupled strategy.17 This work goes further, applying the strategy to the MINLP formulation, as well as to models with a higher number of equations, because of the introduction of nonideal vapor-liquid equilibrium (VLE) definitions or larger columns. The SIMOP approach is introduced to provide a reference comparison. This methodology has several important features that are intended primarily for improving the application of stochastic methods to problems with a large number of equality constraints.6 Its structural analysis capabilities offer an automatic procedure for the choice of design variable sets, unlike that which happens under ASCEND, where the user specifies the decision set (see step 2). This feature, available within SIMOP, allows a straightforward switch between different output sets (viz, solution procedures and, thus, convergence capabilities) of a problem. Besides the possibilities of structural decomposition and partitioning, SIMOP enables the choice of different nonlinear solvers, as well as automatic code generation under FORTRAN 77 for fast and automatic (viz, black-box) connection with the stochastic optimizers in the solution of NLP problems. At the present stage of development, SIMOP does not analyze models from a conceptual point of view. Eventual decisions about the inclusion of initialization procedures are userdependent and such routines must be included manually after generating a simulation. Nevertheless, the user always has the option to choose different decision sets and automatically generate new simulations based on the minimization of their structural solution cost. Moreover, when considering the problems’ structure, it is important to keep two ideas present: (1) for a particular output set, the user has the option to perform an evaluation of the numerical solution cost, using the information on the Jacobian condition number of the nonlinear system of equations;6 (2) within SIMOP, users have the option to decompose the irreducible subsystems in two parts with the purpose to decrease the numerical burden. The procedure, which was first proposed by Luus,30 reformulates the systems as two distinct blocks: a simultaneous block, where only tearing variables are considered as estimates for Newton’s method, and

9613

a sequential block, where the remaining variables are explicitly calculated from the tear estimates. For a more-detailed description of SIMOP, the reader is referred to ref 6. Regarding the stochastic optimizers, both algorithms can handle simulation results in two different ways: (1) as a standard behavior, a feasible path is followed where simulations are evaluated by sudden-death penalty criteria. In this case, the causes of failure by either a failing step in Newton’s method or by infeasibility are considered numerically equivalent and simply discarded; (2) alternatively, with highly constrained problems, the feasible path approach may be unsuitable for obtaining a feasible solution, and both optimizers can be switched to an infeasible path scheme. MSGA1 implements a decreasing and progressive constraints relaxation strategy, while MSIMPSA3 performs a dynamic penalization scheme. In these cases, the constraints violation measure is used as complementary information to find feasible regions within the problem domain. Simulation failures related to the initialization procedure or convergence issues of the Newton-line step continue to be handled by sudden death penalization. GAMS,18 which is a widely used simulation/optimization environment, is also used here as a reference for the NLP optimization of the proposed problem. Case Study: Reactive Distillation The example presented in this paper, shown schematically in Figure 2, is the production of ethylene glycol (EG, component 3) using ethylene oxide (EO, component 1) and water (component 2) via the reaction given by reaction 1. C2H4O + H2O f C2H6O2

(1)

The process is performed in a reactive distillation unit, to promote both the reaction and separation simultaneously. The purpose is to avoid the formation of diethylene glycol (DEG, component 4) in the system, which is an undesired product that appears as the result of the reaction between ethylene glycol and ethylene oxide (reaction 2): C2H6O2 + C2H4O f C4H10O3

(2)

In this case study, the desired production is established and the column feed rate and distribution are unknown. The goal is to minimize the annual cost of operation of the column, and this must be achieved while guaranteeing the desired EG yield (P3) of 25 kmol/h. The presented model was originally proposed in the form of a MINLP formulation,11 where the number of stages (N) is included in the design variable set, leading to 3N + 1 degrees of freedom.4,6 The nonlinear objective function (Fobj) that should be minimized is given as NC

Fobj ) c0 +

N

∑ ∑F

ik + cHQB + cWQC +

ci

i)1

k)1

N

cTD1.55

∑

{

k)1

[

Hmin + 1.27

cSHD H0 +

N

∑ k)1

[

( )] Wk D2

+

Hmin + 1.27

( )]} Wk D2

0.802

(3)

In this equation, the cost is dependent on the reactants’ flow rates (Fik), the costs of steam (cH) and cooling water (cW), the reboiler and the condenser duties (QB and QC), the column diameter (D), the column height (calculated using Hmin and H0), and the liquid holdup in each tray (Wk). Those variables are weighted by the cost coefficients c0, cT, ci, and cSH, to originate the objective function value (Fobj). The problem is also subject

9614

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

Figure 2. General diagram of the reactive distillation system for the production of ethylene glycol (EG).

to inequality constraints that specify bounds for both the liquid and vapor stream rates, as well as for the liquid holdup in each tray. More-complete and detailed information regarding the entire model, its parameters, and the assumptions included in its definition can be found elsewhere.4,6,11 This problem has been extensively studied previously. Besides its nonlinear nature, it combines both reaction and separation phenomena, which can occur in the optimal solutions either simultaneously or in adjacent zones of the column, together with nonisothermal operation and steep temperature profiles.6 These features are usually associated with the presence of important nonconvexities, that originate various local solutions, similar to the ones found in the synthesis of simultaneous reaction/ separation networks. Additionally, by varying the maximum number of equilibrium stages, good control of the size of the optimization problem is possible. These combined characteristics make the present example an interesting problem for testing stochastic optimization algorithms, and an appropriate benchmark for the proposed approach. The presence of multiple optima in this problem was easily assessed with the GAMS model implementation described in more detail in the Discussion and Results section. Here, the OQNLP multistart solver19 was used, together with CONOPT,20 for the solution of the NLPs. The problem considered corresponds exactly to the NLP addressed in the beginning of the Discussion and Results section, switching off the optional distance and merit filters in OQNLP. Because the GAMS

implementation provided a fast and robust NLP model solution (typically requiring ∼1 s of CPU time per optimization run, on a standard 3.4 GHz Pentium IV Linux desktop), the time limit for the OQNLP solver was set at 1 h. Figure 3 shows the results obtained in this multistart experiment. Using default tolerances, a total of 158 and 357 distinct local solutions were identified, for columns with 10 and 20 equilibrium stages (ideal VLE). The flatness of the cost function in the left region of the plots is mainly due to the fact that the cost of the reactants (including the stoichiometric quantities) is included in the definition of the objective; this dominates other contributions, such as the energy costs, which is a characteristic that has been reported previously for the same problem.4 Even on the left side of each of the graphics, significant differences in the compositions and other column profiles could be found for each situation, before the physical dimensions of the columns changed much. These results clearly illustrate the importance of good initialization schemes and global solution methodologies, in reactive separation examples, and, thus, the possible usefulness in using stochastic optimizers. Model Definition. This paper focuses both on NLP and MINLP formulations of the reactive distillation problem presented. Because well-conditioned simulations are desired,4,6 the design variable set chosen was F1k (k ) 1,..., N), F2k (k ) 2,..., N), the fraction of the liquid stream vaporized in the reboiler (denoted as β), and Wk (k ) 1,..., N). This design variable set gives 3N degrees of freedom for the NLP formulation. In the

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

9615

(ii) Reactions rates: these are calculated in each tray by considering an equal distribution of the desired production over the reactive stages. In the other stages, where there is no liquid holdup, their value is set to zero. (iii) Composition profile: (a) The composition in the last tray is assumed to be 95% water and 5% EG (from this point, the EO molar fraction is corrected, considering the occurrence of the reaction); the tray’s composition is finally adjusted by normalization. (b) Next, the entire composition profile is computed, assuming a null fraction of EO in all the trays. (For the other components, a linear interpolation is used in the calculation and, for the reactive stages, the fraction of EO is adjusted by the same procedure used for the last tray of the column.) (iv) Temperature profile: Considering an ideal VLE description, the temperature profile is calculated using the stoichiometric equation of the vapor phase. (v) From this point forward, the alternative approaches follow distinct pathways: (a) Under SIMOP, this initialization also provides the initial guesses of the liquid compositions used by the column FORTRAN model6 (temperatures are used as normalized values for the calculation within the Newton-Raphson step); and (b) when using ASCEND, the scaling capabilities of the simulator were enabled for the model variables (as an equivalent action to variable normalization), and, because of the output set obtained after partitioning, two more assignments were added to the initialization procedure: • A flat profile over the liquid and vapor flow rates, where Lk ) L1 ( k ) 1,..., N + 1) and Vk ) V0 (k ) 0,..., N), was considered over the column and calculated as Figure 3. Local solutions identified with OQNLP: (a) 10 and (b) 20 equilibrium stages (ideal VLE).

MINLP formulation, the discrete variable N also becomes part of the set, leading to 3N + 1 degrees of freedom. It has been verified that atypical column dimensions (namely, the relationship between its diameter and height) could be avoided by constraining both the liquid and vapor streams to a maximum value of 1000 kmol/h.6,21 For that reason, this work only presents results with these constraints on the flow rates. Also, the number of stages was limited to a maximum of 20.21 Details on Model Initialization. Distillation is a classical example of high nonlinearity, because of its recycle complexity, which is a case that could demand an equation-oriented approach that sometimes requires, by itself, robustness improvements.22 This is true even for ideal behavior. Distillation models have nonlinearities that imply a special care over the generation of starting guesses. Good initial estimates are a necessary condition for a successful robust implementation of the Newton or quasiNewton methods in the solution of simultaneous nonlinear systems of equations. Therefore, model initialization is an important aspect for the present optimization study, because, quite often, this step can be critical to the success of the entire optimization study. The initialization procedure followed in this work is mainly based in a previous published scheme.4 A reference initialization procedure for the models used in ASCEND was setup using exactly the same procedure as that proposed under the SIMOP environment:6 (i) A conservative initial guess for F21 (which is not part of the decision set) is obtained considering that the water feed represents, by an excess of 50%, the desired production of EG. Considering both reactants and the stoichiometry of the reactions, a normalized first estimate for the composition of the first tray then is obtained.

N

C Pi 1 L1 ) βNC i)1 xi,1

∑

(4)

V0 ) βL1

(5)

where Pi corresponds to the production of the component i and xi,1 is the fraction of component i in the first column tray. • A first estimate for the ideal equilibrium constants in each tray was computed using the initial temperature profile previously obtained. Discussion and Results Ideal NLP Case. The ideal model written for the NLP formulation has 14N + 23 equations and the ASCEND’s partitioning algorithm generates an output set with 12 + N blocks. A major block holds a system of 13N + 12 equations, which are solved simultaneously, whereas the remaining 11 + N equations are distributed by single equation blocks, which are afterward worked out sequentially. Ten optimization runs were performed for columns with 3, 7, 10, and 15 stages. The best results found for the objective function value are presented in Table 1 and Figure 4. For comparison, the results obtained with two other methodologies (SIMOP6 and GAMS,18 with the CONOPT20,23 solver) are also presented. Despite the equivalent results that have been verified for up to 10 stages, the lowest annual cost obtained under ASCEND was achieved for a column design with seven theoretical trays, which is a case that leads to an annual operation cost of 1.534 × 107 US$/yr and with specifications presented in Table 2. When using SIMOP, 10 optimization runs were performed for each case, and it is noticeable that the solutions from the coupled strategy fall in the weak general trend observed with the number of trays previously determined.6 Nevertheless, when

9616

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

Table 1. Best Results Obtained for the NLP Case with Ideal VLE, Using Various Tools and Methodologies Fobj (× 107 US$/yr) MSGA

MSIMPSA

CONOPT

N

ASCEND

SIMOP

ASCEND

SIMOP

GAMS

3 7 10 15 20

1.655 1.534 1.549 2.625, 1.689a

1.661 1.552 1.550 1.567

1.685 1.546 1.535 8.069, 1.818a

1.662 1.550 1.548 1.574

1.653 1.529 1.518 1.513 1.511

a This result was obtained using the improved initialization approach described in the text.

Figure 4. Comparison of the best objective function values obtained using ASCEND, SIMOP, and GAMS (NLP case, ideal VLE). Table 2. Column Specifications for the Lowest Annual Cost Column Design Achieved (NLP Case with Ideal VLE)a tray number, k

F1k (mol/s)

F2k (mol/s)

Wk (m3)

Tk (K)

0 1 2 3 4 5 6 7 8

0.001 0.001 0.000 1.294 2.060 1.770 2.334

0.141 1.011 0.051 0.000 0.499 4.308 1.259

1.545 0.527 1.213 0.702 1.013 0.553 0.799

467.7 432.2 379.6 373.9 373.7 373.4 373.3

Vk (mol/s) 263.1 263.1 263.1 263.1 265.2 269.0 272.2 277.6

Lk (mol/s) 270.4 270.2 269.2 269.2 271.0 274.3 273.0 277.6

a The design corresponds to a column with seven stages and a Fobj ) 1.534 × 107 US$/yr. Parameters: diameter ) 1.68 m, height ) 10.20 m, β ) 0.973, QB ) 10.5 MW, QC ) 11.1 MW.

using ASCEND, the analysis of the results presented in Table 1 show that an evident impairment of the results occurs at 15 stages for both optimization algorithms. The greater numerical effort associated with the results deterioration can also be seen in Figure 5, which presents the reproducibility of the results, with respect to the ideal NLP cases, through the representation of the relative error distribution (fi is the best value of the final objective function for each run (see Table 1), and ft is the lowest annual cost known for the specified number of trays). Figure 5 shows that, generally, a larger relative error is associated with a larger column dimension (viz, a larger system of equations). As will be discussed later, this deterioration may be explained by the expansion of the number of equations that must be simultaneously solved when greater column dimension models are evaluated. Because the results obtained with the coupled strategy for longer columns were not as robust as those for smaller columns, changes over the initialization procedure used within the

Figure 5. Relative error distribution of the optimum results obtained using ASCEND coupled with MSGA and MSIMPSA for models with 3, 7, 10, and 15 stages (NLP case with ideal VLE).

ASCEND simulations were tested. To obtain better initial guesses for the system of equations that must be solved or, for instance, to provide a different way to deal with some variables, the following modifications were applied to the simulations tested: • The introduction of slack variables to handle both liquid and vapor flow rates in the column was tested for both NLP models (10 and 15 stages), as well as for the MINLP case. The results obtained were equivalent to the previously verified and this type of modification did not seem to be advantageous. • Because the composition profiles in the columns obtained are not linear,4 the initialization procedure was changed. Instead of a linear interpolation, the first estimate for the composition profile was computed using a sigmoid approximation. This ensures that the initial guess for the temperature profile also follows a sigmoid shape, being closer to the real profile obtained at the solution. The modified initialization procedure was tested for NLP models with 10 and 15 stages but, once again, the results were not different from those previously obtained. • To overcome possible difficulties with Newton’s method (QRSlv solver) when dealing with bounds on liquid and vapor flow rates, these were relaxed; however, no improvements were obtained over previous results. • In addition to the standard initialization procedure, two conditions were satisfied in the decision set randomly generated by the stochastic optimizers: (a) normalization of the column’s feeds within a specified stoichiometric excess interval, namely [0; 30%], and tray holdups were normalized between 2 and 4 m3. Both conditions were tested with both optimizers using a column model with 15 stages and the results obtained (also shown in Table 1) represent a significant improvement when compared to the previous ones.

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

The bests results for the NLP approach of this case study were obtained using GAMS. The model was written by careful reformulation of the original model equations, considering the aspects of scaling and nonlinearity. For the initialization, a sigmoid temperature profile for the temperatures and liquid compositions was chosen, after inspection of previously reported solutions. This was propagated throughout the column, considering the vapor compositions, reaction extents, and flow rates respectively. It was found that this scheme was needed to reach a first feasible solution and iterate to local optimality, using both the CONOPT and IPOPT solvers. However, with MINOS and SNOPT, no feasible solution was generally found, especially for larger column sizes. Because the purpose of the test was not a detailed comparison of the relative performance of existing NLP solvers, a more-sophisticated initialization scheme, such as the application of a modified Wang-Henke method,24 was not pursued. Nor was an effort made to compare the results obtained with global solvers such as BARON,25 because the version available could only choose as NLP solvers, other than CONOPT and IPOPT, between MINOS and SNOPT. Although better results were found using GAMS, it should be stated that much simpler initialization schemes were needed (and applied) to find feasible points in both the SIMOP and ASCEND methodologies, because of the use of the stochastic optimizers. The columns obtained with GAMS or the stochastic optimizers are essentially the same, with regard to its diameter, duties, and the value of β. However, because the feeds of the column have the heaviest influence over the objective function value, the differences verified and presented in Table 1 are a consequence of a distinct feed distribution in the solutions achieved. The results obtained with GAMS indicate that the optimum value, despite of the flat shape of the objective function, will be located beyond 20 stages (actually between 30 and 40 stages). The solutions reflect a clear division of the column in two different zones, namely, a reactive zone (near the top) and a distillation zone (near the bottom), with the feeds being distributed only within the reactive zone. Generally, this leads to a lower value of the total feed of the column, as well as to a smaller total liquid holdup (because the holdups are null within the distillation zone). The final result are shorter columns, which also are dependent on less feed, and, as verified by eq 3, lessexpensive columns. This column configuration agrees with the previously presented configuration,11 which also used GAMS to solve this problem. Curiously, these authors present an optimum value of 10 stages verified within a search up to a maximum of 20 trays. The performance of the stochastic optimizers can be improved by tuning their algorithms’ parameters. For example, with MSGA, with results presented in Table 1 having a maximum deviation of 3.6% when compared with GAMS, a maximum deviation of 0.6% can be achieved by performing a more exhaustive search. However, such fine-tuning is not the main purpose of this work, and the results are omitted here, for conciseness. In addition, it should be stressed that not all solutions are strictly interchangeable between the different simulation schemes. For example, the solution obtained under GAMS for seven stages is not feasible under ASCEND or SIMOP, unless the tolerance for solving the system of equations is relaxed to 10-2. This is due to different strategies used by the nonlinear solvers that are applied, but it is a minor problem.

9617

Despite being a simultaneous solution approach (ASCEND and SIMOP are systems that perform partitioning on the systems of equations), GAMS seems to be a successful tool to solve the present case study. This probably means that the local optima given by OQNLOP (Figure 3) are not difficult to surpass. When using stochastic optimizers, somewhat less-sharp results are observed, relative to the feed distribution. The feeds are distributed over the entire column, and a larger total liquid holdup is obtained. Stochastic search-based methods have more difficulties to achieve accurate solutions at the bounds (as, for example, stages with no feed) and that explains some lack of precision when the model dimension increases. MSIMPSA and MSGA provide a search on bounds feature that can help to prevent this effect.1,3 However, because of the dimension of the decision sets for larger columns, the search on bounds was not used to save CPU time. Nonideal NLP Case. The introduction of the Wilson VLE description, switching the column pressure from 1 atm to 15 atm, are conditions that have been proposed in the literature26 and, afterward, tested by other authors.4,6 When using the proposed coupled strategy, changing the VLE model leads to worse results, as presented in Table 3. The results obtained are far away of being close to the previously determined with the ideal VLE models. The NLP result obtained with SIMOP for 10 stages clearly does not present significant deterioration, even when the VLE conditions are modified. The improved initialization scheme was applied within ASCEND IV to models with 10 stages, and, similar to what happened with the ideal models, the new objective function values are significantly better that those obtained with the standard initialization procedure. MINLP Case. The coupled strategy proposed in this work was also applied to the MINLP formulation of the case study, because these optimizers are able to cope with such types of models. The introduction of the discrete variable N into the design variable set is straightforward. Because the optimizers are able to handle discrete variables, only the ASCEND model and the script files need to be changed. An indexed generic column model, with respect to the number of trays, is used. Also, some code lines were added to the management script, to cope with binary expansion calculation (see below) and the column model choice. As a default, the MSIMPSA algorithm only allows a one unit variation (positive or negative) over the discrete variables from one optimization step to the next. This behavior constraints a proper space exploration. MSGA allows greater discrete jumps but, because the evaluation of the optimizers is not a matter of main interest for this work, a common basis was established Table 3. Best Results Obtained for the NLP Case with Nonideal VLE, Using Various Tools and Methodologies Fobj (×107 US$/yr) MSGA N

ASCEND

3 7 10 15

2.357 2.133 2.501,1.865a 3.928

a

MSIMPSA SIMOP

1.579

ASCEND 3.548 3.579 4.513, 1.982a 9.066

SIMOP

1.586

This result was obtained with the improved initialization approach described in the text.

9618

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

Table 4. MINLP Case. Best Results Obtained for the MINLP Case with Both Ideal and Nonideal VLE Models by the Application of the ASCEND Coupled with MSGA and MSIMPSAa Fobj (× 107 US$/yr) VLE ideal Wilson

MSGA

MSIMPSA b

1.753 (5), 1.583 (11) 3.162 (5), 1.871b (13)

Cardoso et al.4

a

1.569 (6), 1.587 (13) 3.317 (8), 1.623a (12)

1.526 (10) 1.512 (10)

a The results obtained by Cardoso et al.4 are presented as a comparative reference. The respective column number of stages is given in parentheses. b Results obtained with the improved initialization approach.

for both tools and wider variations on the discrete variables can be achieved using the binary expansion based on eq 6:27,28 k

y)γ+

∑ (jY ) + (ε - C )Y j

k

k+1

(6)

j)1

This discrete variable calculation was implemented in all the MINLP examples as follows. The discrete variable y is calculated using an expansion based on k binary variables Yk. Because y is constrained by γ e y e ε, k is calculated knowing that Ck e ε and using eq 7. k

Ck ) γ +

∑j

(7)

Table 5. Column Specifications for the Lowest Annual Cost Column Design Achieved Using the MINLP Formulationa tray number, k

F1k (mol/s)

F2k (mol/s)

Wk (m3)

Tk (K)

0 1 2 3 4 5 6 7

0.001 0.000 1.317 1.926 1.744 2.556

1.172 3.571 1.236 0.284 0.769 0.570

1.545 0.527 1.213 0.702 1.013 0.553

454.2 393.6 374.3 373.7 373.5 373.3

Vk (mol/s) 260.2 260.2 260.2 262.3 265.9 267.0 274.8

Lk (mol/s) 267.8 266.6 263.1 263.7 267.0 269.1 274.8

a The design corresponds to a column with six stages and Fobj ) 1.569 × 107 US$/yr. Column details: diameter, 1.67 m; height, 9.34 m; β ) 0.972; QB ) 10.4 MW; QC ) 11.0 MW.

happened with ASCEND and SIMOP, there is no structural analysis and partitioning of the system of equations. The lowest annual cost column design found under ASCEND, using a MINLP model, corresponds to six theoretical trays, an annual cost of 1.569 × 107 US$/yr, and with specifications presented in Table 5. The application of the improved initialization procedure leads generally, one more time, to significant improvements on the objective function values. Moreover, the columns associated with these refined results are also better, from the viewpoint of its height/diameter ratio.29

j)1

In the reactive distillation problem that is presented, the number of trays N was considered to be constrained to 3 e N e 15, giving the binary expansion expression presented in eq 8: C1 ) 4 C2 ) 6 C3 ) 9 C4 ) 13 C5 ) 18

}

k)4 Ck ) 13

}

N ) 3 + Y1 + 2Y2 + 3Y3 + 4Y4 + 2Y5

(8)

As with the NLP formulation, both ideal and Wilson VLE descriptions were implemented with the two optimizers. A total of 30 and 10 optimization runs were respectively performed for the ideal and Wilson models. The best results found for the objective function are presented in Table 4. The MINLP model results obtained by Cardoso et al.4 are given only as a comparative reference, because this approach corresponds to a different simulation strategy, where, contrary to that which

Explaining the Behavior of ASCEND versus SIMOP The results from different approaches on the case study previously presented follow a pattern, where an increase of the column dimension means a deterioration of the associated objective functions. Figure 6 represents the output sets obtained with ASCEND and SIMOP for a column with 15 stages and ideal VLE. The figure shows that, in both cases, there exists a large subsystem of equations (identified as S1). Its relative dimension, when compared with the entire system, is presented in Table 6. The total dimension of the system is different in both approaches, because the equation for the calculation of the objective function is included in the equation set when using ASCEND. The main difference between the two simulation approaches is that, under SIMOP, the subsystem obtained under the partitioning step is divided in two parts.30,31 Using this division allows a different treatment for each part of the subsystem. Subset S1a corresponds to variables that are explicitly calculated, while the remaining ones, within subset S1b, are the tear

Figure 6. Output sets incidence matrices obtained with ASCEND and SIMOP for a column with 15 stages (ideal VLE).

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

9619

Figure 7. Output sets incidence matrices obtained with ASCEND and SIMOP for a column with 10 stages (nonideal VLE). Table 6. Comparison between the Partitioning Step Results When Using ASCEND and SIMOP (Ideal VLE)a ASCEND

SIMOP

N

neq

neqs

neq

neqs

neqt

7 10 15

121 163 233

103 142 207

120 162 232

102 141 206

22 44 68

a The parameters shown are the total number of equations in the system (neq), the total number of equations in the subsystem (neqs), and the number of tear variables of the subsystem (neqt).

Table 7. Comparison between the Partitioning Step Results When Using ASCEND and SIMOP (Nonideal VLE)a ASCEND

SIMOP

N

neq

neqs

neq

neqs

neqt

10

403

382

482

421

50

a

The parameters shown are the total number of equations in the system (neq), the total number of equations in the subsystem (neqs), and the number of tear variables of the subsystem (neqt).

variables of the subsystem variables, which are calculated simultaneously. As presented in Table 6, the division of the subsystem into two subsets leads to a very significant reduction of the number of equations that must be simultaneously solved. In addition, because the variables of the subset S1a are explicitly calculated, the computation of its initial (guess) values is unnecessary. Because SIMOP allows this type of customization, the model can be partitioned into a easier numerical problem. For this particular case, the only variables that must be initialized under SIMOP are the composition and temperature profiles, which are set as tearing variables. Switching the VLE definition from ideal to Wilson has direct consequences over the column model. It follows that new equations are introduced into the problem, changing its dimension from 14N + 23 to 38N + 23 equations (the objective function definition counts as an equation within ASCEND). Because the partitioning algorithm produces an analogous output set (see Figures 6 and 7, as well as Table 7), the new 24N equations are all part of the major block (S1), raising the number of equations to be solved simultaneously and increasing the numerical burden that must be overcome to achieve a satisfactory solution. Therefore, if increasing the column number of trays enlarges the major block set used by ASCEND, switching the VLE model description to nonideal introduces even more equations and that explains the more-severe deterioration of the

results verified with the Wilson model (when using SIMOP, the number of equations presented in Table 7 is greater than 38N + 23, because the equilibrium definition was broken apart in more-explicit terms, to fit data file structures). The deviation in the results is then a consequence that arises from the different dimension of the system of equations that must be simultaneously solved under ASCEND. Nevertheless, despite the deterioration of the results that has been verified for the coupled strategy, the proposed scheme is functional, provided a more-elaborate initialization scheme is used, and a cooperation system was established between MSGA, MSIMPSA, and ASCEND, which is an environment with powerful simulation and analysis capabilities. Moreover, the stochastic optimizers were extended beyond the application of FORTRAN 77-written simulations. ASCEND is not an application devoted for optimization as its main purpose. As previously explained, this equation-oriented environment does not provide tools to implement a choice of design variable sets or any type of feature to customize the output set that comes out of the partitioning step. On the other hand, SIMOP is able to build customizable simulations that could make the difference in this particular problem.6 SIMOP allows one to manage and change the tearing variables of a problem and establish constraints over the output sets generated, avoiding numerical issues related to ill-conditioning. Such types of procedures result in more-stable simulations, which is a necessary condition for the effective use of stochastic optimizers whenever a large set of nonlinear equality constraints occur. Conclusions An efficient and functional communication scheme between the stochastic optimizers MSGA and MSIMPSA and the ASCEND IV environment was successfully built. The stochastic optimizers were extended beyond the scope of FORTRAN 77-written simulations, which seem to be usable with alternative simulation tools. The proposed scheme had been previously applied to small-to-medium scale steady-state and dynamic problems. Although the examined systems are multimodal, the global optimum could be found for a significant number of runs.7,16,17 A proposed reactive distillation case study11 was analyzed for both nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) formulations. The results obtained are comparable to previous solutions.4,6 However, the

9620

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009

results deteriorate for nonideal models or models with more than 10 stages, where a larger number of equations must be simultaneously solved. The introduction of slack variables on flow rates and sigmoid composition and temperature profiles were implemented as alternatives to the standard initialization procedure, with no evident influence on the quality of the achieved results. However, better results are obtained when a modified initialization scheme based in the normalization of the decision set is used. The numerical burden associated with larger systems of equations is probably the reason that leads to the deterioration of the results when using the ASCEND’s proposed output set found by its partitioning algorithm. Because the customizable partitioning step of the SIMOP methodology proved to be robust, despite the NLP model’s size variability, its flexibility and robustness when building simulations suggests an extension of this tool to be applied to the solution of MINLP models. The accuracy of the solutions found by the stochastic optimizers may also be improved by these models when using binary variables to describe limiting situations such as, for example, the existence or absence of feed in each tray. This work is currently underway. Nomenclature c0 ) cost parameter (US$/mol) c1 ) ethylene oxide cost (US$/mol) c2 ) water cost (US$/mol) cH ) cost parameter related with the reboiler (US$/mol) cSH ) cost parameter related with the cost of the column (US$/ mol) cT ) cost parameter related with the cost of the column (US$/mol) cW ) cost parameter related with the reboiler (US$/mol) CK ) binary expansion coefficient D ) column diameter (m) Fik ) molar flow rate of component i feed to tray k (mol/s) Fobj ) objective function fi ) objective function value of the ith optimization run ft ) best-known objective function value Hmin ) minimum space between trays (m) H0 ) extra column height (m) Li ) liquid flow rate of tray k (mol/s) N ) number of trays of the column NC ) number of components Pi ) production rate of component i (mol/s) Pi* ) production rate of component i (value after initialization) (mol/s) QC ) condenser heat duty (W) QB ) reboiler heat duty (W) Tk ) temperature of tray k (K) Vi ) vapor flow rate of tray k (mol/s) WK ) liquid holdup of tray k (m3) * xi,1 ) molar fraction of the component i in the first tray (value after initialization) y ) discrete variable Yi ) binary variable Greek Letters β ) boil-up fraction of L1 vaporized in the reboiler γ ) discrete variable upper bound ε ) discrete variable lower bound Acronyms ASCEND ) Advanced System for Computations in ENgineering Design

DEG ) diethylene glycol EG ) ethylene glycol EO ) ethylene oxide FORTRAN ) FORmula TRANslator GAMS ) general algebraic modeling system MINLP ) mixed-integer nonlinear programming MSGA ) MINLP Salcedo Gonc¸alves Azevedo algorithm MSIMPSA ) MINLP SIMPlex Simulated Annealing algorithm NLP ) nonlinear programming SIMOP ) SIMulation for OPtimization VLE ) vapor-liquid equilibrium

Acknowledgment This work was partially supported by FCT (Portuguese Foundation for Science and Technology), under Contract No. SFRH/BD/10149/2002, employing the computational facilities of Instituto de Sistemas e Robótica-Porto (ISR-P). We are also grateful to the valuable contributions of Professor Arthur Westerberg and Dr. Ben Allan for a deeper understanding of the ASCEND IV environment. The authors further acknowledge the collaboration of Professor Nuno Oliveira (DEQ/ FCTUC, Coimbra, Portugal), who provided the GAMS code and results used in the comparison and also discussed several issues related to the implementation of the model. Literature Cited (1) Salcedo, R. L. Solving Nonconvex Nonlinear-Programming and Mixed-Integer Nonlinear-Programming Problems with Adaptive Random Search. Ind. Eng. Chem. Res. 1992, 31, 262–273. (2) Banga, J. R.; Seider, W. D. Global optimization of chemical processes using stochastic algorithms. In State of the Art in Global Optimization: Computational Methods and Applications; Floudas, C. A., Pardalos, P. M., Eds.; Nonconvex Optimization and Its Applications 7; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 563-583. (3) Cardoso, M. F.; Salcedo, R. L.; de Azevedo, S. F.; Barbosa, D. A simulated annealing approach to the solution of MINLP problems. Comput. Chem. Eng. 1997, 21, 1349–1364. (4) Cardoso, M. F.; Salcedo, R. L.; de Azevedo, S. F.; Barbosa, D. Optimization of reactive distillation processes with simulated annealing. Chem. Eng. Sci. 2000, 55, 5059–5078. (5) Ali, M. M.; Khompatraporn, C.; Zabinsky, Z. B. A numerical evaluation of several stochastic algorithms on selected continuous global optimization test problems. J. Global Optim. 2005, 31, 635–672. (6) Lima, R. M.; Salcedo, R. L.; Barbosa, D. SIMOP: Efficient reactive distillation optimization using stochastic optimizers. Chem. Eng. Sci. 2006, 61, 1718–1739. (7) Salcedo, R. L.; Lima, R. P.; Cardoso, M. F. Simulated Annealing for the Global Optimization of Chemical Processes. Proc. Indian Natl. Sci. Acad., A 2003, 69A, 359–401. (8) Locke, M. H.; Westerberg, A. W. The ASCEND-II SystemsA Flowsheeting Application of a Successive Quadratic Programming Methodology. Comput. Chem. Eng. 1983, 7, 615–630. (9) Piela, P. C.; Epperly, T. G.; Westerberg, K. M.; Westerberg, A. W. ASCEND: An Object-Oriented Computer Environment for Modeling and AnalysissThe Modeling Language. Comput. Chem. Eng. 1991, 15, 53– 72. (10) Allan, B.; Chittur, K.; Epperly, T.; Huss, R.; Pye, J.; Rico-Ramirez, V.; Shao, J.; St. Clair, J.; Thomas, M.;Westerberg, A. The ASCEND Developers Collaboration Project, 2006. Available via the Internet at https:// pse.cheme.cmu.edu/wiki/view/Ascend/AscendProject. (11) Ciric, A. R.; Gu, D. Y. Synthesis of Nonequilibrium Reactive Distillation Processes by MINLP Optimization. AIChE J. 1994, 40, 1479– 1487. (12) ASCEND Development Team, The ASCEND Project Website at Carnegie Mellon University, 2005. Available via the Internet at http:// ascend.cheme.cmu.edu/. (Accessed 2006.) (13) Westerberg, A. A modified least squares algorithm for solVing sparse n×n sets of nonlinear equations; EDRC Technical Report, 1979. (14) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 7, 1st Edition; Cambridge University Press: Oxford, U.K., 1992.

Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 (15) Ousterhout, J. K. Scripting: Higher level programming for the 21st century. IEEE Comput. Mag. 1998, 31, 23–30. (16) Silva, H. G.; Salcedo, R. L. Stochastic Algorithms and ASCEND IV: a coupled equation-oriented strategy for non-linear optimization. In Proceedings of CHEMPOR 2005s9th International Chemical Engineering Conference, Coimbra, Portugal, September 21-23, 2005; Dep. Eng. Química, FCTUC, Univ. Coimbra: Coimbra, Portugal; p 250. (17) Silva, H. G.; Salcedo, R. L. R. Modeling and Optimization of Chemical Processes: ASCEND IV and Stochastic Optimizers. In Proceedings of the 17th IASTED International Conference on Modeling and Simulation; International Association of Science and Technology for Development (IASTED), 2006; pp 493-498. (18) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMSsA User’s Guide; GAMS Development Corporation: Washington, DC, 1998. (19) Lasdon, L.; Plummer, J.; Ugray, Z.; Bussieck, M. ImproVed Filters and Randomized DriVers for Multistart Global Optimization; Technical Report, University of Texas, Austin, TX, 2004. (20) Drud, A. CONOPT; ARKI Consulting and Development A/S: Bagsvaerd, Denmark, 2003. (21) Ciric, A. R. Personal communication, 1995. (22) Morton, W. Equation-Oriented Simulation and Optimization. Proc. Indian Natl. Sci. Acad. 2003, 317–357. (23) Drud, A. S. CONOPTsA Large-Scale GRG Code. ORSA J. Comput. 1994, 6, 207–218. (24) Seader, J. D.; Henley, E. S. Separation Process Principles; John Wiley & Sons: New York, 1998.

9621

(25) Tawarmalani, M.; Sahinidis, N. V. Global optimization of mixedinteger nonlinear programs: A theoretical and computational study. Math. Prog. 2004, 99, 563–591. (26) Okasinski, M. J.; Doherty, M. F. Design method for kinetically controlled, staged reactive distillation columns. Ind. Eng. Chem. Res. 1998, 37, 2821–2834. (27) Cardoso, M. F.; Salcedo, R. L.; de Azevedo, S. F. The simplexsimulated annealing approach to continuous non-linear optimization. Comput. Chem. Eng. 1996, 20, 1065–1080. (28) Floudas,C.A.NonlinearandMixed-IntegerOptimizationsFundamentals and Applications; Oxford University Press: Cambrdige, U.K., 1995. (29) Douglas, J. Conceptual Design of Chemical Processes, 4th Edition; McGraw-Hill Higher Education: New York, 1988. (30) Luus, R. Effective solution procedure for systems of non-linear algebraic equations. Hung. J. Ind. Chem. 1999, 27, 307–310. (31) Shacham, M.; Brauner, N.; Cutlip, M. B. A web-based library for testing performance of numerical software for solving nonlinear algebraic equations. Comput. Chem. Eng. 2002, 26, 547–554.

ReceiVed for reView October 23, 2008 ReVised manuscript receiVed December 17, 2008 Accepted December 18, 2008 IE801613E