A Globally Stable Algorithm for the Integration of High-Index Differential-Algebraic Systems

The problem of constraint stabilization and numerical integration for differential-algebraic systems is addressed using Lyapunov theory. It is observed that the application of stabilization methods which rely on a linear feedback mechanism to nonlinear systems may result in trajectories with finite escape time. To overcome this problem, we propose a method based on a nonlinear stabilization mechanism that guarantees the global existence and convergence of the solutions. Discretization schemes, which preserve the properties of the method, are also presented. The results are illustrated by means of the numerical integration of a slider-crank mechanism.


I. INTRODUCTION
D IFFERENTIAL-ALGEBRAIC systems (also known as DAE systems, descriptor systems or singular systems) are composite systems of differential and algebraic equations. The mathematical modeling of a multitude of engineering systems, such as mechanical systems [1], chemical processes [2] and electrical networks [3], is straightforward in the differentialalgebraic formulation: differential equations can be used to describe each dynamical subsystem independently whereas algebraic equations allow describing phenomena such as environmental and topological constraints, mass and flow conservation, and thermodynamic relations.
Due to the large field of application in which DAE systems are involved, the stability analysis and control of such systems have gained great attention from the research community. In the linear case, DAE systems have been studied in [4], in which several concepts from the theory of dynamical systems have been extended to this class of systems. Then, [5] and [6] have provided a Lyapunov-based analysis for linear DAE systems. In [7], the stability analysis of DAE systems has been performed by means of linear matrix inequalities, yielding necessary and sufficient conditions for the linear case and sufficient conditions for the nonlinear case. One of the first studies in the control of nonlinear DAE systems has been presented in [8], in which global stabilization and tracking problems for constrained mechanical systems have been solved by means of force feedback, while local results have been obtained by means of a linear controller. Many existing control results for DAE systems rely on the derivation of appropriate state-space realizations of the constrained system. For example, McClamroch [9] has studied the feedback stabilization problem for nonlinear DAE systems by means of a change of coordinates that transforms the DAE system in the so-called normal form, and local results have been obtained by the linearization of the transformed system. Recently, the zero dynamics form for DAE systems has been introduced in [10], see also [11].
Although the modeling process is simplified by the use of differential-algebraic equations, numerous complications arise in their numerical integration. 1 An important concept that provides a measure of these difficulties is the concept of index. Loosely speaking, the index indicates the number of time differentiations required to reduce a DAE system to a system of ordinary differential equations (ODEs), see [12] for a precise definition. While numerical methods for the integration of index-1 DAE systems are well established, see, e.g., [13]- [15], for high-index systems direct integration methods are yet to be developed. Good reviews on some of the approaches to the integration of high-index DAE systems are given in [16] and [17].
Coordinates reduction techniques are used to reduce DAE systems to a minimum set of ODEs, which can then be directly integrated by means of classical methods. Such techniques include Maggi's formulation (see [18]- [20]), the Index-1 formulation (see [21]- [24]), the null-space formulation [22], and Udwadia and Kalaba's formulation, see [25] and [26]. However, when numerically integrating the resulting ODEs, the initial conditions must be consistent with the original set of algebraic equations and with a finite number of their time derivatives. In addition, round-off errors introduced by the numerical solver may yield a constraint violation. Since, for DAE systems, the solution manifold 2 is invariant but not attractive, see [27], any constraint violation due to inconsistent initial conditions or numerical approximations may make the solution drift away from the manifold: this is called drift phenomenon or constraint drift. Hence, several methods have been proposed to guarantee that the solution remains in (or close to) the solution manifold. Mattsson and Söderlind [28] have proposed an index reduction technique in combination with the Pantelides's algorithm [29] to find consistent initial values. The result of the reduction process is an augmented, but consistent, index-1 DAE system, which can be directly integrated by means of index-1 methods. Since the required differentiations are carried out analytically rather than numerically, the technique has an accuracy comparable to that of solving a state-space formulation of the problem. However, the reduction comes at the cost of a higher analytical computation due to the higher number of equations and variables involved. Projection methods aim to project the solution onto the solution manifold by means of an iterative process, see [30] and [31]. As a drawback, the adaptive and iterative character of these methods makes them unsuited for real-time simulations. To guarantee that the projection step is successfully completed in a priori fixed sampling interval, Arnold et al. [32] have proposed and discussed a noniterative projection strategy yielding a bound for the constraint error which is independent of time. Other approaches, known as stabilization methods, are based on ideas and tools borrowed from control theory. They aim to avoid drift by adding extra terms that, although vanishing on the solution manifold, have the effect of making this manifold asymptotically attractive. The most widely used method of this kind is the Baumgarte's method [33], which is based upon the principle of feeding back a linear combination of all the constraint errors. The choice of the feedback parameters depends on several factors, such as the model and the numerical solver used. A direct comparison of the aforementioned methods, with a list of advantages and disadvantages and a case study, can be found in [34] and [35]. Some attempts to improve Baumgarte's method have been presented in [36], in which an adaptive computation of the feedback parameters has been proposed, or in [37], in which the integral of the constraint violation has been added to the feedback law. A generalization of the Baumgarte's method for nonlinear systems is presented in [38], in which stability properties of the solution manifold along with discretization schemes that preserve such properties have been derived using classical Lyapunov theory. An interpretation of the Baumgarte's method as an output nullifying feedback control has been provided, in the framework of the Lie algebraic control theory, in [39]. Recently, Chang et al. [40] have proposed a nonlinear stabilization method for dynamical systems on manifolds that preserves first integrals of the system, such as the kinetic energy or the module of the angular momentum for mechanical systems, for any ordinary numerical integrator. However, stabilization methods that rely on a linear feedback mechanism may lead to "closed-loop systems" with solutions having finite escape time. To the best of our knowledge, none of the existing methods addresses this issue. 2 See Remark 1 for a precise definition.
It is important to observe that stabilization methods add extra terms to the equations that affect the solution's trajectory off the solution manifold. Depending on the stabilization method and on the constraint error, stability properties of the equilibrium points of the DAE system may be changed, resulting in possibly unexpected nonphysical behavior. Moreover, additional instabilities may arise whenever a time discretization of the system is performed, i.e., whenever the integration method is replaced by a numerical algorithm. With Baumgarte's method, for instance, the trajectory of the discretized system may drift away from the solution manifold even though the manifold is asymptotically stable for the underlying ODE system in the continuous-time domain. Thus, it is of interest to ensure that stability properties of the solution manifold and of the equilibrium points are preserved when discretization schemes are applied. It is worth stressing, therefore, that Baumgarte's method, and all its variations, rely on notions (and tools) from control theory, hence, improvements can be obtained using ideas and tools borrowed from control theory.
In this paper, the problems of constraint stabilization and numerical integration for DAE systems with high index are addressed by means of concepts borrowed from classical geometric control theory, such as the notions of relative degree and zero dynamics (see [41] for an introduction to these topics), and by means of classical Lyapunov theory. Popular stabilization methods such as Baumgarte's method, or the one in [38], are recast in the proposed framework and are shown to potentially lead to systems exhibiting solutions with finite escape time. On this basis, we propose a nonlinear stabilization method that ensures the existence of the solution for all times. In addition, in the case in which the underlying DAE system has an equilibrium point which is asymptotically stable, we show that the proposed method preserves this property in a neighborhood of the solution manifold. To avoid instabilities due to time discretization, we present two discrete-time schemes and conditions to select the sampling time such that the properties of the method are preserved. The results are illustrated via a step-by-step algorithm for the simulation of DAE systems and validated with an example. Preliminary results have been published in [42]- [44]. The additional contributions of this paper are as follows: the results are presented in a more organized way with formal proofs; two popular stabilization methods have been showed to potentially lead to solutions that may exhibit finite escape time; the simulation algorithm is revisited and complemented by discretization schemes; a new worked-out example with complete calculations and comparisons with two popular stabilization methods is provided.
The rest of this paper is organized as follows. In Section II, we recall some general concepts of the theory of DAE systems and we give some preliminary results. In Section III, we study the problem of constraint stabilization and main results are derived in the continuous-time domain. In Section IV, we present two discretization schemes that preserve the properties highlighted in Section III. Moreover, we summarize the proposed results in an algorithm for the simulation of nonlinear DAE systems. In Section V, the method is validated with an example (a slidercrank mechanism) and simulation results are discussed. Finally, Section VI reports our conclusions.
Notation: We use standard notation. The symbols R ≥0 and N ≥0 indicate the sets of nonnegative real and integer numbers, respectively. The superscript represents the transposition operator. Given a vector x, the symbol x denotes its L 2 norm. Given a map f : R n → R, we use equivalently the symbols ∇ x f and ∂ 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 , provided all partial derivatives exist. The symbol y (k ) (t) indicates the kth time derivative of the function y at time t, provided it exists. Given a manifold Z, the symbol f | Z indicates the restriction of f to Z. All functions and mappings are assumed to be differentiable as many times as needed.

II. PRELIMINARIES
The objective of this section is to introduce some preliminaries on DAE systems and define a formalism to study the stability properties of such systems. The framework presented is useful to highlight some limitations of a class of stabilization methods for the integration of DAE systems based on linear control theory and represents the basis to develop nonlinear discretization schemes that remove the aforementioned limitations.
Consider the DAE system in the semiexplicit 3 forṁ 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 vector fields and mappings. For simplicity of notation, we assume, in the remainder of this paper, that m = 1, if not stated otherwise. We also assume that the DAE system (1) is solvable in the sense of [46, Definition 2.2.1]. Note that the form in (1) contains a general class of DAE systems that includes, for example, mechanical systems with holonomic and nonholonomic constraints [47], and chemical systems described by lumped parameters models or via spatial discretization of the underlying partial differential equations [2]. A concept that plays an important role in the classification and behavior of DAE systems is the one of differentiation index. Definition 1: [46] The differentiation index ν of the DAE system (1) at a point (x • , λ • ) is equal to the minimum number of times it is necessary to differentiate all or part of the equations in (1) with respect to t to determineẋ andλ as continuous functions of x and λ in a neighborhood of (x • , λ • ).
Several notions of index have been defined for DAE systems, see, e.g., [12] for a detailed explanation. Here, we recall the fact that for the DAE system in the semiexplicit form (1), the perturbation index, the tractability index and the geometric index are numerically equal to the differentiation index, although these are defined under different technical (smoothness) assumptions. Finally, the strangeness index, when defined, is equal to ν − 1.
In the remaining of this paper, the term index refers to the differentiation index. Observe that the description (1) represents a class of DAE systems with index ν > 1. This is without loss of generality since index-1 DAE systems are structurally similar to ordinary differential systems and numerical methods for their integration are well established (see [13]- [15]), thus they are not considered in this study.
To study the problem of stabilization of manifolds, it is helpful to consider the DAE system (1) as an input-affine nonlinear system in which the "input" λ(t) and the initial conditions x(t 0 ) = x 0 are consistent with the conditions that the "output" h(x) and its time derivative up to the order ν − 1, i.e., y (i) (t) = h (i) (x(t)), i = 0, . . . , ν − 1, are identically zero for all times. Stabilization of manifolds is a well-established area of research in the framework of differential geometry, see, e.g., [41]. In particular, attempts to study DAE systems in this framework have been made in [9], which has studied the feedback stabilization problem for nonlinear DAE systems by means of a change of coordinates that transforms the DAE system in the so-called normal form. With reference to [41], it is convenient to introduce the concept of relative degree, which is here adapted in the context in which the algebraic equation in (1) is considered as an output.
Definition 2: The following preliminary result shows that there exists a relation between the concepts of index and relative degree.
Lemma 1: Let ν and r be the index and the relative degree of the DAE system (1), respectively, and assume these are well defined at (x • , λ • ) and x • , respectively. Then, ν = r + 1.
Proof: Assume that the system (1) has relative degree r at x • . Differentiating with respect to time the output y = h(x) yields y (k ) (t) = L k f h(x(t)) for all k < r and for all t for which x(t) is defined and close to x • , and Hence, the relative degree r is equal to the number of times the output y has to be differentiated to have the input λ explicitly appearing in (2). Since y (i) (t) = 0, for i ≥ 0 and for all t for which x(t) is defined and close to x • , then, from (2), we obtain and Differentiating once more (3), we obtain the expression ofλ as function of x and λ, proving the claim.
The existence of a well-defined relative degree for the DAE system (1) implies the existence of a local diffeomorphism that induces on the system (1) a special structure, the so-called normal form [41], as described in the following statement. Note also that we remove, for simplicity, reference to the points (x • , λ • ) or x • . All statements are, however, locally defined around a point Proposition 1: Consider the system (1) and assume it has . . , φ n (4) in which the mappings φ i : R n → R, for i = r + 1, . . . , n are such that the Jacobian of Φ is invertible and L g φ i = 0. Then, the system (1) can be written in the coordinates z = (ξ , η ) = Φ(x), as where ξ i (t) ∈ R, for i = 1, . . . , r, ξ = (ξ 1 , . . . , ξ r ) , η(t) ∈ R n −r and a, b, and q are smooth mappings such that a(0, 0) = 0. Proof: Observe that, by definition, the algebraic subsystem of (1) can be written in the coordinates z = Φ(x) as z 1 = ξ 1 = 0. On the other hand, the differential subsystem of (1) can be written as the one in (5), see [41,Proposition 4.1.3]. Therefore, the claim follows by appending the algebraic equation to the differential subsystem.
Remark 1: The solution manifold of the system (5) is and there exists a unique smooth mapping such thatξ r (t) = 0 for all t, which coincides with (3) expressed in the coordinates z = (ξ , η ) defined in Proposition 1. Note also that the η subsystem of (5) restricted to the solution manifold, i.e.,η coincides with the zero dynamics of the auxiliary systemẋ = f (x) + g(x)λ with input λ and output y = h(x). We, therefore, refer to (8) as the zero dynamics of the DAE system (1). A direct consequence of the reformulation in Proposition 1 is that the stability properties of the equilibrium points of (1) can be studied in the special form (5), as shown in the following theorem.
Theorem 1: Suppose that the DAE system (1) has index 1 < ν ≤ n. Stability properties, see [48,Definition 2.3], of the equilibrium points of system (1) are equivalent to the stability properties of the zero equilibrium of the zero dynamics (8).
Proof: From Proposition 1, there exists a change of coordinates that transforms system (1) in the normal form (5). In particular, it is always possible to choose the coordinates Φ(x) = (ξ 1 , . . . , ξ r , η ) such that Φ(x • ) = (ξ • 1 . . . ξ • r (η • ) ) = 0. Note now that for any (ξ(0) , η(0) ) ∈ M, where M is given in (6), ξ i (t) = 0 for all t ≥ 0, for i = 1, . . . , r, while the subsystem of (5) in the variable η reduces to the zero dynamics (8). As a consequence, the stability properties of the equilibrium point x • are equivalent to the stability properties of the equilibrium point η • , proving the claim. (8) has dimension equal to (from Lemma 1) n − r = n − ν + 1, where ν − 1 is the actual 4 number of constraints acting on the system. Thus, (8) is a local representation of the system (1) with the minimum number of independent coordinates.

Remark 2: System
Remark 3: The case ν > n implies the existence of at least n constraints acting on the system. This means that the state belongs (generically) to a zero-dimensional manifold, or the equations have no solution.
Remark 4: In the general case in which m > 1, the existence of a well-defined (vector) relative degree cannot be inferred from the well definiteness of the differential index. However, if x • is a regular point of the zero dynamics algorithm [41], the application of such algorithm leads to the calculation of a maximal zeroing submanifold onto which a zero dynamics can be defined. Moreover, the steps of the zero dynamics algorithm are helpful in defining a local diffeomorphism that transforms the DAE system (1) in the so-called global normal form, which reduces to the form (5) in the case m = 1. Although, in general, the algebraic variable appears explicitly in the η subsystem of (5), there always exists a smooth 5 λ = λ * (ξ, η), solution of the system of equationsξ = 0, which renders the η subsystem function of ξ and η only, namelyη = q(ξ, η). If, in addition, the (vector) relative degree is well defined, then the mapping λ * (ξ, η) is unique.

III. CONSTRAINT STABILIZATION
The objective of this section is to address the problem of integration of DAE systems. Typical problems arising when DAE systems are numerically integrated include inconsistent initial conditions, roundoff errors, and constraint drift. Popular methods to avoid these issues, such as Baumgarte's method [33], rely on ideas borrowed from linear feedback theory. We show that the application of such methods may lead to systems with solutions having finite escape time. Nonlinear "stabilizer" are then designed to cope with this problem. Moreover, in the case in which an equilibrium point of the DAE system is asymptotically stable, the proposed method preserves this property in a neighborhood of the solution manifold. All the aforementioned problems and objectives are recast within the formalism presented in Section II and are studied, initially, in the continuous-time domain. Consider the DAE system (5) with index ν = r + 1, r ∈ R > 0 . Usual index reduction techniques (see [21], [22] or [25]) would normally imposeξ r = 0 and directly integrate the differential subsystem in (5). However, the solution of (5), whenξ r = 0 and the equation for all t ≥t: it is clear that the violation of the constraints would grow unbounded if r > 1. This describes the drift phenomenon or constraint drift and reflects the property of 4 This includes the "hidden" constraints that can be revealed by repeated time differentiations of the original set of algebraic equations. 5 In the general case in which m > 1, we use the symbol ξ to indicate the column vector col(ξ 1 , . . . , ξ m ), in which ξ j = col(ξ j 1 , . . . , ξ j r j ), ξ j i ∈ R for i = 1, . . . , r j and j = 1, . . . , m, r j = ν j − 1 and ν j is the index of the jth component of the vector h(x). However, with the aim of maintaining the notation simple, we drop the superscript j for the case m = 1.
the manifold M of being invariant but not attractive, see [27]. Stabilization methods aim to cope with the constraint drift by feeding back the constraint violation. In general, stabilization methods based on linear feedback directly integrate differential systems described by equations of the forṁ ξ = Aξ,η = q(ξ, η) (9) in place of the system (5), where A ∈ R r ×r is Hurwitz. A popular stabilization technique that is consistent with the aforementioned structure is Baumgarte's [33], which replaces the equationξ r = 0 with the equatioṅ in which the coefficients α j , for j = 1, . . . , r, are such that the roots of the polynomial σ(τ ) = τ r + r −1 j =1 α j τ j , in the indeterminate τ have all negative real part. Despite the conceptual simplicity of such a method, this may be inadequate for nonlinear DAE systems. To see this point, consider the following example.
Example 1: Consider the DAE system with a linear stabilizeŕ a-la Baumgarte described by the equations  [38]. 6 This method can be described, in some cases, by equations of the form (9). To see this, consider the reduced index-1 DAE systeṁ where the mapping f has been obtained from the differential subsystem in (1) with the algebraic variable eliminated, i.e.
with λ * (x) given in (3) and Ascher et al. [38] consider the family of stabilization methods given byẋ where γ > 0 and D(x) chosen such that C(x)D(x) is nonsingular, e.g., is full rank). Note that system (15) can be expressed in the coordinates z = (ξ, η) = Φ(x) asż . (17) Let is such that (7) holds, then system (17) in the coordinates z = (ξ, η) becomesξ where q : R n → R n −r and Note now that in the special case in which the mapping h(x) is linear, which is typical for instance in mechanical systems with linear holonomic constraints (see the example in Section VI), it follows from (16) that the matrix F is constant. Moreover, the matrixJ is constant by definition, thus, as a consequence, the system (18) can be described by equations of the form (9), with A = R − γJF and q(ξ, η) = q(ξ, η) − γ J(Φ −1 (ξ, η))F ξ, proving that the solutions arising from the method proposed in [38] may result in trajectories that may have finite escape time. Remark 6: If ξ(t) = 0 for some t, the trajectory of the system (9) is different from the trajectory of the DAE system (5). Indeed, under the stated smoothness assumptions, the η subsystem of (9) can be rewritten as 7 where q 0 (η) = q(0, η) is the mapping in (8) and q i : R r × R n −r → R n −r , for i = 1, . . . , r are smooth mappings. From (20), it is clear that if ξ(t) = 0 for somet ≥ 0, the solution η(t), for t ≥t, is affected by an external disturbance that depends on the mappings q i for i = 1, . . . , r and on the matrix A in (9).
Observe that not only such disturbances may lead to the finite escape time phenomenon in the worst-case scenario, but they may also "destabilize" trajectories that are supposed to converge to an equilibrium. Based on this observation, for the special case in which the DAE system (1) has an asymptotically stable equilibrium point, an additional objective of the stabilization method is to preserve the convergence of the trajectories initialized in a neighborhood of the solution manifold.
Proof: On the solution manifold M of system (5), the equalities hold, see (7). Note now that with q 0 (η) = q(0, η). Then, by replacing (23) and (24) in (5) and (23) in (21) yields the claim. Since the modified system (21) produces the same trajectory as (5) on the manifold M, we can legitimately modify the dynamics outside M to render M attractive. To this end, we state the main result of this section.
Theorem 2: Consider system (21). Assume there exists a positive definite and radially unbounded function W : for some γ ∈ R and γ 0 ∈ R, and for some W > 0. Then, there existδ > 0 and¯ > 0 such that for all δ ≥δ and ≥¯ , the following statements hold: 1) η(t) and ξ(t) exist for all η(0) ∈ R n −r , ξ(0) ∈ R r and t ≥ 0; 2) lim t→∞ ξ(t) = 0, for all ξ(0) ∈ R r . Proof: To prove claim 1) consider the positive definite function Its time derivative along the trajectories of the system iṡ for any δ > 0, and by replacing the expression of k i (ξ, η) given in (22) Note now that under assumption (26), it is always possible to chooseδ such thatδ where β > 0. Without loss of generality, assume that γ > 0 in (25). Hence, using assumption (25) in (28), for any δ ≥δ and ≥ 0, yieldṡ From the last inequality, it follows that V (ξ, η) ≤ V (ξ(0), η(0))e (β +γ )t + γ 0 t, from which it is clear that ξ(t) and η(t) exist for all t ≥ 0, proving claim 1). To prove claim 2), consider the Lyapunov function candidate and the ξ subsystem of (21). The time derivative of the Lyapunov function along the trajectories of the ξ subsystem iṡ for all ξ i = 0, thus proving claim 2). Note that the ξ subsystem is well defined by claim 1). The global results of Theorem 2 are indeed conceptually appealing but require the explicit knowledge of a function W that satisfies assumptions (25) and (26) globally. This difficulty can be partially overcome if instead of achieving global results, one is interested in local results. For the case in which γ < 0 and γ 0 = 0 in (25), the following result shows that the statements of Theorem 2 may be recast in a local fashion. Moreover, an additional insight on the stability properties of the η subsystem is given.
Remark 8: For the general case in which m > 1 and provided that the DAE system (1) has a well-defined (vector) relative degree, we consider the systeṁ in place of (21), where q j i , for i = 1, . . . , r j and j = 1, . . . , m, are smooth mappings and Observe also that the claims of Theorems 2 and 3 remain valid. The proof follows the same steps as the case m = 1, in which the functions (27) and (31) are replaced by Note that if stability properties of the DAE system (1) are known in advance, these can be exploited to verify the condition (25), as stated in the following corollary.
Corollary 1: If the equilibrium point x • is a globally asymptotically stable equilibrium point for the system (1), then condition (25) of Theorem 2 is satisfied with strict inequality for γ 0 = 0 and any γ ≥ 0.

a(ξ,η )
. Observe also that the previous expression coincides with (7) on the solution manifold M.
The assumptions of Theorems 2 and 3 and the stabilizer (22) require the explicit knowledge of the mappings q 0 (η) and q i (ξ, η), for i = 1, . . . , r, which in turn require the calculation of the inverse mapping x = Φ −1 (z). These steps can be avoided as indicated in the following statements.
Lemma 3: The claims of Theorems 2 and 3 remain valid if the condition (25) is replaced by for some p i : R r × R n −r → R n −r , and the gains (22) are replaced by Proof: Consider the positive definite function (27). Its time derivative along the trajectory of the system (21) can be written asV The rest of the proof follows the same steps of the proof of Theorems 2 and 3, in which q 0 (η) and q i (ξ, η), for i = 1, . . . , r, are replaced by q(ξ, η) and p i (ξ, η), respectively.
The formulation in (21) still requires the inversion x = Φ −1 (ξ, η). In the next result, we show that the results in Theorems 2 and 3 can be derived without inverting any coordinate transformation. To this end, let where ) and ψ(x) = (φ r +1 , . . . , φ n ) , in which the mappings φ i : R n → R, for i = r + 1, . . . , n are such that the Jacobian of Φ is invertible and L g φ i = 0. Consider the nonlinear systeṁ with f (x) given in (13) and in which δ > 0, > 0 and s i : R n → R n , for i = 1, . . . , r to be defined and R given in (19). Then, the following result holds. Theorem 4: Consider the nonlinear system (38). Assume there exist mappings ψ and s i , for i = 1, . . . , r, and a positive definite and radially unbounded function W (ψ(x)) such that the Jacobian of Φ given in (37) is invertible and for some γ ∈ R and γ 0 ∈ R, and sup x∈R n for some W > 0. Then, there existδ > 0 and¯ > 0 such that, for all δ ≥δ and ≥¯ , the following statements hold. 1) x(t) exists for all x(0) ∈ R n and for all t ≥ 0; 2) lim t→∞ h (i) (x(t)) = 0, for i = 0, 1, . . . , r − 1.
In addition, if γ < 0 and γ 0 = 0, then there existδ > 0 and > 0 such that for all δ ≥δ and ≥¯ in (40), the following statement holds. 3 Proof: The proof of the claims is based on the fact that the system (38) is described in the coordinates z = Φ(x) by (21) and (36), for which all the statements hold. Differentiating z = Φ(x) with respect to time, we obtainż = ∇ x Φ(x)ẋ, and by replacing (38) in the previous equation, yieldṡ Note now that by applying the coordinate transformation x = Φ −1 (z) to the first term of (43) yields the differential subsystem in (5), in whichξ r = 0, i.e., . Replacing (39) in the second term of (43) and by applying the inverse coordinate transformation Defining 1 (ξ, η)), then the system (44) is equal to the system (21) with k i (ξ, η) given in (36). By assumption, the positive definite function ≤ γW (ψ(x)) + γ 0 which can be rewritten in the coordinates (ξ, η) as in (35), while (42) can be rewritten as (32). Therefore, by Lemma 3, the claims follow. Remark 11: With arguments similar to those in Theorem 3, the claims of Theorem 4 can be restated in a local fashion.

IV. DISCRETIZATION SCHEMES
Since for a numerical implementation of any integration method, a discretization scheme is required, in this section the problem of the discretization of the system (21) is addressed. In fact, note that the solution trajectory of the difference equations resulting from the discretization of (21) may diverge even though the trajectories of the underlying ODEs converge to the solution manifold. Hence, we now propose two integration schemes for the system (21), which preserve the properties discussed in Section III. In the remainder of this paper, we use the notation x := x(t) and x + := x(t + T ), where t = kT ∀ k ∈ N and T ∈ R > 0 is the sampling period.
Consider the discrete-time system . . , q r ], ϕ T : R n −r → R n −r , q 0 and q i , for i = 1, . . . , r defined in (20), and K : R r × R n −r → R r ×r . Recalling the definition of k i given in (22), let K(ξ, η) be such that Observe that, for T → 0, the system (45) is equal to the system (21), provided the map ϕ T (q 0 ) is such that For instance, if ϕ T (q 0 ) = q 0 (48) then the system (45) and (48) represents the forward Euler discretization of the system (21). Similarly to Theorem 2 for the continuous-time system (21), for the discrete-time system (45) the following result holds. Theorem 5: Consider system (45) and let U ⊆ R n be a compact set that contains the origin. For all T such that there exists a neighborhood of the origin U ⊆ U such that the following claims hold. 1) ξ(t) and η(t) exist for all ξ(0) ∈ R r and η(0) ∈ R n −r ; 2) lim t→+∞ ξ(t) = 0, for all (ξ(0), η(0)) ∈ U . Proof: The proof of claim 1) trivially follows from the fact that the right-hand side of (45) is continuous. To prove claim 2) consider the positive definite and radially unbounded function V (ξ) = ξ ξ. It follows that with (51) for i = 1, . . . , r. Note now that the selection of T as in (49) implies that −2 < T k i (ξ, η) < 0, and since k i (ξ, η) < 0, also d i < 0 for i = 1, . . . , r. Hence, ΔV = T r i=1 d i ξ 2 i < 0, for all ξ = 0, proving claim 2).
Note that different discretization schemes yield different properties. For instance, consider the modified discrete-time system Observe that, since lim T →0 ξ + = lim T →0 ξ(t + T ) = ξ(t) and (47) holds, the system (52) is equal to the system (21) when T → 0. Hence, the following result holds. Theorem 6: Consider the system (52). Then, the following claims hold for any T > 0.
Theorem 6 ensures that the manifold (6) is globally attractive and that the solutions of the system exist for any T > 0. We consider now the stability properties of the zero equilibrium of systems (45) and (52). The following result provides a discretetime equivalent of Theorem 3.
Theorem 7: Let the assumptions of Theorem 3 hold. Let U ⊆ Rn × B be a closed set that contains the origin. Let T be such that in which ϕ T i represents the ith entries of the vector ϕ T , and (52) is used, where T is such that (49) holds. Then, for this selection of T and for δ and selected as in Theorem 3, there exists a neighborhood of the origin U ⊆ U such that the following claims hold: 1) ξ(t) and η(t) exist, for all (ξ(0), η(0)) ∈ U ; 2) lim t→+∞ ξ(t) = 0, for all (ξ(0), η(0)) ∈ U ; 3) lim t→+∞ η(t) = 0, for all (ξ(0), η(0)) ∈ U . Proof: Consider first the discretization scheme (45). The proof of claims 1) and 2) trivially follow from Theorem 5. To prove claim 3), consider the positive definite function V (ξ, η) = W (η) + ξ ξ. The variation of V (ξ, η) along the trajectories of the system in a time interval T is where D T (ξ, η) is defined in (50) and Q i (ξ, η) is the ith row of the matrix Q(ξ, η). Note now that By replacing the previous equation in (58) and by applying the Young inequality, yields where the last inequality is obtained by replacing the expression of (51) in D T (ξ, η). Since T ≤ T 2 , with T satisfying (49), then −1 < T k i (ξ, η) < 0, and as a consequence Note that, under assumption (57), and by selecting δ ≥δ wherē δ is given in (29), for some 0 < β < −γ, yields and, under assumption (56) and by selecting ≥¯ wherē η). Hence, ΔV ≤ T (γ + β)V (ξ, η) < 0, for all (ξ, η) = (0, 0), proving claim 3). Consider now the discretization (52). Claims 1) and 2) follow from Theorem 6. The proof of claim 3) follows the same steps of the previous proof with the difference that D T has the expression in (54) and (55). Hence, the inequality (59) becomes By means of simple calculations, it is easy to show that the inequality T . Therefore, the claim follows by applying the same steps in (60) and (61).
Remark 12: While the discretization scheme in (52) ensures global attractivity of the solution manifold, for the discretization scheme (45), attractivity holds only locally unless q i , for i = 1, . . . , r are bounded, see Theorems 5 and 6. However, this comes at the cost of a smaller sampling time required to preserve the stability properties of the equilibrium points, see the assumptions of Theorem 7.
Remark 13: For the general case in which m > 1, the mappings K(ξ, η) and Q(ξ, η) in (45) and (52) can be written as K(ξ, η) = diag(k 1 1 (ξ, η), . . . , k m r m (ξ, η)) and Q(ξ, η) = [q 1 1 , . . . , q m r m ], respectively, where k j i (ξ, η) is given in (34). Observe also that the claims of Theorems 5-7 remain valid whenever k i (ξ, η) in (49) is replaced by k j i (ξ, η), for j = 1, . . . , m. The proof follows the same steps as the case m = 1, in which d i , q i , and r are replaced by d j i , q j i , and r j , respectively. Remark 14: By using the concepts of relative degree and zero dynamics for sampled nonlinear systems (see [51] and [52]), and with similar arguments as those in Lemma 3 and in Theorem 4, the claims of Theorems 5-7 can be restated in the original coordinates system.
The results presented in Sections II-IV outline a procedure to simulate DAE systems, as described in the following algorithm.

Algorithm
Data: 1) f , g, and h as given in (1) and such that the index ν is well defined; 2) x 0 = x(0) ∈ R n and a function W such that assumptions of Theorem 3 holds. Output: x + such that 1) claims 1) and 2) of Theorem 5 hold if the discretization scheme (45) is used; 2) claims 1) and 2) of Theorem 6 hold if the discretization scheme (52) is used; 3) claims 1)-3) of Theorem 7 hold if x • is locally asymptotically stable.
. . , r and complete the set of φ i (x) for i = r + 1, . . . , n such that the Jacobian of Φ(x) is nonsingular.
Step 5. Compute T such that (49) holds and T such that (56) and (57) hold. Step 6. IF γ 0 = 0 and γ < 0 in Step 3, THEN Select a sampling time T ≤ min( T 2 , T ) OR T ≤ min( Step 7. Accordingly to the choices in Step 6, advance the solution from the approximate state (ξ, η) to the approximate state (ξ + , η + ) according to scheme (45) OR according to scheme (52), in which ϕ T (q) is any one-step scheme such that (47) holds.
Consider now the positive definite and radially unbounded function where is such that (25) holds for any c ≥ 0 and some σ > 0. Note first that for some M > 0. The time derivative of (74) is such thaṫ . Hence, the function W satisfies both assumptions (25) and (26) of Theorem 2. Different simulations have been carried out to evaluate the performances of the proposed stabilization methods. Consider, for instance, the case c = 0 and the initial conditions θ 1 (0) = 0.9999 4 π, θ 2 (0) = 3 4 π, ξ 1 = −0.0001, and ξ 2 = 10. The system (1), (64) and (70) has been integrated in MATLAB by means of the solver ODE45 and default absolute and relative tolerances, using the following three stabilization methods: 1) Baumgarte's method [33] with parameters α 1 = 5 and α 2 = 25 [see (9) and (10)]; 2) the method in [38] with parameter γ = 10 [see (18)]; 3) the proposed method with parameters δ = 1 and = 0. 1 [see (21) and (22)]. Simulation results are plotted in Figs. 2 and 3. In particular, in the subplots (a) and (b) of Fig. 2, it is shown that, for this choice of parameters, the trajectories diverge in finite time for both methods 1) and 2) (dashed-red and dotted-green lines, respectively), while method 3) (solid-blue line) guarantees the existence of the solution for the entire time interval. The time histories of ξ 1 and ξ 2 are shown in the subplots (c) and (d), respectively, of the same figure. We observe that the trajectory obtained with method 3) (solid-blue line) has a faster transient response than the trajectories obtained with methods 1) and 2) (dashed-red and dotted-green lines, respectively). Note that a faster transient response could be achieved using methods 1) and 2) with a different choice of the parameters. However, for any fixed choice of these parameters, there exists a sufficiently large inconsistent initial condition such that the corresponding trajectory diverges in finite time. The trajectories in the Cartesian plane of the masses m 1 and m 2 are shown in the subplots (a) and (b) of Fig. 3, along with their zoomed-in views in (c) and (d), respectively, which reveal a greater numerical accuracy for the method 3).
The system (5) and (71) has been discretized with the schemes (45) and (52) proposed in Section IV in combination with the  forward Euler scheme for the zero dynamics [see (48)]. Simulations have been carried out with the same set of parameters and initial conditions as the previous one. In this setting, it is observed that the condition (49) of Theorem 5 is satisfied for all T ≤ T = 0.0001. A comparison of the trajectories obtained with the two proposed discretization schemes, for different values of the sampling time T , is shown in Fig. 4. Time histories of ξ 1 and ξ 2 obtained with the scheme (45) are shown in the subplots (a) and (b), while the same variables obtained with the scheme (52) are in the subplots (c) and (d) of the same figure. As expected, constraints stabilization is ensured for T ≤ 0.0001 when the scheme (45) is used. For 0.0001 < T ≤ 0.002, stabilization of ξ 2 is still achieved with good accuracy [see Fig. 4(b)], while ξ 1 diverges for T = 0.002 [see Fig. 4(a)]. On the other hand, when scheme (52) is used, constraint stabilization in ensured for any T > 0, see Fig. 4(c) and (d). Consider now the case c = 3, i.e., the equilibrium x • 2 is locally asymptotically stable. Let and consider again the function (75). In this setting, the condition (25) of Theorem 2 is satisfied, for instance, for the parameters γ 0 = 0 and for any γ ≥ −0.001. Moreover, for any selection of the parameter −γ > β > 0, there exists some  σ > 0 in (75) such that M in (32) is arbitrarily small. It is observed that the condition (49) of Theorem 5 is satisfied for any T < T where T = 0.0033. In addition, the condition (56) of Theorem 7 is satisfied for any T < T where T = 0.0045, while for any T > 0, there exists some σ > 0 such that the condition (57) holds for M arbitrarily small. Simulations have been carried out by integrating the system (5) and (71) with the discretization scheme (52) in combination with the forward Euler