Gridding discretization-based multiple stability switching delay search algorithm: The movement of a human being on a controlled swaying bow.

Delay represents a significant phenomenon in the dynamics of many human-related systems-including biological ones. It has i.a. a decisive impact on system stability, and the study of this influence is often mathematically demanding. This paper presents a computationally simple numerical gridding algorithm for the determination of stability margin delay values in multiple-delay linear systems. The characteristic quasi-polynomial-the roots of which decide about stability-is subjected to iterative discretization by means of pre-warped bilinear transformation. Then, a linear and a quadratic interpolation are applied to obtain the associated characteristic polynomial with integer powers. The roots of the associated characteristic polynomial are closely related to the estimation of roots of the original characteristic quasi-polynomial which agrees with the system's eigenvalues. Since the stability border is crossed by the leading one, the switching root locus is enhanced using the Regula Falsi interpolation method. Our methodology is implemented on-and verified by-a numerical bio-cybernetic example of the stabilization of a human-being's movement on a controlled swaying bow. The advantage of the proposed novel algorithm lies in the possibility of the rapid computation of polynomial zeros by means of standard programs for technical computing; in the low level of mathematical knowledge required; and, in the sufficiently high precision of the roots loci estimation. The relationship to the direct search QuasiPolynomial (mapping) Rootfinder algorithm and computational complexity are discussed as well. This algorithm is also applicable for systems with non-commensurate delays.

Introduction

The stability of Time Delay Systems (TDS) has been a challenging and intensively studied research topic over the past few decades [1-4]. It is, i.a., given by the fact that delay appears widely in many real-world systems—e.g. networked, mechanical, electrical, social, economical, etc. [5-7]. Moreover, delay can act as a factor that causes performance degradation and system instability; hence, the latencies and the after-effect phenomena should be considered in the models. Biological systems and human-body models are not unaffected; for instance, a robust estimator of neuronal interaction delays was proposed in [8]; the stabilization mechanism of unstable human quietstanding and motion was studied in [9-11]—where delays are, however, only marginally considered. The human nervous system with reflex delay can also be viewed as the cause of traffic jams [12]. Robust compensation of instability caused by delays for a human standing still was proposed by Asai, et al. [13]. Apart from the above-mentioned habitual studies that take fixed delay values into account; from the robustness point-of-view, it was desirable to also perform research on the influence of particular delay values on stability. This task gives rise to the notion of Delay Dependent Stability—(DDS), which means that the system remains stable in a finite set of delay intervals while taking into consideration fixed or variable system parameters.

Two basic families of DDS methods for computing the delay stability margins are paramount in the literature; namely, Time-domain Indirect and Frequency-domain Direct methods. The former group is based on the Lyapunov-Krasovskii or Lyapunov-Razumikhin approaches, and linear matrix inequalities [14-17]; the Jensen inequality approach [18]; semi-definite programming methods [19]; or the descriptor system approach [20,21]. Such methods suffer from significant computational burdens and provide purely theoretical—and rather conservative results [22]. However, they can deal with both constant and time-varying delays and are even applicable to non-linear systems in most cases. To highlight some important results regarding time-varying delay systems, the maximal allowed upper bounds of multiple time-varying delays expressed in terms of several linear matrix inequalities (LMIs) were investigated e.g. in [23,24]. In [25,26], these bounds were derived from feasible LMIs by employing MATLAB® LMI Toolbox [27] to reach average consensus asymptotically; moreover, the rate of convergence was estimated. Compared to the presented paper, these methods provide only the upper bound, the given inequalities are obtained analytically and they are more mathematically involved; however, the domain of time-varying delays provides more general results. The latter family of frequency domain methods—usually better implementable and applicable—includes approaches based on various ideas, see e.g. [28] for their listing. Namely, from the historical point of view, methods based on the Rekasius substitution, giving rise to a delay-free Associated Characteristic Polynomial (ACP), constitute the group worthy of the highest research interest within this family [29-31]. An interesting solution can be obtained via the argument principle [32]. Besides this, Schur-Cohn methods [3,33]—or those eliminating exponential (delay) terms from the Characteristic QuasiPolynomial (CQP) due to its symmetry [28,34]—have also attracted great attention from the control scientific community. The crucial step of these procedures rests in the determination of all characteristic roots (system poles), i.e. the characteristic quasipolynomial zeros located exactly on the imaginary axis—which expresses the stability margin. However, only the rightmost subset of poles makes the system switch from/to stability/instability, respectively. Hence, let the corresponding delays be called switching delays. The necessary condition for the existence of roots crossing the imaginary axis is that the so-called root tendency is nonzero [30,34]. The frequency-domain direct methods, which are of a low—or even none conservatism and computationally acceptable in many cases, were successfully applied to investigate the stability of communications [28], power [35], electromechanical [36], logistic [37]; and many other systems with constant time delay. Note, however, that they cannot be directly applied to nonlinear TDS, or those with time-varying delays.

In spite of all the favorable features introduced above, a significant disadvantage of many frequency-domain direct methods is that they are only applicable to systems (or models) with commensurate delays where the higher order of commensuracy yields the higher complexity of the analytically derived relations [6]. Motivated by the above deficiency—and also by the endeavor to use easily-accessible standard programs for technical computing (e.g. MATLAB®) by engineers and practitioners, this paper focuses on deriving a computationally simple—yet sufficiently rapid and accurate gridding algorithm for searching for delay margins in both commensurate and non-commensurate TDS with an arbitrary number of independent constant delays. Although the method is not derived analytically in terms of exact analytic formulas, but rather, in a numerical way, the techniques used enable one to retain a very low conservatism level and to reach a high degree of precision; i.e., a discrete time reformulation of the CQP by means of the iterative use of the bilinear transformation with pre-warping [38], followed by the utilization of linear interpolation to obtain the integer-powered ACP, means that the zeros can easily be computed. In addition, a linear, Regula Falsi (RF), interpolation concept is then used for the eventual switching delay estimation.

In this contribution, the designed algorithm is applied to the stability delay margin analysis of a control system with a human on a swaying bow [39]. The skater balances on a hemispheric base and remotely steers the power input to the servo that causes the angle deviation from the bow horizontal position; and consequently, the angle between the skater and the bow symmetry axis is the controlled system output. For fixed controller parameters’ values, switching delays decide about the marginal balance values that cause stability/instability feedback behavior. The model is specified, i.a., by the fact that the delay-free feedback system is unstable—which indicates that the delay margins have upper and lower bounds. Moreover, the unstable controlled plan includes two independent delays that can generally be viewed as non-commensurate ones. Note that a problem of stabilizing the skater moving on a skateboard affected by reflex delay (with a different model of the dynamics) was recently solved in [40].

The efficiency and asymptotical complexity of the proposed algorithm is benchmarked against the possibility of seeking out switching delays by using the direct QuasiPolynomial mapping Rootfinder (QPmR) algorithm [41-43] that proved to be suitable for finding the quasipolynomial zeros within a selected complex plane rectangular region. A detailed example of the use of the algorithm, which provides the reader with real time requirements and the accuracy obtained in comparison to trial-and-error switching delays search procedures, is introduced for illustration as well.

Thorough the paper, , , , denote the set of, respectively, complex numbers, real numbers, positive and non-negative integers. The n-dimensional Euclidean space is denoted as and let . For , Re(s) and Im(s) denote, respectively, the real part and imaginary part of s. The empty set is denoted as ∅. The notation d(s|a) means a polynomial in variable s, the coefficients of which are determined by the parameter a under some particular operation.

Methods

The general TDS model and its spectrum

Let, stand for the Laplace transfer function in of a single-input single-output linear time-invariant TDS with multiple delays as a frequency-domain model where n(s, τ), d(s, τ) are quasipolynomials of the general form: , in which , represent independent delays: and . Let be the associated exponential polynomial related to x(s, τ). If, , then the system defined in Eq (1) is designated as retarded; otherwise, the system is of a neutral type.

In the following section, we assume that there are no common roots of the numerator and denominator in Eq (1). Under this assumption, the system poles (characteristic roots) s agree with the zeros of d(s, τ), i.e. s ≔ {s:d(s, τ) = 0}, and thus: CQP ≡ d(s, τ). Let us denote the spectrum of the system as Σ ≔{s}, the set of all zeros of d(s,) as Σ and introduce some basic spectral properties that are utilizable in this study.

Property 1 [4,6,44]: The system defined in Eq (1) has the following properties:

If nonzero d, λ exist (where d stands for the real-valued coefficient in d(s,) for s with jth delay) for some positive τ and some i, j, k, then the number of poles is infinite.

For a retarded system, for any with |β| < ∞, only a finite number of poles are located in the half-plane Re(s) > β; while infinitely many others are located in Re(s) < β.

For a neutral system, a vertical chain of poles (s) exists at γ ≔ max ReΣ such that lim Re(s) = γ, lim Im(s) = ∞ for |s| < |s| .

Let , then the number of poles with is always finite, and they are isolated.

Isolated poles behave continuously and smoothly with respect to τ on .

If one defines the spectral abscissa as the function: α(τ)≔τ ↦ sup Re Σ, for which the following holds true (Property 2): Function α(τ) may be nonsmooth and hence, not differentiable; e.g. in points with more than one real pole or conjugated pairs with the same maximum real part [45].

It is non-Lipschitz, for instance, at points where the maximum real part has a multiplicity greater than one [46].

Rephrasing Property 2 –especially the second paragraph, abrupt changes may occur in the value of α(τ); however, such cases are rather rare—in at most, a finite number of delay vector points.

TDS stability

Proposition 1 [3,4]: The system as in Eq (1) with a fixed τ is exponentially stable if α(τ) < 0 for retarded systems; and if α(τ) < −ε, ε > 0 for neutral systems.

Note that, if γ < 0 (including infinity); then a TDS is so-called formally stable and Proposition 1 can collapse into the condition α(τ) < 0. The value of γ is, moreover, sensitive to even infinitesimal delay deviations (see item 4) of Property 1). This feature gives rise to the notion of strong stability—which means that a system is preserved from formal instability under small delay changes [4] (i.e. . The necessary and sufficient strong stability condition can be expressed as:

To avoid dealing with a potentially infinite chain of vertical system poles reaching or crossing the imaginary axis, retarded and/or strongly stable neutral systems satisfying Eq (2) are the only ones to be further considered in this study.

Let us now define some other useful notions regarding TDS stability with respect to delay values. A TDS is said to be DDS, if it is exponentially stable for some delay vectors values: , constituting open sets (regions) of stabilizing delays. The root tendency (RT) is defined as:

Crossing delays τ, and the corresponding crossing frequencies, ω, satisfy:

For switching delays , and the corresponding switching frequencies , the condition in Eq (4) holds—and moreover, the corresponding pole () on the imaginary axis is the leading (furthest right) one; i.e. , and the imaginary axis is crossed with a non-zero velocity [17,30,47], i.e.:

Let us denote the set of all switching delays and frequencies as: and , respectively. It is evident from Proposition 1, that stability may only be violated at switching delays, while poles cross the imaginary axis at switching frequencies. This is the key idea of the direct frequency-domain methods cited above. Apart from the cases listed in Property 2, such an imaginary axis crossing is smooth.

Model of a human being on a controlled swaying bow

In this study, we are concerned with a DDS stability analysis of a control system with a model of a (human) skater on a remotely-controlled swaying bow, see Fig 1. The skater controls the power input, P(t) to the servo—giving rise to the angle deviation, u(t) from the horizontal position, and consequently, the angle between the skater and the bow symmetry axis, y(t) emerges.

Fig 1

A skater on a remotely controlled swaying bow.

The skater can send remote signals to the servomechanism (the rectangle labeled as Servo) that exerts power P in order to cause a horizontal-angle input deviation of u. The horizontal asymmetry yields the deviation of angle y from the bow symmetry axis—which is then taken as the system output.

If friction is neglected, the following transfer function: (u(t) → y(t)) can be written: for some particular bow parameters, where delay τ1 expresses the skater’s reaction time, and τ2 means the latency. Nominal controlled system parameters and delay values may be read for instance as: a = −1, b = 0.2, τ1 = 0.3, τ2 = 0.1, as given in [39]. It can be computed that, for such parameters, the controlled TDS is unstable since it has the spectral abscissa value of α = 0.9534.

Let us now assume that the habitual simple negative control feedback loop is applied to a plant given in Eq (6), to restore stability. Let the control loop be equipped with a finite-dimensional linear controller governed by the transfer function: which gives the following retarded CQP:

An optimal spectral abscissa minimizing procedure for the above introduced parameters, see [48], yields the result:

By considering this fixed setting, the nominal spectral abscissa of the control feedback loop is read as follows: α((0.3,0.1)) = −1.4454; while the delay-free case is: α((0,0)) = 0.1323. This means that system stability is delay dependent, and some switching delays might be found for fixed controlled plant and controller parameters.

Discretization tools and techniques

This section intends to provide the reader with the basics of computational discretization techniques and implements useful for the course of the below proposed algorithm.

ACP polynomial approximation

The DDS algorithm introduced further on, in the results section, is designed to be practically implementable by means of most standard programs and utilities for scientific computing—without the necessity of the use of special software. Motivated by this, the first crucial step is to suggest a sufficiently accurate polynomial approximation of the CQP, resulting in the ACP, the leading zeros of which are close to those of the CQP, and can easily be computed.

The concept is based on the successive (iterative) discretization of the CQP, depending on the current leading pole estimation; i.e., a digital-filter-like implementation of the CQP is obtained. The Tustin (or bilinear) transformation is: and represents a standard technique for the corresponding designing of the digital filter to a continuous-time model, where T is the sampling period; and q expresses the shifting operator that agrees with z-1 in the z-transform [9]. Let us recall that the mapping given in Eq (9) stems from the Taylor series expansion of the mutual relation between the s-plane and the z-plane:

The crucial feature of this mapping is that it preserves stability, i.e. the poles in the complex left-half s-plane fall into the interior of the unit circle in the z-plane, and vice versa.

However, it is also known that—for finite-dimensional systems—the transformation given by Eq (9) does not preserve frequencies (i.e. so-called frequency warping occurs). The mutual relations are: where ω, ω stand for the frequency of the continuous-time model and the corresponding digital filter, respectively. Our intention is to make the ACP a digital implementation of the CQP, for which the imaginary part of the leading (dominant) root—after the use of the substitution as per Eq (10)—matches that of the CQP; which, i.a., means that both the corresponding switching imaginary poles might coincide. Hence, the CQP ought to be rescaled so that ω—calculated from Eq (11), equals the desired (original, non-scaled) ω. This idea yields the eventual pre-warping filter design mapping as follows:

The task of the choice of T is not unambiguous. With respect to delta models and derivative discretization, the value of T in Eqs (9) or (12) should be sufficiently small; while, on the contrary—the lower T is, the higher ACP degree is obtained. Note that for the z-transform, a recommendation for periodic systems has been given as: T = [0.2/ω0,0.5/ω0], where ω0 expresses the frequency of undumped oscillations [49]. Adopting the basic knowledge about a second order finite-dimensional model: ω0 = |s0|, where s0 means a pole from the conjugate pair, one can write the interval as:

Note that for infinite-dimensional systems; s0 can be taken as the dominant pole.

The specific feature of TDSs is that the CQP includes exponential terms that are not subjected to either Eqs (9) or (12). In fact, the bilinear transformation acts only on s-powers in the CQP corresponding to output derivatives in the system model. Exponentials are processed by natural shifting:

Switching delay estimation

Once the ACP is found for each grid node given by the delay discretization, i.e. for a particular value of τ, its leading zero is calculated—thus expressing the estimation of the leading zero of the CQP. As the gridding may be very sparse, the imaginary axis can usually be leaped over. A smoother switching delay estimation—and thus, that of the switching pole, can be made by the linear RF interpolation as follows: where is the estimation of the kth element of the switching delay vector; τ, τ stand for the delay elements in the jth and (j + 1)th nodes, respectively; and represents the zeros of the ACP in the corresponding nodes. To rephrase, Eq (15) approximates the zero of the function: . Let us recall that the necessary stability switching condition formulated in Eq (5) must hold for the eventual value of .

QPmR overview

The QPmR mapping algorithm was designed to compute all quasipolynomial zeros located in a given region in . Its leading idea is as follows: Let, s = β + jω, , and thus split the CQP into two real functions R(β,ω)≔Re{d(β + jω,)}, I(β,ω)≔Im{d(β + jω,)}. The region of interest is covered by a rectangular mesh grid: . Then, the zero of the CQP can be found by the intersection of the curves given by the solutions of: R(β,ω) = 0, I(β,ω) = 0, within the gridded region. The accuracy of the zeros is subsequently increased by the Newton method.

The original version (v.1) was designed by Vyhlídal and Zítek [41], and implemented in Matlab based on the use of the Symbolic Math Toolbox. However, the toolbox led to a relatively large computational time. The QPmR was improved in the consecutive work [42], where optimization via the reduction of the scanned regions was implemented. This was achieved by using the argument principle and spectrum distribution analysis [50]. As a result, a significant reduction of computational time was obtained; however, the complexity of the algorithm increased rapidly. The enhanced algorithm was implemented once again in Matlab, as an advanced QPmR (aQPmR). The up-to-date version (v.2) of the algorithm was published in [43]. The authors avoided the use of the Symbolic Math Toolbox—and included the recursive grid density adaptation, which resulted in at least twice as fast computing rates.

Naturally, the QPmR can be utilized for switching delay searches, either as a trial-and-error procedure when directly searching within the given delay range, or as a tool for the evaluation of CQP zeros, instead of the polynomial approximation introduced above. Note however, that the attention of this paper is, i.a., devoted to achieving the goal of avoiding the necessity of using special software tools. Our intention is to employ the QPmR as a benchmark for the algorithm presented herein.

Results and discussion

Prior to the presentation of the proposed DDS algorithm and its implementation on stability margin detection when remote controlling a skater on a swaying bow, some additional methods have to be derived and computational ideas formulated.

Integer-powers polynomial interpolation

The discrete-time formulation of exponentials given by Eq (14) generally results in: ; thus, it is necessary to apply an interpolation: in order to acquire a polynomial with integer powers. A possible way to cope with this task is presented in the following theorem.

Theorem 1: A term , where , 0 < m < 1 can be approximated in the neighborhood of arbitrary as:

Linear interpolation:

Quadratic interpolation: The proof of Theorem 1 is provided in the Appendix.

Corollary 1: If Eq (16) is performed at z0 = 1, which agrees with s0 = 0 in the s-plane, the natural and intuitive result: z−( ≈ (1 − m)z− + mz−( is obtained—which expresses the linear distribution in the accordance with the “closeness” of m to 0 or 1.

In our algorithm, the value z0 is selected as the leading (dominant, right-hand most) pole s0 subjected to Eq (10).

Successive leading pole estimation

Once the ACP in z-1 is found by using Eqs (9) or (12), (14), (16) or (17), its dominant root can be computed simply by means of standard engineering computing software, which is a benefit of the algorithm presented below. In order to obtain a more precise estimation, the iterative procedure rather than a single-step calculation is used. It initially adopts the result obtained from the preceding delay discretization grid node (or more precisely, from the nearest previously solved node)—providing the preliminary estimation of the particular dominant pole: . The choice of ω in Eq (12) and the time period in Eq 13 in its ith iteration is given by: and the interpolation introduced in Eqs (16) or (17) is performed at —resulting from by means of Eq (10), which gives rise to the ACP as where stand for the parameters′ values that determine the polynomial coefficients. In every single iteration step, the updated leading root of the ACP is subjected to Eq (10) so as to consequently update: . In fact, the root closest to the preceding one is evaluated as the leading one for the next iteration step. The iterative procedure for the particular node of the grid is terminated if for a selected ε > 0.

The motivation for the iterative calculation of the ACP in the neighborhood of the currently found leading pole can be explained as follows: The pole furthest right of the CQP spectrum is the only decisive one with respect to TDS stability—as clear from Proposition 1. Let this pole: s0, be isolated and known for a particular τ0 exactly. Due to Property 1, for any ε > 0, there exists δ > 0 such that if: τ1 = τ0 + Δτ,‖Δτ‖ < δ, the corresponding root s1≔{s → min|s −s0|:d(s, τ1) = 0} satisfies |s −s0| < ε. Hence, for a sufficiently small; roots s0 and s1 are close to each other. Let us assume the ACP as: and calculate its leading zero (i.e. that with the maximum modulus). It is supposed that resulting from via the exact mapping given in Eq (10) is located in the vicinity of s0 with . Now, by the iterative calculation of , the corresponding s-plane maps of their zeros, , should satisfy and hence, there exists a descending sequence: which converges for a suitable choice of T. A problem can emerge from Property 2 if for some δ > 0 and any Δτ due to the spectral abscissa discontinuity; however, such cases are rare and can be omitted in practice—since, for instance, the imaginary axis cannot be crossed by a multiple real pole [30,50]. Moreover, only strongly stable systems are taken into consideration.

Switching delay search algorithm

The proposed gridding switching delay search algorithm poses the main goal of this study and, it is summarized below by using ordered steps and a flowchart. The framework of the algorithm consists of the following steps: initialization, ACP computation and seeking its zeros, as well as switching delays and poles′ values enhancement.

Firstly, the mesh grid of delays within the selected ranges must be created, hence: τ = τ + Δτ, τ = 0, k = 1…L, j = 0…N − 1. The discretization step Δτ has to be chosen according to the minimum acceptable estimation error. Note, however, that numerical results have shown that this value may be approximately 103-times higher than eventual obtained accuracy. The counter of found switching delays (poles) are to be initialized as well, i = 1, and select the admissible convergence error while computing the ACP, ε > 0. The initial leading pole, , can easily be calculated exactly for τ = 0 as the rightmost zero of the polynomial d(s,0).

In each node of the grid (implemented as nested loops over τ), the iterative polynomial estimation of the CQP, i.e. the ACP, is computed. Once the final estimations of the dominant poles , in two subsequent grid nodes are found, see the preceding section for details, the test whether the imaginary axis has been crossed is executed. If it passes, the RF as in Eq (15) is successively performed to receive and updated values of , are computed by iterations for each k. The eventual delay vector has to be subjected to the RT condition test formulated in Eq (5). Estimations of switching delays (), switching poles () and their imaginary parts (), i.e. switching frequencies, are obtained as the algorithm outputs. The algorithm summary is as follows (Algorithm 1):

For a neutral TDS—verify the strong stability condition in Eq (2). If it does not pass, terminate the algorithm.

For the given d(s, τ), define the mesh grid: τ = τ + Δτ, τ = 0, k = 1…L, j = 0…N − 1 for a selected delay range, initialize the counter: i = 1 and choose ε > 0. Set: .

Compute:

For: (j1 = 0…N − 1, for (j2 = 0…N − 1… (for j = 0…N − 1 do steps 5–11)).

If: j1 = j2 = … = jL = 0, the inner loop is finished; if not, define: and set:

Compute: according to Eqs (9) or (12), (14), (16) or (17). Find its roots (z) and subject them to Eq (10), which yield s. Define:

While , set ; and then go to step 6.

Set: ; if , the inner loop is finished (see step 4); otherwise, i = i + 1 and go to step 9 (RF).

Set where function is defined in Eq (15) and perform the cycle given in step 10.

For: k = M − 1,…,1 do: If j = 0, set and go back to step 10; or else, set and and then compute the leading root from as in steps 6–7. Update values and find the leading root from according to the cycle in steps 6 and 7 again, and update the value . Then: . When the loop is finished, consolidate and go to step 11.

Update , . Compute the leading zero iteratively: from in accordance with steps 6–7. Set .

Outputs: .

A flowchart of Algorithm 1, beneficial particularly for practitioners, is depicted in Fig 2.

Fig 2

A Flowchart diagram of Algorithm 1.

The flowchart diagram displays the designed DDS algorithm. Its step-by-step formulation is presented in Algorithm 1. The ellipse expresses the starting point, and the oval means the termination point. The rhomb, the trapezium and the square stand for the condition, the output and the execution block, respectively.

Application of the algorithm to the model of a skater

The particular example presented herein, is aimed at Algorithm 1′s accuracy and speed verification, and demonstrates its usability and applicability in bio-cybernetics, by means of the estimation of the stability switching delays for the control loop of a model of a skater on a controlled swaying bow (see Eq (6) and Fig 1) and equipped with the controller given by Eq (7). A set of tests includes the determination of a suitable value of T, the comparison of the use of the linear and quadratic interpolations is according to Eqs (16) and (17) respectively, and the demonstration of the pre-warping asset.

Prior to the presentation of the course of Algorithm 1 itself, let us perform some numerical tests to determine a suitable value of T.

Sampling period determination

Let a particular delay vector be chosen and fixed as: τ = (0.07,0.07), and let us parameterize T as: with respect to the upper bound of the condition given by Eq (13). Let us assume moreover, that the leading zero of d(s,(0.07,0.06)) is known exactly—the use of the QPmR with the precision of 10−9 gives the value: s0 = 0.024805739 ± 3.905723903i—which is used as the initial estimation. An accuracy comparison of leading pole estimations by means of the ACP described in steps 6 and 7 of Algorithm 1 is shown in Table 1. Both the first and second order interpolations as in Eqs (16) and (17) respectively are included—and, moreover, the beneficial impact of pre-warping governed by the Eq (11) is demonstrated. The error expresses the distance of the approximated leading pole from s0 = 0.007985448 ± 3.868912827i; found by the QPmR with the above-mentioned precision. The admissible iteration error was set to: ε = 10−6 (see step 1 of Algorithm 1) and the maximum acceptable number of iteration steps was set to 200. The eventual approximated leading pole loci are given to the reader in S1 Table.

Table 1

Leading pole estimation accuracy benchmark.

Method	kT	Error	itermax
Eqs (9) and (16)	2	7.025e-2	5
	3	3.109e-2	4
	5	1.112e-2	4
	10	2.79e-3	3
	20	6.99e-4	5
	25	4.52e-4	29
	100	-	Inf
Eqs (12) and (16)	2	2.86e-4	6
	3	1.28e-4	4
	5	4.6e-5	4
	10	1.2e-5	3
	20	3e-6	3
	25	1.0e-5	2
	30	3.6e-5	18
	100	-	Inf
Eqs (12) and (17)	2	2.86e-4	5
	3	1.28e-4	4
	5	4.6e-5	4
	10	1.2e-5	3
	20	3e-6	6
	25	4e-6	5
	30	1.9e-5	37
	100	-	Inf

There are the results of a comparison of the use of Eqs (9), (12), (16) and (17) in Algorithm 1, and applied to the control of a skater on a swaying bow for a fixed delay value and several values of T–see Eq (13). Error means the distance of computed leading poles from the exact analytic one found by the QPmR. The value of itermax stands for the number of iterations within the successive polynomial leading root computation until the precision of ε is reached.

As can be seen from Table 1, the degree of precision is one to three orders better with the use of pre-warping than the bilinear transformation without it. Notice moreover, that it is necessary to perform a lower number of iterations to reach the desired ε for higher k—yet, the value is bounded.

Generally, high values of k require enormous number of iterations—or even cause divergence—and, the error thus increases slightly. The growth of the error value can be explained as follows. Algorithm 1 uses the bilinear (Tustin) transformation as in Eqs (9) or (12) to get the associated polynomial, and the (back) transformation from the s-plane to the z-plane via the exponential formula as in Eq (10). Geometrically, the conformal mapping (Eq (9)) maps the closed left-hand side of the complex s-plane inside the circle in the z-plane with the radius of 1/T, i.e., linearly with respect to k. The transformation according to Eq (10) provides the map inside the unit circle, i.e. the radius does not depend on the value of k. From the operator point of view, the Tustin formula represents the trapezoidal (feedback-feedforward) form of the derivative approximation. It is closely related to the delta transform. Hence, intuitively, the shorter the step T is (i.e. the higher k is), the approximation more precise is. However, the latter formula agrees with the transition from the z-transform to the s-transform (or vice-versa). It is well know that the z-transform fails for a very small values of the time period T (for instance, transfer function poles from the transfer function in s collapse to the point 1 in the z-plane). Moreover, the division be a very small number induces numerical problems when computing. The algorithm has to respect a compromise between both the principles. Due the properties of the z-transform and because of the numerical representation of real numbers inside the computer, the value of T must be sufficiently high, i.e. the value of k can not be excessively high.

Only a slight improvement is achieved by means of a quadratic interpolation—which is, however, sacrificed by the necessity of more iteration steps. In fact, two leading roots of the ACP do not constitute a complex conjugate pair due to complex polynomial coefficients. It should be noted that the quadratic interpolation yields better mutual convergence of these roots with ascending k.

The advantage of the iterative calculation of the ACP can be demonstrated for instance, by the evolution of the error —see steps 6 and 7 of Algorithm 1 for ; and s0 stands for the exact dominant system pole—with k = 15 and pre-warping as: {23.7,5.1,5}·10−6. Note that ω in Eq (12) has been set as: .

It seems from Table 1 that the suitable value of T lies approximately in the interval k ∈ [10,25]. Now, let us perform the second test as follows: Let the leading zero of d(s,(0.06,0)) be known exactly and perform Algorithm 1 in a degenerate form for τ1 ∈ [0.06,0.07], τ2 ∈ [0,0.2] with Δτ = 0.01. The numerical results are summarized in Table 2, where expresses the switching delay estimation, s0 is the corresponding leading system pole calculated by the QPmR, the error means the minimum distance of from the one found by the QPmR with the delay resolution of Δτ = 10−7; and the number of iterations (itermax) expresses the mean value rounded to the nearest integer.

Table 2

Switching delay estimation benchmark.

Method	kT	τ¯	s0	Error	itermax
Eqs (9) and (16)	2	(7.00526, 7.77128)e-2	-5.8741e-3 + 3.840208i	3.269e-3	4
	3	(7.00726, 7.58688)e-2	-2.5627e-3 + 3.846954i	1.427e-3	5
	5	(7.00751, 7.49427)e-2	-9.170e-4 + 3.846954i	5.09e-4	3
	10	(7.00746, 7.45718)e-2	-2.324e-4 + 3.851742i	1.27e-4	3
	20	(7.00744, 7.44771)e-2	-6.99e-5 + 3.852092i	3.3e-5	11
	25	(7.00763, 7.44654)e-2	-4.52e-5 + 3.852128i	2.2e-5	14
	100	-	-	-	Inf
Eqs (12) and (16)	2	(7.00754, 7.44444)e-2	-5.7e-6 + 3.8522095i	2.7e-6	5
	3	(7.00747, 7.44445)e-2	-5.6e-6 + 3.8522096i	2.2e-6	4
	5	(7.00744, 7.44454)e-2	-5.6e-6 + 3.8522096i	2.1e-6	4
	10	(7.00743, 7.44456)e-2	-5.6e-6 + 3.8522096i	2.1e-6	3
	20	(7.00743, 7.44454)e-2	-5.4e-5 + 3.8522100i	2.0e-6	5
	25	(7.00744, 7.44455)e-2	-5.7e-5 + 3.8522093i	2.2e-6	34
	30	-	-	-	Inf
Eqs (12) and (17)	2	(7.00754, 7.44444)e-2	-5.7e-6 + 3.8522095i	2.7e-6	5
	3	(7.00747, 7.44445)e-2	-5.6e-6 + 3.8522096i	2.2e-6	4
	5	(7.00744, 7.44454)e-2	-5.6e-6 + 3.8522096i	2.1e-6	4
	10	(7.00743, 7.44456)e-2	-5.6e-6 + 3.8522096i	2.1e-6	3
	20	(7.00744, 7.44454)e-2	-5.8e-5 + 3.8522090i	2.2e-6	4
	25	(7.00732, 7.44455)e-2	-3.4e-5 + 3.8522136i	1.4e-6	18
	30	-	-	-	Inf

Algorithm 1 is tested when the initial CHP root is known exactly for τ ∈ (0.06,0). Delay values go through intervals τ1 ∈ [0.06,0.07], τ2 ∈ [0,0.2] with Δτ = 0.01. The value of stands for the eventual estimated switching delay vector; and s0 is the corresponding leading ACP root. Error means the minimum distance of from the exact one. The number of iterations is rounded.

Again, it can be deduced from Table 2 that the use of pre-warping brings about a significant improvement of the switching delay estimation. The quadratic interpolation however, does not come up to expectations; nevertheless, it can reduce the number of iterations for the ACP in a specific range of k and results in a better mutual approach to the leading complex pair. In order to also include the time consumption factor—given mainly by the number of iterations, the eventual range of k for this example is approximately: k ∈ [3,15] for Eq (9) and k ∈ [3,20] for Eq (13)–which goes beyond the suggestion formulated in Eq (13) and stands somewhere between the z-transform and derivative discretization. Hence, after some additional numerical experiments for the switching delay search herein below, we finally decided to set a uniform k = 8. Let us recall that T is adapted according to Eq (23) in every single iteration step of the ACP calculation with respect to the leading pole estimation ().

Switching delay computation

As the final numerical experiment—the complete course of Algorithm 1 is performed and its results are presented here below. Let the particular region be selected as: R1 ≔ τ1 × τ2 ∈ [0,0.8]×[0,0.8], with Δτ⋅ = 0.01, i.e. N = 80, and ε = 10−6. The linear interpolation formula, according to Eq (16), is used. Entries of the found values are depicted in Fig 3 and compared with the stable/unstable region calculated by the QPmR of a rough delay resolution of Δτ⋅ = 0.01. Supporting files S2A, S2B and S2C Table display signs of the rightmost CQP roots from the QPmR, computed leading roots of the ACP along the stability border lines (with the use of pre-warping), and particular values of and , respectively.

Fig 3

Switching delays found by Algorithm 1.

Delay values computed by using pre-warping are indicated by circles and a joint via the solid line; the ones without pre-warping are represented by diamonds and joint by the dashed line. Results are back-grounded stability/instability regions computed by the QPmR with Δτ⋅ = 0.01.

During computation, we have observed that the convergence of the ACP estimation may be poor when the leading pair of poles found is sufficiently far from the imaginary axis; or, when α(τ) is near to its extreme (minimum). In such cases, the problem was solved by an interim small change in k. Both the plots in Fig 3 displaying results with/without the use of pre-warping, are almost indistinguishable by sight; hence, to provide the reader with results in more detail, select subsets R2 ≔ [0.05,0.1] × [0.05,0.1] and R3 ≔ [0.069,0.81] × [0.063,0.075] of R1 and display found switching delay pairs back-grounded by stable/unstable grid nodes computed by using the QPmR with the gridding Δτ⋅ = 0.005 and Δτ⋅ = 0.01, respectively, see Fig 4A and 4B.

Fig 4

Details of Fig 3 in selected regions.

Two regions were selected, a moderately detailed one R2 (A) and a very detailed one R3 (B). The stable and unstable areas found by the QPmR for Δτ⋅ = 0.005 and Δτ⋅ = 0.001, respectively, are indicated as well.

Table 3 includes the estimated switching delays, the corresponding leading roots found by the QPmR with the precision of 10−9, the error with the same meaning as in Table 2 and the RT vector within the region R2.

Table 3

Highlighted results of Algorithm 1—Inside region R2.

Mapping	τ¯	s0	Error	RT
Eq (9)	(5.00643, 9.68377)e-2	-3.578e-4 + 3.847336i	2.03e-4	(-1.980057–3.902248i,-1.776911–3.712380i)
	(6.00732, 8.57107)e-2	-3.592e-4 + 3.849544i	2.02e-4	(-1.979616–3.901380i,-1.784657–3.703965i)
	(7.00747, 7.46430)e-2	-3.604e-4 + 3.851478i	1.99e-4	(-1.979793–3.901079i,-1.793172–3.696423i)
	(8.00698, 6.36342)e-2	-3.609e-4 + 3.853141i	1.98e-4	(-1.980603–3.901336i,-1.802458–3.689744i)
	(9.00595, 5.26837)e-2	-3.611e-4 + 3.854537i	3.47e-4	(-1.982059–3.902144i,-1.812522–3.683919i)
Eq (12)	(5.00665, 9.66366)e-2	-4.9e-6 + 3.848074i	2.4e-6	(-1.977402–3.901955i,-1.774661–3.711958i)
	(6.00741, 8.55114)e-2	-5.5e-6 + 3.850279i	3.0e-6	(-1.976951–3.901107i,-1.782401–3.703575i)
	(7.00743, 7.44455)e-2	-5.6e-6 + 3.852210i	2.1e-6	(-1.977117–3.900824i,-1.790908–3.696065i)
	(8.00682, 6.34384)e-2	-5.1e-4 + 3.853870i	2.7e-6	(-1.977919–3.901097i,-1.800189–3.689416i)
	(9.00568, 5.24897)e-2	-4.3e-4 + 3.855262i	3.3e-6	(-1.979368–3.901922i,-1.810246–3.683621i)

The results of Algorithm 1 within region R2 = [0.05,0.1] × [0.05,0.1]; Root tendency (RT) was defined in Eq (3). The rest of the quantities have already been defined in Table 2.

As can be seen from Table 3, the results correspond to those in Table 2—and, once again, one can state that the use of pre-warping mapping gives two orders higher precision. The overall closeness to the stability bound computed by the QPmR is very good.

Fig 5A, 5B and 5C display feedback output responses for τ = (0.3,0.1), τ = (0.0700747,0.74643) and τ = (0.04,0.04) respectively, when there is the reference step change Δr = 5° = π/36rad from the zero steady state on the control feedback input at time t = 10s, to verify the result in the time domain. Apparently, delay values from the region determined as stable (namely, the nominal ones) yield a stable step response (5A); while the ones selected from Table 2 as switching delay result in a steady periodic output of a skater’s position angle deviation (see 5B), and the third chosen delays give unstable response (5C). Let us recall that—amazingly, smaller delay values deteriorate stability in this case. Note that S1 File includes simulation MATLAB® (Simulink) scheme of the control loop.

Fig 5

Time domain results verification.

The three plots of the system’s step reference time responses for three different delay vectors are displayed. Whereas the nominal setting τ = (0.3,0.1) (A) gives a stable response asymptotically tracing the reference signal, the second option τ = (0.0700747,0.74643) (B) indicates the stability border by steady oscillations, and the last one: τ = (0.04,0.04) (C), gives unstable oscillations with a rising amplitude.

Computational complexity

Although it is definitely not a difficult or crucial task in the case of the DDS algorithm introduced here, the (asymptotical) computational complexity should be discussed anyway. Parameterized statements are presented first, followed by practically measured data from the above numerical bio-cybernetic example. The basic general facts about Algorithm 1 are going to be presented, as well as the complexity of the QPmR.

The course of Algorithm 1 can be divided into two parts: a cyclic and a non-cyclic one. The non-cyclic part of the program implementation has the time consumption ϑ; consisting of the initialization and the CQP definition. The cyclic one has two subdivisions: with, or without, the RF, ϑ, ϑ, respectively. The former one mainly includes the iterative polynomial approximation and the leading zero calculation, whereas the latter one expresses time for the successive RF interpolation when the imaginary axis is crossed. The following statement can be formulated:

Proposition 2: The computational complexity of Algorithm 1 reads:

The proof of Proposition 2 is given in the Appendix.

As mentioned above, the QPmR can be used in Algorithm 1 instead of steps 6 and 7 to avoid the computation of the ACP; or within the direct procedure, in which for a sufficiently fine grid of discretized delays, the CQP zeros inside a defined region are computed. Regarding the latter case, its complexity can be expressed simply as: where N stands for the number of delay discretization steps for the QPmR; ϑ represents the time that the initialization of the QPmR takes outside the cycle, and ϑ means the time period of the QPmR course over a grid node itself. Note that the duration ϑ significantly depends on the selected region for computation of the zeros, its density, and the desired eventual accuracy.

Real computation time efficiency

Finally, let us provide the reader with the concise real computation-time requirements of Algorithm 1, compared to the use of the QPmR for direct leading pole computation inside a defined region. The particular variables and symbols were defined in Eqs (24) and (25). Computational tests were performed by means of a laptop equipped with 32b Windows® 7 Professional, running on an Intel® Core™2 DuoCPU P8700@2.53GHz, 4GB memory, in the Matlab® 7.11.0.584 (2010b) environment.

The whole course of the final computing test introduced above of Algorithm 1 took 2595 s, which gives the following mean durations of particular operations: ϑ = 0.54s, ϑ = 0.36s, ϑ = 1.64s for found c1 = 1.75; hence, it can be calculated that c2 = 2.66.

Now, let us consider the QPmR and two possible search rectangles in the complex plane and their corresponding computed precisions. The first is chosen as a region of the size with its center in , which expresses the currently found leading pole estimation (in the ith iteration). For such a case, we have measured the mean values ϑ = 0.66s and ϑ = 0.62s (for the selected precision of 10−8). The second rectangle has the fixed size of [1×1] centered in ; again, for which we have obtained ϑ = 0.68s and ϑ = 0.47s (for the reduced precision of 10−7). These results imply that the use of the QPmR for seeking the leading pole instead of the ACP roots computing introduced in Algorithm 1 does not yield a shorter computation time. Moreover, the direct trial-and-error use of the QPmR is hardly applicable; for instance, if the latter region is considered with the discretization step of Δτ⋅ = 0.001 (for a sufficiently smooth result), the overall computation time would approximately be 3 · 105 s; see Eq (25).

Conclusions

The formulation of a novel gridding multiple stability switching delay search algorithm and its verification by a numerical bio-cybernetic example of a human being on a remotely controlled swaying bow were the main objectives of this contribution. The proposed DDS algorithm can be fitted into a group of frequency-domain direct methods, based on an effort to find all characteristic roots (poles) located on the stability border, i.e. on the imaginary axis. Within these methods, it is usually only the upper and lower bounds on the imaginary parts of such poles that can be found by the elimination of delays from the characteristic quasipolynomial [28,33,34]. If the delays are commensurate, a procedure can be relatively fast; however, in the general case of non-commensurate delays, a computationally-lengthy iterative reduction procedure ought to be performed [31]. The algorithm presented herein can deal with non-commensurate delays by more effectively omitting a complex mathematical apparatus. The main advantages of the algorithm reside in the fact that it utilizes basic mathematical operations and standard software tools solely; thus, it is simple to implement by engineers and practitioners.

The effectiveness, simplicity, rapidity and accuracy of the gridding algorithm have been tested by the numerical example. The simple control feedback loop with delays system—applied to a skater on a swaying bow, is controlled by a conventional controller. For such a configuration, the marginal stability values of the delay vector were found by the proposed DDS algorithm; these results can be utilized practically, e.g. in biomechanics. The example presented here, proved that the procedure is computationally more effective than the direct use of the QPmR algorithm [43].

The formulation of the Algorithm 1 introduced above, represents its basic version. There are some natural possibilities regarding how to suggest a modification—for instance, the delays′ discretization can be performed with diverse settings for some k1 ≠ k2 rather than a constant N, or with a non-constant Δτ for various k,j. In addition, it is not necessary to take τ = 0; however, in such a case, a sufficiently accurate estimation of the primal pole has to be made, e.g. by means of the QPmR. These topics might be solved in future research. Its applicability is limited to retarded and strongly stable neutral delay systems—which also yields further natural possible research directions. In addition, the numerical example has shown that the crucial step consists in the determination of a suitable sampling period that influences the accuracy and the convergence rate.

Appendix

Proof of Theorem 1: Assume the linear interpolation first. Then, this can be written: where the left-hand factor has the integer exponent 1 − m ∈ (0,1). Let the intention be to find: such that z1− ≈ a0 + a1z. Now employ the Taylor series expansion at any z0: for some emerging from the remainder. If |z − z0| is sufficiently small, the last term can be neglected. Then, by substituting Eqs (27) into (26) and applying simple algebra on Eq 26 the eventual formula—as in Eq (16) is obtained.

Now consider a proof of Eq (17) and use an analogous idea. Hence: where 2 − m ∈ (1,2). The Taylor series expansion of the latter term can be given as: for some . Again, the cubic term can be cancelled from Eq (29) for some z sufficiently close to z0, and by combining Eqs (29) and (28)—these yield Eq (17), finally.

Proof of Proposition 2: The cyclic part must act in every grid node, the overall number of which is (N +1). Only if the imaginary axis is crossed by the leading pole, is the RT part applied in addition. This happens approximately once or a few times (c1 < N) during the inner-most loop course—and this loop is entered (N +1) times.

The right-hand side of Eq (24) can be obtained by the supposition that the polynomial approximation is the most time consumptive operation. Whereas, outside the RF, there is only one approximation—the number of 2(L − 1) + 1 = 2L − 1 such operations inside the RF is performed, i.e. where assumes other low-time consuming operations within the cycle. Substituting this and into the left-hand side of Eq (24), the right-hand value is obtained. Note that for a particular N, L one can also write C(N, L) = ϑ + c3(2L + N)(N + 1)ϑ, c3 > 0.

The eventual approximated leading pole loci for Table 1.

Leading pole estimations of the ACP received from steps 6 and 7 of Algorithm 1, which gives rise to errors introduced in Table 1.

(XLSX)

Click here for additional data file.

Results of the numerical experiment according to Algorithm 1.

Signs of the rightmost quasipolynomial roots from the QPmR for Fig 3 are given in part A: The data table includes signs of the CQP roots calculated via the QPmR that give rise to the stable and unstable regions in Fig 3. Computed leading roots of the ACP along the stability border line of Fig 3 are provided in part B: The table displays the computed leading roots of the ACP for delay values in the vicinity of the stability border and also their real parts for some other selected delays, as in Fig 3. The pre-warping strategy is used here. Computed switching delays and poles for Fig 3 are introduced in part C: In the table, the reader can find values sets of and for both the strategies (with or without pre-warping) to support data in Fig 3.

(XLSX)

Click here for additional data file.

A simulation scheme of the skater stabilization.

This *.mdl file from MATLAB® enables to test the feedback response for various τ1 and τ2 that can be set inside the "Skater model" block.

(MDL)

Click here for additional data file.