{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Bound on Global Extremum"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: SOSTOOLS' SOSDEMO3 (See Section 4.3 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf))"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(x1, x2)"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x1 x2"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Goldstein-Price function $f(x)$ is defined as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "600 + 720x2 + 720x1 + 3060x2² - 4680x1x2 + 1260x1² + 12288x2³ - 19296x1x2² + 7344x1²x2 - 1072x1³ + 14346x2⁴ - 23616x1x2³ + 7776x1²x2² + 5784x1³x2 - 2454x1⁴ + 1944x2⁵ - 11880x1x2⁴ + 5040x1²x2³ + 9840x1³x2² - 7680x1⁴x2 + 1344x1⁵ - 4428x2⁶ - 1188x1x2⁵ + 8730x1²x2⁴ + 1240x1³x2³ - 5370x1⁴x2² - 168x1⁵x2 + 952x1⁶ - 648x2⁷ + 1944x1x2⁶ + 3672x1²x2⁵ - 3480x1³x2⁴ - 4080x1⁴x2³ + 2592x1⁵x2² + 1344x1⁶x2 - 768x1⁷ + 729x2⁸ + 972x1x2⁷ - 1458x1²x2⁶ - 1836x1³x2⁵ + 1305x1⁴x2⁴ + 1224x1⁵x2³ - 648x1⁶x2² - 288x1⁷x2 + 144x1⁸",
      "text/latex": "$$ 600 + 720x2 + 720x1 + 3060x2^{2} - 4680x1x2 + 1260x1^{2} + 12288x2^{3} - 19296x1x2^{2} + 7344x1^{2}x2 - 1072x1^{3} + 14346x2^{4} - 23616x1x2^{3} + 7776x1^{2}x2^{2} + 5784x1^{3}x2 - 2454x1^{4} + 1944x2^{5} - 11880x1x2^{4} + 5040x1^{2}x2^{3} + 9840x1^{3}x2^{2} - 7680x1^{4}x2 + 1344x1^{5} - 4428x2^{6} - 1188x1x2^{5} + 8730x1^{2}x2^{4} + 1240x1^{3}x2^{3} - 5370x1^{4}x2^{2} - 168x1^{5}x2 + 952x1^{6} - 648x2^{7} + 1944x1x2^{6} + 3672x1^{2}x2^{5} - 3480x1^{3}x2^{4} - 4080x1^{4}x2^{3} + 2592x1^{5}x2^{2} + 1344x1^{6}x2 - 768x1^{7} + 729x2^{8} + 972x1x2^{7} - 1458x1^{2}x2^{6} - 1836x1^{3}x2^{5} + 1305x1^{4}x2^{4} + 1224x1^{5}x2^{3} - 648x1^{6}x2^{2} - 288x1^{7}x2 + 144x1^{8} $$"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "f1 = x1 + x2 + 1\n",
    "f2 = 19 - 14x1 + 3x1^2 - 14x2 + 6x1*x2 + 3x2^2\n",
    "f3 = 2x1 - 3x2\n",
    "f4 = 18 - 32x1 + 12x1^2 + 48x2 - 36x1*x2 + 27x2^2\n",
    "f = (1 + f1^2 * f2) * (30 + f3^2 * f4)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices.\n",
    "We use `SOSModel` instead of `Model` to be able to use the `>=` syntax for Sum-of-Squares constraints."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SumOfSquares\n",
    "using CSDP\n",
    "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n",
    "model = SOSModel(solver);"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "We create the decision variable $\\gamma$ that will be the lower bound to the Goldstein-Price function.\n",
    "We maximize it to have the highest possible lower bound."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "γ",
      "text/latex": "$ γ $"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "@variable(model, γ)\n",
    "@objective(model, Max, γ)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We constrain $\\gamma$ to be a lower bound with the following constraint\n",
    "that ensures that $f(x_1, x_2) \\ge \\gamma$ for all $x_1, x_2$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@constraint(model, f >= γ)\n",
    "\n",
    "JuMP.optimize!(model)"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "We verify that the solver has found a feasible solution:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "NEARLY_FEASIBLE_POINT::ResultStatusCode = 2"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "JuMP.primal_status(model)"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can now obtain the lower bound either with `value(γ)` or `objective_value(model)`:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2.999986072324649"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "objective_value(model)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "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.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}