Stress Constraint Topology Optimization using Backpropagation Method in Design Sensitivity Analysis

Min-Geun Kim; Seok-Chan Kim; Jaeseung Kim; Jai-Kyung Lee; Geun-Ho Lee

doi:10.7734/COSEIK.2022.35.6.367

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 December 2022. 367-374
https://doi.org/10.7734/COSEIK.2022.35.6.367

Stress Constraint Topology Optimization using Backpropagation Method in Design Sensitivity Analysis

설계민감도 해석에서 역전파 방법을 사용한 응력제한조건 위상최적설계

Min-Geun Kim¹^†

Seok-Chan Kim²

Jaeseung Kim³

Jai-Kyung Lee¹

Geun-Ho Lee¹

김 민근¹^†

김 석찬²

김 재승³

이 재경¹

이 근호¹

¹Principal Researcher, Department of Smart Industrial Machine Technology, Korea Institute of Machinery & Materials, Daejeon, 34103, Korea

²Postdoctoral Research Fellow, Department of Smart Industrial Machine Technology, Korea Institute of Machinery & Materials, Daejeon, 34103, Korea

³Researcher, Department of Smart Industrial Machine Technology, Korea Institute of Machinery & Materials, Daejeon, 34103, Korea

¹한국기계연구원 스마트산업기계연구실 책임연구원

²한국기계연구원 스마트산업기계연구실 박사후연구원

³한국기계연구원 스마트산업기계연구실 연구원

^{†Corresponding Author}

ABSTRACT

This papter presents the use of the automatic differential method based on the backpropagation method to obtain the design sensitivity and its application to topology optimization considering the stress constraints. Solving topology optimization problems with stress constraints is difficult owing to singularities, the local nature of stress constraints, and nonlinearity with respect to design variables. To solve the singularity problem, the stress relaxation technique is used, and p-norm for stress constraints is applied instead of local stresses for global stress measures. To overcome the nonlinearity of the design variables in stress constraint problems, it is important to analytically obtain the exact design sensitivity. In conventional topology optimization, design sensitivity is obtained efficiently and accurately using the adjoint variable method; however, obtaining the design sensitivity analytically and additionally solving the adjoint equation is difficult. To address this problem, the design sensitivity is obtained using a backpropagation technique that is used to determine optimal weights and biases in the artificial neural network, and it is applied to the topology optimization with the stress constraints. The backpropagation technique is used in automatic differentiation and can simplify the calculation of the design sensitivity for the objectives or constraint functions without complicated analytical derivations. In addition, the backpropagation process is more computationally efficient than solving adjoint equations in sensitivity calculations.

Keywords

backpropagation method

automatic differentiation

design sensitivity

stress constraint topology optimization

본 논문에서는 역전파 방법 기반 자동미분법을 이용하여 설계민감도를 구하고 이를 응력제한조건을 고려한 위상최적설계에 적용하였다. 응력제한조건이 있는 위상최적화문제는 특이점(singularity)과 응력의 국부성(local nature of stress constraint)문제, 그리고 설계 변수에 대한 비선형성의 문제를 포함하고 최적해를 얻기가 매우 힘들다. 특이점 문제를 해결하기 위해서 응력 완화(stress relaxation) 기법을 사용하였고, 응력의 국부성을 해결하기 위해 p-norm을 이용한 전역 응력치를 제한조건에 사용하였다. 설계 변수에 대한 비선형성을 극복하기 위해 해석적인 방법으로 정확한 설계민감도를 구하는 것이 중요하다. 위상최적설계에서 기존에는 보조변수방법(adjoint variable method)을 사용하여 빠르고 정확한 설계민감도를 구했지만, 설계민감도를 해석적으로 구해야 하고, 보조평형방정식을 추가로 풀어야 하는 어려움이 있다. 이를 해결하기 위해서 인공신경망에서 최적 가중치(weights)와 편차(biases)를 구할 때 쓰이는 역전파 기법을 이용하여 설계민감도를 구하고 이를 응력제한조건을 고려한 위상최적설계에 적용하였다. 역전파 기법은 자동미분에 쓰이는 기법으로 목적함수나 제한조건에 대한 설계민감도를 별도의 수식유도 없이 간단하게 구할 수 있는 장점이 있다. 또한, 미분값을 구하는 역전파의 과정이 보조평형방정식을 푸는 것보다 계산시간이 빠르고 해석적 방법으로 구한 설계민감도와 같은 정확도를 보여준다.

키워드

역전파방법

자동미분법

설계민감도 해석

응력제한조건 위상최적설계

MAIN

1. 서 론
2. 본 론
2.1 응력제한조건을 고려한 위상최적설계 정식화
2.2 응력완화와 전역응력치 계산
2.3 설계민감도 해석
2.4 수치예제 결과 및 고찰
3. 결 론

1. 서 론

유한요소법을 이용한 위상최적설계 기법은 기계구조물의 경량화 설계뿐만 아니라, MEMS, 전자기, 유체 분야에서 원하는 거동을 얻기 위한 최적의 레이아웃을 얻는데 많이 사용된다. Bensøe와 Kikuchi(1988)가 균질화 기법을 이용한 위상최적설계 기법을 개발한 이후 밀도법(SIMP : Solid Isotropic Material with Penalizaation; Bensøe, 1995)을 이용한 위상최적설계 기법은 구현의 간단함으로 인해 현재는 여러 상용프로그램에 적용되어 쉽게 활용되고 있다. Sigmund(2001)은 SIMP기반의 컴플라이언스(compliance) 최소화 문제에 대한 Matlab 코드를 99줄로 구현하여, 약간의 코드 수정만으로 누구나 쉽게 위상최적설계를 활용할 수 있도록 하였다. 하지만, 위상최적화 문제에서 체적조건이 주어진 컴플라이언스 최소화 문제와는 다르게 응력제한조건을 고려한 체적최소화 문제는 최적화 과정에서 특이성(singularity), 응력제한조건의 국부성(local nature of stress constraints), 그리고 설계변수에 대한 높은 비선형성을 포함하고 있기 때문에 최적설계를 수행하는데 큰 어려움이 있다.

Le 등(2010)은 이러한 응력제한조건을 고려한 위상최적설계 문제를 해결하기 위한 방안을 제안하였다. 특이성 문제를 해결하기 위해서 재료밀도의 분포와 마찬가지로 밀도값이 0과 1 사이의 갖도록 하면서 체적에 적용되는 재료 밀도의 차수보다 응력에 적용되는 재료 밀도의 차수를 낮추는 응력완화기법(stress relaxation)을 적용하였다. 응력의 국부성을 해결하기 위해서 p-norm으로 정의하는 전역응력치(global stress measure) $σ_{P N}$ 을 정의하여 최대응력 $σ_{\max}$ 를 근사할 수 있게 하였다. Amir (2021)은 응력제한조건을 가지며 컴플라이언스(compliance)가 최소가 되는 문제의 위상최적설계를 수행하였다.

설계변수에 대한 높은 비선형성을 포함하는 문제에 대해서는 정확한 설계민감도를 구하는 것이 매우 중요하며, 정확한 설계민감도를 얻기 위해 해석적인 수식 유도과정을 동반한다. 이러한 해석적인 방법 중 보조변수방법(adjoint variable method)이 가장 많이 사용되고 있다. 이 방법은 응답해석 시스템 방정식과 유사한 보조방정식을 추가로 풀기만 하면 설계변수의 수에 상관없이 빠르고 정확한 설계민감도를 구하는 방법이다. 하지만, 이러한 방법은 복잡한 수식 유도과정을 동반하게 되며, 제한조건의 수만큼 보조방정식을 풀어야 하는 단점이 있다. 특히 컴플라이언스 문제를 제외한 대부분의 위상최적설계문제에서 해석적 설계민감도 유도과정은 상당히 복잡하기 때문에 많은 노력을 필요로 한다. 최근 Chandrasekhar 등(2021)은 JAX(Bradbury et al., 2018) 라이브러리를 사용하여 설계민감도를 자동미분(automatic differentiation)으로 구하고 이를 적용하여 컴플라이언스 최소화, 변위 컨버터 설계, 그리고 재료설계 문제에 대한 위상최적설계를 수행하였다.

본 논문에서는 역전파 방법(backpropagation method)을 이용한 자동미분법을 이용하여 설계민감도를 구하고 이를 응력제한조건을 고려한 위상최적화 문제에 적용하였다. 응력제한조건에서 발생하는 문제를 해결하기 위해 Le 등(2010)가 제안한 응력완화법과 p-norm 전역응력치를 사용하였다. 설계민감도에서 자동미분을 구현하기 위해 Python의 Pytorch 라이브러리(Paszke et al., 2017)를 사용하였다. Pytorch의 환경에서는 설계변수를 torch변수로 선언하고, 목적함수와 제한조건의 평가는 설계변수의 순방향 전파(forward propatation)을 통해 얻는다. 설계민감도는 이러한 목적함수와 제한조건의 평가 후 역전파 과정을 통해 자동미분의 형태로 구해진다. 자동미분을 이용하여 설계민감도 계산하는 것은 보조변수방법에서 보조평형방정식을 푸는 것 대신에 역전파 과정을 진행하는 것이라 할 수 있다. 본 논문에서는 역전파 기반의 자동미분 방법에 대해서 살펴보고 수치예제를 통해 기존의 보조변수방법으로 얻어진 설계민감도와 역전파 기반 자동미분을 이용한 설계민감도의 정확도를 비교하고 제안된 방법의 효율성을 검증하였다.

2. 본 론

이 장에서는 응력제한조건을 고려한 위상최적설계의 정식화와 함께 기존의 설계민감도 해석 방법과 역전파 기반 자동미분법을 이용한 설계민감도 계산과정을 소개하고자 한다.

2.1 응력제한조건을 고려한 위상최적설계 정식화

밀도법 기반의 응력제한조건을 위상최적설계 문제는 다음과 같이 정식화할 수 있다.

(1a)

\overset{\min}{0 < ρ \leq 1} \sum_{e = 1}^{N_{e}} {\bar{ρ}}_{e} v_{e}

(1b)

subject to : σ_{\max} (\bar{ρ}) \leq σ_{a l l o w}

(1c)

with : K (\bar{ρ}) u = f

여기서, $ρ$ 은 설계변수인 재료 밀도 벡터로 유한요소에서 요소당 0에서 1 사이의 값을 갖는다. $v_{e}$ 는 요소 $e$ 의 체적이며, 유한요소의 $K (\bar{ρ})$ 은 강성행렬이다. 주어진 하중벡터 $f$ 에 대한 변위값은 $u$ 이다. $σ_{\max} (\bar{ρ})$ 는 요소별로 계산되는 von Mises 응력 중 최대값을 나타내며, $σ_{a l l o w}$ 은 허용응력을 나타낸다. $\bar{ρ}$ 은 위상최적설계에서 체커보드(chekerboard) 현상을 완화하고 안정된 최적해를 얻기 위한 필터링된 재료 밀도 벡터로 각 요소별 필터링된 재료 밀도 값은 다음과 같이 계산한다.

(2)

{\bar{ρ}}_{e} = \frac{\sum_{i = 1}^{N_{e}} w (x_{i}) v_{i} ρ_{i}}{\sum_{e = 1}^{N_{e}} w (x_{i}) v_{i}}, w (x_{i}) = r_{\min} - ∥ x_{e} - x_{i} ∥

$w (x_{i})$ 은 $i$ 번째 요소의 중심위치 $x_{i}$ 에서의 가중함수로 필터반경 $r_{\min}$ 내의 $e$ 번째 요소와의 거리에 따라 결정되며 항상 양의 값을 갖는다.

2.2 응력완화와 전역응력치 계산

응력완화는 설계공간에서 응력의 특이성(singularity)를 피하기 위해 응력 설계 공간을 완화하는 방법이다. 응력의 특이성은 하중 지지부와 같이 재료가 필요한 부분에 응력을 낮추기 위해 오히려 재료가 사라지는 것과 같이 우리가 원하지 않는 결과를 주기도 한다. 이 문제를 해결하기 위해서 응력영역에서도 SIMP 방법에서 0과 1 사이의 중간영역(intermediate density value) 재료 밀도를 허용하여 설계 공간을 완화하는 방법을 적용한다. 하지만 이러한 중간영역 재료 밀도는 물리적으로 의미가 없기 때문에 재료 밀도를 분포를 극단적으로(black-and- white design; 즉, 재료 밀도가 0 또는 1에 가까운 값을 갖도록)가져갈 필요가 있다. 이를 위해서 재료의 체적이 $V (\bar{ρ_{e}}) = {\bar{ρ}}_{e} V_{0}$ 와 같을 때 재료상수(또는 강성)에 패널티 값 $p$ 을 고려한 가중치 $η_{c} (\bar{ρ_{e}}) = {\bar{ρ}}_{e}^{p}$ 을 다음과 같이 사용한다.

(3)

E ({\bar{ρ}}_{e}) = E_{\min} + η_{c} ({\bar{ρ}}_{e}) (E_{\max} - E_{\min})

여기서, $E_{\min}$ 과 $E_{\max}$ 은 최대 그리고 최소 재료의 영률(Young’s modulus)이며, $V_{0}$ 은 재료가 모두 채워져 있을 때(solid)의 체적이다. 본 논문에서는 각각 $E_{\min} = 10^{- 6}, E_{\max} = 1$ 의 값을 사용하였다. 응력도 설계영역에서 응력의 특이성을 피하고 응력에 의해 black-and-white 재료 분포를 더 잘 이루어질 수 있도록 응력에 대한 패널티 값 $q$ 을 고려한 응력 가중치 $η_{T} ({\bar{ρ}}_{e}) = {\bar{ρ}}_{e}^{q}$ 을 다음과 같이 사용한다.

(4a)

σ_{e}^{V M} (\bar{ρ}) = (E_{\min} + η_{T} ({\bar{ρ}}_{e}) (E_{\max} - E_{\min})) σ_{1, e}^{V M} (\bar{ρ})

(4b)

σ_{1, e}^{V M} (\bar{ρ}) = (σ_{1, e}^{T} V σ_{1, e})^{1 / 2}

(4c)

σ_{1, e} (\bar{ρ}) = C_{1} B_{e} u_{e} (\bar{ρ})

(4d)

V = [\begin{matrix} 1 & - \frac{1}{2} & 0 \\ - \frac{1}{2} & 1 & 0 \\ 0 & 0 & 3 \end{matrix}]

여기서, $C_{i}$ 은 영률이 $E (\bar{ρ_{e}} = 1) = E_{\max}$ 이고 푸아송 비가 $v$ =0.3일 때의 탄성텐서(elasticity tensor)을 의미하고, $B_{e}$ 는 요소 $e$ 의 변형율-변위 관계 행렬이다. 그리고 $u_{e} (\bar{ρ})$ 은 식 (3)를 요소 $e$ 의 DOF(Degree Of Freedom)에서의 변위이다. $σ_{1, e} (\bar{ρ})$ 은 요소 $e$ 의 응력벡터 성분이며, $σ_{1, e}^{V M} (\bar{ρ})$ 은 von Mises 응력으로 식 (4d)의 보조 벡터 행렬 $V$ 값을 이용해 구할 수 있다. 최종적으로 식 (4b)의 $σ_{e}^{V M} (\bar{ρ})$ 은 요소 $e$ 의 von Mises 응력 $σ_{1, e}^{V M} (\bar{ρ})$ (식 (4c))에 밀도 가중치 $η_{T} (\bar{ρ_{e}}) = {\bar{ρ}}_{e}^{q}$ 을 고려한 완화된 응력치이다.

패널티 값 $p$ 와 $q$ 을 이용하여 응력의 특이성을 완화시키는 방법을 pq-relaxation이라 하며(Bruggi, 2008), 본 논문에서는 $p$ =3, $q$ =1/2의 값을 사용하였다. 이 때, 체적과 강성 그리고 응력 함수에 대한 밀도 근사는 Fig. 1과 같다. 0과 1 사이의 밀도값을 갖는 경우 재료 상수는 체적보다 작은 가중치를 갖도록 하며, 응력은 체적 보다 큰 가중치를 갖도록 한다. 이러한 강제(penalization)는 최종 위상최적설계에서 재료 분포를 나타내는 밀도값이 0 또는 1에 가까운 값을 갖도록 강제할 수 있다(Le et al., 2010).

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F1.jpg

Fig. 1.

Interpolation functions for weights

응력제한조건을 고려한 위상최적설계 문제에서는 다음과 같이 모든 요소에서 평가되는 응력 $σ_{e}^{V M} (\bar{ρ})$ 이 허용응력을 만족해야 한다.

(5)

σ_{e}^{V M} (\bar{ρ}) \leq σ_{a l l o w}, e = 1, 2, . . . N_{e}

식 (5)은 다음과 같이 다시 쓸 수 있다.

(6)

σ_{\max} (\bar{ρ}) = \max (σ_{e}^{V M} (\bar{ρ})) \leq σ_{a l l o w}, e = 1, 2, . . . N_{e}

식 (5)을 고려할 경우, 유한요소의 개수 $N_{e}$ 값이 클수록 제한조건수가 많아지기 때문에 보조변수방법을 이용한 설계민감도 해석시에 계산 비용이 극단적으로 증가하는 단점이 있다. 또한, 식 (6)에서 $σ_{\max}$ 은 설계변수 $\bar{ρ}$ 에 대해 미분가능하지 않다. 이러한 문제를 해결하기 위해 다음과 같이 $σ_{\max}$ 을 근사할 수 있는 전역응력치 $σ_{PN}$ 를 정의하여 $N_{e}$ 개의 국부 응력제한조건 $σ_{e}^{V M}$ 을 갖는 문제를 전역응력치 하나를 제한조건으로 갖는 문제로 다음과 같이 바꾸어 최적화를 진행한다(Le et al., 2010).

(7)

σ_{\max} \approx c σ_{P N} \leq σ_{a l l o w},

(8)

σ_{P N} = {(\sum_{e = 1}^{N_{e}} v_{e} {(σ_{e}^{V M})}^{P_{N}})}^{1 / P_{N}}

여기서, $σ_{P N} (\bar{ρ})$ 은 p-norm 또는 KS(Kreisselmeier-Steinhauser) 함수로 정의된 전역응력치이다. 식 (8)의 전역응력치는 $P_{N} \to \infty$ 일 때, $σ_{\max}$ 값으로 수렴함이 알려져 있으나 너무 큰 값은 수치적으로 안정된 값을 주지 못할 수 있으므로(Le et al., 2010), 본 논문에서는 $P_{N}$ =8을 사용하였다.

실제 $σ_{PN}$ 은 실제 물리적인 값이 아니므로 정규화 파라미터 $c$ 를 식 (7)과 같이 적용하여, 이 값이 $σ_{\max}$ 에 근사화하도록 한다. 파라미터 $c$ 값은 최적화를 통해 업데이트 되며, $I$ 번째( $I \geq 1$ ) 최적화 과정에서의 파라미터 $c^{I}$ 값은 다음과 같이 계산된다.

(9)

c^{I} = α^{I} \frac{σ_{\max}^{I - 1}}{σ_{P N}^{I - 1}} + (1 - α^{I}) c^{I - 1}

$α^{I} \in (0, 1]$ 은 $c^{I}$ 과 $c^{I - 1}$ 사이의 변화량을 조절하는 제어 파라미터이다. 본 논문에서는 전체 최적화 과정에서 $α^{I} = 0.5$ , 정규화 파라미터 초기값은 $c^{0} = 1$ 로 설정하였다.

2.3 설계민감도 해석

위상최적설계에서 최적해를 빠르게 얻기 위해서는 정확한 설계민감도 정보가 필요하다. 위상최적설계에서 설계변수는 각 요소에 정의된 밀도값 $ρ_{e}$ 이며, 유한요소의 수만큼 많은 설계변수를 포함하고 있다. 목적함수의 수보다 설계변수가 많을 경우 보조변수방법을 이용한 설계민감도법이 가장 효율적이지만, 식 (7)과 같이 주어진 응력제한조건에 대한 설계민감도는 이 방법을 이용할 경우 해석적으로 복잡한 수식유도 과정을 거쳐야 한다. 인공신경망에서는 역전파를 이용한 자동미분방법을 이용하여 가중치와 편차에 대해 손실함수의 민감도를 구하고 이러한 손실함수가 최소가 되도록 하는 최적의 가중치와 편차를 얻는다. 자동미분 방법은 목적함수를 설계민감도를 통해 표현하기만 하면 자동으로 손쉽게 설계민감도를 구할 수 있는 장점이 있어 보조변수방법의 복잡한 수식유도과정을 피할 수 있다. 본 장에서는 p-norm으로 정의된 응력에 대한 설계민감도를 보조변수방법으로 구하는 과정을 살펴보고 이러한 설계민감도를 역전파 기반 자동미분으로 구현하는 것을 살펴보도록 하겠다.

2.3.1 보조변수방법

식 (7)에서 p-norm 응력 $σ_{PN}$ 의 밀도 $ρ_{e}$ 에 대한 설계민감도는 식 (2)의 필터를 적용할 때 다음과 같이 구할 수 있다.

(10a)

\frac{\partial σ_{P N}}{\partial ρ_{e}} = \sum_{j = 1}^{N_{e}} \frac{\partial σ_{P N}}{\partial {\bar{ρ}}_{j}} \frac{\partial {\bar{ρ}}_{j}}{\partial ρ_{e}} = \sum_{j = 1}^{N_{e}} \frac{H_{e j}}{\sum_{i = 1}^{N_{e}} H_{i e}} \frac{\partial σ_{P N}}{\partial {\bar{ρ}}_{j}}

(10b)

H_{e j} = w (x_{j}) v_{j} = r_{\min} - ∥ x_{e} - x_{j} ∥

여기서, 식 (10a)의 필터된 밀도 ${\bar{ρ}}_{e}$ 에 대한 $σ_{PN}$ 에 설계민감도는 보조변수법을 이용하여 다음과 같이 구할 수 있다.

(11)

\frac{\partial σ_{P N}}{\partial {\bar{ρ}}_{e}} = {(\sum_{i = 1}^{N_{e}} v_{i} (σ_{i}^{V M})^{P_{N}})}^{\frac{1}{P_{N}} - 1} v_{e} (σ_{e}^{V M})^{P_{N} - 1} q (\bar{ρ})_{e}^{q - 1} σ_{1, e}^{V M} + p (\bar{ρ})^{p - 1} (E_{\max} - E_{\min}) λ_{e}^{T} K_{1, e} u_{e}

식 (10b)는 중심 요소에서 거리 내에 있는 요소들 의 거리 가중함수이다. $λ_{e}$ 는 $e$ 번째 요소 DOF에 해당하는 보조변수 벡터이며 $K_{1, e}$ 은 ${\bar{ρ}}_{e} = 1$ 일 때 요소 $e$ 의 강성행렬이다. 이 보조변수 벡터는 다음과 같은 보조변수방정식을 풀어 얻게 된다.

(12a)

K λ = - {(\frac{\partial σ_{P N}}{\partial u})}^{T}

(12b)

{(\frac{\partial σ_{P N}}{\partial u_{e}})}^{T} = \frac{{(\sum_{i = 1}^{N_{e}} v_{i} (σ_{i}^{V M})^{P_{N}})}^{\frac{1}{P_{N}} - 1} v_{e} (σ_{e}^{V M})^{P_{N} - 1} (\bar{ρ_{e}})^{q}}{σ_{1, e}^{V M}} B_{e} {}_{T}C_{1} V σ_{1, e}

여기서 $K$ 는 식 (1c)에 주어진 선형방정식의 강성행렬과 같다. 식 (11)의 우측항은 ${\bar{ρ}}_{e}$ 대한 $σ_{PN}$ 명시적인 의존항(explicit dependence) 암시적인 의존항(implicit dependence)으로 구성된다. 암시적인 의존항은 식 (12)으로 부터 얻어지는 보조변수로 대체되며, 그 결과 식 (11)의 우측항은 설계변수에 대한 명시적인 의존항만으로 표현된다. 즉, 식 (1c)의 시스템 방정식에서 구성된 강성행렬을 이용하여 식 (12)의 보조평형방정식을 추가적으로 풀기만 하면 설계변수의 수와 상관없이 빠르고 효율적으로 설계민감도를 구할 수 있다.

하지만, 식 (11), (12)에서 알 수 있듯이 설계민감도 계산은 목적함수에 따라 복잡한 유도과정을 거쳐야 하며, 목적함수와 제한조건의 수만큼 보조방정식 (12)를 풀어야 하기때문에 이에 따른 계산비용이 증가할 수 있는 문제가 있다. 본 논문에서는 이러한 복잡한 수식의 유도과정을 필요로 하는 설계민감도법의 단점을 극복하기 위해 역전파 기반의 자동미분법을 이용해 설계민감도를 구현하는 것을 다음과 같이 제안하고자 한다.

2.3.2 역전파 방법

자동미분법은 복잡한 계산과정이 필요한 미분 계산에서 매우 강력한 도구이다. 기본적으로 설계 변수에 대한 미분의 대상이 되는 목적함수를 순방향 전파를 이용해 구하고, 목적함수의 미분값은 역적파 과정을 구하게 된다. 순방향 전파의 개념은 평가하고자 하는 함수 $z = f (x)$ 가 $f (x) = h \circ g (x) = h (g (x))$ 와 같이 함수 $h$ 와 $g$ 의 합성함수로 구성될 때 식 (13a)과 같이 $y$ 와 $z$ 를 순차적으로 평가하는 것으로 설명할 수 있다.

(forward propagation)

(13a)

1) y = g (x) \to z = h (y) \to f (x) = z

(backward propagation)

(13b)

2) J_{1} = \frac{\partial h}{\partial y} = \frac{\partial z}{\partial y}, J_{2} = \frac{\partial g}{\partial x} = \frac{\partial y}{\partial x}

(evaluation of derivative)

(13c)

3) J = \frac{\partial f}{\partial x} = \frac{\partial y}{\partial x} \frac{\partial f}{\partial y} = \frac{\partial g}{\partial x} \frac{\partial h}{\partial y} = J_{2} J_{1},

이때, 함수 $f$ 의 $x$ 에 대한 미분 값은 다음과 같이 구할 수 있다.임의의 함수 $f$ 의 미분값을 구하기 위한 대한 과정은 다음과 같다. 1) 순방향 전파를 이용하여 $f$ 의 값을 구한 후, 2) 식 (13b)에서 함수 $h$ 와 $g$ 의 $y$ 와 $x$ 에 대한 자코비안(Jacobian) $J_{1}$ 와 $J_{2}$ 를 각각 구하고, 3) 식 (13c)을 이용하여 최종적인 자코비안 $J$ 를 계산한다. 이 과정을 Fig. 2에 나타내었다. 2)의 과정을 통해 3)의 미분값을 평가하는 과정을 역전파과장이라 한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F2.jpg

Fig. 2.

Process of automatic differentiation

위상최적설계에서 식 (8)과 같이 주어진 p-norm 응력에 대한 역전파 기반 설계민감도 해석과정을 보여주기 위해 다음과 같이 p-norm응력을 $f_{P N}, σ_{e}^{V M}, u$ , 그리고 $\bar{ρ}$ 의 함성함수 형태로 표현하면 다음과 같다.

(14)

σ_{P N} (ρ) = f_{P N} \circ σ_{e}^{V M} \circ u \circ \bar{ρ} (ρ)

식 (14)에서, $f_{PN}$ 은 식 (8)에서 von-Mises 응력 $σ_{e}^{V M}$ 을 이용하여 계산하는 p-norm을 구하는 KS 함수이며, $\bar{ρ} (ρ)$ 은 식 (2)에서 각 요소의 밀도 $ρ_{e}$ 로 이루어진 밀도 벡터 $ρ$ 의 필터링 함수이다. 밀도 벡터 $ρ$ 를 설계변수라 할 때, p-norm 응력함수의 설계민감도는 다음의 식 (15)와 같이 구할 수 있다.

(15)

\frac{\partial σ_{P N} (ρ)}{\partial ρ} = \frac{\partial f_{P N}}{\partial σ_{e}^{V M}} \frac{\partial σ_{e}^{V M}}{\partial u} ∙ \frac{\partial u}{\partial \bar{ρ}} ∙ \frac{\partial \bar{ρ}}{\partial ρ}

이를 식 (13)과 같이 함수 $σ_{P N} (ρ)$ 을 구하기 위한 순방향 전파 과정과 p-norm 응력의 설계변수 $ρ$ 에 대한 설계민감도를 위한 역방향 전파과정은 다음과 같으며 이러한 과정을 Fig. 3에 나타내었다.

(forward propagation)

(16a)

s 1) \bar{ρ} = \bar{ρ} (ρ) \to u = u (\bar{ρ}) \to σ_{e}^{V M} = σ_{e}^{V M} (u) \to f_{P N} = f_{P N} (σ_{e}^{V M})

(backward propagation)

(16b)

s 2) J_{1} = \frac{\partial f_{P N}}{\partial σ_{e}^{V M}}, {\vec{J}}_{2} = \frac{\partial σ_{e}^{V M}}{\partial u}, J_{3} = \frac{\partial u}{\partial \bar{ρ}}, J_{4} = \frac{\partial \bar{ρ}}{\partial ρ}

(evaluation of derivative)

(16c)

s 3) \vec{J} = \frac{\partial σ_{P N} (ρ)}{\partial ρ} = \frac{\partial \bar{ρ}}{\partial ρ} ∙ \frac{\partial u}{\partial \bar{ρ}} ∙ \frac{\partial σ_{e}^{V M}}{\partial u} \frac{\partial f_{P N}}{\partial σ_{e}^{V M}} = J_{4} ∙ J_{3} ∙ J_{2} J_{1}

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F3.jpg

Fig. 3.

Process of automatic differentiation for p-norm stress

$\vec{J}$ 은 $σ_{P N} (ρ)$ 의 설계변수 $ρ$ 에 대한 설계민감도이다. 식 (16)에서 목적함수 $σ_{P N} (ρ)$ 의 평가부터, $J_{1}, {\vec{J}}_{2}, J_{3}$ 그리고 $J_{4}$ 등을 포함한 자동미분값을 구하는 것은 PyTorch(Paszke et al., 2017)와 같은 python 라이브러리를 이용하여, 별도의 복잡한 수식 유도 없이 독립변수(설계변수) $ρ$ 를 정의하고 미분하고자 하는 함수만 정의하면 쉽게 구현할 수 있다. PyTorch를 이용한 설계변수 $ρ$ (rho 변수)의 선언은 CODE SNIPPET 1 같이 할 수 있다.

CODE SNIPPET 1

rho = torch.tensor(rho_init, dtype=torch.float64, requires_grad=Ture)

CODE SNIPPET 2

def stress_PN(rho):

rho_hat = filter(rho) # density filter(식 (2))

K = assemble(rho) # assemble stiffness matrix K

u = solver(K) # solve Ku=f

sress_PN = stressCal(u) # calculate p-norm stress(식 (8))

return stress_PN

CODE SNIPPET 3

sig_PN = stress_PN(rho) # call the function stress_PN and perform the forward propagation(식 (16a))

sig_PN.backward() # perform back propagation(식 (16b))

sig_PN_sensitivity = rho.grad # evaluation of derivative(식 (16c))

rho_init은 변수 rho의 초기값을 의미하며, requires_grad 옵션을 True로 설정하여 rho가 목적함수의 자동미분을 위한 독립변수임을 선언해야 한다. 구하고자 하는 목적함수 $σ_{P N} (ρ)$ 는 CODE SNIPPET 2와 같이 이미 선언된 독립변수 rho의 함수 stress_PN(rho)로 정의한다. CODE SNIPPET2에서 정의된 함수 stress_PN(rho)을 호출하면 식 (16a)의 순방향전파 과정을 통해 p-norm 응력 $σ_{PN}$ 을 sig_PN로 반환할 수 있다.이 때 반환된 p-norm 응력 sig_PN 값을 통해 내장함수인 backward() 함수를 이용하면 식 (16b)와 같이 설계민감도를 위한 역방향 전파과정이 수행되게 된다. 마지막으로 grad 명령을 통해 최종적으로 식 (16c)와 같이 설계민감도를 sig_PN_sensitivity를 얻을 수 있다. 이 전체 과정은 CODE SNIPPET 3에 나타내었다.

설계변수를 구할 때 가장 어려운 부분은 설계변수에 대해 암시적 의존항을 구하는 것이며, 본 논문에서 이러한 항은 식(16b)의 $J_{3} = \partial u / \partial ρ$ 을 구하는 것이다. 보조변수방법에서는 이러한 암시적 의존항을 직접 구하는 대신에 본래 시스템 방정식과 같은 형태의 보조변수방정식을 세우고 여기서 얻어지는 보조변수로 이 항을 대치한다. 자동미분방법에서는 이러한 보조변수방정식을 푸는 것 대신에 역전파 과정을 이용하여 암시적 의존항을 처리하고 설계민감도를 구하게 된다. 보조변수방법에서 목적함수와 제한조건의 수에 따라 보조변수방정식을 풀어야 하는 것과 마찬가지로 자동미분법에서도 목적함수와 제한조건의 개수에 따라 역전파를 진행하게 된다. 따라서, 두 방법의 효율성을 비교하기 위해서는 보조변수방정식을 푸는 속도와 역전파 과정을 진행하는 속도의 차이를 알아보면 된다. 다음 장에서는 수치예제를 통해 두 설계민감도를 비교하고 정확성과 효율성을 검증하도록 하겠다.

2.4 수치예제 결과 및 고찰

2.4.1 설계민감도 비교 검증

역전파 방법을 이용한 설계민감도의 정확성을 보조변수방법과 비교하여 검증하였다. Fig. 4는 설계민감도 검증을 위한 예제이다. 가로 세로의 길이는 각각 32mm과 20mm이고 요소수는 32×20이다. 작용하는 힘은 $f = 1 N$ 이다. 본 예제에서의 초기 밀도값은 모두 $ρ_{e} = 0.5$ 이며, 필터 반경은 $r_{\min} = 2 m m$ 이다. Fig. 5는 Fig. 4 문제에 대한 von Mises 응력값을 식 (4a)를 이용하여 구하고 이를 요소별로 표현한 그림으로 최대값은 하중이 가해지는 부분에서 $7.417 M P a$ 의 값을 갖는다. p-norm 응력값은 $7.488 M P a$ 로 최대 von Mises 응력값과 유사한 수치를 갖는다. 이 값은 식 (7)에서 조정 파라미터 $c$ 값을 이용하여 최대 von Mises 응력에 근사하게 되며, 식 (9)에 의해 $c$ 는 최적화가 진행되어 업데이트가 된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F4.jpg

Fig. 4.

Cantilever beam problem

Table 1은 Fig. 5에 표시된 5개의 지점에서 밀도값에 대한 p-norm 응력의 설계민감도 값을 보조변수방법(AVM)과 역전파 기반 자동미분방법(AD)을 이용하여 구한 결과를 비교한 것이다. 비교 결과 두 방법이 소수 10자리 수까지 일치하는 것을 알 수 있다.

Table 1.

Design sensitivity comparison

설계변수	AVM(a)	AD(b)	(b)/(a) *100(%)
𝜌₂₂	-3.01682307599E-03	-3.01682307599E-03	100.00
𝜌₁₀₀	-1.77091411254E-04	-1.77091411252E-04	100.00
𝜌₁₃₁	-1.69432425445E-04	-1.69432425445E-04	100.00
𝜌₄₇₅	2.24687016468E-05	2.24687016469E-05	100.00
𝜌₆₃₈	-5.65107034471E-01	-5.65107034471E-01	100.00

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F5.jpg

Fig. 5.

von-Mises stress contour

Table 2에서는 보조변수법(AVM)과 역전파 기반의 자동미분법(AD)을 이용한 설계민감도 계산시간을 비교하였고, 그 결과를 Fig. 6에 나타내었다. 계산 시 사용한 컴퓨터 사양은 Inter(R) Core(TM) i9-9900T CPU 2.11GHz, 6GB RAM이다. 요소수가 늘어남에 따라 설계민감도 해석시간이 늘어나며, 두 방법의 차이는 약 20% 정도로 요소수가 늘어남에 따라 그 차이가 조금씩 늘어나는 경향이 있다. 보조변수방법에서 응답해를 구하는 것과 마찬가지로 자동미분법에서는 순방향 전파시에 응답해를 구하기 때문에 두 방법의 시간차이는 보조변수방정식을 푸는 것과 역전파를 진행하는 시간의 차이라 할 수 있고 역전파를 이용한 민감도 계산이 더 효율적임을 알 수 있다.

Table 2.

Design Sensitivity Computation time (sec)

유한요소 수	AVM(a)	AD(b)	(b)/(a) * 100(%)
640	0.1915	0.0847	226.09
2560	1.9428	1.6127	120.47
5760	15.3071	12.6681	120.83
10240	81.94	66.7675	122.72

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F6.jpg

Fig. 6.

von-Mises stress contour

2.4.2 위상최적설계 결과

다음으로 보조변수법과 역전파 기반 자동미분법을 이용한 설계민감도법을 각각 적용하여 응력제한조건을 고려한 위상최적설계를 수행하였다. Fig. 7은 수행하고자 하는 예제의 하중과 경계조건을 나타낸다. 가로 세로의 길이는 각각 240mm과 40mm이고 요소수는 240×40이다. 작용하는 힘은 $f = 8 N$ 으로 응력집중 현상을 막기 위해 중앙 상단 9개의 노드에 대칭적으로 하중을 가하였다. 본 예제에서의 초기 밀도값은 모두 $ρ_{e} = 0.4$ 이며, 필터 반경은 $r_{\min} = 2 m m$ 이다. 목적함수는 식 (1a)의 체적 $V$ 를 최대체적 $V_{\max}$ 로 정규화 정규화한 $V / V_{\max}$ 이며, 응력허용조건 $σ_{a l l o w} = 5 M P a$ (식 (1b))을 만족시키며 정규화한 체적이 최소화가 되는 위상최적설계를 수행하였다. 사용한 최적설계 알고리즘은 MMA(Svanberg, 1987)이다. 최적화된 재료밀도 $ρ$ 의 분포와 von Mises 응력분포는 Fig. 8과 같고 보조변수방법과 역전파 기반 자동미분 방법을 이용한 최적설계 결과는 Table 3과 같다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F7.jpg

Fig. 7.

MBB beam problem

https://cdn.apub.kr/journalsite/sites/jcoseik/2022-035-06/N0040350606/images/Figure_jcoseik_35_06_06_F8.jpg

Fig. 8.

Topology optimization results

Table 3.

Optimization results comparison

	AVM(a)	AD(b)
$V / V_{m a x}$	0.3957	0.3957
$σ_{\max}$ (MPa)	5.0	5.0
Compliance(Nmm)	4159.315	4159.209
# of iterations	828	828
Gray-level indicator(%)	18.798	18.794

여기서, Gray-level indicator(Sigmund, 2007)는 gray scale 정도, 즉 $\bar{ρ} = 0.5$ 에 가까운 밀도값의 비율을 나타내는 것으로 두 위상최적설계의 재료 밀도의 일치 정도를 살펴보기 위해 비교하였다. 목적함수인 정규화된 최적 체적값을 비교한 결과 모든 수치가 일치함을 알 수 있고, 이는 역전파 기반 자동미분법을 이용한 설계민감도의 값이 해석적으로 정확한 보조변수방법과 일치함을 보여준다.

3. 결 론

본 연구에서는 인공신경망에서 쓰이는 역전파 기법을 이용하여 설계민감도를 구하고 이를 응력제한조건을 고려한 위상최적설계에 적용하였다. 응력제한조건에 대한 보조변수방정식을 이용한 설계민감도법은 해석적으로 정확하고 효율적이나 복잡한 수식유도과정을 거쳐야 하는 단점이 있다. 자동미분법은 이러한 단점을 피할 수 있을 뿐만 아니라 설계민감도 계산시간에서도 약 20% 정도 빠른 것을 확인하였다. 또한, 해석적인 방법으로 구한 설계민감도와 같은 정확도를 보여준다. 몇몇의 수치예제를 통해 제시된 자동미분방법이 매우 정확하며 효율적임을 검증하였다.

Acknowledgements

본 연구는 2022년 정부(방위사업청)의 재원으로 국방기술진흥연구소의 지원을 받아 수행된 연구임(No. KRIT-CT-21- 013).

References

Amir, O. (2021) Efficient Stress-Constrained Topology Optimization using Inexact Design Sensitivities, Int. J. Numer. Methods Eng., 122(13), pp.3241~3272. 10.1002/nme.6662

Bensøe, M.P. (1995) Optimization of Structral Topology, Shape and Material, Springer, New York.

Bensøe, M.P., Kikuchi, N. (1988) Generating Optimal Topologies in Optimal Design using a Homogenization Method, Comput. Methods Appl. Mech. & Eng., 71, pp.197~224. 10.1016/0045-7825(88)90086-2

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M.J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., Zhang, Q. (2018) JAX : Composable Transformations of Python+Numpy Program,Version 0.2. 5.

Bruggi, M. (2008) On an Alternative Approach to Stress Constraints Relaxation in Topology Optimization, Struct. Multidiscip. Optim., 36(2), pp.125~141. 10.1007/s00158-007-0203-6

Chandrasekhar, A., Sridhara, S., Suresh, K. (2021) AuTO:a Framework for Automatic Differentiation in Topology Optimization, Struct. Multidiscip. Optim., 64(6), pp.4355~4365. 10.1007/s00158-021-03025-8

Le, C., Norato, J., Bruns, T., Ha, C., Tortorelli, D. (2010) Stress-based Topology Optimization for Continua, Struct. Multidiscip. Optim., 41(4), pp.605~620. 10.1007/s00158-009-0440-y

Paszke, A., Gross, Sam, Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A. (2017) Automatic Differentation in Pytorch, 31^st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA, pp.1~4.

Sigmund, O. (2001) A 99 line Topology Optimization Code Written in Matlab, Struct. Multidiscip. Optim., 21(2), pp. 120~127. 10.1007/s001580050176

Sigmund, O. (2007) Morphology-based Black and White Filters for topology Optimization, Struct. Multidiscip. Optim., 33(4), pp.401~424. 10.1007/s00158-006-0087-x

Svanberg, K. (1987) The Method of Moving Asymptotes-a New Method for Structural Optimization, Int. J. Numer. Methods Eng., 24(2), pp.359~373. 10.1002/nme.1620240207

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Stress Constraint Topology Optimization using Backpropagation Method in Design Sensitivity Analysis

ABSTRACT

MAIN

(1a)

(1b)

(1c)

(2)

(3)

(4a)

(4b)

(4c)

(4d)

Fig. 1.

Interpolation functions for weights

(5)

(6)

(7)

(8)

(9)

(10a)

(10b)

(11)

(12a)

(12b)

(13a)

(13b)

(13c)

Fig. 2.

Process of automatic differentiation

(14)

(15)

(16a)

(16b)

(16c)

Fig. 3.

Process of automatic differentiation for p-norm stress

Fig. 4.

Cantilever beam problem

Table 1.

Design sensitivity comparison

Fig. 5.

von-Mises stress contour

Table 2.

Design Sensitivity Computation time (sec)

Fig. 6.

von-Mises stress contour

Fig. 7.

MBB beam problem

Fig. 8.

Topology optimization results

Table 3.

Optimization results comparison

Acknowledgements

References