Linear_programming_optimization_with_Python_final.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "68a974bd",
   "metadata": {},
   "source": [
    "# Linear Programming Optimization With Python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97b7eac4",
   "metadata": {},
   "source": [
    "Linear programming is a mathematical optimization technique used to optimize a linear objective function subject to a set of linear constraints. It involves finding the values of decision variables that minimize or maximize the objective function while satisfying the given constraints.\n",
    "\n",
    "The objective function and constraints are represented as linear equations or inequalities, and the decision variables are the unknowns that need to be determined. The goal of linear programming is to find the optimal solution that satisfies all the constraints while optimizing the objective function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c467a40d",
   "metadata": {},
   "source": [
    "Pulp, Pyomo, Scipy, OR-tools are among popular open-source modeling languages used for formulating and solving linear programming problems. They provide powerful tools for modeling and solving linear programming problems. They offer flexible and intuitive syntaxes for specifying optimization models and can be easily integrated into existing Python-based workflows. Additionally, they both support a wide range of solvers (glpk, cbc, gurobi, couenne and so on) making it easy to find the best solver for a particular problem.\n",
    "\n",
    "This notebook contains 8 mini problems from various optimization topics. These problems are purposely designed to provide practice in formulating optimization problems and becoming familiar with popular optimization libraries in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf4f5b7",
   "metadata": {},
   "source": [
    "**1. Problem 1**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ac45466e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GLPK_CMD']"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test available solvers that pulp can use in this machine\n",
    "import pulp\n",
    "pulp.listSolvers(onlyAvailable=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f4d7d9",
   "metadata": {},
   "source": [
    "**Using pyomo with glpk solver**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fde9ed74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value: 960.0\n",
      "Optimal solution: \n",
      "x1 = 1200.0\n",
      "x2 = 0.0\n",
      "x3 = 0.0\n",
      "x4 = 0.0\n"
     ]
    }
   ],
   "source": [
    "# a) Find manufacturing plan to optimize revenue\n",
    "\n",
    "from pyomo.environ import *\n",
    "\n",
    "# Create the model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# Define the decision variables\n",
    "for i in range(1,5):\n",
    "    model.add_component(f\"x{i}\", Var(within=NonNegativeReals))\n",
    "\n",
    "# Define the objective function\n",
    "model.obj = Objective(expr=0.8*model.x1 + 0.3*model.x2 + 0.38*model.x3 + 0.4*model.x4, sense=maximize)\n",
    "\n",
    "# Define the constraints\n",
    "model.con1 = Constraint(expr=0.5*model.x1 + 0.2*model.x2 + 0.3*model.x3 + 0.6*model.x4 <= 600)\n",
    "model.con2 = Constraint(expr=0.1*model.x1 + 0.4*model.x2 + 0.2*model.x3 + 0.5*model.x4 <= 800)\n",
    "\n",
    "# Solve the linear programming problem\n",
    "solver = SolverFactory('glpk', validate=False)\n",
    "results = solver.solve(model)\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Optimal value: {model.obj()}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {getattr(model,f'x{i}')()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3740c554",
   "metadata": {},
   "source": [
    "**Using pyomo with cbc solver**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6331b16e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value: 966.666665\n",
      "Optimal solution: \n",
      "x1 = 1000.0\n",
      "x2 = 0.0\n",
      "x3 = 333.33333\n",
      "x4 = 0.0\n"
     ]
    }
   ],
   "source": [
    "# b) With new constraints and using cbc solver\n",
    "\n",
    "# Create the model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# Define the decision variables\n",
    "for i in range(1,5):\n",
    "    model.add_component(f\"x{i}\", Var(within=NonNegativeReals))\n",
    "\n",
    "# Define the objective function\n",
    "model.obj = Objective(expr=0.8*model.x1 + 0.3*model.x2 + 0.5*model.x3 + 0.4*model.x4, sense=maximize)\n",
    "\n",
    "# Define the constraints\n",
    "model.con1 = Constraint(expr=0.5*model.x1 + 0.2*model.x2 + 0.3*model.x3 + 0.6*model.x4 <= 600) # 600 kg\n",
    "model.con2 = Constraint(expr=0.1*model.x1 + 0.4*model.x2 + 0.2*model.x3 + 0.5*model.x4 <= 800) # 800 kg\n",
    "model.con3 = Constraint(expr=model.x1 + model.x2 >= 1000) \n",
    "\n",
    "# Solve the linear programming problem\n",
    "solver = SolverFactory('cbc', executable='D:/PROJECTS/Data analytics/Solvers/bin/cbc.exe')\n",
    "results = solver.solve(model)\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Optimal value: {model.obj()}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {getattr(model,f'x{i}')()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55f8e42d",
   "metadata": {},
   "source": [
    "**2. Problem 2**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d9097e5",
   "metadata": {},
   "source": [
    "**Using pulp to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9b9dd62e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value: 46.5\n",
      "Optimal solution: \n",
      "x1 = 100.0\n",
      "x2 = 100.0\n",
      "x3 = 150.0\n",
      "x4 = 150.0\n"
     ]
    }
   ],
   "source": [
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_2\", pulp.LpMaximize)\n",
    "\n",
    "# Define the decision variables\n",
    "x1, x2, x3, x4 = [pulp.LpVariable(f\"x{i}\", lowBound=0) for i in range(4)]\n",
    "\n",
    "# Define the objective function\n",
    "prob += 0.07 * x1 + 0.08 * x2 + 0.1 * x3 + 0.11 * x4\n",
    "\n",
    "resource = 500 # unit in million\n",
    "\n",
    "# Define the constraints\n",
    "prob += x1 <= 100\n",
    "prob += x2 <= 100\n",
    "prob += x3 <= 150\n",
    "prob += x4 <= 150\n",
    "prob += x2 + x3 >= 0.5 * resource\n",
    "prob += x1 >= 0.1 * resource\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {eval(f'x{i}.varValue')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3d043e",
   "metadata": {},
   "source": [
    "**Using pyomo with gurobi solver**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5bda80c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value: 46.5\n",
      "Optimal solution: \n",
      "x1 = 100.0\n",
      "x2 = 100.0\n",
      "x3 = 150.0\n",
      "x4 = 150.0\n"
     ]
    }
   ],
   "source": [
    "from pyomo.environ import *\n",
    "\n",
    "# Create the pulp model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# Define the decision variables\n",
    "for i in range(1,5):\n",
    "    model.add_component(f\"x{i}\", Var(within=NonNegativeReals))\n",
    "\n",
    "# Define the objective function\n",
    "model.obj = Objective(expr=0.07*model.x1 + 0.08*model.x2 + 0.1*model.x3 + 0.11*model.x4, sense=maximize)\n",
    "\n",
    "resource = 500 # unit in million\n",
    "\n",
    "# Define the constraints\n",
    "model.con1 = Constraint(expr=model.x1 <= 100) \n",
    "model.con2 = Constraint(expr=model.x2 <= 100) \n",
    "model.con3 = Constraint(expr=model.x3 <= 150) \n",
    "model.con4 = Constraint(expr=model.x4 <= 150) \n",
    "model.con5 = Constraint(expr=model.x1 >= 0.1 * resource) \n",
    "model.con6 = Constraint(expr=model.x2 + model.x3 >= 0.5 * resource) \n",
    "\n",
    "# Solve the linear programming problem using gurobi solver\n",
    "solver = SolverFactory('gurobi', executable='D:/PROJECTS/Data analytics/Solvers/bin/gurobi/bin/gurobi.bat')\n",
    "results = solver.solve(model)\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Optimal value: {model.obj()}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {getattr(model,f'x{i}')()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3c1ec25",
   "metadata": {},
   "source": [
    "**3. Problem 3**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c22bafe1",
   "metadata": {},
   "source": [
    "**Using pulp to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2b4efc7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value (USD): 120500\n",
      "Optimal solution: \n",
      "x1 = 5\n",
      "x2 = 5\n",
      "x3 = 24\n",
      "x4 = 16\n"
     ]
    }
   ],
   "source": [
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_3\", pulp.LpMaximize)\n",
    "\n",
    "# Define the decision variables\n",
    "x1, x2, x3, x4 = [pulp.LpVariable(f\"x{i}\", lowBound=0, cat=\"Integer\") for i in range(4)]\n",
    "\n",
    "# Define the objective function\n",
    "prob += 1000 * x1 + 1500 * x2 + 2500 * x3 + 3000 * x4\n",
    "\n",
    "# Define the constraints\n",
    "prob += x1 >= 5\n",
    "prob += x2 >= 5\n",
    "prob += x3 >= 5\n",
    "prob += x4 >= 5\n",
    "prob += x1 + x2 + x3 + x4 == 50\n",
    "prob += 9 * x1 + 25 * x2 + 25 * x3 + 43 * x4 <= 1460\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value (USD): {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {eval(f'x{i}.varValue')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5daf2a03",
   "metadata": {},
   "source": [
    "**4. Problem 4**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70bbe208",
   "metadata": {},
   "source": [
    "**Using pulp to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "551740f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Infeasible\n",
      "Optimal value: 0\n",
      "Optimal solution: \n",
      "x1 = 0\n",
      "x2 = 0\n",
      "x3 = 0\n",
      "x4 = 0\n"
     ]
    }
   ],
   "source": [
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_4\", pulp.LpMaximize)\n",
    "\n",
    "# Define the decision variables\n",
    "# M=male, F=female, C=cut, M=mill\n",
    "# 1 = MC, 2 = MM, 3 = FC, 4 = FM\n",
    "x1, x2, x3, x4 = [pulp.LpVariable(f\"x{i}\", lowBound=0, cat=\"Integer\") for i in range(4)]\n",
    "\n",
    "# Define the objective function\n",
    "prob += 30*x1 + 25*x2 + 28*x3 + 20*x4\n",
    "\n",
    "# Define the constraints\n",
    "prob += x1 + x2 == 32\n",
    "prob += x3 + x4 == 20\n",
    "prob += 30*x1 + 28*x3 == 25*x2 + 20*x4\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(1,5):\n",
    "    print(f\"x{i} = {eval(f'x{i}.varValue')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b569a6da",
   "metadata": {},
   "source": [
    "The equality constraints are too difficult to meet, so there is no integer point that satisfies all the constraints as defined.\n",
    "\n",
    "Representing x2 in terms of x1, x4 in terms of x3, and substituting into the third equation:\n",
    "\n",
    "30x1 + 28x3 = 25*(32-x1) + 20*(20 - x3)\n",
    "(30+25)x1 + (28+20)x3 = 2532 + 2020\n",
    "\n",
    "If the variables are not required to be integers, then the system of constraint equations can be met.\n",
    "\n",
    "However, when we require x1, x2, x3, and x4 to be integers, there is no feasible point. Geometrically, this means that the line (30+25)x1 + (28+20)x3 = 2532 + 2020 does not intersect any point where both x1 and x3 are integers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fad6359",
   "metadata": {},
   "source": [
    "**5. Problem 5**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed652168",
   "metadata": {},
   "source": [
    "**Using or-tools to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4e813d92",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal objective value: 68581081.08108108\n",
      "xA: 67.56756756756755\n",
      "xB: 0.0\n",
      "xC: 0.0\n",
      "xD: 60.81081081081081\n"
     ]
    }
   ],
   "source": [
    "from ortools.linear_solver import pywraplp\n",
    "\n",
    "# Create the linear solver\n",
    "solver = pywraplp.Solver.CreateSolver('GLOP')\n",
    "\n",
    "# Define the decision variables\n",
    "names = ['xA', 'xB', 'xC', 'xD']\n",
    "variables = [solver.IntVar(0, solver.infinity(), name) for name in names]\n",
    "\n",
    "# Define the objective function\n",
    "coefficients = [250000, 350000, 380000, 850000]\n",
    "obj = solver.Objective()\n",
    "\n",
    "for var, coeff in zip(variables, coefficients):\n",
    "    obj.SetCoefficient(var, coeff)\n",
    "obj.SetMaximization()\n",
    "\n",
    "# Define the constraints\n",
    "constraints = [\n",
    "    {'coefficients': [0.08, 0.12, 0.3, 0.21], 'upper_bound': 350},\n",
    "    {'coefficients': [4, 9, 7, 12], 'upper_bound': 1000},\n",
    "    {'coefficients': [1, 1, 3, 40], 'upper_bound': 2500}\n",
    "]\n",
    "for cons_data in constraints:\n",
    "    constraint = solver.Constraint(-solver.infinity(), cons_data['upper_bound'])\n",
    "    for var, coeff in zip(variables, cons_data['coefficients']):\n",
    "        constraint.SetCoefficient(var, coeff)\n",
    "\n",
    "# Solve the LP problem\n",
    "solver.Solve()\n",
    "\n",
    "# Get the optimal objective value\n",
    "print('Optimal objective value:', solver.Objective().Value())\n",
    "\n",
    "# Get the optimal decision variables\n",
    "for variable in variables:\n",
    "    print(variable.name() + ':', variable.solution_value())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e087cc2c",
   "metadata": {},
   "source": [
    "**6. Problem 6**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b448e0fb",
   "metadata": {},
   "source": [
    "**Using pulp to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "115c5959",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value: 435\n",
      "Optimal solution: \n",
      "x[0][0] = 0\n",
      "x[0][1] = 5\n",
      "x[0][2] = 0\n",
      "x[0][3] = 0\n",
      "x[0][4] = 0\n",
      "x[1][0] = 10\n",
      "x[1][1] = 0\n",
      "x[1][2] = 0\n",
      "x[1][3] = 0\n",
      "x[1][4] = 5\n",
      "x[2][0] = 0\n",
      "x[2][1] = 0\n",
      "x[2][2] = 10\n",
      "x[2][3] = 0\n",
      "x[2][4] = 10\n",
      "x[3][0] = 0\n",
      "x[3][1] = 5\n",
      "x[3][2] = 0\n",
      "x[3][3] = 20\n",
      "x[3][4] = 5\n"
     ]
    }
   ],
   "source": [
    "# a) Problem as described\n",
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_6a\", pulp.LpMinimize)\n",
    "\n",
    "# Define the decision variables\n",
    "factory = [5, 15, 20, 30]\n",
    "warehouse6a = [10, 10, 10, 20, 20]\n",
    "distance = [[5, 1, 4, 6, 7], [3, 4, 2, 7, 8], [4, 3, 1, 7, 9], [6, 5, 4, 9, 11]]\n",
    "\n",
    "dim_i = len(factory)\n",
    "dim_j = len(warehouse6a)\n",
    "\n",
    "# Amount of goods in tons\n",
    "x = [[pulp.LpVariable(\"x_{}_{}\".format(i, j), lowBound=0, cat=\"Integer\") for j in range(dim_j)] for i in range(dim_i)]\n",
    "\n",
    "# Define the objective function as sum across rows and columns\n",
    "prob += pulp.lpSum([[distance[i][j] * x[i][j] for i in range(dim_i)] for j in range(dim_j)])\n",
    "\n",
    "# Define 4 constraints for 4 factories\n",
    "for i in range(dim_i):\n",
    "    prob += pulp.lpSum([x[i][j] for j in range(dim_j)]) == factory[i]\n",
    "# and 5 contraints for warehouses\n",
    "for j in range(dim_j):\n",
    "    prob += pulp.lpSum([x[i][j] for i in range(dim_i)]) == warehouse6a[j]\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(dim_i):\n",
    "    for j in range(dim_j):\n",
    "        print(f\"x[{i}][{j}] = {x[i][j].varValue}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "725698ae",
   "metadata": {},
   "source": [
    "**Using pyomo with couenne solver**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "203230b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: ok\n",
      "Optimal value: 435.0\n",
      "Optimal solution: \n",
      "x[0][0] = 0.0\n",
      "x[0][1] = 3.0\n",
      "x[0][2] = 0.0\n",
      "x[0][3] = 0.0\n",
      "x[0][4] = 2.0\n",
      "x[1][0] = 6.0\n",
      "x[1][1] = 0.0\n",
      "x[1][2] = 0.0\n",
      "x[1][3] = 0.0\n",
      "x[1][4] = 9.0\n",
      "x[2][0] = 2.0\n",
      "x[2][1] = 3.0\n",
      "x[2][2] = 10.0\n",
      "x[2][3] = 5.0\n",
      "x[2][4] = 0.0\n",
      "x[3][0] = 2.0\n",
      "x[3][1] = 4.0\n",
      "x[3][2] = 0.0\n",
      "x[3][3] = 15.0\n",
      "x[3][4] = 9.0\n"
     ]
    }
   ],
   "source": [
    "from pyomo.environ import *\n",
    "\n",
    "# Create the model\n",
    "model = ConcreteModel()\n",
    "\n",
    "# Define the decision variables\n",
    "factory = [5, 15, 20, 30]\n",
    "warehouse6a = [10, 10, 10, 20, 20]\n",
    "distance = [[5, 1, 4, 6, 7], [3, 4, 2, 7, 8], [4, 3, 1, 7, 9], [6, 5, 4, 9, 11]]\n",
    "\n",
    "dim_i = len(factory)\n",
    "dim_j = len(warehouse6a)\n",
    "\n",
    "# Amount of goods in tons\n",
    "model.x = Var(range(dim_i), range(dim_j), within=NonNegativeIntegers)\n",
    "\n",
    "# Define the objective function as sum across rows and columns\n",
    "model.obj = Objective(expr=sum(distance[i][j] * model.x[i, j] for i in range(dim_i) for j in range(dim_j)), sense=minimize)\n",
    "\n",
    "# Define 4 constraints for 4 factories\n",
    "def factory_constraint(model, i):\n",
    "    return sum(model.x[i, j] for j in range(dim_j)) == factory[i]\n",
    "model.factory_constr = Constraint(range(dim_i), rule=factory_constraint)\n",
    "\n",
    "# and 5 contraints for warehouses\n",
    "def warehouse_constraint(model, j):\n",
    "    return sum(model.x[i, j] for i in range(dim_i)) == warehouse6a[j]\n",
    "model.warehouse_constr = Constraint(range(dim_j), rule=warehouse_constraint)\n",
    "\n",
    "# Solve the LP problem using couenne solver\n",
    "solver = SolverFactory('couenne', executable='D:/PROJECTS/Data analytics/Solvers/bin/couenne.exe')\n",
    "results = solver.solve(model)\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {results.solver.status}\")\n",
    "print(f\"Optimal value: {model.obj.expr()}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(dim_i):\n",
    "    for j in range(dim_j):\n",
    "        print(f\"x[{i}][{j}] = {model.x[i, j].value}\") \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "78a2c18b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value: 405\n",
      "Optimal solution: \n",
      "x[0][0] = 0\n",
      "x[0][1] = 5\n",
      "x[0][2] = 0\n",
      "x[0][3] = 0\n",
      "x[0][4] = 0\n",
      "x[1][0] = 10\n",
      "x[1][1] = 0\n",
      "x[1][2] = 0\n",
      "x[1][3] = 0\n",
      "x[1][4] = 5\n",
      "x[2][0] = 0\n",
      "x[2][1] = 0\n",
      "x[2][2] = 15\n",
      "x[2][3] = 0\n",
      "x[2][4] = 5\n",
      "x[3][0] = 0\n",
      "x[3][1] = 5\n",
      "x[3][2] = 0\n",
      "x[3][3] = 15\n",
      "x[3][4] = 10\n"
     ]
    }
   ],
   "source": [
    "# b) Probelem as described but with C3 = C4 = 15 tons\n",
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_6b\", pulp.LpMinimize)\n",
    "\n",
    "# Define the decision variables\n",
    "factory = [5, 15, 20, 30]\n",
    "warehouse6b = [10, 10, 15, 15, 20]\n",
    "distance = [[5, 1, 4, 6, 7], [3, 4, 2, 7, 8], [4, 3, 1, 7, 9], [6, 5, 4, 9, 11]]\n",
    "\n",
    "dim_i = len(factory)\n",
    "dim_j = len(warehouse6b)\n",
    "\n",
    "# Amount of goods in tons\n",
    "x = [[pulp.LpVariable(\"x_{}_{}\".format(i, j), lowBound=0, cat=\"Integer\") for j in range(dim_j)] for i in range(dim_i)]\n",
    "\n",
    "## Define the objective function\n",
    "prob += pulp.lpSum([[distance[i][j] * x[i][j] for i in range(dim_i)] for j in range(dim_j)])\n",
    "\n",
    "## Define the constraints\n",
    "# Define 4 constraints for 4 factories\n",
    "for i in range(dim_i):\n",
    "    prob += pulp.lpSum([x[i][j] for j in range(dim_j)]) == factory[i]\n",
    "# and 5 contraints for warehouses\n",
    "for j in range(dim_j):\n",
    "    prob += pulp.lpSum([x[i][j] for i in range(dim_i)]) == warehouse6b[j]\n",
    "    \n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(dim_i):\n",
    "    for j in range(dim_j):\n",
    "        print(f\"x[{i}][{j}] = {x[i][j].varValue}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28d2d7a2",
   "metadata": {},
   "source": [
    "**7. Problem 7**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d3a4dc6",
   "metadata": {},
   "source": [
    "**Using Scipy to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1c75d454",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimization terminated successfully.\n",
      "Optimal value: 51.0\n",
      "Optimal solution: FI = 1.0, FII = 1.0, FIII = 0.0, MI = 0.0, MII = 0.0, MIII = 3.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<frozen importlib._bootstrap>:228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n",
      "C:\\Users\\huyha\\AppData\\Local\\Temp\\ipykernel_11140\\1732714596.py:26: DeprecationWarning: `method='simplex'` is deprecated and will be removed in SciPy 1.11.0. Please use one of the HiGHS solvers (e.g. `method='highs'`) in new code.\n",
      "  results = linprog(obj_fun, A_ub=lhs_eq, b_ub=rhs_eq, bounds=bounds, method=\"simplex\")\n"
     ]
    }
   ],
   "source": [
    "# 7) Problem 7\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "\n",
    "# Define the objective function\n",
    "obj_fun = [-10, -8, -7, -8, -9, -11]\n",
    "\n",
    "# Define the constraints\n",
    "# Worker: M=male, F=female\n",
    "# Milling_machine: I, II, II\n",
    "\n",
    "lhs_eq = [\n",
    "    [1, 1, 1, 0, 0, 0],\n",
    "    [0, 0, 0, 1, 1, 1],\n",
    "    [1, 0, 0, 1, 0, 0],\n",
    "    [0, 1, 0, 0, 1, 0],\n",
    "    [0, 0, 1, 0, 0, 1]\n",
    "]\n",
    "rhs_eq = [2, 3, 1, 2, 3]\n",
    "\n",
    "# Define the bounds on the variables\n",
    "bounds = [(0, None)] * 6\n",
    "\n",
    "# Solve the LP problem\n",
    "results = linprog(obj_fun, A_ub=lhs_eq, b_ub=rhs_eq, bounds=bounds, method=\"simplex\")\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {results.message}\")\n",
    "print(f\"Optimal value: {-results.fun}\")\n",
    "roundRes = np.around(results.x)\n",
    "print(f\"Optimal solution: FI = {roundRes[0]}, FII = {roundRes[1]}, FIII = {roundRes[2]}, MI = {roundRes[3]}, MII = {roundRes[4]}, MIII = {roundRes[5]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "443dd5df",
   "metadata": {},
   "source": [
    "**8. Problem 8**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0381870d",
   "metadata": {},
   "source": [
    "**Using pulp to solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c39a7fa1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value: 585\n",
      "Optimal solution: \n",
      "x[0][0] = 0\n",
      "x[0][1] = 25\n",
      "x[0][2] = 5\n",
      "x[1][0] = 0\n",
      "x[1][1] = 0\n",
      "x[1][2] = 20\n",
      "x[2][0] = 25\n",
      "x[2][1] = 0\n",
      "x[2][2] = 15\n"
     ]
    }
   ],
   "source": [
    "# a) Solve the problem as described\n",
    "\n",
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_8a\", pulp.LpMaximize)\n",
    "\n",
    "# Define the decision variables\n",
    "rice_type = [30, 20, 40]\n",
    "field_type = [25, 25, 40]\n",
    "productivity8a = [[4, 8, 6], [6, 9, 7], [5, 4, 6]] #tons/hectare\n",
    "\n",
    "dim_i = len(rice_type)\n",
    "dim_j = len(field_type)\n",
    "\n",
    "# Productivity per hectares rice field\n",
    "x = [[pulp.LpVariable(\"x_{}_{}\".format(i, j), lowBound=0, cat=\"Integer\") for j in range(dim_j)] for i in range(dim_i)]\n",
    "\n",
    "# Define the objective function\n",
    "prob += pulp.lpSum([[productivity8a[i][j] * x[i][j] for i in range(dim_i)] for j in range(dim_j)])\n",
    "\n",
    "# Define the constraints\n",
    "for i in range(dim_i):\n",
    "    prob += pulp.lpSum([x[i][j] for j in range(dim_j)]) == rice_type[i]\n",
    "    \n",
    "for j in range(dim_j):\n",
    "    prob += pulp.lpSum([x[i][j] for i in range(dim_i)]) == field_type[j]\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(dim_i):\n",
    "    for j in range(dim_j):\n",
    "        print(f\"x[{i}][{j}] = {x[i][j].varValue}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "854a5cfc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Optimal value: 585\n",
      "Optimal solution: \n",
      "x[0][0] = 0\n",
      "x[0][1] = 25\n",
      "x[0][2] = 5\n",
      "x[1][0] = 0\n",
      "x[1][1] = 0\n",
      "x[1][2] = 20\n",
      "x[2][0] = 25\n",
      "x[2][1] = 0\n",
      "x[2][2] = 15\n"
     ]
    }
   ],
   "source": [
    "# b) Formulate and solve the problem in the case where IR8 rice cannot be grown on type I field\n",
    "\n",
    "import pulp\n",
    "\n",
    "# Create the LP problem\n",
    "prob = pulp.LpProblem(\"LP_Problem_8b\", pulp.LpMaximize)\n",
    "\n",
    "# Define the decision variables\n",
    "rice_type = [30, 20, 40]\n",
    "field_type = [25, 25, 40]\n",
    "productivity8b = [[0, 8, 6], [6, 9, 7], [5, 4, 6]] #tons/hectare\n",
    "\n",
    "dim_i = len(rice_type)\n",
    "dim_j = len(field_type)\n",
    "\n",
    "# Productivity per hectares rice field\n",
    "x = [[pulp.LpVariable(\"x_{}_{}\".format(i, j), lowBound=0, cat=\"Integer\") for j in range(dim_j)] for i in range(dim_i)]\n",
    "\n",
    "# Define the objective function\n",
    "prob += pulp.lpSum([[productivity8b[i][j] * x[i][j] for i in range(dim_i)] for j in range(dim_j)])\n",
    "\n",
    "\n",
    "# Define the constraints\n",
    "for i in range(dim_i):\n",
    "    prob += pulp.lpSum([x[i][j] for j in range(dim_j)]) == rice_type[i]\n",
    "    \n",
    "for j in range(dim_j):\n",
    "    prob += pulp.lpSum([x[i][j] for i in range(dim_i)]) == field_type[j]\n",
    "\n",
    "# Solve the LP problem\n",
    "status = prob.solve()\n",
    "\n",
    "# Print the solution\n",
    "print(f\"Status: {pulp.LpStatus[status]}\")\n",
    "print(f\"Optimal value: {pulp.value(prob.objective)}\")\n",
    "print(\"Optimal solution: \")\n",
    "for i in range(dim_i):\n",
    "    for j in range(dim_j):\n",
    "        print(f\"x[{i}][{j}] = {x[i][j].varValue}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac314398",
   "metadata": {},
   "source": [
    "The condition of IR8 rice cannot be grown on type I field does not affect the total harvest."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}