GA does not solve problems with integer and equality constraints
Show older comments
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)
Accepted Answer
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".
7 Comments
Khlifi Wahid
on 4 Mar 2024
Khlifi Wahid
on 4 Mar 2024
Torsten
on 4 Mar 2024
Does the result make sense at all ? Is it an improvement against the initial parameters ?
Khlifi Wahid
on 4 Mar 2024
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
on 6 Mar 2024
More Answers (2)
John D'Errico
on 3 Mar 2024
Edited: John D'Errico
on 3 Mar 2024
0 votes
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).
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.
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.
Walter Roberson
on 3 Mar 2024
0 votes
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.
5 Comments
John D'Errico
on 3 Mar 2024
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.
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)
capplcsz_opt([1 3 5 7])
[ Ploss,Busvoltage]=power_flow_cap(capplcsz_opt)
[g,h]=constraint_GA(capplcsz_opt, Vmin, Vmax)
idx = find(g>0)
g(idx)
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');
axis([1 33 0 1.5]);
figure(2);
bar(abs(v),'b');
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
on 4 Mar 2024
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
on 4 Mar 2024
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
on 6 Mar 2024
Categories
Find more on Quadratic Programming and Cone Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
