Add .ipynb tests and demo files for eigenvalue extraction

Signed-off-by: Georgii Tishenin <[email protected]>
sogno-platform · May 19, 2024 · ff4b7d9 · ff4b7d9
1 parent e9e2848
commit ff4b7d9
Show file tree

Hide file tree

Showing 2 changed files with 390 additions and 0 deletions.
diff --git a/examples/Notebooks/Eigenvalues/DP_EMT_SinglePhaseRLC_Tests.ipynb b/examples/Notebooks/Eigenvalues/DP_EMT_SinglePhaseRLC_Tests.ipynb
@@ -0,0 +1,226 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Eigenvalue extraction tests in a single-phase RLC circuit in DP & EMT domains"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dpsimpy\n",
+    "import math\n",
+    "import villas.dataprocessing.readtools as rt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# function to assert that two complex numbers are close\n",
+    "def assert_complex_isclose(expected, actual, rel_tol=1e-6):\n",
+    "    assert math.isclose(expected.real, actual.real, rel_tol=rel_tol), \"Real parts not close: {} vs {}\".format(expected.real, actual.real)\n",
+    "    assert math.isclose(expected.imag, actual.imag, rel_tol=rel_tol), \"Imaginary parts not close: {} vs {}\".format(expected.imag, actual.imag)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Circuit parameters\n",
+    "v_volt = 8.5\n",
+    "r_ohm = 100\n",
+    "l_henry = 5\n",
+    "c_farad = 250e-6\n",
+    "deltaT = 1e-4\n",
+    "final_time = 1e-3"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Analytical solution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "eigenvalue0_expected = ((c_farad**2*r_ohm**2 - 4*l_henry*c_farad)**(1/2) - c_farad*r_ohm)/(2*c_farad*l_henry)\n",
+    "eigenvalue1_expected = -((c_farad**2*r_ohm**2 - 4*l_henry*c_farad)**(1/2) + c_farad*r_ohm)/(2*c_farad*l_henry)\n",
+    "\n",
+    "print('Expected eigenvalue 0: ' + str(eigenvalue0_expected))\n",
+    "print('Expected eigenvalue 1: ' + str(eigenvalue1_expected))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Test 1: Eigenvalues extracted from simulation in DP domain match analytical solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# DPsim DP simulation\n",
+    "name = 'DP_SinglePhaseRLC_Test'\n",
+    "log_dir = \"logs/\" + name\n",
+    "dpsimpy.Logger.set_log_dir(log_dir)\n",
+    "\n",
+    "# Create nodes\n",
+    "gnd = dpsimpy.dp.SimNode.gnd\n",
+    "n1 = dpsimpy.dp.SimNode('n1')\n",
+    "n2 = dpsimpy.dp.SimNode('n2')\n",
+    "n3 = dpsimpy.dp.SimNode('n3')\n",
+    "\n",
+    "# Create components\n",
+    "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n",
+    "vs.set_parameters(V_ref=complex(v_volt,0))\n",
+    "r = dpsimpy.dp.ph1.Resistor('r')\n",
+    "r.set_parameters(R=r_ohm)\n",
+    "l = dpsimpy.dp.ph1.Inductor('l')\n",
+    "l.set_parameters(L=l_henry)\n",
+    "c = dpsimpy.dp.ph1.Capacitor('c')\n",
+    "c.set_parameters(C=c_farad)\n",
+    "\n",
+    "# Set connections\n",
+    "vs.connect([gnd, n1])\n",
+    "l.connect([n1, n2])\n",
+    "c.connect([n2, n3])\n",
+    "r.connect([n3, gnd])\n",
+    "\n",
+    "# Create topology\n",
+    "system = dpsimpy.SystemTopology(50, [n1, n2, n3], [vs, r, l, c])\n",
+    "\n",
+    "# Configure and run simulation\n",
+    "sim = dpsimpy.Simulation(name)\n",
+    "sim.set_domain(dpsimpy.Domain.DP)\n",
+    "sim.set_system(system)\n",
+    "sim.do_eigenvalue_extraction(True)\n",
+    "sim.set_time_step(deltaT)\n",
+    "sim.set_final_time(final_time)\n",
+    "sim.run()\n",
+    "\n",
+    "# Read log file\n",
+    "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n",
+    "eigenvalue0 = eigenvalues_timeseries['eigenvalues_0'].values\n",
+    "eigenvalue1 = eigenvalues_timeseries['eigenvalues_1'].values\n",
+    "\n",
+    "# Assert\n",
+    "nSamples = int(final_time/deltaT)\n",
+    "\n",
+    "assert len(eigenvalue0) == nSamples\n",
+    "assert len(eigenvalue1) == nSamples\n",
+    "\n",
+    "for i in range(nSamples):\n",
+    "    assert_complex_isclose(eigenvalue0_expected, eigenvalue0[i], 1e-8)\n",
+    "    assert_complex_isclose(eigenvalue1_expected, eigenvalue1[i], 1e-8)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Test 2: Eigenvalues extracted from simulation in EMT domain match analytical solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# DPsim EMT simulation\n",
+    "name = 'EMT_SinglePhaseRLC_Test'\n",
+    "log_dir = \"logs/\" + name\n",
+    "dpsimpy.Logger.set_log_dir(log_dir)\n",
+    "\n",
+    "# Create nodes\n",
+    "gnd = dpsimpy.emt.SimNode.gnd\n",
+    "n1 = dpsimpy.emt.SimNode('n1')\n",
+    "n2 = dpsimpy.emt.SimNode('n2')\n",
+    "n3 = dpsimpy.emt.SimNode('n3')\n",
+    "\n",
+    "# Create components\n",
+    "vs = dpsimpy.emt.ph1.VoltageSource('vs')\n",
+    "vs.set_parameters(V_ref=complex(v_volt,0))\n",
+    "r = dpsimpy.emt.ph1.Resistor('r')\n",
+    "r.set_parameters(R=r_ohm)\n",
+    "l = dpsimpy.emt.ph1.Inductor('l')\n",
+    "l.set_parameters(L=l_henry)\n",
+    "c = dpsimpy.emt.ph1.Capacitor('c')\n",
+    "c.set_parameters(C=c_farad)\n",
+    "\n",
+    "# Set connections\n",
+    "vs.connect([gnd, n1])\n",
+    "l.connect([n1, n2])\n",
+    "c.connect([n2, n3])\n",
+    "r.connect([n3, gnd])\n",
+    "\n",
+    "# Create topology\n",
+    "system = dpsimpy.SystemTopology(50, [n1, n2, n3], [vs, r, l, c])\n",
+    "\n",
+    "# Configure and run simulation\n",
+    "sim = dpsimpy.Simulation(name)\n",
+    "sim.set_domain(dpsimpy.Domain.EMT)\n",
+    "sim.set_system(system)\n",
+    "sim.do_eigenvalue_extraction(True)\n",
+    "sim.set_time_step(deltaT)\n",
+    "sim.set_final_time(final_time)\n",
+    "sim.run()\n",
+    "\n",
+    "# Read log file\n",
+    "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n",
+    "eigenvalue0 = eigenvalues_timeseries['eigenvalues_0'].values\n",
+    "eigenvalue1 = eigenvalues_timeseries['eigenvalues_1'].values\n",
+    "\n",
+    "# Assert\n",
+    "nSamples = int(final_time/deltaT)\n",
+    "\n",
+    "assert len(eigenvalue0) == nSamples\n",
+    "assert len(eigenvalue1) == nSamples\n",
+    "\n",
+    "for i in range(nSamples):\n",
+    "    assert_complex_isclose(eigenvalue0_expected, eigenvalue0[i], 1e-8)\n",
+    "    assert_complex_isclose(eigenvalue1_expected, eigenvalue1[i], 1e-8)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/Notebooks/Eigenvalues/DP_SinglePhaseRLC_Demo.ipynb b/examples/Notebooks/Eigenvalues/DP_SinglePhaseRLC_Demo.ipynb
@@ -0,0 +1,164 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Eigenvalue extraction during simulation of a single-phase RLC circuit.\n",
+    "\n",
+    "- Simulate RLC circuit in DP domain.\n",
+    "- Simulation includes switch events changing resistance value of the circuit.\n",
+    "- Extract eigenvalues at each time step of simulation.\n",
+    "- Visualize eigenvalue trajectory."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dpsimpy\n",
+    "import villas.dataprocessing.readtools as rt\n",
+    "import villas.dataprocessing.plottools as pt\n",
+    "from villas.dataprocessing.timeseries import TimeSeries as ts\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# DPsim DP simulation\n",
+    "name = 'DP_SinglePhaseRLC_Demo'\n",
+    "log_dir = \"logs/\" + name\n",
+    "dpsimpy.Logger.set_log_dir(log_dir)\n",
+    "\n",
+    "# Create nodes\n",
+    "gnd = dpsimpy.dp.SimNode.gnd\n",
+    "n1 = dpsimpy.dp.SimNode('n1')\n",
+    "n2 = dpsimpy.dp.SimNode('n2')\n",
+    "n3 = dpsimpy.dp.SimNode('n3')\n",
+    "n4 = dpsimpy.dp.SimNode('n4')\n",
+    "n5 = dpsimpy.dp.SimNode('n5')\n",
+    "\n",
+    "# Create components\n",
+    "vs = dpsimpy.dp.ph1.VoltageSource('vs')\n",
+    "vs.set_parameters(V_ref=complex(100,0))\n",
+    "r1 = dpsimpy.dp.ph1.Resistor('r_1')\n",
+    "r1.set_parameters(R=10)\n",
+    "r2 = dpsimpy.dp.ph1.Resistor('r_2')\n",
+    "r2.set_parameters(R=0.1)\n",
+    "l = dpsimpy.dp.ph1.Inductor('l')\n",
+    "l.set_parameters(L=5)\n",
+    "c = dpsimpy.dp.ph1.Capacitor('c')\n",
+    "c.set_parameters(C=250e-6)\n",
+    "s1 = dpsimpy.dp.ph1.Switch('s1')\n",
+    "s1.set_parameters(1e8, 1e-4, True)\n",
+    "s2 = dpsimpy.dp.ph1.Switch('s2')\n",
+    "s2.set_parameters(1e8, 1e-4, False)\n",
+    "\n",
+    "# Create switch events\n",
+    "s1Event = dpsimpy.event.SwitchEvent(0.05, s1, False)\n",
+    "s2Event = dpsimpy.event.SwitchEvent(0.05, s2, True)\n",
+    "\n",
+    "# Set connections\n",
+    "vs.connect([gnd, n1])\n",
+    "r1.connect([n4, gnd])\n",
+    "r2.connect([n5, gnd])\n",
+    "l.connect([n1, n2])\n",
+    "c.connect([n2, n3])\n",
+    "s1.connect([n3, n4])\n",
+    "s2.connect([n3, n5])\n",
+    "\n",
+    "# Create topology\n",
+    "system = dpsimpy.SystemTopology(50, [n1, n2, n3, n4, n5], [vs, r1, r2, l, c, s1, s2])\n",
+    "\n",
+    "# Configure logging\n",
+    "logger = dpsimpy.Logger(name)\n",
+    "logger.log_attribute('current', 'i_intf', l)\n",
+    "logger.log_attribute('voltage', 'v', n3)\n",
+    "\n",
+    "# Configure and run simulation\n",
+    "sim = dpsimpy.Simulation(name)\n",
+    "sim.set_domain(dpsimpy.Domain.DP)\n",
+    "sim.set_system(system)\n",
+    "sim.add_event(s1Event)\n",
+    "sim.add_event(s2Event)\n",
+    "sim.do_eigenvalue_extraction(True)\n",
+    "sim.set_time_step(0.0001)\n",
+    "sim.set_final_time(0.1)\n",
+    "sim.add_logger(logger)\n",
+    "sim.run()\n",
+    "\n",
+    "# Read log file\n",
+    "eigenvalues_timeseries = rt.read_timeseries_dpsim(log_dir + '/eigenvalues.csv')\n",
+    "time = eigenvalues_timeseries['eigenvalues_0'].time\n",
+    "eigenvalue1_real = eigenvalues_timeseries['eigenvalues_0'].values.real\n",
+    "eigenvalue1_imag = eigenvalues_timeseries['eigenvalues_0'].values.imag\n",
+    "eigenvalue2_real = eigenvalues_timeseries['eigenvalues_1'].values.real\n",
+    "eigenvalue2_imag = eigenvalues_timeseries['eigenvalues_1'].values.imag\n",
+    "\n",
+    "simulation_timeseries = rt.read_timeseries_dpsim(log_dir + '/' + name + '.csv')\n",
+    "simulation_timeseries_time_domain = ts.frequency_shift_list(simulation_timeseries, 50)\n",
+    "print(simulation_timeseries_time_domain)\n",
+    "current_timeseries = simulation_timeseries_time_domain['current_shift']\n",
+    "voltage_timeseries = simulation_timeseries_time_domain['voltage_shift']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.title('Current')\n",
+    "plt.plot(current_timeseries.time, current_timeseries.values)\n",
+    "plt.show()\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.title('Voltage across resistor')\n",
+    "plt.plot(voltage_timeseries.time, voltage_timeseries.values)\n",
+    "plt.show()\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.title('First eigenvalue')\n",
+    "plt.plot(time, eigenvalue1_real, label='Real part')\n",
+    "plt.plot(time, eigenvalue1_imag, label='Imaginary part')\n",
+    "plt.legend()\n",
+    "plt.show()\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.title('Second eigenvalue')\n",
+    "plt.plot(time, eigenvalue2_real, label='Real part')\n",
+    "plt.plot(time, eigenvalue2_imag, label='Imaginary part')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}