gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/fun_ccv.m
function [f, g] = fun_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt) %FUN_CCV Evaluates objective function & constraints for OPF. % [f, g] = fun_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt) % MATPOWER Version 2.0 % by Ray Zimmerman, PSERC Cornell 12/12/97 % Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC) % See http://www.pserc.cornell.edu/ for more info. %%----- initialize ----- %% define named indices into gen, branch 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; [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; [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost; %% options alg = mpopt(11); %% constant j = sqrt(-1); %% generator info on = find(gen(:, GEN_STATUS)); %% which generators are on? gbus = gen(on, GEN_BUS); %% what buses are they at? %% 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 %% check for costs for Qg [pcost, qcost] = pqcost(gencost, size(gen, 1), on); if isempty(qcost) %% set number of cost variables ncv = ng; %% only Cp else ncv = 2 * ng; %% Cp and Cq end %% 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 + npv + 1; %% j11:j12 - Cp, cost of Pg %% grab Pg & Qg and their costs Pg = x(j7:j8); %% active generation in p.u. Qg = x(j9:j10); %% reactive generation in p.u. Cp = x(j11:j12); %% active costs in $/hr if ncv > ng %% no free VARs j13 = j12 + 1; j14 = j12 + npv + 1; %% j13:j14 - Cq, cost of Qg Cq = x(j13:j14); %% reactive costs in $/hr end %%----- evaluate objective function ----- %% put Pg & Qg back in gen gen(on, PG) = Pg * baseMVA; %% active generation in MW gen(on, QG) = Qg * baseMVA; %% reactive generation in MVAR %% compute objective value if ncv > ng %% no free VARs f = sum([Cp; Cq]); else %% free VARs f = sum(Cp); end %%----- evaluate constraints ----- if nargout > 1 %% 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); %% rebuild Sbus Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u. %% evaluate power flow equations mis = V .* conj(Ybus * V) - Sbus; %% compute branch power flows Sf = V(branch(:, F_BUS)) .* conj(Yf * V); %% complex power at "from" bus (p.u.) St = V(branch(:, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) %% compute generator cost constraints ( costfcn @ Pg - Cp , etc.) Qcc = []; nsegs = pcost(:, N) - 1; %% number of cost constraints for each gen ncc = sum(nsegs); %% total number of cost constraints Pcc = zeros(ncc, 1); for i = 1:ng xx = pcost(i, COST:2:( COST + 2*(nsegs(i)) ))'; yy = pcost(i, (COST+1):2:( COST + 2*(nsegs(i)) + 1))'; i1 = 1:nsegs(i); i2 = 2:(nsegs(i) + 1); m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1)); b = yy(i1) - m .* xx(i1); Pcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ... m .* gen(on(i), PG) + b - Cp(i); end if ncv > ng %% no free VARs nsegs = qcost(:, N) - 1; %% number of cost constraints for each gen ncc = sum(nsegs); %% total number of cost constraints Qcc = zeros(ncc, 1); for i = 1:ng xx = qcost(i, COST:2:( COST + 2*(nsegs(i)) ))'; yy = qcost(i, (COST+1):2:( COST + 2*(nsegs(i)) + 1))'; i1 = 1:nsegs(i); i2 = 2:(nsegs(i) + 1); m = (yy(i2) - yy(i1)) ./ (xx(i2) - xx(i1)); b = yy(i1) - m .* xx(i1); Qcc(sum(nsegs(1:(i-1))) + [1:nsegs(i)]) = ... m .* gen(on(i), QG) + b - Cq(i); end end %% compute constraint function values g = [ %% equality constraints real(mis); %% active power mismatch for all buses imag(mis); %% reactive power mismatch for all buses %% inequality constraints (variable limits, voltage & generation) bus(:, VMIN) - x(j5:j6); %% lower voltage limit for var V x(j5:j6) - bus(:, VMAX); %% upper voltage limit for var V gen(on, PMIN) / baseMVA - x(j7:j8); %% lower generator P limit x(j7:j8) - gen(on, PMAX) / baseMVA; %% upper generator P limit gen(on, QMIN) / baseMVA - x(j9:j10); %% lower generator Q limit x(j9:j10) - gen(on, QMAX) / baseMVA; %% upper generator Q limit %% inequality constraints (line flow limits) abs(Sf) - branch(:, RATE_A) / baseMVA; %% branch apparent power limits (from bus) abs(St) - branch(:, RATE_A) / baseMVA; %% branch apparent power limits (to bus) %% inequality cost constraints Pcc; Qcc; ]; end return;