{'fraction_non_alphanumeric': 0.0399863391866124, 'fraction_numerical': 0.012977772718217264, 'mean_word_length': 4.220323874610013, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Root Locus for SISO Dead-Time Systems: A Continuation Based Approach

We present a numerical method to plot the root locus of Single-Input-Single-Output (SISO) dead-time systems with respect to the controller gain or the system delay. We compute the trajectories of characteristic roots of the closed-loop system on a prescribed complex right half-plane. We calculate the starting, branch and boundary crossing roots of root-locus branches inside the region. We compute the root locus of each characteristic root based on a predictor-corrector type continuation method. To avoid the high sensitivity of roots with respect to the locus parameter in the neighborhood of branch points, the continuation method relies on a natural parameterization of the root-locus trajectory in terms of a distance in the (characteristic root, locus parameter)-space. The method is numerically stable for high order SISO dead-time systems.


Introduction
The root locus method is an essential tool in modern control engineering for analysis and synthesis problems [15]. This method is successfully implemented for finite dimensional SISO systems and becomes a fundamental tool in control education [6,10].
The closed-loop of the SISO system with a time delay generically has infinitely many poles in the complex plane (see, e.g., [13]). Therefore obtaining the root locus for dead-time systems is a difficult problem. Unlike the finite dimensional case, the root-locus equation contains a time delay term and standard polynomial root-finding algorithms for the root computation are not applicable.
The root locus of SISO dead-time systems with respect to the controller gain is based on two main approaches: methods based on sweeping or gridding parameters on the complex plane and continuation based methods.
The methods in the first group are early methods in the literature. The root locus is obtained on a set of vertical lines [9], horizontal lines [18] or on a rectangular grid [11] in the complex plane by finding points satisfying the root-locus equation. These methods require a large number of grid points for an accurate root-locus plot and may miss the dynamics of some characteristic roots due to the finite number of grid or sweep parameters.
The methods in the second group are continuation based methods. These methods follow characteristic roots of the root locus by a predictor-corrector algorithm inside the root-locus region and detect the ones entering into the region due to the asymptotic root chains of time delay systems. The predicted characteristic root is obtained by computing the slope of the phase equation of the root-locus equation [2], the solution of nonlinear differential equations [17] or by a triangulation method on the complex plane [14] and the predicted values are corrected by a Newton-Raphson iteration. The characteristic roots entering into the root-locus region are detected by checking the sign of function values on the constant grid points of the region's boundary [2] or computing another root locus problem [14,17]. The approach proposed in [2] has two main limitations. First, the root locus is computed on a rectangular region and this is restrictive to observe the overall dynamics of characteristic roots. Second, it requires exhaustive search on a grid for characteristic roots entering into the root-locus region, with a possibility of missing a characteristic root depending on the number of grid points. The articles [14,17] consider general SISO time delay systems, i.e., SISO systems with state delays. The predictor step in both methods, the solution of nonlinear differential equations and the triangulation method, are numerically expensive com-pared to the simple linear predictor. In order to detect asymptotic roots another root locus problem has to be solved and the number of required roots for a complex region is not given.
The root locus of SISO dead-time systems with respect to the time delay allows us to analyze the effects of time delay on system stability and performances and it is less considered in the literature. A computationally expensive brute force approach is to compute the locations of the characteristic roots inside a desired region for a grid over the parameter space using, e.g., a spectral method [4]. A continuation based method in [7] obtains the root locus on the complex right half-plane by solving a nonlinear differential equation and correcting with Newton's method. The predictor step is computationally expensive and the method does not allow to analyze the root locus of stable closed-loop systems since the boundary of the root-locus region is the imaginary axis.
All continuation based methods for the root locus of SISO dead-time systems in the existing literature parameterize root-locus trajectories with respect to the controller gain or the time delay. This parameterization is numerically ill-posed due to the high sensitivity in the neighborhood of intersection points characterized by the presence of multiple roots, as reported in [7]. Both the predictor and corrector steps in all methods require evaluations of functions with time delays. These evaluations are numerically expensive and not numerically stable due to the oscillation and exponential increase of time delay terms in the imaginary axis and the positive real axis directions in the complex plane respectively.
In this article, we compute the root locus of SISO deadtime systems on a given complex right half-plane up to a predefined controller gain or a time delay. Our main contributions are the following: • We compute the trajectories of all characteristic roots entering into the prescribed region with respect to the controller gain or the time delay up to an upper bound value. By choosing the controller gain or time delay sufficiently large, the asymptotic behavior of the roots can be observed. • Our continuation method estimates the next root on the complex plane by a simple linear predictor and corrects this prediction with Newton's method. We parameterize the root-locus trajectories in the (root and gain/delay)-space in terms of the arclength which is numerically robust, as we shall see. • We use an adaptive step size in the prediction step depending on the convergence rate of Newton's method and the distance of the root from the root-locus trajectory. This makes our algorithm scalable since it uses different step sizes for different root-locus trajectories and adapts to the curvature of the trajectories. • We transform the root-locus equation into phase and logarithmic magnitude equations and avoid the oscil-lation and exponential increase problems due to the time delay term in function evaluations, which are only needed in the correction step.
The paper is organized as follows. Section 2 formulates the root-locus problem with respect to the controller gain and the time delay. The critical points of the root locus are computed in Section 3. The predictor-corrector based continuation method is given in Section 4. The overall algorithm for the root locus is presented in Section 5. Section 6 is devoted to numerical examples. In Section 7, some concluding remarks are presented. Notation: C, R, Z : sets of complex, real and integer numbers, ℜ(u) : real part of a complex number u, ℑ(u) : imaginary part of a complex number u, |u|, ∠u : magnitude, phase of a complex number u, ⌈u⌉ : the smallest integer larger than or equal to u, ⌊u⌋ : the largest integer smaller than or equal to u, u T : the transpose of the vector u, A + : the pseudoinverse of the matrix A, sgn(u) : returns +1, −1, 0 given a real number u for u > 0, u < 0, u = 0 respectively, N n2 n1 : the set of integers from n 1 to n 2 . mod(u, v): returns the remainder on division of the dividend u by the divisor v.

Problem Formulation
A SISO dead-time system is a rational, proper SISO plant G with a constant input or output time delay h. The only required input data to compute the root locus of these systems with respect to the controller gain or the time delay are • the poles of G, p i ∈ C for i ∈ N n 1 , • the zeros of G, z r ∈ C for r ∈ N m 1 , • the gain of G, α ∈ R, • the time delay h, h ∈ R and h > 0.
Note that the input data can be obtained from statespace matrices or the transfer function of G.
Let G(s)e −hs be the transfer function representation of the SISO dead-time system where and the real and imaginary parts of system zeros and poles are defined as z r = σ zr + jω zr for r ∈ N m 1 and p i = σ pi + jω pi for i ∈ N n 1 .
The root-locus equation of a SISO dead-time system with respect to the controller gain or the time delay is f (s, λ) = 0 where f (s, λ) = 1 + λG(s)e −hs (2) 1 + G(s)e −λs (3) and λ is the locus parameter for λ ∈ [0, λ max ] where λ max is a given positive number. The root-locus region is a complex right half-plane where σ 0 is a negative real number. The corresponding boundary of the root-locus region is a vertical line through s = σ 0 and parallel to the imaginary axis.
We define root-locus problems with respect to the controller gain or the time delay as follows.
The root-locus trajectory of the characteristic root is a curve on the complex plane on which each point satisfies the root-locus equation for λ ∈ [0, λ max ].The intersection of two or more root-locus trajectories is a branch point. The starting points of the root locus are the characteristic roots inside C σ0 for λ = 0 and the boundary crossing roots are the characteristic roots crossing ℜ(s) = σ 0 for any non-negative λ value smaller then λ max . The critical points are the starting points, branch points and boundary crossing roots of the root locus.
For the root locus w.r.t. the controller gain, we can compute branch points inside C σ0 apriori as we shall see. For the root locus w.r.t. the time delay, we detect the branch points while following the root-locus trajectory. The behavior of a root-locus trajectory around a branch point is given in the following Lemma [17].
Lemma 1 Ifs is a root of the root-locus equation in (2) or (3) with a locus parameterλ satisfying then N root-locus trajectories intersect at s =s. Intersecting trajectories continue straight after the intersection if N is odd and make an angle of − π N angle if N is even.
When a single root s of the root-locus equation crosses the boundary ℜ(s) = σ 0 , we determine whether it enters into or leaves the region C σ0 by computing its crossing direction [13] A root on the boundary enters into or leaves the region C σ0 depending on whether CD(s, λ) > 0 or CD(s, λ) < 0 respectively.
The computation of root-locus trajectories involves two main tasks: computing the critical points and following the root-locus trajectories. The next section focuses on the computation of the critical points of the root locus w.r.t. the controller gain and the time delay. We follow each root trajectory using a continuation method described in Section 4.

Computation of Critical Points
We compute the critical points of two root-locus problems in the next two subsections.

Controller gain as locus parameter
The root-locus equation w.r.t. the controller gain is given in (2). The starting points of the root locus are the poles of G inside the root-locus region C σ0 .
The branch points satisfy the root-locus equation and Thus the branch points are the zeros of the transfer function G ′ (s) − G(s)h inside the region C σ0 . This transfer function can be written as and realized in a state-space representation as a series connection of m + n first-order systems using the poles and zeros of G and the time delay h. Its zeros can be computed by solving an eigenvalue problem constructed from state-space matrices of G. Consequently, the branch points can be determined accurately.
The computation of boundary crossing roots of the rootlocus equation in (2) is a difficult problem. These roots and their crossing directions are computed in the following section.

Computation of boundary crossing roots
A root s on the boundary of the root-locus region ℜ(s) = σ 0 for the controller gain λ satisfies the magnitude and phase equations of the root-locus equation in (2). We first find the intervals on the boundary where the magnitude condition holds for some λ ∈ [0, λ max ]. This is equivalent to finding the intervals on ℜ(s) = σ 0 and ω ∈ [0, ∞) such that Lemma 2 Assume that G has neither poles nor zeros on the boundary of the root-locus region. The functions ω → Λ(ω) and ω → Λ ′ (ω) are continuous and the nonnegative zeros of Λ ′ (ω) are the non-negative real roots of the polynomial Proof. Using the transfer function of G in (1), the function Λ(ω) can be written as ln γ zr (ω) .
(9) The first derivative of the function Λ(ω) is The functions Λ(ω) and Λ ′ (ω) are continuous except the points where γ zr (ω) or γ pi (ω) are equal to zero. These points are the poles or zeros of G on ℜ(s) = σ 0 . The continuity results in Lemma 2 follow from the assumption. The polynomial given in (8) is the numerator of the function Λ ′ (ω) as in equation (10) and the result follows. ✷

Corollary 3
Assume that G has neither poles nor zeros on the boundary of the root-locus region. Then the function ω → Λ(ω) is monotonic on the intervals whose boundary points are consecutive non-negative zeros of Λ ′ (ω), 0 and ∞.
Proof. By Lemma 2, the function Λ(ω) is continuous since σ 0 is chosen such that there are neither poles nor zeros of G on ℜ(s) = σ 0 . Therefore it is monotonic inside the intervals determined by its extremum points and the end points of the boundary of the root-locus region, 0 and ∞. ✷ Since Λ(ω) is monotonic on each interval defined in Corollary 3, we find the subinterval in each interval where Λ(ω) satisfies the inequality in (7). This is done as follows. If both values of Λ(ω) at the interval end points are smaller than ln λ max , then all Λ(ω) values in this interval are smaller than ln λ max because Λ(ω) is monotonic. If one of the values of Λ(ω) at the interval end points is larger and the other one is smaller than ln λ max , we can find the point where Λ(ω) is equal to ln λ max by root-finding algorithms for a bracketed root of monotonic continuous function (such as Brent's method [5]) and take the subinterval satisfying the inequality in (7). If both values of Λ(ω) at the interval end points are larger than ln λ max , we discard that interval since all values of Λ(ω) are larger than ln λ max and the inequality in (7) never holds.
Based on this approach, we can compute the set of intervals I on the boundary of the root-locus region where the magnitude condition (7) is satisfied for some values The boundary crossing roots also satisfy the phase equation of the root-locus equation in (2) (2l over the intervals I on ℜ(s) = σ 0 . Here the function φ(ω) represents the continuous extension of the phase of the transfer function G(s)e −hs . It satisfies the following equation, where ω ∈ [0, ∞). The left-hand side of equation (11) represents constant functions of ω. If we partition the intervals I into the subintervals such that the function φ(ω) is monotonic on each subinterval, we can compute the boundary crossing roots by root-finding algorithms for a bracketed root of monotonic continuous function [5]. The following results allow us to compute the intervals on ℜ(s) = σ 0 where the function φ(ω) is monotonic.
Lemma 4 Assume that G has neither poles nor zeros on the boundary of the root-locus region. Then the functions ω → φ(ω) and ω → φ ′ (ω) are continuous and the nonnegative zeros of φ ′ (ω) are the non-negative real roots of the polynomial where the polynomials Γ r z for r ∈ N m 1 , Γ i p for i ∈ N n 1 , Γ z and Γ p are as defined in Lemma 2.
The first derivative of the function φ(ω) is Following the same arguments in Lemma 2, the functions φ(ω) and φ ′ (ω) are continuous by the assumption. The polynomial given in (12) is the numerator of the function φ ′ (ω) as in equation (15) and the result follows. ✷

Corollary 5
Assume that G has neither poles nor zeros on the boundary of the root-locus region. Denote by I φ the set of intervals whose end points are consecutive nonnegative zeros of φ ′ (ω), 0 and ∞. Then the function ω → φ(ω) is monotonic on each interval in the set I φ .
Proof. The function φ(ω) is continuous. The monotonicity of φ(ω) changes only at the points where φ ′ (ω) = 0. The assertion follows. ✷ The non-negative real roots of the polynomials (8) and (12) are zeros of the rational functions (10) and (15). These rational functions can be realized in a state-space representation as a series connection of second-order systems using the poles and zeros of G and the time delay h and the corresponding zeros can be computed by solving a generalized eigenvalue problem.
The intersection of two sets of intervals I and I φ partitions I into the subintervals as I = ∪ nI i=1 I i where φ(ω) is monotonic on each interval I i . Each intersection of the function φ(ω) and the constant functions in the left-hand side of equation (11) over the intervals I corresponds to a boundary crossing root since any such point on I satisfies both the magnitude condition in (7) and the phase equation in (11) of the root-locus equation on ℜ(s) = σ 0 . Since the function φ(ω) is monotonic on I i , we can compute each intersection point by root-finding algorithms for a bracketed root of φ(ω) over the interval I i . The value of the ω at the intersection point is the imaginary part of the boundary crossing root on the interval I i and the corresponding controller gain is the value of e Λ(ω) for this point. If there is no horizontal line intersecting φ(ω) on I i , we discard this interval since there is no root crossing. Based on the above explanation, the following algorithm computes the boundary crossing roots of the root-locus equation in (2).

Algorithm 1 For each interval in
(1) Compute the maximum and minimum values of the function φ over the interval I i . Since the function is monotonic over the interval, they are the maximum and minimum values of the function φ at the interval ends points, i.e., φ max Consider the integers l for which the constant function y = (2l + 1)π has an intersection with the function y = φ(ω) over the interval I i . Note that the inte- • compute the intersection point between the horizontal line (2l + 1)π and the function φ over the interval I i . This point can be considered as a the bracketed root of a monotonic continuous function φ over I i and accurately computed by root-finding algorithms such as Brent's method [5]. Denote this point by ω cr which is equal to the imaginary part of the boundary crossing root. • compute the corresponding controller gain for the boundary crossing root, i.e., λ cr = e Λ(ωcr ) .
By Algorithm 1, we compute all boundary crossing roots of the root-locus equation in (2) and their corresponding controller gain values for λ ∈ [0, λ max ]. Their crossing directions are determined by the following theorem.

Theorem 6
The crossing direction of a boundary crossing root s cr = σ 0 + jω cr only depends on the imaginary part ω cr on the boundary and is equal to CD(s cr , λ cr ) = −sgn (φ ′ (ω cr )) .
Proof. Using the transfer function G in (1) and (10,15), we obtain By the crossing direction formula in (5) and the equation in (17), the crossing direction of s cr at λ = λ cr is equal to . ✷ By Theorem 6, the crossing directions of the roots are the same when their imaginary parts are inside the same interval of I φ (see Corollary 5). We determine the crossing directions of boundary crossing roots from their imaginary parts. We group them according to their crossing directions and define the sets 1 are the boundary crossing roots entering into or leaving the root-locus region and their controller gains respectively. (5) is well-posed (either +1 or −1) if there are neither poles nor zeros of G or branch points on the boundary of the root-locus region. Remark 2: The closed-loop system is of neutral type (see, e.g., [13]) when G is biproper (i.e., d := G(∞) = 0). Its asymptotic root chains lie outside the root-locus region C σ0 for λ max < e hσ 0 |d| , while the region always has infinitely many roots for larger controller gains. Since it is not numerically possible to follow infinitely many roots, λ max is assumed to be smaller than e hσ 0 |d| .

Time delay as locus parameter
The root-locus equation w.r.t. the time delay is given in (3). The starting points of the root locus are the zeros of the transfer function 1 + G(s) inside the root-locus region C σ0 .
The branch points satisfy the root-locus equation (3) and Since the locus parameter λ = h in equation (18) is unknown, the direct computation of branch points is difficult compared to the previous case. Therefore we detect the branch points of the root locus while following the root-locus trajectory as described in Section 4.5.
The relative stability analysis of the closed-loop of SISO dead-time systems w.r.t. the stability boundary ℜ(s) = σ 0 is given in [8]. This analysis is based on the computation of boundary crossing roots of the root-locus equation in (3) which is similar to the computation for the controller gain case in Section 3.1. Further information can be found in [8]. Note that when the plant is biproper, λ max is assumed to be smaller than max 0, ln |d| |σ0| to have finitely many characteristic roots inside the rootlocus region C σ0 (see Remarks 2 and 3 in [8]).

Computing a Root-Locus Trajectory
We compute the root-locus trajectory of a characteristic root between two critical points. The trajectory satisfies the root-locus equation w.r.t the controller gain in (2) or the time delay in (3), represented by the locus parameter λ.
The starting points and the characteristic roots entering into the root-locus region are computed in Section 3. We follow each root-locus trajectory by a secant-predictor, Newton-corrector continuation method [1]. In the prediction step, a line passing through the last two computed roots and locus parameters is used to estimate the next root and the next locus parameter at a certain distance (step length) in the (characteristic root, locus parameter)-space. This estimate is corrected using Newton's method in the correction step. The next iteration continues in a similar way, though the step length is adaptive. The branch points for the controller gain are computed apriori and those for the time delay are detected while following the trajectory. The new branch direction is determined by Lemma 1 and the trajectory is followed until the upper bound of the locus parameter λ max is reached.

Parameterization of the root-locus trajectory
Continuation based methods follow the characteristic root s and the locus parameter λ based on the parameterization of the root-locus trajectory. In the literature, the characteristic root is parameterized w.r.t. the locus parameter as s = s(λ) for computational purposes [7,17]. In this parameterization, the characteristic roots are highly sensitive w.r.t. changes of the locus parameter in the neighborhood of branch points and relevant numerical problems are reported in [7]. We illustrate this sensitivity on the root locus by means of equation (2), where h = 4π 3 and G(s) = s 2 (s 2 +1)(s 2 +1) . Figure 1 shows the high sensitivity of the characteristic roots w.r.t. changes of λ at λ = 0 where a branch point occurs. As can be seen in Figure 1, this branch point appears as a turning point in the root versus locusparameter space. When a continuation method tries to follow the trajectory 1 or 2, the turning point around λ = 0 causes several numerical problems if λ is used as a continuation parameter. Take for instance trajectory 1. First the sensitivity or derivative of the root w.r.t. λ tends to infinity when λ → 0−. Second when the λ passes zero the same branch cannot be followed anymore because the original trajectory turns back. In the best scenario a point on trajectory 2 will be computed but this is unlikely with a local method because trajectory 2 has another direction than trajectory 1 (see Lemma 1), hence the estimate of the root obtained from the previous points will be bad. On the other hand, trajectory 1 and trajectory 2 appear as smooth curves in the complex plane as shown in Figure 1. In the field of numerical bifurcation analysis (see, e.g., [16]), this was the motivation to parameterize the trajectories in terms of a natural notion of an arclength, i.e., a distance measured along the trajectories.
In our method we parameterize the root-locus trajectories in terms of the arclength γ measured along the trajectory in the combined (characteristic root s, locus parameter λ) space, i.e., s = s(γ) and λ = λ(γ) for γ > 0.   Figure 1 on the complex plane.

Fig. 2. The root-locus trajectories 1 and 2 in
For the example given in (7) the two root-locus trajectories 1 and 2 are shown in Figure 3. The trajectories are smooth and using equal arclength steps for γ (where each step is shown as a point on the trajectory) the trajectories are traced without any numerical problem. Figure 3 also illustrates that a branch point is characterized by an extremum of the function γ → λ(γ) which can be used to detect the presence of a branch point along the trajectory.

Prediction step
The predicted root and the locus parameter computation in the prediction step requires the previous root, the locus parameter, a direction and a step length. For each root-locus trajectory, the starting point s 0 is available. The direction of the prediction step d i is computed as follows: • The initial directions d * 0 ∈ C for the poles of G inside C σ0 and the boundary crossing roots are computed by the root-locus equation, i.e., by solving ∂f ∂s (19) is the mathematical characterization of the tangent vector to the trajectory in the (s, λ) space [16].
• The directions in other iterations are computed using the real and imaginary parts of last two corrected roots and locus parameters, The real and imaginary parts of the predicted root and the locus parameter y p i+1 = [σ p i+1 ω p i+1 λ p i+1 )] T are computed using a line equation with a step length h i The geometric illustration of the prediction step is given in Figure 4 where the predicted point y p i+1 is shown as a gray dot. The initial step length h 0 is fixed. The step lengths in other iterations are calculated adaptively based on previous values as we outline later on.   Figure 1 in the (characteristic root, locus parameter)-space.

Correction step
We use Newton's method to solve a set of nonlinear equations to find the real and imaginary parts of the corrected root and the locus parameter These equations are given by The functions M and P are equivalent representations of the magnitude and phase equations of the root-locus equation. Equation (24) guarantees that the linearized distance of the corrected root and the locus parameter y c i+1 from the predicted root and the locus parameter y p i+1 is equal to the step size h i . Because of the linearization the proposed method can be seen as a pseudo arclenghth continuation method [16].
The function M and P are defined using the root-locus equation w.r.t. the controller gain in (2) or the time delay in (3) as Note that the controller gain k and the time delay h have constant values in the root-locus equation w.r.t. the time delay in (3) and the controller gain in (2) respectively.
We compute the corrected point y c i+1 by Newton's method on (22-24), i.e., by iteratively solving where J k i+1 and f k i+1 are the Jacobian and the gradient vector evaluated atỹ k i+1 . The iteration is initialized by the predicted point such thatỹ 0 i+1 = y p i+1 . We stop the iteration when the prescribed accuracy is reached and set the corrected point y c i+1 to the last point of the iteration. The geometric illustration is given in Figure 4.
Regarding the computational complexity, the above computation requires evaluating the functions on the left-hand side of the equations in (22-24) and the corresponding Jacobian and gradient of these equations at each iteration. This computation is numerically cheap and uses only the poles and zeros of G and the time delay h.

Adaptive step length
The step length computation for the next prediction step depends on two factors [1] • the contraction rate of the first two consecutive Newton steps in the corrector step, i.e., κ i+1 := • the distance to the root-locus trajectory The individual deceleration factors are calculated as κ df = κi+1 κ and δ df = δi+1 δ whereκ andδ are the nominal contraction rate and the distance. The overall deceleration factor h df of the step length is the maximum of individual deceleration factors, κ df and δ df limited to [ 1 2 , 2], i.e., h df := max{min{max{κ df , δ df }, 2}, 1 2 }.
Note that if h df = 2, the predictor step is repeated with a reduced step length. This check is done inside the corrector step to avoid unnecessary Newton iterations (see [1] for further details). The step length for the next prediction step is The adaptive step length selection allows us to follow the root-locus trajectory in a computationally efficient way with a prescribed accuracy. When the trajectory has sharp (resp. wide) curves, the step lengths are smaller (resp. larger). This behavior is illustrated in Figure 4.

Detection of branch points
We directly computed the branch points of the root locus w.r.t. the controller gain in Section 3.1. The branch point detection of the root locus w.r.t. the time delay requires monitoring the time delay parameter along the root-locus trajectory. By Lemma 1, the root-locus trajectory passes through the branch point when the multiplicity of the branch point is odd. In this case, there is no need to detect the branch point. Note that the locus parameter is strictly increasing for this case.
When the multiplicity of the branch point is even, the continuation method follows another root-locus trajectory after passing the branch point. The locus parameter on this trajectory is strictly decreasing since we follow the branch in the opposite direction. When we detect that the current time delay value is smaller than previous time delay values, we know that we have passed a branch point. Since we bracketed the branch point by the current and previous points, we then solve a set of nonlinear equations characterizing the branch point. Based on the multiplicity of the branch point, we decide the branch direction by Lemma 1. The branch point and its direction are added into the list for the next trajectories.

Algorithm
The overall algorithm is as follows.
(1) Compute the critical points of the root locus. Note that due to the symmetry of the spectrum, on the real axis it is sufficient to compute the branch points and continue from the next branch point, i.e., a continuation method is not needed. In our implementation, we set the nominal contraction rate, the nominal distance and the tolerance for the corrector step toκ = 0.5,δ = 10 −3 and 10 −5 .

Examples
We consider three examples with unique characteristics on their root locus. Example 7 has circular trajectories close to each other. Trajectories of Example 8 require different order of accuracy. Example 9 has a branch point, a trajectory going out of the region and trajectories with large complex numbers. We show root-locus plots only on the complex upper half-plane since the root-locus plots are symmetric w.r.t. the real axis due to real-valued coefficients of the plants in Examples.

Example 7
In [12], the stability of equation (3) is considered for the oscillator system, G(s) = ǫs 2 (s 2 +ω 2 1 )(s 2 +ω 2 2 ) . The behavior of the rightmost characteristic roots is analyzed as a function of the delay parameter for small values of a gain parameter. In a particular example, the parameters of the system are set to ǫ = 1, ω 1 = 2, ω 2 = 4 and the oscillatory nature of the rightmost characteristic root and the stability of the closed-loop system are shown.
We plot the root locus in equation (3) of this plant inside ℜ(s) ≥ −1 for λ ∈ [0, 5] in Figure 5. It is well-known that the characteristic roots of time delay systems can only cross the imaginary axis at a finite number of points when the delay is varied. This property can be clearly seen in Figures 5-8. Note in particular in Figure 8 how characteristic roots are passing through same points on the imaginary axis. Root-locus plots in Figures 6-8 illustrate the oscillatory behavior of the characteristic root trajectories centered around ±2j, ±4j and the asymptotic characteristic roots coming from the left. This phenomenon is qualitatively explained in [12]. The closed-loop system is stable for the time delay λ ∈ [0.83, 1.50] ∪ [4.11, 4.50]. The root-locus plot for equation (2) of this plant is given in Figure 9. As shown with cross markers inside gray regions in the figure, the plant has 6 unstable poles. They have all moved to the complex left half-plane at λ = 1.860 and the closed-loop system becomes stable for the con-troller gain λ ∈ [1.860, 4.469]. Note that the root-locus trajectories have lengths in the order of magnitudes, 1, 10 −2 , 10 −3 , 10 −4 respectively as illustrated in Figures 9-12. We follow each trajectory without any problem due to the adaptive step length in the algorithm. Example 9 Finally, we consider the SISO dead-time system where h = 1 and G(s) = s 2 −10s+50 s 3 +4s 2 +4.25s+1.25 . We plot its root-locus trajectories of equation (2) for the controller gain interval λ ∈ [0, 5] inside ℜ(s) ≥ −3.5. Figure 13 illustrates the general behavior of the root-locus trajectories. As the controller gain increases, more characteristic roots cross the boundary ℜ(s) = −3.5 and asymptotic root trajectories get closer to the imaginary axis. Figure 14-15 show the local behavior around the starting points, s = −0.5, −1, −2.5. In Figure 14, we see that the characteristic roots cross the imaginary axis at λ = 0.07 and the closed-loop system becomes unstable. In Figure 15, the root locus has a branch point s = −0.6976 at λ = 0.0009 and the trajectory starting with s = −2.5 leaves the region at λ = 0.0023.  Fig. 11. The root locus inside the bottom gray region in Figure 9 for controller gain λ ∈ [0, 6]. The computation times of root-locus plots for Examples 7-9 are 7.6, 1.4 and 2.5 seconds on a PC with an Intel Core Duo 2.53 GHz processor with 2 GB RAM respectively.

Concluding Remarks
The root locus of SISO dead-time systems with respect to the controller gain or the time delay is computed inside a given complex right half-plane. The method calculates the starting points of root-locus trajectories including the ones crossing the boundary of the root-locus region. Each trajectory is followed by a predictor-corrector continuation method. The characteristic roots on each root-locus trajectory are predicted by a secant method where the step length is adaptive and the predicted values are corrected by Newton's method. This continuation approach, which stems from numerical bifurcation analysis, is general and not restricted to parameterized equations of the form in equations (2,3). The implemen-   tation of the method is numerically stable for high order systems.