In [90]:
import numpy as np
from sympy import *
import itertools
init_printing()
%matplotlib inline
In [107]:
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)
In [67]:
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))
In [83]:
-0.4645+0.3945**2/0.4645
Out[83]:
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}
In [98]:
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))
Out[98]:
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}
In [38]:
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
In [75]:
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
In [84]:
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))
In [77]:
from sympy.solvers.inequalities import solve_rational_inequalities
solve_rational_inequalities([[ ( (Poly(tab[2][0], K),Poly(tab[2][0], K)), '>=') ]])
Out[77]:
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}
In [46]:
a = [Integer(2), -4+h**2*K, 2+h**2*K]
tab= []
juryrec(a,tab)
for line in tab:
print line
In [47]:
line
Out[47]:
In [48]:
e3 = line[0].subs(h,1)
e3
Out[48]:
In [49]:
solve(e3,K)
Out[49]:
In [50]:
e1 = tab[2][0].subs(h,1)
In [51]:
solve(e1,K)
Out[51]:
In [52]:
e1
Out[52]:
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}
In [97]:
2/np.log(2)
Out[97]:
In [105]:
np.arccos(0.75)
Out[105]:
In [106]:
0.72/np.log(2)
Out[106]:
In [109]:
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)
Out[109]:
In [110]:
tab[2][1]/tab[2][0]
Out[110]:
In [ ]:
In [ ]: