Commit ac0e07cc authored by Jonathan_H's avatar Jonathan_H
Browse files

immernoch probleme mit droopequationnode, aber soweit ich festellen konnte...

immernoch probleme mit droopequationnode, aber soweit ich festellen konnte keine Vorzeichenfehler in DroopEquationNode class
parent 47e960b8
This diff is collapsed.
......@@ -51,22 +51,22 @@ class SwingEquationNode(NodeType):
# low_pass_time_const is taken from (Schiffer, Ortega et al, 2014) p.2466
class DroopEquationNode(NodeType):
"""This is a docstring."""
def __init__(self, number_of_variables, v_set = 1., infeed = 1., low_pass_time_const = 0.5, volt_droop_gain = 0.2, freq_droop_gain = 0.4):
def __init__(self, number_of_variables, infeed = 1., v_set = 1., low_pass_time_const = 0.5, volt_droop_gain = 0.2, freq_droop_gain = 0.4):
super(DroopEquationNode, self).__init__(number_of_variables)
self.omega_set = 0
self.v_set = v_set
self.v_set_squared = v_set ** 2
self.infeed = infeed
self.t_inv = 1/low_pass_time_const
self.volt_droop_gain = volt_droop_gain # k_P in (Schiffer, Ortega et al, 2014)
self.freq_droop_gain = freq_droop_gain # k_Q in (Schiffer, Ortega et al, 2014)
self.volt_droop_gain = volt_droop_gain # k_Q in (Schiffer, Ortega et al, 2014)
self.freq_droop_gain = freq_droop_gain # k_P in (Schiffer, Ortega et al, 2014)
# in Schiffer, Ortega et al, 2014 werden k_P und k_Q wie folgt bestimmt:
# k_P = 0.1/S_(N,i) [Hz/pu]
# k_Q = 0.2/S_(N,i) [pu/pu]
# S_(N,i) entspricht dem power rating des nodes in pu - heißt: der hier eingespeisten Energie in pu des Gesamtsystems.
# um eine realistische Größe abzuschätzen muss in diesem Scope bekannt sein wie groß das Netzwerk ist.
#TODO Die gesamtgröße des Netzwerks muss in in diesen Scope, oder sogar schon beim Aufrufen der Funktion node_dynamics bekannt sein
#TODO ist infeed in diesem Fall auch noch konstant?
#done: ist infeed in diesem Fall auch noch konstant? - ja, self.infeed ist infeed desired
def node_dynamics(self, v, omega, i, t, dv, domega, index):
"""
# The current is assumed to be negative when it is flowing away from the node
......@@ -90,10 +90,12 @@ class DroopEquationNode(NodeType):
def node_dynamics_string(self, j="{index}"):
return """
if np.absolute(v[{index}]) == 0:
v_abs = np.absolute(v[{index}])
if v_abs == 0:
dv[{index}] = 0
else:
dv[{index}] = 1.j * v[{index}] * omega[{index}] + {t_inv} * (v[{index}] / np.absolute(v[{index}])) * ( ({v_set} - np.absolute(v[{index}])) - {volt_droop_gain} * ( (v[{index}] * i[{index}].conjugate()).imag - {infeed}.imag ) )
dv[{index}] = 1.j * omega[{index}] * v[{index}] + {t_inv} * (v[{index}] / v_abs) * ( ({v_set} - v_abs)\
- {volt_droop_gain} * ( (v[{index}] * i[{index}].conjugate()).imag - {infeed}.imag ) )
domega[{index}] = {t_inv} * ({omega_set} - omega[{index}]) - {t_inv} * {freq_droop_gain} * ( ( v[{index}] * i[{index}].conjugate()).real - {infeed}.real )""" \
.format(
......@@ -102,7 +104,6 @@ class DroopEquationNode(NodeType):
volt_droop_gain=self.volt_droop_gain,
freq_droop_gain=self.freq_droop_gain,
omega_set=self.omega_set,
v_set_squared=self.v_set_squared,
infeed=self.infeed,
index=j)
......@@ -266,10 +267,9 @@ def approx_jac(x, f, epsilon=0.01):
if __name__ == "__main__":
node_list = list()
node_list.append(DroopEquationNode(3, infeed=1.))
node_list.append(DroopEquationNode(3, infeed=-1.))
# node_list.append(DroopEquationNode(3, infeed=1.))
# node_list.append(DroopEquationNode(3, infeed=-3.))
node_list.append(DroopEquationNode(3, infeed=-1 - 0.000j ))
node_list.append(DroopEquationNode(3, infeed= 1 + 0.000j))
#node_list.append(DroopEquationNode(3, infeed=-3.))
Y = -8.j * np.ones((2, 2), dtype=np.complex128)
Y[0, 0] *= -1.
......@@ -289,7 +289,24 @@ if __name__ == "__main__":
root_rhs = define_root_rhs_omega(2*len(Y), rhs)
guess = 1.* np.ones(7, dtype=np.float64)
# guess = 1.* np.ones(7, dtype=np.float64)
# desired omega, desired voltage and angle difference
omega_des = node_list[0].omega_set
v_des = node_list[0].v_set
phi_diff = 2*np.pi + 0*2.000000000000000*np.pi
guess = np.zeros(7)
#voltage guess node one (real & imag)
#amplitude = v_set & def: phi1 = 0
guess[3] = v_des #+ 0.000000000000001
guess[4] = 0
#voltage guess node two (real & imag)
#calculated from phi1 - phi2 = pi/2
guess[5] = v_des*np.cos(phi_diff)
guess[6] = v_des*np.sin(phi_diff)
#omega guess where omega1 = omega2 = omega_global
guess[1] = omega_des + 1/2 * node_list[0].volt_droop_gain * ( node_list[0].infeed.real + node_list[1].infeed.real)
guess[2] = guess[1]
# print(root_rhs(guess))
......@@ -321,8 +338,19 @@ if __name__ == "__main__":
plt.figure()
plt.plot(times, states[:, 4:])
plt.ylabel('frequency')
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(states[:, 0], states[:, 1], times)
ax.plot(states[:, 2], states[:, 3], times)
plt.ylabel('voltage')
fig = plt.figure()
#plt.plot(times,states[:, 0])
plt.plot(times, (states[:, 0]**2+states[:, 1]**2)**.5)
plt.hold('on')
plt.plot(times, (states[:, 2]** 2 + states[:, 3]** 2) ** .5)
#lt.plot(times,states[:, 2], times, states[:, 3])
#plt.ylim([0,200])
plt.hold('off')
plt.ylabel('voltage')
plt.show()
#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)
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)
res = bao.load_all_fields_from_results(result_file)
if bao.am_master:
if create_test_data:
for key in res.keys():
np.save(os.path.join(sim_dir, key), res[key])
if run_test:
for key in res.keys():
filename = key + ".npy"
with open(os.path.join(sim_dir, filename), mode='rb') as f:
arr = np.load(f)
assert np.allclose(arr, res[key])
else:
# Run Test PYPSD without BAOBAB
times = np.linspace(0., 10., 100)
from scipy.integrate import odeint
brp = bao.BatchRunParameters(number_of_batches=1, simulations_per_batch=1, system_dimension=Y.shape[0])
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()))
#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(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")
ax.set_zlabel(r"$\omega$")
plt.show()
if __name__ == "__main__":
main(flag_baobab=False)
......@@ -14,6 +14,11 @@ except:
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 SwingEquationNode.
:param filename: filename of the dataset
:return: node_list, conductance matrix
"""
data = np.load(filename)
input_power = data['input_power']
Y = data['Y']
......@@ -21,7 +26,15 @@ def load_PYPSA(filename):
node_list = [SwingEquationNode(3, infeed=input_power[i], H=0.1, damping=0.01) 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
......@@ -33,12 +46,25 @@ def load_PyPSA_df(adm, par):
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: I do not know what this is
:return:
"""
root_rhs = define_root_rhs(2 * system_size, rhs)
from scipy.optimize import root
print "start root"
......
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