MATLAB电力系统潮流计算
用MATLAB利用牛顿-拉夫逊法计算电科院22节点系统（使用稀疏技术），计算出各节点的电压相角、各支路的功率等。使用了稀疏计算技术，能计算1000节点以上的系统，且计算时间小于1秒。
计算思路
- 形成节点导纳矩阵，分析各节点的类型，找出待求量；
- 设定各节点电压的初值，并给定迭代误差判据；；
- 计算不平衡量；
- 求修正方程的系数矩阵，即雅可比矩阵的各元素；
- 形成雅可比矩阵，求解修正方程，得到各节点电压变化量；
- 计算各节点电压新值；
- 检查是否收敛，如不收敛，则运用各节点电压新值自第3步开始下一次迭代；
- 计算平衡节点功率和线路功率。
结果展示
22节点网络潮流计算
1047节点潮流计算
部分代码
读取文件
%% 取数据文件和参数
[open_filename] = uigetfile('.TXT');
data=dlmread(open_filename);
node_num=data(1,1); %节点数
branchnum=data(1,2); %支路数
balancenode=data(1,3); %平衡节点号数
v_balancenode=data(1,4); %平衡节点电压
base_power=data(1,6); %基准功率
convergence_accuracy=data(1,7); %运算精度
index=find(data(:,1)==0); %index存放矩阵中第一列为0的元素（区分各块参数）
Line=data(index(1)+1:index(2)-1,:); %线路为第一块
Line_ground=data(index(2)+1:index(3)-1,:); %对地支路为第二块
Transformer=data(index(3)+1:index(4)-1,1:6); %变压器为第三块
PQ=data(index(4)+1:index(5)-1,:); %PQ节点为第四块
PV=data(index(5)+1:index(6)-1,:); %PV节点为第五块
PV_num=PV(:,1); %PV节点数
R=Line(:,4);X=1j*Line(:,5); B=1j*Line(:,6);Z=R+X;Line_Y=1./Z; %Y是Z的逆矩阵,线路参数
Rt=Transformer(:,4); Xt=Transformer(:,5); k=Transformer(:,6); %变压器参数
Gt=Rt./(Rt.*Rt+Xt.*Xt);Bt=-Xt./(Rt.*Rt+Xt.*Xt);
Yt=Gt+1j.*Bt; Yt2=Yt./k; Yt3=Yt./k./k; %变压器等值电路参数
G_ground=1j*Line_ground(:,2); %对地支路参数
Nonzero_num=node_num+2*(size(Line,1)+size(Transformer,1)); %非零元素：对角线+2*（线路+变压器）
形成节点导纳矩阵（应用稀疏矩阵）
tic; %计时
Y=spalloc(node_num,node_num,Nonzero_num); %创建一个节点数*节点数的全零稀疏矩阵Y，该矩阵具有Nonzero_num个非零值的空间
node_i=Line(:,2);node_j=Line(:,3); %节点i、j
Y=Y-sparse(node_i,node_j,Line_Y,node_num,node_num)-sparse(node_j,node_i,Line_Y,node_num,node_num); %把互导纳放进Y的（node_i,node_j）和（node_j,node_i）位置
Y_new=Y;
i=1:node_num; %创建一个由1到22的行向量
Y=Y-sparse(i,i,sum(Y_new,2),node_num,node_num); %把自导纳放进Y的（i,i）位置
node_i=Transformer(:,2);node_j=Transformer(:,3); %将变压器等值参数电路放进矩阵
Y=Y-sparse(node_i,node_j,Yt2,node_num,node_num)-sparse(node_j,node_i,Yt2,node_num,node_num); %互导纳
Y=Y+sparse(node_i,node_i,Yt3,node_num,node_num)+sparse(node_j,node_j,Yt,node_num,node_num); %自导纳
i=Line_ground(:,1);
Y=Y+sparse(i,i,G_ground,node_num,node_num); %自导纳放入对地支路的G
i=Line(:,2);
Y=Y+sparse(i,i,B,node_num,node_num); %自导纳放入i支路的B
i=Line(:,3);
Y=Y+sparse(i,i,B,node_num,node_num); %自导纳放入j支路的B
Y_abs=abs(full(Y)); %形成导纳矩阵的幅值矩阵
Y_angle=angle(full(Y)); %导纳矩阵的相角矩阵
迭代设置
V=ones(node_num,1); %计算启动的电压幅值为1
deta=zeros(node_num,1); %相角为0
%deta=ones(node_num,1)*100*pi/180; %问9：平衡节点相角设为100°
V(PV(:,1),1)=PV(:,2); %PV节点电压设置
V(balancenode)=v_balancenode; %平衡节点电压设置
JDV=[deta;V]; %修正量的向量
P=(PQ(:,2)-PQ(:,4))/base_power; %取标幺值
Q=(PQ(:,3)-PQ(:,5))/base_power;
开始迭代的（部分代码）
count_max=15; %最大迭代次数为count_max
for count=0:count_max
%% 功率方程不平衡量
A=ones(node_num,1);
deta_ij=deta*A'-A*deta'; %deta(ij)=deta(i)-deta(j)
Pi_0=sum(diag(V)*(real(Y).*cos(deta_ij)+imag(Y).*sin(deta_ij))*diag(V),2);%计算Pi
dP=P-Pi_0; %计算dPi
Qi_0=sum(diag(V)*(real(Y).*sin(deta_ij)-imag(Y).*cos(deta_ij))*diag(V),2);%计算Qi
dQ=Q-Qi_0; %计算dQi
dP(balancenode)=0;dQ(balancenode)=0; %平衡节点的dP,Q设为0
dQ(PV(:,1))=0; %pv节点的不平衡量dQ设为0
JDPQ=[dP;dQ]; %修正方程的左边，dp，dq组成的列向量
%% 雅可比矩阵（应用稀疏矩阵）非对角元部分
[Jacobian_i,Jacobian_j]=find(Y);Jacobian_n=size(Jacobian_i,1); %Y中非零元素位置放到Jacobian_i,Jacobian_j中,Jacobian_n为Jacobian_i的行数
判断收敛条件
%% 判断是否收敛
%DV=V.*delta_V;
%Ddelta=[delta_deta;DV];
delta_V=delta(node_num+1:2*node_num).*V;
V=V-delta_V;
deta=deta-delta(1:node_num);
if max(abs(delta_V))<convergence_accuracy %达到收敛精度
break
end
time=toc;
输出结果
%% 运行结束输出结果
disp('迭代次数:'),disp(count);
disp('迭代用时：'),disp(time);
if count>=count_max %不收敛
disp('超过设置的最大迭代次数，结果不收敛！')
else
disp('结果收敛，结果为：');jd_angle=deta.*180/pi; %相角弧度制转换为角度制
disp('各节点电压：（幅值，相角）'),disp([V,jd_angle]);
n_Line=size(Line,1); %线路参数长度
n_trsnsformer=size(Transformer,1); %变压器参数长度
S1=zeros(node_num);
S2=zeros(node_num);
计算平衡节点功率、各支路功率、变压器支路功率等等
%% 计算平衡节点功率
Sb=0;
for t=1:node_num
Ui=V(t,1)*cos(deta(t,1))+1j*V(t,1)*sin(deta(t,1));
Us=V(balancenode,1)*cos(deta(balancenode,1))+1j*V(balancenode,1)*sin(deta(balancenode,1));
Sb=Sb+conj(Y(balancenode,t))*conj(Ui);
end
disp('平衡节点功率：'),Sb=Us*Sb
%% 画图显示结果
subplot(2,2,1),plot(1:node_num,V)
xlabel('节点号');ylabel('电压幅值');grid on;
subplot(2,2,2),plot(1:node_num,jd_angle)
xlabel('节点号');ylabel('电压角度');grid on;
subplot(2,2,3),bar(P,'r');
xlabel('节点号');ylabel('节点注入有功');grid on;
subplot(2,2,4);bar(Q,'b');
xlabel('节点号');ylabel('节点注入无功');grid on;
end
完整源代码可私聊。
可转下载区