电力系统分析潮流计算的计算机算法.docx
《电力系统分析潮流计算的计算机算法.docx》由会员分享,可在线阅读,更多相关《电力系统分析潮流计算的计算机算法.docx(14页珍藏版)》请在冰豆网上搜索。
电力系统分析潮流计算的计算机算法
潮流计算的计算机算法实验报告
姓名:
学号:
班级:
一、实验目的
掌握潮流计算的计算机算法。
熟悉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;
%一:
创建节点导纳矩阵
fori=1:
n1
ifB1(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;
fori=1:
n%对PQ节点的处理
ifi~=isb&B2(i,6)==2
h=h+1;
forj=1:
nOrgS(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
%三:
对PV节点的处理,注意这时不可再将h初始化为0
fori=1:
n
ifi~=isb&B2(i,6)==3
h=h+1;
forj=1:
nOrgS(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;
fori=1:
n
ifB2(i,6)==3
t=t+1;
PVU(t,1)=B2(i,3);
end
end
PVU;
%五:
创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量
h=0;
fori=1:
n%对PQ节点的处理
ifi~=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;
fori=1:
n
%六:
对PV节点的处理,注意这时不可再将h初始化为0
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n%对PQ节点的处理
ifB2(i,6)==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==j%对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
else%非对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
end
ifk==(n-1)%将用于内循环的指针置于初始值,以确保雅可比矩阵换行
k=0;
end
end
end
end
end
k=0;
fori=1:
n%对PV节点的处理
ifB2(i,6)==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==j%对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
Jacbi(2*h,2*k)=2*real(B2(i,3));
else%非对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
Jacbi(2*h,2*k-1)=0;
Jacbi(2*h,2*k)=0;
end
ifk==(n-1)%将用于内循环的指针置于初始值,以确保雅可比矩阵换行
k=0;
end
end
end
end
end
Jacbi;
%九:
求解修正方程,获取节点电压的不平衡量
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
DetaU;
%修正节点电压
j=0;
fori=1:
n%对PQ节点处理
ifB2(i,6)==2
j=j+1;
B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
end
end
fori=1:
n%对PV节点的处理
ifB2(i,6)==3
j=j+1;
B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
end
end
B2;
%十:
开始循环**********************************************************************
whileabs(max(DetaU))>pr
OrgS=zeros(2*n-2,1);%!
!
!
初始功率参数在迭代过程中是不累加的,所以在这里必须将其初始化为零矩阵
h=0;
j=0;
fori=1:
n
ifi~=isb&B2(i,6)==2
h=h+1;
forj=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
fori=1:
n
ifi~=isb&B2(i,6)==3
h=h+1;
forj=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;
%创建DetaS
h=0;
fori=1:
n
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n
ifB2(i,6)==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==j
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));
Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));
else
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);
Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);
end
ifk==(n-1)
k=0;
end
end
end
end
end
k=0;
fori=1:
n
ifB2(i,6)==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==j
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));
Jacbi(2*h,2*k-1)=2*imag(B2(i,3));
Jacbi(2*h,2*k)=2*real(B2(i,3));
else
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));
Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));
Jacbi(2*h,2*k-1)=0;
Jacbi(2*h,2*k)=0;
end
ifk==(n-1)
k=0;
end
end
end
end
end
Jacbi;
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
DetaU;
%修正节点电压
j=0;
fori=1:
n
ifB2(i,6)==2
j=j+1;
B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
end
end
fori=1:
n
ifB2(i,6)==3
j=j+1;
B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);
end
end
B2;
Times=Times+1;%迭代次数加1
end
Times;
五、实验流程
六、实验结果
参数输入:
运行结果:
七、实验体会
通过这次实验,让我第一次接触到了MATLAB,并深切体会到了它的强大之处;潮流计算的计算机算法的实现不仅巩固了我的学过的知识,还让我学到一些MATLAB的编程,虽然在实验的过程中出现了很多的错误,但在老师的细心指导下,问题都解决啦;计算机为我们省去了大量的人工计算,希望在以后的学习中能接触到更多的软件,学习到更多的知识。