Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

ODE の適合、問題ベース

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

問題

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

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

モデル

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

  • C1C2 が 1 つずつ反応して、速度 r1C3 を 1 つ生成

  • C3C4 が 1 つずつ反応して、速度 r2C5 を 1 つ生成

  • C3C4 が 1 つずつ反応して、速度 r3C6 を 1 つ生成

反応速度は、必要な物質の量の積に比例します。したがって、yi が物質 Ci の量を表す場合、C3 を生成する反応速度は r1y1y2 となります。同様に、C5 を生成する反応速度は r2y3y4 となり、C6 を生成する反応速度は r3y3y4 となります。

つまり、この系の変化を制御する微分方程式は、次のとおりです。

dydt=[-r1y1y2-r1y1y2-r2y3y4+r1y1y2-r3y3y4-r2y3y4-r3y3y4r2y3y4r3y3y4].

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

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

実際の反応速度は、r1=2.5r2=1.2r3=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);

show を呼び出して問題を表示します。

show(prob)
  OptimizationProblem : 

	Solve for:
       r

	minimize :
       sum(sum((RtoODE(r, extraParams{1}, extraParams{2})
     - extraParams{3}).^2, 1))

         extraParams{1}:

       Columns 1 through 7

              0    0.0505    0.1010    0.1515    0.2020    0.2525    0.3030

       Columns 8 through 14

         0.3535    0.4040    0.4545    0.5051    0.5556    0.6061    0.6566

       Columns 15 through 21

         0.7071    0.7576    0.8081    0.8586    0.9091    0.9596    1.0101

       Columns 22 through 28

         1.0606    1.1111    1.1616    1.2121    1.2626    1.3131    1.3636

       Columns 29 through 35

         1.4141    1.4646    1.5152    1.5657    1.6162    1.6667    1.7172

       Columns 36 through 42

         1.7677    1.8182    1.8687    1.9192    1.9697    2.0202    2.0707

       Columns 43 through 49

         2.1212    2.1717    2.2222    2.2727    2.3232    2.3737    2.4242

       Columns 50 through 56

         2.4747    2.5253    2.5758    2.6263    2.6768    2.7273    2.7778

       Columns 57 through 63

         2.8283    2.8788    2.9293    2.9798    3.0303    3.0808    3.1313

       Columns 64 through 70

         3.1818    3.2323    3.2828    3.3333    3.3838    3.4343    3.4848

       Columns 71 through 77

         3.5354    3.5859    3.6364    3.6869    3.7374    3.7879    3.8384

       Columns 78 through 84

         3.8889    3.9394    3.9899    4.0404    4.0909    4.1414    4.1919

       Columns 85 through 91

         4.2424    4.2929    4.3434    4.3939    4.4444    4.4949    4.5455

       Columns 92 through 98

         4.5960    4.6465    4.6970    4.7475    4.7980    4.8485    4.8990

       Columns 99 through 100

         4.9495    5.0000

       extraParams{2}:

          1     1     0     1     0     0

       extraParams{3}:

       Columns 1 through 7

         1.0000    0.8879    0.7984    0.7253    0.6644    0.6130    0.5690
         1.0000    0.8879    0.7984    0.7253    0.6644    0.6130    0.5690
              0    0.1074    0.1847    0.2404    0.2805    0.3089    0.3287
         1.0000    0.9953    0.9831    0.9657    0.9449    0.9219    0.8977
              0    0.0034    0.0123    0.0249    0.0401    0.0568    0.0744
              0    0.0013    0.0046    0.0094    0.0150    0.0213    0.0279

       Columns 8 through 14

         0.5308    0.4975    0.4681    0.4420    0.4186    0.3976    0.3786
         0.5308    0.4975    0.4681    0.4420    0.4186    0.3976    0.3786
         0.3421    0.3506    0.3554    0.3574    0.3573    0.3556    0.3527
         0.8729    0.8481    0.8235    0.7994    0.7759    0.7532    0.7313
         0.0924    0.1105    0.1284    0.1459    0.1630    0.1795    0.1954
         0.0347    0.0414    0.0481    0.0547    0.0611    0.0673    0.0733

       Columns 15 through 21

         0.3613    0.3456    0.3311    0.3178    0.3056    0.2942    0.2837
         0.3613    0.3456    0.3311    0.3178    0.3056    0.2942    0.2837
         0.3489    0.3444    0.3395    0.3342    0.3287    0.3230    0.3173
         0.7102    0.6900    0.6706    0.6520    0.6343    0.6173    0.6010
         0.2108    0.2255    0.2396    0.2531    0.2660    0.2783    0.2902
         0.0790    0.0846    0.0898    0.0949    0.0997    0.1044    0.1088

       Columns 22 through 28

         0.2739    0.2647    0.2562    0.2481    0.2406    0.2335    0.2268
         0.2739    0.2647    0.2562    0.2481    0.2406    0.2335    0.2268
         0.3116    0.3059    0.3002    0.2946    0.2891    0.2837    0.2784
         0.5855    0.5706    0.5564    0.5428    0.5297    0.5172    0.5052
         0.3015    0.3123    0.3226    0.3325    0.3420    0.3511    0.3598
         0.1131    0.1171    0.1210    0.1247    0.1283    0.1317    0.1349

       Columns 29 through 35

         0.2205    0.2146    0.2089    0.2035    0.1984    0.1936    0.1890
         0.2205    0.2146    0.2089    0.2035    0.1984    0.1936    0.1890
         0.2732    0.2682    0.2633    0.2585    0.2538    0.2493    0.2449
         0.4938    0.4827    0.4722    0.4620    0.4523    0.4429    0.4339
         0.3682    0.3762    0.3839    0.3913    0.3984    0.4052    0.4117
         0.1381    0.1411    0.1440    0.1467    0.1494    0.1519    0.1544

       Columns 36 through 42

         0.1846    0.1804    0.1763    0.1725    0.1688    0.1653    0.1619
         0.1846    0.1804    0.1763    0.1725    0.1688    0.1653    0.1619
         0.2406    0.2364    0.2324    0.2285    0.2246    0.2209    0.2173
         0.4252    0.4168    0.4087    0.4010    0.3935    0.3862    0.3792
         0.4181    0.4241    0.4300    0.4357    0.4411    0.4464    0.4515
         0.1568    0.1591    0.1613    0.1634    0.1654    0.1674    0.1693

       Columns 43 through 49

         0.1587    0.1556    0.1526    0.1497    0.1469    0.1442    0.1416
         0.1587    0.1556    0.1526    0.1497    0.1469    0.1442    0.1416
         0.2138    0.2104    0.2071    0.2039    0.2007    0.1977    0.1947
         0.3725    0.3660    0.3596    0.3535    0.3476    0.3419    0.3364
         0.4564    0.4611    0.4657    0.4702    0.4744    0.4786    0.4826
         0.1711    0.1729    0.1746    0.1763    0.1779    0.1795    0.1810

       Columns 50 through 56

         0.1392    0.1368    0.1344    0.1322    0.1300    0.1279    0.1259
         0.1392    0.1368    0.1344    0.1322    0.1300    0.1279    0.1259
         0.1918    0.1890    0.1863    0.1836    0.1810    0.1785    0.1761
         0.3310    0.3258    0.3207    0.3158    0.3111    0.3064    0.3019
         0.4866    0.4903    0.4940    0.4976    0.5010    0.5044    0.5077
         0.1825    0.1839    0.1853    0.1866    0.1879    0.1892    0.1904

       Columns 57 through 63

         0.1239    0.1220    0.1202    0.1184    0.1166    0.1149    0.1133
         0.1239    0.1220    0.1202    0.1184    0.1166    0.1149    0.1133
         0.1737    0.1713    0.1690    0.1668    0.1646    0.1625    0.1605
         0.2976    0.2933    0.2892    0.2852    0.2813    0.2775    0.2737
         0.5109    0.5139    0.5169    0.5199    0.5227    0.5255    0.5282
         0.1916    0.1927    0.1939    0.1950    0.1960    0.1971    0.1981

       Columns 64 through 70

         0.1117    0.1101    0.1086    0.1072    0.1057    0.1043    0.1030
         0.1117    0.1101    0.1086    0.1072    0.1057    0.1043    0.1030
         0.1584    0.1565    0.1546    0.1527    0.1508    0.1491    0.1473
         0.2701    0.2666    0.2632    0.2598    0.2566    0.2534    0.2503
         0.5308    0.5334    0.5359    0.5383    0.5407    0.5430    0.5453
         0.1991    0.2000    0.2010    0.2019    0.2028    0.2036    0.2045

       Columns 71 through 77

         0.1017    0.1004    0.0991    0.0979    0.0967    0.0955    0.0944
         0.1017    0.1004    0.0991    0.0979    0.0967    0.0955    0.0944
         0.1456    0.1439    0.1423    0.1407    0.1391    0.1376    0.1361
         0.2472    0.2443    0.2414    0.2385    0.2358    0.2331    0.2304
         0.5475    0.5496    0.5517    0.5538    0.5558    0.5578    0.5597
         0.2053    0.2061    0.2069    0.2077    0.2084    0.2092    0.2099

       Columns 78 through 84

         0.0933    0.0922    0.0911    0.0901    0.0891    0.0881    0.0871
         0.0933    0.0922    0.0911    0.0901    0.0891    0.0881    0.0871
         0.1346    0.1331    0.1317    0.1303    0.1290    0.1277    0.1264
         0.2279    0.2253    0.2229    0.2204    0.2181    0.2157    0.2135
         0.5616    0.5634    0.5652    0.5670    0.5687    0.5704    0.5720
         0.2106    0.2113    0.2119    0.2126    0.2133    0.2139    0.2145

       Columns 85 through 91

         0.0862    0.0852    0.0843    0.0834    0.0826    0.0817    0.0809
         0.0862    0.0852    0.0843    0.0834    0.0826    0.0817    0.0809
         0.1251    0.1238    0.1226    0.1214    0.1202    0.1191    0.1179
         0.2112    0.2091    0.2069    0.2048    0.2028    0.2008    0.1988
         0.5736    0.5752    0.5768    0.5783    0.5798    0.5813    0.5827
         0.2151    0.2157    0.2163    0.2169    0.2174    0.2180    0.2185

       Columns 92 through 98

         0.0801    0.0793    0.0785    0.0777    0.0770    0.0762    0.0755
         0.0801    0.0793    0.0785    0.0777    0.0770    0.0762    0.0755
         0.1168    0.1157    0.1146    0.1136    0.1125    0.1115    0.1105
         0.1969    0.1950    0.1931    0.1913    0.1895    0.1877    0.1860
         0.5841    0.5855    0.5868    0.5882    0.5895    0.5907    0.5920
         0.2190    0.2196    0.2201    0.2206    0.2210    0.2215    0.2220

       Columns 99 through 100

         0.0748    0.0741
         0.0748    0.0741
         0.1095    0.1086
         0.1843    0.1826
         0.5932    0.5944
         0.2225    0.2229




	variable bounds:
       0.1 <= r(1) <= 10
       0.1 <= r(2) <= 10
       0.1 <= r(3) <= 10

問題を解く

最適適合パラメーター 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.8660e-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

新しいパラメーターでは、物質 C1C2 は減りが遅く、物質 C3 はそれほど増えていません。しかし、物質 C4C5C6 は、新しいパラメーターと実際のパラメーターのどちらの場合も、まったく同じ変化をしています。

参考

| |

関連するトピック