fm_opfm.m

电力系统的psat

function fm_opfm
%FM_OPFMR solves the OPF-based  electricity market problem by means of
%         an Interior Point Method with a Merhotra Predictor-Corrector
%         or Newton direction technique.
%
%System equations:
%
%Min:  Cs'*Ps - Cd'Pd
%
%s.t.: f(theta,V,Qg,Ps,Pd) = 0           (PF eq.)
%      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)
%      |Iji(theta,V)|   <= I_max
%      Qg_min  <= Qg  <= Qg_max          (gen. Q limits)
%      V_min  <= V  <= V_max             (bus DAE.V limits)
%
%(* optional constraints)
%
%see also FM_OPFSD FM_PARETO FM_ATC and the OPF structure for settings.
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    05-Mar-2004
%Version:   1.1.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 in order to run OPF routines',2)
  return
end

method = OPF.method;
flow = OPF.flow;

if OPF.show
  fm_disp
  fm_disp('----------------------------------------------------------------')
  fm_disp('Interior Point Method for OPF Computation')
  fm_disp('Social Benefit Objective Function')
  fm_disp(['Data file "',Path.data,File.data,'"'])
  fm_disp('----------------------------------------------------------------')
  fm_disp
end

tic;
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.basepg = ~Settings.ok;
  end
  noDem = 1;
  Demand.bus = 1;
  Demand.con = [1,100,1,zeros(1,12)];
  Demand.con(1,5) = 1e-6;
  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
DPr = [];

qonp = zeros(Demand.n,1);
idx = find(Demand.con(:,3) ~= 0);
if ~isempty(idx)
  qonp(idx) = Demand.con(idx,4)./Demand.con(idx,3);
end
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
%busS = busG(busS);
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

if OPF.enflow
  Iijmax = Line.con(:,12+flow).^2;
  idx_no_I_max = find(Iijmax == 0);
  % no limit for undefined limit
  Iijmax(idx_no_I_max) = 999^2;
else
  Iijmax = 999*999.*ones(Line.n,1);
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+2*n_gen+2*Bus.n+2*Line.n+2*Rsrv.n;
if Rsrv.n, n_s = n_s + 1; end
n_y = Supply.n+Demand.n+n_gen+2*Bus.n+Rsrv.n;
n_a = n_gen + Demand.n + Supply.n;
n2 = 2*Bus.n;
n_b = 2*n_a+n2;
n_c = n_b+2*Line.n;
n_d = Supply.n+Demand.n+n_gen+2*Bus.n;

ro = zeros(2*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
  g_Qg = zeros(n2,n_gen);
  g_Ps = zeros(n2,Supply.n);
  g_Pd = zeros(n2,Demand.n);
  dG_dy = zeros(n_y,1);
  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);
  dF_dy = zeros(n_y,1);
  dG_dy = dG_dy + sparse(n2+n_gen+1:n2+n_a,1,[Csb;-Cdb],n_y,1);
  dH_dy = zeros(n_y,1);
  gy = zeros(n_y,1);
  Jh = zeros(n_s,n_y);
else
  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);
  dF_dy = sparse(n_y,1);
  dG_dy = sparse(n2+n_gen+1:n2+n_a,1,[Csb;-Cdb],n_y,1);
  dH_dy = sparse(n_y,1);
  gy = sparse(n_y,1);
  Jh = sparse(n_s,n_y);
end

g_Pd = g_Pd + sparse(Bus.n+Demand.bus,nD,qonp,n2,Demand.n);
Jh = Jh + sparse(nS,n2+n_gen+nS,-1,n_s,n_y);
Jh = Jh + sparse(Supply.n+nS,n2+n_gen+nS,1,n_s,n_y);
Jh = Jh + sparse(2*Supply.n+nD,n2+n_gen+Supply.n+nD,-1,n_s,n_y);
Jh = Jh + sparse(2*Supply.n+Demand.n+nD,n2+n_gen+Supply.n+nD,1,n_s,n_y);
Jh = Jh + sparse(2*Supply.n+2*Demand.n+nG,n2+nG,-1,n_s,n_y);
Jh = Jh + sparse(2*Supply.n+2*Demand.n+n_gen+nG,n2+nG,1,n_s,n_y);
Jh = Jh + sparse(2*n_a+nB,Bus.n+nB,-1,n_s,n_y);
Jh = Jh + sparse(2*n_a+Bus.n+nB,Bus.n+nB,1,n_s,n_y);

if Rsrv.n, % Power Reserve
  if Settings.octave
    g_Pr  = zeros(n2,Rsrv.n);
  else
    g_Pr  = sparse(n2,Rsrv.n);
  end
  dG_dy(n_d+nR) = Cr;
  Jh = Jh + sparse(n_c+nR,n_d+nR,-1,n_s,n_y);
  Jh = Jh + sparse(n_c+Rsrv.n+nR,n_d+nR,1,n_s,n_y);
  Jh = Jh + sparse(n_c+2*Rsrv.n+1,n_d+nR,1,n_s,n_y);
  Jh = Jh + sparse(n_c+2*Rsrv.n+1,n2+n_gen+Supply.n+nD,-1,n_s,n_y);
end

Jg = sparse(n2,n_y);

% Hessian Matrices but Jlfv
%____________________________________________________________________________

if Settings.octave
  I_smu = eye(n_s);
  Z1 = zeros(n_s,n_y);
  Z2 = zeros(n_s,n_s);
  Z3 = zeros(n2,n2);
  Z4 = zeros(n2,n_s);
  H31 = zeros(n2,n_a+Rsrv.n);
else
  I_smu = speye(n_s);
  Z1 = sparse(n_s,n_y);
  Z2 = sparse(n_s,n_s);
  Z3 = sparse(n2,n2);
  Z4 = sparse(n2,n_s);
  H31 = sparse(n2,n_a+Rsrv.n);
end

% ===========================================================================
% Choosing the Initial Point
% ===========================================================================


% Primal Variables
%____________________________________________________________________________

Ps = Psmin + 0.1*(Psmax-Psmin);
if Rsrv.n, Pdbas = sum(Pr)/Supply.n; else, Pdbas = 0; end
Pd = min(Pdmin + 0.1*(Pdmax-Pdmin) + Pdbas,0.9*Pdmax);
Qg = Bus.Qg(busG);
[Iij,Jac_Iij,Hess_Iij,Iji,Jac_Iji,Hess_Iji] =fm_flows(flow,mu_Iijmax, mu_Ijimax);
a = []; b = [];
a = find(Iij > Iijmax);    b = find(Iji > Iijmax);
if ~isempty(a) | ~isempty(b)
  fm_disp('Actual line flows are greater than imposed limits',1)
  fm_disp(['Check limits of the following lines:'],1)
  fm_disp('          from     to',1)
  for i = 1:length(a)
    fm_disp(['    Line ',fvar(Line.con(a(i),1),5), ...
             ' -> ', fvar(Line.con(a(i),2),5),': ',fvar(sqrt(Iij(a(i))),7), ...
             ' > ', fvar(Line.con(a(i),12+flow),7)],1)
  end
  for i = 1:length(b)
    fm_disp(['    Line ',fvar(Line.con(b(i),2),5), ...
             ' -> ', fvar(Line.con(b(i),1),5),': ',fvar(sqrt(Iji(b(i))),7), ...
             ' > ', fvar(Line.con(b(i),12+flow),7)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
a = find(Qg > Qgmax);
b = find(Qg < Qgmin);
if ~isempty(a),
  fm_disp(['Max Reactive Power limit not respected at buses:'],1)
  for i = 1:length(a)
    fm_disp(['     Bus #',fvar(Bus.con(busG(a(i)),1),4),' -> ',fvar(Qg(a(i)),8), ...
             ' > ', fvar(Qgmax(a(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
if ~isempty(b),
  fm_disp(['Min Reactive Power limit not respected at buses: ',num2str(b')],1)
  for i = 1:length(b)
    fm_disp(['     Bus #',fvar(Bus.con(busG(b(i)),1),4),' -> ',fvar(Qg(b(i)),8), ...
             ' < ', fvar(Qgmin(b(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
a = find(DAE.V > Vmax);
b = find(DAE.V < Vmin);
if ~isempty(a),
  fm_disp(['Max Voltage limit not respected at buses: ',num2str(a')],1)
  for i = 1:length(a)
    fm_disp(['     Bus #',fvar(Bus.con(a(i),1),4),' -> ',fvar(DAE.V(a(i)),8), ...
             ' > ', fvar(Vmax(a(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
if ~isempty(b),
  fm_disp(['Min Voltage limit not respected at buses: ',num2str(b')],1)
  for i = 1:length(b)
    fm_disp(['     Bus #',fvar(Bus.con(b(i),1),4),' -> ',fvar(DAE.V(b(i)),8), ...
             ' < ', fvar(Vmin(b(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end

h_delta_Ps   = Psmax - Psmin;
h_delta_Pd   = Pdmax - Pdmin;
h_delta_Qg   = Qgmax - Qgmin;
h_delta_V    = Vmax  - Vmin;
h_delta_Iij  = Iijmax;

gamma_h = 0.25;

a_Ps  = min(max(gamma_h*h_delta_Ps,Ps-Psmin),(1-gamma_h)*h_delta_Ps);
a_Pd  = min(max(gamma_h*h_delta_Pd,Pd-Pdmin),(1-gamma_h)*h_delta_Pd);