Optimized Proportional-Integral-Derivative Controller for Upper Limb Rehabilitation Robot

: This paper proposes a nature inspired, meta-heuristic optimization technique to tune a proportional-integral-derivative (PID) controller for a robotic arm exoskeleton RAX-1. The RAX-1 is a two-degrees-of-freedom (2-DOFs) upper limb rehabilitation robotic system comprising two joints to facilitate shoulder joint movements. The conventional tuning of PID controllers using Ziegler-Nichols produces large overshoots which is not desirable for rehabilitation applications. To address this issue, nature inspired algorithms have recently been proposed to improve the performance of PID controllers. In this study, a 2-DOF PID control system is optimized o ﬄ ine using particle swarm optimization (PSO) and artiﬁcial bee colony (ABC). To validate the e ﬀ ectiveness of the proposed ABC-PID method, several simulations were carried out comparing the ABC-PID controller with the PSO-PID and a classical PID controller tuned using the Zeigler-Nichols method. Various investigations, such as determining system performance with respect to maximum overshoot, rise and settling time and using maximum sensitivity function under disturbance, were carried out. The results of the investigations show that the ABC-PID is more robust and outperforms other tuning techniques, and demonstrate the e ﬀ ective response of the proposed technique for a robotic manipulator. Furthermore, the ABC-PID controller is implemented on the hardware setup of RAX-1 and the response during exercise showed minute overshoot with lower rise and settling times compared to PSO and Zeigler-Nichols-based controllers.


Introduction
The life expectancy of people, and hence, the number of older adults, is increasing. Elderly people are most vulnerable to strokes. Stroke is among the main causes of limb disabilities and can be fatal [1,2]. According to statistics [3], more than 10 million people suffer from a stroke annually. Consequently, patients become dependent on others for their basic life activities. Physical therapy under the supervision of physiotherapists can help stroke patients to restore the functionality of their disabled limbs. However, traditional physical therapy involves manpower and is highly expensive as the number of patients increases. Rehabilitation robots have been introduced recently to reduce the burden on physiotherapists and increase the number of exercises performed during a therapy session. These robots can provide a more reliable service as they do not face monotony and fatigue failures due to the repetitive nature of the exercises. One of the most widely used control mechanisms employed in such robots is a proportional-integral-derivative (PID) controller. This mechanism is famous for its simple structure and robust performance in a wide range of operating conditions [4]. criteria. Four different objective functions, i.e., integral square error (ISE), integral time square error (ITSE), integral absolute error (IAE) and integral time absolute error (ITAE), have been investigated to find the controller with optimal or near optimal load disturbance response subject to robustness and maximum sensitivity constraints. Maximum sensitivity represents the inverse of the minimum distance on the Nyquist plot between critical point and loop transfer function. Such a method for tuning the controller parameters has proven to be effective for robust performance [21]. Furthermore, performance parameters such as overshoot, rise time, settling time and maximum sensitivity are normalized and the least average error (LAE) is evaluated. The optimal solution found for 2-DOF PID for RAX-1 is then implemented with the RAX-1 hardware for trials with three healthy subjects. The rest of the paper is organized as follows. Section 2 explains the controller design, Section 3 presents the simulation results, Section 4 provides the comparative analysis and discussion. Section 5 concludes the paper.

Methodology
RAX-1 is an exoskeleton device meant to be used for rehabilitation of upper limb extremities. Operating alongside the human arm, exoskeleton devices are required to produce movements similar to those performed by the upper limb. There are nine DOFs in the upper limb, excluding finger joints. This study focuses on the glenohumeral joint in the shoulder, which is a complex ball-and-socket joint that enables the shoulder to perform movements in three DOFs. These movements are commonly referred to as shoulder extension/flexion, abduction/adduction and medial/lateral rotation, also known as internal/external rotation. Figure 1 represents the three movements that can be performed with the shoulder joint. The ranges of motion for the shoulder joint movements performed by a healthy subject are listed in Table 1 [22]. These movement protocols are then implemented on  o find the controller with optimal or near optimal load disturbance response subject to robustness nd maximum sensitivity constraints. Maximum sensitivity represents the inverse of the minimum istance on the Nyquist plot between critical point and loop transfer function. Such a method for uning the controller parameters has proven to be effective for robust performance [21]. Furthermore, erformance parameters such as overshoot, rise time, settling time and maximum sensitivity are ormalized and the least average error (LAE) is evaluated. The optimal solution found for 2-DOF ID for RAX-1 is then implemented with the RAX-1 hardware for trials with three healthy subjects. he rest of the paper is organized as follows. Section II explains the controller design, Section III resents the simulation results, Section IV provides the comparative analysis and discussion. Section concludes the paper.

. Methodology
RAX-1 is an exoskeleton device meant to be used for rehabilitation of upper limb extremities. perating alongside the human arm, exoskeleton devices are required to produce movements similar o those performed by the upper limb. There are nine DOFs in the upper limb, excluding finger joints. his study focuses on the glenohumeral joint in the shoulder, which is a complex ball-and-socket oint that enables the shoulder to perform movements in three DOFs. These movements are ommonly referred to as shoulder extension/flexion, abduction/adduction and medial/lateral otation, also known as internal/external rotation. Figure 1 represents the three movements that can e performed with the shoulder joint. The ranges of motion for the shoulder joint movements erformed by a healthy subject are listed in Table 1 [22]. These movement protocols are then mplemented on RAX-1.

2.1.System Design
The robotic manipulator in the present study comprises two shoulder joints. T manipulator can be calculated by the numbers of links and joints. The 3D model illustrated in Figure 2.
A robot in the planner configuration is defined with three parameters (x, y, θ). H are three dimensional in real-world applications; hence, there is a need for six para yaw, pitch, roll) to describe the position and orientation of a robot in space.

System Design
The robotic manipulator in the present study comprises two shoulder joints. The DOF of the manipulator can be calculated by the numbers of links and joints. The 3D model of the robot is illustrated in Figure 2. The robotic manipulator in the present study comprises two shoulder joints. The DOF of the manipulator can be calculated by the numbers of links and joints. The 3D model of the robot is illustrated in Figure 2.
A robot in the planner configuration is defined with three parameters (x, y, θ). However, robots are three dimensional in real-world applications; hence, there is a need for six parameters (x, y, z, yaw, pitch, roll) to describe the position and orientation of a robot in space.  Figure 3 represents the direct kinematics of the robotic arm. Every rigid body in a serial chain has a label: Link 1 is the rigid body attached to the shoulder joint 1, Link 2 is the rigid body attached to Link 1 and so on. A joint is present between each link. Hence, Joint 1 attaches Link 1 to Link 0 and Joint 2 attaches Link 2 to Link 1. Frame of reference is numbered according to the respective links they are attached to, e.g. Frame 1 is attached to Link1. Eventually, the aim is to calculate the position of Frame 2 relative to that of Frame 0. A robot in the planner configuration is defined with three parameters (x, y, θ). However, robots are three dimensional in real-world applications; hence, there is a need for six parameters (x, y, z, yaw, pitch, roll) to describe the position and orientation of a robot in space. Figure 3 represents the direct kinematics of the robotic arm. Every rigid body in a serial chain has a label: Link 1 is the rigid body attached to the shoulder joint 1, Link 2 is the rigid body attached to Link 1 and so on. A joint is present between each link. Hence, Joint 1 attaches Link 1 to Link 0 and Joint 2 attaches Link 2 to Link 1. Frame of reference is numbered according to the respective links they are attached to, e.g. Frame 1 is attached to Link1. Eventually, the aim is to calculate the position of Frame 2 relative to that of Frame 0.  , and link (i) with their associated joints, joint (i-1), and joint (i). A frame (i) is assigned to link (i) as follows.
1. The zi-1 lies along the axis of motion of the ith joint. 2. The xi axis is normal to the zi-1 axis and pointing away from it. The Denavit-Hartenberg (DH) parameters of a rigid link depends on four geometric parameters (ai, αi, di, θi) [23]. The four parameters describe any revolute joint as follows: 1. ai (Link length) is a distance measured along the xi axis from the point of intersection of xi axis with zi-1 axis to the origin of frame (i). 2. αi (Link Twist) is the angle between the joint axes zi-1 and zi axes measured about xi axis in the right-hand orientation. 3. di (offset) is the distance measured along zi-1 axis from the origin of frame (i-1) to the intersection of xi axis with zi-1 axis. 4. θi (Joint angle) is the angle between xi-1 and xi axes measured about the zi-1 axis in the right-hand sense.  Figure 3a shows a pair of adjacent links which are link (i − 1), and link (i) with their associated joints, joint (i − 1), and joint (i). A frame (i) is assigned to link (i) as follows.

1.
The z i − 1 lies along the axis of motion of the ith joint.

2.
The x i axis is normal to the z i − 1 axis and pointing away from it.
The Denavit-Hartenberg (DH) parameters of a rigid link depends on four geometric parameters (a i , α i, d i , θ i ) [23]. The four parameters describe any revolute joint as follows:

1.
a i (Link length) is a distance measured along the x i axis from the point of intersection of x i axis with z i − 1 axis to the origin of frame (i). 2. α i (Link Twist) is the angle between the joint axes z i − 1 and z i axes measured about x i axis in the right-hand orientation.

3.
d i (offset) is the distance measured along z i − 1 axis from the origin of frame (i − 1) to the intersection of x i axis with z i − 1 axis.

4.
θ i (Joint angle) is the angle between x i − 1 and x i axes measured about the z i − 1 axis in the right-hand sense.
The 2-DOF upper limb robotic manipulator can be calculated from Figure 3b. The transformation matrix from Frame 0 to the end-effector can be defined as: where a 1 = 59.4 cm, a 2 = 105.8 cm and d 1 = 20 cm. The transformations can be used to determine kinematic measurements of the joints. For any joint angle, the position of the end effector can be derived from the transformation matrix.

Dynamic Model
In this study, Euler-Lagrangian approach was applied to calculate the dynamics of the robot manipulator. This approach uses the joint velocities and position to determine the kinetic and potential energies of a system. It generalizes Newtonian mechanics for systems that are subject to a specific class of constraints. These constraints are often expressed in terms of the position or variables describing the system in question.
The Lagrangian equation of motion defined in (1) is written as: where τ j denotes the required torque, L = K − P (Kinetic and potential energies) is the Lagrangian and qj is the generalized coordinate of the jth joint of robot. The inertial matrix D(q) can be determined as follows: where m 1 = 1.5 kg, m 2 = 0.5 kg. Simplifying (2), one can obtain the following (3).
The correction term Christoffel symbols ensures that when the derivatives of the vector field lying in a tangent plane of the configuration manifold are computed, they stay in the same tangent space. The Christoffel symbols in (4) are defined as The potential energy of each joint P i , is the product of the mass of that link m i , position vector to the centre of mass r ci and acceleration due to gravity g: The term φ k is a function of generalized coordinates that does not depend on their derivatives. It is given by the partial derivative of potential energy of the system with respect to the generalized coordinates as follows: Finally, the dynamic equations of the system after substituting various quantities and omitting zero can be expressed as d 11 ..
which, in general can be written in matrix form as: Here, D(q) denotes the inertia matrix of the system and C q, . q gives the Christoffel symbols and g(q) is actually φ k which is determined by taking partial derivative of potential energy with generalized coordinates. τ is a 2 × 1 matrix representing the generalized active forces.

Linearized Model
The linearized state space model for robot exoskeleton (RAX-1) is expressed as follows.

Motor Model
The robot manipulator requires actuators to provide the desired amount of torque at the joints. The actuators convert electrical energy into rotational mechanical energy. DC motors are widely used in robotics as actuators due to their high torque, speed controllability and portability [24]. The internal model of DC motor is illustrated in Figure 4 and can be expressed as follows.

2.4.Motor Model
The robot manipulator requires actuators to provide the desired amount of torque at the joints. The actuators convert electrical energy into rotational mechanical energy. DC motors are widely used in robotics as actuators due to their high torque, speed controllability and portability [24]. The internal model of DC motor is illustrated in Figure 4 and can be expressed as follows.
Here,' ' is the motor inertia, ' ' motor damper, ' ' is the motor constant, ' ' is the proportionality constant between angular velocity of the motor and back emf, whereas ' ' and ' ' inductance and resistance of the armature. If ≪ then, an approximated transfer function of motor is obtained by setting = 0. This converts motor model from a second order system to 1 st order system. DC motors with harmonic gears are used in this system. Table 2 lists the motor parameters.

The Exoskeleton Platform
In this section, the hardware design for the upper limb extremity rehabilitation is described. The platform consists of the design and manufacturing of the mechanical structure, actuators and sensors with hardware implementation of the control algorithm. Here,'J' is the motor inertia, 'B m ' motor damper, 'K τ ' is the motor constant, 'K b ' is the proportionality constant between angular velocity of the motor and back emf, whereas 'L a ' and 'R a ' inductance and resistance of the armature. If L a R a then, an approximated transfer function of motor is obtained by setting L a = 0. This converts motor model from a second order system to 1st order system. DC motors with harmonic gears are used in this system. Table 2 lists the motor parameters.

The Exoskeleton Platform
In this section, the hardware design for the upper limb extremity rehabilitation is described. The platform consists of the design and manufacturing of the mechanical structure, actuators and sensors with hardware implementation of the control algorithm.

Mechanical Structure
The 2-DOF mechanical platform is built using aluminum grade 6061 alloy and weighs approximately 15 kg. The structure is specifically designed to focus on specific exercises for rehabilitation of the shoulder joint ( Figure 5).
The upper extremity exoskeleton device consists mainly of a support frame, height adjustment mechanism and shoulder actuation mechanism. The support frame is attached to wheels enabling the platform to be remote. The manually controlled height regulation mechanism is attached to the support frame. The actuation mechanism is attached to the height regulation mechanism; it adjusts the mechanical frame fixed to human arm to match requirements to the subject height. The human arm is fixed to the exoskeleton with soft wraps both on the forearm and bicep. Other relevant adjustments are made during gait training. Joints at the wrist and elbow are passive and their orientations are fixed. the platform to be remote. The manually controlled height regulation mechanism is attached to the support frame. The actuation mechanism is attached to the height regulation mechanism; it adjusts the mechanical frame fixed to human arm to match requirements to the subject height. The human arm is fixed to the exoskeleton with soft wraps both on the forearm and bicep. Other relevant adjustments are made during gait training. Joints at the wrist and elbow are passive and their orientations are fixed.

3.2.Actuators and Sensors
The two-revolute joint system is equipped with the MAXON EC-90 motor (Maxon Motor, Sachseln, Switzerland) attached to the gear. A brushless flat motor has the maximum angular velocity of 2590 rpm and power rating of 90W which can generate maximum torque of 444 mNm. It is equipped with three Hall effect sensors to measure current and velocity and an internal encoder generating 2040 pulses per rotation. The reducers (harmonic drives) combined with flat motor have speed ratio of 160:1 for the first shoulder joint responsible for the internal/ external movement, while the ratio is 150:1 for the other joint responsible for the extension/ flexion movement. The detailed specification of other motor parameters and gear ratio is provided in Table 3. A force sensor is attached to the wrist handle of the robot, from which external disturbance exerted by the subject is measured while the exoskeleton is working in passive mode. The electrical setup of the system involves an ESCON 50/5 module which acts like a closed-loop speed or current controller for the motor.

3.3.Control Implementation
In the control system, a host computer is used as a graphical user interface. The software is written in visual C#; it allows the user to select one of the possible working modes of the exoskeleton

Actuators and Sensors
The two-revolute joint system is equipped with the MAXON EC-90 motor (Maxon Motor, Sachseln, Switzerland) attached to the gear. A brushless flat motor has the maximum angular velocity of 2590 rpm and power rating of 90W which can generate maximum torque of 444 mNm. It is equipped with three Hall effect sensors to measure current and velocity and an internal encoder generating 2040 pulses per rotation. The reducers (harmonic drives) combined with flat motor have speed ratio of 160:1 for the first shoulder joint responsible for the internal/ external movement, while the ratio is 150:1 for the other joint responsible for the extension/ flexion movement. The detailed specification of other motor parameters and gear ratio is provided in Table 2. A force sensor is attached to the wrist handle of the robot, from which external disturbance exerted by the subject is measured while the exoskeleton is working in passive mode. The electrical setup of the system involves an ESCON 50/5 module which acts like a closed-loop speed or current controller for the motor.

Control Implementation
In the control system, a host computer is used as a graphical user interface. The software is written in visual C#; it allows the user to select one of the possible working modes of the exoskeleton (passive, active or semi-active). The software enables uploading and downloading subject's history over a cloud specifically designed for this device. The user can also define the range of motion at which the exoskeleton should work for several repetitions.
A master/slave network is designed to interconnect the user with the system, where each joint in the system is a slave. CC3200 from Texas Instrument is used as a peripheral device that controls the angular movement. The master is instructed by the user to select an exercise and the number repetitions, as well as set an angular movement via an interface. This information is then communicated to the slaves for further implementation of the task [25] as depicted in Figure 6a,b. The overall system block diagram is illustrated in Figure 6c.
which the exoskeleton should work for several repetitions.
A master/slave network is designed to interconnect the user with the system, where each joint in the system is a slave. CC3200 from Texas Instrument is used as a peripheral device that controls the angular movement. The master is instructed by the user to select an exercise and the number repetitions, as well as set an angular movement via an interface. This information is then communicated to the slaves for further implementation of the task [25] as depicted in Figure 6(a) and 6(b). The overall system block diagram is illustrated in Figure 6(c).

4.1.Proportional Integrator Derivative (PID) Controller
One of the most widely used controllers is the PID controller, which has a very simple structure and is robust in under wide range of operating conditions. The PID error signal, which is the difference between setpoint and measured variable, is calculated and fed back into the system continuously to change the proportional, integral and derivative gains accordingly [26].

Proportional Integrator Derivative (PID) Controller
One of the most widely used controllers is the PID controller, which has a very simple structure and is robust in under wide range of operating conditions. The PID error signal, which is the difference between setpoint and measured variable, is calculated and fed back into the system continuously to change the proportional, integral and derivative gains accordingly [26].
PID is one of the most used control algorithms in industrial control systems. The response of a system is categorized by the rise time, overshoot, settling time and steady state error. In this study, a 2-DOF PID controller is used, and is represented mathematically according to Equation (11). The controller is capable of rejecting disturbances without significant increase of overshoot in setpoint tracking. It includes setpoint weighting on the proportional and derivative terms. A typical 2-DOF PID is composed of feed-forward and feedback compensators. The feed-forward compensator consists of a PD component while the feedback compensator includes PID component as shown in Figure 8.
where u denotes the input given to the plant or system, while K p , K d , K i , denotes the proportional, derivative and integral gains. In this study, 2-DOF PID in MATLAB is used, which has three more parameters to tune; b, c and N, where; b and c denote the setpoint weights, while N denotes a filter coefficient.

4.1.Proportional Integrator Derivative (PID) Controller
One of the most widely used controllers is the PID controller, which has a very simple structure and is robust in under wide range of operating conditions. The PID error signal, which is the difference between setpoint and measured variable, is calculated and fed back into the system continuously to change the proportional, integral and derivative gains accordingly [26].

Particle Swam Optimization (PSO)
PSO is an optimization technique inspired by the social behavior found in nature, such as flocking of birds and fish schooling. In the PSO search space, each solution acts like a flying bird in the search of food, known as a particle. PSO works based on the social behavior of particles in a swarm. This algorithm locates and finds the global best solution by adjusting an objective function. At the end of the process, best solution based on objective function is found for 2-dof pid controller parameters.
First, PSO chooses random solutions to initialize the population. Then, it updates its performance to obtain optimum value. Every particle is characterized by its position and velocity in the swarm. The velocity of a moving particle depends upon the change in position or direction [27]. Each particle updates its new position based on two components, p best and g best , where p best is the best position attained by a particle, while g best is the global best position of the entire swarm. The velocity and position of the particle can be expressed according to Equations (10) and (11), respectively.
where; v ij (t) denotes the particle velocity, x ij (t) denotes the particle position, v ij (t + 1) denotes the particle updated velocity, x ij (t + 1) denotes the particle updated position; p ij (t) denotes the particle best position; g j (t) denotes the global best position of the swarm; w , (w = w max − w min ) denotes the inertia term; r 1 and r 2 are two uniformly distributed random numbers ranging from 0 to 1 and c 1 and c 2 are the acceleration coefficients. The important steps of the PSO are summarized in Figure 9.
According to the PSO algorithm, the swarm size, position, velocity and the constants w, c 1 and c 2 are initialized first. Then, the fitness value of each particle is calculated, P best and g best are defined and the position and velocity of each particle are updated. The algorithm stops when the stopping criterion is met. The optimal solution is chosen according to the latest g best . The PSO parameters and their values used in our study are provided in Table 3.

Artificial Bee Colony (ABC)
ABC is an optimization technique based on the smart behavior of honey bees [16], which are famous for their intelligence and ability to perform complex tasks such as collecting nectar and building nests with a high degree precision [28]. Information about the quality of a food source is communicated within the colony by a particular dance language. The precision of the foraging range of honeybees allows an efficient exploitation of food sources and concentration of foraging on the best patches. There are two main concepts that describe swarm intelligence, namely self-organization and labor division [29]. The self-organizing behavior represents the complex collective behavior that rises from the local interaction between the agents showing a simple self-directed behavior. The mechanism of labor division assigns specific tasks to the agents performing simultaneous activities, which results in a more efficient and time-saving performance. According to the PSO algorithm, the swarm size, position, velocity and the constants , and are initialized first. Then, the fitness value of each particle is calculated, and are defined and the position and velocity of each particle are updated. The algorithm stops when the stopping criterion is met. The optimal solution is chosen according to the latest . The PSO parameters and their values used in our study are provided in Table 4. Table 4. PSO Parameters.

Parameter
Value Number of particles 20 There are three categories of artificial bees in ABC algorithm: employed, scout and onlooker bees. The employed bees go to the food source that has already been visited by them, while unemployed bees look for further sources of food. The number of employed bees is equal to that of food sources. Searching for new sources is performed by the scout bees, whereas the onlooker bees wait for the information about the discovered food sources provided by employed bees via their waggle dance. If the position of a food source does not improve within a number of attempts known as limits, then the employed bees become scout bees. In this manner, the exploitation process is performed by employed and onlooker bees, whereas the scout bees explore the existing solutions. The details of different ABC phases are described in the following subsections. a

. Initialization Phase
Random initialization of the locations of food sources is performed according to the following equation: where SN is the number of food sources taken as half of the bee colony, D is the dimension of the problem, x ij denotes the i th employed bee on j th dimension, wjile x max j and x min j denotes its upper and lower bounds.
b. Employed Bee Phase Every other employed bee is allocated with the food source for further exploitation. Equation (11) represents the process of generating a food source.
where φ = a × rand(0, 1) is a random variable denoting the acceleration coefficients ranging in the interval [−1,1] and v ij is the new solution or a food source. The fitness f it i of the new food source is calculated according to the following equation: where, f i is the objective function of each food source. A selection is made between the original and new food sources to choose the better one by keeping the fitness value in accordance.

c. Probabilistic Selection Phase
An onlooker bee selects a food source with a certain probability calculates as where p i denotes the probability of selecting the ith solution.

d. Onlooker Bee Phase
The employed bees share the information about food sources with the onlooker bees, who select a food source to exploit better solutions according to its selection probability. The fitness values of each exploited food sources is calculated. A greedy selection between original and new food sources is made similar to the employed bees phase.
e. Scout Bee Phase If a food source does not yield better results within the limits L, where L = 0.6 × (Number o f optimization parameters) × (colony size), then this food source is abandoned and the bee associated with it becomes a scout bee. In this case, a new source of food is randomly generated according to Equation (15). All phases will continue until the termination criterion is met. The output is the best food source solution. Figure 10 shows the flowchart of the ABC algorithm which depicts the process of 2-dof pid controller optimization.
According to the ABC algorithm, the colony size, position and the constants L and a are initialized first. Then, the fitness value is calculated, and the best food source is defined. The algorithm stops when the stopping criterion is met. The optimal solution is chosen according to the latest g best . The ABC parameters and their values used in our study are provided in Table 4. According to the ABC algorithm, the colony size, position and the constants and are initialized first. Then, the fitness value is calculated, and the best food source is defined. The algorithm stops when the stopping criterion is met. The optimal solution is chosen according to the latest . The ABC parameters and their values used in our study are provided in Table 5.

Results and Discussions
Based on the simulation, the system is divided into two linearized sub systems representing Joint 1 and Joint 2 [30]. Independent controls for both joints with saturation limits are implemented. A step input is given to the model of the robot and the response is observed. A schematic of the proposed controller's tuning method based on PSO and ABC is shown in Figure 11.

Results and Discussions
Based on the simulation, the system is divided into two linearized sub systems representing Joint 1 and Joint 2 [30]. Independent controls for both joints with saturation limits are implemented. A step input is given to the model of the robot and the response is observed. A schematic of the proposed controller's tuning method based on PSO and ABC is shown in Figure 11. For this research, a comparative study was carried out with four different cost functions to gauge the appropriate objective function for this study. A well-chosen objective function leads to better performance of the system and meets control design expectations. These objective functions include the integral squared error (ISE), integral time squared error (ITSE), integral absolute error (IAE) and integral time absolute error (ITAE), and are defined as follows: = (18) = | | (19) = . (20) = . | |

5.1.Robustness Consideration
Focus on the tuning procedure requires some robustness considerations in the design. This is achieved by using the maximum sensitivity function as a measure of robustness and is given by where | ( )| . The sensitivity function shows the effect of feedback on the output. Disturbances are attenuated if | ( )| < 1 and are amplified if | ( )| 1. The robustness of the closed loop increases with the decrease in . The values for ranging from 1.2 to 2.0 provides reasonable robustness [27].

5.2.Simulation Setup
Before implementing ABC-PID [11] on hardware, simulations are carried out to validate and verify the control algorithm. In this section, a comparative analysis of the PID controller optimized with ABC, PSO, and conventional tuning method based on Zeigler-Nichols is presented. The parametric bounds are provided in Table 6. The Zeigler-Nichols PID parameters and PID parameters For this research, a comparative study was carried out with four different cost functions to gauge the appropriate objective function for this study. A well-chosen objective function leads to better performance of the system and meets control design expectations. These objective functions include the integral squared error (ISE), integral time squared error (ITSE), integral absolute error (IAE) and integral time absolute error (ITAE), and are defined as follows:

Robustness Consideration
Focus on the tuning procedure requires some robustness considerations in the design. This is achieved by using the maximum sensitivity function as a measure of robustness and is given by

Simulation Setup
Before implementing ABC-PID [11] on hardware, simulations are carried out to validate and verify the control algorithm. In this section, a comparative analysis of the PID controller optimized with ABC, PSO, and conventional tuning method based on Zeigler-Nichols is presented. The parametric bounds are provided in Table 5. The Zeigler-Nichols PID parameters and PID parameters tuned using ABC and PSO used for simulation are provided in Tables 6 and 7. The simulations presented here are for internal/external rotation, extension/flexion and abduction/adduction of the shoulder joint.  To evaluate the fitness of the objective function, normalization of the objective function is carried out, which scales objective function within a specified range. The following normalization function is used [31].
where f it i represents the objective to be normalized and f it overall represents the overall fitness. δ is kept 0.99 to avoid any zeros during normalization. Furthermore, normalized average has been evaluated to determine the significance of the objective function. Table 8 shows that IAE and ITAE used as objective functions produce minimal overshoots and suitable rise and settling times for both PSO-PID and ABC-PID. When using ISE and ITSE, the rise time is small, while larger overshoots are detected. However, it is obvious that the performance parameters of ABC-PID for Joints 1 and 2 are much better in terms of the six performance parameters mentioned above compared to that of PSO-PID or Zeigler-Nichols. The results plotted in Figure 12 show the response for Joint 1 for PID controller tuned using the Zeigler-Nichols and optimal parameters found with PSO and ABC for all four-objective functions. The figure shows that the ZN response is quicker with a higher rise time, higher settling time and larger overshoot. The results also show a very low percentage of overshoot in the response with a very low steady state error for all four objective functions. However, the optimal response of the system is found while using IAE as objective function. The results plotted in Figure 12 show the response for Joint 1 for PID controller tuned using the Zeigler-Nichols and optimal parameters found with PSO and ABC for all four-objective functions. The figure shows that the ZN response is quicker with a higher rise time, higher settling time and larger overshoot. The results also show a very low percentage of overshoot in the response with a very low steady state error for all four objective functions. However, the optimal response of the system is found while using IAE as objective function. Figure 12. Step response of Joint 1 using ISE, IAE, ITSE and ITAE. Figure 13 shows the step response for Joint 2. It is evident from the figure that PID tuned Zeigler-Nichols produces large overshoot than the PID optimized with ABC or PSO. However, the overshoot of ABC is 4.4% which is smaller as compared to 22% overshoot of PSO when using IAE. It was found that PID-ABC has the minimum average of the normalized objective function for Joint 1 and Joint 2. Hence, the ABC-optimized PID controller is found to be robust and stable for practical implementation on a robotic arm manipulator. The performance and robustness obtained by calculating the system percent overshoot, rise time, settling time, cost and sensitivity is presented in Table 9. Step response of Joint 1 using ISE, IAE, ITSE and ITAE. Figure 13 shows the step response for Joint 2. It is evident from the figure that PID tuned Zeigler-Nichols produces large overshoot than the PID optimized with ABC or PSO. However, the overshoot of ABC is 4.4% which is smaller as compared to 22% overshoot of PSO when using IAE. It was found that PID-ABC has the minimum average of the normalized objective function for Joint 1 and Joint 2. Hence, the ABC-optimized PID controller is found to be robust and stable for practical implementation on a robotic arm manipulator. The performance and robustness obtained by calculating the system percent overshoot, rise time, settling time, cost and sensitivity is presented in Table 8. Furthermore, an investigation on closed-loop stability analysis is performed by Nyquist plot of the controllers tuned with ABC using IAE for the given systems. The plot for Joint 1 and Joint 2 are shown in Figure 14. It indicates that the both systems are asymptotically stable as no poles lie in the right half of the plane. Furthermore, an investigation on closed-loop stability analysis is performed by Nyquist plot of the controllers tuned with ABC using IAE for the given systems. The plot for Joint 1 and Joint 2 are shown in Figure 14. It indicates that the both systems are asymptotically stable as no poles lie in the right half of the plane. Figure 13. Step response of Joint 2 using ISE, IAE, ITSE and ITAE.
Furthermore, an investigation on closed-loop stability analysis is performed by Nyquist plot of the controllers tuned with ABC using IAE for the given systems. The plot for Joint 1 and Joint 2 are shown in Figure 14. It indicates that the both systems are asymptotically stable as no poles lie in the right half of the plane.

Experimental Evaluation for RAX-1
A preliminary experiment was conducted to test the efficacy and performance of the exoskeleton controlled using the optimized PID controller. Three healthy male subjects participated in the experiment and their properties are listed in Table 10. Subjects were bound with soft wraps and their preliminary tests were performed before and after the exercise. The tests included measuring the body temperature, heart rate and oxygen levels. The evaluation procedure was approved by the Ethics committee of the University of Kuala Lumpur (UniKL), Kuala Lumpur, Malaysia. All subjects were volunteers who signed consent forms before the experiment. The procedure was completed in the presence of physiotherapists from the Royal College of Medicine Perak, Malaysia. The mechanical structure used for this experiment was RAX-1 illustrated in Figure 5. In this experiment, the upper extremity exoskeleton, which is a robot in charge of the rehabilitation protocol, provides passive assistance for the subjects. The two actuators follow the specific trajectories of a pre-defined range of motions. The followed trajectory and joint angles can be

Experimental Evaluation for RAX-1
A preliminary experiment was conducted to test the efficacy and performance of the exoskeleton controlled using the optimized PID controller. Three healthy male subjects participated in the experiment and their properties are listed in Table 9. Subjects were bound with soft wraps and their preliminary tests were performed before and after the exercise. The tests included measuring the body temperature, heart rate and oxygen levels. The evaluation procedure was approved by the Ethics committee of the University of Kuala Lumpur (UniKL), Kuala Lumpur, Malaysia. All subjects were volunteers who signed consent forms before the experiment. The procedure was completed in the presence of physiotherapists from the Royal College of Medicine Perak, Malaysia. The mechanical structure used for this experiment was RAX-1 illustrated in Figure 5. In this experiment, the upper extremity exoskeleton, which is a robot in charge of the rehabilitation protocol, provides passive assistance for the subjects. The two actuators follow the specific trajectories of a pre-defined range of motions. The followed trajectory and joint angles can be measured from the encoders. Figure 15 shows a subject performing the shoulder exercises with the assistance of exoskeleton device. The angular rotation of the joint can be set via the interface, and its current position can be observed on the screen.
A ramp/ sinusoidal repetitive input is provided to the shoulder joint to perform three different exercises, namely shoulder extension/flexion, internal/external rotation and abduction/adduction (Table 1).
For this experiment, a combined shoulder external/internal rotation was performed by the first subject; the upper bound movement during external rotation was limited to 90 • and that during internal rotation was limited to 0 • . Figure 16a illustrates the motion tracking of the robotic arm where a minute overshoot was observed in reaching the target reference. The graph also shows the error representation between the reference and actual response and speed of the motor in terms of rpm. During the experiment, each subject went through the same gait pattern, i.e. the movement was set to a frequency of 0.25 Hz for internal/external rotation. measured from the encoders. Figure 15 shows a subject performing the shoulder exercises with the assistance of exoskeleton device. The angular rotation of the joint can be set via the interface, and its current position can be observed on the screen. A ramp/ sinusoidal repetitive input is provided to the shoulder joint to perform three different exercises, namely shoulder extension/flexion, internal/external rotation and abduction/adduction (Table 1). Figure 15. A subject performing exercise with assistance of exoskeleton device.
For this experiment, a combined shoulder external/internal rotation was performed by the first subject; the upper bound movement during external rotation was limited to 90° and that during internal rotation was limited to 0°. Figure 16(a) illustrates the motion tracking of the robotic arm where a minute overshoot was observed in reaching the target reference. The graph also shows the error representation between the reference and actual response and speed of the motor in terms of rpm. During the experiment, each subject went through the same gait pattern, i.e. the movement was set to a frequency of 0.25 Hz for internal/external rotation. Figure 16(b) shows the current required by the motor to move the shoulder to a desired angle. It also represents the amount of external force exerted by the patient on the system. The external force acts as an external disturbance to the system and its effect can be seen in the Figure 16(b). As the external force to resist the movement increases, the motor current drastically increases to overcome the effect of disturbance. The maximum current drawn by the motor is 5.27 A when the maximum external resisting force of 68 N is applied to the joint. The second exercise, namely shoulder abduction/adduction, was performed by the second subject; the bounds are set from 0° to 90°. Figure 17(a) represents the response of the system during shoulder abduction and adduction. It also represents the difference between the reference signal and response, which demonstrates a minute overshoot in the response. The gait pattern for this exercise is similar to that for the previous exercise, where the movement is set to a frequency of 0.25 Hz.  Figure 16b shows the current required by the motor to move the shoulder to a desired angle. It also represents the amount of external force exerted by the patient on the system. The external force acts as an external disturbance to the system and its effect can be seen in the Figure 16b. As the external force to resist the movement increases, the motor current drastically increases to overcome the effect of disturbance. The maximum current drawn by the motor is 5.27 A when the maximum external resisting force of 68 N is applied to the joint.
The second exercise, namely shoulder abduction/adduction, was performed by the second subject; the bounds are set from 0 • to 90 • . Figure 17a represents the response of the system during shoulder abduction and adduction. It also represents the difference between the reference signal and response, which demonstrates a minute overshoot in the response. The gait pattern for this exercise is similar to that for the previous exercise, where the movement is set to a frequency of 0.25 Hz.  Figure 17(b) represents current drawn by the motor to achieve the target angle and the external force applied by the subject. When external resistance was exerted on the system by the subject, an unusual activity in the current was detected to overcome the external disturbance. A maximum current of 5.14 A was observed with 50N of applied force. The impact of this resistance can also be observed in Figure 17(a), where disruption occurs when external force is applied while achieving the reference.  Figure 17b represents current drawn by the motor to achieve the target angle and the external force applied by the subject. When external resistance was exerted on the system by the subject, an unusual activity in the current was detected to overcome the external disturbance. A maximum current of 5.14 A was observed with 50N of applied force. The impact of this resistance can also be observed in Figure 17a, where disruption occurs when external force is applied while achieving the reference.
The third exercise, namely, shoulder extension/flexion, was performed by the third subject; the range of motion was fixed to oscillate between 0 • to 90 • . In this experiment a pulsating input was fed as a reference (see Figure 18).  Figure 18(a) shows the motion tracking of the shoulder joint while performing extension/flexion. It also shows the error between reference and actual angle with the speed of the motor. In this exercise, each subject went through same gait pattern i.e. the exoskeleton is set to move with the frequency of 0.25 Hz for extension/flexion. Figure 18(b) shows the current driven by the motor to perform the exercise and force exerted by the patient to resist the movement. It can be seen that the motor draws more current than usual where the external force is applied; its impact on the movement  Figure 18a shows the motion tracking of the shoulder joint while performing extension/flexion. It also shows the error between reference and actual angle with the speed of the motor. In this exercise, each subject went through same gait pattern i.e. the exoskeleton is set to move with the frequency of 0.25 Hz for extension/flexion. Figure 18b shows the current driven by the motor to perform the exercise and force exerted by the patient to resist the movement. It can be seen that the motor draws more current than usual where the external force is applied; its impact on the movement can also be observed in Figure 18a. It is evident from the graph that a small amount of overshoot with no steady state was observed. The maximum force applied by the subject was recorded at 55N, while 4.5A of current was drawn by the motor.

Comparison of the PID Controllers Optimized with ABC and the Zeigler-Nichols Method
An experiment was performed to compare PID tuned controllers tuned with ABC and the Zeigler-Nichols method under the effect of disturbance. The parameters set for both techniques are similar to those defined in the simulation setup. The experiment was performed with the extension/flexion movement and the range of motion for both techniques was the same. The response of the system using ABC optimization and the Zeigler-Nichols method is shown in Figure 19. It can be seen from Figure 19 that PID controller tuned with the Ziegler-Nichols has a low settling time and quick response. However, the response of the system is not smooth and has a steady state error. At the same time, ABC-optimized PID response is slower due to a greater settling time. However, it has a smooth response with virtually no steady-state error. The angular rotation of the arm is set to move from 0° to 90°. Figure 19 Figure 19. (a) Response of the system with ABC optimized PID and tuned using Zeigler-Nichols, (b) Error of the system with ABC optimized PID and tuned using Zeigler-Nichols.
It can be seen from Figure 19 that PID controller tuned with the Ziegler-Nichols has a low settling time and quick response. However, the response of the system is not smooth and has a steady state error. At the same time, ABC-optimized PID response is slower due to a greater settling time. However, it has a smooth response with virtually no steady-state error. The angular rotation of the arm is set to move from 0 • to 90 • . Figure 19a illustrates the response of the system tuned with ABC and the Zeigler-Nichols method, while the error signal is shown in Figure 19b. It is evident from the graph that the rise time of ZN-PID is lower than ABC-PID. The overshoot and steady state produced by ZN-PID are also larger than those of ABC-PID. Zeigler-Nichols based PID produces an angular drift of 3.64 • while it is 0.63 • for ABC-PID. A slower response is favorable for the rehabilitation system used in this study. The presence of large overshoot and steady state in the system is dangerous, as it may cause dislocation or fracture of the arm. Hence, it can be concluded that PID controller works better when optimized with ABC compared to PSO or Zeigler-Nichols method as ABC-PID produces very low overshoot, slower rise time and no steady state.

Conclusions
Rehabilitation robotics have been studied for decades, but few researchers have ever considered using optimization techniques for gait training. This paper presents PSO-and ABC-based tuning techniques for a 2-DOF PID controller used in the robot controlling trajectories of different exercise movements performed by the shoulder joints, namely internal/external rotation, abduction/adduction, and extension/flexion. RAX-1 was used as a mechanical experimental platform. The control parameters of the PID controller were tuned using the ABC and PSO algorithms as well as the Zeigler-Nichols method. The simulation results demonstrated better feasibility of the proposed controller in terms of robustness. An ABC-optimized PID controller was also implemented into the hardware for three subjects, with each performing different exercises. It was necessary for the robot to perform steady motion with no steady state error during the process of rehabilitation, as any abrupt movement could dislocate the shoulder joint. The hardware results showed that the controller could trace the desired trajectory with a very minute overshoot in the response and significantly low response times, which are desirable in rehabilitation. Finally, this study compared the performance of PID controllers optimized with ABC and the Zeigler-Nichols method, respectively. The controller tuned using the Zeigler-Nichols method demonstrated a decent rise time but large overshoot in the response, which is dangerous for rehabilitation applications. However, the hardware response of the system had less overshoot and no steady-state error when the ABC optimizer was used to tune the PID parameters. In summary, this study focused on finding optimal parameters of the PID controller used in an upper limb rehabilitation robotic system. In the future, we plan to broaden the spectrum of this study by incorporating various protocols related to the elbow and wrist rehabilitation with an updated mechanical structure. Analysis with other notable optimization techniques such as firefly and ant colony algorithms will be further investigated to compare their parameters and validate their capabilities in rehabilitation applications.