In [1]:
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
import control.matlab as cm
import sympy as sy

omega = 1.0

# h = np.pi/omega
h = 0.6/omega
A = np.array([[0, omega], [-omega, 0]])
B = np.array([[1],[0]])
C = np.array([[1, 0]])
D = np.array([[0]])
sys_c = cm.ss(A,B,C,D)

wh = omega*h
F = np.array([[np.cos(wh), np.sin(wh)], [-np.sin(wh), np.cos(wh)]])
G = 1.0/omega* np.array([[np.sin(wh)], [np.cos(wh)-1]])
sys_d = cm.ss(F,G,C,D, h)

Tc = np.linspace(0,8, 300)
(yc,tc) = cm.step(sys_c, Tc)
Td = h*np.arange(80)
(yd,td) = cm.step(sys_d, Td)
plt.plot(tc,yc)
plt.plot(td,yd[0], '*')
plt.xlim([0,8])


/home/kjartan/.virtualenvs/jupyter_notebook/src/ipython/IPython/kernel/__init__.py:10: UserWarning: The `IPython.kernel` package has been deprecated. You should import from ipython_kernel or jupyter_client instead.
  warn("The `IPython.kernel` package has been deprecated. "
Out[1]:
(0, 8)

In [10]:
sys_d2 = cm.c2d(sys_c, h)
cm.ssdata([sys_d2,sys_d])


Out[10]:
(array([[ 0.]]),
 array([[ 0.,  0.]]),
 array([[ 0.]]),
 matrix([[ A = [[ 0.82533561  0.56464247]
   [-0.56464247  0.82533561]]
 
  B = [[ 0.56464247]
   [-0.17466439]]
 
  C = [[1 0]]
 
  D = [[0]]
 
  dt = 0.6
  ,
          A = [[ 0.82533561  0.56464247]
   [-0.56464247  0.82533561]]
 
  B = [[ 0.56464247]
   [-0.17466439]]
 
  C = [[1 0]]
 
  D = [[0]]
 
  dt = 0.6
  ]], dtype=object))

In [3]:
F = np.array([[np.cos(omega*h), np.sin(omega*h)], [-np.sin(omega*h), np.cos(omega*h)]])
G = 1.0/omega* np.array([[np.sin(omega*h)], [np.cos(omega*h)-1]])
sys_d2 = cm.ss(F,G,C,D, h)

In [35]:
Tc = np.linspace(0,30, 1000)
uc = Tc # Unit ramp
Td = h*np.arange(400)
ud = Td

(yr, tc, xc) = cm.lsim(sys_c, uc, Tc)
(yrd, td, xd) = cm.lsim(sys_d, ud, Td)

plt.plot(tc,yr)
plt.plot(td,yrd[0], '*')
plt.xlim([0,14])


Out[35]:
(0, 14)

In [36]:
lmbda = np.linalg.eigvals(F)
p = np.linalg.eigvals(A)

lmbda - np.exp(p*h)


Out[36]:
array([ 0. -1.38777878e-17j,  0. +1.38777878e-17j])

In [37]:
p


Out[37]:
array([ 0.+1.j,  0.-1.j])

In [38]:
lmbda


Out[38]:
array([ 0.99500417+0.09983342j,  0.99500417-0.09983342j])

In [39]:
import pdb
def mysim(F,G,C,D,u,x0):
    x = x0
    y = []
    for uk in u:
        yk = np.dot(C,x) + np.dot(D,uk)
        y.append(yk[0])
        x = np.dot(F,x) + np.dot(G, uk)
    
    return np.array(y)
        
u = np.ones(10)
x0 = np.zeros((2,1))
y = mysim(F,G,C,D,u,x0)
y-yd[0][:10]


Out[39]:
array([[  0.00000000e+00,  -9.98334166e-02,  -1.98669331e-01,
         -2.95520207e-01,  -3.89418342e-01,  -4.79425539e-01,
         -5.64642473e-01,  -6.44217687e-01,  -7.17356091e-01,
         -7.83326910e-01],
       [  9.98334166e-02,   1.38777878e-17,  -9.88359141e-02,
         -1.95686790e-01,  -2.89584926e-01,  -3.79592122e-01,
         -4.64809057e-01,  -5.44384271e-01,  -6.17522674e-01,
         -6.83493493e-01],
       [  1.98669331e-01,   9.88359141e-02,   2.77555756e-17,
         -9.68508759e-02,  -1.90749012e-01,  -2.80756208e-01,
         -3.65973143e-01,  -4.45548356e-01,  -5.18686760e-01,
         -5.84657579e-01],
       [  2.95520207e-01,   1.95686790e-01,   9.68508759e-02,
          0.00000000e+00,  -9.38981356e-02,  -1.83905332e-01,
         -2.69122267e-01,  -3.48697481e-01,  -4.21835884e-01,
         -4.87806703e-01],
       [  3.89418342e-01,   2.89584926e-01,   1.90749012e-01,
          9.38981356e-02,   5.55111512e-17,  -9.00071963e-02,
         -1.75224131e-01,  -2.54799345e-01,  -3.27937749e-01,
         -3.93908567e-01],
       [  4.79425539e-01,   3.79592122e-01,   2.80756208e-01,
          1.83905332e-01,   9.00071963e-02,   1.11022302e-16,
         -8.52169348e-02,  -1.64792149e-01,  -2.37930552e-01,
         -3.03901371e-01],
       [  5.64642473e-01,   4.64809057e-01,   3.65973143e-01,
          2.69122267e-01,   1.75224131e-01,   8.52169348e-02,
          1.11022302e-16,  -7.95752138e-02,  -1.52713618e-01,
         -2.18684436e-01],
       [  6.44217687e-01,   5.44384271e-01,   4.45548356e-01,
          3.48697481e-01,   2.54799345e-01,   1.64792149e-01,
          7.95752138e-02,   2.22044605e-16,  -7.31384037e-02,
         -1.39109222e-01],
       [  7.17356091e-01,   6.17522674e-01,   5.18686760e-01,
          4.21835884e-01,   3.27937749e-01,   2.37930552e-01,
          1.52713618e-01,   7.31384037e-02,   4.44089210e-16,
         -6.59708187e-02],
       [  7.83326910e-01,   6.83493493e-01,   5.84657579e-01,
          4.87806703e-01,   3.93908567e-01,   3.03901371e-01,
          2.18684436e-01,   1.39109222e-01,   6.59708187e-02,
          4.44089210e-16]])

In [ ]:


In [41]:
pdb.pm()


> /home/kjartan/.virtualenvs/jupyter_notebook/src/ipython/IPython/core/compilerop.py(100)ast_parse()
-> return compile(source, filename, symbol, self.flags | PyCF_ONLY_AST, 1)
(Pdb) 
(Pdb) 
(Pdb) 
(Pdb) q

In [42]:
np.dot(F,x0) + np.dot(D, u[0])


Out[42]:
array([[ 0.],
       [ 0.]])

In [46]:
plt.plot(range(len(y)),y, 'r*')
plt.plot(range(len(yd[0])), yd[0], 'bo')
plt.xlim([0,10])


Out[46]:
(0, 10)

In [1]:
7.5/23.5


Out[1]:
0.3191489361702128

In [ ]: