📄 pq_powerflow.c
字号:
{
i1 = ydz[i];
//ydz[i]是第 i 行第一个元素的首地址，
//既在所有非零非对角元素中的顺序号。
i2 = ydz[i + 1] - 1; //即为第 i 行的非零元素个数。
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[]存放节点电压幅值。
ci = va[i] - va[j];
//va[]存放节点电压相角，ci 即为节点 i 和
//节点 j 之间的相角差(i-j)。
sn = sin(ci); //sin()的近似计算式
cs = cos(ci);
//cos()的近似计算式
bi = bi * vj * sn;
ai = ai * vj * cs;
w[i] = w[i] + ai + bi; //非对角元素 i 的功率
w[j] = w[j] + ai - bi; //非对角元素 j 的功率
}
}
}
}
void bbhl(int kq0)
{
//**** 本函数计算各节点的功率误差，求最大功率误差dsm ****//
//**** 和常数项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[]存放节点电压幅值。
vj = v0[i];
// v0[]存放节点初始电压幅值。v0[]存放的是最后一个节点数据。
// vi[]
// nopt[] = -1: pv节点，对于pv节点，即为该点应维持的电压值。
// nopt[] = 0: pq节点，对于pq节点，即为系统平均电压值。
// vj 此时表示的是节点正常电压的会定值。
if (kq0 == 0)
{
pl0 = pl[i];
pg0 = pg[i];
} // 负荷p,发电机p
else
{
pl0 = ql[i];
pg0 = qg[i];
} // 负荷q,发电机q
if (kk2 == 0)
pl0 = pl0 * ((vi - vj) * af[kq0] / vj + 1.0);
// 考虑负荷静特性。(公式F1-39,P423)
if (nobt[i] == -1 && kq0 == 1)
qg[i] = w[i] - pl0; // pv节点
if(i == mpj && kq0 == 0)
pg[i] = w[i] - pl0; // 平衡节点
b[i] = pg0 + pl0 - w[i];
//pv节点(nobt[] = -1)和平衡节点(mpj)不参与
//求最大功率误差和常数项的运算
if (((kq0 != 1) || (nobt[i] != -1)) && (i != mpj))
{
if (fabs(b[i]) > fabs(dsm))
{
dsm = b[i];
dsd = i;
} // dsm 即为最大功率误差,dsd存放其对应的节点号。
b[i] = b[i] / vi; // 计算修正方程式的常数项。
}
else
b[i]=0.0;
// pv节点(nobt[] = -1)和平衡节点(mpj)不参与
// 求最大功率误差和常数项的运算。
}
}
node_output() //
{
//**** 输出节点数据和最小电压幅值、相角(角度)及其节点号 ****//
//**** (程序框图见p426 F1-16) ****//
double vi, ci;
int i, j, oldnumber;
printc('+', 72);
fprintf(fp2, "\n%5s%8s%10s%11s%11s%11s%11s\n", "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存放其对应的新节点号。
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()
{
//**** 本函数输出支路数据。(程序框图见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; // 统计系统有功网损。
dqh = 0.0; // 统计系统无功网损。
fprintf(fp2, "\n%5s%5s%10s%12s%12s%12s\n ", "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[]: 支路左节点号。
jj = iabs(jzl[l]); // jzl[]: 支路右节点号。
i=nnew[ii];
j=nnew[jj]; // 转换为新节点号。
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];
//支路左、右节点电压值由极坐标转换为直角坐标
ei = vi * cos(ci);
fi = vi * sin(ci);
// ei: 支路左节点电压实部，fi: 支路左节点电压虚部。
ej = vj * cos(cj);
fj = vj * sin(cj);
// ej: 支路右节点电压实部，fj: 支路右节点电压虚部。
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 = ei - ej;
df = fi - fj;
zf = r * r + x * x;
fii = (de * r + df * x) / zf;
fir = (df * r - de * x) / zf;
pij = fii * ei + fir * fi;
qij = fii * fi - fir * ei;
pji = -fii * ej - fir * fj;
qji = -fii * fj + fir * ej;
qij = qij - vi * vi * yk;
qji = qji - vj * vj * yk;
dph = dph + pij + pji;
dqh = dqh + qij + qji;
fprintf(fp2, "\n%5d%5d%12.5f%12.5f%12.5f%12.5f", ii, jj, pij,
qij, pji, qji);
}
fprintf(fp2, "\n\n%8s%19s%18s%13s", "DPH", "DQH", "VMIN", "DSD");
fprintf(fp2, "\n %10.5f%18.5f%18.5f%12d", dph, dqh, vmin, newsort[old[dsd]]);
fprintf(fp2, "\n \n%22sTHE END OF OUTPUT", " ");
}
void iswap(int* m, int* n)
{
//**** 本函数交换m和n值 ****//
int p;
p = *m;
*m = *n;
*n = p;
}
void swap(double* m, double* n)
{
//**** 本函数交换m和n值 ****//
double p;
p = *m;
*m = *n;
*n = p;
}
int iabs(int a)
{
//**** 本函数求绝对值 ****//
a = (a >= 0) ? a : -a;
return a;
}
void printc(char aa, int n)
{
//**** 本函数输出n个aa字符 ****//
int i;
fprintf(fp2, "\n");
for (i = 1; i <= n; i ++)
fprintf(fp2, "%c", aa);
}
void newtoold()
{
int i, j, k, ii1, ii2, zls2, k1, k2, k3, k4, ip;
zls2 = zls + zls;
ii1 = 0;
for (i = 1; i <= zls; i ++)
{
k = izl[i];
ii2 = newsort[(int)fabs(k)];
izl[i] = isgn(ii2, izl[i]);
k = jzl[i];
ii2 = newsort[(int)fabs(k)];
jzl[i] = isgn(ii2, jzl[i]);
printf("izl[%d] = %d, jzl[%d] = %d\n", i, izl[i], i, jzl[i]);
}
for (i = 1; i <= zls - 1; 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] = %d\n", 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] = %d\n", i, nob[i]);
}
void yy1()
{
//**** 本函数形成节点导纳阵(不包括接地支路) ****//
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;
}// 导纳阵对角元(与节点一一对应)先清零。
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)
continue; // 排除左、右节点号相等的情况。
r = zr[l];
x = zx[l];
yk = zyk[l]; // zr[],zx[],zyk[]:支路三参数。
zf = r * r + x * x;
gij = r / zf;
//bij = -x / zf;
bij = -1/x;
yg[ll] = yg[ll] - gij;
yb[ll] = yb[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()
{
//**** 本函数形成节点导纳阵，追加接地支路 ****//
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] - 1.0 / dkk[i];
}//计算接地支路导纳，只影响导纳阵对角元(自导纳)。
ll = 1;
for (l = 1; l <= zls; l ++)
{
i1 = izl[l]; // 支路左节点号。
j1 = jzl[l]; // 支路右节点号。
i = iabs(i1); // 变压器支路有一节点号为负值。
j = iabs(j1);
if (i == j)
continue; // 排除左、右节点号相等的情况。
yk = zyk[l]; // zr[],zx[],zyk[]:支路三参数。
if ((i1 > 0) && (j1 > 0))// 不是变压器支路。是一般支路。
{
bii[i] = bii[i] + bij + yk;
bii[j] = bii[j] + bij + yk;
}
else // 变压器支路。
{
r = zr[l];
x = zx[l];
zf = r * r + x * x;
gij = r / zf;
bij = -x / zf;
if(j1 < 0)
{
i=iabs(j1);
j=iabs(i1);
}
// 若非标准变比在右(j)侧，则左、右互换，保证非标准变比在左侧。
if (kk == 0)
{
gii[i] = 0;
bii[i] = 0;
yg[ll] = 0;
yb[ll] = 0;
}
gii[i] = gii[i] + gij / yk / yk;
bii[i] = bii[i] + bij / yk / yk;
yg[ll] = yg[ll] - 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]);
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -