fm_opfm.m

这是一个很适合研究和学习用的电力系统仿真软件

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-2006 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;

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 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;
PQ = novlim(PQ,'all');
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);

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

busG = double([SW.bus; PV.bus]);
[busS,idxS] = setdiff(busG,Supply.bus);
n_gen = Supply.n+length(busS);
busG = [Supply.bus;busS];
supR = intersect(Supply.bus,Rsrv.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;

if OPF.enreac
  [Qgmax,Qgmin] = fm_qlim(99,-99,'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*ones(Bus.n,1);
  Vmax = OPF.vmax*ones(Bus.n,1);
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
% ===========================================================================

g_Qg = sparse(Bus.n+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(Bus.n+Demand.bus,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,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
  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 Jlfv
% ===========================================================================

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);

% ===========================================================================
% 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 = 0.5*(Qgmax+Qgmin);
[Iij,Jac_Iij,Hess_Iij,Iji,Jac_Iji,Hess_Iji] =fm_flows(flow,mu_Iijmax, mu_Ijimax);
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.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
end
if ~isempty(b)
  DAE.V(Line.to(b)) = DAE.V(Line.from(b));
  DAE.a(Line.to(b)) = DAE.a(Line.from(b));
end

% check voltage limits
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
if ~isempty(a) | ~isempty(b)
  DAE.V = max(DAE.V,Vmin+1e-3);
  DAE.V = min(DAE.V,Vmax-1e-3);
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,DAE.V-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([1:Bus.n]) = -mean([Csb; Cdb]);
warning('off');

% =======================================================================
% 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)
  %________________________________________________________________________