Skip to content

Commit

Permalink
Improvements to Kundur co-simulation Notebook
Browse files Browse the repository at this point in the history 
Signed-off-by: pipeacosta <pipeacosta@gmail.com>
  • Loading branch information
@pipeacosta
pipeacosta committed Jul 25, 2024
1 parent 5d6af8e commit f20ad9f
Showing 1 changed file with 102 additions and 138 deletions.
240 changes: 102 additions & 138 deletions examples/Notebooks/MatDyn/paper/Kundur2Areas_dyn_EMT_cosim_attributes_fault_1e-4.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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)"
]
},
{
Expand All @@ -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",
Expand All @@ -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()"
]
},
{
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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()"
Expand Down Expand Up @@ -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))"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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,
Expand Down

0 comments on commit f20ad9f

Please sign in to comment.