sqppf.m

一种基于MATLAB的多目标最优化仿真工具箱

function [buso, geno, brancho, f, info, et, g, dg] = sqppf(baseMVA, bus, ...
                                    gen, branch, areas, gencost, option, Au, lu, uu)
% SQPPF:  Preconditioner for MINOPF.
%
%Usage:
%
% [bus, gen, branch, f, info] = sqppf(baseMVA, bus, gen, branch, areas, ...
%                                     gencost, mpopt)
%
% [bus, gen, branch, f, info, et, g, dg] = sqppf(baseMVA, bus, gen, ...
%                             branch, areas, gencost, mpopt, Au, lu, uu)

%   MINOPF for MATPOWER
%   $Id: sqppf.m,v 1.7 2006/02/13 19:29:47 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   Copyright (c) 2000-2005 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/minopf/ for more info.


% This function started as an SQP opf but ended being just a
% preconditioner for MINOS. In its current state, it should be
% considered just a quick hack.  Don't look too intensely at the
% code or you will make me blush :-) .


t1 = clock;
if isstr(baseMVA)
  casefile = baseMVA;
  [baseMVA, bus, gen, branch, areas, gencost] = loadcase(casefile);
end

if nargin < 7
  option = mpoption;
end

verbose = option(31);

if nargin < 10
  Au = [];
  lu = [];
  uu = [];
end

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

% If tables do not have multiplier/extra columns, append zero cols.
% Update whenever the data format changes!
if size(bus,2) < MU_VMIN
  bus = [bus zeros(size(bus,1),MU_VMIN-size(bus,2)) ];
end
if size(gen,2) < MU_QMIN
  gen = [ gen zeros(size(gen,1),MU_QMIN-size(gen,2)) ];
end
if size(branch,2) < MU_ANGMAX
  branch = [ branch zeros(size(branch,1),MU_ANGMAX-size(branch,2)) ];
end

% Filter out inactive generators and branches

comgen = find(gen(:,GEN_STATUS) > 0);
offgen = find(gen(:,GEN_STATUS) <= 0);
onbranch  = find(branch(:,BR_STATUS) ~= 0);
offbranch = find(branch(:,BR_STATUS) == 0);

[i2e, bus, gen, branch, areas] = ext2int(bus, gen, branch, areas);
[ref, pv, pq] = bustypes(bus, gen);


c1 = [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
      VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN];
c2 = [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
        GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN];
c3 = [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST];
c4 = [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
        RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST];

col = [ c1 c2 c3 c4];

nb = size(bus, 1);
ng = size(gen, 1);
nl = size(branch, 1);
iycost = find(gencost(:, MODEL) == PW_LINEAR);
ny = size(iycost, 1);
neqc = 2 * nb;
nx = 2*nb + 2*ng;

thbas = 1;                thend    = thbas+nb-1;
vbas     = thend+1;       vend     = vbas+nb-1;
pgbas    = vend+1;        pgend    = pgbas+ng-1;
qgbas    = pgend+1;       qgend    = qgbas+ng-1;
ybas     = qgend + 1;     yend     = ybas + ny - 1;

pmsmbas = 1;              pmsmend = pmsmbas+nb-1;
qmsmbas = pmsmend+1;      qmsmend = qmsmbas+nb-1;
sfbas   = qmsmend+1;      sfend   = sfbas+nl-1;
stbas   = sfend+1;        stend   = stbas+nl-1;

[Ay, by]  = makeAy(baseMVA, ng, gencost, pgbas, qgbas, ybas);

ncy = size(Ay,1);


x=zeros(nx,1);
x(thbas:thend) = bus(:,VA)*pi/180;
x(vbas:vend)   = bus(:,VM);
x(pgbas:pgend) = gen(:,PG)/baseMVA;
x(qgbas:qgend) = gen(:,QG)/baseMVA;

% Some tolerances and other algorithm behavior parameters
maxit  = 15;
dxstop = 0.0001;    % step size termination criterion for dx
gtol   = 0.0005;    % allowable constraint violation
max_dthlp = 0.1;  % max excursion for (ow. unconstrained) theta vars in lp steps
max_dvlp  = 0.01; % as above for voltages,
max_dplp  = 0.05;  % for active power outputs,
max_dqlp  = 0.05;  % and for reactive power outputs.
geweight  = 1000;  % weight on penalty function of equality constraints
W         = geweight * sparse(1:neqc,1:neqc,ones(1,neqc));
KV = 20000;       % penalty constant for voltage violations 0.5*KV*(u(V-Vmax))^2
%KV = 2000;       % penalty constant for voltage violations 0.5*KV*(u(V-Vmax))^2
KPG = 20000;      % penalty constant for Pg limits, 0.5*KPG*(u(PG-PGMAX))^2
%KPG = 2000;      % penalty constant for Pg limits, 0.5*KPG*(u(PG-PGMAX))^2
KQG = 20000;
%KQG = 2000;
KSL  = 500;      % penalty constant for |Sf| and |St| constraints
VMAXEPS = 0.005;   % penalty const. for vmax raised til violations < this number
PMAXEPS = 0.01;
QMAXEPS = 0.01;
alpha0  = geweight/50;
alpha   = 2.7;    % Levenberg-Marquardt Hessian positive-definization parm

A = [ Ay;
      Au ];
l = [ -1e10*ones(ncy,1) ;
       lu ];
u = [  by ;
       uu ];

% Eval constraints and Jacobian at starting point
[dum1, dum2, dum3, f, info, g, dg ] = minopf(baseMVA, bus, gen, ...
   gencost, branch, col, A, l, u, [], 1);
if verbose
  fprintf('\nsqppf: sequential quadratic programming constrained PF solver\n\n');
  fprintf('\n  Itn        F            |dx|          |g|       max|g|  constraint  Step\n');
  fprintf('==========================================================================\n');
end
ge     = g(1:neqc);
dge    = dg(pmsmbas:qmsmend, :);
normge = norm(ge, 2);
[normge1, imax]= max(abs(ge));
gS     = max(g(sfbas:stend), zeros(2*nl,1));
normgs = norm(gS,2);
normgs1= norm(gS, inf);
gV = max(x(vbas:vend)-bus(:,VMAX), zeros(nb,1)) ...
     - max(bus(:,VMIN)-x(vbas:vend), zeros(nb,1));
gP = max(x(pgbas:pgend)-gen(:,PMAX)/baseMVA, zeros(ng,1)) ...
     - max(gen(:,PMIN)/baseMVA-x(pgbas:pgend), zeros(ng,1));
gQ = max(x(qgbas:qgend)-gen(:,QMAX)/baseMVA, zeros(ng,1)) ...
     - max(gen(:,QMIN)/baseMVA-x(qgbas:qgend), zeros(ng,1));
k= 0;

%if (normge/nb > 0.1) | (normge1 > 0.5) | (normgs/nl > 0.5) | (normgs1 > 0.5)
if 1
  next_step = 'obj_ssqr';  % start by minimizing sum of squares of constraints
  fssq = 0.5*ge'*W*ge + 0.5*KSL*gS'*gS + 0.5*KV*gV'*gV + ...
         0.5*KPG*gP'*gP + 0.5*KQG*gQ'*gQ ;
  if verbose
    fprintf(' %3i   %8.6e   %8.6e   %10.6f  %9.5f %5i  ', ...
       k, fssq,  0, normge, normge1, imax);
    fprintf('  none yet\n');
  end
else
  next_step = 'obj_cost';  % Go straight to LP step + powerflow iterations
end

terminate = 0;
k = 1;
iter = 1;
while ~terminate
  if strcmp(next_step, 'obj_ssqr')
    % Do a sum-of-squares-of-constraints reduction.  Assume g(pmsmbas:stend),
    % ge, dge, gs, gV, gP, gQ, fssq are defined.
    gVs  = sign(gV);
    gPs  = sign(gP);
    gQs  = sign(gQ);
    dgother=[sparse(nb,nb)   spdiags(gVs,0,nb,nb)     sparse(nb,ng)   sparse(nb,ng)  ;
     sparse(ng,nb)   sparse(ng,nb)   spdiags(gPs,0,ng,ng)     sparse(ng,ng)  ;
     sparse(ng,nb)   sparse(ng,nb)   sparse(ng,ng)   spdiags(gQs,0,ng,ng)    ];
    Wother= spdiags([ KV*ones(nb,1); KPG*ones(ng,1); KQG*ones(ng,1) ],0,2*ng+nb,2*ng+nb);
    dfssq =   dge' * W * ge ...
              + KSL * dg(sfbas:stend,   : )'  * gS ...
              + dgother' * Wother * [ gV; gP; gQ];
    Q     =   (alpha0*k^0.5) * speye(nx) ...
           + dge' * W * dge ...
           + KSL * dg(sfbas:stend,:)' * spdiags(sign(gS),0,2*nl,2*nl) * dg(sfbas:stend,:) ...
           + spdiags([zeros(nb,1); KV*abs(gVs); KPG*abs(gPs); KQG*abs(gQs)],0,2*nb+2*ng,2*nb+2*ng);
    llist = [ (vbas:vend)';
              (pgbas:pgend)';
              (qgbas:qgend)' ];
    lval  = [ bus(:,VMIN)-x(vbas:vend);
              gen(:,PMIN)/baseMVA-x(pgbas:pgend);
              gen(:,QMIN)/baseMVA-x(qgbas:qgend) ];
    uval  = [ bus(:,VMAX)-x(vbas:vend);
              gen(:,PMAX)/baseMVA-x(pgbas:pgend);
              gen(:,QMAX)/baseMVA-x(qgbas:qgend) ];
    Aqp   = sparse( 1, ref, 1, 1, qgend);
    bqp   = 0;
    e     = 0;
    fnew = fssq + 1;
    i = 0;
    while (fnew > fssq) & (iter <= maxit)
      iter = iter + 1;
      VUB = 1e30 * ones(length(dfssq), 1);
      VLB = -VUB;
      VLB(llist) = lval;
      VUB(llist) = uval;
      [dx, lambdaout, how, success] = mp_qp(Q,dfssq,Aqp,bqp,VLB,VUB,[],length(bqp),-1);

      xnew = x + dx;
      bus(:,VA) = xnew(thbas:thend) * 180/pi;
      bus(:,VM) = xnew(vbas:vend);
      gen(:,PG) = xnew(pgbas:pgend) * baseMVA;
      gen(:,QG) = xnew(qgbas:qgend) * baseMVA;
      [dum1, dum2, dum3, f, info, g, dg ] = minopf(baseMVA, bus, gen, ...
                                   gencost, branch, col, A, l, u, [], 1);
      ge     = g(1:neqc);
      dge    = dg(pmsmbas:qmsmend, :);
      normge = norm(ge, 2);
      [normge1, imax]= max(abs(ge));
      gS = max(g(sfbas:stend), zeros(2*nl,1));
      normgs = norm(gS, 2);
      normgs1= norm(gS, inf);
      gV = max(x(vbas:vend)-bus(:,VMAX), zeros(nb,1)) ...
           - max(bus(:,VMIN)-x(vbas:vend), zeros(nb,1));
      gP = max(x(pgbas:pgend)-gen(:,PMAX)/baseMVA, zeros(ng,1)) ...
           - max(gen(:,PMIN)/baseMVA-x(pgbas:pgend), zeros(ng,1));
      gQ = max(x(qgbas:qgend)-gen(:,QMAX)/baseMVA, zeros(ng,1)) ...
           - max(gen(:,QMIN)/baseMVA-x(qgbas:qgend), zeros(ng,1));
      fnew = 0.5*ge'*W*ge + 0.5*KSL*gS'*gS + 0.5*KV*gV'*gV + ...
           0.5*KPG*gP'*gP + 0.5*KQG*gQ'*gQ ;
      if verbose
        fprintf(' %3i   %8.6e   %8.6e   %10.6f  %9.5f %5i  ', ...
               k, fnew, norm(dx), normge, normge1, imax);
      end
      if fnew > fssq
        i = i + 1;
        Q = Q + (alpha^i) * 100 * speye(nx);
        if verbose
          fprintf('  Q + a*I next\n');
        end
      else
        if verbose
          fprintf('  SSQ step\n');
        end
        k = k + 1;
      end
    end
    fssq = fnew;
    x = xnew;
    if ((normge1 < gtol) &  (normgs1 < gtol)) | (iter > maxit)
      terminate = 1;
    end
  end
  if strcmp(next_step, 'obj_cost')

  % Do a power flow



  % If the power flow did not converge, do QP w/o linearized equality
  % constraints and large  weight on penalty function of equality constraints;
  % otherwise, do QP with equality constraints.

  end

end

t2 = etime(t1,clock);
[dum1, geno, brancho, f, info, g, dg ] = minopf(baseMVA, bus, gen, ...
                                   gencost, branch, col, A, l, u, [], 1);
f = fssq;
info = how;

if nargout
  buso = dum1;
end