I have been trying my hand at a MINLP process synthesis problem that has been running into numerical difficulties due to small non-zero values.

In essence, the model contains continuous flow variables which are supposed to be 0 if their corresponding binary variable is turned off, such as in this example:

m = ConcreteModel()
m.flow = Var(domain = NonNegativeReals)
m.y = Var(domain = Binary)


def FlowExistence(model):
    return model.flow -model.y*F_max <=0

The issue is that because these constraints have a certain tolerance in say IPOPT, the solver returns values for the flow that are not quite zero, but on the order of 1e-9. This causes issues in other constraint in a similar form to the one below, which are needed for mixers and splitters.

def ComponentBalance(model,component):
    return model.flow*model.concentration[component] -model.ComponentFlow[component] <=0

In the above constraint, the low non-zero value for flow creates an extreme entry in the jacobian, causing numerical issues.

Thanks in advance!

I have tried to initialize such variables at 0, yet they do not remain there throughout the solver iteration. Since I have seen many implementations of such constraints without any mentions of difficulties, I was wondering how one could avoid such issues without having to fix the continuous variable between iterations.