Safety Operation of n-DOF Serial Hydraulic Manipulator in Constrained Motion with Consideration of Contact-Loss Fault

Featured Application: The suggested control scheme has widespread applications in ﬁelds of industrial manipulator and can be further expanded to research of automation. Abstract: In consideration of accidental contact-loss due to step-change or accidentally moving out of a constrained framework, this paper focuses on solving this problem during working processes of an n-degree-of-freedom hydraulic manipulator (n-DOF manipulator). In order to overcome this phenomenon, a fault detection methodology-based virtual energy tank is employed with a shaping function to prevent the end-e ﬀ ector from damage or unexpected motion. This technique helps to detect when the contact-loss happens by a virtual energy variable; thus, decoupling a force control regulation. Moreover, a new trajectory for smooth motion after contact-loss detection is also discussed to increase system robustness. Additionally, to enhance tracking performance, adaptive laws are designed to compensate for system uncertainties. Comparative simulations are given on the n-DOF hydraulic manipulator to evaluate e ﬀ ectiveness of the impedance-based energy tank methodology under the sudden step-changed environment. Moreover, inﬂuences of control gains and setup energy parameters to the system behaviors when contact-loss happens are remarkably discussed as indispensable criteria for further development. The simulated results certiﬁed the superior e ﬀ ectiveness and reliability of the suggested methodology over the conventional impedance control for safe operation.


Introduction
The robotic manipulator has currently been essential to widespread industrial applications and services [1][2][3]. Research fields of robotic manipulators are not only restricted to traditional electric manipulators but also extended to pneumatic and hydraulic types owing to advantages of high power-to-weight ratio and stiffness [4]. Different from human self-adaptivity activities, movements of the manipulator are mostly setup based on predetermined typically repeated tasks like loading, positioning. Various studies of algorithm developments have been reported to enhance system behaviors such as advanced adaptive control [5][6][7], optimal control [8], intelligent techniques [9][10][11][12], time-delayed estimation [13][14][15], observers [16][17][18], or even optimizations [19][20][21]. However, those studies intensively focused on the improvement of tracking performance and evaluation of system behaviors in free-motion.
where q, . q, .. q ∈ R n are position, angular velocity, and angular acceleration vectors of each joint, respectively. M(q) ∈ R n×n is the symmetric and positive definite matrix of nominal inertia. C q, . q ∈ R n×n denotes the nominal Coriolis and Centrifugal term matrix. G(q) ∈ R n is the nominal gravity term. ∆ denotes lumped parametric uncertainties, disturbances, and model error. F ext ∈ R n represents external load induced from environment interaction in the end-effector. τ ext ∈ R n is the torque vectors acting on joints. J(q) ∈ R n×n is the Jacobian matrix to transfer from Cartesian to manipulator joint space as x = f(q) . x = ∂f(q) ∂q = J(q) . q (2) Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 23 ( ) A L x (6) where β is the effective bulk modulus,  q X = 0 for any non-zero vector X ∈ R n . Property 2. The inertia matrix M(q) ∈ R n×n is a symmetric and positive definite matrix and is bounded as m l X 2 ≤ X T M(q)X ≤ M u X 2 , in which m l and M u are lower and upper bound positive constants, respectively.
Movements of each link are driven by an electro-hydraulic actuator. Herein, two type of the electro-hydraulic actuator are utilized as depicted in Figure 2. Due to its complicated characteristics, the actuator dynamics should be derived and included in the whole system dynamics. Define x a = x a1 x a2 . . . x an−1 q an ∈ R n is a vector of actuator positions in actuator-space. The actuator position can be relevant calculated from manipulator joint space as x a = f a (q) .
x a = J a (q) . q (3) ( ) where β is the effective bulk modulus,

Assumption 1.
Only internal leakage of the actuator is considered in this paper. The other parametric uncertainties or model errors are ignored. Moreover, the leakage is supposed to be proportional to load pressure as expressed in Equations (5) and (6).
These flows are expressed as [4] ( ) ( ) ω ρ ρ where S P and t P are the supply pressure and the returned-to-tank pressure, respectively. d C is the discharge coefficient. ω is the valve orifice area gradient, ρ is the density of oil, and vi x is the  Without losing generality, let us consider one element of the total robotic system. The actuator force/torque is calculated as where A r , A ai1 , and A ai2 are volumetric area of rotary actuator and areas of piston head part and rod part of cylinders, respectively. P ai1 and P ai2 (i = 1, 2, . . . , n) are pressures of two chambers of each actuator. These parameters are expressed as .
where β is the effective bulk modulus, V ij 0 (q) ( j=1,2) are initial volumes of the two chambers, C leak_i is coefficient of internal leakage, and Q aij (t) ( j=1,2) represents two flows into and out of the two chambers.

Assumption 1.
Only internal leakage of the actuator is considered in this paper. The other parametric uncertainties or model errors are ignored. Moreover, the leakage is supposed to be proportional to load pressure as expressed in Equations (5) and (6).
These flows are expressed as [4] Q ai1 = C d ωx vi where P S and P t are the supply pressure and the returned-to-tank pressure, respectively. C d is the discharge coefficient. ω is the valve orifice area gradient, ρ is the density of oil, and x vi is the motion of the spool of the valve. s(x vi ) = 1 if x vi > 0 and s(x vi ) = 0 when x vi ≤ 0 is a function to determine direction of the spool valve for obtaining flow rate through a proportional valve.

Remark 1.
Ignoring the time-delay phenomenon in valve dynamics and its uncertainties, the spool displacement of the proportional valve can be simplified as a linear function of applied input voltage u vi as x vi = K vi u vi with K vi is the valve gain.

Adaptive Backstepping Sliding Mode Algorithm
From the explanation in Introduction, in this section, the adaptive backstepping sliding mode control is designed to take into account for actuator dynamics to achieve global stability and robustness under the presence of parametric uncertainties or disturbances due to system model error. The overall control scheme is depicted in Figure 3.

Position regulation Torque control
Environment Modify trajectory  A new equivalent variable is defined as New estimated error between the actual and estimated lumped uncertainties is defined as = − ∈ ˆn R Δ Δ Δ in which the estimated lumped uncertainties satisfies: Taking derivative the Lyapunov candidate 1 V with Property 1 and substituting Equations (11) and (12) yields:  To conveniently derive the backstepping sliding mode algorithm, the dynamics of the manipulator including actuators is rewritten from Equations (1), (5)-(8) as y 3 = f(y 1 , y 2 ) + gu + ΨC leak (9) where y = y 1 y 2 y 3 T = q . q (A a1 P a1 − A a2 P a2 ) T is a new state vector of the system. y 1 = y 11 y 12 y 13 . . . y 1n T denotes vector of measurable position of joints (y 1i means angular position of joint ith). y 2 = y 21 y 22 y 23 . . . y 2n T denotes vector of measurable angular velocity (y 2i mean angular velocity of joint ith). f(y 1 , y 2 ), g, Ψ, C leak are defined as the following: Appl. Sci. 2020, 10, 8107

of 22
State error is defined as e 1 = y 1 − y 1d ∈ R n and e 2 = . e 1 = y 2 − y 2d ∈ R n where y 1d , y 2d ∈ R n are the vector of desired angles and angular velocities, respectively.
A new equivalent variable is defined as New estimated error between the actual and estimated lumped uncertainties is defined as ∆ = ∆ −∆ ∈ R n in which the estimated lumped uncertainties satisfies: Assumption 2. The lumped uncertainties are supposed to vary slowly relative to the system dynamics, i.e., . ∆ i 0.
Substituting the adaptive law in Equation (13) yields: .
Remark 2. The aim of using the adaptive law in Equation (13) is to cancel out influences of the lumped uncertainties to force the model as same as an ideal model. In practical simulation, to enhance the convergence rate of estimation, this equation is modified as where K ∆ , η ∆ are diagonal positive definite matrices to increase the convergence rate of the estimated lumped uncertainties.
Remark 3. The lumped uncertainty ∆ is unknown but bounded by ∆ ≤ ∆. Despite using the new adaptive law in Equation (19), an error of estimated lumped uncertainties is still existing, but its estimated error is bounded and predetermined as ∆ = ∆ −∆ ≤ ε ∆ . Then η 1 is chosen to be greater than the error of estimated lumped uncertainties, i.e., η 1i ≥ ε ∆i . This means the controller gains are much smaller instead of designing with large value enough to cover the upper bound of uncertainties in the conventional design.
From Equation (18), the derivative of the Lyapunov function V 1 is semi-negative if the error e 3 converges to zero. Therefore, this error should be tackled as small as possible.
Step 2: Define the sliding surface as Taking the sliding surface s 2 as derivative, one obtains: e 3 = f(y 1 , Choosing the second Lyapunov function V 2 as: where e C leak = C leak −Ĉ leak is an estimated error between a real and estimated leakage coefficient. Taking derivative of the Lyapunov function V 2 and substituting Equation (21) results in: .
To guarantee the system robustness and global stability, the control input signal u v is designed as where K 2 diag K 21 K 22 . . . K 2i . . . K 2n and η 2 diag η 21 η 22 . . . η 2i . . . η 2n are positive definite matrices. σ 2 is the width vector, the term −K 2 s 2 is added to enhance the convergence rate when s 2 is large. A function tanh , . . . , tanh Equations (16) and (24) are used to replace a conventional signum function to achieve chattering-free and facilitate control design. Hence, the derivative of the Lyapunov function V 2 becomes: .
Appl. Sci. 2020, 10, 8107 8 of 22 To achieve system performance as an ideal model, the term e T 3 Ψe leak + e T C leak . e C leak should be cancelled out. Therefore, the following adaptive law is satisfied: Remark 4. The aim of using the adaptive law in Equation (26) is to cancel out influences of the leakage to achieve the ideal model with no uncertainties or disturbances. In practical simulation, to enhance the convergence rate of estimation, this equation is modified as where K C leak and η C leak are a diagonal positive definite matrix to adjust the convergence rate of the estimation. A function sign(• i ) in both Equations (19) and (27) is defined as sign(

Remark 5.
Despite using the adaptive law in Equation (26), an error of estimated leakage coefficient is still existing, but its estimated error is bounded and predetermined as e leak = C leak −Ĉ leak ≤ ε leak . Then η 2 is chosen to be greater than the error of estimated leakage coefficient, i.e., η 2i ≥ ε leak_i . This means the controller gains are much smaller instead of designing with large values enough to cover the upper bound of uncertainties in conventional design.
Consequently, the derivative of Lyapunov function V 2 becomes: .
where ε lumped denotes lumped error of estimated procedures.
Remark 6. In a practical system, the position of the end-effector can be calculated from angular position of each joint. These values can be directly measured by using encoders. However, angular and acceleration are not available. Conventionally, these variables are obtained from taking derivative of joint angles and, thereby, measured noises are also amplified. To effectively estimate unmeasurable variables for designing the controller, a second-order exact differentiation (SOED) is utilized as where z 0 , z 1 , z 2 are estimated values of q, . q, ..

Remark 7.
It should be noted that all estimated errors of the SOED can achieve finite-time convergence no matter what a control input signal is. This implies that the SOED observer can be designed separately to the main controller and, thereby, a full state is available for designing control signals.

Remark 8.
Generally, increasing those observer gains κ 0 , κ 1 , κ 2 results in faster convergence. However, large values of observer gain also result in more noise influence. Too large or too small values of observer gains may lead to overshoot. Therefore, the values of all observer gains κ 0 , κ 1 , κ 2 should be chosen based on careful tuning. Adaptive laws for automatically adjusting observer gains are not considered as the main scope in this paper.

Fault Detection-Based Virtual Energy Tank
The conventional position-based impedance control is redesigned from [49] as where ∆x = x d − x re f and x d , x re f denote the desired and reference trajectories. D m , B m , K m are matrices representing impedance behaviors. u imp is generated signal from force regulation and is defined as with F ext , F d denote vectors of external force induced from interaction with the environment and desired force to be tracked, respectively. K p f , K i f are proportional and integrated controller gains, and the second term η f sign(F ext − F d ) is added to increase the ability of tracking effort. In this section, an accidental contact-loss due to the step-change of environment is considered as a fault. Hence, a fault-detection-based virtual energy tank is suggested to decouple force regulation as prior safety criterion. x t is defined as a state associated with the tank that stores energy. The dynamics of the virtual energy tank variable is expressed by [34]: where u t is the control input of the tank, and u f i is the additive control input signal of the force tracking impedance control scheme. u t and u f i are calculated by: Three switching parameters α, β, γ are defined as the following: where T t is energy stored in the virtual tank. α is a switched parameter when the energy tank reaches its lower bound T l . β is a switched parameter in case the energy tank from the force/impedance controller exceeds an upper bound T u . Then, this parameter is used to detect an occurrence of contact-loss. γ is presented to guarantee the stability of the overall system.

Remark 9.
The value of the lower bound T l is set to be greater than 0 to avoid the singularity. In this paper, the lower bound is set at 0.1, and the parameter α is consequently equal to 1.
The energy stored in the virtual tank T t is calculated as Finally, extended dynamics of the impedance with contact-loss consideration is rewritten as the following: where u imp−cl is the modified control input signal of impedance/contact-loss as As can be seen in Equation (40), the additive term u f i and switching parameters β and γ are supplemented to the force tracking control to decouple the signal from the force regulation when contact-loss happens. To demonstrate the stability of the closed-loop control system, the following candidate Lyapunov function is introduced [50].
Substituting Equation (34) into Equation (44) results in: As introduced in Equation (35), to avoid the singularity, the lower limitation T l should be set to be greater than zero and then, the switching α should be kept at 1 at all times. Therefore, the Equation (24) is equivalent to: Combining with the condition (37), then Equation (46) becomes: Thus, the closed-loop control system with energy tank is passive stable.

Modification of Trajectory
The previous section presented the effectiveness and stability of using the energy tank for contact-loss detection. However, this action only helps to decouple the signal command from force control from the overall controller scheme. Thereby, the end-effector follows the setup desired trajectory. Moreover, due to reaching the upper bound of the energy tank, the end-effector may attempt to regulate force until draining all energy in the virtual tank in some cases. Consequently, unsafe motions may be executed. To completely exhibit safety motion when contact-loss happens, a new trajectory should be designed to smooth the motion of the end-effector. For this purpose, a controlled shaping function ϕ(ζ), where ζ is a shift variable, is introduced as [36] is a shifted variable, ∆p is the error between the modified trajectory and actual position, ∆x si is an arbitrary constant that describes a robust region which denotes the distance between the end-effector to the set-point for shaping motion, and d imax is an arbitrary constant that denotes the width of the region for shaping. It follows that if the desired force F d and the position error ∆x encloses the robust region and passes the frontier, then the controller should be deactivated. The role of designing the shaping function is to handle unexpected motion when contact-loss happens. A shaping function ρ is introduced to generate a smooth transition between the force and impedance control; thus, providing a smooth transition in position performance of the end-effector.
After generating a smooth transition between force and impedance control when contact loss occurs, the next important step is designing a suitable set-point such that the robot can follow without shock or unexpected motion. In order to avoid fast and unexpected motions in case of contact loss, which are caused by the impedance controller for large deflection from the set point, the set-point of the impedance control, called as commanded trajectory x c , is set close to the contact point x contact as x c = x f ade = x contact − ε set , where x contact = x env is the measured position of the end-effector at t 0 and ε set is a small vector (indeed, when contacting with the environment, the measured position of the end-effector is the same as the environment. In the simulation, for illustration of impedance control, x contact is the actual position of the end-effector and this variable is different from the location of environment). A possible fading function is applied to obtain a continuous transition to the new set point. The flowchart for the fading trajectory is illustrated in Figure 4. From this illustration, the setpoint is determined based on the switching parameter β and comparison between the current and last position based time-delay than using predetermined or knowledge of the location of environment as in [35,36].

Configuration of Testing Environments
In this simulation, we verify the proposed methodology on a 6-DOF serial hydraulic KIRO robot. The robot designed is structured with the five first links driven by cylinders and the last link driven by the rotary actuator as shown in Figure 1. The parameters of the 6-DOF KIRO robot is described in Appendix A.
The simulation is constructed using MATLAB/Simulink program. The parameters of the hydraulic specification and hydraulic actuators are described in Tables 1 and 2, respectively. The controller parameters for contact-loss impedance control are described in Table 3. For the position control, the sliding surface is designed with K 1 = 100 × eye(6), λ 1 = 10 × diag (1, 1, 1, 1, 1, 2), and η 1 = 0.1 × eye (6). For the torque control, the controller gains are utilized with K 2 = 100 × eye (6), and η 2 = 0.1 × eye (6). Since the function "tanh(.)" in Equations (16) and (24) takes place of the use of "sign(.)" function, then the values of all σ ij are selected small. In this work, these parameters are set equal to 10 −3 . The sampling time for simulation is set as ts = 0.001s, and the simulation is implemented in the period of T period = 12(s).
Remark 10. In our simulation, the environment stiffness and impedance behavior matrix D m , B m , K m are supposed to be known and set with fixed constants. We aim to demonstrate the response of the system when applying the impedance contact-loss algorithm. Practically, the impedance dynamics directly relates to the design of impedance control gains. Therefore, adaptive gains can be further investigated in future work.

Remark 11.
The gains K P f , K I f are firstly initialized based on system characteristic and trial-and-error method. The technique for identifying the initial value and adaptive laws to adopt with working operation will be developed in future works.

Remark 12.
The value of upper bound T u is chosen based on system characteristics and designer knowledge. Different systems and working conditions may result in different values. Therefore, the choice of a suitable value for the upper bound will be discussed in future works.

Simulation Results
It is noteworthy in this scope that we intensively focus on investigations of system behaviors and stability under contact-loss. The position tracking performance between the conventional and the one with adaptive law algorithms under the presence of uncertainties and disturbances in free-space, and an enhanced force tracking performance can be referred to previous works in [9,10,45].
The first simulation is given in a case of still constrained framework. We suppose the end-effector of the 6-DOF-KIRO robot moves toward along the z-direction as shown in Figure 5 To mimic the contact-loss in two case studies, the environment profile is set as a step function as As can be seen in Equation (50), we suppose that the contact loss happens at the time of t = 7.62 s. Indeed, the contact-loss can suddenly occur if the position along the X-direction suddenly changes. Besides, the end-effector of the KIRO robot must apply a constant force of 100 N to the wall. The desired force induced from contacting command with the environment and moments induced from the command motions such as twisting, bending, etc., at the end-effector. Then, the desired force is set as The manipulator behavior in case of treating impedance-contact-loss based on the virtual energy tank concept is compared with that of conventional impedance control. We suppose that the endeffector is commanded to work along the X-direction. Position tracking performances of conventional impedance and impedance contact-loss in X, Y, and Z directions are displayed in Figures 6-8  To mimic the contact-loss in two case studies, the environment profile is set as a step function as As can be seen in Equation (50), we suppose that the contact loss happens at the time of t = 7.62 s. Indeed, the contact-loss can suddenly occur if the position along the X-direction suddenly changes. Besides, the end-effector of the KIRO robot must apply a constant force of 100 N to the wall. The desired force induced from contacting command with the environment and moments induced from the command motions such as twisting, bending, etc., at the end-effector. Then, the desired force is set as The manipulator behavior in case of treating impedance-contact-loss based on the virtual energy tank concept is compared with that of conventional impedance control. We suppose that the end-effector is commanded to work along the X-direction. Position tracking performances of conventional impedance and impedance contact-loss in X, Y, and Z directions are displayed in Figures 6-8, respectively. The continuous black line represents desired trajectory (•) desired , the black dot-line represents the arbitrary environment (•) env , the blue dot-dash line represents the end-effector position of conventional impedance response (•) imp , and the dash red line represents the end-effector trajectory of impedance contact-loss (•) imp−cl . For the shake of simplicity in the case of force responses, the desired force (F desired ), external impedance force (F ext−imp ), and external impedance contact-loss (F ext−imp−cl ) are presented by italic symbols hereafter since these values are investigated in only X-direction.  Position(mm) X desired X imp X imp-cl X env  At the time of approximately 7.62 s, the sudden change in the environment framework causes contact-loss. At that time, the interacting force immediately decreases to zero, and the tank variable increases until its energy reaches the upper limit u T at the time of approximately 7.63 s. Then, the parameters β and φ simultaneously switch their values as shown in Figure 9. Due to using the energy tank, the end-effector can generally maintain the same level with a very small position shifted. This is because the force control signal is deactivated, and the energy stored in a virtual spring-damper of each joint is dissipated. In the case of conventional impedance control, the end-effector moves until it reaches the constrained framework again at the time of 8.2 s as the requirement of force tracking effort. The force performances in two situations are depicted in Figure 10    At the time of approximately 7.62 s, the sudden change in the environment framework causes contact-loss. At that time, the interacting force immediately decreases to zero, and the tank variable increases until its energy reaches the upper limit u T at the time of approximately 7.63 s. Then, the parameters β and φ simultaneously switch their values as shown in Figure 9. Due to using the energy tank, the end-effector can generally maintain the same level with a very small position shifted. This is because the force control signal is deactivated, and the energy stored in a virtual spring-damper of each joint is dissipated. In the case of conventional impedance control, the end-effector moves until it reaches the constrained framework again at the time of 8.2 s as the requirement of force tracking effort. The force performances in two situations are depicted in Figure 10 in which the continuous black line represents desired force desired F , the blue dot-dash line represents external interactive force , and the dash red line represents the external interactive force of in the X-direction. At the time of approximately 7.62 s, the sudden change in the environment framework causes contact-loss. At that time, the interacting force immediately decreases to zero, and the tank variable increases until its energy reaches the upper limit T u at the time of approximately 7.63 s. Then, the parameters β and ϕ simultaneously switch their values as shown in Figure 9. Due to using the energy tank, the end-effector can generally maintain the same level with a very small position shifted. This is because the force control signal is deactivated, and the energy stored in a virtual spring-damper of each joint is dissipated. In the case of conventional impedance control, the end-effector moves until it reaches the constrained framework again at the time of 8.2 s as the requirement of force tracking effort. The force performances in two situations are depicted in Figure 10 in which the continuous black line represents desired force F desired , the blue dot-dash line represents external interactive force in conventional impedance F ext−imp , and the dash red line represents the external interactive force of impedance-contact-loss F ext−imp−cl in the X-direction. tank, the end-effector can generally maintain the same level with a very small position shifted. This is because the force control signal is deactivated, and the energy stored in a virtual spring-damper of each joint is dissipated. In the case of conventional impedance control, the end-effector moves until it reaches the constrained framework again at the time of 8.2 s as the requirement of force tracking effort. The force performances in two situations are depicted in Figure 10 in which the continuous black line represents desired force desired F , the blue dot-dash line represents external interactive force , and the dash red line represents the external interactive force of impedance-contact-loss in the X-direction.  The parameter β changes its value from 1 to 0 after contact-loss happens along with the change in parameter φ and tank variable t x in Figure 9. The time for switching value of both parameter β is slightly delayed from the time of contact-loss. The reason for this phenomenon is that when contactloss happens at the time of 7.62 s, the value of the tank variable t x increases; thus, increasing the value of energy tank T until reaching the upper limitation u T . When the value of T equals u T , parameter β is automatically passively switched from 1 to 0, and the trajectory is then faded. The switching of parameter φ occurs when satisfying both two conditions: after contact-loss is detected and condition in Equation (27). It is noticed that the problem of overshoot at the beginning of contact motion is not considered in this study. This problem will be considered as the next step of improvements for safety interactions.
In the second simulation, we consider the case of a rippled surface of the constrained framework to evaluate responses of the energy tank variable under several different values of force control gains. The constrained framework is presented as Similar to the first case study, we suppose that the contact-loss happens at the time of 7.62 s. In this case study, only tracking evaluation along the X-direction and response of the virtual energy tank variable are discussed due to similar responses of tracking performances in Y-and Z-directions and other switching parameters. Furthermore, we investigate the end-effector and energy variable behaviors under several different control gains. Herein, we consider nine different control gains and we display them in different line styles. The tracking effort and response of the energy variable are displayed in Figures 11 and 12. The parameter β changes its value from 1 to 0 after contact-loss happens along with the change in parameter ϕ and tank variable x t in Figure 9. The time for switching value of both parameter β is slightly delayed from the time of contact-loss. The reason for this phenomenon is that when contact-loss happens at the time of 7.62 s, the value of the tank variable x t increases; thus, increasing the value of energy tank T until reaching the upper limitation T u . When the value of T equals T u , parameter β is automatically passively switched from 1 to 0, and the trajectory is then faded. The switching of parameter ϕ occurs when satisfying both two conditions: after contact-loss is detected and condition in Equation (27). It is noticed that the problem of overshoot at the beginning of contact motion is not considered in this study. This problem will be considered as the next step of improvements for safety interactions.
In the second simulation, we consider the case of a rippled surface of the constrained framework to evaluate responses of the energy tank variable under several different values of force control gains. The constrained framework is presented as Similar to the first case study, we suppose that the contact-loss happens at the time of 7.62 s. In this case study, only tracking evaluation along the X-direction and response of the virtual energy tank variable are discussed due to similar responses of tracking performances in Y-and Z-directions and other switching parameters. Furthermore, we investigate the end-effector and energy variable behaviors under several different control gains. Herein, we consider nine different control gains and we display them in different line styles. The tracking effort and response of the energy variable are displayed in Figures 11 and 12. Appl. Sci. 2019, 9, x FOR PEER REVIEW 17 of 23 Figure 11. Response of the end-effector in the X-direction under different force control gains. Particularly, green lines, blue lines, and red lines denote tracking efforts when the force proportional gain is set as , respectively. Among that, dashed lines, dasheddot lines, and continuous lines present system behaviors when the force integral gains are , respectively. (Note that since the force control is only regulated in the Xdirection, then in this case, we consider force control gains as a scalar, not matrix as defined above.) As can be seen in Figures 11 and 12, the tracking performance of the end-effector changes under rippled framework no matter how different values of force control gains are. Due to consideration of the rippled surface, the force tracking effort varies around the desired value of 100 (N). The force tracking effort directly affects position tracking performance. Generally, larger values of the force control gains generate more accurate, faster, and smoother response of the end-effector when contactloss happens. However, too-large values of control gains (in this case is when the value of Pf K is greater than 2) may degrade tracking efforts and execute unexpected motions as a trade-off. As revealed in Figure 12, large gains may result in more overshoot phenomena at the beginning. Additionally, from Figure 12, since all lines reach the upper bound at the same time, the switching parameter β changes its state at the time of approximately 7.63 s, the same as the first case study despite different force control gains.
In the third case study, we consider friction exerting in the Z-direction. This friction force is supposed to be proportional to external force exerting at the end-effector in the X-direction as  Particularly, green lines, blue lines, and red lines denote tracking efforts when the force proportional gain is set as , respectively. Among that, dashed lines, dasheddot lines, and continuous lines present system behaviors when the force integral gains are , respectively. (Note that since the force control is only regulated in the Xdirection, then in this case, we consider force control gains as a scalar, not matrix as defined above.) As can be seen in Figures 11 and 12, the tracking performance of the end-effector changes under rippled framework no matter how different values of force control gains are. Due to consideration of the rippled surface, the force tracking effort varies around the desired value of 100 (N). The force tracking effort directly affects position tracking performance. Generally, larger values of the force control gains generate more accurate, faster, and smoother response of the end-effector when contactloss happens. However, too-large values of control gains (in this case is when the value of Pf K is greater than 2) may degrade tracking efforts and execute unexpected motions as a trade-off. As revealed in Figure 12, large gains may result in more overshoot phenomena at the beginning. Additionally, from Figure 12, since all lines reach the upper bound at the same time, the switching parameter β changes its state at the time of approximately 7.63 s, the same as the first case study despite different force control gains.
In the third case study, we consider friction exerting in the Z-direction. This friction force is supposed to be proportional to external force exerting at the end-effector in the X-direction as Particularly, green lines, blue lines, and red lines denote tracking efforts when the force proportional gain is set as K P f = 1, K P f = 1.75, K P f = 2, respectively. Among that, dashed lines, dashed-dot lines, and continuous lines present system behaviors when the force integral gains are K I f = 1, K I f = 2, K I f = 3, respectively. (Note that since the force control is only regulated in the X-direction, then in this case, we consider force control gains as a scalar, not matrix as defined above.) As can be seen in Figures 11 and 12, the tracking performance of the end-effector changes under rippled framework no matter how different values of force control gains are. Due to consideration of the rippled surface, the force tracking effort varies around the desired value of 100 (N). The force tracking effort directly affects position tracking performance. Generally, larger values of the force control gains generate more accurate, faster, and smoother response of the end-effector when contact-loss happens. However, too-large values of control gains (in this case is when the value of K P f is greater than 2) may degrade tracking efforts and execute unexpected motions as a trade-off. As revealed in Figure 12, large gains may result in more overshoot phenomena at the beginning. Additionally, from Figure 12, since all lines reach the upper bound at the same time, the switching parameter β changes its state at the time of approximately 7.63 s, the same as the first case study despite different force control gains.
In the third case study, we consider friction exerting in the Z-direction. This friction force is supposed to be proportional to external force exerting at the end-effector in the X-direction as where µ = 0.001 is ratio parameters. Similar to the second case study, green lines, blue lines, and red lines denote tracking efforts when the force proportional gain is set as K P f = 1, K P f = 1.75, K P f = 2, respectively. Among that, dashed line, dashed-dot lines, and continuous lines present the system behaviors when the force integral gains are K I f = 1, K I f = 2, K I f = 3, respectively.
As can be seen in Figure 13, at the time of contact-loss, different values of force control gains result in different responses of the end-effector. Among the two control gains, the proportional gain K P f dominates the tracking performance. Small gains result in slow response but less overshoot at the beginning. On the contrary, large values of proportional gain generally lead to fast response and smoother movement when contact-loss happens. The integral gains can help to enhance tracking performance; however, these gains cannot be arbitrarily set to large values. Figure 14 shows that responses of the energy variable also vary under the rippled framework and different values of the force control gains. Thereby, one can easily deduce that different times for switching the parameter β are obtained due to different times for reaching the upper bound of the energy variable as shown in Figure 15 (the parameter β switches if and only if the energy variable reaches the upper bound). Therefore, the control gains should be carefully designed in the case of friction as one considerable factor to achieve better performance. As can be seen in Figure 13, at the time of contact-loss, different values of force control gains result in different responses of the end-effector. Among the two control gains, the proportional gain Pf K dominates the tracking performance. Small gains result in slow response but less overshoot at the beginning. On the contrary, large values of proportional gain generally lead to fast response and smoother movement when contact-loss happens. The integral gains can help to enhance tracking performance; however, these gains cannot be arbitrarily set to large values. Figure 14 shows that responses of the energy variable also vary under the rippled framework and different values of the force control gains. Thereby, one can easily deduce that different times for switching the parameter β are obtained due to different times for reaching the upper bound of the energy variable as shown in Figure 15 (the parameter β switches if and only if the energy variable reaches the upper bound). Therefore, the control gains should be carefully designed in the case of friction as one considerable factor to achieve better performance.  Position(mm)  As can be seen in Figure 13, at the time of contact-loss, different values of force control gains result in different responses of the end-effector. Among the two control gains, the proportional gain Pf K dominates the tracking performance. Small gains result in slow response but less overshoot at the beginning. On the contrary, large values of proportional gain generally lead to fast response and smoother movement when contact-loss happens. The integral gains can help to enhance tracking performance; however, these gains cannot be arbitrarily set to large values. Figure 14 shows that responses of the energy variable also vary under the rippled framework and different values of the force control gains. Thereby, one can easily deduce that different times for switching the parameter β are obtained due to different times for reaching the upper bound of the energy variable as shown in Figure 15 (the parameter β switches if and only if the energy variable reaches the upper bound). Therefore, the control gains should be carefully designed in the case of friction as one considerable factor to achieve better performance.  Position(mm) Furthermore, it is noteworthy that at the beginning, when the motion of the manipulator is switched from free-motion to constrained-motion, the free motion of the end-effector is also considered as contact-loss. If the initial point is too far from the constrained framework, more energy is required to drive the end-effector to the target and, thus, the energy may reach the limitation before the system is in constrained-framework and execute undesirable motions. Thereby, the initial point should be set as not far from the contact point. This problem can also be solved by increasing a large value of the upper bound u T ; however, smooth movement may not be executed because more time is required for the energy variable to reach the upper bound. Hence, one should refer step by step to solve this problem by designing adaptive laws to decouple the force control in free motion, or adaptive upper boundary to satisfy different working conditions.

Conclusions
In this paper, we intensively presented the method for solving the problem of contact-loss based on using the energy tank for the n-degree-of-freedom serial manipulator. The contact-loss probably occurs in the case of constrained motion when the end-effector of the manipulator moves out of the constrained framework. In this scenario, the end-effector follows the virtual commanded trajectory resulting from the force tracking effort. When the contact-loss happens, the potential energy stored in the virtual spring/damper mechanism in each joint has instantaneously released. According to this behavior, the fault detection based on the energy tank was designed with several passive switched parameters to decouple the force action. This helped to prevent the system from getting damaged or unexpected motions. Furthermore, to avoid shock at the point of contact-loss, the command trajectory was modified based on time-delay instead of using the location of environment as presented in previous works. This method compared the end-effector positions between the recent step and previous step along with passive switched parameters to smooth the output performance. Additionally, the adaptive laws for the inner position control loop were derived to enhance system performance and facilitate further development of advanced methodology such as force-sensorless using observer or optimization. From the simulations, the end-effector movement after contact-loss happened was maintained moving in same level as in constrained framework and the smooth movement was executed thanks to utilizing the fading trajectory algorithm to generate the new commanded trajectory in comparison with those under conventional impedance control. Consequently, the control methodology showed its ability in both tracking performance and contactloss reaction. Additionally, to comprehensively investigate the force-controlled system behaviors, we discussed the influence of force regulation loop to system response and the end-effector dynamics in which different values of force control gains were examined. Generally, gains of the force control loop, and also boundary parameters of the energy tank, directly decided the smooth motion and Furthermore, it is noteworthy that at the beginning, when the motion of the manipulator is switched from free-motion to constrained-motion, the free motion of the end-effector is also considered as contact-loss. If the initial point is too far from the constrained framework, more energy is required to drive the end-effector to the target and, thus, the energy may reach the limitation before the system is in constrained-framework and execute undesirable motions. Thereby, the initial point should be set as not far from the contact point. This problem can also be solved by increasing a large value of the upper bound T u ; however, smooth movement may not be executed because more time is required for the energy variable to reach the upper bound. Hence, one should refer step by step to solve this problem by designing adaptive laws to decouple the force control in free motion, or adaptive upper boundary to satisfy different working conditions.

Conclusions
In this paper, we intensively presented the method for solving the problem of contact-loss based on using the energy tank for the n-degree-of-freedom serial manipulator. The contact-loss probably occurs in the case of constrained motion when the end-effector of the manipulator moves out of the constrained framework. In this scenario, the end-effector follows the virtual commanded trajectory resulting from the force tracking effort. When the contact-loss happens, the potential energy stored in the virtual spring/damper mechanism in each joint has instantaneously released. According to this behavior, the fault detection based on the energy tank was designed with several passive switched parameters to decouple the force action. This helped to prevent the system from getting damaged or unexpected motions. Furthermore, to avoid shock at the point of contact-loss, the command trajectory was modified based on time-delay instead of using the location of environment as presented in previous works. This method compared the end-effector positions between the recent step and previous step along with passive switched parameters to smooth the output performance. Additionally, the adaptive laws for the inner position control loop were derived to enhance system performance and facilitate further development of advanced methodology such as force-sensorless using observer or optimization. From the simulations, the end-effector movement after contact-loss happened was maintained moving in same level as in constrained framework and the smooth movement was executed thanks to utilizing the fading trajectory algorithm to generate the new commanded trajectory in comparison with those under conventional impedance control. Consequently, the control methodology showed its ability in both tracking performance and contact-loss reaction. Additionally, to comprehensively investigate the force-controlled system behaviors, we discussed the influence of force regulation loop to system response and the end-effector dynamics in which different values of force control gains were examined. Generally, gains of the force control loop, and also boundary parameters of the energy tank, directly decided the smooth motion and stability of the whole system in not only free motion but also constrained framework, and after contact-loss happened. Therefore, depending on working characteristic and requirements, the boundary limitation of the energy tank and controller gains should be dedicatedly designed to obtain optimal performance. Besides, remaining problems such as influences of friction, initial condition, or adaptive boundary evoke interesting topics for further study and development. This algorithm can be expanded to other practical applications for safety such as at the end of mechanical turning, milling, or planning process, manufacture in automation fields, or reaction of the end-effector when suddenly dropping a payload.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
According to [44], Denavit-Hartenberg parameters of the 6-DOF hydraulic manipulator are specified as shown in Table A1.
with c i = cos(q i ), s i = sin(q i ), c ij = c i c j − s i s j , s ij = c i s j + s i c j , (i = 1, 2, . . . , 6; j = 1, 2, . . . , 6; i j). The Jacobian matrix and other matrices of inertia, Coriolis, and Gravity in Cartesian space are calculated from the MATLAB program based on the those in joint space that follow the procedure in [44,51].