GA does not solve problems with integer and equality constraints

Question

Khlifi Wahid 2024-3-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089681-ga-does-not-solve-problems-with-integer-and-equality-constraints

评论： Khlifi Wahid 2024-3-6

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-3-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089681-ga-does-not-solve-problems-with-integer-and-equality-constraints#answer_1420291

编辑：Torsten 2024-3-4

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 个评论
显示 5更早的评论隐藏 5更早的评论

Torsten 2024-3-4

编辑：Torsten 2024-3-4

在 MATLAB Online 中打开

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 2024-3-6

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

请先登录，再进行评论。

Answer 2

John D'Errico 2024-3-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089681-ga-does-not-solve-problems-with-integer-and-equality-constraints#answer_1420211

编辑：John D'Errico 2024-3-3

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 3

Walter Roberson 2024-3-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089681-ga-does-not-solve-problems-with-integer-and-equality-constraints#answer_1420241

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 个评论
显示 3更早的评论隐藏 3更早的评论

Khlifi Wahid 2024-3-3

编辑：Torsten 2024-3-4

在 MATLAB Online 中打开

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

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

%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

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

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

courant(connectedfrom,nodenumber)=conj(complex(BD(nodenumber,2),BD(nodenumber,3))/Busvoltage(nodenumber,1));

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

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

if upperNode==refbus

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

count=0;

p=0;

for p=2:1:Nbus

if abs(voltage(p)-Busvoltage(p))<=epsilon

count=count+1;

end

if count==(Nbus-1)

break

end

iteration=iteration+1;

end

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

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 2024-3-4

在 MATLAB Online 中打开

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 2024-3-6

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

请先登录，再进行评论。