Commit 71cec808 authored by Paul Schultz's avatar Paul Schultz
Browse files

corrected numba code

parent 3198cfc4
......@@ -28,8 +28,8 @@ 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)""" \
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)""" \
.format(
v_set_squared=self.v_set_squared,
H_inv=self.H_inv,
......@@ -55,8 +55,8 @@ 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)""" \
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)""" \
.format(
v_set_squared=self.v_set_squared,
H_inv=self.H_inv,
......@@ -110,9 +110,10 @@ def define_network_rhs(node_list, Y):
return network_rhs
def define_root_rhs(total_length, rhs):
def root_rhs(y):
dydt = rhs(y, 0.)
def define_root_rhs(total_length, func):
rhs = lambda x: func(x, 0)
def root_rhs(y, rhs=rhs):
dydt = rhs(y)
v = y[:total_length].view(np.complex128)
dv = dydt[:total_length].view(np.complex128)
domega = dydt[total_length:]
......@@ -142,8 +143,8 @@ def define_root_rhs_omega(total_length, rhs):
def define_network_rhs_codegen(node_list, Y):
network_rhs_numba = None
assert len(node_list) == len(Y)
length = len(Y)
assert len(node_list) == Y.shape[0]
length = Y.shape[0]
total_length = 2*length
coupling_sp = sps.csr_matrix(Y)
......@@ -154,21 +155,21 @@ def define_network_rhs_codegen(node_list, Y):
indices = coupling_sp.indices
def_network_rhs_string = """
@njit(float64[:](float64[:], float64), cache=True)
def network_rhs_numba(y, t):
dydt = np.empty(total_length + length)
@njit(float64[:](float64[:], float64), cache=True)
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
"""
......
......@@ -24,15 +24,22 @@ def load_PyPSA_df(adm, par):
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)]
return node_list, Y
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"):
def define_gen_rc(brp, rhs):
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):
ic = np.ones(3 * system_size, dtype=np.float64) * 0.5
root_rhs = define_root_rhs(2 * system_size, rhs)
from scipy.optimize import root
result = root(root_rhs, ic)
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())
......@@ -49,11 +56,15 @@ def lorenz_ob(time_series, rc):
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')
node_list, Y = load_PyPSA_df("bus_admittance.npy", "bus_parameters")
rhs = define_network_rhs(node_list, Y)
# nrhs = define_network_rhs_codegen(node_list, Y)
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
......@@ -84,10 +95,7 @@ def main(sim_dir=default_dir, create_test_data=True, run_test=True, flag_baobab=
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)
generate_run_conditions = define_gen_rc(brp, rhs, init=load_flow_sol, method="hybr")
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()))
......
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