روش ESDIRK با قابلیت های تجزیه و تحلیل حساسیت

A new algorithm for numerical sensitivity analysis of ordinary differential equations (ODEs) is presented. The underlying ODE solver belongs to the Runge–Kutta family. The algorithm calculates sensitivities with respect to problem parameters and initial conditions, exploiting the special structure of the sensitivity equations. A key feature is the reuse of information already computed for the state integration, hereby minimizing the extra effort required for sensitivity integration. Through case studies the new algorithm is compared to an extrapolation method and to the more established BDF based approaches. Several advantages of the new approach are demonstrated, especially when frequent discontinuities are present, which renders the new algorithm particularly suitable for dynamic optimization purposes.

Systematic parameter estimation, design and control of systems modelled by ordinary differential equations often require knowledge of sensitivities (Russel, Henriksen, Jørgensen, & Gani, 2000). In addition sensitivity analysis is used for experimental design and model reduction as well as in bifurcation analysis (Jørgensen & Jørgensen, 1998). Enabling efficient and robust sensitivity calculations thus constitutes an important prerequisite for most systematic process and product engineering disciplines (Braatz et al., 2004). The present paper addresses the development of a methodology to achieve a fast, efficient and reliable sensitivity calculation. Nonlinear moving horizon estimation and control also requires efficient and robust algorithms for integration and sensitivity calculations Bauer, 2000, Grötschel et al., 2001 and Jørgensen et al., 2003 Rambau, 2001; Jørgensen, Rawlings, & Jørgensen, 2003). The dual requirement of both efficiency and robustness stems from the online, real-time and unsupervised implementation of nonlinear moving horizon estimation and control in nonlinear model predictive control applications. The sensitivities are needed in the evaluation of gradients required by sequential quadratic programming algorithms. The key contributions of this paper are: • presentation of an ESDIRK numerical integration method that can be applied to ODE as well as index-1 DAE systems. The method can be applied to stiff as well as nonstiff systems. The method is particularly efficient for methods with frequent discontinuities, • construction of a continuous extension which means that the method can be applied to hybrid systems and discrete event systems, • derivation of an efficient method for computation of the state and parameter sensitivities, • description of an implementation (ESDIRK34) of the method that efficiently integrates a system of ODEs along with computation of the sensitivities. In this paper a Runge–Kutta based algorithm for the numerical computation of sensitivities in systems of ODEs is presented. The algorithm has been implemented in the code ESDIRK34 and a number of tests have been conducted comparing the performance of the code with other codes for sensitivity analysis. This comparative study is the topic of Section 4. Section 2 concerns the mathematical basis, while Section 3 covers some important implementation issues including error and convergence control. The considered initial value problem is equation(1) View the MathML sourcex˙=f(t,x,u),x(t0)=x0(u) Turn MathJax on in which View the MathML sourcex∈Rns is the nsns-dimensional state vector and View the MathML sourceu∈Rnp is an npnp-dimensional vector of parameters. In addition to integrating the state equations, ESDIRK34 also computes first order derivative information of the solution with respect to both state variables and parameters. The state sensitivity of the states View the MathML sourcex with respect to the initial state x0ix0i is defined as equation(2) View the MathML sourcess,i=∂x(t,u)∂x0i,i=1,…,ns Turn MathJax on View the MathML sourcess,i satisfies the following sensitivity equations equation(3) View the MathML sources˙s,i=∂f∂xss,i,ss,i(t0)=ei Turn MathJax on where View the MathML sourceei is the ii th column of the nsns-dimensional identity matrix. Furthermore we define the parameter sensitivity of the states View the MathML sourcex with respect to the parameter uiui as equation(4) View the MathML sourcesp,i=∂x(t,u)∂ui,i=1,…,np Turn MathJax on View the MathML sourcesp,i satisfies the following sensitivity equations equation(5) View the MathML sources˙p,i=∂f∂xsp,i+∂f∂ui,sp,i(t0)=∂x0∂ui Turn MathJax on The existence of the sensitivities is given by Gronwall’s theorem, which states that provided the partial derivatives View the MathML source∂f/∂x and View the MathML source∂f/∂u exist and are continuous in a neighbourhood of the solution View the MathML sourcex(t), then state and parameter sensitivities exist, are continuous, and satisfy (3) and (5), respectively (Hairer, Nørsett, & Wanner, 1991). The desire to obtain first order derivative information has extended the original problem of nsns equations to a combined system of nsns state equations and nsns+npnsnsns+npns sensitivity equations. Solving this system directly would be highly inefficient. Instead, we show in this paper how the sensitivity information can be obtained with little extra effort. Solution of (1) in the interval [tk,tk+1]=[tk,tk+Ts][tk,tk+1]=[tk,tk+Ts] with the initial condition, View the MathML sourcexk, and the parameter, View the MathML sourceuk, may be regarded as a discrete map equation(6) View the MathML sourcexk+1=Hk(xk,uk) Turn MathJax on The sensitivities obtained by (3) with respect to the state vectors, View the MathML sourcexk, and (5) with respect to the parameters, View the MathML sourceuk, may in this difference equation interpretation be denoted as equation(7a) View the MathML sourceAk=∇xkHk(xk,uk) Turn MathJax on equation(7b) View the MathML sourceBk=∇ukHk(xk,uk) Turn MathJax on Optimization algorithms for nonlinear model predictive control are developed using this discrete time view of the integration method and the sensitivities Allgöwer et al., 1999, Allgöwer and Zheng, 2000 and Jørgensen et al., 2003 Qin, Rawlings, & Wright, 1999; Allgöwer & Zheng, 2000; Jørgensen et al., 2003). Hence, ESDIRK34 is a method for integration and computation of the sensitivities needed in nonlinear model predictive control applications.

A new algorithm for sensitivity analysis of ODEs has been presented. The algorithm is based on an ESDIRK method of the Runge-Kutta family, and it has been implemented in the code ESDIRK34. The algorithm has order 3 and is A- and L-stable as well as stiffly accurate. The state and sensitivity integrations are performed separately, thereby enabling the linearity of the sensitivity equations to be exploited. A key feature of the new algorithm is the reuse of the Jacobian evaluations for both iteration matrices and sensitivity residuals. In this way the extra effort of calculating sensitivities is minimized. Case studies have shown that clear advantages exist for the use of one-step methods such as ESDIRK34, when frequent discontinuities are present in the solution. The multistep method DASPK proved inefficient compared to ESDIRK34 for this type of problem. This makes ESDIRK34 suitable for use in dynamic optimization and nonlinear model predictive control. Future work should include extension of the solver to handle algebraic equations, which it is already prepared for due to the good stability properties and high internal stage order of the ESDIRK method. Furthermore, inclusion of sparse techniques for large scale problems would be a desirable extension as well as further numerical testing. A more subtle extension could be to experiment with and improve the error controller using control theory as ingeniously pioneered by Gustafsson (1992).