混合感度問題 | 制御系CAD

以下では時間関数 $x(t)$ のラプラス変換を $\displaystyle{\hat{x}(s)=\int_0^\infty x(t)e^{-st}\,dt}$ のように表します。

混合感度問題…Homework

[1] モデルの不確かさがある場合の制御系設計のために次を考えます。

図1 摂動前と摂動後の閉ループ系

上図の閉ループ系について、次の伝達関数（インパルス応答のラプラス変換、または入出力のラブラス変換の比）を定義します。

$\displaystyle{(1)\quad\hat{L}(s)=\hat{G}(s)\hat{K}(s)}$

$\displaystyle{(2)\quad\hat{R}(s)=1+\hat{L}(s)}$

$\displaystyle{(3)\quad\boxed{\hat{S}(s)=\frac{1}{\hat{R}(s)}=\frac{1}{1+\hat{L}(s)}}}$

$\displaystyle{(4)\quad\boxed{\hat{T}(s)=\frac{\hat{L}(s)}{\hat{R}(s)}=\frac{\hat{L}(s)}{1+\hat{L}(s)}}}$

ここで、 $\hat{L}(s)$ , $\hat{R}(s)$ , $\hat{S}(s)$ , $\hat{T}(s)$ はそれぞれループ伝達関数、環送差、感度関数、相補感度関数とよばれ、次の制約式が成り立ちます。

$\displaystyle{(5)\quad\boxed{\hat{S}(s)+\hat{T}(s)=1}}$

制御対象の伝達特性の不確かさを表すために次式を考えます。

$\displaystyle{(6)\quad\hat{G}_{\hat{\Delta}}(s)=(1+\hat{\Delta}(s))\hat{G}(s)}$

ここで、 $\hat{\Delta}(s)$ は乗法的摂動と呼ばれます。このとき摂動後のループ伝達関数、環送差、感度関数、相補感度関数を次のように表します。

$\displaystyle{(7)\quad\hat{L}_{{\Delta}}(s)=\hat{G}_{{\Delta}}(s)\hat{K}(s)=(1+\hat{\Delta}(s))\hat{L}(s)}$

$\displaystyle{(8)\quad\hat{R}_{\Delta}}(s)=1+\hat{L}_{{\Delta}}(s)}$

$\displaystyle{(9)\quad\hat{S}_{{\Delta}}(s)=\frac{1}{\hat{R}_{{\Delta}}(s)}}$

$\displaystyle{(10)\quad\hat{T}_{{\Delta}}(s)=\frac{\hat{L}_{{\Delta}}(s)}{\hat{R}_{{\Delta}}(s)}}$

●以下では、ラプラスの $s$ を虚軸上に限定して（ $s=j\omega$ ）議論を進めます。

まず、乗法的摂動ついては、あとで具体例を示しますが、ここでは次を仮定します。

$\displaystyle{(11)\quad\boxed{|\hat{\Delta}(j\omega)|<|W_T(j\omega)|}}$

一般に、乗法的摂動のノルムバンド $|W_T(j\omega)|$ は高周波域で大きくなります。

また、(5)より、次の制約があり、 $\hat{S}(j\omega)$ と $\hat{T}(j\omega)$ は同時には小さくできません。

$\displaystyle{(12)\quad\hat{S}(j\omega)+\hat{T}(j\omega)=1}$

そこで、低周波域で顕著な重み関数 $|\hat{W}_S|$ と高周波域で顕著な重み関数 $|\hat{W}_T|$ を用いて

$\displaystyle{(13)\quad \boxed{|\hat{W}_S(j\omega)\hat{S}(j\omega)|^2+|\hat{W}_T(j\omega)\hat{T}(j\omega)|^2}} }$

の最小化を図ります。この問題は感度関数 $\hat{S}$ と相補感度関数 $\hat{T}$ を同時に考慮しているので、混合感度問題と呼ばれています。

●これは次のようなループ伝達関数の整形をもたらします。すなわち、低周波域で支配的な $|\hat{W}_S(j\omega)|$ を用いてループゲイン $|\hat{L}(j\omega)|$ を大きくします。また、高周波域で支配的な $|\hat{W}_T(j\omega)|$ を用いてループゲイン $|\hat{L}(j\omega)|$ を小さくすることを意味します。また、安定余裕を確保するために中間周波域ではループゲイン $|\hat{L}(j\omega)|$ をなだらかに落とすことが求められます。

図2 ループ整形

[2] 混合感度問題を解くことによりモデルの不確かさを考慮したロバスト制御系を設計することができるのですが、その妥当性をナイキスト軌跡を用いて考えます。

ノミナル安定性　これはナイキストの安定判別法を表しています。

図3 NS条件

ノミナル性能　これはクリティカルポイント $-1+j0$ を半径 $|\hat{W}_S|$ の円領域に拡大し、ナイキスト軌跡がクリティカルポインから十分離れるようにしたものです。

図4 NP条件

ロバスト安定性　これはモデルの不確かさがある場合のナイキスト軌跡に幅をもたせ場合の安定判別法を表しています。

図5 RS条件

ロバスト性能　これはこれはモデルの不確かさがある場合のナイキスト軌跡に幅をもたせ場合に十分性能を持たせることを意味しています。

図6 RP条件

以上、４つの仕様のうち、ロバスト性能が最も強力で、これが満足されれば、ロバスト安定性とノミナル性能は満足されます。ロバスト性能仕様は絶対値和となり、解析的に扱いにくいので、十分条件を考慮して、混合感度問題の定式化が行われているといえます。

図7 ロバスト性能の十分条件

[3] モデルの不確かさの一例として無駄時間の場合を考えます。ここで、無駄時間 $L$ とは

$\displaystyle{(14)\quad y(t)=u(t-L)}$

のような入出力関係を言います。これをラプラス変換すると

$\displaystyle{(15)\quad \hat{y}(s)=\int_0^\infty u(t-L)e^{-st}\,dt=e^{-sL}\hat{u}(s) }$

したがって、伝達関数は

$\displaystyle{(16)\quad \frac{\hat{u}(s)}{\hat{\delta}(s)}=e^{-Ls} }$

となります。これを有理関数近似またはパデー近似(Pade approximation)することがよく行われます。

$\displaystyle{(17)\quad \begin{array}{l} 1+(-Ls)+\frac{1}{2}(-Ls)^2+\cdots\simeq\frac{1+\beta s}{1+\alpha s}\\ \Rightarrow\ 1+\alpha s+(-Ls)(1+\alpha s)+\frac{1}{2}(-Ls)^2+\cdots\simeq 1+\beta s\\ \Rightarrow\ \alpha-L=\beta,\ \frac{1}{2}L^2-\alpha L=0\\ \Rightarrow\ \alpha=\frac{1}{2}L,\beta=-\frac{1}{2}L\\ \Rightarrow\ \boxed{e^{-Ls}\simeq\frac{1-\frac{1}{2}Ls}{1+\frac{1}{2}Ls}} \end{array} }$

これは１次近似ですが、２次近似の場合は同様にして

$\displaystyle{(18)\quad \boxed{e^{-Ls}\simeq\frac{1-\frac{1}{2}Ls-\frac{1}{12}L^2s^2}{1+\frac{1}{2}Ls+\frac{1}{12}L^2s^2}}}$

となります。ステップ応答を比較してみると次のようになります。

図8 無駄時間Pade 近似の比較

それでは入力端に無駄時間がある場合の乗法的摂動を求めてみます。

$\displaystyle{(19)\quad \underbrace{\hat{G}(s)e^{-Ls}}_{\hat{G}_{\hat{\Delta}(s)}}=(1+\underbrace{e^{-Ls}-1}_{\hat{\Delta}(s)})\hat{G}(s) }$

このノルムバンドは次のように見積ることができます。

$\displaystyle{(20)\quad \begin{array}{lll} \underbrace{|e^{-Ls}-1|}_{|\hat{\Delta}(s)|} &=&|(-Ls)+\frac{1}{2}(-Ls)^2+\cdots|\\ &=&|(-Ls)||1+\frac{1}{2}(-Ls)+\cdots| \simeq\underbrace{(1+\epsilon)\left|\frac{s}{L^{-1}}\right|}_{|\hat{W}_T(s)|} \end{array} }$

次は $L=5$ の場合の例です。

図9 無駄時間のノルムバンド

演習B61…Flipped Classroom
$1^\circ$ 船舶の回頭モデル

$\displaystyle{(21)\quad \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1\\ 0 & -\frac{1}{T} \end{array}\right] \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ \frac{K}{T} \end{array}\right] \delta(t-L) }$

に対して、適当なPIDコントローラ

$\displaystyle{(22)\quad \delta(t)=K_P(\psi_c-\psi(t))-K_Dr(t)+\int_0^t(\psi_c-\psi(\tau))d\tau }$

が与えられているとする。次のプログラムを用いて、モデル誤差のノルムバンド（図9相当）、開ループ整形（図2相当）、閉ループ整形（図7参照）を描き、ロバスト性について検討せよ。

図10 開ループ整形と閉ループ整形

MATLAB

%b61.m 
%-----
 clear all, close all
 T1=118; T2=7.8; T3=18.5; T=T1+T2-T3, K=0.185
 zeta=0.5/sqrt(T*K), wn=sqrt(K/T)
 A=[0 1;0 -2*zeta*wn]; B=[0;wn^2]; C=[1 0];
%-----   
 w=logspace(-3,1,300);
 G(:,:)=freqresp(ss(A,B,C,[]),w); G=G.'; 
 ga_G=20*log10(abs(G));
 ph_G=180/pi*atan2(imag(G),real(G));
 lag=10;
 H=exp(-lag*i*w);
 ga_H=20*log10(abs(H));
 ph_H=180/pi*atan2(imag(H),real(H));
 GH=G.*H;
 ga_GH=20*log10(abs(GH));
 ph_GH=180/pi*atan2(imag(GH),real(GH));
 figure(1)
 subplot(211),semilogx(w,ga_G,w,ga_H,'b',w,ga_GH,'r'),grid,legend('G','H','GH')
 subplot(212),semilogx(w,ph_G,w,ph_H,'b',w,ph_GH,'r'),grid,legend('G','H','GH')
 figure(2)
 plot(real(G),imag(G),'o',real(GH),imag(GH),'+',-1,0,'or'),
 axis([-2 2 -2 2]); grid
%-----  
 ga_DL=20*log10(abs(H-1));
 ga_WT=20*log10(abs(lag*i*w));
 figure(3)
 semilogx(w,ga_DL,w,ga_WT'),grid,legend('DL','WT') 
%-----   
 Kp=2.3145921; Kd=46.886428; Ki=0.0097080;
 K=Kp+Kd*(i*w)+Ki./(i*w);
 ga_K=20*log10(abs(K));
 ph_K=180/pi*atan2(imag(K),real(K));
 L=GH.*K;
 ga_L=20*log10(abs(L));
 ph_L=180/pi*atan2(imag(L),real(L));
 WS=Ki*ones(1,length(w))./(i*w);
 WT=lag*i*w;
 ga_WS=20*log10(abs(WS));
 ga_WT=20*log10(abs(WT));
 figure(4)
 %semilogx(w,ga_GH,w,ga_K,w,ga_L,w,ga_WS,w,-ga_WT),grid,legend('GH','K','L','WS','WTi') 
 semilogx(w,ga_L,w,ga_WS,w,-ga_WT),grid,legend('L','WS','WTi') 
%-----   
 S=ones(1,length(w))./(1+L);
 ga_S=20*log10(abs(S));
 ph_S=180/pi*atan2(imag(S),real(S));
 T=L./(1+L);
 ga_T=20*log10(abs(T));
 ph_T=180/pi*atan2(imag(T),real(T));
 WS=Ki*ones(1,length(w))./(i*w);
 WT=lag*i*w;
 ga_WS=20*log10(abs(WS));
 ga_WT=20*log10(abs(WT));
 mu=20*log10(abs(WS.*S)+abs(WT.*T));
 figure(5)
 %semilogx(w,ga_S,w,ga_T,w,-ga_WS,w,-ga_WT,w,mu),grid,legend('S','T','WSi','WTi','mu') 
 semilogx(w,ga_S,w,ga_T,w,-ga_WS,w,-ga_WT),grid,legend('S','T','WSi','WTi') 
%-----  
%eof