Review of human–robot coordination control for rehabilitation based on motor function evaluation

Abstract

As a wearable and intelligent system, a lower limb exoskeleton rehabilitation robot can provide auxiliary rehabilitation training for patients with lower limb walking impairment/loss and address the existing problem of insufficient medical resources. One of the main elements of such a human–robot coupling system is a control system to ensure human–robot coordination. This review aims to summarise the development of human–robot coordination control and the associated research achievements and provide insight into the research challenges in promoting innovative design in such control systems. The patients’ functional disorders and clinical rehabilitation needs regarding lower limbs are analysed in detail, forming the basis for the human–robot coordination of lower limb rehabilitation robots. Then, human–robot coordination is discussed in terms of three aspects: modelling, perception and control. Based on the reviewed research, the demand for robotic rehabilitation, modelling for human–robot coupling systems with new structures and assessment methods with different etiologies based on multi-mode sensors are discussed in detail, suggesting development directions of human–robot coordination and providing a reference for relevant research.

Graphical abstract

Keywords

human–robot coupling / lower limb rehabilitation / exoskeleton robot / motor assessment / dynamical model / perception

Introduction

Topology optimization (TO), which has been extensively studied over the last decades, is a process of determining optimal layout of materials inside a given design domain. TO has been applied to various structural optimization problems, such as minimum compliance [1,2] vibration [3,4], and thermal issues [5,6], after Bendsøe and Kikuchi [7] proposed the homogenization method. Homogenization is a material distribution method in which a design domain is discretized into small rectangular elements, and each element contains an artificial composite material with microscopic voids. The proposal of the homogenization method was followed by a parallel exploration of the solid isotropic material with penalization (SIMP) method [8,9], which uses an artificial isotropic material whose physical properties are expressed as a function of continuous penalized material density (design variables). The phase-field method [10,11], which is also a material distribution method, is based on the theory of phase transitions. A different type of method called evolutionary structural optimization (ESO) [12,13] has also been proposed. This method eliminates elements with the lowest criterion value on the basis of certain heuristic criteria. ESO is computationally expensive because it requires a much larger number of iterations with an enormous number of intuitively generated solutions compared with material-based methods.

However, these conventional TO methods, which are based on element-wise design variables, suffer from numerical instability problems, such as checkerboards and mesh dependency. Accordingly, several studies have proposed prevention methods. The use of high-order elements has been proven to be an effective means to prevent checkerboards [14,15], but this method entails a considerable increase in computation time. Various filter techniques have been utilized to mitigate checkerboards and mesh dependency because these techniques require only a small amount of extra computation time and are simpler to implement than other methods [16,17]. Notably, filter schemes are purely heuristic. Other prevention schemes, such as perimeter control and gradient constraint, which often make the optimization procedure unstable, are yet to be improved.

A new type of TO approach is the level set-based method (LSM), which was developed by Sethian and Wiegmann [18] to numerically track the motion of structural boundaries. In LSM, boundaries are represented as zero level set contour of an implicit high-dimensional function (level set function or LSF), in which boundary motion, merging, and introduction of new holes are performed. The evolution of a structural interface with respect to time is tracked by solving a Hamilton-Jacobi (H-J) partial differential equation (PDE), where transporting LSF along its outward normal direction is equivalent to moving the boundaries along the descent gradient direction. The conventional level set-based TO approach uses shape derivatives coupled with the original LSM for boundary tracking [19,20]. In this approach, regularization that reinitializes LSF to be signed distance to a zero level set is employed to control the slope of LSF. This conventional approach is updated by solving the H-J equation via an explicit up-wind scheme [21,22]. Variations of the conventional approach include parameterizing LSF using various basis functions, such as finite element method (FEM) basis functions [23], radial basis functions (RBFs) [24,25], and spectral parameterization [26], and corresponding methods for solving the H-J equation. By defining the interfaces between two material phases via the iso-contour of LSF, LSM can handle shape and topology changes during the optimization procedure and provides optimal structures with clear boundaries that are free of checkerboard patterns. Notably, most LSMs rely on finite elements wherein boundaries are still represented by discretized mesh in the analysis field unless alternative techniques are utilized to map the geometry to the analysis model.

Most of these TOs are performed in a fixed domain of finite elements where FEM is used to solve optimization problems. Currently used FEMs are often based on Lagrange polynomials for analysis while the geometrical representation of structures relies on non-uniform rational B-splines (NURBS), which are the criteria in computer aided design (CAD) systems. Thus, conversion of NURBS-based representation into one that is compatible with Lagrange polynomials, that is, mesh generation, is required in structural analysis. The disadvantages of FEM are as follows. First, the geometry approximation inherent in the FEM mesh may generate an approximate error. Second, frequent data interaction between geometry description and the computational mesh, which can be found in several calculations (e.g., fluid, large deformation, and shape optimization problems), is cumbersome and error-prone. An integrating method, namely, isogeometric analysis (IGA) [27], for unifying analysis and CAD processes has been proposed to overcome these disadvantages. This method employs the same basis functions as a technique for describing and analyzing the geometric model, which features the IGA method and CAD-based parameterization of field variables in an isoparametric manner. The first work on isogeometric approximation dates back to 1982 [28]; however, this method is considerably different from the IGA method. Several methods have been devised to help alleviate the difficulties faced by IGA. Special parameterization techniques, such as variational harmonic-based methods [29,30] and analysis-aware parameterization methods for single [31] and multi-domain geometries [32], have been proposed for the computational domain. Alternatives to NURBS, such as T-splines [33,34], polynomial splines over hierarchical T-meshes (PHT-splines) [35–37], and Powell-Sabin splines [38], have been studied for local refinement in IGA due to the limitation of the tensor product form of NURBS in computation refinement. Methods of parameterization of the interior domain while retaining the geometry exactness from the CAD model have been devised [39,40], and isogeometric collocation method is one of the most important among these methods [39]. With regard to interior discretization obstacles, the isogeometric boundary element method is a suitable candidate [41] because only boundary data are required for analysis, and it enables stress analysis [42], fracture analysis [43,44], acoustic analysis [45], and shape optimization [46,47]. Considering that the integral efficiency of IGA is limited by the tensor product structure of NURBS, an efficient quadrature rule, which is more suitable for NURBS-based IGA compared with the Gaussian quadrature rule, has been proposed in Ref. [48]. IGA has been applied to a wide range of problems, such as structural vibrations [49], fluid-structure interaction [50,51], heat conduction analyses [52], shape optimization [53,54], shell analyses [55], TO [56], and electromagnetics [57].

IGA-based shape optimization has been extensively investigated because remeshing is eliminated during the optimization process. IGA has also been recently applied to TO. The most commonly used TO approach is material based, and the most commonly solved TO problem is the minimum compliance case. In Ref. [58], TO was proposed based on isogeometric shape optimization. B-spline curves were introduced to represent the material boundary, and the coordinates of their control points were considered design variables. In Ref. [59], the trimmed spline surface technique was used for spline-based TO. In Refs. [60,61], optimality criteria (OC) and the method of moving asymptotes were implemented in the isogeometric-based SIMP method. In Ref. [62], the TO problem was solved by using a phase-field model, and IGA was utilized for the exact representation of the design domain. In Ref. [56], IGA was introduced to a level set TO method, and NURBS basis functions were used for geometry description and LSF parameterization.

Most studies on TO were restricted to compliance optimization, and the number of studies on TO of dynamic problems is limited. Dynamic problems, such as vibration and noise, are critical in many engineering fields, such as aeronautical and automotive industries. As a typical dynamic problem, structure vibration is controlled by the structure’s dynamic characteristics, which are usually represented by the eigenfrequencies of the structure. Thus, eigenfrequency optimization plays an important role in improving dynamic characteristics.

As an important research topic, eigenfrequency optimization has been studied by many scholars. Díaaz and Kikuchi [63] extended TO to eigenfrequency optimization within the framework of the homogenization approach. Tenek and Hagiwara [64] introduced the SIMP method to structural shape optimization and TO in eigenfrequency problems. Xie and Steven [65] applied the ESO scheme to TO, and the specified eigenfrequency of the structure was used as a constraint. Additional research has also been conducted on frequency optimization [66–70]. However, standard density-based TO methods are unsuitable for eigenfrequency optimization due to localized modes in low-density areas [71]. Low-density areas are much more flexible than areas with full densities; hence, they control the lowest eigenmodes of the entire structure. By changing the penalization of stiffness in the SIMP method, a modified algorithm has been proposed and applied to circumvent this numerical instability problem [3]. LSM that employs a crisp description of structure boundaries has advantages over the density-based approach. The LSM approach can avoid artificial modes localized in the weak phase, which makes LSM a choice for eigenfrequency optimization. Many studies have been performed on level set TO for eigenfrequency problems [72–75].

IGA-based TO has been extensively studied. However, research on the combination of the level set approach and IGA is limited, and that on eigenfrequency problems is absent. In this study, we develop a new optimization method to formulate the TO problem for cases with maximum fundamental eigenfrequency by using LSM under the framework of IGA instead of conventional finite element analysis (FEA). Given that the OC algorithm is particularly efficient for problems with many design variables and few constraints, we consider the OC method for the solution of the optimization problem and conduct a sensitivity analysis. The rest of this paper is organized as follows. In Section 2, IGA and NURBS-based TO are summarized. The eigenvalue optimization problem is described in Section 3, and the TO model is proposed in Section 4. Section 5 presents numerical examples to demonstrate the validity of the proposed approach. The conclusions and discussions are shown in Section 6.

IGA for level set-based TO

In IGA, geometric modeling and analysis are integrated by using NURBS, where the basis function is a bridge of the parameter domain, physical field, and numerical solution. In the proposed method, we consider the NURBS basis function as a bound between IGA and parameterized LSM. We provide a brief review of IGA and NURBS-parameterized LSM.

NURBS-based IGA

We assume that geometrical mapping $\mathrm{\Psi }$ maps the parameter domain $\widehat{\mathrm{\Omega }}$ into the physical domain $\mathrm{\Omega }$. Given a knot vector $\mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s)$ with a non-decreasing sequence of values lying in parameter space, the mapping between two domains can be expressed as

(1)$\mathrm{\Psi }:\widehat{\mathrm{\Omega }}\to \mathrm{\Omega },\mathbf{\Xi }\to \mathrm{\Psi }(\mathbf{\Xi }).$

The NURBS curve is constructed by linear combinations of its basis functions, in which the coefficients are a given set of control points. A NURBS curve of p-degree is defined as

(2)$\mathrm{\Psi }(\xi )={\displaystyle \sum _{i=\mathrm{1}}^n{R}_{i,p}}(\xi )P_i,$

where $n=s-p-\mathrm{1}$ is the number of control points, $P_i\in {ℝ}^d$ is the ith control points in the physical domain, and $R_{i,p}(\xi )$ is the ith univariate NURBS basis function defined in the parameter space $\widehat{\mathrm{\Omega }}$ as follows:

(3)$R_{i,p}(\xi )={\displaystyle \frac{N_{i,p}(\xi ){\omega }_i}{{\displaystyle \sum _{i=\mathrm{1}}^n{N}_{i,p}}(\xi ){\omega }_i}},$

where ${\omega }_i\in (\mathrm{0},\mathrm{1})$ are non-decreasing weights associated with control points and $N_{i,p}(\xi )$ represents the ith B-spline basis functions of p-degree; it is defined by the following Cox-de Boor recursion formula [27].

(4)$\{\begin{array}{l}{N}_{i,\mathrm{0}}(\xi )=\{\begin{array}{lll}\mathrm{1}& & \mathrm{if} {\xi }_i\le \xi <{\xi }_{i+\mathrm{1}}\\ \mathrm{0}& & \mathrm{otherwise}\end{array}\\ N_{i,p}(\xi )={\displaystyle \frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}}{N}_{i,p-\mathrm{1}}(\xi )+{\displaystyle \frac{{\xi }_{i+p+\mathrm{1}}-\xi }{{\xi }_{i+p+\mathrm{1}}-{\xi }_{i+\mathrm{1}}}}{N}_{i+\mathrm{1},p-\mathrm{1}}(\xi )\end{array},$

where basis function $N_{i,p}(\xi )$ has its own support domain $[{\xi }_i,\dots ,{\xi }_{i+p+\mathrm{1}}]$ in which $N_{i,p}(\xi )$ is non-zero. A knot vector is deemed open when the knots are repeated $p+\mathrm{1}$ times at the ends of the vector. In IGA, the open-knot vector is generally used to satisfy the Kronecker delta property at boundary points.

A NURBS surface is a tensor product of bivariate NURBS curves in $\mathbf{\Xi }$ and $\mathrm{H}$ directions with p- and p-degrees, respectively, where a knot vector $\mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t)$ is given in $\mathrm{H}$ direction.

(5)$\mathrm{\Psi }(\xi ,\eta )={\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{R}_{i,j}^{p,q}}}(\xi ,\eta )P_{i,j},$

where $m=t-q-\mathrm{1}$, $P_{i,j}$ are control points and $R_{i,j}^{p,q}$ are bivariate basis functions of the form:

(6)$R_{i,j}^{p,q}={\displaystyle \frac{N_{i,p}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}{{\displaystyle \sum _{i=\mathrm{1}}^n{\displaystyle \sum _{j=\mathrm{1}}^m{N}_{i,p}}}(\xi )N_{j,q}(\eta ){\omega }_{i,j}}},$

where $N_{i,p}(\xi )$ and $N_{j,q}(\eta )$ are B-spline basis functions defined on the knot vector $\mathbf{\Xi }=({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_s)$ and$\mathrm{H}=({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\dots ,{\eta }_t)$, respectively. The interval $\mathbf{\Xi }\times \mathrm{H}$ forms a patch containing all elements, namely, $[{\xi }_k,{\xi }_{k+\mathrm{1}}]\times [{\eta }_l,{\eta }_{l+\mathrm{1}}]$, $\mathrm{1}\le k\le n+p$, and $\mathrm{1}\le l\le m+q$, which are defined by the two knot vectors. The parameter domain and corresponding physical domain for a surface model are depicted in Fig. 1.

Fig.1 Geometrical mapping $\mathrm{\Psi }$maps the common parameter space $(\xi ,\eta )$ onto the physical space

Full size|PPT slide

On the basis of the isoparametric concept, the IGA approach utilizes the same parameters for geometry and analysis models, and the basis functions used for geometry representation are also employed to approximate the numerical solution of PDEs. With Eq. (2), the numerical solution $u$ can be expressed as

(7)$u={\displaystyle \sum _{i=\mathrm{1}}^n{R}_i}(\xi )u_i,$

where $R_i$ is the ith basis function. $u_i$, which is referred to as a control variable at the ith control point, is the coefficient used to approximate the field variable $u$, which plays the same role as the nodal value in FEA. For each element, the shape function and strain-displacement matrix can be expressed as

(8)$\mathit{R}=[R_{\mathrm{1}}{R}_{\mathrm{2}}\mathrm{...}{R}_n],\mathit{B}=\left[\begin{array}{lllll}{\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{0}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial x}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& \mathrm{...}& \mathrm{0}& {\displaystyle \frac{\partial R_n}{\partial y}}\\ {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial y}}& {\displaystyle \frac{\partial R_{\mathrm{1}}}{\partial x}}& \mathrm{...}& {\displaystyle \frac{\partial R_n}{\partial y}}& {\displaystyle \frac{\partial R_n}{\partial x}}\end{array}\right].$

The strain matrix is given by

(9)$\mathit{\epsilon }=\mathit{B}\mathit{u}.$

NURBS-based parameterized LSM

LSM is a TO method that implicitly defines the interfaces between material and void phases by iso-contours of a high-dimensional LSF. Thus, this method allows a crisp description of the material boundary and helps avoid mesh-dependent problems.

The shape of the interpolating functions of LSF directly influences the smoothness of LSF and the material domain. In its most general form, LSF is described by a summation of interpolating functions scaled by their degree of freedom (DOF).

(10)$\mathrm{\Phi }(\mathit{x})={\varphi }^{\mathrm{T}}(\mathit{x})\mathit{s}={\displaystyle \sum _i^N{\varphi }_i(\mathit{x})s_i},$

where $\mathit{x}$ denotes the spatial coordinate, ${\varphi }_i$ comprises interpolating functions associated with $N$spatial points, and $s_i$ are time-dependent optimization variables.

The most commonly used interpolation functions in present LSM are FEM shape functions and RBFs, and their corresponding optimization variables are nodal values and expansion coefficients, respectively. We introduce NURBS basis functions, which can be used to approximate a given set of points with smooth polynomial functions, for parameterizing LSF. Thus, the NURBS-based parameterized LSF is constructed as

(11)$\mathrm{\Phi }(x)={\mathit{R}}^{\mathrm{T}}\mathit{s}={\displaystyle \sum _i^N{R}_i(\mathit{x})s_i},$

where $R_i$ is the ith NURBS basis function and $s_i$ is the ith time-dependent expansion coefficient related to the ith control point.

The evolution of LSF is governed by the following H-J equation.

(12)${\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}+\nabla \mathrm{\Phi }{\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}}=\mathrm{0},$

where $t$ denotes pseudo-time, which represents the evolution of the design in the optimization process. The speed of movement of a point on the level set surface can be expressed by $V=\mathrm{d}\mathit{x}/\mathrm{d}t$. ${\mathit{V}}_{\mathrm{n}}=V\cdot \mathit{n}$ defines the speed of propagation of all level sets along the external normal velocity, where $\mathit{n}=-\nabla \mathrm{\Phi }/|\nabla \mathrm{\Phi }|$. Therefore, Eq. (12) can be rewritten as

(13)${\displaystyle \frac{\partial \mathrm{\Phi }(\mathit{x},t)}{\partial t}}={\mathit{V}}_{\mathrm{n}}|\nabla \mathrm{\Phi }|.$

By substituting Eq. (11) into Eq. (13), the H-J equation can be written as

(14)${\mathit{R}}^{\mathrm{T}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}+{\mathit{V}}_{\mathrm{n}}\cdot \left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|=\mathrm{0.}$

The moving speed of the material free boundary during evolution is related to the time derivative of the expansion coefficient as follows:

(15)${\mathit{V}}_{\mathrm{n}}=-{\displaystyle \frac{{\mathit{R}}^{\mathrm{T}}}{\left|{(\nabla \mathit{R})}^{\mathrm{T}}\mathit{s}\right|}}{\displaystyle \frac{\partial \mathit{s}}{\partial t}}.$

Optimization problems of maximizing eigenvalue

Definition of the eigenvalue problem

We describe the eigenvalue problems in the linear elastic continuum to facilitate the computation of vibration frequencies and modes. A linear elastic continuum structure with a constant mass density is defined in domain $\mathrm{\Omega }\subset {ℝ}^d(d=\mathrm{2}\mathrm{o}\mathrm{r}\mathrm{3})$ with the boundary $\mathrm{\Gamma }=\partial \mathrm{\Omega }$. The weak formulation of the undamped free vibration problem can be expressed as

(16)$a(\mathit{u},\mathit{v})-\lambda b(\mathit{u},\mathit{v})=\mathrm{0},\mathit{v}\in U,$

where eigenfrequency $\lambda$ and corresponding eigenvector $\mathit{u}$, that is, the displacement subdomain in $\mathrm{\Omega }$, are solutions of this eigenvalue problem, $\mathit{v}$ is adjoint displacement, which satisfies the kinematic boundary condition, and $U$ is a space of kinematically admissible displacement fields. In LSM, $a(\cdot ,\cdot )$ and $b(\cdot ,\cdot )$ are respectively defined as

(17)$a(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },$

(18)$b(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega },$

where $D_{ijkl}$ stands for the elasticity tensor component, ${\epsilon }_{ij}$ is the strain tensor component, $\rho$ is the density of the material, and $H(\mathit{\Phi })$ is the Heaviside function, which takes $\mathrm{0}$ for $\mathit{\Phi }<\mathrm{0}$ and $\mathrm{1}$ otherwise.

The eigenvalue problem has a family of solutions ${\lambda }_k$ and $u_k$, $k\ge \mathrm{1}$. The first eigenfrequency and its eigenvector are related to each other as

(19)${\lambda }_{\mathrm{1}}=\min {\displaystyle \frac{{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}}{\epsilon }_{kl}(\mathit{u}){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}{{\displaystyle {\int }_{\mathrm{\Omega }}\rho }\mathit{u}\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}}.$

Optimization model

We consider the TO problem by maximizing the first eigenfrequency under a volume constraint. Under the NURBS-based level set framework, eigenfrequency TO can be expressed as

(20)$\begin{array}{l}\mathrm{Maximize}\mathrm{\hspace{1em}\hspace{1em}}J(\mathit{u},\mathrm{\Phi })={\lambda }_{\mathrm{1}}\\ \mathrm{subject}\hspace{0.17em}\mathrm{to}:\hspace{1em}\hspace{0.17em}a(\mathit{u},\mathit{v},\mathrm{\Phi })={\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi }),\\ {\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }\le V_{\max },\\ s_{\min }\le s_i\le s_{\max },\mathrm{\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}} \end{array}$

where $J(\mathit{u},\mathrm{\Phi })$ is the objective function, $V_{\max }$ represents the maximum admissible volume of the design domain, and $s_{\min }$ and $s_{\max }$ stand for the lower and upper bounds of the design variables, respectively.

However, in the eigenfrequency optimization process, the value of higher-order eigenfrequency may decrease whereas that of lower-order target eigenfrequency may increase, which may possibly lead to the repetition and exchange of mode order number. Given that the objective and constraint functions are typically defined based on a fixed modal order, the sensitivities of these functions are discontinuous in the repeated eigenfrequency case. Approaches are often used to maintain the simplicity of eigenfrequency during the entire optimization process and overcome this ill-posed problem. The modal assurance criterion (MAC) method, an efficient and accurate strategy, is introduced to monitor a single target mode, which is the first mode in this study. The definition of MAC is

(21)$\mathrm{MAC}({\mathit{u}}_a,{\mathit{u}}_b)={\displaystyle \frac{|{\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_b{|}^{\mathrm{2}}}{({\mathit{u}}_a^{\mathrm{T}}{\mathit{u}}_a)({\mathit{u}}_b^{\mathrm{T}}{\mathit{u}}_b)}},$

where ${\mathit{u}}_a$ and ${\mathit{u}}_b$ represent two eigenvectors: One is the reference eigenvector of the current optimization cycle and the other is the objective eigenvector of the previous cycle in the eigenvalue optimization process. The value of MAC varies between 0 and 1. Theoretically, a MAC value of 1 means that the two eigenvectors representing modal shapes are exactly the same. However, this condition is impossible because the structural configuration changes in each iteration, and the modal shapes of adjacent iterations are not orthogonal to each other. By comparing a few reference eigenvectors with an objective eigenvector ${\mathit{u}}_{\mathrm{obj}}^{}$, a new objective eigenvector is obtained in each iteration step of the optimization process, and it can be expressed as

(22)${\mathit{u}}_{\mathrm{obj}}^n={\mathit{u}}_k^n\hspace{1em}\mathrm{that} \underset{k}{\max }[\mathrm{MAC}({\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}},{\mathit{u}}_k^n)], k=\mathrm{1},\mathrm{2},\dots ,N_m,$

where $N_m$ is the number of modal eigenvectors that need to be checked, superscript of displacement indicates the number of iteration step, that is, in the nth iteration, objective eigenvector of previous iteration step ${\mathit{u}}_{\mathrm{obj}}^{n-\mathrm{1}}$ is used to calculate the MAC value.

Sensitivity analysis

Establishing the relationship between the optimization function and design variables by using a sensitivity analysis approach is necessary to solve the optimization problem. According to the material derivative and the adjoint method, the Lagrangian function can be defined as

(23)$L(\mathit{u},\mathrm{\Phi })=J(\mathit{u},\mathrm{\Phi })+a(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }\left[{\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max }\right],$

where $\mathrm{\Lambda }$ is the Lagrangian multiplier. Assuming that $V(\mathit{u},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}H}(\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }-V_{\max }$ is the volume constraint, the shape derivative of Lagrangian function $L(\mathit{u},\mathrm{\Phi })$ is

(24)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={{\lambda }_{\mathrm{1}}}^{\prime }+a^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}{b}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })-{{\lambda }_{\mathrm{1}}}^{\prime }b(\mathit{u},\mathit{v},\mathrm{\Phi })+\mathrm{\Lambda }{V}^{\prime }(\mathit{u},\mathrm{\Phi }) ,$

where the material derivatives of $a(\mathit{u},\mathit{v},\mathrm{\Phi })$ and $b(\mathit{u},\mathit{v},\mathrm{\Phi })$ are respectively given by

(25)$\begin{array}{l}{a}^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}({\mathit{v}}^{\prime })H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ +{\displaystyle {\int }_{\mathrm{\Omega }}{D}_{ijkl}{\epsilon }_{kl}({\mathit{u}}^{\prime }){\epsilon }_{ij}(\mathit{v})\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega ，}}\end{array}$

(26)$b^{\prime }(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_{\mathrm{\Omega }}\rho {\mathit{u}}^{\prime }\mathit{v}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}{\mathit{v}}^{\prime }H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\mathrm{\Omega }}\rho \mathit{u}\mathit{v}\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},$

where ${\mathit{u}}^{\prime }$ and ${\mathit{v}}^{\prime }$ are partial derivatives of $\mathit{u}$ and $\mathit{v}$, respectively, with respect to pseudo-time. $\delta (x)$ is the Dirac function.

The adjoint state equation can be obtained by the Kuhn-Tucker condition.

(27)$a({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b({\mathit{u}}^{\prime },\mathit{v},\mathrm{\Phi })=\mathrm{0},$

(28)$a(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })-{\lambda }_{\mathrm{1}}b(\mathit{u},{\mathit{v}}^{\prime },\mathrm{\Phi })=\mathrm{0},$

(29)$\mathrm{1}-b(\mathit{u},\mathit{v},\mathrm{\Phi })=\mathrm{0.}$

Given that the real mode $\mathit{u}$ is equal to the adjoint mode $\mathit{v}$, Eq. (29) is a normalization condition for the eigenvector. In this case, Eq. (24) can be simplified as

(30)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}={\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\mathrm{\Lambda })\delta (\mathrm{\Phi })|\nabla \mathrm{\Phi }|{\mathit{V}}_{\mathrm{n}}\mathrm{d}\mathrm{\Omega }},$

where $F(\mathit{u})=D_{ijkl}\mathit{\epsilon }(\mathit{u})\mathit{\epsilon }(\mathit{u})-{\lambda }_{\mathrm{1}}\rho \mathit{u}\mathit{u}$. By substituting Eqs. (11) and (15) into the preceding shape derivative equation, we obtain

(31)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial t}}=-{\displaystyle {\int }_{\mathrm{\Omega }}(F(\mathit{u})+\Lambda )\delta (\mathrm{\Phi })\mathit{R}{\displaystyle \frac{\partial s}{\partial t}}\mathrm{d}\mathrm{\Omega }}.$

Given that

(32)${\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}+\mathrm{\Lambda }{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}={\displaystyle \frac{\partial L(\mathit{u},\mathrm{\Phi })}{\partial s}}{\displaystyle \frac{\partial s}{\partial t}},$

the sensitivity of the objective function and volume constraint with respect to the design variables is respectively obtained as follows:

(33)${\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}F(\mathit{u})\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }},$

(34)${\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s}}=-{\displaystyle {\int }_{\mathrm{\Omega }}\delta (\mathrm{\Phi })\mathit{R}\mathrm{d}\mathrm{\Omega }}.$

Numerical implementation

Many design variables, which correspond to large-scale nonlinear equations in the eigenvalue problem, exist in continuum structural TO. Thus, OC is introduced to solve this eigenvalue TO problem. By properly iterating and updating the design variables, this optimization problem is guaranteed to converge to a final solution. Starting from an initial value, the iterative formula of the design variables is expressed as

(35)$s_i^{(k+\mathrm{1})}=c_i^{(k)}{s}_i^{(k)}.$

Theoretically, the iteration coefficients $c_i^{(k)}$ are obtained by setting Eq. (32) equal to 0, which can be written as

(36)$c_i^{(k)}=-{\displaystyle \frac{\partial J(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}/\max \left\{\mu ,{\mathrm{\Lambda }}^{(k)}{\displaystyle \frac{\partial V(\mathit{u},\mathrm{\Phi })}{\partial s_i^{(k)}}}\right\},$

where $\mu$ is a very small number that can avoid singularity when the sensitivity of the volume constraint with respect to the design variables is equal to 0. The Lagrangian multiplier $\mathrm{\Lambda }$ is calculated by the bisection method [76].

A flowchart of the structural TO for maximization of the first eigenfrequency problem is depicted in Fig. 2. Given the condition of constraints, when the relative difference value of the objective function in the current and previous iterations is less than 10⁻³, this optimization process is considered convergent, and the current optimization process is terminated.

Fig.2 Flowchart of the optimization procedure

Full size|PPT slide

Numerical examples

In this section, the proposed IGA-based level set TO framework is applied to two 2D optimization problems. For all examples, the properties of the isotropic material are set as follows: Young’s modulus $E=\mathrm{210}\mathrm{GPa}$and mass density $\rho =\mathrm{7.8}\times {\mathrm{10}}^{\mathrm{3}}{\mathrm{kg/m}}^3$. The properties of the artificial weak material are $E_{\mathrm{0}}=\mathrm{210}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{GPa}$ and mass density $\rho =\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{kg/m}}^3$. Poisson’s ratio $\nu =\mathrm{0.3}$ and plane stress state are assumed for all the materials. The two examples are TO of maximizing the fundamental eigenfrequency of the plane structure with a unit thickness of $\mathrm{0.001}\mathrm{m}$ and a prescribed material volume fraction of $\alpha =\mathrm{50}\%$. In the initial design, the available material is uniformly distributed over the entire admissible design domain. In the following examples, the boundary conditions are imposed by the collocation method, which enforces these conditions to be satisfied pointwise. Given that NURBS basis functions associated with the interior control points vanish at the structural boundary when open-knot vectors are employed, the displacement boundary condition applied on the left and right sides of the beam is imposed by setting the displacement values at left and right boundary control points to zero. All results are produced with programs developed in the MATLAB R2018a environment on a computer with an Intel Core i3-3240 CPU, 3.4 GHz clock speed, and 6 GB RAM. Additional details on the implementation of IGA on the MATLAB platform were presented in Ref. [77].

Cantilever beam

The first numerical example of a short cantilever beam for maximizing the first eigenfrequency optimization problem is shown in Fig. 3. The entire design domain is a rectangle with a size of 0.2 m×0.1 m, a Dirichlet boundary, and fixed displacement at the left edge of the design domain. A concentrated nonstructural mass $M=\mathrm{15.6}$ kg that is one-tenth of the total structural mass of the plate is placed at the center of the right side. Notably, the structure disappears without the nonstructural mass because no structure leads to the highest eigenfrequencies.

Fig.3 Design domain of a cantilever beam structure

Full size|PPT slide

In this example, IGA- and FEA-based TO approaches are applied to a similar problem for comparison. For this 2D example, equally spaced open-knot vectors are used for $x$ and $y$ directions, and the degrees of NURBS basis functions on the two directions are the same, i.e., $p=q=\mathrm{2}$. Nine-node quadratic rectangle elements are used for FEA, and the element number of both methods is the same, which facilitates a comparison under identical conditions.

The computation consumption times of the two methods are shown in Table 1. Notably, for simplicity in the proceeding table and figures, IGA-L and FEM-L represent IGA- and FEM-based LSM, respectively. The solution time of the system equation includes the time spent on assembling the stiffness matrix and solving the equation to obtain an objective function. In this case, when the numbers of Lagrange and NURBS elements are $a\times b$ ($a$ and $b$ represent the number of elements in $x$ and $y$ directions, respectively), the number of DOFs of the IGA-based method is $(a+\mathrm{2})(b+\mathrm{2})=ab+\mathrm{2}a+\mathrm{2}b+\mathrm{4}$ and that of FEM is $(\mathrm{2}a+\mathrm{1})(\mathrm{2}b+\mathrm{1})=\mathrm{4}ab+\mathrm{2}a+\mathrm{2}b+\mathrm{1}$. The DOF of IGA is much less than that of FEM (nearly a quarter of FEM) when the element number is sufficiently large. Table 1 shows that the computation efficiency of IGA-based LSM is higher than that of FEM-based LSM in TO due to the fewer DOFs or smaller size of equations in the IGA-based approach. However, because of the extra calculations of the basis functions and their derivative in the IGA-based approach, the ratio between the iteration time and DOFs of the IGA-based optimization method is higher than that of the FEM-based method.

Tab.1 Comparison of IGA- and FEA-based LSM TO

Method	Number of elements	Number of DOFs	Time of each iteration/s	Time of solution of system equation/s
IGA-L	256×128	67080	873.26	713.68
FEA-L	256×128	263682	1580.42	1237.16
IGA-L	128×64	17160	30.06	25.31
FEA-L	128×64	66306	54.17	43.97
IGA-L	64×32	4488	6.41	5.18
FEA-L	64×32	16770	8.55	6.85

The optimal layouts obtained by using IGA and FEA with 128×64 elements are shown in Fig. 4. In this case, the results are similar, and the crisp boundary is obtained due to the level set method. Adopting B-spline basis functions to parameterize the level set function together with IGA for calculation instead of FEA leads to the rapid convergence of the IGA-based optimization process.

Fig.4 Optimized results obtained by using IGA and FEA methods with 128 × 64 elements

Full size|PPT slide

The convergence history of optimization using IGA and FEA with 128 × 64 elements is shown in Figs. 5 and 6. Figure 5 shows the convergence history of the first eigenfrequency and volume ratio of the structure by using IGA- and FEA-based optimization frameworks, respectively. The initial designs and resultant structures of both methods are basically the same. In optimization with the IGA method, the ${\omega }_i$ of the initial design and the resultant optimum are 173.5 and 188.7, respectively. The volume ratio of the initial design and the resultant optimum are 0.79 and 0.5, respectively. Thus, the fundamental eigenfrequency increased by 8.8% and the volume decreased by 36.7%. Figure 6 shows the iteration history of the first three eigenfrequencies by using both methods. The first eigenfrequency always remains simple, whereas the second and third eigenfrequencies tend to oscillate and overlap. By using the MAC method, problems regarding the orders of eigenfrequency exchange during the optimization process are avoided. The first eigenfrequency generally increases, and the second and third eigenfrequencies decrease as the volume ratio decreases.

Fig.5 Comparison of IGA and FEA in terms of convergence history

Full size|PPT slide

Fig.6 Comparison of IGA and FEA in terms of eigenfrequency history

Full size|PPT slide

Beam with clamped ends

In this section, we present an example of maximizing the fundamental frequency of a clamped beam structure shown in Fig. 7. The working domain has a size of 0.4 m×0.1 m. A fixed displacement boundary condition is imposed on both sides, and a concentrated nonstructural mass M=31.2 kg is placed at the center of structure. In this example, the capability of the proposed method to capture the optimum topology and the effect of the number of elements are studied. For this purpose, a clamped beam is solved with three mesh sizes, namely, 64×16, 128×32, and 256×64. The degree of NURBS basis function is 2 in both directions. The resulting layouts shown in Fig. 8 indicate that the mesh dependency problem is avoided in the proposed method. This figure also shows that an inaccurate optimum topology with a rough boundary is obtained as a consequence of coarse meshes. By refining the mesh, the smoothness of the boundaries of the optimal layout is improved. However, when the number of elements is larger than 128×32, the resulting layout slightly changes. Accounting for the computation cost and precision of results, 128×32 meshes are used in this example.

Fig.7 Design domain of a clamped beam structure

Full size|PPT slide

Fig.8 Optimal layouts obtained by using (a) 64×16 meshes, (b) 128×32 meshes, and (c) 256×64 meshes

Full size|PPT slide

The convergence history of the objective function and volume ratio is given in Fig. 9. The first eigenfrequency of the initial design and the resultant topology are 244.3 and 257.4, respectively. The volume ratio of the initial design and the resultant topology are 0.82 and 0.5, respectively. Thus, the fundamental frequency increases by 5.4%, and the volume decreases by 39%. Figure 10 shows the convergence history of the first three eigenfrequencies. The fundamental frequency remains simple throughout the entire optimization process. The value of first eigenfrequency gradually increases after the fifth iteration, whereas that of the second and third eigenfrequencies decline with the increase in iterations.

Fig.9 Convergence history

Full size|PPT slide

Fig.10 Iteration history of the first three eigenfrequencies

Full size|PPT slide

Conclusions

We solved maximum fundamental eigenfrequency TO problems with a level set model based on the IGA technique. IGA combines the fundamental idea of FEM with the spline technique from a computer-aided geometry design for the integration of CAD and CAE. The IGA method was also introduced to TO due to its superiority over currently used FEM in terms of accuracy and efficiency. The feature of the proposed method is the combination of IGA and LSM in eigenfrequency optimization where the same basis functions (NURBS) are used for geometry representation, dynamic analysis, and parameterization of the implicit LSF. High accuracy and smoothness of LSF were achieved by using smooth NURBS basis functions to approximate LSF. In the case of maximizing the fundamental eigenfrequency, a regularization method for the repetition or exchange of eigenfrequencies was employed to guarantee the simple behavior of structural eigenfrequency.

Two benchmark numerical examples of TO for dynamic problems were applied to verify the validity of the proposed approach. The results obtained from the comparison of FEA- and IGA-based level set TO methods demonstrated that the proposed IGA-based optimization method has better computational efficiency and converges faster than the traditional FEA-based optimization method. The results also showed that solving dynamic TO problems by using IGA together with the level set method is possible. Although we only presented examples of 2D structures, no theoretical difficulties will be encountered if the proposed is extended to the optimization of 3D structures.

References

[1] Ijspeert A J. Biorobotics: using robots to emulate and investigate agile locomotion. Science, 2014, 346( 6206): 196– 203

[2] Shi D, Zhang W X, Zhang W, Ding X L. A review on lower limb rehabilitation exoskeleton robots. Chinese Journal of Mechanical Engineering, 2019, 32( 1): 74

[3] Dollar A M, Herr H. Lower extremity exoskeletons and active orthoses: challenges and state-of-the-art. IEEE Transactions on Robotics, 2008, 24( 1): 144– 158

[4] van Kammen K, Boonstra A M, van der Woude L H V, Visscher C, Reinders-Messelink H A, den Otter R. Lokomat guided gait in hemiparetic stroke patients: the effects of training parameters on muscle activity and temporal symmetry. Disability and Rehabilitation, 2020, 42( 21): 2977– 2985

[5] Hidayah R, Bishop L, Jin X, Chamarthy S, Stein J, Agrawal S K. Gait adaptation using a cable-driven active leg exoskeleton (C-ALEX) with post-stroke participants. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2020, 28( 9): 1984– 1993

[6] Meuleman J, van Asseldonk E, van Oort G, Rietman H, van der Kooij H. LOPES II—design and evaluation of an admittance controlled gait training robot with shadow-leg approach. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2016, 24( 3): 352– 363

[7] Huang R, Cheng H, Qiu J, Zhang J W. Learning physical human–robot interaction with coupled cooperative primitives for a lower exoskeleton. IEEE Transactions on Automation Science and Engineering, 2019, 16( 4): 1566– 1574

[8] Zhou L B, Chen W H, Wang J H, Bai S P, Yu H Y, Zhang Y P. A novel precision measuring parallel mechanism for the closed-loop control of a biologically inspired lower limb exoskeleton. IEEE/ASME Transactions on Mechatronics, 2018, 23( 6): 2693– 2703

[9] Shi D, Zhang W X, Zhang W, Ju L H, Ding X L. Human-centred adaptive control of lower limb rehabilitation robot based on human–robot interaction dynamic model. Mechanism and Machine Theory, 2021, 162: 104340

[10] Long Y, Du Z J, Chen C F, Wang W D, He L, Mao X W, Xu G Q, Zhao G Y, Li X Q, Dong W. Development and analysis of an electrically actuated lower extremity assistive exoskeleton. Journal of Bionics Engineering, 2017, 14( 2): 272– 283

[11] Wei D, Li Z J, Wei Q, Su H, Song B, He W, Li J Q. Human-in-the-loop control strategy of unilateral exoskeleton robots for gait rehabilitation. IEEE Transactions on Cognitive and Developmental Systems, 2021, 13( 1): 57– 66

[12] Ding Y, Kim M, Kuindersma S, Walsh C J. Human-in-the-loop optimization of hip assistance with a soft exosuit during walking. Science Robotics, 2018, 3( 15): eaa r5438

[13] Ding H, Yang X J, Zheng N N, Li M, Lai Y N, Wu H. Tri-co robot: a Chinese robotic research initiative for enhanced robot interaction capabilities. National Science Review, 2018, 5( 6): 799– 801

[14] Meng W, Liu Q, Zhou Z D, Ai Q S, Sheng B, Xie S Q. Recent development of mechanisms and control strategies for robot-assisted lower limb rehabilitation. Mechatronics, 2015, 31: 132– 145

[15] Kalita B, Narayan J, Dwivedy S K. Development of active lower limb robotic-based orthosis and exoskeleton devices: a systematic review. International Journal of Social Robotics, 2021, 13( 4): 775– 793

[16] Zhou J M, Yang S, Xue Q. Lower limb rehabilitation exoskeleton robot: a review. Advances in Mechanical Engineering, 2021, 13( 4): 16878140211011862

[17] Yan T F, Cempini M, Oddo C M, Vitiello N. Review of assistive strategies in powered lower-limb orthoses and exoskeletons. Robotics and Autonomous Systems, 2015, 64: 120– 136

[18] Baud R, Manzoori A R, Ijspeert A, Bouri M. Review of control strategies for lower-limb exoskeletons to assist gait. Journal of NeuroEngineering and Rehabilitation, 2021, 18( 1): 119

[19] Ferris D P, Sawicki G S, Daley M A. A physiologist’s perspective on robotic exoskeletons for human locomotion. International Journal of Humanoid Robotics, 2007, 4( 3): 507– 528

[20] Bohannon R W, Schaubert K. Long-term reliability of the timed up-and-go test among community-dwelling elders. Journal of Physical Therapy Science, 2005, 17( 2): 93– 96

[21] Podsiadlo D, Richardson S. The timed “up & go”: a test of basic functional mobility for frail elderly persons. Journal of the American Geriatrics Society, 1991, 39( 2): 142– 148

[22] Guralnik J M, Simonsick E M, Ferrucci L, Glynn R J, Berkman L F, Blazer D G, Scherr P A, Wallace R B. A short physical performance battery assessing lower extremity function: association with self-reported disability and prediction of mortality and nursing home admission. Journal of Gerontology, 1994, 49( 2): M85– M94

[23] Berg K, Wood-Dauphine S, Williams J I, Gayton D. Measuring balance in the elderly: preliminary development of an instrument. Physiotherapy Canada, 1989, 41( 6): 304– 311

[24] Ganz D A, Bao Y R, Shekelle P G, Rubenstein L Z. Will my patient fall? Journal of the American Medical Association, 2007, 297( 1): 77– 86

[25] Lee T K M, Belkhatir M, Sanei S. A comprehensive review of past and present vision-based techniques for gait recognition. Multimedia Tools and Applications, 2014, 72( 3): 2833– 2869

[26] Scheffer C, Cloete T. Inertial motion capture in conjunction with an artificial neural network can differentiate the gait patterns of hemiparetic stroke patients compared with able-bodied counterparts. Computer Methods in Biomechanics and Biomedical Engineering, 2012, 15( 3): 285– 294

[27] Den Otter A R, Geurts A C H, Mulder T, Duysens J. Abnormalities in the temporal patterning of lower extremity muscle activity in hemiparetic gait. Gait & Posture, 2007, 25( 3): 342– 352

[28] Pickle N T, Shearin S M, Fey N P. Dynamic neural network approach to targeted balance assessment of individuals with and without neurological disease during non-steady-state locomotion. Journal of NeuroEngineering and Rehabilitation, 2019, 16( 1): 88

[29] Swinnen E, Beckwée D, Meeusen R, Baeyens J P, Kerckhofs E. Does robot-assisted gait rehabilitation improve balance in stroke patients? A systematic review. Topics in Stroke Rehabilitation, 2014, 21( 2): 87– 100

[30] Park J H, Kim Y, Lee K J, Yoon Y S, Kang S H, Kim H, Park H S. Artificial neural network learns clinical assessment of spasticity in modified ashworth scale. Archives of Physical Medicine and Rehabilitation, 2019, 100( 10): 1907– 1915

[31] Pinto-Fernandez D, Torricelli D, Sanchez-Villamanan M D C, Aller F, Mombaur K, Conti R, Vitiello N, Moreno J C, Pons J L. Performance evaluation of lower limb exoskeletons: a systematic review. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2020, 28( 7): 1573– 1583

[32] Galen S S, Clarke C J, Allan D B, Conway B A. A portable gait assessment tool to record temporal gait parameters in SCI. Medical Engineering & Physics, 2011, 33( 5): 626– 632

[33] Granat M H, Maxwell D J, Bosch C J, Ferguson A C B, Lees K R, Barbenel J C. A body-worn gait analysis system for evaluating hemiplegic gait. Medical Engineering & Physics, 1995, 17( 5): 390– 394

[34] Neckel N, Pelliccio M, Nichols D, Hidler J. Quantification of functional weakness and abnormal synergy patterns in the lower limb of individuals with chronic stroke. Journal of NeuroEngineering and Rehabilitation, 2006, 3( 1): 17

[35] Neckel N D, Blonien N, Nichols D, Hidler J. Abnormal joint torque patterns exhibited by chronic stroke subjects while walking with a prescribed physiological gait pattern. Journal of NeuroEngineering and Rehabilitation, 2008, 5( 1): 19

[36] Hu B H, Zhang X F, Mu J S, Wu M, Zhu Z J, Liu Z S, Wang Y. Spasticity measurement based on the HHT marginal spectrum entropy of sEMG using a portable system: a preliminary study. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2018, 26( 7): 1424‒ 1434

[37] Dewald J P A, Pope P S, Given J D, Buchanan T S, Rymer W Z. Abnormal muscle coactivation patterns during isometric torque generation at the elbow and shoulder in hemiparetic subjects. Brain, 1995, 118( 2): 495‒ 510

[38] Shestakov M P. Balance of a multijoint biomechanical system in natural and artificial environments: a simulation model. Journal of Physiological Anthropology, 2007, 26( 3): 419– 423

[39] Kasaoka K, Sankai Y. Predictive control estimating operator’s intention for stepping-up motion by exo-skeleton type power assist system HAL. In: Proceedings of 2001 IEEE/RSJ International Conference on Intelligent Robots & Systems. Maui: IEEE, 2001, 3: 1578– 1583

[40] Pei P, Pei Z C, Tang Z Y, Gu H. Position tracking control of PMSM based on fuzzy PID-variable structure adaptive control. Mathematical Problems in Engineering, 2018, ( 1): 5794067

[41] Liu D F, Tang Z Y, Pei Z C. Variable structure compensation PID control of asymmetrical hydraulic cylinder trajectory tracking. Mathematical Problems in Engineering, 2015, ( 1): 890704

[42] Zhang M M, Xie S Q, Li X L, Zhu G L, Meng W, Huang X L, Veale A J. Adaptive patient-cooperative control of a compliant ankle rehabilitation robot (CARR) with enhanced training safety. IEEE Transactions on Industrial Electronics, 2018, 65( 2): 1398– 1407

[43] Shao Y X, Zhang W X, Su Y J, Ding X L. Design and optimisation of load-adaptive actuator with variable stiffness for compact ankle exoskeleton. Mechanism and Machine Theory, 2021, 161: 104323

[44] Yu H Y, Huang S N, Chen G, Pan Y P, Guo Z. Human–robot interaction control of rehabilitation robots with series elastic actuators. IEEE Transactions on Robotics, 2015, 31( 5): 1089– 1100

[45] Zhang W, Zhang W X, Shi D, Ding X L. Design of hip joint assistant asymmetric parallel mechanism and optimization of singularity-free workspace. Mechanism and Machine Theory, 2018, 122: 389– 403

[46] Wang D H, Lee K M, Ji J J. A passive gait-based weight-support lower extremity exoskeleton with compliant joints. IEEE Transactions on Robotics, 2016, 32( 4): 933– 942

[47] Wang Y J, Wu C L, Yu L Q, Mei Y Y. Dynamics of a rolling robot of closed five-arc-shaped-bar linkage. Mechanism and Machine Theory, 2018, 121: 75– 91

[48] Bascetta L, Ferretti G, Scaglioni B. Closed form Newton–Euler dynamic model of flexible manipulators. Robotica, 2017, 35( 5): 1006– 1030

[49] Sun Z B, Li F, Duan X Q, Jin L, Lian Y F, Liu S S, Liu K P. A novel adaptive iterative learning control approach and human-in-the-loop control pattern for lower limb rehabilitation robot in disturbances environment. Autonomous Robots, 2021, 45( 4): 595– 610

[50] Zoss A, Kazerooni H. Design of an electrically actuated lower extremity exoskeleton. Advanced Robotics, 2006, 20( 9): 967– 988

[51] Sun W, Lin J W, Su S F, Wang N, Er M J. Reduced adaptive fuzzy decoupling control for lower limb exoskeleton. IEEE Transactions on Cybernetics, 2021, 51( 3): 1099– 1109

[52] Qiu S Y, Guo W, Caldwell D, Chen F. Exoskeleton online learning and estimation of human walking intention based on dynamical movement primitives. IEEE Transactions on Cognitive and Developmental Systems, 2021, 13( 1): 67– 79

[53] Ghan J, Steger R, Kazerooni H. Control and system identification for the Berkeley lower extremity exoskeleton (BLEEX). Advanced Robotics, 2006, 20( 9): 989– 1014

[54] Huang R, Cheng H, Guo H L, Lin X C, Zhang J W. Hierarchical learning control with physical human-exoskeleton interaction. Information Sciences, 2018, 432: 584– 595

[55] Ruiz Garate V, Parri A, Yan T F, Munih M, Molino Lova R, Vitiello N, Ronsse R. Walking assistance using artificial primitives: a novel bioinspired framework using motor primitives for locomotion assistance through a wearable cooperative exoskeleton. IEEE Robotics & Automation Magazine, 2016, 23( 1): 83– 95

[56] Xu J J, Li Y F, Xu L S, Peng C, Chen S Q, Liu J F, Xu C C, Cheng G X, Xu H, Liu Y, Chen J. A m ulti-mode rehabilitation robot with magnetorheological actuators based on human motion intention estimation. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2019, 27( 10): 2216– 2228

[57] Meijneke C, van Oort G, Sluiter V, van Asseldonk E, Tagliamonte N L, Tamburella F, Pisotta I, Masciullo M, Arquilla M, Molinari M, Wu A R, Dzeladini F, Ijspeert A J, van der Kooij H. Symbitron exoskeleton: design, control, and evaluation of a modular exoskeleton for incomplete and complete spinal cord injured individuals. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2021, 29: 330– 339

[58] Song S, Geyer H. A neural circuitry that emphasizes spinal feedback generates diverse behaviours of human locomotion. The Journal of Physiology, 2015, 593( 16): 3493– 3511

[59] Aguirre-Ollinger G, Nagarajan U, Goswami A. An admittance shaping controller for exoskeleton assistance of the lower extremities. Autonomous Robots, 2016, 40( 4): 701– 728

[60] Kazerooni H, Steger R, Huang L H. Hybrid control of the Berkeley lower extremity exoskeleton (BLEEX). International Journal of Robotics Research, 2006, 25( 5–6): 561– 573

[61] He W, Li Z J, Chen C L P. A survey of human-centered intelligent robots: issues and challenges. IEEE/CAA Journal of Automatica Sinica, 2017, 4( 4): 602– 609

[62] Jamwal P K, Xie S Q, Hussain S, Parsons J G. An adaptive wearable parallel robot for the treatment of ankle injuries. IEEE/ASME Transactions on Mechatronics, 2014, 19( 1): 64– 75

[63] Karavas N, Ajoudani A, Tsagarakis N, Saglia J, Bicchi A, Caldwell D. Tele-impedance based assistive control for a compliant knee exoskeleton. Robotics and Autonomous Systems, 2015, 73: 78– 90

[64] Kao P C, Lewis C L, Ferris D P. Invariant ankle moment patterns when walking with and without a robotic ankle exoskeleton. Journal of Biomechanics, 2010, 43( 2): 203– 209

[65] Kilicarslan A, Grossman R G, Contreras-Vidal J L. A robust adaptive denoising framework for real-time artifact removal in scalp EEG measurements. Journal of Neural Engineering, 2016, 13( 2): 0 26013

[66] Bulea T C, Prasad S, Kilicarslan A, Contreras-Vidal J L. Sitting and standing intention can be decoded from scalp EEG recorded prior to movement execution. Frontiers in Neuroscience, 2014, 8: 376

[67] Liu D, Chen W H, Pei Z C, Wang J H. A brain-controlled lower-limb exoskeleton for human gait training. Review of Scientific Instruments, 2017, 88( 10): 104302

[68] Lyu M X, Chen W H, Ding X L, Wang J H, Pei Z C, Zhang B C. Development of an EMG-controlled knee exoskeleton to assist home rehabilitation in a game context. Frontiers in Neurorobotics, 2019, 13: 67

[69] Huang L P, Zheng J B, Hu H C. Online gait phase detection in complex environment based on distance and multi-sensors information fusion using inertial measurement units. International Journal of Social Robotics, 2022, 14( 2): 413– 428

[70] Kang I, Molinaro D D, Duggal S, Chen Y R, Kunapuli P, Young A J. Real-time gait phase estimation for robotic hip exoskeleton control during multimodal locomotion. IEEE Robotics and Automation Letters, 2021, 6( 2): 3491– 3497

[71] Wang J B, Fei Y Q, Chen W D. Integration, sensing, and control of a modular soft-rigid pneumatic lower limb exoskeleton. Soft Robotics, 2020, 7( 2): 140– 154

[72] Seel T, Raisch J, Schauer T. IMU-based joint angle measurement for gait analysis. Sensors, 2014, 14( 4): 6891– 6909

[73] Beravs T, Reberšek P, Novak D, Podobnik J, Munih M. Development and validation of a wearable inertial measurement system for use with lower limb exoskeletons. In: Proceedings of 2011 the 11th IEEE-RAS International Conference on Humanoid Robots. Bled: IEEE, 2011, 212– 217

[74] Ji J C, Song T, Guo S, Xi F F, Wu H. Robotic-assisted rehabilitation trainer improves balance function in stroke survivors. IEEE Transactions on Cognitive and Developmental Systems, 2020, 12( 1): 43– 53

[75] Chen Z L, Guo Q, Xiong H Y, Jiang D, Yan Y. Control and implementation of 2-DOF lower limb exoskeleton experiment platform. Chinese Journal of Mechanical Engineering, 2021, 34( 1): 22

[76] Chen B J, Zheng E H, Fan X D, Liang T, Wang Q N, Wei K L, Wang L. Locomotion mode classification using a wearable capacitive sensing system. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2013, 21( 5): 744– 755

[77] Casas J, Senft E, Gutiérrez L F, Rincón-Rocancio M, Múnera M, Belpaeme T, Cifuentes C A. Social assistive robots: assessing the impact of a training assistant robot in cardiac rehabilitation. International Journal of Social Robotics, 2021, 13( 6): 1189– 1203

[78] Billinger S A, Arena R, Bernhardt J, Eng J J, Franklin B A, Johnson C M, MacKay-Lyons M, Macko R F, Mead G E, Roth E J, Shaughnessy M, Tang A. Physical activity and exercise recommendations for stroke survivors: a statement for healthcare professionals from the American Heart Association/American Stroke Association. Stroke, 2014, 45( 8): 2532– 2553

[79] Maggioni S, Melendez-Calderon A, van Asseldonk E, Klamroth-Marganska V, Lünenburger L, Riener R, van der Kooij H. Robot-aided assessment of lower extremity functions: a review. Journal of NeuroEngineering and Rehabilitation, 2016, 13( 1): 72

[80] Hussain S, Xie S Q, Jamwal P K. Robust nonlinear control of an intrinsically compliant robotic gait training orthosis. IEEE Transactions on Systems, Man, and Cybernetics. Systems, 2013, 43( 3): 655– 665

[81] Yu X B, Li B, He W, Feng Y H, Cheng L, Silvestre C. Adaptive-constrained impedance control for human ‒ robot co-transportation. IEEE Transactions on Cybernetics, 2021 (in press)

[82] Shi D, Zhang W X, Zhang W, Ding X L. Assist-as-needed attitude control in three-dimensional space for robotic rehabilitation. Mechanism and Machine Theory, 2020, 154: 104044

[83] Shi D, Zhang W X, Zhang W, Ding X L. Force field control for the three-dimensional gait adaptation using a lower limb rehabilitation robot. In: Uhl T, ed. Advances in Mechanism and Machine Science. IFToMM WC 2019. Mechanisms and Machine Science, vol 73. Cham: Springer International Publishing, 2019, 73: 1919– 1928

[84] Wang L T, Wang S Q, van Asseldonk E H F, van der Kooij H. Actively controlled lateral gait assistance in a lower limb exoskeleton. In: Proceedings of 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems. Tokyo: IEEE, 2013, 965– 970

[85] Tsukahara A, Hasegawa Y, Eguchi K, Sankai Y. Restoration of gait for spinal cord injury patients using HAL with intention estimator for preferable swing speed. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2015, 23( 2): 308– 318

[86] Duschau-Wicke A, von Zitzewitz J, Caprez A, Lunenburger L, Riener R. Path control: a method for patient-cooperative robot-aided gait rehabilitation. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2010, 18( 1): 38– 48

[87] Riener R, Lunenburger L, Jezernik S, Anderschitz M, Colombo G, Dietz V. Patient-cooperative strategies for robot-aided treadmill training: first experimental results. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2005, 13( 3): 380– 394

[88] Niu X M, Gao G Q, Liu X J, Fang Z M. Decoupled sliding mode control for a novel 3-DOF parallel manipulator with actuation redundancy. International Journal of Advanced Robotic Systems, 2015, 12( 5): 64

[89] Mohanta J K, Santhakumar M, Kurtenbach S, Corves B, Hüsing M. Augmented PID control of a 2PPR-2PRP planar parallel manipulator for lower limb rehabilitation applications. In: Corves B, Lovasz E C, Hüsing M, Maniu I, Gruescu C, eds. New Advances in Mechanisms, Mechanical Transmissions and Robotics Mechanisms and Machine Science, vol 46. Cham: Springer International Publishing, 2017, 46: 391– 399

[90] Luo L C, Peng L, Wang C, Hou Z G. A greedy assist-as-needed controller for upper limb rehabilitation. IEEE Transactions on Neural Networks and Learning Systems, 2019, 30( 11): 3433– 3443

[91] Li Y, Ge S S. Human‒robot collaboration based on motion intention estimation. IEEE/ASME Transactions on Mechatronics, 2014, 19( 3): 1007– 1014

[92] Emken J L, Bobrow J E, Reinkensmeyer D J. Robotic movement training as an optimization problem: designing a controller that assists only as needed. In: Proceedings of the 9th International Conference on Rehabilitation Robotics (ICORR). Chicago: IEEE, 2005, 307– 312

[93] Zanotto D, Stegall P, Agrawal S K. 2014. Adaptive assist-as-needed controller to improve gait symmetry in robot-assisted gait training. In: Proceedings of 2014 IEEE International Conference on Robotics and Automation (ICRA). Hong Kong: IEEE, 2014, 724– 729

[94] Fineberg D B, Asselin P, Harel N Y, Agranova-Breyter I, Kornfeld S D, Bauman W A, Spungen A M. Vertical ground reaction force-based analysis of powered exoskeleton-assisted walking in persons with motor-complete paraplegia. The Journal of Spinal Cord Medicine, 2013, 36( 4): 313– 321

[95] Long Y, Du Z J, Wang W D, Dong W. Human motion intent learning based motion assistance control for a wearable exoskeleton. Robotics and Computer-Integrated Manufacturing, 2018, 49: 317– 327

[96] Shi D, Zhang W, Wang L D, Zhang W X, Feng Y G, Ding X L. Joint angle adaptive coordination control of a serial parallel lower limb rehabilitation exoskeleton. IEEE Transactions on Medical Robotics and Bionics, 2022, 4(3): 775– 784

[97] Huang R, Cheng H, Chen Y, Chen Q M, Lin X C, Qiu J. Optimisation of reference gait trajectory of a lower limb exoskeleton. International Journal of Social Robotics, 2016, 8( 2): 223– 235

[98] Liu D X, Wu X Y, Du W B, Wang C, Chen C J, Xu T T. Deep spatial-temporal model for rehabilitation gait: optimal trajectory generation for knee joint of lower-limb exoskeleton. Assembly Automation, 2017, 37( 3): 369– 378

[99] Kwakkel G, Kollen B, Lindeman E. Understanding the pattern of functional recovery after stroke: facts and theories. Restorative Neurology and Neuroscience, 2004, 22( 3–5): 281– 299

[100] Langhorne P, Bernhardt J, Kwakkel G. Stroke rehabilitation. The Lancet, 2011, 377( 9778): 1693– 1702

[101] Barbeau H, Ladouceur M, Mirbagheri M M, Kearney R E. The effect of locomotor training combined with functional electrical stimulation in chronic spinal cord injured subjects: walking and reflex studies. Brain Research Reviews, 2002, 40( 1–3): 274– 291

[102] Yang Y R, Wang R Y, Lin K H, Chu M Y, Chan R C. Task-oriented progressive resistance strength training improves muscle strength and functional performance in individuals with stroke. Clinical Rehabilitation, 2006, 20( 10): 860– 870

[103] Morawietz C, Moffat F. Effects of locomotor training after incomplete spinal cord injury: a systematic review. Archives of Physical Medicine and Rehabilitation, 2013, 94( 11): 2297– 2308

[104] Herbert R D, Taylor J L, Lord S R, Gandevia S C. Prevalence of motor impairment in residents of New South Wales, Australia aged 55 years and over: cross-sectional survey of the 45 and Up cohort. BMC Public Health, 2020, 20( 1): 1353

[105] Chen F X, Zhang C, Chen J H, Yang G L. Accurate subdomain model for computing magnetic field of short moving-magnet linear motor with H albach array. IEEE Transactions on Magnetics, 2020, 56( 9): 1– 9

[106] Fang Y, Lerner Z F. Feasibility of augmenting ankle exoskeleton walking performance with step length biofeedback in individuals with cerebral palsy. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2021, 29: 442– 449

[107] Kim Y, Chortos A, Xu W T, Liu Y X, Oh J Y, Son D, Kang J, Foudeh A M, Zhu C X, Lee Y, Niu S M, Liu J, Pfattner R, Bao Z, Lee T W. A bioinspired flexible organic artificial afferent nerve. Science, 2018, 360( 6392): 998– 1003

[108] Huang Y C, He Z X, Liu Y X, Yang R Y, Zhang X F, Cheng G, Yi J G, Ferreira J P, Liu T. Real-time intended knee joint motion prediction by deep-recurrent neural networks. IEEE Sensors Journal, 2019, 19( 23): 11503– 11509

[109] Ugartemendia A, Rosquete D, Gil J J, Diaz I, Borro D. Machine learning for active gravity compensation in robotics: application to neurological rehabilitation systems. IEEE Robotics & Automation Magazine, 2020, 27( 2): 78– 86

[110] Fang W, An Z W. A scalable wearable AR system for manual order picking based on warehouse floor-related navigation. The International Journal of Advanced Manufacturing Technology, 2020, 109( 7–8): 2023– 2037

[111] Oppezzo M, Schwartz D L. Give your ideas some legs: the positive effect of walking on creative thinking. Journal of Experimental Psychology: Learning, Memory, and Cognition, 2014, 40( 4): 1142– 1152

[112] Hidayah R, Chamarthy S, Shah A, Fitzgerald-Maguire M, Agrawal S K. Walking with augmented reality: a preliminary assessment of visual feedback with a cable-driven active leg exoskeleton (C-ALEX). IEEE Robotics and Automation Letters, 2019, 4( 4): 3948– 3954

Acknowledgements

This paper was funded by the National Natural Science Foundation of China (Grant Nos. 91848104, 91748201, and 52105004). The authors thank Yushuang Duan and Hongqian Zhang for their contributions to this study.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format as long as appropriate credit is given to the original author(s) and source, a link to the Creative Commons license is provided, and the changes made are indicated.

The images or other third-party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Visit http://creativecommons.org/licenses/by/4.0/ to view a copy of this license.

RIGHTS & PERMISSIONS

2022 The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn