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.
The problem implementation inputs are:
- the list of nodes
- the list of edges
- the dictionary of lengths (associating a length to each edge)
- the dictionary of sources (associating a volume flow rate to each node)
- 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:
- the diameters of each pipeline
- the pressure at each node
- 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|
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:
The data
junctions_pyo : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
edges_pyo : [(0, 1), (1, 2), (2, 3), (3, 4), (3, 5), (2, 6), (6, 8), (8, 9), (8, 10),(6, 7)]
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}
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}
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}