manvaのエンジニアリング魂

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

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

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

前回までで,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。修正混合感度問題では,外乱応答は改善したが指令応答が悪化した。