特開 2022-168521 最適解の演算方法、装置及びプログラム

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

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2022168521

(43)【公開日】2022-11-08

(54)【発明の名称】最適解の演算方法、装置及びプログラム

(51)【国際特許分類】

   G06N 99/00 20190101AFI20221031BHJP

【ＦＩ】

G06N99/00 180

【審査請求】未請求

【請求項の数】9

【出願形態】ＯＬ

(21)【出願番号】P 2021074043

(22)【出願日】2021-04-26

(71)【出願人】

【識別番号】000003137

【氏名又は名称】マツダ株式会社

(74)【代理人】

【識別番号】100094569

【弁理士】

【氏名又は名称】田中 伸一郎

(74)【代理人】

【識別番号】100059959

【弁理士】

【氏名又は名称】中村 稔

(74)【代理人】

【識別番号】100067013

【弁理士】

【氏名又は名称】大塚 文昭

(74)【代理人】

【識別番号】100130937

【弁理士】

【氏名又は名称】山本 泰史

(74)【代理人】

【識別番号】100162824

【弁理士】

【氏名又は名称】石崎 亮

(72)【発明者】

【氏名】山王丸 将吾

(57)【要約】

【課題】最適解を安定且つ高速に求める最適解の演算方法を提供する。
【解決手段】最適解の演算方法は、或る内部変数において目的関数を線形近似した関数が目標値を満たす解を求める第１ステップと、解に対応する目的関数値が第１ステップの目的関数値未満である場合、この解を次に第１ステップを行うべき内部変数として設定する第２ステップと、解に対応する目的関数値が第１ステップの目的関数値以上である場合、目標値を内部変数に対応する目的関数値に近付ける変更と、線形近似した関数及び変更後の目標値に基づく解の算出とを繰り返す第３ステップと、第３ステップによる解に対応する目的関数値が第１ステップの目的関数値未満になったときに、この解を次に第１ステップを行うべき内部変数として設定する第４ステップと、第２又は第４ステップで設定された内部変数から最適解を決定する第５ステップと、を有する。
【選択図】図１９

【特許請求の範囲】

【請求項1】

内部変数により規定され且つ凸関数である目的関数の最適解を演算するために、コンピュータ装置によって実行される最適解の演算方法であって、
或る内部変数と、当該内部変数に対応する前記目的関数の値である目的関数値とに基づき前記目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１ステップと、
前記第１ステップで求められた前記解に対応する目的関数値が、前記第１ステップで用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に前記第１ステップを行うべき新たな内部変数として設定する第２ステップと、
前記第１ステップで求められた前記解に対応する目的関数値が、前記第１ステップで用いられた内部変数に対応する目的関数値以上である場合に、前記目標値を前記第１ステップで用いられた目的関数値に近付けた値へと変更して、前記線形近似した関数が、変更された前記目標値を満たす解を求める第３ステップであって、前記変更された目標値から求められた前記解に対応する目的関数値が、前記第１ステップで用いられた目的関数値未満になるまで、前記目標値を前記第１ステップで用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、前記線形近似した関数が前記変更された目標値を満たす解を繰り返し求める前記第３ステップと、
前記第３ステップで求められた前記解に対応する目的関数値が、前記第１ステップで用いられた目的関数値未満になったときに、当該解を、次に前記第１ステップを行うべき新たな内部変数として設定する第４ステップと、
前記第２又は第４ステップで設定された内部変数に対応する目的関数値と前記目標値との差分が所定値未満になるまで、前記第１乃至第４ステップを繰り返し行わせ、前記差分が前記所定値未満になったときに、このときに用いられていた内部変数を前記最適解として決定する第５ステップと、
を有することを特徴とする最適解の演算方法。

【請求項2】

前記第３ステップでは、前記第１ステップで用いられた内部変数に対応する目的関数値と前記目標値との中央値を求め、この中央値によって前記目標値を変更する、請求項１に記載の最適解の演算方法。

【請求項3】

更に、前記第１又は第３ステップで求められた前記解に対応する目的関数値が前記目標値未満である場合、当該目標値を当該目的関数値よりも小さい値へと変更する第６ステップを有する、請求項１又は２に記載の最適解の演算方法。

【請求項4】

前記第６ステップでは、前記第１又は第３ステップで求められた前記解に対応する目的関数値が前記目標値未満である場合、当該目的関数値と当該目標値の変更前に用いられていた目標値との間の中央値を求め、この中央値によって前記目標値を変更する、請求項３に記載の最適解の演算方法。

【請求項5】

前記第５ステップでは、前記第２又は第４ステップで設定された内部変数に対応する目的関数値と前記目標値との差分が前記所定値未満で、且つ、前記第２又は第４ステップで設定された内部変数に対応する目的関数値と、前記目標値の変更前に用いられていた目標値との差分が所定値未満となったときに、このときに用いられていた内部変数を前記最適解として決定する、請求項１乃至４のいずれか一項に記載の最適解の演算方法。

【請求項6】

前記最適解の演算方法は、少なくとも一部分が固定され且つ荷重が少なくとも一部分に付与される部材について、当該部材を構成する材料の最適な充填率分布を求めるために用いられ、
前記内部変数は、前記充填率分布に対応するパラメータであり、前記目的関数は、前記荷重による前記部材の変位量に応じたパラメータにより規定される、
請求項１乃至５のいずれか一項に記載の最適解の演算方法。

【請求項7】

前記第１ステップで求められた前記解から前記部材の体積を求め、この体積が所定の制約値以上である場合には、前記体積が前記制約値未満となるまで、前記第２及び第３ステップを行わずに、前記第１ステップを繰り返し行う、請求項６に記載の最適解の演算方法。

【請求項8】

内部変数により規定され且つ凸関数である目的関数の最適解を演算する最適解の演算装置であって、
或る内部変数と、当該内部変数に対応する前記目的関数の値である目的関数値とに基づき前記目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１手段と、
前記第１手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に前記第１手段を行うべき新たな内部変数として設定する第２手段と、
前記第１手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた内部変数に対応する目的関数値以上である場合に、前記目標値を前記第１手段で用いられた目的関数値に近付けた値へと変更して、前記線形近似した関数が、変更された前記目標値を満たす解を求める第３手段であって、前記変更された目標値から求められた前記解に対応する目的関数値が、前記第１手段で用いられた目的関数値未満になるまで、前記目標値を前記第１手段で用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、前記線形近似した関数が前記変更された目標値を満たす解を繰り返し求める前記第３手段と、
前記第３手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた目的関数値未満になったときに、当該解を、次に前記第１手段を行うべき新たな内部変数として設定する第４手段と、
前記第２又は第４手段で設定された内部変数に対応する目的関数値と前記目標値との差分が所定値未満になるまで、前記第１乃至第４手段を繰り返し行わせ、前記差分が前記所定値未満になったときに、このときに用いられていた内部変数を前記最適解として決定する第５手段と、
を有することを特徴とする最適解の演算装置。

【請求項9】

内部変数により規定され且つ凸関数である目的関数の最適解を演算するために、コンピュータ装置によって実行される最適解の演算プログラムであって、
前記コンピュータ装置を、
或る内部変数と、当該内部変数に対応する前記目的関数の値である目的関数値とに基づき前記目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１手段と、
前記第１手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に前記第１手段を行うべき新たな内部変数として設定する第２手段と、
前記第１手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた内部変数に対応する目的関数値以上である場合に、前記目標値を前記第１手段で用いられた目的関数値に近付けた値へと変更して、前記線形近似した関数が、変更された前記目標値を満たす解を求める第３手段であって、前記変更された目標値から求められた前記解に対応する目的関数値が、前記第１手段で用いられた目的関数値未満になるまで、前記目標値を前記第１手段で用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、前記線形近似した関数が前記変更された目標値を満たす解を繰り返し求める前記第３手段と、
前記第３手段で求められた前記解に対応する目的関数値が、前記第１手段で用いられた目的関数値未満になったときに、当該解を、次に前記第１手段を行うべき新たな内部変数として設定する第４手段と、
前記第２又は第４手段で設定された内部変数に対応する目的関数値と前記目標値との差分が所定値未満になるまで、前記第１乃至第４手段を繰り返し行わせ、前記差分が前記所定値未満になったときに、このときに用いられていた内部変数を前記最適解として決定する第５手段と、
として機能させることを特徴とする最適解の演算プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、目的関数の最適解を演算する最適解の演算方法、装置及びプログラムに関する。

【背景技術】

【0002】

従来から、所定の内部変数により規定された目的関数の最適解を演算する技術が知られている。１つの例では、この技術は、与えられた部材の設計領域内（材料分布範囲内）において、部材を構成する材料の最適な充填率分布（密度分布）を求めるための位相最適化（トポロジー最適化）に適用されている。この例では、目的関数は、荷重による部材の変位量に応じたパラメータ（例えばコンプライアンス）により規定され、この目的関数を規定する内部変数には、充填率分布に対応するパラメータが適用される。

【0003】

また、目的関数の最適解を求めるに当たって、最適化問題において目的関数の勾配に関する情報を解の探索に用いる勾配法や、方程式系を数値計算によって解くための反復法によるアルゴリズムであるニュートン法（ニュートン・ラフソン法）などが用いられている。例えば、特許文献１には、勾配法を用いて複数回計算を行うことで、複数個のパラメータの値を同時に選択して、パラメータの最適化の高速化を図る技術が開示されている。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開２０２０－１８１３１８号公報

【発明の概要】

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

【0005】

しかしながら、勾配法は計算時間が膨大にかかるという問題があり、ニュートン法は計算が不安定になる場合があるという問題があった。そこで、本願発明者は、このような問題を解決すべく鋭意研究を重ねた結果、目的関数の最適解を安定且つ高速に求めることができる演算方法を見出した。

【0006】

本発明は、上述した従来技術の問題点を解決するためになされたものであり、内部変数により規定され且つ凸関数である目的関数の最適解を演算する最適解の演算方法、装置及びプログラムにおいて、最適解を安定且つ高速に求めることを目的とする。

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

【0007】

上記の目的を達成するために、本発明は、内部変数により規定され且つ凸関数である目的関数の最適解を演算するために、コンピュータ装置によって実行される最適解の演算方法であって、或る内部変数と、当該内部変数に対応する目的関数の値である目的関数値とに基づき目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１ステップと、第１ステップで求められた解に対応する目的関数値が、第１ステップで用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に第１ステップを行うべき新たな内部変数として設定する第２ステップと、第１ステップで求められた解に対応する目的関数値が、第１ステップで用いられた内部変数に対応する目的関数値以上である場合に、目標値を第１ステップで用いられた目的関数値に近付けた値へと変更して、線形近似した関数が、変更された目標値を満たす解を求める第３ステップであって、変更された目標値から求められた解に対応する目的関数値が、第１ステップで用いられた目的関数値未満になるまで、目標値を第１ステップで用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、線形近似した関数が変更された目標値を満たす解を繰り返し求める第３ステップと、第３ステップで求められた解に対応する目的関数値が、第１ステップで用いられた目的関数値未満になったときに、当該解を、次に第１ステップを行うべき新たな内部変数として設定する第４ステップと、第２又は第４ステップで設定された内部変数に対応する目的関数値と目標値との差分が所定値未満になるまで、第１乃至第４ステップを繰り返し行わせ、差分が所定値未満になったときに、このときに用いられていた内部変数を最適解として決定する第５ステップと、を有することを特徴とする。

【0008】

このように構成された本発明では、第１ステップで求められた解に対応する目的関数値が、第１ステップで用いられた目的関数値以上である場合に、目標値から求められた解に対応する目的関数値が第１ステップで用いられた目的関数値未満になるまで、目標値を第１ステップで用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、線形近似した関数（一次関数）が変更された目標値を満たす解を繰り返し求めて、解に対応する目的関数値が第１ステップで用いられた目的関数値未満になったときに、当該解を、次に第１ステップを行うべき新たな内部変数として設定する。
このような本発明によれば、悪い解も受け入れるしかなかった従来の手法（特にニュートン法）と異なり、最適解の導出のために用いられる制約値としての目標値を変動させことにより、悪い解を切り捨てて良い解を確実に求めることができるので、安定した計算を行うことが可能となる。また、本発明によれば、良い解を見つけるまでの繰り返しの処理において、目的関数を線形近似した関数及び第１ステップで用いられた内部変数に対応する目的関数値を使い回すことができるので、計算ステップ数の増大を効果的に抑制することができる。以上より、本発明によれば、最適解を安定且つ高速に求めることができる。

【0009】

本発明において、好ましくは、第３ステップでは、第１ステップで用いられた内部変数に対応する目的関数値と目標値との中央値を求め、この中央値によって目標値を変更する。
このように構成された本発明によれば、第１ステップで用いられた内部変数に対応する目的関数値及び元の目標値に応じた的確な値へと、目標値を更新することができる。

【0010】

本発明において、好ましくは、更に、第１又は第３ステップで求められた解に対応する目的関数値が目標値未満である場合、当該目標値を当該目的関数値よりも小さい値へと変更する第６ステップを有する。
このように構成された本発明によれば、最適解の演算過程において順次求められる解に対応する目的関数値に応じて、目標値を的確な値へと更新することができる。

【0011】

本発明において、好ましくは、第６ステップでは、第１又は第３ステップで求められた解に対応する目的関数値が目標値未満である場合、当該目的関数値と当該目標値の変更前に用いられていた目標値との間の中央値を求め、この中央値によって目標値を変更する。
このように構成された本発明によれば、解に対応する目的関数値及び前回用いていた目標値に応じた的確な値へと、仮目標値を更新することができる。

【0012】

本発明において、好ましくは、第５ステップでは、第２又は第４ステップで設定された内部変数に対応する目的関数値と目標値との差分が所定値未満で、且つ、第２又は第４ステップで設定された内部変数に対応する目的関数値と、目標値の変更前に用いられていた目標値との差分が所定値未満となったときに、このときに用いられていた内部変数を最適解として決定する。
このように構成された本発明によれば、目的関数値と現在の目標値との差分に加えて、当該目的関数値と前回用いていた目標値との差分を用いることで、解の収束判定を精度良く行うことができる。つまり、より正確な最適解を求めることができる。

【0013】

本発明において、好ましくは、最適解の演算方法は、少なくとも一部分が固定され且つ荷重が少なくとも一部分に付与される部材について、当該部材を構成する材料の最適な充填率分布を求めるために用いられ、内部変数は、充填率分布に対応するパラメータであり、目的関数は、荷重による部材の変位量に応じたパラメータにより規定される。
このように構成された本発明によれば、上述した演算方法を位相最適化（トポロジー最適化）に適用することで、材料の最適な充填率分布（密度分布）を安定且つ高速に求めることができる。

【0014】

本発明において、好ましくは、第１ステップで求められた解から部材の体積を求め、この体積が所定の制約値以上である場合には、体積が制約値未満となるまで、第２及び第３ステップを行わずに、第１ステップを繰り返し行う。
このように構成された本発明によれば、確実に体積制約を満たす解（最適な充填率分布）を求めることができる。

【0015】

他の観点では、上記の目的を達成するために、本発明は、内部変数により規定され且つ凸関数である目的関数の最適解を演算する最適解の演算装置であって、或る内部変数と、当該内部変数に対応する目的関数の値である目的関数値とに基づき目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１手段と、第１手段で求められた解に対応する目的関数値が、第１手段で用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に第１手段を行うべき新たな内部変数として設定する第２手段と、第１手段で求められた解に対応する目的関数値が、第１手段で用いられた内部変数に対応する目的関数値以上である場合に、目標値を第１手段で用いられた目的関数値に近付けた値へと変更して、線形近似した関数が、変更された目標値を満たす解を求める第３手段であって、変更された目標値から求められた解に対応する目的関数値が、第１手段で用いられた目的関数値未満になるまで、目標値を第１手段で用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、線形近似した関数が変更された目標値を満たす解を繰り返し求める第３手段と、第３手段で求められた解に対応する目的関数値が、第１手段で用いられた目的関数値未満になったときに、当該解を、次に第１手段を行うべき新たな内部変数として設定する第４手段と、第２又は第４手段で設定された内部変数に対応する目的関数値と目標値との差分が所定値未満になるまで、第１乃至第４手段を繰り返し行わせ、差分が所定値未満になったときに、このときに用いられていた内部変数を最適解として決定する第５手段と、を有することを特徴とする。
このように構成された本発明によっても、最適解を安定且つ高速に求めることができる。

【0016】

更に他の観点では、上記の目的を達成するために、本発明は、内部変数により規定され且つ凸関数である目的関数の最適解を演算するために、コンピュータ装置によって実行される最適解の演算プログラムであって、コンピュータ装置を、或る内部変数と、当該内部変数に対応する目的関数の値である目的関数値とに基づき目的関数を線形近似した関数が、所定の目標値を満たす解を求める第１手段と、第１手段で求められた解に対応する目的関数値が、第１手段で用いられた内部変数に対応する目的関数値未満である場合に、当該解を、次に第１手段を行うべき新たな内部変数として設定する第２手段と、第１手段で求められた解に対応する目的関数値が、第１手段で用いられた内部変数に対応する目的関数値以上である場合に、目標値を第１手段で用いられた目的関数値に近付けた値へと変更して、線形近似した関数が、変更された目標値を満たす解を求める第３手段であって、変更された目標値から求められた解に対応する目的関数値が、第１手段で用いられた目的関数値未満になるまで、目標値を第１手段で用いられた目的関数値に徐々に近付ける変更を繰り返し行うと共に、線形近似した関数が変更された目標値を満たす解を繰り返し求める第３手段と、第３手段で求められた解に対応する目的関数値が、第１手段で用いられた目的関数値未満になったときに、当該解を、次に第１手段を行うべき新たな内部変数として設定する第４手段と、第２又は第４手段で設定された内部変数に対応する目的関数値と目標値との差分が所定値未満になるまで、第１乃至第４手段を繰り返し行わせ、差分が所定値未満になったときに、このときに用いられていた内部変数を最適解として決定する第５手段と、として機能させることを特徴とする。
このように構成された本発明によっても、最適解を安定且つ高速に求めることができる。

【発明の効果】

【0017】

本発明によれば、内部変数により規定され且つ凸関数である目的関数の最適解を演算する最適解の演算方法、装置及びプログラムにおいて、最適解を安定且つ高速に求めることができる。

【図面の簡単な説明】

【0018】

【図1】本発明の実施形態に係る最適解の演算方法の実行主体の一例であるコンピュータ装置の概略構成図である。

【図2】本発明の実施形態の基本概念についての説明図である。

【図3】位相最適化方法の基本的内容についての説明図である。

【図4】位相最適化方法における最適解についての説明図である。

【図5】位相最適化方法でのコンプライアンスの特性を表す曲線についての説明図である。

【図6】位相最適化方法における最適解を求めるための基本的考え方についての説明図である。

【図7】位相最適化方法において従来の勾配法により最適解を求める方法についての説明図である。

【図8】位相最適化方法において従来の勾配法を適用した場合の結果の一例を示す。

【図9】位相最適化方法において従来のニュートン法により最適解を求める方法についての説明図である。

【図10】位相最適化方法において従来のニュートン法を適用した場合の結果の一例を示す。

【図11】位相最適化方法において従来のニュートン法を適用した場合の結果の他の例を示す。

【図12】ニュートン法による計算が不安定になった１つの原因についての説明図である。

【図13】ニュートン法による計算が不安定になった別の原因についての説明図である。

【図14】本発明の実施形態に係る制約変動法で用いる試行解についての説明図である。

【図15】本発明の実施形態に係る制約変動法において行われる第１処理についての説明図である。

【図16】本発明の実施形態に係る制約変動法において行われる第２処理についての説明図である。

【図17】本発明の実施形態に係る制約変動法において行われる第３処理についての説明図である。

【図18】本発明の実施形態に係る位相最適化方法の全体処理を示すフローチャートであり

【図19】本発明の実施形態に係る制約変動法処理を示すフローチャートである。

【図20】本発明の実施形態に係る制約変動法による結果の一例を示す。

【図21】本発明の実施形態に係る制約変動法による結果の他の例を示す。

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

【0019】

以下、添付図面を参照して、本発明の実施形態による最適解の演算方法、装置及びプログラムについて説明する。

【0020】

＜コンピュータ装置＞
まず、図１を参照して、本発明の実施形態に係る最適解の演算方法の実行主体の一例であるコンピュータ装置について説明する。図１に示すように、コンピュータ装置１０は、利用者などにより情報が入力される入力装置１と、種々の情報を処理する処理装置３と、情報を出力する出力装置５と、を有する。

【0021】

入力装置１は、例えばマウスやキーボードやタッチパネルやマイクなどであり、出力装置５は、表示装置やスピーカなどである。処理装置３は、プログラムを実行する中央演算処理装置（Central Processing Unit：ＣＰＵ）としての１以上のマイクロプロセッサ３ａと、例えばＲＡＭ（Random Access Memory）やＲＯＭ（Read Only Memory）やハードディスクなどにより構成されてプログラム及びデータを格納するメモリ３ｂと、を有する。

【0022】

本実施形態に係る最適解の演算方法は、コンピュータ装置１０の処理装置３によって実行される。具体的には、処理装置３のメモリ３ｂには、この最適解の演算方法に対応するプログラム（本発明の最適解の演算プログラム）が記憶されており、処理装置３のマイクロプロセッサ３ａが、このプログラムをメモリ３ｂから読み出して実行する。これにより、コンピュータ装置１０（特に処理装置３）は、本発明の最適解の演算装置として機能する。具体的には、処理装置３は、本発明における第１乃至第５手段として機能する。

【0023】

＜基本概念＞
次に、図２を参照して、本発明の実施形態の基本概念について説明する。図２は、目的関数の一例を示している。具体的には、図２は、横軸に、目的関数を規定する内部変数（例えば設計変数など）を示し、縦軸に、この内部変数に応じた目的関数（つまり内部変数に対応する目的関数の値（目的関数値））を示している。図２に示すように、目的関数は、下に凸の形状を有しており、最小値を含んでいる。本実施形態による最適解の演算方法では、このような目的関数の最小値に対応する最適解を安定且つ高速に求めることを図っている。

【0024】

＜位相最適化方法＞
次に、本実施形態による最適解の演算方法について具体的に説明する。ここでは、図３乃至図１９を参照して、本実施形態による最適解の演算方法を位相最適化方法に適用した例について説明する。

【0025】

まず、図３は、位相最適化方法の基本的内容についての説明図である。図３（Ａ）に示すように、位相最適化方法は、少なくとも一部分が固定され且つ荷重が少なくとも一部分に付与される部材２０について、当該部材２０を構成する材料の最適な充填率分布を求めるものである。ここでは、部材２０の左辺が固定され、部材２０の右辺の中央部付近に荷重が付与される場合を例示している。位相最適化方法においては、目的関数は、荷重による部材２０の変位量を示すコンプライアンスに相当し、このコンプライアンスは、内部変数としての充填率分布によって規定される。特に、コンプライアンスは、「荷重（外力）×変位」により規定され、この式中の変位が充填率分布に応じて変わる。

【0026】

本実施形態では、このような位相最適化方法に関して、空間の離散化モデル及び境界条件（例えば有限要素法のメッシュ）などを用いて、部材２０の体積制約を課した上で、コンプライアンスが最小となる充填率分布（つまり最適な充填率分布）を最適解として求める。図３（Ｂ）は、３０％の体積制約（部材２０の全部が充填されているときの充填率を１００％とする）が適用されたときに求められた最適な充填率分布の一例を示す。なお、図３（Ｂ）は、グレースケールで表された色によって充填率分布を表現している（以下同様とする）。

【0027】

次いで、図４は、位相最適化方法における最適解の基本的内容についての説明図である。図４は、横軸に充填率分布を示し、縦軸に目的関数としてのコンプライアンスを示している。充填率を「ρ」とし、充填率ρに応じた剛性を「Ｋ（ρ）」とし、変位を「ｕ」とし、荷重を「ｂ」とすると、剛性方程式は「Ｋ（ρ）ｕ＝ｂ」で表され、また、目的関数としてのコンプライアンスを「ｆ₀」とすると、コンプライアンスｆ₀は「ｆ₀＝ｂ・ｕ」と表される。図４に示すように、充填率ρ分布を変化させると、種々のコンプライアンスｆ₀が求められ、体積制約を満たすものの中でコンプライアンスｆ₀が最小になるときの充填率ρ分布が最適解となる。全ての充填率ρ分布のパターンについてコンプライアンスｆ₀を計算する（つまり総当たり計算を行う）ことにより、最適解を求めることができるが、実際は有限要素法の要素数が膨大にあるので、この方法は現実的に不可能である。

【0028】

ここで、１つの充填率ρ分布には、有限要素法などにより部材２０に対して規定された複数の節点数の成分（節点値）が含まれるものとする。また、「ｕ」は、剛性方程式を解いて得られた、部材２０の全体変位ベクトルの節点値（次元×節点数の成分）であり、「ｂ」は、部材２０に付与される全体外力ベクトルの節点値（次元×節点数の成分）であり、「Ｋ（ρ）」は、充填率ρに依存する剛性行列（(次元×節点数）×（次元×節点数）)となる。

【0029】

次いで、図５は、位相最適化方法でのコンプライアンスｆ₀の特性を表す曲線についての説明図である。図５も、横軸に充填率ρ分布を示し、縦軸にコンプライアンスｆ₀を示している。図５に示すように、目的関数としてのコンプライアンスｆ₀の特性を表す曲線がわかれば、上記したような最適解を求めることができ、そして、充填率ρ分布をわずかに変化（微小変化）させて多数の点を計算すれば、このコンプライアンスｆ₀の特性を表す曲線がわかると言える。

【0030】

ところで、充填率ρそのものを設計変数（内部変数）として用いると、ρがマイナスになったり、上限値（１００％）を超えたりするような異常な解が出てしまうことがある。この場合、充填率ρを「０≦ρ≦１．０」に制限するために、ムーブリミットなどの特殊な処理を行うことも可能であるが、そうすると解の性能が悪くなってしまう。したがって、本実施形態では、所定の変数θを使って充填率ρを間接的に表現することとした。つまり、充填率ρの代わりに、変数θを各節点の設計変数とする。具体的には、「ρ（θ）＝１／２×ｔａｎｈθ＋１／２」というシグモイド関数を用いて、充填率ρを設計変数θによって規定する。こうすることで、設計変数θには範囲制限がないので（－∞＜θ＜∞）、ムーブリミットを用いる必要がなくなる。すなわち、ムーブリミットを用いなくても、「０≦ρ（θ）≦１．０」が自動的に成立することとなる。

【0031】

このような設計変数θを用いると、図５に示すように、コンプライアンスｆ₀は、充填率ρ分布の代わりに設計変数θ分布（以下では単に「θ分布」とも呼ぶ。）により規定され、また、コンプライアンスｆ₀に関係する剛性方程式は「Ｋ（ρ（θ））ｕ＝ｂ」と表される。この場合、コンプライアンスｆ₀が最小となる設計変数θを求めて、この設計変数θに対応する充填率ρを求めることにより、部材２０の最適形状がわかる。なお、設計変数θも、充填率ρ分布と同様に、有限要素法などにより規定された複数の節点数の成分（節点値）を含むものである。

【0032】

次いで、図６は、位相最適化方法における最適解を求めるための基本的考え方についての説明図である。図６は、横軸に設計変数θ分布を示し、縦軸にコンプライアンスｆ₀を示している。上述したように、有限要素法の要素数が膨大にあるので、コンプライアンスｆ₀を総当たり計算することは現実的に不可能である。一方で、図６に示すように、或る設計変数θ分布に関する「高さｆ₀」及び「勾配（傾き）ｇ₀」のみを用いて、コンプライアンスｆ₀が改善するようにθ分布を更新していくと、θ分布の最適解を求めることができると言える。

【0033】

次いで、図７は、位相最適化方法において従来の勾配法により最適解を求める方法についての説明図である。図７も、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示している。勾配法は、コンプライアンスｆ₀（目的関数）の特性を山とみなして、この山に沿ってボールを転がして最適解を見つけようとする手法である。具体的には、この勾配法では、ボールは、或る高さｆ₀での勾配ｇ₀に応じた変化ベクトルΔθ（降下ベクトル）だけ転がす。この場合、勾配ｇ₀は「ｇ₀＝∂ｆ₀／∂θ」と表され、変化ベクトルΔθは、この勾配ｇ₀に対して所定の調整定数Ｃ_aを適用することで「Δθ＝－ｇ₀／Ｃ_a」と表される。調整定数Ｃ_aは変化ベクトルΔθを調整するパラメータであり、調整定数Ｃ_aが大きいほど、変化ベクトルΔθの大きさ（絶対値）が小さくなり、調整定数Ｃ_aが小さいほど、変化ベクトルΔθの大きさ（絶対値）が大きくなる。

【0034】

次いで、図８は、位相最適化方法において従来の勾配法を適用した場合の結果の一例を示す。図８は、横軸に計算ステップ（計算時間に対応する）を示し、縦軸にコンプライアンスｆ₀（ｍＪ）を示している。グラフＧ１１は、比較的大きな調整定数Ｃ_aを用いた場合の結果を示し、グラフＧ１３は、比較的小さな調整定数Ｃ_aを用いた場合の結果を示し、グラフＧ１２は、これらグラフＧ１１、Ｇ１３の間の調整定数Ｃ_aを用いた場合の結果を示している。図８に示すように、勾配法によれば、最小のコンプライアンスｆ₀（最適解）が求まるまでに要する計算ステップが多くなることがわかる。特に、比較的大きな調整定数Ｃ_aを用いた場合（グラフＧ１１）には計算ステップ数が多くなり、これよりも小さな調整定数Ｃ_aを用いた場合（グラフＧ１２）には計算ステップ数が少なくなるが、更に小さな調整定数Ｃ_aを用いた場合（グラフＧ１３）には解が発散していることがわかる。

【0035】

次いで、図９は、位相最適化方法において従来のニュートン法により最適解を求める方法についての説明図である。図９も、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示している。ニュートン法は、或るθ分布の点において目的関数（コンプライアンスｆ₀）を二次関数で近似し、この近似した二次関数の極値点によって、次に二次関数による近似を行うべきθ分布の値を更新する、という処理を繰り返すことで、最適解を求める手法である。

【0036】

図９に示すニュートン法の例では、まず、θ分布の値θ１（以下ではθ分布の値を適宜「θ分布値」と呼ぶ。）、及び、当該θ分布値θ１に対応するコンプライアンスｆ₀の値ｆ₀（θ１）（以下ではコンプライアンスｆ₀の値（目的関数値）を適宜「コンプライアンス値」と呼ぶ。）に基づき、目的関数を近似した二次関数ＱＦ１を得ることで、この二次関数ＱＦ１の極値点に対応するθ分布値θ２を求める。次いで、このθ分布値θ２及び当該θ分布値θ２に対応するコンプライアンス値ｆ₀（θ２）に基づき目的関数を近似した二次関数ＱＦ２を得ることで、この二次関数ＱＦ２の極値点に対応するθ分布値θ３を求め、このθ分布値θ３及び当該θ分布値θ３に対応するコンプライアンス値ｆ₀（θ３）に基づき目的関数を近似した二次関数ＱＦ３を得ることで、この二次関数ＱＦ３の極値点に対応するθ分布値θ４を求める、という処理が繰り返される。なお、各コンプライアンス値は、「Ｋ（ρ）ｕ＝ｂ」という剛性方程式を用いて、「ｆ₀＝ｂ・ｕ」というコンプライアンスｆ₀の式から求められる。このようなニュートン法によれば、上述した勾配法よりも、最適解を求めるまでの計算ステップを短縮することができる。

【0037】

次いで、図１０は、位相最適化方法において従来のニュートン法を適用した場合の結果の一例を示す。図１０も、横軸に計算ステップを示し、縦軸にコンプライアンスｆ₀（ｍＪ）を示している。グラフＧ１２は、図８に示したものと同じ勾配法による結果を示し、グラフＧ２は、ニュートン法による結果を示している。図１０に示す例では、ニュートン法によれば、勾配法と比較して、最小のコンプライアンスｆ₀（最適解）が求まるまでに要する計算ステップが少なくなることがわかる。

【0038】

次いで、図１１は、位相最適化方法において従来のニュートン法を適用した場合の結果の他の例を示す。図１１も、横軸に計算ステップを示し、縦軸にコンプライアンスｆ₀（ｍＪ）を示している。図１１に示す例では（グラフＧ３参照）、ニュートン法による計算が非常に不安定になっていることがわかる。この例では解が最終的に収束しているが、解が収束せずに発散する場合もある。

【0039】

次いで、図１２及び図１３を参照して、上記したようにニュートン法による計算が不安定になった原因について説明する。図１２は、ニュートン法による計算が不安定になった１つの原因の説明図である。図１２において、左図は、θ分布とコンプライアンスｆ₀（目的関数）のグラフを示し、右図は、計算のステップとコンプライアンス値のグラフを示している。図１２に示す例でも、図１１で述べたのと同様のニュートン法の手順により計算が順次行われる。この場合、ステップＳ１においては、コンプライアンス値ｆ₀（θ１）が得られ、ステップＳ２においては、このコンプライアンス値ｆ₀（θ１）よりも小さなコンプライアンス値ｆ₀（θ２）が得られ、ステップＳ３において、このコンプライアンス値ｆ₀（θ２）よりも大きなコンプライアンス値ｆ₀（θ３）が得られる。したがって、解の振動が発生していると言える。こうなったのは、コンプライアンスｆ₀（目的関数）の非線形性が強いからであると考えられる。非線形性が強いコンプライアンスｆ₀においては、勾配法のほうが、ニュートン法よりも、解が安定して収束すると考えられる。

【0040】

次いで、図１３は、ニュートン法による計算が不安定になった別の原因の説明図である。図１３において、左図は、θ分布とコンプライアンスｆ₀（目的関数）のグラフを示し、右図は、計算のステップとコンプライアンス値のグラフを示している。図１３に示す例でも、図１１で述べたのと同様のニュートン法の手順により計算が順次行われる。この場合、ステップＳ１においては、コンプライアンス値ｆ₀（θ１）が得られ、ステップＳ２においては、このコンプライアンス値ｆ₀（θ１）よりも小さなコンプライアンス値ｆ₀（θ２）が得られ、ステップＳ３において、このコンプライアンス値ｆ₀（θ２）よりも大きなコンプライアンス値ｆ₀（θ３）が得られ、ステップＳ４においては、このコンプライアンス値ｆ₀（θ３）よりも小さなコンプライアンス値ｆ₀（θ４）が得られる。したがって、解の振動が発生していると言える。こうなったのは、目的関数を近似した二次関数ＱＦ２の凸の向きが変わった（下に凸から上に凸に変わった）からであると考えられる。これは、調整定数（強圧性）を調整することで緩和可能だが、この定数を決めるのに試行錯誤が必要となる。このような場合にも、勾配法のほうが、ニュートン法よりも、解が安定して収束すると考えられる。

【0041】

以上より、ニュートン法では、目的関数の形状が整ったものでないと、計算が安定せずに、最適解が得られにくいと言える。一方、勾配法は、ニュートン法と比べて、計算は安定しているが、計算時間がかかる。したがって、本実施形態では、勾配法及びニュートン法とは別の方法（本実施形態に係る最適解の演算方法を以下では適宜「制約変動法」と呼ぶ。）を位相最適化方法に適用することで、目的関数の最適解を安定且つ高速に求めることを図っている。以下では、図１４乃至図１９を参照して、本実施形態に係る制約変動法の内容について具体的に説明する。

【0042】

図１４は、本実施形態に係る制約変動法で用いる試行解についての説明図である。図１４は、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示している。本実施形態では、コンピュータ装置１０の処理装置３は、或るθ分布値θ_k（或る計算ステップｋでのθ分布値）において目的関数を線形近似した一次関数ＬＦ１が仮目標値ｆ_0cを満たす解を試行解θ_Tryとして求める。具体的には、一次関数ＬＦ１は、θ分布値θ_k及び当該θ分布値θ_kに対応するコンプライアンス値ｆ₀（θ_k）において規定される線形近似関数であり、試行解θ_Tryは、この一次関数ＬＦ１が、「コンプライアンス値ｆ₀＝仮目標値ｆ_0c」となる水平な直線に交わる点のθ分布値である。基本的には、この試行解θ_Tryは、目的関数を線形近似した一次関数が仮目標値ｆ_0cを満たす試行解θ_Tryを次に求めるべき、新たなθ分布値として適用される。したがって、本実施形態では、このような手順で試行解θ_Tryを繰り返し求めて更新していくことで、θ分布値の最適解、つまりコンプライアンス値ｆ₀が最小となるときのθ分布値を最終的に求めるようにしている。なお、仮目標値ｆ_0cは、最適解を求めるために設定されるコンプライアンス値の仮の目標値であり、詳細は後述するが、逐次変更（更新）されるものである。仮目標値ｆ_0cの初期値は、例えば０（変位量０に相当する）に設定される。

【0043】

次に、処理装置３による試行解θ_Tryの具体的な求め方について説明する。まず、θ分布値θ_kと試行解θ_Tryとの差分Δθは「｛Δθ｝＝｛θ_k｝－｛θ_Try｝」と表される。これら｛Δθ｝、｛θ_k｝、｛θ_Try｝は、それぞれ、各パラメータが節点数の成分を含むことを意味する（以下同様とする）。Δθは、制約条件が何もない場合、上記の勾配法と同様に、目的関数ｆ₀の勾配｛ｇ₀｝を用いて、「｛Δθ｝＝－｛ｇ₀｝／Ｃ_a」という式で表される。この式は、関数ｑ（｛Δθ｝）の最小化を示す式「ｑ（｛Δθ｝）＝（｛Δθ｝・｛Δθ｝）／２＋（｛ｇ₀｝・｛Δθ｝）／Ｃ_a」により置き換えられる。しかしながら、本実施形態では体積制約などの制約条件を用いるため、制約条件を満たしながら関数ｑ（｛Δθ｝）を小さくするために、Lagrangeの未定乗数法を用いる。

【0044】

ここで、「ｆ_i（｛θ_k｝＋｛Δθ｝）＝ｆ_ic」を考える（「ｉ」は制約関数の番号であり、「ｆ_ic」は制約値である）。この式の左辺を線形化すると、以下の式が得られる（Taylor展開の第１項までを採用し線形化した制約条件にする）。

この式の第２項は「ｇ_i（｛θ_k｝」となるため、当該式を変形すると、「ｆ_i（｛θ_k｝）＋ｇ_i（｛θ_k｝）・｛Δθ｝－ｆ_ic＝０」となる。

【0045】

他方で、Lagrangeの未定乗数法の公式は、以下の通りである。「λ_i」は、Lagrangeの未定乗数であり、制約条件Ｇ_i（｛ｘ｝）は「＝０」の形を用いるものとする。

【0046】

制約条件を満たす第１項の最小化は、以下の連立方程式を解けばよい。そして、未定乗数λ_iを求めた後、｛ｘ｝を求めればよい。

【0047】

上述したように、最小化したい関数は「（｛Δθ｝・｛Δθ｝）／２＋（｛ｇ₀｝・｛Δθ｝）／Ｃ_a」であり、線形化した制約条件は「ｆ_i（｛θ_k｝）＋ｇ_i（｛θ_k｝）・｛Δθ｝－ｆ_ic＝０」であるので、これらを未定乗数法に適用すると、以下の式となる。

【0048】

よって、｛Δθ｝を求める式は以下の通りである。本実施形態に係る制約変動法では、目的関数ｆ₀を線形化して制約関数化するので（仮目標を満たすようにするので）、下記式中の「｛ｇ₀｝／Ｃ_a」は消えることとなる。

【0049】

このような式を本実施形態に係る制約変動法に適用することで、以下の連立方程式（１）～（４）より、試行解θ_Tryを求めることができる。式（１）は上述した通りで、式（２）～（４）は、上式から得られたものである。式（２）中の「λ₀」、「λ₁」は、Lagrangeの未定乗数である。式（３）は、仮目標値ｆ_0cを用いて目的関数ｆ₀を制約条件化した式であり、特に左辺は線形化したコンプライアンス関数を示す。式（４）は、体積制約を表し、特に、左辺は線形化した体積制約関数を示し、右辺は体積制約値ｆ_1cを示す。

【0050】

具体的には、式（２）～（４）から、Lagrangeの未定乗数λ₀、λ₁及びΔθを求め、このΔθを式（１）に代入することで、試行解θ_Tryが求められる。

【0051】

次に、図１５乃至図１７を参照して、本実施形態に係る制約変動法において行われる第１乃至第３処理について説明する。なお、第１乃至第３処理は、上述した手順により求められた試行解θ_Tryでの体積ｆ₁（θ_Try）が体積制約値ｆ_1c未満であるという体積制約が満たされていることを条件に行われるものとする。換言すると、体積ｆ₁（θ_Try）が体積制約値ｆ_1c以上である場合には、第１乃至第３処理は行われない。

【0052】

図１５は、本実施形態に係る制約変動法において行われる第１処理の説明図である。図１５は、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示している。第１処理では、処理装置３は、上述した手順により求められた試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が、元のθ分布値θ_kに対応するコンプライアンス値ｆ₀（θ_k）未満である場合に（ｆ₀（θ_Try）＜ｆ₀（θ_k））、コンプライアンス値が改善したものと判断する。そして、処理装置３は、この試行解θ_Tryを、次なる試行解θ_Tryを求めるために用いるべき新たなθ分布値として更新する。

【0053】

次いで、図１６は、本実施形態に係る制約変動法において行われる第２処理の説明図である。図１６は、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示すグラフを左右に２つ並べている。本実施形態では、処理装置３は、上述した手順により求められた試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が、θ分布値θ_kに対応するコンプライアンス値ｆ₀（θ_k）以上である場合に（ｆ₀（θ_Try）≧ｆ₀（θ_k））、コンプライアンス値が悪化したものと判断する。この場合、処理装置３は、上記した第１処理を行わずに、つまり当該試行解θ_Tryを次に用いるべき新たなθ分布値として更新せずに、以下のような第２処理を行う。すなわち、第２処理では、処理装置３は、仮目標値ｆ_0cをコンプライアンス値ｆ₀（θ_k）に近付けた値へと変更して、上記のように目的関数を線形近似した一次関数ＬＦ１（これは新たに計算することなく使い回される）が、この変更後の仮目標値ｆ_0cを満たす試行解θ_Try（以下では適宜「仮試行解θ_Try」と呼ぶ。）を新たに求める。具体的には、処理装置３は、変更前の元の仮目標値ｆ_0cを失敗値ｆ_failとして記憶する一方で、この失敗値ｆ_failとコンプライアンス値ｆ₀（θ_k）との中央値に、仮目標値ｆ_0cを更新する。よって、更新後の仮目標値ｆ_0cは「ｆ_0c＝０．５×（ｆ_fail＋ｆ₀（θ_k））」と表される。

【0054】

処理装置３は、更新された仮目標値ｆ_0cと一次関数ＬＦ１から求められた仮試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が、元のθ分布値θ_kに対応するコンプライアンス値ｆ₀（θ_k）未満になるまで、仮目標値ｆ_0cの更新及び仮試行解θ_Tryの算出を繰り返し行う（この場合、失敗値ｆ_failの更新も繰り返す）。そして、処理装置３は、仮試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）がコンプライアンス値ｆ₀（θ_k）未満になったときに、コンプライアンス値が改善したものと判断する。処理装置３は、このときの仮試行解θ_Tryを正式な試行解θ_Tryとして採用し、この試行解θ_Tryを、次なる試行解θ_Tryを求めるために用いるべき新たなθ分布値として更新する。

【0055】

このような第２処理によれば、的確な試行解θ_Tryが見つかるまでの繰り返しの処理において、元のθ分布値θ_kに対応するコンプライアンス値ｆ₀（θ_k）及び目的関数を線形近似した一次関数ＬＦ１を使い回せるので、計算時間に与える影響は少ない。つまり、第２処理によって大きな計算時間の増大を招くことはないと言える。

【0056】

次いで、図１７は、本実施形態に係る制約変動法において行われる第３処理の説明図である。この第３処理は、第１及び第２処理後に行われる処理である。図１７は、横軸にθ分布を示し、縦軸にコンプライアンスｆ₀を示すグラフを左右に２つ並べている。第３処理では、処理装置３は、上述した手順により求められた試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が仮目標値ｆ_0c未満である場合に（ｆ₀（θ_Try）＜ｆ_0c）、コンプライアンス値が過剰改善したものと判断する。そして、処理装置３は、仮目標値ｆ_0cをコンプライアンス値ｆ₀（θ_Try）よりも小さい値へと変更する。具体的には、処理装置３は、コンプライアンス値ｆ₀（θ_Try）と記憶している失敗値ｆ_fail（前回用いていた仮目標値ｆ_0c）との中央値に、仮目標値ｆ_0cを更新する。よって、更新後の仮目標値ｆ_0cは「ｆ_0c＝０．５×（ｆ_fail＋ｆ₀（θ_Try））」と表される。このような第３処理によれば、制約変動法の実行中に、順次求められる試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）に応じて、仮目標値ｆ_0cを的確な値へと更新することができる。

【0057】

なお、上記した第１乃至第３処理は、試行解θ_Tryでの体積ｆ₁（θ_Try）が体積制約値ｆ_1c未満である場合、つまり体積制約を満たす場合に行われものであるが、体積ｆ₁（θ_Try）が体積制約値ｆ_1c以上である場合、つまり体積制約を満たさない場合には、処理装置３は第４処理を行う。この第４処理では、処理装置３は、体積ｆ₁（θ_Try）が体積制約値ｆ_1c以上であるときに求められた試行解θ_Tryによって更新を行うが、体積ｆ₁（θ_Try）が体積制約値ｆ_1c未満になるまで、この試行解θ_Tryの更新を繰り返し行う、つまり試行解θ_Tryを繰り返し求める。この場合、処理装置３は、仮目標値ｆ_0cを更新しない。そして、処理装置３は、体積ｆ₁（θ_Try）が体積制約値ｆ_1c未満になると、このときの試行解θ_Tryを用いて、上記の第１又は第２処理を行う。

【0058】

次に、図１８及び図１９を参照して、本実施形態に係る位相最適化方法及び制約変動法の具体的な処理について説明する。図１８は、本実施形態に係る位相最適化方法の全体処理を示すフローチャートであり、図１９は、本実施形態に係る制約変動法の処理（制約変動法処理）を示すフローチャートである。なお、コンピュータ装置１０の処理装置３（詳しくはマイクロプロセッサ３ａ）が、メモリ３ｂに記憶されたプログラムを読み出して実行することにより、図１８及び図１９に示す処理が実現される。また、これらの処理は所定の周期で繰り返し実行される。

【0059】

位相最適化方法の全体処理が開始されると、ステップＳ１１において、処理装置３は、本処理に必要な種々のデータを取得する。具体的には、処理装置３は、構造解析用の数値計算モデルのデータ、つまり空間の離散化モデル及び境界条件（例えば有限要素法のメッシュ）に関するデータ、設計変数値としてのθ分布値の初期値｛θ₀｝、体積制約値ｆ_1c、仮目標値ｆ_0c、位相最適化方法の収束を判定するための収束判定値ｅｐｓ、を取得する。１つの例では、｛θ₀｝には、充填率ρの５０％に対応する値（つまり０）が適用され、体積制約値ｆ_1cには３０％が適用され、仮目標値ｆ_0cには０が適用される。

【0060】

次いで、ステップＳ１２において、処理装置３は、「［Ｋ（ρ_k（｛θ_k｝））］｛ｕ_k｝＝｛ｂ｝」という剛性方程式を解いて、変位｛ｕ_k｝を求める。ここで、｛θ_k｝は、ステップｋ回目のθ分布値（節点値（節点数の成分））であり、｛ｕ_k｝は、ステップｋ回目での剛性方程式を解いて得られた、部材２０の全体変位ベクトルの節点値（次元×節点数の成分）であり、「ｂ」は、部材２０に付与される全体外力ベクトルの節点値（次元×節点数の成分）であり、「［Ｋ（ρ_k（｛θ_k｝））］」は、ステップｋ回目での充填率ρ_k（｛θ_k｝）に依存する剛性行列（(次元×節点数）×（次元×節点数）の成分)である。

【0061】

次いで、ステップＳ１３において、処理装置３は、ステップＳ１２で求められた｛ｕ_k｝からコンプライアンス値ｆ₀を求めると共に、部材２０の体積ｆ₁を求める。具体的には、処理装置３は、「ｆ₀＝｛ｂ｝・｛ｕ_k｝」を解いてコンプライアンス値ｆ₀を求めると共に、「ρ_k（｛θ_k｝）」を部材２０の空間容量で積分することにより体積ｆ₁を求める。

【0062】

次いで、ステップＳ１４において、処理装置３は、コンプライアンス値ｆ₀を｛θ_k｝で微分することによりコンプライアンス値ｆ₀の勾配｛ｇ₀｝を求め（｛ｇ₀｝＝∂ｆ₀／∂｛θ_k｝）、また、体積ｆ₁を｛θ_k｝で微分することにより体積ｆ₁の勾配｛ｇ₁｝を求める（｛ｇ₁｝＝∂ｆ₁／∂｛θ_k｝）。

【0063】

次いで、ステップＳ１５において、処理装置３は、本実施形態に係る制約変動法処理を実行する。具体的には、処理装置３は、試行解θ_Tryの算出及びこの試行解θ_Tryに応じた第１乃至第４処理を行う。この詳細は、図１９を参照して後述する。

【0064】

次いで、ステップＳ１６において、処理装置３は、解（θ_k）が収束したか否かを判定する。具体的には、現在（ステップｋ回目）のコンプライアンス値ｆ₀（θ_k）と仮目標値ｆ_0cとの差の絶対値が収束判定値ｅｐｓ未満であるか否かを判定すると共に（｜ｆ₀（θ_k）－ｆ_0c｜＜ｅｐｓ）、現在（ステップｋ回目）のコンプライアンス値ｆ₀（θ_k）と失敗値ｆ_failとの差の絶対値が収束判定値ｅｐｓ未満であるか否かを判定する（｜ｆ₀（θ_k）－ｆ_fail｜＜ｅｐｓ）。処理装置３は、これら両方の条件が成立しない場合には、解が収束していないと判定して（ステップＳ１６：Ｎｏ）、ステップＳ１２に戻る。この場合には、処理装置３は、解が収束したと判定されるまで、ステップＳ１２～Ｓ１６の処理を繰り返す。

【0065】

これに対して、処理装置３は、「｜ｆ₀（θ_k）－ｆ_0c｜＜ｅｐｓ」及び「｜ｆ₀（θ_k）－ｆ_fail｜＜ｅｐｓ」の両方の条件が成立する場合には、解が収束したと判定して（ステップＳ１６：Ｙｅｓ）、位相最適化方法の全体処理を終了する。この場合、処理装置３は、現在の解（θ_k）を、θ分布値の最適解（つまりコンプライアンス値ｆ₀が最小となるときのθ分布値）として決定する。そして、処理装置３は、この最適解に対応する、最適な充填率ρ分布（部材２０の最適形状）を求める。例えば、処理装置３は、こうして求められた最適な充填率ρ分布に対応する画像やそれに関連する情報を、出力装置５の表示装置に表示させる。

【0066】

次に、図１９を参照して、本実施形態に係る制約変動法処理について説明する。まず、ステップＳ２１において、処理装置３は、上記の式（１）～（４）の連立方程式を解いて試行解θ_Tryを求める。具体的には、処理装置３は、式（２）～（４）から、Lagrangeの未定乗数λ₀、λ₁及びΔθを求め、このΔθを式（１）に代入することで、試行解θ_Tryを求める。

【0067】

次いで、ステップＳ２２において、処理装置３は、ステップＳ２１で得られた試行解θ_Tryから体積ｆ₁（θ_Try）を求め、この体積ｆ₁（θ_Try）が体積制約値ｆ_1c未満であるか否か（ｆ₁（θ_Try）＜ｆ_1c）、つまり体積制約が満たされているか否かを判定する。処理装置３は、体積制約が満たされている場合（ステップＳ２２：Ｙｅｓ）、ステップＳ２４に進み、体積制約が満たされていない場合（ステップＳ２２：Ｎｏ）、ステップＳ２３に進む。

【0068】

ステップＳ２３において、処理装置３は、体積制約が満たされていないので、第４処理を行う。具体的には、処理装置３は、仮目標値ｆ_0cは更新せずに、現在の試行解θ_Tryを次にステップＳ２１を行うべきθ分布値として更新する。この後、処理装置３は、ステップＳ２１に戻り、更新されたθ分布値を用いて試行解θ_Tryを再度求め、そして、ステップＳ２２に進み、この試行解θ_Tryに対応する体積ｆ₁（θ_Try）に基づき、体積制約が満たされているか否かを判定する。すなわち、処理装置３は、ステップＳ２１～Ｓ２３において、体積制約が満たされるまで、試行解θ_Tryを繰り返し求める。

【0069】

他方で、ステップＳ２４において、処理装置３は、ステップＳ２１で得られた試行解θ_Tryからコンプライアンス値ｆ₀（θ_Try）を求め、このコンプライアンス値ｆ₀（θ_Try）が元のコンプライアンス値ｆ₀（θ_k）未満であるか否か（ｆ₀（θ_Try）＜ｆ₀（θ_k））、つまり試行解θ_Tryが改善したか否かを判定する。処理装置３は、試行解θ_Tryが改善した場合（ステップＳ２４：Ｙｅｓ）、ステップＳ２５に進む。この場合、処理装置３は、第１処理を行う（ステップＳ２５）。具体的には、処理装置３は、現在の試行解θ_Tryを、次なる試行解θ_Tryを求めるために用いるべき新たなθ分布値として更新する。そして、処理装置３は、ステップＳ２７に進む。

【0070】

これに対して、処理装置３は、試行解θ_Tryが改善していない場合（ステップＳ２４：Ｎｏ）、つまり試行解θ_Tryが悪化した場合、ステップＳ２６に進む。この場合、処理装置３は、第２処理を行う（ステップＳ２６）。具体的には、処理装置３は、現在の仮目標値ｆ_0cを失敗値ｆ_failとして記憶する一方で、この失敗値ｆ_failと元のコンプライアンス値ｆ₀（θ_k）との中央値に仮目標値ｆ_0cを更新して（ｆ_0c＝０．５×（ｆ_fail＋ｆ₀（θ_k）））、この更新後の仮目標値ｆ_0cに基づき試行解θ_Try（仮試行解θ_Try）を新たに求める。処理装置３は、この仮試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が、元のコンプライアンス値ｆ₀（θ_k）未満になるまで、仮目標値ｆ_0cの更新及び仮試行解θ_Tryの算出を繰り返し行う（この場合、失敗値ｆ_failの更新も繰り返す）。そして、処理装置３は、仮試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）がコンプライアンス値ｆ₀（θ_k）未満になったときに、コンプライアンス値が改善したものと判断する。処理装置３は、このときの仮試行解θ_Tryを正式な試行解θ_Tryとして採用し、この試行解θ_Tryを、次なる試行解θ_Tryを求めるために用いるべき新たなθ分布値として更新する。そして、処理装置３は、ステップＳ２７に進む。

【0071】

次いで、ステップＳ２７において、処理装置３は、試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）が仮目標値ｆ_0c未満であるか否か（ｆ₀（θ_Try）＜ｆ_0c）、つまり試行解θ_Tryが過剰改善したか否かを判定する。処理装置３は、試行解θ_Tryが過剰改善した場合（ステップＳ２７：Ｙｅｓ）、ステップＳ２８に進み、試行解θ_Tryが過剰改善していない場合（ステップＳ２７：Ｎｏ）、制約変動法処理を終了する。

【0072】

ステップＳ２８において、処理装置３は、試行解θ_Tryが過剰改善したので、第３処理を実行する。具体的には、処理装置３は、試行解θ_Tryに対応するコンプライアンス値ｆ₀（θ_Try）と、記憶している失敗値ｆ_fail（前回用いていた仮目標値ｆ_0c）との中央値に、仮目標値ｆ_0cを更新する（ｆ_0c＝０．５×（ｆ_fail＋ｆ₀（θ_Try）））。そして、処理装置３は、制約変動法処理を終了する。

【0073】

なお、図１９において、ステップＳ２１は、本発明における第１ステップに相当し、ステップＳ２４、Ｓ２５は、本発明における第２ステップに相当し、ステップＳ２４、Ｓ２６は、本発明における第３及び第４ステップに相当し、ステップＳ２７、Ｓ２８は、本発明における第６ステップに相当し、図１８のステップＳ１６は、本発明における第５ステップに相当する。

【0074】

＜作用及び効果＞
次に、図２０及び図２１を参照して、本実施形態に係る最適解の演算方法による作用及び効果について具体的に説明する。

【0075】

図２０は、本実施形態に係る制約変動法による結果の一例を示す。図２０（Ａ）は、勾配法、ニュートン法及び制約変動法のそれぞれを位相最適化方法に適用した場合に得られた最適な充填率分布を示す。ここでは、部材２０の左辺が固定され、部材２０の右辺の中央部に荷重が付与される場合を例示している。図２０（Ｂ）は、勾配法、ニュートン法及び制約変動法のそれぞれによる、最適な充填率分布（つまり最小の（最適な）コンプライアンス）が得られるまでの計算ステップ（計算時間に対応する）を示している。具体的には、グラフＧ１２は、図８に示したものと同じ勾配法による結果を示し、グラフＧ３は、図１１に示したものと同じニュートン法による結果を示し、グラフＧ４は、本実施形態に係る制約変動法による結果を示している。図２０（Ｂ）から、本実施形態に係る制約変動法によれば、勾配法及びニュートン法と比べて、最適解（最小のコンプライアンス）が求まるまでに要する計算ステップ数が非常に少ないこと、及び、ニュートン法と比べて、計算が非常に安定していることがわかる。

【0076】

図２１は、本実施形態に係る制約変動法による結果の他の例を示す。図２１（Ａ）は、勾配法、ニュートン法及び制約変動法のそれぞれを位相最適化方法に適用した場合に得られた最適な充填率分布を示す。ここでは、部材２０の底辺の２箇所が固定され、部材２０内の５箇所に荷重が付与される場合を例示している。図２１（Ｂ）は、勾配法、ニュートン法及び制約変動法のそれぞれによる、最適な充填率分布（つまり最小の（最適な）コンプライアンス）が得られるまでの計算ステップ（計算時間に対応する）を示している。具体的には、グラフＧ５１は、勾配法による結果を示し、グラフＧ５２は、ニュートン法による結果を示し、グラフＧ５３は、本実施形態に係る制約変動法による結果を示している。図２１（Ｂ）からも、本実施形態に係る制約変動法によれば、勾配法及びニュートン法と比べて、最適解（最小のコンプライアンス）が求まるまでに要する計算ステップ数が非常に少ないこと、及び、計算が非常に安定していることがわかる。

【0077】

このようなことから、本実施形態に係る制約変動法によれば、位相最適化方法において最適解を安定且つ高速に求めることができると言える。なお、本実施形態によれば、図２０及び図２１に示した例以外に、部材２０における種々の多点曲げ問題やＬ字曲げ問題やオフセット曲げ問題やなどに適用した場合にも、最適解を安定且つ高速に求めることができる（例えば１５ステップ前後で解を収束させることができる）。

【0078】

以上まとめると、本実施形態に係る最適解の演算方法は、
（１）或るθ分布値と、当該θ分布値に対応するコンプライアンス値に基づき目的関数を線形近似した関数が、仮目標値を満たす試行解を求める第１ステップと、
（２）第１ステップで求められた試行解に対応するコンプライアンス値が、第１ステップで用いられたθ分布値に対応するコンプライアンス値未満である場合に、この試行解を、次に第１ステップを行うべき新たなθ分布値として設定する第２ステップと、
（３）第１ステップで求められた試行解に対応するコンプライアンス値が、第１ステップで用いられたθ分布値に対応するコンプライアンス値以上である場合に、仮目標値を第１ステップで用いられたコンプライアンス値に近付けた値へと変更して、上記の線形近似した関数が、変更された仮目標値を満たす試行解を求める第３ステップであって、変更された仮目標値から求められた試行解に対応するコンプライアンス値が、第１ステップで用いられたコンプライアンス値未満になるまで、仮目標値を第１ステップで用いられたコンプライアンス値に徐々に近付ける変更を繰り返し行うと共に、線形近似した関数が変更された仮目標値を満たす試行解を繰り返し求める第３ステップと、
（４）第３ステップで求められた試行解に対応するコンプライアンス値が、第１ステップで用いられたコンプライアンス値未満になったときに、この試行解を、次に第１ステップを行うべき新たなθ分布値として設定する第４ステップと、
（５）第２又は第４ステップで設定されたθ分布値に対応するコンプライアンス値と仮目標値との差分が所定値未満になるまで、第１乃至第４ステップを繰り返し行わせ、差分が所定値未満になったときに、このときに用いられていたθ分布値を最適解として決定する第５ステップと、
を有する。

【0079】

このような本実施形態によれば、悪い解も受け入れるしかなかった従来の手法（特にニュートン法）と異なり、制約値を変動させることにより、悪い解を切り捨てて良い解を確実に求めることができるので、安定した計算を行うことが可能となる。また、本実施形態に係る最適解の演算方法によれば、良い解を見つけるまでの繰り返しの処理において、目的関数を線形近似した一次関数及び元のθ分布値に対応するコンプライアンス値を使い回すことができるので、計算ステップ数の増大を効果的に抑制することができる。したがって、本実施形態に係る最適解の演算方法によれば、最適解を安定且つ高速に求めることができるのである。

【0080】

また、本実施形態に係る最適解の演算方法によれば、第３ステップでは、第１ステップで用いられたθ分布値に対応するコンプライアンス値と仮目標値との中央値を求め、この中央値によって仮目標値を変更する。これにより、元のコンプライアンス値及び元の仮目標値に応じた的確な値へと仮目標値を更新することができる。

【0081】

また、本実施形態に係る最適解の演算方法によれば、第１又は第３ステップで求められた試行解に対応するコンプライアンス値が仮目標値未満である場合、当該仮目標値を当該コンプライアンス値よりも小さい値へと変更する第６ステップを有する。これにより、制約変動法の実行中に、順次求められる試行解に対応するコンプライアンス値に応じて、仮目標値を的確な値へと更新することができる。

【0082】

また、本実施形態に係る最適解の演算方法によれば、第６ステップでは、第１又は第３ステップで求められた試行解に対応するコンプライアンス値が仮目標値未満である場合、当該コンプライアンス値と当該仮目標値の変更前に用いられていた仮目標値（失敗値）との間の中央値を求め、この中央値によって仮目標値を変更する。これにより、試行解に対応するコンプライアンス値及び前回用いていた仮目標値（失敗値）に応じた的確な値へと仮目標値を更新することができる。

【0083】

また、本実施形態に係る最適解の演算方法によれば、第５ステップでは、第２又は第４ステップで設定されたθ分布値に対応するコンプライアンス値と仮目標値との差分が所定値未満で、且つ、第２又は第４ステップで設定されたθ分布値に対応するコンプライアンス値と、仮目標値の変更前に用いられていた仮目標値（失敗値）との差分が所定値未満となったときに、このときに用いられていたθ分布値を最適解として決定する。これにより、解の収束判定を精度良く行うことができる、つまりより正確な最適解を求めることができる。

【0084】

また、本実施形態に係る最適解の演算方法によれば、第１ステップで求められた試行解から部材２０の体積を求め、この体積が所定の制約値以上である場合には、体積が制約値未満となるまで、第２及び第３ステップを行わずに、第１ステップを繰り返し行う。これにより、体積制約を確実に満たす解を求めることができる。

【0085】

＜本発明の適用例＞
上記した実施形態では、本発明に係る最適解の演算方法を位相最適化方法に適当する例を示したが、本発明に係る最適解の演算方法の適用はこれに限定されない。ここでは、本発明に係る最適解の演算方法の適用例について説明する。

【0086】

本発明に係る最適解の演算方法は、例えば種々の自然現象などを数値計算で求める場合に適用できる。この場合、本発明に係る最適解の演算方法は、種々の自然現象を表すシステム方程式（連立一次方程式）を解くために利用される。１つの例では、本発明に係る最適解の演算方法は、構造力学分野における剛性方程式「［Ｋ］｛ｕ｝＝｛ｂ｝」というシステム方程式（［Ｋ］：剛性行列、｛ｕ｝：変位ベクトル、｛ｂ｝：荷重ベクトル）や、伝熱分野における熱伝導方程式「［Ａ］｛Ｔ｝＝｛Ｆ｝」というシステム方程式（［Ａ］：熱伝導行列、｛Ｔ｝：温度、｛Ｆ｝：熱流束ベクトル）や、流体力学におけるNavier-Stokes方程式「［Ｂ］｛ｐ｝＝｛ｖ｝」というシステム方程式（［Ｂ］：ポアソン行列、｛ｐ｝：圧力、｛ｖ｝：中間流速ベクトル）などに適用できる。

【0087】

ここで、バネ定数Ｋのバネに力ｂが働くことで変位ｕが生じる場合の剛性方程式「Ｋｕ＝ｂ」を例（１自由度の例）に挙げる。この場合、「Ｆ（ｕ）＝Ｋｕ²／２－ｂｕ」という関数を定義する。バネ定数Ｋは正の値なので、「Ｆ（ｕ）」は下に凸な関数である。この「Ｆ（ｕ）」を「ｕ」で微分すると、「∂Ｆ（ｕ）／∂ｕ＝Ｋｕ－ｂ」となる。「∂Ｆ（ｕ）／∂ｕ＝０」のときに「Ｆ（ｕ）」の傾きが０になる。よって、傾きが０になるときは「Ｋｕ－ｂ＝０」となり、これを解くと「Ｋｕ＝ｂ」という解が得られる。したがって、剛性方程式を解くことは、「Ｆ（ｕ）」の最小値を求めることに相当する。

【0088】

次いで、自由度が複数ある剛性方程式「［Ｋ］｛ｕ｝＝｛ｂ｝」を例に挙げる（｛ｕ｝は各節点の変位）。この場合、「Ｆ（｛ｕ｝）＝（［Ｋ］｛ｕ｝・｛ｕ｝）／２－｛ｂ｝・｛ｕ｝」というスカラー関数を定義する。［Ｋ］が正定値行列ならば、Ｆ（｛ｕ｝）は凸関数になる（なお、自然界のシステム方程式の多くは正定値行列により定義される）。ここで、「Ｆ（｛ｕ｝）」を「｛ｕ｝」で微分すると、「∂Ｆ（｛ｕ｝）／∂｛ｕ｝＝［Ｋ］｛ｕ｝－｛ｂ｝」となる。したがって、「∂Ｆ（｛ｕ｝）／∂｛ｕ｝＝０」を解くこと、つまり「Ｆ（｛ｕ｝）」の最小値を求めることは、「［Ｋ］｛ｕ｝＝｛ｂ｝」を解くことと同じである。

【0089】

このようから、従来の手法では、種々の自然現象を表す膨大な数の連立方程式（システム方程式）を解く必要があったが、本発明では、１つのスカラー方程式（上述した「Ｆ（｛ｕ｝）＝（［Ｋ］｛ｕ｝・｛ｕ｝）／２－｛ｂ｝・｛ｕ｝」など）を解けばよい、具体的にはスカラー凸関数の最小値を求めればよい。上述したように、本発明によれば、凸関数である目的関数の最適解を安定且つ高速に求めることができる。したがって、本発明によれば、種々の自然現象を表すシステム方程式を安定且つ高速に解くことができるのである。

【0090】

なお、本発明は、自然現象を数値計算で求める場合に適用することに限定はされず、他の種々の数値計算に適用可能である。要は、本発明は、内部変数により規定され且つ凸関数である目的関数の最適解を演算するために適用される。また、上記した実施形態では、下に凸な目的関数に本発明を適用する例を示したが、本発明は上に凸な目的関数にも適用できる。この場合、目的関数の最大値に対応する最適解を求めればよい。

【符号の説明】

【0091】

１ 入力装置
３ 処理装置
３ａ マイクロプロセッサ
３ｂ メモリ
５ 出力装置
１０ コンピュータ装置
２０ 部材

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】