opf_power_balance_fcn

OPF_POWER_BALANCE_FCN  Evaluates AC power balance constraints and their gradients.
   [G, DG] = OPF_POWER_BALANCE_FCN(X, OM, YBUS, MPOPT)

   Computes the active or reactive power balance equality constraints for
   AC optimal power flow. Computes constraint vectors and their gradients.

   Inputs:
     X : optimization vector
     MPC : MATPOWER case struct
     YBUS : bus admittance matrix
     MPOPT : MATPOWER options struct

   Outputs:
     G  : vector of equality constraint values (active/reactive power balances)
     DG : (optional) equality constraint gradients

   Examples:
       g = opf_power_balance_fcn(x, mpc, Ybus, mpopt);
       [g, dg] = opf_power_balance_fcn(x, mpc, Ybus, mpopt);

   See also OPF_POWER_BALANCE_HESS

0001 function [g, dg] = opf_power_balance_fcn(x, mpc, Ybus, mpopt)
0002 %OPF_POWER_BALANCE_FCN  Evaluates AC power balance constraints and their gradients.
0003 %   [G, DG] = OPF_POWER_BALANCE_FCN(X, OM, YBUS, MPOPT)
0004 %
0005 %   Computes the active or reactive power balance equality constraints for
0006 %   AC optimal power flow. Computes constraint vectors and their gradients.
0007 %
0008 %   Inputs:
0009 %     X : optimization vector
0010 %     MPC : MATPOWER case struct
0011 %     YBUS : bus admittance matrix
0012 %     MPOPT : MATPOWER options struct
0013 %
0014 %   Outputs:
0015 %     G  : vector of equality constraint values (active/reactive power balances)
0016 %     DG : (optional) equality constraint gradients
0017 %
0018 %   Examples:
0019 %       g = opf_power_balance_fcn(x, mpc, Ybus, mpopt);
0020 %       [g, dg] = opf_power_balance_fcn(x, mpc, Ybus, mpopt);
0021 %
0022 %   See also OPF_POWER_BALANCE_HESS
0023 
0024 %   MATPOWER
0025 %   Copyright (c) 1996-2017, Power Systems Engineering Research Center (PSERC)
0026 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0027 %   and Ray Zimmerman, PSERC Cornell
0028 %   and Baljinnyam Sereeter, Delft University of Technology
0029 %
0030 %   This file is part of MATPOWER.
0031 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0032 %   See https://matpower.org for more info.
0033 
0034 %%----- initialize -----
0035 %% define named indices into data matrices
0036 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0037     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0038     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0039 
0040 %% unpack data
0041 [baseMVA, bus, gen] = deal(mpc.baseMVA, mpc.bus, mpc.gen);
0042 if mpopt.opf.v_cartesian
0043     [Vr, Vi, Pg, Qg] = deal(x{:});
0044     V = Vr + 1j * Vi;           %% reconstruct V
0045 else
0046     [Va, Vm, Pg, Qg] = deal(x{:});
0047     V = Vm .* exp(1j * Va);     %% reconstruct V
0048 end
0049 
0050 %% problem dimensions
0051 nb = length(V);             %% number of buses
0052 ng = length(Pg);            %% number of dispatchable injections
0053 
0054 %% ----- evaluate constraints -----
0055 %% put Pg, Qg back in gen
0056 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0057 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0058 
0059 %% rebuild Sbus
0060 if mpopt.opf.v_cartesian
0061     Sbus = makeSbus(baseMVA, bus, gen);             %% net injected power in p.u.
0062 else
0063     Sbus = makeSbus(baseMVA, bus, gen, mpopt, Vm);  %% net injected power in p.u.
0064 end
0065 
0066 %% evaluate complex power balance mismatches
0067 mis = V .* conj(Ybus * V) - Sbus;
0068 
0069 %% assemble active and reactive power balance constraints
0070 g = [ real(mis);    %% active power mismatch
0071       imag(mis) ];  %% reactive power mismatch
0072 
0073 %%----- evaluate constraint gradients -----
0074 if nargout > 1
0075     %% compute partials of injected bus powers
0076     [dSbus_dV1, dSbus_dV2] = dSbus_dV(Ybus, V, mpopt.opf.v_cartesian);  %% w.r.t. V
0077     neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);     %% Pbus w.r.t. Pg
0078     if ~mpopt.opf.v_cartesian
0079         %% adjust for voltage dependent loads
0080         [dummy, neg_dSd_dVm] = makeSbus(baseMVA, bus, gen, mpopt, Vm);
0081         dSbus_dV2 = dSbus_dV2 - neg_dSd_dVm;
0082     end
0083     dg = [
0084         real([dSbus_dV1 dSbus_dV2]) neg_Cg sparse(nb, ng);  %% P mismatch w.r.t V1, V2, Pg, Qg
0085         imag([dSbus_dV1 dSbus_dV2]) sparse(nb, ng) neg_Cg;  %% Q mismatch w.r.t V1, V2, Pg, Qg
0086     ];
0087 end