# Computation of All Stabilizing PID Gain for Second-Order Delay System

Authors: Rihem Farkh; Kaouther Laabidi; Mekki Ksouri
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

The problem of stabilizing a second-order delay system using classical proportional-integral-derivative (PID) controller is considered. An extension of the Hermite-Biehler theorem, which is applicable to quasipolynomials, is used to seek the set of complete stabilizing PID parameters. The range of admissible proportional gains is determined in closed form. For each proportional gain, the stabilizing set in the space of the integral and derivative gains is shown to be either a trapezoid or a triangle.


---

## Body

## 1. Introduction

Dead times are often encountered in various engineering systems and industry processes such as electrical and communication network, chemical process, turbojet engine, nuclear reactor, and hydraulic system. In fact, delays are caused by many phenomena like the time required to transport mass, energy or information, the time processing for sensors, the time needed for controllers to execute a complicated algorithm control, and the accumulation of time lags in a number of simple plants connected in series [1]. Delay makes system analysis and control design much more complex which is explained as follows. First, the effect of the control action takes some time to be felt in the controlled variable. Second, the control action that is applied based on the actual error tries to correct a situation that originated previously. Finally, in the frequency domain, the time delay introduces an extra decrease in the systems phase, which may cause instability [1, 2].Today, proportional-integral (PI) and proportional-integral-derivative (PID) controller types are the most widely used control strategy. It is estimated that over 90% of process control applications employ PID control thanks to its essential functionality and structural simplicity [3]. Since the minimal requirement for PID controller is to guarantee the system stability. It is desirable to know the complete set of stabilizing PID parameters before tuning and design. Surveys have reported that poor tuning, configuration errors, traditional and empirical techniques such as Ziegler Nichols or Cohen-Coon tuning methods are found in the industrial application of PID [4]. A great deal of academic and industrial effort has focused on improving PID control in the areas of tuning rules to decrease the rising gap between engineering practice and control theory. Recent studies have made use of the generalization of Hermite-Biehler theorem to compute the set of all stabilizing PID controllers for a given plant [5]. Since almost all plants encountered in process control contain time delays, so computing the complete set of PID controllers that stabilize time delay system is of considerable importance. A robust PI/PID design via numerical optimization for delay processes was proposed in [6]. Padma Sree et al. [7] propose a PI/PID controllers design for first-order delay system by extracting the coefficients of the numerator and denominator of the closed-loop transfer function. In the work presented in [8], the authors have developed a graphical approach which is based on D-partition method [9], the borders of absolute and relative stability regions are described in the parameter space. A unified approach and Delta operator are used to obtain unified stability boundaries for PI, PD, and PID controller for an arbitrary order delay system; where the boundaries can be found when only the frequency response and not the parameters of the plant are known [10]. Eriksson and Johansson [11] consider a design of PID controller for system with varying time delay using multi-objectives’ optimization. Neural network has been widely applied for nonlinear dynamical system identification via PID controller [12]. Zhang et al. [13] present an intelligent nonlinear PID controller based on the recurrent neural networks (RNNs) to control multivariable plants by using two predictive control schemes. In [14, 15], the authors propose a new PID neural networks (PIDNN) which is a dynamic multilayer network based on P, I, and D neurons. This PIDNN is used to control multivariable plants [15] and has also the ability to control time delay system [14]. In [16], Roy and Iqbal have explored PID tuning of first-order delay system using a first-order Padé approximation and the Hermite-Biehler stabilization framework. The set’s characterization of all stabilizing PI/PID parameters using a version of the Hermit-Biehler theorem for first-order delay system is presented in [5, 17–20]. However, these results are not applicable to the second-order delay system. Although, the PID stabilization for a given time delay plant has been developed [21] using the generalized Nyquist criterion, the stabilizing set of the proportional gain cannot be found. In [22–24], the characterization of the set of all stabilizing P/PI/PID parameters is given by using the Hermit-Biehler Theorem for polynomials for a class of time delay system which verify the interlacing property at high frequencies. This method, however, is complex and there are difficulties to achieve the stabilizing range of proportional parameter. As shown in [23, Example  2], this latter approach does not determine the precise range of proportional gain which is larger than the exact stabilizing range given by [18]. Furthermore, a known stabilizing value of (Kp,Ki,Kd) is first determined using Nyquist criterion, then the entire stability region is established. In our work, the stabilizing problem of PID controller for second order delay system is analyzed using the extended Hermit-Biehler theorem for quasi-polynomials. The exact range of stabilizing proportional gain is first determined. The set of stabilizing integral and derivative constant values are then derived. Our algorithm is more simple and faster in time computing than the other presented in [24].


## 2. Preliminary Results for Analyzing Time Delay System

Several problems in process control engineering are related to the presence of delays. These delays intervene in dynamic models whose characteristic equations are of the following form [5, 25](2.1)δ(s)=d(s)+e-L1sn1(s)+e-L2sn2(s)+⋯+e-Lmsnm(s),
where: d(s) and ni(s) are polynomials with real coefficients and Li represent time delays. These characteristic equations are recognized as quasi-polynomials. Under the following assumptions:(A1)
deg(d(s))=n,deg(ni(s))<n,fori=1,2,…,m,(A2)
L1<L2<⋯<Lm.One can consider the quasi-polynomialsδ*(s) described by(2.2)δ*(s)=esLmδ(s),δ*(s)=esLmd(s)+es(Lm-L1)n1(s)+es(Lm-L2)n2(s)+⋯+nm(s).
The zeros of δ(s) are identical to those of δ*(s)since esLm does not have any finite zeros in the complex plane. However, the quasi-polynomial δ*(s) has a principal term since the coefficient of the term containing the highest powers of s and es is nonzero. If δ*(s)does not have a principal term, then it has an infinity roots with positive real parts [5].The stability of the system with the characteristic equation (2.1) is equivalent to the condition that all the zeros of δ*(s)must be in the open left half of the complex plan. We said that δ*(s)is Hurwitz or is stable. The following theorem gives a necessary and sufficient condition for the stability of δ*(s).Theorem 2.1 (see [5]).
Letδ*(s)be given by (2.2), and write
(2.3)δ*(jω)=δr(ω)+jδi(ω),
where δr(ω) and δi(ω) represent, respectively, the real and imaginary parts of δ*(jω).
Under conditions (A1) and (A2), δ*(s) is stable if and only if(1)
δr(ω) and δi(ω) have only simple, real roots and these interlace,(2)
δ′i(ω0)δr(ω0)-δi(ω0)δ′r(ω0)>0for some ω0in [-∞,+∞],where δ′i(ω) and δ′r(ω) denote the first derivative with respect to ω of δr(ω) and δi(ω), respectively.A crucial stage in the application of the precedent theorem is to make sure thatδr(ω) and δi(ω) have only real roots. Such a property can be checked while using the following theorem.Theorem 2.2 (see [5]).
LetM and N designate the highest powers of s and es which appear in δ*(s). Let η be an appropriate constant such that the coefficient of terms of highest degree in δr(ω) and δi(ω) do not vanish at ω=η. Then a necessary and sufficient condition that δr(ω)and δi(ω) have only real roots is that in each of the intervals -2lπ+η<ω<2lπ+η,l=l0,l0+1,l0+2⋯
δr(ω) or δi(ω) have exactly 4lN+M real roots for a sufficiently large l0.


## 3. PID Control for Second Order Delay System

A second order system with delay can be mathematically expressed by a transfer function having the following form:(3.1)G(s)=Ks2+a1s+a0e-Ls,
where K is the static gain of the plant, L is the time delay, and a0 and a1 are the plant parameter’s which are always positive. The characteristic equation of the closed-loop system is given by(3.2)δ(s)=K(Ki+Kps+Kds2)e-Ls+(s2+a1s+a0)s,
we deduce the quasi-polynomial δ*(s):(3.3)δ*(s)=eLsδ(s)=K(Ki+Kps+Kds2)+s(s2+a1s+a0)eLs,
by replacing s by jω, we get(3.4)δ*(jω)=δr(ω)+jδi(ω),
with:(3.5)δr(ω)=KKi-KKdω+(ω3-a0ω)sin(Lω)-a1ω2cos(Lω),δi(ω)=w[KKp+(a0-ω2)cos(Lω)-a1ωsin(Lω)].Clearly, the parametersKi and Kd only affect the real part of δ*(jω) whereas the parameter Kp affects the imaginary part.Let’s putz=Lω, we get(3.6)δr(z)=KKi-KKdz2L2+sin(z)(z3L3-a0zL)-a1z2L2cos(z),δi(z)=zL(KKp+cos(z)(a0-z2L2)-a1zLsin(z)).The application of the second condition of Theorem2.1, led us to(3.7)E(z0)=δi′(z0)δr(z0)-δi(z0)δr′(z0)>0,
from (3.6) we have(3.8)δi′(z)=KKpL-sin(z)(a0+2a1zL2-z3L3)+cos(z)(a0L-3z2L3-a1z2L2),
for z0=0 (a value that annul δi(z)), and we get:(3.9)E(z0)=δi′(z0)δr(z0)=(KKp+a0L)KKi>0,
which implies Kp>-a0/Ksince K>0 and Ki>0.We pass to the verification of the interlacing condition ofδr(z) and δi(z) roots. For such purpose, we are going to determine the roots from the imaginary part, which is translated by the following relation:(3.10)δi(z)=0⇒z=0orKKp+cos(z)(a0-z2L2)-a1zLsin(z)=0,⇒z=0orKKp+cos(z)(a0-z2L2)=a1zLsin(z),⇒z=0orf(z)=g(z).
We notice the z0=0 is a trivial root of the imaginary part. The others are difficult to solve analytically. However, this can be made graphically. The following two cases are presented.Case 1 (Kp>-a0/K).
In this case, we graph the curves off(z) and of g(z) which are presented in Figure 1. 							(3.11)z1∈[0,π/2],z3∈[3π/2,2π],z5∈[7π/2,4π],⋮z2∈[π/2,π],z4∈[5π/2,3π],z6∈[9π/2,5π],⋮
that is, zj verifies
(3.12)z1∈[0,π2],zj∈[(2j-3)π2,(j-1)π],forj≥2.Figure 1
Representation of the curves off(z) and of g(z) (K=L=a1=1, a0=2, Kp=1.3, and Ku=1.5884) (Case: Kp<Ku).Case 2 (Kp≥Ku).
Figure2 represents the case where the curves of f(z)=KKu+cos(z)(a0-(z2/L2)) and of g(z) are tangent in [0,π], where Ku is the largest number of Kp.
The plot in Figure3 corresponds to the case where Kp>Ku and the plot of f(z) does not intersect g(z) twice in the interval [0,π].Figure 2
Representation of the curves off(z) and of g(z) (K=L=a1=1, a0=2, Kp=Ku=1.5884) (Case: Kp=Ku).Figure 3
Representation of the curves off(z) and of g(z) (K=L=a1=1, a0=2, Kp=2, and Ku=1.5884 (Case: Kp>Ku).Theorem2.2 is used to verify that δi(z) possess only simple roots. By replacing Ls by s1 in (2.2), we rewrite δ*(s) as follows:(3.13)δ*(s)=eLsδ(s)=es1δ(s1)=es1((s1L)3+a1(s1L)2+a0s1L)+K(Kps1L+Ki).
For this new quasi-polynomial, we see that M=3 and N=1 where M and N designate the highest powers of s and es which appear in δ*(s). We choose η=π/4that satisfies the condition giving by theorem 3 as δr(η)≠0 and δi(η)≠0. According to Figure 1, we notice that for -a0/K<Kp<Ku, δi(z) possess four roots in the interval [0,2π-π/4]=[0,7π/4] including the root at origin. As δi(z) is odd function of z, so it possesses seven roots in [-2π+π/4,2π-π/4]=[-7π/4,7π/4]. Hence, we can affirm that δi(z) has exactly 4N+M=7in [-2π+π/4,2π+π/4]=[-7π/4,9π/4]. In addition, it can be shown that δi(z)has two real roots in each of the intervals [2lπ+π/4,2(l+1)π+π/4] and [-2(l+1)π+π/4,-2lπ+π/4] for l=1,2,…. It fallows that δi(z) has exactly 4lN+Mreal roots in [-2lπ+π/4,2lπ+π/4] for -a0/K<Kp<Ku. At the end, according to Theorem 2.2, δi(z) has only real roots for every Kp in [-a0/K,Ku]. For Kp≥Ku, corresponding to Figures 2 and 3, the roots of δi(z)are not real. We pass to determine the superior value of Kp. According to the definition of Ku, if Kp=Ku, then the curves of f(z)and g(z) are tangent in the point α. Which is translated by(3.14)KKu+cos(α)(a0-α2L2)=a1αLsin(α),ddz[KKu+cos(z)(a0-z2L2)]z=α=ddz[a1zLsin(z)]z=α,⇒-2αcos(α)(1+a1L)+sin(α)(α2-a0L2-a1L)=0⇒tan(α)=α(2+a1L)(α2-a0L2-a1L),
once α is determined, the parameter Ku is expressed by the relation:(3.15)Ku=1K(a1αLsin(α)-cos(α)(a0-α2L2)).Theorem 3.1.
Theδi(z) have real roots if and only if:
(3.16)-a0K<Kp<1K(a1αLsin(α)-cos(α)(a0-α2L2)),
where α is the solution of the equation tan(α)=α(2+a1L)/(α2-a1L-a0L2) in the interval [0,π].After determination of the roots of the imaginary partδi(z), we pass to the evaluation of these roots by the real part δr(z)(3.17)δr(z)=KKi-KKdz2L2+sin(z)(z3L3-a0zL)-a1z2L2cos(z)=Kz2L2[-Kd+KiL2z2+LKz[-a1zLcos(z)+sin(z)(z2L2-a0)]]=Kz2L2[-Kd+m(z)Ki+b(z)],
where:(3.18)m(z)=L2z2,b(z)=LKz[-a1zLcos(z)+sin(z)(z2L2-a0)].
From condition 1 of Theorem 2.1, the roots of δr(z) and δi(z) have to interlace for δ*(s) to be stable. We evaluate δr(z) at the roots of the imaginary part δi(z).Forz0=0, we have(3.19)δr(z0)=KKi>0,
for zj≠z0, where j=1,2,3,…, we get:(3.20)δr(zj)=Kzj2L2[-Kd+m(zj)Ki+b(zj)]=Kzj2L2[-Kd+mjKi+bj],
where:(3.21)m(zj)=mj,b(zj)=bj.
Interlacing the roots of δr(z) and δi(z) is equivalent to δr(z0)>0 (since Ki>0), δr(z1)<0,δr(z2)>0… .We can use the interlacing property and do it asδi(z) which has only real roots to reach that δr(z) possess real roots too.From the previous equations we get the following inequalities:(3.22){δr(z0)>0δr(z1)<0δr(z2)>0δr(z3)<0δr(z4)>0⋮⇒{Ki>0Kd>m1Ki+b1Kd<m2Ki+b2Kd>m3Ki+b3Kd<m4Ki+b4⋮
Thus intersecting all these regions in the (Ki,Kd) space, we get the set of (Ki,Kd) values for which the roots of δr(z) and δi(z) interlace for a fixed value of Kp. We notice that all these regions are half planes with their boundaries being lines with positive slopes mj.Example 3.2.
Consider the plant given by relation (3.1) with the following parameters K=L=a1=1 and a0=2:
(3.23)G(s)=1s2+s+2e-s.
The imaginary part δi(z) has only simple real roots if and only if -1<Kp<1.58.
We set the controller parameterKp to 1.3, which is inside the previous range. For this Kp, δi(z) takes the form:
(3.24)z(1.3+cos(z)(2-z2)-zsin(z))=0,
We next compute some of the positive real roots of this equation and arrange them in increasing order of magnitude:
(3.25)z0=0,z1=1.3608,z2=1.8905,z3=4.9829,z4=7.9619,z5=11.0976,z6=14.2017.
Using (3.18) we calculate the parameters mj and bj for j=1,…,6:
(3.26){m1=0.54,m2=0.2798,m3=0.0403,m4=0.0158,{b1=0.54,b2=-0.3150,b3=1.1047,b4=-4.6822.
Interlacing the roots of the real and the imaginary part occurs for Kp=1.3, if and only if the following inequalities are satisfied:
(3.27)Ki>0,Kd>0.54Ki-0.3150,Kd<0.2798Ki+1.1047,Kd>0.0403Ki-4.6822,Kd<0.0158Ki+7.7736.
The boundaries of these regions are illustrated in Figure 4.
In this case, the stability region is defined by only two boundaries:(3.28)Kd=m1Ki+b1,Kd=m2Ki+b2,
because we have the following inequalities:
(3.29)bj<bj+2,forj=2,4,6,…,bj>bj+2,forj=1,3,5,7,…,mj<mj+2,forj≥1.Figure 4
Region boundaries of Example3.2.Example 3.3.
Consider the plant (3.1) with the following parameters K=2,a0=3,a1=1,L=2:(3.30)G(s)=2s2+s+3e-2s.
The imaginary part δi(z) has only simple real roots if and only if -1.5<Kp<0.98.
We now set the controller parameterKp to 0.5, which is inside the previous range. For this Kp, δi(z)takes the form
(3.31)z(0.5+cos(z)(2-z24)-0.5zsin(z))=0.
We next compute some of the positive real roots of this equation and arrange them in increasing order of magnitude
(3.32)z0=0,z1=1.6480,z2=2.9830,z3=5.4560,z4=8.077,z5=11.22,z6=14.26,z7=17.41.
Using (3.18), we calculate the parameters mj and bj for j=1,…,7:(3.33){m1=1.4728,m2=0.4495,m3=0.1344,m4=0.0613,m5=0.0318,m6=0.0197,m7=0.0132,{b1=-1.3656,b2=0.4527,b3=-0.9377,b4=1.7176,b5=-2.5853,b6=3.3906,b7=-4.2097.
Interlacing the roots of the real and imaginary part occurs for Kp=0.5, if and only if the following inequalities are satisfied:
(3.34)Ki>0,Kd>1.4728Ki-1.3656,Kd<0.4495Ki+0.4527,Kd>0.1344Ki-0.9377,Kd<0.0613Ki+1.7176,Kd>0.0318Ki-2.5853,Kd<0.0197Ki+3.3906,Kd>0.0132Ki-4.2097.
The boundaries of these regions are illustrated in Figure 5.
In this case, the stability region is defined by only two boundaries:(3.35)Kd=m1Ki+b1,Kd=m2Ki+b2,Kd=m3Ki+b3,
because we have the following inequalities:
(3.36)bj<bj+2,forj=2,4,6,…,bj>bj+2,forj=3,5,7,…,b1<b3,mj<mj+2,forj≥1.Figure 5
Region boundaries of Example3.3.Figure 6
Plots ofδi(z) (dotted line), (-a1(z/L))cos(z) (thick solid line) and sin(z)(z2/L2-a0) (thin solid line).As pointed out in Examples3.2 and 3.3, the inequalities given by (3.22) represent half planes in the space of Ki and Kd. Their boundaries are given by lines with the following equations:(3.37)Kd=mjKi+bj,forj=1,2,3,….
We notice that the difference between Examples 3.2 and 3.3 resides in change of the sign of (b1-b3) which is explained by the localization of z1 and z3 in different intervals according to (3.12).We now state an important technical lemma that allows us to develop an algorithm for solving the PID stabilization problem. This lemma shows the behavior of the parameterbj, j=2,3,4,… for different values of the parameter Kp inside the range proposed by Theorem 3.1.Lemma 3.4.
If(3.38)-a0K<Kp<1K(a1αLsin(α)-cos(α)(a0-α2L2)),
where α is the solution of the equation tan(α)=α(2+a1L)/(α2-a1L-a0L2) in the interval [0,π].Then:(1)
bj<bj+2forj=2,4,6,…,(2)
bj>bj+2forj=3,5,7,…,(3)
mj>mj+2forj≥1,andlimj→∞mj=0.Proof.
From (3.18) we have
(3.39)m(zj)=L2zj2,0<zj<zj+2⇒L2zj2>L2zj+22,forj≥1,
then
(3.40)m(zj)>m(zj+2),limj→∞m(zj)=0.Remark 3.5.
As we can see from Figure1, the odd roots of δi(z), that is, zj for j=3,5,7,…, verify zj∈[(2j-3)π/2,(j-1)π] and are getting closer to (2j-3)π/2 as j increase, so we have:
(3.41)limj→∞cos(zj)=0,limj→∞sin(zj)=-1.
Moreover, since cosine function and sine function are monotonically increasing between (2j-3)π/2 and (j-1)π, we have
(3.42){cos(z3)>cos(z5)>cos(z7)>⋯sin(z3)>sin(z5)>sin(z7)>⋯⇒{cos(zj)>cos(zj+2)>0,sin(zj+2)<sin(zj)<0.Remark 3.6.
As we can see from Figure1, the even roots of δi(z) that is, zj for j=2,4,6,… verify zj∈[(2j-3)π/2,(j-1)π] and are getting closer to (2j-3)π/2 as j increase, so we have:
(3.43)limj→∞cos(zj)=0,limj→∞sin(zj)=1.
Moreover, since cosine function and sine function are monotonically decreasing between (2j-3)π/2 and (j-1)π, we have:
(3.44){cos(z2)<cos(z4)<cos(z6)<⋯sin(z2)<sin(z4)<sin(z6)<⋯⇒{cos(zj)<cos(zj+2)<0,0<sin(zj)<sin(zj+2).The change ofb(zj) can be found with graphical approach. Given a stabilizing Kp value inside the admissible range by using Theorem 3.7, the curve of δi(z) is plotted to obtain their roots denoted by zj,j=1,2,3,…, the curves of (-a1(z/L))cos(z) and sin(z)(z2/L2-a0) are plotted too in order to understand the behavior of b(zj).As far as the odd roots ofδi(z) are concerned, the corresponding sin(z)(z2/L2-a0) is decreasing by large magnitude, and for the even ones, the corresponding sin(z)(z2/L2-a0) is increasing by large magnitude. However, compared to the change of sin(zj)(zj2/L2-a0), the difference between the values of (-a1(zj+2/L))cos(zj+2) and (-a1(zj/L))cos(zj) is much smaller than both for odd and even j. Thus, the b(zj) have the similar change rules as sin(zj)(zj2/L2-a0),(3.45)bj+2-bj=LK[-a1zj+2Lcos(zj+2)+sin(zj+2)(zj+22L2-a0)]-LK[-a1zjLcos(zj)+sin(zj)(zj2L2-a0)],bj+2-bj=a1L[cos(zj)-cos(zj+2)]+LK[sin(zj+2)zj+2(zj+22L2-a0)-sin(zj)zj(zj2L2-a0)].
Using Remarks 3.5 and 3.6, we can affirm that(i)
bj>bj+2andbj→-∞asj→+∞ for odd values of j,(ii)
bj<bj+2andbj→+∞asj→+∞ for even values of j.We are ready to state the main results of our work.Theorem 3.7.
Under the above assumptions onK, L, a0, and a1, the range of Kp values for which a solution exists to the PID stabilization problem of an open-loop stable plant with transfer function G(s) is given by
(3.46)-a0K<Kp<1K(a1αLsin(α)-cos(α)(a0-α2L2)),
where α is the solution of the equation tan(α)=α(2+a1L)/(α2-a1L-a0L2) in the interval [0,π].ForKp values outside this range, there are no stabilizing PID controllers. The complete stabilizing region is given by(1)
the cross-section of the stabilizing region in the(Ki,Kd) space is the Trapezoid T if b3>b1,(2)
the cross-section of the stabilizing region in the(Ki,Kd) space is the triangle ∆ if b1>b3,The parametersbj,mj,j=1,2,3 necessary for determining the boundaries, can be determined using (3.21) where zj,j=1,2,3 are the positive-real solutions of δi(z) arranged in ascending order of magnitude.In view of Theorem3.7, we propose an algorithm to determine the set of all stabilizing parameters for a second order delay system.An algorithm for determining PID parameters is given as follows
(1)
ChooseKp in the interval suggested by Theorem 3.7 and initialize j=1,(2)
Find the rootszj of δi(z),(3)
Compute the parameterbj,mj associated with the zj for j=1,2,3 founded,(4)
Determine the stability region in the plane(Ki,Kd) using Figure 7 (Theorem 3.7).(5)
Go to step 1.The stabilizing region of(Kp,Ki), (a): b3>b1, (b): b1>b3.
(a)(b)


## 4. Conclusion

In this work, we have proposed an extension of Hermit-Biehler theorem to compute the stability region for second-order delay system controlled by PID controller. The procedure is based first on determining the range of proportional gain valueKp for which a solution to PID stabilization exists. Then, it is shown that for a fixed Kp inside this range, the stabilizing integral Ki and derivative gain Kd values lie inside a region with known shape and boundaries.



---

*Source: 212053-2009-09-27.xml*