Main Content

非線形 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 で評価します。z1z2 と似ていますが、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™ では、リグレッサー セットを明示的に作成する必要はありません。形式 (多項式の次数など)、寄与する変数、ラグを指定するだけで、大規模な集合を暗黙的に生成できます。仕様オブジェクト linearRegressorpolynomialRegressor および 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.InputNamemw2.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) の入力信号と無相関です。ただし、モデル mlinmw3 に若干の自己相関を示しています。モデル 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 モデルの推定には微分可能な非線形関数が必要なため、使用できません。