多変数システムでのデータの解析とモデルの特定

データセット

まず、データセット SteamEng を見てみましょう。

load SteamEng

このデータセットは、ラボ スケール スチーム エンジンから収集しています。制御バルブの後にスチームの入力 Pressure (圧力) (実際には、圧搾空気) と、出力軸に接続されたジェネレーター上に Magnetization voltage (磁気電圧) があります。

出力は、ジェネレーターの Generated voltage (出力電圧) と、ジェネレーターの Rotational speed (回転速度) (出力 AC 電圧の周波数) です。サンプル時間は 50 ms です。

まず、iddata オブジェクトに測定チャネルを収集します。

steam = iddata([GenVolt,Speed],[Pressure,MagVolt],0.05);
steam.InputName  = {'Pressure';'MagVolt'};
steam.OutputName = {'GenVolt';'Speed'};

データをよく見てください。

plot(steam(:,1,1))

plot(steam(:,1,2))

plot(steam(:,2,1))

plot(steam(:,2,2))

ステップ応答とインパルス応答

ダイナミクスを把握する最初のステップでは、データから直接推定された複数のチャネル間のステップ応答を確認します。

mi = impulseest(steam,50);
clf, step(mi)

信頼領域による応答

応答の有意性を評価するには、標準偏差 3 に対応する信頼領域で、インパルス プロットを代わりに使用できます。

showConfidence(impulseplot(mi),3)

非対角影響の方が明らかに大きくなっています (y スケールを比較してください)。つまり、GenVolt は主に (ダイナミクスがあまり大きくない) MagVolt の影響を受け、Speed は主に Pressure に依存しています。一見したところ、MagVolt から Speed への応答の有意性は低くなります。

2 入力/2 出力モデル

最初のクイック テストは、既定の連続時間状態空間予測誤差モデルを確認するためにも使用します。推定には、データの前半のみを使用します。

mp = ssest(steam(1:250))

mp =
  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       x2       x3       x4
   x1   -29.43   -4.561   0.5994     -5.2
   x2   0.4849  -0.8662   -4.101   -2.336
   x3    2.839    5.084   -8.566   -3.855
   x4   -12.13   0.9224    1.818   -34.29
 
  B = 
       Pressure   MagVolt
   x1    0.1033    -1.617
   x2   -0.3027  -0.09415
   x3    -1.566    0.2953
   x4  -0.04477    -2.681
 
  C = 
                 x1       x2       x3       x4
   GenVolt   -16.39   0.3767  -0.7566    2.808
   Speed     -5.623    2.246  -0.5356    3.423
 
  D = 
            Pressure   MagVolt
   GenVolt         0         0
   Speed           0         0
 
  K = 
        GenVolt     Speed
   x1   -0.3555   0.08529
   x2  -0.02311     5.195
   x3     1.526     2.132
   x4     1.787   0.03216
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 40
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Estimated using SSEST on time domain data.              
Fit to estimation data: [86.9;74.84]% (prediction focus)
FPE: 3.897e-05, MSE: 0.01414                            
 

データから直接推定したステップ応答と比較します。

h = stepplot(mi,'b',mp,'r',2); % Blue for direct estimate, red for mp
showConfidence(h)

表示されている信頼限界の範囲内で許容できる変動と適切に一致しています。

状態空間モデルの質をテストするには、推定に使用していないデータ部分でシミュレーションを行ない、出力を比較します。

compare(steam(251:450),mp)

このモデルは、検証データで Generated Voltage (出力電圧) を非常に良好に再現できており、速度でも妥当な結果が得られています (プルダウン メニューを使用すると、各出力の適合を表示できます)。

スペクトル解析

同様に、スペクトル解析推定による mp の周波数応答で比較を行います。

msp = spa(steam);

bode(msp,mp)

clf, bode(msp,'b',mp,'r')

プロットを右クリックすると、別の I/O の組み合わせを選択して詳細を確認することができます。また、[特性]、 [信頼領域] を選択して、ボード線図の信頼性を図示することもできます。

既に述べたように、MagVolt から Speed への応答は有意ではなく、推定は困難です。

単入力/単出力 (SISO) モデル

このデータセットからは、良好なモデルを素早く得ることができます。これを使用しない場合、特定のチャネルでサブモデルを試行し、有意な影響を確認する必要があります。ツールボックス オブジェクトでは、このような作業で必要なデータ記録を完全にサポートしています。入力名と出力名がここでは必須となります。

ステップ応答からは、MagVolt は主に GenVolt に影響を及ぼし、Pressure は主に Speed に影響を及ぼすことがわかります。ここで、2 つのシンプルな SISO モデルを構築します。チャネルを選択するときには、名前と番号の両方を使用できます。

m1 = tfest(steam(1:250,'Speed','Pressure'),2,1); % TF model with 2 poles 1 zero
m2 = tfest(steam(1:250,1,2),1,0) % Simple TF model with 1 pole.

m2 =
  From input "MagVolt" to output "GenVolt":
    18.57
  ---------
  s + 43.53
 
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:                                   
Estimated using TFEST on time domain data.
Fit to estimation data: 73.34%            
FPE: 0.04645, MSE: 0.04535                
 

これらのモデルを MIMO モデル mp と比較します。

compare(steam(251:450),m1,m2,mp)

SISO モデルはフル モデルとほぼ同等です。次に、ナイキスト線図を比較してみましょう。m1 は青色、m2 は緑色、mp は赤色です。並べ替えは自動的に行われます。mp はすべての入力/出力ペアを示し、m1 は Pressure から Speed のみ、m2 は MagVolt から GenVolt のみを含みます。

clf
showConfidence(nyquistplot(m1,'b',m2,'g',mp,'r'),3)

SISO モデルは、各出力を適切に再現しています。

一般的に、出力数が増えると (解析対象が多くなるため)、モデルの適合が難しくなります。反対に、入力数が増えると、容易になります。

2 入力/単出力モデル

出力 GenVolt で良好な結果を得るために、両方の入力を使用できます。

m3 = armax(steam(1:250,'GenVolt',:),'na',4,'nb',[4 4],'nc',2,'nk',[1 1]);
m4 = tfest(steam(1:250,'GenVolt',:),2,1);
compare(steam(251:450),mp,m3,m4,m2)

入力 Pressure をモデル m3 (離散時間) および m4 (連続時間) に含めることにより、入力として MagVolt のみを使用している m2 と比べて約 10% の改善が可能となります。

SISO モデルの結合

必要な場合、零点ダミー モデルを作成し、2 つの SISO モデル m1 および m2 を 1 つの "非対角" モデルにまとめることができます。

mdum = idss(zeros(2,2),zeros(2,2),zeros(2,2),zeros(2,2));
mdum.InputName = steam.InputName;
mdum.OutputName = steam.OutputName;
mdum.ts = 0; % Continuous time model
m12 = [idss(m1),mdum('Speed','MagVolt')];    % Adding Inputs. 
                                             % From both inputs to Speed
m22 = [mdum('GenVolt','Pressure'),idss(m2)]; % Adding Inputs. 
                                             % From both inputs to GenVolt

mm = [m12;m22]; % Adding the outputs to a 2-by-2 model.

compare(steam(251:450),mp,mm)

出力の説明において、"非対角" モデル mm が m1 および m2 と同じように機能することがわかります。