最適化変数を使用して ODE のパラメーターを当てはめる

この例では、最適化変数 (問題ベースのアプローチ) を使用して、常微分方程式 (ODE) を最小二乗法的に最適化するパラメーターを求める方法を示します。

問題

この問題は、複数ステップの反応モデルであり、複数の物質が関与し、そのうちいくつかは相互に反応して、別の物質を生成します。

この問題では、実際の反応速度が不明です。そのため、反応を観測し、速度を推定する必要があります。物質は一連の時間 $$t$$ の間、測定できると仮定します。これらの観測により、最適な一連の反応速度を測定結果に当てはめます。

モデル

モデルには、$$C_1$$ ～ $$C_6$$ の 6 つの物質があり、次のように反応します。

• $$C_1$$ と $$C_2$$ が 1 つずつ反応して、速度 $$r_1$$ で $$C_3$$ を 1 つ生成

• $$C_3$$ と $$C_4$$ が 1 つずつ反応して、速度 $$r_2$$ で $$C_5$$ を 1 つ生成

• $$C_3$$ と $$C_4$$ が 1 つずつ反応して、速度 $$r_3$$ で $$C_6$$ を 1 つ生成

反応速度は、必要な物質の量の積に比例します。したがって、$$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{\mathit{dy}}{\mathit{dt}}=\left[\begin{array}{c}-{\mathit{r}}_1{\mathit{y}}_1{\mathit{y}}_2\\ -{\mathit{r}}_1{\mathit{y}}_1{\mathit{y}}_2\\ -{\mathit{r}}_2{\mathit{y}}_3{\mathit{y}}_4+{\mathit{r}}_1{\mathit{y}}_1{\mathit{y}}_2-{\mathit{r}}_3{\mathit{y}}_3{\mathit{y}}_4\\ -{\mathit{r}}_2{\mathit{y}}_3{\mathit{y}}_4-{\mathit{r}}_3{\mathit{y}}_3{\mathit{y}}_4\\ {\mathit{r}}_2{\mathit{y}}_3{\mathit{y}}_4\\ {\mathit{r}}_3{\mathit{y}}_3{\mathit{y}}_4\end{array}\right].$$

時間 0、点 $$y(0)=[1,1,0,1,0,0]$$ から、微分方程式を開始します。これらの初期値では、すべての物質が完全に反応し、時間の増加とともに $$C_1$$ ～ $$C_4$$ がゼロに近づきます。

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 を使用して、関数を最適化式に変換します。詳細については、非線形関数から最適化式への変換を参照してください。

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

差の二乗和はほぼゼロであり、ODE の解と実際のパラメーターによる解を一致させるパラメーターをソルバーが見つけたことを意味しています。そのため、期待どおり、解には実際のパラメーターが含まれています。

disp(rsol.r)

    2.5000
    1.2000
    0.4500

disp(rtrue)

    2.5000    1.2000    0.4500

限られた観測値

y の成分すべてではなく、y(5) と y(6) の最終出力のみを観測できるとします。この限られた情報に基づいて、すべての反応速度の値を求められるでしょうか。

確認するために、関数 RtoODE を変更して、5 番目と 6 番目の ODE 出力のみを返すようにします。変更した 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 を呼び出し、出力の最後の 2 行を取り出すだけです。

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

ここでも、返される二乗和はほぼゼロです。これは、正しい反応速度をソルバーが見つけたということでしょうか。

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

この問題の正しい反応速度を特定するには、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

関連するトピック

• 問題ベースの最適化ワークフロー
• 離散化された最適軌跡 (問題ベース)