Efficient Numerical Solver for Partially Structured Differential and

Sep 15, 2009 - Given a sparse set of differential and algebraic equations (DAEs), it is ... algebraic equations or (ii) a system of ordinary different...
0 downloads 3 Views 510KB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 9979–9984

9979

Efficient Numerical Solver for Partially Structured Differential and Algebraic Equation Systems Flavio Manenti,*,† Ivan Dones,‡ Guido Buzzi-Ferraris,† and Heinz A. Preisig‡ CMIC Department “Giulio Natta”, Politecnico di Milano, Piazza Leonardo da Vinci 32, I-20133 Milano, Italy, and Chemical Engineering Department, Norwegian UniVersity of Science and Technology (NTNU), Sem Sælands Vei 4, N-7491 Trondheim, Norway

Given a sparse set of differential and algebraic equations (DAEs), it is always recommended to exploit the structure of the system’s sparsity (e.g., tridiagonal blocks matrix, band matrix, and staircase matrix, etc.), thus to use tailored numerical solvers in order to reduce the computation time. Very frequently, though, while highly structured, a couple of elements enter the description which make it difficult for the solvers to reach a solution. They are common in process control applications, where the states added to the plant description by the integral parts of the controllers introduce unstructured elements in the otherwise very structured Jacobian of the mathematical model. Such systems are characterized by a partially structured Jacobian, which inhibits the use of the solvers tailored to fit problems with fully structured matrices. In such cases, one can either use a solver with lower performance, resulting in larger computation times, or alternatively one seeks an approximation for the unstructured points. A solution to the handling of “dirty” Jacobians is presented, which is implemented in a DAE solver package available freely on the Internet. This novel DAE solver fully exploits the overall structure of the system’s sparsity, without compromising CPU computation time and precision of the results. A numerical comparison with different approaches is given by solving a DAE model representing an existing nonequilibrium distillation column. 1. Introduction Modeling of chemical processes is the foundation of any process systems engineering activity including design, control, optimization, planning, and retrofitting. With the model to meet the specifications of the application, models vary with their use and thus a multitude of models may be generated for the same plant each serving a different purpose. Limiting the discussion to the use of first-principle models, the behavior of the desired equipment is represented by either (i) a set of algebraic equations or (ii) a system of ordinary differential equations (ODEs) or differential and algebraic equations (DAEs) or (iii) partial differential and differentialalgebraic equations (PDE and PDAE), including the following: mass, energy, and momentum balances; transfer laws; possible reactions; thermodynamics and geometry, etc.; secondary variables transformations; control laws. Given the size, complexity, and coupling of these equations, solving these equations requires a good and robust numerical solver including an integrator and possibly a local root solver. It is well-known that solvers that exploit the sparse structure of its Jacobian have superior performance in terms of CPU computation time when compared with solvers that do not. Unfortunately, those solvers are not always available or applicable to the specific studied case. For the purpose of discussing the different issues, we consider the case of a generic, binary distillation column. It is well-known that first-principle distillation models (represented by a set of differential and/or algebraic equations) have a sparse structure, well-organized in tridiagonal blocks.1-5 However, this nice structure may be spoiled when control loops are introduced for the purpose of manipulating the dynamic behavior of the plant. * To whom correspondence should be addressed. Tel.: +39 02 2399 3273. Fax: +39 02 7063 8173. E-mail: [email protected]. † Politecnico di Milano. ‡ NTNU.

Control laws containing integral terms, such as proportionalintegral (PI) and proportional-integral-derivative (PID) loops, introduce new states into the representation, thus new equations to be integrated, and consequently also new elements into the Jacobian, which usually lie outside the tridiagonal band. The consequence is that solvers for sparse and well-structured DAE problems cannot be used directly for solving this type of problem. The practical consequences are that one often has to resort to a more general solver, which is less efficient in terms of CPU time, or eventually accept an approximate solution for the unstructured points in order to keep the computational time small. On this subject, this paper shows an efficient approach to precisely solve partially structured DAE systems (section 2), where an overall well-structured Jacobian is spoiled by a few unstructured elements. Some existing classes, belonging to BzzMath library,6-8 are described in section 3. Section 4 proposes a novel approach to solving partially structured DAE systems that exploit at the maximum level the sparsity and structure of the problem. An industrial case study is described in section 5. Finally, the implementation of DAE solvers, numerical results, and a comparison using the new scheme versus using an unstructured integrator are discussed in section 6. 2. Management of System Sparsity and Unstructured Elements Although advanced control and optimization methodologies have been extensively discussed in the scientific literature throughout the past 3 decades,9-18 the conventional control systems are still the most commonly used approaches in the process industry. In this framework, the integral part of each traditional control loop “corrupts” the Jacobian structure by introducing unstructured elements. These considerations motivate the relevance of this research and specifically the develop-

10.1021/ie9007908 CCC: $40.75  2009 American Chemical Society Published on Web 09/15/2009

9980

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

numerical solvers, the reader could be interested in other contributions.25,26

Figure 1. Partially structured Jacobian. X-elements constitute a tridiagonal structure; Y-elements are unstructured entries.

ment of better performing, tailored numerical solvers for partially structured systems. Actually, when the problem has not a fully structured Jacobian (see the example of Figure 1), several very well performing algorithms designed for sparse and well-structured systems cannot be used. Several approaches can be adopted to solve the posed problem: (1) The simplest way out is to use a generic DAE solver for dense systems of equations. Unfortunately, the required CPU time increases exponentially with the size of the system as the algorithm does not take advantage of the sparsity of the system’s matrix. (2) Another approach is to use a DAE solver that accepts the Jacobian existence matrix, thus an incidence matrix only. This provides the solver with the possibility to completely exploit the system sparsity, but not its overall structure. Nevertheless, it is often not possible to provide the Jacobian incidence matrix, especially when the size of the system is very large or when the incidence matrix changes as part of an iterative process. (3) An additional one solves for the unstructured part separately (e.g., quadrature of the integral), while the remaining sparse structure is efficiently calculated with a appropriate solver. While this approach results in a significant improvement in computation time, the approximate solution may not be accurate enough for the application. Particularly for control applications, these unstructured elements are related to the integral part of the control loops and an inaccurate solution may have catastrophic effects mainly if the process is highly dynamic (e.g., surge lines and radicalic systems, etc.), where also the right choice of the discretization time of the control action is critical. (4) Alternatively a novel generalized class of DAE solvers is generated with the objective to solve simultaneously and efficiently the structured and unstructured part. A single DAE solver is then able to handle both the fully structured problem and the corrupted structured problem containing unstructured elements. In the next sections we will focus on the DAE solvers contained in the BzzMath package,6-8 in particular the efficient solver for partially structured problems. 3. DAE Solvers Belonging to BzzMath Library In the past decades, many contributions were published concerning implementations of integration routines. For instance, we point to DDASSL,19,20 RADAU5,21 LSODE,22 STRIDE,23 and DASPK.24 For a more complete survey over the ODE/DAE

In this paper, we focus the attention on the DAE solvers included in the BzzMath library, which is freely downloadable at the Web site www.chem.polimi.it/homes/gbuzzi. 3.1. BzzDae Class. BzzDae is a generic DAE solver that does not need any information about structure and sparsity of the system of equations, which makes it easy to use also in all those cases where unstructured elements are involved. Its weakness may come up to the surface when the solver needs to fill up the Jacobian during the numerical solution process, since no structure and sparsity information is provided. Even more, the linear system solution, necessary for the Newton method, which is used for the convergence of the implicit method, is not quite optimal for the same reason. However, it does automatically detect the possible Jacobian sparsity: classes belonging to the BzzMath library that are designed for solving differential and differential-algebraic systems operate a check to see the sparsity degree, and, sometimes, it is also able to automatically detect the kind of structure in order to select the appropriate method to solve it efficiently. 3.2. BzzDaeSparse Class. This DAE solver requires as input the following: either the incidence matrix of the system of equations, giving so the solver information about the DAE system sparsity and saving time when it comes to populating the Jacobian matrix or information about the kind of structure of the DAE system. In the first case, the full problem (including unstructured elements) can be directly fed to the solver. However, as with BzzDae, no properties of the structure are taken into account and the solver performances are not fully optimized. On the other hand, a significant decrease in terms of CPU time is obtained in the second case, since a smart calculation method can be applied. Unfortunately, the unstructured elements cannot be directly provided to the DAE solver. In the specific case of process control, this implies that the user should evaluate the equations associated with the unstructured elements separately. Possible alternatives include the approximation of the integral part of the controllers for example using a numerical quadrature formula or implementing a discrete controller if appropriate. Both of these approaches will not be considered in this work. 4. DAE Solver for Partially Structured Systems A novel approach has recently been introduced in the BzzMath library. A new DAE solver has been created for greatly sparse systems corrupted with few unstructured elements, as they typically occur in the process control field. The implementation of this algorithm is now part of the BzzDaeSparse class. The matrix reproducing the Jacobian of Figure 1 can be rearranged as reported in Figure 2, so to obtain a matrix A that can be partitioned into four blocks: (1) A1,1 is a structured (i.e., diagonal, tridiagonal, tridiagonal blocks, and staircase, etc..) submatrix. A necessary condition is that this matrix must be well-conditioned. (2) A1,2 and A2,1 are two (usually rectangular) sparse unstructured submatrices. (3) A2,2 is a generally dense submatrix.

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

9981

The resulting four-blocks system is A1,1x1 + A1,2x2 ) b1

(1)

A2,1x1 + A2,2x2 ) b2

(2)

which can be solved readily A1,1x1 ) b1 - A1,2x2

(3)

(A2,2 - A2,1A1,1-1A1,2)x2 ) b2 - A2,1A1,1-1b1

(4)

Even though in eq 4 an inverse matrix appears, it is necessary to solve the problem without explicitly inverting the same matrix.27 To exploit the sparsity of unstructured submatrices, depending on the sparsity degree of the two matrices A1,2 and A2,1, two alternative solving strategies are defined, which are similar in the evaluation of the product z ) A1,1-1b1, but they proceed along two parallel paths: (1) The product W ) A1,1-1A1,2 is evaluated through the solution of the system A1,1W ) A1,2. Once the matrix W is known, the following product is evaluated: A2,1A1,1-1A1,2 ) A2,1W. Then, the product A2,1z is computed, and the system 4 is solved for x2. At last, the vector expression b1 - A1,2x2 is evaluated, and the system 1 is solved for x1. (2) The product V ) A2,1A1,1-1 is evaluated by solving the T T VT ) A2,1 for V. Once the matrix V is known, the equation A1,1 -1 product A2,1A1,1 A1,2 ) VA1,2 is computed. Next, the product A2,1z is calculated, and system 4 is solved for x2. At last, the vector b1 - A1,2x2 is evaluated, and system 1 is solved for x1. The first alternative is used when the matrix A2,1 is sparse and the second alternative when the matrix A1,2 is sparse. The algorithm exploits the sparsity and the structure of matrix A1,1. It is consequently factorized just once.

Figure 2. A1,1, structured (tridiagonal) block; A1,2 and A2,1, sparse unstructured blocks; A2,2, dense block.

Figure 3. Distillation column flow-sheet and control scheme.

In addition, it is important to note that the matrix A2,1A1,1-1A1,2 can be obtained column-by-column in the first case and row-by-row in the second case.2,28 The calculations can thus be done “in-memory” so the matrix must not be stored there by memory allocation. 5. Case Study for Solver Validation Nonequilibrium Distillation Column Model To test the different solvers, a distillation column model was chosen. First principles distillation models yield a tridiagonal block structure, each block representing the properties of a single distillation stage.27,29 A fully dynamic model of an isobutane and normal-butane splitter with control structure is shown in Figure 3. While the details of the model can be taken from the literature,27 here a brief summary: each stage of the column is considered as a nonequilibrium dynamic flash containment, where the liquid phase and gas phase are two uniform lumps exchanging extensive quantities. The main differences with the traditional models are as follows: dynamic mass and energy balance for both phases; gas capacity included; no chemical equilibrium assumption; no thermal equilibrium assumption; no assumption about the liquid density. The discussion of the importance of correct volume calculations can be found elsewhere.27 The dimension of the state is 2NC + 6 states per column stage, where NC stands for the number of components: NC dynamic liquid mass balances; NC dynamic gas mass balances; one dynamic energy balance for the liquid phase; one dynamic energy balance for the gas phase; one algebraic equation for liquid temperature calculation; one algebraic equation for gas temperature calculation; one algebraic equation for liquid volume calculation; one algebraic equation for gas volume calculation. Table 1 summarizes the system of equations given to the solver for the generic stage. The algebraic equations for temperatures and volumes are chosen to be solved directly by the DAE solver, since they are implicit equations and we prefer to use a common numerical solver instead of tuning tolerances for internal and external loop solvers.

9982

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

Table 1. Generic Stage Modela equations

integration states mass balance liquid phase

nLi

mass balance gas phase

nGi

energy balance liquid phase

ULi

energy balance gas phase

UGi

0 ) ULi - UL(TLi, VLi, nLi)

calculate liquid temp

TLi

0 ) UGi - UG(TGi, VGi, nGi)

calculate gas temp

TGi

mech equilib

VLi

vols constraint

VGi

n˙Li ) -nˆLi-Li+1 + nˆLi-1-Li - nˆLi-Gi - nˆLi-env + nˆenv-Li n˙Gi ) -nˆGi-Gi-1 + nˆGi+1-Gi + nˆLi-Gi

˙ L ) -H ˆ L -L + H ˆ L -L - qˆL -G - H ˆ L -B - qˆL -B - pL V˙L U i i i+1 i-1 i i i i i i i i i ˆ L -env + H ˆ env-L - qˆL -env + qˆenv-L H i i i i

˙ G ) -H ˆ G -G + H ˆ G -G + qˆL -G + H ˆ B -G - pG V˙G U i i i-1 i+1 i i i i i i i

0 ) p Li - pGi 0 ) VTOTi - (VLi + VGi)

a Notations: n, NC dimensioned mass vector; U, internal energy; T, temperature; V, volume; p, pressure; q, heat; H, enthalpy. Subscripts: i, stage position (counting from the top); L, liquid phase; G, gas phase; B, boundary between liquid and gas phases; TOT, stage property; env, external environment. Symbols: ∧, stream; · , time derivative.

The secondary variable transformations and all thermodynamic properties are calculated using an external, portable, and easy-to-plug-in thermodynamic package.29 The package accepts temperatures, volumes, and masses, returning Helmholtz’s energy and its derivatives with respect to the inputs. In this way, the liquid (or gas) pressure is calculated as a function of the liquid (or gas) temperature, volume, and mass of that phase, being the pressure the first derivative of Helmholtz’s energy A with respect to temperature for a given volume and mass, such as in ∂ALi ∂VLi TLi,nLi ∂AGi ∂VGi TGi,nGi

) -pLi

(5)

) -pGi

(6)

The pressure drop along the column is not made explicit but is hidden in the resistance that the rising gas streams nˆGi+1-Gi have to face to flow upward in the column. The law defining nˆGi+1-Gi will be a valve-like equation, thus containing the pressure drop in the form of dry and wet tray resistances. This is someway different from the work of Krishnamurthy and Taylor,2 where, in order to express the column pressure drop as a function of tray (or packing) type, for each stage one should consider a hydraulic equation and make the pressure drop of each stage an unknown variable (therefore explicit). Krishnamurthy and Taylor pointed out that many numerical solution methods (such as Newton method for the solution of the an algebraic system) require derivatives, and some of them (with respect to thermodynamic and transport properties) may be unavailable, urging the user to adopt approximations when required. Having available Løvfall’s29,30 stand-alone thermo-

Figure 4. Open-loop DAE structure of the case study: a fully structured tridiagonal block matrix.

modules makes any numerical approximation obsolete. Derivatives are computed from analytically derived expressions internal for those modules, thus providing a world of flexibility for their use in advanced algorithms. Representations of the uncontrolled plant and the plant controlled using proportional controllers only are pure tridiagonal models, as reported in Figure 4 with each block representing the above-mentioned states of a stage. The Jacobian would be highly sparse and correspondingly well-structured, so it would be smart to use structure information to ease the computation of the states.

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

9983

Figure 5. Partially structured tridiagonal block Jacobian. The points in row 941, 942 and in column 941, 942 represent the unstructured elements: precisely, the circles are elements of matrices A1,2 and A2,1, while the diamonds belong to the matrix A2,2.

The implemented controllers are feedback controllers. On the top and bottom the product streams are supervised by proportional controller maintaining the desired liquid level in the respective drums. The column is also equipped with a proportional-integral pressure control, which restricts the gas stream to the condenser to regulate the top pressure of the splitter. Finally there is a proportional-integral temperature controller implemented to control the temperature in the lower section manipulating the reboiler duty. Adding two PI controllers to the set of equations increases the number of the dynamic equations the DAE solver must integrate by two and inserts 14 unstructured elements into the well-structured tridiagonal block Jacobian, as shown in Figure 5. The Jacobian is still highly sparse, even though the few elements in the ellipses corrupt the compact tridiagonal structure. 6. Numerical Results To measure the performance of the new algorithm, two solvers are being used, namely, the full scheme and the advanced sparse solver. The distillation model is excited with an external input. Three cases were generated. (1) Procedure A: BzzDae is used to solve the system. Neither structure nor sparsity information is provided. (2) Procedure B: BzzDaeSparse is used. Information about the overall Boolean Jacobian (incidence matrix) is provided. (3) Procedure C: The new BzzDaeSparse is used providing information about the structure of the partially structured system. The location of the unstructured elements must also be provided. Distillation column models are well-known to be nonlinear. This implies that a +10% step perturbation in the feed composition could display different dynamics and nonsymmetrical responses than a -10% step perturbation. Six different input sequences were studied, and for each case the computational time for computing the process transient was recorded, both for the open-loop and closed-loop configuration. Starting from steady-state conditions, step variations were introduced in the system at time zero: the first 6 s of simulation were in open loop, while from 6 to 24 s the control loops were closed, with a sampling time of 6 s. Only this highly dynamic

Table 2. CPU Time Comparison, Open-Loop (in Seconds) disturbances, %

procedure A

procedure B

procedure C

+10 -10 +20 -20 +50 -50

7.92 7.86 7.91 7.92 7.98 15.78

5.89 1.89 1.92 1.92 2.16 2.16

0.47 0.44 0.47 0.48 0.52 0.94

Table 3. CPU Time Comparison, Closed-Loop (in Seconds) disturbances, %

procedure A

procedure B

procedure C

+10 -10 +20 -20 +50 -50

22.67 22.61 22.69 22.70 22.73 22.81

3.13 3.33 3.33 3.56 3.38 3.59

0.95 0.97 0.97 0.97 1.00 1.,06

regime was chosen, since the solvers perform nearly the same when the process is close to steady-state, and choosing longer simulation times is not adequate for measuring the solvers’ performances. The following six step variations were introduced in the feed composition of isobutane with respect to the nominal value: (10, (20, and (50%. Table 2 and Table 3 summarize the computation times (in seconds) required for the three procedures with the six step perturbations. It is well-known that a method accounting for system sparsity is much faster than the other methods that do not. Looking at the results, this trend is confirmed. Actually, procedure A, which solves the system without any sparsity information, requires about 8 s to simulate 6 s in open loop and about 23 s to simulate the remaining 18 s of the closed-loop scenario on an Intel Core 2 Duo CPU, 2.00 GHz, 2.00 GB of RAM. The integration is rather time-consuming, and it does not allow any application of predictive techniques and moving horizon methodologies, which requires online feasibility and CPU integration times faster than the real time. On the other hand, procedure B, where the user has to define the whole existence Jacobian matrix, is generally four times faster than procedure A. It requires about 2 s to simulate 6 s in open loop and about 3.5 s to simulate the closed-loop scenario.

9984

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

In turn, for all the studied step input signals, procedure B is four times slower than procedure C that does not only exploit the system sparsity and structure but also the location of the unstructured nonzero elements. By adopting the partially structured DAE solver, one needs a half-second to simulate open-loop scenarios and about 1 s for the closed-loop ones. The type of disturbance does not affect the computation time, even though there are three outliers in open-loop calculations, reported in Table 2. A -50% disturbance needs more or less twice the CPU time taken by the other step inputs in procedures A and C, while +10% in procedure C is circa three times slower than the time consumed by the other implemented step disturbances. This is mainly due to the size of the integration steps, which is automatically defined by the stiff solvers and strongly nonlinear process dynamics. When the system is operating in closed loop, controllers are able to stabilize process dynamics and, consequently, they prevent possible less performing solutions in terms of CPU requirements. 7. Conclusions This paper proposes and applies to an industrial case a highperformance algorithm to solve DAE systems with sparse Jacobian, whose elements are partially structured. This is a common case in process modeling and control applications, since the integral term of conventional control loops is usually an unstructured element which may spoil the overall welldefined system structure. On this subject, existing solvers can completely exploit neither the system structure nor, sometimes, the system sparsity, leading to the unavoidable increase of the computational time. On the other hand, the novel approach here developed and implemented in the BzzMath library can exploit the overall structured morphology of the system by requiring only a little extra information about the unstructured elements. It allows a significant improvement in computation time of numerical DAE system integrations, which is comparable to the numerical improvement between dense and sparse solvers in the solution of fully structured systems. Acknowledgment Founding through Mongstad Refinery’s Pilot Project (StatoilHydro) is thankfully acknowledged. We also thank Prof. Eliseo Ranzi (Politecnico di Milano) for valuable suggestions. Literature Cited (1) Gani, R.; Ruiz, C. A.; Cameron, I. T. A generalized model for distillation columnssI: Model description and applications. Comput. Chem. Eng. 1986, 10 (3), 181–198. (2) Krishnamurthy, R.; Taylor, R. Nonequilibrium stage model of multicomponent separation processes. AIChE J. 1985, 31 (12), 1973–1985. (3) Luyben, W. L. Process modeling, simulation and control for chemical engineers, 2nd ed.; Chemical Engineering Series; McGraw-Hill: New York, 1996. (4) Taylor, R.; Kooijman, H. A.; Hung, J. S. A second generation nonequilibrium model for computer simulation of multicomponent separation processes. Comput. Chem. Eng. 1994, 18 (3), 205–217. (5) Venkataraman, S.; Lucia, A. Solving distillation problems by Newton-like methods. Comput. Chem. Eng. 1988, 12 (1), 55–69. (6) Buzzi-Ferraris, G. BzzMath: Numerical library in C++, Politecnico di Milano, http://chem.polimi.it/homes/gbuzzi, 2009.

(7) Buzzi-Ferraris, G.; Manenti, F., Fundamentals and Linear Algebra for the Chemical Engineer SolVing Numerical Problems; WILEY-VCH: Weinheim, Germany, 2009; ISBN 978-3-527-32552-8. (8) Buzzi-Ferraris, G.; Manenti, F. Interpolation and regression models for the chemical engineer solVing numerical problems.; WILEY-VCH: Weinheim, Germany, 2009; ISBN 978-3-527-32652-5. (9) Busch, J.; Oldenburg, J.; Santos, M.; Cruse, A.; Marquardt, W. Dynamic predictive scheduling of operational strategies for continuous processes using mixed-logic dynamic optimization. Comput. Chem. Eng. 2007, 31 (5-6), 574–587. (10) Cutler, C. R.; Ramaker, B. L. Dynamic matrix controlsA computer control algorithm. Proceedings of the Joint Automated Control Conference, 1980. (11) Lang, Y. D.; Biegler, L. T. A software environment for simultaneous dynamic optimization. Comput. Chem. Eng. 2007, 31 (8), 931–942. (12) Manenti, F.; Rovaglio, M. Integrated multilevel optimization in large-scale poly(ethylene terephthalate) plants. Ind. Eng. Chem. Res. 2008, 47 (1), 92–104. (13) Morari, M.; Lee, J. H. Model predictive control: Past, present and future. Comput. Chem. Eng. 1999, 23 (4-5), 667–682. (14) Prett, D. M.; Garcı`a, C. E. Fundamental process control; Butterworth: Boston, MA, 1988. (15) Prett, D. M.; Gillette, R. D. Optimization and constrained multivariable control of a catalytic cracking unit. Proceedings of the Joint Automatic Control Conference, 1980. (16) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control Eng. Practice 2003, 11 (7), 733–764. (17) Rawlings, J. B. Tutorial overview of model predictive control. IEEE Control Syst. Mag. 2000, 20 (3), 38–52. (18) Lima, N. M. N.; Manenti, F.; Maciel Filho, R.; Embiruc¸u, M.; Wolf Maciel, M. R. Fuzzy Model-Based Predictive Hybrid Control of Polymerization Processes. Ind. Eng. Chem. Res. 2009, 48 (18), 8542–8550. (19) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical solution of initial-Value problems in differential-algebraic equations, 2nd ed.; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, 1996. (20) Brenan, K. E.; Petzold, L. R. The numerical-solution of higher index differential algebraic equations by implicit methods. Siam J. Numer. Anal. 1989, 26 (4), 976–996. (21) Hairer, E.; Nørsett, S. P.; Wanner, G., Solving ordinary differential equations IsNonstiff problems. Springer Series in Computational Mathematics, Vol. 8; Springer-Verlag: Berlin, 1987. (22) Hindmarsh, A. C. ODEPACK, a systematized collection of ODE solvers. Scientific Computing, Stepleman ed.; North-Holland: Amsterdam, 1983; pp 55-64. (23) Burrage, K.; Butcher, J. C.; Chipman, F. H. An implementation of singly-implicit Runge-Kutta methods. BIT (Nordisk Tidskrift fo¨r Informationsbahandling) 1980, 20, 326–340. (24) Brown, P. N.; Hindmarsh, A. C.; Petzold, L. R. Using Krylov methods in the solution of large-scale differential-algebraic systems. SIAM J. Sci. Comput. 1994, 1467–1488. (25) Gustafsson, K. An object oriented implementation of software for solving ordinary differential equations. Object oriented numerics conference, Sunriver, OR, 1993. (26) Petzold, L. R. Recent developments in the numerical solution of differential/algebraic systems. Comput. Methods Appl. Mech. Eng. 1989, 75, 77–89. (27) Dones, I.; Preisig, H. A.; de Graaf, S. Towards the dynamic initialisation of C4 splitter models; Proceedings of the 1st Annual Gas Processing Symposium, Doha, Qatar; 2009; p 1. (28) Krishnamurthy, R.; Taylor, R. A nonequilibrium stage model of multicomponent separation processes. Part I: Model description and method of solution. AIChE J. 1985, 31 (3), 449–456. (29) Løvfall, B. T. Computer realisation of thermodynamic models using algebraic objects. Ph.D. thesis, NTNU, Trondheim, Norway, 2008. (30) Preisig, H. A.; Haug-Warberg, T.; Løvfall, B. T.; Pantelides, W. M. a. C. On model portability. Comput.-Aided Chem. Eng. 2006, 1, 21.

ReceiVed for reView May 18, 2009 ReVised manuscript receiVed August 4, 2009 Accepted August 26, 2009 IE9007908