特許 5771893 乱気流抑制制御方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5771893

(24)【登録日】2015年7月10日

(45)【発行日】2015年9月2日

(54)【発明の名称】乱気流抑制制御方法

(51)【国際特許分類】

   B64C 13/16 20060101AFI20150813BHJP

【ＦＩ】

   B64C13/16 A

【請求項の数】2

【全頁数】19

(21)【出願番号】特願2009-264792(P2009-264792)

(22)【出願日】2009年11月20日

(65)【公開番号】特開2011-105248(P2011-105248A)

(43)【公開日】2011年6月2日

【審査請求日】2012年10月26日

【新規性喪失の例外の表示】特許法第３０条第１項適用 日本航空宇宙学会論文集 第５７巻 第６６８号別刷 ２００９年９月号 ｐｐ．３４５−３５３（平成２１年９月５日発行）

(73)【特許権者】

【識別番号】503361400

【氏名又は名称】国立研究開発法人宇宙航空研究開発機構

(73)【特許権者】

【識別番号】504165591

【氏名又は名称】国立大学法人岩手大学

(74)【代理人】

【識別番号】100092200

【弁理士】

【氏名又は名称】大城 重信

(74)【代理人】

【識別番号】100110515

【弁理士】

【氏名又は名称】山田 益男

(72)【発明者】

【氏名】佐藤 昌之

(72)【発明者】

【氏名】横山 信宏

(72)【発明者】

【氏名】佐藤 淳

【審査官】 芦原 康裕

(56)【参考文献】

【文献】 特表２００９−５１８２２２（ＪＰ，Ａ）

【文献】 特開平０１−３０９８９４（ＪＰ，Ａ）

【文献】 特開２００３−１７５８１３（ＪＰ，Ａ）

【文献】 特開平０５−１７０１８４（ＪＰ，Ａ）

【文献】 特表２００８−５００５２５（ＪＰ，Ａ）

【文献】 特開２００２−２５８９０５（ＪＰ，Ａ）

【文献】 特開平０９−１８５４０１（ＪＰ，Ａ）

【文献】 佐藤昌之，横山信宏，佐藤淳，乱気流の事前情報を用いたロバストモデル予測制御によるGust Alleviation制御，第４０期年会講演会講演集，日本，日本航空宇宙学会，２００９年，CD-ROM，A11,p.89-98

【文献】 Hamaki Inokuchi, Hisamichi Tanaka, Toshiyuki Ando，Development of an Onboard Doppler Lidar for Flight Safety，Journal of Aircraft，米国，American Institute of Aeronautics and Astronautics，２００９年 ７月，Ｖｏｌ．４６，Ｎｏ．４，１４１１−１４１５頁

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｂ６４Ｃ １３／１６

(57)【特許請求の範囲】

【請求項1】

連続時間の航空機の線形運動方程式を離散化したシステムを制御対象とし、舵面を動かすアクチュエータや推力変化をもたらすエンジン等の動特性の離散時間システムとしてこれらの結合システムが与えられるものとし、システムのモデリング誤差に対するロバスト性（頑強性）を確保するために制御入力端に[式１]で示される不確かな無駄時間を有する遅れシステムを用い、事前乱気流情報の観測誤差に対するロバスト性を確保するために[式２]に示されるあらかじめ想定された誤差幅を用いた事前乱気流情報を用い、さらに自機の運動の状態の観測誤差に対するロバスト性を確保するために[式３]に示されるあらかじめ想定された誤差幅を用いた機体運動状態を用い、機体の運動に関する状態および入力についての制約条件下において、評価関数を最適化し、得られた最適入力時系列のうち、最初のデータを実際の入力として加え、次のステップでは、新たに得られた乱気流情報および機体運動に関する情報を用いて、再度最適化問題を解き、得られた最適入力時系列の最初のデータを実際の入力として加え、これを毎ステップ繰り返すことを特徴とする乱気流抑制制御方法。

【数1】

ただし、上式において、Ｔ_d は不確かな遅れ時間に相当する遅れステップから構成された遅れステップ数の集合を、d,d+1,‥‥は各要素の遅れ時間を、ｄアンダーバーは時間ステップ数で表現された一番小さい遅れステップの値を、ｄアッパーバーは時間ステップ数で表現された一番大きい遅れステップの値を、ｗ_k+j|k はｋステップにおいてｊステップ先に遭遇する乱気流の観測値、ｗ_k+j はｋステップにおいてｊステップ先に遭遇する乱気流の真値、Ｘ_jΔ_ｗ は乱気流の観測誤差をあらわす。

【請求項2】

機械運動に関する状態量の観測誤差を加え、上記の観測誤差および不確かさのすべての組み合わせに対して共通の入力時系列を求める定式化を行った請求項１に記載の乱気流抑制制御方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、自機の運動の状態を取得できる航空機が、自機の前にある乱気流の風速・風向データをも取得できる場合に、乱気流による垂直加速度を抑制する技術に関する。

【背景技術】

【0002】

最近の日本における大型航空機の事故の大半は乱気流が原因とされ、より安全な航空機運航の実現には、乱気流による突然の揺れを防止する制御技術、すなわち外乱である乱気流の事前情報を何らかの方法で取得し、その情報を有効利用することで乱気流による航空機の揺れを抑制する制御（以下では、これを乱気流抑制制御と称す）が必要である。その制御アルゴリズムは、従来からいくつかの提案がなされている。例えば、特許文献１には「望ましくない横向き運動を減少するための方向舵修正コマンドを生成する方法およびシステム」が提示されている。この発明の目的は、乱気流および突風によって引起こされる横からの荷重および狂いを減少させることによって飛行機の望ましくない横向き運動を減少させることにあり、そのために、乱気流および突風によって引起こされる、飛行機の垂直安定板での有効力を軽減するように飛行機の方向舵コマンドを変化させる方向舵修正コマンドが生成される構成を採っている。より特定的には、垂直安定板の両側で差圧が測定され、ロールレートおよびヨーレート補償される方向舵偏向値を生成するのに用いられる。補償された偏向値は、飛行機のダッチロール周波数の２５％である折点周波数で高域フィルタ処理される。結果は、第２の方向舵偏向値に減算的に合成される第１の方向舵偏向値である。第２の方向舵偏向値は、飛行機のヨーダンパによって生成される慣性横滑り角度比率値を利得調節し、かつ低域フィルタ処理することによって引出される。また、特許文献２である「航空機の上流の乱気流を予測して測定するためのシステム」の発明の目的は、ライダーを使用してシステムに関して前進位相している、気流作動面の、飛行制御システムによる作動を可能にするに十分な距離にある航空機の前方の風速を測定することにある。この発明は、航空機の上流の乱気流を予め測定し、該航空機に搭載されて配置されたシステムに関し、このシステムは、前記航空機の前部に向けて光ビームを送信し、かつ散乱された光ビームを受信するための、例えば、紫外線ライダーのような、ライダーと、制御手段と関連付けられた直接検出デバイスと、第１内部補正アルゴリズムを使用する第１処理素子と、少なくとも１つの航空機制御面のアクチュエータに実行できる命令を送信する第２外部補正アルゴリズムを使用する第２処理素子と、を含むことを特徴とするシステムである。

【0003】

特許文献３には、飛行中、出会う少なくとも１つの垂直乱気流の航空機、特に、輸送機に対する影響を、特に効果的に軽減できる方法を提示する「航空機に対する垂直乱気流の影響を軽減するための方法と装置」が開示されている。この発明が提示する方法は、飛行中、以下の一連の連続する工程、
ａ) 航空機の外側の風の垂直構成要素が航空機の現在の位置で決定され、
ｂ) この風の垂直構成要素により、航空機の外側の垂直乱気流に関する激しさのレベルが上記の現在の位置で決定され、
ｃ) 上記の風の垂直構成要素により、少なくとも１つの制御指令が、航空機の揚力に作用できる少なくとも１つの制御可能な可動部材に対し計算され、この制御指令が、垂直乱気流により航空機に生じる負荷因子の大きさを減少することのできるものであり、
ｄ) 少なくとも上記の激しさのレベルに依存する作動状態が実現しているかどうか確認するため照合がなされ、
ｅ) 上記の作動状態が実現していれば、上記の制御指令が上記の制御可能な可動部材の少なくとも１つの作動器に伝達される、
が自動的および繰り返し実施されるというものである。

【0004】

特許文献４に開示されている「航空機の動的構造負荷の最小化」の目的は、その名称の通り航空機の動的構造負荷の効果的な最小化を提供することである。この発明は、外部励振によって航空機に対して導入される航空機の動的構造負荷を最小化する方法を提供するものであり、前記外部励振を示す信号を生成するステップと、航空機に対して導入される動的構造負荷を低減するように、先行制御ルールに従って、前記励振を示す信号から、航空機の制御要素を駆動するための先行制御信号を導出するステップと、前記先行制御の実績を示す誤差信号を生成するステップと、動的構造負荷を最小化するように、前記誤差信号および/または前記励振を示す信号によって先行制御ルールを最適化するステップとを含むものである。

【0005】

しかし、これら従来の乱気流による航空機の揺れを抑制する制御アルゴリズムは、せいぜいシステムの不確かさのみを考慮した問題に対する制御アルゴリズムまでであって、多くはそれすら考慮しない問題に対する制御アルゴリズムであった。更に進んだ外乱の事前情報データに対する観測誤差およびシステムの状態データに対する観測誤差を考慮した制御アルゴリズムは現在まで提案されていない。しかし、これらの観測値には、必ず観測誤差が含まれ、観測誤差を含まない観測を行うことは不可能である。また、これらの観測誤差を考慮した制御を行わないと、乱気流の影響を悪化させる、すなわち乱気流の揺れを増幅させる制御が加えられる可能性もあり、実システムに適用する場合には、これらの観測誤差およびシステムの不確かさを考慮したアルゴリズムが強く求められるところである。

【発明の概要】

【発明が解決しようとする課題】

【0006】

本発明の課題は、背景技術において指摘した二つの観測誤差、すなわち事前乱気流の観測誤差および自機の運動の状態観測誤差、およびシステムの不確かさを考慮した乱気流抑制制御を行うために、これらの誤差または不確かさの影響をすべて考慮し、これらがもっとも悪く影響し合った場合にも、遭遇する乱気流の影響を最小にする制御を行うアルゴリズムを開発し、提供することにある。

【課題を解決するための手段】

【0007】

本発明の乱気流抑制制御方法は、連続時間の航空機の線形運動方程式を離散化したシステムを制御対象とし、舵面を動かすアクチュエータや推力変化をもたらすエンジン等の動特性の離散時間システムとしてこれらの結合システムが与えられるものとし、システムのモデリング誤差に対するロバスト性（頑強性）を確保するために制御入力端に[式１]で示される不確かな無駄時間を有する遅れシステムを用い、事前乱気流情報の観測誤差に対するロバスト性を確保するために[式２]に示されるあらかじめ想定された誤差幅を用いた事前乱気流情報を用い、さらに自機の運動の状態の観測誤差に対するロバスト性を確保するために[式３]に示されるあらかじめ想定された誤差幅を用いた機体運動状態を用い、機体の運動に関する状態および入力についての制約条件下において、評価関数を最適化し、得られた最適入力時系列のうち、最初のデータを実際の入力として加え、次のステップでは、新たに得られた乱気流情報および機体運動に関する情報を用いて、再度最適化問題を解き、得られた最適入力時系列の最初のデータを実際の入力として加え、これを毎ステップ繰り返すことを特徴とする。

【数1】

ただし、上式において、Ｔ_d は不確かな遅れ時間に相当する遅れステップから構成された遅れステップ数の集合を、d,d+1,‥‥は各要素の遅れ時間を、ｄアンダーバーは時間ステップ数で表現された一番小さい遅れステップの値を、ｄアッパーバーは時間ステップ数で表現された一番大きい遅れステップの値を、ｗ_k+j|k はｋステップにおいてｊステップ先に遭遇する乱気流の観測値、ｗ_k+j はｋステップにおいてｊステップ先に遭遇する乱気流の真値、Ｘ_jΔ_ｗ は乱気流の観測誤差をあらわす。

【0008】

また、機械運動に関する状態量の観測誤差を加え、本発明の乱気流抑制制御方法は、上記の観測誤差および不確かさのすべての組み合わせに対して共通の入力時系列を求める定式化を行うものとした。

【発明の効果】

【0009】

本発明を用いる場合の効果は従来技術と比較して以下のとおりである。
１）事前乱気流の観測誤差があったとしても、想定範囲内ならば実際の乱気流を受けた場合の運動の影響は最適化の結果得られる評価関数値以下である。
２）自機の運動の状態の観測誤差があったとしても、想定範囲内ならば実際の乱気流を受けた場合の運動の影響は最適化の結果得られる評価関数値以下である。
３）搭載システムに不確かな遅れがある場合にも、想定範囲内ならば実際の乱気流を受けた場合の運動の影響は最適化の結果得られる評価関数値以下である。
４）エレベータ等の制御入力の駆動範囲および入力レートリミットを超えない制御入力が得られるため、これらの制御入力の飽和およびレートリミットによる機体運動の不安定化は生じない。

【図面の簡単な説明】

【0010】

【図1】本発明に係る乱気流抑制制御の概念図である。

【図2】制御入力端における不確定無駄時間の影響を示す図である。

【図3】観測誤差を含む事前外乱情報を示す図である。

【図4】本発明の突風抑制性能をシミュレーションした結果をまとめたグラフである。

【図5】本発明の突風抑制性能をシミュレーションした結果の一例を表すグラフである。

【発明を実施するための形態】

【0011】

図１に本発明の概念図を示す。本発明の制御方法では飛行中の航空機は何らかの方法により乱気流の風速・風向データ（Ｖg で表現）が取得でき、また、航空機の運動に関わる情報（Ｖ, θ, etc. で表現）も取得できるものとし、このような場合に、制御を行わないと気流の影響をまともに受けた動きが発生するが、その動きを押さえた動きにするような制御入力コマンドを、乱気流情報の観測誤差、機体運動の観測誤差および機体運動モデルの不確かさに対してロバスト（頑強）であるように生成する制御方法が本発明の基本構想である。

【0012】

本発明の第１の特徴点は、システムのモデリング誤差に対するロバスト性（頑強性）を確保するために制御入力端に不確かな無駄時間を有する遅れシステムを用い、事前乱気流情報の観測誤差に対するロバスト性を確保するためにあらかじめ想定された誤差幅を用いた事前乱気流情報を用い、さらに自機の運動の状態の観測誤差に対するロバスト性を確保するためにあらかじめ想定された誤差幅を用いた機体運動状態を用いた点にある。
第２の特徴点は、上記の観測誤差および不確かさのすべての組み合わせに対して共通の入力時系列を求める定式化を行った点にある。

【0013】

今、航空機の定常飛行状態からの微小運動を考える。その微小運動を表す機体運動モデルが離散時間システムで与えられるとする。ただし、これは通常与えられる連続時間システムを搭載計算機の周期Ｔ_s[s]により零次ホールドによる離散化したシステムでよい。また、エレベータ（昇降舵）等のアクチュエータシステムも同様にして、離散時間システムとして与えられるとする。これも、通常与えられる連続時間システムを搭載計算機の周期Ｔ_s[s]により零次ホールドによる離散化したシステムでよい。このとき、これらを結合したシステムが下記により与えられるとする。

【数2】

ここで、ｘ_k はｋステップにおける結合システムの状態、ｗ_k はｋステップにおける乱気流入力、ｕ_k はｋステップにおけるエレベータ等の制御入力、ｚ_k はｋステップにおける上下加速度の定常値からの変動を表す。さらに，ｕ_k を１ステップ前の入力 ｕ_k-1 とその差分Δｕ_k に分解し、ｕ_k＝ｕ_k-1＋Δｕ_k と表現する。
さらに、エレベータ等の制御入力デバイスによる入力が加えられる制御入力端において、Ｔ_d ∈［Ｔ_dmin，Ｔ_dmax］と不確定な遅れ時間を有する遅れ時間システムが与えられるとする。なお、この不確定な遅れはアクチュエータ等のモデリングの際に発生するモデル化誤差を記述するための仮想的なシステムであって、実際のシステムに遅れが存在しなくても良い。このとき、不確定な遅れ時間Ｔ_d は、先の離散時間システムＰにとって下記のＴ_d による不確かな遅れステップに相当する。

【数3】

そのため、制御入力端にＴ_dの不確かな遅れを有する離散時間システムＰはＴ_d によって表現された有限個の遅れと［数１］のシステムを結合した有限個の要素からなるシステム集合によって表現される。以上で、不確かな遅れを用いたシステムの不確かさ表現が可能となった。

【0014】

次に、乱気流の事前情報の観測誤差を考える。今、現在の時刻である kステップからN-1 ステップまでの乱気流情報が観測されるとする。

【数4】

ここで、ｗ_k+j|k はｋステップにおいてｊステップ先に遭遇する乱気流の観測値、ｗ_k+j はｋステップにおいてｊステップ先に遭遇する乱気流の真値、Ｘ_jΔ_ｗ は乱気流の観測誤差をあらわす。
次に、機体の運動に関する状態の観測誤差を考える。ｋステップの真の状態 x_k がx_k|k と観測され、真の状態と観測値にはｘ_ｋ＝ｘ_ｋ|ｋ＋Δ_ｘの関係が成り立つとする。ただし、観測誤差Δ_ｘはあらかじめ定められるとする。
このとき、不確かな遅れステップ［数２］、乱気流の観測誤差［数３］、状態の観測誤差のすべての組み合わせについて共通する入力時系列を、機体の運動に関する状態および入力についての制約条件下において、あらかじめ定められた評価関数を最適化するよう求める。なお、この問題は二次錐計画（Second Order Cone Programming：SOCP）問題として定式化され、市販のツールを用いて求めることができる。
得られた最適入力時系列のうち、最初のデータを実際の入力として加える。
次のステップでは、新たに得られた乱気流情報および機体運動に関する情報を用いて、再度最適化問題を解き、得られた最適入力時系列の最初のデータを実際の入力として加える。これを毎ステップ繰り返す。

【0015】

以下、本発明の実施の形態について、詳細に説明する。
本明細書で用いる記号の定義を以下に行う。０_n はｎ次元零行列、０_n,m はn×m 次元零行列、Ｉ_n はｎ次元単位行列、１_n はすべての要素が１であるｎ次元ベクトル、０ は適当なサイズの零行列もしくは零ベクトル、Ｒⁿ はｎ次元実ベクトルからなる集合、Ｒ^n×m はｎ×ｍ 次元実行列からなる集合、

【数5】

本明細書を通して、離散時間システム中の変数ｐ について，ｐ_k はｋステップにおける変数ｐ をあらわし、ｐ_k+j|k は、ｋステップにおいて予測されたｋ＋ｊステップにおける予測値をあらわすとする。

【0016】

本明細書では、線形化された航空機の微小縦運動を制御対象とする。また、広く使用されている文献と同様に乱気流の時間変化は機体運動に影響を及ぼさないと仮定する(非特許文献１参照)。現在の機体搭載コンピュータはほとんどディジタルコンピュータであることを考慮して、連続時間の航空機の線形運動方程式を離散化したシステムを制御対象とし、舵面を動かすアクチュエータや推力変化をもたらすエンジン等の動特性も離散時間システムとして与えられると仮定する。このとき、これらの離散時間システムを結合したシステムが（1）式で与えられるとする。

【数6】

ここで、ｘ_k ∈ Ｒ^ｎ, ｗ_k ∈ Ｒ^ｎｗ, ｕ_k ∈Ｒ^ｎｕ,ｚ_k ∈Ｒ^ｎｚはそれぞれ、システムＰの状態変数（ 具体的には速度、角速度、エレベータ舵角等）、外乱入力（具体的には乱気流）、制御入力（具体的には空力舵面変位コマンドなど）、制御出力（上下加速度などの乱気流に起因し、かつ軽減したい運動を特徴づける変数）をあらわすベクトルである。なお、状態ｘ_kは観測誤差を含むものの、すべて観測可能であると仮定する。また、サンプリング周期をＴ_s [s] とした。

【0017】

一般に、システムをモデリング誤差なく求めることは不可能であり、実際のシステムに対する制御系の設計はモデリング誤差等の不確かさを考慮して行われる。このシステムの不確かさの表現方法には、パラメトリックな表現や周波数領域における表現などがあるが、その中の一つの方法として制御入力端における有界な無駄時間を用いた時不変無駄時間システムによる不確かさ表現がある。この方法は、従来技術である非特許文献２ではトラッキング問題、非特許文献３ではモデルマッチング問題と、異なる設計仕様においてもその不確かさ表現としての有効性が示されている。そこで、本明細書においても、制御入力端における無駄時間を用いた不確かさ表現を用いる。具体的には、無駄時間Ｔ_d [s] がＴ_d ∈[Ｔ_dmin,Ｔ_dmax] である時不変無駄時間システムを制御入力端に加えたシステムを考える。ただし、搭載計算機がもたらす１ステップ遅れを考慮して、Ｔ_s≦Ｔ_dmin を満たすと仮定する。いま、制御入力端に無駄時間Ｔ_d が不確定な無駄時間システムを挿入しても、離散時間システムを対象としていることから、制御入力が加えられるステップの数は有限個に限られることがわかる（図２参照）。

【数7】

【0018】

さらに、制御入力ｕ_ｋ を（2）式のように、１ステップ前の値とその差分に分解して表現する。
ｕ_ｋ ＝ ｕ_ｋ−１＋Δｕ_ｋ （2）
これは、入力飽和に対する制約条件のみならず入力レートに対する制約条件も課すためである。
このとき、（1）式のシステムにＴ_d∈[Ｔ_dmin,Ｔ_dmax] を満たす不確かな時不変無駄時間システムを制御入力端に加えたシステムＰ_ｕ は、Δｕ_k を入力とする以下のシステムとして記述される。

【数8】

【0019】

本明細書では、システム（1）の状態量ｘ_ｋ は、観測誤差を含めて観測されるとする。現在時刻をｋステップとし、時刻ｋにおいて観測される状態量ｘ_ｋ|ｋは、観測誤差をあらわす既知の時不変ベクトルΔ_ｘ を用いた以下の式を満たすと仮定する。
ｘ_k ＝ｘ_ｋ|ｋ＋Δ_ｘ (5)
ただし、Δ_ｘの存在範囲は既知の凸多面体で与えられ、その端点から構成された集合をΦ_x とする。すなわち、以下の式が成り立つ。
Δ_ｘ ∈ coΦ_x (6)
上式のcoはコンベックス集合すなわち、coΦ_xは凸多面体の全体及び内部を表す。
ここで、ｘ_k の存在範囲 Ω_x を定義する。

【数9】

【0020】

また、外乱入力である乱気流について、現在からＮ−１ステップ先まで観測誤差を含んだ情報として得られるとする。すなわち、以下の式が成り立つとする。

【数10】

ここで、ｗ_k+j|k はｋステップにおいてｊステップ先に遭遇する乱気流の観測値、ｗ_k+jはｋステップにおいてｊステップ先に遭遇する乱気流の真値、Ｘ_j∈ Ｒ^nw×nw (j ＝0, 1, ・・・N-1) は事前外乱情報の観測誤差をあらわすための与えられた定数行列であり、Δ_ｗ は以下の式を満たす不確かな定数ベクトルとする（図３参照）。

【数11】

【0021】

事前外乱情報がＮ−１ステップ先まで得られることから、Ｎ−１ステップ先まで状態予測を行うとする。

【数12】

上記の仮定のもと、（3）式および（4）式に定義された不確かさを有するシステムＰ_ｕに対して、次に示すような最適化問題（問題１）を考える。すなわち、(5)式のように観測誤差を含む機体運動に関する状態量が観測でき、また(9)式のように観測誤差を含む乱気流情報がえられる場合に、(13),(14),(15)式の制約条件のもと、評価関数J_i(xⁱ, Δ_{u, z} ⁱ) の最大値（具体的には事前乱気流の観測誤差、機械運動の状態量の観測誤差に対して最大化した J_i(xⁱ,Δ_{u, z} ⁱ)のｉについての最大値）を最小にする共通入力時系列 Δ_u を求めるものである。

【数13】

ここで重要なことは、システムＰ_ｕ^ｉのインデックスｉ，事前外乱情報の観測誤差集合Ω_w，状態の観測誤差集合Ω_x に依存しない共通入力差分時系列を要求していることである。その結果、システム集合Ｐ_ｕ によって記述される不確かさ、事前外乱情報の観測誤差および機体運動状態量の観測値に対する観測誤差に対してロバストな入力時系列を求めることが可能となる。

【0022】

従来の乱気流抑制制御である非特許文献５の手法は、「(5)式とは異なり、観測誤差がなく機体運動に関する状態量が観測でき、また、(9)式とは異なり、観測誤差がなく乱気流情報がえられる場合に、(13),(14),(15)式の制約条件のもと、評価関数J_ｉ(xⁱ, Δ_{u, z} ⁱ) の最大値を最小にする共通入力時系列 Δ_u を求めなさい。」という問題を解く方法を提案している。ここでは、インデックスi は、システムの不確かさによるもので、この一種類の不確かさのみを考慮した「最大値」を扱っている。これを一歩進めた乱気流観測誤差のみを考慮した乱気流抑制制御アルゴリズムでは、「(5)式とは異なり、観測誤差がなく機体運動に関する状態量が観測できるが、(9)式と同様に、観測誤差を含む乱気流情報がえられる場合に、(13),(14),(15)式の制約条件のもと、評価関数J_ｉ(xⁱ,Δ_{u, z} ⁱ)の最大値を最小にする共通入力時系列 Δ_u を求めなさい。」という問題を解く方法を提案している。それが、問題２を毎ステップ解くという制御手法である。
すなわち、非特許文献５では、事前外乱情報および機体運動状態量の観測値に観測誤差がない場合の定式化を行っており、乱気流観測誤差のみを考慮した乱気流抑制制御アルゴリズムでは、非特許文献５の結果を発展させ、事前外乱情報のみに観測誤差がある場合の定式化を行っている。よって、本明細書ではその研究の結果のみ示す。以下では、機体運動状態量に観測誤差がないことから、（5）式について，ｘ_ｋ＝ ｘ_ｋ|ｋ が成り立つとする。

【数14】

【0023】

すなわち、非特許文献５では、以上の準備のもと、問題１は次の問題２となる。

【数15】

以上をまとめると、機械運動に関する状態量の観測誤差が無い場合の問題１に対する制御手法は、毎ステップごとに更新される（8）式の拡大状態と（9）式の乱気流情報を用いた問題２ を毎ステップ解き、Δｕ_ｋ＝ Δｕ_ｋ|ｋとした（2）式による入力を用いる制御手法である。

【0024】

【数16】

この問題３は問題２に機体運動に関する状態量の観測誤差を加えた問題となっている。そのため、問題１と問題３は同じものであり、表現を変更しただけである。
このとき、問題３は、(Δ_ｗ,Δ_ｘ)のすべての端点の組み合わせにおいて評価するため、計算複雑度は問題２よりも大きい。しかし、プラントの不確かさが制御入力端の不確定な時不変無駄時間によって記述され、事前外乱情報および機体運動状態量の観測値に対する観測誤差がそれぞれ（9）式および（5）式で表現されるならば、問題３の定式化の際には、何ら近似等の保守性を導入していないことに注意されたい。また、問題３は状態等の制約条件(19)および(20)を満たす制御入力時系列の中から最適な入力時系列を作成する。以上のことから、問題３を用いたシミュレーションを行うことで、突風軽減制御の性能限界評価を行うことが可能である。

【実施例】

【0025】

本節では、非特許文献５と同様に、乱気流による上下加速度の抑制可能性について調べた一例を示す。
非特許文献５と同様に、宇宙航空研究開発機構が有する実験用航空機MuPAL-αの1520 [ｍ]、真対気速度66.5 [ｍ/ｓ]における水平定常飛行状態から線形近似を行った微小縦運動に１/(0.1ｓ＋１) の一次遅れモデルによってモデル化されたエレベータアクチュエータモデルを加えた航空機運動を考える。サンプリング周期Ｔ_sを0.1 [ｓ] とした。なお、状態変数ｘ_ｋ は[ｕ_ｉ ｗ_ｉ q θ δ_e]^Ｔ，外乱入力ｘ_ｋ はｗ_ｇ，制御入力ｕ_ｋ はδ_ｅｃ ，制御出力ｚ_ｋはΔａ_ｚ である。
ここで、ｕ_ｉ，ｗ_ｉ，q，θ，δ_e，ｗ_ｇ，δ_ｅｃ，Δａ_ｚはそれぞれ機体座標系ｘ方向慣性速度[ｍ/ｓ]，機体座標系ｚ方向慣性速度[ｍ/ｓ]，ピッチ角速度[rad/ｓ]，ピッチ角[rad]，エレベータ舵角変位[rad]，慣性座標系ｚ方向突風[ｍ/ｓ]，エレベータ舵角変位コマンド[rad]，慣性座標系ｚ方向加速度の変分[ｍ/ｓ^２] をあらわす。なお，ｗ_ｇは機体下方から上方へ吹く風をプラスの風としている。

【0026】

不確定無駄時間Ｔ_d[s] は、Ｔ_d∈ [0.1, 0.4] とした。このとき、

【数17】

機体運動状態量の観測値に対する観測誤差をあらわすベクトルΔ_x の存在範囲を定める端点集合Φ_x は

【数18】

評価関数Ｊ_ｉ中の重み行列Ｑ, Ｓ は共に零行列と定めた。
なお、（15）式の制御出力制約は課さず、制御入力レートに対する重みであるＲ，および事前乱気流情報の時間長さＴ_ｓ×Ｎはシミュレーションパラメータとした。
また、遭遇する乱気流は次式で与えられると設定した。
ｗ_ｇ(t) = sin(ωt) (21)
ここで，ｔ[s] はシミュレーション開始からの時刻をあらわし、乱気流の周波数であるω[rad/s] はシミュレーションパラメータとした。
事前外乱情報の観測誤差をあらわす行列Ｘ_ｊ は以下に定めた。
Ｘ_ｊ ＝ 0.2 ＋ 0.1 ×(66.5/100 ×Ｔ_ｓ) ｊ (22)
シミュレーションパラメータであるＲ，Ｔａ，ω の集合はそれぞれ以下に定義した。
Ｒ ∈｛ 10¹, 10², 10³,10⁴, 10⁵ ｝
Ｔ_a∈ {1.0, 1.2, 1.4,1.6, 1.8, 2.0, 2.2, 2.4}
ω ∈ {0.1, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0,7.0}
機体運動線形時不変システムおよびエレベータアクチュエータモデルを用いた20 [ｓ] のシミュレーションを行った。なお、ＳＯＣＰ 問題の解は非特許文献４のYALMIPを介した非特許文献６のSeDuMiにより求めた。また、比較対象として事前乱気流情報が得られない条件下における問題３による乱気流抑制制御を行うケースも同一環境下においてシミュレーションを行った。記述上、ケース名を以下に定める。

【0027】

【数19】

この図４はカラー表示で作られたグラフであるが特許図面ではカラー表示ができないため、添付資料としてカラーグラフを別途提出する。上記のシミュレーション回数で２×２倍としているのは乱気流情報誤差について最大と最小を、機体状態量誤差について最大と最小のケースの組み合わせを採るためである。図４に示したものは外乱入力ω(乱気流の周波数)が1.0，2.0，3.0，4.0，5.0，6.0，7.0についての結果である。なお、このグラフでは性能比較のため、ステップ数Ｎをｙ軸に、制御入力レートに対する重みＲをシミュレーションパラメータとしてｘ軸にとり、乱気流の抑制可能性を判断するＪ_Ａ/Ｊ_Ｂと、事前乱気流情報の有用性をあらわすＪ_Ａ/Ｊ_Ｃをｚ軸方向にプロットした。また、同等の性能をあらわすＪ_Ａ/Ｊ_Ｂ ＝ 1.0 およびＪ_Ａ/Ｊ_Ｃ ＝ 1.0 をあらわすメッシュ状の面もプロットした。このとき、メッシュ状の面より下にプロットされているＲとＴ_ａの組み合わせは、Ｊ_Ａのほうが性能が良いことを意味する。ω＝ 0.1, 0.5 の場合では、ほとんどのケースにおいてＪ_Ａ/Ｊ_Ｂ，Ｊ_Ａ/Ｊ_Ｃ共に２を超えていたため、省略した。この結果より、本稿で想定した問題設定においては、周波数が2.0 〜 6.0 [rad/s] の乱気流に対してはその影響を小さくすることが可能であると考えられるが、周波数が1.0[rad/s] より小さい乱気流と7.0 [rads/] より大きい乱気流の影響は、事前に乱気流情報を得てもほとんど抑制できないと考えられる。

【0028】

端的なシミュレーション結果の例として、乱気流の周波数ω＝3.0，Ｒ＝10，Ｔａ＝2.4(Ｎ＝24) の場合のシミュレーション結果を図５に示す。なお、シミュレーションは、事前乱気流の観測誤差が最大もしくは最小の場合、機体運動の状態量の観測誤差が最大もしくは最小の場合、不確かな遅れステップ数が１もしくは２もしくは３もしくは４の場合が存在し、これらすべての組み合わせについて実施している。これらのシミュレーション結果からCaseＡ，Ｂ，Ｃそれぞれの値JA,JB, JC を求め，JA/JB，JA/JC をプロットした点が、それぞれ、図４の ω＝3.0 におけるXおよびYである。図５において実線が上記のCaseＡ（本発明の乱気流抑制制御のデータ）を、破線が上記のCaseＢ（制御なしのデータ）を、点線が上記のCaseＣ（事前乱気流情報ナシのデータ）を示している。グラフに示した１６組すべてについて、CaseＡは制御を行わないCaseＢの場合に比べ、制御幅が顕著に小さいことが分かる。また、点線表示のCaseＣはCaseＢの場合に比べればやや制御幅が小さくなっているが、CaseＡはこれよりはるかに制御幅が小さくなっている。本発明の乱気流抑制制御を施したCaseＡが如何に機体運動を抑制しているかが実線のグラフと破線のグラフ及び点線のグラフとの対比で一目瞭然に確認できる。

【産業上の利用可能性】

【0029】

本発明の乱気流抑制制御方法は、外乱を事前に観測でき、外乱の影響を最小化するなどの最適化仕様が与えられたシステム全般、例えば、船は波浪によって揺れるが、マストなどからの画像によって、今から遭遇する波浪情報が事前に得られると仮定する。そのような場合には、事前に得られた波浪情報を元に、もっとも揺れを抑制する操舵入力を得ることが可能となる。
また、別の応用として、風力発電を行う風車の最適な羽根ピッチ角制御を行うことも可能である。これは，風力発電を行う風車は回転数を羽根のピッチ角によって制御しているが、風車をある一定の回転数で回すことが重要であり、突風などの影響を受けない回転数維持が望まれている。このような場合に、風車の脚もしくは風車の中心などに遠方の乱気流を計測する装置を追加することで、外乱である風速の揺れを観測することが可能となり、事前に得られる風速の揺れを考慮した発電に最適なピッチ角を計算することが可能となる。

【先行技術文献】

【特許文献】

【0030】

【特許文献1】特開平８−３１０４９５号公報 「望ましくない横向き運動を減少するための方向舵修正コマンドを生成する方法およびシステム」 平成８年１１月２６日公開

【特許文献2】特表２００８−５００５２５号公報 「航空機の上流の乱気流を予測して測定するためのシステム」 平成２０年１月１０日公開

【特許文献3】特表２００９−５１１３３９号公報 「航空機に対する垂直乱気流の影響を軽減するための方法と装置」 平成２１年３月１９日公開

【特許文献4】特表２００９−５２４５４７号公報 「航空機の動的構造負荷の最小化」 平成２１年７月２日公開

【非特許文献】

【0031】

【非特許文献1】Stengel, R.E.:”Flight Dynamics”Princeton University Press, 2004,pp.274-297.

【非特許文献2】大野正博，山口恭弘，畑剛，高浜盛雄，宮澤与和，泉達司「ALFLEXのロバスト飛行制御則設計」，計測自動制御学会論文集，Vol. 34, (1998年)，pp. 1905-1912．

【非特許文献3】佐藤昌之，佐藤淳「突風応答と操縦応答の模擬を行う飛行制御器設計− MuPAL-α の横／方向運動に対する設計とその実験結果」日本航空宇宙学会論文集，Vol.54,(2006年)，pp. 71-81.

【非特許文献4】Lofberg, J.:YALMIP: A Toolbox for Modeling and Optimization in MATLAB, Proc. the CACSD Conference, Taipei, Taiwan, 2004,

【非特許文献5】佐藤昌之，横山信宏，佐藤淳「乱気流の事前情報を用いたロバストモデル予測制御によるGust Alleviation制御」日本航空宇宙学会論文集,Vol.57（2009年9月28日），pp. 345-353．

【非特許文献6】Sturm, J.S.: Using SeDuMi1.02, a MATLAB Toolbox for Optimization Over Symmetric Cones, Optimization Methods and Software, Vols. 11-12, (1999), pp. 625-653.

【図1】

【図2】

【図3】

【図5】

【図4】