Data-Driven Variable Impedance Control of a Powered Knee–Ankle Prosthesis for Adaptive Speed and Incline Walking

Most impedance-based walking controllers for powered knee–ankle prostheses use a finite state machine with dozens of user-specific parameters that require manual tuning by technical experts. These parameters are only appropriate near the task (e.g., walking speed and incline) at which they were tuned, necessitating many different parameter sets for variable-task walking. In contrast, this article presents a data-driven, phase-based controller for variable-task walking that uses continuously variable impedance control during stance and kinematic control during swing to enable biomimetic locomotion. After generating a data-driven model of variable joint impedance with convex optimization, we implement a novel task-invariant phase variable and real-time estimates of speed and incline to enable autonomous task adaptation. Experiments with above-knee amputee participants ($N=2$) show that our data-driven controller 1) features highly linear phase estimates and accurate task estimates, 2) produces biomimetic kinematic and kinetic trends as task varies, leading to low errors relative to able-bodied references, and 3) produces biomimetic joint work and cadence trends as task varies. We show that the presented controller meets and often exceeds the performance of a benchmark finite state machine controller for our two participants, without requiring manual impedance tuning.

Impedance control is a common strategy in lower limb wearable robotics because of its simplicity and ability to produce behaviors that are similar to human biology, such as a compliantly controlled interaction with the ground [10] and dynamics similar to what has been observed in skeletal muscles [11]. Further, empirical studies have shown that ankle joint dynamics during walking are well described with an impedance controller [12], [13], [14]. A standard impedance controller calculates joint torque τ based on a joint angle θ and joint velocityθ as where K, B, and θ eq are parameters defining the joint's stiffness, damping, and equilibrium angle, respectively. Traditional methods of impedance control for lower limb prostheses involve segmenting the gait cycle into discrete subphases, where each subphase has its own constant values of K, B, and θ eq . Researchers manually tune the impedance parameters in each subphase until the observed gait is satisfactory [6], [15], [16], [17], [18], [19]. Switching between subphases is controlled by a finite state machine (FSM) with transition criteria based on sensor readings (e.g., elapsed time, leg loading, joint angles, etc.). Like the impedance parameters, these transition criteria are often experimentally tuned for an individual's gait by a technical expert. More elaborate impedance value representations have been suggested [19], [20], [21], [22], but these methods still required manual, expert tuning.
Joint kinematics and kinetics vary based on the ground incline and walking speed [23], [24] (together termed the user's task). Therefore, the necessary impedance parameters and state machine transition criteria also vary. For a standard FSM impedance controller to operate over a wide array of tasks, many tunable parameters are required. For example, one multimodal impedance controller required a total of 140 tunable parameters for five ambulation modes [20]. While only a portion of these parameters were considered necessary to tune, the device's configuration and tuning still required the researchers up to five hours to complete.
In contrast to the standard FSM-based impedance control paradigm, some authors have suggested using continuous functions to define the impedance parameters and how they evolve over the gait cycle [25], [26], [27], [28]. In general, controllers that continually vary a robot's output mechanical impedance with time are known as variable impedance controllers [29]. Biomechanical principles suggest that human joints behave like variable impedance controllers [30] and empirical studies have observed this behavior at the ankle joint during walking [12], [13], [14]. Therefore, variable impedance control may offer a biomimetic solution for controlling powered prosthetic legs. However, how to appropriately define the variable impedance functions to realize walking gaits remains an open question.
A variable impedance controller was suggested in [26] using linear functions for stiffness and damping during stance. The linear functions were hand-tuned and held constant regardless of task. The variable impedance control method in [25] eliminated tuning altogether by using able-bodied kinematic data to generate continuous impedance parameter functions of gait phase. However, this method was limited to the knee joint, did not consider joint kinetics, and was never experimentally validated. Recently, Kumar et al. [27] proposed a similar variable impedance controller, where ankle stiffness and damping were defined as polynomials in gait phase, and the coefficients defining the polynomials were identified using constrained least squares with an able-bodied kinematic and kinetic dataset. The authors utilized piecewise-constant equilibrium angles and demonstrated continuous stiffness and damping expressions that produced satisfactory gait with a postoptimization tuning protocol. This work was later extended to include variable inclines and a phase variable parameterization of stiffness and damping based on the phase portrait of the thigh angle and its integral [28]. However, this phase variable is known to have challenges with nonsteady walking [31], and changes in impedance associated with walking speed were not considered. The authors of [28] also note that their method of identifying the impedance parameters is nonconvex, which does not guarantee a globally optimal solution [32] for their controller.
This article addresses these limitations by presenting a new phase-based, task-adaptive walking controller built on a hybrid combination of continuously variable impedance control during stance and kinematic control during swing (Fig. 1). First, we present a convex, data-driven framework to calculate stance phase joint stiffness, damping, and equilibrium angle as continuous functions of gait phase, walking speed, and incline from an able-bodied dataset [24] (Section III). Paired with an analogous model of swing joint kinematics [23], our hybrid controller adapts behavior across varying tasks based on real-time phase, speed, and incline estimates (Section IV). Next, we present an improved phase variable that avoids kinematic singularities and is robust to the diverse family of thigh trajectories associated with variable-task walking. Then, we perform validation experiments with two above-knee amputee (AKA) participants, demonstrating that the adaptive controller produces biomimetic trends in joint kinematics, kinetics, work, and cadence across varying tasks (Section V). Finally, we show that our presented controller meets or exceeds the performance of a hand-tuned benchmark FSM impedance controller in most tested metrics, suggesting that our optimized kinematic and impedance Fig. 1. A block diagram of the hybrid kinematic impedance controller presented in this work. Real-time estimates of gait phaseŝ and taskχ define desired joint impedance parameters K, B, θ eq , and joint angles θ d using data-driven models. Depending on if the user is in stance or swing, the torque commands τ are calculated using either an impedance controller or a position controller, respectively. models sufficiently capture the key biomechanics of variabletask walking.

II. RELATED WORK
Many researchers have attempted to lessen the manual tuning burden of FSM impedance controllers in previous work. One common approach is to limit impedance control to the stance phase of gait and use kinematic control in swing phase, similar to our proposed architecture. Though many have used this hybrid architecture without a phase variable [6], [33], [34], [35], [36], [37], [38], [39], [40], relatively few have used it with one [28], [41]. Phase variable parameterization can be helpful because it allows continuous regulation of the dynamic interaction between the user and the ground during stance and provides the user with indirect volitional control over foot position during swing [31].
Additionally, some researchers have used biological quasistiffness curves calculated from able-bodied data [33], [37], [38], [41], [42] in lieu of hand-tuned impedance parameters. While the work in [33], [37], [38], and [42] enabled variable-speed walking and the work in [38] enabled obstacle crossing, these approaches were limited to level ground and relied on an FSM to switch between regions of the nonlinear quasi-stiffness curve during stance. Similarly, a quasi-passive ankle prosthesis presented in [43] enabled variable-incline walking with limited tuning by implementing a constant external quasi-stiffness relationship between the global shank angle and ankle torque. This external quasi-stiffness relationship was shown to be invariant across inclines during midstance, obviating the need for real-time incline estimation. However, this invariant relationship was limited to midstance and the controller relied on an FSM with manually tuned behavior for the remainder of the gait cycle. Further, as the control approach was developed for a passive prosthesis, it did not provide a method to increase net ankle work with increasing incline, which is an important characteristic of able-bodied walking [24]. Finally, this method was limited to ankle prostheses, and it is unclear whether the analogous external quasi-stiffness relationship for the knee during midstance is similarly invariant.
Other researchers have used reinforcement learning (RL) to automatically tune the impedance parameters online while a user walks, thus reducing the need for manual expert tuning [44], [45], [46]. Reward functions have been built on knee kinematic similarity to predefined trajectories or the observed contralateral knee's trajectories. However, these approaches were limited to the knee joint only and can require several minutes of walking before the optimal impedance parameters are identified. Further, the RL algorithms focused on kinematic features; the resulting kinetics and overall biomechanics were not investigated.
Non-impedance-based tuning-free controllers have also been developed. In [31], [47], [48], able-bodied kinematic profiles parameterized by a phase variable (i.e., virtual constraints) enabled tuning-free walking. This approach was extended for variable speed and incline walking in [49]. A similar controller was suggested for stair ascent [50]. However, the purely kinematic control paradigm tended to display nonbiomimetic joint torques during stance. In addition, the tuning-free knee-ankle prosthesis controller presented in [22] used an electromyography signal from the biceps femoris to control knee torque. The ankle impedance controller used a constant stiffness and damping with an equilibrium angle calculated from the knee angle. This controller enabled walking, sitting, squatting, and lunging, but was not demonstrated on different slopes.
Finally, our work is most closely related to the phase-varying impedance controller derived from able-bodied data in [28], as discussed in Section I. However, our approach is distinct in multiple important ways. First, our convex optimization formulation provides an approximation of the globally optimal impedance parameter functions. In addition to global optimality, our convex formulation can be solved in polynomial time [32] to facilitate future work on real-time optimization using user data or clinician preference (e.g., [51]). Second, our variable impedance model includes a continuous function for equilibrium angle, mirroring the continuous progression of biological joint dynamics [12], [13], [14]. Third, our variable impedance model is further parameterized by walking speed, which is critical to reproducing normative gait energetics [48]. Fourth, we estimate the task variables in real time, making the system fully autonomous. Fifth, we use a phase variable that is more robust to variable speed and incline behavior than prior phase variable definitions [28], [31], [48], [49]. And sixth, we demonstrate that our approach produces biomimetic trends in joint kinematics, kinetics, work, and cadence for two novel AKA participants over a range of tasks without any manual impedance tuning.

A. Model Framework
To use impedance control for the stance phase of the gait cycle in a continuous, phase-based control framework, we require a model analogous to the kinematic model developed in [23] that describes how the impedance parameters (K, B, and θ eq ) should evolve. Specifically, we require the impedance parameter model to be continuously parameterized by both gait phase s and task χ = (ν, γ), where task is defined by the current walking speed ν and ground incline γ over the ranges 0.8 ≤ ν ≤ 1.2 m/s and −10 ≤ γ ≤ 10 deg.
A model that meets these criteria can be constructed from a linear combination of phase-varying polynomials, where the linear combination weights vary with the task. Polynomial functions of phase are useful to model parameter progression during stance because they are simply parameterized and can represent arbitrary aperiodic signals. We use fourth-order polynomials (d = 4), as they allow sufficient flexibility to model the parameter behavior without overfitting. Once the appropriate polynomial functions are identified for individual tasks in a dataset, bilinear interpolation can be used to create a unified, continuous model with task and phase inputs. First, we define task-specific polynomial functions that represent how the parameters vary during stance for a set of fixed tasks. For convenience, let s st be the stance phase (i.e., s st = s/s TO , where s TO is the phase at toe-off (TO)). Then, the impedance parameters for the pth fixed task χ p are . . , d}} is a set of constant coefficients. Then, the coefficients κ νγ defining the impedance parameter trajectories for an arbitrary task (ν, γ) are calculated through bilinear interpolation of its four nearest neighboring tasks κ ν n ,γ n , where ν n ∈ {ν 1 , ν 2 }, γ n ∈ {γ 1 , γ 2 }. For all j elements in κ νγ , this interpolation is Finally, using κ νγ and (2) evaluated at the current stance phase s st , the impedance parameters are calculated. Therefore, the model is fully defined once each task-specific set of coefficients κ χ p is calculated.

B. Model Fitting
We use an optimization-based approach to fit the model to a dataset of able-bodied walking [24]. The dataset contains kinematic and kinetic joint trajectories recorded from ten participants walking at steady state at 15 distinct points in the task space (i.e., γ ∈ {−10, −5, 0, 5, 10} deg, ν ∈ {0.8, 1.0, 1.2} m/s). Therefore, for each task χ p , we construct an optimization problem to identify the set of impedance parameter coefficients κ * χ p that, when used in (1) and (2), best reproduced the mass-normalized joint torques τ in the dataset, given the dataset kinematics (θ,θ) over all n data points at χ p : 1) Solution Approximation: As written, (4) is difficult to solve, as the product K χ p θ eq,χ p is nonlinear in the unknown parameters, and the overall objective function is nonconvex. To avoid this issue, we solve a similar, convex problem and use its solution to approximate a solution to (4). First, we combine the product of K χ p and θ eq,χ p into a new, higher order polynomial δ χ p with independent coefficients δ ip : By treating the δ ip terms as independent from the k ip terms, the impedance equation forτ becomes linear in the unknown parameters k ip , b ip , and δ ip . We can then write the modified optimization problem as a standard quadratic program (QP), defining a new argument vector x ∈ R 4d+3×1 as x = [k 0p , . . . , k dp , b 0p , . . . , b dp , δ 0p , . . . , δ 2dp ] .
Let α j ∈ R 4d+3×1 be defined for each data point j as Then, the objective function L(κ χ p ) from (4) becomes where 2) Constraints and Regularization: To prevent overfitting, we added a diagonal regularization matrix R = diag(λ) ∈ R 4d+3×4d+3 to H to penalize the L 2 norm of x. The nth diagonal entries in R corresponding to the regularization weights on k i and b i were λ n = 1e −5 while λ n = 1e −2 for the δ i terms. These hyperparameters were chosen prior to the experiments in order to produce a smooth model that captured general behavior instead of overfitting to the training dataset.
Next, we added a constraint matrix A to ensure that K χ p (s) and B χ p (s) remained within ranges that were both physiologically realistic and feasible for the prosthesis to render in a stable manner. Namely, K χ p (s) was constrained above 1.5 Nm/rad/kg and B χ p (s) was constrained between 0.01 and 1.0 Nms/rad/kg. In addition, AKA participants in preliminary experiments noted that a low stiffness at heelstrike (HS) was unsettling, as they were accustomed to a locked knee during early stance with their take-home prostheses. Therefore, a minimum HS stiffness constraint of 3.0 Nm/rad/kg was added to increase participants' confidence that the prosthesis was ready for weight acceptance at HS.
To enforce these constraints, we discretized stance phase into n j points in the range [0,1]. We constructed a constraint matrix A ∈ R 3n j ×4d+3 from submatrices A s ∈ R n j ×d+1 as A column vector b ∈ R 3n j ×1 contained n j copies of the minimum stiffness and damping and maximum damping values, with the first term modified for the HS constraint: Finally, we arrived at the full QP, with the positive offset torque (sum-of-squares) in (8) neglected without loss of generality: We solved this QP for each subject and task χ p combination in the dataset (N = 150) using the MATLAB Optimization Toolbox (R2021b, MathWorks, Natick, MA, USA). Then, we approximated the solution to the original problem (4) by projecting the rational function δ χ p (s st )/K χ p (s st ) = θ eq,χ p (s st ) onto a dth-order polynomial. We assumed the polynomial order was sufficiently high to approximate the rational function δ χ p (s st )/K χ p (s st ) without significant information loss. This assumption was validated by the model's low reconstruction error, detailed in the next section. Then, for each task χ p , the intersubject mean set of coefficientsκ χ p was calculated for use as the final model. Trials that did not well represent the data, measured by a variance accounted for below 75%, were discarded as outliers prior to averaging. Fig. 2 shows the calculated impedance parameter model projected onto a speed of 1 m/s, which was produced by evaluating (2) and (3) withκ χ p . To quantify the impedance parameter model's reconstruction error, we calculatedτ for the knee and ankle over each trial in the dataset using the model:

C. Modeling Results
Then, we calculated the root mean squared error (RMSE) in joint torque over all subjects for each task χ p in the dataset and normalized by the dataset torque's standard deviation for χ p . This dimensionless metric, which we call normalized reconstruction errorĒ, describes how many standard deviationsτ is from the mean dataset torque trajectories, on average. Normalized reconstruction errors below 1.0 indicate that the model is able to predict joint torque to accuracy levels similar to able-bodied intersubject variation. Averaged over all tasks, the knee and ankle normalized reconstruction errors wereĒ k = 0.78 ± 0.11 andĒ a = 0.58 ± 0.09, respectively.

IV. HYBRID KINEMATIC IMPEDANCE CONTROLLER
The proposed hybrid kinematic impedance controller (HKIC, Fig. 1) is an evolution of the purely kinematic controller presented in [49]. In the HKIC, real-time phase and task estimates provide inputs to the impedance model developed in Section III during stance and the kinematic model developed in [23] during swing to provide reference joint behavior. Impedance and position controllers enforce the respective model outputs, described below. Once configured with the user's mass and leg segment lengths, the controller operates autonomously, requiring neither manual impedance tuning nor external knowledge of the terrain. The following sections discuss each component of the HKIC in turn.

A. Task-Invariant Phase Estimation
An estimate of the user's progression through the gait cycle is required in order to synchronize the control outputs with the user's gait. An ideal version of this estimate (termed a phase variable) increases from 0 to 1 at a constant rate between each HS [52]. Similar to [31] and [49], the HKIC's phase variableŝ is calculated using a piecewise-linear mapping of the user's global thigh angle θ th , which has a roughly sinusoidal trajectory (see

Fig. 3(a)
). This angle is measured directly using an Inertial Measurement Unit (IMU, 3DM-CX5-25, LORD Microstrain, Williston, VT, USA) mounted on the proximal end of the prosthesis's knee joint. Mounting the IMU on the prosthesis instead of the person ensures a rigid connection to prevent slipping and vibration, which are commonly associated with soft tissue connections. Proper alignment of the prosthesis by a prosthetist ensures correct alignment of the IMU. The θ th -based method of phase estimation is preferable because it allows the user to start and stop the gait cycle at will and enables nonrhythmic behavior [31].
However, previous iterations of the θ th -based phase variable did not work well for variable-task locomotion because of assumptions made about the shape of the θ th trajectory. For example, [31] and [49] assumed that the θ th trajectory could be divided into two monotonic sections. While this assumption holds fairly well for level-ground and incline walking, it is invalid for steep declines [24] (Fig. 3(a)). Previous methods produced inaccurate, saturated phase estimates for such cases [49]. Further, previous methods did not account for periods of low thigh angular velocity (i.e., when the hip joint is most extended or most flexed), leading to pauses in the phase estimate and subsequent problems in the controller behavior [41], [49]. Therefore, in this work, we relax previous assumptions and add flexibility to the phase variable to better parameterize the gait cycle based on the diverse θ th trajectories observed in variable-task locomotion. First, we introduce short periods of feedforward phase progression that allowŝ to maintain a constant positive rate even when the thigh angular velocity is low, which enables a powerful and biomimetic push-off. Second, we add states to account for thigh trajectories that have more than two monotonic sections (especially common during ramp descent) to prevent excessive phase saturation and gait desynchronization. Third, we introduce a technique to improve the linearity ofŝ, correcting for previous steady-state nonlinearities and thus making it closer to an ideal phase estimate. For brevity, the mathematical details for these improvements are presented in Appendix A.
To illustrate the benefits of the new phase variable over its predecessor [31], [49], we conducted a simulation using thigh kinematic data from [24] (Fig. 3(a)). For each trial of treadmill walking in the dataset, we calculated the phase variable using both the new method (Appendix A) and the previous method described in [49]. For each incline, we averaged the phase trajectories over all strides, participants, and walking speeds, shown in Fig. 3(b) and (c). Notably, the new phase variable eliminated the phase estimate pause associated with maximum hip extension that was observed with the previous phase variable. The new method also reduced the early saturation seen in the previous phase variable, which was particularly prominent at steep ramp declines. Finally, the new method demonstrated improved linearity, particularly during midstance. Compared to an ideal linear phase trajectory, the new method showed 6.25% RMSE with R 2 = 0.990 while the previous method showed 7.48% RMSE and R 2 = 0.976 over all tasks.

B. Task Estimation
In addition to the phase estimate, the HKIC requires an estimate of the user's current taskχ, which is calculated at each TO during steady walking. The estimation methods below are based on [49], with modifications to improve performance. Both estimates update once per stride and are filtered with a moving average over three strides to account for stride-to-stride variation. Although filtering introduces a time delay in the task estimates, experiments in [49] demonstrated that this limitation does not prevent the user from continuing to walk while the estimates converge.
1) Walking Speed: The user's speed is estimated using a three-link leg model, comprising thigh, shank, and foot links, similar to [33], [53], [54], [55], and [49]. Using forward kinematics and inputs from the joint encoders and the thigh IMU, we calculate the Cartesian locations of the heel and toe relative to the hip joint, respectively, given by x heel and x toe . At each TO event, the forward progression of the hip relative to the foot's  [24], where positive angles correspond to hip joint flexion, and (b) and (c) the resulting phase variable trajectories at different inclines. (b) Plot shows the trajectories calculated the previous method described in [49]. (c) Plot shows the trajectories calculated using the new phase variable presented in this work. The new method shows no phase pause near push-off and improved linearity, especially at the point of maximum hip extension. point of contact with the ground during the previous stance phase is calculated as where x HS− heel is the value of x heel from the previous HS. Similarly, by assuming a symmetric gait, the forward progression of the hip relative to the contralateral foot's ground contact point over a swing phase is approximated at each HS as where x TO− toe likewise is x toe from the previous TO. Then, we calculate the total forward progression over the gait cycle as d st + d sw + foot , where foot is a constant accounting for the length of the prosthetic foot. Finally, walking speed is estimated by dividing forward progression by stride time.
2) Incline: The ground inclination is estimated by the global angle of the foot θ f when the foot was flat on the ground, similar to the methods presented in [16], [49], [54], [55], [56]. As it is undesirable to add an extra inertial sensor to the foot, we calculate this angle from the thigh IMU using forward kinematics, along with a correction for foot bending. Prosthetic feet are designed to deflect for energy storage [57], so foot deflection significantly impacts the incline estimate. Offline testing with our prosthesis's foot [39] (Lo Rider, 1E57, Ottobock, Duderstadt, Germany) showed that deflection was correlated with the bending moment in the sagittal plane m y . An on-board six-axis load cell (M3564F, Sunrise Instruments, Nanning, China), located at the distal end of the ankle joint, measures this moment directly. Then, θ f is calculated as where k f is the linear bending coefficient, θ k is the relative knee angle, and θ a is the relative ankle angle. All joint angles are measured positive in flexion and are zero when the user stands upright. The constant offset term θ 0 f accounts for the angular difference between the prosthetic foot, the cosmesis, and the sole of the shoe.
To determine when the foot was flat on the ground, the center of pressure in the foot reference frame cop is calculated using the load cell, similar to [49]. We consider the foot to be flat when 7.5 ≤ cop ≤ 12 cm from the ankle joint, which corresponds to the ground reaction force acting between the middle and the ball of the foot. During this period, θ f is averaged to produce the incline estimate for the stride.
C. Impedance and Kinematic Controllers 1) Stance Impedance Controller: During stance, a variable impedance controller is used to calculate joint torques. First, the stance phase estimateŝ st is calculated bŷ wheres TO is the expected value of the phase variable at TO (see Appendix A for details). Usingŝ st andχ, joint stiffness K, damping B, and equilibrium angle θ eq are calculated using (2) and (3) and the model developed in Section III. Then, the joint torque during stance is calculated with the following impedance control law, scaled by user mass m: 2) Swing Kinematic Controller: A proportional derivative (PD) controller uses constant gains k p and k d to directly track desired joint angle trajectories. This is in contrast to the equilibrium angles of the impedance controller, which do not necessarily align with the normative joint angles. A continuous model of able-bodied joint kinematics [23], generated using data from [24], provides the desired trajectories defined as where b k (s) are Fourier series and c k (χ) are Bernstein basis polynomials. Similar to the impedance model, (19) is evaluated in real time usingŝ andχ. Then, the PD torque command during swing, τ sw , is given by 3) Stance to Swing Transition Smoothing: A time-varying weight w sw ensures a smooth transition from impedance control to position control. Because impedance control may allow the joint angles to vary from their nominal trajectories depending on how the user loads the prosthesis, this smoothing is critical to avoid discrete changes in joint torque. At TO, w sw increases from 0 to 1 over 0.25 s for the knee and 0.05 s for the ankle. The ankle smoothing is faster because close tracking of the ankle kinematics during early swing is important for avoiding toestubbing. The actual output to the joint motors is given by Because the equilibrium angles at HS are close to the kinematic references at the end of the gait cycle, no smoothing is necessary for the swing to stance transition. Close examination of (21) shows that for a brief period just following TO, minimal control action is applied to the joints. This is acceptable because the low-impedance actuators used in our prosthesis [39] allow the joints to continue moving along their current trajectories according to their passive dynamics without control input. Passive early swing knee and ankle dynamics have been shown to produce human-like gait [58], [59], and these passive dynamics may contribute to the biomimetic behavior of the controller.

V. AMPUTEE PARTICIPANT EXPERIMENTS
Experiments with two AKA participants were performed to investigate the ability of the HKIC to produce biomimetic gaits over variable tasks. To benchmark the HKIC's performance against another well-known controller, we also implemented a standard, piecewise-constant FSM impedance controller and tuned it for each participant. The participants completed the experimental protocol once with each controller, detailed below. Photos of the experiment are shown in Fig. 4, and video recordings are available online as supplemental media.

A. Benchmark FSM Impedance Controller
A benchmark finite state machine controller (FSMC) was designed based on the variable-incline FSM impedance controller presented in [16], with an additional stance state and modified transition criteria to improve performance (see Appendix B for details). This controller was chosen as a benchmark because of its simple construction, widespread usage [1], and ability to create biomimetic walking gaits when appropriately tuned [16]. While more sophisticated variants of the FSM impedance control paradigm have shown stronger results, such as those that modulate the impedance parameters based on joint angles or prosthesis axial force [6], [20], [21], [60], [61], the original version from [16] provides a valuable benchmark for comparing novel controllers because its performance and limitations are widely understood [1], [9]. Further, many modern controllers still use FSM impedance control in some if not all sections of the gait cycle [18], [19], [34], [35], [36], [44], [45], [46], [47], [62], so understanding the HKIC's performance relative to the FSMC is scientifically relevant.
The FSMC had five discrete states throughout the gait cycle, each with its own set of constant impedance parameters and transition criteria. Similar to the methods discussed in Section I, these parameters needed to be hand-tuned by an expert researcher in order to produce the desired gait. To enable walking at various inclines, three sets of tunable impedance parameters and transition criteria were instantiated for each joint (i.e., one set for level ground, one set for declines, and one set for inclines). The controller selected between impedance parameter sets based on the estimated inclineγ (Appendix Fig. 15(b)). In total, the FSMC required 96 tunable parameters, including 45 impedance parameters per joint and six FSM transition criteria.

B. Experimental Methods
Two AKA individuals participated in the experiment, with attributes shown in Table I. A third participant was enrolled but was unable to complete the protocol due to excessive swing-phase lateral whipping caused by prosthetic misalignment. Although we worked with the prosthetist to correct the alignment multiple times, the prosthesis would become misaligned again after a short walking bout, possibly due to a combination of his prosthetic socket and weak femur musculature [63]. We suspect that the large distal mass of the robotic prosthesis also exacerbated this problem. Due to this issue, we only present data from the remaining two subjects in this article.
The experimental protocol was approved by the Institutional Review Board of the University of Michigan (HUM00166976), and the participants wore a ceiling-mounted safety harness while walking on the treadmill. For the experiments, the presented HKIC and the comparison FSMC were implemented on a backdrivable, powered knee-ankle prosthesis, shown in Fig. 4 and described in depth in [39]. This prosthesis features quasi-direct drive actuators that enable open-loop joint impedance control.
A licensed prosthetist fits the prosthesis to the participants and ensured proper alignment. The participants were instructed on the expected high-level behavior of both controllers and given time to acclimate to each controller while walking overground within parallel bars. Importantly, the participants were not told which controller was expected to perform better during the experiment. Following this overground acclimation, five trials with each controller were conducted on an in-ground treadmill (Bertec, Columbus, OH, USA). For safety, instrumented handrails were provided on either side of the treadmill. The participants were encouraged to limit body weight support on the handrails to maximize the realism of the experiment, which was verified by handrail force data. Participant P1's mean handrail usage was under 12% bodyweight and participant P2 frequently chose to use only one handrail (Fig. 4(c)-(d)).
The first three trials investigated the performance of the HKIC and the FSMC during steady walking at different speed and incline combinations. Each trial focused on a range of small task deviations (±2 deg, ±0.2 m/s) around one of three baseline tasks: χ = (0 deg, 1 m/s), χ = (5 deg, 1 m/s), and χ = (−5 deg, 1 m/s). We refer to these steady-state task trials as SS-Level, SS-Incline, and SS-Decline, respectively. For the SS-Incline trial, speed was limited to 1.1 m/s to ensure that the participants could safely perform the trial.
The steady-state task trials began with an acclimation period, where the participants walked at the baseline task until feeling comfortable. During this time, the FSMC was tuned by the authors of this work to produce a natural gait, incorporating feedback from the participants and the prosthetist. The authors have significant experience tuning impedance controllers [21], [47], [61]. Tuning continued until the authors, prosthetist, and participant were satisfied with the resulting natural gait (see supplemental video, available online). The time taken to tune the FSMC was recorded. Note that no tuning was done for the HKIC. After tuning and acclimation, the participants walked on the treadmill as it cycled through each of the five tasks near the baseline task, each commanded for 45 s. In these trials, true task feedback was provided to the controllers so that any errors in the task estimates did not affect the results.
The tuning, acclimation, and testing procedure above was repeated for each baseline task. These baseline tasks were chosen to be far apart in the task space in order to sample a wide range of tasks without deviating too far from any one of the FSMC's tuning points. Fig. 5(a) shows the recorded task-space profiles from the treadmill for each trial, where the black dots indicate each commanded task.
The latter two trials consisted of more rapid task changes to investigate each controller's behavior during continuous task variations rather than at steady state and over a wider range of tasks. Also during these trials, the controllers received no real-time knowledge of the task from the treadmill, investigating the autonomous capability of each controller to operate over variable tasks. Both controllers utilized the same task estimation methods (Section IV-B). The FSMC transitioned between the tuned impedance parameters sets based on the estimated incline (Appendix Fig. 15(b)). In these two trials, one with inclines (CV-Incline) and the other with declines (CV-Decline), the treadmill Because the treadmill required time to change task, smooth task trajectories with continuous variations were generated, shown in Fig. 5(b).

1) FSMC Tuning Time:
For the two participants, the FSMC required on average 30 min of tuning to produce normative gaits for the three baseline tasks. On average, the level-ground task required 11 min, the incline task required 15 min, and the decline task required 5 min. Participant-specific tuning times and tuned FSM parameters are listed in Table IV in the Appendix. Trends in the tuned parameters included higher stiffness values during stance than in swing and highly varying knee equilibrium angles across tasks. The observed gait was also noted to be quite sensitive to the tunable FSM transition criteria. Significant variance in the required tuning time for the different tasks was also observed.
2) Steady-State Trials: The kinematic and kinetic trajectories produced by the HKIC during the steady-state trials highlight its ability to reproduce normative biomechanics over variable tasks (Fig. 6). Bilinear interpolation was used to generate the able-bodied reference trajectories for tasks between those reported in the dataset [24]. The observed HKIC trajectories show strong similarity to the able-bodied references, particularly at the ankle. Knee moments are the most different relative to able-bodied for both the HKIC and the FSMC. The separation and trends seen in the HKIC closely resemble those observed the able-bodied data, suggesting appropriate adaptation in response to variable-task walking.
We quantified the similarity between the observed and ablebodied trajectories during stance and swing, showing that the HKIC produced a low RMSE in most metrics (Fig. 7). Stance and swing were treated separately to isolate the performance of the novel impedance parameter model (Section III), as it was only used during stance. The first 15 s at each task were neglected to allow time for the treadmill to reach steady state. Unless otherwise specified, we present interparticipant averages and calculate standard deviations using lumped participant strides.  [24] are also shown for reference. The HKIC produced smooth kinematic variations with incline changes as well as increasing knee flexion and ankle push-off torque with increased speed, resembling the able-bodied trajectories. Fig. 7. Interparticipant RMSE in the observed kinematics (left) and kinetics (right) relative to able-bodied walking data for both the HKIC and FSMC during the steady-state task trials. The error bars represent ±1 standard deviation over lumped participant strides. The HKIC demonstrated lower mean error than the FSMC in seven of eight metrics, with particular improvements at the ankle joint. The high knee kinematic error in swing for the HKIC is the result of intentional early extension to promote user confidence that the prosthesis was ready for weight acceptance. Fig. 8. Interparticipant average cadence for the steady-state task trials as functions of speed for different ramp inclinations: ramp descent (left), level ground (middle), and ramp ascent (right). Error bars represent ±1 standard deviation over lumped participant strides. Both controllers show similar cadence trends as the able-bodied reference (AB) calculated from [24], with increasing step frequency with increasing speed. Overall, the participants preferred longer strides relative to able-bodied, which may be due to the larger mass of the powered prosthesis.
Individual RMSE values for each participant were similar to the interparticipant averages and are available in the Appendix (Table V). The low RMSE values suggest that, in addition to replicating normative trends as task varied, the HKIC produced kinematics and kinetics that were close to the reference values. Further, the HKIC's performance was as good as or better than the hand-tuned FSMC's performance in seven of the eight metrics. The high knee kinematic error during swing can be attributed to the intentional early knee extension meant to improve user confidence (see Appendix A2) and it did not result in adverse gait effects.
Interparticipant spatiotemporal gait metrics also showed similarity to able-bodied data [24]. Both controllers elicited lower cadence gaits (equivalently longer stride length gaits) compared to able-bodied, but show generally similar trends of increasing cadence with walking speed (Fig. 8). In addition, the stance time symmetry ratio r STS was calculated using the ground reaction force data, defined as the ratio between the average prosthetic stance time and the contralateral limb stance time. The mean and standard deviation over the steady-state trials of both participants were r STS = 0.902 ± 0.017 for the HKIC and r STS = 0.892 ± 0.016 for the FSMC. Note that due to a recording error, symmetry data were not available for participant P1's HKIC SS-Decline trial. Both controllers produced  II  TASK ESTIMATE RMSE OBSERVED DURING THE CONTINUOUSLY VARYING  TASK TRIALS AVERAGED OVER LUMPED PARTICIPANT STRIDES a slightly more symmetric gait than average AKA participants with passive prostheses (r STS = 0.784, reported in [64]), but less symmetric gaits than able-bodied people (r STS = 1.02, reported in [65]). The HKIC also produced trends in joint work across variable tasks that were consistent with able-bodied data (Fig. 9). As one of the benefits of impedance control is the ability to control energy exchange with the environment [10], the HKIC should be able to replicate this biological behavior. The HKIC showed similar trends as the able-bodied data, with a linear increase in net work performed with increasing incline, particularly at the ankle (increase of 0.0337 J/kg/deg, R 2 = 0.982). For comparison, able-bodied ankle work increases linearly at 0.0335 J/kg/deg with R 2 = 0.987. The HKIC also increased total work with increasing speed in a manner consistent with able-bodied data, though the work differences between the slow and fast speeds are minor for the able-bodied reference. In contrast, the net work performed by the FSMC decreased with speed and appeared discretized to three levels with respect to inclines, corresponding to its tuned tasks. Interestingly, the HKIC and FSMC showed less energy absorption at the knee during declines, which may reflect the habitual aversion to early stance knee flexion commonly observed in AKA populations [66], [67].

3) Continuously Varying Trials:
The continuously varying trials demonstrated the HKIC's ability to autonomously adapt behavior to the sensed walking speed and ground incline. The kinematic and kinetic errors were calculated in a similar manner for the continuously varying task trials, though this time including strides that occurred during task transients. Fig. 10 shows the interparticipant average error trajectories at both joints for the CV-Incline trial, calculated as the able-bodied references subtracted from the observed values. Table V in the Appendix details the participant-specific stance and swing kinematic and kinetic RMSE for both the CV-Incline and CV-Decline trials. Aside from the late-swing knee kinematics (discussed above), the HKIC shows low errors throughout the gait cycle, particularly at the ankle joint. Further, the magnitude of the FSMC's error is larger than the HKIC's for most of the gait cycle, highlighting the importance of the HKIC's continuously adaptive nature.
As both controllers received no external task input during these trials, the task estimates (and the phase estimate for the HKIC) contributed to the kinematic and kinetic errors. The task estimate RMSE, averaged over each stride and participant, is shown in Table II for each trial. Although the same task estimation algorithms were used with both controllers, FSMC showed higher incline estimate error, suggesting that differences in controller behavior may have impacted the incline estimate's efficacy. In addition, the average phase estimate trajectories produced by HKIC during the CV-Incline and CV-Decline Fig. 9. Interparticipant average prosthesis work per stride over variable (a) inclines and (b) speeds during the steady-state task trials. Error bars represent ±1 standard deviation over lumped participant strides. An able-bodied reference (AB) calculated from [24] shows that the HKIC demonstrated biomimetic energy injection, particularly through a linear increase in ankle work as incline increased, corresponding to 100.6% of the able-bodied rate. Both controllers showed less energy absorption at the knee during steep declines, suggesting that our participants may have had habitual aversions to early stance knee flexion. Fig. 10. Plot of the interparticipant average kinematic and kinetic error trajectories in the continuously varying incline trial, relative to able-bodied data [24]. The knee data are shown in the left column and the ankle in the right. Shaded regions represent ±1 standard deviation over lumped participant strides. Aside from intentional discrepancies in the late-swing knee kinematics (see Appendix A2), the HKIC showed low RMSE across the gait cycle throughout varying tasks, suggesting appropriately adapting biomechanics. trials were highly linear (mean R 2 = 0.989) and accurate (mean RMSE of 6.157%), even as speed and incline varied (see Fig. 11). However, the phase estimate saturated more often for participant P2 than participant P1, suggesting that participant P2's thigh trajectory was less similar to able-bodied trajectories than participant P1's.

A. HKIC Performance
This work presented a data-driven, phase-based walking controller for a powered knee-ankle prosthesis that autonomously adapted its behavior across a continuous range of walking speeds and inclines. To achieve this without manual impedance tuning, Fig. 11. Average phase estimate progression calculated in real time by the HKIC during the continuously varying task trials for participants P1 and P2. Shaded regions represent ±1 standard deviation. The linearity and consistency of the trajectories illustrate the phase variable's ability to adapt to continuous task variations and appropriately parameterize the gait cycle.
we used an able-bodied dataset to optimize for continuous stiffness, damping, and equilibrium angle functions that reproduced biological stance joint torques, given biological kinematics. In an initial offline analysis, we showed that our optimized impedance parameter model produced joint torques with across-task average normalized RMSE values of 0.78 and 0.58 for the knee and ankle, respectively. The low normalized RMSE suggests that the model captures the essential joint dynamics of able-bodied walking.
The subsequent experiments with two AKA participants demonstrated that the identified impedance parameter functions also rendered appropriate stance phase joint mechanics when used for real-time control in the HKIC. Other normative walking features were observed, such as increasing ankle work with increasing incline (Fig. 9) and increasing cadence with walking speed (Fig. 8). Although the kinematic and kinetic profiles produced by the HKIC had small differences relative to able-bodied data (Figs. 6 and 7), the participants exhibited qualitatively normal gait patterns over a wide array of tasks (see supplemental video, available online). Kinematic and kinetic trends emerged with variable speeds and inclines that were consistent with able-bodied data (Fig. 6), including appropriately varying peak ankle moments, stance ankle kinematics, and knee stance kinematics. Knee swing kinematics showed the highest error for the HKIC, which was expected because we intentionally allowed the phase variable to saturate early to ensure full knee extension prior to HS. Pilot testing showed that consistent full knee extension helped eliminate participants' problematic instinctive compensations and promoted confidence that the prosthesis was ready to accept weight (see Appendix A2). Small phase shifts result in large swing kinematic errors due to the large knee range of motion, and although the error values appear large, they did not interfere with the participants' gait or cause toe-stubbing.
Appropriate kinematic and kinetic adaptation are both practically and clinically important for the user. For example, knee swing kinematic adaptations enable the prosthesis to have the proper configuration at HS as incline varies. Without such adaptations, the user may, for example, toe-stub during swing when walking uphill with level-ground kinematics, or vice versa, experience too much flexion to enable HS at the desired time. Further, kinetic adaptation during stance enables increasing peak ankle moments for propulsion as incline and speed increase (Fig. 6). Improper joint kinetics can cause improper ground reaction forces, which can affect user balance [68]. Finally, appropriate kinematic and kinetic co-adaptation enables joint work adaptation, even in cases where both kinematics and kinetics deviate from able-bodied normative trajectories. For example, the HKIC's peak ankle moment at a 7 deg incline is slightly smaller than able-bodied (Fig. 6). However, a corresponding increase in peak plantarflexion angle allows the HKIC to maintain appropriate ankle work (Fig. 9). Biomimetic energy injection is important to prevent compensations from other joints and additional health problems [3], [4], [5], [41].
Qualitative remarks by the participants also testified to the biomimetic task adaptation of the HKIC. Participant P1 remarked while walking at the 7 deg incline that he did not feel like he was walking uphill, suggesting appropriate joint dynamics and energy exchange. Participant P2 remarked that he did not even notice that the treadmill had transitioned to the 2 deg decline and that he could "climb up much easier" while ascending steep inclines. These anecdotal remarks further support the claim that the HKIC adapts to changing tasks to produce normative able-bodied biomechanics, which could result in many practical benefits for the user. For example, the biomimetic energy injection at steep inclines ( Fig. 9) may allow users to walk uphill for longer before fatiguing.
While the FSMC's performance was not drastically worse than the HKIC's in the tested metrics, it required on average 10 min of tuning per tuned task. Although only three tasks were tuned for this study, practical deployment of the FSMC would likely require many more tasks to be tuned. For example, participant P2 noted that the FSMC was "kicking off way too hard" when going uphill at slow speeds, but was happy with its behavior at normal speeds, suggesting that more speedspecific impedance parameter sets could be beneficial. However, adding more tuning points is likely impractical in a clinical setting, especially without specialized equipment such as a variable-incline treadmill. Therefore, the HKIC's potential to produce biomimetic behavior over varying tasks without manual impedance tuning is a significant benefit.
For online implementation of the continuous impedance parameter model, gait phase needed to be estimated in real time. The improved phase variable behavior observed in simulation in Section IV-A was confirmed in the participant experiments. Fig. 11 shows how the phase variable eliminated the previously observed phase pause near push-off. The result of this monotonicity is visible in the kinematics of Fig. 6, as there is not a pause in the kinematic trajectories near push-off, which was observed previously in [49]. Further, the general linearity of the average phase trajectories in Fig. 11 (mean R 2 = 0.989) is improved compared to [49]. Because both the impedance and kinematic models in HKIC assume a perfectly linear phase estimate, the observed linearity keeps the model outputs of the controller synchronized with the user's gait. In addition, Fig. 11 shows that participant P2's phase variable saturates at 1 earlier in the gait cycle than participant P1. This occurs because the methods used to estimate the thigh trajectory features (Appendix A2) prioritize phase saturation (and subsequently full knee extension) to promote participant confidence. This early phase saturation suggests that P2 preferred for the knee to be fully extended earlier in swing, whereas P1 was satisfied with full knee extension occurring right before HS, as it does in able-bodied data.
The task estimates are other critical components required for walking over continuously varying tasks. It is seen in Table II that the error in the speed estimate was fairly consistent over the trials, with RMSE between 0.10 and 0.12 m/s for both controllers. This error is likely due to a slightly asymmetric gait, which violates the assumptions made in the speed estimator's formulation. Gait asymmetries may be the result of our participants' habitual compensations, socket comfort, or the significant mass difference between the robotic prosthesis and participants' passive prostheses. Interestingly, the incline estimate produced lower error with the HKIC (0.61-0.66 deg) than the FSMC (1.03-1.81 deg). We speculate that the higher error in the FSMC is due to a feedback interaction between incline estimate errors and the impedance parameters. Due to the discrete switching behavior of the impedance parameters (Appendix Fig. 15(b)), a small incline estimate error can result in large changes in prosthesis behavior and may affect the θ f and cop progressions. Therefore, the continuous nature of the HKIC may be preferable, as it does not display discrete changes in behavior with small changes in task inputs.

B. Limitations and Future Work
The HKIC and this study were not without limitations. Our experiment provided a somewhat limited view of the HKIC's behavior, as it involved only two participants, each with only one experimental session (although both had prior experience walking with the prosthesis). We expect that the data-driven impedance parameter model identified in Section III will yield similar performance for a wide array of participants, as it was created without a priori knowledge of the participants or their preferences. Preliminary studies of able-bodied users testing the HKIC over varying tasks suggest that this assumption holds [69]. However, this assumption should be validated in future studies with wider AKA participant pools, as the HKIC's ability to generalize to the full AKA population remains unknown. Further, the performance of the HKIC should also be investigated when implemented on different hardware platforms to better validate the framework. While the HKIC could in theory be implemented on any powered prosthetic leg with the ability to render a variable joint impedance, prostheses with nontrivial actuator transmissions (e.g., [6], [48]) or series elasticity (e.g., [21], [61]) would require actuator characterization [70] and, in some cases, closed-loop torque control to accurately render variable joint impedance.
In addition, it is possible that there are users for which the population average impedance parameters are not optimal. There may also exist other impedance parameter functions that produce normative biomechanics, as the human sensorimotor system is highly adaptable [71], [72], [73]. Participants' sensitivity to changes in the impedance parameter model could be investigated in future studies. In addition, there may be factors other than those considered in this work that distinguish the ideal parameter functions, such as user preference. A study investigating users' preferred stiffness in ankle prostheses showed that the preferred joint stiffness varies by user [74], which is likely true for AKA participants as well. For example, participant P2 noted that the knee felt "squishy" when ascending steep slopes and that he would have preferred it to be stiffer. While one of the major advantages of HKIC is that it required no manual impedance parameter tuning, it is currently limited by the lack of an ability to customize to an individual's preferred behavior. Future work will investigate methods to incorporate user preferences in the impedance model, such as weighting the optimization with a single baseline personalization for level-ground walking, as suggested in [75]. This baseline personalization could be gathered using tools in a standard clinic [51], maintaining the minimal-tuning nature of the controller.
Our study also did not investigate discrete changes in ground slope, which may be encountered during daily ambulation (e.g., wheelchair ramps) and should be handled by a variable-task walking controller. While the work in [49] showed that the HKIC's incline estimation algorithm is stable under discrete incline changes, our task estimation methods are limited by their discrete "once-per-step" update nature. Because we detect incline during midstance, the user must, at a minimum, be able to complete the first half of stance phase with the previous stride's task estimate. This is particularly problematic during discrete transitions between steep inclines and level-ground walking because the HS kinematics vary drastically [24]. Future work involving anticipatory algorithms that update the task estimate based on sensed characteristics of the upcoming terrain [76], [77], [78] or user behavior [79], [80] may be necessary to alleviate this limitation.
Further, this work only investigated rhythmic walking over relatively long durations, though almost half of all walking bouts in community ambulation contain less than 12 consecutive steps [81]. One of the unique strengths of the presented phase variable is the ability to intuitively control nonrhythmic tasks [31]. Although this capability of the HKIC was demonstrated at the beginning and end of each trial in this study, it should be explored further and characterized in future studies involving rapid start/stop, lateral movements, and other behaviors that are prominent in agile locomotion. Such studies may also highlight the limitations of using the current impedance parameter model for nonrhythmic tasks. Although the participants were able to achieve start/stop behaviors in this experiment, additional able-bodied data may need to be included in the optimization to produce appropriate impedance parameters for other nonrhythmic tasks.
The continuously adaptive nature of the HKIC may also be a limitation in some circumstances. For example, participant P2 noted that while he appreciated that the HKIC always adapted to the current task, he also preferred the predictability of the FSMC. We plan to investigate methods to preserve the flexibility of the HKIC while increasing predictability in future studies. One way to improve the predictability of the HKIC may be to increase the training duration to allow the participants to better acclimate to and leverage the benefits the HKIC and the powered prosthesis. Perhaps the lack of early stance knee flexion in both controllers (Fig. 10) and the low knee energy dissipation ( Fig. 9) are less due to controller behavior and more due to the participant's habitual compensations developed through years of using a passive prosthesis [66], [67]. Future work may show that as the participants become more comfortable with a powered prosthesis and develop a stronger intuition for the HKIC's behavior, these gait features become more similar to able-bodied data.
Finally, there is much interesting work to be done investigating the relationship between biological joint impedance measured in empirical studies [12], [13], [14] and the impedance parameters used in impedance controllers. Mechanical impedance can only be characterized through perturbation studies, so the impedance parameters found by optimizing over nonperturbed gait data will not necessarily reflect these dynamics. We plan to study the effects of constraining the optimization with known empirical impedance values, as well as to investigate the HKIC's behavior during gait perturbations.

VII. CONCLUSION
This work presented a data-driven walking controller designed to work over a continuum of speeds and inclines. We developed continuous models of joint stiffness, damping, and equilibrium angle for an impedance controller using convex optimization. We also presented an improved phase estimation algorithm, showing increased monotonicity and linearity. Two AKA prosthesis users demonstrated the controller's ability to autonomously produce biomimetic behavior over continuously varying tasks during treadmill experiments. The experiments showed that, when compared with able-bodied data, the presented controller produced biomimetic trends in joint kinematics, kinetics, work, and cadence, indicating its ability to render appropriate joint mechanics as task varied.

A. Task-Invariant Phase Variable Algorithm
The new phase variableŝ is calculated through a series of linear equations with θ th as an input. An FSM controls when each equation is used. Although the FSM contains discrete states, the structure of the linear equations ensures thatŝ is continuous. Each equation is defined by quantitative features of the θ th trajectory, which are measured in real time. Table III lists the features' definitions and notations. First, we give the rationale for each  12. Average global thigh angle trajectory θ th (positive flexion) for 1 m/s 0 deg able-bodied walking, segmented by the phase variable FSM states. The phase variable is defined by linear mappings of θ th during S1, S2, and S4, and by a feedforward phase variable rate during S3 and S5. The feedforward rates for S3 and S5 are given by the average rate of change of the phase estimates during the preceding states, which correspond to periods of constant thigh angular velocity.
FSM state and its corresponding phase variable equation. Then, we present methods to estimate the thigh trajectory features in real time, as well as the steps taken to promote closed-loop stability of the phase estimate.
1) Phase Variable FSM: Consider the average θ th trajectory for an able-bodied individual walking at 1 m/s on level ground, shown in Fig. 12. The pertinent θ th trajectory features used in the phase estimate are labeled, as well as the standard timing of the FSM states. The overall structure of the FSM used to control the phase estimate is shown in Fig. 13.
The FSM begins in S1, occurring just after an HS event. During S1, θ th is linearly scaled as the hip joint extends from θ HS th to θ MHE th such thatŝ increases andŝ = s MHE when θ th = θ MHE th .
Mathematically, this is given bŷ s MHE in S1 and S2.
The FSM transitions to S2 at a phase estimate thresholdŝ 1→2 = 0.1, which typically corresponds to the point in the gait cycle where the θ th trajectory becomes linear. In S2,ŝ is calculated using the same linear relationship as in S1 (22), but is denoted as a distinct state because it represents a portion of the gait cycle where θ th (and thereforeŝ) has constant velocity. The average rate of change ofŝ during S2 (ṡ S2 ) is recorded for use in S3. The FSM transitions to S3 onceŝ 2→3 = 0.9s MHE , which typically corresponds to the end of the linear portion of the thigh trajectory, or ifθ th > 0. This second case rarely occurs during steady walking, but is an important path to S3 in the event of an unusually short stride.
S3 occurs during the section of the gait cycle where θ th reaches its minimum, and thus has a period of low angular velocityθ th . Previous work has shown that sections of lowθ th are problematic because they cause a pause in the phase variable trajectory [31], [41], [49]. This pause violates the assumption thatŝ increases monotonically and at a constant rate, resulting in incorrect kinematic and impedance model outputs. Therefore, during S3, we decoupleŝ from θ th and instead assume that phase continues progressing atṡ S2 : This feedforward phase progression continues until a TO event.
Although this approach limits the user's ability to stop phase progression during S3, such cases are unlikely because stopping would inhibit power delivery from the ankle during pushoff. Moreover, the underactuated dynamics of bipedal walking dictate that once the user's gravity vector passes anterior of the stance foot, the user must continue the gait cycle until the contralateral foot lands [68], [82]. Therefore, we expect the sacrifice in direct control of phase progression during this section of the gait cycle to be negligible. After TO, the FSM transitions to S4, where phase is again estimated via a linear scaling of θ th . This mapping is defined such thatŝ increases fromŝ TO toward s MHF as θ th increases: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
The FSM transitions to S5 when θ th is equivalent to the average of θ HS th and θ MHF th , which typically corresponds to the end of this linear section of the thigh trajectory.
Two problems typically occurred with the previous phase variable methods when θ th ≈ θ MHF th , which occurs during S5 in the new FSM. First, a pause inŝ would occur asθ th slowed and θ th approached θ MHF th , similar to the effect seen in S3. Second, the previous methods assumed that θ MHF th = θ HS th . In cases where θ MHF th > θ HS th , such as the trajectory shown in Fig. 12, the resultinĝ s would saturate prematurely. Excessive saturation in the phase variable can cause desynchronization between the prosthesis and the user, leading to problems such as toe-stubbing. This effect was most exaggerated during declined walking, as the difference between θ MHF th and θ HS th was most pronounced [24]. To avoid both excessive saturation and a phase variable pause, a feedforward phase progression is again enforced based on the average phase rate in S4,ṡ S4 :ŝ This feedforward phase rate continues until either an HS occurs orŝ = 1. If the user is walking consistently and the θ th trajectory feature estimates are correct,ŝ = 1 should occur simultaneously with HS, returning the FSM to S1. Ifŝ = 1 prior to HS, the FSM transitions to S6. S6 is primarily encountered if the user pauses at the end of the gait cycle, so it does not appear in Fig. 12. During S6,ŝ is again calculated using a linear scaling of θ th , giving the user volitional control ofŝ through θ th : This volitional control during S6 is important because it allows movements such as kicking and nonsteady leg swinging [31]. As in S5, an HS event returns the FSM to S1.

2) Thigh Trajectory Feature Prediction:
The θ th features used in (22)-(26) vary from stride to stride with changes in speed, incline, and natural gait variation. Some of these features are used in the phase estimate calculation before they occur in the gait cycle, specifically θ HS th , θ MHE th , θ MHF th , s MHE , and s MHF . For example, θ MHE th is used to calculateŝ during S1 and S2, but it does not typically occur until S3. Therefore, we predict these features in real time based on observations from recent strides. At controller initialization, estimates of the thigh trajectory features are calculated using able-bodied data [24] and updated as new strides became available. Bounds were enforced on all estimated feature values to reject atypical strides and avoid stride-to-stride oscillation in the estimates.
Previous work showed that care must be exercised when predicting features of the thigh trajectory to prevent unwanted interaction between the prediction algorithms and the user's gait progression. For example, Best et al. [49] observed that if a simple moving average was used to calculate θ MHE th , a divergent behavior occurred that resulted in the user taking progressively larger strides. To avoid this behavior, the kinematic features θ HS th , θ MHE th , and θ MHF th were estimated with moving average filters.
These filters recorded the feature values from the previous five strides and averaged the median three for θ HS th and the minimum three for θ MHE th and θ MHF th . These filters were chosen to best reject nonrepresentative strides, and the five-stride window balanced between filter response time and variance rejection.
Another closed-loop interaction was observed during pilot studies regarding the predictions of s MHE and s MHF . In cases when the feature predictors were updating following a rapid change in task, we observed rare strides whereŝ underestimated the true phase at the end of the gait cycle, causing the knee joint to not fully extend before HS. Participants instinctively responded by asymmetrically extending the late swing portion of the gait cycle to try force the knee to full extension. Moving average estimates of s MHE and s MHF , like those used for the kinematic features, caused s MHE and s MHF to decrease, resulting in further underestimation ofŝ on the subsequent stride. We suspect that participants behaved this way because they were accustomed to passive prostheses, which will collapse upon loading if the knee is not fully extended. Therefore, new prediction methods were developed for s MHE and s MHF that favoredŝ saturation over underestimation to combat this instinctive behavior. Let tŝ =1 be the first time during the stride thatŝ = 1. Then, the s MHE and s MHF estimates were calculated as The first quotient in each line of (27) is the true phase where θ MHE th and θ MHF th occurred. The second quotient is an upper bound on this true phase. We average the two so thatŝ favors saturation and full knee extension in late swing, avoiding the potential unstable feedback loop with the user's instinctive compensations. The results of (27) were likewise low-pass filtered with an infinite impulse response (IIR) filter to reject stride-to-stride variation and to prevent step estimate changes.
Finally, to calculate the stance phaseŝ st , the expected value ofŝ at TO,s TO , must be estimated. This was calculated with a minimum moving average filter ofŝ TO observed during previous strides, similar to the thigh trajectory features. The average window was nine strides long, as the TO phase exhibits slow changes with task. Like the thigh trajectory features,s TO was initialized from able-bodied data.
Note that some minor aspects of the thigh trajectory feature estimation algorithms were modified after P1's experiment to better accommodate adaptation for users with thigh kinematics that differ significantly from able-bodied, such as P2. Namely, the feature estimate bounds were added, the θ MHE th and θ MHF th filters were changed from moving average to moving minimum filters, and thes TO filter was changed from an IIR filter to a moving minimum filter. A post hoc simulation of P1's data before and after the minor adjustments showed only a 1.89% mean absolute difference in phase estimate between methods, suggesting that the changes would not have had an appreciable The phase variable is able to parameterize these nonrhythmic motions, allowing the participant to start and stop the gait cycle at will. effect on his results. Further, no distinguishable effects were observed in the able-bodied simulation (Fig. 3(c)).
3) Phase Variable Linearization: The phase variable described above produces a consistent phase estimate trajectory over each stride during steady walking. This consistency allows a linearization map to be formed in order to further improve the phase estimate. Once the θ th feature predictions converged to steady values, the average progression ofŝ was recorded for each steady walking stride and low-pass filtered to produce an average trajectory,s. The time constant of the IIR low-pass filter was chosen to be sufficiently slow (19 strides) such that the transients of the θ th feature predictors were rejected. As a further precaution, any saturated portions ofŝ were discarded prior to averaging, as they diminish as the θ th trajectory feature predictions converge.
The average phase was written as a function of true phase, given bys = σ(s). Although the shape of the thigh trajectory may cause σ(s) to be nonlinear, it is monotonic during normal walking. This implies that an inverse relationship s = σ −1 (s) exists, which can be applied to correct for nonlinearities in s. First, σ(s) was fit with a sixth-order polynomialσ(s) that was constrained with a minimum slope of 0.2. This minimum slope ensured strict monotonicity and numerical stability of the inverse. At each HS event,σ(s) was recalculated to incorporate the previous stride's effect ons. Then, the final, linearized phase estimate was calculated by applying the inverse mapσ −1 to the results of (22)- (26). Fig. 11 highlights the new phase variable algorithm's ability to parameterize variable-incline walking, as consistent phase trajectories were produced for both participants during the continuously varying task trials. The feedforward states S3 and S5 allowed for a positive phase rate, even when thigh velocity was low. In addition, the thigh trajectory features were appropriately estimated, allowing a consistent phase estimates that were independent of the variable thigh trajectories seen with varying inclines. The phase linearization algorithm ensured highly linear estimates, with mean R 2 = 0.997 for participant P1 and mean R 2 = 0.982 for participant P2. Finally, the volitional start/stop behavior of the phase variable originally shown in [31] was preserved. Fig. 14  shows the phase variable trajectory for four nonsteady bouts during the overground acclimation for participant P1, confirming its ability to parameterize nonrhythmic motion.

B. Benchmark FSM Impedance Controller
The FSMC, based on the FSM impedance controller presented in [16], was constructed to provide a benchmark with which to compare the HKIC. The flow of the FSMC's state machine is depicted in Fig. 15. A tunable center of pressure threshold, * cop , controlled the transition from S1 to S2. Then, after a tunable duration, t 2→3 , the FSM transitioned to S3. Next, a TO event triggered the transition to S4. Finally, knee extension (θ k < 0) caused a transition to S5, where the FSM remained until returning to S1 at HS. During transitions, the impedance parameters were rate-limited to prevent step changes in torque. In the FSMC, the torque command was given by (1), where K, B, and θ eq depended on the current FSM state (given in Table IV).
Many methods have been proposed for deciding when to switch between sets of impedance parameters for different tasks, including simple threshold methods [16] and more complex machine learning methods [83], [84]. We employed a strategy similar to [16] where the prosthesis directly estimated the ground incline using the method described in Section IV-B. Then, a secondary FSM was used to select between parameter sets based on the estimated inclineγ. To prevent rapid switching between parameters at the boundaries, overlap was included in the switching thresholds ( Fig. 15(b)).

C. Additional Detailed Results
Table IV lists the results from the impedance parameter tuning for the FSMC, including the tuned impedance parameters, transition thresholds, and tuning times. The stiffness K, damping B, and equilibrium angle θ eq were tuned by the research team for each of the five states (S1, S2, S3, S4, and S5) at three baseline tasks for each participant. Next, Table V shows the kinematic and kinetic RMSE values for each participant during both the steady state and continuously varying trials, each separated by stance and swing.  IV  RESULTS OF THE IMPEDANCE PARAMETER TUNING FOR THE FSMC FOR BOTH PARTICIPANTS P1 AND P2 AT EACH BASELINE TASK. THE FSMC CONSISTS OF FIVE  STATES: S1, S2, S3, S4, AND S5, EACH WITH UNIQUE PARAMETERS   TABLE V  KINEMATIC