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 vector
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 %     H  : 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 %   $Id: opf_consfcn.m,v 1.7 2010/10/12 17:42:47 ray Exp $
0040 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
0041 %   and Ray Zimmerman, PSERC Cornell
0042 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0043 %
0044 %   This file is part of MATPOWER.
0045 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0046 %
0047 %   MATPOWER is free software: you can redistribute it and/or modify
0048 %   it under the terms of the GNU General Public License as published
0049 %   by the Free Software Foundation, either version 3 of the License,
0050 %   or (at your option) any later version.
0051 %
0052 %   MATPOWER is distributed in the hope that it will be useful,
0053 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0054 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0055 %   GNU General Public License for more details.
0056 %
0057 %   You should have received a copy of the GNU General Public License
0058 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0059 %
0060 %   Additional permission under GNU GPL version 3 section 7
0061 %
0062 %   If you modify MATPOWER, or any covered work, to interface with
0063 %   other modules (such as MATLAB code and MEX-files) available in a
0064 %   MATLAB(R) or comparable environment containing parts covered
0065 %   under other licensing terms, the licensors of MATPOWER grant
0066 %   you additional permission to convey the resulting work.
0067 
0068 %%----- initialize -----
0069 %% define named indices into data matrices
0070 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0071     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0072     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0073 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0074     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0075     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0076 
0077 %% unpack data
0078 mpc = get_mpc(om);
0079 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0080 vv = get_idx(om);
0081 
0082 %% problem dimensions
0083 nb = size(bus, 1);          %% number of buses
0084 nl = size(branch, 1);       %% number of branches
0085 ng = size(gen, 1);          %% number of dispatchable injections
0086 nxyz = length(x);           %% total number of control vars of all types
0087 
0088 %% set default constrained lines
0089 if nargin < 7
0090     il = (1:nl);            %% all lines have limits by default
0091 end
0092 nl2 = length(il);           %% number of constrained lines
0093 
0094 %% grab Pg & Qg
0095 Pg = x(vv.i1.Pg:vv.iN.Pg);  %% active generation in p.u.
0096 Qg = x(vv.i1.Qg:vv.iN.Qg);  %% reactive generation in p.u.
0097 
0098 %% put Pg & Qg back in gen
0099 gen(:, PG) = Pg * baseMVA;  %% active generation in MW
0100 gen(:, QG) = Qg * baseMVA;  %% reactive generation in MVAr
0101  
0102 %% rebuild Sbus
0103 Sbus = makeSbus(baseMVA, bus, gen); %% net injected power in p.u.
0104 
0105 %% ----- evaluate constraints -----
0106 %% reconstruct V
0107 Va = x(vv.i1.Va:vv.iN.Va);
0108 Vm = x(vv.i1.Vm:vv.iN.Vm);
0109 V = Vm .* exp(1j * Va);
0110 
0111 %% evaluate power flow equations
0112 mis = V .* conj(Ybus * V) - Sbus;
0113 
0114 %%----- evaluate constraint function values -----
0115 %% first, the equality constraints (power flow)
0116 g = [ real(mis);            %% active power mismatch for all buses
0117       imag(mis) ];          %% reactive power mismatch for all buses
0118 
0119 %% then, the inequality constraints (branch flow limits)
0120 if nl2 > 0
0121   flow_max = (branch(il, RATE_A)/baseMVA).^2;
0122   flow_max(flow_max == 0) = Inf;
0123   if mpopt(24) == 2       %% current magnitude limit, |I|
0124     If = Yf * V;
0125     It = Yt * V;
0126     h = [ If .* conj(If) - flow_max;      %% branch current limits (from bus)
0127           It .* conj(It) - flow_max ];    %% branch current limits (to bus)
0128   else
0129     %% compute branch power flows
0130     Sf = V(branch(il, F_BUS)) .* conj(Yf * V);  %% complex power injected at "from" bus (p.u.)
0131     St = V(branch(il, T_BUS)) .* conj(Yt * V);  %% complex power injected at "to" bus (p.u.)
0132     if mpopt(24) == 1   %% active power limit, P (Pan Wei)
0133       h = [ real(Sf).^2 - flow_max;         %% branch real power limits (from bus)
0134             real(St).^2 - flow_max ];       %% branch real power limits (to bus)
0135     else                %% apparent power limit, |S|
0136       h = [ Sf .* conj(Sf) - flow_max;      %% branch apparent power limits (from bus)
0137             St .* conj(St) - flow_max ];    %% branch apparent power limits (to bus)
0138     end
0139   end
0140 else
0141   h = zeros(0,1);
0142 end
0143 
0144 %%----- evaluate partials of constraints -----
0145 if nargout > 2
0146   %% index ranges
0147   iVa = vv.i1.Va:vv.iN.Va;
0148   iVm = vv.i1.Vm:vv.iN.Vm;
0149   iPg = vv.i1.Pg:vv.iN.Pg;
0150   iQg = vv.i1.Qg:vv.iN.Qg;
0151 
0152   %% compute partials of injected bus powers
0153   [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);           %% w.r.t. V
0154   neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);   %% Pbus w.r.t. Pg
0155                                                         %% Qbus w.r.t. Qg
0156   
0157   %% construct Jacobian of equality (power flow) constraints and transpose it
0158   dg = sparse(2*nb, nxyz);
0159   dg(:, [iVa iVm iPg iQg]) = [
0160     real([dSbus_dVa dSbus_dVm]) neg_Cg sparse(nb, ng);  %% P mismatch w.r.t Va, Vm, Pg, Qg
0161     imag([dSbus_dVa dSbus_dVm]) sparse(nb, ng) neg_Cg;  %% Q mismatch w.r.t Va, Vm, Pg, Qg
0162   ];
0163   dg = dg';
0164 
0165   if nl2 > 0
0166     %% compute partials of Flows w.r.t. V
0167     if mpopt(24) == 2     %% current
0168       [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dIbr_dV(branch(il,:), Yf, Yt, V);
0169     else                  %% power
0170       [dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft] = dSbr_dV(branch(il,:), Yf, Yt, V);
0171     end
0172     if mpopt(24) == 1     %% real part of flow (active power)
0173       dFf_dVa = real(dFf_dVa);
0174       dFf_dVm = real(dFf_dVm);
0175       dFt_dVa = real(dFt_dVa);
0176       dFt_dVm = real(dFt_dVm);
0177       Ff = real(Ff);
0178       Ft = real(Ft);
0179     end
0180   
0181     %% squared magnitude of flow (of complex power or current, or real power)
0182     [df_dVa, df_dVm, dt_dVa, dt_dVm] = ...
0183             dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft);
0184   
0185     %% construct Jacobian of inequality (branch flow) constraints & transpose
0186     dh = sparse(2*nl2, nxyz);
0187     dh(:, [iVa iVm]) = [
0188       df_dVa, df_dVm;                     %% "from" flow limit
0189       dt_dVa, dt_dVm;                     %% "to" flow limit
0190     ];
0191     dh = dh';
0192   else
0193     dh = sparse(nxyz, 0);
0194   end
0195 
0196   %% use non-sparse matrices
0197   if mpopt(51) == 0     %% hijacked SPARSE_QP to indicate MATLAB 6 => full matrices
0198     dg = full(dg);
0199     dh = full(dh);
0200   end
0201 end
opf_consfcn

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
