Matlab及C语言在潮流计算运用.doc

matlab及c语言在潮流计算中的运用（三峡大学电气信息学院20071096班）摘要：潮流计算是研究电力系统稳态运行情况的一种基本电气计算，常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。关键词：电力系统分析 潮流计算 牛顿-拉夫逊法 c语言 matlab一，潮流计算算法原理：牛顿拉夫逊法的基本原理 牛顿-拉夫逊法是一种求解非线性方程的数值解法，由于便于编写程序用计算机求解，应用较广。下面以一元非线性代数方程的求解为例，来说明牛顿-拉夫逊法的基本思想。 设欲求解的非线性代数方程为 f(x)=o设方程的真实解为x*，则必有f(x*)=0。用牛顿-拉夫逊法求方程真实解x*的步骤如下： 首先选取余割合适的初始估值x作为方程f(x)=0的解，若恰巧有f(x)=0，则方程的真实解即为x*= x若f(x)0，则做下一步。 取x=x+x为第一次的修正估值，则 f(x)=f(x+x)其中x为初始估值的增量，即x=x-x。设函数f(x)具有任意阶导数，即可将上式在x的邻域展开为泰勒级数，即： f(x)=f(x+x)=f(x)+f(x)x+f(x)(x)2/2+若所取的|x|足够小，则含(x)的项及其余的一切高阶项均可略去，并使其等于零，即：f(x)f(x)+f(x)x=0故得 x=-f(x)/f(x)从而 x= x-f(x)/f(x) 可见，只要f(x)0，即可根据上式求出第一次的修正估值x，若恰巧有f(x)=0，则方程的真实解即为x*=x。若f(x)0，则用上述方法由x再确定第二次的修正估值x。如此反复叠代下去，直到求得真实解x*为止。二，节点电压用极坐标牛顿拉夫逊法潮流计算 节点的功率方程写成 其中式中，是两节点电压的相角差。把节点功率表示为节点电压的幅值和相角的函数。在有n个节点的系统中，假定第1m号节点为节点，第m+1n-1号节点为节点，第n号节点为平衡节点。un和n是给定的，pv节点的电压幅值um+1un-1也是给定的。因此，只剩下n-1个节点的电压相角1，2，n-1和m个节点的电压幅值u1，u2，um是未知量。实际上，对于每一个pq节点或每一个pv节点都可以列写一个有功功率不平衡量方程式，一共包含了n-1+m个方程式，正好同未知数的数目相同，而比直角坐标形式的方程少了n-1-m个。对于方程式可以写出修正方程式如下：h是(n-1)(n-1)阶方阵，其元素为 ；n是(n-1)m阶矩阵，其元素为 ；k是m(n-1)阶矩阵，其元素为 ；l是mm阶方阵，其元素为 。在这里把节点不平衡功率对节点电压幅值的偏导数都乘以该节点电压，相应地把节点电压的修正量都除以该节点的电压幅值，这样，雅可比矩阵的表达式就具有比较整齐的形式。三，matlab及c语言运算过程的具体实例例：网络接线如图，各支路阻抗和各节点功率均以标幺值标于图中，其中节点2连接的实际是发额定功率的发电厂，设节点1 的电压保持为1.06，用牛顿-拉夫逊法计算系统中的潮流分布。0.02+j0.060.45+j0.150.08+j0.240.01+j0.030.4+j0.050.08+j0.240.04+j0.120.06+j0.180.06+j0.180.2+j0.2gg52134-0.02+j0.16四，程序运行结果：matlab程序运行结果如下：y = 3.7500 -11.2500i -2.5000 + 7.5000i 0 -1.2500 + 3.7500i 0 -2.5000 + 7.5000i 10.8333 -32.5000i -1.6667 + 5.0000i -1.6667 + 5.0000i -5.0000 +15.0000i 0 -1.6667 + 5.0000i 12.9167 -38.7500i -10.0000 +30.0000i -1.2500 + 3.7500i -1.2500 + 3.7500i -1.6667 + 5.0000i -10.0000 +30.0000i 12.9167 -38.7500i 0 0 -5.0000 +15.0000i -1.2500 + 3.7500i 0 6.2500 -18.7500icount = 4u = 1.0555 1.0572 1.0319 1.0331 1.0600angle = 0.0417 0.0055 -0.0261 -0.0236 0error_angle_u = 1.0e-008 * -0.5980 -0.0366 0.1334 0.1379 -0.8207 -0.4688 -0.5584 -0.6216s = 0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 0.0677 + 0.1531iline_s = 0 0.2905 - 0.1333i 0 0.2991 - 0.0163i 0 -0.2956 + 0.1229i 0 0.2185 + 0.0764i 0.2028 + 0.0746i 0.0764 - 0.0760i 0 -0.2109 - 0.0812i 0 -0.0918 - 0.0132i -0.1389 - 0.0789i -0.2931 - 0.0031i -0.1960 - 0.0786i 0.0919 + 0.0130i 0 0 0 -0.0771 + 0.0758i 0.1448 + 0.0773i 0 0 line_i = 0 0.2697 + 0.1377i 0 0.2824 + 0.0273i 0 -0.2790 - 0.1178i 0 0.2071 - 0.0711i 0.1922 - 0.0695i 0.0719 + 0.0723i 0 -0.2022 + 0.0840i 0 -0.0886 + 0.0151i -0.1326 + 0.0799i -0.2836 + 0.0097i -0.1879 + 0.0806i 0.0886 - 0.0146i 0 0 0 -0.0727 - 0.0715i 0.1366 - 0.0729i 0 0 c语言运行结果如下：节点导纳矩阵：3.750000+j-11.250000 -2.500000+j7.500000 -0.000000+j-0.000000 -1.250000+j3.750000 -0.000000+j-0.000000-2.500000+j7.500000 10.833333+j-32.500000 -1.666667+j5.000000 -1.666667+j5.000000 -5.000000+j15.000000-0.000000+j-0.000000 -1.666667+j5.000000 12.916667+j-38.750000 -10.000000+j30.000000 -1.250000+j3.750000-1.250000+j3.750000 -1.666667+j5.000000 -10.000000+j30.000000 12.916667+j-38.750000 -0.000000+j-0.000000-0.000000+j-0.000000 -5.000000+j15.000000 -1.250000+j3.750000 -0.000000+j-0.000000 6.250000+j-18.750000迭代次数：count=4节点相角分别为:=0.041710 0.005480 -0.026114 -0.023594 0.000000节点电压分别为:u=1.055525 1.057178 1.031931 1.033072 1.060000节点功率为:s=0.600000+j-0.100000 0.200000+j0.200000 -0.450000+j-0.150000 -0.400000+j-0.050000 0.067700+j0.153099请按任意键继续. . .经过比较matlab和c语言运行的结果是一样的！最终的潮流分布如下图：0.0764 - 0.0760i0.45+j0.15-0.1389 - 0.0789i-0.0918 - 0.0132i0.4+j0.050.2905 - 0.1333i0.2028 + 0.0746i0.2+j0.2gg125340.2956 - 0.1229i-0.06+0.1i0.2931 + 0.0031i0.0677 + 0.1531i0.2185 + 0.0764i0.2991 - 0.0163i0.2109 + 0.0812i-0.0918 - 0.0132i-0.092 - 0.013i-0.0918 - 0.0132i0.0771 - 0.0758i-0.1448 - 0.0773i0.1960 + 0.0786i五，小结在计算机潮流计算中，使用牛顿-拉夫逊法进行潮流计算，非常快捷。在程序编写过程中，对公式的由来应该非常清楚 ，不然在编写过程中易出错，比如for循环的次数，一些细节问题都应该反复琢磨。在程序调试过程中，应该分步调试，不要一开始就运行完整的程序。对程序中出现的问题应该从运行结果中分析问题的所在，有什么问题解决什么问题，抓住问题的要害。在潮流计算中运用matlab比运用c语言要方便好多，matlab里有大量的数学函数直接调用。用c语言写的话每一个函数都要自己编写，写起来很痛苦,其中求逆矩阵就用了一百多行。最后想说的是，想学没有什么不可以。六，程序代码matlab程序：%writer:chenzhou%time:2009年11月9日-2009年11月15日%编译环境：matlab r2007bclc%说明：为了使节点按照先pq，再pv节点，最后平衡节点的次序编号，以便与公式对照，节点1与节点5对调。%节点阻抗矩阵z=0,0.04+0.12i,0,0.08+0.24i,0;0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.02+0.06i;0,0.06+0.18i,0,0.01+0.03i,0.08+0.24i;0.08+0.24i,0.06+0.18i,0.01+0.03i,0,0;0,0.02+0.06i,0.08+0.24i,0,0;%求互导纳for i=1:5 for j=1:5 if z(i,j)=0 y(i,j)=0; else y(i,j)=1/z(i,j); end endend%求导纳for i=1:5 for j=1:5 if i=j y(i,j)=-y(i,j); else for k=1:5 y(i,j)=sum(y(i,:); end end endend%导纳矩阵g=real(y);b=imag(y);%赋初值u=1,1,1,1,1.06;angle=0,0,0,0,0;%节点功率s=0.6-0.1i,0.2+0.2i,-0.45-0.15i,-0.4-0.05i,0;p=real(s);q=imag(s);error_power=ones(8,1);count=0;%迭代次数while max(error_power)1e-5 for i=1:4 h(i,i)=0;n(i,i)=0;m(i,i)=0;l(i,i)=0;error_p(i)=0;error_q(i)=0; end for i=1:4 for j=1:5 error_p(i)=error_p(i)-u(i)*u(j)*(g(i,j)*cos(angle(i)-angle(j)+b(i,j)*sin(angle(i)-angle(j); error_q(i)=error_q(i)-u(i)*u(j)*(g(i,j)*sin(angle(i)-angle(j)-b(i,j)*cos(angle(i)-angle(j); end error_p(i)=error_p(i)+p(i); error_q(i)=error_q(i)+q(i); enderror_power=error_p,error_q;%求雅可比矩阵%当i=j时h,n,m,l如下 for i=1:4 for j=1:4 if i=j h(i,j)=-u(i)*u(j)*(g(i,j)*sin(angle(i)-angle(j)-b(i,j)*cos(angle(i)-angle(j); n(i,j)=-u(i)*u(j)*(g(i,j)*cos(angle(i)-angle(j)+b(i,j)*sin(angle(i)-angle(j); l(i,j)=-u(i)*u(j)*(g(i,j)*sin(angle(i)-angle(j)-b(i,j)*cos(angle(i)-angle(j); m(i,j)=u(i)*u(j)*(g(i,j)*cos(angle(i)-angle(j)+b(i,j)*sin(angle(i)-angle(j); end endend%当i=j时h，n，m，l如下for i=1:4 for j=1:5 if i=j h(i,i)=h(i,i)+u(i)*u(j)*(g(i,j)*sin(angle(i)-angle(j)-b(i,j)*cos(angle(i)-angle(j); n(i,i)=n(i,i)-u(i)*u(j)*(g(i,j)*cos(angle(i)-angle(j)+b(i,j)*sin(angle(i)-angle(j); m(i,i)=m(i,i)-u(i)*u(j)*(g(i,j)*cos(angle(i)-angle(j)+b(i,j)*sin(angle(i)-angle(j); l(i,i)=l(i,i)-u(i)*u(j)*(g(i,j)*sin(angle(i)-angle(j)-b(i,j)*cos(angle(i)-angle(j); end end n(i,i)=n(i,i)-2*(u(i)2*g(i,i); l(i,i)=l(i,i)+2*(u(i)2*b(i,i);endj=h,n;m,l;error_angle_u=-(inv(j)*error_power);error_angle=zeros(1,4);error_u=zeros(1,4);for i=1:4 error_angle(i)=error_angle_u(i);error_u(i)=error_angle_u(i+4)*u(i);endfor i=1:4 angle(i)=angle(i)+error_angle(i); u(i)=u(i)+error_u(i);endcount=count+1;end%count%平衡节点的注入功率 for j=1:5 p(5)=u(5)*u(j)*(g(5,j)*cos(angle(5)-angle(j)+b(5,j)*sin(angle(5)-angle(j)+p(5); q(5)=u(5)*u(j)*(g(5,j)*sin(angle(5)-angle(j)-b(5,j)*cos(angle(5)-angle(j)+q(5); end s(5)=p(5)+q(5)*sqrt(-1); %节点电压用复数表示 for i=1:5 complex_u(i)=u(i)*exp(angle(i)*sqrt(-1); end %元件端功率 for i=1:5 for j=1:5 line_s(i,j)=conj(complex_u(i)*(conj(complex_u(i)-conj(complex_u(j)*conj(y(i,j); end end %元件端电流 for i=1:5 for j=1:5 line_i(i,j)=conj(line_s(i,j)/conj(complex_u(i); end end y,count,u,angle,error_angle_u,s,line_s,line_i c语言程序：此程序一共三个文件，分别为inversematrix.h，functions.h，chaoliujisuan.c 主函数：/writer:chenzhou/time:2009年11月9日-2009年11月15日/编译环境：visual stdio 2008/说明：为了使节点按照先pq，再pv节点，最后平衡节点的次序编号，以便与公式对照，节点与节点对调。#include#include#includeinversematrix.h#includefunctions.hvoid main()unsigned int i,j,count=0;z_to_y();printf(节点导纳矩阵：n);for (i=0;i5;i+ )for (j=0;j0.00001)|(max_array(error_q)0.00001)one_times();count+;printf(迭代次数：count=%d,count+1);printf(n);printf(节点相角分别为:=);for(i=0;i5;i+)printf(%f ,anglei);printf(n);printf(节点电压分别为:u=);for(i=0;i5;i+)printf(%f ,ui);calculate_s();printf(n节点功率为:s=n);for (i=0;i5;i+ )printf(%f+j%f ,s_pi,s_qi);printf(n); 求逆矩阵的头文件（inversematrix.h文件）int max=8;/=雅可比矩阵=double j88=0;int switchnum28;void switchrow(int r,int k)/*将矩阵中的r,k行交换*/int i;double temp;for (i=0;imax;i+ ) temp=jki; jki=jri; jri=temp;void switchcol(int r,int k)/*将矩阵中的r,k行交换*/int i;double temp;for (i=0;imax;i+ ) temp=jik; jik=jir; jir=temp;int inverse_matrix() int i,j,r,k;double temp;for (k=0;kmax ;k+ ) temp=0; r=k; for (j=k;j=max; j+)/*找主元*/ if (tempfabs(jrj) temp=fabs(jrj); r=j; if (0=temp)/*非奇异矩阵,退出*/ return(0); if(r!=k) switchrow(r,k); switchnum0k=r; switchnum1k=k; else switchnum0k=0; switchnum1k=0; jkk=1/jkk; for (j=0;jmax;j+ ) if (j!=k) jkj=-jkj*jkk; for (i=0;imax ;i+ ) for (j=0;jmax;j+ ) if (i!=k & j!=k) jij=jij+jik*jkj; for (i=0;i=0 ;k-) if(switchnum0k!=0) switchcol(switchnum0k,switchnum1k); 求解过程中的一些功能函数（functions.h文件）/=计算导纳矩阵=/阻抗的实部double z_re55=0,0.04,0,0.08,0,0.04,0,0.06,0.06,0.02,0,0.06,0,0.01,0.08,0.08,0.06,0.01,0,0,0,0.02,0.08,0,0;/阻抗的虚部double z_im55=0,0.12,0,0.24,0,0.12,0,0.18,0.18,0.06,0,0.18,0,0.03,0.24,0.24,0.18,0.03,0,0,0,0.06,0.24,0,0;/导纳矩阵double y_re55=0;double y_im55=0;/=节点电压=double u5=1,1,1,1,1.06;/=节点电压相角=double angle5=0;/=节点电压相角补偿=double error_angle5=0;/=功率不平衡量=double error_p5=0;double error_q5=0;/=节点注入的有功和无功功率=double s_p5=0;double s_q5=0;void z_to_y()/z=m+nj y_re=m/(m*m+n*n) y_im=-n/(m*m+n*n)junsigned char i=0,j=0,k;double y_temp_re55=0;double y_temp_im55=0;for(i=0;i5;i+)for(j=0;j5;j+)if(z_reij=0)&(z_imij=0)y_temp_reij=0;y_temp_imij=0;elsey_temp_reij=z_reij/(z_reij*z_reij+z_imij*z_imij);y_temp_imij=-z_imij/(z_reij*z_reij+z_imij*z_imij);/=导纳矩阵=for(i=0;i5;i+)for(j=0;j5;j+)if(i=j)for(k=0;k5;k+)y_reij=y_reij+y_temp_reik;y_imij=y_imij+y_temp_imik;else y_reij=-y_temp_reij;y_imij=-y_temp_imij;/=计算功率不平衡量=void error_power()unsigned char i,j;double sum_p5=0;double sum_q5=0;double start_p5=0.6,0.2,-0.45,-0.4,0;double start_q5=-0.1,0.2,-0.15,-0.05,0;for(i=0;i5;i+)for(j=0;j5;j+)sum_pi=sum_pi+ui*uj*(y_reij*cos(anglei-anglej)+y_imij*sin(anglei-anglej);sum_qi=sum_qi+ui*uj*(y_reij*sin(anglei-anglej)-y_imij*cos(anglei-anglej);for(i=0;i5;i+)error_pi=start_pi-sum_pi;error_qi=start_qi-sum_qi;/=形成修正方程式=void calculate_j()unsigned char i=0,j=0,k=0;double sum=0;for(i=0;i4;i+)/=h参数=for(j=0;j4;j+)if(i=j)for(k=0;k5;k+)if(i!=k)sum=sum+ui*uk*(y_reik*sin(anglei-anglek)-y_imik*cos(anglei-anglek);else sum=sum+0;jij=sum;sum=0;else jij=-ui*uj*(y_reij*sin(anglei-anglej)-y_imij*cos(anglei-anglej);for(i=0;i4;i+)/=n参数=for(j=0;j4;j+)if(i=j)for(k=0;k5;k+)if(i!=k)sum=sum-ui*uk*(y_reik*cos(anglei-anglek)+y_imik*sin(anglei-anglek);else sum=sum+0;jij+4=sum-2*ui*ui*y_reii;sum=0;else jij+4=-ui*uj*(y_reij*cos(anglei-anglej)+y_imij*sin(anglei-anglej);for(i=0;i4;i+)/=m参数=for(j=0;j4;j+)if(i=j)for(k=0;k5;k+)if(i!=k)sum=sum-ui*uk*(y_reik*cos(anglei-anglek)