tfest

伝達関数モデルの推定

構文

sys = tfest(tt,np)

sys = tfest(u,y,np)

sys = tfest(data,np)

sys = tfest(___,nz)

sys = tfest(___,nz,iodelay)

sys = tfest(___,Name,Value)

sys = tfest(tt,init_sys)

sys = tfest(u,y,init_sys)

sys = tfest(data,init_sys)

sys = tfest(___,opt)

[sys,ic] = tfest(___)

説明

伝達関数モデルの推定

sys = tfest(tt,np) は、timetable tt のすべての入力信号と出力信号を使用して、np 個の極をもつ連続時間の伝達関数 sys を推定します。sys の零点の数は max(np-1,0) です。この構文は SISO および MISO システムに使用できます。この関数では、timetable の最後の変数は単一の出力信号であると仮定します。

tfest を使用して時系列モデルを推定することはできません。これは入力を含まないモデルです。代わりに、時系列モデルには ar、arx、または armax を使用します。

sys = tfest(u,y,np) は、行列 u,y 内の時間領域入力信号と出力信号を使用して連続時間の伝達関数を推定します。データのサンプル時間は 1 秒と仮定されます。この仮定されたサンプル時間は変更できません。サンプル時間が 1 秒以外のデータからモデルを推定する場合は、次の 2 つの代替方法があります。

1. 代わりに、'Ts' の名前と値の引数を使用してサンプル時間を設定することで、離散時間システムを推定します。たとえば、sys = tfest(u,y,np,'Ts',0.1) ではサンプル時間が 0.1 に設定されます。この構文は、SISO、MISO、および MIMO システムで使用できます。

2. 連続時間システムを推定する前に、行列データを timetable オブジェクトまたは iddata オブジェクトに変換します。これらの形式では、サンプル時間の知識をデータに組み込むことができます。詳細については、u,y を参照してください。

行列ベースのデータから連続時間モデルを推定することは推奨されません。

sys = tfest(data,np) は、data の時間領域データまたは周波数領域データを使用します。この構文は特に、周波数領域データまたは周波数応答データを使用して伝達関数を推定する場合、またはデータ オブジェクトが提供する追加情報 (サンプル間動作、データのサンプル時間、実験のラベル付けなど) を利用する場合に使用します。

sys = tfest(___,nz) は、零点の数 nz を指定します。この構文では、前述の任意の入力引数の組み合わせで使用できます。

sys = tfest(___,nz,iodelay) は、入力/出力ペアの伝達遅延 iodelay を指定します。

sys = tfest(___,Name,Value) は、1 つ以上の名前と値のペアの引数で指定された追加のモデル オプションを使用します。たとえば、sys = tfest(um,ym,np,'Ts',0.1) を使用して、0.1 のサンプル時間をもつ行列データから離散時間システムを指定します。sys = tfest(data,np,nz,'InputNames',["u1","u2"],'OutputNames',["y1","y3"]) を使用して、MIMO timetable データに使用する変数と一致する入力および出力信号変数名を指定します。

初期パラメーターの構成

sys = tfest(tt,init_sys) は、線形システム init_sys を使用して、timetable tt を使用した推定用に sys の初期パラメーター化を構成します。

sys = tfest(u,y,init_sys) は、推定に行列データ u,y を使用します。init_sys が連続時間モデルの場合は、行列の代わりに timetable を使用することが推奨されます。

sys = tfest(data,init_sys) は、推定にデータ オブジェクト data を使用します。

追加の推定オプションの指定

sys = tfest(___,opt) は、推定の目的、初期条件の処理、正則化、および推定に使用する数値探索法などのオプションを指定するオプション セット opt を組み込みます。opt は、前述の任意の入力引数の組み合わせの後に指定できます。

推定される初期条件を返す

[sys,ic] = tfest(___) は、推定される初期条件を initialCondition オブジェクトとして返します。この構文は、モデルの応答を同じ推定入力データを使用してシミュレートまたは予測し、その応答を同じ推定出力データと比較する場合に使用します。初期条件を組み込むことで、シミュレーションの初期段階において、測定データと、シミュレーション データまたは予測データとの一致率が高まります。

例

極の数の指定による伝達関数モデルの推定

timetable tt1 内の時間領域システム応答データを読み込みます。

load sdata1.mat tt1;

極の数 np を 2 に設定して、伝達関数を推定します。

np = 2;
sys = tfest(tt1,np);

sys は、推定された 2 極の伝達関数を含む idtf モデルです。

結果の推定モデル sys の分子係数と分母係数を確認します。

sys.Numerator

ans = 1×2

    2.4554  176.9856

sys.Denominator

ans = 1×3

    1.0000    3.1625   23.1631

分子と分母の推定における不確かさおよびその他の情報を表示するには、tfdataを使用します。

推定された伝達関数の極と零点の数の指定

時間領域システムの応答データ z2 を読み込み、そのデータを使用して 2 つの極と 1 つの零点を含む伝達関数を推定します。

load iddata2 z2;
np = 2;
nz = 1;
sys = tfest(z2,np,nz);

sys は、推定された伝達関数を含む idtf モデルです。

既知の伝達遅延を含む伝達関数の推定

データ z2 を読み込みます。これは時間領域システム応答データを含む iddata オブジェクトです。

load iddata2 z2;

2 つの極と 1 つの零点を含み、既知の伝達遅延 iodelay を含む伝達関数モデル sys を推定します。

np = 2;
nz = 1;
iodelay = 0.2;
sys = tfest(z2,np,nz,iodelay);

sys は、推定された伝達関数を含み、IODelay プロパティが 0.2 秒に設定された idtf モデルです。

未知の伝達遅延を含む伝達関数の推定

時間領域システムの応答データ z2 を読み込み、そのデータを使用してシステムの 2 極 1 零点の伝達関数を推定します。iodelay の値を NaN に設定して、伝達関数の未知の伝達遅延を指定します。

load iddata2 z2;
np = 2;
nz = 1;
iodelay = NaN;
sys = tfest(z2,np,nz,iodelay);

sys は、推定された伝達関数を含む idtf モデルで、その IODelay プロパティはそのデータを使用して推定されます。

離散時間の伝達関数の推定

時間領域システム応答データを読み込みます。このデータは、入力行列および出力行列である umat2 と ymat2 に含まれます。

load sdata2.mat umat2 ymat2

2 つの極と 1 つの零点をもつ離散時間の伝達関数を推定します。サンプル時間 Ts を 0.1 秒に指定し、伝達遅延 iodelay を 2 秒に指定します。

np = 2;
nz = 1;
iodelay = 2;
Ts = 0.1;
sysd = tfest(umat2,ymat2,np,nz,iodelay,'Ts',Ts)

sysd =
 
  From input "u1" to output "y1":
                     1.8 z^-1
  z^(-2) * ----------------------------
           1 - 1.418 z^-1 + 0.6613 z^-2
 
Sample time: 0.1 seconds
Discrete-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 3
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Estimated using TFEST on time domain data "umat2,ymat2".
Fit to estimation data: 80.26%                          
FPE: 2.095, MSE: 2.063                                  
 
Model Properties

既定では、モデルに直達はなく、推定された伝達関数の分子多項式は最初の係数 b0 が 0 になります。b0 を推定するには、推定中に Feedthrough プロパティを指定します。

直達を含む離散時間の伝達関数の推定

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

load iddata5 z5

最初に、2 つの極と 1 つの零点を含み、直達がない離散時間の伝達関数モデルを推定します。サンプル時間は z5 の Ts プロパティから取得します。

np = 2;
nz = 1;
sys = tfest(z5,np,nz,'Ts',z5.Ts);

推定された伝達関数は次の形式になります。

$$H(z^{-1})=\frac{b1z^{-1}+b2z^{-2}}{1+a1z^{-1}+a2z^{-2}}$$

既定では、モデルに直達はなく、推定された伝達関数の分子多項式は最初の係数 b0 が 0 になります。b0 を推定するには、推定中に Feedthrough プロパティを指定します。

sys = tfest(z5,np,nz,'Ts',z5.Ts,'Feedthrough',true);

推定された伝達関数の分子多項式は、次のように最初の係数が非ゼロになります。

$$H(z^{-1})=\frac{b0+b1z^{-1}+b2z^{-2}}{1+a1z^{-1}+a2z^{-2}}$$

測定データの遅延の原因の解析

直達および伝達遅延があるモデルとないモデルの 2 つの離散時間モデルを比較します。

測定された入力から出力への遅延がある場合は、直達がないこと、または実際の伝送遅延が原因である可能性があります。離散時間モデルの場合、直達がないことは、入力と出力の間の 1 サンプルのラグに相当します。Feedthrough = false と iodelay = 0 を使用してモデルを推定すると、Feedthrough = true と iodelay = 1 を使用して推定されたシステムと等価の離散時間システムが作成されます。どちらのシステムでも、同じ時間領域応答と周波数領域応答がステップ プロットやボード線図などに示されます。ただし、balred を使用してこれらのモデルを低次元化するか、連続時間表現に変換すると、異なる結果になります。したがって、ベスト プラクティスは、観測された遅延の原因が伝達遅延か、または直達がないことかを確認することです。

直達のない離散時間モデルを推定します。

load iddata1 z1
np = 2; 
nz = 2; 
sys1 = tfest(z1,np,nz,'Ts',z1.Ts);

sys1 には直達がなく、したがって $$z^{-1}$$ で始まる分子多項式があるため、sys1 には 1 サンプルのラグがあります。IODelay プロパティは 0 です。

直達をもち、零点を 2 個から 1 個に減らした、1 サンプルの入出力遅延を発生させる別の離散時間を推定します。

sys2 = tfest(z1,np,nz-1,1,'Ts',z1.Ts,'Feedthrough',true);

モデルのボード応答を比較します。

bode(sys1,sys2);

sys1 と sys2 で使用される離散方程式は等価であるため、ボード応答も等価です。

モデルを連続時間に変換し、これらのモデルのボード応答を比較します。

sys1c = d2c(sys1);
sys2c = d2c(sys2);
bode(sys1c,sys2c);
legend

プロットに示されているように、モデルを連続時間に変換すると 2 つのモデルのボード応答は一致しません。sys1c のように直達がない場合は、必ずラグがあります。sys2c のように直達がある場合は、ラグはあり得ません。連続時間の直達は離散時間の直達にマッピングされます。連続時間のラグは離散時間の遅延にマッピングされます。

個々のチャネルに対して直達と遅延が指定された MISO 離散時間の伝達関数の推定

最初の入力に 2 サンプル、2 番目の入力に 0 サンプルの遅延をもつ、2 入力 1 出力の離散時間の伝達関数を推定します。どちらの入力にも直達はありません。

データを読み込み、データを推定データセットと検証データセットに分けます。

load iddata7 z7
ze = z7(1:300);
zv = z7(200:400);

入力から出力への各伝達関数に対して 2 つの極と 1 つの零点をもつ、2 入力 1 出力の伝達関数を推定します。

Lag = [2;0];
Ft = [false,false];
model = tfest(ze,2,1,'Ts',z7.Ts,'Feedthrough',Ft,'InputDelay',Lag);

選択した Feedthrough 値によって、最初の分子係数が 0 (直達なし) かそれ以外か (非ゼロの直達) かが決まります。遅延は一般的に、InputDelay プロパティまたは IODelay プロパティを使用して個別に表されます。この例では InputDelay のみを使用して遅延を表しています。

推定されたモデルを検証します。データ外れ値を検証対象から除外します。

I = 1:201; 
I(114:118) = [];
opt = compareOptions('Samples',I);
compare(zv,model,opt)

正則化されたインパルス応答モデルを使った伝達関数モデルの推定

正則化されたインパルス応答の推定を使用して、15 次伝達関数モデルを同定します。

データを読み込みます。

load regularizationExampleData m0simdata;

正則化されたインパルス応答 (FIR) モデルを取得します。

opt = impulseestOptions('RegularizationKernel','DC');
m0 = impulseest(m0simdata,70,opt);

次数を 15 に減らした後に、モデルを伝達関数モデルに変換します。

m = idtf(balred(idss(m0),15));

モデルの出力をデータと比較します。

compare(m0simdata,m);

推定オプション セットを使用した伝達関数の推定

初期化と探索法を指定する tfest のオプション セットを作成します。また、各反復の損失関数値が表示されることを指定する表示オプションも設定します。

opt = tfestOptions('InitializeMethod','n4sid','Display','on','SearchMethod','lsqnonlin');

時間領域システム応答データ z2 を読み込み、そのデータを使用して、2 つの極と 1 つの零点をもつ伝達関数を推定します。推定オプションに opt を指定します。

load iddata2 z2;
np = 2;
nz = 1;
iodelay = 0.2;
sys = tfest(z2,np,nz,iodelay,opt);

sys は、推定された伝達関数を含む idtf モデルです。

推定された伝達関数のモデルのプロパティの指定

時間領域システム応答データ z2 を読み込み、そのデータを使用して 2 極 1 零点の伝達関数を推定します。入力遅延を指定します。

load iddata2 z2;
np = 2;
nz = 1;
input_delay = 0.2;
sys = tfest(z2,np,nz,'InputDelay',input_delay);

sys は、入力遅延が 0.2 秒の推定された伝達関数を含む idtf モデルです。

周波数応答データの伝達関数への変換

bode を使用して以下のシステムの振幅と位相応答を取得します。

$$H(s)=\frac{s+0.2}{s^3+2s^2+s+1}$$

0.1 ～ 10 rad/s の範囲の 100 個の周波数点を使用して、周波数応答データを取得します。frd を使用して周波数応答データ オブジェクトを作成します。

freq = logspace(-1,1,100);
[mag,phase] = bode(tf([1 0.2],[1 2 1 1]),freq);
data = frd(mag.*exp(1j*phase*pi/180),freq);

data を使用して 3 極 1 零点の伝達関数を推定します。

np = 3;
nz = 1;
sys = tfest(data,np,nz);

sys は、推定された伝達関数を含む idtf モデルです。

複数入力に対する既知の伝達遅延をもつ伝達関数の推定

時間領域システム応答データ co2data を読み込みます。これには、それぞれ 2 つの入力と 1 つ出力をもつ 2 つの実験からのデータが含まれます。1 番目の実験のデータを、サンプル時間 0.5 秒の iddata オブジェクト data に変換します。

load co2data;
Ts = 0.5;
data = iddata(Output_exp1,Input_exp1,Ts);

探索法、および入力オフセットと出力オフセットについて、推定オプションを指定します。探索の反復最大回数も指定します。

opt = tfestOptions('SearchMethod','gna');
opt.InputOffset = [170;50];
opt.OutputOffset = mean(data.y(1:75));
opt.SearchOptions.MaxIterations = 50;

測定データと推定オプション セット opt を使用して伝達関数を推定します。入力から出力への伝達関数を指定します。

np = 3;
nz = 1;
iodelay = [2 5];
sys = tfest(data,np,nz,iodelay,opt);

iodelay は、1 番目と 2 番目の入力から出力からの、入力から出力への遅延をそれぞれ 2 秒および 5 秒として指定します。

sys は、推定された伝達関数を含む idtf モデルです。

既知および未知の伝達遅延をもつ伝達関数モデルの推定

時間領域システムの応答データを読み込み、そのデータを使用してシステムの伝達関数を推定します。既知および未知の伝達遅延を指定します。

load co2data;
Ts = 0.5;
data = iddata(Output_exp1,Input_exp1,Ts);

data は、2 つの入力チャネルと 1 つの出力チャネルをもち、サンプル レートが 0.5 秒の iddata オブジェクトです。

オプション セット opt を作成します。探索法、および入力オフセットと出力オフセットについて、推定オプションを指定します。探索の反復最大回数も指定します。

opt = tfestOptions('Display','on','SearchMethod','gna');
opt.InputOffset = [170; 50];
opt.OutputOffset = mean(data.y(1:75));
opt.SearchOptions.MaxIterations = 50;

2 秒の既知の遅延を表す 2 と未知の遅延を表す nan を使用して、iodelay の未知および既知の伝達遅延を指定します。iodelay および opt を使用して伝達関数を推定します。

np = 3;
nz = 1;
iodelay = [2 nan];
sys = tfest(data,np,nz,iodelay,opt);

sys は、推定された伝達関数を含む idtf モデルです。

未知の制約付き伝達遅延をもつ伝達関数の推定

予想される分子および分母構造と遅延制約をもつ伝達関数モデルを作成します。

この例では、実験データは 2 つの入力と 1 つの出力で構成されます。どちらの伝達遅延も未知で、同一の上限をもちます。さらに、どちらの入力から出力への伝達関数も構造は同一です。

init_sys = idtf(NaN(1,2),[1,NaN(1,3)],'IODelay',NaN);
init_sys.Structure(1).IODelay.Free = true;
init_sys.Structure(1).IODelay.Maximum = 7;

init_sys は、1 つの入力から出力への伝達関数の構造を記述する idtf モデルです。伝達関数は 1 つの零点、3 つの極、および伝達遅延で構成されます。NaN の使用は、未知の係数を示します。

init_sys.Structure(1).IODelay.Free = true は、伝達遅延が固定ではないことを示します。

init_sys.Structure(1).IODelay.Maximum = 7 は、伝達遅延の上限を 7 秒に設定します。

両方の入力から出力への伝達関数を指定します。

init_sys = [init_sys,init_sys];

時間領域システム応答データを読み込み、そのデータを使用して伝達関数を推定します。tfestOptions オプション セット opt でオプションを指定します。

load co2data;
Ts = 0.5; 
data = iddata(Output_exp1,Input_exp1,Ts);
opt = tfestOptions('Display','on','SearchMethod','gna');
opt.InputOffset = [170;50];
opt.OutputOffset = mean(data.y(1:75));
opt.SearchOptions.MaxIterations = 50;
sys = tfest(data,init_sys,opt);

sys は、推定された伝達関数を含む idtf モデルです。

比較により推定結果を解析します。compareOptions オプション セット opt2 を作成し、入力オフセットと出力オフセットを指定してから、compare を使用します。

opt2 = compareOptions;
opt2.InputOffset = opt.InputOffset;
opt2.OutputOffset = opt.OutputOffset;
compare(data,sys,opt2)

入出力ペアに対して異なる数の極を含む伝達関数の推定

与えられたデータの入出力のペアに対して異なる数の極を含む、複数入力、単出力の伝達関数を推定します。

周波数応答データを取得します。

たとえば、frd を使用して、以下のシステムの周波数応答データ モデルを作成します。

$$G=\left[\begin{array}{c}{e}^{-4s}\frac{s+2}{s^3+2s^2+4s+5}\\ e^{-0.6s}\frac{5}{s^4+2s^3+s^2+s}\end{array}\right]$$

0.01 rad/s ～ 100 rad/s の範囲の 100 個の周波数点を使用して、周波数応答データを取得します。

G = tf({[1 2],[5]},{[1 2 4 5],[1 2 1 1 0]},0,'IODelay',[4 0.6]);
data = frd(G,logspace(-2,2,100));

data は、G の連続時間周波数応答を含む frd オブジェクトです。

data の伝達関数を推定します。

 np = [3 4];
 nz = [1 0];
 iodelay = [4 0.6];
 sys = tfest(data,np,nz,iodelay);

np は、推定された伝達関数の極の数を指定します。np の 1 番目の要素は、1 番目の入力から出力への伝達関数に 3 つの極が含まれることを示します。同様に、np の 2 番目の要素は、2 番目の入力から出力への伝達関数に 4 つの極が含まれていることを示します。

nz は、推定された伝達関数の零点の数を指定します。nz の 1 番目の要素は、1 番目の入力から出力への伝達関数に 1 つの零点が含まれることを示します。同様に、np の 2 番目の要素は、2 番目の入力から出力への伝達関数に零点が含まれていないことを示します。

iodelay は、1 番目の入力から出力への伝達遅延を 4 秒として指定します。2 番目の入力から出力への伝達遅延は 0.6 秒として指定されます。

sys は、推定された伝達関数を含む idtf モデルです。

不安定システムの伝達関数の推定

周波数応答データを使用して、不安定システムを記述する伝達関数を推定します。

idtf を使用して、以下のシステムの伝達関数モデル G を作成します。

$$G=\left[\begin{array}{c}\frac{s+2}{s^3+2s^2+4s+5}\\ \frac{5}{s^4+2s^3+s^2+s+1}\end{array}\right]$$

G = idtf({[1 2], 5},{[1 2 4 5],[1 2 1 1 1]});

idfrd を使用して、G の周波数応答データ モデル data を作成します。0.01 rad/s ～ 100 rad/s の範囲の 100 個の周波数点を指定します。

data = idfrd(G,logspace(-2,2,100));

data は idfrd オブジェクトです。

data の伝達関数を推定します。

np = [3 4];
nz = [1 0];
sys = tfest(data,np,nz);

np は、推定された伝達関数の極の数を指定します。np の 1 番目の要素は、1 番目の入力から出力への伝達関数に 3 つの極が含まれることを示します。同様に、np の 2 番目の要素は、2 番目の入力から出力への伝達関数に 4 つの極が含まれていることを示します。

nz は、推定された伝達関数の零点の数を指定します。nz の 1 番目の要素は、1 番目の入力から出力への伝達関数に 1 つの零点が含まれることを示します。同様に、nz の 2 番目の要素は、2 番目の入力から出力への伝達関数に零点が含まれていないことを示します。

sys は、推定された伝達関数を含む idtf モデルです。

pole(sys)

ans = 7×1 complex

  -1.5260 + 0.0000i
  -0.2370 + 1.7946i
  -0.2370 - 1.7946i
  -1.4656 + 0.0000i
  -1.0000 + 0.0000i
   0.2328 + 0.7926i
   0.2328 - 0.7926i

sys は、極の表示に示されているように、不安定システムです。

モーダル密度が高い周波数応答データを使用した伝達関数の推定

高密度の周波数応答測定データを読み込みます。このデータは、フィードバック制御を使用して平衡に維持される不安定処理に対応します。

load HighModalDensityData FRF f

同定のためにデータを idfrd オブジェクトとしてパッケージ化し、ボード振幅応答を求めます。

G = idfrd(permute(FRF,[2 3 1]),f,0,'FrequencyUnit','Hz');
bodemag(G)

32 個の極と 32 個の零点をもつ伝達関数を推定し、ボード振幅応答を比較します。

sys = tfest(G,32,32);
bodemag(G, sys)
xlim([0.01,2e3])
legend

推定された初期条件の取得と適用

データを読み込んでプロットします。

load iddata1ic z1i
plot(z1i)

出力データ y(1) の初期値を確認します。

ystart = z1i.y(1)

ystart = 
-3.0491

測定された出力は 0 から始まりません。

2 次伝達関数 sys を推定し、推定された初期条件 ic を返します。

[sys,ic] = tfest(z1i,2,1);
ic

ic = 
  initialCondition with properties:

     A: [2×2 double]
    X0: [2×1 double]
     C: [0.2957 5.2441]
    Ts: 0

ic は、sys の自由応答を状態空間形式で X0 の初期状態ベクトルにカプセル化する initialCondition オブジェクトです。

初期条件を組み込まずに、推定データを使用して sys をシミュレートします。シミュレートした出力を測定出力とともにプロットします。

y_no_ic = sim(sys,z1i);
plot(y_no_ic,z1i)
legend('Model Response','Output Data')

測定出力とシミュレートした出力は、シミュレーションの開始時には一致していません。

初期条件を simOptions オプション セットに組み込みます。

opt = simOptions('InitialCondition',ic);
y_ic = sim(sys,z1i,opt);
plot(y_ic,z1i)
legend('Model Response','Output Data')

シミュレーションでは、入力信号に対するモデル応答が、初期条件に対する自由応答と結合されます。測定出力とシミュレートした出力は、シミュレーションの開始時に以前より一致するようになります。この初期条件は、推定データ z1i に対してのみ有効です。

入力引数

tt — timetable ベースの推定データ

timetable | timetables の cell 配列

推定データ。入力信号変数と出力信号変数の両方を含む等間隔でサンプリングされた timetable として、または、複数実験データの場合は timetable の cell 配列として指定します。

timetable 全体の使用

tt 内のすべての入力変数と出力変数を使用し、一連の入力変数の後に一連の出力変数が続くように変数が構成されている場合、次のように指定します。

• SISO システムの場合は、tt を N_s 行 2 列の timetable として指定します。N_s はサンプルの数で、2 つの timetable 変数はそれぞれ測定された入力信号と出力信号を表します。

• MIMO システムの場合は、tt を N_s 行 (N_u+N_y) 列の timetable として指定します。N_u は入力の数、N_y は出力の数です。最初の N_u 個の変数には入力信号が含まれ、残りの N_y 個の変数には出力信号が含まれなければなりません。

  状態空間モデルまたは伝達関数モデルを推定している場合は、次の節で説明されているように、入力チャネルと出力チャネルも明示的に指定しなければなりません。

• 複数実験データの場合は、データを timetable の N_e 行 1 列の cell 配列として指定します。N_e は実験数です。すべての実験のサンプル時間が一致していなければなりません。

timetable から選択した変数の使用

timetable からの変数のサブセットを使用する場合、または入力変数と出力変数が混在している場合は、'InputName' および 'OutputName' 名前と値の引数を使用して、使用する変数を指定します。

たとえば、tt に 6 つの変数 "u1"、"u2"、"u3"、および "y1"、"y2"、"y3" が含まれているとします。変数 "u1" と "u2" を入力、変数 "y1" と "y3" を出力として推定に使用します。推定を実行するには、次のコマンドを使用します。

sys = tfest(tt,__,'InputName',["u1" "u2"],'OutputName',["y1" "y3"])

推定データ型の取り扱いの詳細については、Data Domains and Data Types in System Identification Toolboxを参照してください。

u, y — 行列ベースの推定データ

行列 | 行列の cell 配列

推定データ。SISO システムの場合は、等間隔でサンプリングされた入力および出力時間領域信号値を含む、N_s 行 1 列の実数値の行列のペアとして指定します。ここで、N_s はサンプルの数です。

MIMO システムの場合は、u,y を、以下の次元をもつ入力行列と出力行列のペアとして指定します。

• u — N_s 行 N_u 列の行列。N_u は入力の数です。

• y — N_s 行 N_y 列の行列。N_y は出力の数です。

複数実験データの場合は、u,y を 1 行 N_e 列の cell 配列のペアとして指定します。N_e は実験数です。すべての実験のサンプル時間が一致していなければなりません。

制限

• 行列ベースのデータは、周波数領域データからの推定をサポートしません。iddata オブジェクトや idfrd オブジェクトなどのデータ オブジェクト (data を参照) を使用する必要があります。

• 推定データはサンプル時間を提供しないため、連続時間推定では推定データに行列を使用することは推奨されません。データは 1 Hz でサンプリングされると仮定されます。連続時間推定の場合、入力および出力行列のペアを単一の timetable に変換することが推奨されます。たとえば、1 列の行列 um および ym をサンプル時間が 0.5 分の timetable tt に変換するには、次のコマンドを使用します。

  tt = timetable(um,ym,'rowtimes',minutes(0.5*(1:size(u,1))))

  行列ベースの SISO データを timetable に変換する詳細な例については、timetable への SISO 行列データの変換を参照してください。MIMO 行列のペアを timetable に変換する例については、連続時間モデルを推定するために MIMO 行列データを timetable に変換を参照してください。

  推定データ型の取り扱いの詳細については、Data Domains and Data Types in System Identification Toolboxを参照してください。

data — 推定データ オブジェクト

iddata オブジェクト | frd オブジェクト | idfrd オブジェクト

推定データ オブジェクト。iddata オブジェクト、frd オブジェクト、または等間隔にサンプリングされた入力値と出力値を含む idfrd オブジェクトとして指定します。既定では、モデルのサンプル時間は推定データのサンプル時間に設定されます。

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

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

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

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

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

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

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

  • Domain — 'Frequency'

制限

離散時間周波数領域データを使用して連続時間モデルを推定することはできません。

np — 極の数

非負の整数 | 行列

推定された伝達関数の極の数。非負の整数または行列として指定します。

複数の入力または複数の出力、あるいはその両方を含むシステムの場合は、次のように、np のグローバル値または個々の値を入力/出力のペアに適用できます。

• すべてのペアの極が同数 — np をスカラーとして指定します。

• ペアごとに極の数が異なる — np を N_y 行 N_u 列の行列として指定します。N_y は出力の数、N_u は入力の数です。

例については、極の数の指定による伝達関数モデルの推定を参照してください。

nz — 零点の数

非負の整数 | 行列

推定された伝達関数の零点の数。非負の整数または行列として指定します。

複数の入力または複数の出力、あるいはその両方を含むシステムの場合は、次のように、nz のグローバル値または個々の値を入力/出力のペアに適用できます。

• すべてのペアの極が同数 — nz をスカラーとして指定します。

• ペアごとに極の数が異なる — nz を N_y 行 N_u 列の行列として指定します。N_y は出力の数、N_u は入力の数です。

離散時間データを使用して推定される連続時間モデルの場合は、nz <= np を設定します。

離散時間モデルの推定の場合は、nz を伝達関数の分子多項式の零点の数として指定します。たとえば、tfest(tt,2,1,'Ts',data.Ts) は $$b_1{z}^{-1}/(1+a_1{z}^{-1}+b_2{z}^{-2})$$ 形式の伝達関数を推定しますが、tfest(tt,2,2,'Ts',data.Ts) は $$(b_1{z}^{-1}+b_2{z}^{-2})/(1+a_1{z}^{-1}+b_2{z}^{-2})$$ を推定します。ここで z^-1 は Z 変換ラグ変数です。離散時間の伝達関数の詳細については、離散時間表現を参照してください。例については、離散時間の伝達関数の推定を参照してください。

iodelay — 伝達遅延

[] (既定値) | 非負の整数 | 行列

伝達遅延。非負の整数、NaN スカラー、または行列として指定します。

連続時間システムの場合、data の TimeUnit プロパティに格納された時間単位で伝達遅延を指定します。離散時間システムの場合、サンプル時間 Ts の倍数の遅延を示す整数として伝達遅延を指定します。

N_y 個の出力と N_u 個の入力がある MIMO システムの場合は、iodelay を N_y 行 N_u 列の配列に設定します。この配列の各エントリは、対応する入力/出力のペアの伝達遅延を表す数値です。iodelay をスカラー値に設定して、同一の遅延をすべての入力/出力のペアに適用することもできます。

指定された値は固定遅延として処理されます。未知の伝達遅延を示すには、iodelay 行列で NaN を指定します。

[] または 0 を使用して、伝達遅延がないことを示します。

例については、既知の伝達遅延を含む伝達関数の推定を参照してください。

opt — 推定のオプション

tfestOptions オプション セット

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

• 推定の目的

• 初期条件の処理

• 推定で使用される数値探索法

• サンプル間動作

例については、推定オプション セットを使用した伝達関数の推定を参照してください。

init_sys — sys の初期パラメーター化を構成する線形システム

idtf モデル | 線形モデル | 構造体

sys の初期パラメーター化を構成する線形システム。idtf モデル、線形モデル、または構造体として指定します。init_sys は、測定データを使用した推定を実行して取得することも、直接構築することもできます。

init_sys が idtf モデルの場合、tfest は init_sys のパラメーター値を sys の推定の初期推定として使用します。

init_sys の Structure プロパティを使用して、分子、分母、および伝達ラグの初期パラメーター値と制約を構成します。たとえば次のようにします。

• init_sys の A 行列に対する初期推定を指定するには、init_sys.Structure.Numerator.Value を初期推定に設定します。

• init_sys の B 行列に対する制約を指定するには、次のように設定します。

  • init_sys.Structure.Numerator.Minimum を最小の分子係数値に設定します。

  • init_sys.Structure.Numerator.Maximum を最大の分子係数値に設定します。

  • init_sys.Structure.Numerator.Free を設定して、推定についてどの分子係数が自由であるかを示します。

  例については、未知の制約付き伝達遅延をもつ伝達関数の推定を参照してください。

init_sys が idtf モデルではない場合、最初に init_sys が伝達関数に変換されます。tfest は、結果のモデルのパラメーターを推定の初期推定として使用します。

opt を指定せず、init_sys が作成されたのではなく推定で取得されたものである場合、init_sys.Report.OptionsUsed からの推定オプションが使用されます。

名前と値の引数

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

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

例: sys = tfest(data,np,nz,'Ts',0.1)

InputName — 入力チャネル名

string | 文字ベクトル | string 配列 | 文字ベクトルの cell 配列

入力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、InputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = tfest(tt,__,'InputName',["u1" "u2"]) は、推定で使用する timetable tt からの入力チャネルとして変数 u1 および u2 を選択します。

OutputName — 出力チャネル名

string | 文字ベクトル | string 配列 | 文字ベクトルの cell 配列

出力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、OutputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = tfest(tt,__,'OutputName',["y1" "y3"]) は、推定で使用する timetable tt からの出力チャネルとして変数 y1 および y3 を選択します。

Ts — 推定モデルのサンプル時間

0 (連続時間) (既定値) | 正のスカラー

推定されたモデルのサンプル時間。0 または正のスカラーとして指定します。

• 連続時間モデルの場合、'Ts' を 0 と指定します。

• 離散時間モデルの場合、'Ts' を次のように定義される単位でデータ サンプル時間として指定します。

  • timetable ベースのデータ — timetable Time 列

  • 行列ベースのデータ — 秒

  • データ オブジェクト (iddata オブジェクトなど) — data.TimeUnit プロパティ

  離散の場合、np と nz は、分子多項式と分母多項式の z^-1 の根の数を表します。

  timetable tt のデータのサンプル時間を取得するには、timetable プロパティ tt.Properties.Timestep を使用します。

例については、離散時間の伝達関数の推定を参照してください。

InputDelay — 入力遅延

0 (既定値) | スカラー | ベクトル

各入力チャネルの入力遅延。スカラーまたは数値ベクトルとして指定します。

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

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

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

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

例については、推定された伝達関数のモデルのプロパティの指定を参照してください。

Feedthrough — 離散時間の伝達関数の直達

0 (既定値) | 1 | logical 行列

離散時間の伝達関数の直達。logical スカラーまたは N_y 行 N_u 列の logical 行列として指定します。N_y は出力の数、N_u は入力の数です。すべての入力/出力チャネルに同じ直達を使用するには、Feedthrough をスカラーとして指定します。

2 つの極と 3 つの零点をもつ離散時間モデルについて考えます。

$$H(z^{-1})=\frac{b0+b1z^{-1}+b2z^{-2}+b3z^{-3}}{1+a1z^{-1}+a2z^{-2}}$$

モデルに直達が含まれる場合、b0 は、値が残りのモデル パラメーター b1、b2、b3、a1、および a2 とともに推定される自由パラメーターです。モデルに直達がない場合、b0 はゼロに固定されます。例については、直達を含む離散時間の伝達関数の推定を参照してください。

出力引数

sys — 同定された伝達関数

idtf モデル

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

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

Report のフィールド	説明
Status	モデルのステータスの概要。モデルが構築によって作成されたものか推定によって取得されたものかを示します。
Method	使用された推定コマンド
InitializeMethod	時間領域データを使用した連続時間の伝達関数の推定で、分子と分母を初期化するために使用されたアルゴリズム。次の値のいずれかとして返されます。 'iv' — 操作変数法 'svf' — 状態変数フィルター アプローチ 'gpmf' — 一般化ポアソン モーメント関数アプローチ 'n4sid' — 部分空間状態空間推定アプローチ このフィールドは、推定オプション セットで InitializeMethod オプションが 'all' に設定されている場合に使用されたアルゴリズムを確認するのに特に便利です。
N4Weight	InitializeMethod が 'n4sid' の場合に、特異値分解ステップで使用された重みアルゴリズム。次のいずれの値として返されます。 'MOESP' — Verhaegen による MOESP アルゴリズムを使用。 'CVA' — Larimore による正準変量アルゴリズム (CVA) を使用。 'SSARX' — ARX の推定に基づくアルゴリズムを使用して重み付けを計算する部分空間同定法を使用。 このフィールドは、推定オプション セットで N4Weight オプションが 'auto' に設定されている場合に使用された重み行列を確認するのに便利です。
N4Horizon	InitializeMethod が 'n4sid' である場合に使用された前方および後方の予測範囲。3 要素の行ベクトル [r sy su] として返されます。ここで、r は前方予測範囲、sy は予測に使用された過去の出力の数、su は予測に使用された過去の入力の数です。
InitialCondition	モデル推定時の初期条件の処理。次の値のいずれかとして返されます。 'zero' — 初期条件をゼロに設定。 'estimate' — 初期条件を独立した推定パラメーターとして処理。 'backcast' — 初期条件を最適な最小二乗近似を使用して推定。 このフィールドは、推定オプション セットで InitialCondition オプションが 'auto' に設定されている場合に、初期条件がどのように扱われたかを確認するのに特に便利です。
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	推定に使用されたオプション セット。これは、カスタム オプションを構成していない場合は既定のオプションのセットになります。詳細については、polyestOptions を参照してください。
RandState	推定開始時の乱数ストリームの状態。推定時にランダム化が使用されなかった場合は空 [] になります。詳細については、rng を参照してください。
DataUsed	推定に使用されたデータの属性。次のフィールドをもつ構造体として返されます。 Name — データ セットの名前 Type — データ型 Length — データ サンプルの数 Ts — サンプル時間 InterSample — 入力サンプル間動作。次の値のいずれかとして返されます。 'zoh' — ゼロ次ホールドにより、サンプル間で区分的に一定な入力信号を維持。 'foh' — 1 次ホールドにより、サンプル間で区分的に線形な入力信号を維持。 'bl' — 帯域幅を制限した動作により、ナイキスト周波数を超える連続時間入力信号のパワーがゼロになるように指定。 InputOffset — 推定時に時間領域の入力データから削除されたオフセット。非線形モデルの場合は [] になります。 OutputOffset — 推定時に時間領域の出力データから削除されたオフセット。非線形モデルの場合は [] になります。
Termination	予測誤差の最小化に使用された反復探索の終了条件。次のフィールドを含む構造体として返されます。 WhyStop — 数値探索を終了する理由 Iterations — 推定アルゴリズムで実行された探索の反復回数 FirstOrderOptimality — 探索アルゴリズムが終了したときの勾配探索ベクトルの $$\infty$$ ノルム FcnCount — 目的関数が呼び出された回数 UpdateNorm — 最後の反復における勾配探索ベクトルのノルム。探索法が 'lsqnonlin' または 'fmincon' の場合は省略されます。 LastImprovement — 最後の反復における基準改善。パーセンテージで表されます。探索法が 'lsqnonlin' または 'fmincon' の場合は省略されます。 Algorithm — 'lsqnonlin' または 'fmincon' の探索法で使用されるアルゴリズム。他の探索法が使用されている場合は省略されます。 数値探索の最適化を必要としない推定法の場合、Termination フィールドは省略されます。

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

ic — 初期条件

initialCondition オブジェクト | initialCondition 値のオブジェクト配列

推定される初期条件。initialCondition オブジェクトまたは initialCondition 値のオブジェクト配列として返されます。

• 単一実験データ セットの場合、ic は、推定される初期条件 (x₀) に対する伝達関数モデルの自由応答 (行列 A および C) を状態空間形式で表します。

• N_e 回の実験用の複数実験データ セットの場合、ic は、実験ごとに 1 つの initialCondition 値のセットを含む長さ N_e のオブジェクト配列です。

tfest が 0 の ic 値を返し、非ゼロの初期条件があることがわかっている場合は、tfestOptions の 'InitialCondition' オプションを 'estimate' に設定し、更新されたオプション セットを tfest に渡します。以下に例を示します。

opt = tfestOptions('InitialCondition','estimate')
[sys,ic] = tfest(data,np,nz,opt)

'InitialCondition' の既定の 'auto' 設定では、初期条件が推定誤差最小化プロセス全体に無視できる影響をもたらす場合、'zero' 手法が使用されます。'estimate' を指定すると、ic の値が推定されます。

詳細については、initialCondition を参照してください。この引数の使用例については、推定された初期条件の取得と適用を参照してください。

アルゴリズム

tfest で使用される推定アルゴリズムの詳細は、推定されたモデルおよび推定データのサンプリングを含むさまざまな要因によって異なります。

時間領域データを使用した連続時間の伝達関数の推定

パラメーターの初期化

推定アルゴリズムは、InitializeMethod 推定オプションで指定された方法を使用して推定可能なパラメーターを初期化します。既定の方法は、操作変数 (IV) 法です。

連続時間パラメーター推定の状態変数フィルター (SVF) アプローチと一般化ポアソン モーメント関数 (GPMF) アプローチは、事前にフィルターしたデータを使用します[1][2]。[1]および[2]の定数 $$\frac{1}{\lambda }$$ は、初期化オプション (InitializeOptions) のフィールド FilterTimeConstant に対応します。IV は簡易調整 IV 法で、[3]では SRVIC と呼ばれています。この方法には、SVF で初期化された、現在のモデルの分母となるプレフィルターがあります。モデルの変更が Tolerance より小さくなるまで、このプレフィルターは最大 MaxIterations 回繰り返されます。MaxIterations と Tolerance は、InitializeOptions 構造体を使用して指定できるオプションです。'n4sid' 初期化オプションは、N4SID 推定アルゴリズムを使用して離散モデルを推定し、d2c を使って連続時間に変換します。

tfestOptions を使用して、伝達関数の推定に使用されるオプション セットを作成します。

パラメーターの更新

初期化されたパラメーターは、SearchMethod 推定オプションで指定された非線形最小二乗探索法を使用して更新されます。この探索法の目的は、重み付き予測誤差ノルムを最小化することです。

時間領域データを使用した離散時間の伝達関数の推定

離散時間データの場合、tfest は oe と同じアルゴリズムを使用して、分子および分母多項式係数を決定します。このアルゴリズムでは、arx を使用して初期化が実行され、その後に非線形最小二乗探索ベースの更新で重み付き予測誤差ノルムを最小化します。

連続時間周波数領域データを使用した連続時間の伝達関数の推定

推定アルゴリズムは以下のタスクを実行します。

1. 双一次マッピングを実行して伝達関数の領域 (周波数グリッド) を変換します。連続時間モデルの場合、虚軸は単位円板に変換されます。離散時間モデルの場合、元の領域の単位円板が別の単位円板に変換されます。

2. S-K 反復[4]を実行して、非線形最小二乗問題を解きます。多入力単出力システムを検討します。非線形最小二乗問題は、以下の損失関数を最小化することです。

   $$\underset{D,N_i}{\text{minimize}}{\displaystyle \sum _{k=1}^{n_f}{\left|W({\omega }_k)\left(y({\omega }_k)-{\displaystyle \sum _{i=1}^{n_u}\frac{N_i({\omega }_k)}{D({\omega }_k)}{u}_i({\omega }_k)}\right)\right|}^2}$$

   ここで、W は、指定する周波数依存の重みです。D は推定される伝達関数モデルの分母で、N_i は i 番目の入力に対応する分子です。y および u はそれぞれ測定された出力データと入力データです。n_f および n_u は周波数と入力の数で、w は周波数です。項を再配列すると、次のようになります。

   $$\underset{D,N_i}{\text{minimize}}{\displaystyle \sum _{k=1}^{n_f}{\left|\frac{W({\omega }_k)}{D({\omega }_k)}\left(D({\omega }_k)y({\omega }_k)-{\displaystyle \sum _{i=1}^{n_u}{N}_i({\omega }_k)u_i({\omega }_k)}\right)\right|}^2}$$

   S-K 反復を実行するために、アルゴリズムは反復的に以下を求解します。

   $$\underset{D_m,N_{i,m}}{\text{minimize}}{\displaystyle \sum _{k=1}^{n_f}{\left|\frac{W({\omega }_k)}{D_{m-1}({\omega }_k)}\left(D_m({\omega }_k)y({\omega }_k)-{\displaystyle \sum _{i=1}^{n_u}{N}_{i,m}({\omega }_k)u_i({\omega }_k)}\right)\right|}^2}$$

   ここで、m は現在の反復で、D_m-1(ω) は前の反復で特定された分母応答です。反復の各ステップは線形最小二乗問題です。ここで、特定されたパラメーターは i = 1,2,...n_u に対する応答 D_m(ω) および N_i,m(ω) を取得します。D₀(ω) = 1 を選択することで反復が初期化されます。

   • アルゴリズムの最初の反復によって D₁(ω) が特定されます。D₁(ω) および N_i,1(ω) 多項式は単項式ベースで表されます。

   • 2 回目以降の反復は、単位円板上の直交有理基底関数の項で多項式 D_m(ω) および N_i,m(ω) を表します。これらの基底関数の形式は次のとおりです。

     $$B_{j,m}(\omega )=\left(\frac{\sqrt{1-{\left|{\lambda }_{j,m-1}\right|}^2}}{q-{\lambda }_{j,m-1}}\right){\displaystyle \prod _{r=0}^{j-1}\frac{1-{({\lambda }_{j,m-1})}^{*}q(\omega )}{q(\omega )-{\lambda }_{r,m-1}}}$$

     ここで、λ_j,m-1 は、反復の前のステップ m-1 で特定される j 番目の極です。λ_j,m-1^* は λ_j,m-1 の複素共役で、q は単位円板上の周波数領域変数です。

   • アルゴリズムは最大 20 回の反復に対して実行されます。最後の 3 回の反復で損失関数の値の相対変化が 0.001 より小さい場合、反復は早期に終了します。

   伝達関数係数の範囲を指定する場合、これらの範囲は特定されたパラメーターのアフィン制約に相当します。等式制約 (固定伝達関数係数) のみある場合、対応する等式制約付き最小二乗問題は代数的に求解されます。このために、等式制約行列のヌル空間の直交基底が計算され、このヌル空間内の最小二乗問題が求解されます。伝達関数係数に上限または下限がある場合、対応する不等式制約付きの最小二乗問題は内点法を使用して求解されます。

3. 線形調整を実行します。S-K 反復は、収束する場合でも、常に局所的に最適な解を示すわけではありません。局所的に最適な解を示す可能性がある最適化問題の臨界点を求めるには、2 回目の反復のセットが実行されます。臨界点は、一連の非線形方程式に対する解です。アルゴリズムは、非線形方程式に対する線形近似を連続的に作成し、結果の線形方程式を最小二乗法で解くことで臨界点を探索します。方程式は次のようになります。

   • j 番目の分子パラメーターの方程式:

     $$0=2{\displaystyle \sum _{k=1}^{n_f}\mathrm{Re}\left\{\frac{{\left|W({\omega }_k)\right|}^2{B}_j^{*}({\omega }_k){\displaystyle \sum _{i=1}^{n_u}{N}_{i,m-1}^{*}({\omega }_k)u_i^{*}({\omega }_k)}}{D_{m-1}^{*}({\omega }_k){\left|D_{m-1}({\omega }_k)\right|}^2}\left(D_m({\omega }_k)y({\omega }_k)-{\displaystyle \sum _{i=1}^{n_u}{N}_{i,m}({\omega }_k)u_i({\omega }_k)}\right)\right\}}$$

   • 入力 l に対応する j 番目の分子パラメーターの方程式:

     $$0=-2{\displaystyle \sum _{k=1}^{n_f}\mathrm{Re}\left\{\frac{{\left|W({\omega }_k)\right|}^2{B}_j^{*}({\omega }_k)u_l^{*}({\omega }_k)}{{\left|D_{m-1}({\omega }_k)\right|}^2}\left(D_m({\omega }_k)y({\omega }_k)-{\displaystyle \sum _{i=1}^{n_u}{N}_{i,m}({\omega }_k)u_i({\omega }_k)}\right)\right\}}$$

   最初の反復は、S-K 反復時に分子 N_i パラメーターと分母 D パラメーターで求めた最適解から開始されます。S-K 反復とは異なり、基底関数 B_j(ω) は各反復で変更されません。反復は、S-K 反復で最適解を出した基底関数を使用して実行されます。前回と同様に、アルゴリズムは最大 20 回の反復に対して実行されます。最後の 3 回の反復で損失関数の値の相対変化が 0.001 より小さい場合、反復は早期に終了します。

   伝達関数係数の範囲を指定すると、これらの範囲は一般化ラグランジュ乗数を使用して必要な最適性条件に組み込まれます。結果の制約付き線形最小二乗問題は、S-K 反復ステップで説明した方法と同じ方法で求解されます。

4. 最適解に対応する伝達関数パラメーターを返します。S-K 反復および線形調整反復ステップはどちらも損失関数値の改善を保証しません。アルゴリズムはこれらのステップ中に観測された最適なパラメーター値を追跡し、これらの値を返します。

5. ステップ 1 で実行した双一次マッピングの逆を実行します。

6. SearchMethod 推定オプションで指定された非線形最小二乗探索法を使用して、伝達関数パラメーターの反復調整を実行します。このステップは、以下の状況で実装されます。

   • EnforceStability 推定オプションを true (安定性が要求される) として指定し、このアルゴリズムのステップ 5 の結果が不安定な場合。不安定な極は安定範囲内に反映され、結果のパラメーターは繰り返し調整されます。推定オプションの詳細については、tfestOptions を参照してください。

   • Regularization 推定オプションを使用して損失関数に正則化のペナルティを追加する場合。正則化の例については、動的システムの正則化による同定を参照してください。

   • 離散時間データを使用して連続時間モデルを推定する (離散時間周波数領域データを使用した離散時間の伝達関数の推定を参照)。

   • 周波数領域の入力/出力データを使用して多入力モデルを同定する。

連続時間の周波数領域データを使用して連続時間モデルを推定するために R2016a 以前の推定アルゴリズム (tfest 推定アルゴリズムの更新を参照) を使用する場合、連続時間データと固定遅延については、モデルの推定に出力誤差アルゴリズムが使用されます。連続時間データと自由遅延については、状態空間推定アルゴリズムが使用されます。このアルゴリズムでは、N4SID 推定法を使用してモデル係数が初期化されます。この初期化の後に、非線形最小二乗探索ベースの更新で重み付き予測誤差ノルムが最小化されます。

離散時間周波数領域データを使用した離散時間の伝達関数の推定

推定アルゴリズムは、離散時間データを使用すること以外、連続時間周波数領域データを使用した連続時間の伝達関数の推定と同じです。

R2016a 以前の推定アルゴリズム (tfest 推定アルゴリズムの更新を参照) を使用している場合、アルゴリズムは時間領域データを使用した離散時間の伝達関数の推定用アルゴリズムと同じです。

メモ

連続時間周波数領域データを使用した離散時間の伝達関数の推定はサポートされていません。

離散時間周波数領域データを使用した連続時間の伝達関数の推定

tfest コマンドは、最初に離散時間データから離散時間モデルを推定します。推定されたモデルは、d2c コマンドを使用して連続時間モデルに変換されます。次に、結果の連続時間モデルの周波数応答が、推定データの周波数グリッドで計算されます。次に、目的の (ユーザー指定の) 構造の連続時間モデルはこの周波数応答に適合されます。周波数応答データを使用して連続時間モデルを取得するための推定アルゴリズムは、連続時間データを使用した連続時間の伝達関数の推定のアルゴリズムと同じです。

R2016a 以前の推定アルゴリズム (tfest 推定アルゴリズムの更新を参照) を使用している場合、離散時間データから連続時間モデルを推定するために状態空間推定アルゴリズムが使用されます。このアルゴリズムでは、N4SID 推定法を使用してモデル係数が初期化されます。この初期化の後に、非線形最小二乗探索ベースの更新で重み付き予測誤差ノルムが最小化されます。

遅延の推定

• 遅延値が NaN として指定される場合、delayest を使用して、モデルの分子係数と分母係数から個別に遅延値が推定されます。モデルの残りの部分の反復更新中に、tfest はこれらの遅延値を固定値として処理します。そのため、遅延値が繰り返し更新されることはありません。

• 既定では、離散時間データ (Ts>0) の場合、delayest は遅延の探索を 0 ～ 30 サンプルの範囲に制限します。連続時間モデルの場合、この範囲は 0 ～ 30 Ts 時間単位に変換されます。連続時間データ (Ts = 0) の場合、delayest は探索範囲を 0 ～ 10 時間単位に制限します。これらの範囲を変更するには、最初に idtf を使用してテンプレート モデル init_sys を作成し、次に、init_sys.Structure.IODelay.Minimum および init_sys.Structure(i,j).IODelay.Maximum の値を設定します。

• 以下をもつ初期モデル init_sys の場合:

  • 有限値として指定される init_sys.Structure.IODelay.Value

  • true として指定される init_sys.Structure.IODelay.Free

  初期の遅延値は変更されません。

遅延の推定の多くは困難な問題です。ベスト プラクティスは、遅延の有無と値を評価することです。これを行うには、モデル化されているプロセスの物理的な洞察や、arxstruc、delayest、impulseest などの関数を使用します。入力遅延を決定する例については、モデル構造の選択: モデル次数と入力遅延の決定を参照してください。

参照

[1] Garnier, H., M. Mensler, and A. Richard. “Continuous-Time Model Identification from Sampled Data: Implementation Issues and Performance Evaluation.” International Journal of Control 76, no. 13 (January 2003): 1337–57. https://doi.org/10.1080/0020717031000149636.

[2] Ljung, Lennart. “Experiments with Identification of Continuous Time Models.” IFAC Proceedings Volumes 42, no. 10 (2009): 1175–80. https://doi.org/10.3182/20090706-3-FR-2004.00195.

[3] Young, Peter, and Anthony Jakeman. “Refined Instrumental Variable Methods of Recursive Time-Series Analysis Part III. Extensions.” International Journal of Control 31, no. 4 (April 1980): 741–64. https://doi.org/10.1080/00207178008961080.

[4] Drmač, Z., S. Gugercin, and C. Beattie. “Quadrature-Based Vector Fitting for Discretized H₂ Approximation.” SIAM Journal on Scientific Computing 37, no. 2 (January 2015): A625–52. https://doi.org/10.1137/140961511.

[5] Ozdemir, Ahmet Arda, and Suat Gumussoy. “Transfer Function Estimation in System Identification Toolbox via Vector Fitting.” IFAC-PapersOnLine 50, no. 1 (July 2017): 6232–37. https://doi.org/10.1016/j.ifacol.2017.08.1026.

バージョン履歴

R2012a で導入

R2022b: 時間領域推定データが timetable 形式および行列形式で受け入れられる

ほとんどの推定関数、検証関数、解析関数、およびユーティリティ関数が、時間領域入出力データを、入力と出力両方のデータを含む単一の timetable の形式、または、入力データと出力データを別々に含む行列のペアの形式で受け入れるようになりました。これらの関数は、時間領域データと周波数領域データの両方について、iddata オブジェクトを引き続きデータ ソースとして受け入れます。

R2018a: SearchMethod が 'lsqnonlin' の場合、SearchOptions の詳細オプションは非推奨

Optimization Toolbox™ で lsqnonlin 探索法 (ソルバー) を使用して推定する場合の並列処理を呼び出すオプションを含む、lsqnonlin 関連の詳細オプションの指定は推奨されません。

R2016b: tfest 推定アルゴリズムの更新

R2016b 以降では、周波数領域データから伝達関数の推定を実行するために新しいアルゴリズムが使用されています。新しいアルゴリズムを使うと、特に周波数や振幅の範囲が大きいダイナミクスをもつデータの場合は、迅速かつより正確な結果を得ることができます。ただし、推定結果は前のリリースの結果と一致しない場合があります。以前の推定アルゴリズムを使用して推定を実行するには、'-R2016a' を構文に追加します。

たとえば、周波数領域データ data を使用して、np 個の極をもつ伝達関数モデルを推定しているとします。

sys = tfest(data,np)

以前の推定アルゴリズムを使用するには、次の構文を使用します。

sys = tfest(data,np,'-R2016a')

参考

tfestOptions | idtf | timetable | ssest | procest | ar | arx | oe | bj | polyest | greyest

トピック

• コマンド ラインでの伝達関数モデルの推定
• Estimate Transfer Function Models with Transport Delay to Fit Given Frequency-Response Data
• Estimate Transfer Function Models with Prior Knowledge of Model Structure and Constraints
• Apply Initial Conditions When Simulating Identified Linear Models
• Troubleshoot Frequency-Domain Identification of Transfer Function Models
• 伝達関数モデルとは
• Regularized Estimates of Model Parameters
• Estimating Models Using Frequency-Domain Data