fm_report.m

电力系统的psat

function fm_report
% FM_REPORT write the power 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-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

fm_var

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

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

nseries = Line.n+Ltc.n+Upfc.n+Lines.n+Sssc.n+Tcsc.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
  [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')

  Vmin = 0.8*ones(Bus.n,1);
  Vmax = 1.2*ones(Bus.n,1);
  Qgmin = -99*Settings.mva*ones(Bus.n,1);
  Qgmax =  99*Settings.mva*ones(Bus.n,1);

  if PQ.n
    Vmin(PQ.bus) = PQ.con(:,7);
    Vmax(PQ.bus) = PQ.con(:,6);
  end
  if PV.n
    Vmin(PV.bus) = PV.con(:,9);
    Vmax(PV.bus) = PV.con(:,8);
    Qgmin(PV.bus) = PV.con(:,7);
    Qgmax(PV.bus) = PV.con(:,6);
  end
  if SW.n
    Vmin(SW.bus) = SW.con(:,9);
    Vmax(SW.bus) = SW.con(:,8);
    Qgmin(SW.bus) = SW.con(:,7);
    Qgmax(SW.bus) = SW.con(:,6);
  end
  Vmin(find(Vmin == 0)) = 0.8;
  Vmax(find(Vmax == 0)) = 1.2;

  idx = find(Qgmin == 0 & Qgmax == 0);
  if ~isempty(idx)
    Qgmin(idx) = -99*Settings.mva;
    Qgmax(idx) = 99*Settings.mva;
  end

  vVmax = DAE.V > Vmax;
  vVmin = DAE.V < Vmin;
  vQgmax = Bus.Qg > Qgmax+Settings.lftol;
  vQgmin = Bus.Qg < Qgmin-Settings.lftol;

  Vmax = Vmax.*VB;
  Vmin = Vmin.*VB;
  Qgmax = Qgmax*MVA;
  Qgmin = Qgmin*MVA;

end

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

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 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(Line.to).*conj((VV(Line.to) - ...
         VV(Line.from)./tps).*y + VV(Line.to).*(jay*chrg2));
  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
  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
else
  P_s = [];
  Q_s = [];
  P_r = [];
  Q_r = [];
  from_bus = [];
  to_bus = [];
end

% 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 Upfc.n > 0
  fm_upfc(1);
  from_bus = [from_bus; Upfc.bus1];
  to_bus = [to_bus; Upfc.bus2];
  P_s = [P_s; DAE.gp(Upfc.bus1)];
  Q_s = [Q_s; -DAE.gq(Upfc.bus1)];
  P_r = [P_r; DAE.gp(Upfc.bus2)];
  Q_r = [Q_r; -DAE.gq(Upfc.bus2)];
end

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 Sssc.n > 0
  fm_sssc(1);
  from_bus = [from_bus; Sssc.bus1];
  to_bus = [to_bus; Sssc.bus2];
  P_s = [P_s; DAE.gp(Sssc.bus1)];
  Q_s = [Q_s; -DAE.gq(Sssc.bus1)];
  P_r = [P_r; DAE.gp(Sssc.bus2)];
  Q_r = [Q_r; -DAE.gq(Sssc.bus2)];
end

if Tcsc.n > 0
  from_bus = [from_bus; Tcsc.bus1];
  to_bus = [to_bus; Tcsc.bus2];
  V1 = DAE.V(Tcsc.bus1);
  V2 = DAE.V(Tcsc.bus2);
  t1 = DAE.a(Tcsc.bus1);
  t2 = DAE.a(Tcsc.bus2);
  cc = cos(t1-t2);
  P_s = [P_s; Tcsc.Pe];
  Q_s = [Q_s; V1.*(V1-V2.*cc).*Tcsc.B];
  P_r = [P_r;-Tcsc.Pe];
  Q_r = [Q_r; V2.*(V2-V1.*cc).*Tcsc.B];
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 = (sum(P_s)+sum(P_r))*MVA;
Q_loss = (sum(Q_s)+sum(Q_r))*MVA;

% Inclusion of "Supply" and "Demand" components in case of
% OPF computations
% -----------------------------------------------------------------

Pls = zeros(Bus.n,1);
Pgs = zeros(Bus.n,1);
Qls = zeros(Bus.n,1);

if OPF.init
  if Demand.n > 0
    for i = 1:Demand.n
      k = Demand.bus(i);
      Pls(k) = Demand.con(i,7);
      a = Bus.Pl(k);
      if a == 0
        a = 1;
      end
      Qls(k) = Bus.Ql(k)/a*Demand.con(i,7);
    end
    Pls = Pls(ordbus);
    Qls = Qls(ordbus);
  end
  if Supply.n > 0
    for i = 1:Supply.n
      Pgs(Supply.bus(i)) = Supply.con(i,6);
    end
    Pgs = Pgs(ordbus);
  end
end

% 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-2005';
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} = '';