Payload estimation of weight bearing machinery using multiple model adaptive estimator system and method

Abstract

A method of payload estimation of weight bearing machinery, the method comprising the steps of: (a) creating a series of Kalman filter approximations to the weight bearing machinery for different levels of payload; (b) determining from the series an approximation to the current operating characteristics of the weight bearing machinery; and, (c) utilizing the parameters derived from the step (b) to determine an estimated payload of the weight bearing machinery.

Description

CROSS REFERENCE TO RELATED APPLICATIONS

This application is the U.S. national phase of PCT/AU2007/000515 filed Apr. 19, 2007. PCT/AU2007/000515 claims benefit under the Paris Convention to AU 2006902069 filed Apr. 20, 2006. The disclosures of both of AU 2006902069 and PCT/AU2007/000515 are hereby incorporated herein by reference.

FIELD OF THE INVENTION

The present invention relates to a method and system for estimating the payload of a mining shovel and, in particular, discloses the utilisation of Kalman filtering in the estimation process.

BACKGROUND OF THE INVENTION

Electric mining shovels (EMS) are large electro-mechanical excavators commonly used in open-cut mining to load haul trucks. They are critical production units at most open-cut mine sites and there is an ongoing need to improve their productivity. FIG. 1 illustrates schematically one of these shovels.

It has been recognized since the 1950s that the capability to measure the weight of material collected by a mining shovel (or other excavator) at each dig can enhance safety and facilitate productivity improvement. The present invention is directed to ideas and methods from the field of optimal state and parameter estimation to address the problem of payload weighing for electric mining shovels.

One example of how knowledge of payload can be used to improve safety and productivity is by reducing the frequency of overloading of haul trucks. Truck overloading comes about primarily because mine sites perceive productivity benefits in loading the maximum possible amount of material into each truck. Most mine sites continue to favor the intuitive argument that pushing the limits of rated capacity increases operational efficiency. Moreover, they perceive the benefits are sufficiently great to pursue this objective aggressively. It is, for example, common practice to award operators performance bonuses if sufficient trucks are loaded within an ‘optimal’ band close to the truck capacities.

Truck overloading reduces safety, reduces the life of the truck, and increases down-time and the maintenance costs. To discourage overloading, maintenance and warranty contracts between a mine operator and a truck manufacturer, nowadays, will usually include a penalty clause that makes it mandatory to stop and dump the load where the truck stands if the scales on-board the truck show it to be overloaded. A dumped load results in a direct productivity loss that is compounded by the costs associated with having to rework the dumped material.

Modern trucks are fitted with on-board scales that are capable of monitoring the weight of material in their trays. These systems typically use suspension strut pressure strategies, and can provide weights accurate to plus or minus 3-5% when calibrated. It is generally accepted, however, that these on-board truck scales become accurate only after the truck is in motion, typically when second gear is reached on flat road. The weights reported at the loading point are generally held to be accurate only to plus or minus 20%, due to the influence of mechanical infelicities such as strut seal friction and strut rod bending in combination with other factors which cause the suspension to lock, such as rock spillage around the loading area and uneven distribution of material in the truck tray.

It normally takes four passes of a large modern electric shovel (with a nominal payload of approximately 80 t) to load a haul truck with a 300 t capacity. To manage the loading process, operators currently use information provided by truck scales (which are known to give poor estimates at the load point) in combination with their judgement. The main motivation for the development of on-shovel payload weighing is the belief that the risk of overloading can be better managed if the shovel operator knows, after each digging pass, how much material is in the bucket, together with the remaining capacity, of the truck tray. This information should empower making a ‘go’ or ‘no go’ decision prior to the final dig pass being dumped into the truck.

Other benefits of payload monitoring are envisaged. For example, giving operators real-time payload weights provides a source of feedback that can help them achieve better, more consistent, operation of the machine. There is, moreover, potential value in being able to monitor shovel performance on a cycle-by-cycle basis; the amount of material collected at each pass is widely held as a key performance indicator.

There have been several previous attempts to develop mining shovel payload systems and several commercial systems are currently available. While manufacturers typically quote accuracy to ±2%, the industry widely accepts that none of these systems can reliably achieve this. For a shovel with a nominal 80 t payload, worst case errors are likely to be as large as ±20 t.

Background work relevant to the payload estimation problem includes

(i) previously developed methods specifically applicable to electric mining shovels,

(ii) methods for estimating payload on other mining equipment such as hydraulic excavators and draglines or cognate equipment such as construction cranies, and

(iii) methods proposed in the robotics literature for the identification of the inertial parameters of robotic links and payload.

Methods for determining payloads on mining machines go back to the 1950s. These systems sometimes involved ingenious means for arriving at an estimate of payload, including analog computations using electrical and fluidic circuits. To the inventor's knowledge none of these systems proved successful; certainly none are in current wide spread use today. Most of the approaches in these early patents have been made redundant by the availability, circa 1980, of microprocessors that could be packaged for deployment on mining equipment. A starting assumption for this work is that any future practical system for payload estimation will be based off microprocessor technology and, accordingly, this review is limited to methods intended for computer-based implementation.

Payload Estimation for Electric Mining Shovels

Chang et al. in U.S. Pat. No. 6,225,574 use a static moment balance on the dipper handle assembly to compute payload. The method equates the moment applied to the dipper and handle assembly by gravity to that applied by the hoist ropes. By careful choice of coordinate frames, the equations are formed so as to be functions of dipper position and hoist rope force only, and are independent of the force applied by the crowd drive. This algorithm is known to be the basis for several commercial payload systems.

The hoist rope force is obtained by measuring the armature current and multiplying it by the motor torque constant and transmission ratios. The gravitation force is assumed to act through a known point (the centre-of-mass of the dipper and handle assembly) allowing the determination of an effective force. This force is then divided by the gravitational acceleration constant to give payload mass. Several estimates of payload can be made during the swing phase of the shovel cycle using this approach and these are averaged to determine the final payload estimate.

Chang et al. report that this method performs adequately when the shovel is stationary but gives poor results when in motion. They attempt to overcome this limitation by assigning a ‘fuzzy’ confidence factor between zero and one that weights each payload estimate in the overall sum. The calculation of these confidence factors is based on the observed speed and acceleration of the dipper at the time of measurement. The higher the speed or acceleration, the lower the confidence value; the lower the velocity and acceleration, the higher the confidence value. Confidence limits are used to cull outliers before the final payload estimate is computed as a weighted average of the instantaneous payload estimates.

This weighted-average approach appears to be an ad hoc attempt to deal with the forces resulting from motion of the machine. A mining shovel is a dynamic machine for which the motion on each drive is intermittent and predictable only in its gross character. An open question is the extent to which inertial forces associated with the acceleration and velocity of the dipper influence estimates of payload mass. This includes the centripetal forces acting when the shovel is swinging.

The method of Chang et al. also assumes that the payload centre of mass is invariant and is fixed at the centre of mass of the bucket for all payload sizes. The validity of this assumption is not tested and the sensitivity of payload estimates to errors in payload centre-of-mass is not reported. Moreover, no explicit attempt is made to compensate for drive friction or other sources of loss in arriving at estimates.

The payload estimation method described by Radomilovich in U.S. Pat. No. 4,677,579 uses a dynamic model that accounts for conservative and non-conservative effects in the hoist and crowd systems. The description of the approach includes the discussion of corrections for the inertial forces associated with the rotational inertia of the drive motors and the reduction gear train, the stretch of the hoist ropes, and the inertia of the bucket. Non-conservative losses due to friction and motor inefficiencies are also mentioned. Radomilovich uses a force balance similar to that of Chang et al., from which payload is estimated. The algorithm, as described, ignores several potentially important factors:

• • Swing motions are not included. Centrifugal and gyroscopic forces acting on the system are therefore neglected.
  • The algorithm requires motor accelerations, which are obtained by differentiating motor speed. Motor speed is determined from measurements of armature voltage and current using the steady-state, electrical, motor equations. The differentiation process amplifies noise in the speed signals. The difficulty of obtaining low noise acceleration signals argues for a methodology that does not require inferred measurement of the acceleration.
  • Motor position is calculated by integrating motor speed. As the armature voltage and current are noisy signals, the variance of the inferred position will grow unboundedly with time. In practice this method would require the shovel kinematics to be constantly recalibrated. This limitation can, of course, be overcome by direct position sensing of the motors.

The approach of Blair et al. in U.S. Pat. No. 4,809,794 is based on the observation that, in principle, knowing the co-ordinates of the centre of mass of the dipper allows determination of the strain that the weight of the dipper and its contents induce in any part of the boom or the support structure. Conversely, knowing the strain in a selected part of the support structure and the location of the center-of-mass of the dipper allows the weight of the dipper and its contents to be calculated. They propose the fitting of strain gauges to the support structure (on the A-frame structure that supports the boom) together with the use of experimentally or analytically determined influence coefficients, that are functions of dipper position, to relate measured strain to weight of the dipper, handle, and payload. A large part of this patent is given over to describing logic to establish when the shovel is swinging to the dump location, and there is comparatively little description of how the influence coefficients relating weight to strain might be determined. The algorithm averages the weights obtained at each sample over some part of this swing phase. The approach assumes a location for the payload mass-center and does not account for the dynamic loads associated with motion. The authors note that it is preferable to exclude from the averaging step, a number of the weights calculated at the start and end of the swing when these dynamic loads are often at their highest.

SUMMARY OF THE INVENTION

In accordance with a first aspect of the present invention, there is provided a method of payload estimation of weight bearing machinery, the method comprising the steps of: (a) creating a series of Kalman filter approximations to the weight bearing machinery for different levels of payload; (b) determining from the series an approximation to the current operating characteristics of the digger; (c) utilising the parameters derived from the step (b) determine an estimated payload of the weight bearing machinery.

The weight bearing machinery can comprise an electric mining shovel. The inputs to the Kalman filter can include weight bearing machinery drive control signals for controlling movement of the machinery. The members of the series of Kalman filter approximations simulate increasing payloads on the weight bearing machinery. The step (b) further comprising the steps of: forming a current model probability for each model and then forming a weighted sum of the models to determine an overall model.

BRIEF DESCRIPTION OF THE DRAWINGS

A preferred embodiment will now be described with reference to the accompanying drawings in which:

FIG. 1 illustrates schematically an electric mining shovel;

FIG. 2 is a plot of the hoisted drive signals measured in an electric mining shovel during a typical swing to dump;

FIG. 3 illustrates a plot of the crowd drive signals of an electric mining shovel during a typical swing to dump;

FIG. 4 illustrates a plot of the swing and drive signals of an electric mining shovel during a difficult swing to dump;

FIG. 5 illustrates the shovel configuration at three second time intervals during the swing-to-dump phase of the machine cycle;

FIG. 6 illustrates a block diagram structure of the Kalman filter utilised in the preferred embodiments;

FIG. 7 illustrates an overview of the multiple model algorithm utilised in the preferred embodiment;

FIG. 8 illustrates the steps involved in the preferred embodiment;

FIG. 9 illustrates coordinate systems utilised in the mining shovel geometry;

FIG. 10 illustrates a further coordinate system utilised in the mining shovel system;

DESCRIPTION OF PREFERRED AND OTHER EMBODIMENTS

The preferred embodiment starts from the view that dynamic loads cannot be a priori ignored. Hence we seek to include them. The rigid body dynamics of a mining shovel take the general formM(θ){umlaut over (θ)}+f(θ,{dot over (θ)})=τ  (1)where θ is a vector of configuration variables, e.g. joint angles; {dot over (θ)} is a vector of joint velocities; {umlaut over (θ)} is a vector of joint accelerations; M(θ) is the mass matrix referred to (and dependent on) the relative configuration variables; f (θ,{dot over (θ)}) is the vector of velocity and gravity dependent generalized forces; and τ is a vector of generalized forces applied by the actuators.

A fundamental question is how to incorporate a model such as this into the problem formulation. In reviewing prior art, criticism was made of formulations requiring measurements of positions, velocities, and accelerations of the machine motions, on the grounds that the simultaneous measurement of all three variables for each degree of freedom is, at best, complex and may be impractical. A typical mining shovel will have joint resolvers providing position measurements and could be fitted with tachogenerators to provide velocity measurements, but rotational accelerometers that could be fitted to the motor shaft are not widely available and accelerometers fitted externally, for example to the dipper, would be difficult to maintain. It is reasonable to expect some compliance in the transmission and it is worth noting that, irrespective of where the accelerometers are placed, they would be non-collocated with respect to some system inertias.

It is common practice to overcome the limitation of not having direct velocity and acceleration measurements by numerically differentiating and double differentiating position sensors to get velocity or acceleration, or by integrating and differentiating velocity signals to get position and acceleration. Such schemes are often implemented by high-pass (differentiating) or low-pass (integrating) filters applied to the signals. This approach can become problematic if noise is present on the base signal. This noise is amplified by differentiation and the variance of an integrated signal grows unboundedly.

An alternative approach is to use the causal (or forward) form of the dynamics to predict future motions based on applied forces. This immediately casts the problem into a form suitable for application of optimal state and parameter estimation methods. The advantages of formulating the problem this way are (i) that acceleration measurements are not needed, (ii) that measurement noise and model uncertainty are included naturally within the formulation, and (iii) that the causal structure of the model attenuates noise on driving signals, i.e. ∥M⁻¹∥ is small. The major disadvantage of this approach is that the model retains an ordinary differential equation structure and so the solution procedure is more complicated than that required for the algebraic extraction of payload estimates from the equations of motion.

An optimal state and parameter estimator seeks the best estimates (in a minimum variance sense) of the system dynamic state and model parameters. It provides an implicit, but nevertheless, clean and effective framework for accounting for the inertial loads associated with machine motion.

The preferred embodiment uses the Kalman filter to determine a weight estimate. The Kalman filter is a recursive linear estimator that calculates a minimum variance estimate of the state of a dynamic system using observations that are linearly related to this state. The Kalman filter is optimal in the sense that it minimizes the mean squared estimation error under the assumptions that the state being estimated is that of a linear system driven by Gaussian white noise from which periodic measurements, corrupted by Gaussian white noise, are made.

The Process Model

The process model is used by Kalman filter based method to predict the future state of a system from an estimate of the current state and knowledge of the control inputs. It encodes a priori knowledge about the dynamics of the system under study, and, by enabling the prediction of future system state and state covariance from current estimates of state and state covariance, it propagates information relating to the system state forward in time. In doing so it provides the mechanism for temporal fusion of new information in the form of sensor measurements with knowledge of the current knowledge of the system state.

Kalman filter theory assumes that a linear, stochastic, difference equation can be used to model the system and that measurements made from the system are linearly related to the system states.

The equations describing the system model and the measurements take the formx(t_i)=Φ(t_i,t_i−1)x(t_i−1)+Γ(t_i−1)u(t_i−1)+G(t_i−1)w(t_i−1),  (2)z(t_i)=H(t_i)x(t_i)+v(t_i),t_i=Δt,  (3)where x(t_i) is the n-dimensional state vector; u(t_i) is the l-dimensional vector describing the system inputs; and w(t_i) is the s-dimensional process noise vector representing uncertainty in the inputs and uncertainty resulting from unmodelled dynamics. The matrix Φ(t_i, t_i−1) describes the system dynamics and matrices Γ(t_i−1) and G(t_i−1) describe the distribution of inputs and process noise to the states. The variable z(t_i) is the m-dimensional measurement vector; v(t_i) is the m-dimensional vector of measurement noise; and H(t_i) is the measurement matrix describing how the measurements relate to the states. The sample period, or time between consecutive measurements, is Δt.

The process noise, w(t_i), and the measurement noise, v(t_i), are assumed to be independent zero-mean, white Gaussian noise sequences,E{w(t_i)}=0,  (4)E{v(t_i)}=0,  (5)E{w(t_i)w(t_j)}=Q(t_i)δ_ij,  (6)E{v(t_i)v(t_j)}=R(t_i)δ_ij,  (7)E{w(t_i)v(t_j)}=0.  (8)

The process noise accounts for discrepancies between the process model and the actual behaviour of the plant. The process noise has three components: (i) disturbances that act on the control inputs including measurement errors on the inputs and unmodelled changes in the inputs which are usually assumed to be constant over each prediction step; (ii) uncertainty in model parameters due to both errors in specified values or unmodelled changes in parameters with time; and (iii) stabilizing noise that compensate for plant dynamics not modelled by the process model. It is worth noting that stabilizing noises may not have a well-defined physical basis, but are required to ensure successful filter operation. The process noise covariance, Q(t_i), is positive semidefinite.

The measurement noise describes uncertainties associated with making measurements from sensors. The measurement noise covariance, R(t_i), is positive definite for all i.

The initial condition of the system state is modelled as a Gaussian random vector with known mean and covariance:E{x(t₀)}={circumflex over (x)}₀,E{[x(t₀)−{circumflex over (x)}₀][x(t₀)−{circumflex over (x)}₀]^T}=P(t₀).Derivation of the Kalman Filter

There are several known approaches to deriving the Kalman filter including minimum variance, maximum likelihood, and maximum a posteriori (MAP) approaches. The preferred embodiment uses the MAP formulation because it leads most naturally to the multiple model parameter estimation algorithm used for payload estimation.

Summary of Equations

The optimal state estimate and its covariance are predicted from time t_i−1 to the current time t_i using{circumflex over (x)}(t_i⁻)=Φ(t_i,t_i−1){circumflex over (x)}(t_i−1⁺)+Γu(t_i−1),  (9)P(t_i⁻)=Φ(t_i,t_i−1)P(t_i−1⁺)Φ(t₁,t_i−1)^T+G(t_i)QG(t_i)^T.  (10)

The purpose of the prediction is that it permits the temporal fusion of data.

At time t_i the measurement z(t_i) is incorporated into the state estimate giving the following equations for the updated state estimate and its covariance,A(t_i)=H(t_i)P(t_i⁻)H^T(t_i)+R(t_i),  (11)K(t_i)=P(t_i⁻)H^T(t_i)A(t_i)⁻¹,  (12)r(t_i)=z(t_i)−H(t_i){circumflex over (x)}(t_i⁻),  (13){circumflex over (x)}(t_i⁺)={circumflex over (x)}(t_i⁻)+K(t_i)r(t_i),  (14)P(t_i⁺)=P(t_i⁻)−K(t_i)H(t_i)P(t_i⁻).  (15)

FIG. 6 shows the block diagram structure of the above formulation of a Kalman filter.

The Extended Kalman Filter

The Kalman filter provides the optimal state estimate of a discrete-time, linear system, with discrete, linear measurements. However, all systems are, to some extent, non-linear and in many the non-linearities are significant enough that the linear Kalman filter cannot be applied. Electric mining shovels fall into this category because the equations of motion that govern them are dependent on configuration and velocity. While these dynamics are non-linear, they are smooth in the sense that they can be everywhere differentiated.

The usual approach for solving state estimation problems on non-linear systems with “smooth” dynamics is the extended Kalman filter or EKF. The EKF essentially involves the application a linear methodology to a non-linear problem by linearizing the non-linear state and measurement equations about the state trajectory. The EKF is popular because it is conceptually simple and, for a given non-linear system, it provides best linear estimator measured by minimum mean squared error, and so, in this sense, it is optimal.

The Kalman filter equations given above apply to discrete time system. Where the non-linear dynamics of the system is continuous, it is common to use a hybrid continuous-discrete process model. The non-linear continuous dynamics are used for the prediction step, usually through the application of a standard ordinary differential equations solver such as a Runge-Kutta method. The correction step is then completed using a discretized form of the linearized dynamics. This leads to a variant of the EKF known as continuous-discrete extended Kalman filter (CDEKF).

The derivation of the CDEKF follows from arguments of linearization and the Kalman filter algorithms. Specifically the filter is derived by linearizing the state and observation equations using a Taylor's series expansion.

The state and measurement equations are assumed to have the following form. State and measurement equations:{dot over (x)}(t)=f[x(t),u(t)]+Gw(t),  (16)z(t_i)=h[x(t_i)]+v(t_i).  (17)

Non-linear system equations Non-linear measurement equations

Propagation Equations:

$\begin{array}{cc}\hat{x}\left(t_i^{-}\right)=\hat{x}\left(t_{i-1}^{+}\right)+{\int }_{t_{i-1}}^{t_i}f\left[\hat{x}\left(t_{i-1}^{+}\right),u(t_{i-1})\right]dt,& \left(18\right)\\ P\left(t_i^{-}\right)=\Phi \left(t_i,t_{i-1}\right)P\left(t_{i-1}^{+}\right){\Phi \left(t_i,t_{i-1}\right)}^T+G\left(t_i\right){\mathrm{QG}\left(t_i\right)}^T,& \left(19\right)\end{array}$

Update Equations

$\begin{array}{cc}A\left(t_i\right)=H\left(t_i\right)P\left(t_i^{-}\right)H^T\left(t_i\right)+R\left(t_i\right),& \left(20\right)\\ K\left(t_i\right)=P\left(t_i^{-}\right)H^T\left(t_i\right){A\left(t_i\right)}^{-1},& \left(21\right)\\ r\left(t_i\right)=z\left(t_i\right)-h\left[\hat{x}\left(t_i^{-}\right)\right],& \left(22\right)\\ \hat{x}\left(t_i^{+}\right)=\hat{x}\left(t_i^{-}\right)+K\left(t_i\right)r\left(t_i\right),& \left(23\right)\\ P\left(t_i^{+}\right)=P\left(t_i^{-}\right)-K\left(t_i\right)H\left(t_i\right)P\left(t_i^{-}\right),\mathrm{where}& \left(24\right)\\ \Phi \left(t_i{t}_{i-1}\right)=\exp \left(F\left(t_{i-1}\right)\Delta t\right),& \left(25\right)\\ F\left(t_{i-1}\right)=\frac{\partial f\left(x,u\right)}{\partial x}{❘}_{\underset{u\left(t\right)=u\left(t_{i-1}\right)}{x\left(t\right)=\hat{x}\left(t_{i-1}^{+}\right)}},& \left(26\right)\\ H\left(t_i\right)=\frac{\partial h\left(x\right)}{\partial x}{❘}_{x\left(t\right)=\hat{x}\left(t_i^{-}\right)},& \left(27\right)\\ \Delta t=t_i-t_{i-1}.& \left(28\right)\end{array}$

It is clear that the EKF is very similar to the Kalman filter algorithm with the substitution of transition and measurement matrices by their linearized equivalents.

• • The EKF does not estimate state, but rather perturbations of state from an initial configuration. This is easily accommodated by adding the estimated state perturbation to the initial state.
  • The EKF requires a linearized model that must be computed from an approximate knowledge of state. This imposes a requirement for the filter to be initialized accurately to ensure that linearizations a valid descriptions of the dynamics.
  • Unlike the Kalman filter the, Kalman gain and the state covariance do not converge to a steady state. Moreover it is necessary to continually re-evaluate the Jacobian matrices.
  • If the predicted state trajectory drifts too far from the actual state trajectory, the true covariance will become much larger than the estimated covariance.
  • If the state or measurement equations are highly non-linear and the posterior density is non-Gaussian, the EKF may give a high estimation error.

These points notwithstanding, it is possible to achieve acceptable performance from an EKF in many applications.

The Multiple Model Adaptive Estimator

The Kalman filter and extended Kalman filter provide a method for establishing the dynamic state of a system by balancing the contribution made by a deterministic dynamic process model with that given by the measurement model. The core reason for undertaking state estimation when the objective is to estimate payload mass (and possibly other system parameters) is to use state estimates to implicitly compensate for the inertial loads resulting from machine motion.

The Kalman filter framework extends naturally to the estimation of system states and parameter simultaneously though a multiple model adaptive estimator (MMAE). The approach assumes the parameters to be estimated are constant or piecewise constant.

Generally, a system parameter can vary over a continuous range usually between maximum and a minimum values. The payload mass in the dipper of an electric mining shovel can take any value between zero (dipper empty) and the and a maximum capacity (associated with the dipper being full). The multiple model approach assumes that a parameter of interest, here the payload mass, can take on one of a number, say r, discrete values between these maximum and minimum limits. Denote these byMε{M_j}_j=1^r.  (29)

The idea is to hypothesize r system models, one for each of the r values, to estimate state of each model using a Kalman filter, and to determine the probability that each model is correct from the statistical properties of the innovation sequence generated by the filter for that model. Over time the most accurate model is assigned the largest probability while the less accurate models are assigned lower probabilities. Although payload mass has been used as an example in this discussion, we can in principle use this method on any combination of parameters such as payload mass and payload centre of mass to define a model.

The preferred embodiment uses the MAP approach to obtain an optimal estimate for the system parameters. However, whereas for the state estimation we were trying to find the most likely estimate of a state x(t_i) given the measurement sequence Z(t_i), we are now looking to determine the most likely model M_j, given Z(t_i).

The probability that a model is correct, given the measurement sequence Z(t_i), can expressed recursively as,

$\begin{array}{cc}\begin{array}{c}P\left\{M_j❘Z\left(t_i\right)\right\}=P\left\{M_j❘z\left(t_i\right),Z\left(t_{i-1}\right)\right\},\\ =\frac{f\left(z\left(t_i\right)❘M_j,Z\left(t_{i-1}\right)\right)P\left\{M_j❘Z\left(t_{i-1}\right)\right\}}{f\left(z\left(t_i\right)|Z\left(t_{i-1}\right)\right)},\\ =\frac{f\left(z\left(t_i\right)❘M_j,Z\left(t_{i-1}\right)\right)P\left\{M_j❘Z\left(t_{i-1}\right)\right\}}{\sum _{i=1}^{r}f\left(z\left(t_i\right)❘Z\left(t_{i-1}\right),M_i\right)P\left\{M_i❘Z\left(t_{i-1}\right)\right\}}.\end{array}& \left(30\right)\end{array}$

In the above expression we have made use of Bayes' rule and marginalization. If we use the following notation,μ_j(t_i)=P{M_j|Z(t_i)},μ_j(t_i−1)=P{M_j|Z(t_i−1)},then Eqn. 30 can be written as,

$\begin{array}{cc}{\mu }_j\left(t_i\right)=\frac{f\left(z\left(t_i\right)|M_j,Z\left(t_{i-1}\right)\right){\mu }_j\left(t_{i-1}\right)}{\sum _{k=1}^{r}f\left(z\left(t_i\right)|M_k,Z\left(t_{i-1}\right)\right){\mu }_k\left(t_{i-1}\right)}.& \left(31\right)\end{array}$

Note that if we have r models, we have r of these equations. The probability assigned to each model changes through times as more information becomes available. The initial condition isP(M_j|Z⁰)=μ_j(0),where μ_j(0) is the a priori probability that model j is the correct one. Note that the model probabilities obey the normalizing condition

$\sum _{j=1}^r{\mu }_j\left(t_i\right)=1.$

Equation 31 is a hybrid PDF, with both discrete probabilities and continuous PDFs. Looking at the individual terms in this expression we see that μ_j(t_i−1) is the probability that model j is correct given the measurement sequence up to the previous time step. The probability distribution f(z(t_i)|M_j,Z(t_i−1)) can be expressed for model M_j as,

$\begin{array}{cc}f\left(z\left(t_i\right)|M_j,Z\left(t_{i-1}\right)\right)={2\pi A_j\left(t_i\right)}^{-1/2}\exp \left\{-\frac{1}{2}{r_j\left(t_i\right)}^T{A}_j^{-1}{r}_j\left(t_i\right)\right\}& \left(32\right)\end{array}$

In the above equation, the residual r_j(t_i) and its covariance A_j(t_i) have been generated by assuming that M=M_j.

An estimate {circumflex over (M)} of the parameter M can then be found as a sum of the individual parameters weighted by their probabilities, i.e.

$\begin{array}{cc}\hat{M}\left(t_i\right)=\sum _{j=1}^r{M}_j{\mu }_j\left(t_i\right).& \left(33\right)\end{array}$

Similar expressions can derived for the combined state estimate,

$\begin{array}{cc}\hat{x}\left(t_i^{+}\right)=\sum _{j=1}^r{\mu }_j\left(t_i\right){\hat{x}}_j\left(t_i^{+}\right),& \left(34\right)\end{array}$ and state estimate covariance,

$\begin{array}{cc}\hat{P}\left(t_i^{+}\right)=\sum _{j=1}^r{\mu }_j\left(t_i\right)\left\{P_j\left(t_i^{+}\right)+{\left[{\hat{x}}_j\left(t_i\right)-\hat{x}\left(t_i\right)\right]\left[{\hat{x}}_j\left(t_i\right)-\hat{x}\left(t_i\right)\right]}^T\right\}.& \left(35\right)\end{array}$

This is a sum of probability weighted terms from each model and consists of two components. The variance on the estimate of the j-th model and an outer product due to the fact that {circumflex over (x)}(t_i) is not the same as {circumflex over (x)}_j(t_i).

Essentially, the multiple model approach of the preferred embodiment is as illustrated in 70 in FIG. 7 and involves running a bank of r independent Kalman filters e.g. 71, with each Kalman filter having a distinct model. The residual, residual covariance, and the prior model probability are then used to determine the current model probability 72. The updated probabilities are then used in a weighted sum 73 to determine the current parameter estimate. The MMAE algorithm is detailed in FIG. 8 in more detail.

Algorithm explanatory notes

• • 1. The initialization step sets initial conditions for the state estimate, state covariance. These are typically equal for all of the models. The a priori probability of each parameter is also set. If no prior information is available on parameter probabilities, the probability distribution should be flat and all initial probabilities should be equal to 1/r.
  • 2. The main for-loop begins at the start of the swing and finishes at the end of the swing. In practice start and finish times of the parameter estimator are determined by a cycle decomposition algorithm.
  • 3. In the second for-loop the Kalman filter is run for each for model. It independently updates and state estimate and state covariance for each model.
  • 4. The integration for the state prediction is achieved using a 4th order Runge-Kutta routine or similar ordinary differential equations solver.

When the second for-loop is completed the model probabilities are updated and a new parameter estimate is determined.

Machine Geometry and Coordinate Frames

FIG. 9 shows the parameters (lengths and angles) used to describe the geometry of P&H-class electric mining shovels and the coordinate frames used to describe relative positions of major moving assemblies of these machines. Lengths labeled l and angles labeled φ are fixed by design; length labeled d and angles labeled θ vary under machine motion.

Four frames are needed to describe the relative position of the major shovel assemblies. The world frame is denoted O₀x₀y₀z₀; the O₁x₁y₁z₁-frame is embedded in the machine house; the O₂x₂y₂z₂-frame in the saddle; and the O₃x₃y₃z₃-frame in the dipper. The x- and z-axes of all the body-fixed frames are in the sagittal plane of the machine house, that is the plane parallel to the plane of projection shown in FIG. 9 containing the swing axis. The y-axes of all frames are normal to this plane.

Four-by-four homogeneous transformation matrices can be used to describe the relationship between frames. We denote the matrix describing the transform from O₁x₁y_iz_i to O_jx_jy_jz_j by D_i→j and note that the action of this matrix maps (homogeneous) points in the j-frame to the i-frame. For instance, if p is a point whose location is known in O₃x₃y₃z₃, fixed in the bucket, then the coordinates of that point in frame O₀x₀y₀z₀ can be found from.p′=D_0→3p.

The structure of D_i→j is

$D_{i\to j}=\left(\begin{array}{cc}{R}_{i\to j}& t_{i\to j}\\ 0& 1\end{array}\right)$ where R_i→j is a 3×3 rotation matrix and t_i→j is a 3-dimensional translation vector. Four-by-four homogeneous transformation matrices commute according toD_i→k=D_i→jD_j→k.

The world frame O₀x₀y₀z₀ provides an inertial reference frame when the shovel tracks are locked. The origin of this frame is located at the interface between the upper surface of the tracks and the underside surface of the machine house. The z₀-axis collinear with the swing axis. The x₀-axis points in the direction of forward travel of the crawler tracks, and the y₀-axis completes a right-handed trihedral coordinate frame.

The origin O₁ of frame O₁x₁y₁z₁ is coincident with O₀ and z₁ is collinear with z₀. When θ₁=0, frames O₀x₀y₀z₀ and O₁x₁y₁z₁ coincide. A positive angle θ₁ corresponds to an anti-clockwise rotation of the machine house relative to the tracks when viewed from above. The homogeneous transformation matrix D_0→1 is given by

$\begin{array}{cc}{D}_{0\to 1}=\left(\begin{array}{cccc}\cos {\theta }_1& -\sin {\theta }_1& 0& 0\\ \sin {\theta }_1& \cos {\theta }_1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).& \left(36\right)\end{array}$

Frame O₂x₂y₂z₂ is fixed to the saddle with O₂ at the center of rotation of the saddle. When O₂ is equal to 0, the coordinate directions of O₂x₂y₂z₂ are parallel to those of O₁x₁y₁z₁. The displacement matrix describing the rigid body displacement from Frame 1 to Frame 2 is given by

$\begin{array}{cc}{D}_{1\to 2}=\left(\begin{array}{cccc}\cos {\theta }_2& 0& -\sin {\theta }_2& l_1\cos {\varphi }_1\\ 0& 1& 0& 0\\ \sin {\theta }_2& 0& \cos {\theta }_2& l_1\sin {\varphi }_1\\ 0& 0& 0& 1\end{array}\right),& \left(37\right)\end{array}$ where design parameters l₁ and φ₁ are as shown in FIG. 9

The origin O₃ of O₃x₃y₃z₃ is located as follows. The saddle angle (θ₂) is set equal to 90 degrees so that the handle is horizontal. The handle is then displaced such that the hoist rope falls vertically (θ₅=90 degrees and θ₆=0 degrees). The origin O₃ is located at the intersection of the pitch-line of the handle-rack and the hoist rope; z₃ is set to be collinear with the axis of the hoist rope; x₃ is set parallel to the pitch-line of the handle-rack. Note that axis x₂ is orthogonal to axis x₃. The displacement matrix describing the rigid body displacement D_2→3 is given by

$\begin{array}{cc}{D}_{2\to 3}=\left(\begin{array}{cccc}0& 0& 1& l_2\\ 0& 1& 0& 0\\ -1& 0& 0& -d_3\\ 0& 0& 0& 1\end{array}\right).& \left(38\right)\end{array}$

The multiplication of the above Eqns. gives the location of O₃x₃y₃z₃ in O₀x₀y₀z₀, namely

$\begin{array}{c}{D}_{0\to 3}=D_{0\to 1}{D}_{1\to 2}{D}_{2\to 3}\\ =\left(\begin{array}{cccc}\cos {\theta }_1\sin {\theta }_2& -\sin {\theta }_1& \cos {\theta }_1\cos {\theta }_2& \cos {\theta }_1\left(l_1\cos {\varphi }_1+l_2\cos {\theta }_2+d_3\sin {\theta }_2\right)\\ \sin {\theta }_1\sin {\theta }_2& \cos {\theta }_2& \sin {\theta }_1\cos {\theta }_2& \sin {\theta }_1\left(l_1\cos {\varphi }_1+l_2\cos {\theta }_2+d_3\sin {\theta }_2\right)\\ -\cos {\theta }_2& 0& \sin {\theta }_2& l_1\sin {\varphi }_1+l_2\sin {\theta }_2-d_3\cos {\theta }_2\\ 0& 0& 0& 1\end{array}\right).\end{array}$ Bail Pin Location

The dipper undergoes three parameter motion and it is useful to have an identified point on the dipper whose spatial location serves identify the position of O₃x₃y₃z₃. For this we use the bail-pin, located at b₃=(0,0,l₄,1) in O₃x₃y₃z₃. The coordinates of the bail-pin in O₀x₀y₀z₀ are

$\begin{array}{cc}{b}_0=D_{0\to 3}{b}_3=\left(\begin{array}{c}\cos {\theta }_1\left(l_1\cos {\varphi }_1+\left(l_2+l_4\right)\cos {\theta }_2+d_3\sin {\theta }_2\right)\\ \sin {\theta }_1\left(l_1\cos {\varphi }_1+\left(l_2+l_4\right)\cos {\theta }_2+d_3\sin {\theta }_2\right)\\ l_1\sin {\varphi }_1+\left(l_2+l_4\right)\sin {\theta }_2-d_3\cos {\theta }_2\\ 1\end{array}\right)& \left(39\right)\end{array}$ Constraint Equations

The swing angle, θ₁, the pivot angle, θ₂, and the crowd extension, d₃, parameterize the displacement and rotation of the body fixed frames relative to the world frame. We group these configuration variables as follows:θ=(θ₁,θ₂,d₃)^T.  (40)

The displacements of the swing motor θ_s, crowd motor θ_c, and hoist motor θ_h can be similarly grouped:ψ=(θ_s,θ_c,θ_h)^T.  (41)

The values of θ determine ψ and vice versa. These mappings are not bijective, however, within the physical working range of these variables their correspondence is one-to-one. Note that the specification of either θ or ψ determines the inclination of the hoist rope, a seventh variable, labeled θ₅ in FIG. 9.

To build up constraint equations, start by noting the dependent coordinates θ_s and θ_c are related θ₁ and d₃ by transmission ratios leading to

$\begin{array}{cc}0={\gamma }_1\left(\theta ,\psi \right)=\left[\begin{array}{c}{G}_s{\theta }_1-{\theta }_s\left(42\right)\\ G_c{d}_3-{\theta }_c\end{array}\right],& \left(43\right)\end{array}$ where G_s is the transmission ratio of the swing drive and G_c is the transmission ratio of the crowd drive.

The constraint equations relating θ₂ and θ₅ to θ_c and θ_h can be developed using the vector loop shown in FIG. 10. To simplify the notation, we work in a complex plane mapped to the physical (x₁, z₁)-plane with the real axis collinear with x₁ and imaginary axis with z₁.

Summing vectors in the vector loop of FIG. 10 where the vector z_i is expressed as a complex variable gives0=ξ(θ,ψ,θ₅)=z₁+z₂+z₃+z₄+z₅−z₆z₇,which can be expanded to give,0ξ(θ,ψ,θ₅)=l₁e^iφ1+l₂e^iθ2+d₂e^i(θ2−π2)+l₄e^iθ4+d₅e^iθ6−l₆e^iθ7.  (44)where the variables l_i, d_i, θ_i and φ_i are as defined in FIG. 10. Then using the following relations determined by inspection,

$\begin{array}{cc}{\theta }_6={\theta }_5-\frac{\pi }{2},& \left(45\right)\\ {\theta }_4={\theta }_2,& \left(46\right)\end{array}$

Eqn. 44 can be written as,ξ(θ,ψ,θ₅)=l₁e^iφ1+l₂e^iθ2+d₃e^i(θ2−π2)+l₄e^iθ2+d⁵e^iθ5−l₆e^{i(θ5−π/2)}−l₇e^iθ7.  (47)

To remove d₅ from this equation it is convenient to first introduce the variable d_h. As shown in FIG. 10, d_h represents the distance from the bail pin to the outer quadrant of the hoist sheave when the hoist rope hangs vertically (ie. θ₅=90 deg.). It is related to the angular displacement of the hoist motor by,

$\begin{array}{cc}{d}_h=\frac{{\theta }_h}{G_h},& \left(48\right)\end{array}$ where G_h, is the hoist transmission ratio. The above expression and Equation 45 can be used to relate d₅ and d_h,

$\begin{array}{cc}{d}_5=d_h+l_6{\theta }_6=\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right).& \left(49\right)\end{array}$

The last term on the right hand side of Eqn. 49 accounts for the angle of wrap of the hoist rope around the sheave. Substituting Eqn. 49 into Eqn. 47 gives

$\begin{array}{cc}\zeta \left(\theta ,\psi ,{\theta }_5\right)=l_1{e}^{i{\varphi }_1}+l_2{e}^{i{\theta }_2}+d_3{e}^{i\left({\theta }_2-\mathrm{\pi 2}\right)}+l_4{e}^{i{\theta }_2}+\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]e^{i{\theta }_5}-l_6{e}^{i\left({\theta }_5-\frac{\pi }{2}\right)}-l_7{e}^{i{\varphi }_7}=0.& \left(50\right)\end{array}$

Taking the real and imaginary components of Eqn. 50 relates θ_h and θ₅ to the generalized coordinates,

$\begin{array}{cc}{\gamma }_2\left(\theta ,\psi ,{\theta }_5\right)=\left[\begin{array}{cc}{l}_1\cos {\varphi }_1+l_2\cos {\theta }_2+d_3\sin {\theta }_2+l_4\cos {\theta }_2+\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\cos {\theta }_5-l_6\sin {\theta }_5-l_7\cos {\varphi }_7& \left(51\right)\\ l_1\sin {\varphi }_1+l_2\sin {\theta }_2-d_3\cos {\theta }_2+l_4\sin {\theta }_2+\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\sin {\theta }_5+l_6\cos {\theta }_5-l_7\sin {\varphi }_7& \left(52\right)\end{array}\right].& \left(53\right)\end{array}$

Note that use has been made of the following trigonometric relationships in arriving at Eqn. 51,

$\begin{array}{cc}\cos \left({\theta }_i-\frac{\pi }{2}\right)=\sin {\theta }_i,& \sin \left({\theta }_i-\frac{\pi }{2}\right)=-\cos {\theta }_i.\end{array}$

Concatenating Eqns. 42 and 51 gives

$\begin{array}{cc}\begin{array}{c}0=\Gamma \left(\theta ,\psi ,{\theta }_5\right)\\ =\left[\begin{array}{cc}{G}_s{\theta }_1-{\theta }_s& \left(54\right)\\ G_c{d}_3-{\theta }_c& \left(55\right)\\ l_1\cos {\varphi }_1+l_2\cos {\theta }_2+d_3\sin {\theta }_2+l_4\cos {\theta }_2+\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\cos {\theta }_5-l_6\sin {\theta }_5-l_7\cos {\varphi }_7& \left(56\right)\\ l_1\sin {\varphi }_1+l_2\sin {\theta }_2-d_3\cos {\theta }_2+l_4\sin {\theta }_2+\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\sin {\theta }_5+l_6\cos {\theta }_5-l_7\sin {\varphi }_7& \end{array}\right]\end{array}& \left(57\right)\end{array}$

All valid configurations must satisfy this Eqns. 54. In particular it forms the basis for the solving the kinematics, the statics, and the dynamics of P&H-class mining shovels.

Kinematic Tracking by Newton-Raphson

The kinematic tracking problem is to determine the values of θ and θ₅ given ψ or to determine ψ and θ₅ given θ. We call the first problem forward kinematic tracking and the second inverse kinematic tracking. For notational convenience in distinguishing between these two problems, we writeθ_F=(θ₁,θ₂,θ₃,θ₅)ψ_F=(θ_s,θ_c,θ_h)when working in the domain of the forward kinematic tracking problem andθ_I=(θ₁,θ₂,θ₃)ψ_I=(θ_s,θ_c,θ_h,θ₅)

The difference here is in with the grouping of θ₅. Both problems amount mathematically to solving the non-linear constraint equations. We chose to do this iteratively using a multi-variable Newton's method. The Jacobian matrices developed in expressing a Newton's method solution are used to refer motor inertias to the configuration variables and to solve the statics problem.

Forward Kinematic Tracking

Applying Taylor's series expansion to Eqn 54

$\Gamma \left({\theta }_F+{\mathrm{\Delta \theta }}_F,{\psi }_F+{\mathrm{\Delta \psi }}_F\right)=\Gamma \left({\theta }_F,{\psi }_F\right)+\frac{\mathrm{\delta \Gamma }\left({\theta }_F,{\psi }_F\right)}{{\mathrm{\delta \theta }}_F}{\mathrm{\Delta \theta }}_F+\frac{\mathrm{\delta \Gamma }\left({\theta }_F,{\psi }_F\right)}{{\mathrm{\delta \psi }}_F}{\mathrm{\Delta \psi }}_F+\mathrm{HOT}=0.$

The objective is to find a valid configuration, i.e. Γ(θ_F+Δθ_F,ψ_F′+Δψ_F)=0. It follows that

$\begin{array}{cc}\Gamma \left({\theta }_F,{\psi }_F\right)\approx -\left[\frac{\mathrm{\delta \Gamma }\left({\theta }_F,{\psi }_F\right)}{{\mathrm{\delta \theta }}_F}{\mathrm{\Delta \theta }}_F+\frac{\mathrm{\delta \Gamma }\left({\theta }_F,{\psi }_F\right)}{{\mathrm{\delta \psi }}_F}{\mathrm{\Delta \psi }}_F\right]& \left(58\right)\end{array}$

Equation 58 leads to the iteration equation

$\begin{array}{cc}{\mathrm{\Delta \theta }}_F^k=-{\left(\frac{\partial \Gamma \left({\theta }_F^{k-1},{\psi }_F^k\right)}{\partial {\theta }_F}\right)}^{-1}\left(\Gamma \left({\theta }_F^{k-1},{\psi }_F^k\right)+\frac{\partial \Gamma \left({\theta }_F^{k-1},{\psi }_F^k\right)}{\partial {\psi }_F}{\mathrm{\Delta \psi }}_F^k\right)& \left(59\right)\end{array}$ withθ_F^k=θ_F^k−1+Δθ_F^k  (60)

The algorithm below, gives an algorithm for kinematic tracking for P&H-class shovels. The algorithm takes the current motor positions, ψ_F^k and uses Eqns. 59 and 60 to find the new values of θ_F^k consistent with the constraint equations. For reliable convergence the algorithm requires a good initial values θ_F⁰. In practice this can be achieved by initializing from a well defined configuration such as provided in FIG. 10 where the forward kinematics can be explicitly solved using trigonometry. We note the bail-pin position can be easily found from 0 using Eqn. 39.

Algorithm 3: Forward Kinematic Tracking Using Newton's Method

input: Current motor position ψ_F^k.

output: Values of configuration variables θ_F^k consistent with the constraint equations.

priors: Previous motor and configuration variables: θ_F^k−1, ψ_F^k−1.

Initialization:Δψ_F^k=ψ_K^k−ψ_F^k−1 θ_F^k=θ_F^k−1 γ_F=Γ(θ_F^k,ω_F^k)Iterate Until Converged:

${\Gamma }_F<\mathrm{tol}{\mathrm{\Delta \theta }}_F^k=-\left(\Gamma \left({\theta }_F^k,{\psi }_F^k\right)+\frac{\partial \Gamma \left({\theta }_F^k,{\psi }_F^k\right)}{\partial {\psi }_F}{\mathrm{\Delta \psi }}_F^k\right)$ ${\theta }_F^k={\theta }_F^k+{\mathrm{\Delta \theta }}_F^k$ ${\Gamma }_F=\Gamma \left({\theta }_F^k,{\psi }_F^k\right)$

It is possible to write the Jacobian matrices,

$\begin{array}{cc}\frac{\partial \Gamma }{\partial {\theta }_f}=\left(\begin{array}{cccc}{G}_s& 0& 0& 0\\ 0& 0& G_c& 0\\ 0& -\left(l_2+l_4\right)\sin {\theta }_2+d_3\cos {\theta }_{\mathrm{2`}}& \sin {\theta }_2& -\sin {\theta }_5\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\\ 0& \left(l_2+l_4\right)\cos {\theta }_2+d_3\sin {\theta }_2& -\cos {\theta }_2& \cos {\theta }_5\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\end{array}\right)& \left(61\right)\\ \frac{\partial \Gamma }{\partial {\psi }_F}=\left(\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& \frac{\cos {\theta }_5}{G_h}\\ 0& 0& \frac{\sin {\theta }_5}{G_h}\end{array}\right)& \left(62\right)\end{array}$

Let

$J_F={\left(\frac{\partial \Gamma }{\partial {\theta }_F}\right)}^{-1}\frac{\partial \Gamma }{\partial {\psi }_F}.$ Evaluating expression symbolically we have

$\begin{array}{cc}{J}_F\left(\begin{array}{ccc}\frac{1}{G_s}& 0& 0\\ 0& \frac{\sin \left({\theta }_5-{\theta }_2\right)}{G_c{\gamma }_1}& -\frac{1}{G_h{\gamma }_1}\\ 0& \frac{1}{G_c}& 0\\ 0& \frac{G_h{d}_3}{G_c{\gamma }_1{\gamma }_2}& \frac{\left(l_2+l_4\right)\cos \left({\theta }_5-{\theta }_2\right)-d_3\sin \left({\theta }_5-{\theta }_2\right)}{{\gamma }_1{\gamma }_2}\end{array}\right)\mathrm{with}& \left(63\right)\\ {\gamma }_1=d_3\cos \left({\theta }_5-{\theta }_2\right)+\left(l_2+l_4\right)\sin \left({\theta }_5-{\theta }_2\right){\gamma }_2={\theta }_h+G_h{l}_6\left({\theta }_5-\frac{\pi }{2}\right)& \left(64\right)\end{array}$

Note this matrix relates motor speeds to the speeds of the configuration variables.

Inverse Kinematic Tracking

The same arguments that gave rise to Eqn. 59 and 60 can be used to form an iteration equation for inverse kinematic tracking, namely

$\begin{array}{cc}{\mathrm{\Delta \psi }}_I^k=-{\left(\frac{\partial \Gamma \left({\theta }_I^k,{\psi }_I^{k-1}\right)}{\partial {\psi }_I}\right)}^{-1}\left(\Gamma \left({\theta }_I^k,{\psi }_I^{k-1}\right)+\frac{\partial \Gamma \left({\theta }_I^k,{\psi }_I^{k-1}\right)}{\partial {\theta }_I}{\mathrm{\Delta \theta }}_I^k\right),\mathrm{where}& \left(65\right)\\ {\psi }_I^k={\psi }_I^{k-1}+{\mathrm{\Delta \psi }}_I^k& \left(66\right)\end{array}$

Symbolic forms for

$\frac{\partial \Gamma }{\partial {\psi }_I},\frac{\partial \Gamma }{\partial {\theta }_I},$ and the product

$J_I={\left(\frac{\partial \Gamma }{\partial {\psi }_I}\right)}^{-1}\frac{\partial \Gamma }{\partial {\theta }_I}$ are given as Eqns. 69.Algorithm 4: Inverse kinematic tracking using Newton's method

input: Current configuration variables θ_i^k.

output: Current motor position ψ_I^k.

priors: Previous motor and configuration variables: θ_i^k−1,ψ_i^k−1.

Initialization:Δθ_I^k=θ_I^k−θ_I^k−1 ψ_I^k=ψ_I^k−1 γ_I=γ(θ_i^k,ψ_I^k)Iterate Until Converged:

${\Gamma }_I<\mathrm{tol}{\mathrm{\Delta \psi }}_I^k=-{\left(\frac{\partial \Gamma \left({\theta }_I^k,{\psi }_I^k\right)}{\partial {\psi }_I}\right)}^{-1}\left(\Gamma \left({\theta }_I^k,{\psi }_I^k\right)+\frac{\partial \Gamma \left({\theta }_I^k,{\psi }_I^k\right)}{\partial {\theta }_I}{\mathrm{\Delta \theta }}_I^k\right)$ ${\psi }_I^k={\psi }_I^k+{\mathrm{\Delta \psi }}_I^k$ ${\Gamma }_I=\Gamma \left({\theta }_I^k,{\psi }_I^k\right)$

$\begin{array}{cc}\frac{\partial \Gamma }{\partial {\psi }_i}=\left(\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& \frac{\cos {\theta }_5}{G_h}& -\sin {\theta }_5\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\\ 0& 0& \frac{\sin {\theta }_5}{G_h}& \cos {\theta }_5\left[\frac{{\theta }_h}{G_h}+l_6\left({\theta }_5-\frac{\pi }{2}\right)\right]\end{array}\right),& \left(67\right)\\ \frac{\partial \Gamma }{\partial {\theta }_I}=\left(\begin{array}{ccc}{G}_s& 0& 0\\ 0& 0& G_c\\ 0& -\left(l_2+l_4\right)\sin {\theta }_2+d_3\cos {\theta }_2& \sin {\theta }_2\\ 0& \left(l_2+l_4\right)\cos {\theta }_2+d_3\sin {\theta }_2& -\cos {\theta }_2\end{array}\right),& \left(68\right)\\ J_I=\left(\begin{array}{ccc}{G}_s& 0& 0\\ 0& 0& G_c\\ 0& -G_h\left[\left(l_2+l_4\right)\sin \left({\theta }_5-{\theta }_2\right)+d_3\cos \left({\theta }_5-{\theta }_2\right)\right]& G_h\sin \left({\theta }_5-{\theta }_2\right)\\ 0& \frac{G_h\left[d_3\sin \left({\theta }_5-{\theta }_2\right)-\left(l_2+l_4\right)\cos \left({\theta }_5-{\theta }_2\right)\right]}{{\theta }_h+l_6{G}_h\left({\theta }_5-\pi /2\right)}& \frac{G_h\cos \left({\theta }_5-{\theta }_2\right)}{{\theta }_h+l_6{G}_h\left({\theta }_5-\pi /2\right)}\end{array}\right).& \left(69\right)\end{array}$ Rigid Body Dynamics of P&H-Class Electric Mining Shovels

The preferred embodiment therefore provides a framework for payload estimation on electric mining shovels, is implemented as a bank of Kalman filters running in parallel and estimating shovel state. Each model hypothesizes a value for the payload mass (and possibly other parameters) though a process model, and the probability of each hypothesis being true is evaluated using the statistical properties of the filter residuals. The “true” payload is computed by taking the expected value of the discrete distribution defined by the filter models.

There is a body of received wisdom, well established in the literature, that argues that a process model should be parsimonious—as simple as possible but not simpler—to minimize overall model error. Given the available measurement set, the preferred embodiment uses a process model that models only those dynamics states associated with the rigid body dynamics of the machine. Lagrangian methods are employed to derive the equations of motion describing these rigid body dynamics and develop a framework for computing the linearization of the rigid body dynamics about an arbitrary state trajectory needed for implementation of the extended Kalman filter.

Structure of Process Model

Since the purpose of the process model in the context of payload estimation is to factor out the effect of machine dynamics, it is desirable to measure the states of the model and the inputs directly and accurately. Intuitively, this results in minimal dispersion of information across unmeasured states. A Kalman filter running the correct payload hypothesis should produce state estimates that have greater accuracy than the sensors can measure them. And intuitively, this, in turn, should lead to better payload estimate.

The principle of measuring (rather than estimating) as many of the process model states as possible, argues for two process models: (i) a model comprising the electrical dynamics of the motors and rigid body dynamics of the machine house, saddle, and handle-dipper with the armature voltages serving as inputs, the armature currents and positions serving as measurements, and motor velocities estimated within the Kalman filter structure; or (ii) a model comprising only the rigid body dynamics of the machine house, saddle, and handle-dipper with motor currents serving as inputs and measurements being motor positions and armature voltages. The possible third option of a process model comprising the controller dynamics, the electrical dynamics, and rigid body dynamics is ruled out on the grounds that its more complicated than the two models identified above and requires access to the internal states of the PI and PID controllers, which in practice, are not easily measured.

The second model is chosen because it has fewer states (parsimony) and because the armature voltage and current can be used to obtain an estimate of motor velocity from Eqn. 70. This, of course, is also true for the for Model (i), however the motor dynamics would then be repeated in both the process model and the measurement equations to no apparent advantage.

Given this model there is the decision of which variables to use as states. The obvious candidates with the motor positions and velocities or the configuration variables θ and their derivatives {dot over (θ)}. These variables are kinematically linked via the constraint equations. The choice dictates to what generalized variables the inertias, motions, and forces are all referred. It transpires that the expressions are simpler if we choose the configuration variables and their velocities as states.

Rigid Body Dynamics

The kinematic model represented the shovel as three rigid bodies connected in series. These bodies were (i) machine house and boom, (ii) saddle, and (iii) the dipper and handle. Each body is driven by a motor having its own inertial properties through a transmission. The coupling between the bodies and the motors is assumed to be rigid. The three mass system will be expressed with respect to the configuration coordinates θ and their derivatives as defined previously.

A Lagrangian formulation can be used in conjunction with a screw theoretic description of the kinematic geometry of the shovel to give a high level, description of the dynamics. Lagrange's equations are founded on the principles of virtual work.

The task of deriving the equations describing the rigid body dynamics is complicated by the coupling that exists between the crowd and hoist drives that sees both act directly on the handle-dipper assembly. The constraint equations, derived in above take the form of an implicit relationship between configuration variables and motor position (i.e. γ(θ_i,ψ_Iθ_I))=0). Usually for constrained systems, where the dependent variables cannot be eliminated, a Lagrangian multiplier approach is used to develop the equations of motion. When that approach is applied to this problem, the resulting expressions are not in a form that can be readily applied within a Kalman filtering framework. To overcome this difficulty the equations of motion are developed here using embedded constraint equations. This treatment leads to an unaltered form of the Kalman filtering equations.

Lagrange's equations can be written

$\begin{array}{cc}\frac{d}{dt}\left(\frac{\partial L}{\partial {\stackrel{.}{\theta }}_i}\right)-\frac{\partial L}{\partial {\theta }_i}={\tau }_i,& \left(71\right)\end{array}$

The Lagrangian function (or Lagrangian for short) is given by L=T−V, with T the total kinetic energy, V the total potential energy of the system, and τ_i is the generalized force in the direction of the generalized coordinate θ_i. For this problem, the potential energy of the system is independent of velocity. Thus Eqn. 71 can be written as,

$\begin{array}{cc}\frac{d}{dt}\left(\frac{\partial T}{\partial {\stackrel{.}{\theta }}_i}\right)-\frac{\partial T}{\partial {\theta }_i}={\tau }_i+{\gamma }_i,\mathrm{where},& \left(72\right)\\ {\gamma }_i=\frac{\partial V}{\partial {\theta }_i}.& \left(73\right)\end{array}$

The total kinetic is equal to the sum of the kinetic energy of the rigid bodies comprising the front-end and the kinetic energy of the motor and transmission inertias.

$\begin{array}{cc}T=\frac{1}{2}\left({\theta }_I^T\stackrel{\_}{A}{\theta }_I+{\stackrel{.}{\psi }}_I^T\tilde{A}{\psi }_I\right).& \left(74\right)\end{array}$

In the above equation Ā is a 3×3 matrix describing the inertial properties of the machine house and boom, the saddle, and the handle and dipper assembly referred to the configuration coordinates θ. We refer to the kinetic energy associated with this inertia, the external body kinetic energy because it corresponds to the kinetic energy the external bodies of the machine.

The matrix à is a 4×4 matrix describing the inertias of the motor rotors and transmission elements referred to the dependent variables ψ. We refer to the component of kinetic energy associated with this inertia, the internal kinetic energy.

External Kinetic Energy

The external kinetic energy aggregates energy associate with the mechanical dynamics of the house, boom, saddle, handle, dipper, bail pin, hoist rope, and the swing, crowd and hoist drives. For the development of the equations the house and boom, the saddle, and handle dipper and treated as three rigid bodies. To reduce the complexity of the dynamics, consistent with the principle of parsimony, the following assumptions were made:

1. The hoist rope is inextensible,

2. The hoist rope is massless,

No assumptions are made in the formulation about the direction of the gravity vector in the O₀x₀y₀z₀ frame, although in the absence of inclinometer information, gravity is taken as aligning with the z₀ axis.

To form an expression for the external kinetic energy, a screw theoretic formulation is used. The formulation is based in the observation that the canonical description of an infinitesimal rigid body motion corresponds to a rotation about a defined spatial line together with a translation along that line. That is, the most general instantaneous motion of a rigid body is a spatial twisting. If the motion corresponds, instantaneously, to a translation the line about which the body twists (screw axis) is in the plane at infinity. For all other motions, the corresponding screw axes are lines in affine space. Twisting motions can be concisely described by a 6-dimensional vector of Plücker coordinate, describing the translational and rotational velocities

$\begin{array}{cc}\stackrel{.}{q}=\left(\begin{array}{c}\omega \\ v\end{array}\right),\mathrm{where}\omega ={\left({\omega }_x,{\omega }_y,{\omega }_z\right)}^T,v={\left(v_x,v_y,v_z\right)}^T,& \left(75\right)\end{array}$ are, respectively, the rotational velocity and the translational velocity of the body in the same frame of reference? The angular velocity, ω is a free vector and the translational velocity, v, can be interpreted as the velocity of the (possibly fictitious) point of the body instantaneously at the origin of the frame of reference.

In machines, the relative motion of bodies is commonly constrained by physical joints that allow one degree of relative freedom. The Plücker coordinates describing the relative motion can be written{dot over (q)}_i={dot over (θ)}_is_i where θ_i is the relative speed of the joint and s_i is the normalized Plücker coordinates for the motion known as the joint screw. The convention is to normalize {dot over (q)}_i by dividing each element by the magnitude of the angular speed ω for rotational and helical motions and by magnitude of v for translational motions. In this way the joint screws serve as a basis for describing the instantaneous motion of a linkage. Relative velocities add in the usual way provided all motions are expressed in the same frame of reference.

For example, the velocities of the three external bodies of an EMS are{dot over (q)}₁={dot over (θ)}₁s₁ {dot over (q)}₂={dot over (θ)}₁s₁+{dot over (θ)}₂s₂ {dot over (q)}₃=θ₁s₁+θ₂s₂+d₃s₃  (76)where {dot over (θ)}₁, {dot over (θ)}₂, and {dot over (d)}₃ are the velocities of the configuration variables ands, s₁, s₂, s₃ are the joint screws in for the three freedoms in the O₀x_oy₀z₀-frame whose Plücker coordinates ares₁=(0,0,1;0,0,0)^T s₂=(sin θ₁,−cos θ₁,0;l₁ sin φ₁ cos θ₁,l₁ sin φ₁ sin θ₁,−l₁ cos φ₁)^T  (77)s₃=(0,0,0;cos θ₁ sin θ₂,sin θ₁ sin θ₂,cos θ₂)^T  (77)

The kinetic energy of a rigid body instantaneously twisting with the Plücker vector expressed in an inertial frame is given by

$\begin{array}{cc}T=\frac{1}{2}\left(\begin{array}{cc}\omega & v\end{array}\right)\left(\begin{array}{cc}J& \mathrm{mC}\\ {\mathrm{mC}}^T& {\mathrm{mI}}_3\end{array}\right)\left(\begin{array}{c}\omega \\ v\end{array}\right)& \left(78\right)\end{array}$ where m is the mass of the body, J is the 3×3 rotational inertia matrix of the body in the inertial reference frame

$\begin{array}{cc}J=\left(\begin{array}{ccc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}& I_{\mathrm{xz}}\\ I_{\mathrm{xy}}& I_{\mathrm{yy}}& I_{\mathrm{yz}}\\ I_{\mathrm{xz}}& I_{\mathrm{yz}}& I_{\mathrm{zz}}\end{array}\right),& \left(79\right)\end{array}$

and C is the skew-symmetric form of the centre-of-mass vector of the body, c=(c_x, c_y, c_z)^T, in the inertial frame and

$\begin{array}{cc}C=\left(\begin{array}{ccc}0& -c_z& c_y\\ c_z& 0& -c_x\\ -c_y& c_x& 0\end{array}\right).& \left(80\right)\end{array}$

Note the skew-symmetric form of c takes the cross product with a 3-vector that post multiplies it.

Equation 78 can be written in shorthand notation

$\begin{array}{cc}T=\frac{1}{2}{\stackrel{.}{q}}_b^T{N}^b{\stackrel{.}{q}}_b& \left(81\right)\end{array}$ where the 6×6 matrix N^b known as the inertia matrix. Inertia matrices are most easily expressed in a body fixed frame b. The transformation of an inertia from a (possibly moving) body-fixed frame to another frame, e.g. an inertial frame can be achieved usingN⁰=H_0→b^TN^bH_0→b⁻¹.  (82)where H_0→b is a 6×6 screw transformation matrix having structure.

$\begin{array}{cc}{H}_{0->b}=\left(\begin{array}{cc}{R}_{0->b}& 0\left(19\right)\\ T_{0->b}{R}_{0->b}& R_{0->b}\left(20\right)\end{array}\right)& \left(85\right)\end{array}$

The matrix R_0→b is a 3×3 rotation matrix and T is the skew symmetric form of the translation vector t_0→b. From Eqns. 76 and 81 is hence given by

$\begin{array}{c}T=\frac{1}{2}\left({\stackrel{.}{q}}_1^T{N}_1^0{\stackrel{.}{q}}_1+{\stackrel{.}{q}}_2^T{N}_2^0{\stackrel{.}{q}}_2+{\stackrel{.}{q}}_3^T{N}_3^0{\stackrel{.}{q}}_3\right)\\ =\sum _{i=1}^3\sum _{j=i}^3{\stackrel{.}{\theta }}_i{s}_i^T\left(\sum _{k=j}^3{N}_k^0\right)s_j{\stackrel{.}{\theta }}_j\end{array}$ where for notational convenience {dot over (θ)}₃ has been substituted for {dot over (d)}₃. When written in terms of the derivatives of the configuration variables, θ by

$\begin{array}{cc}T=\frac{1}{2}{\stackrel{.}{\theta }}_I^T\stackrel{\_}{A}{\stackrel{.}{\theta }}_I\mathrm{where}{\stackrel{\_}{a}}_{\mathrm{ij}}=s_i^T\left(\sum _{k=l}^3{N}_k^0\right)s_j,& \left(86\right)\end{array}$ with l=max(i,j). This expression is in a form amenable to expansion via a symbolic manipulation package. The requirements (i) are values for the inertias in body fixed frames of reference; (ii) transformations of the form of Eqn. 81 which, for a P&H-class mining shovel; and (iii) the joint screws which are given by Eqns. 77.Internal Kinetic Energy

The internal kinetic energy is given by

$\begin{array}{cc}{T}_i=\frac{1}{2}{\stackrel{.}{\psi }}_I^T\tilde{A}{\stackrel{.}{\psi }}_I,\mathrm{where}\tilde{A}=\left(\begin{array}{cccc}{J}_s& 0& 0& 0\\ 0& J_c& 0& 0\\ 0& 0& J_h& 0\\ 0& 0& 0& 0\end{array}\right).& \left(87\right)\end{array}$

Here J_s is the inertia of the swing armature and swing transmission referred to the rotor axis; J_c is the inertia of the crowd armature and transmission referred to the rotor axis; and J_h the inertia of the hoist armature and transmission referred to the rotor axis. We desire to express this energy in the configuration variables θ.

The derivative of the constraint equation, with respect to time is

$\begin{array}{c}0=\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{dt},\\ =\frac{\partial \Gamma \left({\theta }_I,{\psi }_I\right)}{\partial {\psi }_I}\frac{d{\psi }_I}{dt}+\frac{\partial \Gamma \left({\theta }_I,{\psi }_I\right)}{\partial {\theta }_I}\frac{d{\theta }_I}{dt},\\ =\frac{\partial \Gamma \left({\theta }_I,{\psi }_I\right)}{\partial {\psi }_I}{\stackrel{.}{\psi }}_I+\frac{\partial \Gamma \left({\theta }_I,{\psi }_I\right)}{\partial {\theta }_I}{\theta }_I.\end{array}$

In the above equation

$\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\psi }_I}$ is the 4×4 matrix and

$\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\theta }_I}$ is a 4×3 matrix. Provided

$\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\psi }_I}$ has full rank, the above expression can be arranged as

$\begin{array}{cc}{\psi }_I=-{\left(\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\psi }_I}\right)}^{-1}\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\theta }_I}{\stackrel{.}{\theta }}_I,=J_i\left({\theta }_I,{\psi }_I\right){\theta }_I,\mathrm{where}{J}_i\left({\theta }_I,{\psi }_I\right)=-{\left(\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\psi }_I}\right)}^{-1}\frac{d\Gamma \left({\theta }_I,{\psi }_I\right)}{d{\theta }_I}& \left(88\right)\end{array}$

We call J_I(θ_I,ψ_I) the constraint Jacobian. It relates the velocities of the configuration variables to those of the dependent variables.

The internal kinetic energy expressed in the configuration variables is thereforeT_i={dot over (θ)}^TJ_I(θ,ψ)^TÃJ_I(θ_I,ω_I){dot over (θ)}_I.  (89)

Again the expression is in a form that readily facilitates symbolic evaluation. The total kinetic energy is

$\begin{array}{cc}T=\frac{1}{2}\left[{\stackrel{.}{\theta }}_I^T\stackrel{\_}{A}{\stackrel{.}{\theta }}_I+{\stackrel{.}{\theta }}_I^T{J_I\left({\theta }_I,{\psi }_I\right)}^T\tilde{A}{J}_I\left({\theta }_I,{\psi }_I\right){\stackrel{.}{\theta }}_I\right]=\frac{1}{2}{\stackrel{.}{\theta }}_I^{T}A{\stackrel{.}{\theta }}_I,\mathrm{where}A=\stackrel{\_}{A}+{J_I\left(\theta ,\psi \right)}^T\tilde{A}{J}_I\left(\theta ,\psi \right).& \left(90\right)\end{array}$

In summation notation we can write Eqn. 90,

$\begin{array}{cc}T=\frac{1}{2}\sum _{i=1}^3\sum _{j=1}^3{a}_{\mathrm{ij}}{\stackrel{.}{\theta }}_i{\stackrel{.}{\theta }}_j.& \left(91\right)\end{array}$

The inertia terms a_ij are functions of both the dependent and independent coordinates due to the introduction of J_I(Q_I,ω_I), but are not explicitly dependent on time. To simplify the notation, the dependence of J_I(θ_I,ω_I) on θ_I and ψ_I will no longer be shown.

Evaluation of Lagrange's Equation

It is necessary to evaluate the derivatives on the left hand side of the Lagrange's equation. As the kinetic and potential energy are independent we shall evaluate their contribution to the Lagrange's equations separately.

Using Eqn. 91, the first term on the left hand side of Eqn. 72 becomes,

$\begin{array}{cc}\begin{array}{c}\frac{d}{dt}\left(\frac{\partial T}{\partial {\stackrel{.}{\theta }}_i}\right)=\frac{d}{dt}\left[\frac{\partial }{\partial {\stackrel{.}{\theta }}_i}\left(\frac{1}{2}\sum _{i=1}^3\sum _{j=1}^3{a}_{\mathrm{ij}}{\stackrel{.}{\theta }}_i{\stackrel{.}{\theta }}_j\right)\right],\\ =\frac{d}{dt}\left(\sum _{j=1}^3{a}_{\mathrm{ij}}{\stackrel{.}{\theta }}_j\right),\\ =\sum _{j=1}^3\left(\frac{d{a}_{\mathrm{ij}}}{dt}{\stackrel{.}{\theta }}_j+a_{\mathrm{ij}}{\ddot{\theta }}_j\right).\end{array}& \left(92\right)\end{array}$

Using the chain rule, the derivative of a_ij with respect to time is

$\begin{array}{c}\frac{d{a}_{\mathrm{ij}}}{dt}=\frac{\partial a_{\mathrm{ij}}}{\partial \psi }\frac{d\psi }{dt}+\frac{\partial a_{\mathrm{ij}}}{\partial \theta }\frac{d\theta }{dt},\\ =\frac{\partial a_{\mathrm{ij}}}{\partial \psi }\stackrel{.}{\psi }+\frac{\partial a_{\mathrm{ij}}}{\partial \theta }\stackrel{.}{\theta }.\end{array}$

Using Eqn. 88 we can remove {dot over (ψ)} from the above expression,

$\begin{array}{cc}\frac{d{a}_{\mathrm{ij}}}{dt}=\left(\frac{\partial a_{\mathrm{ij}}}{\partial \psi }{J}_I+\frac{\partial a_{\mathrm{ij}}}{\partial \theta }\right)\stackrel{.}{\theta }.& \left(93\right)\end{array}$

Alternatively we can write,

$\begin{array}{cc}\frac{d{a}_{\mathrm{ij}}}{dt}=\sum _{k=1}^3\left(\sum _{l=1}^4\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{\partial a_{\mathrm{ij}}}{\partial {\theta }_k}\right){\stackrel{.}{\theta }}_k.& \left(94\right)\end{array}$

Here J_I,lk is the (l,k)-th element of J_I. This allows Eqn. 92 to be expressed as

$\begin{array}{cc}\begin{array}{c}\frac{d}{dt}\left(\frac{\partial T}{\partial {\stackrel{.}{\theta }}_i}\right)=\sum _{j=1}^3\left(\frac{d{a}_{\mathrm{ij}}}{dt}{\stackrel{.}{\theta }}_j+a_{\mathrm{ij}}{\ddot{\theta }}_j\right),\\ =\sum _{j=1}^3\left[\sum _{k=1}^3\left(\sum _{l=1}^4\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{\partial a_{\mathrm{ij}}}{\partial {\theta }_k}\right){\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k+a_{\mathrm{ij}}{\ddot{\theta }}_j\right].\end{array}& \left(95\right)\end{array}$

Considering now the second term on the left hand side of Lagrange's equations, we have

$\begin{array}{cc}\begin{array}{c}\frac{\partial T}{\partial {\theta }_i}=\frac{\partial }{\partial {\theta }_i}\left(\frac{1}{2}\sum _{j=1}^3\sum _{k=1}^3{a}_{\mathrm{jk}}{\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k\right),\\ =\frac{1}{2}\sum _{j=1}^3\sum _{k=1}^3\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}{\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k.\end{array}& \left(96\right)\end{array}$

Applying the chain rule to the partial differential term in this equation gives

$\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}=\frac{\partial a_{\mathrm{jk}}}{\partial \psi }\frac{\partial \psi }{\partial \theta }\frac{\partial \theta }{\partial {\theta }_i}+\frac{\partial a_{\mathrm{jk}}}{\partial \theta }\frac{\partial \theta }{\partial {\theta }_i}.$

The time derivative of the dependent variables, {dot over (ψ)}, can be written

$\begin{array}{cc}\stackrel{.}{\psi }=\frac{\partial \psi }{\partial \theta }\stackrel{.}{\theta }.& \left(97\right)\end{array}$

Equating the above expression with Equation 24 gives

$\begin{array}{cc}\frac{\partial \psi }{\partial \theta }=J_I.& \left(98\right)\end{array}$

Then using Eqn. 98 to eliminate

$\frac{\partial \psi }{\partial \theta }$ from Eqn. 97 gives

$\begin{array}{cc}\begin{array}{c}\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}=\frac{\partial a_{\mathrm{jk}}}{\partial \psi }{J}_I\frac{\partial \theta }{\partial {\theta }_i}+\frac{\partial a_{\mathrm{jk}}}{\partial \theta }\frac{\partial \theta }{\partial {\theta }_i},\\ =\left(\frac{\partial a_{\mathrm{jk}}}{\partial \psi }{J}_I+\frac{\partial a_{\mathrm{jk}}}{\partial \theta }\right)\frac{\partial \theta }{\partial {\theta }_i}.\end{array}& \left(99\right)\end{array}$

As the generalized coordinates are independent, it follows that

$\begin{array}{cc}\frac{\partial \theta }{\partial {\theta }_i}=\{\begin{array}{cc}1& i=j\\ 0& i\ne j.\end{array}& \left(100\right)\end{array}$

Applying this expression Eqn. 99 becomes,

$\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}=\sum _{l=1}^4\frac{\partial a_{\mathrm{jk}}}{\partial {\psi }_l}{J}_{I,\mathrm{li}}+\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i},$ allowing Eqn. 96 can be written

$\begin{array}{cc}\frac{\partial T}{\partial {\theta }_i}=\frac{1}{2}\sum _{j=1}^3\sum _{k=1}^3\left(\sum _{l=1}^4\frac{\partial a_{\mathrm{jk}}}{\partial {\psi }_l}{J}_{I,\mathrm{li}}+\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}\right){\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k.& \left(101\right)\end{array}$

Using Eqn.s 95 and 101 the contribution of the kinetic energy terms to Lagrange's equations becomes,

$\begin{array}{cc}\begin{array}{c}\frac{d}{dt}\left(\frac{\partial T}{\partial {\stackrel{.}{\theta }}_i}\right)-\frac{\partial T}{\partial {\theta }_i}=\sum _{j=1}^3{a}_{\mathrm{ij}}{\ddot{\theta }}_j+\sum _{j=1}^3\sum _{k=1}^3\\ \left[\sum _{l=1}^4\left(\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}-\frac{1}{2}\frac{\partial a_{\mathrm{jk}}}{\partial {\psi }_l}{J}_{I,\mathrm{li}}\right)+\frac{\partial a_{\mathrm{ij}}}{\partial {\theta }_k}-\frac{1}{2}\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}\right]{\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k,\\ =\sum _{j=1}^3{a}_{\mathrm{ij}}{\ddot{\theta }}_j+\sum _{j=1}^3\sum _{k=1}^3{b}_{\mathrm{ijk}}{\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k,\end{array}\mathrm{where}& \left(102\right)\\ b_{\mathrm{ijk}}=\sum _{l=1}^4\left(\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_I}{J}_{I,\mathrm{lk}}-\frac{1}{2}\frac{\partial a_{\mathrm{jk}}}{\partial {\psi }_l}{J}_{I,\mathrm{li}}\right)+\frac{\partial a_{\mathrm{ij}}}{\partial {\theta }_k}-\frac{1}{2}\frac{\partial a_{\mathrm{jk}}}{\partial {\theta }_i}.& \left(103\right)\end{array}$ Gravity Induced Generalized Forces

A screw theoretic approach can also be used to describe the combined action of forces and moments that act on a body. More specifically, any system of forces acting on a body can be reduced to a force acting along a unique line in space together with a moment acting about that line. This couplet of a force F and a moment C is known as a wrench and takes the form of a six dimensional vectorW=(C;F)^T

Just as screws provide a canonical description of infinitesimal rigid body motions, wrenches provide a canonical description of the net force acting on a body.

Gravity imposes a wrench loading on each body shovel model. Suppose N is the inertia matrix for a rigid body expressed in an inertial frame. The wrench acting on this body due to gravity is given byG=Ng  (104)where g is the six-dimensional vector describing gravitational acceleration. In FIG. 9, the z-axis of the inertial frame is aligned with the direction of gravitational acceleration. Hereg=(0,0,0;0,0,g)^T  (105)with g=−9.81 m/s².

Suppose the body is constrained to motion about a joint screw s and let the angle of rotation θ of the body about that screw axis be the generalized coordinate for the body. The component of the generalized force due to gravitational loading, τ_g can be found from the inner product of the gravitational wrench and the joint screw, i.e.τ_g=Ws

The gravitational wrenches acting on each of the rigid bodies of the shovel model areG₁⁰=(N₁⁰+N₂⁰+N₃⁰)g G₂⁰=(N₁⁰+N₂⁰)g G₃⁰=N₃⁰gor in general

$G_i^0=\sum _{j=1}^3{N}_i^{0}g$

The refer G_i to the generalized coordinates, we pair each gravitational wrench with the joint screw, i.e.γ_i=G_i^Ts_i.  (106)

If the shovel is inclined so that the unit vector n=(n_x, n_y, n_z) defines the direction of gravity in the O₀x₀y₀z₀-frame,g=(0,0,0;gn_x,gn_y,gn_z)^T  (107)

The remainder of the formulation remains the same. In a practical implementation, the direction cosines n could be measured by inclinometers. In this instance γ is the 3-dimensional vector

$\begin{array}{cc}\gamma =\left(\begin{array}{c}0\\ g\left[\left(m_p{c}_{\mathrm{xp}}+m_3{c}_{x3}-m_2{c}_{z2}\right)\sin {\theta }_2+\left(m_2{c}_{x2}+m_3{c}_{z3}+m_p{c}_{\mathrm{zp}}\right)\cos {\theta }_2+\left(m_3+m_p\right)d_3\sin {\theta }_2+\left(m_3+m_p\right)l_2\cos {\theta }_2\right]-\\ g\left(m_3+m_p\right)\cos {\theta }_2\end{array}\right)& \left(108\right)\end{array}$ Generalized Actuator Forces

Let τ be the vector of generalized actuator forces, i.e.τ=(τ₁,τ₂,τ₃)  (109)

These generalized forces are generated by net motor torquesT=(T_s,T_c,T_h,0)^T.  (110)where T_s is the net torque applied by the swing drive, T_c the net torque applied by the crowd drive, and T_h net torque applied by the hoist drives. These torques are the electromagnetic torques dependent on the drives less the losses through the transmissions. These torques act on the dependent coordinates ψ. Note the final element of T in Eqn. 110 corresponds to the torque applied on θ₅; this is set to zero as no torque acts in that direction.

It is required to relate T to τ which can be done using the constraint Jacobian and the principle of virtual work. The rate at which the net motor torques do work is given by the inner product of T with {dot over (ψ)}. Since the constraints, by definition, can do no work, it must follow that{dot over (ψ)}^TT={dot over (θ)}^Tτ

Substituting in Eqn. 88 gives the relationship between the net motor torques and the generalized torques,τ=J_I^TT  (111)

It therefore follows that:

$\begin{array}{cc}\tau =\left(\begin{array}{c}{G}_s{T}_s\\ \left(-G_h\left(l_2+l_4\right)\sin \left({\theta }_5-{\theta }_2\right)-G_h{d}_3\cos \left({\theta }_5-{\theta }_2\right)\right)T_h\\ G_c{T}_c+G_h\sin \left({\theta }_5-{\theta }_2\right)T_h\end{array}\right)& \left(112\right)\end{array}$

The sum of Eqns 108 and 112 equated to zero defined the static equations for a P&H-class mining shovel referred to the generalized coordinates θ.

Damping Losses

We can incorporate viscous damping losses into the equations of motion as follows, in terms of the time derivatives of the dependent variables.T_D={tilde over (D)}{dot over (ψ)}where {tilde over (D)} represents the damping due to the motors and the damping on the joints. Thenβ=−^J_I^TT_D=−J_I^T{tilde over (D)}J_{I{dot over (θ)}=D{dot over (θ)}}  (113)whereD=J_I^T{tilde over (D)}J_I Summary of the Rigid Body Dynamics

In summary, the equations governing rigid body motion of the machine can be written

$\begin{array}{cc}\sum _{j=1}^n{a}_{\mathrm{ij}}{\ddot{\theta }}_j+\sum _{j=1}^n\sum _{k=1}^n{b}_{\mathrm{ijk}}{\stackrel{.}{\theta }}_j{\stackrel{.}{\theta }}_k={\tau }_i+{\gamma }_i+{\beta }_i& \left(114\right)\end{array}$ or in matrix form

$\begin{array}{cc}A\ddot{\theta }+B{\stackrel{.}{\theta }}_{\mathrm{ii}}+C{\stackrel{.}{\theta }}_{\mathrm{ij}}=\tau +\gamma +\beta ,\mathrm{where}& \left(115\right)\\ A=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right),& \left(116\right)\\ B=\left(\begin{array}{ccc}{b}_{111}& b_{122}& b_{133}\\ b_{211}& b_{222}& b_{233}\\ b_{311}& b_{322}& b_{333}\end{array}\right),& \left(117\right)\\ C=\left(\begin{array}{ccc}{b}_{112}+b_{121}& b_{113}+b_{131}& b_{123}+b_{132}\\ b_{212}+b_{221}& b_{213}+b_{231}& b_{223}+b_{232}\\ b_{312}+b_{321}& b_{313}+b_{331}& b_{323}+b_{332}\end{array}\right),\mathrm{and}\tau =\left(\begin{array}{c}{\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right),\gamma =\left(\begin{array}{c}{\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\end{array}\right),\beta =\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\end{array}\right),\ddot{\theta }=\left(\begin{array}{c}{\ddot{\theta }}_1\\ {\ddot{\theta }}_2\\ {\ddot{\theta }}_3\end{array}\right),{\stackrel{.}{\theta }}_{\mathrm{ii}}=\left(\begin{array}{c}{\stackrel{.}{\theta }}_1{\stackrel{.}{\theta }}_1\\ {\stackrel{.}{\theta }}_2{\stackrel{.}{\theta }}_2\\ {\stackrel{.}{\theta }}_3{\stackrel{.}{\theta }}_3\end{array}\right),{\stackrel{.}{\theta }}_{\mathrm{ij}}=\left(\begin{array}{c}{\stackrel{.}{\theta }}_1{\stackrel{.}{\theta }}_2\\ {\stackrel{.}{\theta }}_1{\stackrel{.}{\theta }}_3\\ {\stackrel{.}{\theta }}_2{\stackrel{.}{\theta }}_3\end{array}\right).& \left(118\right)\end{array}$

Note the inputs to this model are the net generalized torques c and the states are the generalized coordinates 0 and their derivatives 6.

Linearized Dynamics

The EKF algorithm requires the linearization of the dynamics about state trajectory in discrete-time form. The method for developing this linearization is first to determine the continuous time linearization about a state trajectory{dot over (x)}=Fx+Bu and then convert to the discrete time form required by the CDEKF algorithm, namelyx(t_i)=Φ(t_it_i−1)x(t_i−1)+γ(t_i−1)u(t_i−1)  (119)z(t_i)=H(t_i)x(t_i).  (120)

This section develops this linearization. The starting point is to rearrange Eqn. 115 so we write the rigid body dynamics as

$\begin{array}{cc}\left(\begin{array}{c}\stackrel{.}{\theta }\\ \ddot{\theta }\end{array}\right)=\left(\begin{array}{c}\stackrel{.}{\theta }\\ A^{-1}\left(\tau +\gamma -B{\stackrel{.}{\theta }}_{\mathrm{ii}}-C{\stackrel{.}{\theta }}_{\mathrm{ij}}\right)\end{array}\right).& \left(121\right)\end{array}$

Equation 121 has the same form as previous described, namely{dot over (x)}(t)=f[x(t),u(t)],  (122)andx^T=(θ,{dot over (θ)}),  (123)u^T=(T,g).  (124)

Note that the inputs are T, the motor torques, and g, the gravitational acceleration constant.

Linearized System Matrix

Taking the partial differential of Eqn. 121 with respect to the state vector gives the following block-form for the linearized system matrix F

$\begin{array}{cc}F=\left(\begin{array}{cc}\frac{\partial \stackrel{.}{\theta }}{\partial \theta }& \frac{\partial \stackrel{.}{\theta }}{\partial \stackrel{.}{\theta }}\\ \frac{\partial \ddot{\theta }}{\partial \theta }& \frac{\partial \ddot{\theta }}{\partial \stackrel{.}{\theta }}\end{array}\right)=\left(\begin{array}{cc}0& I_3\\ \frac{\partial \ddot{\theta }}{\partial \theta }& \frac{\partial \ddot{\theta }}{\partial \stackrel{.}{\theta }}\end{array}\right)& \left(125\right)\end{array}$

An expression for the sub-block

$\frac{\partial \ddot{\theta }}{\partial \theta }$ can be found by differentiating Eqn. 121 with respect to the configuration variable θ_k, k=1, 2, 3 gives

$\begin{array}{cc}\begin{array}{c}\frac{\partial \ddot{\theta }}{\partial {\theta }_k}=\frac{\partial }{\partial {\theta }_k}\left[A^{-1}\left(-B{\stackrel{.}{\theta }}_{\mathrm{ii}}-C{\stackrel{.}{\theta }}_{\mathrm{ij}}+\tau +\gamma \right)\right],\\ =-A^{-1}\frac{\partial A}{\partial {\theta }_k}{A}^{-1}\left(-B{\stackrel{.}{\theta }}_{\mathrm{ii}}-C{\stackrel{.}{\theta }}_{\mathrm{ij}}+\tau +\gamma \right)+\\ A^{-1}\left(-\frac{\partial B}{\partial {\theta }_k}{\stackrel{.}{\theta }}_{\mathrm{ii}}-\frac{\partial C}{\partial {\theta }_k}{\stackrel{.}{\theta }}_{\mathrm{ij}}+\frac{\partial \tau }{\partial {\theta }_k}+\frac{\partial \gamma }{\partial {\theta }_k}\right),\end{array}& \left(126\right)\end{array}$ where use has been made of the following relationship,

$\frac{\partial }{\partial {\theta }_k}\left(A^{-1}\right)=-A^{-1}\frac{\partial A}{\partial {\theta }_k}{A}^{-1}.$

Differentiating A, element-by-element, with respect to θ_k gives,

$\begin{array}{cc}\begin{array}{c}\frac{\partial a_{\mathrm{ij}}}{\partial {\theta }_k}=\frac{\partial a_{\mathrm{ij}}}{\partial \psi }\frac{\partial \psi }{\partial {\theta }_k}+\frac{d{a}_{\mathrm{ij}}}{d{\theta }_k},\\ =\sum _{l=1}^4\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_l}\frac{\partial {\psi }_l}{\partial {\theta }_k}+\frac{d{a}_{\mathrm{ij}}}{d{\theta }_k},\\ =\sum _{l=1}^4\frac{\partial a_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{d{a}_{\mathrm{ij}}}{d{\theta }_k}.\end{array}& \left(127\right)\end{array}$

In a similar manner we can now write the partial differential of the components of B, C, τ and γ,

$\begin{array}{cc}\frac{d{b}_{\mathrm{ij}}}{d{\theta }_k}=\sum _{l=1}^4\frac{\partial b_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{\partial b_{\mathrm{ij}}}{\partial {\theta }_k},& \left(128\right)\\ \frac{d{c}_{\mathrm{ij}}}{d{\theta }_k}=\sum _{l=1}^4\frac{\partial c_{\mathrm{ij}}}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{\partial c_{\mathrm{ij}}}{\partial {\theta }_k},& \left(129\right)\\ \frac{d{\tau }_i}{d{\theta }_k}=\sum _{l=1}^4\frac{\partial t_i}{\partial {\psi }_l}{J}_{I,\mathrm{lk}}+\frac{\partial t_i}{\partial {\theta }_k},& \left(130\right)\\ \frac{d{\gamma }_i}{d{\theta }_k}=\frac{\partial g_i}{\partial {\theta }_k}.& \left(131\right)\end{array}$

Since A, B and C do not vary with {dot over (θ)}_i, the partial derivative of {umlaut over (θ)} with respect to {dot over (θ)} is given by,

$\begin{array}{cc}\begin{array}{c}\frac{\partial \ddot{\theta }}{\partial {\stackrel{.}{\theta }}_k}=\frac{\partial }{\partial {\stackrel{.}{\theta }}_k}\left[A^{-1}\left(-B{\stackrel{.}{\theta }}_{\mathrm{ii}}-C{\stackrel{.}{\theta }}_{\mathrm{ij}}+\tau +\gamma \right)\right],\\ =A^{-1}\left(-B\frac{\partial {\stackrel{.}{\theta }}_{\mathrm{ii}}}{\partial {\stackrel{.}{\theta }}_k}-C\frac{\partial {\stackrel{.}{\theta }}_{\mathrm{ij}}}{\partial {\stackrel{.}{\theta }}_k}\right).\end{array}& \left(132\right)\end{array}$

The derivatives

$\frac{\partial {\stackrel{.}{\theta }}_{\mathrm{ii}}}{\partial {\stackrel{.}{\theta }}_k}$ and

$\frac{\partial {\stackrel{.}{\theta }}_{\mathrm{ij}}}{\partial {\stackrel{.}{\theta }}_k}$ follow trivially.

The discrete-time state transition matrix can be expressed as followsΦ(t_i,t_i−1)=exp(F(t_i−1)Δt),  (133)Linearized Distribution Matrix

The linearized input distribution matrix G is similarly given by

$\begin{array}{cc}G=\frac{\partial f}{\partial u}=\left(\begin{array}{cc}\frac{\partial \stackrel{.}{\theta }}{\partial T}& \frac{\partial \stackrel{.}{\theta }}{\partial g}\\ \frac{\partial \ddot{\theta }}{\partial T}& \frac{\partial \ddot{\theta }}{\partial g}\end{array}\right)=\left(\begin{array}{cc}0& 0\\ A^{-1}\frac{\partial \tau }{\partial T}& A^{-1}\frac{\partial \gamma }{\partial g}\end{array}\right)& \left(134\right)\end{array}$

The derivatives

$\frac{\partial \tau }{\partial T}$ and

$\frac{\partial \gamma }{\partial g}$ are easily obtained from Eqns. 112 and 108.The discrete-time inputs satisfy

$\Gamma u\left(t\right)={\int }_{t_{i-1}}^{t_i}B\left(t\right)u\left(t\right)dt,$ from which follows the approximation of Γ byΓ=BΔt.

Note that the distribution matrix is not required for application of CDEKF.

Measurement Equations

The available measurements are motor positions and motor armature voltages and currents. We desire to relate these to state variables x^T=(θ,{umlaut over (θ)}). The relationship is governed by the previous constraint equation.

It is assumed there is no direct measurement of motor speeds, however motor speed can be obtained the equations governing the electrical dynamics of the DC drives as follows. In matrix form these electrical dynamics can be written

$\left(\begin{array}{ccc}{L}_{\mathrm{as}}& 0& 0\\ 0& L_{ac}& 0\\ 0& 0& L_{\mathrm{ah}}\end{array}\right)\left(\begin{array}{c}{\mathrm{di}}_{\mathrm{as}}\mathrm{dt}\\ {\mathrm{di}}_{ac}\mathrm{dt}\\ {\mathrm{di}}_{\mathrm{ah}}\mathrm{dt}\end{array}\right)+\left(\begin{array}{ccc}{R}_{\mathrm{as}}& 0& 0\\ 0& R_{ac}& 0\\ 0& 0& R_{\mathrm{ah}}\end{array}\right)\left(\begin{array}{c}{i}_{\mathrm{as}}\\ i_{ac}\\ i_{\mathrm{ah}}\end{array}\right)+\left(\begin{array}{ccc}{K}_{\mathrm{vs}}\left(i_{\mathrm{fs}}\right)& 0& 0\\ 0& K_{\mathrm{vc}}\left(i_{\mathrm{fc}}\right)& 0\\ 0& 0& K_{\mathrm{vh}}\left(i_{\mathrm{fh}}\right)\end{array}\right)\left(\begin{array}{c}d{\theta }_s\mathrm{dt}\\ d{\theta }_c\mathrm{dt}\\ d{\theta }_h\mathrm{dt}\end{array}\right)=\left(\begin{array}{c}{V}_{\mathrm{as}}\\ V_{ac}\\ V_{\mathrm{ah}}\end{array}\right)$ $\mathrm{or}$ $L\frac{d{I}_a}{dt}+R{I}_a+K_v\frac{d\psi }{dt}=V_a.$

This expression can be rearranged in terms of motor speeds to give

$\begin{array}{cc}\frac{d\psi }{dt}=K_v^{-1}\left(V_a-L\frac{d{I}_a}{dt}-{\mathrm{RI}}_a\right)& \left(135\right)\end{array}$

Because the DC drives have significant inductance, the product

$L\frac{d{I}_a}{dt}$ can be significant at times of high rate of change of current. The rate of change of armature current is not directly measured, but an estimate of it can be obtained by high pass filtering currents. Alternatively motor speeds could be estimated using the motor model within a Kalman filter structure.

The state variables {dot over (θ)} are related to

$\frac{d\psi }{dt}$ by

$\begin{array}{cc}\stackrel{.}{\theta }=J_{\mathrm{dr}}{K}_v^{-1}\left(V_a-L\frac{d{I}_a}{dt}-{\mathrm{RI}}_a\right)& \left(136\right)\end{array}$ where J_dr is obtained as previously defined as:

$J_{\mathrm{dr}}=\left(\begin{array}{ccc}-\frac{1}{G_s}& 0& 0\\ 0& \frac{\sin \left({\theta }_2-{\theta }_5\right)}{G_c{\gamma }_1}& \frac{1}{G_h{\gamma }_1}\\ 0& -\frac{1}{G_c}& 0\end{array}\right),\mathrm{with}$ ${\gamma }_1=d_3\cos \left({\theta }_2-{\theta }_5\right)-\left(l_2+l_4\right)\sin \left({\theta }_2-{\theta }_5\right).$ Measurement Noise Covariance

The measurement errors on motor positions, armature voltages, and armature currents can all be considered independent and will be assumed to be zero mean and Gaussian. This imposes a covariance structure on the effective measurements y^T=(θ,{dot over (θ)}) that is determined by Eqns. 135 and 136.

We begin by considering the covariance on errors in the derived motor speed measurements. The constituent terms are measurement errors in armature current, rate of current, and armature voltage. Let these be given by

${\Sigma }_{Ia}=\left(\begin{array}{ccc}{\sigma }_{Ias}^2& 0& 0\\ 0& {\sigma }_{\mathrm{Iac}}^2& 0\\ 0& 0& {\sigma }_{\mathrm{Iah}}^2\end{array}\right)$ ${\Sigma }_{\stackrel{.}{I}a}=\left(\begin{array}{ccc}{\sigma }_{\stackrel{.}{I}\mathrm{as}}^2& 0& 0\\ 0& {\sigma }_{\stackrel{.}{I}ac}^2& 0\\ 0& 0& {\sigma }_{\stackrel{.}{I}\mathrm{ah}}^2\end{array}\right)$ ${\Sigma }_{\mathrm{Va}}=\left(\begin{array}{ccc}{\sigma }_{\mathrm{Vas}}^2& 0& 0\\ 0& {\sigma }_{\mathrm{Vac}}^2& 0\\ 0& 0& {\sigma }_{\mathrm{Vah}}^2\end{array}\right)$

It is well established that a sum of zero mean Gaussian distributed random variables is a zero mean Gaussian distributed variable with covariance equal to the sum of the covariances on the constituent parts. Hence the covariance of estimated motor speed will beτ_{{dot over (ψ)}}=K_v⁻¹(Σ_Va+LΣ_İaL^T+RΣ_İaR^T)K_v^−T  (137)

Because the constituent matrices are all diagonal, Σ_{{dot over (ψ)}} is diagonal.

Assume measurement errors on motor positions are zero mean and Gaussian with covariance given by

${\Sigma }_{\psi }=\left(\begin{array}{ccc}{\sigma }_{\psi s}^2& 0& 0\\ 0& {\sigma }_{\psi c}^2& 0\\ 0& 0& {\sigma }_{\psi h}^2\end{array}\right)$

The measurement error covariance on the derived state measurements y is hence

$\begin{array}{c}R=\left(\begin{array}{cc}{J}_{\mathrm{dr}}& 0\\ 0& J_{\mathrm{dr}}\end{array}\right)\left(\begin{array}{cc}{\Sigma }_{\psi }& 0\\ 0& {\Sigma }_{\stackrel{.}{\psi }}\end{array}\right)\left(\begin{array}{cc}{J}_{\mathrm{dr}}^T& 0\\ 0& J_{\mathrm{dr}}^T\end{array}\right)\\ =\left(\begin{array}{cc}{J}_{\mathrm{dr}}{\Sigma }_{\psi }{J}_{\mathrm{dr}}^T& 0\\ 0& J_{\mathrm{dr}}^T{\Sigma }_{\stackrel{.}{\psi }}{J}_{\mathrm{dr}}^T\end{array}\right)\end{array}$

This matrix is not diagonal. The off-diagonal terms show a coupling between crowd and hoist. Note that because J_dr is configuration dependent, R is configuration dependent.

Implementation Notes

Practical implementation requires the expansion of both the equations of motion and their linearization into a computable form. In an embodiment practical implementation is achieved using the MATLAB symbolic mathematics toolbox. The approach starts with the transformation matrices and constraint equation and expands them symbolically to form the terms in Eqns. 115 and 125. The resulting expressions are then written out as both MATLAB ‘mfiles’ and ‘C-code’ fragments. The MATLAB code can be used for offline development and debugging while the C-code fragments can be formatted from within MATLAB for insertion directly into the code-base of the real-time system implementation.

The expressions for individual element terms resulting from the symbolic expansions are long and cumbersome and contain repeated expressions, suggesting that a computationally more efficient implementation is possible than that produced directly from the symbolic expansion under MATLAB. While a significant computation load, this does not impinge on the capacity to perform the real-time calculations in modern computer systems, even when running at rates greater than 100 Hz and with a large filter bank (50 models).

The realtime implementation is written in entirely in the C-language and runs under the QNX realtime operating system. The hardware platform is based on the PC-104 format with a 1 GHz Pentium processor. The system is populated with data acquisition cards to acquire signals form the electrical drives as described above. The payload estimator is typically configured to run at a rate 100 Hz with a bank of 50 filters. With this configuration, the CPU is loaded to approximately half its capacity. This load includes data acquisition, data logging, and other necessary activities such as the running of algorithms to detect changes in shovel activity, e.g. the beginning and end of the ‘swing-with-full-bucket’ phase of the shovel cycle when the payload calculations are undertaken.

Both the MATLAB and C versions implement the MMAE algorithm detailed earlier. Matrix calculations lie at the heart of all Kalman filter implementations. The MATLAB implementation uses the matrix capabilities provided within the language. The C implementation uses the MESCHACH public domain matrix library. In both implementations integration of the equation of motion, as needed by the MMAE algorithm, are completed using a 4-th order, fixed step, Runge-Kutta, ODE integrator. A comparison with adaptive step solvers during early simulations found the additional complexity and computational cost associated with such solvers was not warranted. The conversion of the continuous dynamics to discrete form, viz Eqn. 133, is done using the matrix exponential functions provided by MATLAB or MESCHACH as appropriate.

Motors are modelled as torque sources whose output is determined by the field and armature currents according to a characteristic such as that shown previously. These characteristics are represented in software using lookup tables with linear interpolation between grid points. Data for the characteristics was generated experimentally by the shovel manufacturer.

Geometric dimensions, inertias, gear ratios, etc as needed have been supplied from the shovel manufacturer. Dimensional information and gear ratios were sourced from manufacturing drawings. Masses for individual components and inertias have been obtained from CAD models of machines.

Recursive Least Squares Payload Estimation Using Static Force and Moment Balance

Observe that the static force and moment balance previously described can be rearranged in terms of the payload mass m_p as follows

$\begin{array}{cc}\left(\begin{array}{c}g\left(\left(c_{\mathrm{zp}}+l_2\right)\cos {\theta }_2+\left(c_{\mathrm{xp}}+d_3\right)\sin {\theta }_2\right)\\ g\cos {\theta }_2\end{array}\right)m_p=\left(\begin{array}{c}\left[\left(G_h\left(l_2+l_4\right)\sin \left({\theta }_5-{\theta }_2\right)+G_h{d}_3\cos \left({\theta }_5-{\theta }_2\right)\right)T_h-g\left(m_2{c}_{x2}+m_3\left(c_{z2}+l_2\right)\right)\cos {\theta }_2-g\left(m_3\left(c_{x3}+d_3\right)-m_2{c}_{z2}\right)\sin {\theta }_2\right]-\\ m_{3}g\cos {\theta }_2+G_c{T}_c+G_h\sin \left({\theta }_5-{\theta }_2\right)T_h\end{array}\right)\mathrm{summarized}\mathrm{as}{a}_k{m}_p=b_k& \left(138\right)\end{array}$

For n measurements Eqn. 138 provides 2 n equations that can be assembled as follows

${\left(\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_n\end{array}\right)}^{m_p}=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$ $\mathrm{or}$ ${\mathrm{am}}_p=b$

Equation. 138 has the least squares solution{circumflex over (m)}_p=(a^Ta)⁻¹a^Tb.  (139)which can be computed recursively via{circumflex over (m)}_p,k={circumflex over (m)}_p,k−1+K_k(b_k−a_k{circumflex over (m)}_p,k−1)  (140)withK_k=P_ka_k P_k=(a_k^Ta_k)⁻¹  (141)and {circumflex over (m)}_p,k providing a least squares estimate of m_p based on k measurements. In Algorithm 78, Eqn. 140 is used as the basis for the recursive least squares algorithm for payload estimation based on static force and moment balances.Algorithm 5: Payload estimation for static force and moment equations by recursive least squares

Initialization: {circumflex over (m)}_p,0=0, k=0

for duration of swing do

wait (timer interrupt)

read (ψ_F,k, T_k)

Solve for θ_F,k from ψ_F,k using Algorithm 3.3.1

Compute payload estimate {circumflex over (m)}_p,k using Eqn. 140, i.e.P_k=(a_k^Ta_k)⁻¹,K_k=P_k⁻¹a_k,{circumflex over (m)}_p,k=m_p,k−1+K_k(b_k−a_k{circumflex over (m)}_p,k−1)k=k+1A comparison of Static and Dynamic Payload Estimation

The methodology used in comparing Algorithms is as follows:

• • The parameters of the model summarized by Eqn. 115 for a P&H4100 shovel were established from data provided by P&H. This includes dimensions, masses, and inertias from the P&H solid geometry (CAD) models developed during design the machine and from which manufacturing drawings for the machine were generated. Properties of the electrical drives including resistances, inductances, torque and speed characteristics where similarly extracted from experimentally generated characterization data provided by P&H.
  • The parameter-populated dynamic model was augmented with proportional-plus-derivative (PD) loops, closed on each motor position. These loops were tuned to give over-damped responses to step inputs on each axis of motion for the configuration when the payload mass is zero.
  • Measurements were made for sixty-three swing motions of a P&H4100 electric mining shovel. Data for four operators was collected to capture a range of ‘operator styles’. The measured variables are those listed in the table below and included motor positions. This data was collected at 60 Hz.
  • Six configurations of the shovel were randomly chosen from this data to serve as controls for comparing the two algorithms when the machine is stationary. Under these conditions the underlying equations of both algorithms are valid and both algorithms are expected to give accurate results.
  • The recorded swing motions and the static control positions were used as inputs to the motor position loops to generate a time-series of motor-torque data for each ‘loaded swing’ trajectory and the ‘machine stationary’ data. For each loaded swing a fictitious payload, drawn from a uniform distribution in the range 0-80,000 kg, was added to the dipper so that the generated torques included a component required to move dipper payload. The objective is, of course, to determine this payload.
  • The computed motor torques are treated as if they were measured inputs. The motor positions and velocities resulting from these applied torques are treated as measured outputs. Estimates of payload are made from this data using the above Algorithms.

Tables 1 to 5 below give the results of applying the MMAE and static payload estimation errors to the synthesized data. The filter bank for the MMAE estimator was based on 30 Kalman filters. The hypothesized payloads distributed uniformly and regularly spaced over the range 0-80,000 kg.

Table 1 shows the results of the control tests. As expected, both algorithms produce estimates accurate to fractions of a percent of dipper capacity. Each row of Table 1 corresponds to holding the command positions constant to the PD loops on the truth model, so that, apart from some initial motion as the shovel model settled to a steady state under the PD control loops, the computed torques are those required to hold the model at a constant position under gravitational loading. The MMAE algorithm shows lower bias, but the standard deviation across the six control configurations is similar for both methods. The errors resulting from the application of the static payload algorithm are due to initial motion as the model settles to a constant position. If the dynamics are allowed to settle before the algorithm is applied, the estimate is made with zero error. The errors arising from the application of Algorithm arise largely because of the more complicated estimation process which sees the estimate generated as a weighted sum of the computed distribution.

TABLE 1
Performance comparison MMAE and static payload estimates for constant inputs.
				MMAE	MMAE		Static	Static
	True	Swing	MMAE	Payload	error as	Static	Payload	error as
Sw.	Payload	time	Payload	Error	% of full	Payload	Error	% of full
No	(kg)	(s)	(kg)	(kg)	scale	(kg)	(kg)	scale
1	41109.65	18	41056.62	−53.04	−0.07	41133.42	23.76	0.03
2	14058.05	17.83	14084.44	26.4	0.03	14005.46	−52.59	−0.07
3	24690.69	19.3	24707.1	16.4	0.02	24609.21	−81.49	−0.1
4	42762.54	18.5	42751.89	−10.64	−0.01	42688.88	−73.66	−0.09
5	75810.42	19.13	75866.74	56.32	0.07	75730.55	−79.87	−0.1
6	13738.21	20.6	13692.9	−45.31	−0.06	13721.3	−16.91	−0.02
			Average	−1.645		Average	−46.7933
			Std. dev.	42.66397		Std. Dev.	42.2472

Tables 2 to 5 compare the results of both algorithms when applied to the data synthesized by the measured motions from the four different operators. The significant points in these tables are the large errors associated with the static payload estimates. The static payload algorithm is seen to systematically underestimate payload. The mean errors range from approximately 12,000 kg to 20,000 kg. The standard deviations on those estimate errors range from approximately 3600 kg to 8000 kg. These range variations suggest that the static payload algorithm is susceptible to operator style. Inspection of the dipper trajectories reveals that is the collected swing motions the dipper is typically lowered after breaking the face before the material is dumped. The bench conditions during the data collection period were such that the operator was required to move the bucket up the face to achieve reasonable fills. The dipper was then lowered to the appropriate level for dumping to the tray of a haul truck. This net downward motion accounts for the underestimated payload estimates. Under different bench conditions, in particular in fine, loose rilling, dig material, the dipper is often filled earlier, and it is common for the operator to extract the dipper low in the face so that a net upward motion is needed to bring the dipper to truck dump height. In such circumstances, the static payload estimator is expected to overestimate payload.

The dynamic payload estimator achieves significantly better estimates of payload across all motion trajectories. The maximum error for the 63 trajectories is 229 kg, the average errors for each operator are all less than 100 kg, and the standard deviation on the estimates for each operator is less than 60 kg.

These main conclusion drawn here is that the inertial and velocity loads associated with motion during are sufficiently large that they need to be accounted for. The results further suggest that a multiple model adaptive estimator incorporating the rigid body dynamics is able to account for these loads and accurately estimate payload. The MMAE estimator overcomes the need for ad hoc approaches and is expected to out-perform them.

The MMAE estimator is expected to do well in this comparative study because it encodes exact knowledge of the dynamics of the system. Recall that according to the methodology, the same models used to generate torques and measurements treated as if they were measured as is used to process these torques and measurement in the MMAE algorithm. Note that both algorithms performed well under static conditions, showing that under the conditions where underlying system equations are known, both approaches work. In a practical implementation, the MMAE estimator will have only approximate knowledge of the machine dynamics due unmodelled dynamics, parameter errors, and transmission losses.

TABLE 2
Performance comparison MMAE and static payload estimates
for constant inputs for swings by Operator A.
				MMAE	MMAE		Static	Static
	True	Swing	MMAE	Payload	error as	Static	Payload	error as
Sw.	Payload	time	Payload	Error	% of full	Payload	Error	% of full
No	(kg)	(s)	(kg)	(kg)	scale	(kg)	(kg)	scale
1	56178.47	14.82	56198.71	20.24	0.03	45741.73	−10436.7	−13.05
2	18113.35	17.85	18250.68	137.33	0.17	10985.72	−7127.63	−8.91
3	39581.29	15.33	39710.96	129.67	0.16	24381.55	−15199.7	−19
4	9975.89	13.37	10030	54.11	0.07	−1156.42	−11132.3	−13.92
5	6711.63	15.05	6782.09	70.46	0.09	−10645.1	−17356.7	−21.7
6	31170.38	15.88	31236.75	66.37	0.08	14756.78	−16413.6	−20.52
7	22178.41	17.27	22304.63	126.21	0.16	6364.41	−15814	−19.77
			Average	86.34143		Average	−13354.4
			Std. dev.	44.95562		Std. dev.	3809.091

TABLE 3
Performance comparison MMAE and static payload estimates
for constant inputs for swings by Operator B.
				MMAE	MMAE		Static	Static
	True	Swing	MMAE	Payload	error as	Static	Payload	error as
Sw.	Payload	time	Payload	Error	% of full	Payload	Error	% of full
No	(kg)	(s)	(kg)	(kg)	scale	(kg)	(kg)	scale
8	29444.26	18.07	29548.19	103.93	0.13	19540.68	−9903.58	−12.38
9	78676.72	13.72	78712.46	35.74	0.04	60387.14	−18289.6	−22.86
10	42830.9	14.95	42892.69	61.79	0.08	29983.69	−12847.2	−16.06
11	61254.31	13.17	61298.49	44.18	0.06	38865.13	−22389.2	−27.99
12	51717.89	13.8	51806.38	88.49	0.11	35668.45	−16049.4	−20.06
13	61371.5	14.83	61503.33	131.83	0.16	40289.47	−21082	−26.35
14	62418.9	15.53	62500.7	81.8	0.1	44351.95	−18067	−22.58
15	65836.97	15.25	65933.67	96.7	0.12	42099.82	−23737.2	−29.67
16	12153.69	15.15	12203.87	50.18	0.06	−4419.99	−16573.7	−20.72
17	50038.15	11.08	50041.19	3.04	0	37470.25	−12567.9	−15.71
18	25174.11	19.47	25265.62	91.51	0.11	7333.54	−17840.6	−22.3
19	27752.31	13.02	27747.51	−4.8	−0.01	12696.53	−15055.8	−18.82
20	73376.26	13.18	73461.55	85.29	0.11	52077.75	−21298.5	−26.62
21	41580.86	17.05	41621.29	40.43	0.05	23183.51	−18397.4	−23
			Average	65.00786		Average	−17435.6
			Std. dev.	38.99654		Std. dev.	3966.105

TABLE 4
Performance comparison MMAE and static payload estimates
for constant inputs for swings by Operator C.
				MMAE	MMAE		Static	Static
	True	Swing	MMAE	Payload	error as	Static	Payload	error as
Sw.	Payload	time	Payload	Error	% of full	Payload	Error	% of full
No	(kg)	(s)	(kg)	(kg)	scale	(kg)	(kg)	scale
22	32093.26	13.33	32172.86	79.6	0.1	9799.14	−22294.1	−27.87
23	48541.52	15.73	48569.93	28.4	0.04	25054.54	−23487	−29.36
24	62833.95	9.75	63044.27	210.32	0.26	28215.5	−34618.5	−43.27
25	74523.76	13.63	74548.74	24.98	0.03	45171.5	−29352.3	−36.69
26	69594.41	13.63	69667.37	72.96	0.09	47424.97	−22169.4	−27.71
27	69323.4	13.13	69405.6	82.19	0.1	43049.91	−26273.5	−32.84
28	53961.61	12.25	54043.03	81.43	0.1	34332.75	−19628.9	−24.54
29	60673.24	13.05	60717.59	44.35	0.06	31857.58	−28815.7	−36.02
30	46551.71	15.7	46601.94	50.23	0.06	20901.25	−25650.5	−32.06
31	31138.65	13.85	31183.57	44.93	0.06	9117.96	−22020.7	−27.53
32	28450.58	13.87	28558.43	107.85	0.13	10079.75	−18370.8	−22.96
33	16018.56	12.85	16061.77	43.21	0.05	−7296.37	−23314.9	−29.14
34	66154.36	13.2	66269.99	115.63	0.14	44984.01	−21170.4	−26.46
35	33272.5	13.45	33407.59	135.09	0.17	4753.15	−28519.4	−35.65
36	37081.21	14.53	37185.47	104.26	0.13	13530.89	−23550.3	−29.44
37	78334.91	13.77	78428.51	93.6	0.12	46680.9	−31654	−39.57
38	10115.05	14.33	10257.79	142.74	0.18	−6758.72	−16873.8	−21.09
39	17009.8	14.75	17024.2	14.4	0.02	−6409.95	−23419.8	−29.27
40	76677.14	17.28	76748.62	71.48	0.09	58582.81	−18094.3	−22.62
41	58998.38	13.47	59090.99	92.61	0.12	34444.32	−24554.1	−30.69
42	32723.17	14.9	32797.34	74.17	0.09	17699.28	−15023.9	−18.78
43	62409.13	14.62	62407.41	−1.73	0	46070.52	−16338.6	−20.42
44	60631.73	13.6	60757.15	125.42	0.16	43974.49	−16657.2	−20.82
45	76547.75	13.73	76607.19	59.44	0.07	54947.58	−21600.2	−27
			Average	79.065		Average	−23060.5
			Std. dev.	47.07695		Std. dev.	5045.305

TABLE 5
Performance comparison MMAE and static payload estimates
for constant inputs for swings by Operator D.
				MMAE	MMAE		Static	Static
	True	Swing	MMAE	Payload	error as	Static	Payload	error as
Sw.	Payload	time	Payload	Error	% of full	Payload	Error	% of full
No	(kg)	(s)	(kg)	(kg)	scale	(kg)	(kg)	scale
46	2246.16	18	2309.85	63.69	0.08	1908.93	−337.23	−0.42
47	25498.83	14.92	25538.71	39.88	0.05	20661.83	−4837	−6.05
48	60556.05	12.8	60588.2	32.15	0.04	53690.35	−6865.7	−8.58
49	19439.07	17.17	19492.26	53.19	0.07	5352.62	−14086.5	−17.61
50	47164.53	14.73	47079.99	−84.53	−0.11	32649.71	−14514.8	−18.14
51	3471.79	16.73	3480.58	8.79	0.01	−2873.06	−6344.85	−7.93
52	76484.27	16.08	76502.62	18.35	0.02	66165.12	−10319.2	−12.9
53	25530.56	13.88	25523.24	−7.33	−0.01	4250.77	−21279.8	−26.6
54	4748.68	11.32	4857.76	109.08	0.14	−15098.6	−19847.3	−24.81
55	35350.2	17.6	35330.6	−19.6	−0.02	23222.3	−12127.9	−15.16
56	73202.92	13.45	73120.45	−82.46	−0.1	53792.76	−19410.2	−24.26
57	45780.21	17.07	45679.72	−100.48	−0.13	22644.07	−23136.1	−28.92
58	9507.13	19.57	9560.29	53.16	0.07	175.14	−9331.99	−11.66
59	45582.45	16.22	45525.29	−57.16	−0.07	18845.73	−26736.7	−33.42
60	20164.19	17.27	20182.64	18.45	0.02	10072.07	−10092.1	−12.62
61	39669.18	16.73	39625.3	−43.88	−0.05	19270.39	−20398.8	−25.5
62	18938.57	15.9	18947.73	9.16	0.01	−13376.4	−32314.9	−40.39
63	38157.9	18.1	38259.99	102.09	0.13	22706.19	−15451.7	−19.31
			Average	6.252778		Average	−14857.4
			Std. dev.	61.57864		Std. dev.	8267.319

The foregoing describes preferred forms of the present invention. Modifications are obvious to those skilled in the art can be made thereto without departing from the scope of the invention.

Claims

1. A method of payload estimation of weight bearing machinery, the method comprising the steps of: (a) creating a series of Kalman filter approximations to the weight bearing machinery, each approximation representing a system model of the weight bearing machinery having distinct levels of payload;(b) determining from said series an overall approximation of the current operating parameters of the weight bearing machinery; and,(c) utilizing the parameters derived from said step (b) to determine an estimated payload of said weight bearing machinery.

2. The method according to claim 1 wherein said weight bearing machinery comprises an electric mining shovel.

3. The method according to claim 1 wherein the inputs to one or more of the Kalman filters include weight bearing machinery drive control signals for controlling movement of the weight bearing machinery or cognate signals leading to direct or indirect measurement of hoist rope and crowd forces.

4. The method according to claim 1 wherein each Kalman filter approximation simulates increasing payloads on said weight bearing machinery.

5. The method according to claim 1 wherein said step (b) further comprises the steps of forming a current model probability for each system model and then forming a weighted sum of the models to determine the overall approximation of the current operating parameters of the weight bearing machinery.

6. The method according to claim 1 wherein said step (b) further comprises the steps of: identifying one of the system models which most closely approximates the current operating parameters; andsetting the identified model as the overall approximation of the current operating parameters of the weight bearing machinery by assigning the identified model a probability larger than the probability assigned to the other system models.

7. A method of payload estimation of weight bearing machinery, the method including the steps of: (a) utilizing a plurality of adaptive predictive models for predicting a future state of a series of time varying parameters associated with a rigid body dynamic equation of the machinery;(b) determining for each model a difference measure between the projected future state and an actual measured state;(c) utilizing said difference measure to adaptively improve model coefficients of each respective said adaptive predictive model; and,(d) utilizing each said adaptive predictive model to produce said payload estimation.

8. The method according to claim 7 wherein each said adaptive predictive model is directly derived from a system model.

9. The method according to claim 7 wherein said model coefficients are improved such that variance of said difference measure is minimized in a mean squared sense.

10. The method according to claim 7 wherein each said adaptive predictive model is independent and has a distinct model.

11. The method according to claim 10 wherein said payload estimation is produced by a weighted sum of parameters predicted from each said adaptive predictive model.

12. The method according to claim 7 wherein each said adaptive predictive model is a Kalman filter and Kalman filter techniques are utilized to adaptively improve said model coefficients.

13. A method of payload estimation of weight bearing machinery, the method including the steps of: utilizing a multi-model approach to build up a statistical inferencing system directed to the operation of the machinery, andderiving therefrom a payload estimation, the method further including the steps of:(a) creating a plurality of system models, each model representing the weight bearing machinery having distinct levels of payload;(b) determining from said plurality of models an overall approximation of the current operating parameters of the weight bearing machinery; andutilizing said overall approximation to produce said payload estimation.

14. The method according to claim 13 wherein the plurality of system models is approximated by a series of Kalman filters.

15. A system for payload estimation of weight bearing machinery, said system comprising: one or more input sensors configured to measure input coefficients from said machinery; anda processor configured to receive said input coefficients from said input sensors; said processor further configured to:evaluate a plurality of adaptive predictive models for predicting a future state of a series of time varying parameters associated with a rigid body dynamic equation of the machinery,determine for each model a difference measure between said future state and said measured input coefficients,utilize said difference measure to adaptively improve model coefficients of each said adaptive predictive model, andutilize each said adaptive predictive model to produce said payload estimation.

16. The system according to claim 15 wherein said processor is configured to simultaneously evaluate a plurality of adaptive predictive models for predicting said future state of said parameters.

17. The system according to claim 16 wherein each said adaptive predictive model is independent and has a distinct model.

18. The system according to claim 16 wherein said processor is configured to utilize each said adaptive predictive model by computing a weighted sum of parameters predicted from each said adaptive predictive model.

19. The system according to claim 15 wherein each said adaptive predictive model is a Kalman filter and said processor is configured to adaptively improve said filter coefficients using Kalman filter techniques.