利用matlab基础语句实现牛顿法matlab潮流计算
2019-04-10 本文已影响0人
季天泽
问题描述
现代电力系统分析课程留了用计算机进行潮流计算的作业，可以使用各种软件实现，matlab就包含在其中，按老师的意思应该是使用matlab中的某个app。
可是我头铁，app是什么，当然要自己写啊！
题目来自于 电力系统分析教材 华中科技大学版
整体思路
根据Y阵，PQ节点的P、Q值，PV节点的P、V值，平衡节点的δ、V值，利用极坐标的牛顿法，求出每个节点的δ、V值和P、Q值，和ΔP、ΔQ值，根据ΔP、ΔQ值判断是否收敛，同时得出下一个P、Q值，反复迭代求解。
原理及步骤
- 输入Y阵Ymat，PQ、PV节点参数Ps、Qs、Ps、V，平衡节点参数，此为目标值。
- 给出初始V、δ值，其中V节点、平衡节点为1中的值，在迭代过程中不再变化。
- 根据δ、U值求其对应的P、Q值。
- 根据P、Q、δ、V、Ymat求Jacobi矩阵J。
- 将复数矩阵Ymat进行分离，实部组成G阵，虚部组成B阵。根据如下公式计算J阵。
- 根据Ps与P，Qs与Q的差求出修正量DP、DQ。
- 根据J阵与DP、DQ，利用牛顿法公式
其中:
求Δδ、ΔV。
- 根据Δδ、ΔV、δ、V求得下一个δ、V。
- 重复3~8。
程序代码
/chaoliujtz.m
%季天泽
%华中科技大学电气与电子工程学院
%强电磁与新技术国家重点实验室
clear
clc
% 节点必须按照PQ、PV、平衡节点来排序
% PQ节点有m个，共有n+1个节点，即PV节点有n-m个
% 每个PQ节点有1个关于Q的方程式，每个PV或PQ节点有1个关于P的方程式，平衡节点没有方程式
% 所以P有n行，Q有m行，U与δ有n+1行
Ymat=[
1.042093-i*8.242876 -0.588235+i*2.352941 i*3.666667 -0.453858+i*1.891074
-0.588235+i*2.352941 1.069005-i*4.727377 0 -0.480769+i*2.403846
i*3.666667 0 -i*3.333333 0
-0.453858+i*1.891074 -0.480769+i*2.403846 0 0.934627-i*4.261590
];
z=max(size(Ymat));
%Ps、Qs代表实际的P、Q
Ps(1)=-0.3;
Ps(2)=-0.55;
Ps(3)=0.5;
Qs(1)=-0.18;
Qs(2)=-0.13;
n=size(Ps);
m=size(Qs);
n=n(2);
m=m(2);
%delta=[ones(1,10)*-1,0]
delta(z)=0;
%输入U时，注意PV节点和平衡节点直接输入的是目标值的标幺值，而PQ节点输入的仅是迭代初值
U=[1,1,1.1,1.05];
output=[];
error=[];
G=real(Ymat);
B=imag(Ymat);
nn=10%将迭代次数设为10
for ii=1:nn
p1=(U*(G.*cos(delta-delta')+B.*sin(delta-delta')));
q1=(U*(G.*sin(delta-delta')-B.*cos(delta-delta')));
P=U(1:n).*p1(1:n);
Q=U(1:m).*q1(1:m);
clear p1 q1;
DP=Ps-P;
DQ=Qs-Q;
J = Jacobi(P,Q,Ymat,U,delta);
invJ=inv(J);
%修正方程
M=-invJ*[DP';DQ'];
Ddelta=M(1:n)';
DU=(diag(U(1:m))*M(n+1:m+n))';
Ddelta(n+1)=0;
DU(n+1)=0;
delta=delta+Ddelta;
U=U+DU;
error=[error;[DP,DQ,DU]];
output=[output;U,delta];
end
plot(1:nn,output)
/Jacobi.m
%季天泽
%华中科技大学电气与电子工程学院
%强电磁与新技术国家重点实验室
%
function J = Jacobi(P,Q,Y,U,delta)
% U,Y,S都是复矩阵，获取之后进行拆分,m为PQ节点的个数，n为节点总数，n-1-m为PV节点数
% 节点按照PQ、PV、平衡节点来排序
% PQ节点有m个，共有n+1个节点，即PV节点有n-m个
% 每个PQ节点有1个关于Q的方程式，每个PV或PQ节点有1个关于P的方程式，平衡节点没有方程式
% 所以P有n行，Q有m行
n=length(P);
m=length(Q);
G=real(Y);
B=imag(Y);
Q(n+1)=0;
P(n+1)=0;
for i=1:n
for j=1:n
if i~=j
H(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
elseif i==j
H(i,j)=U(i)*U(j)*B(i,j)+Q(i);
end
end
end
for i=1:n
for j=1:m
if i~=j
N(i,j)=-U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
elseif i==j
N(i,j)=-U(i)*U(j)*G(i,j)-P(i);
end
end
end
for i=1:m
for j=1:n
if i~=j
K(i,j)=U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
elseif i==j
K(i,j)=U(i)*U(j)*G(i,j)-P(i);
end
end
end
for i=1:m
for j=1:m
if i~=j
L(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
elseif i==j
L(i,j)=U(i)*U(j)*B(i,j)-Q(i);
end
end
end
J=[
H, N;
K, L
];
程序运行结果
数据迭代结果
数据迭代结果折线图
结论
可以看出，用牛顿法求解电力系统潮流方程，仅需3步即可收敛，收敛速度极快。
