Use of dde23

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

[VBBV]

@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

[VBBV]

@João PedroThe resulting output is below

[João Pedro]

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.

[VBBV]

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.