(57)【特許請求の範囲】
【請求項１】
変数を記憶するための第１及び第２の記憶手段と、プログラムを記憶した第３の記憶手段と、プログラムを実行する実行手段とを備え、前記実行手段が、前記第３の記憶手段に記憶されたプログラムを実行することにより実現される、第１及び第２の求解手段と、設計変数更新手段と、第１及び第２の状態変数更新手段と、定数更新手段と、感度係数更新手段と、残差更新手段と、探索方向更新手段とを有する構造最適設計装置において、
前記第１の求解手段が、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解工程(a)と、前記第２の求解手段が、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解工程(b)とを有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法であって、
最適化すべき構造を縦横に分割して各部分を要素とし、各要素の頂点をノードとして、前記状態変数ベクトルを各ノードにおける変位、前記設計変数ベクトルを各要素における構造部材の存在率とし、
前記第１の求解工程(a)は、
(a-1)前記定数更新手段が、未定定数を更新する定数更新工程と、
(a-2)前記設計変数更新手段が、第１の記憶手段に記憶された設計変数ベクトルおよび状態変数ベクトルを読み出して、前記未定定数及び前記第１の評価汎関数の設計変数ベクトルに対する感度係数に基づいて該設計変数ベクトルを更新し、更新された設計変数ベクトルを該第１の記憶手段に記憶する設計変数更新工程と、
(a-3)前記第１の状態変数更新手段が、第２の記憶手段に記憶された設計変数ベクトルおよび状態変数ベクトルを読み出して状態変数ベクトルを更新し、更新された状態変数ベクトルを該第２の記憶手段に記憶する第１の状態変数更新工程と、
(a-4)前記感度係数更新手段が、前記設計変数ベクトル及び前記設計変数ベクトルに基づいて前記感度係数を更新する感度係数更新工程とを、所定の終了条件を満たすまで繰り返し、
前記第１の状態変数更新工程(a-3)が前記第２の求解工程(b)を有しており、
該第２の求解工程(b)では、
(b-1)前記第２の評価汎関数が、全体剛性行列を作用させた状態変数ベクトルと節点力ベクトルとの差を表す残差ベクトルのノルムにより構成されており、
(b-2-1)前記第２の状態変数更新手段が、探索方向ベクトルと前記残差ベクトルとに基づいて状態変数ベクトルを更新する第２の状態変数更新工程と、
(b-2-2)前記残差更新手段が、前記探索方向ベクトルと前記全体剛性行列とに基づいて前記残差ベクトルを更新する残差更新工程と、
(b-2-3)前記探索方向更新手段が、前記残差ベクトルと前記全体剛性行列とに基づいて前記探索方向ベクトルを更新する探索方向更新工程とを、前記残差ベクトルのノルムの２乗が設定値より小さくなるか所定の反復回数を越えるまで繰り返すことを特徴とする構造最適設計方法。
【請求項２】
前記第１の求解工程は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかにより実現されることを特徴とする請求項１に記載の構造最適設計方法。
【請求項３】
前記第２の求解工程は、共役残差法、GCR法、GCR(k)法、Orthomin(k)法、GMRES(k)法、及びこれらの派生法を含む方法のいずれかにより実現されることを特徴とする請求項１又は２に記載の構造最適設計方法。
【請求項４】
請求項１乃至３のいずれか１項に記載の構造最適設計方法の工程をコンピュータに実行させるためのプログラム。
【請求項５】
請求項４に記載のプログラムを記憶したコンピュータ読み取り可能な記憶媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は設計パラメータの最適化問題の解法である構造物の最適設計、特に構造部材のトポロジーを含む形状最適化のための構造最適設計方法及びプログラム、記憶媒体に関するものである。
【０００２】
【従来の技術】
構造トポロジーの最適設計とは、与えられた条件の下で、構造部材の最適なトポロジー及び形状寸法を決定する問題である。以下、構造部材のトポロジー及び形状寸法を設計変関数といい、上記決定問題を設計変関数の最適化問題という。変関数という理由は、トポロジー及び形状寸法が３次元空間の関数になっているからである。設計変関数の最適化問題では、各設計変関数の値に対して、状態変関数の最適化問題を解かなければならない。この意味から、構造トポロジーの最適設計は、内側に状態変関数の最適化問題を持ち、外側に設計変関数の最適化問題を持つ２重構造最適化問題と捉えることが出来る。
【０００３】
内側の状態変関数最適化問題では、従来技術の蓄積により、空間を有限個の要素に分割するという考え方が採用される。特に構造部材の歪エネルギーを評価汎関数としている問題では、その解析手法として有限要素法が一般的である。有限要素法の解法としては一次方程式に対する直接法が採用されている。
【０００４】
一方、設計変関数の最適化問題に対しては、大別すると、以下に示す３種類の方法が提供されている（非特許文献１、或いは、非特許文献２参照）。
１．Evolutionary法(以下、Ｅ法)
２．Homogenization法(以下、Ｈ法)
３．Material distribution法(以下、ＭＤ法)、或いは密度(Density)法(以下、Ｄ法)
Ｅ法では、空間を分割することによって得られる部分空間のそれぞれをセルと称し、その生成と消滅を適当な規則によって繰り返す。構造部材は、最終的に存在しているセルの集合として与えられる。セルが存在するか否かという２つの状態のみを許すことにより、明確な構造部材が得られる。また、評価汎関数の微分情報を用いないので、局所最適解にトラップされないことから、評価汎関数が多峰性の場合に有効である。特許文献１では、Ｅ法の一種である遺伝的探索法を用いた、骨組構造部材の最適化設計装置が提供されている。この最適化設計装置では、蓄積されたノウハウに基づく試行計算を必要とし、従って膨大な設計変数が存在する実設計問題に対応することが出来なかった従来技術の課題を、以下の方法によって解決している。即ち、骨組部材断面寸法などの離散設計変数データの近似式を使用する近似最適化計算装置と、該設計変数データを使用する詳細最適化計算装置を設け、これら２つの計算装置を結合して骨組構造の最適設計装置を構成している。
【０００５】
Ｈ法は、分割された各部分領域に位置する構造要素に、更に微細な構造を仮定し、連続値を取る設計変関数を新たに導入することによって、感度解析の採用を可能にしている。ここで感度解析とは、設計変関数に関する評価汎関数の微分情報を利用した解析手法のことであり、感度解析が可能になれば、勾配法のような反復解法を用いることが出来、Ｅ法のような総当り的手法に比べて、少なくとも局所最適解の探索に係る計算時間は大幅に短縮される（非特許文献３参照）。
【０００６】
ＭＤ法或いはＤ法は、構造部材の存在率を示す０から１の範囲の実数を各要素に割当てることによって、構造部材のトポロジーや形状寸法変化を表現する方法である。構造部材が存在するか否かという離散的な情報を存在率という連続値に置き換えることによって感度解析を可能にしたという意味でＨ法と同様のものであるが、Ｈ法よりパラメータ数が少ない分、モデル化が容易であり適用範囲も広い。
【０００７】
非特許文献４には、Ｄ法による構造物の位相最適化手法が開示されている。本手法の特徴は以下のとおりである。ボクセル有限要素法（空間を等間隔に分割）を用いているため、あらゆる要素に対する要素剛性マトリクスが同じである。従って、要素剛性マトリクスを予め１度だけ計算しておけば、以後の計算に利用できる。更に、要素が規則正しく配置されているため、各要素の節点番号情報を記憶する必要がない。大規模連立一次方程式を解くために、前処理付共役勾配法とElement-by-Element法を組み合わせて用いたことにより、全体剛性マトリクスを組み立てることなく解が求められるので、必要とするメモリ容量が少なくて済む。
【０００８】
均質化法では、１要素に対して６個の設計変数（３次元の場合）が必要になる。更に設計変数が変化するたびに要素剛性マトリクスを再計算しなければならない。一方、構造部材の存在率を密度比で表現する密度法を採用することによって、１要素に対して１つの設計変数でよい。また設計変数の変化が要素剛性マトリクスに影響を与えない。
【０００９】
【特許文献１】
特開平11-314631号
【非特許文献１】
S. Bulman, J. Sienz, E. Hinton: "Comparisons between algorithms for structural topology optimization using a series of benchmark studies," Computers and Structures, 79, pp.1203-1218 (2001).
【非特許文献２】
Y-L. Hsu, M-S. Hsu, C-T. Chen: "Interpreting results from topology optimization using density contours," Computers and Structures, 79, pp.1049-1058 (2001).
【非特許文献３】
山川宏：“最適化デザイン,”計算力学とCAEシリーズ９,培風館(１９９６)
【非特許文献４】
藤井,鈴木,大坪：“ボクセル有限要素法を用いた構造物の位相最適化,”Transactions of JSCES, Paper No.20000010 (2000).
【００１０】
【発明が解決しようとする課題】
しかしながら、上記従来技術には以下に述べるような問題点があった。
【００１１】
一般に構造最適化問題は、設計変数ベクトルの最適化問題の各反復ステップ毎に、状態変数ベクトルの最適化問題を含む２重の最適化問題として定式化される。設計変数ベクトルの最適化問題を外側の最適化、状態変数ベクトルの最適化問題を内側の最適化と呼ぶと、内側の最適化は、設計変数ベクトルをパラメータとし、即ち設計変数ベクトルを固定し、状態変数ベクトルを求める問題である。これは通常、構造解析と呼ばれ、有限要素法によって線形方程式の解法を用いて解くことが出来る。
【００１２】
ところが、構造が変化し、ある領域の構造部材が存在しなくなる、例えば部材に穴が空いてしまうと、対応する要素の設計変数が０となり、結果的に要素のヤング率が０となる。すると、上記線形方程式の係数行列がフルランクでなくなり、逆行列が計算できないという理由で、直接法では解けなくなる。
【００１３】
そのため従来技術では、設計変数ベクトルが０となる場合に以下のように対処してきた。
【００１４】
係数行列がフルランクになるように、即ちランク数と方程式の数が等しくなるように、方程式を構成し直す。
【００１５】
設計変数ベクトルの値を正確に０とはせずに０に近い微小な値で置き換える。
【００１６】
しかるに、上記前者の方法は、設計変数ベクトルが更新される毎に方程式の更新を行う可能性が生じ、計算時間がかかる。
【００１７】
また、上記後者の方法、即ち要素の設計変数を微小な値に設定するということは、物理的には薄い膜、或いは弱い部材が存在することに相当する。つまり従来法では、物質がない部分を正確に表現していないので、得られた計算結果の精度に疑問が残る。
【００１８】
本発明は、上記従来の問題点に鑑み、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量を軽減する構造最適設計方法及びプログラム、記憶媒体を提供する。
【００１９】
【課題を解決するための手段】
上記課題を解決するために、本発明の構造最適設計方法は、変数を記憶するための第１及び第２の記憶手段と、プログラムを記憶した第３の記憶手段と、プログラムを実行する実行手段とを備え、前記実行手段が、前記第３の記憶手段に記憶されたプログラムを実行することにより実現される、第１及び第２の求解手段と、設計変数更新手段と、第１及び第２の状態変数更新手段と、定数更新手段と、感度係数更新手段と、残差更新手段と、探索方向更新手段とを有する構造最適設計装置において、前記第１の求解手段が、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解工程(a)と、前記第２の求解手段が、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解工程(b)とを有する２重最適化問題として定式化される構造最適設計問題の解を求める構造最適設計方法であって、最適化すべき構造を縦横に分割して各部分を要素とし、各要素の頂点をノードとして、前記状態変数ベクトルを各ノードにおける変位、前記設計変数ベクトルを各要素における構造部材の存在率とし、前記第１の求解工程(a)は、(a-1)前記定数更新手段が、未定定数を更新する定数更新工程と、(a-2)前記設計変数更新手段が、第１の記憶手段に記憶された設計変数ベクトルおよび状態変数ベクトルを読み出して、前記未定定数及び前記第１の評価汎関数の設計変数ベクトルに対する感度係数に基づいて該設計変数ベクトルを更新し、更新された設計変数ベクトルを該第１の記憶手段に記憶する設計変数更新工程と、(a-3)前記第１の状態変数更新手段が、第２の記憶手段に記憶された設計変数ベクトルおよび状態変数ベクトルを読み出して状態変数ベクトルを更新し、更新された状態変数ベクトルを該第２の記憶手段に記憶する第１の状態変数更新工程と、(a-4)前記感度係数更新手段が、前記設計変数ベクトル及び前記設計変数ベクトルに基づいて前記感度係数を更新する感度係数更新工程とを、所定の終了条件を満たすまで繰り返し、前記第１の状態変数更新工程(a-3)が前記第２の求解工程(b)を有しており、該第２の求解工程(b)では、(b-1)前記第２の評価汎関数が、全体剛性行列を作用させた状態変数ベクトルと節点力ベクトルとの差を表す残差ベクトルのノルムにより構成されており、(b-2-1)前記第２の状態変数更新手段が、探索方向ベクトルと前記残差ベクトルとに基づいて状態変数ベクトルを更新する第２の状態変数更新工程と、(b-2-2)前記残差更新手段が、前記探索方向ベクトルと前記全体剛性行列とに基づいて前記残差ベクトルを更新する残差更新工程と、(b-2-3)前記探索方向更新手段が、前記残差ベクトルと前記全体剛性行列とに基づいて前記探索方向ベクトルを更新する探索方向更新工程とを、前記残差ベクトルのノルムの２乗が設定値より小さくなるか所定の反復回数を越えるまで繰り返すことを特徴とするを特徴とする。
【００２０】
【発明の実施の形態】
以下、本発明の構造最適設計法を実現する構造最適設計装置の構成例及び動作例を説明する。
【００２１】
＜本実施形態の構造最適設計装置の構成例＞
図１に、本実施形態の構造最適設計装置の構成例を示す。
【００２２】
図中、２０１は演算・制御用のＣＰＵ、２０２は表示装置、２０３は入力装置、２０４は１次記憶装置、２０５は２次記憶装置、２０６は通信装置、２０７はバスラインである。本実施形態における構成要素は、プログラムとして予め２次記憶装置２０５に格納されており、入力装置２０３或いは通信装置２０６等からのコマンド入力によって１次記憶装置にロードされ、ＣＰＵ２０１により実行されるものである。
【００２３】
図２に、本実施形態の構造最適設計法を実現するために、１次記憶装置２０４及び２次記憶装置２０５に記憶されるデータ及びプログラムの記憶例を示す。
【００２４】
１次記憶装置２０４は、ＲＡＭやＲＯＭからなり、データ記憶領域２１４とプログラム記憶領域２２４を有している。図２には、変更されずに常駐するＯＳやＢＩＯＳのプログラム、変更されないＯＳやＢＩＯＳ用のパラメータなどは図示されていない。データ記憶領域２１４には、本実施形態の構造最適設計を実現するための各種のデータが一時記憶される。この中で再度利用されるデータ／反復して使用されるデータは、一時記憶されることが必須である。保存・保持の必要が無いデータは、一時記憶は必須ではない。
【００２５】
例えば、データ記憶領域２１４には、ｘ（状態変数ベクトル）２１４ａ、ｆ（設計変数ベクトル）２１４ｂ、Ｌ1（状態変数ベクトルの評価関数）２１４ｃ、Ｌ2（設計変数ベクトルの評価関数）２１４ｄ、Ｂ1（状態変数のための境界条件）２１４ｅ、Ｗ0（設計変数のための重量条件）２１４ｆ、Ｘ(0)やｆ(0)（初期値）２１４ｇ、λ（ラグランジェ未定定数）２１８ｈ、Ａ（要素合成行列）２１４ｉ、ｋ（設計変数変更の反復回数）２１４ｊ、Ｃ（感度係数）２１４ｋ、ｒ（残差ベクトル）１２４ｍ、ｐ（探索方向ベクトル）２１５ｎ、ｒ(0)やｐ(0)（初期値）２１４ｐ、ｔ（状態変数変更の反復回数）２１４ｑ、ｂ（節点力ベクトル）２１４ｒ、その他のデータ２１４ｓなどが含まれる。一方、プログラム記憶領域２２４には、２次記憶装置２０５或は通信装置２０６からロードされたプログラムが記憶されて、ＣＰＵ２０１により実行される。
【００２６】
２次記憶装置２０５は、例えば、フロッピーディスクやＣＤなどの外部大容量メモリであり、着脱可能であればなお好ましい。２次記憶装置２０５も、データ記憶領域２１５とプログラム記憶領域２２５を有している。
【００２７】
データ記憶領域２１５には、詳細には示さないが、各種の初期値２１５ａや各種の条件２１５ｂ、あるいは２重最適値問題のデータベース２１５ｃが記憶されている。データベース２１５ｃには、標準化されたマトリクスなどが記憶されていて使用されてもよい。
【００２８】
プログラム記憶領域２２５には、本実施形態で設計変数ベクトルにおける求解法として使用される、最適性規準法モジュール２２５ａ、逐次凸関数近似法モジュール２２５ｂ、逐次線形計画法モジュール２２５ｃ、その他の方法のモジュール２２５ｄが記憶され、状態変数ベクトルにおける求解法として使用される、共役残差法モジュール２２５ｅ、ＧＣＲ法モジュール２２５ｆ、ＧＣＲ（ｋ）法２２５ｇ、Orthomin（ｋ）法２２５ｈ、ＧＭＲＥＧ（ｋ）法２２５ｉ、その他の方法のモジュール２２５ｊが記憶されている。このプログラム記憶領域２２５に記憶されたプログラムから選択された方法に対応するモジュールが１次記憶装置２０４のプログラム記憶領域２２４にロードされて、ＣＰＵ２０１により実行され本実施形態の構造最適設計法が実現する。
【００２９】
すなわち、本構造最適設計装置は、状態変数ベクトルと設計変数ベクトルに対する第１の評価汎関数の最適化問題を解く第１の求解手段と、該状態変数ベクトルと該設計変数ベクトルに対する第２の評価汎関数の最適化問題を解く第２の求解手段とを有し、２重最適化問題として定式化される構造最適設計問題を解く構造最適設計装置であって、前記状態変数ベクトルが各ノードにおける変位、前記設計変数ベクトルが各要素における構造部材の存在率であり、前記第２の求解手段の第２の評価汎関数は残差ベクトルのノルムより構成されており、前記状態変数ベクトルが前記第２の求解工程が開始される時点で初期化されない。前記第１の求解手段は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかを実行する。前記第２の求解手段は、共役残差法、GCR法、GCR(k)法、Orthomin(k)法、GMRES(k)法、及びこれらの派生法を含む方法のいずれかを実行する。
【００３０】
＜本実施形態の構造最適化問題の定式化＞
以下の説明のために、構造最適化問題の定式化を行う。
【００３１】
有限要素法による定式化により変関数は有限次元ベクトルで表現されるものとすると、変関数の評価汎関数は、変数ベクトルの評価関数となる。以下、有限次元ベクトルで表現されているとして記述する。
【００３２】
状態変数ベクトルｘ、設計変数ベクトルｆをそれぞれ列ベクトルとして以下のように書く。
【００３３】
ｘ = (ｘ1,ｘ2,…,ｘm)T  …（式１）
ｆ = (ｆ1,ｆ2,…,ｆn)T  …（式２）
ここで、Ｔは転置を表す。ｘはｍ次元ベクトル、ｆはｎ次元ベクトルである。構造最適化問題において、状態変数ベクトルは各ノードの変位、設計変数ベクトルは各要素内における構造物の存在率より構成される。
【００３４】
例えば、２次元面内での構造最適化問題において、設計領域を縦ｎy、横ｎxに分割すると、
要素数ｎは(ｎx×ｎy)、ノード数は(ｎx＋１)×(ｎy＋１)である。
【００３５】
状態変数ベクトルは各ノードの縦方向及び横方向変位の組なので、ｘの次元ｍは２×(ｎx＋１)×(ｎy＋１)である。
【００３６】
ｘの境界条件は通常、変位拘束条件として与えられる。それをＢ1と記述する。
Ｂ1(ｘ)=0  …（式３）
状態変数ベクトルｘ及び設計変数ベクトルｆに関する評価関数をそれぞれＬ1、Ｌ2とすると、次式のように定義される。
Ｌ1 := Ｌ1(ｘ,ｆ) = (ｂ - Ａｘ)T(ｂ - Ａｘ)  …（式４）
Ｌ2 := Ｌ2(ｆ,ｘ) = (1/2)ｘTＡｘ  …（式５）
ここで、Ａは全体剛性行列、ｂは節点力ベクトルである。ただし、ｂは予め与えられるものであり、Ａは設計変数ベクトルｆの関数である。
【００３７】
全体剛性行列は要素ｊに対する要素剛性行列Ａjの重ね合わせで構成することが出来る。平面歪問題の場合、設計変数ベクトル、即ち構造要素の存在率ｆjを考慮した要素剛性行列を以下に示す。ただし、２次元有限要素として４ノード要素を採用し、要素変位ベクトルｘjを次のように構成する。
ｘj = (ｕ1, ｖ1, ｕ2, ｖ2,… , ｕ4, ｖ4)  …（式６）
ここで、ｕk、ｖkはノードｋの水平方向、及び垂直方向の変位である。
【００３８】
（式６）の要素変位ベクトルに対応する要素剛性行列は、次式のように与えられる。
【００３９】
【数１】
0004095484.tif 000002
【００４０】
上式中
ａi,j = (λ+2μ)＜∂xei,∂xej＞  + μ＜∂yei,∂yej＞  …（式８）
ｂi,j =       λ＜∂xei,∂yej＞  + μ＜∂xej,∂yei＞  …（式９）
ｃi,j =       λ＜∂xej,∂xei＞  + μ＜∂xei,∂yej＞  …（式１０）
ｄi,j = (λ+2μ)＜∂yei,∂yej＞  + μ＜∂xei,∂xej＞  …（式１１）
ここで、ｅjはノードｊに対する基底ベクトル、＜  ,  ＞はベクトルの内積を表す。∂xはｘに関する偏微分、∂yはｙに関する偏微分である。ξは前述の（非特許文献４）に従い、２とする。また、λ、μはラメ定数と呼ばれる材料定数であり、ヤング率Ｅ及びポアソン比νより次式で計算される。
【００４１】
λ = Eν/(1+ν)(1-2ν)  …（式１２）
μ = E/{2(1+ν)}  …（式１３）
（式７）より容易に分かるように、ある要素の存在率が０となると要素剛性行列Ａjのすべての成分が０となる。全体剛性行列は要素剛性行列の重ね合わせで構成される。従って、ノードjを含む要素の設計変数ｆjが０になると、全体剛性行列のノードｊに対応する行及び列の全ての成分が０になる。このように、構造最適化問題において全体剛性行列は一般に特異な行列となる。
【００４２】
設計変数に関する拘束条件は、構造最適設計の場合、総重量の上限値Ｗ0とする。
【００４３】
【数２】
0004095484.tif 000003
【００４４】
また、設計変数には値域に対する拘束条件が以下のように定められている。
【００４５】
0 ≦ ｆj ≦ 1    j = 1, … , n    …（式１５）
以上の表記より、最適設計問題は、拘束条件（式１４）及び（式１５）の下で、（式５）を極小化する問題として定式化される。
【００４６】
また、状態変数ｘについては、拘束条件（式３）の下で（式４）を極小化する解として得られる。
【００４７】
＜本実施形態の構造最適設計装置の動作例＞
上記定式化に基づき、本実施形態の処理手順例を説明する。
【００４８】
図３には、本実施形態の流れ図が示されている。
【００４９】
図中、ステップＳ１０１はシミュレーション対象となる系の諸元を読み込み、変数の初期化を行う処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値としてｘ(0)、f(0)、境界条件Ｂ、評価関数Ｌ1、Ｌ2が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。感度係数の初期値Ｃj(0)は次式で計算する。
【００５０】
Ｃj(0) = - (ｘj(0))T(ξ(ｆj(0))ξ-1Ａj)ｘj(0)    …（式１６）
ステップＳ１０２では、反復回数ｋを１に初期化する。反復回数ｋは、上述の設計変数の求解工程に対する反復処理に対するものであり、予め定められた終了条件（評価関数Ｌ1，Ｌ2の条件）を満たしたら終了する。
【００５１】
ステップＳ１０３では、ラグランジェ未定定数λの更新を行う。ステップＳ１０４では、設計変数ベクトルｆの更新を行う。ステップＳ１０５では、状態変数ベクトルｘの更新を行う。ステップＳ１０６では、感度係数の計算を行う。
【００５２】
上記ステップＳ１０３及びＳ１０４は、拘束条件付最適化問題の解法として採用する方法によって異なる。例えば、この設計変数ベクトルｆ(k)を算出する求解工程は、逐次線形計画法、最適性規準法、逐次凸関数近似法のいずれかにより実現される。以下、解法として最適性規準法を用いる場合と、逐次凸関数近似法を用いる場合について詳述する。
【００５３】
（最適性規準法）
まず、最適性規準法を用いる場合について説明する。構造最適化問題に対する最適性規準法は（非特許文献４）及び（非特許文献５）等によって公知である。
【００５４】
最適性規準法では、（式４）及び（式５）をラグランジェ未定定数λを用いて、１つの式にまとめる。
【００５５】
Ｌ(f,λ) = Ｌ2(f,ｘ) - λ(Ｗ(f) - Ｗ0)  …（式１７）
ステップＳ１０３で実行されるλの更新は、次式を用いて行われる。
【００５６】
λ(k+1) = min[0,[(1/Ｗ0)Ｗ(f(k))] αλ(k)]  …（式１８）
ただし、上付きのｋ，ｋ＋１は反復回数を表す。また、αはべき乗係数で、通常0.85程度の値を設定する。
【００５７】
ステップＳ１０４で実行される設計変数の更新は次式による。
【００５８】
0004095484.tif 000004
ただし、ζは設計変数更新の変動幅の制約値で、0.3程度の値を取る。また、ｓj(k)は次式で与えられる。
【００５９】
ｓj(k) = [(1/λ(k-1))Ｃj(k-1)]αfj(k-1)  …（式２０）
ここで、Ｃj(k)は第２の評価関数Ｌ2の、設計変数ｆj(k)に関する感度係数と呼ばれる量であり、ステップＳ１０６の処理で計算されるものである。
【００６０】
更に（式１９）で計算された設計変数ｆj(k)は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、例えば１０-3程度の値である。
【００６１】
ステップＳ１０５では、状態変数ベクトルｘ(k)を更新する。
【００６２】
状態変数ベクトルｘ(k)は、（式４）の評価関数の極小化の過程で得られる。この解法は共役残差法と呼ばれており、後述する。
【００６３】
ステップＳ１０６では、感度係数Ｃj(k)を次式により計算する。
【００６４】
Ｃj(k) = - (ｘj(k))T(ξ(ｆj(k))ξ-1Ａj)ｘj(k)  …（式２１）
（逐次凸関数近似法）
次に、逐次凸関数近似法を用いた場合の、ステップＳ１０３及びＳ１０４の処理について詳述する。逐次凸関数近似法では設計変数ベクトルの各成分に対してスケーリングを行うことによって、数値計算誤差を避けることができるが、以下では説明を簡単にするため、スケーリングについては記述していない。逐次凸関数近似法の詳細は、（参考文献１）及び（参考文献２）等によって公知である。
【００６５】
（参考文献１）藤井：“パソコンで解く構造デザイン”, 丸善(2002年4月)
（参考文献２）C. Fleury : "CONLIN, an efficient dual optimizer based on convex approximation concepts," Structural Optimization, 1, pp.81-89 (1989)
ステップＳ１０３では、ラグランジェ未定定数を次式により更新する。
【００６６】
λ(k) = λ(k-1) + Δλ(k)  …（式２２）
ただし、初期値λ(0)は例えば0.01程度の値を設定する。
【００６７】
Δλ(k)は次の方程式の解として与えられる。
【００６８】
Δλ(k) = ２×( 1 - Ｗ0/Ｗ(f(k-1)) )λ(k-1)    …（式２３）
ステップＳ１０４では、設計変数ベクトルの更新を行う。
【００６９】
0004095484.tif 000005
更に、（式２４）で計算された設計変数ｆj(k)は、予め設定された値より小さくなると強制的に０にしてもよい。この予め設定された値とは、１０-3程度の値である。
【００７０】
ステップＳ１０５及びＳ１０６の処理は前述のとおりである。
【００７１】
（共役残差法）
ステップＳ１０５において状態変数ベクトルｘ(k)を算出する求解工程は、共役残差法、GCR法、GCR(k)法、Orthomin(k)法、GMRES(k)法、及びこれらの派生法を含む方法のいずれかにより実現される。以下、それらを代表して共役残差法を用いた処理を図４を用いて説明する。
【００７２】
ステップＳ３０１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘ(k)の初期値、設計変数ベクトルｆ(k)の値、変位拘束条件Ｂ、評価関数Ｌ1が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保して、値を設定する。また、状態変数ベクトルｘ、残差ベクトルｒ、探索方向ベクトルｐを初期化する。それぞれの初期値をｘ(0)、ｒ(0)、ｐ(0)と表す。
【００７３】
まず、ｒ(0)は、ｘ(0)と各ノード（節点）に加わる力をベクトル形式で記述した節点力ベクトルｂとを用いて次式で計算される。
【００７４】
ｒ(0) = ｂ - Ａｘ(0)  …（式２６）
探索方向ベクトルｐの初期値ｐ(0)はｒ(0)と等しくする。
【００７５】
ステップＳ３０２では、反復回数ｔを１に初期化する。反復回数ｔは、上述の状態変数の求解工程に対する反復処理に対するものであり、以下のステップＳ３０３からＳ３０７の処理の反復回数が設定値を超えるか、残差ベクトルｒ(t)のノルムの２乗が予め設定された値より小さくなったら、終了する。以下、反復回数ｔ回目の値をｘ(t)のように書く。
【００７６】
ステップＳ３０３では状態ベクトルの更新係数α(t)を計算する。
【００７７】
α(t-1) = (r(t-1),Ａp(t-1))/(Ａp(t-1),Ａp(t-1))  …（式２７）
ステップＳ３０４では、状態変数ベクトルｘを更新する。
【００７８】
ｘ(t) = ｘ(t-1) + α(t-1)ｐ(t-1)  …（式２８）
ステップＳ３０５では、残差ベクトルｒ(t)を次式で計算する。
【００７９】
ｒ(t) = ｒ(t-1) -α(t-1)Ａｐ(t-1)  …（式２９）
ステップＳ３０６では、探索方向ベクトルｐの更新係数βを計算する。
【００８０】
β(t) = - (Ａｒ(t),Ａｐ(t-1))/(Ａｐ(t-1),Ａｐ(t-1))  …（式３０）
ステップＳ３０７では、探索方向ベクトルｐを更新する。
【００８１】
ｐ(t) = ｒ(t) + β(t-1)ｐ(t-1)  …（式３１）
共役残差法及びその様々な派生法は、（参考文献３）、（参考文献４）、（参考文献５）、（参考文献６）、（参考文献７）等によって公知である。これらの方法が、係数行列Ａが特異となる場合においても、安定に近似解に収束することが証明されている。
【００８２】
（参考文献３）  森，杉原，室田：“線形計算，”岩波講座  応用数学，岩波書店(１９９４年)
（参考文献４）  速見：“特異な系に対するGCR(k)法の収束性について，”京都大学数理解析研究所講究録１２６５（２００２年）
（参考文献５）  阿部，緒方，杉原，張，三井：“特異な係数行列をもつ連立一次方程式に対するCR法の収束性，”日本応用数理学会論文誌，Vol.9, No.1, pp.1-13 (1999年)
（参考文献６）  K. Hayami: "On the Behaviour of the Conjugate Residual Method for Singular Systems," NII Technical Report, NII-2001-002E (2002年)
（参考文献７）  S-L. Zhang, Y. Oyanagi, M. Sugihara: "Necessary and sufficient conditions for the convergence of Orthomin(k) on singular and inconsistent systems," Numerische Mathemarik, 87, pp.391-405 (2000年)
（本実施形態の構造最適設計の具体例）
本例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に、本実施形態を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８３】
図５に示すように、構造部材の存在を可能とする設計領域は長方形４０２であり、有限要素法に従って、該領域を縦ｎy、横ｎxに等間隔に分割する。分割された部分領域をセルと呼び、左下及び右上のセルがそれぞれ(１,１)及び(ｎy,ｎx)となるように番号付けを行う。
【００８４】
同様に格子点をノードと呼び、左下及び右上のノードが(１,１)及び(ｎy＋１,ｎx＋１)となるように番号付けを行う。
【００８５】
図中、４０１は支持部材、４０３は荷重ベクトルである。
【００８６】
セル(j,k)には特性関数値ｆ(j,k)が対応する。ここで、特性関数値とはセル(j,k)における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本実施形態における設計変数ベクトルｆの要素である。
ｆ = (f(1,1), f(1,2), …, f(ny,nｘ))T  …（式３２）
同様にノード(j,k)には横方向変位ｕ(j,k)と縦方向変位ｖ(j,k)が対応する。これらは任意の値を取る実数であり、本実施形態における状態変数ベクトルｘの要素である。
ｘ = (u(1,1), v(1,1), u(1,2), v(1,2),…,u(ny+1,nx+1), v(ny+1,nx+1))T…（式３３）
図６に、本具体例に置ける計算結果を示す。図中、黒い領域が構造部材が存在する領域である。本具体施例において、構造部材の設計領域のアスペクト比は縦/横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。図６に示した計算結果は、このような解析解とよく一致していることがわかる。
【００８７】
尚、本実施形態では、図１のように、バスライン２０７に各装置が接続した１つのシステムとして説明しているが、これら各要素をネットワークなどを介して接続して複数の装置のシステムで実現することも、あるいは汎用のコンピュータを複数接続して並列に多重処理あるいは分担処理により処理するようにしても良い。
【００８８】
又、本発明は、上記構造最適設計方法を実現するプログラム、及び該プログラムをコンピュータ読み出し可能に記憶する記憶媒体をも含むものである。
【００８９】
【発明の効果】
以上説明したように、本発明によれば、全体剛性行列が特異になったときでも、特別な処理を施すこと無しに計算を実行することができ、これによってプログラムが簡単になり、更に計算量の軽減が実現できる。
【図面の簡単な説明】
【図１】本実施形態の構造最適設計装置の構成例を示す図である。
【図２】図１の１次記憶装置及び２次記憶装置の記憶構成例を示す図である。
【図３】本実施形態の処理手順例を示すフローチャートである。
【図４】図３のステップＳ１０５を更に詳細に示すの処理手順のフローチャートである。
【図５】本実施形態に適応する具体例の問題設定の説明図である。
【図６】図５に係る本実施形態による計算結果を示す図である。
【図１】
0004095484.tif 000006
【図２】
0004095484.tif 000007
【図３】
0004095484.tif 000008
【図４】
0004095484.tif 000009
【図５】
0004095484.tif 000010
【図６】
0004095484.tif 000011