Linear Model of Robot Joint with Harmonic Drive

Harmonic drive refers to a type of gear mechanism often used in robots, due to its low backlash, high torque transmission, and compact size. It is developed by the Harmonic Drive group companies. It is explained in more detail in [1], as well as other sources from the internet including the manufacturers’ websites. This research paper is also helpful.

The structure of a typical harmonic drive consists of a rigid circular spline, a flexible flexspline, and an elliptical wave generator. The wave generator is coupled to the motor shaft and hence is rotated at high speed. The circular spline has internal teeth. Between the two are the flexspline that has teeth on the outside surface. As the wave generator turns, it deforms the flexspline causing its teeth to mesh with the teeth of circular spline. The effective gear ratio is determined by the difference in the number of teeth of the flexspline and circular spline.

In [1], a linear math model of a harmonic drive from the motor torque to the load angle is derived, without considering the actuator dynamics. The complete model is left as an exercise (problem 6-15). In this article we show how to develop an overall linear model of a DC motor equipped with harmonic drive mechanism, and then perform some computation and simulation with Scilab. Consult [2] for some basic discussion on using Scilab for independent robot joint control.

Overall Transfer Function

A robot joint actuated by DC motor with harmonic drive can be described as in Figure 1. The left side is just an equivalent electrical circuit of DC motor, and on the right is the mechanical subsystem. Note the harmonic drive mechanism is represented as a gear transmission through a flexible shaft with stiffness k. This is the model suggested in [1]. Our task is to include the motor dynamics and gear ratio to the equations already derived there.

Figure 1 a robot joint with DC motor and harmonic drive

From Figure 1, we write down the following dynamical equations

(1)   \begin{equation*}  J_l\ddot{\theta}_l + B_l\dot{\theta}_l + k \left(\theta_l - \theta_{mr}\right) = 0 \end{equation*}

(2)   \begin{equation*}  J_m\ddot{\theta}_{mr} + B_m\dot{\theta}_{mr} - k \left(\theta_l - \theta_{mr}\right) = \tau = k_mi \end{equation*}

(3)   \begin{equation*}  \theta_m = r\theta_{mr} \end{equation*}

(4)   \begin{equation*}  L\dot{i}+Ri = V - k_b\dot{\theta}_m \end{equation*}

where k_m and k_b are torque constant and back EMF constant, respectively. To simplify the notation a bit, we assume here that J_m and B_m are the inertia and friction of motor after transmission through the gear ratio. For real motor parameters, one must adjust these values to achieve a precise model.

Taking Laplace transform and rearranging terms yield

(5)   \begin{equation*}  p_l(s)\theta_l(s) = k\theta_{mr}(s) \end{equation*}

(6)   \begin{equation*}  p_m(s)\theta_{mr}(s) -k\theta_l = k_mi(s) \end{equation*}

(7)   \begin{equation*}  i(s) = \frac{V(s)}{Ls + R} - \frac{k_bs\theta_m(s)}{Ls + R} \end{equation*}

where p_l(s) and p_m(s) are defined as

(8)   \begin{equation*}  p_l(s) = J_ls^2 + B_ls + k \end{equation*}

(9)   \begin{equation*}  p_m(s) = J_ms^2 + B_ms + k \end{equation*}

Substitute (7) into (6)

(10)   \begin{equation*}  p_m(s)\theta_{mr}(S) - k\theta_l = \frac{k_mV(s)}{Ls + R} - \frac{k_mk_bs\theta_m(s)}{Ls + R} \end{equation*}

From (5)

(11)   \begin{equation*}  \theta_{mr}(s) = \frac{p_l(s)}{k}\theta_l(s) \end{equation*}

together with gear relation (3), substitute into (10) and rearrange terms

(12)   \begin{equation*}  \left(\frac{p_m(s)p_l(s)}{k} - k + \frac{rk_mk_bsp_l(s)}{k(Ls + R)}\right)\theta_l(s) = \left(\frac{k_m}{Ls + r}\right)V(s) \end{equation*}

It is straightforward for the reader to verify that the overall transfer function of this system equals

(13)   \begin{equation*}  \frac{\theta_l(s)}{V(s)} = \frac{k_mk}{p_m(s)p_l(s)(Ls + R) - k^2(Ls + R) + rk_mk_bsp_l(s)} \end{equation*}

To study a concrete example, we need to assign parameters to (13). Below we write a Scilab script file hmdrv.sce that creates (13) with some fictitious parameter values. Feel free to change them as desired.

// Dew  JUn 16,2014
// modeling of robot joint with DC motor and harmonic drive
km = 100;  // torque constant
kb = 1;  // back EMF constant
k = 1000;  // torsional stiffness of harmonic drive
r = 10;   // gear ratio
L = 0.1  // armature inductance
R = 1;  // armature resistance
Jm = 1;  //  motor inertia
Bm = 0.01;  // motor shaft friction
Jl = 3;  // load inertia
Bl = 0.05;  // load friction
pnum = km*k;
pl = Jl*s^2+Bl*s+k;
pm = Jm*s^2+Bm*s+k;
pden = (L*s+R)*pm*pl - k^2*(L*s+R) + r*km*kb*s*pl;
P = syslin('c',pnum,pden);
disp("P = ")
disp("Open-loop poles ");

Executing the file to yield the linear transfer function P. The script also computes the open-loop poles.

 P =    
                         2            3        4      5  
    1000060s + 4056.0005s + 3400.0801s + 3.008s + 0.3s   
 Open-loop poles    
  - 4.5241679 + 104.88968i  
  - 4.5241679 - 104.88968i  
  - 0.4891654 + 17.383778i  
  - 0.4891654 - 17.383778i

This is a stable plant with an integrator and two pairs of complex conjugate poles, one is under-damped. Observe the Bode frequency response in Figure 2 by issuing the command


Notice the magnitude peak from the under-damped poles. We will show later how this affected the feedback system.

Figure 2 Bode frequency response of (13) with a set of chosen parameters

System Block Diagram

For simulation purpose, it is more convenient to express this electro-mechanical system (13) as a block diagram in Figure 3.

Figure 3 block diagram of robot joint with DC motor and harmonic drive

This can be constructed as an Xcos model using standard palettes like shown in Figure 4. The parameters are initialized by the script file hmdrv.sce above. Note that the model is contained in a super-block, with input 2 reserved for disturbance. Output 1 and 2 correspond to load and motor angles, respectively. You may want to add a velocity output for a more complicated control scheme such as cascade velocity and angle feedback.

Figure 4 Xcos model of Figure 3

Simple PID Control

To investigate this linear block diagram, we start with PID feedback control schemes. Following the setup in [1], Figure 5 shows an Xcos model pid_hmdrvsim_mf.zcos, a closed-loop system with motor angle used as feedback to the PID controller. The step-response simulation in Figure 6 is the result from PID gains K_p = 2000, K_i = 80, K_d = 1500. Download the model at the bottom of this article (all contained in a zip file) and experiment by yourself to see if you could get a nicer closed-loop response.

It is interesting to note from Figure 6 that, while the motor angle (red) tracks the step command quite nice, the load angle (green) shows oscillatory response that barely dies out. This is reasonable. From the block diagram in Figure 3 or 4 we see that the transfer function from motor angle to load angle is effectively outside the feedback loop. When that under-damped mode is excited, the closed-loop system cannot compensate the undesirable oscillation.

Figure 5 PID control with motor angle feedback pid_hmdrvsim_mf.zcos

Figure 6 step response from PID with motor angle feedback

Exercise: try changing the stiffness constant k. Which system is harder to control, high or low torsional stiffness?

Now, we experiment with an alternate scheme. Figure 7 shows another PID control model pid_hmdrvsim_lf.zcos, with the load angle used as feedback. In this case, everything is inside the feedback loop. It turns out, however, that the added complexity causes difficulty to the PID control. We could observe from the Bode plot in Figure 2 that increasing the gain would cause undesirable effect from the flexible mode. Indeed, the sluggish step response in Figure 8 results from PID gains K_p = 5, K_i = 1, K_d = 2. Cranking up the gains easily causes instability. It seems like the only benefit from this scheme is now we can get rid of the steady-state oscillation.

Figure 7 PID control with load angle feedback pid_hmdrvsim_lf.zcos

Figure 8 step response from PID with load angle feedback


In this article we show how to derive a linear model for robot joint with harmonic drive mechanism. This model may be suitable for basic control system analysis and design. Research control literature exists where more advanced math models of harmonic drive are developed. They may contain nonlinear elements such as hysteresis.

Even for this simple linear model, it is shown that the oscillatory mode resulting from flexible coupling in the harmonic drive can be challenging for simple PID control. In upcoming articles, we will craft higher order controllers to improve stability and performance of this motion control system.

Scilab Script and Model Files

The files used this document are

  • hmdrv.sce : Scilab script file for parameter initialization and transfer function computation
  • pid_hmdrvsim_mf.zcos : Xcos model file for PID with motor angle feedback
  • pid_hmdrvsim_lf.zcos : Xcos model file for PID with load angle feedback

Note: execute the script file first before running simulation

All files are contained in this zip file:


  1. M.W.Spong, S. Hutchinson and M. Vidyasagar, Robot Modeling and Control. John Wiley & Sons. 2006.
  2. V. Toochinda, Robot Analysis and Control with Scilab and RTSX, Mushin Dynamics, 2014.




Linear Model of Robot Joint with Harmonic Drive — 1 Comment

  1. Pingback: Loopshaping Design for Robot Joint with Harmonic Drive | Scilab Ninja