icare

連続時間の代数リカッチ方程式に対する陰的ソルバー

構文

[X,K,L] = icare(A,B,Q,R,S,E,G)

[X,K,L,info] = icare(___)

[___] = icare(___,'noscaling')

[___] = icare(___,'anti')

説明

[X,K,L] = icare(A,B,Q,R,S,E,G) は、次の連続時間の代数リカッチ方程式について、一意の安定化解 X、状態フィードバック ゲイン K、および閉ループ固有値 L を計算します。

$$A^{T}XE+E^{T}XA+E^{T}XGXE-\left(E^{T}XB+S\right)R^{-1}\left(B^{T}XE+S^T\right)+Q=0$$

安定化解 X は、すべての固有値 L を左半平面に配置します。

代数リカッチ方程式は、LQR/LQG 制御、H2∞ および H∞ 制御、カルマン フィルター処理、およびスペクトル分解や既約分解において重要な役割を果たします。

[X,K,L,info] = icare(___) は、連続時間の代数リカッチ方程式の解に関する追加の情報を含む構造体 info も返します。

[___] = icare(___,'noscaling') は、組み込みのスケーリングをオフにしてスケーリング ベクトル info.Sx および info.Sr のすべての要素を 1 に設定します。スケーリングをオフにすると計算は高速化されますが、A,B,Q,R,S,E,G のスケーリングが不適切だと精度に悪影響が及ぶことがあります。

[___] = icare(___,'anti') は、すべての固有値 L を右半平面に配置する、反安定化解 X を計算します。

例

連続時間の代数リカッチ方程式を解く

代数リカッチ方程式 ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-\mathrm{XB}{\mathit{B}}^{\mathit{T}}\mathit{X}+\mathit{C}{\mathit{C}}^{\mathit{T}}=0$ を解くには、次の行列について考えます。

$$A=\left[\begin{array}{ccc}1& -2& 3\\ -4& 5& 6\\ 7& 8& 9\end{array}\right]B=\left[\begin{array}{c}5\\ 6\\ -7\end{array}\right]C=\left[\begin{array}{ccc}7& -8& 9\end{array}\right].$$

一番簡単な方法は、$\mathit{G}=-\mathit{B}{\mathit{B}}^{\mathit{T}}$ と $\mathit{Q}={\mathit{C}}^{\mathit{T}}\mathit{C}$ を使用した後、icare を使用して解を求めることです。

A = [-1,2,3;4,5,-6;7,-8,9];
B = [5;6;-7];
C = [7,-8,9];
G = -B*B';
Q = C'*C;
[X1,K1,L1] = icare(A,[],Q,[],[],[],G)

X1 = 3×3

   15.3201    4.2369   17.0090
    4.2369    2.6252    4.4123
   17.0090    4.4123   19.0374

K1 =

  0x3 empty double matrix

L1 = 3×1

   -3.2139
  -10.1191
  -76.9693

上記の方法では、行列 B と C に大きな要素がある場合、行列 G と Q を計算するためにこれらが正方にされるため、数値が不正確になる可能性があります。数値範囲が限られているため、計算の精度が落ちたり、失敗することもあります。たとえば、norm(B) が 1e6 の場合、norm(G) は 1e12 となり、虚軸の 1e-4 以内のすべての固有値が数値誤差のために '虚数部' として診断される可能性があります。

数値精度を上げるには、代数リカッチ方程式を次のように書き換えます。

${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}–[\mathrm{XB},{\mathit{C}}^{\mathit{T}}]*[I,0;0,-I][{\mathit{B}}^{\mathit{T}}\mathit{X};C]=0$.

上記の方程式は、${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-(\mathrm{XB}+S){\mathit{R}}^{-1}(B^{\mathit{T}}\mathit{X}+S^{\mathit{T}})=0$ の標準形式であり、

ここで、$\mathit{B}=\left[\mathit{B},\mathrm{0}\right],\mathit{S}=\left[0,{\mathit{C}}^{\mathit{T}}\right],\mathrm{and}\mathit{R}=\left[\mathit{I},0;0,-\mathit{I}\right]$ です。

上の値で icare を使用して、解を計算します。

n = size(A,1);
m = size(B,2);
p = size(C,1);
R = blkdiag(eye(m),-eye(p));
BB = [B,zeros(n,p)];
S = [zeros(n,m),C'];
[X2,K2,L2,info] = icare(A,BB,0,R,S,[],[])

X2 = 3×3

   15.3201    4.2369   17.0090
    4.2369    2.6252    4.4123
   17.0090    4.4123   19.0374

K2 = 2×3

  -17.0406    6.0501  -21.7435
   -7.0000    8.0000   -9.0000

L2 = 3×1

   -3.2139
  -10.1191
  -76.9693

info = struct with fields:
        Sx: [3x1 double]
        Sr: [2x1 double]
         U: [3x3 double]
         V: [3x3 double]
         W: [2x3 double]
    Report: 0

ここで、X2 は一意の安定化解であり、K2 には状態フィードバック ゲインが含まれ、L2 には閉ループ固有値が含まれています。

連続時間の代数リカッチ方程式の反安定化解

連続時間の代数リカッチ方程式 ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-\mathrm{XB}{\mathit{B}}^{\mathit{T}}\mathit{X}+\mathit{C}{\mathit{C}}^{\mathit{T}}=0$ の反安定化解を求めるには、次の行列について考えます。

$$A=\left[\begin{array}{ccc}1& -2& 3\\ -4& 5& 6\\ 7& 8& 9\end{array}\right]B=\left[\begin{array}{c}5\\ 6\\ -7\end{array}\right]C=\left[\begin{array}{ccc}7& -8& 9\end{array}\right].$$

数値精度を上げるには、代数リカッチ方程式を次のように書き換えます。

${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}–[\mathrm{XB},{\mathit{C}}^{\mathit{T}}]*[I,0;0,-I][{\mathit{B}}^{\mathit{T}}\mathit{X};C]=0$.

上記の方程式は、${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-(\mathrm{XB}+S){\mathit{R}}^{-1}(B^{\mathit{T}}\mathit{X}+S^{\mathit{T}})=0$ の標準形式であり、

ここで、$\mathit{B}=\left[\mathit{B},\mathrm{0}\right],\mathit{S}=\left[0,{\mathit{C}}^{\mathit{T}}\right],\mathrm{and}\mathit{R}=\left[\mathit{I},0;0,-\mathit{I}\right]$ です。

'anti' オプションを使用して反安定化解を計算します。

A = [-1,2,3;4,5,-6;7,-8,9];
B = [5;6;-7];
C = [7,-8,9];
n = size(A,1);
m = size(B,2);
p = size(C,1);
R = blkdiag(eye(m),-eye(p));
BB = [B,zeros(n,p)];
S = [zeros(n,m),C'];
[X,K,L] = icare(A,BB,0,R,S,[],[],'anti')

X = 3×3

  -18.0978   10.9090   -1.8466
   10.9090   -6.7195    1.4354
   -1.8466    1.4354   -0.9426

K = 2×3

  -12.1085    4.1803    5.9774
   -7.0000    8.0000   -9.0000

L = 3×1

   76.9693
   10.1191
    3.2139

ここで、X は一意の反安定化解であり、K には状態フィードバック ゲインが含まれ、L には閉ループ固有値を含まれています。

入力引数

A,B,Q,R,S,E,G — 入力行列

行列

入力行列。行列として指定します。

行列 Q、R、および G はエルミート行列でなければなりません。正方行列がその複素共役転置と等しい場合、つまり $$a_{i,j}={\overline{a}}_{j,i}$$ である場合、これはエルミート行列です。

エルミート行列の詳細については、ishermitian を参照してください。

行列 R と E は可逆でなければなりません。

行列 B、R、S、E、G が省略されたり、[] に設定された場合、icare は次の既定値を使用します。

• B = 0

• R = I

• S = 0

• E = I

• G = 0

入力 Q、R、G がスカラー値の場合、icare はこれらを単位行列の倍数として解釈します。

'noscaling' — 組み込みのスケーリングをオフにするオプション

'noscaling'

組み込みのスケーリングをオフにするオプション。'noscaling' として指定します。組み込みのスケーリングをオフにした場合、icare はスケーリングベクトル info.Sx および info.Sr のすべての要素を 1 に設定します。スケーリングをオフにすると計算は高速化されますが、A,B,Q,R,S,E,G のスケーリングが不適切だと精度に悪影響が及ぶことがあります。

'anti' — 反安定化解を計算するオプション

'anti'

反安定化解を計算するオプション。'anti' として指定します。このオプションを有効にすると、icare は、(A+G*X*E-B*K,E) のすべての固有値を右半平面に配置する反安定化解 X を計算します。

リカッチ微分方程式の完全な位相図を把握するには、一意の安定化と反安定化の両方が必要です。

出力引数

X — リカッチ方程式の一意の解

行列

連続時間の代数リカッチ方程式の一意の解。行列として返されます。

既定では、X は連続時間の代数リカッチ方程式の安定化解です。'anti' オプションが使用される場合、X は反安定化解になります。

関連するハミルトニアン行列が虚軸上に固有値をもつ場合、icare は X に [] を返します。

ペンシルが特異である場合、つまり [B;S;R] がランク落ちの場合、icare は [] を返します。 (R2024a 以降)

K — 状態フィードバック ゲイン

行列

状態フィードバック ゲイン。行列として返されます。

状態フィードバック ゲイン K は次のように計算されます。

$$K=R^{-1}\left(B^{T}XE+S^T\right).$$

関連するハミルトニアン行列が虚軸上に固有値をもつ場合、icare は K に [] を返します。

ペンシルが特異である場合、つまり [B;S;R] がランク落ちの場合、icare は [] を返します。 (R2024a 以降)

L — 閉ループ固有値

行列

閉ループ固有値。行列として返されます。

閉ループ固有値 L は次のように計算されます。

$$L=eig\left(A+GXE-BK,E\right).$$

関連するハミルトニアン行列が虚軸上に固有値をもつ場合、icare は X と K に [] を返します。言い換えれば、X と K が空の行列であっても、L は空ではありません。

ペンシルが特異である場合、つまり [B;S;R] がランク落ちの場合、icare は L = NaN(n,1) を返します。 (R2024a 以降)

info — 一意の解に関する情報

構造体

一意の解に関する情報。次のフィールドをもつ構造体として返されます。

• Sx — 状態のスケーリングに使用される値のベクトル。

• Sr — 行列 R のスケーリングに使用される値のベクトル。

• U、V、および W — 関連するスケーリングされた行列ペンシルの安定した不変部分空間の基底を表す値のベクトル。詳細については、アルゴリズムを参照してください。

• Report — 以下のいずれかの値をもつスカラー。

  • 0 — 一意の解は正確です。

  • 1 — 解があまり正確でありません。

  • 2 — 解が有限でありません。

  • 3 — [L;-L] で表されるハミルトニアン スペクトルが虚軸上に固有値をもつため、解が見つかりません。

  • 4 — ペンシルが特異 ([B;S;R] がランク落ち) です。 (R2024a 以降)

制限

• 有限安定化解 X が存在し有限であるためには、(A-sE,B) が可安定で、E と R は可逆でなければなりません。一般的にはこれらの条件では不十分ですが、以下の条件が満たされる場合は十分となります。

  • $$\left[\begin{array}{cc}\begin{array}{c}Q\\ S^T\end{array}& \begin{array}{c}S\\ R\end{array}\end{array}\right]\ge 0$$

  • $$\left[\begin{array}{cc}A-B{R}^{-1}{S}^T& Q-S{R}^{-1}{S}^T\end{array}\right]$$ が検出可能

アルゴリズム

不変部分空間の基底

icare は次のペンシルで機能し、このペンシルの安定または反安定の有限固有値に関連付けられた不変部分空間の基底 [U;V;W] を計算します。

$$M-sN=\left[\begin{array}{ccc}A& G& B\\ -Q& -A^T& -S\\ S^T& B^T& R\end{array}\right]-s\left[\begin{array}{ccc}E& 0& 0\\ 0& E^T& 0\\ 0& 0& 0\end{array}\right]$$

データは、虚軸の近傍にある固有値の感度を減らし、安定な不変部分空間と反安定の不変部分空間の間隔を大きくするために、自動的にスケーリングされます。

解、状態フィードバック ゲイン、およびスケーリング ベクトル間の関係

解 X と状態フィードバック ゲイン K は、次の一連の方程式によってスケーリング ベクトルと U,V,W に関係付けられています。

$$\begin{array}{l}X={\text{ D}}_{x}V{U}^{-1}{\text{ D}}_x{E}^{-1}\text{,}\\ K=-D_r{\text{ WU}}^{-1}{\text{ D}}_x,\end{array}$$

ここで、

$$\begin{array}{l}{D}_x=\mathrm{diag}\left(S_x\right),\\ D_r=\mathrm{diag}\left(S_r\right).\end{array}$$

バージョン履歴

R2019a で導入

R2024a: 特異点の検出の改善

ペンシルが特異である場合、つまり [B;S;R] がランク落ちの場合に、関数で X = []、K = []、L = NaN(n,1)、info.Report = 4 が返されるようになりました。

参考

ishermitian | idare | lyap | lqr | lqg | kalman | h2syn (Robust Control Toolbox) | hinfsyn (Robust Control Toolbox) | spectralfact | lncf (Robust Control Toolbox) | rncf (Robust Control Toolbox)