Use of dde23
This topic has been permanently closed and transferred to MATLAB Answers.
Hey guys! I'm undergraduate in aerospace engineering at Brazil and I has done researchs in aeroelasticity fields in the last years.
Recently, I start to use dde23 for solve problems involved time-delay and active control. However, in my last code, I found some problems and I'd like some help for understand what is going wrong. Below, It's the code I'm using.
clc
close all
clear
%% Definicao das variaveis
% Caracteristicas do escoamento
rho = 1; % densidade do escoamento
c0 = 1.0; c1 = 0.165; c2 = 0.0455; c3 = 0.335; c4 = 0.3; % coeficientes de Wagner
omega_c = 17.83; % frequencia critica de flutter
Uc = 26.8903;
U = 1.20 * Uc;
% Caracteristicas do aerofolio
a = -0.15; % distancia entre eixo de referencia e o CE
b = 0.5; % semicorda
m = 20; % massa do aerofolio
H_alpha = 7; % nao-linearidade estrutural - hardening
omega_h = 2*pi; % frequencia natural de plunge
omega_alpha = 6*pi; % frequencia natural de pitch
x_alpha = 0.25; % distancia entre CG e CE
r_alpha = 0.75; % raio de giracao do aerofolio
% Caracteristicas do nonlinear energy sink
delta = 0.75; % posicao adimensional
gamma_n = 5e5; % rigidez adimensional
mi_n = 0.10; % razao de massa
lambda_n_c = (sqrt(3) / 3) * mi_n * omega_c;
lambda_n = 0.5 * lambda_n_c; % amortecimento adimensional
G = 1;
%% Definindo a solucao do sistema
tspan = 0:0.1:120; % vetor temporal
history = [0; deg2rad(1); 0; 0; 0; 0; 0; 0]; % historico para t < 0
delays = [1; 1; 1; 1]; % definicao do delay para as variaveis: xi, alpha, up e xb
ddefun = @(t, y, Z) dydt(t, y, Z, m, a, b, r_alpha, x_alpha, omega_h, omega_alpha, H_alpha, G, mi_n, lambda_n, gamma_n, delta, rho, U, c0, c1, c2, c3, c4);
sol = dde23(ddefun, delays, history, tspan);
%%
function dydt = dydt(t, y, Z, m, a, b, r_alpha, x_alpha, omega_h, omega_alpha, H_alpha, G, mi_n, lambda_n, gamma_n, delta, rho, U, c0, c1, c2, c3, c4)
% transformacao das quatro equacoes que regem a dinamica do sistema
% em oito equacoes - variaveis de estados - com o objetivo de reduzir a
% ordem do sistema:
% xi = y(1); alpha = y(2); up = y(3); xb = y(4);
% xi_d = y(5); alpha_d = y(6); up_d = y(7); xb_d = y(8)
xi_delay = Z(:,1);
alpha_delay = Z(:,2);
up_delay = Z(:,3);
xb_delay = Z(:,4);
y_tau = [xi_delay; alpha_delay; up_delay; xb_delay];
M1 = [1, x_alpha, 0, 0; x_alpha, r_alpha^2, 0, 0; mi_n, -delta*mi_n, -mi_n, 0; 0, 0, 0, 1];
M2 = (pi*rho*b^2/m) * [-1, a, 0, 0; 0, -(1/8+a^2), 0, 0; 0, 0, 0, 0; 0, 0, 0, 0];
M = M1 - M2;
C1 = [0, 0, lambda_n, 0; 0, 0, -delta*lambda_n, 0; 0, 0, -lambda_n, 0; -1, -(1/2-a), 0, (c2+c4)*U/b];
C2 = [2*pi*rho*U*b^2*(c0-c1-c3)*(-1/(m*b)), (pi*rho*b^3*(U/b)+2*pi*rho*U*b^2*(c0-c1-c3)*(1/2-a))*(-1/(m*b)), 0, 2*pi*rho*U*b^2*(U/b)*(c1*c2+c3*c4)*(-1/(m*b));
2*pi*rho*U*b^3*(c0-c1-c3)*(1/2+a)*(1/(m*b^2)), (-pi*rho*b^4*(1/2-a)*(U/b)+2*pi*rho*U*b^3*(1/2+a)*(c0-c1-c3)*(1/2-a))*(1/(m*b^2)), 0, 2*pi*rho*U*b^3*(U/b)*(c1*c2+c3*c4)*(1/2+a)*(1/(m*b^2));
0, 0, 0, 0; 0, 0, 0, 0];
C = C1 - C2;
K1 = [omega_h^2, 0, gamma_n*y(3)^2, 0; 0, (omega_alpha*r_alpha)^2*(1+H_alpha*y(2)^2), -delta*gamma_n*y(3)^2, 0;
0, 0, -gamma_n*y(3)^2, 0; 0, -(U/b), 0, c2*c4*(U/b)];
K2 = [G/m, (2*pi*rho*U*b^2*(c0-c1-c3)*(U/b)*(-1/(m*b)))-(G/m*delta), -G/m, 2*pi*rho*U*b^2*(U/b)^2*(c2*c4)*(c1+c3)*(-1/(m*b));
-G/m*delta, (2*pi*rho*U*b^3*(c0-c1-c3)*(1/2+a)*(U/b)*(1/(m*b^2)))+(G/m*delta^2), G/m*delta, 2*pi*rho*U*b^3*(U/b)^2*(1/2+a)*(c2*c4)*(c1+c3)*(1/(m*b^2));
-G/m, G/m*delta, G/m, 0; 0, 0, 0, 0];
K = K1 - K2;
T = [-G/m, G/m*delta, G/m, 0; G/m*delta, -G/m*delta^2, -G/m*delta, 0; G/m, -G/m*delta, -G/m, 0; 0, 0, 0, 0];
N = zeros(4, 4);
I = eye(4);
A = [N, I; M\-K, M\-C];
B = [N; M\T];
dydt = A*y + B*y_tau;
end
4 Comments
Time Descending
@João Pedro please see the following change that is done in your code
clc
close all
clear
%% Definicao das variaveis
% Caracteristicas do escoamento
rho = 1; % densidade do escoamento
c0 = 1.0; c1 = 0.165; c2 = 0.0455; c3 = 0.335; c4 = 0.3; % coeficientes de Wagner
omega_c = 17.83; % frequencia critica de flutter
Uc = 26.8903;
U = 1.20 * Uc;
% Caracteristicas do aerofolio
a = -0.15; % distancia entre eixo de referencia e o CE
b = 0.5; % semicorda
m = 20; % massa do aerofolio
H_alpha = 7; % nao-linearidade estrutural - hardening
omega_h = 2*pi; % frequencia natural de plunge
omega_alpha = 6*pi; % frequencia natural de pitch
x_alpha = 0.25; % distancia entre CG e CE
r_alpha = 0.75; % raio de giracao do aerofolio
% Caracteristicas do nonlinear energy sink
delta = 0.75; % posicao adimensional
gamma_n = 5e5; % rigidez adimensional
mi_n = 0.10; % razao de massa
lambda_n_c = (sqrt(3) / 3) * mi_n * omega_c;
lambda_n = 0.5 * lambda_n_c; % amortecimento adimensional
G = 1;
%% Definindo a solucao do sistema
tspan = 0:0.1:120; % vetor temporal
history = [0; deg2rad(1); 0; 0; 0; 0; 0; 0]; % historico para t < 0
delays = [1; 1; 1; 1]; % definicao do delay para as variaveis: xi, alpha, up e xb
ddefun = @(t, y, Z) dydt(t, y, Z, m, a, b, r_alpha, x_alpha, omega_h, omega_alpha, H_alpha, G, mi_n, lambda_n, gamma_n, delta, rho, U, c0, c1, c2, c3, c4);
sol = dde23(ddefun, delays, history, tspan)
%%
function dydt = dydt(t, y, Z, m, a, b, r_alpha, x_alpha, omega_h, omega_alpha, H_alpha, G, mi_n, lambda_n, gamma_n, delta, rho, U, c0, c1, c2, c3, c4)
% transformacao das quatro equacoes que regem a dinamica do sistema
% em oito equacoes - variaveis de estados - com o objetivo de reduzir a
% ordem do sistema:
% xi = y(1); alpha = y(2); up = y(3); xb = y(4);
% xi_d = y(5); alpha_d = y(6); up_d = y(7); xb_d = y(8)
xi_delay = Z(:,1);
alpha_delay = Z(:,2);
up_delay = Z(:,3);
xb_delay = Z(:,4);
y_tau = [xi_delay; alpha_delay; up_delay; xb_delay];
M1 = [1, x_alpha, 0, 0; x_alpha, r_alpha^2, 0, 0; mi_n, -delta*mi_n, -mi_n, 0; 0, 0, 0, 1];
M2 = (pi*rho*b^2/m) * [-1, a, 0, 0; 0, -(1/8+a^2), 0, 0; 0, 0, 0, 0; 0, 0, 0, 0];
M = M1 - M2;
C1 = [0, 0, lambda_n, 0; 0, 0, -delta*lambda_n, 0; 0, 0, -lambda_n, 0; -1, -(1/2-a), 0, (c2+c4)*U/b];
C2 = [2*pi*rho*U*b^2*(c0-c1-c3)*(-1/(m*b)), (pi*rho*b^3*(U/b)+2*pi*rho*U*b^2*(c0-c1-c3)*(1/2-a))*(-1/(m*b)), 0, 2*pi*rho*U*b^2*(U/b)*(c1*c2+c3*c4)*(-1/(m*b));
2*pi*rho*U*b^3*(c0-c1-c3)*(1/2+a)*(1/(m*b^2)), (-pi*rho*b^4*(1/2-a)*(U/b)+2*pi*rho*U*b^3*(1/2+a)*(c0-c1-c3)*(1/2-a))*(1/(m*b^2)), 0, 2*pi*rho*U*b^3*(U/b)*(c1*c2+c3*c4)*(1/2+a)*(1/(m*b^2));
0, 0, 0, 0; 0, 0, 0, 0];
C = C1 - C2;
K1 = [omega_h^2, 0, gamma_n*y(3)^2, 0; 0, (omega_alpha*r_alpha)^2*(1+H_alpha*y(2)^2), -delta*gamma_n*y(3)^2, 0;
0, 0, -gamma_n*y(3)^2, 0; 0, -(U/b), 0, c2*c4*(U/b)];
K2 = [G/m, (2*pi*rho*U*b^2*(c0-c1-c3)*(U/b)*(-1/(m*b)))-(G/m*delta), -G/m, 2*pi*rho*U*b^2*(U/b)^2*(c2*c4)*(c1+c3)*(-1/(m*b));
-G/m*delta, (2*pi*rho*U*b^3*(c0-c1-c3)*(1/2+a)*(U/b)*(1/(m*b^2)))+(G/m*delta^2), G/m*delta, 2*pi*rho*U*b^3*(U/b)^2*(1/2+a)*(c2*c4)*(c1+c3)*(1/(m*b^2));
-G/m, G/m*delta, G/m, 0; 0, 0, 0, 0];
K = K1 - K2;
T = [-G/m, G/m*delta, G/m, 0; G/m*delta, -G/m*delta^2, -G/m*delta, 0; G/m, -G/m*delta, -G/m, 0; 0, 0, 0, 0];
N = zeros(4, 4);
I = eye(4);
A = [N, I; M\-K, M\-C];
B = [N; M\T];
% use a loop here ----------------------------
for k = 1:numel(xi_delay)
y_tau = [xi_delay(k); alpha_delay(k); up_delay(k); xb_delay(k)];
dydt{k} = A*y + B*y_tau;
end
dydt = dydt{:}; % return column vector only
end
I appreciate so much your help! I wasn't understand why the function didn't return the column vector.
I will avaliable the results and, if I has another questions, I will be in touch.
Again, thank you so much.
Ok. If you like my solution, don't forget to give a like. The function didn't return the column vector because the dimensions of the matrix multiplication is mismatched.
If you have a matrix of size 8x4 and another matrix of 4x1 then matrix multiplication is possible. But in your case, its not. It is 8x4 and 32x1.
Comments have been disabled.
