Feedback Stabilization of First Order Neutral Delay Systems Using the Lambert W Function

: This paper presents a new approach to stabilize the ﬁrst order neutral delay di ﬀ erential systems with two time delays. First, we provided a few oscillation and non-oscillation criteria for the neutral delay di ﬀ erential equations using spectrum analysis and the Lambert W function. These conditions were explicit and the real roots were analytically expressed in terms of the Lambert W function in the case of non-oscillation. Second, we designed a stabilizing state feedback controller for the neutral delay di ﬀ erential systems with two time delays, wherein the proportional and derivative gains were analytically determined using the results of the non-oscillation criteria. A few examples are given to illustrate the main results.


Introduction
Many practical problems can be described using delay differential equations (DDE), as time delays are intrinsic to many real systems. Additionally, multiple sensors, actuators, controllers, shared communication networks, and haptics introduce multiple time delays [1][2][3][4]. The delays, which are either fixed constants or distributed delays, play an important role in real systems. Neutral delay differential equations (NDDE) are dynamic systems, wherein the delays appear in both the state variables and their time derivatives. Practical applications of NDDE can be seen in circuit theory [5], mechanical systems [6][7][8], neural networks [9], bioengineering [10], economics [11], and control engineering [12][13][14]. Recent research topics in the field of delay differential systems are stability analysis, bifurcation analysis, control theory and application, numerical computation, and periodic solutions analysis.
Stability analysis and stabilization of neutral delay differential systems (NDDS) is a subject of great practical and theoretical importance. This field has received considerable attention over the decades because a time delay might change the performance and stability of a control system considerably. Moreover, the NDDS exhibit drastic differences in their characteristic behavior as compared to the retarded delay differential equation (RDDE), wherein the delay argument does not occur in the highest order derivative term. Comprehensive research on this subject can be found in the literature [14,15]. However, the solution of delay differential systems is complex due to the presence of exponential type transcendental terms in the system characteristic equation. Conventional approaches used for the stability analysis of NDDS are based on the Lyapunov function. These techniques seek to determine the solution of a linear matrix inequality problem. The spectrum approach used to solve the stability problem is mainly based on the solution of an eigenvalue problem. This provides a more accurate stability analysis with less computational burden. The authors in [16], have investigated the stability of the first order RDDE using the Lambert W function, which has been implemented in all major mathematical software packages (for example, MATLAB, Mathematica, Maple, and R).
Blyuss [17] has developed a time-delayed proportional control law that stabilizes the unstable second order neutral delay systems using the parameter space approach. A delayed proportional derivative (PD) controller with a single constant delay for the second-order undamped system has been presented in [18]. In the study described in [19], the authors showed that a non-descriptor type neutral time delay system could be µ-stabilized using derivative feedback under certain conditions. The authors in [20], have designed a proportional state feedback controller for the NDDS with multiple time-delays by the pole placement method. The radius of the essential spectrum of the neural equation solution operator was computed, and the stabilization problem of the neutral equation was reduced to a problem involving only a finite number of characteristic roots. A Lyapunov function-based stability criterion for the scalar NDDE with a single delay has been presented in [21].
In the last few decades, there has been considerable interest in deriving new sufficient conditions for the oscillation or non-oscillation of various types of first-order NDDE solutions [22][23][24][25][26]. This has been the case due to its theoretical interest and important applications. A few sufficient conditions for NDDE with constant coefficients have been presented in [22][23][24]. A number of authors have established a few sufficient conditions for oscillation and non-oscillation using an auxiliary solution [25,26].
In this study, we developed a new stabilization method for first order NDDS with two time-delays. The proposed method consists of two parts: (1) sufficient conditions for the oscillation and non-oscillation of solutions of the NDDE using the geometric approach, and (2) a stabilization method using a state feedback controller. The Lambert W function is well applicable to the DDE with a single delay. Furthermore, the roots of its corresponding characteristic equation can be expressed in terms of the Lambert W function. An outline of the first part of the method is that firstly we define two functions f (λ) and g(λ), which are transcendental equations, such that the characteristic equation p(λ) is the difference of two functions, i.e., p(λ) = f (λ) − g(λ). The solution sets of f (λ) = g(λ) are then equal to the roots of the characteristic equation. Secondly, the geometric properties of the introduced functions are obtained using the Lambert W function. The solutions of the characteristic equations are the intersection points (or tangent points), of the two defined functions. This enables us to easily derive the sufficient oscillation conditions, and the sufficient conditions for the existence of non-oscillations of the characteristic equation using the Lambert W function. Therefore, the solutions of the NDDE are analytically expressed in terms of the Lambert W function. In the second part of the method, we designed a stabilizing state feedback controller for NDDS with two time-delays, where the proportional and derivative gains were analytically determined using the results from part one. This paper is organized as follows: In Section 2, the problem is formulated and a few preliminary concepts are detailed. In Section 3, we present the oscillation and non-oscillation criteria for the NDDE and design a stabilizing state feedback controller. In Section 4, some illustrative examples are given to show the effectiveness of the proposed method. This paper ends with a few concluding remarks.

Preliminary and Problem Formulation
The Lambert W function (hereafter denoted as W(·)), is a multivalued complex function defined as the solution of the following equation: where, z is a complex number. The Lambert W function has infinitely many complex branches when z is a complex number [27]. If z is a real number and z ≥ −e −1 , then W(z) becomes a real function with two possible real branches. By convention, the branch satisfying W(z) ≥ −1 is considered to be the principal branch (W 0 (z)), while that satisfying W(z) < −1 is the negative branch (W −1 (z)). Several function values of W(z) are W(−e −1 ) = −1, W(0) = 0, and W(e) = 1.
In this paper, we consider the first order NDDS with two time delays of the form: . where b, τ 1 , τ 2 are positive constants, a is a negative parameter, and u(t) is a control input. The autonomous NDDE is given by: .
The characteristic equation associated with Equation (3) is given by: where, λ ∈ C are the characteristic roots. All solutions of Equation (4) are oscillatory if and only if the corresponding characteristic equation has no real roots [24].
Here, we define two functions as follows: Using Equations (5) and (6) we can rewrite the characteristic equation as: Thereafter, p(λ) = 0 results in the following equation: Hence, the solution set of Equation (8) is equal to the roots of Equation (4). Throughout this paper, we assume that b, τ 1 , τ 2 are positive constants and Re(λ) denotes the real part of a complex variable λ.

Oscillation and Non-Oscillation Criteria
In this section, we discuss the oscillation and non-oscillation criteria for Equation (4). The non-oscillation criteria indicate the definite existence of some real roots. Lemma 1. If a < 0, then the following are true: Therefore, the global minimum can be found by setting its first derivative to zero.
Proof. The second derivative of f (λ) is given by: Proof. The proof is similar to that of Lemma 2.

Remark 1.
The above lemmas show that f (λ) is not a convex function on Re(λ). As g(λ) is infinitely differentiable and its second derivative is positive, g(λ) is a concave function for all Re(λ).
then every solution of Equation (3) oscillates.
Proof. We obtain from Equation (14), in several steps, the following mutually equivalent inequalities: This completes the proof.
then every solution of Equation (3) oscillates.
Proof. The proof is similar to that of Theorem 1. The maximum value of g(λ) on the interval 0, 1 τ 1 ln(−a) is g( 1 τ 1 ln(−a)) as g(λ) is a concave function. If f (λ m ) is greater than g( 1 τ 1 ln(−a)), then Appl. Sci. 2019, 9, 3539 6 of 12 they do not intersect at any point. We proceed as in the proof of Theorem 1 to immediately obtain the following result: then Equation (3) has two negative real roots.
Proof. First, we determine the maximal number of real roots of the quasi-polynomial function (Equation (4)). From the Pólya-Szegö bound (see Theorem A1), setting ξ = θ = 0 gives a bound on the number of real roots for (4). Here, d 1 = 2 and d 2 = 1, therefore, the maximal number of real roots is 2. Next, two functions f (λ) and g(λ) intersect in the third quadrant of the complex plane if Equation (26) is satisfied. As g(λ) is a strictly increasing function, and f (λ) strictly decreases on the interval 1 τ 1 ln(−a), λ m and strictly increases on the interval (λ m , 0), then g(λ) intersects two points of f (λ) if f (λ m ) < g(λ m ). This completes the proof.
then Equation (3) has two positive real roots.
Proof. The proof is similar to that derived for Theorem 4.

Stabilization of the NDDS
In this section, we consider a state feedback controller of the form: where, k p and k d are the proportional and derivative gains, respectively. Our objective was to analytically determine the set of controller parameters that makes the NDDS (2) exponentially stable. Including the state feedback controller (28), the NDDS (2) can be written as: .
Therefore, the characteristic equation of the NDDS (29) is expressed as follows: where, a f = a − k d and b f = b − k p . Using Theorem 4, we first determine k d such that −1 < a f < 0 and subsequently, we find k p from Equations (8) and (13). The main results of this section are summarized in the following theorem. Theorem 6. If a < k d < a + 1 and where, λ m = 1 τ 1 1 − W(−ea −1 f ) , then the state feedback controller (28) makes the NDDS (2) exponentially stable.

Examples
A few examples and numerical simulation results have been presented here to validate the analytical results of the previous section.
where, a = −2e, b = √ e, τ 1 = 1, τ 2 = 0.5. As a < −1, from Theorem 3 we obtain: Equation (34) does not satisfy the condition of Theorem 3, but satisfies the condition of Theorem 5. Thus, there are two positive real roots. The results are shown in Figure 1c.
The solid red, blue and black lines in Figure 1 represent f (λ), g(λ), and the characteristic equation p(λ), respectively. The intersections of the two functions (red and blue lines), denoted by the blue dots, are the real roots of the corresponding equation. The intersection points between the characteristic equation and the real axis exactly equal the blue dots. Figure 1a shows that there is actually no intersection between p(λ) and the real axis, which implies that every solution of Equation (3) oscillates.   The solid red, blue and black lines in Figure 1 represent ( )  Figure 1a shows that there is actually no intersection between ( ) p λ and the real axis, which implies that every solution of Equation (3) oscillates. Figure 1b,c establishes the existence of real roots.

Example 4. Consider a NDDS with state feedback:
.
and u(t) = k d .
The characteristic equation of Equation (35) in the absence of control becomes a simple quasi-polynomial: Applying Equation (27)      We used a state feedback controller where the controller gains were analytically determined from Theorem 6, as one of the goals of this study was to stabilize the NDDS (35). First, we selected the derivative gain 2.5 Next, we can obtain  We used a state feedback controller where the controller gains were analytically determined from Theorem 6, as one of the goals of this study was to stabilize the NDDS (35). First, we selected the derivative gain k d = −2.5 that satisfies the condition a < k d < a + 1, such that −1 < a f (= a − k d ) < 0. Next, we can obtain λ m = 1 (13), and finally the proportional gain was determined using Equation (31). Figures 2b and 3b show the rightmost spectrum and time response, respectively of the NDDS with the state feedback controller. Applying the QPmR to the NDDS with a feedback system, the four rightmost roots were determined to be −0.2652, −0.5355, −0.9979 + i35.9105, and −0.9979 − i35.9105, which indicates that all the roots of its characteristic equation are located to the left of the imaginary axis. The characteristic equation of Equation (35), with the state feedback controller, becomes λ − 0.5λe −0.7λ + 0.0768e −1.2λ = 0, and the condition described by Equation (26) gives a value of −0.0645668. Based on Theorem 4, the NDDS with the above state feedback controller has two negative real roots. It should be noted that the calculated root λ m equals one of the rightmost roots produced by the QPmR. Therefore, we can assign a pole of the NDDS by adjusting the controller gains.

Example 5.
Consider the following single-species population growth with time delay and feedback control [30]: where N(t) is the population, r is the intrinsic growth rate, τ is the recovering time, c is the growth rate associated with the growth rate at time t − τ, and K is the environmental carrying capacity. Parameters r, τ, c, K are positive. In Equation (37), if q(·) = 0, it is sometimes called a dynamic food-limited population model [31].
Compare Equation (40) with Equation (4), we see in this case, a = rc, b = r, τ 1 = τ 2 = τ. We assume that the feedback controller is in the following form: A dynamic bacteria population model which is represented by an autonomous system of Equation (37), was introduced and discussed in [31]. According to [24], if r, c, τ are positive constants, then every non-oscillation solution of the autonomous system of Equation (37) tends to zero as t → ∞ . However, the real root is not guaranteed to be the rightmost pole. Table 1 shows the five rightmost poles as a function of r, c, and assuming K = 1 and τ = 0.1. The odd columns are without feedback control, and the even columns are with feedback control. In the second and third cases, without feedback control, some complex-conjugated roots are found on the right-hand side of the complex plane and they make the system unstable. With a state feedback controller designed by Theorem 6, all roots are located on the left-hand side of the complex plane. Based on Theorem 4, Equation (40) with the state feedback controller has two negative real roots, and these roots are the rightmost. The simulation results show that the designed state feedback controller makes the NDDS exponentially stable.

Conclusions
In this paper, we presented new oscillation and non-oscillation criteria for the NDDE with two time delays by defining two functions, wherein the difference between the two functions was the corresponding characteristic equation. These criteria provide the number and sign of the real roots of the NDDE. We designed a state feedback controller that makes the NDDS exponentially stable. Moreover, the controller gains were easily and analytically determined using the Lambert W function. Finally, we also included several examples to demonstrate the main results.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Theorem A1. [32,33]. Let τ 1 , · · · , τ N denote real numbers such that τ 1 < · · · < τ N and d 1 , · · · , d N positive integers such that d 1 + · · · + d N = D + N. Let x i,j (s) = s j−1 e τ i s for 1 ≤ i ≤ N and 1 ≤ j ≤ d i . Let # be the number of zeros of the function x(s) = 1≤i≤N,1≤ j≤d i c i,j x i,j (s), that are contained in the horizontal strip ξ ≤ Im(s) ≤ θ. Assuming that This result gives a bound for the number of quasi-polynomial roots in any horizontal strip. Setting ξ = θ = 0, Theorem A1 yields # ≤ D + N − 1, where D stands for the sum of the degrees of the polynomials involved in the quasi-polynomial function x(s), and N designates the associated number of polynomials.