「制御技術」カテゴリーアーカイブ

SMCまとめ

投稿日時: 投稿者:

２次安定性

[1]　非線形系

\displaystyle{(1)\quad \dot{x}(t)=f(x(t))}

に対してリャプノフ関数

\displaystyle{(5)\quad V(x)=x^T(t)Px(t)\quad(P>0)}

を考えます。これを微分して(1)を用いると（(1)の解に沿って微分すると言います）

\displaystyle{(6)\quad %\begin{array}{l} \frac{d}{dt}V(x)=\dot{x}^T(t)Px(t)+x^T(t)P\dot{x}(t)=2x^T(t)P\dot{x}(t)=2x^T(t)Pf(x) %\end{array} }

となります。このとき２次安定とは

\displaystyle{(7)\quad { \boxed{\frac{d}{dt}V(x)=2x^T(t)Pf(x)\le -x^T(t)Qx(t)\quad(Q>0)} }}

が成り立つことを言います。これから次式を示すことができ、漸近安定性が成り立ちます。

\displaystyle{(8)\quad \boxed{||x(t)||\le \beta||x(0)||e^{-\frac{1}{2}\alpha t}} }

ただし

\displaystyle{(9)\quad {\alpha=\sigma_n(P^{-1}Q),\ \beta=\sqrt{\frac{\sigma_1(P)}{\sigma_n(P)} }}

[2]　制御対象が平衡状態にあることは線形状態方程式\dot{x}(t)=Ax(t)+Bu(t)において、x=0, u=0を意味します。そこで平衡状態が乱されてx\ne0となる時刻をt=0にとると、線形状態方程式は次式となります。

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)\qquad(x(0)\ne0) }

これに出力方程式

\displaystyle{(2)\quad y(t)=Cx(t) }

を考慮する場合も含めると、漸近安定性の判定法は次のようにまとめられます。

【漸近安定性の定義とその等価な条件】
定義DA:　\forall x(0)\ne 0: x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow\infty)
条件A0:　\exp(At)\rightarrow 0\quad(t\rightarrow\infty)
条件A1:　{\rm Re}(\lambda_i(A))<0\quad(i=1,\cdots,n)
条件A2:　\exists P>0: PA+A^TP+I_n=0
条件A3:　\exists P>0: PA+A^TP+C^TC=0　ただし、(A,C)は可観測対
条件A4:　\exists P>0: PA+A^TP<0

まず(1)が漸近安定とすると、条件A2のリャプノフ方程式から次式を得ます。

\displaystyle{(3)\quad x^T(t)(PA+A^TP)x(t)=-x(t)^Tx(t) }

次に(1)に対してリャプノフ関数

\displaystyle{(4)\quad V(x)=x^T(t)Px(t) }

を考えるとき、(3)は２次安定性の条件

\displaystyle{(5)\quad \begin{array}{l} \frac{d}{dt}V(x)=2x^T(t)P\dot{x}(t)=2x^T(t)PAx(t)=x^T(t)(PA+A^TP)x(t)\\ =-x^T(t)x(t) \end{array} }

を意味します。これから漸近安定な線形系は２次安定であることが分かります。

●条件A2の両辺からQ^{1/2}\ (Q>0)をかけて

\displaystyle{(6)\quad \underbrace{Q^{1/2}PQ^{1/2}}_{P'}\underbrace{Q^{-1/2}AQ^{1/2}}_{A'}+\underbrace{Q^{1/2}A^TQ^{-1/2}}_{A'^T}\underbrace{PQ^{1/2}Q^{1/2}}_{P'}+\underbrace{Q^{1/2}Q^{1/2}}_{Q}=0 }

を得ます。(A',Q^{1/2})は可観測対なので条件A3を適用して（'をはずして）

条件A2′:　\boxed{\exists P>0: PA+A^TP+Q=0}　ただし、Q>0

も漸近安定性の条件となることに注意します。

●たとえば、安定行列\displaystyle{ A= \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }

が与えられるとき、リャプノフ方程式PA+A^TP+I_n=0を満足するP>0は

\displaystyle{ P=\left[\begin{array}{cc} 2.6250 & 0.1250\\ 0.1250 & 0.6250 \end{array}\right] }

のように求まります。これはMATLABで次のコマンドを用いて計算できます。

MATLAB 
P=lyap(A',eye(2))

ここで、P=lyap(A,Q) はリャプノフ方程式PA^T+AP+Q=0の解を求める関数であることに注意します。

問題設定

[3]　制御対象の状態方程式として次式を考えます。

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n, m< n, {\rm rank}B=m) \end{array} }

ここで、f(t,x,u)はモデル誤差、非線形要素、外乱などの影響を表し、有界かつ未知とします。これに対し、次のスイッチング関数を定義します。

\displaystyle{(2)\quad \boxed{s(t)=Sx(t)}\quad(s(t)\in{\rm\bf R}^m, {\rm rank}S=m) }

いま(1)の解が、ある時刻t_sに対して

\displaystyle{(3)\quad s(t)=0\quad(t\ge t_s) }

を満足するとき、スライディングモードが発生していると言います。

●以上の準備の下で、スライディングモード制御問題は次のように記述されます。

問題１: 動的システムの振る舞いが閉じ込められる超平面 {\cal S}=\{x:Sx=0\} を定義するスイッチング関数を決定せよ。
問題２: 有限時刻 t_s において状態を超平面内に拘束し、引き続き留まらせる(\forall t>t_s: s(t)=0)ことのできるスライディングモード制御則（SM制御則）を設計せよ。

●いま、SM制御則を、２次安定性の条件
を、２次安定性

\displaystyle{(4)\quad %\boxed{ \begin{array}{lll} V(\bar{x})= \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} P_1 & 0\\ 0 & P_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ \Rightarrow \dot{V}(\bar{x})\le - \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} Q_1 & 0\\ 0 & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array}}} }

すなわち

\displaystyle{(4')\quad \begin{array}{lll} V(\bar{x})= \underbrace{x_1^T(t)P_1x_1(t)}_{V(x_1)}+\underbrace{s^T(t)P_2s(t)}_{V(s)}\\ \Rightarrow \dot{V}(\bar{x})\le -x_1^T(t)Q_1x_1(t)-s^T(t)s(t) \end{array}} }

が成り立つように決定します（P_1>0, P_2>0, Q_1>0）。

以下で示すように、このSM制御則は次式で表されます。

\displaystyle{(5)\quad { \boxed{\begin{array}{l} u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }

ここで、SBは正則行列であることを仮定しています。また設計パラメータは、スイッチング関数を決める行列S、安定行列\Phi、スカラー\rho>0です。

●２重積分器に対して、S,\Phi,\rhoを適当に選んだときのSM制御則の計算は次のように行われます。

MATLAB 
%cDINT_smc.m
%-----
clear all, close all
A=[0 1;
   0 0];
B=[0;1];
C=eye(2); 
D=zeros(2,1);
M=1;
S=[M 1];
Phi=-1;
P2=lyap(Phi,1);
P2S=P2*S;
Leq=inv(S*B)*S*A;
LPhi=-inv(S*B)*Phi*S;
L=Leq+LPhi;
rho=0.5;
Ln=inv(S*B)*rho;
x0=B;
sim('DINT_smc')
%-----
%eof

これを実行して、図1aのシミュレーション結果を得ます。


図1a　２重積分器のSM制御

図1b　２重積分器の線形制御（SM制御のスイッチング項を外した場合）

図1bに、SM制御のスイッチング項を外した場合のシミュレーション結果を示します。これらからスイッチング項が速応性の改善に役立っていることが分かります。

●いま(1)においてf(t,x,u)=0、すなわち

\displaystyle{(6)\quad \dot{x}(t)=Ax(t)+Bu(t) }

とし、SBの正則性を仮定します。t\ge t_sにおいて、スライディングモード時の制御則は等価制御と呼ばれ、次式のように求められます。

\displaystyle{(7)\quad \begin{array}{l} s(t)=0\\ \Rightarrow\dot{s}(t)=0\\ \Rightarrow S\dot{x}(t)=SAx(t)+SBu(t)=0\\ \Rightarrow \boxed{u_{eq}(t)=-\underbrace{(SB)^{-1}SA}_{L_{eq}}x(t)} \end{array} }

これによる閉ループ系は次式となります。

\displaystyle{(8)\quad \boxed{\dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)}\quad(t\ge t_s) }

●いま(1)においてf(t,x,u)\ne 0の場合、適合サイズをもつ行列Rを用いて

\displaystyle{(9)\quad　{ \boxed{f(t,x,u)=BR\xi(t,x)}} }

のように表されると仮定します。このとき、等価制御は

\displaystyle{(10)\quad \begin{array}{l} s(t)=0\\ \Rightarrow\dot{s}(t)=0\\ \Rightarrow S\dot{x}(t)=SAx(t)+SBu(t)+SBR\xi(t,x)=0\\ \Rightarrow u(t)=-(SB)^{-1}(SAx(t)+SBR\xi(t,x)) \end{array} }

となります。これによる閉ループ系は次式となります。

\displaystyle{(11)\quad \begin{array}{l} \dot{x}(t)=Ax(t)-B(SB)^{-1}(SAx(t)+SBR\xi(t,x))+BR\xi(t,x)\\ =\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)+ \underbrace{(I_n-B(SB)^{-1}S)B}_0R\xi(t,x)\\ =\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t) \end{array} }

すなわちマッチング条件とよばれる(9)が成り立つとき、スライディングモード時はf(t,x,u)の影響を受けないことが分かります。

●剛体振子の運動方程式は次式で与えられました。

\displaystyle{(12)\quad \ddot{\theta}(t)=-\underbrace{\frac{3g}{4\ell}}_{a}\sin\theta(t)+\underbrace{\frac{1}{J}\tau(t)}_{u(t)} }

x_1=\theta, x_2=\dot{\theta}とおいて、次の非線形状態方程式を得ます。

\displaystyle{(13)\quad \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} x_2(t)\\ -a\sin x_1(t)+u(t) \end{array}\right] }

これは、次式のように、２重積分器に非線形項が加えられたものとみなすことができます。

\displaystyle{(14)\quad \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} 0 & 1\\ 0 & 0 \end{array}\right] \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]+ \left[\begin{array}{c} 0\\ 1 \end{array}\right] u(t)+ \left[\begin{array}{c} 0\\ -a\sin x_1(t) \end{array}\right] }

そこで、図1のSM制御則をそのまま用いたときの様子を調べてみます。図2はa=15とした場合（\ell=0.5）のシミュレーション結果です。


図2a　非線形項を考慮した場合のSM制御

図2b　非線形項を考慮した場合の線形制御（SM制御のスイッチング項を外した場合）

図1のSM制御則は非線形項を考慮していないのですが、マッチング条件が満たされているので、非線形項の影響を抑制していることが分かります。これは図2bのスイッチング項がない場合と比べると明らかです。

SM制御則

[4]　制御対象の状態方程式として次式を考えます。

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n) \end{array} }

以下では、状態方程式とスイッチング関数は次のSMC標準形(regular form)をとるように座標変換されていると仮定します（{\rm rank}B={\rm rank}S=m）。

\displaystyle{(2)\quad \boxed{\underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) + \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)}} }

\displaystyle{(3)\quad \boxed{s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} \ (M=S_2^{-1}S_1)} }

これに対して、座標変換

\displaystyle{(4)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array} }

を行って、次式を得ます。

\displaystyle{(5)\quad \boxed{ \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22}\\ \end{array}\right] }_{\bar{A}=T_s A T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \bar{B}_2 \end{array}\right] }_{\bar{B}=T_sB} u(t)\\ + \underbrace{ \left[\begin{array}{c} f_u(t,x)\\ S_1f_u(t,x)+S_2f_m(t,x,u) \end{array}\right] }_{T_sf(t,x,u)} \end{array}} }

ただし

\displaystyle{(6)\quad \left\{\begin{array}{l} \bar{A}_{11}=A_{11}-A_{12}M\\ \bar{A}_{12}=A_{12}S_2^{-1}\\ \bar{A}_{21}=S_2(M\bar{A}_{11}+A_{21}-A_{22}M)\\ \bar{A}_{22}=S_2(M{A}_{12}+A_{22})S_2^{-1}\\ \bar{B}_{2}=S_2B_2 \end{array}\right. }

以下では、\bar{A}_{11}が安定行列となるようにスイッチング関数が選ばれていると仮定します。

●SMC標準形に対する等価制御は、f_u=0, f_m=0として

\displaystyle{ \begin{array}{cl} & s(t)=0\ \Rightarrow\ \dot{s}(t)=0\\ \Downarrow & by\ (5)\\ &\bar{A}_{21}x_1(t)+\bar{A}_{22}s(t)+\bar{B}_{2}u(t)=0\\ \Downarrow & \\ (7) & u_{eq}(t)=-\underbrace{\bar{B}_{2}^{-1}}_{([0\ I_m]T_sB)^{-1}} \underbrace{\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]}_{[0\ I_m]T_sAT_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ \Downarrow & by\ S=[0\ I_m]T_s, x=T_s^{-1}\bar{x}\\ (8) & \boxed{u_{eq}(t)=-\underbrace{(SB)^{-1}SA}_{L_\ell}x(t)}} \end{array} }

のように得られます。線形制御則u_\ellは、この等価制御をベースして

\displaystyle{ \begin{array}{cl} (9) & u_{eq}(t)=-\underbrace{\bar{B}_{2}^{-1}}_{([0\ I_m]T_sB)^{-1}} (\underbrace{\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]}_{[0\ I_m]T_sAT_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}-\Phi s(t))\\ \Downarrow & by\ S=[0\ I_m]T_s, x=T_s^{-1}\bar{x}\\ (10) & \boxed{u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)}} \end{array} }

のように構成します。ここで、\Phiは適当に設定された安定行列です。

●このとき閉ループ系は次式で与えられます。

\displaystyle{(11)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{s}(t) \end{array}\right] = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ 0 & \Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ + \left[\begin{array}{c} f_u(t,x)\\ \bar{B}_{2}u_n(t)+S_2(Mf_u(t,x)+f_m(t,x,u)) \end{array}\right] \end{array} }

これからスライディングモードs=0を達成するために\Phiの影響を与えることが分かります。簡単な例で確かめると\Phiはダンピングを与え、速くスライディングモードを達成するためのパラメータであることが分かります。

さて、\Phiは安定行列なので

\displaystyle{(12)\quad \begin{array}{l} \boxed{P_2\Phi+\Phi^TP_2=-I} \end{array} }

を満たすP_2>0を選ぶことができます。これを用いて、スイッチング制御則u_nを

\displaystyle{(13)\quad u_n(t)=-\underbrace{\bar{B}_{2}^{-1}}_{(SB)^{-1}}\rho(t,x)\frac{P_2s(t)}{||P_2s(t)||} }

すなわち

\displaystyle{(14)\quad \boxed{u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}} }

と定めます。線形制御則(10)とスイッチング制御則(14)を用いて、２次安定性を示すことができます。

スイッチング関数（１）

[5]　スライディングモードの場合、s(t)=0なので

\displaystyle{(1)\quad s(t)=0 \Rightarrow Sx(t)=S_1x_1(t)+S_2x_2(t)=0\Rightarrow x_2(t)=-\underbrace{S_2^{-1}S_1}_Mx_1(t)}}

となって、x_2の振る舞いはx_1によって決まります。一方、x_1の振る舞いは

\displaystyle{(2)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}M)}_{\bar{A}_{11}}x_1(t)}

によって決まります。これは

\displaystyle{(3)\quad \dot{x}_1(t)=A_{11}x_1(t)+A_{12}x_2(t)}

に対する安定化状態フィードバック

\displaystyle{(4)\quad x_2(t)=-Mx_1(t)}

による閉ループ系とみなすことができます。

したがって、\bar{A}_{11}が安定行列となるようにスイッチング関数が選ぶためには、(3)に対する状態フィードバックによる安定化問題を解けばよいことが分かります。

●そこで、次の２次形式評価関数を最小にするようにMを決定する方法が提案されています。

\displaystyle{(5)\quad \boxed{J=\int_0^\infty \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]^T }_{x^T(t)} \underbrace{ \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} dt} }

これは

\displaystyle{(6)\quad J=\int_0^\infty(x_1^T(t)\underbrace{(Q_{11}-Q_{12}Q_{22}^{-1}Q_{21})}_{\hat{Q}}x_1(t)}
\displaystyle{+\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))^T}_{\hat{u}^T}\underbrace{Q_{22}}_{\hat{R}}\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))}_{\hat{u}})dt}

と書けます。ここで、(3)と(4)をそれぞれ

\displaystyle{(7)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}Q_{22}^{-1}Q_{21})}_{\hat{A}}x_1(t)+\underbrace{A_{12}}_{\hat{B}}\underbrace{(Q_{22}^{-1}Q_{21}x_1(t) +x_2(t))}_{\hat{u}(t)}}

\displaystyle{(8)\quad \underbrace{Q_{22}^{-1}Q_{21}x_1(t) +x_2(t)}_{\hat{u}(t)}=-\underbrace{(M-Q_{22}^{-1}Q_{21})}_{\hat{F}}x_1(t)}

と書き直すと、次のような最適化制御問題の定式化ができていることがわかります。すなわち

\displaystyle{(9)\quad \boxed{\begin{array}{l} \dot{x}_1(t)=\hat{A}x_1(t)+\hat{B}\hat{u}(t)\\ \hat{A}=A_{11}-A_{12}Q_{22}^{-1}Q_{21}\\ \hat{B}=A_{12} \end{array}} }

を２次系式評価関数

\displaystyle{(10)\quad \boxed{\begin{array}{l} J=\int_0^\infty(x_1^T(t)\hat{Q}x_1(t)+\hat{u}^T(t)\hat{R}\hat{u}(t))dt\\ \hat{Q}=Q_{11}-Q_{12}Q_{22}^{-1}Q_{21}\\ \hat{R}=Q_{22} \end{array}} }

を最小にするように

\displaystyle{(11)\quad \boxed{\begin{array}{l} \hat{u}(t)=-\hat{F}x_1(t)\\ \hat{F}=M-Q_{22}^{-1}Q_{21} \end{array}} }

を求める問題です。この解は、リッカチ方程式

\displaystyle{(12)\quad \boxed{\Pi\hat{A}+\hat{A}^T\Pi-\Pi\hat{B}\hat{R}^{-1}\hat{B}^T\Pi+\hat{Q}=0}}

を解いて

\displaystyle{(13)\quad \boxed{\hat{F}=\hat{R}^{-1}\hat{B}^T\Pi}}

を求め、次式からMを決めればよいことがわかります。

\displaystyle{(14)\quad \boxed{M=\hat{F}+Q_{22}^{-1}Q_{21}}}

●２重積分器に対して、関数swflqrによるスイッチング関数Sの決定は次のように行われます。

MATLAB 
%cDINT_smc1.m
%-----
clear all, close all
A=[0 1;
   0 0];
B=[0;1];
C=eye(2); 
D=zeros(2,1);
Mx=0.5; Tx=0.5;
Q=diag([1/Mx,Tx/Mx].^2);
S=swflqr(A,B,Q)
Phi=-1;
P2=lyap(Phi,1);
P2S=P2*S;
Leq=inv(S*B)*S*A;
LPhi=-inv(S*B)*Phi*S;
L=Leq+LPhi;
rho=0.5;
Ln=inv(S*B)*rho;
x0=B;
sim('DINT_smc')
%-----
%eof

これを実行して図3のシミュレーション結果を得ます。


図3　関数swflqrによるスイッチング関数Sの決定

ここで、関数swflqrは、Edwards and Spurgeon によるものです。

スイッチング関数（２）

[6]　等価制御

\displaystyle{(1)\quad u_{eq}(t)=-(SB)^{-1}SAx(t) }

による閉ループ系は次式で表されました。

\displaystyle{(2)\quad \dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)\quad(t\ge t_s) }

ここで

\displaystyle{(3)\quad SA_{eq}=S(I_n-B(SB)^{-1}S)A=0 }

また、A、B、Sは、SMC標準形で与えられるとすると

\displaystyle{(4)\quad \begin{array}{l} A_{eq}= \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & 0 \\ \end{array}\right] }_{I_n-B(SB)^{-1}S} \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right]= \left[\begin{array}{cc} A_{11} & A_{12} \\ -MA_{11} & -MA_{12} \\ \end{array}\right]\\ = \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right] \left[\begin{array}{cc} \bar{A}_{11} & A_{12} \\ 0 & 0 \\ \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right]^{-1} \end{array} }

となって、\bar{A}_{11}が安定行列なので、A_{eq}はm個の零固有値を持ちます。一方、A_{eq}の非零固有値\lambda_i\ (i=1,\cdots,n-m)に対応する固有ベクトルをv_iとすると、

\displaystyle{(5)\quad SA_{eq}=0\Rightarrow SA_{eq}v_i=0\Rightarrow \lambda_iSv_i=0 \Rightarrow Sv_i=0 }

すなわち

\displaystyle{(6)\quad { \boxed{S[v_1\cdots v_{n-m}]=0 %　\Rightarrow [v_1\cdots v_{n-m}]^TS^T=0}} }

が成り立ちます。したがって、まず等価制御による閉ループ系のA_{eq}が望ましい固有値・固有ベクトルをもつように設計して、(5)からSを決定することが考えられます。

●２重積分器に対して、関数swfvplによるスイッチング関数Sの決定は次のように行われます。

MATLAB 
%cDINT_smc2.m
%-----
clear all, close all
A=[0 1;
   0 0];
B=[0;1];
C=eye(2); 
D=zeros(2,1);
lambda=-1;
nocomp=0;
specpos=rand(2,1)
specent=rand(2,1)
S=swfvpl(A,B,lambda,nocomp,specpos,specent)
Phi=-1;
P2=lyap(Phi,1);
P2S=P2*S;
Leq=inv(S*B)*S*A;
LPhi=-inv(S*B)*Phi*S;
L=Leq+LPhi;
rho=0.5;
Ln=inv(S*B)*rho;
x0=B;
sim('DINT_smc')
%-----
%eof

これを実行して図4のシミュレーション結果を得ます。


図4　関数swfvplによるスイッチング関数Sの決定

ここで、関数swfvplは、Edwards and Spurgeon によるものです。

積分動作の導入

[6] 定値外乱を受ける制御対象

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+w\\ y(t)=Cx(t)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, y(t)\in{\rm\bf R}^p, w\in{\rm\bf R}^n, m=p) \end{array} }

の出力を、次のコマンド（定値目標）

\displaystyle{(2)\quad r\in{\rm\bf R}^m %\begin{array}{l} %\dot{r}(t)=\Gamma(r(t)-r_c)\\ %(r(t),r_c\in{\rm\bf R}^m) %\end{array} %\Gammaは安定行列 }

に追従させることを考えます。そのために、積分動作

\displaystyle{(3)\quad \begin{array}{l} \dot{x}_I(t)=y(t)-r\\ %\dot{x}_I(t)=r(t)-y(t)\\ (x_I(t)\in{\rm\bf R}^m) \end{array} }

を考え、次の拡大系を構成します。

偏差系E3：
\displaystyle{(4)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E3}} {\dot u}(t) }

を得ます。この偏差系E3が標準形となっていることに注意して、SMCを設計します。

(4)を、改めて次のように書きます。

\displaystyle{(5)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ I_m \end{array}\right] }_{B_{E3}} {\dot u}(t)\\ (x_1(t)=x(t)-x_\infty, x_2(t)=u(t)-u_\infty) \end{array} }

スイッチング関数として、次式を考えます。

\displaystyle{(6)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} \ (M=S_2^{-1}S_1) }

(15)に対して、座標変換

\displaystyle{(7)\quad \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} }

を行って、次式を得ます。

\displaystyle{(8)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{x}'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right] }_{T_sA_{E3}T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} + \underbrace{ \left[\begin{array}{cc} 0\\ S_2 \end{array}\right] }_{T_sB_{E3}} {\dot u}(t) }

\displaystyle{(9)\quad \left\{\begin{array}{l} \bar{A}_{11}=A-BM\\ \bar{A}_{12}=BS_2^{-1}\\ \bar{A}_{21}=S_1(A-BM)\\ \bar{A}_{22}=S_1BS_2^{-1} \end{array}\right. }

以下では、\bar{A}_{11}が安定行列となるようにスイッチング関数が選ばれていると仮定します。

このとき、SM制御則を、２次安定性

\displaystyle{(10)\quad \begin{array}{lll} V(s)=s(t)^TP_2s(t)& \Rightarrow \dot{V}(s)\le -s^T(t)s(t)&(t\le t_s)\\ V(x_1)=x_1^T(t)P_1x_1(t)& \Rightarrow \dot{V}(x_1)\le -x_1^T(t)Q_1x_1(t)&(t> t_s) \end{array} }

が成り立つように決定します（P_1>0, P_2>0, Q_1>0）。

偏差系E3に対するSMCは次式で与えられます。

\displaystyle{(11)\quad {\dot u}(t)={\dot u}_\ell(t)+{\dot u}_n(t) }

\displaystyle{(12)\quad {\dot u}_\ell(t)=-(SB_{E3})^{-1}(SA_{E3}-\Phi S)x_{E3}(t) }

\displaystyle{(13)\quad {\dot u}_n(t)=-(SB_{E3})^{-1}\rho\frac{P_2s(t)}{||P_2s(t)||}}=-(SB_{E3})^{-1}\rho\,{\rm sgn}(P_2Sx_{E3}(t)) }

これらを積分して、制御対象(1)に対する積分動作をもつSMCを導出します。

\displaystyle{(14)\quad {u}(t)={u}_\ell(t)+{u}_n(t) }

まず(12)は次式のように書けます。

\displaystyle{(15)\quad {\dot u}_\ell(t)=- \boxed{\underbrace{ (SB_{E3})^{-1}(SA_{E3}-\Phi S)S_E^{-1} }_{\left[\begin{array}{cc} F & F_I \end{array}\right]}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)} }

これを積分して

\displaystyle{(16)\quad \boxed{u_\ell(t)=-Fx(t)+F_I\int_0^t(r-y(\tau))d\tau} }

次に(13)は次式のように書けます。

\displaystyle{(17)\quad {\dot u}_n(t) =-(SB_{E3})^{-1}\rho\, {\rm sgn}( \boxed{\underbrace{ P_2SS_E^{-1} }_{\left[\begin{array}{cc} G & G_I \end{array}\right]}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)}) }

これを積分すれば

\displaystyle{(18)\quad \boxed{u_n(t)=(SB_{E3})^{-1}\rho\,\int_0^t\underbrace{{\rm sgn}(-G{\dot x}(\tau)+G_I(r-y(\tau)))}_{0,\pm 1}d\tau} }

●２重積分器に対して、次のように行われます。

MATLAB 
%cDINT_smci.m
%-----
clear all, close all
A=[0 1;
   0 0];
B=[0;1];
C=[1 0]; 
D=0;
n=size(A,1); m=size(B,2); p=size(C,1); 
AE=[A B;zeros(m,n+m)];
BE=[zeros(n,m);eye(m)];
SE=[A B;C D];
QE=eye(n+m);
S=-swflqr(AE,BE,QE)
Phi=-0.5;
Leq=inv(S*BE)*S*AE;
LPhi=-inv(S*BE)*Phi*S;
L=(Leq+LPhi)/SE;
F=L(1:2)
FI=L(3)
P2=lyap(Phi,1)
P2S=P2*S/SE;
G=P2S(1:2)
GI=P2S(3) 
rho=1;
Ln=inv(S*BE)*rho
x0=zeros(2,1);
sim('DINT_smci')
%-----
%eof

これを実行して図5のシミュレーション結果を得ます。


図5　積分動作を導入したSM制御

補遺：2次安定系のロバスト性…Homework

いま漸近安定な線形系

\displaystyle{(1)\quad \dot{x}(t)=Ax(t) }

において次のようなモデルの不確かさがあるとします。

\displaystyle{(2)\quad \dot{x}(t)=(A+\Delta A)x(t)=Ax(t)+\underbrace{\Delta Ax(t)}_{\xi(t,x)}\qquad(x(0)\ne0) }

条件A2’より、適当なP>0とQ>0に対して

\displaystyle{(3)\quad PA+A^TP+Q=0 }

が成り立ちます。このとき

\displaystyle{(4)\quad \boxed{||\xi(t,x)||\le\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||} }

ならば、リャプノフ関数

\displaystyle{(5)\quad V(x)=x^T(t)Px(t) }

に対して

\displaystyle{(6)\quad \frac{d}{dt}V(x)<0 }

となって、リャプノフ安定となります。

●実際、

\displaystyle{(7)\quad \begin{array}{l} \frac{d}{dt}V(x)=x^T(t)P\dot{x}(t)+\dot{x}^T(t)Px(t)\\ =x^T(t)P(Ax(t)+\xi(t,x))+(x^T(t)A^T+\xi^T(t,x))Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)P\xi(t,x)+\xi^T(t,x)Px(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi^T(t,x)Px(t)\\ =-x^T(t)Qx(t)+2\xi^T(t,x)Px(t)\\ \le -x^T(t)Qx(t)+ 2||\xi(t,x)||\underbrace{||Px(t)||}_{\sqrt{x^T(t)P^2x(t)}} \end{array} }

ここで、

\displaystyle{(8)\quad \begin{array}{l} \sigma_n(Q)x^T(t)x(t)\le x^T(t)Qx(t) \Rightarrow \sigma_n(Q)||x(t)||^2\le x^T(t)Qx(t)\\ x^T(t)P^2x(t)\le\sigma_1(P^2)x(t)^Tx(t) \Rightarrow \sqrt{x^T(t)P^2x(t)}\le\sigma_1(P)||x(t)|| \end{array} }

に注意して、(4)の下では

\displaystyle{(9)\quad \begin{array}{l} \frac{d}{dt}V(x) \le -\sigma_n(Q)||x(t)||^2+ 2||\xi(t,x)||\sigma_1(P)x(t)^Tx(t) \\ =(||\xi(t,x)||-\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||)2\sigma_1(P)||x(t)||<0 \end{array} }

となります。

補遺：偏差系の安定化SM制御

投稿日時: 投稿者:


補遺：偏差系の安定化SM制御…Homework

レギュレータ問題は平衡状態の安定化、追値問題は他の平衡状態への安定化を意味します。LQI制御のところでみたように、偏差系を導入すれば、追値問題は偏差系の安定化問題となります。この観点から積分動作導入による追従SM制御を見直してみます。

[1] 定値外乱を受ける制御対象

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+w\\ y(t)=Cx(t)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, y(t)\in{\rm\bf R}^p, w\in{\rm\bf R}^n, m=p) \end{array} }

の出力を、次のコマンド（定値目標）

\displaystyle{(2)\quad r\in{\rm\bf R}^m %\begin{array}{l} %\dot{r}(t)=\Gamma(r(t)-r_c)\\ %(r(t),r_c\in{\rm\bf R}^m) %\end{array} %\Gammaは安定行列 }

に追従させることを考えます。そのために、積分動作

\displaystyle{(3)\quad \begin{array}{l} \dot{x}_I(t)=y(t)-r\\ %\dot{x}_I(t)=r(t)-y(t)\\ (x_I(t)\in{\rm\bf R}^m) \end{array} }

を考え、次の拡大系を構成します。

\displaystyle{(4)\quad %\underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] %}_{\dot{x}_E(t)} = %\underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] %}_{A_E} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] %}_{B_E} u(t) + %\underbrace{ \left[\begin{array}{c} w \\ -r \end{array}\right] %}_{w_E} }

定常状態では

\displaystyle{(5)\quad \left[\begin{array}{c} 0 \\ 0 \end{array}\right] = \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] \left[\begin{array}{c} x_\infty \\ x_{I\infty} \end{array}\right] + \left[\begin{array}{c} B \\ 0 \end{array}\right] u_\infty + \left[\begin{array}{c} w \\ -r \end{array}\right] }

を得ます(x_{\infty},x_{I\infty},u_{\infty}は定数ベクトル)。まず、(4)から(5)を引いて、つぎの偏差系を得ます。

偏差系E1：
\displaystyle{(6)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-x_\infty \\ x_I(t)-x_{I\infty} \end{array}\right] }_{\dot{x}_{E1}(t)} = \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E1}} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ x_I(t)-x_{I\infty} \end{array}\right] }_{x_{E1}(t)} + \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E1}} (u(t)-u_\infty) }

この両辺を微分すれば、状態変数の中の定数ベクトルを除くことができて

偏差系E2：
\displaystyle{(7)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{\dot{x}_{E2}(t)} = \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E2}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)} + \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E2}} {\dot u}(t) }

を得ます。さらに、(1)の状態方程式と観測方程式をまとめた

\displaystyle{(8)\quad \left[\begin{array}{c} {\dot x}(t)-w \\ y(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }

から、(5)すなわち

\displaystyle{(9)\quad \left[\begin{array}{c} -w \\ r \end{array}\right] = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} \left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right] }

を引いて、つぎの関係式が成り立ちます。

\displaystyle{(10)\quad \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{x_{E3}(t)} }

これを用いて、偏差系E2に座標変換を行えば

偏差系E3：
\displaystyle{(11)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E3}} {\dot u}(t) }

を得ます。ここで、つぎの関係式を用いました。

\displaystyle{(12)\quad \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E2}} \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} }

\displaystyle{(13)\quad \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E2}} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E3}} }

\displaystyle{(14)\quad \underbrace{ \left[\begin{array}{cc} 0 & I_m \end{array}\right] }_{C_{E2}} \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S_E} = \underbrace{ \left[\begin{array}{cc} C & 0 \end{array}\right] }_{C_{E3}} }

以下では、この偏差系E3が標準形となっていることに注意して、SMCを設計します。

(11)を、改めて次のように書きます。

\displaystyle{(15)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ I_m \end{array}\right] }_{B_{E3}} {\dot u}(t)\\ (x_1(t)=x(t)-x_\infty, x_2(t)=u(t)-u_\infty) \end{array} }

スイッチング関数として、次式を考えます。

\displaystyle{(16)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} \ (M=S_2^{-1}S_1) }

(15)に対して、座標変換

\displaystyle{(17)\quad \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} }

を行って、次式を得ます。

\displaystyle{(18)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{x}'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right] }_{T_sA_{E3}T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} + \underbrace{ \left[\begin{array}{cc} 0\\ S_2 \end{array}\right] }_{T_sB_{E3}} {\dot u}(t) }

\displaystyle{(19)\quad \begin{array}{l} \bar{A}_{11}=A-BM\\ \bar{A}_{12}=BS_2^{-1}\\ \bar{A}_{21}=S_1(A-BM)\\ \bar{A}_{22}=S_1BS_2^{-1} \end{array} }

以下では、\bar{A}_{11}が安定行列となるようにスイッチング関数が選ばれていると仮定します。

このとき、スライディングモード制御則（SM制御則、SMC則）

\displaystyle{(20)\quad {\dot u}(t)=\underbrace{{\dot u}_\ell(t)}_{linear\ control}+\underbrace{{\dot u}_n(t)}_{switching\ component}} }

を、２次安定性

\displaystyle{(21)\quad \begin{array}{lll} V(\bar{x})= \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} P_1 & 0\\ 0 & P_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)}\\ \Rightarrow \dot{V}(\bar{x})\le - \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} Q_1 & 0\\ 0 & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array}} }

すなわち

\displaystyle{(21')\quad \begin{array}{lll} V(\bar{x})= \underbrace{x_1^T(t)P_1x_1(t)}_{V(x_1)}+\underbrace{s^T(t)P_2s(t)}_{V(s)}\\ \Rightarrow \dot{V}(\bar{x})\le -x_1^T(t)Q_1x_1(t)-s^T(t)s(t) \end{array}} }

が成り立つように決定します（P_1>0, P_2>0, Q_1>0）。

[2] 可到達性の検討

等価制御は

\displaystyle{(22)\quad \begin{array}{l} s(t)=0\Rightarrow\dot{s}(t)=0 \Rightarrow 0=\bar{A}_{21}x_1(t)+\bar{A}_{22}s(t)+S_2{\dot u}(t)\\ \Rightarrow {\dot u}_{eq}(t)=-\underbrace{S_2^{-1}}_{(SB_{E3})^{-1}} \underbrace{\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]x'_{E3}(t)}_{SA_{E3}x_{E3}(t)}} \end{array} }

のように得られます。(20)の第１項{\dot u}_\ellは、この等価制御をベースして

\displaystyle{(23)\quad \begin{array}{l} {\dot u}_\ell(t)=-\underbrace{S_2^{-1}}_{(SB_{E3})^{-1}} \underbrace{(\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]x'_{E3}(t)-\Phi s(t))}_{(SA_{E3}-\Phi S)x_{E3}(t)}} \end{array}} }

のように構成します（\Phiは安定行列）。このとき閉ループ系は次式で与えられます。

\displaystyle{(24)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{s}(t) \end{array}\right] = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ 0 & \Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] + \left[\begin{array}{c} 0\\ S_2{\dot u}_n(t) \end{array}\right] \end{array }

すなわち

\displaystyle{(25)\quad \begin{array}{l} \dot{x}_1(t)=\bar{A}_{11}x_1(t)+\bar{A}_{12}s(t)\\ \dot{s}(t)=\Phi s(t)+S_2{\dot u}_n(t) \end{array} }

ここで、\Phiは安定行列なので

\displaystyle{(26)\quad \begin{array}{l} P_2\Phi+\Phi^TP_2=-I \end{array} }

を満たすP_2>0を選ぶことができます。これを用いて

\displaystyle{(27)\quad {\dot u}_n(t)=-\underbrace{S_2^{-1}}_{(SB_{E3})^{-1}}\rho\frac{P_2s(t)}{||P_2s(t)||}} }

と選びます（\rho>0は定数）。このとき次式が成り立ちます。

\displaystyle{(28)\quad \begin{array}{l} \dot{V}(s)=2s^T(t)P_2\dot{s}(t)\\ =2s^T(t)P_2(\Phi s(t)-\rho\frac{P_2s(t)}{||P_2s(t)||})\\ =s^T(t)(P_2\Phi+\Phi^TP_2)s(t)+2s^T(t)P_2(-\rho\frac{P_2s(t)}{||P_2s(t)||})\\ = -||s(t)||^2-2\rho||P_2s(t)||\\ \le -||s(t)||^2 \end{array} }

[3] スライディングモードの検討

\bar{A}_{11}は安定行列なので

\displaystyle{(29)\quad \begin{array}{l} P_1\bar{A}_{11}+\bar{A}_{11}^TP_1=-Q_1<0 \end{array} }

を満たすP_1>0を選ぶことができます。

\displaystyle{(30)\quad \begin{array}{l} \dot{V}(x_1)=2x_1^T(t)P_1\dot{x}_1(t)\\ =2x_1^T(t)P_1( \bar{A}_{11}x_1(t)+\bar{A}_{12}\underbrace{s(t)}_{0})\\ =x_1^T(t)(P_1\bar{A}_{11}+\bar{A}_{11}^TP_1)x_1(t)\\ =-x_1^T(t)Q_1x_1(t) \end{array} }

[4] 積分動作をもつSMC

上で求めた偏差系E3に対するSMCは次式で与えられました。

\displaystyle{(31)\quad {\dot u}(t)={\dot u}_\ell(t)+{\dot u}_n(t) }

\displaystyle{(32)\quad {\dot u}_\ell(t)=-(SB_{E3})^{-1}(SA_{E3}-\Phi S)x_{E3}(t) }

\displaystyle{(33)\quad {\dot u}_n(t)=-(SB_{E3})^{-1}\rho\frac{P_2s(t)}{||P_2s(t)||}}=-(SB_{E3})^{-1}\rho\,{\rm sgn}(P_2Sx_{E3}(t)) }

これらを積分して、制御対象(1)に対する積分動作をもつSMCを導出します。

\displaystyle{(34)\quad {u}(t)={u}_\ell(t)+{u}_n(t) }

まず(32)は(10)を用いて次式のように書けます。

\displaystyle{(35)\quad {\dot u}_\ell(t)=- \underbrace{ (SB_{E3})^{-1}(SA_{E3}-\Phi S)S_E^{-1} }_{\left[\begin{array}{cc} F & F_I \end{array}\right]} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)} }

これを積分して

\displaystyle{(36)\quad u_\ell(t)=-Fx(t)+F_I\int_0^t(r-y(\tau))d\tau }

次に(33)は(10)を用いて次式のように書けます。

\displaystyle{(37)\quad {\dot u}_n(t) =-S_2^{-1}\rho\, {\rm sgn}( \underbrace{ P_2SS_E^{-1} }_{\left[\begin{array}{cc} G & G_I \end{array}\right]} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)}) }

これを積分すれば

\displaystyle{(38)\quad u_n(t)=S_2^{-1}\rho\,\int_0^t\underbrace{{\rm sgn}(-G{\dot x}(\tau)+G_I(r-y(\tau)))}_{0,\pm 1}d\tau }

ここで、\dot{x}_1(t)=\frac{d}{dt}(x(t)-x_\infty)=\dot{x}(t)に注意し、(25)の第１式を用いて\dot{x}(t)は

\displaystyle{(39)\quad \begin{array}{l} \underbrace{\dot{x}_1(t)}_{\dot{x}(t)}=\underbrace{(A-BM)}_{\bar{A}_{11}}x_1(t)+\underbrace{BS_2^{-1}}_{\bar{A}_{12}}\underbrace{(S_1x_1(t)+S_2x_2(t))}_{s(t)}\\ =(A-BM)x_1(t)+BS_2^{-1}S_2(Mx_1(t)+x_2(t))\\ =Ax_1(t)+Bx_2(t)\\ =Ax(t)+Bu(t)-\left[\begin{array}{cc} A & B \end{array}\right] \left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right]\\ =Ax(t)+Bu(t)-\left[\begin{array}{cc} A & B \end{array}\right]S_E^{-1} \left[\begin{array}{c} -w \\ r \end{array}\right]\\ =Ax(t)+Bu(t)-\left[\begin{array}{cc} I & 0 \end{array}\right] \left[\begin{array}{c} -w \\ r \end{array}\right]\\ =Ax(t)+Bu(t)+w \end{array} }

となって元の状態方程式となりますが、wを無視し、uの近似値を使うことも一手段かもしれません。

[5] 数値例(1)

\displaystyle{(101)\quad \begin{array}{l} \dot{x}(t)=\underbrace{0}_{a}x(t)+\underbrace{1}_{b}u(t)+w\\ y(t)=\underbrace{1}_{c}x(t)\\ \end{array} }

\displaystyle{(109)\quad \left[\begin{array}{c} -w \\ r \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{S_E} \left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right] \Rightarrow \left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right]= \left[\begin{array}{c} r \\ -w \end{array}\right] }

\displaystyle{(111)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_{E3}} {\dot u}(t) }

\displaystyle{(116)\quad s(t)= \underbrace{ \left[\begin{array}{cc} s_1 & s_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{s_2 \left[\begin{array}{cc} m & 1 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} \ (m=s_2^{-1}s_1) }

\displaystyle{(117)\quad \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ s_1 & s_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ -s_2^{-1}s_1 & s_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} }

\displaystyle{(118)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{x}'_{E3}(t)} = \underbrace{ \left[\begin{array}{cc} -m & s_2^{-1} \\ -s_1m & m \\ \end{array}\right] }_{T_sA_{E3}T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} + \underbrace{ \left[\begin{array}{cc} 0\\ s_2 \end{array}\right] }_{T_sB_{E3}} {\dot u}(t) }

\displaystyle{(123)\quad \begin{array}{l} {\dot u}_\ell(t)= -\underbrace{s_2^{-1}}_{(SB_{E3})^{-1}} \underbrace{(\left[\begin{array}{cc} -s_1m & m \\ \end{array}\right]x'_{E3}(t)-\Phi s(t))}_{(SA_{E3}-\Phi S)x_{E3}(t)}}\\ =-s_2^{-1}(\left[\begin{array}{cc} 0 & ms_2 \\ \end{array}\right]- \left[\begin{array}{cc} \Phi s_1 & \Phi s_2 \\ \end{array}\right])x_{E3}(t)\\ =-s_2^{-1}\left[\begin{array}{cc} -\Phi s_1 & ms_2-\Phi s_2 \\ \end{array}\right]x_{E3}(t)\\ =-\left[\begin{array}{cc} -\Phi m & m-\Phi \\ \end{array}\right] \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{S_E^{-1}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)}\\ =-\underbrace{(m-\Phi)}_{f}{\dot x}(t)+\underbrace{m(-\Phi)}_{f_I}(r-y(t))\\ \Rightarrow u_\ell(t)=-fx(t)+f_I\int_0^t (r-y(\tau))d\tau \end{array} }

\displaystyle{(126)\quad \begin{array}{l} p_2\Phi+\Phi^Tp_2=-1\Rightarrow p_2=-\frac{1}{2}\Phi^{-1} \end{array} }

\displaystyle{(127)\quad \begin{array}{l} {\dot u}_n(t) =-s_2^{-1}\rho\, {\rm sgn}(\underbrace{-\frac{1}{2}\Phi^{-1}}_{p_2} \underbrace{ \left[\begin{array}{cc} s_1 & s_2 \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{S_E^{-1}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)})\\ =\rho\, {\rm sgn}(\frac{1}{2}\Phi^{-1} \left[\begin{array}{cc} 1 & m \end{array}\right] \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right])\\ =\rho\, {\rm sgn}(-\underbrace{\frac{1}{2}(-\Phi^{-1})}_{g}{\dot x}(t)+\underbrace{\frac{1}{2}m(-\Phi^{-1})}_{g_I} (r-y(t)))\\ \Rightarrow u_n(t)=\rho\,\int_0^t\underbrace{{\rm sgn}(-g{\dot x}(\tau)+g_I(r-y(\tau)))}_{0,\pm 1}d\tau \end{array} }

m=1, \Phi=-0.5の場合のシミュレーション結果を次に示します。

[5] 数値例(2)

\displaystyle{(201)\quad \begin{array}{l} \dot{x}(t)= \underbrace{\left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]}_{A}x(t) +\underbrace{\left[\begin{array}{cc} 0 \\ 1 \end{array}\right]}_{B}u(t)+ \left[\begin{array}{cc} 0 \\ w \end{array}\right]\\ y(t)=\underbrace{\left[\begin{array}{cc} 1 & 0 \end{array}\right]}_{C}x(t)\\ \end{array} }

\displaystyle{(209)\quad \left[\begin{array}{c} 0 \\ -w \\ r \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{array}\right] }_{S_E} \left[\begin{array}{c} x_{1\infty} \\ x_{2\infty} \\ u_\infty \end{array}\right] \Rightarrow \left[\begin{array}{c} x_{1\infty} \\ x_{2\infty} \\ u_\infty \end{array}\right]= \left[\begin{array}{c} r \\ 0 \\ -w \end{array}\right] }

\displaystyle{(211)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{\dot{x}_{E3}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }_{x_{E3}(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] }_{B_{E3}} {\dot u}(t) }

\displaystyle{(216)\quad \begin{array}{l} s(t)= \underbrace{ \left[\begin{array}{ccc} s_{11} & s_{12} & s_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{s_2 \left[\begin{array}{ccc} m_1 & m_2 & 1 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)}\\ (m_1=s_2^{-1}s_{11},m_2=s_2^{-1}s_{12}) \end{array} }

\displaystyle{(217)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ s_{11} & s_{12} & s_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)}\Leftrightarrow\\ \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x_{E3}(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ -s_2^{-1}s_{11} & -s_2^{-1}s_{12} & s_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} \end{array} }

\displaystyle{(218)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{x}'_{E3}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -m_1 & -m_2 & s_2^{-1} \\ -s_{12}m_1 & -s_{12}m_2 & m_2 \end{array}\right] }_{T_sA_{E3}T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{x'_{E3}(t)} + \underbrace{ \left[\begin{array}{cc} 0\\ 0\\ s_2 \end{array}\right] }_{T_sB_{E3}} {\dot u}(t) }

\displaystyle{(223)\quad \begin{array}{l} {\dot u}_\ell(t)= -\underbrace{s_2^{-1}}_{(SB_{E3})^{-1}} \underbrace{(\left[\begin{array}{ccc} -s_{12}m_1 & -s_{12}m_2 & m_2 \\ \end{array}\right]x'_{E3}(t)-\Phi s(t))}_{(SA_{E3}-\Phi S)x_{E3}(t)}}\\ =-s_2^{-1}(\left[\begin{array}{ccc} 0 & 0 & m_2s_2 \\ \end{array}\right]- \left[\begin{array}{ccc} \Phi s_{11} & \Phi s_{12} & \Phi s_2 \\ \end{array}\right])x_{E3}(t)\\ =-s_2^{-1}\left[\begin{array}{ccc} -\Phi s_{11} & -\Phi s_{12} & m_2s_2-\Phi s_2 \\ \end{array}\right]x_{E3}(t)\\ =-\left[\begin{array}{ccc} -\Phi m_1 & -\Phi m_2 & m_2-\Phi \\ \end{array}\right] \underbrace{ \left[\begin{array}{ccc} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right] }_{S_E^{-1}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)}\\ =-\underbrace{ \left[\begin{array}{cc} -\Phi m_2 & m_2-\Phi \\ \end{array}\right] }_{F}{\dot x}(t)+\underbrace{m_1(-\Phi)}_{F_I}(r-y(t))\\ \Rightarrow u_\ell(t)=-Fx(t)+F_I\int_0^t (r-y(\tau))d\tau \end{array} }

\displaystyle{(226)\quad \begin{array}{l} p_2\Phi+\Phi^Tp_2=-1\Rightarrow p_2=-\frac{1}{2}\Phi^{-1} \end{array} }

\displaystyle{(237)\quad \begin{array}{l} {\dot u}_n(t) =-s_2^{-1}\rho\, {\rm sgn}(\underbrace{-\frac{1}{2}\Phi^{-1}}_{p_2} \underbrace{ \left[\begin{array}{ccc} s_{11} & s_{12} & s_2 \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{ccc} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right] }_{S_E^{-1}} \underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right] }_{x_{E2}(t)})\\ =s_2^{-1}\rho\, {\rm sgn}(\frac{1}{2}\Phi^{-1} \left[\begin{array}{ccc} s_{12} & s_2 & s_{11} \end{array}\right] \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right])\\ =\rho\, {\rm sgn}(\frac{1}{2}\Phi^{-1} \left[\begin{array}{ccc} m_2 & 1 & m_1 \end{array}\right] \left[\begin{array}{c} {\dot x}(t) \\ y(t)-r \end{array}\right])\\ =\rho\, {\rm sgn}(-\underbrace{\frac{1}{2}(-\Phi^{-1}) \left[\begin{array}{cc} m_2 & 1 \end{array}\right] }_{G}{\dot x}(t)+\underbrace{\frac{1}{2}m_1(-\Phi^{-1})}_{G_I} (r-y(t)))\\ \Rightarrow u_n(t)=\rho\,\int_0^t\underbrace{{\rm sgn}(-G{\dot x}(\tau)+G_I(r-y(\tau)))}_{0,\pm 1}d\tau \end{array} }

m_1=m_2=1, \Phi=-0.5の場合のシミュレーション結果を次に示します。

演習…Flipped Classroom

MATLAB 
%cCIP_smci2.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0 %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
%-----
 Tcart=1; Mr=0.5;
 Tpend=1/4*2*pi*sqrt(4/3*ell/g); Mth=3/180*pi;
 Tamp=0.01; Mamp=10; 
 CM=eye(2,4); DM=zeros(2,1);
 CS=[1 0]; C=CS*CM; D=CS*DM;
 n=size(A,1); m=size(C,1); r=size(B,2);
%-----
 AE=[A B;zeros(m,n+m)];
 BE=[zeros(n,m);eye(m)];
 SE=[A B;C D];
 CE=eye(n+m);
 QE=diag([1/Mr 1/Mth Tcart/Mr Tpend/Mth 1/Mamp].^2);
 S=swflqr(AE,BE,QE)
%-----
 Phi=-0.5;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S/SE;
 G=P2S(1:4)
 GI=P2S(5) 
 Leq=inv(S*BE)*S*AE;
 LPhi=-inv(S*BE)*Phi*S;
 L=(Leq+LPhi)/SE;
 F=L(1:4)
 FI=L(5)
 rho=1;
 Ln=inv(S*BE)*rho
%-----
 x0=[0;th0*0;0;0];
 sim('CIP_smci2_2015a.mdl')
%-----
%eof
SCILAB

応用例

投稿日時: 投稿者:


応用例…Homework

[1]　ある航空機の線形状態方程式として、次を考えます。

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \theta & pitch\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 1 & 0 & 0 & 0\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }

A行列の固有値は次のように求められます。

\displaystyle{(2)\quad \lambda(A):\left\{\begin{array}{cl} -1.864\pm j 3.660 & short\ period\ mode\\ 0 & pitch\ attitude\ mode\\ -20 & elevator\ actuator\ mode\\ -20 & flap\ actuator\ mode \end{array}\right. }

状態変数として、pitch angle \thetaの代わりに、fligt path angle \gamma=\theta-\alphaを用いると、次式となります。

\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \gamma=\theta-\alpha & fligt\ path\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 0 & 1.74 & 0.08 & 0.59\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }

また、評価変数は、pitch angle \thetaとfligt path angle \gammaとすると、その出力方程式は次式のように与えられます。

\displaystyle{(4)\quad \begin{array}{l} y(t)=Hx(t)\\ H= \left[\begin{array}{ccccc} 1 & 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 \end{array}\right] \end{array} }

制御目的は、次図のように、pitch angle \thetaとfligt path angle \gammaの間の非干渉化を達成し、独立して設定値に合せることです。


図1 航空機の操縦法

[3] 積分動作導入による追従制御

この場合の制御則は次のように構成されます。

\displaystyle{(20)\quad { \boxed{\begin{array}{l} u(t)=u_\ell(t)+u_n(t)\\ u_\ell(t)=-\underbrace{(SB_E)^{-1}(SA_E-\Phi S)}_{L=L_{eq}+L_\Phi}\left[\begin{array}{c} x_r(t)\\ x(t) \end{array}\right]\\ -\underbrace{(SB_E)^{-1}(\Phi S_r+S_1B_r)}_{L_r} r(t) +\underbrace{(SB_E)^{-1}S_r}_{L_{\dot r}} \dot{r}(t)\\ u_n(t)=-\underbrace{(SB_E)^{-1}\rho(t,x)}_{L_n}\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||} \end{array}}} }

●ここでは、等価制御による閉ループ系のA_{eq}が望ましい固有値・固有ベクトルをもつように設計して、次式からスイッチング関数を決める行列Sを決定します。

\displaystyle{(21)\quad { S[v_1\cdots v_{n-m}]=0 %　\Rightarrow [v_1\cdots v_{n-m}]^TS^T=0} }

A_{eq}の望ましい固有値を

\displaystyle{(22)\quad \lambda(A_{eq}):\left\{\begin{array}{cl} -5.6\pm j 4.2 & short\ period\ mode\\ -1 & pitch\ attitude\ mode\\ -0.4 & elevator\ deflection\ mode\\ -0.7 & flap\ deflextion\ mode \end{array}\right. }

とし、また望ましいモード分布を次により指定します。

\displaystyle{(23)\quad {\tt specpos}= \left[\begin{array}{ccccc} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\\hline 1 & 1 & 1 & 1 & 1\\ 1 & 0 & 1 & 1 & 1\\ 0 & 1 & 1 & 0 & 1\\\hline 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ \end{array}\right],\quad {\tt specent}= \left[\begin{array}{rrrrr} -1 & -1 & -1 & -1 & -1\\ -1 & -1 & -1 & -1 & -1\\\hline 0 & 0 & 1 & 0 & 1\\ 1 & -1 & 0 & 1 & 0\\ -1 & 1 & 0 & -1 & 0\\\hline -1 & -1 & -1 & -1 & -1\\ -1 & -1 & -1 & -1 & -1\\ \end{array}\right] }

このとき次のような行列Sを得ます。

\displaystyle{(24)\quad S=\left[\begin{array}{rrrrr|rr} 0.0531 & 0.0137 & -0.1436 & -0.0260 & -0.1373 & 0.0500 & 0\\ -0.0072 & -0.0612 & 0.1635 & 0.0035 & 0.1661 & 0 & 0.0500 \end{array}\right] }

このときSMC則における行列Lなどは次のように計算されます。

\displaystyle{(25)\quad L=\left[\begin{array}{rrrrrrr} -1.0616 & -0.2743 & 2.9379 & 0.6060 & 2.4605 & -0.4927 & -0.0900\\ 0.1440 & 1.2236 & -3.3390 & -0.2296 & -3.2769 & 0.0671 & 0.0142 \end{array}\right] }

\displaystyle{(26)\quad L_{r}=\left[\begin{array}{rr} -2.9499 & 0.0120\\ 0.4000 & 2.9391 \end{array}\right] }

\displaystyle{(27)\quad L_{\dot{r}}=\left[\begin{array}{rr} -0.1448 & 0.0013\\ 0.0196 & 0.1439 \end{array}\right] }


図2 積分動作導入によるSMC制御系応答

MATLAB 
%cAITCRAFT_smci.m
%-----
 clear all, close all
 A0= [0 0 1.74 0.08 0.59;
     0 -1.99 -13.41 -18.95 -3.6;
     0 1 -1.74 -0.08 -0.59;
     0 0 0 -20 0;
     0 0 0 0 -20];
 B0=[0 0;0 0;0 0;20 0;0 20]; 
 H=[1 0 1 0 0;1 0 0 0 0];
 A=[zeros(2,2) -H ;
    zeros(5,2) A0];
 B=[zeros(2,2); B0];
 Gamma=diag([-0.9,-0.7]);
%----- 
 lambda=[-5.6 4.2 -1 -0.4 -0.7];
 nocomp=1;
 specpos=[zeros(2,5);[1 1 1 1 1;1 0 1 1 1;0 1 1 0 1];zeros(2,5)]
 specent=[-ones(2,5);[0 0 1 0 1;1 -1 0 1 0;-1 1 0 -1 0];-ones(2,5)]
 S=swfvpl(A,B,lambda,nocomp,specpos,specent);
 SS=S(:,6:7)\S*0.05
 F=(SS*B)\(SS*A)
 eig(A-B*F)
%-----
 [nn,mm]=size(B);
 S1=SS(:,1:nn-mm);
 S2=SS(:,nn-mm+1:nn);
 Br=eye(nn-mm,mm);
 B2=B(nn-mm+1:nn,:); 
 Lambda=S2*B2;
 P2=diag([0.025,0.025]);
 Phi=-20/0.025*P2;
 Check=P2*Phi+Phi*P2
 L=inv(Lambda)*(SS*A-Phi*SS)
%-----
%Sr=[-0.1448 0.0013;0.0196 0.1439]
 Sr=zeros(2)
 Lr=inv(Lambda)*(Phi*Sr+S1*Br)
 Ldr=inv(Lambda)*Sr
 Ln=0.1*inv(Lambda)
%-----
 sim('AIRCRAFT_smci_2015a.mdl')
%-----
%eof
SCILAB

追従SMI制御

投稿日時: 投稿者:

積分動作導入による追従SM制御…Homework

[1] 制御対象

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ y(t)=Cx(t)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, y(t)\in{\rm\bf R}^p, m=p) \end{array} }

の出力を、コマンド（次式の解）

\displaystyle{(2)\quad \begin{array}{l} %\dot{r}(t)=\Gamma(r(t)-r_c)\\ \frac{d}{dt}(r(t)-r_c)=\Gamma(r(t)-r_c)\\ (r(t),r_c\in{\rm\bf R}^m) \end{array} }

に追従させることを考えます（\Gammaは安定行列）。そのために、積分動作

\displaystyle{(3)\quad \begin{array}{l} \dot{x}_r(t)=r(t)-y(t)\\ (x_r(t)\in{\rm\bf R}^m) \end{array} }

を導入し、次の拡大系を構成します。ここで、(1)はすでに標準形であるとしています。

\displaystyle{(4)\quad \begin{array}{l} \left[\begin{array}{c} \dot x_r(t)\\ \dot x(t) \end{array}\right] = \left[\begin{array}{c|cc} 0 & -C_1 & -C_2\\\hline 0 & A_{11} & A_{12} \\ 0 & A_{21} & A_{22} \end{array}\right] \left[\begin{array}{c} x_r(t)\\ x(t) \end{array}\right] + \left[\begin{array}{c} 0\\\hline 0\\ B_2 \end{array}\right] u(t) + \left[\begin{array}{c} I_m \\\hline 0 \\ 0 \end{array}\right] r(t)\\ (x_r(t)\in{\rm\bf R}^m, x(t)\in{\rm\bf R}^n) \end{array} }

これを、次のように分割し直しても標準形であることには変わりありません。

\displaystyle{(5a)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc|c} 0 & -C_1 & -C_2\\ 0 & A_{11} & A_{12} \\\hline 0 & A_{21} & A_{22} \end{array}\right] }_{A_E} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0\\ 0\\\hline B_2 \end{array}\right] }_{B_E} u(t) + \left[\begin{array}{c} I_m \\ 0 \\\hline 0 \end{array}\right] r(t)\\ (x_1(t)\in{\rm\bf R}^n, x_2(t)\in{\rm\bf R}^m) \end{array} }

ただし

\displaystyle{(5b)\quad %\begin{array}{l} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] = \left[\begin{array}{c} x_r(t)\\ x(t) \end{array}\right] %\end{array} }

●この積分器による拡大系を安定化できれば、積分器の値x_r(t)は定値となり、被積分項r(t)-y(t)の値は零となり、y(t)はr(t)へ漸近します。そこで、SM制御によって拡大系を安定化し、追従制御系を構成することを考えます。この制御系は特別なr(t)=0の場合を含みますので、まずスイッチング関数として、次式を考えます。

\displaystyle{(6)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} %\underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] %}_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I_m \\ \end{array}\right] }_{S} %\underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] %}_{x(t)} \ (M=S_2^{-1}S_1) }

(5)に対して、座標変換

\displaystyle{(7)\quad \begin{array}{l} \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} I_n & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]\\ \Leftrightarrow \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} I_n & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] \end{array} }

を行って、次式を得ます。

\displaystyle{(8a)\quad \begin{array}{l} %\underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] %}_{\dot{x}'(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right] }_{T_sA_ET_s^{-1}} %\underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] %}_{x'(t)} + \underbrace{ \left[\begin{array}{cc} 0\\ S_2B_2 \end{array}\right] }_{T_sB_E} u(t)\\ + \left[\begin{array}{cc} B_r \\ S_1B_r \end{array}\right] r(t) \end{array} }

ただし

\displaystyle{(8b)\quad \left\{\begin{array}{l} \bar{A}_{11}= \left[\begin{array}{cc} 0 & -C_1 \\ 0 & A_{11} \end{array}\right] -\left[\begin{array}{c} -C_2\\ A_{12} \end{array}\right]M\\ \bar{A}_{12}= \left[\begin{array}{c} -C_2\\ A_{12} \end{array}\right]S_2^{-1}\\ \bar{A}_{21}=S_2(M\bar{A}_{11} + \left[\begin{array}{cc} 0 & A_{21} \end{array}\right] -A_{22}M)\\ \bar{A}_{22}=S_2(M \left[\begin{array}{c} -C_2\\ A_{12} \end{array}\right] +A_{22})S_2^{-1}\\ B_r=\left[\begin{array}{cc} I_m \\ 0 \end{array}\right] \end{array}\right. }

ここで、\bar{A}_{11}が安定行列となるように行列Mが選ばれているとします。

●特別なr(t)=0の場合のスライディングモードはs(t)=0で表されますが、一般のr(t)\ne0の場合のスライディングモードは

\displaystyle{(9)\quad {\boxed{s(t)-S_rr(t)=0} \quad\Rightarrow \dot{s}(t)-S_r\dot{r}(t)=0 } }

で表されるとします。ここで、S_rは適当に選択された行列です。

●以上の準備の下で、SM制御則

\displaystyle{(10)\quad u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component} }

を、２次安定性

\displaystyle{(11)\quad \boxed{ \begin{array}{lll} V(\bar{x})= \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t)-S_rr(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} P_1 & 0\\ 0 & P_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t)-S_rr(t) \end{array}\right] }_{\bar{x}(t)}\\ \Rightarrow \dot{V}(\bar{x})\le - \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t)-S_rr(t) \end{array}\right]^T }_{\bar{x}^T(t)} \underbrace{ \left[\begin{array}{cc} Q_1 & 0\\ 0 & I \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t)-S_rr(t) \end{array}\right] }_{\bar{x}(t)} \end{array}}} }

すなわち

\displaystyle{(11')\quad \begin{array}{lll} V(\bar{x})= \underbrace{x_1^T(t)P_1x_1(t)}_{V(x_1)}+\underbrace{(s(t)-S_rr(t))^TP_2(s(t)-S_rr(t))}_{V(s)}\\ \Rightarrow \dot{V}(\bar{x})\le -x_1^T(t)Q_1x_1(t)-(s(t)-S_rr(t))^T(s(t)-S_rr(t)) \end{array}} }

が成り立つように決定します（P_1>0, P_2>0, Q_1>0）。

[2] 可到達性の検討

ここでは、スライディングモード制御則(10)の具体的な表現を求めます。

●等価制御は、(8)においてr(t)=0の場合

\displaystyle{ \begin{array}{cl} (12.1) & s(t)=0\\ \Downarrow &\\ (12.2) & \dot{s}(t)=0\\ \Downarrow &\\ (12.3) & 0=S_2\bar{A}_{21}x_1(t)+S_2\bar{A}_{22}S_2^{-1}s(t)+S_2B_2u(t)\\ \Downarrow &\\ (12.4) & u_{eq}(t)=-\underbrace{(S_2B_2)^{-1}}_{([0\ I]T_sB_E)^{-1}} \underbrace{\left[\begin{array}{cc} S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right]}_{[0\ I]T_sA_ET_s^{-1}} \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ \Downarrow & by\ S=[0\ I]T_s\\ (12.5) & u_{eq}(t)=-(SB_E)^{-1}SA_ET_s^{-1} \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ \Downarrow & by\ (7),(8b)\\ (12.6) & u_{eq}(t)=-(SB_E)^{-1}SA_E \left[\begin{array}{c} x_r(t)\\ x(t) \end{array}\right]} \end{array} }

のように得られます。(10)の第１項u_\ellは、この等価制御をベースして、(9)を考慮して

\displaystyle{(13)\quad \begin{array}{l} u_\ell(t)=-(S_2B_2)^{-1} (\left[\begin{array}{cc} S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right]\left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]-\Phi s(t))\\ -(S_2B_2)^{-1}(\Phi S_r+S_1B_r)r(t)+(S_2B_2)^{-1}S_r\dot{r}(t)\\ \end{array} }

すなわち

\displaystyle{(13')\quad { \boxed{\begin{array}{l} u_\ell(t) =-\underbrace{(SB_E)^{-1}(SA_E-\Phi S)}_{L=L_{eq}+L_\Phi}\left[\begin{array}{c} x_r(t)\\ x(t) \end{array}\right]\\ -\underbrace{(SB_E)^{-1}(\Phi S_r+S_1B_r)}_{L_r} r(t) +\underbrace{(SB_E)^{-1}S_r}_{L_{\dot r}} \dot{r}(t)\\ \end{array}}} }

のように構成します（\Phiは安定行列）。このとき閉ループ系は次式で与えられます。

\displaystyle{(14)\quad \begin{array}{l} \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ + \left[\begin{array}{cc} 0\\ S_2B_2 \end{array}\right] (-(S_2B_2)^{-1} (\left[\begin{array}{cc} S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right]\left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]-\Phi s(t))\\ -(S_2B_2)^{-1}(\Phi S_r+S_1B_r)r(t)+(S_2B_2)^{-1}S_r\dot{r}(t)+u_n(t)) + \left[\begin{array}{cc} B_r \\ S_1B_r \end{array}\right] r(t)\\ = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ S_2\bar{A}_{21} & S_2\bar{A}_{22}S_2^{-1} \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ + \left[\begin{array}{cc} 0\\ I_m \end{array}\right] \left[\begin{array}{cc} -S_2\bar{A}_{21} & -S_2\bar{A}_{22}S_2^{-1}+\Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]\\ + \left[\begin{array}{cc} 0\\ -(\Phi S_r+S_1B_r)r(t)+S_r\dot{r}(t)+S_2B_2u_n(t) \end{array}\right] + \left[\begin{array}{cc} B_r \\ S_1B_r \end{array}\right] r(t)\\ = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ 0 & \Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right]+ \left[\begin{array}{cc} B_rr(t)\\ -\Phi S_rr(t)+S_r\dot{r}(t)+S_2B_2u_n(t) \end{array}\right] \end{array} }

これを変形して

\displaystyle{(15)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{s}(t)-S_r\dot{r}(t) \end{array}\right] = \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ 0 & \Phi \\ \end{array}\right] \left[\begin{array}{c} x_1(t)\\ s(t)-S_rr(t) \end{array}\right]\\ + \left[\begin{array}{cc} B_rr(t)\\ -\Phi S_rr(t)+S_r\dot{r}(t)+S_2B_2u_n(t) \end{array}\right] \end{array} }

すなわち、次式が成り立ちます。

\displaystyle{ \begin{array}{cl} (15.1) & \dot{x}_1(t)=\bar{A}_{11}x_1(t)+\bar{A}_{12}(s(t)-S_rr(t))+\underbrace{(B_r+\bar{A}_{12}S_r)}_{B_r'}r(t)\\ (15.2) & \dot{s}(t)-S_r\dot{r}(t)=\Phi(s(t)-S_rr(t))+S_2B_2u_n(t) \end{array} }

(15.2)に基づいてスライディングモードs-S_rr=0を達成するu_nを明らかにします。

\Phiは安定行列なので

\displaystyle{(16)\quad \begin{array}{l} P_2\Phi+\Phi^TP_2=-I_m \end{array} }

を満たすP_2>0を選ぶことができます。これを用いて

\displaystyle{(17)\quad \begin{array}{l} u_n(t)=-(S_2B_2)^{-1}\rho(u_\ell,y)\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||} \end{array} }

すなわち

\displaystyle{(17')\quad { \boxed{ \begin{array}{l} u_n(t) =-\underbrace{(SB_E)^{-1}\rho(u_\ell,y)}_{L_n}\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||} \end{array} }} }

と選びます（\rho>0）。このとき次式が成り立ち、(11)が示されます。

\displaystyle{(18)\quad \begin{array}{l} \dot{V}(s)=2(s(t)-S_rr(t))^TP_2(\dot{s}(t)-S_r\dot{r}(t))\\ =2(s(t)-S_rr(t))^TP_2 (\Phi (s(t)-S_rr(t))-\rho(u_\ell,y)\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||})\\ =(s(t)-S_rr(t))^T\underbrace{(P_2\Phi+\Phi^TP_2)}_{-I_m}(s(t)-S_rr(t))\\ +2(s(t)-S_rr(t))^TP_2(-\rho(u_\ell,y)\frac{P_2(s(t)-S_rr(t))}{||P_2(s(t)-S_rr(t))||})\\ \le -||s(t)-S_rr(t)||^2-2\rho(u_\ell,y)||P_2(s(t)-S_rr(t))||\\ \le -||s(t)-S_rr(t)||^2 \end{array} }

これより(15.2)における２次安定性が成り立つことが分かります。２次安定であればs-S_rrは0に漸近します。

[3] スライディングモードの検討

次にs-S_rr=0のもとで、(15.1)からx_1の振舞いについて調べます。

\bar{A}_{11}は安定行列なので

\displaystyle{(19)\quad \begin{array}{l} P_1\bar{A}_{11}+\bar{A}_{11}^TP_1=-Q_1<0 \end{array} }

を満たすP_1>0を選ぶことができます。そこで、(15.1)に対するリャプノフ関数V(x_1)=x_1^T(t)P_1x_1(t)の微分について調べます。

\displaystyle{(20)\quad \begin{array}{l} \dot{V}(x_1)=2x_1^T(t)P_1\dot{x}_1(t)\\ =2x_1^T(t)P_1( \bar{A}_{11}x_1(t)+\bar{A}_{12}\underbrace{(s(t)-S_rr(t))}_{0}+B_r'r(t))\\ =x_1^T(t)\underbrace{(P_1\bar{A}_{11}+\bar{A}_{11}^TP_1)}_{-Q_1}x_1(t)+2x_1^T(t)P_1B_r'r(t)\\ =-x_1^T(t)Q_1x_1(t)+2x_1^T(t)P_1B_r'r(t)\\ \le -\sigma_n(Q_1)x_1^T(t)x_1(t)+2||x_1^T(t)P_1||\,||B_r'r(t)||\\ \le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{x_1^T(t)P_1^2x_1(t)}\,||B_r'r(t)||\\ \le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{\sigma_1(P_1^2)x_1^T(t)x_1(t)}\,||B_r'r(t)||\\ = -\sigma_n(Q_1)||x_1(t)||^2+2\sigma_1(P_1)||x_1(t)||\,||B_r'r(t)||\\ = -\sigma_1(P_1)||x_1(t)||(\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)||-2||B_r'r(t)||) \end{array} }

したがって、次式が成り立ち、リャプノフ安定性が示されます。

\displaystyle{(21)\quad||B_r'r(t)||<\frac{1}{2}\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)|| \Rightarrow \dot{V}(x_1)<0 }

これから、ロバスト性での議論と同様に、x_1が発散することはないことが分かります。

演習C61…Flipped Classroom

MATLAB 
%cCIP_smci.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 A0=A; B0=B; 
 H=[1 0 0 0];
 A=[0 -H ;
    zeros(4,1) A0];
 B=[0; B0];
 CM=[1 0 0 0;
     0 1 0 0];
 Gamma=-5; 
%-----
 lambda=[-0.8 0.5 -1.5 1];
 nocomp=2;
 specpos=rand(4,4);
 specent=rand(4,4);
 S=swfvpl(A,B,lambda,nocomp,specpos,specent)      
 SS=S;
 F=(SS*B)\(SS*A)
 eig(A-B*F)
%-----
[nn,mm]=size(B);
 S1=SS(:,1:nn-mm);
 S2=SS(:,nn-mm+1:nn);
 Br=eye(nn-mm,mm);
 B2=B(nn-mm+1:nn,:); 
 Lambda=S2*B2;
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=2;
 Ln=inv(S*B)*rho;
%-----
 Sr=0;
 Lr=inv(Lambda)*(Phi*Sr+S1*Br)
 Ldr=inv(Lambda)*Sr
%Ln=0.1*inv(Lambda)
%-----
 x0=[0;th0*0;0;0];
 sim('CIP_smci_2015a.mdl')
%-----
%eof
SCILAB

応用例

投稿日時: 投稿者:


応用例…Homework

[1]　ある航空機の線形状態方程式として、次を考えます。

\displaystyle{(1)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \theta & pitch\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 1 & 0 & 0 & 0\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }

A行列の固有値は次のように求められます。

\displaystyle{(2)\quad \lambda(A):\left\{\begin{array}{cl} -1.864\pm j 3.660 & short\ period\ mode\\ 0 & pitch\ attitude\ mode\\ -20 & elevator\ actuator\ mode\\ -20 & flap\ actuator\ mode \end{array}\right. }

状態変数として、pitch angle \thetaの代わりに、fligt path angle \gamma=\theta-\alphaを用いると、次式となります。

\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)\\ x(t)= \left[\begin{array}{ll} \gamma=\theta-\alpha & fligt\ path\ angle\\ q & pitch\ rate\\ \alpha & angle\ of\ atack\\ \eta & elevator\ deflection\\ \delta & flap\ deflection \end{array}\right],\quad u(t)= \left[\begin{array}{ll} \eta_c & elevator\ command\\ \delta_c& flap\ comannd \end{array}\right]\\ A=\left[\begin{array}{rrrrr} 0 & 0 & 1.74 & 0.08 & 0.59\\ 0 &-1.99 & -13.41 & -18.95 & -3.60\\ 0 & 1 & -1.74 & -0.08 & -0.59\\ 0 & 0 & 0 & -20 & 0\\ 0 & 0 & 0 & 0 & -20 \end{array}\right],\quad B=\left[\begin{array}{rr} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 20 & 0\\ 0 & 20 \end{array}\right] \end{array} }

また、評価変数として、pitch angle \thetaとfligt path angle \gammaをとると、その出力方程式は次式のように与えられます。

\displaystyle{(4)\quad \begin{array}{l} y(t)=Hx(t)\\ H= \left[\begin{array}{ccccc} 1 & 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 \end{array}\right] \end{array} }

制御目的は、次図のように、pitch angle \thetaとfligt path angle \gammaの間の非干渉化を達成し、独立して設定値に合せることです。


図1 航空機の操縦法

[2] モデル規範による追従制御

この場合の制御則は次のように構成されます。

\displaystyle{(5)\quad { \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }

●追従すべき規範モデルは

\displaystyle{(6)\quad \begin{array}{l} \dot{x}(t)=A_mx(t)+B_mu(t)\\ y(t)=Hx(t) \end{array}} }

ただし

\displaystyle{(7)\quad \begin{array}{l} A_m=A-BF\\ B_m=BG\\ G=-(H(A-BF)^{-1}B)^{-1} \end{array} }

いまこの規範モデルにもたせるべき望ましい固有値を

\displaystyle{(8)\quad \lambda(A_m):\left\{\begin{array}{cl} -5.6\pm j 4.2 & short\ period\ mode\\ -1 & pitch\ attitude\ mode\\ -20 & elevator\ actuator\ mode\\ -20 & flap\ actuator\ mode \end{array}\right. }

とし、また望ましいモード分布を次により指定します。

\displaystyle{(9)\quad {\tt specpos}= \left[\begin{array}{ccc|cc} 1 & 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 0 & 0\\ 0 & 1 & 1 & 0 & 0\\\hline 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{array}\right],\quad {\tt specent}= \left[\begin{array}{rrr|rr} 0 & 0 & 1 & -1 & -1\\ 1 & -1 & 0 & -1 & -1\\ -1 & 1 & 0 & -1 & -1\\\hline -1 & -1 & -1 & 1 & -1\\ -1 & -1 & -1 & -1 & 1\\ \end{array}\right] }

ここで、{\tt specpos}は望ましい固有ベクトルの要素を指定する場合１、指定しない場合０で表した行列、{\tt specエent}は望ましい固有ベクトルの要素が非零指定の場合１、零の場合０、非零任意の場合－１で表した行列です。

このとき以下に示すプログラムを使って、次が得られます。

\displaystyle{(10)\quad \begin{array}{l} F=WV^{-1} =\left[\begin{array}{rrrrr} -1.9338 & -0.5744 & -2.1188 & 0.4750 & 0.1066\\ 1.9571 & 0.2253 & 3.1273 & -0.0694 & -0.0515 \end{array}\right]\\ W=\left[\begin{array}{rrrrr} 0.8445 & 1.9935 & -0.2742 & 0 & 0\\ -2.7055 & 1.2780 & 0.1767 & 0 & 0 \end{array}\right]\\ V=\left[\begin{array}{rrrrr} 0 & -0.0000 & 0.6667 & 0.0030 & -0.0313\\ 1.0000 & -9.5000 & -0.3333 & 0.9964 & 0.1493\\ -0.9286 & 1.0000 & -0.3333 & -0.0528 & 0.0238\\ -1.8252 & -2.2364 & 0.2886 & 1.0000 & -0.0649\\ 2.9860 & -2.6459 & -0.1860 & -0.0825 & 1.0000 \end{array}\right] \end{array} }

●ここでは、スイッチング関数

\displaystyle{(11)\quad s(t)= \underbrace{ \left[\begin{array}{cc} M & I \end{array}\right] }_{S} x(t) }

を構成するために

\displaystyle{(12)\quad \begin{array}{l} \dot{x}_1(t)=A_{11}x_1(t)+A_{12}\hat{u}(t) \end{array} }

を２次形式評価関数

\displaystyle{(13)\quad J=\int_0^\infty(x_1^T(t)Q_{11}x_1(t)+\hat{u}^T(t)Q_{22}\hat{u}(t))dt }

を最小にする

\displaystyle{(14)\quad \hat{u}(t)=-Mx_1(t) }

を求めます。

\displaystyle{(15)\quad Q_{11}={\rm diag}\{10,5,5\},\ Q_{22}={\rm diag}\{20,20\} }

と選ぶとき、スイッチング関数を決める行列Sが、次のように得られます。

\displaystyle{(16)\quad S=\left[\begin{array}{rrr|rr} 0.7025 & 0.4209 & 0.1894 & -1 & 0\\ -0.0810 & 0.0635 & 0.0267 & 0 & -1\\ \end{array}\right] }

●P_2\Phi+\Phi^TP_2=-Iを満たす\PhiとP_2は、たとえば次のように定めます。

\displaystyle{(17)\quad \underbrace{ \left[\begin{array}{rr} 0.025 & 0\\ 0 & 0.025 \end{array}\right] }_{P_2} \underbrace{ \left[\begin{array}{rr} 20 & 0\\ 0 & 20 \end{array}\right] }_{\Phi} +(*)^T=-I }

このときSMC則における行列Lは次のように計算されます。

\displaystyle{(18)\quad L=\left[\begin{array}{rrrrr} -1.2285 & -0.1858 & -2.1624 & 0.0782 & 0.0460\\ 1.8631 & 0.2827 & 3.0804 & -0.1299 & -0.0060 \end{array}\right] }

図2 モデル規範法によるSMC制御系応答

\displaystyle{{ \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }

MATLAB 
%cAIRCRAFT_smcm.m
%-----
 clear all, close all
 A= [0 0 1.74 0.08 0.59;
     0 -1.99 -13.41 -18.95 -3.6;
     0 1 -1.74 -0.08 -0.59;
     0 0 0 -20 0;
     0 0 0 0 -20];
 B=[0 0;0 0;0 0;20 0;0 20];
%----- 
 lambda=[-5.6 4.2 -1 -20 -20];
 nocomp=1;
%固有ベクトルの非零要素の位置 (１：指定、０：指定せず)
 specpos=[ 1  1  1  0  0; 
           1  0  1  0  0; 
           0  1  1  0  0;
           0  0  0  1  0;
           0  0  0  0  1]
%固有ベクトルの要素の値 (０：零、１：非零指定、－１：非零任意)       
 specent=[ 0  0  1 -1 -1;
           1 -1  0 -1 -1;
          -1  1  0 -1 -1;
          -1 -1 -1  1 -1;
          -1 -1 -1 -1  1]
 F=vplace2(A,B,lambda,nocomp,specpos,specent)       
 pl=eig(A-B*F)
%-----
 H=[1 0 1 0 0;1 0 0 0 0];
 Am=A-B*F;
 G=-inv(H*(Am\B));
 Bm=B*G;
%-----
 Q=diag([10,5,5,20,20]);
 S=swflqr(A,B,Q)
%-----
 P2=diag([0.025,0.025]);
 Phi=-20/0.025*P2;
 Check=P2*Phi+Phi*P2
 L=inv(S*B)*(S*Am-Phi*S)
 Ln=inv(S*B)
%-----
 sim('AIRCRAFT_smcm_2015a.mdl')
%-----
%eof

モデル規範による追従SM制御

投稿日時: 投稿者:

モデル規範による追従SM制御…Homework

[1]　次の状態方程式で表される制御対象を考えます。

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t) }

この状態が理想的なモデル

\displaystyle{(2)\quad \dot{x}_m(t)=A_mx_m(t)+B_mr(t) }

の状態を追従するように、すなわち

\displaystyle{(3)\quad e(t)=x(t)-x_m(t)\rightarrow 0 \quad(t\rightarrow\infty) }

となるように制御則を決定したいとします。

●追従誤差のダイナミックスは

\displaystyle{(4)\quad \begin{array}{l} %\bullet \dot{e}(t)=Ax(t)+Bu(t)-(A_mx_m(t)+B_mr(t))\\ %=Ax(t)-Ax_m(t)+Ax_m(t)-A_mx_m(t)+Bu(t)-B_mr(t)\\ %=Ae(t)+(A-A_m)x_m(t)+Bu(t)-B_mr(t)\\ %\bullet \dot{e}(t)=Ax(t)+Bu(t)-(A_mx_m(t)+B_mr(t))\\ =Ax(t)-A_mx(t)+A_mx(t)-A_mx_m(t)+Bu(t)-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t) \end{array} }

となります。いま条件

\displaystyle{(5)\quad \begin{array}{l} (BB^\dag-I)(A-A_m)=0\Leftrightarrow B\underbrace{B^\dag(A-A_m)}_{F}=A-A_m\\ (BB^\dag-I)B_m=0\Leftrightarrow B\underbrace{B^\dag B_m}_{G}=B_m \end{array} }

を仮定すると（B^\dagはBB^\dag B=Bを満たす疑似逆行列）、制御則

\displaystyle{(6)\quad u(t)=-Ke(t)-Fx(t)+Gr(t) }

の下では、誤差方程式は次式となります。

\displaystyle{(7)\quad \dot{e}(t)=(A_m-BK)e(t) }

実際、

\displaystyle{(8)\quad \begin{array}{l} \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)+B(-Ke(t)-Fx(t)+Gr(t))-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)-BKe(t)-BB^\dag(A-A_m)x(t)\\ +BB^\dag B_mr(t)-B_mr(t)\\ =(A_m-BK)e(t)-\underbrace{(BB^\dag-I)(A-A_m)}_{0}x(t)+\underbrace{(BB^\dag -I)B_m}_{0}r(t)\\ \end{array} }

●以上から、制御則(6)の設計の指針として、望ましいモデル(2)を(1)に対する状態フィードバックと入力変換で得るものとします。すなわち

\displaystyle{(9)\quad \boxed{\begin{array}{l} A_m=A-BF\\ B_m=BG \end{array}} }

ここでGは、rから出力

\displaystyle{(10)\quad y(t)=Hx(t) }

までの定常ゲイン

\displaystyle{(11)\quad %\begin{array}{l} \left\{\begin{array}{l} 0=(A-BF)x(\infty)+BGr(\infty)\\ y(\infty)=Hx(\infty) \end{array}\right. \Rightarrow y(\infty)=\underbrace{-H(A-BF)^{-1}BG}_{I}r(\infty) %\end{array} }

が単位行列となるように、次式で決めます。

\displaystyle{(12)\quad \boxed{G=-(H(A-BF)^{-1}B)^{-1}} }

ちなみにKについては、(7)が漸近安定であることが前提となり、望ましいモデル(2)より、追従誤差は速いダイナミックスをもつことが必要であることに注意します。

[2]　誤差方程式(4)すなわち

\displaystyle{(13)\quad \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t) }

に対して、スライディングモード制御を適用することを考えます。スイッチング関数を

\displaystyle{(14)\quad s(t)=Se(t) }

とするとき、スライディングモード時は

\displaystyle{(15)\quad Se(t)=0\Rightarrow S\dot{e}(t)=0 }

すなわち

\displaystyle{(16)\quad S\dot{e}(t)=S(A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t))=0 }

が成り立ち、等価制御は次式となります。

\displaystyle{(17)\quad u_{eq}(t)=-(SB)^{-1}S(A_me(t)+(A-A_m)x(t)-B_mr(t)) }

このとき追従誤差のダイナミックスは次式となります。

\displaystyle{(18)\quad \begin{array}{l} \dot{e}(t)=A_me(t)+(A-A_m)x(t)+Bu(t)-B_mr(t)\\ =A_me(t)+(A-A_m)x(t)\\ +B(-(SB)^{-1}S(A_me(t)+(A-A_m)x(t)-B_mr(t)))-B_mr(t)\\ =(I-B(SB)^{-1}S)(A_me(t)+(A-A_m)x(t)-B_mr(t))\\ =(I-B(SB)^{-1}S)(A_me(t)+BFx(t)-BGr(t))\\ =(I-B(SB)^{-1}S)A_me(t)+\underbrace{(I-B(SB)^{-1}S)B}_{0}(Fx(t)-Gr(t))\\ =\underbrace{(I-B(SB)^{-1}S)A_m}_{A_{eq}}e(t) \end{array} }

したがって、スライディングモード制御を適用する場合、全体の制御則は、次式で表されます。

\displaystyle{(19)\quad { \boxed{\begin{array}{l} u(t)=-Fx(t)+Gr(t)+\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA_m-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ \displaystyle{u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I)} \end{array}}} }

演習C51…Flipped Classroom

MATLAB 
%cCIP_smcm.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 H=[1 0 0 0];
%-----
 lambda=[-0.8*1.5 0.5 -1.5*1.5 1];
 nocomp=2;
 specpos=rand(4,4);
 specent=rand(4,4);
 F=vplace2(A,B,lambda,nocomp,specpos,specent)       
 pl=eig(A-B*F)
%-----
 H=[1 0 0 0];
 Am=A-B*F;
 G=-inv(H*(Am\B));
 Bm=B*G;
%-----
 Mr=0.5; Mth=3/180*pi; 
 Tc=1; Tpen=0.25*2*pi*sqrt(2*ell/g);
 Q=diag([1/Mr^2,1/Mth^2,Tc^2/Mr^2,Tpen^2/Mth^2]);
 S=swflqr(A,B,Q) 
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=1;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0*0;0;0];
 sim('CIP_smcm_2015a.mdl')
%-----
%eof

固有ベクトル設定問題

投稿日時: 投稿者:

[0] 次の２入力２次系

\displaystyle{(1)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} 0 & 0\\ 0 & -1 \end{array}\right] }_{A} x+ \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{B} u }

に対して２つの異なる状態フィードバック

\displaystyle{(2)\quad %\begin{array}{l} u=- \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{F_1} x,\quad u=- \underbrace{ \left[\begin{array}{cc} 1 & 0\\ 1 & 2 \end{array}\right] }_{F_2} x %\end{array} }

による閉ループ系は、それぞれ

\displaystyle{(3)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} -2 & 0\\ 0 & -3 \end{array}\right] }_{A-BF_1} x,\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} -2 & 2\\ 0 & -3 \end{array}\right] }_{A-BF_2} x }

となって、同じ固有値-2,-3を持ちます。このように多入力系では閉ループ系の固有値を指定しても状態フィードバックは一意に定まりません。これは固有ベクトルの指定に任意性があるからです（ちなみに１入力系では一意に定まります）。そこで、上の２つの閉ループ系の固有ベクトルを求めてみると

\displaystyle{(4)\quad V_1=\left[\begin{array}{cc} 1 & 0\\ 0 & 1 \end{array}\right],\quad V_2=\left[\begin{array}{cc} 1 & 1\\ 0 & 1 \end{array}\right] }

また応答は次式となります。

\displaystyle{(5)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} e^{-2t} & 0\\ 0 & e^{-3t} \end{array}\right] }_{\exp((A-BF_1)t)=V_1\exp(\Lambda t)V_1^{-1}} x(0) }
\displaystyle{(6)\quad \dot{x}= \underbrace{ \left[\begin{array}{cc} 1 & 1\\ 0 & 1 \end{array}\right] \left[\begin{array}{cc} e^{-2t} & 0\\ 0 & e^{-3t} \end{array}\right] \left[\begin{array}{cc} 1 & -1\\ 0 & 1 \end{array}\right] }_{\exp((A-BF_2)t)=V_2\exp(\Lambda t)V_2^{-1}} x(0) }

明らかに(5)ではモードe^{-2t}とe^{-3t}が分離していますが、(6)ではそれらがミックスされています。このように、多入力系の場合の状態フィードバックの設計は、モードをどのように分布させるかがポイントとなります。

[1] 実固有値の場合

\displaystyle{(7)\quad (A-BF)v=\lambda v \Rightarrow \underbrace{ \left[\begin{array}{cc} \lambda I-A & B \end{array}\right] }_{G(\lambda)} \left[\begin{array}{c} v\\ \xi \end{array}\right] =0 \quad(\xi=Fv) }

より、vと\xiは、サイズn\times(n+m)の行列G(\lambda)の零化空間の基底行列に対して、適当なm次元ベクトル\deltaを用いて

\displaystyle{(8)\quad G(\lambda) \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right] =0 \Rightarrow \left[\begin{array}{c} v\\ \xi \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right]\delta }

のように表されます。そこで望ましい固有ベクトルvを与えて

\displaystyle{(9)\quad N(\lambda)\delta=v }

を解いて\deltaを求めます。ただし、望ましい固有ベクトルvのすべての成分を指定できるわけではないので

\displaystyle{(10)\quad \tilde{R}N(\lambda)\delta=\tilde{R}v =\left[\begin{array}{c} \tilde{v}^s\\ \tilde{v}^u\\ \end{array}\right] \Rightarrow \left[\begin{array}{cc} I_q & 0 \end{array}\right] \tilde{R}N(\lambda)\times\delta=\tilde{v}^s }

の最小２乗解として\delta=\delta^*を定めます。ここで、セレクタ行列\tilde{R}により、望ましい固有ベクトルvの指定されたq個の成分だけを集めた\tilde{v}^sが得られるとしています。

このとき、割当可能なv^*と\xi^*は次式で与えられます。

\displaystyle{(11)\quad \left[\begin{array}{c} v^*\\ \xi^* \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ M(\lambda)\\ \end{array}\right]\delta^* }

●もしA-\lambda I_nが正則ならば

\displaystyle{(12)\quad (\lambda I_n-A)\underbrace{(\lambda I_n-A)^{-1}B}_{N(\lambda)}+B\underbrace{(-I_m)}_{M(\lambda)}=0 }

より、N(\lambda)とM(\lambda)が特定できます。(8)より

\displaystyle{(13)\quad \begin{array}{l} v=(\lambda I_n-A)^{-1}B\delta\\ \xi=-\delta \end{array} \quad\Rightarrow\quad \boxed{v=(A-\lambda I_n)^{-1}B\xi} }

となって、指定した固有ベクトルvを近似する\xiを直接求めることができます。

●具体的に、(1)に対して、固有値-2,-3を持たせる状態フィードバックを、対応する固有ベクトルが(4)のV_1となるように求めてみます。

\displaystyle{(14)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} 0-(-2) & 0\\ 0 & -1-(-2) \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{(A-\lambda I_n)^{-1}B} \xi_1 =\underbrace{ \left[\begin{array}{c} 1\\ 0 \end{array}\right]}_{V_1(:,1)}\\ \Rightarrow \xi_1= \left[\begin{array}{cc} 1/2 & 1/2\\ 1 & -1 \end{array}\right]^{-1} \left[\begin{array}{c} 1\\ 0 \end{array}\right] = \left[\begin{array}{cc} 1 & 1/2\\ 1 & -1/2 \end{array}\right] \left[\begin{array}{c} 1\\ 0 \end{array}\right] = \left[\begin{array}{cc} 1\\ 1 \end{array}\right] \end{array} }

\displaystyle{(15)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} 0-(-3) & 0\\ 0 & -1-(-3) \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }_{(A-\lambda I_n)^{-1}B} \xi_2 =\underbrace{ \left[\begin{array}{c} 0\\ 1 \end{array}\right]}_{V_1(:,2)}\\ \Rightarrow \xi_2= \left[\begin{array}{cc} 1/3 & 1/3\\ 1/2 & -1/2 \end{array}\right]^{-1} \left[\begin{array}{c} 0\\ 1 \end{array}\right] = \left[\begin{array}{cc} 3/2 & 1\\ 3/2 & -1 \end{array}\right] \left[\begin{array}{c} 0\\ 1 \end{array}\right] = \left[\begin{array}{cc} 1\\ -1 \end{array}\right] \end{array} }

状態フィードバックゲインFは

\displaystyle{(16)\quad \left[\begin{array}{cc} \xi_1 & \xi_2 \end{array}\right] =FV_1 \quad\Rightarrow\quad F= \left[\begin{array}{cc} \xi_1 & \xi_2 \end{array}\right]V_1^{-1} = \left[\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right] }

[2] 複素共役固有値の場合

\displaystyle{(17)\quad \begin{array}{l} (A-BF)(v_R\pm jv_I)=(\lambda_R\pm j\lambda_I)(v_R\pm jv_I)\\ \Rightarrow (A-BF)v_R=\lambda_Rv_R-\lambda_Iv_I,\ (A-BF)v_I =\lambda_Rv_I+\lambda_Iv_R \\ \Rightarrow \underbrace{ \left[\begin{array}{ccc} \lambda_R I-A & \lambda_I I & B \end{array}\right] }_{G(\lambda)} \left[\begin{array}{cc} v_R & v_I\\ -v_I & v_R\\ \xi_R & \xi_I \end{array}\right] =0 \quad(\xi_R=Fv_R,\xi_I=Fv_I) \end{array} }

より、v_R,v_Iと\xi_R,\xi_Iは、サイズn\times(2n+m)の行列G(\lambda)の零化空間の基底行列に対して、適当なm次元ベクトル\delta_R,\delta_Iを用いて

\displaystyle{(18)\quad G(\lambda) \left[\begin{array}{c} N(\lambda)\\ P(\lambda)\\ M(\lambda) \end{array}\right] =0 \Rightarrow \left[\begin{array}{cc} v_R & v_I\\ -v_I & v_R\\ \xi_R & \xi_I \end{array}\right]= \left[\begin{array}{c} N(\lambda)\\ P(\lambda)\\ M(\lambda) \end{array}\right] \left[\begin{array}{cc} \delta_R & \delta_I \end{array}\right] }

のように表されます。そこで望ましいv_R,v_Iを与えて

\displaystyle{(19)\quad \left[\begin{array}{c} v_R \\ v_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} N(\lambda) \\ -P(\lambda) \end{array}\right] }_{\alpha} \delta_R= \underbrace{ \left[\begin{array}{cc} P(\lambda)\\ N(\lambda) \end{array}\right] }_{\beta} \delta_I }

の解として\delta_R,\delta_Iを求めますが、\alphaと\betaが独立ではないので、その求解には工夫が必要です。ここでは、Edwards and Spurgeon に従って、その手順を述べます。

まず、行列\alphaに対して、\alpha^Tの零化空間の基底行列をK_\alphaで表すと、\alpha^TK_\alpha=0よりK_\alpha^T\alpha=0を得て、K_\alpha^Tの零化空間の基底行列は\alpha、すなわち{\cal N}(K_\alpha^T)={\cal R}(\alpha)となることに注意します。同様に行列\betaに対して、\beta^Tの零化空間の基底行列をK_\betaで表すと、K_\beta^Tの零化空間の基底行列は\beta、すなわち{\cal N}(K_\beta^T)={\cal R}(\beta)となります。そこで、

\displaystyle{(20)\quad \underbrace{ \left[\begin{array}{c} K_\alpha^T\\ K_\beta^T \end{array}\right] }_{\gamma^T} K_\gamma=0 \Rightarrow {\cal R}(K_\gamma)\subset {\cal N}(K_\alpha^T)\cap {\cal N}(K_\beta^T)={\cal R}(\alpha)\cap {\cal R}(\beta) }

に注意して、\gamma=[K_\alpha,K_\beta]の零化空間の基底行列K_\gammaを用いて、望ましいv_R,v_Iを与えて

\displaystyle{(21)\quad \left[\begin{array}{c} v_R \\ v_I \end{array}\right]= K_\gamma \delta }

を解いて\deltaを求めます。ただし、望ましいv_R,v_Iのすべての成分を指定できるわけではないので

\displaystyle{(22)\quad \tilde{R}K_\gamma\delta=\tilde{R}v =\left[\begin{array}{c} \tilde{v}_R^s\\ \tilde{v}_I^s\\ \tilde{v}_R^u\\ \tilde{v}_I^u\\ \end{array}\right] \Rightarrow \left[\begin{array}{cc} I_{2q} & 0 \end{array}\right] \tilde{R}K_\gamma\times\delta= \left[\begin{array}{c} \tilde{v}_R^s\\ \tilde{v}_I^s\\ \end{array}\right] }

の最小２乗解として\delta=\delta^*を定めます。ここで、セレクター行列\tilde{R}により、望ましいv_R,v_Iの指定された2q個の成分だけを集めた\tilde{v}_R^s,\tilde{v}_I^sが得られるとしています。

このとき、割当可能なv_R^*,v_I^*は次式で与えられます。

\displaystyle{(23)\quad \left[\begin{array}{c} v_R^*\\ v_I^* \end{array}\right]= K_\gamma\delta^* }

これは(14)の\delta_R^*,\delta_I^*を、それぞれ\alpha^{\dag}K_\gamma\delta^*と\beta^{\dag}K_\gamma\delta^*で定めれば（\dagは疑似逆行列を表す）、対応する\xi_R^*,\xi_I^*は次式で与えられます。

\displaystyle{(24)\quad \left[\begin{array}{c} \xi_R^*\\ \xi_I^* \end{array}\right]= \left[\begin{array}{c} M(\lambda)\alpha^{\dag}K_\gamma\delta^*\\ M(\lambda)\beta^{\dag}K_\gamma\delta^*\\ \end{array}\right] }

[3] 複素共役固有値の場合（別のアプローチ）

(17)は次式のようにも表すことができます。

\displaystyle{(25)\quad \begin{array}{l} (A-BF)(v_R\pm jv_I)=(\lambda_R\pm j\lambda_I)(v_R\pm jv_I)\\ \Rightarrow (A-BF)v_R=\lambda_Rv_R-\lambda_Iv_I,\ (A-BF)v_I =\lambda_Rv_I+\lambda_Iv_R \\ \Rightarrow %\underbrace{ \left[\begin{array}{cccc} \lambda_R I-A & -\lambda_I I & B & 0\\ \lambda_I I & \lambda_R I-A & 0 & B\\ \end{array}\right] %}_{G(\lambda)} \left[\begin{array}{cc} v_R \\ v_I \\ \xi_R \\ \xi_I \end{array}\right] =0 \quad(\xi_R=Fv_R,\xi_I=Fv_I) \end{array} }

したがって、[2]において

\displaystyle{(26)\quad \begin{array}{l} G(\lambda)\leftarrow \left[\begin{array}{cccc} \lambda_R I-A & -\lambda_I I & B & 0\\ \lambda_I I & \lambda_R I-A & 0 & B \end{array}\right]\\ v\leftarrow\left[\begin{array}{c} v_R\\ v_I \end{array}\right],\ \xi\leftarrow\left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right],\ \delta\leftarrow\left[\begin{array}{c} \delta_R\\ \delta_I \end{array}\right] \end{array} }

と置き換えて同様の議論ができそうです（[2]の議論と等価であるかは検討中です）。

●もしA-\lambda_R I_nが正則ならば、公式

\displaystyle{(27)\quad \begin{array}{l} \left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right]^{-1} = \left[\begin{array}{cc} \Phi^{-1} & -\Phi^{-1}M_{12}M_{22}^{-1} \\ -M_{22}^{-1}M_{21}\Phi^{-1} & M_{22}^{-1}+M_{22}^{-1}M_{21}\Phi^{-1}M_{12}M_{22}^{-1} \end{array}\right]\\ \Phi\stackrel{\rm def}{=}M_{11}-M_{12}M_{22}^{-1}M_{21}\ (M_{22}\ {\rm nonsingular}) \end{array} }

を参照すると

\displaystyle{(28)\quad \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right]^{-1} }

を計算できます。そこで

\displaystyle{(29)\quad \begin{array}{l} \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right] \underbrace{ \left[\begin{array}{cc} \lambda_R I-A & -\lambda_I I\\ \lambda_I I & \lambda_R I-A \end{array}\right]^{-1} \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] }_{N(\lambda)}\\ + \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] \underbrace{(-I_{2m})}_{M(\lambda)}=0 \end{array} }

より、N(\lambda)とM(\lambda)が特定できます。(8)に相当する式より

\displaystyle{(30)\quad \boxed{ \left[\begin{array}{c} v_R\\ v_I \end{array}\right] = \left[\begin{array}{cc} A-\lambda_R I & \lambda_I I\\ -\lambda_I I & A-\lambda_R I \end{array}\right]^{-1} \left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right] \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right]} }

となって、指定した固有ベクトルの実部v_Rと虚部v_Iを近似する\xi_Rと\xi_Iを直接求めることができます。

●具体的に、(1)に対して、固有値-1\pm jを持たせる状態フィードバックを、対応する固有ベクトルが

\displaystyle{(31)\quad \underbrace{ \left[\begin{array}{c} 1\\ 0 \end{array}\right]}_{v_R} \pm j \underbrace{ \left[\begin{array}{c} 0\\ 1 \end{array}\right]}_{v_I} }

となるように求めてみます。

\displaystyle{(32)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cccc} 0-(-1) & 0 & 1 & 0\\ 0 & -1-(-1) & 0 & 1\\ -1 & 0 & 0-(-1) & 0\\ 0 & -1 & 0 & -1-(-1) \end{array}\right]^{-1} }_{\left[\begin{array}{cc} A-\lambda_R I & \lambda_I I\\ -\lambda_I I & A-\lambda_R I \end{array}\right]^{-1}}\\ \times\underbrace{ \left[\begin{array}{cccc} 1 & 1 & 0 & 0\\ 1 & -1& 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 1 & -1 \end{array}\right] }_{\left[\begin{array}{cc} B & 0\\ 0 & B \end{array}\right]} \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right] =\underbrace{ \left[\begin{array}{c} 1\\ 0\\ 0\\ 1 \end{array}\right]}_{\left[\begin{array}{c} v_R\\ v_I \end{array}\right]}\\ \Rightarrow \left[\begin{array}{c} \xi_R\\ \xi_I \end{array}\right]= \left[\begin{array}{cccc} 1/2 & 1/2 & -1/2 & -1/2\\ 0 & 0 & -1 & 1\\ 1/2 & 1/2 & 1/2 & 1/2\\ 1 & -1& 0 & 0 \end{array}\right]^{-1} \left[\begin{array}{c} 1\\ 0\\ 0\\ 1 \end{array}\right]= \left[\begin{array}{c} 1\\ 0\\ -1/2\\ -1/2 \end{array}\right] \end{array} }

状態フィードバックゲインFは

\displaystyle{(33)\quad \begin{array}{l} \left[\begin{array}{cc} \xi_R & \xi_I \end{array}\right] =F \left[\begin{array}{cc} v_R & v_I \end{array}\right]\\ \Rightarrow F= \left[\begin{array}{cc} \xi_R & \xi_I \end{array}\right] \left[\begin{array}{cc} v_R & v_I \end{array}\right]^{-1} = \left[\begin{array}{cc} 1 & -1/2\\ 0 & -1/2 \end{array}\right] \end{array} }

演習C42…Flipped Classroom
1^\circ　(16)と(33)の計算を行う次のプログラムを実行せよ。

MATLAB 
%test_vplace2.m
%-----
 clear all, close all
 A=[0 0;
    0 -1];
 B=[1 1;
    1 -1];
%-----
 lambda=[-2 -3];
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=[ 1  1; 
           1  1];       
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=[ 1  0;
           0  1];
 F=vplace2(A,B,lambda,nocomp,specpos,specent)
 [V,pl]=eig(A-B*F)
 VV=V*diag([1/V(2,1),1/V(1,2)])
 disp("---")
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=[ 1  1; 
           1  1];        
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=[ 1  0;
           0  1];        
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F) 
 VV=V*diag([1/V(1,1),1/V(1,2)])

2^\circ　１入力系の場合は、specposとspecentは乱数で与えてよい、なぜか？

MATLAB 
%test_vplace2a.m
%-----
 clear all, close all
 A=[0 1;
    0 0];
 B=[0;
    1];
%-----
%lambda=[-1 -1];
 lambda=[-2 -3];  
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F)
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(2,2)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(2,2) 
 F=vplace2(A,B,lambda,nocomp,specpos,specent) 
 [V,pl]=eig(A-B*F) 
%-----
%eof

Note C42 固有ベクトル設定プログラム
これはEdwards and Spurgeonによるプログラムspvplを基に作成できます。

スイッチング関数(2)

投稿日時: 投稿者:

スイッチング関数(2)

[1] 状態方程式とスイッチング関数は次の標準形をとるように座標変換されていると仮定します。

\displaystyle{(1)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot x(t)} = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) }

\displaystyle{(2)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} \ (M=S_2^{-1}S_1) }

このとき、等価制御

\displaystyle{(3)\quad u_{eq}(t)=-(SB)^{-1}SAx(t) }

による閉ループ系は次式で表されました。

\displaystyle{(4)\quad \dot{x}(t)=\underbrace{(I_n-B(SB)^{-1}S)A}_{A_{eq}}x(t)\quad(t\ge t_s) }

ここでA、B、Sは、(1)と(2)で与えられるとすると

\displaystyle{(5)\quad \begin{array}{l} A_{eq}= \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & 0 \\ \end{array}\right] }_{I_n-B(SB)^{-1}S} \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right]= \left[\begin{array}{cc} A_{11} & A_{12} \\ -MA_{11} & -MA_{12} \\ \end{array}\right]\\ = \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right] \left[\begin{array}{cc} \bar{A}_{11} & A_{12} \\ 0 & 0 \\ \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -M & I \\ \end{array}\right]^{-1} \end{array} }

となって、\bar{A}_{11}が安定行列なので、A_{eq}はm個の零固有値を持ちます。一方、A_{eq}の非零固有値\lambda_i\ (i=1,\cdots,n-m)に対応する固有ベクトルをv_iとすると、

\displaystyle{(6)\quad SA_{eq}=0\Rightarrow SA_{eq}v_i=0\Rightarrow \lambda_iSv_i=0 \Rightarrow Sv_i=0 }

すなわち

\displaystyle{(7)\quad { \boxed{S[v_1\cdots v_{n-m}]=0 %　\Rightarrow [v_1\cdots v_{n-m}]^TS^T=0}} }

が成り立ちます。したがって、まず等価制御による閉ループ系のA_{eq}が望ましい固有値・固有ベクトルをもつように設計して、(7)からSを決定することが考えられます。

● １入力系の場合は固有ベクトル設定の自由度はありませんが、多入力系の場合は固有ベクトルの向きを設定する必要があり、これは元の状態方程式(状態変数の並び)に対して与えます。そこで(7)が元の状態方程式ではどう変わるかを調べておきます。

元の状態方程式

\displaystyle{(8)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+f(t,x,u)\\ (x(t)\in{\rm\bf R}^n, u(t)\in{\rm\bf R}^m, f(t,x,u)\in{\rm\bf R}^n) \end{array} }

に対する標準形(1)は、適当な直交行列T_r座標による座標変換

\displaystyle{(9)\quad x_{r}(t)=T_rx(t) \Leftrightarrow x(t)=T_r^Tx_r(t) }

を行って（T_r^{-1}=T_r^Tに注意）

\displaystyle{(10)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{x}_{r1}(t)\\ \dot{x}_{r2}(t) \end{array}\right] }_{\dot{x}_r(t)} = \underbrace{ \left[\begin{array}{cc} A_{r11} & A_{r12} \\ A_{r21} & A_{r22} \\ \end{array}\right] }_{A_r=T_rAT_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \Sigma_1V^T \end{array}\right] }_{B_r=T_rB} u(t)\\ + \underbrace{ \left[\begin{array}{l} f_{ru}(t,x)\\ f_{rm}(t,x,u) \end{array}\right] }_{f_r(t,x,u)=T_rf(t,x,u)} \end{array} }

のように得られたものでした。また(2)は次式に

\displaystyle{(11)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_{r1} & S_{r2} \\ \end{array}\right] }_{S_r=ST_r^T} \underbrace{ \left[\begin{array}{c} x_{r1}(t)\\ x_{r2}(t) \end{array}\right] }_{x_r(t)=T_rx(t)} }

に、(7)は次式に相当します。

\displaystyle{(12)\quad S_r[v_{r1}\cdots v_{r,n-m}]=0 }

いま、元の状態空間での固有ベクトルv_1\cdots v_{n-m}は、座標変換後には

\displaystyle{(13)\quad Av_i=\lambda_iv_i\Rightarrow T_rAT_r^TT_rv_i=\lambda_iT_rv_i\Rightarrow v_{ri}=T_rv_i }

に変わることに注意します。(12)から

\displaystyle{(14)\quad S_rT_rT_r^T[v_{r1}\cdots v_{r,n-m}]=S[v_1\cdots v_{n-m}]=0 }

を得ます。これから元の状態方程式に対するスイッチング関数と固有ベクトルについても、(12)と同様の関係が成り立つと言えます。すなわち固有ベクトルv_1\cdots v_{n-m}は元の状態空間で与えておいて、(14)からSを決定してよいことがわかります。

演習C41…Flipped Classroom

1^\circ スイッチング関数決定プログラムswfvplを用いて、図1の台車駆動型倒立振子を安定化する次のSMC則を求めよ。

\displaystyle{{ \boxed{\begin{array}{l} u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }

図1 台車駆動型倒立振子

MATLAB 
%cCIP_smc2.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 C=diag([100 180/pi])*eye(2,4); 
%-----
 lambda=[-0.8 0.5 -3];
 nocomp=1;
 specpos=rand(4,3)
 specent=rand(4,3)
 S=swfvpl(A,B,lambda,nocomp,specpos,specent)
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=0.6;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0;0;0];
 sim('CIP_smc_2015a.mdl')
%-----
%eof

Note C41 スイッチング関数決定プログラム２
Edwards and Spurgeon によって、(26)から、スイッチング関数の行列Sを求めるプログラムswfvplが開発されています。

MATLAB 
%test_swfvpl.m
%-----
 clear all, close all
 A=[0 1 0;
    0 0 1;
    0 0 0];
 B=[0;
    0;
    1];
%-----
 lambda=[-2 -3];
 nocomp=0;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(3,3) 
 S=swfvpl(A,B,lambda,nocomp,specpos,specent) 
 Leq=inv(S*B)*S*A;
 [V,pl]=eig(A-B*Leq) 
%-----
 lambda=[-1 1];
 nocomp=1;
%固有ベクトルの要素の指定の有無 (１：指定、０：任意)
 specpos=rand(3,3)
%固有ベクトルの要素の値 (０：零指定、１：非零指定、－１：非零任意)       
 specent=rand(3,3) 
 S=swfvpl(A,B,lambda,nocomp,specpos,specent) 
 Leq=inv(S*B)*S*A;
 [V,pl]=eig(A-B*Leq) 
%-----
%eof

補遺：LQ制御

投稿日時: 投稿者:

LQ制御…Homework

[1]　可制御な制御対象

\displaystyle{(1)\quad %\left\{\begin{array}{ll} \dot{x}(t)=Ax(t)+Bu(t)\quad (x(t)\in{\bf R}^n,u(t)\in{\bf R}^m)\\ %z(t)=Cx(t)&(z(t)\in{\bf R}^q) %\end{array}\right. }

を安定化する状態フィードバック

\displaystyle{(2)\quad u(t)=-Fx(t) }

の決定法を考えます。一つの方法は，閉ループ系

\displaystyle{(3)\quad %\left\{\begin{array}{l} \dot{x}(t)=(A-BF)x(t) \\ %z(t)=Cx(t) %\end{array}\right. }

の時間応答に関する評価規範として，２次形式評価関数

(4)\quad \boxed{\begin{array}{l} \displaystyle{J=\int_0^\infty (z^T(t)Qz(t)+u^T(t)Ru(t))\,dt\quad (Q>0,R>0)}\\ \displaystyle{z(t)=Cx(t)\quad (z(t)\in{\bf R}^q)} \end{array}}

を設定し，これを最小化する問題を解くことです。ただし、(A,C)は可観測対とします。これによる状態フィードバックのゲイン行列Fは，リッカチ方程式

\displaystyle{(5)\quad {\boxed{\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0}} }

の解\Pi>0を用いて，次式で与えられます。

\displaystyle{(6)\quad {\boxed{F=R^{-1}B^T\Pi}} }

この証明はNoteに示しています。

Note A52-1 行列による微分　

いま、任意の行列Xの(i,j)要素をx_{ij}=[X]_{ij}で表すとき、スカラ関数fを行列変数Xの各要素で微分して得られる行列を\displaystyle{[\frac{\partial f}{\partial X}]_{ij}=\frac{\partial f}{\partial x_{ij}}}で定義します。このとき、行列のトレースについて、次が成り立ちます。

\displaystyle{(1)\ {\rm tr}AB={\rm tr}BA}
\displaystyle{(2)\ \frac{\partial}{\partial X}{\rm tr}AXB=A^TB^T}
\displaystyle{(3)\ \frac{\partial}{\partial X}{\rm tr}AX^TB=BA}
\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T=A^TXB^T+AXB}

実際、

\displaystyle{(1)\ \sum_{i}[AB]_{ii}=\sum_{i}\sum_{j}a_{ij}b_{ji} =\sum_{j}\sum_{i}b_{ji}a_{ij}=\sum_{j}[BA]_{jj}}

\displaystyle{(2)\ [\frac{\partial}{\partial X}{\rm tr}AXB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AXB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ij}b_{jk}}
\displaystyle{=\sum_{k}b_{jk}a_{ki}=[BA]_{ji}=[A^TB^T]_{ij}}

\displaystyle{(3)\ [\frac{\partial}{\partial X}{\rm tr}AX^TB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AX^TB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ji}b_{jk}}
\displaystyle{=\sum_{k}b_{ik}a_{kj}=[BA]_{ij}}

\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T =\frac{\partial}{\partial X}{\rm tr}(AXB)X^T +\frac{\partial}{\partial X}{\rm tr}AX(BX^T)}
\displaystyle{=AXB+A^TXB^T}

Note A52-2 LQ制御問題の解法　

可制御かつ可観測なn次系

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t) }

に対する状態フィードバック

\displaystyle{(2)\quad u(t)=-Fx(t) }

による閉ループ系

\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=(A-BF)x(t),\ z(t)=Cx(t)\\ A_F=A-BF:\ stable\ matrix \end{array} }

に対して、評価関数

\displaystyle{(4)\quad J=\int_0^\infty(z^T(t)Qz(t)+u^T(t)Ru(t))\,dt }

を最小化するようにFを決める問題を考えます。

閉ループ系における状態の振る舞いは次式で与えられます。

\displaystyle{(5)\quad x(t)=\exp({A_Ft})x(0) }

ここで、１次系の場合は初期状態はx(0)\ne0であればよかったのですが、一般の場合はインパルス応答となるようにBの列ベクトルB^{(i)}を考えます。各インパルス応答

\displaystyle{(6)\quad x^{(i)}(t)=\exp({A_Ft})B^{(i)} }

に対する評価関数Jの総和は

(7)\quad \begin{array}{l} \displaystyle{\sum\int_0^\infty x^{(i)}\,^T(t)(C^TQC+F^TRF)x^{(i)}(t)\,dt}\\ \displaystyle{=\sum B^{(i)}\,^T\underbrace{\int_0^\infty \exp(A_F^Tt)(C^TQC+F^TRF)\exp(A_Ft)\,dt}_{\Pi}B^{(i)}}\\ \displaystyle{=\sum B^{(i)}\,^T\Pi B^{(i)}={\rm tr}\ \Pi BB^T} \end{array} }

と書けます。いま\Gammaをラグランジュの未定定数として

\displaystyle{(8)\quad J'={\rm tr}\ \Pi x(0)x^T(0)+{\rm tr}\ ((\Pi A_F+A_F^T\Pi+C^TQC+F^TRF)\Gamma) }

を最小化する問題を考えます。ここで、制約条件は、リャプノフ方程式と呼ばれる

\displaystyle{(9)\quad \Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }

ですが、\Pi>0はA_Fが安定行列を意味することに注意します。

そこで、必要条件として次を得ます。

\displaystyle{(11)\quad \frac{\partial J'}{\partial\Pi}&=&BB^T+\Gamma A_F^T+A_F\Gamma=0\Rightarrow \Gamma>0 }
\displaystyle{(12)\quad \frac{\partial J'}{\partial F}&=&2(RF-B^T\Pi)\Gamma=0\Rightarrow F=R^{-1}B^T\Pi }
\displaystyle{(13)\quad \frac{\partial J'}{\partial\Gamma}&=&\Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }

ここで、第２式から得られるF=R^{-1}B^T\Piを第３式に代入して

\displaystyle{(14)\quad \begin{array}{l} \Pi (A-BR^{-1}B^T\Pi)+(A-BR^{-1}B^T\Pi)^T\Pi+C^TQ\nonumber\\ +(R^{-1}B^T\Pi)^TRR^{-1}B^T\Pi=0 \end{array} }

すなわち、リッカチ方程式と呼ばれる\Piの行列方程式

\displaystyle{(15)\quad \Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0 }

を得ます。これから\Pi>0を求めて、Fは

\displaystyle{(16)\quad F=R^{-1}B^T\Pi }

のように得られます。このような制御方式をLQ制御と呼びます。

一方、十分性の議論は次のように行われます。まず、被積分項は次のように表すことができます。

\displaystyle{(17)\quad \begin{array}{lll} &&y^TQy+u^TRu=(u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\nonumber\\ &&-\underbrace{\frac{d}{dt}x^T\Pi x}_{2x^T\Pi\dot{x}} \end{array} }

実際、右辺に\dot{x}=Ax+Buを代入し、リッカチ方程式を用いると

\displaystyle{(18)\quad \begin{array}{lll} &&{\rm RHS}=u^TRu+2x^T\Pi Bu+x^T\Pi BR^{-1}B^T\Pi x-2x^T\Pi(Ax+Bu)\nonumber\\ &&=u^TRu+x^T(\Pi A+A^T\Pi+C^TQC)x-2x^T\Pi Ax={\rm LHS} \end{array} }

したがって、上記の両辺を積分して

\displaystyle{(19)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt-\left[x^T\Pi x\right]_0^\infty }

を得ます。ここで、x(t)=\exp(A_Ft)x(0)\rightarrow 0\ (t\rightarrow\infty)を前提とするので

\displaystyle{(20)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt+x^T(0)\Pi x(0) }

を得ます。これからu=-R^{-1}B^T\Pi xが評価関数を最小化することが分かります。

スイッチング関数(1)

投稿日時: 投稿者:

スイッチング関数(1)

[1] 状態方程式とスイッチング関数は次の標準形をとるように座標変換されていると仮定します。

\displaystyle{(1)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot x_2(t) \end{array}\right] }_{\dot x(t)} = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t) }

\displaystyle{(2)\quad s(t)= \underbrace{ \left[\begin{array}{cc} S_1 & S_2 \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{S_2 \left[\begin{array}{cc} M & I \\ \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} \ (M=S_2^{-1}S_1) }

これに対して、座標変換

\displaystyle{(3)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)}\\ \Leftrightarrow \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} = \underbrace{ \left[\begin{array}{cc} I & 0 \\ -S_2^{-1}S_1 & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} \end{array} }

を行って、次式を得ます。

\displaystyle{(4)\quad \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)} = \underbrace{ \left[\begin{array}{cc} \bar{A}_{11} & \bar{A}_{12} \\ \bar{A}_{21} & \bar{A}_{22}\\ \end{array}\right] }_{\bar{A}=T_s A T_s^{-1}} \underbrace{ \left[\begin{array}{c} x_1(t)\\ s(t) \end{array}\right] }_{\bar{x}(t)} + \underbrace{ \left[\begin{array}{c} 0\\ \bar{B}_2 \end{array}\right] }_{\bar{B}=T_sB} u(t) }

\displaystyle{(5)\quad \begin{array}{l} \bar{A}_{11}=A_{11}-A_{12}M\\ \bar{A}_{12}=A_{12}S_2^{-1}\\ \bar{A}_{21}=S_2(M\bar{A}_{11}+A_{21}-A_{22}M)\\ \bar{A}_{22}=S_2(M{A}_{12}+A_{22})S_2^{-1}\\ \bar{B}_{2}=S_2B_2 \end{array} }

ここで、s(t)=0、すなわち、状態が超平面内に拘束されているとすると

\displaystyle{(6)\quad s(t)=0 \Rightarrow Sx(t)=S_1x_1(t)+S_2x_2(t)=0\Rightarrow x_2(t)=-\underbrace{S_2^{-1}S_1}_Mx_1(t)}}

となって、x_2の振る舞いはx_1によって決まります。一方、x_1の振る舞いは

\displaystyle{(7)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}M)}_{\bar{A}_{11}}x_1(t)}

によって決まります。これは

\displaystyle{(8)\quad \dot{x}_1(t)=A_{11}x_1(t)+A_{12}x_2(t)}

に対する安定化状態フィードバック

\displaystyle{(9)\quad x_2(t)=-Mx_1(t)}

による閉ループ系とみなすことができます。

したがって、\bar{A}_{11}が安定行列となるようにスイッチング関数を選ぶためには、(8)に対する状態フィードバックによる安定化問題を解けばよいことが分かります。

[2]　そこで、次の２次形式評価関数を最小にするようにMを決定する方法が提案されています。

\displaystyle{(10)\quad \boxed{J=\int_0^\infty \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]^T }_{x^T(t)} \underbrace{ \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} dt} }

これは、恒等式

\displaystyle{(11)\quad \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right]}
\displaystyle{ =\left[\begin{array}{cc} I & Q_{12}Q_{22}^{-1} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} Q_{11}-Q_{12}Q_{22}^{-1}Q_{21} & 0 \\ 0 & Q_{22} \end{array}\right] \left[\begin{array}{cc} I & 0 \\ Q_{22}^{-1}Q_{21} & I \end{array}\right] }

を用いて

\displaystyle{(12)\quad \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right]^T \left[\begin{array}{cc} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{array}\right] \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] =\left[\begin{array}{c} x_1(t)\\ Q_{12}Q_{22}^{-1}x_1(t)+x_2(t) \end{array}\right]^T }
\displaystyle{ \times \left[\begin{array}{cc} Q_{11}-Q_{12}Q_{22}^{-1}Q_{21} & 0 \\ 0 & Q_{22} \end{array}\right] \left[\begin{array}{c} x_1(t)\\ Q_{12}Q_{22}^{-1}x_1(t)+x_2(t) \end{array}\right] }

より

\displaystyle{(13)\quad J=\int_0^\infty(x_1^T(t)\underbrace{(Q_{11}-Q_{12}Q_{22}^{-1}Q_{21})}_{\hat{Q}}x_1(t)}
\displaystyle{+\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))^T}_{\hat{u}^T}\underbrace{Q_{22}}_{\hat{R}}\underbrace{(Q_{12}Q_{22}^{-1}x_1(t)+x_2(t))}_{\hat{u}})dt}

と書けます。ここで、(8)と(6)をそれぞれ

\displaystyle{(14)\quad \dot{x}_1(t)=\underbrace{(A_{11}-A_{12}Q_{22}^{-1}Q_{21})}_{\hat{A}}x_1(t)+\underbrace{A_{12}}_{\hat{B}}\underbrace{(Q_{22}^{-1}Q_{21}x_1(t) +x_2(t))}_{\hat{u}(t)}}

\displaystyle{(15)\quad \underbrace{Q_{22}^{-1}Q_{21}x_1(t) +x_2(t)}_{\hat{u}(t)}=-\underbrace{(M-Q_{22}^{-1}Q_{21})}_{\hat{F}}x_1(t)}

と書き直してみますと、次のような最適化制御問題の定式化ができていることがわかります。すなわち

\displaystyle{(16)\quad \boxed{\begin{array}{l} \dot{x}_1(t)=\hat{A}x_1(t)+\hat{B}\hat{u}(t)\\ \hat{A}=A_{11}-A_{12}Q_{22}^{-1}Q_{21}\\ \hat{B}=A_{12} \end{array}} }

を２次系式評価関数

\displaystyle{(17)\quad \boxed{\begin{array}{l} J=\int_0^\infty(x_1^T(t)\hat{Q}x_1(t)+\hat{u}^T(t)\hat{R}\hat{u}(t))dt\\ \hat{Q}=Q_{11}-Q_{12}Q_{22}^{-1}Q_{21}\\ \hat{R}=Q_{22} \end{array}} }

を最小にするように

\displaystyle{(18)\quad \boxed{\begin{array}{l} \hat{u}(t)=-\hat{F}x_1(t)\\ \hat{F}=M-Q_{22}^{-1}Q_{21} \end{array}} }

を求める問題です。この解は、リッカチ方程式

\displaystyle{(19)\quad \boxed{\Pi\hat{A}+\hat{A}^T\Pi-\Pi\hat{B}\hat{R}^{-1}\hat{B}^T\Pi+\hat{Q}=0}}

を解いて

\displaystyle{(20)\quad \boxed{\hat{F}=\hat{R}^{-1}\hat{B}^T\Pi}}

を求め、次式からMを決めればよいことがわかります。

\displaystyle{(21)\quad \boxed{M=\hat{F}+Q_{22}^{-1}Q_{21}}}

演習C31…Flipped Classroom

1^\circ スイッチング関数決定プログラムswflqrを用いて、図1の台車駆動型倒立振子を安定化する次のSMC則を求めよ。

\displaystyle{{ \boxed{\begin{array}{l} u(t)=\underbrace{u_\ell(t)}_{linear\ control}+\underbrace{u_n(t)}_{switching\ component}\\ u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t) \quad(\Phi:stable\ matrix)\\ u_n(t)=-\underbrace{(SB)^{-1}\rho(t,x)}_{L_n}\frac{P_2s(t)}{||P_2s(t)||}\quad(P_2\Phi+\Phi^TP_2=-I) \end{array}}} }


図1 台車駆動型倒立振子

MATLAB 
%cCIP_smc1.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 ths=0; %input('ths   = <0,180> ')/180*pi;
 th0=3/180*pi; %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
 C=diag([100 180/pi])*eye(2,4); 
%-----
 Mr=0.5; Mth=3/180*pi; 
 Tc=1; Tpen=0.25*2*pi*sqrt(2*ell/g);
 Q=diag([1/Mr^2,1/Mth^2,Tc^2/Mr^2,Tpen^2/Mth^2]);
 S=swflqr(A,B,Q)
%-----
 Phi=-0.1;
 P2=0.5*inv(-Phi);
 Check=P2*Phi+Phi*P2
 P2S=P2*S;
 Leq=inv(S*B)*S*A;
 LPhi=-inv(S*B)*Phi*S;
 L=Leq+LPhi;
 rho=1.5;
 Ln=inv(S*B)*rho;
%-----
 x0=[0;th0;0;0];
 sim('CIP_smc_2015a.mdl')
%-----
%eof

Note C31 スイッチング関数決定プログラム１
Edwards and Spurgeon によって、(19)、(20)、(21)から、スイッチング関数の行列Sを求めるプログラムswflqrが開発されています。

MATLAB 
%test_swflqr.m
%-----
 clear all, close all
 A=[0 1 0;
    0 0 1;
    0 0 0];
 B=[0;
    0;
    1];
%-----
 Q=eye(3,3);
 S=swflqr(A,B,Q)
%-----
%eof