From 29477678d5a22c3cce797c690a8021f6553075a4 Mon Sep 17 00:00:00 2001 From: Marko Hinkkanen Date: Sun, 22 Jun 2025 12:16:00 +0300 Subject: [PATCH] Refactored plotting functions for better clarity and modularity. Enhanced plotting utilities for simulation results and flux maps (optional saving of figures, optional limits and ticks). Separated flux map plotting functions to an own module (_sm_plot_flux_maps.py). Fixed the rated voltage and the rated torque of the 5.6-kW PM-SyRM. --- .vscode/settings.json | 7 +- docs/source/conf.py | 2 +- .../plot_current_vector_pmsm_2kw_diode.py | 3 +- .../plot_current_vector_pmsyrm_thor_sat.py | 3 +- .../plot_flux_vector_ctrl_pmsyrm_thor_sat.py | 5 +- .../plot_flux_vector_pmsyrm_5kw_sat.py | 39 +- .../plot_flux_vector_syrm_7kw_sat.py | 5 +- .../vhz/plot_obs_vhz_pmsm_2kw_two_mass.py | 6 +- examples/drive/vhz/plot_vhz_im_2kw.py | 4 +- examples/drive/vhz/plot_vhz_im_2kw_lc.py | 6 +- .../grid/grid_following/plot_gfl_10kva.py | 4 +- .../grid_following/plot_gfl_dc_bus_10kva.py | 3 +- .../grid/grid_following/plot_gfl_lcl_10kva.py | 5 +- .../grid/grid_forming/plot_gfm_obs_13kva.py | 3 +- .../grid/grid_forming/plot_gfm_rfpsc_13kva.py | 7 +- motulator/common/model/_converter.py | 2 +- motulator/common/utils/__init__.py | 17 +- motulator/common/utils/_plotting.py | 124 +++++ motulator/common/utils/_utils.py | 54 -- motulator/drive/control/_im_current_vector.py | 2 +- motulator/drive/control/_im_flux_vector.py | 2 +- motulator/drive/control/_im_observers.py | 2 +- motulator/drive/control/_sm_observers.py | 2 +- motulator/drive/control/_sm_reference_gen.py | 2 +- motulator/drive/control/_sm_signal_inj.py | 2 +- motulator/drive/model/_lc_filter.py | 2 +- motulator/drive/model/_machine.py | 2 +- motulator/drive/model/_mechanics.py | 2 +- motulator/drive/utils/__init__.py | 14 +- motulator/drive/utils/_plots.py | 493 +++++++++++------- motulator/drive/utils/_sm_flux_maps.py | 283 +--------- .../drive/utils/_sm_plot_control_loci.py | 56 +- motulator/drive/utils/_sm_plot_flux_maps.py | 347 ++++++++++++ motulator/grid/control/_base.py | 2 +- motulator/grid/control/_gfl_current_vector.py | 2 +- motulator/grid/control/_gfm_observer.py | 2 +- motulator/grid/control/_gfm_psc.py | 2 +- motulator/grid/model/_ac_filter.py | 2 +- motulator/grid/model/_ac_source.py | 2 +- motulator/grid/utils/__init__.py | 9 +- motulator/grid/utils/_plots.py | 378 +++++++++----- pyproject.toml | 2 +- 42 files changed, 1194 insertions(+), 717 deletions(-) create mode 100644 motulator/common/utils/_plotting.py create mode 100644 motulator/drive/utils/_sm_plot_flux_maps.py diff --git a/.vscode/settings.json b/.vscode/settings.json index df5a7ac50..058451d56 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -24,6 +24,7 @@ "autoapi", "Awan", "axhline", + "axlim", "Baldor", "Blanco", "Böcker", @@ -36,6 +37,7 @@ "Capecchi", "Capponi", "Claypool", + "cmap", "Degner", "dejavuserif", "docstrings", @@ -154,11 +156,13 @@ "vartheta", "Verghese", "viridis", + "xdata", "xlabel", "xlim", "xtick", "xticklabels", "xticks", + "ydata", "Yepes", "ylabel", "ylabels", @@ -167,6 +171,7 @@ "yticks", "Zigliotto", "zlabel", - "zlim" + "zlim", + "zticks" ], } diff --git a/docs/source/conf.py b/docs/source/conf.py index 2294e6a48..a3610484b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -26,7 +26,7 @@ author = "Aalto Electric Drives" # The full version, including alpha/beta/rc tags -release = "0.6.1" +release = "0.6.2" # -- General configuration ------------------------------------------------------------- diff --git a/examples/drive/current_vector/plot_current_vector_pmsm_2kw_diode.py b/examples/drive/current_vector/plot_current_vector_pmsm_2kw_diode.py index c9691fcaf..9a6e7f8e8 100644 --- a/examples/drive/current_vector/plot_current_vector_pmsm_2kw_diode.py +++ b/examples/drive/current_vector/plot_current_vector_pmsm_2kw_diode.py @@ -55,4 +55,5 @@ # Plot also the stator voltage and currents as well as the DC-bus and grid-side # quantities. -utils.plot_extra(res, base, t_span=(0.8, 0.825)) +utils.plot_stator_waveforms(res, base, t_lims=(0.8, 0.825)) +utils.plot_dc_bus_waveforms(res, base, t_lims=(0.8, 0.825)) diff --git a/examples/drive/current_vector/plot_current_vector_pmsyrm_thor_sat.py b/examples/drive/current_vector/plot_current_vector_pmsyrm_thor_sat.py index fc5b1bb7f..21b30c8dc 100644 --- a/examples/drive/current_vector/plot_current_vector_pmsyrm_thor_sat.py +++ b/examples/drive/current_vector/plot_current_vector_pmsyrm_thor_sat.py @@ -33,7 +33,8 @@ # Get the path of the MATLAB file and load the FEM data p = Path(__file__).resolve().parent if "__file__" in globals() else Path.cwd() fem_flux_map = utils.import_syre_data(str(p / "THOR.mat")) -utils.plot_maps(fem_flux_map, base, x_lims=(-2, 2), y_lims=(-2, 2)) +utils.plot_map(fem_flux_map, "d", base, x_lims=(-2, 2), y_lims=(-2, 2)) +utils.plot_map(fem_flux_map, "q", base, x_lims=(-2, 2), y_lims=(-2, 2)) # %% # Configure the system model. diff --git a/examples/drive/flux_vector/plot_flux_vector_ctrl_pmsyrm_thor_sat.py b/examples/drive/flux_vector/plot_flux_vector_ctrl_pmsyrm_thor_sat.py index a307d79df..2c4d802de 100644 --- a/examples/drive/flux_vector/plot_flux_vector_ctrl_pmsyrm_thor_sat.py +++ b/examples/drive/flux_vector/plot_flux_vector_ctrl_pmsyrm_thor_sat.py @@ -37,10 +37,11 @@ fem_flux_map = utils.import_syre_data(str(p / "THOR.mat")) # %% -# Plot the maps in per-unit values +# Plot the maps in per-unit values. # sphinx_gallery_thumbnail_number = 3 -utils.plot_maps(fem_flux_map, base, x_lims=(-2, 2), y_lims=(-2, 2)) +utils.plot_map(fem_flux_map, "d", base, x_lims=(-2, 2), y_lims=(-2, 2)) +utils.plot_map(fem_flux_map, "q", base, x_lims=(-2, 2), y_lims=(-2, 2)) # %% # Two-dimensional presentation of flux maps. diff --git a/examples/drive/flux_vector/plot_flux_vector_pmsyrm_5kw_sat.py b/examples/drive/flux_vector/plot_flux_vector_pmsyrm_5kw_sat.py index 9fe57d0a4..ce7af1c16 100644 --- a/examples/drive/flux_vector/plot_flux_vector_pmsyrm_5kw_sat.py +++ b/examples/drive/flux_vector/plot_flux_vector_pmsyrm_5kw_sat.py @@ -1,8 +1,8 @@ """ -5.5-kW PM-SyRM, saturated +5.6-kW PM-SyRM, saturated ========================= -This example simulates sensorless stator-flux-vector control of a 5.5-kW PM-SyRM (Baldor +This example simulates sensorless stator-flux-vector control of a 5.6-kW PM-SyRM (Baldor ECS101M0H7EF4) drive. The machine model is parametrized using the flux map data, measured using the constant-speed test. The control system is parametrized using the algebraic saturation model from [#Lel2024]_, fitted to the measured data. This @@ -15,7 +15,6 @@ from pathlib import Path -import matplotlib.pyplot as plt import numpy as np from scipy.io import loadmat @@ -25,7 +24,7 @@ # %% # Compute base values based on the nominal values (just for figures). -nom = utils.NominalValues(U=370, I=8.8, f=60, P=5.5e3, tau=29.2) +nom = utils.NominalValues(U=460, I=8.8, f=60, P=5.6e3, tau=29.7) base = utils.BaseValues.from_nominal(nom, n_p=2) # %% @@ -45,7 +44,7 @@ # Plot the measured data # sphinx_gallery_thumbnail_number = 4 -utils.plot_flux_vs_current(meas_flux_map, base, lims=(-1.5, 1.5)) +utils.plot_flux_vs_current(meas_flux_map, base, x_lims=(-1.5, 1.5)) # %% # Create a saturation model, which will be used in the control system. @@ -72,14 +71,33 @@ # Generate the flux map using the saturation model est_curr_map = i_s_dq_fcn.as_magnetic_model( - d_range=np.linspace(-0.1, 1.3 * base.psi, 256), - q_range=np.linspace(-1.7 * base.psi, 1.7 * base.psi, 256), + d_range=np.linspace(-0.1, base.psi, 256), + q_range=np.linspace(-1.4 * base.psi, 1.4 * base.psi, 256), ) est_flux_map = est_curr_map.invert() # Plot the saturation model (surface) and the measured data (points) -utils.plot_maps( - est_flux_map, base, x_lims=(-1.8, 1.8), y_lims=(-2.15, 2.15), raw_data=meas_flux_map +utils.plot_map( + est_flux_map, + "d", + base, + x_lims=(-2, 2), + y_lims=(-2, 2), + z_lims=(0, 1), + x_ticks=[-2, -1, 0, 1, 2], + y_ticks=[-2, -1, 0, 1, 2], + raw_data=meas_flux_map, +) +utils.plot_map( + est_flux_map, + "q", + base, + x_lims=(-2, 2), + y_lims=(-2, 2), + z_lims=(-1.5, 1.5), + x_ticks=[-2, -1, 0, 1, 2], + y_ticks=[-2, -1, 0, 1, 2], + raw_data=meas_flux_map, ) # %% @@ -113,7 +131,6 @@ mc.plot_current_vs_torque(i_s_vals, base) mc.plot_current_loci(i_s_vals, base) mc.plot_flux_loci(i_s_vals, base) -plt.show() # %% # Set the speed reference and the external load torque. @@ -125,7 +142,7 @@ # Create the simulation object, simulate, and plot the results in per-unit values. sim = model.Simulation(mdl, ctrl) -res = sim.simulate(t_stop=1.5) +res = sim.simulate(t_stop=1.4) utils.plot(res, base) # %% diff --git a/examples/drive/flux_vector/plot_flux_vector_syrm_7kw_sat.py b/examples/drive/flux_vector/plot_flux_vector_syrm_7kw_sat.py index e3ae88dd6..c25bda692 100644 --- a/examples/drive/flux_vector/plot_flux_vector_syrm_7kw_sat.py +++ b/examples/drive/flux_vector/plot_flux_vector_syrm_7kw_sat.py @@ -41,8 +41,9 @@ psi_q_range = np.linspace(-0.5 * base.psi, 0.5 * base.psi, 256) flux_map = curr_map.as_magnetic_model(psi_d_range, psi_q_range).invert() -# Plot the flux map -utils.plot_maps(flux_map, base) +# Plot the flux maps +utils.plot_map(flux_map, "d", base) +utils.plot_map(flux_map, "q", base) # Parameter estimates est_par = control.SaturatedSynchronousMachinePars( diff --git a/examples/drive/vhz/plot_obs_vhz_pmsm_2kw_two_mass.py b/examples/drive/vhz/plot_obs_vhz_pmsm_2kw_two_mass.py index a53b8a620..18a898645 100644 --- a/examples/drive/vhz/plot_obs_vhz_pmsm_2kw_two_mass.py +++ b/examples/drive/vhz/plot_obs_vhz_pmsm_2kw_two_mass.py @@ -62,13 +62,13 @@ # %% # Plot the load speed and the twist angle. -t_span = (0, 1.2) +t_lims = (0, 1.2) _, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 5)) ax1.plot(res.mdl.t, res.mdl.mechanics.w_M, label=r"$\omega_\mathrm{M}$") ax1.plot(res.mdl.t, res.mdl.mechanics.w_L, label=r"$\omega_\mathrm{L}$") ax2.plot(res.mdl.t, res.mdl.mechanics.theta_ML * 180 / np.pi) -ax1.set_xlim(t_span) -ax2.set_xlim(t_span) +ax1.set_xlim(t_lims) +ax2.set_xlim(t_lims) ax1.set_xticklabels([]) ax1.legend() ax1.set_ylabel(r"$\omega_\mathrm{M}$, $\omega_\mathrm{L}$ (rad/s)") diff --git a/examples/drive/vhz/plot_vhz_im_2kw.py b/examples/drive/vhz/plot_vhz_im_2kw.py index 38ad42383..2d58a2f63 100644 --- a/examples/drive/vhz/plot_vhz_im_2kw.py +++ b/examples/drive/vhz/plot_vhz_im_2kw.py @@ -55,7 +55,9 @@ res = sim.simulate(t_stop=1.4) # sphinx_gallery_thumbnail_number = 2 utils.plot(res, base) -utils.plot_extra(res, base, t_span=(1.1, 1.125)) +utils.plot_stator_waveforms(res, base, t_lims=(1.1, 1.125)) +utils.plot_dc_bus_waveforms(res, base, t_lims=(1.1, 1.125)) + # %% # .. note:: diff --git a/examples/drive/vhz/plot_vhz_im_2kw_lc.py b/examples/drive/vhz/plot_vhz_im_2kw_lc.py index 40f7ac496..5f3236faa 100644 --- a/examples/drive/vhz/plot_vhz_im_2kw_lc.py +++ b/examples/drive/vhz/plot_vhz_im_2kw_lc.py @@ -61,20 +61,20 @@ # %% # Plot additional waveforms. -t_span = (1.1, 1.125) # Time span for the zoomed-in plot +t_lims = (1.1, 1.125) # Time span for the zoomed-in plot # Plot the converter and stator voltages (phase a) fig1, (ax1, ax2) = plt.subplots(2, 1) ax1.plot(res.mdl.t, res.mdl.converter.u_c_ab.real / base.u, label=r"$u_\mathrm{ca}$") ax1.plot(res.mdl.t, res.mdl.machine.u_s_ab.real / base.u, label=r"$u_\mathrm{sa}$") -ax1.set_xlim(t_span) +ax1.set_xlim(t_lims) ax1.legend() ax1.set_xticklabels([]) ax1.set_ylabel("Voltage (p.u.)") # Plot the converter and stator currents (phase a) ax2.plot(res.mdl.t, res.mdl.converter.i_c_ab.real / base.i, label=r"$i_\mathrm{ca}$") ax2.plot(res.mdl.t, res.mdl.machine.i_s_ab.real / base.i, label=r"$i_\mathrm{sa}$") -ax2.set_xlim(t_span) +ax2.set_xlim(t_lims) ax2.legend() ax2.set_ylabel("Current (p.u.)") ax2.set_xlabel("Time (s)") diff --git a/examples/grid/grid_following/plot_gfl_10kva.py b/examples/grid/grid_following/plot_gfl_10kva.py index 7bb350f59..225673149 100644 --- a/examples/grid/grid_following/plot_gfl_10kva.py +++ b/examples/grid/grid_following/plot_gfl_10kva.py @@ -50,7 +50,9 @@ sim = model.Simulation(mdl, ctrl) res = sim.simulate(t_stop=0.08) -utils.plot(res, base, plot_pcc_voltage=False) +utils.plot_control_signals(res, base) +utils.plot_grid_waveforms(res, base, plot_pcc_voltage=False) + # Uncomment line below to plot locus of the grid voltage space vector # utils.plot_voltage_vector(res, base) diff --git a/examples/grid/grid_following/plot_gfl_dc_bus_10kva.py b/examples/grid/grid_following/plot_gfl_dc_bus_10kva.py index 90d98452d..a3f4f8e8a 100644 --- a/examples/grid/grid_following/plot_gfl_dc_bus_10kva.py +++ b/examples/grid/grid_following/plot_gfl_dc_bus_10kva.py @@ -52,4 +52,5 @@ sim = model.Simulation(mdl, ctrl) res = sim.simulate(t_stop=0.1) -utils.plot(res, base) +utils.plot_control_signals(res, base) +utils.plot_grid_waveforms(res, base) diff --git a/examples/grid/grid_following/plot_gfl_lcl_10kva.py b/examples/grid/grid_following/plot_gfl_lcl_10kva.py index 87b5d2ed3..8dab6ab1c 100644 --- a/examples/grid/grid_following/plot_gfl_lcl_10kva.py +++ b/examples/grid/grid_following/plot_gfl_lcl_10kva.py @@ -5,7 +5,7 @@ This example simulates a grid-following-controlled converter connected to a strong grid through an LCL filter. The control system includes a phase-locked loop (PLL) to synchronize with the grid, a current reference generator, and a PI-type current -controller. The LCL filter dynamics are not taken into account in the control system. +controller. The LCL-filter dynamics are not taken into account in the control system. """ @@ -53,4 +53,5 @@ # %% # Plot the results. -utils.plot(res, base) +utils.plot_control_signals(res, base) +utils.plot_grid_waveforms(res, base) diff --git a/examples/grid/grid_forming/plot_gfm_obs_13kva.py b/examples/grid/grid_forming/plot_gfm_obs_13kva.py index 93284b02c..718221705 100644 --- a/examples/grid/grid_forming/plot_gfm_obs_13kva.py +++ b/examples/grid/grid_forming/plot_gfm_obs_13kva.py @@ -61,4 +61,5 @@ sim = model.Simulation(mdl, ctrl) res = sim.simulate(t_stop=1.4) -utils.plot(res, base) +utils.plot_control_signals(res, base) +utils.plot_grid_waveforms(res, base) diff --git a/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py b/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py index a431cef14..f7bef3e26 100644 --- a/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py +++ b/examples/grid/grid_forming/plot_gfm_rfpsc_13kva.py @@ -2,8 +2,8 @@ 12.5-kVA converter, RFPSC ========================= -This example simulates reference-feedforward power-synchronization control -(RFPSC) of a converter connected to a weak grid. +This example simulates reference-feedforward power-synchronization control (RFPSC) of a +converter connected to a weak grid. """ @@ -49,4 +49,5 @@ sim = model.Simulation(mdl, ctrl) res = sim.simulate(t_stop=1.4) -utils.plot(res, base) +utils.plot_control_signals(res, base) +utils.plot_grid_waveforms(res, base) diff --git a/motulator/common/model/_converter.py b/motulator/common/model/_converter.py index 7295894f1..a0754e85a 100644 --- a/motulator/common/model/_converter.py +++ b/motulator/common/model/_converter.py @@ -14,7 +14,7 @@ import numpy as np from motulator.common.model._base import Subsystem -from motulator.common.utils import abc2complex, complex2abc, empty_array +from motulator.common.utils._utils import abc2complex, complex2abc, empty_array # %% diff --git a/motulator/common/utils/__init__.py b/motulator/common/utils/__init__.py index e4a835323..3451b630c 100644 --- a/motulator/common/utils/__init__.py +++ b/motulator/common/utils/__init__.py @@ -4,22 +4,7 @@ SequenceGenerator, Step, abc2complex, - clip, complex2abc, - empty_array, - get_value, - sign, - wrap, ) -__all__ = [ - "SequenceGenerator", - "Step", - "abc2complex", - "complex2abc", - "clip", - "empty_array", - "get_value", - "sign", - "wrap", -] +__all__ = ["SequenceGenerator", "Step", "abc2complex", "complex2abc"] diff --git a/motulator/common/utils/_plotting.py b/motulator/common/utils/_plotting.py new file mode 100644 index 000000000..9623685dc --- /dev/null +++ b/motulator/common/utils/_plotting.py @@ -0,0 +1,124 @@ +"""Helper functions for plots.""" + +import matplotlib.pyplot as plt +from cycler import cycler +from numpy.typing import ArrayLike + + +# %% +def set_screen_style() -> None: + """Configure matplotlib for screen viewing.""" + plt.style.use("default") + # plt.style.use("dark_background") + plt.rcParams.update( + { + "axes.prop_cycle": cycler("color", ["b", "r", "m", "c"]), + "lines.linewidth": 1, + "axes.grid": True, + "text.usetex": False, + "axes3d.grid": True, + "figure.constrained_layout.use": True, + } + ) + + +# %% +def set_latex_style() -> None: + """Configure matplotlib for LaTeX documents.""" + plt.style.use("default") + plt.rcParams.update( + { + "axes.prop_cycle": cycler("color", ["b", "r", "m", "c"]), + "lines.linewidth": 1, + "axes.grid": True, + # "axes.autolimit_mode": "round_numbers", + "text.usetex": True, + "text.latex.preamble": (r"\usepackage{newtxtext}\usepackage{newtxmath}"), + "font.family": "serif", + "font.serif": ["Liberation Serif"], + "mathtext.fontset": "dejavuserif", + "font.size": 8, + "axes.labelsize": 8, + "legend.fontsize": 8, + "xtick.labelsize": 8, + "ytick.labelsize": 8, + "axes3d.grid": True, + "figure.figsize": [3.5, 2.5], + # "figure.constrained_layout.use": True, + # "savefig.bbox": "tight", + # "savefig.pad_inches": 0.01, + # Disable automatic adjustment + "figure.constrained_layout.use": False, + "figure.subplot.left": 0.18, # Fixed margins + "figure.subplot.right": 0.95, + "figure.subplot.bottom": 0.18, + "figure.subplot.top": 0.95, + } + ) + + +# %% +def configure_axes( + axes: list, + t_lims: tuple[float, float] | None, + t_ticks: ArrayLike | None, + y_lims: list[tuple[float, float] | None] | None, + y_ticks: list[ArrayLike | None] | None, +) -> None: + """ + Configure time and y-axis settings for all subplots. + + Parameters + ---------- + axes : list + List of matplotlib axes objects to configure. + t_lims : tuple[float, float] | None + Time axis limits. If None, no limits are set. + t_ticks : ArrayLike | None + Time axis tick locations. If None, no ticks are set. + y_lims : list[tuple[float, float] | None] | None + y-axis limits for each subplot. Each element should be a tuple (y_min, y_max) or + None for automatic scaling. + y_ticks : list of ArrayLike or None + y-axis tick locations for each subplot. Each element should be an array of tick + positions or None for automatic ticks. + + """ + # Set time limits and ticks for all subplots + for ax in axes: + if t_lims is not None: + ax.set_xlim(t_lims) + if t_ticks is not None: + ax.set_xticks(t_ticks) + + # Apply y-axis limits + if y_lims is not None: + for ax, ylim in zip(axes, y_lims, strict=False): + if ylim is not None: + ax.set_ylim(ylim) + + # Apply y-axis ticks + if y_ticks is not None: + for ax, ytick in zip(axes, y_ticks, strict=False): + if ytick is not None: + ax.set_yticks(ytick) + + +# %% +def save_and_show(save_path: str | None = None, **savefig_kwargs) -> None: + """ + Helper function to optionally save and show plots. + + Parameters + ---------- + save_path : str | None, optional + Path to save the figure. If None, the figure is not saved. + savefig_kwargs : dict, optional + Additional keyword arguments for `plt.savefig`. + + """ + if save_path is not None: + default_kwargs = {"bbox_inches": "tight"} + default_kwargs.update(savefig_kwargs) + plt.savefig(save_path, **default_kwargs) + plt.show() diff --git a/motulator/common/utils/_utils.py b/motulator/common/utils/_utils.py index 36033a594..a83d75f2f 100644 --- a/motulator/common/utils/_utils.py +++ b/motulator/common/utils/_utils.py @@ -3,9 +3,7 @@ from dataclasses import dataclass from typing import Any, Callable -import matplotlib.pyplot as plt import numpy as np -from cycler import cycler # %% @@ -323,55 +321,3 @@ def clip(value: float, min_value: float, max_value: float) -> float: def sign(x: float) -> int: """Return the sign of x: -1 for negative, 0 for zero, 1 for positive.""" return -1 if x < 0 else (1 if x > 0 else 0) - - -# %% -def set_screen_style() -> None: - """Configure matplotlib for screen viewing.""" - plt.style.use("default") - # plt.style.use("dark_background") - plt.rcParams.update( - { - "axes.prop_cycle": cycler("color", ["b", "r", "m", "c"]), - "lines.linewidth": 1, - "axes.grid": True, - "text.usetex": False, - "axes3d.grid": True, - "figure.constrained_layout.use": True, - } - ) - - -# %% -def set_latex_style() -> None: - """Configure matplotlib for LaTeX documents.""" - plt.style.use("default") - plt.rcParams.update( - { - "axes.prop_cycle": cycler("color", ["b", "r", "m", "c"]), - "lines.linewidth": 1, - "axes.grid": True, - # "axes.autolimit_mode": "round_numbers", - "text.usetex": True, - "text.latex.preamble": (r"\usepackage{newtxtext}\usepackage{newtxmath}"), - "font.family": "serif", - "font.serif": ["Liberation Serif"], - "mathtext.fontset": "dejavuserif", - "font.size": 8, - "axes.labelsize": 8, - "legend.fontsize": 8, - "xtick.labelsize": 8, - "ytick.labelsize": 8, - "axes3d.grid": True, - "figure.figsize": [3.5, 2.5], - # "figure.constrained_layout.use": True, - # "savefig.bbox": "tight", - # "savefig.pad_inches": 0.01, - # Disable automatic adjustment - "figure.constrained_layout.use": False, - "figure.subplot.left": 0.18, # Fixed margins - "figure.subplot.right": 0.95, - "figure.subplot.bottom": 0.18, - "figure.subplot.top": 0.95, - } - ) diff --git a/motulator/drive/control/_im_current_vector.py b/motulator/drive/control/_im_current_vector.py index c10564783..ff75fbba4 100644 --- a/motulator/drive/control/_im_current_vector.py +++ b/motulator/drive/control/_im_current_vector.py @@ -7,7 +7,7 @@ from motulator.common.control import ComplexPIController from motulator.common.control._pwm import PWM -from motulator.common.utils import clip +from motulator.common.utils._utils import clip from motulator.drive.control._base import Measurements from motulator.drive.control._im_observers import ( ObserverOutputs, diff --git a/motulator/drive/control/_im_flux_vector.py b/motulator/drive/control/_im_flux_vector.py index ed76d27a4..e009c9590 100644 --- a/motulator/drive/control/_im_flux_vector.py +++ b/motulator/drive/control/_im_flux_vector.py @@ -6,7 +6,7 @@ from typing import Callable, Sequence from motulator.common.control._pwm import PWM -from motulator.common.utils import sign +from motulator.common.utils._utils import sign from motulator.drive.control._base import Measurements from motulator.drive.control._im_observers import ( ObserverOutputs, diff --git a/motulator/drive/control/_im_observers.py b/motulator/drive/control/_im_observers.py index b914a5581..9c7c5cdba 100644 --- a/motulator/drive/control/_im_observers.py +++ b/motulator/drive/control/_im_observers.py @@ -5,7 +5,7 @@ from math import inf, pi from typing import Callable -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.drive.utils._parameters import InductionMachineInvGammaPars diff --git a/motulator/drive/control/_sm_observers.py b/motulator/drive/control/_sm_observers.py index d0b02b822..54b492afb 100644 --- a/motulator/drive/control/_sm_observers.py +++ b/motulator/drive/control/_sm_observers.py @@ -5,7 +5,7 @@ from math import pi from typing import Callable -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.drive.utils._parameters import ( SaturatedSynchronousMachinePars, SynchronousMachinePars, diff --git a/motulator/drive/control/_sm_reference_gen.py b/motulator/drive/control/_sm_reference_gen.py index 367301a5f..a6ff921f4 100644 --- a/motulator/drive/control/_sm_reference_gen.py +++ b/motulator/drive/control/_sm_reference_gen.py @@ -5,7 +5,7 @@ from scipy.optimize import root_scalar -from motulator.common.utils import clip, sign +from motulator.common.utils._utils import clip, sign from motulator.drive.utils._parameters import ( SaturatedSynchronousMachinePars, SynchronousMachinePars, diff --git a/motulator/drive/control/_sm_signal_inj.py b/motulator/drive/control/_sm_signal_inj.py index 077cf5c1a..65cf91642 100644 --- a/motulator/drive/control/_sm_signal_inj.py +++ b/motulator/drive/control/_sm_signal_inj.py @@ -3,7 +3,7 @@ from cmath import exp, pi from dataclasses import dataclass -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.drive.control._sm_current_vector import ( CurrentVectorController, CurrentVectorControllerCfg, diff --git a/motulator/drive/model/_lc_filter.py b/motulator/drive/model/_lc_filter.py index 16738f66d..064057580 100644 --- a/motulator/drive/model/_lc_filter.py +++ b/motulator/drive/model/_lc_filter.py @@ -6,7 +6,7 @@ import numpy as np from motulator.common.model import Subsystem, SubsystemTimeSeries -from motulator.common.utils import complex2abc, empty_array +from motulator.common.utils._utils import complex2abc, empty_array # %% diff --git a/motulator/drive/model/_machine.py b/motulator/drive/model/_machine.py index 98c8b37b3..8a5308e9f 100644 --- a/motulator/drive/model/_machine.py +++ b/motulator/drive/model/_machine.py @@ -12,7 +12,7 @@ import numpy as np from motulator.common.model import Subsystem, SubsystemTimeSeries -from motulator.common.utils import complex2abc, empty_array, get_value +from motulator.common.utils._utils import complex2abc, empty_array, get_value from motulator.drive.utils._parameters import ( InductionMachineInvGammaPars, InductionMachinePars, diff --git a/motulator/drive/model/_mechanics.py b/motulator/drive/model/_mechanics.py index 4c63086ea..f52cf8b2e 100644 --- a/motulator/drive/model/_mechanics.py +++ b/motulator/drive/model/_mechanics.py @@ -7,7 +7,7 @@ import numpy as np from motulator.common.model import Subsystem, SubsystemTimeSeries -from motulator.common.utils import empty_array, get_value +from motulator.common.utils._utils import empty_array, get_value # %% diff --git a/motulator/drive/utils/__init__.py b/motulator/drive/utils/__init__.py index 4a4a5009d..c3dfa9b79 100644 --- a/motulator/drive/utils/__init__.py +++ b/motulator/drive/utils/__init__.py @@ -6,17 +6,20 @@ SequenceGenerator, Step, ) -from motulator.drive.utils._plots import plot, plot_extra +from motulator.drive.utils._plots import ( + plot, + plot_dc_bus_waveforms, + plot_stator_waveforms, +) from motulator.drive.utils._sm_control_loci import ControlLoci from motulator.drive.utils._sm_flux_maps import ( MagneticModel, SaturationModelPMSyRM, SaturationModelSyRM, import_syre_data, - plot_flux_vs_current, - plot_maps, ) from motulator.drive.utils._sm_plot_control_loci import MachineCharacteristics +from motulator.drive.utils._sm_plot_flux_maps import plot_flux_vs_current, plot_map __all__ = [ "BaseValues", @@ -26,9 +29,10 @@ "MagneticModel", "NominalValues", "plot", - "plot_extra", + "plot_stator_waveforms", + "plot_dc_bus_waveforms", "plot_flux_vs_current", - "plot_maps", + "plot_map", "SaturationModelPMSyRM", "SaturationModelSyRM", "SequenceGenerator", diff --git a/motulator/drive/utils/_plots.py b/motulator/drive/utils/_plots.py index 64dedbe02..023e2d1c9 100644 --- a/motulator/drive/utils/_plots.py +++ b/motulator/drive/utils/_plots.py @@ -4,45 +4,29 @@ import matplotlib.pyplot as plt import numpy as np +from numpy.typing import ArrayLike from motulator.common.model._simulation import SimulationResults -from motulator.common.utils._utils import ( - BaseValues, - complex2abc, +from motulator.common.utils._plotting import ( + configure_axes, + save_and_show, set_latex_style, set_screen_style, ) +from motulator.common.utils._utils import BaseValues, complex2abc - -def _get_machine_type(mdl: Any) -> Literal["im", "sm"]: - return "im" if hasattr(mdl.machine, "psi_R_ab") else "sm" +# LaTeX symbol definition for easier configuration. When per-unit values are used, +# sometimes subscript m is preferred instead of M. +M = r"\mathrm{M}" # %% -def plot( - res: SimulationResults, - base: BaseValues | None = None, - t_span: tuple[float, float] | None = None, - latex: bool = False, -) -> None: - """ - Plot example figures. +def _get_machine_type(mdl: Any) -> Literal["im", "sm"]: + return "im" if hasattr(mdl.machine, "psi_R_ab") else "sm" - Parameters - ---------- - res : SimulationResults - Simulation results. - base : BaseValues, optional - Base values for scaling the waveforms. If not given, the waveforms are plotted - in SI units. - t_span : 2-tuple, optional - Time span. If not given, the whole simulation time is plotted. - latex : bool, optional - Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX - installation, defaults to False. - """ - # ruff: noqa: PLR0912 +def _setup_plot(latex: bool) -> tuple[float, float]: + """Setup plot style and return figure dimensions.""" if latex: set_latex_style() height = plt.rcParams["figure.figsize"][1] * 2.2 @@ -51,151 +35,210 @@ def plot( height = plt.rcParams["figure.figsize"][1] * 1.8 width = plt.rcParams["figure.figsize"][0] - fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, figsize=(width, height)) + return width, height - # For brevity - mdl = res.mdl # Continuous-time data - ctrl = res.ctrl # Discrete-time data - # Check if the time span was given - if t_span is None: - t_span = (0, mdl.t[-1]) - - # Check if the base values were given - if base is None: - pu_vals = False - base = BaseValues.unity() - else: - pu_vals = True - - # Subplot 1: Angular speeds +# %% +def _plot_speeds(ax, ctrl, mdl, base: BaseValues) -> None: + """Plot angular speeds.""" if hasattr(ctrl.ref, "w_M") and np.any(ctrl.ref.w_M): - ax1.plot( + ax.plot( ctrl.t, ctrl.ref.w_M / base.w_M, "--", ds="steps-post", - label=r"$\omega_\mathrm{M,ref}$", + label=rf"$\omega_{M}^\mathrm{{ref}}$", ) - ax1.plot(mdl.t, mdl.machine.w_M / base.w_M, label=r"$\omega_\mathrm{M}$") - ax1.plot( - ctrl.t, - ctrl.fbk.w_M / base.w_M, - label=r"$\hat \omega_\mathrm{M}$", - ds="steps-post", + ax.plot(mdl.t, mdl.machine.w_M / base.w_M, label=rf"$\omega_{M}$") + ax.plot( + ctrl.t, ctrl.fbk.w_M / base.w_M, label=rf"$\hat{{\omega}}_{M}$", ds="steps-post" ) - ax1.legend() - ax1.set_xlim(t_span) - ax1.set_xticklabels([]) + ax.legend() - # Subplot 2: torques + +def _plot_torques(ax, ctrl, mdl, base: BaseValues) -> None: + """Plot torques.""" # Plot torque reference only if it is nonzero if hasattr(ctrl.ref, "tau_M") and np.any(ctrl.ref.tau_M): - ax2.plot( + ax.plot( ctrl.t, ctrl.ref.tau_M / base.tau, "--", - label=r"$\tau_\mathrm{M,ref}$", + label=rf"$\tau_{M}^\mathrm{{ref}}$", ds="steps-post", ) - ax2.plot(mdl.t, mdl.machine.tau_M / base.tau, label=r"$\tau_\mathrm{M}$") + ax.plot(mdl.t, mdl.machine.tau_M / base.tau, label=rf"$\tau_{M}$") # Plot torque estimated only if it is nonzero if hasattr(ctrl.fbk, "tau_M") and np.any(ctrl.fbk.tau_M): - ax2.plot( + ax.plot( ctrl.t, ctrl.fbk.tau_M / base.tau, - label=r"$\hat{\tau}_\mathrm{M}$", + label=rf"$\hat{{\tau}}_{M}$", ds="steps-post", ) if hasattr(mdl.mechanics, "tau_L_tot"): - ax2.plot( + ax.plot( mdl.t, mdl.mechanics.tau_L_tot / base.tau, ":", - label=r"$\tau_\mathrm{L,tot}$", + label=r"$\tau_\mathrm{L}^\mathrm{tot}$", ) - ax2.legend(loc="upper right") - ax2.set_xlim(t_span) - ax2.set_xticklabels([]) + ax.legend(loc="upper right") + - # Subplot 3: Currents +def _plot_currents(ax, ctrl, base: BaseValues) -> None: + """Plot currents.""" if hasattr(ctrl.ref, "i_s") and np.any(ctrl.ref.i_s): - ax3.plot( + ax.plot( ctrl.t, ctrl.ref.i_s.real / base.i, "--", - label=r"$i_\mathrm{d,ref}$", + label=r"$i_\mathrm{d}^\mathrm{ref}$", ds="steps-post", ) - ax3.plot( + ax.plot( ctrl.t, ctrl.fbk.i_s.real / base.i, label=r"$i_\mathrm{d}$", ds="steps-post" ) if hasattr(ctrl.ref, "i_s") and np.any(ctrl.ref.i_s): - ax3.plot( + ax.plot( ctrl.t, ctrl.ref.i_s.imag / base.i, "--", - label=r"$i_\mathrm{q,ref}$", + label=r"$i_\mathrm{q}^\mathrm{ref}$", ds="steps-post", ) - ax3.plot( + ax.plot( ctrl.t, ctrl.fbk.i_s.imag / base.i, label=r"$i_\mathrm{q}$", ds="steps-post" ) - ax3.legend(loc="upper right") - ax3.set_xlim(t_span) - ax3.set_xticklabels([]) + ax.legend(loc="upper right") - # Subplot 4: Voltages - ax4.plot( - ctrl.t, np.abs(ctrl.fbk.u_s) / base.u, label=r"$u_\mathrm{s}$", ds="steps-post" + +def _plot_voltages(ax, ctrl, base: BaseValues) -> None: + """Plot voltages.""" + ax.plot( + ctrl.t, + np.abs(ctrl.fbk.u_s) / base.u, + label=r"$\hat{u}_\mathrm{s}$", + ds="steps-post", ) - ax4.plot( + ax.plot( ctrl.t, ctrl.fbk.u_dc / np.sqrt(3) / base.u, "--", label=r"$u_\mathrm{dc}/\sqrt{3}$", ds="steps-post", ) - ax4.legend(loc="lower right") - ax4.set_xlim(t_span) - _, ymax = ax4.get_ylim() - ax4.set_ylim(0, ymax) - ax4.set_xticklabels([]) + ax.legend(loc="lower right") + _, ymax = ax.get_ylim() + ax.set_ylim(0, ymax) - # Subplot 5: Flux linkages + +def _plot_flux_linkages(ax, ctrl, mdl, base: BaseValues) -> None: + """Plot flux linkages.""" if hasattr(ctrl.ref, "psi_s"): - ax5.plot( + ax.plot( ctrl.t, np.abs(ctrl.ref.psi_s) / base.psi, "--", - label=r"$\psi_\mathrm{s,ref}$", + label=r"$\psi_\mathrm{s}^\mathrm{ref}$", ds="steps-post", ) if _get_machine_type(mdl) == "sm": - ax5.plot( + ax.plot( mdl.t, np.abs(mdl.machine.psi_s_ab) / base.psi, label=r"$\psi_\mathrm{s}$" ) else: - ax5.plot( + ax.plot( mdl.t, np.abs(mdl.machine.psi_s_ab) / base.psi, label=r"$\psi_\mathrm{s}$" ) - ax5.plot( + ax.plot( mdl.t, np.abs(mdl.machine.psi_R_ab) / base.psi, "-.", label=r"$\psi_\mathrm{R}$", ) if hasattr(ctrl.fbk, "psi_s"): - ax5.plot( + ax.plot( ctrl.t, np.abs(ctrl.fbk.psi_s) / base.psi, - label=r"$\hat\psi_\mathrm{s}$", + label=r"$\hat{\psi}_\mathrm{s}$", ds="steps-post", ) - ax5.legend(loc="upper right") - ax5.set_xlim(t_span) - _, ymax = ax5.get_ylim() - ax5.set_ylim(0, ymax) + ax.legend(loc="upper right") + _, ymax = ax.get_ylim() + ax.set_ylim(0, ymax) + + +def plot( + res: SimulationResults, + base: BaseValues | None = None, + t_lims: tuple[float, float] | None = None, + t_ticks: ArrayLike | None = None, + y_lims: list[tuple[float, float] | None] | None = None, + y_ticks: list[ArrayLike | None] | None = None, + latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot example figures. + + Parameters + ---------- + res : SimulationResults + Simulation results. + base : BaseValues, optional + Base values for scaling the waveforms. If not given, the waveforms are plotted + in SI units. + t_lims : tuple[float, float], optional + Time axis limits. If None, uses full time range. + t_ticks : ArrayLike, optional + Time axis tick locations. + y_lims : list[tuple[float, float] | None], optional + y-axis limits for each subplot. + y_ticks : list[ArrayLike | None], optional + y-axis tick locations for each subplot. + latex : bool, optional + Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX + installation, defaults to False. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + # Setup plot style and dimensions + width, height = _setup_plot(latex) + + # Process base values + if base is None: + pu_vals = False + base = BaseValues.unity() + else: + pu_vals = True + + # Set time limits + if t_lims is None: + t_lims = (0, res.mdl.t[-1]) + + # Create figure + fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, figsize=(width, height)) + axes = [ax1, ax2, ax3, ax4, ax5] + + # Plot subplots + _plot_speeds(ax1, res.ctrl, res.mdl, base) + _plot_torques(ax2, res.ctrl, res.mdl, base) + _plot_currents(ax3, res.ctrl, base) + _plot_voltages(ax4, res.ctrl, base) + _plot_flux_linkages(ax5, res.ctrl, res.mdl, base) + + # Configure all axes + configure_axes(axes, t_lims, t_ticks, y_lims, y_ticks) + + # Remove xticklabels for all but the last subplot + for ax in axes[:-1]: + ax.set_xticklabels([]) + ax5.set_xlabel("Time (s)") # Add axis labels if pu_vals: @@ -210,93 +253,111 @@ def plot( ax3.set_ylabel("Current (A)") ax4.set_ylabel("Voltage (V)") ax5.set_ylabel("Flux linkage (Vs)") - ax5.set_xlabel("Time (s)") fig.align_ylabels() - plt.show() + save_and_show(save_path, **savefig_kwargs) # %% -def plot_extra( +def _plot_phase_voltages(ax, mdl, ctrl, base: BaseValues) -> None: + """Plot phase voltages.""" + if _get_machine_type(mdl) == "sm": + theta = ctrl.fbk.theta_m + else: + theta = ctrl.fbk.theta_s + u_s_ab = np.exp(1j * theta) * ctrl.fbk.u_s + + ax.plot(mdl.t, mdl.machine.u_s_ab.real / base.u, label=r"$u_\mathrm{sa}$") + ax.plot( + ctrl.t, u_s_ab.real / base.u, label=r"$\hat{u}_\mathrm{sa}$", ds="steps-post" + ) + ax.legend() + + +def _plot_phase_currents(ax, mdl, ctrl, base: BaseValues) -> None: + """Plot phase currents.""" + if _get_machine_type(mdl) == "sm": + theta = ctrl.fbk.theta_m + else: + theta = ctrl.fbk.theta_s + i_s_ab = np.exp(1j * theta) * ctrl.fbk.i_s + + ax.plot( + mdl.t, + complex2abc(mdl.machine.i_s_ab).T / base.i, + label=[r"$i_\mathrm{sa}$", r"$i_\mathrm{sb}$", r"$i_\mathrm{sc}$"], + ) + ax.plot(ctrl.t, i_s_ab.real / base.i, label=r"$i_\mathrm{sa}(k)$", ds="steps-post") + ax.legend() + + +def plot_stator_waveforms( res: SimulationResults, base: BaseValues | None = None, - t_span: tuple[float, float] | None = None, + t_lims: tuple[float, float] | None = None, + t_ticks: ArrayLike | None = None, + y_lims: list[tuple[float, float] | None] | None = None, + y_ticks: list[ArrayLike | None] | None = None, latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, ) -> None: """ - Plot extra waveforms for a motor drive with a diode bridge. + Plot stator voltage and current waveforms. Parameters ---------- - res : Simulation + res : SimulationResults Should contain the simulated data. base : BaseValues, optional Base values for scaling the waveforms. - t_span : 2-tuple, optional - Time span, defaults to (0, sim.ctrl.t[-1]). + t_lims : tuple[float, float], optional + Time axis limits. If None, uses full time range. + t_ticks : ArrayLike, optional + Time axis tick locations. + y_lims : list[tuple[float, float] | None], optional + y-axis limits for each subplot. + y_ticks : list[ArrayLike | None], optional + y-axis tick locations for each subplot. latex : bool, optional Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ - # ruff: noqa: PLR0915 + # Setup plot style if latex: set_latex_style() else: set_screen_style() - # For brevity - mdl = res.mdl # Continuous-time data - ctrl = res.ctrl # Discrete-time data - - # Check if the time span was given - if t_span is None: - t_span = (0, mdl.t[-1]) - - # Check if the base values were given + # Process base values if base is None: pu_vals = False base = BaseValues.unity() else: pu_vals = True - # Angle of synchronous coordinates - if _get_machine_type(mdl) == "sm": - theta = ctrl.fbk.theta_m - else: - theta = ctrl.fbk.theta_s + # Set time limits + if t_lims is None: + t_lims = (0, res.mdl.t[-1]) - # Quantities in stator coordinates - ctrl.fbk.u_s_ab = np.exp(1j * theta) * ctrl.fbk.u_s - ctrl.fbk.i_s_ab = np.exp(1j * theta) * ctrl.fbk.i_s + # Create figure + fig, (ax1, ax2) = plt.subplots(2, 1) + axes = [ax1, ax2] - fig1, (ax1, ax2) = plt.subplots(2, 1) + # Plot subplots + _plot_phase_voltages(ax1, res.mdl, res.ctrl, base) + _plot_phase_currents(ax2, res.mdl, res.ctrl, base) - # Subplot 1: Voltages - ax1.plot(mdl.t, mdl.machine.u_s_ab.real / base.u, label=r"$u_\mathrm{sa}$") - ax1.plot( - ctrl.t, - ctrl.fbk.u_s_ab.real / base.u, - label=r"$\hat u_\mathrm{sa}$", - ds="steps-post", - ) - ax1.set_xlim(t_span) - ax1.legend() - ax1.set_xticklabels([]) + # Configure all axes + configure_axes(axes, t_lims, t_ticks, y_lims, y_ticks) - # Subplot 2: currents - ax2.plot( - mdl.t, - complex2abc(mdl.machine.i_s_ab).T / base.i, - label=[r"$i_\mathrm{sa}$", r"$i_\mathrm{sb}$", r"$i_\mathrm{sc}$"], - ) - ax2.plot( - ctrl.t, - ctrl.fbk.i_s_ab.real / base.i, - label=r"$\hat i_\mathrm{sa}$", - ds="steps-post", - ) - ax2.set_xlim(t_span) - ax2.legend() + # Remove xticklabels for the first subplot + ax1.set_xticklabels([]) + ax2.set_xlabel("Time (s)") # Add axis labels if pu_vals: @@ -305,38 +366,108 @@ def plot_extra( else: ax1.set_ylabel("Voltage (V)") ax2.set_ylabel("Current (A)") + fig.align_ylabels() + + save_and_show(save_path, **savefig_kwargs) + + +# %% +def _plot_dc_bus_voltages(ax, mdl, base: BaseValues) -> None: + """Plot DC bus and grid voltages.""" + ax.plot(mdl.t, mdl.converter.u_di / base.u, label=r"$u_\mathrm{di}$") + ax.plot(mdl.t, mdl.converter.u_dc / base.u, label=r"$u_\mathrm{dc}$") + ax.plot(mdl.t, mdl.converter.u_g_abc.T[:, 0] / base.u, label=r"$u_\mathrm{ga}$") + ax.legend() + + +def _plot_dc_bus_currents(ax, mdl, base: BaseValues) -> None: + """Plot DC bus and grid currents.""" + ax.plot(mdl.t, mdl.converter.i_L / base.i, label=r"$i_\mathrm{L}$") + ax.plot(mdl.t, mdl.converter.i_dc_int / base.i, label=r"$i_\mathrm{dc}$") + ax.plot(mdl.t, mdl.converter.i_g_ab.real / base.i, label=r"$i_\mathrm{ga}$") + ax.legend() + + +def plot_dc_bus_waveforms( + res: SimulationResults, + base: BaseValues | None = None, + t_lims: tuple[float, float] | None = None, + t_ticks: ArrayLike | None = None, + y_lims: list[tuple[float, float] | None] | None = None, + y_ticks: list[ArrayLike | None] | None = None, + latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot DC bus and grid-side waveforms. + + Parameters + ---------- + res : SimulationResults + Should contain the simulated data. + base : BaseValues, optional + Base values for scaling the waveforms. + t_lims : tuple[float, float], optional + Time axis limits. If None, uses full time range. + t_ticks : ArrayLike, optional + Time axis tick locations. + y_lims : list[tuple[float, float] | None], optional + y-axis limits for each subplot. + y_ticks : list[ArrayLike | None], optional + y-axis tick locations for each subplot. + latex : bool, optional + Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + # Check if DC bus data exists + if not hasattr(res.mdl.converter, "i_L"): + print("DC bus data not available in simulation results.") + return + + # Setup plot style + if latex: + set_latex_style() + else: + set_screen_style() + + # Process base values + if base is None: + pu_vals = False + base = BaseValues.unity() + else: + pu_vals = True + + # Set time limits + if t_lims is None: + t_lims = (0, res.mdl.t[-1]) + + # Create figure + fig, (ax1, ax2) = plt.subplots(2, 1) + axes = [ax1, ax2] + + # Plot subplots + _plot_dc_bus_voltages(ax1, res.mdl, base) + _plot_dc_bus_currents(ax2, res.mdl, base) + + # Configure all axes + configure_axes(axes, t_lims, t_ticks, y_lims, y_ticks) + + # Remove xticklabels for first subplot + ax1.set_xticklabels([]) ax2.set_xlabel("Time (s)") - fig1.align_ylabels() - # Plots the DC bus and grid-side variables (if exist) - if hasattr(mdl.converter, "i_L"): - fig2, (ax1, ax2) = plt.subplots(2, 1) + # Add axis labels + if pu_vals: + ax1.set_ylabel("Voltage (p.u.)") + ax2.set_ylabel("Current (p.u.)") + else: + ax1.set_ylabel("Voltage (V)") + ax2.set_ylabel("Current (A)") + fig.align_ylabels() - # Subplot 1: Voltages - ax1.plot(mdl.t, mdl.converter.u_di / base.u, label=r"$u_\mathrm{di}$") - ax1.plot(mdl.t, mdl.converter.u_dc / base.u, label=r"$u_\mathrm{dc}$") - ax1.plot( - mdl.t, mdl.converter.u_g_abc.T[:, 0] / base.u, label=r"$u_\mathrm{ga}$" - ) - ax1.legend() - ax1.set_xlim(t_span) - ax1.set_xticklabels([]) - - # Subplot 2: Currents - ax2.plot(mdl.t, mdl.converter.i_L / base.i, label=r"$i_\mathrm{L}$") - ax2.plot(mdl.t, mdl.converter.i_dc_int / base.i, label=r"$i_\mathrm{dc}$") - ax2.plot(mdl.t, mdl.converter.i_g_ab.real / base.i, label=r"$i_\mathrm{ga}$") - ax2.legend() - ax2.set_xlim(t_span) - - # Add axis labels - if pu_vals: - ax1.set_ylabel("Voltage (p.u.)") - ax2.set_ylabel("Current (p.u.)") - else: - ax1.set_ylabel("Voltage (V)") - ax2.set_ylabel("Current (A)") - - fig2.align_ylabels() - - plt.show() + save_and_show(save_path, **savefig_kwargs) diff --git a/motulator/drive/utils/_sm_flux_maps.py b/motulator/drive/utils/_sm_flux_maps.py index 43f14e248..c9b2e0dd8 100644 --- a/motulator/drive/utils/_sm_flux_maps.py +++ b/motulator/drive/utils/_sm_flux_maps.py @@ -1,16 +1,14 @@ -"""Manipulate and plot flux linkage maps of synchronous machines.""" +"""Manipulate flux linkage and current maps of synchronous machines.""" from dataclasses import dataclass -from typing import Any, Callable, Literal, cast +from typing import Callable, Literal, cast -import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import RegularGridInterpolator, griddata from scipy.io import loadmat -from motulator.common.utils._utils import BaseValues, set_latex_style, set_screen_style - +# %% @dataclass class MagneticModel: """ @@ -93,10 +91,10 @@ def create_interpolated_model( Parameters ---------- - d_range : Any, optional + d_range : np.ndarray | None, optional Range for the d-axis. If None, the range is determined from the data, defaults to None. - q_range : Any, optional + q_range : np.ndarray | None, optional Range for the q-axis. If None, the range is determined from the data, defaults to None. num : int, optional @@ -129,15 +127,15 @@ def create_grid(d_range: np.ndarray, q_range: np.ndarray) -> np.ndarray: new_type = "current_map" if self.is_flux_map() else "flux_map" # Auto-determine ranges if not provided - psi_d_min = np.max(np.min(np.real(inp), axis=0)) - psi_d_max = np.min(np.max(np.real(inp), axis=0)) - psi_q_min = np.max(np.min(np.imag(inp), axis=1)) - psi_q_max = np.min(np.max(np.imag(inp), axis=1)) + d_min = np.max(np.min(np.real(inp), axis=0)) + d_max = np.min(np.max(np.real(inp), axis=0)) + q_min = np.max(np.min(np.imag(inp), axis=1)) + q_max = np.min(np.max(np.imag(inp), axis=1)) if d_range is None: - d_range = np.linspace(psi_d_min, psi_d_max, num) + d_range = np.linspace(d_min, d_max, num) if q_range is None: - q_range = np.linspace(psi_q_min, psi_q_max, num) + q_range = np.linspace(q_min, q_max, num) # Interpolate the map new_inp = create_grid(d_range, q_range) # type: ignore @@ -206,17 +204,20 @@ def lookup_fcn(dq_input: complex | np.ndarray) -> complex | np.ndarray: ) def invert( - self, d_range: Any = None, q_range: Any = None, num: int | None = None + self, + d_range: np.ndarray | None = None, + q_range: np.ndarray | None = None, + num: int | None = None, ) -> "MagneticModel": """ Invert the map (swap input and output). Parameters ---------- - d_range : Any, optional + d_range : np.ndarray | None, optional Range for the d-axis. If None, the range is determined from the data, defaults to None. - q_range : Any, optional + q_range : np.ndarray | None, optional Range for the q-axis. If None, the range is determined from the data, defaults to None. num : int, optional @@ -432,256 +433,6 @@ def __call__(self, psi_s_dq: complex | np.ndarray) -> complex | np.ndarray: return i_s_dq -# %% -# pyright: reportPrivateImportUsage=false -def plot_maps( - data: MagneticModel, - base: BaseValues | None = None, - x_lims: tuple[float, float] | None = None, - y_lims: tuple[float, float] | None = None, - z_lims: tuple[float, float] | None = None, - raw_data: MagneticModel | None = None, - latex: bool = False, -) -> None: - # ruff: noqa: PLR0912, PLR0913, PLR0915 - """ - Plot flux maps and current maps. - - Parameters - ---------- - data : MagneticModel - Data containing the flux and current information. The coordinates are selected - based on the `type` field, which is either "flux_map" or "current_map" (the - default is "flux_map"). - base : BaseValues, optional - Base values for scaling the maps. - x_lims : tuple[float, float], optional - Range for the x-axis as (min, max). If None, the range is determined from the - data, defaults to None. - y_lims : tuple[float, float], optional - Range for the y-axis as (min, max). If None, the range is determined from the - data, defaults to None. - z_lims : tuple[float, float], optional - Range for the z-axis as (min, max). If None, the range is determined from the - data, defaults to None. - raw_data : MagneticModel, optional - Flux and current information for comparison.. - latex : bool, optional - Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX - installation, defaults to False. - - """ - if latex: - set_latex_style() - width = plt.rcParams["figure.figsize"][0] * 1.4 - height = plt.rcParams["figure.figsize"][1] * 1.4 - size = (width, height) - plt.rcParams.update({"savefig.pad_inches": 0.3}) - else: - set_screen_style() - size = plt.rcParams["figure.figsize"] - - # Check if the base values were given - pu_vals = base is not None - if base is None: - base = BaseValues.unity() - - # Normalize the data - i_s_dq = data.i_s_dq / base.i - psi_s_dq = data.psi_s_dq / base.psi - - fig1, ax1 = plt.subplots(figsize=size, subplot_kw={"projection": "3d"}) - fig2, ax2 = plt.subplots(figsize=size, subplot_kw={"projection": "3d"}) - - if data.type == "flux_map": - x, y, z1, z2 = i_s_dq.real, i_s_dq.imag, psi_s_dq.real, psi_s_dq.imag - - if raw_data is not None: - raw_i_s_dq = raw_data.i_s_dq / base.i - raw_psi_s_dq = raw_data.psi_s_dq / base.psi - ax1.scatter( - raw_i_s_dq.real, - raw_i_s_dq.imag, - raw_psi_s_dq.real, - marker=".", - color="r", - ) - ax2.scatter( - raw_i_s_dq.real, - raw_i_s_dq.imag, - raw_psi_s_dq.imag, - marker=".", - color="r", - ) - - if pu_vals: - xlabel = r"$i_\mathrm{d}$ (p.u.)" - ylabel = r"$i_\mathrm{q}$ (p.u.)" - zlabel1 = r"$\psi_\mathrm{d}$ (p.u.)" - zlabel2 = r"$\psi_\mathrm{q}$ (p.u.)" - else: - xlabel = r"$i_\mathrm{d}$ (A)" - ylabel = r"$i_\mathrm{q}$ (A)" - zlabel1 = r"$\psi_\mathrm{d}$ (Vs)" - zlabel2 = r"$\psi_\mathrm{q}$ (Vs)" - - else: - x, y, z1, z2 = psi_s_dq.real, psi_s_dq.imag, i_s_dq.real, i_s_dq.imag - - if raw_data is not None: - raw_i_s_dq = raw_data.i_s_dq / base.i - raw_psi_s_dq = raw_data.psi_s_dq / base.psi - ax1.scatter( - raw_psi_s_dq.real, - raw_psi_s_dq.imag, - raw_i_s_dq.real, - marker=".", - color="r", - ) - ax2.scatter( - raw_psi_s_dq.real, - raw_psi_s_dq.imag, - raw_i_s_dq.imag, - marker=".", - color="r", - ) - - if pu_vals: - xlabel = r"$\psi_\mathrm{d}$ (p.u.)" - ylabel = r"$\psi_\mathrm{q}$ (p.u.)" - zlabel1 = r"$i_\mathrm{d}$ (p.u.)" - zlabel2 = r"$i_\mathrm{q}$ (p.u.)" - else: - xlabel = r"$\psi_\mathrm{d}$ (Vs)" - ylabel = r"$\psi_\mathrm{q}$ (Vs)" - zlabel1 = r"$i_\mathrm{d}$ (A)" - zlabel2 = r"$i_\mathrm{q}$ (A)" - - ax1.plot_surface( # type: ignore - x, y, z1, cmap="viridis", alpha=0.75, axlim_clip=True - ) - ax2.plot_surface( # type: ignore - x, y, z2, cmap="viridis", alpha=0.75, axlim_clip=True - ) - - ax1.set_xlabel(xlabel) - ax1.set_ylabel(ylabel) - ax1.set_zlabel(zlabel1) # type: ignore - ax2.set_xlabel(xlabel) - ax2.set_ylabel(ylabel) - ax2.set_zlabel(zlabel2) # type: ignore - - if x_lims is not None: - ax1.set_xlim(x_lims) - ax2.set_xlim(x_lims) - if y_lims is not None: - ax1.set_ylim(y_lims) - ax2.set_ylim(y_lims) - if z_lims is not None: - ax1.set_zlim(z_lims) # type: ignore - ax2.set_zlim(z_lims) # type: ignore - - plt.show() - - -# %% -def plot_flux_vs_current( - data: MagneticModel, - base: BaseValues | None = None, - lims: tuple[float, float] | None = None, - latex: bool = False, -) -> None: - """ - Plot the flux vs. current characteristics. - - Parameters - ---------- - data : MagneticModel - Flux map data. The current array should be a rectilinear grid. - base : BaseValues, optional - Base values for scaling the maps. - lims : tuple[float, float], optional - Range for the x-axis as (min, max). If None, determined from the data. - latex : bool, optional - Use LaTeX fonts for the labels, requires a working LaTeX installation. - - """ - if latex: - set_latex_style() - else: - set_screen_style() - - # Check if the base values were given - pu_vals = base is not None - if base is None: - base = BaseValues.unity() - - # Normalize the data - i_s_dq = data.i_s_dq / base.i - psi_s_dq = data.psi_s_dq / base.psi - - # Indices corresponding to i_d = 0 and i_q = 0 - ind_d_0 = np.argmin(np.abs(i_s_dq.real[:, 0])) - ind_q_0 = np.argmin(np.abs(i_s_dq.imag[0, :])) - # Indices corresponding to min(i_d) and max(i_d) - if lims is None: - ind_d_min = np.argmin(i_s_dq.real[:, 0]) - ind_d_max = np.argmax(i_s_dq.real[:, 0]) - else: - ind_d_min = np.argmin(np.abs(i_s_dq.real[:, 0] - lims[0])) - ind_d_max = np.argmin(np.abs(i_s_dq.real[:, 0] - lims[1])) - - fig, ax = plt.subplots() - ax.plot( - i_s_dq.real[:, ind_q_0], - psi_s_dq.real[:, ind_q_0], - color="r", - linestyle="-", - label=r"$\psi_\mathrm{d}(i_\mathrm{d}, 0)$", - ) - ax.plot( - i_s_dq.real[:, -1], - psi_s_dq.real[:, -1], - color="r", - linestyle="--", - label=r"$\psi_\mathrm{d}(i_\mathrm{d}, i_\mathrm{q,max})$", - ) - ax.plot( - i_s_dq.imag[ind_d_0, :], - psi_s_dq.imag[ind_d_0, :], - color="b", - linestyle="-", - label=r"$\psi_\mathrm{q}(0, i_\mathrm{q})$", - ) - ax.plot( - i_s_dq.imag[ind_d_min, :], - psi_s_dq.imag[ind_d_min, :], - color="b", - linestyle=":", - label=r"$\psi_\mathrm{q}(i_\mathrm{d,min}, i_\mathrm{q})$", - ) - ax.plot( - i_s_dq.imag[ind_d_max, :], - psi_s_dq.imag[ind_d_max, :], - color="b", - linestyle="--", - label=r"$\psi_\mathrm{q}(i_\mathrm{d,max}, i_\mathrm{q})$", - ) - - if lims is not None: - ax.set_xlim(lims) - - if pu_vals: - ax.set_xlabel(r"$i_\mathrm{d}$, $i_\mathrm{q}$ (p.u.)") - ax.set_ylabel(r"$\psi_\mathrm{d}$, $\psi_\mathrm{q}$ (p.u.)") - else: - ax.set_xlabel(r"$i_\mathrm{d}$, $i_\mathrm{q}$ (A)") - ax.set_ylabel(r"$\psi_\mathrm{d}$, $\psi_\mathrm{q}$ (Vs)") - ax.legend() - - plt.show() - - # %% def import_syre_data(fname: str, add_negative_q_axis: bool = True) -> MagneticModel: """ diff --git a/motulator/drive/utils/_sm_plot_control_loci.py b/motulator/drive/utils/_sm_plot_control_loci.py index d281f8c3a..08a35c54f 100644 --- a/motulator/drive/utils/_sm_plot_control_loci.py +++ b/motulator/drive/utils/_sm_plot_control_loci.py @@ -6,14 +6,24 @@ import numpy as np from cycler import cycler -from motulator.common.utils._utils import BaseValues, set_latex_style, set_screen_style +from motulator.common.utils._plotting import ( + save_and_show, + set_latex_style, + set_screen_style, +) +from motulator.common.utils._utils import BaseValues from motulator.drive.utils._parameters import ( SaturatedSynchronousMachinePars, SynchronousMachinePars, ) from motulator.drive.utils._sm_control_loci import ControlLoci +# LaTeX symbol definition for easier configuration. When per-unit values are used, +# sometimes subscript m is preferred instead of M. +M = r"\mathrm{M}" + +# %% def _setup_plot( i_s_vals: Any, latex: bool, base: BaseValues | None = None ) -> tuple[np.ndarray, bool, BaseValues]: @@ -59,6 +69,8 @@ def plot_flux_loci( base: BaseValues | None = None, num: int = 16, latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, ) -> None: """ Plot the flux linkage loci. @@ -74,6 +86,10 @@ def plot_flux_loci( Amount of points to be evaluated, defaults to 16. latex : bool, optional Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + If provided, save the figure to this path. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ i_s_vals, pu_vals, base = _setup_plot(i_s_vals, latex, base) @@ -147,7 +163,7 @@ def plot_flux_loci( ax.set_xlabel(r"$\psi_\mathrm{d}$ (Vs)") ax.set_ylabel(r"$\psi_\mathrm{q}$ (Vs)") - plt.show() + save_and_show(save_path, **savefig_kwargs) # %% def plot_current_loci( @@ -156,6 +172,8 @@ def plot_current_loci( base: BaseValues | None = None, num: int = 16, latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, ) -> None: """ Plot the current loci. @@ -171,6 +189,10 @@ def plot_current_loci( Amount of points to be evaluated, defaults to 16. latex : bool, optional Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + If provided, save the figure to this path. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ i_s_vals, pu_vals, base = _setup_plot(i_s_vals, latex, base) @@ -238,7 +260,7 @@ def plot_current_loci( ax.set_aspect("equal") - plt.show() + save_and_show(save_path, **savefig_kwargs) # %% def plot_flux_vs_torque( @@ -247,6 +269,8 @@ def plot_flux_vs_torque( base: BaseValues | None = None, num: int = 16, latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, ) -> None: """ Plot flux magnitude vs. torque characteristics. @@ -262,6 +286,10 @@ def plot_flux_vs_torque( Amount of points to be evaluated, defaults to 16. latex : bool, optional Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + If provided, save the figure to this path. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ i_s_vals, pu_vals, base = _setup_plot(i_s_vals, latex, base) @@ -315,16 +343,16 @@ def plot_flux_vs_torque( ax.legend() if pu_vals: - ax.set_xlabel(r"$\tau_\mathrm{M}$ (p.u.)") + ax.set_xlabel(rf"$\tau_{M}$ (p.u.)") ax.set_ylabel(r"$\psi_\mathrm{s}$ (p.u.)") else: - ax.set_xlabel(r"$\tau_\mathrm{M}$ (Nm)") + ax.set_xlabel(rf"$\tau_{M}$ (Nm)") ax.set_ylabel(r"$\psi_\mathrm{s}$ (Vs)") ax.set_xlim(0, mtpa.tau_M[-1] / base.tau) ax.set_ylim(0, abs(mtpa.psi_s_dq[-1]) / base.psi) - plt.show() + save_and_show(save_path, **savefig_kwargs) # %% def plot_current_vs_torque( @@ -333,6 +361,8 @@ def plot_current_vs_torque( base: BaseValues | None = None, num: int = 16, latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, ) -> None: """ Plot current vs. torque characteristics. @@ -348,6 +378,10 @@ def plot_current_vs_torque( Amount of points to be evaluated, defaults to 16. latex : bool, optional Use LaTeX fonts for the labels, requires a working LaTeX installation. + save_path : str, optional + If provided, save the figure to this path. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ i_s_vals, pu_vals, base = _setup_plot(i_s_vals, latex, base) @@ -380,8 +414,8 @@ def plot_current_vs_torque( # MTPV locus mtpv_i_s_dq_max = self._loci.compute_mtpv_current(i_s_vals[-1]) if not np.isnan(mtpv_i_s_dq_max): # MTPV exists within the current range - max_mtpv_psi_s = float(abs(self.par.psi_s_dq(mtpv_i_s_dq_max))) - mtpv = self._loci.compute_mtpv_locus(max_mtpv_psi_s, num) + mtpv_psi_s_max = float(abs(self.par.psi_s_dq(mtpv_i_s_dq_max))) + mtpv = self._loci.compute_mtpv_locus(mtpv_psi_s_max, num) ax1.plot( mtpv.tau_M / base.tau, @@ -447,10 +481,10 @@ def plot_current_vs_torque( if pu_vals: ax1.set_ylabel(r"$i_\mathrm{d}$ (p.u.)") ax2.set_ylabel(r"$i_\mathrm{q}$ (p.u.)") - ax2.set_xlabel(r"$\tau_\mathrm{M}$ (p.u.)") + ax2.set_xlabel(rf"$\tau_{M}$ (p.u.)") else: ax1.set_ylabel(r"$i_\mathrm{d}$ (A)") ax2.set_ylabel(r"$i_\mathrm{q}$ (A)") - ax2.set_xlabel(r"$\tau_\mathrm{M}$ (Nm)") + ax2.set_xlabel(rf"$\tau_{M}$ (Nm)") - plt.show() + save_and_show(save_path, **savefig_kwargs) diff --git a/motulator/drive/utils/_sm_plot_flux_maps.py b/motulator/drive/utils/_sm_plot_flux_maps.py new file mode 100644 index 000000000..9e4201b70 --- /dev/null +++ b/motulator/drive/utils/_sm_plot_flux_maps.py @@ -0,0 +1,347 @@ +"""Visualize flux linkage and current maps of synchronous machines.""" + +from typing import Literal + +import matplotlib.pyplot as plt +import numpy as np +from numpy.typing import ArrayLike + +from motulator.common.utils._plotting import ( + save_and_show, + set_latex_style, + set_screen_style, +) +from motulator.common.utils._utils import BaseValues +from motulator.drive.utils._sm_flux_maps import MagneticModel + + +# %% +def _setup_plot_style(latex: bool) -> tuple[tuple[float, float], bool]: + """Setup plot style and return figure size and pu_vals flag.""" + if latex: + set_latex_style() + width = plt.rcParams["figure.figsize"][0] * 1.4 + height = plt.rcParams["figure.figsize"][1] * 1.4 + size = (width, height) + plt.rcParams.update({"savefig.pad_inches": 0.3}) + else: + set_screen_style() + size = plt.rcParams["figure.figsize"] + return size + + +# %% +def _filter_and_plot_raw_data( + ax, + raw_data: MagneticModel, + base: BaseValues, + data_type: str, + component: str, + x_lims: tuple[float, float], + y_lims: tuple[float, float], +) -> None: + """Filter and plot raw data as scatter points.""" + raw_i_s_dq = raw_data.i_s_dq / base.i + raw_psi_s_dq = raw_data.psi_s_dq / base.psi + + if data_type == "flux_map": + mask = ( + (raw_i_s_dq.real >= x_lims[0]) + & (raw_i_s_dq.real <= x_lims[1]) + & (raw_i_s_dq.imag >= y_lims[0]) + & (raw_i_s_dq.imag <= y_lims[1]) + ) + x_scatter, y_scatter = raw_i_s_dq[mask].real, raw_i_s_dq[mask].imag + z_scatter = ( + raw_psi_s_dq[mask].real if component == "d" else raw_psi_s_dq[mask].imag + ) + else: + mask = ( + (raw_psi_s_dq.real >= x_lims[0]) + & (raw_psi_s_dq.real <= x_lims[1]) + & (raw_psi_s_dq.imag >= y_lims[0]) + & (raw_psi_s_dq.imag <= y_lims[1]) + ) + x_scatter, y_scatter = raw_psi_s_dq[mask].real, raw_psi_s_dq[mask].imag + z_scatter = raw_i_s_dq[mask].real if component == "d" else raw_i_s_dq[mask].imag + + ax.scatter(x_scatter, y_scatter, z_scatter, marker=".", color="r") + + +# %% +def _get_labels(data_type: str, component: str, pu_vals: bool) -> tuple[str, str, str]: + """Get axis labels based on data type, component, and units.""" + if data_type == "flux_map": + if pu_vals: + xlabel = r"$i_\mathrm{d}$ (p.u.)" + ylabel = r"$i_\mathrm{q}$ (p.u.)" + zlabel = rf"$\psi_\mathrm{{{component}}}$ (p.u.)" + else: + xlabel = r"$i_\mathrm{d}$ (A)" + ylabel = r"$i_\mathrm{q}$ (A)" + zlabel = rf"$\psi_\mathrm{{{component}}}$ (Vs)" + elif pu_vals: + xlabel = r"$\psi_\mathrm{d}$ (p.u.)" + ylabel = r"$\psi_\mathrm{q}$ (p.u.)" + zlabel = rf"$i_\mathrm{{{component}}}$ (p.u.)" + else: + xlabel = r"$\psi_\mathrm{d}$ (Vs)" + ylabel = r"$\psi_\mathrm{q}$ (Vs)" + zlabel = rf"$i_\mathrm{{{component}}}$ (A)" + return xlabel, ylabel, zlabel + + +def _configure_3d_axes( + ax, + x_lims: tuple[float, float], + y_lims: tuple[float, float], + z_lims: tuple[float, float] | None, + x_ticks: ArrayLike | None, + y_ticks: ArrayLike | None, + z_ticks: ArrayLike | None, +) -> None: + """Configure 3D axis limits and ticks.""" + # Set axis limits + ax.set_xlim(x_lims) + ax.set_ylim(y_lims) + if z_lims is not None: + ax.set_zlim(z_lims) # type: ignore + + # Set axis ticks + if x_ticks is not None: + ax.set_xticks(x_ticks) + if y_ticks is not None: + ax.set_yticks(y_ticks) + if z_ticks is not None: + ax.set_zticks(z_ticks) # type: ignore + + +# %% +def plot_map( + data: MagneticModel, + component: Literal["d", "q"], + base: BaseValues | None = None, + x_lims: tuple[float, float] | None = None, + x_ticks: ArrayLike | None = None, + y_lims: tuple[float, float] | None = None, + y_ticks: ArrayLike | None = None, + z_lims: tuple[float, float] | None = None, + z_ticks: ArrayLike | None = None, + raw_data: MagneticModel | None = None, + latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot component (d or q) of flux linkage or current maps. + + Parameters + ---------- + data : MagneticModel + Data containing the flux and current information. + component : {"d", "q"} + Which component to plot: "d" for d-axis, "q" for q-axis. + base : BaseValues, optional + Base values for scaling the maps. If not given, the maps are plotted + in SI units. + x_lims : tuple[float, float], optional + x-axis limits. If None, uses automatic scaling. + x_ticks : ArrayLike, optional + x-axis tick locations. + y_lims : tuple[float, float], optional + y-axis limits. If None, uses automatic scaling. + y_ticks : ArrayLike, optional + y-axis tick locations. + z_lims : tuple[float, float], optional + z-axis limits. If None, uses automatic scaling. + z_ticks : ArrayLike, optional + z-axis tick locations. + raw_data : MagneticModel, optional + Raw data for comparison (shown as scatter points). + latex : bool, optional + Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX + installation, defaults to False. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + # Setup style and base values + size = _setup_plot_style(latex) + pu_vals = base is not None + if base is None: + base = BaseValues.unity() + + # Normalize data + i_s_dq = data.i_s_dq / base.i + psi_s_dq = data.psi_s_dq / base.psi + + # If no limits are provided, determine them from the data + if x_lims is None: + if data.type == "flux_map": + x_lims = (np.min(i_s_dq.real), np.max(i_s_dq.real)) + else: + x_lims = (np.min(psi_s_dq.real), np.max(psi_s_dq.real)) + + if y_lims is None: + if data.type == "flux_map": + y_lims = (np.min(i_s_dq.imag), np.max(i_s_dq.imag)) + else: + y_lims = (np.min(psi_s_dq.imag), np.max(psi_s_dq.imag)) + + # Create figure and axis + fig, ax = plt.subplots(figsize=size, subplot_kw={"projection": "3d"}) + + # Select plot data based on map type + if data.type == "flux_map": + x, y = i_s_dq.real, i_s_dq.imag + z = psi_s_dq.real if component == "d" else psi_s_dq.imag + else: + x, y = psi_s_dq.real, psi_s_dq.imag + z = i_s_dq.real if component == "d" else i_s_dq.imag + + # Plot raw data if provided + if raw_data is not None: + _filter_and_plot_raw_data( + ax, raw_data, base, data.type, component, x_lims, y_lims + ) + + # Plot surface and set labels + ax.plot_surface( # type: ignore + x, y, z, cmap="viridis", alpha=0.75, axlim_clip=True + ) + + xlabel, ylabel, zlabel = _get_labels(data.type, component, pu_vals) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.set_zlabel(zlabel) # type: ignore + + _configure_3d_axes(ax, x_lims, y_lims, z_lims, x_ticks, y_ticks, z_ticks) + + save_and_show(save_path, **savefig_kwargs) + + +# %% +def plot_flux_vs_current( + data: MagneticModel, + base: BaseValues | None = None, + x_lims: tuple[float, float] | None = None, + x_ticks: ArrayLike | None = None, + y_lims: tuple[float, float] | None = None, + y_ticks: ArrayLike | None = None, + latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot the flux vs. current characteristics. + + Parameters + ---------- + data : MagneticModel + Flux map data. The current array should be a rectilinear grid. + base : BaseValues, optional + Base values for scaling the maps. If not given, the maps are plotted + in SI units. + x_lims : tuple[float, float], optional + x-axis limits. If None, uses automatic scaling. + x_ticks : ArrayLike, optional + x-axis tick locations. + y_lims : tuple[float, float], optional + y-axis limits. If None, uses automatic scaling. + y_ticks : ArrayLike, optional + y-axis tick locations. + latex : bool, optional + Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX + installation, defaults to False. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + if latex: + set_latex_style() + else: + set_screen_style() + + # Check if the base values were given + pu_vals = base is not None + if base is None: + base = BaseValues.unity() + + # Normalize the data + i_s_dq = data.i_s_dq / base.i + psi_s_dq = data.psi_s_dq / base.psi + + # Indices corresponding to i_d = 0 and i_q = 0 + ind_d_0 = np.argmin(np.abs(i_s_dq.real[:, 0])) + ind_q_0 = np.argmin(np.abs(i_s_dq.imag[0, :])) + + # Indices corresponding to min(i_d) and max(i_d) + if x_lims is None: + ind_d_min = np.argmin(i_s_dq.real[:, 0]) + ind_d_max = np.argmax(i_s_dq.real[:, 0]) + else: + ind_d_min = np.argmin(np.abs(i_s_dq.real[:, 0] - x_lims[0])) + ind_d_max = np.argmin(np.abs(i_s_dq.real[:, 0] - x_lims[1])) + + fig, ax = plt.subplots() + ax.plot( + i_s_dq.real[:, ind_q_0], + psi_s_dq.real[:, ind_q_0], + color="r", + linestyle="-", + label=r"$\psi_\mathrm{d}(i_\mathrm{d}, 0)$", + ) + ax.plot( + i_s_dq.real[:, -1], + psi_s_dq.real[:, -1], + color="r", + linestyle="--", + label=r"$\psi_\mathrm{d}(i_\mathrm{d}, i_\mathrm{q}^\mathrm{max})$", + ) + ax.plot( + i_s_dq.imag[ind_d_0, :], + psi_s_dq.imag[ind_d_0, :], + color="b", + linestyle="-", + label=r"$\psi_\mathrm{q}(0, i_\mathrm{q})$", + ) + ax.plot( + i_s_dq.imag[ind_d_min, :], + psi_s_dq.imag[ind_d_min, :], + color="b", + linestyle=":", + label=r"$\psi_\mathrm{q}(i_\mathrm{d}^\mathrm{min}, i_\mathrm{q})$", + ) + ax.plot( + i_s_dq.imag[ind_d_max, :], + psi_s_dq.imag[ind_d_max, :], + color="b", + linestyle="--", + label=r"$\psi_\mathrm{q}(i_\mathrm{d}^\mathrm{max}, i_\mathrm{q})$", + ) + + # Set axis limits and ticks + if x_lims is not None: + ax.set_xlim(x_lims) + if x_ticks is not None: + ax.set_xticks(x_ticks) + + if y_lims is not None: + ax.set_ylim(y_lims) + if y_ticks is not None: + ax.set_yticks(y_ticks) + + # Set axis labels + if pu_vals: + ax.set_xlabel(r"$i_\mathrm{d}$, $i_\mathrm{q}$ (p.u.)") + ax.set_ylabel(r"$\psi_\mathrm{d}$, $\psi_\mathrm{q}$ (p.u.)") + else: + ax.set_xlabel(r"$i_\mathrm{d}$, $i_\mathrm{q}$ (A)") + ax.set_ylabel(r"$\psi_\mathrm{d}$, $\psi_\mathrm{q}$ (Vs)") + ax.legend() + + save_and_show(save_path, **savefig_kwargs) diff --git a/motulator/grid/control/_base.py b/motulator/grid/control/_base.py index 6e82585f2..3bf1b5623 100644 --- a/motulator/grid/control/_base.py +++ b/motulator/grid/control/_base.py @@ -4,7 +4,7 @@ from typing import Callable, Protocol, Sequence from motulator.common.control._base import ControlSystem, TimeSeries -from motulator.common.utils import abc2complex, get_value +from motulator.common.utils._utils import abc2complex, get_value from motulator.grid.control._controllers import DCBusVoltageController from motulator.grid.model import GridConverterSystem diff --git a/motulator/grid/control/_gfl_current_vector.py b/motulator/grid/control/_gfl_current_vector.py index 088d430d9..f9e3f1697 100644 --- a/motulator/grid/control/_gfl_current_vector.py +++ b/motulator/grid/control/_gfl_current_vector.py @@ -8,7 +8,7 @@ from motulator.common.control._base import TimeSeries from motulator.common.control._controllers import ComplexPIController from motulator.common.control._pwm import PWM -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.grid.control._base import Measurements from motulator.grid.control._controllers import CurrentLimiter diff --git a/motulator/grid/control/_gfm_observer.py b/motulator/grid/control/_gfm_observer.py index ba8eb917d..e07e90238 100644 --- a/motulator/grid/control/_gfm_observer.py +++ b/motulator/grid/control/_gfm_observer.py @@ -9,7 +9,7 @@ from motulator.common.control._base import TimeSeries from motulator.common.control._pwm import PWM -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.grid.control._base import Measurements from motulator.grid.control._controllers import CurrentLimiter diff --git a/motulator/grid/control/_gfm_psc.py b/motulator/grid/control/_gfm_psc.py index d912f8822..bd52be15e 100644 --- a/motulator/grid/control/_gfm_psc.py +++ b/motulator/grid/control/_gfm_psc.py @@ -7,7 +7,7 @@ from motulator.common.control._base import TimeSeries from motulator.common.control._pwm import PWM -from motulator.common.utils import wrap +from motulator.common.utils._utils import wrap from motulator.grid.control._base import Measurements from motulator.grid.control._controllers import CurrentLimiter diff --git a/motulator/grid/model/_ac_filter.py b/motulator/grid/model/_ac_filter.py index db4c3b51c..e6115a33c 100644 --- a/motulator/grid/model/_ac_filter.py +++ b/motulator/grid/model/_ac_filter.py @@ -13,7 +13,7 @@ import numpy as np from motulator.common.model import Subsystem, SubsystemTimeSeries -from motulator.common.utils import complex2abc, empty_array +from motulator.common.utils._utils import complex2abc, empty_array # %% diff --git a/motulator/grid/model/_ac_source.py b/motulator/grid/model/_ac_source.py index d381ddc08..9da7b11ea 100644 --- a/motulator/grid/model/_ac_source.py +++ b/motulator/grid/model/_ac_source.py @@ -7,7 +7,7 @@ import numpy as np from motulator.common.model import Subsystem, SubsystemTimeSeries -from motulator.common.utils import empty_array, get_value +from motulator.common.utils._utils import empty_array, get_value @dataclass diff --git a/motulator/grid/utils/__init__.py b/motulator/grid/utils/__init__.py index 16785c295..865f410d7 100644 --- a/motulator/grid/utils/__init__.py +++ b/motulator/grid/utils/__init__.py @@ -6,12 +6,17 @@ SequenceGenerator, Step, ) -from motulator.grid.utils._plots import plot, plot_voltage_vector +from motulator.grid.utils._plots import ( + plot_control_signals, + plot_grid_waveforms, + plot_voltage_vector, +) __all__ = [ "BaseValues", "NominalValues", - "plot", + "plot_control_signals", + "plot_grid_waveforms", "plot_voltage_vector", "Step", "SequenceGenerator", diff --git a/motulator/grid/utils/_plots.py b/motulator/grid/utils/_plots.py index 18fbc58eb..f0875f05b 100644 --- a/motulator/grid/utils/_plots.py +++ b/motulator/grid/utils/_plots.py @@ -2,42 +2,21 @@ import matplotlib.pyplot as plt import numpy as np +from numpy.typing import ArrayLike from motulator.common.model._simulation import SimulationResults -from motulator.common.utils._utils import ( - BaseValues, - complex2abc, +from motulator.common.utils._plotting import ( + configure_axes, + save_and_show, set_latex_style, set_screen_style, ) +from motulator.common.utils._utils import BaseValues, complex2abc # %% -def plot( - res: SimulationResults, - base: BaseValues | None, - t_span: tuple[float, float] | None = None, - latex: bool = False, - plot_pcc_voltage: bool = True, -) -> None: - """ - Plot example figures. - - Parameters - ---------- - res : SimulationResults - Should contain the simulated data. - base : BaseValues, optional - Base values for scaling the waveforms. If not given, the waveforms are plotted - in SI units. - t_span : 2-tuple, optional - Time span. If not given, the whole simulation time is plotted. - plot_pcc_voltage : bool, optional - If True, plot the phase voltage waveforms at the point of common coupling (PCC). - Otherwise, plot the grid voltage waveforms, defaults to True. - - """ - # ruff: noqa: PLR0912, PLR0915, PLR0913 +def _setup_plot(latex: bool) -> tuple[float, float]: + """Setup plot style and return figure dimensions.""" if latex: set_latex_style() height = plt.rcParams["figure.figsize"][1] * 2.2 @@ -46,128 +25,184 @@ def plot( height = plt.rcParams["figure.figsize"][1] * 1.8 width = plt.rcParams["figure.figsize"][0] + return width, height - # For brevity - mdl = res.mdl # Continuous-time data - ctrl = res.ctrl # Discrete-time data - - # Check if the time span was given - if t_span is None: - t_span = (0, mdl.t[-1]) # Time span - # Check if the base values were given - if base is None: - pu_vals = False - base = BaseValues.unity() - else: - pu_vals = True - - # Three-phase quantities - i_g_abc = complex2abc(mdl.ac_filter.i_g_ab).T - e_g_abc = complex2abc(mdl.ac_source.e_g_ab).T - u_g_abc = complex2abc(mdl.ac_filter.u_g_ab).T - - # Calculation of active and reactive powers +# %% +def _plot_powers(ax, ctrl, mdl, base: BaseValues) -> None: + """Plot active and reactive powers.""" + # Calculate powers p_g = 1.5 * np.real(mdl.ac_filter.e_g_ab * np.conj(mdl.ac_filter.i_g_ab)) q_g = 1.5 * np.imag(mdl.ac_filter.e_g_ab * np.conj(mdl.ac_filter.i_g_ab)) - # Figure 1 - fig1, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(width, height)) - - # Subplot 1: Active and reactive power - ax1.plot( + ax.plot( ctrl.t, ctrl.ref.p_g / base.p, "--", - label=r"$p_\mathrm{g,ref}$", + label=r"$p_\mathrm{g}^\mathrm{ref}$", ds="steps-post", ) - ax1.plot(mdl.t, p_g / base.p, label=r"$p_\mathrm{g}$") + ax.plot(mdl.t, p_g / base.p, label=r"$p_\mathrm{g}$") + if hasattr(ctrl.ref, "q_g") and np.any(ctrl.ref.q_g): - ax1.plot( + ax.plot( ctrl.t, ctrl.ref.q_g / base.p, "--", - label=r"$q_\mathrm{g,ref}$", + label=r"$q_\mathrm{g}^\mathrm{ref}$", ds="steps-post", ) - ax1.plot(mdl.t, q_g / base.p, label=r"$q_\mathrm{g}$") - ax1.legend() - ax1.set_xlim(t_span) - ax1.set_xticklabels([]) + ax.plot(mdl.t, q_g / base.p, label=r"$q_\mathrm{g}$") + ax.legend() - # Subplot 2: Converter currents - ax2.plot( + +def _plot_currents(ax, ctrl, base: BaseValues) -> None: + """Plot converter currents.""" + ax.plot( ctrl.t, np.real(ctrl.ref.i_c / base.i), "--", - label=r"$i_\mathrm{cd,ref}$", + label=r"$i_\mathrm{cd}^\mathrm{ref}$", ds="steps-post", ) - ax2.plot( + ax.plot( ctrl.t, np.real(ctrl.fbk.i_c / base.i), label=r"$i_\mathrm{cd}$", ds="steps-post", ) - ax2.plot( + ax.plot( ctrl.t, np.imag(ctrl.ref.i_c / base.i), "--", - label=r"$i_\mathrm{cq,ref}$", + label=r"$i_\mathrm{cq}^\mathrm{ref}$", ds="steps-post", ) - ax2.plot( + ax.plot( ctrl.t, np.imag(ctrl.fbk.i_c / base.i), label=r"$i_\mathrm{cq}$", ds="steps-post", ) - ax2.legend() - ax2.set_xlim(t_span) - ax2.set_xticklabels([]) - - # Subplot 3: Converter voltage reference, quasi-static converter voltage, and grid - # voltage magnitudes - ax3.plot( - ctrl.t, np.abs(ctrl.fbk.u_c / base.u), label=r"$u_\mathrm{c}$", ds="steps-post" + ax.legend() + + +def _plot_voltages(ax, ctrl, mdl, base: BaseValues) -> None: + """Plot voltages.""" + ax.plot( + ctrl.t, + np.abs(ctrl.fbk.u_c / base.u), + label=r"$\hat{u}_\mathrm{c}$", + ds="steps-post", ) + if hasattr(ctrl.ref, "v_c"): - ax3.plot( + ax.plot( ctrl.t, np.abs(ctrl.ref.v_c / base.u), - label=r"$v_\mathrm{c,ref}$", + label=r"$v_\mathrm{c}^\mathrm{ref}$", ds="steps-post", ) + if hasattr(ctrl.fbk, "v_c"): - ax3.plot( + ax.plot( ctrl.t, np.abs(ctrl.fbk.v_c / base.u), label=r"$\hat{v}_\mathrm{c}$", ds="steps-post", ) + if np.any(ctrl.ref.u_dc): - ax3.plot( + ax.plot( ctrl.t, ctrl.ref.u_dc / np.sqrt(3) / base.u, "--", - label=r"$u_\mathrm{dc,ref}/\sqrt{3}$", + label=r"$u_\mathrm{dc}^\mathrm{ref}/\sqrt{3}$", ds="steps-post", ) - ax3.plot( + + ax.plot( ctrl.t, ctrl.fbk.u_dc / np.sqrt(3) / base.u, "--", label=r"$u_\mathrm{dc}/\sqrt{3}$", ds="steps-post", ) - ax3.plot( - mdl.t, np.abs(mdl.ac_source.e_g_ab / base.u), "--", label=r"$e_\mathrm{g}$" - ) - _, ymax = ax3.get_ylim() - ax3.set_ylim(0, ymax) - ax3.legend() - ax3.set_xlim(t_span) + ax.plot(mdl.t, np.abs(mdl.ac_source.e_g_ab / base.u), "--", label=r"$e_\mathrm{g}$") + + _, ymax = ax.get_ylim() + ax.set_ylim(0, ymax) + ax.legend() + + +def plot_control_signals( + res: SimulationResults, + base: BaseValues | None = None, + t_lims: tuple[float, float] | None = None, + t_ticks: ArrayLike | None = None, + y_lims: list[tuple[float, float] | None] | None = None, + y_ticks: list[ArrayLike | None] | None = None, + latex: bool = False, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot control signals and converter voltages. + + Parameters + ---------- + res : SimulationResults + Simulation results. + base : BaseValues, optional + Base values for scaling the waveforms. If not given, the waveforms are plotted + in SI units. + t_lims : tuple[float, float], optional + Time axis limits. If None, uses full time range. + t_ticks : ArrayLike, optional + Time axis tick locations. + y_lims : list[tuple[float, float] | None], optional + y-axis limits for each subplot. + y_ticks : list[ArrayLike | None], optional + y-axis tick locations for each subplot. + latex : bool, optional + Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX + installation, defaults to False. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + # Setup plot style and dimensions + width, height = _setup_plot(latex) + + # Process base values + if base is None: + pu_vals = False + base = BaseValues.unity() + else: + pu_vals = True + + # Set time limits + if t_lims is None: + t_lims = (0, res.mdl.t[-1]) + + # Create figure + fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(width, height)) + axes = [ax1, ax2, ax3] + + # Plot subplots + _plot_powers(ax1, res.ctrl, res.mdl, base) + _plot_currents(ax2, res.ctrl, base) + _plot_voltages(ax3, res.ctrl, res.mdl, base) + + # Configure all axes + configure_axes(axes, t_lims, t_ticks, y_lims, y_ticks) + + # Remove xticklabels for all but the last subplot + for ax in axes[:-1]: + ax.set_xticklabels([]) + ax3.set_xlabel("Time (s)") # Add axis labels if pu_vals: @@ -175,60 +210,134 @@ def plot( ax2.set_ylabel("Current (p.u.)") ax3.set_ylabel("Voltage (p.u.)") else: - ax1.set_ylabel("Power (kW, kVAr)") + ax1.set_ylabel("Power (W)") ax2.set_ylabel("Current (A)") ax3.set_ylabel("Voltage (V)") - ax3.set_xlabel("Time (s)") - - fig1.align_ylabels() + fig.align_ylabels() - plt.show() + save_and_show(save_path, **savefig_kwargs) - # %% - # Figure 2 - fig2, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(width, height)) - if not plot_pcc_voltage: - # Subplot 1: Grid voltage - ax1.plot( +# %% +def _plot_three_phase_voltages( + ax, mdl, base: BaseValues, plot_pcc_voltage: bool +) -> None: + """Plot three-phase voltages.""" + if plot_pcc_voltage: + u_g_abc = complex2abc(mdl.ac_filter.u_g_ab).T + ax.plot( mdl.t, - e_g_abc / base.u, - label=[r"$e_\mathrm{ga}$", r"$e_\mathrm{gb}$", r"$e_\mathrm{gc}$"], + u_g_abc / base.u, + label=[r"$u_\mathrm{ga}$", r"$u_\mathrm{gb}$", r"$u_\mathrm{gc}$"], ) else: - # Subplot 1: PCC voltage - ax1.plot( + e_g_abc = complex2abc(mdl.ac_source.e_g_ab).T + ax.plot( mdl.t, - u_g_abc / base.u, - label=[r"$u_\mathrm{ga}$", r"$u_\mathrm{gb}$", r"$u_\mathrm{gc}$"], + e_g_abc / base.u, + label=[r"$e_\mathrm{ga}$", r"$e_\mathrm{gb}$", r"$e_\mathrm{gc}$"], ) - ax1.legend() - ax1.set_xlim(t_span) - ax1.set_xticklabels([]) + ax.legend() - # Subplot 2: Grid currents - ax2.plot( + +def _plot_three_phase_currents(ax, mdl, base: BaseValues) -> None: + """Plot three-phase currents.""" + i_g_abc = complex2abc(mdl.ac_filter.i_g_ab).T + ax.plot( mdl.t, i_g_abc / base.i, label=[r"$i_\mathrm{ga}$", r"$i_\mathrm{gb}$", r"$i_\mathrm{gc}$"], ) - ax2.legend() - ax2.set_xlim(t_span) - ax2.set_xticklabels([]) + ax.legend() - # Subplot 3: Phase angles - ax3.plot(mdl.t, 180 / np.pi * mdl.ac_source.theta_g, label=r"$\theta_\mathrm{g}$") - ax3.plot( + +def _plot_phase_angles(ax, mdl, ctrl) -> None: + """Plot phase angles.""" + ax.plot(mdl.t, 180 / np.pi * mdl.ac_source.theta_g, label=r"$\theta_\mathrm{g}$") + ax.plot( ctrl.t, 180 / np.pi * ctrl.fbk.theta_c, "--", label=r"$\theta_\mathrm{c}$", ds="steps-post", ) - ax3.legend() - ax3.set_xlim(t_span) - ax3.set_ylim(-180, 180) - ax3.set_yticks([-180, -90, 0, 90, 180]) + ax.legend() + ax.set_ylim(-180, 180) + ax.set_yticks([-180, -90, 0, 90, 180]) + + +def plot_grid_waveforms( + res: SimulationResults, + base: BaseValues | None = None, + t_lims: tuple[float, float] | None = None, + t_ticks: ArrayLike | None = None, + y_lims: list[tuple[float, float] | None] | None = None, + y_ticks: list[ArrayLike | None] | None = None, + latex: bool = False, + plot_pcc_voltage: bool = True, + save_path: str | None = None, + **savefig_kwargs, +) -> None: + """ + Plot grid waveforms and phase angles. + + Parameters + ---------- + res : SimulationResults + Simulation results. + base : BaseValues, optional + Base values for scaling the waveforms. If not given, the waveforms are plotted + in SI units. + t_lims : tuple[float, float], optional + Time axis limits. If None, uses full time range. + t_ticks : ArrayLike, optional + Time axis tick locations. + y_lims : list[tuple[float, float] | None], optional + y-axis limits for each subplot. + y_ticks : list[ArrayLike | None], optional + y-axis tick locations for each subplot. + latex : bool, optional + Use LaTeX fonts for the labels. Enabling this option requires a working LaTeX + installation, defaults to False. + plot_pcc_voltage : bool, optional + If True, plot the phase voltage waveforms at the point of common coupling (PCC). + Otherwise, plot the grid voltage waveforms, defaults to True. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). + + """ + # Setup plot style and dimensions + width, height = _setup_plot(latex) + + # Process base values + if base is None: + pu_vals = False + base = BaseValues.unity() + else: + pu_vals = True + + # Set time limits + if t_lims is None: + t_lims = (0, res.mdl.t[-1]) + + # Create figure + fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(width, height)) + axes = [ax1, ax2, ax3] + + # Plot subplots + _plot_three_phase_voltages(ax1, res.mdl, base, plot_pcc_voltage) + _plot_three_phase_currents(ax2, res.mdl, base) + _plot_phase_angles(ax3, res.mdl, res.ctrl) + + # Configure all axes + configure_axes(axes, t_lims, t_ticks, y_lims, y_ticks) + + # Remove xticklabels for all but the last subplot + for ax in axes[:-1]: + ax.set_xticklabels([]) + ax3.set_xlabel("Time (s)") # Add axis labels if pu_vals: @@ -238,15 +347,18 @@ def plot( ax1.set_ylabel("Voltage (V)") ax2.set_ylabel("Current (A)") ax3.set_ylabel("Angle (deg)") - ax3.set_xlabel("Time (s)") + fig.align_ylabels() - fig2.align_ylabels() - - plt.show() + save_and_show(save_path, **savefig_kwargs) # %% -def plot_voltage_vector(res: SimulationResults, base: BaseValues | None) -> None: +def plot_voltage_vector( + res: SimulationResults, + base: BaseValues | None = None, + save_path: str | None = None, + **savefig_kwargs, +) -> None: """ Plot locus of the grid voltage vector. @@ -255,13 +367,17 @@ def plot_voltage_vector(res: SimulationResults, base: BaseValues | None) -> None res : SimulationResults Simulation results. base : BaseValues, optional - Base values for scaling the waveforms. + Base values for scaling the waveforms. If not given, the waveforms are plotted + in SI units. + save_path : str, optional + Path to save the figure. If None, the figure is not saved. + **savefig_kwargs + Additional keyword arguments passed to plt.savefig(). """ - mdl = res.mdl # For brevity - # Check if the base values were given + # Process base values if base is None: pu_vals = False base = BaseValues.unity() @@ -278,10 +394,10 @@ def plot_voltage_vector(res: SimulationResults, base: BaseValues | None) -> None ) ax.axhline(0, color="k") ax.axvline(0, color="k") - ticks = [-1.5, -1, -0.5, 0, 0.5, 1, 1.5] if pu_vals: ax.set_xlabel("Real (p.u.)") ax.set_ylabel("Imaginary (p.u.)") + ticks = [-1.5, -1, -0.5, 0, 0.5, 1, 1.5] ax.set_xticks(ticks) ax.set_yticks(ticks) else: @@ -290,4 +406,4 @@ def plot_voltage_vector(res: SimulationResults, base: BaseValues | None) -> None ax.legend() ax.set_aspect("equal") - plt.show() + save_and_show(save_path, **savefig_kwargs) diff --git a/pyproject.toml b/pyproject.toml index 86ea5aa8f..84a671b43 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "motulator" -version = "0.6.1" +version = "0.6.2" dependencies = [ "numpy", "scipy",