fm_opfsdr.m

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

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-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 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.basepg = ~Settings.ok;
  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
  Bus.Pg(:) = 0;
  Bus.Qg(:) = 0;
  PV = pvzero(PV,'all');
  SW = swzero(SW,'all');
end

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

% Parameters
% ===========================================================================

% Supply parameters
[Csa,Csb,Csc,Dsa,Dsb,Dsc] = costs(Supply);
[Psmax,Psmin] = plim(Supply);
KTBS = tiebreaks(Supply);
ksu = getgamma(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

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

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

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);
DAE.Gl = zeros(n2,1);
DAE.Gk = zeros(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);
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);
Jh = Jh + sparse(2*n_a+n2+nG,1+n4+n_a+nG,-1,n_s,n_y);
Jh = Jh + sparse(2*n_a+n2+n_gen+nG,1+n4+n_a+nG,1,n_s,n_y);
Jh = Jh + sparse(2*n_a+2*n_gen+n2+nB,n3+n_a+nB,-1,n_s,n_y);
Jh = Jh + sparse(2*n_a+2*n_gen+n3+nB,n3+n_a+nB,1,n_s,n_y);
Jh(2*n_a+2*n_gen+n4+1, n_d) = -1;
Jh(2*n_a+2*n_gen+n4+2, n_d) =  1;

if Rsrv.n, % Power Reserve
  g_Pr  = sparse(n2,Rsrv.n);
  dG_dy(n_d+nR) = (1-w)*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(n4,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(n4,n4);
Z4 = sparse(n4,n_s);
Hcolk = sparse(n2,1);
Hrowk = sparse(1,n2+1);
H31 = sparse(n2,n_gen+n_a+n2+2+Rsrv.n);
H41 = sparse(n2+1,n2+n_a);
H42 = sparse(n2+1,n_gen+1+Rsrv.n);
H5  = sparse(n_gen+1+Rsrv.n,n_y);

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

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

Ps = Psmin + 0.5*(Psmax-Psmin);
if Rsrv.n
  Pdbas = sum(Pr)/Supply.n;
else
  Pdbas = 0;
end
Pd = Pdmin + 0.5*(Pdmax-Pdmin) + Pdbas;
Qg = 0.5*(Qgmax+Qgmin);
Qgc = Qg;

[Iij,Jac_Iij,Hess_Iij,Iji,Jac_Iji,Hess_Iji] = fm_flows(OPF.flow,mu_Iijmax, mu_Ijimax);
Iijc = Iij;
Ijic = Iji;
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+OPF.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+OPF.flow),7)],1)
  end
end
if ~isempty(a)
  a1 = Line.to(a);
  a2 = Line.from(a);
  DAE.V(a1) = DAE.V(a2);
  DAE.a(a1) = DAE.a(a2);
  Vc(a1) = DAE.V(a2);
  angc(a1) = DAE.a(a2);
end
if ~isempty(b)
  b1 = Line.to(b);
  b2 = Line.from(b);
  DAE.V(b1) = DAE.V(b2);
  DAE.a(b1) = DAE.a(b2);
  Vc(b1) = DAE.V(b2);
  angc(b1) = DAE.a(b2);
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
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
end
if ~isempty(a) | ~isempty(b)
  DAE.V = max(DAE.V,Vmin+1e-3);
  DAE.V = min(DAE.V,Vmax-1e-3);
  Vc = DAE.V;
end

h_delta_Ps   = Psmax - Psmin;
h_delta_Pd   = Pdmax - Pdmin;
h_delta_Qg   = Qgmax - Qgmin;
h_delta_V    = Vmax  - Vmin;
h_delta_Vc   = Vcmax - Vcmin;
h_delta_lc   = lcmax - lcmin;
h_delta_Iij  = Iijmax;
h_delta_Iji  = Iijmax;
h_delta_Iijc = Iijcmax;
h_delta_Ijic = Iijcmax;

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_Vc  = min(max(gamma_h*h_delta_Vc,Vc-Vcmin),(1-gamma_h)*h_delta_Vc);
a_lc  = min(max(gamma_h*h_delta_lc,lc-lcmin),(1-gamma_h)*h_delta_lc);
a_Iij = min(max(gamma_h*h_delta_Iij,Iij),(1-gamma_h)*h_delta_Iij);
a_Iji = min(max(gamma_h*h_delta_Iji,Iji),(1-gamma_h)*h_delta_Iji);
a_Iijc= min(max(gamma_h*h_delta_Iijc,Iij),(1-gamma_h)*h_delta_Iijc);
a_Ijic= min(max(gamma_h*h_delta_Ijic,Iji),(1-gamma_h)*h_delta_Ijic);

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; ...
     a_Qg; h_delta_Qg - a_Qg; a_Vc; h_delta_Vc - a_Vc; ...
     a_lc; h_delta_lc - a_lc; h_delta_Iij  - a_Iij; ...
     h_delta_Iijc - a_Iijc; h_delta_Iji  - a_Iji; h_delta_Ijic - a_Ijic];

idx = find(s == 0);
if ~isempty(idx),
  s(idx) = 1e-6;
  fm_disp('Problems in initializing slack variables...')
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;
if w < 1, ro([1:Bus.n]) = -mean([Csb; Cdb]); end

% Contingency on selected line
% =====================================================


if OPF.line
  line_orig = Line.con;
  Y_orig = Line.Y;
  Line.con(OPF.line,[8 9 10]) = [0, 1e5, 0];
  line_cont = Line.con;
  fm_y;
  Y_cont = Line.Y;
  Line.con = line_orig;
end

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

  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;
  mu_Qgmin  = mu(idx+nG);     idx = idx + n_gen;
  mu_Qgmax  = mu(idx+nG);     idx = idx + n_gen;
  mu_Vmin   = mu(idx+nB);     idx = idx + Bus.n;
  mu_Vmax   = mu(idx+nB);     idx = idx + Bus.n;
  mu_Qgcmin = mu(idx+nG);     idx = idx + n_gen;
  mu_Qgcmax = mu(idx+nG);     idx = idx + n_gen;
  mu_Vcmin  = mu(idx+nB);     idx = idx + Bus.n;
  mu_Vcmax  = mu(idx+nB);     idx = idx + Bus.n;
  mu_lcmin  = mu(idx+1);      idx = idx + 1;
  mu_lcmax  = mu(idx+1);      idx = idx + 1;
  mu_Iijmax = mu(idx+nL);     idx = idx + Line.n;