电力系统分析潮流计算的计算机算法.docx
《电力系统分析潮流计算的计算机算法.docx》由会员分享，可在线阅读，更多相关《电力系统分析潮流计算的计算机算法.docx（16页珍藏版）》请在冰豆网上搜索。
电力系统分析潮流计算的计算机算法
潮流计算的计算机算法实验报告
：
学号：
班级：
一、实验目的
掌握潮流计算的计算机算法。
熟悉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的编程，虽然在实验的过程中出现了很多的错误，但在老师的细心指导下，问题都解决啦；计算机为我们省去了大量的人工计算，希望在以后的学习中能接触到更多的软件，学习到更多的知识。
升级会员 