电力系统分析潮流计算matlab

1、精选优质文档-倾情为你奉上目录：1、 软件需求说明书.32、 概要设计说明书.41、编写潮流计算程序.4 2、数据的输入测试.4 3、运行得出结果.4 4、进行实验结果验证.43、 详细设计说明书.51、数据导入模块.52、节点导纳矩阵模块.53、编号判断模块.54、收敛条件判定模块.55、雅可比矩阵模块.56、迭代计算模块.57、计算输出参数模块.54、 程序代码.65、 最测试例.151、输入结果.152、输出结果.153、结果验证.151、 软件需求说明书 本次设计利用MATLAB/C+/C（使用MATLAB）编程工具编写潮流计算，实现对节点电压和功率分布的求取。潮流方程的求解基本方法是

2、迭代，包括牛顿-拉夫逊法，以及P-Q分解法，本次设计采用牛顿迭代法。牛顿迭代法（Newton's method）又称为牛顿-拉夫逊方法，它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。 牛顿迭代法是取x0 之后，在这

3、个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知

4、节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。2、 概要设计说明书1、编写潮流计算程序本程序主要分为七个模块：数据导入模块、节点导纳矩阵模块、编号判断模块、收敛条件判定模块、雅可比矩阵模块、迭代计算模块、计算输出参数模块。下图为潮流迭代框图。2、 数据的输入测试 本次设计是将电力系

5、统分析课本P88例题进行潮流计算。3、 运行得出结果 得出电压、有功功率、无功功率、角度的数据。4、进行实验结果验证3、 详细设计说明书 1、数据导入模块： 利用Excel输入已知节点、支路数据，通过“读取”将数据导入MATLAB中。 2、节点导纳矩阵模块： 利用已知的电阻电抗及导纳的值，根据导纳的计算公式，计算出节点的自导纳及节点间的互导纳的值，按照节点编号组成导纳矩阵，利用MATLAB“real”和“imag”调用导纳矩阵中的实部和虚部，分别形成电导和电纳的矩阵。 3、编号判断模块： 当首节点在变压器左侧，设为1，位于变压器右侧，设为2，既非1也非0为不含变压器； 节点类型为PQ时，为1，

6、节点类型为PV时，为2，节点类型为V时，为3。 4、收敛条件判定模块： 根据节点的类型赋初值，并进行失配功率的初步计算，判断是否符合收敛条件maxPi，Qi<。如不符合，则进行后续的计算。 5、雅可比矩阵模块： 根据节点类型确定雅克比矩阵的阶数，然后根据n维非线性方程组的修正方程求出雅克比矩阵。 6、迭代计算模块： 解修正方程，并进行收敛条件判断，如不符合条件则进行下一次迭代，以一直到符合条件为止。 7、计算输出参数模块： 当满足收敛条件maxPi，Qi<时，结束迭代计算，计算电压、有功功率、无功功率、电压角度等参数。4、 程序代码clear %清除变量；filename='

7、;E:MATLABzhilu.xlsx'a=xlsread(filename); %读取支路信息；filename='E:MATLABjiedian.xlsx' c=xlsread(filename); %读取节点信息 b=zeros(a(1,7); %定义节点导纳矩阵G=zeros(a(1,7);B=zeros(a(1,7);for i1=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

8、(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(

9、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(

10、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(

11、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)

12、=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(i

13、1,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); endend %计数各

14、个节点个数PQjd=0;PVjd=0;VOjd=0;for i2=1:a(1,7) if c(i2,6)=1 PQjd=PQjd+1; end if c(i2,6)=2 PVjd=PVjd+1; end if c(i2,6)=3 VOjd=VOjd+1; end end%注入功率的计算zhuru_PQ=zeros(PQjd*2+PVjd,1);k=1;for i3=1:a(1,7) if c(i3,6)=1 zhuru_PQ(k,1)=c(i3,2)-c(i3,4); zhuru_PQ(k+a(1,7)-1,1)=0; for i10=1:a(1,7) zhuru_PQ(k+a(1,7)-1,

15、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 if c(i3,6)=2 zhuru_PQ(k,1)=c(i3,2)-c(i3,4); k=k+1; endendn=a(1,7);%后续矩阵变量定义angle_u=zeros(PQjd*2+PVjd,1); %迭代角度、电压矩阵angle_u(n:PQjd*2+PVjd,1)=1; %电压初值设为1，角度初值设为0delta_au=zeros(PQjd*2+PVjd,1);delta

16、_au(n:PQjd*2+PVjd,1)=0;%角度、电压修正量向量%失配功率的计算 k=1; k_=PQjd+1; delta_PQ= zeros(PQjd*2+PVjd,1); for i4=1:n%PQ节点的计算 if c(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; for i5=1:n if c(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(

17、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 if c(i5,6)=2 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(

18、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 if c(i5,6)=3 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(

19、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节点的计算 if c(i4,6)=2 delta_PQ(k_,1)=zhuru_PQ(k_,1); i6=1; for i5=1:n if c(i5,6)=1 delta_PQ(k_,1)=delta_PQ(k_,1)-c(i4,1)*angle_u

20、(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 if c(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 if c(i5,6)=3 delta_PQ(k_,1)=delta_PQ(k

21、_,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 endzhongjian_delta_PQ=abs(delta_PQ); %取绝对值jingdu=max(zhongjian_delta_PQ); %为迭代精度的判断做赋值准备 while (jingdu>0.0001) %雅各比矩阵的形成 J=zeros(2*PQjd+PVjd); %矩阵定义 U=ones(n,1); %为方便运算将电压取出放入一个新的矩阵 for i=1:PQjd U(i,

22、1)=angle_u(i+n-1,1); end angle=zeros(n,1); %为方便运算将角度取出放入一个新的矩阵 for i=1:n-1 angle(i,1)=angle_u(i,1); end %H% %H部分生成 for i1=1:n-1 for i2=1:n-1 if i1=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 if i1=i2 for i3=1:n if i3=i1 J(i1,i2)=J(i1,i

23、2)+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部分生成 for i1=1:n-1 for i2=1:PQjd if i1=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 if i1=i2 J(i1,i2+n-1)=(-2)*U(i1,1)*U

24、(i1,1)*G(i1,i1); for i3=1:n if i3=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部分生成 for i1=1:PQjd for i2=1:n-1 if i1=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(

25、i1,1)-angle(i2,1); end if i1=i2 for i3=1:n if i3=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部分生成 for i1=1:PQjd for i2=1:PQjd if i1=i2 J(i1+n-1,i2+n-1)=-U(i1,1)*(G(i1,i2)*sin(angle(i1,1)-angle(i2,1

26、)-B(i1,i2)*cos(angle(i1,1)-angle(i2,1); end if i1=i2 J(i1+n-1,i2+n-1)=2*U(i1,1)*U(i1,1)*B(i1,i1); for i3=1:n if i3=i1 J(i1+n-1,i2+n-1)=J(i1+n-1,i2+n-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 delta_au=-inv(J)*delta_PQ; %解修正方程angle_u=ang

27、le_u+delta_au; %形成新的初始值%再次计算失配功率 k=1; k_=PQjd+1; delta_PQ= zeros(PQjd*2+PVjd,1); for i4=1:n% if c(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; for i5=1:n if c(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

28、)+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 if c(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)-angl

29、e_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 if c(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

30、,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 % if c(i4,6)=2 delta_PQ(k_,1)=zhuru_PQ(k_,1); i6=1; for i5=1:n if c(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(i