# 使用优化变量拟合 ODE 参数

### 模型

• 一个 ${C}_{1}$ 和一个 ${C}_{2}$ 以速率 ${r}_{1}$ 反应形成一个 ${C}_{3}$

• 一个 ${C}_{3}$ 和一个 ${C}_{4}$ 以速率 ${r}_{2}$ 反应形成一个 ${C}_{5}$

• 一个 ${C}_{3}$ 和一个 ${C}_{4}$ 以速率 ${r}_{3}$ 反应形成一个 ${C}_{6}$

`$\frac{\mathit{dy}}{\mathit{dt}}=\left[\begin{array}{c}-{\mathit{r}}_{1}{\mathit{y}}_{1}{\mathit{y}}_{2}\\ -{\mathit{r}}_{1}{\mathit{y}}_{1}{\mathit{y}}_{2}\\ -{\mathit{r}}_{2}{\mathit{y}}_{3}{\mathit{y}}_{4}+{\mathit{r}}_{1}{\mathit{y}}_{1}{\mathit{y}}_{2}-{\mathit{r}}_{3}{\mathit{y}}_{3}{\mathit{y}}_{4}\\ -{\mathit{r}}_{2}{\mathit{y}}_{3}{\mathit{y}}_{4}-{\mathit{r}}_{3}{\mathit{y}}_{3}{\mathit{y}}_{4}\\ {\mathit{r}}_{2}{\mathit{y}}_{3}{\mathit{y}}_{4}\\ {\mathit{r}}_{3}{\mathit{y}}_{3}{\mathit{y}}_{4}\end{array}\right].$`

### 在 MATLAB 中表示模型

`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 ```

```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 = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);`

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

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

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

`prob = optimproblem("Objective",obj);`

### 求解问题

```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.8661e-15 ```

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

### 有限的观测值

`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`，然后获取输出的最后两行。

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

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

```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 ```

```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```

```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```