time for parameter estimation using ode 45

Hello everyone,

I used lsqnonlin to find unknown parameters of an ODE. However, when i run the script, it is taking a long time. Until now, it tooks 25 min and it's the first time i encounter this !

The system i'm solvin consists of 3x3 matrix.

Do you think the time for this is reasonable or there is something going wrong ???

I searched a lot on how to know whether the script is stuck somewhere but i didn't succeed.

Based on your experience, is there any method to know where the code is stuck ?

[Sam Chak]

Hi @nado, without seeing the code, we probably cannot provide any meaningful advice.

[nado]

[Ca1,resnorm,residual]=lsqnonlin(@(Ca1)odefit_calib_damping(Ca1),1e+5)

%% lsqnonlin is the function used to optimize the error between experiment and numerical model

function [ err ] = odefit_calib_damping(Ca1 )

time=0:0.1:630;

[t,y]=ode45(@(t,y)calib_damping(t,y,Ca1),time,[10;0;0;0;0;0]);

err=abs(y_exp_surge-y(:,1)).^2; %% y_exp is a vector containing data from experiment

end

function [dydt] = calib_damping( t,y,Ca1)

% this is the function where the system of coupled ODEs are implemented

M(1,1)=8.3163e+7;

M(1,2)=-2.6982e+9;

M(1,3)=1.8107e+6;

M(2,2)=5.6542e+11;

M(2,3)=3.4131e+8;

M(3,3)=1.4249e+6;

M(2,1)=M(1,2);

M(3,1)=M(1,3);

M(3,2)=M(2,3);

K(1,1)=687024.52188410972265443647316524;

K(1,2)=-4.861e+7;

K(1,3)=0;

K(2,1)=-4.861e+7;

K(2,2)=1.4604e+10-1467627741.2664321859450536919855;

K(2,3)=0;

K(3,1)=0;

K(3,2)=0;

K(3,3)=1.4249e+6;

C(1,1)=1.8448e+5 + Ca1;

C(2,2)=1.8286e+10;

C(3,3)=4.0008e+10;

C(1,2)=5.8081e+7;

C(1,3)=-8.5911e+7;

C(2,3)=-2.7048e+10;

C(2,1)=C(1,2);

C(3,1)=C(1,3);

C(3,2)=C(2,3);

M_K=-(M^(-1))*K;

M_C=-(M^(-1))*C;

dydt(1)= y(4);

dydt(2)=y(5);

dydt(3)=y(6);

dydt(4)= M_K(1,1)*y(1)+M_K(1,2)*y(2)+M_K(1,3)*y(3)+M_C(1,1)*y(4)+M_C(1,2)*y(5)+M_C(1,3)*y(6);

dydt(5)= M_K(2,1)*y(1)+M_K(2,2)*y(2)+M_K(2,3)*y(3)+M_C(2,1)*y(4)+M_C(2,2)*y(5)+M_C(2,3)*y(6);

dydt(6)= M_K(3,1)*y(1)+M_K(3,2)*y(2)+M_K(3,3)*y(3)+M_C(3,1)*y(4)+M_C(3,2)*y(5)+M_C(3,3)*y(6);

dydt=[dydt(1);dydt(2);dydt(3);dydt(4);dydt(5);dydt(6)];

end

[Walter Roberson]

Long search time would be expected for stiff systems

[Torsten]

We don't know your data vector, so we cannot run your code.

But at least change the line

err=abs(y_exp_surge-y(:,1)).^2; %% y_exp is a vector containing data from experiment

to

err=y_exp_surge-y(:,1); %% y_exp is a vector containing data from experiment

"lsqnonlin" will add the squares of the vector components of "err" internally.

Working with a mass matrix M to solve

M*y'' = K*y + C*y'

instead of multiplying by M^(-1) might help.

Also changing the ODE solver from ODE45 to ODE15S is worth a try.

[Bjorn Gustavsson]

@nado: In addition to @Torsten's comment you might have a problem if your y_exp_surge variable is a row-array. Then the

err=y_exp_surge-y(:,1); %%

would turn err into a numel(time)-by-numel(time) array - which is not what you want. That would correspond to a "very strange fitting problem". My suggestion is that your modified line should be:

err = y_exp_surge(:) - y(:,1); %% just to make sure that you get the residual for each time-step

[nado]

Thank you @Torsten

I will try what you have say. The reason that i did not put y_exp is because it is so big. Its length is 6301.

[nado]

Thank you @Bjorn Gustavsson

I will do the modifications you suggested

[nado]

Thank you @Walter Roberson

[Sam Chak]

@nado, looks like system identification problem. Your mass matrix and the damping matrix contain some very large values. So, this might explain why the code takes a long time to run.

I have heard of some people used the non-dimensional approach (through normalization). Perhaps, it may reduce the time taken by lsqnonin to solve the problem.

[nado]

Thank you @Sam Chak

[nado]

@Torsten

I tried ode15s and it works !!!

Thanks a lot !!!

Answer 1

[Rajiv Singh]

See also

Estimate Coefficients of ODEs to Fit Given Solution - MATLAB & Simulink (mathworks.com)