Description of opf

0001 function [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il, varargin)
0002 %OPF_CONSFCN  Evaluates nonlinear constraints and their Jacobian for OPF.
0003 %   [H, G, DH, DG] = OPF_CONSFCN(X, OM, YBUS, YF, YT, MPOPT, IL)
0004 %
0005 %   Constraint evaluation function for AC optimal power flow, suitable
0006 %   for use with MIPS or FMINCON. Computes constraint vectors and their
0007 %   gradients.
0008 %
0009 %   Inputs:
0010 %     X : optimization vector
0011 %     OM : OPF model object
0012 %     YBUS : bus admittance matrix
0013 %     YF : admittance matrix for "from" end of constrained branches
0014 %     YT : admittance matrix for "to" end of constrained branches
0015 %     MPOPT : MATPOWER options struct
0016 %     IL : (optional) vector of branch indices corresponding to
0017 %          branches with flow limits (all others are assumed to be
0018 %          unconstrained). The default is [1:nl] (all branches).
0019 %          YF and YT contain only the rows corresponding to IL.
0020 %
0021 %   Outputs:
0022 %     H  : vector of inequality constraint values (flow limits)
0023 %          limit^2 - flow^2, where the flow can be apparent power
0024 %          real power or current, depending on value of
0025 %          opf.flow_lim in MPOPT (only for constrained lines)
0026 %     G  : vector of equality constraint values (power balances)
0027 %     DH : (optional) inequality constraint gradients, column j is
0028 %          gradient of H(j)
0029 %     DG : (optional) equality constraint gradients
0030 %
0031 %   Examples:
0032 %       [h, g] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt);
0033 %       [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt);
0034 %       [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il);
0035 %
0036 %   See also OPF_COSTFCN, OPF_HESSFCN.
0037 
0038 %   MATPOWER
0039 %   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
0040 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0041 %   Ray Zimmerman, PSERC Cornell and Shrirang Abhyankar
0042 %
0043 %   This file is part of MATPOWER.
0044 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0045 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0046 
0047 %%----- initialize -----
0048 %% define named indices into data matrices
0049 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0050     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0051     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0052 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0053     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0054     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0055 
0056 %% unpack data
0057 mpc = get_mpc(om);
0058 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0059 vv = get_idx(om);
0060 
0061 %% problem dimensions
0062 nb = size(bus, 1);          %% number of buses
0063 nl = size(branch, 1);       %% number of branches
0064 ng = size(gen, 1);          %% number of dispatchable injections
0065 nxyz = length(x);           %% total number of control vars of all types
0066 
0067 %% set default constrained lines
0068 if nargin < 7
0069     il = (1:nl);            %% all lines have limits by default
0070 end
0071 nl2 = length(il);           %% number of constrained lines
0072 
0073 %% grab Pg & Qg
0074 Pg = x(vv.i1.Pg:vv.iN.Pg);  %% active generation in p.u.
0075 Qg = x(vv.i1.Qg:vv.iN.Qg);  %% reactive generation in p.u.
0076 
0077 %% put Pg & Qg back in gen
0078 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0079 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0080  
0081 %% ----- evaluate constraints -----
0082 %% reconstruct V
0083 Va = x(vv.i1.Va:vv.iN.Va);
0084 Vm = x(vv.i1.Vm:vv.iN.Vm);
0085 V = Vm .* exp(1j * Va);
0086 
0087 %% rebuild Sbus
0088 Sbus = makeSbus(baseMVA, bus, gen, mpopt, Vm);  %% net injected power in p.u.
0089 
0090 %% evaluate power flow equations
0091 mis = V .* conj(Ybus * V) - Sbus;
0092 
0093 %%----- evaluate constraint function values -----
0094 %% first, the equality constraints (power flow)
0095 g = [ real(mis);            %% active power mismatch for all buses
0096       imag(mis) ];          %% reactive power mismatch for all buses
0097 
0098 %% then, the inequality constraints (branch flow limits)
0099 if nl2 > 0
0100   flow_max = (branch(il, RATE_A)/baseMVA).^2;
0101   flow_max(flow_max == 0) = Inf;
0102   if upper(mpopt.opf.flow_lim(1)) == 'I'    %% current magnitude limit, |I|
0103     If = Yf * V;
0104     It = Yt * V;
0105     h = [ If .* conj(If) - flow_max;    %% branch current limits (from bus)
0106           It .* conj(It) - flow_max ];  %% branch current limits (to bus)
0107   else
0108     %% compute branch power flows
0109     Sf = V(branch(il, F_BUS)) .* conj(Yf * V);  %% complex power injected at "from" bus (p.u.)
0110     St = V(branch(il, T_BUS)) .* conj(Yt * V);  %% complex power injected at "to" bus (p.u.)
0111     if upper(mpopt.opf.flow_lim(1)) == 'P'  %% active power limit, P (Pan Wei)
0112       h = [ real(Sf).^2 - flow_max;         %% branch real power limits (from bus)
0113             real(St).^2 - flow_max ];       %% branch real power limits (to bus)
0114     else                                    %% apparent power limit, |S|
0115       h = [ Sf .* conj(Sf) - flow_max;      %% branch apparent power limits (from bus)
0116             St .* conj(St) - flow_max ];    %% branch apparent power limits (to bus)
0117     end
0118   end
0119 else
0120   h = zeros(0,1);
0121 end
0122 
0123 %%----- evaluate partials of constraints -----
0124 if nargout > 2
0125   %% index ranges
0126   iVa = vv.i1.Va:vv.iN.Va;
0127   iVm = vv.i1.Vm:vv.iN.Vm;
0128   iPg = vv.i1.Pg:vv.iN.Pg;
0129   iQg = vv.i1.Qg:vv.iN.Qg;
0130 
0131   %% compute partials of injected bus powers
0132   [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);           %% w.r.t. V
0133   [dummy, neg_dSd_dVm] = makeSbus(baseMVA, bus, gen, mpopt, Vm);
0134   dSbus_dVm = dSbus_dVm - neg_dSd_dVm;
0135   neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
0136                                                         %% Qbus w.r.t. Qg
0137   
0138   %% construct Jacobian of equality (power flow) constraints and transpose it
0139   dg = sparse(2*nb, nxyz);
0140   dg(:, [iVa iVm iPg iQg]) = [
0141     real([dSbus_dVa dSbus_dVm]) neg_Cg sparse(nb, ng);  %% P mismatch w.r.t Va, Vm, Pg, Qg
0142     imag([dSbus_dVa dSbus_dVm]) sparse(nb, ng) neg_Cg;  %% Q mismatch w.r.t Va, Vm, Pg, Qg
0143   ];
0144   dg = dg';
0145 
0146   if nl2 > 0
0147     %% compute partials of Flows w.r.t. V
0148     if upper(mpopt.opf.flow_lim(1)) == 'I'  %% current
0149       [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V);
0150     else                            %% power
0151       [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V);
0152     end
0153     if upper(mpopt.opf.flow_lim(1)) == 'P'  %% real part of flow (active power)
0154       dFf_dVa = real(dFf_dVa);
0155       dFf_dVm = real(dFf_dVm);
0156       dFt_dVa = real(dFt_dVa);
0157       dFt_dVm = real(dFt_dVm);
0158       Ff = real(Ff);
0159       Ft = real(Ft);
0160     end
0161   
0162     %% squared magnitude of flow (of complex power or current, or real power)
0163     [df_dVa, df_dVm, dt_dVa, dt_dVm] = ...
0164             dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft);
0165   
0166     %% construct Jacobian of inequality (branch flow) constraints & transpose
0167     dh = sparse(2*nl2, nxyz);
0168     dh(:, [iVa iVm]) = [
0169       df_dVa, df_dVm;                     %% "from" flow limit
0170       dt_dVa, dt_dVm;                     %% "to" flow limit
0171     ];
0172     dh = dh';
0173   else
0174     dh = sparse(nxyz, 0);
0175   end
0176 end
opf_consfcn

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
