www.gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/opf.m
function [bus, gen, branch, f, success, et] = opf(baseMVA, bus, gen, gencost, ...
branch, Ybus, Yf, Yt, ref, pv, pq, mpopt)
%OPF Solves an optimal power flow.
% [bus, gen, branch, f, success, et] = opf(baseMVA, bus, gen, gencost, ...
% branch, Ybus, Yf, Yt, ref, pv, pq, mpopt)
% Assumes that Ybus, Yf, Yt, V, ref, pv, pq are consistent with
% (or override) data in bus, gen, branch.
% MATPOWER Version 2.0
% by Ray Zimmerman & David Gan, PSERC Cornell 12/24/97
% Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
% See http://www.pserc.cornell.edu/ for more info.
%%----- initialization -----
%% default arguments
if nargin < 12
mpopt = mpoption; %% use default options
end
%% options
verbose = mpopt(31);
npts = mpopt(14); %% number of points to evaluate when converting
%% polynomials to piece-wise linear
%% define constants
j = sqrt(-1);
%% define named indices into data matrices
[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] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost;
%%----- check/convert costs, set default algorithm -----
%% get cost model, check consistency
model = gencost(:, MODEL);
i_pwln = find(model == PW_LINEAR);
i_poly = find(model == POLYNOMIAL);
if any(i_pwln) & any(i_poly) & verbose
fprintf('not all generators use same cost model, all will be converted to piece-wise linear\n');
end
%% set algorithm
if mpopt(11) == 0
%% use default for this cost model
if find(model ~= PW_LINEAR & model ~= POLYNOMIAL)
error('unknown generator cost model');
elseif any(i_pwln) %% some piece-wise linear, use appropriate alg
mpopt(11) = mpopt(13);
else %% must all be polynomial
mpopt(11) = mpopt(12);
end
end
alg = mpopt(11);
formulation = opf_form(alg);
%% move Pmin and Pmax limits out slightly to avoid problems
%% caused by rounding errors when corner point of cost
%% function lies at exactly Pmin or Pmax
if any(i_pwln)
ng = size(gen, 1);
gen(:, PMIN) = gen(:, PMIN) - 1e-6 * ones(ng, 1);
gen(:, PMAX) = gen(:, PMAX) + 1e-6 * ones(ng, 1);
end
%% check cost model/algorithm consistency
if any( i_pwln ) & formulation == 1
error(sprintf('algorithm %d does not handle piece-wise linear cost functions', alg));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Eventually put code here to fit polynomials to piece-wise linear as needed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% convert polynomials to piece-wise linear
if any(i_poly) & formulation == 2
if verbose
fprintf('converting from polynomial to piece-wise linear cost model\n');
end
[pcost, qcost] = pqcost(gencost, size(gen, 1));
i_poly = find(pcost(:, MODEL) == POLYNOMIAL);
pcost = poly2pwl(pcost(i_poly, :), gen(i_poly, PMIN), gen(i_poly, PMAX), npts);
if ~isempty(qcost)
i_poly = find(qcost(:, MODEL) == POLYNOMIAL);
qcost = poly2pwl(qcost(i_poly, :), gen(i_poly, QMIN), gen(i_poly, QMAX), npts);
end
gencost = [pcost; qcost];
end
%%----- run opf -----
%% start timer
t0 = clock;
%% sizes of things
nb = size(bus, 1);
nl = size(branch, 1);
npv = length(pv);
npq = length(pq);
ng = npv + 1; %% number of generators that are turned on
%% initial state
on = find(gen(:, GEN_STATUS)); %% which generators are on?
V = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
V(gen(on, GEN_BUS)) = gen(on, VG) ./ abs(V(gen(on, GEN_BUS))).* V(gen(on, GEN_BUS));
Pg = gen(on, PG) / baseMVA;
Qg = gen(on, QG) / baseMVA;
%% check for costs for Qg
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);
%% 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 + npv + 1; %% j7:j8 - P of generators
j9 = j8 + 1; j10 = j8 + npv + 1; %% j9:j10 - Q of generators
% j11 = j10 + 1; j12 = j10 + npv + 1; %% j11:j12 - Cp, cost of Pg
% j13 = j12 + 1; j14 = j12 + npv + 1; %% 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];
%% run constrained optimization
[fun, grad] = fg_names(alg);
mpopt(15) = 2 * nb; %% set number of equality constraints
if opf_slvr(alg) == 0 %% use CONSTR
%% 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(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
else %% use LPCONSTR
%% 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);
%% get final objective function value & constraint values
f = feval(fun, x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
end
%%----- calculate return values -----
%% 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.
%% 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);
%% compute elapsed time
et = etime(clock, t0);
return;