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

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

[Drishti]

Can you let me know where you are encountering error in the provided code