Plot frequency responses based on the ODE results

Hello all,

I would like to plot the frequency responses based on the ode results to see the behaviors of the system. I have already established the state-space model for my ode equations. Additionally, I use ss2tf to get the transfer fucntion. Now I am to plot the frequency responses of all degrees of freedom on one figure. However, 'freqz' does not work in my case.

I checked the details from turtorial. unfortunately, I failed. I hope you can help me with this. My codes can be seen in the following. I appreciate your help.

clear; clc;
global all_F
global ex_F
all_F = [];
ex_F = [];

syms z_T

h_f = 150;
h_R = 144.582; % m the height from the MSL to tower top;
H_T = 129.582; % m tower height from the tower bottom;
h_T = 15; % m height from tower base to tower bottom;
h = 29.94;
h_1 = 164.94;
z = 14.94; % m the distance from rotational centre to mooring line
h_t = 92.61;  % m the distance from the centre of gravity of tower to rotational centre
h_p = 14.94; % m the distance from the centre of gravity of platform to rotational centre
g = 9.81;
m = 20093000;  % kg /total mass
m_T = 1.263e5; % 8.6e5; % kg / tower
m_N = 1.017e6; % kg / nacelle
m_p = 1.7838e7;  % kg /platform mass
xi_TA = 0.01;
I_p = 1.2507*10^10;% kg m2 /mass moment of inertia of platform
m_as = 9.64*10^6;  % kg / Added mass for platform surge
I_a = 1.16*10^10;  % kg m2 / Added mass for platform pitch
m_asp = -1.01*10^8; % kg m / Coupling effects of added mass bewteen surge and pitch
c_s = 7.94e4;  % N s2/m2 / damping in surge motion (x-axis translation)
c_sp = -2.25e5; % coupled damping value between surge and platform pitch
c_p = 5.60e8; % damping in pitch motionS 5
k_p = 2.190e9; % rotational stiffness of platform
k_mooringS = 7.965e4;
k_mooringSP = 1.162e6;
k_mooringP = 2.65e8;

tspan = 0:0.025:200;

% definition of TMD parameters


X0 = [0 0 10*pi/180 0 0 0]'; % initial pitch motion

% ==================== definition of the tower properties===============

mu = 0.0084*z_T^3-1.077*z_T^2-171.5*z_T+2.94e04; % mass per length

EI = 1.905e06*z_T^3-2.47e08*z_T^2-5.208e10*z_T+6.851e12; % tower bending stiffness

Phi_TA = 0.9414*(z_T/H_T)^2+0.3468*(z_T/H_T)^3-1.073*(z_T/H_T)^4+1.3139*(z_T/H_T)^5-0.5289*(z_T/H_T)^6; % tower fore-aft first mode shape

% mass component

fun1 = mu*Phi_TA^2; m_TA = double(int(fun1,z_T,0,H_T));

fun2 = mu*Phi_TA; m_1 = double(int(fun2,z_T,0,H_T));

fun3 = mu*(z_T)*Phi_TA; m_2 = double(int(fun3,z_T,0,H_T));

fun4 = mu*(z_T); m_3 = double(int(fun4,z_T,0,H_T));

fun5 = mu*(z_T)^2; m_4 = double(int(fun5,z_T,0,H_T));

% fun6 = mu; m_T = double(int(fun6,z_T,0,H_T));

% Stiffness of the tower

D2y = diff(Phi_TA,z_T,2); Dy = diff(Phi_TA,z_T,1);

fun6 = EI*D2y^2; f1 = int(fun6,0,H_T);
fun7 = mu; f2 = int(fun7,z_T,H_T);
fun8 = g*(m_N+f2)*Dy^2; f3 = int(fun8,0,H_T);
k_TA = double(f1-f3);

% ================= definition of matrices ======================

M = [m_N+m_TA       m_N+m_1                        m_N*h_R+m_2;
     m_N+m_1        m_N+m_p+m_T                   (m_N*h_R-m_p*h_p+m_3+m_T*h_T);
     m_N*h_R+m_2    (m_N*h_R-m_p*h_p+m_3)          m_N*h_R^2+I_p+m_4];
C = [2*xi_TA*sqrt(m_TA*k_TA) 0 0;
     0 0 0;
     0 0 0];
K = [k_TA 0 -(m_N+m_T)*g;
     0 0 0;
     -(m_N+m_T)*g 0 -(m_N*h_R+m_3-m_p*h_p)*g];
M_add = [0 0 0;
         0 m_as m_asp;
         0 m_asp I_a];
C_add = [0 0 0;           % problem_
         0 c_s c_sp;
         0 c_sp c_p];
K_add = [0 0 0;
         0 0 0;
         0 0 k_p];
K_mooring = [0 0 0;      
             0 k_mooringS k_mooringSP;
             0 k_mooringSP k_mooringP];

M_1 = M+M_add;
K_1 = K+K_mooring+K_add;
C_1 = C+C_add;

% state space model--------------

O = zeros(3,3);
O1 = zeros(3,1);
N = ones(3,1);
I = eye(3);
A = [O I;-inv(M_1)*K_1 -inv(M_1)*C_1];
B = [O O; O inv(M_1)];
E = [I O; O O];
D = zeros(6);



% the forces and moments are extracted from Orcaflex 


options = odeset('RelTol',1e-10,'AbsTol',1e-10);

[t,X] = ode45(@(t,X) reducedmodel(t,A,B,X),tspan,X0,options);

[b,a] = ss2tf(A,B,E,D,1);

PtfmPitch_deg = X(:,3)*180/pi;

figure, 

subplot(3,1,1), plot(t,X(:,1)),grid, xlabel('time/ s'), ylabel('TTDspFA/ m')
subplot(3,1,2), plot(t,X(:,2)),grid, xlabel('time/ s'), ylabel('surge/ m')
subplot(3,1,3), plot(t,PtfmPitch_deg),grid, xlabel('time/ s'), ylabel('platform pitch/ deg')

function dXdt = reducedmodel(t,A,B,X)

coes = 9.225e5;
coesp = 8.92e6;
coep = 1.68e10;

F1 = -(coes*X(5)*abs(X(5))-coesp*X(6)*abs(X(6)));
F2 = -(coep*X(6)*abs(X(6))-coesp*X(5)*abs(X(5)));
  
F=[0;0;0;0;F1;F2];

dXdt = A*X+B*F;

disp(t)

end

Best wishes,

Yu

Answer 1

[Star Strider]

I am not certain what result you want. The freqz function only works on Signal Ptocessing Toolbox filter results. To use the Control System Toolbox to get the transfer function plots, do this instead:

sys = ss(A,B,E,D);

figure

bodeplot(sys)

grid

I added that at the end of my code here. (It is not going to be straightforward to interpret those results, at least for me.) I also added my one-sided Fourier transform function ‘FFT1’ since it may be convenient to do element-wise division of its outputs (complex) to determine whatever transfer function you want from the results of integrating your differential equations. I used it to calculate the Fourier transforms of your results, and then plotted them.

I am not certain what result you want, however this may get you started —

clear; clc;
global all_F
global ex_F
all_F = [];
ex_F = [];

syms z_T

h_f = 150;
h_R = 144.582; % m the height from the MSL to tower top;
H_T = 129.582; % m tower height from the tower bottom;
h_T = 15; % m height from tower base to tower bottom;
h = 29.94;
h_1 = 164.94;
z = 14.94; % m the distance from rotational centre to mooring line
h_t = 92.61;  % m the distance from the centre of gravity of tower to rotational centre
h_p = 14.94; % m the distance from the centre of gravity of platform to rotational centre
g = 9.81;
m = 20093000;  % kg /total mass
m_T = 1.263e5; % 8.6e5; % kg / tower
m_N = 1.017e6; % kg / nacelle
m_p = 1.7838e7;  % kg /platform mass
xi_TA = 0.01;
I_p = 1.2507*10^10;% kg m2 /mass moment of inertia of platform
m_as = 9.64*10^6;  % kg / Added mass for platform surge
I_a = 1.16*10^10;  % kg m2 / Added mass for platform pitch
m_asp = -1.01*10^8; % kg m / Coupling effects of added mass bewteen surge and pitch
c_s = 7.94e4;  % N s2/m2 / damping in surge motion (x-axis translation)
c_sp = -2.25e5; % coupled damping value between surge and platform pitch
c_p = 5.60e8; % damping in pitch motionS 5
k_p = 2.190e9; % rotational stiffness of platform
k_mooringS = 7.965e4;
k_mooringSP = 1.162e6;
k_mooringP = 2.65e8;

tspan = 0:0.025:200;

% definition of TMD parameters


X0 = [0 0 10*pi/180 0 0 0]'; % initial pitch motion

% ==================== definition of the tower properties===============

mu = 0.0084*z_T^3-1.077*z_T^2-171.5*z_T+2.94e04; % mass per length

EI = 1.905e06*z_T^3-2.47e08*z_T^2-5.208e10*z_T+6.851e12; % tower bending stiffness

Phi_TA = 0.9414*(z_T/H_T)^2+0.3468*(z_T/H_T)^3-1.073*(z_T/H_T)^4+1.3139*(z_T/H_T)^5-0.5289*(z_T/H_T)^6; % tower fore-aft first mode shape

% mass component

fun1 = mu*Phi_TA^2; m_TA = double(int(fun1,z_T,0,H_T));

fun2 = mu*Phi_TA; m_1 = double(int(fun2,z_T,0,H_T));

fun3 = mu*(z_T)*Phi_TA; m_2 = double(int(fun3,z_T,0,H_T));

fun4 = mu*(z_T); m_3 = double(int(fun4,z_T,0,H_T));

fun5 = mu*(z_T)^2; m_4 = double(int(fun5,z_T,0,H_T));

% fun6 = mu; m_T = double(int(fun6,z_T,0,H_T));

% Stiffness of the tower

D2y = diff(Phi_TA,z_T,2); Dy = diff(Phi_TA,z_T,1);

fun6 = EI*D2y^2; f1 = int(fun6,0,H_T);
fun7 = mu; f2 = int(fun7,z_T,H_T);
fun8 = g*(m_N+f2)*Dy^2; f3 = int(fun8,0,H_T);
k_TA = double(f1-f3);

% ================= definition of matrices ======================

M = [m_N+m_TA       m_N+m_1                        m_N*h_R+m_2;
     m_N+m_1        m_N+m_p+m_T                   (m_N*h_R-m_p*h_p+m_3+m_T*h_T);
     m_N*h_R+m_2    (m_N*h_R-m_p*h_p+m_3)          m_N*h_R^2+I_p+m_4];
C = [2*xi_TA*sqrt(m_TA*k_TA) 0 0;
     0 0 0;
     0 0 0];
K = [k_TA 0 -(m_N+m_T)*g;
     0 0 0;
     -(m_N+m_T)*g 0 -(m_N*h_R+m_3-m_p*h_p)*g];
M_add = [0 0 0;
         0 m_as m_asp;
         0 m_asp I_a];
C_add = [0 0 0;           % problem_
         0 c_s c_sp;
         0 c_sp c_p];
K_add = [0 0 0;
         0 0 0;
         0 0 k_p];
K_mooring = [0 0 0;      
             0 k_mooringS k_mooringSP;
             0 k_mooringSP k_mooringP];

M_1 = M+M_add;
K_1 = K+K_mooring+K_add;
C_1 = C+C_add;

% state space model--------------

O = zeros(3,3);
O1 = zeros(3,1);
N = ones(3,1);
I = eye(3);
A = [O I;-inv(M_1)*K_1 -inv(M_1)*C_1];
B = [O O; O inv(M_1)];
E = [I O; O O];
D = zeros(6);



% the forces and moments are extracted from Orcaflex 


options = odeset('RelTol',1e-10,'AbsTol',1e-10);

[t,X] = ode45(@(t,X) reducedmodel(t,A,B,X),tspan,X0,options);

[b,a] = ss2tf(A,B,E,D,1);

PtfmPitch_deg = X(:,3)*180/pi;

figure, 

subplot(3,1,1), plot(t,X(:,1)),grid, xlabel('time/ s'), ylabel('TTDspFA/ m')
subplot(3,1,2), plot(t,X(:,2)),grid, xlabel('time/ s'), ylabel('surge/ m')
subplot(3,1,3), plot(t,PtfmPitch_deg),grid, xlabel('time/ s'), ylabel('platform pitch/ deg')

[FTx1,Fv1] = FFT1(X(:,1),t);
[FTx2,Fv2] = FFT1(X(:,2),t);
[FTPitch,Fv3] = FFT1(PtfmPitch_deg,t);

figure
tiledlayout(3,1)
nexttile
plot(Fv1, abs(FTx1)*2)
grid
set(gca, 'YScale','log')
ylabel('Amplitude')
title('X_1')
nexttile
plot(Fv1, abs(FTx2)*2)
grid
set(gca, 'YScale','log')
ylabel('Amplitude')
title('X_2')
nexttile
plot(Fv1, abs(FTx1)*2)
grid
set(gca, 'YScale','log')
xlabel('Frequency')
ylabel('Amplitude')
title('PtfmPitch\_deg')
sgtitle('Fourier Transforms')

sys = ss(A,B,E,D,1);

figure
bodeplot(sys)
grid

function dXdt = reducedmodel(t,A,B,X)

coes = 9.225e5;
coesp = 8.92e6;
coep = 1.68e10;

F1 = -(coes*X(5)*abs(X(5))-coesp*X(6)*abs(X(6)));
F2 = -(coep*X(6)*abs(X(6))-coesp*X(5)*abs(X(5)));
  
F=[0;0;0;0;F1;F2];

dXdt = A*X+B*F;

% disp(t)

end


function [FTs1,Fv] = FFT1(s,t)
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

EDIT — (23 May 2024 at 22:02)

Changed ‘Fourier Transform’ figure xlabel to correctly label it as ‘Frequency’.

.

[Yu]

Dear Star,

Many thanks for your quick reply.

Actually, I was trying to plot the frequency responses based on the time series results. However, there are some errors when using 'freqz' function after extracting the transfer function with 'ss2tf' function. Therefore, I used the 'fft' function instead. But the results seems to be not correct. As the frequency of the peak is known already.

Following your codes, the x-axis represents the time. But I prefer frequency which could be helpful to analyze the frequency responses.

Best wishes,

Yu

[Yu]

Dear Star,

If you mean the x-axis represents the frequency. Can the frequency results be this form? All peaks of three degrees of freedom in a same figure.

In addition, I am not sure whether there is the same peak which is near 0.5 Hz in the frequency response plots. If we refering the time-series results, it may be not consistent.

Best wishes,

Yu

[Star Strider]

My pleasure!

Labeling it as ‘Time’ was simply an error (I was not paying sufficiently close attention to what I was doing this morning when I wrote my answer). I edited it to correctly label it as ‘Frequency’.

Rather than extracting the transfer function coefficients from the state space realisation, use the bode or bodeplot functions to plot it directly, as I did. It may be necessary to select the desired input and output vectors to only plot the ones you want. (It plotted the magnitude and phase results for 6 inputs and 6 outputs when I plotted it, so you may want to choose only the ones you are interested in.)

Other than that, you can calculate the transfer functions from the complex Fourier transforms by dividing the output by the input element-wise (using the ./ right division operator). Then plot that result, similarly to the way I plotted the Fourier transform results. I am not certain what those vectors are, so I did not do that division.

[Yu]

Hi Star,

Thank you for your clarification and suggestions.

I still have one question regarding your approach. I used the simple fft to plot the frequency response (seen below) for these three vectors (X1, X2, X3 in your figure). The results are consistent with the physical system. The frequencies of the peaks are obvious different.

I have checked your codes following your explanation. But I am confused about your results now. I found there is a peak existing around 0.5 Hz in all three plots in your results. Do I need to take some action to eliminate this? Or you just put all the peaks in one plot? Additionally, I am thinking what's the meaning and difference of the y-axis between the two approach (complex Fourier transform and FFT). It would also be better if there are some reference you can attach.

Best wishes,

Yu

[Star Strider]

My plesure!

The data are yours. All I did was to provide an alternative approach to calculating and plotting the Fourier transforms, and suggested that using bode or bodeplot is the correct approach to using the identified system to plot the frequency responses.

‘I found there is a peak existing around 0.5 Hz in all three plots in your results. Do I need to take some action to eliminate this?’

That peak exists because it is in your data. I do not understand the reason you would want to eliminate it.

.

[Yu]

Hi Star,

Thank you for your reply.

As I stated before, I believe the results in my last reply (simple results with fft) is convinced. The frequencies of the peak are close to the values in the technical report. In other words, the peaks around 0.5 Hz in second and third plots should not exist.

I am not sure where is the problem as you said the codes you provide is just a approach to get my expected results.

Best wishes,

Yu

[Yu]

Hi Star,

Additionally, could you please add some comments on your codes especially for the calculation of fft using hanning window. It would be better to help me understand.

Best wishes,

Yu

[Star Strider]

I am not certain what you are doing, or the reason the frequencies at 0.5 Hz should not exist. These are your data. If there is an error in your code that creates that peak, you need to find it and correct it. I cannot do that.

[Yu]

Thank you for your reply. I used the hamming window to replace hanning window. The results become more expected.

Besides, my data are from the differential equations. I also tried this process using the data from the software following your approach and direct fft approach. They are similar with the results from my equations.

Yes, I will try to find the reasons. There are always many things to learn.

Best wishes,

Yu

[Star Strider]

My pleasure!

I do not understand what you are doing, however there may be an error in your differential equations. One way to check that would be to use bode or bodeplot and look at each input-output pair individually. That might offer some insight as to what the problem is, or where it is.

The window used for the fft should not make much difference (except for the rectangular window, which is the same as no window).

[Yu]

Dear Star,

Sorry for the late reply. The results from my equations are similar with the results from the FEM software under comparison.

Thank you for your help.

Best wishes,

Yu

[Star Strider]

As always, my pleasure!

No worries! I am happy that everything worked out and you resolved these issues.

[Yu]

Dear Star,

I am not sure if you are available for another question about the conditional statement. I used some conditions in ODE function to generate some external excitations based on the 'real-time' displacement of the system. However, I found it is hard to move to next time step for ode solver.

I am setting a new question abou this issue. It would be good if you can take a look at it.

Best wishes,

Yu

[Star Strider]

@Yu — I have looked for it, however I have not been able to find it yet. Please share the URL for it.