适合歌唱参数使用优化变量

Open Live Script

This example shows how to find parameters that optimize an ordinary differential equation (ODE) in the least-squares sense, using optimization variables (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, $C_{1}$ through $C_{6}$ , that react as follows:

One $C_{1}$ and one $C_{2}$ react to form one $C_{3}$ at rate $r_{1}$
One $C_{3}$ and one $C_{4}$ react to form one $C_{5}$ at rate $r_{2}$
One $C_{3}$ and one $C_{4}$ react to form one $C_{6}$ at rate $r_{3}$

The reaction rate is proportional to the product of the quantities of the required substances. So, if $y_{i}$ represents the quantity of substance $C_{i}$ , then the reaction rate to produce $C_{3}$ is $r_{1} y_{1} y_{2}$ ．Similarly, the reaction rate to produce $C_{5}$ is $r_{2} y_{3} y_{4}$ , and the reaction rate to produce $C_{6}$ is $r_{3} y_{3} y_{4}$ ．

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

$\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}] ．$

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 $C_{1}$ through $C_{4}$ to approach zero as time increases.

Express Model in MATLAB

Thediffunfunction implements the differential equations in a form ready for solution byode45．

typediffun

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 $r_{1} = 2 ． 5$ , $r_{2} = 1 ． 2$ , and $r_{3} = 0 ． 45$ ．Compute the evolution of the system for times zero through five by callingode45．

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);fori = 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 variablerthat has a lower bound of0.1and an upper bound of10．

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 parametersrand the solution with the true parametersyvals．To express this objective function, first write a MATLAB function that computes the ODE solution using parametersr．This function is theRtoODEfunction.

typeRtoODE

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

To useRtoODE在一个目标函数,转换函数an optimization expression by usingfcn2optimexpr．SeeConvert 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 functionobj．

prob = optimproblem("Objective",obj);

Solve Problem

To find the best-fitting parametersr, give an initial guessr0for the solver and callsolve．

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.8659e-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

假设你不能遵守所有的components ofy, but only the final outputsy(5)andy(6)．Can you obtain the values of all the reaction rates based on this limited information?

To find out, modify the functionRtoODEto return only the fifth and sixth ODE outputs. The modified ODE solver is inRtoODE2．

typeRtoODE2

function solpts = RtoODE2(r,tspan,y0) solpts = RtoODE(r,tspan,y0); solpts = solpts([5,6],:); % Just y(5) and y(6) end

TheRtoODE2function simply callsRtoODEand then takes the final two rows of the output.

Create a new optimization expression fromRtoODE2and the optimization variabler, the time span datatspan, and the initial pointy0．

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 expressionmyfcn2and the comparison datayvals2．

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 thaty(5)andy(6)match the true values.

figure plot(tspan,yvals2(1,:),'b-') holdonss2 = 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') holdoff

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent True y(5), New y(5), True y(6), New y(6).

To identify the correct reaction rates for this problem, you must have data for more observations thany(5)andy(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);fori = 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 $C_{1}$ and $C_{2}$ drain more slowly, and substance $C_{3}$ does not accumulate as much. But substances $C_{4}$ , $C_{5}$ , and $C_{6}$ have exactly the same evolution with both the new parameters and the true parameters.