Optimization of trasportation problem in Pyomo

Question

I have an issue as the problem is not optimaze and I am trying to use a binary variable x to identify if the edge is used or not as the number of edges used should be minimized. The issue is that I cannot assign the correct value to x=1 when the Q of the edge is greater then zero. There is surely an issue or in the objective function or in the contraint 7. This is the code I am running:

import matplotlib.pyplot as plt
import networkx as nx
import pyomo.environ as pe
import pyomo.opt as po
import pandas as pd
import numpy as np

d = pd.read_excel("Grafo9Nodi_Bil.xlsx",header=0, index_col=0) #importo le posizioni dal file Excel
#impongo una lunghezza massima del singolo arco
LunMax=8
nodes = range(0,len((d["X"])))
distance = {}
pos={}
edges={}
for i in sorted(nodes):
    pos[i]=(d["X"][i],d["Y"][i])
    for j in sorted(nodes):
        if i!=j and np.round(pe.sqrt((d["X"][i]-d["X"][j])**2+(d["Y"][i]-d["Y"][j])**2),2)<=LunMax:
            distance [i,j]=np.round(pe.sqrt((d["X"][i]-d["X"][j])**2+(d["Y"][i]-d["Y"][j])**2),2)
               
edges={(i,j) for i in sorted(nodes) for j in sorted(nodes) if i!=j and np.round(pe.sqrt((d["X"][i]-d["X"][j])**2+(d["Y"][i]-d["Y"][j])**2),2)<=LunMax}       

model = pe.ConcreteModel()


path_edges = []
model.nodes = pe.Set(initialize=nodes)
model.edges = pe.Set(within=model.nodes*model.nodes, initialize=edges)

model.IN = pe.Var(model.nodes, domain=pe.NonNegativeReals,initialize=0) 
IN=model.IN
model.x = pe.Var(model.nodes*model.nodes, domain=pe.Binary,initialize=0)
x=model.x                 
model.Q = pe.Var(model.nodes*model.nodes, domain=pe.NonNegativeReals) 
Q=model.Q
model.OUT = pe.Var(model.nodes, domain=pe.NonNegativeReals,initialize=0) 
OUT=model.OUT

#funzione obiettivo
def obj_rule(model):
       TCC=sum(d["OpexC"][i]*IN[i] for i in model.nodes) 
       TTC=sum(Q[i,j]*1*distance[i,j] if Q[i,j] is not None and Q[i,j].value != 0 else 0 for [i,j] in model.edges)+sum(distance [i,j]*4*x[i,j] for [i,j] in model.edges if x[i,j].value != None and x[i,j].value != 0)
       TSC=sum(d["OpexS"][i]*OUT[i] for i in model.nodes)
       return TCC+TTC+TSC
model.Objf= pe.Objective (rule=obj_rule, sense = pe.minimize)

#vincoli
def constraint1(model,nodes):
    return IN[nodes]  <= d["IN"][nodes]
model.Const1=pe.Constraint(model.nodes, rule=constraint1)

def constraint2(model,nodes):
    return sum(IN[nodes] for nodes in model.nodes)  >= 1*sum(d["IN"][nodes] for nodes in model.nodes)
model.Const2=pe.Constraint(model.nodes, rule=constraint2)

def constraint4(model,nodes):
    return IN[nodes]+sum(Q[i,nodes] for i in model.nodes if [i,nodes] in model.edges) == sum(Q[nodes,i] for i in model.nodes if [i,nodes] in model.edges) +OUT[nodes]
model.Const4=pe.Constraint(model.nodes, rule=constraint4)

def constraint5(model,nodes):
    return OUT[nodes] <= d["OUT"][nodes]
model.Const5=pe.Constraint(model.nodes, rule=constraint5)

def constraint6(model,nodes):
    return sum(IN[nodes] for nodes in model.nodes)  == sum(OUT[nodes] for nodes in model.nodes)
model.Const6=pe.Constraint(model.nodes, rule=constraint6)

def constraint7(model, i, j):
    if Q[i,j].value != 0:
        return x[i,j] == 1
    else:
        return x[i,j] == 0
model.Const7 = pe.Constraint(model.edges, rule=constraint7)
    
Solver = po.SolverFactory ('gurobi', tee=True)
results = Solver.solve(model)

The result is not optimazed as x[i,j] are all equal to 1 while Q[i,j] have different values. Thanks in advance for the help

Answer 1

Here's a couple things to help you along. This sure looks like a homework project, so I'm reluctant to write the whole thing for you.

Some pyomo errors / tips:

1. You cannot use conditional (if) statements or .value in constructing the model. The values of the variables are unknown when the model is written. That's the solver's job. You will need to re-formulate some of the constraints and objective statements. This can be done for this type of model.

2. Your constraints have "name collision" in the variables. You are over-writing the variable that you take in from the function parameters. For instance in you constraint:

def constraint2(model,nodes):
    return sum(IN[nodes] for nodes in model.nodes)  >= 1*sum(d["IN"][nodes] for nodes in model.nodes)
model.Const2=pe.Constraint(model.nodes, rule=constraint2)

... what does nodes refer to? You are capturing it, but then you are using the same variable names in your summations. If you want to write a constraint for each node, you need to not over-write if you are going to use a sum.

3. Start off with very small data ... maybe just 4 or 5 nodes, and after you write each constraint, print the model and check the math to make sure it makes sense. (This of course assumes that you have the math correct, but it should look as you intend.) Use model.pprint()

Model structure

Your setup for the edges looks good, so why aren't you using that for the flow and usage variables to restrict them to the edges you made? I would expect:

model.q = pe.Var(model.edges, domain=pe.NonNegativeReals, doc='flow')

In order to deal with the balance of flow constraints, you are going to need to organize your edges such that you can find the ones that flow in or out of each node. There are several ways to do that, but perhaps the simplest is just make some dictionaries... Here are a couple ideas:

from collections import defaultdict
import pyomo.environ as pe

nodes = [1, 2, 3, 4, 5]
edges = [(1, 3), (1, 4), (5, 2)]

edges_out_of_node_dict = defaultdict(list)
for edge in edges:
    edges_out_of_node_dict[edge[0]].append(edge)

print(edges_out_of_node_dict)

m   = pe.ConcreteModel('flow')
m.N = pe.Set(initialize=nodes)
m.E = pe.Set(within=m.N*m.N, initialize=edges)

m.q = pe.Var(m.E, domain=pe.NonNegativeReals, doc='flow')

# example constraint:  the total flow out of each connected node is less than 5
def outflow_limit(m, node):
    return sum(m.q[edge] for edge in edges_out_of_node_dict[node]) <= 5
# we only want to apply this to nodes that have outbound arcs or we will get a key error, so...
nodes_with_connections = edges_out_of_node_dict.keys()
m.constraint1 = pe.Constraint(nodes_with_connections, rule=outflow_limit, doc='outflow limit')


m.pprint()

Output:

defaultdict(<class 'list'>, {1: [(1, 3), (1, 4)], 5: [(5, 2)]})
4 Set Declarations
    E : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain   : Size : Members
        None :     2 : E_domain :    3 : {(1, 3), (1, 4), (5, 2)}
    E_domain : Size=1, Index=None, Ordered=True
        Key  : Dimen : Domain : Size : Members
        None :     2 :    N*N :   25 : {(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5)}
    N : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    5 : {1, 2, 3, 4, 5}
    constraint1_index : Size=1, Index=None, Ordered=False
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 5}

1 Var Declarations
    q : flow
        Size=3, Index=E
        Key    : Lower : Value : Upper : Fixed : Stale : Domain
        (1, 3) :     0 :  None :  None : False :  True : NonNegativeReals
        (1, 4) :     0 :  None :  None : False :  True : NonNegativeReals
        (5, 2) :     0 :  None :  None : False :  True : NonNegativeReals

1 Constraint Declarations
    constraint1 : outflow limit
        Size=2, Index=constraint1_index, Active=True
        Key : Lower : Body            : Upper : Active
          1 :  -Inf : q[1,3] + q[1,4] :   5.0 :   True
          5 :  -Inf :          q[5,2] :   5.0 :   True

6 Declarations: N E_domain E q constraint1_index constraint1