MATLAB Answers

Summer
0

"Subscript indices must either be real positive integers or logicals." again! See codes

Asked by Summer
on 10 Jan 2015
Latest activity Answered by Deepak Bhattarai on 9 May 2019
Hi,
I have created two functions in Matlab. One is supposed to produce a concentration profile for a given location and time, but the current code gives this error "Subscript indices must either be real positive integers or logicals.". The second function, called testing, recalls the first function to do the calculation and uses loops instead, but still the same error is appearing.
1. Concentration profile code
function [c,t] = Disp_Conc(x,y,z,t)
%This function calculates the concentration of a an agent at a particular time and location
% Data generated through Dispersion model (This is a 1x61 vector)
[x,cc_x,b_x,betac_x,zc_x,sig_x,t,xc_t,bx_t,betax_t] = Read_Predict('predict');
% Time averaged volume concentration: concentration contour parameters
%erf = error function
sr2 = sqrt(2.0);
xa=(x-xc_t(t)+bx_t(t))/(sr2*betax_t(t));
xb = (x-xc_t(t)-bx_t(t))/(sr2*betax_t(t));
ya = (y+b_x(x))/(sr2*betac_x(x));
yb = (y-b_x(x))/(sr2*betac_x(x));
za = (z-zc_x(x))/(sr2*sig_x(x));
zb = (z+zc_x(x))/(sr2*sig_x(x));
c = (cc_x(x) * (erf(xa)-erf(xb)) * (erf(ya)-erf(yb)) * (exp(-za*za)+exp(-zb*zb)));
time=tm(t);
end
**Running this function, I get an error saying "Subscript indices must either be real positive integers or logicals." due to the fact that values of t and x are both fractional (both x and t are read from a 1x61 vector). To overcome this, I created another .m file that recalls the concentration function and uses loops as below:
2. Function "testing":
function [c, t] = testing( x,y,z,t )
%The values of x and t are already being read from a saved 1x61 vector. z and y can vary and must be specified by the user for a given location. The value of z should be set equal to zero for all cases and doesn't change.
for ts=1:length(t)
for xs=1:length(x)
z=0;
for y=0:100; %(k index?)
X=x(ts);
T=t(xs);
[c(ts),time(ts)]=Disp_Conc(X,y,z,T);
end
end
end
end
**Running "testing" for a particular parameters say testing(5,0,0,8.79), I get the below errors:
"Subscript indices must either be real positive integers or logicals.
Error in Disp_Conc (line 12)
xa=(x-xc_t(t)+bx_t(t))/(sr2*betax_t(t));
Error in testing (line 10)
[c(ts),time(ts)]=Disp_Conc(X,y,z,T);"
I don't know what is wrong this time? The original problem was because of indexing and now I still get the same error with the new function.
Any help would be greatly appreciated. Thx.

  0 Comments

Sign in to comment.

4 Answers

Answer by Star Strider
on 10 Jan 2015
 Accepted Answer

In the line that throws the error:
xa=(x-xc_t(t)+bx_t(t))/(sr2*betax_t(t));
what is the value of ‘t’?

  10 Comments

Typically there are two methods for looping over non-integer values in MATLAB:
  1. derive the values within each loop, using a function or calculating the values from the loop variables.
  2. use the loop variables to look up the values in a predefined array.
Simple examples for each of these, where we want the non-integer value 1/k in each iteration:
for k = 1:4
disp(1/k)
end
A = 1./(1:4);
for k = 1:4
disp(A(k))
end
While the second method looks more complicated, it actually has many advantages in MATLAB: it helps memory management , and makes it easier to write vectorized code. Note how the second example uses indexing, and the index value is an integer!
As you have arrays of data, you should consider using the look-up method. You will probably need these function to help you: size, numel.
"xa(k1) = (x-xc_t(k1)+bx_t(k1))/(sr2*betax_t(k1));"
Worked this time! Thanks :)

Sign in to comment.


Answer by Stephen Cobeldick on 10 Jan 2015
Edited by Stephen Cobeldick on 10 Jan 2015

In MATLAB all array indices must be logical or positive numeric integers. This means that the following is permitted:
>> A = [123,456,789];
>> A(2) % 2 is an integer
ans = 456
But the following produces an error:
>> A(3.14159) % 3.14159 is not an integer
Subscript indices must either be real positive integers or logicals.
Note that negative integers and zero are not permitted indices.
You need to review your code and check all of the array indexing. In particular it seems that you are using some data values (e.g. t, x) as indices.

  4 Comments

Show 1 older comment
Stephen, very good explanation (+ a vote). Since we answer this at least once or twice a day, I don't know why it's not on the FAQ yet. So I used your excellent explanation as the basis for a new FAQ entry: http://matlab.wikia.com/wiki/FAQ#How_do_I_fix_the_error_.22Subscript_indices_must_either_be_real_positive_integers_or_logicals..22.3F
By the way, empty seems to be okay
>> emptyvector = []
emptyvector =
[]
>> A(emptyvector)
ans =
[]
True, [] simply returns an empty array. I have edited my answer to reflect this.

Sign in to comment.


Answer by Luca
on 21 Aug 2017
Edited by Walter Roberson
on 21 Aug 2017

clc; clear;
teta=linspace(0,pi/2,100);
R=1;
g=9.81;
v_quadro = g*R*cos(teta.*2) + R*g(2^0.5-1)-4*R^3*(1-cos(teta))*cos(teta)/((2*R-(2^0.5)*R)^2);
Where is the error, please? Thanks.

  2 Comments

Your g is not a vector, but you are trying to index it at (2^0.5-1).
MATLAB does not have any implicit multiplications. You need to add ".*" or "*" in g(2^0.5-1)

Sign in to comment.


Answer by Deepak Bhattarai on 9 May 2019

Please anyone can help me to fix the following matlab error, it says:
Subscript indices must either be real positive integers or logicals.
Error in withoutSVCbhu>NetPowers (line 396)
PNET(loadbus(ii)) = PNET(loadbus(ii)) - PLOAD(ii);
Error in withoutSVCbhu>NewtonRaphson (line 362)
[PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,QGEN,...
Error in withoutSVCbhu (line 270)
[VM,VA,it] = NewtonRaphson(nmax,tol,itmax,ngn,nld,nbb,bustype,...
Following is the complete matlab Program;
%Bus data
%nbb = number of buses
%bustype = type of bus
%VM = nodal voltage magnitude
%VA = nodal voltage phase angle
nbb = 50 ;
bustype(1) = 3 ; VM(1) = 1 ; VA(1) =0;
bustype(2) = 3 ; VM(2) = 1 ; VA(2) =0 ;
bustype(3) = 3 ; VM(3) = 1 ; VA(3) =0 ;
bustype(4) = 3 ; VM(4) = 1 ; VA(4) =0 ;
bustype(5) = 3 ; VM(5) = 1 ; VA(5) =0 ;
bustype(6) = 3 ; VM(6) = 1 ; VA(6) =0 ;
bustype(7) = 3 ; VM(7) = 1 ; VA(7) =0 ;
bustype(8) = 3 ; VM(8) = 1 ; VA(8) =0 ;
bustype(9) = 3 ; VM(9) = 1 ; VA(9) =0 ;
bustype(10) = 3 ; VM(10) = 1 ; VA(10) =0 ;
bustype(11) = 3 ; VM(11) = 1 ; VA(11) =0 ;
bustype(12) = 3 ; VM(12) = 1 ; VA(12) =0 ;
bustype(13) = 3 ; VM(13) = 1 ; VA(13) =0 ;
bustype(14) = 3 ; VM(14) = 1 ; VA(14) =0 ;
bustype(15) = 3 ; VM(15) = 1 ; VA(15) =0 ;
bustype(16) = 3 ; VM(16) = 1 ; VA(16) =0 ;
bustype(17) = 3 ; VM(17) = 1 ; VA(17) =0 ;
bustype(18) = 3 ; VM(18) = 1 ; VA(18) =0 ;
bustype(19) = 3 ; VM(19) = 1 ; VA(19) =0 ;
bustype(20) = 3 ; VM(20) = 1 ; VA(20) =0 ;
bustype(21) = 3 ; VM(21) = 1 ; VA(21) =0 ;
bustype(22) = 3 ; VM(22) = 1 ; VA(22) =0 ;
bustype(23) = 3 ; VM(23) = 1 ; VA(23) =0 ;
bustype(24) = 3 ; VM(24) = 1 ; VA(24) =0 ;
bustype(25) = 3 ; VM(25) = 1 ; VA(25) =0 ;
bustype(26) = 3 ; VM(26) = 1 ; VA(26) =0 ;
bustype(27) = 3 ; VM(27) = 1 ; VA(27) =0 ;
bustype(28) = 3 ; VM(28) = 1 ; VA(28) =0 ;
bustype(29) = 3 ; VM(29) = 1 ; VA(29) =0 ;
bustype(30) = 3 ; VM(30) = 1 ; VA(30) =0 ;
bustype(31) = 3 ; VM(31) = 1 ; VA(31) =0 ;
bustype(32) = 3 ; VM(32) = 1 ; VA(32) =0 ;
bustype(33) = 3 ; VM(33) = 1 ; VA(33) =0 ;
bustype(34) = 3 ; VM(34) = 1 ; VA(34) =0 ;
bustype(35) = 3 ; VM(35) = 1 ; VA(35) =0 ;
bustype(36) = 3 ; VM(36) = 1 ; VA(36) =0 ;
bustype(37) = 3 ; VM(37) = 1 ; VA(37) =0 ;
bustype(38) = 3 ; VM(38) = 1 ; VA(38) =0 ;
bustype(39) = 3 ; VM(39) = 1 ; VA(39) =0 ;
bustype(40) = 3 ; VM(40) = 1 ; VA(40) =0 ;
bustype(41) = 3 ; VM(41) = 1 ; VA(41) =0 ;
bustype(42) = 3 ; VM(42) = 1 ; VA(42) =0 ;
bustype(43) = 3 ; VM(43) = 1 ; VA(43) =0 ;
bustype(44) = 3 ; VM(44) = 1 ; VA(44) =0 ;
bustype(45) = 1 ; VM(45) = 0.8832 ; VA(45) =0 ;
bustype(46) = 2 ; VM(46) = 0.986 ; VA(46) =0 ;
bustype(47) = 2 ; VM(47) = 0.943 ; VA(47) =0 ;
bustype(48) = 2 ; VM(48) = 0.999 ; VA(48) =0 ;
bustype(49) = 2 ; VM(49) = 1.003 ; VA(49) =0 ;
bustype(50) = 2 ; VM(50) = 0.984 ; VA(50) =0 ;
%Generator data
%ngn = number of generators
%genbus = generator bus number
%PGEN = scheduled active power contributed by the generator
%QGEN = scheduled reactive power contributed by the generator
%QMAX = generator reactive power upper limit
%QMIN = generator reactive power lower limit
ngn = 6 ;
genbus(1) = 45 ; PGEN(1) = 0 ; QGEN(1) = 0 ; QMAX(1) = 5.386 ; QMIN(1) = -2.692;
genbus(2) = 46 ; PGEN(2) =3.409; QGEN(2) =0 ; QMAX(2)=1.774; QMIN(2) = -0.887;
genbus(3) = 47 ; PGEN(3) =0.626; QGEN(3) =0 ; QMAX(3) =0.317; QMIN(3) = -0.158;
genbus(4) = 48 ; PGEN(4) = 0.882; QGEN(4) = 0 ; QMAX(4) =0.605; QMIN(4) = -0.042;
genbus(5) = 49 ; PGEN(5) =0.244; QGEN(5) =0 ; QMAX(5) =0.113; QMIN(5) = -0.0572;
genbus(6) = 50 ; PGEN(6) =0.409 ; QGEN(6) = 0; QMAX(6) =0.192 ; QMIN(6) = -0.096;
%
%Transmission line data
%ntl = number of transmission lines
%tlsend = sending end of transmission line
%tlrec = receiving end of transmission line
%tlresis = series resistance of transmission line
%tlreac = series reactance of transmission line
%tlcond = shunt conductance of transmission line
%tlsuscep = shunt susceptance of transmission line
ntl = 47 ;
tlsend(1) = 2 ; tlrec(1) = 13 ; tlresis(1) = 0.0078 ; tlreac(1) = 0.0444 ;
tlcond(1) = 0 ; tlsuscep(1) = 0.0761 ;
tlsend(2) = 2 ; tlrec(2) = 10 ; tlresis(2) = 0.0043 ; tlreac(2) = 0.0245 ;
tlcond(2) = 0 ; tlsuscep(2) = 0.0421 ;
tlsend(3) = 2 ; tlrec(3) = 9 ; tlresis(3) = 0.0101 ; tlreac(3) = 0.0577 ;
tlcond(3) = 0 ; tlsuscep(3) = 0.0989 ;
tlsend(4) = 2 ; tlrec(4) = 9 ; tlresis(4) = 0.01011 ; tlreac(4) = 0.0577 ;
tlcond(4) = 0 ; tlsuscep(4) = 0.0989 ;
tlsend(5) = 13 ; tlrec(5) = 6 ; tlresis(5) = 0.0065 ; tlreac(5) = 0.037 ;
tlcond(5) = 0 ; tlsuscep(5) = 0.0633 ;
tlsend(6) = 6 ; tlrec(6) = 14 ; tlresis(6) = 0.0067 ; tlreac(6) = 0.0383 ;
tlcond(6) = 0 ; tlsuscep(6) = 0.0657 ;
tlsend(7) = 14 ; tlrec(7) = 4 ; tlresis(7) = 0.0029 ; tlreac(7) = 0.0166 ;
tlcond(7) = 0 ; tlsuscep(7) = 0.0285 ;
tlsend(8) = 14 ; tlrec(8) = 15 ; tlresis(8) = 0.0044 ; tlreac(8) = 0.1921 ;
tlcond(8) = 0 ; tlsuscep(8) = 0.0430 ;
tlsend(9) = 4 ; tlrec(9) = 15 ; tlresis(9) = 0.0079 ; tlreac(9) = 0.0453 ;
tlcond(9) = 0 ; tlsuscep(9) = 0.0777 ;
tlsend(10) = 10 ; tlrec(10) = 11 ; tlresis(10) = 0.0002 ; tlreac(10) = 0.00099 ;
tlcond(10) = 0 ; tlsuscep(10) = 0.0017 ;
tlsend(11) = 10 ; tlrec(11) = 12 ; tlresis(11) = 0.0058 ; tlreac(11) = 0.033 ;
tlcond(11) = 0 ; tlsuscep(11) = 0.0566 ;
tlsend(12) = 10 ; tlrec(12) = 9 ; tlresis(12) = 0.0058 ; tlreac(12) = 0.0336 ;
tlcond(12) = 0 ; tlsuscep(12) = 0.0575 ;
tlsend(13) = 11 ; tlrec(13) = 12 ; tlresis(13) = 0.0059 ; tlreac(13) = 0.0338 ;
tlcond(13) = 0 ; tlsuscep(13) = 0.0580 ;
tlsend(14) = 3 ; tlrec(14) = 17 ; tlresis(14) = 0.0094 ; tlreac(14) = 0.0223 ;
tlcond(14) = 0 ; tlsuscep(14) = 0.0051 ;
tlsend(15) = 3 ; tlrec(15) = 19 ; tlresis(15) = 0.0289 ; tlreac(15) = 0.0687 ;
tlcond(15) = 0 ; tlsuscep(15) = 0.0158 ;
tlsend(16) = 17 ; tlrec(16) = 16 ; tlresis(16) = 0.0276 ; tlreac(16) = 0.0658 ;
tlcond(16) = 0 ; tlsuscep(16) = 0.0151 ;
tlsend(17) = 19 ; tlrec(17) = 25 ; tlresis(17) = 0.0217 ; tlreac(17) = 0.0517 ;
tlcond(17) = 0 ; tlsuscep(17) = 0.0119 ;
tlsend(18) = 19 ; tlrec(18) = 18 ; tlresis(18) = 0.03172 ; tlreac(18) = 0.0755 ;
tlcond(18) = 0 ; tlsuscep(18) = 0.0174 ;
tlsend(19) = 25 ; tlrec(19) = 26 ; tlresis(19) = 0.0098 ; tlreac(19) = 0.0233 ;
tlcond(19) = 0 ; tlsuscep(19) = 0.0054 ;
tlsend(20) = 26 ; tlrec(20) = 27 ; tlresis(20) = 0.0432 ; tlreac(20) = 0.1028 ;
tlcond(20) = 0 ; tlsuscep(20) = 0.0236 ;
tlsend(21) = 18 ; tlrec(21) = 20 ; tlresis(21) = 0.0776 ; tlreac(21) = 0.1847 ;
tlcond(21) = 0 ; tlsuscep(21) = 0.0425 ;
tlsend(22) = 20 ; tlrec(22) = 21 ; tlresis(22) = 0.0304 ; tlreac(22) = 0.0723 ;
tlcond(22) = 0 ; tlsuscep(22) = 0.0166 ;
tlsend(23) = 20 ; tlrec(23) = 22 ; tlresis(23) = 0.0426 ; tlreac(23) = 0.1014 ;
tlcond(23) = 0 ; tlsuscep(23) = 0.0233 ;
tlsend(24) = 22 ; tlrec(24) = 23 ; tlresis(24) = 0.01433 ; tlreac(24) = 0.0341 ;
tlcond(24) = 0 ; tlsuscep(24) = 0.0079 ;
tlsend(25) = 23 ; tlrec(25) = 24 ; tlresis(25) = 0.0466 ; tlreac(25) = 0.1108 ;
tlcond(25) = 0 ; tlsuscep(25) = 0.0255 ;
tlsend(26) = 1 ; tlrec(26) = 7 ; tlresis(26) = 0.00043 ; tlreac(26) = 0.0046 ;
tlcond(26) = 0 ; tlsuscep(26) = 0.1445 ;
tlsend(27) = 1 ; tlrec(27) = 8 ; tlresis(27) = 0.0026 ; tlreac(27) = 0.0280 ;
tlcond(27) = 0 ; tlsuscep(27) = 0.8789 ;
tlsend(28) = 1 ; tlrec(28) = 8 ; tlresis(28) = 0.0026 ; tlreac(28) = 0.0280 ;
tlcond(28) = 0 ; tlsuscep(28) = 0.8789 ;
tlsend(29) = 1 ; tlrec(29) = 8 ; tlresis(29) = 0.0027 ; tlreac(29) = 0.0287 ;
tlcond(29) = 0 ; tlsuscep(29) = 0.8999 ;
tlsend(30) = 7 ; tlrec(30) = 8 ; tlresis(30) = 0.00225 ; tlreac(30) = 0.0240 ;
tlcond(30) = 0 ; tlsuscep(30) = 0.7525 ;
tlsend(31) = 39 ; tlrec(31) = 38 ; tlresis(31) = 0.0815 ; tlreac(31) = 0.1939 ;
tlcond(31) = 0 ; tlsuscep(31) = 0.0028 ;
tlsend(32) = 39 ; tlrec(32) = 40 ; tlresis(32) = 0.0759 ; tlreac(32) = 0.1807 ;
tlcond(32) = 0 ; tlsuscep(32) = 0.0026 ;
tlsend(33) = 38 ; tlrec(33) = 35 ; tlresis(33) = 0.05623 ; tlreac(33) = 0.13379 ;
tlcond(33) = 0 ; tlsuscep(33) = 0.0019 ;
tlsend(34) = 35 ; tlrec(34) = 37 ; tlresis(34) = 0.1244 ; tlreac(34) = 0.2959 ;
tlcond(34) = 0 ; tlsuscep(34) = 0.0043 ;
tlsend(35) = 35 ; tlrec(35) = 36 ; tlresis(35) = 0.0894 ; tlreac(35) = 0.2127 ;
tlcond(35) = 0 ; tlsuscep(35) = 0.0030 ;
tlsend(36) = 35 ; tlrec(36) = 34 ; tlresis(36) = 0.0439 ; tlreac(36) = 0.1044 ;
tlcond(36) = 0 ; tlsuscep(36) = 0.0015 ;
tlsend(37) = 34 ; tlrec(37) = 33 ; tlresis(37) = 0.0657 ; tlreac(37) = 0.1565 ;
tlcond(37) = 0 ; tlsuscep(37) = 0.0023 ;
tlsend(38) = 33 ; tlrec(38) = 31 ; tlresis(38) = 0.0064 ; tlreac(38) = 0.0150 ;
tlcond(38) = 0 ; tlsuscep(38) = 0.00022 ;
tlsend(39) = 31 ; tlrec(39) = 32 ; tlresis(39) = 0.0426 ; tlreac(39) = 0.1015 ;
tlcond(39) = 0 ; tlsuscep(39) = 0.0015 ;
tlsend(40) = 31 ; tlrec(40) = 30 ; tlresis(40) = 0.0906 ; tlreac(40) = 0.2156 ;
tlcond(40) = 0 ; tlsuscep(40) =0.0031 ;
tlsend(41) = 30 ; tlrec(41) = 29 ; tlresis(41) = 0.0812 ; tlreac(41) = 0.1931 ;
tlcond(41) = 0 ; tlsuscep(41) = 0.0028 ;
tlsend(42) = 29 ; tlrec(42) = 28 ; tlresis(42) = 0.0056 ; tlreac(42) = 0.0133 ;
tlcond(42) = 0 ; tlsuscep(42) = 0.00019 ;
tlsend(43) = 28 ; tlrec(43) = 5 ; tlresis(43) = 0.0115 ; tlreac(43) = 0.0273 ;
tlcond(43) = 0 ; tlsuscep(43) = 0.0004 ;
tlsend(44) = 40 ; tlrec(44) = 41 ; tlresis(44) = 0.0622 ; tlreac(44) = 0.1479 ;
tlcond(44) = 0 ; tlsuscep(44) = 0.0021 ;
tlsend(45) = 41 ; tlrec(45) = 42 ; tlresis(45) = 0.03344 ; tlreac(45) = 0.0795 ;
tlcond(45) = 0 ; tlsuscep(45) = 0.0011 ;
tlsend(46) = 41 ; tlrec(46) = 43 ; tlresis(46) = 0.1002 ; tlreac(46)= 0.2383 ;
tlcond(46) = 0 ; tlsuscep(46) = 0.0034 ;
tlsend(47) = 43 ; tlrec(47) = 44 ; tlresis(47) = 0.0551 ; tlreac(47) = 0.1311 ;
tlcond(47) = 0 ; tlsuscep(47) = 0.0019 ;
%
%Shunt data
%nsh = number of shunt elements
%shbus = shunt element bus number
%shresis = resistance of shunt element
%shreac = reactance of shunt element:
%+ve for inductive reactance and –ve for capacitive reactance
nsh = 0 ;
shbus(1) = 0 ; shresis(1) = 0 ; shreac(1) = 0 ;
%
%Load data
%nld = number of load elements
%loadbus = load element bus number
%PLOAD = scheduled active power consumed at the bus
%QLOAD = scheduled reactive power consumed at the bus
nld = 44 ;
loadbus(1) = 1 ; PLOAD(1) = 0 ; QLOAD(1) = 0 ;
loadbus(2) = 2 ; PLOAD(2) = 0.0113 ; QLOAD(2) = 0 ;
loadbus(3) = 3 ; PLOAD(3) = 0.0031 ; QLOAD(3) = 0;
loadbus(4) = 4 ; PLOAD(4) = 0.003 ; QLOAD(4) = 0 ;
loadbus(5) = 5 ; PLOAD(5) = 0.0017 ; QLOAD(5) = 0 ;
loadbus(6) = 6 ; PLOAD(6) = 0.0117 ; QLOAD(6) = 0 ;
loadbus(7) = 7 ; PLOAD(7) = 0 ; QLOAD(7) = 0 ;
loadbus(8) = 8 ; PLOAD(8) = 2.0883 ; QLOAD(8) = -0.752 ;
loadbus(9) = 9 ; PLOAD(9) = 0.643 ; QLOAD(9) = -0.1683 ;
loadbus(10) = 10; PLOAD(10) = 0 ; QLOAD(10) = 0 ;
loadbus(11) = 11 ; PLOAD(11) = 0.582 ; QLOAD(11) = 0.298 ;
loadbus(12) = 12 ; PLOAD(12) = 0 ; QLOAD(12) = 0 ;
loadbus(13) = 13 ; PLOAD(13) = 0 ; QLOAD(13) = 0 ;
loadbus(14) = 14 ; PLOAD(14) = 0.013 ; QLOAD(14) = 0 ;
loadbus(15) = 15 ; PLOAD(15) = 0 ; QLOAD(15) = 0 ;
loadbus(16) = 16 ; PLOAD(16) = 0.0248 ; QLOAD(16) = 0 ;
loadbus(17) = 17 ; PLOAD(17) = 0.0189 ; QLOAD(17) = 0 ;
loadbus(18) = 18 ; PLOAD(18) = 0.0915 ; QLOAD(18) = 0.0383 ;
loadbus(19) = 19 ; PLOAD(19) = 0.0017 ; QLOAD(19) = 0 ;
loadbus(20) = 20 ; PLOAD(20) = 0.0032 ; QLOAD(20) = 0 ;
loadbus(21) = 21 ; PLOAD(21) = 0.0355 ; QLOAD(21) = 0;
loadbus(22) = 22 ; PLOAD(22) = 0.0022 ; QLOAD(22) = 0.0012 ;
loadbus(23) = 23 ; PLOAD(23) = 0.4979 ; QLOAD(23) = -0.023 ;
loadbus(24) = 24 ; PLOAD(24) = 0 ; QLOAD(24) = 0 ;
loadbus(25) = 25 ; PLOAD(25) = 0.0212 ; QLOAD(25) = 0.0018 ;
loadbus(26) = 26 ; PLOAD(26) = 0.148 ; QLOAD(26) = 0.0208 ;
loadbus(27) = 27 ; PLOAD(27) = 0.5108 ; QLOAD(27) = 0.2078 ;
oadbus(28) = 28 ; PLOAD(28) = 0.0117 ; QLOAD(28) = 0 ;
loadbus(29) = 29 ; PLOAD(29) = 0 ; QLOAD(29) = 0 ;
loadbus(30) = 30 ; PLOAD(30) = 0.0432 ; QLOAD(30) = 0.0119 ;
loadbus(31) = 31 ; PLOAD(31) = 0.0541 ; QLOAD(31) = 0 ;
loadbus(32) = 32 ; PLOAD(32) = 0.0432 ; QLOAD(32) = 0.0138 ;
loadbus(33) = 33 ; PLOAD(33) = 0.0585 ; QLOAD(33) = 0.0138 ;
loadbus(34) = 34 ; PLOAD(34) = 0.0348 ; QLOAD(34) = 0.0124 ;
loadbus(35) = 35 ; PLOAD(35) = 0 ; QLOAD(35) = 0 ;
loadbus(36) = 36 ; PLOAD(36) = 0.0223 ; QLOAD(36) = 0.0033 ;
loadbus(37) = 37 ; PLOAD(37) = 0.0102 ; QLOAD(37) = 0 ;
loadbus(38) = 38 ; PLOAD(38) = 0.0057 ; QLOAD(38) = 0.001 ;
loadbus(39) = 39 ; PLOAD(39) = 0.0113 ; QLOAD(39) = 0 ;
loadbus(40) = 40 ; PLOAD(40) = 0.0163 ; QLOAD(40) = 0.012 ;
loadbus(41) = 41 ; PLOAD(41) = 0 ; QLOAD(41) = 0 ;
loadbus(42) = 42 ; PLOAD(42) = 0.7151 ; QLOAD(42) = 0.3051 ;
loadbus(43) = 43 ; PLOAD(43) = 0.0843 ; QLOAD(43) = -0.0383 ;
loadbus(44) = 44 ; PLOAD(44) = 0 ; QLOAD(44) = 0 ;
%General parameters
%itmax = maximum number of iterations permitted before the iterative
%process is terminated – protection against infinite iterative loops
%tol = criterion tolerance to be met before the iterative solution is
%successfully brought to an end
itmax = 100;
tol = 1e-12;
nmax = 2*nbb;
%End of function PowerFlowsData
[YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,tlcond,shbus,...
shresis,shreac,ntl,nbb,nsh);
[VM,VA,it] = NewtonRaphson(nmax,tol,itmax,ngn,nld,nbb,bustype,...
genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA);
[PQsend,PQrec,PQloss,PQbus] = PQflows(nbb,ngn,ntl,nld,genbus,...
loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,...
QLOAD,VM,VA);
it; %Iteration number
VM; %Nodal voltage magnitude (p.u.)
VA = VA*180/pi; %Nodal voltage phase angle(Deg)
YR;
YI;
PQsend'; A=real(PQsend');G=imag(PQsend');
PQrec'; C=real(PQrec'); D=imag(PQrec');
PQloss';E=real(PQloss'); F=imag(PQloss');
totalPL=sum(E);totalQL=sum(F);
%End of MAIN FOR SVC SHUT VARIABLE SUSCEPTANCE
fprintf('Total iteration:%4.0f\n\n',it);
disp('=========================================================');
disp('---------------------------------------------------------');
disp(' Newton Raphson Loadflow Analysis ');
disp('--------------------------------------------------------');
disp(' | Bus | Voltage | Angle | ');
disp(' | No | pu | Degree | ');
for m = 1:nbb
disp('--------------------------------------------------------');
fprintf(' %3g', m);
fprintf(' %8.4f', VM(m));
fprintf(' %8.4f', VA(m));
fprintf('\n');
end
disp('-------------------------------------------------------------');
fprintf('\n\n\n');
disp('==================================================================================');
disp(' Line power flow and losses ');
disp('==================================================================================');
disp(' | From | To | Sending end | Receiving end | Losses | ');
disp(' | bus | bus | MW | MVAR | MW | MVAR | MW | MVAR | ');
for m=1:ntl
S=tlsend(m); R=tlrec(m);
disp('------------------------------------------------------------------------------------');
fprintf(' %4g', S); fprintf(' %4g', R); fprintf(' %8.4f', A(m));
fprintf(' %8.4f', G(m)); fprintf(' %8.4f',C(m));fprintf(' %8.4f',D(m));
fprintf(' %8.4f',E(m));fprintf(' %8.4f',F(m));
fprintf('\n');
end
disp('======================================================================================');
fprintf(' Total active and reactive power losses: ');fprintf(' \t\t%8.4f\t %8.4f\n',[totalPL,totalQL]);
disp('========================================================================================');
fprintf('\n\n');
%End Main Program
%Build up admittance matrix
function [YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,...
tlcond,shbus,shresis,shreac,ntl,nbb,nsh)
YR=zeros(nbb,nbb);
YI=zeros(nbb,nbb);
% Transmission lines contribution
for kk = 1: ntl
ii = tlsend(kk);
jj = tlrec(kk);
denom = tlresis(kk)^2+tlreac(kk)^2;
YR(ii,ii) = YR(ii,ii) + tlresis(kk)/denom + 0.5*tlcond(kk);
YI(ii,ii) = YI(ii,ii) - tlreac(kk)/denom + 0.5*tlsuscep(kk);
YR(ii,jj) = YR(ii,jj) - tlresis(kk)/denom;
YI(ii,jj) = YI(ii,jj) + tlreac(kk)/denom;
YR(jj,ii) = YR(jj,ii) - tlresis(kk)/denom;
YI(jj,ii) = YI(jj,ii) + tlreac(kk)/denom;
YR(jj,jj) = YR(jj,jj) + tlresis(kk)/denom + 0.5*tlcond(kk);
YI(jj,jj) = YI(jj,jj) - tlreac(kk)/denom + 0.5*tlsuscep(kk);
end
% Shunt elements contribution
for kk = 1: nsh
ii = shbus(kk);
denom = shresis(kk)^2+shreac(kk)^2;
YR(ii,ii) = YR(ii,ii) + shresis(kk)/denom;
YI(ii,ii) = YI(ii,ii) - shreac(kk)/denom;
end
end
% End of function YBus
%Carry out iterative solution using the Newton-Raphson method
function [VM,VA,it] = NewtonRaphson(nmax,tol,itmax,ngn,nld,nbb,...
bustype, genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA)
% GENERAL SETTINGS
D = zeros(1,nmax);
flag = 0;
it = 1;
% CALCULATE NET POWERS
[PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,QGEN,...
PLOAD,QLOAD);
while ( it < itmax && flag==0 )
% CALCULATED POWERS
[PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI);
% CHECK FOR POSSIBLE GENERATOR’S REACTIVE POWERS LIMITS VIOLATIONS
[QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,QMAX,...
QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus);
% POWER MISMATCHES
[DPQ,DP,DQ,flag] = PowerMismatches(nmax,nbb,tol,bustype,flag,PNET,...
QNET,PCAL,QCAL);
% JACOBIAN FORMATION
[JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,VM,VA,...
YR,YI);
% SOLVE FOR THE STATE VARIABLES VECTOR
D = JAC\DPQ';
% UPDATE STATE VARIABLES
[VA,VM] = StateVariablesUpdates(nbb,D,VA,VM);
it = it + 1;
end
end
% End function Newton-Raphson
%Function to calculate the net scheduled powers
function [PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,...
QGEN, PLOAD,QLOAD)
% CALCULATE NET POWERS
PNET = zeros(1,nbb);
QNET = zeros(1,nbb);
for ii = 1: ngn
PNET(genbus(ii)) = PNET(genbus(ii)) + PGEN(ii);
QNET(genbus(ii)) = QNET(genbus(ii)) + QGEN(ii);
end
for ii = 1: nld
PNET(loadbus(ii)) = PNET(loadbus(ii)) - PLOAD(ii);
QNET(loadbus(ii)) = QNET(loadbus(ii)) - QLOAD(ii);
end
end
%End function NetPowers
%Function to calculate injected bus powers
function [PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI)
% Include all entries
PCAL = zeros(1,nbb);
QCAL = zeros(1,nbb);
for ii = 1: nbb
PSUM = 0;
QSUM = 0;
for jj = 1: nbb
PSUM = PSUM + VM(ii)*VM(jj)*(YR(ii,jj)*cos(VA(ii)-VA(jj)) +...
YI(ii,jj)*sin(VA(ii)-VA(jj)));
QSUM = QSUM + VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj))-...
YI(ii,jj)*cos(VA(ii)-VA(jj)));
end
PCAL(ii) = PSUM;
QCAL(ii) = QSUM;
end
end
%End of functionCalculatePowers
%Function to check whether or not solution is within generators limits
function [QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,...
QMAX,QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus)
% CHECK FOR POSSIBLE GENERATOR’S REACTIVE POWERS LIMITS VIOLATIONS
if it > 2
flag2 = 0;
for ii = 1: ngn
jj = genbus(ii);
if (bustype(jj) == 2)
if ( QCAL(jj) > QMAX(ii) )
QNET(genbus(ii)) = QMAX(ii);
bustype(jj) = 3;
flag2 = 1;
elseif ( QCAL(jj) < QMIN(ii) )
QNET(genbus(ii)) = QMIN(ii);
bustype(jj) = 3;
flag2 = 1;
end
if flag2 == 1
for ii = 1:nld
if loadbus(ii) == jj
QNET(loadbus(ii))= QNET(loadbus(ii)) - (QLOAD(ii))
end
end
end
end
end
end
end
%End function Generatorslimits
%Function to compute power mismatches
function [DPQ,DP,DQ,flag] = PowerMismatches(nmax,nbb,tol,bustype,...
flag,PNET,QNET,PCAL,QCAL)
% POWER MISMATCHES
DPQ = zeros(1,nmax);
DP = zeros(1,nbb);
DQ = zeros(1,nbb);
DP = PNET - PCAL;
DQ = QNET - QCAL;
% To remove the active and reactive powers contributions of the slack
% bus and reactive power of all PV buses
for ii = 1: nbb
if (bustype(ii) == 1 )
DP(ii) = 0;
DQ(ii) = 0;
elseif (bustype(ii) == 2 )
DQ(ii) = 0;
end
end
% Re-arrange mismatch entries
kk = 1;
for ii = 1: nbb
DPQ(kk) = DP(ii);
DPQ(kk+1) = DQ(ii);
kk = kk + 2;
end
% Check for convergence
for ii = 1: nbb*2
if ( abs(DPQ) < tol)
flag = 1;
end
end
end
%End function PowerMismatches
%Function to built the Jacobian matrix
function [JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,...
VM,VA,YR,YI);
% JACOBIAN FORMATION
% Include all entries
JAC = zeros(nmax,nmax);
iii = 1;
for ii = 1: nbb
jjj = 1;
for jj = 1: nbb
if ii == jj
JAC(iii,jjj) = -QCAL(ii) - VM(ii)^2*YI(ii,ii);
JAC(iii,jjj+1) = PCAL(ii) + VM(ii)^2*YR(ii,ii);
JAC(iii+1,jjj) = PCAL(ii) - VM(ii)^2*YR(ii,ii);
JAC(iii+1,jjj+1) = QCAL(ii) - VM(ii)^2*YI(ii,ii);
else
JAC(iii,jjj) = VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj))...
-YI(ii,jj)*cos(VA(ii)-VA(jj)));
JAC(iii+1,jjj) = -VM(ii)*VM(jj)*(YI(ii,jj)*sin(VA(ii)...
-VA(jj))+YR(ii,jj)*cos(VA(ii)-VA(jj)));
JAC(iii,jjj+1) = -JAC(iii+1,jjj);
JAC(iii+1,jjj+1) = JAC(iii,jjj);
end
jjj = jjj + 2;
end
iii = iii + 2;
end
% Delete the voltage magnitude and phase angle equations of the slack
% bus and voltage magnitude equations corresponding to PV buses
for kk = 1: nbb
if (bustype(kk) == 1)
ii = kk*2-1;
for jj = 1: 2*nbb
if ii == jj
JAC(ii,ii) = 1;
else
JAC(ii,jj) = 0;
JAC(jj,ii) = 0;
end
end
end
if (bustype(kk) == 1) | (bustype(kk) == 2)
ii = kk*2;
for jj = 1: 2*nbb
if ii == jj
JAC(ii,ii) = 1;
else
JAC(ii,jj) = 0;
JAC(jj,ii) = 0;
end
end
end
end
end
%End of function NewtonRaphsonJacobian
%Function to update state variables
function [VA,VM] = StateVariablesUpdates(nbb,D,VA,VM)
iii = 1;
for ii = 1: nbb
VA(ii) = VA(ii) + D(iii);
VM(ii) = VM(ii) + D(iii+1)*VM(ii);
iii = iii + 2;
end
end
%End function StateVariableUpdating
%Function to calculate the power flows
function [PQsend,PQrec,PQloss,PQbus] = PQflows(nbb,ngn,ntl,nld,...
genbus,loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,...
QLOAD,VM,VA);
PQsend = zeros(1,ntl);
PQrec = zeros(1,ntl);
% Calculate active and reactive powers at the sending and receiving
% ends of tranmsission lines
for ii = 1: ntl
Vsend = ( VM(tlsend(ii))*cos(VA(tlsend(ii))) + ...
VM(tlsend(ii))*sin(VA(tlsend(ii)))*i );
Vrec = ( VM(tlrec(ii))*cos(VA(tlrec(ii))) + ...
VM(tlrec(ii))*sin(VA(tlrec(ii)))*i );
tlimped = tlresis(ii) + tlreac(ii)*i;
current =(Vsend - Vrec) / tlimped + Vsend*( tlcond(ii) + ...
tlsuscep(ii)*i )*0.5 ;
PQsend(ii) = Vsend*conj(current);
current =(Vrec - Vsend) / tlimped + Vrec*( tlcond(ii) + ...
tlsuscep(ii)*i )*0.5 ;
PQrec(ii) = Vrec*conj(current);
PQloss(ii) = PQsend(ii) + PQrec(ii);
end
% Calculate active and reactive powers injections at buses
PQbus = zeros(1,nbb);
for ii = 1: ntl
PQbus(tlsend(ii)) = PQbus(tlsend(ii)) + PQsend(ii);
PQbus(tlrec(ii)) = PQbus(tlrec(ii)) + PQrec(ii);
end
% Make corrections at generator buses, where there is load, in order to
% get correct generators contributions
for ii = 1: nld
jj = loadbus(ii);
for kk = 1: ngn
ll = genbus(kk);
if jj == ll
PQbus(jj) = PQbus(jj) + ( PLOAD(ii) + QLOAD(ii)*i );
end
end
end
end

  0 Comments

Sign in to comment.