牛拉法潮流计算

自动化07-1班段佳 07051101

function nl;

%------------------------------------------------------------------------

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

%======================牛顿——拉夫逊法==============================

%===========================潮流计算=================================

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

%-----------------------------------------------------------------------

% % %－－－－－－－－－－－－－－－使用说明部分－－－－－－－－－－－－－－－－－－－－－－－－－－－

display('% %本程序的功能是用牛顿——拉夫逊法进行潮流计算');

display('% %本程序要求用户按照一定的格式将电力系统的参数制成excel表格，系统运行时将从excel中加载这些参数，随后后即可进行潮流计算');

display('% %为了方便运算，用户再给系统节点进行编号时，请按照先PQ节点，再PV节点，最后平衡节点的顺序从小到大编号');

display('% %电力系统潮流计算excel格式——支路参数：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳;5、支路的变比K:1；6、支路首端处于K侧为1，1侧为0'); display('% %电力系统潮流计算excel格式——节点参数：1、节点号；2、电压大小；3、相位角；4、发电机有功；5、发电机无功；6、负载有功；7、负载无功；8、节点类型'); %===================================================================

%==============================数据准备==============================

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

% %---------------------电力系统数据加载部分-----------------------------------------------

clear

x=0;

Branch=0;%支路参数

Note=0;%节点参数

[filename, pathname] = uigetfile('*.xls', 'please choose the excel file with your powersystem parameters ');%从外部excel导入电力系统潮流计算相关参数

try

if filename ~= 0

x=xlsread([pathname,filename],'sheet1', 'A3:F3');

Branch=xlsread([pathname,filename],'sheet1', 'A5:G10');%读支路参数

Note=xlsread([pathname,filename],'sheet1', 'A15:H19');%读节点参数

end

catch

%进行出错处理

errmsg = lasterr;

errordlg(errmsg,'Save as Error');

rethrow(lasterror);

end

% %---------------------支路参数初始化部分-----------------------------------------------

SB=100;UB=220;n=1;m=1;pr=0.0001;

SB=x(5);%功率基准值

UB=x(6);%电压基准值

n=x(1);%节点数

nl=x(2);%支路数

m=x(3);%PQ节点的个数

pr=x(4);;%误差精度

B1(:,1)=Branch(:,1);%1、支路首端号

B1(:,2)=Branch(:,2);%2、末端号

B1(:,3)=Branch(:,3)+Branch(:,4)*i;%3、支路阻抗

B1(:,4)=Branch(:,5)*i;%4、支路对地电纳

B1(:,5)=Branch(:,6);%5、支路的变比K:1；

B1(:,6)=Branch(:,7);%6、支路首端处于K侧为1，1侧为0'

% %

% %---------------------节点参数初始化部分--------------------------------------------------

U=ones(n,1);

a=zeros(n,1);

Ps=zeros(n,1);

Qs=zeros(n,1);

P=zeros(n,1);

Q=zeros(n,1);

detp=zeros(n-1,1);

detq=zeros(m,1);

deta=zeros(n-1,1);

detu=zeros(m,1);

k=0;%迭代次数

U=Note(:,2);%各节点电压初始值(标幺值)

a=Note(:,3);%各节点电压相位初始值（弧度）

Gp=Note(:,4);%各节点发电机有功功率初始值(标幺值)

Gq=Note(:,5);%各节点发电机无功功率初始值(标幺值)

Lp=Note(:,6);%各节点负载有功功率初始值(标幺值)

Lq=Note(:,7);%各节点负载无功功率初始值(标幺值)

type=Note(:,8);%节点类型，PQ节点=1 ，PV节点=2 ，平衡节点=3

for h=1:n

Ps(h)=Gp(h)-Lp(h);%各节点注入的有功功率

Qs(h)=Gq(h)-Lq(h);%各节点注入的无功功率

end

% % %---------------------导纳矩阵计算部分-----------------------------------------------------

Y=zeros(n);

for h=1:nl %支路数

if B1(h,6)==0 %左节点处于低压侧（6、支路首端处于K侧为1，1侧为0）