nv_powerflow.py 3.65 KB
Newer Older
Markus Mirz's avatar
Markus Mirz committed
1
import numpy as np
2
3
from .network import BusType
from .results import Results
Markus Mirz's avatar
Markus Mirz committed
4

Jan Dinkelbach's avatar
Jan Dinkelbach committed
5

6
def solve(system):
Jan Dinkelbach's avatar
Jan Dinkelbach committed
7
    """It performs powerflow by using rectangular node voltage state variables and considering the current mismatch function.
Jan Dinkelbach's avatar
Jan Dinkelbach committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
    
    Solve the non-linear powerflow problem stated by

    r = z-h(state) = 0

    following the Newton-Raphson approach

    delta_state = H^-1 * r
    new_state = old_state + delta_state

    r: residual function (current mismatch) 
    z: expected currents
    state: rectangular voltages (i.e. [V0_re, V1_re, ..., VN_re, V0_im, V1_im, ... , VN_im])
    h: currents calculated from state
    H: Jacobian matrix
    V: same as state but with complex numbers (i.e. [V0_re+j*V0_im, V1_re+j*V1_im, ...])
    """
Jan Dinkelbach's avatar
Jan Dinkelbach committed
25

26
    nodes_num = system.get_nodes_num()
Jan Dinkelbach's avatar
Jan Dinkelbach committed
27
28
29
30
    z = np.zeros(2 * nodes_num)
    h = np.zeros(2 * nodes_num)
    H = np.zeros((2 * nodes_num, 2 * nodes_num))

31
32
33
34
35
36
37
38
39
40
41
42
    for node in system.nodes:
        if node.ideal_connected_with =='':
            i = node.index
            m = 2 * i
            i2 = i + nodes_num
            node_type = node.type
            if node_type == BusType.SLACK:
                z[m] = np.real(node.voltage_pu)
                z[m + 1] = np.imag(node.voltage_pu)
                H[m][i] = 1
                H[m + 1][i2] = 1
            elif node_type is BusType.PQ:
43
44
45
46
                H[m][:node.num] = numpy.real(Ymatrix[i])
                H[m][node.num:] = - numpy.imag(Ymatrix[i])
                H[m+1][:node.num] = numpy.imag(Ymatrix[i])
                H[m+1][node.num:] = numpy.real(Ymatrix[i])
47
48
            elif node_type is BusType.PV:
                z[m + 1] = np.real(node.power)
49
50
                H[m][:node.num] = numpy.real(Ymatrix[i])
                H[m][node.num:] = - numpy.imag(Ymatrix[i])
Jan Dinkelbach's avatar
Jan Dinkelbach committed
51
52

    epsilon = 10 ** (-10)
martin.moraga's avatar
martin.moraga committed
53
    #epsilon = 0.01
Jan Dinkelbach's avatar
Jan Dinkelbach committed
54
55
56
57
58
59
60
    diff = 5
    V = np.ones(nodes_num) + 1j * np.zeros(nodes_num)
    num_iter = 0

    state = np.concatenate((np.ones(nodes_num), np.zeros(nodes_num)), axis=0)

    while diff > epsilon:
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
        for node in system.nodes:
            if node.ideal_connected_with =='':
                i = node.index
                m = 2 * i
                i2 = i + nodes_num
                node_type = node.type
                if node_type is BusType.SLACK:
                    h[m] = np.inner(H[m], state)
                    h[m + 1] = np.inner(H[m + 1], state)
                elif node_type is BusType.PQ:
                    z[m] = (np.real(node.power_pu) * np.real(V[i]) + 
                           np.imag(node.power_pu) * np.imag(V[i])) / (np.abs(V[i]) ** 2)
                    z[m + 1] = (np.real(node.power_pu) * np.imag(V[i]) - 
                           np.imag(node.power_pu) * np.real(V[i])) / (np.abs(V[i]) ** 2)
                    h[m] = np.inner(H[m], state)
                    h[m + 1] = np.inner(H[m + 1], state)
                elif node_type is BusType.PV:
                    z[m] = (np.real(node.power_pu) * np.real(V[i]) + 
79
                            np.imag(node.power_pu) * np.imag(V[i])) / (np.abs(V[i]) ** 2)
80
81
82
83
                    h[m] = np.inner(H[m], state)
                    h[m + 1] = np.abs(V[i])
                    H[m + 1][i] = np.cos(np.angle(V[i]))
                    H[m + 1][i2] = np.sin(np.angle(V[i]))
Jan Dinkelbach's avatar
Jan Dinkelbach committed
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

        r = np.subtract(z, h)
        Hinv = np.linalg.inv(H)
        delta_state = np.inner(Hinv, r)
        state = state + delta_state
        diff = np.amax(np.absolute(delta_state))

        V = state[:nodes_num] + 1j * state[nodes_num:]
        num_iter = num_iter + 1

    # calculate all the other quantities of the grid
    powerflow_results = Results(system)
    powerflow_results.load_voltages(V)
    powerflow_results.calculate_all()

    return powerflow_results, num_iter