Variable type of a vector changes from double to complex double

1 次查看(过去 30 天)
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')

回答(1 个)

Walter Roberson
Walter Roberson 2019-7-25
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.

类别

Help CenterFile Exchange 中查找有关 Wind Power 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by