Add possibility to clear the stored simulation data from the classes

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2021 Aalto Electric Drives
+Copyright (c) 2024 Aalto Electric Drives
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

docs/source/conf.py

@@ -21,7 +21,7 @@
 # -- Project information -----------------------------------------------------
 
 project = "motulator"
-copyright = "2023, Aalto Electric Drives"
+copyright = "2024, Aalto Electric Drives"
 author = "Aalto Electric Drives"
 
 # The full version, including alpha/beta/rc tags

docs/source/control/observers.rst

@@ -79,7 +79,7 @@ Notice that the right-hand side of :eq:`hatws` is independent of :math:`\omega_\
 Gain Selection
 ^^^^^^^^^^^^^^
 
-The estimation-error dynamics are obtained by subtracting :eq:`dhatpsiR` from :eq:`dpsiR`. The resulting system can be linearized for analysis and gain selection purposes. Using the rotor speed as an exmaple, the small-signal deviation about the operating point is :math:`\Delta \omega_\mathrm{m} = \omega_\mathrm{m} - \omega_\mathrm{m0}`, where the subscript 0 refers to the operating point. Linearization of the estimation-error dynamics leads to [#Hin2010]_
+The estimation-error dynamics are obtained by subtracting :eq:`dhatpsiR` from :eq:`dpsiR`. The resulting system can be linearized for analysis and gain selection purposes. Using the rotor speed as an example, the small-signal deviation about the operating point is :math:`\Delta \omega_\mathrm{m} = \omega_\mathrm{m} - \omega_\mathrm{m0}`, where the subscript 0 refers to the operating point. Linearization of the estimation-error dynamics leads to [#Hin2010]_
 
 .. math::
 	\frac{\mathrm{d} \Delta\tilde{\boldsymbol{\psi}}_\mathrm{R}}{\mathrm{d} t} =  \boldsymbol{k}_1\Delta \tilde{\boldsymbol{v}} + \boldsymbol{k}_2\Delta \tilde{\boldsymbol{v}}^* - \mathrm{j}\omega_\mathrm{s0}\Delta\tilde{\boldsymbol{\psi}}_\mathrm{R}

docs/source/model/mechanics.rst

@@ -23,7 +23,7 @@ where :math:`\tau_{\mathrm{L},\omega}` is the speed-dependent load torque and :m
     \tau_{\mathrm{L},\omega} = B\omega_\mathrm{M}
     :label: viscous_friction
 
-where :math:`B` is the viscous damping coefficient. Fiscous friction appears, e.g., due to laminar fluid flow in bearings. Another typical component is quadratic load torque
+where :math:`B` is the viscous damping coefficient. Viscous friction appears, e.g., due to laminar fluid flow in bearings. Another typical component is quadratic load torque
 
 .. math:: 
     \tau_{\mathrm{L},\omega} = k\omega_\mathrm{M}^2\mathrm{sign}(\omega_\mathrm{M})

examples/obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py

@@ -59,7 +59,8 @@
 # on triangularization and allows to use unstructured flux maps.
 
 # Data points for creating the interpolant
-psi_s_data, i_s_data = data.psi_s.ravel(), data.i_s.ravel()
+psi_s_data = np.asarray(data.psi_s).ravel()
+i_s_data = np.asarray(data.i_s).ravel()
 
 # Create the interpolant, i_s = current_dq(psi_s.real, psi_s.imag)
 current_dq = LinearNDInterpolator(

examples/vector/plot_vector_ctrl_pmsm_2kw.py

@@ -11,7 +11,7 @@
 
 import numpy as np
 from motulator import model, control
-from motulator import BaseValues, Sequence, plot
+from motulator import BaseValues, Sequence, plot, plot_extra
 
 # %%
 # Compute base values based on the nominal values (just for figures).
@@ -54,6 +54,16 @@
 # %%
 # Create the simulation object and simulate it.
 
+# Simulate the system without modeling PWM
 sim = model.Simulation(mdl, ctrl, pwm=False)
 sim.simulate(t_stop=4)
 plot(sim, base)  # Plot results in per-unit values.
+
+# Repeat the same simulation with PWM model enabled (takes a bit longer)
+mdl.clear()  # First clear the stored data from the previous simulation run
+ctrl.clear()
+sim = model.Simulation(mdl, ctrl, pwm=True)
+sim.simulate(t_stop=4)
+plot(sim, base)
+# Plot a zoomed view
+plot_extra(sim, t_span=(1.1, 1.125), base=base)

motulator/_plots.py

@@ -176,16 +176,22 @@ def plot_extra(sim, base=None, t_span=None):
     if t_span is None:
         t_span = (0, ctrl.t[-1])
 
-    # Check if the base values were iven
+    # Check if the base values were given
     if base is not None:
         pu_vals = True
     else:
         pu_vals = False
         base = Bunch(w=1, u=1, i=1, psi=1, tau=1)  # Unity base values
 
+    # Angle of synchronous coordinates
+    try:
+        theta = ctrl.theta_s  # Induction machine
+    except AttributeError:
+        theta = ctrl.theta_m  # Synchronous machine
+
     # Quantities in stator coordinates
-    ctrl.u_ss = np.exp(1j*ctrl.theta_s)*ctrl.u_s
-    ctrl.i_ss = np.exp(1j*ctrl.theta_s)*ctrl.i_s
+    ctrl.u_ss = np.exp(1j*theta)*ctrl.u_s
+    ctrl.i_ss = np.exp(1j*theta)*ctrl.i_s
 
     fig1, (ax1, ax2) = plt.subplots(2, 1)
 

motulator/control/_common.py

@@ -207,7 +207,7 @@ class PICtrl:
 
     Notes
     -----
-    This contoller can be used, e.g., as a speed controller. In this case, `y` 
+    This controller can be used, e.g., as a speed controller. In this case, `y` 
     corresponds to the rotor angular speed `w_M` and `u` to the torque reference 
     `tau_M_ref`.
 
@@ -345,8 +345,8 @@ class ComplexPICtrl:
 
     Notes
     -----
-    This contoller can be used, e.g., as a current controller. In this case, `i`
-    corresponds to the stator current and `u` to the stator voltage.
+    This controller can be used, e.g., as a current controller. In this case, 
+    `i` corresponds to the stator current and `u` to the stator voltage.
 
     References
     ----------
@@ -489,16 +489,27 @@ class Ctrl:
     """Base class for the control system."""
 
     def __init__(self):
+        self.clear()
+
+    def clear(self):
+        """
+        Clear the internal data store of the control system.
+
+        This method is automatically run when the instance for the control 
+        system is created. It can also be used in the case of repeated 
+        simulations to clear the data from the previous simulation run.
+
+        """
         self.data = Bunch()  # Data store
         self.clock = Clock()  # Digital clock
 
     def __call__(self, mdl):
         """
         Run the main control loop.
 
-        The main control loop is callable that returns the sampling
-        period `T_s` (float)  and the duty ratios `d_abc` (ndarray, shape (3,)) 
-        for the next sampling period.
+        The main control loop is callable that returns the sampling period `T_s` 
+        (float) and the duty ratios `d_abc` (ndarray, shape (3,)) for the next 
+        sampling period.
 
         Parameters
         ----------
@@ -510,7 +521,7 @@ def __call__(self, mdl):
 
     def save(self, data):
         """
-        Save the internal date of the control system.
+        Save the internal data of the control system.
 
         Parameters
         ----------

motulator/control/im/_observers.py

@@ -1,6 +1,4 @@
-"""
-State observers for induction machines.
-"""
+"""State observers for induction machines."""
 
 import numpy as np
 
@@ -13,7 +11,8 @@ class Observer:
 
     This class implements a reduced-order flux observer for induction machines.
     Both sensored and sensorless operation are supported. The observer structure
-    is similar to [#Hin2010]_. The observer operates in estimated rotor flux coordinates. 
+    is similar to [#Hin2010]_. The observer operates in estimated rotor flux 
+    coordinates. 
 
     Parameters
     ----------
@@ -184,9 +183,9 @@ class FullOrderObserver:
 
     References
     ----------
-    .. [#Tii2023] Tiitinen, Hinkkanen, Harnefors, "Speed-Adaptive Full-Order Observer 
-        Revisited: Closed-Form Design for Induction Motor Drives," Proc. IEEE SLED, 
-        Seoul, South Korea,  Aug. 2023, https://doi.org/10.1109/SLED57582.2023.10261359
+    .. [#Tii2023] Tiitinen, Hinkkanen, Harnefors, "Speed-adaptive full-order 
+       observer revisited: Closed-form design for induction motor drives," Proc. 
+       IEEE SLED, 2023, https://doi.org/10.1109/SLED57582.2023.10261359
 
     """
 
@@ -245,7 +244,10 @@ def _f(self, u_s, i_s, w_m):
             (w_s - w_m)*self.L_sgm*i_err.imag)
 
         # Derivative of the integral term of the speed estimate
-        dw_m = -self.alpha_o*self.alpha_i*self.L_sgm*i_err.imag/self.psi_R if self.psi_R > 0 else 0
+        if self.psi_R > 0:
+            dw_m = -self.alpha_o*self.alpha_i*self.L_sgm*i_err.imag/self.psi_R
+        else:
+            dw_m = 0
 
         return di_s, dpsi_R, dw_m, w_s
 

motulator/control/im/_vector.py

@@ -8,7 +8,7 @@
 import numpy as np
 from motulator._helpers import abc2complex
 from motulator.control._common import PWM, Ctrl, ComplexPICtrl, SpeedCtrl
-from motulator.control.im._observers import Observer, FullOrderObserver
+from motulator.control.im._observers import Observer  # FullOrderObserver
 from motulator._utils import Bunch
 
 

motulator/control/sm/_observers.py

@@ -27,6 +27,8 @@ class Observer:
         Observer gain as a function of the rotor angular speed. The default is 
         ``lambda w_m: 0.25*(R_s*(L_d + L_q)/(L_d*L_q) + 0.2*abs(w_m))`` if
         `sensorless` else ``lambda w_m: 2*pi*15``.
+    sensorless : bool, optional
+        If True, sensorless control is used. The default is True.
 
     Attributes
     ----------
@@ -53,6 +55,7 @@ def __init__(self, par, alpha_o=2*np.pi*40, k=None, sensorless=True):
         self.alpha_o = alpha_o
         self.sensorless = sensorless
         self.k1 = k
+
         if self.sensorless:
             if self.k1 is None:  # If not given, use the default gains
                 sigma0 = .25*par.R_s*(par.L_d + par.L_q)/(par.L_d*par.L_q)
@@ -62,7 +65,7 @@ def __init__(self, par, alpha_o=2*np.pi*40, k=None, sensorless=True):
             if self.k1 is None:
                 self.k1 = lambda w_m: 2*np.pi*15
             self.k2 = lambda w_m: 0
-            
+
         # Initial states
         self.theta_m, self.w_m, self.psi_s = 0, 0, par.psi_f
 
@@ -99,7 +102,7 @@ def update(self, T_s, u_s, i_s, w_m=None):
             eps = -np.imag(e/psi_a) if np.abs(psi_a) > 0 else 0
             w_s = 2*self.alpha_o*eps + self.w_m
         else:
-            k1, k2 = self.k1, 0
+            k1, k2 = self.k1(w_m), 0
             w_s = w_m
             eps = 0
 
@@ -116,11 +119,12 @@ class FluxObserver:
     """
     Sensorless stator flux observer in external coordinates.
 
-    This observer estimates the stator flux linkage and the angle of the coordinate
-    system with respect to the d-axis of the rotor. Speed-estimation is omitted. 
-    The observer gain decouples the electrical and mechanical dynamics and allows placing 
-    the poles of the corresponding linearized estimation error dynamics. This implementation operates 
-    in external coordinates (typically synchronous coordinates defined by reference signals 
+    This observer estimates the stator flux linkage and the angle of the 
+    coordinate system with respect to the d-axis of the rotor. Speed-estimation 
+    is omitted. The observer gain decouples the electrical and mechanical 
+    dynamics and allows placing the poles of the corresponding linearized 
+    estimation error dynamics. This implementation operates in external 
+    coordinates (typically synchronous coordinates defined by reference signals 
     of a control system).
 
     Parameters
@@ -144,8 +148,7 @@ class FluxObserver:
     ----------
     .. [#Tii2022] Tiitinen, Hinkkanen, Kukkola, Routimo, Pellegrino, Harnefors, 
        "Stable and passive observer-based V/Hz control for synchronous Motors," 
-       Proc. IEEE ECCE, Detroit, MI, Oct. 2022, 
-       https://doi.org/10.1109/ECCE50734.2022.9947858
+       Proc. IEEE ECCE, 2022, https://doi.org/10.1109/ECCE50734.2022.9947858
 
     """
 

motulator/model/_simulation.py

@@ -136,7 +136,7 @@ def __call__(self, T_s, d_abc):
 
         # Assume falling edge and compute the normalized switching instants:
         t_n = np.append(0, np.sort(d_abc))
-        # Compute the correponding switching states:
+        # Compute the corresponding switching states:
         q_abc = (t_n[:, np.newaxis] < d_abc).astype(int)
 
         # Durations of switching states

motulator/model/im/_drive.py

@@ -35,7 +35,7 @@ class InductionMachine:
     Notes
     -----
     The Γ model is chosen here since it can be extended with the magnetic
-    saturation model in a staightforward manner. If the magnetic saturation is
+    saturation model in a straightforward manner. If the magnetic saturation is
     omitted, the Γ model is mathematically identical to the inverse-Γ and T
     models [#Sle1989]_.
 
@@ -243,7 +243,7 @@ class Drive:
         Induction machine model.
     mechanics : Mechanics
         Mechanics model.
-    converver : Inverter
+    converter : Inverter
         Inverter model.
 
     """
@@ -253,8 +253,20 @@ def __init__(self, machine=None, mechanics=None, converter=None):
         self.mechanics = mechanics
         self.converter = converter
         self.t0 = 0  # Initial time
-        # Store the solution in these lists
-        self.data = Bunch()  # Stores the solution data
+        self.clear()
+
+    def clear(self):
+        """
+        Clear the simulation data of the system model.
+
+        This method is automatically run when the instance for the system model
+        is created. It can also be used in the case of repeated simulations to 
+        clear the data from the previous simulation run.
+
+        """
+        self.t0 = 0  # Initial time
+        # Solution will be stored in the following lists
+        self.data = Bunch()
         self.data.t, self.data.q = [], []
         self.data.psi_ss, self.data.psi_rs = [], []
         self.data.theta_M, self.data.w_M = [], []
@@ -396,6 +408,11 @@ class DriveWithDiodeBridge(Drive):
     def __init__(self, machine=None, mechanics=None, converter=None):
         super().__init__(machine, mechanics, converter)
         self.converter = converter
+        self.clear()
+
+    def clear(self):
+        """Extend the base class."""
+        super().clear()
         self.data.u_dc, self.data.i_L = [], []
 
     def get_initial_values(self):
@@ -445,7 +462,7 @@ def post_process(self):
         self.data.u_g = abc2complex(u_g_abc)
         # Voltage at the output of the diode bridge
         self.data.u_di = np.amax(u_g_abc, axis=0) - np.amin(u_g_abc, axis=0)
-        # Diode briddge switching states (-1, 0, 1)
+        # Diode bridge switching states (-1, 0, 1)
         q_g_abc = ((np.amax(u_g_abc, axis=0) == u_g_abc).astype(int) -
                    (np.amin(u_g_abc, axis=0) == u_g_abc).astype(int))
         # Grid current space vector
@@ -473,6 +490,11 @@ class DriveTwoMassMechanics(Drive):
 
     def __init__(self, machine=None, mechanics=None, converter=None):
         super().__init__(machine, mechanics, converter)
+        self.clear()
+
+    def clear(self):
+        """Extend the base class."""
+        super().clear()
         self.data.w_L, self.data.theta_ML = [], []
 
     def get_initial_values(self):

motulator/model/sm/_drive.py

@@ -189,11 +189,20 @@ def __init__(self, machine=None, mechanics=None, converter=None):
         self.machine = machine
         self.mechanics = mechanics
         self.converter = converter
-
-        # Initial time
         self.t0 = 0
+        self.clear()
 
-        # Store the solution in these lists
+    def clear(self):
+        """
+        Clear the simulation data of the system model.
+
+        This method is automatically run when the instance for the system model
+        is created. It can also be used in the case of repeated simulations to 
+        clear the data from the previous simulation run.
+
+        """
+        self.t0 = 0
+        # Solution will be stored in the following lists
         self.data = Bunch()
         self.data.t, self.data.q = [], []
         self.data.psi_s, self.data.theta_m = [], []
@@ -329,6 +338,11 @@ class DriveTwoMassMechanics(Drive):
 
     def __init__(self, machine=None, mechanics=None, converter=None):
         super().__init__(machine, mechanics, converter)
+        self.clear()
+
+    def clear(self):
+        """Extend the base class."""
+        super().clear()
         self.data.w_L, self.data.theta_ML = [], []
 
     def get_initial_values(self):