Этот пример показывает, как найти параметры, которые оптимизируют обычное дифференциальное уравнение (ОДУ) в смысле наименьших квадратов, используя подход, основанный на задачах.
Проблема заключается в многоступенчатой модели реакции, включающей несколько веществ, некоторые из которых взаимодействуют друг с другом для получения различных веществ.
Для этой проблемы истинные скорости реакции неизвестны. Итак, нужно наблюдать за реакциями и выводить скорости. Предположим, что вы можете измерить вещества для набора раз T. Из этих наблюдений, подогнать наилучший набор скорости реакции к измерениям.
Модель содержит шесть веществ, через , которые реагируют следующим образом:
Один и один реагируют с образованием одного со скоростью
Один и один реагируют с образованием одного со скоростью
Один и один реагируют с образованием одного со скоростью
Скорость реакции пропорциональна количеству требуемых веществ. Так, если представляет количество вещества , то скорость реакции для получения равна . Аналогично, скорость реакции для получения равна , а скорость реакции для получения равна .
Другими словами, дифференциальное уравнение, управляющее эволюцией системы,
r1y1y2-r3y3y4-r2y3y4-r3y3y4y4r3y4].
Начните дифференциальное уравнение в момент времени 0 в точке 1,1,0,1,0,0]. Эти исходные значения обеспечивают полную реакцию всех веществ, C1 C4 приближаться к нулю по мере увеличения времени.
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
Истинными скоростями реакции являются 2,51,2 = 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
Для подготовки задачи к решению в подходе, основанном на проблеме, создайте трехэлементную переменную оптимизации. r которая имеет нижнюю границу 0.1 и верхняя граница 10.
r = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);
Целевой функцией для этой задачи является сумма квадратов разностей между решением ОДУ с параметрами r и решение с истинными параметрами yvals. Чтобы выразить эту целевую функцию, сначала напишите функцию MATLAB, которая вычисляет решение ODE с помощью параметров 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
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