SciELO - Scientific Electronic Library Online

 
vol.38 número2A review on fault diagnosis of induction machinesAnalysis of heat and mass transfer in a mixed convective diffusion flame attached to a vertical fuel surface índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

  • No hay articulos citadosCitado por SciELO

Links relacionados

  • No hay articulos similaresSimilares en SciELO

Compartir


Latin American applied research

versión impresa ISSN 0327-0793

Lat. Am. appl. res. v.38 n.2 Bahía Blanca abr. 2008

 

Sensitivity analysis for differential algebraic equations

G. Boroni1, P. Lotito1 and A. Clausse1,2

1 CONICET and Universidad Nacional del Centro, Tandil, Argentina.
2 Comisión Nacional de Energía Atómica, Argentina.
{gboroni, plotito, clausse}@exa.unicen.edu.ar

AbstractIn this paper we propose a new method based on adjoint systems for parametric sensitivity analysis of DAE's. This method is employed in a series of experiments, and the results are compared with the estimation of sensitivity by a finite differences method, widely used because of the simplicity of its implementation . In addition, we propose a method extension based on the use of estimation functions that allows transforming the original adjoint system into an algebraic system, making it suitable for many efficient numerical methods.

Keywords — Simulation; Parametric Sensitivity; DAE System.

I. INTRODUCTION

The mathematical models that are currently being used to investigate physical phenomena are becoming more and more realistic. The new characteristic of these models is that they often use parameters whose values may not be accurately known calling for parametric sensitivity analysis. Areas of application include optimization, parameter estimation, optimal control, model simplification, process sensitivity, uncertainty analysis, and experimental design for a wide range of scientific and engineering problems.

Recent work on methods and software for sensitivity analysis of differential algebraic system (DAE) (Cao et al., 2003; Feehery et al., 1997; Li and Petzold, 2000; Li and Petzold, 1999; Li et al., 2000; and Maly and Petzold, 1997), has demonstrated that subsequent sensitivity values can be computed reliably and efficiently via automatic differentiation in combination with DAE solution techniques.

However, there are some difficult issues when trying to make sensitivity analysis of DAE systems with a large number of parameters compared with the number of variables. In these cases, the extension of the original system with the equations required for the sensitive analysis makes the resulting system grow significantly, and therefore, the calculation efficiency decreases.

The first contribution of this paper is the employment of a new method based on adjoint systems for parametric sensitivity analysis, whose theory was recently proposed for differential algebraic equations (Cao et al., 2003). Following this method, it has been proved that it is possible to calculate the sensitivity of parameters using linear adjoint systems, with no need to extend the original system, what allows improving the calculation efficiency considerably.

This method is used in a series of experiments, where the results are compared with the estimation of sensitivity by the finite differences method, widely used for the simplicity of its implementation.

The second contribution of this paper is to propose an extension of the method based on the use of estimation functions that allows transforming the original adjoint system into an algebraic system, making it suitable for the implementation of many efficient numerical methods.

II. SENSITIVITY CALCULATION FOR DIFFERENTIAL ALGEBRAIC SYSTEMS BY ADJOINT METHOD

The numerical solution of initial value problems modeled by DAE has attracted the interest of the numerical community for the last 30 years (Brenan et al., 1996 and Kees and Millen, 1999). Many engineering and scientific problems can be naturally modeled with DAE's, mostly for balance equations of the general form:

(1)

composed by a combination of differential and algebraic equations, where x is the state variable vector, and p is a parameter vector.

The problem of calculating the sensitivity of parameters of a DAE system takes the following form: for a parameter dependent DAE system of type (1), find at time T, for j : 1,..,Np.

Their solution requires the simultaneous solution of the original DAE system with the Np sensitivity systems obtained by differentiating the original DAE with respect to each parameter in turn (Cao et al., 2003).

For large systems this may look like a lot of work, but it can be done efficiently, if Np is relatively small, by exploiting the fact that the sensitivity systems are linear and all share the same Jacobian matrices with the original system (Li et al., 2000).

However, some problems require calculating the sensitivities with respect to a large number of parameters. For these problems, particularly if the number of state variables Nx is also large, the previous sensitivity approach is out of reach.

In Cao et al. (2003) it was demonstrated that these problems can often be handled more efficiently by the adjoint method (Errico 1997). In this approach, we are interested in calculating the sensitivity of an objective function of the perturbed solution and the parameter:

(2)

or alternatively the sensitivity of a function g(x,t,p) defined only at time T. For example, in the mathematical models of nuclear reactors, a common analysis is calculating the sensitivity of speed ui of the cooling liquid, product of the perturbation of the demanded thermal power Q. Thus, the following objective function can be formulated:

(3)

where ū is the average of the speed of the liquid. Next, we describe how to employ the adjoint method to compute the sensitivity.

A. Sensitivity of G(x,p)

To find an expression for , we define the augmented objective function I(x,p) as:

(4)

where λ is a Lagrange multiplier, and F(x,,p,t)dt = 0 by (1). Formally differentiating both members with respect to p, the sensitivity of G(x,p) with respect to p is:

(5)

Then, using integration by parts we obtain an equivalent expression for the term depending on xp instead:

(6)

If the Eq. (6) is replaced in (5) and the terms regrouped, we obtain that:

(7)

and if we require that the following (adjoint) equation be satisfied

(8)

Eq. (7) becomes

(9)

Note that (8) is the adjoint equation, and λ is the adjoint variable of the system. In this work we used the notation of augmented adjoint system corresponding to (8):

(10)

where we have introduced the function α(t).

From (10) onwards, it is observed that the adjoint DAE must be solved backward in time. To calculate the integral term of the Eq. (9), we use a quadrature variable β, expanding again the adjoint system:

(11)

In this way, for t=0,

(12)

yielding the sensitivity equations for :

(13)

B. Sensitivity of g(x,t,p)

Now let us consider the computation of . Interchanging differentiations we have that and from (9), we obtain:

(14)

where the adjoint variable λ(t,T) depends on both t and T, and λT(t,T) denotes . The corresponding adjoint equation is:

(15)

and it is solved in a similar way to system (10).

Again, to calculate the integral term of Eq. (14), we used a quadrature variable β, and we defined the following adjoint system:

(16)

Finally, the expression of the sensitivity equation is given by

(17)

III. NUMERICAL EXPERIMENTS

Nowadays, there is a variety of numerical methods for solving IVP 's given by DAE's, which can be classified into three classes: single-step, multistep, and extrapolation methods. The efficiency of each type of approximation depends on the specific problem characteristics.

In this work the systems (11) and (16) are solved transforming the DAE 's into purely algebraic systems by means of numerical approximations computing through estimating functions as in Boroni and Clausse (2005). In order to do that, let us define the functions and , that represent multistep estimators of the state vector and its derivatives, i.e., λ and . Thus:

(18)

Specifically, in this work we use the adjoint variable and implicit linear multistep estimator (Jackson and Sacks-Davis, 1980) for the derivatives:

(19)

To analyze the application of the adjoint method in the calculation of sensitivity, a series of experiments was made, where the numerical results and errors were examined. In this study two systems were defined: a linear differential system with explicit solution, and a nonlinear conservative oscillating system. So as to simplify the presentation of the results, only the calculation of the sensitivity of the objective function G(x,p) is considered.

A. Linear System

An example of a linear differential equation of the first order is set out, given by:

(20)

This system has only one parameter p, and the explicit solution is:

(21)

The purpose of working with this system was to study the error between the explicit solution and the numerical solution of the sensitivity computed as in (13), for the following objective function

(22)

Using (13) and (22), we obtain that:

(23)

The temporary evolution of the sensitivity can be seen in Fig. 1, where it has the form of a left-biased distribution.


Fig. 1. Temporal evolution of the sensitivity .

In Fig. 2 we can see the absolute error ε between the exact value and the computed value of the sensitivity dG/dp with respect to the explicit solution. In this figure it can be seen that the calculation of the dG/dp generates a error with a maximum close to 1.1x10-4, which stabilizes for a time t ≥ 3.


Fig. 2. Calculation of the absolute error for .

B. Elastic pendulum

The second example is the elastic pendulum (Fig. 3), consisting of a no-linear differential equation system of the fourth order whose natural variables are the string length, the inclination angle with respect to the vertical, and their respective temporal derivatives, that is:

(24)


Fig. 3. Elastic Pendulum.

Since in this case, no explicit solution is available, relative comparisons were made between the explicit and implicit numerical solutions for the following parameters and initial conditions:
k = 7 , L = 1 , g = 8 , r(0) = 1, θ(0) = , (0) = 0 y (0) = 0 .

Figure 4 shows the trajectory of the mass in the ( x, y) plane.


Fig. 4. Trajectory of the mass in the (x, y) plane .

For the calculation of the sensitivity we defined an objective function related with the energy of the pendulum:

(25)

where the term is associated to the harmonic potential energy, mgr(1-cos(θ)) represents the gravitational potential energy and is associated to the kinetic energy. Initially, when the pendulum is standing, the total energy is:

(26)

The purpose of this example is to determine an error 's measure, considering that the objective function has an explicit solution given by the principle of conservation of the energy. Therefore, the following objective function is defined:

(27)

Developing (27) for the pendulum in movement, we obtain that:

(28)

where v2 = 2 + r22 . When the pendulum is standing, the sensitivity is:

(29)

obtaining as a result a constant value equal to 0 (r0 = L). With respect to parameter m, the sensitivity is expressed by:

(30)

For (30) the results indicate that sensitivity grows linearly with time. Using the sensitivity expressions (29) and (30), we made a comparative analysis of the numerical results obtained employing the adjoint method and the finite differences method. In the finite difference method, the system DAE is solved numerically with different values of parameter: a value of perturbed parameter and the rest of the parameters with original values. Thus, it is possible to obtain for example an estimation of the sensitivity equal to:

(32)

The advantage of this method is the simplicity of its implementation, whereas the main disadvantage is that the estimation error is significant.

Figure 5 shows the error defined by the absolute difference εDGDK between the explicit solution and numerical solution. In this figure it can be seen that when adopting the adjoint method the error is close to 1, whereas when adopting the finite differences method, the error is tripled with respect to the adjoint method.


Fig. 5. Absolute error of .

Figure 6 shows the error defined by the absolute difference εDGDM, where the absolute error calculated by the adjoint method is very small if we compare it with the error obtained by finite differences. It can also be observed that for values of time T between 7.7 and 7.9, the sensitivity by finite differences reaches values between 125 and 130, with an εDGDM error between 55 and 60, what practically represents a 50%. This is a clear indicator that the adoption of the finite difference method can produce wrong results.


Fig . 6. Absolute error of .

IV. CONCLUSIONS

It can be seen in the previous examples that the calculation of sensitivity is very often obtained with quite small precision errors. For the linear system (III.A), the proposed method allows obtaining a reasonably accurate solution. When the DAE system is nonlinear (III.B), the calculation of parametric sensitivity by the adjoint method produces a quite acceptable solution.

Finally, the most important advantage of the proposed method is the possibility of calculating the sensitivity of parameters of DAE 's using linear adjoint systems, with no need to extend the original system, thus improving the efficiency in the calculation.

REFERENCES
1. Boroni, G. and A. Clausse, "Una formulación matemática orientada a objetos para simulación continua", Mecánica Computacional, 24, 261 (2005).         [ Links ]
2. Brenan, K., S. Campbell and L. Petzold, Numerical Solution of Initial Value Problems in Differential-Algebraic Equations, SIAM, Philadelphia, PA (1996).         [ Links ]
3. Cao, Y., S. Li, L. Petzold and R. Serban, "Adjoint sensitivity analysis for differential-algebraic equations: The adjoint DAE system and its numerical solution", SIAM J. Sci. Comput., 24, 1076-1089 (2003).         [ Links ]
4. Errico, R.M., "What is an adjoint model?", Bulletin of the American Meteorological Society, 78, 2577-2591 (1997).         [ Links ]
5. Feehery, W.F., J.E. Tolsma and P.I. Barton, "Efficient sensitivity analysis of large-scale differential-algebraic systems", Appl. Numer. Math., 25, 41-54 (1997).         [ Links ]
6. Jackson, R. and R. Sacks-Davis, "An Alternative Implementation of Variable Step-Size Multistep Formulas for Stiff ODEs", ACM Trans. on Mathematical Software, 6, 295-318 (1980).         [ Links ]
7. Kees, C. and C. Miller, "C11 Implementations of Numerical Methods for Solving Differential-Algebraic Equations: Design and Optimization Considerations", ACM Trans. on Mathematical Software, 25, 377-403 (1999).         [ Links ]
8. Li, S. and L.R. Petzold, "Design of New DASPK for Sensitivity Analysis", Technical report, Dept. of Computer Science, UCSB, Santa Barbara (1999).         [ Links ]
9. Li, S. and L.R. Petzold, "Software and algorithms for sensitivity analysis of large-scale differential-algebraic systems", J. Comput. Appl. Math., 125, 131-145 (2000).         [ Links ]
10. Li, S., L.R. Petzold and W. Zhu, "Sensitivity analysis of differential-algebraic equations: A comparison of methods on a special problem", Appl. Numer. Math., 32, 161-174 (2000).         [ Links ]
11. Maly, T. and L.R. Petzold, "Numerical methods and software for sensitivity analysis of differential-algebraic systems", Appl. Numer. Math., 20, 57-79 (1997).
        [ Links ]

Received: April 12, 2007.
Accepted: August 14, 2007.
Recommended by Subject Editor: Jorge Solsona

Creative Commons License Todo el contenido de esta revista, excepto dónde está identificado, está bajo una Licencia Creative Commons