copf.m

matlab语言 编写的电力系统仿真代码

function [busout, genout, branchout, f, success, info, et, g, jac] = copf( ...
                baseMVA, bus, gen, branch, areas, gencost, mpopt)
%COPF  Solves an AC optimal power flow using CONSTR (Opt Tbx 1.x & 2.x).
%
%   Uses the standard/CCV formulations from pre-3.0 versions of MATPOWER, as
%   opposed to the generalized formulation used by fmincopf and MINOPF.
%
%   [bus, gen, branch, f, success] = copf(casefile, mpopt)
%
%   [bus, gen, branch, f, success] = copf(baseMVA, bus, gen, branch, areas, ...
%                                    gencost, mpopt)
%
%   [bus, gen, branch, f, success, info, et, g, jac] = copf(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.
%
%   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).
%
%   Note: the original constr() in The Mathworks' Optimization Toolbox 1.x
%   had a bug in its handling of negative Lagrangian multipliers. For a while,
%   MathWorks distributed a fixed version of constr.m, but it no longer appears
%   to be available on their web and ftp sites.

%   MATPOWER
%   $Id: copf.m,v 1.4 2005/01/21 20:43:19 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

% input arg sorting
t1 = clock;
if isstr(baseMVA) | isstruct(baseMVA)
  casefile = baseMVA;
  if nargin == 1
    mpopt = mpoption;
  elseif nargin == 2
    mpopt = bus;
  else
    error('copf.m: Incorrect input parameter order, number or type');
  end
  [baseMVA, bus, gen, branch, areas, gencost] = loadcase(casefile);
else
  if nargin == 6
    mpopt = mpoption;
  elseif nargin ~= 7
    error('copf.m: Incorrect input parameter order, number or type');
  end
end

% Define constants
j = sqrt(-1);

% options
verbose = mpopt(31);
npts = mpopt(14);

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

%% get bus index lists of each type of bus
[ref, pv, pq] = bustypes(bus, gen);

%% build admittance matrices
[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);

% Check costs
if any( (gencost(:, MODEL) ~= PW_LINEAR) & (gencost(:, MODEL) ~= POLYNOMIAL) )
  error('copf.m: Unrecognized generator cost model in input data');
end

% print warning if there are dispatchable loads
if any(isload(gen))
  fprintf('\ncopf.m: Warning: Found dispatchable load in data (see ''help isload''\n');
  fprintf('        for details on what consitutes a dispatchable load). Only\n');
  fprintf('        the generalized formulation algorithms can maintain the power\n');
  fprintf('        factor of dispatchable loads constant; the Q injection will\n');
  fprintf('        be treated as a free variable.\n');
end

if mpopt(11) == 0 % OPF_ALG not set
  if any( gencost(:, MODEL) == PW_LINEAR )
    mpopt(11) = 200; % CCV formulation, constr
  else  % All are polynomial
    mpopt(11) = 100; % Standard, constr
  end
end

alg = mpopt(11);
formulation = opf_form(alg);
if (alg ~= 100) & (alg ~= 200)
  error('copf: constr.m solver called, but another solver specified in options');
end
if (alg == 200) & any(gencost(:, MODEL) == POLYNOMIAL)  % mixed models
  if verbose
    fprintf('copf.m: some generators use poly cost model, all will be converted\n');
    fprintf('        to piecewise linear.\n');
  end
  [pcost, qcost] = pqcost(gencost, size(gen, 1));
  i_poly = find(pcost(:, MODEL) == POLYNOMIAL);
  tmp = poly2pwl(pcost(i_poly, :), gen(i_poly, PMIN), ...
                     gen(i_poly, PMAX), npts);
  pcost(i_poly, 1:size(tmp,2)) = tmp;
  if ~isempty(qcost)
    i_poly = find(qcost(:, MODEL) == POLYNOMIAL);
    tmp = poly2pwl(qcost(i_poly, :), gen(i_poly, QMIN), ...
                       gen(i_poly, QMAX), npts);
    qcost(i_poly, 1:size(tmp,2)) = tmp;
  end
  gencost = [pcost; qcost];
end

if (alg == 100) & any(gencost(:, MODEL) == PW_LINEAR)
  error('copf.m: Standard formulation requested, but there are piece-wise linear costs in data');
end

% Now go for it.
gbus = gen(:, GEN_BUS);        %% what buses are gens at?

%% sizes of things
nb = size(bus, 1);
nl = size(branch, 1);
npv  = length(pv);
npq  = length(pq);
ng = size(gen, 1);            %% number of generators that are turned on

%% initial state
V  = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
V(gbus) = gen(:, VG) ./ abs(V(gbus)).* V(gbus);
Pg  = gen(:, PG) / baseMVA;
Qg  = gen(:, QG) / baseMVA;

%% check for costs for Qg
[pcost, qcost] = pqcost(gencost, size(gen, 1));

%% set up indexing for x
j1 = 1;      j2  = npv;        %% j1:j2  - V angle of pv buses
j3 = j2 + 1;  j4  = j2 + npq;      %% j3:j4  - V angle of pq buses
j5 = j4 + 1;  j6  = j4 + nb;      %% j5:j6  - V mag of all buses
j7 = j6 + 1;  j8  = j6 + ng;      %% j7:j8  - P of generators
j9 = j8 + 1;  j10  = j8 + ng;      %% j9:j10  - Q of generators
j11 = j10 + 1;  j12  = j10 + ng;      %% j11:j12  - Cp, cost of Pg
j13 = j12 + 1;  j14  = j12 + ng;      %% j13:j14  - Cq, cost of Qg

%% set up x
if formulation == 1
  Cp = [];
  Cq = [];
else
  Cp = totcost(pcost, Pg * baseMVA);
  Cq = totcost(qcost, Qg * baseMVA);  %% empty if qcost is empty
end
x = [angle(V([pv; pq])); abs(V); Pg; Qg; Cp; Cq];

%% objective and constraint function names
[fun, grad]  = fg_names(alg);
mpopt(15)  = 2 * nb;        %% set number of equality constraints

%% set some options
if mpopt(19) == 0
  mpopt(19) = 2 * nb + 150;  %% set max number of iterations for constr
end

%% set up options for Optim Tbx's constr
otopt = foptions;        %% get default options for constr
otopt(1) = (verbose > 0);    %% set verbose flag appropriately
% otopt(9) = 1;          %% check user supplied gradients?
otopt(2)  = mpopt(17);      %% termination tolerance on 'x'
otopt(3)  = mpopt(18);      %% termination tolerance on 'F'
otopt(4)  = mpopt(16);      %% termination tolerance on constraint violation
otopt(13) = mpopt(15);      %% number of equality constraints
otopt(14) = mpopt(19);      %% maximum number of iterations

%% run optimization
[x, otopt, lambda] = constr(fun, x, otopt, [], [], grad, ...
       baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
%% get final objective function value & constraint values
[f, g] = feval(fun, x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ...
               ref, pv, pq, mpopt);

%% check for convergence
if otopt(10) >= otopt(14)  | max(abs(g(1:otopt(13))))      > otopt(4) ...
              | max(g((otopt(13)+1):length(g)))  > otopt(4)
  success = 0;        %% did NOT converge
else
  success = 1;        %% DID converge
end
info = success;

%% reconstruct V
Va = zeros(nb, 1);
Va([ref; pv; pq]) = [angle(V(ref)); x(j1:j2); x(j3:j4)];
Vm = x(j5:j6);
V = Vm .* exp(j * Va);

%% grab Pg & Qg
Sg = x(j7:j8) + j * x(j9:j10);    %% complex power generation in p.u.

%%-----  calculate return values  -----
%% update bus, gen, branch with solution info
[bus, gen, branch] = opfsoln(baseMVA, bus, gen, branch, ...
            Ybus, Yf, Yt, V, Sg, lambda, ref, pv, pq, mpopt);

if nargout > 8  % Do we need to provide the Jacobian?
  [df, dg] = feval(grad, x, baseMVA, bus, gen, gencost, branch, Ybus, ...
                   Yf, Yt, V, ref, pv, pq, mpopt);
  jac = dg';
end

% convert to original external bus ordering
[bus, gen, branch, areas] = int2ext(i2e, bus, gen, branch, areas);

% Now create output matrices with all lines, all generators, committed and
% non-committed
genout = genorg;
branchout = branchorg;
genout(comgen, : ) = gen;
branchout(onbranch, :)  = branch;
% And zero out appropriate fields of non-comitted generators and lines
tmp = zeros(length(offgen), 1);
genout(offgen, PG) = tmp;
genout(offgen, QG) = tmp;
genout(offgen, MU_PMAX) = tmp;
genout(offgen, MU_PMIN) = tmp;
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;

%% compute elapsed time
et = etime(clock, t1);
if  (nargout == 0) & ( info > 0)
  printpf(baseMVA, bus, genout, branchout, f, info, et, 1, mpopt);
end

if nargout, busout = bus; end

return;