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][2][3][4]. It is, i.a., given by the fact that delay appears PLOS ONE | https://doi.org/10.1371/journal.pone.0178950 June 8, 2017 1 / 23 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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 1 ) 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 unstablewhich 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][42][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, C, R, N, N 0 denote the set of, respectively, complex numbers, real numbers, positive and non-negative integers. The n-dimensional Euclidean space is denoted as R n and let R n 0þ :¼ fa ¼ ða 1 ; a 2 ; . . . ; a n Þ 2 R n ja i ! 0; i ¼ 1; 2; . . . ; ng. For s 2 C, 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 . . . ; t L Þ 2 R L þ ; L > 1, represent independent delays: l ij;k 2 N 0 and x ij 2 R. Let x a ðs; τÞ ¼ 1 þ X h n j¼1 x nj exp À s X L k¼1 l nj;k t k be the associated exponential polynomial related to x(s, τ). If, d a ðs; τÞ 2 R, 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 i agree with the zeros of d(s, τ), i.e. s i ≔ {s:d(s, τ) = 0}, and thus: CQP d(s, τ). Let us denote the spectrum of the system as S ≔{s i }, the set of all zeros of d a (s,) as S a 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: 1. If nonzero d ij , λ ij,k exist (where d ij stands for the real-valued coefficient in d(s,) for s i with jth delay) for some positive τ k and some i, j, k, then the number of poles is infinite.
2. For a retarded system, for any b 2 R with |β| < 1, only a finite number of poles are located in the half-plane Re(s) > β; while infinitely many others are located in Re(s) < β.
3. For a neutral system, a vertical chain of poles (s i ) exists at γ ≔ max ReS a such that lim i!1 Re(s i ) = γ, lim i!1 Im(s i ) = 1 for |s i | < |s i+1 | js i j < js iþ1 j.

Let g :¼ sup
Dτ!0 g, then the number of poles with Re s i > g is always finite, and they are isolated.
5. Isolated poles behave continuously and smoothly with respect to τ on C.
1. If one defines the spectral abscissa as the function: α(τ)≔τ 7 ! sup Re S, 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].
2. 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.
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. g < 0Þ. 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 τ c,i , and the corresponding crossing frequencies, ω c,i , satisfy: For switching delays τ i , and the corresponding switching frequencies o i , the condition in Eq (4) holds-and moreover, the corresponding pole (s i ¼ jo i ) on the imaginary axis is the leading (furthest right) one; i.e. aðτ i Þ ¼ Re s i ¼ 0, 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: T and O, 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 frequencydomain 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.
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 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 axiswhich is then taken as the system output. https://doi.org/10.1371/journal.pone.0178950.g001 Stability switching search algorithm: A human on a swaying bow 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 ω, ω d 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 ω dcalculated 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 = |s 0 |, where s 0 means a pole from the conjugate pair, one can write the interval as: Note that for infinite-dimensional systems; s 0 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 t k is the estimation of the kth element of the switching delay vector; τ k,j , τ k,j+1 stand for the delay elements in the jth and (j + 1)th nodes, respectively; andŝ j ;ŝ jþ1 represents the zeros of the ACP in the corresponding nodes. To rephrase, Eq (15) approximates the zero of the function: Reŝðτ k Þ. 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 C. Its leading idea is as follows: Let, s = β + jω, b; o 2 R, 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: W=T = 2 N 0 ; 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 z À ðnþmÞ 2 C, where n 2 N 0 , 0 < m < 1 can be approximated in the neighborhood of arbitrary z 0 2 C as: 1. Linear interpolation: 2. Quadratic interpolation: The proof of Theorem 1 is provided in the Appendix. (16) is performed at z 0 = 1, which agrees with s 0 = 0 in the s-plane, the natural and intuitive result: z −(n+m) % (1 − m)z −n + mz −(n+1) is obtained-which expresses the linear distribution in the accordance with the "closeness" of m to 0 or 1.

Corollary 1: If Eq
In our algorithm, the value z 0 is selected as the leading (dominant, right-hand most) pole s 0 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: s 0;ð1Þ . The choice of ω in Eq (12) and the time period in Eq 13 in its ith iteration is given by:ŝ 0;ðiÞ and the interpolation introduced in Eqs ( 16) or (17) is performed atẑ 0;ðiÞ -resulting fromŝ 0;ðiÞ by means of Eq (10), which gives rise to the ACP asdðz À 1 jτ;ŝ 0;ðiÞ Þ ¼ X i¼0nd i ðτ;ŝ 0;ðiÞ Þz À i where τ;ŝ 0;ðiÞ stand for the parameters 0 values that determine the polynomial coefficients. In every single iteration step, the updated leading rootẑ 0;ðiþ1Þ of the ACP is subjected to Eq (10) so as to consequently update:ŝ 0;ðiþ1Þ . 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 jŝ 0;ðiþ1Þ Àŝ 0;ðiÞ j < ε 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: s 0 , be isolated and known for a particular τ 0 exactly. Due to Property 1, for any ε > 0, there exists δ > 0 such that if: τ 1 = τ 0 + Δτ,kΔτk < δ, the corresponding root s 1 ≔{s ! min|s −s 0 |:d(s, τ 1 ) = 0} satisfies |s 1 −s 0 | < ε. Hence, for a sufficiently small; roots s 0 and s 1 are close to each other. Let us assume the ACP as:dðz À 1 jτ 1 ; s 0 Þ and calculate its leading zeroẑ 1;ð1Þ (i.e. that with the maximum modulus). It is supposed thatŝ 1;ð1Þ resulting fromẑ 1;ð1Þ via the exact mapping given in Eq (10) is located in the vicinity of s 0 with js 1 À s 0 j > js 1 Àŝ 1;ð1Þ j. Now, by the iterative calculation ofd ðiÞ ðz À 1 jτ 1 ; s 1;ðiÞ Þ, the corresponding s-plane maps of their zeros,ŝ 1;ðiþ1Þ , should satisfy js 1 Àŝ 1;ðiÞ j > js 1 Àŝ 1;ðiþ1Þ j and hence, there exists a descending sequence: fjd ðiÞ ðŝ 1;ðiÞ ; τ 1 Þjg i max i¼1 which converges for a suitable choice of T. A problem can emerge from Property 2 if jŝ 1;ðiþ1Þ Àŝ 1;ðiÞ j > d 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 0 values enhancement.
Firstly, the L Â N 2 N 2 mesh grid of delays within the selected ranges must be created, hence: τ k,j+1 = τ k,j + Δτ k,j , τ k,0 = 0, k = 1. . .L, j = 0. . .N − 1. The discretization step Δτ k,j has to be chosen according to the minimum acceptable estimation error. Note, however, that numerical results have shown that this value may be approximately 10 3 -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,ŝ 0 ¼ s 0 , 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 τ k ), the iterative polynomial estimation of the CQP, i.e. the ACP, is computed. Once the final estimations of the dominant polesŝ 0 ,ŝ 1 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 t k and updated values ofŝ 0 ,ŝ 1 are computed by iterations for each k. The eventual delay vectort has to be subjected to the RT condition test formulated in Eq (5). Estimations of switching delays (T), switching poles (Ŝ) and their imaginary parts (Ô), i.e. switching frequencies, are obtained as the algorithm outputs. The algorithm summary is as follows (Algorithm 1): 1. For a neutral TDS-verify the strong stability condition in Eq (2). If it does not pass, terminate the algorithm.

Application of the algorithm to the model of a skater
The particular example presented herein, is aimed at Algorithm 1 0 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: s 0 = 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 s 0 = 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.

-Inf
There are the results of a comparison of the use of Eqs (9), (12), (16) and (17)  Stability switching search algorithm: A human on a swaying bow 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 T -yet, the value is bounded.
Generally, high values of k T 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 T . 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 T . 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 T 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 zplane). 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 T 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 T .
The advantage of the iterative calculation of the ACP can be demonstrated for instance, by the evolution of the error jŝ 1 À s 0 j-see steps 6 and 7 of Algorithm 1 forŝ 1 ; and s 0 stands for the exact dominant system pole-with k T = 15 and pre-warping as: {23.7,5.1,5}Á10 −6 . Note that ω in Eq (12) has been set as: Imŝ 1 .
It seems from Table 1 that the suitable value of T lies approximately in the interval k T 2 [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 2 [0.06,0.07], τ 2 2 [0,0.2] with Δτ = 0.01. The numerical results are summarized in Table 2, where τ expresses the switching delay estimation, s 0 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 (iter max ) expresses the mean value rounded to the nearest integer.
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 T 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 T for this example is approximately: k T 2 [3,15] for Eq (9) and k T 2 [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 T = 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 (ŝ 0 ).
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.
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 T . Both the plots in   Stability switching search algorithm: A human on a swaying bow 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 R 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 1 (Simulink) scheme of the control loop.

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 ϑ nc ; consisting of the initialization and the CQP definition. The cyclic one has two subdivisions: with, or without, the RF, ϑ c,nRF , ϑ c,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 QPmR stands for the number of delay discretization steps for the QPmR; ϑ QPmR,nc represents the time that the initialization of the QPmR takes outside the cycle, and ϑ QPmR,c means the time period of the QPmR course over a grid node itself. Note that the duration ϑ QPmR,c 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 1 7 Professional, running on an Intel 1 Core™2 DuoCPU P8700@2.53GHz, 4GB memory, in the Matlab 1 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: ϑ nc = 0.54s, ϑ c,nRF = 0.36s, ϑ c,RF = 1.64s for found c 1 = 1.75; hence, it can be calculated that c 2 = 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 ½10Refðŝ 0;ðiÞ Àŝ 0;ðiÀ 1Þ Þg Â 10Imfŝ 0;ðiÞ Àŝ 0;ðiÀ 1Þ g with its center inŝ 0;ðiÞ , which expresses the currently found leading pole estimation (in the ith iteration). For such a case, we have measured the mean values ϑ QPmR,nc = 0.66s and ϑ QPmR,nc = 0.62s (for the selected precision of 10 −8 ). The second rectangle has the fixed size of [1×1] centered inŝ 0;ðiÞ ; again, for which we have obtained ϑ QPmR,nc = 0.68s and ϑ QPmR,nc = 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 Á 10 5 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 systemapplied 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 0 discretization can be performed with diverse settings N k 1 6 ¼ N k 2 for some k 1 6 ¼ k 2 rather than a constant N, or with a non-constant Δτ k,j for various k,j. In addition, it is not necessary to take τ k,0 = 0; however, in such a case, a sufficiently accurate estimation of the primal polê s 0;...;0 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 systemswhich 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.
Supporting information S1 Table. 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