i am using matpower in matlab to determine optimal switching schedule for feeder capacitors using GA optimisation to minimise loss with less capacitor switching operations but the best function value is not optimal.

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. 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).

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 in a day. 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. Please suggest if the code has any bugs. I have also attached the the zip file with this post.

if true
%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;
end