Hello,

at our institute we are using OpenSim for ergonomic design studies. Because we are aiming for a completely virtual analysis tool chain, it is not possible to use experimental data in our simulations. Instead motion sequences are computed using constrainted inverse kinematics. However I encountered a problem concerning inverse dynamics when using kinematic constraints. I constructed a simple case study in order to explain the problem:

The picture shows a simple planar system: a rigid body is connected to the ground using a revolute joint and the free end of the body is fixed to the ground via a point constraint. Obviously this closed loop system cannot move since the DOF of the joint is locked by the kinematic constraint. A forward dynamic simulation has been run to prove this behavior. Now when I do inverse dynamics using the "motion" file generated during FD, I get a strange result: the torque of the (actuated joint) does not equal zero! Instead I get exactly the torque necessary to balance the body against gravity. Apparently the constraint forces are ignored during inverse dynamics! I expected a torque of zero since the system is in rest and statically dertermined.

I hope someone can help...

Best regards,

Daniel

## kinematic constraint forces

- Michael Sherman
**Posts:**807**Joined:**Fri Apr 01, 2005 6:05 pm

### Re: kinematic constraint forces

Hi, Daniel. You are correct that OpenSim ignores constraints in ID. Probably it should obey them; we're discussing internally to determine why it was done this way and whether there is a reasonable workaround. Regards, Sherm

- Michael Sherman
**Posts:**807**Joined:**Fri Apr 01, 2005 6:05 pm

### Re: kinematic constraint forces

Hi, Daniel. After some discussion, it is clear that the correct behavior for ID+constraints is not entirely obvious. Although in your posted example the correct result seems clear, in general there are choices to make and we are interested in your thoughts. Here are a few questions:

Thanks and regards,

Sherm

- With the current behavior, any kinematic input is acceptable since the generalized coordinates are independent. What should we do if the given motion does not satisfy the constraints? In general it will be impossible for the motion to satisfy the constraints exactly (especially at the acceleration level where constraint forces must be determined), so we will probably have to accept some slop. But what would that mean? Should there be a "clean up" phase where the motion is made to follow the constraints exactly?
- From your post it is clear that you felt the constraint should carry as much of the load as possible. However the split between joint torques and constraint forces is indeterminate. Is it always preferable to assign the maximum possible forces to the constraints, leaving only the minimum necessary leftovers for the joint torques?

Thanks and regards,

Sherm

- Daniel Krueger
**Posts:**24**Joined:**Fri Jul 09, 2010 12:05 am

### Re: kinematic constraint forces

Hello Sherm,

thank you for your profound reply. Thinking about it again, we came to the same conclusion. I understand that it is hard to find a universally satisfying solution to this problem.

Honestly I have no idea what to anwers to your first question. Perhaps one could track the amount of constraint violation at each time step and abort the operation if this value exceeds a given threshold ?

Concerning the distribution of load between constraints and joints my suggestion is to minimize the load carried by the actuated DOF.

Best regards,

Daniel

thank you for your profound reply. Thinking about it again, we came to the same conclusion. I understand that it is hard to find a universally satisfying solution to this problem.

Honestly I have no idea what to anwers to your first question. Perhaps one could track the amount of constraint violation at each time step and abort the operation if this value exceeds a given threshold ?

Concerning the distribution of load between constraints and joints my suggestion is to minimize the load carried by the actuated DOF.

Best regards,

Daniel

- Michael Sherman
**Posts:**807**Joined:**Fri Apr 01, 2005 6:05 pm

### Re: kinematic constraint forces

OK, thanks. Any other OpenSim users have thoughts on this topic?

- Ton van den Bogert
**Posts:**166**Joined:**Thu Apr 27, 2006 11:37 am

### Re: kinematic constraint forces

I do not use Opensim for inverse dynamics but the problem seems somewhat fundamental and not limited to Opensim. We are working on inverse dynamic analysis of the shoulder mechanism. There is a closed kinematic chain (thorax-clavicle-scapula-thorax) and therefore more unknown joint moments than equations of motion. The way we are dealing with this is to do a combined inverse dynamic analysis and static optimization in one step, rather than sequentially.

The equations of motion are written for the open kinematic chain, and the reaction loads at the loop closure are treated as unknown loads, along with the joint moments. So there is a vector T that contains all these unknown loads. The equations of motion are:

M(q)*qdd + c(q,qd) + B(q)*T = 0

(I use Autolev, which gives me symbolic expressions for M, c, and B)

(d and dd stand for "dot'" and "dotdot", the first and second derivatives with respect to time)

B is a matrix with more columns than rows, so we can't solve T. The classical example would be a planar 4-bar linkage with one grounded link and all joints actuated. The open chain has 3 DOF, so 3 equations of motion, but T has 6 elements: 4 joint torques and X and Y force at the loop closure.

In a musculoskeletal model, T contains joint moments which are related to muscle forces F: T = R(q)*F. R is a matrix with all moment arms. The rows in R that correspond to the reaction loads in T will simply be all zeros. Now we can formulate the static optimization problem:

minimize a cost function f(F)

subject to linear constraints A*F = b, where

A = B(q)*R(q)

b = -M(q)*qdd - c(q,qd)

We're just starting on this, so I can't report results yet but I believe that this is a good and easy way to do it. I would be interested in Sherm's comments.

Alternatively, you could obtain the closed-loop equation of motion, so the reaction loads would not appear at all. You still would have more unknowns than equations of motion. For instance, the closed 4-bar linkage has 1 DOF but 4 unknown driving torques. For the shoulder, we did not want to do it that way because we need to know some of the reaction loads (and because it is much easier to get the open loop equations of motion). Some static optimization approaches include the constraint that the force between scapula and thorax (the loop closure force) must be compressive and normal to the sliding surface (van der Helm, J Biomech 1994). So we needed to bring this force into the equations anyway.

One aspect that is not quite elegant is that we estimate body pose q with an inverse kinematic analysis that does not enforce the kinematic constraint at the loop joint. But if the data is good anyway, that should be OK. In ergonomics, you are probably doing forward kinematics so you will need a way to generate consistent kinematics for the closed loop system. Opensim may be able to do the inverse or forward kinematics with full consideration of the closed loop mechanism.

If your model is torque-driven rather than muscle-driven, it can be done the same way but you then need to define a torque-based cost function f(T). The reaction loads are linearly related to the joint torques that are responsible for the movement (via the equation of motion), but they should be excluded from f(T) because there is no physiological cost. In Daniel's simple 1-segment example, the solution of this static optimization problem would be zero joint torque which is exactly what I would expect a human to do in this situation.

Apologies for the long post, but this is a good opportunity for me to get some early feedback on what we are doing. It may also help solve Daniel's problem or be useful to the Opensim developer team.

Ton van den Bogert

The equations of motion are written for the open kinematic chain, and the reaction loads at the loop closure are treated as unknown loads, along with the joint moments. So there is a vector T that contains all these unknown loads. The equations of motion are:

M(q)*qdd + c(q,qd) + B(q)*T = 0

(I use Autolev, which gives me symbolic expressions for M, c, and B)

(d and dd stand for "dot'" and "dotdot", the first and second derivatives with respect to time)

B is a matrix with more columns than rows, so we can't solve T. The classical example would be a planar 4-bar linkage with one grounded link and all joints actuated. The open chain has 3 DOF, so 3 equations of motion, but T has 6 elements: 4 joint torques and X and Y force at the loop closure.

In a musculoskeletal model, T contains joint moments which are related to muscle forces F: T = R(q)*F. R is a matrix with all moment arms. The rows in R that correspond to the reaction loads in T will simply be all zeros. Now we can formulate the static optimization problem:

minimize a cost function f(F)

subject to linear constraints A*F = b, where

A = B(q)*R(q)

b = -M(q)*qdd - c(q,qd)

We're just starting on this, so I can't report results yet but I believe that this is a good and easy way to do it. I would be interested in Sherm's comments.

Alternatively, you could obtain the closed-loop equation of motion, so the reaction loads would not appear at all. You still would have more unknowns than equations of motion. For instance, the closed 4-bar linkage has 1 DOF but 4 unknown driving torques. For the shoulder, we did not want to do it that way because we need to know some of the reaction loads (and because it is much easier to get the open loop equations of motion). Some static optimization approaches include the constraint that the force between scapula and thorax (the loop closure force) must be compressive and normal to the sliding surface (van der Helm, J Biomech 1994). So we needed to bring this force into the equations anyway.

One aspect that is not quite elegant is that we estimate body pose q with an inverse kinematic analysis that does not enforce the kinematic constraint at the loop joint. But if the data is good anyway, that should be OK. In ergonomics, you are probably doing forward kinematics so you will need a way to generate consistent kinematics for the closed loop system. Opensim may be able to do the inverse or forward kinematics with full consideration of the closed loop mechanism.

If your model is torque-driven rather than muscle-driven, it can be done the same way but you then need to define a torque-based cost function f(T). The reaction loads are linearly related to the joint torques that are responsible for the movement (via the equation of motion), but they should be excluded from f(T) because there is no physiological cost. In Daniel's simple 1-segment example, the solution of this static optimization problem would be zero joint torque which is exactly what I would expect a human to do in this situation.

Apologies for the long post, but this is a good opportunity for me to get some early feedback on what we are doing. It may also help solve Daniel's problem or be useful to the Opensim developer team.

Ton van den Bogert

- Michael Sherman
**Posts:**807**Joined:**Fri Apr 01, 2005 6:05 pm

### Re: kinematic constraint forces

Hi, Ton. I wish we had a "Like" button I could click for your post!

I really like the idea of skipping ID and going straight to the muscle force calculation. It might still be worth separating out the constraint reactions though, if it is true that we would like them to do as much of the work as possible (which is what Daniel was expecting in his example). We have also just started thinking about this, but this is what I had come up with so far, completely untested:

The constrained equations of motion are:

(1) M udot + ~G lambda + T = f

(2) G udot = b

where G is the mXn (m<n) constraint Jacobian, ~G means transpose(G), lambda is the Lagrange multipliers representing the unknown constraint reactions, T is the unknown joint torques and f collects gravity, applied forces, centrifugal forces and other knowns. Here udot (=qdd in your equations) is the given accelerations (cleaned up if necessary so that the constraints in eq. 2 are satisfied).

My thought was to attempt first to satisfy equation (1) using only the constraints. In that case we assume T=0 and then solve the overdetermined linear system

(3) ~G lambda = f - M udot

for lambda that minimizes the residual error. In Daniel's problem this would solve exactly since the constraint alone can prevent motion, but in general there would still be a residual requiring forces that are in the null space of the constraints. If those are to be determined as joint torques, we can just write them down now as

(4) T = f - M udot - ~G lambda

using the lambda we calculated from (3). Or you could solve the muscle activation problem at this point as you proposed.

I think this would have the effect that constraints would carry as much of the load as possible, with joint torques or muscle forces only filling in the remaining gaps. Whether this is a desirable outcome is still a question!

What do you think?

Regards,

Sherm

I really like the idea of skipping ID and going straight to the muscle force calculation. It might still be worth separating out the constraint reactions though, if it is true that we would like them to do as much of the work as possible (which is what Daniel was expecting in his example). We have also just started thinking about this, but this is what I had come up with so far, completely untested:

The constrained equations of motion are:

(1) M udot + ~G lambda + T = f

(2) G udot = b

where G is the mXn (m<n) constraint Jacobian, ~G means transpose(G), lambda is the Lagrange multipliers representing the unknown constraint reactions, T is the unknown joint torques and f collects gravity, applied forces, centrifugal forces and other knowns. Here udot (=qdd in your equations) is the given accelerations (cleaned up if necessary so that the constraints in eq. 2 are satisfied).

My thought was to attempt first to satisfy equation (1) using only the constraints. In that case we assume T=0 and then solve the overdetermined linear system

(3) ~G lambda = f - M udot

for lambda that minimizes the residual error. In Daniel's problem this would solve exactly since the constraint alone can prevent motion, but in general there would still be a residual requiring forces that are in the null space of the constraints. If those are to be determined as joint torques, we can just write them down now as

(4) T = f - M udot - ~G lambda

using the lambda we calculated from (3). Or you could solve the muscle activation problem at this point as you proposed.

I think this would have the effect that constraints would carry as much of the load as possible, with joint torques or muscle forces only filling in the remaining gaps. Whether this is a desirable outcome is still a question!

What do you think?

Regards,

Sherm

- Ton van den Bogert
**Posts:**166**Joined:**Thu Apr 27, 2006 11:37 am

### Re: kinematic constraint forces

Sherm,

I think we have the same idea, just different ways to formalize it.

I did not include your eq. (2), which states that the kinematic constraints must be satisfied for the accelerations. I assumed that this would already have been taken care of in the preceding inverse kinematics or forward kinematics. Either exactly or approximately. Sounds like you would do that also.

My approach is not limited to a muscle-driven model (although that is how I use it), but it can also work with a torque-driven model if you define a torque-based cost function.

In eq. 1, I combined your lambda and T into one vector T. It seems more clear to write the terms separately as you did. On the other hand, my formulation makes it more clear that the lambdas and T's are solved simultaneously from this underdetermined linear system. But the lambdas are not in the cost function.

As far as I can tell, your approach of minimizing the residual in eq. 1 (after eliminating the T term) is exactly equivalent to mine, in the special case that my cost function is the sum of squared joint torques. You wrote eq. (1) without a coefficient matrix in front of the T, which makes the error metric in the linear least-squares approach exactly equivalent to that cost function.

To help justify the approach (so it does not appear to be just a numerical trick), it may be better to present it as linear constraints due to dynamics, plus the minimization of a cost function with some physiological justification. And then you are no longer limited to a cost function that is the sum of squared torques. Someone might want to use cubed torques, weight them unequally, etc.

In the more general formulation, you can also write the equation of motion with a coefficient matrix in front of the T. This may be an advantage for an underactuated system such as gait (fewer T's than DOFs) but I have not really thought that through. For arm movement we have full actuation so we can write (1) as you did.

When the model has muscles, I would definitely suggest to use a muscle-based cost function. That way, any non-muscular contributions to the joint torques are "free", as they should be. (you could also have included them in your term f, to avoid that problem). Also, it seems inelegant to first do a static optimization to solve joint torques, and then another one to solve muscle forces from the torques.

After solving muscle forces directly from the equation of motion, you can still compute the joint torques as T = R(q)*F.

If the closed kinematic chain is inside the body (as in the shoulder), we cannot measure the joint torques or reaction loads to determine how good our answer is. There is only EMG. In Daniel's ergonomics problems, some of these unknown loads (the lambda's) are outside the body, so you can measure them, (and therefore also the joint torques), and validate these theoretical predictions.

Ton

I think we have the same idea, just different ways to formalize it.

I did not include your eq. (2), which states that the kinematic constraints must be satisfied for the accelerations. I assumed that this would already have been taken care of in the preceding inverse kinematics or forward kinematics. Either exactly or approximately. Sounds like you would do that also.

My approach is not limited to a muscle-driven model (although that is how I use it), but it can also work with a torque-driven model if you define a torque-based cost function.

In eq. 1, I combined your lambda and T into one vector T. It seems more clear to write the terms separately as you did. On the other hand, my formulation makes it more clear that the lambdas and T's are solved simultaneously from this underdetermined linear system. But the lambdas are not in the cost function.

As far as I can tell, your approach of minimizing the residual in eq. 1 (after eliminating the T term) is exactly equivalent to mine, in the special case that my cost function is the sum of squared joint torques. You wrote eq. (1) without a coefficient matrix in front of the T, which makes the error metric in the linear least-squares approach exactly equivalent to that cost function.

To help justify the approach (so it does not appear to be just a numerical trick), it may be better to present it as linear constraints due to dynamics, plus the minimization of a cost function with some physiological justification. And then you are no longer limited to a cost function that is the sum of squared torques. Someone might want to use cubed torques, weight them unequally, etc.

In the more general formulation, you can also write the equation of motion with a coefficient matrix in front of the T. This may be an advantage for an underactuated system such as gait (fewer T's than DOFs) but I have not really thought that through. For arm movement we have full actuation so we can write (1) as you did.

When the model has muscles, I would definitely suggest to use a muscle-based cost function. That way, any non-muscular contributions to the joint torques are "free", as they should be. (you could also have included them in your term f, to avoid that problem). Also, it seems inelegant to first do a static optimization to solve joint torques, and then another one to solve muscle forces from the torques.

After solving muscle forces directly from the equation of motion, you can still compute the joint torques as T = R(q)*F.

If the closed kinematic chain is inside the body (as in the shoulder), we cannot measure the joint torques or reaction loads to determine how good our answer is. There is only EMG. In Daniel's ergonomics problems, some of these unknown loads (the lambda's) are outside the body, so you can measure them, (and therefore also the joint torques), and validate these theoretical predictions.

Ton