I've noticed clamp is used four times in solve.jl. The problem with clamp is that it is non-smooth so it won't play nice with anything that utilizes derivatives, like Newton Raphson.

A straightforward clamp_smooth is:

function clamp_smooth(x, lo, hi, m)
    d = (hi - lo)
    @assert m <= (0.5 * d)
    x_norm = (x - lo) / d
    a = 1 / (1 - m)
    
    if x_norm < 0
        y = 0
    elseif x_norm < m
        y = (a / (2 * m)) * x_norm^2
    elseif x_norm < (1 - m)
        y = a * x_norm + 0.5 * (1 - a)
    elseif x_norm < 1
        y = 1 - (a / (2 * m)) * (1 - x_norm)^2
    else
        y = 1
    end
    
    return lo + y * d
end

Where m is the smoothing parameter. This is only C1 continuous which is fine for us now but maybe needs some attention if we need Hessians in the future.

In terms of its behavior, this is with a large m (clamp in orange):

Applied to a pump, it means that the pump is a bit "lazy", then "panicky": it pumps slightly less when at half-capacity, and it pumps a little more when over half-capacity. In general, the value for m is chosen very small, so deviations are negligible.

The reduction_factor is similar in idea:

It uses a cubic approximation to go from 0 to 1 from 0 to threshold (but it doesn't approximate a piecewise linear function like clamp well).

Although it's more costly to evaluate than a basic clamp, it should be helpful for the solver and give some performance benefits.