牛顿——拉夫逊潮流计算的matlab程序

这是我本科毕业设计编写的程序，用matlab语言编写的牛顿-拉夫逊潮流计算程序，程序注释详细。程序从excel中读取原始数据，并将计算结果输出到excel表格中，压缩包内的excel是一个实例的数据，程序可以直接使用。具体的excel和matlab混合编程的资料大家可以自己在网上搜索。

%高斯——赛尔德潮流计算
%[s_nodeu_node]=gaosai（YmnsS_PQP_PVU_PVusQmaxQminexpl）
%返回值：s_node:  节点功率
%       u_node：  节点电压
%输入值：Y：   导纳矩阵
%        m：   PQ节点数包含平衡节点
%        n：   总结点数
%        s：   平衡节点编号，一般定为1号编号
%        S_PQ: PQ节点的视在功率[1:m]其中平衡节点的视在功率暂定为0；
%        P_PV：PV节点的有功功率[m+1:n]
%        U_PV: PV节点的电压给定值只有幅值，没有相位角，是实数[m+1:n]
%        us:   平衡节点的电压给定值，复数。
%        Qmax；PV节点无功功率最大值
%        Qmin：PV节点无功功率最小值
%        epxl: 最大电压变量dumax的最大值，即规定的潮流计算的精度，实数。
function [u_node]=DYgaosai（YmnsS_PQP_PVU_PVusQmaxQminepxl）
%**************************************************************************
%以下为高斯塞尔德潮流计算程序
 fprintf（‘高斯——塞尔德潮流计算开始\n‘）
U=ones（1n）;%节点电压的初始值；
U（s）=us;%将电压初始值的平衡节点的初始电压设为给定值
k=1000;%设置迭代次数
dumax=0;%最大电压变化量，实数
flag=zeros（1n-m）;%PV节点在计算过程中转化为PQ节点的标志，PV转化为PQ节点时，相应的标志位置1
Q_PV=ones（1n-m）;%PV节点的无功功率；
Q_PV=Q_PV.*0.2;
Ss=0;%平衡节点的视在功率
for a=m+1:n
    U（a）=U_PV（a-m）;%第a个节点是PV节点，暂时将PV节点电压给定值赋给节点电压初始值其相角暂时为0
end
while（k>=1）
for i=1:n
    u=U（i）;%保存U（i）的初值，为求du和dumax做准备；
    if（i==s）%第i个节点是否为平衡节点？
       % i=i+1;%d第i个节点是平衡节点，跳过。
    else    %第i节点不是平衡节点，进行以下处理
        A=0;
        for a=1:n
            if（a~=i）
             A=A+Y（ia）*U（a）;
            end
        end
        a=1;
        if（i>=m+1 && flag（i-m）==0） %第i各节点是否为PV节点？
            Q_PV（i-m）=-imag（U（i）‘*Y（ii）*U（i）+U（i）‘*A）;
            if（Q_PV（i-m）<=Qmax && Q_PV（i-m）>=Qmin）
                Upv=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）-A）/Y（ii）;
                ang=angle（Upv）;
                U（i）=abs（U（i））*（cos（ang）+sin（ang）*j）;
%               Q_PV（i-m）=-imag（U（i）‘*Y（ii）*U（i）+U（i）‘*A）;
            else
                if（Q_PV（i-m）>Qmax）
                    Q_PV（i-m）=Qmax;
                    flag（i-m）=1;
                    fprintf（‘PV节点%gQ大于Qmax转化为PQ节点\n‘i）;
                elseif（Q_PV（i-m）                    Q_PV（i-m）=Qmin;
                    flag（i-m）=1;
                    fprintf（‘PV节点%gQ小于Qmin转化为PQ节点\n‘i）;
                end
                U（i）=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）‘-A）/Y（ii）;
            end
        else
            if（i<=m）%第i个节点是PQ节点
                U（i）=（S_PQ（i）‘/U（i）‘-A）/Y（ii）;
            elseif（i>=m+1 && flag（i-m）==1）
                U（i）=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）‘-A）/Y（ii）;
            end
        end
    end
 du=abs（U（i）-u）;
  if（du>dumax）
  dumax=du;
  end
end
i=1;
if（dumax    fprintf（‘迭代次数%g\n‘1000-k+1）
    break;
else
    k=k-1;
    dumax=0;
end
end
if（k==0）
    fprintf（‘高斯——塞尔德潮流计算不能收敛！！！\n‘）
end
for i=1:n %计算平衡节点的视在功率
        Ss=Ss+Y（si）‘*U（i）‘;
end
Ss=Ss*U（s）;%平衡节点的视在功率
    for i=1:n%计算各个节点的视在功率
        if（i==s）
            s_node（i）=Ss;
        else
            if（i<=m）
                s_node（i）=S_PQ（i）;
            elseif（i>=m+1）
                s_node（i）=P_PV（i-m）+Q_PV（i-m）*j;
            end
        end
    end
    s_node
    u_node=U
    fprintf（‘高斯——塞尔德潮流计算结束\n‘）

    
            
            

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3531  2011-05-07 20:29  newton\DYgaosai.m

     文件      21808  2011-06-06 09:25  newton\newton.m

     文件       2597  2011-05-12 08:59  newton\write.m

     文件      16384  2011-06-06 09:25  newton\输出数据.xls

     目录          0  2011-10-18 08:58  newton

----------- ---------  ---------- -----  ----

                44320                    5