on
모델링 - Newton-Raphson method로 모수 최적화
모델링 - 지난 글
목차
모수 추정이란
모델링을 하는 과정에 있어서 모수 추정은 매우 중요한 과정 중 하나이다. 통계학 과점에서 대상의 전체인 모집단을 분석하는 데에 너무나 많은 시간과 비용이 필요하기 때문에 모집단의 일부인 표본을 통해 모집단의 특성을 파악한다고 정의한다. 하지만 표본을 통해 추정하는 것은 100% 일치한다고 보장할 수 없다. 일부에 대한 특성만을 이용하여 전체와 얼마나 일치하는 지를 확인하는 것이기 때문에 추정한 값과 실제 값 사이에는 오차가 존재할 수 밖에 없다.
모델링도 마찬가지이다. 일부 실험을 통해 얻은 자료로 식물의 기작을 설명하고자 하는 것이기 때문에 오차는 필연적으로 발생할 것이다. 하지만 모델러인 우리는 이 오차를 가장 효율적이고 정확하게 줄이는 방법을 알아야한다. 그리고 이 과정은 앞에서 언급했듯 가장 중요한 과정중 하나로 작물 모델링의 교과서인 작물생육 모델링의 이론과 실제에서도 모수 추정의 방법을 가장 먼저 정의하고 설명한다.
이번 글에서는 널리 알려진 모수 추정 방법 중 하나인 Newton 법에 대하여 알아보도록 하겠다.
테일러 급수
어떤 특정한 함수의 성질을 이해할 때 그 함수를 직접적으로 관찰하기보다는 그 함수와 비슷한 함수를 찾아냄으로써 우리가 관심있는 함수의 성징을 이해할 수 있는데 이 때 사용되는 방법이 테일러 전개(Tayler’s Expension)이라 한다.
쉽게 말해 테일러 전개는 우리가 파악하기 어려운 함수들을 우리가 그동안 쉽게 알고 있었던 다양한 함수로 ‘근사’ 시켜서 원래 함수의 성질들을 파악하는 데 사용하는 것이다. 테일러 전개를 위해서는 우선 함수들을 거듭제곱으로 표현할 수 있어야 하는데, 이 때 표현되는 급수를 테일러 급수라 한다. 어떤 함수 f(x)가 무한번 미분가능한 함수라 할 때
위와 같은 다항함수로 표현할 수 있고, 이를 a에서 f(x)의 테일러 급수라 한다.
테일러 급수는 원래의 함수를 ‘근사’ 시킨 것 이므로 오차가 존재할 것이다. 위에서 언급했 듯 테일러 급수는 무한번 미분 가능한 함수이다. 따라서 다음과 같이 표현할 수 있다.
우변의 급수는 무한급수이므로 우리는 부분합의 극한을 생각할 수 있다. 즉, f(x)는 부분합으로 이루어진 수열의 극한을 의미한다. 무한급수
에서는 부분합을 이용하 듯이 테일러 급수에서도 부분합을 이용하는데 이를 n차 테일러 다항식(nth Taylor Polynomial)이라 부르고 함수 f(x)에
대하여 다음과 같이 정의한다.
a에서 f(x)의 n차 테일러 다항식
함수 f(x) = e^x에 대하여 0에서 T1(x), T2(x), T3(x)를 구하면
이를 통해 각각의 다항식을 그래프로 나타내면 n이 증가할 수록 Tn(x)의 그래프는 e^x의 그래프에 근접함을 예측할 수 있는 것이다.
정리하면, 테일러 급수가 필요한 이유는 우리가 잘 모르거나 복잡한 함수를 다루기 쉽고 이해하기 쉬운 다항함수 (ex, 1~3차 다항함수)로 대체시키기 위함이다. 또한 어떤 함수를 테일러 급수로 표현하면 그 함수의 특성을 분석하기가 조금 더 용이해진다. 테일러 급수는 또한 복잡한 함수를 1~3차 다항함수로 근사하여 사용하기 때문에 단순화시킨다고 표현할 수도 있다.
우리는 이러한 테일러 급수의 원리를 적용한 예시인 뉴턴법 (Newton-Raphson method)에 대하여 알아보겠다.
Newton-Raphson method
뉴턴법 (Newton-Raphson method)는 함수 f(x)에 대해 테일러 급수를 통해 근사 함수를 구하고 최소점에서 그 도함수가 0인 사실을 이용하는 것이다. 쉽게 말해 f(x) = 0의 해를 근사적으로 찾을 때 유용하게 사용된다.
뉴턴법은 기본적으로 f’(a)가 x = a에서의 접선의 기울기라는 미분의 기하학적 해석을 이용하며, 현재 x값에서 접선을 그리고 접선이 x축과 만나는 지점으로 x를 이동시켜 가면서 점진적으로 해를 찾는 방법이다.
다음 그림을 보면 이해하기 쉽다. 처음에 x = x1에서 시작했다면, 그 다음 x값은 x2가 될 것이고, x2에서 다시 접선을 그려보면 점차 실제 해에 가까워진다는 것이 원리이다. 그리고 이를 반복하다보면 해를 찾을 수 있는 것이 뉴턴법인 것이다.
이를 수식화하면 다음과 같다. 초기값에서 수렴할 때까지 계속 x를 이동시키며, x 값이 더이상 변화가 없을 경우 종료한다. |xt+1 - xt|이 매우 작은 값이면 뉴턴법을 종료하고 x = xt+1이 해, 즉 f(xt+1) = 0 라고 생각하는 것이다.
Newton-Raphson method의 특징
- f(x)가 연속이고 미분가능해야 한다.
- 해가 여러 개인 경우에는 초기값을 어떻게 주느냐에 따라서 뉴턴법으로 찾아지는 해가 달라질 수 있다. 초기값을 잘 주면 금방 해를 찾을 수 있지만 잘못 주면 시간이 오래 걸리거나 아에 해를 찾지 못할 수도 있다.
광합성 모델에서의 뉴턴법
다음은 광합성 모델 (gasexchange)에서 모수 추정을 위하여 뉴턴법을 사용한 예시이다. 최적화 모델로 배추모델에서 Energy Balance 모델로 부터 예측한 엽 내부 CO2 농도 및 온도 값을 일정 조건에서의 오차를 줄이기 위해 사용되었다.
식물의 잎은 공기, 기공, 세포 내부로 크게 3개의 레이어로 나눌 수 있다. 그 중 우리가 센서로 계측할 수 있는 환경조건(온도, CO2 농도 등)은 공기중 환경에 해당한다. 하지만 식물의 잎에는 기공이 존재하며, 공기 중의 환경에 따라 기공의 열고 닫힘이 조절되며 이에 따라 세포 내부의 온도 및 CO2 농도 또한 바뀌게 된다.
모델링을 하는 과정 에서 잎의 세포 내 온도 및 CO2 농도를 측정하기는 쉽지 않다. 따라서 잎의 기공의 열고 닫음을 컨트롤 하는 기공 전도도와 이에 따른 Energy Balance를 설명하는 수식을 적용하여 입의 세포 내 온도 및 CO2 농도 예측을 하고자 하였고, 이에 대한 정확도를 향상시키기 위해 뉴턴법을 적용하여 오차를 줄이고자 한 것이다.
잎의 레이어 (공기, 기공, 세포 내부)
Energy Balance 모델 수식
파이썬을 통해 구현한 광합성 최적화 모델
def routine(self):
Ca = self.Ca
Tleaf = self.Ta
Ci = self.Ca * 0.9
def func1(Ci, Tleaf):
An = self.leafAssim(Tleaf, Ci)[0]
gbc = self.gbc()
gsco = self.gsc(An, Tleaf)[0]
newCi = self.interCi(An, Tleaf)
return (newCi - Ci)
def func2(Ci, Tleaf):
An = self.leafAssim(Tleaf, Ci)[0]
gbc = self.gbc()
gsco = self.gsc(An, Tleaf)[0]
newTleaf = self.EnergyBal(An, Tleaf)[0]
return (newTleaf - Tleaf)
Ci, Tleaf = self.Newton_2Var(func1, func2, Ci, Tleaf)
func1 함수: 잎의 세포 내부 CO2 예측, func2: 잎의 세포 내부 온도 예측
def Newton_2Var(self, func1, func2, x0, y0, h=0.1, tolerance=1e-8, steps=1000):
for i in range(1, steps):
x0 = 0.01 if x0 == 0 else x0
y0 = 0.01 if y0 == 0 else y0
A = func1(x0, y0)
B = (func1(x0 + h, y0) - func1(x0 - h, y0)) / (2 * h)
C = (func1(x0, y0 + h) - func1(x0, y0 - h)) / (2 * h)
O = func2(x0, y0)
P = (func2(x0 + h, y0) - func2(x0 - h, y0)) / (2 * h)
Q = (func2(x0, y0 + h) - func2(x0, y0 - h)) / (2 * h)
aaa = C * P - B * Q
if aaa == 0:
print("Error: divide by 0")
aaa = 0.001
x = x0 + (C * O - A * Q) / (B * Q - C * P)
y = y0 + (B * O - A * P) / (C * P - B * Q)
error1 = abs((x - x0) / x)
error2 = abs((y - y0) / y)
if (error1 < tolerance) and (error2 < tolerance): break
x0, y0 = x, y
self.newci = x
self.newtleaf = y
return x, y
예측한 CO2, 온도를 각각 x0, y0에 넣어 뉴턴법 적용
작동 원리는 다음과 같다.
- Energy Balance 수식을 통해 예측한 CO2, 온도를 각각 x0, y0 변수에 적용
- 만약 x0 및 y0가 0일 경우 0.01 적용
- x0 및 y0을 x축으로 간주하며 반복할 때 마다 0.1씩 이동하고, x0, y0에서의 미분값을 구하기 위하여 중앙차분법 적용
- x0 와 y0로 2차원이기 때문에 연립방정식 적용을 위하여 자코비안 행렬적용하여 계산
- x-x0 및 y-y0의 오차가 1e-8보다 작으면 종료
자코비안 행렬 적용 원리