(57)【特許請求の範囲】
【請求項１】
プログラムを記憶した記憶手段と、プログラムを実行する実行手段とを備え、前記実行手段が前記記憶手段に記憶されたプログラムを実行することにより実現される、第１及び第２の求解手段と、勾配ベクトル計算手段と、収束判定手段と、第１及び第２の係数計算手段と、探索ベクトル計算手段と、設計変数ベクトル更新手段とを有する構造最適設計装置において、
前記第１の求解手段が、構造部材の存在可能な領域を分割した各領域における構造部材の存在率を表す設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第２の求解手段が、前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計方法において、
前記第２の求解工程は、
前記勾配ベクトル計算手段が、前記設計変数ベクトルに関する前記第２の評価関数の勾配ベクトルを計算する勾配ベクトル計算工程と、
前記収束判定手段が、前記勾配ベクトルのノルムの値が所定値未満か否かを判定する収束判定工程と、
前記第１の係数計算手段が、前記勾配ベクトルのノルムの値に基づいて第１の係数を計算する第１の係数計算工程と、
前記探索ベクトル計算手段が、該第１の係数と前記勾配ベクトルとに基づいて探索ベクトルを計算する探索ベクトル計算工程と、
前記第２の係数計算手段が、前記設計変数ベクトルと前記探索ベクトルとに基づいて第２の係数を計算する第２の係数計算工程と、
前記設計変数ベクトル更新手段が、前記第２の係数と前記探索ベクトルとに基づいて前記設計変数ベクトルを更新する設計変数ベクトル更新工程とを備え、
前記第２の求解工程を、前記収束判定工程において前記第２の評価関数の勾配ベクトルのノルムの値が前記所定値未満となるまで繰り返し実行し、
前記第２の係数計算工程は前記第１の求解工程を含み、該第１の求解工程は勾配ベクトルに基づく反復解法であることを特徴とする最適設計方法。
【請求項２】
前記構造最適設計装置は、前記実行手段が前記記憶手段に記憶されたプログラムを実行することにより実現される勾配ベクトル修正手段を更に有し、前記第２の求解工程は、前記勾配ベクトル修正手段が、前記勾配ベクトル計算工程により計算された勾配ベクトルを前記設計変数ベクトルに対する拘束条件に基づいて修正する勾配ベクトル修正工程を更に備えることを特徴とする請求項１に記載の最適設計方法。
【請求項３】
前記構造最適設計装置は、前記実行手段が前記記憶手段に記憶されたプログラムを実行することにより実現される探索ベクトル修正手段を更に有し、前記第２の求解工程は、前記探索ベクトル修正手段が、前記探索ベクトル計算工程により計算された探索ベクトルを前記設計変数ベクトルに対する拘束条件に基づいて修正する探索ベクトル修正工程を更に備えることを特徴とする請求項１に記載の最適設計方法。
【請求項４】
前記構造最適設計装置は、前記実行手段が前記記憶手段に記憶されたプログラムを実行することにより実現される設計変数ベクトル修正手段を更に有し、前記第２の求解工程は、前記設計変数ベクトル修正手段が、前記設計変数ベクトル更新工程により更新された設計変数ベクトルを前記設計変数ベクトルに対する拘束条件に基づいて修正する設計変数ベクトル修正工程を更に備えることを特徴とする請求項１に記載の最適設計方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、設計パラメタの最適化問題の解法に関し、特に構造部材のトポロジーを含む形状最適化のための自動設計技術に関するものである。
【０００２】
【従来の技術】
構造トポロジー最適設計とは、与えられた条件の下で最適な、構造部材のトポロジーおよび形状寸法を決定する問題である。以下、構造部材のトポロジーおよび形状寸法を設計変関数といい、上記決定問題を設計変関数最適化問題という。変関数という理由は、トポロジーおよび形状寸法が３次元空間の関数になっているからである。
【０００３】
設計変関数最適化問題では、各設計変関数の値に対して、状態変関数の最適化問題を解かなければならない。この意味から、構造トポロジー最適設計は、内側に状態変関数の最適化問題、外側に設計変関数の最適化問題を持つ、２重構造最適化問題と捉える事が出来る。内側の状態変関数最適化問題では、従来技術の蓄積により、空間を有限個の要素に分割するという考え方が採用される。
【０００４】
特に構造部材の歪エネルギーを評価汎関数としている問題では、その解析手法として有限要素法が一般的である。有限要素法の解法としては一次方程式に対する直接法が採用されている。一方、設計変関数の最適化問題に対しては、大別すると、以下に示す３種類の方法が提供されている（〔文献１〕或いは〔文献２〕）：
１．Ｅｖｏｌｕｔｉｏｎａｒｙ法（以下、Ｅ法）
２．Ｈｏｍｏｇｅｎｉｚａｔｉｏｎ法（以下、Ｈ法）
３．Ｍａｔｅｒｉａｌ  ｄｉｓｔｒｉｂｕｔｉｏｎ法（以下、ＭＤ法）
Ｅ法では、空間を分割することによって得られる部分空間のそれぞれをセルと称し、その生成と消滅を適当な規則によって繰り返す。構造部材は、最終的に存在しているセルの集合として与えられる。セルが存在するか否かという２つの状態のみを許すことにより、明確な構造部材が得られる。また、評価汎関数の微分情報を用いないので、局所最適解にトラップされないことから、評価汎関数が多峰性の場合に有効である。
【０００５】
〔文献３〕では、Ｅ法の一種である遺伝的探索法を用いた骨組構造部材の最適化設計装置が提供されている。この最適化設計装置では、骨組部材断面寸法などの離散設計変数データの近似式を使用する近似最適化計算装置と、該設計変数データを使用する詳細最適化計算装置を設け、これら２つの計算装置を結合して、膨大な設計変数が存在する実設計問題に対応することが出来るようにしている。
【０００６】
Ｈ法は、分割された各部分空間に位置する構造要素に、更に微細な構造を仮定し、連続値を取る設計変関数を新たに導入することによって、感度解析の採用を可能にしている。ここで感度解析とは、設計変関数に関する評価汎関数の微分情報を利用した解析手法のことであり、感度解析が可能になれば、勾配法のような反復解法を用いることが出来、Ｅ法のような総当り的手法に比べて、少なくとも局所最適解の探索に掛かる計算時間は大幅に短縮される。
【０００７】
ＭＤ法は、構造部材の存在確率を示す０から１の範囲の実数を各要素に割当てることによって、構造部材のトポロジーや形状寸法変化を表現する方法である。構造部材が存在するか否かという離散的な情報を存在確率という連続値に置き換えることによって感度解析を可能にしたという意味でＨ法と同様のものであるが、Ｈ法よりパラメタ数が少ない分、モデル化が容易であり適用範囲も広い。
【０００８】
〔文献５〕にはＭＤ法による構造物の位相最適化手法が開示されている。本手法の特徴は以下のとおりである：
（１）ボクセル有限要素法（空間を等間隔に分割）を用いているため、あらゆる要素に対する要素剛性マトリクスが同一である。従って、要素剛性マトリクスを予め１度だけ計算しておけば、以後の計算に利用できる。更に、要素が規則正しく配置されているため、各要素の節点番号情報を記憶する必要がない。
【０００９】
（２）大規模連立一次方程式を解くために、前処理付共役勾配法とＥｌｅｍｅｎｔ−ｂｙ−Ｅｌｅｍｅｎｔ法を組み合わせて用いたことにより、全体剛性マトリクスを組み立てることなく解が求められるので、処理に必要なメモリ容量が少なくて済む。
【００１０】
（３）均質化法では、１要素に対して６個の設計変数（３次元の場合）が必要になる。更に設計変数が変化するたびに要素剛性マトリクスを再計算しなければならない。一方、構造部材の存在率を密度比で表現する密度法を採用することによって、１要素に対して１つの設計変数でよい。また設計変数の変化が要素剛性マトリクスに影響を与えない。
〔文献１〕Ｓ．Ｂｕｌｍａｎ，Ｊ．Ｓｉｅｎｚ，Ｅ．Ｈｉｎｔｏｎ：“Ｃｏｍｐａｒｉｓｏｎｓ  ｂｅｔｗｅｅｎ  ａｌｇｏｒｉｔｈｍｓ  ｆｏｒ  ｓｔｒｕｃｔｕｒａｌ  ｔｏｐｏｌｏｇｙ  ｏｐｔｉｍｉｚａｔｉｏｎ  ｕｓｉｎｇ  ａ  ｓｅｒｉｅｓ  ｏｆ  ｂｅｎｃｈｍａｒｋ  ｓｔｕｄｉｅｓ”，  Ｃｏｍｐｕｔｅｒｓ  ａｎｄ  Ｓｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１２０３−１２１８（２００１）．
〔文献２〕Ｙ−Ｌ．Ｈｓｕ，Ｍ−Ｓ．Ｈｓｕ，Ｃ−Ｔ．Ｃｈｅｎ：“Ｉｎｔｅｒｐｒｅｔｉｎｇ  ｒｅｓｕｌｔｓ  ｆｒｏｍ  ｔｏｐｏｌｏｇｙ  ｏｐｔｉｍｉｚａｔｉｏｎ  ｕｓｉｎｇ  ｄｅｎｓｉｔｙ  ｃｏｎｔｏｕｒｓ”，  Ｃｏｍｐｕｔｅｒｓ  ａｎｄ  Ｓｔｒｕｃｔｕｒｅｓ，７９，ｐｐ．１０４９−１０５８（２００１）．
〔文献３〕特開２００１−１３４６２８号公報
〔文献４〕山川宏：“最適化デザイン”，計算力学とＣＡＥシリーズ９，培風館（１９９６）
〔文献５〕藤井，鈴木，大坪：“ボクセル有限要素法を用いた構造物の位相最適化”，Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｆ  ＪＳＣＥＳ，  Ｐａｐｅｒ  Ｎｏ．２０００００１０（２０００）．
【００１１】
【発明が解決しようとする課題】
しかしながら、従来技術には以下に述べるような課題があった。
【００１２】
一般に構造最適化問題は、設計変数ベクトルの最適化問題の各反復ステップ毎に、状態変数ベクトルの最適化問題を含む、２重の最適化問題として定式化される。設計変数ベクトルの最適化問題を外側の最適化、状態変数ベクトルの最適化問題を内側の最適化と呼ぶと、内側の最適化は、設計変数ベクトルをパラメタとし、即ち設計変数ベクトルを固定し、状態変数ベクトルを求める問題である。これは通常、構造解析と呼ばれ、有限要素法によって線形方程式の解法を用いて解く事が出来る。
【００１３】
ところが、構造が変化し、ある領域の構造部材が存在しなくなると、例えば部材に穴が空いてしまうと、対応する要素の設計変数が０となり、結果的に要素のヤング率が０となる。すると、上記線形方程式の係数行列がフルランクでなくなり、逆行列が存在しないという理由で、直接法では解けなくなる。
【００１４】
そのため、殆どの従来技術では、設計変数ベクトルが０となる場合には、正確に０とはせずに０に近い微小な値で置き換えるという方法を取ってきた。〔文献５〕で採用されているボクセル有限要素法でも同様である。
【００１５】
しかるに、要素の設計変数を微小な値に設定すると言うことは、物理的には薄い膜、或いは弱い部材が存在することに相当する。つまり従来法では、物質がない部分を正確に表現していないことになる。
【００１６】
【課題を解決するための手段】
上記課題を解決するために、本発明によれば、プログラムを記憶した記憶手段と、プログラムを実行する実行手段とを備え、前記実行手段が前記記憶手段に記憶されたプログラムを実行することにより実現される、第１及び第２の求解手段と、勾配ベクトル計算手段と、収束判定手段と、第１及び第２の係数計算手段と、探索ベクトル計算手段と、設計変数ベクトル更新手段とを有する構造最適設計装置において、前記第１の求解手段が、構造部材の存在可能な領域を分割した各領域における構造部材の存在率を表す設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、前記第２の求解手段が、前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計方法において、前記第２の求解工程は、前記勾配ベクトル計算手段が、前記設計変数ベクトルに関する前記第２の評価関数の勾配ベクトルを計算する勾配ベクトル計算工程と、前記収束判定手段が、前記勾配ベクトルのノルムの値が所定値未満か否かを判定する収束判定工程と、前記第１の係数計算手段が、前記勾配ベクトルのノルムの値に基づいて第１の係数を計算する第１の係数計算工程と、前記探索ベクトル計算手段が、該第１の係数と前記勾配ベクトルとに基づいて探索ベクトルを計算する探索ベクトル計算工程と、前記第２の係数計算手段が、前記設計変数ベクトルと前記探索ベクトルとに基づいて第２の係数を計算する第２の係数計算工程と、前記設計変数ベクトル更新手段が、前記第２の係数と前記探索ベクトルとに基づいて前記設計変数ベクトルを更新する設計変数ベクトル更新工程とを備え、前記第２の求解工程を、前記収束判定工程において前記第２の評価関数の勾配ベクトルのノルムの値が前記所定値未満となるまで繰り返し実行し、前記第２の係数計算工程は前記第１の求解工程を含み、該第１の求解工程は勾配ベクトルに基づく反復解法であることを特徴とする。
【００１７】
【発明の実施の形態】
本実施形態の構成を図２に示す。図中、２０１はＣＰＵ、２０２は表示装置、２０３は入力装置、２０４は１次記憶装置、２０５は２次記憶装置、２０６は通信装置、２０７はバスラインである。
【００１８】
本実施形態における各処理機能は、プログラムとして予め２次記憶装置２０５に格納されており、入力装置２０３或いは通信装置２０６等からのコマンド入力によって、１次記憶装置にロードされ、実行されるものである。
以下の説明のために対象とする問題の定式化を行う。
【００１９】
有限要素法による定式化により、変関数は有限次元ベクトルで表現されるものとすると、変関数の評価汎関数は、変数ベクトルの評価関数となる。以下、有限次元ベクトルで表現されているとして記述する。
状態変数ベクトルｘ、設計変数ベクトルｆをそれぞれ列ベクトルとして以下のように書く：
（式１）ｘ＝（ｘ１，ｘ２，…，ｘｍ）Ｔ
（式２）ｆ＝（ｆ１，ｆ２，…，ｆｎ）Ｔ
ここでＴは転置を表す。ｘはｍ次元ベクトル、ｆはｎ次元ベクトルである。
ｘおよびｆの境界条件をそれぞれＢ１、Ｂ２とする：
（式３）Ｂ１（ｘ）＝０
（式４）Ｂ２（ｆ）＝０
状態変数ベクトルｘおよび設計変数ベクトルｆに関する評価関数を、それぞれＬ１、Ｌ２とする：
（式５）Ｌ１＝Ｌ１（ｘ；ｆ）
（式６）Ｌ２＝Ｌ２（ｆ，ｘ）
Ｌ１はｘを変数ベクトル、ｆをパラメタとする関数であり、Ｌ２はｘおよびｆを変数ベクトルとする関数である。
【００２０】
一般の最適化問題では、Ｋ１個の等式拘束条件とＫ２個の不等式拘束条件とが課せられることが多い。そこで、設計変数に関するｊ番目の等式拘束条件およびｋ番目の不等式拘束条件を、それぞれＱｊ、Ｒｋとする：
（式７）Ｑｊ（ｆ，ｘ）＝０
（式８）Ｒｋ（ｆ，ｘ）≧０
以上の表記より、最適設計問題は、拘束条件（式４）、（式５）、（式６）、（式７）、および（式８）の下で、次式を満足する解として与えられる：
（式９）ｍｉｎ［Ｌ２］
ただし、状態変数ｘは拘束条件（式３）の下で次式を満足する解として得られる：
（式１０）ｍｉｎ［Ｌ１］
上記定式化に基づき、本実施形態の処理を説明する。
以下では、説明を簡単にするために、（式７）の等式拘束条件を次式の１個とする：
【外１】
0004101047.tif 000002
【００２１】
ここでｃｊ（ｊ＝０，１，・・・，ｎ）は定数である。
【００２２】
また（式８）の不等式拘束条件を次式のように、値の範囲に対するものとする：
（式１２）Ｒｊ−（ｆ）＝ｆｔ（ｊ）−ｆｊｍｉｎ≧０
Ｒｊ＋（ｆ）＝ｆｊｍａｘ−ｆｔ（ｊ）≧０
ただし、ｊ＝１，２，…，ｎ。
【００２３】
従って、不等式拘束条件の数は設計変数の数の２倍となる。
【００２４】
設計変数ベクトルと同次元のフラグベクトルａ＋、およびａ−を以下のように定義する：
ａ＋＝０，      ｉｆ  ０≧Ｒｊ＋（ｆ）
ａ＋＝１，      ｏｔｈｅｒｗｉｓｅ
および
ａ−＝０，      ｉｆ  ０≧Ｒｊ−（ｆ）
ａ−＝１，      ｏｔｈｅｒｗｉｓｅ
また、ａ＋およびａ−が共に１であるｊの集合を∧１、それ以外のインデクスの集合を∧２とする。
【００２５】
図１には、第二の求解工程のフローチャートが示されている。図中、ステップＳ１０１は、シミュレーション対象となる系の諸元を読み込む処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値、境界条件Ｂ１、Ｂ２、評価関数Ｌ１、Ｌ２、拘束条件Ｑｊ、Ｒｋが含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保し、値を設定する。
【００２６】
ステップＳ１０２では、反復回数ｔを１に初期化する。反復回数ｔは、上述の第２の求解工程における第１の反復処理を行なった回数であり、以下のステップＳ１０３から１１１の処理の反復回数が設定値Ｔに達したら終了する。
【００２７】
ステップＳ１０３では勾配ベクトル計算工程を実行する。同工程では、設計変数ベクトルに関する第２の評価関数のグラジェントベクトルｇｔを計算する。ｇｔは次式で計算される：
（式１３）ｇｔ＝（∂Ｌ２／∂ｆ１，∂Ｌ２／∂ｆ２，…，∂Ｌ２／∂ｆｎ）Ｔ
∂Ｌ２／∂ｆｊは解析的に解ける場合とそうでない場合がある。解析解がある場合は、その計算工程を利用して（式１０）を得る。そうでない場合は、自動微分という技術を利用して計算することが出来る。自動微分技術は、以下の〔文献６〕等によって公知である。
〔文献６〕久保田光一，伊理正夫：“アルゴリズムの自動微分と応用，”現代非線形科学シリーズ，コロナ社（１９９８年）
ステップＳ１０４では、（式１３）で得られたｇｔに対して、第一のベクトル修正工程を実行する。同工程は図６を用いて後述する。
【００２８】
ステップＳ１０５では、収束判定工程を実行する。同工程では、ｇｔのノルムの２乗を計算し、その値が０或いは０に近い非常に小さい値ε未満になると処理を終了する。ノルムの２乗は次式で計算する：
（式１４）‖ｇｔ‖２＝ｇｔＴｇｔ
ステップＳ１０６では、第一の係数計算工程を実行する。第一の係数計算工程では、（式１５）により係数βを計算する：
（式１５）β＝‖ｇｔ‖２／‖ｇｔ−１‖２
ただしｔ＝１のときはβ＝１とする。
【００２９】
ステップＳ１０７では探索ベクトル計算工程を実行する。同工程では（式１６）を用いて、反復回数ｔ回目の探索ベクトルｐｔを計算する：
（式１６）ｐｔ←−ｇｔ−１＋βｐｔ−１
ステップＳ１０８では、（式１６）より得られたｐｔに対して、第一のベクトル修正工程を実行する。同工程は図６を用いて後述する。
【００３０】
ステップＳ１０９ではライン探索工程を実行する。同工程では設計変数の更新のための係数δを決定する。同工程は図３を用いて後述する。
【００３１】
ステップＳ１１０では設計変数ベクトル更新工程を実行する。同工程では設計変数ベクトルを更新する：
（式１７）ｆｔ＋１＝ｆｔ＋δｐｔ
ステップＳ１１１では、設計変数ベクトルに対して第二のベクトル修正工程を実行する。同工程は図７を用いて後述する。
【００３２】
ステップＳ１１２ではｔをｔ＋１に更新し、更新後のｔが予め設定された反復回数を超えたら処理を終了する。
【００３３】
以上が第二の求解工程の流れである。以下、各ステップについて詳述する。
【００３４】
図３を参照してライン探索工程を説明する。同工程には、現在の設計変数値ｆｔ、現在の設計変数値における探索ベクトルｐｔが引数として与えられる。更に、評価関数や拘束条件は適宜与えられているとする。
【００３５】
ステップＳ３０１では最大ステップサイズΔを、次式を用いて計算する：
【外２】
0004101047.tif 000003
【００３６】
ただし
0004101047.tif 000004
また、各変数の初期化を次式のように実行する：
（式１９）ｇ１＝（３．０−５．０１／２）／２，
ｇ２＝（５．０１／２−１）／２，
ｄΔ＝Δ／Ｔｍ，
ａ１＝０．０，
ａ２＝ａ１＋ｄΔ，
ただしＴｍは探索範囲を適当な区間に分割するための正整数であり、例えば１０〜１００等の値を設定する。
【００３７】
ステップＳ３０２では探索範囲絞り込み工程を実行する。同工程は図４を用いて後述する。
【００３８】
ステップＳ３０３では探索範囲絞込み工程で設定されたδと最大ステップサイズΔとを比較し、等しければ処理を終了する。そうでなければステップＳ３０４へ進む。
【００３９】
ステップＳ３０４では極小点探索工程を実行する。同工程は図５を用いて後述する。以上でライン探索工程を終了する。
【００４０】
図４を参照して探索範囲絞込み工程を実行する。
【００４１】
ステップＳ４０１では第２の評価関数値を計算する：
（式２０）Ｌ２１＝Ｌ２（ｘ（ｆｔ＋ａ１ｐｔ），ｆｔ＋ａ１ｐｔ）
（式２１）Ｌ２２＝Ｌ２（ｘ（ｆｔ＋ａ２ｐｔ），ｆｔ＋ａ２ｐｔ）
ただしｘ（ｆｔ＋ａ１ｐｔ）、およびｘ（ｆｔ＋ａ２ｐｔ）は図１０を用いて後述する第一の求解工程を実行することによって得られる。
【００４２】
ステップＳ４０２ではＬ２１とＬ２２とを比較し、Ｌ２１が大きければ処理を終了し、そうでなければステップＳ４０３へ進む。
【００４３】
ステップＳ４０３ではｔを０に、μ（ｔ）を０に初期化する。
【００４４】
ステップＳ４０４ではμ（ｔ）を更新する：
（式２２）μ（ｔ）＝μ（ｔ−１）＋ｄΔ
ステップＳ４０５ではＬ２１とＬ２２を更新する：
（式２３）Ｌ２１＝Ｌ２２
（式２４）Ｌ２２＝Ｌ２（ｘ（ｆｘ＋μ（ｔ）ｐｔ），ｆｘ＋μ（ｔ）ｐｔ）ただしｘ（ｆｘ＋μ（ｔ）ｐｔ）は図１２を用いて後述する第一の求解工程を実行することによって得られる。
【００４５】
ステップＳ４０６ではＬ２１とＬ２２とを比較し、Ｌ２１がＬ２２以下ならステップＳ４０８へ、そうでなければステップＳ４０７へ進む。
【００４６】
ステップＳ４０７では、ｔをｔ＋１に更新したあと、ｔとＴｍを比較し、ｔが大きければステップＳ４０９へ、そうでなければステップＳ４０４へ進む。
【００４７】
ステップＳ４０８では、次式を実行する：
（式２５）ａ１＝μ（ｔ）−２ｄΔ
ａ２＝μ（ｔ）
ステップＳ４０９ではμ（ｔ）とΔを比較し、μ（ｔ）が大きければδにΔを代入し、そうでなければ、δにΔ以外の任意の値を代入して、探索範囲絞込み工程を終了する。
【００４８】
図５を用いて極小点探索工程を説明する。
【００４９】
ステップＳ５０１では次式を計算し
（式２６）ｄａ＝ａ２−ａ１
ｄａが予め設定された微小な正の実数以下であればステップＳ５０４へ進み、そうでなければステップＳ５０２へ進む。
【００５０】
ステップＳ５０２では次式を計算する：
（式２７）ｖ１＝ａ１＋ｇ１ｄａ
ｖ２＝ａ１＋ｇ２ｄａ
更に探索範囲の絞込みを行う：
（式２８）ａ２＝ｖ２，ｉｆ  Ｌ２１＜Ｌ２２
ａ１＝ｖ１，  ｏｔｈｅｒｗｉｓｅ
これは、黄金分割法による区間縮小である。
【００５１】
ステップＳ５０３では、以下の式により第二の評価関数の値を計算し、ステップＳ５０１へ進む。
【００５２】
（式２９）Ｌ２１＝Ｌ２（ｘ（ｆｘ＋ａ１ｐｔ），ｆｘ＋ａ１ｐｔ）
Ｌ２２＝Ｌ２（ｘ（ｆｘ＋ａ２ｐｔ），ｆｘ＋ａ２ｐｔ）
ただしｘ（ｆｔ＋ａ１ｐｔ）、およびｘ（ｆｔ＋ａ２ｐｔ）は図１０を用いて後述する第一の求解工程を実行することによって得られる。
【００５３】
ステップＳ５０４では探索区間の中点ａｃと、その位置での第２の評価関数値Ｌ２ｃを計算する：
（式３０）ａｃ＝（ａ１＋ａ２）／２
（式３１）Ｌ２ｃ＝Ｌ２（ｘ（ｆｘ＋ａｃｐｔ），ｆｘ＋ａｃｐｔ）
ステップＳ５０５ではδの値を決定する：
（式３２）δ＝ａｃ，ｉｆ  Ｌ２ｃ＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝ａ１，ｅｌｓｅ  ｉｆ  Ｌ２１＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝ａ２，    ｅｌｓｅ  ｉｆ  Ｌ２２＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝  Δ，    ｅｌｓｅ  ｉｆ  Ｌ２ｅ＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝  ０，  ｏｔｈｅｒｗｉｓｅ
ただしＭＩＮ［］は引数の中の最小値を返す関数である。以上で極小点探索工程を終了する。
【００５４】
ステップＳ１０４およびステップＳ１０８で実施される第一のベクトル修正工程について、図６を用いて説明する。
【００５５】
ステップＳ６０１では第二の係数計算工程を実行する。ここでは、単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、および係数Ｂの値を計算する。同処理は図８を用いて詳述する。
【００５６】
ステップＳ６０２では第一のベクトル射影工程を実行する。この工程は、ステップＳ６０１で計算した単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、及び係数Ｂを使用し、次式に基づいてベクトルＸを修正する：
0004101047.tif 000005
以上で第一のベクトル修正工程を終了する。
【００５７】
ステップＳ１１１で実施される第二のベクトル修正工程について、図７を用いて説明する。
【００５８】
ステップＳ７０１では第二の係数計算工程を実行する。同処理は図８を用いて詳述する。
【００５９】
ステップＳ７０２では第二のベクトル射影工程を実行する。この工程は、ステップＳ７０１で計算した単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、及び係数Ｂを使用し、次式に基づいてベクトルＸを修正する：
0004101047.tif 000006
以上で第二のベクトル修正工程を終了する。
【００６０】
図８を用いて第二の係数計算工程を説明する。
【００６１】
ステップＳ８０１では法線ベクトル計算工程を実行し、ＮＮ、およびＤＳＴを得る。同工程は図９を用いて後述する。
【００６２】
ステップＳ８０２では、次式を用いて係数Ｂを計算する：
【外３】
0004101047.tif 000007
【００６３】
以上で第二の係数計算工程を終了する。
【００６４】
図９を用いてステップＳ８０１の法線ベクトル計算工程を説明する。本工程では、単位法線ベクトルＮＮと、原点から超平面への垂線の長さＤＳＴを計算する。
【００６５】
ステップＳ９０１では次式の変数を計算する：
（式３６）ＤＳＴ＝−ｃ（０）’
ただしｃ（０）’は次式で計算される：
【外４】
0004101047.tif 000008
【００６６】
（式３６）のＤＳＴを次式で修正する：
（式４０）ＤＳＴ＝ＤＳＴ／Ｄ
上式で与えられるＤＳＴが、原点から超平面への垂線の長さとなる。
【００６７】
ステップＳ９０２では法線ベクトルの計算を実行する。
【００６８】
まず、正規化されていない法線ベクトルＮを（式４１）および（式４２）より計算する：
【外５】
0004101047.tif 000009
【００６９】
（式４２）Ｎ（ｊ）＝ｆｔ（ｊ），    ｏｔｈｅｒｗｉｓｅ
（式４１）を用いて、部分空間に対する単位法線ベクトルＮＮを次式で計算する：
【外６】
0004101047.tif 000010
【００７０】
例えば、等式拘束条件として設計変数ベクトルの要素の総和が１である場合は、ｃ（０）＝−１，ｃ（ｊ）＝１，ｊ＝１，２，…，ｎであるので、単位法線ベクトルの成分は次のように簡単になる：
（式４４）ＮＮ（ｊ）＝｜∧１｜−１／２，ｉｆ  ｊ∈∧１
ただし｜∧１｜は∧１の要素数である。
【００７１】
以上で法線ベクトル計算工程を終了する。
【００７２】
図１０を用いて、第一の求解工程を説明する。
【００７３】
構造が変化すると特性関数値が０となる要素が出現する可能性がある。このような構造に対して構想解析を行うと、係数行列が正則でなくなるので、直接法では解けなくなる。そのために、第１の求解工程として、以下の３通りの方法を用いる：
（１）値が０の特性関数値は０に近い微小な値に置き換えて、直接法で連立１次方程式を解く。
【００７４】
（２）値が０を超える要素のみで連立１次方程式を構成し直して、直接法で解く。
【００７５】
（３）逆行列の計算を必要としない、勾配ベクトルを用いた反復法で解く。これには最急降下法、共役勾配法等が含まれる。
【００７６】
以下では共役勾配法を用いた求解工程を説明する。
【００７７】
ステップＳ１００１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘの初期値、設計パラメタｆの値、境界条件Ｂ１、評価関数Ｌ１が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保し、値を設定する。
【００７８】
ステップＳ１００２では反復回数ｔを１に初期化する。反復回数ｔは、上述の第一の求解工程に対する第１の反復処理に対するものであり、以下のステップＳ１００３からステップＳ１０１０の処理の反復回数が設定値を超えたら終了する。
【００７９】
ステップＳ１００３では状態変数ベクトルに関する第一の評価関数のグラジエントベクトルｇ１，ｔを計算する。ｇ１，ｔは次式で計算される：
（式４５）ｇ１，ｔ＝（∂Ｌ１／∂ｘ１，∂Ｌ１／∂ｘ２，…，∂Ｌ１／∂ｘｎ）Ｔ
∂Ｌ１／∂ｘｊは解析的に解ける場合とそうでない場合がある。解析解がある場合は、その計算工程を利用して（式４５）を得る。そうでない場合は、自動微分という技術を利用して計算することが出来る。自動微分技術は〔文献６〕等によって公知である。
【００８０】
ステップＳ１００４ではｇ１，ｔのノルムの２乗を計算し、その値が０或いは０に近い非常に小さい値のときに処理を終了する。ノルムの２乗は次式で計算する：
（式４６）‖ｇ１，ｔ‖２＝ｇ１，ｔＴｇ１，ｔ
ステップＳ１００５では、（式１２３）よりβを計算する：
（式４７）β＝‖ｇ１，ｔ‖２／‖ｇ１，ｔ−１‖２
ただしｔ＝１のときはβ＝１とする。
【００８１】
ステップＳ１００６では（式１２４）を用いて、探索方向ベクトルｐｔを計算する：
（式４８）ｐ１，ｔ＝−ｇ１，ｔ−１＋βｐ１，ｔ−１
ステップＳ１００７ではライン探索工程を実行する。設計変数の更新のための係数αを決定するものである。構造解析問題を有限要素モデルを用いて解く場合には、αの値は次式で与えられる：
（式４９）α＝ｐ１，ｔＴｇ１，ｔ／（ｐ１，ｔＴＡｐ１，ｔ）
ステップＳ１００８では状態変数ベクトルを更新する：
（式５０）ｘｔ＋１＝ｘｔ＋αｐ１，ｔ
ステップＳ１００９ではｔをｔ＋１に更新し、更新後のｔが予め設定された反復回数を超えたら処理を終了する。
【００８２】
以上で第一の求解工程を終了する。
【００８３】
（実施例）
本実施例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に本発明を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８４】
図１１に示すように、構造部材の存在を可能とする設計領域は、長方形１１０２であり、有限要素法に従って、該領域を縦ｎｙ、横ｎｘに等間隔に分割する。分割にされた部分領域をセルと呼び、左下および右上のセルが（１，１）（ｎｙ，ｎｘ）となるように番号付けを行う。同様に格子点をノードと呼び、左下および右上のノードが（１，１）（ｎｙ＋１，ｎｘ＋１）となるように番号付けを行う。
【００８５】
図中１１０１は支持部材、１１０３は荷重ベクトルである。セル（ｊ，ｋ）には特性関数値ｆ（ｊ，ｋ）が対応する。ここで特性関数値とはセル（ｊ，ｋ）における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本発明における設計変数ベクトルｆの要素である：
（式５１）ｆ＝（ｆ（１，１），ｆ（１，２），…，ｆ（ｎｙ，ｎｘ））Ｔ
同様にノード（ｊ，ｋ）には横方向変位ｕ（ｊ，ｋ）と縦方向変位ｖ（ｊ，ｋ）が対応する。これらは任意の値を取る実数であり、本発明における状態変数ベクトルＵの要素である：
（式５２）Ｕ＝（ｕ（１，１），ｖ（１，１），ｕ（１，２），ｖ（１，２），…，ｕ（ｎｙ＋１，ｎｘ＋１），ｖ（ｎｙ＋１，ｎｘ＋１））Ｔ
同様に剛性マトリクスをＡ、加重ベクトルをｂと書けば、有限要素法の良く知られた結果として、状態変数ベクトルＵは、次式で与えられる評価汎関数の最適化問題の解として与えられる：
（式５３）Ｌ１＝（１／２）ＵＴＡＵ−ｂＴＵ
より具体的には、以下の線型方程式の解として与えられることが知られている：
（式５４）ＡＵ−ｂ＝０
第一の求解工程において計算されるｇ１，ｔは（式５４）の左辺として与えられる。（式５４）は連立一次方程式であるので、求解法として、通常、逆行列に基づく直接解法が採用される。しかし構造最適化問題においては、剛性マトリクスＡは設計変数ベクトルｆの関数であり、ｆに０の値を取る要素があれば、剛性マトリクスＡはランク落ちし正則でなくなってしまう。その結果直接解法で解けないという事態に陥るので、本発明で共役勾配法を採用しているのは、このような理由による。
【００８６】
設計変数ベクトルｆに対する評価関数Ｌ２は総歪エネルギーで定義される。
【００８７】
（式５５）Ｌ２＝（１／２）ＵＴＡＵ
第二の求解工程で必要となるｇｔは次式で計算できる：
（式５６）ｇｔ＝∂Ｌ２／∂ｆ（ｅ）＝−（１／２）ＵｅＴＡｅＵｅ
ただしＵｅ、Ａｅは、それぞれ、要素ｅに対応するノードにおける変位ベクトル、該変位ベクトルに対応する要素剛性マトリクスである。
【００８８】
図１２に本実施例に対する計算結果を示す。図中、黒い領域が構造部材が存在する領域である。
【００８９】
本実施例において、構造部材の設計領域のアスペクト比は縦／横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。
【００９０】
図１２に示した計算結果は、このような解析解とよく一致していることがわかる。
【００９１】
尚、本発明は、単一の機器からなる装置に適用しても、複数の機器から構成されるシステムに適用してもよい。また、上述した実施形態の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体を、装置あるいはシステムに供給し、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって達成してもよい。
【００９２】
更に、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって、上述した実施形態の機能を直接実現するばかりでなく、そのプログラムコードの指示に基づいて、コンピュータ上で稼動しているＯＳなどの処理により、上述の機能を実現される場合も含まれる。
【００９３】
これらの場合、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。
【００９４】
以下、上記実施形態に係る本発明の特徴を整理する。
【００９５】
特徴１．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計方法において、
前記第２の求解工程は、
前記設計ベクトルに関する前記第２の評価関数の勾配ベクトルを計算する勾配ベクトル計算工程と、
前記勾配ベクトルのノルムの値に基づいて第１の係数を計算する第１の係数計算工程と、
該第１の係数に基づいて探索ベクトルを計算する探索ベクトル計算工程と、
第２の係数を計算する第２の係数計算工程と、
前記第２の係数に基づいて設計変数ベクトルを更新する設計変数ベクトル更新工程とを備え、
前記第２の係数計算工程は前記第１の求解工程を含み、該第１の求解工程は勾配ベクトルに基づく反復解法であり、当該反復中は前記状態変数ベクトルを初期化しないことを特徴とする最適設計方法。
【００９６】
特徴２．
前記ライン探索工程は、
前記拘束条件に基づいて最大ステップサイズを決定する最大ステップサイズ計算工程と、
前記最大ステップサイズより決定される探索区間の中から極小点を決定する探索範囲絞込み工程と、
極小点における前記第２の評価関数値と探索範囲の終点における前記第２の評価関数値との比較に基づいて極小点を決定する極小点決定工程とを備えることを特徴とする特徴１に記載の最適設計方法。
【００９７】
特徴３．
前記第２の求解工程は、前記設計変数ベクトルより、前記第２の評価関数に寄与しない構造要素に対応する成分を消去する消去工程を含むことを特徴とする特徴１に記載の最適設計方法。
【００９８】
特徴４．
前記消去工程では、前記設計変数ベクトルの各成分の値を増減させたとき、当該設計変数ベクトルに関する前記第２の評価関数の感度の値が常に０である成分を消去することを特徴とする特徴１に記載の最適設計方法。
【００９９】
特徴５．
前記消去工程は、前記第２の求解工程で計算される感度ベクトルが０である成分に対応する要素であって、かつ構造要素の存在率が０でない要素に対して行われることを特徴とする特徴１に記載の最適設計方法。
【０１００】
特徴６．
前記消去工程は、前記第２の求解工程の反復処理の所定回に１回行われることを特徴とする特徴１に記載の最適設計方法。
【０１０１】
特徴７．
前記設計変数ベクトルは、各要素における構造部材の存在率であることを特徴とする特徴１に記載の最適設計方法。
【０１０２】
特徴８．
前記勾配ベクトル計算工程により計算された勾配ベクトルを拘束条件に基づいて修正する勾配ベクトル修正工程を備えることを特徴とする特徴１に記載の最適設計方法。
【０１０３】
特徴９．
前記探索ベクトル計算工程により計算された探索ベクトルを拘束条件に基づいて修正する探索ベクトル修正工程を備えることを特徴とする特徴１に記載の最適設計方法。
【０１０４】
特徴１０．
前記設計変数ベクトル更新工程により更新された設計変数ベクトルを拘束条件に基づいて修正する設計変数ベクトル修正工程を備えることを特徴とする特徴１に記載の最適設計方法。
【０１０５】
特徴１１．
前記勾配ベクトルのノルムの値に基づいて前記第２の求解工程を終了するか否かを判定する収束判定工程を備えることを特徴とする特徴１に記載の最適設計方法。
【０１０６】
特徴１２．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解手段と、
前記第１の求解手段で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解手段とを備える最適設計装置において、
前記第２の求解手段は、
前記設計ベクトルに関する前記第２の評価関数の勾配ベクトルを計算する勾配ベクトル計算手段と、
前記勾配ベクトルのノルムの値に基づいて第１の係数を計算する第１の係数計算手段と、
該第１の係数に基づいて探索ベクトルを計算する探索ベクトル計算手段と、
第２の係数を計算する第２の係数計算手段と、
前記第２の係数に基づいて設計変数ベクトルを更新する設計変数ベクトル更新手段とを備え、
前記第２の係数計算手段は前記第１の求解手段を含み、該第１の求解手段は勾配ベクトルに基づく反復解法を行ない、当該反復中は前記状態変数ベクトルを初期化しないことを特徴とする最適設計装置。
【０１０７】
特徴１３．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計プログラムにおいて、
前記第２の求解工程は、
前記設計ベクトルに関する前記第２の評価関数の勾配ベクトルを計算する勾配ベクトル計算工程と、
前記勾配ベクトルのノルムの値に基づいて第１の係数を計算する第１の係数計算工程と、
該第１の係数に基づいて探索ベクトルを計算する探索ベクトル計算工程と、
第２の係数を計算する第２の係数計算工程と、
前記第２の係数に基づいて設計変数ベクトルを更新する設計変数ベクトル更新工程とを備え、
前記第２の係数計算工程は前記第１の求解工程を含み、該第１の求解工程は勾配ベクトルに基づく反復解法であり、当該反復中は前記状態変数ベクトルを初期化しないことを特徴とする最適設計プログラム。
【０１０８】
【発明の効果】
以上説明したように、本発明によれば、個別の問題に依存する経験則を利用すること無しに、トポロジーを含む構造最適設計を実行することができる。
【図面の簡単な説明】
【図１】第２の求解工程の処理手順を示すフローチャートである。
【図２】本実施形態に係る装置のブロック構成図である。
【図３】ライン探索工程の処理手順を示すフローチャートである。
【図４】探索範囲絞込み工程の処理手順を示すフローチャートである。
【図５】極小点探索工程の処理手順を示すフローチャートである。
【図６】第１のベクトル修正工程の処理手順を示すフローチャートである。
【図７】第２のベクトル修正工程の処理手順を示すフローチャートである。
【図８】係数Ｂの計算工程の処理手順を示すフローチャートである。
【図９】法線ベクトル計算工程の処理手順を示すフローチャートである。
【図１０】第１の求解工程の処理手順を示すフローチャートである。
【図１１】実施例の問題設定の説明図である。
【図１２】実施例の計算結果を示す図である。
【図１】
0004101047.tif 000011
【図２】
0004101047.tif 000012
【図３】
0004101047.tif 000013
【図４】
0004101047.tif 000014
【図５】
0004101047.tif 000015
【図６】
0004101047.tif 000016
【図７】
0004101047.tif 000017
【図８】
0004101047.tif 000018
【図９】
0004101047.tif 000019
【図１０】
0004101047.tif 000020
【図１１】
0004101047.tif 000021
【図１２】
0004101047.tif 000022