Hyfydy User Manual

Introduction

Hyfydy is software for high-performance musculoskeletal simulation, with a focus on biomechanics research and predictive simulation of human and animal motion. It supports bodies, joints, contacts, musculotendon actuators, and torque actuators. Some of its key features include:

• High performance. Hyfydy simulations run approximately 100x faster than OpenSim¹, using the same muscle and contact models*. Hyfydy is also fast enough for use in real-time applications such as games.

• Stable integration. Efficient error-controlled variable-step integration methods ensure the simulation always runs stable at a user-specified accuracy level, without sacrificing performance.

• Accurate models. Hyfydy contains research-grade state-of-the-art models for musculotendon dynamics and contact forces.

• Force-based constraints. Joint constraints in Hyfydy are modeled through forces that mimic the effect of cartilage and ligaments. This allows for customizable joint stiffness and damping, and closed kinematic chains. Hyfydy is optimized to run at high simulation frequencies (>3000Hz) to efficiently handle stiff joint and contact constraint forces.

• Intuitive model format for easy model building and editing. A tool for converting OpenSim models is included (only a subset of OpenSim features is supported).

Hyfydy is currently available as a plugin for the open-source SCONE² simulation software. For more information on SCONE, please visit the SCONE website.

* Performance comparison between Hyfydy and OpenSim was based on a 10 second gait simulation using a reflex-based controller³ with hill-type musculotendon dynamics⁴ and non-linear contact forces⁵. Control parameters were optimized in SCONE² using 1000 generations of CMA-ES⁶ on a 4-core i7 Intel CPU.

Citation

If you use Hyfydy in a research project, please use the following citation:

@misc{Geijtenbeek2021Hyfydy,
   author = {Geijtenbeek, Thomas},
   title = {The {Hyfydy} Simulation Software},
   year = {2021},
   month = {11},
   url = {https://hyfydy.com},
   note = {\url{https://hyfydy.com}}
}

Model

File Format

Hyfydy employs an intuitive text-based model format (.hfd), which is designed to be easily edited or constructed by hand using a text editor. It consists of key-value pairs in the form key = value, curly braces { ... } for groups and square brackets [ ... ] for arrays. Comments are supported through via # (one-line) or /* */ (multi-line).

Example

xxxxxxxxxx
model {
    name = my_model 
    body {
        name = my_body
        mass = 3
        inertia = [ 1.5 1.5 1.5 ]
        pos = [ 0 1 0 ] # single-line comment
    }
    /* multi line
    comment */
}

Data Types

The .hfd file format uses several reoccurring data types, which are explained below.

string

Strings are used to identify components. Quotes are not required, unless an identifier contains whitespace or special characters ({, }, [, ], \ or =):

xxxxxxxxxx
name = my_name
another_name = "Name with special {characters}"

vector3

3D vectors are used to represent position, direction and linear / angular velocity. Values can be entered as an array, or using x, y, or z components:

xxxxxxxxxx
position = [ 0 1 0 ]
velocity { x = 0 y = 1 z = 0 }
earth_gravity { y = -9.81 }

quaternion

Quaternions are used to represent rotations and orientations. They can be initialized using their w, x, y, z components directly, but also via an x-y-z Euler angle (recommended):

xxxxxxxxxx
ori1 { z = 90 } # 90 degree rotation around z
ori2 = [ 45 0 90 ] # 45 degrees around x, then 90 degrees around z
ori3 { w = 1 x = 0 y = 0 z = 0 }

range

Ranges are used to describe any kind boundary, and are written as two numbers separated by two dots (..):

xxxxxxxxxx
joint_limits { x = -10..10 y = 0..0 z = -45..90 }
more_joint_limits = [ -10..10 0..0 -45..90 ]
dof_range = 0..90

Core Components

model

The top-level model component is mostly a container for other components. A typical model looks like this:

xxxxxxxxxx
model {
    name = my_model
    gravity = [ 0 -9.81 0 ]
    material { ... }
    model_options { ... }
    
    body { ... }
    joint { ... }
    geometry { ... }
    ...
}

A model can contain the following properties:

Other components, including material and model_options, are described below.

body

Bodies are specified using the body component. They contain mass properties, position/orientation as well as linear and angular velocity. The body frame-of-reference always has its origin located at the center-of-mass and its axes must be aligned with the principle axes of inertia of the body.

A body component can contain the following properties:

* A body must contain either mass and inertia, or density and shape in order to have valid mass properties.

Example

xxxxxxxxxx
body {
    name = my_body
    mass = 3
    inertia = [ 1.5 1.5 1.5 ]
    pos = [ 0 1 0 ]
    ori = [ 0 0 45 ]
}

When specifying a density and shape and property, the mass and inertia are calculated automatically. Supported shape types are:

• sphere (with radius parameter)

• cylinder and capsule (with radius and height parameters)

• box (with half_dim parameter)

Example

xxxxxxxxxx
body {
    name = my_cube
    density = 1000
    shape {
        type = box
        half_dim = [ 0.1 0.1 0.1 ]
    }
    pos = [ 0 1 0 ]
}

Note: a body component does not include geometry used for contact detection and response. These are defined separately via a geometry component.

joint

Joint components constrain the motion between two bodies. They can contain the following properties:

There are no restrictions to the amount of joints a body can contain, and kinematic loops are allowed. However, adding superfluous joints will impact simulation performance.

Example

xxxxxxxxxx
joint {
    name = knee_r
    parent = femur_r
    child = tibia_r
    pos_in_parent = [ 0 -0.226 0 ]
    pos_in_child = [ 0 0.1867 0 ]
    limits { x = 0..0 y = 0..0 z = -90..0 }
}

Alternatively, a joint component can be specified inside a body component, in which case the child property can be omitted:

xxxxxxxxxx
body {
    name = tibia_r
    mass = 3.7075
    inertia { x = 0.0504 y = 0.0051 z = 0.0511 }
    joint {
        name = knee_r
        parent = femur_r
        child = tibia_r
        pos_in_parent = [ 0 -0.226 0 ]
        pos_in_child = [ 0 0.1867 0 ]
        limits { x = 0..0 y = 0..0 z = -90..0 }
    }
}

geometry

The geometry component defines the properties needed for determining contacts and subsequent contact forces. They can contain the following properties:

Supported shape types for geometry are:

• sphere (with radius parameter)

• capsule (with radius and height parameters)

• box (with half_dim parameter)

• plane (with normal parameter)

An example of a geometry:

xxxxxxxxxx
geometry {
    name = l_heel
    type = sphere
    radius = 0.03
    body = calcn_l
    pos { x = -0.085 y = -0.015 z = 0.005 }
    ori { x = 0 y = 0 z = 0 }
}

Alternatively, geometry components can be specified inside a body component, in which case the body property can be omitted. The shape of the geometry can also be used to automatically calculate the mass properties.

material

The material component describes material properties used to compute contact forces. They contain the following properties:

Note: the interpretation of the friction, stiffness and damping parameters depend on the type of contact force that is used in the simulation.

model_options

The properties defined in model_options are global defaults that apply to all components defined within the model. Their goal is to minimize duplication of common properties. The following properties can be specified within a model_options section:

Actuators

point_path_muscle

The point_path_muscle component specifies a musculotendon unit which path is defined through a series of via points. It contains the following properties:

Note: the way in which muscle force is computed depends on the muscle force that is used in the simulation.

Example

xxxxxxxxxx
point_path_muscle {
    name = hamstrings_r
    tendon_slack_length = 0.31
    optimal_fiber_length = 0.109
    max_isometric_force = 2594
    pennation_angle = 0
    path [
        { body = pelvis pos { x = -0.05526 y = -0.10257 z = 0.06944 } }
        { body = tibia_r pos { x = -0.028 y = 0.1667 z = 0.02943 } }
        { body = tibia_r pos { x = -0.021 y = 0.1467 z = 0.0343 } }
    ]
}

joint_point_path_muscle

The joint_point_path_muscle is identical to a point path muscle, with the exception that with these types of muscles, joint torques are applied instead of forces. The advantage of applying a torque is that it leads to less joint displacement, which in turn can improve performance. Applying torques instead of forces also makes the simulation more similar to reduced coordinate simulation engines, such as OpenSim. If instead accurate computation of joint displacement caused by muscle force is required, this type of component is not recommended.

The properties of the joint_point_path_muscle are identical to those of a point_path_muscle.

joint_motor

Joint motors produce a 3D joint torque based on joint angle and joint velocity. The amount of torque $\tau$$\tau$ is based on base torque ${\tau }_o$$\tau_{o}$, stiffness $k_p$$k_p$, damping $k_d$$k_d$, orientation $q$$q$, target orientation $q_t$$q_t$, angular velocity $\omega$$\omega$, and target velocity ${\omega }_t$$\omega_t$:

The notation $[{]}^{{\tau }_{max}}$$[\space]^{\tau_{max}}$ is used to indicate that the magnitude of the final torque is clamped between $[-{\tau }_{max},{\tau }_{max}]$$[-\tau_{max}, \tau_{max}]$. The function $r:{\mathbb{R}}^4\to {\mathbb{R}}^3$$r: \R^4 \rightarrow \R^3$ converts a quaternion to a 3D rotation vector, which corresponds to the 3D rotation axis scaled by rotation angle in radians.

A joint_motor component can contain the following properties:

Example

xxxxxxxxxx
joint_motor {
    joint = hip_r
    stiffness = 1000
    damping = 10
    max_torque = 100
    target_ori = [ 0 0 90 ]
    target_vel = [ 0 0 0 ]
    torque_offset = [ 0 0 0 ]
}

Modifying joint_motor targets

Joint motors can be modified during the simulation, allowing them to be part of complex control strategies with shifting targets and other properties. In SCONE, joint motor components can be altered through the script interface, via the functions:

• set_motor_target_ori( quaternion )

• set_motor_target_vel( vector3 )

• set_motor_stiffness( number )

• set_motor_damping( number)

• add_motor_torque( vector3 ).

Auxiliary Components

Auxiliary components are not part of the simulation and therefore do not affect the outcome, but can be used by external software for analysis and visualization.

mesh

Components of type mesh can be used by client applications for visualizing bodies. In SCONE, mesh components are visualized in the 3D viewer window. They can contain the following properties:

Example

xxxxxxxxxx
mesh {
    file = example.obj
    pos = [ 0 0.1 0.1 ]
    ori = [ 0 90 0 ]
}

xxxxxxxxxx
mesh {
    shape {
        type = cylinder
        height = 0.5
        radius = 0.05
    }
    pos = [ 0 0 0 ]
    ori = [ 0 0 0 ]
    color = [ 0.8 0.8 0.8 ]
}

xxxxxxxxxx
mesh {
    shape {
        type = box
        dim = [ 0.3 0.1 0.1 ]
    }
    color = [ 0.8 0.8 0.8 ]
}

dof

Components of type dof can be used to describe a specific degree-of-freedom, which can be used by a client application for control and analysis. In SCONE, dof components are used to define model coordinates, which are used for analysis and control.

In Hyfydy, all bodies contain 6 degrees-of-freedom (3 translational and 3 rotational) – specifying specific degrees-of-freedom using a dof component does not affect the simulation.

The dof component can contain the following properties:

In Hyfydy, all bodies contain 6 degrees-of-freedom – 3 translational and 3 rotational. In practice, joints limit their movement…

xxxxxxxxxx
dof {
    name = pelvis_tilt
    source = pelvis_rz
    default = 0
    range = -90..90
}

Forces

The motion of bodies is governed by forces and torques, which cause linear and angular accelerations on individual bodies. In Hyfydy, joint constraints are too enforced through explicit constraint forces. Each force is generated by a specific force component, which can be configured flexibly and at run-time, via a configuration script. It is possible to distinguish between different types of forces, including joint forces, contact forces, actuator forces and external forces.

Example

xxxxxxxxxx
composite_force {
    planar_joint_force_pnld {}
    simple_collision_detection {}
    contact_force_hunt_crossley_sb { transition_velocity = 0.1 }
    muscle_force_m2012fast { xi = 0.1 }
}

Multiple force component specified, actual forces are applied in order of specification. Each individual force component has its own properties.

Joint Forces

Joint forces hold bodies together, like ligaments and cartilage do in human and animal joints. In addition, ligaments and bony structures also limit the range of motion between bodies – these forces become active after a joint rotates beyond a specific threshold. In Hyfydy, the former set of forces is referred to as joint constraint forces, while the latter is referred to as joint limit forces. Both joint constraint and joint limit forces are produces by a single force component.

joint_force_pnld

For each joint $j$$j$, the joint_force_pnld component applies a force $F_j$$F_j$ based on joint displacement $p_j$$p_j$ and displacement velocity $v_j$$v_j$:

The force $F_j$$F_j$ is applied in opposite directions to the child and parent body of the joint, at their respective pos_in_child and pos_in_parent position. See joint for more details.

In addition to a joint constraint force, joint_force_pnld also applies a joint limit torque ${\tau }_j$$\tau_j$ based on joint angle ${\alpha }_j$$\alpha_j$ and angular velocity ${\omega }_j$$\omega_j$. While the limit torque is proportional to joint displacement angle ${\theta }_j$$\theta_j$, the limit damping torque is non-linear. This prevents damping to become active immediately after a joint angle crosses its limit. Instead, Hyfydy follows similar considerations as for the Hunt-Crossley contact restitution force⁵, and scales the limit damping by the limit displacement torque:

in which the joint displacement angle ${\theta }_j$$\theta_j$ is defined as:

The constants $k_p$$k_p$, $k_d$$k_d$, $k_{\alpha }$$k_\alpha$ and $k_{\omega }$$k_\omega$ correspond to the stiffness, damping, limit_stiffness and limit_damping parameters, which can be specified individually per joint, or via defaults based on model_options. As a result, joint_force_pnld does not contain any options for itself and can simply be added as:

xxxxxxxxxx
joint_force_pnld {}

joint_force_pd

This force is similar to joint_force_pnld, with the exception that the limit damping torque is linearly proportional to the angular velocity ${\omega }_j$$\omega_j$:

This damping component is active immediately after a joint angle crosses its limit angle, resulting in an immediate torque in the opposite direction. As a result, using joint_force_pnld is recommended instead.

Planar Joint Forces

All joint forces have planar variants, which generate forces only in the x-y plane, and torques only around the perpendicular z-axis. Planar forces greatly improve simulation performance for planar (2D) models.

The force component names of the planar joint forces are as follows:

Important: planar forces should always be used in combination with a planar integrator.

Contact Forces

Contact forces are applied when the geometries of two bodies intersect. A contact force consists of a restitution and friction component, which are computed based on the penetration depth, the relative velocity of the bodies at the point of contact, and the material properties of the contact. The computation of contact forces consists of two phases:

1. Collision detection. Geometries are checked for intersections, and contacts between geometries are recorded.

2. Collision response. Restitution and friction forces are applied to colliding bodies, based on recorded contacts.

In Hyfydy, each phase is configured using a separate force component, even though technically, the collision detection component only generates contacts and no forces.

Important: both a collision detection and collision response force must be present in order for contact forces to work.

Collision Detection

Collision detection components detect intersections between the geometries of pairs of bodies. Therefore, it is required to have a geometry attached to at least two different bodies in order for collision detection to work.

simple_collision_detection

The simple_collision_detection component detects intersections by going through all relevant pairs of objects. While this is a reasonable strategy when not many geometries are present (as is the case with typical biomechanics simulations), efficiency quickly declines when the number of bodies increases. Therefore, it is important to keep the number of geometries at a minimum when using the simple_collision_detection component.

By default, the component is configured to only detect intersections between a static plane geometry and geometries attached to bodies – and not intersections between different bodies. This behavior can be changed by setting the enable_collision_between_objects property to 1.

Important: setting enable_collision_between_objects = 1 when many geometries (>10) are present severely impacts simulation performance.

The following properties can be set for simple_collision_detection:

Example

xxxxxxxxxx
simple_collision_detection { enable_collision_between_objects = 1 }

Combining Material Properties

When a collision is detected and a contact is generated, the material properties of both geometries are combined into a new set of parameters, which in turn are used by the collision response force.

Collision Response

Collision response components compute actual contact forces, based on the contact found during the collision detection phase.

contact_force_pd

The contact_force_pd produces a linear damped spring contact restitution force, also known as the Kelvin-Voigt contact model. Given the material stiffness coefficient $k$$k$ and dissipation coefficient $c$$c$, the normal force $F_n$$F_n$ is derived from penetration depth $d$$d$ and normal velocity $v_n$$v_n$, which is the contact velocity in the direction of contact normal $\overrightarrow{n}$$\vec{n}$:

The tangent force $\overrightarrow{F_t}$$\vec{F_t}$ depends on the static friction coefficient $\mu$$\mu$, viscosity parameter $\eta$$\eta$ and tangent velocity $\overrightarrow{v_t}$$\vec{v_t}$:

The total contact force $\overrightarrow{F_c}$$\vec{F_c}$ corresponds to:

This force has no extra parameters and can be added by including:

xxxxxxxxxx
contact_force_pd { viscosity = 1000 }

contact_force_hunt_crossley

The non-linear Hunt-Crossley⁵ contact force provides a better model of the dependence of the coefficient of restitution on velocity, based on observed evidence. Given material stiffness $k$$k$, damping $c$$c$, penetration depth $d$$d$ and normal velocity $v_n$$v_n$ the contact restitution force $F_n$$F_n$ corresponds to:

The friction force $\overrightarrow{F_t}$$\vec{F_t}$ and resulting contact force $\overrightarrow{F_c}$$\vec{F_c}$ are defined as in contact_force_pd.

Contact Force Stiffness

Traditionally, the Hunt-Crossley stiffness depends on the radius of the contact sphere. In Hyfydy, the sphere radius is directly incorporated into the stiffness constant. The stiffness constant is updated automatically as part of the conversion from OpenSim to Hyfydy. Given contact sphere radius $r$$r$ and OpenSim stiffness $k_{osim}$$k_{osim}$, the Hyfydy stiffness $k$$k$ corresponds to:

contact_force_hunt_crossley_sb

Similar to contact_force_hunt_crossley, but with a friction model attributed to Michael Hollars and used by Simbody and OpenSim. This force introduces the transition_velocity property, which is described here. Lowering the transition velocity increases accuracy, at the cost of simulation performance as a result of requiring smaller integration time steps.

Example

xxxxxxxxxx
contact_force_hunt_crossley_sb {
    transition_velocity = 0.2
}

Actuator Forces

muscle_force_m2012fast

This is an optimized implementation of the Hill-type muscle model described by Millard et al.⁴, which includes a passive damping term that allows velocity to be determined even when the muscle is deactivated. This is the recommended muscle model to be used in Hyfydy.

In order to increase performance, the Hyfydy implementation of the Millard Equilibrium Muscle Model ⁴ differs from the OpenSim implementation in two significant ways:

1. The curves that describe the force-length and force-velocity relations, as well as the curves for passive tendon and muscle forces, are defined through polynomials instead of splines. The resulting curves in Hyfydy therefore differ slightly from the curves used in the OpenSim implementation.

2. The muscle damping forces are approximated using an explicit method, instead of using an iterative method as suggested in the publication and used in the OpenSim implementation. While this approach greatly increases performance, it causes the resulting muscle damping forces in Hyfydy to differ slightly from the damping the forces produced in the OpenSim implementation.

The muscle_force_m2012fast can include the following properties:

muscle_force_gh2010

This is an implementation of the Hill-type muscle model described by Geyer & Herr³. It includes a passive spring that prevents the muscle from shortening after a specific length threshold. The force-velocity relationship is undefined at zero activation, as a result activation muscle be kept above a threshold (e.g. > 0.01) to ensure stability.

muscle_force_tj2003

This is an implementation of the muscle model described in a document authored by Chand T. John, which describes the implementation in OpenSim of the muscle model published by Thelen⁷. Despite its popularity, this model is best avoided, because its force-velocity relationship is poorly defined at low activation and relies heavily on its ad-hoc extrapolation of the force-velocity curve. We recommend using the muscle_force_m2012fast model instead.

joint_motor_force

The joint_motor_force component produces joint torques defined by joint_motor components, and needs to be included for models that use them. They contain no additional settings.

xxxxxxxxxx
joint_motor_force {}

External Forces

External forces include gravity and perturbation. These forces are applied to automatically when defined in the model.

gravity

Gravity is applied automatically to all non-static bodies in the model. The magnitude and direction is determined by the gravity property defined inside the model. To disable gravity, simply add gravity = [ 0 0 0 ].

perturbation

Perturbation forces are applied automatically based on perturbation components defined inside the model. Since perturbations usually change over time, perturbation components are typically not defined in the model itself, but are added by a client application.

In SCONE, perturbations and perturbation forces are enabled automatically when after specifying enable_external_forces = 1 in the SCONE Model configuration.

Integrators

Overview

Integrators advance the simulation by updating the velocity and position/orientation of each body, based on the linear and angular acceleration that is the result of the sum of all forces. Formally, an integrator $I_h$$I_h$ updates a set of generalized coordinates $q_t$$q_t$ with derivatives ${\dot{q}}_t$$\dot{q}_{t}$ at time $t$$t$ based on accelerations ${\ddot{q}}_t$$\ddot{q}_{t}$ and step size $h$$h$:

The integrator $I_h$$I_h$ is referred to as a fixed-step integrator: during each step, the state is advanced with with a constant time interval. While this approach is common among physics simulation engines, it leaves the difficulty of choosing the step-size $h$$h$ to the user. Step-size that are too small are inefficient, while step-sizes that are too large cause the simulation to become unstable. In practice, the optimal step-size often varies greatly during the course of a simulation, depending on the numerical stiffness at any point in time.

Variable-step integrators circumvent that issue by adapting the step-size to the current state of simulation. A variable-step integrator $I_A$$I_A$ computes the step-size $h$$h$ based on a set of accuracy requirements $A$$A$:

Variable-step integrators are an important part of Hyfydy, and allow for a stable and efficient trade-off between performance and accuracy.

Integration Methods

Forward Euler

The most straightforward integration method is the Forward Euler method, which updates positions $q_t$$q_t$ using velocities ${\dot{q}}_t$$\dot{q}_{t}$, which in turn are updated using accelerations ${\ddot{q}}_t$$\ddot{q}_{t}$:

The downside of this approach is that the positions are updated using a poor estimate of the velocity during interval $h$$h$. This can lead to an accumulation of errors and energy being added to system.

Symplectic Euler

The Symplectic Euler or semi-implicit Euler method improves over the Forward Euler method by updating velocities first and using the updated velocities ${\dot{q}}_{t+h}$$\dot{q}_{t+h}$ to update positions $q_t$$q_t$:

This approach leads to increased stability and energy conservation over the Forward Euler method. However, the position update is still based on a poor estimate of the velocity during the step (even with constant acceleration), and as such this method does not improve accuracy.

Midpoint Euler

The Midpoint Euler method uses the average or midpoint velocity between to update the position, leading to a better position estimate:

With this approach, both position and velocity updates are most accurate if acceleration ${\ddot{q}}_t$$\ddot{q}_{t}$ remains constant over interval $h$$h$. Despite being more accurate, it is typically outperformed by the Symplectic Euler method in terms of stability and energy conservation. However, in combination with a error-controlled variable-step integration strategy, the increased accuracy can be considered more important.

Higher-order Methods

In contrast to the first-order integrators described above, high-order integrators attempt to increase accuracy by evaluating the forces and resulting accelerations at different points in time $t$$t$.

This section is under construction

Muscle Activation Dynamics

In addition to integrating body position and velocity, integrators in Hyfydy also handle muscle activation dynamics. The muscle activation and deactivate are properties that are part of the integrator:

Given muscle activation $a_t$$a_t$ and excitation $u_t$$u_t$ at time $t$$t$, then activation is updated to $a_{t+h}$$a_{t+h}$ according to:

in which activation_rate corresponds to $c_1+c_2$$c_1 + c_2$, while deactivation_rate corresponds to $c_2$$c_2$.

Planar Integrators

Planar integrators only update positions in the x-y plane, and orientations in around the perpendicular z-axis. When using planar models with planar forces, these types of integrators increase simulation performance without loss of accuracy. See Integrator Configuraton for more details.

Integrator Configuration

In Hyfydy, the integrator can be configured via a configuration file. In SCONE, this configuration is directly added to the Model.

This section is under construction

Comparisons

Hyfydy vs OpenSim

Even though Hyfydy is based on the same models as OpenSim¹, there are some differences in implementation that result in slightly different simulation outcomes. As a consequence, results from a Hyfydy simulation are not directly transferable to OpenSim, and vice-versa. While this does not seem ideal, it is important to stress that both simulators are equally valid, and it usually only takes a small number of optimization iterations to convert from one result to another. In the remainder of this section, we provide an overview of all significant differences between Hyfydy and OpenSim.

Joint constraints

OpenSim uses generalized or reduced coordinates to describe the positions and velocities of the articulated bodies in the system, which automatically enforces joint constraints. In Hyfydy, each body has six degrees of freedom, and joint constraints are enforced explicitly through forces that mimic the effect of cartilage, bony structures and tendons. This results in small displacements within joints – just as in biological joints. The amount of displacement is based on the joint stiffness and damping settings, which can be set individually per joint or computed automatically.

Joint limits

Hyfydy uses non-linear damping forces to enforce joint limits, where the amount of damping is proportional to the limit force. In contrast, the OpenSim CoordinateLimitForce uses a constant damping coefficient beyond the limit threshold. The behavior in Hyfydy is similar to damping in Hunt-Crossley contact forces ⁵, and has similar benefits over the OpenSim CoordinateLimitForce.

Friction force

The friction force in the OpenSim HuntCrossleyFroce is based on an unpublished model attributed to Michael Hollars. Hyfydy provides an exact implementation of this model via contact_force_hunt_crossly_sb, which can be used to emulate the behavior of OpenSim friction.

Muscle force

Hyfydy uses an optimized implementation of the Millard Equilibrium Muscle Model ⁴, which differs from the OpenSim implementation in two significant ways:

1. The curves that describe the force-length and force-velocity relations, as well as the curves for passive tendon and muscle forces, are defined through polynomials instead of splines. The resulting curves in Hyfydy therefore differ slightly from the curves used in the OpenSim implementation.

2. The muscle damping forces are computed explicitly instead through the iterative method in the original OpenSim implementation. The resulting muscle damping forces in Hyfydy therefore differ slightly from the forces produced in the OpenSim implementation.

Muscle activation dynamics

Hyfydy uses a different approach for computing muscle activation dynamics than OpenSim. As a result, muscle activation patterns can differ between Hyfydy and OpenSim -- even when excitation patterns are identical. This difference only occurs when the deactivation rate is different from the activation rate.

Numeric Integration

OpenSim and Hyfydy both implement several variable-step integrators with user-configurable error control. In Hyfydy, the accuracy criterium is based on the highest error found across all bodies, while OpenSim uses a weighted sum of errors to determine accuracy. As a result, only in Hyfydy the error of each body is guaranteed to be below the specified accuracy threshold.

Version History

Version history will be documented after the 1.0.0 release.

References