## Abstract

Achieving coordinated motion between transfemoral amputee patients and powered prosthetic joints is of paramount importance for powered prostheses control. In this article, we propose employing an algebraic curve representation of nominal human walking data for a powered knee prosthesis controller design. The proposed algebraic curve representation encodes the desired holonomic relationship between the human and the powered prosthetic joints with no dependence on joint velocities. For an impedance model of the knee joint motion driven by the hip angle signal, we create a continuum of equilibria along the gait cycle using a variable impedance scheme. Our variable impedance-based control law, which is designed using the parameter-dependent Lyapunov function framework, realizes the coordinated hip-knee motion with a family of spring and damper behaviors that continuously change along the human-inspired algebraic curve. In order to accommodate variability in the user's hip motion, we propose a computationally efficient radial projection-based algorithm onto the human-inspired algebraic curve in the hip-knee plane.

## 1 Introduction

Coordinating the motion of transfemoral amputee patients and powered prosthetic joints is of vital importance for powered prostheses control. In order to achieve such coordinated motion patterns, the two important issues of (i) representing the human walking gait and (ii) enforcing the desired walking gait via a proper control scheme need to be addressed.

In the rehabilitation robotics literature, there are two classes of human gait cycle representation. The first class, proposed in Refs. [1]–[6], divides the gait cycle into several periods, each with their own distinct controllers. In order to enforce the desired gait motion patterns, switched impedance schemes are employed in a way that the closed-loop dynamics are forced to behave with a series of *finite states* consisting of passive spring and damper behaviors [1,6]. Due to the passive nature of the closed-loop dynamics, the advantage of such impedance-based schemes is the inherent stability of the interaction between the powered prosthesis and the amputee in each state of the finite state machine.

Despite the advantages of the proposed switched impedance-based strategies in Refs. [1]–[6], one shortcoming of such schemes is that they rely on nonunified representations of walking gait. Consequently, there are dozens of control parameters and transition rules that must be tuned across users and activities [5]. In order to automate the control parameter tuning procedure, techniques such as fuzzy logic-based expert systems [7] and adaptive dynamic programming [8] have been proposed in the literature. A second shortcoming stems from the fact that desynchronizing perturbations increase patient's risk of falling if the nonunified impedance-based controller switches to the wrong state with the wrong control scheme. Indeed, switched impedance-based control techniques that create *isolated equilibria* for the powered prosthetic joint dynamics rely on being able to estimate the precise time for switching from one finite state to another. A potential remedy for the need to determine switching times is based on creating a *continuum of equilibria* along the normal human walking gait and continuously evolving on the continuum, instead of having a series of isolated equilibria.

In contrast to the first class, there is a second class of wearable robot controllers that employ unified representations of the entire human gait cycle [9–16]. Unified gait controllers do not suffer the aforementioned shortcomings of the nonunified gait control schemes. In Ref. [9], the authors employed the center of pressure to unify the stance phase of the gait. Assuming a rocker foot geometry, the center of pressure-based gait in Ref. [9] is a *holonomic* representation of the human stance period. A holonomic gait, which separately unifies the stance phase and swing phase of walking, was proposed in Ref. [17]. The gait in Ref. [17], however, uses the patient's Cartesian coordinates as the phase variables and thus is difficult to measure with onboard sensors. In order to remove the need for measuring the patient's Cartesian coordinates, the authors in Ref. [16] have proposed two different representations using the hip angle measurement. The first one is a piecewise holonomic representation of the human walking gait cycle.

The second representation in Ref. [16], which unifies the stance and swing phases of walking, is a velocity-based (nonholonomic) gait representation using the thigh phase portrait. This second representation method was later experimentally verified in Ref. [18] via a single sensor. As another unified walking gait cycle representation, the authors in Refs. [12] and [13] proposed a thigh angle integral-based representation of walking. However, the integral term in Refs. [12] and [13] needs to be reset at the beginning of every gait cycle to prevent accumulation of drift due to variation in thigh kinematics. Furthermore, the integral-based gait in Refs. [12] and [13] cannot be used for patient's nonrhythmic motion control. In order to automate the phase-based control parameter tuning in real-time, the authors in Ref. [19] proposed using an extremum seeking control, which is run in parallel to the phase-based control scheme in Refs. [12] and [13].

Motivated by the fact that healthy hip-knee angular trajectories create closed and nonintersecting orbits during walking (see Fig. 1), we employed algebraic curves and the 3 L algorithm from the pattern recognition literature [20] to generate a unified holonomic representation of the human gait cycle in our preliminary work [21]. In particular, our human gait cycle representation, which encodes the nominal coordination between the hip and knee angles during normal level walking, is the zero set of bivariate implicit polynomials (IPs). We remark that our proposed algebraic curve fitting in Ref. [21] is closely related to the recent work in Ref. [22], where the authors use an elliptical path in the hip-knee plane to prescribe a normal walking gait for a lower limb exoskeleton during the swing phase. Our algebraic curve representation, however, unifies the entire gait cycle, including both the stance phase and the swing phase. Furthermore, unlike the elliptical path in Ref. [22], which is an open curve converging to infinity, our proposed algebraic curves are closed and bounded, which correspond to periodic human gait cycles. The closed and bounded curve properties cannot be achieved with standard one-dimensional polynomial fitting methods, as the one in Ref. [22]. In an earlier work [23], we employed a symbolic algebraic tool, which is based on computing the resultant of polynomials, for removing phase variables from *autonomous bipedal robot* parametric gaits. However, this approach, similar to Ref. [22], cannot be *practically* used for generating closed algebraic curves from human walking data due to *large degrees* of generated polynomials.

Based on our preliminary work in Ref. [21], we propose employing the algebraic curve representation of nominal walking data for powered knee prosthesis controller design in this article. In order to take into account the variability that is observed during walking, we improve our earlier human-inspired algebraic curve generation algorithm in Ref. [21]. In particular, in order to take into account the fact that hip-knee walking profiles vary during walking, we employ variable geometric contraction and dilation factors, in contrast to constant values in our preliminary work. The variable contraction/dilation factor ensures that the level-sets of the generated implicit polynomials have sufficient distance from each other in a way that step-by-step variability does not significantly perturb the algebraic distance from the normal walking data. Having obtained a unified holonomic representation of walking data via algebraic curves, in the next step, we consider an impedance model of the knee joint motion driven by the hip angle signal. We then create a continuum of equilibria along the human-inspired algebraic curve using a variable impedance scheme. Our variable impedance-based control law, which is designed using the parameter-dependent Lyapunov function framework, realizes the coordinated hip-knee motion with a family of spring and damper behaviors that *continuously* change along the human-inspired algebraic curve. In order to accommodate the variability in user's hip motion, we propose a radial projection-based algorithm onto the human-inspired algebraic curve in the hip-knee plane. The algebraic curve representation makes the radial projection algorithm computationally efficient, with guaranteed convergence in a finite number of steps.

The rest of this paper is organized as follows: first, we briefly review preliminaries from algebraic curves and present the 3 L algorithm for fitting such curves to nominal hip-knee walking data in Sec. 2. Next, in Sec. 3, for an impedance model of the knee joint motion driven by the hip angle signal, we present our control problem formulation and demonstrate that the knee joint closed-loop dynamics take the form of a linear parameter varying (LPV) dynamical system, which is subject to disturbances that change continuously along the human-inspired algebraic curve and are dependent on the hip joint motion speed. Thereafter, in Sec. 4, we propose a parameter-dependent Lyapunov function approach for designing variable stiffness and damping gains along the human-inspired algebraic curve. Additionally, we present a radial projection-based algorithm in order to address variability in the user's hip motion and show the simulation results. In Sec. 5, we present our numerical results. Finally, we conclude this paper with final remarks and outline of possible future research directions in Sec. 6.

*Notation*: Given two vectors (matrices) *a* and *b*, we denote by [*a*; *b*] the vector (matrix) $[a\u22a4,b\u22a4]\u22a4$ where $(\xb7)\u22a4$ is the transpose operator. Given a square symmetric matrix *X*, we denote its maximum (minimum) eigenvalue by $\lambda max(X)$ ($\lambda min(X)$).

## 2 3 L Algorithm for Fitting Human Data to Algebraic Curves

In this section, we review some preliminaries on algebraic curves and use the 3 L algorithm for fitting closed algebraic curves to human data. The 3 L algorithm presented in this section is an extension of the 3 L algorithm in our preliminary work in Ref. [21]. A comprehensive treatment of algebraic curves and their properties may be found in Refs. [24] and [25].

### 2.1 Algebraic Curves.

*n*, an IP

*a*are real numbers. The

_{ij}*degree of polynomial*$h(qH,qK)$ in Eq. (1) is the maximal value of $i+j$ for which $aij\u22600$. Here, we assume that IP

*h*is of degree

*n*. Every IP of degree

*n*has $(n+1)(n+2)/2$ coefficients.

*algebraic distance*of the point $(qH0,qK0)$ to the zero set of $h(qH0,qK0)$. Moreover, the

*zero set*of IP $h(qH,qK)$ given by Eq. (1) is defined to be

A *real algebraic curve* is the zero set of a nonzero real bivariate polynomial *h*. The *degree of an algebraic curve* is defined to be the degree of its associated bivariate implicit polynomial. Algebraic curves of degree 1, 2, 3, 4, $\u2026$, are called *lines*, *conics*, *cubics*, *quartics*, $\cdots $, respectively.

Since we are interested in studying periodic human walking gait profiles, it is desirable to have closed and bounded algebraic curves. The following well-known lemma provides the necessary condition for having a closed and bounded algebraic curve.

Lemma 1. (see Ref. [25]) *An algebraic curve is a closed and bounded plane curve only if it is of even degree*.

### 2.2 3 L Algorithm for Fitting Algebraic Curves to Human Walking Data.

In this section, we fit algebraic curves to the data points obtained from a healthy gait according to Winter's normal cadence walking data [26]. Our fitting algorithm is taken from the pattern recognition literature and is known as the 3 L algorithm [20]. The 3 L algorithm presented in this section is a generalization of the algorithm in our preliminary work in Ref. [21].

#### 2.2.1 Fitting Algebraic Curves to Human Datasets as a Quadratic Optimization Problem.

*N*

_{0}planar data points $(qHl,qKl)$. In the set $H0$, if $(qHl,qKl)$ corresponds to time instant

*t*and $(qHl+1,qKl+1)$ corresponds to time instant $tl+1$, then $tl<tl+1$. Here, the set $H0$ represents the path in the hip-knee plane, which is taken from the Winter's normal cadence walking data [26] (see Fig. 1). The set of ordered data points $(qHl,qKl),\u20091\u2264l\u2264N0$, is the samples of the hip and knee normal cadence walking trajectories corresponding to increasing time instants during a given gait cycle. Since the walking gait profile is periodic, we have $(qH1,qK1)=(qHN0,qKN0)$. The geometric center or centroid of the dataset $H0$ is defined to be the point (see also Fig. 1)

_{l}*p*, such that its zero set $Z(h)$ approximates the set of human hip-knee data points $H0$. This approximation problem is equivalent to minimization of the error functional [20]

*E*in Eq. (4) as a quadratic function of the coefficients of IP $h(qH,qK)$. In order to do so, we rewrite IP $h(qH,qK)=\u2211ijaijqHiqKj$, where $0\u2264i+j\u2264n$, using the inner-product

Therefore, we would like to find the coefficient vector *a* such that the error functional *E* given by Eq. (10) is minimized. The direct minimization of *E*, which is known as the 1 L approach, often fails to generate an acceptable IP fit to a set of data points due to numerical instability problems and/or lack of physically meaningful solutions (see Refs. [20] and [27] for a detailed discussion).

#### 2.2.2 Human Data 3L Fitting Algorithm.

Using the 3 L fitting algorithm developed in Ref. [20], we address the aforementioned numerical shortcomings of the 1 L minimization by introducing two additional datasets, which can be generated from the original human walking dataset $H0$. Being able to vary the two fictitious datasets on either side of the original normal cadence walking dataset $H0$, optimal fitting accuracies of algebraic curves to human walking data can be achieved. Furthermore, it can be shown that under suitable conditions, the algebraic curves generated by the 3 L fitting algorithm are nondegenerate [20,28,29].

Following the experimental results for pose estimation in computer graphics literature [20,25], we chose quartic algebraic curves, i.e., *n *=* *4, to be fitted to the hip-knee normal walking data. However, there is no limitation on the degree of the algebraic curve that can be fitted to human walking datasets. The 3 L algorithm for fitting algebraic curves to human hip-knee walking data can be described in the following three steps.

##### 2.2.2.1 Step 1: Generation of two fictitious datasets.

Given that the dataset $H0$, representing the nominal human walking gait in the hip-knee plane [26], introduces two fictitious datasets close to $H0$. The first set, which is denoted by $H+$ and is located outside the dataset $H0$, has $N+$ points and corresponds to the algebraic distance *c* (see Sec. 2.1 for the definition of algebraic distance), where *c* is a design parameter to be chosen. The second fictitious dataset, which is denoted by $H\u2212$ and is located inside the dataset $H0$, has $N\u2212$ points and corresponds to the algebraic distance—*c*. In this article, we have chosen the design parameter to be *c *=* *1.

where $[qHC,\u2009qKC]\u22a4$ is the centroid of the human dataset $H0$.

In our preliminary work in Ref. [21], we scaled the translated dataset $H0t$ by *constant* scaling factors $\alpha +>1$ and $0<\alpha \u2212<1$ to generate $H+$ and $H\u2212$, respectively. The multiplication of the translated dataset $H0t$ by $\alpha +$ and $\alpha \u2212$ corresponds to *geometric contraction* and *geometric dilation*, respectively. For the presented results in Ref. [21], we chose $\alpha +=1.02$ and $\alpha \u2212=0.98$.

*l*

_{1}and

*l*

_{2}correspond to the walking data at 50% and 75% of the gait cycle, and $\Delta \alpha 1\xb1$ and $\Delta \alpha 2\xb1$ enlarge the distance between the two fictitious datasets with $H0$.

##### 2.2.2.2 Step 2: Defining the three level-set matrix.

##### 2.2.2.3 Step 3: Solving for the unknown implicit polynomials coefficients.

where $m\u22a4(\xb7,\xb7)$ is the function defined by Eq. (6) and $pc=[qHC,\u2009qKC]\u22a4$ is the centroid of the human dataset $H0$. The translation in Eq. (17) corresponds to translation of the geometric center of the obtained algebraic curve from the origin to the human dataset centroid, as explained in Remark 1. The human-inspired IP $h\u22c6(qH,qK)$ represents a desired implicit relationship between the human's hip angle and the wearable robot knee angle. As the human's hip angle $qH(t)$ evolves with time, driving $h\u22c6(qH(t),qK)$ to zero via feedback corresponds to coordinating the motion of the knee with the hip during level walking according to the human walking closed curve in Fig. 1. The algebraic curve given by IP $h\u22c6(qH,qK)$ represents a holonomic relationship as it only depends on the joint angles $qH$ and $qK$, as opposed to the joint angle velocities.

#### 2.2.3 Discussion of the Obtained Algebraic Curve Fits.

Using the 3 L algorithm with two different types of fictitious datasets, we obtained two quartic IPs whose zero level sets are depicted in Fig. 2. Figure 2(a) corresponds to choosing variable geometric factor functions while Fig. 2(b) corresponds to choosing a constant geometric factor coefficient in step 1 of our algorithm. As it can be seen from Fig. 2, the shape of the level sets and their values at different configurations differ from each other in Figs. 2(a) and 2(b). Moreover, the level sets in Fig. 2(a), which is based on the improvement over our preliminary work in Ref. [21], are further from each other during the stance to swing phase transition. Furthermore, if the hip-knee configurations are located further from the nominal Winter's walking data [30], the level set values in Fig. 2(a) in comparison with Fig. 2(b) do not grow very large.

The largest deviation of the quartic algebraic curve fits from Winter's normal cadence walking data, for both of the fits, happens during the stance flexion/extension phase at around the configuration $[0.2677\u2009rad,\u22120.081\u2009rad]\u22a4$, where the knee angle deviates from the nominal value by around $0.04\u2009rad\u22482.3\u2009deg$ for both the cases. The gray band around the zero set of IP $h\u22c6(qH,qK)$ corresponds to the hip and knee configurations that belong to the sublevel set ${(qH,qK):|h\u22c6(qH,qK)|\u2264c0}$ associated with the algebraic distance $c0=0.5$. As it can be seen from Fig. 2, the level sets of the fitted IPs $h\u22c6(qH,qK)$ do not intersect with each other, corresponding to the fact that the generated algebraic curves are nondegenerate. Furthermore, as the algebraic distances from the zero set $Z(h\u22c6)$ change in a small manner, the joint angles do not deviate drastically from the nominal values. This continuity feature is desirable for controller design; since for small values of $|h\u22c6(qH,qK)|$, the hip-knee configurations are still close to the fitted algebraic curve.

## 3 Control Problem Formulation

In this section, we consider an impedance model of the knee joint motion driven by the hip angle signal. Our control problem formulation is motivated by the approach in Refs. [1] and [6], where the impedance-based feedback control is employed to render the closed-loop dynamics of the powered prosthetic joints as a series of passive spring and damper behaviors. Due to the passive nature of the closed-loop dynamics under impedance-based control laws, the interaction between the powered prosthesis and the amputee will be inherently stable. However, we depart from the control framework in Refs. [1] and [6] by creating a continuum of equilibria along the human-inspired algebraic curve instead of having a series of isolated equilibria. As it will be shown in this section, the resulting closed-loop knee joint motion can be described by a LPV dynamical system.

*u*is the torque applied to the powered knee prosthesis. Furthermore, we let the polar representation of the human-inspired algebraic curve in Eq. (17) be

*σ*.

where $ue(\sigma ):=(kkn/Jkn)qK*(\sigma )$. The control input $ue(\sigma )$ renders $xe(\sigma )$ an equilibrium for the knee joint dynamics.

where $Acl(\sigma ):=A\u2212BKc(\sigma ).$

*frozen equilibria*, i.e., whenever the hip comes to rest, we have $\sigma \u02d9=0$. Consequently, from Eq. (25), it follows that $de(\sigma ,\sigma \u02d9)=0$. Therefore, under the frozen equilibria condition, the knee closed-loop error dynamics take the LPV form given by Eq. (29).

## 4 Variable Impedance Control Design

Having formulated our control problem in Sec. 3, we first provide a parameter-dependent Lyapunov function approach for designing variable impedance gains along the human-inspired algebraic curve in Sec. 4.1. The results in Sec. 4.1 are applicable when the human hip evolves according to a nominal gait profile on the algebraic curve. In order to address variability in the user's hip motion, we present a radial projection-based algorithm in Sec. 4.2.

### 4.1 Variable Impedance Control Design Using Parameter Dependent Lyapunov Functions.

In this section, we assume that the hip joint evolves, with a bounded rate of change, along the algebraic curve according to $qH(t)=qH*(\sigma (t))$, where $qH*(\sigma )$ is given by the polar representation in Eq. (19). Under this assumption, we would like to make the tracking error in Eq. (21)*practically s**table* [32]. Practical stability of the error dynamics means that the error can be made arbitrarily small via a proper choice of the state feedback gain $Kc(\sigma )$.

where the variable stiffness $Kp(\sigma )$ and the variable damping $Kd(\sigma )$ are some smooth functions of *σ*, which continuously vary along the human-inspired algebraic curve.

*does not depend on a reference velocity*, e.g., on $q\u02d9K*(\sigma )$. Additionally, under the proposed control law in Eq. (32), the closed-loop state matrix takes the form

Since the stiffness $Kp(\sigma )$ and the damping $Kd(\sigma )$ functions are $2\pi $-periodic by design, it follows that the closed-loop state matrix $Acl(\sigma )$ satisfies the periodicity condition $Acl(\sigma +2\pi )=Acl(\sigma )$, for all values of the angular variable *σ*.

where $X(\sigma )$ is a symmetric and positive definite matrix, which is parameterized by the polar angle *σ*. Furthermore, since *σ* is an angular variable, the matrix $X(\sigma )$ should be chosen to satisfy the periodicity condition $X(\sigma +2\pi )=X(\sigma )$, for all values of *σ*.

*σ*, which determines where the knee joint should be located along the human-inspired algebraic curve. In order to proceed further, we take the derivative of $V(x\delta ,\sigma )$ along the trajectories of Eq. (28) and obtain

*t*. We also assume that the rate of change of the polar angle satisfies

*t*, where ρ is a positive constant. Under the bounded rate of change condition in Eq. (36), the disturbance term given by Eq. (25) satisfies

The values of $|dqK*(\sigma )/d\sigma |$, which will be used for designing the variable impedance gains later, are depicted versus the percentage of stride in Fig. 3.

When the rate of change of the variable *σ* is bounded but an upper bound is not known a priori, the following proposition provides a stability condition, based on a Lyapunov matrix inequality.

*Consider the closed-loop dynamics in Eqs. (28) and (33) and assume that the matrix*$X(\sigma )$

*in Eq. (34) is equal to a symmetric and positive definite constant matrix X. Furthermore, assume that the rate of change of the variable σ is bounded by*$|\sigma \u02d9(t)|\u2264\rho $

*for some unknown constant ρ. If the Lyapunov inequality*

*is satisfied for all σ, then*,

- (i)
*the continuum of equilibria parameterized by σ in**Eq. (22)**is exponentially stable under*$\sigma \u02d9=0$*; and,* - (ii)
*when*$\sigma \u02d9\u22600$*, the error*$x\delta (t)$*is bounded according to*$V(x\delta (t),\sigma (t))\u2264V(x\delta (0),\sigma (0))exp(\u2212\zeta (\sigma )t)+(dqK*(\sigma )/d\sigma )2\zeta (\sigma )\rho 2$(39)*for all time instants t*> 0,*where*$\zeta (\sigma ):=\lambda min(Q(\sigma ))\lambda max(X)\u2212\lambda max(X)$.

*Proof*. We provide a sketch of the proof. From the inequality in Eq. (38) and since the matrix $X(\sigma )=X$ is constant, it follows that:

Therefore, since $de(\sigma ,\sigma \u02d9)=0$ when $\sigma \u02d9=0$, (i) follows. Statement (ii) follows from a direct application of the comparison lemma [34] to the above inequality. $\u25a0$

Proposition 1 provides conditions that can be used for designing the variable stiffness and damping gains along the human-inspired algebraic curve. From standard results in the nonlinear systems literature (see, e.g., Theorem 4.6 in Ref. [34]), the Lyapunov matrix inequality in Eq. (38) is satisfied if and only if for all *σ* the eigenvalues of the matrix $Acl(\sigma )$ given by Eq. (33) are located in the open left-half complex plane. This objective will be achieved if and only if $ap(\sigma )>0$ and $ad(\sigma )>0$ for all values of *σ*. As it will be discussed in Sec. 5, attenuation of the disturbance $de(\xb7)$ in Eq. (23), which exists due to continuous movement of the operating equilibrium point along the algebraic curve during walking, will be achieved if both impedance gains $Kp(\sigma )$ and $Kd(\sigma )$ change continuously with $||de(\xb7)||$ along the algebraic curve.

### 4.2 Extension of Control Design Under Variability in the User's Hip Motion.

Under the nominal hip joint motion, i.e., when $qH(t)=qH*(\sigma (t))$, the reference knee joint angle is given by $qK(t)=qK*(\sigma (t))$ such that the pair $(qH*(\sigma (t)),qK*(\sigma (t)))$ belongs to the human inspired algebraic curve, i.e., $h\u22c6(qH*(\sigma (t)),qK*(\sigma (t)))=0$. However, when the user's hip motion undergoes variability and is no longer given by the nominal function $qH=qH*(\sigma )$, determining the reference knee joint angle for the variable impedance scheme will become nontrivial.

*p*(

*t*) evolves in a neighborhood of $Z(h\u22c6)$, it is possible to project it onto the human-inspired algebraic curve and obtain the reference point

where $\pi Z(h\u22c6):\mathbb{R}2\u2192Z(h\u22c6)$ is a proper projection mapping, to be designed, onto the human-inspired algebraic curve.

A proper projection mapping $\pi Z(h\u22c6)$, which provides the reference points on the human-inspired algebraic curve during walking, should enjoy several properties: (i) for each complete revolution of the pair $p(t)=(qH(t),qK(t))$ on a closed curve in a neighborhood of the human-inspired algebraic curve, the projected point $\pi Z(h\u22c6)(p(t))$ should undergo a complete revolution on the human-inspired algebraic curve $Z(h\u22c6)$; (ii) each point on the human-inspired algebraic curve should be mapped onto itself via the projection mapping, i.e., $\pi Z(h\u22c6)(p)=p$ for all $p\u2208Z(h\u22c6)$; and (iii) the projection mapping computation should be fast enough for control implementation purposes.

*p*and the centroid of the human-inspired algebraic curve $Z(h\u22c6)$. It is immediate to see that the projection mapping $\pi Z(h\u22c6)(\xb7)$ in Eqs. (43) and (44) satisfies properties (i) and (ii). Figure 4 depicts the result of the proposed radial-based projection algorithm for a deviated walking stride onto the human-inspired algebraic curve.

The following proposition states that if the point *p* belongs to a sufficiently small neighborhood of the algebraic curve, then the radial projection mapping can be efficiently computed using the classical bisection root finding algorithm [38].

Proposition 2. *Consider the human-inspired algebraic curve and its centroid. Consider the projection mapping*$\pi Z(h\u22c6)(\xb7)$*in Eqs. (43) and (44). In a sufficiently small neighborhood of*$Z(h\u22c6)$*, the projection mapping*$\pi Z(h\u22c6)(\xb7)$*is well-defined and can be computed via the bisection root finding algorithm.*

*Proof*. Due to the geometry of the level-sets of the IP $h\u22c6(qH,qK)$, the zero set of $h\u22c6(\xb7)$ is an isolated zero level set. In other words, there exists an open neighborhood of the closed curve $Z(h\u22c6)$ denoted by $N$ such that if $q\u2208N$ and $h\u22c6(q)=0$, then $q\u2208Z(h\u22c6)$. Consider an arbitrary point $p\u2208N$ and suppose that $fp(0)\u22600$. Without loss of generality, we assume that $fp(0)>0$. In other words, *p* belongs to the exterior of $Z(h\u22c6)$. Since $Z(h\u22c6)$ is the only zero set of $h\u22c6(\xb7)$ in the neighborhood $N$ and $p\u2208N$, it follows that the line $\u2113p$ in Eq. (45), which passes through *p* and the centroid of the closed curve $Z(h\u22c6)$, intersects the closed curve $Z(h\u22c6)$ at a unique point $p0*$. Suppose that $p0*=p+\tau 0*p\u2212pc||p\u2212pc||$. It follows that $fp(\tau 0*)=0$. It is now possible to choose a point $p\u2032\u2208N$ on line $\u2113p$, associated with the real number $\tau \u2212<0$, in the interior of $Z(h\u22c6)$ such that $fp(\tau \u2212)<0$.

Since $fp(\tau )$ is a smooth function of *τ*, we have $fp(0)>0,\u2009fp(\tau \u2212)<0$. Furthermore, $\tau 0*$ is the only root of the function $fp(\tau )$ in the interval $[\tau \u2212,0]$. From standard results in numerical analysis (see, e.g., Ref. [38]), the proof will be concluded.$\u25a0$

## 5 Simulation Results

We implemented the impedance-based gait control strategy on an impedance model of the knee dynamics of a human subject. In particular, the impedance modeling of knee dynamics was based on the gait data of a healthy 75 kg subject, as derived from body-mass-normalized data from Winter [30]. Using a variable hip joint motion time profile, generated by the method proposed in Ref. [39], we investigated the effectiveness of our continuous impedance-based scheme and the radial-based projection algorithm in Sec. 4.2. In all of our results, we are working with the algebraic curve shown in Fig. 2(a), which is an improvement over our preliminary work in Ref. [21] shown in Fig. 2(b).

### 5.1 Algebraic Distances Under Hip-Knee Variability.

In order to study the variations of the generated IP level sets under hip-knee variability during human walking, we generated 100 walking strides using the methodology in Ref. [39]. These profiles are depicted in Fig. 5(a). The values of the IP $h\u22c6(qH,qK)$ on the generated walking data are depicted in Fig. 5(b). As it can be seen from Fig. 5(b), the algebraic distances given by $h\u22c6(\xb7)$ on the walking strides located at one standard deviation above and below the normal walking data [30] are bounded by the maximum algebraic distance $c0=4$, which happens during the stance to swing transition (around 60% of the walking stride). As it can be seen from Figs. 2(a) and 5(a), this algebraic distance corresponds to a geometric distance of around 5 deg to the nominal hip-knee configuration. Therefore, the generated IP is not overly sensitive to normal variability in human walking data.

### 5.2 Design of Variable Impedance Gains.

where $||\xb7||2$ denotes the *L*_{2}-norm. Moreover, $Np=4$ and *N _{d}* = 2.

Figure 6(a) depicts the normalized variable impedance gains by dividing them by ρ, which is the upper bound on the rate of change of the polar angle (see Eq. (36)), versus the percentage of normal walking stride. The stiffness to damping ratios in our continuous impedance-based scheme are the highest during the stance flexion/extension and swing extension phases of walking and are the lowest during the preswing and swing flexion phases. Interestingly, the stiffness to damping ratios in the conventional switched impedance-based schemes, first proposed in Refs. [1] and [6], also assume their highest and lowest values during the similar phases of walking.

### 5.3 Implementation of the Variable Impedance-Based Scheme.

The block diagram of the overall scheme is depicted in Fig. 7. The projection mapping given by Eqs. (43)–(45) provides the variable impedance-based controller with the reference signal $qK*$ and the polar angle *σ*. In turn, the variable impedance controller, whose stiffness and damping gains are given by Eq. (40), drives the error $qK\u2212qK*(\sigma )$ toward zero, and by doing so, drives the algebraic curve distances to zero.

We tested our variable impedance-based scheme for five walking strides deviated from the Winter's normal motion profile under human hip joint variability. Figure 6(b) depicts the knee time-profiles and their comparison to the normal walking profile (the black curve). Figure 6(c) depicts the resulting traversed paths in the hip-knee plane.

## 6 Conclusions and Future Research Directions

In this article, we employed algebraic curves to represent human walking gait data and achieve coordinated motion between transfemoral amputee patients and powered prosthetic joints. For an impedance model of the knee joint motion driven by the hip angle signal, we proposed using a variable impedance scheme for creating a continuum of equilibria along the human-inspired algebraic curve. Our variable impedance-based control law, which is designed using the parameter-dependent Lyapunov function framework, realized the coordinated hip-knee motion with a family of spring and damper behaviors that continuously change along the human-inspired algebraic curve. Exploiting the geometrical structure of the level-sets of the human-inspired algebraic curve, we proposed a computationally efficient radial projection-based algorithm onto the algebraic curve in the hip-knee plane. We employed the radial projection-based algorithm to accommodate the variability in the human's hip motion in our variable impedance-based scheme. The presented material opens up some potential research directions.

From a theoretical perspective, the 3 L fitting algorithm or its extensions *might* be employed for fitting trivariate (of three variables) implicit polynomials to hip-knee-ankle human walking data. Moreover, the variable impedance-based scheme needs to consider the nonlinear dynamics of a powered knee-ankle prosthesis and use a proper projection algorithm in order to take into account the variability in human hip motion.

From a control system implementation perspective, the proposed methodology has the *potential* to be employed for nonrhythmic motions and thus overcome the limitations of existing unified gait control methods for powered prostheses. As the proposed radial projection algorithm can be implemented using a bisection root finding numerical method, it has the potential to be efficiently implemented on an embedded system in real time. Finally, we can fit the algebraic curve, using the same 3 L numerical algorithm, to a global hip measurement for implementation with an inertial measurement unit.

## Acknowledgment

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. Robert D. Gregg, IV, Ph.D., holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund. Alireza Mohammadi conducted this research when he was with the Locomotor Control Systems Laboratory at the University of Texas at Dallas.

## Funding Data

National Institute of Child Health and Human Development (Grant Nos.: DP2HD080349 and R01HD094772; Funder ID: 10.13039/100000071).

Burroughs Wellcome Fund (Funder ID: 10.13039/100000861).