A method to estimate ground reaction forces (GRFs) in a robot/prosthesis system is presented. The system includes a robot that emulates human hip and thigh motion, along with a powered (active) transfemoral prosthetic leg. We design a continuous-time extended Kalman filter (EKF) and a continuous-time unscented Kalman filter (UKF) to estimate not only the states of the robot/prosthesis system but also the GRFs that act on the foot. It is proven using stochastic Lyapunov functions that the estimation error of the EKF is exponentially bounded if the initial estimation errors and the disturbances are sufficiently small. The performance of the estimators in normal walk, fast walk, and slow walk is studied, when we use four sensors (hip displacement, thigh, knee, and ankle angles), three sensors (thigh, knee, and ankle angles), and two sensors (knee and ankle angles). Simulation results show that when using four sensors, the average root-mean-square (RMS) estimation error of the EKF is 0.0020 rad for the joint angles and 11.85 N for the GRFs. The respective numbers for the UKF are 0.0016 rad and 7.98 N, which are 20% and 33% lower than those of the EKF.

Introduction

Recent advances in microelectronics and robotic technologies have enabled the development of powered prosthetic legs [1] that can help amputees walk up stairs and slopes because of the legs' net power contribution to gait; they are also able to adapt their behavior to various environmental conditions. The company Össur has a lower limb prosthesis called the Power Knee, which is the first commercially available knee to generate power during the gait cycle [2]. Recent developments in powered knee and ankle prostheses are reported in Refs. [38].

Bulky load cells are often employed in robots and prosthetic legs to capture external forces (GRFs) and moments during walking [9]. These data are used as feedback measurements to control the robot or prosthesis with force or impedance controls. Even when not used for feedback, force sensing is important for monitoring and evaluation of prosthesis performance and safety. However, there are several drawbacks to the use of load cells: (i) load cells are expensive; (ii) a 250-lbf load cell weighs about 1 lb, with a length of about 3 in, and thus does not easily fit in a prosthetic leg; (iii) load cell measurements tend to drift, need to be frequently offset, and are noisy and need significant signal conditioning; (iv) load cells can get damaged easily from overloading or off-axis loading; and (v) load cells consume electrical power, which is an important consideration in prosthetics.

The above-mentioned problems do not arise with angle sensors because high-resolution encoders are accurate, reliable, and inexpensive. Nevertheless, they do not measure velocity, and velocity calculated by numerical differentiation is challenging because of the difficult compromise between noise rejection and bandwidth.

There have been several methods to reduce the number of force sensors in robotics. In Ref. [10], a robot compliance controller based on a disturbance observer is presented, where the disturbance observer is used to estimate the external reaction forces. In Ref. [11], a method for force estimation of the end-effector of a Selective Compliance Assembly Robot Arm (SCARA) robot is presented, where servomotor currents and position information are used to estimate forces. In Ref. [12], the external force on the end-effector of a four degrees-of-freedom (DOF) robot manipulator is estimated with a combination of time delay estimation and input estimation. In Ref. [13], a model-based observer is used to estimate the external forces acting on a rigid body. All the aforementioned techniques depend on either the accuracy of the robot model or servomotor current data. Thus, if the accuracy of the robot model or motor current measurement degrades, force estimation can deviate from true force. In Ref. [12], the force estimation does not require a precise robot dynamic model, but the accuracy of estimated force may degrade in the presence of noise, since the external force is considered as an unknown input.

In this paper, we treat GRF as an unknown input with known dynamic properties and known bounds, and we use an EKF to estimate GRFs along with the states of the robot/prosthesis system. Although EKF is efficient in many applications, it has two important potential drawbacks. First, the derivation of the Jacobian matrix for linearization can be complex and can cause numerical implementation difficulties. Second, linearization can lead to cumulative errors which may affect the accuracy of the state estimation. To overcome these limitations, other nonlinear estimators could be used, such as the UKF or particle filter, where the state estimations are obtained without the need for derivatives and Jacobian calculations [1416].

In this research, a powered (active) prosthetic leg is considered for transfemoral amputees. The prosthetic leg attached to a robotic hip/thigh emulator. The combined system includes four degrees-of-freedom: vertical hip displacement, thigh angle, knee angle, and ankle angle. This paper is an extended version of Ref. [17], where we designed an EKF for the online estimation of joint coordinates, velocities, and GRFs.

In this paper, we compare EKF and UKF performance in our 4DOF robot/prosthesis system. We study the EKF because it is the most commonly used nonlinear state estimator due to its versatility and simplicity of implementation. We study the UKF because it generally provides better performance than the EKF. One of the questions we study in this paper is whether the improved performance of the UKF relative to the EKF is worth the increased computational effort.

We also study the effect of various sensor sets on estimation error at different walking speeds. We start with the four main states of the robot/prosthesis system, which are hip displacement, thigh, knee, and ankle angles, as system observations to get the baseline estimation accuracy. Then, we investigate how much estimation accuracy degrades if fewer measurements are used, and we find a tradeoff between the estimation error accuracy and the number of measurements. In general, we prefer to use fewer measurements because fewer sensors mean a simpler system, less complexity, and lower cost.

Finally, we consider the stability of the EKF and boundedness of the estimation error, which are very important for prosthesis control. The stability of the EKF under various values of initial estimation error and noise covariances is explored. We analytically show that the estimation error remains bounded under certain conditions in our 4DOF robot/prosthesis system and confirm the analysis with simulation results.

This research involves the mechanical integration of a prosthesis and a test robot so that we can test the prosthesis without human trials. We will eventually need to use the Kalman filter on a prosthesis that is attached to an amputee. We develop the Kalman filter for the robot/prosthesis system [18] in this paper for later implementation on an amputee/prosthesis system.

The paper is organized as follows: In Sec. 2, the model of the robotic system and prosthetic leg is presented. In Sec. 3, the EKF for state estimation and GRF estimation is discussed. In Sec. 4, the EKF is analyzed by mathematically deriving its stability conditions. Section 5 compares performance between the EKF and UKF when different measurement sensors are used; also the convergence of the EKF is tested with different initial estimation errors and noise magnitudes. Section 6 concludes the paper and suggests future research.

System Model

Robotic testing of transfemoral prostheses is presented here, where motion is limited to the sagittal plane and transverse motion is not considered. Typically, only the sagittal plane is considered in transfemoral prosthesis research. Although the transverse plane is ignored here, sagittal motion captures the essential dynamics of human walking [19]. The model of the test robot/prosthesis is based on the standard robotic framework. Figure 1 shows a diagram of the hip robot and prosthesis combination [20,21]. A general dynamic model for the system is given as follows:

 
D(q)q̈+C(q,q̇)q̇+B(q,q̇)+JeTFe+g(q)=u
(1)

where q=[q1,q2,q3,q4]T is the vector of joint displacements (q1 is vertical hip displacement, q2 is thigh angle, q3 is knee angle, and q4 is ankle angle), D(q) is the inertia matrix, C(q,q̇) is a matrix accounting for centripetal and Coriolis effects, B(q,q̇) is a nonlinear damping vector, Je is the kinematic Jacobian relative to the point of application of the external forces Fe, g(q) is the gravity vector, and u is the four-element vector of control signals [22]. The kinematic and dynamic models of the robot/prosthesis combination are given in Refs. [2325], where a mixed tracking/impedance controller based on passivity methods is designed [22]. The control signals consist of hip force, and thigh, knee, and ankle torques. As Fig. 1 shows, a triangular foot with two points of ground contact is assumed. Horizontal and vertical GRFs are applied to contact points at the toe and heel. The GRFs are denoted as Fxh, Fzh, Fxt, Fzt, which represent the horizontal and vertical GRFs at the heel and toe. Thus, the external force vector Fe in Eq. (1) comprises these four GRFs.

We assume the robot walks along the x-axis. A treadmill provides the walking surface of the robot/prosthesis system. The belt stiffness is modeled to calculate reaction forces during contact between the heel and toe with the treadmill [22,23]. The GRFs are entirely determined by kinematics and are given as follows [23]: 
zh=aHsin(q2+q3+q4+(π2+cos1(ahaH)))+l3sin(q2+q3)+l2sin(q2)+q1
(2)
 
zt=aTsin(q2+q3+q4+(π2cos1(ahaT)))+l3sin(q2+q3)+l2sin(q2)+q1
(3)
 
Fzh=kb(zhsz)(1+sign(zhsz)2)
(4)
 
Fzt=kb(ztsz)(1+sign(ztsz)2)
(5)
 
Fxh=βFzh
(6)
 
Fxt=βFzt
(7)

where kb is the belt stiffness, sz is the treadmill standoff (the distance between the origin of the world coordinate system and the belt when the leg is fully extended), l2 and l3 are the lengths of link 2 (thigh) and link 3 (shank), respectively, aH and aT are the distance from the ankle joint to the heel and the toe, respectively, ah is the height of the ankle joint above the sole of the foot, and β is the friction coefficient between the treadmill belt and the foot. The vertical positions of the toe and heel in the world coordinate system are shown in Fig. 1 as zt and zh, respectively. We thus have four states for the positions and four states for the velocities of the joint displacements.

The main objective of this research is to estimate the GRFs which comprise Fe in Eq. (1). We, therefore, augment the external forces to the state vector. The augmented state vector includes the eight original states of the robot and the four GRFs and is given as follows: 
x=[q1,q2,q3,q4,q1.,q2.,q3.,q4.,Fxt,Fzt,Fxh,Fzh]T
(8)

It should be noted that the forces are not states of the original system model, but are augmented here to the state vector, since we need to estimate them with a state estimator. The 12-element vector of Eq. (8) will be estimated by the state estimator.

Equations (2)(7) are a preliminary approximation to GRF and are accurate only if the amputee (or robot in our case) walks in a highly controlled environment with a known walking surface stiffness. In future work, it will be important to study the robustness of the Kalman filter to modeling errors and particularly to errors in these GRF equations.

We should mention that there are at least two other possible approaches to obtain GRF estimates. In the first alternative approach, GRF could be considered as an unknown input with no prior assumptions about its model. In that case, the EKF would have difficulty estimating GRF, since its effect on the system would be indistinguishable from process noise. In the second alternative approach, state estimation of the robot/prosthesis system could be performed by a nonlinear observer, such as the EKF, and then the estimated states of the robot could be substituted into the GRF Eqs. (2)(7) to obtain the GRF estimates. However, this approach would be less flexible than the approach that we are using. Augmenting the GRF states onto the original state as in Eq. (8) allows us additional tuning flexibility in the artificial process noise that is included in the augmented states, as we will see in Sec. 5.

Extended Kalman Filtering for Robot/Prosthesis State Estimation

The Kalman filter applies directly only to linear systems. However, we can linearize a nonlinear system and then use linear estimation techniques. Over the last few decades, the EKF has become one of the most popular estimation techniques in nonlinear systems. The EKF applies the standard linear Kalman filter methodology to a linearization of the nonlinear system [26]. According to the continuous time-EKF equations [27], the system and measurement equations are assumed as 
ẋ=f(x,u,t)+G(t)w(t)
(9)
 
y=h(x,t)+D(t)v(t)
(10)
where the state x(t), the input u(t), and the output y(t) are in q,p and m, respectively. Moreover, the noise terms w(t) and v(t) are in l and k and are uncorrelated, zero-mean white noise processes with identity covariances. G(t) and D(t) are time-varying matrices of size q × l and m × k. If the nonlinear functions f and h are sufficiently smooth in x, Taylor series expansions can be performed. The following Jacobian matrices are used to linearize the system: 
A(t)=fx|x̂
(11)
 
C(t)=hx|x̂
(12)
The initial state x(0) is a random vector with covariance P(0) 
x̂(0)=E[x(0)]P(0)=E[(x(0)x̂(0))(x(0)x̂(0))T]
and x(0) is assumed to be uncorrelated with w(t) and v(t). The EKF equations are given as 
x̂̇=f(x̂,u,t)+K[yh(x̂,t)]
(13)
 
K(t)=P(t)CT(t)R1(t)
(14)
 
Ṗ=A(t)P(t)+P(t)AT(t)+Q(t)P(t)CT(t)R1(t)C(t)P(t)
(15)
where K(t) and P(t) are the Kalman gain and estimation error covariance matrix. Covariance matrices Q(t) and R(t) are q × q and m × m, respectively, and are defined as follows: 
Q(t)=G(t)GT(t)
(16)
 
R(t)=D(t)DT(t)
(17)
The state-space equations of the robotic/prosthetic system can be obtained from Eq. (1) and the GRF Eqs. (2)(7) must be differentiated to obtain the state-space model of the GRF. However, discontinuities of Eqs. (4) and (5) can be problematic. There are two ways to handle this problem. First, we can use a smooth approximation of the sign functions using hyperbolic tangent functions; but, this approach requires a lot of computational effort for the Jacobian calculation and often results in EKF divergence [16]. In the second approach, which is used in this paper, the sign function is not considered when we take the derivative of Eqs. (2)(7) to include the GRFs in the state-space model. The sign function is re-introduced in the computed Jacobian matrix after taking the derivative. The state-space equations of the GRFs are given as follows: 
Ḟzh=kbzḣ(1+sign(zhsz)2)
(18)
 
Ḟzt=kbzt.(1+sign(ztsz)2)
(19)
 
Ḟxh=βḞzh
(20)
 
Ḟxt=βḞzt
(21)

where zḣ,zṫ represent the derivatives of Eqs. (2) and (3).

The Stability of the Extended Kalman Filter

The analysis of the stability properties of continuous-time extended Kalman filters is complex and has been treated only for especial cases. Safonov and Athans [28] considered the stability of the constant and modified gain EKF. The boundedness of the EKF estimates were investigated when used as a parameter estimator for linear systems [29,30]. Reif et al. [31,32] determined that the estimation error in both the discrete and continuous-time EKF is exponentially bounded under certain conditions. They showed that the estimation error remains bounded if the initial estimation error and disturbances are sufficiently small and the nonlinear system satisfies a detectability rank condition [33,34]. Here, we extend Reif's work [31] to obtain the stability conditions of the state estimator of the robot/prosthesis system. The main purpose of this section is to derive the stability of the estimator in terms of initial condition errors and disturbances.

The nonlinear functions f and h from Eqs. (9) and (10) can be expanded to first-order as follows: 
f(x,u,t)f(x̂,u,t)=A(t)[x(t)x̂(t)]+φ[x,x̂,u,t]
(22)
 
h(x,u)h(x̂,u)=C(t)[x(t)x̂(t)]+X[x,x̂,u]
(23)
where φ and X are the remaining terms of the Taylor series. The estimation error ξ(t)=x(t)x̂(t) can be written recursively as follows: 
dξ(t)=[A(t)K(t)C(t)]ξ(t)+ρ(t)+σ(t)(dw(t)dv(t))
(24)
where σ(t) and ρ(t) are defined by 
σ(t)=[G(t)K(t)D(t)]
(25)
 
ρ(t)=φ[x,x̂,u,t]K(t)X[x,x̂,t]
(26)

In stability theory, there are two commonly used methods to prove boundedness: supermartingales for stochastic differential equations, and Lyapunov functions for deterministic differential equations. Supermartingales can be regarded as the stochastic equivalent of Lyapunov functions in some cases [35]. In order to analyze the error dynamics of the EKF, we first present some preliminary results [31,35]. Note that . represents the spectral norm of a matrix or the Euclidian norm of a vector.

Lemma 4.1. Suppose there is a stochastic processV[ξ(t),t]withξ(t)qgoverned by Eq.(24)and there exist real numbersv1,v2,μ,ϑ0such that 
v1ξ(t)2V[ξ(t),t]v2ξ(t)2
(27)
and 
LV[ξ(t),t]ϑV[ξ(t),t]+μ
(28)
where LV[ξ(t),t] is the differential generator for the stochastic process and is defined as 
LV[ξ(t),t]=Vt(ξ,t)+Vξ(ξ,t)([A(t)K(t)C(t)]ξ(t)+ρ(t))+12i=1qj=1q2V(ξ,t)ξiξj[σ(t)σT(t)]ij
(29)
Then, for every t0 the estimation error, ξ(t) is exponentially bounded in mean square [36,37] 
E{ξ(t)2}v2v1ξ(0)2eϑt+μv1ϑ
(30)

Proof of Lemma 4.1. See Appendix  A.

Theorem 4.1. Consider the continuous-time nonlinear dynamic system Eqs.(9)and(10)and the extended Kalman filter Eqs.(13)(15). Suppose that the following conditions hold:

  • (i)
    There exist p1, p2, c1, c2, q1, q2, r1, r2 such that the following bounds hold for every t0: 
    p1IP(t)p2I
    (31)
     
    c1C(t)c2
    (32)
     
    q1IQ(t)q2I
    (33)
     
    r1IR(t)r2I
    (34)
  • (ii)
    There exist εφ, εx, Kφ,Kx0 such that the nonlinear functions in Eq. (26) are bounded by 
    φ(x,x̂,u)Kφxx̂2
    (35)
     
    X(x,x̂)Kxxx̂2
    (36)
    for x, x̂q and up with xx̂εφ and xx̂εx, respectively.
Then, the estimation error ξ(t) is exponentially bounded in mean square, if the initial estimation error satisfies 
ξ(0)ε
(37)
where ε is a small positive number. Note that the covariance matrices of the noise terms are bounded as in Eqs. (33) and (34), if there is a real number δ0 such that the following bounds are satisfied for every t0: 
G(t)GT(t)δI
(38)
 
D(t)DT(t)δI
(39)

Proof of Theorem 4.1. See Appendix  B.

To calculate Kφ and Kx, we can consider a compact subset Nq; the constants Kφ, Kx can be obtained by 
Kφ=maxisupxNHessfi(x,u)
(40)
 
Kx=maxisupxNHesshi(x)
(41)

If f and h are twice differentiable with respect to x for every xN, the spectral norm of the Hessian matrices of fi and hi are bounded, where fi and hi are the components of f and h, respectively (see Refs. [31] and [38, Chap. 7]). In the given robot/prosthesis system, the Hessian matrix of f is not a function of control input u.

Note that the nonlinear system has to be uniformly detectable in order to satisfy the boundedness condition of Eq. (31) [39,40]. The detectability of the nonlinear system results in the boundedness of the estimation error bounds for the solution P(t) of the Riccati differential equation (15) based on Eq. (31). We present the following definition for the detectability of nonlinear stochastic systems [39].

Definition 4.1. The pair[(f/x)(x,u),(h/x)(x)],xq,upis uniformly detectable if there exists a bounded matrix valued functionΛ(x)such that 
ηT[fx(x,u)+Λ(x)hx(x)]η|λm|η2
(42)

holds for all x, ηq, where λm is the maximum eigenvalue of [(f/x)(x,u)+Λ(x)(h/x)(x)].

Values for ε and δ in Eqs. (B21) and (B24) are obtained in Appendix  B. In Sec. 5.2, the stability of the EKF for the state estimation of the robot/prosthesis system is considered via simulation and related to Theorem 4.1.

Simulation Results

The performance of the robot/prosthesis system during one step of normal walking, which is approximately 1 s, is considered. The reference data are provided by able-bodied research participants at the Motion Studies Laboratory of the Cleveland Department of Veterans Affairs Medical Center (VAMC) [41]. In the prosthetic leg, we have 12 states that are going to be estimated so the A matrix in Eq. (11) is 12 × 12. For our first simulation, we use four states as the measurements: vertical hip position, thigh angle, knee angle, and ankle angle. Therefore, the dimension of the measurement matrix C in Eq. (12) is 4 × 12. The initial value of the state vector and estimate, and the diagonal covariance matrices Q and R for the process and measurement noise are given as 
x(0)=[0.019,1.130,0.090,2.246,0.093,0.770,1.410,2.890,0,0,0,0]T
 
x̂(0)=[0.040,0.901,0.292,2.440,0.293,0.430,1.112,2.620,0,0,0,0]T
 
Q=Qd/Δt=diag{0.0005,0.0005,0.0005,0.0005,0.0002,0.0002,0.0002,0.0002,0,0,0,0}/Δt
 
R=RdΔt=diag{103,103,103,103}Δt

where x̂0 differs from x(0) due to random initial state estimation errors, Δt denotes the discretization step size with Δt=0.5 ms, Q and R represent continuous-time noise covariances, and Qd and Rd represent discrete-time noise covariances [27, Chap. 8]. The process noise represents small values of unmodeled dynamics and parameter uncertainties. The Q matrix is chosen based on prior experience with the accuracy of system dynamic modeling; the last four elements of the diagonal matrix Q are zero since the GRF model is assumed to be perfectly known. The R matrix is usually straightforward to determine on the basis of our knowledge of the accuracy of the measurement system.

The results for the state estimation of the robot/prosthesis system are shown in Figs. 24. Although significant initial estimation errors are imposed on the joint displacements and velocities, the EKF converges to the true states quickly. After toe-off there is no foot contact and the GRF is zero.

Reduced Sensor Sets and Unscented Kalman Filtering.

In a real-world scenario, it would not be practical to use a hip displacement sensor on an amputee because of its invasiveness. The other sensors—thigh angle, knee angle, and ankle angle—can easily be used in the real world because they can be mounted on the prosthesis. In this paper, we use a hip displacement sensor to provide a baseline scenario for the reduced sensor sets, which we study later in this section.

More observations result in more accurate state estimation. However, our goal is to use the fewest possible number of sensors for estimation, which result in a reduction of complexity and cost in the prosthesis hardware. Here, we remove the hip displacement sensor and use only three measurements: the thigh, knee, and ankle angles. For our final test in this section, we will use only two measurements: the knee and ankle angles. Moreover, we will also compare the performance of the EKF and UKF for the state and GRF estimation. We use the same initial values for state and estimated state as used in the EKF for normal walking. We also test the performance of the filters with different gait speeds, and we will see that no extra filter tuning is needed for fast and slow walking.

The original discrete-time form of the UKF has been widely used for state estimation of discrete-time systems. However, the discrete-time UKF cannot be directly applied to continuous-time filtering problems, in which the process and measurement equations are modeled as continuous-time stochastic systems. In Ref. [15], the continuous-time UKF was derived from the discrete-time UKF in matrix form, and this is the algorithm that we implement here. The tuning parameters in the unscented transforms were chosen to be α=0.9, β = 2, and κ=3.

Table 1 compares the accuracy of EKF state estimation in terms of the RMSE in three different gait modes (fast walk, normal walk, slow walk). The EKF works well for different walking speeds with a reduced measurement set. It can be seen, as expected, that for all three walking speeds, the estimation errors with four measurements are smaller than with three or two measurements.

Table 2 shows that the UKF achieves smaller RMS estimation errors than the EKF for almost all measurement sets and walking speeds. Table 3 summarizes the average performance of the EKF and UKF. We see that the UKF achieves 30% improvement in the average RMSE of GRFs with four measurements. Although the UKF performs better than the EKF, its computational effort is between two and three times that of the EKF in terms of the number of multiplications and additions [15].

Stability of the Extended Kalman Filter.

It has been shown that the estimation error in the continuous-time Kalman filter remains bounded under certain conditions. Small initial estimation errors and small noise magnitudes are required conditions to obtain error bounds. In this section, we study the stability of the state estimates of the robot/prosthesis system. First, we need to check that the Jacobian matrices in the robot/prosthesis system satisfy the uniform detectability condition of Definition 4.1. To find an appropriate Λ(x), we assume that [(f/x)(x,u)+Λ(x)(h/x)(x)=Iq×q] and solve for Λ(x)=(I(f/x))(h/x)T((h/x)(h/x)T)1. Note that with our measurement system, ((h/x)(h/x)T)1 is an identity matrix and is thus invertible. Therefore, the uniform detectability condition of Definition 4.1 is satisfied in the robotic/prosthesis system. To test the stability of the robotic/prosthesis EKF, the initial value x̂(0) and the measurement and process noise covariances R, Q are chosen as shown in Table 4, where

 
x(0)=[0.019,1.130,0.090,2.246,0.093,0.770,1.410,2.890,0,0,0,0]T
 
x̂good(0)=[0.027,1.091,0.112,2.350,0.193,0.539,1.3122.797,0,0,0,0]T
 
x̂poor(0)=[0.181,1.991,0.712,2.714,0.393,1.291,2.1123.440,0,0,0,0]T

x̂good(0) is a reasonably accurate initial condition for the state estimate, and we will see that it is good enough to ensure convergence of the EKF. x̂bad(0) is a significantly worse initial condition for the state estimate, and we will see that it is not good enough to ensure convergence of the EKF. The R and Q values in Table 4 have been chosen to demonstrate conditions that, respectively, ensure, or do not ensure, EKF convergence. These values for x̂good(0),x̂bad(0), and R and Q were chosen by trial and error to demonstrate their effect on EKF convergence.

The simulation results are presented in Figs. 57, where the unknown state x6(t) in Eq. (8) (angular velocity of thigh), the estimated state x̂6(t), and its estimation error ξ6(t) are plotted. We can see in Fig. 5 that for small initial estimation error and small noise, the estimation error is bounded. However, Figs. 6 and 7 show that for large initial estimation errors or large noise magnitudes, the estimation error is no longer bounded.

Extensive simulation tests show that boundedness of the estimation error is obtained in the numerical simulation if ξ(0)ε=0.35,G(t)GT(t)δ=2×102 and D(t)DT(t)δ=2×107. The theoretical results for ε and δ via Eqs. (B21) and (B24) from Appendix  B yield smaller bounds for stability: ε6.5×105 and δ8.7×1013 with N={xq,x1000}. Since the estimation error does not diverge in practice even with larger initial condition errors and noise terms than those given by the theorem, we conclude that the bounds are very conservative in this system. However, we note that it is not possible to test all possible conditions that exceed the theoretical error bounds; therefore, we naturally expect the theoretical results to yield more conservative stability bounds than the simulation results.

We want the state to be estimated as accurately as possible for eventual implementation in a state-feedback controller. The controller that we used in this research is robust to estimation errors [22], although more accurate state estimates are always desirable in order to reduce controller errors. We have not explored the relationship between estimation error and controller error in this paper but leave it as an important area for future research.

Conclusion and Future Work

We designed an EKF and an UKF to estimate not only the states of a robot/prosthesis system but also the external forces acting on the prosthetic foot. This approach removes the need for heavy and bulky load cells that are otherwise needed for GRF estimation. We achieved satisfactory estimation errors in various gait speeds for the robot/prosthesis system using four, three, and two measurements. The average RMS estimation errors of the EKF for the thigh, knee, and ankle angles in normal walking with four measurements is 0.0033 rad, in comparison with that the UKF which is 0.0020 rad. Although the UKF outperforms the EKF, it requires more computational effort than the EKF, which will be a consideration for real-time implementation.

We proved mathematically that the estimation error in the EKF is exponentially bounded if the initial estimation errors and disturbances are sufficiently small and if the nonlinear system satisfies a detectability rank condition. In simulation tests, we verified that the estimation errors remained bounded for small initial estimation errors and small disturbances. However, the filter is unstable if the initial estimation errors or disturbances are too large.

As far as we know, the present research is the first time that GRF has been estimated for prostheses using state estimation techniques. As we noted in the introduction, some other methods have been used for external force estimation in robotics. For instance, [12] achieved a 1.72% estimation error for a 4DOF robot, compared to our results, which show a 1.85% error when using the EKF with four sensors during normal walking. Our results are, therefore, quantitatively similar to Ref. [12], although the comparison may not be fair since the robotic systems in the two approaches are much different. Other publications in the area of external force estimation do not include quantitative results [10,11,13].

Future work will include increasing and verifying the robustness of the filters, and experimental implementation and verification of the EKF and UKF, first on robotic hardware [21] and then in human trials.

  • • NSF Grant No. 1344954.

Appendix A

Proof of Lemma 4.1. Let ξ(t) be the unique solution to the stochastic differential in Eq. (24) and let V[ξ(t),t] be a stochastic process which is a non-negative and twice differentiable. Suppose that V[ξ(t),t] satisfies the conditions in Eqs. (27) and (28). By applying Ito's formula [42] we can derive 
dV[ξ(t),t]=LV[ξ(t),t]+V(ξ,t)tσ(t)dw̃(t)
(A1)
where LV[ξ(t),t] is the differential generator given in Eq. (29) and w̃(t) is a Wiener process [35]. We take the integral of both sides of Eq. (A1) 
V[ξ(t),t]V[ξ(0),0]=0tLV[ξ(s),s]ds+0tV(ξ,s)tσ(s)dw̃(s)
(A2)
Since w̃(t) is a Wiener process (Gaussian noise with zero mean), the expectation of the last integral in the preceding equation is zero because of the properties of Ito's formula [28]. Therefore 
E{V[ξ(t),t]}=E{V[ξ(0),0]}+0tE{LV[ξ(s),s]}ds
(A3)
Equation (A3) implies the continuity of E{V[ξ(t),t]}, so differentiating and using Eq. (28) yields 
dE{V[ξ(s),s]}ds=E{LV[ξ(s),s]}ϑE{V[ξ(s),s]}+μ
(A4)
Using Ito's formula, we can write Eq. (A4) as 
dexp(ϑs)E{V[ξ(s),s]}dsμexp(ϑs)
(A5)
Integrating both sides of Eq. (A5) from 0 to t yields 
exp(ϑt)E{V[ξ(t),t]}E{V[ξ(0),0]}μϑ(exp(ϑt)1)
(A6)
Multiplying exp(ϑt) on both sides of Eq. (A6) results in 
E{V[ξ(t),t]}exp(ϑt)E{V[ξ(0),0]}+μϑ(1exp(ϑt))exp(ϑt)E{V[ξ(0),0]}+μϑ
(A7)

Finally, applying the expectation operator to Eq. (27) and substituting the result into Eq. (A7) result in the conclusion that ξ(t) is exponentially bounded in the mean square sense as shown in Eq. (30).

Appendix B

Proof of Theorem 4.1. To prove this theorem, we need to construct a stochastic Lyapunov function to satisfy the conditions of Lemma 4.1. The stochastic Lyapunov function is chosen as follows: 
V[ξ(t),t]=ξT(t)P1(t)ξ(t)
(B1)
For notational convenience, we denote P1(t) as π(t) where P(t) is positive definite with probability one. Therefore, from Eq. (31), it follows that 
1p2Iπ(t)1p1I
(B2)
so that Eq. (B1) we can write 
1p2ξ(t)2V[ξ(t),t]1p1ξ(t)2
(B3)
The first condition of Lemma 4.1, which is Eq. (27), is satisfied. For the second condition, we need to construct the differential generator LV[ξ(t),t] of the stochastic process in Eq. (B1). The differential generator can be defined by Eq. (29). The second part of Eq. (29) can be written as follows [31]: 
i=1qj=1q2V(ξ,t)ξiξj[σ(t)σT(t)]ij=tr(σ(t)σT(t)HessV[ξ(t),t])
(B4)
Therefore, from Eqs. (24), (25), (B1), and (B4), we can write the differential generator as 
LV[ξ(t),t]=ξT(t)dπ(t)dtξ(t)+2ξT(t)π(t)([A(t)K(t)C(t)]ξ(t))+2ξT(t)π(t)ρ(t)+tr([G(t)GT(t)+K(t)D(t)DT(t)KT(t)]π(t))
(B5)
From Eqs. (14), (31), (32), and (34) we have 
K(t)=P(t)CT(t)R1(t)P(t)CT(t)R1(t)p2c2r1
(B6)
We define d=p2c2/r1 which allows us to write Eq. (B6) as 
K(t)d
(B7)
Based on Eqs. (38) and (39), one can obtain 
tr[G(t)GT(t)]δtr[I]=qδ
(B8)
 
tr[D(t)DT(t)]δtr[I]=mδ
(B9)
where q and m are the number of rows of G(t) and D(t), respectively. Moreover, by using Eqs. (31) and (B6)(B9), we obtain 
tr([G(t)GT(t)+K(t)D(t)DT(t)KT(t)]π(t))=tr([G(t)GT(t)π(t)])+tr([D(t)DT(t)K(t)KT(t)π(t)])1p1tr[G(t)GT(t)]+d2p1tr[D(t)DT(t)]Kdδ
(B10)
where Kd is defined as 
Kd=q+d2mp1
(B11)
By substituting for ρ(t) from Eq. (26) into the fourth term of Eq. (B5), it can be shown that 
2ξT(t)π(t)ρ(t)=2ξT(t)π(t)(φ[x,x̂,u,t]K(t)X[x,x̂,t])2ξT(t)π(t)φ[x,x̂,u,t]+2ξT(t)π(t)K(t)X[x,x̂,t]
(B12)
We can use the fact that π(t)P(t)=I and substitute Eq. (14) into Eq. (B12) to obtain 
2ξT(t)π(t)ρ(t)2ξT(t)π(t)φ[x,x̂,u,t]+2ξT(t)CT(t)R1(t)X[x,x̂,t]
(B13)
From Eqs. (31)(36), it follows that 
2ξT(t)π(t)ρ(t)2ξ(t)Kφp1ξ(t)2+2ξ(t)c2Kxr1ξ(t)2
(B14)
We simplify Eq. (B14) to obtain 
2ξT(t)π(t)ρ(t)Knξ(t)3
(B15)
where Kn=(2Kφ/p1)+(2c2Kx/r1). Substituting Eqs. (14), (B10), and (B15) into Eq. (B5) yields 
LV[ξ(t),t]ξT(t){dπdt+AT(t)π(t)+π(t)A(t)2CT(t)R1(t)C(t)}ξ(t)+Knξ(t)3+Kdδ
(B16)
We can find dπ(t) using π(t)P(t)=I 
dπ(t)=π(t)dP(t)π(t)
(B17)
where dP(t) is calculated by the Riccati differential equation which is given in Eq. (15). So, if we substitute Eq. (B17) into Eq. (B16), we obtain 
LV[ξ(t),t]ξT(t){π(t)Q(t)π(t)+CT(t)R1(t)C(t)}ξ(t)+Knξ(t)3+Kdδ
(B18)
We can use Eqs. (31)(34) to write Eq. (B18) as 
LV[ξ(t),t][q1p22+c12r2]ξ(t)2+Knξ(t)3+Kdδ
(B19)
We can also use the upper bound of Eq. (B3) to obtain 
LV[ξ(t),t](q1r2+c12p22p22r2Knξ(t)p22r2)p1V[ξ(t),t]+Kdδ
(B20)
Now, we just need to show that ((q1r2+c12p22p22r2Knξ(t))/p22r2)p10. Since the real numbers p1,p2,c1,r2 are positive, we require ξ(t)((q1r2+c12p22)/p22r2K2) to satisfy this condition. Now, according to Theorem 4.1, we can obtain the initial estimation error bound ε as 
ε=min(εφ,εx,q1r2+c12p222p22r2K2)
(B21)
Finally, it can be seen that the second condition Eq. (28) of Lemma 4.1 is satisfied if 
LV[ξ(t),t](q1p1r2+c12p1p222p22r2)V[ξ(t),t]+Kdδ
(B22)

In summary, the assumptions of Lemma 4.1 Eqs. (27) and (28) are satisfied by Eqs. (B3) and (B22), where (1/p2)=v1,(1/p1)=v2,((q1p1r2+c12p1p22)/2p22r2)=ϑ and Kdδ=μ. We conclude that the estimation error is exponentially bounded in mean square under the conditions in Eqs. (37)(39).

To meet these conditions, we use ε from Eq. (B21) and calculate δ as follows. We need to take care that ε̃ξ(t)ε with some ε̃ε satisfies the inequality in Eq. (B22) 
LV[ξ(t),t](q1p1r2+c12p1p222p22r2)V[ξ(t),t]+Kdδ0
(B23)
So, to guarantee the boundedness of the estimation error, we choose δ on the basis of Eq. (B23) and the lower bound of Eq. (B3) 
δ=(q1p1r2+c12p1p22)ε̃22p23r2Kd
(B24)

References

References
1.
Sub
,
F.
,
Varol
,
H. A.
, and
Goldfarb
,
M.
,
2010
, “
Upslope Walking With a Powered Knee and Ankle Prosthesis: Initial Results With an Amputee Subject
,”
IEEE Trans. Neural Syst. Rehabil. Eng.
,
19
(
1
), pp.
71
78
.
2.
Harvey
,
Z.
,
Potter
,
B. K.
,
Vandersea
,
J.
, and
Wolf
,
E.
,
2011
, “
Prosthetic Advances
,”
J. Surg. Orthop. Adv.
,
21
(
1
), pp.
58
64
http://www.armdynamics.com/caffeine/uploads/files/21-1-9.pdf.
3.
Martinez-Villalpando
,
E. C.
,
Weber
,
J.
,
Elliott
,
G.
, and
Herr
,
H.
,
2008
, “
Design of an Agonist-Antagonist Active Knee Prosthesis
,”
IEEE/RAS-EMBS
International Conference on Biomedical Robotics and Biomechatronics, Scottsdale, AZ, Oct. 19–22, pp.
529
534
.
4.
Au
,
S. K.
,
Herr
,
H.
,
Weber
,
J.
, and
Martinez-Villalpando
,
E. C.
,
2007
, “
Powered Ankle-Foot Prosthesis for the Improvement of Amputee Ambulation
,”
29th Annual IEEE Conference on Engineering in Medicine and Biology Society
(
EMBS
), Lyon, France, Aug. 22–26, pp.
3020
3026
.
5.
Au
,
S. K.
,
Weber
,
J.
, and
Herr
,
H.
,
2007
, “
Biomechanical Design of a Powered Ankle Foot Prosthesis
,”
IEEE International Conference on Rehabilitation Robotics
(ICORR)
, Noordwijk, The Netherlands, June 13–15, pp.
298
303
.
6.
Zhao
,
H.
,
Kolathaya
,
S.
, and
Ames
,
A. D.
,
2014
, “
Quadratic Programming and Impedance Control for Transfemoral Prosthesis
,”
International Conference on Robotics and Automation
(ICRA)
, Hong Kong, China, May 31–June 7, pp.
1341
1347
.
7.
Gregg
,
R. D.
,
Lenzi
,
T.
,
Hargrove
,
L. J.
, and
Sensinger
,
J. W.
,
2014
, “
Virtual Constraint Control of a Powered Prosthetic Leg: From Simulation to Experiments With Transfemoral Amputees
,”
IEEE Trans. Rob.
,
30
(
6
), pp.
1455
1471
.
8.
Young
,
A. J.
,
Simon
,
A. M.
, and
Hargrove
,
L. J.
,
2014
, “
A Training Method for Locomotion Mode Prediction Using Powered Lower Limb Prostheses
,”
IEEE Trans. Neural Syst. Rehabil. Eng.
,
22
(
3
), pp.
671
677
.
9.
Gonzalez
,
M.
,
2014
, “
Biomechanical Analysis of Gait Kinetics Resulting From Use of a Vacuum Socket on a Transtibial Prosthesis
,”
Master's thesis
, University of Nevada, Las Vegas, NVhttp://digitalscholarship.unlv.edu/cgi/viewcontent.cgi?article=1014amp;context=honors_theses.
10.
Murakami
,
T.
,
Nakamura
,
R.
,
Yu
,
F.
, and
Ohnishi
,
K.
,
1993
, “
Force Sensorless Impedance Control by Disturbance Observer
,”
IEEE Power Conversion Conference
, Yokohama, Japan, Apr. 19–21, pp.
352
357
.
11.
Simpson
,
J. W. L.
,
Cook
,
C. D.
, and
Li
,
Z.
,
2002
, “
Sensorless Force Estimation for Robots With Friction
,”
Australasian Conference on Robotics and Automation
, Auckland, New Zealand, Nov. 27–29, pp.
94
99
https://pdfs.semanticscholar.org/17ad/b33f0e6f6ceaa05dbecbe5781232709967b0.pdf.
12.
Phong
,
L. D.
,
Choi
,
J.
,
Lee
,
W.
, and
Kang
,
S.
,
2015
, “
A Novel Method for Estimating External Force: Simulation Study With a 4-DOF Robot Manipulator
,”
Int. J. Precis. Eng. Manuf.
,
16
(
4
), pp.
755
766
.
13.
Hacksel
,
P. J.
, and
Salcudean
,
S. E.
,
1994
, “
Estimation of Environment Forces and Rigid-Body Velocities Using Observers
,”
IEEE International Conference on Robotics and Automation
(ICRA)
, San Diego, CA, May 8–13, pp.
931
936
.
14.
Rigatos
,
G. G.
,
2012
, “
A Derivative-Free Kalman Filtering Approach to State Estimation-Based Control of a Class of Nonlinear Systems
,”
IEEE Trans. Ind. Electron.
,
59
(
10
), pp.
3987
3997
.
15.
Sarkka
,
S.
,
2007
, “
On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems
,”
IEEE Trans. Autom. Control
,
52
(
9
), pp.
1631
1641
.
16.
Julier
,
S. J.
,
Uhlmann
,
J. K.
, and
Durrant-Whyte
,
H. F.
,
1995
, “
A New Approach for Filtering Nonlinear Systems
,”
American Control Conference
(ACC)
, Seattle, WA, June 21–23, pp.
1628
1632
.
17.
Fakoorian
,
S. A.
,
Simon
,
D.
,
Richter
,
H.
, and
Azimi
,
V.
,
2016
, “
Ground Reaction Force Estimation in Prosthetic Legs With an Extended Kalman Filter
,”
IEEE International Systems Conference
(SysCon)
, Orlando, FL, Apr. 18–21, pp.
338
343
.
18.
Richter
,
H.
, and
Simon
,
D.
,
2014
, “
Robust Tracking Control of a Prosthesis Test Robot
,”
ASME J. Dyn. Syst. Meas. Control
,
136
(
3
), p.
031011
.
19.
van den Bogert
,
A. J.
,
Geijtenbeek
,
T.
,
Even-Zohar
,
O.
,
Steenbrink
,
F.
, and
Hardin
,
E. C.
, 2013, “
A Real-Time System for Biomechanical Analysis of Human Movement and Muscle Function
,”
Med. Biol. Eng. Comput.
,
51
(
10
), pp.
1069
1077
.
20.
Azimi
,
V.
,
Simon
,
D.
, and
Richter
,
H.
,
2015
, “
Stable Robust Adaptive Impedance Control of a Prosthetic Leg
,”
ASME
Paper No. DSCC2015-9794.
21.
Richter
,
H.
,
Simon
,
D.
,
Smith
,
W. A.
, and
Samorezov
,
S.
,
2015
, “
Dynamic Modeling, Parameter Estimation and Control of a Leg Prosthesis Test Robot
,”
Appl. Math. Modell.
,
39
(2), pp.
559
573
.
22.
Mohammadi
,
H.
, and
Richter
,
H.
,
2015
, “
Robust Tracking/Impedance Control: Application to Prosthetics
,”
American Control Conference
(ACC)
, Chicago, IL, July 1–3, pp.
2673
2678
.
23.
Warner
,
H.
,
2015
, “
Optimal Design and Control of a Lower-Limb Prosthesis With Energy Regeneration
,”
Master's thesis
, Cleveland State University, Cleveland, OHhttp://engagedscholarship.csuohio.edu/cgi/viewcontent.cgi?article=1678&context=etdarchive.
24.
Azimi
,
V.
,
Simon
,
D.
,
Richter
,
H.
, and
Fakoorian
,
S. A.
,
2016
, “
Robust Composite Adaptive Transfemoral Prosthesis Control With Non-Scalar Boundary Layer Trajectories
,”
American Control Conference
(ACC)
, Boston, MA, July 6–8, pp. 3002–3007.
25.
Warner
,
H.
,
Simon
,
D.
,
Mohammadi
,
H.
, and
Richter
,
H.
,
2016
, “
Switched Robust Tracking/Impedance Control for an Active Transfemoral Prosthesis
,”
American Control Conference
(ACC)
, Boston, MA, July 6–8, pp. 2187–2192.
26.
van der Merwe
,
R.
, and
Wan
,
E. A.
,
2001
, “
The Square-Root Unscented Kalman Filter for State and Parameter-Estimation
,”
IEEE International Conference on Acoustics, Speech, and Signal Processing
(ICASSP'01)
, Salt Lake City, UT, May 7–11, pp.
3461
3464
.
27.
Simon
,
D.
,
2006
,
Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches
,
Wiley
,
Hoboken, NJ
.
28.
Safonov
,
M.
, and
Athans
,
M.
,
1978
, “
Robustness and Computational Aspects of Nonlinear Stochastic Estimators and Regulators
,”
IEEE Trans. Autom. Control
,
23
(
4
), pp.
717
725
.
29.
Ljung
,
L.
,
1979
, “
Asymptotic Behavior of the Extended Kalman Filter as a Parameter Estimator for Linear Systems
,”
IEEE Trans. Autom. Control
,
24
(
1
), pp.
36
50
.
30.
Ursin
,
B.
,
1980
, “
Asymptotic Convergence Properties of the Extended Kalman Filter Using Filtered State Estimates
,”
IEEE Trans. Autom. Control
,
25
(
6
), pp.
1207
1211
.
31.
Reif
,
K.
,
Gunther
,
S.
,
Yaz
,
E.
, and
Unbehauen
,
R.
,
2000
, “
Stochastic Stability of the Continuous-Time Extended Kalman Filter
,”
IEEE Proc. Control Theory Appl.
,
147
(
1
), pp.
45
72
.
32.
Reif
,
K.
,
Gunther
,
S.
,
Yaz
,
E.
, and
Unbehauen
,
R.
,
1999
, “
Stochastic Stability of the Discrete-Time Extended Kalman Filter
,”
IEEE Trans. Autom. Control
,
44
(
4
), pp.
714
728
.
33.
Fitts
,
J. M.
,
1972
, “
On the Observability of Nonlinear Systems With Applications to Nonlinear Regression Analysis
,”
Inf. Sci.
,
4
(
2
), pp.
129
156
.
34.
Sontag
,
E. D.
,
1979
, “
On the Observability of Polynomial Systems—I: Finite-Time Problems
,”
SIAM J. Control Optim.
,
17
(
1
), pp.
139
151
.
35.
Gard
,
T. C.
,
1988
,
Introduction to Stochastic Differential Equations
,
M. Dekker
,
New York
.
36.
Zakai
,
M.
,
1967
, “
On the Ultimate Boundedness of Moments Associated With Solutions of Stochastic Differential Equations
,”
SIAM J. Control
,
5
(
4
), pp.
588
593
.
37.
Shanying
,
Z.
,
Chen
,
C.
,
Li
,
W.
,
Yang
,
B.
, and
Guan
,
X.
,
2013
, “
Distributed Optimal Consensus Filter for Target Tracking in Heterogeneous Sensor Networks
,”
IEEE Trans. Cybern.
,
43
(
6
), pp.
1963
1976
.
38.
Tønne
,
K. K.
,
2007
, “
Stability Analysis of EKF-Based Attitude Determination and Control
,”
Master's thesis
, Norwegian University of Science and Technology, Trondheim, Norwayhttp://folk.ntnu.no/tomgra/Diplomer/Tonne.pdf.
39.
Baras
,
J. S.
,
Bensoussan
,
A.
, and
James
,
M. R.
,
1988
, “
Dynamic Observers as Asymptotic Limits of Recursive Filters: Special Cases
,”
SIAM J. Appl. Math.
,
48
(
5
), pp.
1147
1158
.
40.
Gauthier
,
J. P.
,
Hammouri
,
H.
, and
Othman
,
S.
,
1992
, “
A Simple Observer for Nonlinear Systems Applications to Bioreactors
,”
IEEE Trans. Autom. Control
,
37
(
6
), pp.
875
880
.
41.
Khademi
,
G.
,
Mohammadi
,
H.
,
Hardin
,
E. C.
, and
Simon
,
D.
,
2015
, “
Evolutionary Optimization of User Intent Recognition for Transfemoral Amputees
,”
IEEE Biomedical Circuits and Systems Conference
(BioCAS)
, Atlanta, GA, Oct. 22–24, pp.
1
4
.
42.
Mao
,
X.
,
2007
,
Stochastic Differential Equations and Applications
,
Elsevier
,
Cambridge, UK
.