Estimate and optimize the parameters in EOMs

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

[Yu]

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]

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

[Star Strider]

Your problem appears to be outside my areas of expertise. I thought you wanted to fit differential equations to data, and I can help you with that if I understand the problem well enough. Apparently you are not fitting differential equations, and I do not understand the system you want to fit to your data. There are others here who probably do.

I have moved your comments to my answer to comments to your original question so that others can refer to them.

[Yu]

@Star Strider Sorry, the column headers of Acceleration, Velocity, and Displacement are to distinguish the different data groups.

[Yu]

@Star Strider Hi Star Strider,

Thank you for your help. I thought it is basically same. My equation is Md/dt(dx/dt)+C(dx/dt)+Kx=0. Now, the indenpent variables are known from other software. But some parameters in the matrices are unknown. Therefore, I would like to estimate these parapmeters to let the results of the equation to be 0.

[Sam Chak]

Hi @Yu,

Let's keep it simple. In your system of differential equations, 'odefn', you have a state vector X that comprises 8 state variables. It would be helpful if you could specify which state corresponds to the column data you want the ODEs to fit.

Additionally, I noticed there are 7 unknown parameters in the 'un' vector. Is your task focused on finding the optimal values in 'un' such that the desired state response aligns with one of the experimental data points?

[Sam Chak]

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]

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