Этот пример показывает, как найти параметры, которые оптимизируют обыкновенное дифференциальное уравнение (ОДУ) в смысле наименьших квадратов, используя основанный на проблеме подход.
Проблема заключается в многоступенчатой модели реакции, включающей несколько веществ, некоторые из которых реагируют друг с другом с образованием различных веществ.
Для этой проблемы истинные скорости реакции неизвестны. Итак, нужно наблюдать реакции и выводить скорости. Предположим, что вы можете измерить вещества в течение ряда раз . Из этих наблюдений, соответствует лучшему набору скоростей реакции измерениям.
Модель имеет шесть веществ, через , которые реагируют следующим образом:
Один и один реагировать, чтобы сформировать один по курсу
Один и один реагировать, чтобы сформировать один по курсу
Один и один реагировать, чтобы сформировать один по курсу
Скорость реакции пропорциональна продукту количеств необходимых веществ. Итак, если представляет количество вещества , затем скорость реакции для получения является . Точно так же скорость реакции для получения является , и скорость реакции для получения является .
Другими словами, дифференциальное уравнение, управляющее эволюцией системы,
Запустите дифференциальное уравнение в момент 0 в точке . Эти начальные значения гарантируют, что все вещества реагируют полностью, вызывая через приблизиться к нулю с увеличением времени.
The diffun функция реализует дифференциальные уравнения в форме, готовой к решению ode45.
type diffunfunction 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
Истинные скорости реакции: , , и . Вычислите эволюцию системы на время от нуля до пяти с помощью вызова 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
Чтобы подготовить задачу к решению в основанном на проблеме подходе, создайте трехэлементную переменную оптимизации r который имеет нижнюю границу 0.1 и верхняя граница 10.
r = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);
Целевой функцией для этой задачи является сумма квадратов различий между решением ОДУ с параметрами r и решение с истинными параметрами yvals. Чтобы выразить эту целевую функцию, сначала запишите функцию MATLAB, которая вычисляет решение ОДУ с использованием параметров r. Эта функция является RtoODE функция.
type RtoODEfunction 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 RtoODE2function 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
Чтобы идентифицировать правильные скорости реакции для этой проблемы, вы должны иметь данные для большего количества наблюдений, чем 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
С новыми параметрами, веществами и слить медленнее, и вещество не накапливается так много. Но вещества , , и имеют точно такую же эволюцию как с новыми параметрами, так и с истинными параметрами.
fcn2optimexpr | ode45 | solve