Optimal Control Method for a Lower Limb Exoskeleton with Elastic Elements

Introduction. Lower limb exoskeletons (LLE) are of increasing interest due to the need to address global health issues, such as aging population and increase in neuromuscular injuries [1]. These devices are designed to provide effective solutions to support and improve human motor functions, such as walking assistance [2], rehabilitation and compensation for loss of balance [3], thereby increasing independence and quality of life.

LLE are often designed without elastic elements due to the increased complexity of the stabilization during movement and the influence of additional factors that arise when using elastic elements [4]. On the other hand, these exoskeletons may require elastic elements that improve the ability of the structure to adapt to uneven surfaces. Elastic elements can be installed in the ankle area [5] or used to completely replace it [6].

LLE control is a complex task. Various approaches to its solution have been proposed previously. The development of an effective control method depends on numerous factors, maintaining the constant relevance of research in this area. There are different methods for LLE control: adaptive control method [1][7], robust (stable) methods [8][9], and optimal control method [10][11]. Despite the efficiency of the first two methods, the third is the most successful. The optimal control method takes into account not only the increase in the stability and efficiency of LLE control under dynamic and unpredictable conditions, but also allows for reducing energy costs and the consumption of resources of the control system [10][11]. However, the integration of elastic elements into the LLE feet causes complications in their control. In addition, taking into account external disturbances when controlling the LLE is accompanied by new problems in providing the dynamic stability of the LLE control system.

Based on the above, it can be argued that there is a need to develop a method for optimal control of LLE with elastic elements. Therefore, the objective of this study was to develop a method for optimal control of the lower limb exoskeleton and elastic elements while optimizing energy costs and taking into account external disturbances. This method allows minimizing the quadratic function of energy costs in the presence of elastic elements and external disturbances.

One of the factors that further complicates the solution to the issue of dynamic stability of LLE is the presence of white noise, which is an interference to the control signal [12], which is investigated in this work.

To achieve the stated goal, the following tasks were set:

• development of a mathematical model of an inverted pendulum with elastic elements in the feet;
• development of a method for optimal control of LLE with elastic elements;
• accounting for the impact of external disturbances (Gaussian white noise);
• conducting numerical modeling in the Wolfram Mathematica environment for two types of controls: terminal control and feedback control;
• analysis of the obtained results of the transient processes study on the key parameters of the dynamics of LLE with elastic elements.

Materials. To develop a mathematical model of the LLE with elastic elements in the feet, its kinematic scheme is considered. The foot is represented by elastic elements and is connected to the model by the type of an inverted pendulum (Fig. 1). The model includes inertial properties, kinematic limitations of the joints, and external forces acting on the system.

The characteristics of the elastic elements and the parameters of the model are as follows: k_y, k_z — stiffness coefficients of the horizontal and vertical springs, respectively. Their values — k_y ∈ [ 500, 1500] N/m [13], k_z ∈ [ 7000, 20000] N/m [14]. c_y, c_z — horizontal and vertical damping coefficients, respectively, c_y ∈ [ 30, 200 N∙s/m], c_z ∈ [ 500–2000 N∙s/m] [15].

— moment of exoskeleton inertia relative to the center of mass, kgm²; m = 70 — mass of exoskeleton with patient, kg; h = 1 — length of exoskeleton to the center of mass, m; μ_y = 0.7 and μ_z = 0 — friction coefficients; N_y = 0 and N_z = m g — normal forces, N; g = 9.8 — acceleration of gravity, m/s².

Methods. The dynamics of the LLE is described by the Lagrange equations of the second kind in general form²:

(1)

where L = T – V — Lagrange function; T — kinetic energy of the system; V — potential energy of the system; Q — generalized forces; q = (θ, y_b, z_b)^T — vector of generalized coordinates; — rotary angle of exoskeleton link, measured from horizontal surface (parallel to reference plane) in counterclockwise direction; y_b, z_b — horizontal and vertical displacement of the base.

The coordinates of the center of mass of the pendulum are determined by the following equations:

(2)

The kinetic energy of a pendulum is given by the formula:

(3)

where v — speed of the center of mass, it is equal to .

The potential energy of the systems is described by the following equation:

(4)

The generalized damping forces applied through the Rayleigh dissipation function to simulate linear damping are written as [15]:

(5)

The Coulomb friction model is given by the following equations [15]:

(6)

where α >> 1 — regularization parameter (for approximating the discontinuous function sign(v)), in this paper, α = 100 is selected.

After intermediate calculations, equations (2–6) take the form:

(7)

where — external generalized forces.

In the vicinity of vertical equilibrium (θ ≈ 0), assumptions are valid: sin θ ≈ θ, cos θ = 1, z_b ≈ –mg/k_z. Equations (7) are linearized taking into account smallness and at low speeds:

(8)

Under the assumptions of small angles, the vertical displacement of z becomes dynamically independent of horizontal y and angular θ states. Its behavior is reduced to a harmonic oscillator, which can be analyzed separately [12].

For control synthesis, we present the system in the canonical form of the state space, taking into account external disturbances and measurement noise. Given the state vector , we write the system in the Cauchy form:

(9)

where A — state matrix; B — control matrix; u(t) — vector of control signals; D — matrix of disturbances; w(t) — vector of external disturbances; C — matrix of measurements; v(t) — vector of external disturbances.

The main task of the regulator is to transfer the dynamic system from the initial state x(t₀) to the specified final state x(t₁) in a certain time t₁.

Controllability Gramian W characterizes the ability of the system to reach arbitrary states in a finite time T and is determined by the formula [12]:

(10)

where Ф(t₀, t) = e^{A(t0 – t)} — matrix exponent.

The system is fully controllable on [t₀, t₁] if and only if W(t₁, t₀) is invertible. If the Gramian is invertible, any state x_T ∈ ℝⁿ can be reached using the appropriate control u(t) [12].

The control that minimizes the quadratic function of energy consumption and transfers the system from state x(t₀) to state x(t₁), has the form [12]:

(11)

After intermediate calculations, the quadratic function of energy consumption can be calculated using expression (12):

(12)

The control given by equation (11) is programmatic (open-loop) and time-dependent. To provide asymptotic stabilization of the system at the origin of coordinates under arbitrary initial conditions, it is required to synthesize the closed-loop control law based on real-time feedback. The control task is formulated in the form of minimizing the quadratic function of energy consumption [12]:

(13)

where R — positive definite matrix; Q — positive semi-definite matrix.

Optimal control that minimizes energy costs can be implemented as a negative feedback with variable parameters:

(14)

where k(t) = R^–1B^TP(t); k(t) — feedback coefficient matrix; P(t) — solution to the Riccati equation.

To find P(t) as t → ∞, it is required to solve the algebraic Riccati equation, which is given as follows:

(15)

The system with control signal u(t) is described by the equation:

(16)

In practice, measuring all states is impossible due to the limited number of sensors, noise and measurement errors [12]. To estimate the state vector from the inputs and outputs, an optimal Kalman filter is used.

The system with control signal u(t) and Kalman filter for determining the state vector has the following form [12]:

(17)

where — estimated state; — observer output; L — Kalman gain matrix.

In the presence of external disturbances w(t) and measurement noise v(t) with zero mathematical expectations M[v(t)] = 0, M[w(t)] = 0 and covariance matrices s_w and s_v, respectively, gain matrix L is determined from the solution to the algebraic Riccati equation (18), and, accordingly, stationary Kalman filter is found:

(18)

Control of a linear non-stationary system with external disturbances at the input and output is implemented using linear feedback on the state estimate:

(19)

The type of closed-loop system with feedback for estimating the state vector is determined by the formula:

(20)

After developing the optimal control method for LLE with elastic elements and taking into account the impact of external disturbances, it remains to perform numerical modeling and analysis of the obtained results of the study on transient processes for the angle of deviation and displacement. For this purpose, at the first stage, the optimal control method is focused on restoring equilibrium for a certain stabilization time without taking into account the intermediate trajectory of motion. At the second stage, the optimal feedback control method is used to find the optimal intermediate trajectory of motion both with and without the Kalman filter.

Research Results. The paper investigates the transient processes of the control moment, force, displacement and angle during the transition of the LLE with elastic elements from an unstable position to a vertical equilibrium position. The results of numerical modeling are shown in Figures 2–6 for three initial conditions: Y₁ = [ 1/10, 0, 0 ,0], Y₂ = [ 1/10, 0, 1/10 ,0], Y₃ = [ 1/10 0 1/10 0 0 1/10 0 1/20 0].

The results of terminal control are calculated for a stabilization time of 0.6 s. Figure 2 shows the curves of the control moments of the open-loop system under the initial conditions Y₁ and the target position in the vertical equilibrium state. The blue curve corresponds to control moment M, the orange one — to control force F_y.

The values of the quadratic function of energy consumption for terminal control for different values of stabilization time are calculated using equation (12) and are presented in Table 1. A decrease in the value of the function of energy consumption with an increase in stabilization time can be noticed.

Figure 3 shows the curves of transient processes of the key parameters of the dynamic LLE system (θ, y_b) under the initial conditions Y₁. The blue curve represents the transient process of the deviation angle θ, and the orange curve — change in the coordinate along the ordinate axis y_b during the specified stabilization time of 0.6 s under terminal control.

The simulation results in Figures 4–6 are obtained using the feedback control type, which allows searching for the optimal intermediate trajectory of motion. To study the impact of the elastic elements’ stiffness (changes in spring stiffness coefficient k_y) on the stability of the dynamic LLE system while minimizing energy costs, transient processes of the key parameters of the dynamic LLE system (Fig. 4) were considered under the initial conditions Y₂ and the following parameters of the feedback system: c_y = 100, k_z = 20,000, Q = eye(4) ⋅ 10³, R = eye(2).

The orange line represents dependence, and the green line — . The following values of the stiffness coefficient are selected: k_y = [ 700, 1000, 1500] N/m. The transient processes are shown in Figure 4 a, b, c, respectively.

To study the effect of damping of elastic elements (changes in damping coefficient c_y) on the stability of the dynamic LLE system while minimizing energy costs, transient processes of the key parameters of the dynamic LLE system(θ, y_b) (Fig. 5) are considered under the initial conditions Y₂ and k_y =1000 with feedback. The orange line represents the y_b dependence, and the green line — θ. The damping coefficient values are c_y = [ 0, 50, 100] N∙s/m, and the transient processes are shown in Figure 5 a, b and c, respectively.

This research also takes into account the effect of noise through adding Kalman filter. Figure 6 shows the transient response curves under the initial conditions Y₃ = [ 1/10 0 1/10 0 0 1/10 0 1/20 0] and the covariance matrices s_w = eye(1) and s_v = 0.1 eye(2). The orange line illustrates the transient response y_b, and the green line — transient response θ.

Discussion. To achieve the research objective, the previous section presented the simulation results obtained using the proposed optimal control method to minimize energy costs during the transition of the studied LLE from an unstable position to a vertical equilibrium position.

Table 1 shows that as the time required to reach the target stable state increases, energy costs decrease proportionally. This dependence is consistent with the dynamics of control signals: as the permissible stabilization time increases, the amplitudes of the control moments M and decrease, which, according to formulas (12)–(14), reduces the energy costs. It can be concluded that an increase in the stabilization time (within 0.1–0.4 seconds) reduces energy costs by 98%.

The curves in Figure 3 demonstrate a transient process lasting 0.6 s before reaching the vertical equilibrium position, which is consistent with the principles of terminal control (reaching the target stable state within a certain stabilization time).

The curves in Figure 4 illustrate the effect of changing the stiffness coefficient on the transient processes of the key parameters of the LLE dynamic system (θ, y_b). Changing the coordinate along the ordinate axis y_b at a low value of the stiffness coefficient is characterized by significant overshoot and slow stabilization, and an increase in the stiffness coefficient eliminates overshoot and accelerates the stabilization of the system. On the contrary, changes in angle θ demonstrate an inverse relationship: at a low value of the stiffness coefficient, a smooth transient process without overshoot and a short stabilization time of the system are observed, while an increase in the stiffness coefficient to 1500 causes overshoot in the angle, despite accelerated stabilization. This contradiction emphasizes the competing dynamics between y_b and θ. At the stiffness coefficient k_y =1000 (Fig. 4 b), an optimal compromise is reached: minimization of overshoot in the coordinate y_b, maintaining the stability of angle θ and fast convergence. This mode provides balanced operation of the system.

The reliability of the results obtained is validated by the stability of the closed system, which is confirmed by the negative real parts of its poles for all the studied values of the stiffness coefficient given in Table 2.

The curves in Figure 5 show the effect of the damping coefficient. At a low value of the damping coefficient, significant vibration of the system is observed before reaching the steady state, accompanied by overshoot. Increasing c_y reduces the amplitude of vibrations and eliminates overshoot, but results in a growth of energy costs. To provide a balance between stability and the value of the quadratic function of energy costs, the value of coefficient c_y = 100 is selected, which guarantees the stability of the LLE control system with elastic elements.

From the image in Figure 6, it follows that the Kalman filter provides the convergence of transient processes to zero values within three seconds when exposed to white noise, confirming its robust stabilizing properties.

The results obtained show the following picture. Firstly, the use of open control as a basic stage allows minimizing energy costs at the stage of bringing the system from an unstable state to the target vertical equilibrium at a given stabilization time. Secondly, the transition to closed loop control based on the state feedback provides asymptotic stability and stability of transient modes through solving the Riccati equation and implementing an effective Kalman filter. Thirdly, numerical modeling revealed an important effect of the parameters of elastic elements on the dynamics of the system: increasing the spring stiffness reduces overshoot in angle θ and accelerates stabilization, but it can cause difficulties on the way to a stable trajectory; increasing damping reduces vibrations and reduces overshoot, but increases energy costs. A balance between stability and energy costs is found, which is reached at specific values k_y and c_y (in the examples: k_y = 1000 N/m, c_y = 100 N·s/m).

Based on the previous discussion and analysis of the results obtained in the paper, it can be said that the developed method of optimal control of the LLE with elastic elements has managed to provide the stability of the control system of this exoskeleton for a certain stabilization time and with minimal energy consumption.

Conclusion. A technique for optimal control of a lower limb exoskeleton (LLE) with elastic elements in the feet, taking into account external disturbances and measurement noise, has been formulated and implemented. The main approach is based on representing the LLE dynamics as a system of Lagrange equations, translating it into a canonical form of the state space, and synthesizing the control law through optimization of quadratic functions of energy consumption and system stability. A Kalman filter was used to assess the state, which allowed for correct operation under conditions of a limited number of sensors and the presence of external disturbances.

The practical significance of the results consists in the development of a technique that allows adaptively selecting the parameters of elastic elements and the control mode depending on the conditions of the task and objectives (minimization of energy costs, acceleration of stabilization, minimization of overshoot). In terms of potential applications, this can help improve the efficiency of rehabilitation technologies, reduce energy consumption in prosthetic and orthotic systems, and improve the stability of movements on uneven surfaces. The prospects for the research include experimental verification of the technique and its adaptation to variable loads and complex surfaces, which is a challenge for medical rehabilitation and robotics.