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