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]