Source Code for Module pypower.opf_costfcn

  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 objective function, gradient and Hessian for OPF. 
 18  """ 
 19   
 20  from numpy import array, ones, zeros, arange, r_, dot, flatnonzero as find 
 21  from scipy.sparse import issparse, csr_matrix as sparse 
 22   
 23  from idx_cost import MODEL, POLYNOMIAL 
 24   
 25  from totcost import totcost 
 26  from polycost import polycost 
 27   
 28   
 29 -def opf_costfcn(x, om, return_hessian=False): 
 30      """Evaluates objective function, gradient and Hessian for OPF. 
 31   
 32      Objective function evaluation routine for AC optimal power flow, 
 33      suitable for use with L{pips}. Computes objective function value, 
 34      gradient and Hessian. 
 35   
 36      @param x: optimization vector 
 37      @param om: OPF model object 
 38   
 39      @return: C{F} - value of objective function. C{df} - (optional) gradient 
 40      of objective function (column vector). C{d2f} - (optional) Hessian of 
 41      objective function (sparse matrix). 
 42   
 43      @see: L{opf_consfcn}, L{opf_hessfcn} 
 44   
 45      @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 
 46      Autonoma de Manizales) 
 47      @author: Ray Zimmerman (PSERC Cornell) 
 48      @author: Richard Lincoln 
 49      """ 
 50      ##----- initialize ----- 
 51      ## unpack data 
 52      ppc = om.get_ppc() 
 53      baseMVA, gen, gencost = ppc["baseMVA"], ppc["gen"], ppc["gencost"] 
 54      cp = om.get_cost_params() 
 55      N, Cw, H, dd, rh, kk, mm = \ 
 56          cp["N"], cp["Cw"], cp["H"], cp["dd"], cp["rh"], cp["kk"], cp["mm"] 
 57      vv, _, _, _ = om.get_idx() 
 58   
 59      ## problem dimensions 
 60      ng = gen.shape[0]          ## number of dispatchable injections 
 61      ny = om.getN('var', 'y')   ## number of piece-wise linear costs 
 62      nxyz = len(x)              ## total number of control vars of all types 
 63   
 64      ## grab Pg & Qg 
 65      Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u. 
 66      Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u. 
 67   
 68      ##----- evaluate objective function ----- 
 69      ## polynomial cost of P and Q 
 70      # use totcost only on polynomial cost in the minimization problem 
 71      # formulation, pwl cost is the sum of the y variables. 
 72      ipol = find(gencost[:, MODEL] == POLYNOMIAL)   ## poly MW and MVAr costs 
 73      xx = r_[ Pg, Qg ] * baseMVA 
 74      if any(ipol): 
 75          f = sum( totcost(gencost[ipol, :], xx[ipol]) )  ## cost of poly P or Q 
 76      else: 
 77          f = 0 
 78   
 79      ## piecewise linear cost of P and Q 
 80      if ny > 0: 
 81          ccost = sparse((ones(ny), 
 82                          (zeros(ny), arange(vv["i1"]["y"], vv["iN"]["y"]))), 
 83                         (1, nxyz)).toarray().flatten() 
 84          f = f + dot(ccost, x) 
 85      else: 
 86          ccost = zeros(nxyz) 
 87   
 88      ## generalized cost term 
 89      if issparse(N) and N.nnz > 0: 
 90          nw = N.shape[0] 
 91          r = N * x - rh                   ## Nx - rhat 
 92          iLT = find(r < -kk)              ## below dead zone 
 93          iEQ = find((r == 0) & (kk == 0)) ## dead zone doesn't exist 
 94          iGT = find(r > kk)               ## above dead zone 
 95          iND = r_[iLT, iEQ, iGT]          ## rows that are Not in the Dead region 
 96          iL = find(dd == 1)           ## rows using linear function 
 97          iQ = find(dd == 2)           ## rows using quadratic function 
 98          LL = sparse((ones(len(iL)), (iL, iL)), (nw, nw)) 
 99          QQ = sparse((ones(len(iQ)), (iQ, iQ)), (nw, nw)) 
100          kbar = sparse((r_[ones(len(iLT)), zeros(len(iEQ)), -ones(len(iGT))], 
101                         (iND, iND)), (nw, nw)) * kk 
102          rr = r + kbar                  ## apply non-dead zone shift 
103          M = sparse((mm[iND], (iND, iND)), (nw, nw))  ## dead zone or scale 
104          diagrr = sparse((rr, (arange(nw), arange(nw))), (nw, nw)) 
105   
106          ## linear rows multiplied by rr(i), quadratic rows by rr(i)^2 
107          w = M * (LL + QQ * diagrr) * rr 
108   
109          f = f + dot(w * H, w) / 2 + dot(Cw, w) 
110   
111      ##----- evaluate cost gradient ----- 
112      ## index ranges 
113      iPg = range(vv["i1"]["Pg"], vv["iN"]["Pg"]) 
114      iQg = range(vv["i1"]["Qg"], vv["iN"]["Qg"]) 
115   
116      ## polynomial cost of P and Q 
117      df_dPgQg = zeros(2 * ng)        ## w.r.t p.u. Pg and Qg 
118      df_dPgQg[ipol] = baseMVA * polycost(gencost[ipol, :], xx[ipol], 1) 
119      df = zeros(nxyz) 
120      df[iPg] = df_dPgQg[:ng] 
121      df[iQg] = df_dPgQg[ng:ng + ng] 
122   
123      ## piecewise linear cost of P and Q 
124      df = df + ccost  # The linear cost row is additive wrt any nonlinear cost. 
125   
126      ## generalized cost term 
127      if issparse(N) and N.nnz > 0: 
128          HwC = H * w + Cw 
129          AA = N.T * M * (LL + 2 * QQ * diagrr) 
130          df = df + AA * HwC 
131   
132          ## numerical check 
133          if 0:    ## 1 to check, 0 to skip check 
134              ddff = zeros(df.shape) 
135              step = 1e-7 
136              tol  = 1e-3 
137              for k in range(len(x)): 
138                  xx = x 
139                  xx[k] = xx[k] + step 
140                  ddff[k] = (opf_costfcn(xx, om) - f) / step 
141              if max(abs(ddff - df)) > tol: 
142                  idx = find(abs(ddff - df) == max(abs(ddff - df))) 
143                  print 'Mismatch in gradient' 
144                  print 'idx             df(num)         df              diff' 
145                  print '%4d%16g%16g%16g' % \ 
146                      (range(len(df)), ddff.T, df.T, abs(ddff - df).T) 
147                  print 'MAX' 
148                  print '%4d%16g%16g%16g' % \ 
149                      (idx.T, ddff[idx].T, df[idx].T, abs(ddff[idx] - df[idx]).T) 
150   
151      if not return_hessian: 
152          return f, df 
153   
154      ## ---- evaluate cost Hessian ----- 
155      pcost = gencost[range(ng), :] 
156      if gencost.shape[0] > ng: 
157          qcost = gencost[ng + 1:2 * ng, :] 
158      else: 
159          qcost = array([]) 
160   
161      ## polynomial generator costs 
162      d2f_dPg2 = zeros(ng)               ## w.r.t. p.u. Pg 
163      d2f_dQg2 = zeros(ng)               ## w.r.t. p.u. Qg 
164      ipolp = find(pcost[:, MODEL] == POLYNOMIAL) 
165      d2f_dPg2[ipolp] = \ 
166              baseMVA**2 * polycost(pcost[ipolp, :], Pg[ipolp]*baseMVA, 2) 
167      if any(qcost):          ## Qg is not free 
168          ipolq = find(qcost[:, MODEL] == POLYNOMIAL) 
169          d2f_dQg2[ipolq] = \ 
170                  baseMVA**2 * polycost(qcost[ipolq, :], Qg[ipolq] * baseMVA, 2) 
171      i = r_[iPg, iQg].T 
172      d2f = sparse((r_[d2f_dPg2, d2f_dQg2], (i, i)), (nxyz, nxyz)) 
173   
174      ## generalized cost 
175      if N is not None and issparse(N): 
176          d2f = d2f + AA * H * AA.T + 2 * N.T * M * QQ * \ 
177                  sparse((HwC, (range(nw), range(nw))), (nw, nw)) * N 
178   
179      return f, df, d2f 
180