Main Content

Fit ODE, Problem-Based

This example shows how to find parameters that optimize an ordinary differential equation (ODE) in the least-squares sense, using the problem-based approach.

Problem

The problem is a multistep reaction model involving several substances, some of which react with each other to produce different substances.

For this problem, the true reaction rates are unknown. So, you need to observe the reactions and infer the rates. Assume that you can measure the substances for a set of times t. From these observations, fit the best set of reaction rates to the measurements.

Model

The model has six substances, C1 through C6, that react as follows:

  • One C1 and one C2 react to form one C3 at rate r1

  • One C3 and one C4 react to form one C5 at rate r2

  • One C3 and one C4 react to form one C6 at rate r3

The reaction rate is proportional to the product of the quantities of the required substances. So, if yi represents the quantity of substance Ci, then the reaction rate to produce C3 is r1y1y2. Similarly, the reaction rate to produce C5 is r2y3y4, and the reaction rate to produce C6 is r3y3y4.

In other words, the differential equation controlling the evolution of the system is

dydt=[-r1y1y2-r1y1y2-r2y3y4+r1y1y2-r3y3y4-r2y3y4-r3y3y4r2y3y4r3y3y4].

Start the differential equation at time 0 at the point y(0)=[1,1,0,1,0,0]. These initial values ensure that all of the substances react completely, causing C1 through C4 to approach zero as time increases.

Express Model in MATLAB

The diffun function implements the differential equations in a form ready for solution by 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

The true reaction rates are r1=2.5, r2=1.2, and r3=0.45. Compute the evolution of the system for times zero through five by calling 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

Optimization Problem

To prepare the problem for solution in the problem-based approach, create a three-element optimization variable r that has a lower bound of 0.1 and an upper bound of 10.

r = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);

The objective function for this problem is the sum of squares of the differences between the ODE solution with parameters r and the solution with the true parameters yvals. To express this objective function, first write a MATLAB function that computes the ODE solution using parameters r. This function is the RtoODE function.

type RtoODE
function solpts = RtoODE(r,tspan,y0)
sol = ode45(@(t,y)diffun(t,y,r),tspan,y0);
solpts = deval(sol,tspan);
end

To use RtoODE in an objective function, convert the function to an optimization expression by using fcn2optimexpr. See Convert Nonlinear Function to Optimization Expression.

myfcn = fcn2optimexpr(@RtoODE,r,tspan,y0);

Express the objective function as the sum of squared differences between the ODE solution and the solution with true parameters.

obj = sum(sum((myfcn - yvalstrue).^2));

Create an optimization problem with the objective function obj.

prob = optimproblem("Objective",obj);

View the problem by calling 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

Solve Problem

To find the best-fitting parameters r, give an initial guess r0 for the solver and call 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

The sum of squared differences is essentially zero, meaning the solver found parameters that cause the ODE solution to match the solution with true parameters. So, as expected, the solution contains the true parameters.

disp(rsol.r)
    2.5000
    1.2000
    0.4500
disp(rtrue)
    2.5000    1.2000    0.4500

Limited Observations

Suppose that you cannot observe all the components of y, but only the final outputs y(5) and y(6). Can you obtain the values of all the reaction rates based on this limited information?

To find out, modify the function RtoODE to return only the fifth and sixth ODE outputs. The modified ODE solver is in 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

The RtoODE2 function simply calls RtoODE and then takes the final two rows of the output.

Create a new optimization expression from RtoODE2 and the optimization variable r, the time span data tspan, and the initial point y0.

myfcn2 = fcn2optimexpr(@RtoODE2,r,tspan,y0);

Modify the comparison data to include outputs 5 and 6 only.

yvals2 = yvalstrue([5,6],:);

Create a new objective and new optimization problem from the optimization expression myfcn2 and the comparison data yvals2.

obj2 = sum(sum((myfcn2 - yvals2).^2));
prob2 = optimproblem("Objective",obj2);

Solve the problem based on this limited set of observations.

[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

Once again, the returned sum of squares is essentially zero. Does this mean that the solver found the correct reaction rates?

disp(rsol2.r)
    1.7811
    1.5730
    0.5899
disp(rtrue)
    2.5000    1.2000    0.4500

No; in this case, the new rates are quite different from the true rates. However, a plot of the new ODE solution compared to the true values shows that y(5) and y(6) match the true values.

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

To identify the correct reaction rates for this problem, you must have data for more observations than y(5) and y(6).

Plot all the components of the solution with the new parameters, and plot the solution with the true parameters.

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

With the new parameters, substances C1 and C2 drain more slowly, and substance C3 does not accumulate as much. But substances C4, C5, and C6 have exactly the same evolution with both the new parameters and the true parameters.

See Also

| |

Related Topics