Differential Equation Solving and Constraints

Danske
2006-07-20
2013-04-24
  • Danske
    Danske
    2006-07-20

    I mentioned these papers; here's the link:
    http://www.cs.cmu.edu/~baraff/sigcourse/

    It's a little annoying that they assume you have a equation solving function, but they describe enough in general to allow one to be written, I think.

    The bit about varying the time step relative to the change in values is nice in that it gives an equation for how much to vary the time step by.  In retrospect it's obvious that the error depends on the time derivative of the variables.  If velocities are high then you need a small step size for accurate positions.  Similarly for accelerations and velocity, and d(force)/dt and acceleration.

    It also discusses higher-order Runge-Kutta for solving differential equations.  I was skeptical of the practicability of this in a physical simulation, but if you can express the whole physical system in a (relatively) nice system of equations then that might be fine.  The more accurate integration also means the step size could be higher, but for collisions you'd still want the time step small enough so that the error in objects' positions is small relative the the size of the objects (you don't want objects to miss each other that should have hit in the middle of the time step).

    There's also a good discussion of solving for constraint functions, although they say computationally you just plug them into a generic quadratic program solver. :p

    Another quaternion derivative is also in an appendix.

    They also discuss the Lagranian method, but say it's good for off-line modelling of objects and not within the running simulation.

    It stuff like that which depresses me.  Everything's already been done. ;)

     
    • Martin Baker
      Martin Baker
      2006-07-20

      Thanks, this looks like a really useful reference.

      We may have the theory to simulate most mechanical systems (eventually). I wont consider the job done until I can write a program which allows me to place any virtual component (gears, pulleys, springs, dampers, joints, motors, etc.) into the program and then start it running without any customised programming.

      At the moment its so hard to define the behaviour of these objects, because we have to keep translating between global and different local frames, then we have to keep converting between vectors, tensors, spinors and (in the case of the Lagranian method) individual scalars.

      The methods all have their advantages but no one method suits all situations, some are inaccurate, some are computationally expensive.

      I'm looking for a sort of retro (Newtonian) theory of everything.

      Martin

       
    • Danske
      Danske
      2006-07-21

      A generic black-box object would have to have some inputs and outputs.  It would have some hook function for "compute force and torque" to introduce forces and torques into the system, which would need to query the state of the system.  It would also need to be able to introduce constraints into the system.

      If we consider an jet engine, neglecting its exhaust on other physical objects, it can perform work "on itself", introducing a force in the system that acts only on itself.  Since this is increasing the momentum and kinetic energy of the system this is outside the scope of the general body-body simulation.

      That would be a specific case of the more general case of an object introducting a force on any arbitrary object in the system.  If we wanted the exhaust to affect other objects, then the jet engine (the object responsible for the exhaust behavior) would have to query the system state and be able to introduce forces on other objects as well.  You'd still have to code the way an object introduces forces to the system, but they'd be in a generic "jet engine", "spring", etc. class, and the simulator can just query them for the introduced forces.

      Collision detection and resting contact are also just specific instances of solving for contraints on the system.  More complex behavior of objects rather than being just rigid bodies would introduce new constraints.  If the simulator can introduce dynamically constraints for collisions and contact, then it seems it could also query objects for any constraints they add to the system, or which they change over time.

      Two meshing gears, for example, constrain the velocity of the contact point tangential to their radiuses to zero, i.e. the linear velocity of each gear at the contact point is the same.

      A gear box would constrain the ratio between the angular velocities and angular acceleration of its input and output shafts relative to itself.  It tries to maintain this constraint by reacting with a torque on the two shafts.  This also introduces a torque on itself, which is why you need to bolt it down if you don't want it rotating too. ;)

      Objects can introduce forces and constraints in just the same way as collision detection does.  They query the state of the system and either directly add forces (like a spring), or add constraints (joints, etc.) which need forces/torques added to satisfy them.  If the simulator can dynamically add constraints and forces, then it can ask objects to do the same; you just need a generic interface for doing so.

       
    • Danske
      Danske
      2006-07-24

      Just some thoughts as I've found some other papers to read.  A lot of a simulation is tied up in constraints on the state of the system and solving for forces to maintain those constraints, with those forces subject to their own constraints.

      Non-penetration and joints both put constraints on the position/orientation of objects, as well as their velocities and accelerations (both linear and angular).  The laws of physics in the simulation put constraints on the forces being solved for.

      Non-penetration means that objects may be in contact, but not overlapping.  At the contact point the objects must also not be moving or accelerating toward each other; this is just solving the non-penetration constraint for higher time derivatives of position/orientation.  Joint constraints are similar, just with different positions/orientations, velocities and accelerations that are allowed.

      The main constraint on constraint forces is that we don't want constraint forces to increase the kinetic energy of the system.  They can decrease the KE, simulating loss of KE in a collision or through friction, or maintain KE.  I think in all but collisions this is satisfied by never applying a constraint force in the direction of velocity or acceleration.  If (P=F*v) is negative or zero then no energy is being added to the system.

      However, in collisions we're not solving for force applied over time, but a large transient force applied over an infintesimal amount of time which results in a nearly instantaneous acceleration of the object (impulse force).

      I'm a little opposed to just changing the velocities of the objects instantaneously, but since the contact forces would be so large and over such a small period of time they would "break" the simulation unless we used a very, very small timestep to simulate the collisions with very stiff springs instead.  So I guess I accept the need for applying impulses instead of trying to integrate forces over time.

      These impulses are also subject to constraints.  In addition to using them to satisfy the state constraints on the system, the impulses must not increase kinetic energy.  The restitution coefficient affects how much kinetic energy is lost.  The impulses may also cause other constraints to be violated, such as one from a joint or interpenetration, which will require solving for impulses at those constraint points as well.

      Frictional forces are also subject to constraints.  Static friction at a contact point with a relative velocity of zero between the two objects will resist tangential acceleration with a tangential force determined by the normal force (from the non-interpenetration constraint) and the coefficient of static friction.  Kinetic friction will supply a force in the opposite direction to the relative velocity of the two contacting objects, again determined by the normal force and the coefficient of kinetic friction.  If static friction is not sufficient to prevent acceleration, or kinetic friction causes the relative velocity to drop to zero, then we've transitioned from one to the other.  This transition is another constraint on the frictional forces.

      I'm not entirely sure how impulses interact with frictional constraints.  My first inclination is that static friction will supply a counteracting impulse up to some limit, just as it does for force, if some object being subject to an impulse is in stationary contact with another object and the impulse would introduce a relative velocity between the two objects.  In a collision between two objects with friction taken into account at the point of contact then the collision impulse will also have a tangential component if there is a tangential component to the contact velocity.

      It seems to me that the general method would be to solve for violation of velocity constraints with impulses, as impulse instantaneously fixes velocities, subjecting the impulses to constraints around kinetic energy (as KE is derived from velocity), and then solving for violation of acceleration constraints with forces.

       
    • Martin Baker
      Martin Baker
      2006-07-25

      Thanks, I think I generally agree with the approaches you are suggesting.

      Do you think it would be possible to work out a simulator program for the web site? Not necessary with graphics output just the physics animation. Or possibly just an pseudo program or outline.

      I guess the first thing to work out would be the data structure. In an earlier message you talked about object querying the system state. Do you have any thoughts about how the system state would be encoded?

      I've been coming to the conclusion that such as system could be coded as a scene graph (tree structure) cross linked with a table holding the state information.

      The scene graph would hold for each object:
      * Physical properties, mass, inertia matrix, coefficient of friction, elasticity, etc. in local coordinates.
      * Shape (boundary) in local coordinates.
      * Constrains, for instance a joint, would define limits to movement relative to other objects.
      * Any info required for rendering (textures, bump maps, normals).
      * cross link to state table.

      The state table would hold things that change every frame:
      * linear and angular position (in absolute coordinates).
      * linear and angular velocity (in absolute coordinates).
      * boundary sphere (in absolute coordinates).
      * External forces
      * cross link to scene graph.

      Do you have any thoughts on the data structure requires for such a simulation?

      Martin

       
      • Danske
        Danske
        2006-07-25

        Actually, I'd say the data structure comes out last from how you need to use it.

        On one hand we are solving for the position and velocity (linear and angular) for each object which can be done analytically, but on the other hand the objects are interacting, which requires solving a system of equations.

        For example, with two objects connected by a damped spring, you have a "phantom object" of the spring, connected to the two object.  I call it a phantom object because we may not wish to model it as a deformable object which is being collided with.  For the collision detector it's not really there, hence "phantom".  The damped spring doesn't constraint the motion of the objects, as such, but provides a known force (it doesn't need to be solved for as do constraint forces) based on the relative distance between two fixed points on the objects (spring) and how this distance is changing (damper).

        This requires the damped spring to query the system to determine the location and velocities of the two points on those objects it's connected to.  From that it can add in the two forces on the two objects for our calculation of the state of the system at time + delta_t.

        If we replace the spring with a rigid rod (phantom or not), then we also need to determine the acceleration of the points on the two objects so that they're not accelerating away from each other to satisfy the constraint of the rod's length.  As the force we're applying with the rod will affect the points' accelerations we need to calculate how a unit force along the rod will change those accelerations.  This creates a system of equations we need to solve.

        For just two objects we might be able to hard-code solving a system of just two equations, but the number of objects in the "contact chain" may be greater than just two.  I think this requires inserting all the parameters into a matrix to solve the system of n equations.  If we do model the rod as a real object then there will be at least three equations.

        Since a force applied at one constraint point may acclerate a body such that the constraints at other points are violated, and across several different contact points between bodies in the contact chain, we have to iteratively converge on a solution for all the forces.  At least I think we have to do it iteratively.  With a force between objects A and B affecting the constraints between objects B and C, and any force applied between B and C then affecting the constraints between A and B, we have to numerically converge on a solution rather than solve it analytically.

        We have the set of all objects in the system, some of which may be "phantom" objects and not analyzed by the collision detector.  Each object has its inherent physical properties and its state in global space (other than the phantom objects).  The forces in the system form a set of force pairs between objects.  Some of these forces may directly calculable from the state of the system (springs, gravity, jet exhaust, gyroscopic forces), while others need to be solved for subject to constraints.  The constrained forces would be in the class of contact forces needed to prevent interpenetration or joint forces needed to preserve the constraints of the joints; I think both can be treated together in a general way.  These force pairs are either introduced by the objects themselves (one object dynamically interacts with the other objects in a way not covered by collisions or joints), the collision detector, or through constant joints.

        Oh, and throw in the collision impulse calculations to prevent inter-penetration of objects.  Ah, I guess the same method by which we calculate how a force/torque produces acceleration in an object can be used for how an impulse will affect its velocity.  I think the equations would be identical.

        In general:

        1. Start with a system with no inter-penetrating objects.

        2. Examine collision points and their relative velocities, solving for impulses necessary to prevent velocities which are "toward" interpenetration.  May be messy if this includes friction.  Apply those impulses.

        3. Calculate and apply all forces which do not depend on any other forces (unconstrained), such as springs, dampers, gravity, motors, gyroscopic forces, magnetism, etc.

        4. Examine all the constrained interaction between objects in terms of contact points and force pairs.  Solve numerically for a solution that satisfies all constraints.

        5. Advance the simulation time.  May use the simple Euler for force/acceleration, assuming they're not changing over time, or the modified Euler for a simple average of how forces are changing over time, the Euler Predictor-Corrector to solve for an average force over the time step more accurately, or a method such as fourth-order Runge-Kutta for more initial accuracy.  Halt the simulation time before inter-prenetration has occured, possibly needing to "back up".  I'm not sure how this step should interact with the previous steps.  Anything other than simple Euler evaluates the forces at (t + delta_t) to estimate (dF/dt) and higher derivatives.  Calculating the unconstrained forces is not so time consuming, but re-solving for constrained forces may be so.  On the other hand, the more costly calculations allow us to use a larger step size in the simulation, so the simulation rate may be higher.  On the gripping hand, the previous steps may introduce discontinuities into the forces over time (impulse-forces certainly do) which aren't good for the ODE solving.  However, impulse forces were applied before we got here, and any further collisions necessarily halt this step to re-apply impulses.

        6.  Finish with a new system state with no inter-penetration at (t + delta_t).

        I'm working on exactly how to formulate constraints and constraint solving in equations/code.  If a generic quadratic program-solver were used the geometry of everything would all have to go in one big matrix.  Since we know conceptually the Big Matrix is formed out of several discrete objects and forces, perhaps it would be fine to work on those directly.  We may still need a generic simultaneous equation-solver as the number of force points change dynamically as the simulation runs.

        The calculations involve querying the position, velocity, and acceleration of a point on an object to determine if a constraint is satisfied, and how the object's acceleration responds to force applied at a particular point on the object to solve for forces to satisfy a constraint.  Once I've really worked out how the calculations are performed then I'll know the information needed, which will inform me of how that information should be stored. :)

        The problem is that there's really no simple intermediate case.  Collisions between frictionless rigid bodies is fairly trivial.  If you add static contact forces (which are a subset of general constraint forces) the problem balloons-up to a quadratic program where you have to simultaneously solve all the forces to preserve non-interpenetration.  Throw in fiction and it gets even worse.

        One example of some of the difficulties: imagine a cube sitting on a table.  To prevent interpenetration with the table we must determine forces which provide both zero velocity relative to the table at the four vertexes of the cube in contact with the table.  There is however no unique solution.  As long as the upward forces sum to the force from gravity to give zero linear acceleration to the cube, and the forces from diagonally-opposite points are equal to provide no torque to the cube, the cube will not penetrate the table.  Each point in contact could have a quarter of the weight of the cube on it, or two could have half-each and the other two none.  Depending on how we're attempting to solve for the forces we could either divide by zero somewhere along the way, or circle about never actually converging on a solution.

        Anyway, I'm working on doing the simple pendulum with more sophisticated constraints.  Because of the angular motion I'm afraid that constraints on linear velocities and accelerations at the pivot point will still allow the positional constraint to be violated.  I'm wondering how to fix that.  If it were considered a collision then an impulse would need to be applied just before the positional constraint were violated.  I'm not sure if I should wait for the "slop" value in the position to be violated to apply an impulse, or considered applying an impulse whenever the relative velocity is non-zero.  Even then, with high accelerations the position may integrate to being outside the allowed "slop" in one timestep.

        Oh, I did also play about with varying the time step based on the forces/velocities at play.  In the gyroscope simulation the angular acceleration only prompted a change in time step when the angular velocity was low.  This makes sense, but also points out that if you want to maintain a given accuracy for objects' positions you may still have to scale the time step relative to acceleration even though most of the time it's secondary to velocity.

         
    • Martin Baker
      Martin Baker
      2006-07-27

      It seems that there are a number of options in building a physics simulator, None of these options can cope with all types of physics, I think the ideal system would be able to mix these approaches and dynamically switch between them depending on the accuracy required and computing resources available.

      However, to be practical, I guess we have to start somewhere for example if we are building a 'car racer' which of the following options would you choose:

      Space
      -----
      Options are:
      a) finite steps, e.g.- model each atom, or one point represents say a million or a billion atoms.
      b) solid bodies - each solid body has 6 degrees of freedom, only boundary has to be modelled for collision detection
      c) make solid bodies as a set of simpler shapes, say cubes, or represent shapes and mass distribution by equations.

      Time
      ----
      Options are:
      d) finite steps, represent time varying quantities at step n, in terms of step n-1, n-2, etc.
      e) continuous, represent time varying by equations.

      Forces
      ------
      Options are:
      f) represent all interactions between objects as forces and groups of forces as force + torque.
      g) treat static contact, free movement, constraints (joints and sliding) and collisions separately each with there own algorithms.

      Do you think this represents a comprehensive set of options?
      I've got lots of views on the pros and cons of each of these options. However I cant yet come to a view about a combination which would be simple enough to implement but would produce a useful simulation?

      Martin

       
      • Danske
        Danske
        2006-07-28

        From the reading I've done it is often broken up into discrete concerns:

        Object interaction:

        1. Collisions: constraining objects from interpenetrating by way applying an impulse to keep objects in contact with each other from having velocities toward each other.

        2. Contact: constraining objects in contact with each other but with no velocity toward each other from interpenetrating by applying forces to keep relative acceleration at zero or away from each other.  Joints are also a form of contact constraints.

        3. Dynamic friction: applying a frictional force when two objects are in contact with each other with no normal (toward or away) velocity with a tangential relative velocity.

        4. Static friction: applying a frictional force when two objects are in contact with each other with zero tangential and zero normal velocity, to counteract any tangential acceleration between the objects (maintaining zero relative tangential velocity).

        Object properties:
        1. Rigid, non-deformable.  Constant geometry.

        2. Non-rigid, deformable.  Variable geometry the behavior of which must be modelled.

        For my semi-gamish car simulator I'd pick rigid bodies because I'm not interested in deformation.  Any deformation, e.g. of springs, would be modelled explicitly.

        I wouldn't use collision detection because the parts of the car shouldn't freely collide with each other, their movements being constrained by constant joints instead; no contacts between the rigid-body components should be made or broken.  I'm also not interested in hitting trees or barriers, just the internal behavior of the system of car objects.   Also, the geometry of the component objects isn't important aside from a graphical representation.  Since there's no collision detection I don't need to know the boundaries of the objects.

        The joint constraints would also be specific instances of contact constraints, so no general contact would need to be solved for; again, no contact should be made or lost from the initial joints.

        No friction; joint friction might be interesting to play with, but I'll assume static and dynamic friction are not significant forces.  The joints are well-oiled and don't bind.  Yet again, any contact between component objects is only through joints, and the joints have no friction, so no frictional forces need to be computed.

        The interaction with the road is determined by specific tire/road equations and not a generic friction algorithm; from the simulator's point of view the car never actually touches anything, in collision or contact, so no friction modelling is needed.

        Since I'm not interesting in the simulation rate, I'd go with a fourth-order Runge-Kutta with a small time step.  I think the joint forces would be sufficiently continously absent collisions that this wouldn't present any problems requiring special handling.

        If I wanted to make it more gamish I'd probably add collision detection but not friction.

        Once I'm happy playing with the mechanical behavior of the suspension, I might move on to playing with aerodynamics.  That might involve developing matrices giving lift/drag/torque for a given airflow vector, and possibly some rudimentary computation of ground-effect, but wouldn't model the shape of the aerodynamic surfaces or play with CFD analysis even off-line.

        As another example, if I wanted to do a curling game, I could do everything in limited dimensions.  I'd have to introduce collision detection and friction, but only in two dimensions.  The actual curling of the stones would require specialized abstraction because they don't curl the right way according to the usual Laws of Physics. ;)  Rather than dynamically integrating the force and friction around the base of the stone, a function given the linear and angular velocity would calculate the force/torque the stone is experiencing from friction.

        Yeah, it seems like you can easily take one cross-section of physics you're interested in simulating directly and then either abstract out or ignore the rest, and doing so may actually be required. (Consider a spherical cow of uniform density...)

        Another thing to consider is if the model is accurate.  You can simulate an object as a bunch of connected particles, but if you're trying to have the same behavior as rigid bodies then you're going to need a very small time step, or you'll just simulate spontaneously exploding bodies. ;)

        Even using "F=ma" for the movement of a rigid body is a simplification from directly modelling the individual particles of the body, but it's also very accurate unless you want to simulate the internal forces of the object.

        So you could model contact forces as objects just exerting a very strong, very short-distance repulsive force on each other, but then you need to slow the simulation down to deal with very small distances and very large forces.  If the time step is too high then the forces will be higher than they "should" have been, adding kinetic energy to the system as objects spontaneously fly apart.  If the forces are weak then you can use a larger time step, but objects will interpenetrate more before the repulsive force equals the force pushing the objects together.  If the distance the force operates over is too large then objects supposedly in contact will float apart instead.

        Similarly for joint forces.  They could be modelled as very stiff springs, which is okay if you're running with a very small timestep.  If the natural frequency of the spring and the objects it's attached to is near 1/timestep then there will be problems as it oscillates out of control.

        The "custom" algorithms are just using known equations of behavior that you don't want to simulate directly, integrating everything over very small distances and times.  Pick the chapter of the textbook you want to simulate and model no more and no less than you need to. ;)  Even CFD doesn't model indivdual atoms, just their known mass behavior to a reasonable degree of accuracy.

         
    • Martin Baker
      Martin Baker
      2006-08-03

      I am slowly trying to capture this information and put it on the website.

      I am working on these pages:
      http://www.euclideanspace.com/threed/games/
      for people who are interested in writing games and just want the minimum maths and physics to do that.

      and I'm still working on the physics theory on these pages:
      http://www.euclideanspace.com/physics/dynamics/momentum/

      The idea is to cross link the pages so that people get to the theory they need for their program. Or they could choose to start with the physics theory and then link to the games later.

      I know I've still got a lot of work to do to incorporate these ideas, it could easily be a full time job, I appreciate your help with this.

      Thanks,

      Martin