【有序充电】基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化附Matlab代码
【有序充电】基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化附Matlab代码
TT_Matlab
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，完整matlab代码或者程序定制加qq1575304183。
✅作者简介：热爱科研的Matlab仿真开发者，修心和技 术同步精进，
代码获取、论文复现及科研仿真合作可私信。
个人主页： Matlab科研工作室
个人信条：格物致知。
更多Matlab完整代码及仿真定制内容点击
智能优化算法 神经网络预测 雷达通信 无线传感器 电力系统
信号处理 图像处理 路径规划 元胞自动机 无人机
内容介绍
随着电动汽车的普及，如何高效地进行充电成为了一个亟待解决的问题。在多充电站的情况下，如何实现电动汽车的有序充电，提高充电效率，降低充电成本，成为了一个重要的研究方向。本文将介绍基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化算法流程。
1.问题描述
考虑在IEEE33环境下，有多个充电站和多辆电动汽车。每个充电站都有一定数量的充电桩，每个充电桩的充电速度不同。每辆电动汽车都有一个充电需求，需要在规定时间内完成充电。同时，每个充电站的电力供应有限，需要在保证充电效率的前提下，合理分配电力资源。因此，需要设计一个有序充电方案，使得每辆电动汽车都能在规定时间内完成充电，同时最大化充电效率，最小化充电成本。
2.粒子群算法
粒子群算法是一种群体智能算法，模拟鸟群或鱼群等生物的集体行为，通过模拟群体的协作和竞争，寻找最优解。粒子群算法的基本思想是将每个解看作一个粒子，通过不断迭代，逐步接近最优解。粒子群算法具有收敛速度快、易于实现等特点，适用于解决复杂的优化问题。
3.算法流程
(1)初始化粒子群，包括每辆电动汽车的充电需求、每个充电站的充电桩数量和充电速度等信息。
(2)计算每个充电站的电力供应情况，并根据充电桩的充电速度和电动汽车的充电需求，确定每辆电动汽车的充电站和充电桩。
(3)根据每个充电站的电力供应情况，计算每辆电动汽车的充电时间和充电成本。
(4)根据充电时间和充电成本，更新每个粒子的适应度。
(5)根据粒子的适应度，更新粒子的速度和位置。
(6)重复步骤(2)-(5)，直到达到预设迭代次数或达到最优解。
4.实验结果
通过对IEEE33环境下的多充电站电动汽车群有序充电优化算法进行实验，得到了如下结果：
(1)在充电效率和充电成本方面，粒子群算法优于其他算法。
(2)随着粒子数的增加，算法的收敛速度更快，但计算时间也会增加。
(3)在充电桩数量和充电速度等参数的变化下，算法的性能表现稳定。
5.总结
本文介绍了基于粒子群算法实现IEEE33环境下的多充电站电动汽车群有序充电优化算法流程。该算法具有收敛速度快、易于实现等特点，在充电效率和充电成本方面优于其他算法。未来可以进一步优化算法，提高算法的性能表现。
部分代码
%相关原始数据格式说明如下：
%n——节点个数；n1——支路条数；isb——平衡节点号；H——PQ节点个数（为后面形成PVU存储PV节点初始电压用）；pr——误差精度。
%B1——支路参数矩阵，其中第一列和第二列是起始节点编号和终点节点编号,第三列、第四列、第五列、第六列分别为：支路电阻、电抗、变压器变比、电纳。（不考虑电导）
%B2——节点参数矩阵，其中第一列和第二列为节点编号和节点类型；第三列到第六列分别为：注入有功、注入无功、电压幅值、电压相位。
%节点类型分类如下：“
0
”为平衡节点,“
1
”为PQ，“
2
”为PV节点；“
3
”为PQ（V）节点，“
4
”为PI节点。
function [Ploss ,V_amp]=powerflow(car_1,car_2,car_3,P_flex)
n=
33
;
n=
33
;
n1=
32
;
isb=
1
;
H=
32
;
%%%
%%%
%%%
%%%
%
18
节点加DG PQV处理
pr=
0
.
01
;
v_amp=
0
;
B1=[
1
2
0
.
00
922
0
.
0047
i
1
0
;
2
3
0
.
004
93
0
.
02511
i
1
0
;
3
4
0
.
0366
0
.
01
864i
1
0
;
4
5
0
.
03
811
0
.
01
941i
1
0
;
5
6
0
.0819
0
.
0707
i
1
0
;
6
7
0
.
01
872
0
.
061
88i
1
0
;
7
8
0
.
07114
0
.
02351
i
1
0
;
8
9
0
.
103
0
.
074
i
1
0
;
9
10
0
.
1044
0
.
074
i
1
0
;
10
11
0
.
01
966
0
.
0065
i
1
0
;
11
12
0
.
03744
0
.
0123
8i
1
0
;
12
13
0
.
1468
0
.
1155
i
1
0
;
13
14
0
.
05416
0
.
0712
9i
1
0
;
14
15
0
.
05
91
0
0
.
0526
i
1
0
;
15
16
0
.
07463
0
.
05450
i
1
0
;
16
17
0
.
1289
0
.
1721
i
1
0
;
17
18
0
.
0732
0
.
0574
i
1
0
;
2
19
0
.
0164
0
.
01565
i
1
0
;
19
20
0
.
15042
0
.
13554
i
1
0
;
20
21
0
.
040
95
0
.
047
84i
1
0
;
21
22
0
.
070
89
0
.09373i
1
0
;
3
23
0
.
04512
0
.
030
83i
1
0
;
23
24
0
.0898
0
0
.
070
91i
1
0
;
24
25
0
.0896
0
0
.
07011
i
1
0
;
6
26
0
.
0203
0
.
01034
i
1
0
;
26
27
0
.
02
842
0
.
01447
i
1
0
;
27
28
0
.
1059
0
.09337i
1
0
;
28
29
0
.08042
0
.
07006
i
1
0
;
29
30
0
.
05075
0
.
025
85i
1
0
;
30
31
0
.09744
0
.0963i
1
0
;
31
32
0
.
03105
0
.
0361
9i
1
0
;
32
33
0
.
03410
0
.
05302
i
1
0
];
B2=[
1
0
0
0
1.05
0
;
2
1
-
0
.
01
-
0
.
006
1
0
;
3
1
-
0
.
00
9 -
0
.
004
1
0
;
4
1
-
0
.
012
-
0
.
00
8
1
0
;
5
1
-
0
.
006
-
0
.
003
1
0
;
6
1
-
0
.
006
-
0
.
002
1
0
;
7
1
-
0
.
02
-
0
.
01
1
0
;
8
1
-
0
.
02
-
0
.
01
1
0
;
9
1
-
0
.
006
-
0
.
002
1
0
;
10
1
-
0
.
006
-
0
.
0035
1
0
;
11
1
-
0
.
0045
-
0
.
003
1
0
;
12
1
-
0
.
006
-
0
.
0035
1
0
;
13
1
-
0
.
006
-
0
.
0035
1
0
;
14
1
-
0
.
012
-
0
.
00
8
1
0
;
15
1
-
0
.
006
-
0
.
001
1
0
;
16
1
-
0
.
006
-
0
.
002
1
0
;
17
1
-
0
.
006
-
0
.
002
1
0
;
18
1
-
0
.
00
9 -
0
.
004
1
0
;
19
1
-
0
.
00
9 -
0
.
004
1
0
;
20
1
-
0
.
00
9 -
0
.
004
1
0
;
21
1
-
0
.
00
9 -
0
.
004
1
0
;
22
1
-
0
.
00
9 -
0
.
004
1
0
;
23
1
-
0
.
00
9 -
0
.
005
1
0
;
24
1
-
0
.
042
-
0
.
02
1
0
;
25
1
-
0
.
042
-
0
.
02
1
0
;
26
1
-
0
.
006
-
0
.
0025
1
0
;
27
1
-
0
.
006
-
0
.
0025
1
0
;
28
1
-
0
.
006
-
0
.
002
1
0
;
29
1
-
0
.
012
-
0
.
007
1
0
;
30
1
-
0
.
02
-
0
.
06
1
0
;
31
1
-
0
.
015
-
0
.
007
1
0
;
32
1
-
0
.
021
-
0
.
01
1
0
;
33
1
-
0
.
006
-
0
.
004
1
0
];
for
i=
1
:
33
B2(i,
3
)= B2(i,
3
)*P_flex; %负荷时变系数
B2(i,
4
)= B2(i,
4
)*P_flex;
end
global position
B2(position(
1
),
3
)=B2(position(
1
),
3
)-car_1/
10000
; %电动汽车接入
B2(position(
1
),
4
)=B2(position(
1
),
4
)-
0
.
484
*car_1/
10000
; %功率因素为
0
.
9
B2(position(
2
),
3
)=B2(position(
2
),
3
)-car_2/
10000
; %电动汽车接入
B2(position(
2
),
4
)=B2(position(
2
),
4
)-
0
.
484
*car_2/
10000
; %功率因素为
0
.
9
B2(position(
3
),
3
)=B2(position(
3
),
3
)-car_3/
10000
; %电动汽车接入
B2(position(
3
),
4
)=B2(position(
3
),
4
)-
0
.
484
*car_3/
10000
; %功率因素为
0
.
9
Y=zeros(n); %zeros就是生成一个全
0
的矩阵
Times=
1
; %置迭代次数为初始值
for
i=
1
:n1
p=B1(i,
1
);
q=B1(i,
2
);
Y(p,q)=Y(p,q)-
1
/((B1(i,
3
)+B1(i,
4
))*B1(i,
5
));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+
1
/(B1(i,
3
)+B1(i,
4
))+
0
.
5
*B1(i,
6
);%
Y(q,q)=Y(q,q)+
1
/((B1(i,
3
)+B1(i,
4
))*B1(i,
5
)^
2
)+
0
.
5
*B1(i,
6
);%高压侧阻抗乘以变比平方 输入时注意低压侧在前
end
%disp(
’节点导纳矩阵：’
) ;
Y;
G=real(Y);
B=imag(Y);
OrgS=zeros(
2
*n-
2
,
1
);
DetaS=zeros(
2
*n-
2
,
1
); %将OrgS、DetaS初始化
%创建OrgS，用于存储初始功率参数
Q=
0
;
PQV=
0
;
x=
1.655
;
%%%
x1=
6.7
,x2=
9.85
,x=x1+x2；其中x1为定子漏抗，x2为转子漏抗
xp=
18.8
;
%%%
xc=,xm=,xp=xc*xm/(xc-xm)；其中xc为机端并联电容器电抗，xm为激磁电抗
h=
0
;
for
i=
1
:n
%对PQ(V)节点的处理
h=h+
1
;
if
i~=isb&&B2(i,
2
)==
3
Q(i)=-(B2(i,
5
))^
2
/xp+(-(B2(i,
5
))^
2
+sqrt((B2(i,
5
))^
4
-
4
*(B2(i,
3
))^
2
*x^
2
))/
2
*x;
B2(i,
4
)=Q(i);
B2(i,
2
)=
1
;
PQV=h;
end
end
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%
Ig=
0
.
01
;
Q=
0
;
PI=
0
;
h=
0
;
for
i=
1
:n
%对PI节点的处理
h=h+
1
;
if
i~=isb&&B2(i,
2
)==
4
Q(i)=sqrt(Ig^
2
*((B2(i,
5
))^
2
)-B(i,
3
)^
2
); %e=B2(i,
5
)，f=
0
，e^
2
+f^
2
=B2(i,
5
))^
2
，其中e和f为光伏发电系统接入节点电压的实部和虚部
B2(i,
4
)=Q(i);
B2(i,
2
)=
1
;
PI=h;
end
end
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
h=
0
;
j=
0
;
for
i=
1
:n
%对PQ节点的处理
if
i~=isb&&B2(i,
2
)==
1
h=h+
1
;
for
j=
1
:n
OrgS(
2
*h-
1
,
1
)=OrgS(
2
*h-
1
,
1
)+B2(i,
5
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))+B2(i,
6
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %Pi 书P57页
11
-
45
OrgS(
2
*h,
1
)=OrgS(
2
*h,
1
)+B2(i,
6
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))-B2(i,
5
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %Qi 同上
end
end
end
for
i=
1
:n
%对PV节点的处理
if
i~=isb&&B2(i,
2
)==
2
h=h+
1
;
for
j=
1
:n
OrgS(
2
*h-
1
,
1
)=OrgS(
2
*h-
1
,
1
)+ B2(i,
5
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))+B2(i,
6
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %同上
OrgS(
2
*h,
1
)=OrgS(
2
*h,
1
)+ B2(i,
6
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))-B2(i,
5
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %同上
end
end
end
OrgS;
%创建PVU 用于存储PV节点的初始电压
PVU=zeros(n-H-
1
,
2
);
t=
0
;
for
i=
1
:n
if
B2(i,
2
)==
2
t=t+
1
;
PVU(t,
1
)=B2(i,
5
);
PVU(t,
2
)=B2(i,
6
);
end
end
%disp(
’PV节点初始值:电压、相位ei、fi：’
) ;
PVU;
%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量 书p58
11
-
46
和
11
-
47
h=
0
;
for
i=
1
:n
%对PQ节点的处理
if
i~=isb&&B2(i,
2
)==
1
h=h+
1
;
DetaS(
2
*h-
1
,
1
)=B2(i,
3
)-OrgS(
2
*h-
1
,
1
);
DetaS(
2
*h,
1
)=B2(i,
4
)-OrgS(
2
*h,
1
);
end
end
t=
0
;
for
i=
1
:n
%对PV节点的处理
if
i~=isb&&B2(i,
2
)==
2
h=h+
1
;
t=t+
1
;
DetaS(
2
*h-
1
,
1
)=B2(i,
3
)-OrgS(
2
*h-
1
,
1
);
DetaS(
2
*h,
1
)=PVU(t,
1
)^
2
+PVU(t,
2
)^
2
-B2(i,
5
)^
2
-B2(i,
6
)^
2
;
end
end
%disp(
’P、Q、V不平衡量：’
) ;
DetaS;
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%
%
创建I，用于存储节点电流参数（计算雅克比矩阵用 具体见纸板公式推导）
I=zeros(n-
1
,
1
);
h=
0
;
for
i=
1
:n
if
i~=isb
h=h+
1
;
I(h,
1
)=(OrgS(
2
*h-
1
,
1
)-OrgS(
2
*h,
1
)*sqrt(-
1
))/conj(B2(i,
5
)+B2(i,
6
)*sqrt(-
1
));
end
end
I;
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%
%创建Jacbi(雅可比矩阵) 各元素是书中各元素求法乘个-
1
所以△W=J△V
Jacbi=zeros(
2
*n-
2
);
h=
0
;
k=
0
;
for
i=
1
:n
%对PQ节点的处理
if
B2(i,
2
)==
1
h=h+
1
;
for
j=
1
:n
if
j~=isb
k=k+
1
;
if
i==j %对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)=-B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
)+imag(I(h,
1
));
Jacbi(
2
*h-
1
,
2
*k)=G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
)+real(I(h,
1
));
Jacbi(
2
*h,
2
*k-
1
)=-Jacbi(
2
*h-
1
,
2
*k)+
2
*real(I(h,
1
));
Jacbi(
2
*h,
2
*k)=Jacbi(
2
*h-
1
,
2
*k-
1
)-
2
*imag(I(h,
1
));
else
%非对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)=-B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
);
Jacbi(
2
*h-
1
,
2
*k)=G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
);
Jacbi(
2
*h,
2
*k-
1
)=-Jacbi(
2
*h-
1
,
2
*k);
Jacbi(
2
*h,
2
*k)=Jacbi(
2
*h-
1
,
2
*k-
1
);
end
if
k==(n-
1
) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=
0
;
end
end
end
end
end
k=
0
;
for
i=
1
:n
%对PV节点的处理
if
B2(i,
2
)==
2
h=h+
1
;
for
j=
1
:n
if
j~=isb
k=k+
1
;
if
i==j %对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)= -B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
)+imag(I(h,
1
));
Jacbi(
2
*h-
1
,
2
*k)= G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
)+real(I(h,
1
));
Jacbi(
2
*h,
2
*k-
1
)=
2
*B2(i,
6
);
Jacbi(
2
*h,
2
*k)=
2
*B2(i,
5
);
else
%非对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)= -B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
);
Jacbi(
2
*h-
1
,
2
*k)= G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
);
Jacbi(
2
*h,
2
*k-
1
)=
0
;
Jacbi(
2
*h,
2
*k)=
0
;
end
if
k==(n-
1
) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=
0
;
end
end
end
end
end
%disp(
’雅克比矩阵为：’
) ;
Jacbi; %列参数跟书上的相反 这里面第一列是对电压相位的偏导 第二列是对幅值的偏导
%求解修正方程，获取节点电压的不平衡量
DetaU=zeros(
2
*n-
2
,
1
);
DetaU=JacbiDetaS;
%disp(
’求解修正方程得：’
) ;
DetaU;
%修正节点电压
j=
0
;
for
i=
1
:n
%对PQ节点处理
if
B2(i,
2
)==
1
j=j+
1
;
B2(i,
5
)=B2(i,
5
)+DetaU(
2
*j,
1
);
B2(i,
6
)=B2(i,
6
)+DetaU(
2
*j-
1
,
1
);
end
end
for
i=
1
:n
%对PV节点的处理
if
B2(i,
2
)==
2
j=j+
1
;
B2(i,
5
)=B2(i,
5
)+DetaU(
2
*j,
1
);
B2(i,
6
)=B2(i,
6
)+DetaU(
2
*j-
1
,
1
);
end
end
for
i=
1
:n
%对PQ(V)节点的处理
if
i==PQV
B2(i,
2
)=
3
;
end
end
for
i=
1
:n
%对PI节点的处理
if
i==PI
B2(i,
2
)=
4
;
end
end
%disp(
’修正后的B2阵：’
) ;
B2;
%开始循环**********************************************************************
while
max(abs(DetaU))>pr%先取绝对值在找出最大值来与误差精度比较
for
i=
1
:n
%对PQ(V)节点的处理
if
i~=isb&&B2(i,
2
)==
3
Q(i)=-(B2(i,
5
))^
2
/xp+(-(B2(i,
5
))^
2
+sqrt((B2(i,
5
))^
4
-
4
*(B2(i,
3
))^
2
*x^
2
))/
2
*x;
B2(i,
4
)=Q(i);
end
if
B2(i,
2
)==
3
B2(i,
2
)=
1
;
end
end
for
i=
1
:n
%对PI节点的处理
if
i~=isb&&B2(i,
2
)==
4
Q(i)=sqrt(Ig^
2
*((B2(i,
5
))^
2
)-B(i,
3
)^
2
);
B2(i,
4
)=Q(i);
end
if
B2(i,
2
)==
4
B2(i,
2
)=
1
;
end
end
OrgS=zeros(
2
*n-
2
,
1
);
h=
0
;
j=
0
;
for
i=
1
:n
%对PQ节点的处理
if
i~=isb&&B2(i,
2
)==
1
h=h+
1
;
for
j=
1
:n
OrgS(
2
*h-
1
,
1
)=OrgS(
2
*h-
1
,
1
)+B2(i,
5
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))+B2(i,
6
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %Pi
OrgS(
2
*h,
1
)=OrgS(
2
*h,
1
)+B2(i,
6
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))-B2(i,
5
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
)); %Qi
end
end
end
for
i=
1
:n
%对PV节点的处理
if
i~=isb&&B2(i,
2
)==
2
h=h+
1
;
for
j=
1
:n
OrgS(
2
*h-
1
,
1
)=OrgS(
2
*h-
1
,
1
)+ B2(i,
5
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))+B2(i,
6
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
));
OrgS(
2
*h,
1
)=OrgS(
2
*h,
1
)+ B2(i,
6
)*(G(i,j)*B2(j,
5
)-B(i,j)*B2(j,
6
))-B2(i,
5
)*(G(i,j)*B2(j,
6
)+B(i,j)*B2(j,
5
));
end
end
end
%disp(
’修正后的迭代计算PQ、PV节点参数：’
) ;
OrgS;
%创建DetaS
h=
0
;
for
i=
1
:n
%对PQ节点的处理
if
i~=isb&&B2(i,
2
)==
1
h=h+
1
;
DetaS(
2
*h-
1
,
1
)=B2(i,
3
)-OrgS(
2
*h-
1
,
1
);
DetaS(
2
*h,
1
)=B2(i,
4
)-OrgS(
2
*h,
1
);
end
end
t=
0
;
for
i=
1
:n
%对PV节点的处理
if
i~=isb&&B2(i,
2
)==
2
h=h+
1
;
t=t+
1
;
DetaS(
2
*h-
1
,
1
)=B2(i,
3
)-OrgS(
2
*h-
1
,
1
);
DetaS(
2
*h,
1
)=PVU(t,
1
)^
2
+PVU(t,
2
)^
2
-B2(i,
5
)^
2
-B2(i,
6
)^
2
;
end
end
%disp(
’修正后的迭代计算PQ、PV节点不平衡量：’
) ;
DetaS;
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%创建I
I=zeros(n-
1
,
1
);
h=
0
;
for
i=
1
:n
if
i~=isb
h=h+
1
;
I(h,
1
)=(OrgS(
2
*h-
1
,
1
)-OrgS(
2
*h,
1
)*sqrt(-
1
))/conj(B2(i,
5
)+B2(i,
6
)*sqrt(-
1
));
end
end
I;
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%%%
%创建Jacbi
Jacbi=zeros(
2
*n-
2
);
h=
0
;
k=
0
;
for
i=
1
:n
%对PQ节点的处理
if
B2(i,
2
)==
1
h=h+
1
;
for
j=
1
:n
if
j~=isb
k=k+
1
;
if
i==j %对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)=-B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
)+imag(I(h,
1
));
Jacbi(
2
*h-
1
,
2
*k)=G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
)+real(I(h,
1
));
Jacbi(
2
*h,
2
*k-
1
)=-Jacbi(
2
*h-
1
,
2
*k)+
2
*real(I(h,
1
));
Jacbi(
2
*h,
2
*k)=Jacbi(
2
*h-
1
,
2
*k-
1
)-
2
*imag(I(h,
1
));
else
%非对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)=-B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
);
Jacbi(
2
*h-
1
,
2
*k)=G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
);
Jacbi(
2
*h,
2
*k-
1
)=-Jacbi(
2
*h-
1
,
2
*k);
Jacbi(
2
*h,
2
*k)=Jacbi(
2
*h-
1
,
2
*k-
1
);
end
if
k==(n-
1
) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=
0
;
end
end
end
end
end
k=
0
;
for
i=
1
:n
%对PV节点的处理
if
B2(i,
2
)==
2
h=h+
1
;
for
j=
1
:n
if
j~=isb
k=k+
1
;
if
i==j %对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)= -B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
)+imag(I(h,
1
));
Jacbi(
2
*h-
1
,
2
*k)= G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
)+real(I(h,
1
));
Jacbi(
2
*h,
2
*k-
1
)=
2
*B2(i,
6
);
Jacbi(
2
*h,
2
*k)=
2
*B2(i,
5
);
else
%非对角元素的处理
Jacbi(
2
*h-
1
,
2
*k-
1
)= -B(i,j)*B2(i,
5
)+G(i,j)*B2(i,
6
);
Jacbi(
2
*h-
1
,
2
*k)= G(i,j)*B2(i,
5
)+B(i,j)*B2(i,
6
);
Jacbi(
2
*h,
2
*k-
1
)=
0
;
Jacbi(
2
*h,
2
*k)=
0
;
end
if
k==(n-
1
) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=
0
;
end
end
end
end
end
%disp(
’修正后的雅克比矩阵：’
) ;
Jacbi;
DetaU=zeros(
2
*n-
2
,
1
);
DetaU=JacbiDetaS;
DetaU;
%修正节点电压
j=
0
;
for
i=
1
:n
%对PQ节点处理
if
B2(i,
2
)==
1
j=j+
1
;
B2(i,
5
)=B2(i,
5
)+DetaU(
2
*j,
1
);
B2(i,
6
)=B2(i,
6
)+DetaU(
2
*j-
1
,
1
);
end
end
for
i=
1
:n
%对PV节点的处理
if
B2(i,
2
)==
2
j=j+
1
;
B2(i,
5
)=B2(i,
5
)+DetaU(
2
*j,
1
);
B2(i,
6
)=B2(i,
6
)+DetaU(
2
*j-
1
,
1
);
end
end
for
i=
1
:n
%对PQ(V)节点的处理
if
i==PQV
B2(i,
2
)=
3
;
end
end
for
i=
1
:n
%对PI节点的处理
if
i==PI
B2(i,
2
)=
4
;
end
end
Times=Times+
1
; %迭代次数加
1
end
disp(
’初始潮流计算结果如下：’
);
disp(
’迭代次数：’
);
Times
disp(
’节点电压幅值如下（按节点由小到大的顺序）：’
);
B2(
:
,
5
)
%创建Sb，用于存储平衡节点功率
Sb=zeros(
1
);
for
i=
1
:n
if
i==isb
for
j=
1
:n
Sb=Sb+(B2(i,
5
)+sqrt(-
1
)*B2(i,
6
))*conj(Y(i,j))*conj(B2(j,
5
)+sqrt(-
1
)*B2(j,
6
));
end
end
end
disp(
’初始平衡节点功率：’
)
Sb
%求解各条支路的功率
Sij=zeros(n);
for
i=
1
:n
for
j=
1
:n
Sij(i,j)=Sij(i,j)+(B2(i,
5
)+sqrt(-
1
)*B2(i,
6
))^
2
*conj(B1(
1
,
4
))+(B2(i,
5
)+sqrt(-
1
)*B2(i,
6
))*conj(Y(i,j))*(conj((B2(i,
5
)+sqrt(-
1
)*B2(i,
6
)))-conj((B2(j,
5
)+sqrt(-
1
)*B2(j,
6
))));
end
end
%disp(
’线路功率分布：’
)
Sij;
Sij_P=real(Sij);
Sij_Q=imag(Sij);
for
i=
1
:
32
P_flow(i)= abs( Sij_P(B1(i,
1
),B1(i,
2
)));
end
for
i=
1
:
33
for
j=
1
:
33
if
Sij_P(i,j)<
0
.
0001
Sij_P(i,j)=
0
;
end
if
Sij_Q(i,j)<
0
.
0001
Sij_Q(i,j)=
0
;
end
end
end
%%%
求解系统网损
Ploss=
0
;
for
i=
1
:n
for
j=
1
:n
Ploss=Ploss+B2(i,
5
)*B2(j,
5
)*(G(i,j)*cos(B2(i,
6
)-B2(j,
6
))+B(i,j)*sin(B2(i,
6
)-B2(j,
6
)));
end
end
Ploss=Ploss*
10000
;
V_amp=B2(
:
,
5
);
for
i=
1
:n
%电压罚函数
V_33(i)=abs(
1
-V_amp(i)); %求与额定电压的差值
end
V_flu=sum(V_33);
end
⛳️ 运行结果
参考文献
[1] 龚莉莉.基于利益链的配电网电动汽车充电优化策略研究[J].[2023-11-14].
[2] 李敏,苏小林,阎晓霞,等.多目标分层分区的电动汽车有序充放电优化控制[J].电网技术, 2015, 39(12):7.DOI:10.13335/j.1000-3673.pst.2015.12.033.