From f20ad9f9c6de16bf2281872ea76c60b806187034 Mon Sep 17 00:00:00 2001 From: pipeacosta Date: Thu, 25 Jul 2024 09:01:03 +0200 Subject: [PATCH] Improvements to Kundur co-simulation Notebook Signed-off-by: pipeacosta --- ..._dyn_EMT_cosim_attributes_fault_1e-4.ipynb | 240 ++++++++---------- 1 file changed, 102 insertions(+), 138 deletions(-) diff --git a/examples/Notebooks/MatDyn/paper/Kundur2Areas_dyn_EMT_cosim_attributes_fault_1e-4.ipynb b/examples/Notebooks/MatDyn/paper/Kundur2Areas_dyn_EMT_cosim_attributes_fault_1e-4.ipynb index e8e554a94..48dc51176 100644 --- a/examples/Notebooks/MatDyn/paper/Kundur2Areas_dyn_EMT_cosim_attributes_fault_1e-4.ipynb +++ b/examples/Notebooks/MatDyn/paper/Kundur2Areas_dyn_EMT_cosim_attributes_fault_1e-4.ipynb @@ -81,16 +81,16 @@ "source": [ "time_step = 50e-6\n", "t_f = 1e-3\n", - "# start_fault_time = 0.5\n", - "# fault_clearing_time = 0.6\n", + "start_fault_time = 0.5\n", + "fault_clearing_time = 0.6\n", "\n", - "# switch_closed = 1 # fault resistance\n", - "# node_fault = \"N1\"\n", + "switch_closed = 1 # fault resistance\n", + "node_fault = \"N1\"\n", "\n", "# cosimulation parameters\n", "H = 0.2e-3\n", "interp = 'linear'\n", - "k_map = {'zoh': 0, 'linear': 1}\n", + "k_map = {'none': 0, 'zoh': 0, 'linear': 1}\n", "k = k_map[interp]\n", "\n", "cosim_config = {\n", @@ -186,7 +186,7 @@ " with_avr=False, with_tg=False, with_pss=False)\n", "\n", " ### Extend topology with 3ps fault\n", - " # sw = mpc_reader_dyn.add_3ph_faut(node_fault, switch_closed=switch_closed, switch_open=1e18)\n", + " sw = mpc_reader_dyn.add_3ph_faut(node_fault, switch_closed=switch_closed, switch_open=1e18)\n", " \n", " # create dpsim topology\n", " mpc_reader_dyn.create_dpsim_topology()\n", @@ -222,10 +222,10 @@ " sim.do_system_matrix_recomputation(True)\n", "\n", " # add events\n", - " # sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", - " # sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", - " # sim.add_event(sw_event_1)\n", - " # sim.add_event(sw_event_2)\n", + " sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", + " sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", + " sim.add_event(sw_event_1)\n", + " sim.add_event(sw_event_2)\n", " \n", " sim.run()\n", "\n", @@ -291,7 +291,7 @@ " with_avr=False, with_tg=False, with_pss=False)\n", "\n", "### Extend topology with 3ps fault\n", - "# sw = mpc_reader_dyn.add_3ph_faut(node_fault, switch_closed=switch_closed, switch_open=1e18)\n", + "sw = mpc_reader_dyn.add_3ph_faut(node_fault, switch_closed=switch_closed, switch_open=1e18)\n", "\n", "# create dpsim topology\n", "mpc_reader_dyn.create_dpsim_topology()\n", @@ -345,23 +345,16 @@ "# subdir_name_s1 = interp + '_' + str(H) + '/'\n", "dpsimpy.Logger.set_log_dir('logs/' + sim_name_dyn_s2)\n", "\n", - "init_current_line = sim_pf.get_idobj_attr(\"line5_8-9\", 'current_vector').get()[1][0] + sim_pf.get_idobj_attr(\"line6_8-9\", 'current_vector').get()[1][0]\n", - "init_current_parallel = sim.get_idobj_attr(\"line5_8-9\", 'ParallelCurrentNode1').get()[0][0] + sim.get_idobj_attr(\"line6_8-9\", 'ParallelCurrentNode1').get()[0][0]\n", + "# init_current_line = sim_pf.get_idobj_attr(\"line5_8-9\", 'current_vector').get()[1][0] + sim_pf.get_idobj_attr(\"line6_8-9\", 'current_vector').get()[1][0]\n", + "# init_current_parallel = sim.get_idobj_attr(\"line5_8-9\", 'ParallelCurrentNode1').get()[0][0] + sim.get_idobj_attr(\"line6_8-9\", 'ParallelCurrentNode1').get()[0][0]\n", "\n", "# set initial reference current of current source\n", - "init_current = init_current_line\n", - "\n", - "# print(sim_pf.get_idobj_attr(\"line5_8-9\", 'current_vector').get()[1][0] + sim_pf.get_idobj_attr(\"line6_8-9\", 'current_vector').get()[1][0])\n", + "# init_current = init_current_line\n", + "# init_current = dpsimpy.Math.single_phase_variable_to_three_phase(init_current * np.sqrt(3))\n", "# print(init_current)\n", + "# sys_topo_2.component(\"CS_N8\").set_parameters(init_current, 0)\n", "\n", - "#init_current = sim_pf.get_idobj_attr(\"line5_8-9\", 'current_vector').get()[1][0] + sim_pf.get_idobj_attr(\"line6_8-9\", 'current_vector').get()[1][0]\n", - "init_current = dpsimpy.Math.single_phase_variable_to_three_phase(init_current * np.sqrt(3))\n", - "# print(init_current)\n", - "sys_topo_2.component(\"CS_N8\").set_parameters(init_current, 0)\n", - "\n", - "# set parameters of current source\n", - "# init_voltage = sim_pf.get_idobj_attr(\"N8\", 'v').get()[0]\n", - "# sys_topo_2.component(\"CS_N8\").set_intf_voltage(np.real(init_voltage))\n", + "sys_topo_2.component(\"CS_N8\").set_parameters(dpsimpy.Math.single_phase_variable_to_three_phase(complex(0,0)), 0)\n", "\n", "# log results\n", "logger = dpsimpy.Logger(sim_name_dyn_s2)\n", @@ -388,16 +381,15 @@ "sim2.add_logger(logger)\n", "sim2.do_system_matrix_recomputation(True)\n", "\n", - "# if node_fault in cosim_config[\"nodes\"][1]:\n", - "# sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", - "# sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", - "# sim2.add_event(sw_event_1)\n", - "# sim2.add_event(sw_event_2)\n", + "if node_fault in cosim_config[\"nodes\"][1]:\n", + " sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", + " sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", + " sim2.add_event(sw_event_1)\n", + " sim2.add_event(sw_event_2)\n", " \n", "sim2.start()\n", - "y_2_0 = sim2.get_idobj_attr(\"CS_N8\", \"v_intf\").get()\n", - "print(y_2_0)\n", - "print(init_current_parallel)" + "y_2_0 = sim2.get_idobj_attr(\"CS_N8\", \"v_intf\").derive_coeff(0,0).get()\n", + "print(y_2_0)" ] }, { @@ -420,28 +412,18 @@ "dpsimpy.Logger.set_log_dir('logs/' + sim_name_dyn_s1)\n", "\n", "# set parameters of voltage source\n", - "# init_voltage = sim_pf.get_idobj_attr(\"N8\", 'v').get()[0]\n", - "# sys_topo_1.component(\"VS_N8\").set_parameters(dpsimpy.Math.single_phase_variable_to_three_phase(init_voltage), 0)\n", - "sys_topo_1.component(\"VS_N8\").set_parameters([complex(elem, 0) for elem in y_2_0 * dpsimpy.PEAK1PH_TO_RMS3PH], 0)\n", + "u_1_0 = y_2_0\n", + "sys_topo_1.component(\"VS_N8\").set_parameters(dpsimpy.Math.single_phase_variable_to_three_phase(u_1_0), 0)\n", "\n", "init_current_line = -(sim_pf.get_idobj_attr(\"line3_7-8\", 'current_vector').get()[0][0] + sim_pf.get_idobj_attr(\"line4_7-8\", 'current_vector').get()[0][0])\n", "init_current_parallel = sim.get_idobj_attr(\"line3_7-8\", 'ParallelCurrentNode0').get()[0][0] + sim.get_idobj_attr(\"line4_7-8\", 'ParallelCurrentNode0').get()[0][0]\n", "\n", "# get initial current of voltage source\n", - "init_current = init_current_line - 0*init_current_parallel\n", - " \n", - "# print(-(sim_pf.get_idobj_attr(\"line3_7-8\", 'current_vector').get()[0][0] + sim_pf.get_idobj_attr(\"line4_7-8\", 'current_vector').get()[0][0]))\n", - "# print(init_current)\n", - "\n", + "init_current = init_current_line\n", "init_current = dpsimpy.Math.single_phase_variable_to_three_phase(init_current * dpsimpy.RMS3PH_TO_PEAK1PH * np.sqrt(3))\n", - "# print(init_current)\n", - "sys_topo_1.component(\"VS_N8\").set_intf_current(np.real(init_current))\n", + "print(init_current)\n", "\n", - "# init_current = sim_pf.get_idobj_attr(\"line3_7-8\", 'current_vector').get()[0][0] + sim_pf.get_idobj_attr(\"line4_7-8\", 'current_vector').get()[0][0]\n", - "# init_current = dpsimpy.Math.single_phase_variable_to_three_phase(init_current)\n", - "# sys_topo_1.component(\"VS_N8\").set_intf_current(-np.real(init_current) * np.sqrt(3) * dpsimpy.RMS3PH_TO_PEAK1PH)\n", - "# print(init_current)\n", - "# sys_topo_1.component(\"VS_N8\").set_parameters(init_current, 0)\n", + "sys_topo_1.component(\"VS_N8\").set_intf_current(np.real(init_current))\n", "\n", "# log results\n", "logger = dpsimpy.Logger(sim_name_dyn_s1)\n", @@ -468,18 +450,15 @@ "sim1.do_system_matrix_recomputation(True)\n", "\n", "# add events\n", - "# if node_fault in cosim_config[\"nodes\"][0]:\n", - "# sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", - "# sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", - "# sim1.add_event(sw_event_1)\n", - "# sim1.add_event(sw_event_2)\n", + "if node_fault in cosim_config[\"nodes\"][0]:\n", + " sw_event_1 = dpsimpy.event.SwitchEvent(start_fault_time, sw, True)\n", + " sw_event_2 = dpsimpy.event.SwitchEvent(fault_clearing_time, sw, False)\n", + " sim1.add_event(sw_event_1)\n", + " sim1.add_event(sw_event_2)\n", "\n", "sim1.start()\n", - "sim1.get_idobj_attr(\"VS_N8\", \"V_ref\").set([complex(elem, 0) for elem in y_2_0 * dpsimpy.PEAK1PH_TO_RMS3PH])\n", - "y_1_0 = np.real(init_current)\n", - "print(init_current)\n", - "print(init_current_parallel)\n", - "print(sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get() * dpsimpy.PEAK1PH_TO_RMS3PH)" + "\n", + "y_1_0 = sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get()" ] }, { @@ -497,13 +476,11 @@ }, "outputs": [], "source": [ - "NO_EXTRAP = True\n", - "\n", "t_0 = 0\n", "t_k_v = []\n", - "j_v = []\n", "t_k = 0.0\n", - "y_1 = sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get()\n", + "\n", + "y_1 = y_1_0\n", "\n", "print(y_1)\n", "\n", @@ -517,10 +494,9 @@ "t_m = np.around(np.linspace(t_0, t_f, N + 1), 16)\n", "\n", "# We have to assume the trajectory of y_2 extending its initial value, since we have no prior information\n", - "# y_1_m_prev = np.tile(y_1_0, np.max([k, m]))\n", + "y_1_m_prev = np.tile(y_1_0, np.min([k+1, m]))\n", + "\n", "u_2_m_v = np.zeros((len(y_1_0), n))\n", - "# m = np.min([k+1, m+1])\n", - "y_1_m_prev = np.tile(y_1_0, np.min([k+1, m+1]))\n", "\n", "if interp == 'none':\n", " while t_k <= t_f:\n", @@ -534,63 +510,58 @@ " t_k = sim1.next()\n", " y_1 = sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get()\n", "\n", - "for i in range(0, N):\n", - " y_1_prev = y_1_m_prev[:,-1]\n", - " t_m_i = t[m*i: m*(i+1)]\n", - "\n", - " if interp == 'zoh':\n", - " u_2_m = np.tile(y_1_prev, (m,1)).T\n", - " elif interp == 'linear':\n", - " t_poly = np.array([i*H-H, i*H])\n", - " f_u_2 = np.polynomial.polynomial.polyfit(t_poly, y_1_m_prev[:,-2:].T, 1)\n", - " u_2_m = np.polynomial.polynomial.polyval(t_m_i, f_u_2)\n", - "\n", - " # u_2_m = np.zeros((len(y_1_0), m))\n", - " # for ind_2 in range(0, len(y_1_0)):\n", - " # f_u_2 = np.poly1d(np.polyfit([i*H-H, i*H], y_1_m_prev[ind_2, -2:], 1))\n", - " # u_2_m[ind_2,:] = f_u_2(t_m_i)\n", - "\n", - " u_2_m_v[:,m*i:m*(i+1)] = [np.real(u_2_m_elem) for u_2_m_elem in u_2_m]\n", - " y_1_m = np.zeros((len(y_1_0), m))\n", - " y_2_m = np.zeros((len(y_1_0), m))\n", - "\n", - " # Switch to S_2\n", - " for j in range(0, m):\n", - " u_2 = u_2_m[:,j]\n", - " \n", - " # print(\"u_2({}) = {}\".format(i, u_2))\n", - " sim2.get_idobj_attr(\"CS_N8\", \"I_ref\").set([complex(elem, 0) for elem in u_2 * dpsimpy.PEAK1PH_TO_RMS3PH])\n", - "\n", - " t_k_2 = sim2.next()\n", - " y_2 = sim2.get_idobj_attr(\"CS_N8\", \"v_intf\").get()\n", - " y_2_m[:,j] = y_2.flatten()\n", - "\n", - " # Switch to S_1\n", - " u_1_m = y_2_m \n", - " for j in range(0, m):\n", - " u_1 = u_1_m[:,j]\n", - " \n", - " sim1.get_idobj_attr(\"VS_N8\", \"V_ref\").set([complex(elem, 0) for elem in u_1 * dpsimpy.PEAK1PH_TO_RMS3PH])\n", - " t_k_1 = sim1.next()\n", - " y_1 = sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get()\n", - " # print(y_1)\n", - "\n", - " y_1_m[:,j] = y_1.flatten()\n", - " t_k_v.append(t_k_1)\n", - "\n", - " # Option 1: Take values at macrosteps as samples for extrapolation\n", - " y_1_m_prev = np.hstack((y_1_m_prev, y_1_m[:,-1][np.newaxis].T))\n", - "\n", - " # Option 2: Take values at microsteps as samples for extrapolation\n", - " # y_1_m_prev = y_1_m\n", - "\n", - " # j_v.append(j)\n", - "\n", - " if t_k_1 > t_f:\n", - " print(t_k_1)\n", - " print(\"i: \" + str(i) + \", expected: \" + str(N-1))\n", - " print(\"j: \" + str(j) + \", expected: \" + str(m-1))\n", - " break\n", + "else:\n", + "\n", + " for i in range(0, N):\n", + " y_1_prev = y_1_m_prev[:,-1]\n", + " t_m_i = t[m*i: m*(i+1)]\n", + "\n", + " if interp == 'zoh':\n", + " u_2_m = np.tile(y_1_prev, (m,1)).T\n", + " elif interp == 'linear':\n", + " t_poly = np.array([i*H-H, i*H])\n", + " f_u_2 = np.polynomial.polynomial.polyfit(t_poly, y_1_m_prev[:,-2:].T, 1)\n", + " u_2_m = np.polynomial.polynomial.polyval(t_m_i, f_u_2)\n", + "\n", + " u_2_m_v[:,m*i:m*(i+1)] = [np.real(u_2_m_elem) for u_2_m_elem in u_2_m]\n", + " y_1_m = np.zeros((len(y_1_0), m))\n", + " y_2_m = np.zeros((len(y_1_0), m))\n", + "\n", + " # Switch to S_2\n", + " for j in range(0, m):\n", + " u_2 = u_2_m[:,j]\n", + " \n", + " # print(\"u_2({}) = {}\".format(i, u_2))\n", + " sim2.get_idobj_attr(\"CS_N8\", \"I_ref\").set([complex(elem, 0) for elem in u_2 * dpsimpy.PEAK1PH_TO_RMS3PH])\n", + "\n", + " t_k_2 = sim2.next()\n", + " y_2 = sim2.get_idobj_attr(\"CS_N8\", \"v_intf\").get()\n", + " y_2_m[:,j] = y_2.flatten()\n", + "\n", + " # Switch to S_1\n", + " u_1_m = y_2_m \n", + " for j in range(0, m):\n", + " u_1 = u_1_m[:,j]\n", + " \n", + " sim1.get_idobj_attr(\"VS_N8\", \"V_ref\").set([complex(elem, 0) for elem in u_1 * dpsimpy.PEAK1PH_TO_RMS3PH])\n", + " t_k_1 = sim1.next()\n", + " y_1 = sim1.get_idobj_attr(\"VS_N8\", \"i_intf\").get()\n", + " # print(y_1)\n", + "\n", + " y_1_m[:,j] = y_1.flatten()\n", + " t_k_v.append(t_k_1)\n", + "\n", + " # Option 1: Take values at macrosteps as samples for extrapolation\n", + " y_1_m_prev = np.hstack((y_1_m_prev, y_1_m[:,-1][np.newaxis].T))\n", + "\n", + " # Option 2: Take values at microsteps as samples for extrapolation\n", + " # y_1_m_prev = y_1_m\n", + "\n", + " if t_k_1 > t_f:\n", + " print(t_k_1)\n", + " print(\"i: \" + str(i) + \", expected: \" + str(N-1))\n", + " print(\"j: \" + str(j) + \", expected: \" + str(m-1))\n", + " break\n", "\n", "sim1.stop()\n", "sim2.stop()" @@ -763,11 +734,11 @@ }, "outputs": [], "source": [ - "varname_dpsim = 'N8.V'\n", - "# nominal_voltage = 230000\n", - "nominal_voltage = 1\n", - "plot_node_volt_abs(varname_dpsim, ts_dpsim_fault_emt, ts_dpsim_fault_s1_emt, ts_dpsim_fault_s2_emt, nominal_voltage, timestep_common)\n", - "print('H=' + str(H))" + "# varname_dpsim = 'N8.V'\n", + "# # nominal_voltage = 230000\n", + "# nominal_voltage = 1\n", + "# plot_node_volt_abs(varname_dpsim, ts_dpsim_fault_emt, ts_dpsim_fault_s1_emt, ts_dpsim_fault_s2_emt, nominal_voltage, timestep_common)\n", + "# print('H=' + str(H))" ] }, { @@ -829,9 +800,9 @@ "plt.plot(ts_mono, dpsim_emt_i_intf_values, label='EMT - Monolithic')\n", "plt.plot(ts_s1, dpsim_emt_i_intf_values_s1, '--', label='EMT - S1')\n", "# plt.plot(time_H, u_2_m_v_H, 'o', label='Macro-step')\n", - "plt.plot(ts_s2, dpsim_emt_i_intf_values_s2, label='EMT - S2')\n", - "plt.plot(ts_s2, dpsim_emt_i_ref_values_s2, '-.', label='EMT - S2, ref')\n", - "plt.plot(ts_mono, u_2_m_v[0,:], '.-', label='Linear prediction')\n", + "plt.plot(ts_s2, dpsim_emt_i_intf_values_s2, '-.', label='EMT - S2')\n", + "# plt.plot(ts_s2, dpsim_emt_i_ref_values_s2, '-.', label='EMT - S2, ref')\n", + "plt.plot(ts_mono, u_2_m_v[0,:], ':', label='Linear prediction')\n", "# plt.plot(time_H, dpsim_emt_values_s2_H, 'o', label='Macro-step')\n", "plt.legend(loc='lower right')\n", "plt.xlabel('time (s)')\n", @@ -854,24 +825,17 @@ "plt.plot(ts_mono, dpsim_emt_v_intf_values, label='EMT - Monolithic')\n", "plt.plot(ts_s1, dpsim_emt_v_intf_values_s1, '--', label='EMT - S1')\n", "plt.plot(ts_s2, dpsim_emt_v_intf_values_s2, '-.', label='EMT - S2')\n", - "plt.plot(ts_s1, dpsim_emt_v_ref_values_s1, '--', label='EMT - S1, ref')\n", + "# plt.plot(ts_s1, dpsim_emt_v_ref_values_s1, '--', label='EMT - S1, ref')\n", "# plt.plot(time_H, dpsim_emt_values_s2_H, 'o', label='Macro-step')\n", "plt.legend(loc='lower right')\n", "plt.xlabel('time (s)')\n", "plt.grid()\n", - "plt.ylim([168000, 182000])\n", + "# plt.ylim([168000, 182000])\n", "# plt.xlim([0, 0.6e-3])\n", "plt.xlabel(\"time (s)\")\n", "plt.ylabel(\"Interface voltage\")" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "code", "execution_count": null,