<gh_stars>0
import gurobipy as gp

m = gp.Model("SP")

C_ij = [[999, 1, 1, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 5, 4, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 7, 3, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 2, 1, 999, 999, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 1, 2, 999, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 5, 4, 999, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 3, 999, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 4, 4, 999, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 2, 4, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 2, 999, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 5, 999, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 2, 8, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 4, 999],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 2],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 1],
        [999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999]]

x = m.addVars(16, 16, vtype=gp.GRB.BINARY, name='x')

m.setObjective(gp.quicksum(C_ij[i][j] * x[i, j] for i in range(16) for j in range(16)),
               gp.GRB.MINIMIZE)

m.addConstrs(
    (gp.quicksum(x[i, j] - x[j, i] for j in range(16)) == 1) for i in range(1))
m.addConstrs(
    (gp.quicksum(x[i, j] - x[j, i] for j in range(16)) == -1) for i in range(15, 16))
m.addConstrs(
    (gp.quicksum(x[i, j] - x[j, i] for j in range(16)) == 0) for i in range(1, 15))
m.optimize()
print("\n###########################\nObj:", m.objVal)
