Source Code for Module pypower.t.t_jacobian

  1  # Copyright (C) 2004-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  """Numerical tests of partial derivative code. 
 18  """ 
 19   
 20  from numpy import ones, conj, eye, exp, pi, array 
 21   
 22  from pypower.case30 import case30 
 23  from pypower.ppoption import ppoption 
 24  from pypower.loadcase import loadcase 
 25  from pypower.ext2int import ext2int1 
 26  from pypower.runpf import runpf 
 27  from pypower.makeYbus import makeYbus 
 28  from pypower.dSbus_dV import dSbus_dV 
 29  from pypower.dSbr_dV import dSbr_dV 
 30  from pypower.dAbr_dV import dAbr_dV 
 31  from pypower.dIbr_dV import dIbr_dV 
 32   
 33  from pypower.idx_bus import VM, VA 
 34  from pypower.idx_brch import F_BUS, T_BUS 
 35   
 36  from pypower.t.t_begin import t_begin 
 37  from pypower.t.t_end import t_end 
 38  from pypower.t.t_is import t_is 
 39   
 40   
 41 -def t_jacobian(quiet=False): 
 42      """Numerical tests of partial derivative code. 
 43   
 44      @author: Ray Zimmerman (PSERC Cornell) 
 45      @author: Richard Lincoln 
 46      """ 
 47      t_begin(28, quiet) 
 48   
 49      ## run powerflow to get solved case 
 50      ppopt = ppoption(VERBOSE=0, OUT_ALL=0) 
 51      ppc = loadcase(case30()) 
 52   
 53      results, _ = runpf(ppc, ppopt) 
 54      baseMVA, bus, gen, branch = \ 
 55          results['baseMVA'], results['bus'], results['gen'], results['branch'] 
 56   
 57      ## switch to internal bus numbering and build admittance matrices 
 58      _, bus, gen, branch = ext2int1(bus, gen, branch) 
 59      Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch) 
 60      Ybus_full = Ybus.todense() 
 61      Yf_full   = Yf.todense() 
 62      Yt_full   = Yt.todense() 
 63      Vm = bus[:, VM] 
 64      Va = bus[:, VA] * (pi / 180) 
 65      V = Vm * exp(1j * Va) 
 66      f = branch[:, F_BUS].astype(int)       ## list of "from" buses 
 67      t = branch[:, T_BUS].astype(int)       ## list of "to" buses 
 68      #nl = len(f) 
 69      nb = len(V) 
 70      pert = 1e-8 
 71   
 72      Vm = array([Vm]).T  # column array 
 73      Va = array([Va]).T  # column array 
 74      Vc = array([V]).T   # column array 
 75   
 76      ##-----  check dSbus_dV code  ----- 
 77      ## full matrices 
 78      dSbus_dVm_full, dSbus_dVa_full = dSbus_dV(Ybus_full, V) 
 79   
 80      ## sparse matrices 
 81      dSbus_dVm, dSbus_dVa = dSbus_dV(Ybus, V) 
 82      dSbus_dVm_sp = dSbus_dVm.todense() 
 83      dSbus_dVa_sp = dSbus_dVa.todense() 
 84   
 85      ## compute numerically to compare 
 86      Vmp = (Vm * ones((1, nb)) + pert*eye(nb)) * (exp(1j * Va) * ones((1, nb))) 
 87      Vap = (Vm * ones((1, nb))) * (exp(1j * (Va*ones((1, nb)) + pert*eye(nb)))) 
 88      num_dSbus_dVm = (Vmp * conj(Ybus * Vmp) - Vc * ones((1, nb)) * conj(Ybus * Vc * ones((1, nb)))) / pert 
 89      num_dSbus_dVa = (Vap * conj(Ybus * Vap) - Vc * ones((1, nb)) * conj(Ybus * Vc * ones((1, nb)))) / pert 
 90   
 91      t_is(dSbus_dVm_sp, num_dSbus_dVm, 5, 'dSbus_dVm (sparse)') 
 92      t_is(dSbus_dVa_sp, num_dSbus_dVa, 5, 'dSbus_dVa (sparse)') 
 93      t_is(dSbus_dVm_full, num_dSbus_dVm, 5, 'dSbus_dVm (full)') 
 94      t_is(dSbus_dVa_full, num_dSbus_dVa, 5, 'dSbus_dVa (full)') 
 95   
 96      ##-----  check dSbr_dV code  ----- 
 97      ## full matrices 
 98      dSf_dVa_full, dSf_dVm_full, dSt_dVa_full, dSt_dVm_full, _, _ = \ 
 99              dSbr_dV(branch, Yf_full, Yt_full, V) 
100   
101      ## sparse matrices 
102      dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = dSbr_dV(branch, Yf, Yt, V) 
103      dSf_dVa_sp = dSf_dVa.todense() 
104      dSf_dVm_sp = dSf_dVm.todense() 
105      dSt_dVa_sp = dSt_dVa.todense() 
106      dSt_dVm_sp = dSt_dVm.todense() 
107   
108      ## compute numerically to compare 
109      Vmpf = Vmp[f, :] 
110      Vapf = Vap[f, :] 
111      Vmpt = Vmp[t, :] 
112      Vapt = Vap[t, :] 
113      Sf2 = (Vc[f] * ones((1, nb))) * conj(Yf * Vc * ones((1, nb))) 
114      St2 = (Vc[t] * ones((1, nb))) * conj(Yt * Vc * ones((1, nb))) 
115      Smpf = Vmpf * conj(Yf * Vmp) 
116      Sapf = Vapf * conj(Yf * Vap) 
117      Smpt = Vmpt * conj(Yt * Vmp) 
118      Sapt = Vapt * conj(Yt * Vap) 
119   
120      num_dSf_dVm = (Smpf - Sf2) / pert 
121      num_dSf_dVa = (Sapf - Sf2) / pert 
122      num_dSt_dVm = (Smpt - St2) / pert 
123      num_dSt_dVa = (Sapt - St2) / pert 
124   
125      t_is(dSf_dVm_sp, num_dSf_dVm, 5, 'dSf_dVm (sparse)') 
126      t_is(dSf_dVa_sp, num_dSf_dVa, 5, 'dSf_dVa (sparse)') 
127      t_is(dSt_dVm_sp, num_dSt_dVm, 5, 'dSt_dVm (sparse)') 
128      t_is(dSt_dVa_sp, num_dSt_dVa, 5, 'dSt_dVa (sparse)') 
129      t_is(dSf_dVm_full, num_dSf_dVm, 5, 'dSf_dVm (full)') 
130      t_is(dSf_dVa_full, num_dSf_dVa, 5, 'dSf_dVa (full)') 
131      t_is(dSt_dVm_full, num_dSt_dVm, 5, 'dSt_dVm (full)') 
132      t_is(dSt_dVa_full, num_dSt_dVa, 5, 'dSt_dVa (full)') 
133   
134      ##-----  check dAbr_dV code  ----- 
135      ## full matrices 
136      dAf_dVa_full, dAf_dVm_full, dAt_dVa_full, dAt_dVm_full = \ 
137          dAbr_dV(dSf_dVa_full, dSf_dVm_full, dSt_dVa_full, dSt_dVm_full, Sf, St) 
138      ## sparse matrices 
139      dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm = \ 
140                              dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St) 
141      dAf_dVa_sp = dAf_dVa.todense() 
142      dAf_dVm_sp = dAf_dVm.todense() 
143      dAt_dVa_sp = dAt_dVa.todense() 
144      dAt_dVm_sp = dAt_dVm.todense() 
145   
146      ## compute numerically to compare 
147      num_dAf_dVm = (abs(Smpf)**2 - abs(Sf2)**2) / pert 
148      num_dAf_dVa = (abs(Sapf)**2 - abs(Sf2)**2) / pert 
149      num_dAt_dVm = (abs(Smpt)**2 - abs(St2)**2) / pert 
150      num_dAt_dVa = (abs(Sapt)**2 - abs(St2)**2) / pert 
151   
152      t_is(dAf_dVm_sp, num_dAf_dVm, 4, 'dAf_dVm (sparse)') 
153      t_is(dAf_dVa_sp, num_dAf_dVa, 4, 'dAf_dVa (sparse)') 
154      t_is(dAt_dVm_sp, num_dAt_dVm, 4, 'dAt_dVm (sparse)') 
155      t_is(dAt_dVa_sp, num_dAt_dVa, 4, 'dAt_dVa (sparse)') 
156      t_is(dAf_dVm_full, num_dAf_dVm, 4, 'dAf_dVm (full)') 
157      t_is(dAf_dVa_full, num_dAf_dVa, 4, 'dAf_dVa (full)') 
158      t_is(dAt_dVm_full, num_dAt_dVm, 4, 'dAt_dVm (full)') 
159      t_is(dAt_dVa_full, num_dAt_dVa, 4, 'dAt_dVa (full)') 
160   
161      ##-----  check dIbr_dV code  ----- 
162      ## full matrices 
163      dIf_dVa_full, dIf_dVm_full, dIt_dVa_full, dIt_dVm_full, _, _ = \ 
164              dIbr_dV(branch, Yf_full, Yt_full, V) 
165   
166      ## sparse matrices 
167      dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, _, _ = dIbr_dV(branch, Yf, Yt, V) 
168      dIf_dVa_sp = dIf_dVa.todense() 
169      dIf_dVm_sp = dIf_dVm.todense() 
170      dIt_dVa_sp = dIt_dVa.todense() 
171      dIt_dVm_sp = dIt_dVm.todense() 
172   
173      ## compute numerically to compare 
174      num_dIf_dVm = (Yf * Vmp - Yf * Vc * ones((1, nb))) / pert 
175      num_dIf_dVa = (Yf * Vap - Yf * Vc * ones((1, nb))) / pert 
176      num_dIt_dVm = (Yt * Vmp - Yt * Vc * ones((1, nb))) / pert 
177      num_dIt_dVa = (Yt * Vap - Yt * Vc * ones((1, nb))) / pert 
178   
179      t_is(dIf_dVm_sp, num_dIf_dVm, 5, 'dIf_dVm (sparse)') 
180      t_is(dIf_dVa_sp, num_dIf_dVa, 5, 'dIf_dVa (sparse)') 
181      t_is(dIt_dVm_sp, num_dIt_dVm, 5, 'dIt_dVm (sparse)') 
182      t_is(dIt_dVa_sp, num_dIt_dVa, 5, 'dIt_dVa (sparse)') 
183      t_is(dIf_dVm_full, num_dIf_dVm, 5, 'dIf_dVm (full)') 
184      t_is(dIf_dVa_full, num_dIf_dVa, 5, 'dIf_dVa (full)') 
185      t_is(dIt_dVm_full, num_dIt_dVm, 5, 'dIt_dVm (full)') 
186      t_is(dIt_dVa_full, num_dIt_dVa, 5, 'dIt_dVa (full)') 
187   
188      t_end() 
189   
190   
191  if __name__ == "__main__": 
192      t_jacobian(quiet=False) 
193