lpopf.m

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

function [busout, genout, branchout, f, success, info, et, g, jac] = lpopf( ...
                baseMVA, bus, gen, branch, areas, gencost, mpopt)
%LPOPF  Solves an AC optimal power flow using succesive LPs.
%
%   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] = lpopf(casefile, mpopt)
%
%   [bus, gen, branch, f, success] = lpopf(baseMVA, bus, gen, branch, areas, ...
%                                    gencost, mpopt)
%
%   [bus, gen, branch, f, success, info, et] = lpopf(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).

%   MATPOWER
%   $Id: lpopf.m,v 1.5 2005/01/21 20:43:21 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('lpopf.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('lpopf.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('\nlpopf.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) = 240; % CCV formulation, sparse LP (relaxed)
  else  % All are polynomial
    mpopt(11) = 140; % Standard formulation, sparse LP (relaxed)
  end
end

alg = mpopt(11);
formulation = opf_form(alg);
code = opf_slvr(alg);
if (code < 1) | (code > 3)
  error('lpopf: LP solver called, but another solver specified in options');
end
if (formulation == 2) & any(gencost(:, MODEL) == POLYNOMIAL)  % mixed models
  if verbose
    fprintf('lpopf.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 (formulation == 1) & 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

%% run load flow to get starting point
[x, success_lf] = LPeqslvr(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf,...
                           Yt, V, ref, pv, pq, mpopt);
if success_lf ~= 1
  error('Sorry, cannot find a starting point using power flow, please check data!');
end

%% set step size
cstep = 0;
if ~isempty(Cp)
  cstep = max(abs(Cp));
  if cstep < 1.0e6, cstep = 1.0e6; end
end
step0=[2*ones(nb-1,1);          %% starting stepsize for Vangle
    ones(nb,1);            %% Vmag
    0.6*ones(ng,1);          %% Pg
    0.3*ones(ng,1);          %% Qg
    cstep*ones(length(Cp),1);    %% Cp
    cstep*ones(length(Cq),1)];    %% Cq
idx_xi = [];

%% run optimization
[x, lambda, success] = LPconstr(fun, x, idx_xi, mpopt, step0, [], [], grad, ...
                       'LPeqslvr', baseMVA, bus, gen, gencost, branch, ...
                       Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
info = success;

%% get final objective function value & constraint values
f = feval(fun, x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref,...
          pv, pq, mpopt);

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

g = [];
jac = []; % we do not compute jacobians

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

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;