fm_opfsdr.m

电力系统分析计算程序

  if Rsrv.n
    mu_Prmin  = mu(idx+nR);   idx = idx + Rsrv.n;
    mu_Prmax  = mu(idx+nR);   idx = idx + Rsrv.n;
    mu_sumPrd = mu(idx+1);
  end

  % Computations for the System:  f(thetac,Vc,Qgc,Ps,Pd,lc) = 0
  % =====================================================

  y_snap = DAE.y;
  DAE.y  = yc;
  DAE.Gl = zeros(n2,1);
  DAE.Gk = zeros(n2,1);

  if OPF.line, Line = setstatus(Line,OPF.line,0); end

  Line = gcall(Line);
  gcall(Shunt)
  glambda(SW,1+lc,kg);
  glambda(PV,1+lc,kg);
  glambda(PQ,1+lc);

  Glcall(SW);
  Gkcall(SW);
  Glcall(PV);
  Gkcall(PV);
  Glcall(PQ);

  % generator reactive powers
  DAE.g = DAE.g - sparse(busG+n1,1,Qgc,DAE.m,1);

  % Demand & Supply
  DAE.g = DAE.g + sparse(Demand.bus,1,(1+lc)*Pd,DAE.m,1) ...
          + sparse(Demand.vbus,1,(1+lc)*Pd.*qonp,DAE.m,1);
  DAE.Gl = DAE.Gl + sparse(Demand.bus,1,Pd,DAE.m,1);
  DAE.Gl = DAE.Gl + sparse(Demand.vbus,1,Pd.*qonp,DAE.m,1);
  g_Pdc = sparse(Demand.bus,nD,1+lc,DAE.m,Demand.n);
  g_Pdc = sparse(Demand.vbus,nD,(1+lc)*qonp,DAE.m,Demand.n);
  DAE.g = DAE.g - sparse(Supply.bus,1,(1+lc+kg*ksu).*Ps,DAE.m,1);
  g_Psc = sparse(Supply.bus,nS,-(1+lc+kg*ksu),DAE.m,Supply.n);
  DAE.Gl = DAE.Gl - sparse(Supply.bus,1,Ps,DAE.m,1);
  DAE.Gk = DAE.Gk - sparse(Supply.bus,1,ksu.*Ps,DAE.m,1);

  gc1 = DAE.g;
  gc2p = Line.p;
  gc2q = Line.q;

  Gycall(Line)
  Gycall(Shunt)
  Gyc = DAE.Gy;
  [Iijc,Jijc,Hijc,Ijic,Jjic,Hjic] = fjh2(Line,OPF.flow,mu_Iijcmax,mu_Ijicmax);
  Hess_c = -hessian(Line,ro(n2+1:n4))+Hijc+Hjic-hessian(Shunt,ro(n2+1:n4));

  % Computations for the System:  f(theta,V,Qg,Ps,Pd) = 0
  % =====================================================

  DAE.y = y_snap;

  if OPF.line, Line = setstatus(Line,OPF.line,status); end

  Line = gcall(Line);
  gcall(Shunt)
  gcall(PQ);
  glambda(SW,1,0);
  glambda(PV,1,0);

  % Demand & Supply
  DAE.g = DAE.g + sparse(Demand.bus,1,Pd,DAE.m,1) ...
          + sparse(Demand.vbus,1,Pd.*qonp,DAE.m,1) ...
          - sparse(Supply.bus,1,Ps,DAE.m,1) ...
          - sparse(busG+n1,1,Qg,DAE.m,1);

  Gycall(Line)
  Gycall(Shunt)
  [Iij,Jij,Hij,Iji,Jji,Hji] = fjh2(Line,OPF.flow,mu_Iijmax,mu_Ijimax);
  Hess = -hessian(Line,ro(1:n2))+Hij+Hji-hessian(Shunt,ro(1:n2));

  % Gradient of [s] variables
  % =====================================================

  gs = s.*mu - ms;

  % Gradient of [mu] variables
  % =====================================================

  V = DAE.y(Bus.v);
  Vc = yc(Bus.v);
  gmu = [Psmin-Ps;Ps-Psmax;Pdmin-Pd;Pd-Pdmax;Qgmin-Qg;Qg-Qgmax; ...
         Vmin-V;V-Vmax;Qgmin-Qgc;Qgc-Qgmax;Vcmin-Vc;Vc-Vcmax; ...
         lcmin-lc;lc-lcmax;Iij-Iijmax;Iijc-Iijcmax;Iji-Iijmax; ...
         Ijic-Iijcmax];
  if Rsrv.n
    gmu(Supply.n+supR) = gmu(Supply.n+supR) + Pr;
    gmu = [gmu; Prmin-Pr;Pr-Prmax;sum(Pr)-sum(Pd)];
  end
  gmu = gmu + s;

  % Gradient of [y] = [theta; V; Qg; Ps; Pd] variables
  % =====================================================

  if Rsrv.n,
    Jg = [DAE.Gy,g_Qg,g_Ps,g_Pd,Jz1; ...
          Jz2,g_Psc,g_Pdc,Gyc,DAE.Gk,g_Qgc,DAE.Gl,g_Pr];
  else,
    Jg = [DAE.Gy,g_Qg,g_Ps,g_Pd,Jz1; ...
          Jz2,g_Psc,g_Pdc,Gyc,DAE.Gk,g_Qgc,DAE.Gl];
  end

  dF_dy = Jg'*ro;

  dG_dy(n_d) = -w;                      % max loading factor
  dG_dy(n2+n_gen+1:n2+n_gen+Supply.n) = (1-w)*(Csb + 2*Csc.*Ps + 2*KTBS.*Ps);
  dG_dy(n2+n_gen+Supply.n+1:n2+n_a) =  -(1-w)*(Cdb + 2*Cdc.*Pd + 2*KTBD.*Pd ...
                                               + qonp.*(Ddb + 2*qonp.*Ddc.*Pd));
  dG_dy(n2+nS) = (1-w)*(Dsb + 2*Dsc.*Qg(nS));

  dH_dtV = Jij'*mu_Iijmax + Jji'*mu_Ijimax + [mu_t; mu_Vmax - mu_Vmin];
  dH_dtVc = Jijc'*mu_Iijcmax + Jjic'*mu_Ijicmax + [mu_t; mu_Vcmax - mu_Vcmin];
  if Rsrv.n,
    dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; ...
             mu_Pdmax-mu_Pdmin-mu_sumPrd; ...
             dH_dtVc; 0; mu_Qgcmax - mu_Qgcmin; mu_lcmax - mu_lcmin; ...
             mu_Psmax(idxR)+mu_Prmax-mu_Prmin+mu_sumPrd];
  else
    dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; mu_Pdmax-mu_Pdmin; ...
             dH_dtVc; 0; mu_Qgcmax - mu_Qgcmin; mu_lcmax - mu_lcmin];
  end

  gy = dG_dy - dF_dy + dH_dy;

  Jh(n_b+1:n_b+Line.n,1:n2) = Jij;
  Jh(n_b+1+Line.n:n_b+2*Line.n,1+n2+n_a:n4+n_a) = Jijc;
  Jh(n_b+1+2*Line.n:n_b+3*Line.n,1:n2) = Jji;
  Jh(n_b+1+3*Line.n:n_b+4*Line.n,1+n2+n_a:n4+n_a) = Jjic;

  % Hessian Matrix [D2xLms]
  % --------------------------------------------------------------------

  H3  = sparse(n_a,n_y);
  Hx = -ro(n2+Supply.bus);
  Sidx = 1:Supply.n;
  H3(n_gen+Sidx,n_d) = Hx;
  H3(n_gen+Sidx,n4+1+n_a) = Hx.*ksu;
  H5(n_gen+1,n2+n_gen+Sidx) = Hx';
  H41(end,n2+n_gen+Sidx) = Hx'.*ksu';
  Didx = 1:Demand.n;
  Hx = ro(n2+Demand.bus) + qonp.*ro(n3+Demand.bus);
  H3(n_gen+Supply.n+Didx,n_d) = Hx;
  H5(n_gen+1,n2+n_gen+Supply.n+Didx) = Hx';
  H3 = H3 - sparse(n_gen+nS,n2+n_gen+nS,(1-w)*(2*Csc+2*KTBS),n_a,n_y);
  H3 = H3 - sparse(nS,n2+nS,(1-w)*2*Dsc,n_a,n_y);
  H3 = H3 + sparse(n_gen+Supply.n+nD,n_gen+Supply.n+n2+nD, ...
		   (1-w)*(2*Cdc+2*KTBD+2*Ddc.*qonp.*qonp),n_a,n_y);

  D2xLms = [Hess, H31; -H3; -H41, [Hess_c, Hcolk; Hrowk], H42; -H5];

  switch OPF.method
   case 1 % Newton Directions

    % reduced system
    H_m = sparse(1:n_s,1:n_s,mu./s,n_s,n_s);
    H_s = sparse(1:n_s,1:n_s,1./s,n_s,n_s);
    Jh(:,SW.refbus+n2+n_a) = 0;
    Jh(:,SW.refbus) = 0;
    gy = gy+(Jh.')*(H_m*gmu-H_s*gs);
    Jd = [D2xLms+(Jh.')*(H_m*Jh),-Jg.';-Jg,Z3];
    % reference angle for the actual system
    Jd(SW.refbus,:) = 0;
    Jd(:,SW.refbus) = 0;
    Jd(SW.refbus,SW.refbus) = speye(length(SW.refbus));
    gy(SW.refbus) = 0;
    % reference angle for the critical system
    Jd(:,SW.refbus+n2+n_a) = 0;
    Jd(SW.refbus+n2+n_a,:) = 0;
    Jd(SW.refbus+n2+n_a,SW.refbus+n2+n_a) = speye(length(SW.refbus));
    gy(SW.refbus+n2+n_a) = 0;
    % variable increments
    Dx = -Jd\[gy; -DAE.g; -gc1];
    Ds = -(gmu+Jh*Dx([1:n_y]));
    Dm = -H_s*gs-H_m*Ds;

   case 2 % Mehrotra's Predictor-Corrector

    % -------------------
    % Predictor step
    % -------------------
    % reduced system
    H_m  = sparse(1:n_s,1:n_s,mu./s,n_s,n_s);
    Jh(:,SW.refbus+n2+n_a) = 0;
    Jh(:,SW.refbus) = 0;
    gx = gy+(Jh.')*(H_m*gmu-mu);
    Jd = [D2xLms+(Jh.')*(H_m*Jh),-Jg.';-Jg,Z3];
    % reference angle for the actual system
    Jd(SW.refbus,:) = 0;
    Jd(:,SW.refbus) = 0;
    Jd(SW.refbus,SW.refbus) = speye(length(SW.refbus));
    gx(SW.refbus) = 0;
    % reference angle for the critical system
    Jd(:,SW.refbus+n2+n_a) = 0;
    Jd(SW.refbus+n2+n_a,:) = 0;
    Jd(SW.refbus+n2+n_a,SW.refbus+n2+n_a) = speye(length(SW.refbus));
    gx(SW.refbus+n2+n_a) = 0;
    % LU factorization
    [L,U,P] = lu(Jd);
    % variable increments
    Dx = -U\(L\(P*[gx; -DAE.g; -gc1]));
    Ds = -(gmu+Jh*Dx([1:n_y]));
    Dm = -mu-H_m*Ds;
    % centering correction
    a1 = find(Ds < 0);
    a2 = find(Dm < 0);
    if isempty(a1), ratio1 = 1; else, ratio1 = -s(a1)./Ds(a1);   end
    if isempty(a2), ratio2 = 1; else, ratio2 = -mu(a2)./Dm(a2); end
    alpha_P = min(1,gamma*min(ratio1));
    alpha_D = min(1,gamma*min(ratio2));
    c_gap_af = [s + alpha_P*Ds]'*[mu + alpha_D*Dm];
    c_gap = s'*mu;
    ms = min((c_gap_af/c_gap)^2,0.2)*c_gap_af/n_s;
    gs = mu+(Ds.*Dm-ms)./s;

    % -------------------
    % Corrector Step
    % -------------------
    % new increment for variable y
    gx = gy+(Jh.')*(H_m*gmu-gs);
    gx(SW.refbus) = 0;
    gx(SW.refbus+n2+n_a) = 0;
    % variable increments
    Dx = -U\(L\(P*[gx; -DAE.g; -gc1]));
    Ds = -(gmu+Jh*Dx([1:n_y]));
    Dm = -gs-H_m*Ds;
  end

  % =======================================================================
  % Variable Increments
  % =======================================================================

  Dy  = Dx(nY);          idx = DAE.m;        % curr. sys.
  DQg = Dx(idx+nG);      idx = idx + n_gen;
  DPs = Dx(idx+nS);      idx = idx + Supply.n;
  DPd = Dx(idx+nD);      idx = idx + Demand.n;

  Dyc  = Dx(idx+nY);     idx = idx + DAE.m;  % crit. sys.
  Dkg  = Dx(1+idx);      idx = idx + 1;
  DQgc = Dx(idx+nG);     idx = idx + n_gen;
  Dlc  = Dx(1+idx);      idx = idx + 1;
  if Rsrv.n, DPr = Dx(idx+nR); idx = idx + Rsrv.n;    end

  Dro     = Dx(1+idx:end);                       % Lag. mult.

  % =======================================================================
  % Updating the Variables
  % =======================================================================

  % Step Lenght Parameters [alpha_P & alpha_D]
  %________________________________________________________________________

  a1 = find(Ds  < 0);
  a2 = find(Dm < 0);
  if isempty(a1), ratio1 = 1; else, ratio1 = (-s(a1)./Ds(a1));   end
  if isempty(a2), ratio2 = 1; else, ratio2 = (-mu(a2)./Dm(a2)); end
  alpha_P = min(1,gamma*min(ratio1));
  alpha_D = min(1,gamma*min(ratio2));

  % New primal variables
  %________________________________________________________________________

  DAE.y = DAE.y + alpha_P * Dy;
  Ps = Ps + alpha_P * DPs;
  Pd = Pd + alpha_P * DPd;
  Qg = Qg + alpha_P * DQg;
  if Rsrv.n, Pr = Pr + alpha_P * DPr; end

  yc  = yc  + alpha_P * Dyc;
  kg  = kg  + alpha_P * Dkg;
  Qgc = Qgc + alpha_P * DQgc;
  lc  = lc  + alpha_P * Dlc;
  s   = s   + alpha_P * Ds;

  % New dual variables
  %________________________________________________________________________

  ro = ro + alpha_D*Dro;
  mu = mu + alpha_D*Dm;

  % Objective Function
  %________________________________________________________________________

  s(find(s == 0)) = epsilon_mu;
  Fixd_c = sum(Csa) - sum(Cda) + sum(Dsa) - sum(Dda);
  Prop_c = Csb'*Ps  - Cdb'*Pd + Dsb'*Qg(nS) - Ddb'*(qonp.*Pd);
  TieBreaking = (sum(KTBS.*Ps.*Ps) - sum(KTBD.*Pd.*Pd));
  Quad_c = Csc'*(Ps.*Ps) - Cdc'*(Pd.*Pd) - Ddc'*(qonp.*qonp.*Pd.*Pd);
  Quad_q = Dsc'*(Qg(nS).*Qg(nS));
  if Rsrv.n, Reserve = Cr'*Pr; else, Reserve = 0; end
  G_obj = (1-w)*(Fixd_c + Prop_c + Quad_c + Quad_q + TieBreaking + Reserve) - ...
          ms*sum(log(s)) - w*lc;

  % =======================================================================
  % Reducing the Barrier Parameter
  % =======================================================================

  sigma = max(0.99*sigma, 0.1);     % Centering Parameter
  c_gap = s'*mu;                    % Complementarity Gap
  ms = min(abs(sigma*c_gap/n_s),1); % Evaluation of the Barrier Parameter

  % =======================================================================
  % Testing for Convergence
  % =======================================================================

  test1  = ms <= epsilon_mu;
  norma2 = norm(Dx,inf);
  test2  = norma2 <= epsilon_2;
  norma3 = norm([DAE.g; gc1],inf);
  test3  = norma3 <= epsilon_1;
  norma4 = abs(G_obj-G_obj_k_1)/(1+abs(G_obj));
  test4  = norma4 <= epsilon_2;
  if test1 & test2 & test3 & test4, break, end

  % Displaying Convergence Tests
  %________________________________________________________________________

  iteration = iteration + 1;
  if OPF.show
    fm_disp(['Iter. =',fvar(iteration,5),'  mu =', fvar(ms,8), ...
             '  |dy| =', fvar(norma2,8), '  |f(y)| =', ...
             fvar(norma3,8),'  |dG(y)| =' fvar(norma4,8)])
  end
  fm_status('opf','update',[iteration, ms, norma2, norma3, norma4], ...
            iteration)

  if iteration > iter_max, break, end
end

% Some settings ...
%____________________________________________________________________________

if Settings.matlab, warning('on'); end

Demand = pset(Demand,Pd);
Supply = pset(Supply,Ps);
Rsrv = pset(Rsrv,Pr);

Iij  = sqrt(Iij);
Iji  = sqrt(Iji);
Iijc = sqrt(Iijc);
Ijic = sqrt(Ijic);
Iijmax  = sqrt(Iijmax);
Iijcmax = sqrt(Iijcmax);
MVA = Settings.mva;
Pay = ro(Bus.a).*Line.p*MVA;
ISOPay = sum(Pay);

% Nodal Congestion Prices (NCPs)
%____________________________________________________________________________

fm_setgy(SW.refbus)
dH_dtV(SW.refbus,:) = 0;
OPF.NCP = -DAE.Gy'\dH_dtV;

OPF.obj = G_obj;
OPF.ms = ms;
OPF.dy = norma2;
OPF.dF = norma3;
OPF.dG = norma4;
OPF.iter = iteration;
OPF.gpc = gc2p;
OPF.gqc = gc2q;

SNB.init = 0;
LIB.init = 0;
CPF.init = 0;
OPF.init = 2;

% set Pg, Qg, Pl and Ql
Bus.Pl = OPF.basepl*Snapshot(1).Pl + sparse(Demand.bus,1,Pd,n1,1);
Bus.Ql = OPF.basepl*Snapshot(1).Ql + sparse(Demand.bus,1,Pd.*qonp,n1,1);
Bus.Pg = OPF.basepg*Snapshot(1).Pg + sparse(Supply.bus,1,Ps,n1,1);
Bus.Qg = OPF.basepg*Snapshot(1).Qg + sparse(busG,1,Qg,n1,1);

% Display Results
% --------------------------------------------------

TPQ = totp(PQ);

if (Settings.showlf | OPF.show) & clpsat.showopf

  OPF.report = cell(1,1);
  OPF.report{1,1} = ['Weighting Factor = ',fvar(w,8)];
  OPF.report{2,1} = ['Lambda = ',fvar(lc,8)];
  OPF.report{3,1} = ['Kg = ',fvar(kg,8)];
  OPF.report{4,1} = ['Total Losses = ',fvar(sum(Line.p),8),' [p.u.]'];
  OPF.report{5,1} = ['Bid Losses = ',fvar(sum(Line.p)-Snapshot(1).Ploss,8),' [p.u.]'];

  if ~noDem
    OPF.report{6,1} = ['Total demand = ', ...
                        fvar(sum(Pd),8),' [p.u.]'];
  end
  OPF.report{6+(~noDem),1} = ['TTL = ',fvar(sum(Pd)+TPQ,8),' [p.u.]'];

  fm_disp
  fm_disp('----------------------------------------------------------------')

  Settings.lftime = toc;

  if Fig.stat, fm_stat; end

  if iteration > iter_max
    fm_disp('IPM-OPF: Method did not Converge',2)
  elseif Fig.main
    if ~get(Fig.main,'UserData')
      fm_disp('IPM-OPF: Interrupted',2)
    else
      fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
    end
  else
    fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
    if Settings.showlf == 1
      fm_stat(OPF.report);
    else
      if Settings.beep
        beep
      end
    end
  end
  fm_status('opf','close')
else
  if iteration > iter_max
    fm_disp('IPM-OPF: Method did not Converge',2)
  elseif Fig.main
    if ~get(Fig.main,'UserData')
      fm_disp(['IPM-OPF: Interrupted'],2)
    else
      fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
    end
  else
    fm_disp(['IPM-OPF completed in ',num2str(toc),' s'],1)
  end
end

if iteration > iter_max,
  OPF.conv = 0;
else,
  OPF.conv = 1;
end
if Rsrv.n,
  OPF.guess = [s; mu; DAE.y; Qg; Ps; Pd; ...
               yc; kg; Qgc; lc; Pr; ro];
else
  OPF.guess = [s; mu; DAE.y; Qg; Ps; Pd; ...
               yc; kg; Qgc; lc; ro];
end
OPF.LMP = -ro(1:n1);
OPF.atc = (1+lc)*(sum(Pd)+TPQ)*MVA;
OPF.yc = yc;

if noDem, Demand = restore(Demand); end

Settings.forcepq = forcepq;
if ~OPF.basepl
  PQ = pqreset(PQ,'all');
end
if ~OPF.basepg
  SW = swreset(SW,'all');
  PV = pvreset(PV,'all');
end