damage.md

damage.md

I have a system where something akin to damage occurs over time. The damage goes from level 0 to level 10 in a type of scurve (starting at 0-ish at time=0, and ending at 10-ish after some period). This scurve is slightly stretched and different at different intensities. The whole surface can be modelled with a smooth function using some kind of made up extruded scurve. One simple example being:

def damage_func(duration_at_constant_intensity, intensity, *optimal_parameters):
   return extruded_scurve(duration_at_constant_intensity, intensity, *optimal_parameters) 

With an example of an extruded scurve:

def scurve(x, stretch, mid, peak=10):
    return peak - (peak / (1 + np.exp(-(1 / stretch) * (x - mid))))

def extruded_scurve(x, y, a_stretch, b_stretch, a_mid, b_mid):
    stretch = a_stretch * x + b_stretch
    stretch = np.where(stretch > 0, stretch, 0)
    mid = a_mid * x + b_mid
    return scurve(y, stretch, mid)

The optimal_parameters (in this case a_stretch, b_stretch, a_mid, and b_mid) are found by optimizing real life data to fit this extruded scurve. This real life data is the result of multiple trials where the intensity is kept constant and the damage data points are the resulting damage over a period of time.

The problem we want to tackle is to calculated the running damage for varying intensities. We assume damage is the results of something positively cumulative (like a damage derivative). We envision the underlying damage to behave as follows. Let's say you move from one intensity y[n] to a different intensity y[n+1], with an established damage level of z[n], at a duration of x[n], the damage accumulated from there on is then the same as starting the damage curve at a duration_at_constant_intensity where z[n] = damage_func(duration_at_constant_intensity, y[n+1]) is solved for duration_at_constant_intensity. Note that duration_at_constant_intensity will not be follow the same progression as the true duration x.

With this in mind, to evaluate the unknown damage curve for a given duration and intensity line, we can do a kind of iterative integration. Our current approach is similar to the following.

• Find an inverse function g for z = damage_func(x, y) so that x = g(y, z)
• Find the partial derivative a = (∂damage_func/∂x)(x, y) and call it a = h(x, y). Some kind of accumulation of a will then result in a total damage figure.
• Combine these two functions to find the relationship a = h(g(y, z), y) = a = p(z, y)

It is then easy enough to find a numerical integral estimate over p(z,y). In order to find z for my given x=duration and y=intensity (say x = [0, 0.1, 0.2, ...., 15] and varying y), this involves a loop, starting at z[0] = damage_func(x[0], y[0]) and continuing with z[n] = z[n-1] + p(z[n-1], y[n])*(x[n]-x[n-1]).

Here is a snippet of our current implementation:

d_damage_from_duration_and_intensity = ddx(damage_func)

def iterative_line_integral(intensity, duration, optimal_params):

    delta_duration = duration[1:] - duration[:-1]

    if not np.all(np.abs(delta_duration - delta_duration[0]) < 1e-6):
        raise ValueError(
            "Current implementation requires all time-steps to be of equal length"
        )

    delta_duration = delta_duration[0]

    damages = []
    damage = damage_func(duration[0], intensity[0], *optimal_params)
    damages = [damage]

    dur = duration[0]

    for (i1, i2) in zip(intensity[:-1], intensity[1:]):
        dur = minimize(lambda x: (damage_func(x, i1, *optimal_params) - damage) ** 2, dur).x[0]
        dmid = dur + delta_duration / 2

        damage = (
            damage
            + d_damage_from_duration_and_intensity(dmid, (i1 + i2) / 2, *optimal_params)
            * delta_duration
        )

        damages.append(damage)

    return damages