В этом примере показано, как найти параметры, которые оптимизируют обыкновенное дифференциальное уравнение (ODE) в смысле наименьших квадратов, с помощью подхода, основанного на проблеме.
Проблемой является многоступенчатая модель реакции включение нескольких веществ, некоторые из которых реагируют друг с другом, чтобы произвести различные вещества.
Для этой проблемы истинные скорости реакции неизвестны. Так, необходимо наблюдать реакции и вывести уровни. Примите, что можно измерить вещества для набора времен . От этих наблюдений соответствуйте лучшему набору скоростей реакции к измерениям.
Модель имеет шесть веществ, через , это реагирует можно следующим образом:
Один и один реагируйте, чтобы сформировать тот на уровне
Один и один реагируйте, чтобы сформировать тот на уровне
Один и один реагируйте, чтобы сформировать тот на уровне
Скорость реакции пропорциональна продукту количеств необходимых веществ. Так, если представляет количество вещества , затем скорость реакции, чтобы произвести . Точно так же скорость реакции, чтобы произвести , и скорость реакции, чтобы произвести .
Другими словами, дифференциальное уравнение, управляющее эволюцией системы,
Запустите дифференциальное уравнение во время 0 в точке . Эти начальные значения гарантируют, что все вещества реагируют полностью, вызывая через чтобы приблизиться к нулю как к, времени увеличивается.
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.8668e-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 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