Matlab牛拉法潮流计算程序Word文档下载推荐.docx
- 文档编号:469169
- 上传时间:2023-04-29
- 格式:DOCX
- 页数:7
- 大小:21.95KB
Matlab牛拉法潮流计算程序Word文档下载推荐.docx
《Matlab牛拉法潮流计算程序Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《Matlab牛拉法潮流计算程序Word文档下载推荐.docx(7页珍藏版)》请在冰点文库上搜索。
请输入由各支路参数形成的矩阵:
B1='
B2=input('
请输入各节点参数形成的矩阵:
B2='
Y=zeros(n);
e=zeros(1,n);
f=zeros(1,n);
V=zeros(1,n);
sida=zeros(1,n);
S1=zeros(nl);
%%%---------------------------------------------------
fori=1:
nl %支路数
ifB1(i,6)==0 %左节点处于1侧
p=B1(i,1);
q=B1(i,2);
else%左节点处于K侧
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;
%对角元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;
%迭代次数ICT1、a;
不满足收敛要求的节点数IT2
whileIT2~=0%N0=2*n雅可比矩阵的阶数;
N=N0+1扩展列
IT2=0;
a=a+1;
n
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)
%求i节点有功和无功功率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;
%X3=dQ/deJ(p,N)=DQ节点无功功率差
J(m,q)=X1;
J(m,N)=DP;
q=q+1;
%X1=dP/deJ(m,N)=DP节点有功功率差
J(p,q)=X4;
J(m,q)=X2;
%X4=dQ/dfX2=dp/df
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
p=2*i-1;
J(p,q)=X3;
%扩展列△Q
m=p+1;
J(m,q)=X1;
%扩展列△P
else
%===============下面是针对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;
%PV节点电压误差
J(p,q)=X6;
%PV节点有功误差
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或△U
fork2=k1:
N1 %从k+1列的Jacobi元素到扩展列的△P、△Q或△U
J(k,k2)=J(k,k2)./J(k,k);
%用K行K列对角元素去除K行K列后的非对角元素进行规格化
J(k,k)=1;
%对角元规格化K行K列对角元素赋1
%====================回代运算=======================================
ifk~=3%不是第三行k>
3
k4=k-1;
fork3=3:
k4 %用k3行从第三行开始到当前行的前一行k4行消去
fork2=k1:
N1 %k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
%消去运算(当前行k列元素消为0)
end%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;
%当前行第k列元素已消为0
ifk==N0%若已到最后一行
break;
%==================前代运算==================================
fork3=k1:
N0%从k+1行到2*n最后一行
fork2=k1:
N1%从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
%消去运算
J(k3,k)=0;
else%是第三行k=3
%======================第三行k=3的前代运算========================
N0%从第四行到2n行(最后一行)
N1%从第四列到2n+1列(即扩展列)
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
%消去运算(当前行3列元素消为0)
end%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
%当前行第3列元素已消为0
%====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
fork=3:
2:
N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
%------修改节点电压-----------
N0
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);
%显示各节点的实际电压标幺值E用复数表示
-----------------------------------------------------'
各节点的电压大小V为(节点号从小到大排列):
disp(V);
%显示各节点的电压大小V的模值
各节点的电压相角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=电压X注入电流的共轭值
各节点的功率S为(节点号从小到大排列):
disp(S);
%显示各节点的注入功率
各条支路的首端功率Si为(顺序同您输入B1时一致):
nl
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);
ZF=['
S('
num2str(p),'
'
num2str(q),'
)='
num2str(SSi(p,q))];
disp(ZF);
各条支路的末端功率Sj为(顺序同您输入B1时一致):
p=B1(i,1);
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))...
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
num2str(SSj(q,p))];
各条支路的功率损耗DS为(顺序同您输入B1时一致):
p=B1(i,1);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
DS('
num2str(DDS(i))];
figure
(1);
subplot(1,2,1);
plot(V);
xlabel('
节点号'
ylabel('
电压标幺值'
gridon;
subplot(1,2,2);
plot(sida);
电压角度'
figure
(2);
subplot(2,2,1);
bar(P);
节点注入有功'
subplot(2,2,2);
bar(Q);
节点注入无功'
subplot(2,2,3);
P1=real(Siz);
Q1=imag(Siz);
bar(P1);
支路号'
支路首端注入有功'
subplot(2,2,4);
bar(Q1);
支路首端注入无功'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 牛拉法 潮流 计算 程序