최적화 변수를 사용하여 ODE 파라미터 피팅하기

이 예제에서는 최적화 변수(문제 기반 접근법)를 사용하여, 최소제곱 측면에서 상미분 방정식(ODE)을 최적화하는 파라미터를 찾는 방법을 보여줍니다.

문제

이 문제는 여러 물질이 관련되며 그중 일부가 서로 반응하여 다른 물질을 생성하는 다중 스텝 반응 모델입니다.

이 문제에서 실제 반응 속도는 알려져 있지 않습니다. 따라서 직접 반응을 관측하고 속도를 유추해야 합니다. 시간 집합 $t$ 에 대해 물질을 측정할 수 있다고 가정합니다. 그러한 관측값에서 최적의 반응 속도 집합을 측정값에 피팅하십시오.

모델

모델에는 $C_{1}$ ~ $C_{6}$ 의 6가지 물질이 있고, 이 물질들은 다음과 같이 반응합니다.

한 개의 $C_{1}$ 과 한 개의 $C_{2}$ 가 속도 $r_{1}$ 에서 반응하여 한 개의 $C_{3}$ 을 생성함
한 개의 $C_{3}$ 과 한 개의 $C_{4}$ 가 속도 $r_{2}$ 에서 반응하여 한 개의 $C_{5}$ 를 생성함
한 개의 $C_{3}$ 과 한 개의 $C_{4}$ 가 속도 $r_{3}$ 에서 반응하여 한 개의 $C_{6}$ 을 생성함

반응 속도는 필요한 물질 수량의 곱에 비례합니다. 따라서 $y_{i}$ 가 물질 $C_{i}$ 의 수량을 나타내는 경우 $C_{3}$ 생성을 위한 반응 속도는 $r_{1} y_{1} y_{2}$ 입니다. 이와 유사하게 $C_{5}$ 생성을 위한 반응 속도는 $r_{2} y_{3} y_{4}$ 이고, $C_{6}$ 생성을 위한 반응 속도는 $r_{3} y_{3} y_{4}$ 입니다.

즉, 시스템의 변화를 제어하는 미분 방정식은 다음과 같습니다.

$\frac{dy}{dt} = [\begin{matrix} - r_{1} y_{1} y_{2} \\ - r_{1} y_{1} y_{2} \\ - r_{2} y_{3} y_{4} + r_{1} y_{1} y_{2} - r_{3} y_{3} y_{4} \\ - r_{2} y_{3} y_{4} - r_{3} y_{3} y_{4} \\ r_{2} y_{3} y_{4} \\ r_{3} y_{3} y_{4} \end{matrix}] .$

시간 0에 점 $y (0) = [1, 1, 0, 1, 0, 0]$ 에서 미분 방정식을 시작합니다. 이러한 초기값을 사용하면 모든 물질이 완전히 반응하여 시간이 늘어남에 따라 $C_{1}$ 부터 $C_{4}$ 가 0에 도달하게 됩니다.

MATLAB에서 모델 표현하기

diffun 함수는 ode45로 문제를 풀 수 있는 형식으로 미분 방정식을 구현합니다.

type diffun

function dydt = diffun(~,y,r)
dydt = zeros(6,1);
s12 = y(1)*y(2);
s34 = y(3)*y(4);

dydt(1) = -r(1)*s12;
dydt(2) = -r(1)*s12;
dydt(3) = -r(2)*s34 + r(1)*s12 - r(3)*s34;
dydt(4) = -r(2)*s34 - r(3)*s34;
dydt(5) = r(2)*s34;
dydt(6) = r(3)*s34;
end

실제 반응 속도는 $r_{1} = 2.5$ , $r_{2} = 1.2$ , $r_{3} = 0.45$ 입니다. ode45를 호출하여 시간 0부터 5까지의 시스템 변화를 계산합니다.

rtrue = [2.5 1.2 0.45];
y0 = [1 1 0 1 0 0];
tspan = linspace(0,5);
soltrue = ode45(@(t,y)diffun(t,y,rtrue),tspan,y0);
yvalstrue = deval(soltrue,tspan);
for i = 1:6
    subplot(3,2,i)
    plot(tspan,yvalstrue(i,:))
    title(['y(',num2str(i),')'])
end

최적화 문제

문제 기반 접근법으로 문제를 풀 수 있도록 준비하기 위해 하한이 0.1이고 상한이 10인, 요소를 3개 가진 최적화 변수 r을 만듭니다.

r = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);

이 문제의 목적 함수는 파라미터 r을 사용한 ODE 해와 실제 파라미터 yvals를 사용한 해 사이 차의 제곱합입니다. 이 목적 함수를 표현하기 위해 파라미터 r을 사용하여 ODE 해를 계산하는 MATLAB 함수를 먼저 작성합니다. 이 함수는 RtoODE 함수입니다.

type RtoODE

function solpts = RtoODE(r,tspan,y0)
sol = ode45(@(t,y)diffun(t,y,r),tspan,y0);
solpts = deval(sol,tspan);
end

목적 함수에서 RtoODE를 사용하기 위해 fcn2optimexpr을 사용하여 함수를 최적화 표현식으로 변환합니다. Convert Nonlinear Function to Optimization Expression 항목을 참조하십시오.

myfcn = fcn2optimexpr(@RtoODE,r,tspan,y0);

목적 함수를 ODE 해와 실제 파라미터를 사용한 해 사이 차의 제곱합으로 표현합니다.

obj = sum(sum((myfcn - yvalstrue).^2));

목적 함수 obj를 사용하는 최적화 문제를 만듭니다.

prob = optimproblem("Objective",obj);

문제 풀기

최적의 피팅 파라미터 r을 구하기 위해 솔버의 초기 추측값 r0을 제공하고 solve를 호출합니다.

r0.r = [1 1 1];
[rsol,sumsq] = solve(prob,r0)

Solving problem using lsqnonlin.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

rsol = struct with fields:
    r: [3x1 double]

sumsq = 3.8661e-15

차의 제곱합은 사실상 0입니다. 이는 ODE 해를 실제 파라미터를 사용한 해와 일치시키는 파라미터를 솔버에서 찾았음을 의미합니다. 따라서 예상대로 해에는 실제 파라미터가 포함됩니다.

disp(rsol.r)

    2.5000
    1.2000
    0.4500

disp(rtrue)

    2.5000    1.2000    0.4500

제한된 관측값

y의 모든 성분을 관측할 수는 없고 최종 출력값 y(5)와 y(6)만 관측할 수 있다고 가정해 보겠습니다. 이 제한된 정보에 기반하여 모든 반응 속도 값을 구할 수 있을까요?

이를 확인하기 위해 5번째와 6번째 ODE 출력값만 반환하도록 함수 RtoODE를 수정합니다. 수정된 ODE 솔버는 RtoODE2에 있습니다.

type RtoODE2

function solpts = RtoODE2(r,tspan,y0)
solpts = RtoODE(r,tspan,y0);
solpts = solpts([5,6],:); % Just y(5) and y(6)
end

RtoODE2 함수는 RtoODE를 호출한 다음, 출력값의 마지막 두 개 행만 취합니다.

RtoODE2와 최적화 변수 r, 시간 길이 데이터 tspan, 초기점 y0에서 새 최적화 표현식을 만듭니다.

myfcn2 = fcn2optimexpr(@RtoODE2,r,tspan,y0);

출력값 5와 6만 포함하도록 비교 데이터를 수정합니다.

yvals2 = yvalstrue([5,6],:);

최적화 표현식 myfcn2와 비교 데이터 yvals2에서 새 목적 함수 문제와 새 최적화 문제를 만듭니다.

obj2 = sum(sum((myfcn2 - yvals2).^2));
prob2 = optimproblem("Objective",obj2);

이 제한된 관측값 집합을 기반으로 문제를 풉니다.

[rsol2,sumsq2] = solve(prob2,r0)

Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

rsol2 = struct with fields:
    r: [3x1 double]

sumsq2 = 2.1616e-05

이번에도 반환된 제곱합은 사실상 0입니다. 이는 솔버가 정확한 반응 속도를 찾았다는 의미일까요?

disp(rsol2.r)

    1.7811
    1.5730
    0.5899

disp(rtrue)

    2.5000    1.2000    0.4500

아닙니다. 이 경우 새로 구한 속도는 실제 속도와 상당히 다릅니다. 그러나 실제 값과 비교한 새 ODE 해의 플롯은 y(5)와 y(6)이 실제 값과 일치함을 보여줍니다.

figure
plot(tspan,yvals2(1,:),'b-')
hold on
ss2 = RtoODE2(rsol2.r,tspan,y0);
plot(tspan,ss2(1,:),'r--')
plot(tspan,yvals2(2,:),'c-')
plot(tspan,ss2(2,:),'m--')
legend('True y(5)','New y(5)','True y(6)','New y(6)','Location','northwest')
hold off

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent True y(5), New y(5), True y(6), New y(6).

이 문제의 정확한 반응 속도를 파악하려면 y(5)와 y(6) 외에 더 많은 관측값 데이터가 필요합니다.

새 파라미터를 사용하여 해의 모든 성분을 플로팅하고, 실제 파라미터를 사용하여 해를 플로팅합니다.

figure
yvals2 = RtoODE(rsol2.r,tspan,y0);
for i = 1:6
    subplot(3,2,i)
    plot(tspan,yvalstrue(i,:),'b-',tspan,yvals2(i,:),'r--')
    legend('True','New','Location','best')
    title(['y(',num2str(i),')'])
end

새 파라미터를 사용할 경우 물질 $C_{1}$ 과 $C_{2}$ 는 상대적으로 천천히 고갈되고 물질 $C_{3}$ 은 그만큼 누적되지 않습니다. 하지만 물질 $C_{4}$ , $C_{5}$ , $C_{6}$ 은 새 파라미터와 실제 파라미터를 사용할 때 모두 정확하게 동일한 변화를 보입니다.

참고 항목

solve | ode45 | fcn2optimexpr