The key issue here is mathematical. We review some basic terms.
Kinematics is the study of motion. You can do a lot of simulation using nothing but the two basic definitions
Velocity vx = (dx/dt) the first derivative of position with respect to time. A vector V=(vx, vy, vz) of three similar terms.
Acceleration ax = (dvx/dt) the first derivative of velocity with respect to time, with a similar vector A.
- and their angular equivalents, which are a bit tricker in actual use.
Angular velocity w = (d(theta)/dt) where theta is an orientation angle. This works nicely if theta is oriented around an axis of rotation, but angular velocities don't add together in a nice orthogonal fashion like (x,y,z) linear velocities do. You're best off using quaternions to work with angular orientation and velocity.
Angular acceleration U = dw/dt, and this gets yet messier, as you realize when you play with gyroscopes. Let's leave that topic for a course on physics.
With just kinematic equations you can model bouncing balls, almost all 2d games, and lots of good-looking animations. Special purpose ("ad hoc") tricks are used to handle balls bouncing off walls; e. g. by just reversing the velocity vector in the direction whose coordinates exceed the limits imposed by the wall. But ultimately you need the Real Stuff - dynamics.
Dynamics is the study of force, momentum and energy. Actually the definition is more general than that, with energy being universal to all dynamical systems whereas force and momentum relate specifically to the motion of matter in space. In electronics, for instance, momentum's role is played by inductance and force corresponds to voltage or electromotive force.
We'll stick to spatial motion of massive objects and won't worry even about fluids, for now (though Yann's project involves fluids.) The Navier-Stokes equations are interesting if you want to go learn about them.
F = mA
is the key to it all, where F is a force (vector), m is a mass (scalar) and A is acceleration. There is an analogous equation for rotational energy.
Query 8.1. What is a stateless, or instantaneous model? Give a couple of examples.
ANSWER
It's a model whose output is computed instantaneously
from the input. Pressure of a gas in a closed container of fixed size,
varies instantaneously with temperature. Current in a circuit with fixed
resistance varies instantly with voltage.
/ANSWER
Queries 8.2 and 8.3 are difficult to answer separately, so you can combine your answers if you wish. Just use your own language, not Korn's.
Query 8.2. What is a state variable? How do state variables relate to degrees of freedom for describing a system? This second question is not discussed in the book chapter.
ANSWER
A state variable is the embodiment of a system's memory. It represents
some quantity or attribute which is persistent about a system.
The degrees of freedom of a system usually means how many parameters
are necessary to specify its behavior. These would consist of the effective
number of input variables, together with initial values for the state variables.
By "effective" input variables, I mean that one could always (deliberately
or accidentally) define input variables that control less DOFs than there
are variables. Some inputs might simply not be hooked up to anything. Or
one might allow a user to specify the x,y and z values of a directional
vector, when in fact only two DOFs are needed for a 3space directional
vector after its magnitude is ignored.
/ANSWER
Query 8.3. What is a state transition model? Give an example
not cited in the book.
ANSWER
A state transition model is a model which has memory in it. The output
depends not only on the instantaneous input, but on what happened leading
up to this instant. Consider a water tank, whose input device is a hose
with a gallons-per-second meter on it. Its output is a depth gauge reporting
the number of gallons in the tank. Just looking at the input flow-meter,
you cannot tell me what the gallon-meter will say. You must know its history;
e. g. "5 gal/minute * 60 minutes."
/ANSWER
Read Section 1.1 about Models, Time Histories and State Equations.
Query 8.4. How can a differential equation of order n be reduced to n first-order differential equations? What is a phase variable?
ANSWER
Just cite the text.
/ANSWER
Query 8.5. Consider a mass attached to a spring, in a nearly but not completely friction-free environment. The objects' resting position (when the spring is neither compressed nor stretched) is defined as position x=0. If a graph were drawn with the mass' velocity on the vertical axis and the x-position of the mass' center on the horizontal axis, what trajectory would be traced?
ANSWER
A spiral pattern centered on the origin.
/ANSWER
Query 8.6. Now consider a pendulum that is free to swing in two directions, which is set to swinging in a circle or an ellipse. (My granddaughter Emily does not like to be swung in this fashion!) How would you plot its phase trajectory, and what would it look like?
ANSWER
One way would be to draw a 3d picture with two horizontal axes for
position's x and y components, and the vertical one for speed (magnitude
of velocity). Now you would see a spiral as energy is exchanged between
x and y velocity, with a slow diminution of its radius as the pendulum's
speed drops. Its height is also dropping, so the potential energy is being
lost as well. You get a funnel shaped trajectory.
/ANSWER
Query 8.7. What is a defined variable? Think up an example. Explain the concept in terms of degrees of freedom of the system, if you can.
ANSWER
A defined variable expresses the degrees of freedom of the state and
input variables, without adding another true DOF. It is computed as a function
of the state variables and has no independence of its own. Kinetic energy
is purely a function of velocity and so is not a state variable, in a system
model where velocity is treated as a state variable.
/ANSWER
Please don't make fun of the author's description of 80286 based computers as 'remarkably powerful differential-equation solvers'. We will all be laughing at our wimpy Pentium 600's in 18 months.
A Key Sentence (page 9): "The simulation model is represented in the computer by the derivative procedure which repeatedly evaluates all defined-variable and derivative expressions, to update state-variable time histories." The system performs this updating of time histories by integrating the variable set.
Integration Methods are briefly previewed on Page 11: "multistep formulas" are essentially using curve-fitting techniques to "guess" where the time series x(i), x(i+1) etc. will go, based on previous values of the derivatives. "Runge-Kutta" formulas try to anticipate what the future derivatives will be and use them to update the integral value. We have another paper coming up on integration methods.
The rest of the chapter is about simulation languages, which are not very relevant to realtime simulation. But they are interesting in their own right. The CSSL-type languages take in differential equations and sort them out so that no symbol appears on the right hand side of an assignment before it has received a value by occurring on the left side of some other assignment.
The DESIRE program on page 17 and 18 is interesting. It combines an
interpretive component which sets up the experiment's initial values, specifies
number of passes through the loop, and a DYNAMIC component which is compiled
into an efficient simulation model (and sorts the rules as necessary. If
any loops occur, an error message is generated.) Then the control
system runs the model and plots the results. The main thing that seems
missing is that it doesn't say which integration method to use. On
page 20 the text says the second-orde Runge-Kutta is the default.
We assume the integration step h, and refer to t(n)=nh where n is the number of steps so far taken.
If the function f(t) being integrated is an explicit function of t (that is, given any t you can calculate f without working your way forward from t=0) then life is beautiful and we can do lots of magic with the fact that we can compute f(t+h) any time we want to. But usually f(t) depends on either an input you don't have a time t, or it depends on the system's state at time t. So, f(t+1) would depend on the state at time t+1 - but that's exactly what we're trying to compute.
Query 8.8. Examine figure 1.3. Using good old Euler integration, where dx/dt = f(t) for some function f, if f is rising at time (t) then can you always be sure that Euler integration will underestimate the integral of f from t to t+h?
ANSWER
No, because f(t) might take a nose dive during the integration interval.
We're just GUESSING that f(t+h) will be somewhere in the vicinity of f(t).
Could be higher or lower. So in a funny way, Eule is a better system with
large step sizes than with small ones. In small step sizes you can be pretty
sure that the estimate always lags the true integral, whereas with big
step sizes the function may have turned around - even several times - before
the next sample.
/ANSWER
Query 8.9. CHALLENGE PROBLEM. (A challenge problem is one that is tougher than what I will give on the final exam. If you're really serious about this stuff you'll give it a shot, but if you're skipping stuff to save time, this is one to skip.)
On page 1-7 is discussed the case where f is a linear function of x (and though they don't say it, f has no dependence on input u(t)). Such models are not often needful of realtime simulation, because if they don't depend on the inputs you might as well just pre-simulate them and store the results as an animation. But it's academically interesting to get you used to the trapezoidal integration idea.
Think up an example of a simple physical system where f would be a linear function of x. If you can't think of one, skip down one page to the large word HINT. In fact I'll give you a hint right here: set up the state variable X = (x,v) where v = dx/dt. Then do as they say: solve the equation for X(n+1) in terms of X(n). You'll know you have a solution when everything on the right side of the equations defining X(n+1), contains no reference to anything at time n+1 but only at time n.
What was the crucial step at which the linearity of f(X) mattered?
ANSWER
My example is the classic mass on a spring, where the restoring force
r = -kx for spring constant k. Define v=dx/dt so we'll have X=(x,v) as
our state vector. Then by Newton's famous law of force = mass * acceleration,
m * dv/dt = r = -kx
and of course
dx/dt = v
by definition. So our overall integration formula for trapezoidal integration (1.6) can be broken into two equations, for x and v. The 'f' which must be integrated to compute v, is -kx. The 'f' which must be integrated to compute x, is v. So we can refer to
F=(v, -kx)
and think of the whole thing at once, as
Xn+1 = Xn + h(Fn + Fn+1)/2
or we can do two separate equations for x and v, as follows.
I: xn+1 = xn + h (vn + vn+1)/2
II: vn+1 = vn + h (-kxn - kxn+1)/2
Now it should be clear that we have two linear equations in which xn+1 and vn+1 are the unknowns and everything else is known at time n. If we substitute xn+1 into the second equation, we find pretty quickly that vn+1 is on both sides of the equation and xn+1 is all gone. Solving this for vn+1, we get
vn+1 = (vn (1-(kh**2)/4) - hkxn) / (1-(kh**2)/4)
That looks worse that it is. vn+1 is now easily calculated from stuff that's known at time n, and we can then go back to equation I above and evaluate xn+1 as well. A quick plausibility check shows that if k or h are zero, representing a broken spring and a zero timestep respectively, then the velocity doesn't change and x keeps growing at constant velocity. You can't prove an equation valid by such simple tests but you can sometimes prove them wrong.
/ANSWER
But if you don't have such an easy situation with f, then much more work is required. If f is nonlinear in X but still doesn't depend on inputs, you can generally use iterative methods rather than simply solving a pair (or a larger set of) linear equations.
However in the wide wild world of realtime simulation, most of the time there are live inputs to deal with. So we move on to 'real' integration techniques.
Adams-Bashforth second order (AB-2.)
We might extrapolate fn-1 and fn to guess at fn+1. The integration formula is shown in 1.8. The explanation for this equation is shown in Figure 1.5 and is based on simple geometry. fn - fn-1 represents the height of the triangular portion of the grey area which is sitting on top of a rectangle of width h and height fn. So the triangle area formula of 1/2 * base * height is used.
This game can be continued by using older terms and higher order curves.
Query 8.10: Summarizing the Howe paper.
Make a table of the integral methods discussed in section 1.2. Your purpose in designing this table is to be able to look at it and select the most appropriate integrator for a particular situation in a game.What are your column heads? One of them should be "Compatible with realtime use (yes/no)."
The Lunar Lander Game
In the 1960's and 1970's, we were all very interested in the lunar missions. There was life-or-death drama hanging on correct engineering. So almost every computer scientist with access to any kind of a graphical display programmed or played with the Lunar Lander simulation. And now you will get do do it, too.
Basic facts: the Moon's gravitational acceleration is something like 1 meter per second per second at its surface. But we want the acceleration to vary with altitude, so we'll just use a constant G for the Moon's gravitational constant including the moon's mass. The rocket's non-fuel mass including astronauts is R; the fuel's mass which will decrease as the engine burns is S(t). Total mass M(t)=R+S(t)
The throttle settings range between 0 and 1 and are represented by the input function U(t). The user's task is to adjust the throttle so as to land with velocity of no greater magnitude than 2 meters/second (or the lander will be destroyed.).
Query 8.11: Identify the state variables and any derived variables; then write out the equations of motion for the lunar lander. Then write the recurrence relations for the state variables, under the assumptions of
a) Realtime RK-2 integration, and
b) AB-2 integration, and
c) AB-4 integration.
How would you propose to formally evaluate these three models, if you were charged with coming up with an objective measure of their relative validity for use in a particular realtime simulator? Just building and flying the simulators and evaluating them "by feel" wouldn't be adequate if you were doing the job for NASA.
ANSWER
All scalar variables in this writeup will be shown in lowercase; constants, the input function U and the state vector X will be in uppercase.
Let the thrust of the engine be K1*U(t). The total forces acting on the rocket are therefore
f = K1*U(t) - (K2*m(t)/x2)
where x is the height of the lander above the center of the moon; K1
is a constant representing the amount of thrust delivered when the throttle
is fully open (1.0), and K2 = G*L where G is the universal gravitational
constant and L is the mass of the moon.
The state variables are x, v (velocity) and s the mass of fuel remaining.
The derived variable is m, the total mass.
We now only need to write down derivative equations for the three state variables.
from Force=m*a = m* vdot
(where vdot means the first derivative of velocity v, this being the definition of acceleration)
we get
vdot = Force/m
so we can put it all together and write
vdot = (K1*U(t) - (K2*m(t)/x2)) / m(t)
xdot = v
sdot = K3*U(t)
where K3 is the rate fuel is burned when the throttle is fully open. We assume a linear relationship between throttle value and thrust, and also between throttle value and fuel consumption. These assumptions don't describe how engines really work, but are convenient for our project.
We also need the derived variable equation
m(t) = R + s(t)
Now - how do we use this information to set up the recursion relations for the three integration techniques? From the Hoppe paper we have the unfortunate convention that the derivative expression is called F(Xn, Un) but it's not a force. Also it has three components which we could call Fx, Fv and Fs, whose values are xdot, vdot and sdot as listed above. We have to choose the right one in each case.
The vector X = (x, v, s) is the state vector of the system.
We go back and forth between U(t) and Un whenever convenient, knowing that t=hn.
The X with the hat on it is pronounced "estimator for X" and so we'll type it as estX. We break this one equation out into three components in the following, and so we will be using estx, estv and ests.
Realtime RK-2:
estxn+1/2 = xn + (h/2)*Fx(Xn,Un) = xn + (h/2)*v -- because the derivative expression for x is just v
estvn+1/2 = vn + (h/2)*((K1*Un - (K2*m/x2)) / m) -- the derivative expression for v is messy
estsn+1/2 = sn + (h/2) * K3*Un -- the derivative expression for s is pretty simple.
estmn+1/2 = estsn+1/2 + R -- we need an estimator for m, too, for the next step.
--------------
xn+1 = xn + h* estvn+1/2 -- This is because we are now using the derivative expression of the ESTIMATORS
vn+1 = vn + h* Fv(estXn+1/2,Un+1/2) = vn + h* ((K1*Un+1/2*h) - (K2*estmn+1/2/x2)) / m)
-- So when we say "Fv(estXn+12, Un+1/2)" we really mean to "pick out the useful parts of the state vector X and use them as needed."
sn+1 = sn + h * K3*Un+1/2
So, you can see what one must do to translate the general concept of equations 1.21 and 1.22 into specific equations for an actual simulator.
AB-2 will be very similar but simpler, since we don't have to put in the predictor step.
xn+1 = xn + (h/2)*(3vn - vn-1)
vn+1 = vn + (h/2)*(3* [K1*Un - (K2*mn/xn2)) / mn] - (K1*Un-1 - (K2*mn-1/xn-12))
sn+1 = sn + (h/2) * (K3*(3*Un - Un-1))
AB-4 will be the same technique applied to the messy expression in 1.12. Let's not write it out!