１次系のLQ制御 | 制御系CAD

１次系のLQ制御…Homework

[1] ある平衡状態周りの挙動が状態方程式

$\displaystyle{(1)\quad \dot{x}(t)=ax(t)+bu(t) }$

で表される制御対象に対して、状態フィードバック（state feedback）

$\displaystyle{(2)\quad \boxed{u(t)=-fx(t)} }$

を考えます。ここで、定数 $f$ はゲインと呼ばれ、 $a-bf<0$ を満たすものを選びます。(2)を(1)に代入して、次式を得ます。

$\displaystyle{(3)\quad \boxed{\dot{x}(t)=(a-bf)x(t)\quad (a-bf<0)} }$

これは制御対象とコントローラ（状態フィードバック）が閉路を構成することから（図1参照）、閉ループ系（closed-loop system）と呼ばれます。ここで、物理的に興味深いのは、制御対象が漸近安定でなくても（ $a\ge0$ ）、状態フィードバックを行って閉ループ系を漸近安定（ $a-bf<0$ ）にできることです。そのようなゲイン $f$ を決定する問題を状態フィードバックによる安定化問題と言います。

図1 状態フィードバックによる閉ループ系

[2] １次系(1)に対する安定化状態フィードバック(2)による閉ループ系に対して、評価関数

$\displaystyle{(4)\quad \boxed{J=\int_0^\infty(q^2x^2(t)+r^2u^2(t))\,dt} }$

を最小化するように $f$ を決めるLQ制御問題を考えます。

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

$\displaystyle{(5)\quad x(t)=e^{(a-bf)t}x(0) }$

したがって、評価関数(4)の第１項は、次式のように計算されます。

$(6)\quad \begin{array}{lll} \displaystyle{J_x=\int_0^\infty q^2x^2(t)\,dt}\\ \displaystyle{=\int_0^\infty q^2e^{2(a-bf)t}x^2(0)\,dt}\\ \displaystyle{=q^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty}\\ \displaystyle{=\frac{q^2x^2(0)}{2(a-bf)}\left[\underbrace{e^{2(a-bf)\infty}}_{0}-\underbrace{e^{2(a-bf)0}}_{1}\right]}\\ \displaystyle{=-\frac{q^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0)} \end{array}$

これから $f$ を大きくすれば、 $J_x$ はいくらでも小さくできることが分かります。

●閉ループ系における入力の振舞い $u=-fx$ は次式で与えられます。

$\displaystyle{(7)\quad u(t)=-fe^{(a-bf)t}x(0) }$

したがって、評価関数(4)の第２項は、次式のように計算されます。

$(8)\quad \begin{array}{lll} \displaystyle{J_u=\int_0^\infty r^2u^2(t)\,dt}\\ \displaystyle{=\int_0^\infty r^2f^2e^{2(a-bf)t}x^2(0)\,dt}\\ \displaystyle{=r^2f^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty}\\ \displaystyle{=-\frac{r^2f^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0)} \end{array}$

$J_u$ が下に凸で、 $f=\frac{2b}{a}$ において最小値を取ることは、次式から分かります。

$(9)\quad \begin{array}{lll} \displaystyle{\frac{dJ_u}{df}=\frac{d}{df}\left(-\frac{r^2f^2}{2(a-bf)}x^2(0)\right)}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\left(\frac{2f}{(a-bf)}+(-1)\frac{f^2}{(a-bf)^2}(-b)\right)}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\frac{2(a-bf)f+f^2b}{(a-bf)^2}}\\ \displaystyle{=\frac{-r^2x^2(0)}{2}\frac{(2a-bf)f}{(a-bf)^2}} \end{array}$

$(10)\quad \begin{array}{lll} \displaystyle{\frac{d^2J_u}{df^2}=\frac{-r^2x^2(0)}{2}\left(\frac{2a-2bf}{(a-bf)^2}+(-2)\frac{(2a-bf)f}{2(a-bf)^3}(-b)\right)}\\ \displaystyle{=-r^2x^2(0)\left(\frac{1}{a-bf}+b\frac{2(a-bf)f+bf^2}{(a-bf)^3}\right)}\\ \displaystyle{=-r^2x^2(0)\frac{(a-bf)^2+2(a-bf)bf+b^2f^2}{(a-bf)^3}}r\\ \displaystyle{=-r^2x^2(0)\frac{(a-bf+bf)^2}{2(a-bf)^3}\}\\ \displaystyle{=-r^2x^2(0)\frac{a^2}{2(a-bf)^3}>0} \end{array}$

いま、 $b=1,q=1,r=1$ として、 $a=-1$ と $a=1$ の場合について、 $J_x$ 、 $J_u$ 、 $J=J_x+J_u$ の外形を描くと、次図のようになります。ここで、 $J$ には○で示す最小値が生まれていることに注意してください。

●それでは $J$ の最小値を与える $f$ を求めてみます。評価関数は

$(11)\quad \begin{array}{lll} \displaystyle{J=-\frac{q^2x^2(0)}{2(a-bf)}-\frac{r^2f^2x^2(0)}{2(a-bf)}}\\ \displaystyle{=\underbrace{-\frac{q^2+r^2f^2}{2(a-bf)}}_{\Pi}x^2(0)} \end{array} }$

のように書けますので、 $J$ の最小化は、次の $\Pi$ の最小化と等価です。

$\displaystyle{(12)\quad \boxed{\Pi=-\frac{q^2+r^2f^2}{2(a-bf)}} }$

図2 評価関数の概形

●この極値を求めるために $f$ で微分すると

$(13)\quad \begin{array}{lll} \displaystyle{\frac{d\Pi}{df}=-\frac{2r^2f}{2(a-bf)}-(-1)\frac{q^2+r^2f^2}{2(a-bf)^2}(-b)}\\ \displaystyle{=-\frac{2r^2(a-bf)f+(q^2+r^2f^2)b}{2(a-bf)^2}}\\ \displaystyle{=\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^2}} \end{array}$

これより

$\displaystyle{(14)\quad \frac{d\Pi}{df}=0 \ \Rightarrow \ bf^2-2af-\left(\frac{q}{r}\right)^2b=0 }$

ここで、 $f$ は $a-bf<0$ を満足させるものでなければならないので

$\displaystyle{(15)\quad \boxed{f=\frac{1}{b}\left(a+\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}\right)} }$

となります。そしてこれが一意に定まることは、以下のように $\Pi$ が下に凸でであることから分かります。

$(16)\quad \begin{array}{lll} \displaystyle{\frac{d^2\Pi}{df^2}=\frac{2br^2f-2ar^2}{2(a-bf)^2}+(-2)\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^3}(-b)}\\ \displaystyle{=\frac{-r^2}{(a-bf)}+b\frac{-2r^2(a-bf)f-br^2f^2-q^2b}{(a-bf)^3}}\\ \displaystyle{=\frac{-r^2(a-bf)^2-2r^2(a-bf)bf-r^2b^2f^2-q^2b^2}{(a-bf)^3}}\\ \displaystyle{=\frac{-r^2(a-bf+bf)^2-q^2b^2}{2(a-bf)^3}}\\ \displaystyle{=\frac{-r^2a^2-q^2b^2}{2(a-bf)^3}>0} \end{array}$

[2] LQ制御問題の解は求まったので、今度は評価関数(4)の重み係数 $q$ と $r$ の与え方について考えます。まず被積分項 $q^2x^2(t)+r^2u^2(t)$ は物理量の２乗和になっていますので、予め代表時間と代表長さによる無次元化を行っておく必要があります。たとえば、状態変数と入力変数がとりうる範囲を想定して

$\displaystyle{(17)\quad |x(t)|<M_x, |u(t)|<M_u\Rightarrow q=\frac{1}{M_x}, r=\frac{1}{M_u} }$

とすると、重み係数は物理単位の無次元化の役割を果たすことになります。次に(15)による閉ループ系を考えます。

$\displaystyle{(18)\quad \dot{x}(t)=(a-bf)x(t)=-\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}x(t) }$

ここで、 $a=\frac{1}{T},b=\frac{K}{T}$ の場合を考えると

$\displaystyle{(19)\quad \dot{x}(t)=-\frac{1}{T}\sqrt{1+\left(\frac{q}{r}\right)^2K^2}x(t) =-\frac{1}{\frac{T}{\sqrt{1+\left(\frac{q}{r}\right)^2K^2}}}x(t) }$

となって、新しい時定数は元の時定数 $T$ より小さくなることが分かります。これを決めるパラメータは重み係数の比 $q/r$ で、これが大きいほど時定数 $T$ は小さくなります。

●パラメータ $\rho$ を導入して重み係数を次のように設定します。

$\displaystyle{(20)\quad q=\frac{1}{M_x},\ r=\frac{1}{\rho M_u}\ \Rightarrow\ \frac{q}{r}=\rho\frac{M_u}{M_x}}$

ここで、 $\frac{M_u}{M_x}$ は固定されますので、設計パラメータは $\rho$ のみになります。

$a=-1$ 、 $b=1$ 、 $x(0)=1$ 、 $M_x=1$ 、 $M_u=1$ として、設計パラメータを $\rho=2,4,8$ と変えてみると、次の応答を得ます。

図3 重み係数の比 $\rho=q/r$ による応答の変化

演習A03-1…Flipped Classroom

$1^\circ$ 図2のグラフを描け。またパラメータ $\rho=q/r$ の値による状態と入力の振る舞いを観察し、トレードオフの関係について述べよ。

MATLAB

%a03_1.m
clear all, close all
a=-1; b=1; c=1;
Mx=1; Mu=1; q=1/Mx; r0=1/Mu
x0=1; 
t=0:0.1:4;
%-----
rho=2; r=r0/rho;
f1=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys1=ss(a-b*f1,b,[c;-f1],[]);
%-----
rho=4; r=r0/rho;
f2=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys2=ss(a-b*f2,b,[c;-f2],[]);
%-----
rho=8; r=r0/rho;
f3=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
sys3=ss(a-b*f3,b,[c;-f3],[]);
%-----
initial(sys1,sys2,sys3,x0,t)
grid
legend('rho=2','rho=4','rho=8'); 
title('LQ Control of 1st-order System')
xlabel('time')
ylabel('u(t) and x(t)')
%eof

SCILAB

coming soon

$2^\circ$ ハミルトン行列と呼ばれる

$\displaystyle{(21)\quad M=\left[\begin{array}{cc} a & -r^{-2}b^2 \\ -q^2 & -a \end{array}\right] }$

の安定な固有値

$\displaystyle{(22)\quad \lambda=-\sqrt{a^2+r^{-2}b^2q^2} }$

に対応する固有ベクトルは

$\displaystyle{(23)\quad \left[\begin{array}{cc} v_1 \\ v_2 \end{array}\right] =\left[\begin{array}{cc} 1 \\ \frac{-a-\sqrt{a^2+r^{-2}b^2q^2}}{-r^{-2}b^2} \end{array}\right] }$

となることを確かめ、次式を計算せよ。

$\displaystyle{(24)\quad f=r^{-2}bv_2v_1^{-1} }$

状態オブザーバ…Homework

[3] 状態フィードバックを使うには状態変数の変化をセンサで捉えてフィードバックする必要があるのですが、状態方程式が分かっているのですから、これを計算機内で解いてやれば良いのではないでしょうか？

いま制御対象（１次系）の状態空間表現を

$\displaystyle{(25)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t),\ x(0)\ne0\\ y(t)=cx(t) \end{array}\right. }$

とします。これを計算機で解きますが、その状態と出力を区別して

$\displaystyle{(26)\quad \left\{\begin{array}{lll} \dot{\hat{x}}(t)=a\hat{x}(t)+bu(t),\ \hat{x}(0)=0\\ \hat{y}(t)=c\hat{x}(t) \end{array}\right. }$

で表します。(26)から(25)の第１式どうしを辺々引くと

$\displaystyle{(27)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=a(\hat{x}(t)-x(t)),\ \hat{x}(0)\ne0 }$

を得ます。ここで $a<0$ でない限り

$\displaystyle{(28)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }$

とはなりません。そこで、図4に示すような工夫が行われました。

図4 状態オブザーバの考え方

● いま(2)に次のような補正項を加えます。

$\displaystyle{(29)\quad \dot{\hat{x}}(t)=a\hat{x}(t)+bu(t)-h(c\hat{x}(t)-y(t)),\ \hat{x}(0)=0 }$

すなわち

$\displaystyle{(30)\quad \dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+bu(t)+hy(t),\ \hat{x}(0)=0 }$

ここで、 $h$ は適当な定数です。(30)から(26)の第１式どうしを辺々引くと

$\displaystyle{(31)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=(a-hc)(\hat{x}(t)-x(t)),\ \hat{x}(0)\ne0 }$

を得ます。ここで $a-hc<0$ であるとすると

$\displaystyle{(32)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }$

が成り立ちます。すなわち状態オブザーバ(30)の状態 $\hat{x}(t)$ は、制御対象の状態 $x(t)$ を漸近的に推定することができます。このとき $a-hc<0$ を満たす(30)を状態オブザーバ、 $h$ をオブザーバゲインと呼びます。

● １次系

$\displaystyle{(33)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t),\ x(0)\ne0\\ y(t)=cx(t) \end{array}\right. }$

に対して、状態オブザーバ

$\displaystyle{(34)\quad \boxed{\dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+bu(t)+hy(t),\ \hat{x}(0)=0\quad (a-hc<0)} }$

によって推定された状態 $\hat{x}$ を用いる状態フィードバック

$\displaystyle{(35)\quad u(t)=-f\hat{x}(t) }$

による閉ループ系の状態方程式は、次式となります。

$\displaystyle{(36)\quad \left[\begin{array}{cc} \dot{x}(t)\\ \dot{\hat{x}}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a& -bf\\ hc& a-hc-bf \end{array}\right] }_{A_F} \left[\begin{array}{cc} x(t)\\ \hat{x}(t) \end{array}\right] }$

これから閉ループ系の安定性を判定したいのですが、このままでは見当がつきません。そこで、状態の推定誤差を $e(t)=\hat{x}(t)-x(t)$ とおくと、状態変数を変更した閉ループ系の状態方程式として、次式を得ます。

$\displaystyle{(37)\quad \boxed{ \left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a-bf & -bf\\ 0 & a-hc \end{array}\right] }_{A_F'} \left[\begin{array}{cc} x(t)\\ e(t) \end{array}\right]} }$

実際、状態の推定誤差の定義から、

$\displaystyle{(38)\quad \hat{x}(t)=x(t)+e(t) }$

に注意すると、状態フィードバック(35)は

$\displaystyle{(39)\quad u(t)=-fx(t)-fe(t) }$

これによる閉ループ系は、 $\dot{x}(t)=ax(t)+bu(t)$ に代入して

$\displaystyle{(40)\quad \dot{x}(t)=(a-bf)x(t)-bfe(t) }$

となります。これが(37)の第１式に他なりません。また、第２式は(31)そのもの、すなわち次式となります。

$\displaystyle{(41)\quad \dot{e}(t)=(a-hc)e(t) }$

閉ループ系の状態方程式は(40)と(41)に分離されており、 $a-bf$ と $a-hc$ は共に負となるように $f$ と $h$ を独立に設計してもよいことが分かります。この事実は分離定理と呼ばれています。

なお状態フィードバックを状態オブザーバの出力を用いて実施するコントローラ

$\displaystyle{(42)\quad \boxed{ \left\{\begin{array}{l} \dot{\hat{x}}(t)=(a-bf-hc)\hat{x}(t)+hy(t),\ \hat{x}(0)=0\\ u(t)=-f\hat{x}(t) \end{array}\right.} }$

は、オブザーバベースコントローラと呼ばれています。

演習A03-2…Flipped Classroom

$1^\circ$ １次系 $\dot{x}(t)=ax(t)+bu(t),y(t)=cx(t)$ に対する状態フィードバック $u(t)=-fx(t)$ と，状態オブザーバ $\dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+hy(t)+bu(t)$ （ $a-hc<0$ ）を考える。状態オブザーバによって推定された状態を用いる状態フィードバック $u(t)=-f\hat{x}(t)$ による閉ループ系の状態方程式を求めよ。

$2^\circ$ 閉ループ系の $A$ 行列の固有値が
$\lambda_1=-\frac{1}{T_1}$ ， $\lambda_2=-\frac{1}{T_2}$ （ $T_2>T_1>0$ ）となるように， $f$ と $h$ を定めよ。
$f$ と $h$ に， $\lambda_1$ と $\lambda_2$ をどう対応させたか，およびその理由を述べよ。

１次系のLQG制御…Homework

[5] １次系(33)に対して、オブザーバベースコントローラ(42)による閉ループ系に対して、評価関数

$\displaystyle{(43)\quad J=\int_0^\infty(q^2y\,^2(t)+r^2u\,^2(t))\,dt \label{eq7.2.7}}$

を最小化するように、状態フィードバックゲイン $f$ と状態オブザーバゲイン $h$ を決定するLQG制御問題を考えます。

●評価関数は次のように計算できます。

状態方程式

$\displaystyle{(44)\quad \underbrace{ \left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right] }_{\dot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{cc} a-bf & -bf\\ 0 & a-hc \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{cc} x(t) \\ e(t) \end{array}\right] }_{x_{CL}(t)} }$

と出力方程式

$\displaystyle{(45)\quad \underbrace{ \left[\begin{array}{cc} qx(t)\\ ru(t) \end{array}\right] }_{y_{CL}(t)} = \underbrace{ \left[\begin{array}{cc} q & 0 \\ -rf & -rf \end{array}\right] }_{C_{CL}} \underbrace{ \left[\begin{array}{cc} x(t) \\ e(t) \end{array}\right] }_{x_{CL}(t)} \label{eq7.2.5}}$

で表される閉ループ系の応答は

＜このパートは最初は読み飛ばして結構です＞
$\displaystyle{(46)\quad y_{CL}(t) = C_{CL}\exp(A_{CL}t) x_{CL}(0) \label{eq7.2.6}}$

このとき

$(47)\quad \begin{array}{lll} \displaystyle{J=\int_0^\infty y_{CL}^T(t)y_{CL}(t)\,dt}\\ \displaystyle{=\int_0^\infty x_{CL}^T(0)\exp(A_{CL}^Tt)C_{CL}^TC_{CL}\exp(A_{CL}t)x_{CL}(0)\,dt}\\ \displaystyle{=x_{CL}^T(0)\int_0^\infty \exp(A_{CL}^Tt)\underbrace{C_{CL}^TC_{CL}}_{Q_{CL}}\exp(A_{CL}t)\,dt\,x_{CL}(0)}\\ \displaystyle{=x_{CL}^T(0)\underbrace{\int_0^\infty \exp(A_{CL}^Tt)Q_{CL}\exp(A_{CL}t)\,dt}_{\Pi}\,x_{CL}(0)}\\ \displaystyle{=x_{CL}^T(0)\Pi x_{CL}(0)} \end{array}$

ここで、初期値

$\displaystyle{(48)\quad x_{CL}(0)= \left[\begin{array}{cc} x(0)\\ e(0) \end{array}\right]= \left[\begin{array}{cc} x(0)\\ \hat{x}(0)-x(0) \end{array}\right] \label{eq7.2.9}}$

をどう選ぶかですが、図4に示す閉ループ系を考えます。すなわち、(33)の代わりに、

$\displaystyle{(49)\quad \left\{\begin{array}{lll} \dot{x}(t)=ax(t)+bu(t)+\sigma n_x(t),\ x(0)\ne0\\ y(t)=cx(t)+\rho n_y(t) \end{array}\right. \label{eq7.2.10}}$

を考えます。ここで、新たに加えられた入力 $n_x(t)$ と $n_y(t)$ はそれぞれシステムノイズと観測ノイズとよばれ、定常正規確率過程として扱うのが普通です。特に、 $\sigma^2$ と $\rho^2$ は分散を表しています。ここでは確率的な議論を避けるため、 $n_x(t)$ と $n_y(t)$ は確定的なインパルス入力であると仮定します。

図5 オブザーバベース・コントローラによる閉ループ系の評価

このとき、閉ループ系の状態方程式は、次式となります。

$\displaystyle{(50)\quad &&\left[\begin{array}{cc} \dot{x}(t)\\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a-bf & -bf \\ 0 & a-hc \end{array}\right] }_{A_{CL}} \left[\begin{array}{cc} x(t)\\ e(t) \end{array}\right]\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} \sigma & 0\\ -\sigma & h\rho \end{array}\right] }_{B_{CL}} \left[\begin{array}{cc} n_x(t)\\ n_y(t) \end{array}\right] \label{eq7.2.11}}$

出力方程式(45)の場合のインパルス応答は、インパルス入力を $n_x$ に与える場合と $x_y$ に与える場合の２通り

＜このパートは最初は読み飛ばして結構です＞
$\displaystyle{(51)\quad y^{(1)}_{CL}(t) = C_{CL}\exp(A_{CL}t) \underbrace{ \left[\begin{array}{cc} \sigma \\ -\sigma \end{array}\right] }_{x_{CL}(0)=B^{(1)}_{CL}} \label{eq7.2.12}}$

$\displaystyle{(52)\quad y^{(2)}_{CL}(t) = C_{CL}\exp(A_{CL}t) \underbrace{ \left[\begin{array}{cc} 0 \\ h\rho \end{array}\right] }_{x_{CL}(0)=B^{(2)}_{CL}} \label{eq7.2.13}}$

があるので、評価関数を改めて次式のように両者の２乗和とします。

$\displaystyle{(53)\quad J'=B_{CL}^{(1)T}\Pi B_{CL}^{(1)}+B_{CL}^{(2)T}\Pi B_{CL}^{(2)} \label{eq7.2.14}}$

これは、 $\Pi= \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]$ とおくと

$\displaystyle{(54)\quad J'=\pi_1\sigma^2+\pi_2(\sigma^2+h^2\rho^2)-2\pi_3\sigma^2 \label{eq7.2.15}}$

と計算できます。ここで、 $\pi_1,\pi_2,\pi_3$ を直接計算することは容易ではないのですが、

＜このパートは最初は読み飛ばして結構です＞
$A_{CL}$ が安定行列であることは， $\pi_1,\pi_2,\pi_3$ が，リャプノフ方程式

$\displaystyle{(55)\quad \begin{array}{lll} && \underbrace{ \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] }_{\Pi} \underbrace{ \left[\begin{array}{cc} a-bf & -bf \\ 0 & a-hc \end{array}\right] }_{A_{CL}}\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} a-bf & 0 \\ -bf & a-hc \end{array}\right] }_{A_{CL}^T} \underbrace{ \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] }_{\Pi}\nonumber\\ &&+ \underbrace{ \left[\begin{array}{cc} q^2+r^2f^2 & r^2f^2 \\ r^2f^2 & r^2f^2 \end{array}\right] }_{C_{CL}^TC_{CL}>0} = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array} \label{eq7.2.16}}$

すなわち

$\displaystyle{(56)\quad 2(a-bf)\pi_1+q^2+r^2f^2=0 \label{eq7.2.17}}$
$\displaystyle{(57)\quad -bf\pi_1+(2a-bf-hc)\pi_3+r^2f^2=0 \label{eq7.2.18}}$
$\displaystyle{(58)\quad 2(a-hc)\pi_2-2bf\pi_3+r^2f^2=0 \label{eq7.2.19}}$

の解であることに等しいという結果を用います。

(56)より $\pi_1$ を求めれば，評価関数は

$\displaystyle{(59)\quad J'= \underbrace{ \frac{q^2+r^2f^2}{-2(a-bf)} }_{\pi_1} \sigma^2 +\pi_2(\sigma^2+h^2\rho^2)-2\pi_3\sigma^2 \label{eq7.2.20}}$

となります。ここで，第1項は，(11)で， $x(0)=\sigma$ としたものであることに注意します。いま， $\hat{x}(0)=x(0)$ かつ $w(t)=v(t)=0$ という場合を考えてみると，(50)の第2行目より

$\displaystyle{(60)\quad \dot{e}(t)=(a-hc)e(t),e(0)=0\ \Rightarrow\ e(t)=0\ \Rightarrow\ \hat{x}(t)=x(t) \label{eq7.2.21}}$

だから，(35)は状態フィードバック $u(t)=-fx(t)$ に相当します。この場合の議論はすでに行っており， $\pi_2,\pi_3$ は存在しないので，評価関数(59)は第1項だけとなり，これを最小化する $f$ は，リッカチ方程式

$\displaystyle{(61)\quad -r^{-2}b^2\pi_1^2+2a\pi_1+q^2=0 \label{eq7.2.22}}$

の解 $\pi_1>0$ を用いて次のように与えられます。

$\displaystyle{(62)\quad f=r^{-2}b\pi_1 \label{eq7.2.23}}$

そうすると，評価関数(59)の第2項と第3項は，状態オブザーバを用いるために生じると考えられます。このように，理想的な状態フィードバックに比べて評価関数の値が増えることを評価関数の劣化と呼びます。

つぎに，仮定 $f=r^{-2}b\pi_1$ ， $a-bf<0$ ， $a-hc<0$ のもとでは，(57)より $\pi_3=0$ を得ます。このとき，式(58)は

$\displaystyle{(63)\quad 2(a-hc)\pi_2+r^2f^2=0 \label{eq7.2.24}}$

となり，これから $\pi_2$ が求まります。したがって，評価関数(59)は

$\displaystyle{(64)\quad J'= \underbrace{ \frac{q^2+r^2f^2}{-2(a-bf)} }_{\#1} \sigma^2 + \underbrace{ \frac{\sigma^2+h^2\rho^2}{-2(a-hc)} }_{\#2} r^2f^2 \label{eq7.2.25}}$

と書けます。評価関数の劣化分を最小化するためには，#2を $h$ について最小化すればよいことになります。幸いなことに，#2は#1と同形であるので，これを最小化する $h$ の決定法は $f$ と同様に行ないます。

●以上の議論をまとめると、次のようになります。

アルゴリズム＜1次系のLQG制御＞

入力パラメータ
$a,\,b,\,c,\,q>0,\,r>0,\,\sigma>0,\,\rho>0$
出力パラメータ
$a_K,\,b_K,\,c_K$

ステップ1: 状態フィードバックゲイン $f$ の決定

$\displaystyle{(65)\quad \boxed{{\rm\bf CARE:\ }\displaystyle{-\frac{1}{r^2}b^2\Pi^2+2a\Pi+q^2=0}} }$

を解いて， $\Pi>0$ を求め， $f$ をつぎのように定める。

$\displaystyle{(66)\quad \boxed{f=\frac{1}{r^2}b\Pi} }$

ステップ2: オブザーバゲイン $h$ の決定

$\displaystyle{(67)\quad \boxed{{\rm\bf FARE:\ }\displaystyle{-\frac{1}{\rho^2}c^2\Gamma^2+2a\Gamma+\sigma^2=0}} }$

を解いて， $\Gamma>0$ を求め， $h$ をつぎのように定める。

$\displaystyle{(68)\quad \boxed{h=\frac{1}{\rho^2}c\Gamma} }$

ステップ3: オブザーバベースコントローラの構成

$\displaystyle{(69)\quad \left\{\begin{array}{l} \dot{x}_K(t)=a_Kx_K(t)+b_Ky(t)\\ u(t)=c_Kx_k(t) \end{array}\right. }$

ただし

$\displaystyle{(70)\quad a_K=a-hc-bf,\ b_K=h,\ c_K=-f }$

● $a=-1$ 、 $b=1$ 、 $c=1$ 、 $q=1$ 、 $r=1$ 、 $\sigma=1$ 、 $\rho=0.2$ 、 $x(0)=1$ 、 $\hat{x}(0)=0$ として、状態フィードバックだけとオブザーバベースコントローラを用いた場合とを比較してみると、次の応答を得ます。

図6 状態フィードバックとオブザーバベースコントローラを用いた場合の応答の変化

演習A03-3…Flipped Classroom

$1^\circ$ 図6のグラフを描け。

MATLAB

%a03_2.m
 clear all, close all
 a=-1; b=1; c=1; 
 q=1; r=1; 
 sig=1; rho=0.2;
 f=1/b*(a+sqrt(a^2+(q/r)^2*b^2));
 sys1=ss(a-b*f,[],[c;-f],[]);
 x0=1; 
 t=0:0.1:4;
 initial(sys1,x0,t), hold on
 h=1/c*(a+sqrt(a^2+(sig/rho)^2*c^2));
 AF=[a -b*f;
     h*c a-h*c-b*f];
 sys2=ss(AF,[],diag([c,-f]),[]);
 initial(sys2,[x0;0],t)  
 grid
 legend('sfc','obc')
 title('LQ OBC of 1st-order System')
 xlabel('time')
 ylabel('u(t) and x(t)')
%eof

SCILAB

coming soon

$2^\circ$ 上のシミュレーションにおいて、 $a-hc<a-bc$ となっていることを確かめよ。

Note A03　ラグランジュの未定定数法

１入力１出力１次系に対する最適制御（LQ制御）を考えましたが、一般の多入力多出力高次系の場合も同様な議論ができて、大変有用な制御手段として知られています。その説明はあとで行ないますが、そこではラグランジュの未定定数法によるアプローチが採用されます。

(11)の評価関数 $J$ を最小化する問題を、 $\Gamma$ をラグランジュの未定定数として

$\displaystyle{(1)\quad J=\Pi x^2(0)+\Gamma(2(a-bf)\Pi+q^2+r^2f^2) }$

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

$\displaystyle{(2)\quad \boxed{2(a-bf)\Pi+q^2+r^2f^2=0} }$

ですが、 $\Pi>0$ は $a-bf<0$ を意味することに注意します。そこで、必要条件として次を得ます。

$\displaystyle{(3)\quad \frac{\partial J}{\partial \Pi}&=&x^2(0)+2(a-bf)\Gamma=0\Rightarrow\Gamma>0 }$
$\displaystyle{(4)\quad \frac{\partial J}{\partial f}&=&(-2b\Pi+2r^2f)\Gamma=0\Rightarrow f=r^{-2}b\Pi }$
$\displaystyle{(5)\quad \frac{\partial J}{\partial \Gamma}&=&2(a-bf)\Pi+q^2+r^2f^2=0 }$