幾何学的非線形性を考慮したシェルの形状最適化
Shape Optimization of Shells Considering Geometric Nonlinearities◦ 森友 義弥 ∗1 ，藤田 慎之輔 ∗2
Yoshiya MORITOMO∗1 ，Shinnosuke FUJITA∗2∗1 北九州市立大学大学院 国際環境工学研究科 大学院生
Graduate Student, Faculty of Environmental Engineering, The University of Kitakyushu∗2 北九州市立大学大学院 国際環境工学研究科 准教授 工博
Associate Professor, Faculty of Environmental Engineering, The University of Kitakyushu, Dr. Eng.
キーワード
キーワード：機械学習; 非線形座屈荷重; 幾何学的非線形解析

Keywords：Machine learning; Nonlinear buckling load; Geometric nonlinear analysis

1. はじめに[

近年様々な最適化アルゴリズムによって高い力学的合理性を有するシェルの形状，断面またはトポロジーを獲得する研究が行われている．しかしながらそれらの既往[r= x]

yz

研究の殆どは構造計画の初期段階での形態創生を想定し

x=

y=z=x1[y1[z1·········

]⊤xn
]⊤
yn
]⊤
zn

(4)ているため，機械的性能は線形弾性範囲内でのみ評価されている．一方，シェル構造は他の構造物と比較して幾

この時，節点座標は以下のように離散化される．

何学的非線形性の影響を受けやすいため，既往の研究方法で得られた解形状は非線形座屈荷重が低い場合がある．しかしながら，非線形解析は計算コストが非常に大

xk =

yk =

(5a)

I ∑
J∑

qy,ij Bi,N (uk )Bj,M (vk )(5b)qz,ij Bi,N (uk )Bj,M (vk )(5c)

i=0 j=0

算コストを削減するために，機械学習を使用した非線形解析の応力予測モデルを構築する．従来の構造最適化手

qx,ij Bi,N (uk )Bj,M (vk )i=0 j=0

きいため従来の最適化手法に直接導入することは現実的ではない．従って本研究では非線形性を考慮する際の計I ∑J∑zk =I ∑J∑
i=0 j=0法における繰り返し計算を機械学習による予測システムに置き換えることで，現実的な計算時間での幾何学的非ここで k = 1, · · · , n であり uk と vk は，k 番目の節点線形性を考慮した構造最適化を行う．の節点座標に対応するパラメータである．Bi,N と Bj,M
は B スプライン基底関数 1) であり，N 及び M は関数の次数である．

2. 解析モデルの作成
シェル構造の節点座標は (I + 1) × (J + 1) の制御点ネットを持つ B スプライン曲面で離散化される．各制
御点の x，y ，z 座標のベクトルを qx ，qy ，qz とするとシェル構造物の非線形座屈荷重係数を予測するためにニューラルネットワークを用いる．ニューラルネット

次式のように表される．

[qx = qx,00[qy = qy,00[qz = qz,00

3. ニューラルネットワーク回帰···

qx,0J···

qx,I0······

qy,0J···

qy,I0······

qz,0J···

qz,I0···]⊤qx,IJ (1)]⊤qy,IJ (2)]⊤qz,IJ (3)

節点数を n とすると，節点座標ベクトル r は以下のように表される．

ワークとは図 1 のように，微分可能な変換を繋げて作られた計算グラフである 2) ．
図 1 において円で表されたノードに数値が入っており，このノードの縦方向の集合を層と呼ぶ．また一層目は入力変数を表しており入力層と呼ぶ．最後の層には予測値が出力されるため出力層と呼ぶ．入力層と出力層の間にある層を中間層と呼び，中間層の数によって予測精度は変化する．計算過程の例として，入力層と第一中間

日本建築学会情報システム技術委員会
第 44 回情報･システム･利用･技術シンポジウム論文集， 443-446 ， 2021 年 12 月，京都
Proceedings of the 44th Symposium on Computer Technology of Information ，Systems and Applications ， AIJ ， 443-446 ， Dec. ， 2021 ， Kyoto－ 443 －

層における計算過程を示す．入力変数の個数を n 個，中個とすると，ニューラルネットワークは以下の式によっ間層のノードの個数を k 個とすると，てパラメータベクトルを決定する．

u1 = W 1 h 0[]⊤u1 = u11 · · · u1k
w1,11 w1,12 · · · w1,1n−1 w1,1n .
..
..
.. ..

W1 = .
.
.
. 
 ..
w1,k1 w1,k2 · · · w1,kn−1 w1,kn[]⊤h0 = h01 h01 · · · h0n−1 h0n

1 ∑(tn − yn )2
N n=1[]
y = y1 · · · y n[]t = t1 · · · t nN(6)minimize
W

(7)

(8)

(12)
(13)
(14)

yn：n 番目の予測値, tn：n 番目の目標値(9)

求めたパラメータベクトルと予測モデルの入力変数を入と表すことができる．ここで，h0 は入力層の入力変数力することによって，予測モデルの予測値を求めることベクトル，W1 は一層分のパラメータベクトル，u1 は第ができる．

一中間層の出力ベクトルである．図 1 に示すように，入力と出力の関係が非線形である場合は中間層のノード内

4. 機械学習による予測例題

で出力ベクトルから入力ベクトルへの変換を行う必要が機械学習を行うに当たり，45000 個の学習モデルとある．この変換を行うための関数を活性化関数と呼び，5000 個の予測モデルを作成した．各モデルはそれぞれ

本論文では Reacted Linear Unit(ReLU，図 2 参照) という20m スパンで 5×5 の制御点ネットにより形成された 1/8活性化関数を用いる．ReLU は以下の式で表される．

対称のシェルであり，図 3 は作成したモデルの例である．

[h1 = h11]

···

(10)

h1kn = 1, · · · , kh1n = max(0, u1n )(11)

本研究において式 (12) における目標値は幾何学的非線形解析により求めた非線形座屈荷重係数とする．非線形解析は OpenSeesPy3) を使用し，増分ステップ 0.05 の荷重

ReLU は入力が負の値の場合出力が０で一定となり，入増分法を用いる．入力変数は節点座標 + 制御点座標とす

力が正の値の場合そのまま出力する関数である．以上のる．なお各モデルは 1/8 対称であるため入力変数の各節計算を出力層が１となるまで繰り返すことによって，予点座標及び制御点座標は 1/8 の範囲のものを使用する．

測値 y を求めることができる．教師データの個数を N

またエポック数，バッチサイズは予測精度を左右するパラメータ 2) であるため，それぞれ，[50,60,70,80,90,100]，[4,8,16,32,64] とすることで予測精度の比較検討を行う．
層数はニューラルネットワークにおける中間層の数を示しており，本研究では 3 層とする．バッチサイズとは学習を行うに当たり学習モデルを分割するサイズのことで
20m

20mPerth

図1

ニューラルネットワークのイメージ図:Pin support

4.0

Elevation

3.5
3.0

h

2.5

h=max(0,u)

2.0
1.5

perth(Control point)

1.00.5Control point net

0.043

図2報告 H60

210u

123

Reacted Linear Unit

4I+1=J+1=5
N=M=2

図3

Elevation(Control point)

1/8 対称シェルの例－ 444 －

表120m

入力変数が節点座標 + 制御点座標の予測結果
層数＝ 3エポック数 バッチサイズ
エポック数=50 バッチサイズ=4エポック数=50 バッチサイズ=8エポック数=50 バッチサイズ=16エポック数=50 バッチサイズ=32エポック数=50 バッチサイズ=64エポック数=60 バッチサイズ=4エポック数=60 バッチサイズ=8エポック数=60 バッチサイズ=16エポック数=60 バッチサイズ=32エポック数=60 バッチサイズ=64エポック数=70 バッチサイズ=4エポック数=70 バッチサイズ=8エポック数=70 バッチサイズ=16エポック数=70 バッチサイズ=32エポック数=70 バッチサイズ=64%
53.64
54.80
68.06
67.00
70.64
54.64
54.54
66.68
72.20
68.24
55.28
70.08
63.16
68.16
68.98

エポック数 バッチサイズ
エポック数=80 バッチサイズ=4エポック数=80 バッチサイズ=8
エポック数=80 バッチサイズ=16
エポック数=80 バッチサイズ=32
エポック数=80 バッチサイズ=64
エポック数=90 バッチサイズ=4エポック数=90 バッチサイズ=8
エポック数=90 バッチサイズ=16
エポック数=90 バッチサイズ=32
エポック数=90 バッチサイズ=64
エポック数=100 バッチサイズ=4
エポック数=100 バッチサイズ=8
エポック数=100 バッチサイズ=16
エポック数=100 バッチサイズ=32
エポック数=100 バッチサイズ=64%
52.54
54.12
66.60
70.96
71.50
52.04
56.22
74.32
70.64
68.36
55.62
55.90
57.18
70.44
73.10

20m

:Pin support8m

E=9.87kNm
λp =80.86
λ=116.5λp /λ=0.69
V=46.89m3

19Validation MAE

18I+1=J+1=5
N=M=21716

図515135. 機械学習を用いた形状最適化

12

前節により作成したエポック数=90 バッチサイズ=16110

図4

初期形状の諸量と制御点

142040Epoch6080

収束結果 (層数＝ 3:エポック数＝ 90:バッチサイズ=16)の予測機を用いて非線形座屈荷重係数を最大化する形状最適化を行う．また比較対象として表 1 より二番目に予測精度の高かったエポック数=100 バッチサイズ=64 のある．例えば 100 個の学習モデルに対しバッチサイズを予測機 (予測精度 73.10%)，三番目に予測精度の高かっ4 とすると，バッチ数は 25 個となる．この 25 個のバッたエポック数=60 バッチサイズ=32 の予測機 (予測精度チ学習を終えたとき 1 エポックとなる．エポック数とは72.20%)，四番目に予測精度の高かったエポック数=80

1 エポックごとにデータをシャッフルして組み合わせをバッチサイズ=64 の予測機 (予測精度 71.50%) を用いた

変えながら学習を繰り返すための繰り返し回数のことで形状最適化も行う．以下に最適化問題を定式化する．

ある．
入力変数が節点座標 + 制御点座標の場合の予測結果を度にばらつきがあることが確認できる．バッチサイズに関しては層数やエポック数に関わらずバッチサイズが小さいほど予測精度が低くバッチサイズが大きくなるにつれて予測精度が高くなる傾向がある．
また，最大予測精度はエポック数=90 バッチサイズ=16 のときの 74.32% である．図 4 に最大予測精度時の収束履歴を示す．縦軸は平均絶対許容誤差 (MAE) であり数値が小さいほど予測値と目標値の誤差が少ない事を示す．横軸はエポック数であり，図は繰り返し計算毎の誤差を示している．収束過程については小刻みに振動しているものの繰り返し計算を進めるにつれて誤差が収束している事がわかる．小刻みに振動している理由として, 入力と出力の関係性が複雑でありノイズが発生しているためと推察できる．

報告 H60(15a)subject to V − V ≤ 0(15b)qz

表 1 に示す．表 1 よりエポック数に関しては，エポック数が増減することによる予測精度の変化はなく，予測精maxmize λpλpqz
VV::::

非線形座屈荷重係数の予測値
制御点の z 座標ベクトル
総体積総体積の上限値定式化した最適化問題を解くためのソルバーは pyOpt4)を用いる．
最適化例題として最適化における初期形状を図 5 に示す．初期形状は 20m 四方で高さ 8m の 4 点ピン支持の裁断球殻シェルであり，2 次の B スプライン曲面によって形成されている．荷重条件としては，各節点に鉛直下向きに 1kN の集中荷重を加えている．シェルのヤング率，ポアソン比，板厚はそれぞれ 21GPa，0.2，0.1m であり，λ は幾何学的非線形解析によって算出した真の非線形座屈荷重係数である．λp は予測機を用いて算出した非線形座屈荷重係数の予測値である．λp /λ はシェルの形状－ 445 －

に関わらず非線形座屈荷重係数の予測値 λp は最大化す

E=9.17kNm
λp=131.76
λ＝104.4λp/λ=1.26
V=46.89m3
図6

エポック数=90 バッチサイズ=16 の予測機を用いた最適化後の形状とその諸量ることが可能である．一方，最適化後の真の非線形座屈荷重係数 λ は使用する予測機によって最大化できる場合とできない場合がある．λp /λ では使用する予測機の予測精度が高いほど誤差が大きくなる傾向を示している．
また λ が初期形状より最大化されている結果を示した
ケースは，λp /λ の誤差が約 15% 以内となっている傾向がある．つまり，使用する予測機の予測精度が高ければ

図7

E=14.16kNm
λp=92.09
λ＝113.10λp/λ=0.81
V=46.89m3
エポック数=100 バッチサイズ=64 の予測機を用いた最適化後の形状とその諸量

必ずしも真の非線形座屈荷重係数が最大化できるとは限らず，λp /λ の誤差が約 15% 以内である場合に真の非線形座屈荷重係数は最大化できる傾向がある．

6. おわりに本研究では機械学習を用いることでシェル構造物における非線形座屈荷重の予測モデルの構築し，非線形座屈荷重係数を最大化させるような形状最適化を行った．得

図8E=14.45kNm
λp=116.68
λ＝137.10λp/λ=0.85
V=46.89m3
エポック数=60 バッチサイズ=32 の予測機を用いた最られた知見を以下に示す．

• 入力変数が節点座標 + 制御点座標の場合，層数やエポック数に関わらずバッチサイズが増加するにつれて予測精度が高くなる傾向がある．

• 使用した予測機の予測精度に関わらず非線形座屈適化後の形状とその諸量荷重係数の予測値 λp は最大化することが可能であることがわかる．

E=9.68kNm
λp=131.02
λ＝155.70λp/λ=0.84
V=46.89m3
図9エポック数=80 バッチサイズ=64 の予測機を用いた最適化後の形状とその諸量に対する非線形座屈荷重係数の予測値と真の非線形座屈荷重係数の比であり，1 に近いほど予測精度が高いことを示す．初期形状においては 0.69 と予測値は過小評価されている．

• λp /λ では使用する予測機の予測精度が高いほど誤差が大きくなる傾向を示している．

•使用する予測機の予測精度が高ければ必ずしも真の非線形座屈荷重係数が最大化できるとは限らず，λp /λ の誤差が約 15% 以内である場合に真の非線形座屈荷重係数は最大化できる傾向がある．

今後の展望として，入力変数の見直しや他の予測手法を用いることによって予測機の精度を向上させる事が挙げられる．また，λp /λ の誤差が約 15% 以内になる予測機のパターンを分析する必要がある．

図 6 に予測精度が最大であったエポック数=90 バッチ
サイズ=16 の予測機を用いた最適化後の形状とその諸量を示す．図より λp は 131.76 となっており初期形状に比べて増加している．また総体積も初期形状と同値である
ことから制約を満足した解形状となっていることが確認できる．一方で最適化後の λ は初期形状より低下していることがわかる．また λp /λ も 1.26 となっているため真の値に対して予測値は約 30% 過大評価している事がわ
かる．
図 6，図 7，図 8，図 9 より使用した予測機の予測精度

報告 H60[参考文献]
1） G. Farin, M. S. Kim, and J. Hoschek (eds.). Handbook of Computer
Aided Geometric Design. Elsevier, Amsterdam, 2002.
2） Seth Weidman. Deep Learning from Scratch: Building with python
from First Principles. Oreilly Associates Inc, 2019.09.
3） M. Zhu, F. McKenna, and M. Scott. Openseespy: Python library forthe opensees ﬁnite element framework. SoftwareX, Vol. 7, pp. 6–11,2018.1.
4） R. E. Perez, P. W. Jansen, and J. R. R. A. Martins. pyopt: a pythonbased object-oriented framework for nonlinear constrained optimization. Structural and Multidisciplinary Optimization, Vol. 45, No. 1,pp. 101–118, 2012.1.

－ 446 －