La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

Ajuste de parámetros EDO mediante variables de optimización

Este ejemplo muestra cómo encontrar parámetros que optimicen una ecuación diferencial ordinaria (EDO) mediante mínimos cuadrados, utilizando variables de optimización (enfoque basado en problemas).

Problema

El problema es un modelo de reacción multipaso en el que intervienen varias sustancias, algunas de las cuales reaccionan entre sí para producir sustancias diferentes.

En el caso de este problema, se desconocen las velocidades de reacción verdaderas. Por lo tanto, hay que observar las reacciones e inferir las velocidades. Suponga que puede medir las sustancias para un conjunto de tiempos $$t$$. A partir de estas observaciones, ajuste el mejor conjunto de velocidades de reacción a las mediciones.

Modelo

El modelo tiene seis sustancias, de $$C_1$$ a $$C_6$$, que reaccionan de la siguiente forma:

• Un $$C_1$$ y un $$C_2$$ reaccionan para formar un $$C_3$$ a la velocidad $$r_1$$

• Un $$C_3$$ y un $$C_4$$ reaccionan para formar un $$C_5$$ a la velocidad $$r_2$$

• Un $$C_3$$ y un $$C_4$$ reaccionan para formar un $$C_6$$ a la velocidad $$r_3$$

La velocidad de reacción es proporcional al producto de las cantidades de las sustancias necesarias. Por tanto, si $$y_i$$ representa la cantidad de sustancia $$C_i$$, la velocidad de reacción para producir $$C_3$$ es $$r_1{y}_1{y}_2$$. De forma similar, la velocidad de reacción para producir $$C_5$$ es $$r_2{y}_3{y}_4$$ y la velocidad de reacción para producir $$C_6$$ es $$r_3{y}_3{y}_4$$.

En otras palabras, la ecuación diferencial que controla la evolución del sistema es

$$\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].$$

Inicie la ecuación diferencial en el tiempo 0 y en el punto $$y(0)=[1,1,0,1,0,0]$$. Estos valores iniciales garantizan que todas las sustancias reaccionen completamente, haciendo que $$C_1$$ a través de $$C_4$$ se aproxime a cero a medida que aumenta el tiempo.

Expresar un modelo en MATLAB

La función diffun implementa las ecuaciones diferenciales en una forma que está lista para ser resuelta por 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

Las velocidades de reacción verdaderas son $$r_1=2.5$$, $$r_2=1.2$$ y $$r_3=0.45$$. Calcule la evolución del sistema para los tiempos de cero a cinco llamando a 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

Problema de optimización

Para preparar el problema con el fin de solucionarlo en el enfoque basado en problemas, cree una variable de optimización de tres elementos r que tenga un límite inferior de 0.1 y un límite superior de 10.

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

La función objetivo de este problema es la suma de cuadrados de las diferencias entre la solución de la EDO con los parámetros r y la solución con los parámetros verdaderos yvals. Para expresar esta función objetivo, escriba primero una función de MATLAB que calcule la solución de la EDO utilizando parámetros r. Esta función es la función RtoODE.

type RtoODE

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

Para usar RtoODE en una función objetivo, convierta la función en una expresión de optimización con fcn2optimexpr. Consulte Convertir una función no lineal en una expresión de optimización.

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

Exprese la función objetivo como la suma de las diferencias al cuadrado entre la solución de la EDO y la solución con los parámetros verdaderos.

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

Cree un problema de optimización con la función objetivo obj.

prob = optimproblem("Objective",obj);

Resolver el problema

Para encontrar los parámetros r que mejor se ajustan, dé una estimación inicial r0 para el solver y llame a 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.8659e-15

La suma de las diferencias al cuadrado es esencialmente cero, lo que significa que el solver encontró parámetros que hacen que la solución de la EDO coincida con la solución con parámetros verdaderos. Así que, como era de esperar, la solución contiene los parámetros verdaderos.

disp(rsol.r)

    2.5000
    1.2000
    0.4500

disp(rtrue)

    2.5000    1.2000    0.4500

Observaciones limitadas

Suponga que no se pueden observar todos los componentes de y, sino solo las salidas finales y(5) e y(6). ¿Puede obtener los valores de todas las velocidades de reacción a partir de esta información limitada?

Para averiguarlo, modifique la función RtoODE para que devuelva solo las salidas quinta y sexta de la EDO. El solver de la EDO que se ha modificado está en 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

La función RtoODE2 llama a RtoODE y luego toma las dos últimas filas de la salida.

Cree una nueva expresión de optimización a partir de RtoODE2 y la variable de optimización r, los datos del intervalo de tiempo tspan y el punto inicial y0.

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

Modifique los datos de comparación para incluir solo las salidas 5 y 6.

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

Cree un nuevo objetivo y un nuevo problema de optimización a partir de la expresión de optimización myfcn2 y los datos de comparación yvals2.

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

Resuelva el problema basándose en este conjunto limitado de observaciones.

[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

Una vez más, la suma de cuadrados devuelta es esencialmente cero. ¿Significa esto que el solver ha encontrado las velocidades de reacción correctas?

disp(rsol2.r)

    1.7811
    1.5730
    0.5899

disp(rtrue)

    2.5000    1.2000    0.4500

No, en este caso, las nuevas velocidades son muy diferentes de las verdaderas. Sin embargo, una gráfica que compara la nueva solución de la EDO con los valores reales muestra que y(5) e y(6) coinciden con los valores reales.

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

Para identificar las velocidades de reacción correctas en este problema, debe disponer de datos para más observaciones que y(5) e y(6).

Represente todos los componentes de la solución con los nuevos parámetros y represente la solución con los parámetros reales.

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

Con los nuevos parámetros, las sustancias $$C_1$$ y $$C_2$$ drenan más lentamente y la sustancia $$C_3$$ no se acumula tanto. Sin embargo, las sustancias $$C_4$$, $$C_5$$ y $$C_6$$ tienen exactamente la misma evolución tanto con los nuevos parámetros como con los parámetros verdaderos.

Consulte también

solve | ode45 | fcn2optimexpr

Temas relacionados

• Flujo de trabajo de optimización basada en problemas
• Discretized Optimal Trajectory, Problem-Based