はじめに
皆さんは飛行機に乗るのは好きですか？
筆者は飛行機自体は好きですが、いつも乗るときにこれが人生で最後かもしれないと少し思っています笑
ライト兄弟が有人動力飛行を初めて成功させたのが1903年、そこから航空機は100年余りの間で飛躍的に進化をしています。(wikiより)
特に大きな進展の一つは操縦性を向上させる自動操縦技術の導入です。
自動操縦によってパイロットのワークロードが減った事により、航空機による輸送は大幅に増えました。
そんな自動操縦、最近ではタキシング（空港での地上走行）や離陸、着陸も全部自動で行うという研究がエアバスで行われています。
↑の動画を見てわかるように、カメラからの動画をセマンティックセグメンテーションによって滑走路の端やセンターラインとして各々検出し、逸脱しないように進路や姿勢が制御されています。
ただテストパイロットはちょっと不安気ですね笑
国内においても機体や装備品メーカー、研究機関において信頼性の高い飛行制御システムの研究開発が行われています。
今回記事で紹介するのはそのうちの一つ、突風応答軽減制御（Gust Alleviation : GA）です。
ところで皆さんは飛行機に乗っているときに機体が激しく揺れた経験がありますか？
筆者は何回かあり、簡単には落ちないと分かってはいながらも手に汗を握り、パイロットに祈りを捧げていました笑
機体が巡航中に遭遇する唐突な気流の乱れ、いわゆる乱気流による航空機事故が問題となっています。
乱気流も様々ですが、厄介なものの一つは晴天の中で突如発生する晴天乱気流（Clear-Air Turbulence : CAT）と呼ばれるものです。
CATはパイロットによる目視や気象レーダーでの確認が困難であるため、回避することが難しく、突発的な加速度外乱(G)を機体にもたらすので、対処が非常に困難です。
このような突風外乱をセンサーによって計測して、前もって対処しようということで研究がされているのがGAです。
以下はJAXAにて公開されている技術紹介です。
原理は機体に装備したドップラーライダーと呼ばれるレーザーセンサーを用いて前方の気流の乱れを光学的に検知するという仕組みです。
このセンサー技術を用いると機体がこれから遭遇するであろう突風の事前情報を入手できるため、それに備えて制御を行うことが出来ます。
制御器として研究されているのはフィードフォワードや最適制御などを用いた予見制御・予測制御の技術です。$^{[1][2][3]}$
今回はこの突風応答軽減制御についてシミュレーションで簡単な検証を行ってみたいと思います。
シミュレーションモデル
まずプラントとして機体の運動モデルが必要です。
今回はMATLABの製品ファミリーであるAerospace Blocksetで提供されている3DOF（Degree of Freedomの意）のブロックを利用します。
3DOFは下図の通り機体の縦面内の運動を定義します。
3DOFブロックでは機体に掛かる並進方向の力（N）とピッチングのモーメント（Nm）を入力として受け入れます。
これらの入力には機体に発生する揚力や抗力などの空気力に加え、エンジン等の推進力も加わります。
今回は機体としてB747というジャンボジェットのパラメーターを使用してモデル化をします。
トリムの計算と線形近似
作成したB747の機体モデル（下図）に対して、制御器の設計を行っていきます。
その事前準備として、まず機体のトリムを計算します。
ここで、トリムとは機体が操舵によって運動のつり合いを保っている状態のことです。
（制御工学で言えば、平衡状態のことです）
航空機には機体を制御するため制御舵面、エルロン$\delta_a$、エレベーター$\delta_e$、ラダー$\delta_r$があり、各舵面の角度（舵角）を変えることで発生する空気力やモーメントの向きならびにその大小を変えて機体を制御しています。
参考URL：東京航空計器株式会社
飛行中にパイロットは操縦桿を介してこれらの舵面を動かしていますが、ずっと操縦しているということになるとパイロットに与える負担が大きくなります。
特に離陸して上昇した後の巡行中は、基本一定姿勢で定常飛行を行うので、トリムがとれるように各舵角を一定に保っておく必要があります。
トリムの操作量はこのつり合い状態を保つことが出来る各舵角やエンジン推力のことを指し、自動操縦中はパイロットが操舵せずともトリムが保たれるように制御されます。
それでは、トリムの計算を行っていきます。
Simulinkモデルに対し制御系設計を支援するツールであるSimulink Control Designでは、モデルの平衡状態の計算ならびに線形近似解析をGUIで行う機能モデル線形化器が提供されています。
なお、線形解析機能については過去の記事で紹介をしていますので詳細はコチラ↓。
Simulinkモデルのアプリタブ内からモデル線形化器を起動すると以下のような画面が出てきます。
画面の線形解析タブ内にある操作点のプルダウンメニューを展開し、モデルの平衡化を選択します。
そうすると定常状態マネージャーというアプリの画面が出てきます。
この画面上でモデルから抽出した状態量（積分器など）や入力端子、出力端子に関する情報に対し、各種平衡状態に保つ際の制約を指定することかできます。
今回、機体は高度約6000mを水平定常飛行している状態で一定としたいので、各状態量や入出力に対し定常状態であるかの有無、ならびにその際の値が既知であるか、またその数値はいくつか等を制約として設定していきます。
具体的には、トリムではすべての状態が変化せず一定であるため、定常状態にチェックを入れます。
また高度は初期値の6096mという値で既知にチェック、ピッチ角$\theta$は0[rad]なので0で既知にチェック...という具合に分かっているものについて初期値の指定と計算中に既知として固定化するか、また最大最小といった範囲制約をそれぞれ指定していきます。
そして、最後に「平衡化の開始ボタン」を押すと最適化手法による収束計算が始まり、指定した各プロパティを制約条件として解いていき、解が見つかると結果を返してくれます。
なお、解が見つからない場合はエラーとしてレポートが返ってくるので、設定に誤りがないか（制約が矛盾したり非現実的でないか）を確認し、また少し条件を緩和する等の処置を行ってみて、再度解が収束するかを試みます。
今回は無事に計算が収束し、op_trim1という変数が返って来てるのでクリックして開くと以下のような結果となりました。
入力の1つ目がエレベーター舵角であり、degree換算で0.0184[deg]という値でした。
続いて2つ目がエンジン推力のレベル(0-1[p.u])であり、0.41807ということで、およそエンジンパワーは40%くらい出力するという結果です。
3つ目が機体鉛直方向の突風外乱であり、トリムの解析時では0[m/s]に指定しています。
実際にこれらの入力値を機体のモデルに入力し、トリム状態が実現できているかをシミュレーションで確認してみます。
上図は上から高度、機体鉛直方向の加速度、ピッチ角、ピッチ角速度の2000sのトレンドです。
結果から分かるようにすべての状態量がほぼ一定に保たれていることが分かります。
トリムの計算ができたので、この結果を用いてモデルの線形近似モデルを算出します。
手順としては、モデル線形化器の線形解析タブ内の線形化メニューから表示したい出力を選択します。
今回はボード線図を指定します。
この時、操作点が先ほど求めたop_trim1になっていることに注目してください。
トリム点（平衡点）で非線形モデルを線形近似します。
上図は線形化によるボード線図の結果であり、各入力から出力までのゲインと位相を表示しています。
また線形システムの値はlinsys1という名前の変数に返されており、これは状態空間（A,B,C,D）を格納しています。
そのため、線形化の結果を使って安定性の解析やこの後で行うコントローラーの設計に活用することが出来ます。
GA制御器の設計
今回は突風の事前情報を用いた予測制御を行う目的でモデル予測制御(MPC)を使っていきたいと思います。
MPCの設計にはModel Predictive Control Toolboxを使っていきます。
なお、MPCについては過去の記事↓でも紹介しています。
前の章で取得した線形近似モデルは下図の通りです。
状態方程式は縦面内の状態量に加えて、エレベーターアクチュエーターのダイナミクス（2次遅れ）も加味したモデルとなっています。
このモデルについて、線形MPCを用いて制御器を構成していきます。
MPCでは以下の標準的な2次形式を元にリアルタイムで凸最適化問題を解いていきます。
上式において$y$が出力、$r$が出力の目標、$u$が制御入力のベクトルです。
MPCは現在のステップから予測ホライズン$p$までの応答を予測し、その結果として得られる制御の偏差及び制御入力の総和を小さくするように最適化問題を解いていきます。
今回MPCの計算周期$T_s$は0.1sとし、トリム状態における各出力を目標値として、いかなる突風が発生しても機体がトリム点に留まる、またはズレても素早く復帰することを目的とします。
この際、制御入力であるエレベーター舵角及びその変化率、エンジンのスロットル出力に各々物理的な制限として上下制約を設けます。
%% Design GA controller
% Output order : theta q U W az h
% State order : h theta U W q Elevator
% Input order : de dt wg
sys = linsys1;
sys = setmpcsignals(sys,'MD',3);
Ts = 0.1;
mpcObj = mpc(sys,Ts);
mpcObj.PredictionHorizon = 60;
mpcObj.ControlHorizon = 3;
% Set trim condition
mpcObj.Model.Nominal.U = [de0 dt0 wg0];
mpcObj.Model.Nominal.Y = [theta0 q0 U0 W0 az0 h0];
% Elevator deflection limitation
mpcObj.ManipulatedVariables(1).Max = 20*pi/180;
mpcObj.ManipulatedVariables(1).Min = -20*pi/180;
mpcObj.ManipulatedVariables(1).RateMax = 100*pi/180*Ts;
mpcObj.ManipulatedVariables(1).RateMin = -100*pi/180*Ts;
% Throttle limiation
mpcObj.ManipulatedVariables(2).Max = 1;
mpcObj.ManipulatedVariables(2).Min = 0;
% Scale factor for outputs
mpcObj.OutputVariables(1).ScaleFactor = max(1,abs(theta0));
mpcObj.OutputVariables(2).ScaleFactor = max(1,abs(q0));
mpcObj.OutputVariables(3).ScaleFactor = max(1,abs(U0));
mpcObj.OutputVariables(4).ScaleFactor = max(1,abs(W0));
% Scale factor for inputs
mpcObj.ManipulatedVariables(1).ScaleFactor = 20*pi/180;
% Weights
mpcObj.Weights.OutputVariables = ones(1,6);
weight = 100*[500 500 500 500 5000 5000];
review(mpcObj)
上のコードで設計したmpcオブジェクトmpcObjを下図のSimulinkブロックにて実装します。
MPCブロックの md という入力端子が外乱情報を入力できる端子です。
今回はドップラーライダーによって予測ホライズン先までの時系列情報が取得できるとします。
具体的には、予測ホライズンを60ステップと設定しているので 60*0.1=6s 先までの突風の事前情報が毎制御周期で取得できるようにします。
なお、MPCブロックの ref 及び md ポートは行列も入力できるため、今回のように将来の情報が既知として扱える場合、それを使って先読み（$Look Ahead$）制御を行うことが出来ます。
上図において、ドップラーライダーのブロックはMATLAB Functionブロックにて実装しています。
（↓のサンプルモデル内で使用されているカスタムブロックを流用）
ここで、モデルの全体感は以下の通りです。
突風外乱のシミュレーション
それでは、実際にシミュレーションでGA制御の効果を確認してみたいと思います。
今回考慮するのは、機体が高度6096mで水平定常飛行している際に、下図に示すような機体鉛直方向の突風外乱([m/s])に遭遇するというシナリオです。（突風外乱のトレンドは文献[1]を参考に作成しました）
また性能比較のために、PID及びLQRでそれぞれ飛行制御系を構築し、同様にシミュレーションで評価してみようと思います。
なお、詳細は割愛しますが、PID及びLQRでは突風の先読み情報はありません。
上図は高度に関する結果です。
結果を見ると分かるようにMPCで制御した結果が一番優れており、次いでLQR、PIDという順になっています。
詳細には、MPCによる高度のブレは最大1m程度なのに対し、LQRでは約10m、PIDでは約30mとなっています。
続いて、上図は上からピッチ角、機体の鉛直方向加速度、エレベーター舵角（指令）、突風外乱のトレンドです。
ここで、結果を見る際の留意事項として、
・$\theta$は＋が機首上げ、－が機首下げ
・$az$は＋が下降の加速度、－が上昇の加速度
・$\delta_e$は＋が機首下げ操作、－が機首上げ操作
・$w_g$は＋が機体の上から下側に向かって吹き下ろす風、－が機体下から上側に吹き上げる風
となっています。
結果を見るとMPCは外乱の先読みができるので、将来の風に対し、早期に当て舵を取って対処していることが分かります。
例として、30秒～40秒のトレンドにおいて、吹き上げの大きな突風が発生しますが、この時、MPCでは予め＋の舵角を取っています。
これは機首下げ操作を行っていることを意味し、下から突き上げてくる風に前もってあらがっていることを示しています。
実際この時にピッチ角も－となって機首が下がっています。
LQRも何とか頑張っていますが、MPCと比較すると修正するタイミングが遅れており、なおかつ舵角が上下限いっぱいに動いているため、修正量としては少し過敏な動きです。
ただLQRは突風外乱の事前情報がなく、またMPCと違いアクチュエーターの物理制限を知らないので妥当な動きと言えます。
PIDに至っては、その制御器構造からあまりハイゲインに設定することができませんでした。（安定余裕がない）
そのため、外乱に対する感度が弱く、修正が一番遅い結果となっています。
また加速度の結果を見てもMPCが最も揺れを抑え込んでいることが分かります。
最後に
トレンドで見た結果を3Dの飛行アニメーションとして描画して終わりたいと思います。
Aerospace BlocksetではEpic Games社のUnreal EngineとのCoSim連携機能を提供しており、セスナ機や旅客機などの3Dモデルを用いて航空機の飛行状態を可視化することが出来ます。
※以下のgifアニメーションはすべて5倍速で表示しています。
PIDのアニメーション
LQRのアニメーション
MPCのアニメーション
アニメーションを見るとMPCが格段に良いのが分かります。
これなら安心して乗れそうな気がします。
参考文献
[1] 林 千瑛, 張替 正敏著：最適予見制御による乱気流中の航空機の揺れの制御、日本航空宇宙学会論文集, Vol.58, No.677, pp.164-170, 2010.
[2] 佐藤 昌之, 横山 信宏著：突風の事前情報を用いたGust Alleviation制御器による初期加速度抑制, 日本航空宇宙学会論文集, Vol.56, No.655, pp.353-362, 2008.
[3] 中山 空星, 松田 里香著：航空機の突風応答軽減制御シミュレーションSW開発, MSS技報 Vol.25.