88줄로 쓰는 위상최적화 — 코드가 곧 논문인 경우

위상최적화는 정해진 재료량을 설계영역 안에서 어디에 놓을지를 푸는 문제다. 이 논문의 기여는 새로운 이론이 아니라 코드 자체다 — 88줄짜리 MATLAB 프로그램 하나로 유한요소해석, 감도해석, 필터링, 설계변수 갱신을 모두 수행하고, 전체 코드를 논문 부록에 그대로 실었다.

기준점은 같은 연구실에서 나온 Sigmund의 99줄 코드(2001)다. 이 논문은 거기에 밀도 필터를 추가하고, 배열 사전할당과 루프 벡터화로 속도를 끌어올렸다 — 7,500개 요소 벤치마크에서 100배다.

이 글은 Andreassen, Clausen, Schevenels, Lazarov, Sigmund, "Efficient topology optimization in MATLAB using 88 lines of code", Structural and Multidisciplinary Optimization 43: 1–16 (2011)의 본문과 표를 인용한다.

1. 무엇을 푸는가

대상은 MBB 보(메시 60 × 20, 150 × 50, 300 × 100의 세 가지)이고, 체적 제약은 50%다. 재료 모델은 SIMP다 — 요소의 밀도 x가 영률을 정하고, penalization 지수 p는 통상값 3을 쓴다. 지수를 두는 이유는 중간 밀도를 불리하게 만들어 0과 1 쪽으로 몰기 위해서다.

설계변수 갱신은 최적성 조건(optimality criteria) 방법이고, 체적 제약을 만족시키는 라그랑주 승수는 이분법으로 찾는다. 종료 조건은 명확하다 — "연속한 두 설계의 차이(설계변수 기준)의 L∞ 노름이 1% 미만이면 최적화 루프를 종료한다".

함수 하나로 모든 것이 결정된다.

top88(nelx, nely, volfrac, penal, rmin, ft)

인수는 각각 가로·세로 요소 수, 체적분율, penalization 지수, 필터 반경, 필터 종류다. 논문이 "위상최적화 분야에 처음 들어오는 사람의 학습 곡선을 완만하게 하는 실질적 도구"라고 부른 이유가 이 한 줄에 있다.

2. 필터가 없으면 답이 아니다

위상최적화에서 필터는 부가 기능이 아니라 해가 성립하기 위한 조건이다. 필터가 없으면 두 가지가 나타난다.

체커보드 패턴 — 요소가 흑백으로 번갈아 배열되어 수치적으로만 강해 보이는 가짜 강성이 생긴다.

메시 의존성 — 메시를 세밀하게 할수록 다른 위상이 나온다. 같은 문제의 답이 격자에 따라 달라진다.

88줄 코드는 두 가지 필터를 모두 담는다 — 감도 필터와 밀도 필터이며, 인수 ft로 고른다. 필터 반경은 예제에서 설계영역 폭의 0.04배로, 세 메시에서 각각 2.4, 6, 16이다. 결과에 대한 논문의 문장이 필터의 역할을 정확히 말한다.

감도 필터링과 밀도 필터링 모두 체커보드 패턴을 억제하고 메시에 독립적인 설계로 이끈다. 메시를 세밀하게 하는 것은 해의 정밀화로만 이어지고, 다른 위상으로 이어지지 않는다. — Andreassen et al., §3.4

감도 필터링과 밀도 필터링 모두 체커보드 패턴을 억제하고 메시에 독립적인 설계로 이끈다. 메시를 세밀하게 하는 것은 해의 정밀화로만 이어지고, 다른 위상으로 이어지지 않는다. — Andreassen et al., §3.4

필터 반경을 절대 길이로 준다는 점이 핵심이다. 요소 개수가 아니라 물리적 길이로 최소 부재 크기를 정하기 때문에, 격자를 바꿔도 같은 위상이 나온다.

3. 100배

표 1은 감도 필터를 쓴 MBB 보 최적화의 반복당 계산시간(초)이다. 처음 열 번의 반복 평균이며, 하드웨어까지 명시되어 있다 — Lenovo ThinkPad X301, Intel Core2 Duo U9400, 메모리 2 GB, Windows XP SP3(32비트), MATLAB R2010a.

구현60 × 20150 × 50 (7,500요소)300 × 100

99줄 코드 (2001)0.6575.19—

88줄 코드0.150.721.85

conv2 기반 필터 (71줄)0.130.691.98

PDE 기반 필터0.130.782.18

150 × 50 메시에서 75.19초가 0.72초가 된다. 논문의 표현으로 "100배의 속도 개선"이다. 300 × 100 메시에서 99줄 코드는 아예 시험하지 않았다.

여기서 중요한 것은 알고리즘이 같다는 사실이다. 같은 SIMP, 같은 최적성 조건 갱신, 같은 필터다. 두 자릿수 차이는 배열 사전할당과 루프 벡터화에서 나왔다. 논문이 강조하듯 이 개선은 "가독성을 희생하지 않고" 이뤄졌고, 코드는 오히려 11줄 짧아졌다.

4. 필터를 구현하는 세 가지 방법

논문 후반부는 필터 구현의 대안을 비교한다.

구현성질

기본(88줄)가중치 행렬을 직접 만든다. 비용은 2차원에서 rmin²에, 3차원에서 rmin³에 비례

conv2 기반 (71줄)MATLAB 내장 합성곱을 쓴다. 수학적으로 동등하고 코드가 71줄로 줄지만, 필터링 절차가 가려지고 정규 메시에만 적용된다

PDE 기반헬름홀츠형 편미분방정식으로 필터링해 유한요소 솔버를 재사용한다. 3D에서 크게 빨라지고 병렬 구현이 단순해진다. 필터 반경과의 관계는 R = rmin/(2√3)

PDE 필터의 비용은 필터 반경과 무관하다는 점이 직접 구현과의 결정적 차이다. 반경이 큰 3차원 문제에서 직접 필터는 rmin³로 커지지만 PDE 필터는 그렇지 않다.

다만 2차원의 작은 문제에서는 이 이점이 드러나지 않는다. 밀도 필터를 쓴 표 2에서 300 × 100 메시의 반복당 시간은 88줄 5.67초, conv2 3.30초, PDE 10.08초다. 논문은 그 이유도 적는다 — 밀도 필터는 라그랑주 승수를 찾는 이분법의 매 반복마다 적용되고, PDE 필터에서는 그때마다 후진대입이 상대적으로 비싸기 때문이다.

5. 흑백으로 만드는 필터

밀도 필터는 경계를 흐리게 만든다. 그래서 논문은 Heaviside 투영 필터를 확장으로 제시한다 — 밀도 필터에 계단함수를 붙여 중간 밀도를 0이나 1로 밀어내는 방식이다. 목적은 두 가지다 — 최소 길이 스케일의 확보와 흑백 해의 획득이다.

밀도 필터가 제조 가능한 최소 치수를 보장하고, 투영이 회색 영역을 없앤다. 두 요구는 원래 서로 다른 것인데, 실무에서는 둘 다 필요하다 — 도면으로 옮길 수 없는 회색 밀도장은 설계가 아니기 때문이다.

6. 정리

이 논문의 기여는 코드다 — 88줄 MATLAB 프로그램 전체가 부록에 실려 있고 topopt.dtu.dk에서 내려받을 수 있다.

SIMP 지수는 통상값 p = 3, 예제의 체적 제약은 50%, 갱신은 최적성 조건 방법이다.

수렴은 설계변수 변화의 L∞ 노름이 1% 미만일 때 종료한다.

필터는 선택이 아니라 조건 — 없으면 체커보드와 메시 의존성이 생긴다.

필터 반경은 설계영역 폭의 0.04배로 주어 세 메시에서 2.4 · 6 · 16이 되고, 메시를 세밀하게 해도 위상은 바뀌지 않고 해상도만 올라간다.

150 × 50(7,500요소)에서 99줄 75.19초 대 88줄 0.72초 — 알고리즘이 아니라 사전할당과 벡터화가 만든 100배다.

conv2 구현은 71줄로 줄지만 정규 메시에 한정되고, 필터링 절차가 코드에서 보이지 않게 된다.

PDE(헬름홀츠) 필터의 비용은 필터 반경과 무관하며 R = rmin/(2√3)로 대응된다. 3D와 병렬화에서 유리하다.

Heaviside 투영 필터는 최소 길이 스케일과 흑백 해를 동시에 겨냥한 확장이다.

구조 실무에서 이 논문을 다시 볼 자리는 두 곳이다. 하나는 필터 반경을 무엇으로 정할 것인가이다 — 이 값이 최소 부재 치수를 결정하므로, 제작 공정이 허용하는 최소 두께에서 거꾸로 잡아야 한다. 다른 하나는 격자를 바꿔도 같은 답이 나오는지를 확인하는 습관이다. 이 논문의 예제가 세 가지 메시로 같은 위상을 보이는 것은 자랑이 아니라 필터가 제대로 걸렸다는 증거다. 위상이 메시에 따라 바뀌면 그 결과는 설계가 아니라 격자의 산물이다.