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 .

Loading