optimisation-101-with-or-tools.ipynb

{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Optimisation 101 with OR-Tools\n\n_By Nick Brooks, February 2020_\n\n**Resource:** <br>\n- https://developers.google.com/optimization/lp/lp\n- https://github.com/google/or-tools/tree/stable/examples/python\n- https://docs.mosek.com/modeling-cookbook/intro.html\n\nBig Shoutout to Google Developer Guides for putting together top notch guides! This is not my work.\n\n## Table of Content\n**Types of Optimization:**\n1. Linear Optimisation\n1. Contrained Opimisation\n1. Integer Optimisation\n1. Routing Problems\n\n**Sample Problems:**\n1. Stigler Diet Problem\n1. Sodoku Problem\n1. Clustering Problem\n1. Zebra Problem (Lewis Caroll)\n\n\n## What is Optimization?\nIt is the process of finding the optimal point on a function. This is what fuels most machine learning techniques. However, when people refer to just \"Optimization\", it is often assumed that the problem has contraints, hence **Contrained Optimization**. These introduce a lot of complexity, especially as the problem introduces more dimensions.\n\n\n## Linear Optimisation\nProcess of optimizing a linear fucntion subject to linear contraints. The [Simplex Algorithmn](https://en.wikipedia.org/wiki/Simplex_algorithm) is the most famous of this group. This algorithmn finds an optimal solution by traversing the vertices of the feasible region until there is no more improvement in the objective function\n\nIn this notebook, I will be following the Google Developer Guide for OR-Tools, where I will be adding modifications and thoughts along the way.\n\n\n## Contrainted Optimization\n\n*Contraintes Optimization* is based on feasibility (finding a feasible solution) rather than optimization (finding an optimal solution) and focuses on the constraints and variables rather than the objective function. In fact, a CP problem may not even have an objective function — the goal may simply be to narrow down a vary large set of possible solutions to a more manageable subset by adding constraints to the problem.\n\n\n## Integer Optimization\n\nProblems that require some of the variables to be integers.\n\n\n## Routing Problems\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport time\nimport matplotlib.pyplot as plt\nimport pprint\nfrom collections import Counter\nfrom IPython.display import Image\nfrom IPython.core.display import HTML \n\npd.options.display.max_rows = 999\npd.options.display.width = 300\npd.options.display.max_columns = 500\nstart = time.time()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Problems\n\n## Stigler Diet Problem\n\n[Link to Github Source](https://github.com/google/or-tools/blob/stable/examples/python/stigler_diet.py)\n\nA Resource Allocation Problem. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Define Table..\")\n# Nutrient minimums.\nnutrients = [['Calories (kcal)', 3], ['Protein (g)', 70], [\n    'Calcium (g)', 0.8\n], ['Iron (mg)', 12], ['Vitamin A (KIU)', 5], ['Vitamin B1 (mg)', 1.8],\n             ['Vitamin B2 (mg)', 2.7], ['Niacin (mg)',\n                                        18], ['Vitamin C (mg)', 75]]\n\n# Commodity, Unit, 1939 price (cents), Calories (kcal), Protein (g), Calcium (g), Iron (mg),\n# Vitamin A (KIU), Vitamin B1 (mg), Vitamin B2 (mg), Niacin (mg), Vitamin C (mg)\ndata = [[\n    'Wheat Flour (Enriched)', '10 lb.', 36, 44.7, 1411, 2, 365, 0, 55.4,\n    33.3, 441, 0\n], ['Macaroni', '1 lb.', 14.1, 11.6, 418, 0.7, 54, 0, 3.2, 1.9, 68, 0], [\n    'Wheat Cereal (Enriched)', '28 oz.', 24.2, 11.8, 377, 14.4, 175, 0,\n    14.4, 8.8, 114, 0\n], ['Corn Flakes', '8 oz.', 7.1, 11.4, 252, 0.1, 56, 0, 13.5, 2.3, 68, 0], [\n    'Corn Meal', '1 lb.', 4.6, 36.0, 897, 1.7, 99, 30.9, 17.4, 7.9, 106, 0\n], [\n    'Hominy Grits', '24 oz.', 8.5, 28.6, 680, 0.8, 80, 0, 10.6, 1.6, 110, 0\n], ['Rice', '1 lb.', 7.5, 21.2, 460, 0.6, 41, 0, 2, 4.8, 60, 0], [\n    'Rolled Oats', '1 lb.', 7.1, 25.3, 907, 5.1, 341, 0, 37.1, 8.9, 64, 0\n], [\n    'White Bread (Enriched)', '1 lb.', 7.9, 15.0, 488, 2.5, 115, 0, 13.8,\n    8.5, 126, 0\n], [\n    'Whole Wheat Bread', '1 lb.', 9.1, 12.2, 484, 2.7, 125, 0, 13.9, 6.4,\n    160, 0\n], ['Rye Bread', '1 lb.', 9.1, 12.4, 439, 1.1, 82, 0, 9.9, 3, 66, 0], [\n    'Pound Cake', '1 lb.', 24.8, 8.0, 130, 0.4, 31, 18.9, 2.8, 3, 17, 0\n], ['Soda Crackers', '1 lb.', 15.1, 12.5, 288, 0.5, 50, 0, 0, 0, 0, 0], [\n    'Milk', '1 qt.', 11, 6.1, 310, 10.5, 18, 16.8, 4, 16, 7, 177\n], [\n    'Evaporated Milk (can)', '14.5 oz.', 6.7, 8.4, 422, 15.1, 9, 26, 3,\n    23.5, 11, 60\n], ['Butter', '1 lb.', 30.8, 10.8, 9, 0.2, 3, 44.2, 0, 0.2, 2, 0], [\n    'Oleomargarine', '1 lb.', 16.1, 20.6, 17, 0.6, 6, 55.8, 0.2, 0, 0, 0\n], ['Eggs', '1 doz.', 32.6, 2.9, 238, 1.0, 52, 18.6, 2.8, 6.5, 1, 0], [\n    'Cheese (Cheddar)', '1 lb.', 24.2, 7.4, 448, 16.4, 19, 28.1, 0.8, 10.3,\n    4, 0\n], ['Cream', '1/2 pt.', 14.1, 3.5, 49, 1.7, 3, 16.9, 0.6, 2.5, 0, 17], [\n    'Peanut Butter', '1 lb.', 17.9, 15.7, 661, 1.0, 48, 0, 9.6, 8.1, 471, 0\n], ['Mayonnaise', '1/2 pt.', 16.7, 8.6, 18, 0.2, 8, 2.7, 0.4, 0.5, 0, 0], [\n    'Crisco', '1 lb.', 20.3, 20.1, 0, 0, 0, 0, 0, 0, 0, 0\n], ['Lard', '1 lb.', 9.8, 41.7, 0, 0, 0, 0.2, 0, 0.5, 5, 0], [\n    'Sirloin Steak', '1 lb.', 39.6, 2.9, 166, 0.1, 34, 0.2, 2.1, 2.9, 69, 0\n], ['Round Steak', '1 lb.', 36.4, 2.2, 214, 0.1, 32, 0.4, 2.5, 2.4, 87, 0\n   ], ['Rib Roast', '1 lb.', 29.2, 3.4, 213, 0.1, 33, 0, 0, 2, 0, 0], [\n       'Chuck Roast', '1 lb.', 22.6, 3.6, 309, 0.2, 46, 0.4, 1, 4, 120, 0\n   ], ['Plate', '1 lb.', 14.6, 8.5, 404, 0.2, 62, 0, 0.9, 0, 0, 0], [\n       'Liver (Beef)', '1 lb.', 26.8, 2.2, 333, 0.2, 139, 169.2, 6.4, 50.8,\n       316, 525\n   ], [\n       'Leg of Lamb', '1 lb.', 27.6, 3.1, 245, 0.1, 20, 0, 2.8, 3.9, 86, 0\n   ], [\n       'Lamb Chops (Rib)',\n       '1 lb.', 36.6, 3.3, 140, 0.1, 15, 0, 1.7, 2.7, 54, 0\n   ], [\n       'Pork Chops', '1 lb.', 30.7, 3.5, 196, 0.2, 30, 0, 17.4, 2.7, 60, 0\n   ], [\n       'Pork Loin Roast',\n       '1 lb.', 24.2, 4.4, 249, 0.3, 37, 0, 18.2, 3.6, 79, 0\n   ], ['Bacon', '1 lb.', 25.6, 10.4, 152, 0.2, 23, 0, 1.8, 1.8, 71, 0], [\n       'Ham, smoked', '1 lb.', 27.4, 6.7, 212, 0.2, 31, 0, 9.9, 3.3, 50, 0\n   ], ['Salt Pork', '1 lb.', 16, 18.8, 164, 0.1, 26, 0, 1.4, 1.8, 0, 0], [\n       'Roasting Chicken', '1 lb.', 30.3, 1.8, 184, 0.1, 30, 0.1, 0.9, 1.8,\n       68, 46\n   ], [\n       'Veal Cutlets', '1 lb.', 42.3, 1.7, 156, 0.1, 24, 0, 1.4, 2.4, 57, 0\n   ], [\n       'Salmon, Pink (can)', '16 oz.', 13, 5.8, 705, 6.8, 45, 3.5,\n       1, 4.9, 209, 0\n   ], ['Apples', '1 lb.', 4.4, 5.8, 27, 0.5, 36, 7.3, 3.6, 2.7, 5, 544], [\n       'Bananas', '1 lb.', 6.1, 4.9, 60, 0.4, 30, 17.4, 2.5, 3.5, 28, 498\n   ], ['Lemons', '1 doz.', 26, 1.0, 21, 0.5, 14, 0, 0.5, 0, 4, 952], [\n       'Oranges', '1 doz.', 30.9, 2.2, 40, 1.1, 18, 11.1, 3.6, 1.3, 10, 1998\n   ], [\n       'Green Beans', '1 lb.', 7.1, 2.4, 138, 3.7, 80, 69, 4.3, 5.8, 37, 862\n   ], ['Cabbage', '1 lb.', 3.7, 2.6, 125, 4.0, 36, 7.2, 9, 4.5, 26, 5369], [\n       'Carrots', '1 bunch', 4.7, 2.7, 73, 2.8, 43, 188.5, 6.1, 4.3, 89, 608\n   ], ['Celery', '1 stalk', 7.3, 0.9, 51, 3.0, 23, 0.9, 1.4, 1.4, 9, 313], [\n       'Lettuce', '1 head', 8.2, 0.4, 27, 1.1, 22, 112.4, 1.8, 3.4, 11, 449\n   ], ['Onions', '1 lb.', 3.6, 5.8, 166, 3.8, 59, 16.6, 4.7, 5.9, 21,\n       1184], [\n           'Potatoes', '15 lb.', 34, 14.3, 336, 1.8, 118, 6.7, 29.4, 7.1,\n           198, 2522\n       ], [\n           'Spinach', '1 lb.', 8.1, 1.1, 106, 0, 138, 918.4, 5.7, 13.8, 33,\n           2755\n       ], [\n           'Sweet Potatoes', '1 lb.', 5.1, 9.6, 138, 2.7, 54, 290.7, 8.4,\n           5.4, 83, 1912\n       ], [\n           'Peaches (can)', 'No. 2 1/2', 16.8, 3.7, 20, 0.4, 10, 21.5, 0.5,\n           1, 31, 196\n       ], [\n           'Pears (can)', 'No. 2 1/2', 20.4, 3.0, 8, 0.3, 8, 0.8, 0.8, 0.8,\n           5, 81\n       ], [\n           'Pineapple (can)', 'No. 2 1/2', 21.3, 2.4, 16, 0.4, 8, 2, 2.8,\n           0.8, 7, 399\n       ], [\n           'Asparagus (can)', 'No. 2', 27.7, 0.4, 33, 0.3, 12, 16.3, 1.4,\n           2.1, 17, 272\n       ], [\n           'Green Beans (can)', 'No. 2', 10, 1.0, 54, 2, 65, 53.9, 1.6, 4.3,\n           32, 431\n       ], [\n           'Pork and Beans (can)', '16 oz.', 7.1, 7.5, 364, 4, 134, 3.5,\n           8.3, 7.7, 56, 0\n       ], [\n           'Corn (can)', 'No. 2', 10.4, 5.2, 136, 0.2, 16, 12, 1.6, 2.7, 42,\n           218\n       ], [\n           'Peas (can)', 'No. 2', 13.8, 2.3, 136, 0.6, 45, 34.9, 4.9, 2.5,\n           37, 370\n       ], [\n           'Tomatoes (can)', 'No. 2', 8.6, 1.3, 63, 0.7, 38, 53.2, 3.4, 2.5,\n           36, 1253\n       ], [\n           'Tomato Soup (can)', '10 1/2 oz.', 7.6, 1.6, 71, 0.6, 43, 57.9,\n           3.5, 2.4, 67, 862\n       ], [\n           'Peaches, Dried', '1 lb.', 15.7, 8.5, 87, 1.7, 173, 86.8, 1.2,\n           4.3, 55, 57\n       ], [\n           'Prunes, Dried', '1 lb.', 9, 12.8, 99, 2.5, 154, 85.7, 3.9, 4.3,\n           65, 257\n       ], [\n           'Raisins, Dried', '15 oz.', 9.4, 13.5, 104, 2.5, 136, 4.5, 6.3,\n           1.4, 24, 136\n       ], [\n           'Peas, Dried', '1 lb.', 7.9, 20.0, 1367, 4.2, 345, 2.9, 28.7,\n           18.4, 162, 0\n       ], [\n           'Lima Beans, Dried', '1 lb.', 8.9, 17.4, 1055, 3.7, 459, 5.1,\n           26.9, 38.2, 93, 0\n       ], [\n           'Navy Beans, Dried', '1 lb.', 5.9, 26.9, 1691, 11.4, 792, 0,\n           38.4, 24.6, 217, 0\n       ], ['Coffee', '1 lb.', 22.4, 0, 0, 0, 0, 0, 4, 5.1, 50,\n           0], ['Tea', '1/4 lb.', 17.4, 0, 0, 0, 0, 0, 0, 2.3, 42, 0],\n        ['Cocoa', '8 oz.', 8.6, 8.7, 237, 3, 72, 0, 2, 11.9, 40, 0], [\n            'Chocolate', '8 oz.', 16.2, 8.0, 77, 1.3, 39, 0, 0.9, 3.4, 14, 0\n        ], ['Sugar', '10 lb.', 51.7, 34.9, 0, 0, 0, 0, 0, 0, 0, 0],\n        ['Corn Syrup', '24 oz.', 13.7, 14.7, 0, 0.5, 74, 0, 0, 0, 5, 0], [\n            'Molasses', '18 oz.', 13.6, 9.0, 0, 10.3, 244, 0, 1.9, 7.5, 146,\n            0\n        ], [\n            'Strawberry Preserves', '1 lb.', 20.5, 6.4, 11, 0.4, 7, 0.2,\n            0.2, 0.4, 3, 0\n        ]]\n\nprint(\"Nutritional Contraints Table\")\ndisplay(pd.DataFrame(nutrients, columns = ['Nutrient', 'Daily Recommended Intake']).set_index('Nutrient'))\nprint(\"Food Nutrition Table\")\ndisplay(pd.DataFrame(data,columns = [\n    \"Commodity\",\n    \"Unit\",\n    \"1939 price (cents)\",\n    \"Calories\",\n    \"Protein (g)\",\n    \"Calcium (g)\",\n    \"Iron (mg)\",\n    \"Vitamin A (IU)\",\n    \"Thiamine (mg)\",\n    \"Riboflavin (mg)\",\n    \"Niacin (mg)\",\n    \"Ascorbic Acid (mg)\"]).set_index(\"Commodity\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note on the Food Nutrition Table, it is a bit misleading. While there is the unit and price information, all the other columns are nutrition levels scaled per dollar. This enables the objective function weight to be fixed at 1..\n\n### Framing the Problem\nFraming an optimization problem is a big part of the work. It is here where understanding about the objective and contraints must be materialized mathematically. The nature of these functions will dictate which solution to apply.\n\n#### **Objective Function:** <br>\n\\begin{align}\nf(minimize) & = \\sum_{i} food[i] \\cdot 1\n\\end{align}\n\nwhere Food is the amount of money spent on foodstuff i.\n\n\n#### **Contraints:** <br>\n\\begin{align}\n\\sum_{j} data[j][i+3] \\cdot food[j] \\geq nutrients[i][1]\n\\end{align}\n\nWhere data[j][i+3] is the calories per dollar of food j. and nutrients[i] is the minimal receommended intake value per nutrient i."},{"metadata":{"trusted":true},"cell_type":"code","source":"from six.moves import xrange\nfrom ortools.linear_solver import pywraplp\n\n# Instantiate a Glop solver, naming it LinearExample.\nsolver = pywraplp.Solver('StiglerDietExample',\n                         pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)\n\n# Declare an array to hold our variables.\n# Must be non-negative variables..\nfoods = [solver.NumVar(0.0, solver.infinity(), item[0]) for item in data]\n\n# Objective function: Minimize the sum of (price-normalized) foods.\nobjective = solver.Objective()\n# Here the Objective function looks something like food1 * 1 + food2 * 1 .. \nfor food in foods:\n    objective.SetCoefficient(food, 1)\n# Minimization Problem\nobjective.SetMinimization()\n\n# Create the constraints, one per nutrient.\nconstraints = []\nfor i, nutrient in enumerate(nutrients):\n    # X * Nutrient >= 0\n    constraints.append(solver.Constraint(nutrient[1], solver.infinity()))\n    for j, item in enumerate(data):\n        # Define all the Nutrients and Weights\n        constraints[i].SetCoefficient(foods[j], item[i + 3])\n\nprint('Number of variables =', solver.NumVariables())\nprint('Number of constraints =', solver.NumConstraints())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Solving the Problem"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Solve the system.\nstatus = solver.Solve()\n# Check that the problem has an optimal solution.\nif status != pywraplp.Solver.OPTIMAL:\n    print(\"The problem does not have an optimal solution!\")\n    exit(1)\n\nnutrients_result = [0] * len(nutrients)\nprint('')\nprint('Annual Foods:')\nfor i, food in enumerate(foods):\n    if food.solution_value() > 0.0:\n        print('{}: ${}'.format(data[i][0], 365. * food.solution_value()))\n    for j, nutrient in enumerate(nutrients):\n        nutrients_result[j] += data[i][j + 3] * food.solution_value()\nprint('')\nprint('Optimal annual price: ${:.4f}'.format(365. * objective.Value()))\nprint('')\nprint('Nutrients per day:')\nfor i, nutrient in enumerate(nutrients):\n    print('{}: {:.2f} (min {})'.format(nutrient[0], nutrients_result[i],\n                                       nutrient[1]))\nprint('')\nprint('Advanced usage:')\nprint('Problem solved in ', solver.wall_time(), ' milliseconds')\nprint('Problem solved in ', solver.iterations(), ' iterations')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sodoku Problem\n\n[Link to Github Source](https://github.com/google/or-tools/blob/stable/examples/python/sudoku_sat.py)\n\nI think that this is where optimization problems get really interesting. They are flexible enough to be applied to a wide range of problems, assuming that you are able to find a crafty way to frame the problem correctly.\n\n### Framing the Problem"},{"metadata":{"trusted":true},"cell_type":"code","source":"from ortools.sat.python import cp_model\n\n# Create the model.\nmodel = cp_model.CpModel()\n\ncell_size = 3\nline_size = cell_size**2\nline = list(range(0, line_size))\ncell = list(range(0, cell_size))\n\ninitial_grid = [[0, 6, 0, 0, 5, 0, 0, 2, 0],\n                [0, 0, 0, 3, 0, 0, 0, 9, 0],\n                [7, 0, 0, 6, 0, 0, 0, 1, 0],\n                [0, 0, 6, 0, 3, 0, 4, 0, 0],\n                [0, 0, 4, 0, 7, 0, 1, 0, 0],\n                [0, 0, 5, 0, 9, 0, 8, 0, 0],\n                [0, 4, 0, 0, 0, 1, 0, 0, 6],\n                [0, 3, 0, 0, 0, 8, 0, 0, 0],\n                [0, 2, 0, 0, 4, 0, 0, 5, 0]]\n\nprint(np.matrix(initial_grid))\n\ngrid = {}\nfor i in line:\n    for j in line:\n        grid[(i, j)] = model.NewIntVar(1, line_size, 'grid %i %i' % (i, j))\n\n# AllDifferent on rows.\nfor i in line:\n    model.AddAllDifferent([grid[(i, j)] for j in line])\n\n# AllDifferent on columns.\nfor j in line:\n    model.AddAllDifferent([grid[(i, j)] for i in line])\n\n# AllDifferent on cells.\nfor i in cell:\n    for j in cell:\n        one_cell = []\n        for di in cell:\n            for dj in cell:\n                one_cell.append(grid[(i * cell_size + di,\n                                      j * cell_size + dj)])\n\n        model.AddAllDifferent(one_cell)\n\n# Initial values.\nfor i in line:\n    for j in line:\n        if initial_grid[i][j]:\n            model.Add(grid[(i, j)] == initial_grid[i][j])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Solving the Problem"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Solve and print out the solution.\nsolver = cp_model.CpSolver()\nstatus = solver.Solve(model)\nif status == cp_model.FEASIBLE:\n    for i in line:\n        print([int(solver.Value(grid[(i, j)])) for j in line])\n        \nprint('\\nAdvanced usage:')\nprint(solver.ResponseStats())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Clustering Problem\n\n[Link to Github](https://github.com/google/or-tools/blob/stable/examples/python/clustering_sat.py)"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Define Data\")\n\ndistance_matrix = [\n    [0, 10938, 4542, 2835, 29441, 2171, 1611, 9208, 9528, 11111, 16120, 22606, 22127, 20627, 21246, 23387, 16697, 33609, 26184, 24772, 22644, 20655, 30492, 23296, 32979, 18141, 19248, 17129, 17192, 15645, 12658, 11210, 12094, 13175, 18162, 4968, 12308, 10084, 13026, 15056],\n    [10938, 0, 6422, 9742, 18988, 12974, 11216, 19715, 19004, 18271, 25070, 31971, 31632, 30571, 31578, 33841, 27315, 43964, 36944, 35689, 33569, 31481, 41360, 33760, 43631, 28730, 29976, 27803, 28076, 26408, 23504, 22025, 22000, 13197, 14936, 15146, 23246, 20956, 23963, 25994],\n    [4542, 6422, 0, 3644, 25173, 6552, 5092, 13584, 13372, 13766, 19805, 26537, 26117, 24804, 25590, 27784, 21148, 37981, 30693, 29315, 27148, 25071, 34943, 27472, 37281, 22389, 23592, 21433, 21655, 20011, 17087, 15612, 15872, 11653, 15666, 8842, 16843, 14618, 17563, 19589],\n    [2835, 9742, 3644, 0, 28681, 3851, 4341, 11660, 12294, 13912, 18893, 25283, 24777, 23173, 23636, 25696, 18950, 35927, 28233, 26543, 24127, 21864, 31765, 24018, 33904, 19005, 20295, 18105, 18551, 16763, 13958, 12459, 12296, 10370, 15331, 5430, 14044, 12135, 14771, 16743],\n    [29441, 18988, 25173, 28681, 0, 31590, 29265, 37173, 35501, 32929, 40239, 47006, 46892, 46542, 48112, 50506, 44539, 60103, 54208, 53557, 51878, 50074, 59849, 52645, 62415, 47544, 48689, 46560, 46567, 45086, 42083, 40648, 40971, 29929, 28493, 34015, 41473, 38935, 42160, 44198],\n    [2171, 12974, 6552, 3851, 31590, 0, 3046, 7856, 8864, 11330, 15411, 21597, 21065, 19382, 19791, 21845, 15099, 32076, 24425, 22848, 20600, 18537, 28396, 21125, 30825, 15975, 17101, 14971, 15104, 13503, 10544, 9080, 9983, 13435, 18755, 2947, 10344, 8306, 11069, 13078],\n    [1611, 11216, 5092, 4341, 29265, 3046, 0, 8526, 8368, 9573, 14904, 21529, 21085, 19719, 20504, 22713, 16118, 32898, 25728, 24541, 22631, 20839, 30584, 23755, 33278, 18557, 19545, 17490, 17309, 15936, 12881, 11498, 12944, 14711, 19589, 5993, 12227, 9793, 12925, 14967],\n    [9208, 19715, 13584, 11660, 37173, 7856, 8526, 0, 3248, 7855, 8245, 13843, 13272, 11526, 12038, 14201, 7599, 24411, 17259, 16387, 15050, 13999, 23134, 17899, 26460, 12894, 13251, 11680, 10455, 9997, 7194, 6574, 10678, 20959, 26458, 8180, 5255, 2615, 5730, 7552],\n    [9528, 19004, 13372, 12294, 35501, 8864, 8368, 3248, 0, 4626, 6598, 13168, 12746, 11567, 12731, 15083, 9120, 25037, 18718, 18433, 17590, 16888, 25630, 20976, 29208, 16055, 16300, 14838, 13422, 13165, 10430, 9813, 13777, 22300, 27564, 10126, 8388, 5850, 8778, 10422],\n    [11111, 18271, 13766, 13912, 32929, 11330, 9573, 7855, 4626, 0, 7318, 14185, 14005, 13655, 15438, 17849, 12839, 27179, 21947, 22230, 21814, 21366, 29754, 25555, 33535, 20674, 20872, 19457, 17961, 17787, 15048, 14372, 18115, 24280, 29101, 13400, 13008, 10467, 13375, 14935],\n    [16120, 25070, 19805, 18893, 40239, 15411, 14904, 8245, 6598, 7318, 0, 6939, 6702, 6498, 8610, 10961, 7744, 19889, 15350, 16403, 16975, 17517, 24357, 22176, 28627, 18093, 17672, 16955, 14735, 15510, 13694, 13768, 18317, 28831, 34148, 16326, 11276, 9918, 11235, 11891],\n    [22606, 31971, 26537, 25283, 47006, 21597, 21529, 13843, 13168, 14185, 6939, 0, 793, 3401, 5562, 6839, 8923, 13433, 11264, 13775, 15853, 17629, 21684, 22315, 26411, 19539, 18517, 18636, 16024, 17632, 16948, 17587, 22131, 34799, 40296, 21953, 14739, 14568, 14366, 14002],\n    [22127, 31632, 26117, 24777, 46892, 21065, 21085, 13272, 12746, 14005, 6702, 793, 0, 2608, 4809, 6215, 8151, 13376, 10702, 13094, 15099, 16845, 21039, 21535, 25744, 18746, 17725, 17845, 15232, 16848, 16197, 16859, 21391, 34211, 39731, 21345, 14006, 13907, 13621, 13225],\n    [20627, 30571, 24804, 23173, 46542, 19382, 19719, 11526, 11567, 13655, 6498, 3401, 2608, 0, 2556, 4611, 5630, 13586, 9157, 11005, 12681, 14285, 19044, 18996, 23644, 16138, 15126, 15240, 12625, 14264, 13736, 14482, 18958, 32292, 37879, 19391, 11621, 11803, 11188, 10671],\n    [21246, 31578, 25590, 23636, 48112, 19791, 20504, 12038, 12731, 15438, 8610, 5562, 4809, 2556, 0, 2411, 4917, 12395, 6757, 8451, 10292, 12158, 16488, 16799, 21097, 14374, 13194, 13590, 10943, 12824, 12815, 13779, 18042, 32259, 37918, 19416, 10975, 11750, 10424, 9475],\n    [23387, 33841, 27784, 25696, 50506, 21845, 22713, 14201, 15083, 17849, 10961, 6839, 6215, 4611, 2411, 0, 6760, 10232, 4567, 7010, 9607, 12003, 14846, 16408, 19592, 14727, 13336, 14109, 11507, 13611, 14104, 15222, 19237, 34013, 39703, 21271, 12528, 13657, 11907, 10633],\n    [16697, 27315, 21148, 18950, 44539, 15099, 16118, 7599, 9120, 12839, 7744, 8923, 8151, 5630, 4917, 6760, 0, 16982, 9699, 9400, 9302, 9823, 16998, 14534, 21042, 10911, 10190, 9900, 7397, 8758, 8119, 8948, 13353, 27354, 33023, 14542, 6106, 6901, 5609, 5084],\n    [33609, 43964, 37981, 35927, 60103, 32076, 32898, 24411, 25037, 27179, 19889, 13433, 13376, 13586, 12395, 10232, 16982, 0, 8843, 12398, 16193, 19383, 16423, 22583, 20997, 22888, 21194, 22640, 20334, 22636, 23801, 25065, 28675, 44048, 49756, 31426, 22528, 23862, 21861, 20315],\n    [26184, 36944, 30693, 28233, 54208, 24425, 25728, 17259, 18718, 21947, 15350, 11264, 10702, 9157, 6757, 4567, 9699, 8843, 0, 3842, 7518, 10616, 10666, 14237, 15515, 14053, 12378, 13798, 11537, 13852, 15276, 16632, 19957, 35660, 41373, 23361, 14333, 16125, 13624, 11866],\n    [24772, 35689, 29315, 26543, 53557, 22848, 24541, 16387, 18433, 22230, 16403, 13775, 13094, 11005, 8451, 7010, 9400, 12398, 3842, 0, 3795, 7014, 8053, 10398, 12657, 10633, 8889, 10569, 8646, 10938, 12906, 14366, 17106, 33171, 38858, 21390, 12507, 14748, 11781, 9802],\n    [22644, 33569, 27148, 24127, 51878, 20600, 22631, 15050, 17590, 21814, 16975, 15853, 15099, 12681, 10292, 9607, 9302, 16193, 7518, 3795, 0, 3250, 8084, 6873, 11763, 6949, 5177, 7050, 5619, 7730, 10187, 11689, 13792, 30012, 35654, 18799, 10406, 12981, 9718, 7682],\n    [20655, 31481, 25071, 21864, 50074, 18537, 20839, 13999, 16888, 21366, 17517, 17629, 16845, 14285, 12158, 12003, 9823, 19383, 10616, 7014, 3250, 0, 9901, 4746, 12531, 3737, 1961, 4036, 3588, 5109, 7996, 9459, 10846, 27094, 32690, 16451, 8887, 11624, 8304, 6471],\n    [30492, 41360, 34943, 31765, 59849, 28396, 30584, 23134, 25630, 29754, 24357, 21684, 21039, 19044, 16488, 14846, 16998, 16423, 10666, 8053, 8084, 9901, 0, 9363, 4870, 13117, 11575, 13793, 13300, 15009, 17856, 19337, 20454, 36551, 42017, 26352, 18403, 21033, 17737, 15720],\n    [23296, 33760, 27472, 24018, 52645, 21125, 23755, 17899, 20976, 25555, 22176, 22315, 21535, 18996, 16799, 16408, 14534, 22583, 14237, 10398, 6873, 4746, 9363, 0, 10020, 5211, 4685, 6348, 7636, 8010, 11074, 12315, 11926, 27537, 32880, 18634, 12644, 15358, 12200, 10674],\n    [32979, 43631, 37281, 33904, 62415, 30825, 33278, 26460, 29208, 33535, 28627, 26411, 25744, 23644, 21097, 19592, 21042, 20997, 15515, 12657, 11763, 12531, 4870, 10020, 0, 14901, 13738, 15855, 16118, 17348, 20397, 21793, 21936, 37429, 42654, 28485, 21414, 24144, 20816, 18908],\n    [18141, 28730, 22389, 19005, 47544, 15975, 18557, 12894, 16055, 20674, 18093, 19539, 18746, 16138, 14374, 14727, 10911, 22888, 14053, 10633, 6949, 3737, 13117, 5211, 14901, 0, 1777, 1217, 3528, 2896, 5892, 7104, 7338, 23517, 29068, 13583, 7667, 10304, 7330, 6204],\n    [19248, 29976, 23592, 20295, 48689, 17101, 19545, 13251, 16300, 20872, 17672, 18517, 17725, 15126, 13194, 13336, 10190, 21194, 12378, 8889, 5177, 1961, 11575, 4685, 13738, 1777, 0, 2217, 2976, 3610, 6675, 8055, 8965, 25197, 30774, 14865, 8007, 10742, 7532, 6000],\n    [17129, 27803, 21433, 18105, 46560, 14971, 17490, 11680, 14838, 19457, 16955, 18636, 17845, 15240, 13590, 14109, 9900, 22640, 13798, 10569, 7050, 4036, 13793, 6348, 15855, 1217, 2217, 0, 2647, 1686, 4726, 6000, 6810, 23060, 28665, 12674, 6450, 9094, 6117, 5066],\n    [17192, 28076, 21655, 18551, 46567, 15104, 17309, 10455, 13422, 17961, 14735, 16024, 15232, 12625, 10943, 11507, 7397, 20334, 11537, 8646, 5619, 3588, 13300, 7636, 16118, 3528, 2976, 2647, 0, 2320, 4593, 6093, 8479, 24542, 30219, 13194, 5301, 8042, 4735, 3039],\n    [15645, 26408, 20011, 16763, 45086, 13503, 15936, 9997, 13165, 17787, 15510, 17632, 16848, 14264, 12824, 13611, 8758, 22636, 13852, 10938, 7730, 5109, 15009, 8010, 17348, 2896, 3610, 1686, 2320, 0, 3086, 4444, 6169, 22301, 27963, 11344, 4780, 7408, 4488, 3721],\n    [12658, 23504, 17087, 13958, 42083, 10544, 12881, 7194, 10430, 15048, 13694, 16948, 16197, 13736, 12815, 14104, 8119, 23801, 15276, 12906, 10187, 7996, 17856, 11074, 20397, 5892, 6675, 4726, 4593, 3086, 0, 1501, 5239, 20390, 26101, 8611, 2418, 4580, 2599, 3496],\n    [11210, 22025, 15612, 12459, 40648, 9080, 11498, 6574, 9813, 14372, 13768, 17587, 16859, 14482, 13779, 15222, 8948, 25065, 16632, 14366, 11689, 9459, 19337, 12315, 21793, 7104, 8055, 6000, 6093, 4444, 1501, 0, 4608, 19032, 24747, 7110, 2860, 4072, 3355, 4772],\n    [12094, 22000, 15872, 12296, 40971, 9983, 12944, 10678, 13777, 18115, 18317, 22131, 21391, 18958, 18042, 19237, 13353, 28675, 19957, 17106, 13792, 10846, 20454, 11926, 21936, 7338, 8965, 6810, 8479, 6169, 5239, 4608, 0, 16249, 21866, 7146, 7403, 8446, 7773, 8614],\n    [13175, 13197, 11653, 10370, 29929, 13435, 14711, 20959, 22300, 24280, 28831, 34799, 34211, 32292, 32259, 34013, 27354, 44048, 35660, 33171, 30012, 27094, 36551, 27537, 37429, 23517, 25197, 23060, 24542, 22301, 20390, 19032, 16249, 0, 5714, 12901, 21524, 20543, 22186, 23805],\n    [18162, 14936, 15666, 15331, 28493, 18755, 19589, 26458, 27564, 29101, 34148, 40296, 39731, 37879, 37918, 39703, 33023, 49756, 41373, 38858, 35654, 32690, 42017, 32880, 42654, 29068, 30774, 28665, 30219, 27963, 26101, 24747, 21866, 5714, 0, 18516, 27229, 26181, 27895, 29519],\n    [4968, 15146, 8842, 5430, 34015, 2947, 5993, 8180, 10126, 13400, 16326, 21953, 21345, 19391, 19416, 21271, 14542, 31426, 23361, 21390, 18799, 16451, 26352, 18634, 28485, 13583, 14865, 12674, 13194, 11344, 8611, 7110, 7146, 12901, 18516, 0, 9029, 7668, 9742, 11614],\n    [12308, 23246, 16843, 14044, 41473, 10344, 12227, 5255, 8388, 13008, 11276, 14739, 14006, 11621, 10975, 12528, 6106, 22528, 14333, 12507, 10406, 8887, 18403, 12644, 21414, 7667, 8007, 6450, 5301, 4780, 2418, 2860, 7403, 21524, 27229, 9029, 0, 2747, 726, 2749],\n    [10084, 20956, 14618, 12135, 38935, 8306, 9793, 2615, 5850, 10467, 9918, 14568, 13907, 11803, 11750, 13657, 6901, 23862, 16125, 14748, 12981, 11624, 21033, 15358, 24144, 10304, 10742, 9094, 8042, 7408, 4580, 4072, 8446, 20543, 26181, 7668, 2747, 0, 3330, 5313],\n    [13026, 23963, 17563, 14771, 42160, 11069, 12925, 5730, 8778, 13375, 11235, 14366, 13621, 11188, 10424, 11907, 5609, 21861, 13624, 11781, 9718, 8304, 17737, 12200, 20816, 7330, 7532, 6117, 4735, 4488, 2599, 3355, 7773, 22186, 27895, 9742, 726, 3330, 0, 2042],\n    [15056, 25994, 19589, 16743, 44198, 13078, 14967, 7552, 10422, 14935, 11891, 14002, 13225, 10671, 9475, 10633, 5084, 20315, 11866, 9802, 7682, 6471, 15720, 10674, 18908, 6204, 6000, 5066, 3039, 3721, 3496, 4772, 8614, 23805, 29519, 11614, 2749, 5313, 2042, 0],\n ] # yapf: disable\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from ortools.sat.python import cp_model\n\nnum_nodes = len(distance_matrix)\nprint('Num nodes =', num_nodes)\n\n# Number of groups to split the nodes, must divide num_nodes.\nnum_groups = 4\ngroup_size = num_nodes // num_groups\n\n# Model.\nmodel = cp_model.CpModel()\n\n# Variables.\nneighbors = {}\nobj_vars = []\nobj_coeffs = []\nfor n1 in range(num_nodes - 1):\n    for n2 in range(n1 + 1, num_nodes):\n        same = model.NewBoolVar('neighbors_%i_%i' % (n1, n2))\n        neighbors[n1, n2] = same\n        obj_vars.append(same)\n        obj_coeffs.append(distance_matrix[n1][n2] + distance_matrix[n2][n1])\n\n# Number of neighborss:\nfor n in range(num_nodes):\n    model.Add(sum(neighbors[m, n] for m in range(n)) + \n              sum(neighbors[n, m] for m in range(n + 1, num_nodes)) ==\n              group_size - 1)\n\n# Enforce transivity on all triplets.\nfor n1 in range(num_nodes - 2):\n    for n2 in range(n1 + 1, num_nodes - 1):\n        for n3 in range(n2 + 1, num_nodes):\n            model.Add(\n                neighbors[n1, n3] + neighbors[n2, n3] + neighbors[n1, n2] != 2)\n\n# Redundant constraints on total sum of neighborss.\nmodel.Add(sum(obj_vars) == num_groups * group_size * (group_size - 1) // 2)\n\n# Minimize weighted sum of arcs.\nmodel.Minimize(\n    sum(obj_vars[i] * obj_coeffs[i] for i in range(len(obj_vars))))\n\n# Solve and print out the solution.\nsolver = cp_model.CpSolver()\nsolver.parameters.log_search_progress = True\nsolver.parameters.num_search_workers = 8\n\nstatus = solver.Solve(model)\nprint(solver.ResponseStats())\n\nvisited = set()\nfor g in range(num_groups):\n    for n in range(num_nodes):\n        if not n in visited:\n            visited.add(n)\n            output = str(n)\n            for o in range(n + 1, num_nodes):\n                if solver.BooleanValue(neighbors[n, o]):\n                    visited.add(o)\n                    output += ' ' + str(o)\n            print('Group', g, ':', output)\n            break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Zebra Problem\n\n[Link to Github](https://github.com/google/or-tools/blob/stable/examples/python/zebra_sat.py)\n\n**This is the zebra problem as invented by Lewis Caroll:** <br>\n<br>\nThere are five houses.<br>\nThe Englishman lives in the red house.<br>\nThe Spaniard owns the dog.<br>\nCoffee is drunk in the green house.<br>\nThe Ukrainian drinks tea.<br>\nThe green house is immediately to the right of the ivory house.<br>\nThe Old Gold smoker owns snails.<br>\nKools are smoked in the yellow house.<br>\nMilk is drunk in the middle house.<br>\nThe Norwegian lives in the first house.<br>\nThe man who smokes Chesterfields lives in the house next to the man with the fox.<br>\nKools are smoked in the house next to the house where the horse is kept.<br>\nThe Lucky Strike smoker drinks orange juice.<br>\nThe Japanese smokes Parliaments.<br>\nThe Norwegian lives next to the blue house.<br>\nWho owns a zebra and who drinks water?<br>"},{"metadata":{"trusted":true},"cell_type":"code","source":"from ortools.sat.python import cp_model\n\n# pylint: disable=too-many-statements\n# Create the model.\nmodel = cp_model.CpModel()\n\nred = model.NewIntVar(1, 5, 'red')\ngreen = model.NewIntVar(1, 5, 'green')\nyellow = model.NewIntVar(1, 5, 'yellow')\nblue = model.NewIntVar(1, 5, 'blue')\nivory = model.NewIntVar(1, 5, 'ivory')\n\nenglishman = model.NewIntVar(1, 5, 'englishman')\nspaniard = model.NewIntVar(1, 5, 'spaniard')\njapanese = model.NewIntVar(1, 5, 'japanese')\nukrainian = model.NewIntVar(1, 5, 'ukrainian')\nnorwegian = model.NewIntVar(1, 5, 'norwegian')\n\ndog = model.NewIntVar(1, 5, 'dog')\nsnails = model.NewIntVar(1, 5, 'snails')\nfox = model.NewIntVar(1, 5, 'fox')\nzebra = model.NewIntVar(1, 5, 'zebra')\nhorse = model.NewIntVar(1, 5, 'horse')\n\ntea = model.NewIntVar(1, 5, 'tea')\ncoffee = model.NewIntVar(1, 5, 'coffee')\nwater = model.NewIntVar(1, 5, 'water')\nmilk = model.NewIntVar(1, 5, 'milk')\nfruit_juice = model.NewIntVar(1, 5, 'fruit juice')\n\nold_gold = model.NewIntVar(1, 5, 'old gold')\nkools = model.NewIntVar(1, 5, 'kools')\nchesterfields = model.NewIntVar(1, 5, 'chesterfields')\nlucky_strike = model.NewIntVar(1, 5, 'lucky strike')\nparliaments = model.NewIntVar(1, 5, 'parliaments')\n\nmodel.AddAllDifferent([red, green, yellow, blue, ivory])\nmodel.AddAllDifferent(\n    [englishman, spaniard, japanese, ukrainian, norwegian])\nmodel.AddAllDifferent([dog, snails, fox, zebra, horse])\nmodel.AddAllDifferent([tea, coffee, water, milk, fruit_juice])\nmodel.AddAllDifferent(\n    [parliaments, kools, chesterfields, lucky_strike, old_gold])\n\nmodel.Add(englishman == red)\nmodel.Add(spaniard == dog)\nmodel.Add(coffee == green)\nmodel.Add(ukrainian == tea)\nmodel.Add(green == ivory + 1)\nmodel.Add(old_gold == snails)\nmodel.Add(kools == yellow)\nmodel.Add(milk == 3)\nmodel.Add(norwegian == 1)\n\ndiff_fox_chesterfields = model.NewIntVar(-4, 4, 'diff_fox_chesterfields')\nmodel.Add(diff_fox_chesterfields == fox - chesterfields)\nmodel.AddAbsEquality(1, diff_fox_chesterfields)\n\ndiff_horse_kools = model.NewIntVar(-4, 4, 'diff_horse_kools')\nmodel.Add(diff_horse_kools == horse - kools)\nmodel.AddAbsEquality(1, diff_horse_kools)\n\nmodel.Add(lucky_strike == fruit_juice)\nmodel.Add(japanese == parliaments)\n\ndiff_norwegian_blue = model.NewIntVar(-4, 4, 'diff_norwegian_blue')\nmodel.Add(diff_norwegian_blue == norwegian - blue)\nmodel.AddAbsEquality(1, diff_norwegian_blue)\n\n# Solve and print out the solution.\nsolver = cp_model.CpSolver()\nstatus = solver.Solve(model)\n\nif status == cp_model.FEASIBLE:\n    people = [englishman, spaniard, japanese, ukrainian, norwegian]\n    water_drinker = [\n        p for p in people if solver.Value(p) == solver.Value(water)\n    ][0]\n    zebra_owner = [\n        p for p in people if solver.Value(p) == solver.Value(zebra)\n    ][0]\n    print('The', water_drinker.Name(), 'drinks water.')\n    print('The', zebra_owner.Name(), 'owns the zebra.')\n    \n    print(\"\\n\")\n    print(solver.ResponseStats())\nelse:\n    print('No solutions to the zebra problem, this is unusual!')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Travelling Salesman Problem Between Cities\n\nThis is a classic routing problem, where a vehicle must travel through certain locations in the most efficient way. The framing of the problem below uses a specialised routing since this problem can quickly become intractable. But that does no mean it can't be defined in a classic objective and constraint framework.\n\nThe input data is formated as a `N*N` matrix, where N is the location to be routed through. Routing problems can also be though of as graphs, where vertices are locations, and edges are the paths (weighted with distance or cost) between them. \n\nIn this particular problem, the node is fully connected, which is why there is a distance value between each node. In network theory, this matrix representation is also known as an _adjacency matrix_."},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"Simple travelling salesman problem between cities.\"\"\"\n\nfrom ortools.constraint_solver import routing_enums_pb2\nfrom ortools.constraint_solver import pywrapcp\n\n\ndef create_data_model():\n    \"\"\"Stores the data for the problem.\"\"\"\n    data = {}\n    data['distance_matrix'] = [\n        [0, 2451, 713, 1018, 1631, 1374, 2408, 213, 2571, 875, 1420, 2145, 1972],\n        [2451, 0, 1745, 1524, 831, 1240, 959, 2596, 403, 1589, 1374, 357, 579],\n        [713, 1745, 0, 355, 920, 803, 1737, 851, 1858, 262, 940, 1453, 1260],\n        [1018, 1524, 355, 0, 700, 862, 1395, 1123, 1584, 466, 1056, 1280, 987],\n        [1631, 831, 920, 700, 0, 663, 1021, 1769, 949, 796, 879, 586, 371],\n        [1374, 1240, 803, 862, 663, 0, 1681, 1551, 1765, 547, 225, 887, 999],\n        [2408, 959, 1737, 1395, 1021, 1681, 0, 2493, 678, 1724, 1891, 1114, 701],\n        [213, 2596, 851, 1123, 1769, 1551, 2493, 0, 2699, 1038, 1605, 2300, 2099],\n        [2571, 403, 1858, 1584, 949, 1765, 678, 2699, 0, 1744, 1645, 653, 600],\n        [875, 1589, 262, 466, 796, 547, 1724, 1038, 1744, 0, 679, 1272, 1162],\n        [1420, 1374, 940, 1056, 879, 225, 1891, 1605, 1645, 679, 0, 1017, 1200],\n        [2145, 357, 1453, 1280, 586, 887, 1114, 2300, 653, 1272, 1017, 0, 504],\n        [1972, 579, 1260, 987, 371, 999, 701, 2099, 600, 1162, 1200, 504, 0],\n    ]  # yapf: disable\n    data['num_vehicles'] = 1\n    data['depot'] = 0\n    return data\n\n\n\"\"\"Entry point of the program.\"\"\"\n# Instantiate the data problem.\ndata = create_data_model()\ndisplay(pd.DataFrame(data['distance_matrix']).style.background_gradient(cmap='gray'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def print_solution(manager, routing, assignment):\n    \"\"\"Prints assignment on console.\"\"\"\n    print('Objective: {} miles'.format(assignment.ObjectiveValue()))\n    index = routing.Start(0)\n    plan_output = 'Route for vehicle 0:\\n'\n    route_distance = 0\n    while not routing.IsEnd(index):\n        plan_output += ' {} ->'.format(manager.IndexToNode(index))\n        previous_index = index\n        index = assignment.Value(routing.NextVar(index))\n        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)\n    plan_output += ' {}\\n'.format(manager.IndexToNode(index))\n    print(plan_output)\n    plan_output += 'Route distance: {}miles\\n'.format(route_distance)\n\n\n# Create the routing index manager.\nmanager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),\n                                       data['num_vehicles'], data['depot'])\n\n# Create Routing Model.\nrouting = pywrapcp.RoutingModel(manager)\n\n\ndef distance_callback(from_index, to_index):\n    \"\"\"Returns the distance between the two nodes.\"\"\"\n    # Convert from routing variable Index to distance matrix NodeIndex.\n    from_node = manager.IndexToNode(from_index)\n    to_node = manager.IndexToNode(to_index)\n    return data['distance_matrix'][from_node][to_node]\n\ntransit_callback_index = routing.RegisterTransitCallback(distance_callback)\n\n# Define cost of each arc.\nrouting.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)\n\n# Setting first solution heuristic.\nsearch_parameters = pywrapcp.DefaultRoutingSearchParameters()\nsearch_parameters.first_solution_strategy = (\n    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n\n# Solve the problem.\nassignment = routing.SolveWithParameters(search_parameters)\n\n# Print solution on console.\nif assignment:\n    print_solution(manager, routing, assignment)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#### Circuit Board Drilling Sequence Optimisation"},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"Simple travelling salesman problem on a circuit board.\"\"\"\n\nfrom __future__ import print_function\nimport math\nfrom ortools.constraint_solver import routing_enums_pb2\nfrom ortools.constraint_solver import pywrapcp\n\n\ndef create_data_model():\n    \"\"\"Stores the data for the problem.\"\"\"\n    data = {}\n    # Locations in block units\n    data['locations'] = [\n        (288, 149), (288, 129), (270, 133), (256, 141), (256, 157), (246, 157),\n        (236, 169), (228, 169), (228, 161), (220, 169), (212, 169), (204, 169),\n        (196, 169), (188, 169), (196, 161), (188, 145), (172, 145), (164, 145),\n        (156, 145), (148, 145), (140, 145), (148, 169), (164, 169), (172, 169),\n        (156, 169), (140, 169), (132, 169), (124, 169), (116, 161), (104, 153),\n        (104, 161), (104, 169), (90, 165), (80, 157), (64, 157), (64, 165),\n        (56, 169), (56, 161), (56, 153), (56, 145), (56, 137), (56, 129),\n        (56, 121), (40, 121), (40, 129), (40, 137), (40, 145), (40, 153),\n        (40, 161), (40, 169), (32, 169), (32, 161), (32, 153), (32, 145),\n        (32, 137), (32, 129), (32, 121), (32, 113), (40, 113), (56, 113),\n        (56, 105), (48, 99), (40, 99), (32, 97), (32, 89), (24, 89),\n        (16, 97), (16, 109), (8, 109), (8, 97), (8, 89), (8, 81),\n        (8, 73), (8, 65), (8, 57), (16, 57), (8, 49), (8, 41),\n        (24, 45), (32, 41), (32, 49), (32, 57), (32, 65), (32, 73),\n        (32, 81), (40, 83), (40, 73), (40, 63), (40, 51), (44, 43),\n        (44, 35), (44, 27), (32, 25), (24, 25), (16, 25), (16, 17),\n        (24, 17), (32, 17), (44, 11), (56, 9), (56, 17), (56, 25),\n        (56, 33), (56, 41), (64, 41), (72, 41), (72, 49), (56, 49),\n        (48, 51), (56, 57), (56, 65), (48, 63), (48, 73), (56, 73),\n        (56, 81), (48, 83), (56, 89), (56, 97), (104, 97), (104, 105),\n        (104, 113), (104, 121), (104, 129), (104, 137), (104, 145), (116, 145),\n        (124, 145), (132, 145), (132, 137), (140, 137), (148, 137), (156, 137),\n        (164, 137), (172, 125), (172, 117), (172, 109), (172, 101), (172, 93),\n        (172, 85), (180, 85), (180, 77), (180, 69), (180, 61), (180, 53),\n        (172, 53), (172, 61), (172, 69), (172, 77), (164, 81), (148, 85),\n        (124, 85), (124, 93), (124, 109), (124, 125), (124, 117), (124, 101),\n        (104, 89), (104, 81), (104, 73), (104, 65), (104, 49), (104, 41),\n        (104, 33), (104, 25), (104, 17), (92, 9), (80, 9), (72, 9),\n        (64, 21), (72, 25), (80, 25), (80, 25), (80, 41), (88, 49),\n        (104, 57), (124, 69), (124, 77), (132, 81), (140, 65), (132, 61),\n        (124, 61), (124, 53), (124, 45), (124, 37), (124, 29), (132, 21),\n        (124, 21), (120, 9), (128, 9), (136, 9), (148, 9), (162, 9),\n        (156, 25), (172, 21), (180, 21), (180, 29), (172, 29), (172, 37),\n        (172, 45), (180, 45), (180, 37), (188, 41), (196, 49), (204, 57),\n        (212, 65), (220, 73), (228, 69), (228, 77), (236, 77), (236, 69),\n        (236, 61), (228, 61), (228, 53), (236, 53), (236, 45), (228, 45),\n        (228, 37), (236, 37), (236, 29), (228, 29), (228, 21), (236, 21),\n        (252, 21), (260, 29), (260, 37), (260, 45), (260, 53), (260, 61),\n        (260, 69), (260, 77), (276, 77), (276, 69), (276, 61), (276, 53),\n        (284, 53), (284, 61), (284, 69), (284, 77), (284, 85), (284, 93),\n        (284, 101), (288, 109), (280, 109), (276, 101), (276, 93), (276, 85),\n        (268, 97), (260, 109), (252, 101), (260, 93), (260, 85), (236, 85),\n        (228, 85), (228, 93), (236, 93), (236, 101), (228, 101), (228, 109),\n        (228, 117), (228, 125), (220, 125), (212, 117), (204, 109), (196, 101),\n        (188, 93), (180, 93), (180, 101), (180, 109), (180, 117), (180, 125),\n        (196, 145), (204, 145), (212, 145), (220, 145), (228, 145), (236, 145),\n        (246, 141), (252, 125), (260, 129), (280, 133)\n    ]  # yapf: disable\n    data['num_vehicles'] = 1\n    data['depot'] = 0\n    return data\n\n\ndef compute_euclidean_distance_matrix(locations):\n    \"\"\"Creates callback to return distance between points.\"\"\"\n    distances = {}\n    for from_counter, from_node in enumerate(locations):\n        distances[from_counter] = {}\n        for to_counter, to_node in enumerate(locations):\n            if from_counter == to_counter:\n                distances[from_counter][to_counter] = 0\n            else:\n                # Euclidean distance\n                distances[from_counter][to_counter] = (int(\n                    math.hypot((from_node[0] - to_node[0]),\n                               (from_node[1] - to_node[1]))))\n    return distances\n\n\ndef print_solution(manager, routing, assignment):\n    \"\"\"Prints assignment on console.\"\"\"\n    print('Objective: {}'.format(assignment.ObjectiveValue()))\n    path = []\n    index = routing.Start(0)\n    plan_output = 'Route:\\n'\n    route_distance = 0\n    while not routing.IsEnd(index):\n        next_path = manager.IndexToNode(index)\n        plan_output += ' {} ->'.format(next_path)\n        path.append(next_path)\n        previous_index = index\n        index = assignment.Value(routing.NextVar(index))\n        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)\n    plan_output += ' {}\\n'.format(manager.IndexToNode(index))\n    print(plan_output)\n    plan_output += 'Objective: {}m\\n'.format(route_distance)\n    return path\n\ndef distance_callback(from_index, to_index):\n    \"\"\"Returns the distance between the two nodes.\"\"\"\n    # Convert from routing variable Index to distance matrix NodeIndex.\n    from_node = manager.IndexToNode(from_index)\n    to_node = manager.IndexToNode(to_index)\n    return distance_matrix[from_node][to_node]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Instantiate the data problem.\ndata = create_data_model()\ncircuit = pd.DataFrame(data['locations'])\n\nf, ax = plt.subplots(figsize = [10,8])\nplt.scatter(x = circuit.iloc[:,0], y = circuit.iloc[:,1], c='r')\nplt.title(\"Circuit Board\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create the routing index manager.\nmanager = pywrapcp.RoutingIndexManager(len(data['locations']),\n                                       data['num_vehicles'], data['depot'])\n\n# Create Routing Model.\nrouting = pywrapcp.RoutingModel(manager)\n\ndistance_matrix = compute_euclidean_distance_matrix(data['locations'])\n\ntransit_callback_index = routing.RegisterTransitCallback(distance_callback)\n\n# Define cost of each arc.\nrouting.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)\n\n# Setting first solution heuristic.\nsearch_parameters = pywrapcp.DefaultRoutingSearchParameters()\nsearch_parameters.first_solution_strategy = (\n    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n\n# Solve the problem.\nassignment = routing.SolveWithParameters(search_parameters)\n\n# Print solution on console.\nif assignment:\n    path_solution = print_solution(manager, routing, assignment)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(figsize = [10,8])\nax.plot(circuit.iloc[path_solution,0], circuit.iloc[path_solution,1], 'r-o')\nax.plot(circuit.iloc[path_solution,0][0], circuit.iloc[path_solution,1][0], 'k-v', markersize=15)\nax.plot(circuit.iloc[path_solution[-1],0], circuit.iloc[path_solution[-1],1], 'b-X', markersize=15)\nplt.title(\"Circuit Board - Triangle is Start, X is End\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Multiple Vehicle Routing Problem"},{"metadata":{"trusted":true},"cell_type":"code","source":"coordinates = [\n(456, 320), # location 0 - the depot\n(228, 0),    # location 1\n(912, 0),    # location 2\n(0, 80),     # location 3\n(114, 80),   # location 4\n(570, 160),  # location 5\n(798, 160),  # location 6\n(342, 240),  # location 7\n(684, 240),  # location 8\n(570, 400),  # location 9\n(912, 400),  # location 10\n(114, 480),  # location 11\n(228, 480),  # location 12\n(342, 560),  # location 13\n(684, 560),  # location 14\n(0, 640),    # location 15\n(798, 640)]  # location 16\n\nX = np.array([x[0] for x in coordinates])\nY = np.array([x[1] for x in coordinates])\n\nf, ax = plt.subplots(figsize = [8,6])\nfor i, txt in enumerate(coordinates):\n    ax.text(X[i] + 15, Y[i]+ 15, f\"{i}\")\n\nax.plot(X, Y, 'ko', markersize=8)\nax.plot(X[0], Y[0], 'gX', markersize=15)\nax.set_title(\"Routing Problem Locations\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"Vehicles Routing Problem (VRP).\"\"\"\n\nfrom __future__ import print_function\nfrom ortools.constraint_solver import routing_enums_pb2\nfrom ortools.constraint_solver import pywrapcp\n\n\ndef create_data_model():\n    \"\"\"Stores the data for the problem.\"\"\"\n    data = {}\n    data['distance_matrix'] = [\n        [\n            0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,\n            468, 776, 662\n        ],\n        [\n            548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,\n            1016, 868, 1210\n        ],\n        [\n            776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,\n            1130, 788, 1552, 754\n        ],\n        [\n            696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,\n            1164, 560, 1358\n        ],\n        [\n            582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,\n            1050, 674, 1244\n        ],\n        [\n            274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,\n            514, 1050, 708\n        ],\n        [\n            502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,\n            514, 1278, 480\n        ],\n        [\n            194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,\n            662, 742, 856\n        ],\n        [\n            308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,\n            320, 1084, 514\n        ],\n        [\n            194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,\n            274, 810, 468\n        ],\n        [\n            536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,\n            730, 388, 1152, 354\n        ],\n        [\n            502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,\n            308, 650, 274, 844\n        ],\n        [\n            388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,\n            536, 388, 730\n        ],\n        [\n            354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,\n            342, 422, 536\n        ],\n        [\n            468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,\n            342, 0, 764, 194\n        ],\n        [\n            776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,\n            388, 422, 764, 0, 798\n        ],\n        [\n            662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,\n            536, 194, 798, 0\n        ],\n    ]\n    data['num_vehicles'] = 4\n    data['depot'] = 0\n    return data\n\n\ndef print_solution(data, manager, routing, solution):\n    \"\"\"Prints solution on console.\"\"\"\n    max_route_distance = 0\n    solutions = {}\n    for vehicle_id in range(data['num_vehicles']):\n        index = routing.Start(vehicle_id)\n        plan_output = 'Route for vehicle {}:\\n'.format(vehicle_id)\n        route_distance = 0\n        vehicle_solution = []\n        while not routing.IsEnd(index):\n            vehicle_solution.append(manager.IndexToNode(index))\n            plan_output += ' {} -> '.format(manager.IndexToNode(index))\n            previous_index = index\n            index = solution.Value(routing.NextVar(index))\n            route_distance += routing.GetArcCostForVehicle(\n                previous_index, index, vehicle_id)\n        plan_output += '{}\\n'.format(manager.IndexToNode(index))\n        plan_output += 'Distance of the route: {}m\\n'.format(route_distance)\n        print(plan_output)\n        max_route_distance = max(route_distance, max_route_distance)\n        solutions[vehicle_id] = vehicle_solution\n    print('Maximum of the route distances: {}m'.format(max_route_distance))\n    return solutions\n\n\"\"\"Solve the CVRP problem.\"\"\"\n# Instantiate the data problem.\ndata = create_data_model()\n\n# Create the routing index manager.\nmanager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),\n                                       data['num_vehicles'], data['depot'])\n\n# Create Routing Model.\nrouting = pywrapcp.RoutingModel(manager)\n\n\n# Create and register a transit callback.\ndef distance_callback(from_index, to_index):\n    \"\"\"Returns the distance between the two nodes.\"\"\"\n    # Convert from routing variable Index to distance matrix NodeIndex.\n    from_node = manager.IndexToNode(from_index)\n    to_node = manager.IndexToNode(to_index)\n    return data['distance_matrix'][from_node][to_node]\n\ntransit_callback_index = routing.RegisterTransitCallback(distance_callback)\n\n# Define cost of each arc.\nrouting.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)\n\n# Add Distance constraint.\ndimension_name = 'Distance'\nrouting.AddDimension(\n    transit_callback_index,\n    0,  # no slack\n    3000,  # vehicle maximum travel distance\n    True,  # start cumul to zero\n    dimension_name)\ndistance_dimension = routing.GetDimensionOrDie(dimension_name)\ndistance_dimension.SetGlobalSpanCostCoefficient(100)\n\n# Setting first solution heuristic.\nsearch_parameters = pywrapcp.DefaultRoutingSearchParameters()\nsearch_parameters.first_solution_strategy = (\n    routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n\n# Solve the problem.\nsolution = routing.SolveWithParameters(search_parameters)\n\n# Print solution on console.\nif solution:\n    solution_dict = print_solution(data, manager, routing, solution)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(figsize = [8,6])\n\nax.plot(X, Y, 'ko', markersize=8)\nax.plot(X[0], Y[0], 'gX', markersize=30)\n\nfor i, txt in enumerate(coordinates):\n    ax.text(X[i] + 15, Y[i]+ 15, f\"{i}\")\n\nvehicle_colors = [\"g\",\"k\",\"r\", \"m\"]\nfor vehicle in solution_dict:\n    ax.plot(X[solution_dict[vehicle] + [0]], Y[solution_dict[vehicle] + [0]], f'{vehicle_colors[vehicle]}--')\n    \nax.set_title(\"Four Vehicles Routed for Depot Coverage\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Scheduling Problem\nhttps://developers.google.com/optimization/scheduling/employee_scheduling\n\n- Employee Scheduling\n- The Job Shop Problem"},{"metadata":{"trusted":true},"cell_type":"code","source":"from ortools.sat.python import cp_model\n\n# This program tries to find an optimal assignment of nurses to shifts\n# (3 shifts per day, for 7 days), subject to some constraints (see below).\n# Each nurse can request to be assigned to specific shifts.\n# The optimal assignment maximizes the number of fulfilled shift requests.\nnum_nurses = 5\nnum_shifts = 3\nnum_days = 7\nall_nurses = range(num_nurses)\nall_shifts = range(num_shifts)\nall_days = range(num_days)\nshift_requests = [\n    [[0, 0, 1], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 0, 1]],\n    [[0, 0, 0], [0, 0, 0], [0, 1, 0], [0, 1, 0], [1, 0, 0], [0, 0, 0], [0, 0, 1]],\n    [[0, 1, 0], [0, 1, 0], [0, 0, 0], [1, 0, 0], [0, 0, 0], [0, 1, 0], [0, 0, 0]],\n    [[0, 0, 1], [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 0], [1, 0, 0], [0, 0, 0]],\n    [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 0]]\n]\nrequested_df = pd.DataFrame(shift_requests)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Creates the model.\nmodel = cp_model.CpModel()\n\n# Creates shift variables.\n# shifts[(n, d, s)]: nurse 'n' works shift 's' on day 'd'.\nshifts = {}\nfor n in all_nurses:\n    for d in all_days:\n        for s in all_shifts:\n            shifts[(n, d,\n                    s)] = model.NewBoolVar('shift_n%id%is%i' % (n, d, s))\n\n# Each shift is assigned to exactly one nurse in .\nfor d in all_days:\n    for s in all_shifts:\n        model.Add(sum(shifts[(n, d, s)] for n in all_nurses) == 1)\n\n# Each nurse works at most one shift per day.\nfor n in all_nurses:\n    for d in all_days:\n        model.Add(sum(shifts[(n, d, s)] for s in all_shifts) <= 1)\n\n# Try to distribute the shifts evenly, so that each nurse works\n# min_shifts_per_nurse shifts. If this is not possible, because the total\n# number of shifts is not divisible by the number of nurses, some nurses will\n# be assigned one more shift.\nmin_shifts_per_nurse = (num_shifts * num_days) // num_nurses\nif num_shifts * num_days % num_nurses == 0:\n    max_shifts_per_nurse = min_shifts_per_nurse\nelse:\n    max_shifts_per_nurse = min_shifts_per_nurse + 1\nfor n in all_nurses:\n    num_shifts_worked = 0\n    for d in all_days:\n        for s in all_shifts:\n            num_shifts_worked += shifts[(n, d, s)]\n    model.Add(min_shifts_per_nurse <= num_shifts_worked)\n    model.Add(num_shifts_worked <= max_shifts_per_nurse)\n\n# pylint: disable=g-complex-comprehension\nmodel.Maximize(\n    sum(shift_requests[n][d][s] * shifts[(n, d, s)] for n in all_nurses\n        for d in all_days for s in all_shifts))\n# Creates the solver and solve.\nsolver = cp_model.CpSolver()\nsolver.Solve(model)\n\nresults = {d:{} for d in all_days}\nfor d in all_days:\n    print('Day', d)\n    for n in all_nurses:\n        for s in all_shifts:\n            if solver.Value(shifts[(n, d, s)]) == 1:\n                results[d][s] = n\n                if shift_requests[n][d][s] == 1:\n                    print('Nurse', n, 'works shift', s, '(requested).')\n                else:\n                    print('Nurse', n, 'works shift', s, '(not requested).')\n    print()\n\nprint(\"Results - Format {day:{shift:nurse}}\")\npprint.pprint(results)    \n\nprint(\"\\nShifts per nurse?\")\nnurse_occurence =[]\nfor r in results:\n    for s in results[r]:\n        nurse_occurence.append(results[r][s])\nprint(Counter(nurse_occurence))\n\n# Statistics.\nprint()\nprint('Statistics')\nprint('  - Number of shift requests met = %i' % solver.ObjectiveValue(),\n      '(out of', num_nurses * min_shifts_per_nurse, ')')\nprint('  - wall time       : %f s' % solver.WallTime())\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The Job Shop Problem"},{"metadata":{"trusted":true},"cell_type":"code","source":"Image(url= \"https://developers.google.com/optimization/images/scheduling/schedule1.png\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from __future__ import print_function\n\nimport collections\n\n# Import Python wrapper for or-tools CP-SAT solver.\nfrom ortools.sat.python import cp_model\n\n\ndef MinimalJobshopSat():\n    \"\"\"Minimal jobshop problem.\"\"\"\n    # Create the model.\n    model = cp_model.CpModel()\n\n    jobs_data = [  # task = (machine_id, processing_time).\n        [(0, 3), (1, 2), (2, 2)],  # Job0\n        [(0, 2), (2, 1), (1, 4)],  # Job1\n        [(1, 4), (2, 3)]  # Job2\n    ]\n\n    machines_count = 1 + max(task[0] for job in jobs_data for task in job)\n    all_machines = range(machines_count)\n\n    # Computes horizon dynamically as the sum of all durations.\n    horizon = sum(task[1] for job in jobs_data for task in job)\n\n    # Named tuple to store information about created variables.\n    task_type = collections.namedtuple('task_type', 'start end interval')\n    # Named tuple to manipulate solution information.\n    assigned_task_type = collections.namedtuple('assigned_task_type',\n                                                'start job index duration')\n\n    # Creates job intervals and add to the corresponding machine lists.\n    all_tasks = {}\n    machine_to_intervals = collections.defaultdict(list)\n\n    for job_id, job in enumerate(jobs_data):\n        for task_id, task in enumerate(job):\n            machine = task[0]\n            duration = task[1]\n            suffix = '_%i_%i' % (job_id, task_id)\n            start_var = model.NewIntVar(0, horizon, 'start' + suffix)\n            end_var = model.NewIntVar(0, horizon, 'end' + suffix)\n            interval_var = model.NewIntervalVar(start_var, duration, end_var,\n                                                'interval' + suffix)\n            all_tasks[job_id, task_id] = task_type(\n                start=start_var, end=end_var, interval=interval_var)\n            machine_to_intervals[machine].append(interval_var)\n\n    # Create and add disjunctive constraints.\n    for machine in all_machines:\n        model.AddNoOverlap(machine_to_intervals[machine])\n\n    # Precedences inside a job.\n    for job_id, job in enumerate(jobs_data):\n        for task_id in range(len(job) - 1):\n            model.Add(all_tasks[job_id, task_id +\n                                1].start >= all_tasks[job_id, task_id].end)\n\n    # Makespan objective.\n    obj_var = model.NewIntVar(0, horizon, 'makespan')\n    model.AddMaxEquality(obj_var, [\n        all_tasks[job_id, len(job) - 1].end\n        for job_id, job in enumerate(jobs_data)\n    ])\n    model.Minimize(obj_var)\n\n    # Solve model.\n    solver = cp_model.CpSolver()\n    status = solver.Solve(model)\n\n    if status == cp_model.OPTIMAL:\n        # Create one list of assigned tasks per machine.\n        assigned_jobs = collections.defaultdict(list)\n        for job_id, job in enumerate(jobs_data):\n            for task_id, task in enumerate(job):\n                machine = task[0]\n                assigned_jobs[machine].append(\n                    assigned_task_type(\n                        start=solver.Value(all_tasks[job_id, task_id].start),\n                        job=job_id,\n                        index=task_id,\n                        duration=task[1]))\n\n        # Create per machine output lines.\n        output = ''\n        for machine in all_machines:\n            # Sort by starting time.\n            assigned_jobs[machine].sort()\n            sol_line_tasks = 'Machine ' + str(machine) + ': '\n            sol_line = '           '\n\n            for assigned_task in assigned_jobs[machine]:\n                name = 'job_%i_%i' % (assigned_task.job, assigned_task.index)\n                # Add spaces to output to align columns.\n                sol_line_tasks += '%-10s' % name\n\n                start = assigned_task.start\n                duration = assigned_task.duration\n                sol_tmp = '[%i,%i]' % (start, start + duration)\n                # Add spaces to output to align columns.\n                sol_line += '%-10s' % sol_tmp\n\n            sol_line += '\\n'\n            sol_line_tasks += '\\n'\n            output += sol_line_tasks\n            output += sol_line\n\n        # Finally print the solution found.\n        print('Optimal Schedule Length: %i' % solver.ObjectiveValue())\n        print(output)\n\nMinimalJobshopSat()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Network Flow Problem\n\nhttps://developers.google.com/optimization/flow/maxflow\n\n**Maximum Flows:** <br>\nGet from Node 0 to 4, while maximizing total flow in a unit of time. How to use a networks full capacity.\n\n- At each node, other than the source or the sink, the total flow of all arcs leading in to the node equals the total flow of all arcs leading out of it."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"%%html\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"700\" height=\"450\">\n<!-- arcs -->\n<line x1=\"50\" y1=\"225\" x2=\"250\" y2=\"25\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"50\" y1=\"225\" x2=\"250\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"50\" y1=\"225\" x2=\"250\" y2=\"425\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"250\" y1=\"25\" x2=\"250\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"250\" y1=\"225\" x2=\"450\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"250\" y1=\"25\" x2=\"450\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"250\" y1=\"225\" x2=\"450\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<line x1=\"250\" y1=\"425\" x2=\"450\" y2=\"225\" stroke=\"#762a83\" stroke-width=\"3\"></line>\n<g stroke=\"#762a83\" stroke-width=\"3\" fill=\"none\">\n      <path d=\"M 250 425 a 50 155 0 0 1 0 -200\"></path> <!-- node 3 to 2 -->\n      <path d=\"M 250 425 a 50 155 0 0 0 0 -200\"></path> <!-- node 2 to 3 -->\n</g>\n<!-- weights -->\n<text x=\"120\" y=\"120\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">20</text>\n<text x=\"143\" y=\"215\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">30</text>\n<text x=\"260\" y=\"150\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">40</text>\n<text x=\"355\" y=\"119\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">30</text>\n<text x=\"330\" y=\"215\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">20</text>\n<text x=\"120\" y=\"350\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">10</text>\n<text x=\"355\" y=\"350\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">20</text>\n<text x=\"212\" y=\"329\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">5</text>\n<text x=\"274\" y=\"329\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">10</text>\n<!-- arrows -->\n<g stroke=\"#762a83\" stroke-width=\"1\" fill=\"#762a83\">\n   <path d=\"M 230 40 L 248 27 L 236 46 L 230 40\"></path> <!-- node 0 to 1 -->\n   <path d=\"M 228 221 L 250 225 L 228 229 L 228 220\"></path> <!-- node 0 to 2 -->\n   <path d=\"M 230 410 L 250 425 L 236 404 L 230 410\"></path> <!-- node 0 to 3 -->\n   <path d=\"M 428 221 L 450 225 L 428 229 L 428 220\"></path> <!-- node 2 to 4 -->\n   <path d=\"M 430 212 L 450 225 L 438 206 L 430 212\"></path> <!-- node 1 to 4 -->\n   <path d=\"M 432 235 L 450 225 L 439 242 L 432 235\"></path> <!-- node 3 to 4 -->\n   <path d=\"M 245 203 L 250 220 L 255 203 L 245 203\"></path> <!-- node 1 to 2 -->\n   <path d=\"M 250 403 L 251 420 L 259 406 L 250 403\"></path> <!-- node 2 to 3 -->\n   <path d=\"M 241 245 L 249 228 L 249 248 L 241 245\"></path> <!-- node 3 to 2 -->\n</g>\n<!-- nodes -->\n<g><ellipse fill=\"#e7d4e8\" cx=\"50\" cy=\"225\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"250\" cy=\"25\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"250\" cy=\"225\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"450\" cy=\"225\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"250\" cy=\"425\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<!-- numbers -->\n<text x=\"46\" y=\"229\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">0</text>\n<text x=\"246\" y=\"29\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">1</text>\n<text x=\"246\" y=\"229\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">2</text>\n<text x=\"246\" y=\"429\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">3</text>\n<text x=\"446\" y=\"229\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">4</text>\n</svg>","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"From Taha 'Introduction to Operations Research', example 6.4-2.\"\"\"\n\nfrom __future__ import print_function\nfrom ortools.graph import pywrapgraph\n\ndef main():\n    \"\"\"MaxFlow simple interface example.\"\"\"\n\n    # Define three parallel arrays: start_nodes, end_nodes, and the capacities\n    # between each pair. For instance, the arc from node 0 to node 1 has a\n    # capacity of 20.\n\n    start_nodes = [0, 0, 0, 1, 1, 2, 2, 3, 3]\n    end_nodes = [1, 2, 3, 2, 4, 3, 4, 2, 4]\n    capacities = [20, 30, 10, 40, 30, 10, 20, 5, 20]\n\n    # Instantiate a SimpleMaxFlow solver.\n    max_flow = pywrapgraph.SimpleMaxFlow()\n    # Add each arc.\n    for i in range(0, len(start_nodes)):\n        max_flow.AddArcWithCapacity(\n            start_nodes[i], end_nodes[i], capacities[i])\n\n    # Find the maximum flow between node 0 and node 4.\n    if max_flow.Solve(0, 4) == max_flow.OPTIMAL:\n        print('Max flow:', max_flow.OptimalFlow())\n        print('')\n        print('  Arc    Flow / Capacity')\n        for i in range(max_flow.NumArcs()):\n            print('%1s -> %1s   %3s  / %3s' % (\n                max_flow.Tail(i),\n                max_flow.Head(i),\n                max_flow.Flow(i),\n                max_flow.Capacity(i)))\n        print('Source side min-cut:', max_flow.GetSourceSideMinCut())\n        print('Sink side min-cut:', max_flow.GetSinkSideMinCut())\n    else:\n        print('There was an issue with the max flow input.')\n\nif __name__ == '__main__':\n    main()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Minimum Flows:** <br>\nMinimize the cost of using a network to achieve destination\n\n- At a supply node, a positive amount — the supply — is added to the flow. A supply could represent production at that node, for example.\n- At a demand node, a negative amount — the demand — is taken away from the flow. A demand could represent consumption at that node, for example.\n- The graph below shows a min cost flow problem. The arcs are labeled with pairs of numbers: the first number is the capacity and the second number is the cost. The numbers in parentheses next to the nodes represent supplies or demands. Node 0 is a supply node with supply 20, while nodes 3 and 4 are demand nodes, with demands -5 and -15, respectively."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"%%html\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"700\" height=\"450\">\n<defs>\n    <marker id=\"markerArrow\" markerWidth=\"13\" markerHeight=\"13\" refX=\"2\" refY=\"5\" orient=\"auto\">\n        <path d=\"M2,2 L2,8 L12,5 L2,2\" style=\"fill: #762a83;\"></path>\n    </marker>\n</defs>\n<line x1=\"50\" y1=\"225\" x2=\"250\" y2=\"75\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 0 to 1 -->\n<line x1=\"50\" y1=\"225\" x2=\"250\" y2=\"375\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 0 to 2 -->\n<line x1=\"250\" y1=\"75\" x2=\"250\" y2=\"375\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 1 to 4 -->\n<line x1=\"250\" y1=\"75\" x2=\"550\" y2=\"375\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 1 to 2 -->\n<line x1=\"250\" y1=\"75\" x2=\"550\" y2=\"75\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 1 to 3 -->\n<line x1=\"250\" y1=\"375\" x2=\"550\" y2=\"75\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 2 to 3 -->\n<line x1=\"550\" y1=\"75\" x2=\"550\" y2=\"375\" stroke=\"#762a83\" stroke-width=\"3\"></line> <!-- node 3 to 4 -->\n<g stroke=\"#762a83\" stroke-width=\"3\" fill=\"none\">\n  <path d=\"M 250 375 Q 400 405, 550 375\"></path> <!-- node 3 to 4 -->\n  <path d=\"M 550 375 Q 400 345, 250 375\"></path> <!-- node 4 to 3 -->\n</g>\n<!-- weights -->\n<text x=\"90\" y=\"135\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(15, 4)</text>\n<text x=\"175\" y=\"215\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(20, 2)</text>\n<text x=\"455\" y=\"270\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(10, 6)</text>\n<text x=\"370\" y=\"65\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(4, 2)</text>\n<text x=\"280\" y=\"270\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(15, 1)</text>\n<text x=\"90\" y=\"315\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(8, 4)</text>\n<text x=\"370\" y=\"350\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(5, 3)</text>\n<text x=\"370\" y=\"410\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(4, 3)</text>\n<text x=\"560\" y=\"215\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(20, 2)</text>\n<!-- supplies -->\n<text x=\"0\" y=\"228\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(20)</text>\n<text x=\"236\" y=\"58\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(0)</text>\n<text x=\"236\" y=\"402\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(0)</text>\n<text x=\"560\" y=\"78\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(-5)</text>\n<text x=\"560\" y=\"378\" font-family=\"Verdana\" font-size=\"18\" stroke=\"#762a83\">(-15)</text>\n<!-- arrows -->\n<g stroke=\"#762a83\" stroke-width=\"2\" fill=\"#762a83\">\n   <path d=\"M 50 225 L 230 90\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 0 to 1 -->\n   <path d=\"M 50 225 L 230 360\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 0 to 2 -->\n   <path d=\"M 250 75 L 525 75\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 1 to 3 -->\n   <path d=\"M 250 75 L 533 358\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 1 to 4 -->\n   <path d=\"M 250 375 L 533 92\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 2 to 3 -->\n   <path d=\"M 250 75 L 250 350\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 1 to 2 -->\n   <path d=\"M 550 75 L 550 350\" style=\"marker-end: url(#markerArrow);\"></path> <!-- node 3 to 4 -->\n   <path d=\"M 250 375 Q 400 405, 525 380\" fill=\"none\" style=\"marker-end: url(#markerArrow);\"></path>\n   <!-- node 2 to 4 -->\n   <path d=\"M 550 375 Q 400 345, 275 370\" fill=\"none\" style=\"marker-end: url(#markerArrow);\"></path>\n   <!-- node 4 to 2 -->\n</g>\n<!-- nodes -->\n<g><ellipse fill=\"#e7d4e8\" cx=\"50\" cy=\"225\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"250\" cy=\"75\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"250\" cy=\"375\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"550\" cy=\"75\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<g><ellipse fill=\"#e7d4e8\" cx=\"550\" cy=\"375\" rx=\"8\" ry=\"8\" stroke=\"#762a83\"></ellipse></g>\n<!-- numbers -->\n<text x=\"46\" y=\"229\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">0</text>\n<text x=\"246\" y=\"79\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">1</text>\n<text x=\"246\" y=\"379\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">2</text>\n<text x=\"546\" y=\"79\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">3</text>\n<text x=\"546\" y=\"379\" font-family=\"Verdana\" font-size=\"12\" stroke=\"#762a83\">4</text>\n</svg>","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# \"\"\"From Bradley, Hax, and Magnanti, 'Applied Mathematical Programming', figure 8.1.\"\"\"\n\nfrom __future__ import print_function\nfrom ortools.graph import pywrapgraph\n\n\ndef main():\n    \"\"\"MinCostFlow simple interface example.\"\"\"\n\n    # Define four parallel arrays: start_nodes, end_nodes, capacities, and unit costs\n    # between each pair. For instance, the arc from node 0 to node 1 has a\n    # capacity of 15 and a unit cost of 4.\n\n    start_nodes = [0, 0,  1, 1,  1,  2, 2,  3, 4]\n    end_nodes = [1, 2,  2, 3,  4,  3, 4,  4, 2]\n    capacities = [15, 8, 20, 4, 10, 15, 4, 20, 5]\n    unit_costs = [4, 4,  2, 2,  6,  1, 3,  2, 3]\n\n    # Define an array of supplies at each node.\n\n    supplies = [20, 0, 0, -5, -15]\n\n    # Instantiate a SimpleMinCostFlow solver.\n    min_cost_flow = pywrapgraph.SimpleMinCostFlow()\n\n    # Add each arc.\n    for i in range(0, len(start_nodes)):\n        min_cost_flow.AddArcWithCapacityAndUnitCost(start_nodes[i], end_nodes[i],\n                                                    capacities[i], unit_costs[i])\n\n    # Add node supplies.\n\n    for i in range(0, len(supplies)):\n        min_cost_flow.SetNodeSupply(i, supplies[i])\n\n    # Find the minimum cost flow between node 0 and node 4.\n    if min_cost_flow.Solve() == min_cost_flow.OPTIMAL:\n        print('Minimum cost:', min_cost_flow.OptimalCost())\n        print('')\n        print('  Arc    Flow / Capacity  Cost')\n        for i in range(min_cost_flow.NumArcs()):\n            cost = min_cost_flow.Flow(i) * min_cost_flow.UnitCost(i)\n            print('%1s -> %1s   %3s  / %3s       %3s' % (\n                min_cost_flow.Tail(i),\n                min_cost_flow.Head(i),\n                min_cost_flow.Flow(i),\n                min_cost_flow.Capacity(i),\n                cost))\n    else:\n        print('There was an issue with the min cost flow input.')\n\n\nif __name__ == '__main__':\n    main()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Notebook Runtime: %0.2f Minutes\"%((time.time() - start)/60))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}