このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。
Simulink データを使用した連続時間モデルの推定
この例では、System Identification Toolbox™ を使用して、Simulink® でシミュレートされたモデルを同定する方法を説明します。この例では、連続時間システムおよび遅延を処理する方法と、入力のサンプル間動作の重要性について説明します。
Simulink モデルからのシミュレーション データの取得
以下の Simulink モデルで記述されるシステムを見てみましょう。
open_system('iddemsl1') set_param('iddemsl1/Random Number','seed','0')
赤い部分はシステム、青い部分はコントローラー、基準信号はスイープ正弦波 (チャープ信号) です。データのサンプル時間は 0.5 秒に設定されます。
このシステムは、idpoly 構造を使用して表されます。
m0 = idpoly(1,0.1,1,1,[1 0.5],'Ts',0,'InputDelay',1,'NoiseVariance',0.01)
m0 =
Continuous-time OE model: y(t) = [B(s)/F(s)]u(t) + e(t)
B(s) = 0.1
F(s) = s + 0.5
Input delays (listed by channel): 1
Parameterization:
Polynomial orders: nb=1 nf=1 nk=0
Number of free coefficients: 2
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
iddemsl1 モデルをシミュレートし、そのデータを iddata オブジェクトに保存します。
sim('iddemsl1') dat1e = iddata(y,u,0.5); % The IDDATA object
検証のため、もう一度シミュレーションを実行します。
set_param('iddemsl1/Random Number','seed','13') sim('iddemsl1') dat1v = iddata(y,u,0.5);
最初のシミュレーションで取得した推定データを見てみましょう。
plot(dat1e)
シミュレーション データを使用した離散モデルの推定
まず最初に、既定の次数の離散モデルを評価し、データ特性の予備情報を取得します。
m1 = n4sid(dat1e, 'best') % A default order model
m1 =
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
x1 0.7881 0.1643 0.1116
x2 -0.1214 0.4223 -0.8489
x3 0.155 0.7527 0.2119
B =
u1
x1 -0.0006427
x2 -0.02218
x3 0.07347
C =
x1 x2 x3
y1 -5.591 0.871 1.189
D =
u1
y1 0
K =
y1
x1 -0.001856
x2 0.002363
x3 -0.06805
Sample time: 0.5 seconds
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using N4SID on time domain data "dat1e".
Fit to estimation data: 86.17% (prediction focus)
FPE: 0.01281, MSE: 0.01251
モデルがどれぐらい検証データで再現されているかを確認します。
compare(dat1v,m1)
ご覧のように、検証データは、モデルで正しく予測されています。データ特性をさらに詳しく評価するため、解析の負のラグが自動的に決定される dat1e を使用して算出したノンパラメトリックなインパルス応答を検査します。
ImpModel = impulseest(dat1e,[],'negative');
clf
h = impulseplot(ImpModel);
showConfidence(h,3)
ImpModel は次数 (係数の数) が自動的に決定される FIR モデルです。また、ここでは負の遅延に対するインパルス応答を計算し、フィードバックの影響を解析します。負のラグによる影響は無視できるものとは限りません。これは、制御器によるものです (出力フィードバック)。これは、インパルス応答推定を、時間遅延を決定するときに使用できないことを意味します。その代わりに、複数の低次数 ARX モデルを異なる遅延で構築し、最適な適合を求めます。
V = arxstruc(dat1e,dat1v,struc(1:2,1:2,1:10));
nn = selstruc(V,0) %delay is the third element of nn
nn =
2 2 3
遅延は、3 ラグとなります。(これは正しい値です。むだ時間 1 秒で 2 ラグ遅延となり、残りの 1 つは ZOH ブロックです)。また、対応する ARX モデルは、以下のように算出できます。
m2 = arx(dat1e,nn) compare(dat1v,m1,m2);
m2 =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
A(z) = 1 - 0.2568 z^-1 - 0.3372 z^-2
B(z) = 0.04021 z^-3 + 0.04022 z^-4
Sample time: 0.5 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nk=3
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "dat1e".
Fit to estimation data: 85.73% (prediction focus)
FPE: 0.01346, MSE: 0.0133
推定の改良
2 つのモデル m1 および m2 は、シミュレーションで同様の動作を示します。ここで、次数と遅延を微調整します。遅延を 2 に修正し (直達がないことに加え、最終的な遅延は 3 サンプルになり)、この遅延の値が設定されている既定の次数の状態空間モデルを検出します。
m3 = n4sid(dat1e,'best','InputDelay',2,'Feedthrough',false); % Refinement for prediction error minimization using pem (could also use % |ssest|) m3 = pem(dat1e, m3);
推定されたシステム行列を見てみます。
m3.a % the A-matrix of the resulting model
ans =
0.7332 -0.3784 0.1735
0.4705 0.3137 -0.6955
-0.0267 0.7527 0.6343
3 次ダイナミクスが自動的に選択されます。ここに「追加の」遅延 2 を加えると、5 次状態空間モデルとなります。
自動次数選択に頼りすぎることはお勧めできません。選択は、無作為の誤差に影響されます。最適な方法は、モデルの零点と極を信頼領域と合わせて評価することです。
clf
h = iopzplot(m3);
showConfidence(h,2) % Confidence region corresponding to 2 standard deviations
単位円周での 2 極/零点は明らかにキャンセルされると思われ、1 次のダイナミクスで十分であることがわかります。この情報を使用して、新たに 1 次推定を行います。
m4 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',2,'Ts',dat1e.Ts); compare(dat1v,m4,m3,m1)
compare プロットには、シンプルな 1 次モデル m4 が良好にデータを記述していることが示されています。したがって、このモデルを最終結果として選択します。
離散時間モデルの連続時間モデルへの変換 (LTI)
このモデルを連続時間に変換し、伝達関数の形式で表します。
mc = d2c(m4); idtf(mc)
ans =
From input "u1" to output "y1":
0.09828
exp(-1*s) * ----------
s + 0.4903
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
前述のように、システムの良好な記述が得られます。
連続時間モデルの直接推定
連続時間モデルは直接推定することもできます。離散モデル m4 には 2 サンプル入力遅延があります。これは 1 秒の遅延を表します。ここでは、ssest コマンドを使用してこの推定を実行します。
m5 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1); present(m5)
m5 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1
x1 -0.4903 +/- 0.007948
B =
u1
x1 0.01345 +/- 3.94e+10
C =
x1
y1 7.307 +/- 2.141e+13
D =
u1
y1 0
K =
y1
x1 -0.02227 +/- 6.525e+10
Input delays (seconds): 1
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 4
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 7, Number of function evaluations: 15
Estimated using SSEST on time domain data "dat1e".
Fit to estimation data: 87.34% (prediction focus)
FPE: 0.01054, MSE: 0.01047
More information in model's "Report" property.
不確かさの解析
モデルはデータに 87% 適合していますが、モデル m5 のパラメーターは高水準の不確かさを示しています。これは、モデルが実際に必要なパラメーターより多くのパラメーターを使用しているため、パラメーターの推定の一意性が維持されなくなったことが原因です。モデルの不確かさの実際の影響を表示する方法は 2 通りあります。
不確かさをパラメーターではなくモデルの応答の信頼限界として表示する。
モデルを正準形で推定する。
ここでは両方の方法を試してみましょう。まず、モデルを正準形で推定します。
m5Canon = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); present(m5Canon)
m5Canon =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1
x1 -0.4903 +/- 0.007881
B =
u1
x1 0.09828 +/- 0.001559
C =
x1
y1 1
D =
u1
y1 0
K =
y1
x1 -0.1628 +/- 0.03702
Input delays (seconds): 1
Parameterization:
CANONICAL form with indices: 1.
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 3
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Termination condition: Near (local) minimum, (norm(g) < tol)..
Number of iterations: 7, Number of function evaluations: 15
Estimated using SSEST on time domain data "dat1e".
Fit to estimation data: 87.34% (prediction focus)
FPE: 0.01054, MSE: 0.01047
More information in model's "Report" property.
m5Canon はモデルの正準パラメーター化を使用します。推定データとの適合はモデル m5 と同等です。パラメーターの値の不確かさが小さいことから、信頼性が高いことがわかります。しかし、m5 の場合と同様に、不確かさの値が大きくても "良くない" モデルであるとは限りません。これらのモデルの品質を確実なものとするため、標準偏差 3 に対応する信頼領域を使用して、時間領域および周波数領域での応答を表示します。また、元のシステム m0 を比較のためプロットします。
ボード線図です。
clf opt = bodeoptions; opt.FreqScale = 'linear'; h = bodeplot(m0,m5,m5Canon,opt); showConfidence(h,3) legend show
ステップ プロットです。
clf
showConfidence(stepplot(m0,m5,m5Canon),3)
legend show
2 つのモデルの不確かさの範囲はほぼ同じです。同様に、これらのモデルの信頼領域を使用して、極-零点配置図 (iopzplot) とナイキスト線図 (nyquistplot) を作成することもできます。
idtf(m5)
ans =
From input "u1" to output "y1":
0.09828
exp(-1*s) * ----------
s + 0.4903
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
連続時間推定におけるサンプル間動作の考慮
サンプル データから算出した連続時間モデルを比較する場合、入力信号のサンプル間動作を考慮することが重要です。この例ではここまで、コントローラの Zero-order-Hold (zoh) 回路により、システムへの入力は区分的定数となっています。この回路を取り外し、実際の連続システムを評価します。入力および出力信号は、やはり 2 Hz でサンプリングされますが、その他はすべて同じです。
open_system('iddemsl3') sim('iddemsl3') dat2e = iddata(y,u,0.5);
測定値に合わせて調整されるときに、サンプリング プロパティとサンプル間入力動作 (現在の入力) が盛り込まれるため、離散時間モデルはこれらのデータでも良好です。ただし、連続時間モデルを構築する場合、サンプル間プロパティを把握しておくことが必須となります。まず、ZOH 向けにモデルを構築します。
m6 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); idtf(m6)
ans =
From input "u1" to output "y1":
0.1117
exp(-1*s) * ----------
s + 0.5596
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
推定モデル (m6) を実モデル (m0) と比較します。
step(m6,m0) % Compare with true system
結果は良好に一致しています。ただし、入力に関するデータ オブジェクト情報を含めることができます。近似として、各サンプリングの瞬間の間は、区分的線形 (First-order-hold、FOH) として記述します。この情報は、正規のサンプリングの推定器で使われます。
dat2e.Intersample = 'foh'; m7 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); % new estimation with correct intersample behavior idtf(m7)
ans =
From input "u1" to output "y1":
0.09937
exp(-1*s) * ----------
s + 0.4957
Continuous-time identified transfer function.
Parameterization:
Number of poles: 1 Number of zeros: 0
Number of free coefficients: 2
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by conversion from idss model.
もう一度ステップ応答の比較を見てください。
step(m7,m0) % Compare with true system
このモデル (m7) では、m6 よりもかなり良好な結果が得られています。これでこの例は終了です。
bdclose('iddemsl1'); bdclose('iddemsl3');
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)