pyomo/ipopt: nonlinear network optimization not converging

Question

The question (very short version)

Why can I not decrease the lower boundary for the decision variable model.v_dot (see implementation) below 30 ? As soon as I do so, the solver doesn't return results. After only one iteration I get the solver error: EXIT: Restoration Failed!

The goal is to have the solver return results for a v_dot lower boundary equal to zero.

Long version

I am trying to optimize the diameters of a gas pipeline network.

The topology is given (the node location and hence the pipeline lengths), the edges, the sources' and sinks' location and their throughput.

Network Topology

The problem implementation inputs are:

1. the list of nodes
2. the list of edges
3. the dictionary of lengths (associating a length to each edge)
4. the dictionary of sources (associating a volume flow rate to each node)
5. the dictionary of sinks (same thing) (please see the data below if you want to implement the same instance as mine)

And the outputs should be:

1. the diameters of each pipeline
2. the pressure at each node
3. the (normal) volume flow rate of each pipeline

Important to mention is that:

• there are 2 equality constraints, one physically describing the flow in a pipeline (the pressure drop equation) and one insuring the mass flow conservation in the network
• the pressure is bounded by a maximum and a minimum

Implementation

The implementation on pyomo is as follows:

# declare the model
model = pyo.ConcreteModel()

# declare the sets and parameters

model.junctions = pyo.Set(initialize=junctions_pyo)
model.edges = pyo.Set(initialize=edges_pyo, within= model.junctions*model.junctions)

p_max = 7000  # kpa
p_min = 1600
a_0 = 1  # coefficients for object fct
a_1 = 1
a_2 = 1
a = 0.3  # coefficient for pressure drop equality constraint

model.lengths = pyo.Param(model.edges, initialize=lengths_dict_pyo)

# declare the variables
model.d = pyo.Var(model.edges, bounds=(0.01, 100), within=pyo.NonNegativeReals) 
model.pi = pyo.Var(model.junctions, bounds=(p_min ** 2, p_max ** 2), within=pyo.NonNegativeReals)
model.v_dot = pyo.Var(model.edges, bounds=(30, 1500), within=pyo.Reals)

# define the objective function
def obj_rule(model):
    return sum((a_0 + a_1 * model.d[i, j] + a_2 * model.d[i, j] ** 2) * model.lengths[i, j] for (i, j) in model.edges)
model.obj = pyo.Objective(rule=obj_rule, sense= pyo.minimize)

# define the constraints
# pressure drop equality constraint
def flow_rule(model, i, j):
    return ((model.pi[i] - model.pi[j]) * model.d[i, j] ** 5) == (a * model.v_dot[i, j]**2 * model.lengths[i, j])
model.flow = pyo.Constraint(model.edges, rule=flow_rule)

def conservation_rule(model, i):
    inflow = sources_dict[i] + sum(model.v_dot[k, i] for k in model.junctions if (k, i) in model.edges)
    outflow = sinks_dict[i] + sum(model.v_dot[i, k] for k in model.junctions if (i, k) in model.edges)
    return inflow == outflow
model.conservation = pyo.Constraint(model.junctions, rule=conservation_rule)


# solving
solver = SolverFactory("ipopt")
results = solver.solve(model, tee=True)
model.display()

Results

When attaching a a single source to node 0 and a single sink to node 10 (see the data below) the results should show the fluid flowing through the path 0-1-2-6-4-8-10. However, a lower boundary of 30 has to be assigned to vdot, otherwise the solver doesn't return results(this is discussed further in the section Results discussion). The results I am getting are:

• from the solver, the returned result is a point of local infeasibility

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
WARNING: Loading a SolverResults object with a warning status into
    model.name="unknown";
      - termination condition: infeasible
      - message from solver: Ipopt 3.11.1\x3a Converged to a locally
        infeasible point. Problem may be infeasible.

• objective function value

  Objectives:
obj : Size=1, Index=None, Active=True

     |Key  | Active | Value|
     |-----|--------|------|
     |None |   True | 1091.8088253777207|

• diameters

• (squared) pressures

• (normal) flow rate

• pressure drop constraint

• mass conservation constraint

Results discussion

The diameters for the pipelines/edges (path 0-1-2-4-6-8-10) on which the fluid is supposed to flow are good (consistent to the results of another python package, namely pandapipes) and the pressure values are also what they are supposed to be. However the diameters of the other edges should be reduced to zero by the solver in order to minimize the objective function. This is not the case and the reason might be the lower boundary on the volume flow rate (vdot). however, as soon as I try to lower this bound slightly, the solver returns very bad results. For lower (<= 20) boundary values, the solver doesn't return results at all. After only one iteration I get the solver error: EXIT: Restoration Failed!

The goal is to have the solver return results for v_dot lower boundary equal to 0 (in order to have no fluid flowing on the pipelines where its not supposed to flow)

What I've tried to do:

• I tried to initialize the decision variables at values which are near the solution point.
• I tried to scale down the pressures in the pressure drop equality constraint by multiplying the right hand side of the equation by 10^(-4).

Mathematical formulation

The mathematical problem formulation is as follows:

Problem Formulation

Formulation Table

The data

1. junctions_pyo : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

2. edges_pyo : [(0, 1), (1, 2), (2, 3), (3, 4), (3, 5), (2, 6), (6, 8), (8, 9), (8, 10),(6, 7)]

3. lengths_dict_pyo: {(0, 1): 24.66, (1, 2): 45.58, (2, 3): 94.43, (3, 4): 21.97, (3, 5): 24.69, (2, 6): 77.85, (6, 8): 61.81, (8, 9): 12.30, (8, 10): 27.58, (6, 7): 27.60}

4. sources_dict: {0: 1219.33, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0}

5. sinks_dict : {0: 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 1219.33}

Answer 1

That IPOPT message means that IPOPT could not find a feasible solution to your problem. The reason could be either that:

• Once you set that value below 30, IPOPT can no longer find that basin of attraction (or that basin vanishes). Another feasible solution might exist (unless your problem is convex), but IPOPT can't find it.
• Your problem actually becomes infeasible.

Things you can do:

• Start removing constraints to see if your problem will become feasible at some point.
• Use a deterministic global optimisation solver (e.g. Octeract Engine, BARON, Couenne), which can actually prove whether your problem is infeasible or not. If it's not, they are much more likely to find a solution than standalone IPOPT is.
• Provide a feasible solution in a basin of attraction you want IPOPT to explore.