Optimal Reactive Power Dispatch(ORPD) using Particle Swarm Optimization(PSO)

by

 

12 Apr 2013 (Updated )

I have solved optimal reactive power dispatch problem using Particle Swarm Optimization method

ORPD_30bus()
function [fxmin, xmin] = ORPD_30bus()
clc
clear all
warning off all
global Ybus Yf Yt;

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

[baseMVA, bus, gen, branch, areas, gencost]=case_ieee30; %Input data of IEEE 30 bus test system

% Determine the value of weight change
w_start = 1;   %Initial inertia weight's value
w_end = 0.20;       %Final inertia weight
w_varyfor = floor(0.7*100); 
w_now = w_start;
inertdec = (w_start-w_end)/w_varyfor; %Inertia weight's change per iteration
iter=0;

% Initialize size of Swarm, no. of control varialbes and Velocity
%size of swarm=50 (user may change it), no. of control variables=vg1, vg2,
%vg5,vg8, vg11, vg13 (Generator bus voltage in pu), T1, T2, T3, T4
%(Transformer tap position), QC3, QC10, QC24 (Injected reactive power of
%capacitors in MVAR)

Swarm=[unifrnd(0.90,1.10,50,6),unifrnd(0.95,1.05,50,4),unifrnd(1,20,50,3)]; 
VStep =[unifrnd(0.04,0.04,50,6),unifrnd(0.002,0.002,50,4), unifrnd(4,4,50,3)];
for i=1:50 % no of particles
    global bus gen branch v1 v2 v5 v8 v11 v13
    v1=Swarm(i,1);
    v2=Swarm(i,2);
    v5=Swarm(i,3);
    v8=Swarm(i,4);
    v11=Swarm(i,5);
    v13=Swarm(i,6);
    bus(1,VM)=Swarm(i,1); 
    bus(2,VM)=Swarm(i,2);
    bus(5,VM)=Swarm(i,3);
    bus(8,VM)=Swarm(i,4);
    bus(11,VM)=Swarm(i,5);
    bus(13,VM)=Swarm(i,6);
    gen(1,VG)=Swarm(i,1); % gen1
    gen(2,VG)=Swarm(i,2); %gen2
    gen(3,VG)=Swarm(i,3);  %gen5
    gen(4,VG)=Swarm(i,4);  %gen8
    gen(5,VG)=Swarm(i,5);  %gen11
    gen(6,VG)=Swarm(i,6);  %gen13
    
    branch(11,TAP)=Swarm(i,7);% T1, transformer tap
    branch(12,TAP)=Swarm(i,8); %T2
    branch(15,TAP)=Swarm(i,9); %T3
    branch(36,TAP)=Swarm(i,10); %T4
    bus(3,BS)=Swarm(i,11); % QC3, capacitor
    bus(10,BS)=Swarm(i,12); % QC10, capacitor
    bus(24,BS)=Swarm(i,13); % QC24, capacitor
   
         global bus gen branch Yf Yt Ybus
        [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);% construct Ybus
        [MVAbase,bus,gen,branch,success,et]=runpf; % run NR load flow
         global V
    global Ybus
    bbb=branch(:,[PF,QF]).^2; % ans. is P(from)^2, Q(from)^2
    ccc=bbb'; 
    ddd=(sum(ccc))';
    eee=sqrt(ddd);  %actual values of apparent powers "From bus injection" 
    
    fff=branch(:,[PT,QT]).^2;
    ggg=fff';
    hhh=(sum(ggg))';
    iii=sqrt(hhh); %actual values of apparent powers "To bus injection" 
    jjj=[];
    for count_1=1:41 % jjj will contain actual max. apparent power.("from" bus power is compared with "to" bus power)
        if eee(count_1)>=iii(count_1)
            jjj(count_1)=eee(count_1);
        else
            jjj(count_1)=iii(count_1);
        end
    end
    
   %penalty value calculation for line flow violations
    for count_2=1:41 % max power is compared with thermal limit
      if jjj(count_2)>branch(count_2, RATE_A)
          kkk(count_2)=10000*((jjj(count_2)-branch((count_2), RATE_A))^2); % if limit violates,penalty is imposed
      else
          kkk(count_2)=0; %if no violation, penalty is zero
      end
    end
    penalty_line=sum(kkk);% summation of all penalty of line flow violations
    
  
        %penalty value calculation for bus voltage violation
    lll=bus(:,VM); % voltage mag. of each bus
   for count_3=1:30
        if lll(count_3)>bus(count_3,VMAX)
            penalty_volt(count_3)=10000*(lll(count_3)-bus(count_3,VMAX))^2;
        elseif lll(count_3)<bus(count_3,VMIN)
            penalty_volt(count_3)=10000*(bus(count_3,VMIN)-lll(count_3))^2;
        else
            penalty_volt(count_3)=0;
        end
   end
   penalty_volt=sum(penalty_volt);%summation of penalty for bus voltage violation

    %penalty value calculation for reactive power violation of all
    %generators(PV buses and slack bus)
    mmm=gen(:,QG);% reactive power output of all generators
    for count_4=1:6
        if mmm(count_4)>gen(count_4,QMAX)
            penalty_reactive(count_4)=10000*((mmm(count_4)-gen(count_4,QMAX))^2);
        elseif mmm(count_4)<gen(count_4,QMIN)
            penalty_reactive(count_4)=10000*((gen(count_4,QMIN)-mmm(count_4))^2);
        else
            penalty_reactive(count_4)=0;
        end
    end
    penalty_reactive=sum(penalty_reactive);%summation of penalty for Q limit violation
    
    %penalty value calculation for P violation of slack generator
    nnn=gen(:,PG);
    if nnn(1,1)>gen(1,PMAX)
        penalty_slack=1000*((nnn(1,1)-gen(1,PMAX))^2);
    elseif nnn(1,1)<gen(1,PMIN)
        penalty_slack=1000*((gen(1,PMIN)-nnn(1,1))^2);
    else
        penalty_slack=0;
    end
  %sum of all penalties, penalty function is used to handle inequality
  %constraints
  global penalty_function
 penalty_function=(penalty_slack)+ (penalty_reactive)+(penalty_volt)+(penalty_line);
  %sum of real power losses of all branches
global loss
sum_P_loss=(sum(real(loss)));
global penalty_function
fun_Swarm=(sum_P_loss)+(penalty_function);
% Objective function= sum of active power losses of the transmission lines
fSwarm(i,:)=sum_P_loss;
end 

% Initializing the Best positions matrix and
% the corresponding function values
PBest = Swarm; 
fPBest = fSwarm;
% Finding best particle in initial population
[fGBest, g] = min(fSwarm);
lastbpf = fGBest;
Best = Swarm(g,:); %Used to keep track of the Best particle ever
fBest = fGBest;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  THE  PSO  LOOP                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while( iter < 200 )% no of iterations
    global iter
    iter = iter+1;
    if iter>190
        diary on
    end
    
    % Update the value of the inertia weight w
    if (iter<=w_varyfor) & (iter > 1)
        w_now = w_now - inertdec; %Change inertia weight
    end 
   
    % The PLAIN PSO %
    % Set GBest
    A = repmat(Swarm(g,:), 50, 1); %A = GBest. repmat(X, m, n) repeats the matrix X in m rows by n columns, 50=no. of particles       
    % Generate Random Numbers
    R1 = rand(50, 13); % 50=No. of particles, 13=dimension of the problem(size of swarm=no. of variables)
    R2 = rand(50, 13);
    
    % Calculate Velocity
    VStep = w_now*VStep + 2.1*R1.*(PBest-Swarm) + 2.0*R2.*(A-Swarm);% c1=2.1, c2=2
    
    % Apply Vmax Operator for v > Vmax
    changeRows = VStep > 0.003; 
    VStep(find(changeRows)) = 0.003;
    % Apply Vmax Operator for v < -Vmax
    changeRows = VStep < -0.003; 
    VStep(find(changeRows)) = -0.003;
    
    % ::UPDATE POSITIONS OF PARTICLES::
    Swarm = Swarm + 0.729 * VStep;  % Evaluate new Swarm, Chi=0.729
  
   for j=1:50 %no of particles
    v1=Swarm(j,1);
    v2=Swarm(j,2);
    v5=Swarm(j,3);
    v8=Swarm(j,4);
    v11=Swarm(j,5);
    v13=Swarm(j,6);
    bus(1,VM)=Swarm(j,1); 
    bus(2,VM)=Swarm(j,2);
    bus(5,VM)=Swarm(j,3);
    bus(8,VM)=Swarm(j,4);
    bus(11,VM)=Swarm(j,5);
    bus(13,VM)=Swarm(j,6);
    
    gen(1,VG)=Swarm(j,1); % gen1
    gen(2,VG)=Swarm(j,2); %gen2
    gen(3,VG)=Swarm(j,3);  %gen5
    gen(4,VG)=Swarm(j,4);  %gen8
    gen(5,VG)=Swarm(j,5);  %gen11
    gen(6,VG)=Swarm(j,6);  %gen13
    
    branch(11,TAP)=Swarm(j,7);% T1, transformer tap
    branch(12,TAP)=Swarm(j,8); %T2
    branch(15,TAP)=Swarm(j,9); %T3
    branch(36,TAP)=Swarm(j,10); %T4
    bus(3,BS)=Swarm(j,11); % QC3, capacitor
    bus(10,BS)=Swarm(j,12); % QC10, capacitor
    bus(24,BS)=Swarm(j,13); % QC24, capacitor
   
        [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
        [MVAbase,bus,gen,branch,success,et]=runpf;
         global V
   
    bbb=branch(:,[PF,QF]).^2; % ans. is P(from)^2, Q(from)^2
    ccc=bbb'; 
    ddd=(sum(ccc))';
    eee=sqrt(ddd);  %actual values of apparent powers "From bus injection" 
    
    fff=branch(:,[PT,QT]).^2;
    ggg=fff';
    hhh=(sum(ggg))';
    iii=sqrt(hhh); %actual values of apparent powers "To bus injection" 
    jjj=[];
    for count_1=1:41 % jjj will contain actual max. apparent power.("from" bus power is compared with "to" bus power)
        if eee(count_1)>=iii(count_1)
            jjj(count_1)=eee(count_1);
        else
            jjj(count_1)=iii(count_1);
        end
    end
    
   %penalty value calculation for line flow violations
    for count_2=1:41 % max power is compared with thermal limit
      if jjj(count_2)>branch(count_2, RATE_A)
          kkk(count_2)=10000*((jjj(count_2)-branch((count_2), RATE_A))^2); % if limit violates,penalty is imposed
      else
          kkk(count_2)=0; %if no violation, penalty is zero
      end
    end
    penalty_line=sum(kkk);% summation of all penalty of line flow violations
    
  
        %penalty value calculation for bus voltage violation
    lll=bus(:,VM); % voltage mag. of each bus
   for count_3=1:30
        if lll(count_3)>bus(count_3,VMAX)
            penalty_volt(count_3)=10000*(lll(count_3)-bus(count_3,VMAX))^2;
        elseif lll(count_3)<bus(count_3,VMIN)
            penalty_volt(count_3)=10000*(bus(count_3,VMIN)-lll(count_3))^2;
        else
            penalty_volt(count_3)=0;
        end
   end
   penalty_volt=sum(penalty_volt);%summation of penalty for bus voltage violation

    %penalty value calculation for reactive power violation of all
    %generators(PV buses and slack bus)
    mmm=gen(:,QG);% reactive power output of all generators
    for count_4=1:6
        if mmm(count_4)>gen(count_4,QMAX)
            penalty_reactive(count_4)=10000*((mmm(count_4)-gen(count_4,QMAX))^2);
        elseif mmm(count_4)<gen(count_4,QMIN)
            penalty_reactive(count_4)=10000*((gen(count_4,QMIN)-mmm(count_4))^2);
        else
            penalty_reactive(count_4)=0;
        end
    end
    penalty_reactive=sum(penalty_reactive);%summation of penalty for Q limit violation
    
    %penalty value calculation for P violation of slack generator
    nnn=gen(:,PG);
    if nnn(1,1)>gen(1,PMAX)
        penalty_slack=1000*((nnn(1,1)-gen(1,PMAX))^2);
    elseif nnn(1,1)<gen(1,PMIN)
        penalty_slack=1000*((gen(1,PMIN)-nnn(1,1))^2);
    else
        penalty_slack=0;
    end
  %sum of all penalties  
  global penalty_function
 penalty_function=(penalty_slack)+ (penalty_reactive)+(penalty_volt)+(penalty_line);
  %sum of real power losses of all branches
 global loss
sum_P_loss=(sum(real(loss)))
global penalty_function
fun_Swarm=(sum_P_loss)+(penalty_function);
fSwarm(j,:)=sum_P_loss;
end
   
    % Updating the best position for each particle
    changeRows = fSwarm < fPBest;                      % latest fswarm-previous fswarm
    fPBest(find(changeRows)) = fSwarm(find(changeRows)); %fitness value
    PBest(find(changeRows), :) = Swarm(find(changeRows), :); %position
    lastbpart = PBest(g, :);
    % Updating index g
    [fGBest, g] = min(fPBest);
      
    %Update Best. Only if fitness has improved.
    if fGBest < lastbpf
        [fBest, b] = min(fPBest);
        fbest(iter) = fBest; 
        figure(1);
     
    plot(fbest, '-b' );
    title('Total Active power losses (MW)');
    xlabel('Iteration number');
    ylabel('Total Active power losses (MW)');
    fBest;
    Best = PBest(b,:);
    end
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  END  OF PSO  LOOP                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[fxmin, b] = min(fPBest)
xmin = PBest(b, :)
diary off
return

Contact us