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;