LTI システムの配列を使用した非線形動作の近似

この例では、システムの非線形動作を、相互接続された LTI モデルの配列として近似する方法を説明します。

この例では、操作条件の範囲における機体のピッチ軸のダイナミクスの線形近似について説明します。結果として得られる線形システムの配列は、ダイナミクスの線形パラメーター変動 (LPV) 表現の作成に使用されます。LPV モデルは、非線形ピッチダイナミクスの近似として機能します。

線形パラメーター変動モデル

多くの場合、システムの非線形ダイナミクスは、単純な線形システムを使用して近似する必要があります。単一の線形システムを使用すると、非線形システムの操作点の小範囲の近傍に限定された動作の適切なモデルが得られます。非線形動作を操作条件の範囲で近似する必要がある場合、適切な内挿ルールで相互接続された線形モデルの配列を使用できます。そのようなモデルを LPV モデルといいます。

LPV モデルを生成するには、操作点のグリッド上で非線形モデルを平衡化および線形化します。そのためには、動作空間を少数の "スケジューリングパラメーター" でパラメーター化します。これらのパラメーターは、多くの場合、非線形システムの入力、状態および出力変数のサブセットです。LPV モデルの作成において考慮すべき重要なポイントとして、スケジューリングパラメーターセットを特定し、モデルを線形化するパラメーター値の範囲を選択することがあげられます。

このアプローチを説明するために、この例では、機体のピッチダイナミクスを近似します。

機体のピッチダイナミクス

機体のピッチ軸ダイナミクスの 3 自由度のモデルについて考えます。状態には、地球座標 $(X_e,Z_e)$ 、機体座標 $(u,w)$ 、ピッチ角度 $\theta$ 、ピッチレート $q = \dot\theta$ があります。次の図では、慣性座標系および機体座標系、飛行経路の角度 $\gamma$ 、入射角 $\alpha$ 、ピッチ角 $\theta$ の関係をまとめています。

機体のダイナミクスは非線形であり、空力は速度 $V$ に、モーメントは入射角 $\alpha$ に依存します。scdairframeTRIM モデルは、これらのダイナミクスを示しています。

open_system('scdairframeTRIM')

飛行エンベロープ内のバッチ線形化

速度 $V$ および入射角 $\alpha$ をスケジューリングパラメーターに使用します。つまり、 $\alpha$ 値と $V$ 値のグリッド上で機体モデルを平衡化します。これらの値は、scdairframeTRIM モデルの 5 つの出力の中の 2 つです。

入射角 $\alpha$ は -20 ～ 20 度の間で変動し、速度 $V$ は 700 ～ 1400 m/s の間で変動すると仮定します。スケジューリングに、線形に等間隔な $(\alpha, V)$ ペアの 15 行 12 列のグリッドを使用します。

nA = 15;   % number of alpha values
nV = 12;   % number of V values
alphaRange = linspace(-20,20,nA)*pi/180;
VRange = linspace(700,1400,nV);
[alpha,V] = ndgrid(alphaRange,VRange);

各飛行条件 $(\alpha,V)$ に対して、平衡点における機体のダイナミクスを線形化します (法線加速度とピッチのモーメントが 0)。そのためには、 $w$ および $q$ が定常値となる昇降舵の偏向角 $\delta$ およびピッチレート $q$ を計算する必要があります。

各飛行条件に対して次を行います。

関数 operspec を使用して平衡化条件を指定する。
関数 findop を使用して $\delta$ およびの平衡値を計算する。
関数 linearize を使用して、結果の操作点で機体のダイナミクスを線形化する。

機体座標 $(u,w)$ は、平衡化のための既知の状態です。したがって、これらに適切な値を指定しなければなりません。値は明示的に指定できます。しかし、この例では、モデルが $(\alpha,V)$ の各ペアに基づいてこれらの既知の値を求めます。各飛行条件 $(\alpha,V)$ につき、モデル内の値を更新して操作点の仕様を作成します。以上の手順を、180 の飛行条件すべてについて繰り返します。

clear op report
for ct = 1:nA*nV
   alpha_ini = alpha(ct);      % Incidence [rad]
   v_ini = V(ct);              % Speed [m/s]

   % Specify trim condition
   opspec(ct) = operspec('scdairframeTRIM');

   % Xe,Ze: known, not steady.
   opspec(ct).States(1).Known = [1;1];
   opspec(ct).States(1).SteadyState = [0;0];

   % u,w: known, w steady
   opspec(ct).States(3).Known = [1 1];
   opspec(ct).States(3).SteadyState = [0 1];

   % theta: known, not steady
   opspec(ct).States(2).Known = 1;
   opspec(ct).States(2).SteadyState = 0;

   % q: unknown, steady
   opspec(ct).States(4).Known = 0;
   opspec(ct).States(4).SteadyState = 1;

end
opspec = reshape(opspec,[nA nV]);

すべての操作点仕様でモデルを平衡化します。

opt = findopOptions('DisplayReport','off', 'OptimizerType','lsqnonlin');
opt.OptimizationOptions.Algorithm = 'trust-region-reflective';
[op,report] = findop('scdairframeTRIM',opspec,opt);

配列 op には、線形化に使用される findop で検出された操作点が含まれます。配列 report には、それぞれの点における入力値、出力値、および状態値の記録が含まれます。

モデルを線形化するには、最初に線形化の入力ポイントと出力ポイントを指定します。

io = [linio('scdairframeTRIM/delta',1,'in');            % delta
      linio('scdairframeTRIM/Airframe Model',1,'out');  % alpha
      linio('scdairframeTRIM/Airframe Model',2,'out');  % V
      linio('scdairframeTRIM/Airframe Model',3,'out');  % q
      linio('scdairframeTRIM/Airframe Model',4,'out');  % az
      linio('scdairframeTRIM/Airframe Model',5,'out')]; % gamma

それぞれの平衡化条件でモデルを線形化します。線形化のオフセット情報を構造体 info に保存します。

linOpt = linearizeOptions('StoreOffsets',true);
[G,~,info] = linearize('scdairframeTRIM',op,io,linOpt);
G = reshape(G,[nA nV]);
G.u = 'delta';
G.y = {'alpha','V','q','az','gamma'};
G.SamplingGrid = struct('alpha',alpha,'V',V);

G は、180 の飛行条件 $(\alpha,V)$ における線形化されたプラントモデルの 15 行 12 列の配列です。プラントのダイナミクスは、ローカルダイナミクスが不安定となるスケジューリングの位置を含め、飛行エンベロープ内で大きく変動します。

bodemag(G(3:5,:,:,:))
title('Variations in airframe dynamics')

LPV System ブロック

LPV System ブロックを使用すると、線形パラメーター変動システムのシミュレーションが容易になります。このブロックには、バッチ線形化を使用して生成した状態空間システム配列 G が必要です。また、構造体 info からの入力、出力、状態、および状態微分オフセットを使用してこの情報を拡張します。

オフセット情報を抽出します。

offsets = getOffsetsForLPV(info);
xOffset = offsets.x;
yOffset = offsets.y;
uOffset = offsets.u;
dxOffset = offsets.dx;

LPV モデルのシミュレーション

線形システム配列 G および対応するオフセットを使用する LPV System ブロックが含まれる scdairframeLPV モデル。

このモデルでは、目的の機体の軌跡に基づく入力信号を使用します。この信号 u および対応する時間ベクトル t は、モデルによって読み込まれる scdairframeLPVsimdata.mat ファイルに保存されています。シミュレーションの初期条件を指定します。

alpha_ini = 0;
v_ini = 700;
x0 = [0; 700; 0; 0];

モデルを開いてシミュレートします。

open_system('scdairframeLPV')
sim('scdairframeLPV')

このシミュレーションでは、LPV システムによる機体の応答の最適なエミュレーションが示されます。このシミュレーションにスケジューリング空間の細かいグリッドを指定することで、大量 (180) の線形モデルが導出されました。配列のサイズが大きくなると、実装コストが増加します。一方で、LPV 表現の利点は、以下に基づいてスケジューリンググリッドと、配列の線形システムの個数を調整できることです。

予測された軌跡の範囲内にあるスケジューリング部分空間
アプリケーション内で望まれる精度レベル

前の情報は、スケジューリング変数の範囲の削減に役立ちます。後の情報は、スケジューリング空間のサンプルの最適な分解能 (間隔) の選択に役立ちます。

グリッド化されたスケジューリング空間の背景に対して、前のシミュレーションのスケジューリング変数の実際の軌跡をプロットします。 $(\alpha,V)$ 出力はそれに対応するスコープを使用してログに記録されています (scdairframeLPV の Compare Responses ブロック内に含まれる)。

Stable = false(nA,nV);
for ct = 1:nA*nV
   Stable(ct ) = isstable(G(:,:,ct));
end
alpha_trajectory = Alpha_V_Data.signals(1).values(:,1);
V_trajectory = Alpha_V_Data.signals(2).values(:,1);

plot(alpha(Stable)*180/pi,V(Stable),'g.',...
   alpha(~Stable)*180/pi,V(~Stable),'k.',...
   alpha_trajectory,V_trajectory,'r.')
title('Trajectory of scheduling variables')
xlabel('\alpha')
ylabel('V')
legend('Stable locations','Unstable locations','Actual trajectory')

シミュレーション中に追跡された軌跡は、赤で表示されます。これにより、スケジューリング空間の安定および不安定な領域を調べます。

最小限のメモリを使用して、前のシミュレーションで使用されたものと同様の入力プロファイルのターゲットハードウェア上でこのモデルを実装する必要があるとします。このシミュレーションからは、速度の軌跡が主に 890 ～ 1200 m/s、入射角が主に -15 ～ 12 度の範囲にあることがわかります。この操作領域に対応するスケジューリング空間内のインデックスを求めます。

I1 = find(alphaRange>=-15*pi/180 & alphaRange<=12*pi/180);
I2 = find(VRange>=890 & VRange<=1200);

飛行条件の数を減らすために、サンプリング点の間隔の拡大もできます。たとえば、 $V$ 次元に沿って 2 つおきのサンプルと、 $\alpha$ 次元に沿って 1 つおきのサンプルのインデックスを抽出します。

I1 = I1(1:2:end);
I2 = I2(1:3:end);

LTI システム配列のサブセットを抽出します。

Gr = G(:,:,I1,I2);
size(Gr)

5x2 array of state-space models.
Each model has 5 outputs, 1 inputs, and 4 states.

新しいサンプリンググリッド Gr は、5 行 2 列のよりコンパクトなサイズになります。低次元化されたモデルのシミュレーションを実行し、元の動作を再現する際の忠実度を確認します。

低次元化されたモデルおよび対応するオフセットを使用するように LPV System ブロックを構成します。

lpvblk = 'scdairframeLPV/LPV System';
set_param(lpvblk,...
   'sys','Gr',...
   'uOffset','uOffset(:,:,I1,I2)',...
   'yOffset','yOffset(:,:,I1,I2)',...
   'xOffset','xOffset(:,:,I1,I2)',...
   'dxOffset','dxOffset(:,:,I1,I2)')
sim('scdairframeLPV')

元のモデルの応答とその LPV プロキシとの間のオーバーラップで大きな削減はありません。

LPV モデルは、高速なシミュレーションが必要な状況で、元のシステムの代わりとして機能できます。LPV モデルで使用された線形システムは、システム同定の手法を使用して取得することもできます (配列内で状態の一貫性を維持するには、一層の注意が必要です)。LPV モデルは、Simulink® Design Optimization™ の問題を初期化し、高速なハードウェアインザループのシミュレーションを実行するための代理として最適です。

bdclose('scdairframeLPV')
bdclose('scdairframeTRIM')

参考

LPV System | linearize