1 Introduction

Time delay systems are often encountered in various engineering systems such as electrical and communication network, chemical process, turbojet engine, nuclear reactor and hydraulic system. Delay is frequently a source of instability, oscillation and poor performance in many dynamic systems. Furthermore, delay makes system analysis and control design much more complex[1, 2]. Many processes in industrial and chemical practice are open-loop unstable processes that are known to be difficult to control, such as in the case of polymerization reactors, bioreactors, continuous stirred tank and reactors, exothermic stirred reactors with back mixing, batch reactors, pump with liquid storage tank and combined feed/effluent heat exchanger with adiabatic exothermic reactor, etc., which are inherently open loop unstable by design. Especially, unstable delay processes make control system design a complex task, which has attracted increased attention in the process control community[3]. Recently, many academic researches have been devoted to developing proportional-integral -derivative (PID) control strategies for unstable delay processes. Today, proportional-integral (PI) and PID controller types are the most widely used control strategies. It is estimated that 98% of control loops in the pulp and paper industries are controlled by PI controllers[4] and that in process control applications, more than 95% of the controllers are of PID type[5]. 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[6]. 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. Based on the gain and phase margins criterion, De Paor and O’Malley[7] proposed PID controllers tuning method for the first order delay unstable plant first order delay plant (FODUP). Quinn and Sanathanan[8] investigated a method for the design of controllers for unstable delayed processes based on model matching in the frequency domain. Shafie and Shenton[9] discussed a graphical technique based on the D-partition method for PID controller tuning for stable and unstable processes. Padma Sree et al.[10] presented PI/PID controllers design for first order unstable delay system by extracting the coefficients of the numerator and denominator of the closed-loop transfer function. Huang and Lin[11] developed optimum PID tuning for unstable first and second order delay plant (FODUP and SODUP) based on the minimization of the integral of absolute error (IAE) criterion and using the least-squares method. Poulin and Pomerleau[12] studied a graphical tuning method of PI and PID controllers for integrating processes and unstable processes (FODUP and SODUP) which is based on the analysis of the open-loop frequency response of the process on the Nichols chart. In [13], the authors applied D-partition method to characterize the stability domain in the space of controller parameters. The stability boundary is reduced to a transcendental equation, and the whole stability domain is drawn in two-dimensional plane by sweeping the remaining parameters. However, this result only provides sufficient condition regarding the size of the time delay for stabilization of first-order unstable processes. Using Nyquist criterion, Xiang et al.[14] solved the stabilization problem of second-order unstable delay process by PID controller. An enhanced two-degrees-of-freedom control strategy for second-order unstable process with time-delay is established in [15]. Shamsuzoha introduced an enhanced PID controller for unstable first[16] and second order delay plant[17].Chen et al.[16] investigated the calculation of set point weighting PID parameter for unstable first-order time-delay systems. In [18], the internal model controller equivalent PID tuning rules are synthesized for a low order unstable delay systems in [18]. In [19], the authors presented an adaptive iterative learning control scheme for a class of strict-feedback nonlinear time-delay systems with unknown nonlinearly parameterised and time-varying disturbed functions of known periods. Chen and Li[20] developed an observer-based adaptive iterative learning control scheme for a class of nonlinear systems with unknown time-varying parameters and unknown time-varying delays. In [21], the delay-dependent robust stabilization and an H control for uncertain stochastic Takagi-Sugeno fuzzy systems with discrete interval and distributed time-varying delays are discussed. Kumar and Mittal[22] proposed a parallel fuzzy proportional plus fuzzy integral plus fuzzy derivative (FP+FI+FD) controller for some complex processes, such as first and second-order processes with delay, inverse response process with and without delay and higher order processes.

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[12]. Since almost all plants encountered in process control contain time-delays, computing the complete set of PID controllers that stabilize time delay system is of considerable importance. An analytical approach is used in [2325] to study the characterization of the stability region of delayed systems controlled via PI/PID. Indeed, by using the Hermite-Biehler theorem applicable to the quasi-polynomials [2527], a characterization of all values of the PI/PID stabilization gains for stable and unstable first order delay system is addressed. However, these results are not applicable to the second order delay system. In [28, 29], the stabilizing problem of PI/PID controller for second order delay system is analyzed and then used to obtain all PI and PID gains that stabilize a first and a second order delay interval systems[30, 31].

In this paper, we employed a version of the Hermite-Biehler theorem applicable to quasi-polynomials to investigate the complete set of stabilizing PI/PID parameters for unstable second order process with time delay. However, all these results mentioned above are centered upon the issues of controller design, synthesis and parameter tuning, which do not deal with the determination of the complete set of stabilizing PID controllers based on stability conditions on unstable time delay systems by PID controllers as shown in our paper. We can also improve our approach by searching optimal PI or PID controller inside the stability region for a given performance criteria (integral of square error (ISE), integral of absolute error (IAE), time multiplied by absolute error (ITAE), time multiplied by square error (ITSE), maximum overshoot, rise time and settling time) by using genetic algorithm with regard to the complexity of the optimization problem as indicated in [28].

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[2325]:

$$\matrix{{\delta (s) = d(s) + {{\rm{e}}^{ - {L_1}s}}{n_1}(s) + } \hfill \cr {{e^{ - {L_2}s}}{n_2}(s) + \cdots + {{\rm{e}}^{ - {L_m}s}}{n_m}(s)} \hfill \cr }$$
(1)

where d(s)and n i (s) are polynomials with real coefficients and L i represent time delays. These characteristic equations are recognized as quasi-polynomials. Under the following assumptions

$$\matrix{{({A_1})\,\deg (d(s)) = n\,{\rm{and}}\,\deg ({n_i}(s)) < n} \hfill \cr {{\rm{for}}\,i = 1,2, \cdots, m} \hfill \cr {({A_2})\,0 < {L_1} < {L_2} < \cdots < {L_m}} \hfill \cr }$$
(2)

one can consider the quasi-polynomials δ* (s) described by

$$\matrix{{{\delta ^*}(s) = {{\rm{e}}^{s{L_m}}}\delta (s) = } \hfill \cr {{{\rm{e}}^{s{L_m}}}d(s) + {{\rm{e}}^{s({L_m} - {L_1})}}{n_1}(s) + } \hfill \cr {{{\rm{e}}^{s({L_m} - {L_2})}}{n_2}(s) + \cdots + {n_m}(s).} \hfill \cr }$$
(3)

The zeros of δ(s)are identical to those of δ* (s)since \({{\rm{e}}^{s{L_m}}}\) 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 infinite roots with positive real parts[2325].

The stability of the system with characteristic equation (1) is equivalent to the condition that all the zeros of δ * (s) must be in the open left half of the complex plane. Hence δ * (s) is Hurwitz or is stable. The following theorem gives a necessary and sufficient condition for the stability of δ*(s).

Theorem 1[2325]. Let 5* (s) be given by (3); so

$${\delta ^*}({\rm{j}}\omega) = {\delta _r}(\omega) + {\rm{j}}{\delta _i}(\omega)$$
(4)

where δ r (ω)and δ i (ω) represent the real and imaginary parts of δ* (jω), respectively. Under conditions A 1 and A 2, δ * (s) is stable if and only if:

  1. 1)

    δ r (w)and δ i(ω) have only simple, real roots and these interlace;

  2. 2)

    \(\delta _i^{\prime}({\omega _0}){\delta _r}({\omega _0}) - {\delta _i}({\omega _0})\delta _r^{\prime}({\omega _0}) > 0\)for some ω 0 in [∞,+∞]. where \(\delta _i^{\prime}(\omega)\)and \(\delta _r^{\prime}(\omega)\) denote the first derivative with respect to ω of δ i (ω)and δ r (ω), respectively.

A crucial stage in the application of the preceding 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[2325]. Let M 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 − 2 + η< ω < 2lπ + η, l = l 0,l 0 + 1,l 0 + 2 …, δ r (ω) or δ i (ω)have exactly 4lN + M real roots for a sufficiently large integer l 0.

3 PI control for unstable second order delay system

A second order system with time delay can be mathematically expressed by a transfer function having the following form

$$G(s) = {K \over {{s^2} + {a_1}s + {a_0}}}{{\rm{e}}^{ - Ls}}$$
(5)

where K is the static gain of the plant, L is the time delay, and a 0 and a 1 are the plant parameters. The PI controller is described by the following transfer function

$$C(s) = {K_p} + {{{K_i}} \over s}.$$

According to [25], the closed loop characteristic equation of the system is stable at L = 0, then it will be stable for all L> 0, whereas if it is unstable for L = 0, then it will be unstable for all L> 0. So, a minimal requirement for any control design is that the delay-free closed-loop system be stable.

The characteristic equation of the closed-loop delay-free-system (L =0) is given by

$$\matrix{{\delta (s) = K({K_i} + {K_p}s) + ({s^2} + {a_1}s + {a_0})s = } \hfill \cr {{s^3} + {a_1}{s^2} + ({a_0} + K{K_p})s + K{K_i}.} \hfill \cr }$$
(6)

δ(s) is a third-order polynomial. The closed-loop stability is equivalent to all its coefficients are non zero and have the same sign, we have

$${a_1} > 0,\,{a_0} + K{K_p} > 0,\,\,K{K_i} > 0.$$

For a 0 > 0and a 1 > 0, we obtain an open-loop stable plant. And for a 0 < 0and a 1 > 0, we have an open loop unstable delay plant. In the following, we present the synthesis of PI controller in the case of an unstable second order delay system where K> 0, L> 0, a 0 < 0and a 1 > 0.

We state Theorem 3 determining the range of K p and K i for unstable second order delay system controlled by PI controller.

Theorem 3. Under the assumptions of K> 0, L> 0, a 0 < 0and a 1 > 0, the K p values, for which there is a solution for the stabilization problem of the PI controller of unstable second order delay system, verify:

$$ - {{{a_0}} \over K} < {K_p} < {1 \over K}\left({{a_1}{\alpha \over L}\sin (\alpha) - \cos (\alpha)\left({{a_0} - {{{\alpha ^2}} \over {{L^2}}}} \right)} \right)$$

where α is the solution of the equation: tan(α) = \({{\alpha (2 + {a_1}L)} \over {({\alpha ^2} - {a_1}L - {a_0}{L^2})}}\) in the interval [0, π].

For K p values outside this range, there are no stabilizing PI controllers. The range of K i value is given by

$$0 < {K_i} < \mathop {\min }\limits_{j = 1,3,5, \ldots } \{ {a_j}\} $$

where \({a_j} = a({z_j}) = {{{z_j}} \over {KL}}\left[ {\sin ({z_j})({a_0} - {{z_j^2} \over {{L^2}}}) + {a_1}{{{z_j}} \over L}\cos ({z_j})} \right]\) and Z j ,j = 1, 2, 3, …, are the roots (arranged in ascending order of magnitude) of \({\delta _i}(z) = {\textstyle{z \over L}}(K{K_p} + \cos (z)({a_0} - {{z2} \over {{L^2}}}) - {a_1}{\textstyle{z \over L}}\sin (z))\).

Proof. The characteristic equation of the closed-loop system is given by

$$\delta (s) = K({K_i} + {K_p}s){{\rm{e}}^{ - Ls}} + ({s^2} + {a_1}s + {a_0})s$$
(7)

we deduce the quasi-polynomial δ* (s):

$${\delta ^*}(s) = {{\rm{e}}^{Ls}}\delta (s) = K({K_i} + {K_p}s) + s({s^2} + {a_1}s + {a_0}){{\rm{e}}^{Ls}}$$
(8)

by replacing s by jω,we get

$${\delta ^*}({\rm{j}}\omega) = {\delta _r}(\omega) + {\rm{j}}{\delta _i}(\omega)$$
(9)

where

$$\left\{ {\matrix{{{\delta _r}(\omega) = K{K_i} + ({\omega ^3} - {a_0}\omega)\sin (L\omega) - {a_1}{\omega ^2}\cos (L\omega)} \hfill \cr {{\delta _i}(\omega) = \omega \left[ {K{K_p} + ({a_0} - {\omega ^2})\cos (L\omega) - {a_1}\omega \sin (L\omega)} \right].} \hfill \cr } } \right.$$
(10)

Clearly, the parameter K i only affects the real part of δ* (jω) whereas the parameter K p affects the imaginary part. Let z = Lω,we get:

$$\left\{ {\matrix{{{\delta _r}(z) = K{K_i} + \sin (z)\left({{{{z^3}} \over {{L^3}}} - {a_0}{z \over L}} \right) - {a_1}{{{z^2}} \over {{L^2}}}\cos (z)} \hfill \cr {{\delta _i}(z) = {z \over L}(K{K_p} + \cos (z)\left({{a_0} - {{z2} \over {{L^2}}}) - {a_1}{z \over L}\sin (z)} \right).} \hfill \cr } } \right.$$
(11)

The application of the second condition of Theorem 1 led us to

$$E(z_0) = \delta _i^{\prime} (z_0)\delta _r (z_0) - \delta _i (z_0)\delta _r^{\prime} (z_0) > 0$$

from (11), we have

$$\matrix{{\delta _i^\prime (z) = {{K{K_p}} \over L} - \sin (z)\left({{a_0} + {{2{a_1}z} \over {{L^2}}} - {{{z^3}} \over {{L^3}}}} \right) + } \hfill \cr {\cos (z)\left({{{{a_0}} \over L} - {{3{z^2}} \over {{L^3}}} - {a_1}{{{z^2}} \over {{L^2}}}} \right)} \hfill \cr }$$

for z 0 = 0 (a value that annul δ i (z)), we get

$$E({z_0}) = \delta _i^\prime ({z_0}){\delta _r}({z_0}) = ({{K{K_p} + {a_0}} \over L})K{K_i} > 0$$

which implies \({K_p} > - {\textstyle{{{a_0}} \over K}}\) since K > 0 and K i > 0. □

We consider 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:

where

$$\left\{ {\matrix{{f(z) = K{K_p} + \cos (z)({a_0} - {{{z^2}} \over {{L^2}}})} \hfill \cr {g(z) = {a_1}{z \over L}\sin (z).} \hfill \cr } } \right.$$

We notice z 0 = 0 is a trivial root of the imaginary part. The others are difficult to solve analytically. However, this can be completed graphically. Two cases are presented. In each case, the positive real roots of α i (z) will be denoted by z j ,j =1, 2, 3, …, arranged in increasing order of magnitude.

Case 1.\( - {\textstyle{{{a_0}} \over K}} < {K_p} < {K_u}\)

In this case, we graph the curves of f (z)and g(z)which are presented in Fig. 1.

Fig. 1
figure 1

Representation of the curves of f (z)and g(z)(Case: \( - {\textstyle{{{a_0}} \over K}} < {K_p} < {K_u}\)

K u is defined in second case. We notice that for \( - {\textstyle{{{a_0}} \over K}} < {K_p} < {K_u}\),the curve of f (z) intersects the curve of g(z) twice in the interval [0,π]. Also we can see the following properties:

$$\left\{ {\matrix{{{z_1} \in \left[ {0\,,{\pi \over 2}} \right]} \hfill \cr {{z_3} \in \left[ {{{3\,\pi } \over 2}\,,2\,\pi } \right]} \hfill \cr {{z_5} \in \left[ {{{7\,\pi } \over 2}\,,4\,\pi } \right]} \hfill \cr \vdots \hfill \cr } } \right.\quad {\rm{and}}\quad \left\{ {\matrix{{{z_2} \in \left[ {{\pi \over 2},\pi } \right]} \hfill \cr {{z_4} \in \left[ {{{5\pi } \over 2},3\pi } \right]} \hfill \cr {{z_6} \in \left[ {{{9\pi } \over 2},5\pi } \right]} \hfill \cr \vdots \hfill \cr } } \right.$$

and i.e., Z j verifies:

$$\left\{ {\matrix{{{z_1} \in \left[ {0\,,{\pi \over 2}} \right]} \hfill \cr {{\rm{and}}} \hfill \cr {{z_j} \in \left[ {(2j - 3){\pi \over 2},(j - 1)\pi } \right]\,\,\,{\rm{for}}\,j \geqslant {\rm{ 2}}} \hfill \cr } } \right.$$

and we remark: \({z_j} > z = L\sqrt {{a_0}}, \,\,j = 1,2,3,5,7, \cdots \) …. Case 2. K p K u .

Case 2. K p K u .

Fig. 2 represents the case when K p = K u , and K u is the maximal value of K p , so the plot of f (z) intersects the plot of g(z)twice in the interval [0,⩽] (i.e., f (z) is tangent to g(z)).

Fig. 2
figure 2

Representation of the curves of f (z)and g(z)(Case: K p =K u )

The plot in Fig. 3 corresponds to the case when K p >K u and the plot of f(z)doesnotintersect g(z)twice in the interval [0,π].

Fig. 3
figure 3

Representation of the curves of f(z)and g(z)(Case: K p >K u )

Theorem 2 is used to verify that δ i (z) possesses only simple roots. By replacing Ls by s 1 in (8), we rewrite δ* (s) as

$$\matrix{{{\delta ^*}(s) = {\rm{ }}{{\rm{e}}^{Ls}}\delta (s) = {{\rm{e}}^{{s_1}}}\delta ({s_1}) = } \hfill \cr {{{\rm{e}}^{{s_1}}}\left({{{\left({{{{s_1}} \over L}} \right)}^3} + {a_1}{{\left({{{{s_1}} \over L}} \right)}^2} + {a_0}{{{s_1}} \over L}} \right) + K\left({{K_p}{{{s_1}} \over L} + {K_i}} \right).} \hfill \cr }$$
(12)

For this new quasi-polynomial, we see that M = 3 and N = 1, where M and N designate the highest powers of s 1 and es1 which appear in δ*(s). We choose \(\eta = {\textstyle{\pi \over 4}}\) that satisfies the condition given by Theorem 2 as δ r (η) ≠ 0 and δ i (η) ≠ 0. According to Fig. 1, we notice that for \( - {\textstyle{{{a_0}} \over K}} < {K_p} < {K_u},\,\,{\delta _i}(z)\) possesses four roots in the interval \([0,2\pi - {\textstyle{\pi \over 4}}] = [0,{\textstyle{{7\pi } \over 4}}]\) including the root at origin. As δ i (z) is odd function, so it possesses seven roots in \(\left[ { - 2\pi + {\textstyle{\pi \over 4}},2\pi - {\textstyle{\pi \over 4}}} \right] = \left[ { - {\textstyle{{7\pi } \over 4}},{\textstyle{{7\pi } \over 4}}} \right]\) Hence, we can affirm that δ i (z) has 4N + M = 7 roots in \(\left[ { - 2\pi + {\textstyle{\pi \over 4}},2\pi + {\textstyle{\pi \over 4}}} \right] = \left[ { - {\textstyle{{7\pi } \over 4}},{\textstyle{{9\pi } \over 4}}} \right]\) At the end, according to Theorem 2, δ i (z) has only real roots for every K p in \([ - {\textstyle{{{a_0}} \over K}},{K_u}]\) \({K_p}\geqslant {K_u}\) corresponding Figs. 3 and 4, the roots of δ i (z)are not real. We determine the superior bound of K p, accordingtothe definition of K u , if K p = K u , then the curves of f (z)and g(z)are tangentat the point α, which is expressed by

(13)
Fig. 4
figure 4

Plots of δ i (z) and a(z)

Once α is determined, the parameter K u is expressed by

$${K_u} = {1 \over K}\left({{a_1}{\alpha \over L}\sin (\alpha) - \cos (\alpha)\left({{a_0} - {{{\alpha ^2}} \over {{L^2}}}} \right)} \right).$$
(14)

After the determination of the roots of the imaginary part δ i (z), we move to the evaluation of these roots by the real part δ r (z).

$$\matrix{{{\delta _r}(z) = {\rm{ }}K{K_i} + \sin (z)\left({{{{z^3}} \over {{L^3}}} - {a_0}{z \over L}} \right) - {a_1}{{{z^2}} \over {{L^2}}}\cos (z) = } \hfill \cr {K\left[ {{K_i} - a(z)} \right]} \hfill \cr }$$
(15)

where \(a(z) = {\textstyle{z \over {KL}}}\left[ {\sin (z)\left({{a_0} - {{{z^2}} \over {{L^2}}}} \right) + {a_1}{\textstyle{z \over L}}\cos (z)} \right]\).

Let us denote the roots of δ i (z)by z j , j = 1, 2, 3, ….For z 0,we have

$${\delta _r}({z_0}) = K({K_i} - a(0)) = K{K_i} > 0$$
(16)

for z j z 0,where j = 1, 2, 3, …,we get

$${\delta _r}({z_j}) = K({K_i} - a({z_j})) = K({K_i} - {a_j})$$
(17)

where a(z j)= aj.

Interlacing of the roots of δ r (z)and δ i (z)is equivalent to δ r (z 0) > 0(since K i > 0), δ r (z i ) < 0, δ r (z 2) > 0,….

We can use the interlacing property and the fact that δ i (z) has only real roots to reach that δr(z) possesses real roots, too. From the previous equations, we get the following inequalities.

$$\left\{ {\matrix{{{\delta _r}({z_0}) > 0} \hfill \cr {{\delta _r}({z_1}) < 0} \hfill \cr {{\delta _r}({z_2}) > 0} \hfill \cr {{\delta _r}({z_3}) < 0} \hfill \cr {{\delta _r}({z_4}) > 0} \hfill \cr \vdots \hfill \cr } } \right.\quad \Rightarrow \quad \left\{ {\matrix{{{K_i} > 0} \hfill \cr {{K_i} < {a_1}} \hfill \cr {{K_i} > {a_2}} \hfill \cr {{K_i} < {a_3}} \hfill \cr {{K_i} > {a_4}} \hfill \cr \vdots \hfill \cr } } \right..$$
(18)

From these inequalities, it is clear that the a j odd bounds must be strictly positive. However, a j even bounds are negative in order to find a feasible range of K i ,from which we have:

$$0 < {K_i} < \mathop {\min }\limits_{j = 1,3,5 \cdots } \{ {a_j}\}.$$
(19)

In the following, we are interested to prove that the odd \(a_j^{\prime}{\rm{s}}\) are strictly positive and that the even \(a_j^{\prime}{\rm{s}}\) are negative in order to affirm (19).

The change of a(z j ) can be found with a graphical approach. Given a stabilizing K p value inside the admissible range by using Theorem 3, the curve of δ i (z) is plotted to obtain their roots denoted by z j ,j =1, 2, 3 …,the curve of a(zj) is also plotted in order to understand its behavior (Fig. 4). We note that for every K p in \([ - {\textstyle{{{a_0}} \over K}},{K_u}]\), the parameters a(z j ) verify the following conditions:

$$\left\{ {\matrix{{a({z_j}) > 0,\,{\rm{for}}\,j = 1,3,5,7, \cdots } \hfill \cr {a({z_j}) < 0,\,{\rm{for}}\,j = 2,4,6,8, \cdots.} \hfill \cr } } \right.$$
(20)

Algorithm 1. For determining PI parameters for unstable second order delay system:

  1. 1)

    Choose K p in the interval suggested by Theorem 3;

  2. 2)

    Find the roots z j of δ i (z);

  3. 3)

    Compute the parameter a j associated with the z j previously founded;

  4. 4)

    Determine the lower and the upper bounds for K i as

    $$0 < {K_i} < \mathop {\min }\limits_{j = 1,3,5 \cdots } \{ {a_j}\}.$$
  5. 5)

    Go to Step 1).

Example 1. We consider an unstable second order delay system described by the following transfer function

$$G(s) = {{2{{\rm{e}}^{ - 0.5s}}} \over { - 0.5 + 5s + {s^2}}}.$$

In order to determine K p values, we look for α in interval [0,π] satisfying tan\((\alpha) = {{4.5\alpha } \over {{\alpha ^2} - 2.625}} \Rightarrow \alpha = 1.5617\). K p range is given by 0.25 <K p < 7.85. We remark

$$\left\{ {\matrix{{0.25 < {K_p} < 6.2\,\,\,\,\,\,\, \Rightarrow \,\,\,{K_i} > 0} \hfill \cr {6.2 < {K_p} < 7.85\,\,\,\,\,\,\, \Rightarrow \,\,\,{K_i} < 0.} \hfill \cr } } \right.$$

In consequence, we choose the final range of K p as 0.25 < K p < 6.2.

The controller stability region, obtained in (K p ,K i ) plane, is presented in Fig. 5.

Fig. 5
figure 5

Controller stability domain in (K p ,K i ) for unstable second order delay plant

The step response of the closed-loop system with some PI controllers is shown in Fig. 6.

Fig. 6
figure 6

Time response of the closed-loop system for Example 1

From Fig. 6, we see that the closed-loop system is stable and the output y(t) tracks the step input signal.

4 PID control for unstable second order delay system

The PID controller is described by the following transfer function:

$$C(s) = {K_p} + {{{K_i}} \over s} + {K_d}s.$$

The characteristic equation of the closed-loop delay-free-system (L =0) is given by

$$\matrix{{\delta (s) = K({K_i} + {K_p}s + {K_d}{s^2}) + ({s^2} + {a_1}s + {a_0})s = } \hfill \cr {{s^3} + ({a_1} + K{K_d}){s^2} + ({a_0} + K{K_p})s + K{K_i}.} \hfill \cr }$$

Since this is a third-order polynomial, closed-loop stability is equivalent to all its coefficients are non zero and have the same sign, we have

$${a_1} + K{K_d} > 0,\,\,\,\,{a_0} + K{K_p} > 0,\,\,\,\,K{K_i} > 0.$$

We assume that static gain K of the plant is positive, we obtain the following condition for closed-loop stability of delay-free system

$${K_p} > - {{{a_0}} \over K},\quad {K_i} > 0,\quad {K_d} > - {{{a_1}} \over K}.$$

For a 0 > 0and a 1 > 0, we obtain an open-loop stable plant. And for a 0 < 0 or/and a 1 < 0, we have an open-loop unstable delay plant.

In the following, we present the synthesis of PID controller for the case of an unstable second order delay system where K> 0, L> 0, a 0 < 0 or/and a 1 < 0.

Theorem 4. Under the above assumptions of K> 0, L> 0, a 0 < 0 or/and a 1 < 0, the range of K p values for which a solution exists to the PID stabilization problem of an open-loop unstable plant with transfer function G(s)is given by

$$ - {{{a_0}} \over K} < {K_p} < {1 \over K}\left({{a_1}{\alpha \over L}\sin (\alpha) - \cos (\alpha)\left({{a_0} - {{{\alpha ^2}} \over {{L^2}}}} \right)} \right)$$

where α is the solution of the equation \({\rm{tan}}(\alpha) = {{\alpha (2 + {a_1}L)} \over {({\alpha ^2} - {a_1}L - {a_0}{L^2})}}\,{\rm{in}}\,[0,\pi ]\)

For K p values outside this range, there are no stabilizing PID controllers. The complete stabilizing region given by the cross-section of the stabilizing region in the (K i ,K d ) space is the triangle Δ (Fig. 7).

Fig. 7
figure 7

The stabilizing region of (K i ,K d )

The parameters b j ,m j ,j = 1, 2, necessary for determining the boundaries can be determined using

$$\left\{ {\matrix{{{m_j} = m({z_j}) = {{{L^2}} \over {{z^2}}}} \hfill \cr {{b_j} = b({z_j}) = {L \over {K{z_j}}}\left[ { - {a_1}{{{z_j}} \over L}\cos ({z_j}) + \sin ({z_j})\left({{{z_j^2} \over {{L^2}}} - {a_0}} \right)} \right]} \hfill \cr } } \right.$$

where z j ,j =1, 2, are the positive-real roots of δ i (z)arranged in ascending order of magnitude where δ i (z)is expressed by

$${\delta _i}(z) = {z \over L}\left({K{K_p} + \cos (z)\left({{a_0} - {{z2} \over {{L^2}}}} \right) - {a_1}{z \over L}\sin (z)} \right).$$

Proof. The characteristic equation of the closed-loop system is given by

$$\delta (s) = K({K_i} + {K_p}s + {K_d}{s^2}){{\rm{e}}^{ - Ls}} + ({s^2} + {a_1}s + {a_0})s.$$
(21)

The quasi-polynomial δ *(s)is given by

$$\matrix{{{\delta ^*}(s) = {{\rm{e}}^{Ls}}\delta (s) = } \hfill \cr {K({K_i} + {K_p}s + {K_d}{s^2}) + s({s^2} + {a_1}s + {a_0}){{\rm{e}}^{Ls}}.} \hfill \cr }$$
(22)

By replacing s by ,we get

$${\delta ^*}({\rm{j}}\omega) = {\delta _r}(\omega) + {\rm{j}}{\delta _i}(\omega)$$
(23)

with

$$\left\{ {\matrix{{{\delta _r}(\omega) = K{K_i} - K{K_d}\omega + ({\omega ^3} - {a_0}\omega)\sin (L\omega) - } \hfill \cr {{a_1}{\omega ^2}\cos (L\omega)} \hfill \cr {{\delta _i}(\omega) = \omega \left[ {K{K_p} + ({a_0} - {\omega ^2})\cos (L\omega) - {a_1}\omega \sin (L\omega)} \right].} \hfill \cr } } \right.$$

Clearly, the parameters K i and K d only affect the real part of δ* (jω), whereas the parameter K p affects the imaginary part. Letting z = Lω,we get

$$\left\{ {\matrix{{{\delta _r}(z) = K{K_i} - K{K_d}{{{z^2}} \over {{L^2}}} + \sin (z)\left({{{{z^3}} \over {{L^3}}} - {a_0}{z \over L}} \right) - } \hfill \cr {{a_1}{{{z^2}} \over {{L^2}}}\cos (z)} \hfill \cr {{\delta _i}(z) = {z \over L}\left({K{K_p} + \cos (z)\left({{a_0} - {{z2} \over {{L^2}}}} \right) - {a_1}{z \over L}\sin (z)} \right).} \hfill \cr } } \right.$$
(24)

We notice that δ i (z) is the same in (10) and (24). Consequently, we obtain the same range of K p values, which is given by Theorem 3 that stabilize an unstable second order delay system.

After determination of the roots of the imaginary part δ i (z), we evaluate of these roots by the real part δ r (z).

$$\matrix{{{\delta _r}(z) = K{K_i} - K{K_d}{{{z^2}} \over {{L^2}}} + \sin (z)\left({{{{z^3}} \over {{L^3}}} - } \right.} \hfill \cr {\left. {{a_0}{z \over L}} \right) - {a_1}{{{z^2}} \over {{L^2}}}\cos (z) = } \hfill \cr {K{{{z^2}} \over {{L^2}}}\left\{ { - {K_d} + {K_i}{{{L^2}} \over {{z^2}}} + {L \over {Kz}}\left[ { - {a_1}{z \over L}\cos (z) + } \right.} \right.} \hfill \cr {\left. {\left. {\sin (z)\left({{{{z^2}} \over {{L^2}}} - {a_0}} \right)} \right]} \right\} = } \hfill \cr {K{{{z^2}} \over {{L^2}}}\left[ { - {K_d} + m(z){K_i} + b(z)} \right]} \hfill \cr }$$
(25)

where

$$\left\{ {\matrix{{m(z) = {{{L^2}} \over {{z^2}}}} \hfill \cr {b(z) = {L \over {Kz}}\left[ { - {a_1}{z \over L}\cos (z) + \sin (z)\left({{{{z^2}} \over {{L^2}}} - {a_0}} \right)} \right].} \hfill \cr } } \right.$$
(26)

From condition 1 of Theorem 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).

For z 0 =0, we have

$${\delta _r}({z_0}) = K{K_i} > 0$$
(27)

for z j = z 0,where j =1, 2, 3, …,we get

$$\matrix{{{\delta _r}({z_j}) = K{{z_j^2} \over {{L^2}}}\left[ { - {K_d} + m({z_j}){K_i} + b({z_j})} \right] = } \hfill \cr {K{{z_j^2} \over {{L^2}}}\left[ { - {K_d} + {m_j}{K_i} + {b_j}} \right]} \hfill \cr }$$
(28)

where

$$\left\{ {\matrix{{m({z_j}) = {m_j}} \hfill \cr {b({z_j}) = {b_j}.} \hfill \cr } } \right.$$
(29)

Interlacing the roots of δ r (z)and δ i (z)is equivalent to δ r (z 0) > 0(since K i > 0), δ r (z 1) < 0, δ r (z 2) > 0, ….

We can use the interlacing property and the fact that δ i (z) has only real roots to reach that δ r (z) possesses real roots as well. From the previous equations, we get the following inequalities

(30)

Thus, intersecting all these regions in the (K i K d )space, we get the set of (K i ,K d ) values for which the roots of lowing r (z) and lowing i (z) interlace for a fixed value of K p .We notice that all these regions are half planes with their boundaries being lines with positive slopes m j . □

Example 2. Consider the plant given by (5) with the following parameters K = 2, a 1 =5, a 0 = −0.5 and L = 0.5.

$$G(s) = {{2{{\rm{e}}^{ - 0.5s}}} \over { - 0.5 + 5s + {s^2}}}.$$

The imaginary part δ i (z) has only simple real roots if and only if 0.25 <K p < 7.85. We set the controller parameter K p to 1.5, which is inside the previous range. For this Kp, δ i (z) takes the form:

$$\matrix{{{\delta _i}(z) = {z \over {0.5}}\left({2{K_p} + \cos (z)\left({ - 0.5 - {{{z^2}} \over {0.25}}} \right) - 5{z \over {0.5}}\sin (z)} \right) = } \hfill \cr {2z\left[ {3 + \cos (z)(- 05 - 4{z^2}) - 10z\sin (z)} \right].} \hfill \cr }$$

We next compute some of the positive real roots of this equation and arrange them in increasing order of magnitude:

$$\matrix{{{z_0} = 0,\,\,{z_1} = 0.4375,\,\,{z_2} = 2.292,\,\,{z_3} = 5.185,\,\,{z_4} = 8.141,\,\,} \hfill \cr {{z_5} = 11.22,\,\,{z_6} = 14.31.} \hfill \cr }$$

Using (29), we calculate the parameters m j and b j for j =1,…, 6:

$$\left\{ {\matrix{{{m_1} = 1.3061} \hfill \cr {{m_2} = 0.0476} \hfill \cr {{m_3} = 0.0093} \hfill \cr {{m_4} = 0.0038} \hfill \cr {{m_5} = 0.0020} \hfill \cr {{m_6} = 0.0012} \hfill \cr } } \right.\quad {\rm and}\quad \left\{ {\matrix{{{b_1} = - 1.9581} \hfill \cr {{b_2} = 3.4130} \hfill \cr {{b_3} = - 5.7761} \hfill \cr {{b_4} = 8.5304} \hfill \cr {{b_5} = - 11.5059} \hfill \cr {{b_6} = 14.5353.} \hfill \cr } } \right.$$

Interlacing of the roots of the real and the imaginary part occurs for K p =1.5, if and only if the following inequalities are satisfied:

$$\left\{ {\matrix{{{K_i} > 0} \hfill \cr {{K_d} > {\rm{1}}.{\rm{3061}}{K_i}{\rm{ - 1}}.{\rm{9581}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{0476}}{K_i} + {\rm{3}}.{\rm{413}}} \hfill \cr {{K_d} > {\rm{0}}.{\rm{0093}}{K_i}{\rm{ - 5}}.{\rm{7761}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{0038}}{K_i} + {\rm{8}}.{\rm{5304}}} \hfill \cr {{K_d} > {\rm{0}}.{\rm{002}}{K_i}{\rm{ - 11}}.{\rm{5059}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{0012}}{K_i} + {\rm{14}}.{\rm{5353}}.} \hfill \cr } } \right.$$

The boundaries of these regions are illustrated in Fig. 8.

Fig. 8
figure 8

Region boundaries of Example 2 (K p = 1.5)

Example 3. Consider the plant given by (5) with the following parameters K = 1, a 1 = −0.1, a 0 = 0.5 and L =1.

$$G(s) = {{{{\rm{e}}^{ - s}}} \over {05 - 01s + {s^2}}}.$$

The imaginary part δ i (z) has only simple real roots if and only if −0. 5 <K p < 0. 2314. We set the controller parameter K p to 0.1, which is inside the previous range. For this K p , δ i (z) takes the form

$${\delta _i}(z) = z(01 + \cos (z)(05 - {z^2}) + 01z\sin (z)) = 0.$$

We next compute some of the positive real roots of this equation and arrange them in increasing order of magnitude:

$$\matrix{{{z_0} = 0,\,\,{z_1} = 0.8711,\,\,{z_2} = 1.409,\,{z_3} = 4.695,\,\,{z_4} = 7.839,} \hfill \cr {{z_5} = 10.99,{z_6} = 1413} \hfill \cr }$$

Using (29), we calculate the parameters m j and b j for j = 1, …, 6:

$$\left\{ {\matrix{{{m_1} = 1.3178} \hfill \cr {{m_2} = 0.5037} \hfill \cr {{m_3} = 0.0454} \hfill \cr {{m_4} = 0.0163} \hfill \cr {{m_5} = 0.0083} \hfill \cr {{m_6} = 0.005} \hfill \cr } } \right.\quad {\rm and} \quad \left\{ {\matrix{{{b_1} = 0.2917} \hfill \cr {{b_2} = 1.0565} \hfill \cr {{b_3} = - 4.5895} \hfill \cr {{b_4} = 7.7758} \hfill \cr {{b_5} = - 10.9449} \hfill \cr {{b_6} = 14.095.} \hfill \cr } } \right.$$

Interlacing the roots of the real and the imaginary part occurs for K p = 0.1, if and only if the following inequalities are satisfied:

$$\left\{ {\matrix{{{K_i} > 0} \hfill \cr {{K_d} > {\rm{1}}.{\rm{3178}}{K_i}{\rm{ + 0}}.{\rm{2917}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{5037}}{K_i} + {\rm{1}}.{\rm{0565}}} \hfill \cr {{K_d} > {\rm{0}}.{\rm{0454}}{K_i}{\rm{ - 4}}.{\rm{5895}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{0163}}{K_i} + {\rm{7}}.{\rm{7758}}} \hfill \cr {{K_d} > {\rm{0}}.{\rm{0083}}{K_i}{\rm{ - 10}}.{\rm{9449}}} \hfill \cr {{K_d} < {\rm{0}}.{\rm{005}}{K_i} + {\rm{14}}.{\rm{095}}.} \hfill \cr } } \right.$$

The boundaries of these regions are illustrated in Fig. 9.

Fig. 9
figure 9

Region boundaries of Example 3 (K p = 0.1)

An enlargement of Fig. 9 is shown in Fig. 10.

Fig. 10
figure 10

Stability region of Example 3 (K p =0.1)

Example 4. Consider the plant given by (5) with the following parameters K =1, a 1 = −0.1, a 0 = −0.5and L = 1.

$$G(s) = {{{{\rm{e}}^{ - s}}} \over { - 05 - 01s + {s^2}}}.$$

The imaginary part δ i (z) has only simple real roots if and only if 0. 5 <K p < 0. 7442. We set the controller parameter K p to 0.6, which is inside the previous range. For this K p , δ i (z) takes the form

$$z(- 0.1 + \cos (z)(- 05 - {z^2}) + 0.1z\sin (z)) = 0.$$

We next compute some of the positive real roots of this equation and arrange them in increasing order of magnitude:

$$\matrix{{{z_0} = 0,\,{z_1} = 0.4188,\,{z_2} = 1.1916,\,{z_3} = 4.718,\,{z_4} = 7.8316,} \hfill \cr {{z_5} = 10.9915,\,{z_6} = 14.1271.} \hfill \cr }$$

Using (29), we calculate the parameters m j and b j for j = 1, …, 6:

$$\left\{ {\matrix{{{m_1} = 5.7015} \hfill \cr {{m_2} = 0.7043} \hfill \cr {{m_3} = 0.0449} \hfill \cr {{m_4} = 0.0163} \hfill \cr {{m_5} = 0.0083} \hfill \cr {{m_6} = 0.005} \hfill \cr } } \right.\quad {\rm and}\,\left\{ {\matrix{{{b_1} = 0.7472} \hfill \cr {{b_2} = 1.5338} \hfill \cr {{b_3} = - 4.8233} \hfill \cr {{b_4} = 7.8957} \hfill \cr {{b_5} = - 11.0373} \hfill \cr {{b_6} = 14.1628.} \hfill \cr } } \right.$$

Interlacing the roots of the real and the imaginary part occurs for K p =0. 6, if and only if the following inequalities are satisfied:

$$\left\{ {\matrix{{{K_i} > 0} \hfill \cr {{K_d} > 5.7015{K_i} + 0.7472} \hfill \cr {{K_d} < 0.7043{K_i} + 1.5338} \hfill \cr {{K_d} > 0.0449{K_i} - 4.8233} \hfill \cr {{K_d} < 0.0163{K_i} + 7.8957} \hfill \cr {{K_d} > 0.0083{K_i} - 11.0373} \hfill \cr {{K_d} < 0.005{K_i} + 14.1628.} \hfill \cr } } \right.$$

The boundaries of these regions are illustrated in Fig. 11.

Fig. 11
figure 11

Region boundaries of Example 4 (K p = 0.6)

An enlargement of Fig. 11 is shown in Fig. 12.

Fig. 12
figure 12

Stability region of Example 4 (K p = 0.6)

The stability region is defined by only two boundaries

$$\left\{ {\matrix{{{K_d} = {m_{\rm{1}}}{K_i} + {b_1}} \hfill \cr {{\rm{and}}} \hfill \cr {{K_d} = {m_{\rm{2}}}{K_i} + {b_2}} \hfill \cr } } \right.$$

because we have the following inequalities:

$$\left\{ {\matrix{{{b_j} < {b_{j + 2}},} \hfill & {{\rm{for}}\,\,\,\,j = 2,4,6, \cdots } \hfill \cr {{b_j} > {b_{j + 2}},} \hfill & {{\rm{for}}\,\,\,\,j = 1,3,5,7, \cdots } \hfill \cr {{m_j} < {m_{j + 2}},} \hfill & {{\rm{for}}\,\,\,\,j \geqslant {\rm{1}}{\rm{.}}} \hfill \cr } } \right.$$

As pointed out in Examples 2, 3 and 4, the inequalities given by (30) represent half planes in the space of K i and K d . Their boundaries are given by lines with

$${K_d} = {m_j}{K_i} + {b_j}\,\,\,{\rm{for}}\,\,\,\,j = 1,2,3, \cdots.$$
(31)

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 parameter b j , j = 1, 2, 3, …, for different values of the parameter K p inside the range proposed by Theorem 4.

Lemma 1. If \( - {\textstyle{{{a_0}} \over K}} < {K_p} < {\textstyle{1 \over K}}({a_1}{\textstyle{a \over L}}\sin (\alpha) - \cos (\alpha)({a_0} - {\textstyle{{{a^2}} \over {{L^2}}}}))\), where α is the solution of the equation \(\tan (\alpha) = {{\alpha (2 + {a_1}L)} \over {({\alpha ^2} - {a_1}L - {a_0}{L^2})}}\) in the interval [0,π], then:

$$\left\{ {\matrix{{1 - {b_j} < {b_{j + 2}},\,\,\,\,\,\,{\rm{for}}\,\,\,\,j = 2,4,6, \cdots } \hfill \cr {2 - {b_j} > {b_{j + 2}},\,\,\,\,\,{\rm{for}}\,\,\,\,j = 1,3,5,7, \cdots } \hfill \cr {3 - {m_j} > {m_{j + 2}},\,\,\,\,{\rm{for}}\,\,\,\,j \geqslant 1\,\,\,{\rm{and}}\,\,\,\mathop {\lim }\limits_{j \rightarrow \infty } \,{m_j} = 0} \hfill \cr {4 - {b_1} > - {{{a_1}} \over K}.} \hfill \cr } } \right.$$
(32)

Proof. From equations (26) we have \(m({z_j}) = {{{L^2}} \over {z_j^2}},0 < {z_j} < {z_{j + 2}} \Rightarrow {{{L^2}} \over {z_j^2}} > {{{L^2}} \over {z_{j + 2}^2}}\) for j⩾1 then rn (z j ) > m(zj+2) and lim j→∞ m (z j )=0 □

Remark 1. As we can see from Fig. 1, the odd roots of δ i (z), i.e., for j = 3, 5, 7, …,verify \({z_j} \in [(2j - 3){\textstyle{\pi \over 2}},(j - 1)\pi ]\) and are getting closer to\((2j - 3){\textstyle{\pi \over 2}}\) as j increases, so we have lim j→∞ cos(z j )=0 and lim j→∞ sin(z j )=−1.

Moreover, since cosine function and sine function are monotonically increasing between \((2j - 3){\textstyle{\pi \over 2}}\) and \((j - 1)\pi \), we have

Remark 2. As we can see from Fig. 1, the even roots of δ i (z), i.e., for j = 2, 4, 6, …, verify \((2j - 3){\textstyle{\pi \over 2}}\) and are getting closer to as j increase, so we have lim j→∞ cos(z j ) = 0 and lim j→∞ sin(z j ) = 1.

Moreover, since cosine function and sine function are monotonically decreasing between \((2j - 3){\textstyle{\pi \over 2}}\) and \((j - 1)\pi \), we have

The change of b(z j ) can be found with graphical approach. Given a stabilizing K p value inside the admissible range by using Theorem 4; the curve of δ i (z) is plotted to obtain its roots denoted by z j =1, 2, 3, …,the curves of \( - {a_1}{\textstyle{z \over L}}\cos (z)\) and \(\sin (z)({{{z^2}} \over {{L^2}}} - {a_0})\) are plotted as well in order to understand the behavior of b(z j ) (Fig. 13).

$$\left\{ {\matrix{{{z_1} \in \left[ {0\,,{\pi \over 2}} \right]} \hfill \cr {{z_3} \in \left[ {{{3\pi } \over 2},2\,\pi } \right]} \hfill \cr {{z_5} \in \left[ {{{7\pi } \over 2},4\,\pi } \right]} \hfill \cr \vdots \hfill \cr } } \right.\quad {\rm{and}}\,\,\left\{ {\matrix{{{z_2} \in \left[ {{\pi \over 2},\pi } \right]} \hfill \cr {{z_4} \in \left[ {{{5\pi } \over 2},3\pi } \right]} \hfill \cr {{z_6} \in \left[ {{{9\pi } \over 2},5\pi } \right]} \hfill \cr \vdots \hfill \cr } } \right.$$
Fig. 13
figure 13

Plots of δ i (z) (dotted line), \(\sin (z)({{{z^2}} \over {{L^2}}} - {a_0})\) (thick solid line) and \( - {a_1}{\textstyle{{{z_j}} \over L}}\cos ({z_j})\) (thin solid line)

As far as the odd roots of δ i (z) are concerned, the corresponding \(\sin (z)\left({{{{z^2}} \over {{L^2}}} - {a_0}} \right)\) is decreasing by large magnitude \(\sin ({z_j})({{z_j^2} \over {{L^2}}} - {a_0}) \rightarrow - \infty \) as \(j \rightarrow + \infty \), and for the even ones, the corresponding \(\sin (z)\left({{{{z^2}} \over {{L^2}}} - {a_0}} \right)\) is increasing by large magnitude \(\sin ({z_j})\left({{{z_j^2} \over {{L^2}}} - {a_0}} \right) \rightarrow + \infty \). However, compared to the change of \(j \rightarrow + \infty \) the difference between the values of \(\sin ({z_j})\left({{{z_j^2} \over {{L^2}}} - {a_0}} \right)\) and \( - {a_1}{{{z_{j + 2}}} \over L}\cos ({z_{j + 2}})\) \( - {a_1}{{{z_j}} \over L}\cos ({z_j})\)is much smaller than both odd and even j. Thus, the b(z j ) has the similar change rules as\(\sin ({z_j})({{z_j^2} \over {{L^2}}} - {a_0})\).

We can affirm that:

  1. 1)

    b j >b j+ 2 and b j →−∞ as j→ + ∞ for odd values of j.

  2. 2)

    b j <b j+ 2 and b j → + ∞ as j → + ∞ for even values of j.

Remark 3. Let us prove that \(b({z_j}) > - {\textstyle{{{a_1}} \over K}}\) for j = 1, 3, 5, 7, …,we have i.e., z j verifies

$$\left\{ {\matrix{{{z_1} \in \left[ {0\,{\pi \over 2}} \right]} \hfill \cr {{\rm{and}}} \hfill \cr {{z_j} \in \left[ {(2j - 3){\pi \over 2},(j - 1)\pi } \right],\,\,\,{\rm{for}}\,j \geqslant {\rm{ 2}}{\rm{.}}} \hfill \cr } } \right.$$

For z 1,,we have cos(z 1) > 0 and sin(z 1) > 0. Case 1. a 0 < 0and a 1 > 0.

We suppose

$$\matrix{{b({z_1}) < - {{{a_1}} \over K} \Rightarrow } \hfill \cr {{L \over {Kz}}\left[ { - {a_1}{{{z_1}} \over L}\cos ({z_1}) + \sin ({z_1})({{z_1^2} \over {{L^2}}} - {a_0})} \right] < - {{{a_1}} \over K} \Rightarrow } \hfill \cr { - {a_1}{{{z_1}} \over L}\cos (z) + \sin ({z_1})({{z_1^2} \over {{L^2}}} - {a_0}) < - {a_1}{{{z_1}} \over L} \Rightarrow } \hfill \cr {{a_1}{{{z_1}} \over L}(1 - \cos ({z_1})) + \sin ({z_1})({{z_1^2} \over {{L^2}}} - {a_0}) < 0} \hfill \cr }$$

which is wrong because we have

$$\left\{ {\matrix{{{a_1}{{{z_1}} \over L}(1 - \cos ({z_1})) > 0} \hfill \cr {{\rm{and}}} \hfill \cr {\sin ({z_1})\left({{{z_1^2} \over {{L^2}}} - {a_0}} \right) > 0} \hfill \cr } } \right.$$

so \(b({z_1}) > - {\textstyle{{{a_1}} \over K}}\)

Case 2. a 0 > 0 and a 1 < 0.

For z 1, we have cos(z 1) > 0, sin(z 1) > 0, and in this case, we notice that \({z_1} < L\sqrt {{a_0}} \).

$$\matrix{{\left[ { - {a_1}{{{z_1}} \over L}\cos ({z_1}) + \sin ({z_1})\left({{{z_1^2} \over {{L^2}}} - {a_0}} \right)} \right] > } \hfill \cr { - {a_1}{{{z_1}} \over L}\cos ({z_1}) > - {a_1}{{{z_1}} \over L} \Rightarrow } \hfill \cr {{L \over {K{z_1}}}\left[ { - {a_1}{{{z_1}} \over L}\cos ({z_1}) + \sin ({z_1})\left({{{z_1^2} \over {{L^2}}} - {a_0}} \right)} \right] > - {{{a_1}} \over K} \Rightarrow } \hfill \cr {b({z_1}) = {b_1} > - {{{a_1}} \over K}.} \hfill \cr }$$

Case 3. a 0 < 0 and a 1< 0.

$$\matrix{{\left[ { - {a_1}{{{z_1}} \over L}\cos ({z_1}) + \sin ({z_1})\left({{{z_1^2} \over {{L^2}}} - {a_0}} \right)} \right] > } \hfill \cr { - {a_1}{{{z_1}} \over L}\cos ({z_1}) > - {a_1}{{{z_1}} \over L} \Rightarrow } \hfill \cr {{L \over {K{z_1}}}\left[ { - {a_1}{{{z_1}} \over L}\cos ({z_1}) + \sin ({z_1})\left({{{z_1^2} \over {{L^2}}} - {a_0}} \right)} \right] > - {{{a_1}} \over K} \Rightarrow } \hfill \cr {b({z_1}) = {b_1} > - {{{a_1}} \over K}.} \hfill \cr }$$

We conclude that

$$b({z_1}) > - {{{a_1}} \over K}.$$
(33)

By using (32) and (33), the stability region defined by (30) can be reduced to the following boundaries:

$$\left\{ {\matrix{{{K_d} > {m_{\rm{1}}}{K_i} + {b_1}} \hfill \cr {{K_d} < {m_{\rm{2}}}{K_i} + {b_2}} \hfill \cr {{K_d} > - {{{a_1}} \over K}.} \hfill \cr } } \right.$$
(34)

In view of Theorem 4, we propose an algorithm to determine the set of all stabilizing PID parameters for an unstable second order delay system.

Algorithm 2. For determining PID parameters for unstable second order delay system:

  1. 1)

    Choose K p in the interval suggested by Theorem 4;

  2. 2)

    Find the roots z j of δ i (z);

  3. 3)

    Compute the parameters m j and b j , j =1, 2, associated with the z j previously found;

  4. 4)

    Determine the stability region in the plane (K i ,K d ) using Fig. 7 (Theorem 4);

  5. 5)

    Go to step 1.

Example 5. Finally, let us use the results of this section to determine the entire set of PID parameters to stabilize the unstable delay system which is presented in Example 2. The range of K p values specified by Theorem 4 is given by

$$ - 0.5 < {K_p} < 0.2314.$$

By sweeping over this range and using the algorithm presented earlier, we obtain the stabilizing set of (K p ,K i ,K d ) values sketched in Fig. 14.

Fig. 14
figure 14

The stabilizing region of (K p ,K i ,K d ) values for the PID controller in Example 2

Example 6. Consider an unstable process as given by Huang and Chen[32]

$$G(s) = {{{{\rm{e}}^{ - 0.939s}}} \over {(5s - 1)(207s + 1)}} = {{0.0966} \over {{{\rm{e}}^{ - 0.939s}}}}{s^2} + 2.93s - 0.0966.$$

We present a comparison of different methods for design of PID controllers for unstable second order delay system.

By using our approach, the entire set of PID parameters stabilizing G(s) is obtained by the following figure.

We observe that the PID tuning values of some existing methods in Table 1 are included in the stability region given by Fig. 15.

Fig. 15
figure 15

The stabilizing region of (K p ,K i ,K d ) values for the PID controller

Table 1 Performance comparison

1) Method of Huang and Chen[32]: K p = 6.186,K i = 0.8628, K d =9.1058.

These PID tuning parameters are included in the stability region in the (K i ,K d ) plane for a fixed K p = 6.186 (see Fig. 16).

Fig. 16
figure 16

Stability region in (K i , K d )plane fora K p = 6.186

2) Method of Huang and Lin[33]: K p = 3.954,Ki = 0.7975,K d =8.2006.

These PID tuning parameters are included in the stability region in the (K i ,K d ) plane for a fixed K p = 3.954 (see Fig. 17).

Fig. 17
figure 17

Stability region in (K i , K d )plane fora K p = 3.954

3) Method of Poulin and Pomerleau[12]: K p = 3.050,K i = 0.4036,K d = 6.3135.

These PID tuning parameters are included in the stability region in the (K i ,K d ) plane for a fixed K p = 3.050 (see Fig. 18).

Fig. 18
figure 18

Stability region in (K i , K d )plane fora K p = 3.050

4) Method of Lee et al.[34]: K p = 7.144,K i = 1.0688,K d =11.8233.

These PID tuning parameters are included in the stability region in the (K i , K d ) plane for a fixed K p = 7.144 (see Fig. 19).

Fig. 19
figure 19

Stability region in (K i ,K d )plane fora K p = 7.144

5 Conclusion

In this work, we proposed an extension of Hermite-Biehler theorem to compute the region of stability for unstable second order delay system controlled by PI and PID controllers. Firstly, the procedure is based on determining the range of proportional gain value K p for which a solution to PID stabilization problem exists. Then, it is shown that for a fixed K p inside this range, the stabilizing integral K i and derivative gain K d values lie inside a region with known shape and boundaries.