fm_pareto.m

电力系统的psat

function fm_pareto
% FM_PARETO determine a Pareto set of the multi-objective
%           OPF Problem
%
% FM_PARETO
%
%see also OPF structure
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    25-Feb-2003
%Version:   1.0.1
%
%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 Fig Hdl Varname File Path DAE OPF Settings Varout
global PV SW Bus Line Demand Supply Rsrv History Theme GAMS

tempo1 = clock;

if Fig.main
  hdl = findobj(Fig.main,'Tag','PushClose');
  set(hdl,'String','Stop');
  set(Fig.main,'UserData',1);
end

fm_disp('  ')
fm_disp('PARETO SET COMPUTATIONS')

% open status bar
fm_bar open

% -------------------------------------------------------------------------
%          Indices
% -------------------------------------------------------------------------

idx_NaN = [];
idxo = 0;
ntot = length(OPF.omega);
n_gen = PV.n+SW.n;
muL = ones(Line.n,1);
% resetting time vector in case of previous time simulations
Varout.t = [];

bus_gen = [SW.bus; PV.bus];
if Rsrv.n
  n_s = 2*Supply.n+2*Demand.n+4*n_gen+4*Bus.n+3+4*Line.n+2*Rsrv.n;
  n_y = Supply.n+Demand.n+2*n_gen+4*Bus.n+2+Rsrv.n;
else
  n_s = 2*Supply.n+2*Demand.n+4*n_gen+4*Bus.n+2+4*Line.n;
  n_y = Supply.n+Demand.n+2*n_gen+4*Bus.n+2;
end
n0 = 1:(n_y+Bus.n);
a1 = 2*n_s + 2*Bus.n + Demand.n + Supply.n + n_gen;
a2 = a1 + Bus.n;
n1 = a1+1:a1+Bus.n;
n2 = a2+1:a2+Bus.n;

% -------------------------------------------------------------------------
%          Building OPF.varname
% -------------------------------------------------------------------------

switch OPF.flow
 case 1, flow = 'I_';
 case 2, flow = 'P_';
 case 3, flow = 'S_';
end

Lf = cellstr(num2str(Line.from));
Lt = cellstr(num2str(Line.to));

% -------------------------------------------------------------------------
% unformatted labels
% -------------------------------------------------------------------------
OPF.uname = strcat('theta_',{Varname.bus{:}}');
OPF.uname = [OPF.uname;strcat('V_',{Varname.bus{:}}')];
OPF.uname = [OPF.uname;strcat('Qg_',{Varname.bus{bus_gen}}')];
OPF.uname = [OPF.uname;strcat('PS_',{Varname.bus{Supply.bus}}')];
OPF.uname = [OPF.uname;strcat('PD_',{Varname.bus{Demand.bus}}')];
OPF.uname = [OPF.uname;strcat('thetac_',{Varname.bus{:}}')];
OPF.uname = [OPF.uname;strcat('Vc_',{Varname.bus{:}}')];
OPF.uname = [OPF.uname;{'kg_c'}];
OPF.uname = [OPF.uname;strcat('Qgc_',{Varname.bus{bus_gen}}')];
OPF.uname = [OPF.uname;{'lambda_c'}];
if Rsrv.n
  OPF.uname = [OPF.uname;strcat('PR_',{Varname.bus{Rsrv.bus}}')];
end
OPF.uname = [OPF.uname;strcat('LMP_',{Varname.bus{:}}')];
OPF.uname = [OPF.uname;strcat(flow,Lf,'-',Lt)];
OPF.uname = [OPF.uname;strcat(flow,Lt,'-',Lf)];
OPF.uname = [OPF.uname;strcat(flow,'c',Lf,'-',Lt)];
OPF.uname = [OPF.uname;strcat(flow,'c',Lt,'-',Lf)];
OPF.uname = [OPF.uname; {'TTL'}];

% -------------------------------------------------------------------------
% formatted labels
% -------------------------------------------------------------------------
OPF.fname = strcat('\theta_{',{Varname.bus{:}}','}');
OPF.fname = [OPF.fname;strcat('V_{',{Varname.bus{:}}','}')];
OPF.fname = [OPF.fname;strcat('Q_{g-',{Varname.bus{bus_gen}}','}')];
OPF.fname = [OPF.fname;strcat('P_{S-',{Varname.bus{Supply.bus}}','}')];
OPF.fname = [OPF.fname;strcat('P_{D-',{Varname.bus{Demand.bus}}','}')];
OPF.fname = [OPF.fname;strcat('\theta_{c-',{Varname.bus{:}}','}')];
OPF.fname = [OPF.fname;strcat('V_{c-',{Varname.bus{:}}','}')];
OPF.fname = [OPF.fname;{'k_{gc}'}];
OPF.fname = [OPF.fname;strcat('Q_{gc-',{Varname.bus{bus_gen}}','}')];
OPF.fname = [OPF.fname;{'\lambda_c'}];
if Rsrv.n
  OPF.fname = [OPF.fname;strcat('P_{R-',{Varname.bus{Rsrv.bus}}','}')];
end
OPF.fname = [OPF.fname;strcat('LMP_{',{Varname.bus{:}}','}')];
OPF.fname = [OPF.fname;strcat(flow,'{',Lf,'-',Lt,'}')];
OPF.fname = [OPF.fname;strcat(flow,'{',Lt,'-',Lf,'}')];
OPF.fname = [OPF.fname;strcat(flow,'{c',Lf,'-',Lt,'}')];
OPF.fname = [OPF.fname;strcat(flow,'{c',Lt,'-',Lf,'}')];
OPF.fname = [OPF.fname; {'TTL'}];

% -------------------------------------------------------------------------
%          Running the Pareto Set
% -------------------------------------------------------------------------

OPF.varout = zeros(ntot,n_y+Bus.n+4*Line.n+1);
OPF.wp = OPF.omega;
OPF.show = 0;
OPF.conv = 0;
idxrho = zeros(ntot,1);

lambdamax = OPF.lmax;
lambdamin = OPF.lmin;

for j = 1:ntot
  OPF.w = OPF.omega(j);
  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  OPF.lmax = lambdamax;
  OPF.lmin = lambdamin;
  eval(OPF.fun)

  if OPF.conv == 1

    OPF.lmax = OPF.guess(2*n_s+n_y);
    OPF.lmin = OPF.guess(2*n_s+n_y)-1e-5;
    OPF.w = 0;
    %disp(num2str(OPF.guess(2*n_s+n_y)))
    eval(OPF.fun)
    %fm_opfsdl

    % -------------------------------------------------------------------------
    % Local Marginal Price Computation
    % -------------------------------------------------------------------------

    % Lagrangian Multiplier of Power Flow Equations
    rho = -OPF.guess(2*n_s+n_y+1:end-3*Bus.n)';

    % Alternative Formulas for Lagrangian Multipliers of PF Eqs.
    %Cs = Supply.con(:,8);
    %Cd = Demand.con(:,9);

    %muPsmax = OPF.guess(n_s+Supply.n+1:n_s+2*Supply.n);
    %muPsmin = OPF.guess(n_s+1:n_s+Supply.n);
    %rhoPsc = OPF.guess(2*n_s+n_y+2*Bus.n+Supply.bus);

    %muPdmax = OPF.guess(n_s+2*Supply.n+Demand.n+1:n_s+2*Supply.n+2*Demand.n);
    %muPdmin = OPF.guess(n_s+2*Supply.n+1:n_s+2*Supply.n+Demand.n);
    %rhoPdc = OPF.guess(2*n_s+n_y+2*Bus.n+Demand.bus);
    %rhoQdc = OPF.guess(2*n_s+n_y+3*Bus.n+Demand.bus);
    %rhoQd  = OPF.guess(2*n_s+n_y+Bus.n+Demand.bus);

    %tanphi = Demand.con(:,4)./Demand.con(:,3);
    %lambda = OPF.guess(2*n_s+n_y);
    %kg = OPF.guess(2*n_s+n_y-n_gen-1);

    %rhoS = Cs +(muPsmax - muPsmin + rhoPsc*(1 + lambda + kg));
    %rhoD = Cd +(muPdmin - muPdmax + (rhoPdc + rhoQdc.*tanphi)* ...
    %       (1 + lambda) + rhoQd.*tanphi);

    %rho = [rhoS;rhoD]';
    % -------------------------------------------------------------------------
    % End of Local Marginal Price Computation
    % -------------------------------------------------------------------------

    if abs(rho) < 1e-3, idxrho(j) = j; end
    OPF.varout(j,n0) = [OPF.guess(2*n_s+1:end-4*Bus.n)',rho];
    [Iij, Jij, Hij, Iji, Jji, Hji] = fm_flows(OPF.flow, muL, muL);
    OPF.varout(j,n_y+Bus.n+1:n_y+Bus.n+2*Line.n) = sqrt([Iij', Iji']);
    V_snap = DAE.V;
    ang_snap = DAE.a;
    DAE.V  = OPF.guess(n2);
    DAE.a    = OPF.guess(n1);
    [Iij, Jij, Hij, Iji, Jji, Hji] = fm_flows(OPF.flow, muL, muL);
    OPF.varout(j,n_y+Bus.n+2*Line.n+1:n_y+Bus.n+4*Line.n) = sqrt([Iij', ...
                        Iji']);
    DAE.V = V_snap;
    DAE.a = ang_snap;
    OPF.varout(j,end) = ...
        Settings.mva*sum(OPF.guess(2*n_s+2*Bus.n+n_gen+Supply.n+1:2*n_s+ ...
                                   2*Bus.n+Supply.n+n_gen+Demand.n));
  else
    idx_NaN = [idx_NaN, j];
  end
  idx = j/ntot;
  fm_bar([idxo,idx])
  idxo = idx;
  History.text{end} = [History.text{end}, '  omega = ', ...
                      num2str(OPF.omega(j))];
end

OPF.wp(idx_NaN) = [];
OPF.varout(idx_NaN,:) = [];
idxrho(idx_NaN) = [];

if ~isempty(find(idxrho))
  for i = 2:length(idxrho)
    if idxrho(i) ~= 0 & idxrho(i) ~= length(idxrho)
      OPF.varout(idxrho(i),:) = 0.5*(OPF.varout(idxrho(i)+1,:)+ ...
                                     OPF.varout(idxrho(i)-1,:));
    end
  end
end

% -------------------------------------------------------------------------
%          Graphical settings
% -------------------------------------------------------------------------

% close status bar
fm_bar close
% plot settings
Settings.xlabel = 'weighting factor \omega';
GAMS.hours = [];

if Fig.main
  set(hdl,'String','Close');
  if ~get(Fig.main,'UserData'),
    fm_disp(['Pareto Set Computation interrupted.'],1),
    return
  end
end

fm_disp(['Pareto Set Computation completed in ', ...
         num2str(etime(clock,tempo1)),' s'],1)
OPF.init = 3;
OPF.lmax = lambdamax;
OPF.lmin = lambdamin;