Error in mutation part of artificial immune system code

This is my project based on Artificial Immune System (AIS). I don't know how to fix the error.

clear
clc

%varyQ_all;
tic
SLT=0;
%%%%%%%%%%%%%%%%%%%%%%%%initialization%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pop=0;
num=26;
num2=num*2;
tot_vary=8.2;
crtcl_bus=26;
max_load=18.3;
location=6;

for j=1:40000

      x02=rand*90-40;
      x05=rand*80-40;
      x08=rand*50-10;
      x011=rand*30-6;
      x013=rand*30-6;
      x110=rand*50;
      x124=rand*50;
      nu02=rand*5;
      nu05=rand*5;
      nu08=rand*5;
      nu011=rand*5;
      nu013=rand*5;
      nu110=rand*5;
      nu124=rand*5;
     [basemva,accuracy,maxiter,busdata,linedata]=IEEE30BUSDATA;
     busdata(2,6)=busdata(2,6)-x02;
     busdata(5,6)=busdata(5,6)-x05;
     busdata(8,6)=busdata(8,6)-x08;
     busdata(11,6)=busdata(11,6)-x011;
     busdata(13,6)=busdata(13,6)-x013;
     busdata(10,7)=x110;
     busdata(24,7)=x124;
     lfybus
     lfnewton
     busout
     lineflow
     SLT=0;
     pp=(Vm)';
     loss=real(SLT);
     Pg1=Pg(1);
     if pp>=0.95 & pp<=1.10 & loss<=15.0 & Pg(1)<=500 & Pg(1)>=100
         loop2=0;
         for increment2=[1.0:0.05:1000.0]
             loop2=loop2+1;
              [basemva,accuracy,maxiter,busdata,linedata]=IEEE30BUSDATA;
   busdata(crtcl_bus,location)=busdata(crtcl_bus,location)*tot_vary;
   busdata(crtcl_bus,location)=busdata(crtcl_bus,location)*increment2;
              busdata(2,6)=busdata(2,6)-x02;
              busdata(5,6)=busdata(5,6)-x05;
              busdata(8,6)=busdata(8,6)-x08;
              busdata(11,6)=busdata(11,6)-x011;
              busdata(13,6)=busdata(13,6)-x013;
              busdata(10,7)=x110;
              busdata(24,7)=x124;
     lfybus
     lfnewton
     busout
     lineflow
     SLT=0;
     a=min(Vm);
     b=busdata(crtcl_bus,location);
     if a<0.85
         break
     end
     load2=busdata(crtcl_bus,location);
     Ploss2=real(SLT);
     Vmin2=min(Vm);
     Vmax2=max(Vm);
     Load2a(loop2,:)=busdata(crtcl_bus,location);
     Ploss2a(loop2,:)=real(SLT);
     Vmin2a(loop2,:)=min(Vm);
     Vmax2a(loop2,:)=max(Vm);
         end
         if b>=20.13
             pop=pop+1;
             parents(pop,:)=[x02 x05 x08 x011 x013 x110 x124...
                 nu02 nu05 nu08 nu011 nu013 nu110 nu124...
                 Ploss2 Vmin2 Vmax2]
             if pop==num
                 break
             end
         end
     end
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CLONE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            x02=parents(:,1);
            x05=parents(:,2);
            x08=parents(:,3);
            x011=parents(:,4);
            x013=parents(:,5);
            x110=parents(:,6);
            x124=parents(:,7);
            nu02=parents(:,8);
            nu05=parents(:,9);
            nu08=parents(:,10);
            nu011=parents(:,11);
            nu013=parents(:,12);
            nu110=parents(:,13);
            nu124=parents(:,14);
            load2=parents(:,15);
            Ploss2=parents(:,16);
            Vmin2=parents(:,17);
            Vmax2=parents(:,17);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MUTATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mutatn_iter=0;
for big_loop=1:40000
    mutatn_iter=mutatn_iter+1;
    tau=1/(sqrt(2*sqrt(10)));
    taum=1/(sqrt(2*10));

      for loop3=1:num;
          for j2=1:40000  %loop j2 is to make sure that the mutated value comply to the constraint
        gauss1
              x02_m;
              x05_m;
              x08_m;
              x011_m;
              x013_m;
              x110_m;
              x124_m;
               [basemva,accuracy,maxiter,busdata,linedata]=IEEE30BUSDATA;
                  busdata(2,6)=busdata(2,6)-x02_mutate;
                  busdata(5,6)=busdata(5,6)-x05_mutate;
                  busdata(8,6)=busdata(8,6)-x08_mutate;
                  busdata(11,6)=busdata(11,6)-x011_mutate;
                  busdata(13,6)=busdata(13,6)-x013_mutate;
                  busdata(10,7)=x110_mutate;
                  busdata(24,7)=x124_mutate;
        lfybus
        lfnewton
        busout
        lineflow
        SLT=0;
        qq=(Vm)';
        loss2=real(SLT);
        Pg1a=Pg(1);
  if qq>=0.95 & qq<=1.10 & loss2<=15.0 & Pg(1)<=500 & Pg(1)>=100
                loop4=0;
                for increment4=[1.0:0.05:1000.0]
                    loop4=loop4+1;
              [basemva,accuracy,maxiter,busdata,linedata]=IEEE30BUSDATA;
   busdata(crtcl_bus,location)=busdata(crtcl_bus,location)*tot_vary;
   busdata(crtcl_bus,location)=busdata(crtcl_bus,location)*increment4;
                      busdata(2,6)=busdata(2,6)-x02_mutate;
                  busdata(5,6)=busdata(5,6)-x05_mutate;
                  busdata(8,6)=busdata(8,6)-x08_mutate;
                  busdata(11,6)=busdata(11,6)-x011_mutate;
                  busdata(13,6)=busdata(13,6)-x013_mutate;
                  busdata(10,7)=x110_mutate;
                  busdata(24,7)=x124_mutate;
        lfybus
        lfnewton
        busout
        lineflow
        fvsi
        maxloadability
        c=min(Vm);
        d=busdata(crtcl_bus,location);
        if c<0.85
            break
        end
        Ploss4=real(SLT);
        Vmin4=min(Vm);
        Vmax4=max(Vm);
        load4=busdata(crtcl_bus,loacation);
        load4a(loop4,:)=busdata(crtcl_bus,loacation);
        Ploss4a(loop4,:)=real(SLT);
        Vmin4a(loop4,:)=min(Vm);
        Vmax4a(loop4,:)=max(Vm);
                end
                if d>=20.13
                    offspring(loop3,:)=[x02_mutate x05_mutate x08_mutate...
                        x011_mutate x013_mutate x110_mutate x124_mutate...
                        nu02_mutate nu05_mutate nu08_mutate nu011_mutate...
                        nu013_mutate nu110_mutate nu124_mutate load4...
                        Ploss4 Vmin4 Vmax4]
                    break
                end
  end
          end
      end
x02_mutate=offspring(:,1);
x05_mutate=offspring(:,2);
x08_mutate=offspring(:,3);
x011_mutate=offspring(:,4);
x013_mutate=offspring(:,5);

x110_mutate=offspring(:,6);
x124_mutate=offspring(:,7);

nu02_mutate=offspring(:,8);
nu05_mutate=offspring(:,9);
nu08_mutate=offspring(:,10);
nu011_mutate=offspring(:,11);
nu013_mutate=offspring(:,12);

nu110_mutate=offspring(:,13);
nu124_mutate=offspring(:,14);

load4=offspring(:,15);
Ploss4=offspring(:,16);
Vmin4=offspring(:,17);
Vmax4=offspring(:,18);;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%SELECTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

population=[x02 x05 x08 x011 x013 x110 x124 nu02 nu05 nu08 nu011 nu013...
            nu110 nu124 load2 Ploss2 Vmin2 Vmax2;x02_mutate x05_mutate...
            x08_mutate x011_mutate x013_mutate x110_mutate x124_mutate...
            nu02_mutate nu05_mutate nu08_mutate nu011_mutate...
            nu013_mutate nu110_mutate nu124_mutate...
            load4 Ploss4 Vmin4 Vmax4];

          rank1=sortrows(population,15);
          rank2=flipud(rank1);
          nxt_genrtn=rank2(1:num,:)
          x02=nxt_genrtn(:,1);
          x05=nxt_genrtn(:,2);
          x08=nxt_genrtn(:,3);
          x011=nxt_genrtn(:,4);
          x013=nxt_genrtn(:,5);
          x110=nxt_genrtn(:,6);
          x124=nxt_genrtn(:,7);
          nu02=nxt_genrtn(:,8);
          nu05=nxt_genrtn(:,9);
          nu08=nxt_genrtn(:,10);
          nu011=nxt_genrtn(:,11);
          nu013=nxt_genrtn(:,12);
          nu110=nxt_genrtn(:,13);
          nu124=nxt_genrtn(:,14);
          load2=nxt_genrtn(:,15);
          Ploss2=nxt_genrtn(:,16);
          Vmin2=nxt_genrtn(:,17);
          Vmax2=nxt_genrtn(:,18);
          mx_load=max(load2);
          mn_load=min(load2);
          diff=mx_load-mn_load
          pause(3)
          if diff<0.0001
              break
          end
  end
toc
result1=[x02 x05 x08 x011 x013 x110 x124 nu02 nu05 nu08 nu011 nu013...
         nu110 nu124 load2 Ploss2 Vmin2 Vmax2]
     result=sortrows(result1,15)
     best_reading=result(1,:)
     parents1=parents
     mutatn_iter1=mutatn_iter

Answer 1

[Walter Roberson]

I'm having trouble imagining what the IEEE 30 Bus has to do with Artificial Immune Systems, but anyhow...

I either do not have the toolbox(es) needed to run that myself or else you have omitted some necessary routines or variables.

As I am not able to run the program myself to see what error comes up, you will need to describe the errors you are seeing and where they occur in the code and how the results differ from what you expect.

[Ahmad Firdaus]

thnks for your concern.the error at the mutation part is it says that "Undefined function or variable 'x02_m'".here i attach the IEEE30 bus system and other toolbox incase needed.

for IEEE30 BUS:

function  [basemva,accuracy,maxiter,busdata,linedata]=IEEE30BUSDATA;
clear 
clc

basemva = 100; accuracy = 0.001; acce1 = 1.8;maxiter = 100;
%   IEEE 30-BUS TEST SYSTEM (American Electric Power)
%      Bus Bus Voltage Angle --Load--   ---Generator--- Injected
%       No code Mag. Degree MW    Mvar     MW     Mvar Qmin  Qmax Mvar
busdata=[1  1   1.06    0   0.0    0.0     0.0    0.0  0    0       0
         2  2   1.043   0   21.70  12.7    40.0   0.0  -40  50      0      
         3  0   1.0     0   2.4     1.2     0.0   0.0  0    0       0  
         4  0   1.06    0   7.6     1.6     0.0   0.0  0    0       0
         5  2   1.01    0   94.2    19.0    0.0   0.0   -40 40      0  
         6  0   1.0     0   0.0     0.0     0.0   0.0   0   0       0
         7  0   1.0     0   22.8    10.9    0.0   0.0   0   0       0
         8  2   1.01    0   30.0    30.0    0.0   0.0   -10 40      0    
         9  0   1.0     0   0.0     0.0     0.0   0.0   0   0       0 
         10 0   1.0     0   5.8     2.0     0.0   0.0   0   0       19
         11 2   1.082   0   0.0     0.0     0.0   0.0   -6  24      0
         12 0   1.0     0   11.2    7.5     0     0     0   0       0  
         13 2   1.071   0   0.0     0.0     0     0     -6  24      0
         14 0   1.0     0   6.2     1.6     0     0     0   0       0
         15 0   1.0     0   8.2     2.5     0     0     0   0       0
         16 0   1.0     0   3.5     1.8     0     0     0   0       0
         17 0   1.0     0   9.0     5.8     0     0     0   0       0
         18 0   1.0     0   3.2     0.9     0     0     0   0       0
         19 0   1.0     0   9.5     3.4     0     0     0   0       0
         20 0   1.0     0   2.2     0.7     0     0     0   0       0
         21 0   1.0     0   17.5    11.2    0     0     0   0       0
         22 0   1.0     0   0.0     0.0     0     0     0   0       0
         23 0   1.0     0   3.2     1.6     0     0     0   0       0
         24 0   1.0     0   8.7     6.7     0     0     0   0       4.3
         25 0   1.0     0   0.0     0.0     0     0     0   0       0
         26 0   1.0     0   3.5     2.3     0     0     0   0       0
         27 0   1.0     0   0.0     0.0     0     0     0   0       0
         28 0   1.0     0   0.0     0.0     0     0     0   0       0
         29 0   1.0     0   2.4     0.9     0     0     0   0       0
         30 0   1.0     0   10.6    1.9     0     0     0   0      0];

% Line Data
%
%         Bus     bus     R        X        1/2 B for Line code or
%         n1      nr      pu       pu       pu      tap setting value
linedata=[1       2       0.0192    0.0575  0.02640     1  
          1       3       0.0452    0.1852  0.02040     1  
          2       4       0.0570    0.1737  0.01840     1  
          3       4       0.0132    0.0379  0.00420     1  
          2       5       0.0472    0.1983  0.02090     1  
          2       6       0.0581    0.1763  0.01870     1  
          4       6       0.0119    0.0414  0.00450     1  
          5       7       0.0460    0.1160  0.01020     1  
          6       7       0.0267    0.0820  0.00850     1  
          6       8       0.0120    0.0420  0.00450     1
          6       9       0.0       0.2080  0.0         0.978
          6       10      0.0       0.5560  0.0         0.969  
          9       11      0.0       0.2080  0.0         1  
          9       10      0.0       0.1100  0.0         1      
          4       12      0.0       0.2560  0.0         0.932     
          12      13      0.0       0.1400  0.0         1  
          12      14      0.1231    0.2559  0.0         1  
          12      15      0.0662    0.1304  0.0         1  
          12      16      0.0945    0.1987  0.0         1  
          14      15      0.2210    0.1997  0.0         1
          16      17      0.0824    0.1923  0.0         1  
          15      18      0.1073    0.2185  0.0         1 
          18      19      0.0639    0.1292  0.0         1    
          19      20      0.0340    0.0680  0.0         1
          10      20      0.0936    0.2090  0.0         1  
          10      17      0.0324    0.0845  0.0         1
          10      21      0.0348    0.0749  0.0         1  
          10      22      0.0727    0.1499  0.0         1  
          21      22      0.0116    0.0236  0.0         1
          15      23      0.1000    0.2020  0.0         1  
          22      24      0.1150    0.1790  0.0         1
          23      24      0.1320    0.2700  0.0         1
          24      25      0.1885    0.3292  0.0         1  
          25      26      0.2544    0.3800  0.0         1
          25      27      0.1093    0.2087  0.0         1  
          28      27      0.0000    0.3960  0.0         0.968  
          27      29      0.2198    0.4153  0.0         1
          27      30      0.3202    0.6027  0.0         1  
          29      30      0.2399    0.4533  0.0         1
          8       28      0.0636    0.2000  0.0214      1  
          6       28      0.0169    0.0599  0.065       1];
%
    lfybus                   % Forms the bus admittance matrix
    lfnewton      % Power flow solution by Newton-Raphson method
    busout        % Prints the power flow solution on the screen
    lineflow      % computes and displays the line flow and losses

for lfgauss:

%   Power flow solution by Gauss-Seidel method
%   Copyright (c) 1998 by  H. Saadat
% Revision 1 (Aug. 99) Modified to include two or more parallel lines
Vm=0; delta=0; yload=0; deltad =0;
nbus = length(busdata(:,1));
kb=[];Vm=[]; delta=[]; Pd=[]; Qd=[]; Pg=[]; Qg=[]; Qmin=[]; Qmax=[];  % Added (6-8-00)
Pk=[]; P=[]; Qk=[]; Q=[]; S=[]; V=[]; % Added (6-8-00)
for k=1:nbus
n=busdata(k,1);
kb(n)=busdata(k,2); Vm(n)=busdata(k,3); delta(n)=busdata(k, 4);
Pd(n)=busdata(k,5); Qd(n)=busdata(k,6); Pg(n)=busdata(k,7); Qg(n) = busdata(k,8);
Qmin(n)=busdata(k, 9); Qmax(n)=busdata(k, 10);
Qsh(n)=busdata(k, 11);
    if Vm(n) <= 0  
        Vm(n) = 1.0; V(n) = 1 + j*0;
    else delta(n) = pi/180*delta(n);
         V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));
         P(n)=(Pg(n)-Pd(n))/basemva;
         Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;
         S(n) = P(n) + j*Q(n);
    end
DV(n)=0;
end
num = 0; AcurBus = 0; converge = 1;
Vc = zeros(nbus,1)+j*zeros(nbus,1); Sc = zeros(nbus,1)+j*zeros(nbus,1);

while exist('accel')~=1
   accel = 1.3;
end
while exist('accuracy')~=1
   accuracy = 0.001;
end
while exist('basemva')~=1
   basemva= 100;
end
while exist('maxiter')~=1
   maxiter = 100;
end
%%%%added for parallel lines (Aug. 99)
mline=ones(nbr,1);
for k=1:nbr
      for m=k+1:nbr
         if((nl(k)==nl(m)) && (nr(k)==nr(m)));
            mline(m)=2;
         elseif ((nl(k)==nr(m)) && (nr(k)==nl(m)));
         mline(m)=2;
         else, end
      end
   end

end of statements for parallel lines (Aug. 99)

iter=0;
maxerror=10;
while maxerror >= accuracy && iter <= maxiter
iter=iter+1;
  for n = 1:nbus;
  YV = 0+j*0;
    for L = 1:nbr;
            if (nl(L) == n && mline(L) == 1), k=nr(L);   %modified to handle parallel lines (Aug. 99)
            YV = YV + Ybus(n,k)*V(k);
            elseif (nr(L) == n && mline(L)==1), k=nl(L); %modified to handle parallel lines (Aug. 99)
        YV = YV + Ybus(n,k)*V(k);
            end
    end
       Sc = conj(V(n))*(Ybus(n,n)*V(n) + YV) ;
       Sc = conj(Sc);
       DP(n) = P(n) - real(Sc);
       DQ(n) = Q(n) - imag(Sc);
         if kb(n) == 1
         S(n) =Sc; P(n) = real(Sc); Q(n) = imag(Sc); DP(n) =0; DQ(n)=0;
         Vc(n) = V(n);
         elseif kb(n) == 2
         Q(n) = imag(Sc); S(n) = P(n) + j*Q(n);

             if Qmax(n) ~= 0
               Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
               if abs(DQ(n)) <= .005 && iter >= 10 % After 10 iterations
                 if DV(n) <= 0.045                % the Mvar of generator buses are
                    if Qgc < Qmin(n),             % tested. If not within limits Vm(n)
                    Vm(n) = Vm(n) + 0.005;        % is changed in steps of 0.005 pu
                    DV(n) = DV(n)+.005;           % up to .05  pu in order to bring
                    elseif Qgc > Qmax(n),         % the generator Mvar within the
                    Vm(n) = Vm(n) - 0.005;        % specified limits.
                    DV(n)=DV(n)+.005; end
                 else, end
               else,end
             else,end
           end
         if kb(n) ~= 1
         Vc(n) = (conj(S(n))/conj(V(n)) - YV )/ Ybus(n,n);
         else, end
            if kb(n) == 0
            V(n) = V(n) + accel*(Vc(n)-V(n));
            elseif kb(n) == 2
            VcI = imag(Vc(n));
            VcR = sqrt(Vm(n)^2 - VcI^2);
            Vc(n) = VcR + j*VcI;
             V(n) = V(n) + accel*(Vc(n) -V(n));
            end
     end
    maxerror=max( max(abs(real(DP))), max(abs(imag(DQ))) );
     if iter == maxiter && maxerror > accuracy
     fprintf('\nWARNING: Iterative solution did not converged after ')
     fprintf('%g', iter), fprintf(' iterations.\n\n')
     fprintf('Press Enter to terminate the iterations and print the results \n')
     converge = 0; pause, else, end

end
if converge ~= 1
   tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 
   tech=('                   Power Flow Solution by Gauss-Seidel Method');
end   
k=0;
for n = 1:nbus
  Vm(n) = abs(V(n)); deltad(n) = angle(V(n))*180/pi;
     if kb(n) == 1
     S(n)=P(n)+j*Q(n);
     Pg(n) = P(n)*basemva + Pd(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
     k=k+1;
     Pgg(k)=Pg(n);
     elseif  kb(n) ==2
     k=k+1;
     Pgg(k)=Pg(n);
     S(n)=P(n)+j*Q(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
     end
yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
end
Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);
busdata(:,3)=Vm'; busdata(:,4)=deltad';
clear  AcurBus  DP  DQ  DV  L Sc Vc VcI VcR YV converge delta

for lineflow:

%  This program is used in conjunction with lfgauss or lf Newton
%  for the computation of line flow and line losses.
%
%  Copyright (c) 1998 H. Saadat
SLT = 0;
fprintf('\n')
fprintf('                           Line Flow and Losses \n\n')
fprintf('     --Line--  Power at bus & line flow    --Line loss--  Transformer\n')
fprintf('     from  to    MW      Mvar     MVA       MW      Mvar      tap\n')
for n = 1:nbus
busprt = 0;
   for L = 1:nbr;
       if busprt == 0
       fprintf('   \n'), fprintf('%6g', n), fprintf('      %9.3f', P(n)*basemva)
       fprintf('%9.3f', Q(n)*basemva), fprintf('%9.3f\n', abs(S(n)*basemva))

         busprt = 1;
         else, end
         if nl(L)==n      k = nr(L);
         In = (V(n) - a(L)*V(k))*y(L)/a(L)^2 + Bc(L)/a(L)^2*V(n);
         Ik = (V(k) - V(n)/a(L))*y(L) + Bc(L)*V(k);
         Snk = V(n)*conj(In)*basemva;
         Skn = V(k)*conj(Ik)*basemva;
         SL  = Snk + Skn;
         SLT = SLT + SL;
         elseif nr(L)==n  k = nl(L);
         In = (V(n) - V(k)/a(L))*y(L) + Bc(L)*V(n);
         Ik = (V(k) - a(L)*V(n))*y(L)/a(L)^2 + Bc(L)/a(L)^2*V(k);
         Snk = V(n)*conj(In)*basemva;
         Skn = V(k)*conj(Ik)*basemva;
         SL  = Snk + Skn;
         SLT = SLT + SL;
         else, end
           if nl(L)==n | nr(L)==n
           fprintf('%12g', k),
           fprintf('%9.3f', real(Snk)), fprintf('%9.3f', imag(Snk))
           fprintf('%9.3f', abs(Snk)),
           fprintf('%9.3f', real(SL)),
               if nl(L) ==n & a(L) ~= 1
               fprintf('%9.3f', imag(SL)), fprintf('%9.3f\n', a(L))
               else, fprintf('%9.3f\n', imag(SL))
               end
           else, end
    end
  end
  SLT = SLT/2;
  fprintf('   \n'), fprintf('    Total loss                         ')
  fprintf('%9.3f', real(SLT)), fprintf('%9.3f\n', imag(SLT))
  clear Ik In SL SLT Skn Snk

for lfnewton:

%   Power flow solution by Newton-Raphson method
%   Copyright (c) 1998 by  H. Saadat
%   Revision 1 (Aug. 99) To include two or more parallel lines
ns=0; ng=0; Vm=0; delta=0; yload=0; deltad=0;
nbus = length(busdata(:,1));
kb=[];Vm=[]; delta=[]; Pd=[]; Qd=[]; Pg=[]; Qg=[]; Qmin=[]; Qmax=[];  % Added (6-8-00)
Pk=[]; P=[]; Qk=[]; Q=[]; S=[]; V=[]; % Added (6-8-00)
for k=1:nbus
n=busdata(k,1);
kb(n)=busdata(k,2); Vm(n)=busdata(k,3); delta(n)=busdata(k, 4);
Pd(n)=busdata(k,5); Qd(n)=busdata(k,6); Pg(n)=busdata(k,7); Qg(n) = busdata(k,8);
Qmin(n)=busdata(k, 9); Qmax(n)=busdata(k, 10);
Qsh(n)=busdata(k, 11);
    if Vm(n) <= 0  
        Vm(n) = 1.0; V(n) = 1 + j*0;
    else delta(n) = pi/180*delta(n);
         V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));
         P(n)=(Pg(n)-Pd(n))/basemva;
         Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;
         S(n) = P(n) + j*Q(n);
    end
end
for k=1:nbus
if kb(k) == 1, ns = ns+1; else, end
if kb(k) == 2 
    ng = ng+1; else, end
ngs(k) = ng;
nss(k) = ns;
end
Ym=abs(Ybus); t = angle(Ybus);
m=2*nbus-ng-2*ns;
maxerror = 1; converge=1;
iter = 0;
%%%%added for parallel lines (Aug. 99)
mline=ones(nbr,1);
for k=1:nbr
      for m=k+1:nbr
         if((nl(k)==nl(m)) && (nr(k)==nr(m)));
            mline(m)=2;
         elseif ((nl(k)==nr(m)) && (nr(k)==nl(m)));
         mline(m)=2;
         else, end
      end
   end

end of statements for parallel lines (Aug. 99)

% Start of iterations
clear A  DC   J  DX
while maxerror >= accuracy && iter <= maxiter % Test for max. power mismatch
for ii=1:m
for k=1:m
   A(ii,k)=0;      %Initializing Jacobian matrix
end, end
iter = iter+1;
for n=1:nbus
nn=n-nss(n);
lm=nbus+n-ngs(n)-nss(n)-ns;
J11=0; J22=0; J33=0; J44=0;
   for ii=1:nbr
     if mline(ii)==1   % Added to include parallel lines (Aug. 99)
        if nl(ii) == n || nr(ii) == n
           if nl(ii) == n ,  l = nr(ii); end
           if nr(ii) == n , l = nl(ii); end
         J11=J11+ Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));
         J33=J33+ Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));
            if kb(n)~=1
            J22=J22+ Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));
            J44=J44+ Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));
            else, end
            if kb(n) ~= 1  && kb(l) ~=1
              lk = nbus+l-ngs(l)-nss(l)-ns;
              ll = l -nss(l);
            % off diagonalelements of J1
              A(nn, ll) =-Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));
                  if kb(l) == 0  % off diagonal elements of J2
                  A(nn, lk) =Vm(n)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));end
                  if kb(n) == 0  % off diagonal elements of J3
                  A(lm, ll) =-Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n)+delta(l)); end
                  if kb(n) == 0 && kb(l) == 0  % off diagonal elements of  J4
                  A(lm, lk) =-Vm(n)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));end
             else end
         else , end
      else, end   
   end
   Pk = Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;
   Qk = -Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;
   if kb(n) == 1 
       P(n)=Pk; Q(n) = Qk; end   % Swing bus P
     if kb(n) == 2  
         Q(n)=Qk;
         if Qmax(n) ~= 0
           Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
           if iter <= 7                  % Between the 2th & 6th iterations
              if iter > 2                % the Mvar of generator buses are
                if Qgc  < Qmin(n),       % tested. If not within limits Vm(n)
                Vm(n) = Vm(n) + 0.01;    % is changed in steps of 0.01 pu to
                elseif Qgc  > Qmax(n),   % bring the generator Mvar within
                Vm(n) = Vm(n) - 0.01;end % the specified limits.
              else, end
           else,end
         else,end
     end
   if kb(n) ~= 1
     A(nn,nn) = J11;  %diagonal elements of J1
     DC(nn) = P(n)-Pk;
   end
   if kb(n) == 0
     A(nn,lm) = 2*Vm(n)*Ym(n,n)*cos(t(n,n))+J22;  %diagonal elements of J2
     A(lm,nn)= J33;        %diagonal elements of J3
     A(lm,lm) =-2*Vm(n)*Ym(n,n)*sin(t(n,n))-J44;  %diagonal of elements of J4
     DC(lm) = Q(n)-Qk;
   end
end
DX=A\DC';
for n=1:nbus
  nn=n-nss(n);
  lm=nbus+n-ngs(n)-nss(n)-ns;
    if kb(n) ~= 1
    delta(n) = delta(n)+DX(nn); end
    if kb(n) == 0
    Vm(n)=Vm(n)+DX(lm); end
 end
  maxerror=max(abs(DC));
     if iter == maxiter && maxerror > accuracy 
   fprintf('\nWARNING: Iterative solution did not converged after ')
   fprintf('%g', iter), fprintf(' iterations.\n\n')
   fprintf('Press Enter to terminate the iterations and print the results \n')
   converge = 0; pause, else, end

end

if converge ~= 1
   tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 
   tech=('                   Power Flow Solution by Newton-Raphson Method');
end   
V = Vm.*cos(delta)+j*Vm.*sin(delta);
deltad=180/pi*delta;
i=sqrt(-1);
k=0;
for n = 1:nbus
     if kb(n) == 1
     k=k+1;
     S(n)= P(n)+j*Q(n);
     Pg(n) = P(n)*basemva + Pd(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
     Pgg(k)=Pg(n);
     Qgg(k)=Qg(n);     %june 97
     elseif  kb(n) ==2
     k=k+1;
     S(n)=P(n)+j*Q(n);
     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
     Pgg(k)=Pg(n);
     Qgg(k)=Qg(n);  % June 1997
  end
yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
end
busdata(:,3)=Vm'; busdata(:,4)=deltad';
Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);

%clear A DC DX  J11 J22 J33 J44 Qk delta lk ll lm
%clear A DC DX  J11 J22 J33  Qk delta lk ll lm

for lfybus:

%  This program obtains th Bus Admittance Matrix for power flow solution
%  Copyright (c) 1998 by H. Saadat

j=sqrt(-1); i = sqrt(-1);
nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);
X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);
nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));
Z = R + j*X; y= ones(nbr,1)./Z;        %branch admittance
for n = 1:nbr
if a(n) <= 0  
    a(n) = 1; else end
Ybus=zeros(nbus,nbus);     % initialize Ybus to zero
               % formation of the off diagonal elements
for k=1:nbr;
       Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);
       Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));
    end
end
              % formation of the diagonal elements
for  n=1:nbus
     for k=1:nbr
         if nl(k)==n
         Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);
         elseif nr(k)==n
         Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);
         else, end
     end
end
clear Pgg

that's all.. i also have problem in fvsi part.. i hope u can help me.this is for my project in degree of Electrical Engineering. btw thnks

[Walter Roberson]

Nothing in the code you show assigns anything to x02_m . You do have code that assigns to x02_mutate -- did you rename variables and miss a block of code?

It is difficult to tell for sure looking at the code, but it appears to me that you may be trying to use x02_mutate and related variables before you assign any value to them. You assign values to them after some loops (or perhaps after some "if" statements) but you use the values inside the loops.

[Ahmad Firdaus]

so i should assign the x02_m into some value first?i rename the variables x02_m to indicate that it is the mutation part different from the original which is only x02.

[Walter Roberson]

Variables must be assigned a value before they can be used.

That particular block of code does not do anything useful with those variables: a variable on a line by itself would normally mean to display the variable, but the semi-colon means not to print out the value. Thus what _I_ would do is comment out that section.