clear allclcticiter=0;iteration=200;nvars = 12; %Number of variablesN = 50; %Number of Particles or Swarm size%Accelerat

Business, Finance, Economics, Accounting, Operations Management, Computer Science, Electrical Engineering, Mechanical Engineering, Civil Engineering, Chemical Engineering, Algebra, Precalculus, Statistics and Probabilty, Advanced Math, Physics, Chemistry, Biology, Nursing, Psychology, Certifications, Tests, Prep, and more.
Post Reply
answerhappygod
Site Admin
Posts: 899603
Joined: Mon Aug 02, 2021 8:13 am

clear allclcticiter=0;iteration=200;nvars = 12; %Number of variablesN = 50; %Number of Particles or Swarm size%Accelerat

Post by answerhappygod »

clear allclcticiter=0;iteration=200;nvars = 12; %Number of variablesN = 50; %Number of Particles or Swarm size%Acceleration constantsc1 = 2.05;c2 = 2.05;%Inertia Weightw_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 velocityVelocity =[unifrnd(-0.003,0.003,N,6),unifrnd(-0.003,0.003,N,4), unifrnd(- 0.003,0.003,N,2)];for i=1:Nv1=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%v585bus(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 BSeval(['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:30if 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;elseend endendpenalty_Vl(bus_num)=0;penalty_Vl_violation=sum(penalty_Vl);%Penalty for shunt violationbuus_inf=bus(:,6);for bus_num=1:30if 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;elseend endpenalty_Qc(bus_num)=0;penalty_Qc_violation=sum(penalty_Qc);%Penalty for tap position violationbrch_inf=[branch(11,9); branch(12,9); branch(15,9); branch(36,9)]; for brch_num=1:4if 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;elsepenalty_Tk(brch_num)=0;86endpenalty_Tk_violation=sum(penalty_Tk);%objective function=sum of active power losses of the transmission lineslosses(i)=eval(['initial_losses_',num2str(i)]); %sum of real power losses of all branchesObj_fun_initial(i)=losses(i)+penalty_Vl_violation+penalty_Qc_violation+penalty_T k_violation; %augumented objective function with penalty functionend%% 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 LOOPfigure('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 onfor iter=1:iteration % Update the value of the inertia weight wif iter <= iterationw_temp(iter)=w_max + (((w_min-w_max)/iteration)*iter); %Change inertiaweightend %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:nvarsif 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;elseend end %update positions of particles Swarm=Swarm+0.729*Velocity; for k=1:NOutstep=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 swarm87v8=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 BSeval(['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:30if 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;elseend endpenalty_Vl(bus_num)=0;penalty_Vl_violation=sum(penalty_Vl);%Penalty for shunt violationbuus_inf=bus(:,6);for bus_num=1:30if 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;elseend endpenalty_Qc(bus_num)=0;penalty_Qc_violation=sum(penalty_Qc);%Penalty for tap position violationbrch_inf=[branch(11,9); branch(12,9); branch(15,9); branch(36,9)]; for brch_num=1:4if 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;elseend endpenalty_Tk(brch_num)=0;88penalty_Tk_violation=sum(penalty_Tk);%objective function=sum of active power losses of the transmission lineslosses_temp(k)=eval(['final_losses_',num2str(k)]); %sum of real power losses of all branchesObj_fun_temp(k)=losses_temp(k)+penalty_Vl_violation+penalty_Qc_violation+penalty _Tk_violation; %augumented objective function with penalty function % Final EvaluationVal_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_GbestVal_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; tocend
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
Join a community of subject matter experts. Register for FREE to view solutions, replies, and use search function. Request answer by replying!
Top
Post Reply

Return to “Archived Questions & Answers”