Differential Elimination–Completion Algorithms for DAE and PDAE

Differential–algebraic equations (DAE) and partial differential–algebraic equations (PDAE) are systems of ordinary equations and PDAEs with constraints. They occur frequently in such applications as constrained multibody mechanics, spacecraft control, and incompressible fluid dynamics.
A DAE has differential index r if a minimum of r +1 differentiations of it are required before no new constraints are obtained. Although DAE of low differential index (0 or 1) are generally easier to solve numerically, higher index DAE present severe difficulties.
Reich et al. have presented a geometric theory and an algorithm for reducing DAE of high differential index to DAE of low differential index. Rabier and Rheinboldt also provided an existence and uniqueness theorem for DAE of low differential index. We show that for analytic autonomous first-order DAE, this algorithm is equivalent to the Cartan–Kuranishi algorithm for completing a system of differential equations to involutive form. The Cartan–Kuranishi algorithm has the advantage that it also applies to PDAE and delivers an existence and uniqueness theorem for systems in involutive form. We present an effective algorithm for computing the differential index of polynomially nonlinear DAE. A framework for the algorithmic analysis of perturbed systems of PDAE is introduced and related to the perturbation index of DAE. Examples including singular solutions, the Pendulum, and the Navier–Stokes equations are given. Discussion of computer algebra implementations is also provided.     

The Waveform Relaxation method for systems of differential/algebraic equations

This paper reports efforts towards establishing a parallel numerical algorithm known as Waveform Relaxation (WR) for simulating large systems of differential/algebraic equations. The WR algorithm was established as a relaxation based iterative method for the numerical integration of systems of ODEs over a finite time interval. In the WR approach, the system is broken into subsystems which are solved independently, with each subsystem using the previous iterate waveform as “guesses” about the behavior of the state variables in other subsystems. Waveforms are then exchanged between subsystems, and the subsystems are then resolved repeatedly with this improved information about the other subsystems until convergence is achieved.

In this paper, a WR algorithm is introduced for the simulation of generalized high-index DAE systems. As with ODEs, DAE systems often exhibit a multirate behavior in which the states vary as differing speeds. This can be exploited by partitioning the system into subsystems as in the WR for ODEs. One additional benefit of partitioning the DAE system into subsystems is that some of the resulting subsystems may be of lower index and, therefore, do not suffer from the numerical complications that high-index systems do. These lower index subsystems may therefore be solved by less specialized simulations. This increases the efficiency of the simulation since only a portion of the problem must be solved with specially tailored code. In addition, this paper established solvability requirements and convergence theorems for varying index DAE systems for WR simulation.     

Computed torque control of an under-actuated service robot platform modeled by natural coordinates

The paper investigates the motion planning of a suspended service robot platform equipped with ducted fan actuators. The platform consists of an RRT robot and a cable suspended swinging actuator that form a subsequent parallel kinematic chain and it is equipped with ducted fan actuators. In spite of the complementary ducted fan actuators, the system is under-actuated. The method of computed torques is applied to control the motion of the robot.The under-actuated systems have less control inputs than degrees of freedom. We assume that the investigated under-actuated system has desired outputs of the same number as inputs. In spite of the fact that the inverse dynamical calculation leads to the solution of a system of differential–algebraic equations (DAE), the desired control inputs can be determined uniquely by the method of computed torques.We use natural (Cartesian) coordinates to describe the configuration of the robot, while a set of algebraic equations represents the geometric constraints. In this modeling approach the mathematical model of the dynamical system itself is also a DAE.The paper discusses the inverse dynamics problem of the complex hybrid robotic system. The results include the desired actuator forces as well as the nominal coordinates corresponding to the desired motion of the carried payload. The method of computed torque control with a PD controller is applied to under-actuated systems described by natural coordinates, while the inverse dynamics is solved via the backward Euler discretization of the DAE system for which a general formalism is proposed. The results are compared with the closed form results obtained by simplified models of the system. Numerical simulation and experiments demonstrate the applicability of the presented concepts.     

Parallel optimal control with multiple shooting, constraints aggregation and adjoint methods

In this paper, constraint aggregation is combined with the adjoint and multiple shooting strategies for optimal control of differential algebraic equations (DAE) systems. The approach retains the inherent parallelism of the conventional multiple shooting method, while also being much more efficient for large scale problems. Constraint aggregation is employed to reduce the number of nonlinear continuity constraints in each multiple shooting interval, and its derivatives are computed by the adjoint DAE solver DASPKADJOINT together with ADIFOR and TAMC, the automatic differentiation software for forward and reverse mode, respectively. Numerical experiments demonstrate the effectiveness of the approach.     

Differentiability of Consistency Functions for DAE Systems

In the present paper, parametric initial-value problems for differential-algebraic (DAE) systems are investigated. It is known that the initial values of DAE systems must satisfy not only the original equations in the system but also the derivatives of these equations with respect to time. Whether or not this actually imposes additional constraints on the initial values depends on the particular problem.Often the initial values are not determined uniquely, so that the resulting degrees of freedom can be used to optimize a given performance index. For this purpose, the class of so-called consistency functions is defined. These functions map a set of parameters, which include also those undetermined initial values, to consistent initial values for the DAE system. Because of frequent gradient evaluations of the performance index and the constraints with respect to these system parameters needed by many optimization procedures, we state conditions such that the consistency functions represent differentiable functions with respect to these parameters. Several examples are provided to illustrate the verification of the theoretical assumptions and their differentiability properties.The authors would like to thank the referees for helpful suggestions and comments     

Generic pole assignability,structurally constrained controllers and unimodular completion

In this paper we assume dynamical systems are represented by linear differential-algebraic equations (DAEs) of order possibly higher than one. We consider a structured system of DAEs for both the to-be-controlled plant and the controller. We model the structure of the plant and the controller as an undirected and bipartite graph and formulate necessary and sufficient conditions on this graph for the structured controller to generically achieve arbitrary pole placement. A special case of this problem also gives new equivalent conditions for structural controllability of a plant. Use of results in matching theory, and in particular, ‘admissibility’ of edges and ‘elementary bipartite graphs’, make the problem and the solution very intuitive. Further, our approach requires standard graph algorithms to check the required conditions for generic arbitrary pole placement, thus helping in easily obtaining running time estimates for checking this. When applied to the state space case, for which the literature has running time estimates, our algorithm is faster for sparse state space systems and comparable for general state space systems.     

Some software aspects of DAE simulation & optimization software

Considerable numerical software has been written for simulation and optimization of dynamical systems. From the beginning of their development, differential algebraic equations (DAEs) have often been proposed as a way to make modeling easier. The modeler need only write down equations relating the variables in the model. However, much DAE software requires at least as much user numerical and mathematical expertise as explicit methods. An important aspect of our research has been working toward helping the idea of DAEs achieve its promise in modeling and simulation by both pushing the software to handle more general problems and to also allow for less user expertise. Some recent examples are presented where this research impacts on software and their underlying algorithms. Space necessitates we assume the reader has a rough idea of what a DAE is. The examples are implicit Scicos, and optimization of DAE models. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

A new approach for computing consistent initial values and Taylor coefficients for DAEs using projector-based constrained optimization

This article describes a new algorithm for the computation of consistent initial values for differential-algebraic equations (DAEs). The main idea is to formulate the task as a constrained optimization problem in which, for the differentiated components, the computed consistent values are as close as possible to user-given guesses. The generalization to compute Taylor coefficients results immediately, whereas the amount of consistent coefficients will depend on the size of the derivative array and the index of the DAE. The algorithm can be realized using automatic differentiation (AD) and sequential quadratic programming (SQP). The implementation in Python using AlgoPy and SLSQP has been tested successfully for several higher index problems.     

Minimal comparability completions of arbitrary graphs

A transitive orientation of an undirected graph is an assignment of directions to its edges so that these directed edges represent a transitive relation between the vertices of the graph. Not every graph has a transitive orientation, but every graph can be turned into a graph that has a transitive orientation, by adding edges. We study the problem of adding an inclusion minimal set of edges to an arbitrary graph so that the resulting graph is transitively orientable. We show that this problem can be solved in polynomial time, and we give a surprisingly simple algorithm for it. We use a vertex incremental approach in this algorithm, and we also give a more general result that describes graph classes Π for which Π completion of arbitrary graphs can be achieved through such a vertex incremental approach.     

Sub-exponential graph coloring algorithm for stencil-based Jacobian computations

Partial differential equations can be discretized using a regular Cartesian grid and a stencil-based method to approximate the partial derivatives. The computational effort for determining the associated Jacobian matrix can be reduced. This reduction can be modeled as a (grid) coloring problem. Currently, this problem is solved by using a heuristic approach for general graphs or by developing a formula for every single stencil. We introduce a sub-exponential algorithm using the Lipton–Tarjan separator in a divide-and-conquer approach to compute an optimal coloring. The practical relevance of the algorithm is evaluated when compared with an exponential algorithm and a greedy heuristic.     

Linear-time computation of optimal subgraphs of decomposable graphs

A general problem in computational graph theory is that of finding an optimal subgraph H of a given weighted graph G. The matching problem (which is easy) and the traveling salesman problem (which is not) are well-known examples of this general problem. In the literature one can also find a variety of ad hoc algorithms for solving certain special cases in linear time. We suggest a general approach for constructing linear-time algorithms in the case where the graph G is defined by certain rules of composition (as are trees, series-parallel graphs, and outerplanar graphs) and the desired subgraph H satisfies a property that is “regular” with respect to these rules of composition (as do matchings, dominating sets, and independent sets for all the classes just mentioned). This approach is applied to obtain a linear-time algorithm for computing the irredundance number of a tree, a problem for which no polynomial-time algorithm was previously known.     

Index minimization of differential-algebraic equations in hybrid analysis for circuit simulation

Modern modeling approaches for circuit analysis lead to differential-algebraic equations (DAEs). The index of a DAE is a measure of the degree of numerical difficulty. In general, the higher the index is, the more difficult it is to solve the DAE. The index of the DAE arising from the modified nodal analysis (MNA) is determined uniquely by the structure of the circuit. Instead, we consider a broader class of analysis method called the hybrid analysis. For linear time-invariant electric circuits, we devise a combinatorial algorithm for finding an optimal hybrid analysis in which the index of the DAE to be solved attains the minimum. The optimal hybrid analysis often results in a DAE with lower index than MNA.      

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. AMS subject classification (2000)  34A09, 65L80, 65L05, 41A58     

Mathematical modeling of complex mechanical systems

Mathematical modeling of mechanical systems based on multibody system models is a well tested approach. Generating the equations of motion for complex multibody systems with a large number of degrees of freedom is difficult with paper and pencil. For this reason methods for automatic equation generation have been developed. Most methods result in numerical equations of motion without explicit information about the parameters. In this paper a method is described resulting in symbolic equations of motion. The method allows also the determination of the constraint forces which are important for design purposes. The inverse problem of dynamics is also easily solved.     

On the general nonlinear self-adjoint spectral problem for differential-algebraic systems

We consider a general self-adjoint spectral problem, nonlinear with respect to the spectral parameter, for linear differential-algebraic systems of equations. Under some assumptions, we present a method for reducing such a problem to a general self-adjoint nonlinear spectral problem for a system of differential equations. In turn, this permits one to pass to a problem for a Hamiltonian system of ordinary differential equations. In particular, in this way, one can obtain a method for computing the number of eigenvalues of the original problem lying in a given range of the spectral parameter.     

The Index of an Infinite Dimensional Implicit System

The idea of the index of a differential algebraic equation (DAE) (or implicit differential equation) has played a fundamental role in both the analysis of DAEs and the development of numerical algorithms for DAEs. DAEs frequently arise as partial discretizations of partial differential equations (PDEs). In order to relate properties of the PDE to those of the resulting DAE it is necessary to have a concept of the index of a possibly constrained PDE. Using the finite dimensional theory as motivation, this paper will examine what one appropriate analogue is for infinite dimensional systems. A general definition approach will be given motivated by the desire to consider numerical methods. Specific examples illustrating several kinds of behavior will be considered in some detail. It is seen that our definition differs from purely algebraic definitions. Numerical solutions, and simulation difficulties, can be misinterpreted if this index information is missing.     

On a new approach to asymptotic stabilization problems

A numerical algorithm for solving the asymptotic stabilization problem by the initial data to a fixed hyperbolic point with a given rate is proposed and justified. The stabilization problem is reduced to projecting the resolving operator of the given evolution process on a strongly stable manifold. This approach makes it possible to apply the results to a wide class of semidynamical systems including those corresponding to partial differential equations. By way of example, a numerical solution of the problem of the asymptotic stabilization of unstable trajectories of the two-dimensional Chafee-Infante equation in a circular domain by the boundary conditions is given.     

DAE的Runge-Kutta方法在不可压NS方程求解中的应用

1．引言自然界中的流场通常是非定常复杂流场，要正确模拟和跟踪复杂流场的变化，计算格式的时间精度极为重要．对于常微分方程（＊＊q，一般采用＊K方法及线性多步法来提高格式的时间精度．前者是单步法，在计算过程中可以改变步长，可找到稳定性较好的高精度格式：近年来在发展到偏微分方程的数倩水解中也有很多应用．原始变量的INS方程（二维）为：其中u，v分别是x，y方向速度分量，r是压力，连续方程（1．幻可视为约束条件．从［1］，［2］可见，经空间差分化后（固定空间网格），它可看作带约束的微分方程组，即微分代数方程（DAE－…     

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.