Systems and methods for prosthetic device control

Abstract

Systems and methods for prosthetic device control with a unified virtual constraint that controls an entire gait cycle of the prosthetic device.

Claims

1. A method for controlling a prosthetic device, the method comprising: measuring a plurality of variables with one or more sensors, wherein: the one or more sensors are operatively coupled to a controller; and the controller is configured to control the prosthetic device comprising a plurality of joint elements; determining a value of a monotonic phasing variable based on measurements of the plurality of variables; and adjusting a joint element of the plurality of joint elements in response to the value of the monotonic phasing variable by enforcing a unified virtual constraint that controls an entire gait cycle of the prosthetic device, wherein the monotonic phasing variable is monotonic throughout the entire gait cycle.

2. The method of claim 1 wherein the monotonic phasing variable is a hip position of a user of the prosthetic device.

3. The method of claim 1 wherein the monotonic phasing variable is a hip phase angle of a user of the prosthetic device.

4. The method of claim 1 wherein the monotonic phasing variable is a location of a center of mass of a user of the prosthetic device.

5. The method of claim 1 wherein the entire gait cycle is an ambulatory cycle.

6. The method of claim 5 wherein the ambulatory cycle is a walking cycle.

7. The method of claim 5 wherein the ambulatory cycle is a running cycle.

8. The method of claim 5 wherein the ambulatory cycle is a stair climbing cycle.

9. The method of claim 1 wherein the unified virtual constraint is created utilizing a Discrete Fourier Transform (DFT).

10. The method of claim 1 wherein the unified virtual constraint is created utilizing a polynomial regression fit.

11. The method of claim 1 wherein the joint element of the plurality of joint elements comprises a knee joint.

12. The method of claim 1 wherein the joint element of the plurality of joint elements comprises an ankle joint.

13. The method of claim 1 wherein the plurality of joint elements comprises a knee joint and an ankle joint.

14. The method of claim 1 wherein the plurality of variables comprises an angle of a knee joint.

15. The method of claim 1 wherein the plurality of variables comprises a velocity of a knee joint.

16. The method of claim 1 wherein the plurality of variables comprises an angle of an ankle joint.

17. The method of claim 1 wherein the plurality of variables comprises a velocity of an ankle joint.

18. The method of claim 1 wherein the entire gait cycle comprises a stance phase and a swing phase.

19. The method of claim 18 wherein different control parameters are used for the stance phase and the swing phase to enforce the unified virtual constraint for the entire gait cycle of the prosthetic device.

20. A control system comprising: a plurality of sensors configured to measure a plurality of variables relating to a prosthetic device; and a controller operatively coupled to the plurality of sensors, wherein the controller is configured to: determine a value of a monotonic phasing variable based on measurements of the plurality of variables; and adjust a joint element of the prosthetic device in response to the value of the monotonic phasing variable by enforcing a unified virtual constraint that controls an entire gait cycle of the prosthetic device, wherein the monotonic phasing variable is monotonic throughout the entire gait cycle.

21. The control system of claim 20 wherein the monotonic phasing variable is a hip position of a user of the prosthetic device.

22. The control system of claim 20 wherein the monotonic phasing variable is a hip phase angle of a user of the prosthetic device.

23. The control system of claim 20 wherein the monotonic phasing variable is a location of a center of mass of a user of the prosthetic device.

24. The control system of claim 20 wherein the entire gait cycle is an ambulatory cycle.

25. The control system of claim 24 wherein the ambulatory cycle is a walking cycle.

26. The control system of claim 24 wherein the ambulatory cycle is a running cycle.

27. The control system of claim 24 wherein the ambulatory cycle is a stair climbing cycle.

28. The control system of claim 20 wherein the unified virtual constraint is created utilizing a Discrete Fourier Transform (DFT).

29. The control system of claim 20 wherein the unified virtual constraint is created utilizing a polynomial regression fit.

30. The control system of claim 20 wherein the joint element comprises a knee joint.

31. The control system of claim 20 wherein the joint element comprises an ankle joint.

32. The control system of claim 20 wherein the plurality of joint elements comprises a knee joint and an ankle joint.

33. The control system of claim 20 wherein the plurality of variables comprises an angle of a knee joint.

34. The control system of claim 20 wherein the plurality of variables comprises an angle of an ankle joint.

35. The control system of claim 20 wherein the entire gait cycle comprises a stance phase and a swing phase.

36. The control system of claim 35 wherein the system further comprises different control parameters for the stance phase and the swing phase to enforce the unified virtual constraint for the entire gait cycle of the prosthetic device.

37. A prosthetic device comprising: a plurality of sensors configured to measure a plurality of variables relating to the prosthetic device; and a controller operatively coupled to the plurality of sensors, wherein the controller is configured to: determine a value of a monotonic phasing variable based on measurements of the plurality of variables; and adjust a joint element of the prosthetic device in response to the value of the monotonic phasing variable by enforcing a unified virtual constraint that controls an entire gait cycle of the prosthetic device, wherein the monotonic phasing variable is monotonic throughout the entire gait cycle.

38. The prosthetic device of claim 37 wherein the monotonic phasing variable is a hip position of a user of the prosthetic device.

39. The prosthetic device of claim 37 wherein the monotonic phasing variable is a hip phase angle of a user of the prosthetic device.

40. The prosthetic device of claim 37 wherein the monotonic phasing variable is a location of a center of mass of a user of the prosthetic device.

41. The prosthetic device of claim 37 wherein the entire gait cycle is an ambulatory cycle.

42. The prosthetic device of claim 41 wherein the ambulatory cycle is a walking cycle.

43. The prosthetic device of claim 41 wherein the ambulatory cycle is a running cycle.

44. The control system of claim 41 wherein the ambulatory cycle is a stair climbing cycle.

45. The prosthetic device of claim 37 wherein the unified virtual constraint is created utilizing a Discrete Fourier Transform (DFT).

46. The prosthetic device of claim 37 wherein the unified virtual constraint is created utilizing a polynomial regression fit.

47. The prosthetic device of claim 37 wherein the joint element comprises a knee joint.

48. The prosthetic device of claim 37 wherein the joint element comprises an ankle joint.

49. The prosthetic device of claim 37 wherein the plurality of joint elements comprises a knee joint and an ankle joint.

50. The prosthetic device of claim 37 wherein the plurality of variables comprises an angle of a knee joint.

51. The prosthetic device of claim 37 wherein the plurality of variables comprises an angle of an ankle joint.

52. The prosthetic device of claim 37 wherein the plurality of variables comprises a velocity of a knee joint.

53. The prosthetic device of claim 37 wherein the plurality of variables comprises a velocity of an ankle joint.

54. The prosthetic device of claim 37 wherein the entire gait cycle comprises a stance phase and a swing phase.

55. The prosthetic device of claim 54 wherein the system further comprises different control parameters for the stance phase and the swing phase to enforce the unified virtual constraint for the entire gait cycle of the prosthetic device.

Description

BRIEF DESCRIPTION OF THE FIGURES

(1) FIG. 1 illustrates a schematic view of gait phases used in the prior art to control lower prosthetic limb devices.

(2) FIG. 2 illustrates a schematic of a prosthetic device programmed with the control strategies according to exemplary embodiments.

(3) FIG. 3 illustrates a control strategy flow chart for the operation of the prosthetic device of FIG. 2.

(4) FIG. 4 illustrates a schematic of the unilateral, transfemoral amputee model.

(5) FIG. 5 illustrates a DFT plot for knee and ankle joint trajectories for the model of FIG. 4.

(6) FIG. 6 illustrates graphs of prosthetic leg output functions for ankle and knee joint angle trajectories versus the normalized phase variable for the model of FIG. 4.

(7) FIG. 7 illustrates graphs of angular phase portraits for prosthetic knee and ankle with K=5 and 10 for the output functions for the model of FIG. 4.

(8) FIG. 8 illustrates graphs of simulated joint angles and velocities of prosthetic knee and ankle for the K=10 output function, plotted against the normalized phase variable for the model of FIG. 4.

(9) FIG. 9 illustrates a schematic of a prosthetic device programmed with the control strategies according to exemplary embodiments.

(10) FIG. 10 illustrates a control system block diagram 550 for use with the embodiment of in FIG. 9.

(11) FIG. 11 illustrates plots of the tracking data from walking experiments utilizing exemplary embodiments at a constant speed.

(12) FIG. 12 illustrates a plot of the tracking data from walking experiments utilizing exemplary embodiments at various speeds.

(13) FIG. 13 illustrates empirical human gait data of the ankle and knee joints according to exemplary embodiments.

(14) FIG. 14 illustrates knee and ankle angular position step response according to exemplary embodiments.

(15) FIG. 15 illustrates a phase portrait of hip angle and hip velocity with the hip phase angle according to exemplary embodiments.

(16) FIG. 16 illustrates experimental leg bench top data of a unified controller for a knee joint according to exemplary embodiments.

DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS

(17) Although defined in a similar manner to physical/contact constraints, virtual constraints are enforced by actuator torques rather than external physical forces. The vast majority of virtual constraints used in bipedal robots are holonomic, [13], [23], [32], [12], [33]; therefore, one can translate these design principles into a prosthetic framework by considering virtual holonomic constraints:
h(q)=0(Eq. A)
where q is the vector of generalized coordinates (i.e., all joint angles) of the system, and his a vector-valued function defining one virtual constraint per actuated degree-of-freedom (DOF), e.g., the knee and ankle of the prosthesis. Eq. A is called a unified virtual constraint if it can be used to control the actuated joints through a complete gait cycle (i.e., one stride) without changes to the function h.

(18) Meaningful virtual constraints, i.e., output functions h(q), can be defined in a variety of ways. Reviewing previous work in bipedal robotics (e.g., [23], [13], [12]), virtual constraints can be used to control the actuated joints specified by h.sub.0(q)=(.sub.a, .sub.k).sup.T to a desired (vector-valued) trajectory h.sub.d((q)) as a function of a monotonic quantity (q). This quantity, known as the phase variable or timing variable, provides a unique representation of the gait cycle phase to drive forward the desired pattern in a time-invariant manner. In this case, the output function of (Eq. A) would be defined by h(q)=h.sub.o(q)h.sub.d((q)). In the prior art the desired pattern her is often defined first as a function of time (via boundary-constrained optimization) and then reparameterized into a function of (q) [23].

(19) Given desired virtual constraints (Eq. A), the goal of a virtual constraint controller is to drive output function h(q) to zero using actuator torques. Therefore, the control system output
y:=h(q)(Eq. B)
corresponds to tracking error from the desired constraint (Eq. A).

(20) Some torque control methods are better suited than others at driving this output to zero. If one does not start with a desired joint pattern ha (see discussion above), it is not always possible to solve a virtual constraint (Eq. A) for a unique joint trajectory as needed for joint impedance methods [4], [7], [9], [34]. Bipedal robots typically enforce virtual constraints using partial (i.e., input-output) feedback linearization [23], which has appealing theoretical properties including exponential convergence [11], reduced-order stability analysis [23], and robustness to model errors [13]. However, the application of this method to prosthetics presents unique challenges with human-machine interaction.

(21) FIG. 2 discloses a prosthetic limb 10 that may be programmed according to the exemplary embodiments of control methods described herein. In the embodiment shown, limb 10 comprises a powered ankle 21 and powered knee 31. Each joint can be coupled to an encoder (not shown) that can measure the angle and velocity of the joint. In addition, each joint is coupled to a motor 50 and optionally a transmission 55 that together are capable of producing physiological levels of torque. 50a indicates the motor coupled to the ankle and 50b indicates the motor coupled to the knee. In certain embodiments, e.g. for a user weighing up to 90 kg, at least 80 Nm of torque may be utilized. In particular embodiments, limb 10 and its related components may be the prosthetic limb with powered knee and ankle joints developed by the Vanderbilt Center for Intelligent Mechatronics (Vanderbilt University, TN) and commercialized by Freedom Innovations (Irvine, Calif.) or another powered prosthetic limb. Limb 10 and its components can be powered by a battery 15 in specific embodiments.

(22) One embodiment of the torque control method for a virtual constraint is impedance control (which is equivalent to proportional-derivative control in this context). In such an embodiment, the virtual constraint output would be changing Theta continuously. In the embodiment shown, limb 10 includes controller 100 that is programmed to run two control loops, a joint-level control loop 301 and a motor-level control loop 302 (shown in FIG. 3). In certain embodiments, controller 100 can have an update rate of at least 500 Hz or 1 kHz or more. If controller 100 accepts only impedance commands (i.e. stiffness Kp, viscosity Kd, and equilibrium angle Theta), its impedance controller can be replaced by a desired virtual constraint controller. This can be done for each joint by sending the controller 100 any standard Kp and Kd values and an equilibrium angle given by Theta=JointAngle+(Kd*Joint_Velocity+Torque_Desired)/Kp, where Torque_Desired is the torque calculated from sensor measurements by a torque controller, including for example, one of the torque controllers described herein. In exemplary embodiments, the joint-level loop 301 calculates and re-sends the equilibrium angle for each joint at the same or a different frequency at which the motor-level loop 302 operates.

(23) Controller 100 can be programmed with a control strategy for the control of prosthetic device 10. FIG. 3 discloses a flow chart for the operation of prosthetic device 10 during operation. In 300, joint-level control loop 301 receives measurements of one or more variables relating to elements of prosthetic device 10. In 310 and 320, joint-level control loop 301 calculates a value of a phasing variable and unified virtual constraint to control the entire gait cycle. Joint-level control loop 301 determines the unified virtual constraint error in step 330, and in step 340 calculates torque commands to eliminate the unified virtual constraint error. In step 350, joint-level control loop 301 transmits the torque commands to motor-level control loop 302. In step 360, motor-level control loop 302 then converts the torque commands to electrical current input values that indicate the amount of current to be applied to motors 50. In step 370, control loop 302 causes current to flow to motors 50 on the basis of the current input values. In certain embodiments utilizing closed-loop control, ammeters 45 can measure the current flowing to motors 50 and current flow may be corrected to the desired level if needed. In step 380 motors 50 and transmissions 55 output the appropriate torques to ankle joint 21 and knee joint 31. In particular embodiments, again utilizing closed-loop control, strain gages 35 may measure torques at ankle 21 and knee 31 and current flow may be corrected to achieve intended torque control. Controller 100 may then return to step 300.

(24) In one embodiment, controller 100 may comprise two microcontrollers, including one that implements motor-level control loop 302 and one microcontroller that implements joint-level control loop 301. The PIC32 (Microchip Technologies Inc., Chandler, Ariz.) or a similar microcontroller or a digital signal processor may be used for each microcontroller. The microcontrollers may communicate over the CAN protocol or another digital medium or an analog medium. The motor-level controller can accept torque commands to drive motors 50. Alternately, the motor-level controller may accept impedances instead of torque commands, in which case an additional impedance inversion stage is needed in the joint-level control loop 301. In certain embodiments, controller 100 may comprise a single microcontroller.

(25) As previously mentioned, the unified virtual constraint may also be created utilizing various techniques. In certain embodiments, the unified virtual constraint may be created utilizing a DFT. In other embodiments, the unified virtual constraint may be created utilizing other techniques, including for example, a polynomial regression fit.

(26) In exemplary embodiments the unified virtual constraint error may be calculated in one of various manners. In certain exemplary embodiments, the unified virtual constraint error may be calculated utilizing proportional derivative control. In other embodiments, the unified virtual constraint error may be calculated utilizing feedback linearization.

(27) In exemplary embodiments, the phasing variable may be, for example, the position of the center of mass or the hip position for a person utilizing prosthetic device 10. Any variable that can be measured or calculated and is monotonic throughout the entire gait cycle during a steady gait may be used as the phasing variable.

(28) Referring now to FIG. 4, a schematic of a unilateral, transfemoral amputee model 400 is shown comprising a prosthetic device 410 and a human leg 420. The generalized coordinates used in model 400 are indicated with q terms. In the model shown, angle q.sub.1 is unactuated and for modeling purposes, angles q.sub.2-6 have ideal actuators.

(29) Model 400 shown in FIG. 4 comprises seven leg segments plus a point mass at the hip to represent the upper body [18] and serve as the phasing variable. Calculation of the hip position is described in further detail herein, and the control systems described herein may use hip position as a phase variable for those tasks where hip position is monotonic during the entire gait cycle, including walking on a level surface, walking uphill, walking downhill, heel-to-toe running, or stair climbing.

(30) The full model shown in FIG. 4 is divided into a prosthesis subsystem comprising a prosthetic thigh, shank, and foot, and a human subsystem comprising a contralateral thigh, shank, and foot, a residual thigh on the amputated side, and a point mass at the hip. For modeling purposes, it is assumed that the prosthetic thigh and residual stump are rigidly attached. The thigh and shank segments are modeled using rigid links with both mass and inertia. Rather than model all of the degrees of freedom of the foot, the function of the foot and ankle is modeled using a circular foot [19], [20] plus an ankle joint to capture the stance ankle's positive work [21].

(31) For the purpose of simulating the full model, each stride (a complete gait cycle) can be modeled using four distinct periods. For simulation, a stride starts just after the transition from contralateral stance to prosthesis stance and proceeds through the prosthesis stance period, an impact event, the contralateral stance period and a second impact. The two stance periods can be modeled with continuous, second-order differential equations, and the two impact periods can be modeled using an algebraic mapping that relates the state of the biped at the instant before impact to the state of the biped at the instant after impact.

(32) In this simulation method, the equations of motion during the single support period for a subsystem can be found using the method of Lagrange [18], [22] and written as:
M.sub.ij{umlaut over (q)}.sub.i+C.sub.ij{dot over (q)}.sub.i+N.sub.ij=B.sub.iju.sub.i+J.sub.ij.sup.TF.(1)

(33) The subscript i indicates which subsystem the term is for, with P indicating the prosthesis subsystem and H indicating the human subsystem. The subscript j indicates which leg is in stance, with P indicating that the prosthesis is in stance and C indicating that the contralateral leg is in stance (and that the prosthesis is in swing). Because the biped is assumed to roll without slip, the absolute angle is unactuated. The ideal actuators at the remaining joints generate torque u.sub.i. The interaction forces between the human and prosthesis at the socket are given by F and depend on the motion of both subsystems. Eq. 1 can also be written as a first-order system:

(34) $\begin{array}{cc}{\stackrel{.}{x}}_i=f_{\mathrm{ij}}\left(x_i\right)+g_{\mathrm{ij}}\left(x_i\right)u_i+p_{\mathrm{ij}}\left(x_i\right)F,\mathrm{where}{x}_i=\left[\begin{array}{c}{q}_i\\ {\stackrel{.}{q}}_i\end{array}\right],f_{\mathrm{ij}}\left(x\right)=\left[\begin{array}{c}{\stackrel{.}{q}}_i\\ -M_{\mathrm{ij}}^{-1}\left(C_{\mathrm{ij}}{\stackrel{.}{q}}_i+N_{\mathrm{ij}}\right)\end{array}\right],g_{\mathrm{ij}}\left(x\right)=\left[\begin{array}{c}0\\ M_{\mathrm{ij}}^{-1}{B}_{\mathrm{ij}}\end{array}\right],p_{\mathrm{ij}}\left(x\right)=\left[\begin{array}{c}0\\ M_{\mathrm{ij}}^{-1}{J}_{\mathrm{ij}}^T\end{array}\right].& \left(2\right)\end{array}$

(35) The impacts can be modeled using equations of the form
q.sub.i.sup.+=q.sub.i.sup.,(3)
{dot over (q)}.sub.i.sup.+=A.sub.ij{dot over (q)}.sub.+A.sub.ij,(4)
where is the socket interaction impulse and depends on the pre-impact state of both subsystems [18], the superscript refers to the instant before impact, and the superscript + refers to the instant after impact.

(36) Feedback Linearizing Control

(37) To control both the prosthesis and the human, feedback linearization with virtual constraints is used [11], [23], although the prosthesis controller does not depend on the form of the human controller. Note the human leg is controlled in this manner for purposes of simulation. To perform feedback linearization, the desired motion of the actuated variables is encoded in output functions to be zeroed [23]. These output functions are the virtual constraints:
y.sub.ij=h.sub.ij(q.sub.i)=H.sub.0iq.sub.ih.sub.ij.sup.d(.sub.i(q.sub.i),(5)
where h.sub.ij is a vector-valued function to be zeroed, H.sub.0i is a matrix that maps the generalized coordinates to the actuated angles, h.sub.ij.sup.d is a vector-valued function of the desired joint angles, and .sub.i is the phase variable. For human leg 420, a separate output function is defined for the prosthesis single support period and for the contralateral single support period, and h.sub.Hj.sup.d is encoded using polynomials. For prosthetic device 410, a single, unified output function h.sub.P is defined for the entire stride, i.e., both the stance and swing periods of the prosthesis. Note that h.sub.ij and h.sub.ij.sup.d are functions of configuration only since the motion is parametrized using a kinematic phase variable and not by time explicitly. Differentiating Eq. 5 twice and substituting in the equations of motion (Eq. 2) for {umlaut over (q)}.sub.i gives the output dynamics [16]
.sub.ij=L.sub.f.sub.ij.sup.2h.sub.ij+L.sub.g.sub.ijL.sub.f.sub.ijh.sub.ij.Math.u.sub.i+L.sub.p.sub.ijL.sub.f.sub.ijh.sub.ij.Math.F,(6)
where standard Lie derivative notation [11] has been used and the terms are given by

(38) $L_{f_{\mathrm{ij}}}^2{h}_{\mathrm{ij}}=\frac{}{q_i}\left(\frac{h_{\mathrm{ij}}}{q_i}{\stackrel{.}{q}}_i\right){\stackrel{.}{q}}_i-\frac{h_{\mathrm{ij}}}{q_i}{M}_{\mathrm{ij}}^{-1}\left(C_{\mathrm{ij}}{\stackrel{.}{q}}_i+N_{\mathrm{ij}}\right),L_{g_{\mathrm{ij}}}{L}_{f_{\mathrm{ij}}}{h}_{\mathrm{ij}}=\frac{h_{\mathrm{ij}}}{q_i}{M}_{\mathrm{ij}}^{-1}{B}_{\mathrm{ij}},L_{p_{\mathrm{ij}}}{L}_{f_{\mathrm{ij}}}{h}_{\mathrm{ij}}=\frac{h_{\mathrm{ij}}}{q_i}{M}_{\mathrm{ij}}^{-1}{J}_{\mathrm{ij}}^T.$

(39) To cancel the nonlinearities in the output dynamics, set .sub.ij=v.sub.ij and solve for the input torques:
u.sub.i=.sub.ij+.sub.ij.Math.F,(7)
where v.sub.ij is a stabilizing controller, F is known through measurement, and
.sub.ij=L.sub.g.sub.ijL.sub.f.sub.ijh.sub.ij.sup.1(v.sub.ijL.sub.f.sub.ij.sup.2h.sub.ij),
.sub.ij=L.sub.g.sub.ijL.sub.f.sub.ijh.sub.ij.sup.1.Math.L.sub.p.sub.ijL.sub.f.sub.ijh.sub.ij.

(40) Virtual Constraint Design by DFT

(41) In this embodiment, the objective is to create output functions of virtual constraints that characterize the desired knee and ankle motion over the entire stride of prosthetic device 410. The output function is time-invariant and depends on a phase variable that is strictly monotonic and unactuated [23]. As previously stated, for this embodiment, the phase variable was chosen as the hip x-position q.sub.x, which is measured relative to a coordinate frame created at the impact transition from the contralateral stance to the prosthesis stance period. This phase variable monotonically increases for one complete gait stride, which allows one output function to represent one joint's motion over the entire stride. Taking advantage of the periodic kinematics observed in human gait [24], the method of DFT can be used to parameterize the desired output function.

(42) DFT is a linear transformation of a signal through a sequence of complex numbers to describe discrete frequency components of the signal. In general, the DFT representation is

(43) $\begin{array}{cc}X\left[k\right]=\underset{n=0}{\overset{N-1}{.Math.}}x\left[n\right]W_N^{\mathrm{kn}},k=0,1,\mathrm{.Math.},K& \left(8\right)\end{array}$
where N is the finite number of samples, K is the highest index for the finite sequence of k=0, 1, . . . N1, and W.sub.N=e.sup.j(2.sup..sup.N) is the complex quantity [25]. The original discrete signal in time domain is x[n], which is transformed to a discrete sequence in the frequency domain as X[k]. Because the signal is periodic, there are a finite number of discrete frequencies. As a result, Eq. 8 can be decomposed using a summation of sinusoids. To do so, Euler's relationship e.sup.j=cos j sin , can be used for W.sub.N within the DFT.

(44) Obtaining the frequency content terms X[k] from a signal, the original signal can be reconstructed using Fourier Interpolation, where the basis function is given by:

(45) $\begin{array}{cc}x\left[n\right]=\frac{1}{N}\underset{k=0}{\overset{N-1}{.Math.}}X\left[k\right]W_N^{-\mathrm{kn}},n{}_0& \left(9\right)\end{array}$
and from n=0, 1, . . . , N1.

(46) Note that both x[n] and X[k] may be complex. Defining X[k]=Re {X[k]+j Im{X[k]} and W.sub.N.sup.kn=Re{W.sub.N.sup.kn}+j Im{W.sub.N.sup.kn}, x[n] can be split into its real and imaginary parts:

(47) $\begin{array}{cc}\begin{array}{c}x\left[n\right]=\frac{1}{N}\underset{k=0}{\overset{N-1}{.Math.}}\left(\mathrm{Re}\left\{X\left[k\right]\right\}+j\mathrm{Im}\left\{X\left[k\right]\right\}\right)\\ \left(\mathrm{Re}\left\{W_N^{-\mathrm{kn}}\right\}+j\mathrm{Im}\left\{W_N^{-\mathrm{kn}}\right\}\right)\\ =\frac{1}{N}\underset{k=0}{\overset{N-1}{.Math.}}\left(\mathrm{Re}\left\{X\left[k\right]\right\}\mathrm{Re}\left\{W_N^{-\mathrm{kn}}\right\}-\mathrm{Im}\left\{X\left[k\right]\right\}\mathrm{Im}\left\{W_N^{-\mathrm{kn}}\right\}\right)+\\ j\left(\mathrm{Re}\left\{X\left[k\right]\right\}\mathrm{Im}\left\{W_N^{-\mathrm{kn}}\right\}+\mathrm{Im}\left\{X\left[k\right]\right\}\mathrm{Re}\left\{W_N^{-\mathrm{kn}}\right\}\right)\end{array}& \left(10\right)\end{array}$

(48) For the virtual constraint design, the only portion necessary for the Fourier Interpolation equation x[n] is the real part. It can be validated from the calculated DFT signal that X*[k]=X[k] and X*[Nk]=X[k] holds, which means that the imaginary components of Eq. 10 are zero [26]. The real, conjugate symmetric properties means the magnitude of X[k] for the k.sup.th index is the same for the (Nk).sup.th index [27]. Also, it is clear that the gait kinematic signal is a real-valued sequence. Therefore, if the signal meets the real, conjugate symmetric properties above then this simplifies Eq. 10 to:

(49) $\begin{array}{cc}x\left[n\right]=\frac{1}{N}\underset{k=0}{\overset{N-1}{.Math.}}\mathrm{Re}\left\{X\left[k\right]\right\}\mathrm{Re}\left\{W_N^{-\mathrm{kn}}\right\}-\mathrm{Im}\left\{X\left[k\right]\right\}\mathrm{Im}\left\{W_N^{-\mathrm{kn}}\right\}& \left(11\right)\end{array}$
where the summation terms are only the real part of the Fourier Interpolation. This can be represented by a function with a set of frequency terms from 0 up to N/2, the Nyquist sampling frequency, which is half the number of N samples in the discrete, equally spaced signal. Recall the magnitudes of X[k] and X[Nk] are equal due to the signal being a real-valued sequence, so the interpolation only needs to be computed for half the amount of N samples, as beyond N/2+1 the samples are duplicate values for the DFT. This results in the following exact representation of the original signal

(50) $\begin{array}{cc}x\left[n\right]=\frac{1}{2}{}_0+\underset{k=1}{\overset{\frac{N}{2}-1}{.Math.}}\left[{}_k\mathrm{Re}\left\{W_N^{-\mathrm{kn}}\right\}-{}_k\mathrm{Im}\left\{W_N^{-\mathrm{kn}}\right\}\right]+\frac{1}{2}{}_{\frac{N}{2}}\mathrm{Re}\left\{W_N^{-\frac{N}{2}n}\right\},n=0,1,\mathrm{.Math.},N-1& \left(12\right)\end{array}$
where .sub.k=2 Re{X[k]}R and .sub.k=2 Im{X[k]}R are the scalar coefficients [28]. Note that Eq. 12 is a lower order functional representation of Eq. 11.

(51) The Fourier interpolation presented in Eq. 12 is used to create the unified desired output function, h.sup.d.sub.P, for the entire gait cycle. Because of the unified prostheses controller there is a single output function, h.sub.P, used for both stance and swing period. This is in contrast to the human controller that requires four output functions h.sub.HP and h.sub.HC. To do so, the desired angular trajectories of the knee and ankle are sampled over the entire gait cycle (e.g., from [21], [24], and [29]) at n discrete, equally-spaced values of the phase variable. The phase variable can be normalized using

(52) $\begin{array}{cc}{S}_P\left({}_P\left(q_P\right)\right)=\frac{{}_P-{}_P^{+}}{{}_P^{-}-{}_P^{+}}& \left(13\right)\end{array}$
where the + signifies the start of the stance period for the prosthetic device 410 and the indicates the end of its swing period. Using Eq. 8 to solve for X[k], the real and imaginary parts of X[k] are used to compute the coefficients of Eq. 12. Eq. 12 is used in the output function (Eq. 5) to generate the prosthesis virtual constraints:

(53) $\begin{array}{cc}{h}_P\left(q_P\right)=H_{0P}.Math.q_P-\left(\frac{1}{2}{}_0+\frac{1}{2}{}_{\frac{N}{2}}\cos \left({\mathrm{Ns}}_P\right)+\underset{k=1}{\overset{\frac{N}{2}-1}{.Math.}}\left[{}_k\cos \left(2{\mathrm{ks}}_P\right)-{}_k\sin \left(2{\mathrm{ks}}_P\right)\right]\right).& \left(14\right)\end{array}$

(54) Simulation

(55) Using the method from the above section entitled Virtual Constraint Design by DFT, a unified output function for a prosthetic device 410 can be designed, comprising a knee and ankle joint as shown in FIG. 4. To begin, a desired trajectory that generates a stable gait is needed. For a prosthesis controller, it is generally desirable to mimic human able-bodied motion, so the desired trajectory can be chosen from one of many gait studies [24], [29]. For this simulation work, the trajectory was chosen based on a stable, HZD-based model designed to predict healthy human walking [18]. Using MATLAB to determine the DFT X[k] of the desired gait trajectory, the real and imaginary parts of X[k] and W.sub.N.sup.kn can be computed to create the Fourier Interpolation in Eq. 12. The DFT plot shown in FIG. 5 for the predetermined knee and ankle joint trajectories shows that between the 10.sup.th and the Nyquist sampling frequency of N/2, the magnitude is approximately zero and will not have significant impact in recreating the joint trajectory for h.sup.d.sub.P. As a result, only the first 10 terms in Eq. 12 provide useful information for the Fourier Interpolation, so the series can be truncated to k{1, . . . ,10}. Therefore, the linear transformation for h.sub.i.sup.d will have a minimal number of coefficients, which lessens the computational complexity of the feedback linearization around y.sub.q.

(56) For each joint, desired output functions with K=5, 10 and N/2 were found, where K is the highest index for k in Eq. 14 (see FIG. 6). The lowest order function K=5 is an adequate representation of the original joint trajectory for both the ankle and knee with a coefficient of determination, r=0.995 and r=0.999, respectively. As expected, the output functions for K=10 and N/2 are almost identical, and match the desired trajectory better than K=5. The coefficient of determination for both functions at both joints is r=1.0.

(57) The two lower-order output functions (K=5 and 10) were integrated into the amputee biped simulation as the virtual constraints for prosthetic device 410. For each set of output functions, the human leg 420 plus prosthetic device 410 was simulated until steady-state was reached (see FIG. 7). The two output functions converged to different periodic orbits due to the differences in the desired output functions. The steady-state prosthesis joint angles and velocities for K=10 are shown in FIG. 8. Given the faster convergence to a periodic orbit with K=10 (see FIG. 7) than the output function with K=5, the K=10 desired output function is preferred in this embodiment.

Experiments

(58) FIG. 9 also discloses a prosthetic limb 500 according to exemplary embodiments. This embodiment comprises a joint encoder 505, a motor amplifier 510, a first belt drive 515, a brushless direct current (BLDC) motor 520, a ball screw 525, and a ball nut 530. In addition, this embodiment comprises a motor filter card 535, a journal bearing 540, a second belt drive 545, and a load cell 550. Prosthetic limb 500 further comprises an upper support member 560 and a lower support member 570.

(59) In the embodiment shown in FIG. 3, a linear ball screw actuator with a lever arm at the joint is utilized to satisfy the torque requirements, while limiting the overall weight of prosthetic limb 500. In the illustrated embodiment, the ball screw actuator is a duplicate design between both the knee and ankle actuation system for ease in overall design and manufacturing. In the embodiment shown, BLDC motor 520 is a high power-to-weight ratio BLDC motor, and in specific embodiments is a Maxon 200 Watt three phase BLDC Motor (Model: EC-4pole 30, Maxon Motor, Sachseln, Switzerland) weighing 300 grams. Such a configuration can provide the power to weight ratio needed to meet the actuator torque output at each joint.

(60) For the knee actuator (with a 40 Nm torque requirement), a 2:1 timing belt drive can be used between the motor output and the ball screw. In the embodiment shown, the ankle actuator comprises a 4:1 timing belt drive. Sprockets utilized in prosthetic limb 500 can be machined from 7075 aluminum to reduce weight. In particular embodiments, the ball screw for both the knee and ankle can be a Nook 12 mm diameter, 2 mm lead ball screw (Model: PMBS 12x2r-3vw/0/T10/00/6K/184/0/S, Nook Industries, Cleveland, Ohio, USA). Such a configuration can convert the rotary motion from the motor with belt drive to linear motion of the ball nut driving the lever arm attached to the joint in generating torque output. The ball screw can be supported axially and radially by a double bearing support journal bearing (for example, model EZBK10-SLB, Nook Industries, Cleveland, Ohio, USA). To obtain the full range of motion at each joint with a linear actuator design, a motor mount can be used with a rotational pivot point to eliminate buckling of the ball screw and increase linear motion the ball screw can travel attached to the lever arm.

(61) FIG. 10 illustrates a control system block diagram 550 of the prosthetic limb 500 show in FIG. 9. In this embodiment, a dSPACE microcontroller contains the algorithms in computing the normalized phase variable sp from the IMU measuring hip angle as input to the virtual constraints. The virtual constraints produce knee error .sub.k and ankle error .sub.a to the control law in computing desired current commands for both the knee i.sub.k and ankle i.sub.. The desired current commands are sent to the motor drivers for handling the low-level current control loop to commutate the motors to drive the timing belt and ball screw transmission to actuate the joints.

(62) FIG. 11 and FIG. 12 illustrate plots of the tracking data from walking experiments utilizing the embodiment of FIG. 9. In the illustrated plots, the cmd signal is the desired joint angle produced by the control strategy as a function of the phase variable measurement, which the motors attempt to track. The meas signal is the actual measured joint angle, demonstrating the ability of the prosthesis to track the commanded signal. The plots in FIG. 11 and FIG. 12 show the averaged signals over multiple walking steps on a treadmill. The shaded region represents +/one standard deviation about the average. FIG. 11 illustrates results for walking at approximately 1 mile/hour. FIG. 12 illustrates walking at various speeds from 1 to 3 mile/hour, demonstrating the ability of the controller to automatically adjust walking speed based on measurements of the phase variable.

(63) FIG. 13 illustrates empirical human gait data of the ankle (left column) and knee (right column) joints for able-bodied subjects at normal walking speeds. The left column plots indicate the empirical joint kinematics for the ankle, and the right column plots are for the knee. The top row illustrates the ankle and knee joint angles in radians, the middle row displays the joint torques in Nm per weight in kg, and the bottom row shows the joint mechanical power in Watts per weight in kg. Note, the vertical dashed line indicates the start of the swing period.

(64) To evaluate the performance of the embodiment shown in FIG. 9, a linear proportional-derivative (PD) controller was implemented in position mode to control the leg joints individually. In other embodiments, a different controller may be used in the joint position control loop. The PD controller consist of a control law u.sub.PD=K.sub.P (.sub.d)K.sub.d {dot over ()}, where .sub.d is the desired joint angular position, is the measured joint angular position, {dot over ()} is the measured joint angular velocity, K.sub.p is the proportional gain affecting the stiffness of the joint, and K.sub.d is the derivative gain behaving as a damping term. Step inputs were evaluated at each joint to obtain dynamic response characteristics to evaluate the performance of the system. A position step input of 10 degrees (0.1745 radians) was commanded separately for both knee and ankle (see left column of FIG. 14). Response characteristics were computed as follows (all times in seconds): rise time 0.917 (knee) 0.205 (ankle); peak time 2.57 (knee) 1.64 (ankle); settling time 1.48 (ankle). The ankle step response has better performance with faster response and settling time than the knee. Note the knee's actuation system receives a larger inertia to overcome with the opposing moment due to the mass of the leg during free swing. Thus, the knee has a slower response time compared to the ankle, where the settling time could not be computed based on the undershoot response. Control gains were selected for stable responses. Gains were increased to improve the control response of the knee, but under disturbances it proved to be unstable.

(65) A step duty cycle was commanded to both joints individually to evaluate transient responses during small time durations (see right column of FIG. 14). The desired position commands applied in order were [0, 0.01745, 0, 0.01745, 0, 0.08727, 0.1745, 0.2618, 0.08727, 0.1745, 0.2618, 0.1745, 0.08727, 0, 0.08727, 0 radians] at various time intervals over time length of 30 seconds. Again the ankle actuator was able to meet the desired commanded position, while the knee had steady-state error due to its undershoot behavior response. However, experimentation with the leg attached to amputees, an aided torque will be applied to the knee joint due to the moment created from the thigh motion during swing period. While these results are preliminary bench top testing, future work may entail performing system identification techniques to accurately determine the bandwidth of each actuator. It is important to note it may potentially be a destructive test, particularly for the knee, if a full spectrum analysis is to be performed.

(66) To implement advance control methods for a powered prosthetic leg, it's important to have a fast, real-time control processor to compute the algorithms efficiently and effectively for experimentation purposes with minimal processing delay. Other robotic systems have had success using real-time control for their algorithms, especially on bipedal robots like ATRIAS [32] and MABEL [33]. The real-time control interface chosen for the prosthetic leg is a dSPACE DS1007, where both real-time control and data acquisition are performed. All electronic signals from the leg are sent to the dSPACE processor board for computation and are monitored real-time using a PC with dSPACE ControlDesk Next Generation software. The control algorithms are programmed using MATLAB/Simulink (Version: 2013b, The Mathworks, Natick, Mass., USA). Furthermore, safety conditions were placed with software logic and emergency stop button to safely operate the leg in case a failure was to occur.

(67) The main areas in implementing the control algorithms for the leg include, but not limited to, the virtual constraints, the phase variable, the sensors signal processing/filtering, and the kinematic analysis for torque calculations in performing closed-loop torque control. The virtual constraints were designed offline from empirical able-bodied walking gaits ([24] and FIG. 13) and programmed in MATLAB/Simulink to interface with dSPACE. The overall controller implemented to operate the leg to produce preliminary results was a PD control law using the previously-defined virtual constraints:
u.sub.PD=K.sub.p(H.sub.0PqPh.sub.P.sup.d(s.sub.P))K.sub.d(H.sub.0PqP),(15)

(68) where as before Kp and Kd are gains to control stiffness and damping in the system, respectively. The proposed phase variable sP is based on previous work of finding a robust parameterization of human locomotion from the hip phase angle [35]. Results from human subjects have shown that this approach produces a monotonic and bounded phase variable, even during perturbation. Experiments were performed using motion capture cameras (Vicon, Oxford, UK) in collecting kinematic data of human walking over a walkway with a perturbation mechanism. The phase variable was computed offline by post processing the kinematic data. A major challenge was recreating results from post processed data to real-time computation of the phase variable.

(69) The first step for the control algorithms was producing the phase variable algorithms for real-time measurements. Serial communication for the IMU sensor were developed to output Euler Angles at a sample rate of 500 Hz, and the sensor is mounted along the human's thigh in the sagittal plane measuring the human's hip angle. Eq. 16 presents the normalized phase variable function based on the hip phase angle [35]:

(70) 0$\begin{array}{cc}{s}_P=\frac{\mathrm{atan}2\left(z\left(\stackrel{~}{}+\stackrel{~}{}\right).Math.+\right)}{2},\mathrm{where}z=\frac{.Math.{}_{\max }-{}_{\min }.Math.}{.Math.{\stackrel{~}{}}_{\max }-{\stackrel{~}{}}_{\min }.Math.},=-\left(\frac{{}_{\max }+{}_{\min }}{2}\right),\stackrel{~}{}=-\left(\frac{{\stackrel{~}{}}_{\max }+{\stackrel{~}{}}_{\min }}{2}\right)=+,\stackrel{~}{}={}_0^t\left(\right)d& \left(16\right)\end{array}$

(71) The variable is the hip angle with respect to the vertical/gravity. The variable is chosen such that .sub.0 ()d=.sub.0.sup.T(()+) d=0, where T represents the total amount of time of a gait cycle. The variable z is referred to as the scale factor, and the variables and are what are known as the shift parameters for the hip position and integral, respectively. They are computed from the bounds of the derived phase portrait of the hip angle and integral. The phase variable is normalized from 0 to 2 (based on the orbit of the phase portrait) for a final range of zero to one. The maximum/minimum variables are computed based a four-quadrant window from the phase portrait. Equation 16 is the actual phase variable function as referenced in Eq. 13. FIG. 15 graphically displays the computed maximum/minimum values as the phase angle enters each quadrant (i.e. 0%, 25%, 50%, etc), where an updated maximum or minimum is computed for every passing quadrant of the phase portrait. In FIG. 15, the negative x-axis of quadrant 3 is considered 0% of the phase variable as it increases to a full counter-clockwise revolution to 100%. Each axis crossing of thy phase angle will contain either a maximum or minimum value of the hip angle or hip integral to determining the phase. In other embodiments, the phase may be calculated from the hip angle and hip velocity.

(72) Several preliminary bench top testing were performed to verify the functionality of the unified controller operating on a powered prosthesis. Before operating the leg with the real-time computation of the phase variable, an ideal phase variable signal was constructed to observe how the virtual constraint performs from ideal conditions. The ideal phase variable was a monotonic, increasing linear signal from 0 to 1 to replicate a human walking from 0 to 100% gait cycle. This aided in determining initial control gains to program into the leg for experimentation with human subjects. Bench top testing entailed using a continuous, ideal phase variable for real-time control tuning.

(73) Once response showed low tracking error and adequate transient response, then the controller was switched to use actual real-time measurements of the phase variable. A start/stop bench top experiment was developed where the able-bodied subject began at a rest position, transition to continuous walking, and ended at rest. FIG. 16 illustrates experimental leg bench top results of unified controller for the knee joint (top left) and the ankle joint (top right) using actual phase variable measurements (bottom) from an able-bodied subject at slow walking speed 1.0 miles per hour on a treadmill. In the top two plots. the dotted line is the controller reference command and the solid line is the joint angle response output. In the bottom plot, the dashed line is the phase variable. The subject initiated walking after 6 seconds with stopping after the 18 second mark as illustrated in FIG. 16.

(74) Notice the bottom plot in FIG. 16 that the phase variable is not exactly linear, however it is still monotonic, which is one of the required properties for using a phase variable for a virtual constraint controller. Recall the data in FIGS. 11 and 12 was collected using an exemplary embodiment of the prosethetic device disclosed in FIG. 9. In contrast, the data in FIG. 16 was collected with an ablebodied human walking on a treadmill (Model: Commercial 1750, NordicTrack, Chaska, Minn., USA) at 1.0 miles per hour (slow walking speed) wearing the IMU sensor while the leg is mounted fixed to a table. The human subject started from rest position, then turned on the treadmill to begin walking continuously, and eventually turned off the treadmill to end at rest. Overall, the leg tracked the response well for bench top testing. The leg inherently begins and ends at a stance position, where the knee has small flexion at safe phase variable mode of 0.5. This allows the leg to not collapse with large knee flexion, so the human subject will not fall down.

CONCLUSIONS

(75) The development of a single, unified controller that captures the entire gait cycle eliminates the need to divide the gait into different periods, each with their own controllers. DFT methods were used to produce a unified virtual constraint for a feedback-linearizing controller based on a desired trajectory. While the feasibility of the controller was demonstrated using a walking gait in this embodiment, the DFT method approach can be applied equally well to other locomotion activities with well-characterized joint kinematics. For instance, the control strategy could be applied to different tasks, including for example, stair and slope ascent/descent.

(76) The systems and methods described herein can be used to control other powered lower limb devices besides limb 10. The systems and methods described herein can be used to control a single joint powered lower limb device, such as an ankle prosthesis, or a multi joint powered lower limb device, such as limb 10, which includes a motorized ankle and a motorized knee.

(77) All of the devices, systems and/or methods disclosed and claimed herein can be made and executed without undue experimentation in light of the present disclosure. While the devices, systems and methods of this invention have been described in terms of particular embodiments, it will be apparent to those of skill in the art that variations may be applied to the devices, systems and/or methods in the steps or in the sequence of steps of the method described herein without departing from the concept, spirit and scope of the invention. All such similar substitutes and modifications apparent to those skilled in the art are deemed to be within the spirit, scope and concept of the invention as defined by the appended claims.

REFERENCES

(78) The contents of the following references are incorporated by reference herein: [1] K. Ziegler-Graham, E. J. MacKenzie, P. L. Ephraim, T. G. Travison, and R. Brookmeyer, Estimating the prevalence of limb loss in the united states: 2005 to 2050, Archives of physical medicine and rehabilitation, vol. 89, no. 3, pp. 422-429, 2008. [2] W. C. Miller, A. B. Deathe, M. Speechley, and J. Koval, The influence of falling, fear of falling, and balance confidence on prosthetic mobility and social activity among individuals with a lower extremity amputation, Archives of physical medicine and rehabilitation, vol. 82, no. 9, pp. 1238-1244, 2001. [3] R. S. Galley, M. A. Wenger, M. Raya, N. Kirk, K. Erbs, P. Spyropoulos, and M. S. Nash, Energy expenditure of trans-tibial amputees during ambulation at self-selected pace, Prosthetics & Orthotics Int., vol. 18, no. 2, p. 84, 1994. [4] F. Sup, A. Bohara, and M. Goldfarb, Design and control of a powered transfemoral prosthesis, The International journal of robotics research, vol. 27, no. 2, pp. 263-273, 2008. [5] R. D. Bellman, M. A. Holgate, and T. G. Sugar, Sparky 3: Design of an active robotic ankle prosthesis with two actuated degrees of freedom using regenerative kinetics, in Biomedical Robotics and Biomechatronics, 2008. BioRob 2008. 2nd IEEE RAS & EMBS International Conference on. IEEE, 2008, pp. 511-516. [6] S. K. Au and H. Herr, Powered ankle-foot prosthesis, Robotics & Automation Magazine, IEEE, vol. 15, no. 3, pp. 52-59, 2008. [7] F. Sup, H. A. Varol, and M. Goldfarb, Upslope walking with a powered knee and ankle prosthesis: initial results with an amputee subject, Neural Systems and Rehabilitation Engineering, IEEE Transactions on, vol. 19, no. 1, pp. 71-78, 2011. [8] N. P. Fey, A. Simon, A. J. Young, and L. J. Hargrove, Controlling knee swing initiation and ankle plantarflexion with an active prosthesis on level and inclined surfaces at variable walking speeds, Translational Engineering in Health and Medicine, IEEE Journal of, vol. 2, pp. 1-12, 2014. [9] M. F. Eilenberg, H. Geyer, and H. Herr, Control of a powered ankle-foot prosthesis based on a neuromuscular model, Neural Systems and Rehabilitation Engineering, IEEE Transactions on, vol. 18, no. 2, pp. 164-173, 2010. [10] A. M. Simon, K. A. Ingraham, N. P. Fey, S. B. Finucane, R. D. Lipschutz, A. J. Young, and L. J. Hargrove, Configuring a powered knee and ankle prosthesis for transfemoral amputees within five specific ambulation modes, PloS one, vol. 9, no. 6, p. e99387, 2014. [11] A. Isidori, Nonlinear Control Systems, 3rd ed. London, England: Springer, 1995. [12] A. E. Martin, D. C. Post, and J. P. Schmiedeler, Design and experimental implementation of a hybrid zero dynamics-based controller for planar bipeds with curved feet, The International Journal of Robotics Research, vol. 33, no. 7, pp. 988-1005, 2014. [13] K. Sreenath, H.-W. Park, I. Poulakakis, and J. W. Grizzle, A compliant hybrid zero dynamics controller for stable, efficient and fast bipedal walking on mabel, The International Journal of Robotics Research, vol. 30, no. 9, pp. 1170-1193, 2011. [14] R. D. Gregg, E. J. Rouse, L. J. Hargrove, and J. W. Sensinger, Evidence for a time-invariant phase variable in human ankle control, PLoS ONE, vol. 9, no. 2, p. e89163, 2014, doi: 10.1371/journal.pone.0089163. [15] R. D. Gregg and J. W. Sensinger, Biomimetic virtual constraint control of a transfemoral powered prosthetic leg, in American Control Conference (ACC), 2013. IEEE, 2013, pp. 5702-5708. [16] R. D. Gregg, T. Lenzi, L. J. Hargrove, and J. W. Sensinger, Virtual constraint control of a powered prosthetic leg: From simulation to experiments with transfemoral amputees, IEEE Trans. Robotics, 2014, in press. [17] R. D. Gregg, T. Lenzi, N. P. Fey, L. J. Hargrove, and J. W. Sensinger, Experimental effective shape control of a powered transfemoral prosthesis, in IEEE Int. Conf. Rehab. Robotics, Seattle, Wash., 2013. [18] A. E. Martin and R. D. Gregg, Hybrid Invariance and Stability of a Feedback Linearizing Controller for Powered Prostheses, in American Control Conference, 2015, submitted. [19] A. H. Hansen, D. S. Childress, and E. H. Knox, Roll-over shapes of human locomotor systems: Effects of walking speed, Clinical Biomechanics, vol. 19, no. 4, pp. 407-14, 2004. [20], Prosthetic foot roll-over shapes with implications for alignment of trans-tibial prostheses, Prosthetics and Orthotics International, vol. 24, no. 3, pp. 205-15, 2000. [21] A. E. Martin and J. P. Schmiedeler, Predicting human walking gaits with a simple planar model, Journal of Biomechanics, vol. 47, no. 6, pp. 1416-21, 2014. [22] R. M. Murray, Z. Li, and S. S. Sastry, A Mathematical Introduction to Robotic Manipulation, 1st ed. Boca Raton, Fla.: CRC Press, 1994. [23] E. R. Westervelt, J. W. Grizzle, C. Chevallereau, J. H. Choi, and B. Morris, Feedback Control of Dynamic Bipedal Robot Locomotion. Boca Raton, Fla.: CRC Press, 2007. [24] D. Winter, Biomechancis and Motor Control of Human Movement. Hoboken, N.J.: John Wiley and Sons. Inc, 2009. [25] A. V. Oppenheim, R. W. Schafer, J. R. Buck, et al., Discrete-Time Signal Processing. Upper Saddle River, N.J.: Prentice-Hall Englewood Cliffs, 1989, vol. 2. [26] A. D. Poularikas, Transforms and Applications Handbook. Boca Raton, Fla.: CRC press, 2009. [27] R. G. Lyons, Understanding Digital Signal Processing. Pearson Education, 2010. [28] W. L. Briggs et al., The DFT: An Owners' Manual for the Discrete Fourier Transform. Philadelphia, Pa.: Siam, 1995. [29] J. Perry and J. Burnfield, Gait Analysis: Normal and Pathological Function. Thorofare, N.J.: Slack-Incorporated, 2010. [30] T. Bennett, R. Jafari, and N. Gans, An extended kalman filter to estimate human gait parameters and walking distance, in American Control Conference (ACC), 2013. IEEE, 2013, pp. 752-757. [31] D. J. Villarreal and R. D. Gregg, A survey of phase variable candidates of human locomotion, in IEEE Engineering in Medicine and Biology Conference, Chicago, Ill., August 2014. [32] S. Kolathaya and A. D. Ames, Achieving bipedal locomotion on rough terrain through human-inspired control, presented at the IEEE Int. Symp. Safety Security Rescue Robot. College Station, Tex., USA, 2012 [33] A. Ramezani, J. W. Hurst, K. A. Hamed, and J. W. Grizzle, Performance analysis and feedback control of ATRIAS, a 3D bipedal robot, ASME J. Dyn. Syst. Meas. Control, vol. 136, no. 2, p. 021012, 2013. [34] A. Simon, N. Fey, S. Finucane, R. Lipschutz, and L. Hargrove, Strategies to reduce the configuration time for a powered knee and ankle prosthesis across multiple ambulation modes, presented at the IEEE Int. Conf. Rehab. Robotics, Seattle, Wash., USA, 2013. [35] D. J. Villarreal, H. Poonawala, and R. D. Gregg. A Robust Parameterization of Human Joint Patterns Across Phase-Shifting Perturbations. IEEE Trans Neural Sys and Rehab Eng, 2015, under review.