In [2]:

    

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


    

In [22]:

    

z = sy.symbols('z', real=False)
hh,r1,r2,s0,s1,s2,s3, aa = sy.symbols('h,r1,r2, s0,s1,s2,s3, a')

Bp = 1+aa
Ap = z+aa

pd1 = 0.4
pd2 = 0.4 + 1j*0.3
pd3 = np.conjugate(pd2)
Ac1 = sy.simplify( sy.expand( (z-pd2)*(z-pd3) ) )
Ac = sy.simplify( sy.expand( (z-pd1)*Ac1 ) )  # Desired closed loop poles
print Ac

print sy.simplify( sy.expand( Ap * z*z*(z-1) ) )


    

    




z**3 - 1.2*z**2 + 0.57*z - 0.1
z**2*(a*z - a + z**2 - z)

In [23]:

    

Bp = Bp.subs(aa, -0.8)
Ap = Ap.subs(aa, -0.8)
print Bp
print Ap


    

    




0.200000000000000
z - 0.8

In [33]:

    

Rp = (z-1)*(z**2+r1*z + r2)
Sp = s0*z**3 + s1*z**2 + s2*z + s3
Ao = (z-0.2)**3

diophLHS = sy.simplify( sy.collect( sy.expand( Ap *z**2 * Rp) + Bp*Sp, z ) ) 
diophRHS = sy.simplify( sy.collect( sy.expand( Ac * Ao), z) ) 
diophEq = sy.poly( sy.simplify( sy.expand( diophLHS - diophRHS ) ), z )
dioph=diophEq.all_coeffs()
print sy.latex(diophLHS)
print sy.latex(diophRHS)
print dioph
print sy.expand(Ao)


    

    




0.2 s_{2} z + 0.2 s_{3} + z^{6} + z^{5} \left(r_{1} - 1.8\right) + z^{4} \left(- 1.8 r_{1} + r_{2} + 0.8\right) + z^{3} \left(0.8 r_{1} - 1.8 r_{2} + 0.2 s_{0}\right) + z^{2} \left(0.8 r_{2} + 0.2 s_{1}\right)
z^{6} - 1.8 z^{5} + 1.41 z^{4} - 0.594 z^{3} + 0.138 z^{2} - 0.01656 z + 0.0008
[1.0*r1 + 2.22044604925031e-16, -1.8*r1 + 1.0*r2 - 0.61, 0.8*r1 - 1.8*r2 + 0.2*s0 + 0.594, 0.8*r2 + 0.2*s1 - 0.138, 0.2*s2 + 0.01656, 0.2*s3 - 0.0008]
z**3 - 0.6*z**2 + 0.12*z - 0.008

In [27]:

    

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

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


    

    




-2.22044604925031e-16
0.610000000000000
2.52000000000000
-1.75000000000000
-0.0828000000000000
1.35000000000000

In [ ]: