# 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^{1}, 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^{2} 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 ^{3} with hill-type musculotendon dynamics^{4} and non-linear contact forces^{5}. Control parameters were optimized in SCONE^{2} using 1000 generations of CMA-ES^{6} 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:

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the model | empty |

`gravity` | vector3 [m/s^{2}] | Acceleration 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:

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the body | empty |

`mass` | number [kg] | Body mass | required* |

`inertia` | vector3 | Diagonal components of the inertia matrix | required* |

`density` | number [kg/m^{3}] | Density used to calculate the body mass and inertia, together with shape | required* |

`shape` | Shape | Shape used to calculate the body mass and inertia, together with density | required* |

`pos` | vector3 [m] | Initial position of the body center-of-mass | zero |

`ori` | quaternion | Initial orientation of the body | identity |

`lin_vel` | vector3 [m/s] | Initial linear velocity of the body center-of-mass | zero |

`ang_vel` | vector3 [rad/s] | Initial angular velocity of the body | zero |

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

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the joint | empty |

`parent` | string | Name of the parent body | required |

`child` | string | Name of the child body | required |

`pos_in_parent` | vector3 [m] | Position of the joint in the parent body frame-of-reference | required |

`pos_in_child` | vector3 [m] | Position of the joint in the child body frame-of-reference | required |

`ref_ori` | quaternion | Reference orientation of the child body wrt the parent body | identity |

`stiffness` | number [N/m] | Stiffness property of the joint constraint force | model_options |

`damping` | number | Damping property of the joint constraint force | model_options |

`limits` | range [deg] | Rotational joint limits, expressed as a vector3 of ranges | none |

`limit_stiffness` | number [Nm/rad] | Stiffness property of the joint limit force | model_options |

`limit_damping` | number | Damping property of the joint limit force | model_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**

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

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the model | empty |

`type` | Shape | Shape of the geometry | required |

`body` | string | Name of the body to which the geometry is attached | required |

`material` | string | Name of the material associated with the contact geometry | default |

`pos` | vector3 [m] | Position of the geometry in the body frame-of-reference | `[0 0 0]` |

`ori` | quaternion | Orientation of the geometry relative to the body | identity |

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:

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the model | empty |

`static_friction` | number | Value used to calculate the static friction force | `1` |

`dynamic_friction` | number | Value used to calculate the dynamic friction force | `static_friction` |

`stiffness` | number | Stiffness of the restitution force | `1e6` |

`damping` | number | Damping of the restitution force | `1` |

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

Identifier | Type/Unit | Description | Default |
---|---|---|---|

`mirror` | boolean | Indicate whether components should be mirrored (0 or 1) | `0` |

`scale` | number | Value by which the model is scaled | `1` |

`density` | number [kg/m^{3}] | Default density used for bodies | `1000` |

`joint_stiffness` | number [N/m] | Default stiffness for joints | `0` |

`joint_damping` | number | Default damping for joints | `1` |

`damping_mode` | choice | The way in which default damping is computed | `critical_mass` |

`joint_limit_stiffness` | number [Nm/rad] | Default limit stiffness for joints | `0` |

`joint_limit_damping` | number | Default limit damping for joints | `0` |

`limit_damping_mode` | choice | The way in which default limit damping is computed | `critical_inertia` |

`muscle_force_multiplier` | number | Factor by which muscle `max_isometric_force` is multiplied | `1` |

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

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the model | empty |

`max_isometric_force` | number [N] | Maximum isometric force of the musculotendon unit | required |

`optimal_fiber_length` | number [m] | Optimal fiber length of the musculotendon unit | required |

`tendon_slack_length` | number [m] | Tendon slack length of the musculotendon unit | required |

`pennation_angle` | number [rad] | Pennation angle at optimal fiber length, in radians | `0` |

`stiffness_multiplier` | number | Multiplier applied to passive tendon an muscle elastic forces | `1` |

**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

The notation *rotation vector*, which corresponds to the 3D *rotation axis* scaled by *rotation angle* in radians.

A `joint_motor`

component can contain the following properties:

Identifier | Type | Description | Default |
---|---|---|---|

`joint` | string | name of the joint | required |

`stiffness` | number [N/m] | Stiffness for position-dependent torque ( | `0` |

`damping` | number [Ns/m] | Damping for velocity-dependent torque ( | `0` |

`max_torque` | number [Nm] | Maximum torque magnitude that can be applied by the motor ( | +infinity |

`target_ori` | quaternion | Target orientation for the motor ( | identity |

`target_vel` | vector3 [rad/s] | Target angular velocity ( | `[0 0 0]` |

`torque_offset` | vector3 [Nm] | Base torque ( | `[0 0 0]` |

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

Identifier | Type | Description | Default |
---|---|---|---|

`file` | string | Mesh filename | empty |

`shape` | Shape | Mesh shape | empty |

`pos` | vector3 [m] | Position of the mesh in the body reference frame | `[0 0 0]` |

`ori` | quaternion | Orientation of the mesh in the body reference frame | identity |

`color` | color | Color of the mesh, in format `[r g b a]` | unspecified |

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

Identifier | Type | Description | Default |
---|---|---|---|

`name` | string | Name of the degree-of-freedom | required |

`source` | string | Default 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 |

`default` | number [m or deg] | Default value of the dof | `0` |

`range` | range [m or deg] | Minimum 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…

`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 `joint_force_pnld`

component applies a force

The force `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 ^{5}, and scales the limit damping by the limit displacement torque:

in which the joint displacement angle

The constants `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

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 component | Planar joint force equivalent |
---|---|

`joint_force_pnld` | `planar_joint_force_pnld` |

`joint_force_pd` | `planar_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:

**Collision detection**. Geometries are checked for intersections, and contacts between geometries are recorded.**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`

:

Identifier | Type | Description | Default |
---|---|---|---|

`enable_collision_between_objects` | bool | When set to 1, intersections between each pair of geometries are detected. | 0 |

**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.

Material Property | Hyfydy | Simbody |
---|---|---|

`stiffness` ( | ||

`damping` ( | ||

`static_friction` , `dynamic_friction` ( |

#### 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

The tangent force `viscosity`

parameter

The total contact force

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^{5} contact force provides a better model of the dependence of the coefficient of restitution on velocity, based on observed evidence. Given material stiffness

The friction force

##### 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

##### 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.

Identifier | Type | Description | Default |
---|---|---|---|

`transition_velocity` | number [m/s] | Transition velocity | `0.1` |

**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.^{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:

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.

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:

Identifier | Type | Description | Default |
---|---|---|---|

`xi` | number | Muscle damping factor, in the paper referred to as | 0.1 |

`v_max` | number [optimal_fiber_length/s] | Maximum muscle contraction velocity | 10 |

`use_pennation_during_equilibration` | bool | Use pennation angle during muscle equilibration | 0 |

#### muscle_force_gh2010

This is an implementation of the Hill-type muscle model described by Geyer & Herr^{3}. 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^{7}. 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

The integrator *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

*Variable-step integrators* circumvent that issue by adapting the step-size to the current state of simulation. A variable-step integrator

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

The downside of this approach is that the positions are updated using a poor estimate of the velocity during interval

#### 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

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

#### 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

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

Identifier | Type | Description | Default |
---|---|---|---|

`activation_rate` | number [s^{-1}] | The muscle activation rate | 100 |

`deactivation_rate` | number [s^{-1}] | The muscle deactivation rate | 25 |

Given muscle activation

in which `activation_rate`

corresponds to `deactivation_rate`

corresponds to

### 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.

Integrator | Description | Properties |
---|---|---|

`forward_euler_integrator` | Fixed-step integrator using the Forward Euler method | |

`symplectic_euler_integrator` | Fixed-step integrator using the Symplectic Euler method | |

`midpoint_euler_integrator` | Fixed-step integrator using the Midpoint Euler method | |

`planar_symplectic_euler_integrator` | Planar version of `symplectic_euler_integrator` |

*This section is under construction*

## Comparisons

### Hyfydy vs OpenSim

Even though Hyfydy is based on the same models as OpenSim^{1}, 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:

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.

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.

Version | Date | Description |
---|---|---|

## References