牛顿拉夫讯
谁能帮忙看看这个牛顿拉夫讯怎么改进的，能用于配电网潮流计算？牛顿拉夫讯
double PowerFlow(int iNode, int iBalc, ROAD *pRoad, NodePQ *pNodePQ, bool vOutputVoltage = false)
{
int M=2*(iNode-1);
int flag1; // PV节点标志
double g; // 电导
double b; // 电纳
double dvv; // 电压微分
double angle; // 相角
double p,dp; // 有功及其微分
double q,dq; // 无功及其微分
double e; // 节点电压实部
double f; // 节点电压虚部
double *jcbb = new double;
double *jcb = new double;
int i,j,i1,j1,n,m,t,k1,balc,node/*,flag,npv*/;
double rr,xx,bb,k,com1,com2,sump=0.0,sumq=0.0;
double lossp=0.0, lossq=0.0;
node = iNode; // 节点数
balc = iBalc - 1; // 平衡节点
for(i=0;i<node;i++){
flag1=0;
dp=0.0;
dq=0.0;
dvv=0.0;
}
for(i=0;i<node;i++) {
for(j=0;j<node;j++){
g=0.0;
b=0.0;
}
}
for(n = 0; n < node-1; n++){
i = pRoad.head;
j = pRoad.tail;
rr = pRoad.r;
xx = pRoad.x;
//printf("%d\t%d\t%.3f\t%.3f\n", pRoad.head, pRoad.tail, pRoad.r, pRoad.x);
bb = 0;
k = 1;
g+=-(rr/(rr*rr+xx*xx)/k);
b+=-(-xx/(rr*rr+xx*xx)/k);
g+=-(rr/(rr*rr+xx*xx)/k);
b+=-(-xx/(rr*rr+xx*xx)/k);
g+=rr/((rr*rr+xx*xx)*k*k);
g+=rr/(rr*rr+xx*xx);
b+=bb+(-xx/((rr*rr+xx*xx)*k*k));
b+=bb+(-xx/(rr*rr+xx*xx));
}
for(i=0;i<node;i++){
angle = 0;
v = 1;
p = pNodePQ.p;
q = pNodePQ.q;
//printf("P:%.3f\t Q:%.3f\n", pNodePQ.p, pNodePQ.q);
e=v*cos(angle*PI/180);
f=v*sin(angle*PI/180);
}
// flag1=1;
// printf("共有n=%d条支路\n",n);
// printf("节点导纳矩阵如下：\n");
// for(i=0;i<node;i++){
//for(j=0;j<node;j++)
// printf("%lf+j%lf ", g,b);
//printf("\n\n");}
//
// printf("初始值如下：\n");
// printf("节点号 电压 相角 有功 无功 e f\n");
// for(i=0;i<node;i++)
//printf("%d %lf%lf%lf%lf%lf%lf\n", (i+1), v, angle, p, q, e, f);
for(t=0;t<50;t++){
//printf("\n第%d次迭代开始：\n", t);
for(i=0;i<node;i++){
if(i==balc) /*如果是平衡节点*/
continue;
else{
com1=0.0, com2=0.0;
for(j=0;j<node;j++){
com1+=g*e-b*f;
com2+=g*f+b*e;
}
if(flag1!=1){ /*如果是PQ节点*/
dp=p-e*com1-f*com2;
dq=q-f*com1+e*com2;
}
else{
dp=p-e*com1-f*com2;
dvv=v*v-(e*e+f*f);
}
}
}
for(i=0,j=0;i<node;i++){
if(i==balc) continue;
else if(flag1!=1){
jcbb=dp;
jcbb=dq;
// printf("\ndp[%d]=%lf\ndq[%d]=%lf",(i+1),dp,(i+1),dq);
}
else{
jcbb=dp;
jcbb=dvv;
// printf("\ndp[%d]=%lf\ndvv[%d]=%lf",(i+1),dp,(i+1),dvv);
}
}
//printf("\n");
for(i=0;i<node;i++) {
if(fabs(dp)>=E||fabs(dq)>=E||fabs(dvv)>=E)
break;
}
if(i>=node)
break;
for(i=0,i1=0;i<node;i++){
if(i==balc)
continue;
i1++;
n=2*i1;
for(j=0,j1=0;j<node;j++){
if(j==balc) continue;
j1++;
m=j1*2;
if(i!=j){
jcb=-(g*e+b*f);
jcb=b*e-g*f;
if(flag1!=1)
{
jcb=jcb;
jcb=-jcb;
}
else
jcb=jcb=0.0;
}
else{
com1=0.0, com2=0.0;
for(k1=0;k1<node;k1++)
{
com1+=g*e-b*f;
com2+=g*f+b*e;
}
jcb=-com1-g*e-b*f;
jcb=-com2+b*e-g*f;
if(flag1!=1)
{
jcb=com2+b*e-g*f;
jcb=-com1+g*e+b*f;
}
else
{
jcb=-2*e;
jcb=-2*f;
}
}
}
}
//printf("\n第%d次输出-J矩阵如下：\n", t);
for(i=0;i<2*(node-1);i++){
for(j=0;j<2*(node-1);j++){
jcb*=-1.0;
// printf("%lf ", jcb);
}
// printf("\n\n");
}
/*算线形方程组*/
if (agaus(&jcb,&jcbb,2*(node-1))!=0) {
for (i=0,j=0;i<2*(node-1);i++) {
if(j==balc)
j++;
if(i%2==0){
e+=jcbb;
// printf("de[%d]=%lf e[%d]=%lf\n", (j+1),jcbb,(j+1),e);
}
else{
f+=jcbb;
// printf("df[%d]=%lf f[%d]=%lf\n",(j+1),jcbb,(j+1),f);
j++;
}
}
}
}
// printf("共经过了%d次迭代\n", t);
/*求平衡节点P，Q*/
i=balc;
com1=0.0, com2=0.0;
for(j=0;j<node;j++){
com1+=g*e-b*f;
com2+=g*f+b*e;
}
p=e*com1+f*com2;
q=f*com1-e*com2;
/*求其他各节点P，Q*/
for(i=0;i<node;i++){
if(i==balc) continue;
else
{
com1=0.0, com2=0.0;
for(j=0;j<node;j++)
{
com1+=g*e-b*f;
com2+=g*f+b*e;
}
p=e*com1+f*com2;
q=f*com1-e*com2;
}
}
// printf("最终结果为：\n");
// printf("节点号 电压 相角 有功 无功 e f\n");
for(i=0;i<node;i++)
{
v=sqrt(e*e+f*f);
angle=atan(f/e)*180/PI;
sump+=p;
sumq+=q;
//printf("%d %lf%lf%lf%lf%lf%lf\n", (i+1), v, angle, p, q, e, f);
}
页:
[1]
