From c84f792bff118f79c365ae3f79562b0f770c65e1 Mon Sep 17 00:00:00 2001 From: Marko Hinkkanen <76600872+mhinkkan@users.noreply.github.com> Date: Thu, 11 Apr 2024 20:10:15 +0300 Subject: [PATCH] Refactoring (#115) --- .vscode/settings.json | 2 +- .../flux_vector/plot_flux_vector_pmsm_2kw.py | 2 +- .../plot_flux_vector_pmsyrm_5kw.py | 31 ++-- .../flux_vector/plot_flux_vector_syrm_7kw.py | 2 +- examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py | 6 +- .../obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw.py | 6 +- .../plot_obs_vhz_ctrl_pmsm_2kw_two_mass.py | 6 +- .../obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py | 6 +- .../obs_vhz/plot_obs_vhz_ctrl_syrm_7kw.py | 10 +- .../signal_inj/plot_signal_inj_pmsm_2kw.py | 2 +- .../signal_inj/plot_signal_inj_syrm_7kw.py | 4 +- examples/vector/plot_vector_ctrl_im_2kw.py | 21 +-- examples/vector/plot_vector_ctrl_pmsm_2kw.py | 37 ++--- .../vector/plot_vector_ctrl_pmsm_2kw_diode.py | 58 +++++++ .../vector/plot_vector_ctrl_pmsyrm_thor.py | 8 +- examples/vector/plot_vector_ctrl_syrm_7kw.py | 8 +- examples/vhz/plot_vhz_ctrl_6step_im_2kw.py | 8 +- examples/vhz/plot_vhz_ctrl_im_2kw.py | 13 +- motulator/__init__.py | 2 +- motulator/_helpers.py | 20 +++ motulator/control/_common.py | 36 ++-- motulator/control/im/_obs_vhz.py | 10 +- motulator/control/im/_observers.py | 14 +- motulator/control/im/_vector.py | 8 +- motulator/control/sm/_flux_vector.py | 13 +- motulator/control/sm/_observers.py | 4 +- motulator/control/sm/_torque.py | 14 +- motulator/model/__init__.py | 5 +- motulator/model/_converter.py | 3 - motulator/model/_simulation.py | 131 +++++++++++++-- motulator/model/im/_drive.py | 83 +++------ motulator/model/sm/__init__.py | 2 + motulator/model/sm/_drive.py | 157 ++++++++++++------ 33 files changed, 449 insertions(+), 283 deletions(-) create mode 100644 examples/vector/plot_vector_ctrl_pmsm_2kw_diode.py diff --git a/.vscode/settings.json b/.vscode/settings.json index 2919a875c..90349ae1b 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,6 +1,6 @@ { "editor.rulers": [ - 80 + 79 ], "editor.formatOnType": false, "editor.formatOnSave": true, diff --git a/examples/flux_vector/plot_flux_vector_pmsm_2kw.py b/examples/flux_vector/plot_flux_vector_pmsm_2kw.py index 182b83d89..9ec4a6728 100644 --- a/examples/flux_vector/plot_flux_vector_pmsm_2kw.py +++ b/examples/flux_vector/plot_flux_vector_pmsm_2kw.py @@ -46,7 +46,7 @@ # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=1.6) # %% diff --git a/examples/flux_vector/plot_flux_vector_pmsyrm_5kw.py b/examples/flux_vector/plot_flux_vector_pmsyrm_5kw.py index 660afa457..9a1bf735d 100644 --- a/examples/flux_vector/plot_flux_vector_pmsyrm_5kw.py +++ b/examples/flux_vector/plot_flux_vector_pmsyrm_5kw.py @@ -3,11 +3,11 @@ ========================= This example simulates sensorless stator-flux-vector control of a 5.5-kW -PM-SyRM (Baldor ECS101M0H7EF4) drive. The machine model is parametrized using an -algebraic saturation model, fitted to the flux linkage maps measured using the -constant-speed test. For comparison, the measured data is plotted together with -the model predictions. Notice that the control system used in this example does -not consider the saturation, only the system model does. +PM-SyRM (Baldor ECS101M0H7EF4) drive. The machine model is parametrized using +an algebraic saturation model, fitted to the flux linkage maps measured using +the constant-speed test. For comparison, the measured data is plotted together +with the model predictions. Notice that the control system used in this example +does not consider the saturation, only the system model does. """ @@ -35,14 +35,15 @@ # available later. +# pylint: disable=too-many-locals def i_s(psi_s): """ - Saturation model for a 5.5-kW PM synchronous reluctance machine. + Saturation model for a 5.5-kW PM-SyRM. This model takes into account the bridge saturation in addition to the regular self- and cross-saturation effects of the d- and q-axis. The bridge saturation model is based on a nonlinear reluctance element in parallel - with the Norton-equivalent PM model. + with the Norton-equivalent PM model. Parameters ---------- @@ -56,8 +57,9 @@ def i_s(psi_s): Notes ----- - This model can also be used for other PM synchronous reluctance machines by - changing the model parameters. + For simplicity, the saturation model parameters are hard-coded in the + function below. This model can also be used for other PM-SyRMs by changing + the model parameters. """ # d-axis self-saturation @@ -131,8 +133,8 @@ def i_s(psi_s): plt.show() # %% -# Solve the PM flux linkage for the initial value of the stator flux, which is -# needed in the machine model below. +# Solve the PM flux linkage for the initial value of the stator flux linkage, +# which is needed in the machine model below. res = minimize_scalar( lambda psi_d: np.abs(i_s(psi_d)), bounds=(0, base.psi), method="bounded") @@ -143,12 +145,13 @@ def i_s(psi_s): machine = model.sm.SynchronousMachineSaturated( n_p=2, R_s=.63, current=i_s, psi_s0=psi_s0) -# Magnetically linear PMSyRM model for comparison +# Magnetically linear PM-SyRM model for comparison # machine = model.sm.SynchronousMachine( -# n_p=2, R_s=.63, L_d=18e-3, L_q=110e-3, psi_f=.47) +# n_p=2, R_s=.63, L_d=18e-3, L_q=110e-3, psi_f=.47) mechanics = model.Mechanics(J=.015) converter = model.Inverter(u_dc=540) mdl = model.sm.Drive(machine, mechanics, converter) +# mdl.pwm = model.CarrierComparison() # Enable the PWM model # %% # Configure the control system. @@ -178,7 +181,7 @@ def i_s(psi_s): # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) # %% diff --git a/examples/flux_vector/plot_flux_vector_syrm_7kw.py b/examples/flux_vector/plot_flux_vector_syrm_7kw.py index 46a83fbfd..c53ad01d9 100644 --- a/examples/flux_vector/plot_flux_vector_syrm_7kw.py +++ b/examples/flux_vector/plot_flux_vector_syrm_7kw.py @@ -80,7 +80,7 @@ def i_s(psi_s): # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) # %% 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 44f4ab97e..0cb0502ec 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py +++ b/examples/obs_vhz/plot_obs_vhz_ctrl_im_2kw.py @@ -60,11 +60,9 @@ # mdl.mech.tau_L_w = lambda w_M: np.sign(w_M)*k*w_M**2 # %% -# Create the simulation object and simulate it. You can also enable the PWM -# model (which makes simulation slower). One-sampling-period computational -# delay is modeled. +# Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False, delay=1) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) # %% diff --git a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw.py b/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw.py index 43055efb5..b77b141ef 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw.py +++ b/examples/obs_vhz/plot_obs_vhz_ctrl_pmsm_2kw.py @@ -49,11 +49,9 @@ mdl.mechanics.tau_L_t = Sequence(times, values) # %% -# Create the simulation object and simulate it. You can also enable the PWM -# model (which makes simulation slower). One-sampling-period computational -# delay is modeled. +# Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False, delay=1) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=8) # %% 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 94c9ba2e1..f2b2d9fc2 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 @@ -5,8 +5,8 @@ This example simulates observer-based V/Hz control of a 2.2-kW PMSM drive. The mechanical subsystem is modeled as a two-mass system. The resonance frequency of the mechanics is around 85 Hz. The mechanical parameters correspond to -[#Saa2015]_, except that the torsional damping is set to a smaller value in this -example. +[#Saa2015]_, except that the torsional damping is set to a smaller value in +this example. """ @@ -57,7 +57,7 @@ # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=1.2) # sphinx_gallery_thumbnail_number = 3 plot(sim, base) # Plot results in per-unit values diff --git a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py b/examples/obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py index d4da37ede..b292b3e73 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py +++ b/examples/obs_vhz/plot_obs_vhz_ctrl_pmsyrm_thor.py @@ -118,11 +118,9 @@ def i_s(psi_s): # mdl.mechanics.tau_L_t = Sequence(times, values) # %% -# Create the simulation object and simulate it. You can also enable the PWM -# model (which makes simulation slower). One-sampling-period computational -# delay is modeled. +# Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False, delay=1) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=8) # %% diff --git a/examples/obs_vhz/plot_obs_vhz_ctrl_syrm_7kw.py b/examples/obs_vhz/plot_obs_vhz_ctrl_syrm_7kw.py index c208efbde..5df7aaf28 100644 --- a/examples/obs_vhz/plot_obs_vhz_ctrl_syrm_7kw.py +++ b/examples/obs_vhz/plot_obs_vhz_ctrl_syrm_7kw.py @@ -27,8 +27,8 @@ # but the same model structure can also be used for other synchronous motors. # For PM motors, the magnetomotive force (MMF) of the PMs can be modeled using # constant current source `i_f` on the d-axis [#Awa2018]_, [#Jah1986]_. -# Correspondingly, this approach assumes that the MMFs of the d-axis current and -# of the PMs are in series. This model cannot capture the desaturation +# Correspondingly, this approach assumes that the MMFs of the d-axis current +# and of the PMs are in series. This model cannot capture the desaturation # phenomenon of thin iron ribs, see [#Arm2009]_ for details. For such motors, # look-up tables can be used. @@ -104,11 +104,9 @@ def i_s(psi_s): mdl.mechanics.tau_L_t = Sequence(times, values) # %% -# Create the simulation object and simulate it. You can also enable the PWM -# model (which makes simulation slower). One-sampling-period computational -# delay is modeled. +# Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False, delay=1) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=8) # %% diff --git a/examples/signal_inj/plot_signal_inj_pmsm_2kw.py b/examples/signal_inj/plot_signal_inj_pmsm_2kw.py index 94772d04b..ea8fc904d 100644 --- a/examples/signal_inj/plot_signal_inj_pmsm_2kw.py +++ b/examples/signal_inj/plot_signal_inj_pmsm_2kw.py @@ -54,7 +54,7 @@ # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) # %% diff --git a/examples/signal_inj/plot_signal_inj_syrm_7kw.py b/examples/signal_inj/plot_signal_inj_syrm_7kw.py index 56f46403a..0a4752000 100644 --- a/examples/signal_inj/plot_signal_inj_syrm_7kw.py +++ b/examples/signal_inj/plot_signal_inj_syrm_7kw.py @@ -46,7 +46,7 @@ # Speed reference times = np.array([0, .25, .25, .375, .5, .625, .75, .75, 1])*4 -values = np.array([0, 0, 1, 1, 0, -1, -1, 0, 0])*base.w*.1 +values = np.array([0, 0, 1, 1, 0, -1, -1, 0, 0])*.1*base.w ctrl.w_m_ref = Sequence(times, values) # External load torque times = np.array([0, .125, .125, .875, .875, 1])*4 @@ -56,7 +56,7 @@ # %% # Create the simulation object and simulate it. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) # %% diff --git a/examples/vector/plot_vector_ctrl_im_2kw.py b/examples/vector/plot_vector_ctrl_im_2kw.py index a140da1cb..75d94e3a0 100644 --- a/examples/vector/plot_vector_ctrl_im_2kw.py +++ b/examples/vector/plot_vector_ctrl_im_2kw.py @@ -31,7 +31,7 @@ def L_s(psi): """ - Stator inductance saturation model for a 2.2-kW machine. + Stator inductance saturation model for a 2.2-kW induction machine. Parameters ---------- @@ -59,13 +59,15 @@ def L_s(psi): n_p=2, R_s=3.7, R_r=2.5, L_ell=.023, L_s=L_s) # Unsaturated machine model, using its inverse-Γ parameters (uncomment to try) # machine = model.im.InductionMachineInvGamma( -# R_s=3.7, R_R=2.1, L_sgm=.021, L_M=.224, n_p=2) +# n_p=2, R_s=3.7, R_R=2.1, L_sgm=.021, L_M=.224) # Alternatively, configure the machine model using its Γ parameters # machine = model.im.InductionMachine( -# R_s=3.7, R_r=2.5, L_ell=.023, L_s=.245, n_p=2) +# n_p=2, R_s=3.7, R_r=2.5, L_ell=.023, L_s=.245) mechanics = model.Mechanics(J=.015) converter = model.Inverter(u_dc=540) mdl = model.im.Drive(machine, mechanics, converter) +# mdl.pwm = model.CarrierComparison() # Try to enable the PWM model +# mdl.delay = model.Delay(2) # Try longer computational delay # %% # Configure the control system. @@ -94,22 +96,11 @@ def L_s(psi): # ctrl.w_m_ref = lambda t: (t > .2)*(2*base.w) # mdl.mechanics.tau_L_t = lambda t: 0 -# Speed reversals under the rated load (uncomment to try, change t_stop=8 below) -# import numpy as np -# from motulator.helpers import Sequence -# times = np.array([0, .125, .25, .375, .5, .625, .75, .875, 1])*8 -# values = np.array([0, 0, 1, 1, 0, -1, -1, 0, 0])*base.w -# ctrl.w_m_ref = Sequence(times, values) -# # External load torque -# times = np.array([0, .125, .125, .875, .875, 1])*8 -# values = np.array([0, 0, 1, 1, 0, 0])*base.tau_nom -# mdl.mechanics.tau_L_t = Sequence(times, values) - # %% # Create the simulation object and simulate it. You can also enable the PWM # model (which makes simulation slower). -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=1.5) # %% diff --git a/examples/vector/plot_vector_ctrl_pmsm_2kw.py b/examples/vector/plot_vector_ctrl_pmsm_2kw.py index 6aa5fa923..beefa2313 100644 --- a/examples/vector/plot_vector_ctrl_pmsm_2kw.py +++ b/examples/vector/plot_vector_ctrl_pmsm_2kw.py @@ -1,17 +1,16 @@ """ -2.2-kW PMSM -=========== +2.2-kW PMSM, diode bridge +========================= -This example simulates sensorless vector control of a 2.2-kW PMSM drive. +This example simulates sensorless vector control of a 2.2-kW PMSM drive. """ # %% # Imports. -import numpy as np from motulator import model, control -from motulator import BaseValues, Sequence, plot, plot_extra +from motulator import BaseValues, plot # %% # Compute base values based on the nominal values (just for figures). @@ -27,6 +26,7 @@ mechanics = model.Mechanics(J=.015) converter = model.Inverter(u_dc=540) mdl = model.sm.Drive(machine, mechanics, converter) +# mdl.pwm = model.CarrierComparison() # Enable the PWM model # %% # Configure the control system. @@ -40,30 +40,15 @@ # Set the speed reference and the external load torque. # Speed reference -times = np.array([0, .125, .25, .375, .5, .625, .75, .875, 1])*4 -values = np.array([0, 0, 1, 1, 0, -1, -1, 0, 0])*base.w -ctrl.w_m_ref = Sequence(times, values) -# External load torque -times = np.array([0, .125, .125, .875, .875, 1])*4 -values = np.array([0, 0, 1, 1, 0, 0])*base.tau_nom -mdl.mechanics.tau_L_t = Sequence(times, values) +ctrl.w_m_ref = lambda t: (t > .2)*2*base.w -# mdl.mechanics.tau_L_t = lambda t: (t > .8)*base.tau_nom*.7 -# ctrl.w_m_ref = lambda t: (t > .2)*(2*base.w) +# External load torque +mdl.mechanics.tau_L_t = lambda t: (t > .8)*.7*base.tau_nom # %% # 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) +# Simulate the system and plot results in per-unit values +sim = model.Simulation(mdl, ctrl) +sim.simulate(t_stop=1.4) plot(sim, base) -# Plot a zoomed view -plot_extra(sim, t_span=(1.1, 1.125), base=base) diff --git a/examples/vector/plot_vector_ctrl_pmsm_2kw_diode.py b/examples/vector/plot_vector_ctrl_pmsm_2kw_diode.py new file mode 100644 index 000000000..cf22e3794 --- /dev/null +++ b/examples/vector/plot_vector_ctrl_pmsm_2kw_diode.py @@ -0,0 +1,58 @@ +""" +2.2-kW PMSM, diode bridge +========================= + +This example simulates sensorless vector control of a 2.2-kW PMSM drive, +equipped with a diode bridge rectifier. + +""" + +# %% +# Imports. + +from motulator import model, control +from motulator import BaseValues, plot, plot_extra + +# %% +# Compute base values based on the nominal values (just for figures). + +base = BaseValues( + U_nom=370, I_nom=4.3, f_nom=75, tau_nom=14, P_nom=2.2e3, n_p=3) + +# %% +# Configure the system model. + +machine = model.sm.SynchronousMachine( + n_p=3, R_s=3.6, L_d=.036, L_q=.051, psi_f=.545) +mechanics = model.Mechanics(J=.015) +converter = model.FrequencyConverter(L=2e-3, C=235e-6, U_g=400, f_g=50) +mdl = model.sm.DriveWithDiodeBridge(machine, mechanics, converter) +mdl.pwm = model.CarrierComparison() # Enable the PWM model + +# %% +# Configure the control system. + +par = control.sm.ModelPars( + n_p=3, R_s=3.6, L_d=.036, L_q=.051, psi_f=.545, J=.015) +ref = control.sm.CurrentReferencePars(par, w_m_nom=base.w, i_s_max=1.5*base.i) +ctrl = control.sm.VectorCtrl(par, ref, T_s=250e-6, sensorless=True) + +# %% +# Set the speed reference and the external load torque. + +# Speed reference +ctrl.w_m_ref = lambda t: (t > .2)*base.w + +# External load torque +mdl.mechanics.tau_L_t = lambda t: (t > .6)*base.tau_nom + +# %% +# Create the simulation object and simulate it. + +# Simulate the system +sim = model.Simulation(mdl, ctrl) +sim.simulate(t_stop=1) + +# Plot results in per-unit values +plot(sim, base) +plot_extra(sim, base, t_span=(.8, .825)) diff --git a/examples/vector/plot_vector_ctrl_pmsyrm_thor.py b/examples/vector/plot_vector_ctrl_pmsyrm_thor.py index 3ad2f40f0..816d00a31 100644 --- a/examples/vector/plot_vector_ctrl_pmsyrm_thor.py +++ b/examples/vector/plot_vector_ctrl_pmsyrm_thor.py @@ -62,12 +62,8 @@ mdl.mechanics.tau_L_w = lambda w_M: k*w_M**2*np.sign(w_M) # %% -# Create the simulation object and simulate it. +# Create the simulation object, simulate, and plot results in per-unit values. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=.6) - -# %% -# Plot results in per-unit values. - plot(sim, base) diff --git a/examples/vector/plot_vector_ctrl_syrm_7kw.py b/examples/vector/plot_vector_ctrl_syrm_7kw.py index 66cf41f73..021ac5127 100644 --- a/examples/vector/plot_vector_ctrl_syrm_7kw.py +++ b/examples/vector/plot_vector_ctrl_syrm_7kw.py @@ -50,12 +50,8 @@ mdl.mechanics.tau_L_t = Sequence(times, values) # %% -# Create the simulation object and simulate it. +# Create the simulation object, simulate, and plot results in per-unit values. -sim = model.Simulation(mdl, ctrl, pwm=False) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=4) - -# %% -# Plot results in per-unit values. - plot(sim, base) diff --git a/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py b/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py index 9c8c26869..2ecfc4be5 100644 --- a/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py +++ b/examples/vhz/plot_vhz_ctrl_6step_im_2kw.py @@ -31,6 +31,7 @@ mechanics = model.Mechanics(J=.015) converter = model.Inverter(u_dc=540) mdl = model.im.Drive(machine, mechanics, converter) +mdl.pwm = model.CarrierComparison() # Enable the PWM model # %% # Control system (parametrized as open-loop V/Hz control). @@ -46,7 +47,7 @@ # Speed reference times = np.array([0, .1, .3, 1])*2 -values = np.array([0, 0, 1, 1])*base.w*2 +values = np.array([0, 0, 1, 1])*2*base.w ctrl.w_m_ref = Sequence(times, values) # Quadratic load torque profile (corresponding to pumps and fans) @@ -56,10 +57,9 @@ mdl.mechanics.tau_L_t = lambda t: (t > 1.)*base.tau_nom*0 # %% -# 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, pwm=True) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=2) # %% diff --git a/examples/vhz/plot_vhz_ctrl_im_2kw.py b/examples/vhz/plot_vhz_ctrl_im_2kw.py index f2b01b526..b70afe159 100644 --- a/examples/vhz/plot_vhz_ctrl_im_2kw.py +++ b/examples/vhz/plot_vhz_ctrl_im_2kw.py @@ -30,6 +30,7 @@ # Frequency converter with a diode bridge converter = model.FrequencyConverter(L=2e-3, C=235e-6, U_g=400, f_g=50) mdl = model.im.DriveWithDiodeBridge(machine, mechanics, converter) +mdl.pwm = model.CarrierComparison() # Enable the PWM model # %% # Control system (parametrized as open-loop V/Hz control). @@ -43,20 +44,20 @@ # Set the speed reference and the external load torque. More complicated # signals could be defined as functions. -ctrl.w_m_ref = lambda t: (t > .2)*(1.*base.w) +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) # Stepwise load torque at t = 1 s, 20% of the rated torque -mdl.mechanics.tau_L_t = lambda t: (t > 1.)*base.tau_nom*.2 +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. -sim = model.Simulation(mdl, ctrl, pwm=True) +sim = model.Simulation(mdl, ctrl) sim.simulate(t_stop=1.5) # %% @@ -65,12 +66,12 @@ # .. note:: # The DC link of this particular example is actually unstable at 1-p.u. # speed at the rated load torque, since the inverter looks like a negative -# resistance to the DC link. You could notice this instability if simulating -# a longer period (e.g. set `t_stop=2`). For analysis, see e.g., [#Hin2007]_. +# resistance to the DC link. You can notice this instability if simulating a +# longer period (e.g. set `t_stop=2`). For analysis, see e.g., [#Hin2007]_. # sphinx_gallery_thumbnail_number = 2 plot(sim, base) -plot_extra(sim, t_span=(1.1, 1.125), base=base) +plot_extra(sim, base, t_span=(1.1, 1.125)) # %% # .. rubric:: References diff --git a/motulator/__init__.py b/motulator/__init__.py index a7a5955c1..a150583a7 100644 --- a/motulator/__init__.py +++ b/motulator/__init__.py @@ -3,7 +3,7 @@ This software includes continuous-time simulation models for induction machines and synchronous machines. Furthermore, selected examples of discrete-time -control algortihms are also included as well as various utilities. +control algorithms are also included as well as various utilities. """ diff --git a/motulator/_helpers.py b/motulator/_helpers.py index 088ab1bd8..7acbf3ce3 100644 --- a/motulator/_helpers.py +++ b/motulator/_helpers.py @@ -194,3 +194,23 @@ def __call__(self, t): """ return self.initial_value + (t >= self.step_time)*self.step_value + + +# %% +def wrap(theta): + """ + Limit the angle into the range [-pi, pi). + + Parameters + ---------- + theta : float + Angle (rad). + + Returns + ------- + float + Limited angle. + + """ + + return np.mod(theta + np.pi, 2*np.pi) - np.pi diff --git a/motulator/control/_common.py b/motulator/control/_common.py index ac09d96d7..097034c04 100644 --- a/motulator/control/_common.py +++ b/motulator/control/_common.py @@ -67,9 +67,9 @@ def six_step_overmodulation(u_s_ref, u_dc): References ---------- - .. [#Bol1997] Bolognani, Zigliotto, "Novel digital continuous control of - SVM inverters in the overmodulation range," IEEE Trans. Ind. Appl., - 1997, https://doi.org/10.1109/28.568019 + .. [#Bol1997] Bolognani, Zigliotto, "Novel digital continuous control + of SVM inverters in the overmodulation range," IEEE Trans. Ind. + Appl., 1997, https://doi.org/10.1109/28.568019 """ # Limited magnitude @@ -104,8 +104,8 @@ def duty_ratios(u_s_ref, u_dc): """ Compute the duty ratios for three-phase space-vector PWM. - This computes the duty ratios corresponding to standard space-vector PWM - and minimum-amplitude-error overmodulation [#Hav1999]_. + This computes the duty ratios corresponding to standard space-vector + PWM and minimum-amplitude-error overmodulation [#Hav1999]_. Parameters ---------- @@ -202,14 +202,14 @@ class PICtrl: where `u` is the controller output, `y_ref` is the reference signal, `y` is the feedback signal, and `1/s` refers to integration. The standard PI - controller is obtained by choosing ``k_t = k_p``. The integrator anti-windup - is implemented based on the realized controller output. + controller is obtained by choosing ``k_t = k_p``. The integrator + anti-windup is implemented based on the realized controller output. Notes ----- 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`. + corresponds to the rotor angular speed `w_M` and `u` to the torque + reference `tau_M_ref`. Parameters ---------- @@ -289,7 +289,7 @@ class SpeedCtrl(PICtrl): """ 2DOF PI speed controller. - This provides an interface for a speed controller. The gains are initialized + This is an interface for a speed controller. The gains are initialized based on the desired closed-loop bandwidth and the rotor inertia estimate. Parameters @@ -322,10 +322,10 @@ class ComplexPICtrl: where `u` is the controller output, `i_ref` is the reference signal, `i` is the feedback signal, `w` is the angular speed of synchronous coordinates, - and `1/s` refers to integration. This algorithm is compatible with both real - and complex signals. The 1DOF version is obtained by setting ``k_t = k_p``. - The integrator anti-windup is implemented based on the realized controller - output. + and `1/s` refers to integration. This algorithm is compatible with both + real and complex signals. The 1DOF version is obtained by setting + ``k_t = k_p``. The integrator anti-windup is implemented based on the + realized controller output. Parameters ---------- @@ -486,7 +486,7 @@ def update(self, T_s): # %% class Ctrl: - """Base class for the control system.""" + """Base class for control systems.""" def __init__(self): self.clear() @@ -507,9 +507,9 @@ 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 ---------- diff --git a/motulator/control/im/_obs_vhz.py b/motulator/control/im/_obs_vhz.py index b58f65a01..d943ff7f0 100644 --- a/motulator/control/im/_obs_vhz.py +++ b/motulator/control/im/_obs_vhz.py @@ -1,15 +1,15 @@ """ Observer-based V/Hz control for induction machine drives. -This implements the observer-based V/Hz control method described in [#Tii2022]_. -The state-feedback control law is in the alternative form which uses an +This implements the observer-based V/Hz control method [#Tii2022]_. The +state-feedback control law is in the alternative form which uses an intermediate stator current reference. References ---------- -.. [#Tii2022] Tiitinen, Hinkkanen, Harnefors, "Stable and passive observer-based - V/Hz control for induction motors," Proc. IEEE ECCE, Detroit, MI, Oct. 2022, - https://doi.org/10.1109/ECCE50734.2022.9948057 +.. [#Tii2022] Tiitinen, Hinkkanen, Harnefors, "Stable and passive + observer-based V/Hz control for induction motors," Proc. IEEE ECCE, Detroit, + MI, Oct. 2022, https://doi.org/10.1109/ECCE50734.2022.9948057 """ diff --git a/motulator/control/im/_observers.py b/motulator/control/im/_observers.py index 4012689e4..15bae496f 100644 --- a/motulator/control/im/_observers.py +++ b/motulator/control/im/_observers.py @@ -10,9 +10,9 @@ class Observer: rotor flux coordinates. 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. + Both sensored and sensorless operation are supported. The observer + structure is similar to [#Hin2010]_. The observer operates in estimated + rotor flux coordinates. Parameters ---------- @@ -151,8 +151,8 @@ class FullOrderObserver: Full-order flux observer for induction machines in estimated rotor flux coordinates. - This class implements a full-order flux observer for induction machines. The - observer structure is similar to [#Tii2023]_. The observer operates in + This class implements a full-order flux observer for induction machines. + The observer structure is similar to [#Tii2023]_. The observer operates in estimated rotor flux coordinates. Parameters @@ -184,8 +184,8 @@ class FullOrderObserver: References ---------- .. [#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 + observer revisited: Closed-form design for induction motor drives," + Proc. IEEE SLED, 2023, https://doi.org/10.1109/SLED57582.2023.10261359 """ diff --git a/motulator/control/im/_vector.py b/motulator/control/im/_vector.py index 3490f7878..89c3f7d76 100644 --- a/motulator/control/im/_vector.py +++ b/motulator/control/im/_vector.py @@ -196,8 +196,8 @@ class CurrentReferencePars: Parameters for reference generation. This dataclass stores the nominal and limit values needed for reference - generation. For calculating the rotor flux reference, the machine parameters - are also required. + generation. For calculating the rotor flux reference, the machine + parameters are also required. Parameters ---------- @@ -240,8 +240,8 @@ class CurrentReference: """ Current reference generation. - In the base-speed region, the current reference in rotor-flux coordinates is - given by:: + In the base-speed region, the current reference in rotor-flux coordinates + is given by:: i_s_ref = psi_R_nom/L_M + 1j*tau_M_ref/(1.5*n_p*abs(psi_R)) diff --git a/motulator/control/sm/_flux_vector.py b/motulator/control/sm/_flux_vector.py index 9dafceb11..1c10cc0af 100644 --- a/motulator/control/sm/_flux_vector.py +++ b/motulator/control/sm/_flux_vector.py @@ -17,13 +17,14 @@ class FluxVectorCtrl(Ctrl): This class implements a variant of stator-flux-vector control [#Pel2009]_. Rotor coordinates as well as decoupling between the stator flux and torque - channels are used according to [#Awa2019b]_. Here, the stator flux magnitude - and the electromagnetic torque are selected as controllable variables. + channels are used according to [#Awa2019b]_. Here, the stator flux + magnitude and the electromagnetic torque are selected as controllable + variables. Notes ----- - Proportional controllers are used for simplicity. The magnetic saturation is - not considered in this implementation. + Proportional controllers are used for simplicity. The magnetic saturation + is not considered in this implementation. Parameters ---------- @@ -48,7 +49,7 @@ class FluxVectorCtrl(Ctrl): Flux and torque reference generator. speed_ctrl : SpeedCtrl Speed controller. - w_m_ref : float + w_m_ref : callable Speed reference (electrical rad/s) as a function of time (s). pwm : PWM Pulse-width modulation. @@ -144,7 +145,7 @@ def __call__(self, mdl): # Auxiliary current i_a = psi_s.real/self.L_q + 1j*psi_s.imag/self.L_d - i_s - # Torque-production factor (c_tau = 0 corresponds to the MTPV condition) + # Torque-production factor, c_tau = 0 corresponds to the MTPV condition c_tau = 1.5*self.n_p*np.real(i_a*np.conj(psi_s)) # References for the flux and torque controllers diff --git a/motulator/control/sm/_observers.py b/motulator/control/sm/_observers.py index 55154bfde..a36d2ec34 100644 --- a/motulator/control/sm/_observers.py +++ b/motulator/control/sm/_observers.py @@ -139,8 +139,8 @@ class FluxObserver: Attributes ---------- delta : float - Angle estimate of the coordinate system with respect - to the d-axis of the rotor (electrical rad). + Angle estimate of the coordinate system with respect to the d-axis of + the rotor (electrical rad). psi_s : complex Stator flux estimate (Vs). diff --git a/motulator/control/sm/_torque.py b/motulator/control/sm/_torque.py index 15fde2d32..08066693e 100644 --- a/motulator/control/sm/_torque.py +++ b/motulator/control/sm/_torque.py @@ -1,10 +1,10 @@ """ Torque characteristics for synchronous machines. -This contains computation and plotting of torque characteristics for synchronous -machines, including the MTPA and MTPV loci [#Mor1990]_. The methods can be used -to define look-up tables for control and to analyze the characteristics. This -version omits the magnetic saturation. +This contains computation and plotting of torque characteristics for +synchronous machines, including the MTPA and MTPV loci [#Mor1990]_. The methods +can be used to define look-up tables for control and to analyze the +characteristics. This version omits the magnetic saturation. References ---------- @@ -232,7 +232,8 @@ def mtpa_locus(self, i_s_max, psi_s_min=None, N=20): Parameters ---------- i_s_max : float - Maximum stator current magnitude (A) at which the locus is computed. + Maximum stator current magnitude (A) at which the locus is + computed. psi_s_min : float, optional Minimum stator flux magnitude (Vs) at which the locus is computed. N : int, optional @@ -295,7 +296,8 @@ def mtpv_locus(self, psi_s_max=None, i_s_max=None, N=20): Maximum stator flux magnitude (Vs) at which the locus is computed. Either `psi_s_max` or `i_s_max` must be given. i_s_max : float, optional - Maximum stator current magnitude (A) at which the locus is computed. + Maximum stator current magnitude (A) at which the locus is + computed. N : int, optional Amount of points. The default is 20. diff --git a/motulator/model/__init__.py b/motulator/model/__init__.py index 7a4e247e9..14f0a9d5a 100644 --- a/motulator/model/__init__.py +++ b/motulator/model/__init__.py @@ -2,7 +2,8 @@ from motulator.model._mechanics import Mechanics, MechanicsTwoMass from motulator.model._converter import FrequencyConverter, Inverter -from motulator.model._simulation import CarrierComparison, Simulation +from motulator.model._simulation import ( + CarrierComparison, Simulation, Delay, zoh) import motulator.model.im as im import motulator.model.sm as sm @@ -14,6 +15,8 @@ "Inverter", "CarrierComparison", "Simulation", + "Delay", + "zoh", "im", "sm", ] diff --git a/motulator/model/_converter.py b/motulator/model/_converter.py index ab790fec4..61889743d 100644 --- a/motulator/model/_converter.py +++ b/motulator/model/_converter.py @@ -24,7 +24,6 @@ class Inverter: def __init__(self, u_dc): self.u_dc0 = u_dc - self.q = 0j # Switching state vector @staticmethod def ac_voltage(q, u_dc): @@ -114,8 +113,6 @@ def __init__(self, L, C, U_g, f_g): self.u_g = np.sqrt(2/3)*U_g # Grid angular frequeyncy self.w_g = 2*np.pi*f_g - # Switching state vector - self.q = 0j def grid_voltages(self, t): """ diff --git a/motulator/model/_simulation.py b/motulator/model/_simulation.py index 2b9f3139c..50b3debf6 100644 --- a/motulator/model/_simulation.py +++ b/motulator/model/_simulation.py @@ -4,6 +4,7 @@ from scipy.integrate import solve_ivp from scipy.io import savemat from motulator._helpers import abc2complex +from motulator._utils import Bunch # %% @@ -155,7 +156,7 @@ def __call__(self, T_s, d_abc): # %% -def _zoh(T_s, d_abc): +def zoh(T_s, d_abc): """ Zero-order hold of the duty ratios over the sampling period. @@ -189,25 +190,16 @@ class Simulation: Parameters ---------- - mdl : Drive + mdl : Model Continuous-time system model. ctrl : Ctrl Discrete-time controller. - delay : int, optional - Amount of computational delays. The default is 1. - pwm : bool, optional - Enable carrier comparison. The default is False. """ - def __init__(self, mdl=None, ctrl=None, delay=1, pwm=False): + def __init__(self, mdl=None, ctrl=None): self.mdl = mdl self.ctrl = ctrl - self._delay = Delay(delay) - if pwm: - self._pwm = CarrierComparison() - else: - self._pwm = _zoh def simulate(self, t_stop=1, max_step=np.inf): """ @@ -222,7 +214,7 @@ def simulate(self, t_stop=1, max_step=np.inf): Notes ----- - Other options of `solve_ivp` could be easily changed if needed, but, for + Other options of `solve_ivp` could be easily used if needed, but, for simplicity, only `max_step` is included as an option of this method. """ @@ -243,17 +235,17 @@ def _simulation_loop(self, t_stop, max_step): T_s, d_abc_ref = self.ctrl(self.mdl) # Computational delay model - d_abc = self._delay(d_abc_ref) + d_abc = self.mdl.delay(d_abc_ref) # Carrier comparison - t_steps, q = self._pwm(T_s, d_abc) + t_steps, q = self.mdl.pwm(T_s, d_abc) # Loop over the sampling period T_s for i, t_step in enumerate(t_steps): if t_step > 0: # Update the switching state - self.mdl.converter.q = q[i] + self.mdl.q = q[i] # Get initial values x0 = self.mdl.get_initial_values() @@ -267,7 +259,7 @@ def _simulation_loop(self, t_stop, max_step): self.mdl.set_initial_values(t0_new, x0_new) # Save the solution - sol.q = len(sol.t)*[self.mdl.converter.q] + sol.q = len(sol.t)*[self.mdl.q] self.mdl.save(sol) def save_mat(self, name="sim"): @@ -282,3 +274,108 @@ def save_mat(self, name="sim"): """ savemat(name + "_mdl_data" + ".mat", self.mdl.data) savemat(name + "_ctrl_data" + ".mat", self.ctrl.data) + + +# %% +class Model: + """ + Base class for continuous-time system models. + + This base class is a template for a system model that interconnects the + subsystems and provides an interface to the solver. + + Parameters + ---------- + pwm : zoh | CarrierComparison, optional + Zero-order hold of duty ratios or carrier comparison. If None, the + default is `zoh`. + delay : int, optional + Amount of computational delays. The default is 1. + + """ + + def __init__(self, pwm=None, delay=1): + self.delay = Delay(delay) + self.pwm = zoh if pwm is None else pwm + self.t0 = 0 # Initial time + self.q = 0j # Switching state vector + 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. + + """ + # Initial time + self.t0 = 0 + # Solution will be stored in the following lists + self.data = Bunch() + self.data.t, self.data.q = [], [] + + def get_initial_values(self): + """ + Get the initial values. + + Returns + ------- + x0 : complex list + Initial values of the state variables. + + """ + raise NotImplementedError + + def set_initial_values(self, t0, x0): + """ + Set the initial values. + + Parameters + ---------- + t0 : float + Initial time (s). + x0 : complex ndarray + Initial values of the state variables. + + """ + raise NotImplementedError + + def f(self, t, x): + """ + Compute the complete state derivative list for the solver. + + Parameters + ---------- + t : float + Time (s). + x : complex ndarray + State vector. + + Returns + ------- + complex list + State derivatives. + + """ + raise NotImplementedError + + def save(self, sol): + """ + Save the solution. + + Parameters + ---------- + sol : Bunch + Solution from the solver. + + """ + self.data.t.extend(sol.t) + self.data.q.extend(sol.q) + + def post_process(self): + """Transform the lists to the ndarray format and post-process them.""" + # From lists to the ndarray + for key in self.data: + self.data[key] = np.asarray(self.data[key]) diff --git a/motulator/model/im/_drive.py b/motulator/model/im/_drive.py index 0615b6c6f..27d2c7f1e 100644 --- a/motulator/model/im/_drive.py +++ b/motulator/model/im/_drive.py @@ -6,8 +6,8 @@ """ import numpy as np -from motulator._helpers import abc2complex, complex2abc -from motulator._utils import Bunch +from motulator._helpers import abc2complex, complex2abc, wrap +from motulator.model._simulation import Model # %% @@ -41,8 +41,8 @@ class InductionMachine: References ---------- - .. [#Sle1989] Slemon, "Modelling of induction machines for electric drives," - IEEE Trans. Ind. Appl., 1989, https://doi.org/10.1109/28.44251. + .. [#Sle1989] Slemon, "Modelling of induction machines for electric + drives," IEEE Trans. Ind. Appl., 1989, https://doi.org/10.1109/28.44251. """ @@ -164,8 +164,8 @@ class InductionMachineSaturated(InductionMachine): """ Γ-equivalent model of an induction machine model with main-flux saturation. - This extends the InductionMachine class with a main-flux magnetic saturation - model:: + This extends the InductionMachine class with a main-flux magnetic + saturation model:: L_s = L_s(abs(psi_ss)) @@ -229,13 +229,12 @@ def __init__(self, n_p, R_s, R_R, L_sgm, L_M): # %% -class Drive: +class Drive(Model): """ Continuous-time model for an induction machine drive. - This interconnects the subsystems of an induction machine drive and provides - an interface to the solver. More complicated systems could be modeled using - a similar template. + This interconnects the subsystems of an induction machine drive and + provides an interface to the solver. Parameters ---------- @@ -249,25 +248,14 @@ class Drive: """ def __init__(self, machine=None, mechanics=None, converter=None): + super().__init__() self.machine = machine self.mechanics = mechanics self.converter = converter - self.t0 = 0 # Initial time - 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 = [], [] + """Clear the simulation data of the system model.""" + super().clear() self.data.psi_ss, self.data.psi_rs = [], [] self.data.theta_M, self.data.w_M = [], [] @@ -306,8 +294,7 @@ def set_initial_values(self, t0, x0): self.machine.psi_rs0 = x0[1] # x0[2].imag and x0[3].imag are always zero self.mechanics.w_M0 = x0[2].real - # Limit theta_M0 = x0[3].real into [-pi, pi) - self.mechanics.theta_M0 = np.mod(x0[3].real + np.pi, 2*np.pi) - np.pi + self.mechanics.theta_M0 = wrap(x0[3].real) def f(self, t, x): """ @@ -329,8 +316,7 @@ 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.converter.q, self.converter.u_dc0) + u_ss = 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) mechanics_f = self.mechanics.f(t, w_M, tau_M) @@ -338,17 +324,8 @@ def f(self, t, x): return machine_f + mechanics_f def save(self, sol): - """ - Save the solution. - - Parameters - ---------- - sol : Bunch object - Solution from the solver. - - """ - self.data.t.extend(sol.t) - self.data.q.extend(sol.q) + """Save the solution.""" + super().save(sol) self.data.psi_ss.extend(sol.y[0]) self.data.psi_rs.extend(sol.y[1]) self.data.w_M.extend(sol.y[2].real) @@ -356,17 +333,12 @@ def save(self, sol): def post_process(self): """Transform the lists to the ndarray format and post-process them.""" - # From lists to the ndarray - for key in self.data: - self.data[key] = np.asarray(self.data[key]) - + super().post_process() # Some useful variables self.data.i_ss, _, self.data.tau_M = self.machine.magnetic( self.data.psi_ss, self.data.psi_rs) - self.data.theta_M = np.mod( # Limit into [-pi, pi) - self.data.theta_M + np.pi, 2*np.pi) - np.pi - self.data.theta_m = np.mod( # Limit into [-pi, pi) - self.machine.n_p*self.data.theta_M + np.pi, 2*np.pi) - np.pi + self.data.theta_M = wrap(self.data.theta_M) + self.data.theta_m = wrap(self.machine.n_p*self.data.theta_M) self.data.w_m = self.machine.n_p*self.data.w_M self.data.tau_L = ( self.mechanics.tau_L_t(self.data.t) + @@ -407,7 +379,6 @@ 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): @@ -434,13 +405,13 @@ 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.converter.q, u_dc) - mach_f, i_ss, tau_M = self.machine.f(psi_ss, psi_rs, u_ss, w_M) - i_dc = self.converter.dc_current(self.converter.q, i_ss) + 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) + i_dc = self.converter.dc_current(self.q, i_ss) # Return the list of state derivatives return ( - mach_f + self.mechanics.f(t, w_M, tau_M) + + machine_f + self.mechanics.f(t, w_M, tau_M) + self.converter.f(t, u_dc, i_L, i_dc)) def save(self, sol): @@ -474,9 +445,6 @@ class DriveTwoMassMechanics(Drive): """ Induction machine drive with two-mass mechanics. - This interconnects the subsystems of an induction machine drive and provides - an interface to the solver. - Parameters ---------- machine : InductionMachine | InductionMachineSaturated @@ -508,15 +476,14 @@ def set_initial_values(self, t0, x0): """Extend the base class.""" super().set_initial_values(t0, x0[0:4]) self.mechanics.w_L0 = x0[4].real - self.mechanics.theta_ML0 = np.mod(x0[5].real + np.pi, 2*np.pi) - np.pi + self.mechanics.theta_ML0 = wrap(x0[5].real) def f(self, t, x): """Override the base class.""" # 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.converter.q, self.converter.u_dc0) + u_ss = 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) mechanics_f = self.mechanics.f(t, w_M, w_L, theta_ML, tau_M) diff --git a/motulator/model/sm/__init__.py b/motulator/model/sm/__init__.py index 18cf84cd9..94abf2bad 100644 --- a/motulator/model/sm/__init__.py +++ b/motulator/model/sm/__init__.py @@ -2,6 +2,7 @@ from motulator.model.sm._drive import ( Drive, + DriveWithDiodeBridge, DriveTwoMassMechanics, SynchronousMachine, SynchronousMachineSaturated, @@ -15,6 +16,7 @@ __all__ = [ "Drive", + "DriveWithDiodeBridge", "DriveTwoMassMechanics", "SynchronousMachine", "SynchronousMachineSaturated", diff --git a/motulator/model/sm/_drive.py b/motulator/model/sm/_drive.py index 2fde6a201..8e2ebdcaf 100644 --- a/motulator/model/sm/_drive.py +++ b/motulator/model/sm/_drive.py @@ -5,8 +5,8 @@ """ import numpy as np -from motulator._helpers import complex2abc -from motulator._utils import Bunch +from motulator._helpers import abc2complex, complex2abc, wrap +from motulator.model._simulation import Model # %% @@ -95,7 +95,7 @@ def f(self, psi_s, u_s, w_M): Returns ------- complex list, length 2 - Time derivative of the state vector, [dpsi_s, dtheta_m0] + Time derivative of the state vector, [dpsi_s, dtheta_m] i_s : complex Stator current (A). tau_M : float @@ -104,7 +104,7 @@ def f(self, psi_s, u_s, w_M): Notes ----- In addition to the state derivative, this method also returns the - output signals (stator current `i_ss` and torque `tau_M`) needed for + output signals (stator current `i_s` and torque `tau_M`) needed for interconnection with other subsystems. This avoids overlapping computation in simulation. @@ -166,13 +166,12 @@ def __init__(self, n_p, R_s, current, psi_s0=0j): # %% -class Drive: +class Drive(Model): """ Continuous-time model for a synchronous machine drive. This interconnects the subsystems of a synchronous machine drive and - provides an interface to the solver. More complicated systems could be - modeled using a similar template. + provides an interface to the solver. Parameters ---------- @@ -186,25 +185,14 @@ class Drive: """ def __init__(self, machine=None, mechanics=None, converter=None): + super().__init__() self.machine = machine self.mechanics = mechanics self.converter = converter - self.t0 = 0 - 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 - # Solution will be stored in the following lists - self.data = Bunch() - self.data.t, self.data.q = [], [] + """Clear the simulation data of the system model.""" + super().clear() self.data.psi_s, self.data.theta_m = [], [] self.data.w_M, self.data.theta_M = [], [] @@ -243,9 +231,8 @@ def set_initial_values(self, t0, x0): # x0[1:3].imag are always zero self.machine.theta_m0 = x0[1].real self.mechanics.w_M0 = x0[2].real - # Limit the angles into [-pi, pi) - self.mechanics.theta_m0 = np.mod(x0[1].real + np.pi, 2*np.pi) - np.pi - self.mechanics.theta_M0 = np.mod(x0[3].real + np.pi, 2*np.pi) - np.pi + self.mechanics.theta_m0 = wrap(x0[1].real) + self.mechanics.theta_M0 = wrap(x0[3].real) def f(self, t, x): """ @@ -268,8 +255,7 @@ 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.converter.q, self.converter.u_dc0) + 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 # State derivatives @@ -280,17 +266,8 @@ def f(self, t, x): return machine_f + mechanics_f def save(self, sol): - """ - Save the solution. - - Parameters - ---------- - sol : Bunch - Solution from the solver. - - """ - self.data.t.extend(sol.t) - self.data.q.extend(sol.q) + """Save the solution.""" + super().save(sol) self.data.psi_s.extend(sol.y[0]) self.data.theta_m.extend(sol.y[1].real) self.data.w_M.extend(sol.y[2].real) @@ -298,10 +275,7 @@ def save(self, sol): def post_process(self): """Transform the lists to the ndarray format and post-process them.""" - # From lists to the ndarray - for key in self.data: - self.data[key] = np.asarray(self.data[key]) - + super().post_process() # Compute some useful quantities self.data.i_s, self.data.tau_M = self.machine.magnetic(self.data.psi_s) self.data.w_m = self.machine.n_p*self.data.w_M @@ -310,21 +284,103 @@ def post_process(self): self.mechanics.tau_L_w(self.data.w_M)) self.data.u_ss = self.converter.ac_voltage( self.data.q, self.converter.u_dc0) - self.data.theta_m = np.mod( # Limit into [-pi, pi) - self.data.theta_m + np.pi, 2*np.pi) - np.pi - self.data.theta_M = np.mod( # Limit into [-pi, pi) - self.data.theta_M + np.pi, 2*np.pi) - np.pi + self.data.theta_m = wrap(self.data.theta_m) + self.data.theta_M = wrap(self.data.theta_M) self.data.i_ss = self.data.i_s*np.exp(1j*self.data.theta_m) +# %% +class DriveWithDiodeBridge(Drive): + """ + Synchronous machine drive equipped with a diode bridge. + + This model extends the Drive class with a model for a three-phase diode + bridge fed from stiff supply voltages. The DC bus is modeled as an inductor + and a capacitor. + + Parameters + ---------- + machine : SynchronousMachine | SynchronousMachineSaturated + Induction machine model. + mechanics : Mechanics + Mechanics model. + converter : FrequencyConverter + Frequency converter model. + + """ + + 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): + """Extend the base class.""" + x0 = super().get_initial_values() + [ + self.converter.u_dc0, self.converter.i_L0 + ] + return x0 + + def set_initial_values(self, t0, x0): + """Extend the base class.""" + super().set_initial_values(t0, x0[0:4]) + self.converter.u_dc0 = x0[4].real + self.converter.i_L0 = x0[5].real + + def f(self, t, x): + """Override the base class.""" + # Unpack the states for better readability + 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 + + 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 + i_dc = self.converter.dc_current(self.q, i_ss) + + # Return the list of state derivatives + return ( + machine_f + self.mechanics.f(t, w_M, tau_M) + + self.converter.f(t, u_dc, i_L, i_dc)) + + def save(self, sol): + """Extend the base class.""" + super().save(sol) + self.data.u_dc.extend(sol.y[4].real) + self.data.i_L.extend(sol.y[5].real) + + def post_process(self): + """Extend the base class.""" + super().post_process() + # From lists to the ndarray + 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.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) + # 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 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 + self.data.i_g = abc2complex(q_g_abc)*self.data.i_L + + # %% class DriveTwoMassMechanics(Drive): """ Synchronous machine drive with two-mass mechanics. - This interconnects the subsystems of a synchronous machine drive and - provides an interface to the solver. - Parameters ---------- machine : SynchronousMachine @@ -356,15 +412,14 @@ def set_initial_values(self, t0, x0): """Extend the base class.""" super().set_initial_values(t0, x0[0:4]) self.mechanics.w_L0 = x0[4].real - self.mechanics.theta_ML0 = np.mod(x0[5].real + np.pi, 2*np.pi) - np.pi + self.mechanics.theta_ML0 = wrap(x0[5].real) def f(self, t, x): """Override the base class.""" # 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.converter.q, self.converter.u_dc0) + 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 # State derivatives plus the outputs for interconnections machine_f, _, tau_M = self.machine.f(psi_s, u_s, w_M)