CPG Human Motion Phase Recognition Algorithm for a Hip Exoskeleton with VSA Actuator

Jiaxuan Li¹, Feng Jiang², Longhai Zhang¹, Xun Wang¹, Jinnan Duan¹, Baichun Wei¹, Xiulai Wang^3*, Ningling Ma³, Yutao Zhang³

¹School of Mechatronics Engineering, Harbin Institute of Technology, Harbin, China.
²School of Future Technology, Nanjing University of Information Science & Technology, Nanjing, China.
³Nanjing Jinling Hospital, Affiliated Hospital of Medical School, Nanjing University, Nanjing, China.

DOI: 10.4236/jsip.2024.152002                

Abstract

Due to the dynamic stiffness characteristics of human joints, it is easy to cause impact and disturbance on normal movements during exoskeleton assistance. This not only brings strict requirements for exoskeleton control design, but also makes it difficult to improve assistive level. The Variable Stiffness Actuator (VSA), as a physical variable stiffness mechanism, has the characteristics of dynamic stiffness adjustment and high stiffness control bandwidth, which is in line with the stiffness matching experiment. However, there are still few works exploring the assistive human stiffness matching experiment based on VSA. Therefore, this paper designs a hip exoskeleton based on VSA actuator and studies CPG human motion phase recognition algorithm. Firstly, this paper puts forward the requirements of variable stiffness experimental design and the output torque and variable stiffness dynamic response standards based on human lower limb motion parameters. Plate springs are used as elastic elements to establish the mechanical principle of variable stiffness, and a small variable stiffness actuator is designed based on the plate spring. Then the corresponding theoretical dynamic model is established and analyzed. Starting from the CPG phase recognition algorithm, this paper uses perturbation theory to expand the first-order CPG unit, obtains the phase convergence equation and verifies the phase convergence when using hip joint angle as the input signal with the same frequency, and then expands the second-order CPG unit under the premise of circular limit cycle and analyzes the frequency convergence criterion. Afterwards, this paper extracts the plate spring modal from Abaqus and generates the neutral file of the flexible body model to import into Adams, and conducts torque-stiffness one-way loading and reciprocating loading experiments on the variable stiffness mechanism. After that, Simulink is used to verify the validity of the criterion. Finally, based on the above criterions, the signal mean value is removed using feedback structure to complete the phase recognition algorithm for the human hip joint angle signal, and the convergence is verified using actual human walking data on flat ground.

Keywords

Variable Stiffness Actuator, Plate Spring, CPG Algorithm Convergence Criterion, Human Motion Phase Recognition, Simulink and Adams Co-Simulation

1. Introduction

The research of human exoskeleton and prosthetics is developing rapidly these years. The number of annual published articles in the field has tripled compared to 10 years ago [1]. One of the research hotspots is how to ensure safe, smooth, and effective human-machine interaction in exoskeletons. Classical methods such as the impedance/admittance assistive control algorithm are widely used in assisting the human body [2] [3], and their mechanical structure is mostly driven by direct drive or in-series spring structure design, which rely entirely on the end torque sensor to measure and filter for control. It often shows strong hysteresis to the human body structure, and insufficient bandwidth for disturbance response to human body motion dynamics, which cannot provide an effective experimental platform for further research on stiffness matching issues. Advanced controllers are often identified as showing considerable heterogeneity, each one suited for a specific kind of application, target population and performance target, being not universal enough and lacking the value of comparative study [4].

In recent years, there have been more and more design cases of variable stiffness actuators (VSA), all of which can achieve physical variable stiffness or impedance. Compared with other driving structures, the advantages of VSA structures are their robustness [5], energy efficiency [6], and adaptability to complex tasks [7]. Evidence of these advantages is human joint, which is a typical representative of antagonistic variable stiffness actuators, works far better than current machine joint drives in these properties. According to their working principles, variable stiffness structures can be mainly divided into three types: antagonistic, variable transmission ratio, and material physical properties. One of the typical uses of antagonistic structures is taken out by the Georgia Institute of Technology [8], which uses a nonlinear spring with a force-deformation curve that exhibits a parabolic shape, built with a linear spring and cam structure, and is driven antagonistically with two stepper motors. The typical structures of variable transmission ratio type are Aw AS-I and Aw AS-II developed by the Italian Institute of Technology [9] [10], which uses a torsion spring as an elastic element, and controls stiffness by controlling the force arm r of the torsion spring. This type of structure also has typical representatives such as VS-Joint, QA-Joint, and FSJ developed by the German Aerospace Center [11]-[13]. The advantage of this type of structure is that the output and stiffness control are usually driven by motors separated from each other, with high energy efficiency and a large stiffness adjustment range. However, the volume of the structure is often relatively large, the quality is also relatively high, the internal inertia of the motor is high when outputting, and it is difficult to achieve an arbitrary rotation angle as the motor. Material physical property type structures use the self-physical properties of the elastic element, represented by the Jack spring [14]. The stiffness is changed by controlling the working length of the spring by a rotating central axis extension and contraction of the spring. There are also typical structures that change the external stiffness by changing the working length of the plate spring, such as the arm structure designed by Waseda University, which uses the sliding block to change the working range of the plate spring [13]. There are also methods that use a rotatable plate spring inside a spring to achieve variable stiffness [15] [16], and so on.

Although there are many technical routes and corresponding achievements in structural design, the upper-level control required for VSA still has issues with real-time performance and stability that cannot be achieved. Currently, the upper limb control algorithm methods for calculation principles and starting points are mainly divided into four categories: based on human body models, based on machine learning calculations, based on sensor information feature points, and based on oscillator calculations [17]. The typical control algorithm based on human body models, such as that developed by the Harbin Institute of Technology, obtains the pressure at the bottom of the foot through an optical motion capture system and force platform, and calculates the joint torque of the human lower limb through rigid-body inverse dynamics [18]. Based on this, David G. Lloyd and others from the University of Western Australia changed the parameters of the muscle model by changing the characteristic values of the electromyographic signals to simulate human muscle force, and then solved the joint torque of the human body [19]. Angelos Karatsidis and his colleagues from Delft University of Technology in the Netherlands measured human motion using an inertial measurement unit (IMU) and assisted the optical motion capture system with a force platform to calculate the joint torque of the human body [20]. This method is currently the most widely recognized approach, but it has higher requirements for equipment and sites, expensive equipment costs, and can only be used indoors. The control algorithm based on neural network calculation first began in 1993, when FRANCIS and others from the University of Catholic Chile used electromyography (EMG) to estimate joint angles and joint torque of humans and machines [21], with accuracy up to 73%. In 2019, Fuzhou University only used 10 EMG signals in the lower body and joint signal sources to identify the joint torque of the lower limb of the human body with an accuracy rate of P > 0.9633 [22], and recent research is becoming more and more mature [23]-[25]. However, due to the limitations of EMG sensor fit and neural network transferability, the current practical application scenarios are still relatively narrow. The control algorithm based on sensor information feature points has various types, such as in 2015, when Samsung Industrial University used the periodicity of the vertical acceleration of the human torso during walking to identify human motion phases [26]. The NREL-Exo assistive exoskeleton developed by North Carolina State University determines the walking state of the human body by detecting joint angles and angular velocities [27]. The VariLeg exoskeleton assistive device developed by the Swiss Rehabilitation Engineering Laboratory obtains the current phase of the human body through feedback from the foot pressure sensor and joint angle and angular velocity [28]. In 2016, Harvard University designed a hip exoskeleton that identifies human motion phases using the maximum joint extension [29], and an ankle exoskeleton that uses ankle joint angular velocity zero points as phase identification feature points [30], respectively. L. Bergmann et al. provide an unscented Kalman filter (UKF) to estimate the joint torques the subject gives in to the exoskeleton with VSA [31]. It can be applied over the entire stiffness range of the VSAs and performs better than the common inverse dynamics approach, while requiring reidentification on each subject for accurate fitting. The oscillator-based method is a phase recognition calculation method that uses CPG learning function (rather than the form of a rhythmic signal generated by CPG unit coupling). As a newly emerged control method, although its stability has not yet been determined, it has advantages of simple and accurate use in the recognition of human body movements with changing gait. Compared with methods based on sensor information feature points, CPG is widely present in vertebrates to control periodic movements, so external CPG and human internal CPG can construct influencing paths to make phase recognition reasonably possible. Samsung Polytechnic College designed the world’s first exoskeleton based on CPG assistance [32], which uses an adaptive oscillator to learn human joint movements, estimate the phase with no delay, and finally provide percentage assistance based on the statistical torque curve. Apart from Samsung, Arizona State University has used the Combination of Central Pattern Generator and Human Intent calculation to design a hip exoskeleton assistive device [33]. Currently, there is limited literature on the accuracy of CPG phase calculation, but from the results presented in the current literature [34] [35], it has some practical value, but further exploration is still needed. Also. there is currently no research on using VSA to explore assistive human stiffness matching experiments. So, this paper will start with the human lower limbs that have been more researched, design a hip exoskeleton based on VSA, and study the CPG human motion phase recognition algorithm.

The rest parts of this paper are arranged as follows: In the second section, this paper establishes the performance indicators of the variable stiffness human-assisted experimental apparatus, and then uses two motors to drive a four-linkage mechanism together, and uses a plate spring as the elastic element to establish the principle of variable stiffness and complete the main parameter design. According to the main structure parameters, the mechanism detailed design is completed, the dynamic model of the VSA mechanism is deduced and analyzed, and the nonlinear problems in the variable stiffness mechanism are identified, and the upper-level control structure is divided. In the third section, this paper uses perturbation theory to expand the first-order CPG unit and obtains the main equation for phase convergence under the same frequency and waveform. Then, the convergence of the hip joint angle signal during walking on flat ground is deduced. After determining the main equation, this paper establishes the main equation for CPG convergence of external periodic signal frequency under the premise of circular limit cycle, analyzes the convergence criterion of frequency under circular limit cycle, and identifies the factors affecting the frequency convergence of hip joint angle. Furthermore, this paper extends these criteria to the convergence of CPG units under a complex limit cycle. In the fourth section, this paper uses Abaqus to build a plate spring elastic model, conducts one-way and reciprocating loading experiments on the variable stiffness mechanism in the joint simulation of Simulink and Adams, analyzes the experimental data, and uses the above-derived criteria to establish the CPG phase recognition algorithm, remove the signal mean value through feedback structure and verify it using actual hip joint angle signals during walking.

2. Design and Modeling of Variable Stiffness Actuator

2.1. Definition and Design Requirements of Stiffness

The design of variable stiffness actuators needs to consider the form and parameter requirements of variable stiffness firstly. In the field of human-machine contact, stiffness is divided into two types: intrinsic stiffness, which is manifested as having elastic elements internally and can exhibit displacement according to corresponding stiffness rules when subjected to external disturbances; and apparent stiffness, which depends on the control system to achieve a stiffness rule. When external force disturbances occur, the executing mechanism responds with a displacement according to the control rule, thereby having a certain virtual stiffness from the outside perspective. The typical representative of apparent stiffness is motor impedance control.

The biggest difference between these two types of stiffness is in their external response bandwidth. Intrinsic stiffness is a mechanical property of the structure’s interior, so the external bandwidth and response speed are much greater than the stiffness exhibited by the controlled action. Currently, apparent stiffness is often used to assist the human body in the human-machine coupling direction. Although this method means that the structure and control method are easy to implement, it often faces problems such as sensor information processing, actuator driving bandwidth, and control stability in the face of fast motion. It is difficult to keep stable and effective in various test environments. Therefore, this paper designs variable stiffness starting from actuators with intrinsic stiffness.

In addition, based on statistical parameters such as torque and angular velocity of human lower limb movement, the experiment is designed to be conducted at a walking speed of 5 km/h. At the same time, in combination with the current performance indicators in the field of lower limb exoskeleton design, the remaining performance indicator requirements are specified as follows:

1) The torsional stiffness ranges from 30 N·m/rad to 150 N·m/rad;

2) The change from maximum stiffness to minimum stiffness does not exceed 0.1 s;

3) The torque-bearing capacity at minimum stiffness is not less than 10 N·m;

4) The maximum output torque is not less than 25 N·m;

5) The minimum output angular velocity is not less than 10 rad/s.

6) The volume should be as small as possible, and the internal rotational inertia should be lower than 0.2 kg·m².

2.2. Structural Principles and Mechanism Design

Structural Principles

According to the design requirements, the entire actuator uses a long strip plate spring as the elastic element to control the change of the internal stiffness properties of the actuator by controlling the contact point position of the plate spring. The entire mechanism is divided into front-end mechanism and an elastic output mechanism. The front-end mechanism uses a pair of servo motors to drive a four-bar linkage mechanism together, making it rotate as a whole, controlled by the motor output common mode angle $\frac{{\varphi }_1+{\varphi }_2}{2}$ . The position of the plate spring’s contact point is controlled by the motor differential mode angle $\frac{{\varphi }_1-{\varphi }_2}{2}$ . The elastic output mechanism mainly fixes the plate spring and transmits torque output. The comprehensive design parameters of the mechanism are shown in Table 1.

Table 1. Design of basic mechanism parameters.

Name	Value
Overall shape of the plate spring (Lxbxh)	55 × 10 × 1 (mm)
Effective working part of the plate spring (Lxbxh)	35 × 10 × 1 (mm)
Material of the plate spring	60Si2CrVA, Young’s modulus E = 199e9 Pa, Poisson’s ratio v = 0.3, Yield strength 1.78e9 (Pa)
Outer radius R	55 (mm)
Stiffness range	30 - 216 (N∙m/rad)
Load torque at minimum stiffness	20 (N∙m)
Length of the four-bar linkage	18 mm

Take the differential mode angle of the motor $\frac{{\varphi }_1-{\varphi }_2}{2}={90}^{\circ }$ as the angle balance state, and its deviation angle is used as abscissa. The design deviation ranges from −0.2 rad to 0.8 rad, and the motor differential mode angle changes. The stiffness change curve of the variable stiffness mechanism is shown in Figure 1. The design stiffness ranges from a minimum of 30 N·m/rad to 216 N·m/rad, and the working point radius changes from 24 mm to 33 mm. The torque amplification gain changes from the minimum value of 2.8 at 0 rad to the maximum value of 3.7. It has been verified that the plate spring can withstand a negative load torque of 20 N·m/rad at the minimum stiffness.

Figure 1. Working curve design of variable stiffness actuator.

2.3. Structure Design

In structural design, the core idea is to split the four-bar linkage mechanism, so that the two bars on its both sides are installed on both sides of the plate spring, and a holding bracket is used to offset the torque generated by unbalanced installation, so that the mechanisms do not interfere with each other during operation and the structure is compact.

The core structural design is shown in Figure 2. Two servo motors are respectively connected to gear disk 1 and gear disk 2 with a 2:1 reduction ratio. Gear disk 1 and the driven disk are both mounted on the center axis along the Z direction through the key connection in order to drive the driven disk. In addition, gear disk 2 is also mounted on the center axis via angular contact ball bearings. The driven disk and gear disk 2 both have four mounting threaded holes with the same diameter and evenly distributed around the periphery with a diameter of 36 mm, and are connected to the corresponding swing arm on one end through pins installed in the threaded holes. The other end of the swing arm is connected to the corresponding slider through pins.

Considering that a set of four-bar linkage mechanisms is not in the plane, it produces a torsion tendency perpendicular to the Z direction on the slider during operation. Therefore, a holding bracket is installed on the slider, as shown in Figure 3. The holding bracket is mounted on the center axis via two non-standard angular contact bearings, allowing only rotational freedom around the O-Z direction. The holding bracket is fixed with 8 optical axes through threaded connection. The optical axes restrict the slider to only slide radially in the Z direction through linear bearings installed on the slider, ensuring that the four-bar linkage mechanism works normally.

There are two micro-needle roller bearings installed on the slider for contact with both sides of the plate spring to realize forward and reverse rotation. The plate spring installation and output are shown in Figure 2. To ensure that the plate spring only undergoes deformation by rotating around the Z-axis, the contact surface is made a uniform line contact. First, the plate spring is tightly mounted on the mounting seat with screws, then the mounting seat is fixed on the output bracket through ear holes with screws on both sides, and finally the output bracket outputs torque and displacement outward through the output shaft fixed on it.

Figure 2. Core structure diagram.

Figure 3. Retainer frame diagram.

The overall assembly 3D model of the entire mechanism is shown in Figure 4. The front and rear mounting plates and the output support shell are assembled and fixed with dowels to form the entire mechanism. There are angular contact ball bearings on the centers of the front and rear mounting plates for mounting the entire center shaft, and the majority of the components are restricted in the Z direction by the center shaft.

The overall 3D structure diagram is shown in Figure 5, where the entire mechanism is connected by fixed pins. Two servo motors are mounted on the rear mounting plate, and each motor is equipped with a standard carbon steel pinion with 16 teeth and a modulus of 2, which matches 2:1 reduction ratio to increase torque. An outer 16-bit photoelectric encoder is installed on the outer side of the rear mounting plate to measure the rotational angle output. A rolling needle bearing with an internal diameter of 100 mm (not shown) is installed in the output support shell to install and limit the motion of the output bracket. Since the output force sensor and VSA mechanism are independent of each other, the torque sensor is not shown in this figure.

Figure 4. Output frame diagram.

Figure 5. Three-dimensional schematic diagram of structural design.

2.4. Derivation and Analysis of the Dynamic Model

In the movable components, the main mass rotates on a fixed shaft. The rest are mainly the slider and the rolling needle bearings on it. The linkages with small mass are ignored. The mechanism can be simplified as shown in Figure 6, where $J_i,{\varphi }_i,T_i,i=1,2,4$ respectively represent the corresponding moment of inertia, rotation angle and input torque of the transmission from the motor to the linkage and the output bracket, and T_L represents the external load torque.

The force diagram of the slider of the four-bar linkage mechanism is shown in Figure 7, where F₁ and F₂ represent the forces acting on the slider by the two linkages, F₃, F_s, and F_f represent the radial force, vertical force and friction force generated by the plate spring on the slider. J_h and m₃ represent the moment of inertia of the holding bracket and the mass of the slider, and T₁₁ and T₂₂ represent the torques transferred from the motor output to the linkage, respectively:

Based on Figure 6 and Figure 7, the transmission Equation (1) of motor 1 and motor 2 and the force Equation (2) of the linkage are listed below:

$\begin{array}{l}{J}_1{\ddot{\varphi }}_1=T_1-T_{11}-d_1{\dot{\varphi }}_1\\ J_2{\ddot{\varphi }}_2=T_2-T_{22}-d_2{\dot{\varphi }}_2\\ J_4{\ddot{\varphi }}_4=T_3-T_L-d_4{\dot{\varphi }}_4\end{array}$ (1)

Figure 6. Schematic diagram of mechanism distribution.

Figure 7. Force diagram of the four-bar linkage mechanism.

$\begin{array}{l}{F}_1=\frac{T_{11}}{\sin \left({\varphi }_1-{\varphi }_2\right)l}\\ F_2=\frac{T_{22}}{\sin \left({\varphi }_1-{\varphi }_2\right)l}\end{array}$ (2)

The total mass of the four plate springs is less than 50 g, so the dynamic effects of the plate springs can be ignored. Considering that the slider is rotating as well as moving radially, the following dynamic equations are written for the slider and the support bracket:

$\left(m_3{r}^2+\frac{J_H}{n}\right)\left(\frac{{\ddot{\varphi }}_1+{\ddot{\varphi }}_2}{2}\right)=\frac{T_{\text{11}}+T_{\text{22}}}{n}-2m_{3}r\dot{r}\left(\frac{{\dot{\varphi }}_1+{\dot{\varphi }}_2}{2}\right)-F_{3}r$(3)

$m_3\ddot{r}=m_3{\left(\frac{{\dot{\varphi }}_1+{\dot{\varphi }}_2}{2}\right)}^{2}r+\left(F_2-F_1\right)\cos \left(\frac{{\varphi }_1-{\varphi }_2}{2}\right)-\mu F_3-F_s$ (4)

where n = 4 represents the number of parallel linkages, m₃ represents the mass of the slider, and J_H is the moment of inertia of the support bracket. The radius equation of the contact point is established based on the geometric relationship：

$r=2l\cos \left(\frac{{\varphi }_1-{\varphi }_2}{2}\right)$ (5)

Considering that the stiffness value is actually a function of torque and radius r, the torsional stiffness K_T of the VSA mechanism is defined as follows:

$K_T=f_K\left(T_3,r\right)$ (6)

In Equation (4), F_s represents the downward force of the variable plate spring on the slider. This force is essentially influenced by the load T_L and the stiffness value K_T of the VSA mechanism, but can be indirectly calculated by the energy required for the plate spring to change its stiffness being equal to the energy output change of the external system, without considering the friction. The calculation process is shown in Equation (7):

$\begin{array}{c}{F}_{s}gdr=T_{3}gd{\varphi }_4=T_3\left(\frac{T_3}{f_K\left(T_3,r+dr\right)}-\frac{T_3}{f_K\left(T_3,r\right)}\right)\\ \approx T_3^2\frac{f_K\left(T_3,r\right)-f_K\left(T_3,r+dr\right)}{f_K{\left(T_3,r\right)}^2}=-T_3^2\frac{\frac{\partial f_K\left(T_3,r\right)}{\partial r}}{f_K{\left(T_3,r\right)}^2}dr\end{array}$ (7)

The compensation equation for stiffness variation is obtained as follows:

$F_s=-\frac{T_3^2}{K_T^2}\frac{\partial K_T}{\partial r}$ (8)

Finally, the relationship between output torque, contact point radius, and rotation angle is listed as follows:

$T_{out}=4F_{3}r=K_T\left(\frac{{\varphi }_1+{\varphi }_2}{2}-{\varphi }_4\right)$ (9)

There are two problems in the dynamics of this VSA mechanism. One is that the magnitude of the differential mode torque required for the radial motion control of the slider is significantly different from the common mode torque required for the overall mechanism output. Secondly, the stiffness value of the plate spring shows a load-related characteristic during operation, which will be validated in the subsequent simulation. This section will only analyze the force relationship between various torque components affecting the core slider, which mainly corresponds to Equation (3) and Equation (4).

In Equation (3), the radial rotation of the slider is influenced by inertial force, Coriolis force, plate spring reaction force, and the input moment of the two linkages. Calculated based on the maximum working radius of 32 mm, the moment of inertia is only 15 × 10⁻⁶ kg∙m². Assuming that the maximum stiffness variation to the minimum stiffness requires 0.1 s, and considering that the output speed is under 15 rad/s, the maximum torque produced by the Coriolis force is about 0.03 N.m. The total output torque of the two motors in parallel is 28 N∙m, which is evenly distributed to the four executing mechanisms, totaling 7 N∙m. It can be seen that the effects of inertial force and Coriolis force on the input force are relatively small. Therefore, to avoid loss of generality, and considering the worst case, Equation (3) can be simplified as follows:

$J_3\left(\frac{{\ddot{\varphi }}_1+{\ddot{\varphi }}_2}{2}\right)=\frac{T_r}{n}-2m_{3}r{\dot{r}}_0\left(\frac{{\dot{\varphi }}_1+{\dot{\varphi }}_2}{2}\right)-F_{3}r$ (10)

where $J_3=m_3{r}_0^2+\frac{J_H}{n}=15\times {10}^{-6}\text{kg}\cdot {\text{m}}^2$ ，$T_r=T_{\text{11}}+T_{\text{22}}$ , ${\dot{r}}_0=0.15\text{m}/\text{s}$.

When the slider is moving radially, the load force is mainly composed of centrifugal force, friction force, inertial force, and downward pressure of the plate spring according to Equation (2-2). Based on the maximum speed and maximum radius, the centrifugal force is 0.47 N. Based on the maximum load force, thecontact force is $F_3=\frac{\text{20}\text{N}\cdot \text{m}}{4\times 0.03\text{m}}\simeq 167\text{N}$ . Based on the worst rolling friction coefficient of 0.003, the value is approximately 0.5 N. Since the mass is only $m_3=70\text{g}$ , assuming that the acceleration from rest to the maximum angular velocity requires 0.1 s, the required inertial force is 7 N. The compensation force for variable stiffness depends on the current load torque and the difference between the motor rotation angle and has a numerical change between 0 and several tens of N. Among the four forces, the first three are transformed into the maximum difference in motor output torque according to Equations (2-1) and (2-2), which is 0.025 N∙m. But the downward pressure of the plate spring has different stiffness values, so the differential torque requires a value between 0 and 3 N∙m. The difference between the two values is large, so it is very important to obtain the relationship between the stiffness KT, contact point radius, and load as a function. It cannot be simplified directly.

3. CPG Human Motion Phase Recognition Algorithm

3.1. Convergence Principle of CPG Recognition Algorithm

The convergence characteristics of CPG are based on the adaptive frequency learning method proposed by Ludovic et al. The basic structure under the Cartesian coordinate system is shown in Equation (11):

$\begin{array}{l}\dot{x}=f\left(x,y,\omega \right)+\epsilon F\left(t\right)\\ \dot{y}=f\left(x,y,\omega \right)\\ \dot{\omega }=-\epsilon F\left(t\right)\frac{y}{\sqrt{x^2+y^2}}\end{array}$ (11)

where $x,y$ represents a two-dimensional limit cycle oscillator, and $F\left(t\right)$ represents external disturbances. In the above equation, $\epsilon F\left(t\right)$ represents phase learning, and $-\epsilon F\left(t\right)\frac{y}{\sqrt{x^2+y^2}}$ represents frequency learning behavior of external signals. The learning result is that variable x can basically reproduce the main frequency component of the input signal $F\left(t\right)$ .

This form is widely used in Hopf oscillators. On the one hand, the structure of the Hopf oscillator is simple, and on the other hand, the limit cycle is a uniform circular motion, which can correspond well with the frequency spectrum. Its polar coordinate equation has more obvious mathematical significance, as shown in Equation (12):

$\begin{array}{l}\dot{r}=\left(\mu -r^2\right)r+\epsilon F\cos \phi \\ \dot{\phi }=\omega -\frac{\epsilon }{r}F\sin \phi \\ \dot{\omega }=-\epsilon F\sin \phi \end{array}$ (12)

Here, $r,\phi ,\omega$ represent the radius, phase, and angular velocity of the limit cycle respectively. It can be seen that when the disturbance coefficient $\epsilon =0$ , Equation (12) degenerates into a complete Hopf oscillator. When $\epsilon \ne 0$ , it will produce disturbances to the system and continue to affect the parameters of the oscillator.

The disturbance behavior can be visually explained in Figure 8, where the middle part shows a Hopf oscillator with a circular limit cycle, and the black dot represents the current position. The direction of the limit cycle’s evolution in the next moment under the limit cycle is considered as the vector $P$ . At this time, the external interference $F$ is viewed as a rotating vector around the origin, and the angle and amplitude of the vector represent the phase angle and signal strength of the external interference, respectively. For example, a constant disturbance $F=0$ is represented as a constant vector starting from the origin with a phase angle of 0.

The influence of the disturbance on the limit cycle can be determined by the phase difference between the current external signal and the phase of the current position on the limit cycle, which can be divided into two directions: tangential $e_r$ and normal $e_n$ as shown in Figure 9.

Figure 8. Representation of external disturbance limit cycle.

Figure 9. Tangential and orthogonal impacts of external disturbance on the limit cycle.

To obtain the projection values of $F$ in the tangential and normal directions, the polar coordinate system is transformed back into the Cartesian coordinate system, and the values are obtained as $F\cdot \frac{y}{\sqrt{x^2+y^2}}{e}_r$ and $F\cdot \frac{x}{\sqrt{x^2+y^2}}{e}_n$ , respectively. It can be further written as $F\sin \phi$ and $F\cos \phi$ .

Back to Equation (12), it can be seen that the amplitude term is affected by disturbances along the vertical direction, and once the disturbance is removed, the radius of the limit cycle is still determined by $\mu$ . However, the phase and angular velocity are affected by disturbances along the tangential direction and will change significantly. Therefore, under the condition of only requiring synchronization characteristics, the equation for radius change can be directly omitted, and only the last two terms in Equation (12) that affect the phase are used. It has been proved in Ludovic’s literature that when fixing the radius of the limit cycle to be 1, the ability of synchronous learning has been demonstrated. The following will be generalized to discuss the convergence pattern when the limit cycle is not a circle.

3.2. Convergence Analysis of the CPG Algorithm

The existing literature has shown that the amplitude of the limit cycle does not have a significant impact on the phase. Removing the amplitude and continuing the previous approach in the tangential direction, the following simplified Equation (13) can be obtained:

$\begin{array}{l}\dot{\varphi }=\omega +\epsilon F\left(t\right)g\left(\varphi \right)\\ \dot{\omega }=\epsilon F\left(t\right)g\left(\varphi \right)\end{array}$ (13)

where $\sin \left(\varphi \right)$ representing tangential direction previously is now represented by $g\left(\varphi \right)$ , the derivative of $f\left(\varphi \right)$ with respect to $\varphi$ . In practical applications, $f\left(\varphi \right)$ represents the waveform of a gait for human walking. The time axis of this waveform is normalized in the interval $\left[0,1\right]$ , and the corresponding angle $\varphi$ also changes within the interval $\left[0,1\right]$ . is the derivative of normalized $f\left(\varphi \right)$ with respect to $\varphi$ . $F\left(t\right)$ represents the input signal, which is the joint angle signal.

3.2.1. Explanation of Disturbance Theory Explanation

The disturbance theory has a good estimation effect on the evolution of weakly nonlinear systems. Its basic principle is to perform a Taylor expansion of the nonlinear system at the current moment, and then calculate and estimate based on the expansion term with the strongest influence. Therefore, it is very suitable for studying the learning and convergence behavior of the simple limit cycle CPG in this paper.

As with most linearization methods for nonlinear systems, the estimation is often more accurate in a small neighborhood of the linearization point. Whether the long-term estimation error converges depends on the type of equilibrium point of the system without disturbance. For disturbance theory, the equilibrium point of the oscillator is a non-exponential equilibrium point, which means the long-term estimation based on a certain operating point diverges. As a substitute, observing and studying the current expansion at a time point can determine the evolution trend of the current system. Then, the convergence of all operating points can be concluded, and the global convergence property can be indirectly obtained.

Back to the oscillator model, the theory is mainly used to analyze two points: 1) the convergence characteristics of phase when the frequency and waveform of the limit cycle are consistent with those of the external signal; 2) the convergence characteristics of the limit cycle frequency.

3.2.2. Proof of Convergence of Phase at the Same Frequency

According to the hypothesis, the oscillator equation is rewritten as follows:

$\dot{\varphi }=\omega +\epsilon F\left(t\right)g\left(\varphi \right)$ (14)

According to the theoretical introduction, the first-order estimation equation for the system state variable $\varphi$ is given:

$\varphi ={\varphi }_0+\epsilon {\varphi }_1+{\epsilon }^2{R}_{\varphi }$ (15)

Here, 0-th order estimation, 1st order estimation, and the residual error of the quadratic term of $\varphi$ are represented by ${\varphi }_0,{\varphi }_1,R_{\varphi }$ , respectively. Generally, $f\left(\varphi \right)$ is obtained by filtering signal $F\left(t\right)$ , so $f\left(\varphi \right)$ and $g\left(\varphi \right)$ can be represented as a finite sum of sinusoids:

$\begin{array}{l}f\left(\varphi \right)={\displaystyle \sum }_{n\in Z}{A}_n{\text{e}}^{in\varphi }\\ g\left(\varphi \right)={\displaystyle \sum }_{n\in Z}in{A}_n{\text{e}}^{in\varphi }\end{array}$ (16)

To obtain a multi-order series expansion of $g\left(\varphi \right)$ with respect to $\varphi$ , we first expand the functions $\cos \varphi$ and $\sin \varphi$ as shown in the following Equation (17) and Equation (18):

$\begin{array}{l}\cos \left[kn\left({\varphi }_0+\epsilon {\varphi }_1+{\epsilon }^2{\varphi }_2+{\epsilon }^3{R}_{\varphi }\right)\right]\\ ={\displaystyle \sum }_{p=0}^{\infty }\frac{{\left(-1\right)}^p{\left(kn{\varphi }_0+\epsilon kn{\varphi }_1+{\epsilon }^{2}kn{\varphi }_2+{\epsilon }^{3}kn{R}_{\varphi }\right)}^{2p}}{\left(2k\right)!}\\ =\cos \left(kn{\varphi }_0\right)-\epsilon kn{\varphi }_1\sin \left(kn{\varphi }_0\right)+O\left({\epsilon }^2\right)\end{array}$ (17)

$\begin{array}{l}\sin \left[k\left({\varphi }_0+\epsilon {\varphi }_1+{\epsilon }^2{\varphi }_2+{\epsilon }^3{R}_{\varphi }\right)\right]\\ ={\displaystyle \sum }_{p=0}^{\infty }\frac{{\left(-1\right)}^p{\left(k{\varphi }_0+\epsilon k{\varphi }_1+{\epsilon }^{2}k{\varphi }_2+{\epsilon }^{3}k{R}_{\varphi }\right)}^{2p+1}}{\left(2p+1\right)!}\\ =\sin \left(k{\varphi }_0\right)+\epsilon k{\varphi }_1\cos \left(k{\varphi }_0\right)+O\left({\epsilon }^2\right)\end{array}$ (18)

We can then derive Equation (19):

$\begin{array}{c}{\text{e}}^{ik\varphi }=\cos \left(k\varphi \right)+i\sin \left(k\varphi \right)\\ =\cos \left(k{\varphi }_0\right)-\epsilon k{\varphi }_1\sin \left(k{\varphi }_0\right)+O\left({\epsilon }^2\right)\\ +i\left[\sin \left(kn{\varphi }_0\right)+\epsilon kn{\varphi }_1\cos \left(kn{\varphi }_0\right)+O\left({\epsilon }^2\right)\right]\\ ={\text{e}}^{ik{\varphi }_0}+\epsilon ik{\varphi }_1{\text{e}}^{ik{\varphi }_0}+O\left({\epsilon }^2\right)\end{array}$ (19)

Finally, by substituting formula Equation (19) into Equation (14), we can obtain Equation (20):

${\dot{\varphi }}_0+\epsilon {\dot{\varphi }}_1+{\epsilon }^2{\dot{R}}_{\varphi }={\omega }_0+\epsilon F\left(t\right)g\left({\varphi }_0\right)+{\epsilon }^2{R}_{\varphi }$ (20)

Following the principle of equal corresponding items in disturbance theory, we can obtain the following formula Equation (21):

$\begin{array}{l}{\dot{\varphi }}_0=\omega \\ {\dot{\varphi }}_1=F\left(t\right)g\left({\varphi }_0\right)\end{array}$ (21)

Suppose we integrate from time $t_0$ to obtain the following solution formula Equation (22):

$\begin{array}{l}{\varphi }_0\left(t\right)=\omega t+\omega t_0\\ {\varphi }_1\left(t\right)={\displaystyle {\int }_{t_0}^{t}F\left(t\right)g\left({\varphi }_0\right)\text{d}t}\end{array}$ (22)

Since $f\left(\varphi \right)$ must have the main waveform of the input signal $F\left(t\right)$ for phase synchronization to be meaningful, we can represent the input signal as Equation (23):

$F\left(t\right)=f\left(\omega t\right)+R\left(t\right)$ (23)

Here, $R\left(t\right)$ represents the residual high-frequency signal. The phase difference between $f\left(\varphi \right)$ and $F\left(t\right)$ is already implicit in the starting time point $t_0$ of the limit cycle, which is represented by $\omega t_0$ . At this time, $g\left({\varphi }_0\right)$ can be obtained by Equation (24):

$g\left({\varphi }_0\right)={\displaystyle \sum }_{n\in Z}in\omega A_n{\text{e}}^{in\omega \left(t-t_0\right)}$ (24)

By substituting Equation (24) into Equation (22) and integrating ${\varphi }_1\left(t\right)$ , we can obtain the following Equation (25):

$\begin{array}{c}{\varphi }_1\left(t\right)={\displaystyle {\int }_{t_0}^t{\displaystyle \sum }_{m\in \mathbb{Z}}{A}_m{\text{e}}^{im\omega t}\cdot {\displaystyle \sum }_{n\in \mathbb{Z}}in\omega A_n{\text{e}}^{in\left(\omega t-\omega t_0\right)}+R\left(t\right)g\left({\varphi }_0\right)\text{d}t}\\ ={\displaystyle \sum }_{m\in \mathbb{Z},n\in \mathbb{Z}}{\displaystyle {\int }_{t_0}^{t}in\omega A_m{A}_n{\text{e}}^{im\omega t+in\left(\omega t-\omega t_0\right)}\text{d}t}+\tilde{R}\left(t\right)\end{array}$(25)

At this time, a constant term appears only when $m=-n$ . Otherwise, the effects will be averaged over a sufficiently long time. Assuming the input signal is a real signal, i.e. $A_n$ and $A_{-n}$ are conjugates, the result can be further simplified as Equation (26):

$\begin{array}{c}{\varphi }_1\left(t\right)={\displaystyle \sum }_{m=-n,n\in Z}{\displaystyle {\int }_{t_0}^{t}i\omega n{A}_{-n}{A}_n{\text{e}}^{-i\omega n{t}_0}\text{d}t}+\tilde{R}\left(t\right)\\ ={\displaystyle \sum }_{n\in N}{\displaystyle {\int }_{t_0}^{t}i\omega n{A}_{-n}{A}_n{\text{e}}^{-i\omega n{t}_0}-i\omega n{A}_{-n}{A}_n{\text{e}}^{i\omega n{t}_0}\text{d}t}+\tilde{R}\left(t\right)\\ ={\displaystyle \sum }_{n\in N^{*}}2n\omega {\left|A_n\right|}^2\sin \left(n\omega t_0\right)\left(t-t_0\right)+\tilde{R}\left(t\right)\end{array}$(26)

Note that the constant component in the convergence result does not affect the convergence result, only non-zero frequencies will affect it, which is slightly different from the subsequent frequency convergence derivation. Observing the constant term, when $\omega t_0>0$ , it represents that the oscillator signal phase is lagging behind the external signal. At this time, it will accelerate the phase $\varphi$ to approach the external signal, which will gradually reduce $\omega t_0$ to 0. The opposite is also true. Taking the signal $F\left(t\right)=\sin \left(\pi t\right)+\sin \left(2\pi t\right)$ as an example, the function graph of the constant component $f\left(\varphi \right)=2\pi \sin \left(\varphi \right)+4\pi \sin \left(2\varphi \right)$ is plotted according to Equation (26), as shown in Figure 10:

Figure 10. Schematic diagram of phase convergence.

In the evolution direction, a positive value means that the oscillator will accelerate its own frequency, while a negative value means the opposite. A positive phase difference means that the oscillator phase leads the perturbation signal, while the amplitude is opposite. Looking within the domain P₀, when the internal phase of the oscillator is in phase with the external signal, the oscillator will accelerate its own frequency and reduce the phase lag. The opposite is also true. This drives the oscillator signal to always converge towards P₀, the phase synchronization point.

In fact, each zero point where the derivative of the curve in Figure 10 is negative is a convergence point, while the zero points where the derivative is positive determine the range and boundaries of each convergence interval. As can be seen from Figure 10, the convergence to the zero points between every two blue lines only occurs in the area where the phase difference is approximately ± 2 rad.

For continuous spectrum external signal inputs, the constant component in Eq. 26 can be rewritten in a continuous integration form:

$-{\displaystyle {\int }_0^{\infty }2\omega A_{\omega }^2\sin \left(\omega t_0\right)\text{d}\omega }$ (27)

Taking the hip joint as an example, the amplitude distribution of the hip joint signal after FFT is shown in Figure 11. It can be seen that the spectrum is mainly concentrated in two parts: walking frequency and mean value. Equation (27) indicates that the mean value does not affect the phase convergence, so it can be predicted that the phase will converge when the hip joint angle signal is inputted as shown in Figure 11. Substituting the hip joint FFT amplitude results into Equation (27), we can draw the oscillator evolution trend chart for the hip joint angle signal as an input in Figure 12.

Because the external signal has not been time-normalized, the phase difference in Figure 11 here represents the synchronization time difference. The normal walking speed of the human body is about 0.8 Hz. Obviously, the convergence points are lagging one gait cycle, no lag, and leading one gait cycle. Drawing the

Figure 11. Hip joint angle spectrum during walking at 4.5 Km/h on flat ground.

Figure 12. Convergence trend of hip joint angle curve.

convergence interval, it can be seen that the convergence result is either completely synchronized with the signal, or leading or lagging an integer number of gait cycle times. This phenomenon is determined by the fact that the hip joint angle curve spectrum is highly concentrated and close to a sine curve. In fact, using only Equation (27) to predict the convergence situation is not entirely correct. When there are multiple frequency components, the formula will obtain multiple intervals and different phase-locked points in one signal cycle. Obviously, this is incorrect. The hip joint signal spectrum has an infinite bandwidth, but the fact is that it can still converge in phase. The reason is that other terms with mean effects are not irrelevant during the convergence process. These terms reveal that there are phase back-and-forth oscillations whether in the convergence process or after convergence. These oscillations allow the phase point to jump out of the convergence interval caused by small amplitude frequency components and converge in the interval of the main frequency component. This will be further discussed in the frequency convergence later. In conclusion, it can be inferred that when the hip joint angle signal is input, the phase difference under the same frequency can completely converge to an integer number of gait cycles. Therefore, as long as the frequency converges, the phase must converge, which provides a prerequisite for the next proof.

3.2.3. Frequency Convergence Derivation

Since the disturbance theory is a linearized estimation, the more complex the limit cycle trajectory $f\left(\varphi \right)$ is, the worse the estimation effect and the more serious the error. Integrating the result will be extremely difficult to calculate. Therefore, only the circular limit cycle will be calculated next, and its conclusion will be extended to complex limit cycles. For Equation (13), a third-order estimation formula is listed as below:

$\begin{array}{l}\varphi ={\varphi }_0+\epsilon {\varphi }_1+{\epsilon }^2{\varphi }_2+{\epsilon }^3{R}_{\varphi }\\ \omega ={\omega }_0+\epsilon {\omega }_1+{\epsilon }^2{\omega }_2+{\epsilon }^3{R}_{\omega }\end{array}$ (28)

Here, $\varphi ,\omega ,{\varphi }_i,{\omega }_i,i=1,2$ and $R_{\varphi },R_{\omega }$ are dependent variables of time t. Substituting them into formula (3-18) and combining the above formulas Equation (17), Equation (18) and Equation (19), we derive Equation (29):

$\begin{array}{l}{\dot{\varphi }}_0+\epsilon {\dot{\varphi }}_1+{\epsilon }^2{\dot{\varphi }}_2+{\epsilon }^3{\dot{R}}_{\varphi }={\omega }_0+\epsilon \left[{\omega }_1+F\left(t\right)g\left({\varphi }_0\right)\right]+{\epsilon }^{2}F\left(t\right){\omega }_2{\varphi }_1{g}^{\prime }\left({\varphi }_0\right)\\ {\dot{\omega }}_0+\epsilon {\dot{\omega }}_1+{\epsilon }^2{\dot{\omega }}_2+{\epsilon }^3{\dot{R}}_{\omega }=\epsilon F\left(t\right)g\left({\varphi }_0\right)+{\epsilon }^{2}F\left(t\right){\varphi }_1{g}^{\prime }\left({\varphi }_0\right)\end{array}$ (29)

Similarly, based on the disturbance theory, Equation (30) can be written according to the corresponding equal terms as below:

$\begin{array}{l}{\dot{\varphi }}_0={\omega }_0\\ {\dot{\omega }}_0=0\\ {\dot{\varphi }}_1={\omega }_1+F\left(t\right)g\left({\varphi }_0\right)\\ {\dot{\omega }}_1=F\left(t\right)g\left({\varphi }_0\right)\\ {\dot{\varphi }}_2=F\left(t\right){\omega }_2{\varphi }_1{g}^{\prime }\left({\varphi }_0\right)\\ {\dot{\omega }}_2=F\left(t\right){\varphi }_1{g}^{\prime }\left({\varphi }_0\right)\end{array}$ (30)

Next, we will solve Equation (30), assuming that the input is a finite number of sine waves with a discontinuous spectrum. The signal is represented as formula (31):

$F\left(t\right)={\displaystyle \sum }_{n\in \mathbb{Z}}{A}_n{\text{e}}^{in{\omega }_{F}t}$ (31)

The corresponding limit cycle trajectory of the circle is written as below:

$f\left(\varphi \right)=a\sin \left(\varphi +{\varphi }_0\right)$ (32)

Equation (13) can be written as follows:

$\begin{array}{l}\dot{\varphi }=\omega +\epsilon F\left(t\right)a\cos \left(\varphi +{\varphi }_0\right)\\ \dot{\omega }=\epsilon F\left(t\right)a\cos \left(\varphi +{\varphi }_0\right)\end{array}$ (33)

Firstly, calculate the first-order term as follows:

$\begin{array}{l}{\omega }_0\left(t\right)={\omega }_0\\ {\varphi }_0\left(t\right)={\omega }_0\left(t-t_0\right)\end{array}$ (34)

Assuming $n{\omega }_F\ne {\omega }_0$ , that is, the starting frequency is not synchronized with any frequency of the input signal $F\left(t\right)$ , the second-order term is calculated in sequence as shown in Equation (35) and Equation (36):

$\begin{array}{l}{\omega }_1\left(t\right)={\displaystyle {\int }_{t_0}^{t}F\left(t\right)f^{\prime }\left({\omega }_0\left(t\right)\right)\text{d}t}\\ ={\displaystyle \sum }_{n\in \mathbb{Z}}\frac{A_{n}a{\omega }_0}{2}\left(\frac{{\text{e}}^{i\left(n{\omega }_{F}t+{\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)}-{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{{\text{e}}^{i\left(n{\omega }_{F}t-{\omega }_{0}t+{\omega }_0{t}_0-{\varphi }_0\right)}-{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{i\left(n{\omega }_F-{\omega }_0\right)}\right)\end{array}$ (35)

$\begin{array}{c}{\varphi }_1={\displaystyle {\int }_{t_0}^t{\omega }_1\text{d}t}+{\omega }_1\\ ={\displaystyle \sum }_{n\in Z}\frac{A_{n}a{\omega }_0}{2}[\left(\frac{1}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{1}{{\left(i\left(n{\omega }_F+{\omega }_0\right)\right)}^2}\right){\text{e}}^{i\left(n{\omega }_{F}t+{\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)}\\ -\left(\frac{1}{i\left(n{\omega }_F-{\omega }_0\right)}+\frac{1}{{\left(i\left(n{\omega }_F-{\omega }_0\right)\right)}^2}\right){\text{e}}^{i\left(n{\omega }_{F}t-\left({\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)\right)}\\ -\left(\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{i\left(n{\omega }_F-{\omega }_0\right)}\right)\left(t-t_0\right)\\ -\left(\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{i\left(n{\omega }_F-{\omega }_0\right)}-\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{{\left(n{\omega }_F+{\omega }_0\right)}^2}-\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{{\left(n{\omega }_F-{\omega }_0\right)}^2}\right)]\end{array}$ (36)

As mentioned earlier, frequency convergence is a prerequisite for phase convergence. The first-order estimation result ${\omega }_1\left(t\right)$ only has oscillations with a mean value of 0, which cannot play a decisive role in determining the evolution direction of the frequency. Therefore, it is necessary to further expand the CPG unit and estimate the second-order items for the convergence frequency. Equation (36) is substituted into Equation (30) for further integration of ${\omega }_2\left(t\right)$ , and the calculation result is shown in Equation (37):

$\begin{array}{c}{\dot{\omega }}_2=F\left(t\right)f\left({\varphi }_0\left(t\right)\right){\varphi }_1\\ ={\displaystyle \sum }_{m\in Z}-a{\omega }_0^2{A}_m{\text{e}}^{im{\omega }_{F}t}\sin \left({\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right){\varphi }_1\\ ={\displaystyle \sum }_{m\in Z}\frac{-a{\omega }_0^2{A}_m}{2i}{\varphi }_1\left({\text{e}}^{i\left(m{\omega }_{F}t+{\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)}-{\text{e}}^{i\left(m{\omega }_{F}t-{\omega }_{0}t+{\omega }_0{t}_0-{\varphi }_0\right)}\right)\\ ={\displaystyle \sum }_{m\in Z}{\displaystyle \sum }_{n\in Z}\frac{-a^2{\omega }_0^3{A}_m{A}_n}{4i}\left(E_1+E_2+E_3+E_4\right)\end{array}$(37)

The expressions of E₁ to E₅ are shown in Equation (38) as follows:

$\begin{array}{l}{E}_1=\left(\frac{1}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{1}{i^2{\left(n{\omega }_F+{\omega }_0\right)}^2}\right){\text{e}}^{i\left(\left(n+m\right){\omega }_{F}t+2\left({\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)\right)}\\ E_2=-\left(\frac{1}{i\left(n{\omega }_F-{\omega }_0\right)}+\frac{1}{i^2{\left(n{\omega }_F-{\omega }_0\right)}^2}\right){\text{e}}^{i\left(\left(n+m\right){\omega }_{F}t-2\left({\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)\right)}\\ E_3=\left(\frac{2{\omega }_0}{i\left(n{\omega }_F^2-{\omega }_0^2\right)}+\frac{4n{\omega }_F{\omega }_0}{{\left(n{\omega }_F^2-{\omega }_0^2\right)}^2}\right){\text{e}}^{i\left(n+m\right){\omega }_{F}t}\\ E_4=\left(\frac{i{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{\left(n{\omega }_F+{\omega }_0\right)}+\frac{i{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{\left(n{\omega }_F-{\omega }_0\right)}\right)\left({\text{e}}^{i\left(m{\omega }_{F}t+{\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)}-{\text{e}}^{i\left(m{\omega }_{F}t-{\omega }_{0}t+{\omega }_0{t}_0-{\varphi }_0\right)}\right)\left(t-t_0\right)\\ E_5=-\left({\text{e}}^{i\left(m{\omega }_{F}t+{\omega }_{0}t-{\omega }_0{t}_0+{\varphi }_0\right)}-{\text{e}}^{i\left(m{\omega }_{F}t-{\omega }_{0}t+{\omega }_0{t}_0-{\varphi }_0\right)}\right)(\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{i\left(n{\omega }_F+{\omega }_0\right)}+\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{i\left(n{\omega }_F-{\omega }_0\right)}\\ -\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0+{\varphi }_0\right)}}{{\left(n{\omega }_F+{\omega }_0\right)}^2}-\frac{{\text{e}}^{i\left(n{\omega }_F{t}_0-{\varphi }_0\right)}}{{\left(n{\omega }_F-{\omega }_0\right)}^2})\end{array}$ (38)

$E_1{E}_2{E}_5$ are periodic oscillations with a mean value of 0 after time integration. $E_4$ represents the divergent term, which gradually increases in amplitude over time, but the mean value is still 0 in the integration process. Only the integration mean of $E_3$ when $m=n$ is not 0, leading to a trend of fixed direction change. The simplified expression of ${\omega }_2$ is as follows:

${\omega }_2={\tilde{\omega }}_2+\frac{a^2{\omega }_0^3}{2}\left({\displaystyle \sum }_{n\in N^{*}}\frac{2{\left|A_n\right|}^2{\omega }_0}{\left({\left(n{\omega }_F\right)}^2-{\omega }_0^w\right)}-\frac{A_0^2}{{\omega }_0}\right)\left(t-t_0\right)={\tilde{\omega }}_2+{\omega }_{bias}\left(t-t_0\right)$ (39)

Here, ${\tilde{\omega }}_2$ represents oscillation components with a mean value of 0. At this time, the second-order estimation of the frequency is as follows:

$\omega ={\omega }_0+\epsilon {\omega }_1+{\epsilon }^2{\tilde{\omega }}_2+{\epsilon }^2{\omega }_{bias}\left(t-t_0\right)+{\epsilon }^3{R}_{\omega }$ (40)

The mean value of the remaining components’ convergence contribution to the frequency is 0, so the numerical value of ${\omega }_{bias}$ will mainly determine the direction of system evolution. Next, the frequency convergence property of CPG will be analyzed based on the theory above.

3.2.4. Frequency Convergence Property

The hip joint angle curve is very close to a sine curve. When deriving and using it, it can be directly equivalent to a mean value plus a sine curve for convergence verification. However, for the convenience of parameter debugging or synchronization with other complex signals, this section summarizes the convergence laws, and introduces the basic convergence properties, continuous signal convergence properties, disturbance effects during convergence, and consistent convergence properties of complex trajectory phases, supplemented by simulation verification.

1) Basic Convergence Property

When the initial frequency of the oscillator ${\omega }_0=0$ , ${\omega }_{bias}=0$ , Equation (3-29) indicates that the oscillator cannot learn. If ${\omega }_0\ne 0$ , the convergence direction will be determined by ${\omega }_{bias}$ .

The denominator of the second term $\frac{1}{{\left(n{\omega }_F\right)}^2-{\omega }_0^2}$ in ${\omega }_{bias}$ shows that each frequency component $n{\omega }_F$ of the input signal will have a tendency to pull oscillator frequency ${\omega }_0$ towards it, and the closer it is, the stronger the pull will be, approaching infinity. Therefore, we can draw the following conclusion:

When ${\omega }_0\ne 0$ , there is a neighborhood for each frequency component of the disturbance signal, as long as the angular frequency of the oscillator is close enough to the frequency component, the oscillator frequency component will definitely converge to this frequency. The boundary between each neighborhood is determined by solving Equation (41) below.

$-\frac{A_0^2}{{\omega }_0}+{\displaystyle \sum }_{n\in {\mathbb{N}}^{*}}\frac{2{\left|A_n\right|}^2{\omega }_0}{{\left(n{\omega }_F\right)}^2-{\omega }_0^2}=0$ (41)

Assuming that the input signal is composed of n non-zero frequency real signals, solving Equation (41) with ${\omega }_0$ as the variable yields $n-1$ solutions that divide the positive real axis into n open intervals. Each interval must contain a frequency of the input signal, and as long as the initial frequency of the oscillator is in this interval, it will eventually converge to this frequency.

In practical applications, the influence of the constant component often exceeds what is described in Equation (38). Through experimental testing and adjustment, Equation (42) can be rewritten as:

$-\frac{4A_0^2}{{\omega }_0}+{\displaystyle \sum }_{n\in {\mathbb{N}}^{*}}\frac{2{\left|A_n\right|}^2{\omega }_0}{\left({\left(n{\omega }_F\right)}^2-{\omega }_0^2\right)}=0$ (42)

Simulation verification will be carried out by taking an example. Assuming the input signal is $F\left(t\right)=3+2\sin \left(\pi t\right)+\sin \left(2\pi t\right)$ , substituting it into Equation (42) yields the Equation (43):

$-\frac{4\times 3^2}{{\omega }_0}+\frac{2\times 2^2\times {\omega }_0}{{\pi }^2-{\omega }_0^2}+\frac{2\times 1^2\times {\omega }_0}{{\left(2\pi \right)}^2-{\omega }_0^2}=0$ (43)

By discarding the negative solutions, we obtain 2.8380 and 6.1531. The positive real axis is divided into three intervals: (0, 28380), (2.8380, 6.1531), (6.1531, +∞). In the first interval, the initial frequency converges to the 0 frequency, in the second interval it converges to π, and in the third interval it converges to 2π.

The convergence results with different initial frequencies are shown in Figure 13. It can be clearly seen that there are three convergence results corresponding to the distribution of perturbation signal frequencies, and the boundaries between the convergence results basically match the calculated results. The final convergence frequency results in each convergence interval are close to the lower limit of the interval, indicating that there is a superposition effect among multiple frequencies of the perturbation signal on the oscillator frequency convergence, which is consistent with the description in Equation (43) and verifies the above predicted results.

Figure 13. Convergence conditions with different initial frequencies.

In the previous derivation process, the convergence assumption was made as $n{\omega }_F\ne {\omega }_0$ , and the Equation (38) that determines the convergence direction was obtained. It did not consider the possibility of the average values of E₁, E₂, E₄, and E₅ not being zero at some point during the convergence process. Using the input signal $f\left(\varphi \right)=\sin \left(\varphi \right)$ and learning signal mode $F\left(t\right)=\sin \left(\pi t\right)+\sin \left(2\pi t\right)$ , with an initial angular frequency of 4π rad/s and perturbation gain $\epsilon =0.1$ , a simulation experiment was conducted. The result is shown in Figure 14.

The oscillator finally converged to 3π, which is consistent with the predicted result in Equation (37). It should be noted that the convergence frequency passed through $\omega =3\pi$ rad/s. If perturbation analysis is carried out at this point, E₁ and E₂ are not zero averages. This indicates that E₁ and E₂ do not play a decisive role in convergence, and in fact, E₄ and E₅ cannot have an absolute impact on convergence. The main factor affecting convergence is E₃.

Figure 14. Convergence to the second-order special point chart.

The hip joint angle signal is essentially equivalent to the addition of a constant component and a frequency component. As long as the mean component is removed, and the oscillator for circular limit cycle CPG has a starting frequency that is not 0, according to Equation (3-30), the oscillator will definitely frequency converge to learn the angle signal, and then the phase will inevitably synchronize.

2) Nature of Continuous Perturbation Signal Input

The aforementioned are all discrete disturbance signals input. When the input signal is a continuous signal, i.e., $F\left(t\right)={\displaystyle {\int }_{\omega \in ℝ}F\left(\omega \right){\text{e}}^{i\omega t}\text{d}\omega }$ , the disturbance theory still holds. However, in the derivation process, the spectrum of the input signal ${\omega }_0$ is inevitably involved, which belongs to the input signal condition. At this time, only the contribution of E₃ is considered, and the E₃ term is rewritten in the form of continuous spectrum, as shown in the result Equation (44).

$\frac{a^2{\omega }_0^3}{2}\left[-\frac{4A_0^2}{{\omega }_0}+{\displaystyle {\int }_{0^{+}}^{{\omega }_0^{-}}\frac{2{\left|A_{\omega }\right|}^2{\omega }_0}{{\omega }^2-{\omega }_0^2}\text{d}\omega }+{\displaystyle {\int }_{+\infty }^{{\omega }_0^{+}}\frac{2{\left|A_{\omega }\right|}^2{\omega }_0}{{\omega }^2+{\omega }_0^2}\text{d}\omega }\right]$ (44)

The intuitive physical description of the above equation is shown in Figure 15. The curve represents the spectral amplitude of the disturbance signal. The part of the curve on the right side of ${\omega }_0$ will shift the inherent frequency of the CPG towards a higher frequency, which can be intuitively understood as a force F₂ to the right. Similarly, the spectral on the left side of ${\omega }_0$ will generate a force F₁ to the left.

Figure 15. Schematic diagram of frequency convergence.

The property that the convergence effect is greater for frequencies closer to ${\omega }_0$ is still maintained in Equation (44). In the vicinity of zero, the integral on both sides tends towards infinity in the same form. Assuming there is a small neighborhood ${\omega }_0\pm \Delta \omega$ , the sum of the influences of the two symmetrical points on the left and right on ${\omega }_0$ is given by Equation (44) as follows:

$\frac{2A_{{\omega }_0-\Delta \omega }^2{\omega }_0}{{\left({\omega }_0-\Delta \omega \right)}^2-{\omega }_0^2}+\frac{2A_{{\omega }_0+\Delta \omega }^2{\omega }_0}{{\left({\omega }_0+\Delta \omega \right)}^2-{\omega }_0^2}\frac{2A_{{\omega }_0-\Delta \omega }^2{\omega }_0}{{\left({\omega }_0-\Delta \omega \right)}^2-{\omega }_0^2}+\frac{2A_{{\omega }_0+\Delta \omega }^2{\omega }_0}{{\left({\omega }_0+\Delta \omega \right)}^2-{\omega }_0^2}$ (45)

It can be written as below since $\Delta \omega \ll {\omega }_0$ :

$\frac{\left(A_{{\omega }_0+\Delta \omega }-A_{{\omega }_0-\Delta \omega }\right)\left(A_{{\omega }_0+\Delta \omega }+A_{{\omega }_0-\Delta \omega }\right)}{\Delta \omega }={\frac{2A_{{\omega }_0}}{\Delta \omega }g\frac{\text{d}A}{\text{d}\omega }|}_{\omega ={\omega }_0}$ (46)

It can be seen that in the extremely small convergence region of the oscillator frequency, which is close to ${\omega }_0$ , the convergence effect of the two symmetric frequencies on the oscillator frequency is approximately determined by the amplitude $2A_{{\omega }_0}$ and the derivative of the amplitude with respect to the frequency of the input signal at this point. Moreover, it tends to infinity. When the disturbance signal frequency is not in the frequency domain of ${\omega }_0$ , the attraction to the oscillator frequency ${\omega }_0$ is limited, which means that the direction and speed of convergence are determined by the local properties of ${\omega }_0$ .

Specifically, there are two types of convergence directions. If the spectral amplitude at ${\omega }_0$ is not 0, and the derivative of the spectral amplitude of the input signal is not 0, then the positive or negative value of the derivative will determine the direction of CPG evolution. If the spectral amplitude at this point is 0, or the derivative is 0, then the convergence direction will be completely determined by the calculation result of integral 44. This rule also fully includes the convergence rule described by Equation (42).

The input signal is added with white noise with an energy of 0.01 to form a continuous spectrum. The frequency distribution diagram is shown in Figure 16, and the simulation results are shown in Figure 17.

In terms of convergence results, noise does not affect the convergence results. During the process, when the internal angular frequency of CPG is at 3π rad/s, the derivative of the spectral amplitude of the output signal at this point with

Figure 16. Distribution of initial spectrum.

Figure 17. Learning with noise.

respect to the angular frequency is 0. The input signal shows attraction between frequencies of 2π rad/s and π rad/s, which pulls the CPG frequency towards it. When it reaches 2π rad/s, if it continues to decrease below 2π rad/s, according to Equation (46), 2π rad/s will produce an infinite positive value, which will pull the convergence frequency back to 2π rad/s. From another perspective, this learning behavior has a good filtering effect on noise with uniformly distributed amplitudes.

3) The nature of disturbances during convergence

As mentioned earlier, the convergence of continuous spectrum is related to the position of the spectrum derivative. In computer numerical simulations, the input spectrum is not continuous and depends on the sampling frequency. The simulation experiment results are theoretically very unstable, but in fact, the simulation experiment results are very stable and effective.

Observing the simulated convergence process in Figure 15, the CPG frequency converges during oscillation. When using disturbance theory for estimation, there are sine signal components ${\omega }_1$ ${\tilde{\omega }}_2$ during the convergence process and steady state. These sines force the CPG to oscillate during steady state and convergence, which may cause the frequency component to cross the convergence critical point between the frequency components and towards the convergence of a stronger attraction frequency.

The amplitude of oscillation depends on two factors based on the previous derivation: 1) the amplitude distribution of the spectrum $F\left(t\right)$ ; 2) the disturbance gain $\epsilon$ . Based on the composition of Equation (38), in disturbance estimation, $A_{\omega }$ and a, as well as their corresponding frequencies, directly represent the amplitude of ${\omega }_1$ ${\tilde{\omega }}_2$ . Based on previous simulation results, the disturbance signals are changed to $\sin \left(\pi t\right)+\sin \left(2\pi t\right)$ and $\sin \left(2\pi t\right)$ respectively. The simulation convergence process and steady-state results are shown in Figure 18.

Figure 18. Convergence comparison with different disturbance signals.

The convergence result and convergence speed are determined by the conclusions of the previous section, while the convergence error after convergence increases due to the excess frequency components in the input signal. The gain $\epsilon$ represents the impact of the disturbance quantity $F\left(t\right)$ on the state quantity $\varphi ,\omega$ . Considering the disturbance signal as $\sin \left(\text{1}\text{.8}\pi t\right)+\sin \left(2\pi t\right)$ and the simulation results are shown in Figure 19. The convergence result oscillates back and forth between the two frequency components, causing the final convergence mean to be between the two and greatly increasing the convergence error of the result.

Figure 19. Convergence results with different disturbance coefficients.

5) Convergence characteristics of complex limit cycles

By replacing the limit cycle trajectory with Equation 16, the CPG oscillation signal of this limit cycle trajectory has continuous values on a sufficiently wide spectrum. Therefore, all odd order derivatives of the trajectory signal form a multidimensional orthogonal space with the original trajectory signal. As long as the number of derivative orders is large enough, the input signal $F\left(t\right)$ can be expanded in this space, and it can be expressed as a linear combination of the trajectory signal and multiple odd-order derivatives, similar to the process of Fourier transform. This allows the complex limit cycle trajectory equation to be treated as a single frequency unit, making it theoretically possible to achieve the same derivation and proof method as the circular trajectory.

On the other hand, starting again from E₃ and conducting disturbance theory derivation is equivalent to performing second-order integration on Equation (41), resulting in the following Equation (47).

${\displaystyle {\iint }_{\begin{array}{l}{\omega }_1,{\omega }_2\\ {\omega }_1\ne {\omega }_2\end{array}}\frac{A_{{\omega }_1}^2{\omega }_1^3}{2}}\left(-\frac{4A_0^2}{{\omega }_1}+\frac{2{\left|A_{{\omega }_2}\right|}^2{\omega }_1}{{\omega }_2^2-{\omega }_1^2}\text{d}{\omega }_1\text{d}{\omega }_2\right)$ (47)

Here, ${\omega }_1,{\omega }_2$ represents the frequency components of the limit cycle and the disturbance signal. The convergence direction becomes the superposition effect of the oscillator frequency component and the disturbance signal frequency component.

In disturbance theory derivation of complex limit cycle trajectories, only the term in Equation (39) is independent of the phase components of the limit cycle. Therefore, to prove that the convergence law is still mainly determined by Equation (39), we use the characteristic of no phase information in Equation (39) to conduct the following simulation:

Setting the disturbance signal as:

$F=\sin \left(\pi t\right)+\sin \left(2\pi t\right)+3\cos \left(4\pi t\right)+4\cos \left(5\pi t\right)+5\cos \left(6\pi t\right)$ (48)

The limit cycle trajectory is set to two groups:

$\begin{array}{l}f\left(\varphi \right)=\sin \left(p\right)+\sin \left(2p+\pi /6\right)+3\cos \left(4p+\pi /3\right)\\ +4\cos \left(5p+\pi /2\right)+5\cos \left(6p\right);\end{array}$

$\begin{array}{l}g\left(\varphi \right)=\cos \left(p\right)+2\cos \left(2p+\pi /6\right)-12\sin \left(4p+\pi /3\right)\\ -20\sin \left(5p+\pi /2\right)-30\cos \left(6p\right)\end{array}$

and

$f\left(\varphi \right)=\sin \left(p\right)+\sin \left(2p\right)+3\cos \left(4p\right)+4\cos \left(5p\right)+5\cos \left(6p\right);$

$g\left(\varphi \right)=\cos \left(p\right)+2\cos \left(2p\right)-12\sin \left(4p\right)-20\sin \left(5p\right)-30\cos \left(6p\right)$

The initial angular velocities are both 4π rad/s, and the disturbance coefficient is set to 0.01. The convergence results are shown in Figure 20, where the processes and convergence results of the two are basically the same. This indicates that the consistency of phase does not affect the convergence results at all, which is consistent with the predicted result of Equation (39) and the assumption. Of course, if the phase of the frequency component of the oscillator trajectory signal does not correspond to the phase of the frequency component of the disturbance signal, the final converged phase has no meaning at all.

3.3. CPG Algorithm Construction

The basic convergence property shows that the constant component of the disturbance signal has a strong influence on the frequency convergence of the oscillator. In the hip joint angle spectrum curve, the constant component occupies 50%

Figure 20. Convergence impact on phase consistency.

of the signal energy. Therefore, in order to ensure the stability of the convergence, the average component must be minimized as much as possible. Building feedback structure shown in Figure 21 and the corresponding convergence Equation (49) based on reference:

Figure 21. Algorithm structure.

$\begin{array}{l}\dot{\varphi }=\omega +\epsilon eg\left(\varphi \right)\\ \dot{\omega }=\epsilon eg\left(\varphi \right)\\ \dot{a}=k_{a}ef\left(\varphi \right)\\ {\dot{a}}_0=k_{a}e\\ e=F\left(t\right)-a_0-af\left(\varphi \right)\end{array}$ (49)

where 0 unit is used to learn the signal mean value and corresponds to the fourth line formula in Equation (49). The 1 unit is used to learn the phase of the external periodic signal and is constructed using the learning law in this chapter. Limit cycle $f\left(\varphi \right)$ is a custom trajectory. The corresponding amplitude is learned using the third line formula in Equation (49).

The sum of the external disturbance signal minus the output of the two units generates an error signal, which will drive the two error signals to continue learning the external signal until reaching a convergent stable state. At stable state, the 1 unit has synchronized with the external disturbance signal. Combining its stable angular frequency and current phase, it can obtain the future phase change for a period of time. This information can be transmitted to the lower-level control through the human motion information mapping.

4. Simulation and Discussion

4.1. Simulation Model Establishment and Stiffness Correction

The stiffness value of the VSA mechanism designed in Section 2 does not completely comply with the design rule in the actual working process, which has been pointed out in the stiffness design. Moreover, actual stiffness values are required to correct the variable stiffness compensation force in subsequent control. Therefore, it is necessary to establish a corresponding simulation model to obtain the actual stiffness curve, verify the design results, and conduct subsequent control simulations.

4.1.1. Simulation Model Establishment

Adams is an authoritative software for dynamic simulation, which can effectively simulate a large number of mechanical systems. The VSA mechanism includes plate springs as elastic components. To make the plate springs perform realistically in the Adams model, finite element analysis is needed.

Although Adams internally includes Autoflex, the functions are relatively basic. For this VSA mechanism, the boundary elements and plate spring mounting boundary conditions need to be defined since Autoflex is unable to perform well. Also, there is a lack of reasonable subdivision mesh units for stress deformation. Therefore, it is necessary to define the Modal Neutral File for the plate spring flexible body in an external model and import it. This paper uses Abaqus 6.14 for frequency analysis of the plate spring and imports it into Adams.

Firstly, the material density, elastic modulus, and Poisson ratio of the plate spring are set as 7.85e⁻⁹ (tonne/mm³), 1.97e⁵ (MPa), and 0.3, respectively, in Abaqus. Since the plate spring does not work under extreme conditions, parameters such as yield limit need not be set. Standard linear beam elements are used.

Considering the force curve of the plate spring head in the clamping method, the rigid part that has little effect on the elasticity is removed and replaced with a parabolic shape. The truncated parabolic section is connected to a reference point Rp in the form of a beam element using the MPC multisite constraint method, as shown in Figure 22. At the same time, Rp is taken as the load point, and its 6 degrees of freedom are constrained. The modal analysis is performed in six directions for the RP reference point, and then the Rp point is defined as the sub-model interface marker in Adams. The corresponding mesh is generated with a scale of 1.6.

Plate spring frequencies are extracted every 1 to 15 orders, and the frequency analysis results are used to generate sub-model data, export the mass matrix and stress-strain data, and output the corresponding mnf file according to the unit output described above. The first-order mode extraction result of the plate spring frequency is shown in Figure 23. The first-order frequency of 62Hz indicates that the plate spring has almost no dynamic characteristics in the current application environment, which is consistent with the characteristics of the low mass of the plate spring itself.

In Adams, Rp is defined to be fixedly connected to the mounting seat, and the plate spring is defined to be in contact with the two-side needle roller bearings by a combination of flexible and rigid bodies. In order to prevent the occurrence of piercing during the simulation process, the stiffness of the needle roller bearing is reduced to 1e⁺⁶ (N/m), and the penetration depth is limited to 1e⁻³ (m). The final result of simulation is shown in Figure 24.

Figure 22. Grid partition and MPC constraint settings.

Figure 23. First-order mode of the plate spring.

Figure 24. Schematic diagram of the deformation of the flexible body.

4.1.2. Stiffness Experiment

Stiffness testing was conducted on the VSA mechanism based on the Adams model set up mentioned above. The test was divided into two parts: the one-way loading experiment and the reciprocating loading experiment. The former was mainly used to obtain the stiffness variation trend, stiffness variation range and the reason for the non-linear variation of the stiffness, while the latter aimed to explore the existence of dead zones and hysteresis in the displacement-torque curve during reciprocating loading to obtain a stiffness fitting curve.

1) One-way loading experiment

a) Simulation experiment settings

The experimental device was initially in the form of a square structure composed of four linkages, and the differential angle of the motor was pi/2. In the first two seconds, the two motors were operated in opposite directions to pre-rotate the same angle. When the angle was positive, the contact point of the spring plate was rotated around the center to increase the stiffness, and when the angle was negative, it was decreased. Afterwards, the two motors were fixed to ensure that the position of the slider no longer changed. The load end was uniformly increased from 0 N·m to a maximum of 20 N·m using a ramp function from the 2nd second until the 12th second.

A total of 11 groups of experiments were conducted, and each group achieved different stiffness by adjusting the pre-rotation angles of the two motors. Considering that the four-bar mechanism has a singular position, the pre-rotation angles were set from -0.2rad to 0.8rad with an interval of 0.1rad from low to high. The experimental results are shown in Figure 25.

b) Analysis of experimental results

Two general trends can be observed from Figure 25:

1) The displacement-torque curve becomes steeper as the pre-rotation angle increases, which indicates that the stiffness value is increasing. The displacement at a load of 20 N·m for the curve with the maximum stiffness value is 2.5 times that of the curve with the minimum stiffness value. ii. Under the same pre-rotation angle, as the load torque increases, the increase in displacement decreases and stiffness increases. For example, the slope representing stiffness at point P1 is significantly higher than that at point P0.

In order to obtain the specific stiffness values, a third-order polynomial fitting was performed on the curves in Figure 25. Taking the curve with the most irregular shape and most severe pre-rotation, −0.2 rad, as an example, the fitting result is shown in Figure 26. The residual value of 8.1 indicates that the fitting result is relatively good. The stiffness values were derived based on the fitting curve by taking derivatives, as shown in the stiffness curve in Figure 27, which also conforms to the two trends in Figure 25. The following two conclusions can be drawn: 1) The range of stiffness values varies from the lowest 28 N·m/rad to 275 N·m/rad. The lowest stiffness occurs at a pre-rotation angle of −0.2 rad with almost no load, and the highest stiffness occurs at a pre-rotation of 0.8 rad with a load of 20 N·m. 2) Under the same pre-rotation angle, the trend of stiffness increasing with increasing load torque increases as the pre-rotation angle increases. From the parallel trends between the curves, it can be seen that the stiffness curve of pre-rotation 0.8 rad increases much faster with increasing load than the stiffness curve of pre-rotation −0.2 rad.

The characteristic of the stiffness curve increasing with increasing load torque can be explained by Figure 28 and Figure 29. In the working process, the plate spring is not subject to a single-sided force but is subject to pressure on both sides of the needle roller bearing. When not subjected to torque load, the plate spring is in contact with the needle roller bearings on both sides at the same position, and the stiffness value at this time is close to the design value. As shown in Table 2, where the experimental stiffness represents the stiffness of rotating 0 rad in Figure 22, which is the stiffness under no load. The results show that the design results and experimental results are closer in the low-stiffness range. When subjected to torque load, the contact points on both sides of the plate spring do not have the same length. The two different positions of the needle roller bearings result in a situation similar to a lever for the plate spring. The contact point away from the center of rotation serves as a fulcrum, while the other contact point serves as a force application point, causing the actual working length of the plate spring to become shorter and the stiffness to increase. The influence of the working length of the plate spring on stiffness depends on the current working length, the shorter the working length, the more dramatic the increase in stiffness, which is also the reason for the second point mentioned above.

Figure 25. Displacement curve of unidirectional 0 - 20 N·m uniform loading.

Figure 26. Third-order polynomial fitting of the pre-rotation −0.2 rad displacement-torque curve.

Figure 27. Displacement curve of unidirectional loading.

Figure 28. Pre-rotation 0.8 rad.

Figure 29. Pre-rotation −0.3 rad.

When a load is applied, the lever effect causes the forces on the two needle roller bearings to be unequal. The torques of the two contact pairs on the rotating center axis were monitored as shown in Figure 30. It can be clearly seen that the

Table 2. Comparison between design stiffness and experimental stiffness.

Pre-rotation (rad)	−0.2	−0.1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
Design stiffness (N·m/rad)	30	39	49	61	76	93	113	135	160	187	216
Experimental stiffness (N·m/rad)	26	36	47	58	67	78	87	95	105	115	128

sum of the two torques multiplied by 4 is almost equal to the load-bearing force of 20 N·m. This also indicates that the force on the plate spring is not the same as the ideal design.

Figure 30. Torque of −0.2 rad pair of needle roller bearings on rotating shaft.

2) Reciprocating loading experiment

a) Simulation experiment settings

The simulation process is similar to the one-way loading experiment, except that after 12 seconds, the torque is maintained for 2 seconds, and then the load torque changes from 20 N·m to −20 N·m using a ramp function from 14 seconds to 24 seconds. The torque is then maintained for 2 seconds, and finally the load changes from −20 N·m to 0 N·m from 26 seconds to 36 seconds. The experimental groups were still divided into 11 groups according to the pre-rotation angle from −0.2 rad to 0.8 rad. The experimental results are shown in Figure 31.

b) Analysis of experimental results

The following conclusions can be drawn from the data shown in Figure 31:

1) The mechanism has extremely small creep characteristics, and the creep decreases as the pre-rotation angle increases. At the maximum load point in Figure 31, after the slow loading torque stops increasing and is maintained for 2 seconds, the mechanism still undergoes a small displacement. The displacement decreases gradually from low pre-rotation angles to high pre-rotation angles. In addition, at the pre-rotation angle of −0.2 rad, the creep is only 0.002 rad when the maximum force is maintained, indicating that the creep is minimal. The above statement also implies that creep under low stiffness conditions is more severe than that under high stiffness conditions.

2) The mechanism exhibits a certain hysteresis phenomenon, which is related to the pre-rotation angle and stiffness value. The hysteresis curves in Figure 23 all exhibit loops, and the area of the loops decreases as the pre-rotation angle increases. At −0.2 rad, the highest angle displacement difference can be as large as 0.02 rad (approximately 1 degree) under the same torque.

3) The overall mechanism exhibits a work characteristic of central symmetry, which conforms to the original intention of symmetric structural design.

4) The mechanism as a whole has an extremely small zero error, which can be ignored. Figure 32 shows the curve near zero position of Figure 31, and it can be seen that the maximum zero error is 0.01 rad.

Figure 31. Timing diagram for simulation testing.

Figure 32. Load displacement curves under different stiffness settings.

4.1.3. Stiffness Correction and Working Range

Based on the experimental results mentioned above, it is necessary to simplify the experimental data in the reciprocating experiment to modify the calculation of the torsional stiffness $K_T=f_K\left(T_3,r\right)$ . Since the radius r is a single-value function of the differential angle of the two motors, in order to facilitate calculation and reflect the actual stiffness, corresponding numerical tables and tables of differential angle values are established according to the experimental results.

Firstly, polynomial estimation is applied to the results of the reciprocating experiment in Figure 31, taking the pre-rotation angle of −0.2 rad with the most complex curve as an example. As shown in Figure 33, the fitting curve and experimental data curve fit well, with a residual value of only 62. Then, according to the polynomial fitting results, the slope-stiffness curve is derived from the curve and shown in Figure 34. The torque-stiffness curve is also shown in Figure 35 according to the fitting formula. Figure 34 mainly serves to calculate the current stiffness based on the current rotation angle (equivalent to knowing the difference between the output angle of the two motors and the output angle), and the differential angle of the two motors (equivalent to pre-rotation angle) in subsequent control. Figure 35 mainly serves to calculate the current differential angle of the motors based on the target torque and stiffness in subsequent control. The region A surrounded by the curves in Figure 34 and Figure 35 essentially indicates the working range of the VSA, and any action beyond region A cannot be successfully achieved. Therefore, boundary penalty functions can also be generated based on these two graphs to prevent loss of control.

Figure 33. Load-displacement-torque curve under pre-rotation −0.2 rad.

Figure 34. Reciprocating load angle-stiffness curve.

Figure 35. Reciprocating load torque-stiffness curve.

4.2. Simulation of CPG Phase Recognition Algorithm

To verify the performance of the CPG algorithm constructed in this paper, EMG-based dataset of human walking on a treadmill is collected, and output signals of CPG and real signals are compared. The EMG-based dataset was constructed by our laboratory, with approval from the Chinese Ethics Committee of Registering Clinical Trials. Nine quadrupolar EMG electrodes were mounted on lower-limb muscles (rectus femoris, vastus lateralis, vastus medialis, tibialis anterior, soleus, biceps femoris, semitendinosus, gastrocnemius medial head, and gastrocnemius lateral head) with a sampling frequency of 1111.11 Hz. According to Equation (42), the high-frequency noise contained in the limit cycle will bring additional interference attraction. Therefore, first use a 3rd-order Butterworth filter to filter out noise above 1 Hz from the joint angle signal and angular velocity signal of human walking on a treadmill. Then select the waveform corresponding to one cycle as the and trajectories $f\left(\varphi \right)$ and $g\left(\varphi \right)$ for the 1 unit, and normalize the time.

Set the disturbance coefficient $\epsilon =0.1$ , the amplitude learning gain $k_a=1$ for the 1 unit and 0 unit, and the initial value were set to 20˚ and 10˚ respectively, with an initial phase of 0.5. The experimental results are shown in Figure 36 and Figure 37. It can be clearly seen from the figures that learning can be completed in three steps with significant results. It is foreseeable that the predictive information provided to the lower-level control is quite accurate.

Figure 36. Angle signal learning diagram.

Figure 37. Comparison of normal walking hip joint signals and learning signals.

5. Conclusions

This paper aims to conduct experiments on human-machine interaction variable stiffness and complete the following works:

1) Analyzed the requirements for variable stiffness in human-machine interaction experiments, used plate springs as elastic elements, and four-bar linkages as driving mechanisms to complete the structural design of a small variable stiffness actuator. In the analysis of the dynamic model, it was found that the non-linearity of the mechanism is mainly caused by the non-linear stiffness of the plate spring. Then, Abaqus software was used for finite element analysis of the plate spring to determine that the dynamic characteristics of the plate spring can be ignored. Using the Abaqus software interface to generate a corresponding simulation model in combination with Adams, a unidirectional force and reciprocating force were applied to the variable stiffness actuator under joint simulation in MATLAB and Adams. It was found that the displacement curve of the stiffness of the mechanism was centrally symmetrical and that dead zones and hysteresis phenomena were not significant. Finally, the experimental data was fitted by a third-order polynomial to obtain the corresponding stiffness table.

2) Addressing the lack of research on the convergence characteristics of CPG in human body phase recognition. First, the perturbation theory was used to analyze the convergence characteristics of CPG units with the same frequency and concluded that the hip joint angle signal will converge in phase under the same frequency. Then, the perturbation theory was used to second-order expand the CPG unit, and a judgment formula for frequency convergence was obtained when the limit cycle was a circle. According to the formula, corresponding convergence criteria were obtained and the conclusion was extended to the non-circular case of the limit cycle. The conclusion was that the convergence of the input signal with a negative frequency in the load limit cycle is not affected. Finally, based on the conclusions and criteria, a CPG feedback loop was constructed to eliminate the mean, and the Simulink was used to verify the algorithm’s reliability and convergence.

There are still some shortcomings in this work. Firstly, CPG only optimizes the upper-level control of the hip joint exoskeleton. Secondly, the simulation experiments of the CPG phase recognition algorithm simulate relatively simple types of motion. It is necessary to consider actual application requirements and observe the performance of the algorithm under different motion modes and control accuracy requirements. Finally, the research results of this paper still need further real-machine verification to prove. Future works will focus on improving the lower-level control of the exoskeleton to realize more accurate and stable control of the exoskeleton, especially the feedforward compensation considering the nonlinear stiffness of the VSA. Simulations containing lower-level control should be conducted, and more research and experiments should be taken on the real machine to verify the performance of the whole control system.

Acknowledgments

This work is partly funded by the Technique Program of Jiangsu (No.BE2021086).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bao, G.J., et al. (2019) Academic Review and Perspectives on Robotic Exoskeletons. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 27, 2294-2304. https://doi.org/10.1109/TNSRE.2019.2944655
[2]	Wesley, R., et al. (2017) On the Stiffness Selection for Torque-Controlled Series-Elastic Actuators. IEEE Robotics and Automation Letters, 2, 2255-2262. https://doi.org/10.1109/LRA.2017.2726141
[3]	Zhang, T., Tran, M. and Huang, H. (2019) Admittance Shaping-Based Assistive Control of SEA-Driven Robotic Hip Exoskeleton. IEEE/ASME Transactions on Mechatronics, 24, 1508-1519. https://doi.org/10.1109/TMECH.2019.2916546
[4]	Baud, R., Manzoori, A.R., Ijspeert, A., et al. (2021) Review of Control Strategies for Lower-Limb Exoskeletons to Assist Gait. Journal of Neuro Engineering and Rehabilitation, 18, Article No. 119. https://doi.org/10.1186/s12984-021-00906-3
[5]	Grebenstein, M., Chalon, M., Friedl, W., Haddadin, S., Wimböck, T., Hirzinger, G. and Siegwart, R. (2012) The Hand of the DLR Hand Arm System: Designed for Interaction. The International Journal of Robotics Research, 31, 1531-1555. https://doi.org/10.1177/0278364912459209
[6]	Braun, D.J., Petit, F., Huber, F., Haddadin, S., van der Smagt, P., Albu-Schaffer, A. and Vijayakumar, S. (2013) Robots Driven by Compliant Actuators: Optimal Control under Actuation Constraints. IEEE Transactions on Robotics, 2, 1085-1101. https://doi.org/10.1109/TRO.2013.2271099
[7]	Visser, L.C., Stramigioli, S. and Bicchi, A. (2011) Embodying Desired Behavior in Variable Stiffness Actuators. IFAC Proceedings, 44, 9733-9738. https://doi.org/10.3182/20110828-6-IT-1002.02202
[8]	Migliore, S.A., Brown, E.A. and de Weerth, S.P. (2005) Biologically Inspired Joint Stiffness Control. Proceedings of the 2005 IEEE International Conference on Robotics and Automation, Barcelona, 18-22 April 2005.
[9]	Jafari, A., Tsagarakis, N.G., Vanderborght, B. and Caldwell, D.G. (2010) A Novel Actuator with Adjustable Stiffness (AwAS). 2010 IEEE/RSJ International Conference on Intelligent Robots and Systems, 18-22 October 2010, 4201-4206. https://doi.org/10.1109/iros.2010.5648902
[10]	Jafari, A., Tsagarakis, N.G. and Caldwell, D.G. (2011) AwAS-II: A New Actuator with Adjustable Stiffness Based on the Novel Principle of Adaptable Pivot Point and Variable Lever Ratio. 2011 IEEE International Conference on Robotics and Automation, Shanghai, 9-13 May 2011, 4638-4643. https://doi.org/10.1109/icra.2011.5979994
[11]	Wolf, S., Eiberger, O. and Hirzinger, G. (2011) The DLR FSJ: Energy Based Design of a Variable Stiffness Joint. 2011 IEEE International Conference on Robotics and Automation, Shanghai, 9-13 May 2011, 5082-5089. https://doi.org/10.1109/icra.2011.5980303
[12]	Eiberger, O., Haddadin, S., Weis, M., Albu-Schäffer, A. and Hirzinger, G. (2010) On Joint Design with Intrinsic Variable Compliance: Derivation of the DLR QA-Joint. 2010 IEEE International Conference on Robotics and Automation, Anchorage, 3-7 May 2010, 1687-1694. https://doi.org/10.1109/robot.2010.5509662
[13]	Morita, T. and Sugano, S. (1995) Development of One-DOF Robot Arm Equipped with Mechanical Impedance Adjuster. Proceedings 1995 IEEE/RSJ International Conference on Intelligent Robots and Systems. Human Robot Interaction and Cooperative Robots, Pittsburgh, 5-9 August 1995, 407-412. https://doi.org/10.1109/iros.1995.525828
[14]	Hollander, K.W., Sugar, T.G. and Herring, D.E. (2005) Adjustable Robotic Tendon Using a ‘Jack Spring’. 2005 9th International Conference on Rehabilitation Robotics, Chicago, 28 June-1 July 2005, 113-118.
[15]	Vanderborght, B., Albu-Schaeffer, A., Bicchi, A., Burdet, E., Caldwell, D.G., Carloni, R., et al. (2013) Variable Impedance Actuators: A Review. Robotics and Autonomous Systems, 61, 1601-1614. https://doi.org/10.1016/j.robot.2013.06.009
[16]	Hiroaki, S., et al. (2000) Development of a Robot Joint Mechanism with Variable Compliance by Rotating a Leaf Spring. https://searchworks-lb.stanford.edu/view/4554798
[17]	Tucker, M.R., Olivier, J., Pagel, A., Bleuler, H., Bouri, M., Lambercy, O., et al. (2015) Control Strategies for Active Lower Extremity Prosthetics and Orthotics: A Review. Journal of Neuro Engineering and Rehabilitation, 12, Article No. 1. https://doi.org/10.1186/1743-0003-12-1
[18]	Guo, W., Yang, C.W., et al. (2015) Dynamic Solution of the Lower Extremity Joint Torques in Man-Machine System of Lower Extremity Exoskeleton. Machinery & Electronics, 10, 71-75.
[19]	Farina, D., Merletti, R. and Enoka, R.M. (2004) The Extraction of Neural Strategies from the Surface EMG. Journal of Applied Physiology, 96, 1486-1495. https://doi.org/10.1152/japplphysiol.01070.2003
[20]	Ardestani, M.M., Zhang, X., Wang, L., Lian, Q., Liu, Y., He, J., et al. (2014) Human Lower Extremity Joint Moment Prediction: A Wavelet Neural Network Approach. Expert Systems with Applications, 41, 4422-4433. https://doi.org/10.1016/j.eswa.2013.11.003
[21]	Sepulveda, F., Wells, D.M. and Vaughan, C.L. (1993) A Neural Network Representation of Electromyography and Joint Dynamics in Human Gait. Journal of Biomechanics, 26, 101-109. https://doi.org/10.1016/0021-9290(93)90041-c
[22]	Karatsidis, A., Jung, M., Schepers, H.M., et al. (2018) Predicting Kinetics Using Musculoskeletal Modeling and Inertial Motion Capture. https://doi.org/10.48550/arXiv.1801.01668
[23]	Xiong, B., Zeng, N., Li, H., Yang, Y., Li, Y., Huang, M., et al. (2019) Intelligent Prediction of Human Lower Extremity Joint Moment: An Artificial Neural Network Approach. IEEE Access, 7, 29973-29980. https://doi.org/10.1109/access.2019.2900591
[24]	Goulermas, J.Y., Howard, D., Nester, C.J., Jones, R.K. and Ren, L. (2005) Regression Techniques for the Prediction of Lower Limb Kinematics. Journal of Biomechanical Engineering, 127, 1020-1024. https://doi.org/10.1115/1.2049328
[25]	Zeng, N., Wang, Z., Zineddin, B., Li, Y., Du, M., Xiao, L., et al. (2014) Image-Based Quantitative Analysis of Gold Immunochromatographic Strip via Cellular Neural Network Approach. IEEE Transactions on Medical Imaging, 33, 1129-1136. https://doi.org/10.1109/tmi.2014.2305394
[26]	Souron, R., Bordat, F., Farabet, A., Belli, A., Feasson, L., Nordez, A., et al. (2016) Sex Differences in Active Tibialis Anterior Stiffness Evaluated Using Supersonic Shear Imaging. Journal of Biomechanics, 49, 3534-3537. https://doi.org/10.1016/j.jbiomech.2016.08.008
[27]	Jang, J., Kim, K., Lee, J., Lim, B. and Shim, Y. (2015) Online Gait Task Recognition Algorithm for Hip Exoskeleton. 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems, Hamburg, 28 September-2 October 2015, 5327-5332. https://doi.org/10.1109/iros.2015.7354129
[28]	Zhang, T., Tran, M. and Huang, H.H. (2017) NREL-Exo: A 4-DoFs Wearable Hip Exoskeleton for Walking and Balance Assistance in Locomotion. 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems, Vancouver, 24-28 September 2017, 508-513. https://doi.org/10.1109/iros.2017.8202201
[29]	Schrade, S.O., Dätwyler, K., Stücheli, M., Studer, K., Türk, D., Meboldt, M., et al. (2018) Development of VariLeg, an Exoskeleton with Variable Stiffness Actuation: First Results and User Evaluation from the CYBATHLON 2016. Journal of Neuro Engineering and Rehabilitation, 15, Article No. 18. https://doi.org/10.1186/s12984-018-0360-4
[30]	Ding, Y., Galiana, I., Siviy, C., Panizzolo, F.A. and Walsh, C. (2016) Imu-Based Iterative Control for Hip Extension Assistance with a Soft Exosuit. 2016 IEEE International Conference on Robotics and Automation, Stockholm, 16-21 May 2016, 3501-3508. https://doi.org/10.1109/icra.2016.7487530
[31]	Bergmann, L., Lück, O., Voss, D., Buschermöhle, P., Liu, L., Leonhardt, S., et al. (2023) Lower Limb Exoskeleton with Compliant Actuators: Design, Modeling, and Human Torque Estimation. IEEE/ASME Transactions on Mechatronics, 28, 758-769. https://doi.org/10.1109/tmech.2022.3206530
[32]	Lee, S., Crea, S., Malcolm, P., Galiana, I., Asbeck, A. and Walsh, C. (2016) Controlling Negative and Positive Power at the Ankle with a Soft Exosuit. 2016 IEEE International Conference on Robotics and Automation, Stockholm, 16-21 May 2016, 3509-3515. https://doi.org/10.1109/icra.2016.7487531
[33]	Chinimilli, P.T., Subramanian, S.C., Redkar, S. and Sugar, T. (2019) Human Locomotion Assistance Using Two-Dimensional Features Based Adaptive Oscillator. 2019 Wearable Robotics Association Conference (WearRAcon), Scottsdale, 25-27 March 2019, 92-98. https://doi.org/10.1109/wearracon.2019.8719628
[34]	Seo, K., Lee, J., Lee, Y., Ha, T. and Shim, Y. (2016) Fully Autonomous Hip Exoskeleton Saves Metabolic Cost of Walking. 2016 IEEE International Conference on Robotics and Automation, Stockholm, 16-21 May 2016, 4628-4635. https://doi.org/10.1109/icra.2016.7487663
[35]	Seo, K., Hyung, S.Y., Choi, B.K., Lee, Y. and Shim, Y. (2015) A New Adaptive Frequency Oscillator for Gait Assistance. 2015 IEEE International Conference on Robotics and Automation, Seattle, 26-30 May 2015, 5565-5571. https://doi.org/10.1109/icra.2015.7139977