Sign in to follow this  

Best way to make a physics object hover

This topic is 4106 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hi, If I'm using a physics engine whats the best way to make an object hover? I've looked all over google and its not much help - just abstracted heavily mathematical theory on hovercrafts which isn't useful and such. The approach I've tried is to get the height and at every step, apply an upthrust obeying the inverse square law, assuming that the height will stabilise out. The problem is, when it hits the ground it will suddenly jump high up (K/0.01*0.01 is a large number), so I used clipping. Then another problem is that there won't be enough force above a certain height, so it will float gently down then jump up, when the floors upthrust tips the critical force to outweigh the gravity. Then I tried switching off gravity - no luck. I also tried more force when its falling, which seems a mediocre fix- the object will jut around at a fixed height as it re-adjusts itself. I wouldn't mind if it slowly cycled its height around a certain height, but its the slow falling and sudden jump upward as opposed to moving like a sine wave around its height. I've tried all kinds of laws like exponential, linear .etc And the working solution is very hackish- can anyone impart any advice? If it helps I'm using ODE. Thanks for your time.

Share this post


Link to post
Share on other sites
In these kind fo systems (let it be A.I. sensors, or how much a A.I. car should apply force to correct it's steering), I usually look 'forward' in time, instead of the instantaneous inputs.

Say you use a sort of spring, like you do.

F = k / l(t)^2

instead of looking at l(t), I try to find out what the length will be after say, one second, or 1/2 second, or whatever looks good.

So if the craft is falling looking ahead in time, then you will have to apply a strong force, and has the equilibrium point is reached, the force will counteract the movement gently and naturally.

It's kind of a damping system if you will. There could be a slight sway, but as you said, that's not a problem, and that's tweakable depending on how much you look ahead in time. It's kind of a negative feedback system.

Share this post


Link to post
Share on other sites
Hmm, well, I've implemented a hovercraft (air cushion vehicle) simulation, but you say you aren't interested in the heavy math aspect of it.

Unfortunately, for anything that is realtime, interactive, and true 3D, you are going to have an extremely difficult time implementing this without the difficult theory!

I honestly don't believe that looking ahead in time, or any other time of coherence implementation, will work, except for in extremely idealized situations (i.e. sometime that is limited in what you can do with the hovercraft, or something that isn't particularly interactive).

I implement the air cushion by defining a certain number of sample points arranged in a circular fashion underneath the chassis of the hovercraft (the skirt and plenum chamber area). I run an RK4 solver to predict the new position and orientation, based on the inputs from the air cushion. Each sample point behaves like a spring, but a more detailed fluid flow simulation could be used instead of basic distance finding (which would require more sample points; finite element analysis).

Here's the force and torque function (the 'derivative' function when reading about numerical methods of problem solving). For each sample point, it measures the distance to the ground and applies a force at that point based on the relative velocity between the sample point and the ground, and a dampening factor is also added in there. The force per sample point is the total weight of the hovercraft, divided by the number of thrusters, at a given height (like a spring).

Here's this function:


void Hovertank::Force_and_Torque(Vector3 *output_force,Vector3 *output_torque,Vector3 & position,Quaternion & quat_orientation,
Vector3 & linear_momentum,Vector3 & angular_momentum,Matrix3x3 & mat_orientation,
Vector3 & linear_velocity,Vector3 & angular_velocity)
{
*output_force += Vector3(0.0f,-gpPhysicsManager->mGravAccel*mMass,0.0f); //applies weight

Air_Cushion transformed_cushion = mAirCushion;
Matrix4x4 cushion_sample_transformation;
cushion_sample_transformation.SetTranslation(&position);
cushion_sample_transformation = mat_orientation;//CurrentState.mMat_Orientation;

for(int i = 0; i < transformed_cushion.num_sample_points; i++)
{
transformed_cushion[i] = cushion_sample_transformation * transformed_cushion[i];
}

transformed_cushion.orientation = mat_orientation * transformed_cushion.orientation;

Vector3 start = position;
Vector3 end = position + (transformed_cushion.orientation * transformed_cushion.rest_height);

#if! PRE_FIND_TOUCHED_LEAFS
std::vector<int> leafs;
MoveTouchNode(gpBSP,0,&start,&end,mRadius,leafs); //collect leafs
#endif

Vector3 total_cushion_force(0.0f,0.0f,0.0f);
Vector3 total_cushion_moment(0.0f,0.0f,0.0f);

for(i = 0; i < transformed_cushion.num_sample_points; i++)
{
end = transformed_cushion[i] + (transformed_cushion.orientation *
(transformed_cushion.rest_height)); //FIXME: transformed_cushion.rest_height - radius ???

Vector3 pos = transformed_cushion[i];

HitInfo CGTrace =
#if PRE_FIND_TOUCHED_LEAFS
TraceSphereToBSPFaces(gpBSP,transformed_cushion[i],end,transformed_cushion.sample_radius,gpPhysicsManager->physics_epsilon,mvBSPLeafs);
#else
TraceSphereToBSPFaces(gpBSP,transformed_cushion[i],end,transformed_cushion.sample_radius,gpPhysicsManager->physics_epsilon,leafs);
#endif
if(CGTrace.Time < 1.0f)
{
Vector3 local_pos = transformed_cushion[i] - position; //reflects rotation only
#if USE_REAL_DIST
float dist_to_point = (CGTrace.HitPoint - transformed_cushion[i]).GetLength();
#else
float dist_to_point = (DotProduct(&CGTrace.HitPoint,&transformed_cushion.orientation) -
DotProduct(&transformed_cushion[i],&transformed_cushion.orientation));
#endif
float compression = fabs(dist_to_point - transformed_cushion.rest_height);

float basic_force_magnitude = compression * transformed_cushion.weight_per_thruster;


if(basic_force_magnitude > transformed_cushion.max_sample_force)
basic_force_magnitude = transformed_cushion.max_sample_force;

Vector3 sample_force = transformed_cushion.orientation * basic_force_magnitude * -1.0f;

float total_rel_vel(0.0f),rel_vel_fromlinear(0.0f),rel_vel_fromangular(0.0f);

rel_vel_fromlinear = DotProduct(&transformed_cushion.orientation,&linear_velocity);
rel_vel_fromangular = DotProduct(&transformed_cushion.orientation,
&CrossProduct(&angular_velocity,&local_pos));

total_rel_vel = rel_vel_fromlinear + rel_vel_fromangular;

sample_force += transformed_cushion.orientation * total_rel_vel * -transformed_cushion.dampen_coefficient;

total_cushion_force += sample_force;
total_cushion_moment += CrossProduct(&local_pos,&sample_force);
}
}

*output_force += total_cushion_force;
*output_torque += total_cushion_moment;
}



I dont' know how helpful this is, but seeing as how I've already implemented this it may help to give some ideas.


Share this post


Link to post
Share on other sites
Quote:
Original post by oliii
In these kind fo systems (let it be A.I. sensors, or how much a A.I. car should apply force to correct it's steering), I usually look 'forward' in time, instead of the instantaneous inputs.


What you're describing is known as a "PID controller". Genjix, Googling for that term should turn up plenty of info on what is a very common problem.

Share this post


Link to post
Share on other sites
tthibault:

Didn't say I was interested in heavy math. Just real-life math I cannot simulate ;)

Yours is a very interesting code which I'll have to check out.

Sneftel and oliii:

Looking from the tutorials this looks great. I found this on the net:
typedef struct
{
double dState; // Last position input
double iState; // Integrator state
double iMax, iMin;
// Maximum and minimum allowable integrator state
double iGain, // integral gain
pGain, // proportional gain
dGain; // derivative gain
} SPid;
double UpdatePID(SPid * pid, double error, double position)
{
double pTerm,
dTerm, iTerm;
pTerm = pid->pGain * error;
// calculate the proportional term
// calculate the integral state with appropriate limiting
pid->iState += error;
if (pid->iState > pid->iMax) pid->iState = pid->iMax;
else if (pid->iState
<
pid->iMin) pid->iState = pid->iMin;
iTerm = pid->iGain * iState; // calculate the integral term
dTerm = pid->dGain * (position - pid->dState);
pid->dState = position;
return pTerm + iTerm - dTerm;
}


I assume in my context error would mean the difference in height between my position and the desired position?

Would this roughly translate to:
float celPcHover::UpthrustForce ()
{
float pval, dval, ival;
float error = fabs (WorldToLocal (position - end_pos). Length ())
// calculate the proportional term
pval = error * prop_constant;
// calculate the integral term
istate += error;
ival = istate * int_constant;
// calculate the differential term
dval = pcmechobj->GetLinearVelocity () * vel_constant;
return pval + ival + dval;
}


Am I correct in my reasoning? Would the integrator part be the total amount of excess above the desired value over the past?

Thanks.

Share this post


Link to post
Share on other sites
Quote:

Why do you do this?

cushion_sample_transformation.SetTranslation(&position);

cushion_sample_transformation = mat_orientation;//CurrentState.mMat_Orientation;


Then you loop through all the sample points and transform them by that matrix. And then you again transform the cushion itself (why?)- wouldn't that be a double transform on each of the points?


I've defined the air cushion to have the following properties:
-a bunch of sample points laid out circularly *in object space*.
-An orientation (which points downward).

I realize now there are a few things I've left out which might confuse you, namely that the index operator of the air cushion structure returns the position of the indexed sample point, e.g.:

air_cushion[i] returns the position of sample point at index i.

In order to use the air cushion sample points they must be rotated out to reflect the translation and orientaiton of the hovercraft ( if the hovercraft starts to bank at an angle of 30 degrees, the position of the sample points rotates along with the chassis of the hovercraft). The direction that the sample point fires in is simply rotated.

>>Is weight_per_thruster = weight / n_thruster ?

Yes

>>Then it looks like you clip the force. Whats sample_force ?

Sample force is the force contribution of the current sample force.

Quote:

Then you subtract sample_force from our earlier calculated force. Why is this?

transformed_cushion.orientation * total_rel_vel * - transformed_cushion.dampen_coefficient


Hmm, I'm not entirely sure what you mean by the first part of this question, but the second line is how you model a dampened spring. According to hook's law, a spring exerts a force somewhat proportional to a constant times the compression length (in this case, how far the sample point is being compressed from the rest height). In order to dampen the spring so that the system loses energy and can come to a halt, an added force must be added that acts against the relative velocity (hence the negative). dampen_coefficient is basically just an ad hoc number that works. The larger the dampen coefficient is, the more quickly the system losts energy...the lower the value, the more the hovercraft will 'bounce' (which looks cool).

>>Sorry for all the questions, but this is the first promising thing I've found for weeks. Are you using ODE per chan

Don't apologize for asking questions!!! It took me about a month to get this to actually work. I don't use ODE, I try to write everything by myself (but what I've implemented is far less than any physics library).

Note that I implemented the force_and_torque function, but was using it with an euler solver, and therefore the logic was right but the hovercraft would eventually topple over due to numerical instability. That's what actually took me the longest to get working properly. Here's the RK4 solver that calls force_and_torque. This is essentially taken from Game Physics by David Eberly. Note that force_and_torque is a virtual function defined by the rigid body, and so you can write a new force_and_torque function to simulate cars, boats, airplanes, etc without having to change the structure of the RK4Update function.

If you approach the theory in the correct direction an RK4 solver is actually quite intuitive, but the end result certainly doesn't look it! You should start off with reading about an RK2 solver (which is essentially an improved euler method) and ignore matrices of partial derivatives (they only come into play when you take into account the initial time at the frame, I simply take into account the change in time).



//Only need to normalize quaternions at the last step
void RigidBody::RK4Update(double DT)
{
//part1
double halfdt = 0.5f * DT, sixthdt = DT / 6.0f;


Vector3 XN, PN, LN, VN, WN;
Quaternion QN;
Matrix3x3 RN;

/*
//Step1
*/

Vector3 A1DXDT = CurrentState.mLinearVelocity;
Quaternion A1DQDT = (Quaternion(CurrentState.mAngularVelocity) * CurrentState.mQuat_Orientation) * 0.5f;

Vector3 A1DPDT(0.0f,0.0f,0.0f);
Vector3 A1DLDT(0.0f,0.0f,0.0f);

Force_and_Torque(&A1DPDT,&A1DLDT,CurrentState.mPosition,CurrentState.mQuat_Orientation,CurrentState.mLinearMomentum,
CurrentState.mAngularMomentum,CurrentState.mMat_Orientation,CurrentState.mLinearVelocity,CurrentState.mAngularVelocity);


XN = CurrentState.mPosition + (A1DXDT * halfdt);
QN = CurrentState.mQuat_Orientation + (A1DQDT * halfdt);

PN = CurrentState.mLinearMomentum + (A1DPDT * halfdt);
LN = CurrentState.mAngularMomentum + (A1DLDT * halfdt);
Convert(QN,PN,LN,RN,VN,WN);

/*
//step 2
*/

Vector3 A2DXDT = VN;
Quaternion A2DQDT = (Quaternion(WN) * QN) * .5f;


Vector3 A2DPDT(0.0f,0.0f,0.0f);
Vector3 A2DLDT(0.0f,0.0f,0.0f);

Force_and_Torque(&A2DPDT,&A2DLDT,XN,QN,PN,LN,RN,VN,WN);


XN = CurrentState.mPosition + (A2DXDT * halfdt);
QN = CurrentState.mQuat_Orientation + (A2DQDT * halfdt);

PN = CurrentState.mLinearMomentum + (A2DPDT * halfdt);
LN = CurrentState.mAngularMomentum + (A2DLDT * halfdt);
Convert(QN,PN,LN,RN,VN,WN);

/*
//Step 3
*/

Vector3 A3DXDT = VN;
Quaternion A3DQDT = (Quaternion(WN) * QN) * .5f;


Vector3 A3DPDT(0.0f,0.0f,0.0f);
Vector3 A3DLDT(0.0f,0.0f,0.0f);

Force_and_Torque(&A3DPDT,&A3DLDT,XN,QN,PN,LN,RN,VN,WN);

XN = CurrentState.mPosition + (A3DXDT * DT);
QN = CurrentState.mQuat_Orientation + (A3DQDT * DT);

PN = CurrentState.mLinearMomentum + (A3DPDT * DT);
LN = CurrentState.mAngularMomentum + (A3DLDT * DT);
Convert(QN,PN,LN,RN,VN,WN);

/*
//Step 4
*/


Vector3 A4DXDT = VN;
Quaternion A4DQDT = (Quaternion(WN) * QN) * 0.5f;

Vector3 A4DPDT,A4DLDT;

PredictedStateVector.mPosition = CurrentState.mPosition + (A1DXDT + ((A2DXDT + A3DXDT) * 2.0f) + A4DXDT) * sixthdt;
PredictedStateVector.mQuat_Orientation = CurrentState.mQuat_Orientation + (A1DQDT + ((A2DQDT + A3DQDT) * 2.0f) + A4DQDT) * sixthdt;

PredictedStateVector.mQuat_Orientation.Normalize();

CurrentState.mLinearMomentum = CurrentState.mLinearMomentum + (A1DPDT + ((A2DPDT + A3DPDT) * 2.0f) + A4DPDT) * sixthdt;
CurrentState.mAngularMomentum = CurrentState.mAngularMomentum + (A1DLDT + ((A2DLDT + A3DLDT) * 2.0f) + A4DLDT) * sixthdt;


Convert(PredictedStateVector.mQuat_Orientation,CurrentState.mLinearMomentum,CurrentState.mAngularMomentum,
PredictedStateVector.mMat_Orientation,CurrentState.mLinearVelocity,CurrentState.mAngularVelocity);

PredictedStateVector.mMat_Orientation = PredictedStateVector.mQuat_Orientation.ToMatrix();
PredictedStateVector.mMat_Orientation.ToRadians(&PredictedStateVector.mEuler_Orientation.x,
&PredictedStateVector.mEuler_Orientation.y,&PredictedStateVector.mEuler_Orientation.z);

mJInverseGlobal = CurrentState.mMat_Orientation * mJBodyInverse * CurrentState.mMat_Orientation.Transpose();

//FIXME: this only works when the hulls are always spheres!
CurrentState.mEuler_Orientation = PredictedStateVector.mEuler_Orientation;
CurrentState.mMat_Orientation = PredictedStateVector.mMat_Orientation;
CurrentState.mQuat_Orientation = PredictedStateVector.mQuat_Orientation;
}





Note that this is insanely expensive. The code is geared towards accuracy, not speed. An lower order solver would probably be a better balance. I also apologize for not being particularly eloquent, so feel free to ask me even more questions.

Share this post


Link to post
Share on other sites
Quote:
Original post by Sneftel
Quote:
Original post by oliii
In these kind fo systems (let it be A.I. sensors, or how much a A.I. car should apply force to correct it's steering), I usually look 'forward' in time, instead of the instantaneous inputs.


What you're describing is known as a "PID controller". Genjix, Googling for that term should turn up plenty of info on what is a very common problem.


Ah cool, I'll try to remember the name for the next time [grin]

Share this post


Link to post
Share on other sites

This topic is 4106 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this