manvaのエンジニアリング魂

エンジニアリング・ものづくり・DIYをもっと身近にするためのブログ。インスピレーションを刺激します。

不要な条件なのに「必要条件」という違和感

数学

A\Rightarrow B（AならばB）が成り立つとき、BはAの必要条件であるという。AであるためにはBである「必要」があるから必要条件。うん、BじゃなかったらAじゃないからBである必要があるよな。わかる。
しかし、高校で習ったときからずっと自分の中で何か違う気がしていた。最近ようやくこの違和感の正体がわかった。以下の例を見ればスッキリ理解できる。

\displaystyle{\begin{eqnarray} A:\;x=1 \quad\iff\quad B:\;\left\{\begin{array}{l} x+1=2 \\ x>0 \end{array}\right. \end{eqnarray}}
「AならばB」が明らかに成り立っているので，BはAの必要条件である。（ついでに「BならばA」も成り立っているので互いに必要十分条件。）しかし，お気づきだろうが，x>0という条件はなくても A\Rightarrow B も B\Rightarrow A も成立するので，x>0は「不要」な条件である。不要な条件があるのに「必要条件」と呼ばれることに違和感を感じるのである。
「必要」というのは，AであるためにBが「成り立っていること」が必要なのであり，Bが複数の条件からなる場合，一つ一つの条件自体が必要という意味ではない，ということである。特に「必要十分条件」と言ったときは，言葉尻から，条件自体が過不足ないという意味のような気がしてしまうので違和感を感じたのだろう。
もしも受験生がこの文を読んだなら，数学の試験には一切役に立たないどころか余計混乱するかも知れないので，受験が終わるまでは忘れておくことをお勧めする。

無料で実践ロバスト制御⑥ 重み関数の決め方

制御工学 Python-Control

ほぼ教科書の受け売りになってしまうが，重み関数の決め方や意味を必要最小限だけ説明しておく。

重み関数W_T

無料で実践ロバスト制御④ ーH∞制御問題を解くー - manvaのエンジニアリング魂
の手順３で，重み関数W_Tは乗法的摂動を覆うように決めた。

これは何をしているか。
まず，乗法的摂動\Delta_mを覆うように線を引き，その関数をW_mと置いたとする。そして，乗法的摂動はこれに大きさが1以下の伝達関数\Delta (\|\Delta\|_\infty\leq1)をかけたものである，と考える。

\displaystyle{\begin{eqnarray} \Delta_m=\Delta W_m \end{eqnarray}}
ここで，下図のように摂動を持つ制御対象\tilde{P}を制御器Kで制御したときのロバスト安定性を考える。

図のaからbまでの伝達関数が-T(Tは相補感度関数)となることから，開ループ伝達関数は-\Delta W_m Tであり，\Deltaは1より小さいことと，ノルムを考える際，符号は無視して良いことから，
\displaystyle{\begin{eqnarray} \|W_mT\|_\infty<1 \end{eqnarray}}
を満たせばスモールゲイン定理より安定である。そのため，W_T=W_mとして，乗法的摂動を覆うように決めたW_mをそのまま相補感度関数に対する重み関数W_Tとして使って，
\displaystyle{\begin{eqnarray} \|W_TT\|_\infty<1 \end{eqnarray}}
となることを目指せば，ロバスト安定な制御器が作れることになる。
このように，W_Tは乗法的摂動\Delta_mから半自動的に決められる。

重み関数W_S

重み関数W_Sは感度関数Sの周波数特性を整形するための関数である。感度関数Sはフィードバック制御の性能を表す指標であり，小さいほどよいが，S+T=1の関係があるのでむやみに小さくすることはできない。そこで，低周波が小さくなるような関数を決め，感度関数Sがそれより小さくなるような制御器を設計することを目指す。重み関数W_S(s)はその関数の逆数として定義しておく。つまり，全ての\omegaで

\displaystyle{\begin{eqnarray} \left|S(j\omega)\right|<\frac{1}{|W_S(j\omega)|} \end{eqnarray}}
となることを目指す。W_S(s)を逆数の形で与えておくと，この条件は
\displaystyle{\begin{eqnarray} &\Leftrightarrow& \left|W_S(j\omega)S(j\omega) \right|<1, \quad \forall \omega\\ &\Leftrightarrow& \|W_SS\|_\infty<1 \end{eqnarray}}
と，H∞ノルムで表せる。
教科書では，W_Sは
\displaystyle{\begin{eqnarray} W_S=\frac{15}{s+0.015} \end{eqnarray}}
とほぼ積分にしていた。この0.015を0にして完全に積分にしてしまうと解けないらしく，エラーとなった。

混合感度問題

ここまでで，\|W_SS\|_\infty<1と\|W_TT\|_\infty<1を同時に満たす制御器を求める，という問題に帰着できたが，不等式が2つあると問題が解けないので，この十分条件として

\displaystyle{\begin{eqnarray} \left\| \begin{bmatrix}W_SS \\W_TT \\\end{bmatrix} \right\|_\infty<1 \end{eqnarray}}
を満たす制御器を求める。
「修正混合感度問題」は，
\displaystyle{\begin{eqnarray} W_S=W_{PS}P \end{eqnarray}}
として感度関数の重みに制御対象の極を持たせたものに相当する。

重み\epsilon

これでもまだ解けないので修正が必要。wからyへの直達項を追加する。
無料で実践ロバスト制御④ ーH∞制御問題を解くー - manvaのエンジニアリング魂
の話は，「入力端混合感度問題」と呼ばれるもので，入力端混合感度問題は，仮定A2と仮定A3(教科書参照)を満たさないので標準H∞制御問題となっていない。仮定A2を満たすため小さい重み\epsilonをかけて直達項を追加したので解けた(仮定A3は満たさなくても解けるらしい)。
無料で実践ロバスト制御⑤ ー修正混合感度問題ー - manvaのエンジニアリング魂
の「修正混合感度問題」でも同様に直達項が必要なようだ。

重み関数の決め方の練習

重み関数を自分で決めてみる。
無料で実践ロバスト制御④ ーH∞制御問題を解くー - manvaのエンジニアリング魂
で移植したプログラム4.2を修正し，相補感度関数の重み関数W_Tをもう少しタイトにしてみる。

#%% 重み Wm の定義
Wm = 3*s**2/(s**2+18*s+45**2);
mag, phase, om = bode( Wm, plot=False );
plt.semilogx( om, 20*np.log10(mag) );
Wm = 0.7*s**2/(s**2+5*s+40**2); #もう少しタイトにしてみる
mag, phase, om = bode( Wm, plot=False );
plt.semilogx( om, 20*np.log10(mag) );


このW_Tを使って，↓
無料で実践ロバスト制御⑤ ー修正混合感度問題ー - manvaのエンジニアリング魂
の修正混合感度問題を解いてみると，γが0.546となり，下図のように，感度関数，相補感度関数とそれぞれの重み関数(の逆数)の間にだいぶ隙間が空いた。

γが1を超えない範囲で感度関数をもっと小さくすれば性能が上がるはずである。

#%% 重み関数の定義
#Ws  = 15/(s + 1.5e-2);
#Wps = Ws*0.8;
Wps  = 44/(s + 1.5e-2);  #Wtをタイトにしたら余裕ができたのでこちらを詰める

とすると，γが0.989となり，外乱応答のピークは教科書の図4.17より小さくできた。

感度関数の重みをもっと大きくしてみる。

Wps  = 100/(s + 1.5e-2);

γが1.60となり，1を超えている。これは感度関数，相補感度関数のどこかで重み関数で抑えられない部分ができていることを表している。確認してみると↓，確かに感度関数，相補感度関数共に抑えられておらず，逆転している部分がある。

相補感度関数も抑えられていないため，ロバスト安定条件\|W_TT\|_\infty<1を満たしていないはずである。ところが，外乱応答を見ると↓，持続振動の減衰が悪くなってはいるものの，発散はしていない。

なぜか。スモールゲイン定理は，教科書には「必要十分条件」と書かれているが，「十分条件」だったとすれば理解できる。スモールゲイン定理とは，ボード線図で言えば開ループのゲインが0dBを超えるところがなければ安定，ナイキスト線図で言えば単位円から出るところがなければ安定，という当たり前の話と理解していた。だとすると，かなり保守的な十分条件であるはずである。では教科書が間違っているのかと思って調べてみると，そう単純でもなさそうで，Zhou氏らの教科書(Essentials of Robust ControlやRobust and Optimal Control。pdfがネットに上がっている。著作権とか怪しいのでリンクは貼らないでおく)でも必要十分条件として証明が書かれている。
このあたりの話は、後日調べて別記事に書いたので参考にして欲しい。
manva.hatenablog.com
manva.hatenablog.com

発散するまでゲインを上げてみる。

Wps  = 1000/(s + 1.5e-2);

γは2.60。まだ発散しない。外乱応答↓はかなり改善しちゃってる。

Wps  = 8000/(s + 1.5e-2);

くらいでやっと発散してくれた。もう何やってるかわからん。理論の意味ないな。

無料で実践ロバスト制御⑤ ー修正混合感度問題ー

制御工学 Python-Control

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

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

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

Amazon

前回までで，H∞制御問題の解き方の流れが把握できた。
まずその結果を確認する。

移植したコード「プログラム4.4」
#%% chkperf.m
#%% 閉ループ特性の確認
L = P_nominal*K;
T = feedback(L,1); #% T = L/(1+L)
S = feedback(1,L); #% S = 1/(1+L)
M = feedback(P_nominal,K); #% M = P/(1+L)

plt.figure(5);
#bodemag(T.nominal,'-',1/Wt,':',S.nominal,'--',1/Ws,'-.',w);
#legend('T.nomial','1/Wt','S.nominal','1/Ws',2)
w    = logspace(0,3,100); #% 周波数ベクトルの定義
mag1, phase1, om1 = bode( T, w, plot=False );
plt.semilogx( om1, 20*np.log10(mag1), '-', label='$T$' );
mag2, phase2, om2 = bode( 1/Wt, w, plot=False );
plt.semilogx( om2, 20*np.log10(mag2), ':', label='$1/W_T$' );
mag3, phase3, om3 = bode( S, w, plot=False );
plt.semilogx( om3, 20*np.log10(mag3), '--', label='$S$' );
mag4, phase4, om4 = bode( 1/Ws, w, plot=False );
plt.semilogx( om4, 20*np.log10(mag4), '-.', label='$1/W_S$' );
plt.xlim(1e0, 1e2);
plt.ylim(-30, 30);
plt.legend();

plt.figure(6);
yout, time = step( T, 2 );
plt.plot(time, yout);

plt.figure(7)
plt.subplot(2,1,1);
yout, time = impulse( M, 2 );
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);
    L = P*K;
    T = feedback(L,1); #% T = L/(1+L)
    M = feedback(P,K); #% M = P/(1+L)
    yout, time = step( T, 2 );
    plt.figure(6);
    plt.plot(time, yout);
    yout, time = impulse( M, 2 );
    plt.figure(7);
    plt.plot(time, yout);
plt.subplot(2,1,2);
yout, time = impulse( P_nominal, 2 );
plt.plot(time, yout);
実行結果

f:id:manva:20210711234836p:plain
教科書の図4.9。感度関数Sが重み関数W_Sの逆数で，相補感度関数Tが重み関数W_Tの逆数でそれぞれ抑えられているのがわかる。
f:id:manva:20210711234845p:plain
教科書の図4.10。指令に対するステップ応答の図。良好である。
f:id:manva:20210711234856p:plain
教科書の図4.11。外乱に対するインパルス応答の図。上はH∞制御した結果，下は制御対象(ノミナル)の元々の応答。あまり改善していない。混合感度問題では，極零相殺が起こってしまい，指令応答は良くても外乱応答は良くない，という問題がある。

修正混合感度問題

この問題を改善した，「修正混合感度問題」を解く。
前回と同じく，Python-Controlではsysicが使えないので一般化プラントをMIMOのtfで与える。図4.12(b)より，

\displaystyle{\begin{eqnarray} z_1&=&W_{PS}\left\{P(w_1-u)+\epsilon w_2\right\}\\ z_2&=&W_Tu\\ y&=&P(w_1-u)+\epsilon w_2 \end{eqnarray}}
なので，

WpsP = Wps*Pn;
WpsE = Wps*Weps;
num = [[WpsP.num[0][0], WpsE.num[0][0], (-WpsP).num[0][0]],
       [           [0],            [0],      Wt.num[0][0]], 
       [  Pn.num[0][0],         [Weps],   (-Pn).num[0][0]]];
den = [[WpsP.den[0][0], WpsE.den[0][0], (-WpsP).den[0][0]],
       [           [1],            [1],      Wt.den[0][0]],
       [  Pn.den[0][0],            [1],   (-Pn).den[0][0]]];
G_ = tf(num, den)
G = tf2ss(G_)

とした。
ここから教科書4.3節になるので，プログラムは今までのものを使わずにこれだけで動くように書き直した。

移植したコード「プログラム4.5」
# 実践ロバスト制御4.3節
#プログラム4.5
import numpy as np
import matplotlib.pyplot as plt
import control
from control.matlab import *

# %% defplant.m
# %% パラメータの定義
m1 = 0.8;
m2 = 0.2;
k1 = 100;
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;

#%% defpert.m
#%% 乗法的摂動の見積もり
w    = logspace(0,3,100); #% 周波数ベクトルの定義
k2 = 300;
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_nominal = ss( Ap, Bp, Cp, Dp );
Pn = ss2tf(P_nominal);
#%% 重み Wm の定義
s = tf('s');
Wm = 3*s**2/(s**2+18*s+45**2);

#%% defgp2.m
#%% 重み関数の定義
Ws  = 15/(s + 1.5e-2);
Wps = Ws*0.8;
Wt  = Wm;
Weps = 5e-4;

plt.figure(3)
#bodemag(Ws,':',Wps,Wps*P.nominal,'-.',Wt,'--',w);
#legend('Ws','Wps','Wps*P','Wt',4)
mag,phase,om = bode(Ws, w, plot=False);
plt.semilogx(om, 20*np.log10(mag), ':', label='$W_S$');
mag,phase,om = bode(Wps, w, plot=False);
plt.semilogx(om, 20*np.log10(mag), label='$W_{PS}$');
mag,phase,om = bode(Wps*P_nominal, w, plot=False);
plt.semilogx(om, 20*np.log10(mag), '-.', label='$W_{PS}*P$');
mag,phase,om = bode(Wt, w, plot=False);
plt.semilogx(om, 20*np.log10(mag), '--', label='$W_T$');
plt.xlim()
plt.ylim(-60, 30);
plt.legend();

#%% 一般化プラントの定義
#Pn = P.nominal;
#systemnames   = 'Pn Wps Wt Weps';
#inputvar      = '[w1; w2; u]';
#outputvar     = '[Wps; Wt; Pn+Weps]';
#input_to_Pn   = '[w1 - u]';
#input_to_Wps  = '[Pn + Weps]';
#input_to_Wt   = '[ u ]';
#input_to_Weps = '[ w2 ]';
#G = sysic;
WpsP = Wps*Pn;
WpsE = Wps*Weps;
num = [[WpsP.num[0][0], WpsE.num[0][0], (-WpsP).num[0][0]],
       [           [0],            [0],      Wt.num[0][0]], 
       [  Pn.num[0][0],         [Weps],   (-Pn).num[0][0]]];
den = [[WpsP.den[0][0], WpsE.den[0][0], (-WpsP).den[0][0]],
       [           [1],            [1],      Wt.den[0][0]],
       [  Pn.den[0][0],            [1],   (-Pn).den[0][0]]];
G_ = tf(num, den)
G = tf2ss(G_)

#%% H-inf制御器の計算（実行4.4）
#[K,clp,gamma_min,hinf_info] = hinfsyn(G,1,1,'display','on');
K, CL, gam, rcond = control.hinfsyn(G,1,1)
print(gam)

plt.figure(4)
mag1, phase1, om1 = bode( P_nominal, plot=False );
plt.semilogx( om1, 20*np.log10(mag1), '--', label='$P_{nominal}$' );
mag2, phase2, om2 = bode( K, plot=False );
plt.semilogx( om2, 20*np.log10(mag2), label='$K$' );
plt.ylim(-60, 60);
plt.xlim(1e0, 1e2);
plt.legend();

#%% chkperf2.m
#%% 閉ループ特性の確認
L = P_nominal*K;
T = feedback(L,1); #% T = L/(1+L)
S = feedback(1,L); #% S = 1/(1+L)
plt.figure(5);
#bodemag(T.nominal,'-',1/Wt,':',S.nominal,'--',1/(Wps*P.nominal),'-.',w);
#legend('T.nomial','1/Wt','S.nominal','1/(Wps*P)',2)
mag, phase, om = bode( T, w, plot=False );
plt.semilogx( om, 20*np.log10(mag), '-', label='$T$' );
mag, phase, om = bode( 1/Wt, w, plot=False );
plt.semilogx( om, 20*np.log10(mag), ':', label='$1/W_T$' );
mag, phase, om = bode( S, w, plot=False );
plt.semilogx( om, 20*np.log10(mag), '--', label='$S$' );
mag, phase, om = bode( 1/(Wps*P_nominal), w, plot=False );
plt.semilogx( om, 20*np.log10(mag), '-.', label='$1/(W_{PS}*P)$' );
plt.xlim(1e0, 1e2);
plt.ylim(-30, 30);
plt.legend();

delta = np.arange(-1, 1.01 , 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);
    L = P*K;
    T = feedback(L,1); #% T = L/(1+L)
    M = feedback(P,K); #% M = P/(1+L)
    yout, time = step( T, 2 );
    plt.figure(6);
    plt.plot(time, yout);
    yout, time = impulse( M, 2 );
    plt.figure(7);
    plt.plot(time, yout);
実行結果

f:id:manva:20210711235130p:plain
教科書の図4.13。
f:id:manva:20210716224638p:plain
教科書の図4.14。
f:id:manva:20210714223403p:plain
教科書の図4.15。
f:id:manva:20210711235140p:plain
教科書の図4.16。
f:id:manva:20210711235147p:plain
教科書の図4.17。修正混合感度問題では，外乱応答は改善したが指令応答が悪化した。

無料で実践ロバスト制御④ ーH∞制御問題を解くー

制御工学 Python-Control

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

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

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

Amazon

Python-ControlでH∞制御問題が解けたので，どのようにしたかを書いておく。
Python-ControlにはH∞制御問題を解く関数(hinfsyn)が用意されているので，MATLABのソースを少し修正して移植するだけで簡単にできると思ったが，案外苦労した。

H∞制御問題の概要

ここで，H∞制御問題の概要を説明しておく。H∞制御では，制御対象Gを下図のような「一般化プラント」の形で表す。

H∞制御の目的は，外部入力wから制御量zまでの伝達関数をなるべく小さくするような制御器Kを求めることである。おそらく，最初にこの一般化を考えた人は，wを指令値，zを偏差などとして，偏差をなるべく小さくするような制御器Kを求めればよい，と考えたのではなかろうか。だが，問題はそんなに単純ではなかった。zは何でも良いわけではなく，「混合感度問題」となるように決める必要がある。
また，この問題は，解析的に式では解けなかった。だが，コンピュータを使えば解けたも同然の結果が得られる。それは，「wからzまでの伝達関数のH∞ノルム(大きさ)がある定数γより小さくなる制御器Kが存在するかどうか」というのを判定でき，存在するならその制御器を求められる方法が発見されているからである。それさえ判定できるのなら，コンピュータでγをいろいろ変えて何度も計算して一番小さいγを探せる(γ-イタレーションという)。この解き方はなかなか数学的に難解であるが，MATLABやPython-Controlなどの関数になっているので詳細を理解しなくても道具としては使えるのである。

H∞制御系設計手順

MATLABやPython-Controlで関数hinfsynを使えば，以下のような手順でH∞制御問題を解ける。
１．制御対象の摂動を定義する
２．乗法的摂動のボードゲイン線図を描く
３．それを覆うように重みW_Tを決める
４． 感度関数に対する重みW_Sを決める
５．一般化プラントGを定義
６． 解く！ hinfsyn

重み関数をどのように決めればよいか，などは後日説明することにして，まずは上記の解き方の流れを一通り確認する。
前回↓
無料で実践ロバスト制御③ ー摂動を持つ制御対象ー - manvaのエンジニアリング魂
で上記手順１までできたので，今回はその続きで，手順２から６まで(教科書のプログラム4.2, 4.3, 実行4.1, 4.2)を進めた。

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

手順２． 乗法的摂動のボードゲイン線図を描く

プログラム4.2では，まず前回求めた，摂動のある制御対象の摂動部分のみを乗法的摂動の形で取り出し，周波数応答のゲイン線図を描く。(乗法的摂動とは，摂動のある制御対象\tilde{P}を\tilde{P}=(1+\Delta_m)Pの形で表した\Delta_mの部分)。教科書では，不確かさを持つプラントの周波数応答を計算するためにMATLABのufrdというコマンドを使っており，これがPython-Controlにはないが，これは図を描くためだけなので前回までのようにfor文でパラメータを少しずらしながら何回も描けばよい。

手順３． それを覆うように重みW_Tを決める

その後，それを覆う(全ての周波数で少し上になる)ように重み関数W_m(s)を決める。それは教科書ではもう求められている((4.20)式)ので，それをプロットするだけだ。

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

↓は前回のプログラム実行後に実行する。

#%% 乗法的摂動の見積もり
s = tf('s');
k2 = 300;
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_nominal = ss( Ap, Bp, Cp, Dp );
Pn = ss2tf(P_nominal);

plt.figure(2)
delta_ = np.arange(-1, 1.01 , 0.1)
delta = np.delete(delta_,10)

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);
    P_g = ss2tf(P);
    Dm_g = (P_g - Pn)/Pn; #% 乗法的摂動の計算
    #%% ゲイン線図のプロット
    #bodemag(Dm_g,'--',Wm,'r-',w)
    mag, phase, om = bode( Dm_g, logspace(0,3), plot=False );
    plt.semilogx( om, 20*np.log10(mag), '--' );

#%% 重み Wm の定義
Wm = 3*s**2/(s**2+18*s+45**2);
mag2, phase2, om2 = bode( Wm, plot=False );
plt.semilogx( om2, 20*np.log10(mag2) );
実行結果


教科書の図4.5と同じ図が描けた。(第2共振付近がちょっと怪しいが)

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

次に，「プログラム4.3」に移る。

手順４． 感度関数に対する重みW_Sを決める

教科書のこの例では，簡単に(4.22)式のようにほぼ積分のみの関数に決めている。

手順５． 一般化プラントGを定義

以下では，一般化プラントを定義する。
教科書(MATLAB)では，一般化プラントをsysicというコマンドで定義しているが，これがPython-Controlにはない。教科書のやり方を確認すると，

G = sysic;

の後は，ここで求めたGをhinfsynに↓

[K,clp,gamma_min,hinf_info] = hinfsyn(G,1,1,'display','on');

のように代入するだけなので，これをPython-Controlに移植するためにはsysicが何を出力しているか理解する必要があった。Python-Controlのhinfsynのドキュメントでは，Gは「partitioned lti plant」としか書かれていないので困った。
MATLABのsysicについては教科書p27に説明があり，一般化プラントの「状態空間実現」を求めているとのこと。つまり，

\displaystyle{\begin{eqnarray} \dot{x}&=&Ax+B\begin{bmatrix}w\\u\end{bmatrix}\\ \begin{bmatrix}z\\y\end{bmatrix}&=&Cx+D\begin{bmatrix}w\\u\end{bmatrix} \end{eqnarray}}
という形でGを表す行列A,B,C,Dを求めれば良さそう。
Python-Controlでも同じだろうと考え，これを求めようとしたが，それもなかなか厄介。図4.6より，
\displaystyle{\begin{eqnarray} z_1&=&W_S(w_1-u)\\ z_2&=&W_Tu\\ y&=&P(w_1-u)+\epsilon w_2 \end{eqnarray}}
なので，Ws, Wt, Wepsをコマンドssで状態空間モデルにし，

G1 = [ Ws, 0, -Ws ]
G2 = [  0, 0, Wt ]
G3 = [ P, Weps, -P ]
G = [ G1; G2; G3 ]

みたいにつなげれば良いと思ったが，MATLABでは状態空間モデルで上記のような書き方ができる(↓参照)
モデルの結合 - MATLAB & Simulink Example - MathWorks 日本
が，Python-Controlは，series, parallel, feedback, appendは使えるが，[ , ]や[ ; ]に相当するつなぎ方がないっぽい。
結局，状態空間モデルにせず，MIMO(多入力多出力)のtfを使って伝達関数モデルでGを定義し，後から状態空間モデルにしたらできた。MIMOのtfは，3次元の配列(どの入力からどの出力までかが行列のイメージで２次元配列になっており，その各要素が係数のリストなので計３次元)で与える。上記のG1,G2,G3を分子，分母それぞれの係数のリストにした配列を作れば良い。伝達関数モデルで定義したらSISOでも内部的にはnumやdenがリストの２次元配列になっているようなので，「[0][0]」をつけている。

num = [[Ws.num[0][0],    [0], (-Ws).num[0][0]], 
       [         [0],    [0],    Wt.num[0][0]],
       [Pn.num[0][0], [Weps], (-Pn).num[0][0]]]
den = [[Ws.den[0][0],    [1], (-Ws).den[0][0]], 
       [         [1],    [1],    Wt.den[0][0]],
       [Pn.den[0][0],    [1], (-Pn).den[0][0]]]
G_ = tf(num, den)
G = tf2ss(G_)
手順６． 解く！ hinfsyn

ここまでできれば後は関数hinfsyn ↓で解くだけだ。
https://python-control.readthedocs.io/en/0.8.3/generated/control.hinfsyn.html
関数hinfsynへの引数は，(P, nmeas, ncon)の３つあり，Pに一般化プラントGを与える。一般化プラントGへの入力は，図↓

のイメージのように，wとuを(縦にこの順で)並べたベクトル，出力はzとyを並べたベクトルなので，Gへの入力ベクトルのうち下からいくつ目までがuで，Gからの出力ベクトルのうち下からいくつ目までがyかを引数nmeas, nconで与えてやればよい。

移植したコード「プログラム4.3」「実行4.1」「実行 4.2」
import control
from control.matlab import *

#%% defgp.m
#%% 重み関数の定義
s = tf('s');
Ws = 15/(s + 1.5e-2); #% Ws
Wt = Wm;              #% Wt
Weps = 5e-4;          #% Weps

plt.figure(3)
#bodemag(Ws,Wt,'--',w);
mag, phase, om  = bode(Ws, logspace(0,3), plot=False);
mag2,phase2,om2 = bode(Wt, logspace(0,3), plot=False);
plt.semilogx(om,  20*np.log10(mag), label='$W_S$');
plt.semilogx(om2, 20*np.log10(mag2), '--', label='$W_T$');
#legend('Ws','Wt',4);
plt.legend()

#%% 一般化プラントの定義
#Pn = P.nominal;
#systemnames   = 'Pn Ws Wt Weps';
#inputvar      = '[w1; w2; u]';
#outputvar     = '[Ws; Wt; Pn+Weps]';
#input_to_Pn   = '[w1 - u]';
#input_to_Ws   = '[w1 - u]';
#input_to_Wt   = '[ u ]';
#input_to_Weps = '[ w2 ]';
#G = sysic;
num = [[Ws.num[0][0],    [0], (-Ws).num[0][0]], 
       [         [0],    [0],    Wt.num[0][0]],
       [Pn.num[0][0], [Weps], (-Pn).num[0][0]]]
den = [[Ws.den[0][0],    [1], (-Ws).den[0][0]], 
       [         [1],    [1],    Wt.den[0][0]],
       [Pn.den[0][0],    [1], (-Pn).den[0][0]]]
G_ = tf(num, den)
G = tf2ss(G_)

#%% H-inf制御器の計算（実行4.1）
#[K,clp,gamma_min,hinf_info] = hinfsyn(G,1,1,'display','on');
K, CL, gam, rcond = control.hinfsyn(G,1,1)
print( gam )

plt.figure(4)
mag1, phase1, om1 = bode( P_nominal, plot=False );
plt.semilogx( om1, 20*np.log10(mag1), '--', label='$P_{nominal}$' );
mag2, phase2, om2 = bode( K, plot=False );
plt.semilogx( om2, 20*np.log10(mag2), label='$K$' );
plt.ylim(-60, 60);
plt.xlim(1e0, 1e2);
plt.legend()
実行結果


教科書の図4.7。

教科書の図4.8と同じ図が描けた。
γは0.9675で，MATLABで解いた教科書の値0.9693とほぼ同じ。
まあ解けたんでしょう。

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

制御工学 Python-Control

平田氏の著書「実践ロバスト制御」↓を無料の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で目盛りを書いてしまうので最後に
実行結果

f:id:manva:20210626134845p:plain
図4.4とほぼ同じ図が描けた。
つづく。

無料で実践ロバスト制御② ーまずは古典制御ー

制御工学 Python-Control

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

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

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

Amazon
以下これを「教科書」と呼ぶ。

Python-Controlのインストール

Pythonを使うならAnacondaをインストールするのが簡単そう。
「Pythonによる制御工学」↓

Pythonによる制御工学入門

Pythonによる制御工学入門

Amazon
の著者の南氏のサイト↓にAnacondaやPython-Controlのインストールの仕方が書かれており，これ以上書くこともない。
http://y373.sakura.ne.jp/minami/wp-content/uploads/2019/06/pyctrl_anaconda20190603.pdf

私はWindows10のPCを使っているが，問題なくインストールできた。

Jupyter Labの動作確認

上記でJupyter Labがおすすめのようだったので，Jupyter Labを使うことにする。
私の環境ではAnaconda Navigatorから起動しようとしたらエラーが出たが，Anaconda Navigatorを「管理者として実行」したら起動できた。
とりあえずPython-Control本家のサイトのサンプル↓をJupyter Labにコピーして実行してみた。問題なし。
Examples — Python Control Systems Library dev documentation

アップデート

Anaconda NavigatorやJupyter Labのアップデートができない(途中で何も応答なく止まってしまう)問題があったのだが，↓の「インストール済みパッケージの一括更新など」の通り下記を実行したら更新できた。
https://www.kkaneko.jp/tools/win/anaconda3.html

conda config --remove channels conda-forge
conda upgrade -y --all
conda clean -y --packages

実践ロバスト制御開始(まだロバストではない)

準備ができたので，教科書の内容に移る。
まずは1.3節の例で，剛体(図1.4(a))のPD制御の話。
制御対象は剛体

\displaystyle{ y=\frac{1}{Ms^2}u }
として，PD制御
\displaystyle{ u=(k_p+sk_d)(r-y) }
をしたときの指令rから応答yまでの伝達関数を求めると
\displaystyle{ \frac{y}{r}=\frac{\frac{k_d}{M}s+\frac{k_p}{M}}{s^2+\frac{k_d}{M}s+\frac{k_p}{M}} }
となる。この分母が重根\alphaとなるようにゲインk_p, k_dを決めている。
M=1とし，\alphaが-0.5のときと-2.5のときのステップ応答と相補感度関数を求めている(図1.5(a)と(b))。

教科書の「実行1.1」を移植

上記を求めるためのMATLABのサンプルコードが掲載されている(実行1.1)。これをPython-Controlでやってみた。

変更箇所

・最初に，Pythonで使うライブラリをimport

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

・^　→　**
　　Pythonではべき乗は **
・figure( )　→　plt.figure( )
　　図の描画は，matplotlib.pyplotを使う。
・step(Tn1,'--',Tn2,10)　→
　　(yout, T ) = step( Tn1, np.linspace(0, 10, 100) )
　　(yout2, T2) = step( Tn2, np.linspace(0, 10, 100) )
　　plt.plot( T, yout, '--', T2, yout2 )
　　ステップ応答を求める関数はcontrol.matlabのstep。
・上記のようにplotで２つ同時に書いたらシンプルだが，凡例をつけるには，下記のように２つに分けてlabelを設定しておく必要があるようだ。
　→plt.plot( T, yout, '--', label='$\u03b1 = -0.5$')
　　plt.plot( T2, yout2, label='$\u03b1 = -2.5$' )
　　plt.legend()
　　ギリシャ文字のコードは↓
　　【保存版】ギリシャ文字 Latex, Pythonコマンドリスト | OPTY LIFE(オプティライフ)
・ylim([0 1.4])　→　plt.ylim( 0, 1.4 )
・ボード線図は，Python-Controlにはゲインだけプロットする関数がないので，最初，
　　bodemag(Tn1,'--',Tn2)　→　bode(Tn1,Tn2)
　　としてみたが，これだと位相も表示されてしまうし，細かい設定ができなさそう。bodeではplot=Falseとして，グラフは書かず，いったんmagなどの変数に保存し，semilogxで描けば，ゲインだけプロットしたり破線などの設定ができた。
　→mag,phase,om = bode(Tn1, logspace(-1,4), plot=False);
　　mag2,phase2,om2 = bode(Tn2, logspace(-1,4), plot=False);
　　plt.semilogx(om, 20*np.log10(mag), '--');
　　plt.semilogx(om2, 20*np.log10(mag2));

移植したコード「実行1.1」
import numpy as np
import matplotlib.pyplot as plt
from control.matlab import *
s=tf('s')
M = 1
Pn = 1/(M*s**2)
alpha = -0.5;
K1 = alpha**2*M - 2*alpha*M*s/(0.01*s+1);
alpha = -2.5;
K2 = alpha**2*M - 2*alpha*M*s/(0.01*s+1);
Tn1 = feedback(Pn*K1,1);
Tn2 = feedback(Pn*K2,1);

plt.figure(1)
(yout,  T ) = step( Tn1, np.linspace(0, 10, 100) )
(yout2, T2) = step( Tn2, np.linspace(0, 10, 100) )
plt.plot( T, yout, '--', label='$\u03b1 = -0.5$')
plt.plot( T2, yout2, label='$\u03b1 = -2.5$' )
plt.ylim( 0, 1.4 )
plt.legend()

plt.figure(2)
 mag, phase, om = bode(Tn1, logspace(-1,4), plot=False);
mag2,phase2,om2 = bode(Tn2, logspace(-1,4), plot=False);
plt.semilogx( om, 20*np.log10( mag), '--');
plt.semilogx(om2, 20*np.log10(mag2));
実行結果


教科書の図1.5(a)(b)と(だいたい)同じ図を表示できた。
(ちなみに，JupyterLabのグラフをコピーするには，シフトを押しながら右クリック)

古典制御についてちょっと掘下げる

本筋の話とはそれるが，実応用を考える場合，性能をもっと良くするにはどうすればよいか，どれだけ良くできるか，という観点が重要になる。上図の応答は悪くはないが，良くもない。この結果を見て，古典制御だとこれくらいしかできない，と思わないように補足しておく。
まず，以下のイメージを伝えたい。

フィードバックゲインは大きいほど性能が上がるが，大きすぎると不安定になり，発散する
知っている人にとっては当たり前と思うだろうが，教科書(平田氏の著書のことではなく一般の)で理論だけ学んでいると，意外に明確には書かれていなかったりするので，そのせいで何だかピンと来ないまま勉強している初学者も多いのではないか(私がそうだった)。
なぜ明確に書かれていないかというと，理論的に明確には言えないからである。しかし，実験してみると明らかに上記傾向が経験できる。発振限界付近の挙動については，理論と実際が一致しない場合が多い。実際，今回の例のように単純化した剛体の理論では，ゲインを大きくしても不安定にはならず，ほぼパーフェクトな制御が可能である。ちょっといじってパーフェクトな制御にしてみよう。
上記で極\alphaの絶対値を大きくし，-100と-10000にしてみるとステップ応答は以下のようになる。

\alphaが -100のとき(オレンジ実線)はかなりパーフェクトに近く，-10000(青破線)だと発振してイマイチになっている。周波数特性↓を見ると，

カットオフ周波数あたりでピークができており，ゲインを大きくするほどピークが大きくなっている。「フィードバックゲインが大きいほど性能が上がる」と書いたが，ここでは悪化している。これはなぜかというと，微分sというのが厳密には実現できない(「プロパーではない」という)ので，近似微分
\displaystyle{ \frac{s}{0.01s+1} }
にしているためだ。0.01のところを0.0001にすると，


この様になり，だいぶパーフェクトに近づいた。実際には微分は基本的に前回値と今回値の差分から求める(オイラー微分)ため，この0.01や0.0001という数字は，制御周期によって決まる。制御周期を短くするほど性能が上がりそうであることがわかる。

さて、制御周期が十分に速く，微分が近似ではなくsとみなせるとしても，零点があるので周波数特性は0dBより少し持ち上がる。そこで，もっとわかりやすくするために，上記PD制御に以下のように速度FFを加えてみる。

\displaystyle{ u=(k_p+sk_d)(r-y)-s k_dr }
こうすると指令から応答までの伝達関数は
\displaystyle{ \frac{y}{r}=\frac{\frac{k_p}{M}}{s^2+\frac{k_d}{M}s+\frac{k_p}{M}} }
と，分子が定数項だけになり，分母を重根などとおけば単純な2次のローパスフィルタの特性になる。重根の場合，極\alphaの絶対値がカットオフ周波数[rad/s]になり，この極(すなわちフィードバックゲインk_p, k_d)は理論上いくらでも大きくできるのでパーフェクトな制御ができてしまう。実際には，上記のように微分が厳密には実現できないことや，無駄時間，無視した高周波の共振の影響(スピルオーバ)によって，ゲインを上げていくとどこかで発散する。どのくらいのゲインで発散するのかを，理論的に予測することはなかなか難しく，理論と実際が合わないことが多い。この辺がロバスト制御が生まれた背景とも言えるだろう。

MATLABとPython-Controlの比較。

上記のステップ応答をMATLABでプロットしてみると↓のようになる。

横軸は拡大している。Python-Controlではパーフェクトに見えたが，MATLABではオーバーシュートが出ている。MATLABは計算誤差が大きくなりそうなところの刻み幅を自動的に細かくしてくれたりするので，たぶんMATLABの方が正確なんだろう。やはりMATLABは高くても選ばれるだけの価値があるということだ。

無料で実践ロバスト制御① ー挫折からの再入門ー

制御工学 Python-Control

平田氏の著書「実践ロバスト制御」をPython-Controlを使ってやってみようという話。

 

 今まで，ロバスト制御を試してみたいと思いつつ，なかなか難解なので挫折していた。

なぜ挫折してしまうかというと，以下のような要因がある。

  • 数学的に難解
  • 理解するための前提知識が多い
  • 自分の使いたい用途で，解けるための条件を満たさないかも
  • そんなに難しいのに，結果を見ると古典制御に比べてイマイチそう

だが，平田氏の著書↓が大変わかりやすく，タイトル通り「実践」的である。これならできそうと思わせてくれたので再チャレンジしようとしている。

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

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

  • 作者:平田 光男
  • 発売日: 2017/03/24
  • メディア: 単行本
 

 ただ問題があり，この本はMATLABのRobust Control Toolboxを前提としており，それが使えない。MATLABは制御分野でデファクト・スタンダードとなっているために値段が下がらず，高い。かつてのWindowsのようだ。個人向けのMATLAB Homeというもあるが，心意気次第ではギリギリ出せる，というくらいの価格設定で，Simulinkのブロック数制限などもあり，気持ちよく使える感じではない。(実は私もだいぶ前に買って使っているのだが，Robust Control Toolboxは買っておらず，追加するには，1年間の保守期間を過ぎているのでMATLABやControl System Toolboxももう一度買わないといけないのでやめた)。MATLAB高すぎ問題を解決すべく，同等のフリーソフトを作ろうとする試みは「Octave」「Scilab」「FreeMat」など昔から多数あったが，なかなか牙城を崩せずにいた。が，最近，「Python-Control」がかなり有望である。Simulinkの機能はないが，MATLABの機能はかなり揃っており，関数名などもMATLABと互換性の高いMATLAB Compatibility Module↓というのがあり，上記の本で必要なControl System ToolboxやRobust Control Toolboxの機能もほぼ揃いそうである。

MATLAB Compatibility Module — Python Control Systems Library 0.6d documentation

上記の本のサンプルプログラムを無料のPython-Controlに移植して確認しながら勉強していこうと思う。誰かの参考になるかもしれないので途中経過をブログに書きながらやってみる。

ロバスト制御の実用性について，私は少し懐疑的である。本の通り試すだけでなく古典制御と比較してメリットがあるのかも論じてみたい。

私は制御理論は本で勉強することが多いが，いざ学んだ理論を実際に応用しようとするとギャップを感じることも多い。PID制御などの基礎的な話はすでにネットに情報が溢れているので書かないつもりだが，重要度の割に語られ足りてないと感じる話もあるので，気が向いたらそういう話も織り交ぜてみよう。