In [2]:

    

import numpy as np
import sympy as sy
import control.matlab as cm


    

In [21]:

    

z = sy.symbols('z', real=False)
hh,r1,s0,s1, aa = sy.symbols('h,r1,s0,s1, a')
pc1 = -1.555-1j*1.555
pc2 = np.conjugate(pc1)
Tr = 1
omega0 = 2.2/Tr
#h = 0.2/omega0
h = Tr/10.0
a = -2*omega0
ad = sy.exp(h*a)
#ad = sy.symbols('a_d')
A2p = sy.simplify(sy.expand(sy.poly((z-np.exp(h*pc1))*(z-np.exp(h*pc2)), z))) # Desired closed loop poles
Acp = sy.simplify(sy.expand(sy.poly((z-np.exp(h*pc1))*(z-np.exp(h*pc2))*(z - ad), z))) # Desired charact polynomial 
Ap = sy.poly((z-1)**2, z) # Plant denominator, double integrator
Bp = sy.poly(h**2/2*(z+1), z)
Rp = sy.poly(z+r1, z)
Sp = sy.poly(s0*z + s1, z)
dioph=(Ap*Rp+Bp*Sp-Acp).all_coeffs()
print A2p
print Acp
print dioph


    

    




Poly(1.0*z**2 - 1.69131784560334*z + 0.732713875869332, z, domain='RR')
Poly(1.0*z**3 - 2.33535426668648*z**2 + 1.82198416806576*z - 0.471894422292842, z, domain='RR')
[1.0*r1 + 0.005*s0 + 0.335354266686485, -2.0*r1 + 0.005*s0 + 0.005*s1 - 0.821984168065758, 1.0*r1 + 0.005*s1 + 0.471894422292842]

In [22]:

    

print sy.im(sy.exp(h*pc1))
print z-sy.exp(h*a)
print z-np.exp(h*pc1)


    

    




-0.132570215943283
z - 0.644036421083141
z - 0.845658922801672 + 0.132570215943283*I

In [13]:

    

sol=sy.solve(dioph, (r1,s0,s1))
print sol[r1]
print sol[s0]
print sol[s1]

t0 = A2p.evalf(subs={z:1})/Bp.evalf(subs={z:1,})
print t0


    

    




-0.856007930368168*a_d + 0.143992069631831
-28.7984139263664*a_d + 32.9380169529651
24.6588108997672*a_d - 28.7984139263663
4.13960302659890

In [19]:

    

print A2p
1-1.69+0.737


    

    




Poly(1.0*z**2 - 1.69131784560334*z + 0.732713875869332, z, domain='RR')






    
Out[19]:





0.04700000000000004

In [15]:

    

G = Km * cm.tf([1], [tau, 1, 0])
Gd = Km * cm.tf([tau*(hpt-1+np.exp(-hpt)), tau*(1-(1+hpt)*np.exp(-hpt))], [1, -(1+np.exp(-hpt)), np.exp(-hpt)], h)
Gd2 = cm.c2d(G, h)
print Gd
print Gd2


    

    




 0.001339 z + 0.001247
----------------------
z^2 - 1.807 z + 0.8071

dt = 0.03


 0.001339 z + 0.001247
----------------------
z^2 - 1.807 z + 0.8071

dt = 0.03

In [33]:

    

print A2p
print A2p.evalf(subs={z:1})
print Bp
print Bp.evalf(subs={z:1})


    

    




Poly(z**2 - 1.58555290309441*z - (-0.792776451547206 + 0.168974310731771*I)*(0.792776451547206 + 0.168974310731771*I), z, domain='EX')
0.0714939167206446
Poly(0.00133942860759726*z + 0.00124712240506047, z, domain='RR')
0.00258655101265772

In [4]:

    

0.3/(5*np.sqrt(2))


    

    
Out[4]:





0.042426406871192847

In [6]:

    

np.exp(-0.21)*np.sin(0.21)


    

    
Out[6]:





0.16897431073177133

In [7]:

    

np.exp(0.03*(-14))


    

    
Out[7]:





0.65704681981505675

In [ ]: