Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

n4sid

時間領域または周波数領域のデータを使用した部分空間法による状態空間モデルの推定

説明

sys = n4sid(data,nx) は、時間領域または周波数領域のデータ data を使用して、次数 nx の離散時間状態空間モデル sys を推定します。sys は次の形式のモデルになります。

x˙(t)=Ax(t)+Bu(t)+Ke(t)y(t)=Cx(t)+Du(t)+e(t)

A、B、C、D、および K は状態空間行列です。u(t) は入力、y(t) は出力、e(t) は外乱、x(t) は nx 個の状態のベクトルです。

A、B、C、および K のエントリは、既定ではすべて推定可能な自由パラメーターです。動的システムの場合、D は既定ではゼロに固定されます。これは、システムに直達がないことを意味します。静的システム (nx = 0) の場合、D は既定では推定可能なパラメーターです。

sys = n4sid(data,nx,Name,Value) は、1 つ以上の名前と値のペアの引数で指定された追加のオプションを組み込みます。たとえば、連続時間モデルを推定するには、サンプル時間 'Ts'0 と指定します。A、B、C、D、および K の各行列の既定の動作を変更するには、名前と値のペアの引数 'Form''Feedthrough''DisturbanceModel' を使用します。

sys = n4sid(___,opt) は、推定オプション opt を指定します。これらのオプションには、初期状態、推定の目的、推定に使用する部分空間アルゴリズムに関する選択が含まれます。opt は、前述の任意の入力引数の組み合わせの後に指定します。

[sys,x0] = n4sid(___) は、推定時に計算された初期状態の値を返します。この構文では、前述の任意の入力引数の組み合わせで使用できます。

すべて折りたたむ

状態空間モデルを推定し、その応答を測定出力と比較します。

iddata オブジェクトに格納されている入出力データ z1 を読み込みます。

load iddata1 z1

4 次状態空間モデルを推定します。

nx = 4;
sys = n4sid(z1,nx)
sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
              x1         x2         x3         x4
   x1     0.8392    -0.3129    0.02105    0.03743
   x2     0.4768     0.6671     0.1428   0.003757
   x3   -0.01951    0.08374   -0.09761      1.046
   x4  -0.003885   -0.02914    -0.8796   -0.03171
 
  B = 
              u1
   x1    0.02635
   x2   -0.03301
   x3  7.256e-05
   x4  0.0005861
 
  C = 
            x1       x2       x3       x4
   y1    69.08    26.64   -2.237  -0.5601
 
  D = 
       u1
   y1   0
 
  K = 
              y1
   x1   0.003282
   x2   0.009339
   x3  -0.003232
   x4   0.003809
 
Sample time: 0.1 seconds
  
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 28
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using N4SID on time domain data "z1".  
Fit to estimation data: 76.33% (prediction focus)
FPE: 1.21, MSE: 1.087                            

シミュレートしたモデルの応答を測定出力と比較します。

compare(z1,sys)

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent z1 (y1), sys: 70.44%.

プロットから、シミュレートしたモデルと推定データの間の適合率が 70% を超えていることがわかります。

推定に関する詳しい情報は、idsssys.Report プロパティを調べて確認できます。

sys.Report
ans = 
          Status: 'Estimated using N4SID with prediction focus'
          Method: 'N4SID'
    InitialState: 'estimate'
        N4Weight: 'CVA'
       N4Horizon: [6 10 10]
             Fit: [1x1 struct]
      Parameters: [1x1 struct]
     OptionsUsed: [1x1 idoptions.n4sid]
       RandState: [1x1 struct]
        DataUsed: [1x1 struct]

たとえば、推定された初期状態に関する詳しい情報を調べます。

sys.Report.Parameters.X0
ans = 4×1

   -0.0085
    0.0052
   -0.0193
    0.0282

iddata オブジェクトに格納されている入出力データ z1 を読み込みます。

load iddata1 z1

引数 nx を 1 から 10 までの範囲に指定して、最適なモデル次数を調べます。

nx = 1:10;
sys = n4sid(z1,nx);

自動生成されるプロットに、nx で指定した次数のモデルに対するハンケル特異値が表示されます。

ハンケル特異値が比較的小さい状態については破棄してかまいません。提案される既定の次数の選択は 2 です。

[選択した次数] のリストでモデル次数を選択し、[適用] をクリックします。

推定データを読み込みます。

load iddata2 z2

推定オプションを指定します。重み付けスキーム 'N4Weight''SSARX' に設定し、推定ステータス表示オプション 'Display''on' に設定します。

opt = n4sidOptions('N4Weight','SSARX','Display','on')
Option set for the n4sid command:

          InitialState: 'estimate'
              N4Weight: 'SSARX'
             N4Horizon: 'auto'
               Display: 'on'
           InputOffset: []
          OutputOffset: []
    EstimateCovariance: 1
          OutputWeight: []
                 Focus: 'prediction'
       WeightingFilter: []
      EnforceStability: 0
              Advanced: [1x1 struct]

更新したオプション セットを使用して、3 次状態空間モデルを推定します。

nx = 3;
sys = n4sid(z2,nx,opt);

A、B、および C の各行列の正準形を変更し、D 行列に直達項を含め、K 行列の外乱モデルの推定を排除します。

入出力データを読み込み、n4sid の既定のオプションを使用して 4 次システムを推定します。

load iddata1 z1
sys1 = n4sid(z1,4);

モード形式を指定し、A 行列を既定の A 行列と比較します。

sys2 = n4sid(z1,4,'Form','modal');
A1 = sys1.A
A1 = 4×4

    0.8392   -0.3129    0.0211    0.0374
    0.4768    0.6671    0.1428    0.0038
   -0.0195    0.0837   -0.0976    1.0462
   -0.0039   -0.0291   -0.8796   -0.0317

A2 = sys2.A
A2 = 4×4

    0.7554    0.3779         0         0
   -0.3779    0.7554         0         0
         0         0   -0.0669    0.9542
         0         0   -0.9542   -0.0669

直達項を含め、D 行列を比較します。

sys3 = n4sid(z1,4,'Feedthrough',1);
D1 = sys1.D
D1 = 0
D3 = sys3.D
D3 = 0.0487

外乱のモデル化を排除し、K 行列を比較します。

sys4 = n4sid(z1,4,'DisturbanceModel','none');
K1 = sys1.K
K1 = 4×1

    0.0033
    0.0093
   -0.0032
    0.0038

K4 = sys4.K
K4 = 4×1

     0
     0
     0
     0

連続時間正準形モデルを推定します。

推定データを読み込みます。

load iddata1 z1

モデルを推定します。Ts0 に設定して連続モデルを指定します。

nx = 2;
sys = n4sid(z1,nx,'Ts',0,'Form','canonical');

sys は、正準形の 2 次連続時間状態空間モデルです。

部分空間アルゴリズム SSARX を使用して閉ループ データから状態空間モデルを推定します。このアルゴリズムは、フィードバックの影響を取得するのに他の重み付けアルゴリズムよりも優れています。

ホワイト ノイズで破損した 2 次システムの閉ループ推定データを生成します。

N = 1000; 
K = 0.5;
rng('default');
w = randn(N,1);
z = zeros(N,1); 
u = zeros(N,1); 
y = zeros(N,1);
e = randn(N,1);
v = filter([1 0.5],[1 1.5 0.7],e);
for k = 3:N
   u(k-1) = -K*y(k-2) + w(k);
   u(k-1) = -K*y(k-1) + w(k);
   z(k) = 1.5*z(k-1) - 0.7*z(k-2) + u(k-1) + 0.5*u(k-2);
   y(k) = z(k) + 0.8*v(k);
end
dat = iddata(y, u, 1);

N4SID アルゴリズムで使用される重み付けスキーム 'N4weight' を指定します。2 つのオプション セットを作成します。一方のオプション セットで、'N4weight''CVA' に設定します。もう一方のオプション セットで、'N4weight''SSARX' に設定します。

optCVA = n4sidOptions('N4weight','CVA');
optSSARX = n4sidOptions('N4weight','SSARX');

オプション セットを使用して状態空間モデルを推定します。

sysCVA = n4sid(dat,2,optCVA);
sysSSARX = n4sid(dat,2,optSSARX);

2 つのモデルの推定データとの適合を比較します。

compare(dat,sysCVA,sysSSARX);

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent dat (y1), sysCVA: 71.26%, sysSSARX: 77.1%.

プロットからわかるように、SSARX アルゴリズムを使用して推定したモデルの方が CVA アルゴリズムを使用して推定したモデルよりも適合しています。

入力引数

すべて折りたたむ

推定データ。iddata オブジェクト、frd オブジェクト、または idfrd オブジェクトとして指定します。

時間領域の推定の場合、data は入力信号と出力信号の値を含む iddata オブジェクトでなければなりません。

周波数領域の推定の場合、data は次のいずれかになります。

  • 記録された周波数応答データ (frd (Control System Toolbox) または idfrd)

  • プロパティが次のように指定された iddata オブジェクト。

    • InputData — 入力信号のフーリエ変換

    • OutputData — 出力信号のフーリエ変換

    • Domain'Frequency'

推定データは等間隔サンプルでなければなりません。既定では、モデルのサンプル時間は推定データのサンプル時間に設定されます。

複数実験データの場合、すべての実験のサンプル時間とサンプル間動作が一致していなければなりません。

データの領域に応じて、推定できるモデルのタイプが決まります。

  • 時間領域データまたは離散時間の周波数領域データ — 連続時間モデルと離散時間モデル

  • 連続時間の周波数領域データ — 連続時間モデルのみ

推定モデルの次数。非負の整数、または正の整数の範囲を含むベクトルとして指定します。

  • 推定モデルの目的の次数が既にわかっている場合は、nx をスカラーとして指定します。

  • 推定モデルの最も効果的な次数を選択するために可能性がある次数の範囲を比較する場合は、その範囲を nx に指定します。n4sid により、システムにおけるそれぞれの状態の相対的なエネルギーの寄与を示すハンケル特異値のプロットが作成されます。ハンケル特異値が比較的小さい状態は、モデルの精度にほとんど寄与せず、破棄しても影響はほとんどありません。維持する状態のうちで最も高い状態のインデックスがモデル次数になります。プロット ウィンドウには、使用する次数の提案も示されます。この提案を受け入れることも、別の次数を入力することもできます。例については、推定による最適なモデル次数の特定を参照してください。

    nx を指定しない場合、または nxbest と指定した場合、範囲 1:10 から nx が自動的に選択されます。

  • 静的システムを同定する場合は、nx0 に設定します。

推定オプション。n4sidOptions オプション セットとして指定します。opt で指定するオプションには以下が含まれます。

  • 推定の目的

  • 初期条件の処理

  • 部分空間アルゴリズムに関する選択

オプションを指定する方法の例については、推定オプションの指定および連続時間正準形モデルを参照してください。

名前と値の引数

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

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

例: sys = n4sid(data,nx,'Form','modal')

推定モデルのサンプル時間。'Ts'0 または正のスカラーのいずれかで構成されるコンマ区切りのペアとして指定します。

  • 連続時間モデルの場合、'Ts'0 と指定します。周波数領域では、連続時間の周波数領域データを使用すると連続時間モデルになります。

  • 離散時間モデルの場合、'Ts'TimeUnit プロパティに格納された単位でのデータのサンプル時間に設定されます。

各入力チャネルの入力遅延。'InputDelay' と数値ベクトルで構成されるコンマ区切りのペアとして指定します。

  • 連続時間モデルの場合、TimeUnit プロパティに格納された時間単位で 'InputDelay' を指定します。

  • 離散時間モデルの場合、サンプル時間 Ts の整数倍で 'InputDelay' を指定します。たとえば、'InputDelay'3 に設定すると、3 サンプリング周期の遅延が指定されます。

Nu 個の入力があるシステムの場合、InputDelay を Nu 行 1 列のベクトルに設定します。このベクトルの各エントリは、対応する入力チャネル用の入力遅延を表す数値です。

すべてのチャネルに同じ遅延を適用するには、'InputDelay' をスカラーとして指定します。

sys の正準形のタイプ。'Form' と次の値のいずれかで構成されるコンマ区切りのペアとして指定します。

  • 'free' — 行列 A、B、C、D、および K のすべてのエントリを自由なものとして扱います。

  • 'modal'sys をモード形式で取得します。

  • 'companion'sys をコンパニオン形式で取得します。

  • 'canonical'sys を可観測性の正準形で取得します。

正準形の定義については、正準状態空間実現を参照してください。

これらの形式を同定に使用する方法の詳細については、Estimate State-Space Models with Canonical Parameterizationを参照してください。

例については、形式、直達、外乱モデルの行列の変更を参照してください。

入力から出力への直達。'Feedthrough' と長さ Nu の logical ベクトルで構成されるコンマ区切りのペアとして指定します。Nu は入力の数です。'Feedthrough' を logical スカラーとして指定すると、その値がすべての入力に適用されます。静的システムの場合、'Feedthrough' は常に 1 であると仮定されます。

例については、形式、直達、外乱モデルの行列の変更を参照してください。

K 行列で時間領域のノイズ成分のパラメーターを推定するオプション。'DisturbanceModel' と次の値のいずれかで構成されるコンマ区切りのペアとして指定します。

  • 'estimate' — ノイズ成分を推定します。K 行列は自由パラメーターとして扱われます。時間領域データの場合、'estimate' が既定です。

  • 'none' — ノイズ成分を推定しません。K 行列の要素はゼロで固定されます。周波数領域データの場合、'none' が既定であり、受け入れられる唯一の値です。

例については、形式、直達、外乱モデルの行列の変更を参照してください。

出力引数

すべて折りたたむ

同定された状態空間モデル。idss モデルとして返されます。このモデルは、指定したモデル次数、遅延、および推定オプションを使用して作成されます。

推定結果と使用されたオプションに関する情報は、モデルの Report プロパティに格納されます。Report には次のフィールドがあります。

Report のフィールド説明
Status

モデルのステータスの概要。モデルが構築によって作成されたものか推定によって取得されたものかを示します。

Method

使用された推定コマンド。

InitialState

推定時の初期状態の処理方法。次の値のいずれかとして返されます。

  • 'zero' — 初期状態をゼロに設定。

  • 'estimate' — 初期状態を独立した推定パラメーターとして処理。

このフィールドは、推定オプション セットで 'InitialState' オプションが 'auto' に設定されている場合に特に便利です。

N4Weight

N4SID アルゴリズムで特異値分解に使用された重み付けスキーム。次の値のいずれかとして返されます。

  • 'MOESP' — MOESP アルゴリズムを使用。

  • 'CVA' — 正準変量アルゴリズムを使用。

  • 'SSARX' — 部分空間同定法で ARX の推定に基づくアルゴリズムを使用して重み付けを計算。

このオプションは、推定オプション セットで N4Weight オプションが 'auto' に設定されている場合に特に便利です。

N4Horizon

N4SID アルゴリズムで使用された前方および後方の予測範囲。3 要素の行ベクトル [r sy su] として返されます。ここで、r は前方予測範囲、sy は予測に使用された過去の出力の数、su は予測に使用された過去の入力の数です。

Fit

推定の定量的評価。構造体として返されます。これらの品質メトリクスの詳細については、Loss Function and Model Quality Metricsを参照してください。構造体には、以下のフィールドがあります。

フィールド説明
FitPercent

正規化平方根平均二乗誤差 (NRMSE)。モデルの応答が推定データにどの程度適合するかをパーセンテージで示す尺度で、fitpercent = 100(1-NRMSE) として表されます。

LossFcn

推定完了時の損失関数の値。

MSE

平均二乗誤差 (MSE)。モデルの応答が推定データにどの程度適合するかを示す尺度です。

FPE

モデルの最終予測誤差。

AIC

生の赤池情報量基準 (AIC)。モデルの品質を示す尺度です。

AICc

小さいサンプルサイズの補正された AIC。

nAIC

正規化された AIC。

BIC

ベイズ情報量基準 (BIC)。

Parameters

モデル パラメーターの推定値。

OptionsUsed

推定に使用されたオプション セット。OptionsUsed は、カスタム オプションを構成していない場合は既定のオプションのセットになります。詳細については、n4sidOptions を参照してください。

RandState

推定開始時の乱数ストリームの状態。推定時にランダム化が使用されなかった場合は空 [] になります。詳細については、rng を参照してください。

DataUsed

推定に使用されたデータの属性。次のフィールドをもつ構造体として返されます。

フィールド説明
Name

データ セットの名前。

Type

データ型。

Length

データ サンプルの数。

Ts

サンプル時間。

InterSample

入力サンプル間動作。次の値のいずれかとして返されます。

  • 'zoh' — ゼロ次ホールドにより、サンプル間で区分的に一定な入力信号を維持。

  • 'foh' — 1 次ホールドにより、サンプル間で区分的に線形な入力信号を維持。

  • 'bl' — 帯域幅を制限した動作により、ナイキスト周波数を超える連続時間入力信号のパワーがゼロになるように指定。

InputOffset

推定時に時間領域の入力データから削除されたオフセット。非線形モデルの場合は [] になります。

OutputOffset

推定時に時間領域の出力データから削除されたオフセット。非線形モデルの場合は [] になります。

Report の使用の詳細については、Estimation Reportを参照してください。

推定時に計算された初期状態。各実験に対応する列ベクトルを含む配列として返されます。

この配列は、モデルの Report プロパティの Parameters フィールドにも格納されます。

参照

[1] Ljung, L. System Identification: Theory for the User, Appendix 4A, Second Edition, pp. 132–134. Upper Saddle River, NJ: Prentice Hall PTR, 1999.

[2] van Overschee, P., and B. De Moor. Subspace Identification of Linear Systems: Theory, Implementation, Applications. Springer Publishing: 1996.

[3] Verhaegen, M. "Identification of the deterministic part of MIMO state space models." Automatica, 1994, Vol. 30, pp. 61–74.

[4] Larimore, W.E. "Canonical variate analysis in identification, filtering and adaptive control." Proceedings of the 29th IEEE Conference on Decision and Control, 1990, pp. 596–604.

[5] McKelvey, T., H. Akcay, and L. Ljung. "Subspace-based multivariable system identification from frequency response data." IEEE Transactions on Automatic Control, 1996, Vol. 41, pp. 960–979.

バージョン履歴

R2006a より前に導入