gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/grad_std.m
function [df, dg, d2f] = OTSgra(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt) %OTSGRA Evaluates gradients of objective function & constraints for OPF. % [df, dg] = OTSgra(x, baseMVA, bus, gen, gencost, branch, Ybus, % Yf, Yt, V, ref, pv, pq, mpopt) % Also, if a third output argument is specified, it will compute the 2nd % derivative matrix for the objective function. % 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 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; %% 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 partials of objective function ----- %% generator info on = find(gen(:, GEN_STATUS)); %% which generators are on? gbus = gen(on, GEN_BUS); %% what buses are they at? %% compute values of objective function partials [pcost, qcost] = pqcost(gencost, size(gen, 1), on); df_dPg = zeros(ng, 1); df_dQg = zeros(ng, 1); for i = 1:ng df_dPg(i) = polyval(polyder( pcost( i, COST:(COST+pcost(i, N)-1) ) ), Pg(i)*baseMVA) * baseMVA; %% w.r.t p.u. Pg end if ~isempty(qcost) %% Qg is not free for i = 1:ng df_dQg(i) = polyval(polyder( qcost( i, COST:(COST+qcost(i, N)-1) ) ), Qg(i)*baseMVA) * baseMVA; %% w.r.t p.u. Qg end end df = [ zeros(j6, 1); %% partial w.r.t. Va & Vm df_dPg; %% partial w.r.t. Pg df_dQg ]; %% partial w.r.t. Qg %%----- 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); %% evaluate partials of constraints dg = [ %% equality constraints real(dSbus_dVa(:,[pv;pq])), real(dSbus_dVm), ... real(dSbus_dPg), real(dSbus_dQg); %% P mismatch imag(dSbus_dVa(:,[pv;pq])), imag(dSbus_dVm), ... imag(dSbus_dPg), imag(dSbus_dQg); %% Q mismatch %% inequality constraints (variable limits, voltage & real generation) sparse(nb,j4), -speye(nb,nb), sparse(nb,2*ng); %% Vmin for var V sparse(nb,j4), speye(nb,nb), sparse(nb,2*ng); %% Vmax for var V sparse(ng,j6), -speye(ng,ng), sparse(ng,ng); %% Pmin for generators sparse(ng,j6), speye(ng,ng), sparse(ng,ng); %% Pmax for generators sparse(ng,j8), -speye(ng,ng); %% Qmin for generators sparse(ng,j8), speye(ng,ng); %% Qmax for generators %% inequality constraints (reactive generation & line flow limits) dAf_dVa(:,[pv;pq]), dAf_dVm, sparse(nl,2*ng); %% |Sf| limit dAt_dVa(:,[pv;pq]), dAt_dVm, sparse(nl,2*ng); %% |St| limit ]'; %% make full so optimization toolbox doesn't go wacky dg = full(dg); %% compute 2nd derivative of cost if nargin > 2 d2f_dPg2 = zeros(ng, 1); d2f_dQg2 = zeros(ng, 1); for i = 1:ng d2f_dPg2(i) = polyval(polyder(polyder( pcost( i, COST:(COST+pcost(i, N)-1) ) )), Pg(i)*baseMVA) * baseMVA^2; %% w.r.t p.u. Pg end if ~isempty(qcost) %% Qg is not free for i = 1:ng d2f_dQg2(i) = polyval(polyder(polyder( qcost( i, COST:(COST+qcost(i, N)-1) ) )), Qg(i)*baseMVA) * baseMVA^2; %% w.r.t p.u. Qg end end i = [j7:j10]'; d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2]); end end return;