Смешанное целочисленное программирование - расположение склада (Python + GLPK)

Question

Я относительно новичок в оптимизации, и я пытаюсь оптимизировать проблему (из прошлого курса в Coursera, 2 года назад) о местоположении склада. Проблема в том, что прошло более 6 часов, и он все еще работает на экземпляре со 100 складами и 1000 клиентов.

Проблема в следующем. У меня есть набор складов, которые я могу открыть или нет. Открытие каждого из них имеет стоимость s_w. Кроме того, все они имеют максимальный размер cap_w. С другой стороны, есть группа клиентов, все они должны быть подключены к одному (и только одному) открытому складу. У каждого из них есть спрос d_c, а для каждого клиента существует стоимость транспортировки с каждого склада t_wc. То, что я хочу, это, безусловно, минимизировать общую стоимость.

Итак, у меня есть массив размером, равным общему количеству складов под названием x. Каждое x[w] является целым числом {0,1}, определяющим, открыт склад w или нет. У меня также есть матрица 0 и 1, определяющая, какой склад доставляет каждому клиенту. следовательно, имеется столько строк, сколько клиентов, и столько столбцов, сколько складов. Матрица называется у. y[c][w] равно 1, если товар w доставляет клиенту c, 0 в противном случае.

Все идет нормально. Предполагается, что это проблема MIP. Чтобы сделать это, я делаю это на Python, используя библиотеку PuPL ( https://pythonhosted.org/PuLP/pulp.html) и GLPK для ее решения.

Теперь вот моя модель:

#!/usr/bin/python
# -*- coding: utf-8 -*-

from pulp import *

def solveIt(inputData):

# parse the input
lines = inputData.split('\n')

parts = lines[0].split()
warehouseCount = int(parts[0])
customerCount = int(parts[1])

warehouses = []
for i in range(1, warehouseCount+1):
    line = lines[i]
    parts = line.split()
    warehouses.append((int(parts[0]), float(parts[1])))

customerDemands = []
customerCosts = []

lineIndex = warehouseCount+1
for i in range(0, customerCount):
    customerDemand = int(lines[lineIndex+2*i])
    customerCost = map(float, lines[lineIndex+2*i+1].split())
    customerDemands.append(customerDemand)
    customerCosts.append(customerCost)



x = [LpVariable("x"+str(w),0,1,cat='Integer') for w in range(0,warehouseCount)]
y = [[LpVariable("y"+str(c)+","+str(w),0,1,cat='Integer') for w in range(0,warehouseCount)] for c in range(0,customerCount)]

prob = LpProblem("Warehouse Location",LpMinimize)

#Constraints
# y_cw <= x_w makes sure that no client is delivered by a closed warehouse
for w in range(0,warehouseCount):
    for c in range(0,customerCount):
        prob += y[c][w] <= x[w]

#A client is served by exactly one warehouse
for c in range(0,customerCount):
    affineExpression = []
    for w in range(0,warehouseCount):
        affineExpression.append((y[c][w],1))
    prob += LpAffineExpression(affineExpression) == 1

#For each warehouse, the sum of demand of all the clients it serves is lower than its capacity
for w in range(0,warehouseCount):
    affineExpression = []
    for c in range(0,customerCount):
        affineExpression.append((y[c][w],customerDemands[c]))
    prob += LpAffineExpression(affineExpression) <= warehouses[w][0]

#Objective
#The sum of all the warehouses opening plus the transportation costs has to be minimal
affineExpression = []
for w in range(0,warehouseCount):
    affineExpression.append((x[w],warehouses[w][1]))
    for c in range(0,customerCount):
        affineExpression.append((y[c][w],customerCosts[c][w]))

prob += LpAffineExpression(affineExpression)

print "#######################START SOLVING"
status = prob.solve(GLPK(mip=1,msg = 1))
print LpStatus[status]
for w in range(0,warehouseCount):
    print value(x[w])

solution = []
for c in range(0,customerCount):
    string = ""
    whichOne = -1
    for w in range(0,warehouseCount):
        string += str(value(y[c][w])) + " "
        if value(y[c][w]) == 1:
            whichOne = w
            solution.append(w)
    print string+ "   "+str(whichOne)


# calculate the cost of the solution
obj = sum([warehouses[x][1]*x[w] for x in range(0,warehouseCount)])
for c in range(0, customerCount):
    obj += customerCosts[c][solution[c]]

# prepare the solution in the specified output format
outputData = str(obj) + ' ' + str(0) + '\n'
outputData += ' '.join(map(str, solution))

return outputData

Я знаю, что способ построения матрицы не оптимален, но на самом деле это не занимает много времени. Это начало решаться, и в какой-то момент я достиг точки, где GLPK сказал:

OPTIMAL SOLUTION FOUND
Integer optimization begins...

Я полагаю, что это означает, что он решил LP, и теперь он получил целое число... но прошло около 6 часов или около того, и он прогрессирует, и все еще есть, но не заканчивает. В меньших случаях это работало нормально.

Наверное, мой вопрос... Что-то не так с моей моделью? Некоторые оптимизации я забыл? ИЛИ эта проблема такая огромная?

Кроме того, что касается компьютера, он довольно плохой: только Intel Atom и 1 ГБ оперативной памяти...

Спасибо за помощь!

РЕДАКТИРОВАТЬ: Вот дата: https://github.com/ddeunagomez/DiscreteOptimization/blob/master/04_warehouse_location/warehouse/data/wl_100_1 формат:

NumberOfWarehouses NumberOfCustomers
CapacityWarehouse1 OpeningCostWarehouse1
CapacityWarehouse2 OpeningCostWarehouse2
.....
CapacityWarehouseN OpeningCostWarehouseN
DemandCustomer1
TransportCostW1_C1 TransportCostW2_C1 ....... TransportCostWN_C1
DemandCustomer2
TransportCostW1_C2 TransportCostW2_C2 ....... TransportCostWN_C2
.....
DemandCustomerN
TransportCostW1_CM TransportCostW2_CM ....... TransportCostWN_CM

Answer 1

Это не очень большая проблема в схеме вещей - есть специализированный код для решения гораздо больших экземпляров хранилища, чем этот, и хорошие решения, такие как CPLEX, тоже могут легко его решить. Я не знаю, насколько эффективны GLPK/PuPL, но вполне может быть, что они просто слишком долго используют прямой LP/branch-and-bound (что они и делают).

Одна вещь, которую вы можете попробовать, это позволить переменным y быть непрерывными (0 <= y <= 1) вместо двоичных. Это, скорее всего, ускорит время выполнения, потому что решателю не придется переходить на них. Физическая интерпретация заключается в том, что некоторые клиенты могут разделить свои требования между несколькими складами. На практике очень немногие из них, вероятно, будут разделены на большинство решений. Если емкости достаточно велики, чтобы не связываться, ни одно из требований не будет разделено, и вы всегда получите бинарные решения, даже если разрешите y быть непрерывным.