#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#define M 14 //预计支路数
#define N 14 //预计节点梁拍数
//
struct linestr //支路信息
{
int i,j //和该支路相连的两个节点
float r,x,g,b //该支路的参数，kb为非标准变比
float pij,qij,pji,qji //支路潮流—计算数据
}line[M]
struct pointstr //节点信息
{
int sign //节点类型，0—平衡节点，1—PQ节点，2—PV节点
float voltage,cta,p,q //节点电压幅值和相位角，节点注入有功、无功
}pointpri[N],pointaft[N] //pointpri用于存放节点的初始数据，pointaft用于存放计算后的节点数据
float G[N][N]={0},B[N][N]={0} //节点导纳阵
//
float B1[N][N]={0},B2[N][N]={0}//修正方程系数矩阵
float detp[N]={0},detct[N]={0},detq[N]={0},detv[N]={0}//依次放修正方程的修正量
float pn[N]={0},ctn[N]={0},qn[N]={0},vn[N]={0} //依次放修正量对应的节点号，即变量对应嫌渣的节点号
int nline,npoint //准确的支路数和节点数
int npq,npv,nbal //PQ节点数、PV节点数、平衡节点数
void creat()
{
int i
npq=0npv=0nbal=0
printf("\n请输入总支路数：")scanf("%d",&nline)
for(i=1i<=nlinei++)
{
printf("请输入%d号支路数据：\n",i)
printf("和该支路相连的节点号i,j（请用逗号隔开）：")scanf("%d,%d",&(line[i].i),&(line[i].j))
printf("该支路芹渣悄的电阻和电抗R,X（请用逗号隔开）：")scanf("%f,%f",&(line[i].r),&(line[i].x))
printf("该支路的对地电导和电纳G,B（请用逗号隔开）：")scanf("%f,%f",&(line[i].g),&(line[i].b))
line[i].pij=0line[i].qij=0line[i].pji=0line[i].qji=0
//形成节点导纳阵
G[line[i].i][line[i].j]=-line[i].r/((line[i].r*line[i].r+line[i].x*line[i].x))
G[line[i].j][line[i].i]=G[line[i].i][line[i].j]
B[line[i].i][line[i].j]=line[i].x/((line[i].r*line[i].r+line[i].x*line[i].x))
B[line[i].j][line[i].i]=B[line[i].i][line[i].j]
G[line[i].i][line[i].i]+=line[i].r/((line[i].r*line[i].r+line[i].x*line[i].x))
G[line[i].i][line[i].i]+=line[i].g
B[line[i].i][line[i].i]+=(-line[i].x)/((line[i].r*line[i].r+line[i].x*line[i].x))
B[line[i].i][line[i].i]+=line[i].b
G[line[i].j][line[i].j]+=line[i].r/(line[i].r*line[i].r+line[i].x*line[i].x)
G[line[i].j][line[i].j]+=line[i].g
B[line[i].j][line[i].j]+=(-line[i].x)/(line[i].r*line[i].r+line[i].x*line[i].x)
B[line[i].j][line[i].j]+=line[i].b
}
printf("\n请输入总节点数：")scanf("%d",&npoint)
for(i=1i<=npointi++)
{
printf("请输入节点%d的数据\n",i)
printf("选择节点类型（1－平衡节点，2－PQ节点，3－PV节点）：")
scanf("%d",&(pointpri[i].sign))
if(pointpri[i].sign==1)
{
printf("请输入该节点的电压幅值和相位角（用逗号隔开）：")
scanf("%f,%f",&(pointpri[i].voltage),&(pointpri[i].cta))
pointpri[i].p=0pointpri[i].q=0
nbal++
}
else if(pointpri[i].sign==2)
{
printf("请输入该节点的注入有功和无功（用逗号隔开）：")
scanf("%f,%f",&(pointpri[i].p),&(pointpri[i].q))
pointpri[i].voltage=1pointpri[i].cta=0
npq++
}
else if(pointpri[i].sign==3)
{
printf("请输入该节点的有功注入和电压幅值（用逗号隔开）：")
scanf("%f,%f",&(pointpri[i].p),&(pointpri[i].voltage))
pointpri[i].q=0pointpri[i].cta=0
npv++
}
else
{
printf("选择无效，请重新输入\n")
i--
}
if(nbal>1)
{
printf("出现多个平衡节点，请重新输入\n")
i--nbal--
}
}
}
void view()
{
int i,j
for(i=1i<=npointi++)
{
for(j=1j<=npointj++)
printf("%f+j%f ,",G[i][j],B[i][j])
printf("\n")
}
for(i=1i<=npointi++)
{
printf("节点%d参数：",i)
if(pointpri[i].sign==1)printf("平衡节点 v=%f,theta=%f\n",pointpri[i].voltage,pointpri[i].cta)
else if(pointpri[i].sign==2)printf("PQ节点 P=%f,Q=%f\n",pointpri[i].p,pointpri[i].q)
else printf("PV节点 P=%f,v=%f\n",pointpri[i].p,pointpri[i].voltage)
}
}
void modify()
{
}
void run()
{
int i,j,k,m,n
float E,e//E为题目所要求的迭代精度，由键盘输入，e为每一步迭代的实际计算精度，用以和E比较
int numPQ[N]//PQ节点号
int numPV[N]//PV节点号
int numBAL //平衡节点号
int np_th[N],nq_v[N]//迭代方程中的有功-电压相位角方程、无功－电压幅值方程分别对应的节点号
float B1_fac[N][N]={0},B2_fac[N][N]={0},a
j=1k=1m=1n=1
for(i=1i<=npointi++){B1_fac[i][i]=1B2_fac[i][i]=1}
for(i=1i<=npointi++)
{//将数据存入一个新的结构体数组，保证在迭代过程中即使数据有改变，仍能保留住原始数据
pointaft[i].sign=pointpri[i].sign
pointaft[i].voltage=pointpri[i].voltage
pointaft[i].cta=pointpri[i].cta
pointaft[i].p=pointpri[i].p
pointaft[i].q=pointpri[i].q
//统计PV节点和PQ节点序号，并存入临时数组
if(pointpri[i].sign==1)numBAL=i
else if(pointpri[i].sign==2){numPQ[j++]=inp_th[m++]=inq_v[n++]=i}
else {numPV[k++]=inp_th[m++]=i}
}
//形成有功－相位角的迭代方程系数矩阵B1
for(i=1i<=npoint-1i++)
for(j=1j<=npoint-1j++)
B1[i][j]=B[np_th[i]][np_th[j]]
//形成因子表
for(k=1k<=npoint-1k++)
for(i=k+1i<=npoint-1i++)
{ a=B1[i][k]
for(j=1j<=npoint-1j++)
{
B1[i][j]-=a*B1[k][j]/B1[k][k]
B1_fac[i][j]-=a*B1_fac[k][j]/B1[k][k]
}
}
//形成无功－电压幅值的迭代方程的系数矩阵B2
for(i=1i<=npqi++)
for(j=1j<=npqj++)
B2[i][j]=B[nq_v[i]][nq_v[j]]
//形成因子表
for(k=1k<=npqk++)
for(i=k+1i<=npqi++)
{ a=B2[i][k]
for(j=1j<=npqj++)
{
B2[i][j]-=a*B2[k][j]/B2[k][k]
B2_fac[i][j]-=a*B2_fac[k][j]/B2[k][k]
}
}
printf("请输入迭代精度：")scanf("%f",&E)
k=-1
do
{
e=0
k++
//有功迭代
//首先计算有功不平衡量
for(i=1i<=npoint-1i++)
{
detp[i]=pointaft[np_th[i]].p/pointaft[np_th[i]].voltage
for(j=1j<=npointj++)
detp[i]-=pointaft[j].voltage*(G[np_th[i]][j]*cos(pointaft[np_th[i]].cta-pointaft[j].cta)+B[np_th[i]][j]*sin(pointaft[np_th[i]].cta-pointaft[j].cta))
printf("%f,",detp[i])
}
for(i=1i<=npoint-1i++)
{
detct[i]=0
for(j=1j<=npoint-1j++)
detct[i]+=B1_fac[i][j]*detp[j]
}
for(i=npoint-1i>=1i--)
{
for(j=npoint-1j>ij--)
detct[i]-=B1[i][j]*detct[j]
detct[i]/=B1[i][i]
}
for(i=1i<=npoint-1i++)detct[i]/=(-pointaft[np_th[i]].voltage)//求出相位角的修正量
for(i=1i<=npoint-1i++){pointaft[np_th[i]].cta+=detct[i]e+=detct[i]*detct[i]}//修正角度,并计算精度
printf("第%d次迭代后的相位角：",k)
for(i=1i<=npointi++)printf("%f,",pointaft[i].cta)
printf("\n")
//无功迭代
//首先计算无功不平衡量
for(i=1i<=npqi++)
{
detq[i]=pointaft[nq_v[i]].q/pointaft[nq_v[i]].voltage
for(j=1j<=npointj++)
detq[i]-=pointaft[j].voltage*(G[nq_v[i]][j]*sin(pointaft[nq_v[i]].cta-pointaft[j].cta)-B[nq_v[i]][j]*cos(pointaft[nq_v[i]].cta-pointaft[j].cta))
printf("%f,",detq[i])
}
for(i=1i<=npqi++)
{
detv[i]=0
for(j=1j<=npqj++)
detv[i]+=B2_fac[i][j]*detq[j]
}
for(i=npqi>=1i--)
{
for(j=npqj>ij--)
detv[i]-=B2[i][j]*detv[j]
detv[i]/=B2[i][i]
}
for(i=1i<=npqi++)detv[i]/=(-1)//求出电压的修正量
for(i=1i<=npqi++){pointaft[nq_v[i]].voltage+=detv[i]e+=detv[i]*detv[i]}//修正电压
printf("第%d次迭代后的电压幅值为：",k)
for(i=1i<=npointi++)printf("%f,",pointaft[i].voltage)
printf("\n")
}while(e>E&&k<100)
//输出结果
printf("在满足精度E＝%f的情况下，计算结果为：",E)
for(i=1i<=npointi++)
printf("v%d=%f,cta%d=%f\n",i,pointaft[i].voltage,i,pointaft[i].cta)
//计算平衡节点注入功率
pointaft[numBAL].p=0pointaft[numBAL].q=0
for(i=1i<=npointi++)
{
pointaft[numBAL].p+=pointaft[i].voltage*(G[numBAL][i]*cos(pointaft[numBAL].cta-pointaft[i].cta)+B[numBAL][i]*sin(pointaft[numBAL].cta-pointaft[i].cta))
pointaft[numBAL].q+=pointaft[i].voltage*(G[numBAL][i]*sin(pointaft[numBAL].cta-pointaft[i].cta)-B[numBAL][i]*cos(pointaft[numBAL].cta-pointaft[i].cta))
}
pointaft[numBAL].p*=pointaft[numBAL].voltage
pointaft[numBAL].q*=pointaft[numBAL].voltage
printf("平衡节点注入功率为：%f+j%f\n",pointaft[numBAL].p,pointaft[numBAL].q)
//计算各支路潮流
for(k=1k<=nlinek++)
{
float aa,bb,cc
i=line[k].ij=line[k].j
aa=pointaft[i].voltage*cos(pointaft[i].cta)-pointaft[j].voltage*cos(pointaft[j].cta)
bb=pointaft[i].voltage*sin(pointaft[i].cta)-pointaft[j].voltage*sin(pointaft[j].cta)
cc=line[k].r*line[k].r+line[k].x*line[k].x
line[k].pij=(line[k].r*(pointaft[i].voltage*cos(pointaft[i].cta)*aa+pointaft[i].voltage*sin(pointaft[i].cta*bb))-line[k].x*(pointaft[i].voltage*sin(pointaft[i].cta*aa-pointaft[i].voltage*cos(pointaft[i].cta)*bb)))/cc
line[k].qij=(line[k].x*(pointaft[i].voltage*cos(pointaft[i].cta)*aa+pointaft[i].voltage*sin(pointaft[i].cta*bb))+line[k].r*(pointaft[i].voltage*sin(pointaft[i].cta*aa-pointaft[i].voltage*cos(pointaft[i].cta)*bb)))/cc
line[k].pji=(line[k].r*(pointaft[j].voltage*cos(pointaft[j].cta)*(-aa)+pointaft[j].voltage*sin(pointaft[j].cta*(-bb)))-line[k].x*(pointaft[j].voltage*sin(pointaft[j].cta*(-aa)-pointaft[j].voltage*cos(pointaft[j].cta)*(-bb))))/cc
line[k].qji=(line[k].x*(pointaft[j].voltage*cos(pointaft[j].cta)*(-aa)+pointaft[j].voltage*sin(pointaft[j].cta*(-bb)))+line[k].r*(pointaft[j].voltage*sin(pointaft[j].cta*(-aa)-pointaft[j].voltage*cos(pointaft[j].cta)*(-bb))))/cc
printf("由节点%d流向节点%d的功率为：%f+j%f\n",i,j,line[k].pij,line[k].qij)
printf("由节点%d流向节点%d的功率为：%f+j%f\n",j,i,line[k].pji,line[k].qji)
}
}
main ()
{
int a
printf("\n")
printf("***************电力系统潮流计算**************\n")
printf("\n")
printf(" 设计者： 李 冰 设计：2014 年 5 月\n")
printf("\n")
printf("\n")
printf(" 注：各参数用标么值表示,角度均以弧度制表示\n")
printf("\n")
printf("\n")
do
{
printf("请选择:1-建立电网2-查看数据3-修改数据4-运行计算0-退出程序:")//由键盘输入选项
scanf("%d",&a)
if (a==0) //如果选中0，就表示退出程序
{
exit(0)
}
else if (a==1) //如果选中1，就跳去执行建立电网，即creat处
{
creat()
}
else if (a==2) //如果选中2，就跳去执行查看数据，即view处
{
view()
}
else if (a==3) //如果选中1，就跳去执行修改数据，即modify处
{
modify()
}
else if (a==4) //如果选中1，就跳去执行运行计算，即run处
{
run()
}
else //如果以上都不是，则表示无效输入，重新要求输入选择
{
printf("\n无效输入，请重新选择\n")
}
}while(1)
}
不知道能不能帮到你
function lianxuchaoliuclear
clc
n=9%节点数；
nl=9%支弊铅游路数；
isb=1%平衡节点号；
pr=0.00001%误差精度；
b1=[1 4 0.0576i 0 1.05 14 5 0.017+0.092i 0.158i 1 05 6 0.039+0.17i 0.358i 1 03 6 0.0586i 0 1.05 16 7 0.0119+0.1008i 0.209i 1 07 8 0.0085+0.072i 0.149i 1 08 2 0.0625i 0 1.05 18 9 0.032+0.161i 0.306i 1 09 4 0.01+0.085i 0.176i 1 0]
%依次是支路首端；末端，支路阻抗；对地电纳；支路变比；折算到哪一侧标志（高压侧为1；低压侧为0）；
b2=[0 0 1.05 1.05 0 11.63 0 1.05 1.05 0 30.85 0 1.05 1.05 0 30 0 1 0 0 20 0.9+0.3i 1 0 0 20 0 1 0 0 20 1+0.35i 1 0 0 20 0 1 0 0 20 1.25+0.5i 1 0 0 2]
%依次是节点的发电机功率；负荷功率；节租销点电压初激改值；PV节点电压V给定值；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；
Y=zeros(n)%求导纳阵；
for i=1:nl
if b1(i,6)==0
p=b1(i,1)q=b1(i,2)
else
p=b1(i,2)q=b1(i,1)
end
Y(p,q)=Y(p,q)-1./(b1(i,3)*b1(i,5))
Y(q,p)=Y(p,q)
Y(q,q)=Y(q,q)+1./(b1(i,3)*b1(i,5)^2)+b1(i,4)./2
Y(p,p)=Y(p,p)+1./b1(i,3)+b1(i,4)./2
end
%disp('系统的导纳阵为：')
%disp(Y)
G=real(Y)B=imag(Y)
for i=1:n
e(i)=real(b2(i,3))
f(i)=imag(b2(i,3))
v(i)=b2(i,4)
end
for i=1:n
S(i)=b2(i,1)-b2(i,2)
B(i,i)=B(i,i)+b2(i,5)
end
P=real(S)Q=imag(S)
w=zeros(2*n,1)Jac=zeros(2*n)
t=0
while t==0
for i=1:n
if b2(i,6)~=isb
C=0D=0
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j)
D=D+G(i,j)*f(j)+B(i,j)*e(j)
end
if b2(i,6)==2%P，Q节点；
w(2*i)=P(i)-e(i)*C-f(i)*D
w(2*i-1)=Q(i)-f(i)*C+e(i)*D
else if b2(i,6)==3%P，V节点；
w(2*i)=P(i)-e(i)*C-f(i)*D
w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2)
end
end
else
w(2*i-1)=0
w(2*i)=0
end
end
%disp(w)
w1=w(3:2*n)
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P，Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P，V节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=0
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-2*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
Jac(2*i-1,2*j-1)=0
Jac(2*i,2*j)=0
Jac(2*i-1,2*j)=0
Jac(2*i,2*j-1)=0
end
end
end
%disp(Jac)
Jac1=Jac(3:2*n,3:2*n)
for k=1:2*n-2
m=0
for i=k+1:2*n-2
m=Jac1(i,k)./Jac1(k,k)
w1(i)=w1(i)-m*w1(k)
for j=k+1:2*n-2
Jac1(i,j)=Jac1(i,j)-m*Jac1(k,j)
end
end
end
E(2*n-2)=w1(2*n-2)./Jac1(2*n-2,2*n-2)
for i=2*n-3:-1:1
c=0
for k=i+1:2*n-2
c=c+Jac1(i,k)*E(k)
E(i)=(w1(i)-c)./Jac1(i,i)
end
end
%disp(E)
for i=1:n-1
e(i+1)=e(i+1)-E(2*i-1)
f(i+1)=f(i+1)-E(2*i)
end
%disp(e)
%disp(f)
for i=1:n-1
b(2*i-1)=abs(E(2*i-1))
b(2*i)=abs(E(2*i))
end
KB=max(b)
%disp(KB)
if KB<pr
t=1
else
t=0
end
end
%disp(e)
%disp(f)
for i=1:n
fz(i)=sqrt(e(i)^2+f(i)^2)
end
disp(fz)
JJ1=zeros(2*n-1,2*n-1)Jac1=zeros(2*n-2,2*n-2)ttt=0ccc=1%迭代次数；
L=0.05T=0%L是步长，T是参数；
while ttt==0
%确定扩展矩阵
for i=1:n-1
if b2(i+1,6)==3
JJ1(2*i,2*n-1)=1.5*1.05
JJ1(2*i-1,2*n-1)=0
else if b2(i+1,6)==2
if b2(i+1,2)==0
JJ1(2*i-1,2*n-1)=0
JJ1(2*i,2*n-1)=0
else
JJ1(2*i,2*n-1)=-real(b2(i+1,2))
JJ1(2*i-1,2*n-1)=-imag(b2(i+1,2))
end
end
end
end
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2% P,Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P,V节点，
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
Jac(2*i-1,2*j)=0
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
Jac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
Jac(2*i-1,2*j)=-2*f(i)
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
Jac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
Jac(2*i-1,2*j-1)=0
Jac(2*i,2*j)=0
Jac(2*i-1,2*j)=0
Jac(2*i,2*j-1)=0
end
end
end
Jac1=Jac(3:2*n,3:2*n)
for i=1:2*n-2
for j=1:2*n-2
JJ1(i,j)=Jac1(i,j)
end
end
if ccc==1%以负荷为连续参数；
for j=1:2*n-2
JJ1(2*n-1,j)=0
end
JJ1(2*n-1,2*n-1)=1
w2=zeros(2*n-1,1)
for i=1:2*n-2
w2(i,1)=0
end
w2(2*n-1,1)=1
end
%disp(JJ1)
if ccc>1%以切向量中分量最大（绝对值最大）的状态变量选定为连续参数；
for i=1:2*n-2
Jd(i)=abs(d(i))
end
for i=1:2*n-3
if Jd(i)>=Jd(i+1)
zd=Jd(i)ek=i
else if Jd(i)<=Jd(i+1)
zd=Jd(i+1)ek=i+1
end
end
end
for j=1:2*n-1
JJ1(2*n-1,j)=0
end
JJ1(2*n-1,ek)=1
w2=zeros(2*n-1,1)
for i=1:2*n-1
w2(i,1)=0
end
if d(ek)>0
w2(ek,1)=1
else if d(ek)<0
w2(ek,1)=-1
end
end
end
for k=1:2*n-1
m=0
for i=k+1:2*n-1
m=JJ1(i,k)./JJ1(k,k)
w2(i)=w2(i)-m*w2(k)
for j=k+1:2*n-1
JJ1(i,j)=JJ1(i,j)-m*JJ1(k,j)
end
end
end
d(2*n-1)=w2(2*n-1)./JJ1(2*n-1,2*n-1)
for i=2*n-2:-1:1
c=0
for k=i+1:2*n-1
c=c+JJ1(i,k)*d(k)
d(i)=(w2(i)-c)./JJ1(i,i)
end
end
%disp(d)
for i=1:n-1
e(i+1)=e(i+1)+L*d(2*i-1)
f(i+1)=f(i+1)+L*d(2*i)
end
T=T+L*d(2*n-1)
disp(d(2*n-1))
%disp(e)
%disp(f)
%disp(T)
%对预估的近似解进行校正；
tt=0
while tt==0
for i=1:n
if b2(i,6)~=isb
C=0D=0
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j)
D=D+G(i,j)*f(j)+B(i,j)*e(j)
end
if b2(i,6)==2% P,Q节点；
if b2(i,2)~=0
wj(2*i)=P(i)-T*real(b2(i,2))-e(i)*C-f(i)*D
wj(2*i-1)=Q(i)-T*imag(b2(i,2))-f(i)*C+e(i)*D
else
wj(2*i)=P(i)-e(i)*C-f(i)*D
wj(2*i-1)=Q(i)-f(i)*C+e(i)*D
end
else if b2(i,6)==3%P,V节点；
wj(2*i)=P(i)+T*1.5*1.05-e(i)*C-f(i)*D
wj(2*i-1)=0
end
end
else
wj(2*i)=0
wj(2*i-1)=0
end
end
wj1=wj(3:2*n)
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P,Q节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
JJac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i)
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
JJac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i)
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
JJac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i)
end
end
else if b2(i,6)==3%P,V节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i))
JJac(2*i-1,2*j)=0
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i)
JJac(2*i-1,2*j-1)=0
else if j==i
m=0h=0
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r)
h=h+G(i,r)*f(r)+B(i,r)*e(r)
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i)
JJac(2*i-1,2*j)=-2*f(i)
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i)
JJac(2*i-1,2*j-1)=-2*e(i)
end
end
end
end
else
JJac(2*i-1,2*j-1)=0
JJac(2*i,2*j)=0
JJac(2*i-1,2*j)=0
JJac(2*i,2*j-1)=0
end
end
end
JJac1=JJac(3:2*n,3:2*n)
for k=1:2*n-2
m=0
for i=k+1:2*n-2
m=JJac1(i,k)./JJac1(k,k)
wj1(i)=wj1(i)-m*wj1(k)
for j=k+1:2*n-2
JJac1(i,j)=JJac1(i,j)-m*JJac1(k,j)
end
end
end
E1(2*n-2)=wj1(2*n-2)./JJac1(2*n-2,2*n-2)
for i=2*n-3:-1:1
c=0
for k=i+1:2*n-2
c=c+JJac1(i,k)*E1(k)
E1(i)=(wj1(i)-c)./JJac1(i,i)
end
end
%disp(E1)
for i=1:n-1
e(i+1)=e(i+1)-E1(2*i-1)
f(i+1)=f(i+1)-E1(2*i)
end
for i=1:n-1
bx(2*i-1)=abs(E1(2*i-1))
bx(2*i)=abs(E1(2*i))
end
KB1=max(bx)
if KB1<pr
tt=1
else
tt=0
end
end
%disp(e)
%disp(f)
if d(2*n-1)>0
ttt=0ccc=ccc+1
else
ttt=1
end
end
disp(T)
disp(ccc)
%disp(e)
%disp(f)
for i=1:n
fz1(i)=sqrt(e(i)^2+f(i)^2)
end
disp(fz1)
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）