FOC Simulator: PMSM and Induction Motor Current Control

The amplitude-invariant dq convention uses peak phase quantities. For the PMSM, the plant follows the synchronous-machine voltage equation in Tiitinen et al., IEMDC 2023, attached p. 3, Eqs. (5)–(6):

For the induction machine, transform the linear Γ model in the same paper (attached pp. 2–3 and 6–7) to its equivalent inverse-Γ representation. Rotor-flux orientation gives:

The frame angle integrates $\omega_s$ (or $\omega_e$ for PMSM). An adaptive Dormand–Prince RK45 solver integrates both currents, both PI integrator voltages, rotor flux and frame angle. Its relative tolerance is 10⁻⁶. The 201 output times concentrate around the command step; the phase plot is reconstructed separately over its final two periods.

PI tuning. With $\omega_c=2\pi f_{bw}$, use $K_{pd}=\omega_c L_{ed}$, $K_{pq}=\omega_c L_{eq}$ and $K_i=\omega_c R_{eff}$. PMSM uses its own axis inductances and $R_{eff}=R_s$. The induction model uses $L_{ed}=L_{eq}=L_\sigma$ and $R_{eff}=R_s+R_R$. Feedforward cancels the cross terms, leaving these RL plants. The induction feedforward is $[-R_R\psi_R/L_M-\omega_sL_\sigma i_q,\ \omega_sL_\sigma i_d+\omega_e\psi_R]$. With nominal multipliers and no saturation, pole-zero cancellation yields $i/i^*=\omega_c/(s+\omega_c)$ and a 10–90% rise time $\ln(9)/\omega_c$. These tuning relations are derived here; they are not numerical gains copied from SPRA524.

Voltage limiting and anti-windup. The average voltage vector is limited radially to $V_{dc}/\sqrt3$ for SVPWM or $V_{dc}/2$ for sinusoidal PWM, consistent with the linear modulation geometry described in TI SPRA524 (1999). The integrator back-calculation is $\dot z=K_i(i^*-i)+\omega_c(v_{limited}-v_{raw})$ when enabled. AN1078 (2010), p. 8, describes the equivalent excess-voltage feedback principle.

References and weakening. The PMSM MTPA direction follows stationary torque on a constant-current circle, following Morimoto et al. (1990), DOI:10.1109/28.60058. A scalar root finds the current magnitude for the requested torque.

PMSM weakening uses a continuous feasibility interval rather than a fixed current grid. With conventional saliency, positive speed, $i_d\le0$ and $i_q\ge0$, squared voltage increases with $i_q$. At $i_q=0$:

Complete the square with $a=R_s^2+(\omega_eL_d)^2$, center $d_0=-\omega_e^2L_d\psi_f/a$, and minimum squared voltage $U_{min}^2=(R_s\psi_f\omega_e)^2/a$. Intersect

with $[-I_{max},0]$. A negative radicand or empty intersection proves infeasibility. This interval can be much narrower than any fixed d-current grid.

For each d current, the positive q-current limit is the smaller of the current-circle bound and the resistive voltage-ellipse root. Bounded scalar optimization finds the maximum feasible torque; if the requested torque is lower, bracketed roots locate the d-current interval that can deliver it. Projecting the unconstrained MTPA d current onto that interval minimizes current. A cancellation-resistant quadratic-root expression is used for small q-current limits.

For induction weakening, first restrict d current to

This analytical no-load voltage bound gives a continuously feasible positive-flux domain. Search 180 flux-current levels within that domain, with 20 bisections for each q-current voltage limit. The induction optimum remains discretized, but a narrow feasible domain is not skipped by sampling the entire nominal flux range.

The energy audit independently compares integrated terminal power with stator and rotor copper heat, mechanical conversion, and the change in stored magnetic energy. Phase-current balance, inverse-Park magnitude and the MTPA stationarity equation are checked separately.

Responsiveness before integration. The dimensionless workload estimate is duration multiplied by a conservative sum of rates: the largest $R_{eff}/L$ pole, induction rotor-flux decay and current-limit slip estimate, electrical rotation, and $\omega_c(g_p+\sqrt{g_i})$ for the PI controller. The authoring budget is 1200; it is a responsiveness limit, not an integration-error estimate. The guard branches before calling ode45. Unsupported combinations produce an explanatory input error and withhold results, so a high $R/L$ ratio cannot silently start an expensive explicit solve.

Independent transient reproduction. The original attached harness independent_transients.py (source s7) calls unmodified motulator 0.3.0 SynchronousMachine and stationary-frame Γ InductionMachine classes (s6). Plant states are fluxes; the calc instead integrates dq currents. An external wrapper supplies the same documented continuous PI law and exact measured field angle, while motulator computes plant derivatives, measured currents and torque. It reads no calc file or exported calc trajectory.

Install motulator==0.3.0 and SciPy==1.15.3. Run the attached Python harness; with local dependencies the command is MPLCONFIGDIR=./mplconfig PYTHONPATH=reference_packages python3 independent_transients.py. The PMSM uses p=3, Rs=3.6 Ω, Ld=36 mH, Lq=51 mH, PM flux=0.545 Wb and 1000 rpm. The IM uses p=2, Rs=3.7 Ω, inverse-Γ RR=2.1 Ω, LM=224 mH, Lσ=21 mH, rotor flux=0.9 Wb and 750 rpm; its Γ transformation is performed by motulator. Both use 540 V DC, 200 Hz bandwidth and a q-current step from zero to 0.8 A at 5 ms, with pre-established d current and flux. Equivalent torque inputs are read from the motor models' torque properties, not used as the sole expected output.

SciPy solve_ivp uses rtol=10⁻¹⁰, atol=10⁻¹², max_step=20 μs and separate solves before and after the step. The saved record s8 contains currents, stator and rotor flux, voltage and torque at 5.954929658551372 ms and 20 ms. The additional voltage-limited PMSM record s9 uses the same harness with Q_COMMAND=5 A and run_case("pmsm"); all other settings stay the same. These reproduce published model implementations; the numbers are not claimed to be printed in the IEMDC paper.

• Balanced three-phase sinusoidal machine; constant parameters, ideal rotor angle and speed feedback, fixed DC bus. All fluxes and dq currents use peak values.
• Positive speed and positive motoring torque only. Surface PMSM ($L_d=L_q$) and conventional interior PMSM ($L_q\ge L_d$) are supported in MTPA/weakening mode.
• Shaft speed is imposed by an ideal external drive or dynamometer. There is no mechanical acceleration or outer speed-control loop.
• Before the torque step, $i_d$ is already at its final reference, $i_q=0$, and the PI states hold that equilibrium. For an induction motor $\psi_R=L_Mi_d$ is already established. The simulation does not include initial magnetization or the transition between flux operating points.
• Induction parameters use the inverse-Γ circuit. For a linear Γ model, let $g=L_s/(L_s+L_\ell)$; then $L_M=gL_s$, $L_\sigma=gL_\ell$, $R_R=g^2R_r$.
• Equations, controller implementation, constraint search, sampling layout and drawings are original work derived here from the cited models. No source plots or tables are copied. The motulator project is MIT licensed; its examples provide published external-load commands for tracking checks; independent electrical outputs are reproduced separately with its unmodified machine models.
• Current and voltage constraints apply to the reference generator and averaged inverter respectively. Actual transient current is allowed to exceed the commanded limit and triggers a warning.

Use for current-loop concepts, initial tuning and comparisons of voltage/current constraints. The model omits digital sampling and computational delay, PWM ripple and harmonics, overmodulation, inverter dead time, bus ripple, magnetic saturation and cross-saturation, parameter error, sensorless estimation, iron loss, friction and thermal effects. Gain multipliers change the dynamics away from the nominal bandwidth design.

PMSM field weakening uses continuous feasibility bounds, bracketed roots and bounded scalar optimization. It is not a demagnetization assessment. Induction weakening uses 180 current levels within its analytically feasible positive-flux interval; its optimum is a grid approximation. A point on the voltage boundary may settle slowly because little dynamic voltage remains. In fixed-flux or zero-d mode, the controller can saturate rather than meet the requested torque.

A workload guard includes electrical decay rates as well as speed and controller dynamics, and rejects excessive explicit-solver work before integration. For unsupported high-resistance/low-inductance combinations, shorten the run or use a simulator with a stiff-system integrator. The guard is deliberately conservative and does not imply that a rejected motor operating point is physically impossible.

The 14 N·m PMSM and 14.6 N·m induction cases are tracking checks based on published external-load commands, not independent electrical output benchmarks. Separate cases compare an MTPA point and PMSM/induction current, flux, voltage and torque transients obtained from unmodified motulator 0.3.0 machine models. These use an original matching continuous-PI wrapper, so they do not establish agreement with the complete sensorless sampled-data drive examples. The narrow-ellipse and stiffness cases are numerical regressions supplied by the independent review.

The energy residual uses trapezoidal integration of saved samples, so its small nonzero value is a quadrature error. Phase waveforms are interpolated from dq states and show up to two final electrical periods. Replay changes only interpolated quantities and the drawing; the ODE and verification do not depend on replay time. Current-vector length is normalized by the larger of the current limit and simulated peak, while the d/q arrow lengths are schematic.