电力系统分析潮流计算matlab.docVIP

目录：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）编程工具编写潮流计算，实现对节点电压和功率分布的求取。潮流方程的求解基本方法是迭代，包括牛顿-拉夫逊法，以及P-Q分解法，本次设计采用牛顿迭代法。牛顿迭代法（Newtonsmethod）又称为牛顿-拉夫逊方法，它是牛顿在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、收敛条件判定模块： 根据节点的类型赋初值，并进行失配功率的初步计算，判断是否符合收敛条件maxPi，Qi。如不符合，则进行后续的计算。 5、雅可比矩阵模块： 根据节点类型确定雅克比矩阵的阶数，然后根据n维非线性方程组的修正方程求出雅克比矩阵。 6、迭代计算模块： 解修正方程，并进行收敛条件判断，如不符合条件则进行下一次迭代，以一直到符合条件为止。 7、计算输出参数模块： 当满足收敛条件maxPi，Qi0.0001) %雅各比矩阵的形成 J=zeros(2*PQjd+PVjd); %矩阵定义 U=ones(n,1); %为方便运算将电压取出放入一个新的矩阵 for i=1:PQjd U(i,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,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部分生成 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(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(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)-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=angle_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)+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)-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(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 % 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(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_,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 %再次进行精度判断前的准备 m=0; zhongjian_delta_PQ=delta_PQ;zhongjian_delta_PQ=abs(delta_PQ);jingdu=max(zhongjian_delta_PQ);m=m+1; end %迭代完毕 %电压与角度换算结果输出 U=ones(n,1); for i=1:PQjd U(i,1)=angle_u(i+n-1,1); end angle=zeros(n,1); for i=1:n-1 angle(i,1)=angle_u(i,1); endangle_jiaodu=angle/(2*pi)*360;%计算各个节点功率 P=zeros(n,1); Q=zeros(n,1); for i1=1:n for i2=1:n P(i1,1)=P(i1,1)+U(i1,1)*U(i2,1)*(G(i1,i2)*cos(angle(i1,1)-angle(i