機械系エンジニアの備忘録

20代独身社会人。仕事では機械・機構の研究開発を行っているエンジニアが、自分の専門分野ではないpythonを扱って楽しむブログです。

MENU

【python】【制御工学】pythonで状態フィードバック制御

python-controlで、状態フィードバック制御を行う

■はじめに

前回、python-controlで伝達関数の設定とインパルス・ステップ・初期値応答を計算しました。

stjun.hatenablog.com

前回の伝達関数は極の実軸が負なので元々安定なシステムでした。

ただし世の中のシステムはそのままでは不安定なシステムがあります。 

そこで今回は不安定なシステムに状態フィードバックを使用して安定化させたいと思います。

■コード

from control.matlab import *
import matplotlib.pyplot as plt
import numpy as np

#
A=np.matrix([[0,5],[-2,3]]);
b=np.matrix([[0],[1]]);
C=np.matrix([1,0]);
D=0;
sys=ss(A,b,C,D)

(y0,t0,x0)=initial(sys,X0=[2,3],T=np.arange(0,15,0.01),return_x=True)

print('rank is',np.linalg.matrix_rank(ctrb(A,b)))

print('pole is ',pole(sys))

p=[-3,-4]
k=place(A,b,p)
print('feedback gain is ',k)

sys_1=ss(A-b*k,b,C,D)
(y1,t1,x1)=initial(sys_1,X0=[2,3],T=np.arange(0,15,0.01),return_x=True)

print('pole is ',pole(sys_1))

#グラフの設定
fig=plt.figure(figsize=(10,8))
ax1=plt.subplot2grid((2,1),(0,0))
ax1.plot(t0,x0)
ax2=plt.subplot2grid((2,1),(1,0))
ax2.plot(t1,x1)
plt.show()

実行すると以下の画面がでてきます。上のグラフはフィードバック前のシステムに初期値[2,3]を与えた時の挙動を、下のグラフがフィードバック後です。

f:id:stjun:20191009230150p:plain

上のグラフはy軸が大きすぎるので分かりづらいですが、y軸を-50~50にすると以下のように不安定な挙動を示しているのがよくわかります。

f:id:stjun:20191009230553p:plain
それに対しフィードバック後の下のグラフは3秒あたりで0付近で安定しており、

フィードバック制御によって不安定なシステムを安定化できたことがわかります。

 

なお念のためランクの確認やフィードバック前後の極も出力してみました。

ランクは2なので制御可能であり、またフィードバック前は極が正なのに対し、フィードバック後は極が負になっており、フィードバックによりシステムの極を負に動かして安定化させたことが分かります。

f:id:stjun:20191009230220p:plain

 

■買って良かったもの紹介

今回はipadです。動画視聴や動画編集(Lumafusion)を使用しています。

12.9インチと迷ったのですが、ノートとして使う&持ち運ぶの2点から11インチにしました。満足しています。

ipad pro 12.9インチを買うならsurface pro6の方が良いと思います。windows10なのでpythonもフルに使えます。

 

 

■説明

前回は伝達関数を設定しましたが、下のように状態空間モデルを記述することができます。

A=np.matrix([[0,5],[-2,3]]);
b=np.matrix([[0],[1]]);
C=np.matrix([1,0]);
D=0;
sys=ss(A,b,C,D) 

またフィードバックをしているのは以下の部分になります。

p=[-3,-4]
k=place(A,b,p)
sys_1=ss(A-b*k,b,C,D)

1行目のpはどこに極を動かしたいかを示しています。今は極を[-3,-4]に動かそうとしています。

2行目で極をp=[-3,-4]に動かすためのフィードバック係数kを計算しています。

3行目は2行目で計算したフィードバック係数k込みでの状態空間モデルを記述しています。カッコ内の先頭がAからA-b*kに変わるだけです。

 

また極を調べたい場合は、pole(sys)

ランクを調べたい場合は、np.linalg.matrix_rank(ctrb(A,b))

その他にもdamp(sys)などにすると、固有振動数や減衰比なども計算できます。

 

■最後に

次回はPID制御を気が向いたらやってみます。