このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
非線形 ARX モデルと Hammerstein-Wiener モデルの同定のチュートリアル
この例では、測定した入出力データを使用して、単入力単出力 (SISO) の非線形ブラック ボックス モデルの同定を説明します。ここでは 2 つのタンクからなるシステムの測定データを使用して、さまざまなモデル構造と同定の選択肢について調べます。
入出力データセット
この例では、twotankdata.mat
ファイルに保存されているデータ変数を使用します。これには、0.2 秒のサンプル時間で生成された 2 タンク システムの入出力データ サンプルが 3000 件含まれています。入力 u(t) はポンプに印加された電圧 [V] で、上部タンクへの流入を生成します。この上部タンクの底部にある小さな穴から流出し、下部タンクに流れ込みます。2 タンク システムの出力 y(t) は下部タンクの水位 [m] です。読み込まれたデータをカプセル化する IDDATA オブジェクト z
を作成します。
このデータセットは非線形グレー ボックスの例「2 タンク システム: 連続時間 SISO システムの C MEX ファイル モデリング」でも使用され、そこではシステムを記述する物理方程式が仮定されています。この例は、ブラック ボックス モデルを対象としているため、実際の物理法則に関する知識は想定していません。
load twotankdata z = iddata(y, u, 0.2, 'Name', 'Two tank system');
データの分割とプロット
このデータを、それぞれが 1000 サンプルを含む、同じサイズの 3 つのサブセットに分割します。
ここで、モデルを z1
で推定し、z2
および z3
で評価します。z1
は z2
と似ていますが、z3
には似ていません。したがってある意味では、z2
はモデルの内挿機能の評価に役立ち、z3
は外挿機能の評価に役立ちます。
z1 = z(1:1000); z2 = z(1001:2000); z3 = z(2001:3000); plot(z1,z2,z3) legend('Estimation','Validation 1', 'Validation 2')
リグレッサーの選択
非線形ブラック ボックス モデルを推定する最初のステップは、リグレッサー セットの選択です。リグレッサーは時間遅延入力/出力変数に基づいた単純な式であり、その最も簡単なケースは連続値の小さなセットによって遅延した変数です。たとえば、入力変数の名前が "u"、出力変数の名前が "y" の場合、リグレッサー セットの例は { y(t-1), y(t-2), u(t), u(t-1), u(t-2) } になります。ここで "t" は時間変数を表します。多項式の項を含むもう 1 つの例は、{abs(y(t-2)), y(t-2)*u(t-1), u(t-4)^2} になります。遅延変数には、より複雑な任意の式も可能です。
System Identification Toolbox™ では、リグレッサー セットを明示的に作成する必要はありません。形式 (多項式の次数など)、寄与する変数、ラグを指定するだけで、大規模な集合を暗黙的に生成できます。仕様オブジェクト linearRegressor
、polynomialRegressor
および customRegressor
はこの目的で使用されます。たとえば R = linearRegressor({'y','u'},{[1 2 4], [2 10]}) は、リグレッサー セット {y(t-1), y(t-2), y(t-4), u(t-2), u(t-10)} を意味します。同様に、R = polynomialRegressor({'y'},{[1 2]},3) は、ラグが 1 と 2 である変数 'y' の 3 次多項式リグレッサーを意味します。すなわち、セット {y(t-1)^3, y(t-2)^3} です。絶対値のみを使用する、リグレッサー式にラグの混合を許可するなどの構成もさらに選択できます。
次数行列を使用した線形リグレッサー仕様
連続ラグに基づいた線形リグレッサーのみをモデルが使用する場合、一般的に、linearRegressor
オブジェクトを使用するより、次数行列を使用して指定する方が簡単です。次数行列 [na nb nk] は、リグレッサー式に使用する過去の出力、過去の入力および入力遅延の数を定義します。これは、線形 ARX モデル (ARX、IDPOLY を参照) の推定時に使用した考え方と同じです。たとえば NN = [2 3 4] は、出力変数がラグ (1,2) を使用し、入力変数がラグ (4,5,6) を使用することを意味します。リグレッサー セットは {y(t-1), y(t-2), u(t-4), u(t-5), u(t-6)} となります。
通常、モデル次数は、試行錯誤により選択します。場合によっては、まず線形 ARX モデルを使用して、最適な次数に関するヒントを得ることをお勧めします。まず、関数 arxstruc
および selstruc
を使用して、線形 ARX モデルの最適な次数を求めます。
V = arxstruc(z1,z2,struc(1:5, 1:5,1:5)); % select best order by Akaike's information criterion (AIC) nn = selstruc(V,'aic')
nn = 5 1 3
AIC 条件では、nn = [na nb nk] = [5 1 3] が選択されていますが、これは、選択した ARX モデル構造で、y(t) が 6 つのリグレッサー y(t-1)、y(t-1)、...y(t-5) および u(t-3) で予測されることを意味します。これらの値を非線形モデルの推定時に使用します。
非線形 ARX (IDNLARX) モデル
IDNLARX モデルでは、y(t) = Fcn(y(t-1),y(t-1),...,y(t-5), u(t-3)) のように、モデル出力がリグレッサーの非線形関数となります。
IDNLARX モデルは通常、以下の構文を使用して推定します。
M = NLARX(Data, Orders, Fcn).
これは、線形 ARX コマンドと同じですが、追加の引数 "Fcn" で非線形関数のタイプを指定します。非線形 ARX モデルは 2 段階変換を使用して、出力を生成します。1. 学習データ (入出力信号) をリグレッサーのセットにマッピングします。数値的には、リグレッサー セットは double の行列 R で、データの 1 列が各リグレッサーに対応します。2. リグレッサーを単一出力 y = Fcn® にマッピングします。ここで Fcn() はスカラー非線形関数です。
IDNLARX モデルを推定するため、モデル次数を選択した後、使用する非線形マッピング関数 Fcn() のタイプを選択する必要があります。ウェーブレット ネットワーク関数 (idWaveletNetwork) を最初に試します。
mw1 = nlarx(z1,[5 1 3], idWaveletNetwork);
ウェーブレット ネットワークは、ウェーブレット ユニットを組み合わせて使用し、リグレッサーをモデルの出力にマッピングします。モデル mw1
は、@idnalrx オブジェクトです。その OutputFcn
プロパティに非線形マッピング関数 (ここでは idWaveletNetwork オブジェクト) を格納します。使用するこれらのユニットの数は事前に指定できますが、推定アルゴリズムに任せて最適値を自動的に特定できます。プロパティ NonlinearFcn.NumberOfUnits は、モデルに選択したユニットの数を格納します。
NLFcn = mw1.OutputFcn; NLFcn.NonlinearFcn.NumberOfUnits
ans = 3
モデル mw1
の応答とデータセット z1
の測定出力を比較して、結果を確認します。
compare(z1,mw1);
モデル プロパティの使用
最初のモデルは十分ではないので、このモデルのプロパティの一部を変更します。推定アルゴリズムで関数 idWaveletNetwork のユニット数を自動的に選択させるのではなく、この数値は明示的に指定できます。
mw2 = nlarx(z1,[5 1 3], idWaveletNetwork(8));
既定の InputName および OutputName を使用しています。
mw2.InputName mw2.OutputName
ans = 1×1 cell array {'u1'} ans = 1×1 cell array {'y1'}
IDNLARX モデルのリグレッサーは、GETREG コマンドを使用して表示できます。リグレッサーは、mw2.InputName
、mw2.OutputName
および時間遅延から作成されます。
getreg(mw2)
ans = 6×1 cell array {'y1(t-1)'} {'y1(t-2)'} {'y1(t-3)'} {'y1(t-4)'} {'y1(t-5)'} {'u1(t-3)'}
次数行列は、以下のように、仕様オブジェクトを使用して作成された線形リグレッサー セットと等価である点に注意してください。
R = linearRegressor([z1.OutputName; z1.InputName],{1:5, 3}); mw2_a = nlarx(z1, R, idWaveletNetwork(8)); getreg(mw2_a)
ans = 6×1 cell array {'y1(t-1)'} {'y1(t-2)'} {'y1(t-3)'} {'y1(t-4)'} {'y1(t-5)'} {'u1(t-3)'}
mw2_a
は、基本的に mw2
と同じモデルです。仕様オブジェクトを使用すると柔軟性が向上します。使用するのは、連続ラグを使用しない場合、または出力変数に 1 より大きい最小ラグがない場合です。
マッピング関数へのリグレッサーの割り当て
IDNLARX モデルの出力は、リグレッサーの静的関数となります。通常、このマッピング関数は、線形ブロックと非線形ブロックに加え、バイアス項 (出力のオフセット) を含みます。モデル出力は、このブロックとバイアスの出力の合計となります。モデルの複雑さを緩和するために、リグレッサーのサブセットのみを選択して、非線形ブロックおよび線形ブロックを入力できます。マッピング関数の線形および非線形コンポーネントへのリグレッサーの割り当ては、モデルの RegressorUsage プロパティを使用して実行されます。
RegUseTable = mw2.RegressorUsage
RegUseTable = 6×2 table y1:LinearFcn y1:NonlinearFcn ____________ _______________ y1(t-1) true true y1(t-2) true true y1(t-3) true true y1(t-4) true true y1(t-5) true true u1(t-3) true true
RegUseTable は、どのリグレッサーがウェーブレット ネットワークの線形および非線形コンポーネントへの入力として使用されるかを示す table です。各行は 1 つのリグレッサーに対応します。非線形コンポーネントに入力変数 'u1' の関数であるリグレッサーのみを使用するとします。以下のように、この構成を達成できます。
RegNames = RegUseTable.Properties.RowNames; I = contains(RegNames,'u1'); RegUseTable.("y1:NonlinearFcn")(~I) = false;
例として、mw2
が使用する線形リグレッサーの他に、2 次の多項式リグレッサーを使用するモデルを考えます。最初に、多項式リグレッサー仕様のセットを生成します。
P = polynomialRegressor({'y1','u1'},{1:2, 0:3},2)
P = <strong>Order 2 regressors in variables y1, u1</strong> Order: 2 Variables: {'y1' 'u1'} Lags: {[1 2] [0 1 2 3]} UseAbsolute: [0 0] AllowVariableMix: 0 AllowLagMix: 0 TimeVariable: 't'
次に、仕様 P
をモデルの Regressors プロパティに追加します。
mw3 = mw2; % start with the mw2 structure for new estimation mw3.Regressors = [mw3.Regressors; P]; % add the polynomial set to the model structure
最後にモデルのパラメーターを更新し、1 ステップ先の予測誤差を最小化することでデータを適合させます。
mw3 = nlarx(z1, mw3)
mw3 = <strong>Nonlinear ARX model with 1 output and 1 input</strong> Inputs: u1 Outputs: y1 Regressors: 1. Linear regressors in variables y1, u1 2. Order 2 regressors in variables y1, u1 Output function: Wavelet network with 8 units Sample time: 0.2 seconds Status: Estimated using NLARX on time domain data "Two tank system". Fit to estimation data: 96.76% (prediction focus) FPE: 3.499e-05, MSE: 3.338e-05
推定モデルの評価
線形および非線形のさまざまなモデルを、同じ COMPARE コマンドで一緒に比較できます。
mlin = arx(z1,[5 1 3]); % linear ARX model % compare(z1,mlin,mw2,mw3)
検証データセットで COMPARE コマンドを実行すると、モデル検証を実行できます。最初に、シミュレーション応答が z2
の測定出力信号とどのぐらいよく一致するかを確認することで、モデルを検証します。
compare(z2,mlin,mw2,mw3)
同様に、推定モデルをデータセット z3 に対して検証します。
compare(z3,mlin,mw2,mw3)
モデル残差の解析も、モデル品質を検証するために実行できます。この残差は、ホワイト (ラグ 0 以外では自身と無相関) かつ入力信号と無相関であると予想されます。残差は 1 ステップ先の予測誤差です。RESID コマンドを使用して、残差を有意でない値の境界で可視化できます。検証データセット z2
の例では、3 個の同定されたモデルの残差は以下のとおりです。
resid(z2,mlin,mw2,mw3)
legend show
残差はすべてが z2
(z2.InputData) の入力信号と無相関です。ただし、モデル mlin
と mw3
に若干の自己相関を示しています。モデル mw2
は、両方 (自己相関および相互相関) のテストで効果的に実行されています。
関数 PLOT を使用して、各 IDNLARX モデルの非線形応答を表示できます。このプロットは、本質的には非線形マッピング関数 Fcn() の特性を示しています。
plot(mw2,mw3) % plot nonlinearity response as a function of selected regressors
シグモイド ネットワークを含む非線形 ARX モデル
ウェーブレット ネットワークの代わりに、モデルで他の非線形関数も使用できます。関数 idSigmoidNetwork を試します。この関数は、シグモイド ユニットの和 (線形項およびオフセット項も追加) を使用して、リグレッサーを出力にマッピングします。ネットワークを設定して、(さまざまな値で膨張し、並進された) シグモイドの 8 ユニットを使用します。
ms1 = nlarx(z1,[5 1 3], idSigmoidNetwork(8)); compare(z1,ms1)
新しいモデルをデータセット z2 と
compare(z2,ms1);
データセット z3 で評価します。
compare(z3,ms1);
IDNLARX モデルのその他の非線形マッピング関数
idCustomNetwork は idSigmoidNetwork とよく似ていますが、シグモイド ユニット関数ではなく、独自のユニット関数を MATLAB ファイルで定義できます。この例では関数 gaussunit を使用します。
type gaussunit.m
function [f, g, a] = gaussunit(x) %GAUSSUNIT Unit function for use in idCustomNetwork (example) % %[f, g, a] = GAUSSUNIT(x) % % x: unit function variable % f: unit function value % g: df/dx % a: unit active range (g(x) is significantly non zero in the interval [-a a]) % % The unit function must be "vectorized": for a vector or matrix x, the output % arguments f and g must have the same size as x, computed element-by-element. % Copyright 2005-2021 The MathWorks, Inc. % Author(s): Qinghua Zhang f = exp(-x.*x); if nargout>1 g = - 2*x .* f; a = 0.2; end
customnet ベースのモデルの場合、関数 customnet の非線形コンポーネントへの入力として、3 番目、5 番目、および 6 番目のリグレッサーのみ使用します。
mc1_ini = idnlarx('y1','u1', [5 1 3], idCustomNetwork(@gaussunit)); Use = false(6,1); Use([3 5 6]) = true; mc1_ini.RegressorUsage{:,2} = Use; mc1 = nlarx(z1,mc1_ini);
ツリー パーティション関数は、非線形 ARX モデルに使用できるもう 1 つのマッピング関数です。関数 idTreePartition を使って、モデルを推定します。
mt1 = nlarx(z1, [5 1 3], idTreePartition);
モデル mc1
および mt1
の応答をデータ z1
と比較します。
compare(z1, mc1, mt1)
Statistics and Machine Learning Toolbox の回帰アルゴリズムの使用
Statistics and Machine Learning Toolbox™ にアクセスできる場合、ガウス過程 (GP) 関数および、バギング木またはブースティング木のアンサンブルを使用しても、非線形 ARX モデルを作成できます。GP は idGaussianProcess オブジェクトで表されるのに対して、回帰木アンサンブルは idTreeEnsemble オブジェクトで表されます。GP は既定では二乗指数カーネルを使用します。他のカーネル タイプは、オブジェクト作成時に指定できます。ここに示す例では、'Matern32' カーネルを使用します。
if exist('fitrgp','file')==2 Warn = warning('off','stats:lasso:MaxIterReached'); mp1 = nlarx(z1, [5 1 3], idGaussianProcess('Matern32')) % Create boosted ensemble of 200 regression trees En = idTreeEnsemble; En.EstimationOptions.FitMethod = 'lsboost-resampled'; En.EstimationOptions.NumLearningCycles = 200; En.EstimationOptions.ResampleFraction = 1; mp2 = nlarx(z1, [0 70 3], En) compare(z1, mp1, mp2) % compare model fits to estimation data warning(Warn) end
mp1 = <strong>Nonlinear ARX model with 1 output and 1 input</strong> Inputs: u1 Outputs: y1 Regressors: Linear regressors in variables y1, u1 Output function: Gaussian process function using a Matern32 kernel. Sample time: 0.2 seconds Status: Estimated using NLARX on time domain data "Two tank system". Fit to estimation data: 97.1% (prediction focus) FPE: 2.738e-05, MSE: 2.675e-05 mp2 = <strong>Nonlinear ARX model with 1 output and 1 input</strong> Inputs: u1 Outputs: y1 Regressors: Linear regressors in variables u1 Output function: Bagged Regression Tree Ensemble Sample time: 0.2 seconds Status: Estimated using NLARX on time domain data "Two tank system". Fit to estimation data: 91.14% MSE: 0.0002503
Deep Learning Toolbox からのネットワーク オブジェクトの使用
Deep Learning Toolbox™ が使用可能な場合、ニューラル ネットワークを使用しても、IDNLARX モデルのマッピング関数 Fcn() を定義できます。このネットワークはフィードフォワード ネットワーク (再帰コンポーネントがない) でなければなりません。
ここでは、入力数が不明の単一出力ネットワークを作成します (ネットワーク オブジェクトの零点入力サイズで指定)。入力数は、推定プロセス中にリグレッサーの数に自動的に設定されます。
if exist('feedforwardnet','file')==2 % run only if Deep Learning Toolbox is available ff = feedforwardnet([5 20]); % use a feed-forward network to map regressors to output ff.layers{2}.transferFcn = 'radbas'; ff.trainParam.epochs = 50; % Create a neural network mapping function that wraps the feed-forward network OutputFcn = idFeedforwardNetwork(ff); mn1 = nlarx(z1,[5 1 3], OutputFcn); % estimate network parameters compare(z1, mn1) % compare fit to estimation data end
Hammerstein-Wiener (IDNLHW) モデル
Hammerstein-Wiener モデルは、最大 3 つのブロックで構成されます。これらのブロックとは、非線形静的ブロック、その後に続く別の非線形静的ブロック、さらにその後に続く線形伝達関数ブロックです。最初の非線形静的ブロックがない場合、Wiener モデルと呼ばれ、2 番目の非線形静的ブロックがない場合は、Hammerstein モデルと呼ばれます。
IDNLHW モデルは通常、以下の構文を使用して推定します。
M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
ここで、Orders = [nb bf nk] は、線形伝達関数の次数と遅延を指定します。InputNonlinearity および OutputNonlinearity は 2 つの非線形ブロックの非線形関数を、線形ブロックの入出力端で指定します。M は @idnlhw モデルです。線形出力誤差 (OE) モデルは、自明な (単位ゲイン) 非線形性の場合に対応します。すなわち、IDNLHW モデルの入力および出力非線形関数が両方とも単位ゲインの場合、モデル構造は線形 OE モデルと等価です。
区分的線形非線形関数を含む Hammerstein-Wiener モデル
(区分的線形) 非線形関数 idPiecewiseLinear は、一般的な IDNLHW モデルで有効です。
Orders = [nb nf nk] において、nf は線形伝達関数の極数を指定し、IDNLARX モデルの "na" 次数とある程度関連しています。ゼロの場合、"nb" はその数と関連しています。
mhw1 = nlhw(z1, [1 5 3], idPiecewiseLinear, idPiecewiseLinear);
結果は、前述のように COMPARE で評価できます。
compare(z1,mhw1)
モデルのプロパティは、PLOT コマンドでビジュアル化できます。
ブロック線図をクリックして、入力非線形 (UNL)、線形ブロックまたは出力非線形 (YNL) のビューを選択します。
線形ブロックの場合、ステップ応答、ボード線図、インパルス応答および極-零点図内で、プロットのタイプを選択できます。
plot(mhw1)
2 つの区分的線形関数のブレーク ポイントを確認できます。これは 2 行の行列で、最初の行は入力に、2 番目の行は関数 idPiecewiseLinear の出力に対応します。
mhw1.InputNonlinearity.BreakPoints
ans = Columns 1 through 7 2.6795 3.4262 4.1729 4.9195 5.6324 6.3599 7.0874 -2.0650 -1.6316 -0.2396 1.4602 2.7554 3.6889 4.2632 Columns 8 through 10 7.8148 8.3543 9.1260 4.5074 4.4653 6.2634
飽和および不感帯非線形を含む Hammerstein-Wiener モデル
飽和関数および不感帯関数は、物理からヒントを得た非線形関数です。アクチュエータ/センサーの飽和シナリオおよび乾燥摩擦効果の記述に使用できます。
mhw2 = nlhw(z1, [1 5 3], idDeadZone, idSaturation); compare(z1,mhw2)
非線形が存在しないことは、IDUNITGAIN オブジェクトまたは空の "[]" 値で指定されます。単位ゲインは、変更のない出力への入力信号の直達です。
mhw3 = nlhw(z1, [1 5 3], idDeadZone, idUnitGain); % no output nonlinearity mhw4 = nlhw(z1, [1 5 3], [],idSaturation); % no input nonlinearity
不感帯関数および飽和関数の限界値は、推定後に確認できます。
mhw2.InputNonlinearity.ZeroInterval mhw2.OutputNonlinearity.LinearInterval
ans = -0.9982 4.7565 ans = 0.2601 0.5848
不感帯の場合は ZeroInterval の初期推定、飽和関数の場合は LinearInterval の初期推定を、モデルの推定中に指定できます。
mhw5 = nlhw(z1, [1 5 3], idDeadZone([0.1 0.2]), idSaturation([-1 1]));
Hammerstein-Wiener モデル - プロパティの指定
推定アルゴリズムのオプションは、NLHWOPTIONS コマンドを使用して指定できます。
opt = nlhwOptions();
opt.SearchMethod = 'gna';
opt.SearchOptions.MaxIterations = 7;
mhw6 = nlhw(z1, [1 5 3], idDeadZone, idSaturation, opt);
既に推定済みのモデルを、反復推定によって改善できます。
推定データ z1
で評価された mhw7
は、mhw6
よりも優れた適合を示すはずです。
opt.SearchOptions.MaxIterations = 30; mhw7 = nlhw(z1, mhw6, opt); compare(z1, mhw6, mhw7)
Hammerstein-Wiener モデル - 他の非線形関数を使用
非線形関数である区分線形 (idPiecewiseLinear)、飽和 (idSaturation)、および不感帯 (idDeadZone) は、主として IDNLHW モデルで使用するために設計されてきました。これらは、単一の変数非線形しか推定できません。より一般的な他の非線形関数であるシグモイド ネットワーク (idSigmoidNetwork)、カスタム ネットワーク (idCustomNetwork)、およびウェーブレット ネットワーク (idWaveletNetwork) も、IDNLHW モデル (IDNLARX モデルを除く) で使用できます。ただし、微分不可能な関数であるツリー パーティション (idTreePartition)、複数層ニューラル ネットワーク (idFeedforwardNetwork)、ガウス過程関数 (idGaussianProcess)、および回帰木アンサンブル (idTreeEnsemble) は、IDNLHW モデルの推定には微分可能な非線形関数が必要なため、使用できません。