The dynamics should be separated according to what is calculated on the network links:

complex current
real power
complex power (?)

PSD/classic_model.py

-from __future__ import absolute_import, print_function, division, unicode_literals
-import numpy as np
-import networkx as nx
-import scipy.sparse as sps
-
-
-def define_oscillator_rhs(dim):
-    nn = dim // 2
-
-    def right_hand_side(y, t, coupling, infeed, droop):
-        phases = np.exp(1.j * y[:nn])
-        return np.append(y[nn:], infeed - droop * y[nn:] -
-                         np.imag(phases * coupling.dot(phases).conjugate()))
-
-    return right_hand_side
-
-
-def find_fixpoint(self, args, fp_guess=None):
-    from scipy.optimize import root
-
-    if fp_guess is None:
-        fp_guess = np.zeros(self.size_of_system * 2, dtype=np.float64)
-
-    res = root(lambda x: self._rhs(x, 0., *args), fp_guess)
-
-    if res.success:
-        return res.x
-    else:
-        print("Failed to find fix point")
-        return None
-
-
-def calculate_rhs_args(G):
-    coupling = 8. * sps.csr_matrix(nx.adj_matrix(G))
-    size = len(G)
-    droop = 0.1
-    infeed = np.ones(size, dtype=np.float64)
-    infeed[:size // 2] = -1.
-    infeed[-1] -= np.sum(infeed)
-    return coupling, infeed, droop, size
-
-
-def define_rc_gen(G, brp):
-    coupling, infeed, droop, size = calculate_rhs_args(G)
-    num_batch_inv = 1. / brp.number_of_batches
-
-    def gen(batch, run):
-        ic = np.random.random(brp.system_dimension) - 0.5
-        ic[:size] *= 2. * np.pi
-        ic[size:] *= 1.0
-        return ic, (batch * num_batch_inv * coupling, infeed, droop)

PSD/classic_model_prototype.py

+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+from numba import njit, float64, complex128, void, int32
+# from assimulo.solvers import CVODE
+# from assimulo.problem import Explicit_Problem
+import scipy.sparse as sps
+
+""" The variables at the nodes are theta and omega """
+
+c = sps.csr_matrix(np.ones((2,2)))
+print(c.dot(np.ones(2)))
+
+
+@njit(complex128[:](complex128[:], int32[:], int32[:], complex128[:]))
+def numba_sp_dot(data, indptr, indices, v):
+    res = np.zeros(len(indptr) - 1, dtype=np.complex128)
+    index = 0
+    for row, number_of_entries in enumerate(indptr[1:]):
+        while index < number_of_entries:
+            res[row] += data[index] * v[indices[index]]
+            index += 1
+    return res
+
+
+def numba_sp_powerflows(data, indptr, indices, phases, P):
+    index = 0
+    for row, number_of_entries in enumerate(indptr[1:]):
+        r_temp = 0.
+        while index < number_of_entries:
+            r_temp += data[index] * phases[indices[index]]
+            index += 1
+        P[row] += phases[row] * np.conjugate(r_temp)
+
+
+# noinspection PyPep8Naming
+def define_network_flows(Y, V):
+    coupling = V * (np.conjugate(Y.T * V)).T
+    coupling_sp = sps.csr_matrix(coupling)
+    data = coupling_sp.data
+    indptr = coupling_sp.indptr
+    indices = coupling_sp.indices
+
+    # noinspection PyPep8Naming
+    @njit(void(float64[:], float64[:]))
+    def network_flows(theta, P):
+        phases = np.exp(1.j * theta)
+        # In the next line the slice [:] is necessary to make sure the function writes into the view rather than
+        # replacing it.
+        P[...] += (phases * np.conjugate(numba_sp_dot(data, indptr, indices, phases))).real
+
+    return network_flows
+
+
+# noinspection PyPep8Naming
+def define_synchronous_machine(P_set, alpha, H):
+    # noinspection PyPep8Naming
+    @njit(void(complex128, complex128, complex128, complex128, float64[:], float64[:], float64[:]))
+    def synchronous_machine(P_complex, dP_complex, V_complex, dV_complex, res, omega, domega):
+        res_1 = dV_complex - 1.j * omega[0] * V_complex
+        res[0] = res_1.real
+        res[1] = res_1.imag
+        # Mathematically this keeps V magnitude fixed for this node. Of course we actually need to
+        res[2] = P_complex.real - (P_set - alpha * omega[0] - H * domega[0])
+
+    return synchronous_machine, 1
+
+
+class NodeType(object):
+    def __init__(self):
+        self.name = "passive"
+
+    def node_dynamics(self):
+        pass
+
+
+class SyncMachine(NodeType):
+    def __init__(self, P_set, alpha, H, V):
+        super(SyncMachine, self).__init__()
+        self.name = "SyncMachine"
+        self.P_set = P_set
+        self.alpha = alpha
+        self.H = H
+        self.V_magnitude = V
+
+    def node_dynamics(self):
+        pass
+
+def define_network_rhs(node_list, Y):
+    assert len(node_list) == len(Y)
+    length = len(Y)
+    total_length = 2*length
+
+    infeed = np.array(map(lambda x: x.P_set, node_list), dtype=np.float64)
+    damping = np.array(map(lambda x: x.alpha, node_list), dtype=np.float64)
+    V_mag_array = np.array(map(lambda x: x.V_magnitude, node_list), dtype=np.float64)
+    inv_inertia = np.array(map(lambda x: 1./x.H, node_list), dtype=np.float64)
+
+    infeed *= inv_inertia
+    damping *= inv_inertia
+    Y_over_inertia = np.diag(inv_inertia).dot(Y)
+
+    substract_n_flows = define_network_flows(Y_over_inertia, V_mag_array)
+
+    def network_rhs(y, t):
+        dydt = np.empty(total_length)
+        # For readability we define views into y and dydt:
+        phi = y[:length]
+        omega = y[length:]
+        dphi = dydt[:length]
+        domega = dydt[length:]
+
+        # This syntax makes sure we write into the underlying dydt rather than defining a new dphi variable.
+        dphi[...] = omega
+
+        domega[...] = infeed - damping * omega
+        substract_n_flows(phi, domega)
+
+        return dydt
+
+    return network_rhs
+
+if __name__ == "__main__":
+
+    node_list = list()
+    node_list.append(SyncMachine(1., 0.1, 1., 1.))
+    node_list.append(SyncMachine(-1., 0.1, 1., 1.))
+
+    Y = 8.j * np.ones((2, 2), dtype=np.complex128)
+    Y[0, 0] *= -1.
+    Y[1, 1] *= -1.
+
+    rhs = define_network_rhs(node_list, Y)
+    root_rhs = lambda x: rhs(x, 0.)
+
+    ic = np.ones(4, dtype=np.float64) * 0.5
+
+    from scipy.optimize import root
+
+    result = root(root_rhs, ic)
+    if result.success:
+        print(rhs(result.x, 0.))
+        ic = result.x
+    else:
+        print("failed")
+        exit()
+
+    times = np.linspace(0., 100., 1000)
+
+    from scipy.integrate import odeint
+
+    states = odeint(rhs, ic, times)
+
+    import matplotlib.pyplot as plt
+    plt.figure()
+    plt.plot(times, states[:,2:])
+    plt.show()
\ No newline at end of file

PSD/complex_current_and_nodes.py

+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+from numba import njit, float64, complex128, void, int32
+# from assimulo.solvers import CVODE
+# from assimulo.problem import Explicit_Problem
+import scipy.sparse as sps
+
+# The dynamics should be structured according to what is calculated on the network links.
+
+class NodeType(object):
+    def __init__(self, number_of_variables):
+        self.number_of_variables = number_of_variables
+
+
+class SwingEquationNode(NodeType):
+    def __init__(self, number_of_variables, v_set = 1., infeed = 1., damping = 0.1, H = 1.):
+        super(SwingEquationNode, self).__init__(number_of_variables)
+        self.v_set = v_set
+        self.v_set_squared = v_set ** 2
+        self.infeed = infeed
+        self.damping = damping
+        self.H_inv = 1./H
+
+    def node_dynamics(self, y, dy, i, t):
+        v = np.complex(y[0], y[1])
+        dv = 1.j * y[2] * v - (v.conjugate() * v - self.v_set_squared) * v
+        dy[0] = dv.real
+        dy[1] = dy.imag
+        dy[2] = self.H_inv * (self.infeed - self.damping * y[2] - (v * i.conjugate()).real)
+        return np.array
+
+
+@njit(complex128[:](complex128[:], int32[:], int32[:], complex128[:]))
+def numba_sp_dot(data, indptr, indices, v):
+    res = np.zeros(len(indptr) - 1, dtype=np.complex128)
+    index = 0
+    for row, number_of_entries in enumerate(indptr[1:]):
+        while index < number_of_entries:
+            res[row] += data[index] * v[indices[index]]
+            index += 1
+    return res
+
+@njit
+def numba_sp_complex_currents(data, indptr, indices, v, i):
+    index = 0
+    for row, number_of_entries in enumerate(indptr[1:]):
+        r_temp = 0.
+        while index < number_of_entries:
+            r_temp += data[index] * v[indices[index]]
+            index += 1
+        i[row] = r_temp
+
+
+def define_network_rhs(node_list, Y):
+    assert len(node_list) == len(Y)
+    length = len(Y)
+    total_length = 2*length
+
+    coupling_sp = sps.csr_matrix(Y)
+    data = coupling_sp.data
+    indptr = coupling_sp.indptr
+    indices = coupling_sp.indices
+
+    v_indices = [1, 2, 4, 5] # ...
+
+    def network_rhs(y, t):
+        dydt = np.empty(total_length)
+        i = coupling_sp.dot(v)
+
+        return dydt
+
+    return network_rhs
+
+if __name__ == "__main__":
+
+    node_list = list()
+    node_list.append(SwingEquationNode(3))
+    node_list.append(SwingEquationNode(3))
+
+    Y = 8.j * np.ones((2, 2), dtype=np.complex128)
+    Y[0, 0] *= -1.
+    Y[1, 1] *= -1.
+
+    rhs = define_network_rhs(node_list, Y)
+    root_rhs = lambda x: rhs(x, 0.)
+
+    ic = np.ones(4, dtype=np.float64) * 0.5
+
+    from scipy.optimize import root
+
+    result = root(root_rhs, ic)
+    if result.success:
+        print(rhs(result.x, 0.))
+        ic = result.x
+    else:
+        print("failed")
+        exit()
+
+    times = np.linspace(0., 100., 1000)
+
+    from scipy.integrate import odeint
+
+    states = odeint(rhs, ic, times)
+
+    import matplotlib.pyplot as plt
+    plt.figure()
+    plt.plot(times, states[:,2:])
+    plt.show()
\ No newline at end of file