I can't minimize fval value using this code

3 views (last 30 days)
Khlifi Wahid
Khlifi Wahid on 12 Mar 2024
Commented: Khlifi Wahid on 17 Mar 2024
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
  1 Comment
Torsten
Torsten on 12 Mar 2024
Edited: Torsten on 12 Mar 2024
Do you know a feasible initial vector for "capplcsz_opt" to start with ?

Sign in to comment.

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 12 Mar 2024
You have
constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmax, Vmin);
when you should have
constraints = @(capplcsz_opt) constraint_GA(capplcsz_opt, Vmin, Vmax);
  5 Comments
Show 3 older comments
Torsten
Torsten on 16 Mar 2024
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".

Sign in to comment.

Categories

Find more on Whos in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!