编程计算潮流.docx
《编程计算潮流.docx》由会员分享，可在线阅读，更多相关《编程计算潮流.docx（20页珍藏版）》请在冰豆网上搜索。
编程计算潮流
课程设计任务书
学生姓名：
专业班级：
指导教师：
工作单位：
题目:
简单闭式网络潮流计算程序设计
初始条件：
节点数：
4支路数：
4计算精度：
0.00010
支路1：
13（0.0200+j0.0800）
支路2：
14（0.0400+j0.1200）
支路3：
24（0.0500+j0.1400）
支路4：
34（0.0400+j0.1200）
节点1：
PQ节点，S
（1）=-0.6000-j0.2500
节点2：
PQ节点，S
（2）=-0.8000-j0.3500
节点3：
PV节点，P（3）=0.4000V（3）=0.9500
节点4：
平衡节点，U（4）=1.00000.0000
要求完成的主要任务:
（包括课程设计工作量及其技术要求，以及说明书撰写等具体要求）
（1）编写计算原理。
（2）编写计算程序。
（3）编写程序并计算结果。
（4）结果分析。
时间安排：
第16周，6月6—10号，指导老师下达任务书，审题、查阅相关资料；
第16周，6月11—16号，分析、计算，编写程序，撰写报告；
第17周，6月17号，论文答辩。
指导教师签名：
年月日
系主任（或责任教师）签名：
年月日
摘要
本文运用matlab软件对给定课设题目进行潮流计算。
潮流计算是电力系统课程最基本最常用的计算。
潮流计算是指对电力系统正常运行状况的分析和计算。
通常需要已知系统参数和条件，给定一些初始条件，从而计算出系统运行的电压和功率等；潮流计算方法很多：
高斯-塞德尔法、牛顿-拉夫逊法、PQ分解法、直流潮流法，以及由高斯-塞德尔法、牛顿-拉夫逊法演变的各种潮流计算方法。
通过潮流计算，可以确定各母线的电压幅值和相角，各元件流过的功率，整个系统的功率损耗。
潮流计算是实现安全经济发供电的必要手段和重要工作环节。
因此潮流计算在电力系统的规划计算，生产运行，调度管理及科学计算中都有广泛的运用。
本课程设计采用PQ分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力系统中各节点电压，为进一步进行电力系统分析作准备。
关键词：
PQ分解法matlab潮流计算
1.题目原始数据及化简
原始数据：
编程计算潮流
节点数：
4支路数：
4计算精度：
0.00010
支路1：
13（0.0200+j0..0800）
支路2：
14（0.0400+j0.1200）
支路3：
24（0.0500+j0.1400）
支路4：
34（0.0400+j0.1200）
节点1：
PQ节点，S
（1）=-0.6000-j0.2500
节点2：
PQ节点，S
（2）=-0.8000-j0.3500
节点3：
PV节点，P（3）=0.4000V（3）=0.9500
节点4：
平衡节点，U（4）=1.00000.0000
根据原始数据所画电路简化图如图1
13
42
图1电路简化图
2.PQ分解法
2.1PQ分解法基本思想
PQ分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：
把节点功率表示为电压向量的极坐标形式，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功和无功分开进行迭代其主要特点是以一个（n-1）阶和一个m阶不变的、对称的系数矩阵代替原来的（n+m-1）阶变化的、不对称的系数矩阵M，以此提高计算速度，降低对计算机贮存容量的要求。
牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时，修正方程式为：
或展开为：
（1）
以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特点。
电力系统中有功功率主要与各节点电压向量的角度有关，无功功率则主要受各节点电压幅值的影响。
大量运算经验也告诉我们，矩阵N及J中各元素的数值相对是很小的，因此对牛顿法的第一步简化就是把有功功率和无功功率分开来进行迭代，即将式（4）化简为：
（2）
这样，由于我们把2n阶的线性方程组变成了二个n阶的线性方程组，因而计算量和内存方面都有改善。
但是，H，L在迭代过程中仍然不断变化，而且又都是不对称矩阵。
对牛顿法的第二个化简，也是比较关键的一个化简，即把式（5）中的系数矩阵简化为在迭代过程中不变的对称矩阵。
众所周知，一般线路两端电压的相角差是不大的（通常不超过10~20度），因此可以认为：
（3）
此外，与系统各节点无功功率相应的导纳必定远远小于该节点自导纳的虚部，即：
因此，（4）
考虑到以上关系后，式（5）中系数矩阵中的元素表达式可以化简为：
（5）
这样，式（5）中系数矩阵可以表示为：
（6）
进一步可以把它们表示为以下矩阵的乘积：
（7）
将它代入（5）中，并利用乘法结合率，我们可以把修正方程式变为：
（8）
及
（9）
将以上两式的左右两侧用以下矩阵左乘
就可以得到
（10）
（11）
以上两式就是P-Q分解法达到修正方程式，其中系数矩阵只不过是系统导纳矩阵的虚部，因而是对称矩阵，而且在迭代过程中维持不变。
它们与功率误差方程式
（12）
（13）
构成了P-Q分解法迭代过程中基本计算公式
2.2PQ分解法潮流计算基本步骤
1）形成系数矩阵，并求其逆矩阵。
2）设各节点电压的初值为（i=1,2，…，n，i≠s）和（i=1,2，…，m，i≠s）。
3）通过有功功率的不平衡方程计算有功功率的不平衡量,从而求出（i=1,2，…，n，i≠s）。
4）解修正方程式，求各节点电压相位角的变量（i=1,2，…，n，i≠s）
5）求各节点电压相位角的新值（i=1,2，…，n，i≠s）。
6）通过无功功率的不平衡方程计算无功功率的不平衡量，从而求出（i=1,2，…，m，i≠s）。
7）解修正方程式，求各节点电压大小的变量（i=1,2，…，m，i≠s）。
8）求各节点电压大小的新值（i=1,2，…，m，i≠s）。
9）运用各节点电压的新值自第三步开始进入下一次迭代。
10）计算平衡节点功率和线路功率
3编程及运行
3.1PQ分解法潮流计算程序框图
是
是
否
是
否
是
否
图2PQ分解法潮流计算流程框图
3.2源程序代码
n=input（'请输入节点数：
n='）;
nl=input（'请输入支路数：
nl='）;
isb=input（'请输入平衡母线节点号：
isb='）;
pr=input（'请输入误差精度：
pr='）;
B1=input（'请输入由支路参数形成的矩阵：
B1='）;%输入B1
B2=input（'请输入由支路参数形成的矩阵：
B2='）;%输入B2
X=input（'请输入由节点号及其对地阻抗形成的矩阵：
X='）;%输入X
na=input（'请输入PQ节点数na='）;
Y=zeros（n）;YI=zeros（n）;e=zeros（1,n）;f=zeros（1,n）;V=zeros（1,n）;O=zeros（1,n）;
fori=1:
n
ifX（i,2）~=0;
p=X（i,1）;
Y（p,p）=1./X（i,2）;
end
end
fori=1:
nl
ifB1（i,6）==0
p=B1（i,1）;q=B1（i,2）;
elsep=B1（i,2）;q=B1（i,1）;
end
Y（p,q）=Y（p,q）-1./（B1（i,3）*B1（i,5））;
YI（p,q）=YI（p,q）-1./B1（i,3）;
Y（q,p）=Y（p,q）;
YI（q,p）=YI（p,q）;
Y（q,q）=Y（q,q）+1./（B1（i,3）*B1（i,5）^2）+B1（i,4）./2;
YI（q,q）=YI（q,q）+1./B1（i,3）;
Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）./2;
YI（p,p）=YI（p,p）+1./B1（i,3）;
end%求导纳矩阵
G=real（Y）;B=imag（YI）;BI=imag（Y）;
fori=1:
n
S（i）=B2（i,1）-B2（i,2）;
BI（i,i）=BI（i,i）+B2（i,5）;
end
P=real（S）;Q=imag（S）;
fori=1:
n
e（i）=real（B2（i,3））;
f（i）=imag（B2（i,3））;
V（i）=B2（i,4）;
end
fori=1:
n
ifB2（i,6）==2
V（i）=sqrt（e（i）^2+f（i）^2）;
O（i）=atan（f（i）./e（i））;
end
end
fori=2:
n
ifi==n
B（i,i）=1./B（i,i）;
elseIC1=i+1;
forj1=IC1:
n
B（i,j1）=B（i,j1）./B（i,i）;
end
B（i,i）=1./B（i,i）;
fork=i+1:
n
forj1=i+1:
n
B（k,j1）=B（k,j1）-B（k,i）*B（i,j1）;
end
end
end
end
p=0;q=0;
fori=1:
n
ifB2（i,6）==2
p=p+1;k=0;
forj1=1:
n
ifB2（j1,6）==2
k=k+1;
A（p,k）=BI（i,j1）;
end
end
end
end
fori=1:
na
ifi==na
A（i,i）=1./A（i,i）;
elsek=i+1;
forj1=k:
na
A（i,j1）=A（i,j1）./A（i,i）;
end
A（i,i）=1./A（i,i）;
fork=i+1:
na
forj1=i+1:
na
A（k,j1）=A（k,j1）-A（k,i）*A（i,j1）;
end
end
end
end
ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;
whileICT2~=0|ICT3~=0
ICT2=0;ICT3=0;
fori=1:
n
ifi~=isb
C（i）=0;
fork=1:
n
C（i）=C（i）+V（k）*（G（i,k）*cos（O（i）-O（k））+BI（i,k）*sin（O（i）-O（k）））;
end
DP1（i）=P（i）-V（i）*C（i）;
DP（i）=DP1（i）./V（i）;
DET=abs（DP1（i））;
ifDET>=pr
ICT2=ICT2+1;
end
end
end
Np（K）=ICT2;
ifICT2~=0
fori=2:
n
DP（i）=B（i,i）*DP（i）;
ifi~=n
IC1=i+1;
fork=IC1:
n
DP（k）=DP（k）-B（k,i）*DP（i）;
end
else
forLZ=3:
i
L=i+3-LZ;
IC4=L-1;
forMZ=2:
IC4
I=IC4+2-MZ;
DP（I）=DP（I）-B（I,L）*DP（L）;
end
end
end
end
fori=2:
n
O（i）=O（i）-DP（i）;
end
kq=1;L=0;
fori=1:
n
ifB2（i,6）==2
C（i）=0;L=L+1;
fork=
升级会员 