fm_opfsdr.m

电力系统的psat

  Jh = sparse(n_s,n_y);
end

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
  if Settings.octave
    g_Pr  = zeros(n2,Rsrv.n);
  else
    g_Pr  = sparse(n2,Rsrv.n);
  end
  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

if Settings.octave
  Jg = zeros(n4,n_y);
else
  Jg = sparse(n4,n_y);
end

% Hessian Matrices but Jlfv
%____________________________________________________________________________

if Settings.octave
  I_smu = eye(n_s);
  Z1 = zeros(n_s,n_y);
  Z2 = zeros(n_s,n_s);
  Z3 = zeros(n4,n4);
  Z4 = zeros(n4,n_s);
  Hcolk = zeros(n2,1);
  Hrowk = zeros(1,n2+1);
  H31 = zeros(n2,n_gen+n_a+n2+2+Rsrv.n);
  H41 = zeros(n2+1,n2+n_a);
  H42 = zeros(n2+1,n_gen+1+Rsrv.n);
  H5  = zeros(n_gen+1+Rsrv.n,n_y);
else
  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);
end

% ===========================================================================
% 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 = Bus.Qg(busG);
%Qg = Qgmin + 0.5*(Qgmax-Qgmin);
%DAE.V = Vmin + 0.5*(Vmax-Vmin);
%Vc = DAE.V;
%Iij = 0.1*Iijmax;
%Iji = 0.1*Iijmax;
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 = [];
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
  fm_disp('Optimization routine interrupted.',1)
  return
end
a = find(Qg > Qgmax);
b = find(Qg < Qgmin);
if ~isempty(a),
  fm_disp(['Max Reactive Power limit not respected at buses:'],1)
  for i = 1:length(a)
    fm_disp(['     Bus #',fvar(Bus.con(busG(a(i)),1),4),' -> ',fvar(Qg(a(i)),8), ...
             ' > ', fvar(Qgmax(a(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
if ~isempty(b),
  fm_disp(['Min Reactive Power limit not respected at buses: ',num2str(b')],1)
  for i = 1:length(b)
    fm_disp(['     Bus #',fvar(Bus.con(busG(b(i)),1),4),' -> ',fvar(Qg(b(i)),8), ...
             ' < ', fvar(Qgmin(b(i)),8)],1)
  end
  fm_disp('Optimization routine interrupted.',1)
  return
end
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

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

% =======================================================================
% 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;
  mu_Iijcmax= mu(idx+nL);     idx = idx + Line.n;
  mu_Ijimax = mu(idx+nL);     idx = idx + Line.n;
  mu_Ijicmax= 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(thetac,Vc,Qgc,Ps,Pd,lc) = 0
  %________________________________________________________________________

  V_snap = DAE.V;     ang_snap = DAE.a;
  DAE.V  = Vc;        DAE.a    = angc;
  g_lc  = sparse(n2,1);
  g_kg  = sparse(n2,1);

  if OPF.line, Line.Y = Y_cont; Line.con = line_cont; end
  fm_lf(1);
  for i = 1:PQ.n
    k = PQ.bus(i);
    DAE.gp(k) = (1+lc)*PQ.con(i,4) + DAE.gp(k);
    DAE.gq(k) = (1+lc)*PQ.con(i,5) + DAE.gq(k);
    g_lc(k) = PQ.con(i,4);
    g_lc(Bus.n+k) = PQ.con(i,5);
  end
  for i = 1:SW.n,
    k = SW.bus(i);
    DAE.gp(k) = DAE.gp(k) - (1+lc+kg*ksw(i))*Bus.Pg(k);
    g_lc(k) = -Bus.Pg(k);
    g_kg(k) = -ksw(i)*Bus.Pg(k);
  end
  for i = 1:PV.n,
    k = PV.bus(i);
    DAE.gp(k) = DAE.gp(k) - (1+lc+kg*kpv(i))*PV.con(i,4);
    g_lc(k) = -PV.con(i,4);
    g_kg(k) = -kpv(i)*PV.con(i,4);
  end
  for i = 1:n_gen,
    k = busG(i);
    DAE.gq(k) = DAE.gq(k) - Qgc(i);
  end

  for i = 1:Demand.n,
    k = Demand.bus(i);
    DAE.gp(k) = DAE.gp(k) + (1+lc)*Pd(i);
    DAE.gq(k) = DAE.gq(k) + (1+lc)*Pd(i)*qonp(i);
    g_Pdc(k,i) = 1+lc;
    g_Pdc(Bus.n+k,i) = (1+lc)*qonp(i);
    g_lc(k) = g_lc(k)+Pd(i);
    g_lc(Bus.n+k) = g_lc(Bus.n+k) + qonp(i)*Pd(i);
  end
  for i = 1:Supply.n,
    k = Supply.bus(i);
    DAE.gp(k) = DAE.gp(k) - (1+lc+kg*ksu(i))*Ps(i);
    g_Psc(k,i) = -(1+lc+kg*ksu(i));
    g_lc(k) = g_lc(k)-Ps(i);
    g_kg(k) = g_kg(k)-ksu(i)*Ps(i);
  end
  gc1p = DAE.gp; gc2p = DAE.glfp;
  gc1q = DAE.gq; gc2q = DAE.glfq;
  fm_lf(2);
  Jlfvc = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
  [Iijc, Jijc, Hijc, Ijic, Jjic, Hjic] = fm_flows(OPF.flow,mu_Iijcmax, mu_Ijicmax);

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

  DAE.V = V_snap;
  DAE.a = ang_snap;

  if OPF.line, Line.Y = Y_orig; Line.con = line_orig; end
  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(OPF.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;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; DAE.V; Qg; Ps; Pd] variables
  %________________________________________________________________________

  if Rsrv.n,
    Jg = [Jlfv,g_Qg,g_Ps,g_Pd,Jz1;Jz2,g_Psc,g_Pdc,Jlfvc,g_kg,g_Qgc,g_lc,g_Pr];
  else,
    Jg = [Jlfv,g_Qg,g_Ps,g_Pd,Jz1;Jz2,g_Psc,g_Pdc,Jlfvc,g_kg,g_Qgc,g_lc];
  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+busS) = (1-w)*(Dsb + 2*Dsc.*Qg(busS));

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