Skip to content

Commit

Permalink
Initialization improvements to attirbutes-based co-simulation example.
Browse files Browse the repository at this point in the history 
Signed-off-by: pipeacosta <[email protected]>
  • Loading branch information
@pipeacosta
pipeacosta committed Dec 9, 2024
1 parent 9ca93f3 commit 366c494
Show file tree
Hide file tree
Showing 2 changed files with 143 additions and 38 deletions.
155 changes: 122 additions & 33 deletions examples/Notebooks/Attributes/cosimulation-attributes.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -209,31 +212,77 @@
"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",
" print(m)\n",
"\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",
"\n",
" # 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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -529,9 +623,7 @@
"plt.legend()\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
" \n",
" "
"plt.show()"
]
},
{
Expand All @@ -550,11 +642,8 @@
}
],
"metadata": {
"interpreter": {
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "usr",
"language": "python",
"name": "python3"
},
Expand Down
26 changes: 21 additions & 5 deletions examples/Python/Attributes/emt-cosim-attributes_polynomial.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -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':
Expand All @@ -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

Expand Down

0 comments on commit 366c494

Please sign in to comment.