Avoiding fixed timestep death spiral?

[zorg]

I'm guessing you know the difference between the following two, right?

Code: Select all

local ta = 0.0 -- tick accumulator
local tr = 1/60 -- tick rate (seconds)
function love.update(dt)
    ta = ta + dt
    -- first method
    if ta >= tr then
        -- one full tick elapsed, do stuff using tickrate
        ta = ta - tr
    end
end

Even if the elapsed time (dt) passed is more than the tickrate you choose, if it's an intermittent issue, then the hiccup will dissipate down the line, in an amortized fashion: 0.0 -> 0.005 -> 0.017 -> one tick, subtract tickrate, remainder: 0.0001 -> 0.1 (hiccup) -> one tick, remainder: 0.082 -> 0.0825 -> still one tick, subtract, remainder: 0.061 -> ...
If, however your updates take constantly too long, then this will spiral out of control in the sense that the ticks will be delayed, and the logic will run slower.

Code: Select all

local ta = 0.0 -- tick accumulator
local tr = 1/60 -- tick rate (seconds)
function love.update(dt)
    ta = ta + dt
    -- secondmethod
    while ta >= tr do
        -- any number of full ticks elapsed, do stuff using tickrate
        ta = ta - tr
    end
end

This one's difference is that you're dealing with the problem at once, and you don't wait a whole cycle; this of course means, that some ticks probably won't be shown graphically, due to them being overwritten by the last cycle in that while loop. This is kinda pointless though.

If anything, i'd do this one:

Code: Select all

local ta = 0.0 -- tick accumulator
local tr = 1/60 -- tick rate (seconds)
function love.update(dt)
    ta = ta + dt
    -- third method
    local mult = 0
    while ta >= tr do
        mult = mult + 1
        ta = ta - tr
    end
    if mult > 0 then
       -- do stuff once, but pass everything the tickrate * mult value, instead of just the tickrate.
    end
end

Technically "tick-skipping" (although nothing would be skipped logic-wise), since you update the logic by a multiple of your tickrate, but the drawing isn't really affected (in the sense that it'll be called normally, one update - one draw; but it may present some graphical glitches due to the modified rate the update calculated with; also, whether box2D will like this or not, i don't know; it's technically a fixed rate, just multiplied by an integer... but if it doens't, then it'd be easy to loop only that part times the mult variable, and you'd still save on processing all other logic more times than needed (i.e. only once)

[airstruck]

Right, it looks like we came to pretty much the same conclusion, except you use highest possible multiple of the timestep where I'm just doubling it as many times as possible; not much difference there. I also went for a timestep that was slightly longer than 60 fps and has an exact binary representation, just to be safe, probably no real difference there either. Have you seen an approach like this described anywhere else?

[zorg]

Have you seen an approach like this described anywhere else?

Other than the obligatory one and a few others, haven't seen too many other timestep implementations, but this seemed like the most logical thing to do; i mean, going from the second to the third of my examples, at least.

[airstruck]

I've read that article; it looks more like your second example plus a capped dt. I haven't seen anything like the example I posted at the top of this thread or your third example anywhere, but it seems like it will fit my purposes better. It's just a little worrying that it doesn't seem to be described anywhere. I'm going to go ahead and use it; I'll bump this thread later if there are unanticipated problems with it.

[raidho36]

If you don't discard the time in some way then after its no longer slowed it'll speed up and perform more updates per second until it catches up. You can discard any multiple of timestep but only do one step. Timeline will fall behind, but unless it's an online game that's hardly relevant. I would suggest you postpone some non critical calculations and only do them if there is enough time left. Even better, you can do sensitive calculations every frame, but only do other stuff (like animations) when there is nothing else to compute.

The unity engine uses the catch up method of fixed timestep handling - if more time elapsed since last update than a full timestep, it will perform more steps until it would catch up to real time. Note that it has two different update functions - one with fixed timestep and the other with arbitrary. Keep in mind though that it runs fixed timestep on top of arbitrary, not separately, which is a design oversight. Overall it is the kind of implementation that very well can spiral into extremely low frame rate if fixed step computations are too heavy.

[airstruck]

I would suggest you postpone some non critical calculations and only do them if there is enough time left. Even better, you can do sensitive calculations every frame, but only do other stuff (like animations) when there is nothing else to compute.

I think it's going to have to be something similar to this. Specifically, things that do any kind of integration, from box2d to any kind of acceleration -> velocity -> position scheme, really need a fixed timestep for stability. There's really no way around it, for example you can't get the same result as a normal-length timestep by doing something like splitting acceleration -> velocity into two iterations for a double-length timestep unless you alternated between that and the velocity -> position integration, which would basically be the same as using a fixed timestep anyway.

Most other things probably don't need a fixed timestep, and can just do stuff whenever they have time.

[airstruck]

Hmm, actually the double-length timestep might work if multiple levels of integration happened all in one step. Instead of something like this:

Code: Select all

pos = pos + vel * dt
vel = vel + acc * dt

Something like this should do it for acceleration -> velocity -> position:

Code: Select all

pos = pos + acc * dt ^ 2 / 2 + vel * dt
vel = vel + acc * dt

Not sure if I got that right, not tested yet, but it seems like the basic idea should hold up. Of course this would become more complicated if higher level derivatives need to be modeled (for example "change in acceleration" all the way down to position). I suppose damping or anything else that would affect the outcome would also need to be part of the formula.

Does anyone have any experience with this kind of stuff? Too bad pgimeno hasn't been around, he was great at this stuff.

Edit: Apparently the 3rd and 4th derivative are called "jerk" and "snap." Not sure anything more than that would usually be necessary. I'm thinking if everything in the model is calculated "exactly," timestep won't matter so much for stability.

Here's what I've got so far; these are just notes and haven't been tested but should get the idea across. Would love to hear about it if anyone here has experimented with this stuff before.

Code: Select all

local function i1 (dx1, x, dt)
    return dx1 * dt + x
end

local function i2 (dx2, dx1, x, dt)
    return dx2 * dt ^ 2 / 2 + i1(dx1, x, dt)
end

local function i3 (dx3, dx2, dx1, x, dt)
    return dx3 * dt ^ 3 / 6 + i2(dx2, dx1, x, dt)
end

local function i4 (dx4, dx3, dx2, dx1, x, dt)
    return dx4 * dt ^ 4 / 24 + i3(dx3, dx2, dx1, x, dt)
end

local s = t.snap or 0
local j = t.jerk or 0
local a = t.acceleration or 0
local v = t.velocity or 0
local p = t.position or 0

if s ~= 0 then
    t.position = i4(s, j, a, v, p, dt)
    t.velocity = i3(s, j, a, v, dt)
    t.acceleration = i2(s, j, a, dt)
    t.jerk = i1(s, j, dt)
elseif j ~= 0 then
    t.position = i3(j, a, v, p, dt)
    t.velocity = i2(j, a, v, dt)
    t.acceleration = i1(j, a, dt)
elseif a ~= 0 then
    t.position = i2(a, v, p, dt)
    t.velocity = i1(a, v, dt)
elseif v ~= 0 then
    t.position = i1(v, p, dt)
end

[raidho36]

It's important to keep in mind that this whole thing is computing an integral of a polynomial function in real time with "rectangle" method. It's a reasonably good approximation but it's stability hinges upon sampling frequency, with discrepancy magnitude going to the factor of 2.71828 between sampling reaching infinity and a single sample across entire function, and that's for simple linear function. The simulation will actually behave differently if you use different timestep, which is why Box2D becomes unstable if you feed it variable dt. So while math behind it is sound, it will not behave very nicely.

[airstruck]

it's stability hinges upon sampling frequency, with discrepancy magnitude going to the factor of 2.71828 between sampling reaching infinity and a single sample across entire function

There is some instability, but I thought it might be due to floating point imprecision. At least some of it is due to that, anyway, because it's definitely more stable when all numbers involved are representable in binary (for example timesteps of 1/64 and 1/32 instead of 1/60 and 1/30). It can actually run for a (simulated) hour and not lose any precision at all with certain initial values, but usually after an hour the single- and double-step will have deviated by something like 0.001 (which still seems very stable and ought to be good enough for stuff like demo playback).

Not quite sure what the relevance of Euler's constant is to stability, but I'm not that familiar with this stuff. I thought using polynomials to solve the higher order integrations would get rid of the problem with sampling error (it looks like it has, mostly). Where can I read more about that?

Here's another code dump; this one's actually tested and I think works as intended if anyone wants to play around with it.

Code: Select all

local factorial = setmetatable({ 1, 2, 6, 24, 120 },
    { __index = function (t, i) t[i] = i * t[i - 1] return t[i] end })

local function integrate (dt, d, ...)
    local n = select('#', ...)
    if n == 1 then return d * dt + (...) end
    return d * dt ^ n / factorial[n] + integrate(dt, ...)
end

local function Integrator (keys)
    local vars = {}
    for i = 1, #keys do
        vars[#vars + 1] = ('local %s = t.%s or 0'):format(keys[i], keys[i])
    end
    local lines = {}
    for i = 1, #keys - 1 do
        lines[#lines + 1] = ('elseif %s ~= 0 then'):format(keys[i])
        for j = #keys, i + 1, -1 do
            lines[#lines + 1] = ('t.%s = _integrate(dt, %s)'):format(keys[j],
                table.concat(keys, ',', i, j))
        end
    end
    local source = ([[
        local _integrate = ...
        return function (t, dt) %s if false then %s end end
        ]]):format(table.concat(vars, '\n'), table.concat(lines, '\n'))
    return (loadstring or load)(source)(integrate)
end

local motion = {
    Integrator { 'snapX', 'jerkX', 'accelerationX', 'velocityX', 'x' },
    Integrator { 'snapY', 'jerkY', 'accelerationY', 'velocityY', 'y' },
    Integrator { 'angularSnap', 'angularJerk', 'angularAcceleration',
        'angularVelocity', 'angle' },
    Integrator { 'scaleSnap', 'scaleJerk', 'scaleAcceleration',
        'scaleVelocity', 'scale' },
    apply = function (self, t, dt)
        for i = 1, #self do self[i](t, dt) end
    end
}

local single = 1/64
local double = 1/32
local t1 = { jerkX = -0.125, accelerationX = 4000 }
local t2 = { jerkX = -0.125, accelerationX = 4000 }

for i = 1, 32 * 60 * 60 do
    motion:apply(t1, single)
    motion:apply(t1, single)

    motion:apply(t2, double)
end

for k, v in pairs(t1) do
    print(('%s\t%s | %.17G | %.17G')
        :format(v == t2[k] and ':)' or ':(', k, v, t2[k]))
end

Of course this also might be way too heavy for anything practical, even if stability were perfect; still need to look into that.

[raidho36]

No, by instability I mean the insatbility of the simulation, i.e. identical starting conditions will take you to different end point. This instability is due to rectangle method's inherent behavior. Perhaps I could've used the "chaotic" term to describe it - I didn't because that would be a misnomer.

By carrying out simulation you compute integral of polynomial defined by speed, acceleration, etc. Area of the function is the distance traveled. The way all of that usually simulated, speed is another integral computed with rectangle method over acceleration function, acceleration is another integral and so on. Errors pile up very fast and the only way the system behaves consistently if it's sampling rate is consistent.

The following image is an integral of a function computed with rectangle method. Area is distance traveled, the function is speed. Discrepancy between true value and computed value is the difference between total area of blocks and total area of the function. It's not hard to imagine that area of blocks (distance traveled, speed gained, etc.) will change substantially if you increase or decrease sampling frequency.

integral.png (2.62 KiB) Viewed 4072 times

Additionally, computation errors have magnitude of 10^-15 for doubles, 0.001 is 13 orders of magnitude above basic error, it would take 10 quadrillion cycles to accumulate this much error. And that's assuming exponent of accumulator never goes below 0, in which case it would take even more cycles. So that's not the source of your timing problem.