FIT ODE，基于问题

打开直播脚本

此示例使用基于问题的方法，如何查找如何找到优化最小二乘法的常微分方程（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_{一世}$ 代表物质的数量 $C_{一世}$ ，然后产生反应速度 $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的快递模型

这困境函数以准备好的表格实现差分方程ODE45.。

类型困境

函数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;结尾

真正的反应率是 ${R.}_{1} = 2 。 5.$ 那 ${R.}_{2} = 1 。 2$ ，和 ${R.}_{3.} = 0. 。 4. 5.$ 。通过呼叫来计算系统的演变到零五次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）;为了i = 1：6子图（3,2，i）绘图（tspan，yvalstrue（i，i :)）标题（['y（'，num2str（i），'）']）结尾

图包含6个轴。标题Y（1）的轴1包含类型线的对象。标题Y（2）的轴2包含类型线的对象。带标题Y（3）的轴3包含类型线的对象。标题Y（4）的轴4包含类型线的对象。标题Y（5）的轴5包含类型线的对象。标题Y（6）的轴6包含类型线的对象。

优化问题

要在基于问题的方法中为解决方案做准备问题，请创建三个元素优化变量R.这具有下限0.1和一个上限10.。

r = Optimvar（'r'，3，“下行”，0.1，“上行”，10）;

此问题的目标函数是具有参数的ode解决方案之间的差异的平方和R.以及真正参数的解决方案yvals.。要表达此目标函数，请首先使用参数计算计算ode解决方案的MATLAB功能R.。这个功能是rtoode功能。

类型rtoode

功能解决方案= RTODEE（R，TSPAN，Y0）SOL = ODE45（@（t，y）diffun（t，y，r），tspan，y0）;解决= deval（sol，Tspan）;结尾

使用rtoode在一个客观函数中，通过使用将函数转换为优化表达式FCN2Optimexpr.。看将非线性函数转换为优化表达式。

myfcn = fcn2optimexpr（@ rtoode，r，tspan，y0）;

表达目标函数作为颂歌解决方案与具有真实参数的解决方案之间的平方差异之和。

obj = sum（sum（（myfcn  -  yvalstrue）。^ 2））;

使用目标函数创建优化问题obj.。

prob = OptimProblem（“客观的”，obj）;

通过呼叫来查看问题展示。

展示（prob）

优化问题：求解：R最小化：总和（（Rtooke（r，extraparams {1}） -  extraparams {3}）。^ 2,1））extraparams {1}：列1到7 00.0505 0.1010 0.1515 0.2020 0.2525 0.3030列8至14 0.3535 0.4040 0.4545 0.5051 0.5556 0.6061 0.6566列15至21 0.7071 0.7576 0.8081 0.8586 0.9091 0.9596 1.0101列22至28 1.0606 1.1111 1.1616 1.2121 1.2626 1.3131 1.3636列29至35 1.4141 1.4646 1.5152 1.5657 1.6162 1.6667 1.7172列36至42 1.7677 1.8182 1.8687 1.9192 1.9697 2.0202 2.0707列43〜49分2.1212 2.1717 2.2222 2.2727 2.3232 2.3737 2.4242列50至56分2.4747 2.5253 2.5758 2.6263 2.6768 2.7273 2.7778列57至63分2.8283 2.8788 2.9293 2.9798 3.0303 3.0808 3.1313列64至70 3.1818 3.23233.2828 3.3333 3.3838 3.4343 3.4848 3.4848列71到77 3.5354 3.5859 3.6364 3.6364 3.7374 3.7374 3.7374 3.7879 3.7879 3.8384 3.8384 3.8384 3.8384列78到84 3.8889 3.9394 3.9899 4.049 4.0909 4.1414 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919 4.1919列85H 91 4.2929 4.3939 4.4939 4.4939 4.4949 4.4949 4.5455 4.5455 4.5960 4.6970 4.6970 4.7980 4.8990 4.8990 4.8990 4.8990 4.8990 4.8990 4.8990 4.8990 4.8990柱99到100 4.9495 5.0000 extraparams {2}：1 1 0 1 0 0 Extraparams {3}：列1到7 1.0000 0.7984 0.7984 0.7984 0.7984 0.79840.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列8至14 0.5308 0.4975 0.4681 0.4420 0.4186 0.3976 0.3786 0.53080.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列15至21 0.3613 0.3456 0.3311 0.3178 0.3056 0.2942 0.2837 0.3613 0.3456 0.3311 0.3178 0.30560.2942 0.2837 0.3489 0.3444 0.3395 0.3342 0.3287 0.3230 0.3173 0.7102 0.690.0.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.求解器和呼叫解决。

R0.R = [1 1];[RSOL，SUMSQ] =求解（prob，R0）

使用LSQNONLIN解决问题。发现本地最低限度。优化完成，因为梯度的大小小于最优耐受性的值。

RSOL =结构与字段：R：[3x1双]

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.。

类型rtoode2.

功能解决方案= RTOODEE2（R，TSPAN，Y0）解决方案= RTOOKEE（R，TSPAN，Y0）;解决= solpts（[5,6]，:);％只是y（5）和y（6）结束

这rtoode2.功能只需呼叫rtoode然后采取最终两行的输出。

创建新的优化表达式rtoode2.和优化变量R.，时间跨度数据Tspan.和初始点y0.。

myfcn2 = fcn2optimexpr（@ rtoode2，r，tspan，y0）;

修改比较数据仅包括输出5和6。

YVALS2 = YVALTRUE（[5,6]，:);

从优化表达式创建一个新的目标和新的优化问题myfcn2.和比较数据Yvals2.。

obj2 = sum（sum（（myfcn2  -  yvals2）。^ 2））;prob2 = OptimProblem（“客观的”，obj2）;

根据这一有限的观察来解决问题。

[RSOL2，SUMSQ2] =求解（prob2，R0）

使用LSQNONLIN解决问题。地方最低可能。LSQNONLIN停止，因为相对于其初始值的平方和的最终变化小于功能公差的值。

RSOL2 =结构与字段：R：[3x1双]

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）匹配真实值。

图绘图（TSPAN，YVALS2（1，:)，'b-'） 抓住在ss2 = rtoode2（rsol2.r，tspan，y0）;情节（TSPAN，SS2（1，:)，'r--'）绘图（Tspan，Yvals2（2，:)，'C-'）绘图（TSPAN，SS2（2，:)，'m--'） 传奇（'真正的（5）'那'新y（5）'那'真正的Y（6）'那'新y（6）'那'地点'那'西北'） 抓住离开

图包含轴。轴包含4个类型的4个物体。这些对象表示真y（5），新y（5），true y（6），新y（6）。

要确定此问题的正确反应速率，您必须具有更多观察数据y（5）和y（6）。

用新参数绘制解决方案的所有组件，并用真正的参数绘制解决方案。

图YVALS2 = RTODEE（RSOL2.R，TSPAN，Y0）;为了i = 1：6子图（3,2，i）plot（tspan，yvalstrue（i，:)，'b-'，tspan，yvals2（我，:)，'r--'） 传奇（'真的'那'新的'那'地点'那'最好的事物'） 标题（['y（'，num2str（i），'）']）结尾

图包含6个轴。带有标题Y（1）的轴1包含2个类型的物体。这些对象表示真实，新的。带标题Y（2）的轴2包含2个类型的线。这些对象表示真实，新的。带标题Y（3）的轴3包含2个类型的线。这些对象表示真实，新的。具有标题Y（4）的轴4包含2个类型的型号。这些对象表示真实，新的。标题Y（5）的轴5包含2个类型的2个物体。 These objects represent True, New. Axes 6 with title y(6) contains 2 objects of type line. These objects represent True, New.