Switched electromechanical dynamics for transient phase control of brushed DC servomotor

In typical robot operations, the direct control input to actuators is a varying voltage profile, while the input from the external environment is provided as position/velocity or force/torque feedback. This combination of electromagnetic and mechanical variables, which each exhibits their own characteristic dynamics, requires special consideration in modeling and control. However, typical modeling approaches for electromechanical dynamics consider steady-state and/or periodic phases or use system parameters that are estimated for such conditions.^1,2 As robotic tasks often exceed the scope of steady-state or periodic behavior,^3,4 generally-applicable models are needed specifically for actuators intended to generate transient or aperiodic motion.

Modeling electromechanical systems that cross multiple domains as part of their operation requires the consideration of all their relevant components. Yet, existing electromechanical models of robot actuators typically consider only the DC motor^5–10 because its dynamics is well understood.¹¹ However, a simple DC motor has no inherent mechanism for directional, speed, or torque control. In practice, such control mechanisms are provided by power converter circuits that modulate the voltage applied across the motor terminals. Phase-controlled converters are used to control brushless DC motors, while chopper converters are commonly used to control brushed DC motors through pulse-width modulation (PWM), in which the mean voltage applied at the motor terminals is adjusted through rapid switching of the DC input voltage from the power supply.^12,13

The combination of a motor, gear train, power converter, and feedback control mechanism is referred to as a servomotor. Due to their simple design and dynamic performance capable of rapid acceleration and deceleration, DC servomotors with chopper converters are used as actuators in a broad range of applications.¹¹ The ubiquity of chopper converters and similar switching devices for control in electronics allows for ease of implementation but also introduces non-smooth switched dynamics at the electromagnetic level. Existing electromechanical models of the DC servomotor derived for use in control^1,14,15 typically do not consider all possible switching configurations of the chopper converter during operation. Within mechanical components, such as the gear train, there are the direction-dependent effects of back-drivability¹⁶ and friction,^17,18 which are additional non-smooth aspects.

Analytical approaches to electromechanical modeling typically assume that the switched electromagnetic dynamics are sufficiently fast relative to the mechanical dynamics, in terms of response to varying inputs, such that the electromagnetic and mechanical response can be fully or partially decoupled, as can be demonstrated by singular perturbation theory.¹⁹ The switched dynamics of the power converter is often modeled with an averaging approach (as opposed to data-driven methods, such as model-free reinforcement learning²⁰) that obtains averaged electrical state trajectories from mechanical state trajectories with harmonics of much lower frequencies than the switching frequency. For sufficiently high switching frequencies, such as those used for PWM, the characteristics of the averaged dynamics are similar to those of the full-order switched dynamics. For example, the control designs that stabilize an averaged system will also stabilize the original non-averaged system.^21,22

State-space averaging is often applied to obtain the transfer functions of switching power converters. By taking a weighted average of the state equations of each switching configuration according to their relative time durations, small-signal and steady-state models of the circuit can be obtained under the assumption that the state variables are perturbed around some steady-state operating point.¹² Variations of this method include corrected full-order state-space averaged modeling,²³ which is the most accurate method for both low and high switching frequencies, and parametric average-value modeling,^24,25 which is a computational method that extends state-space averaging to complex systems that are difficult to solve analytically.

In circuit averaging, the input and output variables of each circuit element are averaged across one switching cycle to derive its cycle-averaged input–output relations instead of averaging the state equations directly. The switching power converter is then modeled as an equivalent circuit composed of averaged circuit elements.²⁶ Other circuit-averaging methods include the injected-absorbed-current method, which is mostly suitable for low-frequency, small-signal phenomena.²⁷ 

Averaging approaches are sufficient in applications where only low-frequency behavior is considered, such as typical trajectory tracking control. Modeling averaged rather than full-order switched dynamics also reduces the computational effort because less dense discretization of the solution time interval is required for convergence. However, any high-frequency oscillatory effects in the state variables, such as ripple current, must be recovered from the averaged model with estimation techniques that result in some degree of error.^28,29 For instance, for the circuit-averaging approach, additional relations, such as the energy conservation principle, must be introduced to recover the ripple current from the cycle-averaged current profile in order to model power consumption.^30–32 The use of linearization in deriving averaged dynamics^12,33,34 also results in discrepancies against the full-order nonlinear dynamics. Thus, there remains a need for accurate multi-domain modeling of switched electromechanical systems integrated with averaging approaches to pursue more sophisticated control, especially for transient motion profiles, despite the well-established use of existing models.^1,35

Averaging switched electromechanical dynamics also allows for the modeling of power consumption in general complex tasks. Energy usage is one of the most important criteria for mobile^3,36 or legged robots^{5,8,14,37–43} (e.g., to conserve energy through real-time task management⁴⁴) as well as for manipulators^{4,6,7,10,15,45} to resolve mechanical redundancy. Despite the need for accurate mathematical models to predict and minimize power consumption, broadly applicable (including transient phases and aperiodic tasks), physically-based mathematical models of power consumption are lacking. While bond-graph techniques^37,46,47 provide tools for analyzing interconnections among the multi-domain elements of a system, the specific form of each model term is still required for the quantification of power consumption. Mathematical models of robot power consumption are usually required to be functions of kinematic (joint angles and velocities) and dynamic (joint torques, external forces) variables in the mechanical domain so that they can be used in prediction, design, control, and optimization. The current state of knowledge still relies on ad hoc, empirically-based models calibrated to specific instances, and incomplete proxies like mechanical energy or work,^38,40 square of actuator torques or forces,^3,45,48–50 or their combinations with other terms^36,37,51 are often used to model power consumption instead. While the mechanical power alone may provide some implications about power consumption under periodic or steady-state conditions,³⁸ and to some extent, during positive mechanical-work phases,^43,52 it is different from the true power consumption of a system.¹⁶ 

The aforementioned proxies of power consumption, which are overly simplified, generally cannot differentiate between the significance of positive, zero, and negative output mechanical work or the power consumption of different switching configurations in a controlled system. Actuator torque or its squared value, for instance, lacks the velocity-dependence expected for power consumption and does not imply the actual quantity. When a more complete model of electrical power consumption is considered in the literature,^10,15,39 there is often a condition imposed that omits the energetic effect of negative mechanical or electrical work,^5,16,43 which may be inconsistent with experimental results.^6,14 Data-driven approaches,³⁶ which include, for instance, the use of reinforcement learning for improving locomotion efficiency,⁴² simply sidestep the complexities of multi-domain interactions by collecting more experimental data. While some of the existing approximations may provide acceptable trends for well-defined tasks that are periodic or steady-state, generally applicable models are required for general complex tasks for broader robot applications. Most recent works^47,53–56 along this line are focused on predicting the power consumption of electric vehicles, for which the models are relatively successful in integrating the drive cycle of actual vehicles. However, they present similar issues as described above, i.e., lack of generality and accuracy due to over-simplification.

In this study, a multi-domain framework is established for switched electromechanical dynamics in servomotor systems for their analysis and control, in general, aperiodic tasks, including transient phases. While the approach is applicable to general electromechanical systems with switching-based control, the focus of this work is on permanent-magnet, armature-controlled brushed DC servomotors due to their ubiquity as actuators in robotics, in addition to their simple designs without field windings,⁵⁷ a wide range of torque and speed controllability with identifiable torque-speed characteristics, and compatibility with various control methods.⁵⁸ The complex system-of-systems challenge is addressed in a principled, step-by-step approach where the switched electromechanical dynamics is derived from the individual models of the internal DC motor, gear train, and H-bridge circuit. The coupled models comprehensively integrate all possible distinct switching configurations of the H-bridge and are cycle-averaged to handle non-smoothness and different temporal scales between electrical and mechanical variables. Experiments were conducted to estimate the system parameters and validate the model's prediction accuracy in two distinct tasks—dynamic braking and sinusoidal trajectory following. The model was also used to formulate the servomotor power consumption, which was implemented for optimal control and energy analysis and demonstrated for its merits against other common proxy models.

The main components of the DC servomotor (Fig. 1) are the armature DC motor, gear train, full-bridge inverter (referred to here as the H-bridge), microcontroller, and sensors.^11–13 The system model and the corresponding system parameters are obtained here for the DYNAMIXEL MX-28 servomotor (ROBOTIS Co. Ltd, South Korea),⁵⁹ which is used in many robotic systems,^60–63 but the proposed models apply to general permanent-magnet, armature-controlled brushed DC servomotors due to their common operating principle. To obtain the full state model, the state equations are derived separately for each component and then integrated in terms of electrical and mechanical state variables, where the H-bridge and the DC motor are modeled with equivalent circuits. The electrical power supply is modeled as an ideal voltage source with constant input voltage V_in, where the input current I_in from the power supply is positive in the direction of positive V_in. Due to the constant input voltage, the capacitor typically in parallel with V_in in the power supply does not affect the servomotor dynamics. The electrical current and power consumption drawn by the sensors and microcontroller¹⁶ are not considered here and are assumed to be constant or otherwise motion- and actuation-independent. It is further assumed that the entire system of interest is isolated from any external electromagnetic field effects.

The classic equivalent circuit model of the permanent-magnet, armature-controlled brushed DC motor consists of the inductor, resistor, and back electromotive force (back emf) elements in series^7,8,9,11 (Fig. 1). The voltage V_PWM,H(t) at time t applied by the H-bridge to the DC motor is

where L is the motor inductance, R_a is the armature resistance, $Vemf(t)=Kq˙rot(t)$ (from Faraday's law) is the back emf voltage with motor torque constant K and rotor shaft angle q_rot, I_a is the motor armature current, and the dot indicates time-derivative. The voltage drop of magnitude V_br due to brush-contact resistance associated with slip rings and commutators, which changes sign when the motor current is reversed,¹⁴ is also included in series with the other elements as an additional non-smooth term. Ferromagnetic coupling losses due to eddy current or hysteresis that are usually neglected in control^57,64 and a possible capacitor in parallel with the DC motor in certain designs to suppress electrical noise¹ are not considered here.

The DC motor is attached to a compound reduction gear train of N_G gears, where the ith gear shaft angle q_G,i for i = 1 and i = N_G are the rotor shaft angle q_rot and the servomotor output shaft angle q, respectively. Under a rigid-body assumption, the equivalent moment of inertia inclusive of all rotating mass (rotor inertia, gear train, and output shaft) as reflected on the servomotor output shaft angle is

where $Gi=q˙G,i/q˙$ is the signed gear ratio,¹¹ with $G1=q˙rot/q˙=G$ for brevity, and J_i is the axial moment of inertia of the gear shaft corresponding to q_G,i.

The equation of motion of the combined system of the DC motor and the gear train is

where $τrot(t)=KIa(t)$ (from Lorentz's law) is the rotor torque applied by the stator on the rotor of the motor (where K is assumed to be identical to that used for the back emf), $τf$ is the effective frictional torque (resultant applied by the servomotor housing on gear shaft axles and due to gear contacts) as reflected to the output shaft, and $τout$ is the output torque applied by the servomotor on the external environment through the output shaft. The non-smooth model of frictional torque is $|τf|≤b0$ for $q˙=0$ and $τf=b0sgn(q˙)+b1q˙$ for $q˙≠0$⁠, where b₀ and b₁ are the Coulomb and viscous friction coefficients, respectively, and are negative in sign. Both the load-dependence of $τf$ and the break-away friction are not considered here.

The H-bridge, which consists of four switches (Fig. 2), allows for bidirectional control through PWM. By only closing switches S1 and S4, S2 and S4, and S2 and S3, a voltage of approximately V_in, 0, −V_in is applied to the DC motor, respectively. These three switching configurations are the positive on-state, off-state, and negative on-state, respectively (Table I). By varying the switch timing, a desired variable voltage profile for V_PWM,H(t) within the interval [−V_in, V_in] is approximated with a binary or ternary voltage profile with the values V_in, 0, or −V_in that resembles a series of variable-width pulses occurring at a fixed frequency 1/T_PWM, where T_PWM is the fixed pulse period.⁶⁵ 

The direct control input to the DC servomotor is represented by the signed duty cycle function $D(t)∈[−1,1]$⁠, where the sign and magnitude indicate the direction and the ratio of the on-state duration to T_PWM of each pulse period, respectively¹¹ (Fig. 2). In digital PWM implementations, D(t) is quantized. The microcontroller outputs the PWM voltage $VPWM,M(t)=VMsgn(D(t))H(|D(t)|−c(t))$ to drive the H-bridge, where V_M is the microcontroller supply voltage (constant), H is the Heaviside unit step function defined such that H(0) = 0, and $c(t)∈[0,1]$ is the carrier signal for the selected PWM mode. The signal c(t), which is commonly a sawtooth or triangular wave, has the same frequency as D(t) such that H(|D(t)| − c(t)) is a series of pulses of width |D(t)|T_PWM.¹² The duty cycle D(t) then determines the switching state of the H-bridge as a function of time, regulated by V_PWM,M(t) as communicated through the microcontroller. Although the state equations are derived here specifically for a unipolar leading-edge PWM mode, where an on-state of variable width precedes the off-state with a zero voltage,⁶⁵ the same modeling process can be followed when considering a general servomotor that uses any variations of switching configurations for control.

In addition to the on- and off-states, there exist intermediate switching states where only switch S2 or switch S4 is closed (Table I). These intermediate switching states occur during dead time, also known as blanking time, and are introduced by design to avoid shoot-through when a short-circuit is formed by accidentally closing switches S1 and S2 or S3 and S4 together while toggling between the on- and off-states.¹³ Dead time occurs twice during a pulse period with a fixed dead-time duration T_DT each instance. Thus, for a given pulse beginning at t = t₀, the pulse period interval [t₀, t₄) is divided into the first dead time [t₀, t₁), on-state [t₁, t₂), second dead time [t₂, t₃), and off-state [t₃, t₄), where t₁ = t₀ + T_DT, t₂ = t₀ + |D(t₀)|T_PWM, t₃ = t₀ + |D(t₀)|T_PWM + T_DT, and t₄ = t₀ + T_PWM (Fig. 2).

Ideal instantaneous switching is assumed, where the time to turn the switch on and off, which are the rise and fall time, respectively,¹² is assumed to be zero. Commonly in H-bridges, each switch is a metal–oxide–semiconductor field-effect transistor (MOSFET), which is modeled as a switch with an on-state resistance R_trans in parallel with the body diode of the transistor (Fig. 1). The diode, which is referred to as a freewheel or catch diode in the context of the H-bridge, provides a path for I_a to flow during dead time in order to avoid voltage spikes due to the motor inductance and to protect the transistors and is modeled as an ideal diode in series with the forward voltage V_D and diode resistance R_D.¹³ 

The H-bridge state equations during the on-state, off-state, and dead time are derived from Kirchhoff's voltage and current laws (KVL and KCL) considering the equivalent circuit model of the H-bridge in each switching configuration (Table I and Fig. 3).

For on-state (⁠$t∈[t1,t2)$⁠) and off-state (⁠$t∈[t3,t4)$⁠):

For dead time (⁠$t∈[t0,t1)∪[t2,t3)$⁠):

During dead time, the sign of I_a also determines the path of I_a. For instance, when D(t) > 0 and I_a(t) < 0, current is conducted through switch S4, the diode in parallel with switch S1, and the power supply with non-zero I_in(t) (Fig. 3). If D(t) > 0 and I_a(t) > 0, current is conducted through switch S4 and the diode in parallel with switch S2 but not through the power supply [I_in(t) = 0] because the diode in parallel with switch S3 is reverse-biased.¹³ 

The coupled state equations of the servomotor electromechanical dynamics are the DC motor and gear train equation of motion and the H-bridge circuit equations, with the voltage in terms of V_in and D(t).

On-state (⁠$t∈[t1,t2)$⁠) and off-state (⁠$t∈[t3,t4)$⁠):

Dead time (⁠$t∈[t0,t1)∪[t2,t3)$⁠):

At the boundary of the servomotor system of interest for this study, the mechanical variables are q(t), $τout(t)$⁠, and their time-derivatives, and the electrical variables are D(t), I_in(t) [dependent on I_a(t) and D(t)], and their time-derivatives, while V_in is constant.

In the proposed approach, the instantaneous state variables are cycle-averaged,^12,13,26,57 as in a typical circuit-averaging approach, and their relations are derived. Given a function f(t) of time $t∈[t0,t4)$⁠, which may represent a state or control, its cycle-averaged value $f¯(t)≡1TPWM∫t0t0+TPWMf(t)dt$ can be regarded as a constant function over the given pulse period (Fig. 4). Here, the set of instantaneous and cycle-averaged f(t) values sampled at discretized time points across various pulse periods, which can be non-consecutive or unevenly spaced, are denoted by f and $f¯$⁠, respectively. For sufficiently high 1/T_PWM, the mechanical state trajectories $q≈q¯$⁠, $q˙≈q˙¯$⁠, $q¨≈q¨¯$⁠, and $τout≈τ¯out$ because their variations over the pulse period are negligible relative to the variations of the electrical state variables, which typically have a much smaller time constant than that of the mechanical state variables.^58,64 Due to the quantized nature of D(t) in digital control, $D=D¯$ exactly (Fig. 2).

In this framework, the discretized mechanical state variable trajectories q (and its derivatives) and $τout$ are mapped to the corresponding electrical state variable trajectories $I¯in$ and D. First, the set of instantaneous coupled dynamics [Eqs. (3) and (6)] are cycle-averaged to obtain the state equations for $I¯a$⁠, which are derived without the consideration of dead time due to its short [<5% T_PWM by design in order to minimize the dead-time distortion of V_PWM,H(t) during operation¹³] duration relative to T_PWM:

where the approximation ${I}¯˙a≈I˙¯a$ is used to interchange the order of differentiation and cycle-averaging operations.²⁷ Given $q˙(t)≈q˙(t0)$ and D(t) = D(t₀) for $t∈[t0,t4)$⁠, it can be shown that $τf(t0)≈b0sgn(q˙(t0))+b1q˙(t0)$ and $sgn(D(t))H(|D(t)|−c(t))¯=1TPWM∫t0t0+|D(t0)|TPWMsgn(D(t0))dt=D(t0)$⁠. In these cycle-averaged coupled state equations for $I¯a$⁠, the effect of ripple current appears only in the term $Vbrsgn(Ia)¯$⁠.

The state equation for $I¯in$ is derived in terms of $I¯a$ to include the effects of ripple current. Here, the dead time is included to account for instances of low |D(t) |, in which case T_DT is not negligible with respect to the on-state duration $t∈[t1,t2)$⁠. For the case of V_br = 0, the instantaneous coupled dynamics [Eqs. (6) and (7)] are those of switched RL circuits and the corresponding analytical solution for I_a(t) for $t∈[t0,t4)$ is a piecewise function during the on-state, off-state, and dead time.

For dead time (⁠$t∈[t0,t1)$⁠):

For on-state (⁠$t∈[t1,t2)$⁠):

For dead time (⁠$t∈[t2,t3)$⁠):

For off-state (⁠$t∈[t3,t4)$⁠):

with the values of I_a(t₁), I_a(t₂), and I_a(t₃) satisfying the continuity between the piecewise solutions and written in terms of I_a(t₀).

Due to the brief duration of dead time, the signs of I_a(t) during $t∈[t0,t1)$ and $t∈[t2,t3)$ are assumed to remain constant and equal to their values at t₀ and t₂, respectively, resulting in four possible combinations for the signs of I_a(t₀) and I_a(t₂). The correct combination of signs is identified when the signs of I_a(t₀) and I_a(t₂) resulting from I_a(t) agree with their pre-assumed combination.

Using the analytical forms of I_a(t), $I¯a(t)$ for $t∈[t0,t4)$ can be written as a first-order form with respect to I_a(t₀):

where $ϕ1$ and $ϕ2$ are functions of D(t₀), $q˙(t0)$⁠, and system parameters (G, K, L, R_a, R_D, R_trans, V_in, V_D, T_DT, and T_PWM), obtained through analytical integration depending on pre-assumed signs of I_a(t₀) and I_a(t₂).

From the instantaneous relations [Eqs. (4) and (5)] between I_in(t) and I_a(t) during the on-state and dead time and the analytical solution of I_a(t), for $t∈[t0,t4)$⁠:

where $ϕ3$ and $ϕ4$ are functions of D(t₀), $q˙(t0)$⁠, and the system parameters, and depend on pre-assumed signs of I_a(t₀) and I_a(t₂). At discretized time points across pulse periods, using Eqs. (8) and (14),

For the case of non-zero V_br, the process of obtaining I_a(t) requires solving for the times $t2−∈[t0,t2]$ and/or $t2+∈[t2,t4)$ satisfying I_a(t₂⁻) = 0 and/or I_a(t₂⁺) = 0, respectively, when the sign of I_a(t) changes. If no sign change occurs during $t∈[t0,t2]$ and/or $t∈[t2,t4)$⁠, then these are set as t₂⁻ = t₀ and/or t₂⁺ = t₄, respectively. Note that if the sign of I_a(t) never changes within the period, then t₂⁺−t₂⁻ = T_PWM. For each possible case of t₂⁻ and t₂⁺, the proper piecewise form of I_a(t) is used. For example, if the sign of I_a(t) changes both during the on- and off-states (i.e., $t2−∈[t1,t2)$ and $t2+∈[t3,t4)$⁠), I_a(t) consists of six segments defined over the time intervals [t₀, t₁), [t₁, t₂⁻), [t₂⁻, t₂), [t₂, t₃), [t₃, t₂⁺), and [t₂⁺, t₄) (Fig. 4). Solving for I_a(t) with the assumption of V_br = 0 provides an initial approximation for t₂⁻ and t₂⁺ as the times where I_a(t) = 0. The approximated t₂⁻ and t₂⁺ are used to re-evaluate the cycle-averaged state equations with non-zero V_br [Eqs. (8) and (9)] to update $I¯a(t)$ and D(t₀) in order to obtain I_a(t). Since sgn(I_a(t)) is equal to sgn(I_a(t₂)) for $t∈(t2−,t2+)$⁠, zero for t = t₂⁻ and t = t₂⁺, and −sgn(I_a(t₂)) for the remainder of $t∈[t0,t4)$⁠, $Vbrsgn(Ia(t))¯=Vbrsgn(Ia(t2)){2(t2+−t2−)/TPWM−1}$ for $t∈[t0,t4)$⁠. If needed, t₂⁻ and t₂⁺ of the updated I_a(t) can be used to update $I¯a(t)$ and D(t₀), and the process can be repeated until a desired accuracy is reached.

In contrast to typical circuit-averaging approaches, the cycle-averaging approach here comprehensively considers the effect of all switching configurations, including dead time. In addition, the effect of ripple current is incorporated by analytically obtaining $I¯in$ from $I¯a$ at any set of discretized time points, which requires less computational effort than numerically solving the high-frequency switched dynamics at a denser time-discretization, without requiring additional relations (e.g., energy conservation) or approximate ripple reconstruction algorithms. On the other hand, the analytical solution I_a(t) may not be continuous for consecutive pulse periods as a result of the lack of continuity conditions imposed and the simplifying assumptions applied to obtain the approximated $I¯a$ from the cycle-averaged state equations. Nevertheless, the effect of this potential discontinuity is negligible when the current is cycle-averaged.

The complete switched models of electromechanical dynamics allow for the formulation of servomotor power consumption (SPC) as the input electrical power in terms of mechanical state variables. In the subsequent energy analysis, the servomotor system (excluding the microcontroller and sensors) is identified as the system of interest and SPC is the input electrical power I_inV_in from the power supply. The changes of system energy due to mass transfer across the system boundary are assumed to be negligible during motor operation, and thermal dissipation is the primary form of dissipation considered. From the set of coupled state equations, SPC is derived as

where the thermal dissipation rate $h˙(t)$ depends on the switching configuration:

For on- or off-state (⁠$t∈[t1,t2)∪[t3,t4)$⁠):

For dead time (⁠$t∈[t0,t1)∪[t2,t3)$⁠):

with $Ia(t)=(1/GK){Jq¨(t)−τf(t)+τout(t)}$⁠. Consequently, SPC is a function of only the mechanical state and control of the servomotor. The SPC model is consistent with and can be analyzed with respect to the first law of thermodynamics for a closed system.^66–68 The first law of thermodynamics provides an independent viewpoint through the overall breakdown analysis and identification of the system's various energy terms (potentials, work, and dissipation) in the SPC model and the transfer and transformation between them. Its terms are the rate of change of total kinetic energy $Jq˙q¨=ddt(∑i=1NG12Jiq˙G,i2)$ of the system, including the rotor, gear train, and output shaft; output mechanical power $τoutq˙$ applied by the servomotor on its environment; and the system's internal magnetic potential energy storage rate $LIaI˙a=ddt(12LIa2)=LG2K2(Jq¨−τf+τout)(Jq⃛−τ˙f+τ˙out)$ for the inductor. Any change in the system's potential energy due to external conservative forces, such as gravity or elasticity, will be reflected by the output mechanical power term. The thermal dissipation term is composed of losses due to the (negated) friction work, Joule resistive heat, brush contact, and diode voltage, each of which is always positive to indicate the unidirectional power loss as outbound irreversible dissipation. Accordingly, while all terms in the SPC model are coupled with each other through the mechanical state and control variables, the direct effects of dissipated heat on the other model terms and parameters are assumed to be negligible, which can be ensured by monitoring the internal thermal energy rates through motor temperature in experiments (see below). Unlike existing DC motor electrical power models in the literature,^16,15,43 the inclusion of the entire servomotor components within the system of interest and the incorporation of their dynamic interactions across switching configurations provide more accurate and comprehensive energy terms and their cause-and-effect relationships.

The SPC model can serve as the cost function to resolve the redundancy in robot tasks while minimizing energy usage. The constrained nonlinear optimal control (or trajectory optimization) problem for minimum SPC with the solution time interval $t∈[0,T]$ is formulated as

subject to system dynamics that are cycle-averaged [Eqs. (8), (9), and (16)] and other constraints based on design and task requirements. The optimization variables vector x may include state variables and/or control inputs, depending on the given problem, and is bounded such that the feasible set is realizable by the actual servomotor. Each variable and its time-derivatives, if needed, can be time-parameterized for numerical implementation and solutions.

The system parameters of a typical servomotor can be obtained from a combination of manufacturer specifications, experimental measurements, and estimation techniques. To ensure consistency in the estimated parameter values, experiments were conducted from fixed initial conditions to reduce parameter variation resulting from history dependence on factors, such as current and temperature. For the MX-28 servomotor considered in this study, the parameters obtained from manufacturer specifications include the inductance L of the RE-17 DC motor (Maxon Motor AG, Switzerland; 17 mm diameter, precious metal brushes, 4 W, part number 214897);⁶⁹ the motor torque constant K;⁶⁹ the diode forward voltage V_D and on-state transistor resistance R_trans of the FDS6990A N-Channel Logic Level MOSFETs (Fairchild Semiconductor, Inc., USA) used as switches within the H-bridge;⁷⁰ the dead-time duration T_DT programmed in the IR2104 half-bridge drivers (Infineon Technologies, Inc., USA);⁷¹ and the magnitude of the total gear ratio G,⁵⁹ where G < 0 due to the even number of gears, N_G = 6.

The pulse period T_PWM was measured with an oscilloscope. The diode resistance R_D was chosen as equal to R_trans such that the transistor resistance remains constant throughout its operation, as is typically assumed.¹³ The DC motor's armature resistance R_a, servomotor's capacitance C, and capacitance in parallel with the DC motor, which was confirmed to be negligible (<10 nF), were measured with a digital multimeter. All gear masses, diameters, and tooth counts were measured to obtain G_i and J_i to compute the equivalent moment of inertia J.

For the verification and estimation of the remaining system parameters, experiments were conducted where controlled q(t) or $q˙(t)$ trajectories were externally imposed on the servomotor. A dual-servomotor test hardware platform was developed, similar to an existing design¹⁴ but with modifications, where the MX-28 (tested) servomotor was connected to the DYNAMIXEL XM-430 (driving) servomotor through a shaft coupler (Fig. 5). The driving servomotor was mounted on a linear stage to allow for manual alignment of both servomotor output shafts prior to testing. All experiments were conducted at room temperature, and the motor temperature was monitored with the servomotor's built-in temperature sensor to ensure steady-state thermal dissipation conditions and no overheating of the motor. The test platform was clamped to the table to reduce vibrations such that the imposed driving or tested servomotor trajectories were followed with minimal disturbance.

The realized trajectory q(t), which may deviate from the imposed trajectory due to coupler effects, was measured by the tested servomotor's built-in encoder. The current was measured with an INA219 current sensor (Adafruit Industries, LLC, USA), placed in series between the power supply and the servomotor, connected to an Arduino Uno microcontroller (Arduino AG, Italy). The current of interest I_in(t) was obtained from the measured power supply current by subtracting the idle current that represents the sum of actuation-independent sources of current flow in the system, including the built-in encoder, the temperature sensor, and the microcontroller that were connected to the same power supply as the servomotor. The idle current was measured while no external torque was applied, and the servomotor was stationary; i.e., the servomotor torque was disabled, which is equivalent to D(t) = 0 so that $Ia=τrot=0$⁠. The control input D(t) to the tested servomotor was recorded and synchronized with the other state variable measurements in post-processing. The local proportional-integral-derivative (PID) controller of the servomotor quantized the desired D(t) trajectories in increments of 1/885.

For the estimation of $τout(t)$⁠, the tested servomotor was mounted on a cantilever beam with a pivot point at one end and pressing on a load cell at its other end (Fig. 5). The load cell was connected to a USB-interfaced Wheatstone bridge (Phidgets Inc., Canada). Using the normal force measured by the load cell, $τout(t)$ was obtained from moment equilibrium analysis. All experimental conditions for the tested and driving servomotor trajectories were selected such that there was no lifting of the cantilever beam during testing. To account for the effect of weight and any unavoidable coupler misalignment, a five-term Fourier fitted model of the load cell force, obtained at 0.015 rad increments of q throughout the servomotor's range of motion under static and unactuated conditions, was subtracted from the measured load cell force.

The brush-contact voltage drop V_br was measured by stalling the DC motor within the tested servomotor. The DC motor terminals were desoldered from the H-bridge circuit and directly connected to a benchtop power supply, and the stall current was measured for V_in from 0.3 to 1.5 V. The gear train friction coefficients b₀ and b₁ were obtained for both directions of rotation by disconnecting the DC motor from the gear train [i.e., $τrot(t)=0$] and measuring the average $τout$ exerted by the gear train during constant $q˙(t)$ trajectories imposed by the driving servomotor for various values of $q˙$⁠. In these experiments, the motor temperature sensor was used to ensure that there were no significant thermal effects¹⁸ on the measured brush-contact and gear train friction properties.

To illustrate the proposed approach in a robotic system, in which the mechanical power is transferred from each servomotor to a link through the output shaft, a pendulum was attached to the servomotor output shaft and used for additional prediction and optimal control tasks. The equation of motion of the pendulum is $τout=Jpq¨+Mgdpsin⁡q$⁠, where J_p is the moment of inertia of the pendulum about the servomotor output shaft axis, M is the pendulum mass, g is the gravitational acceleration, d_p is the distance between the center-of-mass of the pendulum and the axis of the output shaft, and q = 0 corresponds to the downward vertical position of the pendulum. For optimal control, the SPC $I¯inVin$ is evaluated directly from $x=[qq˙q¨D]$ and is time-integrated for the cost function. The electrical duty cycle D rather than the mechanical actuator torque $τrot$ is chosen as the control input to solve for D during evaluation of the system dynamics, which requires the inclusion of $q¨$ in x as an additional state variable. Based on the equation of motion and the cycle-averaged equation for D for the case of V_br = 0, the pendulum system dynamics when $q˙$ does not change its sign is

For cases where $q˙$ changes sign, a differentiable sigmoid function approximation of $τf$ should be used in the derivation.

The modeling assumptions were evaluated using the estimated values for the electromechanical system parameters of the MX-28 servomotor (Table II). The PWM frequency 1/T_PWM = 40 kHz (T_PWM = 25 μs) is sufficiently large relative to the no-load angular velocity of the servomotor (5.76 rad/s = 0.92 Hz⁵⁹) such that the use of the cycle-averaged state equations is warranted.¹² According to the manufacturer specifications, the dead-time duration T_DT varies within the range 400–650 ns and has a typical value of 520 ns.⁷¹ That the total dead time 2T_DT is 4.16% of T_PWM warrants its omission when computing D(t) and $I¯a(t)$ from the cycle-averaged state equations. The state equations derived with V_br = 0 are used here, as in common control schemes,⁶⁴ because its estimated value (−0.0076 V from the linear fit of V_in as a function of I_in) is too low to be accurately measured, relative to the power supply input voltage V_in = 12.17 V. The filtering capacitance measured across the DC motor terminals while connected to the H-bridge was minimal, justifying its omission in the model. Using the estimated system parameters, the accuracy and merits of the models are experimentally illustrated and validated through predictions and optimal control tasks. Here, all trajectories in the figures were denoised using a Gaussian filter with a window of 0.16 s.

Selected state variables were predicted for given desired input trajectories and validated against experimental measurements. The first task considered is dynamic braking (by dissipation through armature resistance¹¹) with the pendulum system (M = 0.214 kg; d_p = 0.06928 m; J_p = 1.221 g m²) driven by gravity, where the pendulum falls from rest with initial conditions q(0) = 4.0 rad, $q˙(0)=0rad/s$⁠, and $q¨(0)=0rad/s2$⁠, and the time duration T = 21.7 s (Fig. 6). The H-bridge remained in the off-state without switching. Thus, the trajectory of q(t) was predicted through the models using the given conditions of D(t) = 0 and I_in(t) = 0 at all times. Nevertheless, $I¯a$ was non-zero due to the back emf from $q˙>0$⁠. The magnitude of $τout$ exhibited the same pattern as that of $I¯a$⁠, and its braking effect resulted in $q¨(t)<0$ for t > 7.3 s and q = 5.7 rad at the final time. In comparison, the pendulum would reach the same angle at a much shorter time (t = 0.17 s) if the H-bridge was disconnected such that the DC motor circuit would be open (i.e., $τrot=I¯a=0$⁠), resulting in $q¨>0$ during the falling trajectory due to the positive gravitational torque overcoming the friction (Fig. 6). Since $q˙(t)>0$⁠, $τoutq˙$ was always negative, and its difference from the zero SPC was mostly dissipated as heat due to the non-zero $I¯a$ (discussed in more detail later).

The root-mean-square deviation of the predicted q(t) from the measurement is 0.127 rad or 7.2% of the range of the measured q(t) over the total time duration. Since the dynamic braking task does not involve switching, this discrepancy may be due to the propagation of uncertainties and errors in the estimated system parameters. This may include the typical use of identical K for both V_emf and $τrot$⁠, which is based on the common assumption of ideal power conversion within a DC motor.^57,64 Also, in numerical investigations with a combination of different velocities and friction coefficients, the model predictions showed that, due to the low pendulum angular velocity involved, the q(t) trajectory is highly sensitive to slight variations in the friction coefficients b₀ and b₁, which affect both the acceleration and the final steady-state angle of the pendulum.

The second prediction task is a sinusoidal trajectory following, where a set of sinusoidally varying D(t) trajectories were applied on the tested servomotor as the driving servomotor imposed the sinusoidally varying velocity trajectory corresponding to $q(t)=cos⁡(πt/3+3π/2)−πt/3rad$ with PID position feedback control. The two desired D(t) trajectories tested were $D(t)=±{20sin⁡(πt/3+3π/2)/885+30/885}$⁠, which correspond to net negative and positive SPC, respectively (Fig. 7). The frequencies of both D(t) trajectories were equal to that of $q˙(t)$⁠, while their phases were equal and opposite to that of $q˙(t)$⁠, respectively. These D(t) trajectories were chosen to allow for a periodic continuous testing over a range of accelerations, while avoiding sudden changes in the direction of the force applied to the load cell that may cause vibration or lifting of the cantilever beam. With these desired D(t) trajectories, the trajectories of $I¯in(t)$ were predicted through the models using the measured trajectory of q(t) and the estimated $τout(t)$⁠. Due to the low magnitudes of the D(t) trajectories chosen, the magnitude of I_in(t) was significantly less than that of $I¯a(t)$ (Fig. 7). The root-mean-square deviation of the predicted $I¯in(t)$ from the measured I_in(t) is 13.4% (0.0049 A) and 16.4% (0.0035 A) for the net negative and positive SPC trajectories, respectively, of the range of the corresponding measured I_in(t) over the total time duration. The causes of these discrepancies may include fluctuations in the idle current drawn by the microcontroller and sensors, and the assumption of ideal switching, which neglects switching losses in the transistors of the H-bridge and the microcontroller, nonlinearity of T_DT, and variation of I_a in the DC motor during non-instantaneous switching.¹² Also, since I_in during the experiments was relatively low (<0.05 A) as compared with the stall current 1.4 A of the servomotor,⁵⁹ the assumption of constant parameters without considering the nonlinearity of the DC motor^64,72 and the H-bridge¹¹ across the full operating range of current may be another source of error. In addition, the variation of reaction forces applied by the coupler on the output shaft in the dual-servomotor test platform may have contributed to errors in the $τout(t)$ estimates.

For both tasks, the overall similar trends and the calculated errors between the predicted and measured results confirm the validity of the proposed comprehensive models of servomotor electromechanical dynamics with switching, the cycle-averaging approach, and the estimated system parameters. The results also validate the bi-directional use of the cycle-averaging approach, i.e., mapping from electrical to mechanical (for the dynamic braking) and from mechanical to electrical (for the sinusoidal trajectory following) state trajectories. The use of general electromechanical system parameters, as opposed to custom parameters in specialized models meant to fit a specific state variable (e.g., power consumption) or an operating condition, results in their validity across a wide range of task conditions, including cases of positive, zero, and negative SPC during $τoutq˙<0$⁠, which is often ignored in the literature.^5,43

The optimization variables used for the optimal control tasks with the pendulum system were $x=[qq˙q¨D]$⁠. The initial and final conditions q(0) = 0 rad, $q˙(0)=0rad/s$⁠, $q¨(0)=0rad/s2$⁠, and $q(T)=3π/2rad$ were imposed as constraints, along with $q˙≥0$ (for monotonic motion) and $−1≤D≤1$⁠, over the time duration T = 10 s. The final velocity and acceleration were not specified. The optimal control algorithm was implemented with the open-source trajectory optimization library OptimTraj⁷³ in MATLAB (The MathWorks, Inc., USA). Orthogonal collocation with Chebyshev polynomials was used to discretize the optimization variables, and the cost and constraint function gradients were computed numerically with built-in finite-difference methods.

For the cost function of the time-integrated SPC, the cycle-averaged $I¯in$ [Eq. (16)] in terms of the optimization variables was used as a low-frequency function approximation. In addition, other cost function models that are commonly used as proxies for power consumption^{43,45,48–50,52} are also considered for comparison: the square of the rotor torque (SRT) $τrot2$ and the positive rotor mechanical power (PMP) $max(τrotq˙rot,0)$⁠, which were cycle-averaged as $K2I¯a2$ and $max(GKq˙I¯a,0)$⁠, respectively, and time-integrated. (Here, the PMP was calculated at the rotor, rather than the output shaft as is typical,^43,52 so that the gear train friction directly affects its value.) The resulting optimal q(t) from each optimization was imposed as the desired trajectory for the PID position feedback control in experiments, where the gains were tuned manually such that the measured D(t) trajectory matches that from optimization. The resulting D(t) trajectories were similar in pattern to $I¯in(t)$ and SPC (Fig. 8) and had peak values of 0.20, 0.36, and 0.16 for the SPC, SRT, and PMP optimal trajectories, respectively. The results are analyzed with respect to two perspectives of the use of the proposed model—its merit in optimal control for power consumption minimization and its accuracy for further validation.

The time-integrated model SPC predicted from measured q(t) and D(t) was smaller than those directly measured by 9.1%, 6.0%, and 8.5% for optimal SPC, SRT, and PMP, respectively. The deviation of the measured D(t) from that of the optimization solution, which was due to the tracking error inherent in the PID position feedback control, contributed to the discrepancy between the SPC from experiments [directly measured or from measured q(t) and D(t)] and from optimization. The aforementioned issues of switching losses and other switching nonlinearities may also have contributed to the difference between the predicted and measured SPC, especially for (near) zero I_in(t). For all cost functions, the predicted $I¯in(t)$ reached zero (<1 mA)—for t > 7.4 s for optimal SPC, for 5.4 s < t < 8.8 s for optimal SRT, and near the final time for optimal PMP—while the measured I_in(t) did not (>5.6 mA for optimal SPC, >2.2 mA for optimal SRT, and >5.0 mA for optimal PMP). Overall, as opposed to applying D(t) = 0 for the continuous off-state without switching in the aforementioned dynamic braking task, which resulted in the predicted $I¯in(t)$ and measured I_in(t) to be exactly zero, the proposed models resulted in zero $I¯in$ during falling as regulated by optimization (Figs. 8 and 9). The overall accuracy illustrated serves as another validation of the proposed cycle-averaging approach and the SPC model across different tasks.

The measured SPC integrated over $t∈[0,T]$ was the lowest for the SPC optimal trajectory, which validates that the SPC predicted from the model best reflects the true SPC, as opposed to the other proxies pursued in the literature (Table III). For all cost functions, the measured SPC and $τout$ trajectories peaked approximately when $q(t)=π/2rad$⁠, where the gravitational torque was the largest (Figs. 8 and 9). While $τout(t)≈Mgdpsin⁡q(t)$ depends mainly on the angle (due to the low moment of inertia of the pendulum) and its peak values are almost identical across all cost functions, the variations of the SPC and $τout$ trajectories (including the timing of each peak) reflect each cost function model's distinct behavior. The optimal SRT stayed near the initial equilibrium (within 0.1 rad for t < 1.3 s), then quickly reached the upright equilibrium $q=πrad$ with steep changes of the velocity where $τout$ is near zero, and remained there for most of the time duration (within 0.1 rad for 3.4 s < t < 7.4 s) until around the final time before falling suddenly starts. In this way, the optimal SRT minimized the time spent near angles $q=π/2rad$ and $q=3π/2rad$ where $τout$ peaked (Fig. 9). The lack of velocity-dependence in the SRT cost function resulted in mostly static poses punctuated by sudden motion, as opposed to the smoother trajectories for the other cost function models. As a result, the output mechanical power of the optimal SRT suddenly peaked with relatively large values and then stayed zero during most of the time, which contributed to the distinct time-profile of the measured SPC that is mostly (near) zero with significant increases around $q=π/2rad$ and toward $q=3π/2rad$⁠. The optimal SPC resulted in much larger actuations during most of the time but smaller time-integrated SPC when compared with that of the optimal SRT, which is additional evidence of the distinct behaviors between the two cost function models. The similarity of the SPC and PMP optimal trajectories reflects the similarity in pattern between their SPC and $τoutq˙$⁠, which is approximately proportional to PMP (⁠$Gτrot≈τout$ due to the low moment of inertia and friction in the servomotor) during $τoutq˙>0$ (t < 6.0 s) before falling (Figs. 8 and 9). On the other hand, by construction, PMP is unaffected by variations in SPC when $τrotq˙rot<0$ during the falling phase, and consequently, the time-integrated SPC of the optimal PMP was greater than that of the optimal SPC.

The falling phases (⁠$τoutq˙<0$⁠) demonstrated further notable differences. The measured SPC during falling significantly increased for the optimal SRT, gradually reached (near) zero only at the end of falling for the optimal PMP, and stabilized at zero shortly after falling began for the optimal SPC. That the time t = 7.4 s at which SPC reaches nearly zero after passing the upright position [$q(t)≥π$ and measured SPC < 0.1 W] for the optimal SPC was in between that of the optimal SRT and optimal PMP may reflect the SPC formulation as being mainly the sum of $h˙$ and $τoutq˙$⁠, which are loosely associated with the SRT and PMP models, respectively. Since falling can be achieved solely by gravity given sufficient time, it is unnecessary to have any positive SPC during falling from the perspective of minimizing SPC, which was noticeably demonstrated by the optimal SPC but not by the other two cost function models.

Overall, for optimal SRT, the attributes of the SPC time-profile were significantly distinct from those for optimal SPC (Fig. 8). For optimal PMP, the time-integrated SPC was the largest (Table III), partially due to its lack of capability to take advantage of the gravity during the falling phase. These results indicate that, with respect to their patterns and magnitudes, these existing proxy models are not valid as representations of the actual power consumption. For the purpose of power consumption minimization, the proposed SPC model provides true optimality, both instantaneous and time-integrated, as compared with the proxy models.

Another merit of the proposed mathematical models of the electromechanical system, along with their thermodynamic analysis,⁴⁶ is that they provide the breakdown of SPC into its constituent components (e.g., work vs heat, Joule heating from armature resistance vs frictional heating) and their distribution across time (Fig. 10), which is usually not possible through experiments (e.g., calorimetry⁷⁴). These types of thermodynamic analyses are also not available in other approximation approaches used mainly for control. The inductor's magnetic potential storage energy rate was approximated as $LI¯aI¯˙a$ through a finite-difference method.²⁷ The thermal dissipation, which is always non-negative, was obtained from $h˙≈I¯inVin−Jq˙q¨−τoutq˙−LI¯aI¯˙a$ as the remainder of the SPC after accounting for all non-dissipative components. [Calculating $h˙$ directly from $I¯a$ was avoided since it does not account for the effect of ripple current on the resistor and diode power losses^12,13 and may also introduce further error due to the lack of continuity conditions in deriving the piecewise analytical solutions of I_a(t).] For all cost functions, the relatively low proportion of the kinetic energy change rate is due to the choice of the servomotor as the system of interest, in which the rotor, gear train, and output shaft are the only moving parts. The energy stored or released by the inductor was also minimal, in line with other prior works where the motor inductance is typically neglected.^{5,7,8,16,15,43} The output mechanical work was relatively small (<5%) as a share of the time-integrated SPC, with the remainder mostly dissipated as heat (Table IV). While $τoutq˙$ exhibited a similar trend as that of $h˙$ during the rising phase of the pendulum trajectory, they had opposite trends during the falling phase, which resulted in $h˙$ exceeding the measured SPC. The opposite trend is mainly due to $I¯a$⁠, not $I¯in(t)$⁠, that is (re-)generated by the back emf of the DC motor during the negative mechanical-work phase such that $h˙>0$⁠. The measured SPC was zero or positive (Fig. 8) even when both $I¯a(t)$ and $τout(t)q˙(t)$ were negative during falling (Fig. 9), which implies the distinct roles of $I¯a(t)$ vs $I¯in(t)$⁠. [Note that it is possible to achieve (near) zero $I¯a$ even when I_a(t) can still fluctuate to result in Joule heating from armature resistance or with non-zero D(t) by externally imposing q(t) and $τout(t)$ trajectories to induce (near) zero $τrot$⁠.] It can be seen (Fig. 10) that this aspect was exploited by the optimal SPC during the falling phase when the measured SPC was near zero, which was not noticeably demonstrated by the other cost function models. On the other hand, when zero measured SPC was maintained during the upright equilibrium of the optimal SRT, $τoutq˙$ and $h˙$ were also (near) zero due to zero $I¯a$⁠.

In all of the experiments conducted, the measured power supply current was always positive even during negative SPC due to the measured current being the sum of I_in and the idle current. In general, the case of negative SPC requires special consideration because, while the DC servomotor is inherently a bi-directional energy converter, the power supply is not.^11,16 Power supplies, such as a battery, can be either rechargeable or non-rechargeable, and even in the case of a rechargeable battery, the efficiency of a battery in absorbing power flows depends on battery-specific parameters.

The integration of the coupled models for switched electromechanical dynamics of a typical DC servomotor in the proposed approach, which were experimentally validated for their predictive accuracy and demonstrated for their merits against common proxy models of power consumption, allowed for the analysis and control of general aperiodic tasks, including transient phases. These enhancements of the model were provided by its comprehensiveness in including the switching configurations of on-state, off-state, and dead time, and the combination of cycle averaging with piecewise analytical solutions to handle the non-smooth dynamics and different temporal scales from high-frequency electrical to low-frequency mechanical variables. Future studies will investigate switching losses, idle current, non-ideal power supplies (e.g., batteries), variations of the system parameters, and the scalability of the model for large current/power and multi-degree-of-freedom robotic systems. In addition, validity of the cycle-averaging method under different conditions needs to be further evaluated.

This work was supported in part by the U.S. National Science Foundation (Nos. CMMI-1436636 and IIS-1427193) and a Mitsui USA Foundation scholarship.