A Data-Driven Multi-Scale Online Joint Estimation of States and Parameters for Electro-Hydraulic Actuator in Legged Robot

In order to satisfy the real-time need of model-based controllers for model parameters and full states feedback, this paper has conducted in-depth research on the states and parameters estimation of electro-hydraulic actuator in legged robot with three problems for time-varying parameters estimation (including system parameters and external load force), non-measurable states estimation and measurable states filtering. The first-order trajectory sensitivity method based on the dynamic model is used to determine the parameter set to be estimated, and the parameter fast and slow characteristics are analyzed in detail to obtain the generalized states and slow-varying parameters. Then, the combined algorithm with a fast-varying time scale (composed of a fusion kalman filter and a fast-varying time scale extended kalman filter) and a slow-varying time scale (composed of a slow-varying time scale extended kalman filter) is innovatively proposed to realize the data-driven multi-scale online joint estimation of states and parameters for the actuator system. Finally, the results of three comparative experiments show that the proposed algorithm has better stability, faster convergence speed and more accurate estimation than the dual extended kalman filter algorithm, and the states and parameters estimated by the proposed algorithm accurately reflect the actual characteristics of actuator. Moreover, the algorithm has strong adaptability and robustness in different actuator hardware environment and strong convergence ability for different initial values of states and parameters.


I. INTRODUCTION
Electro-hydraulic actuators are extensively used in heavyduty electromechanical systems and legged robots, for their high load capacity and large power density ratio [1]- [7]. Presently, classical PD control method has been extensively chosen in designing position controller of electro-hydraulic actuators for most legged robots. However, in situations when parameter uncertainties and unknown load disturbances cannot be neglected, the system would exhibit unexpected dynamic behavior and be difficult to meet the performance requirements. In electro-hydraulic actuator system, uncertainties are mostly caused by changes in unknown viscous damping, physical characteristics of the valve, effective bulk modulus and external load force [8].
The associate editor coordinating the review of this manuscript and approving it for publication was Jianyong Yao . Moreover, some parameters and external load changes are significant under different operating conditions [9]. In order to improve tracking performance, several mainstream model-based control schemes have appeared in recent literature, such as adaptive robust control [10], [11], active disturbance rejection control [12], the method combining Nussbaum function and adaptive control [13] and adaptive extended interference observer [14], etc. However, these schemes face some common problems. (1) When the system is disturbed by high frequency and large external load, the adaptive method has the problem of high gain feedback.When system parameters change, the performance of ESO (Extended state observer)-based controller is poor. (2) The impacts of measurement noise during full-state feedback are rarely considered. Practice shows that in some cases, measurement noise has become the core issue in achieving high tracking performance. The measurement signals are usually polluted by heavy noise, which seriously affects the control performance.These schemes are based on the theoretical design of noise-free, however, in practice, various low-pass filters are used to mitigate the effect of noise, causing severe phase lag in the high-frequency range. (3) These schemes are based on full state feedback, which means that in addition to displacement signal, velocity and pressure signals are required. When the actuator cannot install all of these sensors due to cost and/or structural size restrictions, such model-based control methods are difficult to realize directly.
In legged robot applications, the leg performs a large swing motion at a frequency of 3.3 Hz (high frequency). At the mean time, the external load force of actuator changes significantly from −1500(N) to 2000 (N). The system parameters such as equivalent flow coefficient, viscosity coefficient and effective bulk modulus also have large time-varying characteristics. The system parameters and external load force are collectively referred to as parameters. Furthermore, because the robot has set strict upper limit on the structural size and weight of actuator in this paper, only displacement and driving force sensors, and acceleration sensors which must have little volume and weight are configured in the actuators, resulting in the problems that velocity and pressures are not measured and the measurable state (displacement) and signal (driving force) contain noise. Evidently, it is necessary to study the time-varying parameters estimation, unmeasured states estimation and measurable states filtering of the actuator, to lay the foundation for the design of high-precision controllers based on model in the future.
At present, various methods have been proposed to estimate the model parameters of electro-hydraulic system. On the basis of all measurable states, Yuan et al. [15], [16] and Peran [17] have used the open-loop input and output data of electro-hydraulic system collected from the experiment to estimate the parameters of a nonlinear gray box model offline on MATLAB's System Identification Toolbox. Victor et al. [18] has measured the servo spool displacement, flow and input current curves on the experimental platform and estimated the servo valve parameters (port gain) by a fitting method. Moon [19] has adopted the recurrent incremental credit assignment neural network and genetic algorithm to off-line estimate the system parameters of nonlinear electro-hydraulic servo system, such as mass, viscosity coefficient, spring constant and effective bulk modulus. Using the input and output data gathering through experiments, a global solution is obtained.
A common disadvantage of the above methods is that the dynamic changes with varying operating conditions of the model parameters are ignored. Therefore, the reliability and applicability of these estimators are not fully discussed. In order to overcome the disadvantage, online parameter estimation methods are proposed to track the real-time characteristics of actuator. Yao et al. [20] and Ahn et al. [21] have proposed an adaptive robust control, which applies the parameter adaptive estimation algorithm and extended state observer to nonlinear control. Kim et al. [22] has used a high gain PI disturbance observer to estimate and compensate the external load disturbance. Yao et al. [23] have proposed a multilayer neural-networks (NNs) estimator to estimate the disturbance and improve the compensation accuracy of nominal model-based control.The above methods have been fully implemented with the measurable states. Kaddissi et al. [24] has rewriten the system dynamic model in linear parameter form, identified parameters by recursive least square method (RLS) and applied it to the nonlinear backstep algorithm. Sadeghieh et al. [25] has adopted the same parameter estimation method and applied it to the bio-inspired intelligent controller to obtain better control performance than the optimal PID controller. Wos and Dindorf [26] has used the recursive least squares method to identify the nonlinear Hammerstein model of the asymmetric electro-hydraulic servo valve-controlled cylinder system, and the model with nonlinear static block has captured the nonlinear dynamics well. However, the recursive least squares method needs to use the previous data in the recursive calculation, which brings forward higher requirements on the system hardware and software programming, and also requires the estimation model to be linear as well as all states to be measurable. The abundant disadvantages of RLS make it not possible for the joint estimation of states and parameters. Therefore, many researchers have begun to use the optimal filtering methods to deal with the estimation problem, among them the kalman filter is easy to realize and apply in real time on the computer, and can deal with time-varying systems and non-stationary signals. An and Sepehri [27] has applied the extended kalman filter (EKF) to estimate the spool displacement, the chamber pressure of actuator and the velocity of punch, which has a fast and reliable ability to estimate pressure changes. Chinniah et al. [28] has used the extended kalman filter to estimate viscous friction and effective bulk modulus respectively based on sensor data-driven and all states measurability for fault monitoring. Cui et al. [29] has applied the robust kalman filter to indirectly estimate the external leakage of the hydraulic servo motor. However, the above cases are only implemented under the condition that the states are all measurable, but not estimated. Colorado has proposed a novel on-line closed-loop parameter identification algorithm for second order nonlinear systems, designed a cost function based on the optimization method by using a linear combination of the actual and an estimation system, and used algebraic techniques to estimate the velocity and acceleration signals, avoiding noise processing problems. This method converges faster than the online and off-line least squares algorithms, and has strong robustness against disturbance, but without requiring any type of data pre-processing [30], [31]. In addition, he has also proposed a first integrals and adaptive parameter identification method for conservative Hamiltonian systems, and discussed the parameter identification problem as an optimization one [32]. At present, the above two methods are mainly applicable to second order nonlinear systems and conservative Hamiltonian systems.
For state estimation, ESO (Extended state observer) has been used to estimate the unmeasured state [33]. EKF also provides an efficient method for estimating the states of discrete-time nonlinear dynamical system. Nevertheless, the estimation accuracy of states is highly dependent on the predetermined parameter values used in the model. Therefore, based on the three problems in actuator control, the joint estimation of the states and parameters for electrohydraulic actuator is particularly essential. So, in order to make the estimation accuracy of the states more reliable, model parameters need to be constantly calibrated regularly, otherwise the estimation accuracy may be greatly reduced and even lead to invalidation. In other words, when model parameters are not available, a joint estimation algorithm is needed to achieve reliable and accurate estimation of system states and parameters.
In this case, the dual extended kalman filter based on realtime data-driven is proposed, where the states and parameters of the dynamic system are simultaneously estimated by coupling to achieve better prediction accuracy [34]- [37]. The algorithm is an online joint estimation method on the same time scale, however, the estimation errors are relatively much when dealing with the system with a large number of parameters to be estimated and a large difference in the speed of parameter change.
For the three problems of time-varying parameter estimation, unmeasured states estimation, and measurable states filtering of electro-hydraulic actuator under limited sensor configuration, the previous literature review show that it is difficult to obtain good online joint estimation of states and parameters for these methods on this object.
The contribution of this paper is to innovatively propose a data-driven multi-scale online joint estimation algorithm with a fast-varying time scale (composed of a fusion KF (Kalman filter) and a fast-varying time scale EKF (Extended kalman filter)) and a slow-varying time scale (composed of a slow-varying time scale EKF) for states and parameters applied to actuator with three problems. This is the key step and technique to ensure the position accuracy of actuators for model-based controllers. Through comparative analysis, the proposed algorithm has better stability, faster convergence speed and more accurate estimation than the dual extended kalman filter algorithm, and the states and parameters estimated by the algorithm accurately reflect the actual characteristics of actuator. The remainder of this paper is organized as follows. In Section II, the fourth-order gray box state space dynamic model of the actuator is established, the parameter set to be estimated is selected using the first-order trajectory sensitivity, and the change characteristic of parameter set is analyzed in detail. The process and related derivation of data-driven multi-scale online joint estimation algorithm for states and parameters are detailed in Section III. Section IV conducts three comparative experiments on the single-leg experimental platform to fully discuss the estimation performance of the proposed algorithm. Finally, Section V summarizes the main points and future work.

II. DYNAMICES AND PARAMETER ANALYSIS A. DYNAMICES
Compared with the white and black box model, the gray box dynamic model, established based on mechanism, can better describe the model structure of dynamic time-varying parameter system, and also easily reflect the physical meaning of parameters corresponding to the actual system. So the dynamic model of actuator adopts the form of gray box. At the same time, the way to solve the problem that the model structure is too complicated and the computing resources are too large is to find a trade-off between the structural complexity and accuracy. Therefore, the accurate and simple modeling problem is hard to avoid here. In fact, with ensured accuracy, the lower the model order is, the more advantageous it is.
The actuator consists of the servo valve and actuating cylinder. In legged robot application, the servo valve works in the 3.3 Hz frequency band (far less than the natural frequency of servo valve (120 Hz)), so the servo valve is modeled as a proportional link and the higher order dynamics is ignored here. Then the actuator dynamics is composed of the proportional model of servo valve and the fourth-order model of valve-controlled cylinder system. The state variables are set as , and the open-loop state equation is obtained [10], [38]. where x p is actuator output displacement,ẋ p is velocity, P 1 and P 2 are two cavity pressures, K d is equivalent flow coefficient, u is control signal for valve, P s and P 0 respectively are the supply and return pressures, m is the mass of actuator rod, B p is viscosity coefficient, A p1 and A p2 respectively are the piston and rod areas, V 01 and V 02 are the pipeline volumes of piston and rod cavity respectively, L is actuator total stroke, L 0 is piston initial position, β e is effective bulk modulus, c ip is internal leakage coefficient. Under normal operating conditions, the actuator is free from external leakage. F L is the external load force on the piston. The external load force includes inertial force, Coriolis force, gravity, friction force of rigid joint and interference force, which exhibit high dynamic characteristic with the change of cylinder displacement and joint angle. In this paper, the above forces are combined into one as an external load force, which also exhibits high dynamic characteristic. VOLUME 8, 2020

B. PARAMETER SENSITIVITY ANALYSIS
Due to the large number of parameter for electro-hydraulic actuator, the accurate estimation of all parameters is hard to achieve. In fact, only some parameters, called the dominant parameters, exert giant influence on the dynamic performance of position control. Therefore, it is necessary to use the sensitivity analysis method to select the dominant parameters as parameter set, and then employ the estimation algorithm to estimate them. In the paper, the first-order trajectory sensitivity method is used to analyze the influence of parameters on the displacement x p . Under the condition that u and parameters are mutually independent, the first-order trajectory sensitivity equation is as follows [39]- [41].
where λ i n is the first-order trajectory sensitivity function of state vector x to parameter vector α, the mathematical definition is λ i n = [∂x/∂α i ] n , (i = 1, . . . , 12), and its initial value is [∂f /∂x] n is the coefficient term, that is, the partial derivative of state equation to state vector (Jacobian matrix), [∂f /∂α i ] n is the free term, which is the partial derivative of state equation to parameter vector that is Then the first-order approximation function of parameter change α to state change x is It is known from the above equation that to obtain x, it is necessary to first solve the first-order trajectory sensitivity function λ i n corresponding to states, and then multiply it with the parameter vector change α. These equations and variables are calculated on the MATLAB platform.
The actuator position control is applied during the swing period of leg movement. The legged robot performs the desired forward movement with the speed of 1.5 m/s, the duty ratio of 0.5, the gait period of 0.6 second, and the lift height of 0.15 m. Then the operation space planning of swinging leg is performed. Finally the desired trajectory of the hip electro-hydraulic actuator is, as the expected displacement value of dynamic model, obtained according to the inverse kinematics. According to the structure size and measurement, some parameter values such as area are obtained. The empirical values of unknown parameters, such as viscosity coefficient and effective bulk modulus, are determined by referring to other literatures [42]- [45]. The parameter values used in the simulation are shown in Table 1.
The coefficient and free term in Equation (2) contain state variables, so the equation is a differential equation matrix with variable parameters in 4 × 12 dimensions. The calculation process is as follows. Firstly, the proportional controller is used to control the dynamics (Equation (1)), and the numerical solutions of states are obtained. The displacement closed-loop curves are shown in Fig. 1(a). Secondly the states  are introduced into the first-order trajectory sensitivity equation (Equation (2)), jointly calculating the first-order sensitivity function value λ i n of the displacement to parameters. Finally, λ i n is multiplied by the parameter change α to obtain the displacement change x p (Equation (3)). Each parameter is set to change by 10%, and the parameter sensitivity function time-history curves of displacement change x p are shown in Fig. 1(b).  Two indices are used to evaluate the influence of parameter change on displacement. The first one is the maximum displacement change caused by parameter change, and the second one is the integral of displacement change in the sampling period. The histograms are shown in Fig. 2.
In fact, during the operation of actuator, some parameters change very little. Such as, since the constant pressure variable pump is used as the oil source, and the maximum output flow of pump is greater than the maximum flow demand for the simultaneous movement of actuators, the constant pressure setting and flow output ensure that the supply pressure α 2 (P s ) fluctuates very little. It is determined in the design that the change of return pressure α 3 (P 0 ) is small because it is highly related to the design and selection of return pipeline.The parameter fluctuations caused by the machining errors of α 5 (L), α 6 (L 0 ), α 7 (A p1 ) and α 8 (A p2 ) are much less than 1%. So in the absence of structural damage, the structural parameters are fixed by default. Therefore, the new histograms are obtained without considering the sensitivity results of above pressures and structural parameters, as shown in Fig. 3. From Fig. 3(a), it is found that the index I of equivalent flow coefficient α 1 (K d ), effective bulk modulus α 9 (β e ), viscosity coefficient α 11 (B p ) and external load force α 12 (F L ) are all greater than 1.1e −5 (m), wherein the value of equivalent flow coefficient α 1 (K d ) is the largest. From Fig. 3(b), the equivalent flow coefficient α 1 (K d ), effective bulk modulus α 9 (β e ), viscosity coefficient α 11 (B p ) and external load force α 12 (F L ) have sustained influence on the displacement. However, the index II of internal leakage coefficient α 4 (c ip ) and rod mass α 10 (m) are far less than those of above four parameters, so their influences on the system don't last. From the comparison, the four parameters of equivalent flow coefficient α 1 (K d ), effective bulk modulus α 9 (β e ), viscosity coefficient α 11 (B p ) and external load force α 12 (F L ) vary greatly with the load, hydraulic oil temperature and pressure, commutation friction switching and other factors during the control process, and they have an important influence on the displacement accuracy. Therefore, the four parameters are selected as the parameter set (K d , β e , B p , F L ) for online joint estimation.

C. THE FAST AND SLOW CHARACTERISTIC ANALYSIS OF PARAMETER
The parameter set (K d , β e , B p , F L ) obtained so far is distributed in the force balance equation and flow continuity equation of dynamic model (Equation (1)), respectively.
The mathematical expression of equivalent flow coeffi- where C d is the flow coefficient of valve port, w is the area gradient of spool valve and ρ is the density of hydraulic oil. Ideally, the equivalent flow coefficient K d is a fixed value. However in practice, when the valve structure has been determined, K d exhibits a small and slowly varying process with the small change of fluid Reynolds number [46]- [47], so it is a slow-varying parameter. The effective bulk modulus β e has a value range of about 1 × 10 9 − 2 × 10 9 (Par). In high dynamic motion, it is mainly compressed by pressure to show high dynamic fluctuation at the same frequency as pressure, and is a fast-varying parameter.
The viscosity coefficient B p is sensitive to temperature change. After the actuator starts from a cold state, the temperature of hydraulic oil slowly rises to a certain value (the hydraulic radiator controls the temperature within a certain range). At the same time, the viscosity coefficient of hydraulic oil is gradually reduced to a fixed value. When temperature shows a slowly changing characteristic, the viscosity coefficient B p is also a slow-varying parameter. The detailed expression for external load force is F L = J · (M (q)q + C(q,q)q + G(q) + τ f ), where M , C and G are inertial matrix, Coriolis matrix, gravity vector respectively, τ f contains the frictional and interfering forces of rigid joints and unmodeled dynamics. q,q andq are the joint angle, angular velocity and angular acceleration of mechanism movement, and they correspond exactly to the actuator displacement, velocity and acceleration. J is the Jacobian matrix converted between the joint angular velocity and actuator velocity, and the relationship isẋ p = Jq. From the above detailed expression, it is concluded that the external load force is a fast varying parameter at the same frequency as the state displacement and velocity of actuator.
From the above characteristic analysis, the fast varying parameter set is θ fast = [β e , F L ] T , the slow varying parameter set is θ slow = [K d , B p ] T , and the actuator states change rapidly with the high dynamics of operating conditions. Based on the characteristics, the paper sets the fast varying parameter set θ fast = [β e , F L ] T and the states [x p ,ẋ p , P 1 , P 2 ] T on the same fast scale for estimation, and the slow varying parameter set θ slow = [K d , B p ] T on the slow scale for estimation. VOLUME 8, 2020

III. THE DATA-DRIVEN MULTI-SCALE ONLINE JOINT ESTIMATION ALGORITHM FOR STATES AND PARAME-TERS
Under the condition of limited sensor configuration, aiming at the three problems for time-varying parameters estimation, non-measurable states estimation, and measurable states filtering of electro-hydraulic actuators, based on the fast and slow changing characteristics of states and parameters analyzed in the previous section, this paper innovatively proposes a data-driven multi-scale online joint estimation algorithm with a fast-varying time scale (composed of a fusion KF and a fast-varying time scale EKF) and a slow-varying time scale (composed of a slow-varying time scale EKF) to realize the real-time online estimation of actuator states and parameters.

A. SYSTEM DESCRIPTION
Firstly, the two measured values are the acceleration and displacement signal of the axial movement for piston rod (the same direction of two sensors). In this paper, the piston rod is assumed to be a particle, and the axial displacement, velocity and acceleration of piston rod are at the same centroid point. Then the axial linear motion of piston rod is described by the particle motion equation. After obtaining the real-time signal of displacement and acceleration with noise, the particle motion equation is used to estimate the particle velocity as the axial velocity of piston rod, and filter the displacement measurement signal. The measured acceleration is taken as the input signal u 1,k,l , the measured displacement is set as the output signal x p , and the state variables are set as Then the discrete time state equation is as follows. χ 1,k,l+1 = F 1 (χ 1,k,l , u 1,k,l ) + ω 1,k,l y 1,k,l = G 1 (χ 1,k,l , u 1,k,l ) + υ 1,k,l .
where χ 1,k,l is the state vector at time t k,l = t k,0 + l × t(1 ≤ l ≤ L z ), time scale k and l respectively describe slow-varying and fast-varying time scale, L z is the scale conversion limit, that is, one slow-varying time scale is equal to L z fast-varying time scales, t is the calculation time interval. u 1,k,l is the measured acceleration at time t k,l . ω 1,k,l and υ 1,k,l are process and measurement noise matrix respectively, whose corresponding covariance matrices are Q χ 1 and R χ 1 . F 1 (χ 1,k,l , u 1,k,l ) is the transition matrix, and G 1 (χ 1,k,l , u 1,k,l ) is the measurement matrix.
Secondly, the state vector of actuator dynamic equation (Equation (1) In order to distinguish them from the states and parameters with common meanings, χ 2 is collectively referred to as the generalized state vector and θ is referred to as the slow varying parameter set. Then, a multi-scale nonlinear discrete state space model including the generalized states and slow-varying parameters is obtained as follows.
It is important to note here that the displacement and velocity in the measurement matrix are state estimates based on the state equation (4). The specific reason for the usage is detailed in Section III-B. Where χ 2,k,l is the state vector at time t k,l = t k,0 + l × t(1 ≤ l ≤ L z ), u 2,k,l is the input signal at the same time (servo valve control signal). y 2,k,l is the measurement vector at time t k,l . G 2 (χ 1,k,l , χ 2,k,l , θ k , u 2,k,l ) is the measurement matrix. ω 2,k,l and ρ k are the process noise matrix for generalized states and slow-varying parameters respectively, whose covariance matrices are Q χ 2 and Q θ respectively. υ 2,k,l is the measurement noise matrix whose covariance matrices is R χ 2 . Based on the definition of electro-hydraulic actuator system, the goal is to estimate the generalized states χ 2 and slow-varying parameters θ from the measurement data y 1 and y 2 , which contain the acceleration, displacement and driving force with noise. The generalized states refer to the filtered displacement x p , velocityẋ p , two cavity pressures P 1 , P 2 , effective bulk modulus β e and external load force F L . The slow-varying parameters refer to equivalent flow coefficient K d and viscosity coefficient B p . The generalized states are on the fast-varying scale, and the slow-varying parameters are on the slow-varying scale.
The two-scale estimators perform stepwise estimation of the generalized states and slow-varying parameters, and the two are performed alternately with each other as input. Moreover, the estimators use the innovation from same source. The algorithm has a coupling structure that guarantee a stable closed-loop estimation of the final generalized states and slow-varying parameters. And because the state innovation is used, the algorithm adapts the state estimation through the deployment of model parameters on the basis of guaranteeing the state estimation effect. The advantage of the proposed algorithm is that it fixes two slow-varying parameters in the fast-varying time scale, reduces the dimension of generalized states that need to be estimated simultaneously and improves the estimation convergence. The generalized state dimension of algorithm is six, while the generalized state dimension of the ekf and dekf algorithms is eight.
The overall framework of algorithm is as follows. (1) In the fast-varying time scale, since there is no sensor to measure the stateẋ p in actuator's state space model (Equation (5)  The specific calculation steps of the proposed algorithm are summarized as follows and in the flowchart as shown in Fig. 4. (1) Step 1: Initialization, set the initial parameters of fusion KF, fast-varying time scale EKF and slow-varying time scale EKF, respectively.
where χ 1,0,0 , P χ 1 0,0 , Q χ 1 , R χ 1 are respectively the initial state vector, the initial value of state estimation error covariance matrix, the process noise covariance, and the measurement noise covariance of fusion KF. χ 2,0,0 , P χ 2 0,0 , Q χ 2 are VOLUME 8, 2020 respectively the initial generalized state vector, the initial value of state estimation error covariance matrix, and the process noise covariance of fast-varying time scale EKF. θ 0 , P θ 0 , Q θ are the initial slow-varying parameter set, the initial value of parameter estimation error covariance matrix, and the process noise covariance of slow-varying time scale EKF, respectively. R χ 1 and R θ are the measurement noise covariances. Since the same innovation is used, R χ 2 = R θ is satisfied. dχ 2,0,0 /dθ − 1 ∈ R 6×2 is set as a matrix with zero elements for calculation of parameter measurement matrix C θ k . When the estimation starts, the value at time (0) is converted to the value at time (k − 1), and the value at time (0,0) is converted to the value at time (k − 1, l − 1).
For calculation of parameter measurement matrix C θ k . (5)Step 5: Cycle calculation of fast-varying time scale for l = 1 : L z . When the cumulative count is equal to L z , scale conversion is performed to activate the calculation of slowvarying time scale. At this time, make the following switch. χ 2,k,0 =χ 2,k−1,L z , y 2,k,0 = y 2,k−1,L z , (6)Step 6: The posterior estimation of slow-varying time scale EKF where So far, the multi-scale online joint estimation of the generalized states and slow-varying parameters at time k is completed, and then it is ready to enter the cycle at time k +1(end).
The calculation for measurement matrix C θ k of slowvarying parameters is the most important step for ensuring the convergence of coupling estimation, and it is the total derivative of measurement function with respect to the slow-varying parameters. Considering that the generalized states are the functions of slow-varying parameters, the total derivative is decomposed into partial derivatives as follows.
In step 4: The posterior estimation of fusion KF and fast-varying time scale EKF In step 5: Scale conversion is performed to activate the calculation of slow-varying time scale. At this time, make the following switch.

IV. EXPERIMENT AND DISCUSSION
In order to evaluate the performance of the proposed algorithm, the one-legged motion control experimental platform has been built as shown in Fig. 5. The platform consists of the single-legged rigid body, actuators and control system, in which the electro-hydraulic actuators driving the joint VOLUME 8, 2020 mechanism have integrated the acceleration, displacement and force sensor. The control system consists of a computer, PC104 small board, ARM controller, amplifier and 16-bit A/D converter sensor. The QNX operating system with 1000Hz control frequency runs in the controller, where the control and estimation algorithms are implemented.

A. DESCRIPTION OF THE EXPERIMENTAL PROCESS
In the experiment, the hip of single leg is fixed on the bracket to suspend the foot and simulate the leg movement of legged robot in the swing phase.Then the knee and hip actuators are controlled and estimated. Firstly, in the motion controller, the same desired displacement trajectorys of actuators in the sensitivity analysis are used, and then the closed-loop position control is performed using the proportional controller to stabilize actuators. The accuracy of proportional controller is not high, the main purpose of which is to activate the dynamic characteristics of actuator, and to avoid the divergence of open loop control, as well as to prevent damage to the components. Secondly, at each estimation time, the acceleration, displacement and driving force data measured in real time are applied to the proposed algorithm to realize the online joint estimation of generalized states and slow-varying parameters. In order to improve the efficiency of algorithm initial parameter adjustment and performance comparison, and to realize the effective comparison of different estimation results for the same data, the sensor data measured in the experiment are saved to the computer, then the online joint estimation of generalized states and slowvarying parameters is realized based on the matlab / simulink platform.
The algorithm initial parameters are divided into five types as follows. (1) Measurement noise covariances R χ 1 , R χ 2 and R θ are the characteristics related to the sensors. It should be noted that if this values are too large or too small, the filtering effect will be very poor. And the closer the values are to the actual noise, the faster the algorithm converges. Their appropriate values need to be found gradually through experiments. (2) State and parameter process noise covariances Q χ 1 , Q χ 2 and Q θ . Because the state innovation is used, the algorithm adapts the state estimation through the deployment of model parameters on the basis of guaranteeing the state estimation effect. Therefore, the state process noise covariances are set to be small, and only the parameter process noise covariance need to be adjusted through experiments. (3) The estimated error covariances P χ 1 0,0 , P χ 2 0,0 and P θ 0 . As long as they are not zeros, the values have little effect on the filtering effect and converge quickly. So they are set as the certain values without too much analysis. (4) The time scale L z . Too large or too small value can affect the estimate effectiveness. (5) The initial value of states vectors χ 1 , 0, 0 and χ 2,0,0 , the initial value of slow-varying parameter vector θ 0 . Because the true values are unknown, the estimated initial values are generally set based on prior knowledge, but are not accurate. It is necessary to rely on the convergence ability of the proposed algorithm to make their estimated values quickly converge to the true ones.
The adjustment method of initial parameters is as follows. The state process noise covariances(Q χ 1 , the first four dimensions of Q χ 2 ), the estimated error covariances (P χ 1 0,0 , P χ 2 0,0 , P θ 0 ) and the estimated initial values (χ 1,0,0 , χ 2,0,0 , θ 0 ) are set to fixed values, and then the measurement noise (R χ 1 , R χ 2 , R θ )and parameter process noise covariances(the last two dimensions of Q χ 2 , Q θ ), and time scale(L z ) are adjusted by experiments. The above three parameters with a gradually increasing trend are respectively set, and then the optimal values are determined by comparing the effect of estimation.
Evaluation method of algorithm performance: Since generalized state and slow-varying parameters are estimated by coupling with each other, in the case where the estimated state follows the centerline of measurement value, the displacement x p and velocityẋ p are directly related to each other, as are the driving force and the two cavity pressure P 1 , P 2 . Therefore, the comparison between the estimated values of displacement and driving force and the center line of measured ones, combined with the prior knowledge of model parameters and the relationship between the states and parameters are used to judge the estimation effect.
The prior knowledge is that the accuracy error of displacement and driving force sensors is ±0.5%, the generalized state (fast-varying parameter) β e is changed by pressure extrusion in the range of β e ∈ (1e 9 , 2e 9 )(Par), its frequency is similar to pressure and driving force. F L is directly related to the driving force and pressure, and the change frequency is the same. The slow-varying parameter K d has a small nonlinear fluctuation characteristic reflecting the flow capacity of valve port. As the hydraulic oil stabilizes at a certain temperature, B p stabilizes around a certain value in the range of B p ∈ (0, 4e 3 )(N .s/m). In the experiments, the evaluation indices used to estimate performance are mean square error (MSE), maximum error (ME) and average error (AE).

B. ALGORITHM PERFORMANCE VERIFICATION
The section mainly shows the detail of algorithm performance verification. It is divided into three experiments, which are (1) the performance comparison with the combined algorithm (dekf) of fusion kalman filter and dual extended kalman filter, (2) verification of estimated performance on different actuator hardware, and (3) convergence verification of different initial values for generalized states and slow-varying parameters. Experiment (1) and (2) both include comparison with dekf to verify the performance and advantages of the proposed algorithm on different hardware, and experiment (3) shows the ability of the proposed algorithm to converge from different initial values.
Experiment (1) and (3) have used the 50 seconds measurement data of hip actuator, which are the relative acceleration of piston rod to cylinder, the displacement and driving force as shown in Fig. 6. Since the slow-varying parameter B p fully converges near 50 seconds, the data length of 50 seconds  is used. Due to the density and periodicity of measurement data, only the first 10 seconds data is shown in the figures in order to display the data clearly. Experiment (2) has used the 50 seconds measurement data of knee actuator, as shown in Fig. 7 (the length of displayed data is the same as that in Fig. 6). In the data, the original acceleration signal contain certain noise, especially the high-frequency interference caused by vibration makes the signal to have a large deviation. The first-order low-pass filter is used to filter out the high-frequency random signal in the signal, and then the filtered acceleration signal is imported into the proposed algorithm in real time. The displacement and driving force signals are noise-containing measurement values directly measured by sensors, the data detail has been shown in the experiments.
(1) The performance comparison with the combined algorithm (dekf) of fusion kalman filter and dual extended kalman filter In order to verify the estimation effect and performance, the proposed data-driven multi-scale online joint estimation algorithm (mekf) has been implemented on the hip actuator, and compared with the combined algorithm (dekf) based on fusion kalman filter and dual extended kalman filter. The dimension of generalized states in the proposed algorithm is six, while in the dekf algorithm, the eight dimensions of generalized states and slow-varying parameters are in the same fast-varying time scale. The initial parameters of the proposed algorithm (mekf) on hip actuator are obtained by the method described in the fourth paragraph of Section IV.A as follows.
The initial parameters of dekf algorithm are the same as those of mekf. The correlation estimation comparison curves of generalized states χ 2 = [x p ,ẋ p , P 1 , P 2 , β e , F L ] T (Fig. 8,10,11), driving force ( Fig. 9) and slow-varying parameters θ = [K d , B p ] T (Fig. 12) are obtained respectively. As the estimated value of B p converges slowly, the 50 seconds data is shown in Fig. 12 (b). Since other estimated values have converged in about 2 seconds, in order to show the data more clearly, other figures only show the first 10 seconds data. The quantitative performance indices of displacement and driving force under the two algorithms are shown in Table 2.     velocityẋ p and pressures P 1 , P 2 converge to a stable interval quickly after starting from the initial value. The velocitẏ x p changes smoothly between the maximum amplitude of −0.1 − 0.45(m/s), the pressures P 1 and P 2 change smoothly in the range of 5e 6 − 13e 6 (Par) and 8e 6 − 15e 6 (Par) respectively, and their frequencies and phases are consistent with the displacement and driving force respectively. The estimated β e stabilizes quickly in the range of 1e 9 −1.7e 9 (Par), and the frequency is close to the pressures, which is consistent with the change characteristic of effective bulk modulus for hydraulic oil during compression. The estimated F L stabilizes quickly in the range of −2000 − 1500(N ), the frequency and phase are consistent with the driving force, which corresponds to their relationship and also reflects the dynamic characteristic of hip joint dynamic load on the actuator. The estimated K d starts from the initial value, quickly stabilizes and slowly changes in the range of 5.5e −8 − 7e −8 (m 3 /(s.V )) showing the nonlinear slow-varying characteristic. After 20 seconds of adjustment, the estimated B p have gradually converged to around 1200(N .s/m). The parameter is mainly related to the temperature of hydraulic oil. With the system running for a long time, the oil temperature slowly increases and stabilizes to a constant value. Then B p decreases accordingly and finally stabilizes near a fixed value. Therefore, the estimated B p in this experiment is consistent with the actual characteristics. It is clear that K d and B p are updated every 0.2 second, obtained by multiplying the fast-varying scale interval time of one millisecond and L z = 200.
Performance evaluation of the dekf algorithm: The estimated driving force has failed to track the center line of measured value from 3 seconds, and the wave dynamic potential appeared in the overall situation, resulting in similar changes in the pressure P 1 and fast-varying parameter F L . From the direct correlation between them, it is considered that the driving force, pressure P 1 and the fast-varying parameter F L estimated by the dekf algorithm do not conform to the actual dynamic characteristics. The slow-varying parameter B p shows a large fluctuation and does not converge.
It has been concluded from the curves that the estimated values of dekf algorithm are overshooting seriously in the initial stage, while the ones of mekf algorithm have almost no overshooting and converged quickly and smoothly, ensuring the stability of adding the estimated values to the model-based controller in the future. In addition, the quantitative indices of mekf are smaller than those of dekf. The MSE index of driving force is only 31.4% that of dekf, the ME index ratio is 33.1%, and the AE index ratio is 30.9%. Therefore, the mekf algorithm has better stability, faster convergence speed and more accurate estimation than the dekf algorithm. Moreover, the generalized states and slow-varying parameters estimated by mekf algorithm accurately reflect the actual characteristics of actuator.
(2) Verification of estimated performance on different actuator hardware In order to verify the estimated performance in different hardware environments, the proposed algorithm (mekf) has been implemented simultaneously on the knee actuator, and compared with the combined algorithm (dekf) based on fusion kalman filter and dual extended kalman filter. The initial parameter adjustment method of the proposed algorithm for knee actuator is the same as above, and the configuration after adjustment is as follows.   The correlation estimation comparison curves of generalized states χ 2 = [x p ,ẋ p , P 1 , P 2 , β e , F L ] T (Fig. 13, 15, 16), driving force (Fig. 14) and slow-varying parameters θ = [K d , B p ] T (Fig. 17) are obtained respectively. The quantitative performance indices of displacement and driving force under the two algorithms are shown in Table 3. It is also concluded from the curves and quantified indices that the mekf algorithm has better stability, faster convergence speed and more accurate estimation than the dekf algorithm. Moreover, the generalized states and slow-varying parameters estimated by mekf algorithm also accurately reflect the actual characteristics of actuator. According to the estimation results on different actuators, the proposed algorithm has strong adaptability and robustness in different hardware environments.
(3) Convergence verification of different initial values for generalized states and slow-varying parameters In the initial value configuration, the true values of generalized states and slow-varying parameters are unknown. The estimated initial values are generally set based on prior knowledge, but not accurate. Therefore, in the case of inaccurate initial values, it is necessary to rely on the convergence ability of the proposed algorithm to make their estimated    values quickly converge to the true ones. In this experiment, the convergence ability is verified on the hip actuator. Other parameters are the same as experiment (1). Since the actuator is started from the static state, the initial displacement and velocity are determined to be zeros. Five groups initial generalized states and slow-varying parameters are set for estimation and comparison, as follows.
From the curves, it is concluded that after setting five groups initial values with large differences, the estimated    displacement and driving force all have converged to the same curves after 2 seconds and the estimation errors after the time are exactly the same. The estimated velocity is completely recombined into one curve. The generalized states P 1 , P 2 , β e , F L and slow-varying parameters K d all have started from different initial values and quickly converged to the corresponding true values after finite amplitude adjustment of 2 seconds. Only B p has a longer convergence time. It is verified that the proposed algorithm has strong convergence ability for different initial values of generalized states and slow-varying parameters.

V. CONCLUSION
The fourth-order gray box state space model based on mechanism is constructed and the parameter set to be estimated is selected through the first-order trajectory sensitivity analysis. The change characteristics of parameter set are analyzed in detail, and it is divided into the fast-varying parameter set and slow-varying parameter set. And then combined with the fastvarying characteristics of states, the fast-varying parameters and states are divided into fast-varying time scale, and the slow-varying parameters are divided into slow-varying time scale. A data-driven multi-scale online joint estimation algorithm with a fast-varying time scale (composed of a fusion KF and a fast-varying time scale EKF) and a slow-varying time scale (consisted of a slow-varying time scale EKF) is innovatively proposed to realize the efficient and accurate states and parameters estimation of electro-hydraulic actuator system with three problems of time-varying parameter estimation, unmeasur-ed states estimation, and measurable states filtering. Finally, the results of three comparative experiments show that the proposed algorithm has better stability, faster convergence speed, more accurate estimation than the dekf algorithm, and the states and parameters estimated by the proposed algorithm accurately reflect the actual characteristics of actuator. Moreover, the algorithm has strong adaptability and robustness in different actuator hardware environments, and strong convergence ability for different initial values of states and parameters.
In the future, the performance and generalizing ability of the proposed algorithm will be further verified, and it will be used in model-based controller, such as feedback linearization and other model control algorithms, so as to further improve the control accuracy of actuator system under strong time-varying nonlinearity and interference, and to improve the adaptability of control frame.