Design of a Wrist Rehabilitation System with a Novel Mixed Structural Optimization Applying Improved Harmony Search

Abstract

This paper presents the development of a wrist rehabilitation system with a novel approach for structural design, based on the modeling of an optimization problem solved by a metaheuristic algorithm, Improved Harmony Search (ImHS). It is part of a project for developing low-cost rehabilitation systems expressly designed for the population of Latin American countries. A mixed optimization problem is modeled for the design, where the material type is associated with an integer variable and the dimensions of the components are continuous parameters. The novelty is that each element is calculated individually, but considering the combined effect over the structure. The optimization works simultaneously on both the material selection and the meeting of the associated constraints, to guarantee that the system will not fail because of any load, neither it will be unsafe for the patients, since the operation will always be within the limits considered in the modeling. ImHS is a variant of the Harmony Search algorithm, modified to enhance the exploration and exploitation processes. It is a simple yet powerful metaheuristic, implemented in this development with additional modifications to handle constraints and mixed variables. The proposed approach produced quality results, indicating that ImHS can be applied to solve complex engineering problems, facilitating the manufacture and control processes.

1. Introduction

Cerebral strokes are the main cause of neurological impairments worldwide [1], with a survival rate of 60–75%. About 80–90% of the survivors suffer from a reduction or complete loss of their psycho-motor functions [2]. Other factors such as accidents or different diseases also affect these functions. The mobility effect negatively impacts life quality and Daily Living Activities (DLA) such as eating, grooming or dressing. There is a great necessity of rehabilitation for these patients but most classical methods are difficult to implement, because they demand a high physical effort from the therapist [3]. Robotic systems are a good option since they support the high intensity, steady performance and long activity periods required [4]. The wrist is one of the most exposed joints to trauma, injuries or strokes, and diverse devices have been developed for its rehabilitation. An exoskeleton with elastic actuators was developed in [5]. In [6], a rehabilitation system with force sensing was built. A structure for the 3-DOF movements of the wrist by either arm was presented in [7]. In [8], a training glove with pneumatic actuators is designed, and in [9] a 3-DOF exoskeleton is presented.

However, none of the aforementioned works includes the detailed design of the structural component. Most developments are focused on the system control, and generally the devices are over-dimensioned and difficult to produce. The development of computing tools has accelerated the design process, simplifying the manufacture and control while ensuring a safer operation. Few works consider these aspects, highlighting the development of Gandomi et al. [10], with a structural design carried out by metaheuristics, one element at a time. Several methodologies have been proposed, such as Concurrent Design (CD) [11], Design For Control (DFC) [12], Design For A (A = life cycle, manufacture, etc.) [13], and Multidisciplinary Design Optimization (MDO) [14]. Most structural optimization cases have both non-linear objective functions and constraints, with limited search spaces, making it impossible to solve the associated optimization problems by classic techniques, so metaheuristics are a useful tool for their solution. Many works have addressed aspects as the weight optimization to maximize stiffness [15,16], or the topology optimization of beams for different stress types [17,18].

Harmony Search (HS) is a metaheuristic developed by Geem et al. [19] on the basis that a harmony is a certain sound combination and the search for the best harmony is analogous to finding a solution for an optimization problem [20]. It presents a good balance between exploration and exploitation, an easy and low-cost computational implementation with few mathematical requirements, and fast convergence [21]. It has been applied in fields such as engineering [22,23,24,25,26], computing [27,28], image processing [29,30], and electronics [31]. In relation to the scope of this work, HS has been successfully used in different medical applications. In [32], it is applied to select the adequate features for a digital hearing aid device, after classifying the acoustic environment where the user is immerse. A method called MRMR-COA-HS is proposed in [33] for a two-stage gene selection, with a hybrid combination of the Cuckoo Optimization Algorithm (COA) and HS for cancer classification. In [34], HS is used with the entropy of Kapur in a method for image segmentation by thresholding, to detect skin cancer in early stages. Another segmentation technique based on HS is presented in [35], for analyzing magnetic resonance images (MRI). In [36], the authors developed an application for medical physics, using HS to optimize the radioisotope placement of a high dose-rate brachytherapy for prostate cancer. There are also diverse works on HS-based structural design. In [37,38], the efficiency of HS for solving this type of problems is analyzed, and a tool based on HS for the design of a three-span box-girder bridge is presented in [39].

In this work, the structural design of a wrist rehabilitation system is modeled as an optimization problem focused on the weight of the device [40], presenting an approach that simultaneously considers the manufacturing material selection and the meet of the constraints associated to the rehabilitation system. The novelty is that each element is calculated individually, but considering the combined effect over the structure. This paper is organized as follows: Section 2 presents the modeling of the rehabilitation system, describing its functionality and the component analysis while Section 3 explains the optimization method. Section 4 describes the optimization problem and the solution algorithm. The results are given in Section 5 and the final discussion in Section 6.

2. Problem Statement

The best strategy during the design process of mechanical systems is to identify their components [41] and revise them separately, but including the combined effect of all the elements to avoid structural failures. The von-Mises theory is useful for stress analysis when a combined axial-bending-shear load acts on a structural component. The criterion of Goodman is employed to analyze elements in a structure dynamically loaded under combined stress, performing the fatigue analysis by any of the currently accepted failure theories [42,43]. According to the distortion energy theory, the maximum strain is given by Equation (1), where ${\sigma }_f$, ${\sigma }_a$, and ${\tau }_s$ are the stress generated by bending, axial load, and shear, respectively; and $k_f=1$, $k_a=0.85$, and $k_t=0.59$ are the stress concentration factors for bending, axial load, and shear, respectively, as recommended in [42].

$${\sigma }^{\prime }={\left[{\left(k_f{\sigma }_f+k_a\frac{{\sigma }_a}{0.85}\right)}^2+3{\left(k_s{\tau }_s\right)}^2\right]}^{\frac{1}{2}}$$

(1)

Case Study

The wrist joint performs its movements along two axes, each corresponding to a DOF: Flexion/Extension (FL/EX) and Radial Deviation (RD)/Ulnar Deviation (UD). It is oriented by a third DOF, the Pronation (P)/Supination (S) movement of the elbow joint. Figure 1a shows the rehabilitation device proposed by Aponte-Rodriguez et al. [40] to reproduce these 3-DOF. It is part of a join project between the Instituto Politecnico Nacional (Mexico) and the Universidad Militar Nueva Granada (Colombia), to develop low-cost rehabilitation systems expressly designed for the population of Latin American countries. The specifications of the device components were determined from the motion ranges proposed in [1,44] and the hand anthropometry for the Latin American population [45]. The torques required to drive each articulation were considered for the structural analysis, taking into account both the requirements of the DLA and the values proposed in literature [1,3,46]. The RD/UD and FL/EX movements are actuated directly by servomotors connected to the structural components. For this reason, an adequate design is critical for a safe operation.

The development in this paper (Figure 1b) is a continuation of the aforementioned project, and consists in a modification to the extensor arm of the mechanism with two links instead of one to facilitate the wrist movements. The size of each link is determined by its length l, width W, thickness t, and diameter d (Figure 2a), and both of them form an open kinematic chain. The links should be as light as possible, to take advantage of the highest torque in the hand movement during the therapy. The solution method tackles the simultaneous design of more than one structural element, considering the interaction and joint effect while designing. Although this work only analyzes two links, all the mobile elements are designed with the same methodology. A complex mixed-optimization problem is obtained, with both the objective function and the constraints of non-linear type. The free-body diagrams (Figure 2b) consider the bending and torsional external loads for the combined stress. The bending stress is calculated for each link by Equation (2), with Equations (3) and (4) for the bending moments and the second moments of area, respectively. The torque is $T_{xy}=3.571$ Nm, and the forces are $F_{z1}=F_{z2}=25.3N$, and $F_{y1}=F_{y2}+$weight of link 2, and $F_{y2}=30.4N+EF$ [40], with an extra load $EF=49N$ for this case. The shear stress produced by torsion is calculated by Equation (5), with $W>t$.

$${\sigma }_f=\frac{M_y}{I_{yy}}+\frac{M_z}{I_{zz}}$$

(2)

$$\begin{array}{c}\hfill M_y=\frac{l\times W\times F_y}{2},M_z=\frac{l\times t\times F_z}{2}\end{array}$$

(3)

$$\begin{array}{c}\hfill I_{yy}=\left(\frac{1}{12}t\right)(W^3-d^3),I_{zz}=\left(\frac{1}{12}W-d\right)\left(t^3\right)\end{array}$$

(4)

$${\tau }_s=\frac{3T_{xy}}{8a{b}^2}\left(1+0.6095\left(\frac{t}{W}\right)+0.88665{\left(\frac{t}{W}\right)}^3+0.91{\left(\frac{t}{W}\right)}^4\right)$$

(5)

There is no load generating an axial force, so ${\sigma }_a=0$. The maximum total deflection is obtained using Equations (6) and (7), where E is the elastic modulus of the selected material:

$$\begin{array}{c}\hfill {\delta }_{max}=\sqrt{{\left(y_{xy}\right)}^2+{\left(y_{xz}\right)}^2}\end{array}$$

(6)

$$\begin{array}{c}\hfill y_{xy}=\frac{F_y{l}^3}{3E{I}_{yy}},y_{xz}=\frac{F_z{l}^3}{3E{I}_{zz}}\end{array}$$

(7)

3. Optimization Problem

In this work, Improved Harmony Search is used to solve a structural design case, tackling simultaneously the dimensional synthesis and the material selection of more than one element at a time, considering both the interaction and joint effect of the components. A mathematical problem for weight optimization is proposed to obtain the dimensions and material required for each element, meeting the design constraints. This is a mixed problem, because the material selection is a discrete variable, while the dimensions are continuous. The procedure is as follows: given a set of materials, a load state, the maximum allowable deflections and the minimum safety factor allowed, it is desired to find the width, thickness, and appropriate material to build each link.

3.1. Objective Function

It is required to optimize the weight of the links to implicitly optimize the power consumption and the costs associated to the structure. The objective function in Equation (9) corresponds to the total weight, including a sum term for each link, where l and d are the link length and diameter, respectively, and g is the gravitational acceleration. Expression (9) presents the design vector, where $x_1,x_2,$ and $x_3$ are the material type, thickness, and width of link 1, respectively, while $x_4,x_5$, and $x_6$ are the same parameters for link 2.

$$\begin{array}{c}\hfill minf\left(\overrightarrow{x}\right)=x_{1}g{x}_2\left(l{x}_3+\frac{\pi }{4}{x}_3^2-\frac{\pi }{2}{d}^2\right)+x_{4}g{x}_5\left(l{x}_6+\frac{\pi }{4}{x}_6^2-\frac{\pi }{2}{d}^2\right)\end{array}$$

(8)

$$\begin{array}{c}\hfill \overrightarrow{x}={[x_1,x_2,x_3,x_4,x_5,x_6]}^T={[{\rho }_1,t_1,W_1,{\rho }_2,t_2,W_2]}^T\end{array}$$

(9)

The material type is a discrete variable that must guarantee the device will not fail due to loading. This discards materials without a defined fatigue limit, such as polymerics, ceramics or compounds. There are studies establishing the endurance limits for metal alloys [42,47]. Alloys in

Table 1. Physical properties of the selected materials.

		Density $\mathit{\rho }$	Elastic Modulus E	Endurance Limit ${\mathit{S}}_{\mathit{e}}$
Type	Material	[kg/m${}^3$]	[Pa]	[Pa]
1	Aluminum 7075 T6	2810	$71.7\times {10}^9$	$121\times {10}^9$
2	Steel AISI 1045 CD	7870	$205\times {10}^9$	$236\times {10}^9$
3	Steel AISI 304L	8000	$195.4\times {10}^9$	$240\times {10}^9$

were selected because they met the design requirements, are used in medical non invasive applications, and are commercially available. They can be changed in type and/or number to meet specific design requirements.

3.2. Limits and Constraints

The desired values of the security factor and maximum deflection are $n_d=2$ and ${\delta }_d=5\times {10}^{-4}$ m, respectively, and the ranges of the design variables are defined in expressions (10)–(12). Dimensions $d=8\times {10}^{-3}$ m and $l=8.3\times {10}^{-2}$ m are previously set, but can be included as design variables to be calculated if required.

$$\begin{array}{c}\hfill x_1,x_4=\{2810,7870,8000\}\end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{d}{2}\le x_2,x_5\le 6d\end{array}$$

(11)

$$\begin{array}{c}\hfill 1.25d\le x_3,x_6\le 3.3d\end{array}$$

(12)

Equality expressions (13) guarantee a specific stiffness of the links, while inequality expressions (14) guarantee a minimum security factor for each link. Constraints in inequality expressions (15) are derived from the condition $W>t$ in Equation (5):

$$\begin{array}{c}\hfill h_1={\delta }_{max1}-{\delta }_d=0,h_2={\delta }_{max2}-{\delta }_d=0\end{array}$$

(13)

$$\begin{array}{c}\hfill g_1=n_d-\frac{S{e}_1}{{\sigma }_1^{\prime }}\le 0,g_2=n_d-\frac{S{e}_2}{{\sigma }_2^{\prime }}\le 0\end{array}$$

(14)

$$\begin{array}{c}\hfill g_3=(1.1\ast t_1)-W_1\le 0,g_4=(1.1\ast t_2)-W_2\le 0\end{array}$$

(15)

3.3. Problem Complexity

The relation between the feasible zone and the total search space of an optimization problem can be used to express its complexity [48]. This parameter, represented by $C_{\rho }$ and calculated by Equation (16), corresponds to the percentage of feasible vectors F in an arbitrary large number of random solutions S. As can be seen, a very low value of $C_{\rho }$ corresponds to an extremely high complexity, that requires a higher computing effort. One million randomly generated solutions [49] were used to calculate this complexity for the selected problem, and no feasible vectors were found between them, producing $C_{\rho }=0$. The computational characteristics of the design problem are shown in

Table 2. Computational characteristics of the case study.

Design Variables	Equality Constraints	Inequality Constraints	Complexity ${\mathit{C}}_{\mathit{\rho }}$
6	2	4	0

:

$$C_{\rho }=\frac{F}{S}\times 100\%$$

(16)

4. Optimization Method

The first options to consider when solving numerical optimization problems are mathematical programming techniques. Most of those methods require both the objective function and the constraints to be doubly and continuously differentiable with respect to each decision variable. Nevertheless, in many problems, the accomplishment of these conditions cannot be guaranteed, and alternative solution methods are required. Metaheuristics are approximate stochastic methods that have proved to be quite efficient for solving complex engineering design problems, such as the presented in this work.

4.1. Improved Harmony Search

Diverse versions of HS have been developed to improve its stability and the quality of the generated solutions. ImHS [50] is a variant that includes two modifications to enhance the performance of HS without affecting its simplicity. The first one is the change of the parameter bandwidth $bw$, reducing it at each algorithm iteration to improve the exploitation process. Further, $bw$ is personalized for each design variable considering the individual value ranges. The second modification is the inclusion of an additional parameter, the rate of improvement adjustment $r_{ia}$, to give more variety to the process of combining different options to create new solutions, improving the exploration.

4.2. Computational Implementation

Most metaheuristics were developed to solve unconstrained problems. Different schemes have been implemented for constraint handling (a complete description can be found in [51]). In this sense, the feasibility rules of Deb [52] are very popular because of its simplicity and high efficiency. In this work, they are combined with ImHS by applying a Constraint Violation Sum ($CVS$) for comparison-then-selection of solutions. $CVS$ is an accumulator for the absolute excess value of each unmet constraint, to directly measure the feasibility grade of a solution. The rules in ImHS work as follows:

• Between two infeasible harmonies, pick the one with the lowest $CVS$.
• Between a feasible and an infeasible harmonies, select the feasible one.
• Between two feasible harmonies, take the one with the best objective-function value.

ImHS was also modified for processing mixed variables, with a round operation [53] that transforms a real number into the nearest integer to generate the index corresponding to a specific material. Algorithm 1 corresponds to the pseudo-code of the implemented ImHS. It can easily be tuned by trial and error, in this case with $HM=16$ harmonies, $Iterations=5000$ iterations, $R_{accept}=0.85$, and $R_{pa}=0.45$, $r_{ia}=0.66$. However, diverse strategies for a better tuning can be applied, as explained in [54,55].

Algorithm 1:Improved Harmony Search (ImHS)

5. Results

ImHS had a 100% feasible rate in the solution of the case study, since each of their simulations produced one or more feasible solutions. In

Table 3. First five best solutions of ImHS.

		Link 1				Link 2
Rank	Type	Thickness (m)	Width (m)		Type	Thickness (m)	Width (m)	Total Weight (N)
1	1	0.00725	0.02409		1	0.00672	0.01772	0.77042
2	1	0.00720	0.02444		1	0.00780	0.01585	0.77990
3	1	0.00615	0.02382		1	0.00686	0.02149	0.78011
4	1	0.00866	0.01990		1	0.01021	0.01383	0.78143
5	1	0.00951	0.01806		1	0.00563	0.02514	0.81949

, the first five best solutions are presented, while

Table 4. Statistics of the performance of ImHS for the case study.

Best	Worst	Average	Median	Variance	Std Dev
0.77042	1.07628	0.86923	0.87400	0.00745	0.08631

includes a statistical analysis of the performance. There is a similarity between the average value and the median, indicating a homogeneous distribution of the results. In addition, the dispersion of the result set in relation to the arithmetic mean, given by the variance, shows a steady operation of the algorithm. For comparison purposes, the case study was also solved with Differential Evolution (DE) [56], a popular and efficient metaheuristic for solving constrained numerical optimization problems. DE rand/1/bin was selected for this case because of its simple implementation.

Table 5. Best solutions of ImHS and ED/rand/1/bin.

		Link 1				Link 2
Algorithm	Type	Thickness (m)	Width (m)		Type	Thickness (m)	Width (m)	Total Weight (N)
ImHS	1	0.00725	0.02409		1	0.00672	0.01772	0.77042
DE	1	0.01097	0.01593		1	0.01099	0.01590	0.86021

shows the best solution obtained by each algorithm. As can be seen, the solution from ImHS surpassed the best value of DE requiring a considerable lower amount of resources, 30,000 vs. 500,000 evaluations of the objective function.

The best result of ImHS was tested with the Finite Element Method (FEM) using ANSYS Workbench 2015, to verify if it met the problem constraints and boundaries. The simulation includes all the operation loads and stress sources considered in the design (Figure 3). The red end of each link shows the maximum deflection. These deflections were in the order of $1.18\times {10}^{-4}$ m, lower than the established maximum limit of $5\times {10}^{-4}$ m, validating that the links will operate within the operational ranges without failure.

6. Discussion

In this paper, the structural design of a wrist rehabilitation system is carried out by a concurrent methodology, with a novel approach where the optimization process works simultaneously on both the material selection and the constraint meeting. A metaheuristic is applied to the dimensional synthesis and the material selection of more than one structural element at a time, considering the interaction and joint effect while designing. The modeling of the system as an optimization problem adds flexibility to its solution. Once it is modeled, the problem can be adapted to configuration changes. There is no limitation on the materials to be chosen, and the geometries to be analyzed can be as complex as is required by the design, demonstrating the versatility of the method.

The associated optimization problem is a complex one solved with the Improved Harmony Search algorithm, modified for handling constraints and mixed variables. It generates design alternatives that can be evaluated from different perspectives such as cost, ease of manufacture and assembly. The results show that the performance of ImHS is stable and with a low computational cost, making it a good option to solve real-life problems if they are modeled as constrained cases. The quality of the design results guarantees that the physical implementation of the rehabilitation system will be safe for the patients, since they met the safety constraints established in the problem modeling.

As future work, we consider the manufacturing of a prototype of the wrist rehabilitation device designed with the proposed method, in order to test it with real users. For the prototype, it is necessary to develop a control system, but the considerations for the concurrent design will simplify that task. Finally, we suggest considering to implement a strategy for dynamic parameter tuning to improve the performance of ImHS.