fm_opfm.m

电力系统分析计算程序

function fm_opfm
%FM_OPFM 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 voltage limits)
%
%(* optional constraints)
%
%see also FM_OPFSDR 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://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

fm_var

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

if DAE.n
  fm_disp('OPF routine does not support dynamic components',2)
  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;
forcepq = Settings.forcepq;
Settings.forcepq = 1;

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.y(Bus.a) = getzeros(Bus);
  DAE.y(Bus.v) = getones(Bus);
  Bus.Qg = 1e-3*getones(Bus);
else
  length(Snapshot.y);
  DAE.y = Snapshot(1).y;
  Bus.Qg = Snapshot(1).Qg;
end

if ~Demand.n,
  if OPF.basepg & ~clpsat.init
    Settings.ok = 0;
    uiwait(fm_choice(['It is recommended to exclude ' ...
                      'base case powers. Do you want to do so?']))
    OPF.basepg = ~Settings.ok;
    if Fig.opf
      hdl = findobj(Fig.opf,'Tag','CheckboxBaseGen');
      set(hdl,'Value',OPF.basepg)
    end
  end
  noDem = 1;
  Demand = add(Demand,'dummy');
else
  noDem = 0;
end

Bus.Pg = Snapshot(1).Pg;
if ~OPF.basepl
  Bus.Pl(:) = 0;
  Bus.Ql(:) = 0;
  PQ = pqzero(PQ,'all');
end
if ~OPF.basepg
  ploss = Snapshot(1).Ploss;
  Snapshot(1).Ploss = 0;
  Bus.Pg(:) = 0;
  Bus.Qg(:) = 0;
  SW = swzero(SW,'all');
  PV = pvzero(PV,'all');
end

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

% Supply parameters
[Csa,Csb,Csc,Dsa,Dsb,Dsc] = costs(Supply);
[Psmax,Psmin] = plim(Supply);
KTBS = tiebreaks(Supply);

% Demand parameters
[Cda,Cdb,Cdc,Dda,Ddb,Ddc] = costs(Demand);
[Pdmax,Pdmin] = plim(Demand);
KTBD = tiebreaks(Demand);
qonp = tanphi(Demand);

% Reserve parameters
Cr = costs(Rsrv);
[Prmax,Prmin] = plim(Rsrv);
DPr = [];
Pr = [];

busG = double([SW.bus;PV.bus]);
busS = setdiff(busG,Supply.bus);
n_gen = Supply.n+length(busS);
busG = [Supply.bus;busS];
[supR,idxR,dummy] = intersect(Supply.bus,Rsrv.bus);

nS = 1:Supply.n;
nD = 1:Demand.n;
nG = 1:n_gen;
nB = Bus.a';
nV = Bus.v';
nY = 1:DAE.m;
nL = 1:Line.n;
nR = 1:Rsrv.n;

if OPF.enreac
  [Qgmax,Qgmin] = fm_qlim('gen');
else
  Qgmin = -99*Settings.mva*ones(n_gen,1);
  Qgmax =  99*Settings.mva*ones(n_gen,1);
end

if OPF.envolt
  [Vmax,Vmin] = fm_vlim(OPF.vmax,OPF.vmin);
else
  Vmin = OPF.vmin*getones(Bus);
  Vmax = OPF.vmax*getones(Bus);
end

if OPF.enflow
  Iijmax = getflowmax(Line,flow).^2;
else
  Iijmax = 1e6*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_t = getzeros(Bus);

% numbering
n1 = Bus.n;
n2 = 2*n1;
n_a = n_gen + Demand.n + Supply.n;
n_b = 2*n_a+n2;
n_c = n_b+2*Line.n;
n_d = Supply.n+Demand.n+n_gen+n2;
n_s = 2*Supply.n+2*Demand.n+2*n_gen+n2+2*Line.n+2*Rsrv.n;
if Rsrv.n, n_s = n_s + 1; end
n_y = Supply.n+Demand.n+n_gen+n2+Rsrv.n;

ro = zeros(n2,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
% ===========================================================================

g_Qg = sparse(n1+busG,nG,-1,n2,n_gen);
g_Ps = sparse(Supply.bus,nS,-1,n2,Supply.n);
g_Pd = sparse(Demand.bus,nD, 1,n2,Demand.n);
g_Pd = g_Pd + sparse(Demand.vbus,nD,qonp,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);
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,nV,1,n_s,n_y);
Jh = Jh + sparse(2*n_a+nV,nV,1,n_s,n_y);

if Rsrv.n, % Power Reserve
  g_Pr  = sparse(n2,Rsrv.n);
  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 Gy
% ===========================================================================

Z3 = sparse(n2,n2);
H31 = sparse(n2,n_a+Rsrv.n);

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


% Primal Variables
% ===========================================================================

Ps = Psmin + 0.1*(Psmax-Psmin);
if Rsrv.n
  Pr = Prmin + 0.1*(Prmax-Prmin);
  Pdbas = sum(Pr)/Supply.n;
else
  Pdbas = 0;
end
Pd = min(Pdmin + 0.1*(Pdmax-Pdmin) + Pdbas,0.9*Pdmax);
Qg = 0.5*(Qgmax+Qgmin);
[Iij,Iji] = fjh2(Line,flow);
a = []; b = [];

% check flow limits
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.fr(a(i)),5), ...
             ' -> ', fvar(Line.to(a(i)),5),': ',fvar(sqrt(Iij(a(i))),7), ...
             ' > ', fvar(sqrt(Iijmax(a(i))),7)],1)
  end
  for i = 1:length(b)
    fm_disp(['    Line ',fvar(Line.to(b(i)),5), ...
             ' -> ', fvar(Line.fr(b(i)),5),': ',fvar(sqrt(Iji(b(i))),7), ...
             ' > ', fvar(sqrt(Iijmax(a(i))),7)],1)
  end
end
if ~isempty(b)
  DAE.y(Line.vto(b)) = DAE.y(Line.vfr(b));
  DAE.y(Line.to(b))  = DAE.y(Line.fr(b));
end

% check voltage limits
Vbus = DAE.y(Bus.v);
a = find(Vbus > Vmax);
b = find(Vbus < 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(getidx(Bus,a(i)),4), ...
             ' -> ',fvar(Vbus(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(getidx(Bus,b(i)),4), ...
             ' -> ',fvar(Vbus(b(i)),8), ...
             ' < ', fvar(Vmin(b(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
if ~isempty(a) | ~isempty(b)
  Vbus = max(Vbus,Vmin+1e-3);
  Vbus = min(Vbus,Vmax-1e-3);
  DAE.y(Bus.v) = Vbus;
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);
a_Qg  = min(max(gamma_h*h_delta_Qg,Qg-Qgmin),(1-gamma_h)*h_delta_Qg);
a_V   = min(max(gamma_h*h_delta_V,Vbus-Vmin),(1-gamma_h)*h_delta_V);
a_Iij = min(max(gamma_h*h_delta_Iij,Iij),(1-gamma_h)*h_delta_Iij);
a_Iji = min(max(gamma_h*h_delta_Iij,Iji),(1-gamma_h)*h_delta_Iij);

s = [a_Ps; h_delta_Ps - a_Ps; a_Pd; h_delta_Pd - a_Pd; ...
     a_Qg; h_delta_Qg - a_Qg; a_V; h_delta_V  - a_V; ...
     h_delta_Iij  - a_Iij; h_delta_Iij  - a_Iji];

idx = find(s == 0);
if ~isempty(idx),
  s(idx) = 1e-6;
end

if Rsrv.n
  %Pr = Prmin + 0.1*(Prmax-Prmin);
  sumPrd = sum(Pr)-sum(Pd);
  h_delta_Pr   = Prmax - Prmin;
  a_Pr  = min(max(gamma_h*h_delta_Pr,Pr),(1-gamma_h)*h_delta_Pr);
  s = [s; a_Pr; h_delta_Pr - a_Pr; -sumPrd];
end

% Dual Variables
%____________________________________________________________________________

mu = ms./s;
ro(Bus.a) = -mean([Csb; Cdb]);
if Settings.matlab, warning('off'); end

% =======================================================================
% PRIMAL DUAL INTERIOR-POINT METHOD
% =======================================================================

while 1

  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  % =======================================================================
  % KKT optimality necessary condition: [DyLms] = 0
  % =======================================================================

  % Saving Objective Function (k-1)
  %________________________________________________________________________

  G_obj_k_1 = G_obj;

  mu_Psmin  = mu(nS);         idx = Supply.n;
  mu_Psmax  = mu(idx+nS);     idx = idx + Supply.n;
  mu_Pdmin  = mu(idx+nD);     idx = idx + Demand.n;
  mu_Pdmax  = mu(idx+nD);     idx = idx + Demand.n;