time for parameter estimation using ode 45

2022 10 19
1 个回答
更新时间：2022 11 7
2 次查看（30 天）

0 个投票

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 ?

11 个评论
显示 9更早的评论 隐藏 9更早的评论

Sam Chak
Sam Chak 2022-10-19
Hi @nado, without seeing the code, we probably cannot provide any meaningful advice.
nado
nado 2022-10-19
[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
Walter Roberson 2022-10-19
Long search time would be expected for stiff systems
Torsten
Torsten 2022-10-19
编辑：Torsten 2022-10-19
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
Bjorn Gustavsson 2022-10-19
@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
nado 2022-10-19
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
nado 2022-10-19
Thank you @Bjorn Gustavsson
I will do the modifications you suggested
Sam Chak
Sam Chak 2022-10-19
@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
nado 2022-10-19
@Torsten
I tried ode15s and it works !!!
Thanks a lot !!!

请先登录，再进行评论。

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

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

提问:

nado
nado
2022-10-19

回答：

Rajiv Singh
Rajiv Singh
2022-11-7

Community Treasure Hunt

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

Start Hunting!

Translated by