该楼层疑似违规已被系统折叠隐藏此楼查看此楼
n=input('Please input n\n'); %n表示系统的节点数
d0=input('Please input d0\n'); %d0表示系统的误差
m=input('Please inpit m\n'); %m表示系统的pq节点的个数
p0=ones(n-1,1);
q0=ones(m,1);
d=1;
dpv=ones(n-1,1);
dqv=ones(m,1);
dp=zeros(n-1,1);
dq=zeros(m,1);
for a=1:n
for b=1:n
G(a,b)=real(Y(a,b));
B(a,b)=imag(Y(a,b));
end
end %把Y矩阵的实部和虚部分别导出
for a=1:n-1
for b=1:n-1
B1(a,b)=B(a,b);
end
end %表示出B1
for a=1:m
for b=1:m
B2(a,b)=B(a,b);
end
end %表示出B2
while(d>d0)
for a=1:n-1
p0(a)=0;
for b=1:n
p0(a)=p0(a)+v(a)*v(b)*(G(a,b)*cos(w(a)-w(b))+B(a,b)*sin(w(a)-w(b)));
end
end %利用初值求得相应的功率
for a=1:m
q0(a)=0;
for b=1:n
q0(a)=q0(a)+v(a)*v(b)*(G(a,b)*sin(w(a)-w(b))-B(a,b)*cos(w(a)-w(b)));
end
end %利用初值求得相应的无功
for a=1:n-1
dp(a)=p(a)-p0(a);
end %求dp
for a=1:m
dq(a)=q(a)-q0(a);
end %求dq
for a=1:n-1
dpv(a)=dp(a)/v(a);
end %利用初值求出dp/v
for a=1:m
dqv(a)=dq(a)/v(a);
end %利用初值求出dq/v
dw=-1*B1\dpv;
for a=1:n-1
dw(a)=dw(a)/v(a);
end %求出dw的修正值
dv=-1*B2\dqv; %求出v的修正值
d=max(max(abs(dw)),max(abs(dv)));%误差衡量
for a=1:n-1
w(a)=w(a)+dw(a);
end %求修正后的ｗ
for a=1:m
v(a)=v(a)+dv(a);
end %求修正后的v
end
这个是基于PQ分解法解潮流方程，我是按照自己的习惯编写的==。在运行本程序前，需要自行建立p(有功)、q(无功)、Y(节点导纳矩阵)、w(相角)、v(电压幅值)矩阵。
其中，节点编号请按照PQ节点、PV节点、平衡节点的次序编号。为了获得更好的收敛，可以将w统一设置为0，v统一设置为1。PQ解耦法适用于大电力系统，对于小型配电网等电力系统可能会存在较大的计算误差。
