Huge error while using data fitting in optimizing the parameters of EOMs

4 views (last 30 days)
Dear all,
As stated in my last questions, I am doing the parameters identification/estimation for my equations of motion (EOMs).
I am using data fitting method to estimate the unknown parameters in my EOMs.
Firstly, I solve the equations using ODE45 with the guessing values of the unknown parameters and get the acceleration, velocity, and displacement data. Then I set these them as unknown parameters and give the initial values.
Finally, these values are estimated by data fitting method. However, the errors are huge it seems that the cost function can not be converged. Even, I use the same values that are used to solve the differential equations as the initial vaule for data fitting. There are still huge errors. My whole code are below:
Many thanks for your help.
Best wishes,
Yu
clear all; clc, close all
z = 14; % the distance from rotational centre to mooring line
ht = 92.61; % the distance from the centre of gravity of tower to rotational centre
hp = 14.94; % the distance from the centre of gravity of platform to rotational centre
g = 9.81;
m = 20093000; % total mass
mp = 17839000; % platform mass
It = 6.561*10^9; % mass moment of inertia of tower
Ip = 1.251*10^10;% mass moment of inertia of platform
mt = 2254000; % tower+RNA mass
ma_s = 9.64*10^6; % Added mass for platform surge
Ia = 1.16*10^10; % Added mass for platform pitch
Ia_p = -1.01*10^8; % Coupling effects of added mass bewteen surge and pitch
cs = 9.225*10^5; % damping in surge motion (x-axis translation)
ct = 6.9515*10^7; % tower bending damping
cp = 1.16*10^10; % platform pitch damping
cs_p = -8.918*10^6; % coupled damping value between surge and platform pitch
ks = 7.964*10^5; % stiffness in surge motion
kp = 2.19*10^10; % rotational stiffness of platform
kt = 1.5944*10^10; % bending stiffness of tower
tspan = 0:0.1:675;
M = [m mt*ht -mp*hp; mt*ht It+mt*ht.^2 0; -mp*hp 0 Ip+mp*hp.^2];
A = [ma_s 0 Ia_p; 0 0 0; Ia_p 0 Ia;];
C = [cs 0 cs_p; 0 ct -ct; cs_p -ct ct+cp];
K = [ks 0 -ks*z; 0 kt-mt*g*ht -kt; -ks*z -kt kp+ks*z.^2+kt+mp*g*hp];
X0 = [0 5 5 0 0 0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,X] = ode45(@(t,X) reducedmodel(t,X,M,A,C,K),tspan,X0,options);
figure,
subplot(3,1,1), plot(t,X(:,1)),grid, xlabel('time/ s'), ylabel('surge/ m')
subplot(3,1,2), plot(t,X(:,2)),grid, xlabel('time/ s'), ylabel('tower pitch/ deg')
subplot(3,1,3), plot(t,X(:,3)),grid, xlabel('time/ s'), ylabel('platform pitch/ deg')
% save the data as raw vector
x = X(:,1:3)'; % displacement
xdot = X(:,4:6)'; % velocity
xx = (M+A)\(-K*x-C*xdot); % acceleration
% get the results of acceleration
ydata = [x; xdot];
ydata(7:9,:) = xx;
% Optimization
un = zeros(7,1);
% initial value for unknown parameters
un0 = [9.225*10^5; -8.918*10^6; 6.9515*10^7; 1.16*10^10; 7.964*10^3; 1.5944*10^10; 2.19*10^10];
fun = @(un)optfun(un,X);
options = optimset('PlotFcns', @optimplotfval, 'Display', 'iter','MaxFunEvals',3000);
unbest = fminsearch(fun,un0,options);
function dXdt = reducedmodel(t,X,M,A,C,K)
x = X(1:3);
xdot = X(4:6);
xddot = (M+A)\(-K*x-C*xdot);
dXdt = [xdot; xddot];
end
function sse = optfun(un,ydata)
% parameters for turbine model
z = 14;
ht = 92.61;
hp = 14.94;
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_s = 9.64*10^6; % Added mass for platform surge
%ma_h = 2.480*10^7; % Added mass for platform heave
Ia = 1.16*10^10; % Added mass for platform pitch
Ia_p = -1.01*10^8;
%ch = 1.3*10^5;
%kh = 4.470*10^6;
% Divide the vector into displacement, velocity, and acceleration
% xddot = X(:,7:9)'; % Acceleration
% xdot = X(:,4:6)'; % Velocity
% x = X(:,1:3)'; % Displacement
disp = ydata(1:3,:);
vel = ydata(4:6,:);
acc = ydata(7:9,:);
% un(1) = cs, un(2) = cs_p, un(3) = ct, un(4) = cp, un(5) = ks, un(6) =kt
% un(7) = kp
% Mass matrix, added mass martix, damping, stiffness
M = [m mt*ht -mp*hp; mt*ht It+mt*ht.^2 0; -mp*hp 0 Ip+mp*hp^2];
A = [ma_s 0 Ia_p; 0 0 0; Ia_p 0 Ia;];
C = [un(1) 0 un(2) ; 0 un(3) -un(3); un(2) -un(3) un(3)+un(4)];
K = [un(5) 0 -un(5)*z ; 0 un(6)-mt*g*ht -un(6); -un(5)*z -un(6) un(7)+un(5)*z^2+un(6)+mp*g*hp ];
% determine the size of the martrices
% disp(size(M));
% disp(size(A));
% disp(size(C));
% disp(size(K));
% disp(size(xx));
% disp(size(xdot));
% disp(size(x));
% cost function
y_optimize = (M + A) * acc + C * vel + K * disp;
residuals = sum(y_optimize);
sse = sqrt(mean(residuals.^2));
end

Answers (0)

Categories

Find more on Oil, Gas & Petrochemical in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!