Physically Derived Weight for an Exact Torque Penalty in Loss-Minimizing Model Predictive Control of Induction Motor Drives

1. Introduction

An induction motor (IM) driving a pump, fan or conveyor spends most of its operating hours well below rated torque. At rated flux, such a drive keeps paying for magnetizing current and core loss that the load does not need. Loss-minimizing control (LMC) lowers the rotor flux until the extra copper loss of the torque-producing current outweighs the saving in magnetizing and core loss [1,2]. Under the extended product approach of IEC 61800-9-2, the converter losses enter the same balance, because the power drive system is assessed as a whole [3].

Model predictive control (MPC) handles current and voltage limits explicitly and accepts objectives other than tracking errors, which has made it a common framework for AC drives [4,5,6,7]. Economic MPC (EMPC) minimizes the operating cost itself, here the electrical loss [8,9]. Loss-aware predictive torque and current controllers have been reported for traction and linear induction machines [10,11,12], and the efficiency-oriented MPC literature for AC motors has been surveyed recently [13]. In previous work, the authors proposed a thermally adaptive EMPC (TA-EMPC) for a 7.5 kW industrial drive [14]. It minimizes a temperature-coupled loss model of the machine and inverter, adapts the stator and rotor resistances with a sensorless thermal observer, and reserves flux ahead of predicted load rises.

With loss as the objective, the torque reference turns into a constraint, and how it is imposed decides what the drive actually delivers. A hard equality on the predicted torque can make the quadratic program (QP) infeasible while the flux is being rebuilt or the current limit is reached. A quadratic penalty avoids infeasibility but turns the constraint into a trade-off. Each newton-meter costs watts, so the optimizer settles slightly below the reference. In TA-EMPC the torque weight was set so that this error stayed below 0.5% at the rated point. The weight sweep in [14] showed the error falling tenfold for every tenfold increase of the weight (4.77%, 0.50%, 0.05% and 0.005% at 0.1, 1, 10 and 100 times the nominal weight). This scaling was reported but not explained, and no rule was given for other ratings or thermal states.

The choice of weighting factors has a long history in predictive drive control. Proposed remedies include empirical guidelines [15,16], imposed optimized weights [17], ranking and entropy methods [18,19] and sequential formulations without weights [20]. All of them were developed for finite-control-set MPC, where the competing terms are tracking errors of a similar kind. A systematic review of energy-efficient MPC for industrial IM drives lists weighting under an economic objective, where the terms have a physical and monetary meaning, among the open problems [21]. Optimization theory supplies two relevant tools: the exact-penalty property of the ℓ1 norm, used in MPC to soften constraints [22,23,24], and the interpretation of a Lagrange multiplier as the sensitivity of the optimal cost to its constraint [25]. To the authors’ knowledge, they have not yet been used to assign a physical value to the torque weight of a loss-minimizing drive controller.

In this paper the torque term of the EMPC is replaced by a penalty proportional to the absolute torque error. Its weight μ_M is expressed in W/(N·m) and is bounded from below by the marginal loss per unit torque λ_M of the loss-optimal steady state. The contributions are:

• a closed-form expression for the steady torque error of the quadratic penalty, ε_M = λ_M/(2q_MM^*), which accounts for the 1/q_M scaling in [14] and predicts that a fixed weight loses accuracy at light load and at high winding temperature;
• a sufficient condition ${\mu }_M>\Lambda \left(M^{*}\right)={\max }_{0\le M\le M^{*}}{\lambda }_M\left(M\right)$ for zero steady error with the absolute-value penalty, the local condition μ_M > λ_M(M^*), and the form both results take in a controller that holds the flux at its steady-state target, together with an analytic expression and a torque–speed map of λ_M for a 7.5 kW drive at 20 °C and 115 °C;
• a numerical check on the receding-horizon controller, by locating the smallest weight that gives zero steady error in closed loop at 36 operating points, including rated speed with an active voltage limit;
• a simulation study covering current-limited transients with and without the flux reservation of [14], two duty cycles, uncertainty of the loss model, field-orientation error caused by mismatched rotor parameters, current-measurement noise and inverter dead time.

Section 2 describes the models, the controller and the simulation setup. Section 3 presents the results. Section 4 discusses them and states the limitations, and Section 5 concludes. The simulation code and data are provided as described in the Data Availability Statement.

2. Materials and Methods

2.1. Machine Model

The cage IM is modeled in the rotor-flux-oriented d–q frame with the amplitude-invariant transformation. Index 1 refers to the stator and index 2 to the rotor. The rotor flux linkage Ψ₂ lies on the d-axis; the state vector is $\mathbf{x}={\left[i_{1d},\hspace{0.17em}{i}_{1q},\hspace{0.17em}{\Psi }_2\right]}^{\mathrm{T}}$ and the input is $\mathbf{u}={\left[u_{1d},\hspace{0.17em}{u}_{1q}\right]}^{\mathrm{T}}$:

$${\displaystyle \frac{d{\Psi }_2}{\mathit{dt}}}={\displaystyle \frac{R_2}{L_2}}\left(L_m{i}_{1d}-{\Psi }_2\right),$$

(1)

$$\sigma L_1{\displaystyle \frac{d{i}_{1d}}{\mathit{dt}}}=u_{1d}-R_{\sigma }{i}_{1d}+{\omega }_1\sigma L_1{i}_{1q}+k_2{\displaystyle \frac{R_2}{L_2}}{\Psi }_2,$$

(2)

$$\sigma L_1{\displaystyle \frac{d{i}_{1q}}{\mathit{dt}}}=u_{1q}-R_{\sigma }{i}_{1q}-{\omega }_1\sigma L_1{i}_{1d}-k_2{z}_p\omega {\Psi }_2.$$

(3)

Here R₁ and R₂ are the stator and rotor resistances, L₁, L₂ and L_m the stator, rotor and magnetizing inductances, k₂ = L_m/L₂, $\sigma =1-L_m^2/\left(L_1{L}_2\right)$, $R_{\sigma }=R_1+k_2^2{R}_2$, z_p the number of pole pairs, ω the mechanical speed and ω₁ = z_pω + k₂R₂i_1q/Ψ₂ the angular frequency of the stator field. The electromagnetic torque is

$$M={\displaystyle \frac{3}{2}}{z}_p{k}_2{\Psi }_2{i}_{1q}.$$

(4)

By (4), the same torque can be obtained with a high flux and a small i_1q or with a low flux and a large i_1q; LMC exploits this freedom.

2.2. Loss Model and Parameters

The controllable loss of the drive is taken from [14]:

$$\Delta P_{\Sigma }=\Delta P_{\mathit{Cu}1}+\Delta P_{\mathit{Cu}2}+\Delta P_{\mathit{Fe}}+\Delta P_{\mathit{inv}},$$

(5)

$$\Delta P_{\mathit{Cu}1}={\displaystyle \frac{3}{2}}{R}_1\left({\vartheta }_1\right)\left(i_{1d}^2+i_{1q}^2\right), \Delta P_{\mathit{Cu}2}={\displaystyle \frac{3}{2}}{R}_2\left({\vartheta }_2\right)\left(i_{2d}^2+i_{2q}^2\right),$$

(6)

$$\Delta P_{\mathit{Fe}}={\displaystyle \frac{3}{2}}\left(k_h\left|{\omega }_1\right|+k_e{\omega }_1^2\right){\Psi }_2^2, \Delta P_{\mathit{inv}}=a_I{I}_{1m}+b_I{I}_{1m}^2,$$

(7)

where i_2d = (Ψ₂ − L_mi_1d)/L₂ and i_2q = −k₂i_1q are the rotor current components and $I_{1m}=\sqrt{i_{1d}^2+i_{1q}^2}$ is the stator current amplitude. The coefficient a_I = 3(U_CE0 + U_F0)/π + (6/π)f_swE_sw(U_d/U_ref)/I_ref collects the threshold-voltage conduction loss and the switching loss, and b_I = 3(r_CE + r_F)/4 the slope-resistance conduction loss, following Equations (11) and (12) of [14]; with the data of Table 1, a_I = 6.30 W/A and b_I = 0.0338 Ω. The resistances vary as R_j(ϑ_j) = R_j,20[1 + α_j(ϑ_j − 20)], j = 1,2.

Table 1 lists the parameters and their origin. Machine and inverter data are those of Tables 2 and 3 in [14]. The quantities marked “this study” are either not given there or were changed. The switching frequency was raised from 8 kHz in [14] to 10 kHz, so that one PWM period equals the sampling period of 100 µs and the controller updates the voltage reference once per carrier period. The temperature coefficient of the aluminum cage, α₂ = 0.00408 K⁻¹, corresponds to the temperature constant of 225 K for aluminum in IEC 60034-2-1 [26]. The rated rotor flux linkage was set to Ψ_2N = 0.90 Wb; with this value the stator voltage at the rated point, 311.7 V, just fits the linear modulation limit $U_d/\sqrt{3}$ = 311.8 V. The rated torque, speed and current then follow from the model. The input blocking, the weights r and ρ_Ψ and the flux envelope enter the controller of Section 2.4. The parameter set is representative of a 7.5 kW, 400 V, four-pole machine and is not a nameplate record of a specific motor.

2.3. Steady-State Economic Layer and the Marginal Loss per Unit Torque

For a torque reference M^*, speed ω and winding temperatures ϑ = (ϑ₁,ϑ₂), the loss-optimal steady state is the solution of

$$P^{*}\left(M^{*},\omega ,\vartheta \right)=\underset{{\Psi }_2}{\min }\hspace{0.25em}\Delta P_{\Sigma }^{\mathit{ss}}\left({\Psi }_2\right) \mathrm{s}.\mathrm{t}. M=M^{*},\hspace{0.25em}{\mathbf{i}}_1\in {\mathcal{P}}_I,\hspace{0.25em}{U}_{1m}\le U_d/\sqrt{3},\hspace{0.25em}{\Psi }_{\mathit{min}}\le {\Psi }_2\le {\Psi }_{2N},$$

(8)

where ${\mathcal{P}}_I$ is the 16-sided polygon inscribed in the circle of radius I_max that the dynamic controller also uses (Section 2.4), so that every steady-state target is admissible for it; for the voltage the circle is used, which lies inside the voltage polygon of the controller. The steady-state relations are i_1d = Ψ₂/L_m, i_1q = M^*/(1.5z_pk₂Ψ₂), u_1d = R₁i_1d − ω₁σL₁i_1q and u_1q = R₁i_1q + ω₁L₁i_1d, and $U_{1m}=\sqrt{u_{1d}^2+u_{1q}^2}$. The minimizer Ψ_2s serves as the flux target of the dynamic controller. The central quantity of this paper is the derivative of the optimal loss with respect to the torque reference,

$${\lambda }_M\left(M^{*},\omega ,\vartheta \right)={\displaystyle \frac{\partial P^{*}}{\partial M^{*}}} \left[\mathrm{W}/\left(\mathrm{N}\cdot \mathrm{m}\right)\right].$$

(9)

By the sensitivity theorem of nonlinear programming [25], λ_M is the Lagrange multiplier of the equality M = M^* in (8). It states how many watts of loss the last newton-meter of torque costs. While the current and voltage limits are inactive, the first-order conditions of (8) yield

$$\begin{gathered} {\lambda }_M={\displaystyle \frac{3R_{\sigma }{i}_{1q}+\left(a_I/I_{1m}+2b_I\right)i_{1q}+\frac{\partial \Delta P_{\mathit{Fe}}}{\partial {\omega }_1}\frac{k_2{R}_2}{{\Psi }_{2s}}}{\frac{3}{2}{z}_p{k}_2{\Psi }_{2s}}} \\ \approx {\displaystyle \frac{4M^{*}}{9z_p^2{k}_2^2{\Psi }_{2s}^2}}\left(3R_{\sigma }+{\displaystyle \frac{a_I}{I_{1m}}}+2b_I\right). \end{gathered}$$

(10)

The approximation neglects the dependence of the core loss on the slip frequency. Equation (10) shows how λ_M behaves. It rises with the winding resistances, and hence with temperature, and falls with the square of the flux. Under LMC the optimal flux grows roughly as $\sqrt{M^{*}}$ at partial load, so λ_M stays nearly constant there. Once the flux reaches its upper bound, λ_M grows in proportion to the torque. If the current or voltage limit is active, the multiplier also carries the price of that limit, and (10) no longer applies. In general, λ_M is therefore computed by differences of P^* with a step of ±0.05 N·m, central inside the feasible region and one-sided (backward) where M^* + 0.05 N·m is infeasible. The steady-state layer stores Ψ_2s and λ_M in look-up tables (LUTs) on a 53 × 31 grid over 0 ≤ M^* ≤ 1.3M_N and 0.05ω_N ≤ ω ≤ ω_N, one pair of tables per thermal state; beyond M_max the value at the last feasible grid point is held.

2.4. Dynamic Economic MPC

Over the prediction horizon, ω and ω₁ are frozen at their present values, which makes (1)–(3) linear. Exact zero-order-hold discretization with T_s = 100 µs gives x_i+1 = A_dx_i + B_du_i. The input is held over blocks of 1, 1 and 8 steps, so N_c = 3 input vectors $\mathbf{U}={\left[{\mathbf{u}}_0^{\mathrm{T}},{\mathbf{u}}_1^{\mathrm{T}},{\mathbf{u}}_2^{\mathrm{T}}\right]}^{\mathrm{T}}$ cover the horizon N = 10, and condensing gives X = f₀ + GU. At each step the controller solves

$$\underset{\mathbf{U}}{\min }\hspace{0.25em}{\displaystyle \sum _{i=1}^N}\left[\Delta P_{\Sigma }\left({\mathbf{x}}_i\right)T_s+\varphi \left(e_i\right)T_s\right]+r{\displaystyle \sum _{j=0}^{N_c-1}}{∥{\mathbf{u}}_j-{\mathbf{u}}_{j-1}∥}^2+{\rho }_{\Psi }{\left({\Psi }_{2,N}-{\Psi }_{2s}\right)}^2$$

(11)

subject to the predicted dynamics; a 16-sided polygon circumscribing the circle $\parallel {\mathbf{u}}_j\parallel \le U_d/\sqrt{3}$, which admits at most 1/cos(π/16) − 1 ≈ 1.9% overmodulation; a 16-sided polygon inscribed in the circle ||i_1,i||≤I_max at the end of each block; and Ψ_min ≤ Ψ_2,i ≤ Ψ_2N at every step. The predicted torque error is e_i = M^* − M_i, where M_i is linearized about the shifted solution of the previous step, and ϕ(⋅) is the torque term of Section 2.5. The loss terms (6)–(7) are quadratic in x, except a_II_1m, which is replaced by its first-order expansion at the measured current. Since the loss is multiplied by T_s, the cost is an energy and ϕ must be given in watts. In (11), u_j−1 for j = 0 is the input applied in the previous sampling period.

The horizon covers NT_s = 1 ms, whereas the rotor time constant L₂/R₂ is 182 ms at 20 °C and 131 ms at 115 °C. Even with the full current I_max applied to the d-axis, the flux can change within one horizon by no more than 14 mWb, about 5% of the lower flux bound, so the flux trajectory is governed by the terminal term in (11), which pulls Ψ_2,N toward the loss-optimal target. Its weight ρ_Ψ was chosen so that the steady flux stays close to Ψ_2s; at the simulated operating points the deviation was at most 0.23%. As an option, the target can be raised in advance of a load increase by the flux reservation of [14]: Ψ_2s is evaluated at $\max \{M^{*}\left(t\right),{\widehat{M}}_{\mathit{max}}\left(t\right)\}$, where ${\widehat{M}}_{\mathit{max}}$ is the largest reference over a preview window of 100 ms. This is a simplified form of Equation (17) in [14] with the conservativeness factor set to 1. The QP is solved by the dual active-set method of Goldfarb and Idnani [27]. Figure 1 shows the controller structure.

2.5. Torque Term

Four forms of ϕ are compared. In the hard formulation (H), e_i = 0 is imposed at the end of each input block, which gives N_c equality constraints. If this QP is infeasible, the step is solved with the quadratic formulation instead, the usual safeguard in practice. The quadratic formulation (Q) uses ϕ_Q(e) = q_Me² with q_M in W/(N·m)², as in [14]. The absolute-value formulation (L1) uses ϕ_L1(e) = μ_M|e| with μ_M in W/(N·m). It is implemented with one slack variable s_j per input block ${\mathcal{B}}_j$:

$$-s_j\le e_i\le s_j \forall i\in {\mathcal{B}}_j, s_j\ge 0, \mathrm{penalty}{\mu }_M{T}_s{\displaystyle \sum _{j=0}^{N_c-1}}{n}_j{s}_j,$$

(12)

where n_j is the number of steps in block j. Strictly, (12) is a block-wise majorant of the ℓ1 sum ${\mu }_M{T}_s\sum _i\left|e_i\right|$: inside a block it penalizes the largest error, as an ℓ∞ norm would. The two coincide when the error is constant within each block. In steady state the error of the first block, which determines the applied input, is zero; residual errors planned for the later blocks are discussed in Section 3.2. Relative to Q, the QP gains N_c variables and 2N + N_c inequality rows. The weight is either constant, a fixed multiple of ${\bar{\lambda }}_M$, the largest λ_M over the operating envelope in the hot state (Section 2.6), or scheduled, ${\mu }_M\left(k\right)=\kappa {\lambda }_M\left(M^{*}\left(k\right),\omega \left(k\right),\widehat{\vartheta }\left(k\right)\right)$, from the LUT of Section 2.3 (L1-S). The fourth formulation (L1+Q) adds both terms, ϕ(e) = μ_M|e| + q_Me². Since the derivative of q_Me² is zero at e = 0, the quadratic part does not shift the exactness threshold; it only raises the price of large transient errors.

The analysis uses the reduced steady-state problem, in which the flux is re-optimized for the delivered torque M and M is the only decision variable. Let M_max(ω,ϑ) denote the largest torque for which (8) is feasible.

Proposition 1.

(bias of the quadratic penalty). Let P^* be differentiable at the minimizer M_Q of$P^{*}\left(M\right)+q_M{\left(M^{*}-M\right)}^2$

for 0 ≤ M ≤ M_max, and let 0 < M_Q < M_max. Then M^* − M_Q = λ_M(M_Q)/(2q_M), and the relative steady error is ε_M = λ_M(M_Q)/(2q_MM^*).

Proof. 

At an interior minimizer the derivative vanishes: λ_M(M_Q) − 2q_M(M^* − M_Q) = 0. □

Proposition 2.

(zero steady error with the absolute-value penalty). Let 0 < M^* ≤ M_max, let P^* be absolutely continuous for 0 ≤ M ≤ M_max with derivative λ_M almost everywhere, and let λ_M ≥ 0 almost everywhere for M^* ≤ M ≤ M_max. If μ_M > Λ(M^*), where Λ(M^*) is the essential supremum of λ_M for 0 ≤ M ≤ M^*, then M = M^* is the unique global minimizer of P^*(M) + μ_M|M^* − M| for 0 ≤ M ≤ M_max.

Proof. 

For M^* < M ≤ M_max, $P^{*}\left(M\right)-P^{*}\left(M^{*}\right)={\int }_{M^{*}}^M{\lambda }_M\hspace{0.17em}\mathit{dm}\ge 0$

, so the cost exceeds P^*(M^*) by at least μ_M(M − M^*) > 0. For 0 ≤ M < M^*, $P^{*}\left(M^{*}\right)-P^{*}\left(M\right)={\int }_M^{M^{*}}{\lambda }_M\hspace{0.17em}\mathit{dm}\le \Lambda \left(M^{*}\right)\left(M^{*}-M\right)$

, hence P^*(M) + μ_M(M^* − M) ≥ P^*(M^*) + (μ_M − Λ)(M^* − M) > P^*(M^*). □

Proposition 2 is global on the admissible torque interval and needs neither convexity nor monotonicity of λ_M; the only assumptions are that P^* does not decrease above the reference and that λ_M is essentially bounded for torques below it. Both hold on the envelope of this drive, where P^* is continuous and piecewise differentiable, with kinks where a current or voltage limit becomes active. The condition is sufficient, not necessary. Since the cost P^*(M) + μ_M(M^* − M) has the left derivative λ_M(M^*) − μ_M at M^*, the weaker condition μ_M > λ_M(M^*) already makes M^* a strict local minimizer, and μ_M ≥ λ_M(M^*) is necessary for any local minimum at M^*. The two bounds differ where λ_M decreases with torque, which happens at partial load when the optimal flux leaves its lower bound. On a grid with steps of 0.005M_N and 0.05ω_N over the operating envelope defined in Section 2.6, Λ(M^*)/λ_M(M^*) reaches 1.277 (at M^*≈ 0.49M_N in the cold state); the scheduled weight therefore uses κ = 1.35, which satisfies the global condition everywhere.

Both propositions refer to the reduced problem, in which the flux follows the delivered torque. The dynamic controller differs in one respect: its terminal term holds the flux at the target Ψ_2s(M^*) of the reference. The marginal loss that it trades against the torque term is therefore ${\lambda }_M^{\Psi }\left(M\right)=\partial \Delta P_{\Sigma }^{\mathit{ss}}/\partial M$ at Ψ₂ = Ψ_2s(M^*), and the arguments of both propositions apply with λ_M replaced by ${\lambda }_M^{\Psi }$; in particular, Proposition 1 becomes $M^{*}-M_Q={\lambda }_M^{\Psi }\left(M_Q\right)/\left(2q_M\right)$. By the envelope theorem, ${\lambda }_M^{\Psi }\left(M^{*}\right)={\lambda }_M\left(M^{*}\right)$ wherever no current or voltage limit is active in (8). Away from M^* the two differ, because at fixed flux the copper loss makes ${\lambda }_M^{\Psi }$ nearly proportional to the torque, whereas λ_M includes the re-optimization of the flux. Where the voltage limit is active in (8), λ_M also contains the price of that limit, which does not enter ${\lambda }_M^{\Psi }$, so that ${\lambda }_M^{\Psi }\left(M^{*}\right)<{\lambda }_M\left(M^{*}\right)$ there.

The propositions describe the steady state of the reduced problem. Whether the receding-horizon QP, with its 1 ms horizon, linearization, input blocking and terminal term, has the same threshold is not implied by them. For a convex QP, the exact-penalty theorem [23,24] guarantees a finite threshold, equal to a norm of the multipliers of the torque constraints, but gives no value for it. Moreover, the controller optimizes over a short horizon around the present state and can only respond to the local condition. The threshold of the closed loop, denoted ${\mu }_{\mathit{cl}}^{*}$, was therefore determined directly. For each operating point, closed-loop simulations were repeated with different μ_M, and the smallest weight for which the steady error fell below 10⁻⁴% was located by bisection to a relative resolution of 1.4 × 10⁻⁵ (Section 3.2).

If M^* > M_max, the hard formulation has no solution. The absolute-value formulation remains feasible, but for a finite μ_M nothing in the propositions ensures that it delivers M_max; its behavior in this case is examined numerically in Section 3.3. A simple remedy is to limit the reference to $M_{\mathit{max}}\left(\omega ,\widehat{\vartheta }\right)$, which the steady-state layer computes together with Ψ_2s and λ_M; Proposition 2 then applies to the limited reference. This option was not used in the simulations.

2.6. Simulation Setup and Metrics

The plant integrates the nonlinear model (1)–(4), with ω₁ updated continuously, by the fourth-order Runge–Kutta method with five substeps per control period and zero-order-hold inputs; with twenty substeps, the torque in the sequence of Figure 4 changed by less than 10⁻⁸% of M_N and the loss energy by less than 10⁻¹¹%. The inverter is represented by its average model; dead time is added only in scenario (vii). Except in scenarios (vi) and (vii), field orientation is ideal and the controller receives i_1d, i_1q and Ψ₂ without error. The load imposes the speed (dynamometer mode), so the torque controller is tested without a speed loop; this is an idealization, and the interaction with a speed controller is not examined. Two thermal states are used: cold, ϑ₁ = ϑ₂ = 20 °C, and hot, ϑ₁ = ϑ₂ = 115 °C, the hot steady state of [14], at which R₁ is 37.3% and R₂ 38.8% above the cold values. The thermal observer of [14] is represented by the true temperatures.

The quadratic baseline uses the weight q_M,nom = λ_M/(2 ⋅ 0.005M_N) = 52.1 W/(N·m)², evaluated at the cold rated point, i.e., the weight for which the first-order estimate λ_M(M^*)/(2q_MM^*) of Proposition 1 gives 0.5% error at the rated point. This reproduces the selection rule stated in [14]; the absolute weight of [14] is not reported there. The operating envelope is 0 < M^* ≤ min{1.25M_N,M_max(ω)}, 0.05ω_N ≤ ω ≤ ω_N. The constant absolute-value weight is ${\mu }_M=1.2{\bar{\lambda }}_M$ = 68.0 W/(N·m), where ${\bar{\lambda }}_M$ = 56.7 W/(N·m) is the maximum of λ_M over this envelope in the hot state, evaluated on the grid of Section 2.5 and on the boundary M^* = M_max(ω). It is reached at rated speed on the voltage-limit boundary, at M^* = 1.227M_N; at rated speed and 115 °C the voltage limit bounds the torque to M_max = 1.228M_N. The dimensionless factor 1.2 is a margin for revisions of the loss model (Section 3.5). The scheduled weight uses κ = 1.35. Table 2 summarizes the weights used throughout the paper.

Eight scenarios are simulated: (i) steady operation at 18 operating points per thermal state, i.e., 36 points in total, with M^*/M_N ∈ {0.1, 0.25, 0.5, 0.75, 1.0, 1.25} at ω/ω_N ∈ {0.5, 0.9} and M^*/M_N ∈ {0.1, 0.25, 0.5, 0.75, 1.0, 1.2} at ω = ω_N, where the voltage limit is active at the highest torque, and weight sweeps at two further points, 0.5M_N at 0.75ω_N and 0.2M_N at 0.5ω_N; (ii) a torque sequence at 0.75ω_N in the hot state, from 0.2M_N to 0.9M_N, then to 1.4M_N, which exceeds the current limit, and back to 0.5M_N, with and without flux reservation and for horizons N = 10, 20 and 40; (iii) a pump cycle following $M^{*}=M_N{\left(\omega /{\omega }_N\right)}^2$ with speed set-points between 0.45ω_N and 0.9ω_N (12 s); (iv) an intermittent conveyor cycle at 0.75ω_N with ten torque levels between 0.15M_N and 0.9M_N (10 s), with and without flux reservation; (v) a Monte Carlo study of loss-model uncertainty; (vi) a study of field-orientation error, in which the controller and its current-model flux estimator use perturbed values of R₂ and L_m while the plant keeps the nominal ones; (vii) steady operation with current-measurement noise and uncompensated inverter dead time; and (viii) a numerical check of closed-loop convergence, by linearizing the one-step closed-loop map at each equilibrium and by simulating from perturbed initial states.

The metrics are the steady relative torque error ε_M = (M^* − M)/M^*, averaged over each plateau after 0.3 s of settling; the relative integral absolute error ${\mathrm{IAE}}_{\mathit{rel}}=\int \left|M^{*}-M\right|\hspace{0.17em}\mathit{dt}/\int M^{*}\mathit{dt}$; the loss energy $E_{\mathit{loss}}=\int \Delta P_{\Sigma }\hspace{0.17em}\mathit{dt}$; the delivered mechanical work $W=\int M\omega \hspace{0.17em}\mathit{dt}$; and the ratio E_loss/W. The last metric avoids crediting a controller with a saving that comes only from delivering less work. The simulations were written in Python 3.13 with NumPy 2.5 and SciPy 1.18; the QP solver is the quadprog implementation of the Goldfarb–Idnani method. To check reproducibility, the complete computation was repeated from Supplementary Material S1 in a newly created Python environment (NumPy 2.5.4, SciPy 1.18.1, quadprog 0.1.14); all values derived from the results, 30,806 in total, agreed with the stored ones to a relative tolerance of 10⁻⁶, except the measured QP solution times. Generative AI was used in this study as stated in the Acknowledgments.

3. Results

3.1. Marginal Loss Per Unit Torque

Figure 2 maps λ_M over the torque–speed plane together with contours of the loss-optimal flux. For 0.1M_N ≤ M^* ≤ min{1.25M_N,M_max}, λ_M lies between 13.1 and 37.6 W/(N·m) at 20 °C and between 17.1 and 56.7 W/(N·m) at 115 °C. Up to about 0.5M_N the optimal flux is below its rated value and λ_M changes little (13.1–17.4 W/(N·m) cold, 17.1–21.5 W/(N·m) hot), in line with (10). In this range λ_M is largest at 0.1M_N and decreases slightly with torque; the decrease begins where the optimal flux leaves its lower bound Ψ_min, at about 0.06M_N, i.e., below the range considered here. Above about 0.5M_N the flux sits at Ψ_2N, and λ_M grows in proportion to the torque. Near rated speed the voltage limit takes over at high torque; at ω_N it becomes active above 1.01M_N in the cold state and above 0.79M_N in the hot state. The optimal flux then drops below Ψ_2N and λ_M climbs steeply, up to ${\bar{\lambda }}_M$ = 56.7 W/(N·m) in the hot state. At the operating points of Table 4 where no limit is active, heating from 20 °C to 115 °C raises λ_M by 24–35%. The increase follows R_σ in (10) and is diluted by the inverter term, which does not depend on temperature. At the 2 points at rated speed where the voltage limit is active in at least one of the two states, the increase reaches 58%, because the larger resistive voltage drop of the hot machine tightens the voltage limit and forces a lower flux.

Table 3 compares the closed form (10) with differences of P^*. At these points no limit is active, the exact expression reproduces the numerical derivative to within 0.01%, and the simplified form differs by at most 0.4%.

Table 3. Marginal loss per unit torque: numerical derivative of P^* and Equation (10).

ϑ, °C	M*/MN	ω/ωN	Ψ2s, Wb	λM, numerical, W/(N·m)	λM, Eq. (10), exact, W/(N·m)	λM, Eq. (10), approx., W/(N·m)
20	0.1	0.5	0.394	15.98	15.98	15.95
20	0.25	0.5	0.627	14.40	14.40	14.37
20	0.5	0.5	0.891	13.61	13.61	13.57
20	0.5	0.9	0.858	14.69	14.69	14.64
20	1	0.5	0.900	25.61	25.61	25.58
20	1	0.9	0.900	25.63	25.63	25.58
20	1.25	0.5	0.900	31.59	31.59	31.55
115	0.1	0.5	0.398	20.05	20.05	20.01
115	0.25	0.5	0.634	18.47	18.47	18.43
115	0.5	0.5	0.900	17.70	17.70	17.66
115	0.5	0.9	0.874	18.80	18.80	18.72
115	1	0.5	0.900	34.33	34.33	34.28
115	1	0.9	0.900	34.36	34.36	34.28
115	1.25	0.5	0.900	42.48	42.48	42.44

Table 3. Marginal loss per unit torque: numerical derivative of P^* and Equation (10).

ϑ, °C	M*/MN	ω/ωN	Ψ2s, Wb	λM, numerical, W/(N·m)	λM, Eq. (10), exact, W/(N·m)	λM, Eq. (10), approx., W/(N·m)
20	0.1	0.5	0.394	15.98	15.98	15.95
20	0.25	0.5	0.627	14.40	14.40	14.37
20	0.5	0.5	0.891	13.61	13.61	13.57
20	0.5	0.9	0.858	14.69	14.69	14.64
20	1	0.5	0.900	25.61	25.61	25.58
20	1	0.9	0.900	25.63	25.63	25.58
20	1.25	0.5	0.900	31.59	31.59	31.55
115	0.1	0.5	0.398	20.05	20.05	20.01
115	0.25	0.5	0.634	18.47	18.47	18.43
115	0.5	0.5	0.900	17.70	17.70	17.66
115	0.5	0.9	0.874	18.80	18.80	18.72
115	1	0.5	0.900	34.33	34.33	34.28
115	1	0.9	0.900	34.36	34.36	34.28
115	1.25	0.5	0.900	42.48	42.48	42.44

3.2. Steady-State Torque Accuracy

Figure 3a shows the steady error of the quadratic formulation for weights from 0.1 to 100 times q_M,nom. Over three decades the simulated errors follow the 1/q_M law of Proposition 1. Across the 36 operating points of Table 4, the simulated error is 0.82–0.98 times the estimate of Proposition 1 for the reduced problem and 0.888–0.978 times its fixed-flux form (Section 2.5). Two effects account for the difference. First, the two estimates differ at light load, where ${\lambda }_M^{\Psi }\left(M_Q\right)\approx \left(1-{\epsilon }_M\right){\lambda }_M$, and on the voltage-limit boundary at rated speed. Second, input blocking lowers the error: with all inputs free (N = N_c = 10) the ratio to the fixed-flux estimate was 0.997–1.001, and the blocking used here reduced the error by 2.2–10.9%. With the nominal weight the error ranges from 0.46% to 3.88%, a spread by a factor of 8.4 for a weight that does not change; it is largest at 0.1M_N in the hot state. Both trends, 1/M^* at light load and the growth of R_σ with temperature, are those of Proposition 1 and Equation (10).

Figure 3b shows the absolute-value formulation. Below the threshold the controller sheds torque, and the error follows ε_M ≈ 1 − μ_M/λ_M: at fixed flux the marginal loss is nearly proportional to the torque, so the optimizer lowers the torque until its marginal loss equals μ_M. At μ_M = 0.95λ_M the error was 5.1–5.6%, at μ_M = λ_M at most 0.19%, and from μ_M = 1.05λ_M upward it fell below 10⁻⁴% in all four cases.

For the receding-horizon controller, the closed-loop threshold ${\mu }_{\mathit{cl}}^{*}$ (Section 2.5) was 0.999–1.002 times the fixed-flux marginal loss ${\lambda }_M^{\Psi }\left(M^{*}\right)$ at the 36 operating points (Table 4). Relative to λ_M(M^*) it was 0.850–1.002; the 3 points with a ratio clearly below 1 all lie on the voltage-limit boundary at rated speed, where ${\lambda }_M^{\Psi }\left(M^{*}\right)$ is 0.849–0.902 times λ_M(M^*). The 1 ms horizon, the linearization, the input blocking and the terminal flux term thus leave the local threshold unchanged to within 0.2%, and up to this tolerance λ_M(M^*) is an upper bound for it on this drive. At 18 of the points, λ_M decreases somewhere below M^*, and the global bound Λ(M^*) of Proposition 2 is up to 1.266 times λ_M(M^*); the closed loop is nevertheless exact from the local threshold on, because the controller only sees the local trade-off. The constant weight of 68.0 W/(N·m) exceeds ${\mu }_{\mathit{cl}}^{*}$ at all points by a factor of at least 1.47. A stricter requirement, zero predicted error at every step of the horizon, would call for weights of 0.88–3.92 times λ_M, and the constant weight is below this value at one point. It is not needed: the residual predicted errors appear in the later blocks, mostly in the held last block, where the planned drift of the flux toward its target cannot be compensated by a constant input, and they never reach the plant, because only the first input is applied and the predicted error of the first step was at most 10⁻¹⁵ of the reference.

Table 4. Steady relative torque error and loss at 36 operating points: quadratic term with q_M,nom (Q) and absolute-value formulations with μ_M = 68.0 W/(N·m) (L1, L1+Q) and κ = 1.35 (L1-S).

ϑ, °C	M*/MN	ω/ωN	λM, W/(N·m)	${\mathit{\lambda }}_{\mathit{M}}^{\mathit{\Psi }}/{\mathit{\lambda }}_{\mathit{M}}$	Λ/λM	${\mathit{\mu }}_{\mathit{cl}}^{*}/{\mathit{\lambda }}_{\mathit{M}}$	Q: Prop. 1, %	Q: fixed flux, %	Q: sim., %	ΔP, W
20	0.1	0.5	16.0	1.000	1.078	1.002	3.13	3.03	2.96	2.3
20	0.25	0.5	14.4	1.000	1.196	1.002	1.12	1.11	1.06	1.9
20	0.5	0.5	13.6	1.000	1.266	1.001	0.53	0.53	0.49	1.6
20	0.75	0.5	19.6	1.000	1.000	1.002	0.51	0.51	0.48	3.6
20	1	0.5	25.6	1.000	1.000	1.002	0.50	0.50	0.48	6.3
20	1.25	0.5	31.6	1.000	1.000	1.002	0.49	0.49	0.48	9.7
20	0.1	0.9	17.0	1.000	1.060	1.001	3.34	3.23	3.10	2.6
20	0.25	0.9	15.5	1.000	1.168	1.000	1.21	1.19	1.10	2.1
20	0.5	0.9	14.7	1.000	1.231	1.000	0.57	0.57	0.51	1.8
20	0.75	0.9	19.6	1.000	1.000	1.001	0.51	0.51	0.47	3.5
20	1	0.9	25.6	1.000	1.000	1.001	0.50	0.50	0.47	6.1
20	1.25	0.9	31.6	1.000	1.000	1.002	0.49	0.49	0.47	9.5
20	0.1	1	17.4	1.000	1.059	1.001	3.40	3.29	3.14	2.7
20	0.25	1	15.8	1.000	1.164	1.000	1.24	1.22	1.12	2.2
20	0.5	1	15.0	1.000	1.225	1.000	0.59	0.58	0.52	1.9
20	0.75	1	19.6	1.000	1.000	1.000	0.51	0.51	0.46	3.4
20	1	1	25.6	1.000	1.000	1.001	0.50	0.50	0.46	6.1
20	1.2	1	35.4	0.898	1.000	0.899	0.57	0.51	0.49	9.5
115	0.1	0.5	20.1	1.000	1.062	1.001	3.93	3.77	3.69	3.6
115	0.25	0.5	18.5	1.000	1.153	1.001	1.44	1.42	1.36	3.1
115	0.5	0.5	17.7	1.000	1.204	1.001	0.69	0.69	0.65	2.8
115	0.75	0.5	26.1	1.000	1.000	1.001	0.67	0.67	0.65	6.4
115	1	0.5	34.3	1.000	1.000	1.001	0.67	0.67	0.65	11.4
115	1.25	0.5	42.5	1.000	1.000	1.002	0.66	0.66	0.64	17.7
115	0.1	0.9	21.2	1.000	1.053	1.000	4.14	3.98	3.83	3.9
115	0.25	0.9	19.6	1.000	1.138	0.999	1.53	1.51	1.40	3.4
115	0.5	0.9	18.8	1.000	1.186	0.999	0.73	0.73	0.66	3.1
115	0.75	0.9	26.1	1.000	1.000	1.000	0.68	0.67	0.63	6.2
115	1	0.9	34.4	1.000	1.000	1.001	0.67	0.67	0.63	11.1
115	1.25	0.9	42.5	1.000	1.000	1.001	0.66	0.66	0.63	17.3
115	0.1	1	21.5	1.000	1.048	1.000	4.21	4.04	3.88	4.0
115	0.25	1	19.9	1.000	1.130	0.999	1.56	1.53	1.42	3.5
115	0.5	1	19.2	1.000	1.177	0.999	0.75	0.74	0.67	3.1
115	0.75	1	26.1	1.000	1.000	1.000	0.68	0.68	0.63	6.2
115	1	1	40.4	0.902	1.000	0.903	0.78	0.71	0.67	12.5
115	1.2	1	54.4	0.849	1.000	0.850	0.87	0.75	0.71	20.3

${\lambda }_M^{\Psi }$: marginal loss at the flux Ψ_2s(M^*). Λ/λ_M: ratio of the bound of Proposition 2 to the local value. ${\mu }_{\mathit{cl}}^{*}$: closed-loop threshold. Q, Prop. 1: estimate of Proposition 1 with λ_M at the delivered torque; Q, fixed flux: the same with ${\lambda }_M^{\Psi }$. The errors of L1, L1-S and L1+Q were below 10⁻⁴% at every point and are not listed. ΔP: loss of L1 minus loss of Q.

Table 4. Steady relative torque error and loss at 36 operating points: quadratic term with q_M,nom (Q) and absolute-value formulations with μ_M = 68.0 W/(N·m) (L1, L1+Q) and κ = 1.35 (L1-S).

ϑ, °C	M*/MN	ω/ωN	λM, W/(N·m)	${\mathit{\lambda }}_{\mathit{M}}^{\mathit{\Psi }}/{\mathit{\lambda }}_{\mathit{M}}$	Λ/λM	${\mathit{\mu }}_{\mathit{cl}}^{*}/{\mathit{\lambda }}_{\mathit{M}}$	Q: Prop. 1, %	Q: fixed flux, %	Q: sim., %	ΔP, W
20	0.1	0.5	16.0	1.000	1.078	1.002	3.13	3.03	2.96	2.3
20	0.25	0.5	14.4	1.000	1.196	1.002	1.12	1.11	1.06	1.9
20	0.5	0.5	13.6	1.000	1.266	1.001	0.53	0.53	0.49	1.6
20	0.75	0.5	19.6	1.000	1.000	1.002	0.51	0.51	0.48	3.6
20	1	0.5	25.6	1.000	1.000	1.002	0.50	0.50	0.48	6.3
20	1.25	0.5	31.6	1.000	1.000	1.002	0.49	0.49	0.48	9.7
20	0.1	0.9	17.0	1.000	1.060	1.001	3.34	3.23	3.10	2.6
20	0.25	0.9	15.5	1.000	1.168	1.000	1.21	1.19	1.10	2.1
20	0.5	0.9	14.7	1.000	1.231	1.000	0.57	0.57	0.51	1.8
20	0.75	0.9	19.6	1.000	1.000	1.001	0.51	0.51	0.47	3.5
20	1	0.9	25.6	1.000	1.000	1.001	0.50	0.50	0.47	6.1
20	1.25	0.9	31.6	1.000	1.000	1.002	0.49	0.49	0.47	9.5
20	0.1	1	17.4	1.000	1.059	1.001	3.40	3.29	3.14	2.7
20	0.25	1	15.8	1.000	1.164	1.000	1.24	1.22	1.12	2.2
20	0.5	1	15.0	1.000	1.225	1.000	0.59	0.58	0.52	1.9
20	0.75	1	19.6	1.000	1.000	1.000	0.51	0.51	0.46	3.4
20	1	1	25.6	1.000	1.000	1.001	0.50	0.50	0.46	6.1
20	1.2	1	35.4	0.898	1.000	0.899	0.57	0.51	0.49	9.5
115	0.1	0.5	20.1	1.000	1.062	1.001	3.93	3.77	3.69	3.6
115	0.25	0.5	18.5	1.000	1.153	1.001	1.44	1.42	1.36	3.1
115	0.5	0.5	17.7	1.000	1.204	1.001	0.69	0.69	0.65	2.8
115	0.75	0.5	26.1	1.000	1.000	1.001	0.67	0.67	0.65	6.4
115	1	0.5	34.3	1.000	1.000	1.001	0.67	0.67	0.65	11.4
115	1.25	0.5	42.5	1.000	1.000	1.002	0.66	0.66	0.64	17.7
115	0.1	0.9	21.2	1.000	1.053	1.000	4.14	3.98	3.83	3.9
115	0.25	0.9	19.6	1.000	1.138	0.999	1.53	1.51	1.40	3.4
115	0.5	0.9	18.8	1.000	1.186	0.999	0.73	0.73	0.66	3.1
115	0.75	0.9	26.1	1.000	1.000	1.000	0.68	0.67	0.63	6.2
115	1	0.9	34.4	1.000	1.000	1.001	0.67	0.67	0.63	11.1
115	1.25	0.9	42.5	1.000	1.000	1.001	0.66	0.66	0.63	17.3
115	0.1	1	21.5	1.000	1.048	1.000	4.21	4.04	3.88	4.0
115	0.25	1	19.9	1.000	1.130	0.999	1.56	1.53	1.42	3.5
115	0.5	1	19.2	1.000	1.177	0.999	0.75	0.74	0.67	3.1
115	0.75	1	26.1	1.000	1.000	1.000	0.68	0.68	0.63	6.2
115	1	1	40.4	0.902	1.000	0.903	0.78	0.71	0.67	12.5
115	1.2	1	54.4	0.849	1.000	0.850	0.87	0.75	0.71	20.3

${\lambda }_M^{\Psi }$: marginal loss at the flux Ψ_2s(M^*). Λ/λ_M: ratio of the bound of Proposition 2 to the local value. ${\mu }_{\mathit{cl}}^{*}$: closed-loop threshold. Q, Prop. 1: estimate of Proposition 1 with λ_M at the delivered torque; Q, fixed flux: the same with ${\lambda }_M^{\Psi }$. The errors of L1, L1-S and L1+Q were below 10⁻⁴% at every point and are not listed. ΔP: loss of L1 minus loss of Q.

The last column of Table 4 shows the price of the torque error. At every point Q dissipates less than L1, by 1.6–4.0 W up to 0.5M_N and by at most 20.3 W (at 1.2M_N, ω_N, 115 °C), because it delivers less torque. The difference corresponds to the marginal loss of the missing torque, λ_M(M^* − M); their ratio is 0.89–1.05. The lower loss of Q is thus the cost of under-delivery, not a better operating point, and the comparison at equal delivered work is made in Section 3.4.

3.3. Transients, Current Limit and Flux Reservation

Figure 4 shows the torque sequence in the hot state without flux reservation. Before the first step the machine runs at 0.2M_N with the loss-optimal flux of 0.56 Wb. The step to 0.9M_N is therefore limited by the current bound until the flux has been rebuilt. The step to 1.4M_N asks for more than can be delivered within the current limit, M_max = 1.300M_N at this speed. Table 5 lists the metrics.

Figure 4. Torque, rotor flux and stator current amplitude during a torque sequence at 0.75ω_N, 115 °C, without flux reservation: (a) the step from 0.2M_N to 0.9M_N; (b) the whole sequence, with a reference of 1.4M_N beyond the current limit between 350 and 550 ms. L1-S is listed in Table 5.

Figure 4. Torque, rotor flux and stator current amplitude during a torque sequence at 0.75ω_N, 115 °C, without flux reservation: (a) the step from 0.2M_N to 0.9M_N; (b) the whole sequence, with a reference of 1.4M_N beyond the current limit between 350 and 550 ms. L1-S is listed in Table 5.

After the step all formulations first run at the current limit, and the torque rise is set by how fast the flux can be rebuilt. Q with q_M,nom divides the current between i_1d and i_1q and reaches 90% of the step after 36.5 ms. With 10q_M,nom more current goes to i_1q; the 90% point is reached after 15.5 ms, at the price of a slower flux build-up. L1 behaves differently. Its penalty grows only linearly with the error, so the terminal flux term dominates during the rebuild: the current is spent on the flux first. The torque even falls below its initial value, practically to zero, stays below 0.2M_N for 6.7 ms, and reaches half of the step only after 30.2 ms, against 6.0 ms for Q. Its IAE_rel over the step is 7.56%, against 4.07% for Q; L1-S behaves like L1 (7.71%), and L1+Q follows the quadratic transient (3.47%). In steady state Q settles at 0.894M_N and 0.497M_N, whereas all formulations with the absolute-value term deliver 0.900M_N and 0.500M_N. H also reaches the reference in steady state, but its QP was infeasible in 33.5% of the steps of the sequence, which were covered by the quadratic fallback.

Under the 1.4M_N reference, L1 keeps the flux at 0.897 Wb or above and holds 1.297M_N, practically the largest torque admissible under the current polygon, M_max = 1.300M_N. Q with q_M,nom and L1+Q hold 1.267M_N and 1.264M_N with a flux of 0.869 Wb and 0.866 Wb at its lowest. With 10q_M,nom, the controller keeps shifting current from i_1d to i_1q to shrink a large predicted error over a 1 ms horizon; the flux decays to 0.53 Wb, and the torque, 1.05M_N on average over the interval (Table 5), falls to 0.85M_N by its end. This flux depletion results from a large quadratic weight combined with a horizon much shorter than the rotor time constant. It is not a property that the absolute-value term guarantees to avoid, and the behavior of L1 at the current limit is an observation for this scenario.

Figure 5 and the lower part of Table 5 identify the cause of the L1 delay. With flux reservation, the 100 ms preview sees the step in advance and the flux is raised from 0.56 Wb to 0.90 Wb before it, which costs 91 W of additional steady-state loss while the load is still at 0.2M_N. When the target is raised, all three formulations briefly lose torque, down to 0.05M_N for Q and 0.00M_N for L1; averaged over the 50 ms before the step, the torque is 0.197–0.199M_N. The 90% rise time then drops to 1.1 ms for Q, 1.2 ms for L1 and 1.1 ms for L1+Q, and the IAE_rel of the step to 0.78%, 0.15% and 0.14%. With reservation, L1 and L1+Q have the smallest integral errors, because Q keeps its steady bias. Lengthening the horizon without reservation changes little: for L1 the time to half of the step is 30.2, 29.3 and 27.5 ms for N = 10, 20 and 40, since even 4 ms remains short compared with the rotor time constant of 131 ms. The delay of L1 in Figure 4 is therefore a consequence of rebuilding the flux under a short horizon, not of the absolute-value term as such.

3.4. Duty Cycles

Table 6 and Figure 6 summarize the duty cycles. In the hot state, Q with q_M,nom leaves a mean plateau error of 0.97% on the pump cycle and 1.26% on the conveyor cycle, and 0.1q_M,nom raises it to 9.2% and 11.4%. Without flux reservation, all formulations with the absolute-value term, including ${\mu }_M={\bar{\lambda }}_M$ and the scheduled weight, bring the plateau error below 10⁻⁴% in both thermal states.

For the same reservation setting, the loss energy per unit of work hardly depends on the formulation. Without reservation, the absolute-value formulations change E_loss/W relative to Q with q_M,nom by +0.03% to +0.57% of its value; across all hot-state runs without reservation, E_loss/W spans 17.89–18.02% on the pump cycle and 19.78–19.89% on the conveyor cycle. On the pump cycle the reference changes along speed ramps, and the integral error is almost entirely steady bias: IAE_rel is 0.79% for Q and below 0.001% with the absolute-value term. The conveyor cycle contains four steps from a lighter to a heavier load. Here L1 without reservation reaches an IAE_rel of 1.44–1.75%, close to the 1.63% of Q, because the zero plateau error is offset by the delayed rise after each step (Section 3.3). L1+Q lowers it to 0.75%; 10q_M,nom reaches 0.36% but leaves a plateau error of 0.12%. With flux reservation, IAE_rel falls to 0.95% for Q, 0.17% for L1 and 0.08% for L1+Q, while E_loss/W rises to 19.94%, 19.96% and 19.95% because of the higher standby flux; reservation itself thus adds 0.10–0.16 percentage points on the conveyor cycle. The plateau errors of L1 and L1+Q with reservation, 0.17% and 0.05%, arise in the last 100 ms before a load increase, while the flux is being raised; over the rest of the plateaus the mean error of L1 is below 10⁻⁴%. In these windows the error of Q departs from its steady bias by up to 4.2 percentage points and that of L1 reaches 3.4%, so the deviation is caused by the flux ramp rather than by the form of the torque term. H needed its fallback in 0.0% and 1.5% of the steps of the two cycles; on the conveyor cycle its IAE_rel was 0.74%.

3.5. Uncertainty of the Loss Model

The torque equation (4) does not contain the loss coefficients k_h, k_e, a_I and b_I. Errors in these coefficients therefore change the flux chosen by the controller and the loss it incurs, but not the torque it delivers, and exactness is decided by λ_M of the controller’s own loss model, against which μ_M was set. The relevant question for a constant weight is therefore whether it can be fixed in advance when the loss model of the controller is revised later, for example after a new identification of the coefficients or by online adaptation. In 1500 Monte Carlo samples, k_h, k_e, a_I and b_I were scaled by four independent factors drawn uniformly from [0.7, 1.3], and the winding temperature was drawn uniformly from 20–115 °C. For each sample, λ_M was computed at 18 operating points: the 12 points at 0.5ω_N and 0.9ω_N of Table 4 and the six torque levels at rated speed, where the voltage limit becomes active at the highest torques; all 27,000 point evaluations were feasible for the perturbed models. Figure 7a shows the result. λ_M ranged from 13.0 to 55.2 W/(N·m); at the 12 points below rated speed it did not exceed 43.2 W/(N·m), and the larger values all occurred at rated speed. The constant weight μ_M = 68.0 W/(N·m) stayed above the per-sample maximum in all samples; it exceeds the largest sampled value by 23% and the 99th percentile of the per-sample maximum, 54.5 W/(N·m), by 25%. With the scheduled weight the condition holds by construction, since λ_M is taken from the same loss model that the optimizer uses. The temperature range matters as well: λ_M grows with the winding temperature, and the constant weight was set for the hot state of [14]. At 155 °C, the limit of thermal class F, the largest λ_M over the envelope of the nominal model is 65.8 W/(N·m), still below μ_M.

3.6. Field-Orientation Error

The results so far assume ideal field orientation. To show how the penalty-induced error compares with the error caused by wrong rotor parameters, the controller was given values of R₂ and L_m perturbed by uniform factors in [1 − ε,1 + ε] and [1 − ε/2,1 + ε/2], following the perturbation scheme of [14], while the plant kept the nominal values. An error of R₂ by ε also represents an error of the rotor temperature estimate, about 7 K for ε = 2% and 34 K for 10% at 115 °C. The controller estimates the flux with the current model $d{\widehat{\Psi }}_2/\mathit{dt}=\left({\widehat{R}}_2/{\widehat{L}}_2\right)\left({\widehat{L}}_m{i}_{1d}-{\widehat{\Psi }}_2\right)$ and orients its frame with the estimated slip; the plant is simulated with both rotor flux components in that frame, so the torque on the shaft is ${\displaystyle \frac{3}{2}}{z}_p{k}_2\left({\Psi }_{2d}{i}_{1q}-{\Psi }_{2q}{i}_{1d}\right)$. For each ε ∈ {2%,10%}, 100 samples were simulated at two operating points in the hot state (Table 7, Figure 7b). With correct parameters the errors of L1 are zero and those of Q equal the penalty bias, 1.36% and 0.63%. With mismatch, the true error of both formulations is dominated by the orientation error, which reaches several percent at ε = 10%. For the same perturbed parameters, the true error of Q exceeds that of L1 by a median of 1.37 and 0.65 percentage points at ε = 10% for the two operating points (5th–95th percentiles 1.01–1.69 and 0.32–0.97), which is close to the penalty bias of Q at these points. The two contributions thus add up, and the absolute-value term removes one of them. The error of the torque estimated by the controller is also nonzero with mismatch (median absolute value 0.19% for L1 over both values of ε and both points), because the prediction model no longer matches the plant and the controller has no integral action; this offset is of the kind addressed by offset-free MPC [28].

3.7. Measurement Noise and Dead Time

The results above assume exact measurements and an ideal inverter. Two disturbances were therefore added at three operating points in the hot state (Table 8). Gaussian noise with a standard deviation σ_i of 0.1 A or 0.4 A, about 0.45% and 1.8% of I_1N, was added to both measured current components at every sample, and the controller estimated the flux with the current model driven by the noisy i_1d. Dead time was represented by its average effect, a voltage vector of amplitude (4/π)t_df_swU_d = 13.8 V opposite to the current vector for t_d = 2 µs, which the controller did not compensate. Each run lasted 0.6 s, the statistics were taken over its last 0.3 s, and each case with noise was repeated with 10 noise realizations.

With measurement noise the absolute-value term is no longer exact. The mean error of L1 was 0.01–0.15% for σ_i = 0.1 A and 0.48–0.69% for 0.4 A, largest at light load, whereas the error of Q changed by at most 0.01 percentage points for 0.1 A and by −0.03 to +0.23 percentage points for 0.4 A; the 95% confidence intervals are at most ±0.13 percentage points wide. The error of L1 stayed below that of Q at every point, but a weight above λ_M does not ensure zero error once the measured torque fluctuates around the reference. A diagnostic run with the term a_II_1m linearized at the noise-free current instead of the measured one left the error of L1 at 0.1M_N, 0.5ω_N and σ_i = 0.4 A unchanged (0.69% in both cases), so the upward bias of the amplitude of a noisy current is not the cause. The origin of the bias was not identified. The kink of |e| is a plausible contributor for L1, since below the reference the loss term opposes the penalty and above it adds to it, so that apparent excesses of torque are corrected more strongly than apparent shortfalls; however, the error of Q also changes with noise, by −0.03 to +0.23 percentage points at 0.4 A, so part of the effect is common to both formulations. Uncompensated dead time produced a steady offset of 0.93–2.55% for L1 and L1+Q and raised the error of Q by 1.14–3.70 percentage points. Like the orientation error of Section 3.6, this offset comes from a mismatch between the prediction model and the plant, and the absolute-value term does not remove it.

No self-sustained chattering was observed: with dead time but without noise, all formulations settled to a constant input. With noise, L1 passed more of it to the input than Q. The RMS change of the voltage vector between samples was 1.6–2.2 times that of Q for σ_i = 0.1 A and 1.2–1.3 times for 0.4 A, and the standard deviation of the shaft torque was 8–12% higher at 0.4 A. L1+Q behaved similarly, with slightly smaller noise-induced errors at light load.

3.8. Closed-Loop Behavior

The propositions concern the equilibrium, not convergence to it, and no stability proof is available for the nonsmooth absolute-value term. Three numerical checks were therefore made (Table 9). First, the one-step closed-loop map, with the plant state and the previous input as its state, was linearized at the equilibrium of each of the 36 operating points by finite differences, in both directions along each coordinate. The state of the map contains the rotor flux, so the flux dynamics are part of the test. The open-loop rotor pole, $e^{-T_s{R}_2/L_2}$ = 0.9992 at 115 °C, does not appear in the closed-loop spectrum, because the controller acts on the flux through i_1d. The winding temperatures, which change over minutes, are treated as parameters. The largest spectral radius was 0.9084 for Q, 0.9254 for L1 and 0.9254 for L1+Q; the slowest mode of these linearizations has a time constant of at most 1.29 ms. L1 and L1+Q give the same spectral radius to four decimals at 35 of the 36 points (largest difference 1.5 × 10⁻⁴): at the equilibrium the slack constraints hold the predicted error at zero, and the quadratic term, whose gradient vanishes there, hardly changes the local control law. Because the control law of a QP is piecewise affine and, with the absolute-value term, the equilibrium lies where several affine pieces meet, a linearization along the coordinate directions need not represent every piece. Second, therefore, the closed loop was started from 20 small perturbations of each equilibrium in random directions and run for 200 steps. All trajectories returned to the equilibrium; the largest decay rate per step over steps 60–160 was 0.9084 for Q, 0.9255 for L1 and 0.9255 for L1+Q, in agreement with the spectral radii, although the scaled deviation first grew by a factor of up to 397; this transient occurs within the first steps after a perturbation of the flux, which the controller corrects through i_1d and the input. Third, at the 18 operating points in the hot state, six simulations per formulation were started from large perturbations, with the currents scaled by factors between 0.5 and 1.5 and the flux by factors between 0.6 and 1.4. All 108 runs of each formulation returned to the equilibrium, with a final deviation below 10⁻¹² in scaled units. These results show convergence for the tested points and perturbations. They do not establish a region of attraction or stability in general.

3.9. Computational Cost

The absolute-value term enlarges the condensed QP from 6 to 9 decision variables and from 116 to 139 inequality rows (Table 10), i.e., by N_c variables and 2N + N_c rows. The number of active-set iterations grows more than the problem size. During the sequence of Figure 4, the Goldfarb–Idnani method needed a median of 1 and a maximum of 16 iterations for Q, against 15 and 20 for L1 and 13 and 22 for L1+Q: near e = 0 the slack constraints of (12) are active in every block and have to be identified at each step. Each iteration evaluates all constraint rows and updates a factorization of the size of the decision vector, about 1332 multiply–adds for L1 against 732 for Q. With the median iteration counts, the QP itself thus costs roughly 20 × 10³ instead of about 0.7 × 10³ operations, which is no longer small compared with the condensing step of about 1.6 × 10³ operations estimated in [14]. On a desktop processor the median solution time rose only from 38 µs to 66 µs, because the interface overhead dominates there, and the maximum times, 7.1–13.4 ms for the four formulations, are attributed to operating-system delays rather than to the algorithm. A repeated run on the same type of processor (Section 2.6) gave identical iteration counts but median times of 35 µs for Q and 58 µs for L1 and maxima of 0.8–16.6 ms. These measurements therefore cannot decide whether either formulation runs within the 100 µs period on a drive controller. The operation count suggests that the L1 formulation will need a warm-started active set, which is known in steady state, or an explicit solution of the QP to do so; this has not been verified.

4. Discussion

Two effects are mixed in a quadratic torque term, and the results separate them. The first is a steady bias. Under an economic objective it is the optimal answer to the cost the controller was given: every newton-meter costs λ_M watts, and a finite quadratic weight lets the optimizer trade the last fraction of torque for loss. Proposition 1 puts a number on that trade and accounts for the 1/q_M scaling seen in [14]. It also means that a weight tuned at one point cannot hold a given accuracy across the envelope, because λ_M/M^* varies by a factor of 9 over the operating points of Table 4. The absolute-value term removes the bias without an integrator and without an equality constraint that could become infeasible. Its weight has a lower bound in the same unit as the economic cost. On this drive the receding-horizon QP becomes exact at the fixed-flux marginal loss ${\lambda }_M^{\Psi }\left(M^{*}\right)$, which equals λ_M(M^*) except on the voltage-limit boundary, where it is lower (Table 4); the global bound of Proposition 2 adds a margin where λ_M is not monotonic.

The second effect concerns transients. A linear penalty prices a large and a small torque error equally per newton-meter. During a current-limited flux rebuild, L1 therefore gives way to the terminal flux term: without flux reservation it first lowers the torque and reaches half of the step only after 30 ms (Figure 4a). The comparison with the reservation mechanism of [14] and with longer horizons (Figure 5) shows that the lag stems from the flux dynamics, which a 1–4 ms horizon cannot see, and not from the norm of the penalty. With reservation, L1 and L1+Q have the smallest integral errors of all formulations. Without it, L1+Q keeps the quadratic transient and the exact steady state. The quadratic weight of L1+Q is still a tuning parameter, but it now affects only the transient response, not the steady accuracy. Its useful range is bounded on both sides (Table 5). With 0.1q_M,nom, L1+Q behaves like L1, with a 90% rise time of 32.5 ms, and with 10q_M,nom it reproduces the flux depletion of Q at the current limit: the flux falls to 0.53 Wb and the torque to 0.85M_N. The scheduled weight L1-S follows λ_M exactly, but it keeps the transient of L1 (Table 5); L1+Q uses a constant μ_M and obtains its transient behavior from the quadratic part. Both are exact in steady state.

The approach differs from existing weighting methods in what it fixes. Empirical and optimized weights [15,16,17] and ranking methods [18,19] balance tracking terms of the same kind and do not price torque against loss. The sequential evaluation of [20] imposes a strict priority on torque, which behaves like a very large weight and gives no graceful degradation when the reference cannot be met. Offset-free MPC with a disturbance model estimated from the prediction error [28] removes offsets caused by model mismatch, but not the penalty-induced bias: with a correct model the estimated disturbance is zero, and the bias is the optimum of the controller’s own problem. A target calculation that imposes the reference as an equality, or an integrator on the torque error, would remove the bias as well, at the cost of an additional state and of the infeasibility that the penalty was meant to avoid. Section 3.6 shows that both kinds of error occur together and add up, and the absolute-value term does not remove the offset caused by uncompensated dead time either (Section 3.7). With measurement noise it is no longer exact, although its error remained below that of Q, and it passes more of the noise to the input. In an implementation it would therefore be combined with dead-time compensation and filtering of the measured currents, or with a disturbance model.

The practical scope of the result is limited. The penalty-induced error removed here is 0.46–3.88% with the nominal weight. Errors from rotor-parameter mismatch reach several percent at ε = 10% (Table 7) and, according to [14], about 8% for a controller without thermal adaptation, i.e., with the resistances kept at their cold values while the machine heats up. The absolute-value term does not reduce these errors; it removes a bias that the controller would otherwise add on top of them, and it removes it without retuning when temperature or operating point change. The weight also has a meaning beyond tuning: λ_M is the marginal electrical cost of torque, and it could serve as a price signal for an energy-management layer that schedules process throughput. This use is a possibility only and has not been examined here.

The study has limitations of four kinds. First, the plant: all results were obtained in simulation with an average-value inverter model and with the speed imposed by the load; dead time and measurement noise enter only in Section 3.7, in a simplified form. Second, model knowledge: except in Section 3.6 and Section 3.7, field orientation was ideal, and the thermal observer of [14] was replaced by the true temperatures. The R₂ perturbation of Section 3.6 also stands for an error of the rotor temperature estimate, but errors of R₁, Ψ_2N and the dc-link voltage were not studied. The parameter set is representative of a 7.5 kW machine rather than taken from a specific motor, and Ψ_2N, α₂ and f_sw were set in this study. Third, implementation: the QP times were measured on a desktop processor and indicate only the relative cost of the reformulation. Whether the controller, with either torque term, can be executed within the 100 µs sampling period on a drive controller has not been established and requires on-target measurement. The flux reservation assumes that the reference is known 100 ms in advance; where no preview is available, the behavior without reservation in Section 3.3 applies. Fourth, theory: no stability guarantee is claimed for the formulations with the absolute-value term. The stability analysis of [14] was derived for a smooth stage cost and was not extended here, and the evidence of Section 3.8 is numerical and limited to the tested operating points and perturbations. Proposition 2 concerns the reduced steady-state problem; the receding-horizon controller was found to follow the local condition ${\mu }_M>{\lambda }_M^{\Psi }\left(M^{*}\right)$, which was verified only at the simulated operating points. Further work will address validation on a loaded dynamometer with on-target timing, combination with dead-time compensation and an offset-free disturbance model, and longer prediction horizons through move blocking or explicit MPC.

5. Conclusions

The quadratic torque term of a thermally adaptive economic MPC for induction motor drives was replaced by an absolute-value term whose weight is expressed in watts per newton-meter and is bounded from below by the marginal loss per unit torque λ_M of the loss-optimal steady state. An analytic expression and a torque–speed map of λ_M were derived from the machine and inverter loss model. Simulation of a 7.5 kW drive led to the following conclusions.

• The quadratic term causes a steady torque error ε_M = λ_M/(2q_MM^*). With the nominal weight it ranged from 0.46% to 3.88% over the 36 operating points. The simulated errors were 0.82–0.98 times the estimate of Proposition 1 and 0.888–0.978 times its form for a flux held at the steady-state target; without input blocking the latter ratio was 0.997–1.001.
• With the constant weight μ_M = 68.0 W/(N·m), the scheduled weight or L1+Q and exact measurements, the steady error fell below 10⁻⁴% at all 36 operating points. The closed-loop threshold of the receding-horizon controller, located by bisection, was 0.999–1.002 times the fixed-flux marginal loss ${\lambda }_M^{\Psi }\left(M^{*}\right)$, which equals λ_M(M^*) except on the voltage-limit boundary, where it is up to 15% lower; the global bound Λ(M^*) of Proposition 2 is up to 27% higher than λ_M(M^*) at these points.
• On pump and conveyor cycles without flux reservation, the absolute-value term reduced the mean plateau error from 0.97% and 1.26% to below 10⁻⁴%, while the loss energy per unit of work changed by +0.03% to +0.57% of its value.
• Without flux reservation, L1 lowers and then delays the torque during a current-limited flux rebuild; L1+Q avoids this. With the reservation of [14], all formulations reach 90% of the step within 1.2 ms, L1 and L1+Q with the smallest integral errors.
• Under ±30% error in the core and inverter loss coefficients, λ_M stayed within 13.0–55.2 W/(N·m); the constant weight ${\mu }_M=1.2{\bar{\lambda }}_M$ = 68.0 W/(N·m) remained above it in all 1500 samples.
• With mismatched rotor parameters, the orientation error dominates the shaft torque error; the penalty bias adds to it, and the absolute-value term removes that part.
• With current-measurement noise of 1.8% of I_1N, the error of L1 rose to at most 0.69%, still below that of Q, and the input activity increased; uncompensated dead time added an offset of up to 2.55% that the absolute-value term does not remove.
• Linearization at all equilibria gave spectral radii of at most 0.9254, and all simulations from small and large perturbations converged; this is numerical evidence, not a stability proof.
• The reformulated QP has 3 more variables and 23 more inequality rows and needed a median of 15 instead of 1 active-set iterations; execution within the 100 µs period on a drive controller has not been verified.

For practical use, the combined term L1+Q is recommended, preferably together with flux reservation and with a quadratic weight of moderate size, since a large one reintroduces the flux depletion of Q at the current limit. Its absolute-value weight follows from the loss model and need not be retuned when temperature or operating point change, while the quadratic part, which now affects only transients, can be chosen for the desired dynamic response.