第一篇：電力系統(tǒng)潮流計算的C語言實現(xiàn)
////////////////////////////////////////////////////////////////////// //
PQ分解法潮流
// //文件輸入格式：節(jié)點總數(shù)n(包括聯(lián)絡節(jié)點),支路數(shù)zls
//節(jié)點數(shù)(發(fā)電機和負荷)nb,接地電抗數(shù)mdk,迭代精度eps
//考慮負荷靜特性標志kk2(0考慮),平衡節(jié)點號，優(yōu)化標志(0不優(yōu)化)
//最大迭代次數(shù)it1,支路左右節(jié)點號izl[],jzl[],支路電阻zr[],電抗zx[] // //支路容納zyk[],節(jié)點號nob[]及標志nobt[](0-PQ,-1-PV)
//發(fā)電機和負荷有功、無功pg[],qg[],pl[],ql[]
//電壓v0[](pv節(jié)點輸入實際值，PQ節(jié)點任輸入一值)
//電抗節(jié)點號idk[],電抗值dkk[]
//////////////////////////////////////////////////////////////////////
#include “math.h” #include “stdio.h”
#define NS 2000 //最大節(jié)點數(shù) #define NS2 NS * 2 #define NS4 #define ZS #define ZS2 1000 3000 ZS * 2
//NS4、NS必須大于2*zls。//最大支路數(shù)
// // // // // // // #define DKS 200
//最大電抗器數(shù) #define N2 ZS * 4 #define N3 ZS * 8 + NS * 4
FILE *fp1, *fp2;char inname[12], outname[12];// fp1輸入數(shù)據文件指針 fp2輸出文件指針
// inname[]輸入數(shù)據文件名 outname[]輸出數(shù)據文件名 int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;// 節(jié)點總數(shù)n(包括聯(lián)絡節(jié)點)
支路數(shù)(回路數(shù))zls 節(jié)點數(shù)nb(發(fā)電機和負荷)
// 接地電抗數(shù)mdk 精度eps
平衡節(jié)點號mpj
// 節(jié)點優(yōu)化(標志)bnsopton(=0節(jié)點不優(yōu)化，!=0節(jié)點優(yōu)化)
// 最大迭代次數(shù)it1 最低電壓或最大功率誤差節(jié)點號dsd // 負荷靜特性標志(=0考慮負荷靜特性)// 支路數(shù)(雙回線算一條支路)int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];// izl[],jzl[],idk[]:分別存放左、右節(jié)點號和電抗器節(jié)點號。// yds[]存放各行非零非對角元素的個數(shù)。
// ydz[i]是第 i 行第一個非零非對角元素的首地址，// 即在所有非零非對角元素中的次序號 // iy[]存放列足碼。
int nnew[NS4], old[NS], nob[NS], nobt[NS];// nnew[],old[]存放的是新、舊節(jié)點號。// nnew[i]中為i對應的新號 // nob[]存放的是節(jié)點號。nobt[]存放的是節(jié)點類型, 0: pq節(jié)點,-1: pv節(jié)點。double eps, dsm, vmin, dph, dqh, af[3];// eps迭代收斂精度，dsm最大功率誤差
// vmin:系統(tǒng)最低電壓值。dph,dqh:系統(tǒng)有、無功損耗。
// af[0]和af[1]分別是負荷有功功率、無功功率靜態(tài)特性系數(shù)。double v00;
// v00: 系統(tǒng)平均電壓 ci,cj分別作為節(jié)點i,j的電壓相角的臨時存儲單元。double zr[ZS], zx[ZS], zyk[ZS], dkk[DKS], gii[NS], bii[NS], yg[ZS2], yb[ZS2];double pg[NS], qg[NS], pl[NS], ql[NS], v0[NS], v[NS], va[NS];// 支路電阻zr[] 支路電抗zx[] 輸電線路充電容納zyk[](y0/2)// 接地電抗dkk[] 對角元實部gii[] 對角元虛部 // 非對角元實部yg[] 非對角元虛部yb[] // pg[],qg[],pl[],ql[]:發(fā)電機，負荷功率實、虛部 // v[]是電壓幅值,va[]是電壓相角。double w[NS2], kg[3], b[NS2];int newsort[NS4];// newsort[i]存放i對應的老號
void initial();void pqflow();void out();void dataio();void bnsopt();void zlsort(int* nnew);void printo();void printy();void y2();void ya0();void yzb();void jdgl(int kq0);void bbhl(int kq0);void calc();int iabs(int a);void branch_output();void newval(double* aa);void printc();void iswap();void swap();void printf2(double* aa, double* bb, int n);void calc(int* iu, double* u, double* di, int* nfd, double* b);void printi(int* aa, int n);void printf1(double* aa, int n);int find(int k, int a[], int* z);void yzb(int t, int* iu, double* u, double* di, int* nfd);int isgn(int a, int b);void yy1();void y3();void newtoold();
int main(void){ initial();//初始化
pqflow();//pq潮流計算
out();//輸出節(jié)點和支路數(shù)據
return 1;}
int isgn(int a, int b){ //**** 本函數(shù)功能返回值為a的絕對值b的符號 //參數(shù)1提供值，參數(shù)2提供符號//
if(b < 0)
if(a > 0)
a =-a;return a;}
int find(int k, int a[], int* z){
****// //**** 本函數(shù)查找a[]中是否有fabs(k)有則返回0，無則返回1 ****// //參數(shù)1為待查找量，參數(shù)2待搜索數(shù)組，參數(shù)3返回k在a[]中的次序號//
int i;for(i = 1;i <= n;i ++)
if(iabs(k)== a[i])
{
*z = i;
return 1;
} return 0;}
void oldtonew(){ //**** 本函數(shù)將輸入數(shù)據中的節(jié)點號變成從1開始的連續(xù)節(jié)點號 ****//
int i, j, k, ii1, ii2, zls2, k1, k2, k3, k4, ip;zls2 = zls + zls;
for(i = 1;i <= zls2;i ++)
newsort[i] = 0;ii1 = 0;
for(i = 1;i <= zls;i ++){
k = izl[i];
if(!find(k, newsort, &ii2))
{
ii1 ++;
newsort[ii1] = iabs(k);
}
k = jzl[i];
if(!find(k, newsort, &ii2))
{
ii1 ++;
newsort[ii1] = iabs(k);
} } for(i = 1;i <= ii1-1;i ++){
for(j = i+1;j <= ii1;j ++)
{
if(newsort[i] > newsort[j])
{
k = newsort[i];
newsort[i] = newsort[j];
newsort[j] = k;
}
} } for(i = 1;i <= zls;i ++){
k = izl[i];
if(find(k, newsort, &ii2))
{
izl[i] = isgn(ii2, k);
} } else
printf(“error!”);k = jzl[i];if(find(k, newsort, &ii2)){
jzl[i] = isgn(ii2, k);} else
printf(“error!”);printf(“izl[%d] = %d, jzl[%d] = %dn”, i, izl[i], i, jzl[i]);} for(i = 1;i <= nb;i ++){ for(j = 1;j <= n;j ++)
if(nob[i] == newsort[j])
{
nob[i] = j;
break;
} printf(“nob[%d] = %dn”, i, nob[i]);} for(j = 1;j <= n;j ++){ if(mpj == newsort[j]){
mpj = j;
break;} } //電抗器節(jié)點號轉變 for(j = 1;j <= mdk;j ++){ for(i = 1;i <= n;i ++){
if(idk[j] == newsort[i])
{
idk[j] = i;
break;
} } } void initial(){
//**** 本函數(shù)進行初始化工作 ****//
int i, k1;
dataio();//輸入原始數(shù)據
oldtonew();//轉化為新號
if(bnsopton == 0)
//節(jié)點不優(yōu)化，新節(jié)點號即為老節(jié)點號。
for(i = 1;i <= n;i ++)
{
old[i] = i;
nnew[i] = i;
}
else
bnsopt();//節(jié)點優(yōu)化
mpj = nnew[mpj];//mpj:平衡節(jié)點 zlsort(nnew);// sort the r,x and b for(i = 1;i <= mdk;i ++){
k1 = idk[i];
idk[i] = nnew[k1];} for(i = 1;i <= n;i ++)
{
v[i] = v00;
va[i] = 0.0;} // 所有節(jié)點的電壓幅值初值都為1.000(v00)，電壓相角初值都為0。
// exchange the node before and after sort for(i = 1;i <= n;i ++)
yds[i] = 0;
// the immediate for(i = 1;i <= nb;i ++)
{
k1 = nnew[nob[i]];
yds[k1] = nobt[i];} for(i = 1;i <= n;i ++)
nobt[i] = yds[i];newval(pg);newval(qg);newval(pl);newval(ql);newval(v0);
for(i = 1;i <= n;i ++)// nobt[] is type of node if(nobt[i] ==-1)
v[i] = v0[i];
//
nob[] is serials numbe //nobt[] =-1: pv節(jié)點，v0[]存放的是最后一個節(jié)點數(shù)據，//對于pv節(jié)點，即為該點應維持的電壓值。
//nobt[] = 0: pq節(jié)點，v0[]存放的是最后一個節(jié)點數(shù)據，//對于pq節(jié)點，即為系統(tǒng)平均電壓值。printo();//輸出af[]、v00和節(jié)點排序后的支路、節(jié)點和
//接地電抗數(shù)據(僅僅查看中間結果)ya0();//獲得yds[]、ydz[]、列足碼iy[]。(P407)}
void printo(){ //**** 輸出af[]、v00和節(jié)點排序后的支路、節(jié)點和接地電抗數(shù)據 ****//
int i;fprintf(fp2, “n
******AF AND V0 ******n”);fprintf(fp2, “n %7.3f%7.3f%7.3fn”, af[0], af[1], v00);printc('-', 78);fprintf(fp2, “nn
*******ZLB*******n”);for(i = 1;i <= zls;i ++)
{
fprintf(fp2, “n”);
fprintf(fp2, “%8d%8d%8d%8d”, izl[i], jzl[i], old[abs(izl[i])], old[abs(jzl[i])]);
fprintf(fp2, “%9.4f%9.4f%9.4f”, zr[i], zx[i], zyk[i]);
} printc('-', 78);fprintf(fp2, “nn*******BUS*******n”);for(i = 1;i <= nb;i ++)
{
fprintf(fp2, “n”);
fprintf(fp2, “%8d%8d%8d”, nob[i],old[nob[i]], nobt[i]);
fprintf(fp2, “%9.4f%9.4f%9.4f%9.4f%9.4f”, pg[i], qg[i], pl[i], ql[i], v0[i]);
} printc('-', 78);fprintf(fp2,“nn******DKK******n”);for(i = 1;i <= mdk;i ++)
{
fprintf(fp2, “n”);
fprintf(fp2, “%8d%8d%7.4f”, idk[i], old[idk[i]], dkk[i]);
} } void dataio()
{
//**** 系統(tǒng)數(shù)據初始化 ****//
int i;af[0] = 0.6;
af[1] = 2.0;//af[0]和af[1]分別是負荷有功功率、無功功率靜態(tài)特性系數(shù)。
v00 = 1.000;//系統(tǒng)平均電壓
printf(“nplease input the name of data filen”);scanf(“%s”, inname);fp1 = fopen(inname, “r”);printf(“nplease output the name of data filen”);scanf(“%s”, outname);fp2 = fopen(outname, “w”);fscanf(fp1, “%d %d %d %d”, &n, &zls, &nb, &mdk);
// the number of node ,branches, node fscanf(fp1, “%lf %d %d %d %d”, &eps, &kk2, &mpj，&bnsopton, &it1);
//precision, swing node,sort the node,iteration numbers for(i = 1;i <= zls;i ++)
{
fscanf(fp1, “%d %d”, &izl[i], &jzl[i]);
fscanf(fp1, “%lf %lf %lf ”, &zr[i], &zx[i], &zyk[i]);
} for(i = 1;i <= nb;i ++){
fscanf(fp1, “%d %d”, &nob[i], &nobt[i]);
fscanf(fp1, “%lf %lf %lf %lf %lf”, &pg[i], &qg[i], &pl[i]，&ql[i], &v0[i]);} for(i = 1;i <= mdk;i ++)
{
fscanf(fp1, “%d %lf”, &idk[i], &dkk[i]);
} fclose(fp1);}
void pqflow()
{
//**** PQ分解法計算潮流,程序框圖見P164圖3-16(從第 7 步起)****//
int kq0, iu1[N2], nfd1[NS], iu2[N2], nfd2[NS];int i, t;
double u1[N2], u2[N2], di1[NS], di2[NS];yy1();
yzb(0, iu1, u1, di1, nfd1);//form the B matrix of P-0 iteration y2();yzb(1, iu2, u2, di2, nfd2);//form the B matrix of Q-V iteration t = 0;
kq0 = 0;kg[0] = kg[1] = 1;do { jdgl(kq0);// calculating the power bbhl(kq0);// find out the maxi if(kq0 == 0)
printf(“P: %dt%dt%fn”, t, dsd, dsm);else
printf(“Q: %dt%dt%fn”, t, dsd, dsm);
if(fabs(dsm)> eps){
kg[kq0]=1;if(kq0 == 0)
calc(iu1, u1, di1, nfd1, b);if(kq0 == 1)
calc(iu2, u2, di2, nfd2, b);for(i = 1;i <= n;i ++){
if(kq0 == 0)
va[i] = va[i]b[i];} } else kg[kq0] = 0;if(kq0 == 0)
kq0 = 1;else
{
kq0 = 0;
t ++;} if(t > it1)
break;}while((fabs(dsm)> eps)||(kg[kq0]!= 0));
fprintf(fp2, “n%s%d”, “times = ”, t);}
void out()
{
//**** 本函數(shù)輸出節(jié)點和支路數(shù)據 ****//
zlsort(old);
// recover the data if sorted // newtoold();node_output();// node data branch_output();
//branch data printc('-', 78);printc('*', 78);fprintf(fp2, “n”);}
void newval(double* aa){
//**** 本函數(shù)將舊號換成新號 ****//
int i, k1;for(i = 1;i <= n;i ++)
b[i] = 0.0;for(i = 1;i <= nb;i ++)
{
k1 = nnew[nob[i]];
b[k1] = aa[i];
} for(i = 1;i <= n;i ++)
aa[i] = b[i];}
void yzb(int t, int* iu, double* u, double* di, int* nfd){ //**** 本函數(shù)求因子表 ****//
//參數(shù)1為標志(t=0 求B',t=1求B'')// //參數(shù)2因子表上三角矩陣非零非對角元素的列足碼 //參數(shù)3因子表上三角矩陣非零非對角元素的數(shù)值 //參數(shù)4因子表上三角矩陣對角元素 //參數(shù)5因子表上三角各行非零元素個數(shù)
int i, j, k, i1, i2;
int jj, jj1, jj2, im, x, fd[NS];double ai, b[NS];nfd[1] = 1;for(i = 1;i <= n;i ++){
//nobt[] 存放的是節(jié)點類型, 0: pq節(jié)點,-1: pv節(jié)點。if(((t!= 1)||(nobt[i]!=-1))&& i!= mpj)// <---| {
//
| for(j = i + 1;j <= n;j ++)
//
|
b[j] = 0.0;
//
|
b[i] = bii[i];
//
|
if((kk2 == 0)&&(t == 1)&&(nobt[i]!=-1))// 存在(t == 1)的情況，不多余。
b[i] = b[i] + af[1] * ql[i] / v0[i] / v0[i];//af[1] i1 = ydz[i];i2 = ydz[i + 1]1;for(im = 1;im <= i1;im ++){
jj1 = nfd[im];
jj2 = nfd[im + 1]ai * u[k];
}
break;
}
}
}
x = nfd[i];
di[i] = 1.0 / b[i];
ai = di[i];
k = 0;
i1 = i + 1;
for(j = i1;j <= n;j ++)
{
if(fabs(b[j])> 1.0e-15)
{
u[x] = b[j] * ai;
iu[x] = j;
k++;
x++;
}
}
fd[i] = k;
}
else
{
fd[i] = 0;
di[i] = 0.0;
}
nfd[i+1] = nfd[i] + fd[i];} fprintf(fp2, “n********U*********”);for(i = 1;i <= x;i ++){
if(i % 3 == 1)
fprintf(fp2, “n”);
fprintf(fp2, “%10.5f%8i”, u[i], iu[i]);} fprintf(fp2, “n********DI********”);printf1(di, n);}
void printf1(double* aa, int n){ //**** 本函數(shù)輸出aa[i],i=1-n ****// int i;for(i = 1;i <= n;i ++)
{
if(i % 5 == 1)
fprintf(fp2, “n”);
fprintf(fp2, “%9.5f”, aa[i]);
} fprintf(fp2, “nn”);}
void calc(int* iu, double* u, double* di, int* nfd, double* b){
//**** 本函數(shù)利用因子表解線形方程組。(P417圖F1-9)
double bi;int i, j, k, i1, i2;for(i = 1;i <= n;i ++)//前代過程。
{
bi = b[i];
i1 = nfd[i];
i2 = nfd[i + 1];
for(j = i1;j < i2;j ++)
{
k = iu[j];
b[k] = b[k]1;
for(j = i2;j >= i1;j--)
{
k = iu[j];
bi = bi1;i ++)
{
ip = i;
k1 = iabs(izl[i]);
k3 = iabs(jzl[i]);
for(j = i + 1;j <= zls;j ++)
{
k2 = iabs(izl[j]);
k4 = iabs(jzl[j]);
if(k2 < k1 ||(k2 == k1 && k4 < k3))
{
ip = j;
k1 = k2;
k3 = k4;
}
}
if(i!= ip)
{
iswap(&izl[i], &izl[ip]);
}
} }
void bnsopt(){ iswap(&jzl[i], &jzl[ip]);swap(&zr[i], &zr[ip]);swap(&zx[i], &zx[ip]);swap(&zyk[i], &zyk[ip]);//**** 節(jié)點優(yōu)化 ****//
int ii1, ii2, zls2, nomax;int i, j, l, k1, k;
int temp;zls2 = zls + zls;for(i = 1;i <= zls2;i ++)
old[i] = nnew[i] = 0;//先清零。由此可知：NS4、NS必須大于2*zls。
for(i = 1;i <= zls;i ++)
{
old[i] = iabs(izl[i]);
old[i + zls] = iabs(jzl[i]);} //變壓器節(jié)點號由正變負，old[]前zls個為左節(jié)點號，后zls個為右節(jié)點號。
for(i = 1;i <= zls2;i ++)// 冒泡法排序。
{
k1 = i + 1;
for(j = k1;j <= zls2;j ++)
if(old[i] > old[j])
iswap(&old[i], &old[j]);
//交換整數(shù)old[i]、old[j]。小節(jié)點號排在支路左側。
}
nomax = old[zls2];//nomax 即是最大節(jié)點號。Iee30.dat----30 l = 1;for(i = 1;i <= n;i ++)
{
ii1 = old[l];
old[i] = ii1;
for(j = l;j <= zls2;j ++)
{
ii2 = old[j];
if(ii1!= ii2)
{
l = j;
break;
}
nnew[i] ++;
} } for(i = 1;i <= n1;//總支路數(shù)(雙回線算一條支路)ydz[1] = 1;for(i = 1;i <= ngij;yb[ll] = yb[ll]gij / yk;yb[ll] = yb[ll]1;
//即為第 i 行的非零元素個數(shù)。
for(k = i1;k <= i2;k ++)//對第 i 行的所有非零元素進行操作。{
if(kq0!= 0)
{
ai =-yb[k];
bi = yg[k];}//yb[]存放導納陣非對角元的虛部(yg[] + j*yb[])else
{
ai = yg[k];
bi = yb[k];}//yg[]存放導納陣非對角元的實部(yg[] + j*yb[])j = iy[k];//iy[]存放的是列足碼。
vj = vi * v[j];// v[]存放節(jié)點電壓幅值。
ci = va[i]bi;//非對角元素 j 的功率
} }
void bbhl(int kq0)
{
//**** 本函數(shù)計算各節(jié)點的功率誤差，求最大功率誤差dsm ****// //**** 和常數(shù)項b[i]。(程序框圖見P423)****//
int i;double vi, vj;double pl0, pg0;dsm = 0.0;
// dsm 即為最大功率誤差。
for(i = 1;i <= n;i ++){
vi = v[i];// v[]存放節(jié)點電壓幅值。vj = v0[i];// v0[]存放節(jié)點初始電壓幅值。v0[]存放的是最后一個節(jié)點數(shù)據。// vi[] // nopt[] =-1: pv節(jié)點，對于pv節(jié)點，即為該點應維持的電壓值。// nopt[] = 0: pq節(jié)點，對于pq節(jié)點，即為系統(tǒng)平均電壓值。// vj 此時表示的是節(jié)點正常電壓的會定值。if(kq0 == 0){ pl0 = pl[i];pg0 = pg[i];
}
// 負荷p,發(fā)電機p
else
{
pl0 = ql[i];
pg0 = qg[i];
}
// 負荷q,發(fā)電機q if(kk2 == 0)
pl0 = pl0 *((vipl0;
// pv節(jié)點 if(i == mpj && kq0 == 0)
pg[i] = w[i]w[i];//pv節(jié)點(nobt[] =-1)和平衡節(jié)點(mpj)不參與 //求最大功率誤差和常數(shù)項的運算
if(((kq0!= 1)||(nobt[i]!=-1))&&(i!= mpj)){ if(fabs(b[i])> fabs(dsm))
{
dsm = b[i];
dsd = i;}
// dsm 即為最大功率誤差,dsd存放其對應的節(jié)點號。
b[i] = b[i] / vi;
// 計算修正方程式的常數(shù)項。} else
b[i]=0.0;
// pv節(jié)點(nobt[] =-1)和平衡節(jié)點(mpj)不參與
// 求最大功率誤差和常數(shù)項的運算。
} }
node_output()
//
{ //**** 輸出節(jié)點數(shù)據和最小電壓幅值、相角(角度)及其節(jié)點號 ****// //****(程序框圖見p426 F1-16)****//
double vi, ci;int i, j, oldnumber;printc('+', 72);fprintf(fp2, “n%5s%8s%10s%11s%11s%11s%11sn”, “I”，“V”, “CA”, “PL”, “QL”, “PG”, “QG”);vmin = v[1];dsd = 1;for(i = 1;i <= n;i ++){
j = nnew[i];
oldnumber = newsort[i];//轉化為相應舊號
ci = va[j] * 180.0/3.1416;//弧度轉化為角度。
vi = v[j];
if(vi < vmin)
{
vmin = vi;
dsd = j;
}
// vmin即為最小電壓,dsd存放其對應的新節(jié)點號。
fprintf(fp2, “n%5d%11.5f%12.6f”, oldnumber, vi, ci);
fprintf(fp2, “%11.5f%11.5f%11.5f%11.5f”, pl[j], ql[j]，pg[j], qg[j]);} printc('-', 72);} void branch_output()
{
//**** 本函數(shù)輸出支路數(shù)據。(程序框圖見p428 F1-17)
int ii, jj;double r, x, yk, zf, vi, vj, ci, cj;
int i, j, l;double de, df, ei, ej, fi, fj, fii, fir, pij, pji, qij, qji;
****// dph = 0.0;
//
統(tǒng)計系統(tǒng)有功網損。
dqh = 0.0;
//
統(tǒng)計系統(tǒng)無功網損。
fprintf(fp2, “n%5s%5s%10s%12s%12s%12sn ”, “I”, “J”, “PIJ”，“QIJ”, “PJI”, “QJI”);for(i =1;i <= mdk;i ++)
{
j = idk[i];
dkk[i] = v[j] * v[j] / dkk[i];
} for(l = 1;l <= zls;l ++)
{
ii = iabs(izl[l]);
// izl[]: 支路左節(jié)點號。jj = iabs(jzl[l]);
// jzl[]: 支路右節(jié)點號。i=nnew[ii];j=nnew[jj];
//
轉換為新節(jié)點號。ii = newsort[ii];jj = newsort[jj];//轉化為相應舊號 r = zr[l];x = zx[l];yk = zyk[l];vi = v[i];
//
v[]: 電壓幅值。ci = va[i];
//
va[]: 電壓相角。vj = v[j];cj = va[j];//支路左、右節(jié)點電壓值由極坐標轉換為直角坐標 ei = vi * cos(ci);fi = vi * sin(ci);
// ei: 支路左節(jié)點電壓實部，fi: 支路左節(jié)點電壓虛部。ej = vj * cos(cj);fj = vj * sin(cj);
// ej: 支路右節(jié)點電壓實部，fj: 支路右節(jié)點電壓虛部。if((izl[l] < 0)||(jzl[l] < 0))// 變壓器支路。{ if(izl[l] < 0)
{
ei = ei / yk;
fi = fi / yk;
}
// yk=zyk[l]
else
{
ej = ej / yk;
fj = fj / yk;
}
yk = 0.0;
}
de = eifj;
zf = r * r + x * x;
fii =(de * r + df * x)/ zf;
fir =(df * rfir * ei;
pji =-fii * ejvi * vi * yk;
qji = qji1;i ++)
{
ip = i;
k1 = iabs(izl[i]);
k3 = iabs(jzl[i]);
for(j = i + 1;j <= zls;j ++)
{
k2 = iabs(izl[j]);
k4 = iabs(jzl[j]);
if(k2 < k1 ||(k2 == k1 && k4 < k3))
{
ip = j;
k1 = k2;
k3 = k4;
}
}
if(i!= ip)
{
iswap(&izl[i], &izl[ip]);
iswap(&jzl[i], &jzl[ip]);
swap(&zr[i], &zr[ip]);
swap(&zx[i], &zx[ip]);
swap(&zyk[i], &zyk[ip]);
}
} for(i = 1;i <= nb;i ++){
k = nob[i];
nob[i] = newsort[k];
printf(“nob[%d] = %dn”, i, nob[i]);} for(i = 1;i <= nb-1;i ++){
for(j = i+1;j <= nb;j ++)
{
if(nob[i] > nob[j])
{
k = nob[i];
nob[i] = nob[j];
nob[j] = k;
}
} } printf(“nob[%d] = %dn”, i, nob[i]);}
void yy1()
{
//**** 本函數(shù)形成節(jié)點導納陣(不包括接地支路)****// int j1;double r, x, yk, zf, gij, bij;int i, j, i1, l, ll;
for(i = 1;i <= n;i ++)
{
gii[i] = 0.0;
bii[i] = 0.0;
}// 導納陣對角元(與節(jié)點一一對應)先清零。for(i = 1;i <= zls;i ++){
yg[i] = 0.0;
yb[i] = 0.0;} //導納陣非零非對角元(與支路一一對應)先清零。ll = 1;for(l = 1;l <= zls;l ++){
i1 = izl[l];j1 = jzl[l];i = iabs(i1);j = iabs(j1);if(i == j)
// 支路左節(jié)點號。// 支路右節(jié)點號。
// 變壓器支路有一節(jié)點號為負值。
continue;
//
排除左、右節(jié)點號相等的情況。
r = zr[l];x = zx[l];yk = zyk[l];
// zr[],zx[],zyk[]:支路三參數(shù)。
zf = r * r + x * x;gij = r / zf;//bij =-x / zf;bij =-1/x;yg[ll] = yg[ll]bij;gii[i] = gii[i] + gij;bii[i] = bii[i] + bij;gii[j] = gii[j] + gij;bii[j] = bii[j] + bij;
if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1])))
ll++;} // 打印導納矩陣。對角元實部為gii，虛部為bii，// 非零非對角元實部為yb[]，虛部為yb[]，列足碼為iy[]。
fprintf(fp2, “*******GII(1),BII(1)********n”);printf2(gii,bii,n);}
void y3()
{
//**** 本函數(shù)形成節(jié)點導納陣，追加接地支路 ****//
int j1;double r, x, yk, zf, gij, bij;int i, j, i1, l, ll, kk = 0;
for(i = 1;i <= mdk;i ++)
{
j = idk[i];
bii[j] = bii[j]gij / yk;
yb[ll] = yb[ll]-bij / yk;
}
if((i!= iabs(izl[l + 1]))||(j!= iabs(jzl[l + 1])))
{
ll++;
kk = 0;
}
else
kk = 1;} // 打印導納矩陣。對角元實部為gii，虛部為bii，// 非零非對角元實部為yb[]，虛部為yb[]，列足碼為iy[]。
fprintf(fp2, “*******GII,BII********”);printf2(gii,bii,n);
fprintf(fp2, “n*******YYYYY********”);for(i = 1;i <= nzls;i ++)
{
if(i % 2 == 1)
fprintf(fp2, “n”);
fprintf(fp2, “%10.4f%10.4f%8d”, yg[i], yb[i], iy[i]);
} }
第二篇：電力系統(tǒng)通用潮流計算C語言程序
#include
//節(jié)點號
類型
負荷有功
負荷無功母線數(shù)據（類型 1＝PV節(jié)點，2＝PQ節(jié)點，3＝平衡節(jié)點）struct BUS { int busno;int type;float Pd;float Qd;};
//發(fā)電機數(shù)據 節(jié)點號
有功發(fā)電
電壓幅值 struct Generator { int busno;float Pg;float Vg;};
//支路信息 節(jié)點I 節(jié)點J R X B/2 k struct Line { int busi;int busj;float R;float X;float B;float k;};//deltaP deltaQ deltaV^2// void fun1(double YG[][50],double YB[][50],double e[],double f[],int type[],int N,double W[],double P[],double Q[],double V[]){ double dP=0,dQ=0,dV=0;int i,j;for(i=0;i double A=0,B=0; for(j=0;j A+=YG[i][j]*e[j]-YB[i][j]*f[j]; B+=YG[i][j]*f[j]+YB[i][j]*e[j]; } dV=V[i]*V[i]-e[i]*e[i]-f[i]*f[i]; dP=P[i]-e[i]*A-f[i]*B; W[2*i]=dP; dQ=Q[i]-f[i]*A+e[i]*B; if(type[i]==1) W[2*i+1]=dQ; else W[2*i+1]=dV;} } //Jacobi矩陣// void Jacobi(double YG[][50],double YB[][50],double e[50],double f[50],int type[50],int N ,double Ja[100][101]){ int i,j; for(i=0;i for(j=0;j if(i!=j){ if(type[i]==1){ Ja[2*i][2*j]=-(YG[i][j]*e[i]+YB[i][j]*f[i]); Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i]; Ja[2*i+1][2*j]=Ja[2*i][2*j+1]; Ja[2*i+1][2*j+1]=-Ja[2*i][2*j]; } else { Ja[2*i][2*j]=-YG[i][j]*e[i]+YB[i][j]*f[i]; Ja[2*i][2*j+1]=YB[i][j]*e[i]-YG[i][j]*f[i]; Ja[2*i+1][2*j+1]=Ja[2*i+1][2*j]=0; } } else { double a[50]={0},b[50]={0}; for(int k=0;k a[i]+=(YG[i][k]*e[k]-YB[i][k]*f[k]); b[i]+=(YG[i][k]*f[k]+YB[i][k]*e[k]); Ja[2*i][2*j]=-a[i]-YG[i][i]*e[i]-YB[i][i]*f[i]; Ja[2*i][2*j+1]=-b[i]+YB[i][i]*e[i]-YG[i][i]*f[i]; if(type[i]==1){ Ja[2*i+1][2*j]=b[i]+YB[i][i]*e[i]-YG[i][i]*f[i]; Ja[2*i+1][2*j+1]=-a[i]+YG[i][i]*e[i]+YB[i][i]*f[i]; } else { Ja[2*i+1][2*j]=-2*e[i]; Ja[2*i+1][2*j+1]=-2*f[i]; } } } } } } //高斯消元法解方程組函數(shù)// void gauss(double a[][101],int n) { int i,j,k;double c; for(k=0;k c=a[k][k]; for(j=k;j<=n;j++)a[k][j]/=c; for(i=k+1;i c=a[i][k]; for(j=k;j<=n;j++)a[i][j]-=c*a[k][j]; } } a[n-1][n]/=a[n-1][n-1]; for(k=n-2;k>=0;k--) for(j=k+1;j void main(){ ifstream fin;int N=0,GS=0,LD=0,ZLs=0;//節(jié)點數(shù) 發(fā)電機數(shù) 負荷數(shù) BUS *B; Generator *G;Line *L; //從文本中讀入原始數(shù)據到數(shù)組中// fin.open(“C:data.txt”);if(!fin) { cout<<“輸入數(shù)據文件不存在！”< getchar(); } int m1[50]={0},m2[50]={0};float m3[50],m4[50],m5[50],m6[50]; int i,j,l; 支路數(shù)// for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m2[i]>>m3[i]>>m4[i]; N++;} B =new BUS[N];for(i=0;i B[i].busno=m1[i]; B[i].type=m2[i]; B[i].Pd=m3[i]; B[i].Qd=m4[i];} for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m4[i]>>m3[i]; GS++;} G =new Generator[GS];for(i=0;i G[i].busno=m1[i]; G[i].Pg=m4[i]; G[i].Vg=m3[i];} for(i=0;;i++){ fin>>m1[i]; if(m1[i]==0)break; fin>>m2[i]>>m3[i]>>m4[i]>>m5[i]>>m6[i]; ZLs++;} L =new Line[ZLs];for(i=0;i L[i].busi=m1[i]; L[i].busj=m2[i]; L[i].R=m3[i]; L[i].X=m4[i]; } L[i].B=m5[i];L[i].k=m6[i]; LD=N-GS;fin.close(); //節(jié)點導納矩陣形成// double YB[50][50],YG[50][50],BB[50][50],K[50][50];for(i=0;i YB[i][j]=0;YG[i][j]=0;BB[i][j]=0;K[i][j]=1;} } for(l=0;l } for(i=0;i for(j=i;j K[i][j]=K[j][i];K[j][i]=1; } for(j=0;j if(i!=j){ YG[i][i]=YG[i][i]+(YG[i][j]*K[i][j]*K[i][j]); YB[i][i]=YB[i][i]+(YB[i][j]*K[i][j]*K[i][j]+BB[i][j]); } } } //修正后// for(l=0;l // 求A=e+f// double e[50]={0},f[50]={0};double C[100]={0},D[100]={0};for(i=0;i C[2*i]=1;} else C[2*i]=V[i]; } double W[100]={0},Ja[100][101]={0}; //調用Jacobi函數(shù)和高斯函數(shù)// for(int t=1;t<10;t++){ for(i=0;i<2*N-2;i++){ e[i]=C[2*i]; f[i]=C[2*i+1]; } fun1(YG,YB,e,f,type,N,W,P,Q,V); double it=fabs(W[0]); for(i=1;i<2*N-2;i++){ if(it it=fabs(W[i]);j=i; } } //中間迭代過程// cout< cout< if(it<0.00001)break; Jacobi(YG,YB,e,f,type,N,Ja); for(i=0;i<2*N-2;i++){ Ja[i][2*N-2]=W[i]; } //高斯消元法解方程// gauss(Ja,2*N-2); for(i=0;i<2*N-2;i++){ D[i]=-Ja[i][2*(N-1)]; C[i]+=D[i]; } } //平衡節(jié)點// for(i=0;i double a=0,b=0; for(int j=0;j a+=(YG[i][j]*e[j]-YB[i][j]*f[j]); b+=(YB[i][j]*e[j]+YG[i][j]*f[j]); } P[i]=e[i]*a+f[i]*b; Q[i]=f[i]*a-e[i]*b; } //支路// double PZL[100][101]={0},QZL[100][101]={0},pr[100][101]={0},qx[100][101]={0};double x1=0,x2=0,y1=0,y2=0,I2=0;for(int k=0;k i=L[k].busi-1;j=L[k].busj-1;x1=e[i]/L[k].k-e[j];y1=f[i]/L[k].k-f[j];x2=-e[i]*YG[i][j]-f[i]*YB[i][j];y2=-f[i]*YG[i][j]+e[i]*YB[i][j];QZL[i][j]=(x1*y2-x2*y1);PZL[i][j]=(x1*x2+y1*y2);I2=(PZL[i][j]*PZL[i][j]+QZL[i][j]*QZL[i][j])/(e[i]*e[i]+f[i]*f[i]);pr[i][j]=I2*L[k].R;qx[i][j]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[i][j]+=(e[i]*e[i]+f[i]*f[i])*(-L[k].B);x1=e[j]*L[k].k-e[i];y1=f[j]*L[k].k-f[i];x2=-e[j]*YG[j][i]-f[j]*YB[j][i];y2=-f[j]*YG[j][i]+e[j]*YB[j][i];QZL[j][i]=(x1*y2-x2*y1);PZL[j][i]=(x1*x2+y1*y2);I2=(PZL[j][i]*PZL[j][i]+QZL[j][i]*QZL[j][i])/(e[j]*e[j]+f[j]*f[j]);pr[j][i]=I2*L[k].R;qx[j][i]=I2*L[k].X-(e[i]*e[i]+f[i]*f[i]+e[j]*e[j]+f[j]*f[j])*L[k].B;QZL[j][i]+=(e[j]*e[j]+f[j]*f[j])*(-L[k].B); } //全網數(shù)據// int high=1,low=1; double PG=0,PL=0,Prr=0,Vh=sqrt(e[0]*e[0]+f[0]*f[0]),Vl=sqrt(e[0]*e[0]+f[0]*f[0]);for(k=0;k Vh=sqrt(e[k]*e[k]+f[k]*f[k]);high=k+1;} if(sqrt(e[k]*e[k]+f[k]*f[k]) Vl=sqrt(e[k]*e[k]+f[k]*f[k]);low=k+1;} } //輸出數(shù)據到文件databak.txt// ofstream fout; fout.open(“C:databak.txt”); fout<<“節(jié)點”< fout< for(i=0;i fout< } for(j=0;j i=G[j].busno-1; fout< fout<<“支路 ”< fout< for(k=0;k i=L[k].busi-1;j=L[k].busj-1; fout< fout<<“全網數(shù)據”< fout< fout< 南 京 理 工 大 學 《電力系統(tǒng)穩(wěn)態(tài)分析》 課程報告 姓名 XX 學 號： 5*** 自動化學院 電氣工程 基于牛頓-拉夫遜法的潮流計算例題編程報學院(系): 專 業(yè): 題 目: 任課教師 碩士導師 告 楊偉 XX 2015年6月10號 基于牛頓-拉夫遜法的潮流計算例題編程報告 摘要：電力系統(tǒng)潮流計算的目的在于：確定電力系統(tǒng)的運行方式、檢查系統(tǒng)中各元件是否過壓或者過載、為電力系統(tǒng)繼電保護的整定提供依據、為電力系統(tǒng)的穩(wěn)定計算提供初值、為電力系統(tǒng)規(guī)劃和經濟運行提供分析的基礎。潮流計算的計算機算法包含高斯—賽德爾迭代法、牛頓-拉夫遜法和P—Q分解法等，其中牛拉法計算原理較簡單、計算過程也不復雜，而且由于人們引入泰勒級數(shù)和非線性代數(shù)方程等在算法里從而進一步提高了算法的收斂性和計算速度。同時基于MATLAB的計算機算法以雙精度類型進行數(shù)據的存儲和運算, 數(shù)據精確度高，能進行潮流計算中的各種矩陣運算，使得傳統(tǒng)潮流計算方法更加優(yōu)化。 一 研究內容 通過一道例題來認真分析牛頓-拉夫遜法的原理和方法（采用極坐標形式的牛拉法），同時掌握潮流計算計算機算法的相關知識，能看懂并初步使用MATLAB軟件進行編程，培養(yǎng)自己電力系統(tǒng)潮流計算機算法編程能力。 例題如下：用牛頓-拉夫遜法計算下圖所示系統(tǒng)的潮流分布，其中系統(tǒng)中5為平衡節(jié)點，節(jié)點5電壓保持U=1.05為定值，其他四個節(jié)點分別為PQ節(jié)點，給定的注入功率如圖所示。計算精度要求各節(jié)點電壓修正量不大于10-6。 二 牛頓-拉夫遜法潮流計算 1 基本原理 牛頓法是取近似解x(k)之后，在這個基礎上，找到比x(k)更接近的方程的根，一步步地迭代，找到盡可能接近方程根的近似根。牛頓迭代法其最大優(yōu)點是在方程f(x)=0的單根附近時誤差將呈平方減少，而且該法還可以用來求方程的重根、復根。電力系統(tǒng)潮流計算，一般來說，各個母線所供負荷的功率是已知的，各個節(jié)點的電壓是未知的（平衡節(jié)點外）可以根據網絡結構形成節(jié)點導納矩陣，然后由節(jié)點導納矩陣列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計算的問題就轉化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導，得出對應的雅可比矩陣，給未知節(jié)點賦電壓初值，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節(jié)點電壓不平衡量（未知的）構成了誤差方程，解誤差方程，得到節(jié)點電壓不平衡量，節(jié)點電壓加上節(jié)點電壓不平衡量構成節(jié)點電壓新的初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。2 基本步驟和設計流程圖 形成了雅克比矩陣并建立了修正方程式，運用牛頓-拉夫遜法計算潮流的核心問題已經解決，已有可能列出基本計算步驟并編制流程圖。由課本總結基本步驟如下： 1)形成節(jié)點導納矩陣Y； 2)設各節(jié)點電壓的初值，如果是直角坐標的話設電壓的實部e和虛部f；如果是極坐標的話則設電壓的幅值U和相角a； 3)將各個節(jié)點電壓的初值代入公式求修正方程中的不平衡量以及修正方程的系數(shù)矩陣的雅克比矩陣； 4)解修正方程式，求各節(jié)點電壓的變化量，即修正量； 5)計算各個節(jié)點電壓的新值，即修正后的值； 6)利用新值從第（3）步開始進入下一次迭代，直至達到精度退出循環(huán)； 7)計算平衡節(jié)點的功率和線路功率，輸出最后計算結果； ① 公式推導 ② 流程圖 三 matlab編程代碼 clear; % 如圖所示1，2，3，4為PQ節(jié)點,5為平衡節(jié)點 y=0; % 輸入原始數(shù)據，求節(jié)點導納矩陣 y(1,2)=1/(0.07+0.21j); y(4,5)=0;y(1,3)=1/(0.06+0.18j); y(1,4)=1/(0.05+0.10j); y(1,5)=1/(0.04+0.12j); y(2,3)=1/(0.05+0.10j); y(2,5)=1/(0.08+0.24j); y(3,4)=1/(0.06+0.18j); for i=1:5 for j=i:5 y(j,i)=y(i,j); end end Y=0; % 求節(jié)點導納矩陣中互導納 for i=1:5 for j=1:5 if i~=j Y(i,j)=-y(i,j); end end end % 求節(jié)點導納矩陣中自導納 for i=1:5 Y(i,i)=sum(y(i,:)); end Y % Y為導納矩陣 G=real(Y); B=imag(Y);% 輸入原始節(jié)點的給定注入功率 S(1)=0.3+0.3j; S(2)=-0.5-0.15j; S(3)=-0.6-0.25j; S(4)=-0.7-0.2j; S(5)=0; P=real(S); Q=imag(S); % 賦初值，U為節(jié)點電壓的幅值，a為節(jié)點電壓的相位角 U=ones(1,5); U(5)=1.05; a=zeros(1,5); x1=ones(8,1); x2=ones(8,1); k=0; while max(x2)>1e-6 for i=1:4 for j=1:4 H(i,j)=0; N(i,j)=0; M(i,j)=0; L(i,j)=0; oP(i)=0; oQ(i)=0; end end % 求有功、無功功率不平衡量 for i=1:4 for j=1:5 oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j))); end oP(i)=oP(i)+P(i); oQ(i)=oQ(i)+Q(i); end x2=[oP,oQ]'; % x2為不平衡量列向量 % 求雅克比矩陣 % 當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(a(i)-a(j))-B(i,j)*cos(a(i)-a(j))); N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); L(i,j)=H(i,j); M(i,j)=-N(i,j); end end end % 當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(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j))); L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(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); end J=[H,N;M,L] % J為雅克比矩陣 x1=-((inv(J))*x2); % x1為所求△x的列向量 % 求節(jié)點電壓新值，準備下一次迭代 for i=1:4 oa(i)=x1(i); oU(i)=x1(i+4)*U(i); end for i=1:4 a(i)=a(i)+oa(i); U(i)=U(i)+oU(i); end k=k+1; end k,U,a % 求節(jié)點注入功率 i=5; for j=1:5 P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i); Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i); end S(5)=P(5)+Q(5)*sqrt(-1); S % 求節(jié)點注入電流 I=Y*U' 四 運行結果 節(jié)點導納矩陣 經過五次迭代后的雅克比矩陣 迭代次數(shù)以及節(jié)點電壓的幅值和相角（弧度數(shù)） 節(jié)點注入功率和電流 五 結果分析 在這次學習和實際操作過程里：首先，對電力系統(tǒng)分析中潮流計算的部分特別是潮流計算的計算機算法中的牛頓-拉夫遜法進行深入的研讀，弄明白了其原理、計算過程、公式推導以及設計流程。牛頓-拉夫遜法是求解非線性方程的迭代過程，其計算公式為?F?J?X，式中J為所求函數(shù)的雅可比矩陣；?X為需要求的修正值；?F為不平衡的列向量。利用x(*)=x(k+1)+?X(k+1)進行多次迭代，通過迭代判據得到所需要的精度值即準確值x(*)。六 結論 通過這個任務，自己在matlab編程，潮流計算，word文檔的編輯功能等方面均有提高，但也暴漏出一些問題：理論知識儲備不足，對matlab的性能和特點還不能有一個全面的把握，對word軟件也不是很熟練，相信通過以后的學習能彌補這些不足，達到一個新的層次。 電力系統(tǒng)潮流計算發(fā)展史 對潮流計算的要求可以歸納為下面幾點： （1）算法的可靠性或收斂性（2）計算速度和內存占用量（3）計算的方便性和靈活性 電力系統(tǒng)潮流計算屬于穩(wěn)態(tài)分析范疇，不涉及系統(tǒng)元件的動態(tài)特性和過渡過程。因此其數(shù)學模型不包含微分方程，是一組高階非線性方程。非線性代數(shù)方程組的解法離不開迭代，因此，潮流計算方法首先要求它是能可靠的收斂，并給出正確答案。隨著電力系統(tǒng)規(guī)模的不斷擴大，潮流問題的方程式階數(shù)越來越高，目前已達到幾千階甚至上萬階，對這樣規(guī)模的方程式并不是采用任何數(shù)學方法都能保證給出正確答案的。這種情況促使電力系統(tǒng)的研究人員不斷尋求新的更可靠的計算方法。 在用數(shù)字計算機求解電力系統(tǒng)潮流問題的開始階段，人們普遍采用以節(jié)點導納矩陣為基礎的高斯-賽德爾迭代法（一下簡稱導納法）。這個方法的原理比較簡單，要求的數(shù)字計算機的內存量也比較小，適應當時的電子數(shù)字計算機制作水平和電力系統(tǒng)理論水平，于是電力系統(tǒng)計算人員轉向以阻抗矩陣為主的逐次代入法（以下簡稱阻抗法）。 20世紀60年代初，數(shù)字計算機已經發(fā)展到第二代，計算機的內存和計算速度發(fā)生了很大的飛躍，從而為阻抗法的采用創(chuàng)造了條件。阻抗矩陣是滿矩陣，阻抗法要求計算機儲存表征系統(tǒng)接線和參數(shù)的阻抗矩陣。這就需要較大的內存量。而且阻抗法每迭代一次都要求順次取阻抗矩陣中的每一個元素進行計算，因此，每次迭代的計算量很大。 阻抗法改善了電力系統(tǒng)潮流計算問題的收斂性，解決了導納法無法解決的一些系統(tǒng)的潮流計算，在當時獲得了廣泛的應用，曾為我國電力系統(tǒng)設計、運行和研究作出了很大的貢獻。但是，阻抗法的主要缺點就是占用計算機的內存很大，每次迭代的計算量很大。當系統(tǒng)不斷擴大時，這些缺點就更加突出。為了克服阻抗法在內存和速度方面的缺點，后來發(fā)展了以阻抗矩陣為基礎的分塊阻抗法。這個方法把一個大系統(tǒng)分割為幾個小的地區(qū)系統(tǒng)，在計算機內只需存儲各個地區(qū)系統(tǒng)的阻抗矩陣及它們之間的聯(lián)絡線的阻抗，這樣不僅大幅度的節(jié)省了內存容量，同時也提高了節(jié)省速度。 克服阻抗法缺點的另一途徑是采用牛頓-拉夫遜法（以下簡稱牛頓法）。牛頓法是數(shù)學中求解非線性方程式的典型方法，有較好的收斂性。解決電力系統(tǒng)潮流計算問題是以導納矩陣為基礎的，因此，只要在迭代過程中盡可能保持方程式系數(shù)矩陣的稀疏性，就可以大大提高牛頓潮流程序的計算效率。自從20世紀60年代中期采用了最佳順序消去法以后，牛頓法在收斂性、內存要求、計算速度方面都超過了阻抗法，成為直到目前仍被廣泛采用的方法。 在牛頓法的基礎上，根據電力系統(tǒng)的特點，抓住主要矛盾，對純數(shù)學的牛頓法進行了改造，得到了P-Q分解法。P-Q分解法在計算速度方面有顯著的提高，迅速得到了推廣。 牛頓法的特點是將非線性方程線性化。20世紀70年代后期，有人提出采用更精確的模型，即將泰勒級數(shù)的高階項也包括進來，希望以此提高算法的性能，這便產生了保留非線性的潮流算法。另外，為了解決病態(tài)潮流計算，出現(xiàn)了將潮流計算表示為一個無約束非線性規(guī)劃問題的模型，即非線性規(guī)劃潮流算法。 近20多年來，潮流算法的研究仍然非常活躍，但是大多數(shù)研究都是圍繞改進牛頓法和P-Q分解法進行的。此外，隨著人工智能理論的發(fā)展，遺傳算法、人工神經網絡、模糊算法也逐漸被引入潮流計算。但是，到目前為止這些新的模型和算法還不能取代牛頓法和P-Q分解法的地位。由于電力系統(tǒng)規(guī)模的不斷擴大，對計算速度的要求不斷提高，計算機的并行計算技術也將在潮流計算中得到廣泛的應用，成為重要的研究領域。 電力系統(tǒng)潮流計算程序設計 姓名：韋應順 學號：2011021052 電力工程學院 牛頓—拉夫遜潮流計算方法具有能夠將非線性方程線性化的特點，而使用MATLAB語言是由于MATLAB語言的數(shù)學邏輯強，易編譯。 【】【】1.MATLAB程序12 Function tisco %這是一個電力系統(tǒng)潮流計算的程序 n=input（‘n請輸入節(jié)點數(shù)：n=’）； m=input(‘請輸入支路數(shù)：m=’);ph=input(‘n請輸入平衡母線的節(jié)點號：ph=’)； B1=input（‘n請輸入支路信號：B1=’）;%它以矩陣形式存貯支路的情況，每行存貯一條支路 %第一列存貯支路的一個端點 %第二列存貯支路的另一個端點 %第三列存貯支路阻抗 %第四列存貯支路的對地導納 %第五列存貯變壓器的變比，注意支路為1 %第六列存貯支路的序號 B2=input(‘n請輸入節(jié)點信息：B2=’)； %第一列為電源側的功率 %第二列為負荷側的功率 %第三列為該點的電壓值 %第四列為該點的類型：1為PQ，2為PV節(jié)點，3為平衡節(jié)點 A=input（‘n請輸入節(jié)點號及對地阻抗：A=’）； ip=input（‘n請輸入修正值：ip=’）； %ip為修正值);Y=zeros（n）； Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）； Y(p,q)=Y(p,q);f=zeros（1，n）； no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2； End for i=1：n G=real(Y);if A（i2）=0 B=imag(Y);p=A(i1)； Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while max(abs(De))>ip&max(abs(Dfi)>ip t=t+1; e=e+De； f=f+Df; [C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)； J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)； [De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v); jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3)); b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2; c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%該子程序是用來求取Df for i=1:n If i=ph C(i)=0; D(i)=0; For j=i:n C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If B2(i4)=2 p=2*i-1; Df(p)=P(i)-P1;p=p+1;else p=2*i-1; Df(p)=P(i)-P1;p=p+1; Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no?=[];end function [De,Df]=hxf(J,Df,ph,n,no)%該子函數(shù)是為求取De Df DX=JDf;DX1=DX; x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0; For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else DX=[DX1,0,0];End k=0; [x,y]=size(DX);For i=1:2:x K=k+1; Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i); X2=G(i,j)*e(i)+B(i,j)*f(i);%該子程序是用來求取jacci矩陣 for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End 2實例驗證 【例題】設有一系統(tǒng)網絡結線見圖1，各支路阻抗和各節(jié)點功率均已以標幺值標示于圖1中，其中節(jié)點2連接的是發(fā)電廠，設節(jié)點1電壓保持U1＝1.06定值，試計算其中的潮流分布，請輸入節(jié)點數(shù)：n=5 請輸入支路數(shù)：m=7 請輸入平衡母線的節(jié)點號：ph=l 請輸入支路信息： BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 請輸入節(jié)點信息： B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 請輸入節(jié)點號及對地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 請輸入修正值：ip=0．000 0l 參考文獻 [1]陳珩．電力系統(tǒng)穩(wěn)定分析[M]．北京：中國電力出版社，2002：139—187． [2]鄭阿奇．MATLAB實用教程[M]．北京：電子工業(yè)出版社，2005：1-243． [3] 束洪春，孫士云，等．云電送粵交商流混聯(lián)系統(tǒng)全過 程動態(tài)電壓研究[J】．中國電力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn鋤ic Voltage stabil時of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4． [4] 朱新立，湯涌，等．大電網安全分析的全過程動態(tài)仿 真技術[J】．電網技術，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn鋤ic simulation for the stabilhy a眥lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28． [5]Roytelm鋤I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．發(fā)電機勵磁系統(tǒng)數(shù)學模型及參 數(shù)對電網動態(tài)穩(wěn)定性分析結果影響的研究[J】．繼電 器，2007，35(21)：22-27． SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g鋤e翰to璐baScd ∞de詛iled excitation system models柚d parameterS t0 power鏟id dyn鋤ic stabil時【J】．Relay，2007，35(2 1)：22-27． [7] 方思立，朱方．快速勵磁系統(tǒng)對系統(tǒng)穩(wěn)定的影響[J】．中 國電機工程學報，1986，6(1)：20．28． FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞ excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28． [8] 劉?。娏ο到y(tǒng)穩(wěn)定性及發(fā)電機勵磁控制[M】．北京： 中國電力出版社，2007． LIU Qu．Power system S詛bility鋤d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro詛ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996． [10] 陳利芳，陳天祿．淺談自并勵勵磁系統(tǒng)在大容量機組 中的應用【J】．繼電器，2007，35(1)：8l培4． CHEN Li-f抽島CHEN Tian—lIL Application of 辯l仁exci組tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84． [11] 方思立，劉增煌，孟慶和．大型汽輪發(fā)電機自并勵勵 磁系統(tǒng)的應用條件【J]．中國電力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu鋤g，MENG Qin爭hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63． [12]梁小冰，黃方能．利用EMTDC進行長持續(xù)時間過程 的仿真研究【J】．電網技術，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan爭眥ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57． [13]王卉，陳楷，彭哲，等．數(shù)字仿真技術在電力系統(tǒng)中 的應用及常用的幾種數(shù)字仿真工具【J】．繼電器，2004，32(21)：7l一75． wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眥hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75． [14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI時studies【s】．第三篇：電力系統(tǒng)潮流計算
第四篇：電力系統(tǒng)潮流計算發(fā)展史
第五篇：電力系統(tǒng)潮流計算程序設計
點擊咨詢 