src/python/docs/solve_with_cplex_cbc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff0963df",
   "metadata": {},
   "source": [
    "## Solution of the mixed integer linear pump scheduling formulation using IBM's CPLEX solver (commercial) and COINOR's CBC solver (free open source)\n",
    "\n",
    "### T.Janus\n",
    "### 21/05/2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03863d8f",
   "metadata": {},
   "source": [
    "### README\n",
    "\n",
    "The linear programme represented in a `.mps` format is read and solved with two solvers - a commercial **IBM ILOG CPLEX** solver requiring an academic/business license and an open-source free **COINOR CBC** solver. **CPLEX** solver is executd using the IBM's CPLEX Python binding. **CBC** solver is executed with COINOR's [PYTHON-MIP](https://www.python-mip.com/) library for the modeling and solution of Mixed-Integer Linear programs (MIPs) that supports creating large scale linear programmes and comes packaged with CBC binaries.\n",
    "\n",
    "#### TODO:\n",
    "1. Install and test solving the problem with Gurobi using PYTHON-MIP\n",
    "2. Introduction os special ordered sets (SOS) https://docs.python-mip.com/en/latest/sos.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "389b346c",
   "metadata": {},
   "source": [
    "### 1. OPTIMIZE WITH CPLEX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bae418e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import os\n",
    "import cplex\n",
    "import mip\n",
    "# import docplex.mp # Library for creating optimization models in Python\n",
    "import scipy.io\n",
    "import cvxpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9daf9975",
   "metadata": {},
   "source": [
    "TODO : Run via CVXPY\n",
    "https://www.cvxpy.org/install/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24401cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "MODEL_MPS_FILE = pathlib.Path(\"../data/2p1t/2p1t.mps\")\n",
    "MODEL_MATRICES_FILE = pathlib.Path(\"../data/2p1t/2p1t_model.mat\")\n",
    "MODEL_LP_FILE = pathlib.Path(\"../data/2p1t/2p1t.lp\")\n",
    "CPLEX_RESULTS_FILE = pathlib.Path(\"../outputs/2p1t/x_optim_cplex.mat\")\n",
    "CBC_RESULTS_FILE = pathlib.Path(\"../outputs/2p1t/x_optim_cbc.mat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7daab0a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parrarel (multihreaded) execution\n",
    "num_threads = os.cpu_count()\n",
    "print(num_threads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2020a65",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load matrices from Matlab (for testing purposes)\n",
    "matlab_model = scipy.io.loadmat(MODEL_MATRICES_FILE)\n",
    "# print(matlab_model['ub'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc82ffdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert MPS file to LP file\n",
    "instance = mip.Model()\n",
    "instance.read(MODEL_MPS_FILE.as_posix())\n",
    "instance.write(MODEL_LP_FILE.as_posix())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f97b05a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a CPLEX problem object\n",
    "problem = cplex.Cplex()\n",
    "problem.parameters.threads.set(num_threads)\n",
    "problem.parameters.benders.strategy = -1\n",
    "# problem.parameters.lpmethod = 4\n",
    "# Read the MPS file\n",
    "problem.read(MODEL_MPS_FILE.as_posix())\n",
    "# Set the problem type to MILP\n",
    "problem.set_problem_type(cplex.Cplex.problem_type.MILP)\n",
    "# Set the tolerance gap\n",
    "problem.parameters.mip.tolerances.mipgap.set(0.05)\n",
    "# Solve the problem\n",
    "problem.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57a19e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.report()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99282a93",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.print_information()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10ad3e13",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Solution status: {problem.solution.get_status_string()}\")\n",
    "print(f\"Objective value: {problem.solution.get_objective_value()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "357562a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_optim = problem.solution.get_values()\n",
    "scipy.io.savemat(CPLEX_RESULTS_FILE, {\"x\": x_optim, \"obj\": problem.solution.get_objective_value()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a29b689",
   "metadata": {},
   "outputs": [],
   "source": [
    "CPLEX_RESULTS_FILE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b43bae1d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"Solution vector:\")\n",
    "for i, name in enumerate(problem.variables.get_names()):\n",
    "    print(f\"{name} = {problem.solution.get_values(i):.3f}\")\n",
    "    if i > 5:\n",
    "        print(\"...\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dbdc4a6",
   "metadata": {},
   "source": [
    "### 2. OPTIMIZE WITH CBC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efa63052",
   "metadata": {},
   "outputs": [],
   "source": [
    "instance = mip.Model(solver_name=mip.CBC)\n",
    "instance.read(MODEL_MPS_FILE.as_posix())\n",
    "instance.write(MODEL_LP_FILE.as_posix())\n",
    "# Set optimization parameters\n",
    "instance.threads = -1\n",
    "instance.max_gap = 0.15\n",
    "#instance.cut_passes = 10\n",
    "#instance.clique = -1  # -1 means automatic, 0 disables it, 1 enables it and 2 enables more aggressive clique generation\n",
    "#instance.cuts = 2\n",
    "#instance.emphasis = 2\n",
    "#instance.infeas_tol = 1e-4\n",
    "#instance.preprocess = 1\n",
    "#instance.pump_passes = 10\n",
    "instance.lp_method = 1\n",
    "optim_status = instance.optimize(max_seconds=600)\n",
    "if optim_status == mip.OptimizationStatus.OPTIMAL:\n",
    "    print('\\n Optimal solution cost {} found'.format(instance.objective_value))\n",
    "elif optim_status == mip.OptimizationStatus.FEASIBLE:\n",
    "    print('\\n Feasible solution cost {} found, best possible: {}'.format(instance.objective_value, instance.objective_bound))\n",
    "elif optim_status == mip.OptimizationStatus.NO_SOLUTION_FOUND:\n",
    "    print('\\n No feasible solution found, lower bound is: {}'.format(instance.objective_bound))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad3ae441",
   "metadata": {},
   "outputs": [],
   "source": [
    "instance.objective.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db72fdbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "instance.objective_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c899b08f",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_optim_cbc = [var.x for var in instance.vars]\n",
    "scipy.io.savemat(CBC_RESULTS_FILE, {\"x\": x_optim_cbc})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7616688b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#diff = [x1 - x2 for x1, x2 in zip(x_optim, x_optim_cbc)]\n",
    "#diff"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}