Commit 166c0e81 authored by Paul Schultz's avatar Paul Schultz
Browse files

debugging

parent 46b03aae
......@@ -88,9 +88,9 @@ def numba_sp_complex_currents(data, indptr, indices, v, i):
def define_network_rhs(node_list, Y):
assert len(node_list) == len(Y)
length = len(Y)
total_length = 2*length
assert len(node_list) == Y.shape[0]
length = Y.shape[0]
total_length = 2 * length
def network_rhs(y, t):
coupling_sp = sps.csr_matrix(Y)
......
......@@ -20,6 +20,7 @@ def load_PYPSA(filename):
def load_PyPSA_df(adm, par):
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 = [SwingEquationNode(3, infeed=df.p.iloc[i], H=0.1, damping=0.01, v_set=df.v_nom.iloc[i]) for i in range(system_size)]
......@@ -80,20 +81,29 @@ def main(sim_dir=default_dir, create_test_data=True, run_test=True, flag_baobab=
# Run Test PYPSD without BAOBAB
times = np.linspace(0., 10., 100)
from scipy.integrate import odeint
brp = bao.BatchRunParameters(number_of_batches=10, simulations_per_batch=10, system_dimension=100)
brp = bao.BatchRunParameters(number_of_batches=1, simulations_per_batch=1, system_dimension=Y.shape[0])
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()))
generate_run_conditions = define_gen_rc(brp, rhs)
rc = generate_run_conditions(0, 0)
print "generate_run_conditions", "{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[0].infeed += 0.1
rhs = define_network_rhs(node_list, Y)
print "define_network_rhs", "{0.tm_year}-{0.tm_mon}.{0.tm_mday}.-{0.tm_hour}h{0.tm_min}m.timestamp".format(time.localtime(time.time()))
states = odeint(rhs, rc[0], times)
print "integration", "{0.tm_year}-{0.tm_mon}.{0.tm_mday}.-{0.tm_hour}h{0.tm_min}m.timestamp".format(time.localtime(time.time()))
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
for k in range(brp.system_dimension):
for k in range(3):
ax.plot(states[:, k], states[:, brp.system_dimension+k], states[:, 2*brp.system_dimension+k])
ax.set_xlabel(r"$\Re$ V")
ax.set_ylabel(r"$\Im$ V")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment