Stability


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Stability of Open Loop System

  • In order for a system $G(s)={N(s) \over D(s)}$ to be stable all of the roots of the characteristic polynomial need to lie in the left-half plane (LHP).
    • The characteristic equation is the denominator of the transfer function.
    • The roots of the characteristic equation are the exact same as the poles of the transfer function.
    • The eigenvalues of matrix $A$ in the equivalent state space representation are the same as the roots of the characteristic polynomial.
    • In order to have a stable system, roots of $G(s)$ must be in LHP.


$$ G(s) = {1 \over s+a} $$


$$ \mathscr{L}^{-1}(G(s)) = e^{-at} u(t) $$



  • When a pole is negative
    • This root exists in the left half plane
    • Transfer function will ultimately die out
    • The system will eventually be at rest (stable)
  • When a pole is positive
    • This root exists in the right half plane
    • Transfer function will blow up into infinity
    • The system is unstable
  • Transfer function of multiple poles


$$ \begin{align*} G(s) &= {1 \over s+1} \centerdot {1 \over s+3} \centerdot {1 \over s-2}\\\\ G(s) &= {A \over s+1} + {B \over s+3} + {C \over s-2} \\\\ \mathscr{L}^{-1}(G(s)) &= Ae^{-t} + Be^{-3t} + Ce^{2t} \end{align*} $$




  • The last one blows up to infinity to make the whole transfer function unstable
  • Conclusion: a single root in the right half plane makes the whole system unstable

2. Routh-Hurwitz Criterion

Calculating the roots of the system for larger than the second-order polynomial becomes time-consuming and possibly even impossible in a closed-form

  • How can we determine the stability of a higher order polynomial without solving for the roots directly?
    • The great thing about the Routh-Hurwitz criterion is that we do not have to solve for the roots of the characteristic equation


      $$G(s) = {1 \over s^4 + 3s^3-5s^2+s+2}$$

    • If all of the signs are not the same, the system is unstable
    • If you build up a transfer function with a series of poles, then the only way to get a negative coefficient is to have at least one pole exists in right-half plane
  • However, we cannot claim that all positive coefficients are still either stable or unstable


    $$ \begin{align*} G(s) &= \frac{1}{s^4 + 2s^3 + 3s^2 + 10s + 8} \\\\ &= {1 \over s^2-s+4} \centerdot {1 \over s+2} \centerdot {1 \over s+1} \end{align*} $$

2.1. Normal Case

Routh array is a table that can be populated with the coefficients of the polynomial with a few simple rules

  • The number of RHP roots of $D(s)$ is equal to the number of sign changes in the first column of the Routh array


$$As^6 + Bs^5 + Cs^4 + Ds^3 + Es^2 + Fs + G$$




  • Determine the number of roots in RHP by counting the number of sign changes
    • We can determine the number of roots in the right-half plane by looking at this first column
    • For example, it changes sign twice which means that there are two roots in the right half plane


$$G(s) = {1 \over s^4 + 2s^3 + 3s^2 + 10s + 8}$$




2.2. Special Case 1

  • A zero in a row with at least one non-zero appearing later in that same row
    • If you are attempting to access stability of the system, you do not need to complete the rest of the table at this point
    • The system is always unstable because completing Routh array will always result in a sign change of the first column



  • If you are interested in the number of roots located in the right half plane, you can complete the table like below
    • You replace that zero with the Greek symbol epsilon $\epsilon >0$




    • When you finish completing the table, you can take the limit as epsilon $\epsilon$ goes to zero




    • You can see that we still have two unstable roots or two roots in the right half plane

2.3. Special Case 2

  • The second special case is when there is an entire row of zeros, not just a single zero in the row


$$D(s) = s^5 + 2s^4 + 6s^3 + 10s^2 + 8s +12$$




  • Auxiliary polynomial $P(s)$: the row directly above the row of zeros


$$6s^2 + 12s^0 \implies P(s) = s^2 + 2 $$




  • Then $P(s)$ is a factor of the original polynomial $D(s)$



$$ \begin{align*} D(s) &= P(s) \cdot R(s) \\\\\\ s^5 + 2s^4 + 6s^3 + 10s^2 + 8s +12 & = \underbrace{(s^2 + 2)}_{\text{marginally stable}}\;\;\underbrace{(s^3 + 2s^2 + 4s + 6)}_{\text{stable}} \end{align*} $$

  • Apply the Routh-Hurwitz criterion again to $R(s)$



In [1]:
%%html
<center><iframe src="https://www.youtube.com/embed/WBCZBOB3LCA?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
In [2]:
%%html
<center><iframe src="https://www.youtube.com/embed/oMmUPvn6lP8?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
In [3]:
%%html
<center><iframe src="https://www.youtube.com/embed/wGC5C_7Yy-s?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

3. Stability with State Space Representation

  • It is useful to start with scalar systems to get some intuition about what is going on


$$\dot{x} = ax \quad \implies x(t) = e^{at} x(0)$$


$$\begin{cases} a>0: \text{unstable}\\ a<0: \text{asymptotically stable}\\ a=0: \text{critically stable} \end{cases}$$


  • From scalars to matrices?


$$ \dot{x} = Ax + Bu \\ y = Cx + Du $$


  • We cannot say that $A>0$, but we can do the next best thing - eigenvalues !


$$Av = \lambda v$$


  • The eigenvalues tell us how the matrix $A$ 'acts' in different directions (eigenvectors)


$\begin{cases} \text{Re}\,(\lambda)>0: \text{unstable}\\ \text{Re}\,(\lambda)<0: \text{asymptotically stable}\\ \text{Re}\,(\lambda) \leq 0: \text{critically stable} \end{cases}$

4. Stability of Closed Loop System

  • We are interested in the stability of a closed loop system from an open loop system.
    • Stability of a closed loop system can also be determined from the open loop frequency response

4.1. Root Locus (Stability in Time)

  • The closed-loop system is





$$H(s) = \frac{KG}{1+KG}$$


  • A pole exists when the characteristic polynomial in the denominator becomes zero.


$$1+KG(s) = 0 \implies KG(s) = -1 = 1 \angle (2k+1) \pi, \quad k = 0, \pm1, \pm 2, \cdots$$

  • A value of $s^*$ is a closed loop pole if


$$\begin{cases} \lvert KG(s^*) \rvert = 1 \quad \implies \quad K = \frac{1}{\lvert G(s^*)\rvert}\\\\ \angle KG(s^*) = (2k+1)\pi \end{cases}$$


  • Closed-loop poles in the LHP indicate stability

    • The closeness of the poles to the RHP indicate how near to instability the system is

4.2. Relative Stability (Stability in Frequency)

  • Suppose the Bode plot of the open-loop transfer function is given.





  • Question: tell the stability of a closed-loop system from the open-loop frequency response





  • At 180$^o$ of phase lag of the loop, the reference and feedback signal are added.
    • If the magnitude of the loop is greater than 1 the error grows exponentially (unstable)





  • Relative stability is indicated by how close the open-loop frequency response is to the point of 180$^o$ of phase lag and a magnitude of 1
  • More specifically,
    • Gain margin is the distance from a magnitude of 1 (0 dB) at the frequency where $\phi = -180^o$ (phase crossover frequency)
    • Phase margin is the distance from a phase of -180$^o$ at the frequency where $M$ = 0 dB (gain crossover frequency)
  • In order to be stable, both gain and phase margin must be positive
  • Gain and phase margins tell how stable the system would be in closed-loop
    • These quantities can be read from the open-loop data
  • What if $K$ proportional controller is implemented?





  • More intuitively,
    • Gain margin indicates how much you can increase the loop gain $K$ before the system goes unstable
    • Phase margin indicates the amount of phase lag (time delay) you can add before the system goes unstable

by Prof. Richard Hill

In [4]:
%%html
<center><iframe src="https://www.youtube.com/embed/Ky8In4J0v4w?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

4.3. Relative Stability in MATLAB

Open-loop transfer function:

$$G(s) = \frac{1}{s^3+2s^2+s} = \frac{1}{s(s+1)^2}$$
In [5]:
G = tf(1,[1 2 1 0]);





















Gm - Gain margin


Pm - Phase margin


Wcg - Gain crossover frequency


Wcp - Phase crossover frequency













In [6]:



    


margin(G)
[Gm,Pm,Wcg,Wcp] = margin(G)



    















    






Gm =

     2


Pm =

   21.3877


Wcg =

     1


Wcp =

    0.6823












    


























In [7]:



    


20*log10(2)



    















    






ans =

    6.0206























Is stable the closed-loop system with a unity negative feedback?













In [8]:



    


step(feedback(G,1,-1))



    















    















    




























Gain margin

$$CG(s) = 2 \times \frac{1}{s^3+2s^2+s}$$











In [9]:



    


GmG = Gm*G;
margin(GmG)



    















    















    


























In [10]:



    


step(feedback(GmG,1,-1))
xlim([0,50])



    















    















    



























$$CG(s) = 4 \times \frac{1}{s^3+2s^2+s}$$











In [11]:



    


GmG = 4*G;
margin(GmG)



    















    















    


























In [12]:



    


step(feedback(GmG,1,-1))
xlim([0,30])
ylim([-100,100])



    















    















    




























Phase margin


    

  • Add more delay
    • 


$$CG(s) = e^{-j \Phi} \times \frac{1}{s^3+2s^2+s}$$











In [13]:



    


PmG = exp(-j*Pm/180*pi)*G
margin(PmG)



    















    






PmG =
 
  (0.9311-0.3647i)
  ----------------
  s^3 + 2 s^2 + s
 
Continuous-time transfer function.












    


























In [14]:



    


FG = feedback(PmG,1,-1);

t = linspace(0,100,1000);
u = ones(size(t));

[y, tout] = lsim(FG,u,t,0);
plot(tout, abs(y))



    















    















    




























5. Stability in Nyquist Plot



    

  • The gain margin, $K_m = \frac{1}{\lvert G(j\omega) \rvert}$ when $\angle G(j \omega) = 180^o$
      

    • $K_m$ is the maximum stable gain in closed loop
      • 

    • It is easy to find the maximum stable gain from the Nyquist plot
      • 

    

    • 



    

  • The phase margin, $\Phi_m$ is the uniform phase change required to destabilize the system under unitary feedback
    • 















$$G(s) = \frac{20}{(s+1)(s+2)(s+3)}$$











In [15]:



    


s = tf('s');
G = 20/((s+1)*(s+2)*(s+3));



    















    
























In [16]:



    


margin(G)
[Gm,Pm,Wcg,Wcp] = margin(G)



    















    






Gm =

    3.0000


Pm =

   44.4630


Wcg =

    3.3166


Wcp =

    1.8382












    


























In [17]:



    


nyquist(G)

axis([-6,12,-7,7])
axis equal



    















    















    


























In [18]:



    


nyquist(Gm*G)

axis([-6,12,-7,7])
axis equal



    















    















    


























In [19]:



    


% expect instability for large Km

rlocus(G)

axis equal



    















    















    


























In [20]:



    


PmG = exp(-j*Pm/180*pi)*G

margin(PmG)



    















    






PmG =
 
      (14.27-14.01i)
  ----------------------
  s^3 + 6 s^2 + 11 s + 6
 
Continuous-time transfer function.












    


























In [21]:



    


nyquist(PmG)

axis([-6,12,-7,7])
axis equal



    















    















    


























In [22]:



    


%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')