fm_report.m

电力系统分析计算程序

function fm_report
% FM_REPORT write the power flow report.
%
% The report is saved in a text file with the same name of
% the data file followed by "_xx" where xx is a progressive
% number. Finally the report is displayed using the currently
% selected text viewer.
%
% FM_REPORT
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    23-May-2003
%Update:    24-Aug-2003
%Update:    14-Sep-2003
%Version:   2.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

fm_var

if ~Settings.init
  fm_disp('Solve power flow before writing the report file')
end

% General variables
% --------------------------------------------------------------------

nseries = Settings.nseries;
iline = [1:nseries]';
if Fig.stat
  hdlT = findobj(Fig.stat,'Tag','PushSort');
  string = get(hdlT,'UserData');
  switch string
   case '1n'
    [buss,ordbus] = sort(Bus.names(1:end));
    [buss,ordbus] = sort(getidx(Bus,0));
   case 'az'
    [buss,ordbus] = sort(Bus.names(1:end));
   otherwise
    % nothing to do for now
  end
else
  % for command line usage, alphabetical order is used
  [buss,ordbus] = sort(Bus.names(1:end));
end
%[buss,ordbus] = sort(Bus.con(:,1));
nomi_bus = Bus.names;
tab = '****  ';
space = repmat(' ',1,11);
checkabs = Settings.absvalues;
violations = Settings.violations;

switch checkabs
 case 'on'
  MVA  = Settings.mva;
  VB   = getkv(Bus,0,0);
  MW   = '[MW]';
  MVar = '[MVar]';
  kV   = '[kV]';
 otherwise
  MVA  = 1;
  VB   = getones(Bus);
  MW   = '[p.u.]';
  MVar = '[p.u.]';
  kV   = '[p.u.]';
end

% check voltage and generator reactive power limits if necessary
% --------------------------------------------------------------------

if strcmp(violations,'on')

  Vbus = DAE.y(Bus.v);
  [Vmax,Vmin] = fm_vlim(1.2,0.8);
  [Qgmax,Qgmin] = fm_qlim('all');
  buses = getzeros(Bus);
  buses(getbus(PV)) = 1;
  buses(getbus(SW)) = 1;

  vVmax = Vbus > Vmax;
  vVmin = Vbus < Vmin;
  vVminabs = abs(Vbus-Vmin) < Settings.lftol;
  vVmaxabs = abs(Vbus-Vmax) < Settings.lftol;
  vQgmax = Bus.Qg > (Qgmax+Settings.lftol);
  vQgmin = Bus.Qg < (Qgmin-Settings.lftol);
  vQgminabs = (abs(Bus.Qg-Qgmin) < Settings.lftol) & buses;
  vQgmaxabs = (abs(Bus.Qg-Qgmax) < Settings.lftol) & buses;
  Vmax = Vmax.*VB;
  Vmin = Vmin.*VB;
  Qgmax = Qgmax*MVA;
  Qgmin = Qgmin*MVA;

end

% flows in the series components
% --------------------------------------------------------------------

[P_s,Q_s,P_r,Q_r,fr_bus,to_bus] = fm_flows('bus');

line_ffr = [iline, fr_bus, to_bus, P_s*MVA, Q_s*MVA];
line_fto = [iline, to_bus, fr_bus, P_r*MVA, Q_r*MVA];

ntraf = transfno(Line) + Ltc.n + Phs.n;
nline = Hvdc.n + Lines.n + Line.n - ntraf;

% check transmission line limits
% --------------------------------------------------------------------

if strcmp(violations,'on') & Line.n
  Ss = sqrt(P_s(1:Line.n).^2+Q_s(1:Line.n).^2);
  Sr = sqrt(P_r(1:Line.n).^2+Q_r(1:Line.n).^2);
  Is = Ss./DAE.y(Line.vfr);
  Ir = Sr./DAE.y(Line.vto);
  Ps = abs(P_s);
  Pr = abs(P_r);
  Imax = getflowmax(Line,'imax');
  Pmax = getflowmax(Line,'pmax');
  Smax = getflowmax(Line,'smax');
  vIs = Is > Imax;
  vPs = Ps > Pmax;
  vSs = Ss > Smax;
  vIr = Ir > Imax;
  vPr = Pr > Pmax;
  vSr = Sr > Smax;
  Pmax = Pmax*MVA;
  Smax = Smax*MVA;
  Ps = Ps*MVA;
  Pr = Pr*MVA;
  Ss = Ss*MVA;
  Sr = Sr*MVA;
  Imaxs = Imax;
  Imaxr = Imax;
  if strcmp(checkabs,'on')
    Imaxs = sqrt(3)*Imax.*VB(Line.fr)/MVA/1000;
    Imaxr = sqrt(3)*Imax.*VB(Line.to)/MVA/1000;
    Is = sqrt(3)*Is.*VB(Line.fr)/MVA/1000;
    Ir = sqrt(3)*Ir.*VB(Line.to)/MVA/1000;
  end
end

% Power flow report
% --------------------------------------------------------------------

fm_disp
fm_disp('Writing the report file...')

% initialization of report outputs
% --------------------------------------------------------------------

Header = cell(0);
Matrix = cell(0);
Cols = cell(0);
Rows = cell(0);

% general header
% --------------------------------------------------------------------

if OPF.init
  Header{1,1}{1,1} = 'OPTIMAL POWER FLOW REPORT';
elseif CPF.init
  Header{1,1}{1,1} = ['CONTINUATION POWER FLOW REPORT (last ' ...
                      'corrector step solution)'];
elseif LIB.init
  Header{1,1}{1,1} = ['POWER FLOW REPORT (limit-induced bifurcation ' ...
                      'results)'];
elseif SNB.init
  Header{1,1}{1,1} = ['OPTIMAL POWER FLOW REPORT (Saddle-node ' ...
                      'bifurcation results)'];
else
  Header{1,1}{1,1} = 'POWER FLOW REPORT';
end
Header{1,1}{2,1} = ' ';
Header{1,1}{3,1} = ['P S A T  ',Settings.version];
Header{1,1}{4,1} = ' ';
Header{1,1}{5,1} = 'Author:  Federico Milano, (c) 2002-2008';
Header{1,1}{6,1} = 'e-mail:  [e-mail]';
Header{1,1}{7,1} = 'website: http://www.uclm.es/area/gsee/Web/Federico';
Header{1,1}{8,1} = ' ';
Header{1,1}{9,1} = ['File:  ', Path.data,strrep(File.data,'(mdl)','.mdl')];
Header{1,1}{10,1} = ['Date:  ',datestr(now,0)];

Matrix{1,1} = [];
Cols{1,1} = '';
Rows{1,1} = '';

% some useful variables
% --------------------------------------------------------------------

Vs = DAE.y(ordbus+Bus.n).*VB(ordbus);
angs = DAE.y(ordbus);
raddeg = 'rad';
if Fig.stat
  hdlT = findobj(Fig.stat,'Tag','PushAngle');
  string = get(hdlT,'UserData');
  if ~strcmp(string,'rad')
    angs = angs*180/pi;
    raddeg = 'deg';
  end
end

Pgs = Bus.Pg(ordbus)*MVA;
Qgs = Bus.Qg(ordbus)*MVA;
Pls = Bus.Pl(ordbus)*MVA;
Qls = Bus.Ql(ordbus)*MVA;

% losses
% --------------------------------------------------------------------

line_p_losses = line_ffr(:,4) + line_fto(:,4);
line_q_losses = line_ffr(:,5) + line_fto(:,5);

% network statistics
% --------------------------------------------------------------------

Header{2,1} = 'NETWORK STATISTICS';
Matrix{2,1} = Bus.n;
Rows{2,1} = {'Buses:'};
Cols{2,1} = '';
if nline > 0
  Matrix{2,1}(2,1) = nline;
  Rows{2,1}{2,1} = 'Lines:';
  idx = 2;
else
  idx = 1;
end
if ntraf > 0
  idx = idx + 1;
  Matrix{2,1}(idx,1) = ntraf;
  Rows{2,1}{idx,1} = 'Transformers:';
end
Matrix{2,1}(idx+1,1) = SW.n+PV.n+Syn.n;
Rows{2,1}{idx+1,1} = 'Generators:';
Matrix{2,1}(idx+2,1) = PQ.n+Pl.n+Mn.n;
Rows{2,1}{idx+2,1} = 'Loads:';

% statistics of the current solution algorithm
% --------------------------------------------------------------------

Header{3,1} = 'SOLUTION STATISTICS';
Cols{3,1} = '';
Rows{3,1} = {'Number of Iterations:'; ...
	     ['Maximum P mismatch ',MW]; ...
	     ['Maximum Q mismatch ',MVar]};
mp = max(abs(DAE.g(Bus.a)));
mq = max(abs(DAE.g(Bus.v)));
Matrix{3,1} = [Settings.iter; mp*MVA; mq*MVA];
if strcmp(checkabs,'off')
  Matrix{3,1}(4,1) = Settings.mva;
  Rows{3,1}{4,1} = 'Power rate [MVA]';
end

% power flow results
% --------------------------------------------------------------------

Header{4,1} = 'POWER FLOW RESULTS';

if Settings.report % embed line flow in bus report

  nh = 4;
  for i = 1:Bus.n
    % bus report
    nh = nh + 1;
    %Header{nh,1} = nomi_bus{ordbus(i)};
    Header{nh,1} = '';
    bus_name = nomi_bus{ordbus(i)};
    Cols{nh,1} = {bus_name,'V','phase','P gen','Q gen','P load','Q load'; ...
                 ' ', kV, ['[',raddeg,']'],MW,MVar,MW,MVar};
    Rows{nh,1} = {'  '};
    Matrix{nh,1} = [Vs(i),angs(i),Pgs(i),Qgs(i),Pls(i),Qls(i)];

    % violations
    if strcmp(violations,'on')
      hhh = vVmin(i) + vVmax(i) + vVmaxabs(i) + vVminabs(i) + vQgmax(i) ...
            + vQgmin(i) + vQgmaxabs(i) + vQgminabs(i);
      if hhh
        nh = nh + 1;
        Header{nh} = cell(hhh,1);
      end
      hh = 1;
      if vVmin(i)
        Header{nh,1}{hh,1} = sprintf('         *  Minimum voltage limit violation [V_min = %g]',Vmin(i));
        hh = hh + 1;
      end
      if vVmax(i)
        Header{nh,1}{hh,1} = sprintf('         *  Maximum voltage limit violation [V_max = %g]',Vmax(i));
        hh = hh + 1;
      end
      if vVmaxabs(i)
        Header{nh,1}{hh,1} = sprintf('         *  Maximum voltage [V_max = %g]',Vmax(i));
        hh = hh + 1;
      end
      if vVminabs(i)
        Header{nh,1}{hh,1} = sprintf('         *  Minimum voltage [V_min = %g]',Vmin(i));
        hh = hh + 1;
      end
      if vQgmax(i)
        Header{nh,1}{hh,1} = sprintf('         *  Maximum reactive power limit violation [Qg_max = %g]',Qgmax(i));
        hh = hh + 1;
      end
      if vQgmin(i)
        Header{nh,1}{hh,1} = sprintf('         *  Minimum reactive power limit violation [Qg_min = %g]',Qgmin(i));
        hh = hh + 1;
      end
      if vQgmaxabs(i)
        Header{nh,1}{hh,1} = sprintf('         *  Maximum reactive power [Qg_max = %g]',Qgmax(i));
        hh = hh + 1;
      end
      if vQgminabs(i)
        Header{nh,1}{hh,1} = sprintf('         *  Minimum reactive power [Qg_min = %g]',Qgmin(i));
      end
    end
    
    % flows 
    nh = nh + 1;
    Header{nh,1} = '';
    Cols{nh,1} = {'  ','To Bus','Line','P Flow', ...
                  'Q Flow','P Loss','Q Loss'; ...
                  ' ',' ',' ',MW,MVar,MW,MVar};

    idx_to = find(line_ffr(:,2) == ordbus(i));
    idx_fr = find(line_fto(:,2) == ordbus(i));
    idx_ln = [line_ffr(idx_to,3);line_fto(idx_fr,3)];
    Rows{nh,1} = cell(length(idx_ln),2);
    for iii = 1:length(idx_ln)
      Rows{nh,1}{iii,1} = '  ';
      Rows{nh,1}{iii,2} = nomi_bus{idx_ln(iii),1};
    end
    Matrix{nh,1} = [[line_ffr(idx_to,1); line_fto(idx_fr,1)], ...
                    [line_ffr(idx_to,4); line_fto(idx_fr,4)], ...
                    [line_ffr(idx_to,5); line_fto(idx_fr,5)], ...
                    line_p_losses([idx_to;idx_fr]), ...