Solving a non linear ODE with unknown parameter

Hello ! I am working on solving an ODE equation with an unknown kinetic parameter A. I have been using python and deep learning to solve the equation and also determine the value of A , however the loss function is always in the order of 10**4 and the paramter A is wrong , I tried with different hyperparamters but it´s not working. this is the ODE equation : dDP/dt=-k1*([DP]^2) and k1=k= Ae^(1/R(-E/(T+273))) , A is in the order of 10**8, I have DP(t) data.

I am stuck and I would like to know what´s the best way to solve this using matlab ? or is there any examples similar to my problem ?

Any help is highly appreciated !

Answer 1

[Torsten]

%time points
ts=[1 2 3 4 5 6 7 8];
DP=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
p0 = 1e1;
p = fminunc(@(p)fun(p,ts,DP),p0)
E = 111e3;
R = 8.314;
T = 371;
A = p*exp(E/(R*T))
plot(ts,DP)
hold on
plot(ts,1./(1/DP(1)+ A*exp(-E/(R*T))*(ts-ts(1))));

function obj = fun(p,ts,DP)
  DP_model = 1./(1/DP(1)+ p*(ts-ts(1)));   
  obj = sum((DP-DP_model).^2)
end

[khaoula Oueslati]

@Torsten hello ,

I tried with fminunc and this is the output obj =

10.3642

Local minimum found.

Optimization completed because the size of the gradient is less than

the value of the optimality tolerance.

<stopping criteria details>

p =

1.8688e+10

A =

7.9484e+25

the value of p is different from the ground_truth value and also the loss is big. I am confused about certain parts in the code you proposed , how is the objective function calculated exactly ? why in the DP_model , we use the same DP (the measured ones) : DP_model = 1./(1/DP(1)+ p*(ts-ts(1))); and with p in this line of code , you mean A ?

[khaoula Oueslati]

I also dont understand why i keep getting erros when trying to solve the ODE with matlab solver like ODE15s and ODE45 ...

is there an example using Lsqcurvefit please ?

[Torsten]

As I already wrote in a previous response, the solution of your ODE

dDp/dt = -k1*Dp^2, Dp(t0) = Dp0

can explicitly be given as

Dp = 1/(1/Dp0 + k1*(t-t0))

I chose t0 = 1 and Dp(t0) = 1000.

This explains the line

DP_model = 1./(1/DP(1)+ p*(ts-ts(1)));

Since exp(-E/(R*T)) is constant, I chose

p = A*exp(-E/(R*T))

as the parameter to be fitted.

So when fminunc returns p, I have to multiply by exp(E/(R*T)) to get A:

A = p*exp(E/(R*T))

I don't know what code you use, but I get

p = 5.0841e-04
A = 2.1624e+12

[khaoula Oueslati]

Thank you @Torsten for your answer and for your valuable help ! I really appreciate it.

I tried with your code and it gave me the same output as yours , and I tried with another dataset for which I know the A value, this is the code:

% data for regression
DP = [1000	995.6179	991.2741	986.968	982.6991	978.467	974.2712	970.1113	965.9867	961.897	957.8419	953.8207	949.8332	945.8789	941.9574	938.0683	934.2111	930.3855	926.5912	922.8276	919.0945,...
    915.3915	911.7182	908.0743	904.4594	900.8732	897.3152	893.7853	890.283	886.8081	883.3602	879.939	876.5442	873.1755	869.8326	866.5152,...
    863.223	859.9557	856.7131	853.4948	850.3006	847.1302	843.9834	840.8599	837.7594	834.6817	831.6265	,...
    828.5936	825.5827	822.5937	819.6262	816.68	813.755	810.8508	807.9673	805.1042	802.2614	799.4385	796.6354	793.852	791.0879 ,...
    788.343	785.6171	782.9099	780.2214	777.5512	774.8993	772.2654	769.6493	767.0509	764.47	761.9064	759.36	756.8305	754.3178,...
    751.8217	749.3421	746.8788	744.4316	742.0005	739.5851	737.1854	734.8013	732.4325	730.0789	727.7405	725.4169	723.1082,...
    720.814	718.5344	716.2692	714.0182	711.7813	709.5584	707.3494	705.154	702.9722	700.8039	698.6489	696.5072	694.3785	692.2628,...
    690.1599	688.0698	685.9923	683.9273	681.8748	679.8345	677.8063	675.7902	673.7861	671.7939	669.8134	667.8445	665.8872,...
    663.9413	662.0067	660.0834	658.1713	656.2701	654.38	652.5007	650.6321	648.7742	646.9269	645.0901	643.2637	641.4476];
ts=linspace(1,128,128);
E = 111e3;
R = 8.314;
T = 371;
% ts=[1 2 3 4 5 6 7 8];
% DP=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
plot(ts,DP)
ylabel('DP')
xlabel('t')
hold on
p0 = 1e1;
p = fminunc(@(p)fun(p,ts,DP),p0);
disp(['loss value is ' num2str(fun(p,ts,DP))])
disp(['A value is ' num2str(A)])
plot(ts,1./(1/DP(1)+ p*(ts-ts(1))),'bo');
%the value I am looking for
A = p*exp(E/(R*T));
%The objective function
function obj = fun(p,ts,DP)
  DP_model = 1./(1/DP(1)+ p*(ts-ts(1)));   
  obj = sum((DP-DP_model).^2);
end

the loss value is :loss value is 0.35787 and A value is 1.08e10 and the ground_truth A value is 7.8e8

I am not sure the source of this mismatch.

btw , how did you find the DP_model expression ? is it some appoximation ? or after integration we get that expression of the solution ?

[Torsten]

the loss value is :loss value is 0.35787 and A value is 1.08e10 and the ground_truth A value is 7.8e8

I am not sure the source of this mismatch.

I don't know either. Maybe T or E were different. The fit at least is perfect.

btw , how did you find the DP_model expression ? is it some appoximation ? or after integration we get that expression of the solution ?

If you don't trust in my pencil-and-paper solution, here is MATLAB code to solve the differential equation:

syms Dp(t) k1 t0 Dp0
eqn = diff(Dp,t) == -k1*Dp^2;
cond = Dp(t0) == Dp0;
DpSol(t) = dsolve(eqn,cond)

[khaoula Oueslati]

Thanks a lot @Torsten ! it's not that I don't trust in the expression , I just wanted to know how to find it so that I can explain it to my academic supervisor. About the mismatch , it's probably due to some variation in the data or E value , I will further look into that and thanks again !

Answer 2

[Torsten]

Your ODE for D_p gives

D_p = 1/(1/D_p0 + k1*(t-t0))

where D_p0 = D_p(t0).

Now you can apply "lsqcurvefit" to fit the unknown parameter A.

[khaoula Oueslati]

Thanks Torsten for your help! I will try it !

[khaoula Oueslati]

Hello @Torsten , I didnt know how to implement using Lsqcurvefit and since I am not really familiar with matlab , I am struggling to use it for my problem. I tried to use fminunc but erros keep coming up and I didnt know how to tackle them.

this is my code using fminunc

function dDP= Mymodel(t,DP,A)
%known parameters
E=111e3;
R=8.314;
T=371;
dDP= -(A*(exp(-E/R*T)))*(DP^2);

end
function obj = objective(A)
DP0=1000;
%time points
ts=[1 2 3 4 5 6 7 8];
[t,DP]= ode45(@(t,DP)Mymodel(t,DP,A),ts,DP0);
DP_measured=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];


A=(DP-DP_measured).^2;
obj=sum(A);

end
A0=1e8;
fun = @objective;
[A,fval]=fminunc(fun,A0);
disp(['A:' num2str(A)])

These errors come up:
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in objective (line 5)
[t,DP]= ode45(@(t,DP)Mymodel(t,DP,A),ts,DP0);
 
optimization_DP
Error using fminunc
Supplied objective function must return a scalar value.

Error in optimization_DP (line 3)
[A,fval]=fminunc(fun,A0);
 

Answer 3

[Sam Chak]

Hi @khaoula Oueslati

This governing equations are given and you have acquired the data.

The objective is want to find A.

From the data, you can possibly estimate for . Next, can be determined from the differential equation:

Now, if R, E and T are known, then can be determined from the algebraic equation:

Please verify this.

If the data is uniformly distributed, then you can use this method to estimate .

t = -pi:(2*pi/100):pi;
x = sin(t);                     % assume Dp is a sine wave
y = gradient(x)/(2*pi/100);     % estimate dotDp, a cosine wave is expected
plot(t, x, 'linewidth', 1.5, t, y, 'linewidth', 1.5)
grid on
xlabel('t')
ylabel('x(t) and x''(t)')
legend('x(t) = sin(t)', 'x''(t) = cos(t)', 'location', 'northwest')

[khaoula Oueslati]

Thanks a lot Sam ! I will try it, but the thing is , I need to be able to find that parameter using an optimization algorithm and not from an algebraic equation and I tried with python but I didnt figure it out.

Answer 4

[David Willingham]

Hi,

Have you seen this example for solving ODE's using Deep Learning in MATLAB?

Dynamical System Modeling Using Neural ODE

[khaoula Oueslati]

Thanks David ! No I haven't seen it before but I will definitely check it !