I can't minimize fval value using this code

global  ncap
global LDorig BDorig

BDorig = [1    0    0 
2    100  60
3    90   40
4    120  80
5    60   30
6    60   20
7    200  100
8    200  100 
9    60   20
10   60   20
11   45   30
12   60   35
13   60   35
14   120  80
15   60   10
16   60   20
17   60   20
18   90   40
19   90   40
20   90   40
21   90   40
22   90   40
23   90   50
24   420  200
25   420  200
26   60   25
27   60   25
28   60   20
29   120  70
30   200  600
31   150  70
32   210  100
33   60   40];

LDorig = [01   2   0.0922   0.0470
  02   3   0.4930   0.2511
  03   4   0.3661   0.1864 
  04   5   0.3811   0.1941 
  05   6   0.8190   0.7070 
  06   7   0.1872   0.6188
  07   8   0.7115   0.2351
  08   9   1.0299   0.7400
  09   10  1.0440   0.7400 
  10   11  0.1967   0.0651
  11   12  0.3744   0.1298
  12   13  1.4680   1.1549
  13   14  0.5416   0.7129
  14   15  0.5909   0.5260
  15   16  0.7462   0.5449
  16   17  1.2889   1.7210
  17   18  0.7320   0.5739
  02   19  0.1640   0.1565
  19   20  1.5042   1.3555
  20   21  0.4095   0.4784
  21   22  0.7089   0.9373
  3    23  0.4512   0.3084
  23   24  0.8980   0.7091
  24   25  0.8959   0.7071
  6    26  0.2031   0.1034
  26   27  0.2842   0.1447 
  27   28  1.0589   0.9338
  28   29  0.8043   0.7006
  29   30  0.5074   0.2585
  30   31  0.9745   0.9629
  31   32  0.3105   0.3619 
  32   33  0.3411   0.5302];


[sLOSS,v]=power_flow1()
ncap=4;
min=0.01;
max=0.8;
for i=1:2*ncap
k=mod(i,2);
    if k==0
        lb(i,1)=min;
        ub(i,1)=max;
    else
        lb(i,1)=2;
           ub(i,1)=33;
    end
end
Vmin=0.95;
Vmax=1.05;

options = gaoptimset;
options = gaoptimset('PopulationSize', 50,'Generations', 1500,'StallGenLimit',1000,'TimeLimit', 5000,'StallTimeLimit', 1000,'PlotFcn',@gaplotbestf);

constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmax, Vmin);
rng default
[capplcsz_opt,fval,exitflag]=ga(@power_flow_cap,2*ncap,[],[],[],[],lb,ub,constraints,[1 3 5 7  ],options);
capplcsz_opt([1 3 5 7])
[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)
bar(Busvoltage,'r');
axis([1 33 0 1.5]);
grid on
figure(2);
bar(abs(v),'b');
axis([1 33 0 1.5]);
grid on

function [g,h]=constraint_GA(capplcsz_opt, Vmin, Vmax)
g=[];
h=[];
[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)
V=abs(Busvoltage);
 g = [V(:) - Vmax; Vmin - V(:)];

%idx = find(g<0)

end

function [PLOOS,Busvoltage]=power_flow1()
global LDorig BDorig

%display('entrer la puissance de base du system en MVA');
Sbase=10;
%display('entrer la tension de base du system en KV');
Vbase=12.66;
%display('entrer la tolerance epsilon');
epsilon=0.00001;
%display('entrer le nombre d"iteration T');
T=100;
Zbase=Vbase^2/Sbase;
BD=BDorig;%load('bdata33bus.m');
LD=LDorig;%load('ldata33busmod4.m');
cc=sum(BD(:,3))


BD(:,2:3)=BD(:,2:3)/(Sbase*10^3); %%per unit initialization
LD(:,3:4)=LD(:,3:4)/Zbase; %%per unit initialization
dimBD=size(BD);
dimLD=size(LD);
Nbus=dimBD(1,1);%%%% nombre de jeux de barres
Nbr=dimLD(1,1); %%%% nombre de branche
impdance=zeros(Nbus); %%% matrice d'impedance
courant=zeros(Nbus); %%%% matrice courant
terminbus=zeros(Nbus,1); %%% noeuds fineaux
interbus=zeros(Nbus,1); %%%% bus intermedier
for i=1:Nbr
            impdance(LD(i,1),LD(i,2))=complex(LD(i,3),LD(i,4)); %for intermediate
end
for k=1:Nbus 
 co=0; 
 l=0; 
for n=1:Nbr 
if LD(n,1)==k 
 co=co+1; 
end
end
if co==0 
 terminbus(k,1)=k; 
elseif co>= 1
for m=1:Nbr 
 l=l+1; 
if LD(m,2)==k 
 interbus(k,1)=k; 
break
elseif l==Nbr 
 refbus=k; 
end
end
end
end
Busvoltage=ones(Nbus,1);
%Busvoltage(refbus,1)=1.05*Busvoltage(1,1);
iteration=0;
for s=1:T %iteration
                      %%%%%%%%%%%%%%%%%%%
                      %%%backward sweep%%
                      %%%%%%%%%%%%%%%%%%%
for i=1:Nbus
    if terminbus(i,1)~=0

    nodenumber=terminbus(i,1);
    
    for j=1:Nbus
        if impdance(j,nodenumber)~=0
            connectedfrom=j;
        
        end
    end
    courant(connectedfrom,nodenumber)=conj(complex(BD(nodenumber,2),BD(nodenumber,3))/Busvoltage(nodenumber,1));
    end
end

for i=Nbus:-1:1%for intermediate nodes
    if interbus(i,1)~=0
            for k=1:Nbus %find which upper node is connected to current intermediate node
                if impdance(k,interbus(i,1))~=0 %&& k~=IntermediateNodes(i,1)
                    upperNode=k;
                    break;
                end
            end
        courant(upperNode,interbus(i,1))=(conj(complex(BD(interbus(i,1),2),BD(interbus(i,1),3))/Busvoltage(interbus(i,1),1)));%current equal to load or generation connected to intermediate node
            for i1=1:Nbus%current to intermediate node is equal to currents going from inter node to all connected downstream nodes + current of node itself if it exists
                if impdance(interbus(i,1),i1)~=0 && i1~=interbus(i,1)
                    courant(upperNode,interbus(i,1))=courant(upperNode,interbus(i,1))+courant(interbus(i,1),i1);
                end
            end
          if upperNode==refbus
          end
        

    end

end
                     %%%%%%%%%%%%%%%%%
                     %%forward sweep%%
                     %%%%%%%%%%%%%%%%%
                    
 voltage(:,1)=Busvoltage(:,1);

     for i=1:Nbus
    for j=1:Nbus
        if impdance(i,j)~=0
            Busvoltage(j,1)=Busvoltage(i,1)-courant(i,j)*impdance(i,j);
        end
    end
end
count=0;
p=0;
for p=2:1:Nbus 
if abs(voltage(p)-Busvoltage(p))<=epsilon
 count=count+1; 
end
end
if count==(Nbus-1) 
break
end
iteration=iteration+1;
end

                               %%%%%% l'écoulement de puissance%%%%
 bilan=zeros(Nbr,3);
 c=0;
 for i=1:Nbus
    
     for j=1:Nbus
         if courant(i,j)~=0
             c=c+1;
              bilan(c,1)=i;%depart 
              bilan(c,2)=j;%arrive
    bilan(c,3)=Busvoltage(i)*conj(courant(i,j))+Busvoltage(j)*conj(-1*courant(i,j)); %les perts par la methode Sloss=Sij+Sji
         end
     end
 end
sLOSS=sum(bilan(:,3));
Qloss=imag(sLOSS);
 PLOOS=real(sLOSS);
 Voltag=abs(Busvoltage);
% disp(Voltag);
 

end

function [Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)
global  ncap
global LDorig BDorig


%global ldata_o ncap capplc;
%display('entrer la puissance de base du system en MVA');
Sbase=10;
%display('entrer la tension de base du system en KV');
Vbase=12.66;
%display('entrer la tolerance epsilon');
epsilon=0.0001;
%display('entrer le nombre d"iteration T');
T=200;
Zbase=(Vbase^2)/Sbase;
BD=BDorig;%load('bdata33bus.m');
LD=LDorig;%load('ldata33busmod4.m');
%bp=BD(:,2);
%bq=BD(:,3);
BD(:,2:3)=BD(:,2:3)./(Sbase*10^3); %%per unit initialization
LD(:,3:4)=LD(:,3:4)./Zbase; %%per unit initialization
%bppu=bp./Sbase;
%bqpu=bq./Sbase;
%BD(:,2)=bppu;
%BD(:,3)=bqpu;
ncap=4;
%capplc=[18];
%for tt=1:ncap
   for tt=1:ncap
    BD(capplcsz_opt(1,2*tt-1),3)=BD(capplcsz_opt(1,2*tt-1),3)-capplcsz_opt(1,2*tt);
end
%end
dimBD=size(BD);
dimLD=size(LD);
Nbus=dimBD(1,1);%%%% nombre de jeux de barres
Nbr=dimLD(1,1); %%%% nombre de branche
impdance=zeros(Nbus); %%% matrice d'impedance
courant=zeros(Nbus); %%%% matrice courant
terminbus=zeros(Nbus,1); %%% noeuds fineaux
interbus=zeros(Nbus,1); %%%% bus intermedier
for i=1:Nbr
            impdance(LD(i,1),LD(i,2))=complex(LD(i,3),LD(i,4)); %for intermediate
end
for k=1:Nbus 
 co=0; 
 l=0; 
for n=1:Nbr 
if LD(n,1)==k 
 co=co+1; 
end
end
if co==0 
 terminbus(k,1)=k; 
elseif co>= 1
for m=1:Nbr 
 l=l+1; 
if LD(m,2)==k 
 interbus(k,1)=k; 
break
elseif l==Nbr 
 refbus=k; 
end
end
end
end
Busvoltage=ones(Nbus,1);
%Busvoltage(refbus,1)=1.05*Busvoltage(1,1);
%tty=Busvoltage;
iteration=0;
for s=1:T %iteration
                      %%%%%%%%%%%%%%%%%%%
                      %%%backward sweep%%
                      %%%%%%%%%%%%%%%%%%%
for i=1:Nbus
    if terminbus(i,1)~=0

    nodenumber=terminbus(i,1);
    
    for j=1:Nbus
        if impdance(j,nodenumber)~=0
            connectedfrom=j;
        
        end
    end
    courant(connectedfrom,nodenumber)=conj(complex(BD(nodenumber,2),BD(nodenumber,3))/Busvoltage(nodenumber,1));
    end
end

for i=Nbus:-1:1%for intermediate nodes
    if interbus(i,1)~=0
            for k=1:Nbus %find which upper node is connected to current intermediate node
                if impdance(k,interbus(i,1))~=0 %&& k~=IntermediateNodes(i,1)
                    upperNode=k;
                    break;
                end
            end
        courant(upperNode,interbus(i,1))=(conj(complex(BD(interbus(i,1),2),BD(interbus(i,1),3))/Busvoltage(interbus(i,1),1)));%current equal to load or generation connected to intermediate node
            for i1=1:Nbus%current to intermediate node is equal to currents going from inter node to all connected downstream nodes + current of node itself if it exists
                if impdance(interbus(i,1),i1)~=0 && i1~=interbus(i,1)
                    courant(upperNode,interbus(i,1))=courant(upperNode,interbus(i,1))+courant(interbus(i,1),i1);
                end
            end
          if upperNode==refbus
          end
        

    end

end
                     %%%%%%%%%%%%%%%%%
                     %%forward sweep%%
                     %%%%%%%%%%%%%%%%%
                    
 voltage(:,1)=Busvoltage(:,1);

     for i=1:Nbus
    for j=1:Nbus
        if impdance(i,j)~=0
            Busvoltage(j,1)=Busvoltage(i,1)-courant(i,j)*impdance(i,j);
        end
    end
end
count=0;
p=0;
for p=2:1:Nbus 
if abs(voltage(p)-Busvoltage(p))<=epsilon
 count=count+1; 
end
end
if count==(Nbus-1) 
break
end
iteration=iteration+1;
end

                               %%%%%% l'écoulement de puissance%%%%
 bilan=zeros(Nbr,3);
 c=0;
 for i=1:Nbus
    
     for j=1:Nbus
         if courant(i,j)~=0
             c=c+1;
              bilan(c,1)=i;%depart 
              bilan(c,2)=j;%arrive
    bilan(c,3)=Busvoltage(i)*conj(courant(i,j))+Busvoltage(j)*conj(-1*courant(i,j)); %les perts par la methode Sloss=Sij+Sji
         end
     end
 end
sLOSS_cap=sum(bilan(:,3));
Qloss=imag(sLOSS_cap);
Ploss=real(sLOSS_cap);
modul=abs(Busvoltage);

end

Actaully without constraints fval minimized but when I add constraints fval increases, I want to minimize it less than 0.0203

can change some thing to minimize it

[Torsten]

Do you know a feasible initial vector for "capplcsz_opt" to start with ?

Answer 1

[Walter Roberson]

You have

constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmax, Vmin);

when you should have

constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmin, Vmax);

[Khlifi Wahid]

ok,can I make integer variables of these indices distinct from each other?

intcon = [1, 3, 5, 7];

[capplcsz_opt, fval, exitflag] = ga(@power_flow_cap, 2*ncap, [], [], [], [], lb, ub, constraints, intcon, options);

[Torsten]

ok,can I make integer variables of these indices distinct from each other?

Could you explain what you mean ?

[Khlifi Wahid]

ok, I don't want the same number in capplcsz_opt for exemple [4 0.112 4 0.58 5 1.258 9 0.3257 ] number four is repeated

I want to make integer variables of these indices [1 3 5 7 ] distinct

[Torsten]

Add the constraints

(capplcsz_opt(1)-capplcsz_opt(3))^2 >= 1
(capplcsz_opt(1)-capplcsz_opt(5))^2 >= 1
(capplcsz_opt(1)-capplcsz_opt(7))^2 >= 1
(capplcsz_opt(3)-capplcsz_opt(5))^2 >= 1
(capplcsz_opt(3)-capplcsz_opt(7))^2 >= 1
(capplcsz_opt(5)-capplcsz_opt(7))^2 >= 1

in function "nonlcon".

[Khlifi Wahid]

That works,thanks a lot