gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/grad_ccv.m
function [df, dg] = grad_ccv(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt) %GRAD_CCV Evaluates gradients of objective function & constraints for OPF. % [df, dg] = grad_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 data, branch matrices [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 partials of objective function ----- %% compute values of objective function partials df = [ zeros(j10, 1); %% partial w.r.t. Va, Vm, Pg, Qg ones(ncv, 1) ]; %% partial w.r.t. Cp (and Cq) %%----- evaluate partials of 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); %% compute partials of injected bus powers [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); %% w.r.t. V dSbus_dPg = sparse(gbus, 1:ng, -1, nb, ng); %% w.r.t. Pg dSbus_dQg = sparse(gbus, 1:ng, -j, nb, ng); %% w.r.t. Qg %% compute partials of line flows w.r.t. V [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V); %% line limits are w.r.t apparent power, so compute partials of apparent power [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ... dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St); %% cost constraints w.r.t everything ( d(costfcn @ Pg - Cp) , etc.) dQcc_dQg = zeros(0, ng); dQcc_dCq = []; nsegs = pcost(:, N) - 1; %% number of cost constraints for each gen nPcc = sum(nsegs); %% total number of cost constraints dPcc_dPg = sparse([], [], [], nPcc, ng, nPcc); %% nPcc x ng dPcc_dCp = sparse([], [], [], nPcc, ng, nPcc); %% nPcc x ng 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)); ii = sum(nsegs(1:(i-1))) + [1:nsegs(i)]; dPcc_dPg(ii, i) = m * baseMVA; dPcc_dCp(ii, i) = -1 * ones(nsegs(i), 1); end nQcc = 0; if ncv > ng %% no free VARs nsegs = qcost(:, N) - 1; %% number of cost constraints for each gen nQcc = sum(nsegs); %% total number of cost constraints dQcc_dQg = sparse([], [], [], nQcc, ng, nQcc); %% nQcc x ng dQcc_dCq = sparse([], [], [], nQcc, ng, nQcc); %% nQcc x ng 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)); ii = sum(nsegs(1:(i-1))) + [1:nsegs(i)]; dQcc_dQg(ii, i) = m * baseMVA; dQcc_dCq(ii, i) = -1 * ones(nsegs(i), 1); end end %% [ dcc_dV dcc_dPg dcc_dQg dcc_dCp dcc_dCq ] dPcc = [sparse(nPcc,j6), dPcc_dPg, sparse(nPcc,ng), dPcc_dCp, sparse(nPcc,ncv-ng)]; dQcc = [sparse(nQcc,j8), dQcc_dQg, sparse(nQcc,ng), dQcc_dCq]; %% evaluate partials of constraints dg = [ %% equality constraints real(dSbus_dVa(:,[pv;pq])), real(dSbus_dVm), ... real(dSbus_dPg), real(dSbus_dQg), sparse(nb,ncv); %% P mismatch imag(dSbus_dVa(:,[pv;pq])), imag(dSbus_dVm), ... imag(dSbus_dPg), imag(dSbus_dQg), sparse(nb,ncv); %% Q mismatch %% inequality constraints (variable limits, voltage & real generation) sparse(nb,j4), -speye(nb,nb), sparse(nb,2*ng+ncv); %% Vmin for var V sparse(nb,j4), speye(nb,nb), sparse(nb,2*ng+ncv); %% Vmax for var V sparse(ng,j6), -speye(ng,ng), sparse(ng,ng+ncv); %% Pmin for generators sparse(ng,j6), speye(ng,ng), sparse(ng,ng+ncv); %% Pmax for generators sparse(ng,j8), -speye(ng,ng), sparse(ng,ncv); %% Qmin for generators sparse(ng,j8), speye(ng,ng), sparse(ng,ncv); %% Qmax for generators %% inequality constraints (reactive generation & line flow limits) dAf_dVa(:,[pv;pq]), dAf_dVm, sparse(nl,2*ng+ncv); %% |Sf| limit dAt_dVa(:,[pv;pq]), dAt_dVm, sparse(nl,2*ng+ncv); %% |St| limit %% inequality constraints on generator cost dPcc; dQcc; ]'; %% make full so optimization toolbox doesn't go wacky dg = full(dg); end return;