diff --git a/examples/Notebooks/Attributes/cosimulation-attributes.ipynb b/examples/Notebooks/Attributes/cosimulation-attributes.ipynb index ca401bfa7..bc0eba4e8 100644 --- a/examples/Notebooks/Attributes/cosimulation-attributes.ipynb +++ b/examples/Notebooks/Attributes/cosimulation-attributes.ipynb @@ -170,7 +170,10 @@ "outputs": [], "source": [ "time_step = 50e-6\n", - "final_time = 1.0\n", + "final_time = 1.2\n", + "\n", + "# time_step = 0.001\n", + "# final_time = 1.0\n", "\n", "sim = dpsimpy.Simulation(\"EMTCosim\", loglevel=dpsimpy.LogLevel.debug)\n", "sim.set_domain(dpsimpy.Domain.EMT)\n", @@ -209,14 +212,60 @@ "import subprocess\n", "import numpy as np\n", "\n", - "# H_v = np.array([1e-4, 2e-4, 4e-4, 0.001, 0.003, 0.01])\n", - "# H_v_legends = ['1e-4', '2e-4', '4e-4', '1e-3', '3e-3', '1e-2']\n", - "H_v = [0.2e-3]\n", - "H_v_legends = ['0.2e-3']\n", + "H_v = np.array([1e-4, 2e-4, 4e-4, 0.001, 0.003, 0.01])\n", + "H_v_legends = ['1e-4', '2e-4', '4e-4', '1e-3', '3e-3', '1e-2']\n", + "# H_v = [0.2e-3]\n", + "# H_v_legends = ['0.2e-3']\n", + "\n", + "# H_v = [1e-4]\n", + "# H_v_legends = ['1e-4']\n", + "\n", + "# H_v = [0.003]\n", + "# H_v_legends = ['0.003']\n", "\n", - "# H_v = [4e-4]\n", - "# H_v_legends = ['4e-4']\n", + "# Consistent initial conditions for linear extrapolation and num_vs = 0\n", + "# y_1_prev_H = [30.1080, 30.216, 30.432, 31.08, 33.24, 40.8]\n", + "y_1_prev_H = [0, 0, 0, 0, 0, 0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# This here is for evaluating the trapezoidal method error against Matlab\n", + "# i=0\n", + "# for time_step_mono in H_v:\n", + "# sim = dpsimpy.Simulation(\"EMTCosim_\" + H_v_legends[i], loglevel=dpsimpy.LogLevel.debug)\n", + "# sim.set_domain(dpsimpy.Domain.EMT)\n", + "# sim.set_system(sys)\n", + "# sim.set_time_step(time_step_mono)\n", + "# sim.set_final_time(final_time)\n", + "\n", + "# log = dpsimpy.Logger(\"EMTCosim_\" + H_v_legends[i])\n", + "\n", + "# if num_vs > 0:\n", + "# log.log_attribute(\"v0\", \"v\", sys.node(\"n0\"))\n", + "# if num_vs > 1:\n", + "# log.log_attribute(\"v3\", \"v\", sys.node(\"n3\"))\n", + "# log.log_attribute(\"v1\", \"v\", sys.node(\"n1\"))\n", + "# log.log_attribute(\"v2\", \"v\", sys.node(\"n2\"))\n", + "# log.log_attribute(\"r_line.I\", \"i_intf\", sys.component(\"r_line\"))\n", + "\n", + "# sim.add_logger(log)\n", + " \n", + "# sim.run()\n", "\n", + "# i+=1\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "i=0\n", "for H in H_v:\n", " m = int(H/time_step)\n", @@ -224,7 +273,7 @@ "\n", " process1 = subprocess.Popen([\"python3\", \"../../Python/Attributes/emt-cosim-attributes_polynomial.py\", \"-t\", str(time_step), \"-e\", str(final_time), \"-H\", str(H), \"-i\", \"zoh\", \"-p\", \"EMTCosimAttributes_zoh_\" + H_v_legends[i], \"--num-vs\", str(num_vs)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", " output1, error1 = process1.communicate()\n", - " process2 = subprocess.Popen([\"python3\", \"../../Python/Attributes/emt-cosim-attributes_polynomial.py\", \"-t\", str(time_step), \"-e\", str(final_time), \"-H\", str(H), \"-i\", \"linear\", \"-p\", \"EMTCosimAttributes_linear_\" + H_v_legends[i], \"--num-vs\", str(num_vs)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", + " process2 = subprocess.Popen([\"python3\", \"../../Python/Attributes/emt-cosim-attributes_polynomial.py\", \"-t\", str(time_step), \"-e\", str(final_time), \"-H\", str(H), \"-i\", \"linear\", \"--y1-prev\", str(y_1_prev_H[i]), \"-p\", \"EMTCosimAttributes_linear_\" + H_v_legends[i], \"--num-vs\", str(num_vs)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", " output2, error2 = process2.communicate()\n", " process3 = subprocess.Popen([\"python3\", \"../../Python/Attributes/emt-cosim-attributes_polynomial.py\", \"-t\", str(time_step), \"-e\", str(final_time), \"-H\", str(H), \"-i\", \"none\", \"-p\", \"EMTCosimAttributes_no_extrap_\" + H_v_legends[i], \"--num-vs\", str(num_vs)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", " output3, error3 = process3.communicate()\n", @@ -232,8 +281,8 @@ " # Print output and error for debugging\n", " print(\"Output 1:\", output1)\n", " print(\"Error 1:\", error1)\n", - " print(\"Output 2:\", output2)\n", - " print(\"Error 2:\", error2)\n", + " # print(\"Output 2:\", output2)\n", + " # print(\"Error 2:\", error2)\n", " print(\"Output 3:\", output3)\n", " print(\"Error 3:\", error3)\n", "\n", @@ -268,7 +317,7 @@ "import villas.dataprocessing.plottools as pt\n", "import villas.dataprocessing.readtools as rt\n", "import villas.dataprocessing.timeseries as ts\n", - "# %matplotlib widget\n", + "%matplotlib widget\n", "\n", "results_emt = rt.read_timeseries_dpsim('logs/EMTCosim.csv')\n", "\n", @@ -416,8 +465,8 @@ "# pt.set_timeseries_labels(results_emt['v0'], 'v_0 - Monolithic')\n", "# pt.plot_timeseries('Co-simulation results - v_src - no extrap.', results_emt['v0'])\n", "\n", - "pt.set_timeseries_labels(results_emt['v1'], 'v_1 - Monolithic')\n", - "pt.plot_timeseries('Co-simulation results - v_intf - no extrap.', results_emt['v1'])\n", + "# pt.set_timeseries_labels(results_emt['v1'], 'v_1 - Monolithic')\n", + "# pt.plot_timeseries('Co-simulation results - v_intf - no extrap.', results_emt['v1'])\n", "pt.set_timeseries_labels(results_emt['v2'], 'v_2 - Monolithic')\n", "pt.plot_timeseries('Co-simulation results - v_intf - no extrap.', results_emt['v2'])\n", "plt.xlabel('t [s]')\n", @@ -472,19 +521,24 @@ "global_error_zoh = []\n", "global_error_linear = []\n", "\n", - "ts_1_m = results_emt['v1']\n", - "ts_2_m = results_emt['v2']\n", + "ts_v1_m = results_emt['v1']\n", + "ts_v2_m = results_emt['v2']\n", + "ts_i_intf_m = results_emt['r_line.I']\n", + "\n", + "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(7, 6))\n", "\n", "for i in range(0, len(H_v_legends)): \n", " H = H_v[i]\n", " \n", - " m = int(H/time_step)\n", + " m = int(round(H/time_step))\n", "\n", - " ts_1_zoh = results_emt_attributes_zoh_S1[H_v_legends[i]]['v_1']\n", - " ts_2_zoh = results_emt_attributes_zoh_S1[H_v_legends[i]]['v_2']\n", + " ts_v1_zoh = results_emt_attributes_zoh_S1[H_v_legends[i]]['v_1']\n", + " ts_v2_zoh = results_emt_attributes_zoh_S1[H_v_legends[i]]['v_2']\n", + " ts_i_intf_zoh = results_emt_attributes_zoh_S1[H_v_legends[i]]['i_intf']\n", " \n", - " ts_1_linear = results_emt_attributes_linear_S1[H_v_legends[i]]['v_1']\n", - " ts_2_linear = results_emt_attributes_linear_S1[H_v_legends[i]]['v_2']\n", + " ts_v1_linear = results_emt_attributes_linear_S1[H_v_legends[i]]['v_1']\n", + " ts_v2_linear = results_emt_attributes_linear_S1[H_v_legends[i]]['v_2']\n", + " ts_i_intf_linear = results_emt_attributes_linear_S1[H_v_legends[i]]['i_intf']\n", "\n", " # if num_vs == 0:\n", " # ts_1_m = results_emt[1]\n", @@ -505,15 +559,55 @@ " # ts_1_linear = results_emt_attributes_linear_S1[H_v_legends[i]][2]\n", " # ts_2_linear = results_emt_attributes_linear_S1[H_v_legends[i]][3]\n", " \n", - " v_a = np.array([ts_1_m.values[::m], ts_2_m.values[::m]])\n", - " v_zoh = np.array([ts_1_zoh.values[::m], ts_2_zoh.values[::m]])\n", - " v_linear = np.array([ts_1_linear.values[::m], ts_2_linear.values[::m]])\n", + " xy_a = np.array([ts_v1_m.values[::m], ts_v2_m.values[::m], ts_i_intf_m.values[::m], ts_v2_m.values[::m]])\n", + " xy_zoh = np.array([ts_v1_zoh.values[::m], ts_v2_zoh.values[::m], ts_i_intf_zoh.values[::m], ts_v2_zoh.values[::m]])\n", + " xy_linear = np.array([ts_v1_linear.values[::m], ts_v2_linear.values[::m], ts_i_intf_linear.values[::m], ts_v2_linear.values[::m]])\n", " \n", - " gloal_error_zoh_t = np.max(np.linalg.norm(v_a - v_zoh, axis=0))\n", + " # global_error_v1_zoh_t = np.abs(ts_v1_m.values[::m] - ts_v1_zoh.values[::m], axis=0)\n", + " # global_error_v2_zoh_t = np.abs(ts_v2_m.values[::m] - ts_v2_zoh.values[::m], axis=0)\n", + " # global_error_i_intf_zoh_t = np.abs(ts_i_intf_m.values[::m] - ts_i_intf_zoh.values[::m], axis=0)\n", + "\n", + " global_error_v1_linear_t = np.abs(ts_v1_m.values[::m] - ts_v1_linear.values[::m])\n", + " global_error_v2_linear_t = np.abs(ts_v2_m.values[::m] - ts_v2_linear.values[::m])\n", + " global_error_i_intf_linear_t = np.abs(ts_i_intf_m.values[::m] - ts_i_intf_linear.values[::m])\n", + " \n", + " global_error_zoh_t = np.linalg.norm(xy_a - xy_zoh, axis=0)\n", + " global_error_linear_t = np.linalg.norm(xy_a - xy_linear, axis=0)\n", + "\n", + " # plt.plot(ts_1_m.time[::m], gloal_error_zoh_t, label='ZOH - ' + H_v_legends[i])\n", + " # plt.plot(ts_1_m.time[::m], gloal_error_linear_t, label='linear - ' + H_v_legends[i])\n", + "\n", + " axes[0].plot(ts_v1_m.time[::m], global_error_v1_linear_t, label='V1 - linear - ' + H_v_legends[i])\n", + " axes[0].plot(ts_v2_m.time[::m], global_error_v2_linear_t, label='V2 - linear - ' + H_v_legends[i])\n", + "\n", + " axes[1].plot(ts_i_intf_m.time[::m], global_error_i_intf_linear_t, label='i_intf - linear - ' + H_v_legends[i])\n", " \n", - " global_error_zoh.append(gloal_error_zoh_t)\n", - " global_error_linear.append(np.max(np.linalg.norm(v_a - v_linear, axis=0)))\n", + " global_error_zoh.append(np.max(global_error_zoh_t))\n", + " global_error_linear.append(np.max(np.max(global_error_linear_t)))\n", + "\n", + "axes[0].set_xlabel('t [s]')\n", + "axes[0].set_ylabel('Global error')\n", + "# plt.xlim((0, 0.4))\n", + "# plt.ylim((0.0008, 0.0012))\n", + "axes[0].grid()\n", + "axes[0].legend()\n", + "\n", + "axes[1].set_xlabel('t [s]')\n", + "axes[1].set_ylabel('Global error')\n", + "# plt.xlim((0, 0.4))\n", + "# plt.ylim((0.0008, 0.0012))\n", + "axes[1].grid()\n", + "axes[1].legend()\n", " \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "plt.figure()\n", "plt.plot(np.log10(H_v), np.log10(global_error_zoh), 'b', label='ZOH')\n", "plt.plot(np.log10(H_v), np.log10(global_error_zoh), 'bo')\n", @@ -529,9 +623,7 @@ "plt.legend()\n", "\n", "plt.tight_layout()\n", - "plt.show()\n", - " \n", - " " + "plt.show()" ] }, { @@ -550,11 +642,8 @@ } ], "metadata": { - "interpreter": { - "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a" - }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "usr", "language": "python", "name": "python3" }, diff --git a/examples/Python/Attributes/emt-cosim-attributes_polynomial.py b/examples/Python/Attributes/emt-cosim-attributes_polynomial.py index 6a2f82ee1..144681929 100644 --- a/examples/Python/Attributes/emt-cosim-attributes_polynomial.py +++ b/examples/Python/Attributes/emt-cosim-attributes_polynomial.py @@ -141,6 +141,8 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): v_s_2.set_parameters(V_ref=complex(np.sqrt(2)/2, np.sqrt(2)/2), f_src=50) ecs = dpsimpy.emt.ph1.CurrentSource("i_intf", dpsimpy.LogLevel.info) + # ecs.set_parameters(30.0) + c_2 = dpsimpy.emt.ph1.Capacitor("c_2", dpsimpy.LogLevel.info) c_2.set_parameters(c_2_c) r_load = dpsimpy.emt.ph1.Resistor("r_load", dpsimpy.LogLevel.info) @@ -213,6 +215,7 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): parser.add_argument('-e', '--end', type=float, required=True) parser.add_argument('-H', '--macro-step', type=float, required=True) parser.add_argument('-i', '--interp', default='none') + parser.add_argument('--y1-prev', type=float, default=0) parser.add_argument('-p', '--prefix', default='') parser.add_argument('--num-vs', default=0) parser.add_argument('-d', '--debug', type=bool, default=False) @@ -223,10 +226,15 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): t_f = args.end H = args.macro_step interp = args.interp + y_1_prev = args.y1_prev prefix = args.prefix num_vs = int(args.num_vs) debug = args.debug + if interp == 'linear' and y_1_prev is None: + print('Error: Linear extrapolation requires initial previous value!') + exit(1) + if interp == 'none': NO_EXTRAP = True else: @@ -266,15 +274,20 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): y_1 = y_1_0 # We have to assume the trajectory of y_2 extending its initial value, since we have no prior information - # y_1_m_prev = np.tile(y_1_0, m) - y_1_m_prev = np.array([0.0, y_1_0]) + y_1_m_prev = np.tile(y_1_0, m) + # y_1_m_prev = np.array([y_1_prev, y_1_0]) + # y_1_m_prev = np.array([0.0, y_1_0]) if NO_EXTRAP: while t_k <= t_f: u_2 = y_1 sim2.get_idobj_attr("i_intf", "I_ref").set(complex(u_2,0)) - sim2.next() + + if t_k == 0.0: + sim2.start() + else: + sim2.next() y_2 = sim2.get_idobj_attr("i_intf", "v_intf").derive_coeff(0,0).get() u_1 = y_2 @@ -287,7 +300,7 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): for i in range(0, N): y_1_prev = y_1_m_prev[-1] - t_m_i = t[m*i : m*(i+1)] + t_m_i = t[m*(i)+1:m*(i+1)+1] # Extrapolation: Zero order hold if interp == 'zoh': @@ -314,7 +327,10 @@ def set_dpsim2(t_s, t_f, num_vs, logger_prefix): u_2_test = sim2.get_idobj_attr("i_intf", "I_ref").get() print("Input value in S2 after set: {:f}".format(u_2_test)) - sim2.next() + if i == 0 and j == 0: + t_k_2 = sim2.start() + else: + t_k_2 = sim2.next() y_2 = sim2.get_idobj_attr("i_intf", "v_intf").derive_coeff(0,0).get() y_2_m[j] = y_2