JP5019744B2 - 多目的最適化装置、多目的最適化方法および多目的最適化プログラム - Google Patents

多目的最適化装置、多目的最適化方法および多目的最適化プログラム Download PDF

Info

Publication number
JP5019744B2
JP5019744B2 JP2005369267A JP2005369267A JP5019744B2 JP 5019744 B2 JP5019744 B2 JP 5019744B2 JP 2005369267 A JP2005369267 A JP 2005369267A JP 2005369267 A JP2005369267 A JP 2005369267A JP 5019744 B2 JP5019744 B2 JP 5019744B2
Authority
JP
Japan
Prior art keywords
individual
evaluation
fitness
individuals
storage
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2005369267A
Other languages
English (en)
Other versions
JP2007172306A (ja
Inventor
洋隆 梶
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yamaha Motor Co Ltd
Original Assignee
Yamaha Motor Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Yamaha Motor Co Ltd filed Critical Yamaha Motor Co Ltd
Priority to JP2005369267A priority Critical patent/JP5019744B2/ja
Publication of JP2007172306A publication Critical patent/JP2007172306A/ja
Application granted granted Critical
Publication of JP5019744B2 publication Critical patent/JP5019744B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1401—Introducing closed-loop corrections characterised by the control or regulation method
    • F02D41/1406—Introducing closed-loop corrections characterised by the control or regulation method with use of a optimisation method, e.g. iteration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02D—CONTROLLING COMBUSTION ENGINES
    • F02D41/00—Electrical control of supply of combustible mixture or its constituents
    • F02D41/02—Circuit arrangements for generating control signals
    • F02D41/14—Introducing closed-loop corrections
    • F02D41/1401—Introducing closed-loop corrections characterised by the control or regulation method
    • F02D2041/1433—Introducing closed-loop corrections characterised by the control or regulation method using a model or simulation of the system
    • F02D2041/1437—Simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • Genetics & Genomics (AREA)
  • Computational Linguistics (AREA)
  • Mechanical Engineering (AREA)
  • Computer Hardware Design (AREA)
  • Combustion & Propulsion (AREA)
  • Physiology (AREA)
  • Chemical & Material Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、最適化対象のパラメータを最適化する多目的最適化装置、多目的最適化方法および多目的最適化プログラムに関する。
本明細書において、最適化とは、最適解または最適に近い解（準最適解）を求めることをいう。
従来、多目的最適化問題と呼ばれる問題クラスが存在する。例えば、ある製品のコストを最小化し、性能を最大化するという問題を考えた場合、これは２つの目的関数の多目的最適化問題となる。この場合、コストおよび性能が２つの目的関数となる。一般的には、コストを下げると性能が悪化し、性能を上げるとコストがかさむというトレードオフの関係が生じるために、多目的最適化問題の解は一つではない。
図３７は多目的最適化問題をエンジンの最適化に適用した例を示す図である。多目的最適化問題をエンジンの燃費およびトルクの最適化に適用する場合、燃費およびトルクが２つの目的関数ｆ1，ｆ2である。この場合、燃料噴射量、点火時期等のパラメータを調整することにより目的関数ｆ1，ｆ2の値を最適化する。
解Ａは、燃費が解Ｂに比べて優れているが、トルクが解Ｂに比べて劣っている。このように、エンジンの燃費とエンジンのトルクとはトレードオフの関係を有するため、複数の最適解が存在する。使用者は、複数の最適解から目的に合った解を選択することができる。例えば、スポーツ走行に適した自動二輪車に用いるエンジンには解Ａを選択し、ロングツーリングに適した自動二輪車に用いるエンジンには解Ｂを選択する。
一般に多目的最適化問題は、Ｎ個のパラメータについてＭ個の目的関数の値を、各パラメータの制約条件の範囲で最小化する問題と定義される。目的関数の値を最大化する場合は、目的関数に負の符号を付けて目的関数の値を最小化する問題に変換することとする。
このような多目的最適化問題は、一般的に単一の最適解を持たず、パレート最適解と呼ばれる概念で定義される最適解集合を持つ。ここで、パレート最適解とは、ある目的関数の値を改善するためには、少なくとも１つの他の目的関数の値を改悪せざるを得ない解のことをいい、以下のように定義される（例えば、非特許文献１参照）。
〔定義１〕あるｐ個の目的関数ｆk（ｋ＝１，・・・，ｐ）の２つの解ｘ１，ｘ２∈Ｆに関して、ｆk（ｘ１）≦ｆk（ｘ２）（∀ｋ＝１，・・・，ｐ）∧ｆk（ｘ１）＜ｆk（ｘ２）（∃ｋ＝１，・・・，ｐ）のとき、ｘ１はｘ２に優越するという。ここで、Ｆは解の集合である。
〔定義２〕ある解ｘ０が他のすべての解ｘ∈Ｆに優越するとき、ｘ０は最適解である。
〔定義３〕ある解ｘ０に優越する解ｘ∈Ｆが存在しないとき、ｘ０はパレート最適解（または非劣解）である。
パレート最適解集合を求めることは、目的関数のトレードオフに関して最適な解の集合を求めることになる。
図３８はパレート最適解について説明するための図である。図３８は２つの目的関数ｆ1，ｆ2の例を示す。解ａについての目的関数ｆ1の値ｆ1（ａ）は解ｂについての目的関数ｆ1の値ｆ1（ｂ）よりも小さく、解ａについての目的関数ｆ2 の値ｆ2（ａ）は解ｂについての目的関数ｆ2の値ｆ2（ｂ）よりも小さい。したがって、解ａは解ｂに優越する。
同様に、解ａは解ｃ，ｄに優越する。解ａに優越する解は存在しない。同様に、解ｅ，ｆに優越する解も存在しない。したがって、解ａ，ｅ，ｆはパレート最適解である。
なお、解ｇは、弱パレート最適解である。弱パレート最適解とは、ある目的関数についてのみパレート最適解に優越されないパレート解である。弱パレート最適解は、合理的な解ではなく、本来求める必要のない解である。
多目的最適化問題の解法は多数提案されている。最近注目されている方法に多目的進化型アルゴリズム（ＭＯＥＡｓ：Multiple Objective Evolutionary Algorithm）がある。
この方法の最大の特徴は、進化型アルゴリズムの多点探索を利用してパレート最適解集合を一度に求めることである。得られたパレ一ト最適解集合は、その中から目的に合致した解を探す意志決定、またはパレ一ト最適解集合（パレート境界）の形状からの知見の獲得等に用いられる。
進化型アルゴリズムとして遺伝的アルゴリズム（ＧＡ：Genetic Algorithm）を多目的最適化問題に適用する研究が数多く行われている。遺伝的アルゴリズムは、生物の適応進化を模倣した計算手法である。遺伝的アルゴリズムでは、解の候補を個体と呼ぶ。また、目的関数は適応関数と呼び、適応度関数の値を適応度と呼ぶ。
この遺伝的アルゴリズムは、自然進化に見られる過程（染色体の選択、交叉および突然変異）をヒントにして、Ｊ．Ｈｏｌｌａｎｄにより提案されたアルゴリズムである。設計変数を遺伝子とみなして、初期設計の個体集合をランダムに生成し、各個体の適応度を評価する。適応度の良い個体ほど親として選択される可能性が高くなるように親を選択する。そして、交叉（遺伝子の入れ換え）および突然変異（遺伝子のランダムな変化）により子孫を作る。さらに、評価、選択、交叉および突然変異により世代を繰り返し、最適解を探索する。
具体的には、ＦｏｎｓｅｃａらのＭＯＧＡ（Multiple Objective Genetic Algorithm：例えば、非特許文献１参照）、ＤｅｂらのＮＳＧＡ−ＩＩ（Non-Dominated Sorting Genetic Algorithm-II：例えば、非特許文献２参照）、ＺｉｔｚｌｅｒらのＳＰＥＡ２（Strength Pareto Evolutionary Algorithm ２： 例えば、非特許文献３参照）等が提案されている。特に、ＮＳＧＡ−ＩＩおよびＳＰＥＡ２は優秀な多目的進化型アルゴリズムとして知られている。
また、不確実性を伴う適応度関数の最適化のために、探索履歴を用いて真の適応度を推定するＭＦＥＧＡ（Memory-based Fitness Estimation Genetic Algorithm）が佐野、喜多らにより提案されている（例えば、非特許文献４参照）。ここで、ＭＦＥＧＡでは、過去に得られた個体の適応度のサンプル値を探索履歴として保存し、探索履歴を参照して真の適応度を推定する。ＭＦＥＧＡは、不確実性を持つ問題に関して、同一個体の適応度を複数回サンプリングする方法および通常の遺伝的アルゴリズムに比べて優秀な探索性能を持つことが報告されている。
C.M.fonseca，p.J.Flemimg:genetic algorithms for multiobjective optimization：formulation,discussion and generalization，of the 5th international conference on genetic algorithms，pp.416-423(1993) K.Deb，S.Agrawal，A.Pratab，and T.Meyarivan:A Fast Elitist Non-Dominated Sorting Genetic Algorithm for Multi-Objective Optimization:NSGA-II，KanGAL report 20001，Indian Institute of Technology，Kanpur，India(20OO) E.Zitzler,M.Laumanns,L.Thiele:SPEA2:Improving the Performance of the Strength Pareto Evolutionary A1gorithm,Technical Report 103,Computer Engineering and Communication Networks Lab(TIK),Swiss Federal Institute of Technology(ETH)Zurich(2001) 佐野，喜多：探索履歴を利用した遺伝的アルゴリズムによる不確実関数の最適化，電学論Ｃ １２２巻６号，ＰＰ−１００１−１００８（２００２） K.Ikeda, H.Kita, and S.Kobayashi : Failuer of Pareto-Based MOEAs, Does Non-Dominated Really Mean Near to Optimal？ Congress on Evolutionary Computation, pp.957-962(2001) M.D.Berg, et.al. : Computational Geometry : Algorithms and Applications, Springer-Verlag (1997) 今井浩、今井桂子，計算幾何学，情報数学講座１２，共立出版（１９９４） E.Zitzler, K.Deb, L,Thiele : Comparison of Mu1tiobjective Evo1utionary Algorithms : Empirical Results, Evolutionary Computation 8(2), pp.173-195 (2000) Y. Jin, M. Olhofer, B. Sendhoff: On Evolutionary Optimization with Approximate Fitness Functions, Proc. of GECCO2000, pp. 786-792 (2000)
しかしながら、通常の単目的遺伝的アルゴリズムおよび多目的遺伝的アルゴリズムを用いて大きな計算コストを有するシミュレーションまたは実験ベースの最適化を行う場合、適切なパレート最適個体を得るために適応度の推定値の評価回数が非常に多くなり、膨大な時間を要する。そこで、種々の評価回数低減手法が提案されている。
例えば、非特許文献９には、多数の子個体を生成し、それらの子個体の適応度をニューラルネットワークにより推定し、適応度の推定値が高い子個体のみを実環境で評価するという単目的の進化型計算法が提案されている。
一方、適応度の推定値の評価にランキングを用いた多目的遺伝的アルゴリズムにおいては、実環境で評価可能な個体数以上のパレート最適個体が子個体に含まれている場合、実環境で評価すべき子個体をどのようにして選択するかが問題となる。しかしながら、このような場合において子個体の有効な選択方法は提案されていない。
したがって、従来の多目的遺伝的アルゴリズムを用いてパレート最適個体を短時間で効率良く得ることができない。
本発明の目的は、適切なパレート最適個体を短時間で効率良く得ることが可能な多目的最適化装置、多目的最適化方法および多目的最適化プログラムを提供することである。
（１）第１の発明に係る多目的最適化装置は、最適化対象に個体のパラメータの組を与え、複数の目的に対応する複数の適応度関数についての適応度の組を最適化対象から受ける多目的最適化装置であって、個体のパラメータの組および最適化対象から出力される適応度の組を保存用個体として記憶する記憶部と、新たな個体を評価用個体候補として生成し、生成された評価用個体候補のうち評価用個体として選択された個体のパラメータの組を最適化対象および記憶部に与えるとともに、複数組の適応度に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める演算部と、記憶部に記憶された複数の保存用個体に対応する複数組の適応度に基づいて、演算部により生成された評価用個体候補および記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める推定部とを備え、演算部は、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣を比較し、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択するものである。
その多目的最適化装置においては、個体のパラメータの組および最適化対象から出力される適応度の組が保存用個体として記憶部に記憶される。演算部により新たな個体が評価用個体候補として生成され、生成された評価用個体候補のうち評価用個体として推定部により選択された個体のパラメータの組が最適化対象および記憶部に与えられる。また、複数組の適応度に基づいて評価用個体集合が多目的進化型アルゴリズムに従って演算部により評価される。それにより、パレート最適個体集合が求められる。
上記の評価用個体の選択の際には、記憶部に記憶された複数の保存用個体に対応する複数組の適応度に基づいて、評価用個体候補および保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組が推定部により求められる。そして、演算部により複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣が比較され、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けが行われる。次いで、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標が算出され、算出された分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択される。
このように、複数の評価用個体候補が生成され、最上位ランクの複数の評価用個体候補から分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択され、選択された評価用個体が実際に評価される。それにより、不要な評価用個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
（２）第２の発明に係る多目的最適化装置は、最適化対象に個体のパラメータの組を与え、複数の目的に対応する複数の適応度関数についての適応度のサンプル値の組を最適化対象から受ける多目的最適化装置であって、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶する記憶部と、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求める第１の推定部と、第１の推定部により求められた推定値に基づいて新たな個体を評価用個体候補として生成し、生成された評価用個体候補のうち評価用個体として選択された個体のパラメータの組を最適化対象および記憶部に与えるとともに、第１の推定部により求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める演算部と、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、演算部により生成された評価用個体候補および記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める第２の推定部とを備え、演算部は、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣を比較し、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択するものである。
その多目的最適化装置においては、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組が保存用個体として記憶部に記憶される。記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組が第１の推定部により求められる。推定値に基づいて演算部により新たな個体が評価用個体候補として生成され、生成された評価用個体候補のうち評価用個体として第２の推定部により選択された個体のパラメータの組が最適化対象および記憶部に与えられる。また、求められた複数組の推定値に基づいて評価用個体集合が多目的進化型アルゴリズムに従って演算部により評価される。それにより、パレート最適個体集合が求められる。
上記の評価用個体の選択の際には、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、評価用個体候補および保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組が第２の推定部により求められる。そして、演算部により複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣が比較され、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けが行われる。次いで、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標が算出され、算出された分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択される。
このように、複数の評価用個体候補が生成され、最上位ランクの複数の評価用個体候補から分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択され、選択された評価用個体が実際に評価される。それにより、不要な評価用個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
（３）第１の推定部は、記憶部に記憶された複数の保存用個体をｈlとし、注目個体ｘに対応するサンプル値の組をＦ（ｘ）とし、パラメータ空間上で注目個体から距離ｄl離れた個体に対応するサンプル値の組をＦ（ｈl）とし、ｋ’を係数とし、ｌ＝１，…，Ｈとし、ｎを自然数とした場合に、
Figure 0005019744
で表される推定式により注目個体ｘに対応する真の適応度の推定値の組ｆ’（ｘ）を算出してもよい。
この場合、パラメータ空間上で注目個体と他の保存用個体との距離を考慮して真の適応度からの誤差がより小さい推定値を得ることができる。
（４）第２の推定部は、記憶部に記憶された複数の保存用個体をｈlとし、パラメータ空間上で評価用個体候補から距離ｄl離れた保存用個体に対応するサンプル値の組をＦ（ｈl）とし、ｋ’を係数とし、ｌ＝１，…，Ｈとし、ｎを自然数とした場合に、
Figure 0005019744
で表される推定式により評価用個体候補ｘに対応する真の適応度の推定値の組ｆ’（ｘ）を算出してもよい。
この場合、パラメータ空間上で評価用個体候補と保存用個体との距離を考慮して真の適応度からの誤差がより小さい推定値を得ることができる。
（５）ｎは１であってもよい。この場合、推定値の算出において、パラメータ空間上で評価用個体候補から遠く離れた保存用個体の影響が小さくなる。それにより、真の適応度からの誤差が十分に小さい推定値を得ることができる。
（６）ｎは３であってもよい。この場合、推定値の算出において、パラメータ空間上で評価用個体候補から遠く離れた保存用個体の影響が大幅に小さくなる。それにより、真の適応度からの誤差がさらに十分に小さい推定値を得ることができる。
（７）演算部は、適応度関数空間上で最上位ランク内で各評価用個体候補に隣接する２つの保存用個体を検出し、適応度関数ごとに前記２つの保存用個体間の数学的距離を算出し、各評価用個体候補についての複数の適応度関数における数学的距離の合計を前記分布指標として算出し、複数の評価用個体候補から算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択してもよい。
この場合、適応度関数空間上で最上位ランク内で保存用個体が疎らな領域に新たな評価用個体を容易に生成することが可能となる。それにより、多様性の高いパレート最適個体を容易に得ることができる。
（８）演算部は、第１の推定部により求められた推定値の組に基づいてパレート最適個体を表示してもよい。
この場合、使用者は、パレート最適個体を視覚的に認識することができるので、種々の意思決定を容易に行うことができる。
（９）演算部は、多目的進化型アルゴリズムとして遺伝的アルゴリズムを用いて評価用個体集合の個体を評価してもよい。
この場合、遺伝的アルゴリズムに基づいて世代交代を行うことにより、最適なパレート最適個体を容易に得ることができる。
（１０）最適化対象は、機器の複数の性能を評価するための評価システムを含み、パラメータの組は、評価システムのための制御用パラメータの組を含み、複数の適応度関数は評価システムの評価により得られる複数の性能であり、適応度の組は複数の性能の値であってもよい。
この場合、評価システムに制御用パラメータの組が与えられる。評価システムにより制御用パラメータの組に基づいて機器の性能が評価され、複数の性能に対応するサンプル値の組が出力される。この多目的最適化装置によれば、適切な制御用パラメータの組の集合をパレート最適個体として短時間で効率良く得ることができる。
（１１）機器はエンジンであってもよい。この場合、評価システムにエンジン制御用パラメータの組が与えられる。評価システムによりエンジン制御用パラメータの組に基づいてエンジンの性能が評価され、複数の性能に対応するサンプル値の組が出力される。この多目的最適化装置によれば、適切なエンジン制御用パラメータの組の集合をパレート最適個体として短時間で効率良く得ることができる。
（１２）機器はモータであってもよい。この場合、評価システムにモータ制御用パラメータの組が与えられる。評価システムによりモータ制御用パラメータの組に基づいてモータの性能が評価され、複数の性能に対応するサンプル値の組が出力される。この多目的最適化装置によれば、適切なモータ制御用パラメータの組の集合をパレート最適個体として短時間で効率良く得ることができる。
（１３）評価システムは、パラメータの組に基づいて機器を制御するとともに機器の動作により発生される複数の性能の値をサンプル値として出力する機器評価装置であってもよい。
この場合、機器評価装置に制御用パラメータの組が与えられる。機器評価装置により制御用パラメータの組に基づいて機器の性能が評価され、複数の性能に対応するサンプル値の組が出力される。この多目的最適化装置によれば、適切な制御用パラメータの組の集合をパレート最適個体として効率良く短時間で得ることができる。
（１４）評価システムは、パラメータの組に基づいて機器の動作をシミュレーションすることにより複数の性能を評価し、評価された複数の性能の値をサンプル値の組として出力する機器シミュレータであってもよい。
この場合、機器シミュレータに制御用パラメータの組が与えられる。機器シミュレータにより制御用パラメータの組に基づいて機器の動作がシミュレーションされることにより複数の性能が評価され、評価された複数の性能に対応するサンプル値の組が出力される。この多目的最適化装置によれば、適切な制御用パラメータの組の集合をパレート最適個体として短時間で効率良く得ることができる。
（１５）第３の発明に係る多目的最適化方法は、最適化対象に個体のパラメータの組を与え、最適化対象から出力される複数の目的に対応する複数の適応度関数についての適応度のサンプル値の組に基づいてパラメータを最適化する多目的最適化方法であって、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶部に記憶するステップと、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求めるステップと、求められた推定値に基づいて新たな個体を評価用個体候補として生成するステップと、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、生成された評価用個体候補および記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求めるステップと、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣を比較し、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択するステップと、選択された評価用個体のパラメータの組を最適化対象および記憶部に与えるとともに、求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求めるステップとを含むものである。
その多目的最適化方法においては、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組が保存用個体として記憶部に記憶される。記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組が求められる。
次に、推定値に基づいて新たな個体が評価用個体候補として生成される。さらに、複数の保存用個体に対応する複数組のサンプル値に基づいて、評価用個体候補および保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組が求められる。そして、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣が比較され、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けが行われる。次いで、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標が算出され、算出された分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択される。
さらに、評価用個体として選択された個体のパラメータの組が最適化対象および記憶部に与えられる。また、求められた複数組の推定値に基づいて評価用個体集合が多目的進化型アルゴリズムに従って評価される。それにより、パレート最適個体集合が求められる。
このように、複数の評価用個体候補が生成され、最上位ランクの複数の評価用個体候補から分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択され、選択された評価用個体が実際に評価される。それにより、不要な評価用個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
（１６）第４の発明に係る多目的最適化プログラムは、最適化対象に個体のパラメータの組を与え、最適化対象から出力される対応する複数の適応度関数についての適応度のサンプル値の組に基づいてパラメータを最適化するコンピュータにより実行可能な多目的最適化プログラムであって、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶部に記憶する処理と、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求める処理と、求められた推定値に基づいて新たな個体を評価用個体候補として生成する処理と、記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、生成された評価用個体候補および記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める処理と、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣を比較し、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択する処理と、選択された評価用個体のパラメータの組を最適化対象および記憶部に与えるとともに、求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める処理とを、コンピュータに実行させるものである。
その多目的最適化プログラムにおいては、個体のパラメータの組および最適化対象から出力される適応度のサンプル値の組が保存用個体として記憶部に記憶される。記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組が求められる。
次に、推定値に基づいて新たな個体が評価用個体候補として生成される。さらに、複数の保存用個体に対応する複数組のサンプル値に基づいて、評価用個体候補および保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組が求められる。そして、複数の適応度関数の各々について和集合に含まれる複数の個体に対応する推定値の優劣が比較され、複数の適応度関数の各々について優劣の比較結果に基づいて和集合に含まれる複数の個体のランク付けが行われる。次いで、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標が算出され、算出された分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択される。
さらに、評価用個体として選択された個体のパラメータの組が最適化対象および記憶部に与えられる。また、求められた複数組の推定値に基づいて評価用個体集合が多目的進化型アルゴリズムに従って評価される。それにより、パレート最適個体集合が求められる。
このように、複数の評価用個体候補が生成され、最上位ランクの複数の評価用個体候補から分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択され、選択された評価用個体が実際に評価される。それにより、不要な評価用個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
本発明によれば、複数の評価用個体候補が生成され、最上位ランクの複数の評価用個体候補から分布指標が大きい順に所定数の評価用個体候補が評価用個体として選択され、選択された評価用個体が実際に評価される。それにより、不要な評価用個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
（１）第１の実施の形態
まず、本発明の第１の実施の形態に係る多目的最適化装置を図１に基づき説明する。
（ａ）多目的最適化装置の機能的な構成
図１は本発明の第１の実施の形態に係る多目的最適化装置の機能的な構成を示すブロック図である。
図１の多目的最適化装置１は、多目的進化型アルゴリズムとして多目的遺伝的アルゴリズム（ＧＡ）を利用して多目的最適化問題のパレート最適個体集合を算出する。この多目的最適化装置１は、最適化対象６に接続される。
最適化対象６は、機器の性能を評価する評価システムである。評価システムは、実システムを評価する評価装置または確率要素を含むシミュレータである。実システムは、例えばエンジンまたはモータであり、評価装置は、例えばエンジン評価装置またはモータ評価装置である。また、シミュレータは、例えばエンジンシミュレータまたはモータシミュレータである。本実施の形態では、最適化対象６はエンジン評価装置である。
多目的最適化装置１は、多目的進化型アルゴリズム部２、適応度推定部３、出力インタフェース４および入力インタフェース５を含む。適応度推定部３は、適応度推定モジュール３０ａ，３０ｂおよび探索履歴記憶装置３１を含む。
多目的進化型アルゴリズム部２および適応度推定部３の適応度推定モジュール３０ａ，３０ｂは、後述するＣＰＵ１０１（図２）が多目的最適化プログラムを実行することにより実現される。探索履歴記憶装置３１は、後述する外部記憶装置１０６（図２）により構成される。
使用者１０は、多目的進化型アルゴリズム部２に複数の適応度関数（目的関数）を設定する。本実施の形態では、複数の適応度関数としては、燃費、トルク、エンジンの排気ガスに含まれるＣＯ（一酸化炭素）、ＨＣ（炭化水素）、ＮＯx （窒素酸化物）等の成分の濃度等のうち複数が設定される。
ここで、トレードオフの関係としては、トルクと燃費、トルクとＣＯ濃度、トルクとＨＣ濃度、燃費とＮＯx 濃度、ＣＯ濃度とＮＯx 濃度、ＨＣ濃度とＮＯx 濃度等が挙げられる。
また、多目的遺伝的アルゴリズムの個体とは、多目的最適化問題の解の候補であり、複数のパラメータの組および複数の適応度を有する。パラメータは、調整可能な値であり、遺伝的アルゴリズムでは、遺伝子と呼ばれる。パラメータとしては、燃料噴射量、燃料噴射時期、点火時期、スロットル開度等が挙げられる。適応度は、適応度関数の値である。以下、多目的遺伝的アルゴリズムの個体を単に個体と呼ぶ。
多目的進化型アルゴリズム部２は、後述する適応度推定モジュール３０ａ，３０ｂにより算出された真の適応度の推定値を受け、個体のパラメータの組を適応度推定部３の探索履歴記憶装置３１に与えるとともに出力インタフェース４を介して最適化対象６に与える。以下、真の適応度の推定値を単に適応度の推定値または推定値と呼ぶ。
最適化対象６は、多目的最適化装置１から与えられた個体のパラメータの組に基づいて適応度のサンプル値の組を出力する。最適化対象６から出力される各サンプル値は、後述するように、真の適応度およびノイズ成分を含む。サンプル値の詳細については、後述する。
最適化対象６から出力されるサンプル値の組は、入力インタフェース５を介して適応度推定部３の探索履歴記憶装置３１に入力される。探索履歴記憶装置３１は、個体のパラメータの組およびサンプル値の組を探索履歴として記憶する。以下、探索履歴に含まれる個体のパラメータおよびサンプル値の各組を履歴データと呼ぶ。
適応度推定モジュール３０ａ，３０ｂは、探索履歴記憶装置３１に記憶される探索履歴の履歴データに基づいて適応度の推定値の組を算出し、推定値の組を多目的進化型アルゴリズム部２に与える。特に、適応度推定モジュール３０ｂは、後述する子個体の選択時に子個体候補の適応度の推定値を算出するために用いられる。詳細は後述する。適応度推定モジュール３０ａ，３０ｂが１つの適応度推定モジュールにより構成されてもよい。
多目的進化型アルゴリズム部２は、複数組の推定値に基づいて遺伝的アルゴリズムにしたがって複数の個体を発生して多点探索を行い、適応度関数をパレート最適性で評価することにより、パレート最適個体集合を同時に求める。また、多目的進化型アルゴリズム部２は、求められたパレート最適個体集合を使用者１０に提示する。
このように、多目的進化型アルゴリズム部２および適応度推定部３は、協働することにより個体のパラメータの最適化を行う。
多目的進化型アルゴリズム部２および適応度推定部３の詳細な動作については後述する。
（ｂ）多目的最適化装置のハードウエア構成
図２は図１の多目的最適化装置１のハードウエア構成を示すブロック図である。
多目的最適化装置１は、ＣＰＵ（中央演算処理装置）１０１、ＲＯＭ（リードオンリメモリ）１０２、ＲＡＭ（ランダムアクセスメモリ）１０３、入力装置１０４、表示装置１０５、外部記憶装置１０６、記録媒体駆動装置１０７および入出力インタフェース１０８を含む。
入力装置１０４は、キーボード、マウス等からなり、各種指令および各種データを入力するために用いられる。ＲＯＭ１０２にはシステムプログラムが記憶される。記録媒体駆動装置１０７は、ＣＤ（コンパクトディスク）ドライブ、ＤＶＤ（デジタルバーサタイルディスク）ドライブ、フレキシブルディスクドライブ等からなり、ＣＤ、ＤＶＤ、フレキシブルディスク等の記録媒体１０９に対してデータの読み書きを行う。
記録媒体１０９には、多目的最適化プログラムが記録されている。外部記憶装置１０６は、ハードディスク装置等からなり、記録媒体駆動装置１０７を介して記録媒体１０９から読み込まれた多目的最適化プログラムおよび各種データを記憶する。ＣＰＵ１０１は、外部記憶装置１０６に記憶された多目的最適化プログラムをＲＡＭ１０３上で実行する。
表示装置１０５は、液晶表示パネル、ＣＲＴ（陰極線管）等からなり、パレート最適個体集合等の各種画像を表示する。入出力インタフェース１０８は、図１の出力インタフェース４および入力インタフェース５を含む。この入出力インタフェース１０８には最適化対象６が無線通信または有線通信により接続される。入出力インタフェース１０８は、最適化対象６から出力されるサンプル値の組を外部記憶装置１０６に転送するとともに、多目的最適化プログラムにより生成された個体のパラメータの組を最適化対象６に与える。
なお、多目的最適化プログラムを記録する記録媒体１０９として、ＲＯＭ等の半導体メモリ、ハードディスク等の種々の記録媒体を用いることができる。また、多目的最適化プログラムを通信回線等の通信媒体を介して外部記憶装置１０６にダウンロードし、ＲＡＭ１０３上で実行してもよい。
ここで、記録媒体１０９は、コンピュータで読み取り可能な記録媒体であれば、電子的読み取り方式、磁気的読み取り方式、光学的読み取り方式またはその他のあらゆる読み取り方式の記録媒体を含むものである。例えば、上記のＣＤ、ＤＶＤおよびフレキシブルディスクの他、ＣＤＶ（コンパクトディスクビデオ）等の光学的読取方式記録媒体、ＲＡＭ、ＲＯＭ等の半導体記録媒体、ハードディスク等の磁気記録型記録媒体、ＭＯ（光磁気ディスク）等の磁気記憶型／光学的読取方式記録媒体を用いることができる。
（ｃ）最適化対象の構成
図３は最適化対象６の構成の一例を示すブロック図である。図３の最適化対象６はエンジン評価装置である。
最適化対象６は、エンジン６１、ＥＣＵ（エンジン制御ユニット）６２、排気ガス分析計６３、コントローラ６４、スロットルユニット６５およびダイナモ６６を含む。
ＥＣＵ６２は、多目的最適化装置１からシリアル通信によりパラメータの組を受信する。本例では、パラメータの組は、点火時期および燃料噴射時期である。ＥＣＵ６２は、パラメータの組に基づいてエンジン６１の点火時期および燃料噴射時期を制御する。エンジン６１からコントローラ６４に回転数、空燃比等のエンジン情報が与えられる。
コントローラ６４は、エンジン情報に基づいてスロットルユニット６５およびダイナモ６６を制御する。スロットルユニット６５は、エンジン６１の吸入空気量を調整することによりエンジン６１の出力トルクを制御する。ダイナモ６６は、負荷トルクを制御する。
排気ガス分析計６３は、エンジン６１からの排気ガス中の成分を分析し、ＨＣ濃度およびＮＯx濃度をサンプル値の組としてシリアル通信により多目的最適化装置１に出力する。
図４はＨＣ濃度、ＮＯx濃度およびＣＯ濃度と空燃比との関係を示す図である。
燃料が完全に燃焼すれば、排気ガスに二酸化炭素と水とが含まれる。しかし、運転状態が変化すると燃焼状態も変化し、排気ガスにＣＯ、ＨＣおよびＮＯxが含まれる。
図４に示すように、空燃比が小さいと、ＨＣ濃度およびＮＯx 濃度が低くなる。ＮＯx 濃度は空燃比が理論空燃比（１４．７）よりやや小さいときに最大となり、それ以外の領域では減少する。理論空燃比付近では、ＨＣ濃度とＮＯx 濃度はトレードオフの関係を有し、ＣＯ濃度とＮＯx 濃度とはトレードオフの関係を有する。
図３の最適化対象６は、多目的最適化装置１から個体のパラメータの組として点火時期を示す点火時期情報および燃料噴射時期を示す燃料噴射時期情報を受け、サンプル値の組としてＨＣ濃度およびＮＯx 濃度を出力する。
（ｄ）多目的最適化装置の全体処理
図５および図６は図１の多目的最適化装置１の全体処理を示すフローチャートである。
図５に示すように、最適化処理が開始されると、多目的進化型アルゴリズム部２は、初期個体集合として親個体集合Ｐを各パラメータの定められた範囲内でランダムに生成することにより親個体集合Ｐを初期化し、生成された親個体集合Ｐの各個体のパラメータの組を最適化対象６に順次与える（ステップＳ１）。それにより、最適化対象６から適応度のサンプル値の組が順次出力される。
なお、事前知識として知られているパレート最適個体が存在する場合は、そのパレート最適個体を初期個体集合の一部として用いてもよい。それにより、パレート最適個体の探索の収束性が向上する。
適応度推定部３の探索履歴記憶装置３１は、最適化対象６から親個体集合Ｐの各個体に対応するサンプル値の組を取得し、親個体集合Ｐのパラメータの組およびサンプル値の組を探索履歴として記憶する（ステップＳ２）。
次に、適応度推定部３の適応度推定モジュール３０ａは、探索履歴記憶装置３１に記憶された複数の個体に対応するサンプル値の組に基づいて親個体集合Ｐの各個体について適応度の推定値の組を算出する（ステップＳ３）。適応度の推定値の算出方法については後述する。
次に、多目的進化型アルゴリズム部２は、適応度の推定値の組に基づく優劣比較およびパレートランキングにより親個体集合Ｐをランクごとの個体集合に分割する（ステップＳ４）。パレートランキングについては後述する。
次に、多目的進化型アルゴリズム部２は、親個体集合Ｐの各ランクの個体集合に混雑度ソートを行う（ステップＳ５）。それにより、各ランクの個体が混雑度（混雑距離）の大きい順に並べられる。混雑度ソートについては後述する。そして、より上位のランクでより大きな混雑度を有する所定数の個体が選択され、他の個体が削除される。
さらに、多目的進化型アルゴリズム部２および適応度推定モジュール３０ｂは、親個体集合Ｐの最上位ランク（ランク１）の個体集合から特定の１つの親個体を選択し、親個体集合Ｐからランダムに２つの個体を選択し、合計３つの親個体を選択するとともに、３つの親個体に交叉操作を施すことにより子個体集合Ｃを生成する（ステップＳ６ａ）。この場合、多数の子個体候補が生成され、後述する方法で子個体候補から所定数の子個体が選択される。子個体集合Ｃの生成処理の詳細については後述する。
多目的進化型アルゴリズム部２は、生成された子個体集合Ｃの各個体のパラメータの組を最適化対象６に順次与える（ステップＳ６ｂ）。それにより、最適化対象６から適応度のサンプル値の組が順次出力される。
ここで、交叉操作とは、個体の遺伝子を掛け合わせることにより新たな個体を生成することをいう。特定の親個体の選択方法については後述する。
適応度推定部３の探索履歴記憶装置３１は、最適化対象６から出力される子個体集合Ｃの各個体に対応する適応度のサンプル値の組を取得し、各個体に対応するパラメータの組およびサンプル値の組を探索履歴として記憶する（ステップＳ７）。
次いで、多目的進化型アルゴリズム部２は、子個体集合Ｃと親個体集合Ｐとから個体集合Ｆを生成する（ステップＳ８）。
適応度推定部３の適応度推定モジュール３０ａは、探索履歴記憶装置３１に記憶された複数組のサンプル値に基づいて個体集合Ｆの各個体に対応する適応度の推定値の組を算出する（ステップＳ９）。
多目的進化型アルゴリズム部２は、適応度の推定値の組に基づく優劣比較およびパレートランキングにより個体集合Ｆをランクごとの個体集合に分割し、子個体集合Ｃにおいて親個体集合Ｐの個体と重複する個体に最下位ランクを付与する（ステップＳ１０）。
次に、多目的進化型アルゴリズム部２は、個体集合Ｆの各ランクの個体集合に混雑度ソートを行って新たな親個体集合Ｐを生成する（ステップＳ１１）。それにより、各ランクの個体が混雑度（混雑距離）の大きい順に並べられる。混雑度ソートについては後述する。そして、より上位のランクでより大きな混雑度を有する所定数の個体が選択され、他の個体が削除される。
その後、多目的進化型アルゴリズム部２は、世代数が所定の終了条件に到達したか否かを判定する（ステップＳ１２）。
ここで、世代とは、個体集合から親個体を選択する選択ステップ、交叉操作により新たな子個体を生成する交叉ステップおよび親個体と子個体とを入れ替える世代交代ステップから構成される。世代数が所定の終了条件に到達していないと判定した場合には、ステップＳ６ａに移行する。世代数が所定の終了条件に到達したと判定した場合には、多目的進化型アルゴリズム部２は、ステッブＳ１１で生成された親個体集合Ｐをパレート最適個体集合として使用者１０に提示し、処理を終了する。
（ｅ）多目的最適化装置の各処理の具体例
図７〜図１２は多目的最適化装置１の各処理の具体例を示す模式図である。
図７〜図１２には、２つのパラメータｘ1，ｘ2および２つの適応度関数ｆ1，ｆ2の例が示される。図３の最適化対象６の場合には、パラメータｘ1，ｘ2は点火時期および燃料噴射時期であり、適応度関数ｆ1，ｆ2はＨＣ濃度およびＮＯx 濃度である。
（ｅ−１）親個体集合の初期化
図７は初期化により生成される親個体集合を示す模式図であり、（ａ）は適応度関数空間における親個体集合を示し、（ｂ）はパラメータ空間における親個体集合を示す。初期化においては、図７に示すように、適応度関数空間およびパラメータ空間に複数の個体がランダムに生成される。
（ｅ−２）個体の評価方法
多目的最適化問題においては、個体が複数の適応度関数に対応する適応度を有するため、単純な値の大小では個体の優劣を比較できない。本実施の形態では、以下に説明する優劣比較、パレートランキングおよび混雑度ソートを用いて個体を評価する。
（ｅ−２−１）優劣比較
図５のステップＳ４および図６のステップＳ１０における優劣比較について説明する。この優劣比較には、以下に示すα優越戦略（α-domination strategy）が用いられる。なお、α優越戦略の詳紬については、例えば、非特許文献５に掲載されている。
図８はα優越戦略を説明するための模式図である。ここで、一般に、α優越は、次のように定義される。
あるｐ個の目的度関数ｆk（ｋ＝１，・・・，ｐ）の２つの解ｘ１，ｘ２∈Ｆに対して次式（８）で表されるｇk（ｘ１，ｘ２）が次の関係を有する場合、解ｘ１は解ｘ２にα優越する。
ｇk（ｘ１，ｘ２）≦０（∀ｋ＝１，・・・，ｐ）∧ｇk（ｘ１，ｘ２）＜０（∃ｋ＝１，・・・，ｐ）
Figure 0005019744
上式において、ｆk（ｘ１）およびｆj（ｘ１）はそれぞれ解ｘ１に対応する目的関数ｆkおよびｆjの値であり、ｆk（ｘ２）およびｆj（ｘ２）はそれぞれ解ｘ２に対応する目的関数ｆkおよびｆjの値である。多目的遺伝的アルゴリズムでは、解ｘ１および解ｘ２が個体に相当し、目的関数ｆkおよびｆjが適応度関数に相当する。
ここで、図８において、個体Ｉ３に注目する。一般的な優越比較によれば、個体Ｉ３から適応度関数ｆ1に平行に延びる直線Ｌ１１および個体Ｉ３から適応度関数ｆ2に平行に延びる直線Ｌ１２で個体Ｉ３が他の個体に優越する領域が定められる。すなわち、個体Ｉ３は、直線Ｌ１１よりも上でかつ直線Ｌ１２より右の領域にある他の個体Ｉ６，Ｉ７，Ｉ８に優越する。個体Ｉ２，Ｉ４は、個体Ｉ３により優越されない。
個体Ｉ２は、適応度関数ｆ1に関しては個体Ｉ３によりもわずかに優れているが、適応度関数ｆ2に関しては個体Ｉ３よりもかなり劣っている。個体Ｉ４は、適応度関数ｆ2に関しては個体Ｉ３よりもわずかに優れているが、適応度関数ｆ1に関しては個体Ｉ３よりもかなり劣っている。このような個体Ｉ２，Ｉ４の適応度が不確実性（例えばノイズ）を有する場合には、個体Ｉ２，Ｉ４は個体Ｉ３により優越される可能性がある。
これに対して、α優越戦略によれば、個体Ｉ３から適応度関数ｆ1の軸に近づくように傾斜した直線Ｌ１および個体Ｉ３から適応度関数ｆ2の軸に近づくように傾斜した直線Ｌ２により個体Ｉ３が他の個体に優越する領域が定められる。すなわち、個体Ｉ３は、直線Ｌ１よりも上側でかつ直線Ｌ２より右側の領域にある個体Ｉ２，Ｉ４，Ｉ６，Ｉ７，Ｉ８に優越する。α優劣戦略によれば、個体Ｉ２，Ｉ４はパレート最適個体から排除される。
本実施の形態では、α優越戦略により複数の推定値の重み付線形和に基づいて個体の優劣比較が行われる。α優越戦略によれば、ある個体の１つの適応度が１悪くなれば、他の適応度はα悪くなる。すなわち、優劣比較において、１つの適応度の優劣が他の適応度の優劣に影響を与える。それにより、次のように、複数の適応度間の関係を考慮した合理的な解を求めることが可能となる。
弱パレート最適個体は、複数の適応度のうち少なくとも１つが他の個体に優越されない（すなわち少なくとも１つの適応度があるパレート最適解と等しい）解である。このような弱パレート最適個体は、ある適応度関数について最適解を有するが、残りの適応度関数についてはパレート最適個体に劣る。したがって、弱パレート最適個体は、合理的な解とは言えず、本来は求める必要の無い解である。そこで、α優越戦略を導入することにより、弱パレート最適個体を淘汰することができる。
また、適応度が不確実性を伴っている場合、弱パレート最適個体が不確実性によりパレート最適個体となる現象が生じる。弱パレート最適個体が不確実性によりパレート最適個体と判定されると、いつまでも淘汰されることなく個体集合中に存続することとなり、パレート最適個体の探索が停滞する原因となる。そこで、α優越戦略を導入することにより、このようなパレート最適個体と判定される弱パレート最適個体を淘汰することが可能となる。
図９はα優劣戦略による個体の優劣比較を説明するための模式図である。
図９において、個体Ｉ３は個体Ｉ２，Ｉ４，Ｉ５，Ｉ６，Ｉ７に優越しており、個体Ｉ６は個体Ｉ８に優越しており、個体Ｉ７は個体Ｉ８に優越している。個体Ｉ１，Ｉ３，Ｉ５に優越する個体はない。よって、個体Ｉ１，Ｉ３，Ｉ５パレート最適個体である。
（ｅ−２−２）パレートランキング
次に、図５のステップＳ４および図６のステップＳ１０のパレートランキングについて説明する。図１０はパレートランキングを説明するための図である。パレートランキングでは、各個体のランク付けに基づいてパレート最適個体集合を求める。
ｐi個の個体に優越されている個体ｘiのランクｒ（ｘi）は次式で与えられる。
ｒ（ｘi）＝１＋ｐi
ここでは、ランク１を最上位ランクとし、それ以上の数値のランクは数値が大きくなるほど下位のランクとなることにする。
図１０において、個体Ｉ１，Ｉ３，Ｉ５は他の個体に優越されていない。したがって、個体Ｉ１，Ｉ３，Ｉ５のランクは１である。個体Ｉ２，Ｉ４は１つの個体Ｉ３に優越されている。したがって、個体Ｉ２，Ｉ４のランクは２である。同様にして、個体Ｉ６のランクは６であり、個体Ｉ７のランクは５であり、個体Ｉ８のランクは８である。
（ｅ−２−３）混雑度ソート
次に、図５のステップＳ５および図６のステップＳ１１における混雑度ソートについて説明する。図１１は混雑度ソートを説明するための模式図である。
混雑度ソートでは、同じランクの各注目個体について、それに隣接する２つの個体を結ぶ線を対角線とする長方形を想定し、長方形の縦および横の辺の長さの合計で混雑度（混雑距離）を表す。混雑度の値が小さいほど注目個体は混雑した領域に存在する。同じランクの両端の個体には最大の混雑度を与える。
図１１において、個体Ｉ３の混雑度は、隣接する個体Ｉ１，Ｉ５が作る長方形ｓ１の縦および横の辺の合計で表される。個体Ｉ１の混雑度は、隣接する個体Ｉ９，Ｉ３が作る長方形ｓ２の縦および横の辺の合計で表される。個体Ｉ５の混雑度は、隣接する個体Ｉ３，Ｉ１０が作る長方形ｓ３の縦および横の辺の合計で表される。
図１２は多目的進化型アルゴリズム部２による混雑度ソートの処理を示すフローチャートである。
まず、多目的進化型アルゴリズム部２は、個体集合を適応度関数ごとにソートし、適応度関数ごとに同一ランク内で各注目個体に隣接する２つの個体を調べる（ステップＳ３１）。
次に、多目的進化型アルゴリズム部２は、各注目個体に隣接する２つの個体間の数学的距離を適応度関数ごとに算出し、各注目個体についての複数の適応度関数における数学的距離の合計を混雑度として算出する（ステップＳ３２）。ここで、数学的距離としてはユークリッド距離を用いる。
その後、多目的進化型アルゴリズム部２は、各ランクの個体集合の個体を混雑度の値の大きい順にソートする（ステップＳ３３）。
（ｅ−３）探索履歴による推定値の算出
次に、図５のステップＳ３および図６のステップＳ９における推定値の算出について説明する。
図１３は適応度推定部３の適応度推定モジュール３０ａによる推定値の算出を説明するための模式図である。
図１の探索履歴記憶装置３１には、多目的進化型アルゴリズム部２から与えられる個体のパラメータの組および最適化対象６から出力される適応度のサンプル値の組が探索履歴ＨＳとして順次記憶される。図１３においては、個体ごとにパラメータｘ1，ｘ2の組およびサンプル値Ｆ1，Ｆ2の組が探索履歴ＨＳとして記憶されている。
適応度推定モジュール３０ａは、探索履歴ＨＳに基づいて各個体に対応する真の適応度を推定値として算出する。各個体に対応するパラメータの組および推定値の組は、推定結果Ｅとして図１の探索履歴記憶装置３１に記憶される。
図１３に示すように、探索履歴記憶装置３１には、個体ごとにパラメータｘ1，ｘ2の組および推定値ｆ1’，ｆ2'の組が推定結果Ｅとして記憶されている。
また、適応度推定モジュール３０ａは、各個体のパラメータの組および推定値の組に基づいて適応度関数空間およびパラメータ空間上のパレート最適個体集合を図２の表示装置１０５の画面に表示することができる。
図１３においては、適応度関数ｆ1，ｆ2からなる適応度関数空間上およびパラメータｘ1，ｘ2からなるパラメータ空間上にパレート最適個体集合が表示されている。パレート最適個体集合が形成する境界をパレート境界と呼ぶ。
このように、探索履歴記憶装置３１に記憶された探索履歴ＨＳを用いて個体の適応度を推定する方法をメモリベース適応度推定法（Memory-based Fitness Estimated Method：ＭＦＥＭ)と呼ぶ（非特許文献５参照）。
注目個体の推定値を算出する場合、一般に注目個体と探索履歴ＨＳの個体とは異なる探索点である。また、不確実な環境を想定するため、最適化対象６に同じパラメータの組を与えても、異なるサンプル値の組が出力される。したがって、探索履歴ＨＳのサンプル値の組から注目個体の推定値の組を算出するためには適応度関数の性質に何らかの仮定を設ける必要がある。ＭＦＥＭでは、適応度関数が注目個体からの数学的距離に応じてランダムに変動すると考えて不確実な環境をモデル化している。
注目個体をｘとし、その注目個体ｘの真の適応度をｆ（ｘ）とする。パラメータ空間において注目個体ｘから距離ｄだけ離れた個体ｈの適応度ｆ（ｈ）を考える。適応度ｆ（ｈ）の期待値がｆ（ｘ）であり、適応度ｆ（ｈ）の分散が距離ｄに比例して増大する正規分布に従うモデルは次式（１）で表される。
ｆ（ｈ）〜Ｎ（ｆ（ｘ），ｋｄ） …（１）
上式（１）において、ｋは距離による重みを決定する正の定数であり、Ｎ（ｆ（ｘ），ｋｄ）は平均がｆ（ｘ）でかつ分散がｋｄである正規分布の確率密度関数を表す。
ここで、真の適応度ｆ（ｘ）には、平均０および分散σE 2 でかつ個体の位置に無関係な正規分布に従うノイズδが加わるものとする。この場合、個体ｘに対応するサンプル値Ｆ（ｘ）は次式のように定義される。
Ｆ（ｘ）＝ｆ（ｘ）＋δ …（２）
図１４は正規分布に従うノイズδを伴うサンプル値を示す模式図である。ここで、サンプル値Ｆ（ｘ）は適応度関数ｆ1についてのサンプル値Ｆ1（ｘ）および適応度関数ｆ2についてのサンプル値Ｆ2（ｘ）の組であり、真の適応度ｆ（ｘ）は適応度関数ｆ1についての真の適応度ｆ1（ｘ）および適応度関数ｆ2についての真の適応度ｆ2（ｘ）の組である。また、ノイズδは適応度関数ｆ1についてのノイズδ1および適応度関数ｆ2についてのノイズδ2の組である。ノイズδi（ｉ＝１，２）は次式で表される。
δi〜Ｎ（０，σEi 2 ） （ｉ＝１，２）
上式において、Ｎ（０，σEi 2 ）は平均０および分散σE 2である正規分布の確率密度関数を表す。
適応度推定部３は、サンプル値Ｆ（ｘ）の期待値を最小にするパレート最適個体集合を求める。このとき、個体ｈに対応するサンプル値Ｆ（ｈ）は、次式（３．１）および（３．２）としてモデル化される。
Ｆ（ｈ）〜Ｎ（ｆ（ｘ），ｋｄ＋σE 2 ） …（３．１）
ｄ＝｜ｈ−ｘ｜・・・（３．２）
上式（３．１）において、Ｎ（ｆ（ｘ），ｋｄ＋σE 2 ）は平均がｆ（ｘ）でありかつ分散がｋｄ＋σE 2 である正規分布の確率密度関数を表す。
図１５は不確実な適応度関数のモデルを示す模式図である。このモデルでは、注目個体ｘから遠く離れるほどサンプル値Ｆ（ｈ）が大きく不規則に変化するものと仮定している。
このモデルに基づいて探索履歴ＨＳを用いた最尤法により真の適応度ｆ（ｘ）の推定値を算出する。
探索履歴ＨＳに記憶された個体ｈl（ｌ＝１，…，Ｈ）、個体ｈlのサンプル値Ｆ（ｈl ）および注目個体ｘから個体ｈlまでの距離ｄl（ｌ＝１，…，Ｈ）を考えると、サンプル値Ｆ（ｈl），…，Ｆ（ｈH）が得られる確率は次式により表わすことができる。
Figure 0005019744
ここで、ｐ（Ｆ（ｈl），ｄl）は、サンプル値Ｆ（ｈl）が得られる確率を表す確率密度関数であり、次式で表すことができる。
Figure 0005019744
ここで、ｋ’＝ｋ／σE 2 である。本実施の形態においては、定数ｋ’は事前実験にて求めるものとする。非特許文献４には、探索中に定数ｋ’を推定する方法が提案されている。提案された方法により定数ｋ’を求めてもよい。
上記式（４）および（５）を真の適応度ｆ（ｘ）の尤度と考え、最尤法を用いる。それにより、真の適応度ｆ（ｘ）の推定値ｆ’（ｘ）は、次式（６）に示すように、距離ｄl を含む関数により重み付けされた加重平均の式で表すことができる。
Figure 0005019744
図１６は適応度推定部３の適応度推定モジュール３０ａによる推定値の算出処理を示すフローチャートである。
まず、適応度推定モジュール３０ａは、多目的進化型アルゴリズム部２において親個体集合Ｐが初期化されたことを確認する（ステップＳ４１）。次に、適応度推定モジュール３０ａは、探索履歴記憶装置３１の探索履歴ＨＳをすべてクリアする（ステップＳ４２）。
その後、適応度推定モジュール３０ａは、最適化対象６から出力されるサンプル値の組を探索履歴記憶装置３１に探索履歴ＨＳとして記憶する（ステップＳ４３）。次いで、適応度推定モジュール３０ａは、探索履歴記憶装置３１の探索履歴ＨＳに基づいて上記式（６）より各個体に対応する真の適応度の推定値の組を算出する（ステップＳ４４）。
適応度推定モジュール３０ａは、多目的進化型アルゴリズム部２の処理が終了したか否かを判定する（ステップＳ４５）。
多目的進化型アルゴリズム部２の処理が終了していない場合には、ステップＳ４３に戻ってステップＳ４３〜Ｓ４５の処理を繰り返す。多目的進化型アルゴリズム部２の処理が終了した場合には、推定値の算出処理を終了する。
（ｅ−４）親個体の選択から世代交代までの方法
次に、図６のステップＳ６ａの特定の親個体の選択からステップＳ１１の世代交代までの方法を説明する。図１７は特定の親個体の選択から世代交代までの方法を説明するための模式図である。
図１７（ａ）に示すように、パレートランキングにより親個体集合Ｐがランク付けされる。図１７（ｂ）に示すように、親個体集合Ｐのランク１の個体集合について、適応度関数空間における隣接する各２つの個体間のユークリッド距離を分布指標として算出する。本実施の形態では、適応度関数空間は、２つの適応度関数ｆ1およびｆ2からなる。
ユークリッド距離が最大となる２つの個体Ｉａ，Ｉｂのうちいずれか１つの個体を第１の親個体Ｉａとして確率１／２でランダムに選択する。さらに、第２の親個体Ｉｃおよび第３の親個体Ｉｄを親個体集合Ｐからランダム選択で選択する。
本実施の形態においては、分布指標のユークリッド距離Ｌは、図１７（ｆ）に示すように、隣接する２つの個体をｘおよびｙとすると、次式により求められる。
Ｌ＝［｛ｆ1（ｘ）−ｆ1（ｙ）｝2 ＋｛ｆ2（ｘ）−ｆ2（ｙ）｝2 ］1/2 …（７）
次いで、図１７（ｃ）に示すように、第１、第２および第３の親個体Ｉａ，Ｉｃ，Ｉｄから複数の子個体候補を生成し、後述する方法で複数の子個体候補から所定数の個体を含む子個体集合Ｃを生成する。さらに、図１７（ｄ）に示すように、子個体集合Ｃおよび親個体集合Ｐから個体集合Ｆを生成し、個体集合Ｆに上記のα優越戦略を用いたパレートランキングを行う。このとき、子個体集合Ｃにおいて親個体集合Ｐに含まれる個体に重複する個体がある場合には、その個体に最下位ランクを付与する。
その後、図１７（ｅ）に示すように、個体集合Ｆに混雑度ソートを行い、個体のランクおよび各ランク内の混雑度に基づいて所定数の個体を選択し、残りの個体を削除する。それにより、新たな親個体集合Ｐを生成する。このようにして、世代交代が行われる。
図１８は子個体候補の生成処理を説明するための模式図である。図１８（ａ）はパラメータ空間上の親個体集合Ｐを示し、図１８（ｂ）はパラメータ空間上の子個体候補集合を示す。子個体候補集合の個体Ｉ２１のパラメータｘ1，ｘ2が親個体集合Ｐの個体Ｉ１１のパラメータｘ1，ｘ2と一致する場合には、個体Ｉ２１には最下位ランクが付与される。同様に、子個体候補集合の個体Ｉ２２のパラメータｘ1，ｘ2が親個体集合Ｐの個体Ｉ１２のパラメータｘ1，ｘ2と一致する場合には、個体Ｉ２２には最下位ランクが付与される。
なお、特定の親個体の選択方法は、本実施の形態に限定されず、ランク１の個体集合からユークリッド距離Ｌが最大となる２つの個体を第１および第２の親個体として選択し、第３の親個体を親個体集合Ｐからランダム選択、ルーレット選択またはトーナメント選択等の方法により選択してもよい。
図１９は多目的進化型アルゴリズム部２による特定の親個体の選択処理を示すフローチャートである。
まず、多目的進化型アルゴリズム部２は、親個体集合Ｐから選択されたランク１の個体集合を適応度関数ごとにソートする（ステップＳ５１）。
次に、多目的進化型アルゴリズム部２は、ランク１の個体集合において隣接する２つの個体間のユークリッド距離を算出する（ステップＳ５２）。
さらに、多目的進化型アルゴリズム部２は、最大のユークリッド距離を与える２つの個体のうち１つの個体を第１の親個体として確率１／２でランダムに選択し、第２の親個体および第３の親個体を親個体集合Ｐからランダム選択で選択する（ステップＳ５３）。
本実施の形態では、第１、第２および第３の親個体から交叉操作により子個体候補が生成される。交叉操作としては、例えば、ＵＮＤＸ（単峰性正規分布交叉：Unimodal Normal Distribution Crossover）が用いられる。
図２０はＵＮＤＸによる子個体候補の生成処理を示す模式図である。ＵＮＤＸでは、第１の親個体Ｐ１、第２の親個体Ｐ２および第３の親個体Ｐ３の位置関係に基づいて定められる正規分布乱数に従って子個体候補Ｃ１を生成する。この場合、第１の親個体Ｐ１と第２の親個体Ｐ２とを結ぶ軸ＡＸの周辺の正規分布に従って子個体候補Ｃ１が生成されるので、子個体候補Ｃ１が第１〜第３の親個体Ｐ１〜Ｐ３から遠くに離れた位置に生成されることがない。
（ｅ−５）子個体集合Ｃの生成方法
図２１は多目的進化型アルゴリズム部２および適応度推定モジュール３０ｂによる子個体集合Ｃの生成処理を示すフローチャートである。
まず、多目的進化型アルゴリズム部２は、上記の方法で第１、第２および第３の親個体から複数の子個体候補を生成する（ステップＳ６１）。
次に、適応度推定モジュール３０ｂは、全ての親個体および子個体候補についての真の適応度の推定値を算出する（ステップＳ６２）。
子個体候補についての真の適応度ｆ（ｘ）の推定値ｆ’（ｘ）は、次式（６ａ）に示すように、距離ｄlを含む関数により重み付けされた加重平均の式で表すことができる。
Figure 0005019744
ただし、この段階では、子個体候補に対応するサンプル値は得られていないため、子個体候補についての推定値ｆ’（ｘ）の算出式は、上記式（６）とは異なり、その子個体候補のサンプル値Ｆ（ｘ）を含まない。
なお、推定値の算出方法としては、上式（６ａ）以外に、多項式回帰、ニューラルネットワーク、放射基底関数（Radial Basis Function）、クリギング、ガウシアンプロセス、局所重み付き回帰（Locally Weighted Regression）等の他の手法を用いることもできる。
さらに、多目的進化型アルゴリズム部２は、全ての親個体および子個体候補のランキングを行う（ステップＳ６３）。このランキングには、図８のα優越戦略、図９の優劣比較および図１０のパレートランキングを用いる。
次に、多目的進化型アルゴリズム部２は、最上位ランクの１つの子個体候補を親個体集合Ｐに追加して混雑度を算出する（ステップＳ６４）。混雑度の算出方法については、後述する。
次に、多目的進化型アルゴリズム部２は、最上位ランクの全ての子個体の混雑度の算出が終了したか否かを判別する（ステップＳ６５）。最上位ランクの全ての子個体候補の混雑度の算出が終了していない場合には、ステップＳ６４に戻る。
最上位ランクの全ての子個体候補の混雑度の算出が終了した場合には、それらの子個体候補を混雑度の良い順にソートし（ステップＳ６６）、より良い混雑度を有する所定数の子個体候補を子個体集合Ｃとして選択する（ステップＳ６７）。
図２２は子個体集合Ｃの生成処理を示す模式図であり、（ａ）は適応度関数空間における親個体および子個体候補を示し、（ｂ）は適応度空間における混雑度の算出方法を示す。
子個体候補の混雑度の算出方法は、図１１に示した注目個体の混雑度の算出方法と同様である。
図２２（ａ）の適応度関数空間には、親個体Ｐａ〜Ｐｄおよび生成された子個体候補Ｃ１１〜Ｃ２０が示されている。図２２（ｂ）には、親個体Ｐａ〜Ｐｄおよび最上位ランクの１つの子個体候補Ｃ１２が示されている。
図２２（ａ）に示すように、適応度関数空間において親個体から多数の子個体候補が生成される。
混雑度の算出では、１つの子個体候補について、それに隣接する２つの親個体を結ぶ線を対角線とする長方形を想定し、長方形の縦および横の辺の長さの合計で混雑度（混雑距離）を表す。混雑度の値が小さいほど子個体候補は混雑した領域に存在する。同じランクの両端の子個体候補には最大の混雑度を与える。
図２２（ｂ）において、子個体候補Ｃ１２の混雑度は、隣接する親個体Ｐｂ，Ｐｃが作る長方形の縦および横の辺の合計で表される。
（ｆ）第１の実施の形態の効果
本実施の形態に係る多目的最適化装置１においては、多数の子個体候補が生成され、最上位ランクの複数の子個体候補から混雑度が高い順に所定数の子個体候補が子個体として選択され、選択された子個体が実際に評価される。それにより、不要な子個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
また、上式（６）により探索履歴記憶装置３１に記憶された各個体に対応するサンプル値の組に重み付けが行われ、重み付けられた複数組のサンプル値の線形和を用いて注目個体に対応する適応度の推定値の組が求められる。
各個体の重みは、パラメータ空間上で注目個体とその個体との距離を含む関数であるので、真の適応度からの誤差が十分に小さい推定値を得ることができる。したがって、最適化対象から出力されるサンプル値が不確実性を有する場合でも、適切なパレート最適個体集合を得ることができる。
また、α優越戦略により複数の適応度関数の各々について個体集合の複数の個体に対応する推定値の優劣が比較され、複数の適応度関数の各々についての比較結果に重み付けが行われる。そして、複数の適応度関数について重み付けられた複数の比較結果の線形和に基づいて個体集合の複数の個体のランク付けが行われる。それにより、最適化対象が不確実性を有する場合でも、複数の適応度間の関係を考慮した合理的なパレート最適個体を求めることが可能となる。
さらに、適応度関数空間上での最上位ランクの個体の分布において、隣接する個体間の距離を分布指標として用いることにより、疎らな領域に新たな子個体を容易に生成することができる。それにより、最上位ランクの個体を適応度関数空間上の広い領域に偏りなく分布するように生成することが可能となる。したがって、多様性を有するパレート最適個体を容易に得ることができる。
さらに、生成された新たな子個体が親個体集合の個体と重複する場合に、新たな子個体に最下位ランクが付与される。それにより、パレート最適個体の探索初期には、緩やかに悪い個体を減少させることができ、パレート最適個体の探索後期には、パレート最適個体の多様性を維持することができる。
また、適応度推定部３により算出された推定値の組に基づいてパレート最適個体が表示装置１０５の画面上に表示される。それにより、使用者は、パレート最適個体を視覚的に認識することができるので、種々の意思決定を容易に行うことができる。
（２）第２の実施の形態
次に、本発明の第２の実施の形態に係るに多目的最適化装置について説明する。本実施の形態に係る多目的最適化装置は、図１および図２に示した構成を有する。また、本実施の形態に係る多目的最適化装置の全体処理は、図５および図６に示した処理と同様である。
本実施の形態が第１の実施の形態と異なるのは、図５のステップＳ３および図６のステップＳ９における推定値の算出方法および図６のステップＳ６ａにおける特定の親個体の選択方法である。
（ａ）推定値の算出
本実施の形態では、真の適応度ｆ（ｘ）の推定値ｆ’（ｘ）は、改良された推定式により算出される。
Figure 0005019744
上式（９）に示すように、推定値ｆ’（ｘ）はパラメータ空間上の距離ｄlの３乗を含む関数により重み付けされた加重平均の式で表される。
上式（９）によれば、注目個体から探索履歴ＨＳ内の個体までの距離が短いほど重みが大きくなる。一方、注目個体から探索履歴ＨＳ内の個体までの距離が長くなると、極端に重みが小さくなる。したがって、注目個体から離れた個体は推定値ｆ’（ｘ）の算出にほとんど寄与しない。
図２３は探索履歴記憶装置３１に記憶された探索履歴に基づく個体の探索を示す模式図である。図２３（ａ）は単目的最適化における探索初期の個体集合を示し、図２３（ｂ）は多目的最適化における探索初期の個体集合を示し、図２３（ｃ）は単目的最適化における探索後期の個体集合を示し、図２３（ｄ）は多目的最適化における探索後期の個体集合を示す。
図２３（ａ），（ｃ）の縦軸は適応度関数ｆを示し、横軸はパラメータｘを示す。図２３（ｂ），（ｄ）の縦軸は適応度関数ｆ2を示し、横軸は適応度関数ｆ1を示す。
探索初期には、図２３（ａ），（ｂ）に示すように複数の個体が分散している。探索後期には、単目的最適化では、図２３（ｃ）に示すように個体があるパラメータの値の近傍に集中し、多目的最適化では、図２３（ｄ）に示すように複数の個体がパレート最適個体集合を形成する。
このように、多目的最適化では、個体が適応度関数空間の広い範囲に分散するので、式（９）を用いることにより、注目個体から遠く離れた個体の寄与が小さくなるので、推定値の精度が高くなる。
（ｂ）特定の親個体の選択から世代交代までの方法
図２４は親個体の選択から世代交代までの方法を説明するための模式図である。本実施の形態では、３つの適応度関数ｆ1，ｆ2，ｆ3の例が示される。
図２４（ａ）に示すように、パレートランキングにより親個体集合Ｐがランク付けされている。図２４（ｂ）に示すように、親個体集合Ｐのランク１の個体集合について、適応度関数空間における隣接する各３つの個体が形成する三角形の面積を分布指標として算出する。
三角形の形成には、ドローネ三角形分割（Delaunay Triangulation）の方法を用いる（非特許文献６参照）。
ここで、ドローネ三角形分割について簡単に説明する。ドローネ三角形分割は、計算幾何学（Computational Geometry）の中で重要な概念であるボロノイ図（Voronoi Diagram）の双対図形である。平面（空間）上の点集合の三角形分割の中で、種々の意味で最適な三角形分割として知られており、コンピュータグラフィクスにおけるメッシュ生成または有限要素法等にも用いられる。ドローネ三角形分割は、分割された三角形の最小角を最大にする分割方法であり、アルゴリズムとして逐次添加法、分割統治法または幾何変換を用いる方法等がある。なお、ドローネ三角形分割の詳細については、例えば非特許文献７に記載されている。
最大の三角形の面積を与える３つの個体ＩＡ，ＩＢ，ＩＣのうちいずれか１つの個体を第１の親個体ＩＡとして確率１／３でランダムに選択する。さらに、第２の親個体ＩＤおよび第３の親個体ＩＥを親個体集合Ｐからランダム選択で選択する。
次いで、図２４（ｃ）に示すように、第１、第２および第３の親個体ＩＡ，ＩＤ，ＩＥから複数の子個体候補を生成し、上記の方法で複数の子個体候補から所定数の個体を含む子個体集合Ｃを生成する。さらに、図２４（ｄ）に示すように、子個体集合Ｃおよび親個体集合Ｐから個体集合Ｆを生成し、個体集合Ｆに上記のα優越戦略を用いたパレートランキングを行う。このとき、子個体集合Ｃにおいて親個体集合Ｐに含まれる個体に重複する個体がある場合には、その個体に最下位ランクを付与する。
その後、図２４（ｅ）に示すように、個体集合Ｆに混雑度ソートを行い、個体のランクおよび各ランク内の順位に基づいて所定数の個体を選択し、残りの個体を削除する。それにより、新たな親個体集合Ｐを生成する。このようにして、世代交代が行われる。
なお、特定の親個体の選択方法は、本実施の形態に限定されず、最大の三角形の面積を与える３つの個体を第１、第２および第３の親個体として選択してもよく、あるいは最大の三角形の面積を与える３つの個体のうち２つを第１および第２の親個体として選択し、第３の親個体を親個体集合Ｐからランダム選択、ルーレット選択またはトーナメント選択等の方法により選択してもよい。
また、パラメータが３つ以上の場合には、多親拡張された交叉方法、例えばＵＮＤＸ−ｍ等を用いてもよい。
図２５は多目的進化型アルゴリズム部２による特定の親個体の選択処理を示すフローチャートである。
まず、多目的進化型アルゴリズム部２は、親個体集合Ｐから選択されたランク１の個体集合をｆi−ｆj平面に正射影する（ステップＳ６１）。ここで、ｆiおよびｆjは適応度関数である。ｉ，ｊ＝１，２，３であり、かつｉ≠ｊであり、それらの組み合わせは世代ごとに変化させる。
次に、多目的進化型アルゴリズム部２は、正射影した個体集合のドローネ三角形分割を行う（ステップＳ６２）。
さらに、多目的進化型アルゴリズム部２は、ドローネ三角形分割されたランク１の個体集合に適応度関数ｆkの成分を高さ成分として与え、複数の三角形を３次元空間上に展開する（ステップＳ６３）。ここで、ｋ≠ｉ，ｊである。
次いで、多目的進化型アルゴリズム部２は、３次元空間上に展開された複数の三角形の面積をそれぞれ算出する（ステップＳ６４）。
面積最大の三角形を形成する３つの個体のうち１つの個体を第１の親個体として確率１／３でランダムに選択し、親個体集合Ｐから２つの個体を第２の親個体および第３の親個体としてランダム選択で選択する（ステップＳ６５）。
（ｃ）子個体集合Ｃの生成方法
本実施の形態における子個体集合Ｃの生成処理は、図２０に示した処理と同様である。
ただし、子個体候補についての真の適応度ｆ（ｘ）の推定値ｆ’（ｘ）は、改良された推定式により算出される。
Figure 0005019744
上式（９ａ）に示すように、推定値ｆ’（ｘ）はパラメータ空間上の距離ｄlの３乗を含む関数により重み付けされた加重平均の式で表される。
ただし、この段階では、子個体候補に対応するサンプル値は得られていないため、子個体候補についての推定値ｆ’（ｘ）の算出式は、上記式（９）とは異なり、その子個体候補のサンプル値Ｆ（ｘ）を含まない。
（ｄ）第２の実施の形態の効果
本実施の形態に係る多目的最適化装置１においては、多数の子個体候補が生成され、最上位ランクの複数の子個体候補から混雑度が高い順に所定数の子個体候補が子個体として選択され、選択された子個体が実際に評価される。それにより、不要な子個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。したがって、適切なパレート最適個体を短時間で効率良く得ることが可能となる。
また、上式（９）により探索履歴記憶装置３１に記憶された各個体に対応するサンプル値の組に重み付けが行われ、重み付けられた複数組のサンプル値の線形和を用いて注目個体に対応する適応度の推定値の組が求められる。
各個体の重みは、パラメータ空間上で注目個体とその個体との距離の３乗を含む関数であるので、パラメータ空間上で注目個体から遠く離れた他の個体の影響が十分に小さくなる。それにより、真の適応度からの誤差が十分に小さい推定値を得ることができる。したがって、最適化対象から出力されるサンプル値が不確実性を有する場合でも、適切なパレート最適個体集合を得ることができる。
また、適応度関数空間上での最上位ランクの個体の分布において、隣接する３つの個体を頂点とする三角形の面積を分布指標として用いることにより、疎らな領域に新たな子個体を容易に生成することができる。それにより、最上位ランクの個体を適応度関数空間上の広い領域に偏りなく分布するように生成することが可能となる。したがって、多様性を有するパレート最適個体を容易に得ることができる。
（３）他の実施の形態
（ａ）拡張された推定式
上記の推定式（４）および（９）を一般化すると、次式のようになる。
Figure 0005019744
上式（１０）において、ｎは任意の自然数である。第１の実施の形態の推定式（４）はｎ＝１の場合を示し、第２の実施の形態の推定式（９）はｎ＝３の場合を示す。ｎ＝３が好ましいが、ｎが他の自然数であってもよい。
このように、注目個体の真の適応度の推定値は、ノイズ成分を有する適応度のサンプル値を用いずに、探索履歴ＨＳにおける他の個体の推定値の重み付け線形和により算出される。それにより、サンプル値が不確実性を有する場合でも、パレート最適個体を安定に探索することができる。
また、パラメータ空間上で注目個体と他の個体との距離のｎ乗の関数を含む重みが用いられるので、推定値の算出が広範囲に広がる他の個体の影響を大きく受けることが防止される。したがって、推定値を高精度に算出することができる。
（ｂ）親個体の選択
第１の実施の形態で示されるように、２目的最適化問題では、特定の親個体の選択のための分布指標として個体間の距離が用いられる。また、第２の実施の形態で示されるように、３目的最適化問題では、特定の親個体の選択のための分布指標として３つの個体が形成する三角形の面積が用いられる。
特定の親個体の選択のための分布指標をｍ目的最適化問題に拡張した場合、分布指標は、適応度関数空間上で隣接するｍ個の個体が形成する単体(simplex)の大きさである。ｍは２以上の自然数である。上記の単体は、ドローネ三角形分割により形成することができる。
図２６はｍ目的最適化問題に拡張された分布指標を示す図である。図２６に示すように、２目的の場合には、分布指標は隣接する２つの個体間を結ぶ直線の長さであり、３目的の場合には、分布指標は隣接する３つの個体を頂点とする三角形の面積であり、４目的の場合には、分布指標は隣接する４つの個体を頂点とする三角錐の体積である。分布指標は４次元単体の大きさであり、底体積×高さ÷４により算出される。５目的の場合には、分布指標は５次元単体の大きさであり、底４次元面積×高さ÷５により算出される。（ｍ＋１）目的の場合には、分布指標はｍ次元単体の大きさであり、底（ｍ−１）次元面積×高さ÷ｍにより算出される。
このように、分布指標に基づいて特定の親個体を選択することにより、適応度関数空間上で分布が疎らな領域で積極的な個体の探索が行われる。それにより、広い領域で個体の探索が行われるので、推定値を高精度に算出することができるとともに、適応度関数空間上の広い領域で均等にパレート最適個体を見出すことができる。
（ｃ）子個体集合の生成時の拡張された推定式
上記の推定式（４ａ）および（９ａ）を一般化すると、次式のようになる。
Figure 0005019744
上式（１０ａ）において、ｎは任意の自然数である。第１の実施の形態の推定式（４ａ）はｎ＝１の場合を示し、第２の実施の形態の推定式（９ａ）はｎ＝３の場合を示す。ｎ＝３が好ましいが、ｎが他の自然数であってもよい。
なお、上式（１０ａ）の近似モデルとして、多項式回帰、ニューラルネットワーク、放射基底関数（Radial Basis Function）、クリギング、ガウシアンプロセス、局所重み付き回帰（Locally Weighted Regression）等の他の手法を用いることもできる。
（ｄ）エンジンシミュレータへの適応例
図２７は多目的最適化装置をエンジンシミュレータに適用した例を示すブロック図である。
図２７の最適化対象６ａはエンジンシミュレータである。エンジンシミュレータは、例えばパーソナルコンピュータからなる。この最適化対象６ａは、多目的最適化装置１から与えられるパラメータの組に基づいてエンジンの動作のシミュレーションを行い、シミュレーション結果を適応度のサンプル値の組として多目的最適化装置１に出力する。
本実施の形態では、複数の適応度関数としては、燃費、トルク、エンジンの排気ガスに含まれるＣＯ（一酸化炭素）、ＨＣ（炭化水素）、ＮＯx （窒素酸化物）等の成分の濃度等のうち複数が設定される。
また、パラメータとしては、燃料噴射量、燃料噴射時期、点火時期、スロットル開度等が挙げられる。
図２７の多目的最適化装置１によれば、適応度関数の組およびパラメータの組を設定することにより、パレート最適個体集合を短時間で効率良く求めることができる。
（ｅ）モータ評価装置への適応例
図２８は多目的最適化装置をモータ評価装置に適用した例を示すブロック図である。
図２８の最適化対象６ｂはモータ評価装置である。モータ評価装置は、モータ、制御回路および各種検出回路により構成される。この最適化対象６ｂは、多目的最適化装置１から与えられる個体のパラメータの組に基づいてモータを制御するとともにモータの複数の性能項目を測定し、測定結果を適応度のサンプル値の組として多目的最適化装置１に出力する。
複数の適応度関数としては、立ち上がり時間、整定時間、オーバーシュート量、消費電流等のうち複数が設定される。
また、パラメータとしては、ＰＩＤ（比例積分微分：Proportional Integral Derivative）ゲイン、駆動電流等が挙げられる。
ここで、トレードオフの関係としては、立ち上がり時間とオーバーシュート量、立ち上がり時間と消費電流、整定時間とオーバーシュート量等が挙げられる。
図２８の多目的最適化装置１によれば、適応度関数の組およびパラメータの組を設定することにより、パレート最適個体集合を短時間で効率良く求めることができる。
また、パレート最適個体集合をリアルタイムに算出することにより実環境に即したモータのリアルタイム制御を行うことも可能である。
（ｆ）モータシミュレータへの適応例
図２９は多目的最適化装置をモータシミュレータに適用した例を示すブロック図である。
図２９の最適化対象６ｃはモータシミュレータである。モータシミュレータは、例えばパーソナルコンピュータからなる。この最適化対象６ｃは、多目的最適化装置１から与えられるパラメータの組に基づいてモータの動作のシミュレーションを行い、シミュレーション結果を適応度のサンプル値の組として多目的最適化装置１に出力する。
複数の適応度関数としては、立ち上がり時間、整定時間、オーバーシュート量、消費電流等のうち複数が設定される。また、パラメータとしては、ＰＩＤゲイン、駆動電流等が挙げられる。
図２９の多目的最適化装置１によれば、適応度関数の組およびパラメータの組を設定することにより、パレート最適個体集合を短時間で効率良く求めることができる。
（ｇ）多目的進化型アルゴリズムの他の例
上記実施の形態では、多目的進化型アルゴリズムとして遺伝的アルゴリズム（ＧＡ）を用いているが、これに限定されず、遺伝的アルゴリズムの代わりに、進化戦略（ＥＳ：Evolution Strategy）等の同様のアイデアに基づく計算法を用いてもよい。
なお、ＧＡ、ＥＳ等の計算法は、進化アルゴリズム（ＥＡｓ：Evolutionary Algorithms）または進化計算（Evolutionary Computation）と総称される。
（ｈ）４以上の目的への適用
上記実施の形態では、２目的および３目的の最適化を例に挙げて説明したが、本発明は、４以上の目的の最適化にも同様に適用することができる。この場合、トレードオフの関係を有する４以上の適応度関数が設定される。
（ｉ）世代交代方法
上記の図６のステップＳ１０において、子個体集合Ｃのうち親個体集合Ｐの個体と重複しない個体を親個体集合Ｐの下位の個体と入れ換えることにより新たな親個体集合Ｐを生成してもよい。
この場合、パレート最適個体の探索初期には、緩やかに悪い個体を減少させることができるとともに、パレート最適個体の探索後期には、パレート最適個体の多様性を維持することができる。
（ｊ）親個体の再評価
上記実施の形態においては、親個体の再評価を行っているが、実システムまたは大規模シミュレーションにおいて個体の評価回数に制限がある場合には、子個体のみ再評価してもよい。それにより、評価回数を低減することが可能になる。
（ｋ）サンプル値の取得の制限
探索履歴記憶装置３１に記憶されるサンプル値の量が所定の記憶容量に達した時点で探索履歴記憶装置３１へのサンプル値の取得を終了してもよい。それにより、以後は探索履歴記憶装置３１に記憶された探索履歴ＨＳに基づいて推定値を算出し、算出された推定値に基づいてパレート最適個体の探索を進めることが可能になる。
（ｌ）ランク付け
上記実施の形態では、パレートランキングにより複数の個体のランク付けが行われているが、これに限定されず、非優越ソート等の他の方法を用いて複数の個体がランク付けされてもよい。
（ｍ）最適化対象の不確実性を考慮しない多目的最適化装置への適用
上記実施の形態では、最適化対象６が観測ノイズ等の不確実性を有する場合でも多様性を有するパレート最適個体を短時間で得ることができる多目的最適化装置１について説明したが、本発明は、最適化対象の観測ノイズ等の不確実性を考慮しない多目的最適化装置にも適用することができる。この場合にも、不要な子個体候補を評価する時間が省略される。その結果、最適化の速度が向上する。
（ｎ）各部の実現方法
上記実施の形態では、多目的進化型アルゴリズム部２、適応度推定モジュール３０ａ，３０ｂおよび探索履歴記憶部３１がＣＰＵ１０１およびプログラムにより実現されるが、多目的進化型アルゴリズム部２、適応度推定モジュール３０および探索履歴記憶部３１の一部または全てが電子回路等のハードウエアにより実現されてもよい。
（４）実施例および比較例
以下の実施例では、第１および第２の実施の形態に係る多目的最適化装置によりベンチマーク問題を実行した。また、比較例では、子個体集合Ｃの生成方法を除いて第１および第２の実施の形態に係る多目的最適化装置と同様の多目的最適化装置によりベンチマーク問題を実行した。
図３０は実施例および比較例の多目的最適化の条件を示す図である。図３１は実施例および比較例において実行するベンチマーク問題を示す。
図３１に示すように、評価回数を１５００として実施例と比較例とを比較した。なお、参考のため、比較例では、３０００回まで評価を行った。実施例および比較例において、個体集合サイズは１００である。実施例では、第２の実施の形態の子個体集合の生成方法により１００個の子個体候補から３個の子個体を選択し、３個の子個体および２個の親個体を評価した。一方、比較例では、特定の親個体から８個の子個体を生成し、８個の子個体および２個の親個体を評価した。これは、実施例の評価回数が比較例の評価回数の半分であっても良好な結果が得られることを確認するためである。
また、実施例および比較例では、交叉オペレータとしてＵＮＤＸを用い、選択オペレータとして第１および第２の実施の形態の選択方法（図１７および図２４に示した特定の親個体の選択方法）を用いた。さらに、図１７および図２４に示した世代交代方法を用いた。さらに、最上位ランクの決定方法としては、パレートランキングを用いた。
ここで、評価回数とは、（世代数×子個体の数＋初期個体集合の個体の数）で求められる。
実施例および比較例では、ベンチマーク問題として図３１に示すテスト関数ＺＤＴ１、ＺＤＴ２およびＤＴＬＺ２を用いた（非特許文献８参照）。
テスト関数ＺＤＴ１は、パレート最適解集合が凸型のパレート境界を有する２目的最適化問題である。テスト関数ＺＤＴ２は、パレート最適解集合が凹型のパレート境界を有する２目的最適化問題である。ここで、テスト関数ＺＤＴ１およびテスト関数ＺＤＴ２を２変数２目的関数の最小化問題とした。なお、このテスト関数ＺＤＴ１およびテスト関数ＺＤＴ２は弱パレート最適解を持つ。また、それぞれの目的関数には適当なノイズを加える。
テスト関数ＤＴＬＺ２は、パレート最適解集合が凹型のパレート境界を有する３目的最適化問題である。ここで、テスト関数ＤＴＬＺ２を３変数３目的関数の最小化問題とした。また、それぞれの目的関数には適当なノイズを加える。
テスト関数ＺＤＴ１およびテスト関数ＺＤＴ２の定数ｋｉ'（ｉ＝１，２ ）は予備実験にて決定し、ｋｉ'＝３１００００とした。テスト関数ＤＴＬＺ２の定数ｋｉ'（ｉ＝１，２，３）は予備実験にて決定し、ｋｉ'＝３１０００とした。また、α優越戦略におけるαｉｊ（ｉ＝１，…ｎ；ｊ＝１，…，ｎ；ｉ≠ｊ）は０．１とした。ここで、ｎは目的の数である。
ここで、実施例および比較例において、平均絶対誤差Ａｅｒｒｏｒおよび平均推定誤差Ｅｅｒｒｏｒを算出した。
次式はテスト関数ＺＤＴ１およびテスト関数ＺＤＴ２における平均絶対誤差Ａｅｒｒｏｒの算出式である。
Figure 0005019744
次式はテスト関数ＤＴＬＺ２における平均絶対誤差Ａｅｒｒｏｒの算出式である。
Figure 0005019744
上式（１１），（１２）において、｜Ｐ｜は親個体集合Ｐの数を示し、ｍはパラメータ（決定変数）の数を示し、ｘｉｊはパラメータを示す。また、実施例および比較例において、平均推定誤差Ｅｅｒｒｏｒを次式により算出した。
Figure 0005019744
上式（１３）において、｜Ｐ｜は親個体集合Ｐの数を示し、ｎは適応度関数（目的関数）の数を示し、ｆｊ（ｘｉ）は真の適応度の値を示し、ｆ’ｊ（ｘｉ）は真の適応度の推定値を示す。
実施例および比較例では、平均絶対誤差Ａｅｒｒｏｒおよび平均推定誤差Ｅｅｒｒｏｒを３０回の試行で評価した。
図３２（ａ）はテスト関数ＺＤＴ１およびテスト関数ＺＤＴ２における平均絶対誤差Ａｅｒｒｏｒの概念を示す模式図、図３２（ｂ）はテスト関数ＤＴＬＺ２における平均絶対誤差Ａｅｒｒｏｒの概念を示す模式図である。
図３２（ａ）において、ノイズ成分がない場合には、個体ｉ１０１，ｉ１０２，ｉ１０３はパラメータ空間のｘ２＝０の直線上に存在する。最適化により得られた個体Ｉ１０１，Ｉ１０２，Ｉ１０３とノイズ成分がない場合の個体ｉ１０１，ｉ１０２，ｉ１０３との距離の平均値を算出することによりテスト関数ＺＤＴ１およびテスト関数ＺＤＴ２における平均絶対誤差Ａｅｒｒｏｒが得られる。
図３２（ｂ）において、ノイズ成分がない場合には、個体ｉ１１１，ｉ１１２，ｉ１１３はパラメータ空間のｘ３＝０．５の平面上に存在する。最適化により得られた個体Ｉ１１１，Ｉ１１２，Ｉ１１３とノイズ成分がない場合の個体ｉ１１１，ｉ１１２，ｉ１１３との距離の平均値を算出することによりテスト関数ＤＴＬＺ２における平均絶対誤差Ａｅｒｒｏｒが得られる。
図３３は２目的の場合における平均絶対誤差Ａｅｒｒｏｒの概念を示す模式図である。
図３３において、適応度関数空間に真の適応度ｆ（ｘａ），ｆ（ｘｂ），ｆ（ｘｃ）および真の適応度の推定値ｆ’（ｘａ），ｆ’（ｘｂ），ｆ’（ｘｃ）が示されている。
真の適応度ｆ（ｘａ）と推定値ｆ’（ｘａ）との間における適応度関数ｆ１についての距離と適応度関数ｆ２についての距離との合計、真の適応度ｆ（ｘｂ）と推定値ｆ’（ｘｂ）との間における適応度関数ｆ１についての距離と適応度関数ｆ２についての距離との合計、および真の適応度ｆ（ｘｃ）と推定値ｆ’（ｘｃ）との間における適応度関数ｆ１についての距離と適応度関数ｆ２についての距離との合計を算出し、それらの合計の平均を算出することにより平均絶対誤差Ａｅｒｒｏｒが得られる。
図３４は実施例および比較例における平均絶対誤差および平均推定誤差の算出結果を示す図である。
図３４には、平均絶対誤差および平均推定誤差に加えて、参考のために絶対誤差の標準偏差および推定誤差の標準偏差が示されている。また、参考のために、比較例については評価回数が３０００回の場合の平均絶対誤差および平均推定誤差の算出結果も示されている。
図３４に示すように、テスト関数ＺＤＴ１、テスト関数ＺＤＴ２およびテスト関数ＤＴＬＺ２における実施例の平均絶対誤差は、比較例に比べて小さくなり、比較例の評価回数を３０００回とした場合に比べても小さくなった。
また、テスト関数ＺＤＴ１およびテスト関数ＺＤＴ２における実施例の平均推定誤差は、比較例に比べて小さくなり、比較例の評価回数を３０００回とした場合に比べても小さくなった。
これは、実施例では、パレート最適個体の分布が広がりにくいことにより、パレート最適個体のある領域で探索が進みやすいからであると考えられる。なお、テスト関数ＤＴＬＺ２における実施例の平均推定誤差のみが比較例に比べて大きくなった。テスト関数ＤＴＬＺ２の平均推定誤差が悪いのは、探索履歴中の個体が偏って存在していることが原因だと思われる。
図３５は実施例および比較例における平均絶対誤差と評価回数との関係を示す図である。図３５の関係はテスト関数ＺＤＴ１を用いた場合の算出結果を示す。
図３５に示すように、実施例では、比較例の２〜３倍の速度で平均絶対誤差が減少することがわかる。
図３６は実施例および比較例における平均推定誤差と評価回数との関係を示す図である。図３６の関係はテスト関数ＺＤＴ１を用いた場合の算出結果を示す。
図３６に示すように、実施例では、比較例に比べて評価の初期段階から高い推定精度が得られることがわかる。
実施例および比較例から上記第１および第２の実施の形態に係る多目的最適化装置によれば、評価回数を低減することができることがわかった。
（５）請求項の構成要素と実施の形態の各部との対応
以下、請求項の各構成要素と実施の形態の各部との対応の例について説明するが、本発明は下記の例に限定されない。
上記実施の形態では、探索履歴記憶装置３１が記憶部に相当し、適応度推定モジュール３０ａが第１の推定部に相当し、適応度推定モジュール３０ｂが第２の推定部に相当し、多目的進化型アルゴリズム部２が演算部に相当する。
また、親個体が保存用個体に相当し、子個体候補が評価用個体候補に相当し、子個体が評価用個体に相当する。
なお、第１の推定部を用いずに、実システムがシミュレーションから得た適応度を直接用いてもよい。
本発明は、実システムまたは不確実性を有するシミュレーション等の最適化対象のパラメータ、不確実性を有しないシミュレーション等の最適化対象等の種々の最適化対象を最適化するため等に利用することができる。
本発明の第１の実施の形態に係る多目的最適化装置の機能的な構成を示すブロック図である。 図１の多目的最適化装置のハードウエア構成を示すブロック図である。 最適化対象の構成の一例を示すブロック図である。 ＨＣ濃度、ＮＯx濃度およびＣＯ濃度と空燃比との関係を示す図である。 図１の多目的最適化装置の全体処理を示すフローチャートである。 図１の多目的最適化装置の全体処理を示すフローチャートである。 多目的最適化装置の各処理の具体例を示す模式図である。 多目的最適化装置の各処理の具体例を示す模式図である。 多目的最適化装置の各処理の具体例を示す模式図である。 多目的最適化装置の各処理の具体例を示す模式図である。 多目的最適化装置の各処理の具体例を示す模式図である。 多目的最適化装置の各処理の具体例を示す模式図である。 適応度推定部の適応度推定モジュールによる推定値の算出を説明するための模式図である。 正規分布に従うノイズを伴うサンプル値を示す模式図である。 不確実な適応度関数のモデルを示す模式図である。 適応度推定部の適応度推定モジュールによる推定値の算出処理を示すフローチャートである。 特定の親個体の選択から世代交代までの方法を説明するための模式図である。 子個体候補の生成処理を説明するための模式図である。 多目的進化型アルゴリズム部による特定の親個体の選択処理を示すフローチャートである。 ＵＮＤＸによる子個体候補の生成処理を示す模式図である。 多目的進化型アルゴリズム部および適応度推定モジュールによる子個体集合の生成処理を示すフローチャートである。 子個体集合の生成処理を示す模式図である。 探索履歴記憶装置に記憶された探索履歴に基づく個体の探索を示す模式図である。 親個体の選択から世代交代までの方法を説明するための模式図である。 多目的進化型アルゴリズム部による特定の親個体の選択処理を示すフローチャートである。 ｍ目的最適化問題に拡張された分布指標を示す図である。 多目的最適化装置をエンジンシミュレータに適用した例を示すブロック図である。 多目的最適化装置をモータ評価装置に適用した例を示すブロック図である。 多目的最適化装置をモータシミュレータに適用した例を示すブロック図である。 実施例および比較例の多目的最適化の条件を示す図である。 実施例および比較例において実行するベンチマーク問題を示す図である。 テスト関数ＺＤＴ１およびテスト関数ＺＤＴ２における平均絶対誤差の概念を示す模式図およびテスト関数ＤＴＬＺ２における平均絶対誤差の概念を示す模式図である。 ２目的の場合における平均絶対誤差の概念を示す模式図である。 実施例および比較例における平均絶対誤差および平均推定誤差の算出結果を示す図である。 実施例および比較例における平均絶対誤差と評価回数との関係を示す図である。 実施例および比較例における平均推定誤差と評価回数との関係を示す図である。 多目的最適化問題をエンジンの最適化に適用した例を示す図である。 パレート最適解について説明するための図である。
符号の説明
１ 多目的最適化装置
２ 多目的進化型アルゴリズム部
３ 適応度推定部
４ 出力インタフェース
５ 入力インタフェース
６ 最適化対象
１０ 使用者
３０ａ，３０ｂ 適応度推定モジュール
３１ 探索履歴記憶装置
６１ エンジン
６２ ＥＣＵ（エンジン制御ユニット）
６３ 排気ガス分析計
６４ コントローラ
６５ スロットルユニット
６６ ダイナモ
１０１ ＣＰＵ（中央演算処理装置）
１０２ ＲＯＭ（リードオンリメモリ）
１０３ ＲＡＭ（ランダムアクセスメモリ）
１０４ 入力装置
１０５ 表示装置
１０６ 外部記憶装置
１０７ 記録媒体駆動装置
１０８ 入出力インタフェース
１０９ 記録媒体
Ｃ 子個体集合
ｄl 距離
ｆ1，ｆ2 適応度関数
Ｆ 個体集合
ＨＳ 探索履歴
Ｉ１，Ｉ２，Ｉ３，Ｉ４，Ｉ５，Ｉ６，Ｉ７，Ｉ８，Ｉ９，Ｉ１０，Ｉ１１，Ｉ１２，Ｉ２１，Ｉ２２ 個体
Ｌ１，Ｌ２，Ｌ１１，Ｌ１２ 直線
Ｐ 親個体集合
Ｐ１ 第１の親個体
Ｐ２ 第２の親個体
Ｐ３ 第３の親個体
ｓ１，ｓ２，ｓ３ 長方形
ｘ1，ｘ2 パラメータ

Claims (16)

  1. 最適化対象に個体のパラメータの組を与え、複数の目的に対応する複数の適応度関数についての適応度の組を前記最適化対象から受ける多目的最適化装置であって、
    個体のパラメータの組および前記最適化対象から出力される適応度の組を保存用個体として記憶する記憶部と、
    新たな個体を評価用個体候補として生成し、生成された評価用個体候補のうち評価用個体として選択された個体のパラメータの組を前記最適化対象および前記記憶部に与えるとともに、複数組の適応度に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める演算部と、
    前記記憶部に記憶された複数の保存用個体に対応する複数組の適応度に基づいて、前記演算部により生成された評価用個体候補および前記記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める推定部とを備え、
    前記演算部は、前記複数の適応度関数の各々について前記和集合に含まれる複数の個体に対応する推定値の優劣を比較し、前記複数の適応度関数の各々について優劣の比較結果に基づいて前記和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を前記評価用個体として選択することを特徴とする多目的最適化装置。
  2. 最適化対象に個体のパラメータの組を与え、複数の目的に対応する複数の適応度関数についての適応度のサンプル値の組を前記最適化対象から受ける多目的最適化装置であって、
    個体のパラメータの組および前記最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶する記憶部と、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求める第１の推定部と、
    前記第１の推定部により求められた推定値に基づいて新たな個体を評価用個体候補として生成し、生成された評価用個体候補のうち評価用個体として選択された個体のパラメータの組を前記最適化対象および前記記憶部に与えるとともに、前記第１の推定部により求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める演算部と、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、前記演算部により生成された評価用個体候補および前記記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める第２の推定部とを備え、
    前記演算部は、前記複数の適応度関数の各々について前記和集合に含まれる複数の個体に対応する推定値の優劣を比較し、前記複数の適応度関数の各々について優劣の比較結果に基づいて前記和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を前記評価用個体として選択することを特徴とする多目的最適化装置。
  3. 前記第１の推定部は、前記記憶部に記憶された複数の保存用個体をｈlとし、注目個体ｘに対応するサンプル値の組をＦ（ｘ）とし、パラメータ空間上で注目個体から距離ｄl離れた保存用個体に対応するサンプル値の組をＦ（ｈl）とし、ｋ’を係数とし、ｌ＝１，…，Ｈとし、ｎを自然数とした場合に、
    Figure 0005019744
     で表される推定式により注目個体ｘに対応する真の適応度の推定値の組ｆ’（ｘ）を算出することを特徴とする請求項２記載の多目的最適化装置。
  4. 前記第２の推定部は、前記記憶部に記憶された複数の保存用個体をｈlとし、パラメータ空間上で評価用個体候補から距離ｄl離れた保存用個体に対応するサンプル値の組をＦ（ｈl）とし、ｋ’を係数とし、ｌ＝１，…，Ｈとし、ｎを自然数とした場合に、
    Figure 0005019744
     で表される推定式により評価用個体候補ｘに対応する真の適応度の推定値の組ｆ’（ｘ）を算出することを特徴とする請求項２または３のいずれかに記載の多目的最適化装置。
  5. 前記ｎは１であることを特徴とする請求項３または４記載の多目的最適化装置。
  6. 前記ｎは３であることを特徴とする請求項３または４記載の多目的最適化装置。
  7. 前記演算部は、適応度関数空間上で最上位ランク内で各評価用個体候補に隣接する２つの保存用個体を検出し、適応度関数ごとに前記２つの保存用個体間の数学的距離を算出し、各評価用個体候補についての複数の適応度関数における数学的距離の合計を前記分布指標として算出し、複数の評価用個体候補から算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択することを特徴とする請求項２〜６のいずれかに記載の多目的最適化装置。
  8. 前記演算部は、前記第１の推定部により求められた推定値の組に基づいて前記パレート最適個体を表示することを特徴とする請求項２〜７のいずれかに記載の多目的最適化装置。
  9. 前記演算部は、前記多目的進化型アルゴリズムとして遺伝的アルゴリズムを用いて前記評価用個体集合の個体を評価することを特徴とする請求項１〜８のいずれかに記載の多目的最適化装置。
  10. 前記最適化対象は、機器の複数の性能を評価するための評価システムを含み、前記パラメータの組は、前記評価システムのための制御用パラメータの組を含み、前記複数の適応度関数は前記評価システムの評価により得られる前記複数の性能であり、前記適応度の組は前記複数の性能の値であることを特徴とする請求項１〜９のいずれかに記載の多目的最適化装置。
  11. 前記機器はエンジンであることを特徴とする請求項１０記載の最適化装置。
  12. 前記機器はモータであることを特徴とする請求項１０記載の最適化装置。
  13. 前記評価システムは、前記パラメータの組に基づいて前記機器を制御するとともに前記機器の動作により発生される複数の性能の値をサンプル値として出力する機器評価装置であることを特徴とする請求項１０記載の機器最適化装置。
  14. 前記評価システムは、前記パラメータの組に基づいて前記機器の動作をシミュレーションすることにより複数の性能を評価し、評価された複数の性能の値をサンプル値の組として出力する機器シミュレータであることを特徴とする請求項１０記載の機器最適化装置。
  15. 最適化対象に個体のパラメータの組を与え、前記最適化対象から出力される複数の目的に対応する複数の適応度関数についての適応度のサンプル値の組に基づいてパラメータを最適化する多目的最適化方法であって、
    個体のパラメータの組および前記最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶部に記憶するステップと、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求めるステップと、
    求められた前記推定値に基づいて新たな個体を評価用個体候補として生成するステップと、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、前記生成された評価用個体候補および前記記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求めるステップと、
    前記複数の適応度関数の各々について前記和集合に含まれる複数の個体に対応する推定値の優劣を比較し、前記複数の適応度関数の各々について優劣の比較結果に基づいて前記和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択するステップと、
    前記選択された評価用個体のパラメータの組を前記最適化対象および前記記憶部に与えるとともに、求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求めるステップとを含むことを特徴とする多目的最適化方法。
  16. 最適化対象に個体のパラメータの組を与え、前記最適化対象から出力される対応する複数の適応度関数についての適応度のサンプル値の組に基づいてパラメータを最適化するコンピュータにより実行可能な多目的最適化プログラムであって、
    個体のパラメータの組および前記最適化対象から出力される適応度のサンプル値の組を保存用個体として記憶部に記憶する処理と、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて注目個体に対応する真の適応度の推定値の組を求める処理と、
    求められた前記推定値に基づいて新たな個体を評価用個体候補として生成する処理と、
    前記記憶部に記憶された複数の保存用個体に対応する複数組のサンプル値に基づいて、前記生成された評価用個体候補および前記記憶部に記憶された保存用個体の和集合に含まれる各個体に対応する真の適応度の推定値の組を求める処理と、
    前記複数の適応度関数の各々について前記和集合に含まれる複数の個体に対応する推定値の優劣を比較し、前記複数の適応度関数の各々について優劣の比較結果に基づいて前記和集合に含まれる複数の個体のランク付けを行い、最上位ランクの各評価用個体候補について最上位ランクの保存用個体の分布における疎の程度を表す分布指標を算出し、算出された分布指標が大きい順に所定数の評価用個体候補を評価用個体として選択する処理と、
    前記選択された評価用個体のパラメータの組を前記最適化対象および前記記憶部に与えるとともに、求められた複数組の推定値に基づいて評価用個体集合を多目的進化型アルゴリズムに従って評価することによりパレート最適個体集合を求める処理とを、コンピュータに実行させることを特徴とする多目的最適化プログラム。
JP2005369267A 2005-12-22 2005-12-22 多目的最適化装置、多目的最適化方法および多目的最適化プログラム Expired - Fee Related JP5019744B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2005369267A JP5019744B2 (ja) 2005-12-22 2005-12-22 多目的最適化装置、多目的最適化方法および多目的最適化プログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2005369267A JP5019744B2 (ja) 2005-12-22 2005-12-22 多目的最適化装置、多目的最適化方法および多目的最適化プログラム

Publications (2)

Publication Number Publication Date
JP2007172306A JP2007172306A (ja) 2007-07-05
JP5019744B2 true JP5019744B2 (ja) 2012-09-05

Family

ID=38298788

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2005369267A Expired - Fee Related JP5019744B2 (ja) 2005-12-22 2005-12-22 多目的最適化装置、多目的最適化方法および多目的最適化プログラム

Country Status (1)

Country Link
JP (1) JP5019744B2 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20210373888A1 (en) * 2019-08-29 2021-12-02 Huazhong University Of Science And Technology Multi-objective optimization method and system for master production plan of casting parallel workshops

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8098054B2 (en) * 2007-10-10 2012-01-17 John Alexander Verschuur Optimal load controller method and device
JP5124230B2 (ja) * 2007-10-18 2013-01-23 ヤマハ発動機株式会社 パラメトリック多目的最適化装置、パラメトリック多目的最適化方法およびパラメトリック多目的最適化プログラム
JP5156329B2 (ja) * 2007-10-18 2013-03-06 ヤマハ発動機株式会社 パラメトリック多目的最適化装置、パラメトリック多目的最適化方法およびパラメトリック多目的最適化プログラム
JP5235652B2 (ja) * 2008-12-26 2013-07-10 ヤマハ発動機株式会社 多目的最適化装置、多目的最適化方法および多目的最適化プログラム
US7987143B2 (en) * 2009-09-29 2011-07-26 Livermore Software Technology Corporation Methods and systems for multi-objective evolutionary algorithm based engineering desgin optimization
JP2015230702A (ja) * 2014-06-06 2015-12-21 株式会社リコー 処理フロー生成装置、処理フロー生成方法、処理フロー生成プログラム、情報処理装置の製造方法及び情報処理プログラムの製造方法
JP2019095819A (ja) * 2016-03-31 2019-06-20 株式会社インテック 情報処理装置及びプログラム
CN106227988A (zh) * 2016-07-11 2016-12-14 浪潮(北京)电子信息产业有限公司 一种相对误差计算优化方法
WO2019008695A1 (en) * 2017-07-05 2019-01-10 Nec Corporation DEVICE, METHOD AND PROGRAM FOR DETERMINING COST CLASSIFICATION
CN107705110B (zh) * 2017-09-29 2021-03-19 杭州乒乓智能技术股份有限公司 基于多目标的资金池支付序列优化方法及装置
CN107657311B (zh) * 2017-11-03 2019-10-29 电子科技大学 基于多目标粒子群算法的测试优选方法
CN107909194B (zh) * 2017-11-07 2021-07-13 电子科技大学 系统级测试性设计多目标优化方法
CN110728394A (zh) * 2019-04-24 2020-01-24 内蒙古电力（集团）有限责任公司内蒙古电力经济技术研究院分公司 一种考虑潮流分布的高压变电站规划方法
WO2024201834A1 (ja) * 2023-03-29 2024-10-03 富士通株式会社 演算プログラム、演算方法、および情報処理装置

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20210373888A1 (en) * 2019-08-29 2021-12-02 Huazhong University Of Science And Technology Multi-objective optimization method and system for master production plan of casting parallel workshops

Also Published As

Publication number Publication date
JP2007172306A (ja) 2007-07-05

Similar Documents

Publication Publication Date Title
JP5019744B2 (ja) 多目的最適化装置、多目的最適化方法および多目的最適化プログラム
EP1548608A2 (en) Multiobjective optimization
Nam et al. Multiobjective simulated annealing: A comparative study to evolutionary algorithms
Di Pierro et al. An investigation on preference order ranking scheme for multiobjective evolutionary optimization
Zhou et al. Multiobjective evolutionary algorithms: A survey of the state of the art
Yen et al. Dynamic multiobjective evolutionary algorithm: adaptive cell-based rank and density estimation
JP4947903B2 (ja) 最適化方法及び最適化プログラム
Mahmoodabadi et al. A novel combination of Particle Swarm Optimization and Genetic Algorithm for Pareto optimal design of a five-degree of freedom vehicle vibration model
Alajlan et al. Global path planning for mobile robots in large-scale grid environments using genetic algorithms
JP5156329B2 (ja) パラメトリック多目的最適化装置、パラメトリック多目的最適化方法およびパラメトリック多目的最適化プログラム
JP5124230B2 (ja) パラメトリック多目的最適化装置、パラメトリック多目的最適化方法およびパラメトリック多目的最適化プログラム
JP2005285090A (ja) 多目的最適化装置、多目的最適化方法および多目的最適化プログラム
Lu et al. Genetic algorithm modelling and solution of inspection path planning on a coordinate measuring machine (CMM)
US7437336B2 (en) Polyoptimizing genetic algorithm for finding multiple solutions to problems
Benavides et al. Real path planning based on genetic algorithm and Voronoi diagrams
Guo A survey of decomposition based evolutionary algorithms for many-objective optimization problems
Büche Multi-objective evolutionary optimization of gas turbine components
Hamda et al. Multi-objective evolutionary topological optimum design
CN111813669A (zh) 基于多目标群体智能的适应性随机测试用例生成方法
Ott et al. Approximate sequential optimization for informative path planning
Choy et al. Sparse actor-critic: Sparse tsallis entropy regularized reinforcement learning in a continuous action space
Li et al. A surrogate-assisted offspring generation method for expensive multi-objective optimization problems
Zakaria et al. Effects of genetic algorithm parameters on multiobjective optimization algorithm applied to system identification problem
Aliniya et al. Dynamic constrained multi-objective optimization with combination response mechanism
Duan et al. Research on welding robot path planning based on genetic algorithm

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20081104

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20111227

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20120220

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20120410

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20120514

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20120612

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20120612

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

Ref document number: 5019744

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20150622

Year of fee payment: 3

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

LAPS Cancellation because of no payment of annual fees