clc
clear
pnum = 30;
maxit =300 ;
%---------------------initializing DE parameters--------------%

CR=0.9;
F=0.75;
%initfit = 10000000000;
%targfit = 0;
dim = 4;
Pd = [0 0 0 0 0 25 75 450 650 650 650 650 650 650 675 650 625 575 125 0 0 0 0 0 ];
Ph = [0 0 0 0 0 0 100 600 600 600 600 600 600 600 700 600 650 350 0 0 0 0 0 0 ];
M1=min(Pd);
%M2=max(Pd);
%ES_Cap=700;
%Charge_Limit=8*ES_Cap;
for(i=1:24)
xlowt(i,1) = 100;
xhight(i,1) = 325;
xhight(i,3)=180;
xlowt(i,3)=0;
xlowt(i,4)=0;
xhight(i,4)=660;
end
 for i=1:24
      if(Pd(i)>=325)
          xlowt(i,2)=Pd(i)-325;
          xhight(i,2)=xlowt(i,2);
      else
          xhight(i,2)=-150;
          xlowt(i,2)=-50; 
      end
 end

%----------------------Preparation of DE particles---------------%
optimizing=0;
for(i=1:pnum+3)
    for(j=1:24)
        while(3>0)
         sump =0;
         for(k=2:dim-1)
            ind(i,j,k)=xlowt(j,k)+(xhight(j,k)-xlowt(j,k))*rand(1);
            sump=sump+ind(i,j,k);
         end
            ind(i,j,1)=Pd(j)-sump;
            ash=Ph(j)-(0.5*(ind(i,j,1)));
            if(ash>0)
                ind(i,j,4)=ash;
            else
                ind(i,j,4)=0;
            end
                
         %sump = sump + ind(i,j); 
         %end
            if(ind(i,j,1)>=xlowt(j,1) & ind(i,j,1)<=xhight(j,1))                   
                break;
            else
                optimizing=optimizing+1;
            continue;
            end
        end   
    end
end


for(i=1:pnum)
     for(j=1:24)
         for(k=1:dim)
             w(i,j,k)=ind(i+3,j,k)+F*(ind(i+1,j,k)-ind(i+2,j,k));
             if(rand(1)<CR)
                 t(i,j,k)=w(i,j,k);
             else
                 t(i,j,k)=ind(i,j,k);
             end
         end
     end
end



mincost=0;
minind(1,24,4)=0;
fitness(pnum)=0;                                                          
fitness_f(pnum)=0;
p=1 ;
while( p<=maxit+1 )
    for(i=1:pnum)
     fitness(i)=costfunc(ind(i,:,:));
     fitness_f(i)=costfunc(t(i,:,:));
              if(fitness(i)<fitness_f(i))
                ind(i,:,:)=ind(i,:,:);
              else
                ind(i,:,:)=t(i,:,:);
              end
              cost(i)=costfunc(ind(i,:,:));
    end
    
    for(i=2:pnum)
        if(cost(i)<cost(i-1))
            mincost=cost(i);
            minind(1,:,:)=ind(i,:,:);
          
        else
            mincost=cost(i-1);
            minind(1,:,:)=ind(i-1,:,:);
           
        end
    end
    p=p+1;
end

cost=0;
disp('Economic Load sharing among the generators are calculated to be as follows : ');
ans_Generation=minind(1,:,:);
display(ans_Generation)
display(mincost)