直角坐标牛顿-拉夫逊法潮流计算matlab程序(仅供参考)

%该程序仅针对《电力系统分析下》何仰赞P61 的4节点算例。

%节点电压用直角坐标表示时的牛顿-拉夫逊法潮流计算（matlab程序，可能有小错误，仅供学习之用，如果想要通用的程序，可自己动手改或再到pudn、csdn等网站搜索更好的）%南昌大学电力061 李圣涛2009年5月编写，2012年5月上传

clc %清空command windows

clear all %清空workspace

%为了提高可移植性、可读性、通用性，设置以下变量

N=4; %独立节点数

NPQ=2; %PQ节点数

NPV=1; %PV节点数

K=0; %迭代次数

%请输入最大迭代次数Kmax。可从0开始，以观察第Kmax次迭代的结果

Kmax=input('\n\n 请输入最大迭代次数后回车(可从零开始) Kmax=\n');

small=10^(-5); %ε不能太小

%i为节点标号，其中1号……NPQ号为PQ节点，（NPQ+1）

%号……（N-1）号为PV节点，N号节点为平衡节点

%节点导纳矩阵Y的实部

G=[1.042093 -0.588235 0 -0.453858;

-0.588235 1.069005 0 -0.480769;

0 0 0 0 ;

-0.453858 -0.480769 0 0.934627 ];

%节点导纳矩阵Y的虚部

B=[ -8.242876 2.352941 3.666667 1.891074 ;

2.352941 -4.727377 0 2.403846 ;

3.666667 0 -3.3333333 0 ;

1.891074

2.403846 0 -4.261590 ];

%Y矩阵

Y=complex(G,B);

%给定PQ节点的Pnode、Qnode，PV节点的Pnode、Vnode。（Vnode为节点电压的幅值）Pnode=[ -0.3 -0.55 0.5 ];%PQ、PV节点的初值P

Qnode=[ -0.18 -0.13 0 ];%PQ节点的初值Q

Vnode=[ 0 0 1.10 ];%PV节点的初值V

%迭代初值

e=[ 1.0 1.0 1.1 1.05 ];

f=[ 0 0 0 0 ];

%利用for循环来实现多次迭代。

for K=0:Kmax,

%计算△W。△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2】(平衡节点不参与迭代) %用公式（11-46）计算NPQ个PQ节点的△P、△Q。算法：先初始化，再不断地累减

%初始化,得△P=【Pnode(1),Pnode(2),……，Pnode(NPQ)】

for i=1:NPQ,

dP(i)=Pnode(i);

dQ(i)=Qnode(i);

end

%不断地累减,得△P和△Q

for i=1:NPQ,

for j=1:N,

dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);%△P

dQ(i)=dQ(i)-f(i)*G(i,j)*e(j)+f(i)*B(i,j)*f(j)+e(i)*G(i,j)*f(i)+e(i)*B(i,j)*e(j);%△Q

end

end

%用公式（11-47）计算NPV个PV节点的△P、（△V）^2 (其中△P的计算同上)算法：先初始化，再不断地累减

%初始化

for i=(NPQ+1):(N-1),

dP(i)=Pnode(i);

end

%不断地累减,得△P和（△V）^2

for i=(NPQ+1):(N-1),

for j=1:N,

dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);

dV2(i)=Vnode(i)^2-e(i)^2-f(i)^2;

end

end

%组合成△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2】(平衡节点不参与迭代)

%将△P间隔地赋入△W

a=1;

for i=1:(N-1),

dW(a)=dP(i);

a=a+2;

end

%将△Q间隔地赋入△W

a=2;

for i=1:NPQ,

dW(a)=dQ(i);

a=a+2;

end

%将△V^2间隔地赋入△W

a=NPQ*2+2;

for i=(NPQ+1):(NPQ+NPV),

dW(a)=dV2(i);

a=a+2;

end

%判断是否小于ε

if max(dW)

fprintf('\n 迭代是收敛的。第%d次迭代后终止迭代。\n',K-1);

break;

end

%计算雅克比矩阵各元素。算法：先初始化，再不断地累减

J=zeros(N-1);

%初始化

for i=1:N-1,

for j=1:N-1,

if i<=NPQ %求雅克比矩阵中PQ节点的元素

if i==j

J(2*i-1,2*j-1) = -G(i,i)*e(i)-B(i,i)*f(i);

J(2*i-1,2*j ) = B(i,i)*e(i)-G(i,i)*f(i);

J(2*i ,2*j-1) = B(i,i)*e(i)-G(i,i)*f(i);

J(2*i ,2*j ) = G(i,i)*e(i)+B(i,i)*f(i);

elseif i~=j

J(2*i-1,2*j-1) = -G(i,j)*e(i)-B(i,j)*f(i);

J(2*i-1,2*j ) = B(i,j)*e(i)-G(i,j)*f(i);

J(2*i ,2*j-1) = B(i,j)*e(i)-G(i,j)*f(i);

J(2*i ,2*j ) = G(i,j)*e(i)+B(i,j)*f(i);

end

end

if i>NPQ %求雅克比矩阵中PV节点的元素