MATLAB潮流计算仿真Word文档下载推荐.docx
《MATLAB潮流计算仿真Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《MATLAB潮流计算仿真Word文档下载推荐.docx(11页珍藏版)》请在冰豆网上搜索。
163.159e-2+7.479e-2i0.1215i10;
233.68e-3+0.11135i00.9091;
453.68e-3+0.11135i00.9091;
462.808e-2+6.648e-2i0.108i10;
673.0865e-3+0.0833i00.9091;
683.159e-2+7.479e-2i0.1215i10;
6102.457e-2+5.817e-2i0.0945i10;
893.0865e-3+0.0833i00.9091;
8102.808e-2+6.648e-2i0.108i10];
请输入由支路参数形成的矩阵:
B1='
B2=[001.051.0501;
001002;
00.6+0.3718i1002;
00.4+0.247i1002;
00.35+0.2169i1002;
00.5+0.3099i1002;
0.801.051.0503];
请输入各节点参数形成的矩阵:
B2='
Y=zeros(n);
e=zeros(1,n);
f=zeros(1,n);
V=zeros(1,n);
sida=zeros(1,n);
S1=zeros(nl);
%-------修改部分------------
ym=0;
SB=100;
UB=220;
%ym=input('
您输入的参数是标么值?
(若不是则输入一个不为零的数值)'
ifym~=0
%SB=input('
请输入功率基准值:
SB='
%UB=input('
请输入电压基准值:
UB='
YB=SB./UB./UB;
BB1=B1;
BB2=B2;
fori=1:
nl
B1(i,3)=B1(i,3)*YB;
B1(i,4)=B1(i,4)./YB;
end
disp('
B1矩阵B1='
disp(B1)
n
B2(i,1)=B2(i,1)./SB;
B2(i,2)=B2(i,2)./SB;
B2(i,3)=B2(i,3)./UB;
B2(i,4)=B2(i,4)./UB;
B2(i,5)=B2(i,5)./SB;
B2矩阵B2='
disp(B2)
end
%%%---------------------------------------------------
fori=1:
nl%支路数
ifB1(i,6)==0%左节点处于低压侧
p=B1(i,1);
q=B1(i,2);
else
p=B1(i,2);
q=B1(i,1);
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;
%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
%求导纳矩阵
disp('
导纳矩阵Y='
disp(Y)
%----------------------------------------------------------
G=real(Y);
B=imag(Y);
%分解出导纳阵的实部和虚部
n%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
n%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);
%i节点无功补偿量
%===================================================================
P=real(S);
Q=imag(S);
ICT1=0;
IT2=1;
N0=2*n;
N=N0+1;
a=0;
whileIT2~=0
IT2=0;
a=a+1;
ifi~=isb%非平衡节点
C(i)=0;
D(i)=0;
forj1=1:
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);
%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);
%Σ(Gij*fj+Bij*ej)
P1=C(i)*e(i)+f(i)*D(i);
%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);
%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)
%求P'
Q'
V2=e(i)^2+f(i)^2;
%电压模平方
%=========以下针对非PV节点来求取功率差及Jacobi矩阵元素=========
ifB2(i,6)~=3%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%===============以上为除平衡节点外其它节点的功率计算=================
%=================求取Jacobi矩阵===================
ifj1~=isb&
j1~=i%非平衡节点&
非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
%dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
%dP/df=dQ/de
X3=X2;
%X2=dp/dfX3=dQ/de
X4=-X1;
%X1=dP/deX4=dQ/df
p=2*i-1;
q=2*j1-1;
J(p,q)=X3;
J(p,N)=DQ;
m=p+1;
J(m,q)=X1;
J(m,N)=DP;
q=q+1;
J(p,q)=X4;
J(m,q)=X2;
elseifj1==i&
j1~=isb%非平衡节点&
对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);
%dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);
%dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);
%dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);
%dQ/df
J(p,q)=X3;
%扩展列△Q
m=p+1;
J(p,q)=X4;
%扩展列△P
J(m,q)=X2;
%===============下面是针对PV节点来求取Jacobi矩阵的元素===========
%PV节点有功误差
DV=V(i)^2-V2;
%PV节点电压误差
%dP/de
%dP/df
X5=0;
X6=0;
J(p,q)=X5;
J(p,N)=DV;
J(p,q)=X6;
X5=-2*e(i);
X6=-2*f(i);
%=========以上为求雅可比矩阵的各个元素=====================
fork=3:
N0%N0=2*n(从第三行开始,第一、二行是平衡节点)
k1=k+1;
N1=N;
%N=N0+1即N=2*n+1扩展列△P、△Q
fork2=k1:
N1%扩展列△P、△Q
J(k,k2)=J(k,k2)./J(k,k);
%非对角元规格化
J(k,k)=1;
%对角元规格化
ifk~=3%不是第三行
%============================================================
k4=k-1;
fork3=3:
k4%用k3行从第三行开始到当前行前的k4行消去
N1%k3行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
%消去运算
J(k3,k)=0;
ifk==N0
break;
%==========================================
fork3=k1:
N0
N1
%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====
2:
N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
%------修改节点电压-----------
DET=abs(J(k,N));
ifDET>
=pr%电压偏差量是否满足要求
IT2=IT2+1;
%不满足要求的节点数加1
ICT2(a)=IT2;
ICT1=ICT1+1;
%用高斯消去法解"
w=-J*V"
迭代次数:
'
disp(ICT1);
没有达到精度要求的个数:
disp(ICT2);
fork=1:
V(k)=sqrt(e(k)^2+f(k)^2);
sida(k)=atan(f(k)./e(k))*180./pi;
E(k)=e(k)+f(k)*j;
%===============计算各输出量===========================
各节点的实际电压标幺值E为(节点号从小到大排列):
disp(E);
EE=E*UB;
disp(EE);
-----------------------------------------------------'
各节点的电压大小V为(节点号从小到大排列):
disp(V);
VV=V*UB;
disp(VV);
各节点的电压相角sida为(节点号从小到大排列):
disp(sida);
forp=1:
C(p)=0;
forq=1:
C(p)=C(p)+conj(Y(p,q))*conj(E(q));
S(p)=E(p)*C(p);
各节点的功率S为(节点号从小到大排列):
disp(S);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
SS=S*SB;
disp(SS);
各条支路的首端功率Si为(顺序同您输入B1时一致):
ifB1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Si(p,q));
SSi(p,q)=Si(p,q)*SB;
ZF=['
S('
num2str(p),'
'
num2str(q),'
)='
num2str(SSi(p,q))];
disp(ZF);
%disp(SSi(p,q));
各条支路的末端功率Sj为(顺序同您输入B1时一致):
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Sj(q,p));
SSj(q,p)=Sj(q,p)*SB;
num2str(SSj(q,p))];
%disp(SSj(q,p));
各条支路的功率损耗DS为(顺序同您输入B1时一致):
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i)*SB;
DS('
num2str(DDS(i))];
%disp(DDS(i));
figure
(1);
subplot(2,2,1);
plot(V);
xlabel('
节点号'
ylabel('
电压标幺值'
gridon;
subplot(2,2,2);
plot(sida);
电压角度'
subplot(2,2,3);
bar(real(S));
节点注入有功'
subplot(2,2,4);
bar(Siz);
支路首端无功'
(注:
文档可能无法思考全面,请浏览后下载,供参考。
可复制、编制,期待你的好评与关注)