mopf.m

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

else
  iang = find((branch(:, ANGMIN) & branch(:, ANGMIN) > -360) | ...
              (branch(:, ANGMAX) & branch(:, ANGMAX) < 360));
  iangl = find(branch(iang, ANGMIN));
  iangh = find(branch(iang, ANGMAX));
  nang = length(iang);
end

% Find out problem dimensions
nb = size(bus, 1);                              % buses
ng = size(gen, 1);                              % variable injections
nl = size(branch, 1);                           % branches
iycost = find(gencost(:, MODEL) == PW_LINEAR);  % y variables for pwl cost
ny    = size(iycost, 1);
neqc  = 2 * nb;                                 % nonlinear equalities
nusr  = size(Au, 1);                            % # linear user constraints
nx    = 2*nb + 2*ng;                            % control variables
nvl   = size(vload, 1);                         % dispatchable loads
npqh  = size(ipqh, 1);                          % general pq capability curves
npql  = size(ipql, 1);
if isempty(Au)
  nz = 0;
  Au = sparse(0,nx);
  if ~isempty(N)        % still need to check number of columns of N
    if size(N,2) ~= nx;
      error(sprintf('mopf.m: user supplied N matrix must have %d columns.', nx));
    end
  end
else
  nz = size(Au,2) - nx;                       % additional linear variables
  if nz < 0
    error(sprintf('mopf.m: user supplied A matrix must have at least %d columns.', nx));
  end
end
nxyz = nx+ny+nz;                                % total # of vars of all types

% Form vector of indexes into table columns; this vector is needed
% by minopf.mex$(PLATFORM)
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];

% Definition of indexes into optimization variable vector and constraint
% vector.
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;
zbas     = yend + 1;      zend     = zbas + nz - 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;
usrbas  = stend+1;        usrend  = usrbas+nusr-1; % warning: nusr could be 0
pqhbas  = usrend+1;       pqhend  = pqhbas+npqh-1; % warning: npqh could be 0
pqlbas  = pqhend+1;       pqlend  = pqlbas+npql-1; % warning: npql could be 0
vlbas   = pqlend+1;       vlend   = vlbas+nvl-1;   % warning: nvl could be 0
angbas  = vlend+1;        angend  = angbas+nang-1; % not done yet, need # of
                                                 % Ay constraints.

% Let makeAy deal with any y-variable for piecewise-linear convex costs.
% note that if there are z variables then Ay doesn't have the columns
% that would span the z variables, so we append them.
if ny > 0
  [Ay, by]  = makeAy(baseMVA, ng, gencost, pgbas, qgbas, ybas);
  if nz > 0
    Ay = [ Ay  sparse(size(Ay,1), nz) ];
  end
else
  Ay = [];
  by =[];
end
ncony = size(Ay,1);
yconbas = angend+1;       yconend = yconbas+ncony-1; % finally done with
                                                     % constraint indexing

% Make Aang, lang, uang for branch angle difference limits
if nang > 0
  ii = [(1:nang)'; (1:nang)'];
  jj = [branch(iang, F_BUS); branch(iang, T_BUS)];
  Aang = sparse(ii, jj, [ones(nang, 1); -ones(nang, 1)], nang, nxyz);
  uang = 1e10*ones(nang,1);
  lang = -uang;
  lang(iangl) = branch(iang(iangl), ANGMIN) * pi/180;
  uang(iangh) = branch(iang(iangh), ANGMAX) * pi/180;
else
  Aang =[];
  lang =[];
  uang =[];
end

% Make Avl, lvl, uvl in case there is a need for dispatchable loads
if nvl > 0
  xx = gen(vload, PMIN);
  yy = Qlim;
  pftheta = atan2(yy, xx);
  pc = sin(pftheta);
  qc = -cos(pftheta);
  ii = [ (1:nvl)'; (1:nvl)' ];
  jj = [ pgbas+vload-1; qgbas+vload-1 ];
  Avl = sparse(ii, jj, [pc; qc], nvl, nxyz);
  lvl = zeros(nvl, 1);
  uvl = lvl;
else
  Avl =[];
  lvl =[];
  uvl =[];
end

% Make Apqh if there is a need to add general PQ capability curves;
% use normalized coefficient rows so multipliers have right scaling
% in $$/pu
if npqh > 0
  Apqhdata = [gen(ipqh,QC1MAX)-gen(ipqh,QC2MAX), gen(ipqh,PC2)-gen(ipqh,PC1)];
  ubpqh = (gen(ipqh,QC1MAX)-gen(ipqh,QC2MAX)) .* gen(ipqh,PC1) ...
         + (gen(ipqh,PC2)-gen(ipqh,PC1)) .* gen(ipqh,QC1MAX);
  for i=1:npqh,
    tmp = norm(Apqhdata(i,:));
    Apqhdata(i,:) = Apqhdata(i, :) / tmp;
    ubpqh(i) = ubpqh(i) / tmp;
  end
  Apqh = sparse([1:npqh, 1:npqh]', [(pgbas-1)+ipqh;(qgbas-1)+ipqh], ...
                Apqhdata(:), npqh, nxyz);
  ubpqh = ubpqh / baseMVA;
  lbpqh = -1e10*ones(npqh,1);
else
  Apqh = [];
  ubpqh = [];
  lbpqh = [];
end

% similarly Apql
if npql > 0
  Apqldata = [gen(ipql,QC2MIN)-gen(ipql,QC1MIN), gen(ipql,PC1)-gen(ipql,PC2)];
  ubpql= (gen(ipql,QC2MIN)-gen(ipql,QC1MIN)) .* gen(ipql,PC1) ...
         - (gen(ipql,PC2)-gen(ipql,PC1)) .* gen(ipql,QC1MIN) ;
  for i=1:npql,
    tmp = norm(Apqldata(i, : ));
    Apqldata(i, :) = Apqldata(i, :) / tmp;
    ubpql(i) = ubpql(i) / tmp;
  end
  Apql = sparse([1:npql, 1:npql]', [(pgbas-1)+ipql;(qgbas-1)+ipql], ...
                Apqldata(:), npql, nxyz);
  ubpql = ubpql / baseMVA;
  lbpql = -1e10*ones(npql,1);
else
  Apql = [];
  ubpql = [];
  lbpql = [];
end

% Insert y columns in Au and N as necessary
if ny > 0
  if nz > 0
    Au = [ Au(:,1:qgend) sparse(nusr, ny) Au(:, qgend+(1:nz)) ];
    if ~isempty(N)
        N = [ N(:,1:qgend) sparse(size(N,1), ny) N(:, qgend+(1:nz)) ];
    end
  else
    Au = [ Au sparse(nusr, ny) ];
    if ~isempty(N)
        N = [ N sparse(size(N,1), ny) ];
    end
  end
end

% Now form the overall linear restriction matrix;
% note the order of the constraints.

if (ncony > 0 )
  A = [ Au;
        Apqh;
        Apql;
        Avl;
        Aang;
        Ay;
        sparse(ones(1,ny), ybas:yend, ones(1,ny), 1, nxyz ) ];  % "linear" cost
  l = [ lbu;
        lbpqh;
        lbpql;
        lvl;
        lang;
       -1e10*ones(ncony+1, 1) ];
  u = [ ubu;
        ubpqh;
        ubpql;
        uvl;
        uang;
        by;
        1e10];
else
  A = [ Au; Apqh; Apql; Avl; Aang ];
  l = [ lbu; lbpqh; lbpql; lvl; lang ];
  u = [ ubu; ubpqh; ubpql; uvl; uang ];
end


% Problem data for MINOPF is ready now.

% So, can we do anything good about lambda initialization?
if all(bus(:, LAM_P) == 0)
  bus(:, LAM_P) = (10)*ones(nb, 1);
end

% initialize some default options
if mpopt(61) == 0           % MNS_FEASTOL
  mpopt(61) = mpopt(16);    % = OPF_VIOLATION by default
end
if mpopt(62) == 0           % MNS_ROWTOL
  mpopt(62) = mpopt(16);    % = OPF_VIOLATION by default
end
if mpopt(63) == 0           % MNS_XTOL
  mpopt(63) = mpopt(17);    % = CONSTR_TOL_X by default
end

% Call minopf and hopefully MINOS will solve it on the first try.
[buso, geno, brancho, f, info, g, jac, xr, pimul] = minopf(baseMVA, bus,gen, ...
                      gencost, branch, col, A, l, u, mpopt, N, fparm, H, Cw, ...
                      z0, zl ,zu);

success = (info == 0);  % | (info == 9) | (info == 13);
if ~success & have_fcn('bpmpd')
  % oh oh, MINOS barfed; let's try to precondition the case a little bit.
  [bus1, gen1, branch1, f, info1 ] = sqppf(baseMVA, bus, gen, branch, ...
                                  areas, gencost, mpopt);
  [buso, geno, brancho, f, info, g, jac, xr, pimul] = minopf(baseMVA, bus1, gen1, ...
                      gencost, branch1, col, A, l, u, mpopt, N, fparm, H, Cw, ...
                      z0, zl, zu);
  success = (info == 0);    % | (info == 9) | (info == 13);
end

% If we succeeded and there were generators with general pq curve
% characteristics, this is the time to re-compute the multipliers,
% splitting any nonzero multiplier on one of the linear bounds among the
% Pmax, Pmin, Qmax or Qmin limits, producing one multiplier for a P limit and
% another for a Q limit. For upper Q limit, if we are neither at Pmin nor at 
% Pmax, the limit is taken as Pmin if the Qmax line's normal has a negative P
% component, Pmax if it has a positive P component. Messy but there really
% are many cases.  Remember multipliers in pimul() are negative.
muPmax = geno(:, MU_PMAX);
muPmin = geno(:, MU_PMIN);
if success & (npqh > 0)
  k = 1;
  for i = ipqh'
    if muPmax(i) > 0
      geno(i,MU_PMAX)=geno(i,MU_PMAX)-pimul(pqhbas+k-1)*Apqhdata(k,1)/baseMVA;
    elseif muPmin(i) > 0
      geno(i,MU_PMIN)=geno(i,MU_PMIN)+pimul(pqhbas+k-1)*Apqhdata(k,1)/baseMVA;
    else
      if Apqhdata(k, 1) >= 0
         geno(i,MU_PMAX)=geno(i,MU_PMAX)-pimul(pqhbas+k-1)*Apqhdata(k,1)/baseMVA;
      else
         geno(i,MU_PMIN)=geno(i,MU_PMIN)+pimul(pqhbas+k-1)*Apqhdata(k,1)/baseMVA;
      end
    end
    geno(i,MU_QMAX)=geno(i,MU_QMAX)-pimul(pqhbas+k-1)*Apqhdata(k,2)/baseMVA;
    k = k + 1;
  end
end

if success & (npql > 0)
  k = 1;
  for i = ipql'
    if muPmax(i) > 0
      geno(i,MU_PMAX)=geno(i,MU_PMAX)-pimul(pqlbas+k-1)*Apqldata(k,1)/baseMVA;
    elseif muPmin(i) > 0
      geno(i,MU_PMIN)=geno(i,MU_PMIN)+pimul(pqlbas+k-1)*Apqldata(k,1)/baseMVA;
    else
      if Apqldata(k,1) >= 0
        geno(i,MU_PMAX)=geno(i,MU_PMAX)-pimul(pqlbas+k-1)*Apqldata(k,1)/baseMVA;
      else
        geno(i,MU_PMIN)=geno(i,MU_PMIN)+pimul(pqlbas+k-1)*Apqldata(k,1)/baseMVA;
      end
    end
    geno(i,MU_QMIN)=geno(i,MU_QMIN)+pimul(pqlbas+k-1)*Apqldata(k,2)/baseMVA;
    k = k + 1;
  end
end

% angle limit constraints
if success & (nang > 0)
  tmp = [angbas:angend];
  ii = find(pimul(tmp) > 0);
  brancho(iang(ii), MU_ANGMIN) = pimul(tmp(ii)) * pi/180;
  ii = find(pimul(tmp) < 0);
  brancho(iang(ii), MU_ANGMAX) = -pimul(tmp(ii)) * pi/180;
end

% With these modifications, printpf must then look for multipliers
% if available in order to determine if a limit has been hit.

% Start preparation of output data.

% first reorder generators
geno = geno(inv_gen_ord, :);

% then copy gen bus voltages back to gen matrix
geno(:, VG) = buso(geno(:, GEN_BUS), VM);

% convert to original external bus ordering
[buso, geno, brancho, areas] = int2ext(i2e, buso, geno, brancho, areas);

% Now create output matrices with all lines, all generators, committed and
% non-committed
genout = genorg;
branchout = branchorg;
genout(comgen, : ) = geno;
branchout(onbranch, :)  = brancho;
% And zero out appropriate fields of non-comitted generators and lines
if ~isempty(offgen)
  tmp = zeros(length(offgen), 1);
  genout(offgen, PG) = tmp;
  genout(offgen, QG) = tmp;
  genout(offgen, MU_PMAX) = tmp;
  genout(offgen, MU_PMIN) = tmp;
end
if ~isempty(offbranch)
  tmp = zeros(length(offbranch), 1);
  branchout(offbranch, PF) = tmp;
  branchout(offbranch, QF) = tmp;
  branchout(offbranch, PT) = tmp;
  branchout(offbranch, QT) = tmp;
  branchout(offbranch, MU_SF) = tmp;
  branchout(offbranch, MU_ST) = tmp;
end

et = etime(clock,t1);
if (nargout == 0) & ( (info == 0) | (info == 9) | (info == 13) )
  printpf(baseMVA, buso, genout, branchout, f, success, et, 1, mpopt);
end

if nargout, busout = buso; end

return;