Unsteady simulation problem

[Julio Jerez]

Oh you can do more than 5 rag doll, the demo has 50, and when every thing is active it take 12 MS.
but the is in single core, with tow or three core you can cut that by half.

Can you post the exact code of you power ragdoll, so that I am make a test in the sandbox demo?

[JoeJ]

Oh you can do more than 5 rag doll, the demo has 50, and when every thing is active it take 12 MS.
but the is in single core, with tow or three core you can cut that by half.

Can you post the exact code of you power ragdoll, so that I am make a test in the sandbox demo?

i'll copy past my code from other thread and see ift it reproduces...

[JoeJ]

I changed StandartJoins.cpp and added my ragdoll joint.
Only difference is about avoiding GrammSchmidt, see comment, that's not important.

It adds 3 powered "ragdolls"

One without animation and 5 bodies - that's fine.
Same with animation - lay it with "tail" on ground to see some jumpy motion.
And one with 20 bodies - touch it a little to see the bug.

[Julio Jerez]

Yes I see the problem, there is one thing that I tweak in the joint
when the angle is small you use an identity matrix and the base for rotation, I change to the
rotate parent frame, the does not ne the extra matrix , like this

Code: Select all

struct JoesRagdollJoint
{
   NewtonBody* m_body0;
   NewtonBody* m_body1; // parent
   dMatrix m_localMatrix0;
   dMatrix m_localMatrix1;
   NewtonJoint* m_joint;

   dQuaternion m_target; // relative target rotation to reach at next timestep

   float m_reduceError;
   float m_pin_length;
   float m_angularFriction;
   float m_stiffness;

   float m_anim_speed;
   float m_anim_offset;
   float m_anim_time;

   JoesRagdollJoint (NewtonBody* child, NewtonBody* parent, const dMatrix &localMatrix0, const dMatrix &localMatrix1, NewtonWorld *world)
   {
      m_body0 = child;
      m_body1 = parent;
      m_localMatrix0 = localMatrix0;
      m_localMatrix1 = localMatrix1;
      m_target = dQuaternion (dVector(1.0f,0,0), 0.0f);
      m_joint = NewtonConstraintCreateUserJoint (world, 6, JoesRagdollJoint::Callback, 0, child, parent);
      NewtonJointSetUserData (m_joint, (void*) this);

      m_reduceError = 0.2f; // amount of error to reduce per timestep (more -> oszillation)
      m_angularFriction = 300.0f;
      m_stiffness = 0.9f;

      m_anim_speed = 0.0f;
      m_anim_offset = 0.0f;
      m_anim_time = 0.0f;

//     m_basis = dGetIdentityMatrix();
   }

   float CalculateAngle (const dVector& dir, const dVector& cosDir, const dVector& sinDir) const
   {
      float cosAngle = dir % cosDir;
      float sinAngle = (dir * cosDir) % sinDir;
      return atan2(sinAngle, cosAngle);
   }

   void SubmitConstraints(dFloat timestep, int threadIndex)
   {
      float invTimestep = 1.0f / timestep;

      dMatrix body0Matrix; NewtonBodyGetMatrix(m_body0, &body0Matrix[0][0]);
      dMatrix body1Matrix; NewtonBodyGetMatrix(m_body1, &body1Matrix[0][0]);
      dMatrix matrix0 = m_localMatrix0 * body0Matrix;
      dMatrix matrix1 = m_localMatrix1 * body1Matrix;

      if (m_anim_speed != 0.0f) // some animation to illustrate purpose
      {
         m_anim_time += timestep * m_anim_speed;
         float a0 = sin(m_anim_time);
         float a1 = m_anim_offset * 3.14f;
         dVector axis(sin(a1), 0.0f, cos(a1));
         //dVector axis (1,0,0);
         m_target = dQuaternion(axis, a0 * 0.5f);
      }

      // point to point constraint

      dVector p0(matrix0.m_posit);
      dVector p1(matrix1.m_posit);
      NewtonUserJointAddLinearRow(m_joint, &p0[0], &p1[0], &matrix1.m_front[0]);
      NewtonUserJointAddLinearRow(m_joint, &p0[0], &p1[0], &matrix1.m_up[0]);
      NewtonUserJointAddLinearRow(m_joint, &p0[0], &p1[0], &matrix1.m_right[0]);

      // measure error
      dQuaternion q0(matrix0);
      dQuaternion q1(matrix1);
      dQuaternion qt0 = m_target * q1;
      dQuaternion qErr = (q0.DotProduct(qt0) < 0.0f
         ? dQuaternion(-q0.m_q0, q0.m_q1, q0.m_q2, q0.m_q3)
         : dQuaternion(q0.m_q0, -q0.m_q1, -q0.m_q2, -q0.m_q3)) * qt0;

      dQuaternion qErr1 = ((q0.DotProduct(qt0) < 0.0f) ? dQuaternion (-q0.m_q0,  q0.m_q1,  q0.m_q2,  q0.m_q3) * qt0 : dQuaternion ( q0.m_q0, -q0.m_q1, -q0.m_q2, -q0.m_q3)) * qt0;

      float errorAngle = 2.0f * acos(max(-1.0f, min(1.0f, qErr.m_q0)));
      dVector errorAngVel(0, 0, 0);

      dMatrix basis;
      if (errorAngle > 1.0e-10f) {
         dVector errorAxis(qErr.m_q1, qErr.m_q2, qErr.m_q3, 0.0f);
         errorAxis = errorAxis.Scale(1.0f / dSqrt(errorAxis % errorAxis));
         errorAngVel = errorAxis.Scale(errorAngle * invTimestep);

         basis = dGrammSchmidt(errorAxis);
      }
      else
      {
         basis = dMatrix (qt0, dVector (0.0f, 0.0f, 0.0f, 1.0f));
      }

      dVector angVel0, angVel1;
      NewtonBodyGetOmega(m_body0, (float*)&angVel0);
      NewtonBodyGetOmega(m_body1, (float*)&angVel1);

      dVector angAcc = (errorAngVel.Scale(m_reduceError) - (angVel0 - angVel1)).Scale(invTimestep);

      // motor
      for (int n = 0; n < 3; n ++)
      {
         // calculate the desired acceleration
         dVector &axis = basis[n];
         float relAccel = angAcc % axis;

         NewtonUserJointAddAngularRow(m_joint, 0.0f, &axis[0]);
         NewtonUserJointSetRowAcceleration(m_joint, relAccel);
         NewtonUserJointSetRowMinimumFriction(m_joint, -m_angularFriction);
         NewtonUserJointSetRowMaximumFriction(m_joint, m_angularFriction);
         NewtonUserJointSetRowStiffness(m_joint, m_stiffness);
      }
   }

   static void Callback (const NewtonJoint* joint, float timestep, int threadIndex)
   {
      JoesRagdollJoint *custom = (JoesRagdollJoint*) NewtonJointGetUserData(joint);
      custom->SubmitConstraints (timestep, threadIndex);
   }
};


let me see what is wrong with solver, the pops are ver bad.

[Julio Jerez]

here is a image to explain what is going on

spline.png (199.06 KiB) Viewed 5752 times

I listed wrong, the cut start form the bottom which the root note.
what happens is that the code see that second joint form the root violate the limit and removed that for the active set and res tart the solver,
the it as find the nest nest force break the limit and do the same, then as if goes up the stack, it keep finding more force breaking the limit.

If instead it did a reach for the worse offender, it worse find that the force in the middle break the limit the worse, and after removing that the force on the bottom does not need to be that strong anymore, but that requires two passed, one for estimating the worse force, and one the update them. I tried to update while removing and that is a big mistake


I thought, that the force calculate in of importance start from the root of the array, but the is no true, because the force is not only a function of the acieration require to enforce the constrain, when there are motors, the acceleration of the right size is arbitrary of the I can no assume any order, instead I need to do tow passes, on pass to fin the worse forces and another to update that node and possible bring in force that maybe dead.
this dopes no seems good new because it makes the algorithm square.

I nee to implement Lempel again to make comparison of the correctness , I have a version floating around form newton 1.5 but I do not know where is was I though is was on the general vectors class but it is not. In any case I will be much better that a general Danzig or Lemkel, but it will no be liner,
instead I will be proportion to the square of the unilateral constrains.
Then I can add heuristic like flagging motors

This is also the bug with the Tank demo too, I did not see before because I did no has motors but as soon as I added then the jerk motion shoued up.

[JoeJ]

Interesting.
One idea to fight the performance drop coming from fixing this would be a sparse net of hard joints over a dense net of soft joints.
For the snake from the picture that would maen a hard joint from body 1 to body 5, from 5 to 9 and so on.

For a ragdoll it would mean a hard joint from foot to hip, hip to torso, torso to hand.
That might not look 100% realistic in some extreme conditions, but the main problem would be the requirement af a 'super joint', which can power both orientation and position.
I tried that a while ago and it was not stable, so i came up with a joint that only powers rotation about one axis and keeps a minimum distance (can push but not pull).
This works and cuts the 'soft joints' problem to a half with very few additional rows.

Tricks like that might become handy here too...
But we'll see. I'm happy enough if you can get it to work

[Julio Jerez]

Hey Joe,
I have being thinking about a powered ragdoll controller for a few days.
and one of the problem I have is how to relate the inverted pendulum to the dynamic articulated body skeleton.

I have being reading paper but I have no found a single one that do the method that I thing about, almost every one use the proportional derivative controller which I believe is probably the worse method I can thing off. I have tried using pdf mumbo-jumbo in the pass and each and is always happen he you calibrate a set of going for a fix set of controller and all it take is add or remove one parameter for all the calibration to go to the toilet.

However this morning is just dawn on me that this is quite easy, My idea is to make the Inverted Pendulum itself a control joint.
This method complete avoid all of the calculation, that an analytical approach has to do, all I nee to do is to control limits and velocity, the physic solve will do all the calculation for me.
It is so easy that I am I can not even believe it, I am going to try to do it this weekend.

I have the feeling that s going to be very cool because as it usually is when and Idea is a real solution for a problem, the ideas itself is surpricilly trivial to the point that you say how did not I think of that before.

By having the inverted pendulum be part of the skeleton, the doing the IK form animation is very trivial, veasicall it all comes to extracting the kinematic information form the animation and convert the information to joint angles.
Blending animation is quite trivial because the system is base of solving fro the difference between the target angle common form an animation clip, of a procedural animation, and the current skeleton pose.
Now that we have the skeleton solver that is fast and strong I believe this is going to be very nice feature of the character controller.

[JoeJ]

one of the problem I have is how to relate the inverted pendulum to the dynamic articulated body skeleton

To me this is one of the easier problems.
E.g. i use inverted pendulum to get the desired ankle angle:
Divide ragdoll bodies in 2 groups - the bodies above and the bodies below the ankle joint.
For each group build one virtual single body with it's properly summed up physical properties.
Create Inverted Pendulum from the 2 virtual bodies and calculate ankle torque / velocity from control target.
Use that results to build a pose for the real ragdoll, and try to conserve physical quantities from the IP at control target.

Sounds easy, but it involves the IK solvers, some proper interpolation from two real ankles to one IP ankle, and the problem to calculate inertia for a group of bodies (i still loop over all bodies because a inertia matrix can't be build in my opinion as discussed in some other thread).
The entire mapping between ragdoll and IP is not based on analytical math and will introduce some error to me, but i don't worry it's relevant.

The only place whre i'm exactly analytical is the IP itself.
As the upper body moves, the distance to the edge of the foot changes, so does the max acceleration that can be used without tipping over.
So the IP is a system with position dependent acceleration, and the equations to control such a system become very very complex. I've used tools like Mathematica to get them.
(Allthough i think the win over the approximate equations i've used before may not be worth it - i'll check that back when it works)

I'm back to global illumination right now and take a break from physics.
The final plroblem i have with the ragdoll is that i don't know how to predict the small error or temporal lag that Newton gives me.
This problem is there, no matter if i use a 20 bodies ragdoll, or a simple 2 boxes IP model directly.
Seems no big problem, but i have no idea at the moment.

You need to bulid a two body IP model yourself and try to rotate the upper body to target angles as fast as possible without loosing balance.
As fast as possible means the center of pressure, visualized by summed foot contacts weighted by contact force, has to be almost at the edge of the foot -> best performance AND highest risk to loose balance at the same time.
You will be surprised how hard this is, and humans do this all the time.
I'm sure this is the key problem to be solved, no matter on any other details.


For the rest of your post i'm not sure if i understand but we may be similar.
E.g. i could use animation too and the IP controller / ik solver would care to maintain balance automatically (but my goal is to replace animkation completely - no time or mony for animation ).

About the hard joint skeleton:
I think the reason i did not get it to work for balancing ragdoll is simple:
The contact forces are not part of the skeleton, but they are as important as anything inside the skeleton.
Because the skeleton is so strong, there is additional stress on the contacts and so they fail to keep accurate enough.

A solution might be: The possibility to get a good inital guess for contact forces by callback from the user, which the solver can refine during it's iterations.
This way i could calculate the guess from a virtual single body of the whole ragdoll.
This would also make it possible to run at 60 Hz. I need at least 90 - at 60 the contact forces are to bad anyways.
If that sounds a bad idea to you, maybe it could be done with the user contact stuff you plan to add.
(let me know if something similar or better is already happening inside Newton)

What i do at the moment is this trick: Do not use hard joints, instead use distance joints (1 row) from the feet to ALL other bodies. 20 Joints to one body And it is stable without issues!
Also it is rock solid stiff against gravity, just like using hard joints.

[Julio Jerez]

To me this is one of the easier problems.
E.g. i use inverted pendulum to get the desired ankle angle:
Divide ragdoll bodies in 2 groups - the bodies above and the bodies below the ankle joint.
For each group build one virtual single body with it's properly summed up physical properties.
Create Inverted Pendulum from the 2 virtual bodies and calculate ankle torque / velocity from control target.
Use that results to build a pose for the real ragdoll, and try to conserve physical quantities from the IP at control target.

ha but the is the p[art that is confusing to me, and I believe is ultimately wrong.

The idea that you build tow bodes one with a fake center of mass made of all the bodies part, is frat out wrong and has many problems.
1-all the body part from and single mass. there are connected by joint and each joint accent forces on each other, pretending the form a single body is just wrong
2-the com is always changin position so the equation is hard for formulate
3- I read some paper and what the do about that is the come of with a mapping function for the IP, to the bodies and vise verse. but they is not a clear and each paper hind how the go about is totally different and obscure way. some systems go to the extreme hat after the do all that nonsense the system is no stable and they fake weight the gain of the PDF the mass and inertia of the skeleton. This is never explained and made no sence physically.

My method is different, all of the part are physical and do not changes.
the inverted pendulum is a joint that place the pelvis on a stick just like you do, by the pelvis is always on the stick.
this joint is more like a inverts rope joint the top body rest on the pendulum length and is free to fall on any direction horizontally and it keep the distance between the pelvis and the ground constant as long as the pelvis assert force on the ground, but If the force point upwork then it act like a contact.
This joint is part of eh skeleton at all time.

The joint of the controller now is to balance the pendulum using the other body parts.
I explain the first controoler that I have in mind. The Balance condoner

The balance controller will work like this.
with will assume that some body part are rigid and the one that are no rigid are ignored.

wit this information it will calculate the center of mass of all connect parts, (I may ignore he low body part, but this is not necessary)

this new center and the center of the inverted pendulum for a pair of point particles that allow the control to calculate the direction in witch the fake center has to move to prevent the pelvis center to keep falling.
this direction is not map to the different bodies parts using biomechanical heuristics.
if the angle if the direction is capable to stop or river the direction of motion of the pelvis, the body maintain the equilibrium. if it can not the this is a signal that the center of the IP must be changed discreetly, (this is like stepping to keep the balance)
This is no different that the balancing of a pendulum of a car excepts that the card is now the compound COM of the body parts which in general will be above the pelvis center.
all control are base on control parametric joint angle, and all signal come for the geometry of the skeleton.
My objective is not to make the skeleton to general motion by physical forces, like you are try to do. although I believe that if I subside on my first step, then you goal is just another controller.
My goal is to control the skeleton with poses give by a frame, the frame will be given by an extern source like an animation, procedural, and the environment.
so in my method the physic solve the legs just to keep them in position, the knees and leg joint are no rigid. the friction with the found is only use to determine the character speed but they can slide.
an example will be a running animation, but I am getting ahead of me now.

Let me get the first step in, maybe we can joint effort in getting a cool physical base animation system.

[Julio Jerez]

The cool thing of all of this is that if I am correct, the Character controller will not need the stupid huge capsule is will be the skeleton.
I can envision a contrail that simple push and the pelvis and the skeleton moves in the direction that the player want obeying all obstacles, and adapting to the evirimnet naturally.

the player could move robotically if it is procedurally, but is can be very organic if the delta position are take from a motion capture sequence, and the adaption to the complex evirimnet is 100% a emerging behavior for eh physical simulation.

The only concern is performance since the player is not loge a simple body, but that is the price the system will have to pay.
In he sandbox demo, in a single thread the engine can do 100 ragdoll made of 19 body part in 17 milliseconds per frame. and my computer is about 8 years old now. and it was a entry level system even when I bought it. so I believe this will not be a problem.

on the good side it allow for the unification of the entire scene physically which is in my opinion the biggest problem today.

[Julio Jerez]

For the rest of your post i'm not sure if i understand but we may be similar.
E.g. i could use animation too and the IP controller / ik solver would care to maintain balance automatically (but my goal is to replace animkation completely - no time or mony for animation ).

yes but I believe that this is a mistake, it is better to work with the animation instead of against animation.
The biggest problem with animation is transition blending that always look stupid.
The system I am proposing, I believe that reduce the number animation to a minimum.

for example say you wan the play to wall, in my system the walk animation uses the joint angle o the body part to get the target joint angles, that only the goal, the physics solver will minimize the error between the goal, the constraint imposed by the controller and the environment constraints.
what is different now is that the solve is capable to really find the solution and pose to all other system that the phsy engine is no capable to drive the skeleton with high fidelity.

Imaging a transition form wall to run, any animation system can be use to do that, as along as the transition points are math (must system do that) the final pose is just a pose not different that any other, the controller see that as a new target pose that it need to move too.
It can even deal with missing part, for example say you cut the arm, and the arm will simple dangle.

Basically this system will use the conservation of momentum to do the blending and transitions, so as long as the animation is decent them the final animation should look as natural a can be possible.

so the animation is just a way to get a sequence of natural looking poses.
and is use the strength of the physic solver to it advantage.

the reason I thong all real time animation system *, is that no physics engine has a solver the can solve convert the key frame form an animation to drive the skeleton, at leas no in real time. if that was available I believe that animations system will be a dime a dozen and will be every where.
if this works, maybe we should write a paper.

[JoeJ]

1-all the body part from and single mass. there are connected by joint and each joint accent forces on each other, pretending the form a single body is just wrong
2-the com is always changin position so the equation is hard for formulate

It's correct only under the assumption that all joints are stiff and prevent relative movement between the unique bodies. Because the mapping is rebuild every frame and relative movement is small it's at least good enough. Com IP == com ragdoll and momentum also equals between the according parts.
So the IP has all information necessary to feed the controller.
The controller result is then used to predict the IP positions at next frame, and mapping this back to the ragdoll may introduce some error that i'm willing to ignore for now.
(e. g. the IP can not change height - only the ragdoll can do this)

However, im curious to see what you'll come up with...

What you say sounds a bit like Havoks character stuff, did you try this?
Don't know if it still exists in the form i remember from the time i came to Newton,
but they had functionality to map animation to ragdolls (also allowing simplified ragdoll skeleton and more detailed animation skeletion).
I'm pretty sure there was nothing about balancing or walking, but i found it pretty impressive.

[Julio Jerez]

What you say sounds a bit like Havoks character stuff, did you try this?
Don't know if it still exists in the form i remember from the time i came to Newton,
I'm pretty sure there was nothing about balancing or walking, but i found it pretty impressive.

I heard of that and also the euphoria and natural motion, I have no seen havoc by I had seen natural motion, but the do no use the pendulum as far as I know.
In any case the Havok and natural motion is really awful. They show goo looking videos but all I have seen do no come near what that videos show.

Havok has always trumpeted there soft keyframe monbojumpbo, but in reality what it is is the use the softness of the physic solve to fake elastic joints.
I work on a game before "All start vs Legends" and also Madagascar, and I can tell you all animator hate the system. because when the map the animation to the skeleton the solve can not really play the animation they way the animation was created. basically the soft key frame filet out all of the high frequencies from the animation and the animation look featureless.

The animator always start vey excited, and it take a process of six moth for three reality to sink in and it is delegate to death, animation and in some cases reaction to punch of bullet.
This is a very frustration process, because those technologies are expensive, and the technical are very supported dishonest and simple blame the shortcoming of the technology on the programmer,

I crates very awkward situations, because usually the person who license that kind of technology and even the core team of programmer know very lither about applied math, and the believe what technical support says. What I want is a system capable to play a skeleton physically based 100% o the time. I know this is hard but that is my goal.

do you have any site with info about the havoc animation I seem some of the Euphoria and Natural motion but the do no say anything about just videos. I must say the videos are quite impressive.

[JoeJ]

do you have any site with info about the havoc animation

No. But Havok did not come close to Euphoria. It just tried to play back anmiation (i remember the low pass effect you mention). Looking at their site they don't mention character control anywhere now. I guess it's still there but not promoted.

I don't know much about Euphoria. I tried their animation tool ("Endorphin"?) years ago and found it's purpose is mainly to make automatic transitions between animation clips.
But their realtime stuff may be very different.
Someone wrote that their system needs to learn desired behaviours, and the game studio needs to work with Natural Motion people for some time to accomplish that.
Their videos are awesome, but because they keep everything secret and it's rarely used i would be very suspiciously.

[Julio Jerez]

I will have to add a new feature to newton 3.14.
basically newton until newton 3.14 the equation that the engine solve is

Code: Select all

1) M * a = Fe + Jt * f
2) J dot a  + a0 = 0

from 1 you can get:
a = inv (M) * Fe + inv(M) * Jt * f

and replacing in two
J * inv (M) * Fe + J * inv(M) * Jr * f + a0 = 0

and
J * inv(M) * Jt * f = - J * inv (M) * Fe - a0

which is the familiar rigid body equation.

I never bothered adding soft joint constraints since the solver was already soft and the last thing I wanted to have was lots of people complaining about tuning parameters like we saw on the early days with engines like ODE. To me adding softness to softness makes a fussy system and I did not want to deal wit that.
Now that we do have a linear time solver, we can use the more general version of the contraint equation: this

Code: Select all

J dot (a)  + C * f +  K * x + a0 = 0

here K and C are a diagonal matrices of spring stiffness and damper that allows the solver to violate a constraint forces. These technique is called matrix regularization and is used for avoiding singular values when solving very large matrix systems. These are two new vectors that are part of the joint. as K and C set to zero by default. then a new joint function can written to values to generate the joint stiffness and dampening property.
In fact the generalization of the constraint equation applied to tetrahedral elements is the basics of many Finite Element Solvers. Any way what the method does for now is that we can have
stiff springs damper joint that can be solved as part of the general system solver.
The general idea is that we make the skeletal joint with joint that are 100% rigid where they have to be and soft where they need not to.

example, say we what a hinge joint to control the arm and the lower are:
this joint will have 5 dof that are 100% rigid, an one dof which is powered soft joint with high stiffness. This controller them only need to be concerned with the powered joint parametric value.
This will be the target dof vector.

I believe these are the two elements needed for a full controllable physically based system.
the IP and soft joints that actually work on the desired DOF.