fm_opfsdr.m

电力系统的psat

function fm_opfsdr
%FM_OPFSDR solve the OPF-based  electricity market problem by means of
%          an Interior Point Method with a Merhotra Predictor-Corrector
%          or Newton direction technique.
%          Voltage stability limits are included (see equations below).
%
%System equations:
%
%Min:  (1-w)*(Cs'*Ps - Cd'Pd + Cr*Pr) - w*lc
%
%s.t.: f(theta,V,Qg,Ps,Pd) = 0           (PF eq.)
%      f(thetac,Vc,Qgc,Ps,Pd,kg,lc) = 0  (PF eq. crit.)
%      lc_max >= lc >= lc_min            (min. stab. margin)
%      Ps_min <= Ps <= Ps_max            (supply bid blocks)
%      Pd_min <= Pd <= Pd_max            (demand bid blocks)
%  *   Ps + Pr <= Ps_max                 (reserve blocks)
%  *   Pr_min <= Pr <= Pr_max
%  *   sum(Pr) <= sum(Pd)
%      |Iij(theta,V)|   <= I_max         (thermal or power limits)
%      |Iij(thetac,Vc)| <= I_max
%      |Iji(theta,V)|   <= I_max
%      |Iji(thetac,Vc)| <= I_max
%      Qg_min  <= Qg  <= Qg_max          (gen. Q limits)
%      Qgc_min <= Qgc <= Qgc_max
%      V_min  <= V  <= V_max             (bus DAE.V limits)
%      Vc_min <= Vc <= Vc_max
%
%(* optional constraints)
%
%see also FM_OPFM FM_PARETO FM_ATC and the OPF structure for settings.
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%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.

fm_var

if ~autorun('Optimal Power Flow')
  return
end

if ~Supply.n,
  fm_disp('Supply data have to be specified before in order to run OPF routines',2)
  return
end

lcmin = OPF.lmin;
lcmax = OPF.lmax;
w = OPF.w;
if w > 1 | w < 0
  fm_disp('The weigthing factor range is [0,1]')
  return
end

if OPF.show
  fm_disp
  fm_disp('----------------------------------------------------------------')
  fm_disp('Interior Point Method for OPF Computation')
  fm_disp('Spot Price of Security (including loading parameter)')
  if w > 0 & w < 1,
    fm_disp('Multi-Objective Function')
  elseif w == 0,
    fm_disp('Social Benefit Objective Function')
  elseif w == 1,
    fm_disp('Maximization of the Distance to Collapse')
  end
  fm_disp(['Minimum Loading Parameter = ',num2str(lcmin)])
  fm_disp(['Weighting Factor = ',num2str(w)])
  fm_disp(['Data file "',Path.data,File.data,'"'])
  fm_disp('----------------------------------------------------------------')
  fm_disp
end

tic;
if w == 0, w = 1e-5; end
if OPF.flatstart == 1
  DAE.V = ones(Bus.n,1);
  DAE.a = zeros(Bus.n,1);
  Bus.Qg = 1e-3*ones(Bus.n,1);
else
  length(Snapshot.V);
  DAE.V = Snapshot(1).V;
  DAE.a = Snapshot(1).ang;
  Bus.Qg = Snapshot(1).Qg;
end
if ~Demand.n,
  if OPF.basepg & ~clpsat.init
    Settings.ok = 0;
    uiwait(fm_choice(['It is strongly recommended to exclude ' ...
                      'base case powers. Do you want to do so?']))
    OPF.bsaegen = ~Settings.ok;
  end
  noDem = 1;
  Demand.bus = 1;
  Demand.con = [1,100,1,zeros(1,12)];
  Demand.con(1,5) = 1e-5;
  Demand.n = 1;
else
  noDem = 0;
end

Bus.Pg = Snapshot(1).Pg;
if ~OPF.basepl
  buspl = Snapshot(1).Pl;
  busql = Snapshot(1).Ql;
  Bus.Pl(:) = 0;
  Bus.Ql(:) = 0;
  pqcon = PQ.con(:,[4 5]);
  PQ.con(:,[4 5]) = 0;
end
if ~OPF.basepg
  ploss = Snapshot(1).Ploss;
  Snapshot(1).Ploss = 0;
  buspg = Snapshot(1).Pg;
  busqg = Snapshot(1).Qg;
  Bus.Pg(:) = 0;
  Bus.Qg(:) = 0;
  pvcon = PV.con(:,4);
  PV.con(:,4) = 0;
end

% ===========================================================================
% Definition of vectors: parameters, variables and Jacobian matrices
% ===========================================================================

% Parameters
%____________________________________________________________________________

Csa = Supply.con(:,7)/100;
Csb = Supply.con(:,8);
Csc = 100*Supply.con(:,9);

Dsa = Supply.con(:,10)/100;
Dsb = Supply.con(:,11);
Dsc = 100*Supply.con(:,12);

Cda = Demand.con(:,8)/100;
Cdb = Demand.con(:,9);
Cdc = 100*Demand.con(:,10);

Dda = Demand.con(:,11)/100;
Ddb = Demand.con(:,12);
Ddc = 100*Demand.con(:,13);

Psmax = Supply.con(:,4);
Pdmax = Demand.con(:,5);
Psmin = Supply.con(:,5);
Pdmin = Demand.con(:,6);

% Tiebreaks status
if ~OPF.tiebreak,
  KTBD = zeros(Demand.n,1);
  KTBS = zeros(Supply.n,1);
else
  % Demand Tiebreaks
  if length(Demand.con(1,:)) < 15,
    KTBD = zeros(Demand.n,1);
  else,
    KTBD = Demand.con(:,15);
  end
  idx = find(Pdmax == 0);
  if ~isempty(idx)
    KTBD(idx) = 0;
  end
  idx = find(Pdmax);
  if ~isempty(idx)
    KTBD(idx) = KTBD(idx)./Pdmax(idx);
  end
  % Supply Tiebreaks
  if length(Supply.con(1,:)) < 14,
    KTBS = zeros(Supply.n,1);
  else,
    KTBS = Supply.con(:,14);
  end
  idx = find(Psmax == 0);
  if ~isempty(idx)
    KTBS(idx) = 0;
  end
  idx = find(Psmax);
  if ~isempty(idx)
    KTBS(idx) = KTBS(idx)./Psmax(idx);
  end
end

if Rsrv.n,
  Cr  = Rsrv.con(:,5);
  Prmax = Rsrv.con(:,3);
  Prmin = Rsrv.con(:,4);
end

%idx = find(Psmax == 0);
%if~isempty(Psmax)
%  Psmax(idx) = 0.00001;
%end
%idx = find(Pdmax == 0);
%if~isempty(Pdmax)
%  Pdmax(idx) = 0.00001;
%end

idx = find(Demand.con(:,3) ~= 0);
qonp = zeros(Demand.n,1);
qonp(idx) = Demand.con(idx,4)./Demand.con(idx,3);

n_gen = PV.n+SW.n;
busG = [SW.bus; PV.bus];
busS = [];
for i = 1:Supply.n,
  busS = [busS; find(busG == Supply.bus(i))];
end

busR = [];
supR = [];
for i = 1:Rsrv.n,
  busR = [busR; find(busG == Rsrv.bus(i))];
  supR = [supR; find(Supply.bus == Rsrv.bus(i))];
end

nS = 1:Supply.n;
nD = 1:Demand.n;
nG = 1:n_gen;
nB = 1:Bus.n;
nL = 1:Line.n;
nR = 1:Rsrv.n;

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*Settings.mva;
    Qgmax(idx) = 99*Settings.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*Settings.mva*ones(n_gen,1);
  Qgmax =  99*Settings.mva*ones(n_gen,1);
end

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
Vcmin = Vmin;
Vcmax = Vmax;

if OPF.enflow
  Iijmax = Line.con(:,12+OPF.flow).^2;
  idx_no_I_max = find(Iijmax == 0);
  % no limits if Iijmax is zero
  Iijmax(idx_no_I_max) = 999^2;
else
  Iijmax = 999*999.*ones(Line.n,1);
end
Iijcmax = Iijmax;

% distributed slack bus coefficients
if SW.n, ksw = SW.con(:,11); end
if PV.n, kpv = PV.con(:,10); end
ksu = zeros(Supply.n,1);
for i = 1:PV.n
  idx = find(Supply.bus == PV.bus(i));
  if idx, ksu(idx) = PV.con(i,10); end
end
for i = 1:SW.n
  idx = find(Supply.bus == SW.bus(i));
  if idx, ksu(idx) = SW.con(i,11); end
end

iteration = 0;
iter_max = Settings.lfmit;

% Graphical Settings
%____________________________________________________________________________

fm_status('opf','init',iter_max,{Theme.color08,'b','g','y'}, ...
          {'-','-','-','-'},{Theme.color11, Theme.color11, ...
                    Theme.color11,Theme.color11})

% Variables
%____________________________________________________________________________

% Lagrangian multipliers [mu]
mu_Iijmax  = zeros(Line.n,1);
mu_Ijimax  = zeros(Line.n,1);
mu_t = zeros(Bus.n,1);

% numbering
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
n_y = Supply.n+Demand.n+2*n_gen+4*Bus.n+2+Rsrv.n;
n_a = n_gen + Demand.n + Supply.n;
n2 = 2*Bus.n;
n3 = 3*Bus.n;
n4 = 4*Bus.n;
n_b = 2*n_a+2*n_gen+n4+2;
n_c = n_b+4*Line.n;
n_d = Supply.n+Demand.n+2*n_gen+4*Bus.n+2;

% y variables
Vc = DAE.V;
angc = DAE.a;
lc = lcmin + 0.001*(lcmax-lcmin);
kg = 0.0001;

ro = zeros(4*Bus.n,1);     % Dual Variables [ro]
sigma = OPF.sigma;         % Centering Parameter [sigma]
ms = sigma/n_s;            % Barrier Parameter [ms]
gamma = OPF.gamma;         % Safety Factor [gamma]
epsilon_mu = OPF.eps_mu;   % Convergence Tolerances
epsilon_1  = OPF.eps1;
epsilon_2  = OPF.eps2;
G_obj = 1;                 % Objective Function

% Jacobian Matrices
%____________________________________________________________________________

if Settings.octave
  Jz1 = zeros(n2,n2+n_gen+2+Rsrv.n);
  Jz2 = zeros(n2,n2+n_gen);
  g_Qg = zeros(n2,n_gen);
  g_Ps = zeros(n2,Supply.n);
  g_Pd = zeros(n2,Demand.n);
  g_Qg = g_Qg + sparse(Bus.n+busG,nG,-ones(n_gen,1),n2,n_gen);
  g_Ps = g_Ps + sparse(Supply.bus,nS,-ones(Supply.n,1),n2,Supply.n);
  g_Pd = g_Pd + sparse(Demand.bus,nD,ones(Demand.n,1),n2,Demand.n);
  g_Pd = g_Pd + sparse(Bus.n+Demand.bus,nD,qonp,n2,Demand.n);
  g_Qgc = zeros(n2,n_gen);
  g_Qgc = g_Qgc + sparse(Bus.n+busG,nG,-ones(n_gen,1),n2,n_gen);
  g_Psc = zeros(n2,Supply.n);
  g_Pdc = zeros(n2,Demand.n);
  g_lc  = zeros(n2,1);
  g_kg  = zeros(n2,1);
  dF_dy = zeros(n_y,1);
  dG_dy = zeros(n_y,1);
  dG_dy = dG_dy + sparse(n2+n_gen+1:n2+n_a,1,(1-w)*[Csb;-Cdb],n_y,1);
  dH_dy = zeros(n_y,1);
  gy = zeros(n_y,1);
  Jh = zeros(n_s,n_y);
else
  Jz1 = sparse(n2,n2+n_gen+2+Rsrv.n);
  Jz2 = sparse(n2,n2+n_gen);
  g_Qg = sparse(Bus.n+busG,nG,-ones(n_gen,1),n2,n_gen);
  g_Ps = sparse(Supply.bus,nS,-ones(Supply.n,1),n2,Supply.n);
  g_Pd = sparse(Demand.bus,nD,ones(Demand.n,1),n2,Demand.n);
  g_Pd = g_Pd + sparse(Bus.n+Demand.bus,nD,qonp,n2,Demand.n);
  g_Qgc = sparse(Bus.n+busG,nG,-ones(n_gen,1),n2,n_gen);
  g_Psc = sparse(n2,Supply.n);
  g_Pdc = sparse(n2,Demand.n);
  g_lc  = sparse(n2,1);
  g_kg  = sparse(n2,1);
  dF_dy = sparse(n_y,1);
  dG_dy = sparse(n2+n_gen+1:n2+n_a,1,(1-w)*[Csb;-Cdb],n_y,1);
  dH_dy = sparse(n_y,1);
  gy = sparse(n_y,1);