fm_opfm.m

电力系统的psat

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

% =======================================================================
% 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_Iijmax = mu(idx+nL);     idx = idx + Line.n;
  mu_Ijimax = mu(idx+nL);     idx = idx + Line.n;
  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(theta,V,Qg,Ps,Pd) = 0
  %________________________________________________________________________

  fm_lf(1); fm_pq(1);
  for i = 1:SW.n,  k = SW.bus(i);  DAE.gp(k) = DAE.gp(k) - Bus.Pg(k); end
  for i = 1:PV.n,  k = PV.bus(i);  DAE.gp(k) = DAE.gp(k) - PV.con(i,4); end
  for i = 1:n_gen, k = busG(i);    DAE.gq(k) = DAE.gq(k) - Qg(i); end

  for i = 1:Demand.n,
    k = Demand.bus(i);
    DAE.gp(k) = DAE.gp(k) + Pd(i);
    DAE.gq(k) = DAE.gq(k) + Pd(i)*qonp(i);
  end
  for i = 1:Supply.n,
    k = Supply.bus(i);
    DAE.gp(k) = DAE.gp(k) - Ps(i);
  end

  fm_lf(2);
  Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
  [Iij, Jij, Hij, Iji, Jji, Hji] = fm_flows(flow,mu_Iijmax, mu_Ijimax);

  % Gradient of [s] variables
  %________________________________________________________________________

  gs = s.*mu - ms;

  % Gradient of [mu] variables
  %________________________________________________________________________

  gmu = [Psmin-Ps;Ps-Psmax;Pdmin-Pd;Pd-Pdmax;Qgmin-Qg;Qg-Qgmax; ...
         Vmin-DAE.V;DAE.V-Vmax;Iij-Iijmax;Iji-Iijmax];
  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; DAE.V; Qg; Ps; Pd] variables
  %________________________________________________________________________

  Jg = [Jlfv, g_Qg, g_Ps, g_Pd];
  if Rsrv.n, Jg = [Jg, g_Pr]; end

  dF_dy = (Jg.')*ro;
  dG_dy(n2+n_gen+nS) = (Csb + 2*Csc.*Ps + 2*KTBS.*Ps);
  dG_dy(n2+n_gen+Supply.n+nD) = -(~noDem)*(Cdb+2*Cdc.*Pd+2*KTBD.* ...
                                           Pd+qonp.*(Ddb+2*qonp.*Ddc.*Pd));
  dG_dy(n2+busS) = (Dsb + 2*Dsc.*Qg(busS));
  dH_dtV = (Jij.')*mu_Iijmax + (Jji.')*mu_Ijimax + [mu_t; mu_Vmax-mu_Vmin];
  if Rsrv.n,
    dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; mu_Pdmax-mu_Pdmin-mu_sumPrd; ...
             mu_Psmax+mu_Prmax-mu_Prmin+mu_sumPrd];
  else
    dH_dy = [dH_dtV; mu_Qgmax-mu_Qgmin; mu_Psmax-mu_Psmin; mu_Pdmax-mu_Pdmin];
  end
  gy = dG_dy - dF_dy + dH_dy;

  Jh(n_b+nL,1:n2) = Jij;
  Jh(n_b+Line.n+nL,1:n2) = Jji;

  % Hessian Matrix [D2xLms]
  %________________________________________________________________________

  if Settings.octave
    H3 = zeros(n_a+Rsrv.n,n_y);
  else
    H3 = sparse(n_a+Rsrv.n,n_y);
  end
  H3 = H3 - sparse(n_gen+nS,n2+n_gen+nS,(2*Csc+2*KTBS),n_a+Rsrv.n,n_y);
  H3 = H3 - sparse(busS,n2+busS,2*Dsc,n_a+Rsrv.n,n_y);
  H3 = H3 + sparse(n_gen+Supply.n+nD,n_gen+Supply.n+n2+nD, ...
                   (2*Cdc+2*KTBD+2*Ddc.*qonp.*qonp),n_a+Rsrv.n,n_y);

  Hess = -fm_hessian(ro(1:n2))+Hij+Hji;
  D2xLms = [Hess, H31; -H3];

  % Complete Hessian Matrix [D2yLms]
  %________________________________________________________________________

  switch method
   case 1 % Newton Directions

    % reduced system
    if Settings.octave
      H_m  = diag(mu./s);
      H_s = diag(1./s);
    else
      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);
    end
    Jh(:,SW.bus) = 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.bus,:) = 0;
    Jd(:,SW.bus) = 0;
    if Settings.octave
      Jd(SW.bus,SW.bus) = eye(SW.n);
    else
      Jd(SW.bus,SW.bus) = speye(SW.n);
    end
    gy(SW.bus) = 0;
    % variable increments
    Dx = -Jd\[gy; -DAE.gp; -DAE.gq];
    Ds = -(gmu+Jh*Dx([1:n_y]));
    Dm = -H_s*gs-H_m*Ds;

   case 2 % Mehrotra's Predictor-Corrector

    % -------------------
    % Predictor step
    % -------------------
    % reduced system
    if Settings.octave
      H_m  = diag(mu./s);
    else
      H_m  = sparse(1:n_s,1:n_s,mu./s,n_s,n_s);
    end
    Jh(:,SW.bus) = 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.bus,:) = 0;
    Jd(:,SW.bus) = 0;
    if Settings.octave
      Jd(SW.bus,SW.bus) = eye(SW.n);
    else
      Jd(SW.bus,SW.bus) = speye(SW.n);
    end
    gx(SW.bus) = 0;
    % LU factorization
    [L,U,P] = lu(Jd);
    % variable increments
    Dx = -U\(L\(P*[gx; -DAE.gp; -DAE.gq]));
    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.bus) = 0;
    % variable increments
    Dx = -U\(L\(P*[gx; -DAE.gp; -DAE.gq]));
    Ds = -(gmu+Jh*Dx([1:n_y]));
    Dm = -gs-H_m*Ds;
  end

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

  Dtheta = Dx(nB);                idx = Bus.n;            % curr. sys.
  DV     = Dx(idx+nB);            idx = idx + Bus.n;
  DQg    = Dx(idx+nG);            idx = idx + n_gen;
  DPs    = Dx(idx+nS);            idx = idx + Supply.n;
  DPd    = Dx(idx+nD);            idx = idx + Demand.n;
  if Rsrv.n, DPr = Dx(idx+nR);    idx = idx + Rsrv.n;    end
  Dro    = Dx(1+idx:end);                                 % Lag. mult.

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

  % Step Length 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.a = DAE.a + alpha_P*Dtheta;
  DAE.V = DAE.V + alpha_P*DV;
  Ps = Ps + alpha_P*DPs;
  Pd = Pd + alpha_P*DPd;
  Qg = Qg + alpha_P*DQg;
  if Rsrv.n, Pr = Pr + alpha_P*DPr; end
  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  - (~noDem)*Cdb'*Pd + Dsb'*Qg(busS) - 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(busS).*Qg(busS));
  if Rsrv.n, Reserve = Cr'*Pr; else, Reserve = 0; end
  G_obj = (Fixd_c+Prop_c+Quad_c+Quad_q+TieBreaking+Reserve)-ms*sum(log(s));

  % =======================================================================
  % 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.gp; DAE.gq],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;

  fm_status('opf','update',[iteration, ms, norma2, norma3, norma4], ...
            iteration)

  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
  if iteration > iter_max, break, end
end

% Updating Demand.con & Supply.con
%____________________________________________________________________________

Demand.con(:,7) = Pd;
Supply.con(:,6) = Ps;
if Rsrv.n, Rsrv.con(:,10) = Pr; end
MVA = Settings.mva;
Pay = ro(1:Bus.n).*DAE.glfp*MVA;
ISOPay = sum(Pay);
Iij  = sqrt(Iij);
Iji  = sqrt(Iji);
Iijmax  = sqrt(Iijmax);

% Computation of Nodal Congestion Prices (NCPs)
%____________________________________________________________________________

Jlfv(SW.bus,:) = [];
Jlfv(:,SW.bus) = [];
dH_dtV(SW.bus,:) = [];
NCP = -Jlfv'\dH_dtV;
OPF.NCP = [NCP(1:SW.bus-1);0;NCP(SW.bus:end)];

OPF.obj = G_obj;
OPF.ms = ms;
OPF.dy = norma2;
OPF.dF = norma3;
OPF.dG = norma4;
OPF.iter = iteration;

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

% set Pg, Qg, Pl and Ql
for i = 1:Demand.n,
  k = Demand.bus(i);
  Bus.Pl(k) = Snapshot(1).Pl(k)+Pd(i);
  Bus.Ql(k) = Snapshot(1).Ql(k)+Pd(i)*qonp(i);
end

for i = 1:Supply.n
  k = Supply.bus(i);
  Bus.Pg(k) = Snapshot(1).Pg(k)+Ps(i);
end

for i = 1:n_gen
  Bus.Qg(nG(i)) = Qg(i);
end

% Display Results
%____________________________________________________________________________

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

  OPF.report = cell(1,1);
  OPF.report{1,1} = ['TTL = ', ...
                     fvar(sum(Pd)+sum(PQ.con(:,4)),8), ' [p.u.]'];
  if ~noDem
    OPF.report{2,1} = ['Total demand = ',fvar(sum(Pd),8), ' [p.u.]'];
  end
  OPF.report{2+(~noDem),1} = ['Bid Losses = ', ...
                      fvar(sum(DAE.glfp)-Snapshot(1).Ploss,8), ' [p.u.]'];
  OPF.report{3+(~noDem),1} = ['Total Losses = ', ...
                      fvar(sum(DAE.glfp),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.a; DAE.V; Qg; Ps; Pd; Pr; ro];
else
  OPF.guess = [s; mu; DAE.a; DAE.V; Qg; Ps; Pd; ro];
end

if noDem
  Demand.con = [];
  Demand.n = 0;
  Demand.bus = [];
end

if ~OPF.basepl
  Bus.Pl = buspl;
  Bus.Ql = busql;
  PQ.con(:,[4 5]) = pqcon;
end
if ~OPF.basepg
  Snapshot(1).Ploss = ploss;
  Bus.Pg = buspg;
  Bus.Qg = busqg;
  PV.con(:,4) = pvcon;
end