clear all
clc
tic
iter=0;
iteration=200;
nvars = 12; %Number of variables
N = 50; %Number of Particles or Swarm size
%Acceleration constants
c1 = 2.05;
c2 = 2.05;
%Inertia Weight
w_max=0.9;
w_min=0.4;
w_temp(1)=w_max;
%Load IEEE 30-bus data
[baseMVA, bus, gen, branch]=loadcase(case_ieee30);
%% Initialization of Swarm & velocity Swarm=[unifrnd(0.95,1.10,N,6),unifrnd(0.90,1.10,N,4),unifrnd(0.00,0.20,N,2)]; %Initialize velocity
Velocity =[unifrnd(-0.003,0.003,N,6),unifrnd(-0.003,0.003,N,4), unifrnd(- 0.003,0.003,N,2)];
for i=1:N
v1=Swarm(i,1);
bus(1,8)=v1;
gen(1,6)=v1;
v2=Swarm(i,2);
bus(2,8)=v2;
gen(2,6)=v2;
v5=Swarm(i,3);
%v1
%Vm, column 8 is voltage magnitude (p.u.)
%Vg, column 6 is voltage magnitude setpoint (p.u.) %v2
%v5
85
bus(5,8)=v5;
gen(3,6)=v5;
v8=Swarm(i,4);
bus(8,8)=v8;
gen(4,6)=v8;
v11=Swarm(i,5);
bus(11,8)=v11;
gen(5,6)=v11;
v13=Swarm(i,6);
bus(13,8)=v13;
gen(6,6)=v13;
t1=Swarm(i,7);
branch(11,9)=t1;
t2=Swarm(i,8);
branch(12,9)=t2;
t3=Swarm(i,9);
branch(15,9)=t3;
t4=Swarm(i,10);
branch(36,9)=t4;
qc10=Swarm(i,11);
bus(10,6)=qc10;
qc24=Swarm(i,12);
bus(24,6)=qc24;
%v8
%v11
%v13
%tp1 6-9, column 9 is tap position
%tp2 6-10
%tp3 4-12
%tp4 27-28
%Shunt capacitor 10, column 6 is BS
%Shunt capacitor 10, column 6 is BS
eval(['savecase (''case_ieee30_test', num2str(i), '.mat'', baseMVA, bus, gen, branch)']);
eval(['initial_results_',num2str(i),'=runpf(''case_ieee30_test', num2str(i) '.mat'')']);
eval(['initial_losses_',num2str(i),'=sum(real(get_losses(initial_results_', num2str(i),')))']);
%Penalty for bus voltage violation
bus_inf=bus(:,8);
for bus_num=1:30
if bus_inf(bus_num)>1.10 penalty_Vl(bus_num)=10000*(bus_inf(bus_num)-1.10)^2;
elseif bus_inf(bus_num)<0.95 penalty_Vl(bus_num)=10000*(bus_inf(bus_num)-0.95)^2;
else
end end
end
penalty_Vl(bus_num)=0;
penalty_Vl_violation=sum(penalty_Vl);
%Penalty for shunt violation
buus_inf=bus(:,6);
for bus_num=1:30
if buus_inf(bus_num)>0.20 penalty_Qc(bus_num)=10000*(buus_inf(bus_num)-0.20)^2;
elseif buus_inf(bus_num)<0.00 penalty_Qc(bus_num)=10000*(buus_inf(bus_num)-0.00)^2;
else
end end
penalty_Qc(bus_num)=0;
penalty_Qc_violation=sum(penalty_Qc);
%Penalty for tap position violation
brch_inf=[branch(11,9); branch(12,9); branch(15,9); branch(36,9)]; for brch_num=1:4
if brch_inf(brch_num)>1.10 penalty_Tk(brch_num)=10000*(brch_inf(brch_num)-1.10)^2;
elseif brch_inf(brch_num)<0.90 penalty_Tk(brch_num)=10000*(brch_inf(brch_num)-0.90)^2;
else
penalty_Tk(brch_num)=0;
86
end
penalty_Tk_violation=sum(penalty_Tk);
%objective function=sum of active power losses of the transmission lines
losses(i)=eval(['initial_losses_',num2str(i)]); %sum of real power losses of all branches
Obj_fun_initial(i)=losses(i)+penalty_Vl_violation+penalty_Qc_violation+penalty_T k_violation; %augumented objective function with penalty function
end
%% Initialize best position (Pbest) and global best postion (Gbest) matrix Pbest=Swarm;
Val_Pbest=Obj_fun_initial;
%finding best particle in initial population
[Val_Gbest,m]=min(Val_Pbest);
Gbest=Swarm(m,:); %used to keep track of the best particle ever Gbest_calc=repmat(Swarm(m,:),N,1);
%% PSO LOOP
figure('NumberTitle', 'off', 'Name', 'PSO Algorithm Based Optimal Reactive Power Dispatch');
title('ACTIVE POWER LOSS MINIMIZATION');
ylabel('Total Active Power Loss (MW)');
xlabel('Iteration Number');
grid on;
hold on
for iter=1:iteration
% Update the value of the inertia weight w
if iter <= iteration
w_temp(iter)=w_max + (((w_min-w_max)/iteration)*iter); %Change inertia
weight
end
%generate random numbers
R1=rand(N,nvars);
R2=rand(N,nvars);
% constriction factor
%2.05+2.05=4.1;
%2/abs(2-4.1-sqrt(4.1*4.1-4*4.1))=0.729
%calculate velocity Velocity=(w_temp(iter)*Velocity+c1*R1.*(Pbest-Swarm)+c2*R2.*(Gbest_calc-
Swarm));
for v_iter=1:nvars
if v_iter==nvars Outstep=Velocity(:,v_iter)>0.003; Velocity(find(Outstep),v_iter)=0.003; Outstep=Velocity(:,v_iter)<-0.003; Velocity(find(Outstep),v_iter)=-0.003;
else
end end
%update positions of particles
Swarm=Swarm+0.729*Velocity;
for k=1:N
Outstep=Velocity(:,v_iter)>0.01; Velocity(find(Outstep),v_iter)=0.01; Outstep=Velocity(:,v_iter)<-0.01; Velocity(find(Outstep),v_iter)=-0.01;
v1=Swarm(k,1);
bus(1,8)=v1;
gen(1,6)=v1;
v2=Swarm(k,2);
bus(2,8)=v2;
gen(2,6)=v2;
v5=Swarm(k,3);
bus(5,8)=v5;
gen(3,6)=v5;
%v1
%Vm, column 8 is voltage magnitude (p.u.)
%Vg, column 6 is voltage magnitude setpoint (p.u.) %v2
%v5
%evaluate a new swarm
87
v8=Swarm(k,4);
bus(8,8)=v8;
gen(4,6)=v8;
v11=Swarm(k,5);
bus(11,8)=v11;
gen(5,6)=v11;
v13=Swarm(k,6);
bus(13,8)=v13;
gen(6,6)=v13;
t1=Swarm(k,7);
branch(11,9)=t1;
t2=Swarm(k,8);
branch(12,9)=t2;
t3=Swarm(k,9);
branch(15,9)=t3;
t4=Swarm(k,10);
branch(36,9)=t4;
qc10=Swarm(k,11);
bus(10,6)=qc10;
qc24=Swarm(k,12);
bus(24,6)=qc24;
%v8
%v11
%v13
%tp1 6-9, column 9 is tap position
%tp2 6-10
%tp3 4-12
%tp4 27-28
%Shunt capacitor 10, column 6 is BS
%Shunt capacitor 10, column 6 is BS
eval(['savecase (''case_ieee30_test', num2str(k), '.mat'', baseMVA, bus, gen, branch)']);
eval(['final_results_',num2str(k),'=runpf(''case_ieee30_test', num2str(k) '.mat'')']);
eval(['final_losses_',num2str(k),'=sum(real(get_losses(final_results_',num2str(k ),')))']);
%Penalty for bus voltage violation
bus_inf=bus(:,8);
for bus_num=1:30
if bus_inf(bus_num)>1.10 penalty_Vl(bus_num)=10000*(bus_inf(bus_num)-1.10)^2;
elseif bus_inf(bus_num)<0.95 penalty_Vl(bus_num)=10000*(bus_inf(bus_num)-0.95)^2;
else
end end
penalty_Vl(bus_num)=0;
penalty_Vl_violation=sum(penalty_Vl);
%Penalty for shunt violation
buus_inf=bus(:,6);
for bus_num=1:30
if buus_inf(bus_num)>0.20 penalty_Qc(bus_num)=10000*(buus_inf(bus_num)-0.20)^2;
elseif buus_inf(bus_num)<0.00 penalty_Qc(bus_num)=10000*(buus_inf(bus_num)-0.00)^2;
else
end end
penalty_Qc(bus_num)=0;
penalty_Qc_violation=sum(penalty_Qc);
%Penalty for tap position violation
brch_inf=[branch(11,9); branch(12,9); branch(15,9); branch(36,9)]; for brch_num=1:4
if brch_inf(brch_num)>1.10 penalty_Tk(brch_num)=10000*(brch_inf(brch_num)-1.10)^2;
elseif brch_inf(brch_num)<0.90 penalty_Tk(brch_num)=10000*(brch_inf(brch_num)-0.90)^2;
else
end end
penalty_Tk(brch_num)=0;
88
penalty_Tk_violation=sum(penalty_Tk);
%objective function=sum of active power losses of the transmission lines
losses_temp(k)=eval(['final_losses_',num2str(k)]); %sum of real power losses of all branches
Obj_fun_temp(k)=losses_temp(k)+penalty_Vl_violation+penalty_Qc_violation+penalty _Tk_violation; %augumented objective function with penalty function
% Final Evaluation
Val_Pbest_temp(k)=Obj_fun_temp(k);
end
if Val_Pbest_temp<Val_Pbest
losses=losses_temp;
Val_Pbest=Val_Pbest_temp;
Pbest=Swarm;
end
[Val_Gbest_temp,n]=min(Val_Pbest); if Val_Gbest_temp<Val_Gbest
Val_Gbest=Val_Gbest_temp; Gbest=Swarm(n,:); Gbest_calc=repmat(Swarm(n,:),N,1);
end
Val_Gbest_rec(iter)=Val_Gbest;
plot(Val_Gbest_rec);
drawnow;
toc
end
- Hello If You Notice In Equation 1 Page 13 The Objective Function Is Only Minimizing Power Loss While We Want Also To 1 (48.86 KiB) Viewed 40 times
- Hello If You Notice In Equation 1 Page 13 The Objective Function Is Only Minimizing Power Loss While We Want Also To 2 (48.86 KiB) Viewed 40 times
The Honey Bee Optimization (HBO) algorithm [60], which is a nature inspired algorithm that mimics the mating behaviour of the bee in the exploration and exploitation search, is also employed to solve ORPD. There are three kinds of bees in the colony, the queen, the workers and the drones. This technique was realised by sorting of all drones based on their fitness function. Crossover and mutation operators were applied to this technique to solve the ORPD problem. 2.2 THE ORPD PROBLEM FORMULATION The objective of the ORPD is to minimize the active power loss in the transmission network, which can be described as follows: Minimize f (x, u) (1) while satisfying (2) g(x, u)-0 h (x, u) ≤0 (3) 13 where f(x, u) is the objective function to be optimized, g(x, u) and h(x, u) are the set of equality and inequality constraints respectively. x is a vector of state variables, and u is the vector of control variables. The state variables are the load bus (PQ bus) voltages, phase angles, generator bus voltages and the slack active generation power. The control variables are the generator bus voltages, the shunt capacitors/reactors and the transformers tap settings. The objective function of the ORPD is to minimize the active power losses in the transmission lines/network, which can be defined as follows: Ni F = Min. PL=G₁ [V? + V-2V/V, cos(8) - 8)] where k refers to the branch between buses i and j: PL... and G₂ are the active loss and mutual conductance of branch k respectively; 6, and dy are the voltage angles at bus i and j; N₁ is the total number of transmission lines. The above minimization objective function is subjected to the both equality and inequality