fmincopf.m

最优潮流计算

function [busout, genout, branchout, f, success, info, et, g, jac] = ...
      fmincopf(baseMVA, bus, gen, branch, areas, gencost, Au, lbu, ubu, mpopt)
%FMINCOPF  Solves an AC optimal power flow using FMINCON (Opt Tbx 2.x & later).
%
%   [bus, gen, branch, f, success] = fmincopf(casefile, mpopt)
%
%   [bus, gen, branch, f, success] = fmincopf(casefile, A, l, u, mpopt)
%
%   [bus, gen, branch, f, success] = fmincopf(baseMVA, bus, gen, branch,...
%                                    areas,  gencost, mpopt)
%
%   [bus, gen, branch, f, success] = fmincopf(baseMVA, bus, gen, branch,...
%                                    areas, gencost, A, l, u, mpopt)
%
%   [bus, gen, branch, f, success, info, et, g, jac] = fmincopf(casefile)
%
%   The data for the problem can be specified in one of 3 ways: (1) the name of
%   a case file which defines the data matrices baseMVA, bus, gen, branch,
%   areas and gencost, (2) a struct containing the data matrices as fields, or
%   (3) the data matrices themselves.
%
%   When specified, A, l, u represent additional linear constraints on the
%   optimization variables, l <= A*x <= u. For an explanation of the
%   formulation used and instructions for forming the A matrix, type
%   'help genform'.
%
%   The optional mpopt vector specifies MATPOWER options. Type 'help mpoption'
%   for details and default values.
%
%   The solved case is returned in the data matrices, bus, gen and branch. Also,
%   returned are the final objective function value (f) and a flag which is
%   true if the algorithm was successful in finding a solution (success).
%   Additional optional return values are an algorithm specific return status
%   (info), elapsed time in seconds (et), the constraint vector (g) and the
%   Jacobian matrix (jac).
% 
%   Rules for A matrix: If the user specifies an A matrix that has more columns
%   than the combined number of "x" (OPF) and "y" (pwl cost) variables, thus
%   allowing for extra linearly costed and linearly constrained "z" variables,
%   then it is the user's responsibility to fill the overall linear cost in the
%   last row of A, including that reflected in the "y" variables.
%
%   NOTE: The shadow prices (lambda's and mu's) produced by fmincon appear to
%         be slightly inaccurate.

%   MATPOWER
%   $Id: fmincopf.m,v 1.9 2005/01/25 14:39:46 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/matpower/ for more info.

% Sort out input arguments
t1 = clock;
if isstr(baseMVA) | isstruct(baseMVA)
  casefile = baseMVA;
  if nargin == 5
    Au = bus;
    lbu = gen;
    ubu = branch;
    mpopt = areas;
  elseif nargin == 4
    Au = bus;
    lbu = gen;
    ubu = branch;
    mpopt = mpoption;
  elseif nargin == 2
    Au = sparse(0,0);
    lbu = [];
    ubu = [];
    mpopt = bus;
  elseif nargin == 1
    Au = sparse(0,0);
    lbu = [];
    ubu = [];
    mpopt = mpoption;
  else
    error('fmincopf.m: Incorrect input parameter order, number or type');
  end;
  [baseMVA, bus, gen, branch, areas, gencost] = loadcase(casefile);
else
  if nargin == 9 
    mpopt = mpoption;
  elseif nargin == 7
    mpopt = Au;
    Au = sparse(0,0);
    lbu = [];
    ubu = [];
  elseif nargin == 6
    mpopt = mpoption;
    Au = sparse(0,0);
    lbu = [];
    ubu = [];
  elseif nargin ~= 10
    error('fmincopf.m: Incorrect input parameter order, number or type');
  end
end

% Load column indexes for case tables.
[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;
[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] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
        GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

% If tables do not have multiplier/extra columns, append zero cols
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_ST
  branch = [ branch zeros(size(branch,1),MU_ST-size(branch,2)) ];
end

% Filter out inactive generators and branches; save original bus & branch

comgen = find(gen(:,GEN_STATUS) > 0);
offgen = find(gen(:,GEN_STATUS) <= 0);
onbranch  = find(branch(:,BR_STATUS) ~= 0);
offbranch = find(branch(:,BR_STATUS) == 0);
genorg = gen;
branchorg = branch;
ng = size(gen,1);         % original size(gen), at least temporally
gen   = gen(comgen, :);
branch = branch(onbranch, :);
if size(gencost,1) == ng
  gencost = gencost(comgen, :);
else
  gencost = gencost( [comgen; comgen+ng], :);
end

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

% Sort generators in order of increasing bus number;
ng = size(gen,1);
[tmp, igen] = sort(gen(:, GEN_BUS));
[tmp, inv_gen_ord] = sort(igen);  % save for inverse reordering at the end
gen  = gen(igen, :);
if ng == size(gencost,1)
  gencost = gencost(igen, :);
else
  gencost = gencost( [igen; igen+ng], :);
end

% Find out if any of these "generators" are actually dispatchable loads.
% (see 'help isload' for details on what consitutes a dispatchable load)
% Dispatchable loads are modeled as generators with an added constant
% power factor constraint. The power factor is derived from the original
% value of Pmin and either Qmin (for inductive loads) or Qmax (for capacitive
% loads). If both Qmin and Qmax are zero, this implies a unity power factor
% without the need for an additional constraint.
vload = find( isload(gen) & (gen(:, QMIN) ~= 0 | gen(:, QMAX) ~= 0) );
% At least one of the Q limits must be zero (corresponding to Pmax == 0)
if any( gen(vload, QMIN) ~= 0 & gen(vload, QMAX) ~= 0 )
	error('Either Qmin or Qmax must be equal to zero for each dispatchable load.');
end
% Initial values of PG and QG must be consistent with specified power factor
% This is to prevent a user from unknowingly using a case file which would
% have defined a different power factor constraint under a previous version
% which used PG and QG to define the power factor.
Qlim = (gen(vload, QMIN) == 0) .* gen(vload, QMAX) + ...
    (gen(vload, QMAX) == 0) .* gen(vload, QMIN);
if any( abs( gen(vload, QG) - gen(vload, PG) .* Qlim ./ gen(vload, PMIN) ) > 1e-4 )
    errstr = sprintf('%s\n', ...
        'For a dispatchable load, PG and QG must be consistent', ...
        'with the power factor defined by PMIN and the Q limits.' );
    error(errstr);
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
nx = 2*nb + 2*ng;                              % control variables
nvl = size(vload, 1);                          % dispatchable loads
nz = size(Au,2) - 2*nb - 2*ng - ny;            % number of extra z variables
nz = max(nz,0);

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

% Let makeAy deal with any y-variable for piecewise-linear convex costs.
[Ay, by]  = makeAy(baseMVA, ng, gencost, pgbas, qgbas, ybas);
ncony = size(Ay,1);

% 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, yend);
  lvl = zeros(nvl, 1);
  uvl = lvl;
else
  Avl =[];
  lvl =[];
  uvl =[];
end

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

if (nz > 0)
  % user defined z variables thus becoming responsible for
  % defining any Ay needed as well as the cost row (whether there are
  % y variables or not) in the last row of Au.