牛拉潮流程序

来源：年旅网

clear; clc;

%%%+++++++++++++++++++++++++++++++++++++输入带有变压器的支路矩阵中各节点对应各变比 % function[]

%========================================================================== [NODE,Branch]=OpDF_; %打开矩阵 (test.m)文件 Node=NODE;

N=Node(:,1); %节点号 Type=Node(:,2); %节点类型

BR=Branch %将支路信息保存在BR中%%

K=Branch(:,6); %支路变压器变比，0代表没有变压器 n=length(N); %节点数

nbr=length(K); %支路数

Total_of_Bus1=size(NODE);%%%取节点矩阵的行和列

Total_of_Bus=Total_of_Bus1(1,1)%%bus矩阵的行数即节点数 Total_of_Branch1=size(Branch);%%%%取支路矩阵的行和列

Total_of_Branch=Total_of_Branch1(1,1);%% 支路branch矩阵行数即支路数

Z=zeros(Total_of_Bus1); %将节点排序重新存储节点信息%%%%定义为节点数的方阵 format short

b=1; %排序标志位 pq=0; %PQ节点标志位 pv=0; %PV节点标志位 ph=0; %平衡节点标志位

%---------------------按照PQ,PV,平衡节点的次序排序各种节点

%-------------------------------------统计PQ节点数 0代表是pq节点 for a=1:Total_of_Bus if NODE(a,2) == 0 Z(b,:)=NODE(a,:); b=b+1;

pq=pq+1; end end

%-------------------------------------统计PV节点数2代表pv节点 for a=1:Total_of_Bus if NODE(a,2) == 2 Z(b,:)=NODE(a,:); b=b+1;

pv=pv+1; end end

%-------------------------------------统计平衡节点数 3代表平衡节点 for a=1:Total_of_Bus if NODE(a,2) == 3 Z(b,:)=NODE(a,:); b=b+1;

ph=ph+1; end end Z

Z2=Z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将节点进行重新排序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mm=zeros(n,1); for i=1:n

mm(i,1)=i; end

Z1(:,1)=mm(:,1);

Branch1=zeros(nbr,2); for i=1:n

if Z(i,1)~=Z1(i) for j=1:nbr

if Branch(j,1)==Z(i,1) Branch1(j,1)=Z1(i); end

if Branch(j,2)==Z(i,1) Branch1(j,2)=Z1(i); end end

else

for j=1:nbr

if Branch(j,1)==Z(i,1) Branch1(j,1)=Z(i,1); end

if Branch(j,2)==Z(i,1) Branch1(j,2)=Z(i,1); end end end end

Branch(:,1)=Branch1(:,1); Branch(:,2)=Branch1(:,2); Z(:,1)=Z1(:,1);

j=sqrt(-1);

%--------矩阵已经完成按照PQ，PV，平衡节点的顺序排列起来------ YSNODE=Z; %保存排序后的原始节点数据

%%%%=======================================================================

Y=zeros(n,n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求互导纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:n

for t=1:nbr

if (Branch(t,1)==i||Branch(t,2)==i)&& Branch(t,6)==0 %%非变压器支路%%%%%%% Y(Branch(t,1),Branch(t,2))=-1/(Branch(t,3)+j*Branch(t,4)); Y(Branch(t,2),Branch(t,1))=Y(Branch(t,1),Branch(t,2));

else if (Branch(t,1)==i||Branch(t,2)==i)&&Branch(t,6)~=0 %%变压器支路%%%%%% Y(Branch(t,1),Branch(t,2))=(-1/(j*Branch(t,4)))/Branch(t,6);

Y(Branch(t,2),Branch(t,1))=Y(Branch(t,1),Branch(t,2)); end end end end

%%%%%%%%%%%%%%%%%%%%%%%%%求自导纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:n

for t=1:nbr

if (Branch(t,1)==i||Branch(t,2)==i)&& Branch(t,6)==0 %%非变压器支路%%%%%%% Y(i,i)=Y(i,i)+1/(Branch(t,3)+j*Branch(t,4))+(1/2)*j*Branch(t,5);

else if Branch(t,1)==i&&Branch(t,6)~=0 %%%%%%变压器支路——————————且i为首节点%%%%%%%%

Y(i,i)=Y(i,i)+1/(j*Branch(t,4));

else if Branch(t,2)==i&&Branch(t,6)~=0%%%%%%%%%%%__变压器支路————且i为末节点%%%%%%%%

Y(i,i)=Y(i,i)+(1/(j*Branch(t,4)))/(Branch(t,6)*Branch(t,6)); end end end end end

%%%%%%%%%若有并联电容器组，则自导纳要加上并联电容器的导纳%%%%%%%%%%%%%%%%%%%%%% for i=1:n

if NODE(i,13)~=0

Y(i,i)=Y(i,i)+j*NODE(i,13) end end Y

n=length(N);

G=real(Y); %实部，即电导 B=imag(Y); %虚部，即电纳

%%%%%%%%%%%%%%%%%%%%%%%%%%%给定初始的电压值与相位值%%%%%%%%%%%%%%%%%%%%%% U_first=Z(:,3); %初始电压幅值 phase_first=Z(:,4); %初始相位值 e=U_first.*cos(phase_first); f=U_first.*sin(phase_first);

%%%%%%%%%%%%%%%%%计算Delta_P初始功率量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P=Z(:,5); %节点负荷有功分量 Q=Z(:,6); %节点负荷无功分量

PG=Z(:,7); %发电机发出的有功 QG=Z(:,8); %发电机发出的无功 U0=Z(:,9); %节点电压都的初始值 Delta_P=zeros(1,n-1); for i=1:n-1 for j=1:n

Delta_P(i)=Delta_P(i)-e(i)*(G(i,j)*e(j)-B(i,j)*f(j))+f(j)*(G(i,j)*f(j)+B(i,j)*e(j)); end end

for i=1:n-1

Delta_P(i)=Delta_P(i)-(P(i)-PG(i)); end Delta_P

%%%%%%%%%%%%%%%%%%计算Delta_Q初始功率量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% m=0;

for i=1:n;

if Type(i)==2; %计算PV节点的个数 m=m+1; %m代表pv节点个数 end end

Delta_Q=zeros(1,n-m-1); for i=1:n-m-1 for j=1:n

Delta_Q(i)=Delta_Q(i)-f(i)*(G(i,j)*e(j)-B(i,j)*f(j))+e(i)*(G(i,j)*f(j)+B(i,j)*e(j)); end end

for i=1:n-m-1

Delta_Q(i)=Delta_Q(i)-(Q(i)-QG(i)); end Delta_Q

Delta_V=zeros(1,m); for i=1:m for j=1:n

if Type(j)==2

Delta_V(i)=U0(j)^2-(e(j)^2+f(i)^2); end end end Delta_V num=0;

disp(['第',num2str(num),'次时的Delta总的失配量为：'])

%--------------------------进入循环体判断是否满足条件--------------------------------% %--------------------先算出最大值，作为判断是否收敛的依据----------------------------% DEL=[Delta_P Delta_Q]; %Delta_P Delta_Q______________%%%%%% MAX =max(abs(DEL)); MAX

Theta_first=zeros(1,n);

U_f=U_first';

Delta_F_E1=[Theta_first(1:n-1) U_f(1:n-m-1)]; Delta_F=Delta_F_E1';

Delta_Cor=Delta_F_E1; %_______________Delta_the Delta_u________________%%

disp(['第一次最大失配量误差：',num2str(MAX)]) %-------------循环判断-----------%

if MAX>1e-004 % 判断依据

disp('-----------------------下面开始下一次迭代过程！-----------------------') end

while MAX>1e-004 num=num+1;

%%%%%%%%%%%%%%%%%%%形成雅克比矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%% %--------------------------先求非对角元素--(H)-------------% Hik=zeros(n-1,n-1); for i=1:n-1 for k=1:n-1 if i~=k