2019年12月26日 | イノマタの趣味部屋

「pythonによる制御工学入門」を読ませて頂きました。面白かったのでご紹介します。

サポートページ =>『Pythonによる制御工学入門』サポートページ

期間:12日間程度(1日2~3時間)

特徴

内容

pythonの基礎
制御モデル
- 伝達関数モデル
- 状態空間モデル
制御対称のふるまい
- 時間応答
- 周波数応答
制御系設計
- PID制御
- 2自由度制御
- ゲインチューニング(限界感度法)
- ゲインチューニング(モデルマッチング)
- 状態フィードバック制御(極配置法)
- 状態フィードバック制御(最適レギュレータ)
安定性
- 閉ループ
- 開ループ
アドバンスドな制御系設計
- オブザーバ
- ロバスト制御

感想

私の事前知識:

古典制御の伝達関数やPID制御は授業で出てきたので何となくイメージ出来ていました(わかるとは言ってない)。現代制御は全くわかりませんでした。

pythonに関してはゆるふわで常用しています。

全体の感想

表題からわかるようにpythonですぐ動かせます。パラメータを弄って理解を深めたり、自分で設計するうえで真似すればに使えるため便利です。

実用や、制御を勉強するための入門として全体を俯瞰する目的で書かれていると思われます。そのため、後半は特にあっさりした説明になっています。

pythonの記法に関しては、(私は読み飛ばしましたがパッと見た感じ)必要な部分は書いてありそうです。pythonを使ったことのない方には多少大変かもしれません。

モジュールはNumpy,Matplotli,Scipy,Sympy,Python-control(, slycotもインストールする必要あり)。 Python-control (リンクは2019/12/25時点で最新バージョンのマニュアルへ)はMATLABっぽく書けるパッケージなので、MATLABを使ったことがある方に良いかもしれません。

モデル

最適レギュレータを下で紹介していますので、状態空間モデルに関して少し書きます。

伝達関数モデルではフィードバック系全体について、出力y(t)、入力u(t)に対してP(t)を

$$\frac{y}{u}=P(t)$$

の形で表します。それに対して状態空間モデルでは、内部の状態を全て観測して、

$$\dot{ \textbf{x}}(t)= \textbf{Ax}(t)+ \textbf{Bu}(t)$$ $$\textbf{y}(t)= \textbf{Cx}(t)+ \textbf{Du}(t)$$

の形で表します。例えば、このようなアームを考えてみます。

トルクを$$\tau$$

粘性摩擦による力を$$\mu \dot{\theta}$$

重力加速度をgとすると、運動方程式は

$$\tau(t)=J\ddot{\theta}+\mu \dot{\theta}+Mglsin\theta$$

と書けます。入力u(t)、出力y(t)をそれぞれ

$$u(t)=\tau(t)$$

$$y(t)=\theta(t)$$

とし、角度 $\theta$ を十分小さいとして $sin(\theta)$ を

$$sin(y(t))=y(t)$$

と近似すると、

$$u(t)=J\ddot{y}(t)+\mu \dot{y}(t)+Mgly(t)$$

と書けます。これをもとに伝達関数が作れます。

さらに、状態空間モデルでは、

$$\textbf{x}(t)= \begin{bmatrix} \theta & \dot{\theta} \end{bmatrix}^T$$

と置いて

$$\dot{\textbf{x}}(t)= \begin{bmatrix} \dot{\theta} \cr \ddot{\theta} \end{bmatrix} =\begin{bmatrix} 0 & 1\cr -\frac{Mgl}{J} & -\frac{\mu}{J} \end{bmatrix}\textbf{x}(t) + \begin{bmatrix} 0 \cr \frac{1}{J} \end{bmatrix}u(t) $$

$$y(t)=\theta(t)=\begin{bmatrix} 1 & 0 \end{bmatrix} \textbf{x}(t) $$

となります。状態空間モデルの形にまとめると、

$$\dot{ \textbf{x}}(t)= \textbf{Ax}(t)+ \textbf{Bu}(t)$$ $$\textbf{y}(t)= \textbf{Cx}(t)+ \textbf{Du}(t)$$

$$\textbf{A}=\begin{bmatrix} 0 & 1\cr -\frac{Mgl}{J} & -\frac{\mu}{J} \end{bmatrix},\textbf{B}=\begin{bmatrix} 0 \cr \frac{1}{J} \end{bmatrix}, \textbf{C}=\begin{bmatrix} 1 & 0 \end{bmatrix},\textbf{D}=0$$

のように表せます。

ただし、状態空間モデルでは各パラメータが観測出来る必要があります。

序盤はかなりわかりやすく書かれています。古典制御と現代制御の基本的な部分はわかった気がします。

制御系設計

pythonで自分で動かせるという点で最もよかったのが設計の部分でした。教科書のグラフをサラッと見るよりわかった感がありました。

授業で使った教科書ではPID制御の理論をやったあとで、「なんで制御できるのかはわかったけど、kp,ki,kdの具体的な調べ方がわからん」と思ってたので、ゲインチューニングがコードと一緒に書いてあるのはかなり有り難かったです。

MATLABのようにマクローリン展開などが一瞬で出来ます。

(ここでちょっと思ったのが、MATLABのような有料ツールに近いことが無料で出来るのだから計算速度等が落ちているのではということ。実際のところはどうなのかは調べてません。)

個人的ハイライト:最適レギュレータ

p.176~状態空間モデルでのフィードバックゲインFを下のように与えます。

$$u=\textbf{Fx}$$

この F を"いい感じ"(雑)に選ぶ為に使われる手法の1つが最適レギュレータらしいです。

最適レギュレータでは、最適な値を求めるために、評価関数を最小にすることを考えます。

$$\textbf{Q}=\textbf{Q}^T>0,\textbf{R}=\textbf{R}^T>0$$

に対して評価関数

$$J=\int^\infty_0\textbf{x}(t)^T\textbf{Qx}(t)+u^T(t)\textbf{R}u(t)dt$$

を最小化するコントローラは

$$u=\textbf{F}_{opt}\textbf{x}$$

$$\textbf{F}_{opt}=\textbf{R}^{-1}\textbf{B}^T\textbf{P}$$

の形で得られます。ただし、

$$\textbf{P}=\textbf{P}^T>0$$

はリカッチ方程式

$$\textbf{A}^T\textbf{P}+\textbf{A}\textbf{P}+\textbf{P}\textbf{B}\textbf{R}^{-1}\textbf{B}\textbf{P}+\textbf{Q}=0$$

を満たす唯一の正定対象行列です。 ~~なんでこうなるのかは詳しくは書いてませんが今後調べます。~~ また、Jの最小値は

$$\textbf{x}(0)^T\textbf{Px}(0)$$

です。

python上でMATLABの関数に似たものを使用できるpython_controlというものがあります。これで、一見難しそうな計算が一瞬の内に出来ます！！これはやばい。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
from control.matlab import *
# グラフの設定
def plot_set(fig_ax,*args):
    # x軸のラベルを１つ目を引数に
    fig_ax.set_xlabel(args[0])
    # y軸のラベルを2つ目を引数に
    fig_ax.set_ylabel(args[1])
    fig_ax.grid(ls=':')
    if len(args)==3:
        fig_ax.legend(loc=args[2])

A = [[0,1],[-4,-5]]
B = [[0],[1]]
C= np.eye(2)
D = np.zeros([2,1])
P= ss(A,B,C,D)#状態空間モデル
#最適レギュレータによるFの計算
Q = [[100,0],[0,1]]
R=1
X,E,F=care(P.A, P.B, Q, R)#フィードバックゲインF,リカッチ方程式の解X,閉ループ極Eを返す。
F=-F
#状態空間モデルでの応答
Acl = P.A+P.B*F
Pfb = ss(Acl, P.B, P.C, P.D)#状態空間モデル
Td = np.arange(0, 5, 0.01)#時間
Ud = 1*(Td>0)#時間Tdが正の範囲では1,負では0となるステップ信号。
X0 = [-0.3,0.4]#初期値
xst,t=step(Pfb,Td)#ステップ応答
xin,_=initial(Pfb,Td,X0)
x,_,_=lsim(Pfb,Ud,Td,X0)

fig, ax = plt.subplots(1,2,figsize=(6,2.3))
for i in [0,1]:
    ax[i].plot(t,x[:,i])
    ax[i].plot(t,xst[:,i], ls='-')
    ax[i].plot(t,xin[:,i], ls='-.')

plot_set(ax[0],'t','$x_1$')
plot_set(ax[1],'t','$x_2$')
plt.savefig('最適レギュレータによるフィードバック制御.png')