Python's PuLP module can't solve my integer linear programming model. However, I'm convinced there are integer solutions to my problem. I definitely don't understand why PuLP isn't returning any solutions. Does anyone have any ideas? Here's my code:
from pulp import *
import networkx as nx
def solver(k,G,safe,Qmax):
n = len(G)
mu = [1.0/n]*n #Temporary
unsafe = [i for i in range(n) if not i in safe]
maet = LpProblem("ResolveMAET", LpMinimize)
q = {i:(str(i)) for i in G}
p = {(i,j):(str(i)+","+str(j)) for i in unsafe for j in G.neighbors(i)}
q_var = LpVariable.dicts("q",[qi for qi in q.values()],lowBound = 0, upBound = Qmax, cat = 'Continuous')
p_var = LpVariable.dicts("p",[pij for pij in p.values()],lowBound = 0, upBound= 1, cat = 'Binary')
### Objective function
maet += lpSum([ mu[i] * q_var[q[i]] for i in q])
### Constraints
#Markov Chain (Constraint n°1)
for i in unsafe:
Mij = 1.0 / G.degree(i) #Temporary
maet += q_var[q[i]] - lpSum([Mij*q_var[q[j]] for j in G.neighbors(i)]) + Qmax * lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) >= lpSum([ Mij*G[i][j]["weight"] for j in G.neighbors(i) ])
#Sign (Constraint n°2)
for i,j in p:
maet += q_var[q[i]] - q_var[q[j]] - Qmax*p_var[p[i,j]] >= G[i][j]["weight"] - Qmax
#At Most One Sign Per Node (Constraint n°3)
for i in unsafe:
maet += lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) <= 1
#At Most k Signs (Constraint n°4)
maet += lpSum([p_var[p[i,j]] for i,j in p]) <= k
#Shelter (Constraint n°5)
for i in safe:
maet += q_var[q[i]] <= 0 #=0 car q_var[q[i]] >= 0
### Résolution
maet.solve()
return [(i,j) for i,j in p if p_var[p[i,j]].varValue==1]
######################
v= 4.0 #km/h
# Graph initialization
G = nx.Graph()
file_density = open("densite.txt", 'r', encoding="UTF8")
mu = dict()
id_to_pos = {}
pos_to_id = {}
id = 0
for l in file_density:
data = l.split(",")
for i in range(0,len(data),3):
vertex = (float(data[i][1:]),float(data[i+1]))
number = int(data[i+2][6:])
id_to_pos[id] = vertex
pos_to_id[vertex] = id
mu[id] = number
id+=1
file_density.close()
# Ajout des noeuds
G.add_nodes_from([i for i in range(id)])
file_graph = open("graph.txt", 'r', encoding="UTF8")
file_graph.readline()
for l in file_graph:
data = l.split(" ")
weight = float(data[1])
extremity1 = pos_to_id[(float(data[2].split(",")[0][1:]),float(data[2].split(",")[1][0:]))]
extremity2 = pos_to_id[(float(data[3].split(",")[0][1:]),float(data[3].split(",")[1][0:]))]
G.add_edge(extremity1, extremity2, weight=v*weight)
file_graph.close()
from random import randint
safe = [randint(0,len(G)) for i in range(10)]
print(solver(100,G,safe,60000.0))
I've checked the variables and coefficients, and they all seem correct. I think the problem comes from the Qmax constant.
Remark : This constant transforms conditional constraints into linear constraints. For example, the constraint $q_i = t_{i,j} + q_j$ if $p_{i,j} = 1$ is replaced by the linear constraint $q_i \geq t_{i,j} + q_j - Q_{max} (1-p_{i,j})$ (contraint n°2)
I just made up a small set of fake data and ran this and it returns optimal with some kind of solution. I have no idea what it is doing or what the meaning of the result is, but the model appears to be generating integer results. Perhaps this will help you in troubleshooting.
Also note that I dropped in a print statement to print out the model, which you should be doing for troubleshooting.
from pulp import *
import networkx as nx
def solver(k,G,safe,Qmax):
n = len(G)
mu = [1.0/n]*n #Temporary
unsafe = [i for i in range(n) if not i in safe]
maet = LpProblem("ResolveMAET", LpMinimize)
q = {i:(str(i)) for i in G}
p = {(i,j):(str(i)+","+str(j)) for i in unsafe for j in G.neighbors(i)}
q_var = LpVariable.dicts("q",[qi for qi in q.values()],lowBound = 0, upBound = Qmax, cat = 'Continuous')
p_var = LpVariable.dicts("p",[pij for pij in p.values()],lowBound = 0, upBound= 1, cat = 'Binary')
### Objective function
maet += lpSum([ mu[i] * q_var[q[i]] for i in q])
### Constraints
#Markov Chain (Constraint n°1)
for i in unsafe:
Mij = 1.0 / G.degree(i) #Temporary
maet += q_var[q[i]] - lpSum([Mij*q_var[q[j]] for j in G.neighbors(i)]) + Qmax * lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) >= lpSum([ Mij*G[i][j]["weight"] for j in G.neighbors(i) ])
#Sign (Constraint n°2)
for i,j in p:
maet += q_var[q[i]] - q_var[q[j]] - Qmax*p_var[p[i,j]] >= G[i][j]["weight"] - Qmax
#At Most One Sign Per Node (Constraint n°3)
for i in unsafe:
maet += lpSum([p_var[p[i,j]] for j in G.neighbors(i)]) <= 1
#At Most k Signs (Constraint n°4)
maet += lpSum([p_var[p[i,j]] for i,j in p]) <= k
#Shelter (Constraint n°5)
for i in safe:
maet += q_var[q[i]] <= 0 #=0 car q_var[q[i]] >= 0
### Résolution
print(maet) # <-- take a look at the math and make sure it matches your math model before moving to larger dataset
maet.solve()
return [(i,j) for i,j in p if p_var[p[i,j]].varValue==1]
######################
v= 4.0 #km/h
# file_density = open("densite.txt", 'r', encoding="UTF8")
# mu = dict()
# id_to_pos = {}
# pos_to_id = {}
# id = 0
# for l in file_density:
# data = l.split(",")
# for i in range(0,len(data),3):
# vertex = (float(data[i][1:]),float(data[i+1]))
# number = int(data[i+2][6:])
# id_to_pos[id] = vertex
# pos_to_id[vertex] = id
# mu[id] = number
# id+=1
# file_density.close()
# # Ajout des noeuds
# G.add_nodes_from([i for i in range(id)])
# file_graph = open("graph.txt", 'r', encoding="UTF8")
# file_graph.readline()
# for l in file_graph:
# data = l.split(" ")
# weight = float(data[1])
# extremity1 = pos_to_id[(float(data[2].split(",")[0][1:]),float(data[2].split(",")[1][0:]))]
# extremity2 = pos_to_id[(float(data[3].split(",")[0][1:]),float(data[3].split(",")[1][0:]))]
# G.add_edge(extremity1, extremity2, weight=v*weight)
# file_graph.close()
connections = [
(0, 1, 0.5),
(1, 2, 2.0),
(1, 3, 2.5),
(2, 4, 3.0),
(3, 5, 6.4),
(4, 5, 1.0),]
# Graph initialization
G = nx.Graph()
for c in connections:
G.add_edge(c[0], c[1], weight=c[2])
print(list(G.neighbors(1)))
from random import randint
# safe = [randint(0,len(G)) for i in range(10)]
safe = [2, ]
print(solver(3,G,safe,60.0))