Proposal of Novel Model for a 2 DOF Exoskeleton for Lower-Limb Rehabilitation

Nowadays, engineering is working side by side with medical sciences to design and create devices which could help to improve medical processes. Physiotherapy is one of the areas of medicine in which engineering is working. There, several devices aimed to enhance and assist therapies are being studied and developed. Mechanics and electronics engineering together with physiotherapy are developing exoskeletons, which are electromechanical devices attached to limbs which could help the user to move or correct the movement of the given limbs, providing automatic therapies with flexible and configurable programs to improve the autonomy and fit the needs of each patient. Exoskeletons can enhance the effectiveness of physiotherapy and reduce patient rehabilitation time. As a contribution, this paper proposes a dynamic model for two degrees of freedom (2 DOF) leg exoskeleton acting over the knee and ankle to treat people with partial disability in lower limbs. This model has the advantage that it can be adapted for any person using the variables of mass and height, converting it into a flexible alternative for calculating the exoskeleton dynamics very quickly and adapting them easily for a child’s or young adult’s body. In addition, this paper includes the linearization of the model and an analysis of its respective observability and controllability, as preliminary study for control strategies applications.


Introduction
Currently a new concept of a human-machine robot is being studied and developed; different from traditional robots, exoskeletons are intelligent robot systems used widely around the world in many areas.Exoskeleton applications involve using them to increase human strength, speed, and performance in areas such as the military [1]; in medicine, they are used to train humans on surgery devices [2], to help humans perform surgeries [3], to increase the effectiveness of gait rehabilitation therapies [4] and for general physical therapies [5,6].The entertainment industry employs them in video games as controllers, and manufacture industry applications include teaching movements to robot arms tele-operated by exoskeletons [7], heavy load lifting, and many others.In this paper, we propose a dynamic model to design a two degrees of freedom leg exoskeleton for the rehabilitation of knee and ankle injuries.Several methods for exoskeleton modeling have been employed, some of which independently model the human leg and the exoskeleton frame [8].Among the different proposals, there are, for instance, developments that use Denavit-Hartenberg parameters to model the exoskeleton leg and the use of PID (Proportional-Integral and Derivative) controllers to control the angular position of each joint [2,9].Others use parallel robots instead of rigid structures coupled to the body [10], while still other kinds of models, like the one presented in [11], use basic mechanical and dynamical data to create a simple representation of the exoskeleton aiming to control purposes or virtual models to implement a dynamic and cinematic model using software such as OpenSim [12].Most of these proposed models for exoskeletons use Euler-Lagrange equations to determine the system dynamics, as this allows them to design robust controllers with gravity compensation and virtual joint torque control [13,14].However, in most of these cases some variables are assumed, which are not described by their respective studies.Other works consider the mathematical modelling of elements to be added to existent exoskeletons, for instance, the development of a hand crank generator for powering lower-limb exoskeletons [15].In this paper, a model for a two degrees of freedom (2 DOF) exoskeleton is developed using Euler-Lagrange mechanics where the exoskeleton structure and the leg are considered as one body, which is modeled as a two-link pendulum combined with a biomechanics model in order to include leg kinematics and dynamics as described in [16].This paper is organized as follows: Section 2 describes the model development and its linear approximation; Section 3 shows the simulation process as well as the simulation results of the nonlinear and linear model with a controllability and observability analysis; and finally, Section 4 shows the description of the exoskeleton considered.Conclusions are discussed with considerations for future research at the end of this paper.

Methods
This section shows the development of the proposed model to describe the dynamics of the exoskeleton attached to the leg acting over two joints (knee and ankle); it involves double pendulum dynamics combined with a popular biomechanical model.

Model Description
The proposed model for the exoskeleton is a variation of a two-link pendulum dynamic model, where biomechanical and anthropometric models are added to create a generalized, versatile, and flexible dynamic model which can describe the dynamics of a human leg performing movements in the sagittal plane.In addition, the model can be calculated and altered for use with a patient's individual characteristics.
For this purpose, a Hanavan biomechanical model was used, which describes human limbs as geometrical forms, specifically, truncated cones for legs and feet [17].To make sure that the model can be calculated for every kind of person, a generalized Hanavan model was calculated by a geometrical analysis of the truncated cone (see Figure 1) and its kinematics were combined with the anthropometric model of Drillis and Contini [18].To calculate the lower radius R l of the truncated cone, according to Figure 1, consider that: replacing it into the cone mass equation: Thus, the resulting expression for the lower truncated cone radius is: where R l and R h are the lower and higher radius of the truncated section, l the truncated cone height (all in meters), ρ is the density of the piece, and m its mass.As the lower radius is function of itself, the ratio R h /R l was calculated using direct measurements with a sample of 50 Colombians aged between 18 and 45 years.The length of each section of the leg was calculated using the Drillis and Contini model illustrated in Figure 2. From Equation (3), density (ρ) is calculated using the Drillis and Contini equation for average density (human body density), which is shown in Equation ( 4) [19].Average density is expressed in terms of the c parameter, which is called the ponderal index or index of corpulence; it is a measure of the contexture calculated as the ratio between mass and body  Average density is expressed in terms of the c parameter, which is called the ponderal index or index of corpulence; it is a measure of the contexture calculated as the ratio between mass and body Mass and distance to the center of mass are calculated using tabulated data shown in [19]; to simplify calculations, the center of mass distance of each limb segment was considered to be located at 50% of each cone length, resulting in the generalized Hanavan model of the human leg illustrated in Table 1.

Limb Section Parameters
Thigh where r t , r l , r f are the corresponding lower radii, R t , R l , R f are the higher radii, L is the limb length, and H T is the height of the person (all in meters); ρ is the density of the body in Kg/m 3 , m i is the segment mass, and M T is the person mass in [Kg].
Once the Hanavan model is generalized, it is necessary to calculate its kinematics.To do that, the moment of inertia across the "X" axes (see Figure 1) is calculated and is expressed as follows: This generalized biomechanical model is combined with the double link pendulum model to enhance it and describe the behavior of the exoskeleton coupled to the human leg.According to this, each link of the pendulum is considered as a truncated cone; the kinematics for the proposed model is described in Figure 3.The displacement of the center of mass is defined as: where , is the angular position of each link, , represents each link length.The angular position of the second link is referenced to the first link, so it will be zero if it is aligned with it.
The speed of center of mass is defined as: The moment of inertia across the center of mass of each segment is defined using Equation (5): The displacement of the center of mass is defined as: where θ 1,2 is the angular position of each link, L 1,2 represents each link length.The angular position of the second link is referenced to the first link, so it will be zero if it is aligned with it.
The speed of center of mass is defined as: The moment of inertia across the center of mass of each segment is defined using Equation ( 5): To find the dynamic model of the double pendulum enhanced with the biomechanical model, Lagrange equations of the second kind or Euler-Lagrange equations were used.Those are a reformulation of classical mechanics where the constraints are incorporated directly by the selection of generalized coordinates for the system being studied.Instead of forces, Lagrangian mechanics use the energies in the system, described using the Lagrangian; the nonrelativistic Lagrangian for a system of particles can be defined by: where T and V are the kinematic energy and potential energy of the system, respectively.The Euler-Lagrange equation is obtained from Equation ( 10) and, establishing that, if not all the forces acting on a system are derivable from a potential, those can be written as [20]: where L contains the potential of conservative forces (Equation ( 10)), q j represents the generalized coordinates of the system, Q j denotes the generalized forces acting outside the system, and j represents each particle or body analyzed.According to this, external forces acting over the exoskeleton coupled to the leg are considered; those external forces are torque and viscous friction at each joint, and the Euler-Lagrange equation for the system will be: where j represents links 1 and 2 of the pendulum, which denotes the leg and foot.The angle θ represents the position of each link, τ j is the applied torque, and b j is the viscous friction coefficient.
To calculate the Lagrangian, it is necessary to define the kinetic and potential energy of the system; kinetic energy is defined as follows [21]: where E CT is the translational kinetic energy and E CR is the rotational kinetic energy; v is the linear speed and ω is the angular speed of the center of mass, m represents bodies mass and I CM is the moment of inertia at the center of mass.Since there are two joints to be considered, the total kinetic energy T TS of the system will be expressed as: Potential energy is defined as: so total potential energy V TS for the system will be: where m 1,2 is the mass of each segment, g is the gravity used as 9.8 m/s 2 , and y 1,2 is the center of mass altitude calculated using the y coordinate of Equations ( 6) and (7).Using Matlab-MuPAD, dynamic equations which describe the exoskeleton are calculated by replacing Equations ( 8) and ( 9) into Equation ( 14) to calculate the total kinetic energy of the system; after calculating the kinetic and potential energy, the Lagrangian ( 10) is calculated and replaced into Equation (12).Once all the parameters are established, Equation ( 12) is solved and simplified for .. θ 1 and .. θ 2 , resulting in two nonlinear equations, one for each link describing the dynamic behavior of the proposed system.The resulting dynamic model is: (18) where: Those constants are calculated by using Table 1 parameters where R l is the leg higher radius, r l the leg lower radius, R f the foot higher radius, R f the foot lower radius, L 1 the leg length, and L 2 the foot length.
This model can be represented as the general form of the Euler-Lagrange equations, such that: where M is the inertial properties matrix, V the coriolis and centripetal vector, and G the gravity vector of the system.According to the generalized coordinates, the generalized equation for our system is: where the parameters of Equation ( 20) are: 2 ) 2 2 ) 2

Linear Approximation of the Model
To perform a comparison between the nonlinear model and a linear model suitable to design controllers, the system is linearized.
Since Equations ( 17) and ( 18) are second order, it is necessary to transform them into a pair of first-order differential equations, therefore the following variables are introduced: Those equations are linearized around an equilibrium point based on the Taylor series expansion [22].Thus, the Equations in (22) can be represented as: so, the linear approximation of the system will be calculated as follows: .
The equilibrium point was selected at the point where the system has lost all its energy; according to the reference settled in Figure 3 it is: Finally, the linear approximation of the system in state-space is: .
where P 2 and P 3 are the same as in Equations ( 17) and ( 18), so:

Simulation Results
The aim of the experiments is to validate the developed model by simulation.In this sense, the dynamic models (nonlinear and linear) are implemented to validate the results in the linearization process and the equivalence with the nonlinear model, after which the controllability and observability analysis are performed to the linearized model as a further step to verify the feasibility of control implementation.

Dynamic Model Simulation
This section describes the developed process to simulate the linear and nonlinear dynamics of the proposed model; for this purpose, Matlab-Simulink was used.Proposed linear and nonlinear models were simulated, calculating its parameters by considering a person with a height of 1.90 m and weight of 90 Kg.Viscous friction at each joint was assumed as b 1 = b 2 = 0.3 and applied torque as τ 1 = τ 2 = 0, to obtain a lightly damped response of the system.To simulate the exoskeleton, motors and structure weights were added to each link mass; the motors to be used are flat brushless DC motors coupled to a reduction gear capable of generating a maximum torque of 15 Nm.Those motors coupled with the reduction gear box weigh around 0.57 Kg, and the exoskeleton structure was assumed to weigh 0.13 Kg, amounting to a total of 0.7 Kg.These values were added to the second link mass; as the motor acting over the knee at the first link is not going to lift its own weight over the first link, only the structure weight was added.

Nonlinear Model Simulation
For simulating the nonlinear model, the ordinary differential equation (ODE) system composed by the generalized Euler-Lagrange equation was solved by using its explicit form: ..
where M is the inertial properties matrix, V the coriolis and centripetal vector, and G the gravity vector of the system.According to the generalized coordinates, the generalized equation for our system is: ..
where each matrix parameter is shown in Equation ( 21).After calculations with selected personal characteristics, they are defined as: Once the system of ordinary differential equations (ODEs) is defined, it is necessary to reduce it from two second-order equations to four first-order equations; this is done by applying the same process employed for linearizing the model.This time, for the nonlinear functions described by Equation ( 27), the following variables must be introduced: where f i (x, u) corresponds to each .. θ i expression given by solving Equation (27): where each matrix parameter is shown in Equation ( 21).After calculations with selected personal characteristics, they are defined as: ( , ) = 0.2706 cos( ) + 0.8133 Once the system of ordinary differential equations (ODEs) is defined, it is necessary to reduce it from two second-order equations to four first-order equations; this is done by applying the same process employed for linearizing the model.This time, for the nonlinear functions described by Equation ( 27), the following variables must be introduced: where ( , ) corresponds to each expression given by solving Equation ( 27): Simulink "ode15s" solver was selected to solve this ODE system, and as a nonlinear model reference uses horizontal axes as 0° degrees and the gravity vector (G(q)) is used, no initial conditions were needed.The simulation block diagram is illustrated in Figure 4 and the dynamic behavior of the model is shown in Figure 5.  Simulink "ode15s" solver was selected to solve this ODE system, and as a nonlinear model reference uses horizontal axes as 0 • degrees and the gravity vector (G(q)) is used, no initial conditions were needed.The simulation block diagram is illustrated in Figure 4 and the dynamic behavior of the model is shown in Figure 5.The model behavior was as expected; the first link falls and moves describing a damped sine wave, which loses energy due to friction causing it to stabilize in −90°, and the second link movement is induced by the movement of the first link-as its angle is referenced to the first link, it stabilizes on line with it, which is why it reaches 0°.The model speed response is illustrated in Figure 6.
Using Figure 7, it is easy to see that the nonlinear system has two equilibrium points at (90, 0) and (0, 0); this is what allowed us to define and apply the linearization method mentioned before.
As the purpose of the proposed model is to be controlled, to check if the input variable of the system affects its stabilizing point for a possible position control, different torque values were applied to each link.Tests were done by applying torques to the first link in between 0 and 18 Nm while the applied torque to the second link was 0 Nm; for the second link, due to its low mass, torques between 0 and 2.5 Nm were applied while the applied torque to the first link was 0 Nm.Results are shown in Figures 8 and 9.
For the first link, the test torque vector is defined as Equation (30); each element in the vector corresponds to an applied torque at each joint: = [0 2 4 6 8 10 12 14 16 18] = 0 (29) After applying this vector to the system, is clear that this input modifies the model steady state value, which is as expected for control purposes.Furthermore, it also affects the oscillation behavior of the system, making it slower and decreasing its amplitude (see Figure 8).The model behavior was as expected; the first link falls and moves describing a damped sine wave, which loses energy due to friction causing it to stabilize in −90 • , and the second link movement is induced by the movement of the first link-as its angle is referenced to the first link, it stabilizes on line with it, which is why it reaches 0 • .The model speed response is illustrated in Figure 6.
Using Figure 7, it is easy to see that the nonlinear system has two equilibrium points at (90, 0) and (0, 0); this is what allowed us to define and apply the linearization method mentioned before.
As the purpose of the proposed model is to be controlled, to check if the input variable of the system affects its stabilizing point for a possible position control, different torque values were applied to each link.Tests were done by applying torques to the first link in between 0 and 18 Nm while the applied torque to the second link was 0 Nm; for the second link, due to its low mass, torques between 0 and 2.5 Nm were applied while the applied torque to the first link was 0 Nm.Results are shown in Figures 8 and 9.For the first link, the test torque vector is defined as Equation (31); each element in the vector corresponds to an applied torque at each joint: τ 1 = 0 2 4 6 8 10 12 14 16 18 After applying this vector to the system, is clear that this input modifies the model steady state value, which is as expected for control purposes.Furthermore, it also affects the oscillation behavior of the system, making it slower and decreasing its amplitude (see Figure 8).For the second link, the test torque vector is defined as: = 0 = [0 0.1 0.3 0.6 0.8 1 1.2 1.6 1. 8 2] (30) Making some tests and applying different torque values to the second link, some instability points were encountered for input ( = ).The model response for the second link is shown in Figure 9, due to its nonlinearities, it is possible to observe that different answers are obtained with different behaviors, but after 1.2 Nm the system becomes completely unstable (see Figure 9).For the second link, the test torque vector is defined as: Making some tests and applying different torque values to the second link, some instability points were encountered for input u 2 (u 2 = τ 2 ).The model response for the second link is shown in Figure 9, due to its nonlinearities, it is possible to observe that different answers are obtained with different behaviors, but after 1.2 Nm the system becomes completely unstable (see Figure 9).
Making some tests and applying different torque values to the second link, some instability points were encountered for input ( = ).The model response for the second link is shown in Figure 9, due to its it is possible to observe that different answers are obtained with different behaviors, but after 1.2 Nm the system becomes completely unstable (see Figure 9).To check the behavior of the model when user characteristics vary, several tests were performed following the simulation scheme used before; a few of those tests are shown below.For this purpose, we refer to them as test 2 and test 3; test 2 was conducted by simulating a person of 65 Kg, 1.70 m tall, and test 3 simulated a person of 50 Kg, 1.50 m tall (See Figures 10-13).

Test 2:
Robotics 2017, 6, 20 14 of 25 To check the behavior of the model when user characteristics vary, several tests were performed following the simulation scheme used before; a few of those tests are shown below.For this purpose, we refer to them as test 2 and test 3; test 2 was conducted by simulating a person of 65 Kg, 1.70 meters tall, and test 3 simulated a person of 50 Kg, 1.50 m tall (See Figures 10-13).In the above-mentioned tests, many more nonlinearities were encountered when varying the user characteristics, which could be mitigated by using robust nonlinear controllers.Besides the nonlinearities, the system behavior was as expected and followed the results shown before.

Linear Model Simulation
Calculating all parameters and applying them to Equation (25), it is possible to obtain the following state-space representation of the linear approximation of the model: Matlab-Simulink was used to obtain the system time response; due to the linear approximation method used, made around the equilibrium point (Equation 24), the system reference changed.Thus, the system 0 • degrees are located at the negative vertical axes so initial conditions were needed, which were taken as θ 1 (0) = θ 2 (0) = 90 • .Results are shown in Figure 14.
Matlab-Simulink was used to obtain the system time response; due to the linear approximation method used, made around the equilibrium point (22), the system reference changed.Thus, the system 0° degrees are located at the negative vertical axes so initial conditions were needed, which were taken as (0) = (0) = 90°.Results are shown in Figure 14.The linear model differs from the nonlinear model response, but behaves as expected following the dynamics of a double pendulum stabilizing at 0° degrees.This because the linear model is merely an approximation made to work with the nonlinearities present in the original model.This model is affected in a higher way for the friction factor than the dynamical one; therefore, it loses energy faster and reaches the steady state faster.The model speed response is illustrated in Figure 15.The linear model differs from the nonlinear model response, but behaves as expected following the dynamics of a double pendulum stabilizing at 0 • degrees.This because the linear model is merely an approximation made to work with the nonlinearities present in the original model.This model is affected in a higher way for the friction factor than the dynamical one; therefore, it loses energy faster and reaches the steady state faster.The model speed response is illustrated in Figure 15.With this linear approximation, the phase diagram softens and converges into the steady state value more smoothly than the nonlinear one (see Figure 16).With this linear approximation, the phase diagram softens and converges into the steady state value more smoothly than the nonlinear one (see Figure 16).With this linear approximation, the phase diagram softens and converges into the steady state value more smoothly than the nonlinear one (see Figure 16).To check if the input variable of the system affects its stabilizing point, the same vectors of Equations ( 31) and (32) were applied to the linear model, the results of which are shown in Figures 17 and 18.To check if the input variable of the system affects its stabilizing point, the same vectors of Equations ( 31) and (32) were applied to the linear model, the results of which are shown in Figures 17 and 18.For both links, each input produces a significant change in the steady state of the model.It also reveals that the nonlinearity which produced the uncontrolled increase in the response of the system was eliminated, making the system proper for control design and tuning.
The results from tests 2 and 3 are shown below (Figures 19-22).For both links, each input produces a significant change in the steady state of the model.It also reveals that the nonlinearity which produced the uncontrolled increase in the response of the system was eliminated, making the system proper for control design and tuning.
The results from tests 2 and 3 are shown below (Figures 19-22).As expected, the linearized model behavior stabilizes at a different rate according to the mass and length of each segment of the leg.This speed change requires a robust speed control in order to avoid injuring a person with different characteristics when using the same device; it also could be configurable to enhance the therapy.As expected, the linearized model behavior stabilizes at a different rate according to the mass and length of each segment of the leg.This speed change requires a robust speed control in order to avoid injuring a person with different characteristics when using the same device; it also could be configurable to enhance the therapy.

Observability and Controllability Analysis
The linear model is planned to be used to design and tune linear control strategies, which is why it must be analyzed to ensure that it is controllable and/or observable.Defining the controllability matrix (Co) and observability matrix (O) as: where n is the range of A matrix.Using the linear state-space of Equation (33), which is in the form: .
it is easy to define that n = 4, so our analysis will be defined as: replacing the matrices of Equation (33) into Equations ( 34) and (35) will result in a controllability matrix of dimensions (4 × 8) and an observability matrix of dimensions (8 × 4), as the model is a multi-input multi-output (MIMO) system with two inputs and two outputs.To obtain the controllability matrix for each input-output pair, it is necessary to take columns (1-7) of Co to build the controllability matrix for the first pair (u 1 − y 1 ) and to take columns (2-8) to build the controllability matrix for the second pair (u 2 − y 2 ); in the case of observability, is the same process but selecting rows.The controllability and observability matrices will look like:  Nonsingular matrices indicate that the system is fully observable and controllable.

Exoskeleton Design Proposal
The model previously described is oriented to describe a simple exoskeleton with 2 DOF.Figures 23 and 24 show the general idea for the implementation; the exoskeleton consists of active and passive elements which are attached to the legs in a non-invasive way.The main idea is that the device can be wearable and adapted to the features of individual patients.The structure is intended to work as Master-Slave, where one leg controls the other one.This enables the user to use the non-injured leg to perform rehabilitation exercises that can be replicated in the injured one.The employed actuators are planned to be brushless DC motors coupled to a reduction gear (see Figures 25 and 26).The structure is intended to work as Master-Slave, where one leg controls the other one.This enables the user to use the non-injured leg to perform rehabilitation exercises that can be replicated in the injured one.The employed actuators are planned to be brushless DC motors coupled to a reduction gear (see Figures 25 and 26

Conclusions
The contribution of this paper focuses on extending the exoskeleton line of research by proposing a new way of modeling such that the followed objective is to develop a dynamic model using a combination of different models which mathematically describe the human body.The proposed model is a combination of three models, which, when united, can describe the physical behavior of the exoskeleton frame coupled to the human leg.This model has the advantage that it can be calculated for any kind of person given their its mass and height, converting it into a flexible alternative for calculating the exoskeleton dynamics very quickly and adapting them easily for a child's or young adult's body.As a secondary goal, the paper focuses on showing that the model can be controlled in different ways to satisfy the different characteristics of each possible rehabilitation therapy for which the exoskeleton could be used.As the resulting model is a highly nonlinear model,

Conclusions
The contribution of this paper focuses on extending the exoskeleton line of research by proposing a new way of modeling such that the followed objective is to develop a dynamic model using a combination of different models which mathematically describe the human body.The proposed model is a combination of three models, which, when united, can describe the physical behavior of the exoskeleton frame coupled to the human leg.This model has the advantage that it can be calculated for any kind of person given their its mass and height, converting it into a flexible alternative for calculating the exoskeleton dynamics very quickly and adapting them easily for a child's or young adult's body.As a secondary goal, the paper focuses on showing that the model can be controlled in different ways to satisfy the different characteristics of each possible rehabilitation therapy for which the exoskeleton could be used.As the resulting model is a highly nonlinear model, it was linearized and tested using simulation environments to ensure that it works as expected and its able to design and test control strategies.The linear representation results were fully controllable and observable, which fulfills the secondary objective proposed.Since this work was focused in the mathematical model description, and its validation was performed by software simulation, it is necessary to perform a validation with the physical device.As future work, tests with the physical device will be conducted and tests with human patients will then be considered so as to define a medical protocol to evaluate the device under the supervision of physical therapists.This model is also planned to be extended to three links, to be used to describe all leg movements in the sagittal plane.Also, control strategies including force feedback and/or joint virtual torque control with gravity compensation will be tested with the real equipment.

Figure 2 .
Figure 2. Body segment length according to D&C model.

Figure 2 .
Figure 2. Body segment length according to D&C model.

Figure 2 .
Figure 2. Body segment length according to D&C model.

Table 1 .
Generalized Hanavan model of the human leg.