Mathematical Model of Lego EV3 Motor

Some time ago I calculated Lego NXT motor parameters to use them in mathematical modeling of Lego motor. Also I used them to calculate controller constants for Segway-like robot.

When I bought Lego EV3 set, I wanted to port my Segway program to EV3. But I need to build EV3 motor model to calculate balance controller constants. In this post I described how I measured and estimated the motor's parameters.

Motor's mathematical model can be described by following equations:

$$\begin{equation} \left\{\begin{aligned} \displaystyle \dot\omega&=\frac{K_\tau\cdot{I}-B\cdot\omega-A_r-\tau_d}{J}\\ \displaystyle \dot{I}&=\frac{U-R_a\cdot{I}-K_b\cdot{\omega}}{L_a} \end{aligned}\right. \label{eq:eq1} \end{equation}$$
where:

I U	(A) (V)	armature current armature voltage
ω τd	(rad/sec) (N·m)	shaft's angular velocity shaft's load torque
La Ra Kτ Kb	(H) (Ω) (N·m/A) (V/(rad/s))	armature inductance armature resistance torque constant back electromotive force coefficient
J B Ar	(kg·m²) (N·m/(rad/s)) (N·m)	rotor's moment of inertia viscosity resistance coefficient dry friction force

Some of these parameters can be measured directly. I could measure $L_a$ and $R_a$ using RLC meter.
$R_a$ = 4.196 Ohm
$L_a$ = 4.94 mH
But $R_a$ in the equation $\eqref{eq:eq1}$ includes the full circuit resistance, for example internal resistance of batteries. It must be calculated from other measurements.

We will use following scenarios to calculate motor's parameters:

1. Analyzing deceleration curve, when the motor rotates freely without external load and stops due to internal friction. This measurement allow us to find relations between $A_r, B$ and $J$.
2. Measure rotation speed and the motor's current under load. In this case the motor rotates with a constant speed. This measurement allows us to estimate $B, K_b, R_a$ and $K_\tau$. 
3. Analyzing acceleration curve, when the motor starts rotating without external load. This measurement allows us to adjust parameters calculated on the previous step.

1. Deceleration curve 

 In this case the motor circuit is open, external load is missing, so the motor model becomes simpler ($I=0, \tau_d=0$):
$$\begin{equation} \left\{\begin{aligned} \displaystyle \dot\omega&=\frac{-B\cdot\omega-A_r}{J}\\ \displaystyle U&=K_b\cdot{\omega} \end{aligned}\right. \label{eq:eq2} \end{equation}$$
The second expression allows us to find $K_b$ if we know rotation speed. But in this case we cannot measure rotation speed and we will use only the first equation:

$$\displaystyle \dot\omega=\frac{-B\cdot\omega-A_r}{J}$$

The solution of this equation can be found using the initial conditions ($\omega(0)=\omega_0$):

$$\omega(t) = -\frac{A_r}{B}+ e^{\textstyle-\frac{B \cdot t}{J}}(\omega_0+\frac{A_r}{B})$$
Because we can measure motor position but not the rotation speed, we need to find motor position function. Integration of rotation speed function with initial condition $\phi(0)=0$ gives us the following expression:
$$\begin{equation} \phi(t) = -\frac{A_r\cdot t}{B}+ (1-e^{\textstyle-\frac{B \cdot t}{J}})(\omega_0+\frac{A_r}{B})\frac{J}{B} \label{eq:eq_phi} \end{equation}$$
Now we can use least square method to find constants, but some of them are depended on others.
So we should introduce additional independent constants:
$$\left\{\begin{aligned} \displaystyle T_1&=\frac{B}{J}\\ \displaystyle T_2&=\frac{A_r}{B} \end{aligned}\right.$$
After that we can rewrite the function $\phi(t)$ from $\eqref{eq:eq_phi}$:
$$\phi(t)=-T_2\cdot t+\frac{\omega0+T_2}{T_1}(1-e^{\textstyle-T_1 \cdot t})$$
I wrote a program that starts the motor and turning the power off after one second of rotation. It measures the motor position every 5 ms and writes them to a log-file. Then a wrote a SciLab script that uses Least Squares method to find dependencies between $A_r, B$ and $J$. The calculation results are following relations:
$$\begin{equation} \frac{B}{J}=0.4837581433546762\\ \frac{A_r}{B}=10.697523425732065 \label{eq:eq_decel1} \end{equation}$$
and graph that shows measured and calculated motor position.

Motor position graph: measured and calculated values

Now we have relations between $A_r, B$ and $J$, and we can use them in the next steps.

2. Measurement under load

In case that the motor is steady state ($\omega=const$ and $I=const$), then these differential equations are represented by the following linear equations:
$$\begin{equation} \left\{\begin{aligned} \displaystyle &K_\tau\cdot{I}-B\cdot\omega=A_r+\tau_d\\ \displaystyle &R_a\cdot{I}+K_b\cdot{\omega}=U \end{aligned}\right. \label{eq:eq_load} \end{equation}$$
The solution of this system can be found in form:
$$\left\{\begin{aligned} I &= \frac{U\cdot B + K_b (A_r+\tau_d)}{R_a \cdot B + K_b \cdot K_\tau}\\ \omega &= \frac{R_a (A_r+\tau_d) - K_\tau \cdot U}{R_a \cdot B + K_b \cdot K_\tau} \end{aligned}\right.$$
 We can measure $I, \omega, U$ and $\tau_d$. Both equations are linear and have $\tau_d$ as an argument. So, if we know $A_r$ we can find $B, R_a, K_\tau, K_b$ using results of measurements in two points.
These measurement gives us three parameters for each point: ($\tau_{d_1}, I_1, \omega_1$) and ($\tau_{d_2}, I_2, \omega_2$).
The system that contains results of two measurements:
$$\left\{\begin{aligned} I_1 &= \frac{U\cdot B + K_b (A_r+\tau_{d_1})}{R_a \cdot B + K_b \cdot K_\tau}\\ \omega_1 &= \frac{R_a (A_r+\tau_{d_1}) - K_\tau \cdot U}{R_a \cdot B + K_b \cdot K_\tau}\\ I_2 &= \frac{U\cdot B + K_b (A_r+\tau_{d_2})}{R_a \cdot B + K_b \cdot K_\tau}\\ \omega_2 &= \frac{R_a (A_r+\tau_{d_2}) - K_\tau \cdot U}{R_a \cdot B + K_b \cdot K_\tau} \end{aligned}\right.$$

Solution of this system can be found in following form:
$$\begin{equation} \left\{\begin{aligned} B &= \frac{I_1 (A_r+\tau_{d_2}) - I_2 (A_r+\tau_{d_1}) }{\omega_2\cdot I_1-\omega_1\cdot I_2}\\ K_b &= \frac{U(I_1- I_2) }{\omega_2\cdot I_1-\omega_1\cdot I_2}\\ R_a &= \frac{U(\omega_1- \omega_2) }{\omega_2\cdot I_1-\omega_1\cdot I_2}\\ K_\tau &= \frac{\omega_1 (A_r+\tau_{d_2}) - \omega_2 (A_r+\tau_{d_1}) }{\omega_2\cdot I_1-\omega_1\cdot I_2} \end{aligned}\right. \label{eq:eq_load2} \end{equation}$$

Measurement assembly

I used a special assembly (LDD model) to measure rotation speed and the current under load. The assembly lifts objects with known weight and recorded the motor positions. From these positions I calculate the rotation speed. Also I used a circuit to measure the motor current using a multimeter.  

Rendered assembly

The real assembly photo

NXT based assembly

This photo was taken on the first attempt of EV3 motor parameters measurement, when I used NXT brick with LeJOS to record motor positions because NXT LeJOS is some kind of "real-time" OS. But LeJOS for EV3 uses embedded Java that requires some time to warm-up and does not allow me to organize stable and precise measurement time intervals. Unfortunately, I found that EV3 and NXT have different electrical parameters so I had to repeat my measurements on EV3 using C++ instead of Java for my program.

I used a wheel from Lego 8051 MotorBike set. The groove when I put the thread has diameter 68 mm. I used scales to measure weight of lifted objects. Then the torque can be calculated:
$$\tau_d=g \cdot m\cdot r$$
Where $g$ is gravitational acceleration, $m$ is object mass, and $r$ is the wheel radius.

The measurement results

I made a number of measurements and put the results into a table:

Load torque, N*cm	Motor current, A	Rotation speed, rad/sec
0.00	0.054	15.8825
6.04	0.24	13.1598
9.74	0.36	11.5017
15.28	0.54	9.0583
19.01	0.66	7.1035

The graphical representation of measured results:

Motor parameters estimation

Now we can find the motor parameters. We will use the first and the last row of the table with
measurement results. So:
$$\left\{\begin{aligned} \tau_{d_1}&=0 ~ N\cdot m\\ I_1&=0.054 ~ A\\ \omega_1&=15.8825 ~ rad/sec \end{aligned}\right.$$
$$\left\{\begin{aligned} \tau_{d_2}&=0.1901 ~ N\cdot m\\ I_2&=0.66 ~ A\\ \omega_2&=7.1035 ~ rad/sec \end{aligned}\right.$$
Using these reference data we can find the solution of the system $\eqref{eq:eq_load2}$:
$$\left\{\begin{aligned} B&=-0.06000677881 \cdot A_r+0.001016586247 \\ K_b &= 0.4716532815\\ R_a &= 6.832750917\\ K_\tau &= 0.8693067325 \cdot A_r+0.2989986520 \end{aligned}\right.$$
Using the expression for B and the value for T2 found during the deceleration curve measurement $\eqref{eq:eq_decel1}$ we can find $A_r$:
$$-0.06000677881 \cdot A_r+0.001016586247=\frac{A_r}{10.697523425732065}$$
From this equation we can find $A_r=0.006623300293$. Then we can calculate other parameters:
$$\begin{equation} \begin{aligned} B&=0.0006191433314 \\ K_b &= 0.4716532815\\ R_a &= 6.832750917\\ K_\tau &= 0.3047563315 \end{aligned} \label{eq:eq_step2_params} \end{equation}$$

Because this method has relatively low accuracy (because I used a simple multimeter), I used acceleration curve measurement to find motor parameters with better precision.

3. Acceleration curve

In this case we will find the solution of the motor equation system $\eqref{eq:eq1}$ for $\omega(t)$. It has following form:
$$\begin{equation} \omega(t)=C_1 e^{\textstyle (\varkappa_1+\varkappa_3)t} + C_2 e^{\textstyle (\varkappa_1-\varkappa_3)t}+\frac{K_\tau U-R_a (A_r+\tau_d)}{C_5} \label{eq:eq_omega} \end{equation}$$
where
$$C_1=\frac{C_4}{C_5}(\frac{\varkappa_2}{\varkappa_3} - 1)  (\varkappa_1 - \varkappa_3)\\ C_2= - \frac{C_3}{C_5}(\frac{\varkappa_2}{\varkappa_3} + 1) (\varkappa_1 + \varkappa_3)\\ C_3=\frac{L_a}{2}(A_r + \tau_d + \frac{J U}{K_b}(\varkappa_2 - \varkappa_3))\\ C_4=\frac{L_a}{2}(A_r + \tau_d+ \frac{J U}{K_b}(\varkappa_2 + \varkappa_3))\\ C_5=B R_a + K_b K_\tau$$
and where
$$\varkappa_1 = -\frac{B L_a + J R_a} {2 J L_a}\\ \varkappa_2 = \frac{B L_a - J R_a} {2 J L_a}\\ \varkappa_3 = \sqrt{\varkappa_2^2 - \frac{K_b K_\tau}{J L_a}}$$

Because we can get only motor rotation angle from the motor's encoder, we should find expression for $\phi(t)$ by integrating $\eqref{eq:eq_omega}$ with respect of initial condition  $\phi(0) = 0$:
$$\phi(t)=C_6 e^{\textstyle (\varkappa_1+\varkappa_3)t} + C_7 e^{\textstyle (\varkappa_1-\varkappa_3)t}+\frac{(K_\tau U-R_a (A_r+\tau_d))t}{C_5}-(C_6+C_7)$$
where
$$C_6=\frac{C_1}{\varkappa_1+\varkappa_3}\\ C_7=\frac{C_2}{\varkappa_1-\varkappa_3}$$
I used the SciLab script to calculate motor parameters. It uses values $\eqref{eq:eq_step2_params}$ as initial values and uses $\eqref{eq:eq_decel1}$ to calculate $J$ and $A_r$. The source data is the same as for deceleration curve measurement. The script produces following motor parameters:
$$\begin{equation} \begin{aligned} R_a&=6.832749059810827\\ B&=0.000726962269165\\ K_\tau&=0.304766706036738\\ K_b&=0.459965726538748\\ J&=0.001502739083882\\ A_r&=0.007776695904018\\ L_a&=0.00494 \end{aligned} \label{eq:eq_step3_params} \end{equation}$$
The comparison of experimental and calculated motor positions:

Experimental and calculated motor angle

Conclusion

The mathematical model of Lego EV3 motor allows us to calculate controllers using modern techniques like this and check them in a virtual environment.