Source Code for Module pypower.opf_consfcn

  1  # Copyright (C) 1996-2011 Power System Engineering Research Center (PSERC) 
  2  # Copyright (C) 2011 Richard Lincoln 
  3  # 
  4  # PYPOWER is free software: you can redistribute it and/or modify 
  5  # it under the terms of the GNU General Public License as published 
  6  # by the Free Software Foundation, either version 3 of the License, 
  7  # or (at your option) any later version. 
  8  # 
  9  # PYPOWER is distributed in the hope that it will be useful, 
 10  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 11  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 12  # GNU General Public License for more details. 
 13  # 
 14  # You should have received a copy of the GNU General Public License 
 15  # along with PYPOWER. If not, see <http://www.gnu.org/licenses/>. 
 16   
 17  """Evaluates nonlinear constraints and their Jacobian for OPF. 
 18  """ 
 19   
 20  from numpy import zeros, ones, conj, exp, r_, Inf, arange 
 21   
 22  from scipy.sparse import lil_matrix, vstack, hstack, csr_matrix as sparse 
 23   
 24  from idx_gen import GEN_BUS, PG, QG 
 25  from idx_brch import F_BUS, T_BUS, RATE_A 
 26   
 27  from makeSbus import makeSbus 
 28  from dSbus_dV import dSbus_dV 
 29  from dIbr_dV import dIbr_dV 
 30  from dSbr_dV import dSbr_dV 
 31  from dAbr_dV import dAbr_dV 
 32   
 33   
 34 -def opf_consfcn(x, om, Ybus, Yf, Yt, ppopt, il=None, *args): 
 35      """Evaluates nonlinear constraints and their Jacobian for OPF. 
 36   
 37      Constraint evaluation function for AC optimal power flow, suitable 
 38      for use with L{pips}. Computes constraint vectors and their gradients. 
 39   
 40      @param x: optimization vector 
 41      @param om: OPF model object 
 42      @param Ybus: bus admittance matrix 
 43      @param Yf: admittance matrix for "from" end of constrained branches 
 44      @param Yt: admittance matrix for "to" end of constrained branches 
 45      @param ppopt: PYPOWER options vector 
 46      @param il: (optional) vector of branch indices corresponding to 
 47      branches with flow limits (all others are assumed to be 
 48      unconstrained). The default is C{range(nl)} (all branches). 
 49      C{Yf} and C{Yt} contain only the rows corresponding to C{il}. 
 50   
 51      @return: C{h} - vector of inequality constraint values (flow limits) 
 52      limit^2 - flow^2, where the flow can be apparent power real power or 
 53      current, depending on value of C{OPF_FLOW_LIM} in C{ppopt} (only for 
 54      constrained lines). C{g} - vector of equality constraint values (power 
 55      balances). C{dh} - (optional) inequality constraint gradients, column 
 56      j is gradient of h(j). C{dg} - (optional) equality constraint gradients. 
 57   
 58      @see: L{opf_costfcn}, L{opf_hessfcn} 
 59   
 60      @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 
 61      Autonoma de Manizales) 
 62      @author: Ray Zimmerman (PSERC Cornell) 
 63      @author: Richard Lincoln 
 64      """ 
 65      ##----- initialize ----- 
 66   
 67      ## unpack data 
 68      ppc = om.get_ppc() 
 69      baseMVA, bus, gen, branch = \ 
 70          ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"] 
 71      vv, _, _, _ = om.get_idx() 
 72   
 73      ## problem dimensions 
 74      nb = bus.shape[0]          ## number of buses 
 75      nl = branch.shape[0]       ## number of branches 
 76      ng = gen.shape[0]          ## number of dispatchable injections 
 77      nxyz = len(x)              ## total number of control vars of all types 
 78   
 79      ## set default constrained lines 
 80      if il is None: 
 81          il = arange(nl)         ## all lines have limits by default 
 82      nl2 = len(il)              ## number of constrained lines 
 83   
 84      ## grab Pg & Qg 
 85      Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u. 
 86      Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u. 
 87   
 88      ## put Pg & Qg back in gen 
 89      gen[:, PG] = Pg * baseMVA  ## active generation in MW 
 90      gen[:, QG] = Qg * baseMVA  ## reactive generation in MVAr 
 91   
 92      ## rebuild Sbus 
 93      Sbus = makeSbus(baseMVA, bus, gen) ## net injected power in p.u. 
 94   
 95      ## ----- evaluate constraints ----- 
 96      ## reconstruct V 
 97      Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]] 
 98      Vm = x[vv["i1"]["Vm"]:vv["iN"]["Vm"]] 
 99      V = Vm * exp(1j * Va) 
100   
101      ## evaluate power flow equations 
102      mis = V * conj(Ybus * V) - Sbus 
103   
104      ##----- evaluate constraint function values ----- 
105      ## first, the equality constraints (power flow) 
106      g = r_[ mis.real,            ## active power mismatch for all buses 
107              mis.imag ]           ## reactive power mismatch for all buses 
108   
109      ## then, the inequality constraints (branch flow limits) 
110      if nl2 > 0: 
111          flow_max = (branch[il, RATE_A] / baseMVA)**2 
112          flow_max[flow_max == 0] = Inf 
113          if ppopt['OPF_FLOW_LIM'] == 2:       ## current magnitude limit, |I| 
114              If = Yf * V 
115              It = Yt * V 
116              h = r_[ If * conj(If) - flow_max,     ## branch I limits (from bus) 
117                      It * conj(It) - flow_max ].real    ## branch I limits (to bus) 
118          else: 
119              ## compute branch power flows 
120              ## complex power injected at "from" bus (p.u.) 
121              Sf = V[ branch[il, F_BUS].astype(int) ] * conj(Yf * V) 
122              ## complex power injected at "to" bus (p.u.) 
123              St = V[ branch[il, T_BUS].astype(int) ] * conj(Yt * V) 
124              if ppopt['OPF_FLOW_LIM'] == 1:   ## active power limit, P (Pan Wei) 
125                  h = r_[ Sf.real**2 - flow_max,   ## branch P limits (from bus) 
126                          St.real**2 - flow_max ]  ## branch P limits (to bus) 
127              else:                ## apparent power limit, |S| 
128                  h = r_[ Sf * conj(Sf) - flow_max, ## branch S limits (from bus) 
129                          St * conj(St) - flow_max ].real  ## branch S limits (to bus) 
130      else: 
131          h = zeros((0,1)) 
132   
133      ##----- evaluate partials of constraints ----- 
134      ## index ranges 
135      iVa = arange(vv["i1"]["Va"], vv["iN"]["Va"]) 
136      iVm = arange(vv["i1"]["Vm"], vv["iN"]["Vm"]) 
137      iPg = arange(vv["i1"]["Pg"], vv["iN"]["Pg"]) 
138      iQg = arange(vv["i1"]["Qg"], vv["iN"]["Qg"]) 
139      iVaVmPgQg = r_[iVa, iVm, iPg, iQg].T 
140   
141      ## compute partials of injected bus powers 
142      dSbus_dVm, dSbus_dVa = dSbus_dV(Ybus, V)           ## w.r.t. V 
143      ## Pbus w.r.t. Pg, Qbus w.r.t. Qg 
144      neg_Cg = sparse((-ones(ng), (gen[:, GEN_BUS], range(ng))), (nb, ng)) 
145   
146      ## construct Jacobian of equality constraints (power flow) and transpose it 
147      dg = lil_matrix((2 * nb, nxyz)) 
148      blank = sparse((nb, ng)) 
149      dg[:, iVaVmPgQg] = vstack([ 
150              ## P mismatch w.r.t Va, Vm, Pg, Qg 
151              hstack([dSbus_dVa.real, dSbus_dVm.real, neg_Cg, blank]), 
152              ## Q mismatch w.r.t Va, Vm, Pg, Qg 
153              hstack([dSbus_dVa.imag, dSbus_dVm.imag, blank, neg_Cg]) 
154          ], "csr") 
155      dg = dg.T 
156   
157      if nl2 > 0: 
158          ## compute partials of Flows w.r.t. V 
159          if ppopt['OPF_FLOW_LIM'] == 2:     ## current 
160              dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \ 
161                      dIbr_dV(branch[il, :], Yf, Yt, V) 
162          else:                  ## power 
163              dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft = \ 
164                      dSbr_dV(branch[il, :], Yf, Yt, V) 
165          if ppopt['OPF_FLOW_LIM'] == 1:     ## real part of flow (active power) 
166              dFf_dVa = dFf_dVa.real 
167              dFf_dVm = dFf_dVm.real 
168              dFt_dVa = dFt_dVa.real 
169              dFt_dVm = dFt_dVm.real 
170              Ff = Ff.real 
171              Ft = Ft.real 
172   
173          ## squared magnitude of flow (of complex power or current, or real power) 
174          df_dVa, df_dVm, dt_dVa, dt_dVm = \ 
175                  dAbr_dV(dFf_dVa, dFf_dVm, dFt_dVa, dFt_dVm, Ff, Ft) 
176   
177          ## construct Jacobian of inequality constraints (branch limits) 
178          ## and transpose it. 
179          dh = lil_matrix((2 * nl2, nxyz)) 
180          dh[:, r_[iVa, iVm].T] = vstack([ 
181                  hstack([df_dVa, df_dVm]),    ## "from" flow limit 
182                  hstack([dt_dVa, dt_dVm])     ## "to" flow limit 
183              ], "csr") 
184          dh = dh.T 
185      else: 
186          dh = None 
187   
188      return h, g, dh, dg 
189