Matlab直角坐标牛拉法机算潮流程序代码

解决问题

以何仰赞电力系统分析（下）例11-5为例，编写直角坐标系下牛顿拉夫逊法潮流计算的Matlab程序。对于类似模型的计算机潮流计算问题具有较好的适用性。

简介

参照何仰赞电力系统分析（下）的牛拉法程序流程框图，将程序分为节点导纳矩阵的形成、判敛、不平衡量计算、雅可比矩阵计算、节点电压修正、计算功率等几个部分。

程序代码

clear                                               %清除工作空间变量
clc                                                 %清屏
Accuracy=1.0e-10;                                    %判敛精度
Renew=10;                                           %迭代次数上限
Node=[-0.3-0.4i 0 2;                               %行：节点
      0.2 0 2;                              %列：电源注入-负荷吸收、给定电压、节点类型
      -0.15-0.4i 0 2;
      0 1 1];
Line=[1 2 0 0.2i 0 1.05 1.05 1;                       %行：线路
      1 4 0 0.01 0 1 1 1;         %列：1-a侧节点   2-b侧节点  3-PI型导纳  4-PI型阻抗  5-PI型导纳  
      2 3 0 0.04 0 1 1 1;        %列：6-变比  7-阻抗在K侧或1侧  8-变压器1侧节点标号  
      3 4 0 0.04 0 1 1 1]; 
Cap=[0 0 0 0];                                
N_Node=size(Node,1);
N_Line=size(Line,1);
Num=0;
Renew_OK=0;

Y=zeros(N_Node);                                    %创建节点导纳矩阵
for i=1:N_Line
    Line(i,4)=1/Line(i,4);                          %矩阵中阻抗变导纳
    if Line(i,6)~=1                                 %变压器判断并PI等值变换
        K=Line(i,6);                                %变比
        if Line(i,1)==Line(i,8)                     %变压器1侧接母线a
            Line(i,3)=Line(i,4)*(K-1)/K;            %两种展开PI型等效电路方式为两支腿不同
            Line(i,5)=Line(i,4)*(1-K)/K^2;
        elseif Line(i,2)==Line(i,8)                 %变压器1侧接母线b
            Line(i,3)=Line(i,4)*(1-K)/K^2;
            Line(i,5)=Line(i,4)*(K-1)/K;
        else 
            disp('错误！请核对矩阵元素！');          %核对矩阵变压