{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Bilinear terms"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: [Floudas1999; Section 3.1](@cite) and [Lasserre2009; Table 5.1](@cite)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Introduction"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Consider the polynomial optimization problem from [Floudas1999; Section 3.1](@cite)."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Basic semialgebraic Set defined by no equality\n22 inequalities\n 1.0 - 0.0025*x[6] - 0.0025*x[4] ≥ 0\n 1.0 - 0.0025*x[7] - 0.0025*x[5] + 0.0025*x[4] ≥ 0\n 1.0 - 0.01*x[8] + 0.01*x[5] ≥ 0\n 83333.33333333333 - 8333.33252*x[4] - 100.0*x[1] + x[1]*x[6] ≥ 0\n -1250.0*x[5] + 1250.0*x[4] + x[2]*x[7] - x[2]*x[4] ≥ 0\n -1.25e6 + 2500.0*x[5] + x[3]*x[8] - x[3]*x[5] ≥ 0\n -100.0 + x[1] ≥ 0\n 10000.0 - x[1] ≥ 0\n -1000.0 + x[2] ≥ 0\n 10000.0 - x[2] ≥ 0\n -1000.0 + x[3] ≥ 0\n 10000.0 - x[3] ≥ 0\n -10.0 + x[4] ≥ 0\n 1000.0 - x[4] ≥ 0\n -10.0 + x[5] ≥ 0\n 1000.0 - x[5] ≥ 0\n -10.0 + x[6] ≥ 0\n 1000.0 - x[6] ≥ 0\n -10.0 + x[7] ≥ 0\n 1000.0 - x[7] ≥ 0\n -10.0 + x[8] ≥ 0\n 1000.0 - x[8] ≥ 0\n"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x[1:8]\n",
    "p = sum(x[1:3])\n",
    "using SumOfSquares\n",
    "K = @set 0.0025 * (x[4] + x[6]) <= 1 &&\n",
    "    0.0025 * (-x[4] + x[5] + x[7]) <= 1 &&\n",
    "    0.01 * (-x[5] + x[8]) <= 1 &&\n",
    "    100x[1] - x[1] * x[6] + 8333.33252x[4] <= 250000/3 &&\n",
    "    x[2] * x[4] - x[2] * x[7] - 1250x[4] + 1250x[5] <= 0 &&\n",
    "    x[3] * x[5] - x[3] * x[8] - 2500x[5] + 1250000 <= 0 &&\n",
    "    100 <= x[1] && x[1] <= 10000 &&\n",
    "    1000 <= x[2] && x[2] <= 10000 &&\n",
    "    1000 <= x[3] && x[3] <= 10000 &&\n",
    "    10 <= x[4] && x[4] <= 1000 &&\n",
    "    10 <= x[5] && x[5] <= 1000 &&\n",
    "    10 <= x[6] && x[6] <= 1000 &&\n",
    "    10 <= x[7] && x[7] <= 1000 &&\n",
    "    10 <= x[8] && x[8] <= 1000"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will now see how to find the optimal solution using Sum of Squares Programming.\n",
    "We first 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."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Clarabel.MOIwrapper.Optimizer"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "import Clarabel\n",
    "solver = Clarabel.Optimizer"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "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 function searches for the largest lower bound and finds zero using the `d`th level of the hierarchy`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "solve (generic function with 1 method)"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "function solve(d)\n",
    "    model = SOSModel(solver)\n",
    "    @variable(model, α)\n",
    "    @objective(model, Max, α)\n",
    "    @constraint(model, c, p >= α, domain = K, maxdegree = d)\n",
    "    optimize!(model)\n",
    "    println(solution_summary(model))\n",
    "    return model\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "The first level of the hierarchy gives a lower bound of `2100`"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------------\n",
      "           Clarabel.jl v0.7.1  -  Clever Acronym              \n",
      "                   (c) Paul Goulart                          \n",
      "                University of Oxford, 2022                   \n",
      "-------------------------------------------------------------\n",
      "\n",
      "problem:\n",
      "  variables     = 21\n",
      "  constraints   = 29\n",
      "  nnz(P)        = 0\n",
      "  nnz(A)        = 64\n",
      "  cones (total) = 2\n",
      "    : Zero        = 1,  numel = 9\n",
      "    : Nonnegative = 1,  numel = 20\n",
      "\n",
      "settings:\n",
      "  linear algebra: direct / qdldl, precision: Float64\n",
      "  max iter = 200, time limit = Inf,  max step = 0.990\n",
      "  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,\n",
      "  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32\n",
      "  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07\n",
      "  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, \n",
      "               max iter = 10, stop ratio = 5.0\n",
      "  equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05\n",
      "               max iter = 10\n",
      "\n",
      "iter    pcost        dcost       gap       pres      dres      k/t        μ       step      \n",
      "---------------------------------------------------------------------------------------------\n",
      "  0  -2.3512e+03  -2.3512e+03  1.93e-16  2.04e-02  5.29e-02  1.00e+00  9.80e+01   ------   \n",
      "  1  -2.1101e+03  -2.1100e+03  9.17e-06  5.10e-04  9.42e-04  3.70e-02  1.86e+00  9.81e-01  \n",
      "  2  -2.1001e+03  -2.1001e+03  9.66e-08  5.27e-06  9.66e-06  3.84e-04  1.93e-02  9.90e-01  \n",
      "  3  -2.1000e+03  -2.1000e+03  9.66e-10  5.27e-08  9.66e-08  3.84e-06  1.93e-04  9.90e-01  \n",
      "  4  -2.1000e+03  -2.1000e+03  9.66e-12  5.27e-10  9.66e-10  3.84e-08  1.93e-06  9.90e-01  \n",
      "---------------------------------------------------------------------------------------------\n",
      "Terminated with status = solved\n",
      "solve time =  115ms\n",
      "* Solver : Clarabel\n",
      "\n",
      "* Status\n",
      "  Result count       : 1\n",
      "  Termination status : OPTIMAL\n",
      "  Message from the solver:\n",
      "  \"SOLVED\"\n",
      "\n",
      "* Candidate solution (result #1)\n",
      "  Primal status      : FEASIBLE_POINT\n",
      "  Dual status        : FEASIBLE_POINT\n",
      "  Objective value    : 2.10000e+03\n",
      "  Dual objective value : 2.10000e+03\n",
      "\n",
      "* Work counters\n",
      "  Solve time (sec)   : 1.14878e-01\n",
      "  Barrier iterations : 4\n",
      "\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "model2 = solve(2)\n",
    "nothing # hide"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "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
}