Prerequisites

본 포스트를 잘 이해하기 위해선 아래의 내용에 대해 알고 오는 것이 좋습니다.

• 방향장과 오일러 방법
• 위상 평면
• ODE 암시적 솔버 (implicit solvers)
• 라그랑주 승수에 대한 이해 (SVM 편의 챕터 2)

들어가면서

PSPICE나 Simscape 같은 물리 시뮬레이션 소프트웨어를 사용해본 적이 있는가? 공대생이라면 학부 시절 한번쯤은 써봤을 수도 있는 소프트웨어가 아닐까 한다만, 생긴건 아래와 같이 블록으로 표현되는 소자들을 선으로 연결해주면 물리 시뮬레이션을 수행해주는 기능을 한다. 아래 그림은 Simscape의 예시인데, 왼쪽에 보이는 Mass-Spring-Damper 모델을 가운데에 보이는 모델처럼 표현할 수 있다.

그림 1. Mass-Spring-Damper 문제를 Simscape 모델로 푸는 예시

이런 물리 모델링/시뮬레이션은 볼수록 신기하다. 그저 블록을 연결할 뿐인데 물리 모델링을 해준다니? Mass-Spring-Damper 모델이야 학부 저학년때도 배우는 단순한 모델이니까 그렇다쳐도, 꽤 복잡한 물리 현상들까지도 척척 풀어주는 걸 보면 이 개발자들은 신(God)인가 싶다.

알아보니 이런 물리 모델링/시뮬레이션 소프트웨어의 뒷면에는 대수방정식(Differential Algebraic Systems of Equation)이라는 수학적 이론이 포함되어 있다고 했다. 몇 년동안 이게 뭔지 이해하지 못하다가 드디어 약간 알것만 같아 휘뚜루마뚜루 정리해본다.

ODE vs. DAE

비유적으로 설명하면 ODE는 방향장 혹은 위상 평면을 따라 달리는 자유로운 강아지라면, DAE는 제약 사항을 가지고 움직이는 기찻길 위의 기차같은 것이다. (이 비유는 “제약사항”에 주목해야 함.) 즉, DAE가 말하는 것은

• D (Differential, 미분): “어떻게 변화하는가?” 시간에 따른 상태의 변화량($\dot{x}$)을 나타낸다. 가속도, 속도, 전류의 변화 등 시스템의 동특성(Dynamics)을 담당한다.
• A (Algebraic, 대수): “어떤 규칙을 따르는가?” 미분항이 없는 방정식으로, 시스템이 반드시 지켜야 하는 제약 조건(Constraints)이나 보존 법칙을 의미한다. “실의 길이는 일정해야 한다”, “한 노드에서 전류의 합은 0이다” 같은 물리적 약속이다.

결국 DAE는 “물리적인 변화(D)를 하되, 정해진 규칙(A)을 반드시 지키며 움직여라”라는 명령을 담은 방정식 체계라고 할 수 있다.

좋은 설명이지만 정확한 이해를 위해 우선 두 방정식의 정의부터 다시 한번 출발해보자. 상미분방정식(Ordinary Differential Equation, ODE)과 대수방정식 (Differential Algebraic Systems of Equations, DAE)의 수학적 정의부터 한번 파헤쳐보자.

ODE (상미분 방정식)

수학적으로 Explicit(명시적) 형태를 띈다. $n$차 ODE는 항상 1차 연립 방정식으로 변환 가능하며, 벡터 필드 위에서 상태가 어떻게 변할지 직접 가르쳐준다.

• 특징: 모든 변수의 변화율($\dot{\mathbf{x}}$)이 정의되어 있다.
• 기하학적 의미: 상태 공간의 모든 지점에서 ‘어디로 가야 할지’ 화살표가 그려져 있는 상태이다.

이는 방향장과 오일러 방법 편이나 위상 평면 편에서도 다룬 것 처럼, 미분 방정식을 푼다는 것은 (쉽게 설명하자면) 방향장 혹은 위상 평면 위에 한 점이 있고, 그 점이 미분 계수로 표현되는 주변과의 관계를 통해 다음 점이 어디로 이동할지를 결정하는 것과 같다는 것이다.

그림 2. (2, -1)이라는 점으로부터 시작해 위상평면 상의 기울기를 따라 스텝사이즈 $\Delta t = 0.5$ 만큼씩 이동하여 나타내는 솔루션 (출처: 위상평면 편)

DAE (미분 대수 방정식)

미분 방정식에 Implicit(암시적) 제약 조건이 추가된 형태이다. 가장 일반적인 형태는 아래와 같은 형태이다.

• 특징: $\mathbf{x}$는 미분 변수(Differential variable)이고, $\mathbf{z}$는 대수 변수(Algebraic variable)이다. $\mathbf{z}$는 미분식이 없어서 어떻게 변할지 바로 알 수 없고, 아래의 제약 조건 $\mathbf{g}=0$을 만족하는 값으로 결정되어야 한다.
• 기하학적 의미: 상태 공간 전체를 움직이는 게 아니라, $\mathbf{g}(\mathbf{x}, \mathbf{z})=0$이라는 매니폴드(Manifold, 곡면) 위에서만 움직여야 한다.

처음보면 좀 당황스러울 수 있지만 차근히 생각해보자. 우선 식 (2)의 첫번째 식 $\dot{\mathbf{x}} = \mathbf{f}(t, \mathbf{x}, \mathbf{z})$은 아무튼간에 변수 $z$가 추가된 미분방정식이라 생각할 수 있다. 그리고 두 번째 식 $\mathbf{0} = \mathbf{g}(t, \mathbf{x}, \mathbf{z})$은 $x$를 미분한 미분계수가 없는 식이라는 것이 핵심이다. 미분계수가 없는 식을 “Algebraic Equation”이라고 부르는 것이고, 이것을 “대수 규칙”이라고 볼 수 있다. 물리적인 제약사항을 표시하기 위해 존재한다.

이런게 왜 필요할까? 예를 들어서 생각해보자.

RC 회로의 ODE 풀이와 DAE 풀이

RC 회로를 푸는 것을 보면 ODE와 DAE의 차이를 쉽게 알 수 있다. 뿐만 아니라 왜 DAE가 SPICE 소프트웨어나 Simscape 같은 블록 조합의 모델링에 유용한지도 이해할 수 있다.

여기서는 저항 하나에 캐패시터 하나가 직렬로 나열된 회로를 생각하고자 한다.

ODE는 흔히 사용하는 explicit solver인 ode45를 사용하고 DAE는 제약 조건(미분계수가 없는 식)을 포함하여 ode15s로 풀어볼 수 있다.

ODE 풀이

ODE의 풀이는 구성요소 전체를 아우르는 모델을 수식으로 합치는 과정을 포함한다. 키르히호프 전압 법칙을 써서 $V_{in}=IR+\frac{Q}{C}$를 만들고, $I=\frac{dQ}{dt}$임을 이용하면 하나의 미분방정식으로 모델링할 수 있다.

%% --- [Case 1] ODE 방식 (커패시터 전압 Vc 기반) ---
% C * dVc/dt = I = (Vin - Vc) / R
% dVc/dt = (Vin - Vc) / (R*C)
f_ode = @(t, Vc) (Vin - Vc) / (R * C);
[t_ode, V_ode] = ode45(f_ode, tspan, 0); % 초기 전압 0V

DAE 풀이

DAE는 각 소자들(저항, 캐패시터)의 물리 법칙을 나열하고, 제약사항인 대수 법칙을 작성해준다. 저항의 법칙은 $V=IR$이고 캐패시터는 $I=C\frac{dV_C}{dt}$이다. 두 소자 사이를 연결하는 법칙은 키르히호프의 전압 법칙 $V_{in}-V_R-V_C=0$이다. 즉, SPICE나 Simscape에서 블록을 연결하는 것 처럼 각 소자들의 지배방정식을 나열해주고 전체를 관통하는 대수방정식으로 (domain의) 제약사항을 기술해준다. 대수방정식은 SPICE나 Simscape에서는 물리 모델링 블록들을 선으로 연결하는 행위에 해당한다고 볼 수 있다.

%% --- [Case 2] DAE 방식 (소자 식 나열: 비인과적 모델링) ---
% 상태 변수 Y = [Vr; Vc; I]
% 1. Vr - I*R = 0          (저항의 법칙 - 대수)
% 2. C * Vc' = I           (커패시터의 법칙 - 미분)
% 3. Vin - Vr - Vc = 0     (KVL 연결 법칙 - 대수)

% Mass Matrix: 미분항이 있는 Vc(2번 변수)만 1, 나머지는 0
M = [0 0 0; 
     0 C 0; 
     0 0 0];

f_dae = @(t, Y) [
    Y(1) - Y(3)*R;         % 식 1: Vr = I * R
    Y(3);                  % 식 2: C*Vc' = I (M 행렬에 의해 C가 곱해짐)
    Vin - Y(1) - Y(2)      % 식 3: Vin - Vr - Vc = 0
];

opts = odeset('Mass', M, 'MassSingular', 'yes');
% 초기값: [Vr; Vc; I] = [5; 0; 5/100]
y0 = [Vin; 0; Vin/R];

[t_dae, y_dae] = ode15s(f_dae, tspan, y0, opts);

그림 3. RC회로에 대한 ODE 모델링과 DAE 모델링

진자 운동 (Pendulum)의 ODE 풀이와 DAE 풀이

위 RC 회로 모델링으로 ODE와 DAE의 모델링 철학을 이해했다면 이번에는 좀 더 어려운 예제를 풀어보고자 한다. 이 예제는 앞선 RC 회로 모델링에 비해서는 좀 더 어려운데, 굳이 가져온 이유는 DAE를 풀 때 필요한 index의 개념과 라그랑즈 승수($\lambda$)에 대해 설명하기 위함이다. 따라서, index나 제약 사항을 기술하는 라그랑주 승수대해 따라올 필요가 없는 경우 과감히 스킵해도 괜찮다.

DAE의 “Index (인덱스)” 개념에 대해서

아래의 ODE와 DAE의 예제를 풀 때 인덱스 (Index)의 개념이 필요해 개념 정리를 할 필요가 있다. Index는 DAE를 ODE로 바꾸기 위해 제약 조건을 몇 번 미분해야 하는지를 나타내는 척도이다. Index가 0인 방정식은 일반 상미분방정식을 나타내고, 대수 방정식을 한 번 미분하면 상미분방정식으로 표현되는 상태를 Index 1, 최소 두 번 미분해야 상미분 방정식으로 나타낼 수 있는 대수 방정식의 index는 2가 되는 것이다.

ODE 풀이

우선 각도 $\theta$만 이용해서 ODE 방식으로 푸는 경우이다. 상태 변수는 $[\theta, \dot{\theta}]$로 놓고 풀 수 있을 것이다.

% --- ODE Solver (ode45) ---
L = 1; g = 9.81;
f_ode = @(t, y) [y(2); -(g/L)*sin(y(1))]; % y(1)=theta, y(2)=dot_theta

[t_ode, y_ode] = ode45(f_ode, [0 5], [pi/4; 0]); % 45도에서 시작

DAE 풀이

이번엔 같은 문제를 DAE로 풀어보자. 제약 조건으로 줄 대수 방정식(algebraic equation)은

으로 줄 수 있다. 이 조건을 풀어서 설명해보자면 진자를 달고 있는 끈이 “강체”(rigid body)이고 그 길이는 항상 $L$이라는 제약조건인 것이다.

그런데, MATLAB에서 ode15s로 DAE를 풀기 위해서는 식(3)과 같은 형태로는 index 3의 수준에서는 바로 풀 수 없고, 제약 조건을 딱 한번만 미분하면 미분방정식 형태로 바꿀 수 있는 수준인 index 1까지 미분해서 제공해야한다. 즉, 위치에 관한 물리적 제약을 가속도 수준까지 풀어서 줘야한다. 그러기 위해선 이 제약 조건을 “장력”의 관점에서 다시 풀어서 생각해볼 필요가 있다. 식 (3)의 물리적 의미에 따라 무슨 일이 있어도 끈의 길이를 “L”로 맞추려면 pendulum이 $(x, y)$ 위치에 있을 때 $(x, y)$ 벡터의 방향에 $\lambda$만큼 비례하는 힘 (장력)을 준다는 것이다.

그림 4. 진자에 작용하는 장력(tension)과 중력(gravity)

따라서, 장력 벡터 $T$를 이렇게 정의할 수 있다.

여기서 $\lambda$는 라그랑주 승수 (링크의 챕터 2)이며 시간에 따라 변하는 변수이다. 우리의 솔버 ode15s는 이 $\lambda$ 변수를 포함해서 문제를 풀어주게 된다.

이제 식 (3)을 미분해가면서 최소 한번만 미분해주면 대수 방정식을 미분방정식으로 만들 수 있게 바꿔보자. (미분이 필요한 최소한의 값을 지표로 index 라고 하며, 이 경우 index를 3을 1로 만드는 것이다.)

우선 식 (3)의 위치 제약 식을 미분하면 아래와 같다.

여기서 $\dot{x}=u, \dot{y}=v$이다. 식 (7)을 한번 더 미분하면,

그림 4에서 볼 수 있는 장력과 중력의 힘들을 이용해 x 방향에 작용하는 힘과 y 방향으로 작용하는 힘들을 각각 계산해서 $F=ma$ 형태로 묶어보면,

식 (11), (12)를 식 (10)에 대입하면 아래와 같은 식을 얻을 수 있다.

상태변수는 $[x, y, u, v, \lambda]$ 총 5개이고 여기서 $\lambda$는 제약 조건을 유지하기 위한 힘 혹은 장력에 관한 라그랑주 승수이다.

% --- DAE Solver (ode15s) ---
m = 1; L = 1; g = 9.81;

% Mass Matrix
M = diag([1 1 1 1 0]); 

f_dae = @(t, Y) [
    Y(3);                               % x' = u
    Y(4);                               % y' = v
    Y(5)*Y(1)/m;                        % u' = (lambda * x)/m
    (Y(5)*Y(2) - m*g)/m;                % v' = (lambda * y - mg)/m
    % Index-1 제약 조건: 가속도 수준 (x*u' + y*v' + u^2 + v^2 = 0)
    % 여기에 u'과 v' 자리에 위의 식들을 대입하여 정리한 형태입니다.
    Y(3)^2 + (Y(1)*(Y(5)*Y(1)/m)+ Y(4)^2 + Y(2)*((Y(5)*Y(2) - m*g)/m) )
];

opts = odeset('Mass', M, 'MassSingular', 'yes');
y0 = [sin(pi/4); -cos(pi/4); 0; 0; 0]; 

[t_dae, y_dae] = ode15s(f_dae, [0 5], y0, opts);

그림 5. ODE와 DAE의 풀이는 거의 일치하는 것을 알 수 있다.

Numerical Drift

그런데, 그림 5의 진자 운동의 ODE와 DAE의 풀이를 time span을 늘려서 5초까지가 아니라 10초까지로 늘려서 시뮬레이션해보면 ODE와 DAE의 차이가 벌어지는 것을 알 수 있다.

그림 6. 대수방정식의 index를 낮춰 해를 구하면 numerical drift가 발생할 수 있다.

그림 6에서와 같이 5초 이후로는 계속해서 에러가 누적되는데, 이러한 에러가 발생하는 이유는 대수방정식 제약사항을 index 3에서 index 1로 낮춰 관리했기 때문이다. 즉, 우리가 궁극적으로 보고 싶은 것은 진자가 L이라는 길이의 실 위서 그 길이를 지키며 잘 있는지 궁금했지만 ode15s 솔버가 지키고자 하는 조건은 가속도가 특정 식을 만족하도록 $\lambda$에 적절한 값을 찾는 것이었기 때문이다.

그래서, Simscape나 PSPICE에서는?

지금까지 우리는 물리 현상을 모델링 및 시뮬레이션 하기 위해 미분방정식을 세우거나, 대수방정식이 포함된 연립 미분방정식을 풀면서 의미롤 고찰했다. 하지만 실제 산업 현장에서 수천 개의 부품이 얽힌 시스템을 매번 이렇게 수동으로 모델링하는 것은 불가능에 가깝다. 여기서 물리 모델링 소프트웨어의 진가가 드러난다.

1. 비인과적 모델링 (Acausal Modeling)의 실현

우리가 MATLAB 코드를 짤 때 f_dae에 식을 나열했던 것처럼, PSPICE나 Simscape는 사용자가 블록을 선으로 연결하는 행위 자체를 대수 방정식($\mathbf{0} = \mathbf{g}$)의 생성으로 간주한다.

• 저항 블록을 갖다 놓으면 $V=IR$이 추가된다.
• 두 블록을 선으로 이으면 그 지점의 전압이 같다는 제약($V_1 = V_2$)과 전류의 합이 0이라는 제약($\sum I = 0$)이 자동으로 수식 시스템에 편입된다. 이것이 바로 입력과 출력을 미리 정하지 않아도 되는 비인과적 모델링의 핵심이다.

비인과적 모델링이 중요한 것은 물리현상을 모델링할 때에는 특정 현상이 인과적인 경우가 드물기 때문이다. 예를 들어, 두 기어가 맞물려 돌아가는 것은 하나의 기어가 다른 기어를 움직이는 인과 관계를 가지기 때문인가? 아니면 그 반대로 작동하는 것인가? 기어는 서로 맞물려 돌아가기 때문에 비인과적 관계를 가진다.

2. 자동 Index Reduction

우리가 진자 문제에서 Index 3 식을 Index 1으로 바꾸기 위해 손으로 두 번 미분했던 과정을 기억하는가? Simscape 엔진은 (정확한 작동 원리는 모르지만, 유저 입장에서 봤을 때) 이 과정을 자동으로 수행해준다. 사용자가 어떤 복잡한 구속 조건을 걸더라도, 솔버가 풀 수 있는 최적의 형태로 수식을 알아서 변형(여기선 Reduction) 해준다.

3. 수치적 오차의 강건한 보정

우리가 직접 짠 DAE 코드에서 주기가 미세하게 어긋났던 이유는 가속도 제약만 지키다 보니 위치 제약이 조금씩 풀렸기 때문이다. 구체적으로 공개된 정보는 없지만 Simscape는 이와 같은 오차를 막기 위한 transient initialization이나 transient solving 알고리즘을 이용해 결과를 보정한다 [2]. 덕분에 장시간 시뮬레이션에서도 에너지가 보존되고 물리적 일관성이 유지된다.

좋아. 지금까지의 논의를 종합해서, Simscape 시뮬레이션의 안정성과 솔버의 철학을 다루는 핵심 챕터를 작성해 보았다. 이 내용은 단순한 사용법을 넘어 “왜 내 모델이 터지는가?”에 대한 본질적인 답을 제시해 줄 것이다.

시뮬레이션의 사투: Stiff 문제와 솔버의 선택

지금까지 DAE와 물리 시뮬레이션 소프트웨어에서의 사용에 대해 약간 알아보았다. 이와 관련하여 예를 들어 Simscape를 사용할 때 자주 발생하는 솔버가 “터져버리는” 문제에 대해서도 잠깐 생각해보자.

Simscape 모델을 돌리다 보면 “솔버가 스텝 사이즈를 줄이다가 멈췄습니다”라거나 “수렴에 실패했습니다”라는 에러를 자주 만난다. 이는 우리가 푸는 물리 모델이 수학적으로 Stiff하기 때문에 발생하는 현상이다. (Stiff를 “경직성”이라고도 번역하는 것 같은데 한글 용어가 익숙치 않아서 그냥 Stiff로 쓰겠다.)

Stiff 문제: “억울하게” 천천히 가야 하는 상태

Simscape 모델을 돌리다 보면 왜 자꾸 “스텝 사이즈를 줄이다가 멈췄다”는 에러가 날까? 보통 ‘Stiff(경직성) 문제’라고 하면 해가 급격히 변하는 구간을 떠올리기 쉽다. 그런데 이렇게 생각해보자. 해가 급격히 변할 때 스텝을 줄이는 건 너무나 당연한 일이다. 변화가 빠르니까 세밀하게 관찰해야 하기 때문이다. 그걸 ‘경직되었다’고 표현하는 건 좀 이상하지 않은가?

실제로 MathWorks의 ode15s 솔버의 문서를 뒤져보고 나서야 무릎을 탁 쳤다. 진짜 Stiff의 본질은 “해는 전반적으로 완만하게 변하고 있는데, 시스템 내부에 숨어 있는 민감한 성분 때문에 억지로 기어가야 하는 상태”를 말하는 것이었다. 다시 말해 Stiff하다는 건, 시스템의 ‘가장 빠른 성분’에 발목이 잡힌 상태다.

• 사고의 전환: 겉으로는 평온해 보이는 평지인데, 사실 발밑에 아주 예민한 지뢰(빠른 고유값)가 깔려 있는 상황이다. ode45 같은 일반적인 Explicit(명시적) 솔버는 이 지뢰를 밟지 않으려고 아주 조심조심(작은 스텝) 이동하다가 결국 지쳐서 포기해버린다. 아니면 잠깐 방심해서 스텝을 크게 딛는 순간, 지뢰가 터지면서 시스템의 에너지가 무한대로 발산하며 시뮬레이션이 ‘폭발’해버리기도 한다.
• Simscape는 왜 더 심할까? 우리가 앞에서 다룬 식 (2)의 대수 방정식($g=0$) 성분 때문이다. “그 즉시” 평형을 맞춰야 한다는 제약은 물리적으로 변화 속도가 ‘무한대’인 성분과 같다. 이 녀석들이 시스템 도처에 지뢰를 깔아놓으니, Simscape 모델은 태생적으로 Stiff할 수밖에 없다.

그래서 우리는 ode15s 같은 Implicit 솔버를 써야 한다. 이런 종류의 솔버는 단순히 현재 기울기로 앞을 예측하는 게 아니라, 미래의 지점이 평형을 이루는지 미리 방정식을 풀어서 확인하고 이동한다. 덕분에 발밑에 지뢰가 깔려 있어도 수치적으로 폭발하지 않고 안정적으로 큰 걸음을 내디딜 수 있는 것이다.

수렴 실패: “네 모델은 물리적으로 말이 안 돼!”

자, 이제 솔버까지 잘 골랐는데 이번엔 “수렴 실패(Convergence Error)”라는 메시지가 우리를 괴롭힌다. 이건 또 왜 그럴까? Implicit 솔버가 다음 지점을 찾는 방식은 일종의 ‘정답 맞히기 게임’이다. 이때 사용하는 도구가 바로 Newton-Raphson(NR) 반복법이다.

솔버는 스스로에게 이런 질문을 던진다. “이 지점에서 전류의 합이 0인가? 알짜힘은 0인가?” 이 질문에 대한 오차를 잔차(Residual)라고 부르는데, 이 잔차가 0이 될 때까지 계속해서 계산을 반복한다. 즉, 물리 법칙을 완벽히 만족하는 평형 상태를 찾아가는 과정이다.

그런데 여기서 시뮬레이션이 터지는 진짜 이유가 나온다.

• 물리적 불연속성과 에너지 점프: 만약 내가 모델을 설계할 때, 물체가 딱딱한 벽을 순식간에 뚫고 지나가야 하거나 에너지가 찰나의 순간에 점프해야 하도록 만들었다면? 혹은 인덕터 두 개를 직렬로 연결해서 전류가 부딪히게 만들었다면 어떻게 될까?
• 솔버의 절규: 솔버는 NR 반복을 수천 번 돌려도 잔차를 0으로 만드는 해를 찾을 수 없다. 물리적으로 그런 상태는 존재할 수 없기 때문이다! 잔차가 줄어들지 않고 요동치거나 엉뚱한 곳으로 튀어버린다.

결국 수렴 실패라는 건, 솔버가 나에게 보내는 마지막 경고이다. “아무리 계산해봐도 네가 만든 상황은 물리 법칙상 해가 존재할 수 없어!”라고 말이다. 결국 에러를 고치는 건 솔버 옵션을 만지는 게 아니라, 내 모델의 물리적 개연성에 문제가 있지는 않은지, 에너지 흐름이 너무 급격하게 변하는 구간은 없는지 다시 검토하는 것에서 시작해야 한다는 것을 깨달았다.

“DAE는 안 풀리는 문제가 있는데, ODE로는 풀릴까요?”

한 구독자분께서 아주 날카로운 질문을 던져주셨다.

“DAE 솔버가 물리 법칙을 만족하지 못해서 해를 못 찾는다면… 수식을 잘 정리해서 ODE로 만들면 돌아갈까요? DAE로는 안 풀리는 문제가 ODE로는 풀릴 수 있나요?”

이 질문을 듣고 나도 잠시 멈춰서 생각에 빠졌다. “돌아가는 것”과 “정답인 것”의 차이를 고민해 볼 수 있는 아주 본질적인 질문이기 때문이다. 꼭 이 답변이 정답은 아닐 수 있겠지만 내 사고의 과정을 따라 답을 적어본다.

그림 7. 철저한 원칙주의로 물리적 모순에 멈춰서는 DAE와, 수치적 오차를 감수하며 유연하게 해를 찾아가는 ODE의 대비 (Google 나노바나나 이용하여 그림 제작함.)

1) DAE 솔버의 ‘결벽증’

우선 Simscape 같은 DAE 솔버의 성격을 이해해야 한다. 앞에서 말했듯 얘네는 Newton-Raphson 반복을 통해 제약 조건을 엄격하게 체크한다.

이렇게 생각해보자. DAE 솔버는 “기찻길 위를 달리는 기차”다. 철로가 0.1mm라도 어긋나 있으면 기차는 탈선할 위험을 감수하지 않고 그 자리에 즉시 멈춰버린다. “전류의 합이 0이 아니잖아? 난 이대로는 단 한 발자국도 못 가!” 이런 결벽증 때문에 물리적으로 아주 미세하게 모순되는 지점(불연속적인 에너지 점프 등)만 만나도 솔버는 “수렴 실패”를 선언하고 파업에 들어간다.

2) ODE 솔버의 ‘무던함’

반면, 수식을 열심히 정리해서 ODE 형태로 만들면 상황이 묘하게 달라진다. ODE는 제약 조건을 이미 수식 안에 녹여버렸기 때문에, 솔버 입장에서는 실시간으로 감시해야 할 ‘대수 방정식’이라는 감옥이 사라진 셈이다.

ODE 솔버는 앞서 비유한 것 처럼 “지도를 들고 초원을 달리는 강아지”와 같다. 강아지는 지도가 좀 잘못 그려져 있어도, 혹은 앞에 작은 웅덩이가 있어도 일단 앞발을 내디딘다. 수치적 오차가 발생하더라도 미분 계수라는 ‘경향성’을 따라 적분하며 다음 지점으로 어떻게든 넘어간다.

3) 돌아는 가는데, 그게 ‘진짜’일까?

그래서 질문에 대한 답을 내리자면 이렇다. “ODE로 바꾸면 돌아갈 수는 있다. 하지만 그 결과가 물리적으로 맞다는 보장은 없다.”

DAE 솔버가 멈췄던 지점에서 ODE 솔버는 수치적 오차를 내며 억지로 뚫고 지나갈 것이다. 시뮬레이션은 끝까지 완료되겠지만, 데이터를 뜯어보면 에너지가 보존되지 않거나 물리 법칙을 무시하고 값이 튀어버린 ‘가짜 결과’일 가능성이 높다.

4) 결론: DAE의 멈춤은 ‘경고’다

결국 “DAE로는 안 풀리는데 ODE로는 풀리는 상황”은 솔버의 성능 차이가 아니라, 솔버가 물리 법칙을 대하는 태도의 차이에서 온다.

DAE 솔버가 까다롭게 구는 이유는 우리가 만든 가상 세계가 현실의 물리 법칙을 완벽히 대변하게 만들기 위해서다. 만약 DAE가 해를 찾지 못하고 멈췄다면, ODE로 우회할 수도 있겠지만 “내 모델의 물리적 구조 어디에 모순이 숨어 있는가”를 들여다볼 필요도 있다.

마치며

대수방정식에 대해서 오랫동안 궁금했는데, 이번 포스팅을 계기로 약간의 궁금증이 해소되었다. 물리 모델링을 수행할 때 블록을 불러들이는 것과 블록들을 선으로 연결하는 행위의 뒷면에는 이처럼 놀라운 계산들이 숨어있었다는게 알아보면서도 놀라울 따름이다. 물리 모델링의 백스테이지에 대해 궁금했던 나같은 사람들에게 도움이 되었으면 좋겠다.

참고 문헌

• [1] Solve Differential Algebraic Equations (DAEs)
• [2] How Simscape Simulation Works
• [3] ode15s - 경직성(stiff) ODE 풀기 예제