無料で実践ロバスト制御③ ー摂動を持つ制御対象ー

平田氏の著書「実践ロバスト制御」↓を無料のPython-Controlで勉強する試み③。

実践ロバスト制御 (システム制御工学シリーズ)

作者:平田光男
コロナ社

Amazon

まずはロバスト制御系設計の流れを体験してみたい。

MATLABやPython-Controlを使う前提で
理論は後回し
簡単な例で

とすれば，だいぶ楽なはず。

例として，上記教科書4.2節の「2自由度振動系に対する設計例」が良さそうなのでここから始める。

摂動を持つ制御対象のモデル

下図(教科書の図4.3)のような制御対象を扱う。
f:id:manva:20210626224852p:plain
制御対象への入力 $u$ は，リニアモータのイメージで電圧[V]としており，推力定数 $K_s$ [N/V]をかけて力[N]にしている。( $K_s=1$ として，入力 $u$ は質量 $m_1$ に加えられる力と考えた方がわかりやすいかもしれない)。観測できる出力 $y$ は質量 $m_2$ の位置 $p_2$ である。
この制御対象の運動方程式を状態方程式の形で表す。状態変数を

$\displaystyle{ x=\begin{bmatrix}p_1\\p_2\\\dot{p_1}\\\dot{p_2} \end{bmatrix} }$

とし，

$\displaystyle{\begin{eqnarray} M&=&\begin{bmatrix}m_1 & 0 \\0 & m_2 \end{bmatrix}\\ C&=&\begin{bmatrix}c_1+c_2 & -c_2 \\-c_2 & c_2 \end{bmatrix}\\ K&=&\begin{bmatrix}k_1+k_2 & -k_2 \\-k_2 & k_2 \end{bmatrix}\\ F&=&\begin{bmatrix}K_s \\0 \end{bmatrix} \end{eqnarray}}$

とおくと，状態方程式を $\dot{x}=A_px+B_pu$ ，出力方程式を $y=C_px$ とすると，

$\displaystyle{\begin{eqnarray} A_p&=&\begin{bmatrix}O & I \\-M^{-1}K & -M^{-1}C \end{bmatrix}\\ B_p&=&\begin{bmatrix}O\\M^{-1}F\end{bmatrix}\\ C_p&=&\begin{bmatrix}0&1&0&0\end{bmatrix} \end{eqnarray}}$

$k_2$ に20%の摂動(不確かさ，ばらつき)があるとしている。

教科書の「プログラム4.1」を移植

プログラム4.1は，上記制御対象の入力 $u$ から出力 $y$ までのボード線図を描くものである。

変更箇所

・摂動の表現　最初の問題は，摂動を持つ制御対象をどう表すか。MATLABのRobust Control Toolboxでは，urealという便利な関数があるのだが，Python-Controlでは使えない。調べてみると，南氏の著書のソースコード↓があり，摂動を表現できていた。
https://github.com/373yuki/pyctrl/blob/master/chap7.ipynb
どうやっているのか確認。まず，

delta = np.arange(-1, 1 , 0.1)

として，-1から1まで0.1刻みの等差数列(arangeは等差数列)を用意し，for文で

P = (1 + WT*delta[i])*Pn

のようにずらして繰返し計算していた。これを真似してやってみる。

・行列の定義

M = np.matrix([ [ m1,  0 ],
                [  0, m2 ] ]);

のように，matrix関数を使わないと配列(array)になってしまう。matrix関数を使って定義しておけば「*」が行列の掛け算になる。

・inv(M) →　np.linalg.inv(M)
・zeros(2,2)　→　np.zeros( (2, 2) )
・eye(2,2)　→　np.eye(2)

・行列の結合　　MATLABのように

Ap = [ [ np.zeros((2,2)),  np.eye(2) ],
       [           -iM*K,      -iM*C ] ];

と書いてみたが，これじゃだめらしい。concatenate関数を使う。axis=1を指定すれば横にくっつける，デフォルトでは縦にくっつける。

Ap = np.concatenate([
     np.concatenate([np.zeros((2,2)),  np.eye(2)], axis=1),
     np.concatenate([          -iM*K,     -iM*C ], axis=1) ]);

・ボード線図のゲインの方のｙ軸範囲を変更する方法は↓を参考に，下記のようにsubplotの軸ハンドラを取得しておくことで解決。
https://jckantor.github.io/CBE30338/05.03-Creating-Bode-Plots.html

plt.subplot(2,1,1);
plt.subplot(2,1,2);
fig1, fig2 = plt.gcf().axes

・bode関数で、plot=Trueでプロットすると位相はdeg単位だが，mag,phase,om = bode( );のように戻値として受取ったphaseは，rad単位になるらしい。deg=Trueとしても同じ。なのでプロットするとき*180/np.piとした。
・目盛りを20dB, 90deg刻みにする

plt.yticks(np.arange(-60, 41, 20))
plt.yticks(np.arange(-360, 1, 90))

移植したコード「プログラム4.1」

# 実践ロバスト制御4.2節
import numpy as np
import matplotlib.pyplot as plt
from control.matlab import *

# %% defplant.m
# %% パラメータの定義
m1 = 0.8;
m2 = 0.2;
k1 = 100;

# k2の摂動によらない部分だけ先に計算
c1 = 1;
c2 = 0.3;
Ks = 100;
# %% 運動方程式のM,K,C行列を定義
M = np.matrix([ [ m1,  0 ],
                [  0, m2 ] ]);
C = np.matrix([ [ c1+c2, -c2 ],
                [ -c2,    c2 ] ]);
F = np.matrix([ [ Ks ],
                [  0 ] ]);
# %% 状態空間実現
iM = np.linalg.inv(M);
Bp = np.concatenate([ np.zeros((2,1)), iM*F ]);
Cp = [ 0, 1, 0, 0 ];
Dp = 0;

plt.figure(1);
plt.subplot(2,1,1);
plt.ylim(-60, 40);
plt.yticks(np.arange(-60, 41, 20))
plt.xlim(1e0, 1e2);
plt.subplot(2,1,2);
plt.ylim(-360, 0);
plt.yticks(np.arange(-360, 1, 90))
plt.xlim(1e0, 1e2);
fig1, fig2 = plt.gcf().axes # subplotのハンドラを取得

# k2 = ureal('k2',300,'percent',20);
delta = np.arange(-1, 1 , 0.1)
for i in range(len(delta)):
    k2 = 300+300*0.2*delta[i];
    K = np.matrix([ [ k1+k2, -k2 ],
                    [ -k2,    k2 ] ]);
    Ap = np.concatenate([
         np.concatenate([np.zeros((2,2)),  np.eye(2)], axis=1),
         np.concatenate([          -iM*K,     -iM*C ], axis=1) ]);
    # %% 制御対象の定義
    P = ss(Ap,Bp,Cp,Dp);
    #bode(P,{1e0,1e2}) # % ボード線図
    mag,phase,om = bode(P, plot=False);
    plt.sca(fig1);
    plt.semilogx(om, 20*np.log10(mag));
    plt.sca(fig2)
    plt.semilogx(om, phase*180/np.pi);
plt.sca(fig1);
plt.xticks([]); # ゲインの方の横軸目盛りを消す。semilogxで目盛りを書いてしまうので最後に