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

Скрипт Open Live Script

В этом примере показано, как найти параметры, которые оптимизируют обыкновенное дифференциальное уравнение (ODE) в смысле наименьших квадратов, с помощью подхода, основанного на проблеме.

Проблема

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

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

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 возвратить только пятый и шестой 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

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

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

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