Just throwing another tip out there: if you're finding yourself having to model anything mildly nonlinear, but find yourself shackled by the linearity constraints of a Mixed Integer Program (MIP), you can use a statistical/ML technique called MARS [1] to regress a piecewise-linear model from your nonlinear data to approximate said nonlinear function. MARS models are essentially models with constant or max functions as basis functions, e.g.
y = a₀ + a₁ max(0, x − b₁) + a₂ max(0, x − b₂) + ...
max(0,x) functions are very easy to model in MIPs.
y = a₀ + a₁ max(0, x − b₁) + a₂ max(0, x − b₂) + ...
max(0,x) functions are very easy to model in MIPs.
[1] https://en.wikipedia.org/wiki/Multivariate_adaptive_regressi...