1 Introduction

The inverse kinematics problem for the robots involves determining the joints values that match the input information of the end effector’s position and direction (Köker and Çakar 2016). These matched values will ensure that subsequent robot control process will follow the desired trajectory. This is one of very important issues in the robotic field because it is related to other aspects such as motion planning, dynamic analysis and control (Huang et al. 2012). Traditionally, there are some solutions to solve the inverse kinematics for robots such as: geometry method (Shanda et al. 2019; Fan et al. 2019) is the method using geometric and trigonometric relationships to resolve, the iteration method is the often required inversion of a Jacobian matrix. However, when applying them to solve the inverse kinetics of robots, especially redundant driving robots, it is often much more complicated and time-consuming. The reason is the nonlinearity of the formulas and the geometry transformation complexity between the working space (task space) and the joint space. In addition, the difficult point is in the singularity, the multiplicity of these formulas as well as the necessary variation of the formulas corresponding to the changes of different robot structures (Rubio et al. 2014; Kou et al. 2014; Lopez-Franco et al. 2018).

For upper exoskeleton rehabilitation robot type which normally has the same redundancy as the human arm, configuring the joint angles of the robot in accordance with those of human’s arm is one of the crucial control mechanisms to minimize the energy exchange between people and robot. In Kim and Rosen (2015), the authors proposed the redundancy resolution algorithm uses a wrist position and orientation algorithm that allows the robot to form the natural human arm configuration as the operator changes the position and orientation of the end effector. In this way, the redundancy of the arm is mathematically expressed by the swivel angle. In Wang et al. (2019), the resolution of the redundancy through computation of the swivel angle made it possible to obtain unique solutions for the joint space. However, in this way knowing or predicting the direction of the endpoint is a prerequisite. Authors in Li et al. (2019) used the reverse coordinate method to complete the inverse kinematics solutions and also proposed a new multicubic polynomial-interpolation method for planning joint trajectories.

In addition to those existing methods of solving the IK problems, in recent years, the application of meta-heuristic optimization algorithm in solving inverse kinematics has become increasingly common. The optimization-based approach has represented suitable ways to overcome the aforementioned problems. Indeed, these methods are generally based on the ways to choose objective function, constraints, and the chosen solution approach as well. There are many forms for objective functions based on some criteria such as singularity avoidance, obstacle avoidance, and physical constraints (Hwang and Jeon 2015; Whitney 2010; Runarsson and Yao 2005). In Kenwright (2017), the author used a Differential Evolutionary algorithm to generate approximated solutions for a range of IK problems such as positions, orientation, physical factors like overall center-of-mass location. In Karpinska et al. (2012), the redundant IK is formulated as the minimization of an infinity norm and translated to a dual optimization problem. In Laitenberger et al. (2014), a method based on optimization was proposed to implicitly resolve inverse kinematics of human upper limb. Authors (Dereli and Köker 2020) compared the position deviation and calculation time of two meta-heuristic techniques of PSO and ABC. The results showed that both algorithms successfully find the angles of the joints of a redundant manipulator robot with 7-DoF. Urrea and Saa (2020) presented a graphics simulator that allows for characterizing the kinematics and dynamics behavior of redundant planar manipulator robot. They calculated the IK based on the probabilistic method called Simulated Annealing.

In this study, we proposed to improve our algorithm (Bui et al. 2013) by improving the searching area in joint space. Previously, we tested this improvement into the DE algorithm and have proven effectiveness (Thanh-Trung et al. 2020). Therefore, in this study we proposed to use algorithm of self-adaptive control parameters in Differential Evolution with search space improvement (Pro-ISADE) to resolve the inverse kinematics for the redundant driven model. And the research model towards applying is the human arm model in some ADL movements and sports activities. This is a very important research task for the later stage of designing and controlling the rehabilitation robot. In addition, it can also be a prerequisite to better understand the natural and optimal movement of the human arm’s joints. From which we can apply to more widely later in industrial robots, serve robots, etc. After experimenting and simulating, we found the advantage of the proposed algorithm is to significantly reduce the complexity and volume of calculations compared with the other methods. The results of the predicted joints' angles from the application of the algorithm were compared with those measured by the Exoskeleton type Human Motion Capture System (E-HMCS) showing the similarity.

The remainder of the paper is divided into the following sections: Sect. 2 describes the model of the human arm. Next, the theory of ISADE algorithm and ISADE algorithm with improvement search space (Pr-ISADE) will be presented in Sect. 3. Section 4 covers the setup process to make experiment and run the proposed algorithm. The predicted results after applying the algorithm and comparing with the measured results are shown in Sect. 5. Finally, the conclusions are outlined in Sect. 6.

2 Human arm and E-HMCS model

As described above, in order to determine the endpoint trajectories of the human arm in some daily activities, research needs to have equipment to measure these movements. With the difficult access to modern and expensive kinetic measurement systems such as Advanced Realtime Tracking (ART) (https://ar-tracking.com/en), at the same time, for the purpose of evaluating the proposed algorithm, this study has suggested and manufactured an Exoskeleton type Human Motion capture system (E-HMCS) (Thanh-Trung et al. 2021). In addition, to convert between the angular values measured from the device to the corresponding angular joints value of the human arm, the study also proposed a 3D model of the human arm (Trung et al. 2020). This section will present a summary of the design of these models as well as the conversion of corresponding values between measured values and human hand angles.

2.1 Human arm and E-HMCS model

When ignoring the movement of scapular and clavicle motions, the Human arm is considered composed of links connected to each other by three anatomical joints (shoulder joint, elbow joint and wrist joint) (Thanh-Trung et al. 2020). With these joints, in this study, 7 DoFs were analyzed including shoulder flexion/ extension (Shld. Flex./Ext.), shoulder abduction/ adduction (Shld. Abd./Add.), shoulder internal/External rotation (Shld. Int./Ext. Rot.), elbow flexion/ extension (Elb. Flex./Ext.), wrist pronation/ supination (Wrist Pro./Sup.), wrist flexion/ extension (Wrist Flex./Ext.) and wrist Radial/Unal deviation (Wrist Rad./Unal. Dev.). In order to support the research and manufacturing process of upper limb rehabilitation robot, the authors also proposed a human arm 3D model that ensures the basic anthropometric parameters of Asians (Trung et al. 2020). The human arm anthropometric parameters and Range of Motion (RoM) are shown in the Table 1 and the reason will be explained in the calculation below Table 2, respectively.

Table 1 Some anthropometric parameters of human arm and human arm model
Full size table
Table 2 RoM of human arm and human arm model
Full size table

The structure model of the joints is shown in detail as Fig. 1. The shoulder joint is composed of 3 axes that are perpendicular to and intersect as Fig. 1a–c. The remaining 4 joints including (Elb. Flex./Ext.), (Wrist Pro./Sup.), (Wrist Flex./Ext.) And (Wrist Rad./Unal. Dev.) are represented in the Fig. 1d–f. The typical movements of human arm and corresponding of proposed model are presented in Fig. 2. The segments' length of model is chosen based on the segments' length of the person wearing the E-HMCS.

Fig. 1
figure 1

Joints’ structure of human arm model (Thanh-Trung et al. 2021)

Full size image
Fig. 2
figure 2

Movements of human arm and corresponding proposed human arm model (Thanh-Trung et al. 2021)

Full size image

For the purpose of evaluating the efficiency of the proposed algorithm, and determining the endpoint trajectory in ADL activities and sport movements, the research also developed an Exoskeleton type Human Motion Capture System (E-HMCS) (Thanh-Trung et al. 2021). The device uses 5 Potentiometers (POT) of CPP-35B (Midori Precision Co., Ltd.) and 2 Encoders of E6B2-CWZ6C (Omron Corp.), those are placed at joints as in Fig. 3a. The electric setups were presented as Fig. 3b. This is a measuring device, which has structure like an exoskeleton type, was fabricated by 3D printing technology and has properties of cheap, simple calculation process. With these features, this device is also easily accessible to developing countries, limiting resources for expensive systems such as Video motion capture system. The E-HMCS was produced and tested as in the Fig. 4. The overall specification of the system is summarized in the Table 3.

Fig. 3
figure 3

Setup sensors for E-HMCS device

Full size image
Fig. 4
figure 4

Exoskeleton type human motion capture system (E-HMCS)

Full size image
Table 3 Basic parameters of E-HMCS
Full size table

The human arm model wearing the E-HMCS device is shown in Fig. 5. The goal of this study is to use this measuring device to obtain the endpoint trajectory as well as the measured joints’ angle values corresponding. The endpoint trajectories are then put into proposed algorithm (Pro-ISADE algorithm) to predict the joints angle values and compare with the above measured values. With the aim of evaluating the coordinates of the end point position, the model will measure as well as predict for 5 measured angle values including \({q}_{1},{q}_{2},{q}_{3},{q}_{4},{q}_{5}\). The reason will be explained in the calculation below.

Fig. 5
figure 5

Human arm model wearing the E-HMCS

Full size image

2.2 Forward kinematic of the E-HMCS system

Forward kinematics analysis uses the kinematic equations of a robot to compute the position and orientation of the end-effector from specified values for the joint parameters. The homogeneous transformation matrix can be used to obtain the forward kinematics of the robot manipulator, using the DH parameters in Eq. 1:

$$ T_{{i - 1}}^{i} = \left[ {\begin{array}{*{20}c} {C\theta _{i} } & { - S\theta _{i} } & 0 & {a_{i} } \\ {S\theta _{i} Ca_{i} } & {C\theta _{i} Ca_{i} } & { - Sa_{i} } & { - d_{i} Sa_{i} } \\ {S\theta _{i} Sa_{i} } & {S\theta _{i} Sa_{i} } & { - Ca_{i} } & { - d_{i} Ca_{i} } \\ 0 & 0 & 0 & 1 \\ \end{array} } \right] $$
(1)

where: S and C denote the sine and cosine functions (Table 4)

$$ T_{0}^{5} = {\text{~}}T_{0}^{1} \times T_{1}^{2} \times T_{2}^{3} \times T_{3}^{4} \times T_{4}^{5} $$
(2)
Table 4 D-H parameters
Full size table

The position of end-effector coordinates in the working space are determined by:

$$ \begin{aligned} xE = & \left( {2869 \times cos \left( {q1} \right) \times sin \left( {q2~ - ~pi/2} \right)} \right)/10~ - ~336 \times sin \left( {q4} \right) \\ & \times \left( {sin \left( {q1} \right) \times sin \left( {q3} \right)~ - ~cos \left( {q1} \right) \times cos \left( {q3} \right) \times cos \left( {q2~ - ~pi/2} \right)} \right)~ \\ & + ~336 \times cos \left( {q1} \right) \times cos \left( {q4} \right) \times sin \left( {q2~ - ~pi/2} \right) \\ \end{aligned} $$
$$ \begin{aligned} yE = & \left( {2869 \times sin \left( {q1} \right) \times sin \left( {q2~ - ~pi/2} \right)} \right)/10~ + ~336 \times sin \left( {q4} \right) \\ & \quad \times \left( {cos \left( {q1} \right) \times sin \left( {q3} \right)~ + ~cos \left( {q3} \right) \times cos \left( {q2~ - ~pi/2} \right) \times sin \left( {q1} \right)} \right)~ \\ & \quad + ~336 \times cos \left( {q4} \right) \times sin \left( {q1} \right) \times sin \left( {q2~ - ~pi/2} \right) \\ \end{aligned} $$
$$ \begin{aligned} zE = & \left( {2869 \times cos \left( {q2~ - ~pi/2} \right)} \right)/10~ + ~336 \times cos \left( {q4} \right) \times cos \left( {q2~ - ~pi/2} \right) \\ & \quad - ~336 \times cos \left( {q3} \right) \times sin \left( {q4} \right) \times sin \left( {q2~ - ~pi/2} \right) + 143 \\ \end{aligned} $$
(3)

We can see that position of end-effector coordinates depends on joint variables \({q}_{1},{q}_{2},{q}_{3},{q}_{4}\). Joint variables \({q}_{5},{q}_{6},{q}_{7}\) determine the orientaion of the end-effector. To prove the effectiveness of the proposed Pro-ISADE algorithm, the study added one more joint of \({q}_{5}\) in the process of dealing with the inverse kinetics problem.

2.3 Convert the measured angle to the human arm angles

Measuring model E-HMCS and human arm model have different order of the first two joints {\({q}_{1},{q}_{2}\)} as well as different origin coordinate system (Fig. 2). Therefore, they will have different D-H parameters and the matching angle 1, 2, 3 of the 2 models will be different. From angular \({q}_{4}\) onwards, the angle of the measured model will coincide with the angle of the arm model. To convert angles from gauge to human arm, it is needed to find \(\left\{{q}_{1}^{\prime},{q}_{2}^{\prime},{q}_{3}^{\prime}\right\}\) angle of human arm so that the direction of link 2 (elbow point) of 2 models are the same. After having the results of the joint angle and the measured position, we will calculate the matrix for determining the elbow direction of the measurement model with the form:

$$ A_{3} = \left( {\begin{array}{*{20}c} {B_{{11}} } & {B_{{12}} } & {B_{{13}} } \\ {B_{{21}} } & {B_{{22}} } & {B_{{23}} } \\ {B_{{31}} } & {B_{{32}} } & {B_{{33}} } \\ \end{array} } \right) $$
(4)

where

$$B11=cos(q1) \times cos(q3) \times cos(q2 - pi/2) - sin(q1) \times sin(q3)$$
$$B12= cos(q1) \times sin(q2 - pi/2)$$
$$B13= - cos(q3) \times sin(q1) -cos(q1) \times cos(q2-pi/2) \times sin(q3)$$
$$B21=cos(q1) \times sin(q3)+cos(q3) \times cos(q2-pi/2) \times sin(q1)$$
$$ B22 = - sin \left( {q1} \right) \times sin \left( {q2~ - \frac{{pi}}{2}} \right),\,cos \left( {q1} \right) \times cos \left( {q3} \right) - cos \left( {q2~ - \frac{{pi}}{2}} \right) \times sin \left( {q1} \right) \times sin \left( {q3} \right) $$
(5)
$$B23=cos(q1) \times cos(q3)-cos(q2-pi/2) \times sin(q1) \times sin(q3)$$
$$B31=-cos(q3) \times sin(q2-pi/2)$$
$$B32=-cos(q2-pi/2)$$
$$B33=sin(q3) \times sin(q2-pi/2)$$

Matrix determines the elbow direction of the human arm model (compared with the original coordination of the measuring model) has the form:

$$ C_{3} = \left( {\begin{array}{*{20}c} {C_{{11}} } & {C_{{12}} } & {C_{{13}} } \\ {C_{{21}} } & {C_{{22}} } & {C_{{23}} } \\ {C_{{31}} } & {C_{{32}} } & {C_{{33}} } \\ \end{array} } \right) $$
(6)

where

$$C11=cos(q^{\prime}1) \times cos(q^{\prime}2) \times cos(q^{\prime}3)-sin(q^{\prime}1) \times sin(q^{\prime}3)$$
$$C12=cos(q^{\prime}1) \times sin(q^{\prime}2)$$
$$C13=-cos(q^{\prime}3) \times sin(q^{\prime}1)-cos(q^{\prime}1) \times cos(q^{\prime}2) \times sin(q^{\prime}3)$$
(7)
$$C21=-cos(q^{\prime}3) \times sin(q^{\prime}2)$$
$$C22=cos(q^{\prime}2)$$
$$C23=sin(q^{\prime}2) \times sin(q^{\prime}3)$$
$$C31=cos(q^{\prime}1) \times sin(q^{\prime}3)+cos(q^{\prime}2) \times cos(q^{\prime}3) \times sin(q^{\prime}1)$$
$$C32=sin(q^{\prime}1) \times sin(q^{\prime}2)$$
$$C33=cos(q^{\prime}1) \times cos(q^{\prime}3)-cos(q^{\prime}2) \times sin(q^{\prime}1) \times sin(q^{\prime}3)$$

With condition: B3 = C3 we can find out the values of \(\{{q}_{1}^{\prime},{q}_{2}^{\prime},{q}_{3}^{\prime}\}\) for human arm as following:

$$ \begin{aligned} q^{\prime}1{\text{~}} = & - 2 \times atan((2 \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right){\text{~}} - {\text{~}}2 \times tan\left( {q1/2} \right)^{2} \\ & \quad \times tan(q2/2 - pi + {\text{~}}((tan\left( {q1/2} \right)^{2} + tan(q2/2 - {\text{~}}pi/4)^{2} {\text{~}} \\ & \quad + tan\left( {q1/2} \right)^{2} \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right)^{2} - {\text{~}}4 \times tan\left( {q1/2} \right) \\ & \quad \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right){\text{~}} + {\text{~}}1) \times (tan\left( {q1/2} \right)^{2} + tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right)^{2} \\ & \quad + tan\left( {q1/2} \right)^{2} \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right)^{2} + {\text{~}}4 \times tan\left( {q1/2} \right) \\ & \quad \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right) + {\text{~}}1))^{{1/2}} )/(tan\left( {q1/2} \right)^{2} {\text{~}} - {\text{~}}tan(q2/2 - {\text{~}}pi/4)^{2} \\ & \quad - tan\left( {q1/2} \right)^{2} \times tan\left( {q2/2{\text{~}} - {\text{~}}pi/4} \right)^{2} + {\text{~}}1)){\text{~}} \\ \end{aligned} $$
$$ \begin{array}{lllll} {q^\prime }2 = & - 2 \times atan\left( {\left( {\left( {cos\left( {q1} \right) + 1} \right) \times \left( {sin\left( {q2} \right) + 1} \right)} \right.} \right.\\ & \quad \times {\left( { - \left( {cos{{\left( {q2} \right)}^2} \times sin{{\left( {q1} \right)}^2} - 1} \right)/\left( {cos{{\left( {q1/2} \right)}^4} \times {{\left( {sin\left( {q2} \right) + \;1} \right)}^2}} \right)} \right)^{1/2}}\\ & \quad /\left. {\left( {2 \times \left( {cos\left( {q2} \right) \times sin\left( {q1} \right) + \;1} \right)} \right)} \right) \end{array} $$
$$ \begin{aligned} {q^\prime }3 = & - 2 \times atan\left( {\left( {{{\left( {tan\left( {q2/2} \right)\; + \;1} \right)}^2} \times (2 \times tan\left( {q3/2} \right)} \right.} \right.\\ & \quad + {\left( { - \left( {cos{{\left( {q2} \right)}^2} \times sin{{\left( {q1} \right)}^2} - \;1} \right)/(cos{{\left( {q1/2} \right)}^4} \times {{(sin\left( {q2} \right) + 1)}^2})} \right)^{(1/2)}}\\ & \quad + tan{\left( {q3/2} \right)^2} \times {\left( { - \left( {cos{{\left( {q2} \right)}^2} \times sin{{\left( {q1} \right)}^2} - \;1} \right)/\left( {cos{{\left( {q1/2} \right)}^4} \times {{\left( {sin\left( {q2} \right)\; + \;1} \right)}^2}} \right)} \right)^{(1/2)}}\\ & \quad - \;2 \times tan{\left( {q1/2} \right)^2} \times tan\left( {q3/2} \right) + sin\left( {q2} \right) \times \left( { - \left( {cos{{\left( {q2} \right)}^2} \times sin{{\left( {q1} \right)}^2} - \;1} \right)/} \right.\left( {cos{{\left( {q1/2} \right)}^4}} \right.\\ & \quad \times {(sin\left( {q2} \right) + \;1)^2}){)^{(1/2)}} + 2 \times tan\left( {q1/2} \right) \times sin\left( {q2} \right) + tan{\left( {q3/2} \right)^2}\\ & \quad \times sin\left( {q2} \right) \times {\left( { - \left( {cos{{\left( {q2} \right)}^2} \times sin{{\left( {q1} \right)}^2} - \;1} \right)/\left( {cos{{\left( {q1/2} \right)}^4} \times {{\left( {sin\left( {q2} \right) + 1} \right)}^2}} \right)} \right)^{(1/2)}}\\ & \quad - 2 \times tan\left( {q1/2} \right) \times tan{\left( {q3/2} \right)^2} \times sin\left( {q2} \right)))/\left( {\left( {sin\left( {q2} \right)\; + \;1} \right)} \right. \times \left( {tan{{\left( {q1/2} \right)}^2}} \right.\\ & \quad - tan{\left( {q2/2} \right)^2} + tan{\left( {q3/2} \right)^2} + tan{\left( {q1/2} \right)^2} \times tan{\left( {q2/2} \right)^2} - tan{\left( {q1/2} \right)^2}\\ & \quad \times tan{\left( {q3/2} \right)^2} + tan{\left( {q2/2} \right)^2} \times tan{\left( {q3/2} \right)^2} + \;8 \times tan\left( {q1/2} \right) \times tan\left( {q2/2} \right)\\ & \quad \left. { \times tan\left( {q3/2} \right)\; - \;tan{{\left( {q1/2} \right)}^2} \times tan{{\left( {q2/2} \right)}^2} \times tan{{\left( {q3/2} \right)}^2} - 1} \right) \end{aligned} $$
(8)

3 Proposed optimization algorithm

In the Bui et al. (2013) we suggested to develop a new version of DE algorithm that can automatically adapt the learning strategies and the parameters settings during evolution. The main ideas of the ISADE algorithm are summarized below.

3.1 Mutation operator

ISADE probabilistically selects one out of several available learning strategies in the mutation operator for each individual in the current population. In this research, we select three learning strategies in the mutation operator as candidates: “DE/best/1/bin”, “DE/best/2/bin” and “DE/rand to best/1/bin” that are respectively expressed as (Fig. 6):

$$ DE/best/1{\text{:}}\,\,\,V_{{i,j}}^{G} = X_{{best,j}}^{G} + F \times (X_{{r1,j}}^{G} - X_{{r2,j}}^{G} ) $$
(9)
$$ DE/best/2{\text{:}}\,\,V_{{i,j}}^{G} = X_{{best,j}}^{G} + F \times (X_{{r1,j}}^{G} - X_{{r2,j}}^{G} ) + F \times (X_{{r3,j}}^{G} - X_{{r4,j}}^{G} ) $$
(10)
$$ DE/rand\,to\,best/1{\text{:}}\,\,V_{{i,j}}^{G} = X_{{r1,j}}^{G} + F \times (X_{{best,j}}^{G} - X_{{r1,j}}^{G} ) + F \times (X_{{r2,j}}^{G} - X_{{r3,j}}^{G} ) $$
(11)

where \(i = \left\{ {1,2, \ldots ,NP} \right\};j = \left\{ {1, \ldots ,D} \right\}\) are current population and design variable, respectively.

Fig. 6
figure 6

ISADE flowchart

Full size image

\( DE/Rand{\text{~}}to{\text{~}}best/1/bin\) strategy usually demonstrates good diversity while the \( DE/{\text{~}}best/1/bin\) and \( DE/best/2/bin\) strategy show good convergence property, which we also observe in our trial experiments.

3.2 Adaptive scaling factor F and crossover control parameter CR

In the ISADE algorithm, the author suggested to use the sigmoid function to control neighborhood parameter. we sort the particles by estimating their fitness. A ranked particle is labeled with ranked number and assigned F that corresponds with its number. The formula for F by sigmoid function as following:

$$ F_{i} = \frac{1}{{1 + \exp \left(\alpha \times \frac{{i - \frac{{NP}}{2}}}{{NP}}\right)}} $$
(12)

where \(\alpha ,i{\text{~}}\) denote the gain of the sigmoid function, particle of the \({\text{~}}i_{{th}} {\text{~}}\) in \(NP\), respectively.

For better performance of ISADE, the scale factor F should be high in the beginning to have much exploration and after curtain generation F needs to be small for proper exploitation. Thus, we proposed to calculate the F as follow:

$$ F_{{iter}}^{{mean}} = F_{{\min }} + (F_{{\max }} - F_{{\min }} )\left( {\frac{{iter_{{\max }} - iter}}{{iter_{{\max }} }}} \right)^{{n{}_{{iter}}}} $$
(13)

where \(F_{{max}} ,F_{{min}} ,iter,iter_{{max}} {\text{~}}and{\text{~}}n_{{iter}}\) are the lower boundary condition of \(F\), upper boundary condition of \(F\), current generation, maximum generation and nonlinear modulation index, respectively.

The author introduced a novel approach of scale factor \({\text{~}}F_{i}\) of each particle with their fitness in Eq.13. Thus, in one generation the value of \(F_{i}^{{iter}} \left( {i = 1, \ldots ,NP} \right)\) are not the same for all particles in the population rather they are changed in each generation. The final value of scale factor for each generation is calculated as follow:

$$ F_{{iter}}^{i} = \frac{{F_{i} - F_{{iter}}^{{mean}} }}{2} $$
(14)

where \(iter = 1, \ldots ,{\text{~}}iter_{{max}} {\text{~}}and{\text{~}}i = 1, \ldots ,NP.\)

The control parameter \(CR{\text{~}}\) is adapted as following:

$$ CR_i^{G + 1} = {\rm{ }}\left\{ {\begin{array}{*{20}{l}} {ran{d_2}}&{if{\mkern 1mu} {\mkern 1mu} ran{d_1} \le \tau }\\ {CR_i^G}&{othewise} \end{array}} \right.{\mkern 1mu} $$
(15)

where rand1 and rand2 are uniform random values ∈ [0, 1], τ represents probabilities to adjust CR, same as(6) we assign τ = 0.10.

3.3 Proposed ISADE algorithms with searching space improvement (Pro-ISADE)

As mention in the introduction part, the disadvantage of many studies using optimization algorithms to solve the IK problem of redundant robots is almost only to focus on the results related to the optimal running process such as execution time, number of generation, etc., but have not yet considered the feasibility of the joints’ variable values. In order to overcome these drawbacks, the author of this research (Thanh-Trung et al. 2020) proposed this algorithm that is explained as following. The solution to improve the continuity of joints’ values is limiting the initialization domain of X. This help the program to achieve the dual goal of increasing calculation speed, accuracy and ensuring continuity for the value of joints’ variables. In this algorithm, firstly the robot from any position moves to the first point of the trajectory. With this first point, the initialization values for the particles are randomly selected in the full Range of Motion \(\left( {ROM} \right)\) of joints. In addition, the target function in this case has the form (Thanh-Trung et al. 2020):

$$ Func.1 = a \times \sqrt {\sum\limits_{k = 1}^5 {{{(q_1^k - q_0^k)}^2}} } + b \times \sqrt {\begin{array}{*{20}{l}} {{{({x_i} - {x_{ei}})}^2} + {{({y_i} - {y_{ei}})}^2} + {{({z_i} - {z_{ei}})}^2}}\\ { + {{(R{x_i} - R{x_e})}^2} + {{(Rx - R{x_e})}^2} + {{(R{y_i} - R{y_e})}^2} + {{(R{z_i} - R{z_e})}^2}} \end{array}} $$
(16)

where \(i~\left( {1/n} \right)\) are the points in the trajectories; \(k~\left( {1/5} \right)\) are the joints of the robot.

Values \(q_{1}^{k} {\text{~}}and{\text{~}}q_{0}^{k}\) (i = 1) are the joints’ variable values at the original position and 1st point on the trajectory, respectively; \((x_{i}\),\(~y_{i} ,z_{i} )\) and \((x_{{ei}}\),\(~y_{{ei}} ,z_{{ei}} )\) are the End-effector coordinates for the i-point found by the algorithm and the desired End-effector coordinates; \(\left( {R_{{xi}} ,R_{{yi}} ,R_{{zi}} } \right)\) and \(\left( {R_{{xei}} ,R_{{yei}} ,R_{{zei}} } \right)\) are corresponding rotation cosine angles performing orientation of the end-effector which are found by Algorithm and orientation of the desired end-effector; \(a,~b\) are penalty coefficients. Cost function as Eq. (16) ensures the energy spent in the joints to reach the 1st desired position is minimized. Besides, it also minimizes the distance error between the actual and desired end-effector position. The condition to stop algorithm is that the Cost Func.1 value is less than value of e or the number of iterations reaches \(iter_{{max}} ~.\)

In the case of solving the IK for consecutive point, after calculating for 1st point of the trajectory, the remaining points are calculated with a search limitation around the previous optimal joints’ values. By limiting this, the program's search space will be limited while ensuring the continuity of the joint variables. In this case, the target function will still be same as the function of 1st point, but it has coefficient a = 0.

Applied in this study, the research has initialized the first values of the joints (Zero positions of the joints) were equal to the measured values for the first point on the trajectory. This means \(q_{0}^{k} = q_{1}^{k} \,(k = 1 / 5)\).

4 Setup for experiments and proposed algorithm

4.1 Setup for measuring experiments of ADL and sport activities

In exercise activities, the research did experiments with two ADL activities of:

  • Drinking water task.

  • Brushing teeth task.

    In addition, the study also tested 2 sport exercises:

  • Playing ball task.

  • Catching ball task.

As in the Fig. 7, in this study, the measuring device E-HMCS was worn on a person with height 1.72 (m) and 22 years old. This person performed four groups of the complex ADL activities as aforementioned. The Figure presents two common ADL tasks of drinking water and brushing teeth. The drinking water task is divided into two stages of reaching and lifting. In order to perform these tasks, the experimental system should note a few points:

  • Adjust the segments of the measuring system to match the length of the wearer's segments.

  • The measuring device is firmly linked to the human hand by the strap band at upper arm.

  • Testers need to wear and become familiar with the equipment. It is important to note that when performing the activities, the examiner needs to do it naturally and comfortably. If this person feels hindered by the measuring equipment, the adjustment should be made to get the most natural feeling in the joint position.

Fig. 7
figure 7

Experiment tasks

Full size image

The experimenter wore the device and perform the ADL and sport activities. Each activity was implemented 20 experiments and with confirmation of comfort level after each session. Figure 8 below are the endpoint orbital results for all the missions. In addition to these results, we also got the results of the matching joints’ variables values. They will be used to compare with the predicted results when using the Pro-ISADE algorithm to solve the inverse kinetics problem.

Fig. 8
figure 8

Experiment results of endpoint trajectories of different tasks

Full size image

4.2 Setup for Proposed Pro-ISADE algorithm

When solving the IK problem for the 7-DOF serial robot manipulator of human arm, the study focused on three main aspects. The first of these is the sensitivity of the solution, or in the other word, the amount distance error of end effector. The second criterion was the execution time. In order to avoid the endless loop, the maximum numbers of generation \((iter_{{max}} ){\text{~}}\) was set as 600. And the final aspect is the searching space of joints’ variables. Normally, almost studies have been using the Range of Motion (\(RoM\)) of joints for this searching space. Our algorithm (Thanh-Trung et al. 2020) proposed to use the searching space of current generation is around previous optimal joints’ values. In the Thanh-Trung et al. (2020) the \(ubs_{{i + 1}} {\text{~}}and{\text{~}}lbs_{{i + 1}}\) are the joints’ upper and lower boundary of the current generation. when setting up the maximum distance error by the fitness value setting for the end effector position, the study set the value of \(1e - 5{\text{~}}\left( {mm} \right)\) for Pro-ISADE algorithm (Table 5).

Table 5 Optimization parameters used for proposed algorithm
Full size table

Running the IK problem in order to predict the 5 joints of \(q_{1} ,q_{2} ,q_{3} ,q_{4} ,q_{5}\) were coded by Matlab version 2019a and run on the computer equipped with an Intel Core i5-4258U @2.4 GHz processor and 8 GB Ram memory.

5 Results and discussion

The endpoint trajectories obtained from measuring experiments shown in Fig. 8 were inputted in the proposed Pro-ISADE algorithm to handle the inverse kinetics problem and obtain the predicted joints’ values. The results of the of inverse kinetics for the first phase of the drinking water task, the second stage for the drinking task, the brushing task, the playing ball task and the catching ball task are shown as shown in Figs. 9, 10, 11, 12 and 13, respectively. The evaluation results for each task include the endpoint deviation and comparison between the predicted joints values by the algorithm and the measured values by the experiment.

Fig. 9
figure 9

Results of Stage 1 Drinking water task (Reaching stage)

Full size image
Fig. 10
figure 10

Results of Stage 2 Drinking water task (lifting stage)

Full size image
Fig. 11
figure 11

Results of brushing teeth task

Full size image
Fig. 12
figure 12

Results of playing ball task

Full size image
Fig. 13
figure 13

Results of catching ball task

Full size image

For the drinking water task, the measured endpoint movement trajectories, i.e. the wrist joint, of the two movements of reaching a cup and elevating the cup to mouth are shown in Fig. 8a and b, respectively. After applying the proposed Pro-ISADE the Fig. 9 show the results for reaching stage. Meanwhile, the results for the lifting cup stage are shown in the Fig. 10. In these figures, the Fig. 9a and 10a partly show the quality of the Pro ISADE algorithm when solving the IK problem for these orbits. This quality is demonstrated by the very close grip of the simulated trajectories compared to reference trajectories that generated from measured experiments. Specifically, it can be seen as in the Table 6 the average distance error between the end point orbits for stage 1 and stage 2 are \(1.9551e - 6~\left( {mm} \right)\) and \(2.3364e - 6~\left( {mm} \right)\), respectively. Another very important result of this study is the quality of the predicted joints' angle values, that are gotten from solving inverse kinematics by using proposed Pro-ISADE algorithm, shown in Figs. 9b, c and 10b, c. According to these results, the predicted values of joint variables are continuous and very close to the measured values. For this drinking motion, in both phases the direction of end effector point, i.e. the wrist joint, is virtually unchanged. This is shown through the constant angle \(q_{3} ,q_{5}\) corresponding to the upper arm internal/external rotation and forearm pronation/supination movements. Deviations between the predicted values and the measured values of the angles are shown in Fig. 10c. These deviations are relatively small and are all less than 1 degree.

Table 6 Summarized results
Full size table

To demonstrate the effectiveness of the Pro-ISADE algorithm, the study tested inverse kinetics for stage 1 of drinking water task by using ISADE algorithm. The results of the solutions of the predicted joint values were shown in Fig. 14. It is easy to see that the algorithm can still determine the values of these joints’value, but these values fluctuate dramatically. It is completely different from the measured joint values shown in Fig. 9b. At the same time, these values are not controllable in terms of control.

Fig. 14
figure 14

Results in stage 1 of drinking water task when using ISADE algorithm

Full size image

In the second ADL activity of brushing teeth, the experimenter lifted the toothbrush from one point below the body to the mouth and perform the repetitive action of bringing the toothbrush to the left and right. Similar to the drinking water task, distance error between simulation and measured reference trajectories, results of predicted and measured joints values are demonstrated in the Fig. 11. Results after applying the algorithm Pro-ISADE showed that the simulated end effector trajectory closely followed the reference endpoint orbit. The average deviation between these 2 orbits is \(1.844e - 6~\left( {mm} \right)\). In addition, the predicted joints variable values are also very similar to the measured joints variable values. In these values the variables \(q_{3}, q_{5}\) are constant while the last stage of angle \(q_{2}\) is sinusoidal. The reason for this phenomenon is because in the brushing task as described, the experimenters hardly used the shoulder Internal- External rotation and forearm Pronation- Supination movement. Meanwhile, in order for the arm to perform the repetitive movement left to right, the \(q_{2}\) joint will move forward and reverse again and again. Figure 11b and c show the quality of the Pro-ISADE algorithm application to predict the values of the joints’ values. Similar to the task of drinking water, the proposed algorithm has handled the problem of inverse kinetics very well to ensure that the predicted angle results closely follow the joints angle values of the human’s normal movements.

To increase the convincing of the proposed algorithm, in this study we also chose 2 activities other than ADL. These two practice tasks are Playing ball and Catching ball tasks. The process of implementing, collecting and processing data are completely similar to the ADL tasks mentioned above. The results for these two tasks in turn are shown in Figs. 12 and 13. It can be easily seen that: similar to the two ADL activities mentioned above, in these two sports operations the proposed Pro-ISADE algorithm still handles the inverse kinematics problem very well. This quality is shown through the simulation endpoint orbital deviation compared to the reference trajectory and especially the value of the predicted joints’ variables compared to the measured joints’ values from the experiment.

After using the proposed Pro-ISADE algorithm to calculate the inverse kinetics problem for the measuring equipment system, the aforementioned comparison process has confirmed the high feasibility of the method, the study then transferred the values of the measuring device into the value of the human arm joints. Applying formula Equa.8 we can get the human joints’ values for the different tasks shown in Fig. 15.

Fig. 15
figure 15

Results of Solving IK problem for human arm joints in all tasks

Full size image

6 Conclusion

In this study, proposed algorithm of self-adaptive control parameters in Differential Evolution with search space improvement (Pro-ISADE) was used to solve the Inverse Kinematic problems of human arm in four different activities including 2 ADL and 2 sport activities. One of very important results is the predicted values of joint variables of the human arm. To evaluate this proposed algorithm, firstly an Exoskeleton type Human Motion Capture System was used to measure the wrist (end effector point) trajectories as well as the joints ‘value corresponding for all the above activities. In order to ensure the integrity of the operation, the experimenter should be familiar with measuring the equipment and perform the activities as naturally as possible. The measured endpoint orbits were then put into the suggested Pro-ISADE algorithm to predict the joints variable values. The experiment as well as the results after running the algorithm showed that the predicted endpoints and especially the predicted joint angle values of the IK problem closely follow the measured end effector trajectories and the joint angle values. Deviations between the predicted angle and the measured angle in the results can be caused by unnaturalness when wearing a rigid measuring frame to obtain the results. The results of the measuring system’s joints were then converted to values of the human joints to complete the inverse kinetics problem for human upper limb. Through the results obtained in the study, it can be concluded that the pro-ISADE proposed algorithm can be used to solve problems of inverse kinetics for natural human movements.