コンテンツへスキップ

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

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軸のラベルを1つ目を引数に
    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')
最適レギュレータのステップ応答。
( 青が応答、オレンジがゼロ状態応答、点線がゼロ入力応答。)

最適レギュレータの関数:

現代文明に感動しました。伝達関数モデルでのゲインチューニングの方法もすぐできて衝撃的でした。

次に読みたい

最後に参考文献としておすすめ本が載っているのもポイント高い。

短期集中:振動論と制御理論 工学系の数学入門 の感想です。

振動を例にモデル化しグラフを描くまで、制御の一番はじめの体験ができる本。

短期集中:振動論と制御理論 工学系の数学入門 日本評論社 吉田 勝俊

期間: ~2019年2月11日の3週間程度

内容(個人的まとめ)

  1. 力学の運動方程式
  2. 単振動,減衰振動,強制振動
  3. スケール変換
  4. ハーモニックバランス
  5. 複素数,ラプラス変換,伝達関数
  6. 解析力学入門
  7. 棒立て制御(倒立振子)
  8. 線形化

特徴

単振り子や倒立振子を取り上げて、ソフト「MATX」で計算を体験しながら理解できる本でした。全体的に簡単に書いているので、なんとなく理解できます。 キーワードが載っているので、次につなげることができます。

感想

(参考)読む前: 電気科3年後期のテスト期間でした。
力学は入試+学部1年でやったきりで、さっぱりわかりませんでした。 制御工学の初歩的知識は持っており、なんとなく伝達関数を記述できました。 また、私はpythonで計算しましたが、pythonはこれで初めてさわりました(最下部のリンク)。

全体

この本の内容は

うちの卒研生が配属直後の1カ月以内に本書の全課題をこなし、

本書 「はじめに」より

と書かれている通りの難易度です。

力学、振動の運動方程式の解き方など、非常に基本的な部分からスタートして倒立振子の姿勢制御までを紹介しています。 数値計算ソフトMATXを使用しながら演習できます。

私は MATX ではなくPythonでやりました。pythonは方程式を解いたりグラフを描いたりすることが簡単にできるのでおすすめです。MATXのコードが載っているので、コーディングの参考にできました。

力学がさっぱりわからない状態で、倒立振子の制御がやりたかった私に最適な本でした。参考にした本、この後に読むといい本が載っているのもかなり有難かったです。

1. 力学の運動方程式

よくある内容です。

2. 単振動,減衰振動,強制振動

電気科の私にはLRC回路が馴染みがありましたが、ばねとダンパーで復習できました。

3. スケール変換

2次方程式を解くのは面倒なので、パラメータを減らして1次方程式にします。このような操作をスケール変換と呼びます。 1次方程式にすると、プログラムで解くのが簡単になります。

4. ハーモニックバランス法

解き方です。

5. 複素数,ラプラス変換,伝達関数

基本的な内容です。

6. 解析力学入門

この本で、解析力学というものを初めて知りました。それだけでもこの本を読んで良かったと思いました。

解析力学では、運動の座標を変化させたり一般化して表現したりするみたいです。ラグランジュの運動方程式の知識があると、一般化された座標の運動方程式が比較的簡単に立てられます。

ラグランジュの運動方程式をこの本だけで理解するのは厳しかったです。

7. 棒立て制御(倒立振子)

ラグランジュの運動方程式は使わずに倒立振子の運動方程式を立てて、グラフに書きました(下記の関連記事)。途中式は(よくあることですが)書かれていないので大変でした。それっぽい波形が得られたので満足です。

8. 線形化

あまりよくわからなかったのでさらっと読みました。

次に読みたい

解析力学キャンパス・ゼミ 改訂2 マセマ出版社 馬場 敬之 (著)

力学 裳華房 原島 鮮 (著)

Pythonによる制御工学入門 オーム社 南 裕樹 (著)

関連記事

Jupyter Notebookで倒立振子の振動波形を表示する

Scilab/Xcosが必要になったので入れて使ってみます。動作の確認として、MATLABでやったのと同じ解析をしてみました。

Scilab/Xcosが必要になったので入れて使ってみます。windows10,Scilab6.0.2です。動作の確認として、MATLABでやったのと同じ解析をしてみました。

目次

  • 導入
    • scilabとは
    • Xcosとは
  • インストール
  • 使ってみる
    • ブロックを並べる
    • 設定
    • シミュレーション
    • 構文
      • polyとは
      • syslinとは
      • bodeとは
      • dampとは
  • MATLAB/SimlinkとScilab/Xocsの比較
    • MATLAB/Simlink
    • Scilab/Xcos

導入

Scilabとは

Scilabはフリーのオープンソフトウェアで、端的に言ってしまえばMATLABのフリー版です。様々な数値解析ができます。

詳しくは公式にて

MATLABに関してはこちらの記事をご覧ください。=>MATLAB使い方個人的まとめ

Xcos

動的な連続/離散時間システムのモデリングやシミュレーションができます。 MATLABにとってのSimlinkです。グラフィカルにブロックを並べることでモデリングでき、コンパイルやシミュレーションもまとめてできる優れものです。

X2C(https://www.lcm.at/virtuelle-entwicklung/x2c/)

インストール

公式ページ: Download Scilab 6.0.2 -scilab

言語設定でjapaneseを選択し、デフォルトのままfullバージョンでインストールします。

インストールページのスクショ

使ってみる

MATLABと同様に、コマンドウィンドウ、Xcosウィンドウ、SimNoteのウィンドウをそれぞれ使って解析を行います。アプリケーション/Xcosで起動するとXcosウィンドウとライブラリが現れます。

Scilabの画面
Xcos

ブロックを並べてブロック線図を作り、sin波を表示させてみます。

ブロックを並べる

まずは信号源、出力/表示などからsin波発生器、CSCOPE、CLOCKを見つけて右クリックかドラッグでXcos上に配置し、矢印をドラッグして配線します。 スコープは赤い矢印があるものを選びましょう。これはイベントシミュレーションを駆動するイベント信号になります。

ブロック図の作成

設定

それぞれのブロックで右クリックすると項目が現れて、その中からブロックパラメータを開きます。 CFSCOPEを設定します。[参考1]

  • Color or mark vector: スコープで表示するグラフ線の色またはマークの連続によるグラフを指定する。入力信号数にかかわらず8セットの数字が必要。数字が正のとき、連続線で色が指定でき、負のとき、マーク形状が指定できる。(plot2dコマンドに準ずる)
  • Output window number: グラフ・ウィンドウ番号の指定。
  • Output window position: グラフ・ウィンドウの表示位置指定。左上の位置をx,yのドットで示す。
  • Output window size: グラフ・ウィンドウの表示サイズをドット数で指定。
  • Ymin: グラフ・ウィンドウのY軸最小値
  • Ymax: グラフ・ウィンドウのY軸最大値
  • Refresh period: グラフ・ウィンドウのX軸(時間)最大値 シミュレーション設定時間がこの値より大きい場合、グラフをクリアし再度続きの時間から描画する。基本的使い方では、シミュレーション時間と一致させる。
  • Buffer size: 入力データを一時的に蓄えておくメモリ領域のサイズを指定。グラフ表示速度に大きく影響するので、変化を観察したい場合以外は、少なくともRefresh period以上に設定する。最小は2。
  • Accept herited events: (inherited?) 継承の無効(0)・有効(1)を指定。1にするとイベント・ポート赤が削除される。アクティベーションが上流で行われる場合、継承が可能。

SCOPEの積算時間を10に、最小値~最大値を-10~10にしました。 クロックとsin発生器のパラメータも設定します。初期化時間を0.0に、ゲイン(振幅)5にします。

シミュレーションの設定

シミュレーション

シミュレーションを開始します。シミュレーション/パラメータ設定から、[参考1]

  • Final integration time: シミュレーション終了時間[単位は任意]
  • Real time scaling: 正の実数nを入力すると、シミュレーション設定単位時間tがn×t秒としてシミュレーションが実行される。n=1とした場合、シミュレーションの設定単位時間が1秒として解析表示される。
  • Integrator absolute tolerance: ソルバーの絶対誤差公差
  • Integrator relative tolerance: ソルバーの相対誤差公差
  • Tolerance on time: ODE(常微分方程式)ソルバーの最小時間公差
  • Max integration time interval: ソルバー呼び出しの最大時間間隔。"too many calls"というエラーメッセージが出た場合、値を小さくする。(解析精度を落とす。)
  • Solver 0(CVODE) - 100(IDA): 数値ソルバーの選択。
  • maximum step size: ソルバーがとる最大時間ステップ

最終積算時間を1.0E1.0、すなわち10sにしました。 AEBならA*10^Bになります。

document/test_0620/image.scgはこのようになりました。指定した通りに動いていることが確認できます。

document/ test_0620/image.scg の出力

構文

ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 3. 付録: Scilab言語の基本構文 -SimCircuit

Scilab コンソールの使い方の基本

「;」を後ろにつけると結果を表示しないみたいですね。

polyとは

[参考2]によれば、ploy関数は以下のように書きます。

p=poly(実数もしくは行列,"変数")
-->a=[-3 -4];
-->p=poly(a,"x")
 p  = 

            2
   12 +7x +x 

コンソールに上のように入力すると、12+7x+x^2を返します。((x+3)(x+4)=0, のときx=-3,-4) pを多項式変数として扱う事ができます。

-->s=poly(0,'s')
s =

    s
polyをコンソールに入力したときの出力

syslinとは

[参考3]によれば、syslin関数は以下のように書きます。

q=syslin(連続時間システムc or離散時間システムd,多項式)
--> s=poly(0,'s');

--> P=syslin('c',1/(s+1))
 P  = 

           
     1     
   ------  
           
   1 + s   

bodeとは

[参考4]によれば、syslin関数は以下のように書きます。

bode(syslinで定義された線形システムの多項式,x軸の最小値,x軸の最大値,x軸の刻み,['コメント'])

コメントにタイトルなどをかきます。 こだわらないなら関数のみ与えればOKです。

dampとは

[参考5]によれば、syslin関数は以下のように書きます。

[wn,z] = damp(連続時間のsys関数)

MATLAB/SimlinkとScilab/Xcosの比較

某授業のレポート1でやった問題の伝達関数を使って確認していきます。

問題: オープンループの伝達関数L(s)が以下の時の、回ループと閉ループでのボーデ線図をそれぞれ表示します。

$$L(s)=\frac{0.1}{s(s+0.5)(s+1)}$$

MATLAB

  • オープンループ

$$L(s)=\frac{0.1}{s^3+1.5s^2+0.5s}$$

MATLABオープンループでのボーデ線図

位相余裕:59° 利得余裕:17.5dB

  • フィードバック $$L(s)=\frac{0.1}{s^3 + 1.5s^2 + 0.5s + 0.1}$$

A. 周波数応答

MATLABフィードバックでのボーデ線図

B. ステップ応答

設定がうまくいかず、ぎこちないですがこんな感じ。最終値は1です。

MATLABフィードバックでのステップ応答

Scilab

  • オープンループ

複数行のコマンドを何度も使用する場合、アプリケーション/SciNoteに書くとスクリプトのようにまとめて実行できます。 見やすさのために出力なしで実行しました。

SciNoteを使用している様子

A. 周波数応答

まずはコマンドで実行してみます。

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
bode(p)
show_margins(p)
[gm,fp]=g_margin(p) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(p) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 

出力

ゲイン余裕[dB]

   17.501225

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   59.289834

 ゲイン交点角周波数[Hz]

   0.0293667

ゲイン余裕:位相が-180°と交差する周波数(位相交差角周波数)においてゲインが0dBより何dB小さいか

位相余裕:ゲインが0dBと交差する周波数(ゲイン交差角周波数)

オープンループボーデ線図は

Scilab オープンループでのボーデ線図

MATLABと一致してますね。

B. ステップ応答

Scilab オープンループでのステップ応答
  • フィードバック

A. 周波数応答

コマンドからやります。フィードバックした場合の伝達関数G(s)は以下のようになるので

$$G(s)=\frac{L(s)}{1+L(s)}$$

コード

s=poly(0,'s');
p=syslin('c',0.1/(s*(s+0.5)*(s+1)))
g=p/(1+p)
disp(g)
bode(g)
show_margins(g)
[gm,fp]=g_margin(g) disp("ゲイン余裕[dB]") disp(gm) disp("位相交差角周波数[Hz]") disp(fp) 
[pm,fg]=p_margin(g) disp("位相余裕[°]") disp(pm) disp("ゲイン交差角周波数[Hz]") disp(fg) 
[wn,z] = damp(g); disp("固有周波数[Hz]") disp(wn) disp("減衰係数") disp(z) 

出力

                           
            0.1            
   ----------------------  
                    2   3  
   0.1 + 0.5s + 1.5s + s   

 ゲイン余裕[dB]

   16.258267

 位相交差角周波数[Hz]

   0.1125395

 位相余裕[°]

   114.71226

 ゲイン交差角周波数[Hz]

   0.0313482

 固有周波数[Hz]

   0.2964606
   0.2964606
   1.1378001

 減衰係数

   0.6108736
   0.6108736
   1.

同じ曲線を描いていることを確認しました。

Scilab フィードバックでのボーデ線図

B. ステップ応答

最終値1で、100秒まで計算しています。

Scilab フィードバックでのステップ応答

コードは[参考6]を参考にさせて頂きました。

雑記

  • Xcosのライブラリで固まったらEscキーを押しましょう!ドラッグより右クリックの方がいいです。
  • 周波数応答を調べるとき、関数に伝達関数を与えています。なのでXcosでブロック線図を描こうと思ったらブロックを並べて伝達関数を計算させて、それを取得して関数に入れる、という手順を踏むと思うので回りくどいです。ここでやった例はすぐ伝達関数を計算できるので、コマンドでやった方が早いと思いました。
  • 何気に変数ブラウザ便利ですね。

参考

1.ビギナーのための科学技術計算ソフト(Scilab)の使い方講座 2.3 Xcosによるシミュレーション・ツールとしての使い方 (GUI型ブロック線図入力モジュール) -SimCircuit

2. Scilab関数の説明:polyとは

3. Scilab関数の説明:syslinとは

4. Scilabでのボード線図の描き方

5. damp -scilab

6. scilabで制御系設計