Parameters Sensitivity Analysis of Position-Based Impedance Control for Bionic Legged Robots ’ HDU

Kaixian Ba 1, Bin Yu 1,2,*, Zhengjie Gao 1, Wenfeng Li 1, Guoliang Ma 1 and Xiangdong Kong 1,2,3 1 School of Mechanical Engineering, Yanshan University, Qinhuangdao 066004, China; bkx@ysu.edu.cn (K.B.); gzj@stumail.ysu.edu.cn (Z.G.); xrlzlwf@stumail.ysu.edu.cn (W.L.); mgl@stumail.ysu.edu.cn (G.M.); xdkong@ysu.edu.cn (X.K.) 2 Hebei Provincial Key Laboratory of Heavy Machinery Fluid Power Transmission and Control, Qinhuangdao, China 3 National Engineering Research Center for Local Joint of Advanced Manufacturing Technology and Equipment, Yanshan University, Qinhuangdao 066004, China * Correspondence: yb@ysu.edu.cn; Tel.: +86-0335-807-4618


Introduction
Bionic legged robots are better at adapting to unknown and unstructured environments.Their unique advantages, such as overcoming obstacles and executing tasks in the wild, have made them a major focus of research in the robotic domain [1][2][3][4].For the hydraulic drive robot, the highly integrated valve-controlled cylinder composes the drive component, which is called the hydraulic drive unit (HDU) [5,6].During the robotic motion process, the robotic feet interact with the ground frequently.This means that the demand for HDU not only includes characteristics of response ability and high control accuracy, but also dynamic compliance.Thus, the impact on the hydraulic system can be obviously reduced, which can protect the mechanical structure and components, and improve the moving stability of the robots.
The selection of control methods directly affects the compliance of the HDU.As a commonly used control method for compliance, the impedance control method has been widely applied to motor-driving legged robots such as the Tekken [7], Scout [8] and MIT cheetah robot [9].In recent years, as the hydraulic-driven legged robot became the focus of increased research, impedance control was also applied to this kind of robot, as exemplified by robots such as Bigdog [10], HyQ [11], and Scalf-1 [12].Force-based and position-based impedance control methods are often used for dynamic compliance control.Their basic principle can be expressed as follows.Firstly, the hydraulic control system is taken as the control inner loop, and a dynamic control outer loop is attached to the system.When an external disturbance acts on the system, the input signal of the control inner

Introduction of the HDU and Its Performance Test Platform
As the driver of the leg joint on bionic legged robots, the HDU is a highly integrated system of servo valve-controlled symmetrical cylinder.The author's institute participates in the design of the hydraulic quadruped robot.The quadruped robot prototype, single leg, and HDU performance test platform are shown in Figure 1a-c, respectively.The performance of the HDU directly affects the performance of the whole robot.Thus, a special performance test platform is built to study the methods for HDU.The schematic of the test platform is shown in Figure 2. The principle of the electro-hydraulic load simulator in Figure 2 is widely used in many fields such as aviation, aerospace, and vessel and construction machining [5].In Figure 2, the left part is a HDU-adopted position closed-loop control that contains a small servo valve, servo cylinder, position sensor, and force sensor.The right part is another HDU-adopted force closed-loop control that contains the same type of servo valve and servo cylinder.Two parts' cylinder rods are jointed rigidly by the thread of a force sensor.The HDU performance test platform is showed in Figure 3a.The controller adopted is dSPACE, a semi-physical simulation platform, which is showed in Figure 3b.The performance of the HDU directly affects the performance of the whole robot.Thus, a special performance test platform is built to study the methods for HDU.The schematic of the test platform is shown in Figure 2. The performance of the HDU directly affects the performance of the whole robot.Thus, a special performance test platform is built to study the methods for HDU.The schematic of the test platform is shown in Figure 2. The principle of the electro-hydraulic load simulator in Figure 2 is widely used in many fields such as aviation, aerospace, and vessel and construction machining [5].In Figure 2, the left part is a HDU-adopted position closed-loop control that contains a small servo valve, servo cylinder, position sensor, and force sensor.The right part is another HDU-adopted force closed-loop control that contains the same type of servo valve and servo cylinder.Two parts' cylinder rods are jointed rigidly by the thread of a force sensor.The HDU performance test platform is showed in Figure 3a.The controller adopted is dSPACE, a semi-physical simulation platform, which is showed in Figure 3b.The principle of the electro-hydraulic load simulator in Figure 2 is widely used in many fields such as aviation, aerospace, and vessel and construction machining [5].In Figure 2, the left part is a HDU-adopted position closed-loop control that contains a small servo valve, servo cylinder, position sensor, and force sensor.The right part is another HDU-adopted force closed-loop control that contains the same type of servo valve and servo cylinder.Two parts' cylinder rods are jointed rigidly by the thread of a force sensor.The HDU performance test platform is showed in Figure 3a.The controller adopted is dSPACE, a semi-physical simulation platform, which is showed in Figure 3b.The performance of the HDU directly affects the performance of the whole robot.Thus, a special performance test platform is built to study the methods for HDU.The schematic of the test platform is shown in Figure 2. The principle of the electro-hydraulic load simulator in Figure 2 is widely used in many fields such as aviation, aerospace, and vessel and construction machining [5].In Figure 2, the left part is a HDU-adopted position closed-loop control that contains a small servo valve, servo cylinder, position sensor, and force sensor.The right part is another HDU-adopted force closed-loop control that contains the same type of servo valve and servo cylinder.Two parts' cylinder rods are jointed rigidly by the thread of a force sensor.The HDU performance test platform is showed in Figure 3a.The controller adopted is dSPACE, a semi-physical simulation platform, which is showed in Figure 3b.

Mathematic Model of HDU Position-Based Impedance Control
Impedance control is one type of active compliance control.In particular, it refers to an active compliance control that applies the system equivalent to the second-order mass-spring-damping system.By adopting impedance control, a system can be equipped with the dynamic compliance of a second-order mass-spring-damping system when a disturbance force is applied to the system.
The impedance control is composed by an impedance control inner loop and impedance control outer loop.The impedance control inner loop refers to the closed-loop control, which is realized in the inner loop of the hydraulic position closed-loop control system during the impedance control.In the position-based impedance control of this paper, the impedance control inner loop refers to the position closed-loop control.The impedance control outer loop refers to the open-loop control where the external disturbance signal is transferred to the input signal of impedance control inner loop during the impedance control.

Principles of Impedance Control
The inner loop of the position-based impedance control is a closed-loop control.When the inner loop is affected by a disturbance force, the impedance control outer loop should be added to the system to equip the system with impedance control characteristics.The function of this outer loop is to transform the disturbance force into position error.Then, the desired stiffness K D can be obtained, which causes an elastic force.In the same way, desired damping C D and desired mass m D are obtained, which can cause viscous force and inertia force, respectively.The HDU force schematics with impedance control outer loop are shown respectively in Figures 4 and 5.
During the robot's walking process, the load, such as grounds and steps, provides the disturbance force to the HDU, because the force sensor is mounted on the piston end of the HDU.In this paper, the force control system of the simulated load provides the disturbance force to the performance test platform.

Mathematic Model of HDU Position-Based Impedance Control
Impedance control is one type of active compliance control.In particular, it refers to an active compliance control that applies the system equivalent to the second-order mass-spring-damping system.By adopting impedance control, a system can be equipped with the dynamic compliance of a second-order mass-spring-damping system when a disturbance force is applied to the system.
The impedance control is composed by an impedance control inner loop and impedance control outer loop.The impedance control inner loop refers to the closed-loop control, which is realized in the inner loop of the hydraulic position closed-loop control system during the impedance control.In the position-based impedance control of this paper, the impedance control inner loop refers to the position closed-loop control.The impedance control outer loop refers to the open-loop control where the external disturbance signal is transferred to the input signal of impedance control inner loop during the impedance control.

Principles of Impedance Control
The inner loop of the position-based impedance control is a closed-loop control.When the inner loop is affected by a disturbance force, the impedance control outer loop should be added to the system to equip the system with impedance control characteristics.The function of this outer loop is to transform the disturbance force into position error.Then, the desired stiffness D K can be obtained, which causes an elastic force.In the same way, desired damping D C and desired mass D m are obtained, which can cause viscous force and inertia force, respectively.The HDU force schematics with impedance control outer loop are shown respectively in Figures 4 and 5.
During the robot's walking process, the load, such as grounds and steps, provides the disturbance force to the HDU, because the force sensor is mounted on the piston end of the HDU.In this paper, the force control system of the simulated load provides the disturbance force to the performance test platform.As it is shown in Figure 4, the load is pressed to position X Δ , 2 p B refers to damping coefficient at load, 2 t m refers to equivalent mass at load, and 2 f F refers to friction at load.The force acting on the piston, which is provided by the sensor, is the disturbance force of the HDU position control system, and is defined as Δ a F .Further, the force acting on the load, which is provided by the sensor, is defined as Δ b F .The force balance relation in Figure 4 can be expressed as follows: Third, the force schematic of the HDU is showed in Figure 5.As it is shown in Figure 4, the load is pressed to position ∆X, B p2 refers to damping coefficient at load, m t2 refers to equivalent mass at load, and F f 2 refers to friction at load.The force acting on the piston, which is provided by the sensor, is the disturbance force of the HDU position control system, and is defined as ∆F a .Further, the force acting on the load, which is provided by the sensor, is defined as ∆F b .The force balance relation in Figure 4 can be expressed as follows: Third, the force schematic of the HDU is showed in Figure 5.   F refers to friction in the cylinder.The position error can be expressed as follows: Based on the theoretical analysis above, the schematic of the position-based impedance control is shown in Figure 6.As can be seen in Figure 6, when the disturbance force tested by the force sensor acts on the HDU, the impedance control outer loop generates a corresponding calculated position that disturbs the input position.Then, the new input to the position control system is formed.Thus, the final input signal enters the position control inner loop, and a new output position is formed to equip the system with impedance characteristics.

The Block Diagram and State Space Presentation of Impedance Control System
The block diagram of HDU position-based impedance control is shown in Figure 7, where the detailed deduction and performance analysis of the inner loop is presented in previous research [5,15].In Figure 7, t m is conversion mass (including the piston, the displacement sensor, the force sensor, the connecting pipe, and the oil in the servo cylinder), r X is input position, x K is position sensor gain, PID K is proportion-integration-differentiation (PID) controller gain including proportional gain P K , integral gain I K and differential gain D K , K is load stiffness, p B is load and 2 Q is outlet oil flow.Where B p1 refers to the viscous damping coefficient, m t1 refers to the equivalent mass at piston, and F f 1 refers to friction in the cylinder.The position error can be expressed as follows: Based on the theoretical analysis above, the schematic of the position-based impedance control is shown in Figure 6.F refers to friction in the cylinder.The position error can be expressed as follows: Based on the theoretical analysis above, the schematic of the position-based impedance control is shown in Figure 6.As can be seen in Figure 6, when the disturbance force tested by the force sensor acts on the HDU, the impedance control outer loop generates a corresponding calculated position that disturbs the input position.Then, the new input to the position control system is formed.Thus, the final input signal enters the position control inner loop, and a new output position is formed to equip the system with impedance characteristics.

The Block Diagram and State Space Presentation of Impedance Control System
The block diagram of HDU position-based impedance control is shown in Figure 7, where the detailed deduction and performance analysis of the inner loop is presented in previous research [5,15].In Figure 7, t m is conversion mass (including the piston, the displacement sensor, the force sensor, the connecting pipe, and the oil in the servo cylinder), r X is input position, x K is position sensor gain, PID K is proportion-integration-differentiation (PID) controller gain including and 2 Q is outlet oil flow.As can be seen in Figure 6, when the disturbance force tested by the force sensor acts on the HDU, the impedance control outer loop generates a corresponding calculated position that disturbs the input position.Then, the new input to the position control system is formed.Thus, the final input signal enters the position control inner loop, and a new output position is formed to equip the system with impedance characteristics.

The Block Diagram and State Space Presentation of Impedance Control System
The block diagram of HDU position-based impedance control is shown in Figure 7, where the detailed deduction and performance analysis of the inner loop is presented in previous research [5,15].
In Figure 7, m t is conversion mass (including the piston, the displacement sensor, the force sensor, the connecting pipe, and the oil in the servo cylinder), K x is input position, is position sensor gain, K PID is proportion-integration-differentiation (PID) controller gain including proportional gain K P , integral gain K I and differential gain  Denote, where servo K is the transfer function of the servo valve, ω is the natural frequency of the servo valve, ζ is the damping ratio of the servo valve, and axv K is the servo valve gain.
( ) ( ) ( ) where, 1 K and 2 K express the transfer function of nonlinear pressure flow, defined as conversion coefficient in this paper), d C is the orifice flow coefficient of the spool valve, W is the area gradient of the spool valve, ρ is the density of hydraulic oil, s p is the system supply's oil pressure, 1 p is the left cavity pressure of the servo cylinder, 2 p is the right cavity pressure of the servo cylinder, and 0 p is the system return oil pressure.
( ) where, 3 K and 4 K express the transfer function of flow continuity, L represents the total piston stroke of the servo cylinder, 0 L is the initial piston position of the servo cylinder, ip C is the internal leakage coefficient of the servo cylinder, ep C is the external leakage coefficient of the servo cylinder, p A is the effective piston area of the servo cylinder, and e β is the effective bulk modulus.
The state variables in Figure 7 are expressed as follows: where the input variables are expressed as follows: Denote, where K servo is the transfer function of the servo valve, ω is the natural frequency of the servo valve, ζ is the damping ratio of the servo valve, and K axv is the servo valve gain.
where, K 1 and K 2 express the transfer function of nonlinear pressure flow, K d = C d W 2/ρ (K d is defined as conversion coefficient in this paper), C d is the orifice flow coefficient of the spool valve, W is the area gradient of the spool valve, ρ is the density of hydraulic oil, p s is the system supply's oil pressure, p 1 is the left cavity pressure of the servo cylinder, p 2 is the right cavity pressure of the servo cylinder, and p 0 is the system return oil pressure.
where, K 3 and K 4 express the transfer function of flow continuity, L represents the total piston stroke of the servo cylinder, L 0 is the initial piston position of the servo cylinder, C ip is the internal leakage coefficient of the servo cylinder, C ep is the external leakage coefficient of the servo cylinder, A p is the effective piston area of the servo cylinder, and β e is the effective bulk modulus.The state variables in Figure 7 are expressed as follows: where the input variables are expressed as follows: Disturbance variables are expressed as follows: The state space of the system can be expressed as follows: .
The physical meanings and initial values of the parameters in the control system block diagram are shown in Table 1.

Experiment of Position-Based Impedance Control
As a typical signal, the sinusoidal response is able to evaluate the impedance control performance under the input of different frequencies and amplitudes.In this paper, a sinusoidal signal is adopted to analyze the system performance and the sensitivity of the main parameters.To study their variation patterns in different conditions, four working conditions are tested in this paper.The details of the working conditions are shown in Table 2.In a position-based impedance control system, the load acts as an external disturbance force that is simulated by the force-based control system.Thus, the desired position can be defined as the ratio of the actual force to the desired impedance characteristic Z D , i.e., the ratio of the value of the force sensor to Z D .The actual position is the output of the system to be tested, i.e., the value of the position sensor.The experimental and simulation curves of the sinusoidal response are shown in  in sequence of working conditions [17].In a position-based impedance control system, the load acts as an external disturbance force that is simulated by the force-based control system.Thus, the desired position can be defined as the ratio of the actual force to the desired impedance characteristic D Z , i.e., the ratio of the value of the force sensor to D Z .The actual position is the output of the system to be tested, i.e., the value of the position sensor.The experimental and simulation curves of the sinusoidal response are shown in  in sequence of working conditions [17].In a position-based impedance control system, the load acts as an external disturbance force that is simulated by the force-based control system.Thus, the desired position can be defined as the ratio of the actual force to the desired impedance characteristic D Z , i.e., the ratio of the value of the force sensor to D Z .The actual position is the output of the system to be tested, i.e., the value of the position sensor.The experimental and simulation curves of the sinusoidal response are shown in  in sequence of working conditions [17].The mean values of performance indexes in different conditions are shown in Table 3.As it can be seen in Table 3, the values of the two performance indexes are close in experiment and simulation, which indicates that the experimental curves fit the simulation curves well.As the position-based impedance control theory in Figure 8 shows, the actual position is greater than the desired one.Thus, the amplitude attenuation is a negative value.It increases with the increase of the disturbance force's amplitude, but has few relationships with the frequency.In contrast, the phase angle delay increases with the frequency, but has few relationships with the amplitude.The mean values of performance indexes in different conditions are shown in Table 3.As it can be seen in Table 3, the values of the two performance indexes are close in experiment and simulation, which indicates that the experimental curves fit the simulation curves well.As the position-based impedance control theory in Figure 8 shows, the actual position is greater than the desired one.Thus, the amplitude attenuation is a negative value.It increases with the increase of the disturbance force's amplitude, but has few relationships with the frequency.In contrast, the phase angle delay increases with the frequency, but has few relationships with the amplitude.

Contrast between First-Order and Second-Order Sensitivity Analysis
Desired stiffness K D and desired damping C D are control parameters of the impedance control outer loop.Proportional gain K P and integral gain K i are control parameters of the position control inner loop.They affect the impedance control performance in different ways.So, in this paper, the system output position is mainly discussed, which is influenced by the variation of the four parameters.Because of space limitations, only one working condition (bias: 1500 N, amplitude: 1000 N, frequency: 2 Hz) is studied.The first-order and second-order trajectory sensitivity analysis methods researched in our previous works [15,16] are adopted to analyze the position variation when the four parameters increase 10% to 20%.The contrast curves are shown in Figures 12-15.
control inner loop.They affect the impedance control performance in different ways.So, in this paper, the system output position is mainly discussed, which is influenced by the variation of the four parameters.Because of space limitations, only one working condition (bias: 1500 N, amplitude: 1000 N, frequency: 2 Hz) is studied.The first-order and second-order trajectory sensitivity analysis methods researched in our previous works [15,16] are adopted to analyze the position variation when the four parameters increase 10% to 20%.The contrast curves are shown in Figures 12-15.outer loop.Proportional gain P K and integral gain i K are control parameters of the position control inner loop.They affect the impedance control performance in different ways.So, in this paper, the system output position is mainly discussed, which is influenced by the variation of the four parameters.Because of space limitations, only one working condition (bias: 1500 N, amplitude: 1000 N, frequency: 2 Hz) is studied.The first-order and second-order trajectory sensitivity analysis methods researched in our previous works [15,16] are adopted to analyze the position variation when the four parameters increase 10% to 20%.The contrast curves are shown in Figures 12-15.outer loop.Proportional gain P K and integral gain i K are control parameters of the position control inner loop.They affect the impedance control performance in different ways.So, in this paper, the system output position is mainly discussed, which is influenced by the variation of the four parameters.Because of space limitations, only one working condition (bias: 1500 N, amplitude: 1000 N, frequency: 2 Hz) is studied.The first-order and second-order trajectory sensitivity analysis methods researched in our previous works [15,16] are adopted to analyze the position variation when the four parameters increase 10% to 20%.The contrast curves are shown in Figures 12-15.The following conclusions can be reached from the four sets of curves.There is little difference between the first-order position variation and the second-order variation when the values of the four parameters increases 10%.Specifically, the curves of desired damping C D and integral gain i K almost overlap.In the other two sets, the maximum error is no more than 10% of the amplitude.When the values of the four parameters increase 20%, the variations of first-order and second-order are still The following conclusions can be reached from the four sets of curves.There is little difference between the first-order position variation and the second-order variation when the values of the four parameters increases 10%.Specifically, the curves of desired damping C D and integral gain K i almost overlap.In the other two sets, the maximum error is no more than 10% of the amplitude.When the values of the four parameters increase 20%, the variations of first-order and second-order are still close, although there is some increase of their error.Specifically, for desired damping C D and integral gain K i , the variation curves of first-order and second-order are still close.As for the proportional gain K P , its maximum error is no more than 10% of the amplitude.When it comes to the desired stiffness K D , its maximum error is no more than 20% of the amplitude.

Position variation Time t/s
The method of second-order trajectory sensitivity analysis has a very high accuracy, while its calculation is complicated and demands a lot of hard work.Particularly in this paper, when the four parameters increase less than 20%, the corresponding results of the first-order and second order sensitivity analysis method are close.So, in order to ensure the simplicity of calculation and application, the method of first-order sensitivity analysis is adopted to analyze the sensitivity of the four parameters under different working conditions.
The method studied previously can precisely analyze the sensitivity of the parameters.Compared with the second-order method, the calculation has been largely simplified under the method of first-order trajectory sensitivity analysis, but it also requires solving first-order linear non-homogeneous differential equations with time-varying factors, which makes the program complicated.So, in order to ensure high solving accuracy, a new method with an easier calculation is proposed.That is the first-order matrix sensitivity analysis.

Deduction of First-Order Matrix Sensitivity Analysis Theory
Combined with the position-based impedance control method, equation of the HDU system can be expressed as follows: g(x, u, α, t) = 0 (11) where x is m−dimensional state vector, u is r−dimensional vector unrelated to α, α is p−dimensional vector, and t is time.
The initial value of the state vector x 0 can be obtained by giving the initial value of the input vector u 0 and the initial value of parameter vector α 0 , and the initial state of the equations is: g(x 0 , u 0 , α 0 , t) = 0 (12) where the variation of parameter vector ∆α and input vector ∆u can change the value of ∆x, the error of the state variable x, which is expressed as follows: Expanded in the form of the first-order Taylor Series, Equation ( 9) can be expressed as follows: If we bring Equation (12) into Equation ( 14) and ignore the higher-order terms, then we can get: Equation ( 15) can also be expressed as follows: In Equation ( 17), supposing: where S u is m × r order matrix, and the n-th line indicates the relation among the n-th state variable x n .
In Equation ( 17), supposing: where, S α is the m × p order matrix, the n-th line indicates the relation among the n-th state variable and p parameter vectors.Bring Equations ( 17) and (18) into Equation ( 16), then: Equation ( 19) is an approximate expression of ∆x resulted from the change of parameter vector ∆α and input vector ∆u, in which S α indicates an m × p order parameter sensitivity matrix of parameter vector α with time-varying factors.S u indicates an m × r order input sensitivity matrix of input vector u with time-varying factors.
When taking no account of the variation of input vector, Equation ( 19) can be simplified as follows: The output equation of the system can be expressed as follows: where C and D are matrices of output equation factors.The change of output variable ∆Y resulting from the parameter variation can be reached after solving parameter sensitivity matrix S α .

Contrast between Two Analysis Method of First-Order Sensitivity
The state vectors' initial values of the servo-cylinder's position, velocity, and pressure of two chambers, the servo-valve's position, velocity and acceleration are zero.So, the initial value of parameter sensitivity matrix S α can be expressed as follows: Solve the sensitivity matrices in MATLAB, and then compare the results from the trajectory sensitivity in Section 3.1 and −S α , the inverse value of the control parameter's sensitivity matrix of the system output position.Due to space limitations, only one working condition (1500 N bias, 1000 N amplitude, 2 Hz frequency) is shown.The curves of parameter dynamic sensitivity in this situation is shown in Figure 16.It can be seen that the dynamic compliance curves of first-order matrix sensitivity deviate little from the curves of first-order trajectory sensitivity.Particularly, by comparing Figures 12-15 in Section 3.1, it can be found that the value calculated by first-order matrix sensitivity is more approximate to the result of second-order trajectory sensitivity and has higher precision than firstorder trajectory sensitivity.Moreover, only a two-dimension matrix calculation is needed, which avoids solving complicated differential equations with time-varying factors.So, in the research field of this paper, the first-order matrix sensitivity analysis method is more adapted than the first-order trajectory sensitivity analysis method.Whether calculating a high-order and multi-dimensional matrix is easier than solving differential equations with time-varying factors cannot be determined, so it requires further research to find which is better between the high-order matrix sensitivity analysis method and the high-order trajectory sensitivity analysis method.However, due to space constraints, this will not be discussed in this paper.

Contrast of Dynamic Sensitivity Analysis in Each Working Condition
For the convenience of contrast, the variation of system position response when each parameter increases 10% is calculated according to Equation (21).The curves of position variation with time is shown in Figure 17.It can be seen that the dynamic compliance curves of first-order matrix sensitivity deviate little from the curves of first-order trajectory sensitivity.Particularly, by comparing Figures 12-15 in Section 3.1, it can be found that the value calculated by first-order matrix sensitivity is more approximate to the result of second-order trajectory sensitivity and has higher precision than first-order trajectory sensitivity.Moreover, only a two-dimension matrix calculation is needed, which avoids solving complicated differential equations with time-varying factors.So, in the research field of this paper, the first-order matrix sensitivity analysis method is more adapted than the first-order trajectory sensitivity analysis method.Whether calculating a high-order and multi-dimensional matrix is easier than solving differential equations with time-varying factors cannot be determined, so it requires further research to find which is better between the high-order matrix sensitivity analysis method and the high-order trajectory sensitivity analysis method.However, due to space constraints, this will not be discussed in this paper.

Contrast of Dynamic Sensitivity Analysis in Each Working Condition
For the convenience of contrast, the variation of system position response when each parameter increases 10% is calculated according to Equation (21).The curves of position variation with time is shown in Figure 17.As it can be seen in Figure 17: 1.The variation of each parameter affects impedance control position output.The position varies periodically with the sinusoidal disturbance force.Among the parameters, the desired stiffness 2. The order of magnitudes of output position increases with the increase of sinusoidal disturbance force.However, there isn't an obvious relationship between the effects on output position and the frequency of the disturbance force.
Using the method of dynamic sensitivity analysis, the qualitative effects on impedance control performance is analyzed.In order to analyze the effects on main system performance indexes and the varying patterns of parameter sensitivity under different working conditions, a quantificational analysis is needed.

Sensitivity Indexes
Two measurement indexes are introduced to analyze the effects on main performance indexes resulting from the variation of parameters in different working conditions.As it can be seen in Figure 17: 1.
The variation of each parameter affects impedance control position output.The position varies periodically with the sinusoidal disturbance force.Among the parameters, the desired stiffness K D affects the output position much more than the others.The influence from integral gain K i is the most irrelevant.Desired damping C D and proportional gain K p have similar influences on the output position.With the disturbance force in sinusoidal variation, C D varies by following the curve of minus cosine, and K p varies by following the curve of cosine.

2.
The order of magnitudes of output position increases with the increase of sinusoidal disturbance force.However, there isn't an obvious relationship between the effects on output position and the frequency of the disturbance force.
Using the method of dynamic sensitivity analysis, the qualitative effects on impedance control performance is analyzed.In order to analyze the effects on main system performance indexes and the varying patterns of parameter sensitivity under different working conditions, a quantificational analysis is needed.

Sensitivity Indexes
Two measurement indexes are introduced to analyze the effects on main performance indexes resulting from the variation of parameters in different working conditions.
For sinusoidal response, in a stable sinusoidal period, the variation of parameters results in the change of the output position amplitude.The mean of amplitude attenuation is defined as the first sensitivity measurement index s 1 , which is expressed as follows: where, Similarly, for phase angle delay, another important index, its mean is defined as the second sensitivity measurement index s 2 , which is expressed as follows: where, Using the above two indexes, s 1 and s 2 , the effect on output position resulting from the variation of parameters can be quantificationally analyzed.

Sensitivity Histograms in Different Working Conditions
According to Equations ( 23)-( 27), the two sensitivity histograms are shown in Figures 18-21, when the four control parameters increase 10% in four working conditions.For sinusoidal response, in a stable sinusoidal period, the variation of parameters results in the change of the output position amplitude.The mean of amplitude attenuation is defined as the first sensitivity measurement index 1  s , which is expressed as follows: where, Similarly, for phase angle delay, another important index, its mean is defined as the second sensitivity measurement index 2 s , which is expressed as follows: where, Using the above two indexes, 1 s and 2 s , the effect on output position resulting from the variation of parameters can be quantificationally analyzed.

Sensitivity Histograms in Different Working Conditions
According to Equations ( 23)-( 27), the two sensitivity histograms are shown in Figures 18-21, when the four control parameters increase 10% in four working conditions.As it can be seen in the above Figures 18-21: 1.
In all conditions, S 1 is a positive value, so the increase of the four parameters will result in the reduction of the position amplitude attenuation.The S 1 of desired stiffness K D is much greater than the other parameters, which illustrates that the variation of K D has a remarkable influence on output position amplitude attenuation, and more relation with disturbance force amplitude than the disturbance force frequency.Under the first working condition, amplitude attenuation reduced about 0.2 mm when K D increased 10%.The S 1 of proportional gain K p is greater than that of C D and K i .The S 1 of K D and K p have some correlation with disturbance force amplitude, but have little correlation with disturbance force frequency.

2.
The S 2 of the four control parameters has positive or negative values in different conditions.An increase of the parameters has different influences on the phase angle delay of the output position.Specifically, having more relation with disturbance force frequency than the amplitude, the S 2 of K D and C D have nearly the same absolute but an opposite sign.Their increase has totally different influences on phase angle delay.An approximate 10% increase of them generates about 0.18 • of phase angle delay error.The S 2 of K p has correlation with both disturbance force amplitude and frequency.The larger the frequency, the more remarkable the influence on phase delay.The S 2 of K p is close to the S 2 of K D .The S 2 of K i is about 0.001 • on magnitude order, much less than the other control parameters.

Experiment
Compared with structural parameters and some working parameters, control parameters can change and be measured during the working process.The two sensitivity indexes of four control parameters are studied by experiments in this section.The mean of several samples is included to ensure the accuracy of the experiment.By comparing the measured data and the simulation result of the first-order sensitivity matrix, the contrast histograms of experiment value and simulation value for the two sensitivity indexes when the four parameters increase 10% under four working conditions are shown in Figures 22-25  As it can be seen in the above Figures 18-21: 1.In all conditions, 1 S is a positive value, so the increase of the four parameters will result in the reduction of the position amplitude attenuation.

Experiment
Compared with structural parameters and some working parameters, control parameters can change and be measured during the working process.The two sensitivity indexes of four control parameters are studied by experiments in this section.The mean of several samples is included to ensure the accuracy of the experiment.By comparing the measured data and the simulation result of the first-order sensitivity matrix, the contrast histograms of experiment value and simulation value for the two sensitivity indexes when the four parameters increase 10% under four working conditions are shown in Figures 22-25  In Figures 22-25, the maximum and mean errors between the experimental value and the simulation value are shown in Table 4.In Figures 22-25, the maximum and mean errors between the experimental value and the simulation value are shown in Table 4.In Figures 22-25, the maximum and mean errors between the experimental value and the simulation value are shown in Table 4.In Figures 22-25, the maximum and mean errors between the experimental value and the simulation value are shown in Table 4.

Figure 2 .
Figure 2. Schematic of HDU performance test platform.

Figure 2 .
Figure 2. Schematic of HDU performance test platform.

Figure 2 .
Figure 2. Schematic of HDU performance test platform.

Figure 2 .
Figure 2. Schematic of HDU performance test platform.

Figure 4 .
Figure 4. Force schematic of the load.

1 pB
refers to the viscous damping coefficient, 1 t m refers to the equivalent mass at piston, and 1 f

1 pB
refers to the viscous damping coefficient,

Figure 7 .
Figure 7. Block diagram of HDU position-based impedance control.

Figure 7 .
Figure 7. Block diagram of HDU position-based impedance control.

Figure 16 .
Figure 16.The curves of dynamic sensitivity between two sensitivity methods.(a) First-order dynamic sensitivity of D K ; (b) First-order dynamic sensitivity of D C ; (c) First-order dynamic sensitivity of

Figure 16 .
Figure 16.The curves of dynamic sensitivity between two sensitivity methods.(a) First-order dynamic sensitivity of K D ; (b) First-order dynamic sensitivity of C D ; (c) First-order dynamic sensitivity of K p ; (d) First-order dynamic sensitivity of K i .

Figure 17 .
Figure 17.Variation of position response when each parameter increases 10%.(a) The first working condition; (b) The second working condition; (c) The third working condition; (d) The fourth working condition.

DK
affects the output position much more than the others.The influence from integral gain i K is the most irrelevant.Desired damping D C and proportional gain p K have similar influences on the output position.With the disturbance force in sinusoidal variation, D C varies by following the curve of minus cosine, and p K varies by following the curve of cosine.

Figure 17 .
Figure 17.Variation of position response when each parameter increases 10%.(a) The first working condition; (b) The second working condition; (c) The third working condition; (d) The fourth working condition.
and Q 2 is outlet oil flow.

Table 1 .
Parameters and initial values of the simulation model.

Table 2 .
Four working conditions researched in this paper.

Table 2 .
Four working conditions researched in this paper.

Table 2 .
Four working conditions researched in this paper.

Table 3 .
Mean values of performance indexes in different conditions.

Table 3 .
Mean values of performance indexes in different conditions.
have some correlation with disturbance amplitude, but have little correlation with disturbance force frequency.2. The 2 S of the four control parameters has positive or negative values in different conditions.An increase of the parameters has different influences on the phase angle delay of the output position.Specifically, having more relation with disturbance force frequency than the amplitude, the 2 S of D K and D C have nearly the same absolute value, but an opposite sign.Their increase has totally different influences on phase angle delay.An approximate 10% increase of them generates about 0.18° of phase angle delay error.The 2 S of p K has correlation with both disturbance force amplitude and frequency.The larger the frequency, the more remarkable the influence on phase delay.The 2 S of p K is close to the 2 S of D K .The 2 S of i K is about 0.001° on magnitude order, much less than the other control parameters.
D K and p K

Table 4 .
Maximum and mean errors between experimental value and simulation value.

Table 4 .
Maximum and mean errors between experimental value and simulation value.

Table 4 .
Maximum and mean errors between experimental value and simulation value.

Table 4 .
Maximum and mean errors between experimental value and simulation value.