There is one thing I can't quite wrap my head around with sequential impulse solvers:

Let's say a body has collided with the ground, the solver has adjusted the linear and angular velocities without restitution and it's been integrated such that the body is no longer colliding. At the start of the next frame, the body still has the linear and angular velocities required to resolve the collision in the last frame. This results in a tiny bounce. The bounce does not come from restitution, it's a result of the previous velocities required to resolve the collision, carrying over to the next frame.

I haven't found any online resources discuss this. Is there a good way to prevent this bounce or is it unavoidable with sequential impulse solvers?

The way to solve this is to have two different velocities. One is the usual velocity, and the other is an instantaneous “bias” velocity which is used to resolve penetrations and other position constraint errors (e.g. in joints). When integrating velocity to position, you add both velocities together before multiplying by the time step. After integration, the bias velocity gets reset to 0 so that it doesn't carry over to the next frame.

When you are solving the constraints, you apply different types of impulses to one of the two velocity vectors. Constraints on velocity (e.g. relative velocity along normal must be ≥ 0) modify the usual velocity, and constraints on position (e.g. distance along normal must be ≥ 0) modify the bias velocity.

Ah that makes sense! Thank you, I'll give that a try.

@aressera

Thinking through the implementation details of this brought up a question. When I'm calculating the magnitude of the impulse applied to the bias velocity (the bias lambda), I need the contact velocity. See the line of code below marked HERE.

for constraint_set in &fixed_constraint_sets {
	rigid_body := constraint_set.rigid_body;

	for constraint in &constraint_set.constraints {
		{ // Normal velocity correction
			contact_velocity := rigid_body.velocity + linalg.cross(rigid_body.angular_velocity, constraint.r);
			velocity_error_n := linalg.dot(contact_velocity, constraint_set.n);
			lambda_n := -velocity_error_n / constraint.effective_mass_inv_n;
					
			prev_total_impulse_n := constraint.total_impulse_n;
			constraint.total_impulse_n = max(constraint.total_impulse_n + lambda_n, 0.0);
			total_impulse_delta_n := constraint.total_impulse_n - prev_total_impulse_n;

			rigid_body.velocity += constraint_set.n * total_impulse_delta_n / rigid_body.mass;
			rigid_body.angular_velocity += rigid_body.tentative_inv_global_inertia_tensor * constraint.rxn * total_impulse_delta_n;
		}

		{ // Normal position correction
			bias_contact_velocity := rigid_body.velocity + rigid_body.bias_velocity + linalg.cross(rigid_body.angular_velocity + rigid_body.bias_angular_velocity, constraint.r); // HERE
			bias_velocity := linalg.dot(bias_contact_velocity, constraint_set.n);
			bias_lambda_n := -(bias_velocity + constraint.bias) / constraint.effective_mass_inv_n;

			prev_total_bias_impulse_n := constraint.total_bias_impulse_n;
			constraint.total_bias_impulse_n = max(constraint.total_bias_impulse_n + bias_lambda_n, 0.0);
			total_bias_impulse_delta_n := constraint.total_bias_impulse_n - prev_total_bias_impulse_n;

			rigid_body.bias_velocity += constraint_set.n * total_bias_impulse_delta_n / rigid_body.mass;
			rigid_body.bias_angular_velocity += rigid_body.tentative_inv_global_inertia_tensor * constraint.rxn * total_bias_impulse_delta_n;
		}
	}
}

To get the bias lambda n variable, I need the bias contact velocity. My question is about the correctness of this line:

bias_contact_velocity := rigid_body.velocity + rigid_body.bias_velocity + linalg.cross(rigid_body.angular_velocity + rigid_body.bias_angular_velocity, constraint.r);

Was I right to include both the regular and bias velocities here? This feels like the right thing to do. My thinking is, if the regular velocity isn't quite zeroed out, the bias impulse should account for that. It may need to be stronger than if the regular velocity is zeroed out. I could calculate the bias impulse in complete isolation of the regular impulse but it seems like it should consider the regular velocity since it is not guaranteed to be zero along the collision normal.

Am I thinking about this correctly? Does what I've done here look reasonable?

The bias impulse/velocity is totally separate from the regular velocity. So, you should calculate a relative bias velocity along the normal using only the bias linear and angular velocities.

I think it's easier to think about this in terms of the two types of error you're correcting: position and velocity.

Velocity error gets corrected via impulses as usual: with no restitution, there should be no bounce, because it should be zeroing the inward-along-normal velocity at the contact point and that's it. (If not, you might have a bug wrt how you accumulate/clamp the impulse, Erin Catto has several great presentations about this.)

You can use the same impulse-based approach to solve position errors ("bias velocity") but note that this doesn't work that well because during your iterations, your constraint gradients are fixed (since positions aren't changing). AFAIK Catto moved away from using bias and directly solves the position error (see SolvePositionConstraints() here https://github.com/erincatto/box2d/blob/main/src/dynamics/b2_contact_solver.cpp

This

sigh.. gamedev.net is such buggy garbage lately. Anyway what I was saying is just: correcting the pos/orn directly is equivalent to, after each bias-impulse iteration, applying the bias-impulses and then clearing them. But it's simpler to just directly update pos/orn, see line 735 of the file I linked to.

Solving position directly converges better since you're updating the constraint gradients at each iteration rather than keeping them fixed (as they are with impulses); of course this means you have to recalculate them at each iteration so it's a bit more expensive.

IMO you want to consider “solve velocity error” and “solve position error” to be two separate processes which don't occur simultaneously, you do one and then the other. Box2D is definitely a great reference for this. You might also be interested in Position-Based Dynamics which gets rid of impulses and just does everything at the position level – I think it's a lot simpler than impulses, and for joints it's much more stable/strong (the down side is that it makes collision detection more complex/expensive since things are moving around during iteration rather than being fixed in place).

@aressera @raigan
Okay I'll remove the regular velocity from the bias velocity calculation since they're separate things.

It's interesting that in Box2D, the positional correction is not done at the velocity level but rather, the positions and orientations are directly updated. I think my current solver is stable enough for now. Maybe I'll switch to that approach if I need more stability in the future given it's more expensive. So the extra calculations you'd have to do each iteration of the solver would be essentially be related to the effective mass and the inertia tensor? And you'd have to recalculate those because you'd actually be moving the object each iteration and those two things depend on the object's position/orientation right?

@CasualKyle the constraint gradient direction itself (ie the position of the contact point and direction of the contact normal) also change when you change the pos/orn.

@raigan Right, yes that makes sense. Thanks for the great information.