Optimisation PSO pour systèmes énergétiques
Optimisation PSO pour systèmes énergétiques
fermertout
clc
%% Définition du problème
%hade bala
C=0,9;%prix de l'électricité
W=0,2 ; %probabilité de perte de charge%
K=0,99;%facteur d'énergie renouvelable%
nvars=1;%un seul système
LB=[0 0 1 0];% Limite inférieure du problème
UB=[45 8 20 10];% Limite supérieure du problème
%%
20;%100;% Nombre maximum d'itérations
NPOP=5;%20;% Nombre de la population
% Déterminer le pas maximum pour la vitesse
1:4
siLB(d)>-1e20 && UB(d)<1e20
(UB(d)-LB(d))/NPOP
sinon
inf
fin
fin
%% Paramètres initiaux du PSO
% w=0.5; % Poids d'inertie
% wdamp=0.99; % Rapport d'amortissement du poids d'inertie
% c1=2; % Coefficient d'apprentissage personnel
% c2=2; % Coefficient d'apprentissage global
2.05;
2.05;
phi=phi1+phi2;
chi=2/(phi-2+sqrt(phi^2-4*phi));
chi; % Poids d'inertie
1
c1=chi*phi1; % Coefficient d'apprentissage personnel
c2=chi*phi2; % Coefficient d'apprentissage global
%% Initialisation
tic
empty_particle.position=[];
empty_particle.velocity=[];
empty_particle.cost=[];
empty_particle.[Link]=[];
empty_particle.[Link]=[];
%---------
repmat(empty_particle, NPOP, 1);
inf
[Link]=[];
inf
pour i=1:NPOP
cc=1;%une valeur pour le coût
ww=0.3;% une valeur pour la probabilité de perte de charge%
kkk=2;%renewable energy factor%
0
%C=0.9;%prix de l'électricité
%W=0.2;%probabilité de perte de charge%
%K=0.99;%facteur d'énergie renouvelable%
whileww>=W | kkk>=K% LOLP/REfactor
particle(i).position(1,:)=unifrnd(0,45,1,nvars);%pv kW
particule(i).position(2,:)=unifrnd(0,8,1,nvars; % jours d'autonomie
particle(i).position(3,:)=unifrnd(1,20,1,nvars);%nombre de maisons
particle(i).position(4,:)=unifrnd(0,10,1,nvars);%nombre d'éoliennes
1:4
particle(i).velocity(g,:)=rand(1,nvars);
fin
%----convert------------
p_npv=particule(i).position(1);
ad=particule(i).position(2);
maisons=arrondi(particule(i).position(3));
nwt=arrondi(particule(i).position(4));
%-----------------------
LPSP
ff=ff+1;
ww(i)=LPSP;
kkk(i)=facteur_renouvelable;
fin
prix_electricité
particle(i).[Link]=particle(i).position;
particle(i).[Link]=particle(i).cost;
si particule(i).[Link]ût < [Link]ût
globalbest=particule(i).meilleur;
fin
fin
Fminn= zéro(max_it,1);
%% Boucle principale PSO
% disp('Itération Fiabilité');
% disp('-----------------------------');
foru=1:max_it
vv=0;
pour i=1:NPOP
cc=1;% une valeur pour le coût
LPSP=0,3;% une valeur pour la probabilité de perte de charge%
renewable_factor=2;%
bb=0;
%C=0.9;%prix de l'électricité
%W=0.2;%perte de probabilité de charge%
%K=0.99;%facteur d'énergie renouvelable%
% ww=0.3*100;% une valeur pour la probabilité de perte de charge%
% kkk=2*100;%facteur d'énergie renouvelable%
tandis que (LPSP) >= W | (facteur_renouvelable) >= K
fory=1:4
particle(i).velocity(y,:)=w1*particle(i).velocity(y,:)+c1*rand*...
(particule(i).[Link](y,:)-particule(i).position(y,:))...
+c2*rand*([Link](y,:)-particle(i).position(y,:));
%particule(i).vitesse(y,:)=min(max(particule(i).vitesse(y,:),-velmax),velmax);
particle(i).position(y,:)=particle(i).position(y,:)+particle(i).velocity(y,:);
% drapeau=(particule(i).position(kk,:)<LB(kk) | particule(i).position(kk,:)>UB(kk));
% particule(i).vitesse(drapeau)=-particule(i).vitesse(drapeau);
particle(i).position(y,:)=min(max(particle(i).position(y,:),LB(y)),UB(y));
fin
%p_npv(i,:)=arrondi(particle(i).position(1,:));
oo=0;
%ad(i,:)=arrondi(particule(i).position(2,:));
p_npv=arrondi(particule(i).position(1));
ad=arrondi(particule(i).position(2));
maisons=arrondi(particule(i).position(3));
nwt=arrondi(particule(i).position(4));
%[cc ww]=contraintes(c,w,nnn(i,:),zz(i,:),zmax,nmax,N);
LPSP
bb=bb+1;
fin
%----convert------------
% az(i,:)=arrondir(particle(i).position(1,:));
% z(i,:)=arrondir(particle(i).position(2,:));
% n(i,:) = arrondi(particle(i).position(3,:));
% sol(i).pos(1,:)=arrondir(particle(i).position(1,:));
sol(i).pos(2,:)=arrondi(particule(i).position(2,:));
% sol(i).pos(3,:)=arrondir(particle(i).position(3,:));
%-----------------------
%[LPSP
prix_electricité
facteur_renouvelable
objectif_fonction(sol(i).pos);
vv=vv+1;
si particule(i).coût < particule(i).[Link]ût
particle(i).[Link]=particle(i).cost;
particle(i).[Link]=particle(i).position;
si la meilleure coût de la particule(i) < coût_global_meilleur %& rnwfct < rnwfct_meilleur
globalbest=particule(i).meilleur;
% rnwfct_best=rnwfct;
fin
fin
fin
Fminn(u)=[Link];
Xmin=[Link];
p_npv=arrondi([Link](1));
ad=arrondi([Link](2));
maisons=round([Link](3));
nwt=arrondi([Link](4));
% w=wdamp*w;
% aff_disp(['Itération ',num2str(u),': Meilleur coût= ',num2str(Fminn(u))]);
fin
temps=toc;
%% Résultats
fermerTout;
% figure;
% tracer(Fminn,'Largeur de ligne',3);
% ylabel('prix de l'électricité');
% xlabel('Nombre d’itérations');
% ylim([0 0.4])
% xlim([1 max_it ])
Fmin=min(Fminn);
Xmin
LPSP
formatlong
result=[p_npv,ad,houses,nwt,LPSP,price_electricity,renewable_factor,b(1),b(2),b(3),b(4),b(5)];
% barre(resultat)
% légende('pv(kW)','jours d'autonomie','nombre de maisons','nombre de vent
turbines
électricité($/kW)
figure,subplot(2,1,1),x=1:168;bar(x,ali(:,1:4),'groupe');colormap
jet,subplot(2,1,2),x=1:168;bar(x,ali(:,5:6),'groupe');colormapjet
p_npv
annonce
maisons
nwt
LPSP
prix_de_l'électricité
facteur renouvelable
fonction
LPSP
nPng)
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% %% 3) charger les entrées RAFSANJAN/KERMAN
% charger('[Link]');
% charger('[Link]');
% vitesse_du_vent=rafsanjan(:,1);[vitesse_du_vent]=min10tohourly(vitesse_du_vent);%horaire
température=rafsanjan(:,2);[température]=min10àhoraire(température);%horaire
% radiation_solaire=rafsanjan(:,3)/1000;[radiation_solaire]=min10tohourly(radiation_solaire);%horaire%kw
% dégager rafsanjan
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% %% 2) charger les entrées KHASH/SISTAN BALOCHESTAN
% charger('sistan_khash.mat');
% charger('[Link]');
% vitesse_du_vent=sistan_khash(:,1);[vitesse_du_vent]=min10tohourly(vitesse_du_vent);%horaire
% température=sistan_khash(:,2);[température]=min10tohourly(température);%horaire
solar_radiation=sistan_khash(:,3)/1000;[solar_radiation]=min10tohourly(solar_radiation);%horaire%kw
% effacer sistan_khash
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%% 1)charger les entrées NAHAVAND
charger('[Link]');nahavand = décaler(nahavand,60);
charger('[Link]');
nahavand(:,1);[wind_speed]=min10tohourly(wind_speed);%horaire
température=nahavand(:,2);[température]=min10tohourly(température);%horaire
nahavand(:,3)/1000;[solar_radiation]=min10tohourly(solar_radiation);%horaire%kw
clearnahavand
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%% énergie solaire
%#######entrées##############entrées####################entrées##
%température ambiante-hararate mohit
température
7,3
g=rayonnement_solaire;%rayonnement_solaire_horaire_nahavand';%kW
%############################################################
gref=1 ;%1000kW/m^2
25;%température à la condition de référence
kt=-3.7e-3;% coefficient de température de la puissance maximale (1/c0)
tc=tamb2+(0.0256).*g;
upv=0.986;%efficacité du pv avec un angle incliné>>98,6%
p_pvout_hourly=upv*(p_npv.*(g/gref)).*(1+kt.*(tc-tref));%puissance de sortie(kw)_horaire
cleartc
cleartamb2
clairg
%% batterie
%#######entrées##############entrées####################entrées##
%charger la demande >> profil de charge typique horaire d'un foyer rural en kW
load1=[1.5 1 0.5 0.5 1 2 2 2.5 2.5 3 3 5 4.4 4 3.43 3 1.91 2.48 3 3 3.42 3.44 2.51 2];
charge1=charge1/5; charge1=charge1*2;% le maximum serait de 2kW, la moyenne est de 1kW
%la courbe de charge d'une journée complète de consommation typique, étude de cas palestinienne
%load1=[1 1 1 1 1.5 2 2.5 2.5 1.5 1.5 1.5 2 2 2.5 2.5 3 3.5 4 4 3.5 2.5 1.5 1 1];
%houses=1;%number of houses in a village
load2=maisons.*load1;%charge totale en une journée pour tout le village
%données de charge horaire pour une année
a=0;
pour i = 1:1:360
a=[a,charger2];
fin
a(1)= [];
Pl1=a;
%ad=3;%autonomie quotidienne
0,92;
ub=0.85;
0,8 ; % de décharge 0,5 0,7 dans l'article 80%
el=mean(load2);
capacité de batterie 40 kWh
% ############################################################
cwh=(el*ad)/(uinv*ub*dod);%capacité de stockage pour batterie,bmax,kW
Éolienne
%#######entrées##############entrées####################entrées##
4;
rw=cell2mat(WindTurbines(MCH,2));% rw=4;%diamètre des pales(m)
aw=cell2mat(WindTurbines(MCH,3));% aw=pi*(rw)^2;% Zone Balayée >> pi x Rayon = Zone Balayée par les Lames
uw=cell2mat(WindTurbines(MCH,4));% uw=0.95;%
vco=cell2mat(WindTurbines(MCH,5));% vco=25;% arrêt
vci=cell2mat(WindTurbines(MCH,6));% vci=3;% seuil d'entrée
vr=cell2mat(Eoliennes(MCH,7));% vr=8;%vitesse nominale(m/s)
pr=cell2mat(HélicesÉoliennes(MCH,8));% pr=2;% puissance nominale(kW)
pmax=cell2mat(WindTurbines(MCH,10));% % pmax=2.5;% puissance de sortie maximale (kW)
pfurl=cell2mat(WindTurbines(MCH,9));% pfurl=2.5;%puissance de sortie à la vitesse de coupure (9kW)
% ############################################################
vitesse_du_vent;
groupe électrogène diesel
Png=4;Png=nPng*Png;%kW puissance de sortie du générateur diesel
Bg=0.08145;%1/kW
Ag=0,246;%1/kW
Pg=4;Pg=nPng*Pg;%puissance nominale kW
consommation de carburant du générateur diesel
Fg=Bg*Pg+Ag*Png;
%% PROGRAMME PRINCIPAL
contribution=zeros(5,8640);% pv, éolien, batterie, contribution diesel à chaque heure
Ebmax=cwh;%40kWh%capacité de la batterie 40 kWh
Ebmin=cwh*(1-dod);%40kWh
SOCb=0,2 ; %état de charge de la batterie >> 20%
Eb=zeros(1,8640);
time1=zeros(1,8640);
zeros(1,8640);
Edump=zeros(1,8640);
Edch=zeros(1,8640);
Ech=zeros(1,8640);
Eb(1,1)=SOCb*Ebmax;%état de charge pour le temps de démarrage
%^^^^^^^^^^^^^^DÉBUT^^^^^^^^^^^^^^^^^^^^^^^^
Pl=Pl1;
clearPl1;
%^^^^^^^^^^Calcul de la puissance de sortie^^^^^^^^
%calcul de l'énergie solaire
Pp=p_pvout_hourly;%puissance de sortie(kw)_horaire
pouri=1:1:8640
siPp(i)>p_npv
Pp(i)=p_npv;%si la puissance de sortie du pv dépasse le maximum
fin
fin
% calcul de l'énergie éolienne
1:1:8640
%pr *((v2(t)-vci)/(vr-vci))^3pr+(((pfurl-pr)/(vco-vr))*(v2(t)-vr));
si v2(t) < vci % v2 >> vitesse_du_vent_horaire ;
0;
sinon si vci <= v2(t) && v2(t) <= vr
pwtg(t)=(pr/(vr^3-vci^3))*(v2(t))^3-(vci^3/(vr^3-vci^3))*(pr);
elseif vr <= v2(t) && v2(t) <= vco
pwtg(t)=pr;
sinon
pwtg(t)=0;
fin
Pw(t)=pwtg(t)*uw*nwt;%puissance électrique d'une éolienne
fin
2:1:8640
%^^^^^^^^^^^^^^LIRE LES ENTRÉES^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^COMPARAISON^^^^^^^^^^^^^^^^^^^
siPw(t)+Pp(t)>=(Pl(t)/uinv)
%^^^^^^EXÉCUTER LE CHARGEMENT AVEC ÉOLIENNE ET PV^^^^^^
siPw(t)+Pp(t)>Pl(t)
%^^^^^^^^^^^^^^CHARGE^^^^^^^^^^^^^^^^^^^^^^^^^^
[Edump,Eb,Ech] = charge(Pw,Pp,Eb,Ebmax,uinv,Pl,t,Edump,Ech);
time1(t)=1;
Pp(t)
contribution(5,t)=Edump(t);
sinon
Eb(t)=Eb(t-1);
retourner
fin
sinon
Décharge
[Eb, Edump, Edch, diesel, time1, t] =
décharge(Pw,Pp,Eb,Ebmax,uinv,Pl,t,Pg,Ebmin,Edump,Edch,Ech,diesel,time1);
Pp(t)
contribution(5,t)=Pl(t);
fin
fin
%% traçage
% figure
a=contribution';
b=somme(a);
b(4)/(b(1)+b(2));
% h=diagramme(b);
% colormap jet;
% légende('PV','ÉOLIEN','BATTERIE','DIESEL');
%fiabilité
% de probabilité de perte de charge = somme(chargement - pv - vent + batterie) / somme(chargement)
total_loss=0;
1:1:8640
% aa(t)=Pl(t)-Pp(t)-Pw(t)+Eb(t);
siPl(t)>(Pp(t)+Pw(t)+(Eb(t)-Ebmin)+diesel(t))
total_loss=total_loss+(Pl(t)-(Pp(t)+Pw(t)+(Eb(t)-Ebmin)+diesel(t)));
fin
fin
LPSP=perte_totale/(somme(Pl));
fiabilité = somme(aa) / somme(Pl);
économique(diesel, Pl, Fg, cwh);
économie_rapide_iran(diesel,Pl,Fg,cwh,p_npv,nwt,maisons);
ali=[Pp(1:168)',Pw(1:168)',Eb(1:168)',diesel(1:168)',Pl(1:168)',Edump(1:168)'];
Edump=somme(Edump);
fin