Variable type of a vector changes from double to complex double

Hello everyone. I am trying to solve an ode using the forward Euler method. I have a strange problem! When I try to save the solution of ode in a vector through a loop, the type of this vector changes from double to complex double. I do not have any complex numbers in my code. My code is as follows:

clear all
close all
clc
tf=5;
dt=0.001;
t=0:dt:tf;
N=length(t);
% System Parameters %
Rs=2.9; Rr=1.52; Ls=0.223; Lr=0.229; M=0.217; P=2; Kt=8.29e-5; roh=1.225; R=2.5;
Jt=0.0048; ng=4.86; C_pmax=0.45; landa_opt=7; wn=1050*pi/30; ws=1500*pi/30; Vs=220; phi_s=Vs/ws; s=(ws-wn)/ws;
sigma=1-(M^2/(Ls*Lr)); K_opt=(0.5*pi*roh*(R^5)*C_pmax)/(landa_opt^3);
% Controller Parameters %
k1_d=1; k2_d=1; epsilon_d=1; h1_d=1; h2_d=1; delta_d=10;ld=5/3;
k1_q=1; k2_q=1; epsilon_q=1; h1_q=1; h2_q=1; delta_q=10;lq=5/3;
%========== Initializing ==========%
% Wind Speed %
v=0*ones(N,1);
% Actual Signals %
Ird=zeros(N+1,1);
Irq=zeros(N+1,1); 
% Refernce Signals %
Ird_ref=zeros(N,1); Irdref_dot=zeros(N,1);
Irq_ref=zeros(N,1); Irqref_dot=zeros(N,1);
Tem_ref=zeros(N,1);
wmr_ref=zeros(N,1);
% Tracking Errors %
e_Ird=zeros(N,1);e_Irq=zeros(N,1);
%
int_rd_dot=zeros(N,1); int_rd=zeros(N+1,1);
int_rq_dot=zeros(N,1); int_rq=zeros(N+1,1);
% Sliding Surface %
s_Ird=zeros(N,1);s_Irq=zeros(N,1);
% Control Inputs %
vrd=zeros(N,1);vrq=zeros(N,1);
% Torques %
Tem=zeros(N,1); % Electro Magnetic Torque %
Tg=zeros(N,1); % Generator Torque %
Ta=zeros(N,1); % Aerodynamic Torque %
% Tip Speed Ratio %
landa=zeros(N,1);
Gama=zeros(N,1);
% Pitch Angle %
Betta=zeros(N,1);
% Power Coefficient %
Cp=zeros(N,1);
% Mechanical Angualar speed (Turbine Rotor Speed)%
wmr=zeros(N+1,1);
wmr(1,1)=100;
% Main %
for k=1:N
    
  % Wind Speed %
  v(k,1)=12;
% Refernce Signals %
Ird_ref(k,1)=(Vs/(ws*M)); Irdref_dot(k,1)=0;
Irq_ref(k,1)=(1/ng)*(-Ls/(P*M*phi_s))*(0.5*pi*roh*(R^3)*(C_pmax/landa_opt)*((v(k,1))^2)); Irqref_dot(k,1)=0;
Tem_ref(k,1)=(-P*M*phi_s/(Ls))*Irq_ref(k,1);
wmr_ref(k,1)=landa_opt*v(k,1)/R;
% Functions in System Equations %
fd=(1/(sigma*Lr))*(-Rr*Ird(k,1)+s*ws*sigma*Lr*Irq(k,1)); gd=(1/(sigma*Lr));
fq=(1/(sigma*Lr))*(-Rr*Irq(k,1)-s*ws*sigma*Lr*Ird(k,1)-s*ws*(M/Ls)*phi_s); gq=gd;

% Tracking Errors %
e_Ird(k,1)=Ird(k,1)-Ird_ref(k,1); 
e_Irq(k,1)=Irq(k,1)-Irq_ref(k,1); 

% Sliding Surfaces %
int_rd_dot(k,1)=k1_d*((e_Ird(k,1))^ld)+k2_d*(e_Ird(k,1))+epsilon_d*(tanh(e_Ird(k,1)));
int_rd(k+1,1)=dt*int_rd_dot(k,1)+int_rd(k,1);
s_Ird(k,1)=e_Ird(k,1)+int_rd(k,1);
%
int_rq_dot(k,1)=k1_q*((e_Irq(k,1))^lq)+k2_q*(e_Irq(k,1))+epsilon_q*(tanh(e_Irq(k,1)));
int_rq(k+1,1)=dt*int_rq_dot(k,1)+int_rq(k,1);
s_Irq(k,1)=e_Irq(k,1)+int_rq(k,1);

% Control Inputs %
vrd(k,1)=(-1/gd)*(fd-Irdref_dot(k,1)+k1_d*((e_Ird(k,1))^ld)+k2_d*(e_Ird(k,1))+epsilon_d*(tanh(e_Ird(k,1)))+h1_d*((s_Ird(k,1))^ld)+h2_d*(s_Ird(k,1))+delta_d*(tanh(s_Ird(k,1))));
%
vrq(k,1)=(-1/gq)*(fq-Irqref_dot(k,1)+k1_q*((e_Irq(k,1))^lq)+k2_q*(e_Irq(k,1))+epsilon_q*(tanh(e_Irq(k,1)))+h1_q*((s_Irq(k,1))^lq)+h2_q*(s_Irq(k,1))+delta_q*(tanh(s_Irq(k,1))));

System Equations %
Ird_dot=fd+gd*vrd(k,1);
Irq_dot=fq+gq*vrq(k,1);
Ird(k+1,1)=dt*Ird_dot+Ird(k,1);
Irq(k+1,1)=dt*Irq_dot+Irq(k,1);
Tem(k,1)=-P*(M/Ls)*phi_s*(Irq(k,1)); Tg(k,1)=ng*Tem(k,1);
landa(k,1)=R*(wmr(k,1))/(v(k,1)); Betta(k,1)=0;
Gama(k,1)=(1/(landa(k,1)+0.08*Betta(k,1)))-(0.035/((Betta(k,1))^3+1));
Cp(k,1)=(0.5176*(116*Gama(k,1)-0.4*Betta(k,1)-5)*(exp(-21*Gama(k,1)))+0.0068*landa(k,1));
Ta(k,1)=0.5*pi*roh*(R^3)*(Cp(k,1)/landa(k,1))*((v(k,1))^2);
wmr_dot=(1/Jt)*(Ta(k,1)-Kt*wmr(k,1)-Tg(k,1));
wmr(k+1,1)=dt*wmr_dot+wmr(k,1);
%
end
figure(1)
plot(t,Ird_ref,'b',t,Ird(1:end-1,1),'--r','linewidth',2)
xlabel('time(s)')
legend('Ird_r_e_f','Ird')
figure(2)
plot(t,Irq_ref,'b',t,Irq(1:end-1,1),'--r','linewidth',2)
xlabel('time(s)')
legend('Irq_r_e_f','Irq')
% figure(3)
% plot(t,wmr_ref,'b',t,wmr(1:end-1,1),'--r','linewidth',2)
% xlabel('time(s)')
% legend('wmr_r_e_f','wmr')
% figure(4)
% plot(t,Tem_ref,'b',t,Tem,'--r','linewidth',2)
% xlabel('time(s)')
% legend('Tem_r_e_f','Tem')
% figure(5)
% plot(t,Ta,'b',t,Tg,'--r','linewidth',2)
% xlabel('time(s)')
% legend('T_a','T_g')

Answer 1

[Walter Roberson]

You have

ld=5/3;              %fractional
[...]
Ird=zeros(N+1,1);   %so initialized to 0
[...]
Ird_ref(k,1)=(Vs/(ws*M));  %0 or positive
[...]
e_Ird(k,1)=Ird(k,1)-Ird_ref(k,1);   %0 minus 0 or positive, so 0 or negative
[...]
int_rd_dot(k,1)=k1_d*((e_Ird(k,1))^ld)+k2_d*(e_Ird(k,1))+epsilon_d*(tanh(e_Ird(k,1)));

We established that e_Ird is 0 or negative, and it is being raised to the fractional value ld (=5/3) . That is going to give you a result that is either 0 or complex valued.

[hossein dastres]

Thank you so much.