HIL Implementation of Harmonic Drive Motor (Part I)

Part I: Discretization and Simulation

One approach quite common in industrial control nowadays, when doing analysis and design with a real plant is expensive, is to implement the plant model on a computer, usually an embedded system. Some advantages of this so-called Hardware-In-the-Loop (HIL) approach, in addition to lowering development cost and time, are flexibility in adjusting the parameters, and the ability to perform extreme experiments which might damage the real plant.

From our previous articles on robot joint driven by DC motor with harmonic drive mechanism (from now I will abbreviate it as HDM for short) , we obtain a linear model and simulate with PID and custom-designed controllers. Here we continue our discussion on how to implement this math model on a real-time embedded system, to be used as an HIL plant for feedback control analysis and design.

To be honest, this is my first HIL project that is based on a math model of a physical plant, except for some simple memoryless systems. Certain factors, both hardware and software related, such as ADC and PWM accuracy, delay, fixed-point math and other numerical issues, could affect the performance of this embedded system. With limited resources, at the end I am not so sure it could resemble the real plant good enough to warrant its place in a control lab. Nevertheless, it is worth a try. If things work out properly, it could be an interesting topic for my control engineering students for the upcoming semester. If not, I believe we learn something together in such attempt. 😉

The proposed procedure is as follows

  1. Discretize the continuous-time transfer function of the HDM
  2. Simulate the discrete-time plant and compare the responses with continuous-time model
  3. Adjust overall plant gain to suit the input and output limits of embedded system
  4. Implement on the selected microcontroller (in this case, a PIC24EP product from Microchip)
  5. Experiment with the HIL model

Part I of this material covers step 1 – 3. The rest is left to future article(s).


From the HDM block diagram derived from our previous article and displayed again in Figure 1, we can manipulate to achieve a reduced block diagram in Figure 2, consisting of only 2 sub-blocks

(1)   \begin{equation*}  P_1(s) = \frac{k_mk}{a_{c5}s^5 + a_{c4}s^4 + a_{c3}s^3 + a_{c2}s^2 + a_{c1}s} \end{equation*}

(2)   \begin{equation*}  P_2(s) = \frac{rP_l(s)}{k} \end{equation*}


(3)   \begin{eqnarray*}  a_{c5} &=& LJ_mJ_l  \nonumber \\ a_{c4} &=& RJ_mJ_l + L(J_mB_l + J_lB_m)  \nonumber \\ a_{c3} &=& L(J_mk + B_mB_l + J_lk) + R(J_mB_l + J_lB_m) + rk_mk_bJ_l  \nonumber \\ a_{c2} &=& Lk(B_m + B_l) +R(J_mk + B_mB_l + J_lk) + rk_mk_bB_l  \nonumber \\ a_{c1} &=& k(Lk + R(B_m + B_l) - kL + rk_mk_b) \end{eqnarray*}


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

Figure 1 block diagram of HDM plant

Figure 2 reduced block diagram of HDM plant

The reduction process is straightforward but quite tedious. We leave the details to the reader. Note that P_2(s) is not proper, so this diagram is not suitable for continuous-time simulation. However, we can use this structure conveniently in the discrete-time domain.

There are a couple of discretization methods. In our experiment we exploit the so-called bilinear, or Tustin method. The easiest way to implement is to rely on Scilab or MATLAB to discretize the plant for a particular set of parameters, and then hard-code the results in the C-program (or write some online coefficient updating function). The drawback is you need to redo the discretization offline each time a parameter is changed.

So here we take the pain to derive discrete-time transfer function coefficients as functions of physical parameters. By using the bilinear relationship

(5)   \begin{equation*}  s = \frac{2}{T}\frac{(z - 1)}{(z + 1)} \end{equation*}

with sampling period T, we eventually have the discrete-time representation of P_1 and P_2 as

(6)   \begin{equation*}  P_1(z) = \frac{k_mk(b_{1z5}z^5 + b_{1z4}z^4 + b_{1z3}z^3 + b_{1z2}z^2 + b_{1z1}z + b_{1z0})}{a_{1z5}z^5 + a_{1z4}z^4 + a_{1z3}z^3 + a_{1z2}z^2 + a_{1z1}z + a_{1z0}} \end{equation*}

(7)   \begin{equation*}  P_2(z) = \frac{r(b_{lz2}z^2 + b_{lz1}z + b_{lz0})}{k(a_{lz2}z^2 + a_{lz1}z + a_{lz0})} \end{equation*}


(8)   \begin{eqnarray*}  a_{1z5} &=& 32a_{c5} + 16Ta_{c4} + 8T^2a_{c3} + 4T^3a_{c2} + 2T^4a_{c1} \nonumber \\ a_{1z4} &=& -160a_{c5} - 48Ta_{c4} - 8T^2a_{c3} + 4T^3a_{c2} + 6T^4a_{c1} \nonumber \\ a_{1z3} &=& 320a_{c5} + 32Ta_{c4} - 16T^2a_{c3} - 8T^3a_{c2} + 4T^4a_{c1} \nonumber \\ a_{1z2} &=& -320a_{c5} + 32Ta_{c4} + 16T^2a_{c3} - 8T^3a_{c2} - 4T^4a_{c1} \nonumber \\ a_{1z1} &=& 160a_{c5} - 48Ta_{c4} + 8T^2a_{c3} + 4T^3a_{c2} - 6T^4a_{c1} \nonumber \\ a_{1z0} &=& -32a_{c5} + 16Ta_{c4} - 8T^2a_{c3} + 4T^3a_{c2} - 2T^4a_{c1} \nonumber \\ b_{1z5} &=& T^5 \nonumber \\ b_{1z4} &=&5 T^5 \nonumber \\ b_{1z3} &=& 10T^5 \nonumber \\ b_{1z2} &=& 10T^5 \nonumber \\ b_{1z1} &=& 5T^5 \nonumber \\ b_{1z0} &=& T^5 \nonumber \\ a_{lz2} &=& T^2 \nonumber \\ a_{lz1} &=& 2T^2 \nonumber \\ a_{lz0} &=& T^2 \nonumber \\ b_{lz2} &=& 4J_l + 2B_lT + kT^2 \nonumber \\ b_{lz1} &=& 2kT^2 - 8J_l\nonumber \\ b_{lz0} &=& 4J_l - 2B_lT + kT^2  \end{eqnarray*}

Once again, the routine work is left as an exercise.

Verification with Scilab

Before doing anything else, it makes sense to verify that our derivation of the discrete-time transfer functions (6), (7) and their coefficients (8) are correct. Here we show how to verify (6) using discretization command in Scilab. You can use similar process for (7).

Note :all Scilab commands listed below are included in the script file hdminit.sce

Given a set of parameters, we construct the continuous-time plant as before.

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
// plant denominator coefficient
ac5 = L*Jm*Jl;
ac4 = R*Jm*Jl+L*(Jm*Bl + Jl*Bm);
ac3 = L*(Jm*k + Bm*Bl + Jl*k) + R*(Jm*Bl + Jl*Bm) + r*km*kb*Jl;
ac2 = L*k*(Bm+Bl)+R*(Jm*k + Bm*Bl + Jl*k) + r*km*kb*Bl;
ac1 = k*(L*k + R*(Bm+Bl) - k*L + r*km*kb);
pl = Jl*s^2+Bl*s+k;
pm = Jm*s^2+Bm*s+k;
// transfer function from V to load angle
pnum1 = km*k;
pden1 = ac5*s^5 + ac4*s^4 + ac3*s^3 + ac2*s^2 + ac1*s;;
P1 = syslin('c',pnum1,pden1);

Executing the above script yields the continuous-time transfer function

(9)   \begin{equation*}  P(s) = \frac{33333.33}{s(s^2 - 0.0483s + 11022.313)(s^2 -0.9783s + 302.435)}  \end{equation*}

Specify the sampling period of discrete-time system, say, 10 milliseconds.

T = 0.01;

To use the Scilab command cls2dls to discretize the plant, first it must be described in state-space form

P1ss = tf2ss(P1);

Then put this plant and sampling time as arguments to cls2dls, which outputs a discrete-time state space system using Bilinear transformation.

Pd1ss = cls2dls(P1ss, T);

For comparison with our hand computation, we prefer a discrete-time transfer function. There is some other way to do this, but here we simply use the control textbook relationship P(z) = C(zI-A)^{-1}B + D .This is done in Scilab as follows

z = poly(0,'z');
Pd1 = Cd*inv(z*eye(size(Ad,1),size(Ad,2)) - Ad)*Bd + Dd;

Running the script in Scilab yields

 Pd1  =
                                       2            3            4            5  
    0.0000008 + 0.0000039z + 0.0000078z + 0.0000078z + 0.0000039z + 0.0000008z   
                                           2            3            4   5       
      - 0.9224925 + 3.8350285z - 6.9849487z + 7.1298534z - 4.0574407z + z

Now we construct another discrete-time plant from our hand calculation.

// numerator coefficients of p1(z) 
b1z5 = T^5;
b1z4 = 5*T^5;
b1z3 = 10*T^5;
b1z2 = 10*T^5;
b1z1 = 5*T^5;
b1z0 = T^5;
// denominator coefficients of p1(z) 
a1z5 = 32*ac5 + 16*T*ac4 + 8*T^2*ac3 + 4*T^3*ac2 + 2*T^4*ac1;
a1z4 = -160*ac5 - 48*T*ac4 - 8*T^2*ac3 + 4*T^3*ac2 + 6*T^4*ac1;
a1z3 = 320*ac5 + 32*T*ac4 - 16*T^2*ac3 - 8*T^3*ac2 + 4*T^4*ac1;
a1z2 = -320*ac5 + 32*T*ac4 + 16*T^2*ac3 - 8*T^3*ac2 - 4*T^4*ac1;
a1z1 = 160*ac5 - 48*T*ac4 + 8*T^2*ac3 + 4*T^3*ac2 - 6*T^4*ac1;
a1z0 = -32*ac5 + 16*T*ac4 - 8*T^2*ac3 + 4*T^3*ac2 - 2*T^4*ac1;
// -----  p1(z) -------------------
p1numz = km*k*(b1z5*z^5 + b1z4*z^4 + b1z3*z^3 + b1z2*z^2 + b1z1*z + b1z0);
p1denz = a1z5*z^5 + a1z4*z^4 + a1z3*z^3 + a1z2*z^2 +a1z1*z + a1z0;

To make a comparison, notice that the coefficient of highest z in the denominator of Pd1 is normalized to 1, but in p1denz it equals alz5. So we have to normalize both p1numz and p1denz with this quantity. Using syslin command for this discrete-time plant yields the same transfer function as the one computed by Scilab cls2dls command.

-->Pd1h = syslin('d',p1numz/a1z5,p1denz/a1z5)
 Pd1h  =
                                       2            3            4            5  
    0.0000008 + 0.0000039z + 0.0000078z + 0.0000078z + 0.0000039z + 0.0000008z   
                                           2            3            4   5       
      - 0.9224925 + 3.8350285z - 6.9849487z + 7.1298534z - 4.0574407z + z

and this should make us happy for the moment.

Step Response Simulation

To see how well the resulting discrete-time transfer functions (6) and (7) represent their continuous-time counterparts, we must simulate the two and compare their responses. Before running any simulation model below, first you need to execute the script hdminit.sce to initialize all the variables used in the diagram.

Figure 3 shows an Xcos diagram pid_lf_cdcompare.zcos, a PID control comparison with load angle as feedback. The continuous time and discrete time plants are placed in the upper and lower feedback loops, respectively. So the upper loop is purely continuous system, while the lower loop is a hybrid system (continuous-time PID controller, discrete-time plant).

Figure 3 [pid_lf_cdcompare.zcos] PID control with load angle feedback

Note also that an LPF transfer function, with pole well above closed-loop bandwidth, is inserted in the feedback path of the hybrid system. This effectively prevents algebraic loop error during simulation.

This LPF does have its usage in our real implementation later on. The output variables may be sent out as PWM signals, which require simple RC circuits as LPF before feeding to an A/D of the controller.

Running the simulation with PID gains: k_p = 4, k_i = 1, k_d = 0.5, yields the step responses in Figure 4. The responses from two feedback systems closely matches, though that of the discrete-time plant (red) tends to oscillate more. This oscillation is more pronounced when the controller gains are increased.

Figure 4 step responses from model in Figure 3

We also construct a comparison model for the case of using motor angle feedback, as shown in Figure 5. The responses from PID gains: k_p = 5, k_i = 2, k_d = 1 are shown in Figure 6.

5 [pid_mf_cdcompare.zcos] PID control with motor angle feedback

Figure 6 step responses from model in Figure 5

Adjusting Plant Gain to Embedded Hardware

So far we have developed a discrete-time model quite ready to be implemented on a microcontroller. There is one issue left to be considered, and is fundamentally related to the properties of hardware used.

We ask a simple question. How can you be sure that the plant gain is comparable with the range of input and output of your hardware? The extreme cases are too high gain that easily causes output saturation, or too small gain that requires too excessive input to drive the system than the controller could deliver. So we need to calibrate the plant gain to take advantage of the full input and output ranges.

For the plant transfer functions described by (6),(7), the plant gain is a complicated function of many parameters. In fact, due to the presence of integrator the DC gain is infinity. So if we simulate an open-loop step response, the motor does not settle at a value but continues to rotate. This makes gain adjustment somewhat tricky. Nevertheless, it is not too difficult when we realize that, in a real motor, a property that is already available as standard specs is the velocity, commonly specified in Revolution Per Minute (RPM).

In a simulation diagram, motor velocity can be computed by adding a derivative block as shown in Figure 7. Note in the diagram that we “measure” both \theta_l and \theta_m and their derivatives. It is however only one variable, in this case \theta_m, the motor shaft velocity, we want to calibrate, since two angle outputs are proportionally related by the gear ratio.

Figure 7 [speedtuning.zcos] Xcos diagram for open-loop gain adjustment

For calibration purpose, we add an amplifier gain block kv to the HDM subblock as in Figure 8. This in effect makes the overall plant gain independently adjustable.

Figure 8 adding an amplifier block k_v for overall gain adjustment

To show how to compute k_v, let us add some concrete parameter values. Suppose the resolutions of input ADC and output PWM are 10 and 12 bits, respectively. The maximum velocity rating for this motor is represented by S_{max} (RPM).

From Figure 2, with new gain block k_v added at the input, the transfer function from V to \theta_m equals

(10)   \begin{equation*}  \frac{\theta_m(s)}{V(s)} = \frac{1}{s}\left(\frac{k_vk_mr(J_ls^2 + B_ls + k)}{a_{c5}s^4 + a_{c4}s^3 + a_{c3}s^2 + a_{c2}s + a_{c1}}\right) \end{equation*}

Note that the term in parenthesis represents the transfer function from input voltage to motor shaft velocity. Let \delta_m represent output rate of change per second, then at DC we have

(11)   \begin{equation*}  \delta_m = \frac{k_vk_mrk}{a_{c1}}V \end{equation*}

At maximum input, we want the motor to rotate at maximum speed, and translate this to maximum output. We also have to take into consideration the motor sense, clockwise or counter-clockwise. So the range for each direction is half of the input range, say, for ADC input value of 511, the motor rotates at S_{max} RPM, and the PWM maximum output is 4095.

Substituting all these quantities and do some math, we have

(12)   \begin{equation*}  k_v = 0.1333\frac{a_{c1}S_{max}}{k_mrk} \end{equation*}

For example, a motor with maximum speed 3,000 RPM yields k_v = 400. With this gain and step input of 511, the output changes at rate 204,350 units/sec. The equation (12) is coded in hdminit.sce and used in Xcos model speedtuning.zcos. Figure 9 shows the simulation result of the output rate, when the step input is applied at t = 1 second. We see the fluctuation around the expected value.

Figure 9 simulation result showing the rate of motor angle output

Even more interesting open-loop load angle rate of change is shown in Figure 10. We see a large sinusoid fluctuation that dies out quite slowly.

Figure 10 simulation result showing the rate of load angle output

Let us perform further calculation from this setup. Suppose we assign the maximum PWM value of 4095 to represent the load angle of 360 degrees (full turn). With motor turned at maximum velocity 3,000 RPM, the load turns at 300 RPM, one tenth of the motor velocity. The PWM value for load angle will reach its maximum in about 200 milliseconds.

We will see in Part II whether this setup is suitable for our embedded system. Stay tuned.

Scilab Script and Model Files

The files used this document are

Note :execute the script file first before running simulation

Download all as zip file: hilhdm1.zip


  1. Show derivation details of continuous-time transfer functions and their coefficients (1) – (4).
  2. Show derivation details of discrete-time transfer functions and their coefficients (6) – (8).
  3. Derive (12).
  4. Experiment with Xcos model comparisons. Try changing sampling time, and PID gains. Can you explain their effects on the responses?



Comments are closed.