Mass estimation of a simple hydraulic crane using discrete extended Kalman filter and inverse dynamics for online identification

Automatization of hydraulic machinery requires accurate information of the current dynamic state of the machinery but also information of the underlying dynamic model characterized by a set of parameters. Some of the parameters can be considered static and well defined, such as machinery dimensions, whereas a part of the parameter set is time varying and needs to be identified based on observations. Particularly, difficult parameters to estimate are the ones, from which no prior knowledge is available. Consequently, the parameter corrections cannot be assumed to be small, which is commonly required for the existing parameter estimation algorithms. This study creates an online capable identification algorithm for estimation of a load mass operated by a hydraulic crane. In the case of load mass estimation, the unknown parameter can be practically any positive value, which implies the parameter corrections to be large. In this study, the estimation problem is divided in two parts: First, the dynamical states of the system are estimated based on the system kinematic relationships and dynamics of the hydraulic circuit. Secondly, the unknown load mass is estimated based on the known hydraulic forces and kinematics using the inverse dynamics of the mechanical structure. The proposed algorithm is tested with both artificially created measurements and with an experimental setup. The results show that both the kinematics of the structure and hydraulic pressures can be accurately estimated using the proposed method. Moreover, the method can be used to further estimate the payload mass. A drawback related to inverse dynamics is that it produces biased estimates in static equilibrium because of the discontinuous nature of static friction force. However, this drawback can be avoided, in part, by not updating the payload estimate in the low-velocity region. The proposed estimation methodology is capable for online identification, and as such, it can be used to adapt the control laws of automated machinery. Moreover, the methodology can be useful to record and document the amount of payload being handled during a work cycle.

based on the known hydraulic forces and kinematics using the inverse dynamics of the mechanical structure.The proposed algorithm is tested with both artificially created measurements and with an experimental setup.The results show that both the kinematics of the structure and hydraulic pressures can be accurately estimated using the proposed method.Moreover, the method can be used to further estimate the payload mass.A drawback related to inverse dynamics is that it produces biased estimates in static equilibrium because of the discontinuous nature of static friction force.However, this drawback can be avoided, in part, by not updating the payload estimate in the lowvelocity region.The proposed estimation methodology is capable for online identification, and as such, it can be used to adapt the control laws of automated machinery.Moreover, the methodology can be useful to record and document the amount of payload being handled during a work cycle.

Introduction
The great ambition to increase the level of automation in mechatronic machinery aims to produce more effective, efficient and safe machines.However, for the effective control and automation of the machinery, the unknown parameters of the system need to be identified.A parameter that is probably one of the most inherit to vary is the mass of the payload of a machinery.The payload mass has a great significance on system dynamics and its identification is needed to tune the controller gains and safety limits accordingly.Moreover, the weight of the handled objects is also under great interest, since the information can be used, e.g., to document the total mass of loaded material.
To identify a parameter affecting the dynamical behavior of a system, such as the load mass, the system dynamical state needs to be determined.Commonly, the dynamical state of the system cannot be directly measured, but a state estimation algorithm can be applied.Kalman filter is a commonly used method to estimate states of a system described by a set of linear differential equations [8].However, systems including high nonlinearities require its nonlinear variant to be used [30].In estimation of nonlinear mechanisms, extended Kalman filter [29] and unscented Kalman filter [16] are commonly applied as in [24].However, the drawback of unscented Kalman filter is its higher computational cost compared to extended Kalman filter [24], which limits its applicability in online estimation applications.
The application of extended Kalman filter into multibody dynamics was pioneered over a decade ago by Cuadrado et al in [4,6], where the extended Kalman filter was applied in its continuous time form.Later, the extended Kalman filter has been applied in multibody systems in its discrete time form in [28,31] and in error state form [27].The error state extended Kalman filter has been further expanded to a mechanical system including coupling with hydraulics [12].Unscented Kalman filter has been proposed as an alternative for linearization-based extended Kalman filter [24], however, with the disadvantage of increased computational burden due to the computation of sigma points [24] [21].All the aforementioned studies incorporate the equations of motion of the system to the estimator structure to predict the evolvement of the states and, as such, use the equation of motion as an a priori knowledge for the velocities at next time step.
Another approach in the state estimation of mechanisms is the use of so-called kinematic Kalman filter (KKF) [22].In this approach, the forces affecting the system are assumed unknown and the estimation is based on purely the measurements of kinematic quantities and known constraints of the motion.This methodology has been further applied in estimation of unknown forces by pairing the kinematic estimator with a dynamic estimator, which incorporates the equation of motion for the estimation of the unknown forces [23].
Parameter estimation methods can be classified by multiple factors, such as in online or offline methods, in parametric or nonparametric methods [19] and in frequency-domain or time-domain methods [20].The selection between these methods is case specific and depends highly on the purpose of the identification.In general, online identification is more demanding than offline identification because of the data available and computational requirements.However, the online identification is required in many applications, where the system needs to be observed in real time.
The nonparametric methods differ from the parametric methods by the coefficients of the identified model having no clear physical interpretation but producing similar input-output behavior with the physical system.For this reason, parametric methods are more suitable in applications where also the physical background of the system dynamics needs to be known.
Frequency-domain methods are based on the idea of modeling a system as a signal amplifier for specific input and output signals.Consequently, the system dynamics can be characterized by comparing the frequency spectra of system inputs and outputs.However, frequency-domain methods often lead nonparametric system identification [20], where the physical quantities are difficult to determine.Moreover, the frequencydomain methods require a proper excitation signal to be fed to the system, which, in the case of hydraulic machinery, is often challenging.
Another classification for state and parameter estimation is related to whether they are performed jointly (joint estimation) or separately (dual estimation).In joint parameter and state estimation, the unknown parameter vector is fused with the original state vector and the resultant augmented state vector is estimated [32].In dual estimation, the states of the system are estimated separately from the parameters and the resultant state estimates are used as an input for the parameter estimation [32].Correspondingly, the estimated parameters can be fed back to the state estimation algorithm.
The research problem of this study was motivated by a practical problem.The problem was to identify a payload mass in an experimental setup of simple hydraulic crane by using only position sensors and pressure sen-sors of the hydraulics.A similar problem was presented in [23] where excavator digging forces were estimated based on measured positions, accelerations and actuator forces.However, that methodology was not applicable for the problem of the present study, because of lack of acceleration data and due to payload mass appearing in the equation of motion in both external forces and inertial forces of the system.
The objective of this study is to create an online identification algorithm for an unknown payload of a hydraulically actuated machinery.The algorithm is to be capable of online identification, meaning that only the information from the past shall be used.The estimator follows a similar two stage estimation structure as in [23].However, the missing acceleration data are replaced by incorporating the measurements from the hydraulic pressures and the known dynamics of the hydraulic circuit.Moreover, the dynamic observer is reformulated so that instead of estimating the unknown forces, the payload mass is to be estimated.Additionally, the state transition model to be used utilizes linearization and exponential integration, originally proposed in [25], to formulate a system of linear differential equations at each time step, which is integrated with a single step exponential integration.It is important to note that the proposed method differs from the dual estimation method [32], by not incorporating the estimated parameter to the state estimation algorithm.The reason for this is to avoid closed-loop coupling between parameter estimation and state estimation which may lead stability problems.Moreover, the methodology differs from the estimators coupling multibody dynamics with hydraulics [12] [17] by not including equations of motion to the Kalman filter algorithm.

Modeling of hydraulic machinery
A typical construction of a hydraulic machinery consists of passive structural elements, or bodies, linked together with joints, and hydraulic actuators determining the movement between these bodies.The dynamics of these systems are coupled together since hydraulic pressures produce the external forces for the body structure and correspondingly the movement of the bodies affects the hydraulic pressures through the actuator displacements.In this study, lumped fluid theory is used to evaluate the dynamics of the hydraulic circuit and that model is later used for the estimation of hydraulic Fig. 1 Illustration of a closed-loop topology system forces and kinematics.Furthermore, multibody dynamics is used to model the structure of mutually connected bodies, when the external forces from the hydraulics are known.This multibody model is later used to solve the unknown payload mass of the system.

Modeling of mechanics: multibody dynamics
The mechanics of hydraulic machinery can be modeled using a topological multibody dynamics method [15].The dynamics in this method are first formulated in the Cartesian coordinates and then expressed in the relative joint coordinates.The reference point in this study is a point on the moving body that instantaneously coincides with the origin of the inertial reference frame.Closed-loop topologies are temporarily made open-loop topologies by introducing cut joints, and later, the corresponding loop-closure constraints are imposed accordingly.In this study, the loop-closure constraints are incorporated using the coordinate partitioning method and this approach is also known as the double-step semi-recursive method.The methodology is comprehensively presented in [5], [14] and [11].Therefore, the actual derivation of equation of motion is excluded from this study.However, the main points of the methodology are presented here for the convenience of the reader.
Consider a tree-structure open-loop topology system of N b rigid bodies with any number of branches as shown in Fig. 1.The open-loop kinematics can be recursively expressed using the classical kinematic relations [14,15].Accordingly, the Cartesian velocities Z ∈ R 6N b ×1 and accelerations Ż ∈ R 6N b ×1 can be transformed into a set of relative joint velocities ż ∈ R N b ×1 and accelerations z ∈ R N b ×1 using a velocity transformation matrix R ∈ R 6N b ×N b as [5,14]: where T ∈ R 6N b ×6N b is a constant path matrix and R d ∈ R 6N b ×N b is a block diagonal matrix containing joint specific vectors b [5].Using the previous definitions, the equation of motion for the open-loop system can be written as: where M and Q are a mass matrix and force vector, respectively, with respect to aforementioned special type of reference point coordinates, while M and Q are the mass matrix and force vector, respectively, in the joint space.
To incorporate cut joints for a closed-loop topology system with N f degrees of freedom, a set of N m loop-closure constraints (z) = 0 are imposed in Eq. (3) using the coordinate partitioning method [3,14], where z ∈ R N b ×1 is the vector of relative joint coordinates.The constraints here are assumed to be holonomic and scleronomic.The relative joint coordinates can be transformed into a minimum set of relative joint coordinates using a velocity transformation matrix R z ∈ R N b ×N f as [14]: where żd ∈ R N m and żi ∈ R N f are the respective dependent and independent relative joint velocity vectors, zi ∈ R N f is the independent relative joint acceleration vector and d z ∈ R N m ×N m and i z ∈ R N m ×N f are the dependent and independent columns of the Jacobian matrix z , respectively.The inverse of d z is assured by ignoring redundant constraints and singular configurations.By substituting Eqs. ( 4) and (5) in Eq. ( 3), the equations of motion for the closed-loop topology system can be written as: where M and Q are the mass matrix and force vector, respectively, in the independent joint coordinate space and the term Ṙz żi represents the Cartesian accelerations under true relative joint velocities and zeroindependent relative joint accelerations.

Lumped fluid theory
Lumped fluid theory is a common approach in modeling of fluid power systems [33].In lumped fluid theory, the fluid power system is divided into finite control volumes with uniform pressures and zero volume elements carrying the pressure losses of the circuit.In the case of a simple hydraulic circuit, the control volumes are composed of the hydraulic lines and cylinders, while the pressure losses are concentrated in the hydraulic valves, representing the borders between the control volumes.This assumption allows the mass density inside a control volume to be expressed using a single scalar value, and consequently, the dynamics of the fluid can be described by applying the theory of conservation of mass as: where m and ρ are the total fluid mass inside the control volume and its density, respectively, Q tot is the total flow rate into the control volume and V is the control volume volumetric size.In this approach, it is assumed that the deformation of the pressurized structure does not create any change in the control volume size, but the effect of deformation is considered as a part of the effective bulk modulus [33].The effective bulk modulus B e describes the relationship between the fluid density ρ and the pressure p as: which can be applied to Eq. ( 7) resulting: The effective bulk modulus can be evaluated using the definition of compressibility being the inverse of bulk modulus and the property of summability of compressibilities of the hydraulic fluid and its surrounding mechanical structure elements.The contribution of different structural elements (hoses, cylinders, etc.) is taken into account by their volumetric weighted average as: where κ e is the effective compressibility of the hydraulic oil and the structure, B oil is the bulk modulus of the hydraulic oil, n c is the number of subvolumes and V k and B k are their respective volumes and bulk modules.
Because the deformation of the structure is taken into account in the effective bulk modulus and it can be assumed to be small compared to the total volume [33], the only variable volume in the structure is the volume of the hydraulic cylinder: where V cyl is the volume of a cylinder chamber, A cyl is the piston area and s is the distance between the cylinder end and the piston.This distance can be defined as s = z in piston side and s = z 0 − z in piston rod side, where z is the piston position and z 0 is the length of the cylinder stroke.The total distance between the cylinder attachment points to the surroundings can be defined as l = l 0 + z, where l 0 represents the cylinder length, when the cylinder is at its shortest.If all the leakages in the hydraulic system are ignored, only the valves are allowed to introduce flow to the system.In hydraulic machinery, the flow can be assumed to be turbulent through the valves and therefore proportional to the square root of the pressure loss [26].In this study, a semi-empirical approach is used, where the flow through a valve can be calculated as [9]: where K is a coefficient dependent on fluid properties, discharge coefficient and cross-sectional area of the valve, Q is the flow through the valve, p is the pressure drop in the direction of Q and sign is the signum function.When considering a voltage controlled critically overlapped linear directional control valve, the value of K can be expressed as [13]: where U is the spool position of directional control valve and C v is a semi-empirical flow rate coefficient.
It is important to note that the structure of connected volumes will change when the sign of the spool position changes.The value of C v can be solved from manufacturers data sheet values as: where Q nom , U nom and p nom are flow, valve position and pressure difference in a nominal operational point.
Although the relationship of Eq. ( 12) is commonly used in the literature and can be found in valve data sheets, it has a significant computational drawback of the derivative of flow with respect to pressure loss being infinite at p = 0.This computational drawback can be solved by assuming the pressure loss-flow rate relationship to be linear with small p [12].The assumption of linear pressure drop with small pressure loss can be also justified by physics of the flow, since with the small pressure loss the flow can be assumed to be laminar, which implies linear pressure loss-flow rate behavior [26].To form a continuous pressure loss-flow function, this is realized by linear interpolation between the origin and the pressure loss of which the flow is turbulent, as in [12]: where p lim is the pressure level, where the flow is assumed to be turbulent.One of the numerical advantages of this formulation is its property of finite derivatives at p = 0. On top of the flows producing change to control volume pressures, also the volumetric changes of the control volumes themselves change their pressures as described in Eq. (7).As stated earlier, the only variable sub-volume considered is the cylinder chamber volume, which is determined by Eq. (11).By differentiating Eq. ( 11) with respect to time, the rate of change of a single cylinder chamber volume can be determined as: where ṡ = ±ż depending on the side of the hydraulic cylinder.
The spool position U can be modeled by first-order model: where U in is the valve control signal and τ is time constant of the valve.The scalar equations for single hydraulic system elements, presented above, can be expanded into more general matrix form, which derivation is presented in the following paragraphs.This matrix form is a basis of the later discussed kinematic observer.Consider a system consisting of N vol control volumes, N flow connections between the control volumes, N u control signals, and N cyl hydraulic cylinders.If both inertial and external forces affecting to the hydraulic circuit are unknown, the hydraulic circuit can be modeled with random walk piston accelerations: where X is the state vector, z i , żi and zi ∈ R N cyl ×1 are the vectors of respective cylinder positions, velocities and accelerations, p ∈ R N vol ×1 is vector of pressures of control volumes, U in and U ∈ R N u ×1 are vectors of respective spool positions and control signals and τ ∈ R N u ×N u is a diagonal matrix of valve time constants.The matrix T p,Q ∈ R N vol ×N flow is a path matrix with values −1, 0 and 1 defining the connections of flows.The matrix T V,cyl ∈ R N vol ×2N cyl is a Boolean matrix defining which cylinder chamber is part of which control volume.Additionally, the flows Q are organized in a vector Q ∈ R N flow ×1 , effective bulk modules divided by the corresponding control volume in vector B e,rel ∈ R N vol ×1 and time derivatives of the cylinder chamber volumes in vector Vcyl ∈ R 2N cyl ×1 .Example of a hydraulic system topology is depicted in Fig. 2.
It is important to note the fundamental difference between the previous studies applying state estimation for hydraulically actuated multibody systems [12], [17] and the present study in terms of Eq. ( 18).In the aforementioned studies, acceleration is not a part of the state vector and correspondingly the derivative of velocity is evaluated by means of equation of motion.
On the contrary, Eq. ( 18) utilizes, in part, the idea introduced in [23], where accelerations are included as a part of the state vector and correspondingly the velocity derivatives are not based on equation of motion but instead the values of acceleration state variables.That simplification essentially ignores the prior information of equation of motion, which, if accurately known, could be used to compute more accurate estimates.However, as the purpose of the algorithm is to provide more accurate mass estimates, the assumption is that the equation of motion is not all the time accurate enough to be used to predict acceleration levels.Moreover, using equation of motion (including the payload mass) to predict accelerations and consecutively computing values of the payload mass based on computed Fig. 2 Example of a hydraulic system accelerations would result in a circular argument kind of algorithm, which might show undesirable dynamics and would probably require extensive tuning.

Linearization of the system model
The discrete extended Kalman filter employs the tangent space of the originally nonlinear state-space to evaluate the system response.The nonlinear state-space model described in Eq. ( 18) can be augmented with a null-valued virtual control vector U v being the same size as the state vector [25]: where f aug is the augmented state-space model.It can be shown, that a point U v = −f(X, U in ) is an equilibrium point for the augmented state-space model regardless of the selection of X and U in and the augmented state-space model is linear as a function of U v .For this reason, the system can be analyzed near an equilibrium point by introducing small deviations re f to state vector, control vector and virtual control vector, respectively: where respective reference state vector and reference virtual control vector, X ref and U in,ref , can be arbitrarily selected.Because of the definition of U v being a null vector and the property of U v,re f = −f(X, U in ) in the equilibrium point, the deviation u v with respect to the equilibrium point can be written as: When applying the equation ( 20) to (18), it follows: where where s, ṡ ∈ R 2N cyl ×1 are the respective vectors of cylinder chamber heights and rate of changes and p ∈ R N flow ×1 is vector of pressure differences across the corresponding flows.In following derivations, a convention is followed, where n is referring to the row index and m to column index of an element in a Jacobian matrix.The partial derivatives in matrix A related to piston movement can be presented as: and ∂s/∂z = ∂ ṡ/∂ ż ≡ T s,z is a matrix with values -1, 0 and 1 depending on the definitions of positive directions of piston movement and in which cylinder specific cylinder chamber height s is related to.The partial derivatives related to the term A 44 can be expressed as: where δ is Kronecker delta operator and The Jacobian of flow vector Q with respect to spool position vector U can be calculated as: where an additional condition needs to be applied that It should be noted that Eq. ( 32) imposes a discontinuity to the system as the configuration changes when U changes the sign.As a consequence, the local linearization underestimates the sensitivity of the system near the point U = 0, which leads underestimated uncertainties in the estimation.
Proper scaling of the state-space models is an important aspect for the numerical accuracy of the models.When using SI units, the dimensions of positions, velocities and accelerations in hydraulic systems are typically in the range of 10 −3...1 while the dimensions of pressures are in the range of 10 5...7 .The issue of different scales is overcome in this study by using MPa unit for the pressures and m, m/s and m/s 2 for positions, velocities and accelerations, respectively.

Estimator structure
The estimator is constructed using dual estimation methodology, where the estimation problem is separated into estimation of system kinematics and state of the hydraulic system and this information is further used to estimate the mass of a payload.The aim of this approach is to allow robust estimation of system kinematics and hydraulics and to consequently address observability issues related to mass estimation using a separate algorithm.The overall structure of the estimator is depicted in Fig. 3.The estimation of system kinematics and dynamics of hydraulic circuit is named here as kinematic observer, because it does not use any information of the system inertial or external forces.On the other hand, based on the estimated motion and state of the hydraulic circuit, the mass of the payload can be solved from the equation of motion using inverse dynamics.The mass estimate is further filtered to remove noise that the solution of inverse dynamics includes.Because of the coarse state-transition model of the estimator, the estimator has poor ability to perform the estimation without measurement signal available.Moreover, the numerical stiffness of the hydraulic system may produce an unstable estimator in the case of infrequent and large corrections.The mismatch between the measurement frequency and the estimator frequency is in this study solved by linearly interpolating the measurement signal between consecutive data points.In real-time applications, this creates an additional delay equal to time period between measurements.However, because the size of a measurement time step is considered insignificant compared to the dynamics of the estimation algorithm, this can be considered as acceptable.

Kinematic observer
The state-space model presented in Eq. ( 22) can be augmented with the statistical uncertainties of the system.In practice, this is performed by assuming additional zero mean disturbance signals to be affecting the system states [7]: where W ∈ R N w × 1 is the vector of N w disturbance signals with any length and matrix B 1 maps the disturbance signals to derivatives of state variables.
The disturbance signals are defined to have the respective power spectral densities of R wpsd,1 , R wpsd,2 , ..., R wpsd,N w , respectively, and zero cross-spectral densities.When defining the reference state to be the current state in Eq. ( 24), the process can be discretized using the zero-order hold for the control as [18]: ) where where the matrix G 1 k is the discrete equivalent for B 1 .The a priori estimate Xk+1 can be evaluated using the state-transition model of Eqs.(34) as [7]: The a priori estimate can be corrected using the information of sensors and consequently the posterior estimate Xk+1 can be calculated as [7]: where K k is the Kalman gain, Y k is the vector of measurements at a single time step, and h( Xk+1 ) is the measurement function evaluated in a priori state estimate.The evolution of the Kalman gain is evaluated by following set of equations [29]: where P − k+1 represents the estimate covariance before the measurements, P + k+1 is the estimate covariance after a measurement, P k is plant noise covariance matrix, R v is measurement noise covariance matrix and H is Jacobian of measurement function h() with respect to state vector X evaluated at X = X ref .
The measurement noise covariance matrix and plant noise covariance matrix are the main tuning parameters of the filter.The measurement noise covariance matrix can be often tuned directly using the assumed measurement noise properties of each sensor.The formation of plant noise covariance matrix is not equally easy, because of the dynamical relationships between the states.For this reason, the plant noise covariance matrix is evaluated by means of continuous time model of the system, for which the noise components definitions are more intuitive.
For the evaluation of plant noise covariance, assumptions of plant noise effect to the system is needed.It is reasonable to assume, that position-velocityacceleration relationship is accurately modeled, but the uncertainties are related to accelerations, pressures and valve positions of the system.For this reason, B 1 and R wpsd are defined as: where R wpsd,z ∈ R N cyl ×N cyl , R wpsd, p ∈ R N vol ×N vol and R wpsd,U ∈ R N u ×N u are diagonal matrices with power spectral densities of individual noise components on the diagonal.By using the definition of Eqs. ( 43) and ( 44), the plant noise covariance matrix can be formed as [7]: where The detailed derivation of Eqs.(43-49) is presented in appendix A.
Instead of tuning each element of the plant noise covariance matrix, the tuning can be performed using fewer parameters, which are the power spectral densities of plant noises.

Inverse dynamics
In inverse dynamics, the basis assumption is the motion of the system described by positions, velocities and accelerations to be known [10].In the traditional approach, the external forces of the system are solved based on this known or desired motion.However, in this study, the inverse dynamics is used to solve the unknown mass property of the system.When the pressures acting on the hydraulic cylinders are known, the vector of hydraulic forces can be solved as: where • is element-wise vector multiplication and A cyl is vector of piston areas.The force produced by the pressures is corrected with the friction force of the cylinder to get the effective force the hydraulic actuator is applying on the system.For simplicity, the friction force is assumed to be a combination of pressuredependent terms and velocity-dependent terms as in [1].For velocity-dependent terms, Brown-McPhee friction model is applied as it includes static, dynamic and viscous friction [2] and it has been earlier applied to describe the cylinder friction [13].The computational benefit of Brown-McPhee model is that it does not involve any discontinuities [2].Here, a variant of Brown-McPhee model is being used, where the normal force is assumed constant and sufficient to enable viscous friction [13].Consequently, the friction force can be expressed as [13]: where F μB−M is the Brown-McPhee friction force, F s is the static friction force, F c is the Coulomb friction force, v s is the transition velocity and σ is the viscous friction coefficient.To consider the pressure dependency of the friction force, Brown-McPhee friction model can be enhanced with pressure-dependent terms.
As in [1] the pressure-dependent friction is assumed to be linearly dependent on both the pressure difference between cylinder chambers and the piston rod side pressure: where F μp is the pressure-dependent friction component, p piston and p rod are the respective pressures on piston side and piston rod side, and C μp1 and C μp2 are the friction coefficients related to pressure-dependent friction force component.The total friction force can be calculated as a combination of Eqs.(51, 52) as: Consequently, the effective actuator force is: where F μ consists of respective scalar components produced by Eq. ( 53).
When the independent coordinates are defined to be the piston positions, the actuator force vector is acting in the independent joint space.Therefore, the equation of motion (Eq. 3) can be written with composite external force vector without the actuator forces and the actuator forces can be summed separately.The residual for such written equation of motion can be written as: Assuming that the number of degrees of freedom of the system is the same as the number of unknown masses, and both the mass properties and degrees of freedom are correctly defined, an exact solution for Eq.(55) can be found.It is important to note that if the number of degrees of freedom is larger than number of masses, Eq. (55) becomes optimization problem for finding min{|g|}.Obviously, in specific points of general case, the solution may not be found, such as if 1. zi j = 0 and Q j = 0 or 2. M zi j = Q j , where j is referring to a respective row of a column vector.However, nonexistence of the solution can be considered acceptable point-wise, when using proper filtering to limit the effect of resultant solution outliers.Considering a case, where the number of unknown masses matches with the number of degrees of freedom of the system, Eq. ( 55) can be solved by applying Newton-Raphson iteration as: where h is the index of iteration algorithm.The final value from the previous time step is always used as an initial value of the iteration.In this study, the Jacobian matrix of g with respect to m is calculated numerically using forward difference: where the value for is selected to be 1 × 10 −2 kg.

Mass update algorithm
When the evolution of load mass as a function of time is modeled using random walk model, both the statetransition matrix F and the measurement matrix H become identity matrices.For this kind of system, the application of Kalman filter would produce an estimator that, after a transient behavior at start, becomes equal to first-order filter, when the Kalman gain K becomes constant.For this reason, the mass is estimated using first-order low-pass filter: where K f is a filter coefficient, m l is the value of estimated load mass produced by inverse dynamics and m l,L P is its respective value after low-pass filtering.To maximize the filter convergence rate, low-pass filter is combined with additional weighted moving average filter with a window size N avg .The weight distribution of the filter is chosen to be symmetrical triangular, where a new value enters the filter window with close to zero weight and correspondingly the oldest value exists the window with the same weight as: where y (k+1) ≡ m l,W M A,(k+1) is the output of the filter, x ≡ m l,L P is the input of the filter and the window size is equal to The true nature of static friction force is that it behaves like sign function as a function of velocity.As a consequence, the true value of the friction force is difficult to determine with zero velocity.This leads to biased estimates as the friction force can resist movement of the cylinder in static equilibrium different to what is indicated by Eq.( 53).For this reason, the value of filtered mass estimate is not updated when the piston velocity is below the transition velocity v s .

Numerical example and experimental setup
A simple hydraulic crane, illustrated in Fig. 4, was studied as numerical example.The only unknown parameter of the system is the mass of the payload the crane is lifting.The proposed estimation algorithm was studied with two methods, which differ by the origin of the measurement data.In the first method, the data were created artificially by simulating the system and augmenting the simulation results with additional measurement noise.The original simulation results without noise are here referred to as precise values, the same results with additive white noise as artificial measurements and the simulation model as artificial plant.
In the second method, an experimental setup was used producing real measurement data for the estimation algorithm.These measurement data are here referred to as real measurements and the measurement setup as experimental plant.In the case of artificial plant, the control signal of the hydraulic control valve and the initial condition of the system were taken from the experimental plant to achieve results that are mutually comparable.The simulation model and experimental setup are depicted in Figs. 4 and 5, respectively.The hydraulic system in Fig. 4 consists of a hydraulic cylinder, two control volumes V 1 and V 2 , a voltage controlled directional control valve and external pressures p P (pump pressure) and p T (tank pressure).The measurands of the system were the piston position and pressures at V 1 and V 2 .The frequency of the measurements was 100 Hz.
The mechanical structure is modeled using four bodies: hydraulic cylinder, hydraulic piston, arm of the crane and the payload.The hydraulic cylinder is assumed to be massless, while the mass of the hydraulic piston is considered.It should be noted that this is equivalent to modeling the cylinder and piston as a composite body attached to the pillar with a pin-inslot type of joint.The open loop structure of bodies is built using the path (cylinder-piston-arm-load mass).Because the load mass is rigidly attached to the lift arm, it can be modeled using prismatic joint and an additional constraint equation defining the related coordinate to be always zero.The other cut joint of the system is located at the joint between the pillar and lift arm, and a corresponding constraint equation is added to form the closed-loop structure.
When producing the artificial measurements, the pressure-dependent component of the friction model (Eq.52) is being smoothed by replacing the sign(ż) function with tanh(4ż/v s ) to avoid numerical problems.The friction model for the system is tuned using a training measurement data set from which the cylinder velocity and the hydraulic force are estimated using the kinematic estimator.The fit of the friction model to the hydraulic force data is depicted in Fig. 6 where the colors represent the number of data points per interval in logarithmic scale.Also, the median friction force, calculated from data points, is depicted in Fig. 6.The parameters of the system are described in Table 1.

Results and discussion
For the mass estimation algorithm performance, the most important indicator is the accuracy and convergence rate of the estimated mass.For this reason, both the filtered mass value and its original value are analyzed.Another important aspect is to understand the performance of the kinematic observer, as it is used as an input for the calculation of inverse dynamics.From the system with artificial measurements, estimated and true values of the states are presented and the measured values are presented when applicable.In the case of experimental setup, the information of real values of the states cannot be accessed as in the case of artificial measurements and for this reason only measurement data and estimates are presented.Moreover, for the estimates produced by the kinematic observer, the confidence intervals are presented in the figures with 95% confidence bounds, which correspond to ±1.96σ and are calculated from the diagonal terms of P + k shown in Eq. (41).

Work cycle presentation
Both the experimental plant and artificial plant have the same control signal applied to the directional control valve.Also, the initial states are assumed to be the same and are based on the first samples of the measurement data set.The control signal U in is depicted in Fig. 7a with the corresponding simulated valve position.It is important to note that here ± 1 represents fully open valve position, where positive value represents Pump-V 1 , Tank-V 2 connection and correspondingly negative value represents Pump-V 2 , Tank-V 1 connection.The resultant values of cylinder position from both the experimental plant and artificial plant are presented in Fig. 7b.As b shows, the motion of the simulation model is close to the experimental setup, but due to modeling inaccuracies they are diverging from each other.

Results with artificial measurements
Figure 8 shows the estimated, measured and precise values of piston position for the artificial plant.In the larger scale, all these values seem to be coinciding.However, the detailed view shows that a small difference is present between the estimated and precise value and the measured value differs from the precise value by a zero mean noise.Moreover, the measurement noise prevents recognizing any fast velocity dynamics directly from the measurement data of piston position.The estimation confidence interval can be considered narrow while the confidence interval does not show any growing trend, which indicate both good observability properties and relatively low process noise effects assigned on the position level in estimator tuning.In Fig. 9, the estimated and precise values of piston velocities of artificial plant are depicted.Here, the velocity estimate follows accurately its precise value and, as the detailed view shows, is able to track the rapid fluctuations having frequency of around 10 Hz.The confidence interval on velocity level is already broader than on position level (Fig. 8), which can be partly explained by the lack of direct velocity measurements.However, the confidence interval can be still considered relatively narrow and not showing significant trend of growing, which indicates the assumed stochastic sys- tem having good observability properties also on the velocity level.This can be explained, in part, by the coupling between hydraulic pressures and piston velocity and availability of those pressure measurements for the estimation algorithm.What is also notable here is the similarity of the shape of the velocity curve to earlier presented spool position curve (Fig. 7a).This behavior is due to the friction force (velocity-dependent component) and gravity force (static component) dominating the force balance over inertial forces (acceleration-dependent component).
Figure 10 shows the precise and estimated values of acceleration of the piston in artificial measurement configuration.In this application, the accurate estimation of acceleration is challenging due to lack of acceleration measurement but also due to a priori assumption of zero time derivative.As can be seen from Fig. 10, the estimated piston acceleration follows the precise piston acceleration even in the fluctuations of around 10 Hz in frequency.However, the acceleration data show clear phase shift between the precise and estimated value of piston acceleration.The main reason for this phase shift is the lack of a priori knowledge of the acceleration behavior and as such inability to make reliable predictions of it.Consequently, the acceleration estimate is The lack of a priori knowledge of acceleration is also visible in confidence interval of acceleration estimate.The confidence interval is significantly broader than in the case of position and velocity estimates.However, to have acceleration estimate with high bandwidth, the process noise for the acceleration, modeled with random walk, needs to be assumed relatively high, which automatically results in broad confidence intervals.Moreover, the confidence interval shows relatively stable, non-growing behavior, which indicates acceptable observability properties at the acceleration level.
Estimated, measured and precise values of pressures p 1 and p 2 with artificial measurements are shown in Fig. 11a and 11b, respectively.The both pressure estimates show good agreement with the corresponding precise values and the phase shift can also be considered small.The confidence intervals for both of the pressure estimates show non-growing behavior, which indicate their good observability properties.The good observability properties of the pressures were also expected, because the pressure measurements are directly provided for the Kalman filter algorithm.
Based on the kinematic estimates, the payload mass is estimated using inverse dynamics as described in Sect.3.2.Both the resultant mass estimate from inverse dynamics and its value after low-pass filtering and moving weighted average filtering are presented in Fig. 12.Because the mass estimate before filtering is purely based on the solution of equations of motion, it convergences instantly to the actual value of the load mass.However, the unfiltered estimate includes clear outliers.A closer look to velocity estimate shows that these peaks occur specifically in situations where the abso-Fig.11 Hydraulic pressures with artificial measurements Fig. 12 Mass estimation with artificial measurements, m load = 50 kg lute value of velocity rapidly decreases to zero making the velocity estimate to overshoot.Because the friction force is sensitive for small velocity changes when the velocity magnitude is small, the aforementioned behavior creates error in the mass estimation.However, the effect of wrong estimates with small velocity magnitudes is here solved by filtering the estimate and also by ignoring the new mass data in the low-velocity region.

Results with real measurements
Figure 13 shows both the estimated and measured piston positions in the experimental setup.Similarly, as in the case of artificial measurements, the position esti-Fig.13 Piston position z in experimental setup mate follows the measured value accurately.However, the difference between the measured and estimated values is observed to be larger than in the case of artificial measurements.The confidence interval for the position estimate shows similar non-growing behavior as in the case of artificial measurements.
For the velocity estimate (Fig. 14a) and acceleration estimate (Fig. 14b), a formal verification could not be performed with the experimental setup.However, the estimated velocity and acceleration levels are similar to the corresponding simulation results, which indicates successful estimation.
Similarly as in the case of artificial measurements, the confidence intervals grow toward the higher time derivatives of position (velocity and acceleration).However, in the both cases the confidence do not show significant trend of growing, which indicates acceptable observability properties.
The estimated and measured values of pressures p 1 and p 2 are shown in Figs.15a and b, respectively.The pressure graphs show higher pressure fluctuations compared to the results from the artificial plant.The higher fluctuation level may be a consequence of system flexibility which is not included in the artificial plant.Moreover, the pressure estimates are not able to follow the measured pressures as accurately as in the case of artificial measurements.The reason for this is probably the higher level of fluctuations, which makes estimation more difficult.As in the case of artificial measurements, the directly available information from pressure measurements prevents the confidence intervals of estimated pressures from growing.
In Fig. 16 the values of mass produced by inverse dynamics and its corresponding filtered value are presented for the experimental measurement configuration.As the graphs show, the unfiltered estimate is biased in the beginning and at the end of the data sam- ple when the mechanism is not moving.Moreover, clear outliers exist in the unfiltered estimate, and those outliers occur mainly in the low-velocity region.How-Fig.16 Mass estimation with real measurements, m load = 50 kg ever, because the filtered value is not updated during the periods of low velocities, the filtered estimate does not suffer as severely for the friction force bias.
Compared to earlier presented estimation results with artificial measurements, the mass estimate with real measurements fluctuates significantly more.The reason for increased fluctuations might be in the performance of the kinematic estimator, but also in the accuracy of inverse dynamics.The most probable reason in inverse dynamics can be considered to be the rigid body assumption, which is not valid for the lifting arm, that bends significantly in the experimental setup.

Numerical efficiency
The numerical efficiency of algorithm is analyzed in MATLAB implementation using watchdog timer (tictoc function).The parts kinematic observer and inverse dynamics combined with post-filtering are computed separately, and the total calculation time is divided by the length of the time series to get the ratio between computational time and real time.The kinematic observer used in average 0.84 s computational time for 23.5 s data which means ratio of 3.6 % between computational time and real time.Correspondingly, inverse dynamics and post-filtering consumed in average 5.3 s computational time for 23.5 s data making the ratio between computational time and real time 22.6 %.Consequently, the total ratio between computational time and real time is 26 % meaning the algorithm is clearly real time capable with current implementation and level of complexity of the structure.

Conclusion
This work proposes an estimation algorithm for evaluating the amount of payload mass attached to a hydraulically driven crane.The estimation was divided in two major parts: 1. estimation of the system state based on kinematics of the mechanism and dynamics of the hydraulic system and 2. solving the unknown mass by substituting the previously estimated states to the equation of motion.The kinematic observer was implemented using a novel approach incorporating the pressure dynamics of the hydraulic circuit to the statetransition model.
The results show that when using artificial measurements the mean value of the unfiltered mass estimate converges rapidly to the actual value of mass.However, the unfiltered estimate includes a few outliers which occur at the time instances when the velocity is low.For this reason, the mass value received from inverse dynamics was filtered using low-pass filter and moving weighted average filter and the fluctuations in the resultant estimate were significantly reduced.The experimental tests show similar results.However, fluctuations in the mass estimate are higher, which is assumed to be due to the assumption of undeformable bodies.
As indicated by the friction parameters, the friction has a significant role in the evaluation of total actuator force.Moreover, for the accurate evaluation of the friction force, the velocities need to be known.The results from artificial measurements show that the kinematic observer is able to evaluate the actuator velocity accurately.Moreover, the velocity estimate with real measurements has similar values to the respective estimate with artificial measurements, which indicates successful estimation of velocity also with the experimental setup.
In future work, the proposed method could be applied in more complicated systems, such as excavators or loaders to investigate the usability of the method under real work cycle conditions.Moreover, incorporating the structure flexibility to equation of motion could decrease the fluctuations in the unfiltered estimate.Consequently, the low-pass filter bandwidth could be increased leading to faster convergence.
Funding Open Access funding provided by LUT University (previously Lappeenranta University of Technology (LUT)).This study is, in part, supported by the Academy of Finland (project # 316106) and by Business Finland within the project Next Generation Digital Twin Solutions -XGen-DigiTwin.

Availability of data and material
The datasets generated during and/or analyzed during the current study are not publicly available but are available from the corresponding author on reasonable request.

Conflicts of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Derivation of plant noise covariance matrix
The kinematic observer used in this work requires knowledge of the statistical properties of discrete time plant noise and measurement noise, where the plant noise covariance components can be considered difficult to tune manually.However, the plant noise properties may be easier to obtain in continuous time frame by the spectral properties of the signal and convert the spectral properties of continuous time signals to statistical properties of the discrete time signals.In this study, this conversion is performed using Eq. ( 45), where the state-transition matrix F as well as matrices B 1 and R wpsd defining the noise influence to system is needed.When using the approximation order of n = 0 in zeroorder hold discretization (Eq.35), the state-transition matrix can be written as: As described in Eq. ( 43), matrix B 1 can be defined as: where the noises are assumed to be applied on derivatives of accelerations, hydraulic pressures and spool positions.
The power spectral density matrix for the noise signals can be formed as: where a diagonal element represents power spectral density of a noise signal and an off-diagonal element represents cross-spectral density between two noise signals.Because of the assumption of independent noise signals, all the off-diagonal elements are zeros.
When performing the matrix multiplication of Eq. (45), it results: where

Fig. 4 Fig. 5
Fig. 4 Construction of the simulation model of the simple hydraulic crane

Fig. 10
Fig. 10 Piston acceleration z with artificial measurements purely based on observations, which creates the phase shift between the precise and estimated value.The lack of a priori knowledge of acceleration is also visible in confidence interval of acceleration estimate.The confidence interval is significantly broader than in the case of position and velocity estimates.However, to have acceleration estimate with high bandwidth, the process noise for the acceleration, modeled with random walk, needs to be assumed relatively high, which automatically results in broad confidence intervals.Moreover, the confidence interval shows relatively stable, non-growing behavior, which indicates acceptable observability properties at the acceleration level.Estimated, measured and precise values of pressures p 1 and p 2 with artificial measurements are shown in Fig.11a and 11b, respectively.The both pressure estimates show good agreement with the corresponding precise values and the phase shift can also be considered small.The confidence intervals for both of the pressure estimates show non-growing behavior, which indicate their good observability properties.The good observability properties of the pressures were also expected, because the pressure measurements are directly provided for the Kalman filter algorithm.Based on the kinematic estimates, the payload mass is estimated using inverse dynamics as described in Sect.3.2.Both the resultant mass estimate from inverse dynamics and its value after low-pass filtering and moving weighted average filtering are presented in Fig.12.Because the mass estimate before filtering is purely based on the solution of equations of motion, it convergences instantly to the actual value of the load mass.However, the unfiltered estimate includes clear outliers.A closer look to velocity estimate shows that these peaks occur specifically in situations where the abso-

Fig. 14 Fig. 15
Fig. 14 Estimated piston velocity and acceleration in experimental setup

Table 1
Parameters used in the numerical example