from IPython.display import Image
fig = Image(url='http://mathopt.sakura.ne.jp/LP.png')
fig 

import pulp

prob = pulp.LpProblem("myProblem",pulp.LpMinimize)

z1 = pulp.LpVariable("z1",0,3)

z2 = pulp.LpVariable("z2",0,1)

prob += z1 + z2 <= 2

print prob

prob += -1*z1 - 2*z2 


print prob

prob+=20*z1+30*z2

print prob

prob += -z1-2*z2

status = prob.solve()

print pulp.LpStatus[status]

print pulp.value(prob.objective)

print pulp.value(z1), pulp.value(z2)

prob_d = pulp.LpProblem("myProblem",pulp.LpMinimize)
x1=pulp.LpVariable("x1",0,3)
x2=pulp.LpVariable("x2",0,1)
prob_d+=-x1-2*x2
prob_d+=x1+x2<=2,"First_Constraint"
print prob_d

status = prob_d.solve()
print pulp.LpStatus[status]

print prob_d.constraints.keys()

print prob_d.constraints["First_Constraint"].pi

fig = Image(url='http://mathopt.sakura.ne.jp/VRP.png')
fig 

routing_model=pulp.LpProblem("VRP",pulp.LpMinimize)

pow_s=['s1','s2','s3']
pow_s_vars=pulp.LpVariable.dicts('x', pow_s, lowBound=0, upBound=1, cat=pulp.LpBinary)
print pow_s_vars

print pow_s_vars['s1']

cargos=["c0","c1","c2","c3"];  print cargos


v_feas_rt={}; v_feas_rt['v0']=[]
v_feas_rt['v0']+=[('v0',('c0'))]
v_feas_rt['v0']+=[('v0',('c0','c1'))]
v_feas_rt['v0']+=[('v0',('c1','c3'))]
print v_feas_rt['v0']

route_cost={}

route_cost[('v0',('c0'))]=100
route_cost[('v0',('c0','c1'))]=200
route_cost[('v0',('c1','c3'))]=300


v_feas_rt['v1']=[]
v_feas_rt['v1']+=[('v1',('c3'))]
v_feas_rt['v1']+=[('v1',('c0','c2'))]
v_feas_rt['v1']+=[('v1',('c0','c3'))]
print v_feas_rt['v1']

route_cost[('v1',('c3'))]=150
route_cost[('v1',('c0','c2'))]=220
route_cost[('v1',('c0','c3'))]=180
print route_cost

vehicles,feas_rt=['v0','v1'],[]
for v in vehicles:
        feas_rt+=v_feas_rt[v]
print feas_rt

x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
print x

y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
print y

[x[route] for route in feas_rt if 'c0' in route[1] ]

for c in cargos:
    routing_model+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c

print routing_model

for v in vehicles:
    routing_model += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
print routing_model

routing_model+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print routing_model

dual_cargo,dual_vehicle={},{}

status = routing_model.solve()
print pulp.LpStatus[status]

for c in cargos:
    dual_cargo[c]=routing_model.constraints["Must_assign_%s"%c].pi

for v in vehicles:
    dual_vehicle[v]=routing_model.constraints["Vehicle_assign_%s"%v].pi

print dual_cargo,dual_vehicle

import networkx as nx
G={}
G['v0']=nx.DiGraph()
G['v0'].add_weighted_edges_from([('s','c0',100),('s','c1',50),('c0','c1',20),('c0','c2',40),('c0','c3',30),('c1','c3',35),('c2','t',70),('c3','c2',15),('c3','t',45)])
G['v1']=nx.DiGraph()
G['v1'].add_weighted_edges_from([('s','c0',150),('s','c1',10),('c0','c2',40),('c0','c3',30),('c1','c3',35),('c2','t',70),('c3','t',45)])



import networkx as nx
pos=nx.shell_layout(G['v0']);
nx.draw_shell(G['v0'],node_size=1000,node_color='w');
nx.draw_networkx_edge_labels(G['v0'],pos);

import networkx as nx
pos=nx.shell_layout(G['v1']);
nx.draw_shell(G['v1'],node_size=1000,node_color='w');
nx.draw_networkx_edge_labels(G['v1'],pos);

import pulp
cargos=["c0","c1","c2","c3"]
vehicles,feas_rt=['v0','v1'],[]
v_feas_rt={}
v_feas_rt['v0'],v_feas_rt['v1']=[],[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]

master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
#for v in vehicles:
#    master_p += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p


status=master_p.solve()
print pulp.LpStatus[status]

dual_cargo={}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
print dual_cargo


print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
	

for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
	

route_cost={}
v_feas_rt['v0']+=[('v0',('c0','c1','c2','c3'))]
v_feas_rt['v1']+=[('v1',('c1','c3'))]
route_cost[('v0',('c0','c1','c2','c3'))]=240
route_cost[('v1',('c1','c3'))]=90

import pulp
cargos=["c0","c1","c2","c3"]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p += sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p


status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle

print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	

v_feas_rt['v1']+=[('v1',('c0','c2'))]
route_cost[('v1',('c0','c2'))]=260

import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
print feas_rt
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p

status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle

print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	

v_feas_rt['v0']+=[('v0',('c1','c2','c3'))]
route_cost[('v0',('c1','c2','c3'))]=170

import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p

status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle

print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	

v_feas_rt['v0']+=[('v0',('c0','c2','c3'))]
v_feas_rt['v0']+=[('v0',('c0','c2'))]
route_cost[('v0',('c0','c2','c3'))]=215
route_cost[('v0',('c0','c2'))]=210

import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
master_p=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpContinuous)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpContinuous)
for c in cargos:
    master_p+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    master_p+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
master_p+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
print master_p

status=master_p.solve()
print pulp.LpStatus[status]
dual_cargo,dual_vehicle={},{}
for c in cargos:
	dual_cargo[c]=master_p.constraints["Must_assign_%s"%c].pi
for v in vehicles:
    dual_vehicle[v]=master_p.constraints["Vehicle_assign_%s"%v].pi
print dual_cargo,dual_vehicle

print "s-t paths for v0"
for path in nx.all_simple_paths(G['v0'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v0'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	
print "s-t paths for v1"
for path in nx.all_simple_paths(G['v1'],source='s',target='t'):
	tcost=0
	for i in range(len(path)-1):
		tcost+=G['v1'][path[i]][path[i+1]]['weight']
		if path[i]!='s':
			tcost-=dual_cargo[path[i]]
	print path,tcost	

	

import pulp
feas_rt=[]
for v in vehicles:
	feas_rt+=v_feas_rt[v]
spp=pulp.LpProblem("MasterProblem",pulp.LpMinimize)
x=pulp.LpVariable.dicts('x',feas_rt,lowBound=0,upBound=1,cat=pulp.LpBinary)
y=pulp.LpVariable.dicts('y',cargos,lowBound=0,upBound=1,cat=pulp.LpBinary)
for c in cargos:
    spp+=sum([x[route] for route in feas_rt if c in route[1]])+y[c]==1,"Must_assign_%s"%c
for v in vehicles:
    spp+= sum( [x[route] for route in v_feas_rt[v]])<=1,"Vehicle_assign_%s"%v
spp+=sum([route_cost[route]*x[route] for route in feas_rt])+sum([1000*y[c] for c in cargos])
status=spp.solve()
print pulp.LpStatus[status]

for v in spp.variables():
    if v.varValue>0.1:
        print v.name,v.varValue