gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/fun_std.m
function [f, g] = fun_std(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt) %FUN_STD Evaluates objective function & constraints for OPF. % [f, g] = fun_std(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; %% constant j = sqrt(-1); %% 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 %% 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 %% grab Pg & Qg Pg = x(j7:j8); %% active generation in p.u. Qg = x(j9:j10); %% reactive generation in p.u. %%----- evaluate objective function ----- %% generator info on = find(gen(:, GEN_STATUS)); %% which generators are on? gbus = gen(on, GEN_BUS); %% what buses are they at? %% 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 [pcost, qcost] = pqcost(gencost, size(gen, 1), on); f = sum( [totcost(pcost, gen(on, PG)); ... %% cost of Pg totcost(qcost, gen(on, QG)) ] ); %% cost of Qg, empty if no qcost %%----- 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 injected at "from" bus (p.u.) St = V(branch(:, T_BUS)) .* conj(Yt * V); %% complex power injected at "to" bus (p.u.) %% 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) ]; end return;