Source Code for Module pypower.opf_hessfcn

  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 Hessian of Lagrangian for AC OPF. 
 18  """ 
 19   
 20  from numpy import array, zeros, ones, exp, arange, r_, flatnonzero as find 
 21  from scipy.sparse import vstack, hstack, issparse, csr_matrix as sparse 
 22   
 23  from idx_gen import PG, QG 
 24  from idx_brch import F_BUS, T_BUS 
 25  from idx_cost import MODEL, POLYNOMIAL 
 26   
 27  from polycost import polycost 
 28  from d2Sbus_dV2 import d2Sbus_dV2 
 29  from dSbr_dV import dSbr_dV 
 30  from dIbr_dV import dIbr_dV 
 31  from d2AIbr_dV2 import d2AIbr_dV2 
 32  from d2ASbr_dV2 import d2ASbr_dV2 
 33  from opf_costfcn import opf_costfcn 
 34  from opf_consfcn import opf_consfcn 
 35   
 36   
 37 -def opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il=None, cost_mult=1.0): 
 38      """Evaluates Hessian of Lagrangian for AC OPF. 
 39   
 40      Hessian evaluation function for AC optimal power flow, suitable 
 41      for use with L{pips}. 
 42   
 43      Examples:: 
 44          Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt) 
 45          Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il) 
 46          Lxx = opf_hessfcn(x, lmbda, om, Ybus, Yf, Yt, ppopt, il, cost_mult) 
 47   
 48      @param x: optimization vector 
 49      @param lmbda: C{eqnonlin} - Lagrange multipliers on power balance 
 50      equations. C{ineqnonlin} - Kuhn-Tucker multipliers on constrained 
 51      branch flows. 
 52      @param om: OPF model object 
 53      @param Ybus: bus admittance matrix 
 54      @param Yf: admittance matrix for "from" end of constrained branches 
 55      @param Yt: admittance matrix for "to" end of constrained branches 
 56      @param ppopt: PYPOWER options vector 
 57      @param il: (optional) vector of branch indices corresponding to 
 58      branches with flow limits (all others are assumed to be unconstrained). 
 59      The default is C{range(nl)} (all branches). C{Yf} and C{Yt} contain 
 60      only the rows corresponding to C{il}. 
 61      @param cost_mult: (optional) Scale factor to be applied to the cost 
 62      (default = 1). 
 63   
 64      @return: Hessian of the Lagrangian. 
 65   
 66      @see: L{opf_costfcn}, L{opf_consfcn} 
 67   
 68      @author: Ray Zimmerman (PSERC Cornell) 
 69      @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad 
 70      Autonoma de Manizales) 
 71      @author: Richard Lincoln 
 72      """ 
 73      ##----- initialize ----- 
 74      ## unpack data 
 75      ppc = om.get_ppc() 
 76      baseMVA, bus, gen, branch, gencost = \ 
 77          ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"] 
 78      cp = om.get_cost_params() 
 79      N, Cw, H, dd, rh, kk, mm = \ 
 80          cp["N"], cp["Cw"], cp["H"], cp["dd"], cp["rh"], cp["kk"], cp["mm"] 
 81      vv, _, _, _ = om.get_idx() 
 82   
 83      ## unpack needed parameters 
 84      nb = bus.shape[0]          ## number of buses 
 85      nl = branch.shape[0]       ## number of branches 
 86      ng = gen.shape[0]          ## number of dispatchable injections 
 87      nxyz = len(x)              ## total number of control vars of all types 
 88   
 89      ## set default constrained lines 
 90      if il is None: 
 91          il = arange(nl)            ## all lines have limits by default 
 92      nl2 = len(il)           ## number of constrained lines 
 93   
 94      ## grab Pg & Qg 
 95      Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u. 
 96      Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u. 
 97   
 98      ## put Pg & Qg back in gen 
 99      gen[:, PG] = Pg * baseMVA  ## active generation in MW 
100      gen[:, QG] = Qg * baseMVA  ## reactive generation in MVAr 
101   
102      ## reconstruct V 
103      Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]] 
104      Vm = x[vv["i1"]["Vm"]:vv["iN"]["Vm"]] 
105      V = Vm * exp(1j * Va) 
106      nxtra = nxyz - 2 * nb 
107      pcost = gencost[arange(ng), :] 
108      if gencost.shape[0] > ng: 
109          qcost = gencost[arange(ng, 2 * ng), :] 
110      else: 
111          qcost = array([]) 
112   
113      ## ----- evaluate d2f ----- 
114      d2f_dPg2 = zeros(ng)#sparse((ng, 1))               ## w.r.t. p.u. Pg 
115      d2f_dQg2 = zeros(ng)#sparse((ng, 1))               ## w.r.t. p.u. Qg 
116      ipolp = find(pcost[:, MODEL] == POLYNOMIAL) 
117      d2f_dPg2[ipolp] = \ 
118              baseMVA**2 * polycost(pcost[ipolp, :], Pg[ipolp] * baseMVA, 2) 
119      if any(qcost):          ## Qg is not free 
120          ipolq = find(qcost[:, MODEL] == POLYNOMIAL) 
121          d2f_dQg2[ipolq] = \ 
122                  baseMVA**2 * polycost(qcost[ipolq, :], Qg[ipolq] * baseMVA, 2) 
123      i = r_[arange(vv["i1"]["Pg"], vv["iN"]["Pg"]), 
124             arange(vv["i1"]["Qg"], vv["iN"]["Qg"])] 
125  #    d2f = sparse((vstack([d2f_dPg2, d2f_dQg2]).toarray().flatten(), 
126  #                  (i, i)), shape=(nxyz, nxyz)) 
127      d2f = sparse((r_[d2f_dPg2, d2f_dQg2], (i, i)), (nxyz, nxyz)) 
128   
129      ## generalized cost 
130      if issparse(N) and N.nnz > 0: 
131          nw = N.shape[0] 
132          r = N * x - rh                    ## Nx - rhat 
133          iLT = find(r < -kk)               ## below dead zone 
134          iEQ = find((r == 0) & (kk == 0))  ## dead zone doesn't exist 
135          iGT = find(r > kk)                ## above dead zone 
136          iND = r_[iLT, iEQ, iGT]           ## rows that are Not in the Dead region 
137          iL = find(dd == 1)                ## rows using linear function 
138          iQ = find(dd == 2)                ## rows using quadratic function 
139          LL = sparse((ones(len(iL)), (iL, iL)), (nw, nw)) 
140          QQ = sparse((ones(len(iQ)), (iQ, iQ)), (nw, nw)) 
141          kbar = sparse((r_[ones(len(iLT)), zeros(len(iEQ)), -ones(len(iGT))], 
142                         (iND, iND)), (nw, nw)) * kk 
143          rr = r + kbar                  ## apply non-dead zone shift 
144          M = sparse((mm[iND], (iND, iND)), (nw, nw))  ## dead zone or scale 
145          diagrr = sparse((rr, (arange(nw), arange(nw))), (nw, nw)) 
146   
147          ## linear rows multiplied by rr(i), quadratic rows by rr(i)^2 
148          w = M * (LL + QQ * diagrr) * rr 
149          HwC = H * w + Cw 
150          AA = N.T * M * (LL + 2 * QQ * diagrr) 
151   
152          d2f = d2f + AA * H * AA.T + 2 * N.T * M * QQ * \ 
153                  sparse((HwC, (arange(nw), arange(nw))), (nw, nw)) * N 
154      d2f = d2f * cost_mult 
155   
156      ##----- evaluate Hessian of power balance constraints ----- 
157      nlam = len(lmbda["eqnonlin"]) / 2 
158      lamP = lmbda["eqnonlin"][:nlam] 
159      lamQ = lmbda["eqnonlin"][nlam:nlam + nlam] 
160      Gpaa, Gpav, Gpva, Gpvv = d2Sbus_dV2(Ybus, V, lamP) 
161      Gqaa, Gqav, Gqva, Gqvv = d2Sbus_dV2(Ybus, V, lamQ) 
162   
163      d2G = vstack([ 
164              hstack([ 
165                  vstack([hstack([Gpaa, Gpav]), 
166                          hstack([Gpva, Gpvv])]).real + 
167                  vstack([hstack([Gqaa, Gqav]), 
168                          hstack([Gqva, Gqvv])]).imag, 
169                  sparse((2 * nb, nxtra))]), 
170              hstack([ 
171                  sparse((nxtra, 2 * nb)), 
172                  sparse((nxtra, nxtra)) 
173              ]) 
174          ], "csr") 
175   
176      ##----- evaluate Hessian of flow constraints ----- 
177      nmu = len(lmbda["ineqnonlin"]) / 2 
178      muF = lmbda["ineqnonlin"][:nmu] 
179      muT = lmbda["ineqnonlin"][nmu:nmu + nmu] 
180      if ppopt['OPF_FLOW_LIM'] == 2:       ## current 
181          dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = dIbr_dV(Yf, Yt, V) 
182          Hfaa, Hfav, Hfva, Hfvv = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, muF) 
183          Htaa, Htav, Htva, Htvv = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, muT) 
184      else: 
185          f = branch[il, F_BUS].astype(int)    ## list of "from" buses 
186          t = branch[il, T_BUS].astype(int)    ## list of "to" buses 
187          ## connection matrix for line & from buses 
188          Cf = sparse((ones(nl2), (arange(nl2), f)), (nl2, nb)) 
189          ## connection matrix for line & to buses 
190          Ct = sparse((ones(nl2), (arange(nl2), t)), (nl2, nb)) 
191          dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \ 
192                  dSbr_dV(branch[il,:], Yf, Yt, V) 
193          if ppopt['OPF_FLOW_LIM'] == 1:     ## real power 
194              Hfaa, Hfav, Hfva, Hfvv = d2ASbr_dV2(dSf_dVa.real, dSf_dVm.real, 
195                                                  Sf.real, Cf, Yf, V, muF) 
196              Htaa, Htav, Htva, Htvv = d2ASbr_dV2(dSt_dVa.real, dSt_dVm.real, 
197                                                  St.real, Ct, Yt, V, muT) 
198          else:                  ## apparent power 
199              Hfaa, Hfav, Hfva, Hfvv = \ 
200                      d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, muF) 
201              Htaa, Htav, Htva, Htvv = \ 
202                      d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, muT) 
203   
204      d2H = vstack([ 
205              hstack([ 
206                  vstack([hstack([Hfaa, Hfav]), 
207                          hstack([Hfva, Hfvv])]) + 
208                  vstack([hstack([Htaa, Htav]), 
209                          hstack([Htva, Htvv])]), 
210                  sparse((2 * nb, nxtra)) 
211              ]), 
212              hstack([ 
213                  sparse((nxtra, 2 * nb)), 
214                  sparse((nxtra, nxtra)) 
215              ]) 
216          ], "csr") 
217   
218      ##-----  do numerical check using (central) finite differences  ----- 
219      if 0: 
220          nx = len(x) 
221          step = 1e-5 
222          num_d2f = sparse((nx, nx)) 
223          num_d2G = sparse((nx, nx)) 
224          num_d2H = sparse((nx, nx)) 
225          for i in range(nx): 
226              xp = x 
227              xm = x 
228              xp[i] = x[i] + step / 2 
229              xm[i] = x[i] - step / 2 
230              # evaluate cost & gradients 
231              _, dfp = opf_costfcn(xp, om) 
232              _, dfm = opf_costfcn(xm, om) 
233              # evaluate constraints & gradients 
234              _, _, dHp, dGp = opf_consfcn(xp, om, Ybus, Yf, Yt, ppopt, il) 
235              _, _, dHm, dGm = opf_consfcn(xm, om, Ybus, Yf, Yt, ppopt, il) 
236              num_d2f[:, i] = cost_mult * (dfp - dfm) / step 
237              num_d2G[:, i] = (dGp - dGm) * lmbda["eqnonlin"]   / step 
238              num_d2H[:, i] = (dHp - dHm) * lmbda["ineqnonlin"] / step 
239          d2f_err = max(max(abs(d2f - num_d2f))) 
240          d2G_err = max(max(abs(d2G - num_d2G))) 
241          d2H_err = max(max(abs(d2H - num_d2H))) 
242          if d2f_err > 1e-6: 
243              print 'Max difference in d2f: %g' % d2f_err 
244          if d2G_err > 1e-5: 
245              print 'Max difference in d2G: %g' % d2G_err 
246          if d2H_err > 1e-6: 
247              print 'Max difference in d2H: %g' % d2H_err 
248   
249      return d2f + d2G + d2H 
250