3.1. Hybrid-dynamics engine

The hybrid-dynamics engine performs dynamics calculations in which a selection of joints is configured as a forward- and the others as an inverse-dynamics joint (Featherstone, Reference Featherstone2008). In case of the former, joint torque $ {\tau}_i $ is provided and the hybrid-dynamics engine calculates joint acceleration $ {\ddot{q}}_i $ , in which case joint position $ {q}_i $ and velocity $ {\dot{q}}_i $ are obtained through numerical integration. In case of the latter, $ {q}_i $ , $ {\dot{q}}_i $ and $ {\ddot{q}}_i $ are provided to calculate $ {\tau}_i $ . Joint configurations are as follows:

• • Joints 5 and 8 are that of the knee and ankle and controlled by the prosthesis, for which they are forward-dynamics joints ( $ {\tau}_5\hskip0.35em =\hskip0.35em {\tau}_{\mathrm{knee}} $ and $ {\tau}_8\hskip0.35em =\hskip0.35em {\tau}_{\mathrm{ankle}} $ ).

• • Joints 6 and 7 are introduced for the purpose of obtaining force and moment readings, for which they are configured as inverse-dynamics joints, and in fact locked, so that $ {q}_6\hskip0.35em =\hskip0.35em {\dot{q}}_6\hskip0.35em =\hskip0.35em {\ddot{q}}_6\hskip0.35em =\hskip0.35em 0 $ and $ \hskip2em {q}_7\hskip0.35em =\hskip0.35em {\dot{q}}_7\hskip0.35em =\hskip0.35em {\ddot{q}}_7\hskip0.35em =\hskip0.35em 0 $ .Footnote ¹ $ {\tau}_6 $ corresponds to vertical force and $ {\tau}_7 $ to the bending moment.

• • Joints 1 to 4 are configured differently depending on the scope of the simulation. If ground contact is involved (stance), at least two independent joints must be configured as forward dynamic joints to assure motion freedom, typically trunk translation (1 and 2). If also ground impact is involved, like in walking, it is advised to configure also the hip joint as forward dynamic joint in order to reduce impact force and slip. The case studies presented in this manuscript use this configuration.

3.2. Symbol conventions

For readability and consistency with literature, additional symbols are introduced to represent variables with units different from the SI system and signs that could be opposite from their definition in the dynamic model in Section 3.1. Angular positions and velocities $ \theta $ and $ \omega $ are expressed in $ {}^{\circ } $ and rpm, respectively. Corresponding rotational joints are defined so that a positive value corresponds to flexion: we have trunk angle $ {\theta}_{\mathrm{t}}\hskip0.35em =\hskip0.35em {q}_3 $ , knee angle $ {\theta}_{\mathrm{k}}\hskip0.35em =\hskip0.35em -{q}_5 $ and ankle angle $ {\theta}_{\mathrm{a}}\hskip0.35em =\hskip0.35em -{q}_8 $ . In addition, hip angle $ {\theta}_{\mathrm{h}} $ is expressed with relation to the absolute vertical, so that $ {\theta}_{\mathrm{h}}\hskip0.35em =\hskip0.35em {q}_3+{q}_4 $ . Forces $ F $ are expressed per body weight (N N⁻¹) and moments $ M $ and torques $ T $ per body mass (N m kg⁻¹). Furthermore, several sign changes apply: we have trunk horizontal force $ {F}_{\mathrm{x}}\hskip0.35em =\hskip0.35em {\tau}_1 $ , trunk vertical force $ {F}_{\mathrm{y}}\hskip0.35em =\hskip0.35em {\tau}_2 $ , trunk torque $ {T}_{\mathrm{t}}\hskip0.35em =\hskip0.35em {\tau}_3 $ , hip torque $ {T}_{\mathrm{h}}\hskip0.35em =\hskip0.35em {\tau}_4 $ , knee torque $ {T}_{\mathrm{k}}\hskip0.35em =\hskip0.35em {\tau}_5 $ , tube force $ {F}_{\mathrm{u}}\hskip0.35em =\hskip0.35em -{\tau}_6 $ , tube moment $ {M}_{\mathrm{u}}\hskip0.35em =\hskip0.35em -{\tau}_7 $ and ankle torque $ {T}_{\mathrm{a}}\hskip0.35em =\hskip0.35em -{\tau}_8 $ . Values for associated controller parameters are expressed in corresponding units.

3.3. Compliant ground model

Ground contact is modelled by a compliant ground model (Azad and Featherstone, Reference Azad and Featherstone2010), instead of by employing non-continuous impulsive dynamics and switching constraint functions. This approach allows for impact forces to remain realistically finite, and can be used to take an extent of compliance into account that is found in a real foot prosthesis, socket and hip joint, which is otherwise not modelled by the multi-rigid-body mechanism. The compliant ground model implements a non-linear spring damper element between defined contact points and the ground, so that the normal force $ {F}_{\mathrm{n},j} $ of a point $ j $ equals

$$ {F}_{\mathrm{n},j}\hskip0.35em =\hskip0.35em -\sqrt{\max \left(-{y}_j,0\right)}\left({Ky}_j+D{\dot{y}}_j\right), $$

where $ {y}_j $ is the absolute y-coordinate of point $ j $ , and $ y\hskip0.35em =\hskip0.35em 0 $ defines the ground plane. At least two contact points must be defined in order to exercise a moment: one in the heel and one in the ball of the foot, respectively, defined by the coordinates $ [{x}_{\mathrm{heel}},{-y}_{\mathrm{foot}}] $ and $ [{x}_{\mathrm{ball}},{-y}_{\mathrm{foot}}] $ in the coordinate frame of the foot, as shown in Figure 2. Kinematic parameter choices are listed in Table 1. Note that the distance between the heel and ball essentially represents the effective foot length, which is found to be between 0.63 to 0.81 times the total length of a prosthetic foot (Hansen et al., Reference Hansen, Sam and Childress2004). Stiffness and damping constants for the two contact points are set to $ K\hskip0.35em =\hskip0.35em 2\times {10}^5\hskip0.35em {\mathrm{Nm}}^{-1.5} $ and $ D\hskip0.35em =\hskip0.35em 2\times {10}^4\hskip0.35em {\mathrm{Nsm}}^{-1.5} $ , respectively. This stiffness implies that the foot sinks approximately 1.5 cm in the ground when the total body weight (ca. 750 N) is applied on the foot, which might suggest that the floor is rather soft. However, the stiffness also models vertical or longitudinal compliance found in a real leg and foot prosthesis, which is otherwise not taken account for by the rigid bodies and its rotational joints of the presented leg model. An increased stiffness leads to higher internal peak forces and torques at impact, but was anyway found to not significantly affect system dynamics: increasing it tenfold did not destabilize the system. Alternatively, a higher stiffness can be combined with the tube extension joint (joint 6) configured as a stiff forward dynamics joint, as mentioned in Section 3.1, so that it too acts as a spring. The damping value was tuned heuristically to achieve an approximate critical damping behavior for typical impacts that occur during walking. The model furthermore implements the ability to both stick and slip, with a friction coefficient of $ \mu \hskip0.35em =\hskip0.35em 0.5 $ . Corresponding tangential forces are calculated using the same spring-damper constants.

4.1. Trajectory generation

Trajectory generation is realized so that smooth trajectories are obtained for applicability as inverse dynamics references, while not being predetermined to allow for interactive change of user intent, ambulation velocity and varying step times. This is achieved by a two-step procedure: rough and smooth joint reference generation.

Rough joint references

Firstly, rough joint position reference templates are predefined for every user action, which can be discontinuous. These include a reference trajectory for standing still (i.e., a constant zero joint angle), an initiation step, a termination step, and a cyclic walking step, the latter of which normalized for step time, as shown in Figure 4. Through the FSM of the human controller, these templates are scaled in time, cut, and stitched together to define the rough joint reference $ {\boldsymbol{q}}_{\mathrm{rr}} $ , based on user intent (walking or standing still), change of intent (starting or stopping a walk), intended ambulation speed, and step cycle completion.

Figure 4. Predefined rough reference trajectories $ {\theta}_{\mathrm{trr}} $ and $ {\theta}_{\mathrm{hrr}} $ for the trunk and hip, respectively, for each of the four main states of the human controller (see Figure 3), normalized for the anticipated stride cycle duration $ {\tilde{t}}_{\mathrm{cyc}} $ (see x-axis). The resulting smooth reference trajectories $ {\theta}_{\mathrm{tr}} $ and $ {\theta}_{\mathrm{hr}} $ for this particular chain of reference trajectories are depicted as well, as obtained by passing the rough reference trajectories through the filter in equation 1. While walking, the trunk is sent into a forward lean of $ {\theta}_{\mathrm{trr}}\hskip0.35em =\hskip0.35em -4\deg $ , and the hip roughly follows a trajectory between $ {\theta}_{\mathrm{h},\max}\hskip0.35em =\hskip0.35em {25}^{\circ } $ and $ {\theta}_{\mathrm{h},\min}\hskip0.35em =\hskip0.35em -{10}^{\circ } $ .

If the human intent changes from stand to walk, the human controller will pass once through a transition state to initiate gait (ReStepStart), and vice versa (ReStepStop). In reality, gait initiation and termination can be executed by both legs and in different ways (van Keeken et al., Reference van Keeken, Vrieling, Hof and Postema2008). However, in the presented model, the initiation step is implemented such that the intact leg makes the first step (executed by the SLIP model). Forward thrust is generated with the residual leg by flexing the hip while standing. During walking, the reference position template of the next step is defined at the heel strike (step initiation reflex) of the residual leg (i.e., when transitioning to ReStepCycle or ReStepStop) according to the registered intent and intended ambulation speed at that time instance, at which point the reference trajectory of the previous step is terminated and that of the new step initiated. In particular, the reference trajectory is scaled in time as function of the intended ambulation speed. As such, no different trajectories are defined for different ambulation velocities, though such an implementation would be possible too. The scaling is done with an estimate of the duration of the next step cycle $ {\tilde{t}}_{\mathrm{cyc}} $ according to $ {\tilde{t}}_{\mathrm{cyc}}\hskip0.35em =\hskip0.35em {\tilde{x}}_{\mathrm{cyc}}/{v}_{\iota } $ , where the estimate of stride length $ {\tilde{x}}_{\mathrm{cyc}} $ is the stride length of the previous step cycle $ {x}_{\mathrm{cyc}} $ , measured as the distance between the heel of the residual leg at previous and current heel strike (initial estimate $ \tilde{x}\hskip0.35em =\hskip0.35em 1.33\;\mathrm{m} $ ).

Smooth joint references

Secondly, smooth position reference trajectories $ {\boldsymbol{q}}_{\mathrm{r}} $ , continuous velocity reference trajectories $ {\dot{\boldsymbol{q}}}_{\mathrm{r}} $ and discontinuous but finite acceleration reference trajectories $ {\ddot{\boldsymbol{q}}}_{\mathrm{r}} $ are obtained through a two-dimensional filter in order to make them eligible as inputs for the hybrid dynamics engine for joints configured as inverse-dynamics joints. Because of the filter, the output remains smooth even during state switches of the FSM, which is crucial as state changes are not timed and can occur unexpectedly due to, for example, the step initiation reflex.

The filter to obtain $ {\boldsymbol{q}}_{\mathrm{r}} $ and its derivatives is implemented as follows:

(1) $$ 4{\zeta}^2{\boldsymbol{q}}_{\mathrm{r}\mathrm{r}}\hskip0.35em =\hskip0.35em {\ddot{\boldsymbol{q}}}_{\mathrm{r}}{t}_{\mathrm{delay}}^2+4{\zeta}^2{\dot{\boldsymbol{q}}}_{\mathrm{r}}{t}_{\mathrm{delay}}+4{\zeta}^2{\boldsymbol{q}}_{\mathrm{r}} $$

with a damping ratio $ \zeta \hskip0.35em =\hskip0.35em 1/\sqrt{2} $ and a variable time delay $ {t}_{\mathrm{delay}}\hskip0.35em =\hskip0.35em 0.1{\tilde{t}}_{\mathrm{cyc}} $ (during standing, $ {t}_{\mathrm{delay}}\hskip0.35em =\hskip0.35em 0.2 $ ) so that the shapes of the reference trajectories $ {\boldsymbol{q}}_{\mathrm{r}} $ ( $ {\theta}_{\mathrm{t},\mathrm{r}}\hskip0.35em =\hskip0.35em {q}_{3\mathrm{r}} $ and $ {\theta}_{\mathrm{h},\mathrm{r}}\hskip0.35em =\hskip0.35em {q}_{3\mathrm{r}}+{q}_{4\mathrm{r}} $ for the trunk and hip) become invariant of the commanded ambulation velocity. The result of the filter output is also shown in Figure 4.

An additional advantage of the two-step procedure for trajectory generation is that it requires the definition of only compact and tangible or descriptive data structures that describe the approximate motion of the leg (i.e., those of the rough reference trajectories), rather than having to rely on long lookup tables that are difficult to customize. As shown in Figure 4, the hip trajectory for a cyclic walking step can be described as gradually moving backwards from a starting flexion angle $ {\theta}_{\mathrm{h},\max } $ to an extension angle $ {\theta}_{\mathrm{h},\min } $ during stance, after which it accelerates back toward its starting angle during the beginning of swing, while anticipating the next foot strike.

Joint commands

As mentioned in Section 3.1, both the trunk and hip rotation joints are eligible to be configured as inverse-dynamics joints without over-constraining the body kinematics during stance. This can be realized by directly taking the generated smooth reference trajectory $ {q}_{\mathrm{r}} $ and derivatives of the corresponding joint as inputs to the hybrid-dynamics engine. By doing so for the trunk, the system dynamics are equivalent for assuming the trunk to be a point mass, which has previously shown to be a fair approach in simple dynamic walking models (McGeer, Reference McGeer1990; Geyer et al., Reference Geyer, Seyfarth and Blickhan2006; Blickhan et al., Reference Blickhan, Seyfarth, Geyer, Grimmer, Wagner and Günther2007). When the trunk is a point mass, it is implied that—as long as no joint torque limits are imposed on the trunk joint, and as long as hip rotation directly compensates for trunk rotation (which it does, because hip rotation trajectories are defined with respect to the absolute vertical, i.e., $ {\theta}_{\mathrm{h}}\hskip0.35em =\hskip0.35em {q}_3+{q}_4 $ )—the trunk orientation does not affect the outcome of the dynamics. Nevertheless, for purpose of demonstration and aesthetics, the trunk is sent in a forward lean of $ {4}^{\circ } $ while walking, as shown in Figure 4.

The hip joint follows a more dynamic trajectory. Whereas it too can be configured as inverse-dynamics joint, it is instead configured as a forward-dynamics joint to follow its reference trajectory with a stiff joint-impedance controller, in order to moderate impact forces and to eliminate occurrences of (micro-)slip throughout stance:

$$ {\tau}_4\hskip0.35em =\hskip0.35em {T}_{\mathrm{h}}\hskip0.35em =\hskip0.35em {P}_{\mathrm{h}}\left({q}_4-{q}_{4\mathrm{r}}\right)+{D}_{\mathrm{h}}\left({\dot{q}}_4-{\dot{q}}_{4\mathrm{r}}\right), $$

with a manually tuned stiffness gain $ {P}_{\mathrm{h}}\hskip0.35em =\hskip0.35em 2\;\mathrm{kN}\;\mathrm{m}\hskip0.35em {\mathrm{rad}}^{-1} $ while walking and $ {P}_{\mathrm{h}}\hskip0.35em =\hskip0.35em 500\;\mathrm{N}\;\mathrm{m}\hskip0.35em {\mathrm{rad}}^{-1} $ while standing, and damping gain $ {D}_{\mathrm{h}}\hskip0.35em =\hskip0.35em 100\;\mathrm{N}\;\mathrm{m}\;\mathrm{s}\hskip0.35em {\mathrm{rad}}^{-1} $ . The stiffness is manually tuned to be maximal for eliminating visible foot slip during stance.

4.2. Walking: SLIP trunk forces

During stance of the intact leg (HeStance and DoSupport), the trunk is subject to a force from the SLIP model $ {F}_{\mathrm{p}} $ (as displayed also in Figure 2), a forward thrust $ {F}_{\mathrm{thrust}} $ and optionally a vertical lift force $ {F}_{\mathrm{lift}} $ :

The resulting applied trunk forces are calculated as follows:

$$ {F}_{\mathrm{x}}\hskip0.35em =\hskip0.35em {\tau}_1\hskip0.35em =\hskip0.35em {F}_{\mathrm{p}}\frac{x_{\mathrm{p}}}{l_{\mathrm{p}}}+{F}_{\mathrm{thrust}}\hskip5em {F}_{\mathrm{y}}\hskip0.35em =\hskip0.35em {\tau}_2\hskip0.35em =\hskip0.35em {F}_{\mathrm{p}}\frac{y_{\mathrm{p}}}{l_{\mathrm{p}}}+{F}_{\mathrm{lift}}, $$

in which $ {l}_{\mathrm{p}}\hskip0.35em =\hskip0.35em \sqrt{x_{\mathrm{p}}^2+{y}_{\mathrm{p}}^2} $ , $ {x}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {q}_1-{x}_{\mathrm{p},0} $ and $ {y}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {q}_2 $ are the pendulum’s actual length and its corresponding x- and y-components, respectively. $ {x}_{\mathrm{p},0} $ is the pendulum’s floor coordinate. $ {F}_{\mathrm{p}} $ is calculated as $ {F}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {k}_{\mathrm{p}}\left({l}_{\mathrm{p}}-{l}_{\mathrm{p},0}\right)+{c}_{\mathrm{p}}{\dot{l}}_{\mathrm{p}} $ , in which the parameters $ {k}_{\mathrm{p}} $ and $ {c}_{\mathrm{p}} $ are linear stiffness and damping coefficients, respectively, and $ {l}_{\mathrm{p},0} $ the pendulum’s rest length.

Parameters are chosen so as to reproduce the double ‘camel back’ bounce that is typically observed in the ground reaction force of healthy human walking. Damping $ {c}_{\mathrm{p}} $ has been added in order to account for natural losses at ground impact and to modulate an excess of bouncy behavior. Accordingly, $ {k}_{\mathrm{p}} $ and $ {c}_{\mathrm{p}} $ are calculated as follows:

$$ {k}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {m}_{\mathrm{p}}{\left(2\pi /{\tilde{t}}_{\mathrm{p}}{\sqrt{1-{\zeta}^2}}^{-1}\right)}^2\hskip6em {c}_{\mathrm{p}}\hskip0.35em =\hskip0.35em 2\zeta \sqrt{k_{\mathrm{p}}{m}_{\mathrm{p}}} $$

with the estimated pendulum ground contact duration $ {\tilde{t}}_{\mathrm{p}} $ , an effective pendulum inertia $ {m}_{\mathrm{p}} $ and a fixed damping ratio of $ \zeta \hskip0.35em =\hskip0.35em 0.2 $ . Per estimate, $ {m}_{\mathrm{p}} $ is chosen equal to the trunk mass $ {m}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {m}_3 $ , and $ {\tilde{t}}_{\mathrm{p}} $ is estimated to be half that of $ {\tilde{t}}_{\mathrm{cyc}} $ .

The pendulum’s rest length $ {l}_{\mathrm{p},0} $ and floor coordinate $ {x}_{\mathrm{p},0} $ are chosen with an aim to produce a symmetric gait. For this purpose, $ {l}_{\mathrm{p},0} $ is chosen equal to the total length of the residual leg: $ {l}_{\mathrm{p},0}\hskip0.35em =\hskip0.35em {l}_{\mathrm{leg}} $ . $ {x}_{\mathrm{p},0} $ is set at the transition to HeStance to achieve rest length at impact: $ {x}_{\mathrm{p},0}\hskip0.35em =\hskip0.35em {q}_1+\sqrt{l_{\mathrm{p},0}^2-{y}_{\mathrm{p}}^2} $ . This transition is triggered when $ {y}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {q}_2 $ falls below a threshold height $ {y}_{\mathrm{p},0} $ , which is chosen to correspond to the anticipated trunk height at heel strike of the residual leg, see Figure 2. The SLIP’s corresponding angle of attack can be calculated as $ {\alpha}_{\mathrm{p}}\hskip0.35em =\hskip0.35em {\mathit{\cos}}^{-1}\left({y}_{\mathrm{p},0}/{l}_{\mathrm{p},0}\right)\hskip0.35em =\hskip0.35em {20.6}^{\circ } $ .

4.3. Walking: additional trunk support forces

In addition to forces resulting from the SLIP model, support forces can be applied to the trunk, which can have different functions. The application of any such force is optional. However, it is necessary to include a forward force to enable a gait that is energetically symmetric, since the SLIP model itself is purely dissipative, for which it can only generate asymmetric (though plausible) gaits without the injection of positive power.

In this work, we introduce a forward thrust force, which for purpose of simplicity is constant throughout stance of the intact leg (i.e., during activation of the SLIP model), but in addition, it is adjusted from step to step to regulate for the intended ambulation velocity $ {v}_{\iota } $ , as further explained below. Note that whereas this force is constant, the total horizontal force on the trunk is not, because it includes also the horizontal force components that is imposed by the SLIP model. It is not claimed that the resulting force profile is exactly similar to that exerted by a human, but it imposes possible and even realistic (due to the balancing of power injection) dynamics on the lower limb prosthesis, which is what matters for its functionality and feedback control.

In addition to the forward thrust force, other forces can be applied to model any type of realistic disturbance, caused by the human or environment, which could either help or perturb the prosthesis system and its controller. One example of common user behavior is lifting the hip of the residual leg to create additional ground clearance, enforcing an asymmetric gait. Especially for low ambulation velocities, it was found that the corresponding low inertial forces do not allow for sufficient knee flexion to prevent scuffing or tripping when using a passive knee and ankle prosthesis. Since the SLIP model parameters are defined to realize a symmetric gait, a vertical support force can be added to create additional hip lift. Note that whereas such asymmetric user behaviors should not be promoted, the prosthesis should also work for them. The possible implementation of a conceptual lifting force to enable walking with lower ambulation velocities is further explained below.

Also other forces or impulses can be applied, such as those resulting from a push, gust of wind or uneven terrain. These are not further investigated in this work, but are valid sources as disturbances for improving the robustness of a controller.

Horizontal thrust force for power injection and speed management

A support forward thrust force is applied to inject positive power and to regulate the intended ambulation velocity $ {v}_{\iota } $ . The applied thrust force is constant throughout the stance phase of the intact leg and applied using the quadratic relation $ {F}_{\mathrm{thrust}}\hskip0.35em =\hskip0.35em {K}_{\mathrm{thrust}}{v}_{\iota}^2 $ with gain $ {K}_{\mathrm{thrust}} $ . This relation is used because the losses it aims to compensate for are strongly positively correlated with velocity, typically approximately quadratic or of higher order, depending on many factors of the non-linear system, including the employed controller, internal friction and impact conditions. Since $ {K}_{\mathrm{thrust}} $ is not algebraically derivable and poorly predictable, it is updated once every step cycle at initiation of the stance phase of the residual leg, with the aim to converge to the desired ambulation velocity $ {v}_{\iota } $ , requiring only an initial estimate (chosen as $ {K}_{\mathrm{thrust}}\hskip0.35em =\hskip0.35em 50 $ ). The update is based on an energy balance that compares the difference of the translational kinetic energy of the system for the actual and target ambulation velocity $ {v}_{\mathrm{target}} $ and a prediction of work done by the thrust force given an estimation of stride length:

(2) $$ {K}_{\mathrm{thrust}}\leftarrow {K}_{\mathrm{thrust}}+\frac{\Delta {E}_{\mathrm{kx}}}{\mathrm{1}/\mathrm{2}{\tilde{x}}_{\mathrm{cyc}}{v}_{\mathrm{target}}^2}, $$

where $ \Delta {E}_{\mathrm{kx}}\hskip0.35em =\hskip0.35em \mathrm{1}/\mathrm{2}m\left({v}_{\mathrm{target}}^2-{v}_{\mathrm{x}}^2\right) $ denotes the difference of translational kinetic energy for the measured ambulation velocity $ {v}_{\mathrm{x}} $ and target velocity $ {v}_{\mathrm{target}} $ . $ {v}_{\mathrm{target}} $ is defined as the average of the intended ambulation velocity $ {v}_{\iota } $ registered at initiation of the current and previous step cycle of the residual leg, as an implementation of a finite impulse response filter to prevent updates of $ {K}_{\mathrm{thrust}} $ to overshoot. $ {v}_{\mathrm{x}} $ is calculated as the average velocity of the previous step cycle of the residual leg: $ {v}_{\mathrm{x}}\hskip0.35em =\hskip0.35em {x}_{\mathrm{cyc}}/{t}_{\mathrm{cyc}} $ , with $ {t}_{\mathrm{cyc}} $ the duration between the last two heel strikes of the residual leg.

In addition, when the target velocity changes, an additional thrust force is added to accelerate or brake the system according to a prediction of the work required to realize a corresponding change of translational kinetic energy, employing a similar energy balance as that used for updating $ {K}_{\mathrm{thrust}} $ as in equation 2.

Vertical lift to emulate an asymmetric gait and prevent foot scuff

A support vertical lift force can be applied to realize additional lift of the hip during stance of the residual leg for lower ambulation velocities, to prevent tripping with a passive prosthesis. In the presented model, this vertical support force $ {F}_{\mathrm{lift}} $ is synchronized with the SLIP rotation and defined as follows:

$$ {F}_{\mathrm{lift}}\hskip0.35em =\hskip0.35em \left\{\begin{array}{ll}0& \mathrm{if}\;{v}_{\mathrm{target}}\hskip0.35em \ge \hskip0.35em 1.1\;{\mathrm{ms}}^{-1}\\ {}0.1\;{k}_{\mathrm{p}}\left(1.1-{v}_{\mathrm{target}}\right)\cos \left(\frac{\pi }{2}\frac{\theta_{\mathrm{p}}}{\alpha_{\mathrm{p}}}\right)& \mathrm{if}\;{v}_{\mathrm{target}}<1.1\;{\mathrm{ms}}^{-1}.\end{array}\right. $$

This implementation lifts the leg up to approximately an additional 1 cm for every 0.1 ms $ {}^{-1} $ that the target velocity is lower than 1.1 ms $ {}^{-1} $ , which was heuristically found to avoid scuffing issues for ambulation speeds of at least 0.8 ms $ {}^{-1} $ for the reference passive prosthetic leg that is showcased in Section 5. Additional lift is added during gait termination to prevent tripping. Naturally, depending on the goal of the simulation study, the implementation of an intentional asymmetric gait through a lift force can be increased, decreased, changed or removed entirely.

5.1. Tuning procedure

Following a preliminary implementation of the prosthesis device and controller, initial parameter estimates can mostly be interpreted from reference material. Tuning of parameters or further refinements of the control system are done directly in forward dynamic simulations, which in our work are performed in Matlab Simulink (version 2020b). For new designs, the primary objective of a tuning procedure is to realize a feasible and stable gait with the closed-loop system. Secondarily, the objective is to improve (bio)mechanical characteristics for that gait, such as reducing internal peak forces, or to improve its robustness to deviations from the gait, such as different and changing ambulation speeds, starting and stopping behavior, but also other possible sources of disturbances and changes of human actions.

Generally, the holistic design and tuning of a multidimensional and intrinsically unstable forward dynamic system is a daunting task, for which typically is resorted to automated optimization techniques rather than manual heuristic processes. Especially for more complex systems like OpenSim, this is a sensible approach. However, disadvantages of such holistic optimization problems are that they are difficult to define and time consuming, and they provide lack of insight into the system. An initial heuristic approach can be favorable to set the scope of design solutions and possible future optimizations, for as long as the system is comprehensible. We found that a major advantage of the presented simplistic framework is that achieving a feasible solution (stabilization) and basic improvements can be achieved heuristically with a comprehensible and straightforward tuning approach. In addition, a simulation lasts on the order of only seconds,Footnote ² allowing for fast progression.

Whereas tuning remains an iterative process, the general approach is to move from the device to the controller, and from simple actions (slow walk) to extended behaviors (change of speed, starting and stopping, disturbances, etc.). In particular, the approach is to first obtain a feasible result for a stable unperturbed walking gait with a slow target ambulation speed, yet one that is fast enough to allow for passive knee flexion as a result of inertial forces. For this purpose, we found $ {v}_{\iota}\hskip0.35em =\hskip0.35em 1.1\;{\mathrm{ms}}^{-1} $ to be a good target velocity. Initial conditions of state variables for a cyclic (rather than starting or stopping) step are chosen as $ {\theta}_{\mathrm{k}}\hskip0.35em =\hskip0.35em {6}^{\circ } $ , $ {\theta}_{\mathrm{a}}\hskip0.35em =\hskip0.35em {0}^{\circ } $ , the hip and trunk angles equal to their initial angle of their rough references for a walk cycle (see Figure 4, that is, $ {\theta}_{\mathrm{t}}\hskip0.35em =\hskip0.35em -{4}^{\circ } $ and $ {\theta}_{\mathrm{h}}\hskip0.35em =\hskip0.35em {25}^{\circ } $ ), the trunk height so that ground contact with the heel is realized, and the trunk forward velocity so that it is 20% higher than $ {v}_{\iota } $ (a higher than average velocity prior to impact of the heel strike of the residual leg). Other initial state variables are zero.

An overview of the proposed tuning procedure is visualized in Figure 5. To reduce system complexity, it is advisable to firstly assume the use of a passive ankle prosthesis, like in the presented case studies. For the knee prosthesis, the procedure can then be summarized as follows:

1. 1. Remove the actuator(s) from the joint,Footnote ³ leaving only the physical end stops but an otherwise transparent knee joint, and observe limit cycle instability for cyclic walking at $ {v}_{\iota}\hskip0.35em =\hskip0.35em 1.1\;{\mathrm{ms}}^{-1} $ . The overly transparent knee joint typically results in a fall after several steps, but not immediately.

2. 2. Re-introduce the actuator (friction, inertia, springs), and observe if cyclic behavior has improved or deteriorated. A minimum amount of friction (energy dissipation) generally improves the behavior, leading to a later fall or even a stable limit cycle.Footnote ⁴ If instead behavior deteriorates, likely sources of friction and inertia are too high, leading to scuffing or tripping during swing as a result of insufficient knee flexion.

3. 3. If the behavior has deteriorated, either directly attempt the implementation of a feasible controller that can counter some of these actuator or device deficiencies (e.g., friction, inertia) (step 5), such as a transparency, minimum jerk trajectory or impedance controller if the prosthesis is powered, or firstly relax device deficiencies until walking behavior is improved. Depending on the scope of the design study, these relaxations can be the foundation of an actual design iteration if the mechanical design has not yet been decided on. Otherwise, these changes can be temporal, and be reversed in a subsequent tuning stage when the controller is introduced.

4. 4. If a high level walking controller is employed, such as an FSM, verify that it functions correctly, despite the lack of controller outputs. In case of an FSM, improve transitions where needed by analyzing sensor data, until it cycles through walking states according to expectations, without error.

5. 5. Search to improve walking behavior by including and tuning controller outputs, and possibly by changing device parameters if the stage of design permits this (e.g., changing a spring constant). If device deficiencies had previously been relaxed, they can be reintroduced and compensated for by controller actions. (Bio)mechanical features can be used and pursued as reference for improving walking behavior (Grimmer and Seyfarth, Reference Grimmer and Seyfarth2014), such as:

   • • A maximum knee flexion angle of approximately $ {60}^{\circ } $ during swing. A significantly lower angle can lead to insufficient ground clearance, but an angle that is significantly higher might cause the knee to extend insufficiently when entering stance. Ground clearance itself can be monitored as y-coordinates of foot contact points, to understand if the model is (close to) scuffing or tripping.

   • • Limited knee flexion during stance. In particular, at ground impact and lift-off, the knee is ideally nearly extended, and during mid-stance the knee flexes maximally approximately $ {15}^{\circ } $ . Additional flexion increases the likelihood of buckling, especially when braking or terminating gait. Additional stance flexion resistance could be introduced to lower the maximum knee angle.

   • • Reduction of internal peak forces as a result of impacts with physical end stops. For example, (angle-based) leg resistance during swing extension can be increased to lower impact forces with the knee extension end stop.

   • • Reduction of prosthesis energy consumption and motor saturation. This is of particular interest for powered prostheses.

6. 6. Once a satisfactory walking behavior is established for an ambulation speed of $ {v}_{\iota}\hskip0.35em =\hskip0.35em 1.1\;{\mathrm{ms}}^{-1} $ , the system can be challenged and tuned for different and changing intended ambulation speeds and gait initiation and termination, which can be seen as disturbances to the control system. Especially braking challenges the controller during stance, as it results in additional flexion load on the knee that could result in the system buckling to the ground. At this stage, also high level intent recognition and speed adaptation control modules can be experimented with.

The successful completion of these steps results in a prosthesis and controller system that is feasible for a certain repertoire of possible human behaviors, which formulate a strong basis for design realization. However, it should be noted that the converged solution is likely not unique, nor is it claimed to be optimal in any sense, especially because the human system is simplified, and not exercising optimal actions itself.

Figure 5. A proposed tuning procedure: whereas one could resort to optimizations to holistically tune the dynamic human-prosthesis system, the relative simplicity of the system also allows for manual tuning.

If the design project also includes the development of an active ankle prosthesis, one could proceed to its introduction and repeat the tuning process. In addition, it is certainly possible to repeat the tuning process for a broader range of disturbances and human actions (e.g., change or removal of support forces), trajectories (e.g., different hip reference trajectory) and parameters (e.g., different body sizing), in order to improve controller robustness for a broader repertoire of possible human behaviors. However, it should be noted that the effectiveness of manual tuning might crawl to a slow when increasing the dimensionality of human deviations, and that it might instead be more productive to implement an optimization procedure at this stage.

5.2. Variable damping commercial leg

The implementation of the commercial leg is inspired from a simplification and interpretation of material presented in a patent of Össur (Bisbee III et al., Reference Bisbee, Elliott and Oddsson2016). Its FSM and state-dependent damping formulae are depicted in Figure 6. Parameter choices of the FSM follow largely from direct interpretations, and are partly tuned according to the procedure described in Section 5.1, as further detailed in the subsequent Section 5.2.2.

Figure 6. Prosthesis FSM and state-based variable damping control for walking (states 1 to 5) and standing (states 1 and 2), largely inspired from Bisbee III et al. (Reference Bisbee, Elliott and Oddsson2016), implemented for proofing the simulation framework. Dashed transitions incorporate standing, non-standard walking and corrective behavior, which are not triggered in this demonstration.

The controller runs at a frequency of 200 Hz,Footnote ⁵ and the encoder is quantized for a resolution of $ {2}^{12} $ ticks per revolution, as realized by the prosthesis sensors module. The controller commands the quadratic damping coefficient $ c $ , which for purposes of demonstration is assumed to be a perfect damping source, so that the applied torque by the damper on the knee equals $ {T}_{\mathrm{D}}\hskip0.35em =\hskip0.35em -c{\dot{q}}_5\mid {\dot{q}}_5\mid $ , as implemented by the prosthesis device module. This module also implements sources of friction to define the knee torque, and implements a passive ankle prosthesis to define the ankle torque, as further detailed below. Kinematic and dynamic parameters of the prosthesis are according to those listed in Table 1.

5.2.1. Prosthesis device

In addition to the damper torque $ {T}_{\mathrm{D}} $ , a realistic amount of joint viscous friction, stiction, soft parallel extension spring and an extension end stop has been added to the prosthesis device model, so that

$$ {T}_{\mathrm{k}}\hskip0.35em =\hskip0.35em {\tau}_5\hskip0.35em =\hskip0.35em {T}_{\mathrm{D}}-{k}_{\mathrm{k}}{q}_5-{c}_{\mathrm{k}}{\dot{q}}_5-{T}_{\mathrm{stic}}\left({\dot{q}}_5\right)-{T}_{\mathrm{extlim}}\left({q}_5,{\dot{q}}_5\right)-{T}_{\mathrm{flxlim}}\left({q}_5,{\dot{q}}_5\right), $$

with stiffness $ {k}_{\mathrm{k}}\hskip0.35em =\hskip0.35em 1\;\mathrm{Nm}/\mathrm{rad} $ , damping $ {c}_{\mathrm{k}}\hskip0.35em =\hskip0.35em 0.2\;\mathrm{Nms}/\mathrm{rad} $ , stiction $ {T}_{\mathrm{stic}}\hskip0.35em =\hskip0.35em 0.1\;\tanh \left(100{\dot{q}}_5\right) $ and extension and flexion end stop torques $ {T}_{\mathrm{extlim}} $ and $ {T}_{\mathrm{flxlim}} $ , respectively.

To accurately incorporate the behavior of an end stop of finite thickness, it is implemented using an explosively progressive non-linear compression spring element. For this purpose, the one-dimensional Neo-Hookean element $ {f}_{\mathrm{neo}}\hskip0.35em =\hskip0.35em \lambda -{\lambda}^{-2} $ is employed, where $ {f}_{\mathrm{neo}} $ is proportional to the output force and $ \lambda $ the extension ratio of the compressed element, which accurately describes the behavior of materials with a Poisson ratio near 0.5, such as rubbers (Treloar, Reference Treloar1975). By including also a damping term, the extension limit torque $ {T}_{\mathrm{extlim}} $ is computed as follows:

$$ {T}_{\mathrm{extlim}}\hskip0.35em =\hskip0.35em \left\{\begin{array}{ll}\frac{2}{7}{T}_{\mathrm{half}}\left({f}_{\mathrm{neo}}\left(\lambda \right)+{t}_{\mathrm{c}}\dot{\lambda}\hskip0.1em \frac{\partial {f}_{\mathrm{neo}}}{\partial \lambda}\right)& \mathrm{if}\hskip0.35em {\theta}_{\mathrm{k}}<{\theta}_{\mathrm{extlim}}\\ {}0& \mathrm{if}\hskip0.35em {\theta}_{\mathrm{k}}\hskip0.35em \ge \hskip0.35em {\theta}_{\mathrm{extlim}}\end{array}\right. $$

with the extension ratio $ \lambda \hskip0.35em =\hskip0.35em \max \left({q}_5+{\theta}_{\mathrm{extlim}},0\right)/d+1 $ , starting angle $ {\theta}_{\mathrm{extlim}}\hskip0.35em =\hskip0.35em {3}^{\circ } $ , thickness $ d\hskip0.35em =\hskip0.35em {2}^{\circ }, $ torque at half compression (i.e., $ \lambda \hskip0.35em =\hskip0.35em 0.5 $ , or $ {1}^{\circ } $ of compression) $ {T}_{\mathrm{half}}\hskip0.35em =\hskip0.35em 100\hskip0.35em \mathrm{Nm} $ and damping time constant $ {t}_{\mathrm{c}}\hskip0.35em =\hskip0.35em 0.01\;\mathrm{s} $ . An identical end stop has been implemented for flexion, which operates in the opposite direction at a starting angle of $ {\theta}_{\mathrm{flxlim}}\hskip0.35em =\hskip0.35em {115}^{\circ } $ , producing $ {T}_{\mathrm{flxlim}} $ .

The prosthesis system also includes the ankle, which is modelled to be passive, and as such does not require a controller. It is implemented as a linear spring-damper system: $ {\tau}_8\hskip0.35em =\hskip0.35em {k}_{\mathrm{a}}\left({q}_{8,0}-{q}_8\right)-{c}_{\mathrm{a}}{\dot{q}}_8 $ , with rest angle $ {q}_{8,0}\hskip0.35em =\hskip0.35em {4}^{\circ } $ , stiffness $ {k}_{\mathrm{a}}\hskip0.35em =\hskip0.35em 6\hskip0.35em \mathrm{N}\hskip0.35em \mathrm{m}\hskip0.35em {\mathrm{rad}}^{-1} $ and damping $ {c}_{\mathrm{a}}\hskip0.35em =\hskip0.35em 1\hskip0.35em \mathrm{N}\hskip0.35em \mathrm{m}\hskip0.35em \mathrm{s}\hskip0.35em {\mathrm{rad}}^{-1} $ .

5.2.2. Results

Intermediate results

The procedure described in Section 5.1 has been used as a general guide to tune a number of controller parameters. Figure 7 shows several of the intermediate results of this process. The top graphs depict results for walking with $ {v}_{\iota}\hskip0.35em =\hskip0.35em 1.1\hskip0.35em {\mathrm{ms}}^{-1} $ and setting all joint friction and stiffness values and mid-layer controller parameters to zero, leaving only the end stops to generate a knee torque, so that $ {\tau}_5\hskip0.35em =\hskip0.35em -{T}_{\mathrm{extlim}}\left({q}_5,{\dot{q}}_5\right)-{T}_{\mathrm{flxlim}}\left({q}_5,{\dot{q}}_5\right) $ . The system ends up falling due to an overly transparent knee joint, but still manages to complete several steps.

Figure 7. Intermediate results of the tuning procedure, showing angular trajectories of the knee $ {\theta}_{\mathrm{k}} $ and hip $ {\theta}_{\mathrm{h}} $ (including also the hip rough and smooth reference trajectories $ {\theta}_{\mathrm{hrr}} $ and $ {\theta}_{\mathrm{hr}} $ , respectively), and torque values of the knee $ {T}_{\mathrm{k}} $ , damper $ {T}_{\mathrm{D}} $ and tube sensor $ {M}_{\mathrm{u}} $ . The top graphs show instability as a result of removing the actuator (both its friction and spring and controller output). By introducing the actuator without any controlled output, a limit cycle gait can already be achieved (step 2), but it is subject to internal impacts and sensitive to disturbance. The introduction of angle-based damping (soft end stops) rids these internal impacts and improves stability of the gait cycle (step 3).

The bottom left graphs show results for reintroduction of basic joint damping and stiffness values (i.e., $ {\tau}_5\hskip0.35em =\hskip0.35em -{k}_{\mathrm{k}}{q}_5-{c}_{\mathrm{k}}{\dot{q}}_5-{T}_{\mathrm{stic}}\left({\dot{q}}_5\right)-{T}_{\mathrm{extlim}}\left({q}_5,{\dot{q}}_5\right)-{T}_{\mathrm{flxlim}}\left({q}_5,{\dot{q}}_5\right) $ ). The introduction of minimal damping caused the system to no longer fall, but neither did it converge to a single-step limit cycle, possibly due to the lack of angle-based damping (soft end stop) in swing flexion. In addition, the knee joint collided with the physical extension end stop at the end of swing extension due to the lack of damping, resulting in high impacts, which are visible in the bottom left figure. As regards for the high level controller, we found that the FSM—directly interpreted from Bisbee III et al. (Reference Bisbee, Elliott and Oddsson2016)—correctly identified transitions, except for occasional preliminary transitions to late stance (pre-swing). This was caused because of a different moment reading. Where Bisbee III et al. (Reference Bisbee, Elliott and Oddsson2016) uses a moment measure closer to the knee joint to detect a transition to late stance with a threshold of 0.2 N m kg $ {}^{-1} $ , we employ a moment measure closer to the ankle (i.e., the tube sensor), which is higher than 0.2 N m kg $ {}^{-1} $ during most of stance. To resolve this, we have increased the threshold to 0.5 N m kg $ {}^{-1} $ , as listed in Figure 6. The framework has proven exceptionally useful to understand sensor data that are customary to specific design solutions, of which little or no data are available in the literature.

The bottom right graphs in Figure 7 show results after introduction of angle-based mid-layer damping commands in swing flexion and swing extension (soft end stops), as listed in Figure 6, resulting in a more stable gait and no more knee joint impacts with the physical knee extension end stop. However, due to the lack of stance support, the model was still found to exhibit high knee stance flexion angles, and to fall down when lowering the ambulation speed. This was resolved by introducing and tuning additional damping commands in stance, as listed in Figure 6.

Main results

All final controller parameter values are depicted in Figure 6. Simulations have demonstrated stable walking for a range of ambulation velocities of at least $ 0.8<{v}_{\iota }<1.4\hskip0.35em {\mathrm{ms}}^{-1} $ . Figure 8 shows results from a 40 s trial, of which various joint angles, forces and torques are plotted for the first 8 s. During the 40 s trial, the system starts from a vertical standing orientation, is commanded to start walking ( $ \iota \hskip0.35em =\hskip0.35em 1 $ ) at $ t\hskip0.35em =\hskip0.35em 1\hskip0.35em s $ with an ambulation velocity of $ {v}_{\iota}\hskip0.35em =\hskip0.35em 1.1\hskip0.35em {\mathrm{ms}}^{-1} $ , to accelerate and decelerate with both step and ramp commands, up to $ {v}_{\iota }<1.4\hskip0.35em {\mathrm{ms}}^{-1} $ and down to $ {v}_{\iota }<0.8\hskip0.35em {\mathrm{ms}}^{-1} $ , and to stop walking ( $ \iota \hskip0.35em =\hskip0.35em 0 $ ) at $ t\hskip0.35em =\hskip0.35em 40\hskip0.35em \mathrm{s} $ .

Figure 8. Simulation results of a walking trial, starting from and ending in standstill. The first plot shows animation frames, sampled at 10 Hz, with the SLIP drawn in red, and where the hip, toe and heel trajectories are drawn to demonstrate gait symmetry and foot clearance. The second plot shows the ambulation velocity $ {v}_{\mathrm{cyc}} $ (updated at heel strike of both the intact and residual leg) of the full simulation, which shows that it correctly tracks the intended velocity $ {v}_{\iota } $ with a delay of several steps. In this graph, transparent lines correspond to results of the simulation for complete removal of the forward thrust force, leading to an asymmetric “limping” gait and poorly tracked intended ambulation speed. The last three plots show various joint positions (global hip angle $ {\theta}_{\mathrm{h}} $ , its rough and smooth reference trajectories $ {\theta}_{\mathrm{hrr}} $ and $ {\theta}_{\mathrm{hr}} $ , knee angle $ {\theta}_{\mathrm{k}} $ , and ankle angle $ {\theta}_{\mathrm{a}} $ ), forces (horizontal $ {F}_{\mathrm{x}} $ and vertical trunk force $ {F}_{\mathrm{y}} $ and tube force $ {F}_{\mathrm{u}} $ ), and torques (knee torque $ {T}_{\mathrm{k}} $ , of which torque contributed by the damper $ {T}_{\mathrm{D}} $ , tube moment $ {M}_{\mathrm{u}} $ and ankle torque $ {T}_{\mathrm{a}} $ ), respectively, for the first 8 s of the simulation.

The results plotted in Figure 8 show high resemblance with human data recordings found in literature. The knee trajectory shows a healthy stance bounce and peak flexion values similar to those found in healthy user data, a comparison of which is shown in Figure 8 with data from Grimmer and Seyfarth (Reference Grimmer and Seyfarth2014). The ankle trajectory instead is very similar specifically to that of users wearing a passive ankle prosthesis, a comparison of which too is shown in Figure 8 with data from Grimmer and Seyfarth (Reference Grimmer and Seyfarth2014). In addition to kinematic data, also kinetic data follow expected patterns. Notable is the observed ‘camel back’ shape of the vertical force. This shape is also present in the tube force $ {F}_{\mathrm{u}} $ , but it appears with additional bounces since the feet have been modelled with only two ground contact points. In addition, it is sensitive to peaking at impact due to the high rigidity of the rigid-body system. Also peak values of torque measures correspond well to those found in literature, namely ca. 1 N m kg $ {}^{-1} $ and 1.5 N m kg $ {}^{-1} $ for the knee and ankle torque, respectively, at an ambulation velocity of $ 1.1\hskip0.35em {\mathrm{ms}}^{-1} $ .

Symmetry

The SLIP model is observed to successfully generate a symmetric gait regarding hip height, except for low ambulation velocities ( $ <1\hskip0.35em {\mathrm{ms}}^{-1} $ ), for which we have intentionally introduced an asymmetric lift force (Section 4.3) to generate sufficient ground clearance. Nevertheless, the asymmetry of the hip motion is minimal and barely visible in the animation plot of Figure 8. While walking with an ambulation speed of $ 1.1\hskip0.35em {\mathrm{ms}}^{-1} $ (i.e., when the lift force is not yet active), the hip peaks at a height of 902 mm and 900 mm from the ground for the stance phase of the residual and intact limb, respectively. This height corresponds with the leg length $ \left({l}_{\mathrm{leg}}\hskip0.35em =\hskip0.35em 920\hskip0.35em \mathrm{mm}\right) $ in combination with ground softness and minimal knee flexion. In the last step prior to terminating the walking gait, the hip peaks at a height of 911 mm during stance of the residual leg as a result of the lift force, which prevents it from falling. If instead the lift force is entirely removed from the simulation, the model is observed to start scuffing around $ t\hskip0.35em =\hskip0.35em 32\hskip0.35em \mathrm{s} $ , when decelerating to an ambulation speed of approximately $ 1\hskip0.35em \mathrm{m}{\mathrm{s}}^{-1} $ , and falls two steps later, before finishing the trial.

Regarding symmetry of forward velocity, it can be observed that $ {\dot{q}}_1 $ oscillates slightly more during stance of the intact (SLIP) limb than during that of the residual limb. Such an asymmetry is to be expected given the different modeling approach for the residual and intact leg: one having inertia, the other being a mass-less SLIP. However, the average stance phase (ignoring initiation and termination steps) of the residual and intact leg last 53% and 57% of the stride cycle, respectively (i.e., the double support phases amount to 10% of the stride cycle). This is only a small difference, which falls within the range of asymmetries found in recorded data (Grimmer and Seyfarth, Reference Grimmer and Seyfarth2014). The forward thrust force (Section 4.3) during stance of the intact limb allows for this near-symmetry, as it enables power injection during both stance phases. Instead, when this thrust force is completely removed, the only positive power injection occurs during stance phase of the residual leg. As a result, the stance phase of the intact leg lasts longer, and the gait becomes reminiscent of limping. In fact, by doing so, the average stance phase of the residual and intact leg were measured to last 65% and 45% of the stride cycle, respectively. Interestingly, despite this strong asymmetry, the model did not fall. However, the intended ambulation speed was neither tracked properly, as this was also regulated by the thrust force. The resulting ambulation speed and $ {\dot{q}}_1 $ of this simulation are also shown in the velocity plot of Figure 8, in transparent colors.

5.2.3. Introducing socket dynamics

As mentioned in Section 3, it is possible to introduce additional joints and bodies to pursue more realistic results, at the cost of increased simulation time. One such addition concerns the interaction between the patient’s stump and socket, which is currently modelled to be rigid, given the utilization of a single rigid body to model both (i.e., the thigh), but in reality is reported to be significantly soft (LaPrè et al., Reference LaPrè, Price, Wedge, Umberger and Sup2018). Prosthesis dynamics and control mostly act in the sagittal plane, for which the most relevant degrees of freedom are pistoning (longitudinal or vertical displacement) and anterior–posterior tilt (LaPrè et al., Reference LaPrè, Price, Wedge, Umberger and Sup2018). Instead, anterior–posterior displacement (forward horizontal displacement) is less relevant.

Two joints and bodies have been introduced between the hip and knee joints to model stump-socket softness for vertical displacement and anterior–posterior tilt, respectively. As a result, the dynamic system now consists of 10 instead of 8 joints: the stump-socket joints are joints 5 and 6, and the knee, tube, and ankle joints have moved down in the serial kinematic chain from joints 5 to 8 to joints 7 to 10. Per demonstration, The location of the two new joints are exactly halfway the thigh, and inertial properties of bodies have not been altered (the corresponding two new bodies have no mass). They are furthermore configured as forward dynamics joints, and per demonstration, their softness is modelled with translational and rotational linear spring-damper systems, so that

$$ {\displaystyle \begin{array}{c}{F}_{\mathrm{s}}\hskip0.35em =\hskip0.35em -{\tau}_5\hskip0.35em =\hskip0.35em {k}_5{q}_5+{c}_5{\dot{q}}_5\\ {}{M}_{\mathrm{s}}\hskip0.35em =\hskip0.35em -{\tau}_6\hskip0.35em =\hskip0.35em {k}_6{q}_6+{c}_6{\dot{q}}_6.\end{array}} $$

Whereas exact stiffness and damping gains are unknown and patient specific, here stiffness gains are chosen to realize any displacement difference on the order of $ d\hskip0.35em =\hskip0.35em 1\hskip0.35em \mathrm{cm} $ , in line with calculated displacement errors of markers in inverse kinematic studies (LaPrè et al., Reference LaPrè, Price, Wedge, Umberger and Sup2018). Damping gains are chosen to realize a system that is at least critically damped. As such, for the longitudinal displacement, we have selected the translational stiffness as $ {k}_5\hskip0.35em =\hskip0.35em mg/d\approx 74\hskip0.35em \mathrm{k}\hskip0.35em {\mathrm{Nm}}^{-1} $ and the translational damping as $ 2\sqrt{mk}\approx 4.7\hskip0.35em \mathrm{kN}\hskip0.35em \mathrm{s}\hskip0.35em {\mathrm{m}}^{-1} $ . For the rotational joint, we have selected $ {k}_6\hskip0.35em =\hskip0.35em 0.5m/d\approx 4.7\hskip0.35em \mathrm{kN}\hskip0.35em \mathrm{m}\hskip0.35em {\mathrm{rad}}^{-1} $ , to result in approximately 0.01 rad of displacement, and $ {c}_6\hskip0.35em =\hskip0.35em 2\sqrt{\frac{1}{4} mk}\approx 530\hskip0.35em \mathrm{N}\hskip0.35em \mathrm{m}\hskip0.35em \mathrm{s}\hskip0.35em {\mathrm{rad}}^{-1} $ , where $ \frac{1}{4}m $ is an upper estimate of effective inertia at the joint.

Without applying any further adjustments to the system, the controller or human intent, the simulation has been rerun. Generally, odds for instability are higher with increased softness, especially when the prosthesis and human controller are not further refined. However, results show that the model did not fall for the full duration of the 40 s trial. The corresponding simulation time has approximately doubled, but was still under 1 min. Nevertheless, when rerunning the simulation for stiffness values half that of the originally selected ones, a fall is observed around 11 s, after the increase of ambulation speed.

Results of socket displacement, forces and torques of the first 8 s for the originally selected stiffness values are shown in Figure 9. For comparison with previous results, the tube force and torque in Figure 8 have been plotted too. As expected, the force is very similar, but the moment is different. Nevertheless, the force signal is visibly smoother and its peak values have reduced with up to 20%, as a result of the new softness in the system.

Figure 9. Results demonstrating the effect of socket dynamics on the system, using otherwise similar properties and inputs as that of the simulation whose results are shown in Figure 7, for purposes of comparison.

5.3. Powered research leg

The second case study is that of an existing powered research prosthesis (Guercini et al., Reference Guercini, Tessari, Driessen, Buccelli, Pace, De Giuseppe, Traverso, De Michieli and Laffranchi2022). In particular, this research leg employs a dual motor architecture to assist multiple types of activities. It has a small, fast, yet low torque actuator, which is permanently engaged with the knee joint. In addition, it has a large, slow, but high torque actuator, which can optionally engage to provide additional torque in extension for activities that benefit from this, such as stair climbing. The use of only its small actuator should suffice to realize walking, despite its rated peak torque limit of approximately 30 Nm, whereas the knee torque of healthy subjects is observed to peak at approximately 1 NN $ {}^{-1} $ during walking (75 Nm for a subject of 75 kg). The presented simulation framework can be used to demonstrate this. Toward this purpose, an FSM-based controller inspired from that presented in Section 5.2 has been created, as shown in Figure 10. This controller has been used as the basis for the control system that has been presented for experimental validation by Guercini et al. (Reference Guercini, Tessari, Driessen, Buccelli, Pace, De Giuseppe, Traverso, De Michieli and Laffranchi2022), as realized with a Speedgoat interface in combination with a subject wearing an able-bodied adapter.

Figure 10. Overview of the FSM-impedance-based prosthesis controller module and electrically actuated research prosthesis device module, as designed and tuned to enable walking, starting and stopping. The FSM states are the same as those in Figure 6.

The prosthesis device is mainly different from that of the commercial quasi-passive hydraulic knee prosthesis (Section 5.2) in that the simplified hydraulic damping actuation model has been replaced for an electric drive train, featuring a model of a current-controlled DC motor and the gearing, which takes the input of a commanded current rather than a damping coefficient, as shown in Figure 10. The purpose of this implementation is to incorporate (1) realistic effects and limitations as a result of gearbox friction, fast shaft inertia and supply voltage and current limitations, (2) the possibility to analyze efficiency and heating as a result of joule and frictional losses. The end stops described in Section 5.2 is also included in the system of the research prosthesis.

To reduce simulation time, high-frequency dynamics are avoided by modeling the field-oriented controlled brushless Maxon EC 45 flat 24 V motor in continuous time and as a perfect current-commanded but saturated torque source that takes supply voltage and current limitations into account according to the equivalent brushed DC-motor model. With a torque constant of $ {k}_{\tau}\hskip0.35em =\hskip0.35em 36.9\hskip0.35em \mathrm{mN}\hskip0.35em {\mathrm{A}}^{-1} $ , a terminal resistance of $ R\hskip0.35em =\hskip0.35em 0.608 $ , a maximum permissible current of $ {I}_{\mathrm{max}}\hskip0.35em =\hskip0.35em 6.4\hskip0.35em \mathrm{A} $ (ca. double that of the specified nominal current) and a supply voltage of $ {V}_{\mathrm{max}}\hskip0.35em =\hskip0.35em 24\hskip0.35em \mathrm{V} $ , the theoretical maximum torque and no-load speed are $ {k}_{\tau }{I}_{\mathrm{max}}\hskip0.35em =\hskip0.35em 0.24\hskip0.35em \mathrm{N}\hskip0.35em \mathrm{m} $ and $ {U}_{\mathrm{max}}/{k}_{\tau}\hskip0.35em =\hskip0.35em 650\hskip0.35em \mathrm{rad}\hskip0.35em {\mathrm{s}}^{-1} $ _, respectively, on the fast shaft. The gearbox model incorporates drive-direction dependent (forward or backward) stiction, viscous friction and coulomb friction of the Conedrive CBC140100 strain wave gear with a gear ratio of $ {n}_{\mathrm{gear}}\hskip0.35em =\hskip0.35em 100 $ , and incorporates inertia of both itself and the motor rotor. The gearbox-motor unit connects to the joint via a relatively stiff linear spring element with a stiffness of 1000 N m rad $ {}^{-1} $ and a damping of 30 N m s rad $ {}^{-1} $ , to account for gearbox and joint flexibility. The knee torque $ {T}_{\mathrm{k}} $ is computed as the sum of the output torque of the gearbox $ {T}_{\mathrm{A}} $ and the end stop torque $ {T}_{\mathrm{extlim}} $ .

5.3.1. Results

At the time of controller development for the work by Guercini et al. (Reference Guercini, Tessari, Driessen, Buccelli, Pace, De Giuseppe, Traverso, De Michieli and Laffranchi2022), the focus was on realizing a single feasible walking gait for experimentation with an able-bodied user, without taking other actions into account or prosthesis energy consumption.

Following the procedure in Section 5.1, a main difference between the variable damper knee and powered knee is that inclusion of the electric actuator without any control output (step 2) leads to deterioration rather than improvement of walking behavior: not being able to complete a single step cycle, as shown in Figure 11. The presence of inertia and significantly more friction of the electric actuator cause the knee to extend too late and slow during swing, causing the system to trip. A fixed extension torque has been applied for knee angles $ {\theta}_{\mathrm{k}}>{55}^{\circ } $ to reliably overcome gearbox stiction and realize a repeatable transition from swing flexion to swing extension.

Figure 11. In contrast to the variable damping leg, the powered research leg shows impeded walking behavior when its actuator is introduced without a control output.

In addition, an impedance controller has been introduced that includes non-zero stiffness gains. The purpose in swing is countering friction to accelerate knee flexion and extension, especially at initiation of the walking gait, for the system would otherwise trip and fall. In stance, the purpose of these terms is providing additional support for braking. The introduction of additional damping gains was not necessary, as the gearbox friction already realizes sufficient damping on its own. Nevertheless, similar to the previous case study, damping has been introduced to brake the knee near extension ( $ {\theta}_{\mathrm{k}}<{15}^{\circ } $ ) in swing extension, to lower the impact torque with the extension end stop.

Results of the simulation for walking at a target velocity of $ 1.1\hskip0.35em {\mathrm{ms}}^{-1} $ starting from standstill are shown in Figure 12, corresponding to the designed and tuned controller shown in Figure 10. In this figure, $ {T}_{\mathrm{cmd}} $ corresponds to the commanded current $ {I}_{\mathrm{cmd}} $ , saturated by its absolute maximum of $ 6.4\hskip0.35em \mathrm{A} $ , multiplied by the torque constant $ {k}_{\tau } $ and gear ratio $ {n}_{\mathrm{gear}} $ , and as such represents a hypothetical or ideal torque output considering a drive train that is free of friction and inertia, which peaks at ca. 24 Nm. Impact torques with the extension end stop are still visible (e.g., at $ t\hskip0.35em =\hskip0.35em 3.1\hskip0.35em \mathrm{s} $ , 4.3 s, 5.6 s, etc.), but have been reduced by approximately 40% due to the inclusion of additional damping in swing extension.

Figure 12. Simulation results of an electrically powered research prosthesis, which demonstrate the possibility of realizing a walking gait including a knee bounce with a commanded knee torque $ {T}_{\mathrm{cmd}} $ that has a peak value (ca. 24 Nm) significantly lower than that found in a biological knee (ca. 75 Nm for a person weighing 75 kg). Due to gearbox friction and inertia, the actuator output torque $ {T}_{\mathrm{A}} $ is significantly different from the commanded $ {T}_{\mathrm{cmd}} $ . Spikes correspond to collisions with the knee extension end stop.

Foremost, the results illustrate the possibility to realize a possible walking gait despite the use of a weak actuator with limited torque output. Due to friction, the actual torque output $ {T}_{\mathrm{A}} $ is lower when forward driven, but higher when backward driven, or during rapid accelerations as caused by drive train inertia, for which it is significantly different from $ {T}_{\mathrm{cmd}} $ . The latter is visible during the aforementioned collisions with the extension end stop, which—in combination with the relatively stiff drive train—result in a high reaction torque on the actuator ( $ {T}_{\mathrm{A}} $ ). This provides an argument against the use of stiff (non-series elastic actuated) electric actuation with high inertia in knee prostheses.

An analysis of energy flows of the actuator can be performed to illustrate the sub-optimal performance of the present hardware and control solution. These can be calculated by integrating the mechanical and electrical power flows from its various components over time, which is preferably done for a whole number of steps. Average power can then be calculated by dividing the resulting energy figures by the total duration $ \Delta t $ of these steps. The actuator was shown to damp the system by doing negative work, calculated as $ \int {\dot{q}}_5{T}_{\mathrm{A}}\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em -9.9\hskip0.35em \mathrm{W} $ , averaging five whole steps from the third step onward, from $ t\hskip0.35em =\hskip0.35em 3.52\hskip0.35em \mathrm{s} $ to $ t\hskip0.35em =\hskip0.35em 9.24\hskip0.35em \mathrm{s} $ ( $ \Delta t\hskip0.35em =\hskip0.35em 5.72\hskip0.35em \mathrm{s} $ ). Yet, the motor provides a net mechanical power of $ \int {P}_{\mathrm{rotor}}\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em 2.6\hskip0.35em \mathrm{W} $ to the actuator, where the motor rotor mechanical power $ {P}_{\mathrm{rotor}}\hskip0.35em =\hskip0.35em {\omega}_{\mathrm{rotor}}{T}_{\mathrm{rotor}} $ , with $ {T}_{\mathrm{rotor}}\hskip0.35em =\hskip0.35em {k}_{\tau }I $ the applied torque on the motor rotor (fast shaft), and $ {\omega}_{\mathrm{rotor}}\hskip0.35em =\hskip0.35em {n}_{\mathrm{gear}}{\dot{q}}_5 $ its rotational velocity. This suggests $ 2.6+9.9\hskip0.35em =\hskip0.35em 12.5\hskip0.35em \mathrm{W} $ of heat losses in the transmission alone. Nevertheless, the power flows in the motor are significantly higher than the net average figure of 2.6 W. In particular, it is the sum of positive work of $ \int \max \left({P}_{\mathrm{rotor}},0\right)\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em 10.9\hskip0.35em \mathrm{W} $ and negative work of $ \int \min \left({P}_{\mathrm{rotor}},0\right)\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em -8.3\hskip0.35em \mathrm{W} $ . However, the main issue is with electric energy consumption, calculated as $ \int \max \left({P}_{\mathrm{elec}},0\right)\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em 15.8\hskip0.35em \mathrm{W} $ , of which up to $ -\int \min \left({P}_{\mathrm{elec}},0\right)\hskip0.35em \mathrm{d}t/\Delta t\hskip0.35em =\hskip0.35em 6.3\hskip0.35em \mathrm{W} $ is possibly eligible for regeneration,Footnote ⁶ where the electrical supply power $ {P}_{\mathrm{elec}} $ is estimated according to the DC-motor equivalent $ {P}_{\mathrm{elec}}\hskip0.35em =\hskip0.35em {I}^2R+{IU}_{\mathrm{emf}} $ , with the counter-electromotive force $ {U}_{\mathrm{emf}}\hskip0.35em =\hskip0.35em {\omega}_{\mathrm{rotor}}{k}_{\tau } $ , and $ \int {I}^2R\hskip0.35em \mathrm{d}t\hskip0.35em =\hskip0.35em 6.9\hskip0.35em \mathrm{W} $ is permanently lost to joule heating in the windings of the motor. Depending on energy regenerative capabilities of the electronics, this suggests an efficiency for the motor alone of only 16% to 27%. Other control strategies can likely improve this figure. Nevertheless, also regarding hardware, it is noteworthy that—despite the specific design choice of using two motors of which only the smaller and faster one is used for passive activities such as walking—it still suffers from high inertial reaction torques and friction for which the motor partly needs to compensate, in order to realize a (transparent) walking gait. An optimization of not only the controller, but also the hardware—featuring a faster and more efficient drive train still—could lead to significant reduction of energy consumption.

The simulation provides a basis for an initial feasible control scheme to perform walking experiments of the research leg (Guercini et al., Reference Guercini, Tessari, Driessen, Buccelli, Pace, De Giuseppe, Traverso, De Michieli and Laffranchi2022). However, after its implementation, whereas most parameters are kept unchanged, changes to the controller are nonetheless made on the real device. After all, the presented model remains a simplification of reality, in particular regarding human behavior. Whereas the model produces feasible results for a certain scope of possible human behaviors, it is not guaranteed that the human will exhibit exactly those behaviors in reality, nor does the model take human psychology into account. Other reasons for controller changes include unknown or imperfectly modelled sources of friction or device imperfections, leading to the removal or inclusion of additional friction compensation terms to the implemented controller.