Differential-algebraic Equation Solver Research Articles

A Sparse Differential Algebraic Equation (DAE) and Stiff Ordinary Differential Equation (ODE) Solver in Maple

This paper implements efficient numerical methods in Maple to solve index-1 nonlinear Differential Algebraic Equations (DAEs) and stiff Ordinary Differential Equations (ODEs) systems. Single-step methods (like Trapezoid (TR), Implicit-mid point (IMP), Euler-backward (EB), Radau IIA (Rad) methods, TRBDF2, TRX2) and backward-difference formula of order 2 are implemented with adaptive time-stepping methods in Maple to solve index-1 nonlinear DAEs. Maple’s robust and efficient ability to search within a list/set is exploited to identify the sparsity pattern and automatically calculate the analytic Jacobian. The algorithm and implementation are robust and efficient for index-1 DAE problems and scale well for finite difference/finite element discretization of two-dimensional models with system size up to 10,000 nonlinear DAEs and solve the same in a few seconds. The computational efficiency of the proposed algorithm (provided as an open-access code) compares favorably with the commercial solver gPROMs, one of the most commonly used sparse DAE solvers in the industry.

BioModME for building and simulating dynamic computational models of complex biological systems.

Molecular mechanisms of biological functions and disease processes are exceptionally complex, and our ability to interrogate and understand relationships is becoming increasingly dependent on the use of computational modeling. We have developed "BioModME," a standalone R-based web application package, providing an intuitive and comprehensive graphical user interface to help investigators build, solve, visualize, and analyze computational models of complex biological systems. Some important features of the application package include multi-region system modeling, custom reaction rate laws and equations, unit conversion, model parameter estimation utilizing experimental data, and import and export of model information in the Systems Biology Matkup Language format. The users can also export models to MATLAB, R, and Python languages and the equations to LaTeX and Mathematical Markup Language formats. Other important features include an online model development platform, multi-modality visualization tool, and efficient numerical solvers for differential-algebraic equations and optimization. All relevant software information including documentation and tutorials can be found at https://mcw.marquette.edu/biomedical-engineering/computational-systems-biology-lab/biomodme.php. Deployed software can be accessed at https://biomodme.ctsi.mcw.edu/. Source code is freely available for download at https://github.com/MCWComputationalBiologyLab/BioModME.

A Fundamental Study for Numerical Integration Procedures for the ANCF Consistent Rotation-Based Formulation Beam Elements

This study presents a numerical integration technique for a shear deformable beam element based on the absolute nodal coordinate formulation (ANCF). In particular, the consistent rotation-based formulation (CRBF) is used to develop a plane ANCF beam element, namely, the ANCF/CRBF shear deformable element. Even though the ANCF/CRBF beam element introduces rotation parameters as nodal coordinates, it does not need for interpolating the rotation field. Therefore, the ANCF/CRBF formulation defines a unique rotation field, and does not suffer from a coordinate redundancy problem encountered when using the large rotation vector formulation (LRVF). On the other hand, ANCF/CRBF formulation leads to highly nonlinear inertia forces including non-zero Coriolis and centrifugal forces. The objective of this study is to develop a constrained formulation and a numerical integration technique for the ANCF/CRBF shear deformable beam element. In particular, this study focuses on the fact that the global position gradients can be expressed in terms of the finite rotation parameters. In the present method, such algebraic relations are treated as constraints, and the ANCF/CRBF beam element is formulated as the differential algebraic equations (DAEs). As a result, the present formulation ensures the constant mass matrix, however, DAE solvers have to be introduced for solving the resulting equations. In addtion, since the present method can be regarded as a kind of augmented formulation, the number of degrees of freedom is significant from a standpoint of numerical efficiency. In order to cope with increase of the degrees of freedom, this study develops size reduction techniques focusing on the algebraic expressions of constraints. Then, the results obtained from the proposed method are compared with the results obtained using the fully parameterized ANCF shear deformable element.

Simulation of differential-algebraic equation systems with optimization criteria embedded in Modelica

Differential-algebraic equations with embedded optimization criteria (DAEO) are a class of mathematical models for underdetermined differential-algebraic equation (DAE) systems with less algebraic equations than algebraic variables. The algebraic variables may be calculated as the solution of an embedded (non)linear program, yielding a DAEO system. An example for DAEOs is the dynamic flux balance analysis (DFBA) approach, where the formulation of metabolic reaction networks leads to an underdetermined equation system for the intracellular fluxes that are assumed to behave optimally with respect to some cell-specific optimization criterion.We present a toolbox that allows formulation of DAEOs in the object-oriented Modelica modeling language. The solution method is based on substituting the embedded optimization problem with its first-order Karush-Kuhn-Tucker conditions to obtain a nonsmooth DAE system that can be simulated by a root-finding DAE solver. One nonlinear example and two examples based on DFBA demonstrate the performance of the toolbox.

Application of the DAE approach to the nonlinear sloshing problem

The nonlinear sloshing problem is an important issue for design of liquid storage tanks, liquid cargo transportations, tuned liquid dampers and so on. This paper is concerned with development of a novel approach for the nonlinear sloshing problem based on the Hamiltonian mechanics. In particular, this study is aimed at developing the method available to analyze the nonlinear liquid surface behavior like a traveling wave observed in the small liquid depth. In the present formulation, the fluid is assumed to be inviscid incompressible and irrotational flow. Then, the liquid surface motion is described by nonlinear multimodal models. However, since the sloshing problem based on such assumptions yields an irregular Lagrangian, it makes formulation difficult in a straightforward way. Therefore, the present approach employs the constrained Hamiltonian mechanics with the Lagrange’s method of undetermined multipliers to derive equations of motion. The resulting system is comprised of differential equations and algebraic equations, referred to as differential algebraic equations (DAEs). In addition, the present method takes full account of the nonlinear mode-to-mode interactions without reduction methods focusing on the predominant sloshing modes. However, the multimodal models without such reduction methods suffer from severe numerical stiff problem. Therefore, the numerical integration techniques based on implicit schemes (DAE solver) are incorporated as remedies for the stiff problem. Specifically, discrete forms of the equations of motions are derived by employing the Galerkin method and a discrete derivative. The proposed approach is validated by comparisons with an existing model and an experiment in time domain analysis.

Optimisation of dynamic reactive power sources using mesh adaptive direct search

Dynamic reactive power sources can be used to effectively mitigate the fault-induced delayed voltage recovery and transient voltage instability issues. When many var sources need to be installed at planned locations, optimisation of their sizes is a complicated non-linear optimisation problem due to its non-convexity and the dependence of the constraint on time-series trajectories of post-fault voltages. Solving this optimal sizing problem needs to utilise both a non-linear optimisation solver and a power system differential-algebraic equation solver. This study proposes a new approach for solving this problem under both a single contingency and multiple contingencies by using the mesh adaptive direct search algorithm interfaced with a power system simulator. The proposed approach is validated by case studies on a North American Eastern Interconnection model to optimise the sizes of planned STATCOMs against critical contingencies.

Comparison of several formulations and integration methods for the resolution of DAEs formulations in event-driven simulation of nonsmooth frictionless multibody dynamics

This article is devoted to the comparison of numerical integration methods for nonsmooth multibody dynamics with joints, unilateral contacts and impacts in an industrial context. With an event–driven strategy, the smooth dynamics, which is integrated between two events, may be equivalently formulated as a Differential Algebraic Equation (DAE) of index 1, 2 or 3. It is well-known that these reformulations are no longer equivalent when a numerical time–integration technique is used. The drift-off effect and the stability of the numerical scheme strongly depend on the index of the formulation. But, besides the standard properties of accuracy and stability of the DAE solvers, the event–driven context imposes some further requirements that are crucial for a robust and efficient event-driven strategy. In this article, several state–of–the–art numerical time integration methods for each formulation are compared: the generalized-α scheme for index-3 formulation and stabilized index-2 formulation, (Partitioned) Runge–Kutta Half–Explicit Method of order 5 (HEM5 and PHEM56) for index-2 DAEs with projection techniques, and Runge-Kutta explicit scheme of order 5, the Dormand-Prince scheme (DOPRI5), for index-1 DAEs with projection techniques (MDOP5). We compare these schemes in terms of efficiency, violation of the constraints and the way they handle stiff dynamics on numerous industrial benchmarks, where a CAD software is in this loop. One of the major conclusions is that the index-2 DAEs solvers prove to be better than other schemes to maintain low violations at position and acceleration levels. The best compromise allows us to design efficient event–driven solvers. When the dynamics is stiff, implicit schemes outperform explicit and half-explicit methods whichare sometimes unable to compute the dynamics when the system’s frequency range is wide. Furthermore, in industrial context, some solvers fail to reproduce the properties that they enjoy in theory. This is particularly true for half-explicit schemes when the Jacobian of the constraints has not full rank.

Energy Transfer Efficiency Optimization in an Electromagnetic Railgun

Electromagnetic railguns (EMRs) launch projectiles by the systematic triggering of current pulse forming units (PFUs), which release stored capacitor energy. Past work has focused on maximizing projectile velocity at the muzzle for given PFU capacitor charge voltages, which maximizes efficiency. However, operators in naval applications need to be able to select the muzzle velocity for a projectile to achieve a desired ballistic path. Thus, there is a demand to maximize efficiency for a given muzzle velocity. As such, this paper investigates maximizing the efficiency of transferring the stored PFU capacitor electrical energy to projectile kinetic energy while minimizing the current value when the projectile exits the rails to reduce damage to them from arcing. This is accomplished by optimizing the choice of the initial capacitor voltage in each PFU and the trigger timing subject to an electrical-mechanical dynamic model of the EMR and projectile and any additional scenario constraints (described shortly). This optimization is solved using numerical programming in the context of a relaxation with the initial capacitor voltages and triggering times as the only explicit variables to the optimizer; additional needed values are obtained from a differential algebraic equation solver. The optimization is successfully demonstrated for three scenarios: operation for a projectile mass and muzzle velocity associated with published operating parameters, operation for several projectile mass-muzzle velocity pairs with both equal and distinct initial capacitor voltages, and operation for several projectile mass-muzzle velocity pairs with a load current bound and equal initial capacitor voltages. Results show improved efficiency in comparison to the use of the published parameters, the independence of efficiency from projectile mass when there is no current bound, virtually identical results for equal and distinct initial capacitor voltages, and decreased efficiency with a current bound.

Safety Verification Utilizing Model-based Development for Safety Critical Cyber-Physical Systems

The application of cyber-physical systems (CPSs) in safety-critical application domain requires rigorous verification of their functional correctness and safety-relevant properties. We propose a practical verification process which enables to conduct safety verification of safety critical CPSs. The verification process consists of (a) a system model construction method, which generates a system model by combining software described in C and plant model code reused from model-based development, (b) a model transformation method, which transforms the plant models including differential algebraic equations (DAE) to approximate models without DAE to reduce verification complexity induced by DAE solver execution, (c) a model simplification framework, which enables the simplification of bond-graph plant models using domain-knowledge-based replacement of complex model components for further verification overhead reductions, and (d) a formal verification based on symbolic execution. We implemented the proposed methods and framework, and successfully applied the proposed verification process for safety verification of automotive brake control systems. The results of the study demonstrate that the verification detects a complex failure condition in a real-world brake control system from the generated system model and that the automated model transformations of the CPS models yield significant verification complexity reductions without impairing the ability to detect unsafe behavior. © 2017 Information Processing Society of Japan.

Innovative process simulation software not only for electromembrane processes

An effective tool is required for a successful and effective development and design of technologies. It is often necessary to propose several conceptual designs for the particular application. To assess and compare all designs, the process simulation software is usually used. Especially, for the cases of complex technologies which consist of several unit operations. There are several commercial products allowing comprehensive design, simulation, and optimization, such as Aspen HYSYS® AspenTech (commercial product) or DWSim (open source). Nevertheless, none of them includes directly electromembrane unit operations such as electrodialysis and electrodeionization. Custom unit operations have to be added in the form of user models implemented in various program or scripting languages e.g. Fortran or Python. Our objective is to develop a useful process simulator consisting of: (i) user-friendly graphic interface; (ii) set of solvers for differential algebraic equations; (iii) program for the calculation of dissociation equilibria including the database of chemical compounds and dissociation constants; and (iv) set of unit operations including electrodialysis, electrodeionization, reverse osmosis, different filtrations, simple mixers, splitters, tanks, heat exchangers, pumps, etc.

Model Reformulation and Efficient Simulation of Two Dimensional Electrochemical Thermal Models

Simulations of a two-dimensional battery model that includes temperature considerations will be performed. The two dimensional model will be based on the porous electrode model developed by Doyle et al.1 using concentrated solution theory to describe the internal behavior of a lithium-ion sandwich consisting of positive and negative porous electrodes, a separator, and current collectors. This model has proven versatile enough to incorporate additional advancements allowing for the development of similar models.2-7However, even with only a single linear dimension, solving this model is difficult and computationally expensive.Several reformulation and order reduction techniques have been presented in the past to improve the computational efficiency to allow for parameter estimation and optimization.8-13However, most of these models consider only a single macro dimension in addition to the pseudo dimension for the solid electrode particles. This is sufficient in many cases. However, when temperature is considered, a second dimension (parallel to the electrodes) becomes more important. Since the height and width of a battery are often orders of magnitude larger than the thickness of the cell, a significant temperature gradient can develop parallel to the electrodes. This can affect battery health and utilization as well as safety. Therefore, an effective method to solve the battery equations in two linear dimensions would allow for better analysis of battery operation when temperature is a concern.The model will be solved using techniques previously developed13to simulate a one dimensional coupled thermal electrochemical model. That technique will be expanded to incorporate 2D effects. Since a large number of equations arise even for the 1D case, reformulation and orthogonal collocation are essential to reduce the computational cost and maintain accuracy in the two spatial dimensions when nonlinear temperature dependent parameters are used.The reformulation will be performed using a coordinate transformation to rescale each region in the x-direction (perpendicular to the electrodes) to be solved in the interval from 0 to 1. This allows orthogonal collocation techniques to be used in the solution while utilizing simple trial functions. This has been used successfully in the past for 1D models. The reformulated system of equations will be solved using a fast differential algebraic equation solver, and the results will be compared with those obtained by a rigorous COMSOL simulation of the non-reformulated 2D thermal model. Acknowledgements The authors are thankful for the financial support from the United States Government, Advanced Research Projects Agency – Energy (ARPA-E), and the U.S. Department of Energy, under award number DE-AR0000275.

B-Spline Gaussian Collocation Software for Two-Dimensional Parabolic PDEs

AbstractIn this paper we describe newB-spline Gaussian collocation software for solving two-dimensional parabolic partial differential equations (PDEs) defined over a rectangular region. The numerical solution is represented as a bi-variate piecewise polynomial (using a tensor productB-spline basis) with time-dependent unknown coefficients. These coefficients are determined by imposing collocation conditions: the numerical solution is required to satisfy the PDE and boundary conditions at images of the Gauss points mapped onto certain subregions of the spatial domain. This leads to a large system of time-dependent differential algebraic equations (DAEs) which is solved using the DAE solver, DASPK. We provide numerical results in which we use the new software, called BACOL2D, to solve three test problems.

Time Efficient Solution of Phase Equilibria in Dynamic and Distributed Systems with Differential Algebraic Equation Solvers

Solution of phase equilibria with flash calculations is central in many processes. During the integrating of the conservation equations in these systems, flash calculations are traditionally solved in inner loops at each integration step. Some of these systems can be solved more efficiently using modern DAE solvers, where the differential equations and the algebraic equations describing the phase equilibrium are solved simultaneously. In this paper we present a framework, called the Thermodynamic Differential Algebraic Equation (TDAE) method, which handles most two-phase flash variants. The phase boundary is tracked, enabling a robust solution also with phase changes. The time consumption of the TDAE method has been compared to the traditional approach in several examples implemented in both Fortran and Matlab. In some cases, the TDAE method is more than 800 times faster, and in other cases the traditional methodology should be used. We will give insight into how and when DAE solvers can be used to speed up phase equilibrium calculations.

Simulation of Constrained Mechanical Systems—Part II: Explicit Numerical Integration

This paper presents an explicit to integrate differential algebraic equations (DAEs) method for simulations of constrained mechanical systems modeled with holonomic and nonholonomic constraints. The proposed DAE integrator is based on the equation of constrained motion developed in Part I of this work, which is discretized here using explicit ordinary differential equation schemes and applied to solve two nontrivial examples. The obtained results show that this integrator allows one to precisely solve constrained mechanical systems through long time periods. Unlike many other implicit DAE solvers which utilize iterative constraint correction, the presented DAE integrator is explicit, and it does not use any iteration. As a direct consequence, the present formulation is simple to implement, and is also well suited for real-time applications.

One-step 9-stage Hermite–Birkhoff–Taylor DAE solver of order 10

The HBT(10)9 method for ODEs is expanded into HBT(10)9DAE for solving nonstiff and moderately stiff systems of fully implicit differential algebraic equations (DAEs) of arbitrarily high fixed index. A scheme to generate first-order derivatives at off-step points is combined with Pryce scheme which generates high order derivatives at step points. The stepsize is controlled by a local error estimator. HBT(10)9DAE uses only the first four derivatives of y instead of the first 10 required by Taylor’s series method T10DAE of order 10. Dormand–Prince’s DP(8,7)13M for ODEs is extended to DP(8,7)DAE for DAEs. HBT(10)9DAE wins over DP(8,7)DAE on several test problems on the basis of CPU time as a function of relative error at the end of the interval of integration. An index-5 problem is equally well solved by HBT(10)9DAE and T10DAE. On this problem, the error in the solution by DP(8,7)DAE increases as time increases.

A Modular Modeling Approach to Simulate Interactively Multibody Systems With a Baumgarte/Uzawa Formulation

In this paper, a modular modeling approach of multibody systems adapted to interactive simulation is presented. This work is based on the study of the stability of two differential algebraic equation solvers. The first one is based on the acceleration-based augmented Lagrangian formulation and the second one on the Baumgarte formulation. We show that these two solvers give the same results and have to satisfy the same criteria to stabilize the algebraic constraint acceleration error. For a modular modeling approach, we propose to use the Baumgarte formulation and an iterative Uzawa algorithm to solve external constraint forces. This work is also the first step to validate the concept of two types of numerical components for object-oriented programming.

A quick and efficient method for consistent initialization of battery models

Secondary batteries are usually modeled as a system of coupled nonlinear partial differential equations. These models are typically solved by applying finite differences or other discretization techniques in the spatial directions and solving the resulting system of differential algebraic equations (DAEs) numerically in time. These DAEs are very difficult to solve even using popular DAE solvers. The complications arise partly due to the difficulty in obtaining consistent, or closely consistent, initial conditions for the DAEs. In this paper, a shooting method is proposed as an effective and rapid technique for the initialization of battery models. This method is built on a regionwise shooting approach with initial guess at one end of the electrode and physics based shooting criterion on the other end that can ultimately satisfy all the required conditions in a battery unit. Notably, the computation time required for the proposed method is only milliseconds in a FORTRAN environment for the case of initializing a standard physics based lithium-ion battery model. Also the initial values obtained are exact and can readily be fed into any DAE solver for achieving accurate solutions without solver failure. This rapid method will help in simulating batteries in hybrid environments in real-time (milliseconds).

Modeling reactive transport using exact solutions for first-order reaction networks

Operator splitting is often used for solving advection-dispersion-reaction (ADR) equations. Each operator can be solved separately using an algorithm appropriate to its mathematical behavior. Although a lot of research has been done in operator splitting for solving ADR equations, numerical approaches for the reaction operator are computationally expensive. To meet the convergence criteria of ODE (ordinary differential equation) or DAE (differential algebraic equation) solvers, a transport time step has to be subdivided into a large number of reaction time steps. Additional computation effort is also required to reduce the splitting error. In this paper, we develop exact solutions of various first-order reaction networks for the reaction operator and couple those solutions with numerical solutions of the transport operator. The reactions are treated as local phenomena and simulated using exact solutions that we develop, while advection and dispersion are treated as global processes and simulated numerically. The proposed method avoids the numerical error from the reaction operator and requires a single-step calculation to solve the reaction operator. Compared to conventional operator-splitting methods, the proposed method offers both computational efficiency and simulation accuracy.

Solving differential-algebraic equations by Taylor series (II): Computing the System Jacobian

The authors have developed a Taylor series method for solving numerically an initial-value problem differential-algebraic equation (DAE) that can be of high index, high order, nonlinear, and fully implicit, BIT, 45 (2005), pp. 561–592. Numerical results have shown that this method is efficient and very accurate. Moreover, it is particularly suitable for problems that are of too high an index for present DAE solvers. This paper develops an effective method for computing a DAE’s System Jacobian, which is needed in the structural analysis of the DAE and computation of Taylor coefficients. Our method involves preprocessing of the DAE and code generation employing automatic differentiation. Theory and algorithms for preprocessing and code generation are presented. An operator-overloading approach to computing the System Jacobian is also discussed.

Arbitrary order Krylov deferred correction methods for differential algebraic equations

In this paper, a new framework for the construction of accurate and efficient numerical methods for differential algebraic equation (DAE) initial value problems is presented. The methods are based on applying spectral deferred correction techniques as preconditioners to a Picard integral collocation formulation for the solution. The resulting preconditioned nonlinear system is solved using Newton–Krylov schemes such as the Newton-GMRES method. Least squares based orthogonal polynomial approximations are computed using Gaussian type quadratures, and spectral integration is used to avoid the numerically unstable differentiation operator. The resulting Krylov deferred correction (KDC) methods are of arbitrary order of accuracy and very stable. Preliminary results show that these new methods are very competitive with existing DAE solvers, particularly when high precision is desired.