{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Extended formulation"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In this notebook, we show how to work with the extended formulation of a polyhedron.\n",
    "The convex hull of the union of polyhedra that are H-represented can be obtained\n",
    "as the projection of a H-representation [Theorem 3.3, B85].\n",
    "In order to use the resulting polyhedron as a constraint set in an optimization problem,\n",
    "there is no need to compute the resulting H-representation of this projection.\n",
    "Moreover, other operations such as intersection are also implemented between extended H-representations.\n",
    "We illustrate this with a simple example. We start by defining the H-representation of a square with JuMP."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "[B85] Balas, E., 1985.\n",
    "*Disjunctive programming and a hierarchy of relaxations for discrete optimization problems*.\n",
    "SIAM Journal on Algebraic Discrete Methods, 6(3), pp.466-486."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "H-representation LPHRep{Float64}:\n4-element iterator of HalfSpace{Float64, SparseArrays.SparseVector{Float64, Int64}}:\n HalfSpace(  [1]  =  -1.0, 1.0)\n HalfSpace(  [2]  =  -1.0, 1.0)\n HalfSpace(  [1]  =  1.0, 1.0)\n HalfSpace(  [2]  =  1.0, 1.0)"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using JuMP\n",
    "model = Model()\n",
    "@variable(model, -1 <= x <= 1)\n",
    "@variable(model, -1 <= y <= 1)\n",
    "using Polyhedra\n",
    "square = hrep(model)"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the names of the JuMP variables are used as the names of the corresponding dimensions for\n",
    "the polyhedron `square`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Vector{String}:\n \"x\"\n \"y\""
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "dimension_names(square)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "In the following, `diagonal` and `antidiag` are projections of extended H-representations of dimension 7\n",
    "hence `diamond` is a projection of an extended H-representation of dimensions 12."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Polyhedra.Projection{LPHRep{Float64}, Base.OneTo{Int64}}(HyperPlane(  [1]  =  1.0\n  [3]  =  -1.0\n  [5]  =  -1.0, 0.0) ∩ HyperPlane(  [2]  =  1.0\n  [4]  =  -1.0\n  [6]  =  -1.0, 0.0) ∩ HyperPlane(  [1 ]  =  1.0\n  [8 ]  =  -1.0\n  [10]  =  -1.0, 0.0) ∩ HyperPlane(  [2 ]  =  1.0\n  [9 ]  =  -1.0\n  [11]  =  -1.0, 0.0) ∩ HalfSpace(  [3]  =  -1.0\n  [7]  =  -3.0, 0.0) ∩ HalfSpace(  [4]  =  -1.0\n  [7]  =  -3.0, 0.0) ∩ HalfSpace(  [3]  =  1.0\n  [7]  =  1.0, 0.0) ∩ HalfSpace(  [4]  =  1.0\n  [7]  =  1.0, 0.0) ∩ HalfSpace(  [5]  =  -1.0\n  [7]  =  -1.0, -1.0) ∩ HalfSpace(  [6]  =  -1.0\n  [7]  =  -1.0, -1.0) ∩ HalfSpace(  [5]  =  1.0\n  [7]  =  3.0, 3.0) ∩ HalfSpace(  [6]  =  1.0\n  [7]  =  3.0, 3.0) ∩ HalfSpace(  [7]  =  -1.0, 0.0) ∩ HalfSpace(  [7]  =  1.0, 1.0) ∩ HalfSpace(  [8 ]  =  -1.0\n  [12]  =  -3.0, 0.0) ∩ HalfSpace(  [9 ]  =  -1.0\n  [12]  =  1.0, 0.0) ∩ HalfSpace(  [8 ]  =  1.0\n  [12]  =  1.0, 0.0) ∩ HalfSpace(  [9 ]  =  1.0\n  [12]  =  -3.0, 0.0) ∩ HalfSpace(  [10]  =  -1.0\n  [12]  =  -1.0, -1.0) ∩ HalfSpace(  [11]  =  -1.0\n  [12]  =  3.0, 3.0) ∩ HalfSpace(  [10]  =  1.0\n  [12]  =  3.0, 3.0) ∩ HalfSpace(  [11]  =  1.0\n  [12]  =  -1.0, -1.0) ∩ HalfSpace(  [12]  =  -1.0, 0.0) ∩ HalfSpace(  [12]  =  1.0, 1.0), Base.OneTo(2))"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "diagonal = convexhull(translate(square, [-2, -2]), translate(square, [2, 2]))\n",
    "antidiag = convexhull(translate(square, [-2,  2]), translate(square, [2, -2]))\n",
    "diamond = diagonal ∩ antidiag"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the names the first two dimensions are still identical to the names of the JuMP variables\n",
    "and the auxiliary variables have no name."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "12-element Vector{String}:\n \"x\"\n \"y\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\"\n \"\""
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "dimension_names(diamond.set)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We don't need to compute the result of the projection to solve an optimization problem\n",
    "over `diamond`. For instance, to compute the maximal value that `y` can take over this\n",
    "polytope with the GLPK solver, we can do as follows.\n",
    "Note that if we use anonymous JuMP variables, the name of the JuMP variables will be\n",
    "the names of the corresponding dimensions of the polyhedron.\n",
    "Therefore, we can retrieve the JuMP variable according to the corresponding dimension name\n",
    "with `variable_by_name`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(0.0, 2.0)"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "import GLPK\n",
    "model = Model(GLPK.Optimizer)\n",
    "@variable(model, [1:2] in diamond)\n",
    "x = variable_by_name(model, \"x\")\n",
    "y = variable_by_name(model, \"y\")\n",
    "@objective(model, Max, y)\n",
    "optimize!(model)\n",
    "value(x), value(y)"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "In the optimization problem, above, the auxiliary variables of the extended formulation\n",
    "are transparently added inside a bridge.\n",
    "To manipulate the auxiliary variables, one can use the extended H-representation directly instead of its projection.\n",
    "Note that as the auxiliary dimensions have no name, we cannot use `variable_by_name` to retrieve the corresponding\n",
    "JuMP variables. We can instead catch the returned value of `@variable` in some variable `v` in order to use anonymous JuMP\n",
    "variables while still assigning the created JuMP variables to `v`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "12-element Vector{Float64}:\n  0.0\n  2.0\n -0.75\n -0.25\n  0.75\n  2.25\n  0.25\n -0.75\n  2.25\n  0.75\n -0.25\n  0.75"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "import GLPK\n",
    "model = Model(GLPK.Optimizer)\n",
    "v = @variable(model, [1:12] in diamond.set)\n",
    "y = variable_by_name(model, \"y\")\n",
    "@objective(model, Max, y)\n",
    "optimize!(model)\n",
    "value.(v)"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.1"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.1",
   "language": "julia"
  }
 },
 "nbformat": 4
}