balreal

平衡化状態空間実現

構文

[sysb,g] = balreal(sys)

[sysb,g] = balreal(sys,Name=Value)

[sysb,g,TL,TR] = balreal(sys,___)

説明

[sysb,g] = balreal(sys) は、LTI モデル sys の平衡化状態空間実現を計算します。安定なモデルの場合、sys は、可制御グラミアンと可観測グラミアンが等しく、対角行列であり、対角要素がハンケル特異値のベクトル g で構成される等価実現です。これにより、入力から状態および状態から出力へのエネルギー伝達が平衡化されます。g の小さいエントリは、xelim で削除できる状態を示します。

sys が不安定な極をもつ場合、関数は安定部分を分離して平衡化してから、不安定部分に戻して sysb を形成します。不安定モードに対応する g の要素は、Inf に設定されます。

[sysb,g] = balreal(sys,Name=Value) は、1 つ以上の名前と値の引数として指定したオプションを使用して平衡化実現を計算します。 (R2023b 以降)

[sysb,g,TL,TR] = balreal(sys,___) は、平衡化状態変換行列 TL および TR も返します。 (R2023b 以降)

例

安定したシステムの平衡実現

次のほぼ相殺された極-零点のペアをもつ零点-極-ゲイン モデルを考えます。

sys = zpk([-10 -20.01],[-5 -9.9 -20.1],1)

sys =
 
     (s+10) (s+20.01)
  ----------------------
  (s+5) (s+9.9) (s+20.1)
 
Continuous-time zero/pole/gain model.

平衡化されたグラミアンをもつ状態空間実現を取得します。

[sysb,g] = balreal(sys);

結合グラミアンの対角要素を確認します。

g'

ans = 1×3

    0.1006    0.0001    0.0000

これは、sysb の最後の 2 つの状態と、入力および出力との結合が弱いことを示しています。次に、xelim を使用して、これらの状態を削除できます。

sysr = xelim(sysb,[2 3]);

これにより、元のシステムの 1 次近似を取得することができます。

zpk(sysr)

ans =
 
  -0.00011597 (s-8638)
  --------------------
       (s+4.981)
 
Continuous-time zero/pole/gain model.

元のモデルと低次元化されたモデルのボード応答を比較します。

opt = bodeoptions;
opt.PhaseMatching = "on";
bodeplot(sys,sysr,'r--',opt)

このプロットは、2 番目と 3 番目の状態を削除してもシステム ダイナミクスにあまり影響がないことを示します。

高速で正確な結果を得るには、モデルの低次元化ワークフローに reducespec を使用します。

不安定なシステムの平衡実現

この例では、不安定な極をもつシステムの平衡化実現を実行する方法を示します。

以下のシステムを考えます。

$$\mathit{H}\left(\mathit{s}\right)=\frac{0.1}{\left(\mathit{s}+2\right)\left(\mathit{s}+1\right)\left(\mathit{s}-1\right)\left(\mathit{s}+0.001\right)}$$

システムには 3 つの安定した極と 1 つの不安定な極が含まれています。1 つの安定した極は安定境界のごく近くにあります。

モデルを作成します。

sys = ss(zpk(1,[-2 -1 1 -0.001],0.1));

balreal を適用して、平衡化実現を作成します。

[sysbal,g,TL,TR] = balreal(sys);
g

g = 4×1

       Inf
   25.0374
    0.0391
    0.0018

sys が不安定な極をもつ場合、balreal は安定部分を分離して平衡化してから、不安定部分を戻して平衡化実現を形成します。不安定な極が、ベクトル g に Inf として現れます。

オフセット オプションを使用して安定境界をシフトし、名目上は安定している極を除外することもできます。

[sysbal,g,TL,TR] = balreal(sys,Offset=1e-2);
g

g = 4×1

       Inf
       Inf
    0.0393
    0.0018

この実現では、balreal は s = –0.001 の極を不安定として扱います。

入力引数

sys — 動的システム

数値 LTI モデル

分解する動的システム。ss または tf モデルなどの数値 LTI モデルとして指定します。

名前と値の引数

引数のオプションのペアを Name1=Value1,...,NameN=ValueN として指定します。Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後に表示されなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: [sys,g2] = balreal(sys,Goal="relative",FreqIntervals=[10,100]);

Goal — 平衡化の目標

"absolute" (既定値) | "relative"

平衡化の目標。"absolute" または "relative" として指定します。このオプションを使用して平衡化アルゴリズムを選択し、低次元化されたモデル G_r のタイプを制御します。

• "absolute" — 絶対誤差 $${\Vert G-G_r\Vert }_{\infty }$$ を最小化します。

• "relative" — 相対誤差 $${\Vert G^{-1}\left(G-G_r\right)\Vert }_{\infty }$$ を最小化します。

相対誤差が周波数全体のより良好な一致を表すのに対し、絶対誤差は大部分のゲインをもつ領域を強調します。

Regularization — 正則化レベル

"auto" (既定値) | 非負のスカラー

正則化レベル。"auto" または非負のスカラー値として指定します。Goal を "relative" に設定すると、balreal は sys ではなくモデル [sys,α*I] を平衡化します。このオプションを "auto" に設定すると、balreal は適切な正則化レベルの α 値を選択します。値 α ≥ 0 を指定すると、この既定の設定をオーバーライドします。

この正則化値は、すべての周波数において適切に定義された相対誤差を提供するのに役立ちます。

FreqIntervals — グラミアンの周波数範囲の制限

[] (既定値) | 2 列の行列

周波数を制限したハンケル特異値を計算するための周波数範囲。2 つの列をもつ行列として指定します。各行は周波数範囲 [fmin,fmax] を指定します。ここで、fmin および fmax は非負の周波数であり、モデルの周波数単位で表されます。切り捨てる低エネルギー状態を特定する場合、ソフトウェアはこれらの周波数範囲のみでのシステム動作に対する状態の寄与を計算します。以下に例を示します。

• モデルの周波数単位を rad/s と仮定して、計算の範囲を 3 ～ 15 rad/s に制限するには、FreqIntervals を [3,15] に設定します。

• 計算を 2 つの周波数範囲 (3 ～ 15 rad/s と 40 ～ 60 rad/s) に制限するには、[3,15;40,60] を使用します。

• カットオフ周波数 fcut より低いすべての周波数を指定するには、[0,fcut] を使用します。

• カットオフを超えるすべての周波数を指定するには、連続時間では [fcut,Inf]、離散時間では [fcut,pi/Ts] を使用します。ここで、Ts はモデルのサンプル時間を表します。

既定値 [] では周波数制限が適用されません。この既定値は、連続時間の [0,Inf] または離散時間の [0,pi/Ts] と等価です。ただし、[] 以外の TimeIntervals 値を指定する場合、この制限によって FreqIntervals = [] がオーバーライドされます。TimeIntervals 値と FreqIntervals 値の両方を指定する場合、計算ではこれらの間隔の和集合が使用されます。

TimeIntervals — グラミアンの時間間隔の制限

[] (既定値) | 2 列の行列

時間を制限したハンケル特異値を計算するための時間間隔。2 つの列をもつ行列として指定します。各行は時間間隔 [tmin,tmax] を指定します。ここで、tmin および tmax は非負の時間であり、モデルの時間単位で表されます。切り捨てる低エネルギー状態を特定する場合、ソフトウェアはこれらの時間間隔にあるシステムのインパルス応答に対する状態の寄与のみを計算します。以下に例を示します。

• モデルの時間単位を秒と仮定して、計算の範囲を 3 ～ 15 秒に制限するには、TimeIntervals を [3,15] に設定します。

• 計算を 2 つの時間間隔 (3 ～ 15 秒と 40 ～ 60 秒) に制限するには、[3,15; 40,60] を使用します。

• ゼロから最大でカットオフ時間 tcut までのすべての時間を指定するには、[0,tcut] を使用します。カットオフ後のすべての時間を指定するには、[tcut,Inf] を使用します。

既定値 [] では時間制限が適用されません。この既定値は [0,Inf] と等価です。ただし、[] 以外の FreqIntervals 値を指定する場合、この制限によって Timeintervals = [] がオーバーライドされます。TimeIntervals 値と FreqIntervals 値の両方を指定する場合、計算ではこれらの間隔の和集合が使用されます。

SepTol — 安定/不安定分解での精度の最大損失

10 (既定値) | 正のスカラー

安定/不安定分解の最大精度損失係数。正のスカラーとして指定します。不安定な極をもつモデルの場合、balreal はまず変換 TL、TR を使用して安定ダイナミクスを抽出します。SepTol を使用して分解の精度を制御します。SepTol を増やすことは、精度を犠牲にして、近傍の安定モードと不安定モードを分離する場合には役立ちます。詳細については、stabsep を参照してください。

Offset — 安定/不安定境界のオフセット

1e-8 (既定値) | 正のスカラー

安定/不安定境界のオフセット。正のスカラー値として指定します。安定/不安定分解では、安定項は次の条件を満たす極のみを含みます。

• Re(s) < -Offset × max(1,|Im(s)|) (連続時間)

• |z| < 1 - Offset (離散時間)

安定境界の近くの極を不安定として扱うには、Offset の値を増加します。

出力引数

sysb — 平衡化状態空間実現

状態空間モデル

平衡化状態空間実現。状態空間モデルとして返されます。

g — ハンケル特異値

ベクトル

ハンケル特異値。ベクトルとして返されます。

不安定モードに対応する g の要素は、Inf に設定されます。

TL — 変換の左辺の行列

行列

R2023b 以降

平衡化状態変換の左辺の行列。Nx 行 Nx 列の次元をもつ行列として返されます。ここで Nx はモデル sys 内の状態の数です。

アルゴリズムはモデルの状態空間実現 (A, B, C, D, E) を (A_b, B_b, C_b, D_b, E_b) に変換します。これは以下によって与えられます。

• 陽的な状態空間モデルの場合

  $$\begin{array}{l}\begin{array}{cc}{A}_b=T_{L}A{T}_R,& B_b=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_b=C{T}_R,& D_b=D,& E_b=T_L{T}_R=I\end{array}\end{array}$$

• 記述子状態空間モデルの場合

  $$\begin{array}{l}\begin{array}{cc}{A}_b=T_{L}A{T}_R,& B_b=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_b=C{T}_R,& D_b=D,& E_b=T_{L}E{T}_R\end{array}\end{array}$$

入力モデル sys が状態空間モデルでない場合、関数はこの引数に対して空の値 [] を返します。

TR — 変換の右辺の行列

行列

R2023b 以降

平衡化状態変換の右側の行列。Nx 行 Nx 列の次元をもつ行列として返されます。ここで、Nx はモデル sys 内の状態の数です。アルゴリズムはモデルの状態空間実現 (A, B, C, D, E) を (A_b, B_b, C_b, D_b, E_b) に変換します。これは以下によって与えられます。

• 陽的な状態空間モデルの場合

  $$\begin{array}{l}\begin{array}{cc}{A}_b=T_{L}A{T}_R,& B_b=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_b=C{T}_R,& D_b=D,& E_b=T_L{T}_R=I\end{array}\end{array}$$

• 記述子状態空間モデルの場合

  $$\begin{array}{l}\begin{array}{cc}{A}_b=T_{L}A{T}_R,& B_b=T_{L}B,\end{array}\\ \begin{array}{ccc}{C}_b=C{T}_R,& D_b=D,& E_b=T_{L}E{T}_R\end{array}\end{array}$$

入力モデル sys が状態空間モデルでない場合、関数はこの引数に対して空の値 [] を返します。

ヒント

モデル次数の低次元化のために、reducespec を使用します。

アルゴリズム

安定した状態空間モデルを考えます。

$$\begin{array}{l}E\dot{x}=Ax+Bu\\ y=Cx+Du\end{array}$$

これには、可制御グラミアン W_c と可観測グラミアン W_o があります。

状態座標変換 $$x=T_R\overline{x}$$ は、次の等価モデルを生成します。

$$\begin{array}{l}{T}_{L}E{T}_R\dot{\overline{x}}=T_{L}A{T}_R\overline{x}+T_{L}Bu\\ y=C{T}_R\overline{x}+Du\end{array}$$

balreal はグラミアンを次のものに変換します。

$$\begin{array}{cc}{\overline{W}}_c=T_L{W}_c{T}_L^T,& {\overline{W}}_o=T_R^T{W}_o\end{array}{T}_R$$

その際、以下が満たされるようにします。

$${\overline{W}}_c={\overline{W}}_o=diag(g)$$

モデルが不安定な極をもつ場合、関数は安定部分を分離して平衡化してから、不安定部分に戻して実現を形成します。不安定モードに対応する g の要素は、Inf に設定されます。

アルゴリズムの詳細については、[1] および [2] を参照してください。

TimeIntervals オプションまたは FreqIntervals オプションを使用した場合、balreal は時間または周波数を制限した可制御グラミアンと可観測グラミアンを基に平衡実現を行います。時間および周波数を制限したグラミアンの計算の詳細については、gram および [4] を参照してください。

参照

[1] Laub, A., M. Heath, C. Paige, and R. Ward. “Computation of System Balancing Transformations and Other Applications of Simultaneous Diagonalization Algorithms.” IEEE Transactions on Automatic Control 32, no. 2 (February 1987): 115–22. https://doi.org/10.1109/TAC.1987.1104549.

[2] Moore, B. “Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction.” IEEE Transactions on Automatic Control 26, no. 1 (February 1981): 17–32. https://doi.org/10.1109/TAC.1981.1102568.

[3] Laub, Alan J. “Computation of ‘Balancing’ Transformations.” Joint Automatic Control Conference, no. 17 (1980): 84. https://doi.org/10.1109/JACC.1980.4232165.

[4] Gawronski, Wodek, and Jer-Nan Juang. “Model Reduction in Limited Time and Frequency Intervals.” International Journal of Systems Science 21, no. 2 (February 1990): 349–76. https://doi.org/10.1080/00207729008910366.

バージョン履歴

R2006a より前に導入

R2023b: 簡易構文の使用

R2023b 以降では、1 つ以上の名前と値の引数としてオプションを直接指定できます。以下に例を示します。

[sysb,g] = balreal(sys,Goal="relative",...
            Offset=1e-3,FreqIntervals=[10,100]);

この変更の結果として、構文 [sysb,g] = balreal(sys,opts) は推奨されません。ここでは、balredOptions を使用して opts を作成しました。

参考

reducespec | gram | xelim | ssequiv