备份 多产品库存路径问题python+docplex代码

from docplex.mp.conflict_refiner import ConflictRefiner
from docplex.mp.conflict_refiner import ConflictRefinerResult
import numpy as np
import csv
import pandas as pd
import math
# 导入库
from docplex.mp.model import Model
from docplex.util.environment import get_environment


class Node():
    def __init__(self):
        self.node_no = 0
        self.x_coord = 0.
        self.y_coord = 0.
        self.init_inv = []
        self.upper = 0
        self.lower = 0.
        self.demands = []
        self.hi = 0.


def calDistance(node_list, depot):
    '''
    计算节点之间距离
    输入：node_list-node列表以及depot
    输出：距离矩阵-dis_matrix 用dis_matrix.iloc[i][j]读 0是depot 1-n 是零售商
    '''
    dis_matrix = pd.DataFrame(data=None, columns=range(len(node_list) + 1), index=range(len(node_list) + 1))
    x0, y0 = depot.x_coord, depot.y_coord
    dis_matrix.iloc[0, 0] = 0
    for i in range(len(node_list)):
        xi, yi = node_list[i].x_coord, node_list[i].y_coord
        dis_matrix.iloc[i + 1, 0] = np.around(math.sqrt((xi - x0) ** 2 + (yi - y0) ** 2), 2)
        dis_matrix.iloc[0, i + 1] = dis_matrix.iloc[i + 1, 0]
    for i in range(len(node_list)):
        xi, yi = node_list[i].x_coord, node_list[i].y_coord
        for j in range(len(node_list)):
            xj, yj = node_list[j].x_coord, node_list[j].y_coord
            dis_matrix.iloc[i + 1, j + 1] = np.around(math.sqrt((xi - xj) ** 2 + (yi - yj) ** 2), 2)
    return dis_matrix


def read_csv_file(filepath):
    node_list = []
    with open(filepath) as absXnY:
        f_csv = csv.reader(absXnY)
        # header = next(f_csv)
        i = -2
        for line in f_csv:
            if len(line) == 3:
                retailer_num = int(line[0]) - 1
                period = int(line[1])
                veh_cap = float(line[2])
            elif int(line[0]) == 1:
                depot = Node()
                depot.node_no = int(line[0]) - 1
                depot.x_coord = float(line[1])
                depot.y_coord = float(line[2])
                depot.init_inv = eval(line[3])
                depot.upper = float(line[4])
                depot.hi = float(line[5])
                depot.lower = None
                depot.demands = []
            else:
                node_list.append(Node())
                node_list[i].node_no = int(line[0]) - 1
                node_list[i].x_coord = float(line[1])
                node_list[i].y_coord = float(line[2])
                node_list[i].init_inv = eval(line[3])
                node_list[i].upper = float(line[4])
                node_list[i].lower = float(line[5])
                node_list[i].demands = eval(line[6])
                node_list[i].hi = float(line[7])
            i += 1
    # print(len(node_list))
    # print(node_list[1].node_no)
    # print(node_list[1].init_inv)
    # print(node_list[3].demands)
    return node_list, depot, retailer_num, period, veh_cap


def build_MIRP(mdl, node_list, dis_matrix, retailer_num, period, veh_cap, veh_num, prod_num, **kwargs):
    """ Takes as input:
        - a list of product tuples (name, demand, inside, outside)
        - a list of resource tuples (name, capacity)
        - a list of consumption tuples (product_name, resource_named, consumed)
    """
    # --- decision variables ---
    x = []
    y = []
    v = []
    q = []
    z = []
    w = []
    for t in range(period):
        x.append([])
        for k in range(veh_num):
            x[t].append([])
            for j in range(retailer_num):
                x[t][k].append([])
                x[t][k][j].append(mdl.integer_var(lb=0, ub=2, name='x{0}{1}{2}{3}'.format(0, j, k, t)))
                for i in range(1, retailer_num):
                    x[t][k][j].append(mdl.binary_var(name='x{0}{1}{2}{3}'.format(i, j, k, t)))
    for t in range(period):
        y.append([])
        for k in range(veh_num):
            y[t].append([])
            for i in range(retailer_num):
                y[t][k].append(mdl.binary_var(name='y{0}{1}{2}'.format(i, k, t)))
    for t in range(period):
        v.append([])
        for m in range(prod_num):
            v[t].append([])
            for i in range(retailer_num):
                v[t][m].append(mdl.integer_var(lb=0, name='v{0}{1}{2}'.format(i, m, t)))
    for t in range(period):
        q.append([])
        for k in range(veh_num):
            q[t].append([])
            for m in range(prod_num):
                q[t][k].append([])
                for i in range(retailer_num):
                    q[t][k][m].append(mdl.integer_var(lb=0, name='q{0}{1}{2}{3}'.format(i, m, k, t)))
    for t in range(period):
        z.append([])
        for j in range(1, retailer_num):
            if j == 1:
                z[t].append([])
            z[t].append([])
            for i in range(retailer_num):
                z[t][j].append(mdl.binary_var(name='z{0}{1}{2}'.format(i, j, t)))
    for t in range(period):
        w.append([])
        for m in range(prod_num):
            w[t].append([])
            for j in range(retailer_num):
                w[t][m].append([])
                w[t][m][j].append(mdl.continuous_var(lb=0, name='x{0}{1}{2}{3}'.format(0, j, m, t)))
                for i in range(1, retailer_num):
                    if i == 1:
                        w[t][m][j].append([])
                    w[t][m][j].append(mdl.binary_var(name='x{0}{1}{2}{3}'.format(i, j, m, t)))
    # mdl.x_vars = x
    # mdl.y_vars = y
    # mdl.v_vars = v
    # mdl.q_vars = q
    # mdl.z_vars = z
    # mdl.w_vars = w
    # --- constraints ---
    # demand satisfaction
    # (2)
    # mdl.add_constraints((mdl.v_vars[t][m][0] == mdl.v_vars[t - 1][m][0] + M - mdl.sum(
    #     mdl.q_vars[t][k][m][i] for k in range(veh_num) for i in range(1, retailer_num)) - mdl.sum(
    #     mdl.w_vars[t][m][i][0]) for i in range(1, retailer_num)) for m in range(prod_num) for t in range(1, period))

    mdl.add_constraints(v[t][m][0] == v[t - 1][m][0] + M - mdl.sum(
        q[t][k][m][i] for k in range(veh_num) for i in range(1, retailer_num)) - mdl.sum(
        w[t][m][i][0]) for i in range(1, retailer_num) for m in range(prod_num) for t in range(1, period))


    # (3)
    mdl.add_constraints((mdl.v_vars[t][m][i] == mdl.v_vars[t - 1][m][i] + mdl.sum(
        mdl.q_vars[t][k][m][i] for k in range(veh_num)) + mdl.sum(
        mdl.w_vars[t][m][i][j] for j in range(retailer_num)) - mdl.sum(
        mdl.w_vars[t][m][j][i] for j in range(1, retailer_num)) - node_list[i].demands[m]) for m in range(prod_num) for
                        i in range(1, retailer_num) for t in range(1, period))
    # (4)
    mdl.add_constraints((mdl.sum(mdl.v_vars[t][m][i] for m in range(prod_num)) <= node_list[i].upper
                         for i in range(1, retailer_num)) for t in range(period))
    # (5)
    mdl.add_constraints((mdl.sum(mdl.q_vars[t][k][m][i] for k in range(veh_num) for m in range(prod_num)) <=
                         node_list[i].upper - mdl.sum(mdl.v_vars[t - 1][m][i] for m in range(prod_num)) - mdl.sum(
                mdl.w_vars[t][m][j][i] for j in range(1, retailer_num) for m in range(prod_num))) for i in
                        range(1, retailer_num) for t in range(1, period))
    # (6)
    mdl.add_constraints(
        (mdl.sum(mdl.q_vars[t][k][m][i] for m in range(prod_num)) <= mdl.y_vars[t][k][i] * node_list[i].upper) for i in
        range(1, retailer_num) for m in range(prod_num) for k in range(veh_num) for t in range(period))
    # (7)
    mdl.add_constraints((mdl.sum(
        mdl.q_vars[t][k][m][i] for m in range(prod_num) for i in range(1, retailer_num)) <= veh_cap * mdl.y_vars[t][k][
                             0]) for k in range(veh_num) for t in range(period))
    # (8) ×
    # (9)
    mdl.add_constraints(
        (mdl.v_vars[t - 1][i] + mdl.sum(mdl.q_vars[t][k][m][i] for k in range(veh_num)) <= node_list[i].upper) for m in
        range(prod_num) for i in range(1, retailer_num) for t in range(1, period))
    # (10) ×
    # (11)
    mdl.add_constraints((node_list[i].demands[m] - mdl.v_vars[t - 1][m][i] - mdl.sum(
        mdl.q_vars[t][k][m][i] for k in range(veh_num)) - M * (1 - mdl.z_vars[t][j][i]) < mdl.v_vars[t - 1][m][
                             j] + mdl.sum(mdl.q_vars[t][k][m][j] for k in range(veh_num)) - node_list[j].demands[
                             m] <= M * (1 - mdl.z_vars[t][j][i])) for m in range(prod_num) for i in
                        range(1, retailer_num) for j in range(1, retailer_num) for t in range(1, period))
    # (12)
    mdl.add_constraints((mdl.v_vars[t - 1][m][j] + mdl.sum(mdl.q_vars[t][k][m][j] for k in range(veh_num)) -
                         node_list[j].demands[m] <= M * (1 - mdl.z_vars[t][j][0])) for m in range(prod_num) for j in
                        range(1, retailer_num) for t in range(period))
    # (13)
    mdl.add_constraints((mdl.w_vars[t][m][j][i] <= mdl.v_vars[t - 1][m][i] + mdl.sum(
        mdl.q_vars[t][k][m][i] for k in range(veh_num)) - node_list[i].demands[m] + M * (1 - mdl.z_vars[t][j][i])) for m
                        in range(prod_num) for i in range(1, retailer_num) for j in range(1, retailer_num) for t in
                        range(period))
    # (14) --
    # (15) --
    # (16) --

    # --- objective ---
    mdl.total_routing_cost = mdl.sum(
        dis_matrix.iloc[i, j] * mdl.x_vars[t][k][j][i] for t in range(period) for k in range(veh_num) for i in
        range(retailer_num) for j in range(retailer_num)) + mdl.sum(
        dis_matrix.iloc[i, j] * mdl.z_vars[t][j][z] for t in range(period) for i in range(retailer_num) for j in
        range(retailer_num))
    mdl.add_kpi(mdl.total_routing_cost, "routing cost")
    mdl.total_shipping_quantity = mdl.sum(
        mdl.q_vars[t][k][m][i] for t in range(period) for k in range(veh_num) for m in range(prod_num) for i in
        range(1, retailer_num)) + mdl.sum(
        mdl.w_vars[t][m][j][i] for t in range(period) for m in range(prod_num) for j in range(retailer_num) for i in
        range(retailer_num))
    mdl.add_kpi(mdl.total_shipping_quantity, "shipping quantity")
    mdl.minimize(mdl.total_routing_cost)

    return mdl


def print_production_solution(mdl, retailer_num, veh_num, period):
    obj = mdl.objective_value
    print("* 目标函数值: {:g}".format(obj))
    print("* 总路径成本{0}".format(mdl.total_routing_cost))
    for i in range(retailer_num):
        for j in range(retailer_num):
            for k in range(veh_num):
                for t in range(period):
                    print("x{0}{1}{2}{3}:".format(i, j, k, t) + str(mdl.x_vars[t][k][j][i]))


if __name__ == '__main__':
    # Build the model
    node_list, depot, retailer_num, period, veh_cap = read_csv_file('G1H3n5lessP.csv')
    veh_num = 2
    prod_num = 5
    M = 9999
    dis_matrix = calDistance(node_list, depot)
    with Model(name='MIRPLR') as model:
        model = build_MIRP(model, node_list, dis_matrix, retailer_num, period, veh_cap, veh_num, prod_num)
        model.print_information()
        # Solve the model.
        if model.solve():
            print_production_solution(model, retailer_num, veh_num, period)
            # Save the CPLEX solution as "solution.json" program output
            with get_environment().get_output_stream("solution.json") as fp:
                model.solution.export(fp, "json")
        else:
            print("Problem has no solution")