e2nIEE / pandapower
 1 ```# -*- coding: utf-8 -*- ``` 2 3 ```# Copyright 1996-2015 PSERC. All rights reserved. ``` 4 ```# Use of this source code is governed by a BSD-style ``` 5 ```# license that can be found in the LICENSE file. ``` 6 7 ```# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics ``` 8 ```# and Energy System Technology (IEE), Kassel. All rights reserved. ``` 9 10 11 1 ```"""Solves the power flow using a full Newton's method. ``` 12 ```""" ``` 13 14 1 ```from numpy import angle, exp, linalg, conj, r_, Inf, arange, zeros, max, zeros_like, column_stack ``` 15 1 ```from scipy.sparse.linalg import spsolve ``` 16 17 1 ```from pandapower.pf.iwamoto_multiplier import _iwamoto_step ``` 18 1 ```from pandapower.pypower.makeSbus import makeSbus ``` 19 1 ```from pandapower.pf.create_jacobian import create_jacobian_matrix, get_fastest_jacobian_function ``` 20 21 22 1 ```def newtonpf(Ybus, Sbus, V0, pv, pq, ppci, options): ``` 23 ``` """Solves the power flow using a full Newton's method. ``` 24 ``` Solves for bus voltages given the full system admittance matrix (for ``` 25 ``` all buses), the complex bus power injection vector (for all buses), ``` 26 ``` the initial vector of complex bus voltages, and column vectors with ``` 27 ``` the lists of bus indices for the swing bus, PV buses, and PQ buses, ``` 28 ``` respectively. The bus voltage vector contains the set point for ``` 29 ``` generator (including ref bus) buses, and the reference angle of the ``` 30 ``` swing bus, as well as an initial guess for remaining magnitudes and ``` 31 ``` angles. ``` 32 ``` @see: L{runpf} ``` 33 ``` @author: Ray Zimmerman (PSERC Cornell) ``` 34 ``` @author: Richard Lincoln ``` 35 ``` Modified by University of Kassel (Florian Schaefer) to use numba ``` 36 ``` """ ``` 37 38 ``` # options ``` 39 1 ``` tol = options['tolerance_mva'] ``` 40 1 ``` max_it = options["max_iteration"] ``` 41 1 ``` numba = options["numba"] ``` 42 1 ``` iwamoto = options["algorithm"] == "iwamoto_nr" ``` 43 1 ``` voltage_depend_loads = options["voltage_depend_loads"] ``` 44 1 ``` v_debug = options["v_debug"] ``` 45 1 ``` use_umfpack = options["use_umfpack"] ``` 46 1 ``` permc_spec = options["permc_spec"] ``` 47 48 1 ``` baseMVA = ppci['baseMVA'] ``` 49 1 ``` bus = ppci['bus'] ``` 50 1 ``` gen = ppci['gen'] ``` 51 52 ``` # initialize ``` 53 1 ``` i = 0 ``` 54 1 ``` V = V0 ``` 55 1 ``` Va = angle(V) ``` 56 1 ``` Vm = abs(V) ``` 57 1 ``` dVa, dVm = None, None ``` 58 1 ``` if iwamoto: ``` 59 0 ``` dVm, dVa = zeros_like(Vm), zeros_like(Va) ``` 60 61 1 ``` if v_debug: ``` 62 0 ``` Vm_it = Vm.copy() ``` 63 0 ``` Va_it = Va.copy() ``` 64 ``` else: ``` 65 1 ``` Vm_it = None ``` 66 1 ``` Va_it = None ``` 67 68 ``` # set up indexing for updating V ``` 69 1 ``` pvpq = r_[pv, pq] ``` 70 ``` # generate lookup pvpq -> index pvpq (used in createJ) ``` 71 1 ``` pvpq_lookup = zeros(max(Ybus.indices) + 1, dtype=int) ``` 72 1 ``` pvpq_lookup[pvpq] = arange(len(pvpq)) ``` 73 74 ``` # get jacobian function ``` 75 1 ``` createJ = get_fastest_jacobian_function(pvpq, pq, numba) ``` 76 77 1 ``` npv = len(pv) ``` 78 1 ``` npq = len(pq) ``` 79 1 ``` j1 = 0 ``` 80 1 ``` j2 = npv # j1:j2 - V angle of pv buses ``` 81 1 ``` j3 = j2 ``` 82 1 ``` j4 = j2 + npq # j3:j4 - V angle of pq buses ``` 83 1 ``` j5 = j4 ``` 84 1 ``` j6 = j4 + npq # j5:j6 - V mag of pq buses ``` 85 86 ``` # evaluate F(x0) ``` 87 1 ``` F = _evaluate_Fx(Ybus, V, Sbus, pv, pq) ``` 88 1 ``` converged = _check_for_convergence(F, tol) ``` 89 90 1 ``` Ybus = Ybus.tocsr() ``` 91 1 ``` J = None ``` 92 93 ``` # do Newton iterations ``` 94 1 ``` while (not converged and i < max_it): ``` 95 ``` # update iteration counter ``` 96 1 ``` i = i + 1 ``` 97 98 1 ``` J = create_jacobian_matrix(Ybus, V, pvpq, pq, createJ, pvpq_lookup, npv, npq, numba) ``` 99 100 1 ``` dx = -1 * spsolve(J, F, permc_spec=permc_spec, use_umfpack=use_umfpack) ``` 101 ``` # update voltage ``` 102 1 ``` if npv and not iwamoto: ``` 103 1 ``` Va[pv] = Va[pv] + dx[j1:j2] ``` 104 1 ``` if npq and not iwamoto: ``` 105 1 ``` Va[pq] = Va[pq] + dx[j3:j4] ``` 106 1 ``` Vm[pq] = Vm[pq] + dx[j5:j6] ``` 107 108 ``` # iwamoto multiplier to increase convergence ``` 109 1 ``` if iwamoto: ``` 110 0 ``` Vm, Va = _iwamoto_step(Ybus, J, F, dx, pq, npv, npq, dVa, dVm, Vm, Va, pv, j1, j2, j3, j4, j5, j6) ``` 111 112 1 ``` V = Vm * exp(1j * Va) ``` 113 1 ``` Vm = abs(V) # update Vm and Va again in case ``` 114 1 ``` Va = angle(V) # we wrapped around with a negative Vm ``` 115 116 1 ``` if v_debug: ``` 117 0 ``` Vm_it = column_stack((Vm_it, Vm)) ``` 118 0 ``` Va_it = column_stack((Va_it, Va)) ``` 119 120 1 ``` if voltage_depend_loads: ``` 121 1 ``` Sbus = makeSbus(baseMVA, bus, gen, vm=Vm) ``` 122 123 1 ``` F = _evaluate_Fx(Ybus, V, Sbus, pv, pq) ``` 124 125 1 ``` converged = _check_for_convergence(F, tol) ``` 126 127 1 ``` return V, converged, i, J, Vm_it, Va_it ``` 128 129 130 1 ```def _evaluate_Fx(Ybus, V, Sbus, pv, pq): ``` 131 ``` # evalute F(x) ``` 132 1 ``` mis = V * conj(Ybus * V) - Sbus ``` 133 1 ``` F = r_[mis[pv].real, ``` 134 ``` mis[pq].real, ``` 135 ``` mis[pq].imag] ``` 136 1 ``` return F ``` 137 138 139 1 ```def _check_for_convergence(F, tol): ``` 140 ``` # calc infinity norm ``` 141 1 ``` return linalg.norm(F, Inf) < tol ```

Read our documentation on viewing source code .