Reliable computation of PID gain space for general second-order time-delay systems

ABSTRACT This paper addresses the problem of determining the stability gain space of a PID controller for general second-order time-delay systems. First, a review of existing results and the associated drawbacks is presented. Subsequently, a new algorithm to compute the entire PID stability gain space is developed. The new algorithm is based upon existing results on the relationship between the stability of a quasi-polynomial and its derivatives, an extended version of the Hermit–Biehler theorem, and also the Nyquist criterion. The algorithm entails extraction of an admissible range for the PID parameter Kp, and then based on this range, a stability region in the (Ki − Kd) plane is computed. Well-known examples are studied to demonstrate the reliability and accuracy of the results.


Introduction
Despite a great deal of development in advanced controller design methods during the past decades, proportional-integral-derivative (PID) controllers are still widely used in industrial process systems. By various accounts, over 95% of the currently working controllers are of the PID type (Brisk, 2004). This widespread application is based on several advantages such as simplicity in control structure, control stability, effective performance, model-free tuning, robust behaviour, and easy maintenance (Astrom & Hagglund, 2006;Panda, 2008). Even where more advanced controllers are employed, the PID controller yet remains in charge of regulating the first control automation level (Visioli, 2006).
Modern digital industrial control systems are required to operate much more efficiently than their previous and older analogue counter parts. As a result, the current trend in tuning PID controllers is to move away from traditional design heuristics such as the Ziegler-Nichols method (Ziegler & Nichols, 1942), and towards powerful optimisation-based synthesis methods (Ho, 2003). In case of the latter, knowing the closed-loop stability region with respect to the (K i , K d , K p ) parameters is critical and necessary to be used as constraints of the tuning optimisation problem. In recent years, a number of methods to determine this space have been proposed. For example, Ho, Datta, and Bhattacharyya (1997) developed a linear programming method to determine the range of stabilising PI and PID controller parameters for a given system CONTACT M. Firouzbahrami fbahrami@ee.sharif.ir model. Based on the Nyquist criterion, Söylemez, Munro, and Baki (2003) developed a fast method for calculating stabilising PID controller parameters, and Tan, Kaya, and Atherton (2003) presented a new method for estimation of stabilising PI and PID controllers by plotting the stability boundary locus. The algorithms cited previously are only applicable to finite-dimensional systems. This is unfortunate since time-delay is an omnipresent feature of all industrial processes. More recently, attention has been shifted to algorithms which deal directly with non-finite dimensional time-delay systems. Initial works, such as Bhattacharyya (2001, 2002) and Martelli (2005), considered first-order time-delay systems. Second-order integrating systems with time-delay were studied in Ou, Tang, Gu, and Zhang (2005) and Ou, Zhang, and Gu (2006). Integral time-delay systems have also been studied in  and Padula and Visioli (2012). In addition to special system models, some solutions for general time-delay systems have also been proposed which are complex to evaluate (for example, see Almodaresi, Bozorg, & Taghirad, 2016;Bozorg & Termeh, 2011;Hohenbichler, 2009;Ou, Zhang, & Yu, 2009) or have a large computational cost (for example, see Yu, Le, Xian, & Wang, 2013). There are also seminal works in the domain of fractional order PID controllers some of which are Caponetton and Dongola (2013), Badri and Tavazoei (2013), and Hamamci (2008).
Second-order process models are of special interest since they are the lowest order model which can exhibit principal characteristics of real-life processes. In Wang (2007), a new graphical approach is proposed for the stabilisation of second-order time-delay systems with two real poles. However, according to Martelli (2009), 'the stability conditions were not explicit and no finite number of required computational steps is specified' . Moreover, as demonstrated in this paper, in some instances, the algorithm is not reliable. A second-order time-delay system with an additional finite zero has been considered in Martelli (2009). In this work, a new algorithm is developed to compute the entire stability gain space. While the algorithm is precise, it is complex and difficult to implement. Second-order time-delay systems with a pair of complex poles have been studied in Wang and Zhang (2010) by using a similar approach as presented in Wang (2007). Computation of stabilising PI and PID parameters for second-order time-delay systems with real unstable poles has also been considered in Farkh, Laabidi, and Ksouri (2014).
The contributions of this paper are two-fold. First, we demonstrate some inaccuracies in some of the existing results cited above and highlight what appears to be the source of these inaccuracies. Second, we propose an algorithm for computation of regions for PID stabilising gain space of general (real or complex poles) second-order time-delay systems. The algorithm computes efficiently and is based on an extended version of the Hermite-Biehler theorem and the Nyquist criterion. We also demonstrate that it is consistent with fundamental theorems on stability of quasi-polynomials and therefore returns a reliable solution.

Review on some previous results
Consider the standard closed-loop system as shown in Figure 1. Let the plant model G(s) be a second-order time-delay system defined as, where k > 0 is the steady-state gain, θ > 0 is the timedelay, and T 1 and T 2 are the time constants of the plant with arbitrary sign. System (1), as considered in Wang (2007), is a special class of the general second-order timedelay systems and contains two real poles only. An ideal PID-type controller is considered in the following form: (2) From system model (1) and controller (2), the characteristic equation for the closed-loop system in Figure 1 is Equations such as (3) are known as quasi-polynomials.
One of the earliest studies on the stability of quasipolynomials was carried out by Pontryagain (Pontryagain, 1955). Based on Pontryagain's results, a significant extension of the Hermite-Biehler theorem has been developed to study the stability of quasi-polynomials (see Bellman & Cook, 1963). The algorithm outlined in Wang (2007) first determines an admissible range for K p , and then for a specific value of K p within this range, an associated stability region for (K i − K d ) is computed by graphical means. Central to the methodology is the following result which sets out a necessary and sufficient condition for the existence of stabilising PID controllers.
Lemma 2.1 (Wang, 2007): A necessary and sufficient condition for existence of a stabilising PID controller with stability degree σ > 0 for the plant described in (1) is The value of σ determines the line s = −σ which will constrain all the closed-loop poles to its left.
From Lemma 2.1 for a given plant and stability degree σ > 0, Inequality (4) becomes a function of the controller parameter K p . It therefore appears that there always exists some K p which will satisfy (4). Accordingly, there must exist at least one PID controller for any second-order system of the form specified in (1) which will place the closed-loop poles of the system to the left of an arbitrary line on the negative real axis. Note that from Lemma 2.1, a necessary and sufficient condition for existence of stabilising PID controller is 6(n + kK p e θσ ) + σ 2 m m = 6(n + kK p e θσ ) m + σ 2 > σ 2 , or equivalently, n + kK p e θσ m > 0.
It is clear that for any values of m and n, a K p can be found (positive or negative) such that (7) is satisfied. This assertion is correct for delay-free second-order systems, but not generally correct, as shown in Example 2.1.
Example 2.1: Consider the following second-order timedelay system: Take σ = 0.1. From Lemma 2.1, 6.9458 6(1.3158 which indicates that for K p > −0.9747, there exists at least one stabilising PID controller. In addition, Remark 2 in Wang (2007) emphasises that the stability region is nonempty since for K p > −0.9747 there exists some ε > 0 such that, < 6.9458 6(1.3158 + K p e 0.3 ) + 0.069458 To see why this cannot be true, consider that if there exists a PID controller which places the closed-loop poles to the left of s = −0.1, then there must exist at least one PID controller which places the poles to the left of the origin. However, according to Theorem 1 in Wang (2007), there does not exist any stabilising PID controller for this plant. This theorem states that for a stabilising PID controller to exist, (11) must be satisfied in the interval (0, π), However, from Figure 2 , it is clear that for the plant specified by (8), Equation (11) does not have any roots in the said interval.
There is also a secondary contradiction which appears to stem from an earlier proposition in Silva et al. (2002)  on the relationship between the stability of a quasipolynomial and its delay-free version. It states 'a minimal requirement for any control design is that the delayfree closed-loop system be stable' (also see, Remark 2 and Theorem 1 in Wang, 2009, and Condition 1 in Wang, 2012, andalso Farkh et al., 2014;Marquez-Rubio, del Muro-Cuéllar, and Ramírez, 2014;Wang & Zhang, 2010). This statement is also quoted elsewhere in the literature (for example, see Silva et al., 2001;Silva, Datta, and Bhattacharyya, 2004) and forms the basis for several recent works in this topic (Farkh et al., 2014;Keel & Bhattacharyya, 2008;Marquez-Rubio et al., 2014;Silva et al., 2004;Wang, 2007Wang, , 2009Wang, , 2012Xu, Datta, & Bhattacharyya, 2003).
Despite extensive referrals, a formal proof of the proposition appears to be missing. In addition, there are counterexamples which show that the statement might not be correct. For example, Kharitonov, Niculescu, Moreno, and Michiels (2005) show that for the proportional-only control, the stability gain space of delayed and delay-free versions of a plant can have an empty intersection.
Example 2.2 (Wang, 2007): Consider the following second-order time-delay system: From, Equation (15a) in Wang (2007), the delay-free version of the plant is stable if and only if, This result can be simply obtained from Routh Hurwitz table. Now, consider the following two PID controllers: Since C 1 and C 2 invalidate the the second and the third conditions of (13), respectively, they will not stabilise the delay-free plant, and by implication of the statement should not also stabilise the delayed version either. However, both of these controllers will actually stabilise the delayed version of the system. The step responses of the closed-loop system with C 1 and C 2 are shown in Figure 3.

Main results
In this section, we consider a general form of secondorder time-delay systems with real or complex poles. For this general form, we first develop new necessary conditions to compute the stability gain space of the PID controller. We then make direct use of the Hermit-Biheler theorem and the Nyquist criterion to develop a simple and reliable algorithm for computation of the PID stabilising gain space.

Problem formulation and preliminary results
Consider the closed-loop system as shown in Figure 1 with the general second-order time-delay system G(s), Without loss of generality, we assume k > 0 in the remainder of this paper. The controller C(s) we consider is in the form of (2). From (15) and (2), the characteristic equation of the system in Figure 1 becomes Throughout the paper, we consider only the PID parameter K i to be positive. Alternatively, the results can be simply extended to the region of K i < 0 using methods presented in this work. Now, multiplying both sides of (16) by e θs gives Since e θs has no finite zeros, the quasi-polynomial (16) is Hurwitz stable if and only if (17) is Hurwitz stable. That is, the unity feedback loop in Figure 1 is stable if and only if all zeros of (17) are in the left half complex plane. The following theorem provides a necessary and sufficient condition for stability of quasi-polynomials such as f(s) in (17).
Theorem 3.1 (Bellman & Cook, 1963): Consider the following quasi-polynomial: where θ 1 < θ 1 < < θ r , θ r + θ 1 > 0 and δ 0r ࣔ 0. The quasi-polynomial (18) is stable if and only if, (1) f r (w) and f i (w) have only real roots and these roots interlace, where f r (w), f i (w), f r (w), and f i (w) denote the real and imaginary parts of f(jw) and f (jw) (the first derivative of f(jw) with respect to w), respectively.
Theorem 3.1 is known as the extended Hermite-Biehler theorem. Other important theorems which have been used in the development of our proposed conditions are the following two theorems.
Theorem 3.2 (Bellman & Cook, 1963): Consider the quasi-polynomial (18) in Theorem 3.1. Let M and N denote the highest powers of s and e θs in f(s). Also, suppose that f r (w) and f i (w) are the real and imaginary parts of f(jw). Let η be a constant such that the coefficients of the highest degree in f r (w), and f i (w) do not vanish at w = η. Then, a necessary and sufficient condition for f r (w) or f i (w) to have only real roots is that in the interval −2lπ + η ࣘ w ࣘ 2lπ + η, f r (w) or f i (w) have exactly 4lN + M real roots for a sufficiently large l. Theorem 3.3 (Kharitonov et al., 2005)

: Consider the quasi-polynomial f(s) in (18). If f(s) is stable, f (s) (derivation of f(s) with respect to s) is also a stable quasipolynomial.
A straightforward conclusion of Theorem 3.3 is that if f (s) is an unstable quasi-polynomial, then f(s) is also unstable.

Necessary conditions for PID stabilisation
In this section, we use Theorem 3.3 to derive tight parameter bounds which are then used as a set of necessary conditions for existence of a PID stabilising controller for the plant considered in (15). (17).

Lemma 3.1: Consider the quasi-polynomial f(s) in
Proof: Computing the third derivation of f(s) in (17) with respect to s yields where f a (s) = (θ 3 s 3 + θ 2 (αθ + 9)s 2 + θ (βθ 2 + 6αθ + 18)s + (βθ 2 + 2αθ + 2)). The conditions proposed in Lemma 3.1 are consistent with existing result by Lee, Wang, and Xiang (2010). Using their results, it can be easily shown that the third statement of (19) is a necessary and sufficient condition for existence of PID stabilising controllers for the special class of plant (15) in which there are only one stable and one unstable pole. The plant considered in Lee et al. (2010) can contain more than two poles; however, unlike the system considered in this work, it does not permit two unstable (or two stable) poles at the same time and also oscillatory transients. Moreover, the entire stability area was not considered in Lee et al. (2010). Lemma 3.2 is another result which can be obtained using Theorem 3.3. It is used as a necessary condition to compute the stabilising gain space of the PID parameters, in the final result of the paper. (17). Let

Lemma 3.2: Consider the quasi-polynomial f(s) in
and n be the number of unstable poles of f c (s) in (22). Then, f(s) is stable only if the Nyquist plot of H (s) = 2ke −θ s f c (s) encircles the ( − 1, 0) point −n times. (17),

Proof: From
where f c (s) is defined in (22). Evidently, f c (s) in (23) where GM 1 and GM 2 are the gain margin of −H(s) and H(s), respectively.

An admissible range for K p
Substituting s = jw into (17) yields Taking z = θw, the real and imaginary parts of the characteristic equation (25) can be separated as follows: In (27), only the proportional term of the PID controller appears. Accordingly, Theorem 3.2 can be used to obtain a range for K p under which f i (z) has only real roots. By rewriting f i (z), we have wherē Consider the following lemmas which are used in deriving an admissible range for the PID parameter K p in the next section. (27) and (29). Then, f i (z) has only simple real roots if and only iff i (z) + kK p has 2l + 1 zeros in the interval 0 < z ࣘ 2lπ for a sufficiently large value of l.

Lemma 3.3: Consider f i (z) andf i (z) as defined by
Proof: From Theorem 3.2, we have M = 3, and N = 1. Hence, f i (z) only has simple real roots if and only if (1) the number of its roots in the interval −2lπ + η ࣘ z ࣘ 2lπ + η for sufficiently large values of l is equal to 4l + 3, and (2) the term cos (z) which is the coefficient of the highest power in f i (z) does not vanish for z = η.
To satisfy the second condition, let η = 0. Then, f i (z) must have 4l + 3 zeros in −2lπ ࣘ z ࣘ 2lπ for sufficiently large values of l to satisfy the first condition. Considering the relationship between f i (z) andf i (z) + kK p in (28), and thatf i (z) + kK p is an even function of z, implies f i (z) has 4l + 3 zeros in −2lπ ࣘ z ࣘ 2lπ if and only iff i (z) + kK p has 2l + 1 zeros in the interval 0 < z ࣘ 2lπ.

Remark 3.2:
Since the number of roots is a natural number, the value of l in Lemma 3.3 must be from the set of half natural numbers L = 0.5, 1, 1.5, 2, 2.5, ….
It is known that for every A, B, and z, From (29) and (31), we havē where A = −α z θ , and B = β − z θ 2 . From Lemma 3.4, for β ࣙ α 2 /2, the amplitude off i (z) in (32) has a minimum point which is used in computing the admissible range of K p . Specifically, for the roots interlacing property that appears in the next section, all local maximum and minimum points off i (z) + kK p must be positive and negative, respectively. We summarise this condition in Lemma 3.5. (29). Define E 1 and E 2 as follows:

Therefore, the extremum points off i (z) + kK p do not have any change of sign if and only if
Proof: The proof is trivial by considering the definitions of E 1 and E 2 .

Remark 3.3: One of the local maximum points off i (z)
is at z = 0. As a result, from Lemma 3.5,f i (z) + kK p for z = 0 in the K p admissible range must have positive sign. Then, from (27), a necessary condition for Lemma 3.5 to be satisfied isf i (0) + kK p > 0 or, equivalently, K p > − β k . It is now possible to propose a new algorithm which can be effectively used to derive an admissible range for the PID parameter K p .
Step 2: From Lemma 3.5, determine the values of E 1 and E 2 and compute the interval (33).
Step 3: Determine some values of l ࢠ L (L is defined in Remark 3.2) and for each one a relevant interval (K p1 , K p2 ) such that for all values of K p in the interval, the condition of Lemma 3.3 is satisfied. Step 4: The admissible range of K p is the intersection region of (E 1 , E 2 ) and the union of the intervals (K p1 , K p2 ).

Remark 3.4:
For the special case β ࣘ α 2 /2, it is possible to propose some direct equations to compute the admissible range of K p (for example, see Wang, 2007). However, in the general case, using such results can lead to an inaccurate range, especially when β α 2 /2.
(34) Figure 4 shows the functionf i (z) along with two limit values of E 1 = 3 and E 2 = −2.984 for G 1 (s). Therefore, for Step 2 in Algorithm 3.1, we have Now, let l = 1, K p1 = −1.5 and K p2 = 1.474 in the third step of the algorithm. These values are simply realised by considering the two points p 1 and p 2 in Figure 4 as the maximum values for whichf i (z) − p 1 andf i (z) − p 2 retain its three roots in the interval 0 < z ࣘ 2π. Therefore, Step 3 also leads to (35). Consequently, the combination of the two achieved intervals from Steps 2 and 3 leads to an admissible range for the PID parameter K p which is the same as (35). It was also possible to let l = 0.5 in Step 3. However, in this case, (K p1 , K p2 ) = ( − 0.695, 1.474), which is clearly a reduced subset of the obtained result. For G 2 (s), the functionf i (z) and its constraints E 1 = 9.279 and E 2 = −9.5 are shown in Figure 5.
Step 2 in Algorithm 3.1 gives For the third step, first let l = 2, K p1 = −0.4639, and K p2 = −0.1245 which are found by considering the points p 11 and p 21 in Figure 5 as the maximum possible shift off i (z) (downward and upward) for which it keeps its five roots in the interval 0 < z ࣘ 4π. Alternatively, let l = 2.5, K p1 = −0.3854, and K p2 = 0.475 which are determined by p 12 and p 22 in Figure 5. By combining the results, we have K p ࢠ ( − 0.4639, 0.475) for Step 3. Finally, from Steps 2 and 3, the admissible range for the PID parameter K p is found to be as (36).

Stability region in the (K i − K d ) plane
In this section, for any given value of K p within its admissible range, we derive a region for (K i − K d ) in which the interlacing property for the roots of f i (z) and f r (z) in (26) and (27)   where z h for h = 1, 2, … are the positive roots of f i (z), and Proof: From (26) and (27), where From (31), (39) can be rewritten as Suppose the parameters K i , K p , and K d are equal to zero.
The two functions f i (z) and f r (z) will then have a relative phase shift equal to π/2 but with equal amplitude and period. Therefore, the roots of f i (z) and f r (z) interlace. If any of the PID parameters are nonzero, a dynamic (non-periodic) y-intercept is added to each of the functions f i (z) and f r (z) and one of the following cases must occur: Case1: Some roots of f i (z) do not interlace any roots of f r (z) (see Figure 6 as an example). Case2: The roots of f r (z) and f i (z) interlace (see Figure 7 as an example). Case3: Some roots of f r (z) and f i (z) interlace pairwise (see Figure 8 as an example).
Since K i > 0, then f r (z 0 ) = kK i > 0 (since z 0 = 0). Therefore, from (39), f r (z) can be restated as follows: where a(z) and b(z) are defined in (38). Since kz 2 θ 2 has a positive value, condition (37) is equivalent to Satisfying (43) implies that between every two consecutive roots of f i (z), there exists at least one root of f r (z) and therefore only Case 2 can occur. This completes the proof.
Summary: For any given value of K p within its admissible range resulted by Algorithm 3.1, the intersection region of the infinite number of inequalities in (37) determines a set of (K i , K d ) values under which the interlacing property of Theorem 3.1 for the roots of f r (z) in (26) and the roots of (27)  Step 3: Using Algorithm 3.1, find an admissible range for the proportional gain K p . In addition, compute the boundary values K − d and K + d as the lower and the upper boundaries of PID parameter K d using Lemma 3.2.
Step 4: Select a fixed value for K p within its admissible range.
Step 5: Determine e and o as defined in Theorem 3.4.
Step 6: Compute the intersection region of the three boundaries K − id from below, K + id from above, and K i = 0 from the left-hand side in the (K i − K d ) plane.
Step 7: Go to Step 4. Example 3.2: Consider again the plant G 1 (s) in (34) from Example 3.1. This system satisfies the necessary conditions outlined in Lemma 3.1. From Example 3.1, the admissible range for the PID parameter K p is (− 1.5, 1.4740). In addition, from Lemma 3.2, the admissible range for K d is ( − 3.25, 5.7116). Taking K p = 1 for illustration purposes, we find o = 3 and e = 8. The two lines (44) for h = 1, 3 and K d = K − d determine K − id , and the lines (44) for h = 2, 4, 6 and K d = K + d determine K + id . All of the stated lines and the resulting intersection region with the vertical line K i = 0 are plotted in Figure 9. The steps are then repeated for different values of K p . The resulting intersection regions are shown in Figure 10. This system satisfies the necessary conditions obtained in Lemma 3.1. Moreover, from Algorithm 3.1, the range for K p is determined as follows: which, in this instance, is similar to the interval obtained for K p in Wang (2007). The transfer function H(s) defined in Lemma 3.2 becomes H(s) = 4 e −0.2s 0.04s 3 + 1.08s 2 + 3.68s − 5.2 Computing the poles of H(s) shows it has an unstable pole at s = 1.0662. The Nyquist plot of H(s) is shown in Figure 11. The interval between P 1 and P 2 refers to where the Nyquist diagram encircles the point (−1, 0) once anticlockwise. Accordingly, from the Nyquist theorem and  Lemma 3.2, the admissible range for the PID parameter K d can be determined as follows: where P 1 = −0.2262 and P 2 = −0.7692. For illustration purposes, we select two different values −0.95 and −0.98 for K p . This leads to the two stability regions in the (K i − K d ) plane shown in Figure 12. Both controllers C 1 (s) and C 2 (s) which were given in (14) are within this region. Accordingly, the proposed algorithm is reliable and consistent with the results obtained in Example 2.2.

Conclusion
This paper studied the problem of determining the stability gain space of PID controllers for general secondorder time-delay systems. Based on an extension of the Hermite-Biehler theorem, and the Nyquist criterion, a new simple algorithm for determining the stability gain space has been proposed. The algorithm not only applies to both stable and unstable systems, but can also be utilised for systems with real or complex poles. To investigate the validity of the proposed algorithm, three different types of second-order time-delay systems have been studied.

Disclosure statement
No potential conflict of interest was reported by the authors.