Подгонка ОДУ, основанная на проблеме

Открыть Live Script

Этот пример показывает, как найти параметры, которые оптимизируют обыкновенное дифференциальное уравнение (ОДУ) в смысле наименьших квадратов, используя основанный на проблеме подход.

Проблема

Проблема заключается в многоступенчатой модели реакции, включающей несколько веществ, некоторые из которых реагируют друг с другом с образованием различных веществ.

Для этой проблемы истинные скорости реакции неизвестны. Итак, нужно наблюдать реакции и выводить скорости. Предположим, что вы можете измерить вещества в течение ряда раз $t$ . Из этих наблюдений, соответствует лучшему набору скоростей реакции измерениям.

Модель

Модель имеет шесть веществ, $C_{1}$ через $C_{6}$ , которые реагируют следующим образом:

Один $C_{1}$ и один $C_{2}$ реагировать, чтобы сформировать один $C_{3}$ по курсу $r_{1}$
Один $C_{3}$ и один $C_{4}$ реагировать, чтобы сформировать один $C_{5}$ по курсу $r_{2}$
Один $C_{3}$ и один $C_{4}$ реагировать, чтобы сформировать один $C_{6}$ по курсу $r_{3}$

Скорость реакции пропорциональна продукту количеств необходимых веществ. Итак, если $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}$ приблизиться к нулю с увеличением времени.

Экспресс- Модель в MATLAB

The 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.

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

Figure contains 6 axes. Axes 1 with title y(1) contains an object of type line. Axes 2 with title y(2) contains an object of type line. Axes 3 with title y(3) contains an object of type line. Axes 4 with title y(4) contains an object of type line. Axes 5 with title y(5) contains an object of type line. Axes 6 with title y(6) contains an object of type line.

Задача оптимизации

Чтобы подготовить задачу к решению в основанном на проблеме подходе, создайте трехэлементную переменную оптимизации r который имеет нижнюю границу 0.1 и верхняя граница 10.

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

Целевой функцией для этой задачи является сумма квадратов различий между решением ОДУ с параметрами r и решение с истинными параметрами yvals. Чтобы выразить эту целевую функцию, сначала запишите функцию MATLAB, которая вычисляет решение ОДУ с использованием параметров r. Эта функция является 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);

Выразите целевую функцию как сумму квадратов различий между решением ОДУ и решением с истинными параметрами.

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

Сумма квадратов различий по существу равна нулю, что означает, что решатель нашел параметры, которые заставляют решение ОДУ соответствовать решению с истинными параметрами. Итак, как и ожидалось, решение содержит истинные параметры.

disp(rsol.r)

    2.5000
    1.2000
    0.4500

disp(rtrue)

    2.5000    1.2000    0.4500

Ограниченные наблюдения

Предположим, что вы не можете наблюдать все компоненты y, но только конечные выходы y(5) и y(6). Можете ли вы получить значения всех скоростей реакции на основе этой ограниченной информации?

Чтобы выяснить, измените функцию RtoODE для возврата только пятого и шестого выходов ОДУ. Модифицированный решатель ОДУ находится в 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

The 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

Снова, возвращенная сумма квадратов по существу равна нулю. Значит ли это, что решатель нашел правильные скорости реакции?

disp(rsol2.r)

    1.7811
    1.5730
    0.5899

disp(rtrue)

    2.5000    1.2000    0.4500

Нет; в этом случае новые ставки довольно сильно отличаются от истинных. Однако график нового решения ОДУ по сравнению с истинными значениями показывает, что 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. The axes 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}$ имеют точно такую же эволюцию как с новыми параметрами, так и с истинными параметрами.

См. также

fcn2optimexpr | ode45 | solve

Документация