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://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

fm_var

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

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

nseries = Line.n+Ltc.n+Lines.n+Hvdc.n+Phs.n;
iline = [1:1:nseries]';
if Fig.stat
  hdlT = findobj(Fig.stat,'Tag','PushSort');
  string = get(hdlT,'UserData');
  switch string
   case '1n'
    [buss,ordbus] = sort(Varname.bus(1:end));
    [buss,ordbus] = sort(Bus.con(:,1));
   case 'az'
    [buss,ordbus] = sort(Varname.bus(1:end));
   otherwise
    % nothing to do for now
  end
else
  % for command line usage, alphabetical order is used
  [buss,ordbus] = sort(Varname.bus(1:end));
end
%[buss,ordbus] = sort(Bus.con(:,1));
nomi_bus = Varname.bus;
tab = '****  ';
space = repmat(' ',1,11);
mgp = max(abs(DAE.gp));
mgq = max(abs(DAE.gq));

checkabs = Settings.absvalues;
checkshunt = Settings.shuntvalues;
violations = Settings.violations;

switch checkabs
 case 'on'
  MVA  = Settings.mva;
  VB   = Bus.con(:,2);
  MW   = '[MW]';
  MVar = '[MVar]';
  kV   = '[kV]';
 otherwise
  MVA  = 1;
  VB   = ones(Bus.n,1);
  MW   = '[p.u.]';
  MVar = '[p.u.]';
  kV   = '[p.u.]';
end

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

if strcmp(violations,'on')

  [Vmax,Vmin] = fm_vlim(1.2,0.8);
  [Qgmax,Qgmin] = fm_qlim(99,-99,'all');

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

end

% flows in the "Line" components
% -----------------------------------------------------------------

P_s = [];
Q_s = [];
P_r = [];
Q_r = [];
from_bus = [];
to_bus = [];

if Line.n  
  tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
  VV = DAE.V.*exp(jay*DAE.a);
  r = Line.con(:,8);
  rx = Line.con(:,9);
  chrg = Line.con(:,10)/2;
  if strcmp(checkshunt,'on')
    chrg1 = chrg - jay*Shunt.g(Line.from) + Shunt.b(Line.from);
    chrg2 = chrg - jay*Shunt.g(Line.to) + Shunt.b(Line.to);
  else
    chrg1 = chrg;
    chrg2 = chrg;
  end
  z = r + jay*rx;
  y = ones(Line.n,1)./z;
  MW_s = VV(Line.from).*conj((VV(Line.from) - tps.*VV(Line.to)).*y ...
         + VV(Line.from).*(jay*chrg1))./(tps.*conj(tps));
  P_s = real(MW_s); % active power P_sr
  Q_s = imag(MW_s); % reactive power Q_sr
  MW_r = VV(Line.to).*conj((VV(Line.to) - ...
         VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg2));
  P_r = real(MW_r); % active power P_rs 
  Q_r = imag(MW_r); % reactive power Q_rs
  from_bus = Line.from;
  to_bus = Line.to;
  if strcmp(violations,'on')
    Ss = abs(MW_s);
    Sr = abs(MW_r);
    Is = Ss./DAE.V(Line.from);
    Ir = Sr./DAE.V(Line.to);
    Ps = abs(P_s);
    Pr = abs(P_r);
    Imax = Line.con(:,13);
    idx = find(Imax == 0);
    if ~isempty(idx)
      Imax(idx) = 1e6*Settings.mva;
    end
    Pmax = Line.con(:,14);
    idx = find(Pmax == 0);
    if ~isempty(idx)
      Pmax(idx) = 1e6*Settings.mva;
    end
    Smax = Line.con(:,15);
    idx = find(Smax == 0);
    if ~isempty(idx)
      Smax(idx) = 1e6*Settings.mva;
    end
    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.from)/MVA/1000;
      Imaxr = sqrt(3)*Imax.*VB(Line.to)/MVA/1000;
      Is = sqrt(3)*Is.*VB(Line.from)/MVA/1000;
      Ir = sqrt(3)*Ir.*VB(Line.to)/MVA/1000;
    end
  end
end

[P_r,P_s,Q_r,Q_s] = flows(Tcsc,P_r,P_s,Q_r,Q_s);
[P_r,P_s,Q_r,Q_s] = flows(Sssc,P_r,P_s,Q_r,Q_s);
[P_r,P_s,Q_r,Q_s] = flows(Upfc,P_r,P_s,Q_r,Q_s);

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

if Ltc.n > 0
  from_bus = [Line.from; Bus.int(Ltc.con(:,1))];
  to_bus = [Line.to; Bus.int(Ltc.con(:,2))];
  P_s = [P_s; -real(Ltc.dat(:,5))];
  Q_s = [Q_s; imag(Ltc.dat(:,5))];
  P_r = [P_r; -real(Ltc.dat(:,6))];
  Q_r = [Q_r; imag(Ltc.dat(:,6))];
end

gp_old = DAE.gp;
gq_old = DAE.gq;

DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
if Lines.n > 0
  fm_lines(1);
  from_bus = [from_bus; Lines.bus1];
  to_bus = [to_bus; Lines.bus2];
  P_s = [P_s; DAE.gp(Lines.bus1)];
  Q_s = [Q_s; -DAE.gq(Lines.bus1)];
  P_r = [P_r; DAE.gp(Lines.bus2)];
  Q_r = [Q_r; -DAE.gq(Lines.bus2)];
end

DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
if Hvdc.n > 0
  fm_hvdc(1);
  from_bus = [from_bus; Hvdc.bus1];
  to_bus = [to_bus; Hvdc.bus2];
  P_s = [P_s; DAE.gp(Hvdc.bus1)];
  Q_s = [Q_s; -DAE.gq(Hvdc.bus1)];
  P_r = [P_r; DAE.gp(Hvdc.bus2)];
  Q_r = [Q_r; -DAE.gq(Hvdc.bus2)];
end

DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
if Phs.n > 0
  fm_phs(1);
  from_bus = [from_bus; Phs.bus1];
  to_bus = [to_bus; Phs.bus2];
  P_s = [P_s; DAE.gp(Phs.bus1)];
  Q_s = [Q_s; -DAE.gq(Phs.bus1)];
  P_r = [P_r; DAE.gp(Phs.bus2)];
  Q_r = [Q_r; -DAE.gq(Phs.bus2)];
end

DAE.gp = gp_old;
DAE.gq = gq_old;

line_ffrom = [iline from_bus to_bus P_s*MVA Q_s*MVA];
line_fto   = [iline to_bus from_bus P_r*MVA Q_r*MVA];

% Computation of total losses
% -----------------------------------------------------------------

P_loss = MVA*sum(P_s+P_r);
Q_loss = MVA*sum(Q_s+Q_r);

% Creation of the text file and report writing
% ---------------------------------------------------------------

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-2006';
Header{1,1}{6,1} = 'e-mail:  [e-mail]';
Header{1,1}{7,1} = 'website: http://thunderbox.uwaterloo.ca/~fmilano';
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.V(ordbus).*VB(ordbus);
angs = DAE.a(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;
if Line.n,
  ntrasf = length(find(Line.con(:,7) ~= 0));
else
  ntrasf = 0;
end
nline = Lines.n + Line.n-ntrasf;

% 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 ntrasf > 0
  idx = idx + 1;
  Matrix{2,1}(idx,1) = ntrasf;
  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]};
mgp = max(abs(DAE.gp));
mgq = max(abs(DAE.gq));
Matrix{3,1} = [Settings.iter; mgp*MVA; mgq*MVA];
if strcmp(checkabs,'off')
  Matrix{3,1}(4,1) = Settings.mva;
  Rows{3,1}{4,1} = 'Power rate [MVA]';