Fuzzy forward dynamics of distinct gait phases with a prosthetic foot

This work considers epistemic uncertainty in the form of fuzzy parameters in a multibody forward dynamics simulation of the human leg with a prosthetic foot. The thigh and shank are modelled as rigid bodies while the prosthetic foot, modelled after a carbon spring prosthesis, is represented by a predeformed geometrically exact beam model. A variational integrator is used to solve the equations of motion in the forward dynamics simulation and the Graph Follower algorithm is used to include epistemic uncertainty. Two cases are considered. Large movements are examined using a pendulum swing, similar to the swing phase during human gait. To validate the deformation of the prosthesis, a second case is examined, where the prosthesis is fixed in space and deforms under the weight of the body segments. Both cases consider a fuzzy Young’s modulus and determine the envelopes of the resulting uncertain target output. The Graph Follower algorithm was modified to substantially increase computational efficiency, enabling the propagation of fuzzy uncertainty through the complex multibody model with rigid and flexible bodies.


Introduction
The ability to walk upright is a key aspect of human locomotion. Losing this ability due to injury or illness has a significant impact on the life of the unfortunate human. For instance, Waters and Mulroy [1] examines the effect of various illnesses and gait impairments on the effort it takes to walk at different speeds and [2] summarises the causes for amputation and examines different lower limb prosthesis designs and examines future trends. Studying the effects of impairments on gait and improving the design of prosthetic feet is an ever current field of research [3], with the oldest prosthetic limbs dated around 2700 BC [4]. The main goal of lower limb prosthetic design is to mimic the biological counterpart in functionality and behaviour, while at the same time, respecting weight, size and visual appearance constraints. Apart from providing the same functionality, the prosthetic limb should also reduce the psychological effect of the prosthesis on the patient [2]. This makes prosthesis B Eduard S. Scheiterer eduard.s.scheiterer@fau.de 1 Institute of Applied Dynamics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Immerwahrstrasse 1, 90158 Erlangen, Germany design a complex and interdisciplinary endeavour. As with most engineering design processes, accurate simulation provides a powerful tool in this design process.
To accurately simulate a given system, it is necessary for the model to represent the key aspects of the system. This also requires precise knowledge of the parameters used to describe the model. However, precise determination of the model's parameters requires effort and may not always be possible. How much effort is put into the model and the precise determination of parameters, of course, depends on the goal of the simulation. This requires for a balance to be found, between accuracy and available computing power. In case of the human leg with a prosthetic foot during gait, the model has to represent the movement of the different body parts while also capturing the deformation energy stored in the prosthetic foot. In this work, thigh and shank are modelled as rigid bodies in the director formulation, see [5], while the prosthetic foot is modelled as a predeformed geometrically exact beam, as seen in [6]. However, most simulations rely on deterministic inputs which cannot always be obtained or require enormous measurement effort to approximate [7][8][9][10]. This means that when simulating, uncertainties have to be taken into account. These uncertainties can be subdivided into model and data uncertainty and can further be classified by the cause for the uncertainty. If the uncertainty is caused by incompleteness or impreciseness, it is classified as epistemic uncertainty and can be represented with fuzzy approaches. If the uncertainty is caused by variance it is classified as aleatoric uncertainty and can be represented with stochastic approaches [11]. To account for epistemic uncertainty, the Graph Follower algorithm developed in [12] is used in this work. In [12][13][14] the Graph Follower algorithm is applied to a multi-story frame with three stories. This system, a model commonly used in structural analysis, see for instance [15], has three degrees of freedom. In this work, the human leg is modelled as a multibody system with a geometrically exact beam prosthetic foot which results in a system with many degrees of freedom. This increase in degrees of freedom and the addition of predeformed geometrically exact beam theory into the multibody system result in a highly non-linear and complex system. To compensate for this increase in the complexity and lower the computational cost of propagating fuzzy uncertainty through this system, modifications to the Graph Follower algorithm are proposed to greatly increase its efficiency.
The main contribution of this work are the efficiency increasing modifications to the Graph Follower algorithm, which enable its use to propagate epistemic uncertainty through much more complex systems. This is then demonstrated for a highly non-linear forward dynamics simulation of a complex multibody model with rigid and flexible bodies. This leads to the following structure of this work. The multibody model with rigid and flexible bodies of the human leg and a prosthetic foot is detailed before summarising the Graph Follower algorithm. Then, the modifications made to the Graph Follower algorithm to substantially increase its efficiency such that it can be applied to the complex system of the human leg with a geometrically exact beam prosthetic foot model are described, followed by the examination of two configurations of the model and the discussion of the results.

Modelling the human leg with a prosthetic foot
This work focuses on the propagation of epistemic parameter uncertainty in a forward dynamics simulation of the human leg with a prosthetic carbon spring foot. To this end, a suitable model of the human leg with a deformable prosthetic foot is developed first. This model has to allow for the leg's movement and has to represent the deformation energy of the prosthesis as accurately as possible while not incurring a too high computational cost. The leg is modelled as a kinematic chain, i.e. a series of rigid bodies (thigh, shank) and flexible bodies (prosthetic foot) connected by joints using the director formulation from [16].

The leg model
The leg consists of the thigh, the shank and the prosthetic foot connected by joints. The thigh and shank are assumed to be rigid bodies while the prosthesis can deform. The hip is modelled as a spherical joint between the thigh and a fixed reference frame. The knee joint is reduced to a single degree of freedom, the rotation about one axis, and is thus modelled as a revolute joint allowing for the primary movement occurring during normal gait [17] in the sagittal plane. A rigid body, as well as a cross-section of the beam are represented by a placement ϕ ∈ R 3 and an orthonormal director triad d 1 , d 2 , d 3 ∈ R 3 , see the middle of Fig. 1.
Orthonormaltiy of the directors is enforced via so called internal constraints, representing the rigidity of the bodies or of the cross-sections respectively. The interconnection between the rigid and flexible bodies and with the fixed reference frame gives rise to external constraints representing the joints. Joint constraint formulations are explained in detail in [5,16]. In this work, parameters based on [18] are used for mass, inertia and size of the thigh and shank and are summarised in Table 1. The prosthetic foot is modelled after the Össur Vari-Flex ® carbon spring prosthesis [19]. The prosthetic foot is rigidly attached to the shank, mimicking a passive prosthetic foot. This means that there is no actuator inputting energy into the prosthetic foot during gait, opposed to for instance the "Empower Prothesenfuss" from [20], which actively controls the prosthesis during gait and can reduce metabolic cost for the wearer [2,21]. However, such an active prosthesis requires external energy. Here, a passive prosthetic foot is examined, which only stores and releases energy provided by the patient.

The carbon spring prosthesis model
In this work, the passive prosthetic foot Össur Vari-Flex ® from [19], shown in the top of Fig. 1, is modelled in [22,23] with geometrically exact beam theory from [6,24]. This formulation allows for accurate representation of complex strain states in the prosthesis while still keeping the computational demand relatively low compared to full 3D finite element models. This is crucial for the fuzzy forward dynamics simulation, since the model has to accurately represent the prosthesis' deformation energy while still keeping the computational cost low enough to make the many forward dynamic evaluations required in the uncertainty propagation algorithm feasible. The location of the centerline of the beam is denoted by ϕ(s, t) ∈ R 3 , where s runs along the beam's length and t denotes time. The orientation of a given cross-section is represented by an orthonormal director triad d 1 (s, t), d 2 (s, t), d 3 (s, t) ∈ R 3 . Due to the orthonormality, the cross-sections of the beam are modelled as rigid sections, omitting warping and distortion effects in the   [25,26] which describes an ideally elastic material behaviour.
The material parameters and geometry of the beam are included in the two matrices D Γ = diag(G A, G A, E A) and D K = diag(E I 1 , E I 2 , G J ) consisting of the Young's modulus E, the shear modulus G and geometry specific parameters e.g. the cross-section area A, the area moments of inertia I 1 , I 2 and the polar moment of inertia J . The values derived via reverse engineering in [23] are summarised in Table 2.
The strain measures Γ (q) and K (q) quantify shear, elongation, flexion and torsion with where i, j, k = 1, 2, 3. Figure 1 shows the prosthesis in the top, a geometrically exact beam with the centerline and directors in the middle, while the resulting model of the prosthesis from [22,23] in an undeformed and a deformed state, using the predeformed geometrically exact beam formulation for the prosthesis geometry, is shown in the bottom. To complete the model, viscous damping is included in the beam model. Following [25][26][27], the effective material constitutive equations are extended by a viscous term, proportional to the material strain ratesΓ i ,K i , with viscosity parameters for extension η E and shear η G . Similar to the energy function above, the Kelvin-Voigt viscous dissipative potential V is defined as combine material and geometric parameters, as for instance in [26]. Since no information on the viscous damping parameters η E , and η G for the prosthesis are currently available, this work uses η E = η G = 1.0, to counteract high frequency vibrations inside the beam. This continuous formulation of the beam can be discretised using one-dimensional finite elements, see for instance [28,29], and solved numerically. In this work, the prosthesis is discretised with 20 nodes resulting in 19 one dimensional finite elements, striking a balance between accuracy and computational cost. The choice of 20 spatial nodes is based on a convergence study building on the work detailed in [30] for the conservative case. The viscous damping introduced to the model is included in the convergence study here and the results are shown in Fig. 2. The prosthesis' highest total deformation energy W int during a 1 second simulation in a gravitational field is used as a comparison value when varying the number of nodes n used to discretise the prosthesis. Figure 2 shows the convergence of the internal deformation energy error n relative to a reference solution (computed on a very fine grid of n max = 250 nodes), when increasing the number of nodes used to discretise the prosthesis. The error is calculated for each considered number of nodes n.
Increasing the number of nodes leads to an increase in the forward dynamics computation time. Due to the impact the computation time of the forward dynamics has on the fuzzy forward dynamics simulation, 20 nodes are chosen as a balance between accuracy and computational cost in this work. This examination is based on the forward dynamics simulation, detailed in the next section. Together with the two rigid bodies, the choice of 20 nodes for the prosthesis' discretisation leads to a total of 264 degrees of freedom for the model, constrained by 132 internal constraints and 14 exter- nal constraints (three for the hip, five for the knee and six for the ankle joint).

Forward dynamics simulation
To examine the behaviour of the leg with the prosthesis before including parametric uncertainty, the model is placed in a gravitational field orientated in the negative z-direction and held in place at the hip with a spherical joint. The resulting motion of the simulation is a swinging motion of the human leg model with prosthetic foot, similar to the swing phase during the human gait cycle. A variational integrator is used to approximate the dynamics, derived from the discrete action principle for the discretised augmented Lagrangian, see [5]. All derivatives and gradients in this work are calculated automatically with Casadi [31]. For the simulation shown in Fig. 3, a timestep of 0.001 s is used. Furthermore, the Young's modulus of the prosthesis is reduced and it's density is increased, compared to carbon fibre laminate, in order to emphasise its deformation due to inertia. However, the simulation also runs stably with the original material parameters. As can be seen in Fig. 4 the energy is exchanged between the potential, kinetic and internal deformation energy. Due to the variational integrator, which is structure preserving, the total energy of the system is not artificially reduced or increased by numerical integration. The viscous damping is present, however, it's influence is so small that it is not visible in the energy evolution plot. dynamics simulation of the human leg with a prosthetic foot. A modified version of the Graph Follower algorithm developed in [12] is used to compute the uncertainty propagation through the forward dynamics simulation. Before the modifications to the Graph Follower algorithm performed in this work can be described, a basic understanding of the Graph Follower algorithm and epistemic uncertainty in the form of fuzzy numbers is necessary. A fuzzy numberp is characterised by an interval that defines its possible values and its membership function μ(p), see for instance [32].

Summary of the Graph Follower algorithm
Uncertain parametersp modelled by a triangular fuzzy number, as required for the Graph Follower algorithm, are defined by an intervalp ∈ [ p min , p max ] and their associated membership function μ p (p) ∈ [0, 1]. Only one value in the interval has the function value 1, shown in the middle of Fig. 5, leading to a convex fuzzy number. To get a triangular fuzzy number, the membership function μ p is assumed to be linear from the extremes of the interval with the value μ p ( p min ) = μ p ( p max ) = 0 to the deterministic parameter μ p ( p det ) = 1, resulting in a triangular shape for the membership function. Propagating this uncertain parameter through a forward dynamics simulation has the goal to calculate both the resulting interval of a target output quantity f d (p) as well as the corresponding membership function μ f at any time, resulting in the fuzzy target outputf d (p). This results in an optimisation problem to find the extremes of the target output given the uncertain input interval. The target output f d represents a scalar quantity of interest that is calculated from the models trajectory q(p, t) and the uncertain parametersp. The uncertain parametersp either directly affect the target output or indirectly affect it, by only influencing the forward dynamics q(p, t).
To compute this optimisation on a computer, a discretisation of the problem is necessary. Using α-level cuts, it is possible to discretise the membership function and intervals, shown in the middle of Fig. 5 for the the α-level cuts α = 0.0 and α = 0.5. For α = 1.0 the interval of the fuzzy parameters reduce to a single parameter, p det . For more details on this discretisation see for instance [7,8]. The Graph Follower algorithm computes the target output envelopes using α-level optimisation. This means the intervals of the target output function are calculated with an optimisation based on the input's fuzzy parameter interval, while the membership function is inferred from the chosen α-level of the uncertain input parameter, see [8]. This can be done for a singular uncertain parameter or multiple, however, depending on the interaction of these parameters, further unintentional uncertainty may be introduced [8]. The method of propagating fuzzy uncertainty with α-level optimisation is based on the extension principle for fuzzy numbers [32] and Nguyen's note on the extension principle [33] and is detailed in e.g. [7,8] and only summarised briefly here. The goal of propagating a fuzzy uncertain parameter through a model is to calculate the fuzzy target outputf d , which means both the numerical values f d as well as the associated membership function values μ f . These can be calculated with α-level optimisation which discretises the input fuzzy numbers with α-level cuts and then performs an optimisation to find the exterme values of the target output f d,α k for the given α-level cut α k . When the two extreme values of the target output at a given timestep are known for a given α-level cut, two points of the target output membership function μ f are known since the values are inherited from the membership function value of the α-level cut due to Nguyen's principle. The target output fuzzy numberf d can be computed by solving the following optimisation problem for each timestep and multiple α-level cuts as described for instance in [12][13][14].
Solving this optimisation problem includes the evaluation of the deterministic mapping of the target function f d and thus for this work the forward dynamics simulation. The Graph Follower algorithm combines the α-level optimisation with an approximation of the forward dynamics simulation via Taylor expansion and post-processing steps to efficiently calculate the target output envelope and the target output's membership function over time, resulting in an approximation of the fuzzy target output. The process of forward propagation of epistemic uncertainty to calculate a target output envelope is visualised in Fig. 5.
The Graph Follower algorithm from [12] employs two main methods to make α-level optimisation feasible with forward dynamics. The first is the linearisation of the forward dynamics with a Taylor expansion with respect to the last known optimal parameter, allowing for quick computations of the target output at the cost of accuracy. The second method, is the combination of storage and post-processing. When a new parameter is found by the optimisation, the forward dynamics are calculated with it and the target output is computed for the entire simulation time. Both the parameter as well as its associated target output evolution are stored in a library within one α-level. This library of previously computed parameters and target output trajectories is used to chose a parameter for the linearisation of the forward dynamics based on the currently known extreme target output trajectory. The parameter p that is associated with the currently known extreme target output trajectory is used for the linearisation to approximate the forward dynamics and allow an efficient optimisation. Furthermore, at the end of an α-level optimisation for a given α-level, all target output trajectories are compared and only the extreme values are transferred to the target output envelope. This is repeated with several α-levels, to define the membership function value for the envelopes.

Fuzzy parameters of the prosthesis model
The prosthesis' Young's modulus is assumed to be subject to uncertainty. Alternative uncertain parameters could be the body segment length, joint locations, other material parameters for the prosthesis or any combination of these parameters. However, since this work's goal is to demonstrate the propagation of epistemic uncertainty through a complex nonlinear model and to allow for an easy interpretation of the results, only the Young's modulus is assumed to be uncertain here before examining multiple uncertain parameters. The prosthesis' Young's modulus deterministic value E det is subjected to epistemic uncertainty by modelling it as a triangular fuzzy number as [0.9, 1.0, 1.1]· E det , shown exemplary in the middle of Fig. 5. The fuzzy Young's modulus is discretised with three α-levels: [1.0, 0.5, 0.0], where α = 1.0 is expected to yield the same results as the deterministic solution for triangular input fuzzy numbers, while α = 0.0 allows for the largest variance in the parameter. It should be noted, that the parameter values in the interval have to lead to solutions for the forward dynamics simulation before the Graph Follower algorithm can be applied.

Speeding up the Graph Follower algorithm
This section details the modifications made to the Graph Follower algorithm in order to achieve the shown results. These modifications are necessary to keep it computationally feasible for the complexity of the analysed system.

Reduction of evaluation points
To accomplish this, the timestep of the forward dynamics and the Graph Follower algorithm are separated. To compute stably, the forward dynamics of the leg with a geometrically exact prosthetic foot requires a timestep of 0.001 s. Using the same timestep for the Graph Follower algorithm for a simulation time of 2.0 s would result in 2000 α-level optimisations for each discrete α-level envelope, meaning once for the upper and once for the lower boundary, for every examined α-level. Additionally, for every new optimal parameter found by the Graph Follower algorithm, a new forward dynamics simulation is necessary to calculate the associated target output. In the worst case, this would lead to n α−levels * n timesteps * 2 forward dynamics simulations, which is computationally expensive.
To reduce this, the α-level optimisation is only performed every 100 steps of the forward dynamics simulation. This leads to greatly reduced computational cost. The resulting envelopes are still calculated with the forward dynamics timestep of 0.001 s, resulting in smooth trajectories. The number of required optimisations is greatly reduced at the cost of possibly not finding all extreme trajectories. Figure 6 shows the resulting envelopes for two different optimisation timesteps to compare the effects of this reduced optimisation timestep. Table 3 lists the maximum difference between the envelopes. Both simulations were performed using the parameters for the swing configuration, described in detail in the next section. The difference between the two compared optimisation timesteps is negligible for this case.

Tolerance distinguishing optimal parameters
The Graph Follower algorithm maximises or minimises the target output function within the parameter interval. If an optimal parameter is found, it is stored in a library and used for the approximation of the forward dynamics in the optimisation. The last parameter associated with the extreme target output for the current examined timestep is used to linearise the forward dynamics during the optimisation. This parameter library is also used to check whether the new optimal parameter is indeed new, by comparing it to the parameter used for linearisation. A new forward dynamics simulation is only performed if the difference of the optimised parameter to known parameters exceeds a percentual tolerance. This tolerance can be estimated by a sensitivity analysis of the forward dynamics with respect to the uncertain input parameter. To further reduce the computational cost, the forward dynamics simulation is recalculated and the new target output is calculated from the resulting trajectory only if the found optimal parameter is sufficiently different from all previously calculated parameters. Otherwise, the stored values are used. Ensuring the novelty of a parameter within an α-level greatly reduces the computational cost of the Graph Follower algorithm. Of course, this tolerance directly affects the accuracy of the target output envelopes. This means it has to be dimensioned to fit to the examined problem. Also, multiple uncertain parameters make it more difficult to apply, since multiple parameters do not affect the target output without interdependencies. In this work, new forward dynamics calculations are performed if the relative difference p diff,rel between the current fuzzy parameter and the closest stored parameter exceeds p tol of the deterministic value. In the case of a simulation with multiple fuzzy parameters, a new simulation is performed if any individual parameter difference exceeds p tol of the respective deterministic parameter.
The tolerance value p tol = 0.01% was determined experimentally. For more precision in the choice of this tolerance, a sensitivity analysis of the target output with respect to the fuzzy parameter can be performed.

Limitations of the Graph Follower algorithm
Before examining the scenarios simulated in this work, the limitations of the Graph Follower algorithm are briefly summarised. The goal of propagating fuzzy input parameters through a model, is to calculate the target output fuzzy number. This includes the interval of the target output f d and the associated membership values μ f . As detailed above, the Graph Follower algorithm uses α-level optimisation and a linearisation of the forward dynamics with respect to the fuzzy input parameter to efficiently calculate the target output envelopes for multiple α-levels, where the membership function is derived from the current α-level, resulting in an approximation of the target output fuzzy number. The α-level optimisation with forward dynamics imposes limitations on the accuracy of the resulting envelopes. To work precisely, the α-level optimisation requires a global optimiser and precise calculation of the target output, see [7]. However, global optimisation is computationally expensive and the accurate calculations of the forward dynamics required for the exact target output calculation are also computationally expensive with such a complex model. Thus, to be feasible, the Graph Follower algorithm only approximates global optimisation, by using Matlab's MultiStart function with ten starting points. This significantly increases the result quality with an acceptable increase in computation time, compared to using fmincon with just one starting point. The linearisation of the forward dynamics with respect to the fuzzy parameter introduces more inaccuracy. Both of these limitations can lead to an overlap of α-levels. The modifications necessary for the Graph Follower algorithm to work with the complex model of the human leg with a deformable prosthesis, also require a trade-off between accuracy and computational cost. The separation of the timestep of the forward dynamics simulation and the α-level optimisation greatly reduces the amount of necessary optimisations but also reduces the chances to find extremising parameters. Also, the tolerance when comparing new parameters to stored parameters can reduce accuracy if it is chosen to large, however, it greatly reduces the computational cost of the modified Graph Follower algorithm. Nonetheless, the modified Graph Follower algorithm, as is used here, produces an approximation of the target output fuzzy number for the complex multibody system of the human leg with a predeformed geometrically exact beam prosthetic foot. Furthermore, these limitations can be reduced or even become negligible by increasing the available computing power. Also, if the precise membership function value of the target output is not required, the trajectories of the target output of the individual α-levels can be combined to get an accurate global envelope of the target out- put, providing information on the extreme values of the target output based on the fuzzy input parameter. Thus, even with limitations, the modified Graph Follower algorithm provides valuable information about how epistemic model parameters affect a desired target output for a highly complex multibody system with flexible and rigid bodies, allowing for a quicker examination of the uncertainties effects than simulating many different sets of parameters.

Simulation scenarios
In the two exemplary scenarios analysed in this work, the leg moves under its weight due to being placed in a gravitational field oriented in the negative z-direction. The first examined scenario is a pendulum swing of the multibody system, where the prosthesis deforms due to inertia, as was already used to check the models implementation with forward dynamics, shown as a timelapse in Fig. 3. The pendulum movement is similar to the swing phase occurring during human gait (for more information on the biomechanics of the human leg and the human gait cycle see for instance [17]). The parameters used for the simulation are summarised in Table 4.
The second scenario simulates the leg in a squatting position, that is the thigh is horizontal with the shank vertical. The prosthetic foot is held in place at two nodes. This is done by constraining two nodes to remain in their initial position. However, rotation of the cross-sections at these nodes is permitted. The prosthesis deforms under the leg's weight periodically. Figure 7 shows the initial configuration of the leg alongside a deformed configuration of the resulting movement. The deformation is due to the mass of the leg in a gravitational field. To emulate the process of realistic load on the prosthesis, the mass of the shank is increased to represent the influence of body weight, while the prosthesis uses the parameters derived with reverse engineering. The parameters used for the simulation are summarised in Table 5.

Emulation of the swing phase
In the simulation of the leg's swing movement, the chosen target output f d is the total deformation energy stored in the prosthesis at a certain time t n .
The total deformation energy stored and released during gait in the prosthetic foot is related to the burden of gait for a patient, see for instance [1]. This means it is directly related to the walking comfort of the patient with the prosthetic foot and is therefore an important measure for the prosthesis' design and it's acceptance by the patient. Predicting the prosthesis' behaviour with respect to the target output is important during the design process. However, predicting the movement of a double pendulum is difficult, given its chaotic nature. Adding a flexible body with uncertain stiffness to the end of the double pendulum further complicates the prediction. Including fuzzy uncertainty and a target output that depends on both the behaviour of the multibody model and the fuzzy uncertainty increases the need for algorithmic computation of the desired target output. Figure 8 shows the calculated trajectories for the minimal (0.9 E det ), deterministic (E det ) and maximum (1.1 E det ) values of the prosthesis' Young's modulus. The Graph Follower algorithm is used to calculate the target output envelopes. Analysing the uncertainty propagation allows for an efficient calculation of the correlation between the parametric input uncertainty, in this case the Young's modulus, and the target output, the stored internal deformation energy of the prosthesis. The results are shown in Fig. 9. The target output envelopes for the different alpha levels are shown alongside the deterministic simulation's result. The α = 1-level envelope matches the deterministic solution, as expected. The total deformation energy of the prosthesis oscillates which is represented well by the variational integrator used in the forward dynamics which avoids artificial energy loss due to numerical integration. As shown in the examination of the forward dynamics of this model, the energy drain of the viscous damping is negligible for the chosen parameters. While the input uncertainty is triangular and symmetric, the output envelopes shown in Fig. 9 do not inherit these properties. The limitations of the Graph Follower algorithm can be seen at t = 1.4s when the upper envelope from α = 0.0-level intersects the envelope from α = 0.5-level.
However, when considering the entire evolution of the envelopes, the inaccuracy is very localised to parts of the tar-

Emulation of the prosthesis during a squatting exercise
The second target output analyses the maximal local deformation energy over all finite elements during the simulation at a certain time t n . This is related to the structural integrity of the prosthesis during a load scenario.
The calculations use the squatting pose configuration of the model and more realistic parameters for the material of the prosthesis. Simulations were performed with realistic parameters. However, due to the prosthesis being lightweight and stiff the Courant-Friedrichs-Lewy-condition [34] requires the simulation to run with a very small timestep which leads to very computationally intensive forward dynamics simulations. The load applied to the prosthesis emulates the load of the body during a squat exercise. The envelopes of the target function are shown in Fig. 11. As expected, the prosthesis acts as a spring under the load of the shank. The leg moves up and down periodically, while the uncertain Young's modulus affects both the amplitude and the period of the oscillation. Since the amplitude and frequency of the oscillating movement of the leg is directly related to the prosthesis' deformation, the oscillating movement can also be seen in the resulting envelopes for the chosen target function. Artefacts from the simplifications necessary for the Graph Follower algorithm to run with the high number of degrees of freedom of the model are present, but can be reduced by increasing computational power.

Multiple fuzzy parameters
After demonstrating the application of the Graph Follower algorithm on the complex multibody system of the human leg with a geometrically exact beam prosthesis with a single fuzzy uncertain parameter, the next step is to examine the structure with multiple uncertain parameters. For this, two aspects are of interest. How does the method perform with regards to the increased number of optimisation parameters and which parameters can be affected with uncertainty.

Numerical efficiency of Graph Follower algorithm with multiple uncertain parameters
Here a brief evaluation of the Graph Follower algorithm with respect to its numerical efficiency is discussed. A comprehensive evaluation of the Graph Follower algorithm and comparison to other methods can be found in [12]. There, it is shown that the Graph Follower algorithm is suitable for propagating multiple fuzzy parameters since its computing time scales much less with the amount of uncertain parameters than the benchmark algorithm Qua.Si.III. To compare The relative change is displayed on the outside, with the increased complexity runtime divided by the simpler simulation runtime the modified Graph Follower algorithm used in this work, the runtime for four simulations are compared, one with a single fuzzy parameter and one with 11 fuzzy parameters, each for 1 and 2 s with 1000 and 2000 timsteps respectively.
The resulting runtimes are shown in Table 6. This shows, that an increase in the amount of uncertain parameters does not affect the computation time of the simulation strongly. This also means, that focus should be put on reducing the computational time of the forward dynamics simulation, to reduce the computational cost when using the Graph Follower algorithm to propagate fuzzy uncertainty.

Choice of fuzzy parameters
The main limiting factor when choosing which parameters to model as uncertain, is the requirement, that the parameters are independent of each other and without interaction [7]. Interdependence or interaction between fuzzy numbers can negatively affect the α-level optimisation by introducing unwanted uncertainty. Possible uncertain parameters are the leg's segment length, the segment's mass and moments of inertia, the segment's initial orientation and the prosthesis material parameters, namely Young's modulus, density and Poisson's ratio. Here, affecting the segment length or initial position with uncertainty would lead to fuzzy configurations of the leg. However, the configuration is not independent, since it has to satisfy the internal constraints of the rigid bodies of the segments and the joint constraints and can therefore not be affected by uncertainty without further investigation. The prosthesis' Young's modulus, Poisson's ratio and density parameters can safely be assumed to be independent and are affected with uncertainty. For the leg, its mass is affected with uncertainty, including the moments of inertia. These values are not independent of each other, however the exact distribution of the mass along the leg is unknown for an individual and difficult to measure in vivo, which allows for an independent treatment in this work. This leads to a simulation with a total of 11 uncertain parameters. All fuzzy parameters are discretised with the same α-level cuts and for comparability also have the same parameters and interval of p det ± 0.1 p det as the squatting simulation.

Results
As in the previous simulation the movement is an oscillation. With the increased number of fuzzy parameters, the period and magnitude of the oscillations differ more than in the simulation with only the fuzzy Young's modulus leading to larger envelopes, shown in Fig. 12. Due to the multiple uncertain parameters, visualisation of the parameter distribution and their correlation to the envelopes is no longer trivial. Figure 13 visualises the parameter distribution for the α = 0 level with the individual trajectories calculated from the corresponding parameter set and the resulting upper and lower envelope. The left side shows the relative difference with respect to the upper (green) and lower (red) bounds for the parameters with the deterministic parameters in the middle (blue). As can be seen, all of the parameters vary and lead to individually extreme target outputs for a given timestep, shown on the right of Fig. 13. It can also be seen, that for instance the Poisson's ratio (ν) leading to the envelope varies a lot less than for instance the mass of the thigh segment (L 1,M ), giving an indication of which parameters have a strong effect on the target output. Also of note, the moments of inertia for the z-direction (L i,T 3 ) do not have an effect on the forward dynamics, since the movement has no rotation around the z-axis. This means, that these parameters can vary without affecting the target output. This results in most trajectories of the simulation shown in Fig. 13 being based on the deterministic value.

Conclusions
In this work, the forward propagation of epistemic uncertainty on the parameter level is investigated for two cases of a human leg with a geometrically exact beam prosthetic foot in a fuzzy forward dynamics simulation. The leg is modelled as a multibody system with rigid and flexible bodies connected by joints. The prosthesis' Young's modulus, subject to fuzzy uncertainty, is propagated forward through the model and multiple target outputs are evaluated. Multiple fuzzy input parameters are also successfully propagated through the forward dynamics model. Remarkably, increasing the number of uncertain parameters does not significantly increase affect the computation cost. The main contribution of this work is the adaptation of the Graph Follower algorithm to a highly non-linear and complex multibody dynamics setting with rigid and flexible bodies. To achieve this, the efficiency of the Graph Follower algorithm is increased substantially. The computational cost of the algorithm is reduced by coarsening the evaluation grid and adding a tolerance to identify new extreme parameters during the optimisation. This slightly reduces the accuracy in the resulting envelopes, but allows for the use of the Graph Follower algorithm to propagate an uncertain Young's modulus through to the leg model with a prosthetic foot. This provides an approximation of the fuzzy number of the target output and it's evolution during the simulation time.
The forward propagation of fuzzy uncertainty through a multibody model with rigid and flexible bodies, can be extended in multiple directions. The model can be expanded from a forward dynamics model, mimicking gait phases and loads, to an optimal control simulation of the human gait. Furthermore, the uncertainty aspect can be extended. As briefly mentioned in the introduction, there are multiple ways to model uncertainty, divided into aleatoric and epistemic uncertainty, or the combination of the two, polymorphic uncertainty. Developing the forward propagation of polymorphic uncertainty is planned through either forward dynamics and also optimal control.