2. Hopping robot dynamics

This section discusses the dynamic model design of SLHR using the EL equation. For balancing the robot, an arm is pivoted to the central pivot and coerces the pitching motion of the robot body. For many-legged machines, having dynamic characteristics be captured by the robot and do not affect by the constraint of the pitching motion. The degree of freedom (DOF) of the robot leg is two. So because of one, the robot is rotating about the axis which is parallel to the support arm, that DOF is called a revolute DOF and because of another one, the robot leg compresses or extends, that DOF is called a prismatic DOF [Reference Lynch and Park21]. The revolute DOF is actuated or active while the prismatic DOF is passive and contains a linear passive spring.

Because we want to study the robot’s planar motion that is why we use 2D robot model dynamics which is shown in Fig. 1. For the same model, we have obtained the following [Reference Cherouvim and Papadopoulos22, Reference Yu, Fu and Chen23].

m = Mass due to the robot mass and partly because of the arm of support

L = Leg’s length when it is not compressed

k = Stiffness of the linear spring

$\tau$ = Torque applied by the hip actuator

r = Length of the leg at any time

b = Co-efficient of viscous friction with which energy loss in the leg prismatic DOF is modeled

Effects due to the impacts with the ground are included in the lumped parameters k, b.

For the phase stance, by the Lagrangian approach [Reference Morin24], we find the dynamics with Cartesian coordinate x and y. The body inertial force and spring force are much more than the force which is by the leg mass. So the effect of leg mass is neglected. Because of complexity, we find the dynamics in the function of r and $\phi$ . We also neglect the effect of gravity because the leg is already of small length and during stance, it will be further reduced. Now we apply the Lagrangian approach, consider Fig. 1.

The total kinetic energy (T)

\begin{equation}\nonumber T={{\frac{1}{2}}m\dot{x}^2}+{{\frac{1}{2}}m\dot{y}^2}\end{equation}

And the total potential energy is

\begin{equation*}V={{\frac{1}{2}}k{x}^{2}}+{{\frac{1}{2}}k{y}^{2}}\end{equation*}

So the Lagrangian ( $\mathcal{L}_a$ ) is

\begin{equation*}{\mathcal{L}}_{a}={T-V}\end{equation*}

\begin{equation}\mathcal{L}_a=\left({{\frac{1}{2}}m\dot{x}^{2}}+{{\frac{1}{2}}m\dot{y}^{2}}\right)-\left({{\frac{1}{2}}k{x}^{2}}+{{\frac{1}{2}}k{y}^{2}}\right)\nonumber\end{equation}

The Lagrangian equation for x coordinate of form is

\begin{equation}\frac{d}{dt}\left(\frac{\mathcal{L}_a}{\dot{x}}\right)-\frac{\mathcal{L}_a}{x}=F_x\nonumber\end{equation}

The Lagrangian equation for y coordinate of form is

\begin{equation}\frac{d}{dt}\left(\frac{\mathcal{L}_a}{\dot{y}}\right)-\frac{\mathcal{L}_a}{y}=F_y\nonumber\end{equation}

where $F_x$ and $F_y$ are the non-conservative components of the force by viscous friction and actuator torque. The Lagrangian equation for x and y direction is

(1) \begin{equation}m\ddot{x}+k(L-r)sin(\phi) -b\dot{r}\;sin(\phi)=-\frac{\tau}{r}cos(\phi)\end{equation}

(2) \begin{equation}m\ddot{y}-k(L-r)cos(\phi) +b\dot{r}\;cos(\phi)=-\frac{\tau}{r}sin(\phi)\end{equation}

For the flight phase, only the gravitational force is acting so the speed in the horizontal direction is constant and because of gravitation, a constant acceleration is acted in the vertical direction.

(3) \begin{equation}\ddot{x}=0\end{equation}

(4) \begin{equation}\ddot{y}=-g\end{equation}

3. The apex height

As per the [Reference Cherouvim and Papadopoulos17] the apex height is the height of the robot at the topmost point of its body when it is in any of the two phase (namely, stance and flight). Hence, the correct measure of the apex height is very important when it comes to overall balanced performance of the hopping robot. If the terrain/surface on which the robot is moving is uneven, then also it should be able to figure out the necessary distance from the surface to stay in stable operating mode. The fixed distance from the ground to the top point of the robot during its operation is the apex height.

Now we will compute the rate of changes $\dot{\phi}$ and $\dot{r}$ which is used for the determination of vertical dynamics. From Fig. 1 of the stance phase, we get

(5) \begin{equation} x-x_{fp}=-r\;sin{\phi}\end{equation}

(6) \begin{equation} y=rcos{\phi}\end{equation}

where $x_{fp}$ is the position during stance of the foot on the ground.

By differentiating Eqs. (5) and (6), we get $\dot{x}$ and $\dot{y}$ of the robot body.

(7) \begin{equation} \dot{x}=-\dot{r} sin{\phi}-r\dot{\phi}cos{\phi}\end{equation}

(8) \begin{equation} \dot{y}=\dot{r}cos{\phi}-r\dot{\phi}sin{\phi}\end{equation}

For small leg angle, Eq. (8) becomes

(9) \begin{equation} \dot{y}=\dot{r}cos{\phi}\end{equation}

We initiate from Eq. (2) of dynamics, which is about the robot’s vertical motion. Substituting the leg’s length r and the leg’s angle $\phi$ as expressions of the coordinates x, y of the robot body, and by the use of small-angle approximations of trigonometry, the vertical dynamics becomes:

(10) \begin{equation}m\ddot{y}+b\dot{y}+ky=kLcos{\phi}\end{equation}

4. State observer design

Practically for feedback some of the state variables are not available. This is why the estimation of the unavailable state variables becomes important. For the unmeasurable state variable estimation is done. The process of estimation is called observation.

State Observer: State variables are observed or estimated by a device called a state observer or simply an observer [Reference Chen25, Reference Ogata26]. The block diagram of state observer is given in Fig. 2.

• Based upon the output and control variables measurement, the observation of state variables is done by state observer

• Only after satisfaction and fulfillment of the observability condition design of state observer is possible

The vertical dynamics is

\begin{equation}m\ddot{y}+b \dot{y}+k y = k L\cos{\phi}\nonumber\end{equation}

Now,

(11) \begin{equation} \ddot{y} =-\frac{b}{m}\dot{y}-\frac{k}{m}y+\frac{k L}{m}\cos{\phi} \end{equation}

Figure 2. Block diagram of state observer (uncolored portion is the actual system and colored portion is the estimated system).

Let

\begin{align*} x_1 & = y\\ \dot{x_1} & = x_{2} = \dot{y}\\ \dot{x_2} &= \ddot{y} \end{align*}

So,

(12) \begin{equation} \dot{x_1}=x_2 \end{equation}

and

\begin{equation} \dot{x_2}=-\frac{b}{m}x_2- \frac{k}{m}x_1+\frac{k L}{m}\cos{\phi} \nonumber \end{equation}

Let

\begin{equation*}u=\cos{\phi}\end{equation*}

Therefore,

(13) \begin{equation} \dot{x_2}=-\frac{b}{m}x_2- \frac{k}{m}x_1+\frac{k L}{m}u \end{equation}

Using Eqs. (12) and (13), the state equation $\dot{x}=Ax+Bu $ is

(14) \begin{equation} \begin{bmatrix}\dot{x_1} \\ \dot{x_2} \end{bmatrix}= \begin{bmatrix} 0\;\;\;\; & 1 \\-\frac{k}{m}\;\;\;\; & -\frac{b}{m} \end{bmatrix} \begin{bmatrix} x_1\\x_2 \end{bmatrix} +\begin{bmatrix} 0 \\ \frac{k L}{m} \end{bmatrix} u\end{equation}

where $x_1,x_2 $ are the state vector, and u is the input vector.

Hence,

\begin{align*}A=\begin{bmatrix} 0\;\;\;\; & 1 \\-\frac{k}{m}\;\;\;\; & -\frac{b}{m} \end{bmatrix} \quad and \quad B=\begin{bmatrix} 0 \\ \frac{k L}{m} \end{bmatrix} \end{align*}

Since

\begin{equation*}y=x_1\end{equation*}

Therefore, the output equation is

(15) \begin{equation}y=\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix} \begin{bmatrix} x_1\\x_2 \end{bmatrix}\end{equation}

Hence,

\begin{align*} C=\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix} \quad and \quad D=0\end{align*}

Let us first check the observability condition by computing the observability matrix.

(16) \begin{equation}\mathcal{O}=[C^*\;:\;A^*C^*]=\begin{bmatrix} 0\;\;\;\; & 1\\1\;\;\;\; & 0\end{bmatrix}\end{equation}

It can be observed that the observability matrix $\mathcal{O}$ has a rank 2. Since the system is completely observable therefore it is possible to determine the required observer gain matrix.

Let us assume

\begin{align*}u=-k\tilde{x}\end{align*}

For the system to be stable, we assume the eigenvalues of the observer matrix are

\begin{align*} \mu_{1} &=-5 \\ \mu_{2} & = -5 \end{align*}

Hence, the characteristics equation for assumed eigenvalues is

(17) \begin{equation} s^2+10s+25=0\end{equation}

The model and the values of L, m and b are inspired by the article of Cherouvim [Reference Cherouvim and Papadopoulos17]. Here the variation in the spring constant for three different values like $k_1,k_2$ , and $k_3$ has been done. The changes in the apex height and vertical velocity are presented and shown in the section of simulation results. The models and analysis are as follows:

\begin{align*}L=0.275\;{\rm m}, m = 4\;{\rm kg},\; and\quad b=6\Big(\frac{Ns}{m} \Big)\end{align*}

For ${{k}}_{{1}}\;\textbf{=}\;\textbf{4000}$ $\Big(\frac{N}{m}\Big)$

So

\begin{align*}\frac{k_1}{m}=1000,\frac{b}{m}=1.5\quad and\quad \frac{k_1 L}{m}=275\end{align*}

Therefore for the numeric value the state equation is

\begin{align*} \begin{bmatrix}\dot{x_1} \\ \dot{x_2} \end{bmatrix}&= \begin{bmatrix} 0\;\;\;\; & 1 \\-1000\;\;\;\; & -1.5 \end{bmatrix} \begin{bmatrix} x_1\\x_2 \end{bmatrix} +\begin{bmatrix} 0 \\ 275 \end{bmatrix} u \end{align*}

Now we write $\dot{x_1}$ and $\dot{x_2}$ separately

(18) \begin{align} \dot{x_1}=x_2 \end{align}

(19) \begin{align} \dot{x_2} = -1000{x_1}-1.5{x_2}+275{u} \end{align}

For observer design of given system:

\begin{align*}A_1=\begin{bmatrix} 0\;\;\;\; & 1 \\-1000\;\;\;\; & -1.5 \end{bmatrix},B_1=\begin{bmatrix} 0 \\ 275 \end{bmatrix} \text{and C be the same.}\end{align*}

The state observer design is reduced to the evaluation of appropriate observer gain matrix $k_e$ . The observer error equation is

\begin{equation*} \dot{e}=(A_1-K_eC)e\end{equation*}

The characteristics equation of the observer becomes

\begin{equation*}|sI-A_1+K_eC|=0\end{equation*}

where

\begin{equation*}k_e=\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\end{equation*}

Then the characteristics equation becomes

\begin{equation*}\Bigg| \begin{bmatrix} s\;\;\;\; & 0\\0\;\;\;\; & s\end{bmatrix}-\begin{bmatrix} 0\;\;\;\; & 1 \\-1000\;\;\;\; & -1.5 \end{bmatrix}+\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg |=0\\\end{equation*}

(20) \begin{equation} s^2+(K_{e1}+1.5)s+1.5K_{e1}+K_{e2}+1000=0 \end{equation}

By comparing Eqs. (20) and (17), we get the value of $k_e$ .

\begin{align*} K_{e1} &=8.5\\ K_{e2} &=-987.75\end{align*}

So the equation of full order state observer is

\begin{equation*} \dot{\tilde{x}}=(A_1-K_e{C})\tilde{x}+B_1u+K_e{y}\end{equation*}

By putting $y=C{x}$ , we get

(21) \begin{equation} \dot{\tilde{x}}=(A_1-K_e{C})\tilde{x}+B_1u+K_e{C{x}}\end{equation}

\begin{align*} \begin{bmatrix} \dot{\tilde{x}}_{1} \\ \dot{\tilde{x}}_{2} \end{bmatrix}=\Bigg( \begin{bmatrix} 0\;\;\;\; & 1 \\-1000\;\;\;\; & -1.5 \end{bmatrix}-\begin{bmatrix} 8.5 \\ -987.75\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg )\begin{bmatrix} \tilde{x_1} \\ \tilde{x_2} \end{bmatrix}+\begin{bmatrix} 0 \\ 275 \end{bmatrix}u+\Bigg [ \begin{bmatrix} 8.5 \\ -987.75\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg ]\begin{bmatrix} {x_1} \\ {x_2} \end{bmatrix}\end{align*}

Therefore, the $\dot{\tilde{x}}_{1}$ and $\dot{\tilde{x}}_{2}$ are

(22) \begin{align} \dot{\tilde{x}}_{1} =-8.5\tilde{x_1}+\tilde{x_2}+8.5x_1 \end{align}

(23) \begin{align} \dot{\tilde{x}}_{2} =-12.25\tilde{x_1}-1.5\tilde{x_2}-987.75x_1+275{u} \end{align}

For $k_{2}\,\textbf{=}\,\textbf{4800}$ $\Big(\frac{N}{m}\Big)$

So

\begin{align*}\frac{k_2}{m}=1200,\frac{b}{m}=1.5\quad and\quad \frac{k_2 L}{m}=330\end{align*}

Therefore for the numeric value the state equation is

\begin{align*} \begin{bmatrix}\dot{x_1} \\ \dot{x_2} \end{bmatrix}&= \begin{bmatrix} 0\;\;\;\; & 1 \\-1200\;\;\;\; & -1.5 \end{bmatrix} \begin{bmatrix} x_1\\x_2 \end{bmatrix} +\begin{bmatrix} 0 \\ 330 \end{bmatrix} u \end{align*}

Now we write $\dot{x_1}$ and $\dot{x_2}$ separately

(24) \begin{align} \dot{x_1} = x_2 \end{align}

(25) \begin{align} \dot{x_2} = -1200{x_1}-1.5{x_2}+330{u} \end{align}

For observer design of given system:

\begin{align*}A_2=\begin{bmatrix} 0\;\;\;\; & 1 \\-1200\;\;\;\; & -1.5 \end{bmatrix},B_2=\begin{bmatrix} 0 \\ 330 \end{bmatrix} \text{and C be the same.}\end{align*}

The observer error equation is

\begin{equation*} \dot{e}=(A_2-K_eC)e\end{equation*}

The characteristics equation of the observer becomes

\begin{equation*}|sI-A_2+K_eC|=0\end{equation*}

where

\begin{equation*}k_e=\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\end{equation*}

Then the characteristics equation becomes

\begin{equation*}\Bigg| \begin{bmatrix} s\;\;\;\; & 0\\0\;\;\;\; & s\end{bmatrix}-\begin{bmatrix} 0\;\;\;\; & 1 \\-1200\;\;\;\; & -1.5 \end{bmatrix}+\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg |=0\\\end{equation*}

(26) \begin{equation} s^2+(K_{e1}+1.5)s+1.5K_{e1}+K_{e2}+1200=0 \end{equation}

By comparing Eqs. (26) and (17), we get the value of $k_e$ .

\begin{align*} K_{e1} &=8.5\\ K_{e2} &=-1187.75\end{align*}

So the equation of full order state observer is

\begin{equation*} \dot{\tilde{x}}=(A_2-K_e{C})\tilde{x}+B_2u+K_e{y} \end{equation*}

By putting $y=C{x}$ , we get

(27) \begin{equation} \dot{\tilde{x}}=(A_2-K_e{C})\tilde{x}+B_2u+K_e{C{x}}\end{equation}

\begin{align*} \begin{bmatrix} \dot{\tilde{x}}_{1} \\ \dot{\tilde{x}}_{2} \end{bmatrix}=\Bigg( \begin{bmatrix} 0\;\;\;\; & 1 \\-1200\;\;\;\; & -1.5 \end{bmatrix}-\begin{bmatrix} 8.5 \\ -1187.75\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg )\begin{bmatrix} \tilde{x_1} \\ \tilde{x_2} \end{bmatrix}+\begin{bmatrix} 0 \\ 330 \end{bmatrix}u+\Bigg [ \begin{bmatrix} 8.5 \\ -1187.75 \end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg ]\begin{bmatrix} {x_1} \\ {x_2} \end{bmatrix}\end{align*}

Therefore, the $\dot{\tilde{x}}_{1}$ and $\dot{\tilde{x}}_{2}$ are

(28) \begin{align} \dot{\tilde{x}}_{1} =-8.5\tilde{x_1}+\tilde{x_2}+8.5x_{1} \end{align}

(29) \begin{align} \dot{\tilde{x}}_{2} =-12.25\tilde{x_1}-1.5\tilde{x_2}-1187.75x_1+330{u} \end{align}

For $k_3\;=\;5600$ $\Big(\frac{N}{m}\Big)$

So

\begin{align*}\frac{k_3}{m}=1400,\frac{b}{m}=1.5\quad and\quad \frac{k_3 L}{m}=385\end{align*}

Therefore for the numeric value the state equation is

\begin{align*} \begin{bmatrix}\dot{x_1} \\ \dot{x_2} \end{bmatrix}&= \begin{bmatrix} 0\;\;\;\; & 1 \\-1400\;\;\;\; & -1.5 \end{bmatrix} \begin{bmatrix} x_1\\x_2 \end{bmatrix} +\begin{bmatrix} 0 \\ 385 \end{bmatrix} u \end{align*}

Now we write $\dot{x_1}$ and $\dot{x_2}$ separately

(30) \begin{align} \dot{x_1}=x_2 \end{align}

(31) \begin{align} \dot{x_2}=-1400{x_1}-1.5{x_2}+385{u} \end{align}

For observer design of given system:

\begin{align*}A_3=\begin{bmatrix} 0\;\;\;\; & 1 \\-1400\;\;\;\; & -1.5 \end{bmatrix},B_3=\begin{bmatrix} 0 \\ 385 \end{bmatrix} \text{and C be the same.}\end{align*}

The observer error equation is

\begin{equation*} \dot{e}=(A_3-K_eC)e\end{equation*}

The characteristics equation of the observer becomes

\begin{equation*}|sI-A_3+K_eC|=0\end{equation*}

where

\begin{equation*}k_e=\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\end{equation*}

Then the characteristics equation becomes

\begin{equation*}\Bigg| \begin{bmatrix} s\;\;\;\; & 0\\0\;\;\;\; & s\end{bmatrix}-\begin{bmatrix} 0\;\;\;\; & 1 \\-1400\;\;\;\; & -1.5 \end{bmatrix}+\begin{bmatrix} K_{e1} \\ K_{e2}\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg |=0\\\end{equation*}

(32) \begin{equation} s^2+(K_{e1}+1.5)s+1.5K_{e1}+K_{e2}+1400=0 \end{equation}

By comparing Eqs. (32) and (17), we get the value of $k_e$ .

\begin{align*} K_{e1} &=8.5\\ K_{e2} &=-1387.75\end{align*}

So the equation of full order state observer is

\begin{equation*} \dot{\tilde{x}}=(A_3-K_e{C})\tilde{x}+B_3u+K_e{y} \end{equation*}

By putting $y=C{x}$ , we get

(33) \begin{equation} \dot{\tilde{x}}=(A_3-K_e{C})\tilde{x}+B_3u+K_e{C{x}}\end{equation}

\begin{align*} \begin{bmatrix} \dot{\tilde{x}}_{1} \\ \dot{\tilde{x}}_{2} \end{bmatrix}=\Bigg( \begin{bmatrix} 0\;\;\;\; & 1 \\-1400\;\;\;\; & -1.5 \end{bmatrix}-\begin{bmatrix} 8.5 \\ -1387.75\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg )\begin{bmatrix} \tilde{x}_{1} \\ \tilde{x}_{2} \end{bmatrix}+\begin{bmatrix} 0 \\ 385 \end{bmatrix}u+\Bigg [ \begin{bmatrix} 8.5 \\ -1387.75\end{bmatrix}\begin{bmatrix} 1\;\;\;\; & 0 \end{bmatrix}\Bigg ]\begin{bmatrix} {x_1} \\ {x_2} \end{bmatrix}\end{align*}

Therefore, the $\dot{\tilde{x}}_{1}$ and $\dot{\tilde{x}}_{2}$ are

(34) \begin{align} \dot{\tilde{x}}_{1} =-8.5\tilde{x_1}+\tilde{x_2}+8.5x_{1} \end{align}

(35) \begin{align} \dot{\tilde{x}}_{2} =-12.25\tilde{x_1}-1.5\tilde{x_2}-1387.75x_1+385{u} \end{align}

This is the required estimated dynamics with $k_3=5600$ . In the next section, the simulation results for all the above models have been presented.

5. Simulation results and discussion

In this section, we discuss the simulation of observer design of the apex height and vertical velocity of the SLHR. For the same, we use the equation of vertical dynamics during the stance phase of SLHR using Matlab. After simulation of the vertical height and the vertical velocity of the actual system, we can compare it with the estimated system. Here we have obtained six-state trajectories which are further segregated into three parts for three different spring constants. Each part is having vertical velocity and an apex height of the actual and estimated system. As it is observed in Figs. 3(a), 4(a), and 5(a), the actual height is the same as the estimated one; this implies that our estimation is correct. It can be seen in Figs. 3(a), 4(a), and 5(a) that the robot is launched from the initial height $0.5$ m. It performs hopping and a few seconds later the robot acquires its original height of $0.275$ m because of damping. The robot apex height during hopping is large for large spring constant which we can be observed in Figs. 3(a), 4(a), and 5(a).

Figure 3. Estimated and actual height and velocity of SLHR for spring constant $k_1=4000$ . (a) Estimated and actual height and (b) estimated and actual velocity.

Figure 4. Estimated and actual height and velocity of SLHR for spring constant 4800. (a) Estimated and actual height and (b) estimated and actual velocity.

Figure 5. Estimated and actual height and velocity of SLHR for spring constant 5600. (a) Estimated and actual height and (b) estimated and actual velocity.

The vertical velocity of SLHR is following the actual system’s vertical velocity; this also implies that our estimation is correct and after 7–8 s, the vertical velocity approaches zero that means the robot is attaining rest position due to damping. Velocity is positive and negative because of the direction of movement of the robot either upward or downward. For the value of spring constant as 4000 the velocity is $1.99$ m/s, which can be noticed in the first spike of the oscillatory trajectory shown in Fig. 3(b). Whereas vertical velocity for spring constant 4800 and 4800 stays close to 2 m/s in the first spike and then decreases with time shown in Figs. 4(b) and 5(b). For vertical velocity, the trajectory is almost the same for all the three spring constants. There are a few limitations of the present contribution is that the span of choices of spring constant values is limited and due to physical constraints of the observer design, it is difficult for the case of flight phase to make an optimal choice of spring constant values.