loadflow.m

潮流计算

function [bus_sol,line_sol,line_flow] = ...
   loadflow(bus,line,tol,iter_max,acc,display,flag)
% Syntax:    [bus_sol,line_sol,line_flow] =
% loadflow(bus,line,tol,iter_max,acc,display,flag)
% 8/12/97
% Purpose:   solve the load-flow equations of power systems
%            modified to eliminate do loops and improve the use 
%            sparse matices解电力系统的潮流方程，用稀疏矩阵的特点减小循环次数，改进算法
% Input:    bus       - bus data 节点数据
% 输入          line      - line data 线路数据
%           tol       - tolerance for convergence 收敛精度
%           iter_max  - maximum number of iterations 最大迭代次数
%           acc       - acceleration factor 加速因素
%           display   - 'y', generate load-flow study report
%                        else, no load-flow study report
%                        显示：0输入y，显示发电机潮流报告，否则不显示
%           flag      - 1, form new Jacobian every iteration
%                       2, form new Jacobian every other 
%                           iteration
%                       标志位：1，每次迭代形成新的雅克比矩阵，2，隔一次形成一个新的雅克比矩阵
% Output:   bus_sol   - bus solution (see report for the
% 输出                        solution format)节点的潮流解
%           line_sol  - modified line matrix线路修正矩阵
%           line_flow - line flow solution (see report)线路潮流解
%
% See also:  
%
% Algorithm: Newton-Raphson method using the polar form of 
%   the equations for P(real power) and Q(reactive power).
%算法：有功和无功极坐标形式的方程的牛顿－拉夫逊法
% Calls:     Y_sparse, calc, form_jac chq_lim
%调用
%
% (c) Copyright 1991 Joe H. Chow - All Rights Reserved
%
% History (in reverse chronological order)
% Modification to correct generator var error on output 
% Graham Rogers November 1997
% Version:   2.1
% Author:     Graham Rogers
% Date:      October 1996
% Purpose:   To add generator var limits and on-load tap changers
% Version:   2.0
% Author:    Graham Rogers
% Date:      March 1994
% Version:   1.0
% Authors:   Kwok W. Cheung, Joe H. Chow
% Date:      March 1991
%
% ***********************************************************
global bus_int                                         %定义全局变量
global Qg bus_type g_bno PQV_no PQ_no ang_red volt_red 
global Q Ql
global gen_chg_idx
global ac_line n_dcl
tt = clock;     % start the total time clock开始计时
jay = sqrt(-1); %-1的平方根
load_bus = 3;   %赋值给load_bus，3是负荷节点
gen_bus = 2;    %2是发电机节点
swing_bus = 1;  %1是平衡节点
if exist('flag') == 0 % exist：检查函数或变量是否被定义，为0是没有定义
   flag = 1;  %那么flag置1
end
lf_flag = 1; %置为1


% set solution defaults设置缺省值：
if isempty(tol);tol = 1e-9;end  %如果收敛精度为空，则设为1e-9
if isempty(iter_max);iter_max = 30;end  %如果最大迭代次数为空，设为30
if isempty(acc);acc = 1.0; end; %如果加速因素为空，设为1.0
if isempty(display);display = 'n';end; %如果显示为空，设为n

if flag <1 | flag > 2 %如果flag<1或者>2
   error('LOADFLOW: flag not recognized')%显示错误
end
[nline nlc] = size(line);    % number of lines and no of line colsb把线路矩阵的行数和列数赋给nline和nlc
[nbus ncol] = size(bus);     % number of buses and number of col把节点矩阵的行数和列数赋给nbus和ncol


% set defaults设置缺省值：

% bus data defaults节点数据缺省值：

%如果节点矩阵列数小于15，当矩阵不完整时，设置默认值
if ncol<15 
   % set generator var limits设置发电机无功的限制
   if ncol<12 %如果列数小于12，把矩阵bus第11列的元素赋值为一个nbus*1的全1阵元素*9999
                              %把12列的元素赋值为nbus*1的全1阵元素*-9999
                              
      bus(:,11) = 9999*ones(nbus,1);%即：把11、12列的元素都赋为9999和-9999（11列是最大无功，12列是最小无功，无功的上下限）
      bus(:,12) = -9999*ones(nbus,1);
   end 
  if ncol<13;bus(:,13) = ones(nbus,1);end
   %如果节点矩阵列数小于13，bus第13列元素为1（13列是额定电压）
   bus(:,14) = 1.5*ones(nbus,1);%bus第14列元素为1.5（最大电压标幺值）
   bus(:,15) = 0.5*ones(nbus,1);%bus第15列元素为0.5（最小电压标幺值）
   volt_min = bus(:,15);%电压的最小限制是bus矩阵的15列元素
   volt_max = bus(:,14);%电压的最大限制是bus矩阵的14列元素
else
   volt_min = bus(:,15);%电压最小限制是第15列元素
   volt_max = bus(:,14);%电压最大限制是第14列元素
end

%如果bus矩阵是完整的，但有些节点的数据没有设置，设置这些数据
no_vmin_idx = find(volt_min==0);% 查找最小电压限制为零的母线号码对应的列数
if ~isempty(no_vmin_idx)%如果找到这个位置
   volt_min(no_vmin_idx) = 0.5*ones(length(no_vmin_idx),1);%这个位置的元素赋为0.5
end
no_vmax_idx = find(volt_max==0);%最大电压的默认值是1.5
if ~isempty(no_vmax_idx)
   volt_max(no_vmax_idx) = 1.5*ones(length(no_vmax_idx),1);
end
no_mxv = find(bus(:,11)==0);%找11列元素为0的位置
no_mnv = find(bus(:,12)==0);%找12列元素为0的位置
if ~isempty(no_mxv);  bus(no_mxv,11)=9999*ones(length(no_mxv),1);end %最大无功的默认值是9999,列矩阵
if ~isempty(no_mnv);bus(no_mnv,12) = -9999*ones(length(no_mnv),1);end%最小无功的默认值-9999，列矩阵

no_vrate = find(bus(:,13)==0);%查找额定电压为零的列
if ~isempty(no_vrate);bus(no_vrate,13) = ones(length(no_vrate),1);end%额定电压的默认值1

tap_it = 0;
tap_it_max = 10;%分接头最大动作次数
no_taps = 0;
% line data defaults, sets all tap ranges to zero - this fixes taps设置线路默认数据
if nlc < 10%如果线路矩阵的列数小于10
   line(:,7:10) = zeros(nline,4);%将line的7到10列的所以元素赋为0
   no_taps = 1;  % disable tap changing不考虑分接头的影响,如果程序没有给出分接头调整，则无需调整
 
end

% outer loop for on-load tap changers有载分接头的外循环

%mm_chk=1;
%while (tap_it<tap_it_max&mm_chk)%迭代次数没达到最大
%   tap_it = tap_it+1;%继续迭代
   % build admittance matrix Y建立稀疏导纳矩阵
   [Y,nSW,nPV,nPQ,SB] = y_sparse(bus,line);
   % process bus data
   bus_no = bus(:,1);%节点序号，一个数组
   V = bus(:,2);
   ang = bus(:,3)*pi/180;%角度转化为弧度
   Pg = bus(:,4);
   Qg = bus(:,5);%全局变量
   Pl = bus(:,6);
   Ql = bus(:,7);
   Gb = bus(:,8);
   Bb = bus(:,9);
   bus_type = round(bus(:,10));% round：取最接近的整数
   qg_max = bus(:,11);
   qg_min = bus(:,12);
   sw_bno=ones(nbus,1);%全1阵  一列的矩阵
   g_bno=sw_bno;
   % set up index for Jacobian calculation
   %% form PQV_no and PQ_no
   bus_zeros=zeros(nbus,1);%空阵  一列
   swing_index=find(bus_type==1); %找到平衡节点的序号
   sw_bno(swing_index)=bus_zeros(swing_index);%这个矩阵是一列的矩阵，平衡节点处的元素是0，其他元素是1
 
   PQV_no=find(bus_type >=2);%非平衡节点
   PQ_no=find(bus_type==3);%负荷节点
   gen_index=find(bus_type==2);%发电机节点
   g_bno(gen_index)=bus_zeros(gen_index); % g_bno is a vector having ones everywhere but the generator bus locations
                                          %一列的矩阵，发电机节点处元素为0，其他元素为1
     
   % construct sparse angle reduction matrix
   il = length(PQV_no);
   ii = [1:1:il]';%1列
   ang_red = sparse(ii,PQV_no,ones(il,1),il,nbus);%生成一个稀疏矩阵，在(ii,PQV_no)的位置是1，其他是0
   % 0 1 0 0 0 0 0 0 0 0 0 0 0
   % 0 0 1 0 0 0 0 0 0 0 0 0 0
   % 0 0 0 1 0 0 0 0 0 0 0 0 0
   % 0 0 0 0 1 0 0 0 0 0 0 0 0
   % 0 0 0 0 0 1 0 0 0 0 0 0 0
   % 0 0 0 0 0 0 1 0 0 0 0 0 0
   % 0 0 0 0 0 0 0 1 0 0 0 0 0
   % 0 0 0 0 0 0 0 0 1 0 0 0 0 
   % 0 0 0 0 0 0 0 0 0 1 0 0 0
   % 0 0 0 0 0 0 0 0 0 0 1 0 0
   % 0 0 0 0 0 0 0 0 0 0 0 1 0
   % 0 0 0 0 0 0 0 0 0 0 0 0 1
   
   % construct sparse voltage reduction matrix
   il = length(PQ_no);
   ii = [1:1:il]';
   volt_red = sparse(ii,PQ_no,ones(il,1),il,nbus);%稀疏矩阵，在(ii,PQ_no)的位置是1，其他是0
%volt_red =
%   (1,3)        1