Code covered by the BSD License

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

### kartik pandya (view profile)

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

```