diff --git a/.gitignore b/.gitignore index 5fc9f9c..37f56d9 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ /examples/models/custom_multi_machine/data/plots /examples/models/ieee_9bus/data/plots /examples/models/ibb_model/data/plots +/examples/models/k2a_manual/data/plots diff --git a/examples/models/custom_multi_machine/custom_multi_machine_opt.py b/examples/models/custom_multi_machine/custom_multi_machine_opt.py index 6832a28..745daea 100644 --- a/examples/models/custom_multi_machine/custom_multi_machine_opt.py +++ b/examples/models/custom_multi_machine/custom_multi_machine_opt.py @@ -5,8 +5,8 @@ import torch import examples.models.custom_multi_machine.custom_multi_machine_model as mdl -from src.diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation +from diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation np.random.seed(0) diff --git a/examples/models/custom_multi_machine/custom_multi_machine_sim.py b/examples/models/custom_multi_machine/custom_multi_machine_sim.py index 10f94cd..48f53c5 100644 --- a/examples/models/custom_multi_machine/custom_multi_machine_sim.py +++ b/examples/models/custom_multi_machine/custom_multi_machine_sim.py @@ -4,7 +4,7 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.custom_multi_machine.custom_multi_machine_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss def record_desired_parameters(simulation): diff --git a/examples/models/ibb_manual/data/original_data.npy b/examples/models/ibb_manual/data/original_data.npy new file mode 100644 index 0000000..3bfee4e Binary files /dev/null and b/examples/models/ibb_manual/data/original_data.npy differ diff --git a/examples/models/ibb_manual/ibb_manual_sim.py b/examples/models/ibb_manual/ibb_manual_sim.py new file mode 100644 index 0000000..aef016d --- /dev/null +++ b/examples/models/ibb_manual/ibb_manual_sim.py @@ -0,0 +1,83 @@ +""" +This example shows how to manually create and simulate the IBB model. +""" +import os + +import numpy as np +import matplotlib.pyplot as plt +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss + +from diffpssi.power_sim_lib.models.synchronous_machine import SynchMachine +from diffpssi.power_sim_lib.models.static_models import * + + +def record_desired_parameters(simulation): + """ + Records the desired parameters of the simulation. + Args: + simulation: The simulation to record the parameters from. + + Returns: A list of the recorded parameters. + + """ + # Record the desired parameters + record_list = [ + simulation.busses[1].models[0].omega.real, + simulation.busses[1].models[0].e_q_st.real, + simulation.busses[1].models[0].e_d_st.real, + ] + return record_list + + +def main(): + """ + This function simulates the IBB model. + """ + parallel_sims = 1 + + sim = Pss(parallel_sims=parallel_sims, + sim_time=10, + time_step=0.005, + solver='heun', + ) + + sim.fn = 60 + sim.base_mva = 2200 + sim.base_voltage = 24 + + sim.add_bus(Bus(name='Bus 0', v_n=24)) + sim.add_bus(Bus(name='Bus 1', v_n=24)) + + sim.add_line(Line(name='L1', from_bus='Bus 0', to_bus='Bus 1', length=1, s_n=2200, v_n=24, unit='p.u.', + r=0, x=0.65, b=0, s_n_sys=2200, v_n_sys=24)) + + sim.add_generator(SynchMachine(name='IBB', bus='Bus 0', s_n=22000, v_n=24, p=-1998, v=0.995, h=3.5e7, d=0, + x_d=1.81, x_q=1.76, x_d_t=0.3, x_q_t=0.65, x_d_st=0.23, x_q_st=0.23, t_d0_t=8.0, + t_q0_t=1, t_d0_st=0.03, t_q0_st=0.07, f_n_sys=60, s_n_sys=2200, v_n_sys=24)) + sim.add_generator(SynchMachine(name='Gen 1', bus='Bus 1', s_n=2200, v_n=24, p=1998, v=1, h=3.5, d=0, x_d=1.81, + x_q=1.76, x_d_t=0.3, x_q_t=0.65, x_d_st=0.23, x_q_st=0.23, t_d0_t=8.0, t_q0_t=1, + t_d0_st=0.03, t_q0_st=0.07, f_n_sys=60, s_n_sys=2200, v_n_sys=24)) + + sim.set_slack_bus('Bus 0') + + sim.add_sc_event(1, 1.05, 'Bus 1') + sim.set_record_function(record_desired_parameters) + + # Run the simulation. Recorder format shall be [batch, timestep, value] + t, recorder = sim.run() + + # Plot the results + plt.figure() + for i in range(len(recorder[0, 0, :])): + plt.subplot(len(recorder[0, 0, :]), 1, i + 1) + plt.plot(t, recorder[0, :, i].real) + plt.ylabel('Parameter {}'.format(i)) + plt.xlabel('Time [s]') + plt.show() + + if os.environ.get('DIFFPSSI_TESTING') == 'True': + np.save('./data/original_data.npy', recorder[0].real) + + +if __name__ == '__main__': + main() diff --git a/examples/models/ibb_model/ibb_opt.py b/examples/models/ibb_model/ibb_opt.py index d4a6303..07bbb66 100644 --- a/examples/models/ibb_model/ibb_opt.py +++ b/examples/models/ibb_model/ibb_opt.py @@ -4,8 +4,8 @@ import numpy as np import torch import examples.models.ibb_model.ibb_model as mdl -from src.diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss np.random.seed(0) @@ -95,3 +95,7 @@ def main(parallel_sims=100): ) opt.run() + + +if __name__ == '__main__': + main() diff --git a/examples/models/ibb_model/ibb_sim.py b/examples/models/ibb_model/ibb_sim.py index 4d02c4b..6fea67a 100644 --- a/examples/models/ibb_model/ibb_sim.py +++ b/examples/models/ibb_model/ibb_sim.py @@ -4,7 +4,7 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.ibb_model.ibb_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss def record_desired_parameters(simulation): diff --git a/examples/models/ibb_transformer/ibb_trans_sim.py b/examples/models/ibb_transformer/ibb_trans_sim.py index ed69327..eeb4b5f 100644 --- a/examples/models/ibb_transformer/ibb_trans_sim.py +++ b/examples/models/ibb_transformer/ibb_trans_sim.py @@ -4,7 +4,7 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.ibb_transformer.ibb_trans_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss def record_desired_parameters(simulation): diff --git a/examples/models/ibb_with_controllers/ibb_wc_sim.py b/examples/models/ibb_with_controllers/ibb_wc_sim.py index 7a57879..53100bf 100644 --- a/examples/models/ibb_with_controllers/ibb_wc_sim.py +++ b/examples/models/ibb_with_controllers/ibb_wc_sim.py @@ -4,8 +4,8 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.ibb_with_controllers.ibb_wc_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss -from src.diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.backend import * def record_desired_parameters(simulation): diff --git a/examples/models/ieee_9bus/ieee_9bus_opt.py b/examples/models/ieee_9bus/ieee_9bus_opt.py index 35d5994..c99ce61 100644 --- a/examples/models/ieee_9bus/ieee_9bus_opt.py +++ b/examples/models/ieee_9bus/ieee_9bus_opt.py @@ -6,8 +6,8 @@ import torch import examples.models.ieee_9bus.ieee_9bus_model as mdl -from src.diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation +from diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation np.random.seed(0) diff --git a/examples/models/ieee_9bus/ieee_9bus_sim.py b/examples/models/ieee_9bus/ieee_9bus_sim.py index c0b09b9..143aa62 100644 --- a/examples/models/ieee_9bus/ieee_9bus_sim.py +++ b/examples/models/ieee_9bus/ieee_9bus_sim.py @@ -4,7 +4,7 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.ieee_9bus.ieee_9bus_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss def record_desired_parameters(simulation): diff --git a/examples/models/k2a/k2a_opt.py b/examples/models/k2a/k2a_opt.py index 878122c..4a8f6ed 100644 --- a/examples/models/k2a/k2a_opt.py +++ b/examples/models/k2a/k2a_opt.py @@ -6,8 +6,8 @@ import torch import examples.models.k2a.k2a_model as mdl -from src.diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation +from diffpssi.optimization_lib.ps_optimization import PowerSystemOptimization +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation np.random.seed(0) diff --git a/examples/models/k2a/k2a_sim.py b/examples/models/k2a/k2a_sim.py index 8e5e9d8..b97cd07 100644 --- a/examples/models/k2a/k2a_sim.py +++ b/examples/models/k2a/k2a_sim.py @@ -4,7 +4,7 @@ import numpy as np import matplotlib.pyplot as plt import examples.models.k2a.k2a_model as mdl -from src.diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss def record_desired_parameters(simulation): diff --git a/examples/models/k2a_manual/data/original_data.npy b/examples/models/k2a_manual/data/original_data.npy new file mode 100644 index 0000000..b2f1baf Binary files /dev/null and b/examples/models/k2a_manual/data/original_data.npy differ diff --git a/examples/models/k2a_manual/data/original_data_t.npy b/examples/models/k2a_manual/data/original_data_t.npy new file mode 100644 index 0000000..f732a23 Binary files /dev/null and b/examples/models/k2a_manual/data/original_data_t.npy differ diff --git a/examples/models/k2a_manual/k2a_manual_sim.py b/examples/models/k2a_manual/k2a_manual_sim.py new file mode 100644 index 0000000..fb167e3 --- /dev/null +++ b/examples/models/k2a_manual/k2a_manual_sim.py @@ -0,0 +1,167 @@ +""" +File contains an example of how to simulate the K2A model. +""" +import numpy as np +import matplotlib.pyplot as plt +from diffpssi.power_sim_lib.simulator import PowerSystemSimulation as Pss + +from diffpssi.power_sim_lib.models.synchronous_machine import SynchMachine +from diffpssi.power_sim_lib.models.static_models import * +from diffpssi.power_sim_lib.models.exciters import SEXS +from diffpssi.power_sim_lib.models.governors import TGOV1 +from diffpssi.power_sim_lib.models.stabilizers import STAB1 + + +def record_desired_parameters(simulation): + """ + Records the desired parameters of the simulation. + Args: + simulation: The simulation to record the parameters from. + + Returns: A list of the recorded parameters. + + """ + # Record the desired parameters + record_list = [ + simulation.busses[0].models[0].delta.real, + simulation.busses[1].models[0].delta.real, + simulation.busses[2].models[0].delta.real, + simulation.busses[3].models[0].delta.real, + + simulation.busses[0].models[0].omega.real, + simulation.busses[1].models[0].omega.real, + simulation.busses[2].models[0].omega.real, + simulation.busses[3].models[0].omega.real, + + simulation.busses[0].models[0].e_q_st.real, + simulation.busses[1].models[0].e_q_st.real, + simulation.busses[2].models[0].e_q_st.real, + simulation.busses[3].models[0].e_q_st.real, + + simulation.busses[0].models[0].e_d_st.real, + simulation.busses[1].models[0].e_d_st.real, + simulation.busses[2].models[0].e_d_st.real, + simulation.busses[3].models[0].e_d_st.real, + ] + return record_list + + +def main(): + """ + This function simulates the K2A model. + """ + parallel_sims = 1 + + sim = Pss(parallel_sims=parallel_sims, + sim_time=10, + time_step=0.005, + solver='heun', + ) + + sim.fn = 50 + sim.base_mva = 900 + sim.base_voltage = 230 + + sim.add_bus(Bus(name='B1', v_n=20)) + sim.add_bus(Bus(name='B2', v_n=20)) + sim.add_bus(Bus(name='B3', v_n=20)) + sim.add_bus(Bus(name='B4', v_n=20)) + sim.add_bus(Bus(name='B5', v_n=230)) + sim.add_bus(Bus(name='B6', v_n=230)) + sim.add_bus(Bus(name='B7', v_n=230)) + sim.add_bus(Bus(name='B8', v_n=230)) + sim.add_bus(Bus(name='B9', v_n=230)) + sim.add_bus(Bus(name='B10', v_n=230)) + sim.add_bus(Bus(name='B11', v_n=230)) + + sim.add_line(Line(name='L5-6', from_bus='B5', to_bus='B6', length=25, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L6-7', from_bus='B6', to_bus='B7', length=10, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L7-8-1', from_bus='B7', to_bus='B8', length=110, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L7-8-2', from_bus='B7', to_bus='B8', length=110, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L8-9-1', from_bus='B8', to_bus='B9', length=110, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L8-9-2', from_bus='B8', to_bus='B9', length=110, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L9-10', from_bus='B9', to_bus='B10', length=10, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + sim.add_line(Line(name='L10-11', from_bus='B10', to_bus='B11', length=25, s_n=100, v_n=230, unit='p.u.', + r=1e-4, x=1e-3, b=1.75e-3, s_n_sys=900, v_n_sys=230)) + + sim.add_transformer(Transformer(name='T1', from_bus='B1', to_bus='B5', s_n=900, v_n_from=20, v_n_to=230, r=0, + x=0.15, s_n_sys=900)) + sim.add_transformer(Transformer(name='T2', from_bus='B2', to_bus='B6', s_n=900, v_n_from=20, v_n_to=230, r=0, + x=0.15, s_n_sys=900)) + sim.add_transformer(Transformer(name='T3', from_bus='B3', to_bus='B11', s_n=900, v_n_from=20, v_n_to=230, r=0, + x=0.15, s_n_sys=900)) + sim.add_transformer(Transformer(name='T4', from_bus='B4', to_bus='B10', s_n=900, v_n_from=20, v_n_to=230, + r=0, x=0.15, s_n_sys=900)) + + sim.add_load(Load(name='L1', bus='B7', p=967, q=100, model='Z', s_n_sys=900)) + sim.add_load(Load(name='L2', bus='B9', p=1767, q=100, model='Z', s_n_sys=900)) + + sim.add_shunt(Shunt(name='C1', bus='B7', v_n=230, q=200, model='Z', s_n_sys=900)) + sim.add_shunt(Shunt(name='C2', bus='B9', v_n=230, q=350, model='Z', s_n_sys=900)) + + sim.add_generator(SynchMachine(name='G1', bus='B1', s_n=900, v_n=20, p=700, v=1.03, h=6.5, d=0, x_d=1.8, + x_q=1.7, x_d_t=0.3, x_q_t=0.55, x_d_st=0.25, x_q_st=0.25, t_d0_t=8.0, t_q0_t=0.4, + t_d0_st=0.03, t_q0_st=0.05, f_n_sys=50, s_n_sys=900, v_n_sys=230)) + sim.add_generator(SynchMachine(name='G2', bus='B2', s_n=900, v_n=20, p=700, v=1.01, h=6.5, d=0, x_d=1.8, + x_q=1.7, x_d_t=0.3, x_q_t=0.55, x_d_st=0.25, x_q_st=0.25, t_d0_t=8.0, t_q0_t=0.4, + t_d0_st=0.03, t_q0_st=0.05, f_n_sys=50, s_n_sys=900, v_n_sys=230)) + sim.add_generator(SynchMachine(name='G3', bus='B3', s_n=900, v_n=20, p=719, v=1.03, h=6.175, d=0, x_d=1.8, + x_q=1.7, x_d_t=0.3, x_q_t=0.55, x_d_st=0.25, x_q_st=0.25, t_d0_t=8.0, t_q0_t=0.4, + t_d0_st=0.03, t_q0_st=0.05, f_n_sys=50, s_n_sys=900, v_n_sys=230)) + sim.add_generator(SynchMachine(name='G4', bus='B4', s_n=900, v_n=20, p=700, v=1.01, h=6.175, d=0, x_d=1.8, + x_q=1.7, x_d_t=0.3, x_q_t=0.55, x_d_st=0.25, x_q_st=0.25, t_d0_t=8.0, t_q0_t=0.4, + t_d0_st=0.03, t_q0_st=0.05, f_n_sys=50, s_n_sys=900, v_n_sys=230)) + + sim.add_governor(TGOV1(name='GOV1', gen='G1', r=0.05, d_t=0.02, v_min=0, v_max=1, t_1=0.5, t_2=1, t_3=2)) + sim.add_governor(TGOV1(name='GOV2', gen='G2', r=0.05, d_t=0.02, v_min=0, v_max=1, t_1=0.5, t_2=1, t_3=2)) + sim.add_governor(TGOV1(name='GOV3', gen='G3', r=0.05, d_t=0.02, v_min=0, v_max=1, t_1=0.5, t_2=1, t_3=2)) + sim.add_governor(TGOV1(name='GOV4', gen='G4', r=0.05, d_t=0.02, v_min=0, v_max=1, t_1=0.5, t_2=1, t_3=2)) + + sim.add_exciter(SEXS(name='AVR1', gen='G1', k=100, t_a=2.0, t_b=10.0, t_e=0.1, e_min=-3, e_max=3)) + sim.add_exciter(SEXS(name='AVR2', gen='G2', k=100, t_a=2.0, t_b=10.0, t_e=0.1, e_min=-3, e_max=3)) + sim.add_exciter(SEXS(name='AVR3', gen='G3', k=100, t_a=2.0, t_b=10.0, t_e=0.1, e_min=-3, e_max=3)) + sim.add_exciter(SEXS(name='AVR4', gen='G4', k=100, t_a=2.0, t_b=10.0, t_e=0.1, e_min=-3, e_max=3)) + + sim.add_pss(STAB1(name='PSS1', gen='G1', k=50, t=10.0, t_1=0.5, t_2=0.5, t_3=0.05, t_4=0.05, h_lim=0.03)) + sim.add_pss(STAB1(name='PSS2', gen='G2', k=50, t=10.0, t_1=0.5, t_2=0.5, t_3=0.05, t_4=0.05, h_lim=0.03)) + sim.add_pss(STAB1(name='PSS3', gen='G3', k=50, t=10.0, t_1=0.5, t_2=0.5, t_3=0.05, t_4=0.05, h_lim=0.03)) + sim.add_pss(STAB1(name='PSS4', gen='G4', k=50, t=10.0, t_1=0.5, t_2=0.5, t_3=0.05, t_4=0.05, h_lim=0.03)) + + sim.set_slack_bus('B3') + + sim.add_sc_event(1, 1.1, 'B1') + + sim.set_record_function(record_desired_parameters) + t, recorder = sim.run() + + # Format shall be [batch, timestep, value] + # create a new subplot for each parameter + plt.figure() + + for i in range(len(recorder[0, 0, :])): + plt.subplot(len(recorder[0, 0, :]), 1, i + 1) + plt.plot(t, recorder[0, :, i].real) + plt.ylabel('Parameter {}'.format(i)) + plt.xlabel('Time [s]') + + plt.savefig('data/plots/original.png'.format()) + plt.show() + + # add time to the first column + saver = np.zeros((len(t), len(recorder[0, 0, :]) + 1)) + saver[:, 0] = t + saver[:, 1:] = recorder[0, :, :].real + + np.save('data/original_data_t.npy', saver) + np.save('data/original_data.npy', recorder[0].real) + + +if __name__ == '__main__': + main() diff --git a/pyproject.toml b/pyproject.toml index 9b25ab9..43ee8e4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,10 +4,17 @@ build-backend = "setuptools.build_meta" [project] name = "diffpssi" -version = "0.0.1" +version = "0.0.3" authors = [ { name="Georg Kordowich", email="georg.kordowich@fau.de" }, ] +dependencies = [ + "matplotlib", + "numpy", + "torch", + "tqdm", +] + description = "A differentiable power system simulation tool" readme = "README.md" requires-python = ">=3.8" diff --git a/readme.md b/readme.md index acd8032..95dc4a5 100644 --- a/readme.md +++ b/readme.md @@ -1,4 +1,18 @@ # DiffPSSi: A framework for differentiable power system simulations +## Quickstart +1. Install Package: + ``` + pip install diffpssi + ``` +2. Copy-Paste example simulation from: + ``` + examples/models/ibb_manual/ibb_sim_manual.py + ``` +3. Run the simulation: + ``` + python ibb_sim_manual.py + ``` + ## Overview DiffPSSi contains a framework designed for simulating and optimizing the dynamic behavior of power systems. The core idea of the framework is to allow the use of automatic differentiation for dynamic @@ -14,13 +28,15 @@ The code is strongly based on [this](https://github.com/hallvar-h/DynPSSimPy) re enable the gradient calculation for optimization purposes. The code is still under development and will be extended in the future. +**Note: This repository is still under development and will be extended in the future. Use at your own risk.** + ## Features -- **Inherently Parallel Implementation:** A unique and most important feature of this simulation framework, as it allows the execution of multiple simulations in parallel by using vectors of parameters for every element. +- **Inherently Parallel Implementation:** A unique and important feature of this simulation framework, as it allows the execution of multiple simulations in parallel by using vectors of parameters for every element. - **Dynamic Simulation:** Allows for detailed dynamic simulations of power systems, including interactions between various components. - **Optimization Library:** Features optimizers based on BFGS and automatic differentiation for efficient gradient computation. This part can be used for parameter optimization or identification purposes. - **Extensible Model Library:** Contains models of AVRs, governors, stabilizers, static models like lines, loads, transformers, and more. - **Backend Flexibility:** Choose between `torch` and `numpy` as backend for computations. -- **Solver Options:** Includes Euler and Heun methods for numerical integration. +- **Solver Options:** Includes Euler and Runge Kutta methods for numerical integration. ## Installation @@ -28,9 +44,9 @@ in the future. ``` git clone git@github.com:georgkordowich/diffpssi.git ``` -2. Install the required dependencies: +2. Install the package: ``` - pip install -r requirements.txt + pip install diffpssi ``` 3. Run example simulations: ``` @@ -38,10 +54,43 @@ in the future. ``` ## Usage and Examples -Detailed usage instructions and examples can be found in the `examples` directory. For a guided introduction, check out the IBB model simulation example under `examples/models/ibb_model/ibb_sim.py` and the corresponding model in `examples/models/ibb_model.py`. +Detailed usage instructions and examples can be found in the `examples` directory. -To define a model, you first need to creat the power system simulation (Pss) itself. Afterward, you can add busses, -generators, and lines as desired. +Generally, there are two options to create a simulation. One option is to create the model as a dictionary and pass it +to the simulation. This is the recommended way. For an example, check out +the IBB model simulation example under `examples/models/ibb_model/ibb_sim.py` and the corresponding model +in `examples/models/ibb_model.py`. + +The other option is to create the model manually in the simulation file. This can be seen in the example under +`examples/models/ibb_model/ibb_sim_manual.py`. For this option, first a simulation must be created and afterward, +busses, generators, and lines can be added to the simulation. The following code snippet shows how to create +a simulation and add busses, generators, and lines to it. +```python +sim = Pss(parallel_sims=parallel_sims, + sim_time=10, + time_step=0.005, + solver='heun', + ) + +sim.fn = 60 +sim.base_mva = 2200 +sim.base_voltage = 24 + +sim.add_bus(Bus(name='Bus 0', v_n=24)) +sim.add_bus(Bus(name='Bus 1', v_n=24)) + +sim.add_line(Line(name='L1', from_bus='Bus 0', to_bus='Bus 1', length=1, s_n=2200, v_n=24, unit='p.u.', + r=0, x=0.65, b=0, s_n_sys=2200, v_n_sys=24)) + +sim.add_generator(SynchMachine(name='IBB', bus='Bus 0', s_n=22000, v_n=24, p=-1998, v=0.995, h=3.5e7, d=0, + x_d=1.81, x_q=1.76, x_d_t=0.3, x_q_t=0.65, x_d_st=0.23, x_q_st=0.23, t_d0_t=8.0, + t_q0_t=1, t_d0_st=0.03, t_q0_st=0.07, f_n_sys=60, s_n_sys=2200, v_n_sys=24)) +sim.add_generator(SynchMachine(name='Gen 1', bus='Bus 1', s_n=2200, v_n=24, p=1998, v=1, h=3.5, d=0, x_d=1.81, + x_q=1.76, x_d_t=0.3, x_q_t=0.65, x_d_st=0.23, x_q_st=0.23, t_d0_t=8.0, t_q0_t=1, + t_d0_st=0.03, t_q0_st=0.07, f_n_sys=60, s_n_sys=2200, v_n_sys=24)) + +sim.set_slack_bus('Bus 0') +``` Once the model is defined, you can run a simulation. One unique feature of this framework is that you can define the number of parallel simulations to run. This is useful for parameter optimization, where you can run multiple diff --git a/src/diffpssi/optimization_lib/ps_optimization.py b/src/diffpssi/optimization_lib/ps_optimization.py index ce19fad..cfc74d1 100644 --- a/src/diffpssi/optimization_lib/ps_optimization.py +++ b/src/diffpssi/optimization_lib/ps_optimization.py @@ -6,7 +6,7 @@ import torch from matplotlib import pyplot as plt -from src.diffpssi.optimization_lib.optimizers import CustomBFGSREALOptimizer +from diffpssi.optimization_lib.optimizers import CustomBFGSREALOptimizer # currently only bfgs is supported, as it works best by far optimizer_dict = { @@ -52,7 +52,7 @@ def __init__(self, self.sim = sim if sim.backend == 'numpy': - raise NotImplementedError('Optimization is only supported for the PyTorch backend.' + raise NotImplementedError('Optimization is only supported for the PyTorch backend. ' 'Please set the backend to PyTorch in power_sim_lib/backend.py') self.target_data = original_data @@ -84,6 +84,7 @@ def default_loss_function(sim_result, target_data): Returns: A vector of the mean absolute error for each batch element of the size (batch-size) """ + # noinspection PyArgumentList return torch.mean(torch.sum(torch.abs(target_data - sim_result), dim=2), axis=1) self.loss_function = default_loss_function @@ -151,11 +152,17 @@ def run(self, max_steps=100): Args: max_steps: The maximum number of optimization steps that should be performed. """ - if os.environ['DIFFPSSI_FORCE_OPT_ITERS'] is not None: - max_steps = int(os.environ['DIFFPSSI_FORCE_OPT_ITERS']) + if os.environ.get('DIFFPSSI_FORCE_OPT_ITERS') is not None: + max_steps = int(os.environ.get('DIFFPSSI_FORCE_OPT_ITERS')) print('WARNING: FORCING THE USE OF {} OPTIMIZATION ITERATION.' - 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ['DIFFPSSI_FORCE_OPT_ITERS'])) + 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ.get('DIFFPSSI_FORCE_OPT_ITERS'))) opt_start_time = time.time() + + min_loss_idx = None # the index of the current best batch element + results = None # the simulation results + t = None # the timesteps + opt_step = None # the current optimization step + for opt_step in range(max_steps): opt_step_start = time.time() # set the gradients to zero in order to accumulate the diff --git a/src/diffpssi/power_sim_lib/backend.py b/src/diffpssi/power_sim_lib/backend.py index b0d1f9a..bed4972 100644 --- a/src/diffpssi/power_sim_lib/backend.py +++ b/src/diffpssi/power_sim_lib/backend.py @@ -8,12 +8,10 @@ BACKEND = 'numpy' # 'torch' or 'numpy' DEFAULT_DEVICE = 'cpu' # 'cuda' or 'cpu' -try: - BACKEND = os.environ['DIFFPSSI_FORCE_SIM_BACKEND'] +if os.environ.get('DIFFPSSI_FORCE_SIM_BACKEND') is not None: + BACKEND = os.environ.get('DIFFPSSI_FORCE_SIM_BACKEND') print('WARNING: FORCING THE USE OF {} AS BACKEND.' - 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ['DIFFPSSI_FORCE_SIM_BACKEND'])) -except KeyError: - pass + 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ.get('DIFFPSSI_FORCE_SIM_BACKEND'))) if BACKEND == 'torch': import torch diff --git a/src/diffpssi/power_sim_lib/load_flow.py b/src/diffpssi/power_sim_lib/load_flow.py index dcdab3a..8166c4c 100644 --- a/src/diffpssi/power_sim_lib/load_flow.py +++ b/src/diffpssi/power_sim_lib/load_flow.py @@ -3,7 +3,7 @@ """ from itertools import count import time -from src.diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.backend import * def do_load_flow(ps_sim): diff --git a/src/diffpssi/power_sim_lib/models/blocks.py b/src/diffpssi/power_sim_lib/models/blocks.py index 7c9b097..2f4260e 100644 --- a/src/diffpssi/power_sim_lib/models/blocks.py +++ b/src/diffpssi/power_sim_lib/models/blocks.py @@ -1,7 +1,7 @@ """ File contains implemented controller blocks for power system simulations. """ -from src.diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.backend import * class PT1Limited(object): @@ -20,7 +20,7 @@ class PT1Limited(object): state_1 (float or torch.Tensor): Internal state of the PT1Limited block. """ - def __init__(self, t_pt1, gain_pt1, lim_min, lim_max, parallel_sims=None): + def __init__(self, t_pt1, gain_pt1, lim_min, lim_max): """ Initializes the PT1Limited block with specified parameters. @@ -29,13 +29,11 @@ def __init__(self, t_pt1, gain_pt1, lim_min, lim_max, parallel_sims=None): gain_pt1 (float): Gain of the PT1 transfer function. lim_min (float): Minimum limit for the output. lim_max (float): Maximum limit for the output. - parallel_sims (int, optional): Number of parallel simulations to enable. """ self.t_pt1 = t_pt1 self.gain_pt1 = gain_pt1 self.lim_min = lim_min self.lim_max = lim_max - self.enable_parallel_simulation(parallel_sims) self.input = 0 @@ -124,7 +122,7 @@ class Limiter(object): limit (float or torch.Tensor): The limit value for both positive and negative sides. """ - def __init__(self, limit, parallel_sims=None): + def __init__(self, limit): """ Initializes the Limiter block with a specified limit. @@ -133,7 +131,6 @@ def __init__(self, limit, parallel_sims=None): parallel_sims (int, optional): Number of parallel simulations to enable. """ self.limit = limit - self.enable_parallel_simulation(parallel_sims) def get_output(self, input_var): """ @@ -173,7 +170,7 @@ class LeadLag(object): state_1 (float or torch.Tensor): Internal state of the LeadLag block. """ - def __init__(self, t_1, t_2, parallel_sims=None): + def __init__(self, t_1, t_2): """ Initializes the LeadLag block with specified parameters. @@ -184,7 +181,6 @@ def __init__(self, t_1, t_2, parallel_sims=None): """ self.t_1 = t_1 self.t_2 = t_2 - self.enable_parallel_simulation(parallel_sims) self.input = 0 self.state_1 = 0 @@ -265,7 +261,7 @@ class Washout(object): state_1 (float or torch.Tensor): Internal state of the Washout block. """ - def __init__(self, k_w, t_w, parallel_sims=None): + def __init__(self, k_w, t_w): """ Initializes the Washout filter with specified parameters. @@ -276,7 +272,6 @@ def __init__(self, k_w, t_w, parallel_sims=None): """ self.k_w = k_w self.t_w = t_w - self.enable_parallel_simulation(parallel_sims) self.input = 0 self.state_1 = 0 diff --git a/src/diffpssi/power_sim_lib/models/exciters.py b/src/diffpssi/power_sim_lib/models/exciters.py index d65c66a..8e93990 100644 --- a/src/diffpssi/power_sim_lib/models/exciters.py +++ b/src/diffpssi/power_sim_lib/models/exciters.py @@ -2,8 +2,8 @@ Contains the SEXS (Simplified Excitation System) model for power system simulations. Other models can be added here as well. """ -from src.diffpssi.power_sim_lib.backend import * -from src.diffpssi.power_sim_lib.models.blocks import LeadLag, PT1Limited +from diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.models.blocks import LeadLag, PT1Limited class SEXS(object): @@ -20,16 +20,44 @@ class SEXS(object): v_setpoint (float or torch.Tensor): The setpoint for the system voltage. bias (float or torch.Tensor): A bias value for voltage control. """ - def __init__(self, param_dict, parallel_sims, v_setpoint): + + def __init__(self, param_dict=None, + name=None, + gen=None, + t_a=None, + t_b=None, + t_e=None, + k=None, + e_min=None, + e_max=None, + ): """ Initializes the SEXS model with specified parameters. Args: param_dict (dict, optional): A dictionary of parameters for the model. - parallel_sims (int, optional): Number of parallel simulations to enable. - v_setpoint (float, optional): The setpoint for the system voltage. Default is 1.0. + name (str, optional): The name of the model. + gen (str, optional): The name of the generator the model is connected to. + t_a (float, optional): The time constant t1 of the lead-lag control block. + t_b (float, optional): The time constant t2 of the lead-lag control block. + t_e (float, optional): The time constant t of the PT1 transfer function block. + k (float, optional): The gain of the PT1 transfer function block. + e_min (float, optional): The minimum value of the output of the PT1 transfer function block. + e_max (float, optional): The maximum value of the output of the PT1 transfer function block. """ + if param_dict is None: + param_dict = { + 'name': name, + 'gen': gen, + 'T_a': t_a, + 'T_b': t_b, + 'T_e': t_e, + 'K': k, + 'E_min': e_min, + 'E_max': e_max, + } self.name = param_dict['name'] + self.gen = param_dict['gen'] self.t_a = param_dict['T_a'] self.t_b = param_dict['T_b'] self.t_e = param_dict['T_e'] @@ -37,11 +65,10 @@ def __init__(self, param_dict, parallel_sims, v_setpoint): self.e_min = param_dict['E_min'] self.e_max = param_dict['E_max'] - self.lead_lag = LeadLag(t_1=self.t_a, t_2=self.t_b, parallel_sims=parallel_sims) - self.pt1_lim = PT1Limited(t_pt1=self.t_e, gain_pt1=self.k, lim_min=self.e_min, lim_max=self.e_max, - parallel_sims=parallel_sims) + self.lead_lag = LeadLag(t_1=self.t_a, t_2=self.t_b) + self.pt1_lim = PT1Limited(t_pt1=self.t_e, gain_pt1=self.k, lim_min=self.e_min, lim_max=self.e_max) - self.v_setpoint = v_setpoint + self.v_setpoint = 0.0 self.bias = 0.0 def differential(self): @@ -96,6 +123,8 @@ def enable_parallel_simulation(self, parallel_sims): """ self.v_setpoint = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.v_setpoint self.bias = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.bias + self.lead_lag.enable_parallel_simulation(parallel_sims) + self.pt1_lim.enable_parallel_simulation(parallel_sims) def initialize(self, e_fd): """ diff --git a/src/diffpssi/power_sim_lib/models/governors.py b/src/diffpssi/power_sim_lib/models/governors.py index c95d9a6..f327633 100644 --- a/src/diffpssi/power_sim_lib/models/governors.py +++ b/src/diffpssi/power_sim_lib/models/governors.py @@ -1,8 +1,8 @@ """ File containing the TGOV1 model. Other governor models can be added here. """ -from src.diffpssi.power_sim_lib.backend import * -from src.diffpssi.power_sim_lib.models.blocks import LeadLag, PT1Limited +from diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.models.blocks import LeadLag, PT1Limited class TGOV1(object): @@ -18,16 +18,48 @@ class TGOV1(object): lead_lag (LeadLag): The LeadLag block representing the compensator. p_ref (float or torch.Tensor): The reference power setpoint for the governor. """ - def __init__(self, param_dict, parallel_sims): + + def __init__(self, param_dict=None, + name=None, + gen=None, + r=None, + d_t=None, + t_1=None, + t_2=None, + t_3=None, + v_min=None, + v_max=None, + ): """ Initializes the TGOV1 model with specified parameters. Args: param_dict (dict, optional): A dictionary of parameters for the model. - parallel_sims (int, optional): Number of parallel simulations to enable. + name (str, optional): The name of the model. + gen (str, optional): The name of the generator the model is connected to. + r (float, optional): The droop of the governor. + d_t (float, optional): The damping coefficient of the governor. + t_1 (float, optional): The time constant t1 of the PT1Limited block. + t_2 (float, optional): The time constant t1 of the first lead-lag compensator. + t_3 (float, optional): The time constant t2 of the first lead-lag compensator. + v_min (float, optional): The minimum value of the governor output. + v_max (float, optional): The maximum value of the governor output. """ + if param_dict is None: + param_dict = { + 'name': name, + 'gen': gen, + 'R': r, + 'D_t': d_t, + 'T_1': t_1, + 'T_2': t_2, + 'T_3': t_3, + 'V_min': v_min, + 'V_max': v_max, + } self.name = param_dict['name'] + self.gen = param_dict['gen'] self.r = param_dict['R'] self.d_t = param_dict['D_t'] self.t_1 = param_dict['T_1'] @@ -37,9 +69,8 @@ def __init__(self, param_dict, parallel_sims): self.v_max = param_dict['V_max'] droop = 1 / self.r - self.pt1_lim = PT1Limited(t_pt1=self.t_1, gain_pt1=droop, lim_min=self.v_min, lim_max=self.v_max, - parallel_sims=parallel_sims) - self.lead_lag = LeadLag(t_1=self.t_2, t_2=self.t_3, parallel_sims=parallel_sims) + self.pt1_lim = PT1Limited(t_pt1=self.t_1, gain_pt1=droop, lim_min=self.v_min, lim_max=self.v_max) + self.lead_lag = LeadLag(t_1=self.t_2, t_2=self.t_3) self.p_ref = 0.0 @@ -99,6 +130,8 @@ def enable_parallel_simulation(self, parallel_sims): parallel_sims (int): Number of parallel simulations. """ self.p_ref = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.p_ref + self.pt1_lim.enable_parallel_simulation(parallel_sims) + self.lead_lag.enable_parallel_simulation(parallel_sims) def initialize(self, p_mech): """ diff --git a/src/diffpssi/power_sim_lib/models/stabilizers.py b/src/diffpssi/power_sim_lib/models/stabilizers.py index c660cbf..773e113 100644 --- a/src/diffpssi/power_sim_lib/models/stabilizers.py +++ b/src/diffpssi/power_sim_lib/models/stabilizers.py @@ -1,8 +1,8 @@ """ Contains the STAB1 model for power system stabilizers. Other models can be added here as well. """ -from src.diffpssi.power_sim_lib.backend import * -from src.diffpssi.power_sim_lib.models.blocks import LeadLag, Washout, Limiter +from diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.models.blocks import LeadLag, Washout, Limiter class STAB1(object): @@ -20,15 +20,45 @@ class STAB1(object): limiter (Limiter): The Limiter block to restrict the output within a specific range. """ - def __init__(self, param_dict, parallel_sims): + def __init__(self, param_dict=None, + name=None, + gen=None, + k=None, + t=None, + t_1=None, + t_2=None, + t_3=None, + t_4=None, + h_lim=None, + ): """ Initializes the STAB1 model with specified parameters. Args: param_dict (dict, optional): A dictionary of parameters for the model. - parallel_sims (int, optional): Number of parallel simulations to enable. + name (str, optional): The name of the model. + gen (str, optional): The name of the generator the model is connected to. + k (float, optional): The gain of the washout filter. + t (float, optional): The time constant of the washout filter. + t_1 (float, optional): The time constant t1 of the first lead-lag compensator. + t_2 (float, optional): The time constant t1 of the second lead-lag compensator. + t_3 (float, optional): The time constant t2 of the first lead-lag compensator. + t_4 (float, optional): The time constant t2 of the second lead-lag compensator. """ + if param_dict is None: + param_dict = { + 'name': name, + 'gen': gen, + 'K': k, + 'T': t, + 'T_1': t_1, + 'T_2': t_2, + 'T_3': t_3, + 'T_4': t_4, + 'H_lim': h_lim, + } self.name = param_dict['name'] + self.gen = param_dict['gen'] self.k_w = param_dict['K'] self.t_w = param_dict['T'] self.t_1 = param_dict['T_1'] @@ -37,10 +67,10 @@ def __init__(self, param_dict, parallel_sims): self.t_4 = param_dict['T_4'] self.h_lim = param_dict['H_lim'] - self.washout = Washout(k_w=self.k_w, t_w=self.t_w, parallel_sims=parallel_sims) - self.lead_lag1 = LeadLag(t_1=self.t_1, t_2=self.t_3, parallel_sims=parallel_sims) - self.lead_lag2 = LeadLag(t_1=self.t_2, t_2=self.t_4, parallel_sims=parallel_sims) - self.limiter = Limiter(limit=self.h_lim, parallel_sims=parallel_sims) + self.washout = Washout(k_w=self.k_w, t_w=self.t_w) + self.lead_lag1 = LeadLag(t_1=self.t_1, t_2=self.t_3) + self.lead_lag2 = LeadLag(t_1=self.t_2, t_2=self.t_4) + self.limiter = Limiter(limit=self.h_lim) def differential(self): """ @@ -97,7 +127,10 @@ def enable_parallel_simulation(self, parallel_sims): Args: parallel_sims (int): Number of parallel simulations. """ - pass + self.washout.enable_parallel_simulation(parallel_sims) + self.lead_lag1.enable_parallel_simulation(parallel_sims) + self.lead_lag2.enable_parallel_simulation(parallel_sims) + self.limiter.enable_parallel_simulation(parallel_sims) def initialize(self, v_pss): """ diff --git a/src/diffpssi/power_sim_lib/models/static_models.py b/src/diffpssi/power_sim_lib/models/static_models.py index 035edcc..7e10032 100644 --- a/src/diffpssi/power_sim_lib/models/static_models.py +++ b/src/diffpssi/power_sim_lib/models/static_models.py @@ -5,7 +5,7 @@ """ import torch -from src.diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.backend import * class ScEvent(object): @@ -20,6 +20,7 @@ class ScEvent(object): end_time (float): The end time of the short circuit event. bus (int): The index of the bus where the short circuit occurs. """ + def __init__(self, start_time, end_time, bus): """ Initializes the ScEvent object with the start time, end time, and bus index. @@ -63,23 +64,31 @@ class Bus(object): models (list): List of models (generators, loads, etc.) connected to the bus. voltage (torch.Tensor): Current voltage at the bus. """ - def __init__(self, param_dict, parallel_sims): + + def __init__(self, param_dict=None, + name=None, + v_n=None): """ Initializes the Bus object with the given name, load flow type, and nominal voltage. Args: param_dict (dict): Dictionary of parameters for the bus. - parallel_sims (int): Number of parallel simulations to run. + name (str): The name of the bus. + v_n (float): Nominal voltage at the bus in kV. """ - self.parallel_sims = parallel_sims + self.parallel_sims = None + if param_dict is None: + param_dict = { + 'name': name, + 'V_n': v_n, + } + self.name = param_dict['name'] self.v_n = param_dict['V_n'] self.models = [] self.lf_type = 'PQ' - self.voltage = torch.ones((parallel_sims, 1), dtype=torch.complex128) - - self.enable_parallel_simulation(parallel_sims) + self.voltage = 1.0 def get_current_injections(self): """ @@ -144,6 +153,8 @@ def enable_parallel_simulation(self, parallel_sims): parallel_sims (int): Number of parallel simulations. """ self.v_n = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.v_n + self.voltage = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.voltage + self.parallel_sims = parallel_sims class Line(object): @@ -154,40 +165,69 @@ class Line(object): reactance, and susceptance, between two buses in the power system. Attributes: - from_bus (int): The index of the starting bus of the line. - to_bus (int): The index of the ending bus of the line. + from_bus_id (int): The index of the starting bus of the line. + to_bus_id (int): The index of the ending bus of the line. r (torch.Tensor): Resistance of the line. x (torch.Tensor): Reactance of the line. b (torch.Tensor): Susceptance of the line. """ - def __init__(self, param_dict, bus_from, bus_to, sys_mva, sys_voltage, parallel_sims): + def __init__(self, s_n_sys, v_n_sys, + param_dict=None, + name=None, + from_bus=None, + to_bus=None, + length=None, + s_n=None, + v_n=None, + unit=None, + r=None, + x=None, + b=None): """ Initializes the Line object with the specified electrical parameters and connected buses. Args: - param_dict (dict): Dictionary of parameters for the line. - bus_from (int): The index of the starting bus of the line. - bus_to (int): The index of the ending bus of the line. - sys_mva (float): Base power in MVA for the system. - sys_voltage (float): Base voltage in kV for the system. - parallel_sims (int): Number of parallel simulations to run. - """ + param_dict (dict, optional): Dictionary of parameters for the line. + name (str, optional): The name of the line. + from_bus (str, optional): The name of the starting bus of the line. + to_bus (str, optional): The name of the ending bus of the line. + length (float, optional): The length of the line in km. + s_n (float, optional): The nominal power of the line in MVA. + v_n (float, optional): The nominal voltage of the line in kV. + unit (str, optional): The unit of the line parameters. Either 'Ohm' or 'p.u.'. + r (float, optional): The resistance of the line (either in Ohm or p.u.)/length. + x (float, optional): The reactance of the line (either in Ohm or p.u.)/length. + b (float, optional): The susceptance of the line (either in Ohm or p.u.)/length. + """ + if param_dict is None: + param_dict = { + 'name': name, + 'from_bus': from_bus, + 'to_bus': to_bus, + 'length': length, + 'S_n': s_n, + 'V_n': v_n, + 'unit': unit, + 'R': r, + 'X': x, + 'B': b, + } self.name = param_dict['name'] - self.from_bus = bus_from - self.to_bus = bus_to + self.from_bus_name = param_dict['from_bus'] + self.to_bus_name = param_dict['to_bus'] + + self.from_bus_id = None + self.to_bus_id = None length = param_dict['length'] - s_n = param_dict.get('S_n', sys_mva) - v_n = param_dict.get('V_n', sys_voltage) + s_n = param_dict.get('S_n', s_n_sys) + v_n = param_dict.get('V_n', v_n_sys) unit = param_dict['unit'] r = param_dict['R'] x = param_dict['X'] b = param_dict['B'] - s_n_sys = sys_mva - v_n_sys = sys_voltage - z_n_sys = v_n_sys ** 2 / s_n_sys z_n = v_n ** 2 / s_n @@ -202,8 +242,6 @@ def __init__(self, param_dict, bus_from, bus_to, sys_mva, sys_voltage, parallel_ else: raise ValueError('Unit not supported') - self.enable_parallel_simulation(parallel_sims) - def get_admittance_diagonal(self): """ Calculates the diagonal admittance value of the line. @@ -229,14 +267,14 @@ def enable_parallel_simulation(self, parallel_sims): Args: parallel_sims (int): Number of parallel simulations. """ - self.from_bus = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.from_bus - self.to_bus = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.to_bus + self.from_bus_id = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.from_bus_id + self.to_bus_id = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.to_bus_id self.r = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.r self.x = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.x self.b = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.b -class Transfomer(object): +class Transformer(object): """ Represents a transformer in the power system simulation. @@ -244,36 +282,66 @@ class Transfomer(object): reactance, and connection between two buses in the power system. Attributes: - from_bus (int): The index of the primary side bus of the transformer. - to_bus (int): The index of the secondary side bus of the transformer. - r (torch.Tensor): Resistance of the transformer. - x (torch.Tensor): Reactance of the transformer. - b (torch.Tensor): Susceptance of the transformer (usually zero). + from_bus_id (int): The index of the primary side bus of the transformer. + to_bus_id (int): The index of the secondary side bus of the transformer. + r (torch.Tensor): Resistance of the transformer in p.u. + x (torch.Tensor): Reactance of the transformer in p.u. + b (torch.Tensor): Susceptance of the transformer in p.u. """ - def __init__(self, param_dict, bus_from, bus_to, s_n_sys, parallel_sims): + + def __init__(self, s_n_sys, + param_dict=None, + name=None, + from_bus=None, + to_bus=None, + s_n=None, + r=None, + x=None, + v_n_from=None, + v_n_to=None, + b=0): """ Initializes the Transformer object with the specified electrical parameters and connected buses. Args: - param_dict (dict): Dictionary of parameters for the transformer. - bus_from (int): The index of the primary side bus of the transformer. - bus_to (int): The index of the secondary side bus of the transformer. - s_n_sys (float): Base power in MVA for the system. - parallel_sims (int): Number of parallel simulations to run. - """ + param_dict (dict, optional): Dictionary of parameters for the transformer. + name (str, optional): The name of the transformer. + from_bus (str, optional): The name of the primary side bus of the transformer. + to_bus (str, optional): The name of the secondary side bus of the transformer. + s_n (float, optional): The nominal power of the transformer in MVA. + r (float, optional): The resistance of the transformer in p.u. + x (float, optional): The reactance of the transformer in p.u. + v_n_from (float, optional): The nominal voltage of the primary side of the transformer in kV. + v_n_to (float, optional): The nominal voltage of the secondary side of the transformer in kV. + b (float, optional): The susceptance of the transformer in p.u. + """ + if param_dict is None: + param_dict = { + 'name': name, + 'from_bus': from_bus, + 'to_bus': to_bus, + 'S_n': s_n, + 'R': r, + 'X': x, + 'V_n_from': v_n_from, + 'V_n_to': v_n_to, + 'B': b, + } self.name = param_dict['name'] - self.from_bus = bus_from - self.to_bus = bus_to + self.from_bus_name = param_dict['from_bus'] + self.to_bus_name = param_dict['to_bus'] + + self.from_bus_id = None + self.to_bus_id = None self.s_n = param_dict['S_n'] self.r = param_dict['R'] self.x = param_dict['X'] self.v_n_from = param_dict['V_n_from'] self.v_n_to = param_dict['V_n_to'] - self.s_n_sys = s_n_sys self.b = param_dict.get('B', 0) - self.enable_parallel_simulation(parallel_sims) + self.s_n_sys = s_n_sys def get_admittance_diagonal(self): """ @@ -282,7 +350,7 @@ def get_admittance_diagonal(self): Returns: torch.Tensor: Diagonal admittance of the transformer. """ - return (1 / (self.r + 1j * self.x) + 1j * self.b / 2)*self.s_n/self.s_n_sys + return (1 / (self.r + 1j * self.x) + 1j * self.b / 2) * self.s_n / self.s_n_sys def get_admittance_off_diagonal(self): """ @@ -291,7 +359,7 @@ def get_admittance_off_diagonal(self): Returns: torch.Tensor: Off-diagonal admittance of the transformer. """ - return -1 / (self.r + 1j * self.x)*self.s_n/self.s_n_sys + return -1 / (self.r + 1j * self.x) * self.s_n / self.s_n_sys def enable_parallel_simulation(self, parallel_sims): """ @@ -300,8 +368,8 @@ def enable_parallel_simulation(self, parallel_sims): Args: parallel_sims (int): Number of parallel simulations. """ - self.from_bus = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.from_bus - self.to_bus = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.to_bus + self.from_bus_id = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.from_bus_id + self.to_bus_id = torch.ones((parallel_sims, 1), dtype=torch.int32) * self.to_bus_id self.r = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.r self.x = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.x self.v_n_from = torch.ones((parallel_sims, 1), dtype=torch.float64) * self.v_n_from @@ -319,17 +387,35 @@ class Load(object): p_soll_mw (float): Desired active power (MW) of the load. q_soll_mvar (float): Desired reactive power (MVAR) of the load. """ - def __init__(self, param_dict, s_n_sys, v_n_sys, parallel_sims): + + def __init__(self, s_n_sys, + param_dict=None, + name=None, + bus=None, + p=None, + q=None, + model=None): """ Initializes the Load object with specified active and reactive power demands. Args: - param_dict (dict): Dictionary of parameters for the load. - s_n_sys (float): Base power in MVA for the system. - v_n_sys (float): Base voltage in kV for the system. - parallel_sims (int): Number of parallel simulations to run. - """ + param_dict (dict, optional): Dictionary of parameters for the load. + name (str, optional): The name of the load. + bus (str, optional): The name of the bus where the load is connected. + p (float, optional): The active power demand of the load in MW. + q (float, optional): The reactive power demand of the load in MVAR. + model (str, optional): The model of the load. Currently only 'Z' is supported. + """ + if param_dict is None: + param_dict = { + 'name': name, + 'bus': bus, + 'P': p, + 'Q': q, + 'model': model, + } self.name = param_dict['name'] + self.bus = param_dict['bus'] self.p_soll_mw = param_dict['P'] self.q_soll_mvar = param_dict['Q'] self.model = param_dict['model'] @@ -341,8 +427,6 @@ def __init__(self, param_dict, s_n_sys, v_n_sys, parallel_sims): self.y_load = None - self.enable_parallel_simulation(parallel_sims) - def get_lf_power(self): """ Calculates and returns the load flow power of the Load. @@ -354,7 +438,7 @@ def get_lf_power(self): complex: The calculated load flow power. """ - return -(self.p_soll_mw + 1j * self.q_soll_mvar)/self.s_n_sys + return -(self.p_soll_mw + 1j * self.q_soll_mvar) / self.s_n_sys def get_admittance(self, dyn): """ @@ -376,6 +460,7 @@ def get_admittance(self, dyn): else: return torch.zeros_like(self.p_soll_mw) + # noinspection PyUnusedLocal def initialize(self, s_calc, v_bb): """ Initializes the Load by calculating its admittance. @@ -386,7 +471,7 @@ def initialize(self, s_calc, v_bb): s_calc (float): The calculated complex power. v_bb (torch.Tensor): Busbar voltage value. """ - s_load = (self.p_soll_mw + 1j * self.q_soll_mvar)/self.s_n_sys + s_load = (self.p_soll_mw + 1j * self.q_soll_mvar) / self.s_n_sys z_load = torch.conj(torch.abs(v_bb) ** 2 / s_load) self.y_load = 1 / z_load @@ -440,23 +525,41 @@ class Shunt(object): Attributes: q_soll_mvar (float): Desired reactive power (MVAR) of the load. """ - def __init__(self, param_dict, s_n_sys, v_n_sys, parallel_sims): + + def __init__(self, s_n_sys, + param_dict=None, + name=None, + bus=None, + v_n=None, + q=None, + model=None): """ Initializes the Load object with specified active and reactive power demands. Args: - param_dict (dict): Dictionary of parameters for the load. - s_n_sys (float): Base power in MVA for the system. - v_n_sys (float): Base voltage in kV for the system. - parallel_sims (int): Number of parallel simulations to run. - """ + param_dict (dict, optional): Dictionary of parameters for the load. + name (str, optional): The name of the load. + bus (str, optional): The name of the bus where the load is connected. + v_n (float, optional): The nominal voltage of the load in kV. + q (float, optional): The reactive power demand of the load in MVAR. + model (str, optional): The model of the load. Currently only 'Z' is supported. + """ + if param_dict is None: + param_dict = { + 'name': name, + 'bus': bus, + 'V_n': v_n, + 'Q': q, + 'model': model, + } self.name = param_dict['name'] + self.bus = param_dict['bus'] self.v_n = param_dict['V_n'] self.q_soll_mvar = param_dict['Q'] self.model = param_dict['model'] if not self.model == 'Z': - raise ValueError('Only Z model is supported for loads') + raise ValueError('Only Z model is supported for shunts') self.s_n_sys = s_n_sys @@ -464,8 +567,6 @@ def __init__(self, param_dict, s_n_sys, v_n_sys, parallel_sims): z = torch.conj(1 / s_shunt) self.y_shunt = 1 / z - self.enable_parallel_simulation(parallel_sims) - def get_lf_power(self): """ Calculates and returns the load flow power of the Load. @@ -479,6 +580,7 @@ def get_lf_power(self): return torch.zeros_like(self.q_soll_mvar) + # noinspection PyUnusedLocal def get_admittance(self, dyn): """ Computes and returns the admittance of the Load. diff --git a/src/diffpssi/power_sim_lib/models/synchronous_machine.py b/src/diffpssi/power_sim_lib/models/synchronous_machine.py index 4597bfc..b51d6e7 100644 --- a/src/diffpssi/power_sim_lib/models/synchronous_machine.py +++ b/src/diffpssi/power_sim_lib/models/synchronous_machine.py @@ -2,12 +2,7 @@ The model of a generator (Synchronous Machine Model) is defined in this file. The model is a 6th order differential equation, which is solved during the simulation. """ -import torch - -from src.diffpssi.power_sim_lib.models.exciters import SEXS -from src.diffpssi.power_sim_lib.backend import * -from src.diffpssi.power_sim_lib.models.governors import TGOV1 -from src.diffpssi.power_sim_lib.models.stabilizers import STAB1 +from diffpssi.power_sim_lib.backend import * class SynchMachine(object): @@ -24,16 +19,78 @@ class SynchMachine(object): governor (TGOV1): The governor model associated with the machine. stabilizer (STAB1): The power system stabilizer model associated with the machine. """ - def __init__(self, param_dict, f_n_sys, s_n_sys, v_n_sys, parallel_sims): + + def __init__(self, f_n_sys, s_n_sys, v_n_sys, + param_dict=None, + name=None, + bus=None, + s_n=None, + v_n=None, + p=None, + v=None, + h=None, + d=None, + x_d=None, + x_q=None, + x_d_t=None, + x_q_t=None, + x_d_st=None, + x_q_st=None, + t_d0_t=None, + t_q0_t=None, + t_d0_st=None, + t_q0_st=None, + ): """ Initializes the SynchMachine object with specified parameters. Args: param_dict (dict, optional): A dictionary of parameters for the machine. - s_n, v_n, p_soll_mw, v_soll, h, d, x_d, x_q, x_d_t, x_q_t, x_d_st, x_q_st, t_d0_t, t_q0_t, t_d0_st, t_q0_st: - Electrical and mechanical characteristics of the machine. - parallel_sims (int, optional): Number of parallel simulations to enable. + name (str, optional): The name of the machine. + bus (Bus, optional): The bus the machine is connected to. + s_n (float, optional): The apparent power rating of the machine in MVA. + v_n (float, optional): The rated voltage of the machine in kV. + p (float, optional): The active power of the machine in MW. + v (float, optional): The voltage of the machine in per unit. + h (float, optional): The inertia constant of the machine in seconds. + d (float, optional): The damping coefficient of the machine in pu. + x_d (float, optional): The direct-axis transient reactance of the machine in pu. + x_q (float, optional): The quadrature-axis transient reactance of the machine in pu. + x_d_t (float, optional): The direct-axis transient reactance of the machine in pu. + x_q_t (float, optional): The quadrature-axis transient reactance of the machine in pu. + x_d_st (float, optional): The direct-axis subtransient reactance of the machine in pu. + x_q_st (float, optional): The quadrature-axis subtransient reactance of the machine in pu. + t_d0_t (float, optional): The direct-axis open-circuit time constant of the machine in seconds. + t_q0_t (float, optional): The quadrature-axis open-circuit time constant of the machine in seconds. + t_d0_st (float, optional): The direct-axis open-circuit subtransient time constant of the machine in + seconds. + t_q0_st (float, optional): The quadrature-axis open-circuit subtransient time constant of the machine in + seconds. """ + if param_dict is None: + param_dict = { + 'name': name, + 'bus': bus, + 'S_n': s_n, + 'V_n': v_n, + 'P': p, + 'V': v, + 'H': h, + 'D': d, + 'X_d': x_d, + 'X_q': x_q, + 'X_d_t': x_d_t, + 'X_q_t': x_q_t, + 'X_d_st': x_d_st, + 'X_q_st': x_q_st, + 'T_d0_t': t_d0_t, + 'T_q0_t': t_q0_t, + 'T_d0_st': t_d0_st, + 'T_q0_st': t_q0_st, + } + + self.parallel_sims = None + self.name = param_dict['name'] self.bus = param_dict['bus'] self.s_n = param_dict['S_n'] @@ -78,10 +135,6 @@ def __init__(self, param_dict, f_n_sys, s_n_sys, v_n_sys, parallel_sims): self.governor = None self.stabilizer = None - self.parallel_sims = parallel_sims - - self.enable_parallel_simulation(parallel_sims) - def differential(self): """ Computes the differential equations governing the dynamics of the synchronous machine. @@ -120,36 +173,33 @@ def differential(self): return deriv_vec - def add_exciter(self, exciter_dict, parallel_sims): + def add_exciter(self, exciter_model): """ Associates an exciter model with the synchronous machine. Args: - exciter_dict (dict): A dictionary containing parameters for initializing the exciter model. - parallel_sims (int): Number of parallel simulations to enable. + exciter_model: The exciter model to associate with the machine. """ - self.exciter = SEXS(param_dict=exciter_dict, parallel_sims=parallel_sims, v_setpoint=self.v_soll) + self.exciter = exciter_model - def add_governor(self, governor_dict, parallel_sims): + def add_governor(self, governor_model): """ Associates a governor model with the synchronous machine. Args: - governor_dict (dict): A dictionary containing parameters for initializing the governor model. - parallel_sims (int): Number of parallel simulations to enable. + governor_model: The governor model to associate with the machine. """ - self.governor = TGOV1(param_dict=governor_dict, parallel_sims=parallel_sims) + self.governor = governor_model - def add_pss(self, pss_dict, parallel_sims): + def add_pss(self, pss_model): """ Associates a power system stabilizer model with the synchronous machine. Args: - pss_dict (dict): A dictionary containing parameters for initializing the power system stabilizer model. - parallel_sims (int): Number of parallel simulations to enable. + pss_model: The power system stabilizer model to associate with the machine. """ - self.stabilizer = STAB1(param_dict=pss_dict, parallel_sims=parallel_sims) + self.stabilizer = pss_model def set_state_vector(self, x): """ @@ -207,7 +257,7 @@ def calc_current_injections(self): i_q = self.e_d_st / (1j * self.x_q_st) * torch.exp(1j * (self.delta - torch.pi / 2)) # transform it to the base system. sn/vn is local per unit and sys_s_n/sys_v_n is global per unit - i_inj = (i_d + i_q)*(self.s_n / self.s_n_sys) + i_inj = (i_d + i_q) * (self.s_n / self.s_n_sys) return i_inj def update_internal_vars(self, v_bb): @@ -237,7 +287,7 @@ def initialize(self, s_calc, v_bb): # First calculate the currents of the generator at the busbar. # Those currents can then be used to calculate all internal voltages. - i_gen = torch.conj(s_calc / v_bb)/(self.s_n / self.s_n_sys) + i_gen = torch.conj(s_calc / v_bb) / (self.s_n / self.s_n_sys) # Calculate the internal voltages and angle of the generator # Basically this is always U2 = U1 + jX * I @@ -260,7 +310,7 @@ def initialize(self, s_calc, v_bb): self.e_q_st = v_q + self.x_d_st * i_d self.e_d_st = v_d - self.x_q_st * i_q - self.p_m = s_calc.real/(self.s_n / self.s_n_sys) + self.p_m = s_calc.real / (self.s_n / self.s_n_sys) self.e_fd = self.e_q_t + i_d * (self.x_d - self.x_d_t) if self.exciter: @@ -297,7 +347,7 @@ def get_admittance(self, dyn): torch.Tensor: Admittance of the synchronous machine. """ if dyn: - return -1j/self.x_d_st * (self.s_n / self.s_n_sys) + return -1j / self.x_d_st * (self.s_n / self.s_n_sys) else: return torch.zeros((self.parallel_sims, 1), dtype=torch.complex128) @@ -308,7 +358,7 @@ def get_lf_power(self): Returns: torch.Tensor: The load flow power. """ - return (self.p_soll_mw + 1j * 0)/self.s_n_sys + return (self.p_soll_mw + 1j * 0) / self.s_n_sys def enable_parallel_simulation(self, parallel_sims): """ @@ -317,35 +367,36 @@ def enable_parallel_simulation(self, parallel_sims): Args: parallel_sims (int): Number of parallel simulations. """ + self.parallel_sims = parallel_sims - self.s_n = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.s_n - self.v_n = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.v_n + self.s_n = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.s_n + self.v_n = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.v_n self.p_soll_mw = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.p_soll_mw - self.v_soll = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.v_soll - - self.h = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.h - self.d = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.d - self.x_d = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_d - self.x_q = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_q - self.x_d_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_d_t - self.x_q_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_q_t - self.x_d_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_d_st - self.x_q_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.x_q_st - self.t_d0_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.t_d0_t - self.t_q0_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.t_q0_t - self.t_d0_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.t_d0_st - self.t_q0_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.t_q0_st - - self.omega = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.omega - self.delta = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.delta - self.e_q_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.e_q_t - self.e_d_t = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.e_d_t - self.e_q_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.e_q_st - self.e_d_st = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.e_d_st - - self.p_m = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.p_m - self.p_e = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.p_e - self.e_fd = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.e_fd - self.i_d = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.i_d - self.i_q = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.i_q - self.v_bb = torch.ones((parallel_sims, 1), dtype=torch.complex128)*self.v_bb + self.v_soll = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.v_soll + + self.h = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.h + self.d = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.d + self.x_d = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_d + self.x_q = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_q + self.x_d_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_d_t + self.x_q_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_q_t + self.x_d_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_d_st + self.x_q_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.x_q_st + self.t_d0_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.t_d0_t + self.t_q0_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.t_q0_t + self.t_d0_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.t_d0_st + self.t_q0_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.t_q0_st + + self.omega = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.omega + self.delta = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.delta + self.e_q_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.e_q_t + self.e_d_t = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.e_d_t + self.e_q_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.e_q_st + self.e_d_st = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.e_d_st + + self.p_m = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.p_m + self.p_e = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.p_e + self.e_fd = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.e_fd + self.i_d = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.i_d + self.i_q = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.i_q + self.v_bb = torch.ones((parallel_sims, 1), dtype=torch.complex128) * self.v_bb diff --git a/src/diffpssi/power_sim_lib/simulator.py b/src/diffpssi/power_sim_lib/simulator.py index f57babe..4bf020f 100644 --- a/src/diffpssi/power_sim_lib/simulator.py +++ b/src/diffpssi/power_sim_lib/simulator.py @@ -8,11 +8,14 @@ import numpy as np from tqdm import tqdm -from src.diffpssi.power_sim_lib.load_flow import do_load_flow -from src.diffpssi.power_sim_lib.models.synchronous_machine import SynchMachine -from src.diffpssi.power_sim_lib.models.static_models import * -from src.diffpssi.power_sim_lib.backend import * -from src.diffpssi.power_sim_lib.solvers import solver_dict +from diffpssi.power_sim_lib.load_flow import do_load_flow +from diffpssi.power_sim_lib.models.synchronous_machine import SynchMachine +from diffpssi.power_sim_lib.models.static_models import * +from diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.solvers import solver_dict +from diffpssi.power_sim_lib.models.governors import TGOV1 +from diffpssi.power_sim_lib.models.stabilizers import STAB1 +from diffpssi.power_sim_lib.models.exciters import SEXS class PowerSystemSimulation(object): @@ -52,12 +55,13 @@ class PowerSystemSimulation(object): reset: Resets the simulation to its initial state. run: Runs the simulation. """ + def __init__(self, time_step, sim_time, parallel_sims, solver, - grid_data, + grid_data=None, verbose=True, ): """ @@ -88,18 +92,20 @@ def __init__(self, self.record_func = None self.verbose = verbose - try: + if os.environ.get('DIFFPSSI_FORCE_INTEGRATOR') is not None: # this should only be used for integration tests - self.solver = solver_dict[os.environ['DIFFPSSI_FORCE_INTEGRATOR']]() + self.solver = solver_dict[os.environ.get('DIFFPSSI_FORCE_INTEGRATOR')]() print('WARNING: FORCING THE USE OF THE {} INTEGRATOR. ' - 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ['DIFFPSSI_FORCE_INTEGRATOR'])) - except KeyError: + 'THIS SHOULD ONLY HAPPEN FOR UNITTESTS'.format(os.environ.get('DIFFPSSI_FORCE_INTEGRATOR'))) + else: self.solver = solver_dict[solver]() - self.fn = grid_data['f'] - self.base_mva = grid_data['base_mva'] - self.base_voltage = grid_data['base_voltage'] - self.create_grid(grid_data) + self.base_voltage = None + self.base_mva = None + self.fn = None + + if grid_data is not None: + self.create_grid(grid_data) self.backend = BACKEND @@ -125,6 +131,10 @@ def create_grid(self, grid_data): Args: grid_data (dict): Dictionary containing the grid data. """ + self.fn = grid_data['f'] + self.base_mva = grid_data['base_mva'] + self.base_voltage = grid_data['base_voltage'] + transformed_data = {} # first transform the data to a more convenient format for key, value in grid_data.items(): @@ -142,117 +152,192 @@ def create_grid(self, grid_data): transformed_data[key] = value grid_data = transformed_data - for bus in grid_data['busses']: - self.add_bus(bus) + for bus_dict in grid_data['busses']: + bus_model = Bus(param_dict=bus_dict) + bus_model.enable_parallel_simulation(self.parallel_sims) + self.add_bus(bus_model) generators = grid_data.get('generators', []) if isinstance(generators, dict): - for model in generators['GEN']: - self.add_generator(model) - - for load in grid_data.get('loads', []): - self.add_load(load) - - for shunt in grid_data.get('shunts', []): - self.add_shunt(shunt) - - for line in grid_data.get('lines', []): - self.add_line(line) - - for transformer in grid_data.get('transformers', []): - self.add_transformer(transformer) + for gen_dict in generators['GEN']: + generator_model = SynchMachine(param_dict=gen_dict, + s_n_sys=self.base_mva, + v_n_sys=self.base_voltage, + f_n_sys=self.fn) + self.add_generator(generator_model) + + for load_dict in grid_data.get('loads', []): + load_model = Load(param_dict=load_dict, + s_n_sys=self.base_mva) + self.add_load(load_model) + + for shunt_dict in grid_data.get('shunts', []): + shunt_model = Shunt(param_dict=shunt_dict, + s_n_sys=self.base_mva) + self.add_shunt(shunt_model) + + for line_dict in grid_data.get('lines', []): + line_model = Line(param_dict=line_dict, + s_n_sys=self.base_mva, + v_n_sys=self.base_voltage) + self.add_line(line_model) + + for transformer_dict in grid_data.get('transformers', []): + transformer_model = Transformer(param_dict=transformer_dict, + s_n_sys=self.base_mva) + self.add_transformer(transformer_model) exciters = grid_data.get('avr', []) if isinstance(exciters, dict): - for model in exciters['SEXS']: - self.get_generator_by_name(model['gen']).add_exciter(model, self.parallel_sims) + for sexs_dict in exciters['SEXS']: + exciter_model = SEXS(param_dict=sexs_dict) + self.add_exciter(exciter_model) governors = grid_data.get('gov', []) if isinstance(governors, dict): - for model in governors['TGOV1']: - self.get_generator_by_name(model['gen']).add_governor(model, self.parallel_sims) + for gov_dict in governors['TGOV1']: + gov_model = TGOV1(param_dict=gov_dict) + self.add_governor(gov_model) psss = grid_data.get('pss', []) if isinstance(psss, dict): - for model in psss['STAB1']: - self.get_generator_by_name(model['gen']).add_pss(model, self.parallel_sims) + for pss_dict in psss['STAB1']: + pss_model = STAB1(param_dict=pss_dict) + self.add_pss(pss_model) - self.busses[self.bus_idxs[grid_data['slack_bus']]].lf_type = 'SL' + self.set_slack_bus(grid_data['slack_bus']) - def add_bus(self, param_dict): + def set_slack_bus(self, slack_bus): + """ + Sets the slack bus of the system. + + Args: + slack_bus (str): The name of the slack bus. + """ + slack_bus_idx = self.bus_idxs[slack_bus] + self.busses[slack_bus_idx].lf_type = 'SL' + + def add_bus(self, bus_model): """ Adds a bus to the system. Args: - param_dict (dict): Dictionary containing the parameters of the bus. + bus_model (Bus): A bus model to add to the grid. """ # Assign an index to the bus and add it to the list of buses - self.bus_idxs[param_dict['name']] = len(self.busses) - self.busses.append(Bus(param_dict, self.parallel_sims)) + self.bus_idxs[bus_model.name] = len(self.busses) + bus_model.enable_parallel_simulation(self.parallel_sims) + self.busses.append(bus_model) - def add_generator(self, param_dict): + def add_generator(self, generator_model): """ Adds a generator to a specified bus in the system. Args: - param_dict (dict): Dictionary containing the parameters of the generator. + generator_model (SynchMachine): A generator model to add to the grid. """ - bus = self.busses[self.bus_idxs[param_dict['bus']]] - bus.add_model(SynchMachine(param_dict, self.fn, self.base_mva, bus.v_n, self.parallel_sims)) + bus = self.busses[self.bus_idxs[generator_model.bus]] + generator_model.enable_parallel_simulation(self.parallel_sims) + bus.add_model(generator_model) # fit the voltage of the bus to the generator bus.update_voltages(bus.models[-1].v_soll) bus.lf_type = 'PV' - def add_load(self, param_dict): + def add_load(self, load_model): """ Adds a load to a specified bus in the system. Args: - param_dict (dict): Dictionary containing the parameters of the load. + load_model (Load): A load model to add to the grid. """ - bus = self.busses[self.bus_idxs[param_dict['bus']]] - bus.add_model(Load(param_dict, self.base_mva, bus.v_n, self.parallel_sims)) + bus = self.busses[self.bus_idxs[load_model.bus]] + load_model.enable_parallel_simulation(self.parallel_sims) + bus.add_model(load_model) - def add_shunt(self, param_dict): + def add_shunt(self, shunt_model): """ - Adds a load to a specified bus in the system. + Adds a shunt to a specified bus in the system. Args: - param_dict (dict): Dictionary containing the parameters of the shunt element. + shunt_model (Shunt): A shunt model to add to the grid. """ - bus = self.busses[self.bus_idxs[param_dict['bus']]] - bus.add_model(Shunt(param_dict, self.base_mva, bus.v_n, self.parallel_sims)) + bus = self.busses[self.bus_idxs[shunt_model.bus]] + shunt_model.enable_parallel_simulation(self.parallel_sims) + bus.add_model(shunt_model) - def add_line(self, param_dict): + def add_line(self, line_model): """ Adds a transmission line between two buses in the system. Args: - param_dict (dict): Dictionary containing the parameters of the line. + line_model (Line): A line model to add to the grid. """ - bus_from = self.bus_idxs[param_dict['from_bus']] - bus_to = self.bus_idxs[param_dict['to_bus']] - self.lines.append(Line(param_dict, - bus_from, - bus_to, - self.base_mva, - self.base_voltage, - parallel_sims=self.parallel_sims)) + bus_from = self.bus_idxs[line_model.from_bus_name] + bus_to = self.bus_idxs[line_model.to_bus_name] + + line_model.from_bus_id = bus_from + line_model.to_bus_id = bus_to + + line_model.enable_parallel_simulation(self.parallel_sims) - def add_transformer(self, param_dict): + self.lines.append(line_model) + + def add_transformer(self, transformer_model): """ Adds a transformer between two buses in the system. Args: - param_dict (dict): Dictionary containing the parameters of the transformer. + transformer_model (Transformer): A transformer model to add to the grid. + """ + bus_from = self.bus_idxs[transformer_model.from_bus_name] + bus_to = self.bus_idxs[transformer_model.to_bus_name] + + transformer_model.from_bus_id = bus_from + transformer_model.to_bus_id = bus_to + + transformer_model.enable_parallel_simulation(self.parallel_sims) + + self.trafos.append(transformer_model) + + def add_exciter(self, exciter_model): + """ + Adds an exciter to a specified generator in the system. + Different exciter models work, for example the SEXS. + + Args: + exciter_model (object): An exciter model to add to the grid. + """ + generator = self.get_generator_by_name(exciter_model.gen) + exciter_model.v_setpoint = generator.v_soll + exciter_model.enable_parallel_simulation(self.parallel_sims) + + generator.add_exciter(exciter_model) + + def add_governor(self, governor_model): + """ + Adds a governor to a specified generator in the system. + Different governor models work, for example the TGOV1. + + Args: + governor_model (object): A governor model to add to the grid. """ - bus_from = self.bus_idxs[param_dict['from_bus']] - bus_to = self.bus_idxs[param_dict['to_bus']] - self.trafos.append(Transfomer(param_dict, - bus_from, - bus_to, - s_n_sys=self.base_mva, - parallel_sims=self.parallel_sims)) + generator = self.get_generator_by_name(governor_model.gen) + governor_model.enable_parallel_simulation(self.parallel_sims) + + generator.add_governor(governor_model) + + def add_pss(self, pss_model): + """ + Adds a PSS to a specified generator in the system. + Different PSS models work, for example the STAB1. + + Args: + pss_model (object): A PSS model to add to the grid. + """ + generator = self.get_generator_by_name(pss_model.gen) + pss_model.enable_parallel_simulation(self.parallel_sims) + + generator.add_pss(pss_model) def admittance_matrix(self, dynamic): """ @@ -277,19 +362,19 @@ def admittance_matrix(self, dynamic): len(self.busses)), dtype=torch.complex128) for line in self.lines: - self.dynamic_y_matrix[:, line.from_bus, line.to_bus] += line.get_admittance_off_diagonal() - self.dynamic_y_matrix[:, line.to_bus, line.from_bus] += line.get_admittance_off_diagonal() - self.dynamic_y_matrix[:, line.from_bus, line.from_bus] += line.get_admittance_diagonal() - self.dynamic_y_matrix[:, line.to_bus, line.to_bus] += line.get_admittance_diagonal() + self.dynamic_y_matrix[:, line.from_bus_id, line.to_bus_id] += line.get_admittance_off_diagonal() + self.dynamic_y_matrix[:, line.to_bus_id, line.from_bus_id] += line.get_admittance_off_diagonal() + self.dynamic_y_matrix[:, line.from_bus_id, line.from_bus_id] += line.get_admittance_diagonal() + self.dynamic_y_matrix[:, line.to_bus_id, line.to_bus_id] += line.get_admittance_diagonal() for transformer in self.trafos: - self.dynamic_y_matrix[:, transformer.from_bus, transformer.to_bus] += ( + self.dynamic_y_matrix[:, transformer.from_bus_id, transformer.to_bus_id] += ( transformer.get_admittance_off_diagonal()) - self.dynamic_y_matrix[:, transformer.to_bus, transformer.from_bus] += ( + self.dynamic_y_matrix[:, transformer.to_bus_id, transformer.from_bus_id] += ( transformer.get_admittance_off_diagonal()) - self.dynamic_y_matrix[:, transformer.from_bus, transformer.from_bus] += ( + self.dynamic_y_matrix[:, transformer.from_bus_id, transformer.from_bus_id] += ( transformer.get_admittance_diagonal()) - self.dynamic_y_matrix[:, transformer.to_bus, transformer.to_bus] += ( + self.dynamic_y_matrix[:, transformer.to_bus_id, transformer.to_bus_id] += ( transformer.get_admittance_diagonal()) for i, bus in enumerate(self.busses): @@ -303,19 +388,19 @@ def admittance_matrix(self, dynamic): len(self.busses)), dtype=torch.complex128) for line in self.lines: - self.static_y_matrix[:, line.from_bus, line.to_bus] += line.get_admittance_off_diagonal() - self.static_y_matrix[:, line.to_bus, line.from_bus] += line.get_admittance_off_diagonal() - self.static_y_matrix[:, line.from_bus, line.from_bus] += line.get_admittance_diagonal() - self.static_y_matrix[:, line.to_bus, line.to_bus] += line.get_admittance_diagonal() + self.static_y_matrix[:, line.from_bus_id, line.to_bus_id] += line.get_admittance_off_diagonal() + self.static_y_matrix[:, line.to_bus_id, line.from_bus_id] += line.get_admittance_off_diagonal() + self.static_y_matrix[:, line.from_bus_id, line.from_bus_id] += line.get_admittance_diagonal() + self.static_y_matrix[:, line.to_bus_id, line.to_bus_id] += line.get_admittance_diagonal() for transformer in self.trafos: - self.static_y_matrix[:, transformer.from_bus, transformer.to_bus] += ( + self.static_y_matrix[:, transformer.from_bus_id, transformer.to_bus_id] += ( transformer.get_admittance_off_diagonal()) - self.static_y_matrix[:, transformer.to_bus, transformer.from_bus] += ( + self.static_y_matrix[:, transformer.to_bus_id, transformer.from_bus_id] += ( transformer.get_admittance_off_diagonal()) - self.static_y_matrix[:, transformer.from_bus, transformer.from_bus] += ( + self.static_y_matrix[:, transformer.from_bus_id, transformer.from_bus_id] += ( transformer.get_admittance_diagonal()) - self.static_y_matrix[:, transformer.to_bus, transformer.to_bus] += ( + self.static_y_matrix[:, transformer.to_bus_id, transformer.to_bus_id] += ( transformer.get_admittance_diagonal()) for i, bus in enumerate(self.busses): diff --git a/src/diffpssi/power_sim_lib/solvers.py b/src/diffpssi/power_sim_lib/solvers.py index c1f21a1..904685e 100644 --- a/src/diffpssi/power_sim_lib/solvers.py +++ b/src/diffpssi/power_sim_lib/solvers.py @@ -2,7 +2,7 @@ File contains all the solvers for the power system simulation. The solvers are used to integrate the differential equations of the power system simulation. More solvers can be added here """ -from src.diffpssi.power_sim_lib.backend import * +from diffpssi.power_sim_lib.backend import * class Euler(object): diff --git a/tests/integration/test_correct_backend.py b/tests/integration/test_correct_backend.py new file mode 100644 index 0000000..31aa218 --- /dev/null +++ b/tests/integration/test_correct_backend.py @@ -0,0 +1,14 @@ +""" +Unit tests for the simulation examples in the examples folder. +""" +import unittest + + +class TestCorrectBackend(unittest.TestCase): + def test_correct_backend(self): + """ + This test checks if the correct backend is used. + """ + import diffpssi.power_sim_lib.backend as backend + + self.assertEqual(backend.BACKEND, 'numpy') diff --git a/tests/integration/test_example_opts.py b/tests/integration/test_example_opts.py index 699359d..7ce80f3 100644 --- a/tests/integration/test_example_opts.py +++ b/tests/integration/test_example_opts.py @@ -5,15 +5,22 @@ import unittest import numpy as np -# set environment variable that forces the simulator to always use the heun integrator -os.environ['DIFFPSSI_FORCE_INTEGRATOR'] = 'heun' -os.environ['DIFFPSSI_FORCE_SIM_BACKEND'] = 'torch' -os.environ['DIFFPSSI_FORCE_OPT_ITERS'] = '1' -os.environ['DIFFPSSI_FORCE_PARALLEL_SIMS'] = '2' - -class TestCustomMultiMachine(unittest.TestCase): - def test_simulation_output(self): +def set_env_vars(): + """ + Set the environment variables that are used to force the use of the heun integrator, the use of torch, and the + number of parallel simulations. + """ + # set environment variable that forces the simulator to always use the heun integrator + os.environ['DIFFPSSI_FORCE_INTEGRATOR'] = 'heun' + os.environ['DIFFPSSI_FORCE_SIM_BACKEND'] = 'torch' + os.environ['DIFFPSSI_FORCE_OPT_ITERS'] = '1' + os.environ['DIFFPSSI_FORCE_PARALLEL_SIMS'] = '2' + + +class TestOptRunAbility(unittest.TestCase): + def test_custom_multi_machine_opt(self): + set_env_vars() import examples.models.custom_multi_machine.custom_multi_machine_opt as custom_multi_machine_opt # get the path of the custom_multi_machine_sim.py file @@ -24,9 +31,8 @@ def test_simulation_output(self): # Run the optimization custom_multi_machine_opt.main(parallel_sims=2) - -class TestIbbModel(unittest.TestCase): - def test_simulation_output(self): + def test_ibb_opt(self): + set_env_vars() import examples.models.ibb_model.ibb_opt as ibb_opt # get the path of the ibb_sim.py file @@ -37,9 +43,8 @@ def test_simulation_output(self): # Run the optimization ibb_opt.main(parallel_sims=2) - -class TestIeee9Bus(unittest.TestCase): - def test_simulation_output(self): + def test_ieee_9bus_opt(self): + set_env_vars() import examples.models.ieee_9bus.ieee_9bus_opt as ieee_9bus_opt # get the path of the ibb_sim.py file @@ -50,9 +55,7 @@ def test_simulation_output(self): # Run the optimization ieee_9bus_opt.main(parallel_sims=2) - -class TestK2A(unittest.TestCase): - def test_simulation_output(self): + def test_k2a_opt(self): import examples.models.k2a.k2a_opt as k2a_opt # get the path of the ibb_sim.py file diff --git a/tests/integration/test_example_sims.py b/tests/integration/test_example_sims.py index 36d510e..365ecb4 100644 --- a/tests/integration/test_example_sims.py +++ b/tests/integration/test_example_sims.py @@ -5,15 +5,21 @@ import unittest import numpy as np -# set environment variable that forces the simulator to always use the heun integrator -os.environ['DIFFPSSI_FORCE_INTEGRATOR'] = 'heun' -atol = 1e-8 -rtol = 1e-8 +class TestSimulationOutput(unittest.TestCase): + def set_env_vars(self): + """ + Set the environment variables for the tests. + """ + # set environment variable that forces the simulator to always use the heun integrator + os.environ['DIFFPSSI_FORCE_INTEGRATOR'] = 'heun' + os.environ['DIFFPSSI_TESTING'] = 'True' + self.atol = 1e-8 + self.rtol = 1e-8 -class TestCustomMultiMachine(unittest.TestCase): - def test_simulation_output(self): + def test_custom_multi_machine_sim(self): + self.set_env_vars() import examples.models.custom_multi_machine.custom_multi_machine_sim as custom_multi_machine_sim # get the path of the custom_multi_machine_sim.py file @@ -34,11 +40,34 @@ def test_simulation_output(self): expected_output = np.load('./data/custom_multi_machine.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) + def test_ibb_manual_sim(self): + self.set_env_vars() + import examples.models.ibb_manual.ibb_manual_sim as ibb_manual_sim -class TestIbbModel(unittest.TestCase): - def test_simulation_output(self): + # get the path of the ibb_manual_sim.py file + file_path = ibb_manual_sim.__file__ + + os.chdir(os.path.dirname(file_path)) + + # Run the simulation + ibb_manual_sim.main() + + # Load the output file and compare it to the expected output + output_sim = ibb_manual_sim.np.load('./data/original_data.npy') + + # change the directory back to where this file here is + os.chdir(os.path.dirname(os.path.abspath(__file__))) + + # Load the expected output in the tests folder + expected_output = np.load('./data/ibb_model.npy') + + # compare the two outputs + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) + + def test_ibb_model_sim(self): + self.set_env_vars() import examples.models.ibb_model.ibb_sim as ibb_sim # get the path of the ibb_sim.py file @@ -59,11 +88,10 @@ def test_simulation_output(self): expected_output = np.load('./data/ibb_model.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) - + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) -class TestIbbTransformer(unittest.TestCase): - def test_simulation_output(self): + def test_ibb_transformer_sim(self): + self.set_env_vars() import examples.models.ibb_transformer.ibb_trans_sim as ibb_trans_sim # get the path of the ibb_trans_sim.py file @@ -84,11 +112,10 @@ def test_simulation_output(self): expected_output = np.load('./data/ibb_transformer.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) - -class TestIbbWithControllers(unittest.TestCase): - def test_simulation_output(self): + def test_ibb_with_controllers_sim(self): + self.set_env_vars() import examples.models.ibb_with_controllers.ibb_wc_sim as ibb_wc_sim # get the path of the ibb_wc_sim.py file @@ -109,11 +136,10 @@ def test_simulation_output(self): expected_output = np.load('./data/ibb_with_controllers.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) - + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) -class TestIEEE9Bus(unittest.TestCase): - def test_simulation_output(self): + def test_ieee_9bus_sim(self): + self.set_env_vars() import examples.models.ieee_9bus.ieee_9bus_sim as ieee_9bus_sim # get the path of the ieee_9bus_sim.py file @@ -134,11 +160,10 @@ def test_simulation_output(self): expected_output = np.load('./data/ieee_9bus.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) - -class TestK2A(unittest.TestCase): - def test_simulation_output(self): + def test_k2a_sim(self): + self.set_env_vars() import examples.models.k2a.k2a_sim as k2a_sim # get the path of the k2a_sim.py file @@ -159,4 +184,28 @@ def test_simulation_output(self): expected_output = np.load('./data/k2a.npy') # compare the two outputs - self.assertTrue(np.allclose(output_sim, expected_output, atol=atol, rtol=rtol)) + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol)) + + def test_k2a_manual_sim(self): + self.set_env_vars() + import examples.models.k2a_manual.k2a_manual_sim as k2a_manual_sim + + # get the path of the k2a_manual_sim.py file + file_path = k2a_manual_sim.__file__ + + os.chdir(os.path.dirname(file_path)) + + # Run the simulation + k2a_manual_sim.main() + + # Load the output file and compare it to the expected output + output_sim = k2a_manual_sim.np.load('./data/original_data.npy') + + # change the directory back to where this file here is + os.chdir(os.path.dirname(os.path.abspath(__file__))) + + # Load the expected output in the tests folder + expected_output = np.load('./data/k2a.npy') + + # compare the two outputs + self.assertTrue(np.allclose(output_sim, expected_output, atol=self.atol, rtol=self.rtol))