Critically Damped PD controller

[raigan2]

Hi,
I'm having trouble trying to add a critically damped harmonic oscillator (spring/PD-controller) constraint to my simulator, and the references I've found are somewhat conflicting. Currently I'm using regular PD controllers outside of the solver (i.e the controllers are driving motor torques) and it's proving quite difficult to get well-behaving results, especially when the environment varies.

The original reference is here, which describes using CFM/ERP terms to simulate given P and D values: http://ode.org/ode-0.5-userguide.html#sec_3_8_2

My solver is based on the PBD/Kalman Filter as presented by Antonio on this board (R and S terms correspond to CFM and ERP respectively): http://www.bulletphysics.com/Bullet/php ... =cfm#p6001

Basically I'm unable to figure out what the correct formula is for calculating CFM/ERP values which will produce a critically damped spring; I've also got no clue as to where to begin deriving my own, as it seems to depend on the integrator as well as the constraint solver..

ODE gives (h is timestep, P and D are spring and damping coefficients):

Code: Select all

CFM =   1 / (D + h*P)
ERP = h*P / (D + h*P)

Poking inside Box2D's distance joint gives:

Code: Select all

CFM =       1 / (h*D + h*h*P) 
ERP = C * h*P / (h*D + h*h*P)

David Wu (PhysicsAndAnimationGDC2003.doc) gives:

Code: Select all

kp  =         P / (1 + h*D + h*h*P)
kd  = (D + h*P) / (1 + h*D + h*h*P)

Box2D and ODE should be similar enough to use the same terms, yet the former has everything scaled by an additional 1/h and includes the current constraint error (C) while the latter doesn't.. I'm not clear on why.

The above would allow the simulation of a harmonic oscillator; a critically damped oscillator requires additional calculation of appropriate P and D terms. An alternate parameterization of PD controllers (frequency+damping rather than spring coefficients ks+kd) is presented in Box2D and Wu's presentation:

Box2D:

Code: Select all

P = 4 * pi*pi * mass * frequency*frequency
D = 4 * pi * mass * damping * frequency

Wu:

Code: Select all

P = 9 * frequency*frequency
D = 4.5 * damping * frequency

I'd like to know how I might go about deriving such formulas for my own solver, and/or which of the above formulas are correct. At this point I'm quite confused

[Erin Catto]

Box2D's soft distance constraint is using impulses rather than forces, hence the 1/h term.

With a frequency/damping ratio model you just set the damping ratio to 1 to get critical damping.

To use frequency + damping ratio, you need to know the effective mass.

Here's how I derived the relationship between ERP/CFM and a mass-spring-damper system.

First, model a 1D particle with a 1D weld constraint. Solve for the final velocity symbolically in terms of the initial velocity and the ERP/CFM parameters.

Second, model a 1D MSD system and compute the velocity update using implicit Euler.

You now have two velocity update formulas. The first has ERP/CFM and the second has mass, spring, and damping coefficients. It turns out you can determine a correspondence between these formulas and parameters by inspection.

[raigan2]

Box2D's soft distance constraint is using impulses rather than forces, hence the 1/h term.

Of course, I'm sorry -- I forgot that ODE was forces/LegrangeMultiplier... I've been away from this stuff for about a year

With a frequency/damping ratio model you just set the damping ratio to 1 to get critical damping.

This is about the only thing I did learn from the various resources

To use frequency + damping ratio, you need to know the effective mass.

Ah.. but for most purposes the effective mass is extremely non-trivial and changes dynamically based on the configuration of the world. I was hoping there was a way to get the solver to figure out and use the stiffest frequency possible (with damping always set to 1). I may be confused but it seems like the solver should sort of "know" the effective mass at each joint, or at least that information is implicitly available, and an approximation could be built up over iterations based on the corrections that have been applied and the current remaining error.

Here's how I derived the relationship between ERP/CFM and a mass-spring-damper system. <..>

Thanks, I'll give this a shot. Frankly I have a very poor intuition regarding calculus, which means that implicit integration and analytical MSD systems are still somewhat over my head

[Erin Catto]

Here are some more details (check my work).

MSD: http://en.wikipedia.org/wiki/Damping

m * vDot + c * v + k * x = 0

w = sqrt(k/m), z = c / (2 * sqrt(k * m))

vDot = -2 * z * w * v - w * w * x
xDot = v

Implicit Euler:

v2 = v1 - h * (2 * z * w * v2 - w * w * x2)
x2 = x1 + h * v2

Substitute second into first:

v2 = v1 - h * (2 * z * w * v2 - w * w * (x1 + h * v2))

Now solve for v2 (left to you).


Constrained system:

Newton:
m * vDot = lambda (force)

Velocity constraint:
v + (beta/h) * x + gamma * lambda = 0

Symplectic Euler:
v2 = v1 + h * lambda / m
x2 = x1 + h * v2

Constraint:
v2 + (beta/h) * x1 + gamma * lambda = 0 (ERP + CFM)

Unknowns: x2, v2, lambda
Solve for v2 and compare with IE.