A modified collocation method for solving differential-algebraic equations

https://doi.org/10.1016/S0096-3003(99)00138-1Get rights and content

Abstract

In this study, a modified collocation method is introduced for solving nonlinear initial value problems and their sensitivity profiles of mixed differential and algebraic equations. The weighted residual equation in the modified method is defined based on the integral expression of the problem. The additional equality constraints for enforcing continuity of the solution at each element boundary are noticeably not required in this study. As a result, the collocation equations in global and finite element time domain are able to unify into a compact formula so that coefficients of the solutions to differential-algebraic equations can be obtained, element by element. The proposed method can be straightforwardly applied for handling sensitivity matrix equations of differential-algebraic problems. The solution and sensitivity profiles of the differential-algebraic equations are simultaneously evaluated by this approach.

Introduction

The mathematical models of process engineering systems operating under transient conditions are usually described by mixed sets of differential and algebraic equations (DAEs). Differential equations in such systems arise from the unsteady state material and energy balances, whereas the algebraic equations commonly result from pseudo-steady-state assumptions, equilibrium relations and process constraints, and so on. In realistic processes, these equations are usually nonlinear, and one must resort to appropriate numerical methods to obtain the desired time profiles. Many, but not all, of DAEs can be solved conveniently using numerical codes of ordinary differential equations (ODEs) [1]. Historically, a set of DAEs was solved as ODEs often destroying the natural structure of the system 2, 3. This was done primarily to put the problem in the form required by the particular integration method or software package being used. Today, it is becoming more common to deal with original DAEs. One important reason for this is that the variables in the original DAEs typically have some physical significance, whereas those that result after manipulation into ODEs may not [4]. Thus by solving the original DAEs directly, the effect of parameter changes or model changes is more easily studied.

Several papers have appeared in the past decade which are related to the solution of DAEs. Brenan et al. [4] review more than 200 papers on the solution of DAEs and provide a detailed description of methods for the solution of DAEs, including convergence and solvability analysis. Popular DAE solution codes such as DASSL [5], are based on the use of the well-known backward differentiation formulas, as first suggested by Gear [6]. Newton-like Runge–Kutta methods are the other common algorithms for solving DAEs 7, 8. Collocation methods have been used for more than thirty years to solve ordinary differential equations [9]. In the literature, such conventional collocation methods are rarely employed to solve DAEs. In cases of dynamic optimization problems for DAEs, collocation methods are used only to transform dynamic equations into algebraic ones 10, 11, 12. A nonlinear programming method is then applied to determine an extreme solution of the approximate problem.

In this study, a modified collocation method is developed to solve DAEs. This method can be straightforwardly applied to solve sensitivity equations of DAEs. The weighted residual equation is defined based on the integral expression of DAEs. As a result, the approximate equations in the global and local time domain can be unified into a compact formula. Therefore, the solution and sensitivity profiles of DAEs can be simultaneously computed, element-by-element.

Section snippets

Collocation on global time domain

Consider a dynamic system, described by the following set of differential and algebraic equations:dxdt=f(x(t),y(t),t;θ),t0⩽t⩽tf,h(x,y,t;θ)=0,where x(t)≡[x1,x2,…,xn]TRn,y(t)≡[y1,y2,…,yr]TRr are referred to as “dynamic-state” and “steady-state” variables. The superscript T denotes the transpose. The parameter vector, θ∈Rp, is considered to be held constant in evaluating the solution of this problem. The notation of this parameter vector is dropped in the following sections for notational

A unified formula

The accuracy and efficiency of collocation methods on global time domain are largely dependent upon the degree of shape polynomials. More computational efforts required to solve the approximate problem arising by the higher-degree polynomials. Although there is no theoretical limit to the degree of the polynomial one could employ as a shape function. An alternative to collocation methods on global time domain uses piecewise polynomial approximations towards reducing the degree.

The whole time

Numerical computational algorithm

The remained task in this study is to find the expansion coefficients of the state variables in , , , . These equations can be solved by Newton–Raphson method or its modified methods. There are two phases to find these coefficients. In the first phase, the initial conditions of DAEs are determined to obtain the consistent initial values. In the second phase, , are then used to find expansion coefficients of the approximate problem, element-by element.

The Newton–Raphson method with step length

Sensitivity analysis

Sensitivity analysis provides a systematic framework in which to study the accuracy and robustness of a mathematical model. Sensitivity matrices, which are partial derivatives of the state variables respect to the system parameters, play an important role in uncertainty analysis, parameter estimation and optimization, experimental data analysis, and model discrimination. Several methods exist to solve sensitivity equations derived from systems of ODEs. Extensive reviews can be found in Rabitz

Applications

The computational schemes were tested on a wide variety of problems. Two reactor examples were presented here. MATLAB and its toolboxes were used to implement the proposed algorithms and to test it with these reactor problems.

Example 1. The free radical polymerization of styrene in a batch reactor catalyzed by mixed initiators is considered in this example. This system was modeled by a set of implicit formulas of DAEs discussed in Butala [18]. The modeling equations are expressed as follows:dx1d

Conclusions

A modified collocation method have proposed for solving DAEs in this study. The weighted residual equation is defined based on the integral expression of DAEs. The proposed method compared with the conventional collocation method is also discussed in detail. The additional equality constraints for enforcing continuity of the solution at each element boundary are noticeably not required in this study. As a result, the approximate equations for DAEs in the global and local time domain can be

Acknowledgements

The author would like to thank the National Science Council of the ROC for financially supporting this research under Contract No. NSC84-2214-E194-005.

References (20)

There are more references available in the full text version of this article.

Cited by (19)

  • An efficient sparse approach to sensitivity generation for large-scale dynamic optimization

    2011, Computers and Chemical Engineering
    Citation Excerpt :

    In chemical engineering, general DAEs often correspond to the natural statement of the underlying physical models. A reformulation or index reduction is often not convenient and for large nonlinear problems may even not be possible (Cao et al., 2002; Logsdon & Biegler, 1989; Wang, 2000). Thus, solvers able to handle general DAEs present basically some advantages.

  • Inverse problems of biological systems using multi-objective optimization

    2008, Journal of the Chinese Institute of Chemical Engineers
  • Kinetic modeling using S-systems and lin-log approaches

    2007, Biochemical Engineering Journal
View all citing articles on Scopus
View full text