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

SM制御則

投稿日時: 投稿者:

SM制御則…Homework

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

\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, 1\le m< n, {\rm rank}B=m) \end{array} }

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

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

以下では、状態方程式とスイッチング関数は次のSM標準形(regular form)をとるように座標変換されているとします(Note C22-1参照)。

\displaystyle{(3)\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{(4)\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_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)} }

ここで、B_2とS_2はともにサイズm\times mの正則行列とします。したがって、f_mはマッチング条件を満たしますが、f_uはマッチング条件を満たさないことに留意します。

これに対して、座標変換

\displaystyle{(5)\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} }

を行うために、まず(5)を(3)に代入して

\displaystyle{(6)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & S_2^{-1} \\ \end{array}\right] }_{T_s^{-1}} \underbrace{ \left[\begin{array}{c} \dot x_1(t)\\ \dot s(t) \end{array}\right] }_{\dot{\bar{x}}(t)}\\ = \underbrace{ \left[\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{cc} I & 0 \\ -M & 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)} + \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)} \end{array} }

左からT_sをかけて

\displaystyle{(7)\quad \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} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{cc} A_{11}-A_{12}M & A_{12}S_2^{-1} \\ A_{21}-A_{22}M & A_{22}S_2^{-1} \\ \end{array}\right] }_{AT_s^{-1}} \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} 0\\ B_2 \end{array}\right] }_{B} u(t) + \underbrace{ \left[\begin{array}{cc} I & 0 \\ S_1 & S_2 \\ \end{array}\right] }_{T_s} \underbrace{ \left[\begin{array}{l} f_u(t,x)\\ f_m(t,x,u) \end{array}\right] }_{f(t,x,u)} \end{array} }

すなわち、次式を得ます。

\displaystyle{(8)\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{(8')\quad \left\{\begin{array}{l} \bar{A}_{11}=A_{11}-A_{12}M\quad(M=S_2^{-1}S_1)\\ \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}がサイズn-m\times n-mの安定行列となるようにスイッチング関数が選ばれていると仮定します。

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

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

を、２次安定性

\displaystyle{(10)\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{(11)\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] 可到達性の検討

ここでは、スライディングモード制御則(9)の具体的な表現を求め、スライディングモードに到達する時刻t_sの見積もりを行います。

●等価制御は、(7)においてf_u=0, f_m=0として

\displaystyle{ \begin{array}{cl} (12.1) & s(t)=0\\ \Downarrow &\\ (12.2) & \dot{s}(t)=0\\ \Downarrow &\\ (12.3) & \bar{A}_{21}x_1(t)+\bar{A}_{22}s(t)+\bar{B}_{2}u(t)=0\\ \Downarrow &\\ (12.4) & u_{eq}(t)=-\underbrace{\bar{B}_{2}^{-1}}_{([0\ I]T_sB)^{-1}} \underbrace{\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]}_{[0\ I]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]T_s, x=T_s^{-1}\bar{x}\\ (12.5) & u_{eq}(t)=-(SB)^{-1}SAx(t)} \end{array} }

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

\displaystyle{(13.1)\quad u_\ell(t)=-\bar{B}_{2}^{-1} (\left[\begin{array}{cc} \bar{A}_{21} & \bar{A}_{22} \\ \end{array}\right]\bar{x}(t)-\Phi s(t)) }

すなわち

\displaystyle{(13.2)\quad \boxed{u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)}} }

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

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

\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} \\ 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} }

すなわち

\displaystyle{ \begin{array}{l} (15.1)\quad \dot{x}_1(t)=\bar{A}_{11}x_1(t)+\bar{A}_{12}s(t)+f_u(t,x)\\ (15.2)\quad \dot{s}(t)=\Phi s(t)+\bar{B}_{2}u_n(t)+\underline{S_2(Mf_u(t,x)+f_m(t,x,u))} \end{array} }

閉ループ系の状態はx_1とsから成りますが、(15.2)はsだけに関係しているので、まずその下線部に抗してスライディングモードs=0を達成するu_nを明らかにします。

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

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

を満たすP_2>0を選ぶことができます。これを用いて、(7)のスイッチング項u_nを

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

すなわち

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

と定めます。

●\rho(t,x)は、仮定

{\displaystyle{ \begin{array}{l} (18.1)\quad ||f_u(t,x)|| \le k_1||x(t)||+k_2\\ (18.2)\quad ||f_m(t,x,u)|| \le k_3||u(t)||+\alpha(t,x)\\ (18.3)\quad k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}|| < 1 \end{array} }

のもとで

\displaystyle{(19)\quad \boxed{ %\begin{array}{l} \rho(t,x)\ge \frac{||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}{1-k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||} %\end{array} }}

のように選びます。このとき

\displaystyle{ \begin{array}{l} (20.1)\quad S_2=S_2B_2 B_2^{-1} \Rightarrow ||S_2||<||\bar{B}_{2}||\,||B_2^{-1}||\\ (20.2)\quad ||u_n(t)||\le\rho(t,x)||\bar{B}_{2}^{-1}||\\ (20.3)\quad ||u(t)||\le ||u_\ell(t)||+||u_n(t)|| \end{array} }

に注意して

\displaystyle{ \begin{array}{cl} (21.1) & \rho(t,x)\ge k_3 \rho(t,x)||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||\\ & +||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (20.1)\\ (21.2) & \ge k_3 \rho(t,x)||S_2||\,||\bar{B}_{2}^{-1}||\\ & + ||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & \\ (21.3) & = ||S_2|| \{||M||(k_1||x(t)||+k_2)\\ & +k_3(||u_\ell(t)||+ \rho(t,x)||\bar{B}_{2}^{-1}||)+\alpha(t,x)\}+\gamma_2}\\ \Downarrow & by\ (20.2) \\ (21.4) & \ge ||S_2|| (||M||(k_1||x(t)||+k_2) +k_3(||u_\ell(t)||+||u_n(t)||)+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (20.3) \\ (21.5) & \ge ||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u(t)||+\alpha(t,x))+\gamma_2}\\ \Downarrow & by\ (18.1),(18.2) \\ (21.6) & \ge ||S_2|| (||M||\,||f_u(t,x)||+||f_m(t,x,u)||)+\gamma_2} \end{array} }

を得て、次式が成り立ちます。(15.2)はsだけに関係しているので、リャプノフ関数V(s)=s^T(t)P_2s(t)の微分について調べます。

\displaystyle{ \begin{array}{cl} (22.1) & \dot{V}(s)=2s^T(t)P_2\dot{s}(t)\\ \Downarrow & by\ (15.2) \\ (22.2) & =2s^T(t)P_2 (\Phi s(t)-\rho(t,x)\frac{P_2s(t)} {||P_2s(t)||}+S_2(Mf_u(t,x)+f_m(t,x,u)))\\ \Downarrow & \\ (22.3) & =s^T(t)(P_2\Phi+\Phi^TP_2)s(t)\\ & +2s^T(t)P_2(-\rho(t,x)\frac{P_2s(t)} {||P_2s(t)||}+S_2(Mf_u(t,x)+f_m(t,x,u)))\\ \Downarrow & by\ (16) \\ (22.4) & \le-||s(t)||^2\\ & -2||P_2s(t)||(\rho(t,x)-||S_2||(||M||||f_u(t,x)||+||f_m(t,x,u)||)\\ \Downarrow & by\ (21.6) \\ (22.5) & \le -||s(t)||^2-2\gamma_2||P_2s(t)||\\ (22.6) & \le -||s(t)||^2 \end{array} }

これより(15.2)における２次安定性が成り立つことが分かります。２次安定であれば状態sは0に漸近します。しかし興味深いのは、ある時刻t_s以降はs=0を保持するような振舞いをすることです。

●その到達時刻t_sを見積もるために、(22.5)から

\displaystyle{(23)\quad \dot{V}(s)　\le -2\gamma_2||P_2s(t)|| }

も成り立つことと

\displaystyle{(24)\quad ||P_2s(t)||^2=(P_2^{1/2}s(t))^TP_2(P_2^{1/2}s(t)) \ge \sigma_n(P_2)||P_2^{1/2}s(t)||^2=\sigma_n(P_2)V(s) }

に注意して

\displaystyle{(25)\quad \dot{V}(s)\le -2\gamma_2\sigma_n^{1/2}(P_2)V^{1/2}(s) }

よって

\displaystyle{(26)\quad \begin{array}{l} \dot{V}(s)=\frac{d}{dt}(V^{1/2}(s))^2=2V^{1/2}(s)\frac{d}{dt}V^{1/2}(s)< -2\gamma_2\sigma_n^{1/2}(P_2) V^{1/2}(s)\\ \Rightarrow \frac{d}{dt}V^{1/2}(s)< -\gamma_2\sigma_n^{1/2}(P_2) \end{array} }

この両辺を積分すると

\displaystyle{(27)\quad \begin{array}{l} \underbrace{V^{1/2}(s(t_s))}_0-V^{1/2}(s(0))< -\gamma_2\sigma_n^{1/2}(P_2) t_s\\ \Leftrightarrow \boxed{t_s\le \frac{V^{1/2}(s(0))}{\gamma_2\sigma_n^{1/2}(P_2)}} \end{array} }

を得て、s=0に到達する時刻t_sを見積もることができます。

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

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

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

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

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

\displaystyle{ \begin{array}{cl} (29.1)&\dot{V}(x_1)=2x_1^T(t)P_1\dot{x}_1(t)\\ \Downarrow & by\ (7) \\ (29.2)&=2x_1^T(t)P_1(\bar{A}_{11}x_1(t)+\bar{A}_{12}\underbrace{s(t)}_{0}+f_u(t,x))\\ \Downarrow & \\ (29.3)&=x_1^T(t)(P_1\bar{A}_{11}+\bar{A}_{11}^TP_1)x_1(t)+2x_1^T(t)P_1f_u(t,x)\\ \Downarrow & by\ (28) \\ (29.4)&=-x_1^T(t)Q_1x_1(t)+2x_1^T(t)P_1f_u(t,x)\\ \Downarrow &by\ \sigma_n(Q_1)x_1^Tx_1\le x_1^TQ_1x_1\\ (29.5)&\le -\sigma_n(Q_1)x_1^T(t)x_1(t)+2||x_1^T(t)P_1||\,||f_u(t,x)||\\ \Downarrow &\\ (29.6)&\le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{x_1^T(t)P_1^2x_1(t)}\,||f_u(t,x)||\\ \Downarrow &by\ x_1^TP_1^2x_1\le\sigma_1(P_1^2)x_1^Tx_1\\ (29.7)&\le -\sigma_n(Q_1)||x_1(t)||^2+2\sqrt{\sigma_1(P_1^2)x_1^T(t)x_1(t)}\,||f_u(t,x)||\\ \Downarrow &by\ \sigma_1(P_1^2)=\sigma_1^2(P_1)\\ (29.8)&= -\sigma_n(Q_1)||x_1(t)||^2+2\sigma_1(P_1)||x_1(t)||\,||f_u(t,x)||\\ \Downarrow &\\ (29.9)&= -\sigma_1(P_1)||x_1(t)||(\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)||-2||f_u(t,x)||) \end{array} }

ここで

\displaystyle{(30)\quad \boxed{ ||f_u(t,x)||<\frac{1}{2}\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)|| }}

のとき

\displaystyle{(31)\quad \dot{V}(x_1)\le -\sigma_1(P_1)||x_1(t)||(\frac{\sigma_n(Q_1)}{\sigma_1(P_1)}||x_1(t)||-2||f_u(t,x)||)<0 }

となって、リヤプノフ安定性が成り立ちます。したがって、ロバスト性での議論と同様に、x_1が発散することはないことが分かります。

Note C22-1 SM標準形について

列フルランクをもつBは、サイズn\times nの直交行列Uとサイズm\times mの直交行列Vを用いて

\displaystyle{(1)\quad B=U\Sigma V^T}

のように特異値分解できます。ただし

\displaystyle{(2)\quad \Sigma= \left[\begin{array}{cc} \Sigma_1 \\ 0_{n-m\times m} \end{array}\right] }

\displaystyle{(3)\quad \Sigma_1={\rm diag}\{\sigma_1,\cdots,\sigma_m\}\quad(\sigma_1\ge\cdots\ge\sigma_m>0) }

いま

\displaystyle{(4)\quad U= \left[\begin{array}{cc} U_1(n\times m) & U_2(n\times n-m) \end{array}\right] }

と分割するとき、本論(1)に対する座標変換

\displaystyle{(5)\quad x_{r}(t)= \underbrace{ \left[\begin{array}{cc} U_2 & U_1 \end{array}\right] }_{T_r} x(t) \Leftrightarrow x(t)=T_r^Tx_r(t) }

を行うと（T_r^{-1}=T_r^Tに注意）

\displaystyle{(6)\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} }

のように本論(3)が得られます。ちなみに本論(4)は次式に相当します。

\displaystyle{(7)\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)} }

これから、標準形に対して求めたスイッチング関数を元の状態方程式に関するものに戻すには

\displaystyle{(8)\quad S=S_rT_r }

とすればよいことが分かります。

Note C22-2 条件(18)について
次式で表されるDCモータを考えます。

\displaystyle{(1)\quad \left\{\begin{array}{l} J_0\ddot{\theta}=K_ti_a\\ L_0\dot{i}_a+Ri_a=v_a-K_e\dot{\theta} \end{array}\right. }

これより、次の状態方程式を得ます。

\displaystyle{(2)\quad \underbrace{ \left[\begin{array}{c} \dot{\theta}\\ \dot{\omega}\\ \dot{i}_a \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cc|c} 0 & 1 & 0\\ 0 & 0 & \frac{K_t}{J_0}\\\hline 0 & -\frac{K_e}{L_0} & -\frac{R}{L_0} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \theta\\ \omega\\ i_a \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{c} 0\\ 0\\\hline \frac{1}{L_0} \end{array}\right] }_{B} \underbrace{ v_a }_{u} }

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

\displaystyle{(3)\quad s= \underbrace{ \left[\begin{array}{ccc} \frac{J_0}{K_t}\omega_n^2 & 2\frac{J_0}{K_t}\zeta\omega_n & 1 \end{array}\right] }_{S}x \quad(M= \left[\begin{array}{ccc} \frac{J_0}{K_t}\omega_n^2 & 2\frac{J_0}{K_t}\zeta\omega_n \end{array}\right]) }

さらに次式のようなパラメータの不確かさを考えます。

\displaystyle{(4)\quad \left\{\begin{array}{l} \underbrace{(1+\xi_J)J_0}_{J}\ddot{\theta}=K_ti_a\\ \underbrace{(1+\xi_L)L_0}_{L}\dot{i}_a+Ri_a=v_a-K_e\dot{\theta} \end{array}\right. \Rightarrow \left\{\begin{array}{l} \ddot{\theta}=\frac{K_t}{J_0}i_a\underbrace{-\xi_J\ddot{\theta}}_{f_{u}}\\ \dot{i}_a=-\frac{R}{L_0}i_a-\frac{K_e}{L_0}\dot{\theta}+\frac{1}{L_0}v_a\underbrace{-\xi_L\dot{i}_a}_{f_m} \end{array}\right. }

ただし、次を仮定します。

\displaystyle{(5)\quad |\xi_J|<0.5,\ |\xi_L|<0.1 }

これにより状態方程式は次式となります。

\displaystyle{(6)\quad \dot{x}=Ax+Bu+ \left[\begin{array}{c} 0\\ f_{u}\\ f_{m} \end{array}\right]%} }

まずf_mは

\displaystyle{(7)\quad f_m=-\xi_L\dot{i}_a=-\xi_L(-\frac{R}{L_0}i_a-\frac{K_e}{L_0}\dot{\theta}+\frac{1}{L_0}v_a) }

と表されるので、仮定(18)の第２式におけるk_3と\alphaは次のように特定できます。

\displaystyle{(8)\quad |f_m|<\frac{\xi_L}{L_0}|v_a|+ \frac{\xi_L}{L_0}(R|i_a|+K_e|\dot{\theta}|) <\underbrace{\frac{1}{10L_0}}_{k_3}|v_a|+ \underbrace{\frac{1}{10L_0}(R|i_a|+K_e|\dot{\theta}|)}_{\alpha} }

またf_uについては、閉ループ系では

\displaystyle{(9)\quad \begin{array}{l} f_u=-\xi_J\dot{\omega} =-\xi_J\frac{K_t}{J_0}(i_a+\frac{J_0}{K_t}\omega_n^2\theta+2\frac{J_0}{K_t}\zeta\omega_n\dot{\theta}) \end{array} }

と考えられるので、仮定(18)の第１式におけるk_1とk_2は次のように特定できます。

\displaystyle{(10)\quad |f_u|<\xi_J||M||||x|| <\underbrace{0.5||M||}_{k_1}||x||+\underbrace{0}_{k_2} }

さて、(19)を見積もってみます。

\displaystyle{(11)\quad \begin{array}{l} \rho(t,x)\ge \frac{||S_2|| (||M||(k_1||x(t)||+k_2)+k_3||u_\ell(t)||+\alpha(t,x))+\gamma_2}{1-k_3 ||\bar{B}_{2}||\,||\bar{B}_{2}^{-1}||\,||B_2^{-1}||}\\ =\frac{0.5||M||^2||x(t)||+\frac{1}{10L_0}||u_\ell(t)||+\frac{1}{10L_0}(R|i_a|+K_e|\dot{\theta}|)+\gamma_2}{1-\frac{1}{10L_0} L_0}\\ =\frac{5L_0||M||^2||x(t)||+||u_\ell(t)||+R|i_a|+K_e|\dot{\theta}|+10L_0\gamma_2}{9L_0} \end{array} }

以上に基づいて、DCモータに対するSMCのシミュレーションを行ってみます。

MATLAB 
%cDCM_smc.m
%-----
 clear all, close all
 R=1.2; L0=0.05; Ke=0.6; Kt=0.6; J=0.1352;
 A=[0 1 0;
    0 0 Kt/J;
    0 -Ke/L0 -R/L0];
 B=[0;0;1/L0];
 C=eye(3); D=zeros(3,1);
 A11=A(1:2,1:2); A12=A(1:2,3);
 A21=A(3,1:2);   A22=A(3,3);
%
 L=0.046;
 Ad=[0 1 0;
     0 0 6;
     0 -Ke/L -R/L];
 Bd=[0;0;1/L];
%-----
 wn=2; zeta=0.95;
 M=[wn^2 2*zeta*wn]/A12(2,1);
 S=[M 1];
 Leq=inv(S*B)*S*A
%
 Phi=-2
 LPhi=-inv(S*B)*Phi*S
% 
 P2=-1/(Phi*2);
 P2S=P2*S;
 gam2=0.01
 rho=(1+Ke*1+R*1+5*L0*norm(M)^2*1+10*L0*gam2)/(9*L0)
 Ln=inv(S*B)*rho
%-----
 x0=[1;0;0];
 sim('DCM_smc')
%-----
%eof


図1 DCモータに対するSMCのシミュレーション例

Note C22-4 スライディングモード制御系の安定性について

スライディングモード制御系の振舞いは、スライディングモードへの突入動作と、そこでの滑り動作から成ります。前者は、s=0を達成することですから、２次安定性

\displaystyle{(1)\quad V(s)=s^T(t)P_2s(t)&\quad\Rightarrow \dot{V}(s)\le -s^T(t)s(t)&}\quad(t\le t_s) }

を検討します。これが可到達条件となります。ただ、これではまだ原点x=0への収束は保証されないので、２次安定性

\displaystyle{(2)\quad V(x_1)=x_1^T(t)P_1x_1(t)&\quad\Rightarrow \dot{V}(x_1)\le -x_1^T(t)Q_1x_1(t)&}\quad(t> t_s) }

を検討します。

Note C22-5 安定行列\bar{A}_{11}=A_{11}-A_{12}M\ (M=S_2^{-1}S_1)について

●「(A,B)が可制御対 \Leftrightarrow (A_{11},A_{12})も可制御対」

実際、
\displaystyle{(1)\quad \begin{array}{l} {\rm rank}\left[\begin{array}{cc} B & \lambda I_n-A \end{array}\right]\\ = {\rm rank}\left[\begin{array}{ccc} 0 & A_{11}-\lambda I_m & A_{12}\\ B_{2} & A_{21} & A_{22}-\lambda I_m \end{array}\right]\\ = {\rm rank}\left[\begin{array}{ccc} A_{11}-\lambda I_m & A_{12} \end{array}\right]+m \end{array} }

これより

\displaystyle{(2)\quad {\rm rank}\left[\begin{array}{cc} B & \lambda I_n-A \end{array}\right]=n \Leftrightarrow {\rm rank}\left[\begin{array}{ccc} A_{12} & A_{11}-\lambda I_m \end{array}\right]=n-m }

●「\bar{A}_{11}の固有値はスイッチング関数を出力方程式としたときの不変零点に等しい」

まず、線形系

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

の不変零点とは、行列\left[\begin{array}{cc} A & B\\ C & 0 \end{array}\right]の一般化固有値すなわち

\displaystyle{(3)\quad \left[\begin{array}{cc} A & B\\ C & 0 \end{array}\right] \left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right] =\lambda \left[\begin{array}{cc} I_n & 0\\ 0 & 0 \end{array}\right] \left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right]\quad (\left[\begin{array}{cc} v_\lambda\\ u_\lambda \end{array}\right]\ne0) }

を満足する\lambdaとして定義されます。これから次式が成り立ちます。

\displaystyle{(4)\quad \left\{\begin{array}{l} v_\lambda=(\lambda I_n-A)^{-1}Bu_\lambda\\ Cv_\lambda=0 \end{array}\right. \Rightarrow C(\lambda I_n-A)^{-1}Bu_\lambda=0 }

詳しい説明は省きますが、これは指数関数入力e^{\lambda t}u_\lambdaに対する応答が零になることを意味しています。

(3)が非自明な解をもつためには次式が必要十分となります。

\displaystyle{(5)\quad {\rm det}\left[\begin{array}{cc} \lambda I_n-A & -B\\ -C & 0 \end{array}\right]=0 }

さて、次の線形系を考えます。

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

この不変零点を求めます。

\displaystyle{(7)\quad \begin{array}{l} {\rm det}\left[\begin{array}{cc} \lambda I_n-A & -B\\ -S & 0 \end{array}\right]\\ = {\rm det}\left[\begin{array}{ccc} \lambda I_{n-m}-A_{11} & -A_{12} & 0\\ -A_{21} & \lambda I_m-A_{22} & -B_{2}\\ -S_1 & -S_2 & 0 \end{array}\right]=0\\ \Leftrightarrow {\rm det}\left[\begin{array}{ccc} \lambda I_{n-m}-A_{11} & -A_{12} \\ -S_1 & -S_2 \end{array}\right]=0\\ \Leftrightarrow{\rm det} \left[\begin{array}{cc} I_{n-m} & -A_{12}S_2^{-1} \\ 0 & I_m \end{array}\right]\left[\begin{array}{cc} \lambda I_{n-m}-\bar{A}_{11} & 0 \\ 0 & -S_2 \end{array}\right] \left[\begin{array}{cc} I_{n-m} & 0 \\ M & I_m \end{array}\right]\\ ={\rm det}(\lambda I_{n-m}-\bar{A}_{11}){\rm det}(-S_2)=0\\ \Leftrightarrow {\rm det}(\lambda I_{n-m}-\bar{A}_{11})=0 \end{array} }

すなわち、スイッチング関数を出力方程式とするシステム(6)の不変零点は、行列\bar{A}_{11}の固有値と一致します。この行列\bar{A}_{11}はスライディングモードの振舞いを決めるので、\bar{A}_{11}は安定行列とすることが前提となります。したがって、どのようなスイッチング関数が望ましいかは、どのような超平面上をどのようにスライディングさせるかに関わっていると言えます。

問題設定

投稿日時: 投稿者:

問題設定…Homework

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

\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, 1\le 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)ことのできるスライディングモード制御則（SMC則）を設計せよ。

この問題の解、SMC則は、次のように表されます。

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

\displaystyle{(5)\quad { \boxed{u_\ell(t)=-\underbrace{(SB)^{-1}(SA-\Phi S)}_{L=L_{eq}+L_{phi}}x(t)}} }

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

ここで、SBは正則行列であることを仮定しています。また設計パラメータは、安定行列\Phi、スカラー\rho>0です。P_2は次式から定まる正定行列です。

\displaystyle{(7)\quad P_2\Phi+\Phi^TP_2=-I_m }

以下ではこの解の妥当性を２次安定化の観点から検討します。

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

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

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

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

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

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

ちなみに、(4)の第１項u_\ellは、この等価制御をベースして

\displaystyle{(11)\quad u_\ell(t)=-\underbrace{(SB)^{-1}SA}_{L_{eq}}x(t)-\underbrace{(SB)^{-1}(-\Phi) S}_{L_{phi}}x(t) }

のように構成されています。この\Phiの役割についてはあとで述べます。

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

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

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

\displaystyle{(13)\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{(14)\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} }

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


演習…Flipped Classroom
1^\circ 次の剛体振子の状態方程式

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

に対して、スイッチング関数を

\displaystyle{ s(t)= \underbrace{ \left[\begin{array}{cc} m & 1 \end{array}\right] }_{S} \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }

と選ぶとき、SMC則の線形制御の部分は次式となることを示せ。

\displaystyle{ u_\ell(t)=- \left[\begin{array}{cc} -m\Phi & m-\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] }

※-\Phiを改めて\Phiとおけば、可到達条件のところで述べた次式(14)となります。

\displaystyle{(14)\quad u(t)=- \left[\begin{array}{cc} m\Phi & m+\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\frac{s(t)}{|s(t)|} }

2^\circ 上の剛体振子の状態方程式についてマッチング条件(12)を調べよ。

Note C21１入力系の場合

●ここでは次の１入力系を考えます。

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

\displaystyle{(2)\quad A= \left[\begin{array}{cccc|c} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & 0 & \vdots \\ \vdots & & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 1 \\\hline -a_1 & -a_2 & \cdots & \cdots & -a_n \end{array}\right],\quad B=\left[\begin{array}{c} 0 \\ 0 \\ \vdots \\ 0 \\\hline 1 \end{array}\right] }

\displaystyle{(3)\quad \begin{array}{r} \Delta A(t)= \left[\begin{array}{cccc|c} 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 0 & \vdots \\ \vdots & & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 0 \\\hline -\delta_1(t) & -\delta_2(t) & \cdots & \cdots & -\delta_n(t) \end{array}\right]\\ k_i^-\le\delta_i(t)\le k_i^+\quad(i=1,\cdots,n) \end{array} }

また次のスイッチング関数を考えます。

\displaystyle{(4)\quad s(t)=Sx(t) }

\displaystyle{(5)\quad S=\left[\begin{array}{ccccc} m_1 & m_2 & \cdots & m_{n-1} & 1 \end{array}\right] }

このとき、SMC則

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

を可到達条件を満たすように構成することを考えます。

●線形制御u_\ellは、

(7)\quad \begin{array}{l} \displaystyle{\dot{s}(t)=\sum_{i=1}^{n-1}m_i\dot{x}_i(t)+\dot{x}_n(t)}\\ \displaystyle{=\sum_{i=1}^{n-1}m_i x_{i+1}(t)-\sum_{i=1}^{n}(a_i+\delta_i(t))x_i(t)+u(t)}\\ \displaystyle{=-a_1x_1(t)+\sum_{i=2}^{n}(m_{i-1}-a_i) x_i(t)- \sum_{i=1}^{n}\delta_i(t)x_i(t)+u(t)} \end{array}

において、\delta_i=0とおいて

\displaystyle{(8)\quad \boxed{u_\ell(t)=a_1x_1(t)+\sum_{i=2}^{n}(a_i-m_{i-1}) x_i(t)} }

のように定めます。

●u_nの候補として次の２つを考えます。

\displaystyle{(9)\quad \boxed{u_n(t)=-\rho(t,x)\frac{s(t)}{|s(t)|}} }

\displaystyle{(10)\quad \boxed{u_n(t)=\sum_{i=1}^{n}k_ix_i(t)-\eta\frac{s(t)}{|s(t)|}} }

以下では、これらが適当な条件の下で可到達条件s\dot{s}<0を満足しSMC則になることを示します。

(9)の場合：

\displaystyle{(11)\quad \begin{array}{l} \rho(t,x)=\sum_{i=1}^{n}\max\{|k_i^+|,|k_i^-|\}|x_i(t)|+\eta\\ \ge |\sum_{i=1}^{n}\delta_i(t)x_i(t)|+\eta \end{array} }

を満足する\rho(t,x)を選ぶと

\displaystyle{(12)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+u_n(t))\\ =s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)-\rho(t,x)\frac{s(t)}{|s(t)|})\\ =-s(t)\sum_{i=1}^{n}\delta_i(t)x_i(t)-|s(t)|\rho(t,x)\\ \le |s(t)|(|\sum_{i=1}^{n}\delta_i(t)x_i(t)|-\rho(t,x))\\ < -\eta|s(t)| \end{array} }

(10)の場合：

\displaystyle{(13)\quad k_i= \left\{\begin{array}{ll} k_i^- & (sx_i>0)\\ k_i^+ & (sx_i<0) \end{array}\right. }

のようにk_iを選ぶと

\displaystyle{(14)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+u_n(t))\\ =s(t)(-\sum_{i=1}^{n}\delta_i(t)x_i(t)+\sum_{i=1}^{n}k_ix_i(t)-\eta\frac{s(t)}{|s(t)|})\\ =\sum_{i=1}^{n}\underbrace{(k_i-\delta_i(t))s(t)x_i(t)}_{<0}-\eta|s(t)|\\ < -\eta|s(t)| \end{array} }

ロバスト性

投稿日時: 投稿者:

線形系の２次安定性…Homework

対象が平衡状態にあることは線形状態方程式\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

条件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 \frac{d}{dt}V(x)=2x^T(t)P\dot{x}(t)=2x^T(t)PAx(t)=-x^T(t)Qx(t)\quad (Q=I) }

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

●条件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{Q^{1/2}PQ^{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

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

モデルの不確かさがある場合…Homework

いま漸近安定な線形系(1)において次のようなモデルの不確かさがあるとします。

\displaystyle{(7)\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{(8)\quad PA+A^TP+Q=0 }

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

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

ならば、リャプノフ関数

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

に対して

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

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

●実際、

\displaystyle{(12)\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{(13)\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} }

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

\displaystyle{(14)\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)|| \\ =(||\xi(t,x)||-\frac{1}{2}\frac{\sigma_n(Q)}{\sigma_1(P)}||x(t)||)2\sigma_1(P)||x(t)||<0 \end{array} }

となります。

●次の数値例を考えてみます。

\displaystyle{(15)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -4 & -1 \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{\left[\begin{array}{c} 0 \\ 1 \end{array}\right]}_{D}\underbrace{\sin(2t)}_{\xi(t)} }

Aは安定行列なので、リャプノフ方程式PA+A^TP+I_n=0を満足するP>0は

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

のように求まります。

このPをもつリャプノフ関数(10)に対して

\displaystyle{(17)\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)+D\xi(t))+(x^T(t)A^T+\xi^T(t)D^T)Px(t)\\ =x^T(t)PAx(t)+x^T(t)A^TPx(t)+x^T(t)PD\xi(t)+\xi(t)D^TPx(t)\\ =x^T(t)(PA+A^TP)x(t)+2\xi(t)D^TPx(t)\\ =-x^T(t)Qx(t)+2\xi(t)D^TPx(t)\\ \le -x^T(t)x(t)+ 2|\xi(t)|||D^TP||||x(t)||\\ \le -||x(t)||^2+ 2\sigma_1(PD)||x(t)|| \\ =(2\sigma_1(PD)-||x(t)||)||x(t)||<0 \end{array} }

これから

\displaystyle{(18)\quad ||x(t)||>2\sigma_1(PD)=1.2748 }

のときは\frac{d}{dt}V(x)<0となって原点に漸近していきますが、||x(t)||<1.2748になるとそれは保証されないことが分かります(図1参照)。

図1 モデルの不確かさがある場合の解軌道

演習C12…Flipped Classroom
1^\circ　smc17.slxを用いて、図1のシミュレーションを行え。

Note C12 振り子の安定性

振り子の運動方程式

\displaystyle{(1)\quad \underbrace{\frac{4}{3}m\ell^2}_J \dot{\omega}(t)=-mg\ell\sin\theta(t)-c\,\omega(t)\quad(\omega(t)=\dot{\theta}(t)) }

から得られる非線形状態方程式

\displaystyle{(2)\quad \frac{d}{dt} \underbrace{\left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]}_{x} = \underbrace{ \left[\begin{array}{c} \omega(t) \\ -\frac{mg\ell}{J}\sin\theta(t)-\frac{c}{J}\omega(t) \end{array}\right] }_{f(x)} }

に対して、次のリャプノフ関数の候補を考えます。

\displaystyle{(3)\quad V=\frac{1}{2} \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right]^T \left[\begin{array}{cc} \frac{1}{2}\frac{c^2}{J} & \frac{1}{2}c \\ \frac{1}{2}c & J \end{array}\right] \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] }_{x^T(t)Px(t)} +mg\ell(1-\cos\theta(t))>0 }

この時間微分を計算すると

\displaystyle{(4)\quad \begin{array}{l} \dot{V}(t)=\frac{1}{2}(x^TP\dot{x}+\dot{x}^TPx)+mg\ell\omega\sin\theta =x^TP\dot{x}+mg\ell\omega\sin\theta\\ =\left[\begin{array}{c} \theta \\ \omega \end{array}\right]^T \left[\begin{array}{cc} p_1 & p_3 \\ p_3 & p_2 \end{array}\right] \left[\begin{array}{c} \omega \\ -\frac{mg\ell}{J}\sin\theta-\frac{c}{J}\omega \end{array}\right]+mg\ell\omega\sin\theta\\ =\left[\begin{array}{cc} \theta &\omega \end{array}\right] \left[\begin{array}{cc} p_1\omega -p_3(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \\ p_3\omega -p_2(\frac{mg\ell}{J}\sin\theta+\frac{c}{J}\omega) \end{array}\right]+mg\ell\omega\sin\theta\\ =p_1\theta\omega -p_3(\frac{mg\ell}{J}\theta\sin\theta+\frac{c}{J}\theta\omega)+ p_3\omega^2 -p_2(\frac{mg\ell}{J}\omega\sin\theta+\frac{c}{J}\omega^2)+mg\ell\omega\sin\theta\\ =-\frac{mg\ell}{J}\sin\theta(p_3\theta+p_2\omega-J\omega) +(p_3-p_2\frac{c}{J})\omega^2+(p_1-p_3\frac{c}{J})\theta\omega \end{array} }

ここで、p_1=\frac{1}{2}\frac{c^2}{J}、p_2=J、p_3=\frac{1}{2}cを代入して

\displaystyle{(5)\quad \dot{V}(t)=-\frac{1}{2}\frac{mg\ell}{J}c\,\theta(t)\sin\theta(t) -\frac{1}{2}c\,\omega^2(t) }

となります。ここで

\displaystyle{(6)\quad \theta(t)\sin\theta(t)>0\quad (|\theta(t)|<\pi,\theta(t)\ne0) }

だから、(5)は負値をとります。したがって、(3)はリャプノフ関数となっています。

●ちなみに(2)を次のように書き直してみます。

\displaystyle{(7)\quad \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{mg\ell}{J} & -\frac{c}{J} \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ \frac{mg\ell}{J} \end{array}\right] }_{D} \underbrace{ (\theta(t)-\sin\theta(t)) }_{\xi(x)} }

我々は軸摩擦をもつ振り子の角度は十分時間がたてば零に収束することを知っています。上の議論に沿って解釈すれば、モデルの不確かさ\xi(x)自体が限りなく零に近づき（\sin\theta\simeq\theta）、線形系となると考えられます。

Note C12 　MATLABによるリャプノフ方程式の求解

たとえば、安定行列\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で次のコマンドを用いました。

P=lyap(A’,eye(2))

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

Loading

２次安定性

投稿日時: 投稿者:

リャプノフの安定定理

非線形システム理論に関する文献：
Hassan K.Khalil: Nonlinear Systems, 3rd ed., Prentice Hall, 2002
から次の定理（Theorem 4.1, p.114）を引用します。

リャプノフの安定定理　状態方程式

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

に対して、その平衡状態 x^*=0 を含む領域 {\cal D} を考える。いま {\cal D} で連続微分可能な実数値関数 V に対して

\displaystyle{(2)\quad \left\{\begin{array}{l} V(0)=0 \\ V(x)>0\quad(x\in{\cal D},x\not=0) \end{array}\right.}

このとき

\displaystyle{(3)\quad \frac{d}{dt}V(x)\le0 \quad(x\in{\cal D})}

ならば、平衡状態 x^*=0 は安定、また

\displaystyle{(4)\quad \frac{d}{dt}V(x)<0 \quad(x\in{\cal D},x\not=0)}

ならば、平衡状態 x^*=0 は漸近安定である。

すなわち、(2)と(4)を満足するリャプノフ関数Vの存在を示して（発見して）、非線形系(1)の漸近安定性を示すことができます。

２次安定性

非線形系(1)に対してリャプノフ関数

\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 \boxed{\alpha=\frac{\sigma_n(Q)}{\sigma_1(P)},\ \beta=\sqrt{\frac{\sigma_1(P)}{\sigma_n(P)}} %\alpha=\frac{\sigma_n(Q)}{\sigma_1(P)}\le\sigma_n(P^{-1}Q),\ }

以下では、これを導出します。

●予備知識として、任意の対称行列R>0に対して、特異値（固有値に等しい）を

\displaystyle{(10)\quad \sigma_1(R)\ge\cdots\ge\sigma_n(R) }

と表記し、また任意のv\ne0に対して

\displaystyle{(11)\quad \sigma_n(R)\le\frac{v^TRv}{v^Tv}\le\sigma_1(R) }

および

\displaystyle{(12)\quad \sigma_1(R^{-1})=\frac{1}{\sigma_n(R)} }

が成り立つこと使います。

●そこで

\displaystyle{(13)\quad \begin{array}{l} \sigma_n(Q)x^T(t)x(t)\le x^T(t)Qx(t)\\ x^T(t)Px(t)\le\sigma_1(P)x(t)^Tx(t) \end{array} }

に注意して、次式を得ます。

\displaystyle{(14)\quad \frac{d}{dt}V(x)\le -x^T(t)Qx(t) \le -\sigma_n(Q)x^T(t)x(t) \le -\sigma_n(Q) \frac{x^T(t)Px(t)}{\sigma_1(P)} \le -\alpha V(x) }

これから

\displaystyle{(15)\quad V(x)\le V(0)e^{-\alpha t} \Leftrightarrow x^T(t)Px(t)\le x^T(0)Px(0)e^{-\alpha t} }

ここで

\displaystyle{(16)\quad {\begin{array}{l} %\sigma_n^2(P)x^Tx\le x^TPx\le\sigma_1^2(P)x^Tx\\ %\sigma_n^2(P)x^T(0)x(0)\le x^T(0)Px(0)\le\sigma_1^2(P)x(0)^Tx(0) \sigma_n(P)x^T(t)x(t)\le x^T(t)Px(t)\\ x^T(0)Px(0)\le\sigma_1(P)x(0)^Tx(0) \end{array} }}

に注意して

\displaystyle{(17)\quad {\sigma_n(P)x^T(t)x(t) \le \sigma_1(P) x^T(0)x(0)e^{-\alpha t} }}

すなわち

\displaystyle{(18)\quad {\underbrace{x^T(t)x(t)}_{||x(t)||^2} \le \underbrace{\frac{\sigma_1(P)}{\sigma_n(P)}}_{\beta^2} \underbrace{x(0)^Tx(0)}_{||x(0)||^2} e^{-\alpha t} }}

を得ます。これから(8)が成り立ちます。なお一般には

\displaystyle{(19)\quad \underbrace{\frac{\sigma_n(Q)}{\sigma_1(P)}}_{\alpha} =\frac{1}{\sigma_1(Q^{-1})\sigma_1(P)} \le \frac{1}{\sigma_1(Q^{-1}P)} =\underbrace{\sigma_n(P^{-1}Q)}_{\alpha'} }

ですが、PとQを対角行列に選ぶときは等号が成り立ち、(8)において\alphaではなく\alpha'を用いることができます。もう少し詳しくは、対称行列X=PとY=Q^{-1}が可換（XY=YX）で、同時対角化可能（TXT^{-1}とTYT^{-1}を対角行列とするTが存在）の場合には等号が成り立ち、(8)において\alphaではなく\alpha'を用いることができます。

Note C11 (11),(12)の導出

●任意の正定行列R>0（サイズn\times n）に対して、次の特異値分解を考えます。

\displaystyle{(1)\quad \begin{array}{l} R=U\Sigma U^T\\ \Sigma={\rm diag}\{\sigma_1,\cdots,\sigma_n\}\quad(\sigma_1\ge\cdots\ge\sigma_n)\\ UU^T=U^TU=I_n \end{array} }

ここで、Rは対称行列なので固有値は実数、しかも正定なのでそれらはすべて正であること、したがって特異値\sigma_1はRの最大固有値、特異値\sigma_nはRの最小固有値であることに注意してください。このとき、R=R^{1/2}R^{1/2}を満たすRの平方根行列は

\displaystyle{(2)\quad \begin{array}{l} R^{1/2}=U\Sigma^{1/2} U^T\\ \Sigma^{1/2}={\rm diag}\{\sqrt{\sigma_1},\cdots,\sqrt{\sigma_n}\} \end{array} }

で与えられます。また、Rの逆行列は

\displaystyle{(3)\quad \begin{array}{l} R^{-1}=U\Sigma^{-1} U^T\\ \Sigma^{-1}={\rm diag}\{\frac{1}{\sigma_1},\cdots,\frac{1}{\sigma_n}\} \quad(\frac{1}{\sigma_n}\ge\cdots\ge\frac{1}{\sigma_1}) \end{array} }

のように表されます。これから(12)は自明です。

いま、R^{1/2}を用いて線形写像y=R^{1/2}xを考えると、これは３つの線形写像に分解されます。

\displaystyle{(4)\quad \begin{array}{l} x'=U^Tx\\ y'=\Sigma^{1/2} x'\\ y=Uy' \end{array} }

任意のn次元ベクトルx\in{\rm\bf R}^nのノルムとして、次の２ノルムを考えます。

\displaystyle{(5)\quad ||x||_2=\sqrt{x_1^2+\cdots+x_n^2}&\Rightarrow&||x||_2^2=x^Tx\\ }

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

\displaystyle{(6)\quad ||x'||_2^2=||U^Tx||_2^2=x^TUU^Tx=x^Tx=||x||_2^2\ \Rightarrow\ \frac{||x'||_2}{||x||_2}=1 }

\displaystyle{(7)\quad \begin{array}{l} ||y'||_2^2=||\Sigma^{1/2} x'||_2^2=\sigma_1x_1'\,^2+\cdots+\sigma_nx_n'\,^2\\ \le \sigma_1(x_1'\,^2+\cdots+x_n'\,^2)=\sigma_1x'\,^Tx'=\sigma_1||x'||_2^2 \ \Rightarrow\ \frac{||y'||_2}{||x'||_2}\le\sqrt{\sigma_1} \end{array} }

\displaystyle{(8)\quad ||y||_2^2=||Uy'||_2^2=y'\,^TU^TUy'=y'\,^Ty'=||y'||_2^2\ \Rightarrow\ \frac{||y||_2}{||y'||_2}=1 }

すなわち

\displaystyle{(9)\quad \frac{||x'||_2}{||x||_2}\times\frac{||y'||_2}{||x'||_2}\times\frac{||y||_2}{||y'||_2}\le\sqrt{\sigma_1} \ \Rightarrow\ \frac{||y||_2}{||x||_2}=\frac{||R^{1/2}x||_2}{||x||_2}\le\sqrt{\sigma_1} }

一方

\displaystyle{(10)\quad \begin{array}{l} ||y'||_2^2=||\Sigma^{1/2} x'||_2^2=\sigma_1x_1'\,^2+\cdots+\sigma_nx_n'\,^2\\ \ge \sigma_n(x_1'\,^2+\cdots+x_n'\,^2)=\sigma_nx'\,^Tx'=\sigma_n||x'||_2^2 \ \Rightarrow\ \frac{||y'||_2}{||x'||_2}\ge\sqrt{\sigma_n} \end{array} }

より、次式が成り立ちます。

\displaystyle{(11)\quad \frac{||x'||_2}{||x||_2}\times\frac{||y'||_2}{||x'||_2}\times\frac{||y||_2}{||y'||_2}\ge\sqrt{\sigma_n} \ \Rightarrow\ \frac{||y||_2}{||x||_2}=\frac{||R^{1/2}x||_2}{||x||_2}\ge\sqrt{\sigma_n} }

したがって、(7)と合わせて、次式を得ます。

\displaystyle{(12)\quad \sqrt{\sigma_n} \le \frac{||R^{1/2}x||_2}{||x||_2}=\frac{\sqrt{x^TR^{1/2}R^{1/2}x}}{\sqrt{x^Tx}} \le\sqrt{\sigma_1} }

これを２乗して次が得られます。

\displaystyle{(13)\quad \sigma_n \le \frac{x^TRx}{x^Tx} \le\sigma_1 }

すなわち(11)が導出されました。

可到達条件

投稿日時: 投稿者:

可到達条件…Homework

[0] 剛体振り子の状態方程式として次式を考えます。

\displaystyle{(1)\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] }

これに対して次のSMCを考えていました。

\displaystyle{(2)\quad u(t)=-{\rm sgn}(s(t)) %-\frac{s(t)}{|s(t)|} }

\displaystyle{(3)\quad s(t)=mx_1(t)+x_2(t) }

振り子を離す初期角度が小さければ容易にスライディング直線に拘束されますが、大きな初期角度をつけて話すと一旦スライディング直線に到達してもこれを突き抜けてしまう可能性があります（前節の図１と図２に示す状況です）。そこで可到達条件

\displaystyle{(4)\quad %\begin{array}{l} \boxed{s(t)\dot{s}(t)< 0} %\end{array} }

をどう満足させるかについて検討します。

[1]　制御則(2)に対して、減衰力を与えるために状態フィードバックを加えてみます。

\displaystyle{(5)\quad u(t)=- \left[\begin{array}{cc} f_1 & f_2 \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} }

ここで、\rhoというパラメータも導入されており、高速スイッチングする操作入力の振幅を表しています。このとき、次式を得ます。

\displaystyle{(6)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(m\dot{x}_1(t)+\dot{x}_2(t))\\ =s(t)(mx_2(t)-a\sin x_1(t)+u(t))\\ =s(t)(mx_2(t)-a\sin x_1(t)-f_1x_1(t)-f_2x_2(t)-\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} ) \end{array} }

ここで、f_1=0, f_2=mと選ぶと、-s(t)a\sin x_1(t)\le|s(t)|aに注意して

\displaystyle{(7)\quad \begin{array}{l} s(t)\dot{s}(t) =s(t)(-a\sin x_1(t)-\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} )\\ =-s(t)a\sin x_1(t)-\rho|s(t)|\\ <|s(t)|(a-\rho) \end{array} }

したがって、ある定数\eta>0に対して

\displaystyle{(8)\quad a-\rho < -\eta \quad \Leftrightarrow \quad \boxed{\rho > a+\eta} }

を満たす\rhoを選べば

\displaystyle{(9)\quad \boxed{s(t)\dot{s}(t)< -\eta|s(t)|} }

が成り立ちます。

結局、制御則

\displaystyle{(10)\quad { u(t)=- \left[\begin{array}{cc} 0 & m \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} }}

を用いれば、スライディング直線に一旦到達すれば、以後は留まるようにできます。


図1 x_1(0)=3, x_2(0)=0のときの状態軌道

図2 x_1(0)=3, x_2(0)=0のとき、m=1,\rho=1の制御則(5)の下での状態軌道

●その到達時刻を見積もってみます。

\displaystyle{(11)\quad s(t)\dot{s}(t)< -\eta|s(t)| \quad\Leftrightarrow\quad \left\{\begin{array}{ll} \dot{s}(t)< -\eta & (s(t)>0)\\ \dot{s}(t)> -\eta & (s(t)<0) \end{array} \right. }

から

\displaystyle{(12)\quad \frac{1}{2}\frac{d}{dt}s^2(t)< -\eta|s(t)| \Rightarrow \frac{1}{2}\frac{d}{dt}|s(t)|^2=|s(t)|\frac{d}{dt}|s(t)|< -\eta|s(t)| \Rightarrow \frac{d}{dt}|s(t)|< -\eta }

両辺を積分すると

\displaystyle{(13)\quad \underbrace{|s(t_s)|}_0-|s(0)|< -\eta t_s \quad\Leftrightarrow\quad \boxed{t_s\le \frac{|s(0)|}{\eta}} }

を得て、スライディング直線に到達する時刻t_sを見積もることができます。

高速スイッチングする操作入力の振幅\rhoは小さい方が望ましいといえます。ただ、(8)より\etaも小さくせざる得ず、このことはスライディング直線に到達する速さに影響を与えます。

[2]　もう一つ別のパラメータ\Phi>0を導入して、制御則(10)を次のように変更してみます。

\displaystyle{(14)\quad { u(t)=- \left[\begin{array}{cc} m\Phi & m+\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} }}

このとき、(3)と(8)に注意して、次式を得ます。

\displaystyle{(15)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(m\dot{x}_1(t)+\dot{x}_2(t))\\ =s(t)(mx_2(t)-a\sin x_1(t)+u(t))\\ =s(t)(mx_2(t)-a\sin x_1(t)-m\Phi x_1(t)-(m+\Phi)x_2(t)-\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} )\\ =-\Phi s(t)(m x_1(t)+x_2(t))-s(t)(a\sin x_1(t)+\rho\,{\rm sgn}(s(t)) %\frac{s(t)}{|s(t)|} )\\ \le -\Phi s^2(t)+|s(t)|(a-\rho)\\ < -\Phi s^2(t)-\eta |s(t)|\\ \end{array} }

これより

\displaystyle{(16)\quad \boxed{s(t)\dot{s}(t)< -\Phi s^2(t)} \quad \Leftrightarrow \quad \left\{\begin{array}{ll} \dot{s}(t)< -\Phi s(t) & (s(t)>0)\\ \dot{s}(t)> -\Phi s(t) & (s(t)<0) \end{array}\right. }

これは\Phiが|s(t)|\rightarrow 0となる速さを決めていることを意味しています。


図3 x_1(0)=3, x_2(0)=0のとき、m=1,\Phi=1,\rho=0.3の制御則(14)の下での状態軌道

●ちなみに、計算の都合上、(14)の代わりにつぎのスイッチング方式が用いられます。

\displaystyle{(17)\quad { u(t)=- \left[\begin{array}{cc} m\Phi & m+\Phi \end{array}\right] \left[\begin{array}{cc} x_1(t)\\ x_2(t) \end{array}\right] -\rho\frac{s(t)}{|s(t)|+\delta}} }

ここで\delta>0は適切に選ばれた定数です。

演習C02…Flipped Classroom
1^\circ　smc15.slxを用いて、図1、図2のシミュレーションを行え。
2^\circ　smc16.slxを用いて、図3のシミュレーションを行え。

SM制御

投稿日時: 投稿者:

SM制御…Homework

[1] 次の無定位系(２重積分器)を考えます。

\displaystyle{(1)\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) }

これに対して、スイッチング関数と呼ばれる

\displaystyle{(2)\quad s(t)=mx_1(t)+x_2(t) }

を用いて、制御則

\displaystyle{(3)\quad u(t)=\left\{\begin{array}{ll} 1 & (s(t)<0) \\ 0 & (s(t)=0) \\ -1 & (s(t)>0) \end{array}\right. }

すなわち

\displaystyle{(4)\quad \boxed{u(t)=-{\rm sgn}(s(t)) %-\frac{s(t)}{|s(t)|}} }

を考えます。\dot{x}(t)=Ax(t)+B(\pm 1)の解は(\tau'=t-\tau\Rightarrow d\tau=-d\tau'に注意して)

\displaystyle{(5)\quad x(t)=\pm\int_0^t\exp(A(t-\tau))Bd\tau =\mp\int_t^0\exp(A\tau')Bd\tau' =\pm\int_0^t\exp(A\tau)Bd\tau }

また

\displaystyle{(6)\quad \exp( \left[\begin{array}{cc} 0 & 1\\ 0 & 0 \end{array}\right]\tau) = \left[\begin{array}{cc} 1 & \tau \\ 0 & 1 \end{array}\right] \Rightarrow \exp(A\tau)B= \left[\begin{array}{cc} 1 & \tau \\ 0 & 1 \end{array}\right] \left[\begin{array}{c} 0\\ 1 \end{array}\right] = \left[\begin{array}{c} \tau \\ 1 \end{array}\right] }

となることから

\displaystyle{(7)\quad \left\{\begin{array}{l} x_1(t)=\pm\int_0^t \tau d\tau=\pm\frac{1}{2}t^2+x_1(0)\\ x_2(t)=\pm\int_0^t 1 d\tau=\pm t+x_2(0) \end{array}\right. \Rightarrow x_1(t)-x_1(0)=\pm \frac{1}{2}(x_2(t)-x_2(0))^2 }

と計算できます。ここで、状態軌道は放物線となっていることに注意します。

●m=2、x_1(0)=x_2(0)=4のときの状態軌道を図1に示します。


図1 x_1(0)=x_2(0)=4のときの状態軌道

これから直線x_2=-mx_1の上側と下側にある（横に寝た）放物線を切り替えて状態軌道が構成されていることが分かります。ところが、原点に近くなると奇妙な現象が生じます。

●m=2、x_1(0)=x_2(0)=0.4のときの状態軌道を図2に示します。


図2 x_1(0)=x_2(0)=0.4のときの状態軌道

これから状態軌道が直線x_2=-mx_1に到達すると、この直線上に拘束されて原点まで滑っていることが分かります。一体、何が起こっているのでしょう？

●これを調べるために

\displaystyle{(9)\quad \begin{array}{l} s(t)\dot{s}(t)=s(t)(m\dot{x}_1(t)+\dot{x}_2(t))=s(t)(mx_2(t)+u(t))\\ =s(t)(mx_2(t)-\frac{s(t)}{|s(t)|})=mx_2(t)s(t)-|s(t)|\\ \Rightarrow s(t)\dot{s}(t)\le (m|x_2(t)|-1)|s(t)| \end{array} }

から、次が成り立つことがわかります。

\displaystyle{(10)\quad \begin{array}{l} m|x_2(t)|<1 \Rightarrow s(t)\dot{s}(t)< 0 \end{array} }

いま(2)で定義したsはx_1-x_2平面におけるx_2軸の切片であることに注意します。(10)より、s>0のとき\dot{s}<0なので切片sはx_2軸正の方から0に漸近します。一方、s<0のとき\dot{s}>0なので切片sはx_2軸負の方から0に漸近します。このことは状態は直線

\displaystyle{(11)\quad mx_1(t)+x_2(t)=0 }

の上に拘束されてしまうことを意味します。これから

\displaystyle{(12)\quad \dot{x}_1(t)=-mx_1(t) }

となるので、つぎを得ます。

\displaystyle{(13)\quad x_1(t)=e^{-mt}x_1(0)\rightarrow 0\quad(t\rightarrow\infty) }

これは直線(2)を滑るように原点に向かうことを意味します。

●以上のように、適当なスイッチング関数s(t)に対して、可到達条件

\displaystyle{(14)\quad　{ \begin{array}{l} s(t)\dot{s}(t)< 0 \end{array}} }

を満足させ、スライディング直線に乗せたあと、原点に向かって滑らせる制御方式をスライディングモード制御と呼びます。

[2]　いまスライディング直線に乗る時刻をt_sとすると、次が成り立ちます。

\displaystyle{(15)\quad s(t)=0\quad(t>t_s) \quad\Rightarrow\quad \dot{s}(t)=0\quad(t>t_s) }

また(9)では次式を用いました。

\displaystyle{(16)\quad \dot{s}(t)=m\dot{x}_1(t)+\dot{x}_2(t)=mx_2(t)+u(t) }

そうすると(15)と(16)を合わせて、スライディング直線を滑らせる制御則は

\displaystyle{(17)\quad 0=mx_2(t)+u_{eq}(t)\Rightarrow u_{eq}(t)=-mx_2(t)\quad(t>t_s) }

のように考えられ、等価制御と呼ばれています。

●等価制御はスライディングモード制御そのものではないのですが、どのように違うか、次のシミュレーション例を見てください。これからスライディングモード制御則を１次システム1/(0.04s+1)を通したものが等価制御とよく合っていることがわかります。この１次システムはアクチュエータのダイナミックス（による帯域制限）を表しているとみなすことができます。


図3 SMCとその等価制御（m=1, x_1(0)=1, x_2(0)=0）

そうであれば最初から等価制御を用いればよいではないかと思うかもしれませんが、等価制御はあくまでスライディングモードに突入後の制御則です。したがって、まずは可到達条件(14)を満足させて、スライディングモードを実現する必要があることに注意してください。

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

\displaystyle{(18)\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{(19)\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{(20)\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] }

そこで、これまでのスライディングモード制御(4)をそのまま用いたときの様子を調べてみます。

次は仮にa=0.25とした場合のシミュレーション結果です。一見、図3と変わらないように見えます。


図4 振り子のSMCとその等価制御（m=1, x_1(0)=1, x_2(0)=0）

ここには、次の等価制御もシミュレーションされています。

\displaystyle{(21)\quad 0=mx_2(t)-a\sin x_1(t)+u_{eq}(t)\Rightarrow u_{eq}(t)=a\sin x_1(t)-mx_2(t)\quad(t>t_s) }

実は、無定位系そのものは非線形項を考慮していないのに、これに非線形項が加わった場合の影響を抑制していることが分かります。この事実がスライディングモード制御を大変魅力的なものとしています。

演習C01…Flipped Classroom
1^\circ　smc11.slxとsmc12.slxを用いて、図1、図2のシミュレーションを行え。
2^\circ　smc13.slxとsmc14.slxを用いて、図3、図4のシミュレーションを行え。

柔軟ビーム(続)

投稿日時: 投稿者:

[1] 柔軟構造物の制振問題の例として

のような柔軟ビームを、一定角度だけ振動を抑制しながら回転させる制御系設計を考えていきます。そのために柔軟ビームに関する運動方程式

\displaystyle{(1a)\quad \underbrace{ \left[\begin{array}{c|ccc} \alpha & m_1 & \dots &m_N \\ \hline m_1 & 1 & & 0 \\ \vdots & & \ddots & \\ m_N & 0 & & 1 \\ \end{array}\right] }_{M=\left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right]} \underbrace{ \left[\begin{array}{c} \ddot{\theta} \\ \hline \ddot{r}_1 \\ \vdots \\ \ddot{r}_N \end{array}\right] }_{\ddot{\xi}} + \underbrace{ \left[\begin{array}{c|ccc} 0 & 0 & \dots &0 \\ \hline 0 & \Omega_1^2 & & 0 \\ \vdots & & \ddots & \\ 0 & 0 & & \Omega_N^2 \\ \end{array}\right] }_{K=\left[\begin{array}{cc} K_{11} & K_{12} \\ K_{21} & K_{22} \end{array}\right]} \underbrace{ \left[\begin{array}{c} \theta \\ \hline r_1 \\ \vdots \\ r_N \end{array}\right] }_{\xi} = \underbrace{ \left[\begin{array}{c} 1 \\ \hline 0 \\ \vdots \\ 0 \end{array}\right] }_{B_2} u' }
\displaystyle{(1b)\quad \eta(\xi)= \underbrace{ \left[\begin{array}{c|ccc} 0 & \phi_1(\xi) & \dots &\phi_N(\xi) \end{array}\right] }_{C_1} \underbrace{ \left[\begin{array}{c} \theta \\ \hline r_1 \\ \vdots \\ r_N \end{array}\right] }_{\xi} }

から次の状態空間表現を得ていました。

\displaystyle{(2a)\quad \underbrace{ \left[\begin{array}{c} \dot{\xi} \\ \ddot{\xi} \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cc} 0_{N+1\times N+1} & I_{N+1} \\ M^{-1}K & 0_{N+1\times N+1} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \xi \\ \dot{\xi} \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{cc} 0_{N+1\times 1} \\ M^{-1}B_2 \end{array}\right] }_{B} u'}
\displaystyle{(2b)\quad \eta(\xi)= \underbrace{ \left[\begin{array}{cc} C_1 & 0_{1\times N} \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \xi \\ \dot{\xi} \end{array}\right] }_{x} }

このモデルでは、柔軟ビームを回転させるモータのトルクの時系列を与えたときの弾性変位の計算はできます。したがって、トルク制御を行うことが考えられますが、一般にモータの購入時には、回転角制御または回転数制御のどちらかの目的に応じて、ドライバも一緒に購入します。ここでは、回転数制御用のドライバが準備されていると仮定します。すなわち、回転数指令値を操作入力として考えます。そうすると、上のモデルを書き替える必要が出てきます。

一般にモノを動かすときは、増速、定速、減速の速度パターンを与えて、面積が望ましい変位となるようにします。いま、モータの回転数制御系のダイナミックスを、十分小さな時定数T_aを用いて

\displaystyle{(3)\quad\frac{d}{dt}\dot{\theta}=-\frac{1}{T_a}\dot{\theta}+\frac{1}{T_a}\dot{\theta}_c}

で表しておきます。次に、上の運動方程式の第２ブロック行は

\displaystyle{(4)\quad \underbrace{ \left[\begin{array}{c} m_1 \\ \vdots \\ m_N \\ \end{array}\right] }_{M_{21}} \ddot{\theta} + \underbrace{ \left[\begin{array}{c} \ddot{r}_1 \\ \vdots \\ \ddot{r}_N \end{array}\right] }_{\ddot{r}} + \underbrace{ \left[\begin{array}{ccc} \Omega_1^2 & & 0 \\ & \ddots & \\ 0 & & \Omega_N^2 \\ \end{array}\right] }_{K_{22}} \underbrace{ \left[\begin{array}{c} r_1 \\ \vdots \\ r_N \end{array}\right] }_{r} = \left[\begin{array}{c} 0 \\ \vdots \\ 0 \end{array}\right] }

となっていますが、これを

\displaystyle{(5)\quad z=M_{21}\theta+r}

を用いて、次のように書き替えます。

\displaystyle{(6)\quad\underbrace{\frac{d^2}{dt^2}(M_{21}\theta+r)}_{\ddot{z}}+K_{22}(z-M_{21}\theta)=0}

以上から、制御系設計用のモデルとして次式を得ます。

\displaystyle{(7a)\quad \underbrace{ \left[\begin{array}{c} \dot{z} \\ \dot{\theta} \\\hline \ddot{z} \\ \ddot{\theta} \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cc|cc} 0_{N\times N} & 0_{N\times 1} & I_N & 0_{N\times 1} \\ 0_{1\times N} & 0 & 0_{1\times N} & 1 \\\hline -K_{22} & K_{22}M_{21} & 0_{N\times N} & 0_{N\times 1} \\ 0_{1\times N} & 0 & 0_{1\times N} & -\frac{1}{T_a} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {z} \\ {\theta} \\\hline \dot{z} \\ \dot{\theta} \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{cc} 0_{N\times 1} \\ 0 \\\hline 0_{N\times 1} \\ \frac{1}{T_a} \end{array}\right] }_{B} \dot{\theta}_c}
\displaystyle{(7b)\quad \eta(\xi)= \underbrace{ \left[\begin{array}{ccc} \phi_1(\xi) & \cdots & \phi_N(\xi) \end{array}\right] \left[\begin{array}{cc|cc} I_N & -M_{21} & 0_{N\times N} & 0_{N\times 1} \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {z} \\ {\theta} \\\hline \dot{z} \\ \dot{\theta} \end{array}\right] }_{x} }

演習B64b…Flipped Classroom

\displaystyle{L=1} m
\displaystyle{D=0.005} m
\displaystyle{\rho=7980} kg/m^3
\displaystyle{E=1950000} N/mm^2
\displaystyle{I=\frac{\pi}{64}D^4 } m^4
\displaystyle{A=\frac{\pi}{4}D^2} m^2
\displaystyle{J_0(=\frac{1}{2}M_mr_m^2) } kgm^2
\displaystyle{M_e=0.1} kg

MATLAB 
%beam1.m
%-----
 clear all, close all
 L=1; D=0.005; A=pi/4*D^2;
 rho=7980; E=195000e6; I=pi/64*D^4;
 J0=0.5*1*0.01; Me=0.1;
 J=(J0+1/3*rho*A*L^3+Me*L^2)/(rho*A*L^3); beta=Me/(rho*A*L);
 T=1/(L^2*sqrt(rho*A/(E*I)));
%-----
 x=0:0.01:30;
 y=err(x,beta);
 figure(1)
 plot(x,sign(y),[0 30],[0 0]); axis([0 30 -2 2]);
 N=3;
% w=locate(N);
% om=[];
% for i=1:N
%   x=fsolve(w(1,i),err); om=[om x];
% end
 om=[1.3604444, 4.0789285, 7.1670195]
%-----
 n_d=100; h=L/n_d; x=0:h:L; x=x';
 y1=phi(1,x,om); y2=phi(2,x,om); y3=phi(3,x,om);
 figure(2)
 plot(x,y1,x,y2,x,y3)
%-----
 y1=dphi(1,x,om); y2=dphi(2,x,om); y3=dphi(3,x,om);
 figure(2)
 plot(x,y1,x,y2,x,y3)
%-----
 mode_fun=[]; 
 for i=1:N, mode_fun =[mode_fun  phi(i,x,om) ]; end
 mode_dfun=[]; 
 for i=1:N, mode_dfun=[mode_dfun dphi(i,x,om)]; end
%-----
 t_simp=zeros(1, n_d+1);
 t_simp(1,1)=1/3*h;
 for i=1:(n_d+1)/2, t_simp(1,2*i)=2/3*h; t_simp(1,2*i+1)=4/3*h; end
 t_simp(1,n_d+1)=1/3*h;
 c_simp=t_simp(1,:);
%-----
 for i=1:N
   m(i) = c_simp*(x.*mode_fun(:,i))+beta*mode_fun(n_d+1,i);
   OM2(i)=om(i)^4;
 end
% M=[J m';m eye(N,N)];
% K=zeros(N+1,N+1); K(2:N+1,2:N+1)=diag(OM2);
% AA=[zeros(N+1,N+1) eye(N+1,N+1);
% -M\K zeros(N+1,N+1)];
% B =[zeros(N+1,1); M\eye(N+1,1)];
 M21=m'; K22=diag(OM2);
 AA=[zeros(N+1,N+1) eye(N+1,N+1);
     [-K22 K22*m';zeros(1,N+1)] zeros(N+1,N+1)];
 Ta=0.1;
 AA(2*(N+1),2*(N+1))=-1/Ta; 
 B =[zeros(2*N+1,1); 1/Ta];
 C=[eye(N,N) -M21 zeros(N,N+1)];
 CC=[mode_fun*C;zeros(1,N),1,zeros(1,N+1)];
 P=ss(AA,B,CC,[]);
 G=tf(P);
 w=logspace(-3,3,100);
 figure(3)
 for i=[1,21,41,61,81,101]
   bode(G(i,1),w),hold on
 end
%-----
 sys=ss(AA,B,CC([1,21,41,61,81,101],:),[])
 state0=(L/(E*I))*B; 
%state0=zeros(8,1)
 t0=0; t=0:T/1000:T;
 figure(4)
 initial(sys,state0,t)
%-----
 [F,p]=opt(AA,B,eye(8,8),eye(8,8),1)
 ACL=AA-B*F;
 sysCL=ss(ACL,B,CC([1,21,41,61,81,101],:),[]) 
 figure(5)
 initial(sysCL,state0,t)
%====
 function errval=err(x,beta)
   errval=beta*x.*(cosh(x).*sin(x)-sinh(x).*cos(x))-(1+cosh(x).*cos(x));
 end
%-----
 function y=phi(i,x,om)
   y=(sinh(om(i))+sin(om(i))).*(cosh(om(i)*x)-cos(om(i)*x))...
   -(cosh(om(i))+cos(om(i))).*(sinh(om(i)*x)-sin(om(i)*x));
   y=y/max(abs(y));
 end
%-----
 function y=dphi(i,x,om)
   y=(sinh(om(i))+sin(om(i))).*(sinh(om(i)*x)+sin(om(i)*x))...
   -(cosh(om(i))+cos(om(i))).*(cosh(om(i)*x)-cos(om(i)*x));
   y=om(i)*y;
   y=y/max(abs(y));
 end
%-----
%eof
SCILAB

柔軟ビーム

投稿日時: 投稿者:

[1] 柔軟構造物の制振問題をどう取り扱うかを考えて行きます。

この写真は米国の某大学で製作された制御実験装置です。柔軟ビームの一端をロボットアームで把持し、振動を抑制しながら回転させることが制御目的です。この柔軟ビームは水平方向ばかりでなく、垂直方向にも振動するのですが、ここでは水平方向だけの振動抑制問題を考えます。参考にしたのは次の文献です。

阿部・児島著：「無駄時間・分布定数系の制御」、コロナ社、2007

この本の６章「振動系」では、次のような制御対象「柔軟ビーム」を扱っています。

ここで、ビームの長さをL[m]、断面積をA[m^2]、密度をA[kg/m^3]、縦弾性係数をE[N/m^2]、断面２次モーメントをI[m^4]、ハブの回転慣性モーメントをJ_0[kgm^2]、ペイロードをM_e[kg]とします。また、時刻t[sec]における、ハブの回転角を\theta(t)[rad]、ハブの支持点から距離x[m]のビーム上の点の弾性変位をy(x,t)[m]、その点の座標を(X(t),Y(t))とします。

このときハブの支持点x=0において次が成り立ちます。

\displaystyle{(1)\quad y(0,t)=\frac{\partial y(0,t)}{\partial x}=0}

いま、\sin(\cdot)、\cos(\cdot)をそれぞれS_{(\cdot)}、C_{(\cdot)}と略記すると、次式が成り立ちます。

\displaystyle{(2a)\quad X=xC_\theta-yS_\theta\ \Rightarrow\ \dot{X}=-xS_\theta\dot{\theta}-\dot{y}S_\theta-yC_\theta\dot{\theta}=-(xS_\theta+yC_\theta)\dot{\theta}-\dot{y}S_\theta }
\displaystyle{ \Rightarrow\ \dot{X}^2=(xS_\theta+yC_\theta)^2\dot{\theta}^2+2(xS_\theta+yC_\theta)\dot{\theta}\dot{y}S_\theta+\dot{y}^2S_\theta^2 }
\displaystyle{(2b)\quad Y=xS_\theta+yC_\theta\ \Rightarrow\ \dot{Y}=xC_\theta\dot{\theta}+\dot{y}C_\theta-yS_\theta\dot{\theta}=(xC_\theta-yS_\theta)\dot{\theta}+\dot{y}C_\theta}
\displaystyle{ \Rightarrow\ \dot{Y}^2=(xC_\theta-yS_\theta)^2\dot{\theta}^2+2(xC_\theta-yS_\theta)\dot{\theta}\dot{y}C_\theta+\dot{y}^2C_\theta^2 }
\displaystyle{(2c)\quad \dot{X}^2+\dot{Y}^2=(x^2+y^2)\dot{\theta}^2+2x\dot{\theta}\dot{y}+\dot{y}^2 }

以下では、双曲線関数\sinh(\cdot)、\cosh(\cdot)をそれぞれS^h_{(\cdot)}、C^h_{(\cdot)}と略記し、次の基本式を多用します。

\displaystyle{(3a)\quad C_x^h=\frac{1}{2}(e^x+e^{-x}),S_x^h=\frac{1}{2}(e^x-e^{-x}) },
\displaystyle{(3b)\quad C_x=\frac{1}{2}(e^{jx}+e^{-jx}),S_x=\frac{1}{2}(e^{jx}-e^{-jx}) }
\displaystyle{(3c)\quad (C_x^h)'=S_x^h,\ (S_x^h)'=C_x^h}, \displaystyle{C_x'=-S_x,\ S_x'=C_x}
\displaystyle{(3d)\quad C_x^{h2}-S_x^{h2}=1}, \displaystyle{C_x^2+S_x^2=1}

また、次のような変分を取ることを頻繁に行います。

\displaystyle{(4)\quad \begin{array}{l} y=f(x) \Rightarrow \\ \delta y=\epsilon(\left.\frac{f(x+\epsilon x_\epsilon)}{d\epsilon}\right|_{\epsilon=0})=\epsilon(\left.f'(x+\epsilon x_\epsilon)x_\epsilon\right|_{\epsilon=0})=f'(x)\epsilon x_\epsilon=f'(x)\delta x \end{array} }

●柔軟ビームの運動方程式を、Hamilton Principleに基づいて導出します。

\displaystyle{(5a)\quad \int_{t_0}^{t_1} ( \delta L + \delta W ) dt =\int_{t_0}^{t_1} ( \delta T-\delta V + \delta W ) dt=0 }

\displaystyle{(5b)\quad \delta T=\delta T_0+ \delta T_1+ \delta T_2 }

●\displaystyle{(6a)\quad T_0 = \frac{1}{2} J_0\dot{\theta}^2 }

\displaystyle{(6b)\quad T_0(\epsilon) = \frac{1}{2} J_0(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)^2 }

\displaystyle{(6c)\quad \frac{dT_0(\epsilon) }{d\epsilon}=J_0(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)\dot{\theta}_\epsilon }

\displaystyle{(6d)\quad \delta T_0=\epsilon\left(\left.\frac{dT_0(\epsilon)}{d\epsilon}\right|_{\epsilon=0}\right)=J_0\dot{\theta}\epsilon\dot{\theta}_\epsilon=J_0\dot{\theta}\delta\dot{\theta} }

\displaystyle{(6e)\quad \int_{t_0}^{t_1}\delta T_0 dt= \int_{t_0}^{t_1} J_0\dot{\theta}\delta\dot{\theta} dt = \left[J_0\dot{\theta}\delta{\theta} \right]_{t_0}^{t_1} -\int_{t_0}^{t_1} J_0\ddot{\theta}\delta{\theta} dt = -\int_{t_0}^{t_1} J_0\ddot{\theta}\delta{\theta} dt }

●\displaystyle{(7a)\quad T_1 = \frac{1}{2} \rho A \int_0^{L} (\dot{X}^2+\dot{Y}^2)dx }
\displaystyle{= \frac{1}{2} \rho A \int_0^{L} ((x^2+y^2)\dot{\theta}^2+2x\dot{\theta}\dot{y}+\dot{y}^2)dx }

\displaystyle{(7b)\quad T_1(\epsilon)}
\displaystyle{ = \frac{1}{2} \rho A \int_0^{L} ((x^2+({y}+\epsilon{y}_\epsilon)^2)(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)^2+2x(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)(\dot{y}+\epsilon\dot{y}_\epsilon)+(\dot{y}+\epsilon\dot{y}_\epsilon)^2 )dx }

\displaystyle{(7c)\quad \frac{dT_1(\epsilon) }{d\epsilon}=\rho A \int_0^{L} ( ({y}+\epsilon{y}_\epsilon){y}_\epsilon(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)^2+(x^2+({y}+\epsilon{y}_\epsilon)^2)(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)\dot{\theta}_\epsilon }
\displaystyle{ +x\dot{\theta}_\epsilon(\dot{y}+\epsilon\dot{y}_\epsilon)+x(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)\dot{y}_\epsilon +(\dot{y}+\epsilon\dot{y}_\epsilon)\dot{y}_\epsilon )dx }

\displaystyle{(7d)\quad \delta T_1=\epsilon\left(\left.\frac{dT_1(\epsilon)}{d\epsilon}\right|_{\epsilon=0}\right)}
\displaystyle{=\rho A \int_0^{L} ( {y}\epsilon{y}_\epsilon\dot{\theta}^2+(x^2+{y}^2)\dot{\theta\epsilon}\dot{\theta}_\epsilon +x\epsilon\dot{\theta}_\epsilon\dot{y}+x\dot{\theta}\epsilon\dot{y}_\epsilon+\dot{y}\epsilon\dot{y}_\epsilon )dx }
\displaystyle{=\rho A \int_0^{L} ({y}\dot{\theta}^2\epsilon{y}_\epsilon+(x^2+{y}^2)\dot{\theta\epsilon}\dot{\theta}_\epsilon+x\dot{y}\epsilon\dot{\theta}_\epsilon+x\dot{\theta}\epsilon\dot{y}_\epsilon+\dot{y}\epsilon\dot{y}_\epsilon )dx }
\displaystyle{=\rho A \int_0^{L} ({y}\dot{\theta}^2\delta{y}+((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta\dot{\theta}+(x\dot{\theta}+\dot{y})\delta\dot{y} )dx }

\displaystyle{(7e)\quad \int_{t_0}^{t_1}\delta T_1 dt=\int_{t_0}^{t_1}\rho A \int_0^{L} ({y}\dot{\theta}^2\delta{y}+((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta\dot{\theta}+(x\dot{\theta}+\dot{y})\delta\dot{y} )dxdt }
\displaystyle{ =\int_{t_0}^{t_1}\rho A \int_0^{L} {y}\dot{\theta}^2\delta{y}dxdt+\int_{t_0}^{t_1}\rho A \int_0^{L} ((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta\dot{\theta}dxdt+\int_{t_0}^{t_1}\rho A \int_0^{L} (x\dot{\theta}+\dot{y})\delta\dot{y} dxdt }
\displaystyle{ =\int_{t_0}^{t_1}\rho A \int_0^{L} {y}\dot{\theta}^2\delta{y}dxdt }
\displaystyle{ +\left[\rho A \int_0^{L} (x^2+{y}^2)\dot{\theta}+x\dot{y})\delta{\theta}dx\right]_{t_0}^{t_1}-\int_{t_0}^{t_1}\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta{\theta}dxdt }
\displaystyle{ +\left[\rho A \int_0^{L} (x\dot{\theta}+\dot{y}) \delta{y}dx\right]_{t_0}^{t_1}-\int_{t_0}^{t_1}\rho A \int_0^{L}\frac{d}{dt}(x\dot{\theta}+\dot{y})\delta{y}dxdt }
\displaystyle{ =-\int_{t_0}^{t_1}\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta{\theta}dxdt -\int_{t_0}^{t_1}\rho A \int_0^{L} ({y}\dot{\theta}^2+\frac{d}{dt}(x\dot{\theta}+\dot{y}))\delta{y}dxdt }

●\displaystyle{(8a)\quad T_2 = \frac{1}{2} M_e (\dot{X}_L^2+\dot{Y}_L^2) }
\displaystyle{= \frac{1}{2} M_e ((L^2+y^2(L))\dot{\theta}^2+2L\dot{\theta}\dot{y}(L)+\dot{y}^2(L)) }

\displaystyle{(8b)\quad T_2(\epsilon) = \frac{1}{2} M_e ((L^2+({y}(L)+\epsilon{y}_\epsilon(L))^2)(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)^2 }
\displaystyle{+2L(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)(\dot{y}(L)+\epsilon\dot{y}_\epsilon(L))+(\dot{y}(L)+\epsilon\dot{y}_\epsilon(L))^2 ) }

\displaystyle{(8c)\quad \frac{dT_2(\epsilon)}{d\epsilon}}
\displaystyle{=M_e (({y}(L)+\epsilon{y}_\epsilon(L)){y}_\epsilon(L)(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)^2+(L^2+({y}(L)+\epsilon{y}_\epsilon(L))^2)(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)\dot{\theta}_\epsilon }
\displaystyle{+L\dot{\theta}_\epsilon(\dot{y}(L)+\epsilon\dot{y}_\epsilon(L))+L(\dot{\theta}+\epsilon\dot{\theta}_\epsilon)\dot{y}_\epsilon(L)+(\dot{y(L)}+\epsilon\dot{y}_\epsilon(L))\dot{y}_\epsilon(L) ) }

\displaystyle{(8d)\quad \delta T_2=M_e ({y}(L)\dot{\theta}^2\delta{y}(L)+((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L))\delta\dot{\theta}+(L\dot{\theta}+\dot{y}(L))\delta\dot{y}(L) ) }

\displaystyle{(8e)\quad \int_{t_0}^{t_1}\delta T_2 dt }
\displaystyle{=-\int_{t_0}^{t_1}M_e \frac{d}{dt}((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L))\delta{\theta}dt -\int_{t_0}^{t_1}M_e ({y}(L)\dot{\theta}^2+\frac{d}{dt}(L\dot{\theta}+\dot{y}(L)))\delta{y}(L)dt }

●\displaystyle{(9a)\quad V=\frac{1}{2} EI \int_0^{L} y''^2 dx\quad(y''=\frac{\partial^2 y}{\partial x^2}) }

\displaystyle{(9b)\quad V(\epsilon) = \frac{1}{2} EI \int_0^{L} (y'' +\epsilon{y_\epsilon''} )^2 dx }

\displaystyle{(9c)\quad frac{dV(\epsilon) }{d\epsilon}=EI \int_0^{L} (y'' +\epsilon{y_\epsilon''} ){y_\epsilon''}dx }

\displaystyle{(9d)\quad \delta V=\epsilon\left(\left.\frac{dV(\epsilon)}{d\epsilon}\right|_{\epsilon=0}\right) = EI \int_0^{L} y'' \epsilon{y_\epsilon''} dx = EI \int_0^{L} y'' \delta y'' dx }
\displaystyle{ = EI \left[y'' \delta y'\right]_0^{L} - EI \int_0^{L} y''' \delta y' dx }
\displaystyle{ = EI y''(L) \delta y'(L)-EI \left[y''' \delta y\right]_0^{L} + EI \int_0^{L} y'''' \delta y dx }
\displaystyle{ = EI y''(L) \delta y'(L)-EI y'''(L) \delta y(L) + EI \int_0^{L} y'''' \delta y dx }

\displaystyle{(9e)\quad \int_{t_0}^{t_1}\delta V dt=\int_{t_0}^{t_1} (EI y''(L) \delta y'(L)-EI y'''(L) \delta y(L) + EI \int_0^{L} y'''' \delta y dx )dt }

●\displaystyle{(10a)\quad W = u\theta }

\displaystyle{(10b)\quad W(\epsilon) = u(\theta+\epsilon\theta_\epsilon) }

\displaystyle{(10c)\quad \frac{dW(\epsilon) }{d\epsilon}= u\theta_\epsilon }

\displaystyle{(10d)\quad \delta W=\epsilon\left(\left.\frac{dW(\epsilon)}{d\epsilon}\right|_{\epsilon=0}\right)=u\epsilon\theta_\epsilon=u\delta\theta }

\displaystyle{(10e)\quad \int_{t_0}^{t_1}\delta W dt= \int_{t_0}^{t_1} u\delta\theta dt }

●以上から、ハミルトン原理を表す上式の各項を

\displaystyle{(11a)\quad \int_{t_0}^{t_1}\delta T_0 dt=-\int_{t_0}^{t_1} J_0\ddot{\theta}\delta{\theta} dt }
\displaystyle{(11b)\quad \int_{t_0}^{t_1}\delta T_1 dt=　}
\displaystyle{-\int_{t_0}^{t_1}\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta{\theta}dxdt -\int_{t_0}^{t_1}\rho A \int_0^{L} ({y}\dot{\theta}^2+\frac{d}{dt}(x\dot{\theta}+\dot{y}))\delta{y}dxdt }
\displaystyle{(11c)\quad \int_{t_0}^{t_1}\delta T_2 dt= }
\displaystyle{=-\int_{t_0}^{t_1}M_e \frac{d}{dt}((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L))\delta{\theta}dt -\int_{t_0}^{t_1}M_e ({y}(L)\dot{\theta}^2+\frac{d}{dt}(L\dot{\theta}+\dot{y}(L)))\delta{y}(L)dt }
\displaystyle{(11d)\quad \int_{t_0}^{t_1}\delta V dt=\int_{t_0}^{t_1} (EI y''(L) \delta y'(L)-EI y'''(L) \delta y(L) + EI \int_0^{L} y'''' \delta y dx )dt }
\displaystyle{(11e)\quad \int_{t_0}^{t_1}\delta W dt= \int_{t_0}^{t_1} u\delta\theta dt }

のように得たので、次式が成り立ちます。

\displaystyle{(12)\quad \int_{t_0}^{t_1} ( \delta T_0+ \delta T_1+ \delta T_2-\delta V + \delta W ) dt= }
\displaystyle{ -\int_{t_0}^{t_1} J_0\ddot{\theta}\delta{\theta} dt }
\displaystyle{ -\int_{t_0}^{t_1}\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})\delta{\theta}dxdt -\int_{t_0}^{t_1}\rho A \int_0^{L} ({y}\dot{\theta}^2+\frac{d}{dt}(x\dot{\theta}+\dot{y}))\delta{y}dxdt }
\displaystyle{ -\int_{t_0}^{t_1}M_e \frac{d}{dt}((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L))\delta{\theta}dt -\int_{t_0}^{t_1}M_e ({y}(L)\dot{\theta}^2+\frac{d}{dt}(L\dot{\theta}+\dot{y}(L)))\delta{y}(L)dt }
\displaystyle{ -\int_{t_0}^{t_1} (EI y''(L) \delta y'(L)-EI y'''(L) \delta y(L) + EI \int_0^{L} y'''' \delta y dx )dt }
\displaystyle{ +\int_{t_0}^{t_1} u\delta\theta dt }
\displaystyle{=-\int_{t_0}^{t_1}\underline{ (J_0\ddot{\theta}+M_e \frac{d}{dt}((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L))+\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})dx -u)} \delta\theta dt }
\displaystyle{ -\int_{t_0}^{t_1}\underline{ (\rho A \int_0^{L} ({y}\dot{\theta}^2+\frac{d}{dt}(x\dot{\theta}+\dot{y})) dx + EI \int_0^{L} y'''' dx)}\delta y dt }
\displaystyle{ -\int_{t_0}^{t_1}\underline{(M_e ({y}(L)\dot{\theta}^2+\frac{d}{dt}(L\dot{\theta}+\dot{y}(L)))-EI y'''(L) )}\delta{y}(L)dt }
\displaystyle{ -\int_{t_0}^{t_1} \underline{EI y''(L) }\delta y'(L) dt }=0

これから各変分の係数（下線部）を0と置いて、柔軟ビームの運動方程式として次式を得ます。

\displaystyle{(13a)\quad J_0\ddot{\theta} +M_e \frac{d}{dt}((L^2+{y}^2(L))\dot{\theta}+L\dot{y}(L)) +\rho A \int_0^{L} \frac{d}{dt}((x^2+{y}^2)\dot{\theta}+x\dot{y})dx -u=0 }
\displaystyle{(13b)\quad \rho A \int_0^{L} ({y}\dot{\theta}^2+\frac{d}{dt}(x\dot{\theta}+\dot{y})) dx + EI \int_0^{L} y'''' dx =0 }
\displaystyle{(13c)\quad M_e ({y}(L)\dot{\theta}^2+\frac{d}{dt}(L\dot{\theta}+\dot{y}(L)))-EI y'''(L) =0 }
\displaystyle{(13d)\quad EI y''(L) }=0

●高次項や高階微分項は微小であるとして、これらの近似を行います。まず第１式(13a)は

\displaystyle{(14a)\quad J_0\ddot{\theta} +M_e L^2\ddot{\theta}+M_e\frac{d}{dt}({y}^2(L)\dot{\theta})+M_eL\ddot{y}(L) }
\displaystyle{+\rho A \int_0^{L} (x^2\ddot{\theta}+\frac{d}{dt}({y}^2\dot{\theta}))dx +\rho A \int_0^{L}x\ddot{y}dx -u=0 }

においてy^2が微小であるとして、次のように近似します。

\displaystyle{(14a')\quad \underbrace{(J_0+\frac{1}{3}\rho AL^3+M_eL^2)}_J \ddot{\theta}+M_e L\ddot{y}(L)+\rho A \int_0^{L} x\ddot{y}dx -u=0 }

また、第２式(13b)と第３式(13c)は、\dot{\theta}^2が微小であるとして、次のように近似します。

\displaystyle{(14b)\quad \rho A \int_0^{L} (x\ddot{\theta}+\ddot{y}) dx + EI \int_0^{L} y'''' dx =0 }
\displaystyle{(14c)\quad M_e (L\ddot{\theta}+\ddot{y}(L))-EI y'''(L) =0 }

これらを次のように無次元化します。

(15a)\quad \frac{L}{EI}\times第１式(14a’)：
\displaystyle{\underbrace{\frac{J}{\rho AL^3}}_{\alpha} \underbrace{\frac{d^2\theta}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\theta}} +\underbrace{\frac{M_e}{\rho AL}}_{\beta} \underbrace{\frac{L}{L}}_1 \underbrace{\frac{d^2\frac{y(L)}{L}}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\eta}(1)} +\int_0^{\frac{L}{L}} \underbrace{\frac{x}{L}}_\xi\underbrace{\frac{d^2\frac{y}{L}}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\eta}} \underbrace{d\frac{x}{L}}_{d\xi} -\underbrace{\frac{u}{\frac{EI}{L}}}_{u'}=0 }
(15b)\quad \frac{L^2}{EI}\times第２式(14b)：
\displaystyle{\int_0^{\frac{L}{L}} (\underbrace{\frac{x}{L}}_\xi\underbrace{\frac{d^2\theta}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\theta}} +\underbrace{\frac{d^2\frac{y}{L}}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\eta}}) \underbrace{d\frac{x}{L}}_{d\xi} + \int_0^{\frac{L}{L}} \underbrace{\frac{d^4\frac{y}{L}}{d(\frac{x}{L})^4}}_{\eta''''} \underbrace{d\frac{x}{L}}_{d\xi} =0 }
(15c)\quad \frac{L^2}{EI}\times第３式(14c)：
\displaystyle{\underbrace{\frac{M_e}{\rho AL}}_{\beta} (\underbrace{\frac{L}{L}}_1\underbrace{\frac{d^2\theta}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\theta}}+\underbrace{\frac{d^2\frac{y(L)}{L}}{d(\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}}_{\ddot{\eta}(1)})-\underbrace{\frac{d^3\frac{y(L)}{L}}{d(\frac{x}{L})^3}}_{\eta'''(1)} =0}
(15d)\quad \frac{L}{EI}\times第４式(13d)：
\displaystyle{\underbrace{\frac{d^2\frac{y}{L}}{d(\frac{x}{L})^2}}_{\eta''(1)} }=0}

すなわち、代表長さL、代表時間T=L^2\sqrt{\frac{\rho A}{EI}}、代表トルク\frac{EI}{L}として、無次元化した柔軟ビームの運動方程式は次式となります。

\displaystyle{(16a)\quad \alpha\ddot{\theta}+\beta\ddot{\eta}(1)+\int_0^1\xi\ddot{\eta}d\xi-u'=0}
\displaystyle{(16b)\quad \int_0^1(\xi\ddot{\theta}+\ddot{\eta}+\eta'''')d\xi=0}
\displaystyle{(16c)\quad \beta (\ddot{\theta}+\ddot{\eta}(1))-\eta'''(1)=0}
\displaystyle{(16d)\quad \eta''(1)=0}

ただし

\displaystyle{(17a)\quad T=L^2\sqrt{\frac{\rho A}{EI}},\ \tau=\frac{t}{T}}
\displaystyle{(17b)\quad \xi=\frac{x}{L},\ \eta=\frac{y}{L},\ u'=\frac{u}{\frac{EI}{L}}}
\displaystyle{(17c)\quad \alpha=\frac{J}{\rho AL^3}=\frac{J_0+\frac{1}{3}\rho AL^3+M_eL^2}{\rho AL^3}=\frac{J_0}{\rho AL^3}+\frac{1}{3}+\beta }
\displaystyle{(17d)\quad \beta=\frac{M_e}{\rho AL}}

[2] 以下では、モード法による求解方法を考えていきます。

そのために、\ddot{\theta}=0を仮定すると

\displaystyle{(18a)\quad \eta''''+\ddot{\eta}=0}
\displaystyle{(18b)\quad \eta'''(1)=\beta\ddot{\eta}(1)}
\displaystyle{(18c)\quad \eta(0)=\eta'(0)=\eta''(1)=0}

を得ます。ここで、変数分離\eta(\xi,\tau)=\phi(\xi)r(\tau)を行います。

\displaystyle{(19a)\quad \phi''''r +\phi\ddot{r} =0 }
\displaystyle{(19b)\quad \phi'''(1)r =\beta\phi(1)\ddot{r} }
\displaystyle{(19c)\quad \phi(0)r=\phi'(0)r=\phi''(1)r =0 }

さらに、\displaystyle{\Omega^2=-\frac{\ddot{r}}{r} },　\displaystyle{ \ddot{r}+\Omega^2r =0 }を仮定しますと

\displaystyle{(20a)\quad \phi'''' -\Omega^2\phi =0 }
\displaystyle{(20b)\quad \phi'''(1) =-\beta\Omega^2\phi(1) }
\displaystyle{(20c)\quad \phi(0)=\phi'(0)=\phi''(1) =0 }

を得ます。これから

\displaystyle{(21a)\quad \phi'''' -\Omega^2\phi =0 }
\displaystyle{\cdot\qquad\Downarrow\qquad \omega^4=\Omega^2}
\displaystyle{(21b)\quad \phi'''' -\omega^4\phi =0 }
\displaystyle{\cdot\qquad\Downarrow\qquad \phi(\xi)=c_1e^{\omega\xi}+c_2e^{-\omega\xi}+c_3e^{j\omega\xi}+c_4e^{-j\omega\xi}}
\displaystyle{(21c)\quad \phi(\xi)=c_1C_{\omega\xi}+c_2S_{\omega\xi}+c_3C^h_{\omega\xi}+c_4S^h_{\omega\xi}}
\displaystyle{\cdot\qquad\Downarrow\qquad\phi(0)=\phi'(0)=\phi''(1) =0}
\displaystyle{(21d)\quad c_1+c_3=0,c_2+c_4=0}
\displaystyle{c_3(C^h_{\omega}+C_{\omega})+c_4(S^h_{\omega}+S_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow\qquad for\ arbitrary\ \gamma}
\displaystyle{(21f)\quad c_3=-c_1=\alpha(S^h_{\omega}+S_{\omega}),c_4=-c_2=-\gamma(C^h_{\omega}+C_{\omega})}
\displaystyle{\cdot\qquad\Downarrow\qquad }
\displaystyle{(21g)\quad \phi(x)=\gamma((S^h_{\omega}+S_{\omega})(C^h_{\omega\xi}-C_{\omega\xi})-(C^h_{\omega}+C_{\omega})(S^h_{\omega\xi}-S_{\omega\xi}))}

を得ます。さらに\omegaが満足すべき制約式を次のように得ます。

\displaystyle{(22a)\quad \phi=\gamma((S^h_{\omega}+S_{\omega})(C^h_{\omega\xi}-C_{\omega\xi})-(C^h_{\omega}+C_{\omega})(S^h_{\omega\xi}-S_{\omega\xi}))}
\displaystyle{\phi'=\gamma\omega((S^h_{\omega}+S_{\omega})(S^h_{\omega\xi}+S_{\omega\xi})-(C^h_{\omega}+C_{\omega})(C^h_{\omega\xi}-C_{\omega\xi}))}
\displaystyle{\phi''=\gamma\omega^2((S^h_{\omega}+S_{\omega})(C^h_{\omega\xi}+C_{\omega\xi})-(C^h_{\omega}+C_{\omega})(S^h_{\omega\xi}+S_{\omega\xi}))}
\displaystyle{\phi'''=\gamma\omega^3((S^h_{\omega}+S_{\omega})(S^h_{\omega\xi}-S_{\omega\xi})-(C^h_{\omega}+C_{\omega})(C^h_{\omega\xi}+C_{\omega\xi}))}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(22b)\quad \phi(1)=2\gamma(C^h_{\omega}S_{\omega}-S^h_{\omega}C_{\omega})}
\displaystyle{\phi'''(1)=-2\gamma\omega^3(1+C^h_{\omega}C_{\omega})}
\displaystyle{\cdot\qquad\Downarrow\qquad \phi'''(1) =-\beta\Omega^2\phi(1)}
\displaystyle{(22c)\quad -2\gamma\omega^3(1+C^h_{\omega}C_{\omega}) =-\beta\omega^42\gamma(C^h_{\omega}S_{\omega}-S^h_{\omega}C_{\omega})}
\displaystyle{\cdot\qquad\Downarrow\qquad }
\displaystyle{(22d)\quad \beta\omega(C^h_{\omega}S_{\omega}-S^h_{\omega}C_{\omega})-(1+C^h_{\omega}C_{\omega})=0}

●これを満足する\omegaを、\omega_1,\omega_2,\cdotsとしますと、次のモード関数群を得たことになります。

\displaystyle{(23a)\quad \phi_i(x)=\gamma_i((S^h_{\omega_i}+S_{\omega_i})(C^h_{\omega_i\xi}-C_{\omega_i\xi})-(C^h_{\omega_i}+C_{\omega_i})(S^h_{\omega_i\xi}-S_{\omega_i\xi}))}

ただし

\displaystyle{(23b)\quad \beta\omega_i(C^h_{\omega_i}S_{\omega_i}-S^h_{\omega_i}C_{\omega_i})-(1+C^h_{\omega_i}C_{\omega_i})=0}

これらのモード関数はお互いに直交することが次のようにして示されます（i,j=1,2,\cdots）。

\displaystyle{(24a)\quad \int_0^1 \phi''_i\phi''_jd\xi=[ \phi''_i\phi'_j]_0^1-\int_0^1 \phi'''_i\phi'_jd\xi}
\displaystyle{=\underbrace{\phi''_i(1)}_0\phi'_j(1)-\phi''_i(0)\underbrace{\phi'_j(0)}_0-\int_0^1 \phi'''_i\phi'_jd\xi}
\displaystyle{=-[ \phi'''_i\phi_j]_0^1+\int_0^1 \phi''''_i\phi_jd\xi}
\displaystyle{=\underbrace{-\phi'''_i(1)}_{\beta\Omega_i^2\phi_i(1)}\phi_j(1)+\phi'''_i(0)\underbrace{\phi_j(0)}_0+\int_0^1 \Omega_i^2\phi_i\phi_jd\xi}
\displaystyle{=\Omega_i^2(\int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1))}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(24b)\quad \int_0^1 \phi''_i\phi''_jdx=\Omega_i^2(\int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1))}
\displaystyle{\int_0^1 \phi''_j\phi''_idx=\Omega_j^2(\int_0^1 \phi_j\phi_id\xi+\beta\phi_j(1)\phi_i(1))}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(24c)\quad 0=(\Omega_i^2-\Omega_j^2)(\int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1))}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(24d)\quad \int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1)=0\quad(i\ne j)}
\displaystyle{\int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1)\ne0\quad(i=j)}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(24e)\quad \int_0^1 \phi_i\phi_jd\xi+\beta\phi_i(1)\phi_j(1)=\delta_{ij}}
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(24f)\quad \int_0^1 \phi''_i\phi''_jd\xi=\Omega_i^2\delta_{ij}}

●結局、弾性変位\eta(\xi,\tau)を、モード関数\phi_i(\xi)を時間関数r_i(\tau)で重み付けて

\displaystyle{(25)\quad \eta(\xi,\tau)=\sum_{i=0}^\infty\phi_i(\xi)r_i(\tau) }

のように表します。この時間関数r_i(\tau)の支配方程式は

\displaystyle{(26a)\quad \alpha\ddot{\theta}+\int_0^1\xi\ddot{\eta}d\xi+\beta\ddot{\eta}(1)=u'}
\displaystyle{(26b)\quad \int_0^1(\xi\ddot{\theta}+\ddot{\eta}+\eta'''')d\xi=0}

から得られます。まず、第１式(26a)から

\displaystyle{(27)\quad \alpha\ddot{\theta}+\int_0^1\xi\ddot{\eta}d\xi+\beta\ddot{\eta}(1) }
\displaystyle{=\alpha\ddot{\theta}+\int_0^1 \xi\sum_{i=0}^\infty\phi_i\ddot{r}_i d\xi+ \beta\sum_{i=0}^\infty\phi_i(1)\ddot{r}_i }
\displaystyle{=\alpha\ddot{\theta}+\sum_{i=0}^\infty\underbrace{(\int_0^1 \xi\phi_i d\xi+ \beta\phi_i(1))}_{m_i}\ddot{r}_i=u' }

を得ます。次に

\displaystyle{(28)\quad \int_0^{1} \phi_j \eta'''' d\xi=[ \phi_j\eta''']_0^1-\int_0^1 \phi'_j\eta'''d\xi}
\displaystyle{=\phi_j(1)\eta'''(1)-\underbrace{\phi_j(0)}_0\eta'''(0)-\int_0^1 \phi'_j\eta'''d\xi}
\displaystyle{=\phi_j(1)\eta'''(1)-[\phi'_j\eta'']_0^1+\int_0^1 \phi''_j\eta''d\xi}
\displaystyle{=\phi_j(1)\eta'''(1)-\phi'_j(1)\underbrace{\eta''(1)}_0+\underbrace{\phi'_j(0)}_0\eta''(0)+\int_0^1 \phi''_j\eta''d\xi}
\displaystyle{=\phi_j(1)\beta(\ddot{\theta}+\ddot{y}(1))+\int_0^1 \phi''_j\sum_{i=0}^\infty\phi''_ir_id\xi}
\displaystyle{=\phi_j(1)\beta(\ddot{\theta}+\sum_{i=0}^\infty\phi_i(1)\ddot{r}_i)+\Omega^2_jr_j}

を考慮して、第２式(26b)から

\displaystyle{(29)\quad \int_0^1\phi_j(\xi\ddot{\theta}+\ddot{\eta}+\eta'''')d\xi }
\displaystyle{=\int_0^{1}\phi_j(\xi\ddot{\theta}+\sum_{i=0}^\infty\phi_i\ddot{r}_i)d\xi +\phi_j(1)\beta(\ddot{\theta}+\sum_{i=0}^\infty\phi_i(1)\ddot{r}_i)+\Omega^2_jr_j }
\displaystyle{=\underbrace{(\int_0^1 \xi\phi_jd\xi+\beta\phi_j(1))}_{m_j}\ddot{\theta} +\sum_{i=0}^\infty\underbrace{(\int_0^1 \phi_j\phi_id\xi+\beta\phi_j(1)\phi_i(1))}_{\delta_{ij}}\ddot{r}_i+\Omega^2_jr_j}
\displaystyle{=m_j\ddot{\theta}+\ddot{r}_j+\Omega_j^2r_j=0}

を得ます。すなわち時間関数r_i(\tau)は次式を満足すべきことがわかります。

\displaystyle{(30a)\quad \alpha\ddot{\theta}+\sum_{i=0}^\infty m_i\ddot{r}_i=u' }
\displaystyle{(30b)\quad m_i\ddot{\theta}+\ddot{r}_i+\Omega_i^2r_i=0\quad(i=1,0,\cdots)}

●これから、柔軟ビームに関するもう一つの運動方程式として、次を得ます。

\displaystyle{(31a)\quad \underbrace{ \left[\begin{array}{c|ccc} \alpha & m_1 & \dots &m_N \\ \hline m_1 & 1 & & 0 \\ \vdots & & \ddots & \\ m_N & 0 & & 1 \\ \end{array}\right] }_{M=\left[\begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array}\right]} \underbrace{ \left[\begin{array}{c} \ddot{\theta} \\ \hline \ddot{r}_1 \\ \vdots \\ \ddot{r}_N \end{array}\right] }_{\ddot{\xi}} + \underbrace{ \left[\begin{array}{c|ccc} 0 & 0 & \dots &0 \\ \hline 0 & \Omega_1^2 & & 0 \\ \vdots & & \ddots & \\ 0 & 0 & & \Omega_N^2 \\ \end{array}\right] }_{K=\left[\begin{array}{cc} K_{11} & K_{12} \\ K_{21} & M_{22} \end{array}\right]} \underbrace{ \left[\begin{array}{c} \theta \\ \hline r_1 \\ \vdots \\ r_N \end{array}\right] }_{\xi} = \underbrace{ \left[\begin{array}{c} 1 \\ \hline 0 \\ \vdots \\ 0 \end{array}\right] }_{B_2} u' }
\displaystyle{(31b)\quad \eta(\xi)= \underbrace{ \left[\begin{array}{c|ccc} 0 & \phi_1(\xi) & \dots &\phi_N(\xi) \end{array}\right] }_{C_1} \underbrace{ \left[\begin{array}{c} \theta \\ \hline r_1 \\ \vdots \\ r_N \end{array}\right] }_{\xi} }

ただし

\displaystyle{(32a)\quad T=L^2\sqrt{\frac{\rho A}{EI}},\ \tau=\frac{t}{T}}
\displaystyle{(32b)\quad \xi=\frac{x}{L},\ \eta=\frac{y}{L},\ u'=\frac{u}{\frac{EI}{L}}}
\displaystyle{(32c)\quad \alpha=\frac{J}{\rho AL^3},\ J=J_0+\frac{1}{3}\rho AL^3+M_eL^2 }
\displaystyle{(32d)\quad \beta=\frac{M_e}{\rho AL}}
\displaystyle{(32e)\quad m_i=\int_0^1 \xi\phi_jd\xi+\beta\phi_j(1)}
\displaystyle{(32f)\quad \Omega_i^2=\omega_i^4}
\displaystyle{(32g)\quad \beta\omega_i(C^h_{\omega_i}S_{\omega_i}-S^h_{\omega_i}C_{\omega_i})-(1+C^h_{\omega_i}C_{\omega_i})=0}

これから次の状態空間表現を得ます。

\displaystyle{(33a)\quad \underbrace{ \left[\begin{array}{c} \dot{\xi} \\ \ddot{\xi} \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cc} 0_{N+1\times N+1} & I_{N+1} \\ M^{-1}K & 0_{N+1\times N+1} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \xi \\ \dot{\xi} \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{cc} 0_{N+1\times 1} \\ M^{-1}B_2 \end{array}\right] }_{B} u'}
\displaystyle{(33b)\quad \eta(\xi)= \underbrace{ \left[\begin{array}{cc} C_1 & 0_{1\times N} \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \xi \\ \dot{\xi} \end{array}\right] }_{x} }

Note B64　モード関数

●ビーム（梁）の境界条件を、自由支持・ピン支持・固定支持と変えた場合のモード関数と振動数方程式を示します。ここで、Ax=0 が解 x\ne0 を持つための条件 {\rm det}A=0 が使われていることに注意してください。


（藤田勝久：振動工学、森北出版より）

\displaystyle{(1a)\quad M=-EI\frac{\partial^2 y}{\partial x^2} }
\displaystyle{(1b)\quad S=\frac{\partial M}{\partial x}=-EI\frac{\partial^3 y}{\partial x^3} }
\displaystyle{(1c)\quad \frac{\partial S}{\partial x}=-EI\frac{\partial^4 y}{\partial x^4} }

\displaystyle{(2a)\quad \rho A\frac{\partial^2 y}{\partial t^2}=\frac{\partial S}{\partial x} }, \displaystyle{m\frac{\partial^2 y(L)}{\partial t^2}=-S(L) }
\displaystyle{\cdot\qquad\Downarrow }
\displaystyle{(2b)\quad {\rho A}\frac{\partial^2 y}{\partial t^2}+{EI}\frac{\partial^4 y}{\partial x^4}=0 }, \displaystyle{m\frac{\partial^2 y(L)}{\partial t^2}-EI\frac{\partial^3 y(L)}{\partial x^3}=0 }
\displaystyle{\cdot\qquad\Downarrow \times \frac{L^3}{EI}, \times \frac{L^2}{EI} }
\displaystyle{(2c)\quad \frac{\partial^2 \frac{y}{L}}{\partial (\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}+\frac{\partial^4 \frac{y}{L}}{\partial (\frac{x}{L})^4}=0 }, \displaystyle{\frac{m}{\rho AL}\frac{\partial^2 \frac{y(L)}{L}}{\partial (\frac{t}{L^2\sqrt{\frac{\rho A}{EI}}})^2}-\frac{\partial^3 \frac{y(L)}{L}}{\partial (\frac{x}{L})^3}=0 }
\displaystyle{\cdot\qquad\Downarrow T=L^2\sqrt{\frac{\rho A}{EI}},\ \tau=\frac{t}{T},\ \xi=\frac{x}{L},\ \eta=\frac{y}{L},\ \beta=\frac{m}{\rho AL} }
\displaystyle{(2d)\quad \frac{\partial^2 \eta}{\partial \tau^2}+\frac{\partial^4 \eta}{\partial \xi^4}=0 }, \displaystyle{\beta\frac{\partial^2 \eta(1)}{\partial \tau^2}-\frac{\partial^3 \eta(1)}{\partial \xi^3}=0}
\displaystyle{\cdot\qquad\Downarrow \eta(\xi,\tau)=\phi(\xi)r(\tau) }
\displaystyle{(2e)\quad \phi(\xi)\ddot{r}(\tau)+\frac{\partial^4 \phi(\xi)}{\partial \xi^4}r(\tau)=0 }, \displaystyle{\beta\phi(1)\ddot{r}(\tau)-\frac{\partial^3 \phi(1)}{\partial \xi^3}r(\tau)=0 }
\displaystyle{\cdot\qquad\Downarrow \ddot{r}(\tau)=-\Omega^2r(\tau) }
\displaystyle{(2f)\quad -\phi(\xi)\Omega^2r(\tau)+\frac{\partial^4 \phi(\xi)}{\partial \xi^4}r(\tau)=0 }, \displaystyle{-\phi(1)\Omega^2r(\tau)-\frac{\partial^3 \phi(1)}{\partial \xi^3}r(\tau)=0 }
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(2g)\quad \frac{\partial^4 \phi}{\partial \xi^4}-\Omega^2\phi(\xi)=0 }, \displaystyle{\frac{\partial^3 \phi(1)}{\partial \xi^3}+\beta\Omega^2\phi(1)=0 }
\displaystyle{\cdot\qquad\Downarrow \omega^4=\Omega^2}
\displaystyle{(2h)\quad \phi''''(\xi)-\omega^4\phi(\xi)=0 }, \displaystyle{\phi'''(1)+\beta\omega^4\phi(1)=0 }
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(2i)\quad \phi(\xi)=c_1C_{\omega\xi}+c_2S_{\omega\xi}+c_3C^h_{\omega\xi}+c_4S^h_{\omega\xi}}
\displaystyle{\phi'(\xi)=\omega(-c_1S_{\omega\xi}+c_2C_{\omega\xi}+c_3S^h_{\omega\xi}+c_4C^h_{\omega\xi})}
\displaystyle{\phi''(\xi)=\omega^2(-c_1C_{\omega\xi}-c_2S_{\omega\xi}+c_3C^h_{\omega\xi}+c_4S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(c_1S_{\omega\xi}-c_2C_{\omega\xi}+c_3S^h_{\omega\xi}+c_4C^h_{\omega\xi})}

\displaystyle{ \begin{array}{|c|c|c|c|}\hline & free & pinned & fixed \\\hline free & 1^\circ & & \\\hline pinned & 3^\circ,3'^\circ & 2^\circ & \\\hline fixed & 6^\circ,6'^\circ & 5^\circ & 4^\circ \\\hline \end{array} }

1^\circ 自由・自由支持の場合 （\beta=0）

\displaystyle{\phi''(0)=\phi'''(0)=0, \phi''(1)=\phi'''(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{-c_1+c_3=0, -c_2+c_4=0 \Rightarrow c_3=c_1, c_4=c_2}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi''(\xi)=\omega^2(-c_1C_{\omega\xi}-c_2S_{\omega\xi}+c_1C^h_{\omega\xi}+c_2S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(c_1S_{\omega\xi}-c_2C_{\omega\xi}+c_1S^h_{\omega\xi}+c_2C^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1(C_{\omega}-C^h_{\omega})+c_2(S_{\omega}-S^h_{\omega})=0}
\displaystyle{c_1(S_{\omega}+S^h_{\omega})-c_2(C_{\omega}-C^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(C_{\omega}-C^h_{\omega})^2+(S^2_{\omega}-S^{h2}_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{1-C_{\omega}C^h_{\omega}=0}}

2^\circ ピン・ピン支持の場合

\displaystyle{\phi(0)=\phi''(0)=0,\ \phi(1)=\phi''(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, -c_1+c_3=0 \Rightarrow c_1=c_3=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi(\xi)=c_2S_{\omega\xi}+c_4S^h_{\omega\xi}}
\displaystyle{\phi''(\xi)=\omega^2(-c_2S_{\omega\xi}+c_4S^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_2S_{\omega}+c_4S^h_{\omega}=0, -c_2S_{\omega}+c_4S^h_{\omega}=0 \Rightarrow c_4=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{S_{\omega}=0}}

3^\circ ピン・自由支持の場合 （\beta=0）

\displaystyle{\phi(0)=\phi''(0)=0, \phi''(1)=\phi'''(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, -c_1+c_3=0 \Rightarrow c_1=c_3=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi''(\xi)=\omega^2(-c_2S_{\omega\xi}+c_4S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(-c_2C_{\omega\xi}+c_4C^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{-c_2S_{\omega}+c_4S^h_{\omega}=0}
\displaystyle{-c_2C_{\omega}+c_4C^h_{\omega}=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{S_{\omega}C^h_{\omega}-C_{\omega}S^h_{\omega}=0}}

3'^\circ ピン・質点付自由支持の場合

\displaystyle{\phi(0)=\phi''(0)=0, \phi''(1) =0,\ \phi'''(1) =-\beta\omega^4\phi(1)}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, -c_1+c_3=0 \Rightarrow c_1=c_3=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi''(\xi)=\omega^2(-c_2S_{\omega\xi}+c_4S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(-c_2C_{\omega\xi}+c_4C^h_{\omega\xi})}
\displaystyle{\phi(\xi)=c_2S_{\omega\xi}+c_4S^h_{\omega\xi}}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{-c_2S_{\omega}+c_4S^h_{\omega}=0}
\displaystyle{c_2(-C_{\omega}+\beta\omega S_{\omega})+c_4(C^h_{\omega}+\beta\omega S^h_{\omega})=0}
\displaystyle{\cdot\Downarrow}
\displaystyle{-S_{\omega}C^h_{\omega}-\beta\omega S_{\omega}S^h_{\omega}+C_{\omega}S^h_{\omega}-\beta\omega S_{\omega}S^h_{\omega}=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{(S_{\omega}C^h_{\omega}-C_{\omega}S^h_{\omega})+2\beta\omega S_{\omega}S^h_{\omega}=0}}

4^\circ 固定・固定支持の場合

\displaystyle{\phi(0)=\phi'(0)=0,\ \phi(1)=\phi'(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, c_2+c_4=0 \Rightarrow c_3=-c_1, c_4=-c_2}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi(\xi)=c_1C_{\omega\xi}+c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi}}
\displaystyle{\phi'(\xi)=\omega(-c_1S_{\omega\xi}+c_2C_{\omega\xi}-c_1S^h_{\omega\xi}-c_2C^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1(C_{\omega}-C^h_{\omega})+c_2(S_{\omega}-S^h_{\omega})=0}
\displaystyle{-c_1(S_{\omega}+S^h_{\omega})+c_2(C_{\omega}-C^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(C_{\omega}-C^h_{\omega})^2+(S^2_{\omega}-S^{h2}_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{1-C_{\omega}C^h_{\omega}=0}}

5^\circ 固定・ピン支持の場合

\displaystyle{\phi(0)=\phi'(0)=0,\ \phi(1)=\phi''(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, c_2+c_4=0 \Rightarrow c_3=-c_1, c_4=-c_2}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi(\xi)=c_1C_{\omega\xi}+c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi}}
\displaystyle{\phi''(\xi)=\omega^2(-c_1C_{\omega\xi}-c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1(C_{\omega}-C^h_{\omega})+c_2(S_{\omega}-S^h_{\omega})=0}
\displaystyle{c_1(C_{\omega}+C^h_{\omega})+c_2(S_{\omega}+S^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(C_{\omega}-C^h_{\omega})(S_{\omega}+S^h_{\omega})-(C_{\omega}+C^h_{\omega})(S_{\omega}-S^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{S_{\omega}C^h_{\omega}-C_{\omega}S^h_{\omega}=0}}

6^\circ 固定・自由支持の場合 （\beta=0）

\displaystyle{\phi(0)=\phi'(0)=0,\ \phi''(1) =\phi'''(1) =0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, c_2+c_4=0 \Rightarrow c_3=-c_1, c_4=-c_2}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi''(\xi)=\omega^2(-c_1C_{\omega\xi}-c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(c_1S_{\omega\xi}-c_2C_{\omega\xi}-c_1S^h_{\omega\xi}-c_2C^h_{\omega\xi})}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1(C_{\omega}+C^h_{\omega})+c_2(S_{\omega}+S^h_{\omega})=0}
\displaystyle{c_1(S_{\omega}-S^h_{\omega})-c_2(C_{\omega}+C^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(C_{\omega}+C^h_{\omega})^2+(S^2_{\omega}-S^{h2}_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{1+C^h_{\omega}C_{\omega}=0}}

6'^\circ 固定・質点付自由支持の場合

\displaystyle{\phi(0)=\phi'(0)=0,\ \phi''(1) =0,\ \phi'''(1) =-\beta\omega^4\phi(1)}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1+c_3=0, c_2+c_4=0 \Rightarrow c_3=-c_1, c_4=-c_2}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\phi''(\xi)=\omega^2(-c_1C_{\omega\xi}-c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi})}
\displaystyle{\phi'''(\xi)=\omega^3(c_1S_{\omega\xi}-c_2C_{\omega\xi}-c_1S^h_{\omega\xi}-c_2C^h_{\omega\xi})}
\displaystyle{\phi(\xi)=c_1C_{\omega\xi}+c_2S_{\omega\xi}-c_1C^h_{\omega\xi}-c_2S^h_{\omega\xi}}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{c_1(C_{\omega}+C^h_{\omega})+c_2(S_{\omega}+S^h_{\omega})=0}
\displaystyle{c_1((S_{\omega}-S^h_{\omega})+\beta\omega(C_{\omega}-C^h_{\omega}))+c_2(-(C_{\omega}+C^h_{\omega})+\beta\omega(S_{\omega}-S^h_{\omega}))=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{(C_{\omega}+C^h_{\omega})^2+(S^2_{\omega}-S^{h2}_{\omega})-\beta\omega(C_{\omega}+C^h_{\omega})(S_{\omega}-S^h_{\omega})+\beta\omega(C_{\omega}-C^h_{\omega})(S_{\omega}+S^h_{\omega})=0}
\displaystyle{\cdot\qquad\Downarrow}
\displaystyle{\boxed{(1+C^h_{\omega}C_{\omega})+\beta\omega(C_{\omega}S^h_{\omega}-S_{\omega}C^h_{\omega})=0}}

船舶のロバスト制御

投稿日時: 投稿者:

船舶のロバスト制御…Homework


図1 制御対象（船舶の回頭運動制御）

●船舶の針路制御のために、次のNomotoモデルを考えます。

\displaystyle{(1)\quad \dot{r}(t)=-\frac{1}{T(t)}r(t)+\frac{K(t)}{T(t)}\delta(t-L)+w(t)} }

ただし、r(t)=\dot{\psi}(t)かつ

\displaystyle{(2)\quad T(t)=\frac{L}{U(t)}T',\ K(t)=\frac{U(t)}{L}K' \quad (U_1\le U(t)\le U_2) }

ここで、舵の効果が出るまでの無駄時間Lを想定しています。混同の恐れはないので船長Lと同じ表記を用います。また前進速度Uの変動により、時定数Tと定常ゲインKが変動しますが、次式から舵は前進速度の２乗で効いてくることがわかります。

\displaystyle{(3)\quad \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{L}\right)\frac{1}{T'}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{L}\right)^2\frac{K'}{T'}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }

いま、ある前進速度U^*に注目しますと、次の表現もできます。

\displaystyle{(4a)\quad \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{U^*}\right)\frac{1}{T^*}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{U^*}\right)^2\frac{K^*}{T^*}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }

ただし

\displaystyle{(4b)\quad T^*=\frac{L}{U^*}T',\ K^*=\frac{U^*}{L}K' \quad (U_1\le U^*\le U_2) }

いま、次のような前進速度の変動を考えます。このとき、単位フィードバックによる閉ループ系は大きく変動します。


図2 前進速度が変動する場合の回頭運動

●ある前進速度U^*を基準にしたNOMOTOモデルを考えます。

\displaystyle{(5)\quad \left\{\begin{array}{l} \dot{\psi}(t)=r(t) \\ \displaystyle{\dot{r}(t)=-\left(\frac{U}{U^*}\right)\frac{1}{T^*}r(t) +\left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*}\delta(t)} \end{array}\right. }

これに、舵のダイナミクス

\displaystyle{(6)\quad \dot{\delta}(t)=-\frac{1}{T_a}\delta(t)+\frac{K_a}{T_a}u(t) }

を接続した状態空間表現を次のように表します。

\displaystyle{(7a)\quad \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{\delta}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & -\left(\frac{U}{U^*}\right)\frac{1}{T^*} & \left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*} \\ 0 & 0 & -\frac{1}{T_a} \end{array}\right] }_{A(U,U^2)} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \frac{K_a}{T_a} \end{array}\right] }_{B} u(t) }
\displaystyle{(7b)\quad \underbrace{ \psi(t) }_{y(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} }

ここで、次のようなポリトピック表現ができることに注目します（本式は未発表です）。

\displaystyle{(8a)\quad A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3 }

ただし

\displaystyle{(8b)\quad \begin{array}{l} A_1=A(U_1,U_1^2)\\ A_2=A(U_2,U_2^2)\\ A_3=A(\frac{U_1+U_2}{2},U_1U_2) \end{array} }

また、U_3=\frac{U_1+U_2}{2}を定義して

\displaystyle{(8c)\quad \begin{array}{l} p_1(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U-U_3 & U_2-U_3 \\ U^2-U_1U_2 & U_2^2-U_1U_2 \\ \end{array}\right]\\ p_2(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_3 & U-U_3 \\ U_1^2-U_1U_2 & U^2-U_1U_2 \\ \end{array}\right]\\ p_3(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_2 & U_2-U \\ U_1^2-U_2^2 & U_2^2-U^2 \\ \end{array}\right]\\ p_0=\det \left[\begin{array}{cc} U_1-U_2 & U_2-U_3 \\ U_1^2-U_2^2 & U_2^2-U_1U_2 \\ \end{array}\right]\\ p_1(U,U^2)+p_2(U,U^2)+p_3(U,U^2)=1 \end{array} }

これは端点(U_1,U_1^2)、(U_2,U_2^2)と仮想的な端点(U_3,U_1U_2)が作る３角形領域における変動をカバーすることができます。


図3 船舶のパラメータ変動凸領域

●これに基づいて、次のようなLPV制御の仕組みを考えます。


図4 船舶のLPV制御の仕組み

まず次の２ポートシステムを構成します。

\displaystyle{(9a)\quad P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A(U,U^2)& 0 \\ -C & 0 \end{array}\right] }_{A(U,U^2)} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_1} r + \underbrace{\left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc} 0 &\omega_I\\ \omega_DCA(U,U^2) & 0 \end{array}\right] }_{C_1} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{11}} r + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_DCB \end{array}\right] }_{D_{12}} u\\ \underbrace{ \left[\begin{array}{c} y \\ x_I \end{array}\right] }_{y_2} = \underbrace{ \left[\begin{array}{cc} C & 0\\ 0 & 1 \end{array}\right] }_{C_2} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{21}} r \end{array}\right.}

\displaystyle{(9b)\quad A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3 }

次に、これに対して、次のコントローラを設計します。

\displaystyle{(10a)\quad K_0: \left\{\begin{array}{l} \dot{x}_K=A_K(U,U^2)x_K+ \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & B_K^{(2)}(U,U^2) \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \\ u=C_K(U,U^2)x_K + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \end{array}\right.}

\displaystyle{(10b)\quad \begin{array}{l} A_K(U,U^2)=p_1(U,U^2)A_{K1}+p_2(U,U^2)A_{K2}+p_3(U,U^2)A_{K3}\\ B_K(U,U^2)=p_1(U,U^2)B_{K1}+p_2(U,U^2)B_{K2}+p_3(U,U^2)B_{K3}\\ C_K(U,U^2)=p_1(U,U^2)C_{K1}+p_2(U,U^2)C_{K2}+p_3(U,U^2)C_{K3}\\ D_K(U,U^2)=p_1(U,U^2)D_{K1}+p_2(U,U^2)D_{K2}+p_3(U,U^2)D_{K3} \end{array} }

これから、次の出力フィードバックを得ます。

\displaystyle{(11)\quad K: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_K \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A_K(U,U^2) & B_K^{(2)}(U,U^2) \\ 0 & 0 \end{array}\right] }_{A_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right]\\ + \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & 0\\ -1& 1 \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \\ u= \underbrace{ \left[\begin{array}{cc} C_K(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{C_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right]\\ + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & 0 \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \end{array}\right. }

●いま、U=U^*として設計したH_\infty制御による閉ループ系のシミュレーションを次に示します。３種類の速度変動に応じてばらつきか見られます。



図5 船舶のLTI制御のシミュレーションと性能評価

●これに対して、速度変動に応じたLPV制御による閉ループ系のシミュレーションを次に示します。H_\infty制御の場合のばらつきが抑制されていることが分かります。




図6 船舶のLOP制御のシミュレーションと性能評価

演習B63…Flipped Classroom
1^\circ 図2は次のプログラムを用いて得られている。これを参考にして、MATLAB用プログラムを作成しなさい。

MATLAB
SCILAB 
//ship7.sce
//-----
T1=118; T2=7.8; T3=18.5; Tship=T1+T2-T3; Kship=0.185; tL=6;
Tdelta=1; Kdelta=1;
//
Us=7.7; U1=Us*0.5; U2=Us*1.5; U3=(U1+U2)/2; U32=U1*U2;
a0=[ 0  1          0                 ;
     0  0          0                 ;
     0  0         -1/Tdelta          ];
a1=[ 0  0          0                 ;
     0 -1/Us/Tship 0                 ;
     0  0          0                 ];
a2=[ 0  0          0                 ;
     0  0          1/Us^2*Kship/Tship;
     0  0          0                 ];
A1=a0+U1*a1+U1^2*a2; 
A2=a0+U2*a1+U2^2*a2;  
A3=a0+U3*a1+U32*a2;
B0=[0;0;Kdelta/Tdelta];
//-----
function Q=interp3(P1,P2,P3,P)
 x1=P1(1); x2=P2(1); x3=P3(1); x=P(1);
 y1=P1(2); y2=P2(2); y3=P3(2); y=P(2);
 alpha=((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
 Q(1) =((x -x2)*(y2-y3)-(x2-x3)*(y -y2))/alpha;
 Q(2) =((x1-x3)*(y -y3)-(x -x3)*(y1-y3))/alpha;
 Q(3) =((x1-x2)*(y2-y )-(x2-x )*(y1-y2))/alpha;
endfunction
//
function Ut=U(t,ID)
 if ID==0, Ut=Us,
 else if t<=ID*100, Ut=Us-(Us-U1)/(ID*100)*t, else Ut=U1, end, end
endfunction
//
function dxG=fG(t,xG1), dxG=AG(t)*xG1+B0*ut, endfunction
//
function AGt=AG(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 AGt=p1*A1+p2*A2+p3*A3;
endfunction
//-----
t0=0; t1=300; nt=1500; td=(t1-t0)/nt; t=t0:td:t1; iL=tL/td;
v0=0; v1=0; v2=0; v3=0; 
for i=1:nt
 v0=[v0 U(t(i),0)]; //constant velocity
 v1=[v1 U(t(i),1)]; //rapid velocity variation
 v2=[v2 U(t(i),2)]; //medium velocity variation
 v3=[v3 U(t(i),3)]; //slow velocity variation
end
clf(0)
scf(0);subplot(211),plot(t,v0,'b',t,v1,'r',t,v2,'m',t,v3,'k'),mtlb_grid,mtlb_axis([t0 t1 0 10]);
legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
//
clf(1)
xG=zeros(3,1); y=0; u=0; ID=0
for i=1:nt
 if i<=iL, yt=0;         ut=0;
 else      yt=y(:,i-iL); ut=1-yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(0);subplot(212),plot(t,y,'b');
//
xG=zeros(3,1); y=0; u=0; ID=1
for i=1:nt
 if i<=iL, yt=0;         ut=0;
 else      yt=y(:,i-iL); ut=1-yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(0);subplot(212),plot(t,y,'r')
//
xG=zeros(3,1); y=0; u=0; ID=2
for i=1:nt
 if i<=iL, yt=0;         ut=0;
 else      yt=y(:,i-iL); ut=1-yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(0);subplot(212),plot(t,y,'m')
//
xG=zeros(3,1); y=0; u=0; ID=3
for i=1:nt
 if i<=iL, yt=0;         ut=0;
 else      yt=y(:,i-iL); ut=1-yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(0);subplot(212),plot(t,y,'k'),mtlb_grid,mtlb_axis([t0 t1 -1 3]);
legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
//-----
//eof

2^\circ 図5は次のプログラムを用いて得られている。これを参考にして、MATLAB用プログラムを作成しなさい。

MATLAB
SCILAB 
//ship8.sce
//-----
function [LME,LMI,OBJ]=synlmiof5(YLIST)
 [gam,R,S,AK,BK,CK,DK]=YLIST(:);
 LME1=R-R';
 LME2=S-S';
 LME=list(LME1,LME2);
 LMI0=[R eye(A);eye(A) S];
 AW=[A*R+B2*CK A+B2*DK*C2
     AK S*A+BK*C2];
 LMI1=-(AW+AW'+2*alpha*LMI0);
 LMI2=-[-r*LMI0  AW;
        AW'      -r*LMI0];
 LMI3=-[ sin(th)*AW   cos(th)*AW;
        -cos(th)*AW   sin(th)*AW];
 LMI3=LMI3+LMI3';
 BW=[B1+B2*DK*D21
     S*B1+BK*D21];
 CW=[C1*R+D12*CK C1+D12*DK*C2];
 [p,m]=size(D11);
 LMI4=-[AW+AW' BW CW'
        BW' -gam*eye(m,m) (D11+D12*DK*D21)'
        CW D11+D12*DK*D21 -gam*eye(p,p)];
 [p,m]=size(DK);
 LMI5=-[-1e0*eye(m,m) DK'
        DK -1e0*eye(p,p)];
 LMI=list(LMI0,LMI1,LMI2,LMI3,LMI4,LMI5);
 OBJ=gam;
endfunction
//-----
function [AK,BK,CK,DK]=syncont(R,S,ak,bk,ck,dk)
 [u,sd,v]=svd(eye()-S*R); Ni=sqrt(sd)\u'; Mti=v/sqrt(sd); 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti;
 DK=dk; 
endfunction
//=====
T1=118; T2=7.8; T3=18.5; Tship=T1+T2-T3; Kship=0.185;
Tdelta=1; Kdelta=1; tL=9; wD=tL; wI=0.01;
//
A0=[0 1 0;0 -1/Tship Kship/Tship;0 0 -1/Tdelta]; B0=[0;0;Kdelta/Tdelta]; 
C0=[1 0 0];
A=[A0 zeros(3,1);-C0 0]; B1=[zeros(3,1);1]; B2=[B0;0];
C1=[zeros(1,3) wI;wD*C0*A0 0;zeros(1,4)]; D11=zeros(3,1); D12=[0;wD*C0*B0;1];
C2=[C0 0;zeros(1,3) 1]; D21=[0;0]; D22=[0;0];
//-----
alpha=0.01; r=10; th=%pi/16;
gam0=100; R0=eye(4,4); S0=eye(4,4);
AK0=-eye(4,4); BK0=ones(4,2); CK0=ones(1,4); DK0=ones(1,2);
YLIST0=list(gam0,R0,S0,AK0,BK0,CK0,DK0);
YLIST=lmisolver(YLIST0,synlmiof5);
[gam,R,S,ak,bk,ck,dk]=YLIST(:); 
//
[AK,BK,CK,DK]=syncont(R,S,ak,bk,ck,dk)
plK=spec(AK)
ACL=[A+B2*DK*C2 B2*CK;BK*C2 AK];
BCL=[B1+B2*DK*D21;
     BK*D21];
CCL=[C2(1,:) zeros(1,4)];
plCL=spec(ACL)
//-----
AK=[AK BK(:,2);zeros(1,5)];
BK=[BK(:,1) zeros(4,1); -1 1];
CK=[CK DK(:,2)];
DK=[DK(:,1) 0];
AL=[A B2*CK;zeros(5,4) AK];
BL=[B2*DK; BK];
CL=[C2 D22*CK];
DL=D22*DK;
//-----
clf(0),clf(1)
w=logspace(-3,0,100); nw=length(w);
g=freq(A0,B0,C0,%i*w);
for i=1:nw, ga(i)=20*log10(norm(g(:,i))); end
scf(0);plot2d(w,ga,logflag='ln')
gws=20*log10(abs(wI./(%i*w))); gwt=20*log10(abs(wD^(-1)./(%i*w)));
scf(0);plot2d(w,gws,logflag='ln')
scf(0);plot2d(w,gwt,logflag='ln')
g=freq(AL,BL(:,1),CL(1,:),DL(1,1),%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(0);plot2d(w,ga,logflag='ln'),mtlb_grid
g=freq(ACL,BCL,CCL,%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(1);plot2d(w,ga,logflag='ln')
for i=1:nw, ga(i)=20*log10(abs(1-g(:,i))); end
scf(1);plot2d(w,ga,logflag='ln')
scf(1);plot2d(w,-gws,logflag='ln')
scf(1);plot2d(w,gwt,logflag='ln'),mtlb_grid
//-----
clf(2)
AG=[0 1;0 -1/Tship]; BG=[0 0;Kship/Tship 1]; CG=[1 0];
function dxG=fG(t,xG),dxG=AG*xG+BG*ut, endfunction
function dxK=fK(t,xK),dxK=AK*xK+BK*yt, endfunction
t0=0; t1=300; nt=1500; td=(t1-t0)/nt; t=t0:td:t1; iL=tL/td;
//
xG=zeros(2,1); y=0; xK=zeros(5,1); u=0; w=0
for i=1:nt
 if i<=iL, yt=[0;0];         ut=[0;w];
 else      yt=[y(:,i-iL);1]; ut=[CK*xK(:,i)+DK*yt;w]; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)]; 
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y CG*xG(:,i+1)]; u=[u ut(1)];
end
scf(2);subplot(211),plot(t,y,'r')
scf(2);subplot(212),plot(t,u,'r')
//
xG=zeros(2,1); y=0; xK=zeros(5,1); u=0; w=0.5*1e-3
for i=1:nt
 if i<=iL, yt=[0;0];         ut=[0;w];
 else      yt=[y(:,i-iL);1]; ut=[CK*xK(:,i)+DK*yt;w]; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)]; 
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y CG*xG(:,i+1)]; u=[u ut(1)];
end
scf(2);subplot(211),plot(t,y,'r--'),mtlb_grid,mtlb_axis([t0 t1 -0.5 1.5]);
scf(2);subplot(212),plot(t,u,'r--'),mtlb_grid,mtlb_axis([t0 t1 -1 1]);
//=====
Us=7.7; U1=Us*0.5; U2=Us*1.5; U3=(U1+U2)/2; U32=U1*U2;
a0=[ 0  1          0                 ;
     0  0          0                 ;
     0  0         -1/Tdelta          ];
a1=[ 0  0          0                 ;
     0 -1/Us/Tship 0                 ;
     0  0          0                 ];
a2=[ 0  0          0                 ;
     0  0          1/Us^2*Kship/Tship;
     0  0          0                 ];
A1=a0+U1*a1+U1^2*a2; 
A2=a0+U2*a1+U2^2*a2;  
A3=a0+U3*a1+U32*a2;
//-----
function Q=interp3(P1,P2,P3,P)
 x1=P1(1); x2=P2(1); x3=P3(1); x=P(1);
 y1=P1(2); y2=P2(2); y3=P3(2); y=P(2);
 alpha=((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
 Q(1) =((x -x2)*(y2-y3)-(x2-x3)*(y -y2))/alpha;
 Q(2) =((x1-x3)*(y -y3)-(x -x3)*(y1-y3))/alpha;
 Q(3) =((x1-x2)*(y2-y )-(x2-x )*(y1-y2))/alpha;
endfunction
//
function Ut=U(t,ID)
 if ID==0, Ut=Us,
 else if t<=ID*100, Ut=Us-(Us-U1)/(ID*100)*t, else Ut=U1, end, end
endfunction
//
function dxG=fG(t,xG), dxG=AG(t)*xG+B0*ut, endfunction
//
function AGt=AG(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 AGt=p1*A1+p2*A2+p3*A3;
endfunction
//-----
clf(3)
t0=0; t1=300; nt=1500; td=(t1-t0)/nt; t=t0:td:t1; iL=tL/td;
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=0
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK*xK(:,i)+DK*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(3);subplot(211),plot(t,y,'b');
scf(3);subplot(212),plot(t,u,'b');
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=1
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK*xK(:,i)+DK*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(3);subplot(211),plot(t,y,'r')
scf(3);subplot(212),plot(t,u,'r')
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=2
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK*xK(:,i)+DK*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(3);subplot(211),plot(t,y,'m')
scf(3);subplot(212),plot(t,u,'m')
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=3
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK*xK(:,i)+DK*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(3);subplot(211),plot(t,y,'k'),mtlb_grid,mtlb_axis([t0 t1 -0.5 1.5]);
//legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
scf(3);subplot(212),plot(t,u,'k'),mtlb_grid,mtlb_axis([t0 t1 -1 1]);
legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
//-----
//eof

3^\circ 図6は次のプログラムを用いて得られている。これを参考にして、MATLAB用プログラムを作成しなさい。

MATLAB
SCILAB 
//ship9.sce
//-----
function [LME,LMI,OBJ]=synlmi(YLIST)
 [gam,R,S,AK1,BK1,CK1,DK1,AK2,BK2,CK2,DK2,AK3,BK3,CK3,DK3]=YLIST(:);
 LME1=R-R'; 
 LME2=S-S'; 
 LME=list(LME1,LME2);
 LMI0=[R eye(A1);eye(A1) S];
 AW1=[A1*R+B2*CK1 A1+B2*DK1*C2;
      AK1 S*A1+BK1*C2];
 LMI11=-(AW1+AW1'+2*alpha*LMI0);
 LMI21=-[-r*LMI0 AW1;AW1 -r*LMI0];
 LMI31=-[sin(th)*AW1 cos(th)*AW1;-cos(th)*AW1 sin(th)*AW1];
 LMI31=LMI31+LMI31';
 BW=[B1+B2*DK1*D21;
     S*B1+BK1*D21];
 CW1=[C11*R+D121*CK1 C11+D121*DK1*C2];
 [p,m]=size(D11);
 LMI41=-[AW1+AW1' BW CW1'
         BW' -gam*eye(m,m) (D11+D121*DK1*D21)'
         CW1  D11+D121*DK1*D21 -gam*eye(p,p)];
 [p,m]=size(DK1);
 LMI51=-[-1e2*eye(m,m) DK1';DK1 -1e2*eye(p,p)];
//
 AW2=[A2*R+B2*CK2 A2+B2*DK2*C2;
      AK2 S*A2+BK2*C2];
 LMI12=-(AW2+AW2'+2*alpha*LMI0);
 LMI22=-[-r*LMI0 AW2;AW2 -r*LMI0];
 LMI32=-[sin(th)*AW2 cos(th)*AW2;-cos(th)*AW2 sin(th)*AW2];
 LMI32=LMI32+LMI32';
 BW=[B1+B2*DK2*D21;
     S*B1+BK2*D21];
 CW2=[C12*R+D122*CK2 C12+D122*DK2*C2];
 [p,m]=size(D11);
 LMI42=-[AW2+AW2' BW CW2'
        BW' -gam*eye(m,m) (D11+D122*DK2*D21)'
        CW2  D11+D122*DK2*D21 -gam*eye(p,p)];
 [p,m]=size(DK2);
 LMI52=-[-1e2*eye(m,m) DK2';DK2 -1e2*eye(p,p)];
//
 AW3=[A3*R+B2*CK3 A3+B2*DK3*C2;
      AK3 S*A3+BK3*C2];
 LMI13=-(AW3+AW3'+2*alpha*LMI0);
 LMI23=-[-r*LMI0 AW3;AW3 -r*LMI0];
 LMI33=-[sin(th)*AW3 cos(th)*AW3;-cos(th)*AW3 sin(th)*AW3];
 LMI33=LMI33+LMI33';
 BW=[B1+B2*DK3*D21;
     S*B1+BK3*D21];
 CW3=[C13*R+D123*CK3 C13+D123*DK3*C2];
 [p,m]=size(D11);
 LMI43=-[AW3+AW3' BW CW3'
        BW' -gam*eye(m,m) (D11+D123*DK3*D21)'
        CW3  D11+D123*DK3*D21 -gam*eye(p,p)];
 [p,m]=size(DK3);
 LMI53=-[-1e2*eye(m,m) DK3';DK3 -1e2*eye(p,p)];
 LMI=list(LMI0,LMI11,LMI21,LMI31,LMI41,LMI51,LMI12,LMI22,LMI32,LMI42,LMI52,LMI13,LMI23,LMI33,LMI43,LMI53);
 OBJ=gam;
endfunction
//-----
function [AK,BK,CK,DK]=syncont(R,S,ak,bk,ck,dk)
 [u,sd,v]=svd(eye()-S*R); Ni=sqrt(sd)\u'; Mti=v/sqrt(sd); 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti;
 DK=dk; 
endfunction
//=====
T1=118; T2=7.8; T3=18.5; Tship=T1+T2-T3; Kship=0.185;
Tdelta=10; Kdelta=1; tL=9; 
wD1=tL; wD2=tL; wD3=tL; wI1=0.01; wI2=0.01; wI3=0.01;
//
A0=[0 1 0;0 -1/Tship Kship/Tship;0 0 -1/Tdelta]; B0=[0;0;Kdelta/Tdelta]; 
C0=[1 0 0];
//-----
Us=7.7; U1=Us*0.5; U2=Us*1.5; U3=(U1+U2)/2; U32=U1*U2;
a0=[ 0  1          0                  0;
     0  0          0                  0;
     0  0         -1/Tdelta           0;
    -1  0          0                  0];
a1=[ 0  0          0                  0;
     0 -1/Us/Tship 0                  0;
     0  0          0                  0;
     0  0          0                  0];
a2=[ 0  0          0                  0;
     0  0          1/Us^2*Kship/Tship 0;
     0  0          0                  0;
     0  0          0                  0];
A1=a0+U1*a1+U1^2*a2; 
A2=a0+U2*a1+U2^2*a2;  
A3=a0+U3*a1+U32*a2;
B1=[zeros(3,1);1]; 
B2=[B0;0];
C11=[zeros(1,3) wI1;wD1*C0*A0 0;zeros(1,4)]; 
C12=[zeros(1,3) wI2;wD2*C0*A0 0;zeros(1,4)]; 
C13=[zeros(1,3) wI3;wD3*C0*A0 0;zeros(1,4)]; 
D11=zeros(3,1); 
D121=[0;wD1*C0*B0;1];
D122=[0;wD2*C0*B0;1];
D123=[0;wD3*C0*B0;1];
C2=[C0 0;zeros(1,3) 1]; 
D21=[0;0]; 
D22=[0;0];
//-----
alpha=0.01; r=15; th=%pi/4;
gam0=100; R0=eye(4,4); S0=eye(4,4);
AK0=-eye(4,4); BK0=ones(4,2); CK0=ones(1,4); DK0=ones(1,2);
YLIST0=list(gam0,R0,S0,AK0,BK0,CK0,DK0,AK0,BK0,CK0,DK0,AK0,BK0,CK0,DK0);
YLIST=lmisolver(YLIST0,synlmi);
[gam,R,S,ak1,bk1,ck1,dk1,ak2,bk2,ck2,dk2,ak3,bk3,ck3,dk3]=YLIST(:); 
//
A=A1; [AK1,BK1,CK1,DK1]=syncont(R,S,ak1,bk1,ck1,dk1);
plK1=spec(AK1),
ACL1=[A1+B2*DK1*C2 B2*CK1;BK1*C2 AK1];
BCL1=[B1+B2*DK1*D21;BK1*D21];
CCL1=[C2(1,:) zeros(1,4)];
plCL1=spec(ACL1)
//
A=A2; [AK2,BK2,CK2,DK2]=syncont(R,S,ak2,bk2,ck2,dk2);
plK2=spec(AK2),
ACL2=[A2+B2*DK2*C2 B2*CK2;BK2*C2 AK2];
BCL2=[B2+B2*DK2*D21;BK2*D21];
CCL2=[C2(1,:) zeros(1,4)];
plCL2=spec(ACL2)
//
A=A3; [AK3,BK3,CK3,DK3]=syncont(R,S,ak3,bk3,ck3,dk3);
plK3=spec(AK3),
ACL3=[A3+B2*DK3*C2 B2*CK3;BK3*C2 AK3];
BCL3=[B1+B2*DK3*D21;BK3*D21];
CCL3=[C2(1,:) zeros(1,4)];
plCL3=spec(ACL3)
//-----
AK1=[AK1 BK1(:,2);zeros(1,5)];
BK1=[BK1(:,1) zeros(4,1); -1 1];
CK1=[CK1 DK1(:,2)];
DK1=[DK1(:,1) 0];
AL1=[A1 B2*CK1;zeros(5,4) AK1];
BL1=[B2*DK1; BK1];
CL1=[C2 D22*CK1];
DL1=D22*DK1;
clf(0),clf(1)
w=logspace(-3,0,100); nw=length(w);
g=freq(A0,B0,C0,%i*w);
for i=1:nw, ga(i)=20*log10(norm(g(:,i))); end
scf(0);plot2d(w,ga,logflag='ln')
gws=20*log10(abs(wI1./(%i*w))); gwt=20*log10(abs(wD1^(-1)./(%i*w)));
scf(0);plot2d(w,gws,logflag='ln')
scf(0);plot2d(w,gwt,logflag='ln')
g=freq(AL1,BL1(:,1),CL1(1,:),DL1(1,1),%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(0);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
g=freq(ACL1,BCL1,CCL1,%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(1);plot2d(w,ga,logflag='ln')
for i=1:nw, ga(i)=20*log10(abs(1-g(:,i))); end
scf(1);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
scf(1);plot2d(w,-gws,logflag='ln')
scf(1);plot2d(w,gwt,logflag='ln')
//
clf(2),clf(3)
AK2=[AK2 BK2(:,2);zeros(1,5)];
BK2=[BK2(:,1) zeros(4,1); -1 1];
CK2=[CK2 DK2(:,2)];
DK2=[DK2(:,1) 0];
AL2=[A2 B2*CK1;zeros(5,4) AK1];
BL2=[B2*DK1; BK1];
CL2=[C2 D22*CK1];
DL2=D22*DK1;
w=logspace(-3,0,100); nw=length(w);
g=freq(A0,B0,C0,%i*w);
for i=1:nw, ga(i)=20*log10(norm(g(:,i))); end
scf(2);plot2d(w,ga,logflag='ln')
gws=20*log10(abs(wI2./(%i*w))); gwt=20*log10(abs(wD2^(-1)./(%i*w)));
scf(2);plot2d(w,gws,logflag='ln')
scf(2);plot2d(w,gwt,logflag='ln')
g=freq(AL2,BL2(:,1),CL2(1,:),DL2(1,1),%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(2);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
g=freq(ACL2,BCL2,CCL2,%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(3);plot2d(w,ga,logflag='ln')
for i=1:nw, ga(i)=20*log10(abs(1-g(:,i))); end
scf(3);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
scf(3);plot2d(w,-gws,logflag='ln')
scf(3);plot2d(w,gwt,logflag='ln')
//
clf(4),clf(5)
AK3=[AK3 BK3(:,2);zeros(1,5)];
BK3=[BK3(:,1) zeros(4,1); -1 1];
CK3=[CK3 DK3(:,2)];
DK3=[DK3(:,1) 0];
AL3=[A3 B2*CK1;zeros(5,4) AK1];
BL3=[B2*DK1; BK1];
CL3=[C2 D22*CK1];
DL3=D22*DK1;
w=logspace(-3,0,100); nw=length(w);
g=freq(A0,B0,C0,%i*w);
for i=1:nw, ga(i)=20*log10(norm(g(:,i))); end
scf(4);plot2d(w,ga,logflag='ln')
gws=20*log10(abs(wI3./(%i*w))); gwt=20*log10(abs(wD3^(-1)./(%i*w)));
scf(4);plot2d(w,gws,logflag='ln')
scf(4);plot2d(w,gwt,logflag='ln')
g=freq(AL3,BL3(:,1),CL3(1,:),DL3(1,1),%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(4);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
g=freq(ACL3,BCL3,CCL3,%i*w); 
for i=1:nw, ga(i)=20*log10(abs(g(:,i))); end
scf(5);plot2d(w,ga,logflag='ln')
for i=1:nw, ga(i)=20*log10(abs(1-g(:,i))); end
scf(5);plot2d(w,ga,logflag='ln'),mtlb_grid,mtlb_axis([10^(-3) 10^0 -80 60])
scf(5);plot2d(w,-gws,logflag='ln')
scf(5);plot2d(w,gwt,logflag='ln')
//return
//=====
a0=[ 0  1          0                 ;
     0  0          0                 ;
     0  0         -1/Tdelta          ];
a1=[ 0  0          0                 ;
     0 -1/Us/Tship 0                 ;
     0  0          0                 ];
a2=[ 0  0          0                 ;
     0  0          1/Us^2*Kship/Tship;
     0  0          0                 ];
A1=a0+U1*a1+U1^2*a2; 
A2=a0+U2*a1+U2^2*a2;  
A3=a0+U3*a1+U32*a2;
//-----
function Q=interp3(P1,P2,P3,P)
 x1=P1(1); x2=P2(1); x3=P3(1); x=P(1);
 y1=P1(2); y2=P2(2); y3=P3(2); y=P(2);
 alpha=((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
 Q(1) =((x -x2)*(y2-y3)-(x2-x3)*(y -y2))/alpha;
 Q(2) =((x1-x3)*(y -y3)-(x -x3)*(y1-y3))/alpha;
 Q(3) =((x1-x2)*(y2-y )-(x2-x )*(y1-y2))/alpha;
endfunction
//
function Ut=U(t,ID)
 if ID==0, Ut=Us,
 else
  if t<=ID*100, Ut=Us-(Us-U1)/(ID*100)*t, else Ut=U1, end
 end
endfunction
//
function dxG=fG(t,xG), dxG=AG(t)*xG+B0*ut, endfunction
//
function AGt=AG(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 AGt=p1*A1+p2*A2+p3*A3;
endfunction
//
function dxK=fK(t,xK), 
 dxK=AK(t)*xK+BK(t)*yt, 
endfunction
//
function AKt=AK(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 AKt=p1*AK1+p2*AK2+p3*AK3;
endfunction
//
function BKt=BK(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 BKt=p1*BK1+p2*BK2+p3*BK3;
endfunction
//
function CKt=CK(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 CKt=p1*CK1+p2*CK2+p3*CK3;
endfunction
//
function DKt=DK(t), 
 P1=[U1;U1^2]; P2=[U2;U2^2]; P3=[U3;U3^2]; Ut=U(t,ID); P=[Ut;Ut^2]; 
 Q=interp3(P1,P2,P3,P); p1=Q(1); p2=Q(2); p3=Q(3); 
 DKt=p1*DK1+p2*DK2+p3*DK3;
endfunction
//-----
clf(6)
t0=0; t1=300; nt=1500; td=(t1-t0)/nt; t=t0:td:t1; iL=tL/td;
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=0;
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK(t(i))*xK(:,i)+DK(t(i))*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(6);subplot(211),plot(t,y,'b');
scf(6);subplot(212),plot(t,u,'b');
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=1;
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK(t(i))*xK(:,i)+DK(t(i))*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(6);subplot(211),plot(t,y,'r');
scf(6);subplot(212),plot(t,u,'r');
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=2;
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK(t(i))*xK(:,i)+DK(t(i))*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(6);subplot(211),plot(t,y,'m');
scf(6);subplot(212),plot(t,u,'m');
//
xG=zeros(3,1); y=0; xK=zeros(5,1); u=0; ID=3;
for i=1:nt
 if i<=iL, yt=[0;1];         ut=0;
 else      yt=[y(:,i-iL);1]; ut=CK(t(i))*xK(:,i)+DK(t(i))*yt; end
 xG=[xG ode(xG(:,i),t(i),t(i+1),fG)];
 xK=[xK ode(xK(:,i),t(i),t(i+1),fK)];
 y=[y xG(1,i+1)]; u=[u ut];
end
scf(6);subplot(211),plot(t,y,'k'),mtlb_grid,mtlb_axis([t0 t1 -0.5 1.5]);
//legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
scf(6);subplot(212),plot(t,u,'k'),mtlb_grid,mtlb_axis([t0 t1 -1 1]);
legend(['const.vel.';'rapid vel.var.';'medium vel.var.';'slow vel.var.'])
//-----
//eof

問題設定

投稿日時: 投稿者:

問題設定…Homework

[1] H_\infty制御系を設計する際の２ポートシステムの構成法について述べます。


図1 ２ポートシステム

次の制御対象を考えます。

\displaystyle{(1)\quad G: \left\{\begin{array}{l} \dot{x}=Ax+Bu\\ y=Cx \end{array}\right.}

これに対する評価出力として、次を考えます。

\displaystyle{(2)\quad W_S: \left\{\begin{array}{l} \dot{x}_I=r-y\\ y_{11}=\omega_Ix_I \end{array}\right. }

\displaystyle{(3)\quad W_T: y_{12}=\omega_D\dot{y}=\omega_DCAx+\omega_DCBu}


図2 積分器を含める場合の２ポートシステム

このとき、次の２ポートシステムを構成します。

\displaystyle{(4)\quad \boxed{P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A& 0 \\ -C & 0 \end{array}\right] }_{A} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_1} r + \underbrace{\left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc} 0 &\omega_I\\ \omega_DCA & 0 \end{array}\right] }_{C_1} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{11}} r + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_DCB \end{array}\right] }_{D_{12}} u\\ \underbrace{ \left[\begin{array}{c} y \\ x_I \end{array}\right] }_{y_2} = \underbrace{ \left[\begin{array}{cc} C & 0\\ 0 & 1 \end{array}\right] }_{C_2} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{21}} r \end{array}\right.} }

これに対して、次のコントローラを設計します。

\displaystyle{(5)\quad K_0: \left\{\begin{array}{l} \dot{x}_K=A_Kx_K+ \underbrace{ \left[\begin{array}{cc} B_{K1} & B_{K2} \end{array}\right] }_{B_K} \left[\begin{array}{c} y \\ x_I \end{array}\right] \\ u=C_Kx_K + \underbrace{ \left[\begin{array}{cc} D_{K1} & D_{K2} \end{array}\right] }_{D_K} \left[\begin{array}{c} y \\ x_I \end{array}\right] \end{array}\right.}

これによる閉ループ系

\displaystyle{(6)\quad P_{CL}: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_K \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} A+B_2D_KC_2 & B_2C_K \\ B_KC_2 & A_K \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ \left[\begin{array}{ccc} B_1+B_2D_KD_{21} \\ B_KD_{21} \end{array}\right] }_{B_{CL}} r\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{ccc} C_1+D_{12}D_KC_2 & D_{12}C_K \end{array}\right] }_{C_{CL}} \left[\begin{array}{c} x \\ x_K \end{array}\right] + \underbrace{ (D_{11}+D_{12}D_KD_{21}) }_{D_{CL}} r \end{array}\right. }

において、そのH_\inftyノルムが\gammaより小、すなわち

\displaystyle{(7)\quad \sup_{\omega\in{\rm\bf R}} \underbrace{||\hat{P}_{CL}(j\omega)||_2}_{\sigma_1(\hat{P}_{CL}(j\omega))} <\gamma }

ただし

\displaystyle{(8)\quad \hat{P}_{CL}(s)= \left[\begin{array}{c} \hat{G}_{y_{11}}(s)\\ \hat{G}_{y_{12}}(s) \end{array}\right]= C_{CL}(sI-A_{CL})^{-1}B_{CL}+D_{CL} }

となるように出力フィードバックのゲインA_K,B_K,C_K,D_Kを求める問題を考えます。

これをLMIベース設計問題として定式化すると次のようになります。

Minimize \gamma on R=R^T,S=S^T,{\cal A}_K,{\cal B}_K,{\cal C}_K,D_K subject to

\displaystyle{(9)\quad \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] >0 }

\displaystyle{(10)\quad \left[\begin{array}{ccc} \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right]+(*)^T & \left[\begin{array}{c} B_1+B_2D_KD_{21} \\ SB_1+{\cal B}_KD_{21} \end{array}\right] & (*)^T \\ (*)^T & -\gamma^2 I & (*)^T \\ \left[\begin{array}{cc} C_1R+D_{12}{\cal C}_K & C_1+D_{12}D_KC_2 \end{array}\right] & D_{11} & -I \end{array}\right]<0 }

Determine the output feedback controller A_K,B_K,C_K　by

\displaystyle{(11)\quad \begin{array}{lll} A_K&=&N^{-1}({\cal A}_K-S(A-B_2D_KC_2)R-{\cal B}_KC_2R-SB_2{\cal C}_K)M^{-T} \\ B_K&=&N^{-1}({\cal B}_K-SB_2D_K) \\ C_K&=&({\cal C}_K-D_KC_2R)M^{-T} \\ &&(I-SR=NM^T) \end{array} }

実際には、次のＬＭＩも追加することが行われます。

\displaystyle{(12)\quad \begin{array}{l} \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] +(*)^T %\nonumber\\&& +\alpha \left[\begin{array}{cc} R & I \\ I & S \end{array}\right]<0}\\ \left[\begin{array}{cc} -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] & \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right] \\ (*)^T & -r \left[\begin{array}{cc} R & I \\ I & S \end{array}\right] \end{array}\right] <0\\ \left[\begin{array}{cc} \sin\theta & \cos\theta \\ -\cos\theta & \sin\theta \end{array}\right]\otimes \left[\begin{array}{cc} AR+B_2{\cal C}_K & A+B_2D_KC_2\\ {\cal A}_K & SA+{\cal B}_KC_2 \end{array}\right]+ (*)^T <0 \end{array} }

ちなみに、コントローラK_0に積分器

\displaystyle{(13)\quad\dot{x}_I=r-y}

を接続して、制御対象に対するコントローラとして次を得ます。

\displaystyle{(14)\quad \boxed{K: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_K \\ \dot{x}_I \end{array}\right]= \left[\begin{array}{cc} A_K & B_{K2} \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \left[\begin{array}{cc} B_{K1} & 0\\ -1& 1 \end{array}\right] \left[\begin{array}{c} y \\ r \end{array}\right] \\ u= \left[\begin{array}{cc} C_K & D_{K2} \end{array}\right] \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \left[\begin{array}{cc} D_{K1} & 0 \end{array}\right] \left[\begin{array}{c} y \\ r \end{array}\right] \end{array}\right.} }

[2] 上の制御系設計は混合感度問題として定式化されていることを示すために、次式に注目します。

\displaystyle{(15)\quad \hat{P}_{CL}(s)= \left[\begin{array}{c} \hat{G}_{y_{11}}(s)\\ \hat{G}_{y_{12}}(s) \end{array}\right]= \left[\begin{array}{c} \hat{W}_S(s)\hat{S}(s)\\ \hat{W}_T(s)\hat{T}(s) \end{array}\right] }

ただし

\displaystyle{(16)\quad \hat{G}_{y_{11}}(s) =\omega_I\frac{\frac{1}{s}}{1+\hat{G}(s)\hat{K}_0(s)\frac{1}{s}} =\underbrace{\frac{\omega_I}{s}}_{\hat{W}_S(s)} \underbrace{\frac{1}{1+\hat{G}(s)\hat{K}(s)}}_{\hat{S}(s)} }
\displaystyle{(17)\quad \hat{G}_{y_{12}}(s) =\omega_Ds\frac{\hat{G}(s)\hat{K}_0(s)\frac{1}{s}}{1+\hat{G}(s)\hat{K}_0(s)\frac{1}{s}} =\underbrace{\omega_Ds}_{\hat{W}_T(s)} \underbrace{\frac{\hat{G}(s)\hat{K}(s)}{1+\hat{G}(s)\hat{K}(s)}}_{\hat{T}(s)} }

したがって、

\displaystyle{(18)\quad ||\hat{P}_{CL}(j\omega)||_2=\sigma_1(\hat{P}_{CL}(j\omega)) =\sqrt{|\hat{W}_S(j\omega)\hat{S}(j\omega)|^2+|\hat{W}_T(j\omega)\hat{T}(j\omega)|^2} }

に注意して

\displaystyle{(19)\quad \sup_{\omega\in{\rm\bf R}}||\hat{P}_{CL}(j\omega)||_2<\gamma }

を最小化することは

\displaystyle{(20)\quad \forall\omega: |\hat{W}_S(j\omega)\hat{S}(j\omega)|^2+|\hat{W}_T(j\omega)\hat{T}(j\omega)|^2<\gamma^2 }

を最小化していることになります。

演習B62…Flipped Classroom
1^\circ 水タンクにおいて、バルブに無駄時間がある場合の次の状態方程式を考えます。

\displaystyle{(21)\quad\underbrace{\frac{d}{dt}(h(t)-h^*)}_{\dot{x}(t)}=\underbrace{-\frac{k}{2S\sqrt{h^*}}}_{a}\underbrace{(h(t)-h^*)}_{x(t)}+\underbrace{\frac{1}{S}}_{b}\underbrace{(q_i(t-L)-q_i^*)}_{u(t-L)}}

これに対するH_\infty制御系を設計する次のプログラムを実行して得られる３枚の図の意味について説明せよ。

MATLAB 
%b62.m 
%-----
 clear all, close all
 S=1; k=0.1; hs=1; L=10;
 a=-k/(2*sqrt(hs)*S); b=1/S; c=1;
 wD=L; wI=0.02;
 A=[a 0;-c 0]; B1=[0;1]; B2=[b;0];
 C1=[0 wI;wD*c*a 0]; D11=[0;0]; D12=[0;wD*c*b];
 C2=[c 0;0 1]; D21=[0;0]; D22=[0;0];
%-----
 alpha=0.001; r=0.1; th=pi/3;
 LMIs=of_synlmi6(A,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th);
 cobj=zeros(1,decnbr(LMIs)); 
 cobj(1)=1; 
 [cost,xopt]=mincx(LMIs,cobj);  
 gopt=dec2mat(LMIs,xopt,1)
 R=dec2mat(LMIs,xopt,2); 
 S=dec2mat(LMIs,xopt,3);
 ak=dec2mat(LMIs,xopt,4); 
 bk=dec2mat(LMIs,xopt,5); 
 ck=dec2mat(LMIs,xopt,6); 
 dk=dec2mat(LMIs,xopt,7); 
 [u,sd,v]=svd(eye(size(A,1))-R*S); 
 sd=diag(sqrt(1./diag(sd))); 
 Ni=sd*v'; Mti=u*sd; 
 AK=Ni*(ak-S*(A-B2*dk*C2)*R-bk*C2*R-S*B2*ck)*Mti; 
 BK=Ni*(bk-S*B2*dk); 
 CK=(ck-dk*C2*R)*Mti; 
 DK=dk;
%----- 
 pl=eig(A);
 ACL=[A+B2*DK*C2 B2*CK;BK*C2 AK];
 plCL=eig(ACL)
 figure(1)
 dregion(-alpha,0,r,th,r*[-2,1,-1,1]) 
 plot(real(pl),imag(pl),'o',real(plCL),imag(plCL),'x')
%-----  
 ACL=[A+B2*DK*C2 B2*CK;BK*C2 AK];
 BCL=[B1+B2*DK*D21;BK*D21];
 CCL=[C2(1,:) zeros(1,2)];
 AK=[AK BK(:,2);zeros(1,3)];
 BK=[BK(:,1) zeros(2,1); -1 1];
 CK=[CK DK(:,2)];
 DK=[DK(:,1) 0];
%-----   
 w=logspace(-2,1,100); 
 H=exp(-L*i*w); ga_H=20*log10(abs(H));
 ga_DL=20*log10(abs(H-1));
 ga_WT=20*log10(abs(L*i*w));
 figure(2)
 w=logspace(-3,0,100); 
 semilogx(w,ga_DL,w,ga_WT),grid, legend('Delta','WT') 
 %-----  
 ga_WS=20*log10(abs(wI./(i*w))); 
 ga_WT=20*log10(abs(wD*(i*w)));
 G(:,:)=freqresp(ss(a,b,c,[]),w); G=G.'; ga_G=20*log10(abs(G));
 K(:,:)=freqresp(ss(AK,BK(:,1),CK,DK(:,1)),w); K=K.'; ga_K=20*log10(abs(K)); 
 T(:,:)=freqresp(ss(ACL,BCL(:,1),CCL(1,:),[]),w); T=T.'; ga_T=20*log10(abs(T)); 
 ga_S=20*log10(abs(1-T)); 
 figure(3)
 subplot(121),semilogx(w,ga_G+ga_K,w,ga_WS,'b',w,-ga_WT,'r'),grid,legend('L=GK','WS','1/WT')
 subplot(122),semilogx(w,ga_T,'r',w,-ga_WT,'r',w,ga_S,'b',w,-ga_WS,'b'),grid,legend('T','1/WT','S','1/WS')
%-----
function LMIs=of_synlmi6(A,B1,B2,C1,C2,D11,D12,D21,D22,alpha,r,th)
 [n,m]=size(B2);  [p,n]=size(C2);
 setlmis([]); 
 gam=lmivar(1,[1 0]);  
 [R,xxx,Rdec]=lmivar(1,[n 1]); 
 [S,xxx,Sdec]=lmivar(1,[n 1]); 
 Ak=lmivar(2,[n n]); 
 Bk=lmivar(2,[n p]); 
 Ck=lmivar(2,[m n]); 
 Dk=lmivar(2,[m p]); 
% 
 lmiRS=newlmi;
 lmiterm([lmiRS 1 1 R],A,1,'s');       %#1:R*A'+AR 
 lmiterm([lmiRS 1 1 Ck],B2,1,'s');     %#1:Ck'*B2'+B2*Ck 
 lmiterm([lmiRS 2 1 0],A');            %#1:A' 
 lmiterm([lmiRS 2 1 Ak],1,1);          %#1:Ak 
 lmiterm([lmiRS 2 1 -Dk],C2',B2');     %#1:C2'*Dk'*B2' 
 lmiterm([lmiRS 2 2 S],1,A,'s');       %#1:A'*S+S*A 
 lmiterm([lmiRS 2 2 Bk],1,C2,'s');     %#1:C2'*Bk'+Bk*C2 
 lmiterm([lmiRS 1 3 0],B1);            %#1:B1 
 lmiterm([lmiRS 1 3 Dk],B2,D21);       %#1:B2*Dk*D21 
 lmiterm([lmiRS 2 3 S],1,B1);          %#1:S*B1 
 lmiterm([lmiRS 2 3 Bk],1,D21);        %#1:Bk*D21 
 lmiterm([lmiRS 3 3 gam],-1,1);        %#1:-gam 
 lmiterm([lmiRS 4 1 R],C1,1);          %#1:C1*R 
 lmiterm([lmiRS 4 1 Ck],D12,1);        %#1:D12*Ck 
 lmiterm([lmiRS 4 2 0],C1);            %#1:C1 
 lmiterm([lmiRS 4 2 Dk],D12,C2);       %#1:D12*Dk*C2 
 lmiterm([lmiRS 4 3 0],D11);           %#1:D11 
 lmiterm([lmiRS 4 3 Dk],D12,D21);      %#1:D12*Dk*D21 
 lmiterm([lmiRS 4 4 gam],-1,1);        %#1:-gam 
%-----
 lmiPL1=newlmi; 
 lmiterm([lmiPL1 1 1 R],A,1,'s');      %#2:R*A'+A*R 
 lmiterm([lmiPL1 1 1 Ck],B2,1,'s');    %#2:Ck'*B2'+B2*Ck 
 lmiterm([lmiPL1 2 1 Ak],1,1);         %#2:Ak 
 lmiterm([lmiPL1 1 2 0],A);            %#2:A 
 lmiterm([lmiPL1 1 2 Dk],B2,C2);       %#2:B2*Dk*C2 
 lmiterm([lmiPL1 2 2 S],1,A,'s');      %#2:A'*S+S*A 
 lmiterm([lmiPL1 2 2 Bk],1,C2,'s');    %#2:C2'*Bk'+Bk*C2 
% 
 lmiterm([lmiPL1 1 1 R],2*alpha,1);    %#2:2*alpha*R
 lmiterm([lmiPL1 2 1 0],2*alpha);      %#2:2*alpha*I 
 lmiterm([lmiPL1 2 2 S],2*alpha,1);    %#2:2*alpha*S 
%-----
 lmiPL2=newlmi; 
 lmiterm([lmiPL2 1 1 R],-r,1);         %#3:-r*R 
 lmiterm([lmiPL2 2 1 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 2 2 S],-r,1);         %#3:-r*S 
% 
 lmiterm([lmiPL2 1 3 R],A,1);          %#3:A*R
 lmiterm([lmiPL2 1 3 Ck],B2,1);        %#3:B2*Ck 
 lmiterm([lmiPL2 2 3 Ak],1,1);         %#3:Ak 
 lmiterm([lmiPL2 1 4 0],A);            %#3:A 
 lmiterm([lmiPL2 1 4 Dk],B2,C2);       %#3:B2*Dk*C2 
 lmiterm([lmiPL2 2 4 S],1,A);          %#3:S*A 
 lmiterm([lmiPL2 2 4 Bk],1,C2);        %#3:Bk*C2 
% 
 lmiterm([lmiPL2 3 3 R],-r,1);         %#3:-r*R
 lmiterm([lmiPL2 4 3 0],-r);           %#3:-r*I 
 lmiterm([lmiPL2 4 4 S],-r,1);         %#3:-r*S 
%-----
 sth=sin(th); cth=cos(th);
 lmiPL3=newlmi; 
 lmiterm([lmiPL3 1 1 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R) 
 lmiterm([lmiPL3 1 1 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 2 1 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 1 2 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 1 2 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 2 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 2 2 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
% 
 lmiterm([lmiPL3 1 3 R],cth*A,1);      %#1:cth*(A*R) 
 lmiterm([lmiPL3 1 3 R],1,-cth*A');    %#1:cth*(-R*A')
 lmiterm([lmiPL3 1 3 Ck],cth*B2,1);     %#1:cth*(B*Ck) 
 lmiterm([lmiPL3 1 3 -Ck],-cth*B2',1);  %#1:cth*(-Ck'*B') 
 lmiterm([lmiPL3 2 3 Ak],cth,1);       %#4:cth*(Ak) 
 lmiterm([lmiPL3 1 4 -Ak],-cth,1);     %#4:cth*(-Ak') 
 lmiterm([lmiPL3 1 4 0],A);            %#4:cth*(A) 
 lmiterm([lmiPL3 2 3 0],-A');          %#4:cth*(-A') 
 lmiterm([lmiPL3 1 4 Dk],cth*B2,C2);   %#4:cth*(B2*Dk*C2) 
 lmiterm([lmiPL3 2 3 -Dk],-cth*C2',B2');%#4:cth*(-C2'*Dk'*B2') 
 lmiterm([lmiPL3 2 4 S],1,cth*A);      %#4:cth*(S*A) 
 lmiterm([lmiPL3 2 4 S],-cth*A',1);    %#4:cth*(-A'*S) 
 lmiterm([lmiPL3 2 4 Bk],1,cth*C2);    %#4:cth*(Bk*C2) 
 lmiterm([lmiPL3 2 4 -Bk],-cth*C2',1); %#4:cth*(-C2'*Bk') 
% 
 lmiterm([lmiPL3 3 3 R],sth*A,1,'s');  %#4:sth*(R*A'+A*R)
 lmiterm([lmiPL3 3 3 Ck],sth*B2,1,'s');%#4:sth*(Ck'*B2'+B2*Ck) 
 lmiterm([lmiPL3 4 3 Ak],sth,1);       %#4:sth*(Ak) 
 lmiterm([lmiPL3 3 4 0],sth*A);        %#4:sth*(A) 
 lmiterm([lmiPL3 3 4 Dk],sth*B2,C2);   %#4:sth*(B2*Dk*C2) 
 lmiterm([lmiPL3 4 4 S],1,sth*A,'s');  %#4:sth*(A'*S+S*A) 
 lmiterm([lmiPL3 4 4 Bk],1,sth*C2,'s');%#4:sth*(C2'*Bk'+Bk*C2) 
%-----
 posX=-newlmi;
 lmiterm([posX 1 1 R],1,1);            %#5:R 
 lmiterm([posX 2 1 0],1);              %#5:I 
 lmiterm([posX 2 2 S],1,1);            %#5:S 
%-----
 lmiDk=-newlmi;                        
 lmiterm([lmiDk 1 1 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 2 0],1e2);           %#6:1e2 
 lmiterm([lmiDk 2 1 Dk],1,1);          %#6:Dk 
%-----
 lmig=newlmi;                          
 lmiterm([lmig,1,1,gam],1,1);          %#7:gam 
 lmiterm([lmig,1,1,0],-1e3);           %#7:1e3 
 LMIs=getlmis;
end
%-----
%eof