weber_soco.py

Go to the documentation of this file.

 
"""
Copyright (c) by Joao Pedro PEDROSO, Masahiro MURAMATSU and Mikio KUBO, 2012
"""
from pyscipopt import Model, quicksum
 
 
def weber(I, x, y, w):
    """weber: model for solving the single source weber problem using soco.
    Parameters:
        - I: set of customers
        - x[i]: x position of customer i
        - y[i]: y position of customer i
        - w[i]: weight of customer i
    Returns a model, ready to be solved.
    """
 
    model = Model("weber")
 
    X, Y, z, xaux, yaux = {}, {}, {}, {}, {}
    X = model.addVar(lb=-model.infinity(), vtype="C", name="X")
    Y = model.addVar(lb=-model.infinity(), vtype="C", name="Y")
    for i in I:
        z[i] = model.addVar(vtype="C", name="z(%s)" % (i))
        xaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s)" % (i))
        yaux[i] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s)" % (i))
 
    for i in I:
        model.addCons(xaux[i] * xaux[i] + yaux[i] * yaux[i] <= z[i] * z[i], "MinDist(%s)" % (i))
        model.addCons(xaux[i] == (x[i] - X), "xAux(%s)" % (i))
        model.addCons(yaux[i] == (y[i] - Y), "yAux(%s)" % (i))
 
    model.setObjective(quicksum(w[i] * z[i] for i in I), "minimize")
 
    model.data = X, Y, z
    return model
 
 
import random
 
 
def make_data(n, m):
    """creates example data set"""
    I = range(1, n + 1)
    J = range(1, m + 1)
    x, y, w = {}, {}, {}
    for i in I:
        x[i] = random.randint(0, 100)
        y[i] = random.randint(0, 100)
        w[i] = random.randint(1, 5)
    return I, J, x, y, w
 
 
if __name__ == "__main__":
 
    random.seed(3)
    n = 7
    m = 1
    I, J, x, y, w = make_data(n, m)
    print("data:")
    print("%s\t%8s\t%8s\t%8s" % ("i", "x[i]", "y[i]", "w[i]"))
    for i in I:
        print("%s\t%8g\t%8g\t%8g" % (i, x[i], y[i], w[i]))
    print
 
    model = weber(I, x, y, w)
    model.optimize()
    X, Y, z = model.data
    print("Optimal value=", model.getObjVal())
    print("Selected position:", )
    print("\t", (round(model.getVal(X)), round(model.getVal(Y))))
    print
    print("Solution:")
    print("%s\t%8s" % ("i", "z[i]"))
    for i in I:
        print("%s\t%8g" % (i, model.getVal(z[i])))
    print
    try:  # plot the result using networkx and matplotlib
        import networkx as NX
        import matplotlib.pyplot as P
 
        P.clf()
        G = NX.Graph()
 
        G.add_nodes_from(I)
        G.add_nodes_from(["D"])
 
        position = {}
        for i in I:
            position[i] = (x[i], y[i])
        position["D"] = (round(model.getVal(X)), round(model.getVal(Y)))
 
        NX.draw(G, pos=position, node_size=200, node_color="g", nodelist=I)
        NX.draw(G, pos=position, node_size=400, node_color="w", nodelist=["D"], alpha=0.5)
        # P.savefig("weber.pdf",format="pdf",dpi=300)
        P.show()
 
    except ImportError:
        print("install 'networkx' and 'matplotlib' for plotting")
 
 
def weber_MS(I, J, x, y, w):
    """weber -- model for solving the weber problem using soco (multiple source version).
    Parameters:
        - I: set of customers
        - J: set of potential facilities
        - x[i]: x position of customer i
        - y[i]: y position of customer i
        - w[i]: weight of customer i
    Returns a model, ready to be solved.
    """
    M = max([((x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2) for i in I for j in I])
    model = Model("weber - multiple source")
    X, Y, v, u = {}, {}, {}, {}
    xaux, yaux, uaux = {}, {}, {}
    for j in J:
        X[j] = model.addVar(lb=-model.infinity(), vtype="C", name="X(%s)" % j)
        Y[j] = model.addVar(lb=-model.infinity(), vtype="C", name="Y(%s)" % j)
        for i in I:
            v[i, j] = model.addVar(vtype="C", name="v(%s,%s)" % (i, j))
            u[i, j] = model.addVar(vtype="B", name="u(%s,%s)" % (i, j))
            xaux[i, j] = model.addVar(lb=-model.infinity(), vtype="C", name="xaux(%s,%s)" % (i, j))
            yaux[i, j] = model.addVar(lb=-model.infinity(), vtype="C", name="yaux(%s,%s)" % (i, j))
            uaux[i, j] = model.addVar(vtype="C", name="uaux(%s,%s)" % (i, j))
 
    for i in I:
        model.addCons(quicksum(u[i, j] for j in J) == 1, "Assign(%s)" % i)
        for j in J:
            model.addCons(xaux[i, j] * xaux[i, j] + yaux[i, j] * yaux[i, j] <= v[i, j] * v[i, j],
                          "MinDist(%s,%s)" % (i, j))
            model.addCons(xaux[i, j] == (x[i] - X[j]), "xAux(%s,%s)" % (i, j))
            model.addCons(yaux[i, j] == (y[i] - Y[j]), "yAux(%s,%s)" % (i, j))
            model.addCons(uaux[i, j] >= v[i, j] - M * (1 - u[i, j]), "uAux(%s,%s)" % (i, j))
 
    model.setObjective(quicksum(w[i] * uaux[i, j] for i in I for j in J), "minimize")
 
    model.data = X, Y, v, u
    return model
 
 
if __name__ == "__main__":
    random.seed(3)
    n = 7
    m = 1
    I, J, x, y, w = make_data(n, m)
    model = weber_MS(I, J, x, y, w)
    model.optimize()
    X, Y, w, z = model.data
    print("Optimal value:", model.getObjVal())
    print("Selected positions:")
    for j in J:
        print("\t", (model.getVal(X[j]), model.getVal(Y[j])))
    for (i, j) in sorted(w.keys()):
        print("\t", (i, j), model.getVal(w[i, j]), model.getVal(z[i, j]))
 
    EPS = 1.e-4
    edges = [(i, j) for (i, j) in z if model.getVal(z[i, j]) > EPS]
 
    try:  # plot the result using networkx and matplotlib
        import networkx as NX
        import matplotlib.pyplot as P
 
        P.clf()
        G = NX.Graph()
 
        G.add_nodes_from(I)
        G.add_nodes_from("%s" % j for j in J)  # for distinguishing J from I, make nodes as strings
        for (i, j) in edges:
            G.add_edge(i, "%s" % j)
 
        position = {}
        for i in I:
            position[i] = (x[i], y[i])
        for j in J:
            position["%s" % j] = (model.getVal(X[j]), model.getVal(Y[j]))
        print(position)
 
        NX.draw(G, position, node_color="g", nodelist=I)
        NX.draw(G, position, node_color="w", nodelist=["%s" % j for j in J])
        P.show()
    except ImportError:
        print("install 'networkx' and 'matplotlib' for plotting")