潮流计算结果有错误的原因:这个错误和matlab的版本有关,有些版本可以运行无错误,有些版本则不行。
6次左右收敛,PQ分解法是12次左右收敛。首先要确定是你读取的数据正确,这个可以找一个正确的潮流程序进行验证,也可以将真值作为初始值在你的潮流程序中进行验证。你提到的迭代次数给人感觉好像是PQ分解法。
解的多值性和存在性:
对于非线性方程组的求解,从数学的观点来看,应该有多组解。根据程序中所设定的初值,一般都能收敛到合理解。但也有收敛到不合理解(电压过低或过高)的特殊情况。这些解是数学解(因为它们满足节点平衡方程式)而不是实际解。
为此需改变运行条件后再重新计算。此外,对于潮流计算问题所要求的节点电压的分量(幅值和角度或实部和虚部)。只有当其为实数时才有意义。如果所给的运行条件中无实数解,则认为该问题无解。
因此,当迭代不收敛时,可能有两种情况:一是解(指实数解)不存在,此时需修改运行方式;另一是计算方法不收敛,此时需更换计算方法。
t=0;
s=0;
r=0;
w=0;
number=input('How many node are there=');
% Convert Pq to a new array
for ii=1:number
if data(ii,4)==1
t=t+1;
for jj=1:14
new_data1(t,jj)=data(ii,jj);
end;
a(1,t)=ii;
s=s+1; %record the number of the PQ node
end;
end;%Convert pv to a new array
for ii=1:number
if data(ii,4)==2
t=t+1;
for jj=1:14
new_data1(t,jj)=data(ii,jj);
end;
a(1,t)=ii;
r=r+1; %record the number of the PV node
end;
end;
%Convert set_v to a new array
for ii=1:number
if data(ii,4)==3
t=t+1;
for jj=1:14
new_data1(t,jj)=data(ii,jj);
end;
a(1,t)=ii;
w=w+1;
end;
end;%creat a new_data2
[x,y]=size(data2)
for ii=1:x
for jj=1:2
for mm=1:number
if data2(ii,jj)==a(1,mm)
new_data2(ii,jj)=mm;
end;
end;
end;
end;for ii=1:x
for jj=3:14
new_data2(ii,jj)=data2(ii,jj);
end;
end;%creat a Y
Y=zeros(number,number);
YY=zeros(number,number);
yy=zeros(number,number);for ii=1:x
% for jj=1:14
iii=new_data2(ii,1);
jjj=new_data2(ii,2);
if new_data2(ii,5)==2
sub=new_data2(ii,6)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))-new_data2(ii,7)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))i;
Y(iii,jjj)=-sub/new_data2(ii,14);
YY(iii,jjj)=sub/new_data2(ii,14);
Y(jjj,iii)=-sub/new_data2(ii,14);
YY(jjj,iii)=sub/new_data2(ii,14);
yy(iii,jjj)=(1-new_data2(ii,14))/(new_data2(ii,14)new_data2(ii,14))sub;
yy(jjj,iii)=(new_data2(ii,14)-1)/(new_data2(ii,14))sub;
else
Y(iii,jjj)=-new_data2(ii,6)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))+new_data2(ii,7)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))i;
YY(iii,jjj)=new_data2(ii,6)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))-new_data2(ii,7)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))i;
Y(jjj,iii)=-new_data2(ii,6)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))+new_data2(ii,7)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))i;
YY(jjj,iii)=new_data2(ii,6)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))-new_data2(ii,7)/(new_data2(ii,7)new_data2(ii,7)+new_data2(ii,6)new_data2(ii,6))i;
yy(iii,jjj)=new_data2(ii,8)/2i;
yy(jjj,iii)=new_data2(ii,8)/2i;
end;
%end;
end;for iii=1:number
Y(iii,iii)=0;
end;%for ii=1:x
% for jj=1:14
for iii=1:number
for jj=1:number
% if iii~=jj
Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);
% end;
end;
end;%creat B, Gfor ii=1:number
for jj=1:number
G(ii,jj)= real(Y(ii,jj));
B(ii,jj)= imag(Y(ii,jj));
end;
end;%creat Initial_P Initial_Q Initial_V
for ii=1:(s+r)
set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7))/100;
end;
for ii=1:s;
set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8))/100;
end;for ii=1:r
set_V(ii,1)=new_data1(ii+s,12)new_data1(ii+s,12);%try to modify for sike of correcting
end;Initial_p_q_v=[set_P;set_Q;set_V];
disp(Initial_p_q_v);
%creat Initial_e,Initial_f
for ii=1:number-1
e(ii,1)=1;
f(ii,1)=00;%change f to test used to be 10
end; e(number,1)=new_data1(number,12);
f(number,1)=0;
% e(64,1)=088;%test 118ieee
% f(64,1)=039395826829394;
% f(14,1)=0;
% e(10,1)=1045;
%e(11,1)=101;
%e(12,1)=107;
%e(13,1)=109;
%//////////////////////////////////////////////////////////////////////////
%/////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
% Start NEWTOWN CALULATION
for try_time=1:25
%Creat every node consume P Q and U
n=s;
m=r;
for ii=1:(n+m)
sum1=0;
for jj=1:(n+m+1)
sum1=sum1+e(ii,1)(G(ii,jj)e(jj,1)-B(ii,jj)f(jj,1))+f(ii,1)(G(ii,jj)f(jj,1)+B(ii,jj)e(jj,1));
end;
p(ii,1)=sum1;
end;for ii=1:n
sum2=0;
for jj=1:(n+m+1)
sum2=sum2+f(ii,1)(G(ii,jj)e(jj,1)-B(ii,jj)f(jj,1))-e(ii,1)(G(ii,jj)f(jj,1)+B(ii,jj)e(jj,1));
end;
q(ii,1)=sum2;
end;disp('q=');
disp(q);u=zeros((n+m),1);
for ii=(n+1):(n+m)
u(ii,1)=e(ii,1)e(ii,1)+f(ii,1)f(ii,1);
end;for ii=n+1:(n+m)
extra_u((ii-n),1)=u(ii,1);
end;
disp('extra_u=');
disp(extra_u);sum=[p;q;extra_u];
disp(sum)
disp(s);
disp(p);
%creat Jacobian
disp(n);
disp(m);
for ii=1:(n+m)
for jj=1:(n+m)
if (ii~=jj)
PF(ii,jj)=B(ii,jj)e(ii,1)-G(ii,jj)f(ii,1);
PE(ii,jj)=-G(ii,jj)e(ii,1)-B(ii,jj)f(ii,1);
else
ss=0;
qq=0;
for num=1:(n+m+1)
ss=ss+G(ii,num)f(num,1)+B(ii,num)e(num,1);
qq=qq+G(ii,num)e(num,1)-B(ii,num)f(num,1);
end;
PF(ii,jj)=-ss+B(ii,jj)e(ii,1)-G(ii,jj)f(ii,1);%TEST+1
PE(ii,jj)=-qq-G(ii,jj)e(ii,1)-B(ii,jj)f(ii,1);%TEST+1
end;
end;
end;copy=314159;disp('================copy================')
for ii=1:n
for jj=1:m+n
if (ii~=jj)
QE(ii,jj)=B(ii,jj)e(ii,1)-G(ii,jj)f(ii,1);%TEST+1
QF(ii,jj)=G(ii,jj)e(ii,1)+B(ii,jj)f(ii,1);%TEST+1
else
ss=0;
qq=0;
for num=1:(n+m+1)
ss=ss+G(ii,num)f(num,1)+B(ii,num)e(num,1);
qq=qq+G(ii,num)e(num,1)-B(ii,num)f(num,1);
end;
QF(ii,jj)=-qq+G(ii,jj)e(ii,1)+B(ii,jj)f(ii,1);%TEST+1
QE(ii,jj)=ss+B(ii,jj)e(ii,1)-G(ii,jj)f(ii,1);%TEST+1
end;
end;
end;
%disp('QF');
%disp(QF);
%disp('QE');
%disp(QE);
UE=zeros((n+m),(n+m));
UF=zeros((n+m),(n+m));for ii=n+1:n+m
for jj=1:(n+m)
if (ii~=jj)
UE(ii,jj)=0;
UF(ii,jj)=0;
else
ss=0;
qq=0;
for num=1:(n+m+1)
ss=ss+G(ii,num)f(num,1)+B(ii,num)e(num,1);
qq=qq+G(ii,num)e(num,1)-B(ii,num)f(num,1);
end;
UF(ii,jj)=-2f(ii,1);
UE(ii,jj)=-2e(ii,1);
end;
end;
end;for ii=(n+1):(n+m)
for jj=1:(n+m)
extra_UE((ii-n),jj)=UE(ii,jj);
extra_UF((ii-n),jj)=UF(ii,jj);
end;
end;
%disp('extra_UE');
%disp(extra_UE);
%disp('extra_Uf');
%disp(extra_UF);
Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];
%disp('Jacobian=');
%disp(Jacobian);%creat substract result
substract_result=Initial_p_q_v-sum;
%disp('substract_result');
%disp(substract_result);
%calculate delta_f_e
delta_f_e=-inv(Jacobian)substract_result;
%disp(delta_f_e);
for ii=1:number-1;
f(ii,1)=f(ii,1)+delta_f_e(ii,1);
e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1);
end;if max(substract_result)<1e-4
break;
end ;
end;%disp('substract_result');
%disp(substract_result);
%disp('e=');%disp(e);
%disp('f=');
%disp(f);for ii=1:number
uuu(ii,1)= e(ii,1)e(ii,1)+f(ii,1)f(ii,1);
U_RESULT(ii,1)=sqrt(uuu(ii,1));
end;for ii=1:number
for jj=1:number
if ii==a(1,jj)
Old_Uresult(ii,1)=U_RESULT(jj,1)
end;
end;
end;for ii=1:number
Old_Uresult(ii,2)=ii;
end;
%disp('U_result');%disp(U_RESULT);
disp('=====================================');
disp('The last result is :')disp('===========U===================BUS-NO');
disp('U=')
disp(Old_Uresult);
%calculate the angle
PI=3141592
for ii=1:number
Angle(ii,1)=atan(f(ii,1)/e(ii,1))/PI180;
end;for ii=1:number
for jj=1:number
if ii==a(1,jj)
Old_Angle(ii,1)=Angle(jj,1);
Old_Angle(ii,2)=ii;
end;
end;
end;
disp('=======Angle===================BUS-NO');
disp('Angle=');
disp(Old_Angle);
disp('=====Try-times=======================')
disp('Try-times=')
disp(try_time);%disp('p====================');
%disp(p); % for jj=1:number
% if a(1,jj)==118
% man=jj
% end;%end;
%disp('man=========');
%disp(man)sum4=0;
for jj=1:number
Y_conj(number,jj)=conj(Y(number,jj));
sum4=sum4+Y_conj(number,jj)(e(jj,1)-f(jj,1)i);
end;
%sum4=sum4e(number,1);
disp('===============Balance P Q=========BUS-NO');
%disp(sum4);Blance_Q(1,1)=imag(sum4)100;
Blance_Q(1,2)=a(1,number);
Blance_P(1,1)=real(sum4)100;
Blance_P(1,2)=a(1,number);
disp('Q Of the Balance node= ');disp(Blance_Q);
disp('P Of the Balance node= ')
disp(Blance_P);
disp('=================================BUS-NO');
%calculate the Q of the P-V nodeQ_PV_node=zeros(number,2);
Y_conj=conj(Y);
for ii=(s+1):(s+r)
for jj=1:number
Q_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)i)(Y_conj(ii,jj)(e(jj,1)-f(jj,1)i));
end;
end;
for ii=(s+1):(s+r);
Q_PV_node(ii,1)=Q_PV_node(ii,1)100+new_data1(ii,8)i;
end;disp('This program is from blogsinacomcn/breadwinner') ; for ii=1:number
old_number=a(1,ii);
Q_PV_node_old(old_number,1)=Q_PV_node(ii,1);
end;
for ii=1:number
Q_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));
end;
for ii=1:number
Q_PV_node_old(ii,2)=ii;
end;disp('Q gen=');disp(Q_PV_node_old); 1、矩阵的稀疏性没有体现,delta_f_e=-inv(Jacobian)substract_result;最好用一个函数来求稀疏线性方程组的解,在该函数中可以考虑Jacobian阵的结构特点(用all函数),然后进行合适的节点编号,使得注入元达到最小,而后进行前代和回代过程;
2、可以考虑用极坐标来形成Jacobian矩阵,这样的话,矩阵在结构上是位子对称的,很容易用矩阵分块技术来形成结构对称的分块对角矩阵,而后用递归或者是并行求解该Jacobian方程,这样的方法比较适用于大规模的电力系统仿真计算中。
PQ分解法是一种用于计算电力系统潮流的方法,它的计算速度较快且占用的内存比较小,应用较为广泛。
41P-Q分解法的基本原理:P-Q分解法是从简化一极坐标表示的牛顿-拉夫逊法潮流修正方程基础上派生出来的,考虑到了电力系统本身的特点。
潮流计算的一般提法是:已知电力网络的结构和参数,已知各负荷点、电源点吸取或发出的有功功率和无功功率(PQ节点),给定电压控制点的电压幅值和有功功率(PV节点),对指定的一个平衡节点给定其电压幅值和相位角(Vθ点),求解全网各节点电压幅值和相位角,并进一步算出各支路的功率分布和网络损耗。求解潮流问题的基本方程式是节点功率平衡方程。若全网有n个节点,对其中任一节点,可写出其节点功率平衡方程式i=1,2,…,n
式中Pi、Qi分别为节点注入有功功率和无功功率,妭i为节点电压相量,Yik为节点导纳矩阵元素。这一方程描述了节点电压同功率之间的非线性关系,是潮流计算的基本方程式。对潮流计算的数字计算机求解方法提出的基本要求是:①计算速度快;②占用存储量少;③收敛性好;④方法简单。
数值解法潮流计算在数学上是求解一组非线性方程,基本的方法是迭代法。首先发展的潮流问题数字解法是导纳矩阵迭代法。它占用计算机存储量少,适合于计算机发展初期阶段的实际条件,其缺点是收敛性较差。其后发展了阻抗矩阵迭代法,克服了导纳矩阵迭代法收敛性差的缺点,但对大电力系统的计算,占用计算机存储量大。
60年代末期出现了以导纳矩阵为基础、采用稀疏矩阵和节点编号优化技术的牛顿-拉夫森法。该法以其在收敛性、存储量和计算时间方面的优越性逐渐取代了其他方法,在当今的潮流计算中应用得最为广泛。在数学上,牛顿-拉夫森法是求解非线性方程的有效方法。它把非线性方程的求解变成反复对相应的线性方程迭代求解的过程:设非线性方程组为F(X)=0,求解X的第t步迭代格式是 F'(X(t))ΔX(t)=-F(X(t))
X(t+1)=X(t)+ΔX(t)
式中F┡(X)是非线性函数矢量F(X)对变量矢量X的一阶偏导数矩阵,称为雅可比矩阵;ΔX为X的偏差矢量;上标t表示第t次迭代值。当X的相邻两次迭代值之差ΔX小于给定误差时,迭代收敛。一般潮流问题,经6~7次迭代即可收敛。牛顿-拉夫森法的迭代收敛性与初值X(10)的选取有关。当X(10)与其解X接近时,极易收敛。对电力系统潮流问题,当选用标幺值计算时,节点电压通常在1附近,给定各点电压初值为1,可以有相当好的收敛性。 在牛顿-拉夫森的基础上,发展了一种更加简化的方法,即PQ分解法。这种方法利用高压电网电抗值远大于电阻值、有功功率变化主要与电压相角有关、无功功率变化主要与电压幅值有关的特点,将有功功率和无功功率的迭代分开进行,使占用计算机存储量和计算量进一步减少。 潮流计算方法的进一步发展将使潮流计算更加快速,收敛性更好,适应各种系统运行状况的能力更强,并在一些专门领域中如优化潮流和在线潮流中得到实际的应用。 参考书目 东北电业管理局调度局编:《电力系统运行 *** 作和计算》,水利电力出版社,北京,1977。
这是牛顿法原理
把非线性函数f(x)在x = 0处展开成泰勒级数
牛顿法
取其线性部分,作为非线性方程f(x)=0的近似方程,则有
f(0 )+(x-0 ) f′(0 )=0
设f′(0 )≠0,则其解为x = - xf(1)
再把f(x)在x 处展开为泰勒级数,取其线性部分为f(x)=0的近似方程,若f′(x ) ≠0,则得x = - 如此继续下去,得到牛顿法的迭代公式:x = - (n=0,1,2,…) (2)
例1 用牛顿法求方程f(x)=x +4x -10=0在[1,2]内一个实根,取初始近似值x =15。 解 f′(x)=3x +8x所以迭代公式为:
x = - n=0,1, 2,
列表计算如下:
n
0
1
2
3
15
13733333
136526201
136523001
这是不一定的,要看情况,只是因为现在电力系统都比较复杂,才总体上表现为高斯-赛德尔法迭代次数比较多。
高斯-赛德尔法与PQ分解法、牛拉法所用的迭代矩阵不一样,收敛的快慢就是要看迭代矩阵的谱半径。谱半径小于1说明收敛,否则不收敛。谱半径越小,收敛速度越快。
对于Ax=b,PQ分解法和牛拉法的迭代方程为:
x=B1x+f1, B1=I-inv(D)A,f1=inv(D)b,D为对角矩阵。
对于Ax=b,高斯-赛德尔法的迭代方程为:
x=B2x+f2, B2=inv(D-L)U,f2=inv(D-L)b,L、U为下三角和上三角矩阵。
对于病态电网,例如重负荷、很多长线路、负电抗变压器等等,都是病态电网的表现。病态电网的特点使得迭代方程中的B1、B2的谱半径不一样,B2往往大于1(矩阵近似奇异或者奇异),而B1往往小于1,此时高斯-赛德尔法表现为发散,所以迭代很可能不收敛。
当B2的谱半径小于B1时,高斯-赛德尔法的迭代次数是要少于PQ分解法、牛拉法的。
是的。
PQ分解法无法应用于配电网络的潮流计算。
PQ分解法潮流计算:基于MATLAB软件的P-Q分解法潮流计算。电力系统潮流计算是研究电力系统稳态运行情况的一种重要的分析计算。
电力系统中的潮流计算归结为求解一个用功率作为注入的非线性方程组。牛顿拉夫逊法和PQ分解法都是求解潮流方程的典型算法。这方面的软件有很多,国外、国内都有。在国内应用较广的如电科院的PSASP,中国版BPA都有潮流计算的功能。不同的软件在输入方式、数据文件格式和输入内容上都有差别。但直接与功率方程相关的几个重要参数都是一样的。比如,不同类型节点(PQ,PV)的已知量,平衡节点定义,线路阻抗、变压器阻抗和变比等。
以上就是关于潮流计算结果有错误的原因全部的内容,包括:潮流计算结果有错误的原因、谁知道风电场并网的潮流计算MATLAB程序、写出牛顿拉夫逊法和PQ分解法的原理和优缺点及应用范围等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)