Variable type of a vector changes from double to complex double
1 view (last 30 days)
Show older comments
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')
0 Comments
Answers (1)
Walter Roberson
on 25 Jul 2019
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.
See Also
Categories
Find more on Wind Power in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!