JP4101048B2 - Optimal design method - Google Patents

Description

【００２１】
ここでｃ_ｊ（ｊ＝０，１，・・・，ｎ）は定数である。
【００２２】
また（式８）の不等式拘束条件を次式のように、値の範囲に対するものとする：
（式１２）Ｒ_ｊ ⁻（ｆ）＝ｆ_ｔ（ｊ）−ｆ_ｊｍｉｎ≧０
Ｒ_ｊ ^＋（ｆ）＝ｆ_ｊｍａｘ−ｆ_ｔ（ｊ）≧０
ただし、ｊ＝１，２，…，ｎ。
【００２３】
従って、不等式拘束条件の数は設計変数の数の２倍となる。
【００２４】
設計変数ベクトルと同次元のフラグベクトルａ^＋、およびａ⁻を以下のように定義する：
ａ^＋＝０， ｉｆ ０≧Ｒ_ｊ ^＋（ｆ）
ａ^＋＝１， ｏｔｈｅｒｗｉｓｅ
および
ａ⁻＝０， ｉｆ ０≧Ｒ_ｊ ⁻（ｆ）
ａ⁻＝１， ｏｔｈｅｒｗｉｓｅ
また、ａ^＋およびａ⁻が共に１であるｊの集合を∧_１、それ以外のインデクスの集合を∧_２とする。
【００２５】
図１には、第二の求解工程のフローチャートが示されている。図中、ステップＳ１０１は、シミュレーション対象となる系の諸元を読み込む処理である。諸元の読み込みは、入力装置２０３或いは通信装置２０６からの入力データで行っても良いし、予め２次記憶装置２０５にファイルとして格納しておいたデータを読み出して利用することもできる。系の諸元にはｘ、ｆの初期値、境界条件Ｂ_１、Ｂ_２、評価関数Ｌ_１、Ｌ_２、拘束条件Ｑ_ｊ、Ｒ_ｋが含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保し、値を設定する。
【００２６】
ステップＳ１０２では、反復回数ｔを１に初期化する。反復回数ｔは、上述の第２の求解工程における第１の反復処理を行なった回数であり、以下のステップＳ１０３から１１１の処理の反復回数が設定値Ｔに達したら終了する。
【００２７】
ステップＳ１０３では勾配ベクトル計算工程を実行する。同工程では、設計変数ベクトルに関する第２の評価関数のグラジェントベクトルｇ_ｔを計算する。ｇ_ｔは次式で計算される：
（式１３）ｇ_ｔ＝（∂Ｌ_２／∂ｆ_１，∂Ｌ_２／∂ｆ_２，…，∂Ｌ_２／∂ｆ_ｎ）^Ｔ
∂Ｌ_２／∂ｆ_ｊは解析的に解ける場合とそうでない場合がある。解析解がある場合は、その計算工程を利用して（式１０）を得る。そうでない場合は、自動微分という技術を利用して計算することが出来る。自動微分技術は、以下の〔文献６〕等によって公知である。
【００２８】
〔文献６〕久保田光一，伊理正夫：“アルゴリズムの自動微分と応用，”現代非線形科学シリーズ，コロナ社（１９９８年）
ステップＳ１０４では、（式１３）で得られたｇ_ｔに対して、第一のベクトル修正工程を実行する。同工程は図６を用いて後述する。
【００２９】
ステップＳ１０５では、収束判定工程を実行する。同工程では、ｇ_ｔのノルムの２乗を計算し、その値が０或いは０に近い非常に小さい値ε未満になると処理を終了する。ノルムの２乗は次式で計算する：
（式１４）‖ｇ_ｔ‖^２＝ｇ_ｔ ^Ｔｇ_ｔ
ステップＳ１０６では、第一の係数計算工程を実行する。第一の係数計算工程では、（式１５）により係数βを計算する：
（式１５）β＝‖ｇ_ｔ‖^２／‖ｇ_ｔ−１‖^２
ただしｔ＝１のときはβ＝１とする。
【００３０】
ステップＳ１０７では探索ベクトル計算工程を実行する。同工程では（式１６）を用いて、反復回数ｔ回目の探索ベクトルｐ_ｔを計算する：
（式１６）ｐ_ｔ←−ｇ_ｔ−１＋βｐ_ｔ−１
ステップＳ１０８では、（式１６）より得られたｐ_ｔに対して、第一のベクトル修正工程を実行する。同工程は図６を用いて後述する。
【００３１】
ステップＳ１０９ではライン探索工程を実行する。同工程では設計変数の更新のための係数δを決定する。同工程は図３を用いて後述する。
【００３２】
ステップＳ１１０では設計変数ベクトル更新工程を実行する。同工程では設計変数ベクトルを更新する：
（式１７）ｆ_ｔ＋１＝ｆ_ｔ＋δｐ_ｔ
ステップＳ１１１では、設計変数ベクトルに対して第二のベクトル修正工程を実行する。同工程は図７を用いて後述する。
【００３３】
ステップＳ１１２ではｔをｔ＋１に更新し、更新後のｔが予め設定された反復回数を超えたら処理を終了する。
【００３４】
以上が第二の求解工程の流れである。以下、各ステップについて詳述する。
【００３５】
図３を参照してライン探索工程を説明する。同工程には、現在の設計変数値ｆ_ｔ、現在の設計変数値における探索ベクトルｐ_ｔが引数として与えられる。更に、評価関数や拘束条件は適宜与えられているとする。
【００３６】
ステップＳ３０１では最大ステップサイズΔを、次式を用いて計算する：
【外２】

【００３７】
ただし

また、各変数の初期化を次式のように実行する：
（式１９）ｇ１＝（３．０−５．０^１／２）／２，
ｇ２＝（５．０^１／２−１）／２，
ｄΔ＝Δ／Ｔ_ｍ，
ａ１＝０．０，
ａ２＝ａ１＋ｄΔ，
ただしＴ_ｍは探索範囲を適当な区間に分割するための正整数であり、例えば１０〜１００等の値を設定する。
【００３８】
ステップＳ３０２では探索範囲絞り込み工程を実行する。同工程は図４を用いて後述する。
【００３９】
ステップＳ３０３では探索範囲絞込み工程で設定されたδと最大ステップサイズΔとを比較し、等しければ処理を終了する。そうでなければステップＳ３０４へ進む。
【００４０】
ステップＳ３０４では極小点探索工程を実行する。同工程は図５を用いて後述する。以上でライン探索工程を終了する。
【００４１】
図４を参照して探索範囲絞込み工程を実行する。
【００４２】
ステップＳ４０１では第２の評価関数値を計算する：
（式２０）Ｌ２１＝Ｌ２（ｘ（ｆ_ｔ＋ａ１ｐ_ｔ），ｆ_ｔ＋ａ１ｐ_ｔ）
（式２１）Ｌ２２＝Ｌ２（ｘ（ｆ_ｔ＋ａ２ｐ_ｔ），ｆ_ｔ＋ａ２ｐ_ｔ）
ただしｘ（ｆ_ｔ＋ａ１ｐ_ｔ）、およびｘ（ｆ_ｔ＋ａ２ｐ_ｔ）は図１０を用いて後述する第一の求解工程を実行することによって得られる。
【００４３】
ステップＳ４０２ではＬ２１とＬ２２とを比較し、Ｌ２１が大きければ処理を終了し、そうでなければステップＳ４０３へ進む。
【００４４】
ステップＳ４０３ではｔを０に、μ（ｔ）を０に初期化する。
【００４５】
ステップＳ４０４ではμ（ｔ）を更新する：
（式２２）μ（ｔ）＝μ（ｔ−１）＋ｄΔ
ステップＳ４０５ではＬ２１とＬ２２を更新する：
（式２３）Ｌ２１＝Ｌ２２
（式２４）Ｌ２２＝Ｌ_２（ｘ（ｆ_ｘ＋μ（ｔ）ｐ_ｔ），ｆ_ｘ＋μ（ｔ）ｐ_ｔ）
ただしｘ（ｆ_ｘ＋μ（ｔ）ｐ_ｔ）は図１２を用いて後述する第一の求解工程を実行することによって得られる。
【００４６】
ステップＳ４０６ではＬ２１とＬ２２とを比較し、Ｌ２１がＬ２２以下ならステップＳ４０８へ、そうでなければステップＳ４０７へ進む。
【００４７】
ステップＳ４０７では、ｔをｔ＋１に更新したあと、ｔとＴ_ｍを比較し、ｔが大きければステップＳ４０９へ、そうでなければステップＳ４０４へ進む。
【００４８】
ステップＳ４０８では、次式を実行する：
（式２５）ａ１＝μ（ｔ）−２ｄΔ
ａ２＝μ（ｔ）
ステップＳ４０９ではμ（ｔ）とΔを比較し、μ（ｔ）が大きければδにΔを代入し、そうでなければ、δにΔ以外の任意の値を代入して、探索範囲絞込み工程を終了する。
【００４９】
図５を用いて極小点探索工程を説明する。
【００５０】
ステップＳ５０１では次式を計算し
（式２６）ｄａ＝ａ２−ａ１
ｄａが予め設定された微小な正の実数以下であればステップＳ５０４へ進み、そうでなければステップＳ５０２へ進む。
【００５１】
ステップＳ５０２では次式を計算する：
（式２７）ｖ１＝ａ１＋ｇ１ｄａ
ｖ２＝ａ１＋ｇ２ｄａ
更に探索範囲の絞込みを行う：
（式２８）ａ２＝ｖ２，ｉｆ Ｌ２１＜Ｌ２２
ａ１＝ｖ１， ｏｔｈｅｒｗｉｓｅ
これは、黄金分割法による区間縮小である。
【００５２】
ステップＳ５０３では、以下の式により第二の評価関数の値を計算し、ステップＳ５０１へ進む。
【００５３】
（式２９）Ｌ２１＝Ｌ２（ｘ（ｆ_ｘ＋ａ１ｐ_ｔ），ｆ_ｘ＋ａ１ｐ_ｔ）
Ｌ２２＝Ｌ２（ｘ（ｆ_ｘ＋ａ２ｐ_ｔ），ｆ_ｘ＋ａ２ｐ_ｔ）
ただしｘ（ｆ_ｔ＋ａ１ｐ_ｔ）、およびｘ（ｆ_ｔ＋ａ２ｐ_ｔ）は図１０を用いて後述する第一の求解工程を実行することによって得られる。
【００５４】
ステップＳ５０４では探索区間の中点ａｃと、その位置での第２の評価関数値Ｌ２ｃを計算する：
（式３０）ａｃ＝（ａ１＋ａ２）／２
（式３１）Ｌ２ｃ＝Ｌ２（ｘ（ｆ_ｘ＋ａｃｐ_ｔ），ｆ_ｘ＋ａｃｐ_ｔ）
ステップＳ５０５ではδの値を決定する：
（式３２）δ＝ａｃ，ｉｆ Ｌ２ｃ＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝ａ１， ｅｌｓｅ ｉｆ Ｌ２１＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝ａ２， ｅｌｓｅ ｉｆ Ｌ２２＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝Δ， ｅｌｓｅ ｉｆ Ｌ２ｅ＝ＭＩＮ［Ｌ２１，Ｌ２２，Ｌ２ｃ，Ｌ２ｅ］
＝０， ｏｔｈｅｒｗｉｓｅ
ただしＭＩＮ［］は引数の中の最小値を返す関数である。以上で極小点探索工程を終了する。
【００５５】
ステップＳ１０４およびステップＳ１０８で実施される第一のベクトル修正工程について、図６を用いて説明する。
【００５６】
ステップＳ６０１では第二の係数計算工程を実行する。ここでは、単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、および係数Ｂの値を計算する。同処理は図８を用いて詳述する。
【００５７】
ステップＳ６０２では第一のベクトル射影工程を実行する。この工程は、ステップＳ６０１で計算した単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、及び係数Ｂを使用し、次式に基づいてベクトルＸを修正する：

以上で第一のベクトル修正工程を終了する。
【００５８】
ステップＳ１１１で実施される第二のベクトル修正工程について、図７を用いて説明する。
【００５９】
ステップＳ７０１では第二の係数計算工程を実行する。同処理は図８を用いて詳述する。
【００６０】
ステップＳ７０２では第二のベクトル射影工程を実行する。この工程は、ステップＳ７０１で計算した単位法線ベクトルＮＮ、原点から超平面への垂線の長さＤＳＴ、及び係数Ｂを使用し、次式に基づいてベクトルＸを修正する：

以上で第二のベクトル修正工程を終了する。
【００６１】
図８を用いて第二の係数計算工程を説明する。
【００６２】
ステップＳ８０１では法線ベクトル計算工程を実行し、ＮＮ、およびＤＳＴを得る。同工程は図９を用いて後述する。
【００６３】
ステップＳ８０２では、次式を用いて係数Ｂを計算する：
【外３】

【００６４】
以上で第二の係数計算工程を終了する。
【００６５】
図９を用いてステップＳ８０１の法線ベクトル計算工程を説明する。本工程では、単位法線ベクトルＮＮと、原点から超平面への垂線の長さＤＳＴを計算する。
【００６６】
ステップＳ９０１では次式の変数を計算する：
（式３６）ＤＳＴ＝−ｃ（０）’
ただしｃ（０）’は次式で計算される：
【外４】

【００６７】
（式３６）のＤＳＴを次式で修正する：
（式４０）ＤＳＴ＝ＤＳＴ／Ｄ
上式で与えられるＤＳＴが、原点から超平面への垂線の長さとなる。
【００６８】
ステップＳ９０２では法線ベクトルの計算を実行する。
【００６９】
まず、正規化されていない法線ベクトルＮを（式４１）および（式４２）より計算する：
【外５】

【００７０】
（式４２）Ｎ（ｊ）＝ｆ_ｔ（ｊ）， ｏｔｈｅｒｗｉｓｅ
（式４１）を用いて、部分空間に対する単位法線ベクトルＮＮを次式で計算する：
【外６】

【００７１】
例えば、等式拘束条件として設計変数ベクトルの要素の総和が１である場合は、ｃ（０）＝−１， ｃ（ｊ）＝１，ｊ＝１，２，…，ｎであるので、単位法線ベクトルの成分は次のように簡単になる：
（式４４）ＮＮ（ｊ）＝｜∧_１｜^−１／２， ｉｆ ｊ∈∧_１
ただし｜∧_１｜は∧_１の要素数である。
【００７２】
以上で法線ベクトル計算工程を終了する。
【００７３】
図１０を用いて、第一の求解工程を説明する。
【００７４】
構造が変化すると特性関数値が０となる要素が出現する可能性がある。このような構造に対して構想解析を行うと、係数行列が正則でなくなるので、直接法では解けなくなる。そのために、第１の求解工程として、以下の３通りの方法を用いる：
（１）値が０の特性関数値は０に近い微小な値に置き換えて、直接法で連立１次方程式を解く。
【００７５】
（２）値が０を超える要素のみで連立１次方程式を構成し直して、直接法で解く。
【００７６】
（３）逆行列の計算を必要としない、勾配ベクトルを用いた反復法で解く。これには最急降下法、共役勾配法等が含まれる。
【００７７】
以下では共役勾配法を用いた求解工程を説明する。
【００７８】
ステップＳ１００１では、シミュレーション対象となる系の諸元を読み込む。系の諸元にはｘの初期値、設計パラメタｆの値、境界条件Ｂ_１、評価関数Ｌ_１が含まれている。プログラムは、この情報によって、必要な変数領域を１次記憶装置２０４に確保し、値を設定する。
【００７９】
ステップＳ１００２では反復回数ｔを１に初期化する。反復回数ｔは、上述の第一の求解工程に対する第１の反復処理に対するものであり、以下のステップＳ１００３からステップＳ１０１０の処理の反復回数が設定値を超えたら終了する。
【００８０】
ステップＳ１００３では状態変数ベクトルに関する第一の評価関数のグラジエントベクトルｇ_１，ｔを計算する。ｇ_１，ｔは次式で計算される：
（式４５）ｇ_１，ｔ＝（∂Ｌ_１／∂ｘ_１，∂Ｌ_１／∂ｘ_２，…，∂Ｌ_１／∂ｘ_ｎ）^Ｔ
∂Ｌ_１／∂ｘ_ｊは解析的に解ける場合とそうでない場合がある。解析解がある場合は、その計算工程を利用して（式４５）を得る。そうでない場合は、自動微分という技術を利用して計算することが出来る。自動微分技術は〔文献６〕等によって公知である。
【００８１】
ステップＳ１００４ではｇ_１，ｔのノルムの２乗を計算し、その値が０或いは０に近い非常に小さい値のときに処理を終了する。ノルムの２乗は次式で計算する：
（式４６）‖ｇ_１，ｔ‖^２＝ｇ_１，ｔ ^Ｔｇ_１，ｔ
ステップＳ１００５では、（式１２３）よりβを計算する：
（式４７）β＝‖ｇ_１，ｔ‖^２／‖ｇ_{１，ｔ−１}‖^２
ただしｔ＝１のときはβ＝１とする。
【００８２】
ステップＳ１００６では（式１２４）を用いて、探索方向ベクトルｐ_ｔを計算する：
（式４８）ｐ_１，ｔ＝−ｇ_{１，ｔ−１}＋βｐ_{１，ｔ−１}
ステップＳ１００７ではライン探索工程を実行する。設計変数の更新のための係数αを決定するものである。構造解析問題を有限要素モデルを用いて解く場合には、αの値は次式で与えられる：
（式４９）α＝ｐ_１，ｔ ^Ｔｇ_１，ｔ／（ｐ_１，ｔ ^ＴＡｐ_１，ｔ）
ステップＳ１００８では状態変数ベクトルを更新する：
（式５０）ｘ_ｔ＋１＝ｘ_ｔ＋αｐ_１，ｔ
ステップＳ１００９ではｔをｔ＋１に更新し、更新後のｔが予め設定された反復回数を超えたら処理を終了する。
【００８３】
以上で第一の求解工程を終了する。
【００８４】
（実施例）
本実施例では、任意の位置に加重を受ける片持ち梁の最適形状自動設計に本発明を適用するものである。説明を簡単にするために平面歪問題に限定する。
【００８５】
図１１に示すように、構造部材の存在を可能とする設計領域は、長方形１１０２であり、有限要素法に従って、該領域を縦ｎ_ｙ、横ｎ_ｘに等間隔に分割する。分割にされた部分領域をセルと呼び、左下および右上のセルが（１，１）（ｎ_ｙ，ｎ_ｘ）となるように番号付けを行う。同様に格子点をノードと呼び、左下および右上のノードが（１，１）（ｎ_ｙ＋１，ｎ_ｘ＋１）となるように番号付けを行う。
【００８６】
図中１１０１は支持部材、１１０３は荷重ベクトルである。セル（ｊ，ｋ）には特性関数値ｆ（ｊ，ｋ）が対応する。ここで特性関数値とはセル（ｊ，ｋ）における構造部材の存在確率を示す０から１の正の実数値をとる変数であり、本発明における設計変数ベクトルｆの要素である：
（式５１）ｆ＝（ｆ（１，１），ｆ（１，２），…，ｆ（ｎ_ｙ，ｎ_ｘ））^Ｔ
同様にノード（ｊ，ｋ）には横方向変位ｕ（ｊ，ｋ）と縦方向変位ｖ（ｊ，ｋ）が対応する。これらは任意の値を取る実数であり、本発明における状態変数ベクトルＵの要素である：
（式５２）Ｕ＝（ｕ（１，１），ｖ（１，１），ｕ（１，２），ｖ（１，２），…，ｕ（ｎ_ｙ＋１，ｎ_ｘ＋１），ｖ（ｎ_ｙ＋１，ｎ_ｘ＋１））^Ｔ
同様に剛性マトリクスをＡ、加重ベクトルをｂと書けば、有限要素法の良く知られた結果として、状態変数ベクトルＵは、次式で与えられる評価汎関数の最適化問題の解として与えられる：
（式５３）Ｌ_１＝（１／２）Ｕ^ＴＡＵ−ｂ^ＴＵ
より具体的には、以下の線型方程式の解として与えられることが知られている：
（式５４）ＡＵ−ｂ＝０
第一の求解工程において計算されるｇ_１，ｔは（式５４）の左辺として与えられる。（式５４）は連立一次方程式であるので、求解法として、通常、逆行列に基づく直接解法が採用される。しかし構造最適化問題においては、剛性マトリクスＡは設計変数ベクトルｆの関数であり、ｆに０の値を取る要素があれば、剛性マトリクスＡはランク落ちし正則でなくなってしまう。その結果直接解法で解けないという事態に陥るので、本発明で共役勾配法を採用しているのは、このような理由による。
【００８７】
設計変数ベクトルｆに対する評価関数Ｌ_２は総歪エネルギーで定義される。
【００８８】
（式５５）Ｌ_２＝（１／２）Ｕ^ＴＡＵ
第二の求解工程で必要となるｇ_ｔは次式で計算できる：
（式５６）ｇ_ｔ＝∂Ｌ_２／∂ｆ（ｅ）＝−（１／２）Ｕ_ｅ ^ＴＡ_ｅＵ_ｅ
ただしＵ_ｅ、Ａ_ｅは、それぞれ、要素ｅに対応するノードにおける変位ベクトル、該変位ベクトルに対応する要素剛性マトリクスである。
【００８９】
図１２に本実施例に対する計算結果を示す。図中、黒い領域が構造部材が存在する領域である。
【００９０】
本実施例において、構造部材の設計領域のアスペクト比は縦／横が２対１であり、その解析解は、水平方向に対して±４５度の梁が組み合わせられたものであることが知られている。
【００９１】
図１２に示した計算結果は、このような解析解とよく一致していることがわかる。
【００９２】
尚、本発明は、単一の機器からなる装置に適用しても、複数の機器から構成されるシステムに適用してもよい。また、上述した実施形態の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体を、装置あるいはシステムに供給し、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって達成してもよい。
【００９３】
更に、装置あるいはシステム内のコンピュータが記憶媒体に格納されたプログラムコードを読み出して実行することによって、上述した実施形態の機能を直接実現するばかりでなく、そのプログラムコードの指示に基づいて、コンピュータ上で稼動しているＯＳなどの処理により、上述の機能を実現される場合も含まれる。
【００９４】
これらの場合、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。
【００９５】
以下、上記実施形態に係わる本発明の特徴を整理する。
【００９６】
特徴１．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計方法において、
前記第２の求解工程が、前記設計変数ベクトルを逐次更新する設計変数ベクトル更新工程を含み、該設計変数ベクトル更新工程は、
始点から探索して極小点を得る極小点探索工程と、
該極小点における前記第２の評価関数の値と、終点における前記第２の評価関数値とに基づいて最適点を決定する端点評価工程とを備えることを特徴とする最適設計方法。
【００９７】
特徴２．
前記極小点探索工程は、黄金分割を用いた区間縮小法を用いることを特徴とする特徴１に記載の最適設計方法。
【００９８】
特徴３．
前記極小点探索工程は、２次曲線近似による推定法を用いることを特徴とする特徴１に記載の最適設計方法。
【００９９】
特徴４．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解手段と、
前記第１の求解手段で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解手段とを備える最適設計方法において、
前記第２の求解手段が、前記設計変数ベクトルを逐次更新する設計変数ベクトル更新手段を含み、該設計変数ベクトル更新手段は、
始点から探索して極小点を得る極小点探索手段と、
該極小点における前記第２の評価関数の値と、終点における前記第２の評価関数値とに基づいて最適点を決定する端点評価手段とを備えることを特徴とする最適設計方法。
【０１００】
特徴５．
設計変数ベクトルをパラメタとして状態変数ベクトルに対する第１の評価関数の最適化問題を解く第１の求解工程と、
前記第１の求解工程で求められた状態変数ベクトルと、前記設計変数ベクトルとに対する第２の評価関数の最適化問題を解く第２の求解工程とを備える最適設計プログラムにおいて、
前記第２の求解工程が、前記設計変数ベクトルを逐次更新する設計変数ベクトル更新工程を含み、該設計変数ベクトル更新工程は、
始点から探索して極小点を得る極小点探索工程と、
該極小点における前記第２の評価関数の値と、終点における前記第２の評価関数値とに基づいて最適点を決定する端点評価工程とを備えることを特徴とする最適設計プログラム。
【０１０１】
【発明の効果】
以上説明したように、本発明によれば、個別の問題に依存する経験則を利用すること無しに、トポロジーを含む構造最適設計を実行することができる。
【図面の簡単な説明】
【図１】第２の求解工程の処理手順を示すフローチャートである。
【図２】本実施形態に係る装置のブロック構成図である。
【図３】ライン探索工程の処理手順を示すフローチャートである。
【図４】探索範囲絞込み工程の処理手順を示すフローチャートである。
【図５】極小点探索工程の処理手順を示すフローチャートである。
【図６】第１のベクトル修正工程の処理手順を示すフローチャートである。
【図７】第２のベクトル修正工程の処理手順を示すフローチャートである。
【図８】係数Ｂの計算工程の処理手順を示すフローチャートである。
【図９】法線ベクトル計算工程の処理手順を示すフローチャートである。
【図１０】第１の求解工程の処理手順を示すフローチャートである。
【図１１】実施例の問題設定の説明図である。
【図１２】実施例の計算結果を示す図である。[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a solution of a design parameter optimization problem, and more particularly to an automatic design technique for shape optimization including the topology of a structural member.
[0002]
[Prior art]
Structural topology optimal design is a problem that determines the optimal topology and geometry of a structural member under given conditions. Hereinafter, the topology and shape dimensions of the structural member are referred to as a design variable function, and the determination problem is referred to as a design variable function optimization problem. The reason for the variable function is that the topology and shape dimensions are functions of a three-dimensional space.
[0003]
In the design variable function optimization problem, the state variable function optimization problem must be solved for each design variable function value. In this sense, the structural topology optimum design can be regarded as a double structure optimization problem having a state variable function optimization problem inside and a design variable function optimization problem outside. The inner state variable function optimization problem employs the idea of dividing the space into a finite number of elements by accumulating the prior art.
[0004]
In particular, in the problem where the strain energy of the structural member is an evaluation functional, the finite element method is generally used as the analysis method. A direct method for linear equations is adopted as a solution method of the finite element method. On the other hand, the design variable function optimization problem is roughly classified into the following three methods ([Document 1] or [Document 2]):
1. Evolutionary method (hereinafter referred to as E method)
2. Homogenization method (H method)
3. Material distribution method (MD method)
In the E method, each of the partial spaces obtained by dividing the space is referred to as a cell, and generation and disappearance are repeated according to appropriate rules. The structural member is given as a collection of cells that are finally present. By allowing only two states, whether or not a cell is present, a clear structural member is obtained. In addition, since the differential information of the evaluation functional is not used, it is not trapped by the local optimum solution, so that it is effective when the evaluation functional is multimodal. [Document 3] provides an optimization design apparatus for a frame structure member using a genetic search method which is a kind of E method. In this optimization design device, trial calculation based on accumulated know-how is required, and therefore the problems of the prior art that could not cope with the actual design problem with a large number of design variables are solved by the following method. ing. That is, an approximate optimization calculation device that uses an approximate expression of discrete design variable data such as a cross-sectional dimension of a frame member and a detailed optimization calculation device that uses the design variable data are provided, and these two calculation devices are combined to provide a framework. It constitutes an optimal design device for the structure.
[0005]
The H method makes it possible to adopt sensitivity analysis by newly introducing a design variable function that assumes a finer structure and introduces a continuous value in a structural element located in each divided subspace. Here, sensitivity analysis is an analysis method that uses differential information of an evaluation functional related to a design variable function. If sensitivity analysis becomes possible, an iterative solution method such as a gradient method can be used. Compared with the brute force method as described above, at least the calculation time for searching for the local optimum solution is greatly shortened.
[0006]
The MD method is a method of expressing the topology and shape dimension change of a structural member by assigning each element a real number in the range of 0 to 1 indicating the existence probability of the structural member. It is the same as the H method in the sense that it enables sensitivity analysis by replacing discrete information on whether or not there is a structural member with a continuous value of existence probability, but with a smaller number of parameters than the H method. Modeling is easy and the application range is wide.
[0007]
[Document 5] discloses a phase optimization method for a structure by the MD method. The features of this method are as follows:
(1) Since the voxel finite element method (space is divided at equal intervals) is used, the element stiffness matrix for all elements is the same. Therefore, if the element stiffness matrix is calculated once in advance, it can be used for subsequent calculations. Furthermore, since the elements are regularly arranged, it is not necessary to store node number information of each element.
[0008]
(2) In order to solve a large-scale simultaneous linear equation, a solution can be obtained without assembling the entire stiffness matrix by using a combination of the preconditioned conjugate gradient method and the Element-by-Element method. Less memory capacity.
[0009]
(3) The homogenization method requires six design variables (one-dimensional case) for one element. In addition, the element stiffness matrix must be recalculated each time the design variable changes. On the other hand, one design variable may be used for one element by adopting a density method that expresses the abundance of the structural member as a density ratio. Also, changes in design variables do not affect the element stiffness matrix.
[Reference 1] S. Bulman, J .; Sienz, E .; Hinton: “Comparisons between algorithms for structural topology optimization using a series of benchmark studies,” Computers and Structures. 79. 1203-1218 (2001).
[Reference 2] YL. Hsu, MS. Hsu, C-T. Chen: “Interpreting results from topology optimization using density controls,” Computers and Structures, 79, pp. 1049-1058 (2001).
[Document 3] JP 2001-134628 [Document 4] Hiroshi Yamakawa: “Optimization Design,” Computational Mechanics and CAE Series 9, Baifukan (1996)
[Document 5] Fujii, Suzuki, Otsubo: “Topology optimization of structures using voxel finite element method,” Transactions of JSCES, Paper No. 20000010 (2000).
[0010]
[Problems to be solved by the invention]
However, the prior art has the following problems.
[0011]
In general, the structural optimization problem is formulated as a double optimization problem including a state variable vector optimization problem at each iteration step of the design variable vector optimization problem. When the optimization problem of the design variable vector is called outer optimization and the optimization problem of the state variable vector is called inner optimization, the inner optimization uses the design variable vector as a parameter, that is, the design variable vector is fixed. This is a problem of obtaining a state variable vector. This is usually called structural analysis, and can be solved by solving a linear equation by the finite element method.
[0012]
However, if the structure is changed and there is no structural member in a certain region, for example, if a hole is formed in the member, the design variable of the corresponding element becomes 0, and as a result, the Young's modulus of the element becomes 0. Then, the coefficient matrix of the linear equation is not full rank and cannot be solved by the direct method because there is no inverse matrix.
[0013]
For this reason, in most conventional techniques, when the design variable vector is 0, a method of replacing it with a minute value close to 0 instead of being exactly 0 has been taken. The same applies to the voxel finite element method employed in [Document 5].
[0014]
However, setting the element design variable to a minute value is equivalent to the presence of a physically thin film or a weak member. In other words, the conventional method does not accurately represent a portion without a substance.
[0015]
[Means for Solving the Problems]
In order to solve the above-described problem, according to the present invention, a storage unit that stores a program and an execution unit that executes the program are provided, and the execution unit executes the program stored in the storage unit. First and second solving means, gradient vector calculating means, convergence determining means, first and second coefficient calculating means, search vector calculating means, design variable vector updating means, and maximum step In the structure optimum design apparatus having a size calculating means, a narrowing means, and a coefficient determining means, the first solution finding means is a design that represents the existence ratio of the structural member in each area obtained by dividing the area where the structural member can exist. A first solution step for solving the optimization problem of the first evaluation function for the state variable vector with the variable vector as a parameter, and the second solution means include the first solution In the optimal design method, comprising the state variable vector obtained in the step and a second solving step for solving the optimization problem of the second evaluation function for the design variable vector, the second solving step includes the gradient A vector calculation means calculates a gradient vector of the second evaluation function related to the design variable vector, and the convergence determination means determines whether the norm value of the gradient vector is less than a predetermined value. A convergence determination step, a first coefficient calculation step in which the first coefficient calculation means calculates a first coefficient based on a norm value of the gradient vector, and the search vector calculation means A search vector calculation step of calculating a search vector based on a coefficient of the vector and the gradient vector, and the second coefficient calculation means includes the design variable vector and the search vector A second coefficient calculation step of calculating a second coefficient based on the design variable vector, and a design variable vector in which the design variable vector update means sequentially updates the design variable vector based on the second coefficient and the search vector The second coefficient calculation step is repeatedly executed until the norm value of the gradient vector of the second evaluation function becomes less than the predetermined value in the convergence determination step. The step includes a maximum step size calculating step in which the maximum step size calculating means determines a maximum step size based on a constraint condition for the design variable vector, and a narrowing means is a search section determined by the maximum step size. The narrowing-down process of narrowing down the search range for searching for the minimum point from the inside, and the coefficient determination means, the end point and the midpoint of the search range And a coefficient determination step of determining a point that gives a smaller evaluation function value as the second coefficient based on the comparison of the second evaluation functions.
[0016]
DETAILED DESCRIPTION OF THE INVENTION
The configuration of the present embodiment is shown in FIG. In the figure, 201 is a CPU, 202 is a display device, 203 is an input device, 204 is a primary storage device, 205 is a secondary storage device, 206 is a communication device, and 207 is a bus line.
[0017]
Each processing function in the present embodiment is stored in advance in the secondary storage device 205 as a program, and is loaded into the primary storage device and executed by a command input from the input device 203 or the communication device 206 or the like. is there.
[0018]
For the following explanation, we formulate the target problem.
[0019]
If the variable function is expressed by a finite-dimensional vector by the formulation by the finite element method, the evaluation functional of the variable function becomes an evaluation function of the variable vector. Hereinafter, it is described as being expressed by a finite dimensional vector.
The state variable vector x and the design variable vector f are respectively written as column vectors as follows:
(Formula 1) x = (x ₁ , x ₂ ,..., X _m ) ^T
(Expression 2) f = (f ₁ , f ₂ ,..., F _n ) ^T
Here, T represents transposition. x is an m-dimensional vector, and f is an n-dimensional vector.
Let the boundary conditions of x and f be B ₁ and B ₂ respectively:
(Formula 3) B ₁ (x) = 0
(Formula 4) B ₂ (f) = 0
Assume that the evaluation functions related to the state variable vector x and the design variable vector f are L ₁ and L ₂ , respectively:
(Formula 5) L ₁ = L ₁ (x; f)
(Expression 6) L ₂ = L ₂ (f, x)
L ₁ is a function having x as a variable vector and f as a parameter, and L ₂ is a function having x and f as variable vectors.
[0020]
In general optimization problems, K ₁ equality constraints and K ₂ inequality constraints are often imposed. Therefore, the j-th equality constraint condition and the k-th inequality constraint condition for the design variable are respectively Q _j and R _k :
(Expression 7) Q _j (f, x) = 0
(Equation 8) R _k (f, x) ≧ 0
From the above notation, the optimum design problem is given as a solution satisfying the following expression under the constraint conditions (Expression 4), (Expression 5), (Expression 6), (Expression 7), and (Expression 8). :
(Formula 9) min [L ₂ ]
However, the state variable x is obtained as a solution satisfying the following equation under the constraint condition (Equation 3):
(Formula 10) min [L ₁ ]
Based on the above formulation, the processing of this embodiment will be described.
In the following, for ease of explanation, the equality constraint in (Equation 7) is one of the following:
[Outside 1]

[0021]
Here, c _j (j = 0, 1,..., N) is a constant.
[0022]
Also, the inequality constraint condition of (Equation 8) is for the range of values as follows:
(Equation ^{_{12) R j - (f)}} = f t (j) -f jmin ≧ 0
R _j ⁺ (f) = f _jmax −f _t (j) ≧ 0
However, j = 1, 2,..., N.
[0023]
Therefore, the number of inequality constraints is twice the number of design variables.
[0024]
Define flag vectors a ⁺ and a ⁻ of the same dimension as the design variable vector as follows:
a ⁺ = 0, if 0 ≧ R _j ⁺ (f)
a ⁺ = 1, otherwise
And a ⁻ = 0, if 0 ≧ R _j ⁻ (f)
a ⁻ = 1, otherwise
Further, a ⁺ and a ^- is ∧ a set of j are both 1 _1, the set of the other indexes and ∧ _2.
[0025]
FIG. 1 shows a flowchart of the second solution solving process. In the figure, step S101 is a process of reading the specifications of the system to be simulated. The specification can be read by input data from the input device 203 or the communication device 206, or data stored as a file in the secondary storage device 205 in advance can be read and used. The system specifications include initial values of x and f, boundary conditions B ₁ and B ₂ , evaluation functions L ₁ and L ₂ , and constraint conditions Q _j and R _k . Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value.
[0026]
In step S102, the number of iterations t is initialized to 1. The number of iterations t is the number of times the first iteration process has been performed in the above-described second solving step, and ends when the number of iterations of the following steps S103 to 111 reaches the set value T.
[0027]
In step S103, a gradient vector calculation step is executed. In the process, calculating a gradient vector g _t of the second evaluation functions for the design variable vector. g _t is calculated as:
(Formula 13) g _t = (∂L ₂ / ∂f ₁ , ２L ₂ / ， f ₂ ,..., L ₂ / ∂f _n ) ^T
∂L ₂ / ∂f _j may or may not be solved analytically. When there is an analytical solution, (Equation 10) is obtained using the calculation process. Otherwise, it can be calculated using a technique called automatic differentiation. The automatic differentiation technique is known from [Document 6] below.
[0028]
[Reference 6] Koichi Kubota, Masao Iri: "Automatic differentiation and application of algorithms," Modern Nonlinear Science Series, Corona (1998)
In step S104, it executes the obtained _{g t} (Equation 13), the first vector correction process. This process will be described later with reference to FIG.
[0029]
In step S105, a convergence determination step is executed. In the process, calculates the square of the norm of g _t, the process ends and the value is very less than small value ε close to 0 or 0. The square of the norm is calculated as:
(Expression 14) ‖g _t ‖ ² = g _t ^T g _t
In step S106, a first coefficient calculation step is executed. In the first coefficient calculation step, the coefficient β is calculated by (Equation 15):
(Formula 15) β = ‖g _t ‖ ² / ‖g _t-1 ‖ ²
However, when t = 1, β = 1.
[0030]
In step S107, a search vector calculation step is executed. In the same process, the search vector pt of the _{tth iteration} is calculated using (Equation 16):
(Expression 16) p _t ← −g _t−1 + βp _t−1
In step S108, the obtained _{p t} from the equation (16), performing a first vector correction process. This process will be described later with reference to FIG.
[0031]
In step S109, a line search process is executed. In the same process, a coefficient δ for updating the design variable is determined. This process will be described later with reference to FIG.
[0032]
In step S110, a design variable vector update process is executed. The process updates the design variable vector:
(Expression 17) f _{t + 1} = f _t + δp _t
In step S111, a second vector correction step is performed on the design variable vector. This process will be described later with reference to FIG.
[0033]
In step S112, t is updated to t + 1, and the process ends when the updated t exceeds the preset number of iterations.
[0034]
The above is the flow of the second solving process. Hereinafter, each step will be described in detail.
[0035]
The line search process will be described with reference to FIG. In the same process, the current design variable value f _t and the search vector p _t in the current design variable value are given as arguments. Furthermore, it is assumed that the evaluation function and constraint conditions are given as appropriate.
[0036]
In step S301, the maximum step size Δ is calculated using the following equation:
[Outside 2]

[0037]
However,

Also, initialize each variable as follows:
(Formula 19) g1 = (3.0−5.0 ^1/2 ) / 2
g2 = (5.0 ^1/2 -1) / 2
dΔ = Δ / T _m ,
a1 = 0.0,
a2 = a1 + dΔ,
However, _Tm is a positive integer for dividing the search range into appropriate sections, and is set to a value such as 10 to 100, for example.
[0038]
In step S302, a search range narrowing step is executed. This process will be described later with reference to FIG.
[0039]
In step S303, δ set in the search range narrowing step is compared with the maximum step size Δ, and if they are equal, the process ends. Otherwise, the process proceeds to step S304.
[0040]
In step S304, a minimum point search step is executed. This process will be described later with reference to FIG. This completes the line search process.
[0041]
The search range narrowing step is executed with reference to FIG.
[0042]
In step S401, a second evaluation function value is calculated:
(Equation _{20) L21 = L2 (x (} f t + a1p t), f t + a1p t)
(Equation _{21) L22 = L2 (x (} f t + a2p t), f t + a2p t)
However _{x (f t + a1p t)} , and _x (f _t + a2p t) is obtained by performing a first solution finding process to be described later with reference to FIG. 10.
[0043]
In step S402, L21 and L22 are compared. If L21 is larger, the process ends. If not, the process proceeds to step S403.
[0044]
In step S403, t is initialized to 0 and μ (t) is initialized to 0.
[0045]
In step S404, μ (t) is updated:
(Expression 22) μ (t) = μ (t−1) + dΔ
In step S405, L21 and L22 are updated:
(Formula 23) L21 = L22
(Expression 24) L22 = L ₂ (x (f _x + μ (t) p _t ), f _x + μ (t) p _t )
However, x (f _x + μ (t) p _t ) is obtained by executing a first solution solving process to be described later with reference to FIG.
[0046]
In step S406, L21 and L22 are compared. If L21 is equal to or smaller than L22, the process proceeds to step S408, and if not, the process proceeds to step S407.
[0047]
In step S407, t is updated to t + 1, t and _Tm are compared, and if t is large, the process proceeds to step S409, and if not, the process proceeds to step S404.
[0048]
In step S408, the following equation is executed:
(Equation 25) a1 = μ (t) −2dΔ
a2 = μ (t)
In step S409, μ (t) and Δ are compared. If μ (t) is large, Δ is substituted for δ, and if not, any value other than Δ is substituted for δ, and the search range narrowing step is performed. finish.
[0049]
The minimum point search process will be described with reference to FIG.
[0050]
In step S501, the following expression is calculated (Expression 26) da = a2-a1.
If da is less than or equal to a preset small positive real number, the process proceeds to step S504, and if not, the process proceeds to step S502.
[0051]
In step S502, the following equation is calculated:
(Formula 27) v1 = a1 + g1da
v2 = a1 + g2da
Further refine the search range:
(Formula 28) a2 = v2, if L21 <L22
a1 = v1, otherwise
This is a section reduction by the golden section method.
[0052]
In step S503, the value of the second evaluation function is calculated by the following formula, and the process proceeds to step S501.
[0053]
(Equation _{29) L21 = L2 (x (} f x + a1p t), f x + a1p t)
_{L22 = L2 (x (f x} + a2p t), f x + a2p t)
However _{x (f t + a1p t)} , and _x (f _t + a2p t) is obtained by performing a first solution finding process to be described later with reference to FIG. 10.
[0054]
In step S504, the midpoint ac of the search section and the second evaluation function value L2c at that position are calculated:
(Formula 30) ac = (a1 + a2) / 2
(Equation _{31) L2c = L2 (x (} f x + acp t), f x + acp t)
In step S505, the value of δ is determined:
(Expression 32) δ = ac, if L2c = MIN [L21, L22, L2c, L2e]
= A1, else if L21 = MIN [L21, L22, L2c, L2e]
= A2, else if L22 = MIN [L21, L22, L2c, L2e]
= Δ, else if L2e = MIN [L21, L22, L2c, L2e]
= 0, otherwise
However, MIN [] is a function that returns the minimum value in the argument. This is the end of the minimum point search process.
[0055]
The first vector correction process performed in step S104 and step S108 will be described with reference to FIG.
[0056]
In step S601, the second coefficient calculation step is executed. Here, the unit normal vector NN, the length DST of the perpendicular from the origin to the hyperplane, and the value of the coefficient B are calculated. This process will be described in detail with reference to FIG.
[0057]
In step S602, the first vector projection process is executed. This process uses the unit normal vector NN calculated in step S601, the perpendicular length DST from the origin to the hyperplane, and the coefficient B, and corrects the vector X based on the following equation:

This completes the first vector correction process.
[0058]
The second vector correction process performed in step S111 will be described with reference to FIG.
[0059]
In step S701, a second coefficient calculation step is executed. This process will be described in detail with reference to FIG.
[0060]
In step S702, a second vector projection process is executed. This process uses the unit normal vector NN calculated in step S701, the perpendicular length DST from the origin to the hyperplane, and the coefficient B, and corrects the vector X based on the following equation:

This completes the second vector correction process.
[0061]
The second coefficient calculation step will be described with reference to FIG.
[0062]
In step S801, a normal vector calculation step is executed to obtain NN and DST. This process will be described later with reference to FIG.
[0063]
In step S802, the coefficient B is calculated using the following equation:
[Outside 3]

[0064]
The second coefficient calculation step is thus completed.
[0065]
The normal vector calculation process in step S801 will be described with reference to FIG. In this step, a unit normal vector NN and a perpendicular length DST from the origin to the hyperplane are calculated.
[0066]
In step S901, the following variable is calculated:
(Expression 36) DST = −c (0) ′
Where c (0) ′ is calculated as:
[Outside 4]

[0067]
Modify the DST in (Equation 36) with the following equation:
(Formula 40) DST = DST / D
The DST given by the above equation is the length of the perpendicular from the origin to the hyperplane.
[0068]
In step S902, a normal vector is calculated.
[0069]
First, an unnormalized normal vector N is calculated from (Equation 41) and (Equation 42):
[Outside 5]

[0070]
(Equation _{42) N (j) = f} t (j), otherwise
Using (Equation 41), the unit normal vector NN for the subspace is calculated by the following equation:
[Outside 6]

[0071]
For example, when the sum of the elements of the design variable vector is 1 as an equality constraint condition, c (0) = − 1, c (j) = 1, j = 1, 2,. The normal vector components are simplified as follows:
(Formula 44) NN (j) = | ∧ ₁ | ^−1/2 , if j∈∧ ₁
However, | ∧ ₁ | is the number of elements of ∧ ₁ .
[0072]
This completes the normal vector calculation step.
[0073]
The first solution finding process will be described with reference to FIG.
[0074]
When the structure changes, an element having a characteristic function value of 0 may appear. When concept analysis is performed on such a structure, the coefficient matrix is not regular and cannot be solved by the direct method. Therefore, the following three methods are used as the first solution process:
(1) The characteristic function value of 0 is replaced with a minute value close to 0, and the simultaneous linear equations are solved by the direct method.
[0075]
(2) Reconstruct a simultaneous linear equation with only elements whose value exceeds 0 and solve it by the direct method.
[0076]
(3) Solve by an iterative method using gradient vectors, which does not require calculation of the inverse matrix. This includes the steepest descent method, the conjugate gradient method, and the like.
[0077]
Hereinafter, a solution finding process using the conjugate gradient method will be described.
[0078]
In step S1001, the specifications of the system to be simulated are read. The system specifications include an initial value of x, a value of the design parameter f, a boundary condition B ₁ , and an evaluation function L ₁ . Based on this information, the program secures a necessary variable area in the primary storage device 204 and sets a value.
[0079]
In step S1002, the number of iterations t is initialized to 1. The number of iterations t is for the first iteration process for the first solution step described above, and ends when the number of iterations of the processing from step S1003 to step S1010 below exceeds a set value.
[0080]
In step S1003, a gradient vector g _{1, t} of the first evaluation function related to the state variable vector is calculated. g _{1, t} is calculated by the following formula:
(Equation 45) g _{1, t} = (∂L ₁ / ∂x ₁ , ∂L ₁ / ∂x ₂ ,..., L ₁ / ｎx _n ) ^T
∂L ₁ / ∂x _j may or may not be solved analytically. When there is an analytical solution, (Equation 45) is obtained using the calculation process. Otherwise, it can be calculated using a technique called automatic differentiation. The automatic differentiation technique is known from [Document 6] and the like.
[0081]
In step S1004, the norm square of g1 _{, t} is calculated, and the process ends when the value is 0 or a very small value close to 0. The square of the norm is calculated as:
(Equation 46) ‖g _{1, t} ‖ ^{_{2 = g 1, t T g}} 1, t
In step S1005, β is calculated from (Equation 123):
(Formula 47) β = ‖g _{1, t} ‖ ² / ‖g _{1, t-} 1−１ ²
However, when t = 1, β = 1.
[0082]
In step S1006, the search direction vector p _t is calculated using (Equation 124):
(Formula 48) p1 _{, t} = -g1 _{, t-1} + (beta) p1 _{, t-1}
In step S1007, a line search process is executed. The coefficient α for updating the design variable is determined. When solving a structural analysis problem using a finite element model, the value of α is given by:
(Formula 49) α = p _{1, t} ^T g _{1, t} / (p _{1, t} ^T Ap _{1, t} )
In step S1008, the state variable vector is updated:
(Expression 50) x _{t + 1} = x _t + αp _{1, t}
In step S1009, t is updated to t + 1, and the process ends when the updated t exceeds the preset number of iterations.
[0083]
This completes the first solution process.
[0084]
(Example)
In this embodiment, the present invention is applied to the optimum shape automatic design of a cantilever beam that receives a load at an arbitrary position. For simplicity of explanation, we will limit it to the plane distortion problem.
[0085]
As shown in FIG. 11, the design area that allows the existence of the structural member is a rectangle 1102, and the area is divided into vertical n _y and horizontal _nx at equal intervals according to the finite element method. Called partial region which is to split the cell, performs numbered as lower left and upper right of the cell becomes _{(1,1) (n y, n} x). Similarly called grid points and node performs numbered as the lower-left and upper-right node is _{(1,1) (n y + 1} , n x +1).
[0086]
In the figure, 1101 is a support member, and 1103 is a load vector. The characteristic function value f (j, k) corresponds to the cell (j, k). Here, the characteristic function value is a variable taking a positive real value from 0 to 1 indicating the existence probability of the structural member in the cell (j, k), and is an element of the design variable vector f in the present invention:
(Equation 51) f = (f (1,1), f (1,2),..., F ( _ny , _nx )) ^T
Similarly, the node (j, k) corresponds to the lateral displacement u (j, k) and the longitudinal displacement v (j, k). These are real numbers taking arbitrary values, and are elements of the state variable vector U in the present invention:
(Equation 52) U = (u (1,1 ), v (1,1), u (1,2), v (1,2), ..., u (n y + 1, n x +1), v ( n _y +1, n _x +1)) ^T
Similarly, if the stiffness matrix is written as A and the weighting vector is written as b, as a well-known result of the finite element method, the state variable vector U is given as a solution to the optimization problem of the evaluation functional given by the following equation:
(Formula 53) L ₁ = (1/2) U ^T AU−b ^T U
More specifically, it is known to be given as a solution to the following linear equation:
(Formula 54) AU−b = 0
G _{1, t} calculated in the first solving process is given as the left side of (Formula 54). Since (Formula 54) is a simultaneous linear equation, a direct solution based on an inverse matrix is usually adopted as a solution. However, in the structure optimization problem, the stiffness matrix A is a function of the design variable vector f, and if there is an element in which f takes a value of 0, the stiffness matrix A drops in rank and becomes non-regular. As a result, the solution cannot be solved by the direct solution, and the conjugate gradient method is adopted in the present invention for the above reason.
[0087]
Evaluation function L ₂ with respect to the design variable vector f is defined by the total strain energy.
[0088]
(Formula 55) L ₂ = (1/2) U ^T AU
The g _t required in the second solution finding process can be calculated by the following formula:
(Expression 56) g _t = ∂L ₂ / ∂f (e) =-(1/2) U _e ^TA _e U _e
Here, U _e and A _e are a displacement vector at a node corresponding to the element e and an element stiffness matrix corresponding to the displacement vector, respectively.
[0089]
FIG. 12 shows the calculation results for this example. In the figure, black areas are areas where structural members are present.
[0090]
In this embodiment, the aspect ratio of the design area of the structural member is 2/1 in the vertical / horizontal direction, and the analytical solution is known to be a combination of beams of ± 45 degrees with respect to the horizontal direction. ing.
[0091]
It can be seen that the calculation result shown in FIG. 12 is in good agreement with such an analytical solution.
[0092]
Note that the present invention may be applied to an apparatus composed of a single device or a system composed of a plurality of devices. Also, a storage medium storing software program codes for realizing the functions of the above-described embodiments is supplied to the apparatus or system, and a computer in the apparatus or system reads out and executes the program code stored in the storage medium. May be achieved.
[0093]
Further, the computer in the apparatus or system reads out and executes the program code stored in the storage medium, thereby not only directly realizing the functions of the above-described embodiments but also on the computer based on the instruction of the program code. The case where the above-described functions are realized by processing of an OS or the like that is running on is also included.
[0094]
In these cases, the storage medium storing the program code constitutes the present invention.
[0095]
The features of the present invention according to the above embodiment will be summarized below.
[0096]
Features 1.
A first solution step for solving the optimization problem of the first evaluation function for the state variable vector with the design variable vector as a parameter;
In an optimal design method comprising: a state variable vector obtained in the first solution solving step; and a second solution solving step that solves an optimization problem of a second evaluation function for the design variable vector,
The second solving step includes a design variable vector update step of sequentially updating the design variable vector, and the design variable vector update step includes:
Searching for a minimum point by searching from the start point and obtaining a minimum point;
An optimal design method comprising: an endpoint evaluation step for determining an optimal point based on the value of the second evaluation function at the minimum point and the second evaluation function value at the end point.
[0097]
Feature 2.
The optimal design method according to claim 1, wherein the local minimum point searching step uses a section reduction method using a golden section.
[0098]
Feature 3.
2. The optimum design method according to claim 1, wherein the minimum point search step uses an estimation method based on quadratic curve approximation.
[0099]
Feature 4.
First solution means for solving the optimization problem of the first evaluation function for the state variable vector using the design variable vector as a parameter;
In an optimum design method comprising: a state variable vector obtained by the first solution finding means; and a second solution finding means for solving an optimization problem of a second evaluation function for the design variable vector,
The second solving means includes design variable vector update means for sequentially updating the design variable vector, and the design variable vector update means includes:
A minimum point search means for searching from the start point to obtain a minimum point;
An optimum design method comprising: an endpoint evaluation means for determining an optimum point based on the value of the second evaluation function at the minimum point and the second evaluation function value at the end point.
[0100]
Feature 5.
A first solution step for solving the optimization problem of the first evaluation function for the state variable vector with the design variable vector as a parameter;
In an optimal design program comprising: a state variable vector obtained in the first solving step; and a second solving step for solving an optimization problem of a second evaluation function for the design variable vector,
The second solving step includes a design variable vector update step of sequentially updating the design variable vector, and the design variable vector update step includes:
Searching for a minimum point by searching from the start point and obtaining a minimum point;
An optimum design program comprising: an endpoint evaluation step for determining an optimum point based on the value of the second evaluation function at the minimum point and the second evaluation function value at the end point.
[0101]
【The invention's effect】
As described above, according to the present invention, it is possible to execute a structure optimum design including a topology without using an empirical rule depending on an individual problem.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a processing procedure of a second solution solving process.
FIG. 2 is a block configuration diagram of an apparatus according to the present embodiment.
FIG. 3 is a flowchart showing a processing procedure of a line search process.
FIG. 4 is a flowchart showing a processing procedure of a search range narrowing step.
FIG. 5 is a flowchart showing a processing procedure of a minimum point search process.
FIG. 6 is a flowchart showing a processing procedure of a first vector correction step.
FIG. 7 is a flowchart showing a processing procedure of a second vector correction step.
FIG. 8 is a flowchart showing a processing procedure of a coefficient B calculation process.
FIG. 9 is a flowchart showing a processing procedure of a normal vector calculation step.
FIG. 10 is a flowchart showing a processing procedure of a first solution finding process.
FIG. 11 is an explanatory diagram of problem setting according to the embodiment.
FIG. 12 is a diagram showing a calculation result of an example.

Claims (2)

Storage means for storing a program, and a execution means for executing the program, said execution unit is realized by executing a program stored in the storage means, the first and second solving means, the gradient Vector calculation means, convergence determination means, first and second coefficient calculation means, search vector calculation means, design variable vector update means, maximum step size calculation means, narrowing means, and coefficient determination means In the structural optimum design apparatus having
The first solving means solves the optimization problem of the first evaluation function for the state variable vector using the design variable vector representing the existence rate of the structural member in each region obtained by dividing the region where the structural member can exist as a parameter. 1 solution process,
Optimum design, wherein the second solution solving means includes a state variable vector obtained in the first solution solving process and a second solution solving process for solving an optimization problem of a second evaluation function for the design variable vector. In the method
The second solving step includes
A gradient vector calculating step in which the gradient vector calculating means calculates a gradient vector of the second evaluation function with respect to the design variable vector;
A convergence determination step in which the convergence determination means determines whether the norm value of the gradient vector is less than a predetermined value;
A first coefficient calculating step in which the first coefficient calculating means calculates a first coefficient based on a norm value of the gradient vector;
A search vector calculation step in which the search vector calculation means calculates a search vector based on the first coefficient and the gradient vector;
A second coefficient calculating step in which the second coefficient calculating means calculates a second coefficient based on the design variable vector and the search vector;
Said design variable vector updating means, and a design variable vector updating step of sequentially updating the design variable vector based on said search vector and the second coefficient,
The second solving step is repeatedly executed until the norm value of the gradient vector of the second evaluation function is less than the predetermined value in the convergence determining step,
The second coefficient calculation step includes:
A maximum step size calculating step in which the maximum step size calculating means determines a maximum step size based on a constraint condition on the design variable vector ;
The narrowing-down process narrows down the search range for searching for a minimum point from the search section determined from the maximum step size;
The coefficient determination means includes a coefficient determination step of determining a point that gives a smaller evaluation function value as the second coefficient based on the comparison of the second evaluation function at the end point and the middle point of the search range. An optimal design method characterized by that.   The optimal design method according to claim 1, wherein the narrowing-down step uses a section reduction method using a golden section.

Images