loadflow.m

潮流计算

%   (2,4)        1
%   (3,5)        1
%   (4,8)        1
 %  (5,9)        1
%   (6,10)       1
 %  (7,12)       1
  % (8,13)       1
   
   iter = 0;     % initialize iteration counter迭代的初始值
   
   % **************calculate the power mismatch and check convergence第一次计算功率偏差，检查收敛，利用calc.m
      [delP,delQ,P,Q,conv_flag] =...
      calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol);
  %计算得到是否收敛，并得到conv_flag,delP,delQ,P,Q
   
     st = clock;     % start the iteration time clock迭代开始计时
   %% start iteration process for main Newton_Raphson solution用牛拉法进行迭代
   while (conv_flag == 1 & iter < iter_max)%conv_flag=1表示没有收敛，没有超过最大迭代次数
      iter = iter + 1;
     
      %*******************Form the Jacobean matrix形成雅可比矩阵 见form_jac.m
      clear Jac
      Jac=form_jac(V,ang,Y,ang_red,volt_red);
      
     % reduced real and reactive power mismatch vectors简化的功率偏差矢量
      red_delP = ang_red*delP;% 列矩阵，length=(PQV_no)除去了平衡节点，（12行）
      red_delQ = volt_red*delQ;% 列矩阵， 除去了平衡节点和发电机节点，（8行）
      clear delP delQ
      % solve for voltage magnitude and phase angle increments电压幅值和相角增量
      temp = Jac\[red_delP; red_delQ];%Jac*temp=[red_delp;red_delQ]——修正方程；[]中是一列矢量；解修正方程；temp是列矢量（20）
      % expand solution vectors to all buses将解矢量扩展为所有节点
      delAng = ang_red'*temp(1:length(PQV_no),:);%temp(1:length(PQV_no),:)表示temp的前12个元素，即相角的增量。delAng是在平衡节点的位置上置0
      delV = volt_red'*temp(length(PQV_no)+1:length(PQV_no)+length(PQ_no),:);%temp的后8个元素，即电压增量，
                                                                             %恢复13个元素的矩阵，在平衡节点和发电机节点处为0
      
      % update voltage magnitude and phase angle修正电压幅值和相角
      V = V + acc*delV;     %修正各点电压
      V = max(V,volt_min);  % voltage higher than minimum比最小电压高的电压
      V = min(V,volt_max);  % voltage lower than maximum比最高电压低的电压，即没有超出范围
      ang = ang + acc*delAng;%修正各点相角
     
      %************************ calculate the power mismatch and check convergence重新计算功率偏差，检查收敛情况
      [delP,delQ,P,Q,conv_flag] =...调用calc.m
         calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol);
     
     % check if Qg is outside limits检查无功是否在范围内
      gen_index=find(bus_type==2);%找出发电机节点
      Qg(gen_index) = Q(gen_index) + Ql(gen_index);%发电机节点的无功＝计算出的无功＋负荷的无功
 %进行发电机结点类型转换

    %******************调用chq_lim可能改变节点类型
      lim_flag = chq_lim(qg_max,qg_min);%检查发电机无功有无超出范围，见chq_lim.m
        if lim_flag == 1; %lim_flag＝1表示节点有改变
           disp('Qg at var limit');%显示（）
        end
   end  %如果收敛（由calc得到），检查完无功越界，或达到最大迭代次数30，就停止迭代

   if iter == iter_max%如果迭代次数达到了最大次数
      imstr = int2str(iter_max);%转化为字符串
      disp(['inner ac load flow failed to converge after ', imstr,' iterations'])%显示交流潮流在第多少次迭代后没有收敛
      %  tistr = int2str(tap_it);%转化为字符串
    %  disp(['at tap iteration number ' tistr]) %显示第几次选择分接头的迭代
   else
     disp('inner load flow iterations')%否则显示（内部潮流迭代）
    disp(iter)%显示次数
   end   
   
   % if no_taps == 0
   %   lftap   %调用lftap，每次潮流都会调节，如果电压越限制，调节抽头大小如果还未调节成功，则不收敛
   % else
   %    mm_chk = 0;
   %  end
   %
%end  %如果计算分接头达到最大次数10，mm_chk=0，停止迭代


ste = clock;     % end the iteration time clock结束迭代计时

vmx_idx = find(V==volt_max);%找出最大电压
vmn_idx = find(V==volt_min);%找出最小电压

if ~isempty(vmx_idx)%如果找到
   disp('voltages at') %显示
   bus(vmx_idx,1)' 
   disp('are at the max limit')%这个节点的电压达到了最大极限
end
if ~isempty(vmn_idx)%如果找到最小的
   disp('voltages at')% 显示这个节点的电压达到了最小极限
   bus(vmn_idx,1)' 
   disp('are at the min limit');
end 
gen_index=find(bus_type==2);%chq_lim会把无功越界的PV节点改成PQ节点。  新的发电机节点
load_index = find(bus_type==3);%负荷节点
Pg(gen_index) = P(gen_index) + Pl(gen_index);%发电机有功＝计算的有功＋负荷有功
Qg(gen_index) = Q(gen_index) + Ql(gen_index);%发电机无功＝计算的无功＋负荷无功
gend_idx = find((bus(:,10)==2)&(bus_type~=2));%找出迭代后无功超出极限的发电机节点
if ~isempty(gend_idx)%如果找到
   disp('the following generators are at their var limits')%显示这些发电机节点达到了无功极限，即从PV变成了PQ的节点
   disp('    bus#    Qg')
   disp([bus(gend_idx,1)  Q(gend_idx)])
   Qlg = Ql(gend_idx)-bus(gend_idx,7);% the generator var part of the load负荷的发电机无功部分
   Qg(gend_idx)=Qg(gend_idx)-Qlg;% restore the generator vars存储发电机无功
   Ql(gend_idx)=bus(gend_idx,7);% restore the original load vars存储原始的负荷无功
end
Pl(load_index) = Pg(load_index) - P(load_index);%负荷有功＝发电机有功-计算的有功
Ql(load_index) = Qg(load_index) - Q(load_index);%负荷无功＝发电机无功-计算的无功

Pg(SB) = P(SB) + Pl(SB); Qg(SB) = Q(SB) + Ql(SB);%******************************************？？
VV = V.*exp(jay*ang);  % solution voltage 电压解
% calculate the line flows and power losses计算线路潮流和功率损耗
tap_index = find(abs(line(:,6))>0);%找出分接头变比绝对值>0的
tap_ratio = ones(nline,1);%全1阵，一列矢量
tap_ratio(tap_index)=line(tap_index,6);%设置了分抽头的都按原来的值，没有设置的为1
phase_shift(:,1) = line(:,7);%分抽头相角
tps = tap_ratio.*exp(jay*phase_shift*pi/180);%变比写成复数形式
from_bus = line(:,1);%起始节点                   % 把线路矩阵的各列赋给各个矩阵
from_int = bus_int(round(from_bus));%**********************************************？bus_int定义为global
to_bus = line(:,2);%中止节点
to_int = bus_int(round(to_bus));%***************************************************？？
r = line(:,3);%电阻
rx = line(:,4);%电抗
chrg = line(:,5);%充电电容电压
z = r + jay*rx;%阻抗
y = ones(nline,1)./z;%导纳


MW_s = VV(from_int).*conj((VV(from_int) - tps.*VV(to_int)).*y ...|VV|^2-变比×末端电压×支路导纳+起始电压×(充电电容电压/2)/|变比|^2
   + VV(from_int).*(jay*chrg/2))./(tps.*conj(tps));
P_s = real(MW_s);     % active power sent out by from_bus to to_bus从起始节点发出的注入末端节点有功
Q_s = imag(MW_s);     % reactive power sent out by from_bus to to_bus从起始节点发出的注入末端节点的无功
MW_r = VV(to_int).*conj((VV(to_int) ...%**************************************************？？？？
   - VV(from_int)./tps).*y ...
   + VV(to_int).*(jay*chrg/2));
P_r = real(MW_r);     % active power received by to_bus from from_bus末端节点收到的从起始节点来的有功
Q_r = imag(MW_r);     % reactive power received by to_bus from from_bus末端节点收到的从起始节点来的无功
iline = [1:1:nline]';%列阵，从1到nline
line_ffrom = [iline from_bus to_bus P_s Q_s];%把这次矩阵的元素全部排列起来的一个数组
line_fto   = [iline to_bus from_bus P_r Q_r];%把这次矩阵的元素全部排列起来的一个数组
% keyboard
P_loss = sum(P_s) + sum(P_r) ;%有功损耗
Q_loss = sum(Q_s) + sum(Q_r) ;%无功损耗

bus_sol=[bus_no  V  ang*180/pi Pg Qg Pl Ql Gb Bb...节点的潮流解
      bus_type qg_max qg_min bus(:,13) volt_max volt_min];
line_sol = line;%修正的线路矩阵
line_flow(1:nline, :)  =[iline from_bus to_bus P_s Q_s];
line_flow(1+nline:2*nline,:) = [iline to_bus from_bus P_r Q_r];%线路的潮流解

% Give warning of non-convergence没有收敛，给出警告
if conv_flag == 1%如果还没收敛，显示交流潮流收敛失败
   disp('ac load flow failed to converge')
   error('stop')
end

% display results显示结果
if display == 'y',
   clc
   disp('                             LOAD-FLOW STUDY')
   disp('                    REPORT OF POWER FLOW CALCULATIONS ')
   disp(' ')
   disp(date)
   fprintf('SWING BUS                  : BUS %g \n', SB)
   fprintf('NUMBER OF ITERATIONS       : %g \n', iter)
   fprintf('SOLUTION TIME              : %g sec.\n',etime(ste,st))
   fprintf('TOTAL TIME                 : %g sec.\n',etime(clock,tt))
   fprintf('TOTAL REAL POWER LOSSES    : %g.\n',P_loss)
   fprintf('TOTAL REACTIVE POWER LOSSES: %g.\n\n',Q_loss)
   if conv_flag == 0,
      disp('                                      GENERATION             LOAD')
      disp('       BUS     VOLTS     ANGLE      REAL  REACTIVE      REAL  REACTIVE ')
      disp(bus_sol(:,1:7))
      
      disp('                      LINE FLOWS                     ')
      disp('      LINE  FROM BUS    TO BUS      REAL  REACTIVE   ')
      disp(line_ffrom)
      disp(line_fto)
   end
end; 
if iter > iter_max,
   disp('Note: Solution did not converge in %g iterations.\n', iter_max)
   lf_flag = 0
end

return