%本程序為電動汽車充電的模擬退火PSO算法。并且加入了高斯變異++考慮用戶側滿意度、電動汽車建模的概率性
%為了使節點按先PQPV最后平衡節點次序編號，以便與公式對照，節點1與節點10對調。
clear;
num=xlsread（‘網絡參數.xls‘）;
e=xlsread（‘節點功率參數.xls‘）;
num（:6）=num（:4）+num（:5）*sqrt（-1）;
Y=zeros（10）;???%生成空導納矩陣
%創建節點導納矩陣
for?i=1:9;???%不考慮變壓器和線路對地導納
????p=num（i2）;q=num（i3）;
????Y（pq）=Y（pq）-1/num（i6）;
????Y（qp）=Y（pq）;
????Y（pp）=Y（pp）+1/num（i6）;
????Y（qq）=Y（qq）+1/num（i6）;
end????%根據支路參數算導納矩陣
temp=Y（102）;%保持原來的數據
P=zeros（948）;?%節點基本負荷，9個PQ節點一天48個點的負荷
Q=zeros（948）;
Pv2r=zeros（948）;?%汽車用來行駛的電功率
Pev_leave=zeros（948）;%電動汽車離開電網的時間
%Pv2r1=zeros（948）;
Xc=zeros（1148）;?%Xc（1:D:）為一天48個點的控制變量，Xc（D+1:）為一天48個點的目標函數值
%獲取節點日負荷曲線數據
for?i=1:9?%平衡節點10不計算，標幺值
???for?j=1:48
????????P（ij）=e（26j）*e（i+44）;?Q（ij）=e（26j）*e（i+45）;
????????%Pv2r（ij）=e（18j）*e（i+49）/1000;?%汽車行駛功率，標幺值
????????Pv2r（ij）=e（20j）*e（i+49）/1000;
????????Pev_leave（ij）=e（27j）;%改變了汽車離開電網的分布。
????????
???end
end
for?i=1:48
????[UAngleline_Ploss_P_1（i）P0_10（i）]=niulafa（-P（:i）-Q（:i）1Ytemp）;
end
meanP=mean（P0_10）;
%粒子群優化算法
N=50;?%種群大小
D=48;?%搜索空間維數，即控制變量個數
Tmax=1200;?%最大迭代次數
c1=1.49445;c2=1.49445;?%學習因子
wmax=0.9;wmin=0.4;?%慣性權重
T=1e6;b=0.9;?%模擬退火的溫度控制參數
scope=[00.0312??%Pev范圍，均為標幺值
???????0.951.05?%分接頭t范圍
???????0.931.07?%節點電壓范圍
???????0.0240.24
???????-0.03120.0312
???????02];?%電池儲存能量范圍
???%初始化粒子群
?xv_swarm=zeros（112*DN）;?%生成粒子群位置、速度、當前最優值矩陣
?%xv_swarm（1:101:DN）為N個粒子的位置，xv_swarm（1:10D+1:2*DN）為N個粒子的速度，xv_swarm（111:DN）為N個粒
?%子各個時刻的適應值（目標函數值），xv_swarm（即前面48個值的和1149N）為第N個粒子總的適應值，。?
?pg_swarm=zeros（11D+1N+1）;?%生成個體最優位置、群體最優位置和相應最優值矩陣
?%pg_swarm（1:101:D，N）為N個粒子的各自歷史最優位置，pg_swarm（111:DN）為N個粒子的各自歷史每時刻的最小適應值，
?%pg_swarm（1149N）為N個粒子的各自歷史總的最小適應值，
?%pg_swarm（1:101:DN+1）為N個粒子的群體最優位置，pg_swarm（111:DN+1）為N個粒子的群體各個時刻最小適應值，
?%pg_swarm（1149N+1）為N個粒子的群體總的最小適應值
?ev_soc_re=zeros（949）;?%汽車SOC第一列為初始值，ev_soc（：2:49）為9個節點一天48個點的soc值
ev_soc=zeros（949N）;
%ev_soc_re（:1）=rand（91）;
%ev_soc_re（:1）=ev_soc_re（:1）*（scope（42）-scope（41））+scope（41）;?%soc初始化
ev_soc_re（:1）=rand*（0.12-0.05）+0.05;?%soc初始化0.2-0.5
for?n=1:N
ev_soc（:1n）=ev_soc_re（:1）;
end
%ev_soc（1:91）=0.72;
%Pv2r=zeros（948）;
%對N個粒子初始化
for?n=1:N?
????sum_P=0;
?for?j=1:48
?????U=2;?%使進入while循環
?????count1=0;
??while??max（U）>1.07||min（U）<0.93||max（max（abs（line_P）））>3.932||max（ev_soc（:j+1n））>scope（42）||min（ev_soc（:j+1n））