人力飛行機空力設計と粒子群最適化

ご無沙汰しておりました。Yuukivelです。
直近の前刷り原稿締め切りでバタバタして、お盆に入って実家でダラダラしております。
一人暮らしだとなかなか好物の梨とか食べられないし、やはり食の充実という意味で実家はいいですね!

さて、、

今回は以前記事にした人力飛行機設計について、OBチームにて実際に作ることになった機体の空力設計について解説する。
前記事はこちら→ http://d.hatena.ne.jp/yuukivel/20140101/1388592540

先に三面図を公開しておきましょう。ちょうど良い時期ですし、様々に参考にしていただければ。

一応PDFとXFLR5ファイルについては、まだ変更の可能性があるので今後公開するかも、としておく。

最適化について

本題の空力設計の話に入る前に、前提の知識として、「最適化」について解説しておく。
TLを見ていると、最適化を研究している者としてかなり気になる用法でこの言葉が使われていることがある。
あまり詳しく立ち入ると記事が長くなってしまうので、ここでは簡単に。

最適化とは、「ある評価関数があって、その評価関数を実行可能な範囲で最大・最小化すること」だ。
一応、最大化でも最小化でも、評価関数の正負を逆にすれば最大化にも最小化にもできるので、ここでは全て最小化、とする。

実行可能な範囲全てをくまなく調べて、本当にここより小さな解は存在しない、という解を大域的最適解という。
また、実行可能な範囲で、その解の周りだけみると最小となっている解を局所的最適解という。

ある解が大域的最適解だと言い切るためには、評価関数が以下である必要がある。
・実行可能な範囲で、関数の形が分かっている
・実行可能な範囲で、凸関数であることが証明できる

凸関数とは、例えば二次関数のように、谷が1つしか無いような関数のこと。


上のグラフはy = 0.5x^4-0.2x^3-3x^2+4x+2を表している。
このうち赤い○で囲ったのが「xが-3から3の範囲で」大域的最適解といえる。ちなみにこの関数であれば-3から0の範囲で凸関数。2階微分がこの範囲は正だからだ。
しかし、もしかしたらxが遙か彼方で谷を持つかもしれない。つまりは真の意味での大域的最適解は「神のみぞ知る」解といえる。

人間は厳密な意味での大域的最適解を得るのは不可能であり、ましてや関数の形が分からないし凸関数であることも証明できない問題には局所最適解で甘んじるしかないのだ。
局所最適解であることの必要十分条件としてKKT条件(Karush-Kuhn-Tucker condition)がある。これを満たせば、一応局所最適解ですよということは示される。
KKT条件について詳しくは立ち入らないが、実行可能かどうか、その点でのラグランジアンの傾きが0かどうか、といった条件をひとまとめにしたもの。

言いたかったことはなにかというと
・「最適化」とは評価関数があって、実行可能な範囲で局所最適解を求めること、ということ。それ「改善」じゃない?と言いたくなることが多々。
・最適性の条件を満たしていればそれはその評価関数において最適解であるということ。あとこの部分性能悪いから、それ最適解じゃないよとかいうのは、ちょっとこまる。だいたい評価関数の選び方はコンセプトで決まってしまうのだから。

最適解を求める方法には大きく分けて二種類あって
・勾配法(ニュートン法とか逐次二次計画法(SQP)とか) 
関数の傾きの情報を使って(局所)最適解を探索する。関数が微分可能でないと使えない。速いと言われるけどヤコビアンやヘッシアンが分かっていないとそれらを評価するのに時間がかかる。
・メタヒューリスティクス(遺伝的アルゴリズムとか粒子群最適化とか)
ランダムサーチの強化版みたいなもの。理論的根拠は希薄だが上手くいく。厳密な局所最適解は得られないが、関数が連続ですらなくても実行可能でロバスト。ちなみにXFOILを用いた最適化はFortranの小数点の問題(丸め誤差)のせいでメタヒューリスティクスでしか動かせません。

これらのうち長所短所をみながら問題に合わせて使っていく。

で、今回はメタヒューリスティクスの中でも粒子群最適化(Particle Swarm Optimizer:PSO)をつかって人力飛行機を設計しましたよというお話

空力設計(コンセプト)

この機体の空力設計におけるコンセプトは、以前の記事で作ったものと変わっていなくて、
・構造優先設計
・抵抗の徹底的な削減(薄翼化・高いテーパー比の尾翼等)
ということ。抵抗の削減であって出力の最小化ではないことも注意。取り回しの関係で翼幅は33.75mに固定。また設計速度は風に負けないように少々早めの7.5m/s。こうした制約の中で、そして存在する桁を使って、できるかぎり抵抗の小さい機体を設計しましょう、というコンセプト。

前者については以前の記事の通りで、桁を(上手い人に)設計してもらって、実際に作って重量を実測した上で、循環分布を設計した。
桁の剛性は重さの割に結構高く、たわみ2.5mを実現する循環分布を求めると楕円循環分布よりほんの少し揚力を真ん中によせた形になった。
そうこうして、各翼の端での循環の値が求まったとして、次の「α制約分散PSOによる主翼自動設計」の項に行きましょう。

空力設計(α制約分散PSOによる主翼自動設計)

各翼での循環値が定まったら、その循環値と指定した翼厚をみたす翼型の中から最小抵抗の翼型を設計する。
最適化問題風に記述すると以下のようになる。

翼型については4つの翼型の線形和と追加キャンバーで表現する。

1.使用する4翼型の翼型座標ファイルを用意しx座標の点の位置をそろえておく
2.係数に沿って翼型座標ファイルのy座標を単純に足し算する(引き算も可)
3.係数に従ってNACA4digit系列のキャンバーラインを作成する。
4.線形和で得た翼型のy座標に3.で得たキャンバーを足し算する。

以上の手順で翼型をパラメーターで表現する。以前GAでやっていたころは追加キャンバーをPARSECの応用でやっていたが、今回はNACA4digit系のキャンバーを利用した。
こうすることで得られる翼型がぐんと滑らかになった。

で、こうしてパラメータ表現された翼型に対しPSOを適用するのだが、そもそもPSOとは、、
Wikipediaの記事がシンプルでよいのでリンクを張っておきます。 粒子群最適化
PSOは魚群や鳥群の餌に群がる際の情報伝達をモデルとしたメタヒューリスティクスであり、同じメタヒューリスティクスである遺伝的アルゴリズム(GA)に比べて、実数問題や多峰性問題に向いていると言われる。
今回はPSOを並列実行できるように拡張し、octaveを7プロセス並列で一世代あたり490個体で計算させた。
PSOは各世代の位置と速度を、以下の式に従って更新する。

このときlbestは各プロセスでの過去から現在までの最良個体値、gbestは全プロセスを通した最良個体値とする。
ωやC,σは適当な定数。挙動をみながら調整する必要有り。

またgbestやlbestを決める上で最良個体を選ばなければならないが、今回は制約条件が存在するため、αレベル比較を実装した。
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1457-30.pdf
αレベル比較とは、制約条件の「満足度」を定義し、比較する両個体の満足度がともにα以上であれば目的関数の比較を、どちらかでもα以下であれば満足度の大小で比較を行う方法。
詳しい実装方法は論文に載っているので省略します。

このようにPSOにαレベル比較を実装することで、「欲しい循環を満たしながら、抵抗が最小となる翼型」を探索することが可能となった。
従って、
1.循環分布を各セクションに対して求める。
2.循環γ=0.5UCCLを満たす翼型を順次求める
3.各翼の中身を線形補間で翼型を求める。
のループを自動で回すことで、可能な限り抵抗の小さい主翼を求めることができる。
もちろんメタヒューリスティクスは厳密な最適解を求めることはできないから、実行毎に解が変わる可能性がある。
そこは自分が満足するのが出てきたらおわり、ということで(笑)
主翼最適設計とせず、主翼自動設計と言っているのは、こうしたところに理由がある。
1⇒2⇒3のループを繰り返し実行して、解を求めるようにすれば主翼最適設計といっても口がむずがゆくならない。

実際に設計を行った様子を動画にあげているので、もしよろしければ。

(見返したら翼自動最適設計といっていてむずがゆい、、確かに循環最適⇒翼型1最適⇒翼型2最適・・・となっているはずなので最適設計ではあるのだけれど)

おわりに

尾翼やペラの設計についても話したかったのだけれど、書く力が尽きてしまったので、いったんここで記事を締めくくろうと思う。
実際XFLR5の安定性解析結果をぱっとUnityで作ったフライトシミュに読み込んで安定性を議論・・・とかやっててそっちの話もしたいのだけれど、また今度ということで。。。
もうとにかく最適設計、最適設計、最適設計とやってきた機体なので実際にどう飛ぶか非常に楽しみ。

机上の空論となるか、計算は人間を追い越すか。

みたいな煽り文句をつければ、話題のタネになるかな(笑)



追記)三面図の厨二版あるんだけど、あまりに評判が悪いので悲しんでいます。上手く表現したつもりだったんだけどなぁ。。