using list with join

PSD/complex_current_and_nodes.py

@@ -29,7 +29,7 @@ class SwingEquationNode(NodeType):
     def node_dynamics_string(self, j="{index}"):
         return """
     dv[{index}] = 1.j * omega[{index}] * v[{index}] - (v[{index}].real * v[{index}].real + v[{index}].imag * v[{index}].imag - {v_set_squared}) * v[{index}]
-    domega[{index}] = {H_inv:.64f} * ({infeed} - {damping} * omega[{index}] - (v[{index}] * i[{index}].conjugate()).real)"""\
+    domega[{index}] = {H_inv:.64f} * ({infeed} - {damping} * omega[{index}] - (v[{index}] * i[{index}].conjugate()).real)""" \
             .format(
             v_set_squared=self.v_set_squared,
             H_inv=self.H_inv,
@@ -56,7 +56,7 @@ class DroopEquationNode(NodeType):
     def node_dynamics_string(self, j="{index}"):
         return """
     dv[{index}] = 1.j * omega[{index}] * v[{index}] - (v[{index}].real * v[{index}].real + v[{index}].imag * v[{index}].imag - {v_set_squared}) * v[{index}]
-    domega[{index}] = {H_inv:.64f} * ({infeed} - {damping} * omega[{index}] - (v[{index}] * i[{index}].conjugate()).real)"""\
+    domega[{index}] = {H_inv:.64f} * ({infeed} - {damping} * omega[{index}] - (v[{index}] * i[{index}].conjugate()).real)""" \
             .format(
             v_set_squared=self.v_set_squared,
             H_inv=self.H_inv,
@@ -154,30 +154,29 @@ def define_network_rhs_codegen(node_list, Y):
     indices = coupling_sp.indices
 
     def_network_rhs_string = """
-@njit(float64[:](float64[:], float64))
-def network_rhs_numba(y, t):
-    dydt = np.empty(total_length + length)
+    @njit(float64[:](float64[:], float64))
+    def network_rhs_numba(y, t):
+        dydt = np.empty(total_length + length)
 
-    v = y[:total_length].view(np.complex128)
-    omega = y[total_length:]
-    dv = dydt[:total_length].view(np.complex128)
-    domega = dydt[total_length:]
+        v = y[:total_length].view(np.complex128)
+        omega = y[total_length:]
+        dv = dydt[:total_length].view(np.complex128)
+        domega = dydt[total_length:]
 
-    i = np.zeros(l_indptr - 1, dtype=np.complex128)
-    index = 0
-    for row, number_of_entries in enumerate(indptr[1:]):
-        while index < number_of_entries:
-            i[row] += data[index] * v[indices[index]]
-            index += 1
+        i = np.zeros(l_indptr - 1, dtype=np.complex128)
+        index = 0
+        for row, number_of_entries in enumerate(indptr[1:]):
+            while index < number_of_entries:
+                i[row] += data[index] * v[indices[index]]
+                index += 1
 
-"""
+    """
 
-    for j, node in enumerate(node_list):
-        def_network_rhs_string += node.node_dynamics_string(j)
+    def_network_rhs_string = ''.join([node.node_dynamics_string(j) for j, node in enumerate(node_list)])
 
     def_network_rhs_string += """
     return dydt
-"""
+    """
 
     context = globals()
     context.update(locals())

PSD/microgrid_frequency_dynamics.py

@@ -26,7 +26,7 @@ def define_gen_rc(brp, rhs):
         result = root(root_rhs, ic)
         if result.success:
             ic = result.x
-            ic[2 * system_size + batch+1] += 0.1
+            ic[2 * system_size + batch + 1] += .1 * (1. - 2. * np.random.random())
         else:
             print("failed")
         return ic, ()
@@ -39,7 +39,7 @@ 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):
+def main(sim_dir=default_dir, create_test_data=True, run_test=True, flag_baobab=True):
 
     node_list, Y = load_PYPSA('microgrid_testcase.npz')
     rhs = define_network_rhs(node_list, Y)
@@ -53,7 +53,7 @@ def main(sim_dir=default_dir, create_test_data=True, run_test=True,flag_baobab=T
         times = np.linspace(0, 10, 100)
         brp = bao.BatchRunParameters(number_of_batches=10, simulations_per_batch=10, system_dimension=100)
         ob = bao.combine_obs(bao.run_condition_observer, lorenz_ob)
-        bao.run_experiment(result_file, brp, None, define_gen_rc(brp,rhs), ob, times, figures=(True, 1), rhs=rhs)
+        bao.run_experiment(result_file, brp, None, define_gen_rc(brp, rhs), ob, times, figures=(True, 1), rhs=rhs)
         res = bao.load_all_fields_from_results(result_file)
 
         if bao.am_master:
@@ -73,8 +73,8 @@ def main(sim_dir=default_dir, create_test_data=True, run_test=True,flag_baobab=T
         from scipy.integrate import odeint
         brp = bao.BatchRunParameters(number_of_batches=10, simulations_per_batch=10, system_dimension=100)
 
-        generate_run_conditions = define_gen_rc(brp,rhs)
-        rc = generate_run_conditions(0,0)
+        generate_run_conditions = define_gen_rc(brp, rhs)
+        rc = generate_run_conditions(0, 0)
 
         #node_list[0].infeed += 0.1
         rhs = define_network_rhs(node_list, Y)
@@ -82,14 +82,15 @@ def main(sim_dir=default_dir, create_test_data=True, run_test=True,flag_baobab=T
         import matplotlib.pyplot as plt
         from mpl_toolkits.mplot3d import Axes3D
 
-        plt.figure()
-        plt.plot(times, states[:, 2*brp.system_dimension:])
         fig = plt.figure()
         ax = fig.gca(projection='3d')
-        ax.plot(states[:, 0], states[:, 1], times)
-        ax.plot(states[:, 2], states[:, 3], times)
+        for k in range(brp.system_dimension):
+            ax.plot(states[:, k], states[:, brp.system_dimension+k], states[:, 2*brp.system_dimension+k])
+        ax.set_xlabel("real V")
+        ax.set_ylabel("imag V")
+        ax.set_zlabel(r"$\omega$")
         plt.show()
 
 
 if __name__ == "__main__":
-    main()
+    main(flag_baobab=True)