Introduction

Stroke and spinal cord injury are the two leading causes of damage to the central nervous system, often leading to hemiplegia or paraplegia1. Stroke is a clinical manifestation of acute cerebrovascular disease characterized by high morbidity and disability2. Hemiplegia can significantly reduce an individual’s mobility, usually resulting in altered metabolic function and body composition3,4. A damage to the spinal cord can cause either whole or partial paralysis, affecting neuronal signaling and regulation both above and below the site of injury. This can lead to impairments in autonomic, sensory, neuromusculoskeletal, and motor functions5. Functional electrical stimulation (FES) is a technique that uses low-frequency electrical pulses to artificially activate paralyzed muscles to produce movements such as walking, grasping and standing6,7. FES seeks to assist these individuals in improving or restoring the functional ability of their paralyzed muscles8. FES-assisted rehabilitation programs have many benefits, such as pain relief, cardiopulmonary recovery, bone development and metabolic function enhancement9,10,11.

FES-assisted cycling is a rehabilitation method that enables people with limb weakness to train on tricycles or stationary cycling trainers12. In general, FES-assisted cycling promotes coordinated alternating contractions of specific lower limb muscles in a sequence that varies according to the speed and angle of the ride13. There are some FES pedal efficiency issues for traditional FES cycling stimulation mode. Trial-and-error tests are a typical way to find each patient’s stimulation pattern during cyclic therapy14,15. But the physiotherapist’s experience determine show effective the stimulation pattern is. Not only is it a waste of the patient’s recovery time, but it also places a burden on the patient and the healthcare provider.

The closed-loop control strategy of FES pedal movement has the potential to provide services with high intensity, low staffing costs and good control accuracy16,17,18,19,20. Griffis et al.21proposed a closed-loop adaptive based on Deep Neural Networks control scheme for FES cycling to ensure the asymptotic convergence of tracking errors. However, this control method has not been extended to FES control inputs. In previous studies that adopted adaptive closed-loop control schemes based on proportional-integral-derivative (PID) control algorithms, the fractional-order PID controller proposed by Podlubny demonstrated excellent performance in handling complex dynamic systems22. In FES systems, fractional-order controllers have also been attempted to improve control accuracy and system robustness23. On the other hand, Wen et al. proposed an observer-based PID controller that exhibited stable performance in controlling random errors and noise influences24. For fractional-order controllers, their design and implementation are relatively complex, requiring in-depth mathematical background and significant computational resources, which may be limited in real-time control. Observer-based methods demand high accuracy of the system model; model uncertainties and parameter variations may affect the performance of the observer. However, there are still several problems in the current research phase of using electrical stimulators to help patients with bicycle training. Firstly, Traditional stimulation modes are not able to maintain a steady training speed for long periods of time. Secondly, some studies have used real-time electromyography (EMG) to drive FES system, but the FES generated artifact contaminated EMG signals25,26. Thirdly, the PID control algorithm for the adaptive closed-loop control scheme has defects and,has the disadvantages of easy to cause local optimal solutions and slow convergence speed.

Hence, in this study, IMU sensors were used to obtain data to avoid the artifacts generated by FES from polluting the EMG signals, and the stimulation pattern was obtained from the motor rehabilitation of incomplete spinal cord injury patients (T9-T12) using pedaling motion hill modeling and optimizing the target with maximum moment efficiency. We proposed a PSO-BP-PID controller. It has the adaptive ability to adjust control parameters in real-time according to system dynamics, adapting to individual differences among patients. The particle swarm optimization algorithm improves convergence speed, avoiding the problem of traditional PID algorithms easily falling into local optima, enhances computational efficiency, and is easy to deploy and implement in practical systems.

Methods

Establishment of mathematical model

In order to study in depth, the principles of human-system interaction during stirrups movement of patients, a dynamic model of the system was developed. The model consists of a lower limb dynamic model for cycling and a dynamic muscle model.

Lower limb dynamic model for cycling

In the process of FES cycling training, joint torque is generated by muscle contraction to drive the relative movement of each part of the lower extremity. The study demonstrates that electrical stimulation activates muscle fibers, leading to contractions muscle joint angles and angular velocities in knee joint movements. Therefore, joint torque can be directly obtained through the activation of muscle fibers and joint angle. Meanwhile, the transfer of joint torque to crank torque is determined by the position of lower limb movement. To get the right stimulation pattern need to know how the lower limb positions change during the pedaling process. Physically, the pedaling motion can be regarded as a single degree of freedom system, and the kinematic model of pedaling can be established, and then the relative positions of lower limb parts can be obtained through the crank angle.

Fig. 1
figure 1

(a) Schematic diagram of FES Cycling training. (b) Schematic diagram of unilateral pedal motion structure.

Full size image

The following will select the left lower limb movement for analysis. As shown in part (b) of Figure 1, \(L_{t}\), \(L_{s}\) and \(L_{cr}\) represent the thigh, calf and crank length(cm) respectively. The angle between the horizontal line and the thigh, calf, and crank are denoted by \(\theta _{1}\), \(\theta _{2}\)and \(\theta _{3}\), respectively. The hip joint ___location served as the origin of the coordinate system, while the crank’s center was used as (\(x_c\),\(y_c\)). Set up a system of equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} \textrm{L}_\textrm{t}\cos \theta _1+\textrm{L}_\textrm{s}\cos \theta _2+\textrm{L}_\mathrm {{cr}}\cos \theta _3=\textrm{x}_\textrm{c}\\ \textrm{L}_\textrm{t}\sin \theta _1+\textrm{L}_\textrm{s}\sin \theta _2+\textrm{L}_\mathrm {{cr}}\sin \theta _3=\textrm{y}_\textrm{c}& \end{array}\right. }\end{aligned}$$
(1)

As shown in part (b) of Figure 1, the relationship between \(\theta _{3}\) and \(\theta _{cr}\) is \(\theta _3=-\theta _{\textrm{cr}}\), substituted in formula (1) to get:

$$\begin{aligned} {\left\{ \begin{array}{ll}\textrm{L}_\textrm{t}\cos \theta _1+\textrm{L}_\textrm{s}\cos \theta _2+\textrm{L}_\mathrm {{cr}}\cos \theta _{cr}=\textrm{x}_\textrm{c}\\ \textrm{L}_\textrm{t}\sin \theta _1+\textrm{L}_\textrm{s}\sin \theta _2-\textrm{L}_\mathrm {{cr}}\sin \theta _{cr}=\textrm{y}_\textrm{c}& \end{array}\right. }\end{aligned}$$
(2)

Suppose \(a=\frac{L_t}{L_s}\), \(b=\frac{y_c+L_{cr}\sin \theta _{cr}}{L_s}\), \(c=\frac{x_c-L_{cr}\sin \theta _{cr}}{L_s}\), substituted in formula (2) to get \({\left\{ \begin{array}{ll}\cos \theta _2=c-a\cos \theta _1\\ \sin \theta _2=b-a\sin \theta _1\end{array}\right. }\). According to the trigonometric formula:

$$\begin{aligned} & l=b^2-2ab\sin \theta _1+a^2\sin ^2\theta _1+c^2-2ac\cos \theta _1+a^2\cos ^2\theta _1 \end{aligned}$$
(3)
$$\begin{aligned} & \theta _{1}={\left\{ \begin{array}{ll}\pi -\arcsin \biggl (\frac{a^2+b^2+c^2-1}{2a\sqrt{b^2+c^2}}\biggr )-\arctan \frac{c}{b}& \theta _{cr}\in [0,\pi )\\[2ex]-\arcsin \biggl (\frac{a^2+b^2+c^2-1}{2a\sqrt{b^2+c^2}}\biggr )-\arctan \frac{c}{b}& \theta _{\text {cr}}\in [\pi ,2\pi )\end{array}\right. }\end{aligned}$$
(4)
$$\begin{aligned} & \theta _1-\theta _2=\arccos \left( \frac{b^2+c^2-a^2-1}{2a}\right) \end{aligned}$$
(5)

Through the above calculation, the relation of each bar angle is obtained. Part (b) of figure 1 shows the formula \(\theta _h=\theta _1\). According to the above functions, \(\theta _h\) and \(\theta _k\) are divided into functions of \(\theta _cr\), The meanings of \(\theta _h\) and \(\theta _k\) are: \(\theta _\textrm{h}=\mathrm {f_h}(\theta _\textrm{cr})\), \(\theta _\textrm{k}=\mathrm {f_k}(\theta _\textrm{cr})\). Taking the time derivative of formula (2), we get the angular velocity \(W_{1}\) and \(W_{2}\) of \(\theta _{1}\) and \(\theta _{2}\), \(\mathrm {w_1=-\frac{L_{cr}w\sin (\theta +\theta _1)}{L_t\sin (\theta _1-\theta _2)}}\), \(\mathrm {w_2=-\frac{L_{cr}w\sin (\theta +\theta _2)}{L_t\sin (\theta _2-\theta _1)}}\). where w is the angular velocity of the crank. Taking the second derivative of formula (2) in time, we get the acceleration \(a_{1}\) and \(a_{1}\) of \(\theta _{1}\) and \(\theta _{2}\), where a is the acceleration of the crank.

$$\begin{aligned} & \textrm{a}_\textrm{1}=\frac{\theta _1^\textrm{2}\textrm{L}_t\cos (\theta _1+\theta _2)+\textrm{w}_2^\textrm{2}\textrm{L}_\textrm{s}+\textrm{a}\textrm{L}_{\textrm{cr}}\cos (\theta -\theta _2)}{\textrm{L}_\textrm{t}\sin (\theta _1-\theta _2)} \end{aligned}$$
(6)
$$\begin{aligned} & \textrm{a}_\textrm{2}=\frac{\theta _2^\textrm{2}\textrm{L}_t\cos (\theta _1+\theta _2)+\textrm{w}_1^\textrm{2}\textrm{L}_\textrm{s}+\textrm{a}\textrm{L}_\mathrm {{cr}}\cos (\theta -\theta _1)}{L_t\sin (\theta _2-\theta _1)} \end{aligned}$$
(7)

The crank torque is calculated by Lagrange method and virtual work principle. First calculate the kinetic energy of the left leg \(T_{l}\), The kinetic energy of the leg is composed of moving kinetic energy and rotational kinetic energy.

$$\begin{aligned} T_l=\frac{1}{2}\textrm{m}_\textrm{t}\textrm{L}_\textrm{t}^\textrm{2}\theta _h^2+\frac{1}{2}\textrm{I}_{s}\theta _h^2+\frac{1}{2}m_s\left[ \textrm{L}_\textrm{t}^\textrm{2}\theta _h^2+\textrm{2}\textrm{L}_\textrm{t}\textrm{L}_\textrm{s}\theta _h\theta _2\cos (\theta _h-\theta _2)+\textrm{L}_\textrm{s}^2\theta _2^2\right] +\frac{1}{2}\textrm{I}_\textrm{s}\theta _2^2\end{aligned}$$
(8)

Express \(T_{l}\) relationship with \(\theta _{cr}\) as a function of \(T_l=T(\theta _{cr}){\theta _{cr}}^2\). The gravitational potential energy is formulated as follows:

$$\begin{aligned} T_l=\mathrm {L_1}\sin \theta _\textrm{h}\times \mathrm {m_tg}+(\mathrm {L_1}\sin \theta _\textrm{h}+\mathrm {L_2}\sin \theta _2)\times \mathrm {m_sg}\end{aligned}$$
(9)

Express \(V_{l}\) relationship with \(\theta _{cr}\) as a function of \(V_l=V(\theta _{cr})\). According to the dynamic analysis of the right leg lags \(\pi\) bar relative to the left leg \(T_{r}\), can be obtained \(T_r=T(\theta _{cr}-\pi )\theta _{cr}^2\) , \(V_{r}=V(\theta _{cr})\). The generalized force simplified by the Lagrangian formula is shown below27 :

$$\begin{aligned} \textrm{Q}_\textrm{cr}=\textrm{2}[\textrm{T}(\theta _{cr})\theta _{cr}+\textrm{T}(\theta _{cr}-\pi )\theta _{cr}]+\frac{\textrm{dT}(\theta _{cr})}{\textrm{dt}}\theta _{cr}+\frac{\textrm{dT}(\theta _{cr}-\pi )}{\textrm{dt}}\theta _{cr}+\frac{\textrm{dV}(\theta _{cr})}{\textrm{d}\theta _{cr}}+\frac{\textrm{dV}(\theta _{cr}-\pi )}{\textrm{d}\theta _{cr}}\end{aligned}$$
(10)

Muscle dynamic model

Fig. 2
figure 2

Schematic of the Hill-Huxley model.

Full size image

The electrical stimulation causes the muscle contraction process to activate the corresponding muscle fibers, and the degree to which the muscle fibers are activated is the degree of muscle activation. Schematic of the Hill-Huxley model is shown in Figure 2. The formula used to calculate muscle activation is:

$$\begin{aligned} \text {Act} = {\left\{ \begin{array}{ll} 0, & \text {pw}< \text {Thres} \\ \text {Sf} \frac{\text {pw} - \text {Thres}}{\text {Sat} - \text {Thres}}, & \text {Thres} \le \text {pw} < \text {Sat} \\ \text {Sf}, & \text {pw} \ge \text {Sat} \end{array}\right. } \end{aligned}$$
(11)

Where \(\text {pw}\) is the pulse width, \(\text {Thres}\) and \(\text {Sat}\) are threshold pulse width and saturation pulse width, respectively. The \(\text {Act}\) is muscle activation in the range of 0–1, \(\text {sf}\)is scale factor, takes the value of 1. In the second section of the muscle contraction module, muscle force is a function of muscle activation, muscle length, and velocity28. To simplify the calculation, the equations for angle and angular velocity are used as relations29.

$$\begin{aligned} & \textrm{M}_{{ma}}=\exp \left\{ -\left( \frac{\phi -\textrm{K}_1}{\textrm{K}2}\right) ^2\right\} \end{aligned}$$
(12)
$$\begin{aligned} & \textrm{f}_\textrm{v}=1+\frac{1}{\textrm{v}_\textrm{m}}\frac{\textrm{dl}}{\textrm{dt}} \end{aligned}$$
(13)
$$\begin{aligned} & \textrm{M}_{\textrm{mv}}=1+\frac{1}{\textrm{v}_\textrm{m}}\frac{\textrm{dl}}{\textrm{d}\phi }\frac{\textrm{d}\phi }{\textrm{dt}}=1-\frac{1}{\textrm{v}_\textrm{m}}\textrm{ma}\dot{\phi }(t)=1-\textrm{K}_3\dot{\phi } \end{aligned}$$
(14)

Where \(\phi\) is the knee angle \((0\le \phi \le 180)\), \(\textrm{M}_{{ma}}\) is the moment angle coefficient of the knee joint. \(\textrm{M}_{{mv}}\) is the angular velocity coefficient of the moment in the knee joint. l is the length of the quadriceps. \(\textrm{V}_{\textrm{m}}\) is the rate of contraction of the quadriceps. \(ma\phi\) is the power of the quadriceps muscle in the knee joint Rule; \(\textrm{k}_{\textrm{1}}\)and \(\textrm{k}_{\textrm{2}}\) are two parameters to be identified. Since \(\textrm{k}_{\textrm{3}}\) is difficult to determine, it is calculated from the available data described in the literature as \(\textrm{k}_{\textrm{3}}\)= 0.0430,31. The active torque of muscle activation is \(\mathrm {M_{act}=Act\cdot M_{ma}\cdot M_{mv}}\). In this model, active and passive muscle attributes are separate. The active moments belong to the active muscle properties. Passive muscle properties were calculated in the third part of muscle contraction dynamics. Passive muscle properties include passive damping moments \(\textrm{M}_{\textrm{vis}}={\mu \phi }^{\cdot }\), and passive stiffness moments \(\textrm{M}_{\textrm{cla}}=\lambda (\phi -\phi _0)\). Where \(\mu ,\lambda\) and \(\phi _{0}\) represent the hysteresis coefficient, elasticity coefficient and reference angle, respectively. The three parameters \(\mu ,\lambda\) and \(\phi _{0}\)are to be identified. According to the study of Ferrarin32, The kinetic equation provides the total external moment acting on the knee joint: \(\textrm{J}\ddot{\phi }=\textrm{M}_\textrm{act}+\textrm{M}_\textrm{grav}+\textrm{M}_\textrm{vis}+\textrm{M}_{cla}\).

Stimulation pattern

By referring to the existing FES cycling stimulation pattern, it was found that this study used a quantitative method to maximize the torque conversion rate from the knee joint to the crank by using Thres and Sat as the critical constraint values of the controlled object not only coordinates the movement of the leg muscles but also reduces muscle fatigue and maximizes muscle output. The stimulation pattern is related to the equipment model, the subject’s training posture and the subject’s height and weight information. In this study, the stimulation mode was optimized for the human body work efficiency as the optimization goal.

Firstly, input some basic parameters necessary for the design (human body size, experimental device size, etc.). Secondly, import these parameters into the already built lower limb dynamic model for cycling to obtain the motion state of each part of the lower limb. Thirdly, the effective moment and knee joint moment calculation was performed by combining the motion state with the muscle dynamics model. After the above calculation, it is easy to get the muscle output torque and average muscle energy expenditure during the movement. Finally, enter the calculation of energy efficiency and controller adjusted the intensity of the stimulation to the optimal stimulation mode so that the training in this mode can achieve the maximum average efficiency.

Calculation formula of instantaneous working efficiency in FES cycling: \(\eta =\frac{1}{\textrm{Angle}}\int _{\mathrm {A_1}}^{\mathrm {A_2}}\frac{\mathrm {Q_{cr}}}{\mathrm {M_{act}}}\). Where \(\textrm{Q}_{\textrm{cr}}\) is torque at the crank, \(\textrm{M}_{\textrm{act}}\) is torque at the knee joint, Angle is the angular interval, \(\textrm{A}_{\textrm{1}}\) and \(\textrm{A}_{\textrm{2}}\)is the initial angle and the termination angle of the stimulus, respectively. The stimulation angle interval was selected after the references as 10016,33,34.

Controller design

In order to achieve stable cycling in FES cycling, the PID controller was used to adjust the stimulus signal PW. The particle swarm optimization approach was proposed to optimize the BP neural network’s connection weight matrix in order to accelerate the learning convergence speed and keep it from entering local minima.

PID controller

Figure 3 depicts the PID control system’s fundamental architecture. The PID controller consists of three links: proportional, integral and differential, which are described mathematically as follows:

Fig. 3
figure 3

Basic structure of PID control system.

Full size image

Where \(K_{p}\) is a scale factor, \(T_{i}\) is integration time constant, and \(T_{d}\) is differential time constant. The proportional link instantly reflects the deviation signal e(t) proportionally, and the controller immediately generates a control effect to reduce it through negative feedback. The integral link is used to eliminate the static difference to improve the non-differential degree of the system. In contrast, the system’s performance indices acquired by the conventional tuning procedure are frequently subpar when dealing with complicated real-world controlled objects, such significant overshoot and lengthy correction times. The popular BP neural network offers high self-adaptability and self-learning capabilities for complicated and unpredictable applications. It has infinite accuracy for approximating any continuous nonlinear function.

PSO-BP-PID controller

The basic structure of the PSO-BP-PID control system is shown in Figure 4, where u(t) is the expected riding speed given by the system; v is the pulse width of the electrical stimulator in the system; Y is the actual speed at which the system expects the cycling training; e(t) is the error between the actual velocity and the expected velocity.

Fig. 4
figure 4

Basic structure of PSO-BP-PID control system.

Full size image

The basic idea behind PSO is to start with a collection of particles that have no mass or volume. Then, each particle is viewed as a potential solution to the particular kind of optimization issue, and the set fitness function is used to assess the particle quality35. The particle swarm optimization algorithm can be mathematically described as follows36:

$$\begin{aligned} \mathrm {v_{id}^{t+1}=wv_{id}^t+c1^*rand()^*(p_{id}-x_{id}^t)+c2^*rand()^*(p_{gd}-x_{id}^t)}\end{aligned}$$
(15)
$$\begin{aligned} \mathrm {x_{id}^{t+1}=x_{id}^t+v_{id}^{t+1}}\end{aligned}$$
(16)

Where \(r_{1}\) and \(r_{2}\) are random numbers between 0 and 1; \(c_{1}\) and \(c_{2}\) are acceleration constants, which represent the degree to which particles are influenced by social knowledge and individual cognition and are usually set to the same value. The w is the inertia weight factor, which is used to balance the global exploration and local development ability of particles. \(v_{id}^t\) and \(x_{id}^t\) are the velocity and position of the d dimension of particle i in the t iteration, respectively. \(p_{id}\) is the individual extreme position of particle i in the d dimension and \(p_{gd}\)is the global extreme position of the population in the d dimension. The BP neural network is widely used as a feedback network37,38. The basic idea of the BP-PID algorithm is to make the three neurons of output layer correspond to the three parameters of proportion, integral and differential in the PID structure. With the help of the self-learning capacity for the neural network, the output of the PID coefficient is controlled by adjusting the weight coefficient between layers. The structure of the BP neural network is shown in Figure 5.

Fig. 5
figure 5

The structure of the BP neural network.

Full size image

Experimental

Apparatus and system

To prevent subjects from having excessive lower limb movements, we designed an ankle-foot orthosis and fixed subjects’ lower limbs to the pedals. Using a modified pedal bike as an experimental platform, each lower leg was inserted into a fixation device fixed to the pedal so that leg movement was restricted to the sagittal plane. An inertial sensor was used to obtain the crank angle, and the stimulation pattern was implemented according to the crank angle. A personal computer was used to collect the crank angle data through a frequency of 200 Hz, and the individual controlled the electrical stimulator through a USB interface via the computer USB interface. Pulse width modulation (threshold-saturation) was used with balanced bipolar stimulation pulses with a constant frequency (20 Hz) and amplitude. The pedaling rhythm was maintained at a constant rate by the intensity of the stimulation delivered to the muscle. Schematic diagram of the experimental scene is shown in Figure 6.

Fig. 6
figure 6

Schematic diagram of the experimental scene.

Full size image

Subjects

We recruited six subjects from July to September 2022 at Shandong Provincial Third Hospital. The study has been approved by the Medical Ethics Committee of Shandong Provincial Third Hospital (KYLL-2020023). All subjects signed the informed consent form. All methods were performed in accordance with the relevant guidelines and regulations. The characteristics of the subjects are shown in Table 1. All subjects had no severe cognitive or communication problems. From Table 1, we can see that the study covered patients with different lower limb movement disorders. Subject 1 was a healthy person; subject 2, subject 4, and subject 6 were spinal cord injury patients; subject 3, and subject 5 were stroke patients. We selected stroke patients who had The NIH stroke scale (NIHSS) scores of 1–6. The NIHSS scoring system are: 0 = no stroke;1–4 = minor stroke; 5–15 = moderate stroke; 15–20 = moderate/severe stroke; 21–42= severe stroke17. The Hospital Ethics Committee approved the study, and all subjects signed informed consent. All subjects had no severe cognitive or communication problems.

Table 1 Characteristics of subjects (N=7).
Full size table

We collected all subjects’ characteristics, including information on height, weight, electrical stimulation pulse width threshold and saturation value. Each subject underwent three trials: (1) experiment without an electrical stimulator (traditional training), (2) experiment with a fixed electrical stimulation pulse width (training with fixed pulse width), and (3) experiment with the designed controller controlling the stimulation pulse width (training with self-adaptation control). Each experiment lasted 30 min with an interval of one day between tests to avoid muscle fatigue. And the experiments were executed sequentially, at a fixed point in time each day.

Performance indexes

An indicator of tracking accuracy is the root mean square error(RMSE).

$$\begin{aligned} \mathrm {RMSE=\sqrt{\frac{1}{T}\sum _{t=1}^T(x(t)-x_d(t))^2}}\end{aligned}$$
(17)

Where x denotes cycling speed, \(x_{dis}\) the desired speed, and T is the duration of cycling training. The absolute value of error (AVE)fluctuations was calculated as a measure of velocity fluctuations.The formula for AV is: \(\mathrm {AVE=max(x_n)-min(x_n)}\) Where \(\max (\cdot )\) is to take the maximum value, \(\min (\cdot )\) is to take the minimum value. The training duration and number of stops were calculated as a facilitation of FES.

Results

Figure 7 shows the results of a thirty-minute cycling experiment with two subjects are shown: Traditionally trained cycling speed, training with fixed pulse width cycling speed and training with self-adaptation control cycling speed. In the experiment, it was found that the change in pulse width in a short period of time had little effect on the running speed, so the average running speed of three revolutions was used as the running speed of the FES rehabilitation vehicle. We can see from Figure 7 that the subjects using the adaptive control scheme had a more stable cycling speed.

Fig. 7
figure 7

Different training modes cycling speed for subject 4 and subject 6.

Full size image

From Table 2, it is evident that under the conditions of electrical stimulation, the mean squared deviations are rather minimal. Table 3 shows that under the condition of FES participation, the speed mean square deviation of cycling training is relatively small, and in the three conditions using adaptive control were reduced in order. For example, Subject 6 had a self-adaptation control reduction of 0.644 and 0.7402 units compared to the traditional training mode and the fixed electrical stimulation mode, respectively. However, for subject 1, the velocity fluctuations in the conventional condition were smaller. Subjects’ actual cycling time (minus the intervening pauses for rest) during the thirty-minute subject time can be obtained from Table 4. It can be seen that the effective cycling time in the fixed pulse width electrical stimulation condition was greater than the conventional cycling time. The effective cycling time in the adaptive control electrical stimulation condition was greater than the cycling time in the fixed pulse width electrical stimulation condition. For example, the Cycling time of subject 2 in training with self-adaptation control was 578 s and 152 s longer than that of the traditional training mode and the fixed electrical stimulation mode, respectively. However, for subject 5, the cycling experiment in the adaptive control condition ended early without waiting for the end of the training condition, so there is a short ride time in the adaptive controller mode.

The experimental results show that the FES control system is more stable than both the traditional cycling mode and the fixed pulse width mode. Specifically,

$$\begin{aligned} \text {RMSE}_{\text {Fixed Pulse Width}}&< \text {RMSE}_{\text {Self-adaptation Control}} \\ \text {AVE}_{\text {Self-adaptation Control}}&< \text {AVE}_{\text {Fixed Pulse Width}}< \text {AVE}_{\text {Tranditional Traing Mode}} \\ \text {Training Time}_{\text {Traditional Training Mode}}&< \text {Training Time}_{\text {Fixed Pulse Width}}< \text {Training Time}_{\text {Self-adaptation Control}} \\ \text {Number of Stops}_{\text {Self-adaptation Control}}&< \text {Number of Stops}_{\text {Fixed Pulse Width}} < \text {Number of Stops}_{\text {Traditional Training Mode}} \end{aligned}$$

Therefore, the FES rehabilitation vehicle control system proposed in this study has good control performance for rehabilitation cycling training.

Table 2 RMSE of subjects speed in different training modes.
Full size table
Table 3 The Ave of subjects speed in different training modes.
Full size table
Table 4 Cycling time of subjects and the number of cycling pauses in different training modes.
Full size table

Discussion

The experimental results clearly show that by studying patients with different levels of severity, we have found that the worse the lower limb movement is for a given patient type, the better the control of the FES adaptive controller is. The traditional FES cycling stimulation pattern is still based on the experience of physiotherapist and is a fixed pattern16,39. This not only makes the rehabilitation of patients less efficient but also requires a large number of experienced physiotherapists. Our research is more stable than the optimization of other control algorithms, The RMSE achieved in our study was 4 units lower than that reported by NBoudville et al. (2018)40, 1.1 units lower than the results of Rouhani et al. (2017)41, and 0.16 units lower than the findings of de Sousa et al. (2016)42.

One of the main contributions of our research work is the proposed automatic stimulation timing based on the patient’s characteristics. The stimulation pattern can be obtained after the first measurement of the patient’s body parameters and can be used later in the training process. The stimulation pattern can be obtained dynamically from each patient to suit his or her needs. Another major contribution of this study is the proposed method of adaptive control of the FES. In this study, We proved this method can easily deploy on other systems and not require manual adjustment of the stimulation intensity and allows adaptive adjustment of the assisted stimulation pulse width according to the rehabilitation training speed of each patient. But we also faced the common limitation about may not reflect the true effect in larger or more heterogeneous populations due to the small sample size, resulting in limited external validity. The proposed modeling method requires more parameters to be identified, and therefore more preparatory experiments are needed.

There are some limitations in this study. We only recruited few participants. Future research should strengthen collaboration with healthcare organizations to access larger and more diverse datasets, enhancing the external validity of the findings and addressing potential underlying. Additionally, in terms of application, the developed model should be further optimized to achieve lightweight implementation without compromising accuracy, making it more feasible for deployment on endpoint devices.

Conclusion

In this study, pedal hill modeling was used to find the optimal stimulus mode with the maximum torque efficiency optimization objective, and the PSO-BP-PID control algorithm was used to construct the adaptive controller system. We found that the designed adaptive controller can effectively make the subjects cycling steadily at the expected speed. These results provide the possibility that subjects have longer continuous training times and more stable circuit training speeds during training.