Max function in differential equation

[Ge0rges]

Hello all,

I have a set of differential equations I’m trying to solve. One of them includes a max() function. I am wondering if that is the reason why I am sometimes unable to solve this equation. Can the solver “integrate” max()? Do I need to specifically handle that discontinuity somehow?

Here is the model code:

# Load state conditions
    pOC_content = u[1]
    dOC_content = u[2]
    inorganic_carbon_content = u[3]
    cell_count = u[4]

    ## CELL COUNT
    # Growth
    growth = mu_max * (dOC_content / (Ks + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

    # Organic carbon requirement
    required_dOC_per_cell = maintenance_per_cell
    required_dOC = required_dOC_per_cell * cell_count

    # Starvation Deaths
    dOC_missing = max(required_dOC - dOC_content, 0)
    starvation_deaths = dOC_missing == 0 ? 0 : dOC_missing / required_dOC_per_cell

    # Total Deaths
    deaths =  starvation_deaths

    ## CARBON
    dOC_consumption = required_dOC_per_cell * (cell_count - deaths)
    fixed_carbon = inorganic_carbon_fixing_rate * inorganic_carbon_content

    # Particulate Organic carbon
    du[1] = pOC_input_rate - eea_rate*cell_count

    # Dissolved Organic carbon
    du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count + fixed_carbon

    # Inorganic carbon
    du[3] = inorganic_carbon_input_rate + inorganic_carbon_per_cell * (deaths - growth) + required_dOC - fixed_carbon

    # Net cell count change
    du[4] = growth - deaths

[Oscar_Smith]

the discontinuity is a problem. check out the documentation for callbacks.

[ChrisRackauckas]

Handling Instability When Solving ODE Problems

Thanks for the example. I will share as little as possible but just a snippet in order to explain what’s going on to others. This is a good teaching opportunity. I was able to find out that the problem is that, in the statement k9*X3/(k10*k11/(k12*(1.-min(1.,k13+k14*X5))) + X3) when X3 and X5 are dual numbers, the value of k11/(k12*(1.-min(1.,k13+k14*X5))) is NaN for the derivative components when the condition makes it the constant 1.0 (since the derivative of the minimum of in the denominat…

[Ge0rges]

Thanks for that. I started by trying a soft max approach as follows.

    # Define a soft max() function.
    function softmax(x, y, k)
        return log(exp(k*x) + exp(k*y))/k
    end
    dOC_missing = softmax(required_dOC - dOC_content, 0, 10) # max(required_dOC - dOC_content, 0)

However the yielded instability in previously stable solutions, even for larger or smaller K (1, 10, 100). I then opted to implement the callback solution, however this runs out of iterations only in this case.

    ...
    condition(u, t, integrator) = p[2] * u[4] - u[2]
    affect!(integrator) = nothing
    max_cb = ContinuousCallback(condition, affect!, save_positions=(false,false))

    # Callback list
    cbs = CallbackSet(carbon_add_cb, max_cb)

    # Build the ODE Problem and Solve
    sol = solve(prob, Rosenbrock23(), callback=cbs, tstops=stops, isoutofdomain=is_invalid_domain, maxiters=1e7)

Tried with autodiff=False and with abstol=1e-10, reltol=1e-10 as well.

Given that my max is max(X, 0), the condition is simply X. I believe the reason this is happening is because the values in X approach equality but not exactly, therefore X tends to 0 rather is equal to it. What is the right way to handle this? I have seen this throughout my model where some variables which are analytically equal to 0 in certain cases tend to 0 instead (e.g. 1e-324).

Reproducible callback failing:

function solve_model(p, u0)
    p = [2.767317774014066e-7, 0.7433220714328409, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0]
    u0 = [1.6382099999999998e13, 3.4061512877350002e10, 0.0, 100000.0, 1.461e7]

    # Callback for the max() function in the model - causes a discontinuity.
    # If max equation changes, this condition will have to change.
    condition(u, t, integrator) = p[2] * u[4] - u[2]
    affect!(integrator) = nothing
    max_cb = ContinuousCallback(condition, affect!, save_positions=(false,false))

    # Build the ODE Problem and Solve
    prob = ODEProblem(model, u0[1:end-1], (0.0, last(u0)), p)
    sol = solve(prob, Rosenbrock23(), callback=max_cb, tstops=stops, isoutofdomain=is_invalid_domain, maxiters=1e7)

    return sol
end


# Implements the differential equations that define the model
function model(du,u,p,t)
    # Paramaters
    mu_max = p[1]
    maintenance_per_cell = p[2]
    dOC_per_cell = p[3]

    carrying_capacity = p[4]
    pOC_input_rate = p[5]
    dOC_input_rate = p[6]
    inorganic_carbon_input_rate = p[7]
    inorganic_carbon_fixing_rate = p[8]
    inorganic_carbon_per_cell = p[9]
    eea_rate = p[10]
    Ks = p[11]

    # Load state conditions
    pOC_content = u[1]
    dOC_content = u[2]
    inorganic_carbon_content = u[3]
    cell_count = u[4]


    growth = mu_max * (dOC_content / (Ks + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

    required_dOC = maintenance_per_cell * cell_count

    dOC_missing = max(required_dOC - dOC_content, 0)
    starvation_deaths = dOC_missing == 0 ? 0 : dOC_missing / required_dOC_per_cell

    deaths =  starvation_deaths

    dOC_consumption = maintenance_per_cell * (cell_count - deaths)
    
    eea_rate = eea_rate * (1 / (1 + exp(-pOC_content))) * 2 - eea_rate
    
    du[1] = pOC_input_rate - eea_rate*cell_count

    du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count

    du[4] = growth - deaths
end

Example of working parameters:

p = [0.06, 0.1806273414373398, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0]