# -*- coding: utf-8 -*-
"""
Bicriteria Transportation Problem
Y. P. Aneja, K. P. K. Nair
Management Science, Vol. 25, No. 1 (Jan., 1979), pp. 73-78

Solved with epsilon contraint for demonstration purpose. Here a correct 
solution is known so the potential appoximation of the epsilon constraint
method can be demonstated.
"""
from gurobipy import Model, GRB, quicksum
import matplotlib.pyplot as plt


supply = [8, 19, 17]

demand = [11,3,14,16]

cost = [[1,2,7,7],
        [1,9,3,4],
        [8,9,4,6]]

distance = [[4,4,3,4],
            [5,8,9,10],
            [6,2,5,1]]

M = range(3)
N = range(4)

def createModel(is_f1):
    model= Model("Bi-Objective Transport Problem -%r" % is_f1)
    
    # Decide from which supplier to deliver to which customer
    x = model.addVars(M, N, lb=0, vtype=GRB.CONTINUOUS, name="x")
    
    model.update()
    # All supply has to be used
    model.addConstrs(quicksum(x[i,j] for j in N) == supply[i] for i in M)
    # All demand has to be satisfied
    model.addConstrs(quicksum(x[i,j] for i in M) == demand[j] for j in N)
    
    if is_f1:
        # f1: Optimise for cost
        model.setObjective(quicksum(cost[i][j] * x[i, j] for i in M for j in N), GRB.MINIMIZE)
    else:
        # f1: Optimise for distance
        model.setObjective(quicksum(distance[i][j] * x[i, j] for i in M for j in N), GRB.MINIMIZE)
        
    return model, x

r1, x1 = createModel(True)
r1.optimize()
r1.getVars()
r2, x2 = createModel(False)
r2.optimize()
r2.getVars()

other_target = sum(distance[i][j] * x1[i, j].x for i in M for j in N) # f2
opt_f1 = (r1.objVal, other_target)
print("Optimal cost: %s" % r1.objVal)
print("With distance of: %s" % other_target)   


other_target = sum(cost[i][j] * x2[i, j].x for i in M for j in N) # f1
opt_f2 = (other_target, r2.objVal)
print("Optimal distance: %s" % r2.objVal)
print("With cost of: %s" % other_target)  


fig, ax = plt.subplots()
plt.xlabel('cost')
plt.ylabel('distance')
plt.scatter([opt_f1[0], opt_f2[0]], [opt_f1[1], opt_f2[1]])
offset = 1; 
plt.text(opt_f1[0] + offset, opt_f1[1] + offset, "f1 opt")
plt.text(opt_f2[0] + offset, opt_f2[1] + offset, "f2 opt")

  
"""
Epsilon test
""" 
r = opt_f1[1] - opt_f2[1] # Takes value for f2 from both
x_values_e = []
y_values_e = []
maxIter = 5
for i in range(maxIter + 1):
    epsilon = opt_f2[1] + (r / maxIter) * i 
    model, x = createModel(True)
    
    # Epsilon constraint
    model.addConstr(quicksum(distance[i][j] * x[i, j] for i in M for j in N) <= epsilon)
    model.optimize()
    model.getVars()
    
    # Save values
    x_value = model.objVal
    y_value = sum(distance[i][j] * x[i, j].x for i in M for j in N) # f2
    x_values_e.append(x_value)
    y_values_e.append(y_value)
    
    
    # Write the epsilon values next to the points
    if i != 0 and i != maxIter:
        plt.text(x_value + offset, y_value + offset, "e: %.1f" % epsilon)

# Plot the points that were used to create the front
plt.scatter(x_values_e, y_values_e)

# Plot the front
x_values_e.reverse()
y_values_e.reverse()    
plt.plot(x_values_e, y_values_e) 


"""
Complete pareto front for show purpose 
Plot to see how epsilon constraint error develops
"""
x_values = [143,156,176,186,208]
y_values = [265,200,175,171,167]    
plt.plot(x_values, y_values)