Error using ^ Incorrect dimensions for raising a matrix to a power when trying to simulate a steady state simulation of DFIG

4 views (last 30 days)
I'm Trying to simulate this code to check the steady state operation of DFIG, but it gives this Error. Not having idea of how to resolve it. Here's the Script and error:
Error using ^
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To perform elementwise matrix powers, use '.^'.
Error in steady_state_operation (line 82)
Is(i) = sqrt(ids^2+iqs^2);
%% DFIG Parameters -> Rotor Parameters referred to the Stator Side:
f = 50; % Stator Frequency (Hz)
Ps = 2E6; % Rated Stator Power (W)
N = 1500; % Rated Rotation Speed (rpm)
Vs = 690; % Rated Stator Voltage (V)
Is = 1760; % Rated Stator Current (A)
Tem = 12732; % Rated Torque (A)
p = 2; % Pole pair
u = 1/3; % Stator/Rotor Turns Ratio
Vr = 2070; % Rated Rotor Voltage (non reached)(V)
Smax = 1/3; % Maximum Slip
Vr_Stator = (Vr*Smax)*u; % Rated rotor voltage referred to stator (V)
Rs = 2.6E-3; % Stator Resistance (Ohm)
Lsi = 0.087E-3; % Leakage inductance (Stator & Rotor) (H)
Lm = 2.5E-3; % Magnetising inductance (H)
Rr = 2.9E-3; % Rotor Resistance referred to stator (Ohm)
Ls = Lm + Lsi; % Stator Inductance (H)
Lr = Lm + Lsi; % Rotor Inductance (H)
Vbus = Vr_Stator*sqrt(2); % DC Bus Voltage referred to the Stator (V)
Sigma = 1 - Lm^2/(Ls*Lr);
Vs = 690*sqrt(2/3);
ws = 2*pi*f;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum and Minimum Speed @ 2.4 MW:
wt_nom = 18;
wt_min = 9;
% Gear Box Ratio:
GB = 100;
% Torque and Speed arrays (INPUTS)
wt = [ 0.9425, 0.9425, 0.9425, 0.9524, 1.0476, 1.1429, 1.2381,...
1.3333, 1.4286, 1.5238, 1.6190, 1.7143, 1.8095, 1.8850, 1.8850,...
1.8850, 1.8850, 1.8850, 1.8850, 1.8850, 1.8850];
Torque = -1E6*[ 0.0501, 0.1164, 0.1925, 0.2696, 0.3262, 0.3883, 0.4557,...
0.5285, 0.6067, 0.6902, 0.7792, 0.8736, 0.9733, 1.0894, 1.2474,...
1.2745, 1.2745, 1.2745, 1.2745, 1.2745, 1.2745];
Temm = Torque/GB; % Convert to Machine Side
wmm = wt*(60/2/pi)*GB; % Convert to Machine Side (in rpm)
sss = (1500-wmm)/1500 % Slip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Qs = 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ref_Qs_ = 0*ones(size(sss));
i = 0;
for ss =sss,
i =i+1;
Tem = Temm(i);
Ref_Qs = Ref_Qs_(i);
Wm = ws*(1-ss);
s(i) = ss;
A = ws*ws*Lm*Lm;
B = 4*Rs*Tem*ws/3-((Vs)^2);
C = 4*Rs*Rs/(9*Lm*Lm)*((Ref_Qs^2)/(ws^2)+(Tem.^2)/(p^2));
X1 = (-B+sqrt(B.^2-4*A.*C))./(2*A);
X2 = (-B-sqrt(B.^2-4*A.*C))./(2*A);
Fs = Lm*sqrt(X1);
Fs(i) = Fs;
% Stator Currents:
ids = 2*Ref_Qs./(3*ws*Fs);
iqs = 2*Tem./(3*p*Fs);
Is(i) = sqrt(ids^2+iqs^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%Error Line %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Staor Voltages:
Uds = Rs*ids;
Uqs = Rs*iqs+ws*Fs;
Us(i) = sqrt(Uds^2+Uqs^2);
% Rotor Currents:
idr = -Ls*ids/Lm+Fs/Lm;
iqr = -Ls*iqs/Lm;
Ir(i) = sqrt(idr^2+iqr^2);
% wr:
wr = ws*ss;
% Rotor Voltage:
Udr = Rr*idr-Lr*wr.*iqr-Lm*wr.*iqs;
Uqr = Rr*iqr+Lr*wr.*idr+Lm*wr.*ids;
Vr(i) = sqrt(Udr^2+Uqr^2);
% Rotor Power:
Pr(i) = 1.5*(Udr.*idr+Uqr.*iqr);
Qr(i) = 1.5*(Uqr.*idr-Uqr.*iqr);
Mod_Sr = sqrt(Pr.^2+Qr.^2);
% Stator Flux:
Fsd = Ls*ids+Lm*idr;
Fsq = Ls*iqs+Lm*iqr;
Fs(i) = sqrt(Fsd*Fsd+Fsq*Fsq);
% Rotor Flux:
Frd = Lr*idr+Lm*ids;
Frq = Lr*iqr+Lm*iqs;
Fr(i) = sqrt(Frd*Frd+Frq*Frq);
% Mechanical Power:
Pmec(i) = Tem*Wm/p;
if Pmec(i) >= 0
R(i) = Pmec(i)/(Ps(i)+Pr(i)); % Motor Mode
else
R(i) = (Ps(i)+Pr(i))/Pmec(i); % Generator Mode
end
end
Figure(1)
Subplot(3,3,1)
hold on
plot(wmm,Temm,'r',wmm,Temm,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Tem(Nm)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,2)
hold on
plot(wmm,Ps+Pr,'r',wmm,Ps+Pr,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Pt(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,3)
hold on
plot(wmm,Pr,'r',wmm,Pr,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Pr(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,4)
hold on
plot(wmm,Ps,'r',wmm,Ps,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Ps & Pr(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,5)
hold on
plot(wmm,Is,'r',wmm,Is,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Is(A)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,6)
hold on
plot(wmm,Ir,'r',wmm,Ir,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Ir(A)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,7)
hold on
plot(wmm,Vr,'r',wmm,Vr,'or',wmm,Us,'r',wmm,Us,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Vr & Vs(V)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,8)
hold on
plot(wmm,Qs,'r',wmm,Qs,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Qs(VAr)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,9)
hold on
plot(wmm,Qr,'r',wmm,Qr,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Qr(VAr)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,10)
hold on
plot(wmm,R,'r',wmm,R,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Efficiency(pu)','FontSize',14);
xlim([700 , 2000]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Idr = 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = 0;
for ss =sss,
i =i+1;
Tem = Temm(i);
Wm = ws*(1-ss);
s(i) = ss;
A = ws*ws*Lm*Lm + (Rs*Lm/Ls)^2;
B = 4*Rs*Tem*ws/3/p-((Vs)^2);
C = 4*Rs*Rs*Tem*Tem/(9*p*p*Lm*Lm);
X1 = (-B+sqrt(B.^2-4*A.*C))./(2*A);
X2 = (-B-sqrt(B.^2-4*A.*C))./(2*A);
ims = sqrt(X1);
ims(i)=ims;
Fs = ims*Lm;
% Rotor Currents:
Idr = 0;
Iqr = -2*Tem*Ls./(3*p*Lm*Fs);
Ir(i) = sqrt(Idr^2+Iqr^2);
% Stator Currents:
Ids = Fs/Ls;
Iqs = -Iqr*Lm/Ls;
Is(i) = sqrt(Ids^2+Iqs^2);
% Stator Voltages:
Uds = Rs*Ids;
Uqs = Rs*Iqs+ws*Fs;
Us(i) = sqrt(Uds^2+Uqs^2);
% Stator Power:
Ps(i) = 1.5*(Rs*Ids.^2+Rs*Iqs.^2+ws*Fs.*Iqs);
Qs(i) = 1.5*(Uqs.*Ids-Uds.*Iqs);
Mod_Ss = sqrt(Ps.^2+Qs.^2);
% wr:
wr = ws*ss;
% Rotor Voltage:
Udr = Rr*Idr-Lr*wr.*Iqr-Lm*wr.*Iqs;
Uqr = Rr*Iqr+Lr*wr.*Idr+Lm*wr.*Ids;
Vr(i) = sqrt(Udr^2+Uqr^2);
% Rotor Power:
Pr(i) = 1.5*(Udr.*Idr+Uqr.*Iqr);
Qr(i) = 1.5*(Uqr.*Idr-Udr.*Iqr);
Mod_Sr = sqrt(Pr.^2+Qr.^2);
% Stator Flux:
Fsd = Ls*Ids+Lm*Idr;
Fsq = Ls*Iqs+Lm*Iqr;
Fs(i) = sqrt(Fsd*Fsd+Fsq*Fsq);
% Rotor Flux:
Frd = Lr*Idr+Lm*Ids;
Frq = Lr*Iqr+Lm*Iqs;
Fr(i) = sqrt(Frd*Frd+Frq*Frq);
% Mechanical Power:
Pmec(i) = Tem*Wm/p;
if Pmec(i) >= 0
R(i) = Pmec(i)/(Ps(i)+Pr(i)); % Motor Mode
else
R(i) = (Ps(i)+Pr(i))/Pmec(i); % Generator Mode
end
end
Figure(1)
Subplot(3,3,1)
hold on
plot(wmm,Temm,'g',wmm,Temm,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Tem(Nm)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,2)
title('red: Qs = 0, Green: Idr = 0','fontsize',12)
hold on
plot(wmm,Ps+Pr,'g',wmm,Ps+Pr,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Pt(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,3)
hold on
plot(wmm,Pr,'r',wmm,Pr,'or','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Pr(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,4)
hold on
plot(wmm,Ps,'g',wmm,Ps,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Ps & Pr(W)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,5)
hold on
plot(wmm,Is,'g',wmm,Is,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Is(A)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,6)
hold on
plot(wmm,Ir,'g',wmm,Ir,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Ir(A)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,7)
hold on
plot(wmm,Vr,'g',wmm,Vr,'og',wmm,Us,'g',wmm,Us,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Vr & Vs(V)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,8)
hold on
plot(wmm,Qs,'g',wmm,Qs,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Qs(VAr)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,9)
hold on
plot(wmm,Qr,'g',wmm,Qr,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Qr(VAr)','FontSize',14);
xlim([700 , 2000]);
Subplot(3,3,10)
hold on
plot(wmm,R,'g',wmm,R,'og','linewidth',1.5);
xlabel('n(rpm)','FontSize',14);
ylabel('Efficiency(pu)','FontSize',14);
xlim([700 , 2000]);

Answers (1)

Star Strider
Star Strider on 23 Mar 2019
I’m not certain what you want, or what your code does.
One option is to address only the latest elements of the vectors (even though you don’t specifically subscript them earlier in your code when you create them):
Is(i) = sqrt(ids(i)^2+iqs(i)^2);
If that provides the result you want, make similar changes to the other lines that will throw the same error.
If it doesn’t, I have no other suggestions.

Community Treasure Hunt

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

Start Hunting!