Add an LC filter model, refactoring (#116)

* New LC filter model, refactoring

* Update _lc_filter.py

docs/source/installation.rst

@@ -32,7 +32,7 @@ Several powerful open-source IDEs are available for Python. The following instru
 1)	Install VS Code, Python, and ``git`` on your computer. Install also the recommended Python extensions in VS Code.
 2) Clone the project::
     
-      git clone https://github.com/AaltoElectricDrives/motulator
+      git clone https://github.com/Aalto-Electric-Drives/motulator
 
    This will create a folder called *motulator* in your current directory. 
 

examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py

@@ -57,7 +57,7 @@
 
 # Quadratic load torque profile, e.g. pumps and fans (uncomment to enable)
 # k = 1.1*base.tau_nom/(base.w/base.p)**2
-# mdl.mech.tau_L_w = lambda w_M: np.sign(w_M)*k*w_M**2
+# mdl.mechanics.tau_L_w = lambda w_M: np.sign(w_M)*k*w_M**2
 
 # %%
 # Create the simulation object and simulate it.

examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py

@@ -73,6 +73,7 @@
 ax1.set_xlim(t_span)
 ax2.set_xlim(t_span)
 ax1.set_xticklabels([])
+ax1.legend()
 ax1.set_ylabel(r"$\omega_\mathrm{M}$, $\omega_\mathrm{L}$ (rad/s)")
 ax2.set_ylabel(r"$\vartheta_\mathrm{ML}$ (deg)")
 ax2.set_xlabel("Time (s)")

examples/vhz/plot_vhz_ctrl_6step_im_2kw.py

@@ -2,11 +2,11 @@
 2.2-kW induction motor, 6-step mode
 ===================================
 
-This example simulates V/Hz of a 2.2-kW induction motor drive. The six-step
-overmodulation is enabled, which increases the fundamental voltage as well as
-the harmonics. Since the PWM is not synchronized with the stator frequency, the
-harmonic content also depends on the ratio between the stator frequency and the 
-sampling frequency.
+This example simulates V/Hz control of a 2.2-kW induction motor drive. The 
+six-step overmodulation is enabled, which increases the fundamental voltage as 
+well as the harmonics. Since the PWM is not synchronized with the stator 
+frequency, the harmonic content also depends on the ratio between the stator 
+frequency and the sampling frequency.
 
 """
 # %%
@@ -42,8 +42,7 @@
 ctrl.rate_limiter = control.RateLimiter(2*np.pi*120)
 
 # %%
-# Set the speed reference and the external load torque. More complicated
-# signals could be defined as functions.
+# Set the speed reference and the external load torque.
 
 # Speed reference
 times = np.array([0, .1, .3, 1])*2

examples/vhz/plot_vhz_ctrl_im_2kw.py

@@ -41,8 +41,7 @@
 ctrl.rate_limiter = control.RateLimiter(2*np.pi*120)
 
 # %%
-# Set the speed reference and the external load torque. More complicated
-# signals could be defined as functions.
+# Set the speed reference and the external load torque.
 
 ctrl.w_m_ref = lambda t: (t > .2)*base.w
 
@@ -54,8 +53,7 @@
 mdl.mechanics.tau_L_t = lambda t: (t > 1.)*.2*base.tau_nom
 
 # %%
-# Create the simulation object and simulate it. The option `pwm=True` enables
-# the model for the carrier comparison.
+# Create the simulation object and simulate it.
 
 sim = model.Simulation(mdl, ctrl)
 sim.simulate(t_stop=1.5)

examples/vhz/plot_vhz_ctrl_im_2kw_lc.py

@@ -0,0 +1,92 @@
+"""
+2.2-kW induction motor, LC filter
+=================================
+
+This example simulates open-loop V/Hz control of a 2.2-kW induction machine
+drive equipped with an LC filter. 
+
+"""
+# %%
+# Imports.
+
+import numpy as np
+import matplotlib.pyplot as plt
+from motulator import model, control
+from motulator import BaseValues, plot
+
+# %%
+# Compute base values based on the nominal values (just for figures).
+
+base = BaseValues(
+    U_nom=400, I_nom=5, f_nom=50, tau_nom=14.6, P_nom=2.2e3, n_p=2)
+
+# %%
+# Create the system model. The filter parameters correspond to [#Sal2006]_.
+
+machine = model.im.InductionMachineInvGamma(  # Inverse-Γ parameters
+    R_s=3.7, R_R=2.1, L_sgm=.021, L_M=.224, n_p=2)
+mechanics = model.Mechanics(J=.015)
+converter = model.Inverter(u_dc=540)
+lc_filter = model.LCFilter(L=8e-3, C=9.9e-6, R_L=.1)
+mdl = model.im.DriveWithLCFilter(machine, mechanics, converter, lc_filter)
+mdl.pwm = model.CarrierComparison()  # Enable the PWM model
+
+# %%
+# Control system (parametrized as open-loop V/Hz control).
+
+# Inverse-Γ model parameter estimates
+par = control.im.ModelPars(R_s=0*3.7, R_R=0*2.1, L_sgm=.021, L_M=.224)
+ctrl = control.im.VHzCtrl(250e-6, par, psi_s_nom=base.psi, k_u=0, k_w=0)
+ctrl.rate_limiter = control.RateLimiter(2*np.pi*120)
+
+# %%
+# Set the speed reference and the external load torque.
+
+ctrl.w_m_ref = lambda t: (t > .2)*base.w
+
+# Quadratic load torque profile (corresponding to pumps and fans)
+k = 1.1*base.tau_nom/(base.w/base.n_p)**2
+mdl.mechanics.tau_L_w = lambda w_M: k*w_M**2*np.sign(w_M)
+
+# %%
+# Create the simulation object and simulate it.
+
+sim = model.Simulation(mdl, ctrl)
+sim.simulate(t_stop=1.5)
+
+# %%
+# Plot results in per-unit values.
+
+# sphinx_gallery_thumbnail_number = 2
+plot(sim, base)
+
+# %%
+# Plot additional waveforms.
+
+t_span = (1.1, 1.125)  # Time span for the zoomed-in plot
+mdl = sim.mdl.data  # Continuous-time data
+# Plot the converter and stator voltages (phase a)
+fig1, (ax1, ax2) = plt.subplots(2, 1)
+ax1.plot(mdl.t, mdl.u_cs.real/base.u, label=r"$u_\mathrm{ca}$")
+ax1.plot(mdl.t, mdl.u_ss.real/base.u, label=r"$u_\mathrm{sa}$")
+ax1.set_xlim(t_span)
+ax1.legend()
+ax1.set_xticklabels([])
+ax1.set_ylabel("Voltage (p.u.)")
+# Plot he converter and stator currents (phase a)
+ax2.plot(mdl.t, mdl.i_cs.real/base.i, label=r"$i_\mathrm{ca}$")
+ax2.plot(mdl.t, mdl.i_ss.real/base.i, label=r"$i_\mathrm{sa}$")
+ax2.set_xlim(t_span)
+ax2.legend()
+ax2.set_ylabel("Current (p.u.)")
+ax2.set_xlabel("Time (s)")
+
+plt.tight_layout()
+plt.show()
+
+# %%
+# .. rubric:: References
+#
+# .. [#Sal2006] Salomäki, Hinkkanen, Luomi, "Sensorless control of induction
+#    motor drives equipped with inverter output filter," IEEE Trans. Ind.
+#    Electron., 2006, https://doi.org/10.1109/TIE.2006.878314

motulator/_plots.py

@@ -196,7 +196,7 @@ def plot_extra(sim, base=None, t_span=None):
     fig1, (ax1, ax2) = plt.subplots(2, 1)
 
     # Subplot 1: voltages
-    ax1.plot(mdl.t, mdl.u_ss.real/base.u)
+    ax1.plot(mdl.t, mdl.u_cs.real/base.u)
     ax1.plot(ctrl.t, ctrl.u_ss.real/base.u)
     ax1.set_xlim(t_span)
     ax1.legend([r"$u_\mathrm{sa}$", r"$\hat u_\mathrm{sa}$"])

motulator/control/im/_obs_vhz.py

@@ -18,7 +18,7 @@
 import numpy as np
 from motulator.control._common import Ctrl, PWM, RateLimiter
 from motulator.control.im._observers import FluxObserver
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator._utils import Bunch
 
 
@@ -110,6 +110,7 @@ def __init__(self, par, ctrl_par, T_s=250e-6):
         # States
         self.theta_s, self.tau_M_ref, self.w_r_ref = 0, 0, 0
 
+    # pylint: disable=too-many-locals
     def __call__(self, mdl):
         """
         Run the main control loop.
@@ -185,8 +186,8 @@ def __call__(self, mdl):
         self.observer.update(self.T_s, u_s, i_s, w_s)
         self.w_r_ref += self.T_s*self.alpha_r*(w_r - self.w_r_ref)
         self.tau_M_ref += self.T_s*self.alpha_f*(tau_M - self.tau_M_ref)
-        self.theta_s += self.T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_s = np.mod(self.theta_s + np.pi, 2*np.pi) - np.pi
+        self.theta_s += self.T_s*w_s
+        self.theta_s = wrap(self.theta_s)
         self.clock.update(self.T_s)
 
         # Calculate the duty ratios and update the voltage estimate state

motulator/control/im/_observers.py

@@ -1,6 +1,7 @@
 """State observers for induction machines."""
 
 import numpy as np
+from motulator._helpers import wrap
 
 
 # %%
@@ -103,8 +104,8 @@ def _f(self, T_s, u_s, i_s, w_m):
     def _update(self, T_s, i_s, dpsi_R, w_s, w_m):
         # Update the states
         self.psi_R += T_s*dpsi_R
-        self.theta_s += T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_s = np.mod(self.theta_s + np.pi, 2*np.pi) - np.pi
+        self.theta_s += T_s*w_s
+        self.theta_s = wrap(self.theta_s)
         self.w_m += T_s*self.alpha_o*(w_m - self.w_m)
         # Store the old current
         self._i_s_old = i_s
@@ -215,7 +216,7 @@ def __init__(
                     self.alpha - 1j*w_m)
         else:
             raise NotImplementedError(
-                'Sensored mode not implemented for full-order observer')
+                "Sensored mode not implemented for full-order observer")
 
         # States
         self.psi_R, self.i_s, self.theta_s, self.w_m = 0, 0, 0, 0
@@ -257,8 +258,8 @@ def _update(self, T_s, di_s, dpsi_R, dw_m, w_s):
         self.i_s += T_s*di_s
         self.psi_R += T_s*dpsi_R
         self.w_m += T_s*dw_m
-        self.theta_s += T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_s = np.mod(self.theta_s + np.pi, 2*np.pi) - np.pi
+        self.theta_s += T_s*w_s
+        self.theta_s = wrap(self.theta_s)
 
     def __call__(self, T_s, u_s, i_s, w_m=None):
         """

motulator/control/im/_vhz.py

@@ -17,7 +17,7 @@
 # %%
 import numpy as np
 from motulator.control._common import Ctrl, PWM
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator._utils import Bunch
 
 
@@ -33,7 +33,7 @@ class VHzCtrl(Ctrl):
 
     """
 
-    # pylint: disable=too-many-instance-attributes
+    # pylint: disable=too-many-instance-attributes, too-many-arguments
     def __init__(self, T_s, par, psi_s_nom, k_u=1., k_w=4., six_step=False):
         super().__init__()
         self.T_s = T_s
@@ -50,7 +50,7 @@ def __init__(self, T_s, par, psi_s_nom, k_u=1., k_w=4., six_step=False):
         # Instantiate classes
         self.pwm = PWM(six_step=six_step)
         self.rate_limiter = callable
-        self.six_step = six_step
+        # self.six_step = six_step
         # States
         self.i_s_ref, self.theta_s, self.w_r_ref = 0j, 0, 0
         # self.T_s0 = T_s
@@ -120,8 +120,8 @@ def __call__(self, mdl):
         # with ref at the end of the variable name.
         self.i_s_ref += self.T_s*self.alpha_i*(i_s - self.i_s_ref)
         self.w_r_ref += self.T_s*self.alpha_f*(w_r - self.w_r_ref)
-        self.theta_s += self.T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_s = np.mod(self.theta_s + np.pi, 2*np.pi) - np.pi
+        self.theta_s += self.T_s*w_s
+        self.theta_s = wrap(self.theta_s)
         self.clock.update(self.T_s)
 
         # Compute the duty ratios

motulator/control/sm/_flux_vector.py

@@ -2,7 +2,7 @@
 
 from dataclasses import dataclass, InitVar
 import numpy as np
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator.control._common import Ctrl, PWM, SpeedCtrl
 from motulator.control.sm._torque import TorqueCharacteristics
 from motulator.control.sm._vector import ModelPars
@@ -126,10 +126,8 @@ def __call__(self, mdl):
         else:
             # Measure the rotor speed
             w_m = self.n_p*mdl.mechanics.meas_speed()
-            # Limit the electrical rotor position into [-pi, pi)
-            theta_m = np.mod(
-                self.n_p*mdl.mechanics.meas_position() + np.pi,
-                2*np.pi) - np.pi
+            # Measure the electrical rotor position
+            theta_m = wrap(self.n_p*mdl.mechanics.meas_position())
 
         # Current vector in estimated rotor coordinates
         i_s = np.exp(-1j*theta_m)*abc2complex(i_s_abc)

motulator/control/sm/_obs_vhz.py

@@ -2,7 +2,7 @@
 
 from dataclasses import dataclass, InitVar
 import numpy as np
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator.control._common import Ctrl, PWM, RateLimiter
 from motulator.control.sm._flux_vector import (
     FluxTorqueReference, FluxTorqueReferencePars)
@@ -166,8 +166,8 @@ def __call__(self, mdl):
         # Update the states
         self.observer.update(self.T_s, u_s, i_s, w_s)
         self.tau_M_ref += self.T_s*self.alpha_f*(tau_M - self.tau_M_ref)
-        self.theta_s += self.T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_s = np.mod(self.theta_s + np.pi, 2*np.pi) - np.pi
+        self.theta_s += self.T_s*w_s
+        self.theta_s = wrap(self.theta_s)
         self.clock.update(self.T_s)
 
         # PWM output

motulator/control/sm/_observers.py

@@ -2,6 +2,7 @@
 State observers for synchronous machines.
 """
 import numpy as np
+from motulator._helpers import wrap
 
 
 # %%
@@ -110,8 +111,8 @@ def update(self, T_s, u_s, i_s, w_m=None):
         self.psi_s += T_s*(
             u_s - self.R_s*i_s - 1j*w_s*self.psi_s + k1*e + k2*np.conj(e))
         self.w_m += T_s*self.alpha_o**2*eps
-        self.theta_m += T_s*w_s  # Next line: limit into [-pi, pi)
-        self.theta_m = np.mod(self.theta_m + np.pi, 2*np.pi) - np.pi
+        self.theta_m += T_s*w_s
+        self.theta_m = wrap(self.theta_m)
 
 
 # %%

motulator/control/sm/_signal_inj.py

@@ -1,7 +1,7 @@
 """Sensorless control with signal injection for synchronous machine drives."""
 
 import numpy as np
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator.control._common import Ctrl, PWM, SpeedCtrl
 from motulator.control.sm._vector import CurrentCtrl, CurrentReference
 from motulator._utils import Bunch
@@ -261,5 +261,5 @@ def update(self, T_s, err):
 
         # Update the states
         self.w_m += T_s*self.k_i*err
-        self.theta_m += T_s*w_m  # Next line: limit into [-pi, pi)
-        self.theta_m = np.mod(self.theta_m + np.pi, 2*np.pi) - np.pi
+        self.theta_m += T_s*w_m
+        self.theta_m = wrap(self.theta_m)

motulator/control/sm/_vector.py

@@ -2,7 +2,7 @@
 
 from dataclasses import dataclass, InitVar
 import numpy as np
-from motulator._helpers import abc2complex
+from motulator._helpers import abc2complex, wrap
 from motulator.control._common import Ctrl, ComplexPICtrl, PWM, SpeedCtrl
 from motulator.control.sm._torque import TorqueCharacteristics
 from motulator.control.sm._observers import Observer
@@ -122,8 +122,7 @@ def __call__(self, mdl):
             w_m, theta_m = self.observer.w_m, self.observer.theta_m
         else:
             w_m = self.n_p*mdl.mechanics.meas_speed()
-            theta_m = self.n_p*mdl.mechanics.meas_position()
-            theta_m = np.mod(theta_m + np.pi, 2*np.pi) - np.pi
+            theta_m = wrap(self.n_p*mdl.mechanics.meas_position())
 
         # Current vector in estimated rotor coordinates
         i_s = np.exp(-1j*theta_m)*abc2complex(i_s_abc)

motulator/model/__init__.py

@@ -2,6 +2,7 @@
 
 from motulator.model._mechanics import Mechanics, MechanicsTwoMass
 from motulator.model._converter import FrequencyConverter, Inverter
+from motulator.model._lc_filter import LCFilter
 from motulator.model._simulation import (
     CarrierComparison, Simulation, Delay, zoh)
 
@@ -13,6 +14,7 @@
     "MechanicsTwoMass",
     "FrequencyConverter",
     "Inverter",
+    "LCFilter",
     "CarrierComparison",
     "Simulation",
     "Delay",