Merge remote-tracking branch 'origin/jonathan' into jonathan

# Conflicts:
#	PSD/complex_current_and_nodes.py

PSD/__init__.py

-from complex_current_and_nodes import *
-
-__all__ = ['complex_current_and_nodes']
\ No newline at end of file
+# from complex_current_and_nodes import *
+#
+# __all__ = ['complex_current_and_nodes_new']

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_test.py

+import numpy as np
+import complex_current_and_nodes as ccn
+from scipy.optimize import root
+
+
+def approx_jacobian(x, f, epsilon=0.01):
+    identity = np.diag(np.ones_like(x))
+    # This makes use of numpy array broadcasting rules, do not expect to understand it straight away.
+    return (np.apply_along_axis(f,1,epsilon * identity + x) - f(x))/epsilon
+
+# This assumes we are using a node type with exactly one internal variable = omega.
+# The lines where this is implicitly assumed are marked with # omega
+if __name__ == "__main__":
+
+    # create network nodes
+    nodeType = ccn.EnumNodeType.SWING
+    node_list = list()
+    node_list.append(ccn.SwingEquationNode(infeed= 1))
+    node_list.append(ccn.SwingEquationNode(infeed= 1))
+    node_list.append(ccn.SwingEquationNode(infeed= 1))
+    node_list.append(ccn.SwingEquationNode(infeed=-1))
+    node_list.append(ccn.SwingEquationNode(infeed=-1))
+    node_list.append(ccn.SwingEquationNode(infeed=-1))
+
+    # create network vertices - lossless network
+    # sum along a row/column must be 0!
+    # ToDo: document in .tex
+    # coupling of 8 is chosen at random but well above minimum for system stability for 3 node system with total power
+    # input/output of /pm 2.
+    if len(node_list) > 3:
+        assert (len(node_list)%2 == 0)
+        import networkx as nx
+        #
+        graph = nx.random_regular_graph(5, len(node_list), np.random.rand(1)[0])
+        Y = np.zeros((len(node_list), len(node_list)), dtype=np.complex128)
+        Y += nx.to_numpy_matrix(graph, dtype=np.complex128) * -8 * 1.j
+        for j in range(len(node_list)):
+            Y[j, j] = -1 * Y[:, j].sum(0)
+    else:
+        Y = np.ones((len(node_list), len(node_list)), dtype=np.complex128) * -8 * 1.j
+        for i in range(len(node_list)):
+            Y[i, i] *= -1
+
+    # define system of gde
+    rhs = ccn.define_network_rhs(node_list,Y)
+    # define system of algebraic equations for usage in root (fixed point) finding algorithm
+    root_rhs = ccn.define_root_rhs_omega(len(node_list), rhs) # omega
+    # find fixed point using scipy.root
+    # Attention! For these systems the .success flag fails to accurately reflect the algorithms success use instead:
+    # the indicator used here.
+    guess = np.random.rand(len(node_list) * 3 + 1)            # omega
+    guess = np.ones(len(node_list) * 3 + 1, dtype=np.float64) * 0.5
+    root_result = root(root_rhs,guess)
+
+    # print the result quantitatively and qualitatively
+    # quantitatively
+    if np.all(root_rhs(root_result.x) < 1e-6):
+        print("\nFound a fixed point of the system! At:\n")
+        print(root_result.x)
+    else:
+        print("\nNO fixed point found!\n")
+    print("\nrhs(result_result.x[:-1], 0) =\n{}".format(rhs(root_result.x[:-1], 0)))
+    # qualitatively
+    jacobian = approx_jacobian(root_result.x[:-1], lambda x: rhs(x, 0))
+    w, _ = np.linalg.eig(jacobian)
+    print("\nThe eigenvalues of the jacobian are:\n{}".format(w[::-1]))
+
+    # perturb system to
+    # - check stability of fixed point (i.e. see if it will return to fixed point)
+    # - check dynamic behaviour of Nodes (i.e. check for error in calculations or bug in code)
+    perturbation = (np.random.rand(len(node_list) * 3) - 0.5 ) / 100 * 0
+    print('\nCausing perturbation:\n{}'.format(perturbation))
+
+    # solve system dynamics
+    from scipy.integrate import odeint
+
+    times = np.linspace(0., 100., 100000)
+    states = odeint(rhs, root_result.x[:-1] + perturbation, times)
+
+    # plot results
+    import matplotlib.pyplot as plt
+    from mpl_toolkits.mplot3d import Axes3D
+
+    plt.figure()
+    plt.plot(times, states[:, len(node_list) * 2:])
+    plt.xlabel('time')
+    plt.ylabel('frequency')
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+    for i in range(len(node_list)):
+        ax.plot(states[:, i], states[:, i + 1], times)
+    plt.ylabel('voltage')
+    fig = plt.figure()
+    for i in range(len(node_list)):
+        plt.plot(times, (states[:, i] ** 2 + states[:, i + 1] ** 2) ** .5)
+    plt.ylabel('voltage')
+    plt.show()
+
+

PSD/complex_current_and_nodes_test_fixed_point_finding.py

+import numpy as np
+from scipy.optimize import root
+from scipy.integrate import odeint
+import complex_current_and_nodes as ccn
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+def approx_jac(x, f, epsilon=0.01):
+    identity = np.diag(np.ones_like(x))
+    # This makes use of numpy array broadcasting rules, do not expect to understand it straight away.
+    return (np.apply_along_axis(f,1,epsilon * identity + x) - f(x))/epsilon
+
+
+def test_fp_finding(no_of_runs, node_list, Y):
+    rhs = ccn.define_network_rhs(node_list, Y)
+    root_rhs = ccn.define_root_rhs_omega(2 * len(Y), rhs)
+
+    xRange = np.arange(0, len(node_list) * 3 + 1, 1)
+    yValues = np.zeros((len(node_list) * 3 + 1, 2), dtype=np.float64)
+
+    for i in range(no_of_runs):
+        guess = np.random.rand(len(node_list) * 3 + 1)
+        result = root(root_rhs, guess)
+
+        jac = approx_jac(result.x[:-1], lambda x: rhs(x, 0.))
+        w, _ = np.linalg.eig(jac)
+
+        no_pos_ev = 0
+        for eigen in w:
+            if eigen > 0:
+                no_pos_ev += 1
+
+        if result.success:
+            yValues[no_pos_ev, 0] += 1
+        else:
+            yValues[no_pos_ev, 1] += 1
+
+    fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, sharey=True, figsize=(6, 6))
+    if node_list[0].nodeType == ccn.EnumNodeType.SWING:
+        fig.suptitle('{} node network; swingeq.; fully interconnected'.format(len(node_list)))
+    elif node_list[0].nodeType == ccn.EnumNodeType.DROOP:
+        fig.suptitle('{}-node network; volt. droopeq.; fully interconnected'.format(len(node_list)))
+
+    fig.text(0.5, 0.04, 'no. of positive eigenvalues of jacobian', ha='center')
+
+    plt.subplot(211)
+    if yValues[:,0].sum(0) == 0:
+        plt.gca().set_ylim([0,1])
+    plt.ylabel('N out of {}'.format(int(yValues[:,0].sum(0))))
+    fig.text(2. * len(yValues[:,0]) / 3., 2. * yValues[:,0].max(0) / 3., 'root search success')
+    plt.bar(xRange, yValues[:,0], 1./2, color='green')
+
+    plt.subplot(212)
+    if yValues[:,1].sum(0) == 0:
+        plt.gca().set_ylim([0,1])
+    plt.ylabel('N out of {}'.format(int(yValues[:,1].sum(0))))
+    fig.text(2. * len(yValues[:,0]) / 3., 2. * yValues[:, 1].max(0) / 3., 'root search failures')
+    plt.bar(xRange, yValues[:,1], 1./2, color='red')
+
+    plt.show()
+
+def test_fp_finding_random_network(no_of_runs, no_of_feeding_nodes, no_of_using_nodes, node_type):
+    no_of_nodes = no_of_feeding_nodes + no_of_using_nodes
+    assert (no_of_nodes%2 == 0)
+
+    success = 0.
+    for i in range(no_of_runs):
+        # create new power distribution
+        power_input = np.random.rand(no_of_feeding_nodes)
+        power_input /= power_input.sum(0)
+        power_usage = np.random.rand(no_of_using_nodes) * -1
+        power_usage /= np.absolute(power_usage.sum(0))
+        power_input = np.concatenate((power_input, power_usage))
+        np.set_printoptions(precision=3)
+
+        # take power distribution and make node_list
+        node_list = list()
+        if node_type == ccn.EnumNodeType.SWING:
+            for j in range(len(power_input)):
+                node_list.append(ccn.SwingEquationNode(3, power_input[j]))
+        elif node_type == ccn.EnumNodeType.DROOP:
+            for j in range(len(power_input)):
+                node_list.append(ccn.DroopEquationNode(3, power_input[j]))
+
+        # create network admittance matrix
+        if no_of_nodes > 3:
+            import networkx as nx
+            #
+            graph = nx.random_regular_graph(3, len(node_list), np.random.rand(1)[0])
+            Y = np.zeros((len(node_list), len(node_list)), dtype=np.complex128)
+            Y += nx.to_numpy_matrix(graph, dtype=np.complex128) * -8 * 1.j
+            for j in range(len(node_list)):
+                Y[j, j] = -1 * Y[:, j].sum(0)
+        else:
+            Y = np.ones((2, 2), dtype=np.complex128) * -8 * 1.j
+            Y[0, 0] *= -1
+            Y[1, 1] *= -1
+
+        def root_rhs_2(y, phase_coupling=Y, input_power=power_input):
+            phases = np.exp(1.j * y)
+            return input_power - np.real(phases * phase_coupling.dot(phases).conjugate())
+
+
+        # print exemplary network
+        if i == 0:
+            print("sum({})=\n{}".format(power_input, power_input.sum(0)))
+            print(Y)
+        #
+        # create function determining the network dynamics
+        rhs = ccn.define_network_rhs(node_list, Y)
+        # look for fixed points
+        root_rhs = ccn.define_root_rhs_omega(node_list, rhs)
+        guess = np.random.rand(len(node_list) * 3 + 1)
+        # root_rhs = root_rhs_2
+        # guess = np.zeros(len(Y))
+
+        we_are_done_here = False
+        it_counter = 0
+        while not we_are_done_here:
+            if it_counter > 1:
+                print("it_counter={}, done={}".format(it_counter, we_are_done_here))
+            try:
+                result = root(root_rhs, guess)#, options={'maxiter': int(1e7)})#, method="Krylov", tol = 10e-6, )
+            except ValueError, e:
+                if str(e) == "Jacobian inversion yielded zero vector. This indicates a bug in the Jacobian approximation.":
+                    print("            Bad initial guess. Try again")
+                    continue
+                else:
+                    raise e
+
+            result.success = np.all(root_rhs(result.x) < 1e-3)
+            we_are_done_here = True
+
+            if result.message == "The iteration is not making good progress, as measured by the \
+  improvement from the last five Jacobian evaluations.":
+                print("#############")
+                we_are_done_here = False
+                it_counter += 1
+                guess = np.random.rand(len(node_list) * 3 + 1)
+
+
+        if result.success:
+            print("success    sum(root_rhs(result.x))/no_of_nodes = {:9.6f}".format(sum(np.abs(root_rhs(result.x)))/no_of_nodes))
+            success += 1.
+        else:
+            print("failure    sum(root_rhs(result.x))/no_of_nodes = {:9.6f}".format(sum(np.abs(root_rhs(result.x)))/no_of_nodes))
+            print(result.message)
+    print("success rate = {}".format(success/no_of_runs))
+    return success / (no_of_runs)
+
+
+if __name__ == "__main__":
+
+    node_type = ccn.EnumNodeType.DROOP
+    no_of_max_nodes = 40
+    no_of_runs = 20
+    fp_success_rate_array = np.zeros((no_of_max_nodes, no_of_max_nodes), dtype=np.float64)
+    # avg_neg_jac_success = np.zeros((no_of_max_nodes, no_of_max_nodes), dtype=np.float64)
+    # avg_neg_jac_failure = np.zeros((no_of_max_nodes, no_of_max_nodes), dtype=np.float64)
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    ax.set_xlabel('no. of generating nodes')
+    ax.set_ylabel('no. of using nodes')
+    ax.set_zlabel('fp finding success rate')
+
+    fig_jac_s = plt.figure()
+    ax_jac_s = fig_jac_s.add_subplot(111, projection='3d')
+    ax_jac_s.set_xlabel('no. of generating nodes')
+    ax_jac_s.set_ylabel('no. of using nodes')
+    ax_jac_s.set_zlabel('avg no. of neg jac EV - fp success')
+
+    fig_jac_f = plt.figure()
+    ax_jac_f = fig_jac_f.add_subplot(111, projection='3d')
+    ax_jac_f.set_xlabel('no. of generating nodes')
+    ax_jac_f.set_ylabel('no. of using nodes')
+    ax_jac_f.set_zlabel('avg no. of neg jac EV - fp failure')
+
+
+    nodes = 2
+    while nodes <= no_of_max_nodes:
+        producer = 1
+        while producer < nodes:
+            user = nodes - producer
+            fp_success_rate = test_fp_finding_random_network(
+                                                                    no_of_runs, producer, user, node_type)
+            fp_success_rate_array[producer, user] = fp_success_rate
+
+            plt.figure(fig.number)
+            ax.bar3d(producer-0.5, user-0.5, 0, 1, 1, fp_success_rate, color='green')
+            plt.figure(fig_jac_s.number)
+            # ax_jac_s.bar3d(producer-0.5, user-0.5, 0, 1, 1, avg_neg_jac_success, color='green')
+            # plt.figure(fig_jac_f.number)
+            # ax_jac_f.bar3d(producer - 0.5, user - 0.5, 0, 1, 1, avg_neg_jac_failure, color='red')
+            producer += 1
+        nodes += 2
+
+    print(fp_success_rate_array)
+
+    plt.show()
\ No newline at end of file

PSD/lib.py

+# the tutorial on how to create sphinx told me to create this file. Not sure what it does, or wether it does anything. \
+# But keeping it to be on the safe side.
\ No newline at end of file

PSD/microgrid_droopEq_dynamics.py

-#from __future__ import division, print_function, absolute_import, unicode_literals
-
-import numpy as np
-from numba import njit
-import os
-from complex_current_and_nodes import DroopEquationNode, SwingEquationNode, define_network_rhs, define_network_rhs_codegen, define_root_rhs
-
-try:
-    import baobap as bao
-except:
-    ImportWarning("Cannot find pyBAOBAP")
-
-
-default_dir = os.path.join("simulation_data", "microgrid")
-
-def load_PYPSA(filename):
-    """
-    Loads the specified dataset and returns a (possibly sparse) matrix determining the network connections and their conductances, as well as a node_list of appropriate length. The nodes are of type DroopEquationNodes.
-    :param filename: filename of the dataset
-    :return: node_list, conductance matrix
-    """
-    #TODO replace constants with proper values. Attention: k_P,i != k_P,j & k_Q,i != k_Q,j - !!!remember: S^N = 1!!!
-    data = np.load(filename)
-    input_power = data['input_power']
-    for p in input_power:
-        print "{}".format(p)
-    volt_droop_gain = 0.2 * input_power
-    Y = data['Y']
-    system_size = len(input_power)
-    node_list = [DroopEquationNode(3, infeed=input_power[i], volt_droop_gain=0.1*input_power[i], freq_droop_gain=0.2*input_power[i]) for i in range(system_size)]
-    return node_list,Y
-
-
-def load_PyPSA_df(adm, par):
-    #TODO: documentation
-    """
-
-    :param adm:
-    :param par:
-    :return:
-    """
-    from pandas import read_csv
-    Y = np.load(adm)
-    Y = Y[()] # this is import to recover a csr_matrix
-    df = read_csv(par, index_col=0)
-    system_size = df.shape[0]
-    node_list = [DroopEquationNode(3, infeed=input_power[i], volt_droop_gain=0.1*input_power[i], freq_droop_gain=0.2*input_power[i]) for i in range(system_size)]
-    load_flow_sol = np.zeros(3 * system_size, dtype=np.float64)
-    load_flow_sol[:2*system_size] = np.array(df.v_pu.values * np.exp(1.j * df.phi.values), dtype=np.complex128).view(np.float64)
-    return node_list, Y, load_flow_sol
-
-def define_gen_rc(brp, rhs, init=None, method="krylov"):
-    """
-
-    :param brp: Variable holding the baobap batch run parameters(see: baobap.BatchRunParameters(...)).
-    :param rhs: Function handle of function determining the behaviour of the network nodes.
-    :param init: A vector of length 3 * brb.system_dimension, containing the initial guesses for the root finding algorithim.
-    :param method: String detailing which root finding method to use. Must be one of those implemented in scipy.optimize.root(...).
-    :return: returns a function handle to the function generate_run_conditions(batch, run)
-    """
-    system_size = brp.system_dimension
-    if init is None:
-        init = np.ones(3 * system_size, dtype=np.float64) * 0.5
-
-    def generate_run_conditions(batch, run):
-        """
-Takes returns a vector of length
-        :param batch:
-        :param run: This seems to ne non-essential.
-        :return:
-        """
-        root_rhs = define_root_rhs(2 * system_size, rhs)
-        from scipy.optimize import root
-        print "start root"
-        result = root(root_rhs, init, method="{}".format(method))
-        print "end root"
-        if result.success:
-            ic = result.x
-            ic[2 * system_size + batch + 1] += .1 * (1. - 2. * np.random.random())
-        else:
-            print("failed")
-            ic = init
-        return ic, ()
-    return generate_run_conditions
-
-
-# This function observes the results of the simulation run and returns its observations in a dict
-# noinspection PyUnusedLocal
-def lorenz_ob(time_series, rc):
-    return{"max_values": np.max(time_series, axis=0)}
-
-
-def main(sim_dir=default_dir, create_test_data=True, run_test=True, flag_baobab=True):
-    import time
-    print "start", "{0.tm_year}-{0.tm_mon}.{0.tm_mday}.-{0.tm_hour}h{0.tm_min}m.timestamp".format(time.localtime(time.time()))
-
-    node_list, Y = load_PYPSA('microgrid_testcase.npz')
-    load_flow_sol = None
-    # node_list, Y, load_flow_sol = load_PyPSA_df("bus_admittance.npy", "bus_parameters")
-    # rhs = define_network_rhs(node_list, Y)
-    rhs = define_network_rhs_codegen(node_list, Y)
-
-    print "compilation finished", "{0.tm_year}-{0.tm_mon}.{0.tm_mday}.-{0.tm_hour}h{0.tm_min}m.timestamp".format(time.localtime(time.time()))
-
-    if flag_baobab:
-        # Run PYPSD code with BAOBAB
-        np.random.seed(0)
-
-        result_file, batch_dir = bao.prep_dir(sim_dir)
-        times = np.linspace(0, 10, 100)