Skip to content

Instantly share code, notes, and snippets.

@ljwolf

ljwolf/ortools_toy_optimization_models.ipynb

Forked from jGaboardi/ortools_toy_optimization_models.ipynb
Created August 13, 2018 18:37
Show Gist options
  • Save ljwolf/eef3cf1617becae24e606e684787ffdc to your computer and use it in GitHub Desktop.
Save ljwolf/eef3cf1617becae24e606e684787ffdc to your computer and use it in GitHub Desktop.
Download ZIP
7x5 OD facility location models with ortools
Display the source blob
Display the rendered blob
Raw
ortools_toy_optimization_models.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 7x5 OD facility location models with [`ortools`](https://developers.google.com/optimization/)\n",
"\n",
"#### James Gaboardi <[email protected]>, 08/2018"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from ortools.linear_solver import pywraplp\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate synthetic data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(1991)\n",
"\n",
"# n clients and n facilities\n",
"client_count, facility_count = 7, 5\n",
"\n",
"# client demand (weights)\n",
"weights = np.random.randint(1, 10, (client_count,1))\n",
"\n",
"# distance matrix\n",
"cost_matrix = np.random.uniform(0.1, 4.,\n",
" (client_count,facility_count))\n",
"cost_matrix = np.round(cost_matrix, decimals=2)\n",
"\n",
"# candidate facilites to site\n",
"facilities = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---------------------------"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Location Set Covering Problem"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def lscp(name, cij=None, s=None):\n",
" \"\"\"\n",
" Location Set Covering Problem\n",
" \n",
" Originally Published:\n",
" Toregas, C. and ReVelle, Charles. 1972.\n",
" Optimal Location Under Time or Distance Constraints.\n",
" Papers of the Regional Science Association. 28(1):133 - 144.\n",
" \n",
" Parameters\n",
" ----------\n",
" name : str\n",
" model name\n",
" cij : numpy.ndarray\n",
" cost matrix\n",
" s : float\n",
" service radius\n",
" \"\"\"\n",
" # client count and range\n",
" n_cli, r_cli = cij.shape[0], xrange(cij.shape[0])\n",
" \n",
" # facility count and range\n",
" n_fac, r_fac = cij.shape[1], xrange(cij.shape[1])\n",
" \n",
" # binary coverage matrix from cij\n",
" aij = np.zeros(cij.shape)\n",
" aij[cij <= s] = 1.\n",
" \n",
" # create a solver instance\n",
" solver_instance = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING\n",
" \n",
" # instantiate a model\n",
" model = pywraplp.Solver(name, solver_instance)\n",
" \n",
" # facility decision variables\n",
" fac_vars = {j: model.IntVar(0,1, 'y[%i]' % (j)) for j in r_fac}\n",
" \n",
" # coverage constraints\n",
" for i in r_cli:\n",
" model.Add(model.Sum([aij[i,j] * fac_vars[j]\\\n",
" for j in r_fac]) >= 1)\n",
" \n",
" # objective function\n",
" model.Minimize(model.Sum([fac_vars[j] for j in r_fac]))\n",
" \n",
" # solve\n",
" solve = model.Solve()\n",
" for j in r_fac:\n",
" if fac_vars[j].solution_value() > 0:\n",
" print 'Facility', fac_vars[j], 'open'\n",
" print \"\"\n",
" print 'Best objective = ', model.Objective().Value()\n",
" print \"Time = \", model.WallTime(), \"milliseconds\"\n",
" \n",
" # linear programming formulation\n",
" lp_formulation = model.ExportModelAsLpFormat(True)\n",
" print \"\"\n",
" print lp_formulation"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Facility y[1] open\n",
"Facility y[2] open\n",
"\n",
"Best objective = 2.0\n",
"Time = 13 milliseconds\n",
"\n",
"\\ Generated by MPModelProtoExporter\n",
"\\ Name : lscp\n",
"\\ Format : Free\n",
"\\ Constraints : 7\n",
"\\ Variables : 5\n",
"\\ Binary : 5\n",
"\\ Integer : 0\n",
"\\ Continuous : 0\n",
"Minimize\n",
" Obj: +1.000000 V0 +1.000000 V1 +1.000000 V2 +1.000000 V3 +1.000000 V4 \n",
"Subject to\n",
" C0: +1.000000 V0 +1.000000 V1 +1.000000 V2 +1.000000 V3 +1.000000 V4 >= 1.000000\n",
" C1: +1.000000 V1 +1.000000 V2 >= 1.000000\n",
" C2: +1.000000 V1 +1.000000 V2 +1.000000 V3 +1.000000 V4 >= 1.000000\n",
" C3: +1.000000 V0 +1.000000 V1 +1.000000 V2 +1.000000 V4 >= 1.000000\n",
" C4: +1.000000 V0 +1.000000 V2 +1.000000 V3 +1.000000 V4 >= 1.000000\n",
" C5: +1.000000 V0 +1.000000 V2 +1.000000 V3 +1.000000 V4 >= 1.000000\n",
" C6: +1.000000 V0 +1.000000 V1 +1.000000 V3 >= 1.000000\n",
"Bounds\n",
" 0 <= V0 <= 1\n",
" 0 <= V1 <= 1\n",
" 0 <= V2 <= 1\n",
" 0 <= V3 <= 1\n",
" 0 <= V4 <= 1\n",
"Binaries\n",
" V0\n",
" V1\n",
" V2\n",
" V3\n",
" V4\n",
"End\n",
"\n"
]
}
],
"source": [
"coverage = 3.\n",
"lscp(\"lscp\", cij=cost_matrix, s=coverage)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"--------------------"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## *p*-median Problem"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def pmp(name, ai=None, cij=None, p=None):\n",
" \"\"\"\n",
" p-median Problem\n",
" \n",
" Originally Published:\n",
" S. L. Hakimi. 1964. Optimum Locations of Switching Centers and \n",
" the Absolute Centers and Medians of a Graph. Operations Research.\n",
" 12 (3):450-459.\n",
" Adapted from:\n",
" -1-\n",
" ReVelle, C.S. and Swain, R.W. 1970. Central facilities location.\n",
" Geographical Analysis. 2(1), 30-42.\n",
" -2-\n",
" Toregas, C., Swain, R., ReVelle, C., Bergman, L. 1971. The Location\n",
" of Emergency Service Facilities. Operations Research. 19 (6),\n",
" 1363-1373.\n",
" - 3 -\n",
" Daskin, M. (1995). Network and discrete location: Models, algorithms,\n",
" and applications. New York: John Wiley and Sons, Inc.\n",
" \n",
" Parameters\n",
" ----------\n",
" name : str\n",
" model name\n",
" ai : numpy.ndarray\n",
" client weight vector\n",
" cij : numpy.ndarray\n",
" cost matrix\n",
" p : int\n",
" density of facilities to site\n",
" \"\"\"\n",
" # client count and range\n",
" n_cli, r_cli = cij.shape[0], xrange(cij.shape[0])\n",
" \n",
" # facility count and range\n",
" n_fac, r_fac = cij.shape[1], xrange(cij.shape[1])\n",
" \n",
" # weighted demand\n",
" sij = ai * cij\n",
" \n",
" # create a solver instance\n",
" solver_instance = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING\n",
" \n",
" # instantiate a model\n",
" model = pywraplp.Solver(name, solver_instance)\n",
" \n",
" # client decision variables\n",
" cli_vars = {(i,j): model.IntVar(0, 1, 'x[%i,%i]' % (i,j))\\\n",
" for i in r_cli for j in r_fac}\n",
" # facility decision variables\n",
" fac_vars = {j: model.IntVar(0, 1, 'y[%i]' % (j)) for j in r_fac}\n",
" \n",
" # assignment constraints\n",
" for i in r_cli:\n",
" model.Add(model.Sum([cli_vars[i,j] for j in r_fac]) == 1)\n",
" \n",
" # facilty constraint\n",
" model.Add(model.Sum([fac_vars[j] for j in r_fac]) == p)\n",
"\n",
" # opening constraints\n",
" for i in r_cli:\n",
" for j in r_fac:\n",
" model.Add(fac_vars[j] - cli_vars[i,j] >= 0)\n",
" \n",
" # objective function\n",
" obj = [sij[i,j] * cli_vars[i,j]\\\n",
" for i in r_cli for j in r_fac]\n",
" model.Minimize(model.Sum(obj))\n",
" \n",
" # solve\n",
" solve = model.Solve()\n",
" for j in r_fac:\n",
" if fac_vars[j].solution_value() > 0:\n",
" print 'Facility', fac_vars[j], 'open'\n",
" for i in r_cli:\n",
" if cli_vars[i,j].solution_value() > 0: \n",
" print '\\tServing', cli_vars[i,j]\n",
" print \"\"\n",
" print 'Best objective = ', model.Objective().Value()\n",
" print \"Time = \", model.WallTime(), \"milliseconds\"\n",
" \n",
" # linear programming formulation\n",
" lp_formulation = model.ExportModelAsLpFormat(True)\n",
" print \"\"\n",
" print lp_formulation"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Facility y[0] open\n",
"\tServing x[0,0]\n",
"\tServing x[1,0]\n",
"\tServing x[6,0]\n",
"Facility y[4] open\n",
"\tServing x[2,4]\n",
"\tServing x[3,4]\n",
"\tServing x[4,4]\n",
"\tServing x[5,4]\n",
"\n",
"Best objective = 41.07\n",
"Time = 14 milliseconds\n",
"\n",
"\\ Generated by MPModelProtoExporter\n",
"\\ Name : pmp\n",
"\\ Format : Free\n",
"\\ Constraints : 43\n",
"\\ Variables : 40\n",
"\\ Binary : 40\n",
"\\ Integer : 0\n",
"\\ Continuous : 0\n",
"Minimize\n",
" Obj: +1.420000 V00 +1.770000 V01 +2.610000 V02 +1.560000 V03 +2.790000 V04 +12.560000 V05 +4.920000 V06 +11.840000 V07 +13.520000 V08 +12.840000 V09 +26.240000 V10 +15.440000 V11 +22.720000 V12 +9.280000 V13 +10.640000 V14 +12.240000 V15 +8.960000 V16 +1.280000 V17 +31.120000 V18 +6.400000 V19 +4.900000 V20 +18.700000 V21 +8.950000 V22 +10.700000 V23 +4.500000 V24 +6.270000 V25 +10.920000 V26 +7.110000 V27 +1.050000 V28 +3.090000 V29 +2.460000 V30 +12.780000 V31 +19.320000 V32 +12.180000 V33 +20.040000 V34 \n",
"Subject to\n",
" C00: +1.000000 V00 +1.000000 V01 +1.000000 V02 +1.000000 V03 +1.000000 V04 = 1.000000\n",
" C01: +1.000000 V05 +1.000000 V06 +1.000000 V07 +1.000000 V08 +1.000000 V09 = 1.000000\n",
" C02: +1.000000 V10 +1.000000 V11 +1.000000 V12 +1.000000 V13 +1.000000 V14 = 1.000000\n",
" C03: +1.000000 V15 +1.000000 V16 +1.000000 V17 +1.000000 V18 +1.000000 V19 = 1.000000\n",
" C04: +1.000000 V20 +1.000000 V21 +1.000000 V22 +1.000000 V23 +1.000000 V24 = 1.000000\n",
" C05: +1.000000 V25 +1.000000 V26 +1.000000 V27 +1.000000 V28 +1.000000 V29 = 1.000000\n",
" C06: +1.000000 V30 +1.000000 V31 +1.000000 V32 +1.000000 V33 +1.000000 V34 = 1.000000\n",
" C07: +1.000000 V35 +1.000000 V36 +1.000000 V37 +1.000000 V38 +1.000000 V39 = 2.000000\n",
" C08: -1.000000 V00 +1.000000 V35 >= 0.000000\n",
" C09: -1.000000 V01 +1.000000 V36 >= 0.000000\n",
" C10: -1.000000 V02 +1.000000 V37 >= 0.000000\n",
" C11: -1.000000 V03 +1.000000 V38 >= 0.000000\n",
" C12: -1.000000 V04 +1.000000 V39 >= 0.000000\n",
" C13: -1.000000 V05 +1.000000 V35 >= 0.000000\n",
" C14: -1.000000 V06 +1.000000 V36 >= 0.000000\n",
" C15: -1.000000 V07 +1.000000 V37 >= 0.000000\n",
" C16: -1.000000 V08 +1.000000 V38 >= 0.000000\n",
" C17: -1.000000 V09 +1.000000 V39 >= 0.000000\n",
" C18: -1.000000 V10 +1.000000 V35 >= 0.000000\n",
" C19: -1.000000 V11 +1.000000 V36 >= 0.000000\n",
" C20: -1.000000 V12 +1.000000 V37 >= 0.000000\n",
" C21: -1.000000 V13 +1.000000 V38 >= 0.000000\n",
" C22: -1.000000 V14 +1.000000 V39 >= 0.000000\n",
" C23: -1.000000 V15 +1.000000 V35 >= 0.000000\n",
" C24: -1.000000 V16 +1.000000 V36 >= 0.000000\n",
" C25: -1.000000 V17 +1.000000 V37 >= 0.000000\n",
" C26: -1.000000 V18 +1.000000 V38 >= 0.000000\n",
" C27: -1.000000 V19 +1.000000 V39 >= 0.000000\n",
" C28: -1.000000 V20 +1.000000 V35 >= 0.000000\n",
" C29: -1.000000 V21 +1.000000 V36 >= 0.000000\n",
" C30: -1.000000 V22 +1.000000 V37 >= 0.000000\n",
" C31: -1.000000 V23 +1.000000 V38 >= 0.000000\n",
" C32: -1.000000 V24 +1.000000 V39 >= 0.000000\n",
" C33: -1.000000 V25 +1.000000 V35 >= 0.000000\n",
" C34: -1.000000 V26 +1.000000 V36 >= 0.000000\n",
" C35: -1.000000 V27 +1.000000 V37 >= 0.000000\n",
" C36: -1.000000 V28 +1.000000 V38 >= 0.000000\n",
" C37: -1.000000 V29 +1.000000 V39 >= 0.000000\n",
" C38: -1.000000 V30 +1.000000 V35 >= 0.000000\n",
" C39: -1.000000 V31 +1.000000 V36 >= 0.000000\n",
" C40: -1.000000 V32 +1.000000 V37 >= 0.000000\n",
" C41: -1.000000 V33 +1.000000 V38 >= 0.000000\n",
" C42: -1.000000 V34 +1.000000 V39 >= 0.000000\n",
"Bounds\n",
" 0 <= V00 <= 1\n",
" 0 <= V01 <= 1\n",
" 0 <= V02 <= 1\n",
" 0 <= V03 <= 1\n",
" 0 <= V04 <= 1\n",
" 0 <= V05 <= 1\n",
" 0 <= V06 <= 1\n",
" 0 <= V07 <= 1\n",
" 0 <= V08 <= 1\n",
" 0 <= V09 <= 1\n",
" 0 <= V10 <= 1\n",
" 0 <= V11 <= 1\n",
" 0 <= V12 <= 1\n",
" 0 <= V13 <= 1\n",
" 0 <= V14 <= 1\n",
" 0 <= V15 <= 1\n",
" 0 <= V16 <= 1\n",
" 0 <= V17 <= 1\n",
" 0 <= V18 <= 1\n",
" 0 <= V19 <= 1\n",
" 0 <= V20 <= 1\n",
" 0 <= V21 <= 1\n",
" 0 <= V22 <= 1\n",
" 0 <= V23 <= 1\n",
" 0 <= V24 <= 1\n",
" 0 <= V25 <= 1\n",
" 0 <= V26 <= 1\n",
" 0 <= V27 <= 1\n",
" 0 <= V28 <= 1\n",
" 0 <= V29 <= 1\n",
" 0 <= V30 <= 1\n",
" 0 <= V31 <= 1\n",
" 0 <= V32 <= 1\n",
" 0 <= V33 <= 1\n",
" 0 <= V34 <= 1\n",
" 0 <= V35 <= 1\n",
" 0 <= V36 <= 1\n",
" 0 <= V37 <= 1\n",
" 0 <= V38 <= 1\n",
" 0 <= V39 <= 1\n",
"Binaries\n",
" V00\n",
" V01\n",
" V02\n",
" V03\n",
" V04\n",
" V05\n",
" V06\n",
" V07\n",
" V08\n",
" V09\n",
" V10\n",
" V11\n",
" V12\n",
" V13\n",
" V14\n",
" V15\n",
" V16\n",
" V17\n",
" V18\n",
" V19\n",
" V20\n",
" V21\n",
" V22\n",
" V23\n",
" V24\n",
" V25\n",
" V26\n",
" V27\n",
" V28\n",
" V29\n",
" V30\n",
" V31\n",
" V32\n",
" V33\n",
" V34\n",
" V35\n",
" V36\n",
" V37\n",
" V38\n",
" V39\n",
"End\n",
"\n"
]
}
],
"source": [
"pmp(\"pmp\", ai=weights, cij=cost_matrix, p=facilities)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"--------------------------------------------"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## *p*-center Problem"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def pcp(name, cij=None, p=None):\n",
" \"\"\"\n",
" p-center Problem\n",
" \n",
" Originally Published:\n",
" S. L. Hakimi. 1964. Optimum Locations of Switching Centers and \n",
" the Absolute Centers and Medians of a Graph. Operations Research. \n",
" 12 (3):450-459.\n",
" Adapted from:\n",
" Daskin, M. (1995). Network and discrete location: Models, algorithms,\n",
" and applications. New York: John Wiley and Sons, Inc.\n",
" \n",
" Parameters\n",
" ----------\n",
" name : str\n",
" model name\n",
" cij : numpy.ndarray\n",
" cost matrix\n",
" p : int\n",
" density of facilities to site\n",
" \"\"\"\n",
" # client count and range\n",
" n_cli, r_cli = cij.shape[0], xrange(cij.shape[0])\n",
" \n",
" # facility count and range\n",
" n_fac, r_fac = cij.shape[1], xrange(cij.shape[1])\n",
" \n",
" # create a solver instance\n",
" solver_instance = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING\n",
" \n",
" # instantiate a model\n",
" model = pywraplp.Solver(name, solver_instance)\n",
" \n",
" # client decision variables\n",
" cli_vars = {(i,j): model.IntVar(0, 1, 'x[%i,%i]' % (i,j))\\\n",
" for i in r_cli for j in r_fac}\n",
" # facility decision variables\n",
" fac_vars = {j: model.IntVar(0, 1, 'y[%i]' % (j)) for j in r_fac}\n",
" \n",
" # minimize maximum variable\n",
" W = model.NumVar(0, model.infinity(), 'W')\n",
" \n",
" # assignment constraints\n",
" for i in r_cli:\n",
" model.Add(model.Sum([cli_vars[i,j] for j in r_fac]) == 1)\n",
" \n",
" # facilty constraint\n",
" model.Add(model.Sum([fac_vars[j] for j in r_fac]) == p)\n",
"\n",
" # opening constraints\n",
" for i in r_cli:\n",
" for j in r_fac:\n",
" model.Add(fac_vars[j] - cli_vars[i,j] >= 0)\n",
" \n",
" # minimize maximum constraints\n",
" for j in r_fac:\n",
" model.Add(model.Sum([cij[i,j] * cli_vars[i,j] for i in r_cli]) <= W)\n",
" \n",
" # objective function\n",
" model.Minimize(W)\n",
" \n",
" # solve\n",
" solve = model.Solve()\n",
" print \"\"\n",
" for j in r_fac:\n",
" if fac_vars[j].solution_value() > 0:\n",
" print 'Facility', fac_vars[j], 'open'\n",
" for i in r_cli:\n",
" if cli_vars[i,j].solution_value() > 0: \n",
" print '\\tServing', cli_vars[i,j]\n",
" print \"\"\n",
" print 'Best objective = ', model.Objective().Value()\n",
" print \"Time = \", model.WallTime(), \"milliseconds\"\n",
" \n",
" # linear programming formulation\n",
" lp_formulation = model.ExportModelAsLpFormat(True)\n",
" print \"\"\n",
" print lp_formulation"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Facility y[0] open\n",
"\tServing x[0,0]\n",
"\tServing x[3,0]\n",
"\tServing x[4,0]\n",
"\tServing x[6,0]\n",
"Facility y[3] open\n",
"\tServing x[1,3]\n",
"\tServing x[2,3]\n",
"\tServing x[5,3]\n",
"\n",
"Best objective = 4.89\n",
"Time = 587 milliseconds\n",
"\n",
"\\ Generated by MPModelProtoExporter\n",
"\\ Name : pcp\n",
"\\ Format : Free\n",
"\\ Constraints : 48\n",
"\\ Variables : 41\n",
"\\ Binary : 40\n",
"\\ Integer : 0\n",
"\\ Continuous : 1\n",
"Minimize\n",
" Obj: +1.000000 V40 \n",
"Subject to\n",
" C00: +1.000000 V00 +1.000000 V01 +1.000000 V02 +1.000000 V03 +1.000000 V04 = 1.000000\n",
" C01: +1.000000 V05 +1.000000 V06 +1.000000 V07 +1.000000 V08 +1.000000 V09 = 1.000000\n",
" C02: +1.000000 V10 +1.000000 V11 +1.000000 V12 +1.000000 V13 +1.000000 V14 = 1.000000\n",
" C03: +1.000000 V15 +1.000000 V16 +1.000000 V17 +1.000000 V18 +1.000000 V19 = 1.000000\n",
" C04: +1.000000 V20 +1.000000 V21 +1.000000 V22 +1.000000 V23 +1.000000 V24 = 1.000000\n",
" C05: +1.000000 V25 +1.000000 V26 +1.000000 V27 +1.000000 V28 +1.000000 V29 = 1.000000\n",
" C06: +1.000000 V30 +1.000000 V31 +1.000000 V32 +1.000000 V33 +1.000000 V34 = 1.000000\n",
" C07: +1.000000 V35 +1.000000 V36 +1.000000 V37 +1.000000 V38 +1.000000 V39 = 2.000000\n",
" C08: -1.000000 V00 +1.000000 V35 >= 0.000000\n",
" C09: -1.000000 V01 +1.000000 V36 >= 0.000000\n",
" C10: -1.000000 V02 +1.000000 V37 >= 0.000000\n",
" C11: -1.000000 V03 +1.000000 V38 >= 0.000000\n",
" C12: -1.000000 V04 +1.000000 V39 >= 0.000000\n",
" C13: -1.000000 V05 +1.000000 V35 >= 0.000000\n",
" C14: -1.000000 V06 +1.000000 V36 >= 0.000000\n",
" C15: -1.000000 V07 +1.000000 V37 >= 0.000000\n",
" C16: -1.000000 V08 +1.000000 V38 >= 0.000000\n",
" C17: -1.000000 V09 +1.000000 V39 >= 0.000000\n",
" C18: -1.000000 V10 +1.000000 V35 >= 0.000000\n",
" C19: -1.000000 V11 +1.000000 V36 >= 0.000000\n",
" C20: -1.000000 V12 +1.000000 V37 >= 0.000000\n",
" C21: -1.000000 V13 +1.000000 V38 >= 0.000000\n",
" C22: -1.000000 V14 +1.000000 V39 >= 0.000000\n",
" C23: -1.000000 V15 +1.000000 V35 >= 0.000000\n",
" C24: -1.000000 V16 +1.000000 V36 >= 0.000000\n",
" C25: -1.000000 V17 +1.000000 V37 >= 0.000000\n",
" C26: -1.000000 V18 +1.000000 V38 >= 0.000000\n",
" C27: -1.000000 V19 +1.000000 V39 >= 0.000000\n",
" C28: -1.000000 V20 +1.000000 V35 >= 0.000000\n",
" C29: -1.000000 V21 +1.000000 V36 >= 0.000000\n",
" C30: -1.000000 V22 +1.000000 V37 >= 0.000000\n",
" C31: -1.000000 V23 +1.000000 V38 >= 0.000000\n",
" C32: -1.000000 V24 +1.000000 V39 >= 0.000000\n",
" C33: -1.000000 V25 +1.000000 V35 >= 0.000000\n",
" C34: -1.000000 V26 +1.000000 V36 >= 0.000000\n",
" C35: -1.000000 V27 +1.000000 V37 >= 0.000000\n",
" C36: -1.000000 V28 +1.000000 V38 >= 0.000000\n",
" C37: -1.000000 V29 +1.000000 V39 >= 0.000000\n",
" C38: -1.000000 V30 +1.000000 V35 >= 0.000000\n",
" C39: -1.000000 V31 +1.000000 V36 >= 0.000000\n",
" C40: -1.000000 V32 +1.000000 V37 >= 0.000000\n",
" C41: -1.000000 V33 +1.000000 V38 >= 0.000000\n",
" C42: -1.000000 V34 +1.000000 V39 >= 0.000000\n",
" C43: +1.420000 V00 +3.140000 V05 +3.280000 V10 +1.530000 V15 +0.980000 V20 +2.090000 V25 +0.410000 V30 -1.000000 V40 <= 0.000000\n",
" C44: +1.770000 V01 +1.230000 V06 +1.930000 V11 +1.120000 V16 +3.740000 V21 +3.640000 V26 +2.130000 V31 -1.000000 V40 <= 0.000000\n",
" C45: +2.610000 V02 +2.960000 V07 +2.840000 V12 +0.160000 V17 +1.790000 V22 +2.370000 V27 +3.220000 V32 -1.000000 V40 <= 0.000000\n",
" C46: +1.560000 V03 +3.380000 V08 +1.160000 V13 +3.890000 V18 +2.140000 V23 +0.350000 V28 +2.030000 V33 -1.000000 V40 <= 0.000000\n",
" C47: +2.790000 V04 +3.210000 V09 +1.330000 V14 +0.800000 V19 +0.900000 V24 +1.030000 V29 +3.340000 V34 -1.000000 V40 <= 0.000000\n",
"Bounds\n",
" 0 <= V00 <= 1\n",
" 0 <= V01 <= 1\n",
" 0 <= V02 <= 1\n",
" 0 <= V03 <= 1\n",
" 0 <= V04 <= 1\n",
" 0 <= V05 <= 1\n",
" 0 <= V06 <= 1\n",
" 0 <= V07 <= 1\n",
" 0 <= V08 <= 1\n",
" 0 <= V09 <= 1\n",
" 0 <= V10 <= 1\n",
" 0 <= V11 <= 1\n",
" 0 <= V12 <= 1\n",
" 0 <= V13 <= 1\n",
" 0 <= V14 <= 1\n",
" 0 <= V15 <= 1\n",
" 0 <= V16 <= 1\n",
" 0 <= V17 <= 1\n",
" 0 <= V18 <= 1\n",
" 0 <= V19 <= 1\n",
" 0 <= V20 <= 1\n",
" 0 <= V21 <= 1\n",
" 0 <= V22 <= 1\n",
" 0 <= V23 <= 1\n",
" 0 <= V24 <= 1\n",
" 0 <= V25 <= 1\n",
" 0 <= V26 <= 1\n",
" 0 <= V27 <= 1\n",
" 0 <= V28 <= 1\n",
" 0 <= V29 <= 1\n",
" 0 <= V30 <= 1\n",
" 0 <= V31 <= 1\n",
" 0 <= V32 <= 1\n",
" 0 <= V33 <= 1\n",
" 0 <= V34 <= 1\n",
" 0 <= V35 <= 1\n",
" 0 <= V36 <= 1\n",
" 0 <= V37 <= 1\n",
" 0 <= V38 <= 1\n",
" 0 <= V39 <= 1\n",
" 0.000000 <= V40\n",
"Binaries\n",
" V00\n",
" V01\n",
" V02\n",
" V03\n",
" V04\n",
" V05\n",
" V06\n",
" V07\n",
" V08\n",
" V09\n",
" V10\n",
" V11\n",
" V12\n",
" V13\n",
" V14\n",
" V15\n",
" V16\n",
" V17\n",
" V18\n",
" V19\n",
" V20\n",
" V21\n",
" V22\n",
" V23\n",
" V24\n",
" V25\n",
" V26\n",
" V27\n",
" V28\n",
" V29\n",
" V30\n",
" V31\n",
" V32\n",
" V33\n",
" V34\n",
" V35\n",
" V36\n",
" V37\n",
" V38\n",
" V39\n",
"End\n",
"\n"
]
}
],
"source": [
"pcp(\"pcp\", cij=cost_matrix, p=facilities)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"--------------------------------------"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Maximal Covering Location Problem"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def mclp(name, ai=None, cij=None, s=None, p=None):\n",
" \"\"\"\n",
" Originally Published:\n",
" Church, R. L and C. ReVelle. 1974. The Maximal Covering Location\n",
" Problem. Papers of the Regional Science Association. 32:101-18.\n",
" \n",
" Parameters\n",
" ----------\n",
" name : str\n",
" model name\n",
" ai : numpy.ndarray\n",
" client weight vector\n",
" cij : numpy.ndarray\n",
" cost matrix\n",
" s : float\n",
" service radius\n",
" p : int\n",
" density of facilities to site\n",
" \"\"\"\n",
" # client count and range\n",
" n_cli, r_cli = cij.shape[0], xrange(cij.shape[0])\n",
" \n",
" # facility count and range\n",
" n_fac, r_fac = cij.shape[1], xrange(cij.shape[1])\n",
" \n",
" # binary coverage matrix from cij\n",
" aij = np.zeros(cij.shape)\n",
" aij[cij <= s] = 1.\n",
" \n",
" # create a solver instance\n",
" solver_instance = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING\n",
" # instantiate a model\n",
" model = pywraplp.Solver(name, solver_instance)\n",
" \n",
" # client decision variables\n",
" cli_vars = {i: model.IntVar(0, 1, 'x[%i]' % (i)) for i in r_cli}\n",
" \n",
" # facility decision variables\n",
" fac_vars = {j: model.IntVar(0, 1, 'y[%i]' % (j)) for j in r_fac}\n",
" \n",
" # facilty constraint\n",
" model.Add(model.Sum([fac_vars[j] for j in r_fac]) == p)\n",
" \n",
" # maximal coverage constraints\n",
" for i in r_cli:\n",
" model.Add(model.Sum([aij[i,j] * fac_vars[j]\\\n",
" for j in r_fac]) >= cli_vars[i])\n",
"\n",
" # objective\n",
" obj = [ai.flatten()[i] * cli_vars[i] for i in r_cli]\n",
" model.Maximize(model.Sum(obj))\n",
" \n",
" # solve\n",
" solve = model.Solve()\n",
" print \"\"\n",
" for j in r_fac:\n",
" if fac_vars[j].solution_value() > 0:\n",
" print 'Facility', fac_vars[j], 'open'\n",
" print \"\"\n",
" print 'Best objective = ', model.Objective().Value()\n",
" print \"Time = \", model.WallTime(), \"milliseconds\"\n",
" \n",
" # linear programming formulation\n",
" lp_formulation = model.ExportModelAsLpFormat(True)\n",
" print \"\"\n",
" print lp_formulation "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Facility y[0] open\n",
"Facility y[4] open\n",
"\n",
"Best objective = 31.0\n",
"Time = 48 milliseconds\n",
"\n",
"\\ Generated by MPModelProtoExporter\n",
"\\ Name : mclp\n",
"\\ Format : Free\n",
"\\ Constraints : 8\n",
"\\ Variables : 12\n",
"\\ Binary : 12\n",
"\\ Integer : 0\n",
"\\ Continuous : 0\n",
"Maximize\n",
" Obj: +1.000000 V00 +4.000000 V01 +8.000000 V02 +8.000000 V03 +5.000000 V04 +3.000000 V05 +6.000000 V06 \n",
"Subject to\n",
" C0: +1.000000 V07 +1.000000 V08 +1.000000 V09 +1.000000 V10 +1.000000 V11 = 2.000000\n",
" C1: -1.000000 V00 +1.000000 V07 >= 0.000000\n",
" C2: -1.000000 V01 +1.000000 V08 >= 0.000000\n",
" C3: -1.000000 V02 +1.000000 V10 +1.000000 V11 >= 0.000000\n",
" C4: -1.000000 V03 +1.000000 V08 +1.000000 V09 +1.000000 V11 >= 0.000000\n",
" C5: -1.000000 V04 +1.000000 V07 +1.000000 V11 >= 0.000000\n",
" C6: -1.000000 V05 +1.000000 V10 +1.000000 V11 >= 0.000000\n",
" C7: -1.000000 V06 +1.000000 V07 >= 0.000000\n",
"Bounds\n",
" 0 <= V00 <= 1\n",
" 0 <= V01 <= 1\n",
" 0 <= V02 <= 1\n",
" 0 <= V03 <= 1\n",
" 0 <= V04 <= 1\n",
" 0 <= V05 <= 1\n",
" 0 <= V06 <= 1\n",
" 0 <= V07 <= 1\n",
" 0 <= V08 <= 1\n",
" 0 <= V09 <= 1\n",
" 0 <= V10 <= 1\n",
" 0 <= V11 <= 1\n",
"Binaries\n",
" V00\n",
" V01\n",
" V02\n",
" V03\n",
" V04\n",
" V05\n",
" V06\n",
" V07\n",
" V08\n",
" V09\n",
" V10\n",
" V11\n",
"End\n",
"\n"
]
}
],
"source": [
"coverage = 1.5\n",
"mclp(\"mclp\", ai=weights, cij=cost_matrix, s=coverage, p=facilities)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"-------------------------------------"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:py2_at1866]",
"language": "python",
"name": "conda-env-py2_at1866-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment