System framework needs to accommodate non-analytical potential energy

sherm1 commented 4 years ago

Background

One of the few definitive validations we can perform of complicated models is to verify that energy is conserved. This is particularly important within MultibodyPlant where a wide variety of joints, force elements, and constraints must be perfectly implemented.

The conserved quantity of interest is E = U(t) + KE(t) − Wɴᴄ(t) where

U(q) is the potential energy (PE) from conservative force elements (e.g. gravity, springs, stiffness-dependent contact forces),
KE(q,v) is the kinetic energy of moving masses, and
Wɴᴄ(t) = ∫ Pɴᴄ(x,u) dt is the work done (signed) on the system by non-conservative model elements (e.g. dissipation from dampers, injected power from actuators). Pɴᴄ is instantaneous power from those elements.
(Noting that q=q(t), v=v(t), x=x(t), u=u(t)).

PE can also be defined in terms of "conservative power" Pᴄ such that U(t) = U(0) − Wᴄ(t) where Wᴄ(t) = ∫ Pᴄ(q,v) dt is the work done on the system (signed) by conservative elements. Wᴄ is not as useful for validation as is a separately-computed U(q) because Pᴄ is typically computed using forces produced by the elements and it is normally those forces we would like to validate.

The System API currently supports the above with these methods, implemented by each LeafSystem and summed up by each Diagram:

U   = CalcPotentialEnergy(context);
Pᴄ  = CalcConservativePower(context)
Pɴᴄ = CalcNonConservativePower(context)
KE  = CalcKineticEnergy(context);

The problem

The current System API assumes that we can always calculate PE analytically given the system kinematic variables q. That's true for simple force elements, but in general all or part of the PE may have to be computed by integrating a power (force x velocity) calculation. For example, hydroelastic contact has a potential energy (as proved in section X.A of this paper). However, that is an integral over the displaced volume and we have no closed form for that. It is easily calculated by integrating power, though. Another example is the symmetric bushing under development in #11946, where most of the PE has a closed form but where under some circumstances there is an additional term for which a closed form has not yet been discovered.

So for implementation purposes we must break our PE calculation into analytical and work integral portions: U(t) = Uᴀ(q) + Uɪ(t) where Uᴀ(q) is a known function but Uɪ(t) = −∫ Pᴄɪ(q,v) dt must be integrated numerically. We have to split up the conservative power computation similarly so that Pᴄ = Pᴄᴀ(q,v) + Pᴄɪ(q,v).

The problem is that we don't have a way to express this split computation in the current System API.

Proposal

The most straightforward change would be to modify the overrideable APIs to this set:

Uᴀ  = CalcPotentialEnergyAnalytical(context);
Pᴄᴀ = CalcConservativePowerAnalytical(context);
Pᴄɪ = CalcConservativePowerNumerical(context);
Pɴᴄ = CalcNonConservativePower(context)
KE  = CalcKineticEnergy(context);

then define these non-virtual methods:

Uɪ  (= −∫Pᴄɪ)        = GetPotentialEnergyNumerical(t);
Wɴᴄ (= ∫Pɴᴄ)         = GetNonConservativeWork(t);
U   (= Uᴀ + Uɪ)      = CalcPotentialEnergy(context);
Pᴄ  (= Pᴄᴀ + Pᴄɪ)    = CalcConservativePower(context);
E   (= U + KE − Wɴᴄ) = CalcEnergy(context);

While we're changing this API, let's consider switching to user-supplied lambdas with specified prerequisites as we do for output ports and cache entries.

Workaround

Currently for the symmetric bushing in #11946 we are misusing the current API like this:

Uᴀ        = CalcPotentialEnergy(context);
Pᴄᴀ       = CalcConservativePower(context)
Pɴᴄ + Pᴄɪ = CalcNonConservativePower(context)

cc @mitiguy @amcastro-tri

mitiguy commented 4 years ago

This looks really good -- very clear. Nit-picking: Perhaps we can find a better than w in Pᴄᴡ: Pᴄᴡ = CalcConservativePowerNumerical(context); I think w is too close to "work", "angular velocity" and "world". It is unclear how "w" associates with "conservative" or "numerical". Perhaps Pᴄʙ = CalcConservativePowerNumerical(context); Pᴄɪ = CalcConservativePowerIntegrate(context);

sherm1 commented 4 years ago

Thanks -- switched to Pᴄɪ and Uɪ per your suggestion, PTAL. I'm thinking the method names should still be Analytical vs. Numerical, which is the important distinction.

amcastro-tri commented 4 years ago

Probably asking for potential energies (whether analytical or not) is too much. Why not to just get away with "conservative power" and "non-conservative power"? then our integrators can just integrate those for energy budgets. If a particular system "knows" an analytical form, it would be reported at that system's level API for comparisons, but probably too much asking all systems to be able to provide that.

sherm1 commented 4 years ago

Well, no one has to provide the analytical PE since they can just provide the non-analytical conservative power. But there is little verification value in that since it is derived directly from the forces. Also the analytical PE can be correct to all decimal places while the integrated power will generally have large errors. When working on complicated elements like some constraints it is easy to forget cross terms which can be very small -- there is nothing like seeing them conserve energy to 15 decimal places to provide a warm feeling about correctness.

So I think we should continue to try to get analytical PE whenever possible. This will mostly affect elements of the MultibodyPlant where we should attempt to do the best energy bookkeeping we can.

RobotLocomotion / drake