GA does not solve problems with integer and equality constraints

When I test this code which is shown below. I got an error of " GA does not solve problems with integer and equality constraints."

what should I change to solve this problem.

[ploss,v]=power_flow1()
ncap=4;
min=0.1;
max=1;
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', 500,'StallGenLimit',100,'TimeLimit', 500,'StallTimeLimit', 50,'PlotFcn',@gaplotbestf);
constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmax, Vmin);
[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)

Answer 1

[Torsten]

Your code from above works for me. Test it.

Because you only have inequality constraints, you must set h = [] instead of h = zeros(60,1) in function "constraint_ga".

[Khlifi Wahid]

I found a value of -4 for the exitflag. do you know what should i change to solve this problem

[Torsten]

The code you posted gives exitflag = -2 (No feasible point found) (see above). A remedy would be to start with a feasible point.

I don't know which actual code with exitflag = - 4 you are talking about.

[Khlifi Wahid]

can i change something to make it 1

[Torsten]

Does the result make sense at all ? Is it an improvement against the initial parameters ?

[Khlifi Wahid]

Unfortunately, the values are not within the maximum and minimum value, and there is no improvement in fval.

[Torsten]

You want to have V=abs(Busvoltage) within Vmin and Vmax, and this is the case (see the code below).

[g,h]=constraint_GA(capplcsz_opt, Vmin, Vmax)
idx = find(g>0)

gives idx = [].

I don't know how you see that there is no improvement for fval.

[Khlifi Wahid]

Actually, I want to minimize the fval value; it must be less than 0.0203.

Answer 2

[John D'Errico]

The coupling of equality constraints with integer constraints makes a significant problem. For example, suppose I wanted to solve a problem where x and y were constrained to be integer, together with a rather simple equality constraint? Thus just require that

x^2 - 13*y^2 = 1

This is in fact known as a Pell equation. When you multiply y^2 there by 13, there are not too many integer solutions. The first integer solution does not arise until arise until you get to (x,y) == (649,180). And if we change the constraint to

x^2 - 109*y^2 = 1

the smallest positive integer solution apparently lives at (158070671986249,15140424455100).

https://en.wikipedia.org/wiki/Pell%27s_equation

Worse, if the factor applied to y is the square of an integer, then no solutions will ever exist. Should GA be able to know these facts about Pell equations?

So why am I pointing out Pell equations? Because they are a simple class of equation form where the solutions are not trivial to identify, AND the solutions often lie few and far between. Should GA know the necessary number theory to be able to recognize a Pell equation for what it is, then to solve a Pell equation as part of the constraints? OF COURSE NOT!

Even simpler, suppose the constraints were purely linear? Even there, things get tricky. Solving a system of linear Diophantine equations just to identify potential feasible solutions again requires an understanding of number theory, though much simpler in that case.

Maybe you think that GA should be able to look at sets of infeasible parameter values, and from that, know how to find a feasible solution? NOPE. Again, look at the Pell equation example. Given the constraint

x^2 - 109*y^2 = 1

and some set of sample values for x and y that fail to satisfy the constraint, can you simply identify a pair that does satisfy the constraint? Again, unless you understand Pell equations well enough to know how to solve for a general solution to that equation, nothing will easily help you.

I can give a huge variety of examples. Um, how about this simple constraint on the integer x:

isprime(x*2^x + 1) == true

that is, require that x is an integer, AND that x*2^x+1 is a prime number. This family of numbers are known as the Cullen numbers.

https://en.wikipedia.org/wiki/Cullen_number

They are rather rarely prime. Sequence A005849 in the OEIS (Online Encyclopedia of Integer Sequences) lists the complete known set of values for x, such that a prime number arises is:

{1, 141, 4713, 5795, 6611, 18496, 32292, 32469, 59656, 90825, 262419, 361275, 481899, 1354828, 6328548, 6679881}

The point being, GA would need to be an expert in number theory to handle these problems. It is not so. The combination of equality constraints with integer variables simply makes the problem too complex.

So, to try to answer your actual question as to what should you change to resolve it, you could just try leaving the variables to be non-integer, only rounding them inside the objective function. GA might survive the discontinuities introduced in the objective function. Or, you might just test all possible sets of integer variables, taking the best option. That would completely avoid any optimizer.

Answer 3

[Walter Roberson]

What you have to do is not declare [1 3 5 7] as being integer valued, and instead declare your own MutationFcn and CrossoverFcn and InitialPopulation, with those functions being defined in such a way that they "just happen" to always make [1 3 5 7] integers

ga() cannot handle the combination of integer valued variables together with nonlinear constraints, so you have to declare all the variables to be continuous but arrange the MutationFcn and CrossOverFcn to enforce integer values for the appropriate variables.

[John D'Errico]

As I suggested, this is equivalent to rounding the values inside, but not enforcing the variables to be integer. However you do it though, this will leave you with a piecewise constant objective function. It will probably result in poor convergence, doing a lot more work than it really needs, or it may get stuck.

[Khlifi Wahid]

actually, I want to place 4 batteries in 4 nodes out of 33 nodes. The location is identified based on the lowest node voltages (Busvoltag). Each node has a voltage value calculated by the power_flow_cap function. We have two variables for each battery (so 8 variabls): its dimension and its location. The location ranges from 1 to 33, and the dimension ranges from 0.1 to 1. Additionally, the location must be an integer. After placement, the voltage must be within the range of 0.95 to 1.05. this objective function is to minimize the losses (ploss) and Busvoltage is the voltage of 33 nodes

You can find attached the code files.

global LDD BDD
LDD=[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];

BDD=[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];

%clear all
%clc
%[ploss,v]=power_flow1()
ncap=4;
min=0.1;
max=1;
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', 500,'StallGenLimit',100,'TimeLimit', 500,'StallTimeLimit', 50,'PlotFcn',@gaplotbestf);
constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmax, Vmin);
[capplcsz_opt,fval,exitflag]=ga(@power_flow_cap,2*ncap,[],[],[],[],lb,ub,constraints,[1 3 5 7],options)

ga stopped because the average change in the penalty function value is less than options.FunctionTolerance but 
constraints are not satisfied to within options.ConstraintTolerance.

capplcsz_opt = 1×8

15.0000 1.0000 17.0000 1.0000 14.0000 1.0000 33.0000 1.0000

fval = 0.4831

exitflag = -2

capplcsz_opt([1 3 5 7])

ans = 1×4

15 17 14 33

[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)

Ploss = 0.4831

Busvoltage = 
   1.0000 + 0.0000i
   0.9980 - 0.0020i
   0.9889 - 0.0127i
   0.9852 - 0.0205i
   0.9816 - 0.0288i
   0.9784 - 0.0511i
   0.9853 - 0.0595i
   0.9830 - 0.0719i
   0.9872 - 0.0928i
   0.9918 - 0.1139i
   0.9916 - 0.1174i
   0.9915 - 0.1241i
   1.0020 - 0.1544i
   1.0108 - 0.1670i
   1.0145 - 0.1755i
   1.0156 - 0.1806i
   1.0224 - 0.1908i
   1.0219 - 0.1908i
   0.9974 - 0.0022i
   0.9939 - 0.0033i
   0.9932 - 0.0037i
   0.9925 - 0.0040i
   0.9854 - 0.0131i
   0.9787 - 0.0145i
   0.9754 - 0.0152i
   0.9772 - 0.0517i
   0.9756 - 0.0526i
   0.9703 - 0.0577i
   0.9667 - 0.0614i
   0.9648 - 0.0630i
   0.9664 - 0.0707i
   0.9677 - 0.0731i
   0.9707 - 0.0756i

[g,h]=constraint_GA(capplcsz_opt, Vmin, Vmax)

g = 60×1

-0.0500 -0.0480 -0.0390 -0.0354 -0.0321 -0.0298 -0.0371 -0.0356 -0.0416 -0.0484

h =

     []

idx = find(g>0)

idx =

  0×1 empty double column vector

g(idx)

ans =

  0×1 empty double column vector

a = [];
b = [];
for i = 1:length(capplcsz_opt)
    if mod(i, 2) == 0
        b(end + 1) = capplcsz_opt(i);
    else
        a(end + 1) = capplcsz_opt(i);
    end
end
%constraints = @(capsz_opt) constraint_GA(capsz_opt, Vmax, Vmin);

%[capplcsz_opt,fval]=ga(@power_flow_cap1,2*ncap,[],[],[],[],lb,ub,constraints,options);

%[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)

figure(1);
%hold on

bar(Busvoltage,'r');

Warning: Using only the real component of complex data.

axis([1 33 0 1.5]);

figure(2);

bar(abs(v),'b');

Unrecognized function or variable 'v'.

axis([1 33 0 1.5]);
grid on
%disp('entrer le nombre des condensateurs');
%ncap=input('entrer le nombre des c\n n= ');
%lb=zeros(2*ncap,1);
%ub=zeros(2*ncap,1);
%min = input('entrer la val min des c\n Cmin= ')
%max=input('entrer la val max des c\n Cmax= ')

function [ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)
global LDD BDD
%global ldata_o ncap capplc;
%display('entrer la puissance de base du system en MVA');
Sbase=1;
%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;
LD = LDD;
BD = BDD;
%BD=load('bdata33bus.m');
%LD=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);
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);
Ploss=real(sLOSS);
modul=abs(Busvoltage);

end

function [g,h]=constraint_GA(capplcsz_opt, Vmin, Vmax)
g=zeros(60,1);
h=[];
[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt);
V=abs(Busvoltage);
g(1)=V(1)-Vmax ;
 g(12)=   V(2)-Vmax;
g(3) = V(3)-Vmax;
 g(4) = V(4)-Vmax;
  g(5) =V(5)-Vmax;
  g(6) =V(6)-Vmax;
 g(7)=  V(7)-Vmax;
 g(8) = V(8)-Vmax;
 g(9)  =V(9)-Vmax;
  g(10) =V(10)-Vmax;
  g(11) =V(11)-Vmax;
  g(12) =V(12)-Vmax;
  g(13) =V(13)-Vmax;
  g(14) =V(14)-Vmax;
  g(15) =V(15)-Vmax;
  g(16) =V(16)-Vmax;
  g(17) =V(17)-Vmax;
   g(18)=V(18)-Vmax;
 g(19)  =V(19)-Vmax;
  g(20) =V(20)-Vmax;
 g(21)  =V(21)-Vmax;
 g(22)  =V(22)-Vmax;
 g(23)  =V(23)-Vmax;
g(24)   =V(24)-Vmax;
 g(25)  =V(25)-Vmax;
   g(26)=V(26)-Vmax;
   g(27)=V(27)-Vmax;
  g(28) =V(28)-Vmax;
  g(29) =V(29)-Vmax;
  g(30)= V(30)-Vmax;
  g(31) =V(31)-Vmax;
  g(32) =V(32)-Vmax;
  g(33) =V(33)-Vmax;
 g(1) =Vmin-V(1);
 g(2) =Vmin- V(2);
 g(3)=Vmin-V(3);
 g(4)=Vmin-V(4);
 g(5)=Vmin-V(5);
 g(6)=Vmin-V(6);
 g(7)=Vmin-V(7);
 g(8)=Vmin-V(8);
 g(9)=Vmin-V(9);
 g(10)=Vmin-V(10);
g(11)=Vmin-V(11);
g(12)=Vmin-V(12);
g(13)=Vmin-V(13);
g(14)=Vmin-V(14);
g(15)=Vmin-V(15);
g(16)=Vmin-V(16);
g(17)=Vmin-V(17);
g(18)=Vmin-V(18);
g(19)=Vmin-V(19);
g(20)=Vmin-V(20);
g(21)=Vmin-V(21);
g(22)=Vmin-V(22);
g(23)=Vmin-V(23);
g(24)=Vmin-V(24);
g(25)=Vmin-V(25);
g(26)=Vmin-V(26);
g(27)=Vmin-V(27);
g(28)=Vmin-V(28);
g(29)=Vmin-V(29);
g(30)=Vmin-V(30);
g(31)=Vmin-V(31);
g(32)=Vmin-V(32);
g(33)=Vmin-V(33);
end

[Walter Roberson]

You still need to write your own MutationFcn and CrossoverFcn that accept double precision inputs and "just happen" to give integer results for particular variables. And you need the InitialPopulation large enough to completely fill one generation.

[Walter Roberson]

g(1)=V(1)-Vmax ;
 g(12)=   V(2)-Vmax;
g(3) = V(3)-Vmax;

You output to g(12) there instead of to g(2). You later output to g(12) without having outout to g(2)

 g(1) =Vmin-V(1);

You overwrite what you have written to g(1) and the rest.

You could simply code

 g = [V(:) - Vmax; Vmin - V(:)];

[Khlifi Wahid]

Actually, the problem of limiting the values ​​between vmax and vmin has been solved but I did not find an optimal faval value and exitflag=-5, I want to minimize it less than 0.0203 can i change something to do it