MILP_optimization_with_Python_P8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "68a974bd",
   "metadata": {},
   "source": [
    "# Linear Programming and MILP Optimization With Python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c467a40d",
   "metadata": {},
   "source": [
    "This notebook demonstrates the use of Mixed Integer Linear Programming (MILP) to solve this optimization problem. The goal is to allocate different rice varieties to different field types to maximize the total harvest yield. Pyomo and OR-tools are used together with the linear and PLEX solvers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf4f5b7",
   "metadata": {},
   "source": [
    "**1. Using Google OR-tools solver**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64840e8e-6f9e-42f9-aad3-d53875e646ec",
   "metadata": {},
   "source": [
    "**a) Formulate and solve the optimization problem.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ac45466e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up OR-tools library\n",
    "from ortools.linear_solver import pywraplp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "17101f7a-c6a3-4b86-b33a-d87e422de550",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the MILP solver with the SCIP backend.\n",
    "solver = pywraplp.Solver.CreateSolver('SCIP')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cddb6fe7-8966-42aa-b300-81f080b0062f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define decision variables\n",
    "fields = ['I', 'II', 'III']\n",
    "rice_types = ['IR8', 'CBC', 'IR132']\n",
    "x = {}\n",
    "for f in fields:\n",
    "    for r in rice_types:\n",
    "        x[(f, r)] = solver.IntVar(0, solver.infinity(), f'x_{f}_{r}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fde9ed74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ortools.linear_solver.pywraplp.Constraint; proxy of <Swig Object of type 'operations_research::MPConstraint *' at 0x000001AC70485890> >"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Objective function\n",
    "yield_matrix = {'I': {'IR8': 4, 'CBC': 6, 'IR132': 5},\n",
    "                'II': {'IR8': 8, 'CBC': 9, 'IR132': 4},\n",
    "                'III': {'IR8': 6, 'CBC': 7, 'IR132': 6}}\n",
    "\n",
    "solver.Maximize(solver.Sum(x[(f, r)] * yield_matrix[f][r] for f in fields for r in rice_types))\n",
    "\n",
    "# Constraints for total hectares for each rice type\n",
    "solver.Add(solver.Sum(x[(f, 'IR8')] for f in fields) == 30)\n",
    "solver.Add(solver.Sum(x[(f, 'CBC')] for f in fields) == 20)\n",
    "solver.Add(solver.Sum(x[(f, 'IR132')] for f in fields) == 40)\n",
    "\n",
    "# Constraints for total hectares for each field type\n",
    "solver.Add(solver.Sum(x[('I', r)] for r in rice_types) == 25)\n",
    "solver.Add(solver.Sum(x[('II', r)] for r in rice_types) == 25)\n",
    "solver.Add(solver.Sum(x[('III', r)] for r in rice_types) == 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "90c6512e-827b-47d0-8b14-5896daa9916a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Field I, Rice IR8: 0.0 hectares\n",
      "Field I, Rice CBC: 20.0 hectares\n",
      "Field I, Rice IR132: 5.0 hectares\n",
      "Field II, Rice IR8: 25.0 hectares\n",
      "Field II, Rice CBC: 0.0 hectares\n",
      "Field II, Rice IR132: 0.0 hectares\n",
      "Field III, Rice IR8: 5.0 hectares\n",
      "Field III, Rice CBC: 0.0 hectares\n",
      "Field III, Rice IR132: 35.0 hectares\n"
     ]
    }
   ],
   "source": [
    "# Solve the problem\n",
    "status = solver.Solve()\n",
    "\n",
    "# Output results\n",
    "if status == pywraplp.Solver.OPTIMAL:\n",
    "    print(\"Status:\", \"Optimal\")\n",
    "    for f in fields:\n",
    "        for r in rice_types:\n",
    "            print(f\"Field {f}, Rice {r}: {x[(f, r)].solution_value()} hectares\")\n",
    "else:\n",
    "    print(\"The problem does not have an optimal solution.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3740c554",
   "metadata": {},
   "source": [
    "**b) Additional scenario where IR8 cannot be grown on Field Type I**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6331b16e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n",
      "Field I, Rice IR8: 0.0 hectares\n",
      "Field I, Rice CBC: 20.0 hectares\n",
      "Field I, Rice IR132: 5.0 hectares\n",
      "Field II, Rice IR8: 25.0 hectares\n",
      "Field II, Rice CBC: 0.0 hectares\n",
      "Field II, Rice IR132: 0.0 hectares\n",
      "Field III, Rice IR8: 5.0 hectares\n",
      "Field III, Rice CBC: 0.0 hectares\n",
      "Field III, Rice IR132: 35.0 hectares\n"
     ]
    }
   ],
   "source": [
    "solver.Add(x[('I', 'IR8')] == 0)\n",
    "\n",
    "# Solve the problem\n",
    "status = solver.Solve()\n",
    "\n",
    "# Output results\n",
    "if status == pywraplp.Solver.OPTIMAL:\n",
    "    print(\"Status:\", \"Optimal\")\n",
    "    for f in fields:\n",
    "        for r in rice_types:\n",
    "            print(f\"Field {f}, Rice {r}: {x[(f, r)].solution_value()} hectares\")\n",
    "else:\n",
    "    print(\"The problem does not have an optimal solution.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55f8e42d",
   "metadata": {},
   "source": [
    "**2. Using CPLEX solver**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c36c9613-c955-4478-9155-8337a053cc3b",
   "metadata": {},
   "source": [
    "**a) Formulate and solve the optimization problem.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fa51182e-85eb-4b8f-b6fe-1e01dac0af6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import ConcreteModel, Var, Objective, Constraint, NonNegativeIntegers, SolverFactory, maximize\n",
    "\n",
    "# Create a new model\n",
    "model = ConcreteModel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4db6da82-cba3-4d80-a569-0d6c3239fcb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sets\n",
    "fields = ['I', 'II', 'III']\n",
    "rice_types = ['IR8', 'CBC', 'IR132']\n",
    "\n",
    "# Decision variables\n",
    "model.x = Var(fields, rice_types, domain=NonNegativeIntegers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a6ebf59a-36ab-431a-8ab4-9b7f39c0067a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective: Maximize yield\n",
    "yield_coefficients = {\n",
    "    ('I', 'IR8'): 4, ('I', 'CBC'): 6, ('I', 'IR132'): 5,\n",
    "    ('II', 'IR8'): 8, ('II', 'CBC'): 9, ('II', 'IR132'): 4,\n",
    "    ('III', 'IR8'): 6, ('III', 'CBC'): 7, ('III', 'IR132'): 6\n",
    "}\n",
    "model.objective = Objective(expr=sum(model.x[f, r] * yield_coefficients[(f, r)] for f in fields for r in rice_types), sense=maximize)\n",
    "\n",
    "# Constraints\n",
    "# Total hectares for each rice type\n",
    "model.total_hectares_ir8 = Constraint(expr=sum(model.x[f, 'IR8'] for f in fields) == 30)\n",
    "model.total_hectares_cbc = Constraint(expr=sum(model.x[f, 'CBC'] for f in fields) == 20)\n",
    "model.total_hectares_ir132 = Constraint(expr=sum(model.x[f, 'IR132'] for f in fields) == 40)\n",
    "\n",
    "# Total hectares for each field type\n",
    "model.total_field_i = Constraint(expr=sum(model.x['I', r] for r in rice_types) == 25)\n",
    "model.total_field_ii = Constraint(expr=sum(model.x['II', r] for r in rice_types) == 25)\n",
    "model.total_field_iii = Constraint(expr=sum(model.x['III', r] for r in rice_types) == 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9b9dd62e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer Community Edition 22.1.1.0\n",
      "  with Simplex, Mixed Integer & Barrier Optimizers\n",
      "5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21\n",
      "Copyright IBM Corp. 1988, 2022.  All Rights Reserved.\n",
      "\n",
      "Type 'help' for a list of available commands.\n",
      "Type 'help' followed by a command name for more\n",
      "information on commands.\n",
      "\n",
      "CPLEX> Logfile 'cplex.log' closed.\n",
      "Logfile 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmprivcz9xy.cplex.log' open.\n",
      "CPLEX> Problem 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmppxseebgu.pyomo.lp' read.\n",
      "Read time = 0.02 sec. (0.00 ticks)\n",
      "CPLEX> Problem name         : C:\\Users\\huyha\\AppData\\Local\\Temp\\tmppxseebgu.pyomo.lp\n",
      "Objective sense      : Maximize\n",
      "Variables            :       9  [General Integer: 9]\n",
      "Objective nonzeros   :       9\n",
      "Linear constraints   :       6  [Equal: 6]\n",
      "  Nonzeros           :      18\n",
      "  RHS nonzeros       :       6\n",
      "\n",
      "Variables            : Min LB: 0.000000         Max UB: all infinite   \n",
      "Objective nonzeros   : Min   : 4.000000         Max   : 9.000000       \n",
      "Linear constraints   :\n",
      "  Nonzeros           : Min   : 1.000000         Max   : 1.000000       \n",
      "  RHS nonzeros       : Min   : 20.00000         Max   : 40.00000       \n",
      "CPLEX> Version identifier: 22.1.1.0 | 2022-11-26 | 9160aff4d\n",
      "Found incumbent of value 560.000000 after 0.00 sec. (0.00 ticks)\n",
      "Tried aggregator 1 time.\n",
      "Reduced MIP has 6 rows, 9 columns, and 18 nonzeros.\n",
      "Reduced MIP has 0 binaries, 9 generals, 0 SOSs, and 0 indicators.\n",
      "Presolve time = 0.00 sec. (0.01 ticks)\n",
      "Tried aggregator 1 time.\n",
      "Detecting symmetries...\n",
      "Reduced MIP has 6 rows, 9 columns, and 18 nonzeros.\n",
      "Reduced MIP has 0 binaries, 9 generals, 0 SOSs, and 0 indicators.\n",
      "Presolve time = 0.00 sec. (0.01 ticks)\n",
      "MIP emphasis: balance optimality and feasibility.\n",
      "MIP search method: dynamic search.\n",
      "Parallel mode: deterministic, using up to 8 threads.\n",
      "Root relaxation solution time = 0.00 sec. (0.01 ticks)\n",
      "\n",
      "        Nodes                                         Cuts/\n",
      "   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap\n",
      "\n",
      "*     0+    0                          560.0000     1385.0000           147.32%\n",
      "*     0     0      integral     0      585.0000      585.0000        2    0.00%\n",
      "Elapsed time = 0.02 sec. (0.04 ticks, tree = 0.00 MB, solutions = 2)\n",
      "\n",
      "Root node processing (before b&c):\n",
      "  Real time             =    0.02 sec. (0.04 ticks)\n",
      "Parallel b&c, 8 threads:\n",
      "  Real time             =    0.00 sec. (0.00 ticks)\n",
      "  Sync time (average)   =    0.00 sec.\n",
      "  Wait time (average)   =    0.00 sec.\n",
      "                          ------------\n",
      "Total (root+branch&cut) =    0.02 sec. (0.04 ticks)\n",
      "\n",
      "Solution pool: 2 solutions saved.\n",
      "\n",
      "MIP - Integer optimal solution:  Objective =  5.8500000000e+02\n",
      "Solution time =    0.03 sec.  Iterations = 2  Nodes = 0\n",
      "Deterministic time = 0.04 ticks  (1.25 ticks/sec)\n",
      "\n",
      "CPLEX> Incumbent solution written to file 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmp6v2us8jd.cplex.sol'.\n",
      "CPLEX> Solution status: ok\n",
      "Field I, Rice IR8: -0.0 hectares\n",
      "Field I, Rice CBC: -0.0 hectares\n",
      "Field I, Rice IR132: 25.0 hectares\n",
      "Field II, Rice IR8: 5.0 hectares\n",
      "Field II, Rice CBC: 20.0 hectares\n",
      "Field II, Rice IR132: -0.0 hectares\n",
      "Field III, Rice IR8: 25.0 hectares\n",
      "Field III, Rice CBC: -0.0 hectares\n",
      "Field III, Rice IR132: 15.0 hectares\n"
     ]
    }
   ],
   "source": [
    "# Solve the model using CPLEX\n",
    "solver = SolverFactory('cplex', executable ='C:/Program Files/IBM/ILOG/CPLEX_Studio_Community2211/cplex/bin/x64_win64/cplex.exe')\n",
    "result = solver.solve(model, tee=True)\n",
    "\n",
    "# Print the results\n",
    "if result.solver.status == 'ok' and result.solver.termination_condition == 'optimal':\n",
    "    print(\"Solution status:\", result.solver.status)\n",
    "    for f in fields:\n",
    "        for r in rice_types:\n",
    "            print(f\"Field {f}, Rice {r}: {model.x[f, r].value} hectares\")\n",
    "elif result.solver.termination_condition == 'infeasible':\n",
    "    print(\"No feasible solution found.\")\n",
    "else:\n",
    "    print(\"Solver Status: \", result.solver.status)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caa53636-f1ff-49cf-8a6e-c3383beb6778",
   "metadata": {},
   "source": [
    "**b) Additional scenario where IR8 cannot be grown on Field Type I**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cf0df3d7-a461-40cf-9578-41c92124824d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer Community Edition 22.1.1.0\n",
      "  with Simplex, Mixed Integer & Barrier Optimizers\n",
      "5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21\n",
      "Copyright IBM Corp. 1988, 2022.  All Rights Reserved.\n",
      "\n",
      "Type 'help' for a list of available commands.\n",
      "Type 'help' followed by a command name for more\n",
      "information on commands.\n",
      "\n",
      "CPLEX> Logfile 'cplex.log' closed.\n",
      "Logfile 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmpt62jq_du.cplex.log' open.\n",
      "CPLEX> Problem 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmplvvdyrdt.pyomo.lp' read.\n",
      "Read time = 0.00 sec. (0.00 ticks)\n",
      "CPLEX> Problem name         : C:\\Users\\huyha\\AppData\\Local\\Temp\\tmplvvdyrdt.pyomo.lp\n",
      "Objective sense      : Maximize\n",
      "Variables            :       9  [General Integer: 9]\n",
      "Objective nonzeros   :       9\n",
      "Linear constraints   :       7  [Equal: 7]\n",
      "  Nonzeros           :      19\n",
      "  RHS nonzeros       :       6\n",
      "\n",
      "Variables            : Min LB: 0.000000         Max UB: all infinite   \n",
      "Objective nonzeros   : Min   : 4.000000         Max   : 9.000000       \n",
      "Linear constraints   :\n",
      "  Nonzeros           : Min   : 1.000000         Max   : 1.000000       \n",
      "  RHS nonzeros       : Min   : 20.00000         Max   : 40.00000       \n",
      "CPLEX> Version identifier: 22.1.1.0 | 2022-11-26 | 9160aff4d\n",
      "Found incumbent of value 585.000000 after 0.00 sec. (0.00 ticks)\n",
      "Tried aggregator 2 times.\n",
      "MIP Presolve eliminated 1 rows and 1 columns.\n",
      "Aggregator did 2 substitutions.\n",
      "Reduced MIP has 4 rows, 6 columns, and 12 nonzeros.\n",
      "Reduced MIP has 0 binaries, 6 generals, 0 SOSs, and 0 indicators.\n",
      "Presolve time = 0.02 sec. (0.01 ticks)\n",
      "Tried aggregator 1 time.\n",
      "Reduced MIP has 4 rows, 6 columns, and 12 nonzeros.\n",
      "Reduced MIP has 0 binaries, 6 generals, 0 SOSs, and 0 indicators.\n",
      "Presolve time = 0.00 sec. (0.01 ticks)\n",
      "MIP emphasis: balance optimality and feasibility.\n",
      "MIP search method: dynamic search.\n",
      "Parallel mode: deterministic, using up to 8 threads.\n",
      "Root relaxation solution time = 0.00 sec. (0.00 ticks)\n",
      "\n",
      "        Nodes                                         Cuts/\n",
      "   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap\n",
      "\n",
      "*     0+    0                          585.0000     1005.0000            71.79%\n",
      "      0     0        cutoff            585.0000                      2     --- \n",
      "\n",
      "Root node processing (before b&c):\n",
      "  Real time             =    0.02 sec. (0.03 ticks)\n",
      "Parallel b&c, 8 threads:\n",
      "  Real time             =    0.00 sec. (0.00 ticks)\n",
      "  Sync time (average)   =    0.00 sec.\n",
      "  Wait time (average)   =    0.00 sec.\n",
      "                          ------------\n",
      "Total (root+branch&cut) =    0.02 sec. (0.03 ticks)\n",
      "\n",
      "Solution pool: 1 solution saved.\n",
      "\n",
      "MIP - Integer optimal solution:  Objective =  5.8500000000e+02\n",
      "Solution time =    0.02 sec.  Iterations = 2  Nodes = 0\n",
      "Deterministic time = 0.03 ticks  (2.21 ticks/sec)\n",
      "\n",
      "CPLEX> Incumbent solution written to file 'C:\\Users\\huyha\\AppData\\Local\\Temp\\tmpo4i7p94e.cplex.sol'.\n",
      "CPLEX> Solution status: ok\n",
      "Field I, Rice IR8: 0.0 hectares\n",
      "Field I, Rice CBC: 20.0 hectares\n",
      "Field I, Rice IR132: 5.0 hectares\n",
      "Field II, Rice IR8: 25.0 hectares\n",
      "Field II, Rice CBC: 0.0 hectares\n",
      "Field II, Rice IR132: 0.0 hectares\n",
      "Field III, Rice IR8: 5.0 hectares\n",
      "Field III, Rice CBC: 0.0 hectares\n",
      "Field III, Rice IR132: 35.0 hectares\n"
     ]
    }
   ],
   "source": [
    "# Additional scenario where IR8 cannot be grown on Field Type I\n",
    "model.no_ir8_field_i = Constraint(expr=model.x['I', 'IR8'] == 0)\n",
    "\n",
    "# Solve the model using CPLEX\n",
    "solver = SolverFactory('cplex', executable ='C:/Program Files/IBM/ILOG/CPLEX_Studio_Community2211/cplex/bin/x64_win64/cplex.exe')\n",
    "result = solver.solve(model, tee=True)\n",
    "\n",
    "# Print the results\n",
    "if result.solver.status == 'ok' and result.solver.termination_condition == 'optimal':\n",
    "    print(\"Solution status:\", result.solver.status)\n",
    "    for f in fields:\n",
    "        for r in rice_types:\n",
    "            print(f\"Field {f}, Rice {r}: {model.x[f, r].value} hectares\")  # Changed from .value() to .value\n",
    "elif result.solver.termination_condition == 'infeasible':\n",
    "    print(\"No feasible solution found.\")\n",
    "else:\n",
    "    print(\"Solver Status: \", result.solver.status)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "443dd5df",
   "metadata": {},
   "source": [
    "**3. Original Problem #8 from tstran155**"
   ]
  },
  {
   "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.0\n",
      "Optimal solution: \n",
      "x[0][0] = 0.0\n",
      "x[0][1] = 25.0\n",
      "x[0][2] = 5.0\n",
      "x[1][0] = 20.0\n",
      "x[1][1] = 0.0\n",
      "x[1][2] = 0.0\n",
      "x[2][0] = 5.0\n",
      "x[2][1] = 0.0\n",
      "x[2][2] = 35.0\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.0\n",
      "Optimal solution: \n",
      "x[0][0] = 0.0\n",
      "x[0][1] = 25.0\n",
      "x[0][2] = 5.0\n",
      "x[1][0] = 20.0\n",
      "x[1][1] = 0.0\n",
      "x[1][2] = 0.0\n",
      "x[2][0] = 5.0\n",
      "x[2][1] = 0.0\n",
      "x[2][2] = 35.0\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": "d38f1a17-ffbc-421c-8ecd-cf9b2700ac08",
   "metadata": {},
   "source": [
    "**3. Discussions**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac314398",
   "metadata": {},
   "source": [
    "a. OR-tools solver:\n",
    "\n",
    "Results of the first and second scenario where IR8 cannot be grown on Field Type I are indeed identical. The unchanged solution when adding the constraint that IR8 cannot be grown in Field I suggests two things: firstly, that the restriction aligns with the optimal planting strategy derived without it, and secondly, it could indicate a potential area of focus for further investigation like planting strategy, costs, or yields could change this result.\n",
    "\n",
    "The solutions provided by the LP and MILP techniques are different. Here is the insights:\n",
    "\n",
    "Allocation of Rice Types:\n",
    "\n",
    "- IR8: The LP solution does not utilize IR8 in Field II and uses all 25 hectares in Field III. In contrast, the MILP solution fully utilizes Field II for IR8 planting, which suggests different optimization dynamics between integer constraints (MILP) and continuous variables (LP).\n",
    "    \n",
    "- IR132: The LP solution has a significant allocation of IR132 in Field II (20 hectares), which is not present in the MILP solution. This again reflects the flexibility in continuous allocation in LP that might not be as effective or feasible when restricted to integer values in MILP.\n",
    "\n",
    "Field Utilization:\n",
    "\n",
    "- Field II in LP has only IR132 planted and utilizes only 20 out of 25 possible hectares, indicating a potentially suboptimal use of available land under the constraints and yield assumptions of LP.\n",
    "    \n",
    "- Field III shows differing allocations between IR8 and IR132 across the solutions, which might indicate that the profitability or yield per hectare changes dramatically based on how the constraints are treated (integer vs. continuous).\n",
    "\n",
    "b. CPLEX solver:\n",
    "\n",
    "- The optimal value reached is 585.0, which represents the maximum total yield based on your model's constraints and objective function.\n",
    "  \n",
    "- Rice IR132 is heavily allocated to Field I, which might indicate it's the most profitable or effective use of that field, while rice CBC is primarily allocated to Field II, with Rice IR8 also placed significantly in Field III and minimally in Field II.\n",
    "  \n",
    "- The results are reasonable and align with typical outcomes in agricultural optimization, where certain crops are prioritized in specific fields to optimize the use of available resources and meet demand constraints effectively.\n",
    "  \n",
    "c. LP vs. MILP:\n",
    "\n",
    "- LP offers more flexibility in how variables can be allocated, potentially allowing for more \"fine-tuned\" solutions at a fractional level, which might not always be practical or feasible in real-world scenarios.\n",
    "  \n",
    "- MILP provides solutions that are more practical for real-world applications where the decision variables must be whole numbers (like the number of hectares that can realistically be planted).\n",
    "  \n",
    "- Depending on the exact nature of the constraints and objective function, MILP might yield a less optimal solution in terms of the raw numbers but a more feasible and practically implementable solution.\n"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}