Piotr Banaszkiewicz
Najpierw zaimportuję wszystkie potrzebne moduły Pythona do załadowania plików Matlaba oraz wygenerowania odpowiednich wykresów.
In [1]:
%matplotlib inline
import numpy as np
import scipy.io # do odczytania pliku z Matlaba
import scipy.signal # do obiektów transmitancji
import scipy.interpolate # do interpolacji i porównania wyników symulacji
import sympy # do obliczeń symbolicznych
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font="DejaVu Sans")
sns.set_context('poster')
s = sympy.symbols('s') # zdefiniowana tutaj, aby nie tracić czasu inicjalizacji za każdym razem gdy jest wywoływana RMSE3
In [2]:
matlab = scipy.io.loadmat("obiekt.mat")
In [3]:
obiekt = matlab['y']
matlab['y'].shape
Out[3]:
Jak widać, obiekt został zaimportowany poprawnie. Możemy teraz łatwo wyświetlić zapisaną odpowiedź skokową tego obiektu:
In [4]:
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
plt.legend()
Out[4]:
Odpowiedź tego rzeczywistego obiektu regulacji odbiega znacząco od "ładnych", akademickich obiektów regulacji. Wynika to z faktu, że na obiekt oddziaływały zakłócenia. Pomimo tego możemy powiedzieć o tym obiekcie, że albo jest niskiego rzędu wraz z opóźnieniem, albo wysokiego rzędu, gdyż tylko takie obiekty mają podobne odpowiedzi skokowe.
Transmitancja takiego obiektu wynosi: $$ G(s) = \frac{ke^{-s\theta}}{Ts+1} $$
Parametrami $k$, $T$ oraz $\theta$ możemy manipulować. Jednym z zastosowanych podczas zajęć podejść była identyfikacja metodą prób i błędów.
Własnoręcznie zidentyfikowałem takie parametry obiektu:
In [5]:
k = 2.135
T = 15.57
theta = 7.1
Zobaczmy odpowiedź skokową obiektu o ww. transmitancji i tychże parametrach:
In [6]:
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
obiekt1_pb = scipy.signal.lti([k], [T, 1])
time_space = np.linspace(0, 59, 1000)
obiekt1_pb_y = obiekt1_pb.output(U=[0 if t < theta else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt1_pb_y[1], label=u'Obiekt przybliżony metodą prób i błędów')
plt.legend()
Out[6]:
Niestety z powodu braku osbługi outputdelay w pakietach naukowych Pythona, musiałem opóźnienie zasymulować inaczej: odpowiedź skokowa (1) zaczyna się dopiero od $t = \theta$, wcześniej $t = 0$.
Możemy zobaczyć zbliżoną postać odpowiedzi do rzeczywistego obiektu. Aby sprawdzić jak "blisko" jest zasymulowany obiekt do obiektu rzeczywistego użyłem pomiaru pierwiastka błędu średniokwadratowego:
In [7]:
def RMSE_interp(y1, t1, y2, t2):
"1-D interpolate first series of data, compute Mean Square Error for both data"
interp = scipy.interpolate.InterpolatedUnivariateSpline(t1, y1)
y1_interpolated = interp(t2)
error = np.subtract(y2, y1_interpolated)
return np.sum(error ** 2) / len(error)
def RMSE1(args):
global obiekt
k, T, theta = args
obiekt2 = scipy.signal.lti([k], [T, 1])
time = np.linspace(0, 59, 1000)
obiekt2_y = obiekt2.output(U=[0 if t < theta else 1 for t in time], T=time)
return RMSE_interp(obiekt[:,0], range(60), obiekt2_y[1], time)
In [8]:
RMSE_interp(obiekt[:,0], range(60), obiekt1_pb_y[1], time_space)
Out[8]:
In [9]:
RMSE1([k, T, theta])
Out[9]:
Jak widać przy okazji zdefiniowałem funkcję, którą będzie później można wykorzystać w metodzie automatycznej optymalizacji, opartej na szukaniu minimum lokalnego. Funkcja ta, RMSE1, interpoluje odpowiedź obiektu rzeczywistego, gdyż jest on oparty tylko na 60 pomiarach. Aby uzyskać większą dokładność w znalezieniu $\theta$, postanowiłem wykorzystać 1000 punktów.
Jak wspomniałem wyżej, metoda ta opiera się na szukaniu minimum lokalnego w pobliżu pewnego punktu. W pakiecie scipy Pythona znajduje się wiele funkcji optymalizujących, ja skorzystam z najprostszej fmin.
In [10]:
obiekt1_args = k2, T2, theta2 = scipy.optimize.fmin(RMSE1, (k, T, theta))
Jak widać algorytmowi udało się znaleźć lepszą (mniejszą) wartość MSE:
In [11]:
RMSE1([k, T, theta]) - RMSE1([k2, T2, theta2])
Out[11]:
Parametry, które uzyskał algorytm, różnią się niewiele od tych odkrytych przeze mnie "ręcznie":
In [12]:
k2, T2, theta2
Out[12]:
In [13]:
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
obiekt1_pb = scipy.signal.lti([k], [T, 1])
time_space = np.linspace(0, 59, 1000)
obiekt1_pb_y = obiekt1_pb.output(U=[0 if t < theta else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt1_pb_y[1], label=u'Obiekt przybliżony metodą prób i błędów')
obiekt1_opt = scipy.signal.lti([k2], [T2, 1])
time_space = np.linspace(0, 59, 1000)
obiekt1_opt_y = obiekt1_opt.output(U=[0 if t < theta2 else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt1_opt_y[1], label=u'Obiekt przybliżony metodą optymalizacji')
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.legend(loc=2)
Out[13]:
Różnica jest tak niewielka, że musiałem powiększyć wyraźnie wykres, aby była chociaż trochę zauważalna. Świadczy to bardzo dobrze o ręcznym wyniku identyfikacji obiektu. Pozwolę sobie przypomnieć wyniki błędu średniokwadratowego. Metoda prób i błędów:
In [14]:
RMSE1((k, T, theta))
Out[14]:
Metoda optymalizacji:
In [15]:
RMSE1((k2, T2, theta2))
Out[15]:
Transmitancja takiego obiektu wynosi: $$ H(s) = \frac{ke^{-s\theta}}{(T_1s+1)(T_2s+1)} = \frac{ke^{-s\theta}}{T_1T_2s^2 + (T_1+T_2)s + 1}$$
Parametrami $k$, $T_1$, $T_2$ oraz $\theta$ możemy, jak w poprzednim przykładzie, manipulować. W przypadku tego obiektu będziemy posługiwać się tylko metodą automatycznego szukania minimalnego błędu średniokwadratowego.
Na początek zdefinuję funkcję RMSE2 zwracającą błąd średniokwadratowy dla nowego obiektu inercyjnego drugiego rzędu. Jest ona bardzo zbliżona do funkcji RMSE1 dla poprzedniego obiektu (nadal wykorzystujemy interpolację dla obiektu rzeczywistego):
In [16]:
def RMSE2(args):
global obiekt
k, T1, T2, theta = args
obiekt2 = scipy.signal.lti([k], [T1 * T2, T1 + T2, 1])
time = np.linspace(0, 59, 100)
obiekt2_y = obiekt2.output(U=[0 if t < theta else 1 for t in time], T=time)
return RMSE_interp(obiekt[:,0], range(60), obiekt2_y[1], time)
Zobaczmy, jak zadziała ta funkcja dla przykładowych danych $k = 2$, $T_1 = 2$, $T_2 = 13$, $\theta = 6$:
In [17]:
k = 2
T1 = 2
T2 = 13
theta = 6
print RMSE2((k, T1, T2, theta))
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
obiekt2 = scipy.signal.lti([k], [T1 * T2, T1 + T2, 1])
time_space = np.linspace(0, 59, 100)
obiekt2_y = obiekt2.output(U=[0 if t < theta else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt2_y[1], label=u'Obiekt przybliżony metodą optymalizacji')
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.legend(loc=2)
Out[17]:
Sprawdźmy, czy damy radę zoptymalizować ten wynik:
In [18]:
obiekt2_args = k_2, T1_2, T2_2, theta_2 = scipy.optimize.fmin(RMSE2, (k, T1, T2, theta))
Metoda optymalizacji otrzymała dużo lepszy wynik, chociaż trochę gorszy niż w przypadku obiektu inercyjnego 1. rzędu.
In [19]:
k_2, T1_2, T2_2, theta_2, RMSE2((k_2, T1_2, T2_2, theta_2))
Out[19]:
In [20]:
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
obiekt2 = scipy.signal.lti([k_2], [T1_2 * T2_2, T1_2 + T2_2, 1])
time_space = np.linspace(0, 59, 100)
obiekt2_y = obiekt2.output(U=[0 if t < theta_2 else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt2_y[1], label=u'Obiekt przybliżony metodą optymalizacji')
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.legend(loc=2)
Out[20]:
Podczas moich testów okazało się, że funkcja scipy.optimize.fmin jest bardzo czuła na punkt początkowy w przypadku tego obiektu i łatwo może znaleźć nieoptymalny wynik.
Transmitancja takiego obiektu wynosi: $$ J(s) = \frac{k}{(Ts+1)^n} $$
Parametrami $k$, $T$, oraz $n$ możemy, jak w poprzednich przykładach, manipulować. W przypadku tego obiektu będziemy ponownie posługiwać się tylko metodą automatycznego szukania minimalnego błędu średniokwadratowego.
W tym przykładzie użyjemy pakietu do obliczeń symbolicznych Pythona, Sympy. Wykorzystamy go, aby utworzyć wielomian $n$-tego stopnia, który jest mianownikiem szukanej transmitancji.
In [21]:
def RMSE3(args):
global obiekt
global s
k, T, n = args
n = int(n) # musi być całkowity wykładnik
poly = sympy.expand((T * s + 1) ** n).as_poly()
coeffs = poly.all_coeffs()
#print poly, coeffs, k
obiekt2 = scipy.signal.lti([k], [float(x) for x in coeffs])
time = np.linspace(0, 59, 100)
obiekt2_y = obiekt2.output(U=[1 for t in time], T=time)
return RMSE_interp(obiekt[:,0], range(60), obiekt2_y[1], time)
In [22]:
obiekt3_args = k, T, n = scipy.optimize.fmin(RMSE3, (2, 4, 3))
Już teraz możemy ocenić, że dla tych danych wejściowych ten obiekt jest najgorzej przybliżony. Zobaczmy to na wykresie:
In [23]:
n = int(n)
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
poly = sympy.expand((T * s + 1) ** n).as_poly()
coeffs = poly.all_coeffs()
obiekt3 = scipy.signal.lti([k], [float(x) for x in coeffs])
time_space = np.linspace(0, 59, 100)
obiekt3_y = obiekt3.output(U=[1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt3_y[1], label=u'Obiekt przybliżony metodą optymalizacji')
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.legend(loc=2)
Out[23]:
Przeprowadziłem analizę dla różnych typów obiektów. W przypadku każdego z nich identyfikacja przebiegła pomyślnie: obiekt rzeczywisty został dosyć dokładnie przybliżony numerycznie:
In [24]:
plt.plot(range(60), obiekt[:,0], label=u'Obiekt rzeczywisty')
k, T, theta = obiekt1_args
obiekt1 = scipy.signal.lti([k], [T, 1])
time_space = np.linspace(0, 59, 1000)
obiekt1_y = obiekt1.output(U=[0 if t < theta else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt1_y[1], label=u'Obiekt inercyjny 1. rzędu')
error1 = RMSE1(obiekt1_args)
k, T1, T2, theta = obiekt2_args
obiekt2 = scipy.signal.lti([k], [T1 * T2, T1 + T2, 1])
time_space = np.linspace(0, 59, 100)
obiekt2_y = obiekt2.output(U=[0 if t < theta else 1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt2_y[1], label=u'Obiekt inercyjny 2. rzędu')
error2 = RMSE2(obiekt2_args)
k, T, n = obiekt3_args
error3 = RMSE3(obiekt3_args)
n = int(n)
poly = sympy.expand((T * s + 1) ** n).as_poly()
coeffs = poly.all_coeffs()
obiekt3 = scipy.signal.lti([k], [float(x) for x in coeffs])
obiekt3_y = obiekt3.output(U=[1 for t in time_space], T=time_space)
plt.plot(time_space, obiekt3_y[1], label=u'Obiekt wieloinercyjny bez opóźnienia')
plt.xlabel(u'Czas')
plt.ylabel(u'Odpowiedź')
plt.legend(loc=4)
Out[24]:
Poniżej przedstawiam wyliczone błędy średniokwadratowe dla poszczególnych obiektów:
In [25]:
print "Obiekt inercyjny 1. rzędu: ", error1
print "Obiekt inercyjny 2. rzędu: ", error2
print "Obiekt wieloinercyjny bez opóźnienia: ", error3
Najlepiej dopasowany został obiekt inercyjny 2. rzędu z opóźnieniem.