grad_std.m

可进行电力系统多节点系统的优化潮流计算

function [df, dg, d2f] = grad_std(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%GRAD_STD  Evaluates gradients of objective function & constraints for OPF.
%   [df, dg] = grad_std(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
%   $Id: grad_std.m,v 1.8 2004/09/17 14:37:59 ray Exp $
%   by Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ 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);

%% generator info
on = find(gen(:, GEN_STATUS) > 0);      %% 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 = length(on);                        %% 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 -----
%% 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);
    br = find(branch(:, BR_STATUS));
    nbr = length(br);
    
    %% compute partials of line flow constraints
    if mpopt(24) == 1   %% branch active power limits
        dFlow_dV = [
            real(dSf_dVa(br,[pv;pq])), real(dSf_dVm(br,:)); %% from bus
            real(dSt_dVa(br,[pv;pq])), real(dSt_dVm(br,:)); %% to bus
        ];
    else                %% branch apparent power limits
        dFlow_dV = [
            dAf_dVa(br,[pv;pq]), dAf_dVm(br,:);             %% from bus
            dAt_dVa(br,[pv;pq]), dAt_dVm(br,:);             %% to bus
        ];
    end
    
    %% 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 & 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 (line flow limits)
        dFlow_dV, sparse(2*nbr,2*ng);
    ]';

    %% make full if using constr or non-sparse QP solver
    if opf_slvr(mpopt(11)) == 0 | ~have_fcn('sparse_qp') | mpopt(51) == 0
        dg = full(dg);
    end

    %% compute 2nd derivative of cost
    if nargout > 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;