珍藏多年的matlab潮流计算程序源代码集合，包含多个潮流计算程序

// power flow2.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include"math.h" #include"stdio.h" #include"stdlib.h" #include "Conio.h" #define N 4 /*定义节点数*/ #define M 2 /*定义PQ节点数*/ #define V N-M-1 /*定义PV节点数*/ #define N 4 /*定义节点数*/ #define M 2 /*定义PQ节点数*/ #define V N-M-1 /*定义PV节点数*/ #define n 6 #define precision 1e-16 int gauss(double a[][n+2],double x[n+1]) {int i,j,k,r;double c; for(k=1;k<=n-1;k++) {r=k; for(i=k;i<=n;i++) if(fabs(a[i][k])>fabs(a[r][k]))r=i; if(fabs(a[r][k])<precision) {printf("\n det A=0.failed!");exit(0);} if(r!=k) {for(j=k;j<=n+1;j++) {c=a[k][j];a[k][j]=a[r][j];a[r][j]=c;} } for(i=k+1;i<=n;i++) {c=a[i][k]/a[k][k]; for(j=k+1;j<=n+1;j++) a[i][j]=a[i][j]-c*a[k][j]; } } if(fabs(a[n][n])<precision) {printf("\n det A=0.failed!");exit(0);} for(k=n;k>=1;k--) {x[k]=a[k][n+1]; for(j=k+1;j<=n;j++) x[k]=x[k]-a[k][j]*x[j]; x[k]=x[k]/a[k][k]; } return(1); } int _tmain(int argc, _TCHAR* argv[]) { int k=0,i,j; /*k为迭代次数,n为节点个数*/ double temp1,temp2; double g[N][N]={{1.0421,-0.5882,0.0,-0.4539},{-0.5822,1.069,0.0,-0.4808}, {0.0,0.0,0.0,0.0},{-0.4539,-0.4808,0.0,0.9346}},b[N][N]={{-8.2429,2.3529,3.6666,1.8911}, {2.3529,-4.7274,0.0,2.4038},{3.6666,0.0,-3.3333,0.0},{1.8911,2.4038,0.0,-4.2616}}; double e[N]={1.0,1.0,1.0,1.05},f[N]={0.0,0.0,0.0,0.0}; double p[N]={-0.3,-0.55,0.5,0.0},q[N]={-0.18,-0.13,0.0,0.0}; double u[N]={0.0,0.0,1.1,1.05}; double dp[N]={0.0},dq[N]={0.0}; double du2[N]={0.0}; double jac[2*(N-1)][2*(N-1)]; double h[N-1][N-1],nn[N-1][N-1],jj[M][N-1],l[M][N-1],r[N-1][N-1],s[N-1][N-1]; double dpq[2*(N-1)],dfe[2*(N-1)]; double fe[2*(N-1)]={0.0,1.0,0.0,1.0,0.0,1.0}; int re,det; double aa[n][n+1],a[n+1][n+2],x[n+1]; do { for(i=0;i<N-1;i++) /*计算节点功率,电压不平衡量*/ { temp1=0; temp2=0; for(j=0;j<N;j++) { temp1+=(g[i][j]*e[j]-b[i][j]*f[j]); temp2+=(g[i][j]*f[j]+b[i][j]*e[j]); } dp[i]=p[i]-e[i]*temp1-f[i]*temp2; } for(i=0;i<N-2;i++) { temp1=0; temp2=0; for(j=0;j<N;j++) { temp1+=(g[i][j]*e[j]-b[i][j]*f[j]); temp2+=(g[i][j]*f[j]+b[i][j]*e[j]); } dq[i]=q[i]-f[i]*temp1+e[i]*temp2; } for(i=M;i<N-1;i++) { du2[i]=u[i]*u[i]-(e[i]*e[i]+f[i]*f[i]); } for(i=0;i<N-1;i++) //显示dp printf("dp(%d)=%f\n",i+1,dp[i]); for(i=0;i<M;i++) //显示dq printf("dq(%d)=%f\n",i+1,dq[i]); for(i=M;i<(M+V);i++) //显示du2 printf("du2(%d)=%f\n",i+1,du2[i]); for(i=0;i<N-1;i++) /*计算雅克比矩阵各元素*/ { temp1=0; temp2=0; for(j=0;j<N;j++) {temp1+=(g[i][j]*f[j]+b[i][j]*e[j]); temp2+=(g[i][j]*e[j]-b[i][j]*f[j]); } for(j=0;j<N-1;j++) /*计算H,N元素*/ {if(j!=i) { h[i][j]=(-1)*b[i][j]*e[i]+g[i][j]*f[i]; nn[i][j]=g[i][j]*e[i]+b[i][j]*f[i]; } if(j==i) { h[i][i]=(-1)*b[i][i]*e[i]+g[i][i]*f[i]+temp1; nn[i][i]=g[i][i]*e[i]+b[i][i]*f[i]+temp2; } } } for(i=0;i<M;i++) /*计算J,L元素*/ { temp1=0; temp2=0; for(j=0;j<N;j++) {temp1+=(g[i][j]*f[j]+b[i][j]*e[j]); temp2+=(g[i][j]*e[j]-b[i][j]*f[j]); } for(j=0;j<N-1;j++) {if(j!=i) { jj[i][j]=(-1)*nn[i][j]; l[i][j]=h[i][j]; } if(j==i) { jj[i][i]=(-1)*g[i][i]*e[i]-b[i][i]*f[i]+temp2; l[i][i]=(-1)*b[i][i]*e[i]+g[i][i]*f[i]-temp1; } } } for(i=M;i<N-1;i++) //计算R,S元素 {for(j=0;j<N-1;j++) {if(j!=i) { r[i][j]=0; s[i][j]=0; } if(j==i) { r[i][i]=2*f[i]; s[i][i]=2*e[i]; } } } for(i=0;i<N-1;i++) /*将H,N元素放入雅克比矩阵*/ { for(j=0;j<N-1;j++) {jac[2*i][2*j]=h[i][j]; jac[2*i][2*j+1]=nn[i][j]; } } for(i=0;i<M;i++) /*将J,L元素放入雅克比矩阵*/ { for(j=0;j<N-1;j++) { jac[2*i+1][2*j]=jj[i][j]; jac[2*i+1][2*j+1]=l[i][j]; } } for(i=M;i<N-1;i++) /*将R,S元素放入雅克比矩阵*/ { for(j=0;j<N-1;j++) { jac[2*M+1][2*j]=r[i][j]; jac[2*M+1][2*j+1]=s[i][j]; } } for(i=0;i<N-1;i++) /*得到功率,电压不平衡向量*/ { dpq[2*i]=dp[i]; } for(i=0;i<M;i++) { dpq[2*i+1]=dq[i]; } for(i=0;i<V;i++) { dpq[2*M+1+i]=du2[M]; } for(i=0;i<n;i++) {for(j=0;j<n;j++) aa[i][j]=jac[i][j]; aa[i][j]=dpq[i]; } for(i=1;i<=n;i++)for(j=1;j<=n+1;j++) a[i][j]=aa[i-1][j-1]; det=gauss(a,x);/*用高斯消元法解修正方程*/ for(i=0;i<2*(N-1);i++) /*得到df与de的解*/ dfe[i]=x[i+1]; k=k+1; /*迭代次数加1*/ if(k>=10) /*迭代次数过多,出错*/ { printf("error"); return(0); } for(i=0;i<2*(N-1);i++) /*得到下一轮迭代f,e的初值*/ { fe[i]+=dfe[i]; } for(i=0;i<(N-1);i++) { f[i]=fe[2*i]; e[i]=fe[2*i+1]; } for(i=0;i<2*(N-1);i++) dfe[i]=fabs(dfe[i]); for(i=0;i<2*(N-1);i++) /*判断收敛条件*/ { if(dfe[i]<=0.00001) re=0; else {re=1; break;} } } while(re); for(i=0;i<N;i++) /*显示节点电压*/ {printf("e(%d)=%f\n",i+1,e[i]); printf("f(%d)=%f\n",i+1,f[i]); } return 0; }