(57)【要約】【課題】全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量を軽減する構造物の最適設計方法を提供する。【解決手段】第１及び第２の求解工程を有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法で、状態変数ベクトルを各ノードにおける変位、設計変数ベクトルを各要素における構造部材の存在率とし、第１の求解工程は、その１ステップとして第２の求解工程を含み、記憶された設計変数ベクトル及び状態変数ベクトルを読み出して設計変数ベクトルを更新し、第２の求解工程は、記憶された設計変数ベクトル及び状態変数ベクトルを読み出して状態変数ベクトルを更新するＳ３０４。第２の求解工程の第２の評価汎関数が残差ベクトルのノルムより構成されておりＳ３０５、状態変数ベクトルが第２の求解工程が開始される時点で初期化されない。【選択図】    図４
2004310375.tif 000002
【特許請求の範囲】
【請求項１】
状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解工程と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解工程とを有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法において、
前記状態変数ベクトルを各ノードにおける変位、前記設計変数ベクトルを各要素における構造部材の存在率とし、
前記第１の求解工程は、その１ステップとして前記第２の求解工程を含み、第１の記憶手段に記憶された設計変数ベクトル及び状態変数ベクトルを読み出して設計変数ベクトルを更新し、更新された設計変数ベクトルを該第１の記憶手段に記憶する設計変数更新工程を含み、
前記第２の求解工程は、第２の記憶手段に記憶された設計変数ベクトル及び状態変数ベクトルを読み出して状態変数ベクトルを更新し、更新された状態変数ベクトルを該第２の記憶手段に記憶する状態変数更新工程を含み、前記第２の求解工程の第２の評価汎関数が残差ベクトルのノルムより構成されており、前記状態変数ベクトルが前記第２の求解工程が開始される時点で初期化されないことを特徴とする構造最適設計方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は設計パラメータの最適化問題の解法である構造物の最適設計、特に構造部材のトポロジーを含む形状最適化のための自動設計技術に関するものである。
【０００２】
【従来の技術】
構造トポロジーの最適設計とは、与えられた条件の下で、構造部材の最適なトポロジー及び形状寸法を決定する問題である。以下、構造部材のトポロジー及び形状寸法を設計変関数といい、上記決定問題を設計変関数の最適化問題という。変関数という理由は、トポロジー及び形状寸法が３次元空間の関数になっているからである。設計変関数の最適化問題では、各設計変関数の値に対して、状態変関数の最適化問題を解かなければならない。この意味から、構造トポロジーの最適設計は、内側に状態変関数の最適化問題を持ち、外側に設計変関数の最適化問題を持つ２重構造最適化問題と捉えることが出来る。
【０００３】
内側の状態変関数最適化問題では、従来技術の蓄積により、空間を有限個の要素に分割するという考え方が採用される。特に構造部材の歪エネルギーを評価汎関数としている問題では、その解析手法として有限要素法が一般的である。有限要素法の解法としては一次方程式に対する直接法が採用されている。
【０００４】
一方、設計変関数の最適化問題に対しては、大別すると、以下に示す３種類の方法が提供されている（非特許文献１、或いは、非特許文献２参照）。
１．Ｅｖｏｌｕｔｉｏｎａｒｙ法（以下、Ｅ法）
２．Ｈｏｍｏｇｅｎｉｚａｔｉｏｎ法（以下、Ｈ法）
３．Ｍａｔｅｒｉａｌ  ｄｉｓｔｒｉｂｕｔｉｏｎ法（以下、ＭＤ法）、或いは密度（Ｄｅｎｓｉｔｙ）法（以下、Ｄ法）
Ｅ法では、空間を分割することによって得られる部分空間のそれぞれをセルと称し、その生成と消滅を適当な規則によって繰り返す。構造部材は、最終的に存在しているセルの集合として与えられる。セルが存在するか否かという２つの状態のみを許すことにより、明確な構造部材が得られる。また、評価汎関数の微分情報を用いないので、局所最適解にトラップされないことから、評価汎関数が多峰性の場合に有効である。特許文献１では、Ｅ法の一種である遺伝的探索法を用いた、骨組構造部材の最適化設計装置が提供されている。この最適化設計装置では、蓄積されたノウハウに基づく試行計算を必要とし、従って膨大な設計変数が存在する実設計問題に対応することが出来なかった従来技術の課題を、以下の方法によって解決している。即ち、骨組部材断面寸法などの離散設計変数データの近似式を使用する近似最適化計算装置と、該設計変数データを使用する詳細最適化計算装置を設け、これら２つの計算装置を結合して骨組構造の最適設計装置を構成している。
【０００５】
Ｈ法は、分割された各部分領域に位置する構造要素に、更に微細な構造を仮定し、連続値を取る設計変関数を新たに導入することによって、感度解析の採用を可能にしている。ここで感度解析とは、設計変関数に関する評価汎関数の微分情報を利用した解析手法のことであり、感度解析が可能になれば、勾配法のような反復解法を用いることが出来、Ｅ法のような総当り的手法に比べて、少なくとも局所最適解の探索に係る計算時間は大幅に短縮される（非特許文献３参照）。
【０００６】
ＭＤ法或いはＤ法は、構造部材の存在率を示す０から１の範囲の実数を各要素に割当てることによって、構造部材のトポロジーや形状寸法変化を表現する方法である。構造部材が存在するか否かという離散的な情報を存在率という連続値に置き換えることによって感度解析を可能にしたという意味でＨ法と同様のものであるが、Ｈ法よりパラメータ数が少ない分、モデル化が容易であり適用範囲も広い。
【０００７】
非特許文献４には、Ｄ法による構造物の位相最適化手法が開示されている。本手法の特徴は以下のとおりである。ボクセル有限要素法（空間を等間隔に分割）を用いているため、あらゆる要素に対する要素剛性マトリクスが同じである。従って、要素剛性マトリクスを予め１度だけ計算しておけば、以後の計算に利用できる。更に、要素が規則正しく配置されているため、各要素の節点番号情報を記憶する必要がない。大規模連立一次方程式を解くために、前処理付共役勾配法とＥｌｅｍｅｎｔ−ｂｙ−Ｅｌｅｍｅｎｔ法を組み合わせて用いたことにより、全体剛性マトリクスを組み立てることなく解が求められるので、必要とするメモリ容量が少なくて済む。
【０００８】
均質化法では、１要素に対して６個の設計変数（３次元の場合）が必要になる。更に設計変数が変化するたびに要素剛性マトリクスを再計算しなければならない。一方、構造部材の存在率を密度比で表現する密度法を採用することによって、１要素に対して１つの設計変数でよい。また設計変数の変化が要素剛性マトリクスに影響を与えない。
【０００９】
【特許文献１】
特開平１１−３１４６３１号
【非特許文献１】
Ｓ．  Ｂｕｌｍａｎ，  Ｊ．  Ｓｉｅｎｚ，  Ｅ．  Ｈｉｎｔｏｎ：  ”Ｃｏｍｐａｒｉｓｏｎｓ  ｂｅｔｗｅｅｎ  ａｌｇｏｒｉｔｈｍｓ  ｆｏｒ  ｓｔｒｕｃｔｕｒａｌ  ｔｏｐｏｌｏｇｙ  ｏｐｔｉｍｉｚａｔｉｏｎ  ｕｓｉｎｇ  ａ  ｓｅｒｉｅｓ  ｏｆ  ｂｅｎｃｈｍａｒｋ  ｓｔｕｄｉｅｓ，”  Ｃｏｍｐｕｔｅｒｓ  ａｎｄ  Ｓｔｒｕｃｔｕｒｅｓ，  ７９，  ｐｐ．１２０３−１２１８  （２００１）．
【非特許文献２】
Ｙ−Ｌ．  Ｈｓｕ，  Ｍ−Ｓ．  Ｈｓｕ，  Ｃ−Ｔ．  Ｃｈｅｎ：  ”Ｉｎｔｅｒｐｒｅｔｉｎｇ  ｒｅｓｕｌｔｓ  ｆｒｏｍ  ｔｏｐｏｌｏｇｙ  ｏｐｔｉｍｉｚａｔｉｏｎ  ｕｓｉｎｇ  ｄｅｎｓｉｔｙ  ｃｏｎｔｏｕｒｓ，”  Ｃｏｍｐｕｔｅｒｓ  ａｎｄ  Ｓｔｒｕｃｔｕｒｅｓ，  ７９，  ｐｐ．１０４９−１０５８  （２００１）．
【非特許文献３】
山川宏：“最適化デザイン，”計算力学とＣＡＥシリーズ９，培風館（１９９６）
【非特許文献４】
藤井，鈴木，大坪：“ボクセル有限要素法を用いた構造物の位相最適化，”Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｆ  ＪＳＣＥＳ，  Ｐａｐｅｒ  Ｎｏ．２０００００１０  （２０００）．
【００１０】
【発明が解決しようとする課題】
しかしながら、上記従来技術には以下に述べるような問題点があった。
【００１１】
一般に構造最適化問題は、設計変数ベクトルの最適化問題の各反復ステップ毎に、状態変数ベクトルの最適化問題を含む２重の最適化問題として定式化される。設計変数ベクトルの最適化問題を外側の最適化、状態変数ベクトルの最適化問題を内側の最適化と呼ぶと、内側の最適化は、設計変数ベクトルをパラメータとし、即ち設計変数ベクトルを固定し、状態変数ベクトルを求める問題である。これは通常、構造解析と呼ばれ、有限要素法によって線形方程式の解法を用いて解くことが出来る。
【００１２】
ところが、構造が変化し、ある領域の構造部材が存在しなくなる、例えば部材に穴が空いてしまうと、対応する要素の設計変数が０となり、結果的に要素のヤング率が０となる。すると、上記線形方程式の係数行列がフルランクでなくなり、逆行列が計算できないという理由で、直接法では解けなくなる。
【００１３】
そのため従来技術では、設計変数ベクトルが０となる場合に以下のように対処してきた。
【００１４】
係数行列がフルランクになるように、即ちランク数と方程式の数が等しくなるように、方程式を構成し直す。
【００１５】
設計変数ベクトルの値を正確に０とはせずに０に近い微小な値で置き換える。
【００１６】
しかるに、上記前者の方法は、設計変数ベクトルが更新される毎に方程式の更新を行う可能性が生じ、計算時間がかかる。
【００１７】
また、上記後者の方法、即ち要素の設計変数を微小な値に設定するということは、物理的には薄い膜、或いは弱い部材が存在することに相当する。つまり従来法では、物質がない部分を正確に表現していないので、得られた計算結果の精度に疑問が残る。
【００１８】
本発明は、上記従来の問題点に鑑み、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量を軽減する構造物の最適設計方法を提供する。
【００１９】
【課題を解決するための手段】
上記課題を解決するために、本発明の構造最適設計方法は、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解工程と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解工程とを有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法において、前記状態変数ベクトルを各ノードにおける変位、前記設計変数ベクトルを各要素における構造部材の存在率とし、前記第１の求解工程は、その１ステップとして前記第２の求解工程を含み、第１の記憶手段に記憶された設計変数ベクトル及び状態変数ベクトルを読み出して設計変数ベクトルを更新し、更新された設計変数ベクトルを該第１の記憶手段に記憶する設計変数更新工程を含み、前記第２の求解工程は、第２の記憶手段に記憶された設計変数ベクトル及び状態変数ベクトルを読み出して状態変数ベクトルを更新し、更新された状態変数ベクトルを該第２の記憶手段に記憶する状態変数更新工程を含み、前記第２の求解工程の第２の評価汎関数が残差ベクトルのノルムより構成されており、前記状態変数ベクトルが前記第２の求解工程が開始される時点で初期化されないことを特徴とする。
【００２０】
【発明の実施の形態】
以下、本発明の構造最適設計法を実現する構造最適設計装置の構成例及び動作例を説明する。
【００２１】
＜本実施形態の構造最適設計装置の構成例＞
図１に、本実施形態の構造最適設計装置の構成例を示す。
【００２２】
図中、２０１は演算・制御用のＣＰＵ、２０２は表示装置、２０３は入力装置、２０４は１次記憶装置、２０５は２次記憶装置、２０６は通信装置、２０７はバスラインである。本実施形態における構成要素は、プログラムとして予め２次記憶装置２０５に格納されており、入力装置２０３或いは通信装置２０６等からのコマンド入力によって１次記憶装置にロードされ、ＣＰＵ２０１により実行されるものである。
【００２３】
図２に、本実施形態の構造最適設計法を実現するために、１次記憶装置２０４及び２次記憶装置２０５に記憶されるデータ及びプログラムの記憶例を示す。
【００２４】
１次記憶装置２０４は、ＲＡＭやＲＯＭからなり、データ記憶領域２１４とプログラム記憶領域２２４を有している。図２には、変更されずに常駐するＯＳやＢＩＯＳのプログラム、変更されないＯＳやＢＩＯＳ用のパラメータなどは図示されていない。データ記憶領域２１４には、本実施形態の構造最適設計を実現するための各種のデータが一時記憶される。この中で再度利用されるデータ／反復して使用されるデータは、一時記憶されることが必須である。保存・保持の必要が無いデータは、一時記憶は必須ではない。
【００２５】
例えば、データ記憶領域２１４には、ｘ（状態変数ベクトル）２１４ａ、ｆ（設計変数ベクトル）２１４ｂ、Ｌ１（状態変数ベクトルの評価関数）２１４ｃ、Ｌ２（設計変数ベクトルの評価関数）２１４ｄ、Ｂ１（状態変数のための境界条件）２１４ｅ、Ｗ０（設計変数のための重量条件）２１４ｆ、Ｘ（０）やｆ（０）（初期値）２１４ｇ、λ（ラグランジェ未定定数）２１８ｈ、Ａ（要素合成行列）２１４ｉ、ｋ（設計変数変更の反復回数）２１４ｊ、Ｃ（感度係数）２１４ｋ、ｒ（残差ベクトル）１２４ｍ、ｐ（探索方向ベクトル）２１５ｎ、ｒ（０）やｐ（０）（初期値）２１４ｐ、ｔ（状態変数変更の反復回数）２１４ｑ、ｂ（節点力ベクトル）２１４ｒ、その他のデータ２１４ｓなどが含まれる。一方、プログラム記憶領域２２４には、２次記憶装置２０５或は通信装置２０６からロードされたプログラムが記憶されて、ＣＰＵ２０１により実行される。
【００２６】
２次記憶装置２０５は、例えば、フロッピーディスクやＣＤなどの外部大容量メモリであり、着脱可能であればなお好ましい。２次記憶装置２０５も、データ記憶領域２１５とプログラム記憶領域２２５を有している。
【００２７】
データ記憶領域２１５には、詳細には示さないが、各種の初期値２１５ａや各種の条件２１５ｂ、あるいは２重最適値問題のデータベース２１５ｃが記憶されている。データベース２１５ｃには、標準化されたマトリクスなどが記憶されていて使用されてもよい。
【００２８】
プログラム記憶領域２２５には、本実施形態で設計変数ベクトルにおける求解法として使用される、最適性規準法モジュール２２５ａ、逐次凸関数近似法モジュール２２５ｂ、逐次線形計画法モジュール２２５ｃ、その他の方法のモジュール２２５ｄが記憶され、状態変数ベクトルにおける求解法として使用される、共役残差法モジュール２２５ｅ、ＧＣＲ法モジュール２２５ｆ、ＧＣＲ（ｋ）法２２５ｇ、Ｏｒｔｈｏｍｉｎ（ｋ）法２２５ｈ、ＧＭＲＥＧ（ｋ）法２２５ｉ、その他の方法のモジュール２２５ｊが記憶されている。このプログラム記憶領域２２５に記憶されたプログラムから選択された方法に対応するモジュールが１次記憶装置２０４のプログラム記憶領域２２４にロードされて、ＣＰＵ２０１により実行され本実施形態の構造最適設計法が実現する。
【００２９】
すなわち、本構造最適設計装置は、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解手段と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解手段とを有し、２重最適化問題として定式化される構造最適設計問題を解く構造最適設計装置であって、前記状態変数ベクトルが各ノードにおける変位、前記設計変数ベクトルが各要素における構造部材の存在率であり、前記第２の求解手段の第２の評価汎関数は残差ベクトルのノルムより構成されており、前記状態変数ベクトルが前記第２の求解工程が開始される時点で初期化されない。前記第１の求解手段は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかを実行する。前記第２の求解手段は、共役残差法、ＧＣＲ法、ＧＣＲ（ｋ）法、Ｏｒｔｈｏｍｉｎ（ｋ）法、ＧＭＲＥＳ（ｋ）法、及びこれらの派生法を含む方法のいずれかを実行する。
【００３０】
＜本実施形態の構造最適化問題の定式化＞
以下の説明のために、構造最適化問題の定式化を行う。
【００３１】
有限要素法による定式化により変関数は有限次元ベクトルで表現されるものとすると、変関数の評価汎関数は、変数ベクトルの評価関数となる。以下、有限次元ベクトルで表現されているとして記述する。
【００３２】
状態変数ベクトルｘ、設計変数ベクトルｆをそれぞれ列ベクトルとして以下のように書く。
【００３３】
ｘ  ＝  （ｘ１，ｘ２，…，ｘｍ）Ｔ  …（式１）
ｆ  ＝  （ｆ１，ｆ２，…，ｆｎ）Ｔ  …（式２）
ここで、Ｔは転置を表す。ｘはｍ次元ベクトル、ｆはｎ次元ベクトルである。構造最適化問題において、状態変数ベクトルは各ノードの変位、設計変数ベクトルは各要素内における構造物の存在率より構成される。
【００３４】
例えば、２次元面内での構造最適化問題において、設計領域を縦ｎｙ、横ｎｘに分割すると、
要素数ｎは（ｎｘ×ｎｙ）、ノード数は（ｎｘ＋１）×（ｎｙ＋１）である。
【００３５】
状態変数ベクトルは各ノードの縦方向及び横方向変位の組なので、ｘの次元ｍは２×（ｎｘ＋１）×（ｎｙ＋１）である。
【００３６】
ｘの境界条件は通常、変位拘束条件として与えられる。それをＢ１と記述する。
Ｂ１（ｘ）＝０  …（式３）
状態変数ベクトルｘ及び設計変数ベクトルｆに関する評価関数をそれぞれＬ１、Ｌ２とすると、次式のように定義される。
Ｌ１  ：＝  Ｌ１（ｘ，ｆ）  ＝  （ｂ  −  Ａｘ）Ｔ（ｂ  −  Ａｘ）  …（式４）
Ｌ２  ：＝  Ｌ２（ｆ，ｘ）  ＝  （１／２）ｘＴＡｘ  …（式５）
ここで、Ａは全体剛性行列、ｂは節点力ベクトルである。ただし、ｂは予め与えられるものであり、Ａは設計変数ベクトルｆの関数である。
【００３７】
全体剛性行列は要素ｊに対する要素剛性行列Ａｊの重ね合わせで構成することが出来る。平面歪問題の場合、設計変数ベクトル、即ち構造要素の存在率ｆｊを考慮した要素剛性行列を以下に示す。ただし、２次元有限要素として４ノード要素を採用し、要素変位ベクトルｘｊを次のように構成する。
ｘｊ  ＝  （ｕ１，  ｖ１，  ｕ２，  ｖ２，…  ，  ｕ４，  ｖ４）  …（式６）
ここで、ｕｋ、ｖｋはノードｋの水平方向、及び垂直方向の変位である。
【００３８】
（式６）の要素変位ベクトルに対応する要素剛性行列は、次式のように与えられる。
【００３９】
【数１】
2004310375.tif 000003
【００４０】
上式中
ａｉ，ｊ  ＝  （λ＋２μ）＜∂ｘｅｉ，∂ｘｅｊ＞    ＋  μ＜∂ｙｅｉ，∂ｙｅｊ＞  …（式８）
ｂｉ，ｊ  ＝              λ＜∂ｘｅｉ，∂ｙｅｊ＞    ＋  μ＜∂ｘｅｊ，∂ｙｅｉ＞  …（式９）
ｃｉ，ｊ  ＝              λ＜∂ｘｅｊ，∂ｘｅｉ＞    ＋  μ＜∂ｘｅｉ，∂ｙｅｊ＞  …（式１０）
ｄｉ，ｊ  ＝  （λ＋２μ）＜∂ｙｅｉ，∂ｙｅｊ＞    ＋  μ＜∂ｘｅｉ，∂ｘｅｊ＞  …（式１１）
ここで、ｅｊはノードｊに対する基底ベクトル、＜  ，  ＞はベクトルの内積を表す。∂ｘはｘに関する偏微分、∂ｙはｙに関する偏微分である。ξは前述の（非特許文献４）に従い、２とする。また、λ、μはラメ定数と呼ばれる材料定数であり、ヤング率Ｅ及びポアソン比νより次式で計算される。
【００４１】
λ  ＝  Ｅν／（１＋ν）（１−２ν）  …（式１２）
μ  ＝  Ｅ／｛２（１＋ν）｝  …（式１３）
（式７）より容易に分かるように、ある要素の存在率が０となると要素剛性行列Ａｊのすべての成分が０となる。全体剛性行列は要素剛性行列の重ね合わせで構成される。従って、ノードｊを含む要素の設計変数ｆｊが０になると、全体剛性行列のノードｊに対応する行及び列の全ての成分が０になる。このように、構造最適化問題において全体剛性行列は一般に特異な行列となる。
【００４２】
設計変数に関する拘束条件は、構造最適設計の場合、総重量の上限値Ｗ０とする。
【００４３】
【数２】
2004310375.tif 000004
【００４４】
また、設計変数には値域に対する拘束条件が以下のように定められている。
【００４５】
０  ≦  ｆｊ  ≦  １    ｊ  ＝  １，  …  ，  ｎ    …（式１５）
以上の表記より、最適設計問題は、拘束条件（式１４）及び（式１５）の下で、（式５）を極小化する問題として定式化される。
【００４６】
また、状態変数ｘについては、拘束条件（式３）の下で（式４）を極小化する解として得られる。
【００４７】
＜本実施形態の構造最適設計装置の動作例＞
上記定式化に基づき、本実施形態の処理手順例を説明する。
【００４８】
図３には、本実施形態の流れ図が示されている。
【００４９】
図中、ステップＳ１０１はシミュレーション対象となる系の諸元を読み込み、変数の初期化を行う処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値としてｘ（０）、ｆ（０）、境界条件Ｂ、評価関数Ｌ１、Ｌ２が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。感度係数の初期値Ｃｊ（０）は次式で計算する。
【００５０】
Ｃｊ（０）  ＝  −  （ｘｊ（０））Ｔ（ξ（ｆｊ（０））ξ−１Ａｊ）ｘｊ（０）    …（式１６）
ステップＳ１０２では、反復回数ｋを１に初期化する。反復回数ｋは、上述の設計変数の求解工程に対する反復処理に対するものであり、予め定められた終了条件（評価関数Ｌ１，Ｌ２の条件）を満たしたら終了する。
【００５１】
ステップＳ１０３では、ラグランジェ未定定数λの更新を行う。ステップＳ１０４では、設計変数ベクトルｆの更新を行う。ステップＳ１０５では、状態変数ベクトルｘの更新を行う。ステップＳ１０６では、感度係数の計算を行う。
【００５２】
上記ステップＳ１０３及びＳ１０４は、拘束条件付最適化問題の解法として採用する方法によって異なる。例えば、この設計変数ベクトルｆ（ｋ）を算出する求解工程は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかにより実現される。以下、解法として最適性規準法を用いる場合と、逐次凸関数近似法を用いる場合について詳述する。
【００５３】
（最適性規準法）
まず、最適性規準法を用いる場合について説明する。構造最適化問題に対する最適性規準法は（非特許文献４）及び（非特許文献５）等によって公知である。
【００５４】
最適性規準法では、（式４）及び（式５）をラグランジェ未定定数λを用いて、１つの式にまとめる。
【００５５】
Ｌ（ｆ，λ）  ＝  Ｌ２（ｆ，ｘ）  −  λ（Ｗ（ｆ）  −  Ｗ０）  …（式１７）
ステップＳ１０３で実行されるλの更新は、次式を用いて行われる。
【００５６】
λ（ｋ＋１）  ＝  ｍｉｎ［０，［（１／Ｗ０）Ｗ（ｆ（ｋ））］  αλ（ｋ）］  …（式１８）
ただし、上付きのｋ，ｋ＋１は反復回数を表す。また、αはべき乗係数で、通常０．８５程度の値を設定する。
【００５７】
ステップＳ１０４で実行される設計変数の更新は次式による。
【００５８】
2004310375.tif 000005
ただし、ζは設計変数更新の変動幅の制約値で、０．３程度の値を取る。また、ｓｊ（ｋ）は次式で与えられる。
【００５９】
ｓｊ（ｋ）  ＝  ［（１／λ（ｋ−１））Ｃｊ（ｋ−１）］αｆｊ（ｋ−１）  …（式２０）
ここで、Ｃｊ（ｋ）は第２の評価関数Ｌ２の、設計変数ｆｊ（ｋ）に関する感度係数と呼ばれる量であり、ステップＳ１０６の処理で計算されるものである。
【００６０】
更に（式１９）で計算された設計変数ｆｊ（ｋ）は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、例えば１０−３程度の値である。
【００６１】
ステップＳ１０５では、状態変数ベクトルｘ（ｋ）を更新する。
【００６２】
状態変数ベクトルｘ（ｋ）は、（式４）の評価関数の極小化の過程で得られる。この解法は共役残差法と呼ばれており、後述する。
【００６３】
ステップＳ１０６では、感度係数Ｃｊ（ｋ）を次式により計算する。
【００６４】
Ｃｊ（ｋ）  ＝  −  （ｘｊ（ｋ））Ｔ（ξ（ｆｊ（ｋ））ξ−１Ａｊ）ｘｊ（ｋ）  …（式２１）
（逐次凸関数近似法）
次に、逐次凸関数近似法を用いた場合の、ステップＳ１０３及びＳ１０４の処理について詳述する。逐次凸関数近似法では設計変数ベクトルの各成分に対してスケーリングを行うことによって、数値計算誤差を避けることができるが、以下では説明を簡単にするため、スケーリングについては記述していない。逐次凸関数近似法の詳細は、（参考文献１）及び（参考文献２）等によって公知である。
【００６５】
（参考文献１）藤井：“パソコンで解く構造デザイン”，  丸善（２００２年４月）
（参考文献２）Ｃ．  Ｆｌｅｕｒｙ  ：  ”ＣＯＮＬＩＮ，  ａｎ  ｅｆｆｉｃｉｅｎｔ  ｄｕａｌ  ｏｐｔｉｍｉｚｅｒ  ｂａｓｅｄ  ｏｎ  ｃｏｎｖｅｘ  ａｐｐｒｏｘｉｍａｔｉｏｎ  ｃｏｎｃｅｐｔｓ，”  Ｓｔｒｕｃｔｕｒａｌ  Ｏｐｔｉｍｉｚａｔｉｏｎ，  １，  ｐｐ．８１−８９  （１９８９）
ステップＳ１０３では、ラグランジェ未定定数を次式により更新する。
【００６６】
λ（ｋ）  ＝  λ（ｋ−１）  ＋  Δλ（ｋ）  …（式２２）
ただし、初期値λ（０）は例えば０．０１程度の値を設定する。
【００６７】
Δλ（ｋ）は次の方程式の解として与えられる。
【００６８】
Δλ（ｋ）  ＝  ２×（  １  −  Ｗ０／Ｗ（ｆ（ｋ−１））  ）λ（ｋ−１）    …（式２３）
ステップＳ１０４では、設計変数ベクトルの更新を行う。
【００６９】
2004310375.tif 000006
ただし、ｇ（ｋ）  ＝  −（ｆｊ（ｋ−１））２Ｃｊ（ｋ−１）  …（式２５）
更に、（式２４）で計算された設計変数ｆｊ（ｋ）は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、１０−３程度の値である。
【００７０】
ステップＳ１０５及びＳ１０６の処理は前述のとおりである。
【００７１】
（共役残差法）
ステップＳ１０５において状態変数ベクトルｘ（ｋ）を算出する求解工程は、共役残差法、ＧＣＲ法、ＧＣＲ（ｋ）法、Ｏｒｔｈｏｍｉｎ（ｋ）法、ＧＭＲＥＳ（ｋ）法、及びこれらの派生法を含む方法のいずれかにより実現される。以下、それらを代表して共役残差法を用いた処理を図４を用いて説明する。
【００７２】
ステップＳ３０１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘ（ｋ）の初期値、設計変数ベクトルｆ（ｋ）の値、変位拘束条件Ｂ、評価関数Ｌ１が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。また、状態変数ベクトルｘ、残差ベクトルｒ、探索方向ベクトルｐを初期化する。それぞれの初期値をｘ（０）、ｒ（０）、ｐ（０）と表す。
【００７３】
まず、ｒ（０）は、ｘ（０）と各ノード（節点）に加わる力をベクトル形式で記述した節点力ベクトルｂとを用いて次式で計算される。
【００７４】
ｒ（０）  ＝  ｂ  −  Ａｘ（０）  …（式２６）
探索方向ベクトルｐの初期値ｐ（０）はｒ（０）と等しくする。
【００７５】
ステップＳ３０２では、反復回数ｔを１に初期化する。反復回数ｔは、上述の状態変数の求解工程に対する反復処理に対するものであり、以下のステップＳ３０３からＳ３０７の処理の反復回数が設定値を超えるか、残差ベクトルｒ（ｔ）のノルムの２乗が予め設定された値より小さくなったら、終了する。以下、反復回数ｔ回目の値をｘ（ｔ）のように書く。
【００７６】
ステップＳ３０３では状態ベクトルの更新係数α（ｔ）を計算する。
【００７７】
α（ｔ−１）  ＝  （ｒ（ｔ−１），Ａｐ（ｔ−１））／（Ａｐ（ｔ−１），Ａｐ（ｔ−１））  …（式２７）
ステップＳ３０４では、状態変数ベクトルｘを更新する。
【００７８】
ｘ（ｔ）  ＝  ｘ（ｔ−１）  ＋  α（ｔ−１）ｐ（ｔ−１）  …（式２８）
ステップＳ３０５では、残差ベクトルｒ（ｔ）を次式で計算する。
【００７９】
ｒ（ｔ）  ＝  ｒ（ｔ−１）  −α（ｔ−１）Ａｐ（ｔ−１）  …（式２９）
ステップＳ３０６では、探索方向ベクトルｐの更新係数βを計算する。
【００８０】
β（ｔ）  ＝  −  （Ａｒ（ｔ），Ａｐ（ｔ−１））／（Ａｐ（ｔ−１），Ａｐ（ｔ−１））  …（式３０）
ステップＳ３０７では、探索方向ベクトルｐを更新する。
【００８１】
ｐ（ｔ）  ＝  ｒ（ｔ）  ＋  β（ｔ−１）ｐ（ｔ−１）  …（式３１）
共役残差法及びその様々な派生法は、（参考文献３）、（参考文献４）、（参考文献５）、（参考文献６）、（参考文献７）等によって公知である。これらの方法が、係数行列Ａが特異となる場合においても、安定に近似解に収束することが証明されている。
【００８２】
（参考文献３）  森，杉原，室田：“線形計算，”岩波講座  応用数学，岩波書店（１９９４年）
（参考文献４）  速見：“特異な系に対するＧＣＲ（ｋ）法の収束性について，”京都大学数理解析研究所講究録１２６５（２００２年）
（参考文献５）  阿部，緒方，杉原，張，三井：“特異な係数行列をもつ連立一次方程式に対するＣＲ法の収束性，”日本応用数理学会論文誌，Ｖｏｌ．９，  Ｎｏ．１，  ｐｐ．１−１３  （１９９９年）
（参考文献６）  Ｋ．  Ｈａｙａｍｉ：  ”Ｏｎ  ｔｈｅ  Ｂｅｈａｖｉｏｕｒ  ｏｆ  ｔｈｅ  Ｃｏｎｊｕｇａｔｅ  Ｒｅｓｉｄｕａｌ  Ｍｅｔｈｏｄ  ｆｏｒ  Ｓｉｎｇｕｌａｒ  Ｓｙｓｔｅｍｓ，”  ＮＩＩ  Ｔｅｃｈｎｉｃａｌ  Ｒｅｐｏｒｔ，  ＮＩＩ−２００１−００２Ｅ  （２００２年）
（参考文献７）  Ｓ−Ｌ．  Ｚｈａｎｇ，  Ｙ．  Ｏｙａｎａｇｉ，  Ｍ．  Ｓｕｇｉｈａｒａ：  ”Ｎｅｃｅｓｓａｒｙ  ａｎｄ  ｓｕｆｆｉｃｉｅｎｔ  ｃｏｎｄｉｔｉｏｎｓ  ｆｏｒ  ｔｈｅ  ｃｏｎｖｅｒｇｅｎｃｅ  ｏｆ  Ｏｒｔｈｏｍｉｎ（ｋ）  ｏｎ  ｓｉｎｇｕｌａｒ  ａｎｄ  ｉｎｃｏｎｓｉｓｔｅｎｔ  ｓｙｓｔｅｍｓ，”  Ｎｕｍｅｒｉｓｃｈｅ  Ｍａｔｈｅｍａｒｉｋ，  ８７，  ｐｐ．３９１−４０５  （２０００年）（本実施形態の構造最適設計の具体例）
本例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に、本実施形態を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８３】
図５に示すように、構造部材の存在を可能とする設計領域は長方形４０２であり、有限要素法に従って、該領域を縦ｎｙ、横ｎｘに等間隔に分割する。分割された部分領域をセルと呼び、左下及び右上のセルがそれぞれ（１，１）及び（ｎｙ，ｎｘ）となるように番号付けを行う。
【００８４】
同様に格子点をノードと呼び、左下及び右上のノードが（１，１）及び（ｎｙ＋１，ｎｘ＋１）となるように番号付けを行う。
【００８５】
図中、４０１は支持部材、４０３は荷重ベクトルである。
【００８６】
セル（ｊ，ｋ）には特性関数値ｆ（ｊ，ｋ）が対応する。ここで、特性関数値とはセル（ｊ，ｋ）における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本実施形態における設計変数ベクトルｆの要素である。
ｆ  ＝  （ｆ（１，１），  ｆ（１，２），  …，  ｆ（ｎｙ，ｎｘ））Ｔ  …（式３２）
同様にノード（ｊ，ｋ）には横方向変位ｕ（ｊ，ｋ）と縦方向変位ｖ（ｊ，ｋ）が対応する。これらは任意の値を取る実数であり、本実施形態における状態変数ベクトルｘの要素である。
ｘ  ＝  （ｕ（１，１），  ｖ（１，１），  ｕ（１，２），  ｖ（１，２），…，ｕ（ｎｙ＋１，ｎｘ＋１），  ｖ（ｎｙ＋１，ｎｘ＋１））Ｔ…（式３３）
図６に、本具体例に置ける計算結果を示す。図中、黒い領域が構造部材が存在する領域である。本具体施例において、構造部材の設計領域のアスペクト比は縦／横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。図６に示した計算結果は、このような解析解とよく一致していることがわかる。
【００８７】
尚、本実施形態では、図１のように、バスライン２０７に各装置が接続した１つのシステムとして説明しているが、これら各要素をネットワークなどを介して接続して複数の装置のシステムで実現することも、あるいは汎用のコンピュータを複数接続して並列に多重処理あるいは分担処理により処理するようにしても良い。
【００８８】
又、本発明は、上記構造最適設計方法を実現するプログラム、及び該プログラムをコンピュータ読み出し可能に記憶する記憶媒体をも含むものである。
【００８９】
【発明の効果】
以上説明したように、本発明によれば、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量の軽減が実現できる。
【図面の簡単な説明】
【図１】本実施形態の構造最適設計装置の構成例を示す図である。
【図２】図１の１次記憶装置及び２次記憶装置の記憶構成例を示す図である。
【図３】本実施形態の処理手順例を示すフローチャートである。
【図４】図３のステップＳ１０５を更に詳細に示すの処理手順のフローチャートである。
【図５】本実施形態に適応する具体例の問題設定の説明図である。
【図６】図５に係る本実施形態による計算結果を示す図である。
【図１】
2004310375.tif 000007
【図２】
2004310375.tif 000008
【図３】
2004310375.tif 000009
【図４】
2004310375.tif 000010
【図５】
2004310375.tif 000011
【図６】
2004310375.tif 000012