multi-modal transportation.py

# -*- coding: utf-8 -*-
"""
Created on Thu Mar 15 09:32:04 2018

@author: Ken Huang
"""
from docplex.mp.model import Model
from itertools import product
import numpy as np
import cvxpy as cp
import pandas as pd
import json


class MMT:
    '''a Model class that solves the multi-model transportation optimization problem.'''

    def __init__(self, framework='DOCPLEX'):
        # parameters
        self.portSpace = None
        self.dateSpace = None
        self.goods = None
        self.indexPort = None
        self.portIndex = None
        self.maxDate = None
        self.minDate = None
        self.tranCost = None
        self.tranFixedCost = None
        self.tranTime = None
        self.ctnVol = None
        self.whCost = None
        self.kVol = None
        self.kValue = None
        self.kDDL = None
        self.kStartPort = None
        self.kEndPort = None
        self.kStartTime = None
        self.taxPct = None
        self.transitDuty = None
        self.route_num = None
        self.available_routes = None
        # decision variables
        self.var = None
        self.x = None
        self.var_2 = None
        self.y = None
        self.var_3 = None
        self.z = None
        # result & solution
        self.xs = None
        self.ys = None
        self.zs = None
        self.whCostFinal = None
        self.transportCost = None
        self.taxCost = None
        self.solution_ = None
        self.arrTime_ = None
        self.objective_value = None
        # helping variables
        self.var_location = None
        self.var_2_location = None
        self.var_3_location = None

        if framework not in ['CVXPY', 'DOCPLEX']:
            raise ValueError('Framework not supported, the model only supports CVXPY and DOCPLEX')
        else:
            self.framework = framework

    def set_param(self, route, order):
        '''set model parameters based on the read-in route and order information.'''

        bigM = 100000
        route = route[route['Feasibility'] == 1]
        route['Warehouse Cost'][route['Warehouse Cost'].isnull()] = bigM
        route = route.reset_index()

        portSet = set(route['Source']) | set(route['Destination'])

        self.portSpace = len(portSet)
        self.portIndex = dict(zip(range(len(portSet)), portSet))
        self.indexPort = dict(zip(self.portIndex.values(), self.portIndex.keys()))

        self.maxDate = np.max(order['Required Delivery Date'])
        self.minDate = np.min(order['Order Date'])
        self.dateSpace = (self.maxDate - self.minDate).days
        startWeekday = self.minDate.weekday() + 1
        weekday = np.mod((np.arange(self.dateSpace) + startWeekday), 7)
        weekday[weekday == 0] = 7
        weekdayDateList = {i: [] for i in range(1, 8)}
        for i in range(len(weekday)):
            weekdayDateList[weekday[i]].append(i)
        for i in weekdayDateList:
            weekdayDateList[i] = json.dumps(weekdayDateList[i])

        source = list(route['Source'].replace(self.indexPort))
        destination = list(route['Destination'].replace(self.indexPort))
        DateList = list(route['Weekday'].replace(weekdayDateList).apply(json.loads))

        self.goods = order.shape[0]
        self.tranCost = np.ones([self.portSpace, self.portSpace, self.dateSpace]) * bigM
        self.tranFixedCost = np.ones([self.portSpace, self.portSpace, self.dateSpace]) * bigM
        self.tranTime = np.ones([self.portSpace, self.portSpace, self.dateSpace]) * bigM

        for i in range(route.shape[0]):
            self.tranCost[source[i], destination[i], DateList[i]] = route['Cost'][i]
            self.tranFixedCost[source[i], destination[i], DateList[i]] = route['Fixed Freight Cost'][i]
            self.tranTime[source[i], destination[i], DateList[i]] = route['Time'][i]

        self.transitDuty = np.ones([self.portSpace, self.portSpace]) * bigM
        self.transitDuty[source, destination] = route['Transit Duty']

        # make the container size of infeasible routes to be small enough, similar to bigM
        self.ctnVol = np.ones([self.portSpace, self.portSpace]) * 0.1
        self.ctnVol[source, destination] = route['Container Size']
        self.ctnVol = self.ctnVol.reshape(self.portSpace, self.portSpace, 1)
        self.whCost = route[['Source', 'Warehouse Cost']].drop_duplicates()
        self.whCost['index'] = self.whCost['Source'].replace(self.indexPort)
        self.whCost = np.array(self.whCost.sort_values(by='index')['Warehouse Cost'])
        self.kVol = np.array(order['Volume'])
        self.kValue = np.array(order['Order Value'])
        self.kDDL = np.array((order['Required Delivery Date'] - self.minDate).dt.days)
        self.kStartPort = np.array(order['Ship From'].replace(self.indexPort))
        self.kEndPort = np.array(order['Ship To'].replace(self.indexPort))
        self.kStartTime = np.array((order['Order Date'] - self.minDate).dt.days)
        self.taxPct = np.array(order['Tax Percentage'])

        # add available route indexes
        self.route_num = route[['Source', 'Destination']].drop_duplicates().shape[0]
        routes = route[['Source', 'Destination']].drop_duplicates().replace(self.indexPort)
        self.available_routes = list(zip(routes['Source'], routes['Destination']))
        # localization variables of decision variables in the matrix
        var_location = product(self.available_routes, range(self.dateSpace), range(self.goods))
        var_location = [(i[0][0], i[0][1], i[1], i[2]) for i in var_location]
        self.var_location = tuple(zip(*var_location))

        var_2_location = product(self.available_routes, range(self.dateSpace))
        var_2_location = [(i[0][0], i[0][1], i[1]) for i in var_2_location]
        self.var_2_location = tuple(zip(*var_2_location))

        self.var_3_location = self.var_2_location

    def build_model(self):
        '''overall function to build up model objective and constraints'''
        if self.framework == 'CVXPY':
            self.cvxpy_build_model()
        elif self.framework == 'DOCPLEX':
            self.cplex_build_model()

    def cvxpy_build_model(self):
        '''build up the mathematical programming model's objective and constraints using CVXPY framework.'''

        # 4 dimensional binary decision variable matrix
        self.var = cp.Variable(self.route_num * self.dateSpace * self.goods, boolean=True, name='x')
        self.x = np.zeros((self.portSpace, self.portSpace, self.dateSpace, self.goods)).astype('object')
        self.x[self.var_location] = list(self.var)
        # 3 dimensional container number matrix
        self.var_2 = cp.Variable(self.route_num * self.dateSpace, integer=True, name='y')
        self.y = np.zeros((self.portSpace, self.portSpace, self.dateSpace)).astype('object')
        self.y[self.var_2_location] = list(self.var_2)
        # 3 dimensional route usage matrix
        self.var_3 = cp.Variable(self.route_num * self.dateSpace, boolean=True, name='z')
        self.z = np.zeros((self.portSpace, self.portSpace, self.dateSpace)).astype('object')
        self.z[self.var_3_location] = list(self.var_3)
        # warehouse related cost
        warehouseCost, arrTime, stayTime = self.warehouse_fee(self.x)
        ###objective###
        transportCost = np.sum(self.y * self.tranCost) + np.sum(self.z * self.tranFixedCost)
        transitDutyCost = np.sum(np.sum(np.dot(self.x, self.kValue), axis=2) * self.transitDuty)
        taxCost = np.sum(self.taxPct * self.kValue) + transitDutyCost
        objective = cp.Minimize(transportCost + warehouseCost + taxCost)
        ###constraint###
        constraints = []
        # 1.Goods must be shipped out from its origin to another node and shipped to its destination.
        constraints += [np.sum(self.x[self.kStartPort[k], :, :, k]) == 1 for k in range(self.goods)]
        constraints += [np.sum(self.x[:, self.kEndPort[k], :, k]) == 1 for k in range(self.goods)]
        # 2.For each goods k, it couldn't be shipped out from its destination or shipped to its origin.
        constraints += [np.sum(self.x[:, self.kStartPort[k], :, k]) == 0 for k in range(self.goods)]
        constraints += [np.sum(self.x[self.kEndPort[k], :, :, k]) == 0 for k in range(self.goods)]
        # 3.constraint for transition point
        for k in range(self.goods):
            for j in range(self.portSpace):
                if (j != self.kStartPort[k]) & (j != self.kEndPort[k]):
                    constraints.append(np.sum(self.x[:, j, :, k]) == np.sum(self.x[j, :, :, k]))
        # 4.each goods can only be transitioned in or out of a port for at most once
        constraints += [np.sum(self.x[i, :, :, k]) <= 1 for k in range(self.goods) for i in range(self.portSpace)]
        constraints += [np.sum(self.x[:, j, :, k]) <= 1 for k in range(self.goods) for j in range(self.portSpace)]
        # 5.transition-out should be after transition-in
        constraints += [stayTime[j, k] >= 0 for j in range(self.portSpace) for k in range(self.goods)]
        # 6.constraint for number of containers used
        numCtn = np.dot(self.x, self.kVol) / self.ctnVol
        constraints += [self.y[i, j, t] - numCtn[i, j, t] >= 0 \
                        for i in range(self.portSpace) for j in range(self.portSpace) for t in
                        range(self.dateSpace) if not isinstance(self.y[i, j, t] - numCtn[i, j, t] >= 0, bool)]
        # 7. constraint to check whether a route is used
        constraints += [self.z[i, j, t] >= (np.sum(self.x[i, j, t, :]) * 10e-5) \
                        for i in range(self.portSpace) for j in range(self.portSpace) for t in
                        range(self.dateSpace) if
                        not isinstance(self.z[i, j, t] >= (np.sum(self.x[i, j, t, :]) * 10e-5), bool)]
        # 8.time limitation constraint for each goods
        constraints += [np.sum(arrTime[:, self.kEndPort[k], :, k]) <= self.kDDL[k] for k in range(self.goods)
                        if not isinstance(np.sum(arrTime[:, self.kEndPort[k], :, k]) <= self.kDDL[k], bool)]
        model = cp.Problem(objective, constraints)

        self.objective = objective
        self.constraints = constraints
        self.model = model

    def cplex_build_model(self):
        '''build up the mathematical programming model's objective and constraints using DOCPLEX framework.'''
        model = Model()
        # 4 dimensional binary decision variable matrix
        self.var = model.binary_var_list(self.route_num * self.dateSpace * self.goods, name='x')
        self.x = np.zeros((self.portSpace, self.portSpace, self.dateSpace, self.goods)).astype('object')
        self.x[self.var_location] = self.var
        # 3 dimensional container number matrix
        self.var_2 = model.integer_var_list(self.route_num * self.dateSpace, name='y')
        self.y = np.zeros((self.portSpace, self.portSpace, self.dateSpace)).astype('object')
        self.y[self.var_2_location] = self.var_2
        # 3 dimensional route usage matrix
        self.var_3 = model.binary_var_list(self.route_num * self.dateSpace, name='z')
        self.z = np.zeros((self.portSpace, self.portSpace, self.dateSpace)).astype('object')
        self.z[self.var_3_location] = self.var_3
        # warehouse related cost
        warehouseCost, arrTime, stayTime = self.warehouse_fee(self.x)
        ###objective###
        transportCost = np.sum(self.y * self.tranCost) + np.sum(self.z * self.tranFixedCost)
        transitDutyCost = np.sum(np.sum(np.dot(self.x, self.kValue), axis=2) * self.transitDuty)
        taxCost = np.sum(self.taxPct * self.kValue) + transitDutyCost
        model.minimize(transportCost + warehouseCost + taxCost)
        ###constraint###
        # 1.Goods must be shipped out from its origin to another node and shipped to its destination.
        model.add_constraints(np.sum(self.x[self.kStartPort[k], :, :, k]) == 1 for k in range(self.goods))
        model.add_constraints(np.sum(self.x[:, self.kEndPort[k], :, k]) == 1 for k in range(self.goods))
        # 2.For each goods k, it couldn't be shipped out from its destination or shipped to its origin.
        model.add_constraints(np.sum(self.x[:, self.kStartPort[k], :, k]) == 0 for k in range(self.goods))
        model.add_constraints(np.sum(self.x[self.kEndPort[k], :, :, k]) == 0 for k in range(self.goods))
        # 3.constraint for transition point
        for k in range(self.goods):
            for j in range(self.portSpace):
                if (j != self.kStartPort[k]) & (j != self.kEndPort[k]):
                    model.add_constraint(np.sum(self.x[:, j, :, k]) == np.sum(self.x[j, :, :, k]))
        # 4.each goods can only be transitioned in or out of a port for at most once
        model.add_constraints(np.sum(self.x[i, :, :, k]) <= 1 for k in range(self.goods) for i in range(self.portSpace))
        model.add_constraints(np.sum(self.x[:, j, :, k]) <= 1 for k in range(self.goods) for j in range(self.portSpace))
        # 5.transition-out should be after transition-in
        model.add_constraints(stayTime[j, k] >= 0 for j in range(self.portSpace) for k in range(self.goods))
        # 6.constraint for number of containers used
        numCtn = np.dot(self.x, self.kVol) / self.ctnVol
        model.add_constraints(self.y[i, j, t] - numCtn[i, j, t] >= 0 \
                              for i in range(self.portSpace) for j in range(self.portSpace) for t in
                              range(self.dateSpace) if not isinstance(self.y[i, j, t] - numCtn[i, j, t] >= 0, bool))
        # 7. constraint to check whether a route is used
        model.add_constraints(self.z[i, j, t] >= (np.sum(self.x[i, j, t, :]) * 10e-5) \
                              for i in range(self.portSpace) for j in range(self.portSpace) for t in
                              range(self.dateSpace) if
                              not isinstance(self.z[i, j, t] >= (np.sum(self.x[i, j, t, :]) * 10e-5), bool))
        # 8.time limitation constraint for each goods
        model.add_constraints(np.sum(arrTime[:, self.kEndPort[k], :, k]) <= self.kDDL[k] for k in range(self.goods)
                              if not isinstance(np.sum(arrTime[:, self.kEndPort[k], :, k]) <= self.kDDL[k], bool))

        self.objective = model.objective_expr
        self.constraints = list(model.iter_constraints())
        self.model = model

    def solve_model(self, solver=cp.CBC):
        '''
        solve the optimization model & cache the optimized objective value, route and arrival time for each goods.
        :param solver: the solver to use to solve the LP problem when framework is CVXPY, has no effect to the model
        when framework is DOCPLEX. Default solver is cvxpy.CBC, other open source solvers do not perform that well.
        :return: None
        '''
        try:
            if self.framework == 'CVXPY':
                self.objective_value = self.model.solve(solver)
                self.xs = np.zeros((self.portSpace, self.portSpace, self.dateSpace, self.goods))
                self.xs[self.var_location] = self.var.value
                self.ys = np.zeros((self.portSpace, self.portSpace, self.dateSpace))
                self.ys[self.var_2_location] = self.var_2.value
                self.zs = np.zeros((self.portSpace, self.portSpace, self.dateSpace))
                self.zs[self.var_3_location] = self.var_3.value

            elif self.framework == 'DOCPLEX':
                ms = self.model.solve()
                self.objective_value = self.model.objective_value
                self.xs = np.zeros((self.portSpace, self.portSpace, self.dateSpace, self.goods))
                self.xs[self.var_location] = ms.get_values(self.var)
                self.ys = np.zeros((self.portSpace, self.portSpace, self.dateSpace))
                self.ys[self.var_2_location] = ms.get_values(self.var_2)
                self.zs = np.zeros((self.portSpace, self.portSpace, self.dateSpace))
                self.zs[self.var_3_location] = ms.get_values(self.var_3)

        except:
            raise Exception('Model is not solvable, no solution will be provided')

        nonzeroX = list(zip(*np.nonzero(self.xs)))
        nonzeroX = sorted(nonzeroX, key=lambda x: x[2])
        nonzeroX = sorted(nonzeroX, key=lambda x: x[3])
        nonzeroX = list(map(lambda x: (self.portIndex[x[0]], self.portIndex[x[1]], \
                                       (self.minDate + pd.to_timedelta(x[2], unit='days')).date().isoformat(),
                                       x[3]), nonzeroX))

        self.whCostFinal, arrTime, _ = self.warehouse_fee(self.xs)
        self.transportCost = np.sum(self.ys * self.tranCost) + np.sum(self.zs * self.tranFixedCost)
        self.taxCost = np.sum(self.taxPct * self.kValue) + \
                       np.sum(np.sum(np.dot(self.xs, self.kValue), axis=2) * self.transitDuty)
        self.solution_ = {}
        self.arrTime_ = {}
        for i in range(self.goods):
            self.solution_['goods-' + str(i + 1)] = list(filter(lambda x: x[3] == i, nonzeroX))
            self.arrTime_['goods-' + str(i + 1)] = (self.minDate + pd.to_timedelta \
                (np.sum(arrTime[:, self.kEndPort[i], :, i]), unit='days')).date().isoformat()

    def get_output_(self):
        '''After the model is solved, return total cost, final solution and arrival
        time for each of the goods'''

        return self.objective_value, self.solution_, self.arrTime_

    def warehouse_fee(self, x):
        '''return warehouse fee, arrival time and stay time for each port.'''

        startTime = np.arange(self.dateSpace).reshape(1, 1, self.dateSpace, 1) * x
        arrTimeMtrx = startTime + self.tranTime.reshape(self.portSpace, \
                                                        self.portSpace, self.dateSpace, 1) * x
        arrTime = arrTimeMtrx.copy()
        arrTimeMtrx[:, self.kEndPort.tolist(), :, range(self.goods)] = 0
        stayTime = np.sum(startTime, axis=(1, 2)) - np.sum(arrTimeMtrx, axis=(0, 2))
        stayTime[self.kStartPort.tolist(), range(self.goods)] -= self.kStartTime
        warehouseCost = np.sum(np.sum(stayTime * self.kVol, axis=1) * self.whCost)

        return warehouseCost, arrTime, stayTime

    def txt_solution(self, route, order):
        '''transform the cached results to text.'''

        travelMode = dict(zip(zip(route['Source'], route['Destination']), route['Travel Mode']))
        txt = "Solution"
        txt += "\nNumber of goods: " + str(order['Order Number'].count())
        txt += "\nTotal cost: " + str(self.transportCost + self.whCostFinal + self.taxCost)
        txt += "\nTransportation cost: " + str(self.transportCost)
        txt += "\nWarehouse cost: " + str(self.whCostFinal)
        txt += "\nTax cost: " + str(self.taxCost)

        for i in range(order.shape[0]):
            txt += "\n------------------------------------"
            txt += "\nGoods-" + str(i + 1) + "  Category: " + order['Commodity'][i]
            txt += "\nStart date: " + pd.to_datetime(order['Order Date']) \
                .iloc[i].date().isoformat()
            txt += "\nArrival date: " + str(self.arrTime_['goods-' + str(i + 1)])
            txt += "\nRoute:"
            solution = self.solution_['goods-' + str(i + 1)]
            route_txt = ''
            a = 1
            for j in solution:
                route_txt += "\n(" + str(a) + ")Date: " + j[2]
                route_txt += "  From: " + j[0]
                route_txt += "  To: " + j[1]
                route_txt += "  By: " + travelMode[(j[0], j[1])]
                a += 1
            txt += route_txt

        return txt


def transform(filePath):
    '''Read in order and route data, transform the data into a form that can
    be processed by the operation research model.'''
    order = pd.read_excel(filePath, sheet_name='Order Information')
    route = pd.read_excel(filePath, sheet_name='Route Information')
    order['Tax Percentage'][order['Journey Type'] == 'Domestic'] = 0
    route['Cost'] = route[route.columns[7:12]].sum(axis=1)
    route['Time'] = np.ceil(route[route.columns[14:18]].sum(axis=1) / 24)
    route = route[list(route.columns[0:4]) + ['Fixed Freight Cost', 'Time', \
                                              'Cost', 'Warehouse Cost', 'Travel Mode', 'Transit Duty'] + list(
        route.columns[-9:-2])]
    route = pd.melt(route, id_vars=route.columns[0:10], value_vars=route.columns[-7:] \
                    , var_name='Weekday', value_name='Feasibility')
    route['Weekday'] = route['Weekday'].replace({'Monday': 1, 'Tuesday': 2, 'Wednesday': 3, \
                                                 'Thursday': 4, 'Friday': 5, 'Saturday': 6, 'Sunday': 7})

    return order, route


if __name__ == '__main__':
    order, route = transform("model data.xlsx")
    m = MMT()
    # m = MMT(framework='CVXPY') # for open source framework
    m.set_param(route, order)
    m.build_model()
    m.solve_model()
    txt = m.txt_solution(route, order)
    with open("Solution.txt", "w") as text_file:
        text_file.write(txt)