Source Code for Module pypower.t.t_hessian

  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 2nd derivative code. 
 18  """ 
 19   
 20  from numpy import pi, random, ones, zeros, exp 
 21  from scipy.sparse import csr_matrix as sparse 
 22   
 23  from pypower.case30 import case30 
 24  from pypower.ppoption import ppoption 
 25  from pypower.runpf import runpf 
 26  from pypower.ext2int import ext2int1 
 27  from pypower.makeYbus import makeYbus 
 28  from pypower.dSbus_dV import dSbus_dV 
 29  from pypower.d2Sbus_dV2 import d2Sbus_dV2 
 30  from pypower.dSbr_dV import dSbr_dV 
 31  from pypower.d2Sbr_dV2 import d2Sbr_dV2 
 32  from pypower.dIbr_dV import dIbr_dV 
 33  from pypower.d2Ibr_dV2 import d2Ibr_dV2 
 34  from pypower.dAbr_dV import dAbr_dV 
 35  from pypower.d2ASbr_dV2 import d2ASbr_dV2 
 36  from pypower.d2AIbr_dV2 import d2AIbr_dV2 
 37   
 38  from pypower.idx_bus import VM, VA 
 39  from pypower.idx_brch import T_BUS, F_BUS 
 40   
 41  from pypower.t.t_begin import t_begin 
 42  from pypower.t.t_is import t_is 
 43  from pypower.t.t_end import t_end 
 44   
 45   
 46 -def t_hessian(quiet=False): 
 47      """Numerical tests of 2nd derivative code. 
 48   
 49      @author: Ray Zimmerman (PSERC Cornell) 
 50      @author: Richard Lincoln 
 51      """ 
 52      t_begin(44, quiet) 
 53   
 54      ## run powerflow to get solved case 
 55      ppopt = ppoption(VERBOSE=0, OUT_ALL=0) 
 56      results, _ = runpf(case30(), ppopt) 
 57      baseMVA, bus, gen, branch = \ 
 58          results['baseMVA'], results['bus'], results['gen'], results['branch'] 
 59   
 60      ## switch to internal bus numbering and build admittance matrices 
 61      _, bus, gen, branch = ext2int1(bus, gen, branch) 
 62      Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch) 
 63      Vm = bus[:, VM] 
 64      Va = bus[:, VA] * (pi / 180) 
 65      V = Vm * exp(1j * Va) 
 66      f = branch[:, F_BUS]       ## list of "from" buses 
 67      t = branch[:, T_BUS]       ## list of "to" buses 
 68      nl = len(f) 
 69      nb = len(V) 
 70      Cf = sparse((ones(nl), (range(nl), f)), (nl, nb))  ## connection matrix for line & from buses 
 71      Ct = sparse((ones(nl), (range(nl), t)), (nl, nb))  ## connection matrix for line & to buses 
 72      pert = 1e-8 
 73   
 74      ##-----  check d2Sbus_dV2 code  ----- 
 75      t = ' - d2Sbus_dV2 (complex power injections)' 
 76      lam = 10 * random.rand(nb) 
 77      num_Haa = zeros((nb, nb), complex) 
 78      num_Hav = zeros((nb, nb), complex) 
 79      num_Hva = zeros((nb, nb), complex) 
 80      num_Hvv = zeros((nb, nb), complex) 
 81      dSbus_dVm, dSbus_dVa = dSbus_dV(Ybus, V) 
 82      Haa, Hav, Hva, Hvv = d2Sbus_dV2(Ybus, V, lam) 
 83      for i in range(nb): 
 84          Vap = V.copy() 
 85          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
 86          dSbus_dVm_ap, dSbus_dVa_ap = dSbus_dV(Ybus, Vap) 
 87          num_Haa[:, i] = (dSbus_dVa_ap - dSbus_dVa).T * lam / pert 
 88          num_Hva[:, i] = (dSbus_dVm_ap - dSbus_dVm).T * lam / pert 
 89   
 90          Vmp = V.copy() 
 91          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
 92          dSbus_dVm_mp, dSbus_dVa_mp = dSbus_dV(Ybus, Vmp) 
 93          num_Hav[:, i] = (dSbus_dVa_mp - dSbus_dVa).T * lam / pert 
 94          num_Hvv[:, i] = (dSbus_dVm_mp - dSbus_dVm).T * lam / pert 
 95   
 96      t_is(Haa.todense(), num_Haa, 4, ['Haa', t]) 
 97      t_is(Hav.todense(), num_Hav, 4, ['Hav', t]) 
 98      t_is(Hva.todense(), num_Hva, 4, ['Hva', t]) 
 99      t_is(Hvv.todense(), num_Hvv, 4, ['Hvv', t]) 
100   
101      ##-----  check d2Sbr_dV2 code  ----- 
102      t = ' - d2Sbr_dV2 (complex power flows)' 
103      lam = 10 * random.rand(nl) 
104      # lam = [1 zeros(nl-1, 1)] 
105      num_Gfaa = zeros((nb, nb), complex) 
106      num_Gfav = zeros((nb, nb), complex) 
107      num_Gfva = zeros((nb, nb), complex) 
108      num_Gfvv = zeros((nb, nb), complex) 
109      num_Gtaa = zeros((nb, nb), complex) 
110      num_Gtav = zeros((nb, nb), complex) 
111      num_Gtva = zeros((nb, nb), complex) 
112      num_Gtvv = zeros((nb, nb), complex) 
113      dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, _, _ = dSbr_dV(branch, Yf, Yt, V) 
114      Gfaa, Gfav, Gfva, Gfvv = d2Sbr_dV2(Cf, Yf, V, lam) 
115      Gtaa, Gtav, Gtva, Gtvv = d2Sbr_dV2(Ct, Yt, V, lam) 
116      for i in range(nb): 
117          Vap = V.copy() 
118          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
119          dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap = \ 
120              dSbr_dV(branch, Yf, Yt, Vap) 
121          num_Gfaa[:, i] = (dSf_dVa_ap - dSf_dVa).T * lam / pert 
122          num_Gfva[:, i] = (dSf_dVm_ap - dSf_dVm).T * lam / pert 
123          num_Gtaa[:, i] = (dSt_dVa_ap - dSt_dVa).T * lam / pert 
124          num_Gtva[:, i] = (dSt_dVm_ap - dSt_dVm).T * lam / pert 
125   
126          Vmp = V.copy() 
127          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
128          dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp = \ 
129              dSbr_dV(branch, Yf, Yt, Vmp) 
130          num_Gfav[:, i] = (dSf_dVa_mp - dSf_dVa).T * lam / pert 
131          num_Gfvv[:, i] = (dSf_dVm_mp - dSf_dVm).T * lam / pert 
132          num_Gtav[:, i] = (dSt_dVa_mp - dSt_dVa).T * lam / pert 
133          num_Gtvv[:, i] = (dSt_dVm_mp - dSt_dVm).T * lam / pert 
134   
135      t_is(Gfaa.todense(), num_Gfaa, 4, ['Gfaa', t]) 
136      t_is(Gfav.todense(), num_Gfav, 4, ['Gfav', t]) 
137      t_is(Gfva.todense(), num_Gfva, 4, ['Gfva', t]) 
138      t_is(Gfvv.todense(), num_Gfvv, 4, ['Gfvv', t]) 
139   
140      t_is(Gtaa.todense(), num_Gtaa, 4, ['Gtaa', t]) 
141      t_is(Gtav.todense(), num_Gtav, 4, ['Gtav', t]) 
142      t_is(Gtva.todense(), num_Gtva, 4, ['Gtva', t]) 
143      t_is(Gtvv.todense(), num_Gtvv, 4, ['Gtvv', t]) 
144   
145      ##-----  check d2Ibr_dV2 code  ----- 
146      t = ' - d2Ibr_dV2 (complex currents)' 
147      lam = 10 * random.rand(nl) 
148      # lam = [1, zeros(nl-1)] 
149      num_Gfaa = zeros((nb, nb), complex) 
150      num_Gfav = zeros((nb, nb), complex) 
151      num_Gfva = zeros((nb, nb), complex) 
152      num_Gfvv = zeros((nb, nb), complex) 
153      num_Gtaa = zeros((nb, nb), complex) 
154      num_Gtav = zeros((nb, nb), complex) 
155      num_Gtva = zeros((nb, nb), complex) 
156      num_Gtvv = zeros((nb, nb), complex) 
157      dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, _, _ = dIbr_dV(branch, Yf, Yt, V) 
158      Gfaa, Gfav, Gfva, Gfvv = d2Ibr_dV2(Yf, V, lam) 
159   
160      Gtaa, Gtav, Gtva, Gtvv = d2Ibr_dV2(Yt, V, lam) 
161      for i in range(nb): 
162          Vap = V.copy() 
163          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
164          dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap = \ 
165              dIbr_dV(branch, Yf, Yt, Vap) 
166          num_Gfaa[:, i] = (dIf_dVa_ap - dIf_dVa).T * lam / pert 
167          num_Gfva[:, i] = (dIf_dVm_ap - dIf_dVm).T * lam / pert 
168          num_Gtaa[:, i] = (dIt_dVa_ap - dIt_dVa).T * lam / pert 
169          num_Gtva[:, i] = (dIt_dVm_ap - dIt_dVm).T * lam / pert 
170   
171          Vmp = V.copy() 
172          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
173          dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp = \ 
174              dIbr_dV(branch, Yf, Yt, Vmp) 
175          num_Gfav[:, i] = (dIf_dVa_mp - dIf_dVa).T * lam / pert 
176          num_Gfvv[:, i] = (dIf_dVm_mp - dIf_dVm).T * lam / pert 
177          num_Gtav[:, i] = (dIt_dVa_mp - dIt_dVa).T * lam / pert 
178          num_Gtvv[:, i] = (dIt_dVm_mp - dIt_dVm).T * lam / pert 
179   
180      t_is(Gfaa.todense(), num_Gfaa, 4, ['Gfaa', t]) 
181      t_is(Gfav.todense(), num_Gfav, 4, ['Gfav', t]) 
182      t_is(Gfva.todense(), num_Gfva, 4, ['Gfva', t]) 
183      t_is(Gfvv.todense(), num_Gfvv, 4, ['Gfvv', t]) 
184   
185      t_is(Gtaa.todense(), num_Gtaa, 4, ['Gtaa', t]) 
186      t_is(Gtav.todense(), num_Gtav, 4, ['Gtav', t]) 
187      t_is(Gtva.todense(), num_Gtva, 4, ['Gtva', t]) 
188      t_is(Gtvv.todense(), num_Gtvv, 4, ['Gtvv', t]) 
189   
190      ##-----  check d2ASbr_dV2 code  ----- 
191      t = ' - d2ASbr_dV2 (squared apparent power flows)' 
192      lam = 10 * random.rand(nl) 
193      # lam = [1 zeros(nl-1, 1)] 
194      num_Gfaa = zeros((nb, nb), complex) 
195      num_Gfav = zeros((nb, nb), complex) 
196      num_Gfva = zeros((nb, nb), complex) 
197      num_Gfvv = zeros((nb, nb), complex) 
198      num_Gtaa = zeros((nb, nb), complex) 
199      num_Gtav = zeros((nb, nb), complex) 
200      num_Gtva = zeros((nb, nb), complex) 
201      num_Gtvv = zeros((nb, nb), complex) 
202      dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = dSbr_dV(branch, Yf, Yt, V) 
203      dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm = \ 
204                              dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St) 
205      Gfaa, Gfav, Gfva, Gfvv = d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, lam) 
206      Gtaa, Gtav, Gtva, Gtvv = d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, lam) 
207      for i in range(nb): 
208          Vap = V.copy() 
209          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
210          dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap = \ 
211              dSbr_dV(branch, Yf, Yt, Vap) 
212          dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap = \ 
213              dAbr_dV(dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap) 
214          num_Gfaa[:, i] = (dAf_dVa_ap - dAf_dVa).T * lam / pert 
215          num_Gfva[:, i] = (dAf_dVm_ap - dAf_dVm).T * lam / pert 
216          num_Gtaa[:, i] = (dAt_dVa_ap - dAt_dVa).T * lam / pert 
217          num_Gtva[:, i] = (dAt_dVm_ap - dAt_dVm).T * lam / pert 
218   
219          Vmp = V.copy() 
220          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
221          dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp = \ 
222              dSbr_dV(branch, Yf, Yt, Vmp) 
223          dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp = \ 
224              dAbr_dV(dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp) 
225          num_Gfav[:, i] = (dAf_dVa_mp - dAf_dVa).T * lam / pert 
226          num_Gfvv[:, i] = (dAf_dVm_mp - dAf_dVm).T * lam / pert 
227          num_Gtav[:, i] = (dAt_dVa_mp - dAt_dVa).T * lam / pert 
228          num_Gtvv[:, i] = (dAt_dVm_mp - dAt_dVm).T * lam / pert 
229   
230      t_is(Gfaa.todense(), num_Gfaa, 2, ['Gfaa', t]) 
231      t_is(Gfav.todense(), num_Gfav, 2, ['Gfav', t]) 
232      t_is(Gfva.todense(), num_Gfva, 2, ['Gfva', t]) 
233      t_is(Gfvv.todense(), num_Gfvv, 2, ['Gfvv', t]) 
234   
235      t_is(Gtaa.todense(), num_Gtaa, 2, ['Gtaa', t]) 
236      t_is(Gtav.todense(), num_Gtav, 2, ['Gtav', t]) 
237      t_is(Gtva.todense(), num_Gtva, 2, ['Gtva', t]) 
238      t_is(Gtvv.todense(), num_Gtvv, 2, ['Gtvv', t]) 
239   
240      ##-----  check d2ASbr_dV2 code  ----- 
241      t = ' - d2ASbr_dV2 (squared real power flows)' 
242      lam = 10 * random.rand(nl) 
243      # lam = [1 zeros(nl-1, 1)] 
244      num_Gfaa = zeros((nb, nb), complex) 
245      num_Gfav = zeros((nb, nb), complex) 
246      num_Gfva = zeros((nb, nb), complex) 
247      num_Gfvv = zeros((nb, nb), complex) 
248      num_Gtaa = zeros((nb, nb), complex) 
249      num_Gtav = zeros((nb, nb), complex) 
250      num_Gtva = zeros((nb, nb), complex) 
251      num_Gtvv = zeros((nb, nb), complex) 
252      dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = dSbr_dV(branch, Yf, Yt, V) 
253      dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm = \ 
254             dAbr_dV(dSf_dVa.real, dSf_dVm.real, dSt_dVa.real, dSt_dVm.real, Sf.real, St.real) 
255      Gfaa, Gfav, Gfva, Gfvv = d2ASbr_dV2(dSf_dVa.real, dSf_dVm.real, Sf.real, Cf, Yf, V, lam) 
256      Gtaa, Gtav, Gtva, Gtvv = d2ASbr_dV2(dSt_dVa.real, dSt_dVm.real, St.real, Ct, Yt, V, lam) 
257      for i in range(nb): 
258          Vap = V.copy() 
259          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
260          dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap = \ 
261              dSbr_dV(branch, Yf, Yt, Vap) 
262          dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap = \ 
263              dAbr_dV(dSf_dVa_ap.real, dSf_dVm_ap.real, dSt_dVa_ap.real, dSt_dVm_ap.real, Sf_ap.real, St_ap.real) 
264          num_Gfaa[:, i] = (dAf_dVa_ap - dAf_dVa).T * lam / pert 
265          num_Gfva[:, i] = (dAf_dVm_ap - dAf_dVm).T * lam / pert 
266          num_Gtaa[:, i] = (dAt_dVa_ap - dAt_dVa).T * lam / pert 
267          num_Gtva[:, i] = (dAt_dVm_ap - dAt_dVm).T * lam / pert 
268   
269          Vmp = V.copy() 
270          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
271          dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp = \ 
272              dSbr_dV(branch, Yf, Yt, Vmp) 
273          dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp = \ 
274              dAbr_dV(dSf_dVa_mp.real, dSf_dVm_mp.real, dSt_dVa_mp.real, dSt_dVm_mp.real, Sf_mp.real, St_mp.real) 
275          num_Gfav[:, i] = (dAf_dVa_mp - dAf_dVa).T * lam / pert 
276          num_Gfvv[:, i] = (dAf_dVm_mp - dAf_dVm).T * lam / pert 
277          num_Gtav[:, i] = (dAt_dVa_mp - dAt_dVa).T * lam / pert 
278          num_Gtvv[:, i] = (dAt_dVm_mp - dAt_dVm).T * lam / pert 
279   
280      t_is(Gfaa.todense(), num_Gfaa, 2, ['Gfaa', t]) 
281      t_is(Gfav.todense(), num_Gfav, 2, ['Gfav', t]) 
282      t_is(Gfva.todense(), num_Gfva, 2, ['Gfva', t]) 
283      t_is(Gfvv.todense(), num_Gfvv, 2, ['Gfvv', t]) 
284   
285      t_is(Gtaa.todense(), num_Gtaa, 2, ['Gtaa', t]) 
286      t_is(Gtav.todense(), num_Gtav, 2, ['Gtav', t]) 
287      t_is(Gtva.todense(), num_Gtva, 2, ['Gtva', t]) 
288      t_is(Gtvv.todense(), num_Gtvv, 2, ['Gtvv', t]) 
289   
290      ##-----  check d2AIbr_dV2 code  ----- 
291      t = ' - d2AIbr_dV2 (squared current magnitudes)' 
292      lam = 10 * random.rand(nl) 
293      # lam = [1 zeros(nl-1, 1)] 
294      num_Gfaa = zeros((nb, nb), complex) 
295      num_Gfav = zeros((nb, nb), complex) 
296      num_Gfva = zeros((nb, nb), complex) 
297      num_Gfvv = zeros((nb, nb), complex) 
298      num_Gtaa = zeros((nb, nb), complex) 
299      num_Gtav = zeros((nb, nb), complex) 
300      num_Gtva = zeros((nb, nb), complex) 
301      num_Gtvv = zeros((nb, nb), complex) 
302      dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = dIbr_dV(branch, Yf, Yt, V) 
303      dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm = \ 
304                              dAbr_dV(dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It) 
305      Gfaa, Gfav, Gfva, Gfvv = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, lam) 
306      Gtaa, Gtav, Gtva, Gtvv = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, lam) 
307      for i in range(nb): 
308          Vap = V.copy() 
309          Vap[i] = Vm[i] * exp(1j * (Va[i] + pert)) 
310          dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap = \ 
311              dIbr_dV(branch, Yf, Yt, Vap) 
312          dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap = \ 
313              dAbr_dV(dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap) 
314          num_Gfaa[:, i] = (dAf_dVa_ap - dAf_dVa).T * lam / pert 
315          num_Gfva[:, i] = (dAf_dVm_ap - dAf_dVm).T * lam / pert 
316          num_Gtaa[:, i] = (dAt_dVa_ap - dAt_dVa).T * lam / pert 
317          num_Gtva[:, i] = (dAt_dVm_ap - dAt_dVm).T * lam / pert 
318   
319          Vmp = V.copy() 
320          Vmp[i] = (Vm[i] + pert) * exp(1j * Va[i]) 
321          dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp = \ 
322              dIbr_dV(branch, Yf, Yt, Vmp) 
323          dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp = \ 
324              dAbr_dV(dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp) 
325          num_Gfav[:, i] = (dAf_dVa_mp - dAf_dVa).T * lam / pert 
326          num_Gfvv[:, i] = (dAf_dVm_mp - dAf_dVm).T * lam / pert 
327          num_Gtav[:, i] = (dAt_dVa_mp - dAt_dVa).T * lam / pert 
328          num_Gtvv[:, i] = (dAt_dVm_mp - dAt_dVm).T * lam / pert 
329   
330      t_is(Gfaa.todense(), num_Gfaa, 3, ['Gfaa', t]) 
331      t_is(Gfav.todense(), num_Gfav, 3, ['Gfav', t]) 
332      t_is(Gfva.todense(), num_Gfva, 3, ['Gfva', t]) 
333      t_is(Gfvv.todense(), num_Gfvv, 2, ['Gfvv', t]) 
334   
335      t_is(Gtaa.todense(), num_Gtaa, 3, ['Gtaa', t]) 
336      t_is(Gtav.todense(), num_Gtav, 3, ['Gtav', t]) 
337      t_is(Gtva.todense(), num_Gtva, 3, ['Gtva', t]) 
338      t_is(Gtvv.todense(), num_Gtvv, 2, ['Gtvv', t]) 
339   
340      t_end() 
341   
342   
343  if __name__ == '__main__': 
344      t_hessian(quiet=False) 
345