航空機の非線形ピッチ ダイナミクスのゲイン スケジュール コントローラーの設計と検証

この例では、線形パラメーター変動状態空間モデルを使用して非線形システムのゲイン スケジュール コントローラーを設計して検証する方法を示します。この例では、LTI システムの配列を使用した非線形動作の近似 (Simulink Control Design)の例を基に、コマンド ラインで非線形動作を近似する方法を示します。

この例では、次を行います。

1. 平衡点のグリッドで非線形の Simulink® モデルを線形化します。

2. LTI モデルと平衡オフセットのグリッド配列から LPV モデルを作成します。

3. 平衡点のグリッドでコントローラーを設計します。このグリッドは、線形化に使用されるグリッドとは異なる場合があります。

4. コマンド ラインの LTI シミュレーション (固定の操作点) と LPV シミュレーション (パラメーターの軌跡を指定) を使用してコントローラーをテストします。

5. Simulink の非線形プラントでコントローラーをテストします。

航空機モデル

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

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

open_system("scdairframeOpenLoop")

操作条件の範囲

入射角 $\alpha$ は –20 ～ 20 度の間で変動し、速度 $\mathit{V}$ は 700 ～ 1400 m/s の間で変動すると仮定します。スケジューリングに、線形に等間隔な ($\alpha$,$\mathit{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$,$\mathit{V}$) に対して、平衡点における機体のダイナミクスを線形化します (法線加速度とピッチのモーメントが 0)。そのためには、$\mathit{w}$ および $\mathit{q}$ が定常値となる昇降舵の偏向角 $\delta$ およびピッチ レート $\mathit{q}$ を計算する必要があります。

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

   % Specify trim condition
   opspec(ct) = operspec("scdairframeOpenLoop");

   % 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("scdairframeOpenLoop",opspec,opt);

バッチ線形化

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

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

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

linOpt = linearizeOptions("StoreOffsets",true);
[G,~,info] = linearize("scdairframeOpenLoop",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$,$\mathit{V}$) における線形化されたプラント モデルの 15 行 12 列の配列です。

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

ssInterpolant を使用して、これらの線形化されたモデルをグリッド点間に内挿する LPV モデルを作成します。

offsets = info.Offsets;
Glpv = ssInterpolant(G,offsets);

制御設計: ピッチ軸の安定増大

制御設計の ($\alpha$,$\mathit{V}$) グリッドは、プラント モデルのグリッドに比べて粗く、範囲が小さくなっています。

nAcd = 5;   % number of alpha values
nVcd = 3;   % number of V values
alphacd = linspace(-10,10,nAcd)*pi/180;
Vcd = linspace(900,1200,nVcd);
[alphacd,Vcd] = ndgrid(alphacd,Vcd);
Ncd = nAcd*nVcd;

ピッチ レート ループには 2-DOF のゲイン スケジュール アーキテクチャを使用します。これにより、誤差信号 $\mathit{e}={\mathit{q}}_{\mathrm{ref}}-\mathit{q}$ の I のみの項が $\mathit{q}$ の比例項に結合されます。

$$\delta =\frac{{\mathit{K}}_{\mathit{i}}}{\mathit{s}}\left({\mathit{q}}_{\mathrm{ref}}-\mathit{q}\right)+{\mathit{K}}_{\mathit{p}}\mathit{q}$$

調整可能なゲイン ${\mathit{K}}_{\mathit{p}}$ および ${\mathit{K}}_{\mathit{i}}$ は、($\alpha$,$\mathit{V}$) に対する線形依存性があるゲイン曲面として定義されます。

s = tf("s");
Domain = struct("alpha",alphacd,"V",Vcd);
ShapeFcn = @(x,y) [x y];
Ki = tunableSurface("Ki",-1,Domain,ShapeFcn);
Kp = tunableSurface("Kp",1,Domain,ShapeFcn);
K = [Ki/s Kp] * [1 -1;0 1];
K.InputName = {'qref';'q'};
K.OutputName = {'delta'};

機体の LPV モデルを調整グリッド (alphacd,Vcd) でサンプリングし、閉ループ モデルを作成します。

[Ga,Goffsets] = sample(Glpv,[],alphacd,Vcd);
T0 = connect(Ga,K,'qref','q','delta');

Warning: The following block outputs are not used: alpha,V,az,gamma.

ステップ応答とループ整形の要件に従ってゲイン曲面を調整します。

R1 = TuningGoal.StepTracking("qref","q",1/2);
R2 = TuningGoal.MaxLoopGain("delta",50,1);
R3 = TuningGoal.MaxGain("qref","q",tf(15,[1 0.01]));
T = systune(T0,[R1 R2],R3);

Final: Soft = 4.94, Hard = 0.19413, Iterations = 30

Warning: StepTracking goal: Feedback configuration has fixed integrators that cannot be stabilized with available tuning parameters. Make sure these are modeling artifacts rather than physical instabilities.

showTunable(T)

Ki =
   -4.4209   -0.0712    1.2825
-----------------------------------
Kp =
    1.1167   -0.0000   -0.2609

調整結果を表示します。

viewGoal([R1 R2 R3],T)

調整後のゲイン曲面は、$\alpha$ に対する依存性が弱く、$\mathit{V}$ に対する依存性が強くなっています。

Ki = setBlockValue(Ki,T);
clf
subplot(211)
viewSurf(Ki)
Kp = setBlockValue(Kp,T);
subplot(212)
viewSurf(Kp)

LPV コントローラーを作成します。コントローラーで平衡性についての段階的な補正を行うように、平衡偏向角を出力のオフセットとして含めます。

Ka = ss(setBlockValue(K,T));
Koffsets = reshape(struct("y",{Goffsets.u}),size(Goffsets));
Klpv = ssInterpolant(Ka,Koffsets);

閉ループ LPV モデルを作成します。プラントとコントローラーのサンプリングと内挿を 2 つの異なる ($\alpha$,$\mathit{V}$) グリッドで行っていることに注意してください。

Tlpv = feedback(Glpv*Klpv,1,2,3,+1);
Tlpv = Tlpv(:,1); % from qref to plant outputs

設計グリッドでの ${\mathit{q}}_{\mathrm{ref}}$ から $\mathit{q}$ までの閉ループ応答をプロットします。

clf
step(sample(Tlpv(3,:),[],alphacd,Vcd),6);
grid
title("Closed-loop response qref to q on design grid")

非線形シミュレーション

非線形と LPV のシミュレーションを比較するために、いずれかの平衡操作点から初期化されるピッチ レート ループのステップ応答をシミュレートします。この平衡化条件は密なグリッドで選択され、設計グリッドには属しません。

aidx = 9;
Vidx = 6;
alpha0 = alphaRange(aidx);
V0 = VRange(Vidx);
q0 = offsets(aidx,Vidx).x(4);

LPV コントローラーの初期状態 ${\mathit{x}}_{\mathrm{K0}}$ を選択して、選択された操作点に対する平衡偏向角 ${\delta }_0$ を提供します。LPV コントローラーの出力は次のとおりです。

$$\delta ={\delta }_0+{\mathit{C}}_{\mathit{K}}{\mathit{x}}_{\mathit{K}}+{\mathit{D}}_{\mathit{K},2}\mathit{q},$$

${\mathit{x}}_{\mathrm{K0}}$ は次を満たさなければなりません。

$${\mathit{C}}_{\mathit{K}}{\mathit{x}}_{\mathrm{K0}}+{\mathit{D}}_{\mathit{K},2}{\mathit{q}}_0=0.$$

K0 = sample(Klpv,[],alpha0,V0);
xK0 = -K0.C\K0.D(2)*q0;

(alpha0,V0) の平衡化条件から開始して、t=0 における q0 から q0+0.05 までのステップ変化を適用します。

tstep = 0;
rstep0 = q0;
rStepAmp = 0.05;
rstepf = rstep0+rStepAmp;
Tf = 5;

閉ループの Simulink モデルを開き、非線形シミュレーションを初期化します。

open_system("scdairframeClosedLoop")

alpha_ini = alpha0;
v_ini = V0;
q_ini = q0;
yK0 = reshape([Koffsets.y],[1 1 nAcd nVcd]);

非線形シミュレーションの結果を取得します。

sim("scdairframeClosedLoop",[0 tstep+Tf]);
tsim = y.Time;
ysim = y.Data;

比較のために、最初に次の操作条件で LTI の応答を計算します。

t = linspace(0,Tf,250);
y = step(sample(Tlpv,[],alpha0,V0),t);
ylin = [alpha0,V0,q0] + rStepAmp * y(:,1:3);

LPV シミュレーション

非線形の機体モデルの LPV 近似 Glpv を使用して応答を計算します。これを行うには、$\mathit{p}=\left(\alpha ,\mathit{V}\right)$ の軌跡が必要です。1 つ目の選択肢として、上記の LTI シミュレーションから得られるおおよその軌跡を使用できます ($\alpha$ と $\mathit{V}$ が Tlpv の最初の 2 つの出力であることを思い出してください)。

p_ideal = ylin(:,1:2);

プラントの初期状態は平衡化解析から取得できます。Klpv の初期状態 xK0 を加えることで全体の初期状態 xinit が得られます。

xinit = [offsets(aidx,Vidx).x ; xK0];

おおよそのパラメーターの軌跡 p_ideal を使用して LPV の応答をシミュレートします。

qref = rstepf*ones(1,numel(t));
y1 = lsim(Tlpv,qref,t,xinit,p_ideal);

実際のパラメーターの軌跡は、$\alpha$,$\mathit{V}$ が状態 $\mathit{u}$,$\mathit{w}$ に次のように依存するため内因的なものになります。

$$\tan \left(\alpha \right)=\frac{\mathit{w}}{\mathit{u}},$$

$$\mathit{V}=\sqrt{{\mathit{u}}^2+{\mathit{w}}^2}.$$

より正確な結果を得るために、この依存関係を関数 $\mathit{p}=\mathit{F}\left(\mathit{t},\mathit{x},\mathit{u}\right)$ を使用して記述できます。

pFcn = @(t,x,u) [atan2(x(3),x(2)) ; sqrt(x(2)^2+x(3)^2)];
[y2,~,~,p2] = lsim(Tlpv,qref,t,xinit,pFcn);

シミュレーション結果を比較します。

figure(1)
clf
subplot(411)
plot(t,y1(:,1),t,y2(:,1),t,ylin(:,1),tsim,ysim(:,1),"k--")
ylabel("alpha")
grid on
title("Linear, LPV, and Nonlinear Simulations")
legend("p ideal","p true","Linear","Nonlinear","location","southeast")
subplot(412)
plot(t,y1(:,2),t,y2(:,2),t,ylin(:,2),tsim,ysim(:,2),"k--")
ylabel("V")
grid on
subplot(413)
plot(t,y1(:,3),t,y2(:,3),t,ylin(:,3),tsim,ysim(:,3),"k--")
ylabel("q")
grid on
subplot(414)
plot(t,y1(:,5),t,y2(:,5),tsim,ysim(:,5),"k--")
ylabel("gamma")
grid on

LPV シミュレーションは非線形シミュレーションに非常に近く、機体の LPV モデルが有効な代理となり、ゲイン スケジュール LPV コントローラーの性能が十分であることがわかります。予想どおり、実際のパラメーターの軌跡を使用した LPV シミュレーションの方が代理で p_ideal を使用した場合よりも精度が少し高くなっています。

参考

ssInterpolant | sample | linearize (Simulink Control Design) | systune | ndgrid

関連するトピック

• LTV モデリングと LPV モデリング
• MATLAB および Simulink での LTV モデルと LPV モデルの使用
• 昇圧コンバーター モデルの LPV 近似
• 磁気浮上システムのバッチ線形化の結果からの LPV モデル
• 風力タービンの制御設計