{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Consider the polynomial optimization problem of\n",
    "minimizing the polynomial $x^3 - x^2 + 2xy -y^2 + y^3$\n",
    "over the polyhedron defined by the inequalities $x \\ge 0, y \\ge 0$ and $x + y \\geq 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1//4, 0)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x y\n",
    "p = x^3 - x^2 + 2x*y -y^2 + y^3\n",
    "using SemialgebraicSets\n",
    "S = @set x >= 0 && y >= 0 && x + y >= 1\n",
    "p(x=>1, y=>0), p(x=>1//2, y=>1//2), p(x=>0, y=>1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but [Ipopt](https://github.com/JuliaOpt/Ipopt.jl/) only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "******************************************************************************\n",
      "\n",
      "termination_status(model) = LOCALLY_SOLVED::TerminationStatusCode = 4\n",
      "value(a) = 0.4999999966705758\n",
      "value(b) = 0.49999999667062867\n",
      "objective_value(model) = 0.24999999500590342\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.24999999500590342"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuMP\n",
    "using Ipopt\n",
    "model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))\n",
    "@variable(model, a >= 0)\n",
    "@variable(model, b >= 0)\n",
    "@constraint(model, a + b >= 1)\n",
    "@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)\n",
    "optimize!(model)\n",
    "\n",
    "@show termination_status(model)\n",
    "@show value(a)\n",
    "@show value(b)\n",
    "@show objective_value(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the problem can be written equivalently as follows using [registered functions](http://www.juliaopt.org/JuMP.jl/v0.19.0/nlp/#User-defined-Functions-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = LOCALLY_SOLVED::TerminationStatusCode = 4\n",
      "value(a) = 0.499999995002878\n",
      "value(b) = 0.4999999950028773\n",
      "objective_value(model) = 0.24999999250431656\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.24999999250431656"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuMP\n",
    "using Ipopt\n",
    "model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))\n",
    "@variable(model, a >= 0)\n",
    "@variable(model, b >= 0)\n",
    "@constraint(model, a + b >= 1)\n",
    "peval(a, b) = p(x=>a, y=>b)\n",
    "register(model, :peval, 2, peval, autodiff=true)\n",
    "@NLobjective(model, Min, peval(a, b))\n",
    "optimize!(model)\n",
    "\n",
    "@show termination_status(model)\n",
    "@show value(a)\n",
    "@show value(b)\n",
    "@show objective_value(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sum-of-Squares approach\n",
    "\n",
    "We will now see how to find the optimal solution using Sum of Squares Programming."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using SumOfSquares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first need to pick an SDP solver, see [here](http://www.juliaopt.org/JuMP.jl/dev/installation/#Getting-Solvers-1) for a list of the available choices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using CSDP\n",
    "optimizer = with_optimizer(CSDP.Optimizer, printlevel=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using MosekTools\n",
    "optimizer = with_optimizer(Mosek.Optimizer, QUIET=true);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Sum-of-Squares certificate that $p \\ge \\alpha$ over the domain `S`, ensures that $\\alpha$ is a lower bound to the polynomial optimization problem.\n",
    "The following program searches for the largest upper bound and find zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = OPTIMAL::TerminationStatusCode = 1\n",
      "objective_value(model) = -2.0092710828478744e-10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-2.0092710828478744e-10"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = SOSModel(optimizer)\n",
    "@variable(model, α)\n",
    "@objective(model, Max, α)\n",
    "@constraint(model, c3, p >= α, domain = S)\n",
    "\n",
    "optimize!(model)\n",
    "@show termination_status(model)\n",
    "@show objective_value(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the solution $(1/2, 1/2)$ found by Ipopt of objective value $1/4$\n",
    "and this certificate of lower bound $0$ we know that the optimal objective value is in the interval $[0, 1/4]$\n",
    "but we still do not know what it is (if we consider that we did not try the solutions $(1, 0)$ and $(0, 1)$ as done in the introduction).\n",
    "If the dual of the constraint `c3` was atomic, its atoms would have given optimal solutions of objective value $0$ but that is not the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using MultivariateMoments\n",
    "ν3 = moment_matrix(c3)\n",
    "extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortunately, there is a hierarchy of programs with increasingly better programs that can be solved until we get one with atom dual variables.\n",
    "This comes from the way the Sum-of-Squares constraint with domain `S` is formulated.\n",
    "The polynomial $p - \\alpha$ is guaranteed to be nonnegative over the domain `S` if there exists Sum-of-Squares polynomials $s_0$, $s_1$, $s_2$, $s_3$ such that\n",
    "$$ p - \\alpha = s_0 + s_1 x + s_2 y + s_3 (x + y - 1). $$\n",
    "Indeed, in the domain `S`, $x$, $y$ and $x + y - 1$ are nonnegative so the right-hand side is a sum of squares hence is nonnegative.\n",
    "Once the degrees of $s_1$, $s_2$ and $s_3$ have been decided, the degree needed for $s_0$ will be determined but we have a freesom in choosing the degrees of $s_1$, $s_2$ and $s_3$.\n",
    "By default, they are chosen so that the degrees of $s_1 x$, $s_2 y$ and $s_3 (x + y - 1)$ match those of $p - \\alpha$ but this can be overwritten using the $maxdegree$ keyword argument.\n",
    "\n",
    "### The maxdegree keyword argument\n",
    "\n",
    "The maximum total degree (i.e. maximum sum of the exponents of $x$ and $y$) of the monomials of $p$ is 3 so the constraint in the program above is equivalent to `@constraint m p >= α domain = S maxdegree = 3`..\n",
    "That is, since $x$, $y$ and $x + y - 1$ have total degree 1, the sum of squares polynomials $s_1$, $s_2$ and $s_3$ have been chosen with maximum total degree $2$.\n",
    "Since these polynomials are sums of squares, their degree must be even so the next maximum total degree to try is 4.\n",
    "For this reason, the keywords `maxdegree = 4` and `maxdegree = 5` have the same effect in this example.\n",
    "In general, if the polynomials in the domain are not all odd or all even, each value of `maxdegree` has different effect in the choice of the maximum total degree of $s_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = OPTIMAL::TerminationStatusCode = 1\n",
      "objective_value(model) = 8.700006062680539e-8\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8.700006062680539e-8"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = SOSModel(optimizer)\n",
    "@variable(model, α)\n",
    "@objective(model, Max, α)\n",
    "@constraint(model, c5, p >= α, domain = S, maxdegree = 5)\n",
    "\n",
    "optimize!(model)\n",
    "@show termination_status(model)\n",
    "@show objective_value(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, the dual variable is atomic as it is the moments of the measure\n",
    "$$0.5 \\delta(x-1, y) + 0.5 \\delta(x, y-1)$$\n",
    "where $\\delta(x, y)$ is the dirac measure centered at $(0, 0)$.\n",
    "Therefore the program provides both a certificate that $0$ is a lower bound and a certificate that it is also an upper bound since it is attained at the global minimizers $(1, 0)$ and $(0, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atomic measure on the variables x, y with 2 atoms:\n",
       " at [-0.000242397, 1.00024] with weight 0.4936307922624579\n",
       " at [0.999972, 2.76679e-5] with weight 0.5063248377957446"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using MultivariateMoments\n",
    "μ5 = moment_matrix(c5)\n",
    "extractatoms(μ5, 1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A deeper look into atom extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `extractatoms` function requires a `ranktol` argument that we have set to `1e-3` in the preceding section.\n",
    "This argument is used to transform the dual variable into a system of polynomials equations whose solutions give the atoms.\n",
    "This transformation uses the SVD decomposition of the matrix of moments and discard the equations corresponding to a singular value lower than `ranktol`.\n",
    "When this system of equation has an infinite number of solutions, `extractatoms` concludes that the measure is not atomic.\n",
    "For instance, with `maxdegree = 3`, we obtain the system\n",
    "$$x + y = 1$$\n",
    "which contains a whole line of solution.\n",
    "This explains `extractatoms` returned `nothing`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Algebraic Set defined by 1 equality\n",
       " -x - 0.9999999999999998*y + 1.0000000000569824 == 0\n"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ν3 = moment_matrix(c3)\n",
    "MultivariateMoments.computesupport!(ν3, 1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `maxdegree = 5`, we obtain the system\n",
    "\\begin{align}\n",
    "  x + y & = 1\\\\\n",
    "  y^2 & = y\\\\\n",
    "  xy & = 0\\\\\n",
    "  x^2 + y & = 1\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Algebraic Set defined by 4 equalities\n",
       " -x - 0.9999999962812296*y + 1.0000000580930948 == 0\n",
       " -y^2 + 1.0002701267987304*y - 2.7674560303336815e-5 == 0\n",
       " -x*y - 1.293349230893601e-6*y + 3.3938404378219955e-5 == 0\n",
       " -x^2 - 1.0002601608230963*y + 1.0002330609868753 == 0\n"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ν5 = moment_matrix(c5)\n",
    "MultivariateMoments.computesupport!(ν5, 1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system can be reduced to the equivalent system\n",
    "\\begin{align}\n",
    "  x + y & = 1\\\\\n",
    "  y^2 & = y\n",
    "\\end{align}\n",
    "which has the solutions $(0, 1)$ and $(1, 0)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Algebraic Set defined by 2 equalities\n",
       " x + 0.9999999962812296*y - 1.0000000580930948 == 0\n",
       " y^2 - 1.0002701267987304*y + 2.7674560303336815e-5 == 0\n"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))\n",
    "ν5.support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `extractatoms` then reuse the matrix of moments to find the weights $1/2$, $1/2$ corresponding to the diracs centered respectively at $(0, 1)$ and $(1, 0)$.\n",
    "This details the how the function obtained the result\n",
    "$$0.5 \\delta(x-1, y) + 0.5 \\delta(x, y-1)$$\n",
    "given in the previous section."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}