# Nyquist



In [3]:

import numpy as np
import sympy as sy
import itertools
import matplotlib.pyplot as plt
import control.matlab as cm
init_printing()
%matplotlib inline



## 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 $$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}}$$

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

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



In [4]:

z,h = sy.symbols('z,h')
eh = sy.exp(-h)
H = ( (h-1+eh)*z + (1-eh-h*eh) )/( z*z - (1+eh)*z + 0.5)

B,A = sy.fraction(H)
print A




z**2 - z*(1 + exp(-h)) + 0.5




In [28]:

Hs = {}
Hs2 = {}
H0 = cm.tf([1.0], [1, 1, 0])
Hs[0] = (H0,[1.0],[1.0, 1.0, 0])
Hs2[0] = (H0,[1.0],[1.0, 1.0, 0])
for hh in [0.01, 0.1, 0.2, 0.4, 0.8, 1]:
Bp = sy.Poly(B.subs(h,hh))
Ap = sy.Poly(A.subs(h,hh))
a = []
for el in Ap.coeffs():
a.append(float(sy.N(el)))
b = []
for el in Bp.coeffs():
b.append(float(sy.N(el)))
Hs[hh] = (cm.tf(b,a,hh), b,a)
Hs2[hh] = (cm.c2d(H0, hh), b,a)




In [33]:

z = np.exp(1j*np.linspace(0.02,np.pi, 800))
s = 1j*np.linspace(0.3,100,800)

zz = Hs[0][0](s);

xy = np.column_stack((np.real(zz),np.imag(zz)))
np.savetxt('nyquistH0150908.out', xy, fmt='%10.5f')

for (hh,HH) in Hs2.iteritems():
if hh > 0:
zz = HH[0](z);
xy = np.column_stack((np.real(zz),np.imag(zz)))
np.savetxt('nyquistHzz%d150908.out' % (100*hh), xy, fmt='%10.5f')

zz = Hs2[0][0](s);
zz2 = Hs2[0.01][0](z);
zz3 = Hs2[0.1][0](z);
zz4 = Hs2[0.2][0](z);
zz5 = Hs2[0.8][0](z);

plt.figure()
plt.plot(np.real(zz), np.imag(zz))
plt.plot(np.real(zz2), np.imag(zz2))
plt.plot(np.real(zz3), np.imag(zz3))
plt.plot(np.real(zz4), np.imag(zz4))
plt.plot(np.real(zz5), np.imag(zz5))
plt.ylim(-4,1)
plt.xlim(-2,1)




Out[33]:

(-2, 1)




In [38]:

Hln2 = H.subs(h, np.log(2))
Hln2 = cm.tf([0.19, 0.15], [1, -1.5, 0.5], np.log(2))
Hln2




Out[38]:

0.19 z + 0.15
-----------------
z^2 - 1.5 z + 0.5

dt = 0.69314718056




In [41]:

cm.rlocus(Hln2)
plt.ylim(-3,3)
plt.xlim(-3,3)




Out[41]:

(-3, 3)




In [ ]:



### Apply Jury's criterion



In [ ]:




In [ ]:




In [ ]:




In [ ]:




In [ ]: