図11.6, 11.7

In [1]:
import numpy as np
import sympy as sym
from control import matlab
from matplotlib import pyplot as plt
from scipy import arange
In [2]:
#システムパラメータを与える
A = np.array([[0.0, 1.0], [-6.0,-5.0]]) #行列A
b = np.array([[0.0], [1.0]]) #ベクトルb
c = np.array([[1.0, 0.0]]) #ベクトルc
d = 0 #スカラーd
In [3]:
#配置するオブザーバの極を与える
op = [-5+6j, -5-6j] #-4，-5をオブザーバの極とする

#オブザーバゲインを求める．双対システムA^Tとc^Tの状態フィードバックゲインを求め，転置を取るとオブザーバゲインとなる
h = matlab.place(A.T, c.T, op).T #双対システムで極配置し，転置を求める
print("observer gain =", h) #オブザーバゲインの表示

observer gain = [[ 5.]
 [30.]]
In [4]:
#Ah = A-hcとオブザーバの極の計算と表示
Ah = A - h*c #Ah = A-hcの計算
(eig_obs, eig_vec_obs) = np.linalg.eig(Ah) #Ahの固有値を求める
print('Ah = ',Ah) #Ahを表示
print('observer pole =',eig_obs) #オブザーバの極を表示

Ah =  [[ -5.   1.]
 [-36.  -5.]]
observer pole = [-5.+6.j -5.-6.j]
In [5]:
#システムとオブザーバ（誤差システム（11.15)式）の状態空間表現を求める
cc = np.array([[1.0, 0.0],[0.0,1.0]])#状態変数をプロットするためにcを単位行列にする
dd = np.array([[0.0],[0.0]])#cに合わせてdを零ベクトルとする

sys = matlab.ss(A, b, cc, dd) #システムの状態空間表現
e_sys = matlab.ss(Ah, b, cc, dd) #誤差システムの状態空間表現
In [6]:
#時間変数の定義
t = arange(0, 10, 0.01) #0から10まで0.01刻み

#システムとオブザーバの初期値を与える
x0 = np.array([[1.0],[1.0]]) #x(0) = [1;1]'
x0_ob = np.array([[0.0],[0.0]])  #xob(0) = [0;0]'

#初期値応答の計算
y, t = matlab.initial(sys, t, x0) #システムの初期値応答
e, t = matlab.initial(e_sys, t, x0_ob-x0) #誤差システムの初期値応答

#オブザーバの状態を求める
h_x = e + y #e(t) = hat{x}(t) - x(t)よりhat{x}を求める
In [7]:
# 図11.6(a)のプロット
plt.plot(t, y[:,0], label = "x_1(t)")#x_{1}を抽出してplotする
plt.plot(t, y[:,1], label = "x_2(t)")#x_{2}を抽出してplotする
plt.plot(t, h_x[:,0], label = "hat{x}_1(t)")#hat{x}_{1}を抽出してplotする
plt.plot(t, h_x[:,1], label = "hat{x}_2(t)")#hat{x}_{2}を抽出してplotする
plt.xlim([0,10]) #横軸（時間軸）の範囲の指定
plt.ylim([-1.0, 3.0]) #縦軸の範囲の設定
plt.grid(color='gray') #罫線を表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("state") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示
In [8]:
# 図11.6(b)の描画
plt.plot(t, e[:,0], label = "e_1(t)")#e_{1}を抽出してplotする
plt.plot(t, e[:,1], label = "e_2(t)")#e_{2}を抽出してplotする
plt.xlim([0,10]) #横軸（時間軸）の範囲の指定
plt.ylim([-2.0, 4.0]) #縦軸の範囲の設定
plt.grid(color='gray') #罫線を表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("state") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示
In [9]:
observer_poles2=[-4,-5] 
h2 = matlab.place(A.T, c.T, observer_poles2).T
Ah2 = A - h2*c
e_sys2 = matlab.ss(Ah2, b, cc, dd)

#時間変数の定義
t = arange(0, 10, 0.01) #0から10まで0.01刻み

#システムとオブザーバの初期値を与える
x0 = np.array([[1.0],[1.0]]) #x(0) = [1;1]'
x0_ob = np.array([[0.0],[0.0]])  #xob(0) = [0;0]'

(e2, t) = matlab.initial(e_sys2, t, x0_ob - x0) #誤差システムの初期値応答
#オブザーバの状態を求める
h_x2 = e2 + y #e(t) = hat{x}(t) - x(t)よりhat{x}を求める

# 図11.7(a)の描画
plt.plot(t, y[:,0], label = "x_1(t)")#x_{1}を抽出してplotする
plt.plot(t, y[:,1], label = "x_2(t)")#x_{2}を抽出してplotする
plt.plot(t, h_x2[:,0], label = "hat{x}_1(t)")#hat{x}_{1}を抽出してplotする
plt.plot(t, h_x2[:,1], label = "hat{x}_2(t)")#hat{x}_{2}を抽出してplotする
plt.xlim([0,4]) #横軸（時間軸）の範囲の指定
plt.ylim([-1.0, 3.0]) #縦軸の範囲の設定
plt.grid(color='gray') #罫線を表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("output") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示
In [10]:
# 図11.7(b)の描画
plt.plot(t, y[:,0], label = "x_1(t)")#x_{1}を抽出してplotする
plt.plot(t, y[:,1], label = "x_2(t)")#x_{2}を抽出してplotする
plt.plot(t, h_x2[:,0], label = "hat{x}_1(t)")#hat{x}_{1}を抽出してplotする
plt.plot(t, h_x2[:,1], label = "hat{x}_2(t)")#hat{x}_{2}を抽出してplotする
plt.xlim([0,4]) #横軸（時間軸）の範囲の指定
plt.ylim([-1.0, 3.0]) #縦軸の範囲の設定
plt.grid(color='gray') #罫線を表示
plt.xlabel("time   t[s]") #横軸のラベル表示
plt.ylabel("output") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示