hyfydy-user-manual

Hyfydy User Manual

Introduction

About

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 50-100x faster than OpenSim1, 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 SCONE2 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 controller3 with hill-type musculotendon dynamics4 and non-linear contact forces5. Control parameters were optimized in SCONE2 using 1000 generations of CMA-ES6 on a 4-core i7 Intel CPU.

Citation

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

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

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 =):

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:

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):

range

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

Core Components

model

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

A model can contain the following properties:

IdentifierTypeDescriptionDefault
namestringName of the modelempty
gravityvector3Acceleration due to gravity[0 -9.81 0]

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:

IdentifierTypeDescriptionDefault
namestringName of the bodyempty
massnumberBody massrequired*
inertiavector3Diagonal components of the inertia matrixrequired*
densitynumberDensity used to calculate the body mass and inertia, together with shaperequired*
shapeShapeShape used to calculate the body mass and inertia, together with densityrequired*
posvector3Initial position of the body center-of-masszero
oriquaternionInitial orientation of the bodyidentity
lin_velvector3Initial linear velocity of the body center-of-masszero
ang_velvector3Initial angular velocity of the bodyzero

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

Example

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

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:

IdentifierTypeDescriptionDefault
namestringName of the jointempty
parentstringName of the parent bodyrequired
childstringName of the child bodyrequired
pos_in_parentvector3Position of the joint in the parent body frame-of-referencerequired
pos_in_childvector3Position of the joint in the child body frame-of-referencerequired
ref_oriquaternionReference orientation of the child body wrt the parent bodyidentity
stiffnessnumberStiffness property of the joint constraint forcemodel_options
dampingnumberDamping property of the joint constraint forcemodel_options
limitsrangeRotational joint limits, expressed as a vector3 of rangesnone
limit_stiffnessnumberStiffness property of the joint limit forcemodel_options
limit_dampingnumberDamping property of the joint limit forcemodel_options

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

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

geometry

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

IdentifierTypeDescriptionDefault
namestringName of the modelempty
typeShapeShape of the geometryrequired
bodystringName of the body to which the geometry is attachedrequired
materialstringName of the material associated with the contact geometrydefault
posvector3Position of the geometry in the body frame-of-reference[0 0 0]
oriquaternionOrientation of the geometry relative to the bodyidentity

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:

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:

IdentifierTypeDescriptionDefault
namestringName of the modelempty
static_frictionnumberValue used to calculate the static friction force1
dynamic_frictionnumberValue used to calculate the dynamic friction forcestatic_friction
stiffnessnumberStiffness of the restitution force1e6
dampingnumberDamping of the restitution force1

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:

IdentifierTypeDescriptionDefault
mirrorbooleanIndicate whether components should be mirrored (0 or 1)0
scalenumberValue by which the model is scaled1
densitynumberDefault density used for bodies1000
joint_stiffnessnumberDefault stiffness for joints0
joint_dampingnumberDefault damping for joints1
damping_modechoiceThe way in which default damping is computedcritical_mass
joint_limit_stiffnessnumberDefault limit stiffness for joints0
joint_limit_dampingnumberDefault limit damping for joints0
joint_stiffnessnumberDefault stiffness for joints0
limit_damping_modechoiceThe way in which default limit damping is computedcritical_inertia
muscle_force_multipliernumberFactor by which muscle max_isometric_force is multiplied1

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:

IdentifierTypeDescriptionDefault
namestringName of the modelempty
max_isometric_forcenumberMaximum isometric force of the musculotendon unitrequired
optimal_fiber_lengthnumberOptimal fiber length of the musculotendon unitrequired
tendon_slack_lengthnumberTendon slack length of the musculotendon unitrequired
pennation_anglenumberPennation angle at optimal fiber length, in radians0
stiffness_multipliernumberMultiplier applied to passive tendon an muscle elastic forces1

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

Example

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 τ is based on base torque τo, stiffness kp, damping kd, orientation q, target orientation qt, angular velocity ω, and target velocity ωt:

τ=[τo+kpr(q1qt)+kd(vtv)]τmax

The notation [ ]τmax is used to indicate that the magnitude of the final torque is clamped between [τmax,τmax]. The function r:R4R3 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:

IdentifierTypeDescriptionDefault
jointstringname of the jointrequired
stiffnessnumberStiffness for position-dependent torque (kp)0
dampingnumberDamping for velocity-dependent torque (kd)0
max_torquenumberMaximum torque magnitude that can be applied by the motor (τmax)+infinity
target_oriquaternionTarget orientation for the motor (qt)identity
target_velvector3Target angular velocity (vt)[0 0 0]
torque_offsetvector3Base torque (τo)[0 0 0]

Example

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:

IdentifierTypeDescriptionDefault
filestringMesh filenameempty
shapeShapeMesh shapeempty
posvector3Position of the mesh in the body reference frame[0 0 0]
oriquaternionOrientation of the mesh in the body reference frameidentity
colorcolorColor of the mesh, in format [r g b a]unspecified

Example

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:

IdentifierTypeDescriptionDefault
namestringName of the degree-of-freedomrequired
sourcestringDefault name of the dof. This consists of the name of the body appended by _rx, _ry or _rz for rotation relative to its parent, or _tx, _ty and _tz for translation of a root body. Dof values can be inverted by prepending the source with a minus sign (-).required
defaultnumberDefault value of the dof0
rangerangeMinimum and maximum values of the dof. Note that these values are not enforced during the simulation, use a joint limit force for that instead.-90..90

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

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

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, the joint_force_pnld component applies a force Fj based on joint displacement pj and displacement velocity vj:

Fj=kppjkvvj

The force Fj 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 τj based on joint angle αj and angular velocity ωj. While the limit torque is proportional to joint displacement angle θ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 force5, and scales the limit damping by the limit displacement torque:

τj=kαθj(αj,αmin,αmax)(1+kωωj)

in which the joint displacement angle θj is defined as:

θj(αj,αmin,αmax)={αjαmin,if αj<αminαjαmax,if αj>αmax0,otherwise

The constants kp, kd, kα and kω 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:

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 ωj:

τj=kαθj(αj,αmin,αmax)kωω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:

Joint force componentPlanar joint force equivalent
joint_force_pnldplanar_joint_force_pnld
joint_force_pdplanar_joint_force_pd

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:

IdentifierTypeDescriptionDefault
enable_collision_between_objectsboolWhen set to 1, intersections between each pair of geometries are detected.0

Example

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.

Material PropertyHyfydySimbody
stiffness (k1, k2)k1k2k1+k2k1k2k1+k2
damping (c1, c2)c1+c22k1k1+k2c1+k2k1+k2c2
static_friction, dynamic_friction (μ1, μ2)μ1μ22μ1μ2μ1+μ2

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 and dissipation coefficient c, the normal force Fn is derived from penetration depth d and normal velocity vn, which is the contact velocity in the direction of contact normal n:

Fn=kdcvn

The tangent force Ft depends on the static friction coefficient μ, viscosity parameter η and tangent velocity vt:

Ft=vtvtmin(ηvt,μFn)

The total contact force Fc corresponds to:

Fc=Fnn+Ft

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

contact_force_hunt_crossley

The non-linear Hunt-Crossley5 contact force provides a better model of the dependence of the coefficient of restitution on velocity, based on observed evidence. Given material stiffness k, penetration depth d and normal velocity vn the contact restitution force Fn corresponds to:

Fn=kd32(132cvn)

The friction force Ft and resulting contact force Fc are defined as in contact_force_pd.

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.

IdentifierTypeDescriptionDefault
transition_velocitynumberTransition velocity0.1

Example

Actuator Forces

muscle_force_m2012fast

This is an optimized implementation of the Hill-type muscle model described by Millard et al.4, 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 4 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:

IdentifierTypeDescriptionDefault
xinumberMuscle damping factor, in the paper referred to as β0.1
v_maxnumberMaximum muscle contraction velocity [m/s]10
use_pennation_during_equilibrationbooleanUse pennation angle during muscle equilibration0

muscle_force_gh2010

This is an implementation of the Hill-type muscle model described by Geyer & Herr3. 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 Thelen7. 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.

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 Ih updates a set of generalized coordinates qt with derivatives q˙t at time t based on accelerations q¨t and step size h:

{qt+h,q˙t+h}=Ih(qt,q˙t,q¨t,h)

The integrator Ih 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 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 IA computes the step-size h based on a set of accuracy requirements A:

{h,qt+h,q˙t+h}=IA(qt,q˙t,q¨t,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 qt using velocities q˙t, which in turn are updated using accelerations q¨t:

qt+h=qt+hq˙tq˙t+h=q˙t+hq¨t

The downside of this approach is that the positions are updated using a poor estimate of the velocity during interval 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 q˙t+h to update positions qt:

q˙t+h=q˙t+hq¨tqt+h=qt+hq˙t+h

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:

q˙t+h=q˙t+hq¨tqt+h=pt+h2(q˙t+q˙t+h)

With this approach, both position and velocity updates are most accurate if acceleration q¨t remains constant over interval 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.

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:

IdentifierTypeDescriptionDefault
activation_ratenumberThe muscle activation rate [s1]100
deactivation_ratenumberThe muscle deactivation rate [s1]25

Given muscle activation at and excitation ut at time t, then activation is updated to at+h according to:

at+h=at+h(utat)(c1ut+c2)

in which activation_rate corresponds to c1+c2, while deactivation_rate corresponds to c2.

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.

IntegratorDescriptionProperties
forward_euler_integratorFixed-step integrator using the Forward Euler method 
symplectic_euler_integratorFixed-step integrator using the Symplectic Euler method 
midpoint_euler_integratorFixed-step integrator using the Midpoint Euler method 
planar_symplectic_euler_integratorPlanar version of symplectic_euler_integrator 

This section is under construction

Comparisons

Hyfydy vs OpenSim

Even though Hyfydy is based on the same models as OpenSim1, 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 5, 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 4, 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.

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.

VersionDateDescription
   

References

 


1 Seth, A., Hicks, J. L., Uchida, T. K., Habib, A., Dembia, C. L., Dunne, J. J., … Delp, S. L. (2018). OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. PLoS Computational Biology, 14(7). https://doi.org/10.1371/journal.pcbi.1006223
2 Geijtenbeek, T. (2019). SCONE: Open Source Software for Predictive Simulation of Biological Motion. Journal of Open Source Software, 4(38), 1421. https://doi.org/10.21105/joss.01421
3 Geyer, H., & Herr, H. (2010). A muscle-reflex model that encodes principles of legged mechanics produces human walking dynamics and muscle activities. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 18(3), 263–273. **https://doi.org/10.1109/TNSRE.2010.2047592
4 Millard, M., Uchida, T., Seth, A., & Delp, S. L. (2013). Flexing computational muscle: modeling and simulation of musculotendon dynamics. Journal of Biomechanical Engineering, 135(2), 021005. https://doi.org/10.1115/1.4023390
5 Hunt, K. H., & Crossley, F. R. E. (1975). Coefficient of Restitution Interpreted as Damping in Vibroimpact. Journal of Applied Mechanics, 42(2), 440.
6 Hansen, N. (2006). The CMA evolution strategy: a comparing review. Towards a New Evolutionary Computation, 75–102.
7 Thelen, D. G. (2003). Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. Journal of Biomechanical Engineering, 125(1), 70–77. https://doi.org/10.1115/1.1531112