fm_opfrep.m

一个较好的MATLAB潮流程序

function fm_opfrep
%FM_OPFREP writes OPF report file
%
%see also FM_OPFM FM_OPFSD and OPF structure for settings.
%
%Author:    Federico Milano
%Date:      09-Mar-2005
%Version:   1.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.

global OPF DAE Bus Varname Settings File Path
global Line SW PV PQ Demand Supply Rsrv Snapshot

if OPF.init == 0
  fm_disp('Solve OPF before writing the OPF report.',2)
end

if OPF.init > 2
  fm_disp('OPF report is not defined for the current OPF solution.',2)
end

% some useful data and indexes
MVA = Settings.mva;

Psmax = Supply.con(:,4);
Psmin = Supply.con(:,5);
if Demand.n 
  Pdmax = Demand.con(:,5); 
  Pdmin = Demand.con(:,6); 
  noDem = 0;
else
  noDem = 1;
  Demand.n = 1;
end
if Rsrv.n
  Prmax = Rsrv.con(:,3);
  Prmin = Rsrv.con(:,4);
end

n_gen = PV.n+SW.n;
busG = [SW.bus; PV.bus];
nS = 1:Supply.n;
nD = 1:Demand.n;
nG = 1:n_gen;
nB = 1:Bus.n;
nL = 1:Line.n;
nR = 1:Rsrv.n;

Vmin = OPF.vmin*ones(Bus.n,1);
Vmax = OPF.vmax*ones(Bus.n,1);
if OPF.envolt
  Vmin(PQ.bus) = PQ.con(:,7);
  Vmax(PQ.bus) = PQ.con(:,6);
  Vmin(PV.bus) = PV.con(:,9);
  Vmax(PV.bus) = PV.con(:,8);
  Vmin(SW.bus) = SW.con(:,9);
  Vmax(SW.bus) = SW.con(:,8);
  Vmin(find(Vmin == 0)) = OPF.vmin;
  Vmax(find(Vmax == 0)) = OPF.vmax;
end

if OPF.enflow
  Iijmax = Line.con(:,12+OPF.flow);
  idx_no_I_max = find(Iijmax == 0);
  % no limit for undefined limit
  Iijmax(idx_no_I_max) = 999;
else
  Iijmax = 999*ones(Line.n,1);
end

if OPF.enreac
  Qgmin = [SW.con(:,7); PV.con(:,7)];
  Qgmax = [SW.con(:,6); PV.con(:,6)];
  idx = find(Qgmin == 0 & Qgmax == 0);
  if ~isempty(idx)
    Qgmin(idx) = -99*MVA;
    Qgmax(idx) = 99*MVA;
  end
  if OPF.flatstart == 1
    Bus.Qg(busG) = max(Bus.Qg(busG),Qgmin+1e-3);
    Bus.Qg(busG) = min(Bus.Qg(busG),Qgmax-1e-3);
  end
else
  Qgmin = -99*MVA*ones(n_gen,1);
  Qgmax =  99*MVA*ones(n_gen,1);
end

n2 = 2*Bus.n;
n3 = 3*Bus.n;
n4 = 4*Bus.n;
n5 = 5*Bus.n;
n6 = 6*Bus.n;
n7 = 7*Bus.n;

% numbering
if OPF.init == 1

  n_s = 2*Supply.n+2*Demand.n+2*n_gen+2*Bus.n+2*Line.n+2*Rsrv.n;
  if Rsrv.n, n_s = n_s + 1; end

  mu_Psmin  = OPF.guess(n_s+nS);       idx = Supply.n;
  mu_Psmax  = OPF.guess(n_s+idx+nS);   idx = idx + Supply.n;
  mu_Pdmin  = OPF.guess(n_s+idx+nD);   idx = idx + Demand.n;
  mu_Pdmax  = OPF.guess(n_s+idx+nD);   idx = idx + Demand.n;
  mu_Qgmin  = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Qgmax  = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Vmin   = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_Vmax   = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_Iijmax = OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  mu_Ijimax = OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  if Rsrv.n
    mu_Prmin  = OPF.guess(n_s+idx+nR); idx = idx + Rsrv.n;
    mu_Prmax  = OPF.guess(n_s+idx+nR); idx = idx + Rsrv.n;
    mu_sumPrd = OPF.guess(n_s+idx+1);
  end

  a = OPF.guess(2*n_s+nB);
  V = OPF.guess(2*n_s+Bus.n+nB);
  Qg = OPF.guess(2*n_s+n2+nG);
  Ps = OPF.guess(2*n_s+n2+n_gen+nS);
  Pd = OPF.guess(2*n_s+n2+n_gen+Supply.n+nD);
  if Rsrv.n
    Pr = OPF.guess(2*n_s+n2+n_gen+Supply.n+Demand.n+nR);
  end
  rhop = OPF.guess(2*n_s+n2+n_gen+Supply.n+Demand.n+Rsrv.n+nB);
  rhoq = OPF.guess(2*n_s+n3+n_gen+Supply.n+Demand.n+Rsrv.n+nB);
  
elseif OPF.init == 2

  n_s = 2*Supply.n+2*Demand.n+4*n_gen+4*Bus.n+2+4*Line.n+2*Rsrv.n;
  if Rsrv.n, n_s = n_s + 1; end
  
  mu_Psmin  = OPF.guess(n_s+nS);       idx = Supply.n;
  mu_Psmax  = OPF.guess(n_s+idx+nS);   idx = idx + Supply.n;
  mu_Pdmin  = OPF.guess(n_s+idx+nD);   idx = idx + Demand.n;
  mu_Pdmax  = OPF.guess(n_s+idx+nD);   idx = idx + Demand.n;
  mu_Qgmin  = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Qgmax  = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Vmin   = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_Vmax   = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_Qgcmin = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Qgcmax = OPF.guess(n_s+idx+nG);   idx = idx + n_gen;
  mu_Vcmin  = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_Vcmax  = OPF.guess(n_s+idx+nB);   idx = idx + Bus.n;
  mu_lcmin  = OPF.guess(n_s+idx+1);    idx = idx + 1;
  mu_lcmax  = OPF.guess(n_s+idx+1);    idx = idx + 1;
  mu_Iijmax = OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  mu_Iijcmax= OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  mu_Ijimax = OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  mu_Ijicmax= OPF.guess(n_s+idx+nL);   idx = idx + Line.n;
  if Rsrv.n
    mu_Prmin  = OPF.guess(n_s+idx+nR); idx = idx + Rsrv.n;
    mu_Prmax  = OPF.guess(n_s+idx+nR); idx = idx + Rsrv.n;
    mu_sumPrd = OPF.guess(n_s+idx+1);
  end

  a = OPF.guess(2*n_s+nB);
  V = OPF.guess(2*n_s+Bus.n+nB);
  Qg = OPF.guess(2*n_s+n2+nG);
  Ps = OPF.guess(2*n_s+n2+n_gen+nS);
  Pd = OPF.guess(2*n_s+n2+n_gen+Supply.n+nD);
  ac = OPF.guess(2*n_s+n2+n_gen+Supply.n+Demand.n+nB);
  Vc = OPF.guess(2*n_s+n3+n_gen+Supply.n+Demand.n+nB);
  kg = OPF.guess(2*n_s+n4+n_gen+Supply.n+Demand.n+1);
  Qgc = OPF.guess(2*n_s+n4+n_gen+Supply.n+Demand.n+1+nG);
  lambda = OPF.guess(2*n_s+n4+n_gen+Supply.n+Demand.n+2+n_gen);
  if Rsrv.n
    Pr = OPF.guess(2*n_s+n4+2*n_gen+Supply.n+Demand.n+2+nR);
  end
  rhop = OPF.guess(2*n_s+n4+2*n_gen+Supply.n+Demand.n+2+Rsrv.n+nB);
  rhoq = OPF.guess(2*n_s+n5+2*n_gen+Supply.n+Demand.n+2+Rsrv.n+nB);
  rhopc = OPF.guess(2*n_s+n6+2*n_gen+Supply.n+Demand.n+2+Rsrv.n+nB);
  rhoqc = OPF.guess(2*n_s+n7+2*n_gen+Supply.n+Demand.n+2+Rsrv.n+nB);
  
end

Pay = rhop.*DAE.glfp*MVA;
ISOPay = sum(Pay);
[Iij, Jij, Hij, Iji, Jji, Hji] = fm_flows(OPF.flow,mu_Iijmax, mu_Ijimax);
Iij  = sqrt(Iij);
Iji  = sqrt(Iji);

if OPF.init == 2
  vold = DAE.V;
  aold = DAE.a;
  DAE.V = Vc;
  DAE.a = ac;
  if OPF.line
    line_orig = Line.con;
    Y_orig = Line.Y;
    Line.con(OPF.line,[8 9 10]) = [0, 1e5, 0];
    line_cont = Line.con;
    fm_y;
  end
  [Iijc, Jijc, Hijc, Ijic, Jjic, Hjic] = fm_flows(OPF.flow,mu_Iijcmax, mu_Ijicmax);
  if OPF.line
    Line.Y = Y_orig; 
    Line.con = line_orig; 
  end  
  Iijc = sqrt(Iijc);
  Ijic = sqrt(Ijic);
end

% initialize report structures
Header{1,1}{1,1} = 'OPTIMAL POWER FLOW REPORT';
switch OPF.init
 case 1, Header{1,1}{2,1} = '(Standard OPF)';
 case 2, Header{1,1}{2,1} = '(Multiobjective OPF)';
 otherwise, % nothing to do ...
end
Header{1,1}{3,1} = ' ';
Header{1,1}{4,1} = ['P S A T  ',Settings.version];
Header{1,1}{5,1} = ' ';
Header{1,1}{6,1} = 'Author:  Federico Milano, (c) 2002-2005';
Header{1,1}{7,1} = 'e-mail:  [e-mail]';
Header{1,1}{8,1} = 'website: http://thunderbox.uwaterloo.ca/~fmilano';
Header{1,1}{9,1} = ' ';
Header{1,1}{10,1} = ['File:  ', Path.data,strrep(File.data,'(mdl)','.mdl')];
Header{1,1}{11,1} = ['Date:  ',datestr(now,0)];

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

Header{2,1} = 'NETWORK STATISTICS';
Matrix{2,1} = Bus.n;
Rows{2,1} = {'Buses:'};
Cols{2,1} = '';

if Line.n,
  ntrasf = length(find(Line.con(:,7) ~= 0));
else
  ntrasf = 0;
end
nline = Line.n-ntrasf;

if nline > 0
  Matrix{2,1}(2,1) = nline;
  Rows{2,1}{2,1} = 'Lines:';
  idx = 2;