Main Content

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 C1 a C6, que reaccionan de la siguiente forma:

  • Un C1 y un C2 reaccionan para formar un C3 a la velocidad r1

  • Un C3 y un C4 reaccionan para formar un C5 a la velocidad r2

  • Un C3 y un C4 reaccionan para formar un C6 a la velocidad r3

La velocidad de reacción es proporcional al producto de las cantidades de las sustancias necesarias. Por tanto, si yi representa la cantidad de sustancia Ci, la velocidad de reacción para producir C3 es r1y1y2. De forma similar, la velocidad de reacción para producir C5 es r2y3y4 y la velocidad de reacción para producir C6 es r3y3y4.

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

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

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 C1 a través de C4 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 r1=2.5, r2=1.2 y r3=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

Figure contains 6 axes objects. Axes object 1 with title y(1) contains an object of type line. Axes object 2 with title y(2) contains an object of type line. Axes object 3 with title y(3) contains an object of type line. Axes object 4 with title y(4) contains an object of type line. Axes object 5 with title y(5) contains an object of type line. Axes object 6 with title y(6) contains an object of type line.

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

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

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

Figure contains 6 axes objects. Axes object 1 with title y(1) contains 2 objects of type line. These objects represent True, New. Axes object 2 with title y(2) contains 2 objects of type line. These objects represent True, New. Axes object 3 with title y(3) contains 2 objects of type line. These objects represent True, New. Axes object 4 with title y(4) contains 2 objects of type line. These objects represent True, New. Axes object 5 with title y(5) contains 2 objects of type line. These objects represent True, New. Axes object 6 with title y(6) contains 2 objects of type line. These objects represent True, New.

Con los nuevos parámetros, las sustancias C1 y C2 drenan más lentamente y la sustancia C3 no se acumula tanto. Sin embargo, las sustancias C4, C5 y C6 tienen exactamente la misma evolución tanto con los nuevos parámetros como con los parámetros verdaderos.

Consulte también

| |

Temas relacionados