Problen with optimization of two parameters, and with differential equations

1 次查看（过去 30 天）
Jorge Armando Vazquez
编辑： Torsten 2023-4-17
I made this program for a chemistry class, but for some reason it shows that the objective function enters optimal values but there really isn't a good similarity between the real data and the optimized data, I don't know if it's a problem of how I wrote the optimizer or because of the resolution of differentials
the data is:
0 0
60 0.28
120 0.36
180 0.42
240 0.44
300 0.445
360 0.47
420 0.5
clc
clear all
%Valores iniciales de parametros
par0 = [690,30 ];
data = load('acido.txt');
%Entrada de prueba de presion
fun_objetivo = @(par)FunObjetivo(par,data(:,:));
t = data(:,1);
xe= data(:,2);
%Argumentos de entrada
A =[]; b =[];
Aeq =[]; beq =[];
lb =[500 0]; ub =[1000 100];
IntCon =[]; nonlcon =[];
nvar =2;
%Nelder_Mead
options = optimset('MaxIter',5000,'MaxFunEvals',...
3000,'FunValCheck','off','Display','iter');
%We use here the Nelder-Mead method
par_optimos = fminsearchbnd(fun_objetivo,par0,lb,ub, options);
function FO =FunObjetivo(par,data)
%Vector tiempo
t = data(:,1);
%Vector conversion real
xe = data(:,2);
%Constantes
cc=0.365851537;
ca0=7.718472259;
keq=3.6021;
T=323.15;
k0=par(1);
Ea=par(2);
%Simulation
y0=0;
dxadt=@(t,xa) k0*exp(-Ea/(8310*T))*cc*ca0*((1-xa^2)-((xa^2)/keq));
tspan=[0:60:420];
[t,xa]=ode45(dxadt,tspan,y0);
%Valor de la Funcion Objetivo
FO = sum((xe-xa).^2);
end

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2023-4-17
移动：Torsten 2023-4-17
The parameters you want to fit are not independent.
Since T remains constant, the complete expression
k0*exp(-Ea/(8310*T))*cc*ca0
can be merged to one parameter to be fitted.
  3 个评论
显示 1更早的评论
Torsten
Torsten 2023-4-17
编辑：Torsten 2023-4-17
No. They cannot be estimated independently if T does not change during your integration.
Jorge Armando Vazquez
Jorge Armando Vazquez 2023-4-17
Ok I change it, but the problem remains, maybe is other thing?
clc
clear all
%Valores iniciales de parametros
par0 = [0.0008];
data = load('acido.txt');
%Entrada de prueba de presion
fun_objetivo = @(par)FunObjetivo(par,data(:,:));
t = data(:,1);
xe= data(:,2);
%Argumentos de entrada
A =[]; b =[];
Aeq =[]; beq =[];
lb =[0]; ub =[1];
IntCon =[]; nonlcon =[];
nvar =1;
%Nelder_Mead
options = optimset('MaxIter',5000,'MaxFunEvals',...
3000,'FunValCheck','off','Display','iter');
%We use here the Nelder-Mead method
par_optimos = fminsearchbnd(fun_objetivo,par0,lb,ub, options);
function FO =FunObjetivo(par,data)
%Vector tiempo
t = data(:,1);
%Vector conversion real
xe = data(:,2);
%Constantes
cc=0.365851537;
ca0=7.718472259;
keq=3.6021;
Te=323.15;
k1=par(1);
%Simulation
y0=0;
dxadt=@(t,xa) k1*cc*ca0*((1-xa^2)-((xa^2)/keq));
%dxadt=@(t,xa) k0*exp(-Ea/(8310*Te))*cc*ca0*((1-xa^2)-((xa^2)/keq));
tspan=[0:60:420];
[t,xa]=ode45(dxadt,tspan,y0);
%Valor de la Funcion Objetivo
FO = sum((xe-xa).^2);
end

请先登录，再进行评论。

Torsten
Torsten 2023-4-17
编辑：Torsten 2023-4-17
Ran in:
All curves that stem from your model tend to sqrt(keq/(1+keq)). Since your measurement data tend to 0.5, my guess is that you have to choose keq such that 0.5 = sqrt(keq/(1+keq)) to get a proper fit. The constant you try to fit gives the velocity with which this equilibrium value of 0.5 is reached.
syms c t keq xa(t)
eqn = diff(xa,t)==c*((1-xa^2)-xa^2/keq);
cond = xa(0)==0;
sol = dsolve(eqn,cond)
sol = 
keqnum = double(solve(sqrt(keq/(keq+1))==0.5,keq))
keqnum = 0.3333
%keqnum=3.6021;
sol = matlabFunction(subs(sol,keq,keqnum));
data=[0 0
60 0.28
120 0.36
180 0.42
240 0.44
300 0.445
360 0.47
420 0.5];
t = data(:,1);
xe = data(:,2);
c = 0:0.001:0.01;
hold on
plot(t,xe,'o')
for i = 1:numel(c)
xa = sol(c(i),t);
plot(t,xa)
end
hold off

类别

在 Help Center 和 File Exchange 中查找有关 Quadratic Programming and Cone Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by