Stability of neutral delay diﬀerential equations with applications in a model of human balancing

. In this paper the exponential stability of linear neutral second order diﬀerential equations is studied. In contrast with many other works, coeﬃcients and delays in our equations can be variable. The neutral term makes this object essentially more complicated for the study. A new method for the study of stability of neutral equation based on an idea of the Azbelev W-transform has been proposed. An application to stabilization in a model of human balancing has been described. New stability tests in explicit form are proposed.


Introduction
The main object of this paper is the second order neutral delay differential equation with the corresponding initial functions x(ξ) = ϕ(ξ), x (ξ) = ψ(ξ), x (ξ) = η(ξ) for ξ < 0, the following directions. Although the study of delay equations historically started with equations with constant delays θ i , τ i (i=1,...,m), this case does not look naturally in modeling control of technical devices. The first development is in the study of equation with variable delays in the frame of applications. In real applications the measurements of the process at discrete moments of times are used. It means that the deviating arguments are actually of the form t − θ i (t) = g i ,t − τ i (t) = h i , t ∈ (t i , t i+1 ), i = 1, 2, .., m, where g i ,h i are constants (i = 1, 2, ..., m). Simulation of control of smart prostheses, which is one of most important applications of the human balancing model, discussed below in Section 4, is based in fact on the values of angles and speeds of their changes is also based on delays of this sort. In many cases variable delays appear naturally in mathematical models. For example, control can depend on geometry of electromechanical device, for instance on the distance between sensors, mechanical and control blocks. If the distances can change we come to variable delays in the model. The second direction of our development of the model of human balancing [15,23] is in considering variable coefficients. The mechanical stiffness in biological systems is rather periodic than constant and this in a corresponding sense explains considering variable coefficients. But there is also another reason. Below, discussing the model of human balancing in Section 4, we explain additional ways of appearing variable coefficients in this model. It will be demonstrate, for instance, that a combination of small external unpredictable influences and nonlinearities can imply variable coefficients in the model. A "small" external force can change the position of the stationary point, and this could lead to the fact that the coefficients of the system linearized in the neighborhood of this new stationary point change because they depend on the solution. Example 4.1 demonstrates that a system without a right-hand side linearized in the neighborhood of the solution x = 0 can be exponentially stable, and a system linearized in the neighborhood of x = δ , where δ is a "small" constant, may be Lyapunov unstable. It is natural in these cases to consider equations with variable coefficients instead of equations with constant ones. If a person falls, then the theoretical fact about stability gives nothing. This demonstrated the importance of considering non-homogeneous equations. In practical tasks, the influence of the right-hand side on the solution should be evaluated. It is necessary to obtain a dependence of the solution on the right-hand side bounded by corresponding constants. We will get such estimates in Section 3. This is the third direction in the development of human balancing model. The principal development of the results presented in [3] is in considering cases of oscillation behavior of solutions to neutral equation (1.1), where positivity-based approach does not work (see Remark 3.3). The results of [3] cannot be used for the model of human balancing, see, for example, equation (4.10), in which the coefficient of ϕ(t) is negative and the assumption of the main result of [3] about exponential stability of the ordinary part ϕ + Aϕ + Bϕ = 0 is not fulfilled. Even in the case of θ 2 = 0 in equation (4.10), the main result of [3] does not work in the cases of such small q i and large coefficients of x and x .
In this paper we propose a method for stability analysis of neutral equations. This method leads us to conditions on smallness of delays that is naturally from the physical point of view. Our paper is built as follows. In Section 2 we describe known results which are used in the proofs. In Section 3 we formulate results on stability of equation (1.1). In Section 4 we describe the known model of human balancing and our developments concerning variable delays and coefficients from mechanical point of view. In Section 5 applications of our results to estimates of delays in the model of human balancing are obtained. In Section 6 examples for estimate of acceleration coefficient for the human balancing are presented. In Section 7 we prove main results. In Section 8 discussion and open questions are presented.

Preliminary
We understand a solution of equation ( where x 1 (t), x 2 (t) are two solutions of homogeneous equation (2.2), (1.2), where 2) satisfying the conditions respectively. The kernel C(t, s) in representation (2.1) is called the Cauchy function of equation (1.1).
Note that that in the classical books on delay differential equations [9,19], the homogeneous equations are considered as equation (2.2) with absolutely continuous initial function It is known that for equation (

Formulations of Main Results
Let us consider the following ordinary differential equation with constant positive coefficients A and B, where z is an essentially bounded measurable function.
Denote by W (t, s) the Cauchy function of equation (3.1). It is known that for every fixed s the function W (t, s), as a function of the variable t, satisfies the homogeneous equation (3.2) and the initial conditions We have to construct W (t, s) in every one of the cases: if A 2 > 4B then and if A 2 < 4B then It is known that the solution of the equation (3.1) which satisfies the initial conditions can be written in the form Its derivatives are the following Let us denote Denote by x(t) the solution of homogeneous equation (2.2) and by X(t) the solution of nonhomogeneous equation (1.1) with the same initial conditions, and f * = esssup t≥0 |f (t)|.
Theorem 3.1. Let A > 0, B > 0 and the following inequality be fullfilled: Remark 3.1 Denoting P = ||W ||, Q = ||W t ||, R = ||W tt ||, we have a simple geometrical interpretation of this result. Let us define the coordinates: is fulfilled for every internal point (X, Y, Z) of the pyramid bounded by the planes Remark 3.2 It is clear from inequality (3.13) that in the case, when the coefficients a i (t) and b i (t) (i = 1, ..., m) are close to constants, the second and fourth terms are small, and in the case of small delays θ i (t) and τ i (t) (i = 1, ..., m), the first and third terms are small. We can make a conclusion that equation (2.2) in this case preserves the property of the exponential stability of equation (3.1).
Let us formulate several tests of the exponential stability based on estimates of W , W t and W tt . Proposition 3.2. Assume that A 2 > 4B and: Then equation (2.2) is uniformly exponentially stable, and the difference of solution x(t) of homogeneous equation (2.2) and the solution X(t) of non-homogeneous equation (1.1) with the same initial conditions satisfies the inequality Assume that A 2 = 4B and: then equation (2.2) is uniformly exponentially stable, and the difference of solution x(t) of homogeneous equation (2.2) and the solution X(t) of non-homogeneous equation (1.1) with the same initial conditions satisfies the inequality Proposition 3.4. Assume that A 2 < 4B and: then equation (2.2) is uniformly exponentially stable, and the difference of solution x(t) of homogeneous equation (2.2) and the solution X(t) if non-homogeneous equation (1.1) with the same initial conditions satisfies the inequality Remark 3.3 In the paper [3], the exponential stability of neutral equations was obtained under the assumption A 2 ≥ 4B. In Proposition 3.4 in contrast with this assumption, the case of A 2 < 4B is considered. This means that we do not limit our analysis only by the case of nonoscillatory equation x + Ax + Bx = 0 and consider also the case of oscillating model equation.

Description of neural-mechanical model of human balancing
Let us explain how second order delay equation appears in the model of balancing. The neural-mechanical model of human balancing in the sagittal plane is depicted in Figure  1 which was constructed in [23]. If we consider point A as fixed point we can model the human body as an inverted pendulum with mass m, moment of inertia J A with respect to pivot A, while l stands for the distance between the center of gravity C and pivot A. b is the torsional dashpot of damping and k is the torsional spring of stiffness. The torque Q is regulated by the central nervous system based on the sensory signals about rotation angle ϕ, angular velocityφ and angular accelerationφ of the human body. By using the Newton second law for the moment at fixed point A we will obtain the following equation where

2)
Q p is a passive torque, Q a -an active torque, f (t) defines corresponding torques which is implied by actions which seem negligible and is not usually taken into account in the construction of the model (for example, a person standing on a slowly moving belt or overboard, a weak wind that increases resistance to motion, etc). The passive torque is defined as In physiological systems, the coefficients b(t), k(t) are rather periodic than constant and small delays θ 1 (t), θ 2 (t) exist. In [23] θ 1 and θ 2 are considered as negligible and the coefficients b(t) = b, k(t) = k are constants. The active torque Q a (t) is assumed in the nonlinear form where α denotes the limit of the active torque and its linear part is with K p , K d and K a being the positive constants describing controls. In physiological systems they are rather almost periodic functions than constants. τ 1 ,τ 2 ,τ 3 are delays of these torques respectively. The delayed signals of angle ϕ and angular velocityφ are provided by the vestibular system and proprioceptors, while the angular acceleration ϕ signal is related to the information coming from the mechanoreceptors according to Newton's Second Law [15]. The saturation effect is shown in Figure 2 which is taken from [23], where the active control torque Q a tends to saturation torque limit α as Q c increases. The dynamics of the quiet standing process is governed by the following secondorder nonlinear delay differential equation: These nonlinearities can be presented by variable coefficients in the model i(t, ϕ) = 1 + ∆i(t, y), K p (t, ϕ) = K p +∆K p (t, y), K d (t, ϕ) = K d +∆K d (t, y), K a (t, ϕ) = K a +∆K a (t, y), where ∆i(t, y) → 0, ∆K p (t, y) → 0, ∆K d (t, y) → 0, ∆K a (t, y) → 0 when y → 0. Note that tanh y y Making the classical linearization we obtain Formulas (4.8) allow us to obtain estimates of ∆i(t, ϕ), ∆K p (t, ϕ), ∆K d (t, ϕ) and ∆K a (t, ϕ), if we have a priori estimates of the solution. Analysis of stability can be reduced to one of the homogeneous equation If we assume that there are no delays in reactions, then in the case of K p + k > mgl, the second order ordinary differential equation is exponentially stable. In practice, it is known that the mechanical stiffness is insufficient for stability of human balancing [14]. The source of instability is the delay in physiological reactions. To estimate a smallness of the delays which does not break stability is the main goal in this model. We make this estimation in Section 5. Let us take a closer look at the standard transition from equation (4.6) to linearized equations (4.7) and (4.10). Suppose the function x * (t) is a solution to the nonlinear equation (4.6). Is it physically achievable? It is clear that for this it must be at least Lyapunov stable. Does the exponential stability of equations (4.10), which was obtained according to the standard linearization scheme, imply exponential stability, or at least Lyapunov stability of this solution x * (t)? Generally speaking, no. taking into account two first terms in every one of presentations (4.9) and not one term only as it was done in the classical linearization (4.10). We obtain Let us set f (t) = −60(ϕ * ) 3 , where ϕ * will be defined below. Assume for simplicity that After substitution of the coefficients we obtain where coefficients are defined by (4.12). The constant ϕ * = √ 38 √ 14399997 is a solution of this equation. The coefficient of ϕ(t) is negative. In this case one of the roots of the characterictic polynomial for the homogeneous equation is positive and there are no even the Lyapunov stability of the solution ϕ * = √ 38 √ 14399997 . Instead of assumptions (4.13) we can suggest that the delays are small enough.
So we get the following equation In the case of θ 1 = 0 and θ 2 = 0 we obtain This equation is a particular case of equation (1.1). Thus we come after the classical linearization to the following homogeneous equation for stability studies

Application to model of human balancing: estimates of delays
Let us study the stability of (4.17) and (4.18) on the basis of Theorem 3.1 and Propositions 3.2-3.4. The trivial solution ϕ = 0 of (4.19) corresponds to the equilibrium of balancing. Let us describe the numerical values of several physical parameters [23]: human weight mg = 600, length from pivot point A to center of mass C l = 1, moment of inertia J A = 60 and passive stiffness k = 471. The contribution of the passive dumping is usually small, it is assumed to be negligible in further analysis: b ≈ 0. In our terms we get that a 1 = 0 and b 1 = −2.15. We will give range for the proportional, derivative and acceleration gains K p , K d and K a as follows [23]: To study the balancing of inactive people, we will choose the maximum values for the demonstration. We will also choose the reaction delays to be fixed at τ 1 ∈ [0, 1] and τ 2 ∈ [0, 1].
Below we use A and B defined by (3.12) in the case m = 2 According to Theorem 3.1, the inequality that implies exponential stability for equation Substituting all the terms described in the end of the previous section, namely mgl = 600, J A = 60, k = 471, b = 0, K p = 3000, K d = 600 and K a = 48, we obtain the case of A 2 < 4B. We will evaluate W t and W tt by Proposition 3.4. We can get the area of stability below the blue line at Figure 3.
However, we can also consider the cases where θ 1 and θ 2 are not zeros (this case can be important, for example, for people with disabilities). If we use the same steps of the proof like in Theorem 3.1 (see Section 7) we get that the following inequality Figure 3: Domain of stability is below the blue line: 1.728X + 19.125Y < 1 Figure 4: Domain of stability with θ 2 is inside the pyramid: 1.728X+19.125Y +0.271Z < 1 implies the exponential stability for equation (4.17). Note that from the mechanical point of view a 1 is sufficiently small.
We can see that for small θ 2 we get similar results to inequality (5.3). If θ 2 is variable, for example 0 ≤ θ 2 (t) ≤ 1, we can see at Figure 4 that the area of stability is inside the pyramid.

Examples on estimates of acceleration coefficient q
Let us demonstrate few examples that allows us to estimate q in the human balancing model.
In this case we will need to use Proposition 3.2 for the stability analysis. Let us take numbers that match this case: A = 8 and B = 7.85. Now if we substitute the physical parameters (see Section 4 for details) in Proposition 3.2 we will get the following inequality Now by substituting mgl = 600, J a = 60, k = 471, b = 0, K d = 480, K p = 600 and τ * 1 = τ * 2 = 0.05 (see Section 4 for details) we can evaluate the range of q: In this case we will need to use Proposition 3.3 for the stability analysis. Let us take numbers that match this case: A = 4 and B = 4. Now if we substitute the physical parameters (see Section 4 for details) in Proposition 3.3 we will get the following inequality Now by substituting mgl = 600, J a = 60, k = 471, b = 0, K d = 240, K p = 369 and τ * 1 = τ * 2 = 0.05 (see Section 4 for details) we can evaluate the range of q: |q| < 0.1343 (6.8) Example 6.3. The inequality A 2 < 4B (6.9) In this case we will need to use Proposition 3.4 for the stability analysis. Let us take numbers that match this case: A = 10 and B = 40. Now if we substitute the physical parameters (see Section 4 for details) in Proposition 3.4 we will get the following inequality (6.11) Now by substituting mgl = 600, J a = 60, k = 471, b = 0, K d = 600, K p = 2529 and τ * 1 = τ * 2 = 0.05 (see Section 4 for details) we can evaluate the range of q: |q| < 0.4886 (6.12)

Proofs
In the begining of Section 3 we defined the Cauchy function W (t, s) of ordinary differential equation ( Proof of Theorem 3.1: Consider the following equation where r i (t) = t − φ i (t), g i (t) = t − τ i (t) and h i (t) = t − θ i (t) with the zero initial functions x(ξ) = x (ξ) = x (ξ) = 0 for ξ < 0. (7.11) It is known [1] that in the analysis of stability, we can consider only the zero initial values Let us start with the case of g i (t) ≥ 0, h i (t) ≥ 0. We can rewrite equation (7.10): and Let us make the Azbelev W -transform [1], substituting where z ∈ L ∞ (L ∞ in the space of essentially bounded functions z : [0, ∞) → (−∞, ∞)), in equation (7.14). We obtain the following equation: where the operator K : L ∞ → L ∞ is defined by the equality where Define the operator |K| : t 0 |W tt (t, s)|z(s)ds + z(t) . In the general case (i.e. without the assumption g i (t) ≥ 0, h i (t) ≥ 0 for t ≥ 0) we can write equation (7.14) in the form ∆b i (t)x(h i (t)) + q(t)x (r i (t)) =f (t), t ∈ [0, +∞), (7.19) wherẽ x (s)ds. (7.20) It is clear thatf (t) is bounded. The exponential estimates of the Cauchy function in the first case imply the boundedness of solution x(t) for t ∈ [0, ∞) for every bounded on the semiaxis [0, ∞) right-hand side f (t). The Bohl-Perron theorem implies the exponential estimate of the Cauchy function and nontrivial solutions. It follows from the fact ||K|| < 1 that |z| ≤ 1 1−Q and |x(t)| ≤

Open questions
It would be interesting to develop our approach to the study of stability of the neutral delay differential equation of the order n: a i n−1 x (n−1) (t−θ i n−1 (t))+...+ m i=1 a i 0 x(t−θ i 0 (t)) = 0 (8.1) Systems of second order equations could be studied on the basis of the construction of the Cauchy matrix. Such systems describe model in which a motion of point objects are studied. In practice, usually 3-dimensional systems have to be considered.
It would be interesting to consider the neutral equation without damping terms (i.e. without the first derivative x (t − τ i (t))). Asymptotic properties of equation (8.2) without damping term (i.e. in the case of a i (t) ≡ 0 for t ∈ [0, ∞)) were studied in ( [19], Chapter III, Section 16, pp. 105-106), were instability of the equation x (t) + bx(t − θ) = 0, for every pair of positive constants b and θ was obtained.
Conditions of stability of second order equations with variable coefficients and delays were obtained in [7].
Results about boundedness of solutions for vanishing delays (θ i (t) → 0 for t → ∞) and about asymptotic representations of solutions were obtained in [16,20], see also ([19], Chapter III, Section 16). Boundedness of solutions for equations with advanced arguments (θ i (t) ≤ 0) was studied in [8]. Results on the exponential stability of the equation x (t)+ax(t)−bx(t−θ) = 0 with constant coefficients and delay were obtained in [5]. First results based on nonoscillation for the exponential stability of the second order equation without damping term and with variable coefficients and delays were obtained in the paper [6]. thesis.

Formatting of funding sources
Non-relevant.

Funding
Non-relevant.