MATLAB简单潮流计算程序如下：
function lianxuchaoliu
clear
clc
n=9%节点数；
nl=9%支路数；
isb=1%平衡节点号；
pr=0.00001%误猛桐局差精度；
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首轮歼屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的枝让程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
function lianxuchaoliuclear
clc
n=9%节点数；
nl=9%支弊铅游路数；
isb=1%平衡节点号；
pr=0.00001%误差精度；
b1=[1 4 0.0576i 0 1.05 14 5 0.017+0.092i 0.158i 1 05 6 0.039+0.17i 0.358i 1 03 6 0.0586i 0 1.05 16 7 0.0119+0.1008i 0.209i 1 07 8 0.0085+0.072i 0.149i 1 08 2 0.0625i 0 1.05 18 9 0.032+0.161i 0.306i 1 09 4 0.01+0.085i 0.176i 1 0]
%依次是支路首端；末端，支路阻抗；对地电纳；支路变比；折算到哪一侧标志（高压侧为1；低压侧为0）；
b2=[0 0 1.05 1.05 0 11.63 0 1.05 1.05 0 30.85 0 1.05 1.05 0 30 0 1 0 0 20 0.9+0.3i 1 0 0 20 0 1 0 0 20 1+0.35i 1 0 0 20 0 1 0 0 20 1.25+0.5i 1 0 0 2]
%依次是节点的发电机功率；负荷功率；节租销点电压初激改值；PV节点电压V给定值；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；
Y=zeros(n)%求导纳阵；
for i=1:nl
if b1(i,6)==0
p=b1(i,1)q=b1(i,2)
else
p=b1(i,2)q=b1(i,1)
end
Y(p,q)=Y(p,q)-1./(b1(i,3)*b1(i,5))
Y(q,p)=Y(p,q)
Y(q,q)=Y(q,q)+1./(b1(i,3)*b1(i,5)^2)+b1(i,4)./2
Y(p,p)=Y(p,p)+1./b1(i,3)+b1(i,4)./2
end
%disp('系统的导纳阵为：')
%disp(Y)
G=real(Y)B=imag(Y)
for i=1:n
e(i)=real(b2(i,3))
f(i)=imag(b2(i,3))
v(i)=b2(i,4)
end
for i=1:n
S(i)=b2(i,1)-b2(i,2)
B(i,i)=B(i,i)+b2(i,5)
end
P=real(S)Q=imag(S)
w=zeros(2*n,1)Jac=zeros(2*n)
t=0
while t==0
for i=1:n
if b2(i,6)~=isb
C=0D=0
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j)
D=D+G(i,j)*f(j)+B(i,j)*e(j)
end
if b2(i,6)==2%P，Q节点；
w(2*i)=P(i)-e(i)*C-f(i)*D
w(2*i-1)=Q(i)-f(i)*C+e(i)*D
else if b2(i,6)==3%P，V节点；
w(2*i)=P(i)-e(i)*C-f(i)*D
w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2)
end
end
else
w(2*i-1)=0
w(2*i)=0
end
end
%disp(w)
w1=w(3:2*n)
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P，Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P，V节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=0
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-2*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
Jac(2*i-1,2*j-1)=0
Jac(2*i,2*j)=0
Jac(2*i-1,2*j)=0
Jac(2*i,2*j-1)=0
end
end
end
%disp(Jac)
Jac1=Jac(3:2*n,3:2*n)
for k=1:2*n-2
m=0
for i=k+1:2*n-2
m=Jac1(i,k)./Jac1(k,k)
w1(i)=w1(i)-m*w1(k)
for j=k+1:2*n-2
Jac1(i,j)=Jac1(i,j)-m*Jac1(k,j)
end
end
end
E(2*n-2)=w1(2*n-2)./Jac1(2*n-2,2*n-2)
for i=2*n-3:-1:1
c=0
for k=i+1:2*n-2
c=c+Jac1(i,k)*E(k)
E(i)=(w1(i)-c)./Jac1(i,i)
end
end
%disp(E)
for i=1:n-1
e(i+1)=e(i+1)-E(2*i-1)
f(i+1)=f(i+1)-E(2*i)
end
%disp(e)
%disp(f)
for i=1:n-1
b(2*i-1)=abs(E(2*i-1))
b(2*i)=abs(E(2*i))
end
KB=max(b)
%disp(KB)
if KB<pr
t=1
else
t=0
end
end
%disp(e)
%disp(f)
for i=1:n
fz(i)=sqrt(e(i)^2+f(i)^2)
end
disp(fz)
JJ1=zeros(2*n-1,2*n-1)Jac1=zeros(2*n-2,2*n-2)ttt=0ccc=1%迭代次数；
L=0.05T=0%L是步长，T是参数；
while ttt==0
%确定扩展矩阵
for i=1:n-1
if b2(i+1,6)==3
JJ1(2*i,2*n-1)=1.5*1.05
JJ1(2*i-1,2*n-1)=0
else if b2(i+1,6)==2
if b2(i+1,2)==0
JJ1(2*i-1,2*n-1)=0
JJ1(2*i,2*n-1)=0
else
JJ1(2*i,2*n-1)=-real(b2(i+1,2))
JJ1(2*i-1,2*n-1)=-imag(b2(i+1,2))
end
end
end
end
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2% P,Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P,V节点，
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=0
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-2*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
Jac(2*i-1,2*j-1)=0
Jac(2*i,2*j)=0
Jac(2*i-1,2*j)=0
Jac(2*i,2*j-1)=0
end
end
end
Jac1=Jac(3:2*n,3:2*n)
for i=1:2*n-2
for j=1:2*n-2
JJ1(i,j)=Jac1(i,j)
end
end
if ccc==1%以负荷为连续参数；
for j=1:2*n-2
JJ1(2*n-1,j)=0
end
JJ1(2*n-1,2*n-1)=1
w2=zeros(2*n-1,1)
for i=1:2*n-2
w2(i,1)=0
end
w2(2*n-1,1)=1
end
%disp(JJ1)
if ccc>1%以切向量中分量最大（绝对值最大）的状态变量选定为连续参数；
for i=1:2*n-2
Jd(i)=abs(d(i))
end
for i=1:2*n-3
if Jd(i)>=Jd(i+1)
zd=Jd(i)ek=i
else if Jd(i)<=Jd(i+1)
zd=Jd(i+1)ek=i+1
end
end
end
for j=1:2*n-1
JJ1(2*n-1,j)=0
end
JJ1(2*n-1,ek)=1
w2=zeros(2*n-1,1)
for i=1:2*n-1
w2(i,1)=0
end
if d(ek)>0
w2(ek,1)=1
else if d(ek)<0
w2(ek,1)=-1
end
end
end
for k=1:2*n-1
m=0
for i=k+1:2*n-1
m=JJ1(i,k)./JJ1(k,k)
w2(i)=w2(i)-m*w2(k)
for j=k+1:2*n-1
JJ1(i,j)=JJ1(i,j)-m*JJ1(k,j)
end
end
end
d(2*n-1)=w2(2*n-1)./JJ1(2*n-1,2*n-1)
for i=2*n-2:-1:1
c=0
for k=i+1:2*n-1
c=c+JJ1(i,k)*d(k)
d(i)=(w2(i)-c)./JJ1(i,i)
end
end
%disp(d)
for i=1:n-1
e(i+1)=e(i+1)+L*d(2*i-1)
f(i+1)=f(i+1)+L*d(2*i)
end
T=T+L*d(2*n-1)
disp(d(2*n-1))
%disp(e)
%disp(f)
%disp(T)
%对预估的近似解进行校正；
tt=0
while tt==0
for i=1:n
if b2(i,6)~=isb
C=0D=0
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j)
D=D+G(i,j)*f(j)+B(i,j)*e(j)
end
if b2(i,6)==2% P,Q节点；
if b2(i,2)~=0
wj(2*i)=P(i)-T*real(b2(i,2))-e(i)*C-f(i)*D
wj(2*i-1)=Q(i)-T*imag(b2(i,2))-f(i)*C+e(i)*D
else
wj(2*i)=P(i)-e(i)*C-f(i)*D
wj(2*i-1)=Q(i)-f(i)*C+e(i)*D
end
else if b2(i,6)==3%P,V节点；
wj(2*i)=P(i)+T*1.5*1.05-e(i)*C-f(i)*D
wj(2*i-1)=0
end
end
else
wj(2*i)=0
wj(2*i-1)=0
end
end
wj1=wj(3:2*n)
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P,Q节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
JJac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
JJac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
JJac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P,V节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i-1,2*j)=0
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
JJac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
JJac(2*i-1,2*j)=-2*f(i)
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
JJac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
JJac(2*i-1,2*j-1)=0
JJac(2*i,2*j)=0
JJac(2*i-1,2*j)=0
JJac(2*i,2*j-1)=0
end
end
end
JJac1=JJac(3:2*n,3:2*n)
for k=1:2*n-2
m=0
for i=k+1:2*n-2
m=JJac1(i,k)./JJac1(k,k)
wj1(i)=wj1(i)-m*wj1(k)
for j=k+1:2*n-2
JJac1(i,j)=JJac1(i,j)-m*JJac1(k,j)
end
end
end
E1(2*n-2)=wj1(2*n-2)./JJac1(2*n-2,2*n-2)
for i=2*n-3:-1:1
c=0
for k=i+1:2*n-2
c=c+JJac1(i,k)*E1(k)
E1(i)=(wj1(i)-c)./JJac1(i,i)
end
end
%disp(E1)
for i=1:n-1
e(i+1)=e(i+1)-E1(2*i-1)
f(i+1)=f(i+1)-E1(2*i)
end
for i=1:n-1
bx(2*i-1)=abs(E1(2*i-1))
bx(2*i)=abs(E1(2*i))
end
KB1=max(bx)
if KB1<pr
tt=1
else
tt=0
end
end
%disp(e)
%disp(f)
if d(2*n-1)>0
ttt=0ccc=ccc+1
else
ttt=1
end
end
disp(T)
disp(ccc)
%disp(e)
%disp(f)
for i=1:n
fz1(i)=sqrt(e(i)^2+f(i)^2)
end
disp(fz1)
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）