C语言报告-潮流.docx

节点电压用直角坐标表示时的牛顿拉夫逊法潮流计算的C语言实现一、实验目的 根据所给的电力系统，编制牛顿拉夫逊法潮流计算程序，最后完成一个切实可行的电力系统计算应用程序。通过实验的过程加深对于原课程内容的掌握，并了解将其实现的手段。二、原理介绍采用直角坐标时，节点电压可表示为 导纳矩阵元素则表示为将上述表示式带入， i=1,2,3.,m并展开分出实部和虚部便可得 i=1,2,3.,n假定系统中的第1,2.，m号节点为PQ节点，第个节点的给定功率设为和，对该节点可列写方程i=1,2,3.,m假定系统中的第m+1,m+2,.,n-1号节点为PV节点，则对其中每一个节点可以列写方程i= m+1,m+2,.,n-1第n号节点为平衡节点，其电压是给定的，故不参加迭代由上不难写出如下的修正方程式式中上述方程中雅克比矩阵的各元素，可以求偏导数获得，当时当时修正方程式还可以写成分块矩阵的形式其中对于PQ节点对于PV节点用牛顿拉夫逊法计算潮流，首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值和，置迭代计数。然后开始进入牛顿法的迭代过程。在进行第次迭代时，其计算步骤如下：（1） 按上一次迭代算出的节点电压值和（当时即为给定的初值），利用公式计算各类节点的不平衡量、和。（2） 按条件校验收敛，即如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。（3） 利用公式计算雅克比矩阵的个元素。（4） 解修正方程式求节点电压的修正量和。（5） 修正各节点的电压（6） 迭代计数加1，返回第一步继续迭代过程。迭代结束后，还要算出平衡节点的功率和网络中的功率分布。输电线路功率的计算公式如下三、源代码与注释/ Y-matrix.cpp : 定义控制台应用程序的入口点。/#include stdafx.h#include #define array_size 10int _tmain(int argc, _TCHAR* argv)int node_num;int branch_num;float rxbarray_size5;FILE *fp;/读入数据fp = fopen(branchdata.txt,r);fscanf(fp,%d%d,&node_num,&branch_num);int i = 0;int j = 0;int k = 0;for (i = 0;i = branch_num-1; i+)for (j = 0;j = 4; j+)fscanf(fp,%f,&rxbij);fclose(fp);float ratioarray_size;/分离变比for (i=0;i=branch_num-1;i+)if (rxbi4=0)ratioi = 1;rxbi4 = -rxbi4;elseratioi = rxbi4;rxbi4 = 0;float y_rearray_size;float y_imarray_size;float Ymatrix_rearray_sizearray_size;float Ymatrix_imarray_sizearray_size;for (i=0;iarray_size;i+)y_rei=0;y_imi=0;for (j=0;jarray_size;j+)Ymatrix_reij=0;Ymatrix_imij=0;for (i=0;ibranch_num;i+)/计算节点间导纳y_rei=rxbi2/(rxbi2*rxbi2+rxbi3*rxbi3);y_imi=-rxbi3/(rxbi2*rxbi2+rxbi3*rxbi3);for (i=0;ibranch_num;i+)/导纳阵非对角线元素Ymatrix_re(int)rxbi0-1(int)rxbi1-1=-y_rei/ratioi;Ymatrix_re(int)rxbi1-1(int)rxbi0-1=Ymatrix_re(int)rxbi0-1(int)rxbi1-1;Ymatrix_im(int)rxbi0-1(int)rxbi1-1=-y_imi/ratioi;Ymatrix_im(int)rxbi1-1(int)rxbi0-1=Ymatrix_im(int)rxbi0-1(int)rxbi1-1;for (i=0;inode_num;i+)/导纳阵对角线元素for (j=0;jbranch_num;j+)if (int)rxbj1=(i+1)Ymatrix_reii=Ymatrix_reii+y_rej/(ratioj*ratioj);Ymatrix_imii=Ymatrix_imii+y_imj/(ratioj*ratioj)+rxbj4; else if (int)rxbj0=(i+1)Ymatrix_reii=Ymatrix_reii+y_rej;Ymatrix_imii=Ymatrix_imii+y_imj+rxbj4;float nodedataarray_size4;/ 0 - 节点号1-节点类型(0pq 1 pv 2 平衡) 2 3 -节点电压初始值int pq_num = 0;int pv_num = 0;int slack_num = 0;fp = fopen(nodedata.txt,r);for (i = 0;i node_num; i+)for (j = 0;j 4; j+)fscanf(fp,%f,&nodedataij);if (nodedatai1 = 0)pq_num+;else if (nodedatai1=1)pv_num+;elseslack_num+;fclose(fp);/设定变量，设初值float Ps50;float Qs50;float Vs50;float max_error = pow(10.0,-5);Ps0 = -0.30;Qs0 = -0.18;Ps1 = -0.55;Qs1 = -0.13;Ps2= 0.5;Vs2 = 1.10;float P_delta50;float Q_delta50;float V2_delta50;float delta_max = 0;float e50;float e_delta50;float f50;float f_delta50;float Jacobi5050;float invJacobi5050;float Jacobi_temp5050;float L5050;float R5050; float X5050;float P_slack = 0;float Q_slack = 0;for (i=0;i50;i+)P_deltai=0;Q_deltai=0;V2_deltai=0;ei=0;e_deltai=0;fi=0;f_deltai=0;for (j=0;j50;j+)Jacobiij=0;invJacobiij=0;Lij=0;Rij=0; Xij=0;for (i = 0;i0 )delta_max = 0;for (i=0;ipq_num;i+)/计算deltaP PQP_deltai = Psi;for (j = 0;jnode_num;j+)P_deltai = P_deltai -ei*(Ymatrix_reij*ej-Ymatrix_imij*fj);for (j = 0;j delta_max) delta_max = abs(P_deltai);for (i=0;ipq_num;i+)/计算deltaQ PQQ_deltai = Qsi;for (j = 0;jnode_num;j+)Q_deltai = Q_deltai -fi*(Ymatrix_reij*ej-Ymatrix_imij*fj);for (j = 0;j delta_max)delta_max = abs(Q_deltai);for (i = pq_num;ipq_num+pv_num;i+)/计算deltaP PVP_deltai = Psi;for (j = 0;jnode_num;j+)P_deltai = P_deltai -ei*(Ymatrix_reij*ej-Ymatrix_imij*fj);for (j = 0;j delta_max)delta_max = abs(P_deltai);for (i = pq_num;i delta_max)delta_max = abs(V2_deltai);if (delta_max max_error)/判断是否收敛calFlag = -1;break;/雅克比矩阵Pfor (i =0;ipq_num+pv_num;i+)for (j=0;jpq_num+pv_num;j+)if (i = j)Jacobi2*i+12*j = -Ymatrix_reii*ei-Ymatrix_imii*fi;/p efor (k = 0;knode_num;k+)Jacobi2*i+12*j = Jacobi2*i+12*j - (Ymatrix_reik*ek - Ymatrix_imik*fk);Jacobi2*i+12*j+1 = Ymatrix_imii*ei-Ymatrix_reii*fi;/p ffor (k = 0;knode_num;k+)Jacobi2*i+12*j+1 = Jacobi2*i+12*j+1 - (Ymatrix_reik*fk + Ymatrix_imik*ek);elseJacobi2*i+12*j = -(Ymatrix_reij*ei+Ymatrix_imij*fi);/p eJacobi2*i+12*j+1 = Ymatrix_imij*ei - Ymatrix_reij*fi; /p f/雅克比矩阵pq Qfor (i=0;ipq_num;i+)for (j=0;jpq_num+pv_num;j+)if (i = j)Jacobi2*i2*j = Ymatrix_imii*ei-Ymatrix_reii*fi;/q efor (k = 0;knode_num;k+)Jacobi2*i2*j = Jacobi2*i2*j + (Ymatrix_reik*fk + Ymatrix_imik*ek);Jacobi2*i2*j+1 = Ymatrix_reii*ei+Ymatrix_imii*fi;/q ffor (k = 0;knode_num;k+)Jacobi2*i2*j+1 = Jacobi2*i2*j+1 - (Ymatrix_reik*ek - Ymatrix_imik*fk); elseJacobi2*i2*j = Ymatrix_imij*ei - Ymatrix_reij*fi;/q eJacobi2*i2*j+1 = Ymatrix_reij*ei+Ymatrix_imij*fi;/ 雅克比矩阵pv vfor (i = pq_num;ipq_num+pv_num;i+)for (j=0;jpq_num+pv_num;j+)if (i =j)Jacobi2*i2*i = -2*ei;Jacobi2*i2*i+1 = -2*fi; elseJacobi2*i2*i = 0;/Jacobi2*i2*i+1 = 0;/求雅克比矩阵的逆for (i = 0; i (node_num-1)*2;i+ )for (j = 0; j (node_num-1)*2;j+ )if (i = j)invJacobiij = 1;elseinvJacobiij = 0;for (i = 0; i 50 ;i+ )for (j = 0; j 50;j+ )Jacobi_tempij = Jacobiij;/高斯消元法float dia_element_temp = 0;float row_element0_temp = 0;for (i = 0; i (node_num-1)*2;i+)dia_element_temp = Jacobi_tempii;for (j = 0; j(node_num-1)*2;j+)Jacobi_tempij = Jacobi_tempij/dia_element_temp;invJacobiij = invJacobiij/dia_element_temp;for (j = 0; j(node_num-1)*2;j+)if (i != j)row_element0_temp = Jacobi_tempji;for (k = 0; k(node_num-1)*2;k+)Jacobi_tempjk = Jacobi_tempjk - Jacobi_tempik*row_element0_temp;invJacobijk = invJacobijk - invJacobiik*row_element0_temp;for (i = 0; i (node_num-1);i+)e_deltai = 0;f_deltai = 0;for (j = 0; j pq_num ;j+)e_deltai = e_deltai+invJacobi2*i2*j+1*P_deltaj+invJacobi2*i2*j*Q_deltaj;f_deltai = f_deltai + invJacobi2*i+12*j+1*P_deltaj + invJacobi2*i+12*j*Q_deltaj;for (j = pq_num; j pq_num+pv_num ;j+)e_deltai = e_deltai + invJacobi2*i2*j+1*P_deltaj + invJacobi2*i2*j*V2_deltaj;f_deltai = f_deltai + invJacobi2*i+12*j+1*P_deltaj + invJacobi2*i+12*j*V2_deltaj;e_deltai = -e_deltai;f_deltai = -f_deltai;for (i = 0; i = (node_num-1);i+)ei = ei + e_deltai;fi = fi + f_deltai;calFlag-;caltime+;for (i = 0; i node_num;i+)P_slack = P_slack +e