%% 程序初始化
clc,clear;
%% 读取电力系统节点和线路参数
eps=power13;
baseMVA=eps.baseMVA;
basekV=eps.basekV;
bus=eps.bus;
generator=eps.generator;
branch=eps.branch;
%% 形成节点导纳矩阵
%统计节点、发电机、线路数量
[Number_bus,~]=size(bus); %记为n
[Number_generator,~]=size(generator);
[Number_branch,~]=size(branch);
%统计PQ节点、PV节点、平衡节点数量
Number_PQ=0; %记为m
Number_PV=0; %记为n-1-m
Number_reference=0;
for i=1:Number_bus
if bus(i,2)==1
Number_PQ=Number_PQ+1;
elseif bus(i,2)==2
Number_PV=Number_PV+1;
else
Number_reference=Number_reference+1;
end
end
%形成支路导纳矩阵Yb
Yb=zeros(Number_branch,1);
for i=1:Number_branch
Yb(i)=1/(branch(i,4)+1j*branch(i,5));
end
%形成节点导纳矩阵Y
Y=zeros(Number_bus,Number_bus);
for i=1:Number_branch
p=branch(i,2);
q=branch(i,3);
Y(p,p)=Y(p,p)+Yb(i);
Y(q,q)=Y(q,q)+Yb(i);
Y(p,q)=-Yb(i);
Y(q,p)=-Yb(i);
end
Y=Y*baseMVA/(basekV^2); %换算成有名值
G=real(Y);
B=imag(Y);
%% 设定状态量初值
%设定电压相角初值
angle=bus(1:Number_bus,6)*pi/180; %换算成弧度
%设定电压幅值初值
V=bus(1:Number_bus,5);
V=V*basekV; %换算成有名值
%给定节点注入有功列向量
Psp=-bus(1:Number_PQ+Number_PV,3);
Psp(Number_PQ+1:Number_PQ+Number_PV)=Psp(Number_PQ+1:Number_PQ+Number_PV)+generator(Number_reference:Number_reference+1,3);
%给定节点注入无功列向量
Qsp=-bus(1:Number_PQ,4);
%% 设置迭代次数
for i=1:50
%% 计算功率偏差
%计算有功功率偏差
dP=Psp;
for p=1:Number_PQ+Number_PV %1到PQ+PV节点
for q=1:Number_bus %1到所有节点
dP(p)=dP(p)-V(p)*V(q)*(G(p,q)*cos(angle(p)-angle(q))+B(p,q)*sin(angle(p)-angle(q)));
end
end
%计算无功功率偏差
dQ=Qsp;
for p=1:Number_PQ %1到PQ节点
for q=1:Number_bus %1到所有节点
dQ(p)=dQ(p)-V(p)*V(q)*(G(p,q)*sin(angle(p)-angle(q))-B(p,q)*cos(angle(p)-angle(q)));
end
end
%形成功率偏差
dF=[dP;dQ];
%% 计算雅可比矩阵
H=zeros(Number_bus-1,Number_bus-1); %(n-1)维
N=zeros(Number_bus-1,Number_PQ); %(n-1)行m列
M=zeros(Number_PQ,Number_bus-1); %m行(n-1)列
L=zeros(Number_PQ,Number_PQ); %m维
%计算H
for p=1:Number_bus-1
for q=1:Number_bus-1
if p==q
for s=1:Number_bus
if s~=p
H(p,p)=H(p,p)+V(p)*V(s)*(G(p,s)*sin(angle(p)-angle(s))-B(p,s)*cos(angle(p)-angle(s)));
end
end
else
H(p,q)=-V(p)*V(q)*(G(p,q)*sin(angle(p)-angle(q))-B(p,q)*cos(angle(p)-angle(q)));
end
end
end
%计算N
for p=1:Number_bus-1
for q=1:Number_PQ
if p==q
for s=1:Number_bus
if s~=p
N(p,p)=N(p,p)-V(p)*V(s)*(G(p,s)*cos(angle(p)-angle(s))+B(p,s)*sin(angle(p)-angle(s)));
end
end
N(p,p)=N(p,p)-2*V(p).^2*G(p,p);
else
N(p,q)=-V(p)*V(q)*(G(p,q)*cos(angle(p)-angle(q))+B(p,q)*sin(angle(p)-angle(q)));
end
end
end
%计算M
for p=1:Number_PQ
for q=1:Number_bus-1
if p==q
for s=1:Number_bus
if s~=p
M(p,p)=M(p,p)-V(p)*V(s)*(G(p,s)*cos(angle(p)-angle(s))+B(p,s)*sin(angle(p)-angle(s)));
end
end
else
M(p,q)=V(p)*V(q)*(G(p,q)*cos(angle(p)-angle(q))+B(p,q)*sin(angle(p)-angle(q)));
end
end
end
%计算L
for p=1:Number_PQ
for q=1:Number_PQ
if p==q
for s=1:Number_bus
if s~=p
L(p,p)=L(p,p)-V(p)*V(s)*(G(p,s)*sin(angle(p)-angle(s))-B(p,s)*cos(angle(p)-angle(s)));
end
end
L(p,p)=L(p,p)+2*V(p).^2*B(p,p);
else
L(p,q)=-V(p)*V(q)*(G(p,q)*sin(angle(p)-angle(q))-B(p,q)*cos(angle(p)-angle(q)));
end
end
end
%形成雅可比矩阵
J=[H,N;M,L]; %(n-1+m)维
%% 计算修正量
dX=-(J\dF); %(n-1+m)行1列
dangle=dX(1:Number_bus-1);
dV=dX(Number_bus:Number_bus-1+Number_PQ).*V(1:Number_PQ); %由于dX=dV/V，所以dV=dX*V
%% 更新状态量
angle=angle+[dangle;0];
V=V+[dV;zeros(Number_bus-Number_PQ,1)];
%% 收敛判据判断
if max(abs(dX))<=1e-4
break;
end
end
%% 输出潮流计算结果
P=zeros(Number_bus,1);
Q=zeros(Number_bus,1);
for p=1:Number_bus
for q=1:Number_bus
%节点有功功率
P(p)=P(p)+real(V(p)*exp(angle(p)*1j)*conj(Y(p,q)*V(q)*exp(angle(q)*1j)));
%节点无功功率
Q(p)=Q(p)+imag(V(p)*exp(angle(p)*1j)*conj(Y(p,q)*V(q)*exp(angle(q)*1j)));
end
end
%支路功率
SS=zeros(Number_bus,Number_bus);
for p=1:Number_bus
for q=1:Number_bus
SS(p,q)=V(p)*exp(angle(p)*1j)*conj((V(p)*exp(angle(p)*1j)-V(q)*exp(angle(q)*1j))*Y(p,q));
end
end
S=[SS(13,1);SS(1,2);SS(1,5);SS(1,3);SS(5,6);SS(5,11);SS(5,9);SS(2,12);SS(3,4);SS(6,7);SS(6,8);SS(9,10)];
%支路功率损耗
ddSS=zeros(Number_bus,Number_bus);
for p=1:Number_bus
for q=1:Number_bus
ddSS(p,q)=SS(p,q)+SS(q,p);
end
end
dS=[ddSS(13,1);ddSS(1,2);ddSS(1,5);ddSS(1,3);ddSS(5,6);ddSS(5,11);ddSS(5,9);ddSS(2,12);ddSS(3,4);ddSS(6,7);ddSS(6,8);ddSS(9,10)];
%节点电压幅值
V=V/basekV; %换算成标幺值
%节点电压相角
angle=angle*180/pi; %换算成角度
评论0