i am using matpower in matlab to determine optimal switching schedule for feeder capacitors using GA optimisation based on 24 hours dayahead forecasted load profile for every bus load. the objective fxn is minimisation of I2R loss, subjected to voltage constraint (ie bus voltage should lie within a dead band) and the switching operation should be less than 7 times a day for each capacitor.load profile is a 30*48 (bus hourly active,reactive load) read from xlsx file. my network has 6 capacitors,hence my decision variable should be a string of 6*24=144 switch status, representing 24 status of switch for each capacitor in each hour of a day.
To test the code i am supposing 100% loading of all the bus loads
throughout the day thus imposing a constant load on the substation
transformer all day. this means the result of optimisation should be
theoritically constant switch status of all feeder capacitors all day (as
load hasn't changed). after running the GA optimisation the result shows
change in switch status of capacitors (though within the constraint) which
i beieve should be same all day as the loading is constant at every bus
throughout the day.
I even tried all loads operating at pf 1 all day; still it shows several
switchings of capacitors which i beleive should have been at status '0' for
pf 1.
Has the caseformat i used contained any bug? Please suggest. I have also
attached the the zip file with this post.
%%case30bus
function mpc = case30original(hr, cap1, cap2, cap3, cap4, cap5, cap6)
if cap1==1
c1=0.6;
else c1=0;
end
if cap2==1
c2=0.6;
else c2=0;
end
if cap3==1
c3=0.6;
else c3=0;
end
if cap4==1
c4=0.3;
else c4=0;
end
if cap5==1
c5=0.9;
else c5=0;
end
if cap6==1
c6=0.9;
else c6=0;
end
%% Read Load Profile
global lp;
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 100;
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va
baseKV zone Vmax Vmin
mpc.bus = [
1 2 0 0 0 0 1 1 0 23 1 1.2 0.6;
2 1 lp(2,2*hr-1) lp(2,2*hr)-c1 0 0 1 1 0 23 1 1.2 0.6;
3 1 0 0 0 0 1 1 0 23 1 1.2 0.6;
4 1 lp(4,2*hr-1) lp(4,2*hr) 0 0 1 1 0 23 1 1.2 0.6;
5 1 0 0 0 0 1 1 0 23 1 1.2 0.6;
6 1 0 0 0 0 1 1 0 23 1 1.2 0.6;
7 1 0 0 0 0 1 1 0 23 1 1.2 0.6;
8 1 0 0 0 0 1 1 0 23 1 1.2 0.6;
9 1 lp(9,2*hr-1) lp(9,2*hr) 0 0 1 1 0 23 1 1.2 0.6;
10 1 0 0 0 0 3 1 0 23 1 1.2 0.6;
11 1 lp(11,2*hr-1) lp(11,2*hr) 0 0 1 1 0 23 1 1.2 0.6;
12 1 lp(12,2*hr-1) lp(12,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
13 1 lp(13,2*hr-1) lp(13,2*hr)-c2 0 0 2 1 0 23 1 1.2 0.6;
14 1 lp(14,2*hr-1) lp(14,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
15 1 lp(15,2*hr-1) lp(15,2*hr)-c3 0 0 2 1 0 23 1 1.2 0.6;
16 1 lp(16,2*hr-1) lp(16,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
17 1 lp(17,2*hr-1) lp(17,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
18 1 lp(18,2*hr-1) lp(18,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
19 1 lp(19,2*hr-1) lp(19,2*hr)-c4 0 0 2 1 0 23 1 1.2 0.6;
20 1 lp(20,2*hr-1) lp(20,2*hr) 0 0 2 1 0 23 1 1.2 0.6;
21 1 lp(21,2*hr-1) lp(21,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
22 1 lp(22,2*hr-1) lp(22,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
23 1 0 0-c5 0 0 2 1 0 23 1 1.2 0.6;
24 1 lp(24,2*hr-1) lp(24,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
25 1 lp(25,2*hr-1) lp(25,2*hr)-c6 0 0 3 1 0 23 1 1.2 0.6;
26 1 lp(26,2*hr-1) lp(26,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
27 1 lp(27,2*hr-1) lp(27,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
28 1 lp(28,2*hr-1) lp(28,2*hr) 0 0 1 1 0 23 1 1.2 0.6;
29 1 lp(29,2*hr-1) lp(29,2*hr) 0 0 3 1 0 23 1 1.2 0.6;
30 1 lp(30,2*hr-1) lp(30,2*hr) 0 0 3 1 0 23 1 1.2 0.6
];
%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax
Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10
ramp_30 ramp_q apf
mpc.gen=[
1 15 0 0 0 1 100 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
%22 0 0 0 0 0 100 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%28 0 0 0 0 0 100 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
%27 26.91 0 48.7 -15 1 100 0 55 0 0 0 0 0 0 0 0 0 0 0 0;
%23 19.2 0 40 -10 1 100 0 30 0 0 0 0 0 0 0 0 0 0 0 0;
%13 37 0 44.7 -15 1 100 0 40 0 0 0 0 0 0 0 0 0 0 0 0;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio
angle status angmin angmax
mpc.branch=[
1 2 0.11 0.32 0 130 130 130 1 0 1 -360 360;
2 3 0.07 0.07 0 130 130 130 0 0 1 -360 360;
3 4 0.22 0.19 0 65 65 65 0 0 1 -360 360;
4 5 0.10 0.09 0 130 130 130 0 0 1 -360 360;
5 6 0.31 0.18 0 130 130 130 0 0 1 -360 360;
6 7 0.26 0.14 0 65 65 65 0 0 1 -360 360;
7 8 0.26 0.14 0 90 90 90 0 0 1 -360 360;
8 9 0.25 0.14 0 70 70 70 0 0 1 -360 360;
9 10 0.25 0.14 0 130 130 130 0 0 1 -360 360;
10 11 0.75 0.42 0 32 32 32 0 0 1 -360 360;
11 12 0.35 0.2 0 65 65 65 0 0 1 -360 360;
12 13 0.14 0.08 0 32 32 32 0 0 1 -360 360;
13 14 0.29 0.16 0 65 65 65 0 0 1 -360 360;
8 15 0.09 0.08 0 65 65 65 0 0 1 -360 360;
15 16 0.14 0.08 0 65 65 65 0 0 1 -360 360;
16 17 0.25 0.14 0 65 65 65 0 0 1 -360 360;
6 18 0.09 0.14 0 32 32 32 0 0 1 -360 360;
18 19 0.30 0.08 0 32 32 32 0 0 1 -360 360;
19 20 0.29 0.26 0 32 32 32 0 0 1 -360 360;
6 21 0.11 0.16 0 16 16 16 0 0 1 -360 360;
3 22 0.11 0.10 0 16 16 16 0 0 1 -360 360;
22 23 0.06 0.06 0 16 16 16 0 0 1 -360 360;
23 24 0.11 0.09 0 16 16 16 0 0 1 -360 360;
24 25 0.28 0.24 0 32 32 32 0 0 1 -360 360;
25 26 0.20 0.17 0 32 32 32 0 0 1 -360 360;
26 27 0.29 0.16 0 32 32 32 0 0 1 -360 360;
2 28 0.09 0.00 0 32 32 32 0 0 1 -360 360;
28 29 0.31 0.17 0 32 32 32 0 0 1 -360 360;
29 30 0.21 0.12 0 32 32 32 0 0 1 -360 360;
];
%% fitness_function
function y=fitness(x)
y=0;
global Simu_Hour;
%Simu_Hour=24;
for j=1:Simu_Hour;
temp = case30original(j,x(j),x(24+j),x(48+j),x(72+j),x(96+j),x(120+j));
T=runpf(temp);
for b=1:29
y=y+T.branch(b,14)+T.branch(b,16);
end
end
%%CONSTRAINTS FUNCTION
function [c,ceq] = constraints(x)
global Simu_Hour;
global bus_vol_hr;
c1_tap=0;c2_tap=0;c3_tap=0;c4_tap=0;c5_tap=0;c6_tap=0;
for j=1:Simu_Hour-1
if (x(j)~=x(j+1))
c1_tap = c1_tap+1;
else
end
if (x(24+j)~=x(25+j))
c2_tap=c2_tap+1;
else
end
if (x(48+j)~=x(49+j))
c3_tap= c3_tap+1;
else
end
if (x(72+j)~=x(73+j))
c4_tap=c4_tap+1;
else
end
if (x(96+j)~=x(97+j))
c5_tap=c5_tap+1;
else
end
if (x(120+j)~=x(121+j))
c6_tap=c6_tap+1;
else
end
end %for loop end
V_day=zeros(1,24);
Vmax=zeros(1,24);
Vmin=zeros(1,24);
for j=1:Simu_Hour
temp=case30original(j,x(j),x(24+j),x(48+j),x(72+j),x(96+j),x(120+j));
T=runpf(temp);
for b=1:30
V_day(b)=T.bus(b,8);
bus_vol_hr(j,b)= V_day(b);
end
Vmax(j)=max(V_day);
Vmin(j)=min(V_day);
end
vol=[
-1*Vmin(1)-(-0.9);Vmax(1)-1.1;
-1*Vmin(2)-(-0.9);Vmax(2)-1.1;
-1*Vmin(3)-(-0.9);Vmax(3)-1.1;
-1*Vmin(4)-(-0.9);Vmax(4)-1.1;
-1*Vmin(5)-(-0.9);Vmax(5)-1.1;
-1*Vmin(6)-(-0.9);Vmax(6)-1.1;
-1*Vmin(7)-(-0.9);Vmax(7)-1.1;
-1*Vmin(8)-(-0.9);Vmax(8)-1.1;
-1*Vmin(9)-(-0.9);Vmax(9)-1.1;
-1*Vmin(10)-(-0.9);Vmax(10)-1.1;
-1*Vmin(11)-(-0.9);Vmax(11)-1.1;
-1*Vmin(12)-(-0.9);Vmax(12)-1.1;
-1*Vmin(13)-(-0.9);Vmax(13)-1.1;
-1*Vmin(14)-(-0.9);Vmax(14)-1.1;
-1*Vmin(15)-(-0.9);Vmax(15)-1.1;
-1*Vmin(16)-(-0.9);Vmax(16)-1.1;
-1*Vmin(17)-(-0.9);Vmax(17)-1.1;
-1*Vmin(18)-(-0.9);Vmax(18)-1.1;
-1*Vmin(19)-(-0.9);Vmax(19)-1.1;
-1*Vmin(20)-(-0.9);Vmax(20)-1.1;
-1*Vmin(21)-(-0.9);Vmax(21)-1.1;
-1*Vmin(22)-(-0.9);Vmax(22)-1.1;
-1*Vmin(23)-(-0.9);Vmax(23)-1.1;
-1*Vmin(24)-(-0.9);Vmax(24)-1.1
];
capacitor_switching=[
c1_tap-7;...
c2_tap-7;...
c3_tap-7;...
c4_tap-7;...
c5_tap-7;...
c6_tap-7
];
%%All NON-LINEAR CONSTRAINTS
c = [vol;capacitor_switching];
%%No equality constraints
ceq=[];
%% main_body
tic;
clear all;
clc;
nvars=144;
LB=zeros(1,144);
UB=ones(1,144);
global Simu_Hour;
global lp;
global bus_vol_hr;
bus_vol_hr = zeros(24,30);
Simu_Hour=24;
lp=xlsread('networkdata.xlsx','busdata','F2:BA31');% Read Load Profile
ObjectiveFunction = @fitness;
ConstraintFunction = @constraints;
opts= gaoptimset(...
'PopulationSize',300,...
'Generations',600,...
'EliteCount',20,...
'StallGenLimit',30,...
'TolFun',1e-8,...
'PlotFcns',@gaplotbestf);
[x,fbest,exitflag]= ga(ObjectiveFunction,nvars,[],[],[],[],...
LB,UB,ConstraintFunction,1:144,opts);
display(x);
toc;
load100%_activeonly.rar
Description: application/rar
