From a05b8c23276b2841e67a9b13652d9761b0581fdf Mon Sep 17 00:00:00 2001 From: Marko Hinkkanen <76600872+mhinkkan@users.noreply.github.com> Date: Mon, 15 Apr 2024 07:43:46 +0300 Subject: [PATCH] Add an LC filter model, refactoring (#116) * New LC filter model, refactoring * Update _lc_filter.py --- docs/source/installation.rst | 2 +- examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py | 2 +- .../plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py | 1 + examples/vhz/plot_vhz_ctrl_6step_im_2kw.py | 13 ++- examples/vhz/plot_vhz_ctrl_im_2kw.py | 6 +- examples/vhz/plot_vhz_ctrl_im_2kw_lc.py | 92 +++++++++++++++++ motulator/_plots.py | 2 +- motulator/control/im/_obs_vhz.py | 7 +- motulator/control/im/_observers.py | 11 ++- motulator/control/im/_vhz.py | 10 +- motulator/control/sm/_flux_vector.py | 8 +- motulator/control/sm/_obs_vhz.py | 6 +- motulator/control/sm/_observers.py | 5 +- motulator/control/sm/_signal_inj.py | 6 +- motulator/control/sm/_vector.py | 5 +- motulator/model/__init__.py | 2 + motulator/model/_converter.py | 18 ++-- motulator/model/_lc_filter.py | 83 ++++++++++++++++ motulator/model/im/__init__.py | 2 + motulator/model/im/_drive.py | 89 +++++++++++++++-- motulator/model/sm/__init__.py | 2 + motulator/model/sm/_drive.py | 99 ++++++++++++++++--- 22 files changed, 400 insertions(+), 71 deletions(-) create mode 100644 examples/vhz/plot_vhz_ctrl_im_2kw_lc.py create mode 100644 motulator/model/_lc_filter.py diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 987e4e201..288a8d54a 100644 --- a/docs/source/installation.rst +++ b/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. diff --git a/examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py b/examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py index 0cb0502ec..101ba9f38 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py +++ b/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. diff --git a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py b/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py index f2b2d9fc2..8233d44f4 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py +++ b/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)") diff --git a/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py b/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py index 2ecfc4be5..463754372 100644 --- a/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py +++ b/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 diff --git a/examples/vhz/plot_vhz_ctrl_im_2kw.py b/examples/vhz/plot_vhz_ctrl_im_2kw.py index b70afe159..79e1b1f6f 100644 --- a/examples/vhz/plot_vhz_ctrl_im_2kw.py +++ b/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) diff --git a/examples/vhz/plot_vhz_ctrl_im_2kw_lc.py b/examples/vhz/plot_vhz_ctrl_im_2kw_lc.py new file mode 100644 index 000000000..bd1ff616f --- /dev/null +++ b/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 diff --git a/motulator/_plots.py b/motulator/_plots.py index 6f6ce3d09..529642d7e 100644 --- a/motulator/_plots.py +++ b/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}$"]) diff --git a/motulator/control/im/_obs_vhz.py b/motulator/control/im/_obs_vhz.py index d943ff7f0..a8434f33d 100644 --- a/motulator/control/im/_obs_vhz.py +++ b/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 diff --git a/motulator/control/im/_observers.py b/motulator/control/im/_observers.py index 15bae496f..1d7040cef 100644 --- a/motulator/control/im/_observers.py +++ b/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): """ diff --git a/motulator/control/im/_vhz.py b/motulator/control/im/_vhz.py index 24e3f9f5c..22fe6ad30 100644 --- a/motulator/control/im/_vhz.py +++ b/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 diff --git a/motulator/control/sm/_flux_vector.py b/motulator/control/sm/_flux_vector.py index 1c10cc0af..bd50c98b7 100644 --- a/motulator/control/sm/_flux_vector.py +++ b/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) diff --git a/motulator/control/sm/_obs_vhz.py b/motulator/control/sm/_obs_vhz.py index 1892194b6..8fa37683a 100644 --- a/motulator/control/sm/_obs_vhz.py +++ b/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 diff --git a/motulator/control/sm/_observers.py b/motulator/control/sm/_observers.py index a36d2ec34..84d9fe2d7 100644 --- a/motulator/control/sm/_observers.py +++ b/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) # %% diff --git a/motulator/control/sm/_signal_inj.py b/motulator/control/sm/_signal_inj.py index 41d928cda..db0f7f1ae 100644 --- a/motulator/control/sm/_signal_inj.py +++ b/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) diff --git a/motulator/control/sm/_vector.py b/motulator/control/sm/_vector.py index 3f35f3d3f..d4e46f412 100644 --- a/motulator/control/sm/_vector.py +++ b/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) diff --git a/motulator/model/__init__.py b/motulator/model/__init__.py index 14f0a9d5a..03b7e094d 100644 --- a/motulator/model/__init__.py +++ b/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", diff --git a/motulator/model/_converter.py b/motulator/model/_converter.py index 61889743d..af78817d8 100644 --- a/motulator/model/_converter.py +++ b/motulator/model/_converter.py @@ -39,15 +39,15 @@ def ac_voltage(q, u_dc): Returns ------- - u_ac : complex - AC-side voltage (V). + u_c : complex + AC-side converter voltage (V). """ - u_ac = q*u_dc - return u_ac + u_c = q*u_dc + return u_c @staticmethod - def dc_current(q, i_ac): + def dc_current(q, i_c): """ Compute the DC-side current of a lossless inverter. @@ -55,8 +55,8 @@ def dc_current(q, i_ac): ---------- q : complex Switching state vector. - i_ac : complex - AC-side current (A). + i_c : complex + AC-side converter current (A). Returns ------- @@ -64,7 +64,7 @@ def dc_current(q, i_ac): DC-side current (A). """ - i_dc = 1.5*np.real(q*np.conj(i_ac)) + i_dc = 1.5*np.real(q*np.conj(i_c)) return i_dc def meas_dc_voltage(self): @@ -111,7 +111,7 @@ def __init__(self, L, C, U_g, f_g): self.u_dc0 = np.sqrt(2)*U_g # Peak-valued line-neutral grid voltage self.u_g = np.sqrt(2/3)*U_g - # Grid angular frequeyncy + # Grid angular frequency self.w_g = 2*np.pi*f_g def grid_voltages(self, t): diff --git a/motulator/model/_lc_filter.py b/motulator/model/_lc_filter.py new file mode 100644 index 000000000..9f6038e1a --- /dev/null +++ b/motulator/model/_lc_filter.py @@ -0,0 +1,83 @@ +""" +Continuous-time model for an output LC filter. + +The space vector model is implemented in stator coordinates. + +""" +from motulator._helpers import complex2abc + + +# %% +class LCFilter: + """ + LC-filter model. + + Parameters + ---------- + L : float + Inductance (H). + C : float + Capacitance (F). + R_L : float, optional + Series resistance (Ω) of the inductor. The default is 0. + G_C : float, optional + Parallel conductance (S) of the capacitor. The default is 0. + + """ + + def __init__(self, L, C, R_L=0, G_C=0): + self.L, self.C, self.R_L, self.G_C = L, C, R_L, G_C + self.i_cs0, self.u_ss0 = 0j, 0j + + def f(self, i_cs, u_ss, u_cs, i_ss): + """ + Compute the state derivative. + + Parameters + ---------- + i_cs : complex + Converter output current (A). + u_ss : complex + Capacitor (stator) voltage (V). + u_cs : complex + Converter output voltage (V). + i_ss : complex + Stator current (A). + + Returns + ------- + complex list, length 2 + Time derivative of the state vector, [di_cs, du_ss] + + """ + + di_cs = (u_cs - u_ss - self.R_L*i_cs)/self.L + du_ss = (i_cs - i_ss - self.G_C*u_ss)/self.C + return [di_cs, du_ss] + + def meas_currents(self): + """ + Returns the converter phase currents at the end of the sampling period. + + Returns + ------- + i_c_abc : 3-tuple of floats + Phase currents (A). + + """ + i_c_abc = complex2abc(self.i_cs0) + return i_c_abc + + def meas_voltages(self): + """ + Returns the capacitor (stator) phase voltages at the end of the + sampling period. + + Returns + ------- + u_s_abc : 3-tuple of floats + Phase voltages (V). + + """ + u_s_abc = complex2abc(self.u_ss0) + return u_s_abc diff --git a/motulator/model/im/__init__.py b/motulator/model/im/__init__.py index 3c64b643e..0e7225c84 100644 --- a/motulator/model/im/__init__.py +++ b/motulator/model/im/__init__.py @@ -4,6 +4,7 @@ Drive, DriveWithDiodeBridge, DriveTwoMassMechanics, + DriveWithLCFilter, InductionMachine, InductionMachineSaturated, InductionMachineInvGamma, @@ -13,6 +14,7 @@ "Drive", "DriveWithDiodeBridge", "DriveTwoMassMechanics", + "DriveWithLCFilter", "InductionMachine", "InductionMachineSaturated", "InductionMachineInvGamma", diff --git a/motulator/model/im/_drive.py b/motulator/model/im/_drive.py index 27d2c7f1e..5ec7a9f2d 100644 --- a/motulator/model/im/_drive.py +++ b/motulator/model/im/_drive.py @@ -316,9 +316,9 @@ def f(self, t, x): # Unpack the states psi_ss, psi_rs, w_M, _ = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, self.converter.u_dc0) + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) # State derivatives plus the outputs for interconnections - machine_f, _, tau_M = self.machine.f(psi_ss, psi_rs, u_ss, w_M) + machine_f, _, tau_M = self.machine.f(psi_ss, psi_rs, u_cs, w_M) mechanics_f = self.mechanics.f(t, w_M, tau_M) # List of state derivatives return machine_f + mechanics_f @@ -343,7 +343,7 @@ def post_process(self): self.data.tau_L = ( self.mechanics.tau_L_t(self.data.t) + self.mechanics.tau_L_w(self.data.w_M)) - self.data.u_ss = self.converter.ac_voltage( + self.data.u_cs = self.converter.ac_voltage( self.data.q, self.converter.u_dc0) # Compute the inverse-Γ rotor flux @@ -405,8 +405,8 @@ def f(self, t, x): psi_ss, psi_rs, w_M, _, u_dc, i_L = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, u_dc) - machine_f, i_ss, tau_M = self.machine.f(psi_ss, psi_rs, u_ss, w_M) + u_cs = self.converter.ac_voltage(self.q, u_dc) + machine_f, i_ss, tau_M = self.machine.f(psi_ss, psi_rs, u_cs, w_M) i_dc = self.converter.dc_current(self.q, i_ss) # Return the list of state derivatives @@ -427,7 +427,7 @@ def post_process(self): self.data.u_dc = np.asarray(self.data.u_dc) self.data.i_L = np.asarray(self.data.i_L) # Some useful variables - self.data.u_ss = self.converter.ac_voltage(self.data.q, self.data.u_dc) + self.data.u_cs = self.converter.ac_voltage(self.data.q, self.data.u_dc) self.data.i_dc = self.converter.dc_current(self.data.q, self.data.i_ss) u_g_abc = self.converter.grid_voltages(self.data.t) self.data.u_g = abc2complex(u_g_abc) @@ -483,9 +483,9 @@ def f(self, t, x): # Unpack the states psi_ss, psi_rs, w_M, _, w_L, theta_ML = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, self.converter.u_dc0) + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) # State derivatives plus the outputs for interconnections - machine_f, _, tau_M = self.machine.f(psi_ss, psi_rs, u_ss, w_M) + machine_f, _, tau_M = self.machine.f(psi_ss, psi_rs, u_cs, w_M) mechanics_f = self.mechanics.f(t, w_M, w_L, theta_ML, tau_M) # List of state derivatives return machine_f + mechanics_f @@ -502,3 +502,76 @@ def post_process(self): # From lists to the ndarray self.data.w_L = np.asarray(self.data.w_L) self.data.theta_ML = np.asarray(self.data.theta_ML) + + +# %% +class DriveWithLCFilter(Drive): + """ + Induction machine drive with an output LC filter. + + Parameters + ---------- + machine : InductionMachine | InductionMachineSaturated + Induction machine model. + mechanics : Mechanics + Mechanics model. + converter : Inverter + Inverter model. + lc_filter : LCFilter + LC-filter model. + + """ + + def __init__( + self, + machine=None, + mechanics=None, + converter=None, + lc_filter=None): + super().__init__(machine, mechanics, converter) + self.lc_filter = lc_filter + self.clear() + + def clear(self): + """Extend the base class.""" + super().clear() + self.data.i_cs, self.data.u_ss = [], [] + + def get_initial_values(self): + """Extend the base class.""" + x0 = super().get_initial_values() + [ + self.lc_filter.i_cs0, self.lc_filter.u_ss0 + ] + return x0 + + def set_initial_values(self, t0, x0): + """Extend the base class.""" + super().set_initial_values(t0, x0[0:4]) + self.lc_filter.i_cs0 = x0[4] + self.lc_filter.u_ss0 = x0[5] + + def f(self, t, x): + """Override the base class.""" + # Unpack the states + psi_ss, psi_rs, w_M, _, i_cs, u_ss = x + # Interconnections: outputs for computing the state derivatives + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) + # State derivatives plus the outputs for interconnections + machine_f, i_ss, tau_M = self.machine.f(psi_ss, psi_rs, u_ss, w_M) + mechanics_f = self.mechanics.f(t, w_M, tau_M) + lc_filter_f = self.lc_filter.f(i_cs, u_ss, u_cs, i_ss) + # List of state derivatives + return machine_f + mechanics_f + lc_filter_f + + def save(self, sol): + """Extend the base class.""" + super().save(sol) + self.data.i_cs.extend(sol.y[4]) + self.data.u_ss.extend(sol.y[5]) + + def post_process(self): + """Extend the base class.""" + super().post_process() + # From lists to the ndarray + self.data.i_cs = np.asarray(self.data.i_cs) + self.data.u_ss = np.asarray(self.data.u_ss) diff --git a/motulator/model/sm/__init__.py b/motulator/model/sm/__init__.py index 94abf2bad..fbc46a6d0 100644 --- a/motulator/model/sm/__init__.py +++ b/motulator/model/sm/__init__.py @@ -4,6 +4,7 @@ Drive, DriveWithDiodeBridge, DriveTwoMassMechanics, + DriveWithLCFilter, SynchronousMachine, SynchronousMachineSaturated, ) @@ -18,6 +19,7 @@ "Drive", "DriveWithDiodeBridge", "DriveTwoMassMechanics", + "DriveWithLCFilter", "SynchronousMachine", "SynchronousMachineSaturated", "import_syre_data", diff --git a/motulator/model/sm/_drive.py b/motulator/model/sm/_drive.py index 8e2ebdcaf..e2ffdeb11 100644 --- a/motulator/model/sm/_drive.py +++ b/motulator/model/sm/_drive.py @@ -255,11 +255,12 @@ def f(self, t, x): psi_s, theta_m, w_M, _ = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, self.converter.u_dc0) - u_s = np.exp(-1j*theta_m)*u_ss # Stator voltage in rotor coordinates + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) + # Converter (stator) voltage in rotor coordinates + u_c = np.exp(-1j*theta_m)*u_cs # State derivatives - machine_f, _, tau_M = self.machine.f(psi_s, u_s, w_M) + machine_f, _, tau_M = self.machine.f(psi_s, u_c, w_M) mechanics_f = self.mechanics.f(t, w_M, tau_M) # List of state derivatives @@ -282,7 +283,7 @@ def post_process(self): self.data.tau_L = ( self.mechanics.tau_L_t(self.data.t) + self.mechanics.tau_L_w(self.data.w_M)) - self.data.u_ss = self.converter.ac_voltage( + self.data.u_cs = self.converter.ac_voltage( self.data.q, self.converter.u_dc0) self.data.theta_m = wrap(self.data.theta_m) self.data.theta_M = wrap(self.data.theta_M) @@ -338,10 +339,10 @@ def f(self, t, x): psi_s, theta_m, w_M, _, u_dc, i_L = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, u_dc) - u_s = np.exp(-1j*theta_m)*u_ss # Stator voltage in rotor coordinates + u_cs = self.converter.ac_voltage(self.q, u_dc) + u_c = np.exp(-1j*theta_m)*u_cs - machine_f, i_s, tau_M = self.machine.f(psi_s, u_s, w_M) + machine_f, i_s, tau_M = self.machine.f(psi_s, u_c, w_M) i_ss = np.exp(1j*theta_m)*i_s # Stator current in stator coordinates i_dc = self.converter.dc_current(self.q, i_ss) @@ -363,7 +364,7 @@ def post_process(self): self.data.u_dc = np.asarray(self.data.u_dc) self.data.i_L = np.asarray(self.data.i_L) # Some useful variables - self.data.u_ss = self.converter.ac_voltage(self.data.q, self.data.u_dc) + self.data.u_cs = self.converter.ac_voltage(self.data.q, self.data.u_dc) self.data.i_dc = self.converter.dc_current(self.data.q, self.data.i_ss) u_g_abc = self.converter.grid_voltages(self.data.t) self.data.u_g = abc2complex(u_g_abc) @@ -419,10 +420,10 @@ def f(self, t, x): # Unpack the states psi_s, theta_m, w_M, _, w_L, theta_ML = x # Interconnections: outputs for computing the state derivatives - u_ss = self.converter.ac_voltage(self.q, self.converter.u_dc0) - u_s = np.exp(-1j*theta_m)*u_ss # Stator voltage in rotor coordinates + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) + u_c = np.exp(-1j*theta_m)*u_cs # State derivatives plus the outputs for interconnections - machine_f, _, tau_M = self.machine.f(psi_s, u_s, w_M) + machine_f, _, tau_M = self.machine.f(psi_s, u_c, w_M) mechanics_f = self.mechanics.f(t, w_M, w_L, theta_ML, tau_M) # List of state derivatives return machine_f + mechanics_f @@ -439,3 +440,79 @@ def post_process(self): # From lists to the ndarray self.data.w_L = np.asarray(self.data.w_L) self.data.theta_ML = np.asarray(self.data.theta_ML) + + +# %% +class DriveWithLCFilter(Drive): + """ + Synchronous machine drive with an output LC filter. + + Parameters + ---------- + machine : SynchronousMachine | SynchronousMachineSaturated + Synchronous machine model. + mechanics : Mechanics + Mechanics model. + converter : Inverter + Inverter model. + lc_filter : LCFilter + LC-filter model. + + """ + + def __init__( + self, + machine=None, + mechanics=None, + converter=None, + lc_filter=None): + super().__init__(machine, mechanics, converter) + self.lc_filter = lc_filter + self.clear() + + def clear(self): + """Extend the base class.""" + super().clear() + self.data.i_cs, self.data.u_ss = [], [] + + def get_initial_values(self): + """Extend the base class.""" + x0 = super().get_initial_values() + [ + self.lc_filter.i_cs0, self.lc_filter.u_ss0 + ] + return x0 + + def set_initial_values(self, t0, x0): + """Extend the base class.""" + super().set_initial_values(t0, x0[0:4]) + self.lc_filter.i_cs0 = x0[4] + self.lc_filter.u_ss0 = x0[5] + + # pylint: disable=too-many-locals + def f(self, t, x): + """Override the base class.""" + # Unpack the states + psi_s, theta_m, w_M, _, i_cs, u_ss = x + # Interconnections: outputs for computing the state derivatives + u_cs = self.converter.ac_voltage(self.q, self.converter.u_dc0) + u_s = np.exp(-1j*theta_m)*u_ss # Stator voltage in rotor coordinates + # State derivatives plus the outputs for interconnections + machine_f, i_s, tau_M = self.machine.f(psi_s, u_s, w_M) + i_ss = np.exp(1j*theta_m)*i_s # Stator current in stator coordinates + mechanics_f = self.mechanics.f(t, w_M, tau_M) + lc_filter_f = self.lc_filter.f(i_cs, u_ss, u_cs, i_ss) + # List of state derivatives + return machine_f + mechanics_f + lc_filter_f + + def save(self, sol): + """Extend the base class.""" + super().save(sol) + self.data.i_cs.extend(sol.y[4]) + self.data.u_ss.extend(sol.y[5]) + + def post_process(self): + """Extend the base class.""" + super().post_process() + # From lists to the ndarray + self.data.i_cs = np.asarray(self.data.i_cs) + self.data.u_ss = np.asarray(self.data.u_ss)