Home > matpower5.1 > opf_consfcn.m

opf_consfcn

PURPOSE ^

OPF_CONSFCN Evaluates nonlinear constraints and their Jacobian for OPF.

SYNOPSIS ^

function [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il, varargin)

DESCRIPTION ^

OPF_CONSFCN  Evaluates nonlinear constraints and their Jacobian for OPF.
   [H, G, DH, DG] = OPF_CONSFCN(X, OM, YBUS, YF, YT, MPOPT, IL)

   Constraint evaluation function for AC optimal power flow, suitable
   for use with MIPS or FMINCON. Computes constraint vectors and their
   gradients.

   Inputs:
     X : optimization vector
     OM : OPF model object
     YBUS : bus admittance matrix
     YF : admittance matrix for "from" end of constrained branches
     YT : admittance matrix for "to" end of constrained branches
     MPOPT : MATPOWER options struct
     IL : (optional) vector of branch indices corresponding to
          branches with flow limits (all others are assumed to be
          unconstrained). The default is [1:nl] (all branches).
          YF and YT contain only the rows corresponding to IL.

   Outputs:
     H  : vector of inequality constraint values (flow limits)
          limit^2 - flow^2, where the flow can be apparent power
          real power or current, depending on value of
          opf.flow_lim in MPOPT (only for constrained lines)
     G  : vector of equality constraint values (power balances)
     DH : (optional) inequality constraint gradients, column j is
          gradient of H(j)
     DG : (optional) equality constraint gradients

   Examples:
       [h, g] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt);
       [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt);
       [h, g, dh, dg] = opf_consfcn(x, om, Ybus, Yf, Yt, mpopt, il);

   See also OPF_COSTFCN, OPF_HESSFCN.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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-2015 by Power System Engineering Research Center (PSERC)
0040 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   and Ray Zimmerman, PSERC Cornell
0042 %
0043 %   $Id: opf_consfcn.m 2644 2015-03-11 19:34:22Z ray $
0044 %
0045 %   This file is part of MATPOWER.
0046 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0047 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0048 
0049 %%----- initialize -----
0050 %% define named indices into data matrices
0051 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0052     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0053     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0054 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0055     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0056     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0057 
0058 %% unpack data
0059 mpc = get_mpc(om);
0060 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0061 vv = get_idx(om);
0062 
0063 %% problem dimensions
0064 nb = size(bus, 1);          %% number of buses
0065 nl = size(branch, 1);       %% number of branches
0066 ng = size(gen, 1);          %% number of dispatchable injections
0067 nxyz = length(x);           %% total number of control vars of all types
0068 
0069 %% set default constrained lines
0070 if nargin < 7
0071     il = (1:nl);            %% all lines have limits by default
0072 end
0073 nl2 = length(il);           %% number of constrained lines
0074 
0075 %% grab Pg & Qg
0076 Pg = x(vv.i1.Pg:vv.iN.Pg);  %% active generation in p.u.
0077 Qg = x(vv.i1.Qg:vv.iN.Qg);  %% reactive generation in p.u.
0078 
0079 %% put Pg & Qg back in gen
0080 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0081 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0082  
0083 %% rebuild Sbus
0084 Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u.
0085 
0086 %% ----- evaluate constraints -----
0087 %% reconstruct V
0088 Va = x(vv.i1.Va:vv.iN.Va);
0089 Vm = x(vv.i1.Vm:vv.iN.Vm);
0090 V = Vm .* exp(1j * Va);
0091 
0092 %% evaluate power flow equations
0093 mis = V .* conj(Ybus * V) - Sbus;
0094 
0095 %%----- evaluate constraint function values -----
0096 %% first, the equality constraints (power flow)
0097 g = [ real(mis);            %% active power mismatch for all buses
0098       imag(mis) ];          %% reactive power mismatch for all buses
0099 
0100 %% then, the inequality constraints (branch flow limits)
0101 if nl2 > 0
0102   flow_max = (branch(il, RATE_A)/baseMVA).^2;
0103   flow_max(flow_max == 0) = Inf;
0104   if upper(mpopt.opf.flow_lim(1)) == 'I'    %% current magnitude limit, |I|
0105     If = Yf * V;
0106     It = Yt * V;
0107     h = [ If .* conj(If) - flow_max;    %% branch current limits (from bus)
0108           It .* conj(It) - flow_max ];  %% branch current limits (to bus)
0109   else
0110     %% compute branch power flows
0111     Sf = V(branch(il, F_BUS)) .* conj(Yf * V);  %% complex power injected at "from" bus (p.u.)
0112     St = V(branch(il, T_BUS)) .* conj(Yt * V);  %% complex power injected at "to" bus (p.u.)
0113     if upper(mpopt.opf.flow_lim(1)) == 'P'  %% active power limit, P (Pan Wei)
0114       h = [ real(Sf).^2 - flow_max;         %% branch real power limits (from bus)
0115             real(St).^2 - flow_max ];       %% branch real power limits (to bus)
0116     else                                    %% apparent power limit, |S|
0117       h = [ Sf .* conj(Sf) - flow_max;      %% branch apparent power limits (from bus)
0118             St .* conj(St) - flow_max ];    %% branch apparent power limits (to bus)
0119     end
0120   end
0121 else
0122   h = zeros(0,1);
0123 end
0124 
0125 %%----- evaluate partials of constraints -----
0126 if nargout > 2
0127   %% index ranges
0128   iVa = vv.i1.Va:vv.iN.Va;
0129   iVm = vv.i1.Vm:vv.iN.Vm;
0130   iPg = vv.i1.Pg:vv.iN.Pg;
0131   iQg = vv.i1.Qg:vv.iN.Qg;
0132 
0133   %% compute partials of injected bus powers
0134   [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);           %% w.r.t. V
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

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005