c语言编写的PQ分解法解潮流程序

1、这里要提一下，本人使用的数据出自鞠平主编电力工程，PQ分解法可能由于求逆程序的精确度问题而比matlab迭代次数稍多，可将float改为dobule float一试，若有不当之处，还望包涵（另外，本人使用c语言编译器为tc3.0）/PQ分解法解潮流程序#include#include#define N 4 /节点数#define n_PQ 2 /PQ节点数#define n_PV 1 /PV节点数#define n_br 5 /串联支路数#define PI 3.1415void main()void disp_matrix(float *disp_p,int disp_m,int disp_

2、n); /矩阵显示函数float Us2*N=1.0,0,1.0,0,1.05,0,1.05,0; /电压初值float PsN=0,-0.5,0.2; /有功初值float QsN=0,-0.3; /无功初值 float GNN,BNN,B1N-1N-1,B2n_PQn_PQ,invB1N-1N-1,invB2n_PQn_PQ; struct /阻抗参数 int nl; /左节点 int nr; /右节点 float R; /串联电阻值 float X; /串联电抗值 float Bl; /左节点并联电纳 float Br; /右节点并联电纳 ydatan_br= 1,2,0,0.1880,

3、-0.6815,0.6040, 1,3,0.1302,0.2479,0.0129,0.0129, 1,4,0.1736,0.3306,0.0172,0.0172, 3,4,0.2603,0.4959,0.0259,0.0259, 2,2,0,0.05,0,0 ;float Z2; /Z2=R2+X2 各串联阻抗值的平方float U_ampN,U_angN; /e,f存储电压的x轴分量和y轴分量,dfe存储电压修正值float dU_ampN,dU_angN; /电压幅值和相角修正值float mid1N,mid2N,dPN-1,dQn_PQ; /mid1、mid2存储计算雅克比行列式对角线元

4、素的中间值,dS存储PQU的不平衡量float dPQU=1.0; /PQU不平衡量最大值int kk=0; /迭代次数int i,j,k; float t; /形成导纳矩阵-for(i=0;iN;i+) for(j=0;jN;j+) Gij=0; Bij=0; for(i=0;in_br;i+) if(ydatai.nl!=ydatai.nr) Z2=(ydatai.R)*(ydatai.R)+(ydatai.X)*(ydatai.X); /串联阻抗等效导纳值 /非对角元素 Gydatai.nl-1ydatai.nr-1=(-ydatai.R)/Z2; Bydatai.nl-1ydatai.

5、nr-1=ydatai.X/Z2; Gydatai.nr-1ydatai.nl-1=(-ydatai.R)/Z2; Bydatai.nr-1ydatai.nl-1=ydatai.X/Z2; /对角元素 Gydatai.nl-1ydatai.nl-1+=ydatai.R/Z2; Gydatai.nr-1ydatai.nr-1+=ydatai.R/Z2; Bydatai.nl-1ydatai.nl-1+=(-ydatai.X/Z2); Bydatai.nr-1ydatai.nr-1+=(-ydatai.X/Z2); /并联导纳等效导纳值 Bydatai.nl-1ydatai.nl-1+=ydata

6、i.Bl; Bydatai.nr-1ydatai.nr-1+=ydatai.Br; else Gydatai.nl-1ydatai.nr-1+=ydatai.R; Bydatai.nl-1ydatai.nr-1+=ydatai.X; printf(G=n);disp_matrix(*G,N,N);printf(B=n);disp_matrix(*B,N,N);/分离U,deltafor(i=0;i0.00001) for(i=0;iN-1;i+) mid1i=0; mid2i=0; for(j=0;jN;j+) mid1i=mid1i+Gij*U_ampj*cos(U_angj)-Bij*U_

7、ampj*sin(U_angj); mid2i=mid2i+Gij*U_ampj*sin(U_angj)+Bij*U_ampj*cos(U_angj); dPi=Psi-(U_ampi*cos(U_angi)*mid1i+U_ampi*sin(U_angi)*mid2i); if(in_PQ) dQi=Qsi-(U_ampi*sin(U_angi)*mid1i-U_ampi*cos(U_angi)*mid2i); dPQU=0; for(i=0;iN-1;i+) if(dPi0&dPQU0&dPQUdPi) dPQU=dPi; for(i=0;in_PQ;i+) if(dQi0&dPQU0&d

8、PQU0.00001) kk+; /形成B,B矩阵 for(i=0;iN-1;i+) for(j=0;jN-1;j+) B1ij=Bij; if(in_PQ&jn_PQ) B2ij=Bij; /求P,delta修正值 for(i=0;iN-1;i+) dPi=dPi/U_ampi; if(in_PQ) dQi=dQi/U_ampi; /B求逆- for(i=0;iN-1;i+) for(j=0;jN-1;j+) if(i!=j) invB1ij=0; else invB1ij=1; for(i=0;iN-1;i+) for(j=0;jN-1;j+) if(i!=j) t=B1ji/B1ii;

9、for(k=0;kN-1;k+) B1jk-=B1ik*t; invB1jk-=invB1ik*t; for(i=0;iN-1;i+) if(B1ii!=1) t=B1ii; for(j=0;jN-1;j+) invB1ij=invB1ij/t; for(i=0;iN-1;i+) dU_angi=0; for(i=0;iN-1;i+) for(j=0;jN-1;j+) dU_angi+=invB1ij*dPj; dU_angi=dU_angi/U_ampi; /Q,U_amp修正值 /B求逆 for(i=0;in_PQ;i+) for(j=0;jn_PQ;j+) if(i!=j) invB2i

10、j=0; else invB2ij=1; for(i=0;in_PQ;i+) for(j=0;jn_PQ;j+) if(i!=j) t=B2ji/B2ii; for(k=0;kn_PQ;k+) B2jk-=B2ik*t; invB2jk-=invB2ik*t; for(i=0;in_PQ;i+) if(B2ii!=1) t=B2ii; for(j=0;jn_PQ;j+) invB2ij=invB2ij/t; for(i=0;in_PQ;i+) dQi=dQi/U_ampi; for(i=0;in_PQ;i+) dU_ampi=0; for(i=0;in_PQ;i+) for(j=0;jn_PQ

11、;j+) dU_ampi+=invB2ij*dQj; for(i=0;iN;i+) U_ampi-=dU_ampi; U_angi-=dU_angi; else break;/循环结束-/求平衡节点功率-mid1N-1=0;mid2N-1=0;for(j=0;jN;j+) mid1N-1=mid1N-1+GN-1j*U_ampj*cos(U_angj)-BN-1j*U_ampj*sin(U_angj); mid2N-1=mid2N-1+GN-1j*U_ampj*sin(U_angj)+BN-1j*U_ampj*cos(U_angj);PsN-1=U_ampN-1*cos(U_angN-1)*m

12、id1N-1+U_ampN-1*sin(U_angN-1)*mid2N-1;QsN-1=U_ampN-1*sin(U_angN-1)*mid1N-1-U_ampN-1*cos(U_angN-1)*mid2N-1;for(i=n_PQ;iN-1;i+) Qsi=U_ampi*sin(U_angi)*mid1i-U_ampi*cos(U_angi)*mid2i;/-/ 显示输出结果printf(kk=%dn,kk);printf(P=);for(i=0;iN;i+) printf(%9.4f,Psi);printf(nQ=);for(i=0;iN;i+) printf(%9.4f,Qsi);printf(nU_amp=);for(i=0;iN;i+) printf(%9.4f,U_ampi);printf(nU_ang(deg)=);for(i=0;iN;i+) pri