电力系统分析潮流计算matlab.docx
《电力系统分析潮流计算matlab.docx》由会员分享，可在线阅读，更多相关《电力系统分析潮流计算matlab.docx（20页珍藏版）》请在冰点文库上搜索。
电力系统分析潮流计算matlab
目录：
1、软件需求说明书...................................................................................3
2、概要设计说明书...................................................................................4
1、编写潮流计算程序.........................................................................4
2、数据的输入测试.............................................................................4
3、运行得出结果.................................................................................4
4、进行实验结果验证.........................................................................4
3、详细设计说明书...................................................................................5
1、数据导入模块.................................................................................5
2、节点导纳矩阵模块.........................................................................5
3、编号判断模块.................................................................................5
4、收敛条件判定模块..........................................................................5
5、雅可比矩阵模块..............................................................................5
6、迭代计算模块..................................................................................5
7、计算输出参数模块..........................................................................5
4、程序代码...............................................................................................6
5、最测试例.............................................................................................15
1、输入结果.......................................................................................15
2、输出结果.......................................................................................15
3、结果验证.......................................................................................15
1、软件需求说明书
本次设计利用MATLAB/C++/C（使用MATLAB）编程工具编写潮流计算，实现对节点电压和功率分布的求取。
潮流方程的求解基本方法是迭代，包括牛顿-拉夫逊法，以及P-Q分解法，本次设计采用牛顿迭代法。
牛顿迭代法（Newton's method）又称为牛顿-拉夫逊方法，它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。
多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。
方法使用函数f（x）的泰勒级数的前面几项来寻找方程f（x） = 0的根。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f（x） = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
牛顿迭代法是取x0之后，在这个基础上，找到比x0更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f（x）=0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。
为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
2、概要设计说明书
1、编写潮流计算程序
本程序主要分为七个模块：
数据导入模块、节点导纳矩阵模块、编号判断模块、收敛条件判定模块、雅可比矩阵模块、迭代计算模块、计算输出参数模块。
下图为潮流迭代框图。
2、数据的输入测试
本次设计是将《电力系统分析》课本P88例题进行潮流计算。
3、运行得出结果
得出电压、有功功率、无功功率、角度的数据。
4、进行实验结果验证
3、详细设计说明书
1、数据导入模块：
利用Excel输入已知节点、支路数据，通过“读取”将数据导入MATLAB中。
2、节点导纳矩阵模块：
利用已知的电阻电抗及导纳的值，根据导纳的计算公式，计算出节点的自导纳及节点间的互导纳的值，按照节点编号组成导纳矩阵，利用MATLAB“real”和“imag”调用导纳矩阵中的实部和虚部，分别形成电导和电纳的矩阵。
3、编号判断模块：
当首节点在变压器左侧，设为1，位于变压器右侧，设为2，既非1也非0为不含变压器；
节点类型为PQ时，为1，节点类型为PV时，为2，节点类型为VƟ时，为3。
4、收敛条件判定模块：
根据节点的类型赋初值，并进行失配功率的初步计算，判断是否符合收敛条件max{ΔPi，ΔQi}<
。
如不符合，则进行后续的计算。
5、雅可比矩阵模块：
根据节点类型确定雅克比矩阵的阶数，然后根据n维非线性方程组的修正方程求出雅克比矩阵。
6、迭代计算模块：
解修正方程，并进行收敛条件判断，如不符合条件则进行下一次迭代，以一直到符合条件为止。
7、计算输出参数模块：
当满足收敛条件max{ΔPi，ΔQi}<
时，结束迭代计算，计算电压、有功功率、无功功率、电压角度等参数。
4、程序代码
clear%清除变量；
filename='E:
\MATLAB\zhilu.xlsx';
a=xlsread（filename）;%读取支路信息；
filename='E:
\MATLAB\jiedian.xlsx';
c=xlsread（filename）;%读取节点信息
b=zeros（a（1,7））;%定义节点导纳矩阵
G=zeros（a（1,7））;
B=zeros（a（1,7））;
fori1=1:
（a（1,7））
if（a（i1,5）==0）%%%首节点在变压器左侧
b（a（i1,1）,a（i1,1））=b（a（i1,1）,a（i1,1））+1/（a（i1,3）+1j*a（i1,4））;
b（a（i1,2）,a（i1,2））=b（a（i1,2）,a（i1,2））+1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）^2）;
b（a（i1,1）,a（i1,2））=b（a（i1,1）,a（i1,2））+（-1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）））;
b（a（i1,2）,a（i1,1））=b（a（i1,2）,a（i1,1））+（-1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）））;%%%进行导纳计算
G（a（i1,1）,a（i1,1））=real（b（a（i1,1）,a（i1,1）））;
G（a（i1,2）,a（i1,2））=real（b（a（i1,2）,a（i1,2）））;
G（a（i1,1）,a（i1,2））=real（b（a（i1,1）,a（i1,2）））;
G（a（i1,2）,a（i1,1））=real（b（a（i1,2）,a（i1,1）））;
B（a（i1,1）,a（i1,1））=imag（b（a（i1,1）,a（i1,1）））;
B（a（i1,2）,a（i1,2））=imag（b（a（i1,2）,a（i1,2）））;
B（a（i1,1）,a（i1,2））=imag（b（a（i1,1）,a（i1,2）））;
B（a（i1,2）,a（i1,1））=imag（b（a（i1,2）,a（i1,1）））;
elseif（a（i1,5）==1）%%%首节点在变压器右侧
b（a（i1,1）,a（i1,1））=b（a（i1,1）,a（i1,1））+1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）^2）;
b（a（i1,2）,a（i1,2））=b（a（i1,2）,a（i1,2））+1/（a（i1,3）+1j*a（i1,4））;
b（a（i1,1）,a（i1,2））=b（a（i1,1）,a（i1,2））+（-1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）））;
b（a（i1,2）,a（i1,1））=b（a（i1,2）,a（i1,1））+（-1/（（a（i1,3）+1j*a（i1,4））*a（i1,6）））;%%%进行导纳计算
G（a（i1,1）,a（i1,1））=real（b（a（i1,1）,a（i1,1）））;
G（a（i1,2）,a（i1,2））=real（b（a（i1,2）,a（i1,2）））;
G（a（i1,1）,a（i1,2））=real（b（a（i1,1）,a（i1,2）））;
G（a（i1,2）,a（i1,1））=real（b（a（i1,2）,a（i1,1）））;
B（a（i1,1）,a（i1,1））=imag（b（a（i1,1）,a（i1,1）））;
B（a（i1,2）,a（i1,2））=imag（b（a（i1,2）,a（i1,2）））;
B（a（i1,1）,a（i1,2））=imag（b（a（i1,1）,a（i1,2）））;
B（a（i1,2）,a（i1,1））=imag（b（a（i1,2）,a（i1,1）））;
else%%%不含变压器支路
b（a（i1,1）,a（i1,1））=b（a（i1,1）,a（i1,1））+1/（a（i1,3）+1j*a（i1,4））+1j*a（i1,5）;
b（a（i1,2）,a（i1,2））=b（a（i1,2）,a（i1,2））+1/（a（i1,3）+1j*a（i1,4））+1j*a（i1,5）;
b（a（i1,1）,a（i1,2））=b（a（i1,1）,a（i1,2））+（-1/（（a（i1,3）+1j*a（i1,4））））;
b（a（i1,2）,a（i1,1））=b（a（i1,2）,a（i1,1））+（-1/（（a（i1,3）+1j*a（i1,4））））;%%%进行导纳计算
G（a（i1,1）,a（i1,1））=real（b（a（i1,1）,a（i1,1）））;
G（a（i1,2）,a（i1,2））=real（b（a（i1,2）,a（i1,2）））;
G（a（i1,1）,a（i1,2））=real（b（a（i1,1）,a（i1,2）））;
G（a（i1,2）,a（i1,1））=real（b（a（i1,2）,a（i1,1）））;
B（a（i1,1）,a（i1,1））=imag（b（a（i1,1）,a（i1,1）））;
B（a（i1,2）,a（i1,2））=imag（b（a（i1,2）,a（i1,2）））;
B（a（i1,1）,a（i1,2））=imag（b（a（i1,1）,a（i1,2）））;
B（a（i1,2）,a（i1,1））=imag（b（a（i1,2）,a（i1,1）））;
end
end
%%%%计数各个节点个数
PQjd=0;
PVjd=0;
VOjd=0;
fori2=1:
a（1,7）
ifc（i2,6）==1
PQjd=PQjd+1;
end
ifc（i2,6）==2
PVjd=PVjd+1;
end
ifc（i2,6）==3
VOjd=VOjd+1;
end
end
%注入功率的计算
zhuru_PQ=zeros（PQjd*2+PVjd,1）;
k=1;
fori3=1:
a（1,7）
ifc（i3,6）==1
zhuru_PQ（k,1）=c（i3,2）-c（i3,4）;
zhuru_PQ（k+a（1,7）-1,1）=0;
fori10=1:
a（1,7）
zhuru_PQ（k+a（1,7）-1,1）=zhuru_PQ（k+a（1,7）-1,1）+B（i3,i10）;
end
zhuru_PQ（k+a（1,7）-1,1）=zhuru_PQ（k+a（1,7）-1,1）-c（i3,5）
k=k+1;
end
ifc（i3,6）==2
zhuru_PQ（k,1）=c（i3,2）-c（i3,4）;
k=k+1;
end
end
n=a（1,7）;
%%%%%后续矩阵变量定义
angle_u=zeros（PQjd*2+PVjd,1）;%迭代角度、电压矩阵
angle_u（n:
PQjd*2+PVjd,1）=1;%电压初值设为1，角度初值设为0
delta_au=zeros（PQjd*2+PVjd,1）;
delta_au（n:
PQjd*2+PVjd,1）=0;
%角度、电压修正量向量
%%%%%失配功率的计算
k=1;
k_=PQjd+1;
delta_PQ=zeros（PQjd*2+PVjd,1）;
fori4=1:
n
%%%%%%%%%%%%%%%%%%%PQ节点的计算
ifc（i4,6）==1
delta_PQ（k,1）=zhuru_PQ（k,1）;
delta_PQ（k+n-1,1）=zhuru_PQ（k+n-1,1）;
i6=1;
fori5=1:
n
ifc（i5,6）==1
delta_PQ（k,1）=delta_PQ（k,1）-angle_u（n-1+k）*angle_u（n-1+i6）*（G（i4,i5）*cos（angle_u（k,1）-angle_u（i6,1））+B（i4,i5）*sin（angle_u（k,1）-angle_u（i6,1）））;
delta_PQ（k+n-1,1）=delta_PQ（k+n-1,1）-angle_u（n-1+k）*angle_u（n-1+i6）*（G（i4,i5）*sin（angle_u（k,1）-angle_u（i6,1））-B（i4,i5）*cos（angle_u（k,1）-angle_u（i6,1）））;
i6=i6+1;
end
ifc（i5,6）==2
delta_PQ（k,1）=delta_PQ（k,1）-angle_u（n-1+k）*c（i5,1）*（G（i4,i5）*cos（angle_u（k,1）-angle_u（i6,1））+B（i4,i5）*sin（angle_u（k,1）-angle_u（i6,1）））;
delta_PQ（k+n-1,1）=delta_PQ（k+n-1,1）-angle_u（n-1+k）*c（i5,1）*（G（i4,i5）*sin（angle_u（k,1）-angle_u（i6,1））-B（i4,i5）*cos（angle_u（k,1）-angle_u（i6,1）））;
i6=i6+1;
end
ifc（i5,6）==3
delta_PQ（k,1）=delta_PQ（k,1）-angle_u（n-1+k）*c（i5,1）*（G（i4,i5）*cos（angle_u（k,1））+B（i4,i5）*sin（angle_u（k,1）））;
delta_PQ（k+n-1,1）=delta_PQ（k+n-1,1）-angle_u（n-1+k）*c（i5,1）*（G（i4,i5）*sin（angle_u（k,1））-B（i4,i5）*cos（angle_u（k,1）））;
end
end
k=k+1;
end
%%%%%%%%%%%%%%%%%%%%%PV节点的计算
ifc（i4,6）==2
delta_PQ（k_,1）=zhuru_PQ（k_,1）;
i6=1;
fori5=1:
n
ifc（i5,6）==1
delta_PQ（k_,1）=delta_PQ（k_,1）-c（i4,1）*angle_u（n-1+i6）*（G（i4,i5）*cos（angle_u（k_,1）-angle_u（i6,1））+B（i4,i5）*sin（angle_u（k_,1）-angle_u（i6,1）））;
i6=i6+1;
end
ifc（i5,6）==2
delta_PQ（k_,1）=delta_PQ（k_,1）-c（i4,1）*c（i5,1）*（G（i4,i5）*cos（angle_u（k_,1）-angle_u（i6,1））+B（i4,i5）*sin（angle_u（k_,1）-angle_u（i6,1）））;
i6=i6+1;
end
ifc（i5,6）==3
delta_PQ（k_,1）=delta_PQ（k_,1）-c（i4,1）*c（i5,1）*（G（i4,i5）*cos（angle_u（k_,1））+B（i4,i5）*sin（angle_u（k_,1）））;
end
end
k_=k_+1;
end
end
zhongjian_delta_PQ=abs（delta_PQ）;%%%%取绝对值
jingdu=max（zhongjian_delta_PQ）;%%%%为迭代精度的判断做赋值准备
while（jingdu>0.0001）
%%%%%%%雅各比矩阵的形成
J=zeros（2*PQjd+PVjd）;%%%%矩阵定义
U=ones（n,1）;%%%%为方便运算将电压取出放入一个新的矩阵
fori=1:
PQjd
U（i,1）=angle_u（i+n-1,1）;
end
angle=zeros（n,1）;%%%%为方便运算将角度取出放入一个新的矩阵
fori=1:
n-1
angle（i,1）=angle_u（i,1）;
end
%%%%%H%%%%%%%%%%H部分生成
fori1=1:
n-1
fori2=1:
n-1
ifi1~=i2
J（i1,i2）=-U（i1,1）*U（i2,1）*（G（i1,i2）*sin（angle（i1,1）-angle（i2,1））-B（i1,i2）*cos（angle（i1,1）-angle（i2,1）））;
end
ifi1==i2
fori3=1:
n
ifi3~=i1
J（i1,i2）=J（i1,i2）+U（i1,1）*U（i3,1）*（G（i1,i3）*sin（angle（i1,1）-angle（i3,1））-B（i1,i3）*cos（angle（i1,1）-angle（i3,1）））;
end
end
end
end
end
%%%%%N%%%%%%%%%%N部分生成
fori1=1:
n-1
fori2=1:
PQjd
ifi1~=i2
J（i1,i2+n-1）=-U（i1,1）*（G（i1,i2）*cos（angle（i1,1）-angle（i2,1））-B（i1,i2）*sin（angle（i1,1）-angle（i2,1）））;
end
ifi1==i2
J（i1,i2+n-1）=（-2）*U（i1,1）*U（i1,1）*G（i1,i1）;
fori3=1:
n
ifi3~=i1
J（i1,i2+n-1）=J（i1,i2+n-1）-U（i3,1）*（G（i1,i3）*cos（angle（i1,1）-angle（i3,1））+B（i1,i3）*sin（angle（i1,1）-angle（i3,1）））;
end
end
end
end
end
%%%%%K%%%%%%%%%%K部分生成
fori1=1:
PQjd
fori2=1:
n-1
ifi1~=i2
J（i1+n-1,i2）=U（i1,1）*U（i2,1）*（G（i1,i2）*cos（angle（i1,1）-angle（i2,1））-B（i1,i2）*sin（angle（i1,1）-angle（i2,1）））;
end
ifi1==i2
fori3=1:
n
ifi3~=i1
J（i1+n-1,i2）=J（i1+n-1,i2）-U（i1,1）*U（i3,1）*（G（i1,i3）*cos（angle（i1,1）-angle（i3,1））+B（i1,i3）*sin（angle（i1,1）-angle（i3,1）））;
end
end
end
end
end
%%%%%L%%%%%%%%%%L部分生成
fori1=1:
PQjd
fori2=1:
PQjd
ifi1~=i2
J（i1+n-1,i2+n-1）=-U（i1,1）*（G（i1,i2）*sin（angle（i1,1）-angle（i2,1））-B（i1,i2）*cos（angle（i1,1）-angle（i2,1）））;
end
ifi1==i2
J（i1+n-1,i2+n-1）=2*U（
升级会员 