Discretization schemes for constraint stabilization in nonlinear differential-algebraic systems

In this paper the problem of simulation of differential-algebraic systems is addressed. In modelling me- chanical systems the use of redundant coordinates and con- straints results in differential-algebraic equations, the integra- tion of which can lead to numerical instabilities, such as the so-called drift phenomenon. In [1] the authors have proposed a globally convergent conceptual continuous-time algorithm for the integration of constrained mechanical systems which ensures the existence of solutions and global attractivity of the solution manifold. The objective of this paper is to study the numerical implementation of the algorithm presented in [1]. In addition, the stability properties of the constrained system in the manifold are studied in both the continuous and discrete time cases. The proposed technique is illustrated by means of a simple example.


I. INTRODUCTION
The numerical analysis of differential-algebraic (DA) systems has gained great relevance in the last decades.The DA formulation is more flexible than that given by ordinary differential equations (ODEs) in modelling largescale systems, often composed by the interconnection of multiple differential subsystems and in which the relations among the subsystems are described by algebraic constraints.For instance, constrained mechanical systems are inherently simpler to model in the framework of differential-algebraic equations (DAEs).The interest raised in the simulation of mechanical systems, together with the increasing computational power of computers, have motivated the development of different techniques for the integration of DA systems, see for example [2], [3], [4], [5] and [6].
One of the most relevant problems in integrating DAEs is the drift phenomenon.Since the solution manifold is invariant but not attractive, see [7], any constraint violation due to numerical approximations may make the solution drift away from the manifold.Stabilization methods provide a possible solution to this problem within a control theory perspective.The basic idea of such methods is to stabilize the constraint manifold by feeding back the constraint violation.Among all the methods of this class, Baumgarte's stabilization has gained prominence because of its simplicity, see [8].However, this method does not solve all possible numerical instabilities as recently pointed out in [1].
In [1] the authors have proposed a globally convergent conceptual continuous-time algorithm for the integration of constrained mechanical systems based on the geometric interpretation of DA systems given in [9].This algorithm ensures global convergence of the constraint violation onto the manifold and existence of solutions.However, a discretization scheme is required for any numerical implementation of the algorithm.In addition, no results have been provided on the preservation of stability properties of the system in the manifold.
Based on the results presented in [9] and [1], the objective of this study is to investigate both the aforementioned problems.In Section II we introduce some preliminaries and recall results from [1] which are instrumental for the rest of the paper.The section is concluded with a new result regarding the stability properties for the "stabilized system" in continuous-time.In Section III we provide the main theoretical results of the paper: two discretization schemes for the algorithm proposed in [1] are presented and conditions for the preservation of the continuous-time stability properties given in [1] and in Section II are provided.In Section IV the method is validated via an example in which a simple nonlinear oscillator is integrated with one of the proposed schemes.Finally in Section V we report our conclusions.
Notation.We use standard notation.The set of real and natural numbers (including 0) are denoted, respectively, by R and N. The symbols R ≥0 and N ≥0 indicate the sets of non-negative real numbers and non-negative integer numbers, respectively.The superscripts and − represent the transposition operator and the transposition of the inverse operator, respectively.I represents the identity matrix.Given a matrix A, A ≺ 0 denotes that A is negative definite and symmetric.Given a vector x, the symbol ||x|| denotes its L 2 norm while the symbol x i denotes the i-th element of x.Given a map f : R n → R, we use the symbol ∂f ∂x to denote the row vector of partial derivatives of f with respect to the vector x ∈ R n .The symbol L f h denotes the Lie derivative of the function h along the vector field f , while

II. PRELIMINARIES
In this section we introduce some definitions and give some preliminary results which are instrumental for the remaining of the paper.To begin with we recall the relation between differential-algebraic systems and a class of ordinary differential equations presented in [1].At the end of the section a new result is given regarding the stability properties of the solutions of the system.
Consider a nonlinear differential-algebraic system of the form ẋ = f (x) + g(x)λ, where x(t) ∈ R n is the state vector, λ(t) ∈ R m is the algebraic variable, and f : R n → R n , g : R n → R n×m and h : R n → R m are smooth mappings.For simplicity we assume that m = 1 and that the origin is the only equilibrium point of the system.All statements are valid locally around the origin if not otherwise stated.We now recall the definition of index for a DA system.
[10] The differentiation index ν of the differential-algebraic system (1) is equal to the minimum number of times it is necessary to differentiate (1) to find an explicit expression of ẋ and λ as continuous functions of x and λ.
For the DA system (1) the following proposition holds.
[9] Consider system (1).Assume it has index (2) Assume that the Jacobian of Φ(x) is invertible and that φ i : R n → R, for i = r + 1, . . ., n, are such that L g φ i (x) = 0. Then system (1) can be written as where and q are smooth mappings and a(0, 0) = 0.Moreover, the solution manifold of system Usual index reduction techniques (see [11], [12] or [13]) would normally impose ξr = 0 and directly integrate the differential subsystem in (3).However, since the solution of (3) is it is clear that violation of the constraints would grow unbounded, if r > 1.This phenomenon is known as constraint "drift" and many stabilization methods have been proposed in the literature to deal with this issue, see for example [8], [14], [15] and [16].In addition, as pointed out in [1], stabilization methods which rely on linear feedback control theory do not guarantee existence of solutions for all t ≥ 0. Thus, in [1] a globally stable convergent continuous-time algorithm has been presented.The main properties of the algorithm therein are recalled in Proposition 2 and Theorem 1.Consider the following differential system where . ., ξ r ) , q 0 and q i , for i = 1, . . ., r, are smooth mappings and for i = 1, . . ., r, with > 0 and δ > 0. For such a system the following results hold.
[1] Consider system (5).Assume there exists a positive definite radially unbounded function W (µ) such that ∂W (µ) ∂µ q 0 (µ) ≤ γW (µ) + γ 0 , for some γ ∈ R and γ 0 ∈ R, and for some W > 0. Then there exists δ > 0 such that for all δ ≥ δ in (6) the following statements hold. 1) µ(t) and ξ(t) exist for all µ(0) ∈ R n−r and ξ(0 Theorem 1 ensures that the manifold M is globally attractive and that the solutions of the system exist for all t ≥ 0. However, it does not give any information on the stability properties of the constrained system.The first objective of this paper is to study the behaviour of the solutions in the manifold M. To this end we first recall a preliminary result.Let μ = q 0 (µ), be the constrained dynamics, also called zero dynamics, of system (3) (equivalently of system ( 5)) in the manifold M, i.e. for ξ(t) ≡ 0. Then the following result holds.
Proposition 3. [9] Suppose that the differential-algebraic system (1) has differentiation index ν ≤ n.Stability properties of the zero equilibrium of system (1) are equivalent to stability properties of the zero equilibrium of the zero dynamics of system (3).
The problem we want to address in the remainder of the section is to find out under which conditions the stability properties of the zero dynamics of system (3), and according to Proposition 3 of the DA system (1), are preserved by the stabilization method in equation (5).In particular, the following result holds.
1) µ(t) and ξ(t) exist for all µ(0) ∈ R n−r and ξ(0 Remark 1. Observe that assumption (7) with γ < 0 and γ 0 = 0 in Theorem 3 can be rewritten as where the function W (µ) evolves along the trajectory of system (9), i.e. the zero dynamics.According to Proposition 3 and Theorem 2, if the equilibrium point of system ( 1) is globally asymptotically stable then also the equilibrium point of the "stabilized" system ( 5) is globally asymptotically stable.

III. MAIN RESULT
In this section the problem of the discretization of system (5) is addressed.Since for a numerical implementation of any integration method a discretization scheme is required, we now study under which conditions the results obtained in [1] and in Section II of this paper for the continuous-time case are mantained when a discretization is performed.In the remainder of the paper we use the notation x := x(t) and x + := x(t + T ), ∀t ∈ N, where T ∈ R >0 is the sample period.
Consider the discrete-time system where . . ., q r ], ϕ T : R n−r → R n−r , ϕ T , q 0 and q i , for i = 1, . . ., r are smooth mappings, and K : R r × R n−r → R r×r .Recalling the definition of k i given in (6), let K(ξ, µ) be such that Observe that system (11) is a discretization scheme for system (5), provided the map ϕ T (q 0 ) is such that see for example [17] for conserving time-integration schemes for mechanical systems.Similarly to Theorem 1 for the continous-time system (5), for the discrete-time system (11) the following result holds.
Note that different discretization schemes provide different properties.For instance, consider the modified discrete-time system Observe that, since lim and since (13) holds, system ( 15) is a discretization of system (5).The following result holds.
Theorem 4 ensures that the manifold is globally attractive and that the solutions of the system exist for all t ∈ N ≥0 and for any T > 0. We consider now the stability properties of the zero equilibrium of systems (11) and (15).The following result provides a discrete-time equivalent of Theorem 2.
Remark 2. Observe that assumption ( 16) can be rewritten as where W (µ) is a positive definite function which evolves along the trajectories of the system i.e. the constrained system.Moreover, note that for T → 0 equation ( 16) becomes lim Hence, according to Remark 1, condition (16) is the discretetime counterpart of (7).Similarly, condition (17) is the discrete-time counterpart of (8).This observation implies that for T small enough, ( 7) and ( 8) imply ( 16) and (17).
Observe that while the discretization scheme in (15) ensures global attractivity of the constraint manifold, for the discretization scheme (11) attractivity holds only locally, see Theorems 3 and 4.However, this comes at the cost of a smaller sample time required to preserve the stability properties of the system, see the assumptions of Theorem 5.The results so far presented, together with the algorithm presented in [1], outline a procedure to simulate DA systems, as shown in the following example.

IV. EXAMPLE
Consider the nonlinear oscillator depicted in Fig. 1 with mass m and spring constant k.Let (X, Y ) and ( Ẋ, Ẏ ) be the position and the velocity coordinates of the mass m, respectively.Assume that the only external forces acting on the system are the gravitational force (with g a the gravitational acceleration), the constraint force λ of the floor, the friction force F d = −d Ẋ and the nonlinear spring force F k = −kX 3 .Such a system can be described by a set of differential-algebraic equations in the state space form (1). m k X Y Fig. 1: A nonlinear oscillator.
To this end, let and set calculated according to Proposition 1.The system in the new coordinates is described by the equations For simulation purposes we consider the "stabilized" and discretized system where and where has been selected according to [17] to preserve the total mechanical energy.
Note that the assumptions of Theorem 5 are satisfied for T = 0.001, with and γ = 0.00001.Simulation results are shown in Figs. 2, 3, 4 and 5.As expected from Theorems 4 and 5, we observe in Fig. 2 that the constraint violation converges to zero, while from Fig. 3 we observe that the origin is an asymptotically stable equilibrium for the constrained system.Fig. 4 displays the trajectory of the system in the X-Y plane, while Fig. 5 shows the total mechanical energy which is dissipated because of the friction force acting on the system.
Consider now the case where d = 0.The results of the simulation for the same set of parameters are shown in Figs. 6, 7, 8 and 9.As expected from Theorem 4, we observe from Fig. 7 that the constraint violation converges to zero.In addition, observe from Fig. 9 that the mechanical energy is preserved at steady state, i.e. once solutions reach the manifold M.This observation suggests that the method proposed can be used also in the case in which the equilibrium is marginal stable, even though this has not been proved.

V. CONCLUSIONS
We have studied the problem of integration of DA systems.In Section II we have extended the theory developed in [1] by studying the stability properties of the "stabilized" system in continuous-time.We have concluded that if the equilibrium  of the constrained system is globally asymptotically stable then the stabilization method preserves this property.In Section III we have extended the results to discrete-time.We have concluded that the solution manifold remains globally attractive for one of the discretization scheme proposed for any sampling time.Moreover asymptotic stability of the equilibrium is preserved (locally).Finally, in Section IV we have demonstrated the technique with an example and discussed the results.

Fig. 2 :
Fig. 2: Time history of the coordinate Y for d = 0.03.

Fig. 5 :
Fig. 5: Time history of the total mechanical energy for d = 0.03.

Fig. 8 :Fig. 9 :
Fig. 8: Path of the mass m on the X-Y plane for d = 0.