最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

同定による複雑なシステムの線形近似

この例では、線形モデルの同定により、複雑な非線形システムの線形近似を取得する方法を説明します。この方法は、システムを励起する入力信号の選択に基づいています。線形近似は、選択した入力信号の非線形モデルをシミュレーションした応答に線形モデルを適合することで得られます。

この例では Simulink®、Control System Toolbox™ および Simulink Control Design™ を使用します。

はじめに

多くの場合、線形モデルは、複数の局所的条件により複雑になっている非線形システムを単純化することにより得ます。たとえば、航空機のダイナミクスの忠実度の高いモデルは詳細な Simulink モデルで表すことができます。このようなシステムのシミュレーションの高速化、操作点での局所的動作の調査、補償器の設計を行うための一般的な方法として、システムの線形化が用いられます。操作点に関して元のモデルを解析的に線形化すると、通常は元のモデルの状態と同じ数または近い数の次数のモデルが得られます。解析または制御システム設計で使用すると、この次数は入力のクラスに対して必要以上に高くなる場合があります。そのため、システムのシミュレーションから入出力データを収集し、そのデータを使用して目的の次数の線形モデルを生成するための代替方式を検討することが必要になる場合もあります。

F14 モデルの解析的な線形化

F14 モデルを考えます。このモデルは既に線形モデルですが、Derivative ブロックと外乱の原因が含まれているため、出力の特性に影響が生じる場合があります。次のように、1 つの入力端子と 2 つの出力端子の間でこのモデルを "線形化" することができます。

open_system('idF14Model')
IO = getlinio('idF14Model')
syslin = linearize('idF14Model',IO)
3x1 vector of Linearization IOs: 
--------------------------
1. Linearization input perturbation located at the following signal:
- Block: idF14Model/Pilot
- Port: 1
- Signal Name: Stick Input
2. Linearization output measurement located at the following signal:
- Block: idF14Model/Gain5
- Port: 1
- Signal Name: Angle of Attack
3. Linearization output measurement located at the following signal:
- Block: idF14Model/Pilot G force
- Port: 1
- Signal Name: Pilot G force

syslin =
 
  A = 
                 Transfer Fcn    Derivative  Transfer Fcn   Derivative1
   Transfer Fcn       -0.6385             0         689.4             0
   Derivative               1        -1e+05             0             0
   Transfer Fcn      -0.00592             0       -0.6571             0
   Derivative1              0             0             1        -1e+05
   Actuator Mod             0             0         1.424             0
   Alpha-sensor      0.001451             0             0             0
   Stick Prefil             0             0             0             0
   Pitch Rate L             0             0             1             0
   Proportional             0             0       -0.8156             0
 
                 Actuator Mod  Alpha-sensor  Stick Prefil  Pitch Rate L
   Transfer Fcn         -1280             0             0             0
   Derivative               0             0             0             0
   Transfer Fcn        -137.7             0             0             0
   Derivative1              0             0             0             0
   Actuator Mod           -20         2.986        -39.32         -1.67
   Alpha-sensor             0        -2.526             0             0
   Stick Prefil             0             0        -22.52             0
   Pitch Rate L             0             0             0        -4.144
   Proportional             0         -1.71         22.52        0.9567
 
                 Proportional
   Transfer Fcn             0
   Derivative               0
   Transfer Fcn             0
   Derivative1              0
   Actuator Mod        -3.864
   Alpha-sensor             0
   Stick Prefil             0
   Pitch Rate L             0
   Proportional             0
 
  B = 
                 Stick Input
   Transfer Fcn            0
   Derivative              0
   Transfer Fcn            0
   Derivative1             0
   Actuator Mod            0
   Alpha-sensor            0
   Stick Prefil            1
   Pitch Rate L            0
   Proportional            0
 
  C = 
                 Transfer Fcn    Derivative  Transfer Fcn   Derivative1
   Angle of Att      0.001451             0             0             0
   Pilot G forc         -3106     3.106e+08     7.083e+04    -7.081e+09
 
                 Actuator Mod  Alpha-sensor  Stick Prefil  Pitch Rate L
   Angle of Att             0             0             0             0
   Pilot G forc             0             0             0             0
 
                 Proportional
   Angle of Att             0
   Pilot G forc             0
 
  D = 
                 Stick Input
   Angle of Att            0
   Pilot G forc            0
 
Continuous-time state-space model.

syslin は 2 つの出力、1 つの入力、9 つの状態があるモデルです。これは、"Stick Input" 入力から 2 つの出力への線形化パスに元のシステムの 9 つの状態が含まれているためです。これは operpoint を使用して検証できます。

operpoint('idF14Model')
 Operating point for the Model idF14Model.
 (Time-Varying Components Evaluated at time t=0)

States: 
----------
(1.) idF14Model/Actuator Model
      x: 0            
(2.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.1
      x: 0            
(3.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.2
      x: 0            
(4.) idF14Model/Controller/Alpha-sensor Low-pass Filter
      x: 0            
(5.) idF14Model/Controller/Pitch Rate Lead Filter
      x: 0            
(6.) idF14Model/Controller/Proportional plus integral compensator
      x: 0            
(7.) idF14Model/Controller/Stick Prefilter
      x: 0            
(8.) idF14Model/Dryden Wind Gust Models/Q-gust model
      x: 0            
(9.) idF14Model/Dryden Wind Gust Models/W-gust model
      x: 0            
      x: 0            

Inputs: None 
----------

この次数を低くし、一方で選択した矩形波 ("Stick input") 入力に対する応答の忠実度を維持することは可能でしょうか。

同定データの準備

モデルをシミュレーションし、0:30 秒間の入力矩形波 (u) と出力 "Angle of attack" (y1) および "Pilot G force" (y2) のログを生成します。このデータは、等間隔の時間ベクトル (サンプル時間 0.0444 秒) への内挿の後に、"idF14SimData.mat" ファイルに保存されます。

load idF14SimData
Z = iddata([y1, y2],u,Ts,'Tstart',0);
Z.InputName = 'Stick input'; Z.InputUnit = 'rad/s';
Z.OutputName = {'Angle of attack', 'Pilot G force'};
Z.OutputUnit = {'rad', 'g'};
t = Z.SamplingInstants;

subplot(311)
plot(t,Z.y(:,1)), ylabel('Angle of attack (rad)')
title('Logged Input-Output Data')

subplot(312)
plot(t,Z.y(:,2)), ylabel('Pilot G force (g)')

subplot(313)
plot(t,Z.u), ylabel('Stick input (rad/s)')
axis([0 30 -1.2 1.2])
xlabel('Time (seconds)')

状態空間モデルの推定

ssest コマンドを使用して、次数 2 ~ 4 の状態空間モデルを推定します。ここでは、推定で "シミュレーション" フォーカスを使用するよう設定し、モデルの外乱コンポーネントを推定しないよう選択します。

opt = ssestOptions('Focus','simulation');
syslin2 = ssest(Z, 2, 'DisturbanceModel', 'none', opt);
syslin3 = ssest(Z, 3, 'DisturbanceModel', 'none', opt);
syslin4 = ssest(Z, 4, 'DisturbanceModel', 'none', opt);

線形化されたモデル syslin と 3 つの同定されたモデルのパフォーマンスをデータと比較します。syslin は SS モデルで、syslin2syslin3syslin4 は IDSS モデルです。

syslin.InputName = Z.InputName;
syslin.OutputName = Z.OutputName; % reconcile names to facilitate comparison
clf
compare(Z, syslin, syslin2, syslin3, syslin4)

プロットからは、3 次モデル (syslin3) が既定 (t=0) の操作条件に関する航空機のダイナミクスの線形近似として十分に機能していることがわかります。解析的な線形化により得られたもの (syslin) と比べて、データの適合はやや優れたものとなっています。元の idF14Model が線形の場合、線形化の結果 (syslin) がデータと 100% 適合しないのはなぜでしょうか。これには 2 つの理由があります。

  1. 測定された出力は突風の影響を受けています。そのため、ログに記録された出力は Stick input の単なる関数ではありません。外乱が影響を及ぼしています。

  2. "Pilot G force" ブロックでは Derivative ブロックが使用されていて、線形化は時定数 "c" の値に依存しています。"c" の値は小さくしなければなりませんが (ここでは 1e-5 を使用します)、0 にはならないようにします。c の値が 0 以外の場合、線形化で近似誤差が発生します。

モデル syslin3 のパラメーターを見てみましょう。応答は十分に再現されているようです。

syslin3
syslin3 =
  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
   x1   -1.006   -2.029  -0.5842
   x2    8.284   -19.39    5.611
   x3   -2.784    12.63   -6.956
 
  B = 
       Stick input
   x1       0.2614
   x2        5.512
   x3       -3.606
 
  C = 
                      x1       x2       x3
   Angle of att   -8.841  -0.5347   -1.402
   Pilot G forc   -86.42   -15.85   -66.12
 
  D = 
                 Stick input
   Angle of att            0
   Pilot G forc            0
 
  K = 
       Angle of att  Pilot G forc
   x1             0             0
   x2             0             0
   x3             0             0
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 18
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                       
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.4;97.02]%         
FPE: 2.367e-05, MSE: 0.1103                   

低次元化と推定によるモデルの単純化

また、線形化されたモデル syslin を低次元化し、この低次元化されたモデルのパラメーターがデータ Z に最もよく適合するように調整するという方法もあります。低次元化のための適切な値を求めるには、hsvd を使用します。

[S, BalData] = hsvd(syslin);
clf; bar(S)

棒グラフからは、状態 4 以降では特異値が非常に小さいことがわかります。そのため、低次元化して 3 次にするのが最適と考えられます。

sysr = balred(syslin,3,BalData)
opt2 = bodeoptions; opt2.PhaseMatching = 'on';
clf; bodeplot(sysr,syslin,opt2)
sysr =
 
  A = 
             x1        x2        x3
   x1    -1.205     -5.28     45.51
   x2   -0.5497     4.546    -27.37
   x3  -0.06248     7.327    -27.93
 
  B = 
       Stick input
   x1       -328.5
   x2        -1048
   x3        682.7
 
  C = 
                         x1          x2          x3
   Angle of att  -0.0006108  -0.0005806    0.000794
   Pilot G forc   -0.007392    -0.03346      0.1282
 
  D = 
                 Stick input
   Angle of att     -0.03784
   Pilot G forc       -1.617
 
Continuous-time state-space model.

ボード線図は 10 ラジアン/秒の高い忠実度を示しています。sysrcompare プロットが示すように、状態が 9 つある元のモデルと同じ程度に応答をエミュレートできます。

compare(Z, sysr, syslin)

sysr のパラメーターを調整して、データとの適合を改善しましょう。この推定では、"Levenberg-Marquardt" 検索方式を選択し、許容できる反復の最大数を 10 に変更します。これらの選択は試行錯誤で行われます。また、推定の進捗状況の表示もオンにします。

opt.Display = 'on';
opt.SearchMethod = 'lm';
opt.SearchOptions.MaxIterations = 10;
sysr2 = ssest(Z, sysr, opt)
compare(Z, sysr2)
sysr2 =
  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
   x1   -2.227   -5.344     45.5
   x2  0.03107    3.495   -27.19
   x3   0.6561    6.546   -29.82
 
  B = 
       Stick input
   x1       -328.5
   x2        -1048
   x3        682.7
 
  C = 
                         x1          x2          x3
   Angle of att  -0.0006251  -0.0003531   -0.001759
   Pilot G forc     -0.0112    -0.02443      0.1133
 
  D = 
                 Stick input
   Angle of att     0.003051
   Pilot G forc       0.6397
 
  K = 
       Angle of att  Pilot G forc
   x1             0             0
   x2             0             0
   x3             0             0
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: yes
   Disturbance component: none
   Number of free coefficients: 20
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                       
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.78;97.03]%        
FPE: 1.434e-05, MSE: 0.1099                   

改善されたモデル sysr2 は F14 モデルの応答と非常によく一致しています (最初の出力では約 99%、次の出力では 97%)。

まとめ

ここでは、複雑なシステムの線形近似を取得するための、解析的な線形化にかわる代替方法について説明しました。結果は特定の入力信号から生成され、厳密に言えばその入力に対してしか適用できません。結果をさまざまな入力プロファイルに適用できるようにするには、各種の入力を使用して複数のシミュレーションを実行します。その後、生成されたデータセットを 1 つの複数実験データセットに結合し (iddata/merge を参照)、これを使用して推定を実行します。この例では、わかりやすいように、複雑ではありますが線形のシステムを使用しました。この方式の実際の利点は、非線形システムで見られます。

また、線形システムを低次元化し、その一方で低次元化されたモデルが元の Simulink モデルのシミュレーションされた応答に忠実になるようにする方法についても説明しました。低次元化されたモデル sysr の役割は、推定されたモデル sysr2 の初期推定を生成することです。この手法では、線形システムであれば、クラスが異なっていても、推定の初期モデルとして使用できます。

bdclose('idF14Model')

その他の情報

System Identification Toolbox を使った動的システムの同定の詳細については、System Identification Toolbox 製品の情報ページを参照してください。