潮流计算的计算机算法实验报告
姓名：学号：班级：
一、实验目的
掌握潮流计算的计算机算法。
熟悉MATLAB，并掌握MATLAB程序的基本调试方法。
二、实验准备
根据课程内容，熟悉MATLAB软件的使用方法，自行学习MATLAB程序的基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法(或PQ分解法) 的计算程序，用相应的算例在MATLAB上进行计算、调试和验证。
三、实验要求
每人一组，在实验课时内，用MATLAB调试和修改运行程序，用算例计算输出潮流结果。
四、实验程序
clear; %清空内存
n=input('请输入节点数:n=');
n1=input('请输入支路数:n1=');
isb=input('请输入平衡节点号:isb=');
pr=input('请输入误差精度:pr=');
B1=input('请输入支路参数:B1=');
B2=input('请输入节点参数:B2=');
X=input('节点号和对地参数:X=');
Y=zeros(n);
Times=1;
%一：创建节点导纳矩阵
for i=1:n1
if B1(i,6)==0 %不含变压器的支路
p=B1(i,1);
q=B1(i,2);
Y(p,q)=Y(p,q)-1/B1(i,3);
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);
Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);
else %含有变压器的支路
p=B1(i,1);
q=B1(i,2);
Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1/B1(i,3);
Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));
end
end
Y;
%将OrgS、DetaS初始化
OrgS=zeros(2*n-2,1);
DetaS=zeros(2*n-2,1);
%二：创建OrgS，用于存储初始功率参数
h=0;
j=0;
for i=1:n %对PQ节点的处理
if i~=isb&B2(i,6)==2
h=h+1;
for j=1:n
OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j ,3))-imag(Y(i,j))*I
mag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i, j))*real(B2(j,3)));
OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3)) -imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2( j,3))+imag(Y(i,j))*real(B2(j,3)));
end
end
end
%三：对PV节点的处理，注意这时不可再将h初始化为0
for i=1:n
if i~=isb&B2(i,6)==3
h=h+1;
for j=1:n
OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j ,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag (B2(j,3))+imag(Y(i,j))*real(B2(j,3)));
OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3)) -imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2( j,3))+imag(Y(i,j))*real(B2(j,3)));
end
end
end
OrgS;
%四：创建PVU 用于存储PV节点的初始电压
PVU=zeros(n-h-1,1);
t=0;
for i=1:n
if B2(i,6)==3
t=t+1;
PVU(t,1)=B2(i,3);
end
end
PVU;
%五：创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;
for i=1:n %对PQ节点的处理
if i~=isb&B2(i,6)==2
h=h+1;
DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);
end
end
t=0;
for i=1:n
%六：对PV节点的处理，注意这时不可再将h初始化为0
if i~=isb&B2(i,6)==3
h=h+1;
t=t+1;
DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);
DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,
3))^2-imag( B2(i,3))^2;
end
end
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,3));
end
end
I;
%八：创建Jacbi(雅可比矩阵)
Jacbi=zeros(2*n-2);
h=0;
k=0;
for i=1:n %对PQ节点的处理
if B2(i,6)==2
h=h+1;
for j=1:n
if j~=isb
k=k+1;
if i==j %对角元素的处理
