Clear Filters
Clear Filters

Plot frequency responses based on the ODE results

17 views (last 30 days)
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

Accepted Answer

Star Strider
Star Strider on 23 May 2024
Edited: Star Strider on 23 May 2024 at 22:02
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’.
.
  14 Comments
Yu
Yu on 28 May 2024 at 11:17
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
Star Strider on 28 May 2024 at 13:22
@Yu — I have looked for it, however I have not been able to find it yet. Please share the URL for it.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!