# Jury's stability criterion


```python
import numpy as np
from sympy import *
import itertools
init_printing()
%matplotlib inline
```


```python
def jury(a):
    n = len(a)
    tab = []
    for i in range(n):
        line1 = a[:n-i]
        line2 = line1[::-1]
        tab.append(line1)
        tab.append(line2)
        alpha = line1[-1]/line2[-1]
        a[:n-i] = line1-alpha*line2
    for line in tab:
        print line

def juryrec(a,tab):
    n = len(a)
    if n==1:
        tab.append(a)
    else:
        line1 = a
        line2 = line1[::-1]
        tab.append(line1)
        tab.append(line2)
        alpha = line1[-1]/line2[-1]
        aa = [el1 - alpha*el2 for (el1,el2) in itertools.izip(line1,line2)]
        juryrec(aa[:-1],tab)

        
```


```python
aa1 = np.array([1,-2,1.4,-0.1])
aa2 = np.array([4,-4,-7,2])
a2 = [Integer(4),Integer(-4),Integer(-7),Integer(2)]

tab=[]
juryrec(a2, tab)
print "Coefficients"
print aa2
for line in tab:
    print line
np.roots(aa2)
tab=[]
juryrec(aa1, tab)
print "Coefficients"
print aa1
for line in tab:
    print line
print "Roots"
print np.roots(aa1)
print "Magnitudes"
print np.abs(np.roots(aa1))
```

    Coefficients
    [ 4 -4 -7  2]
    [4, -4, -7, 2]
    [2, -7, -4, 4]
    [3, -1/2, -5]
    [-5, -1/2, 3]
    [-16/3, -4/3]
    [-4/3, -16/3]
    [-5]
    Coefficients
    [ 1.  -2.   1.4 -0.1]
    [ 1.  -2.   1.4 -0.1]
    [-0.1  1.4 -2.   1. ]
    [0.98999999999999999, -1.8600000000000001, 1.2]
    [1.2, -1.8600000000000001, 0.98999999999999999]
    [-0.46454545454545459, 0.39454545454545475]
    [0.39454545454545475, -0.46454545454545459]
    [-0.12945205479452027]
    Roots
    [ 0.95986894+0.56970991j  0.95986894-0.56970991j  0.08026212+0.j        ]
    Magnitudes
    [ 1.11620687  1.11620687  0.08026212]



```python
-0.4645+0.3945**2/0.4645
```




    -0.12945102260495156



## Proportional control of second order system (DC-motor with zero-cancellation)
The transfer function is given by
\begin{equation} 
H(z) = \frac{1}{z^2 - 1.5z + 0.5} = \frac{B(z)}{A(z)}.
\end{equation}
The closed loop pulse-transfer function becomes
\begin{equation} 
H_c(z) = \frac{K H(z)}{KH(z) + 1} = \frac{K B(z)}{A(z) + KB(z)}.
\end{equation}
The characteristic equation of the closed loop system is thus
\begin{equation}
z^2  -1.5z + 0.5+K 
\end{equation}


```python
K = symbols('K')
a = [Integer(1), -1.5, 0.5+K]
tab=[]
juryrec(a,tab)
for line in tab:
    print line
crit1 = simplify(tab[2][0])
crit2 = simplify(tab[4][0])
print crit1
print solve(crit1, K)
print crit2
print solve(crit2, K)
plot(crit2, xlim=(-1.0, 5.0), ylim=(-1,1))
plot(crit1, xlim=(-1.0, 5.0), ylim=(-1,1))
```

## Proportional control of the normalized DC-motor
Zero-order hold sampling of the DC motor with transfer function $G(s)=\frac{1}{s(s+1)}$ gives the discrete time system
\begin{equation} H(z) = \frac{\big(h-1+e^{-h}\big)z + \big(1-e^{-h}-he^{-h}\big)}{z^2 -\big(1+e^{-h}\big)z + e^{-h}} \end{equation}

Let $h=\ln 2 \approx 0.693$. This gives the pulse-transfer function
\begin{equation} H(z) = \frac{B(z)}{A(z)} = \frac{0.19z + 0.15}{z^2 - 1.5z + 0.5} \end{equation}

Proportional control gives the closed loop system 
\begin{equation} 
H_c(z) = \frac{K H(z)}{KH(z) + 1} = \frac{K B(z)}{A(z) + KB(z)}.
\end{equation}
The characteristic equation of the closed loop system is
\begin{equation}
z^2 + (-1.5+0.19K)z + 0.5+0.15K 
\end{equation}




```python
z,h = symbols('z,h')
eh = exp(-h)
H = ( (h-1+eh)*z + 1-eh-h*eh ) / (z**2 -(1+eh)*z + eh )
print H
H1 = H.subs(h,log(2))
print H1
H2 = N(H1)
print H2
```

    (-h*exp(-h) + z*(h - 1 + exp(-h)) + 1 - exp(-h))/(z**2 - z*(1 + exp(-h)) + exp(-h))
    (z*(-1/2 + log(2)) - log(2)/2 + 1/2)/(z**2 - 3*z/2 + 1/2)
    (0.193147180559945*z + 0.153426409720027)/(z**2 - 1.5*z + 0.5)


### Apply Jury's criterion


```python
K = symbols('K')
a = [Integer(1), -1.5+0.19*K, 0.5+0.15*K]
tab=[]
juryrec(a,tab)
for line in tab:
    print line

```

    [1, 0.19*K - 1.5, 0.15*K + 0.5]
    [0.15*K + 0.5, 0.19*K - 1.5, 1]
    [-(0.15*K + 0.5)**2 + 1, 0.19*K - (0.15*K + 0.5)*(0.19*K - 1.5) - 1.5]
    [0.19*K - (0.15*K + 0.5)*(0.19*K - 1.5) - 1.5, -(0.15*K + 0.5)**2 + 1]
    [-(0.15*K + 0.5)**2 + 1 - (0.19*K - (0.15*K + 0.5)*(0.19*K - 1.5) - 1.5)**2/(-(0.15*K + 0.5)**2 + 1)]



```python
print tab[2][0]
print solve(tab[2][0], K)
print solve(tab[4][0], K)
print N(tab[4][0].subs(K,3))
print N(tab[2][0].subs(K,3))
```

    -(0.15*K + 0.5)**2 + 1
    [-10.0000000000000, 3.33333333333333]
    [0.0, 75.0000000000000]
    0.0753230769230770
    0.0975000000000000



```python
from sympy.solvers.inequalities import solve_rational_inequalities
solve_rational_inequalities([[ ( (Poly(tab[2][0], K),Poly(tab[2][0], K)), '>=') ]])

```




    (-oo, -10) U (-10, 10/3) U (10/3, oo)



## Proportional control of the double integrator
\begin{equation}
H(z) = \frac{B(z)}{A(z)} = \frac{h^2(z+1)}{2(z-1)^2} = \frac{h^2z + h^2}{2z^2 -4z + 2}
\end{equation}

Proportional control gives closed-loop system
\begin{equation}
H_c(z) = \frac{KH(z)}{1+KH(z)} = \frac{KB(z)}{A(z) + KB(z)} 
\end{equation}
with denominator polynomial
\begin{equation}
A(z) + KB(z) = 2z^2 +(-4+Kh^2)z + 2+Kh^2
\end{equation}



```python
a = [Integer(2), -4+h**2*K, 2+h**2*K]
tab= []
juryrec(a,tab)
for line in tab:
    print line

```

    [2, K*h**2 - 4, K*h**2 + 2]
    [K*h**2 + 2, K*h**2 - 4, 2]
    [-(K*h**2/2 + 1)*(K*h**2 + 2) + 2, K*h**2 - (K*h**2/2 + 1)*(K*h**2 - 4) - 4]
    [K*h**2 - (K*h**2/2 + 1)*(K*h**2 - 4) - 4, -(K*h**2/2 + 1)*(K*h**2 + 2) + 2]
    [-(K*h**2/2 + 1)*(K*h**2 + 2) + 2 - (K*h**2 - (K*h**2/2 + 1)*(K*h**2 - 4) - 4)**2/(-(K*h**2/2 + 1)*(K*h**2 + 2) + 2)]



```python
line
```




    [-(K*h**2/2 + 1)*(K*h**2 + 2) + 2 - (K*h**2 - (K*h**2/2 + 1)*(K*h**2 - 4) - 4)**2/(-(K*h**2/2 + 1)*(K*h**2 + 2) + 2)]




```python
e3 = line[0].subs(h,1)
e3
```




    -(K/2 + 1)*(K + 2) + 2 - (K - (K/2 + 1)*(K - 4) - 4)**2/(-(K/2 + 1)*(K + 2) + 2)




```python
solve(e3,K)
```




    []




```python
e1 = tab[2][0].subs(h,1)
```


```python
solve(e1,K)
```




    [-4, 0]




```python
e1
```




    -(K/2 + 1)*(K + 2) + 2



## PD control of the double integrator
The pulse-transfer function for the system obtained by zero-order hold sampling of $G(s)=\frac{1}{s^2}$: 
\begin{equation}
H(z) = \frac{B(z)}{A(z)} = \frac{h^2/2(z+1)}{(z-1)^2} = \frac{h^2}{2}\frac{z + 1}{z^2 -2z + 1}
\end{equation}

PD control (simple form) can be implemented using the filter
\begin{equation}
u(kh) = k_1 e(kh) + k_2 e(kh-h)
\end{equation}
with z-transform
\begin{equation}
R(z)U(z) = zU(z) = zk_1 E(z) + k_2 E(z) = S(z) E(z).
\end{equation}

This gives the closed loop system

\begin{equation}
H_c(z) = \frac{\frac{S(z)}{R(z)} \frac{B(z)}{A(z)} } {1 + \frac{S(z)}{R(z)} \frac{B(z)}{A(z)}} = \frac{S(z)B(z)}{R(z)A(z) + S(z)B(z)} 
\end{equation}
with denominator polynomial
\begin{equation}
R(z)A(z) + S(z)B(z) = z(z^2 -2z + 1) + h^2/2(k_1z + k_2)(z+1) = z^3 -2z^2 + z + h^2/2\big(k_1z^2 + (k_1+k_2)z + k_2\big) = z^3 + (h^2/2k_1 - 2)z^2 + (1 + h^2/2(k_1+k_2)z
\end{equation}



```python
2/np.log(2)
```




    2.8853900817779268




```python
np.arccos(0.75)
```




    0.72273424781341566




```python
0.72/np.log(2)
```




    1.0387404294400537




```python
aa1 = np.array([9,-11.81,3.15])
tab=[]
juryrec(aa1, tab)
print "Coefficients"
print aa1
for line in tab:
    print line
np.roots(aa1)

```

    Coefficients
    [  9.   -11.81   3.15]
    [  9.   -11.81   3.15]
    [  3.15 -11.81   9.  ]
    [7.8975, -7.6765000000000008]
    [-7.6765000000000008, 7.8975]
    [0.43581563786008104]





    array([ 0.93980424,  0.37241798])




```python
tab[2][1]/tab[2][0]
```




    -0.97201646090534988




```python

```


```python

```

eb08f7a4726d07a880e9dad0eeae30579ca07759
ipynb
Jupyter Notebook
discrete-time-systems/notebooks/Jurys-criterion-example.ipynb
alfkjartan/control-computarizado
5b9a3ae67602d131adf0b306f3ffce7a4914bf8e
MIT
2020-12-22T09:57:05.000Z
2020-12-22T09:57:05.000Z
discrete-time-systems/notebooks/Jurys-criterion-example.ipynb
alfkjartan/control-computarizado
5b9a3ae67602d131adf0b306f3ffce7a4914bf8e
MIT
discrete-time-systems/notebooks/Jurys-criterion-example.ipynb
alfkjartan/control-computarizado
5b9a3ae67602d131adf0b306f3ffce7a4914bf8e
MIT
2019-09-25T20:02:23.000Z
2019-09-25T20:02:23.000Z
Qwen/Qwen-72B
1. YES
2. YES