In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%pylab inline
In [ ]:
import control as con
In [ ]:
K = 1
d = 0.5
T = 10
delay = 10
a0 = 1
a1 = (2 * d * T) #16
a2 = (T**2) #100
b0 = K
tf_1 = con.matlab.tf(K, [a2, a1, a0])
#print tf_1
ss_1a = con.matlab.tf2ss(tf_1)
#print ss_1a
d_num, d_den = con.pade(delay, 1)
tf_delay = con.tf(d_num, d_den)
ss_delay = con.series(tf_delay, tf_1)
#print con.matlab.tf2ss(ss_delay)
In [ ]:
d_yout, d_T = con.matlab.step(ss_delay)
yout, T = con.matlab.step(tf_1) # step without delay
plt.plot(d_T, d_yout, 'r-', label='poly_est')
plt.plot(np.add(d_T, delay), yout, 'g-', label='idealized') #delay in timeaxis!
In [ ]:
In [ ]:
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import control
from control.matlab import *
x = [[1],[2],[3]]
#x_dot = np.dot(A, x) + np.dot(b, .0) # x° = Ax + bu
#show(x_dot)
#y = np.dot(c, x)# + np.dot(d, .0) # y = cx + du
#show(y)
#uboot_ss = control.ss(A, b, c, d)
#control.pole(uboot_ss)
#control.bode_plot(control.feedback(uboot_ss))
A = [[0, 1, 0], [0, 0, 1], [0, 0, -0.005]]
B = [[1], [0], [0]]
C = [[1., 0, 0]]
D = [[0.]]
sys1 = control.ss(A, B, C, D)
#mag, phase, omega = bode(sys1)
arr1, arr2 = control.step(sys1)
control.pole(sys1)
#mag, phase, omega = control.bode(control.feedback(sys1))
#sys = ss("0. 1 0; 0. 0 1; 0 0 -0.005", "0; 0; 1", "1 0 0", "0")
#mag, phase, omega = bode(sys)
#real, imag, freq = control.nyquist_plot(control.feedback(sys1))
control.matlab.pzmap(control.feedback(sys1))
In [ ]:
import numpy as np
import sympy as sp
from sympy import init_printing
from sympy import pprint as show
from __future__ import division
def sat(v, u_max):
return np.clip(v, -u_max, u_max)
#Test
#print sat([2,-3,4,1], 2)
#print sat([[2, 3],[-3, -4],[4, 0],[1, 0.5]], 3)
In [ ]:
# uboot
A = np.matrix([[0., 1., 0.],
[0., 0., 1.],
[0., 0., -0.005]])
a = A[np.array(A).shape[1]-1] #last line of Regelungsnormalform (characteristical polynom coefficients)
b = np.matrix([[0], [0], [1.]])
d = 0
c = np.matrix([[1], [0], [0]])
u_max = 2.5e-5
#simulation time
T = 1.0/100.0
#Introduced for fun
u_max_sys = 2.5e-5
In [ ]:
from IPython.display import clear_output
# Some other regulator
max_iter = 150000
# Initial state
x0 = np.array([[0],[0],[-0.004]])
# Timeloop
x_t = x0
y_t = np.zeros(max_iter)
u_t = np.zeros(max_iter)
k_t = np.matrix([[ 1.67656250e-06],
[ 3.07437500e-04],
[ 3.10195000e-02]]
)
for t in range(1, max_iter):
## Calc u
#print k_t
#print x_t
u = np.dot(-k_t.T, x_t)
#print u
if t%10000 == 1:
clear_output()
print t*T, "seconds done ->", t/max_iter*100, "%"
#u = sat(u[0][0], u_max)
#print u
## System response
# System saturation u (trivial if u_max and u_max_sys are identical)
# Calc x_dot
#show(A)
#show(x_t)
x_dot = np.dot(A, x_t) + b * sat(u, u_max_sys)
#print x_dot
#print b*u
#print x_dot + b*u
x_t = x_t + x_dot*T
#print x_t
y_t[t] = x_t[0]
u_t[t] = u
#print y_t
import matplotlib.pyplot as plt
plt.plot(np.arange(0, len(y_t))*T, y_t, 'b')
plt.xlabel('t [s]')
plt.ylabel('y(t)')
In [ ]:
from IPython.display import Image
from IPython.core.display import HTML
Image(url = "https://www.eit.hs-karlsruhe.de/mesysto/fileadmin/images/Skript_SYS_V_10_0_2/Kapitel_10_3/Grafik_10_3_7_HQ.png")
In [ ]:
from numpy import linalg as LA
# Berechnung der inversen Steuerbarkeitsmatrix
n = A.shape[0]
Q = b #
for i in range(1, n):
Q = np.hstack([Q, LA.matrix_power(A,i)*b])
Q_inv = LA.inv(Q)
print Q_inv
In [ ]:
#Zeilenvektor t_1.T entspricht der letzten Zeile der inversen Steuerbarkeitsmatrix
Image(url="https://www.eit.hs-karlsruhe.de/mesysto/fileadmin/images/Skript_SYS_V_10_0_2/Kapitel_10_3/Formel_10_3_51_HQ.png")
In [ ]:
t1 = Q_inv[-1,:]
print t1
In [ ]:
# Berechnung der Transformationsmatrix
Image(url="https://www.eit.hs-karlsruhe.de/mesysto/fileadmin/images/Skript_SYS_V_10_0_2/Kapitel_10_3/Grafik_10_3_8_HQ.png")
In [ ]:
n = A.shape[0]
T = t1
for i in range(1, n):
T = np.vstack([T, t1*LA.matrix_power(A,i)])
print T
In [ ]:
#Bestimmung der Zustandsraumdarstellung in Regelungsnormalform
Image(url="https://www.eit.hs-karlsruhe.de/mesysto/fileadmin/images/Skript_SYS_V_10_0_2/Kapitel_10_3/Grafik_10_3_9_HQ.png")
In [ ]:
Image(url="https://www.eit.hs-karlsruhe.de/mesysto/fileadmin/images/Skript_SYS_V_10_0_2/Kapitel_10_3/Grafik_10_3_10_HQ.png")
In [ ]:
A0 = T*A*LA.inv(T)
b0 = T*b
c0 = (c.T * LA.inv(T)).T
print A0, "\n", b0, "\n", c0
In [ ]:
import optim_tools as tools
A = np.matrix([[-4, -2], [1, -1]])
b = np.matrix([[1], [0]])
c = np.matrix([[4],[0]])
d = 0
(A0, b0, c0, d0), Q = tools.get_Steuerungsnormalform(A, b, c, d)
print A0
print b0
print c0
print d0
In [ ]:
In [ ]:
In [ ]:
# uboot
A = np.matrix([[0., 1., 0.],
[0., 0., 1.],
[0., 0., -0.005]])
a = -A[-1][:].T #!!!!
b = np.matrix([[0], [0], [1.]])
d = 0
c = np.matrix([[1], [0], [0]])
sys = control.ss(A, b, c.T, d)
#mag, phase, omega = bode(sys1)
#arr1, arr2 = control.step(sys)
#plt.plot(arr2, arr1)
#plt.show()
plt.axis([-2, .01, -1.1, 1.1])
control.matlab.pzmap(sys)
roots = [-1, -1+1j, -1-1j]
k = control.place(A, b, roots)
print k
a_tilde = np.matrix(np.poly(roots)[:0:-1]).T
#print a_tilde
k1 = a_tilde - a
print "a~", a_tilde
print "a", a
print "k1", k1
control.matlab.pzmap(control.ss(A-b*k, b, c.T, d))
In [ ]:
xx = np.array([1,2,3,4,5,6])
xx[:0:-1]
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: