JP2005316614A - 最適化方法及び最適化プログラム - Google Patents

Abstract

【課題】 従来方法では最適解の探索が困難な最適化問題について、確率過程を用いて目的関数を最適化する最適解の高精度近似解を求める最適化手法及びプログラムを提供する。
【解決手段】 機体１に対して飛行経路１〜Ｋ（ｕ，ｍａｘ）における経路積分が実行され、機体１に対する期待値経路が求められる。機体２、・・・機体ｋ、・・・機体Ｋｍａｘについても同様にそれぞれ、飛行経路にしたがって経路積分が実行され、各機体について期待値経路が求められる。全機を使って、期待値・機体を算出することで、機体の最適解に近似する期待値が求められる。また、全機体の期待値経路を使って、期待値機体・経路を算出する。同様に、最適解がパレート集合となる静的な多目的最適化問題についても、複数の目的関数を同時に最適化する最適解の高精度近似解が求められる。
【選択図】 図１２

Description

この発明は、目的関数が多峰性関数であっても最適解の近似解を求めることを可能にする単目的又は多目的の最適化方法及び最適化プログラムに関する。

システムを設計・構築する場合、工学においては常に何らかの最適化が要求される。最適化手法としては、従来、多くの手法が提案されてきている。最適化手法の一例として、山登り法がある。山登り法は、まず初めに適当な解Ｓを生成し、次に、解Ｓの一部を手直しすることにより解Ｓ′を生成する。ここでもしＳ′の目的関数値がＳの目的関数値よりも良ければＳをＳ′に更新し、もし良くなければ別のＳ′を生成する。このような操作をＳよりも目的関数がよくなるようなＳ′を作ることが不可能になるまで続け、そのときのＳを最適解として出力する。多くの局所解が存在するような複雑な解空間の場合には、山登り法で求められる解は殆どの場合に局所解である。

変分原理を応用した手法に代表される確定論に基づく最適化手法は、厳密解が得られるという利点があるが、解空間が複雑である場合には、最適解への効率的な探索をするために多くの工夫や知識・経験が必要である。更に、目的関数が多峰性関数である場合には、いずれの従来の方法も、最適解を求めるのに長い計算時間と計算機内に大容量のメモリを必要とするなどの取り扱い上の制約があることが多い。目的関数の解空間が複雑な場合には、局所解に陥る危険性が高く、様々な工夫をしても、場合によっては解を求めることが事実上不可能になることもある。

山登り法を改良し、厳密な最適解を得るための手法として、統計力学的な機構を形式的に最適化問題の解法に取り入れたシミュレーテッドアニーリング法が提案されている。シミュレーテッドアニーリング法は、解Ｓから解Ｓ′に移動するかどうかを判定する際に、最適値で確率最大となる確率分布を導入し、パラメータを操作しながら最適値（確率分布の頂点）を数値的に探索する手法であって、Ｓ′に対する目的関数値が悪い場合でもある確率でＳ′に移動することを許容する（良くなる場合にはもちろん移動する）。広い解空間を探索することができるが、目的関数値を悪化させる解へ移動（遷移）する確率は、探索が進むに連れて小さくされる。この確率の変化の仕方を適当にコントロールすることにより、局所解にとらわれずに厳密な最適解を得ることができ、解空間が複雑な問題では山登り法よりもかなり良い解を得ることができるとされている。しかしながら、クーリングスケジュールと呼ばれるパラメータ操作が不可欠であるとともに、問題によっては厳密な最適解に収束するには莫大な時間が必要となることがある。

量子的な振る舞いをする粒子が採る運動経路は古典軌道（古典力学で求められる軌道）上で最も確率が高くなる。即ち、古典力学における最小作用の原理に基づく式（１）に示すラグランジェ関数Ｌの作用積分Ｉが最小となる運動において最大となる。古典力学における粒子の運動経路はこの最も出現確率の高い運動経路と一致するので、これの近似解を期待値として経路積分を用いて計算することが可能である。更に、確定論として扱われる他の様々な問題についても、経路積分を用いて扱うことが可能である。与えられた境界条件を満たす無数の経路を用意したとき、粒子がそれぞれの経路を採った場合の確率（厳密には確率振幅）は式（１）によって定まり、更にその経路の期待値は多重積分として表される（１９４８年、ファインマン（Ｒ．Ｆｅｙｎｍａｎ）が発表）。したがって、最も出現確率の高い運動経路の近似解は、そうした経路積分による期待値（確率平均）として知ることができる。以上のように、経路積分により求められる期待値は、作用積分Ｉを最小とする経路の高精度近似解となる。

そこで、本発明者らは、量子力学で使われている経路積分の考えに基づいて、最適化問題を確率過程に置き換え、厳密解ではないが期待値という形で最適解の近似解を求める手法を提案している（非特許文献１〜３）。この近似解を初期解として確定論に基づく手法を適用すれば、局所解に捕らわれることなく厳密解を得ることができる。この手法は、アルゴリズムが簡潔であり、人為的に設定するべきパラメータが唯一つであることが特徴であり、解空間が比較的滑らかであることが期待できる工学的諸問題においては、大域解についての良好な近似解を得ることができる。

上記の確率過程を用いた最適化法においては、経路積分の発想を基にして、最適化すべき変数の発生確率Ｐを式（２）の形で定義する。式（２）は、式（１）に示す目的関数Ｉが最小のとき最大の発生確率をもち、Ｉが大きくなるほど、発生確率が急速に小さくなることを示している。

ここで、ｈは量子力学ではプランクの定数であり、ｈ→０の極限で古典力学系に一致するものであるが、ここでは、解の揺らぎの大きさを与える任意のパラメータとして用いられる。ｈが大であるならば広く緩やかな確率分布を、またｈが小であるならば狭く急峻な確率分布を与える。Ａは、確率の総和を１にするための正規化定数である。目的関数Ｉの大きさが大きいほど確率Ｐは小さい値となり、目的関数Ｉの値が小さいほど確率Ｐは大きい値となる。

一次元の運動を例に取ると、図１９に示すように、時間軸をＮ個に分割して、２点（ｘ0 ，ｘN ）間を結ぶ経路（一例）を描くことができる。確定論からすれば、ニュートンの運動方程式から定まるように厳密解（作用積分Ｉを最小化する解）がただ一つ存在する。一方、量子力学的には、Ｎ個に分割された時間軸（横軸）上の各時点ｔi （０≦ｉ≦Ｎ）における各点（縦軸）ｘi をそれぞれ直線で結んでできる無数の経路が存在する。ランダムに経路ｘ（ｔ）を発生させて式（１）からその経路に対応した作用積分Ｉの値を求める。即ち、ランダムに発生される非常に多数の経路ｘ（ｔ）については、式（１）によって求められた作用積分Ｉを用いて、式（２）に従って求められる発生確率が与えられる。本最適化法では、作用積分Ｉを最適化（最小化）すべき目的関数で形式的に置き換え、目的関数の最小値（最適値）を与える最適解ｘ（ｔ）の近似解として、式（３）で得られる期待値を求めようとするものである。

ここでＤｘ＝ｄｘ1 ｄｘ2 ・・・ｄｘN は、経路についての多重積分の測度であり、図１９に示す一次元の運動を例に採れば、時点ｔi における位置ｘi の期待値は、ｘi にその発生確率を掛けて総和を取ったものであり、連続量としてとらえれば位置ｘにその位置での発生確率Ｐを乗じつつ積分（多重積分）したものである。なお、ｘは、最適化すべき一般の変数であり、必ずしも経路である必要はない。

多重積分で計算される期待値は、モンテカルロ法を用いた数値計算に載せやすいという特徴を有している。上記のように、期待値の計算には多重積分が必要となるが、確率分布に従って設計変数を自動生成し、それらの平均値を取ることで、期待値が求められる。図２０は、一様重力場中の投げ上げ問題において、初期条件（ｘ0 ＝０，ｔ0 ＝０）、終端条件（ｘN ＝１，ｔN ＝１）、質量ｍ＝２ｋｇとしたとき、解析解（太線）と、確率分布に従って自動生成された経路（細線）とが示されており、解析解付近に評価値の高い経路が集中しているのが理解される。このように、任意の目的関数Ｉを用いて確率Ｐを定義し、これに基づいて発生させた多数の経路について期待値を求めれば、目的関数を最も小さくする変数（最適解）の近似解が得られる。

図２１は上記の計算の具体的手順を示すフローチャートである。スタート（ステップ４０、「Ｓ４０」と略す。以下同じ）して初期化（Ｓ４１）の後、乱数を発生して適当な解を作成して初期解とする（Ｓ４２）。解の評価値を計算して、その評価値に基づく確率分布を計算する（Ｓ４３）。Ｓ４３で求められた確率に従って解が発生するように、例えばメトロポリス法によって解の採否を決定する（Ｓ４４）。次に、期待値の計算をする（Ｓ４５）。解の生成過程はマルコフ過程と見なせるので、計算の終了は、計算の繰り返し回数Ｋが所定値Ｋｍａｘを超えたか否かで判定（Ｓ４６）する。繰り返し回数Ｋが所定値Ｋｍａｘを超えない場合にはＫの値を更新してＳ４２に戻り、繰り返し回数Ｋが所定値Ｋｍａｘを超えた場合には終了する（Ｓ４７）。

この最適化手法によれば、確定論に基づく最適化が初期値をどこに設定するかに応じて局所的な解に捕らわれるのに対して、原理的に初期値依存性がない。また、確率的平均値を積分によって得るだけであるので、最適化計算を複数の計算機で並列計算が可能である。また、設定すべきパラメータは揺らぎの程度を示すｈだけであり、しかも解の発散や振動を引き起こすこともない。また、拘束条件によって最適化すべきシステムの計算領域が限定されているほど計算の実行が容易になる。更に、状態量拘束条件への対応や知見の得られていないシステムへの最適化にも適用可能であるという優れた特徴を備えている。

工学問題においては、対象となるシステムに要求される目的は一つではなく、複数存在することが多い。複数の目的が存在する場合の目的関数を最適化する問題は多目的最適化と称されている。多目的最適化においては、目的間でトレードオフの関係、即ち、ある一つの目的のみを最適化すると、他の目的は望ましくない状態になるという関係が生じることが多く、その結果、解はパレート集合となる。そのため、多目的問題を形式的に単一目的問題に変形して取り扱うことで目的間のバランスをとり、そうすることで最適化した選好解を得ることが行われている。単一目的化のためには、主に、各目的関数の線形重み和等の線形の選好関数が提案されている。非線形な選好関数の下でニューラルネットワークを用いた多目的最適化の手法が提案されている（特許文献１参照）。
特開２００２−２８８６２５号（段落［００５５］〜［００６３］、図１） 吉田洋明 外３名、「経路積分の概念に基づいた工学的最適化問題における新たな近似解法」、第５回最適化シンポジウム講演論文集、日本機械学会、２００２年１０月１２日〜１３日、ｐ１８９−１９４ 近藤理良 外３名、「確率過程を用いた最適化問題の近似解法」、第１５回計算力学講演会講演論文集、日本機械学会、２００２年１１月２日〜４日、ｐ４６３−４６４ 吉田洋明 外４名、「工学的問題における確率過程を用いた最適化法」、第１４回計算力学講演会講演論文集、日本機械学会、２００１年１１月２８日〜３０日、ｐ３１５−３１６

目的関数が複数個存在するような場合、従来の手法においては、それら複数個の目的関数に適当な重みを付けて線形結合することにより１つの仮想的な目的関数を作り、その目的関数を最適化するという作業が行われている。しかしながら、この仮想的な目的関数に関する最適化は必ずしもユーザが求めたい解を得ることに対応しているとは限らず、線形結合における重みの設定法によっては、例えば１つの目的関数の値はとても良くなるが、他の目的関数の値はとても悪くなるというような解が得られる可能性もあるという問題点がある。

そこで、従来方法では最適解の探索が困難な最適化問題について、確率過程を用いて、目的関数を最適化する最適解の高精度近似解を求める点、及び、上記の単一目的関数を対象とした確率論的な最適化手法を、互いの間にトレードオフの関係がある複数の目的関数が存在することに起因して最適解が点の集合となるような場合、多目的最適化問題に拡張して適用し、厳密解ではないが、複数の目的関数をそれぞれできるだけ最適化する解として十分実用に耐える近似解を得ることを可能にする点に解決すべき課題がある。

この発明の目的は、経路積分の概念を適用した単目的最適化法及びそのプログラム、及び互いにトレードオフの関係にある複数の目的関数を統合的に最適化を図ることができる多目的最適化法及びそのプログラムを提供することである。

上記の課題を解決するため、この発明による最適化方法は、最適化すべき変数の関数である目的関数について、前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与え、前記確率分布をもって分布する前記変数の多重積分として定義される期待値を前記目的関数の最適解に対する近似解として求めることから成っている。

この最適化方法によれば、最小化すべき目的関数ついて、目的関数の値の関数であって各目的関数の値が小さいほど大きい値を示す確率分布が定められる。確率分布を被積分関数の中に組み込んだ多重積分として期待値が計算されるので、計算機による計算が容易になる。なお、この最適化方法は、動的な単目的における最適化方法に適している。

この発明による最適化プログラムは、最適化すべき変数の関数である目的関数を設定するステップ、前記目的関数の値によって定まり前記各目的関数の値が小さいほど大きい値を示す確率分布を与えるステップ、前記目的関数の最適解に対する近似解を求めるため、前記確率分布をもって分布する前記変数の多重積分として定義される期待値を計算するステップをコンピュータに実行させることから成っている。

コンピュータがこの最適化プログラムを実行するときには、目的関数を設定するステップでは、最適化すべき変数の関数である目的関数を設定し、確率分布を設定するステップでは、その目的関数の値の関数であって目的関数の値が小さいほど大きい値を示す確率分布を定める。期待値を計算するステップでは、その確率分布をもって分布する変数の期待値が多重積分として計算される。計算された期待値は、目的関数を最小化する最適解の十分良い近似解を与える。

この最適化方法及び最適化プログラムにおいて、前記変数をハングライダーのパイロットの操作による主翼の迎角及び前記ハングライダーの機体の構造変数とし、前記目的関数を前記ハングライダーの滑空距離の逆数と飛行拘束条件を破るときに与えられるペナルティとの和とすることができる。機体形状と構造を固定した場合、ハングライダーの最適滑空問題は空気力学を考慮した運動方程式に関連する変分問題として定式化可能である。しかしながら、機体形状と構造にも設計の自由度がある場合に従来のように変分問題として扱うことは、設計可能な機体のすべてに対して最適滑空問題を解く必要があるため、非常に困難である。本発明では、上記のように変数と目的関数とを定めることにより、飛行拘束条件の下で、最も長く滑空することができる迎角の時間経過とその飛行経路及び機体の形状・構造を求めることができる。

また、この発明による最適化方法は、同時に最適化すべき変数の関数である複数の目的関数のそれぞれについて、前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与え、前記複数の目的関数に対応する前記確率分布の積として統合確率分布が表されるように前記複数の目的関数の線形和として統合目的関数を定義し、前記統合確率分布をもって分布する前記変数の多重積分として定義される期待値を、前記統合目的関数の最適解に対する近似解として求めることから成っている。

この最適化方法によれば、同時に最小化すべき複数の目的関数のそれぞれについて、各目的関数の値の関数であって各目的関数の値が小さいほど大きい値を示す確率分布が定められる。複数の目的関数の線形和として統合目的関数を定義するが、統合目的関数に対応する統合確率分布が各目的関数の各確率分布の積として表されるように定義される。その結果、期待値を求める際の多重積分被積分関数の中に組み込まれる確率分布が、各目的関数の確率分布の積の形となるので、計算機による積分計算が容易になる。

また、この発明による最適化プログラムは、同時に最適化すべき変数の関数である複数の目的関数を設定するステップ、前記各目的関数について前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与えるステップ、前記複数の目的関数に対応する前記確率分布の積として統合確率分布が表されるように前記複数の目的関数の線形和として統合目的関数を定義するステップ、及び前記統合目的関数の最適解に対する近似解を求めるため、前記統合確率分布をもって分布する前記変数の多重積分として定義される期待値を計算するステップをコンピュータに実行させることから成っている。

コンピュータがこの最適化プログラムを実行するとき、目的関数を設定するステップでは、最適化すべき変数の関数であり且つ同時に最小化すべき複数の目的関数をそれぞれ設定し、確率分布を設定するステップでは、各目的関数の値の関数であって各目的関数の値が小さいほど大きい値を示す確率分布を定める。統合目的関数を定義するステップでは、複数の目的関数の線形和として統合目的関数が定義される。この際、統合目的関数に対応する統合確率分布が、各目的関数の各確率分布の積として定められる。期待値を計算するステップでは、統合確率分布をもって分布する変数の期待値が多重積分として計算される。計算された期待値は、統合目的関数を最小化する最適解の十分良い近似解を与える。

この最適化方法又は最適化プログラムにおいて、前記変数は構造物を構築するために必要な材料使用量又はそれを誘導する長さ、面積、体積又は重量であり、前記目的関数は前記構造物の変位又は前記構造物が提供する容積に関する量とすることができる。このように変数と目的関数とを定めることにより、一般には別途寸法上の制約がある中で、構造物を構築するために必要な材料使用量を最も少なくしながら、構造物が提供する剛性としての変位又は内部の利用空間としての容積を最大にする最適解の集合を得ることができる。

前記構造物は棒材を組み合わせたトラス又はラーメン構造であり、前記変数は前記棒材の断面積と長さであり、前記目的関数は前記トラス又はラーメン構造に作用する荷重が印加される位置における変位であるとすることができる。このように変数と目的関数とを定めることにより、一般には別途寸法上の制約がある中で、構造物としてのトラス又はラーメン構造を構築するために必要な材料使用量を最も少なくしながら、荷重の作用点における変位量を小さくし、トラス又はラーメン構造の全体体積のような規模を小さくする最適解の集合を得ることができる。

前記構造物は殻構造を有するタンクであり、前記変数は前記タンクに使用する金属薄板材の肉厚と表面積を定める寸法とであり、前記目的関数は前記タンクの体積の逆数であるとすることができる。このように変数と目的関数とを定めることにより、一般には別途寸法上の制約がある中で、タンク内部容量としての最大体積、使用材料量としての最小重量という最適解の集合を得ることができる。

この発明は、上記のように構成されているので、動的な問題を含む単目的最適化問題の場合には、従来方法では局所解に陥るなどのために最適解の探索が困難な最適化問題であっても、確率過程を用いて目的関数を最適化する最適解の高精度近似解を求めることができる。また、多目的最適化問題の場合には、同時に最小化すべき複数の目的関数のそれぞれについて各目的関数の値の関数としての確率分布が定められ、複数の目的関数の線形和として統合目的関数を定義するが、その際に、統合目的関数に対応する統合確率分布が各目的関数の各確率分布の積として表されるように定義する。その結果、期待値を求める際の多重積分被積分関数の中に組み込まれる確率分布が、各目的関数の確率分布の積の形となるので、単一目的における最適化法の場合と同様に、計算機による多重積分計算を行うことができる。この発明によれば、従来方法において、多目的最適化問題における目的関数を１つしか設定できないのに対して、互いにトレードオフの関係にある複数の目的関数を統合的に最適化を図ることができる多目的最適化法及びそのプログラムを提供することができる。変数としては、時間に依存するもの、及び時間に依存しないもののいずれでも対応して統合的に最適化を図ることができる。

以下、添付した図面に基づいて、この発明による単目的・多目的最適化法及びそのためのプログラムの実施例を説明する。この発明による単目的最適化法を、ハングライダーの機体設計と飛行経路の統合的最適化について説明する。図１は、ハングライダーの滑空距離を飛行拘束条件の下で最大にする最適化問題の概要を示す図である。図１に示すように、ハングライダー１が高度１２ｍから滑空を行い、高度２ｍに達するまでの滑空距離Ｌが最大となる機体形状（構造設計を含む）と飛行経路Ｐとを求める最適化問題を考える。この問題の特徴は、ハングライダー１の機体設計という時間に依存しない設計上の最適化と、ハングライダー１の滑空距離Ｌに結びつくパイロットの操縦という時間に依存する変数の最適化との両方を最適化する点にある。飛行経路Ｐは、パイロットの操縦操作によって決まるため、実際にはパイロットの操作で生じる機体の主翼の迎角αa の時間履歴を求めることになる。この最適化問題は、荷重倍数の制限という飛行拘束条件が付いており、そうした状態量拘束条件下での最適化問題である。

機体を質点として扱い、その運動を縦平面内に限ったとき、ハングライダーの運動方程式は、次の式（４）〜式（７）のように表される。

ここで、Ｖ：速度、γ：経路角、ｘ：距離、ｚ：高度、ｍ：質量（機体＋パイロット）、ρ：空気密度、Ｓ^m ：主翼面積、ＣL ：揚力係数 Ｃ’D ：抗力係数、ｇ：重力加速度、ｔ：時間である。また、Ｃ’D は、地面効果を含んだ係数であり、有害抗力係数ＣD 0 と誘導抗力係数Ｃ Di とを用いて式（８）及び式（９）のように表される。

ここで、Ａ^m geは、地面効果を考慮した主翼縦横比、ｂは主翼幅である。
また、境界条件は、以下のように定められる。
初期条件

終端条件
終端高度：ｚf ＝２［ｍ］

ハングライダー１の機体の設計変数が図２に示されている。ハングライダー１の機体は、概略、主翼２と尾翼３から成ると捉えることができ、図２に示す設計変数はそれぞれ、Ｓ^m ：主翼面積、Ａ^m ：主翼アスペクト比、λ^m ：主翼テーパ比、Ｓ^h ：水平尾翼面積、Ａ^h ：水平尾翼アスペクト比、λ^h ：水平尾翼テーパ比、Ｙ：パイプ長、Ｄ：パイプ径、Ｌt ：機体長である。

図３は、翼平面形状のバリエーションを説明する図である。翼の平面形状は、図３に示すように、翼弦（横軸）と翼幅（縦軸）との変化によって変化する。翼の形状は、更に、翼端を切り落とすことにより、翼面積、テーパ比、及びアスペクト比が決定される。したがって、翼の平面形状だけを取り上げても、図３に示す平面形状に加えて、垂直方向（ｚ軸方向）に翼端を切り落とした形状が並んでいる。この三次元空間に分布する翼平面形状から一つの翼平面形状を選択することにより、翼（翼面積、テーパー比、アスペクト比）が決定される。

一つの翼の平面形状が選択されると、その選択された翼の平面形状に対して、桁の設計変数が次のように決定される。主翼２の桁は、図２に示すように、３本のパイプ４，５，６で構成されているので、各パイプ４，５，６の長さ、即ち、主翼２の付け根からパイプ５の先端までの長さＹ1 と、主翼２の付け根からパイプ６の先端までの長さＹ2 の比率を変化させると、図４に示すようなパイプ長さの組合せが存在する。この中から一つの組合せを選択することによって、それぞれのパイプ長さＹ1 ，Ｙ2 が決定される。

最後に、それぞれのパイプ４，５，６の径を変化させると、３本のパイプ径の変化による三次元空間が形成される。図５には、主翼２の付け根側の２本のパイプ５，６の径Ｄ2 とＤ3 とを変化させたときの二次元平面を示してある。実際には、図５に示す平面に垂直な方向にＤ1 の変化による桁の組合せが存在している。このパイプ径の組合せによる三次元区間の中から一つのパイプ径の組合せを選ぶことで、パイプ径が決定される。以上の操作によって、主翼２の設計変数が決定される。しかしながら、機体全体では、更に尾翼３について、水平尾翼の平面形状とパイプ径を決定する必要があるが、水平尾翼については、主翼２の場合と同様であるので詳細については省略する。これらの組合せによって、機体形状の取り得る自由度は大きな次元を持つことが判る。

地上待機中に自重によって機体が破壊しないという条件を満たすこと以外に、滑空中に機体が空気力や慣性力で破壊しないように、翼桁にかかる曲げ応力σと剪断応力τについて、その最大値（σmax 、τmax ）に式（１０）と式（１１）に示すような拘束条件が設定される。

ここでσmax とτmax とは、それぞれ、０．２％引張・圧縮耐力と０．２％剪断耐力であり、１．５は航空機の安全率である。この拘束は滑空中に満たされるべきものであるから、設計した機体が耐えられる飛行を行うように、パイロットの操縦操作、即ち、迎角αa の時間履歴に対してこの拘束がかけられる。なお、空気力によって発生する翼への荷重は揚力線理論によって求められる。なお、揚力線は翼桁に一致するとして、翼には曲げと剪断がのみが作用するとして計算が行われる。

本最適化問題の目的は、滑空距離Ｌが最も長くなるように最適化を図ることである。そこで、式（１２）に示すように、目的関数Ｉとして高度２ｍまで降下するまでの滑空距離ｘf （以下、Ｌの代わりに変数としてｘf を用いる）の逆数を採用することで、最大化問題は最小化問題に変換される。また、滑空中には、曲げ応力σと剪断応力τについて拘束がかけられているので、この拘束を破るような操縦操作には目的関数Ｉにおいて２番目の項に示すペナルティ項を加えることで、評価値が悪化するように図られている。ペナルティ項は、経路上の各的における値の和で表されている。

上記の条件での計算結果が、図６〜図１０に示されている。図６は、最適解、即ち、滑空距離ｘf を最も長くするときのパイロットの操縦操作である迎角αa （度）の飛行距離に対する経過を示すグラフである。図６から、最初の大きな迎角αa から短時間で迎角αa を小さくし、その後徐々に迎角αa を大きくしていることが最適な操縦操作であることが判る。図７は、最適解における飛行距離と速度との関係を示すグラフである。図６とも合わせて考察すると、最適解における飛行態様は、まず迎角αa を減少させて増速することで揚力を大きくし、その後は、滑空を続けると共に速度が次第に低下することに対応して迎角αa を次第に大きくすることで揚力を保つ飛行をしているものと考えられる。図８は、最適解における飛行距離と高度との関係を示すグラフである。最初、迎角αa を急激に小さくして増速を図っているために大きく高度を下げているが、ある程度速度が得られた後は徐々に高度を下げていることが判る。最終的には、３５０ｍ以上の滑空距離に到達していることが判る。

図９は、滑空中における翼桁の曲げ応力と飛行距離との関係を示すグラフである。翼の付け根に当たる第３パイプ６の曲げ応力比が、滑空開始直後の引き起し動作の際に、拘束条件に近づく以外は、拘束条件を満たしていることが判る。したがって、拘束条件を満たすような操縦操作と機体設計とが行われており、最適化が達成されていることが判る。なお、図示しないが、剪断応力については、滑空中に拘束条件を十分満たしており、剪断応力よりも曲げ応力の方が厳しい拘束条件を示すことも判明した。

図１０は、最適化によって得られた機体の形状を示す平面図である。図１０に示すように、最適解を与える機体形状は、アスペクト比が大きな機体であり、常識的に予想されるところの滑空に向いた機体であることが判る。

上記の滑空距離の逆数を目的関数とする本発明による単目的最適化のフローチャートの一例が図１１に示されている。
ステップ１（「Ｓ１」と略す。以下同じ）において、任意の初期機体設計変数ｘ（ｂ，ｏｌｄ）を生成する。なお、ここでｘはベクトル量とする。
ステップ２では、機体設計変数ｘ（ｂ，ｏｌｄ）に対して飛行経路の最適化が行われる。ステップ２における飛行経路の最適化は、図の右側に示す以下の小ステップ（２−１）〜（２−８）に細分化されている。 即ち、
ステップ（２−１）では、任意の初期飛行経路（ここでは、迎角αa のような操縦変数とする）ｘ（ｕ，ｏｌｄ）を生成する。
ステップ（２−２）では、操縦変数ｘ（ｕ，ｏｌｄ）に対する目的関数値Ｉ（ｏｌｄ）を計算する。
ステップ（２−３）では、乱数により新たな操縦変数ｘ（ｕ，ｎｅｗ）を生成する。 ステップ（２−４）では、新たな操縦変数ｘ（ｕ，ｎｅｗ）に対する目的関数値Ｉ（ｎｅｗ）を計算する。
ステップ（２−５）では、メトロポリス法を用いて、目的関数値から定まる確率分布に従うように、ｘ（ｏｌｄ）かｘ（ｎｅｗ）を選択する。具体的には、ΔＩ＝Ｉ（ｎｅｗ）−Ｉ（ｏｌｄ）を計算し、ΔＩ≦０であればｘ（ｕ，ｎｅｗ）を採用するが、ΔＩ＞０であってもｅｘｐ（−ΔＩ／ｈ）の確率でｘ（ｕ，ｎｅｗ）を採用する。それ以外の場合には、ｘ（ｕ，ｏｌｄ）を採用する。そして、採用した操縦変数と目的関数値をそれぞれｘ（ｕ，ｎｅｗ），Ｉ（ｎｅｗ）とする。
ステップ（２−６）では、採用された操縦変数ｘ（ｕ，ｎｅｗ）を使って、操縦変数の期待値＜ｘ（ｕ）＞を計算する。
ステップ（２−７）では、採用された操縦変数ｘ（ｕ，ｎｅｗ）を新たにｘ（ｕ，ｏｌｄ）とする。
ステップ（２−８）では、計算終了回数に達していれば計算を終了し、そうでなければ、ステップ（２−３）に戻る。
ステップ３では、機体設計変数ｘ（ｂ，ｏｌｄ）と、この機体に対して行った飛行経路の最適化によって求めた操縦変数の期待値＜ｘ（ｕ）＞を合わせて、
ｘ（ｏｌｄ）＝｛ｘ（ｂ，ｏｌｄ），＜ｘ（ｕ）＞｝とし、この機体に対する目的関数値Ｉ（ｏｌｄ）を計算する。
ステップ４では、乱数により、新たな機体設計変数ｘ（ｂ，ｎｅｗ）を生成する。
ステップ５では、新たな機体設計変数ｘ（ｂ，ｎｅｗ）に対してステップ２と同様に飛行経路の最適化を行う。
ステップ６では、機体設計変数ｘ（ｂ，ｎｅｗ）と、この機体に対して行った飛行経路の最適化によって求めた操縦変数の期待値＜ｘ（ｕ）＞を合わせて、
ｘ（ｎｅｗ）＝｛ｘ（ｂ，ｎｅｗ），＜ｘ（ｕ）＞｝とし、この機体に対する目的関数値Ｉ（ｎｅｗ）を計算する。
ステップ７では、メトロポリス法を用いて、目的関数値Ｉ（ｏｌｄ）とＩ（ｎｅｗ）から定まる確率分布に従うように、ｘ（ｏｌｄ）かｘ（ｎｅｗ）を選択する。具体的な方法は、ステップ（２−５）と同様であり、採用された変数をｘ（ｎｅｗ）とする。
ステップ８では、採用された変数ｘ（ｎｅｗ）を使って、操縦変数の期待値＜ｘ＞を計算する。
ステップ９では、採用された変数の組ｘ（ｎｅｗ）をｘ（ｏｌｄ）とする。
ステップ１０では、計算終了回数に達していれば計算を終了し、そうでなければ、ステップ４に戻る。計算が終了したとき、期待値＜ｘ＞は、機体設計変数と操縦変数の最適解の近似解となっている。

図１２は、図１１に示すフローチャートで実行される計算のイメージを示す図である。機体１に対して飛行経路１〜Ｋ（ｕ，ｍａｘ）における経路積分が実行され、機体１に対する期待値としての経路が求められる。機体２、・・・機体ｋ、・・・機体Ｋｍａｘについても同様にそれぞれ、飛行経路にしたがって経路積分が実行され、各機体について期待値としての経路が求められる。全機を使って、期待値・機体を算出することで、機体の最適解のよい近似である期待値が求められる。

２つの目的関数を同時に統合化して最小化する多目的最適化問題について説明する。２つの目的関数のそれぞれに、式（２）に対応して、次の式（１４）及び式（１５）に示す確率分布を適用する。一般的な多目的最適化問題では、唯一の最適解が存在することはなく、目的関数同士にトレードオフの関係があるパレート最適解が存在する。

ここでは、２つの目的関数の線形和である統合目的関数Ｉt を、α，βを目的関数Ｉ1 ，Ｉ2 についての重みとして、式（１６）のように定める。

式（１４）と式（１５）におけるｈ1 ，ｈ2 を式（１７）のように定めると、統合確率分布Ｐt は、式（１８）のように求めることができる。式（１８）に定める確率分布を用いて最適化を図ることによって、単一目的から多目的へと拡張が可能になる。

凸関数においては、パレート最適解上のある一点の確率が最も高くなり、最終的には
その確率が最も高い１点が解になるが、計算過程においては、パレート最適解にほぼ一様に解が集まることが判った。

図１３は、この発明を構造物の最適化問題へ適用した場合の一例を示す図である。図１３は、二次元トラス構造の一例を示している。設計変数は各部材の断面積ａi と左端から垂直部材までの距離ｂであり、拘束条件は下記に掲げるとおりであり、目的関数は、荷重Ｗ1 ，Ｗ2 が印加される点の変位の２乗とトラス構造体全体の体積である。また、本題において用いた定数は図１４に示すテーブルに記載のとおりである。

四つのコンセプトにおけるパレート最適解は、図１５に示されているように点の集合として表される。

次に、図１６は、この発明による多目的最適化が適用される圧力容器の一例を示す斜視図である。設計変数は容器の半径Ｒ（０．１≦Ｒ≦３６．０［ｉｎ］），板厚Ｔ（０．５≦Ｔ≦６．０［ｉｎ］）、長さＬｖ（０．１≦Ｌｖ≦１４０．０［ｉｎ］）であって、下記に示すように、拘束条件の下で、目的関数は、容器の容積と重量である。本容器に用いる定数は図１７に示すテーブルに記載のとおりである。また、最適化の結果は、図１８に示すとおりであり、点の集合として表されていることが判る。

ハングライダーの滑空距離が最大となる最適化問題の概要を示す図である。 機体の設計変数を説明する図である。 翼平面形状のバリエーションを説明する図である。 翼桁を構成するパイプの長さの組合せを示す図である。 翼の付け根側の２本のパイプの径を変化させたときの二次元平面を示す図である。 最適解におけるパイロットの操縦操作である仰角の時間履歴を示すグラフである。 最適解における滑空距離と速度との関係を示すグラフである。 最適解における滑空距離と高度との関係を示すグラフである。 滑空中における翼桁の曲げ応力と滑空距離との関係を示すグラフである。 最適解を与える機体形状の平面図である。 ハングライダーの最大滑空距離問題を解くためのフローチャートである。 図１１に示すフローチャートで実行される計算のイメージ図である。 この発明による多目的最適化が適用される構造物の一例を示す構造図である。 図１３に示す構造物の最適化問題において用いた定数のテーブルである。 図１３に示す構造物の最適化問題における四つのコンセプトのパレート最適解を示す図である。 この発明による多目的最適化が適用される圧力容器の一例を示す斜視図である。 図１６に示す圧力容器の最適化問題において用いた定数のテーブルである。 図１６に示す圧力容器の最適化問題におけるパレート最適解を示す図である。 ２点間の一次元の運動について、経路の時間経過の一例を示す図である。 一様重力場中での投げ上げ問題の解析解（太線）と、単一目的最適化においてメトロポリス法により自動生成された経路（細線）とを示す図である。 単一目的最適化における計算の具体的手順を示すフローチャートである。

符号の説明

１ ハングライダー
２ 主翼
３ 尾翼
４，５，６ パイプ
Ｐ 飛行経路
Ｌ 滑空距離
αa 迎角
Ｉ 目的関数

Claims (12)

1. 最適化すべき変数の関数である目的関数について、前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与え、前記確率分布をもって分布する前記変数の多重積分として定義される期待値を前記目的関数の最適解に対する近似解として求めることから成る最適化方法。
2. 前記変数をハングライダーのパイロットの操作による主翼の迎角及び前記グハンライダーの機体の構造変数とし、前記目的関数を前記ハングライダーの滑空距離の逆数と飛行拘束条件を破るときに与えられるペナルティとの和とすることから成る請求項１に記載の最適化方法。
3. 同時に最適化すべき変数の関数である複数の目的関数のそれぞれについて、前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与え、前記複数の目的関数に対応する前記確率分布の積として統合確率分布が表されるように前記複数の目的関数の線形和として統合目的関数を定義し、前記統合確率分布をもって分布する前記変数の多重積分として定義される期待値を、前記統合目的関数の最適解に対する近似解として求めることから成る最適化方法。
4. 前記変数は構造物を構築するために必要な材料使用量又はそれを誘導する長さ、面積、体積又は重量であり、前記目的関数は前記構造物の変位又は前記構造物が提供する容積に関する量であることから成る請求項３に記載の最適化方法。
5. 前記構造物は棒材を組み合わせたトラス又はラーメン構造であり、前記変数は前記棒材の断面積と長さであり、前記目的関数は前記トラス又はラーメン構造に作用する荷重が印加される位置における変位であることから成る請求項４に記載の最適化方法。
6. 前記構造物は殻構造を有するタンクであり、前記変数は前記タンクに使用する金属薄板材の肉厚と表面積を定める寸法とであり、前記目的関数は前記タンクの体積の逆数であることから成る請求項４に記載の最適化方法。
7. 最適化すべき変数の関数である目的関数を設定するステップ、前記目的関数の値によって定まり前記各目的関数の値が小さいほど大きい値を示す確率分布を与えるステップ、前記目的関数の最適解に対する近似解を求めるため、前記確率分布をもって分布する前記変数の多重積分として定義される期待値を計算するステップをコンピュータに実行させる最適化プログラム。
8. 前記変数をハングライダーのパイロットの操作による主翼の迎角及び前記グハンライダーの機体の構造変数とし、前記目的関数を前記ハングライダーの滑空距離の逆数と飛行拘束条件を破るときに与えられるペナルティとの和とすることから成る請求項７に記載の最適化プログラム。
9. 同時に最適化すべき変数の関数である複数の目的関数を設定するステップ、前記各目的関数について前記目的関数の値によって定まり前記目的関数の値が小さいほど大きい値を示す確率分布を与えるステップ、前記複数の目的関数に対応する前記確率分布の積として統合確率分布が表されるように前記複数の目的関数の線形和として統合目的関数を定義するステップ、及び前記統合目的関数の最適解に対する近似解を求めるため、前記統合確率分布をもって分布する前記変数の多重積分として定義される期待値を計算するステップをコンピュータに実行させる最適化プログラム。
10. 前記変数は構造物を構築するために必要な材料使用量又はそれを誘導する長さ、面積、体積又は重量であり、前記目的関数は前記構造物の変位又は前記構造物が提供する容積に関する量であることから成る請求項９に記載の最適化プログラム。
11. 前記構造物は棒材を組み合わせたトラス又はラーメン構造であり、前記変数は前記棒材の断面積と長さであり、前記目的関数は前記トラス又はラーメン構造に作用する荷重が印加される位置における変位であることから成る請求項１０に記載の最適化プログラム。
12. 前記構造物は殻構造を有するタンクであり、前記変数は前記タンクに使用する金属薄板材の肉厚と表面積を定める寸法とであり、前記目的関数は前記タンクの体積の逆数であることから成る請求項１０に記載の最適化プログラム。

Images