ガウス過程回帰を用いたラチスシェルの形状最適化に関する基礎的研究
Fundamental Study on Shape Optimization of Reticulated Shell
Using Gaussian Process Regression○松原 大悟*1，瀧内 雄二*2，中澤 祥二*3
Daigo MATSUBARA＊1, Yuji TAKIUCHI＊2 and Shoji NAKAZAWA＊3*1 豊橋技術科学大学

大学院生

Graduate Student, Toyohashi University of Technology*2 豊橋技術科学大学

助教

博士(工学)

Assistant Professor, Toyohashi University of Technology, Dr.Eng.
*3 豊橋技術科学大学教授

博士(工学)

Professor, Toyohashi University of Technology, Dr.Eng.
キーワード：ガウス過程回帰; 構造最適化; ラチスシェル; ベイズ最適化; NURBS 曲線; 機械学習
Keywords: Gaussian process regression; structural optimization; reticulated shell; bayesian optimization;
NURBS curve; machine learning.

1. はじめにそこで本研究では，少ない解析回数で最適解を探索す

近年，施工技術の発達により，意匠性に富んだ複雑な形るベイズ最適化と，大量のデータから近似形状（応答曲状のラチスシェルが多く計画されている。ラチスシェルに面）を生成して最適解を探索する手法を構造最適化問題

は，外力を主に面内力で支持部まで伝達する力学特性があに適用し，その有用性を検討することを目的とする。本り，構造的に合理的な形状の探索が重要となる。これまで，研究では，一例としてラチスシェルの形状を設計変数とシェル構造の最適化に関する様々な研究 1)2)が行われておし，固定荷重によるひずみエネルギの最小化を目的関数

り，力学的な合理性のみならず，意匠性や施工性を目的関とする。

数に取った最適化手法も提案されている。
力学的な合理性を測る指標に注目すると，ひずみエネル

2. ベイズ最適化の概要

ギがよく利用されている。一方で，ラチスシェルは座屈荷

2.1. ベイズ最適化の流れ重が設計上重要となるため，非線形座屈荷重の最大化や弾ベイズ最適化 5)は少ないデータから目的関数 f を近似す塑性座屈荷重の推定値の最大化も提案されている。大規模る関数 f *をモデル化し効率的に最適値を求める手法であなラチスシェルの構造設計では弾塑性座屈解析や弾塑性る。ベイズ最適化では以下の手続きを行う。

地震応答解析を利用して安全性を確認する。このため，非step1 : データを基に f をモデル化し近似関数 f *を作成す

線形解析の結果を目的関数に利用することができれば有る。

用であると考えられる。しかし，非線形解析は線形解析にstep2 : f *を基に，次に探索する点 x*を計算する。

比べて解析時間がかかる場合が多い。最適化計算では目的step3 : f (x*)を計算し，f *を更新する。

関数の値を繰り返し計算するため，1 回の目的関数の計算

step4 : step2，3 を所定の回数繰り返し，終了する。に時間がかかる場合は，最適な形態を獲得するまでに膨大

本研究では目的関数のモデル化には次節で説明するガな時間を要するという問題点がある。このため，非線形解ウス過程回帰を用いる。ガウス過程回帰では任意の x に対析の結果を利用する最適化は，少ない解析回数で行うことする目的関数の期待値 f *と標準偏差 σ*が計算される。このが望まれる。

性質を利用し，2.3 節に示す獲得関数が最大となる点を次一方で，近年では機械学習の研究が建築分野でも進めらに探索する候補点として決定する。

れ大量のデータを有した研究もある。例えば Zhu ら 3)は大

ここでは例として，図 1 の 1 変数関数 f をベイズ最適化

量の座屈解析の結果を集めて，機械学習により弾塑性座屈により最小化する問題を考える。図 1(a)の●印は f により荷重の推定を行っている。また，瀧内ら 4)はシェルの振動求めた観測値である。x = 0 から 0.5 付近は σ*が小さく，回性状のデータベースを構築し分析している。このようなデ帰関数の精度が良い。一方で x < -1 や x > 1 の範囲は σ*がータベースが社会的に整備されれば，最適化を実施する上大きく，十分に回帰が行えていない。ベイズ最適化では f *で，データベースを利用した最適化手法も考えられる。の不確実性を考慮して獲得関数を計算し，獲得関数が最大

日本建築学会情報システム技術委員会
第 44 回情報･システム･利用･技術シンポジウム論文集， 491-494 ， 2021 年 12 月，京都
Proceedings of the 44th Symposium on Computer Technology of Information ，Systems and Applications ， AIJ ， 491-494 ， Dec. ， 2021 ， Kyoto－ 491 －となる図 1(a)の□の点を次のステップの観測点 x*とする。fあり，f +はその解析 step までの目的関数の最小値である。

(x*)が計算され，f *, σ*が更新される。(a)と(b)を比較するとベイズ最適化はガウス過程により f *の事後確率の期待値x = 1 付近の σ*が小さくなることがわかる。同様のことをと標準偏差を考慮し，次のステップで計算する点を探索す

繰り返し，f *の精度を上昇させながら最良解を探索する。

る。

f,f*

f,f*

aEI

3. 解析条件
3.1. 解析モデル

aEI

解析モデルを図 2(a)に示す。解析モデルはスパン 50m x50m の矩形平面を有するラチスシェルであり，支持条件は4 隅の 3 節点をピン支持とする。部材は鋼管を想定し，接f (x*)xx(a)

合部は剛接合とする。
部材のヤング係数 E を 205,000N/mm2，*

f *±σ*f*f

図1x(b)
獲得関数

降伏応力度 σy=235N/mm2 とする。X 軸，Y 軸と平行の部材を格子材，その他の部材を斜材とし，断面諸元を表 1 にまベイズ最適化の流れ

とめる。
荷重条件は比較的扁平なシェルを想定し，水平投影面積2.2. ガウス過程回帰による目的関数のモデル化
ガウス過程回帰は M 個のデータ(xi, yi)，(i = 1,2,…,M)が与えられたとき，新しい M+1 個目のデータ xM+1 に対して，事前分布から事後確率を求めるベイズ確率を利用してyM+1 を求める手法である。X = {x1, x2,…, xM}T に対し，y ={y1, y2,…, yM}T の同時分布 p(y|X)が式(1)のように多変量ガ当たり 1.0kN/m2 の荷重を固定荷重として作用させる。
本研究ではより少ない変数で多様なシェル形状を表現するために，3x3 次の NURBS 曲面を利用する。NURBS 曲面上の任意の点 C(u, v)は二つの媒介変数 u, v(0≦u，v≦1)により式(7)で表される。
nウス分布に従うと仮定する。

C ( u, v ) =p( y X ) = N (  , K )(1)

ここで N(μ, K)は平均ベクトル μ，共分散行列 K の多変量ガウス分布を表す。K の i 行 j 列要素はカーネル関数 k(xi,xj)を用いて計算する。任意の xM+1 に対する yM+1 の値は事後確率として式(2), (3)により，その期待値 f *と分散 σ2*を求めることができる。k*は X と xM+1 から計算される K の i
行 M＋1 列要素であり，M 次のベクトルである。

f * = E  yM +1 xM +1 , ( X , y ) = k *T K −1 y

2*=  ( yM +1 xM +1 , ( X , y ) )(2)

2= k ( xM +1 , xM +1 ) − k K k*T−1(3)* ( N (u )M (v)w P )mi =0 j =0nmi,kj ,lijij(7) ( Ni , k (u )M j ,l (v)wij )i =0 j =0

式中の Ni, k(u)，Mj, l(v)は NURBS 基底関数である。Pij = {PXij,PYij, PZij}T は制御点座標，wij は重み係数であり本研究ではすべての wij は 1.0 とする。
形状最適化解析では制御点 Pij の高さ方向成分 PZij を設計変数とする。なお，形状の対称性を考慮し，図 2(b)の●印の 9 個の PZij を設計変数として考慮する。形状最適化解
析においてライズが極めて高くなることを避けるために，
本研究では設計変数の範囲は 0  PZij  10 m とした。
Yk ( x1 , x M +1 ) , k ( x2 , x M +1 ) ,k =
..., k ( x M , x M +1 ) T*

2.3. 獲得関数

X P03 P13 P23(1)50m

P02 P12 P
22(2)本研究では式(4)~(6)に示すベイズ最適化でよく利用される，期待改善度 EI，改善確率 PI，信頼性下限関数 LCB

P33

P01 P11(3)の 3 種類の獲得関数を用いる。

()+**
EI: aEI = f − f ( x )  ( Z ) +  ( x ) ( Z )Z=f + − f * ( x)

 ( x)

図2

解析モデル

表1

部材断面諸元*

 f + − f * ( x) PI: aPI = P ( f ( x )  f ) =  *  ( x) +(5)log R(6)
R:データ数R
ここで Φ,  は正規分布の累積分布関数と確率密度関数でLCB: aLCB = f * ( x ) −  * ( x ) , =報告 H854.2m 8.3m 12.5m(b)制御点の位置

50m(a)平面形状

(4)

部材

λdt

DAIy

格子材

40353.6

7.9361.587.75

13710

斜材

70285.79.0294.780.788243

λ:細長比 d:管厚中心径[mm] t:管厚[mm] D:外径[mm]A:断面積[cm2] Iy:断面 2 次モーメント[cm4]－ 492 －

3.2. 最適化解析の条件

図 3 にベイズ最適化(LCB)における目的関数の推移を示

本研究は，ベイズ最適化と応答曲面を用いた最適化によto minimize f ( x ) =

1 T
D SD2

す。
グラフの横軸は解析 step である。
250step までに log10( f )が 7.0 程度に収束していることがわかる。って式(8)に示すひずみエネルギを最小化する。

表2(8)

ここで，D は変位ベクトル，S は剛性マトリクスである。

手法

ひずみエネルギの最小化では勾配法を用いた研究事例が数多く報告されている。本研究は基礎的段階として，ベイ

ひずみエネルギの値[kNm]
AB

CEIPILCB

乱数 1

20.08115.936

9.547

乱数 2

31.84221.3638.911

ズ最適化と応答曲面を用いた最適化の有効性を確認する

乱数 3

25.41922.5219.283

ために，勾配法の 1 種である L-BFGS-B 法と結果を比較す

乱数 4

22.14720.04010.430

乱数 517.339

17.3398.362

μ23.365

19.449.307-7.346

σ4.9932.459

0.688-0.143る。なお，ベイズ最適化と応答曲面を用いた最適化は目的関数の勾配を必要としない。非線形解析の結果など高コス

7.565
7.240
9.193

7.469
7.216
7.240

トな目的関数を用いた最適化は今後の課題とする。
ガウス過程回帰は数値解析ライブラリの Scikit-learn を8.5
8.0
7.5
7.06.5用いて行う。なお，ガウス過程回帰ではひずみエネルギの対数を取って正規化したデータを利用する。また，カーネル関数は予備解析の結果より精度が最もよかったMatern5/2 カーネル 5)を用いる。

イズ最適化は与えられたデータによって回帰関数の形が変わるため，初期値を変えて各 5 回解析を行う。初期値は乱数により 10 個与え，解析 step 数は 1000 とする。
手法 B（応答曲面を用いた最適化）：一様分布の乱数より生成した 5000 個の構造データから応答曲面を回帰し，回帰値の最小値を求める。この方法はベイズ最適化とは異なり回帰関数の更新を行わない。また，蓄積された構造物のデータを有効に利用できる手法である。
手法 C（L-BFGS-B 法）：手法 A と手法 B の解析結果を評価するため，変数の上下限の制約条件を考慮できる準ニュートン法である L-BFGS-B 法を利用する。計算には数値解析ライブラリの Scipy を用いる。探索開始点はベイズ最適化の初期値と同様に乱数により与える。収束判定定数はstep0
図3

手法 A（ベイズ最適化）：ベイズ最適化では，獲得関数として EI，PI，LCB の 3 種類を考慮しそれぞれ解析を行う。ベ

log10(f )5001000目的関数の推移(LCB, 乱数 5)4.2. 最良形状各手法で得られた最良形状を図 4 に示す。
どの手法でも，中心部が高く，外周部と斜材方向にアーチ構造が形成される形状が探索される。図 2(a)に示す断面での，各手法の最良形状の節点の高さ分布を図 5 に示す。断面 1 ではシェル
中央の X = 0 付近で各形状の高さの違いは小さいが，外周部に向かうにつれて節点高さの差が大きくなる。断面 2 ではどの形状も同様の高さ分布である。断面 3 での最高高さは，L-BFGS-B 法によって求まった形状が最も高く，次にベイズ最適化(LCB)，応答曲面を用いた最適化によって求められた形状の順である。この順番はひずみエネルギが小さい順番であり，外周部の高さの差がひずみエネルギの差に影響していると考えられる。

ε=1.0×10-3 とする。
4. 解析結果
4.1. ひずみエネルギの値
表 2 に各手法で得られたひずみエネルギの最良解の結果と 5 回の解析結果の平均値 μ と標準偏差 σ を示す。ベイズ最適化(LCB)と応答曲面を用いた最適化によって得られたひずみエネルギの値は，L-BFGS-B 法によって得られた値よりも平均値で 1.3 倍程度大きい。ベイズ最適化での獲得関数ごとの解析結果の違いに注目すると，どの初期値でもLCB が最小値を探索することがわかる。また，標準偏差も(a) ベイズ最適化 (LCB)(b) 応答曲面を用いた最適化ZYX(c) L-BFGS-B 法

図4

各手法の最良形状

LCB の結果が最も小さい。ベイズ最適化(LCB)と応答曲面

図 6 に最良形状の各断面での固定荷重作用時における軸を用いた最適化の結果を比較すると，平均値では応答曲面応力分布を示す。3 個の形状の軸応力分布は似た性状となを用いた最適化が優れているが，最良解はベイズ最適化っていることがわかる。また，どの形状も曲げモーメント(LCB)のほうが優れている。の発生が小さいことを確認している。

報告 H85－ 493 －108X[m]6-30121086420-20-104 1(a) 断面h[m]-3010010210ベイズ最適化(LCB)
系列380系列18 -50応答曲面を用いた最適化
系列166系列2L-BFGS-B
4系列24X[m]2-20 2 -100系列3系列30 3(c) 断面0 -50
図-505 各最良形状の節点高さ
0σN[N/mm2]0X[m]-10σN[N/mm2]

-30

00

DBO[m]
2.0

4.0

6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0(a) ベイズ最適化(LCB)

0.2
0.0
0.0

DRS[m]
2.0

4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0(b) 応答曲面を用いた最適化

図7

最良解からの距離

これより，ベイズ最適化は応答曲面を用いた最適化よりも最良解に近い範囲のデータが多く，効率よく最良解を探索していることがわかる。

本研究ではベイズ最適化と応答曲面を用いた最適化をラチスシェルの形状最適化問題に適用した。以下に得られた結論を示す。

X[m]-25 -20 -15 -10 -5(a) 断面4 1
σN[N/mm2]0.00.0-20-40 系列1-25 -20 -15 6-10 -50.40.2

5. まとめ-3080.4-1010-20データ数の密度12h[m]
108642 系列1X[m]0-30-20-100系列2(b) 断面 2h[m]データ数の密度

121086420(1) ベイズ最適化と応答曲面を用いた最適化によるひずみ0(b) 断面 2系列2エネルギの最小値は L-BFGS-B 法に比べて 1.2~1.4 倍程度となった。得られた形状は概ねすべて同じ形状となり，ベ10210系列3ベイズ最適化(LCB)-1080
系列18 -50応答曲面を用いた最適化
系列1-2066-30L-BFGS-B 法
系列24系列24 X[m]-402-25 -20 -15 2-10 -5 0系列3系列30 3(c) 断面0 -50
図 6 各最良形状の軸応力-50

イズ最適化と応答曲面を用いた最適化の有効性が確認できる。
(2) ベイズ最適化は獲得関数や初期値によってばらつきがあるが，応答曲面を直接最小化するより，少ない解析回数で，ひずみエネルギが小さい結果が得られる。
今後は解析 step や初期データの個数と最適化の精度の関係性を分析することに加え，弾塑性地震応答解析など 1回の計算コストが大きい問題に対しても検討を行う。

4.3. 最良解とデータの関係
ベイズ最適化による最良解を xBO，ベイズ最適化で探索されたデータを ixBO，応答曲面を用いた最適化による最良解を xRS，構造データの j 番目のデータを jx とし式(9)のように書く。最良解とすべてのデータの距離 iDBO, jDRS を式(10), (11)のように計算する。
xBO = PBO1z , PBO 2 z ,..., PBO 9 z  , i xBO =  i P1z , i P2 z ,..., i P9 z T

TxRS = PRS1z , PRS 2 z ,..., PRS 9 z  , j x =  j P1z , j P2 z ,..., j P9 z TTi

DBO ( xBO , i xBO )=j( PBO1z − i P1z ) + ( PBO 2 z − i P2 z ) +22+ ( PBO 9 z − i P9 z )DRS ( xRS , j x )=(P − P ) + (P − P ) +21z

j 1z

22z

j2z+ ( P9 z − j P9 z )22(9)(10)(11)

図 7 に DBO，DRS の分布のヒストグラムを示す。DRS は構造データを一様分布の乱数により生成しているため，正規分布になっている。DBO は最良解から 0.6m-0.8m の範囲が最［参考文献］
1） 浜田英明, 大森博司：設計者の選好と力学的合理性を勘案した自由曲面シェル構造の構造形態創生法の提案 その 1 多目的遺伝的アルゴリズムによる発見的手法, 日本建築学会構造系論文集 Vol.71, No.609, pp.105-111, 2006.112） 山本憲司, 皆川洋一, 大森博司：座屈荷重を目的関数とする空間構造の形状最適化に関する研究, 日本建築学会構造系論文集, No.564, pp.95-102, 2003.23） S. Zhu, M. Ohsaki, and X. Guo. Prediction of non-linear bucklingload of imperfect reticulated shell using modified consistentimperfection and machine learning. Engineering Structures,
Volume 226, 2021.14） 瀧内雄二, 中澤祥二, 東山侑生, 加藤史郎：振動データベースの分析に基づいた単層ラチスシェルの静的地震荷重に関する基礎的研究 その 1 振動データベースの構築, 日本建築学会大会学術講演梗概集, pp.653-654, 2021.95） J. Snoek, H. Larochelle, and R. Adams. Practical Bayesian
optimization of machine learning algorithms. Advances in Neural
Information Processing Systems 25, pp.2960–2968, 2012

も多く，DRS は最良解から 1.2m-1.4m の範囲が最も多い。

報告 H85－ 494 －