docplex官方tutorial搬运1 - 选址问题

描述问题

一家虚构的Coffee Company计划在不久的将来开设N家商店，并且需要了解以下几点，以确定他们应位于的位置：
- 这个咖啡屋的大多数客户喜欢阅读和借书，因此目标是要以使所有城市公共图书馆都在最小步行距离之内的方式来定位这些商店。
在此示例中，我们使用芝加哥公开数据。
我们采用K-Median模型来获得未来商店的最佳位置。

Use decision optimization

Step 1: Import the library¶

Run the following code to import the Decision Optimization CPLEX Modeling library. The DOcplex 库包括两个模型包，分别是数学规划和约束规划（我这里不完全确定Mathematical Programming and Constraint Programming的意思）

import sys
try:
    import docplex.mp
except:
    raise Exception('Please install docplex. See https://pypi.org/project/docplex/')

(对于这个问题用的是docplex.mp）

Step 2: Model the data

这里用到的数据十分简单，包括公共图书馆的列表以及他们的地理位置信息。
数据源来自 Chicago open data ，数据是json格式的，如下: data" : [ [ 1, "13BFA4C7-78CE-4D83-B53D-B57C60B701CF", 1, 1441918880, "885709", 1441918880, "885709", null, "Albany Park", "M, W: 10AM-6PM; TU, TH: 12PM-8PM; F, SA: 9AM-5PM; SU: Closed", "Yes", "Yes ", "3401 W. Foster Avenue", "CHICAGO", "IL", "60625", "(773) 539-5450", [ "http://www.chipublib.org/locations/1/", null ], [ null, "41.975456", "-87.71409", null, false ] ]

Step 3: Prepare the data

我们需要收集各个地点的经纬度数据

# Store longitude, latitude and street crossing name of each public library location.
class XPoint(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __str__(self):
        return "P(%g_%g)" % (self.x, self.y)

class NamedPoint(XPoint):
    def __init__(self, name, x, y):
        XPoint.__init__(self, x, y)
        self.name = name
    def __str__(self):
        return self.name

定义计算距离的函数

为了简化计算，直接使用python的第三方package geopy

try:
    import geopy.distance
except:
    if hasattr(sys, 'real_prefix'):
        #we are in a virtual env.
        !pip install geopy 
    else:
        !pip install --user geopy

# Simple distance computation between 2 locations.
from geopy.distance import great_circle
 
def get_distance(p1, p2):
    return great_circle((p1.y, p1.x), (p2.y, p2.x)).miles

声明图书馆数据列表

将json数据进行处理

def build_libraries_from_url(url):
    import requests
    import json
    from six import iteritems

    r = requests.get(url)
    myjson = json.loads(r.text, parse_constant='utf-8')
    
    # find columns for name and location
    columns = myjson['meta']['view']['columns']
    name_col = -1
    location_col = -1
    for (i, col) in enumerate(columns):
        if col['name'].strip().lower() == 'name':
            name_col = i
        if col['name'].strip().lower() == 'location':
            location_col = i
    if (name_col == -1 or location_col == -1):
        raise RuntimeError("Could not find name and location columns in data. Maybe format of %s changed?" % url)
    
    # get library list
    data = myjson['data']

    libraries = []
    k = 1
    for location in data:
        uname = location[name_col]
        try:
            latitude = float(location[location_col][1])
            longitude = float(location[location_col][2])
        except TypeError:
            latitude = longitude = None
        try:
            name = str(uname)
        except:
            name = "???"
        name = "P_%s_%d" % (name, k)
        if latitude and longitude:
            cp = NamedPoint(name, longitude, latitude)
            libraries.append(cp)
            k += 1
    return libraries

libraries = build_libraries_from_url('https://data.cityofchicago.org/api/views/x8fc-8rcq/rows.json?accessType=DOWNLOAD')
print("There are %d public libraries in Chicago" % (len(libraries)))

定义咖啡屋的open数量

定义一个常量表示我们需要开几个咖啡屋
nb_shops = 5
print("We would like to open %d coffee shops" % nb_shops)

验证一下数据

通过folium这个python第三方工具

try:
    import folium
except:
    if hasattr(sys, 'real_prefix'):
        #we are in a virtual env.
        !pip install folium 
    else:
        !pip install --user folium

import folium
map_osm = folium.Map(location=[41.878, -87.629], zoom_start=11)
for library in libraries:
    lt = library.y
    lg = library.x
    folium.Marker([lt, lg]).add_to(map_osm)
map_osm

运行上述代码后，将显示数据，但是仅通过查看地图就无法确定理想情况下在哪里开设咖啡店。
让我们设置DOcplex来编写和解决一个优化模型，该模型将帮助我们确定以最佳方式在哪里放置咖啡店。
所以接下来才是重点...

Step 4: 定义模型

from docplex.mp.environment import Environment
env = Environment()
env.print_information()

创建一个docplex模型

(以下废话太多了，我简单合并一下)

from docplex.mp.model import Model

mdl = Model("coffee shops")
BIGNUM = 999999999

# 保证点唯一
libraries = set(libraries)
# 考虑简单的情况，咖啡屋的候选点和图书馆的数量一致
# 也就是说，每个图书馆的位置都是我们的候选点
coffeeshop_locations = libraries

# 决策变量
# 0-1决策变量表示那个候选点将被选择
coffeeshop_vars = mdl.binary_var_dict(coffeeshop_locations, name="is_coffeeshop")
# 这个决策变量表示两点是否被连接
link_vars = mdl.binary_var_matrix(coffeeshop_locations, libraries, "link")

# First constraint: if the distance is suspect, it needs to be excluded from the problem.
# 第一个约束：如果距离超过最大值，排除
for c_loc in coffeeshop_locations:
    for b in libraries:
        if get_distance(c_loc, b) >= BIGNUM:
            mdl.add_constraint(link_vars[c_loc, b] == 0, "ct_forbid_{0!s}_{1!s}".format(c_loc, b))

# Second constraint: each library must be linked to a coffee shop that is open.
# 第二个约束：每一个图书馆必须连接一个咖啡屋
mdl.add_constraints(link_vars[c_loc, b] <= coffeeshop_vars[c_loc]
                   for b in libraries
                   for c_loc in coffeeshop_locations)
mdl.print_information()

# Third constraint: each library is linked to exactly one coffee shop.
# 第三个约束：每一个图书馆只能连接一个咖啡屋
mdl.add_constraints(mdl.sum(link_vars[c_loc, b] for c_loc in coffeeshop_locations) == 1
                   for b in libraries)
mdl.print_information()

# Fourth constraint: there is a fixed number of coffee shops to open.
# 第四个约束：咖啡屋的数量是确定的
# Total nb of open coffee shops
mdl.add_constraint(mdl.sum(coffeeshop_vars[c_loc] for c_loc in coffeeshop_locations) == nb_shops)

# Print model information
mdl.print_information()

# 目标函数：最小化距离
# Minimize total distance from points to hubs
total_distance = mdl.sum(link_vars[c_loc, b] * get_distance(c_loc, b) for c_loc in coffeeshop_locations for b in libraries)
mdl.minimize(total_distance)

# 求解
print("# coffee shops locations = %d" % len(coffeeshop_locations))
print("# coffee shops           = %d" % nb_shops)

assert mdl.solve(), "!!! Solve of the model fails"

第五步：查看求解结果，进行案例分析

The solution can be analyzed by displaying the location of the coffee shops on a map.

total_distance = mdl.objective_value
open_coffeeshops = [c_loc for c_loc in coffeeshop_locations if coffeeshop_vars[c_loc].solution_value == 1]
not_coffeeshops = [c_loc for c_loc in coffeeshop_locations if c_loc not in open_coffeeshops]
edges = [(c_loc, b) for b in libraries for c_loc in coffeeshop_locations if int(link_vars[c_loc, b]) == 1]

print("Total distance = %g" % total_distance)
print("# coffee shops  = {0}".format(len(open_coffeeshops)))
for c in open_coffeeshops:
    print("new coffee shop: {0!s}".format(c))

# 展示结果，咖啡屋将被标红
import folium
map_osm = folium.Map(location=[41.878, -87.629], zoom_start=11)
for coffeeshop in open_coffeeshops:
    lt = coffeeshop.y
    lg = coffeeshop.x
    folium.Marker([lt, lg], icon=folium.Icon(color='red',icon='info-sign')).add_to(map_osm)
    
for b in libraries:
    if b not in open_coffeeshops:
        lt = b.y
        lg = b.x
        folium.Marker([lt, lg]).add_to(map_osm)
    

for (c, b) in edges:
    coordinates = [[c.y, c.x], [b.y, b.x]]
    map_osm.add_child(folium.PolyLine(coordinates, color='#FF0000', weight=5))

map_osm