Estimate and optimize the parameters in EOMs

Question

Yu 2024-3-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2093556-estimate-and-optimize-the-parameters-in-eoms

编辑： Yu 2024-3-13

Dear Sir or Madam,

I hope you are all doing well. I am trying to estimate and optimize the parameters with the numerical data (Displacement, Velocity, Acceleration) in my equations (Matrixes). And there are totally 7 parameters needed to estimate.

I have the structure's motion responses from the software (could be assumed as experimental data). I would like to use these data to estimate and optimize some parameters in the equations as some parameters are nonliner and hard to express mathematically.

I was refering to the 'Non-linear data fit' of help center. However, I have no idea about how to apply the optimization approach to my case. From my perspectives, I need to transfer my equations of motion to first-order system. I assume the external function is 0 (free-decay test). Therefore, the cost function should be 'MyFunction' - 0 = 0. My function could be seen as follows:

I would much appreciate your support.

Thank you.

Best wishes,

Yu Gao

% X1 = Surge; X2 = PlatformPitch; X3 = RelativeAngularDiplacement; X4=Heave
% X5 = SurgeVelocity; X6 = PlatformPitchVelocity; 
% X7 = RelativeAugularVelocity; X8 = HeaveVelocity
function yu_gao
function kineticsmodel(un,t)
X0 = [0;10;10;0;0;0;0;0]; % [x10;x20;x30;x40; dx1dt0;dx2dt0;dx3dt0;dx4dt0]
[t, Xv] = ode45(@odefn,t,X0);
    function dxdt = odefn(t,X)
    
    x = X(1:4);
    xdot = X(5:8);
    
    M = [m mt*ht Iac 0;...
         mt*ht It 0 0;...
         Iac 0 Ip+mp*hp^2 0;...
         0 0 0 m];
    
    A = [ma 0 -mp*hp 0;...
        0 mt*ht^2 0 0;...
        -mp*hp 0 Ia 0;...
        0 0 0 mh];
    
    C = [un(4) 0 un(7) 0;...
         0 un(5) -un(5) 0;...
         un(7) -un(5) un(5)+un(6) 0;...
         0 0 0 ch];
    
    K = [un(1) 0 -un(1)*z 0;...
         0 un(2)-mt*ht -un(2) 0;...
         -un(1)*z -un(2) un(3)+un(1)*z^2+un(2)+mp*g*hp 0;...
         0 0 0 kh];
    xddot = (M+A)\(-K*x-C*xdot);
    dXdt = [xdot xddot];
    end
X=Xv;
end
tspan = 0:0.1:600;
X0 = [0;10;10;0;0;0;0;0];
X = [Surge PlatfromPitch VarName17 Heave SurgeVelocity PlatformPitchsVelocity VarName11 HeaveVelocity];
% X1 = Surge; X2 = PlatformPitch; X3 = RelativeAngularDiplacement; X4=Heave
% X5 = SurgeVelocity; X6 = PlatformPitchVelocity; 
% X7 = RelativeAugularVelocity; X8 = HeaveVelocity
z = 14; 
ht = 56.50;
hp = 14.94;
height_t = 129.13;
g = 9.81;
m = 20093000; % total mass
mp = 17839000; % platform mass
It = 6.561*10^9;
Ip = 1.251*10^10;
mt = 2254000; % tower+RNA mass
ma = 9.64*10^6; % Added mass for platform surge
mh = 2.480*10^7; % Added mass for platform heave
Ia = 1.16*10^10; % Added mass for platform pitch
Iac = -1.01*10^8;
ch = 1.3*10^5;
kh = 4.470*10^6;
% initial value of the unknown parameters.
un0 = [7.964*10^3 1.5944*10^10 4.453*10^6 9.225*10^5 6.9515*10^7 1.676*10^9 -8.918*10^6];
[un,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kineticsmodel,un0,t,X);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, un(k1))
end
tv = linspace(min(tspan), max(tspan));
unfit = kineticsmodel(un, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, unfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'un_1(t)', 'un_2(t)', 'un_3(t)', 'un_4(t)', 'un_5(t)','un_6(t)','un_7(t)','Location','N')
end

8 个评论
显示 6更早的评论隐藏 6更早的评论

Yu 2024-3-13

移动：Star Strider 2024-3-13

在 MATLAB Online 中打开

Optimization Data.xlsx

Dear Star Strider,

I appreciate your detailed explanation. I was trying to use your script to reproduce my case. However, I found my data can not be read when the 'X' is within the 'function'. Therefore, I can not check the potential problem based on your example. My data and code have been attached.

The error is:

Unrecognized function or variable 'Surge'.

Error in yu_gao (line 44)

X = [Surge PlatfromPitch VarName17 Heave SurgeVelocity PlatformPitchsVelocity VarName11 HeaveVelocity];

Best wishes,

Yu

function yu_gao
function kineticsmodel(un,t)
X0 = [0;10;10;0;0;0;0;0]; % [x10;x20;x30;x40; dx1dt0;dx2dt0;dx3dt0;dx4dt0]
[t, Xv] = ode45(@odefn,t,X0);
    function dxdt = odefn(t,X)
    
    x = X(1:4);
    xdot = X(5:8);
    
    M = [m mt*ht Iac 0;...
         mt*ht It 0 0;...
         Iac 0 Ip+mp*hp^2 0;...
         0 0 0 m];
    
    A = [ma 0 -mp*hp 0;...
        0 mt*ht^2 0 0;...
        -mp*hp 0 Ia 0;...
        0 0 0 mh];
    
    C = [un(4) 0 un(7) 0;...
         0 un(5) -un(5) 0;...
         un(7) -un(5) un(5)+un(6) 0;...
         0 0 0 ch];
    
    K = [un(1) 0 -un(1)*z 0;...
         0 un(2)-mt*ht -un(2) 0;...
         -un(1)*z -un(2) un(3)+un(1)*z^2+un(2)+mp*g*hp 0;...
         0 0 0 kh];
    xddot = (M+A)\(-K*x-C*xdot);
    dXdt = [xdot xddot];
    end
X=Xv;
end
tspan = 0:0.1:675;
X0 = [0;10;10;0;0;0;0;0];
X = [Surge PlatfromPitch VarName17 Heave SurgeVelocity PlatformPitchsVelocity VarName11 HeaveVelocity];
z = 14; 
ht = 56.50;
hp = 14.94;
height_t = 129.13;
g = 9.81;
m = 20093000; % total mass
mp = 17839000; % platform mass
It = 6.561*10^9;
Ip = 1.251*10^10;
mt = 2254000; % tower+RNA mass
ma = 9.64*10^6; % Added mass for platform surge
mh = 2.480*10^7; % Added mass for platform heave
Ia = 1.16*10^10; % Added mass for platform pitch
Iac = -1.01*10^8;
ch = 1.3*10^5;
kh = 4.470*10^6;
un0 = [7.964*10^3 1.5944*10^10 4.453*10^6 9.225*10^5 6.9515*10^7 1.676*10^9 -8.918*10^6];
[un,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kineticsmodel,un0,t,X);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, un(k1))
end
tv = linspace(min(tspan), max(tspan));
unfit = kineticsmodel(un, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, unfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'un_1(t)', 'un_2(t)', 'un_3(t)', 'un_4(t)', 'un_5(t)','un_6(t)','un_7(t)','Location','N')
end

Yu 2024-3-13

移动：Star Strider 2024-3-13

Optimization Data.xlsx

I am sorry for the misleading message. There are three parts of this excel.

The firts part is the acceleration of different degrees of freedom of the structure.

The second part is the velocity of the structure.

The third part is the displacement of the structure.

This is the wind turbine mathematical model. Surge motion is the translational motion along X-axis. Heave motion is the translational motion along Z-axis. I assume there is a relative movement between the tower and platform. Therefore, the PlatformPitch represents the platform's rotational angles around the Y-axis. Also, the TowerPitch is the rotational angles around the Y-axis. However, there is a difference between these two motions. So the relative motion can be represented by using Nacelle's motion (connected to the tower) minus Platform motion as shown in the table.

In conclusion, you can only focus on the RelativePitchAngle, Surge, Heave (Coluimn R, S, T), RelativePitchAngularVelocity, SurgeVelocity, HeaveVelociy (Column K, L, M)...

What I want to do is that I would like use these known results from the software (numerical results) to estimate the unknown parameters. I suppose it should be called 'offline' method as I do not need to calculate the matrixes. By the way, the equations can probably not be transfered to the differential form as there are many coupled terms. Then, I could got the true intial vaules of these parameters. After that, I could use 'Online' method to get the accurate values.

I reupload the file. Please refer to it. Many thanks for your help.

Best wishes,

Yu

Sam Chak 2024-3-13

Also, instead of labeling technical names in the Excel, try labeling them to indicate which column data corresponds to X1, X2, X3, and so on up to X8.

By the way, it seems that the x-axis data {-75 to 600} is not aligned with the rest of the columns, likely due to a hasty copy-paste. It is uncommon for software-generated data to have this kind of mistake.

Yu 2024-3-13

编辑：Yu 2024-3-13

Optimization Data.xlsx

Hi @Sam Chak Thank you for your suggestions. I have added some comments in the codes. Please see the updated codes at the top.

Yes, there are 7 unknown parameters (possible more). And I want to find the parameters to satisfy my equations by using numerical or experimental results.

The vector X should comprise velocity, and displacment. I am not sure whether I should use acceleration as it is the second-order differential. As you can see in my kinetic model. They are multipled by mass, damping, and stiffness respectively. Yes, I would use these parameters from experiment or numerical simulation to find the optimal parameters. To be honest, it should be basic step. This is the 'offline' method (I was told). Because I know the variables' values. I can introduce these values into the matrices (Md/dt(dx/dt)+C(dx/dt)+Kx=0). Let F=Md/dt(dx/dt)+C(dx/dt)+Kx. So the finally target is find the 7 unknow parameters to let F=0.

Thank you for your sugestions. The minus time is becasue it is static analysis which is distinguished by it. I have modified it. Please refer the file attached.

Best wishes,

Yu

请先登录，再进行评论。

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

Estimate and optimize the parameters in EOMs

8 个评论
显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Estimate and optimize the parameters in EOMs

8 个评论 显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论