{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: Conic Programming\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Originally Contributed by**: Arpit Bhatia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is a Cone?\n",
    "A subset $C$ of a vector space $V$ is a cone if $\\forall x \\in C$ and positive scalars $\\alpha$,\n",
    "the product $\\alpha x \\in C$. A cone C is a convex cone if $\\alpha x + \\beta y \\in C$,\n",
    "for any positive scalars $\\alpha, \\beta$, and any $x, y \\in C$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conic Programming\n",
    "Conic programming problems are convex optimization problems in which a convex function is minimized\n",
    "over the intersection of an affine subspace and a convex cone.\n",
    "An example of a conic-form minimization problems, in the primal form is:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "& \\min_{x \\in \\mathbb{R}^n} & a_0^T x + b_0 \\\\\n",
    "& \\;\\;\\text{s.t.} & A_i x + b_i & \\in \\mathcal{C}_i & i = 1 \\ldots m\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding dual problem is:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "& \\max_{y_1, \\ldots, y_m} & -\\sum_{i=1}^m b_i^T y_i + b_0 \\\\\n",
    "& \\;\\;\\text{s.t.} & a_0 - \\sum_{i=1}^m A_i^T y_i & = 0 \\\\\n",
    "& & y_i & \\in \\mathcal{C}_i^* & i = 1 \\ldots m\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where each $\\mathcal{C}_i$ is a closed convex cone and $\\mathcal{C}_i^*$ is its dual cone."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some of the Types of Cones Supported by JuMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP\n",
    "using ECOS\n",
    "using LinearAlgebra\n",
    "using Random\n",
    "\n",
    "Random.seed!(1234);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By this point we have used quite a few different solvers.\n",
    "To find out all the different solvers and their supported problem types, check out the\n",
    "[solver table](http://www.juliaopt.org/JuMP.jl/v0.19.0/installation/#Getting-Solvers-1) in the docs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second-Order Cone\n",
    "The Second-Order Cone (or Lorentz Cone) of dimension $n$ is of the form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "Q^n = \\{ (t,x) \\in \\mathbb{R}^\\mbox{n} : t \\ge ||x||_2 \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Second-Order Cone rotated by $\\pi/4$ in the $(x_1,x_2)$ plane is called a Rotated Second-Order Cone.\n",
    "It is of the form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "Q_r^n = \\{ (t,u,x) \\in \\mathbb{R}^\\mbox{n} : 2tu \\ge ||x||_2^2, t,u \\ge 0 \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These cones are represented in JuMP using the MOI sets `SecondOrderCone` and `RotatedSecondOrderCone`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Euclidean Projection on a hyperplane\n",
    "For a given point $u_{0}$ and a set $K$, we refer to any point $u \\in K$\n",
    "which is closest to $u_{0}$ as a projection of $u_{0}$ on $K$.\n",
    "The projection of a point $u_{0}$ on a hyperplane $K = \\{u | p' \\cdot u = q\\}$ is given by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\min & ||u - u_{0}|| \\\\\n",
    "& \\;\\;\\text{s.t.} & p' \\cdot u = q  \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = rand(10)\n",
    "p = rand(10)\n",
    "q = rand();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can model the above problem as the following conic program:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\min & t \\\\\n",
    "& \\;\\;\\text{s.t.} & p' \\cdot u = q \\\\\n",
    "& & (t, u - u_{0}) \\in Q^{n+1}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On comparing this with the primal form of a conic problem we saw above,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& x = (t , u) &\\\\\n",
    "& a_0 = e_1 &\\\\\n",
    "& b_0 = 0 &\\\\\n",
    "& A_1 = (0, p) &\\\\\n",
    "& b_1 = -q &\\\\\n",
    "& C_1 = \\mathbb{R}_- &\\\\\n",
    "& A_2 = 1 &\\\\\n",
    "& b_2 = -(0, u_0) &\\\\\n",
    "& C_2 = Q^{n+1} &\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we can obtain the dual problem as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\max & y_1 + (0, u_0)^T y_2 \\\\\n",
    "& \\;\\;\\text{s.t.} & e_1 - (0,p)^T y_1 - y_2 = 0 \\\\\n",
    "& & y_1 \\in \\mathbb{R}_- \\\\\n",
    "& & y_2 \\in Q^{n+1}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +0.000e+00  -0.000e+00  +1e+01  4e-01  2e-04  1e+00  6e+00    ---    ---    1  1  - |  -  - \n",
      " 1  +9.352e-02  +9.474e-02  +7e-01  4e-02  2e-05  7e-02  4e-01  0.9308  3e-04   2  2  2 |  0  0\n",
      " 2  +7.130e-01  +7.349e-01  +7e-02  5e-03  2e-06  3e-02  5e-02  0.9427  7e-02   2  2  2 |  0  0\n",
      " 3  +7.410e-01  +7.413e-01  +8e-04  6e-05  2e-08  3e-04  6e-04  0.9890  1e-04   2  2  2 |  0  0\n",
      " 4  +7.413e-01  +7.413e-01  +8e-06  6e-07  2e-10  4e-06  6e-06  0.9890  1e-04   1  1  1 |  0  0\n",
      " 5  +7.413e-01  +7.413e-01  +9e-08  7e-09  2e-12  4e-08  7e-08  0.9890  1e-04   1  1  1 |  0  0\n",
      " 6  +7.413e-01  +7.413e-01  +1e-09  8e-11  2e-14  5e-10  8e-10  0.9890  1e-04   1  1  1 |  0  0\n",
      "\n",
      "OPTIMAL (within feastol=7.8e-11, reltol=1.4e-09, abstol=1.0e-09).\n",
      "Runtime: 0.000137 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model = Model(optimizer_with_attributes(ECOS.Optimizer, \"printlevel\" => 0))\n",
    "@variable(model, u[1:10])\n",
    "@variable(model, t)\n",
    "@objective(model, Min, t)\n",
    "@constraint(model, [t, (u - u0)...] in SecondOrderCone())\n",
    "@constraint(model, u' * p == q)\n",
    "optimize!(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective_value(model) = 0.7413476010760025\n",
      "value.(u) = [0.3124638243936739, 0.7621182340845496, 0.5377408631527775, 0.049622730071534914, 0.516584856932142, 0.8058883387736702, 0.08216859863776912, 0.019039831335664947, 0.22253665470728054, 0.21813451035590134]\n"
     ]
    }
   ],
   "source": [
    "@show objective_value(model);\n",
    "@show value.(u);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +0.000e+00  -0.000e+00  +3e+00  4e-01  4e-01  1e+00  1e+00    ---    ---    1  1  - |  -  - \n",
      " 1  -4.087e+00  -1.050e+00  +4e-01  2e-01  3e-01  4e+00  2e-01  0.8614  3e-03   1  1  1 |  0  0\n",
      " 2  -1.100e+00  -1.173e+00  +1e-01  7e-02  3e-02  3e-02  5e-02  0.8459  1e-01   2  2  2 |  0  0\n",
      " 3  -7.558e-01  -7.447e-01  +5e-03  3e-03  1e-03  1e-02  3e-03  0.9794  4e-02   2  2  2 |  0  0\n",
      " 4  -7.415e-01  -7.414e-01  +5e-05  3e-05  1e-05  2e-04  3e-05  0.9890  1e-04   2  1  1 |  0  0\n",
      " 5  -7.413e-01  -7.413e-01  +6e-07  3e-07  1e-07  2e-06  4e-07  0.9890  1e-04   2  1  1 |  0  0\n",
      " 6  -7.413e-01  -7.413e-01  +7e-09  3e-09  2e-09  2e-08  4e-09  0.9890  1e-04   3  1  1 |  0  0\n",
      "\n",
      "OPTIMAL (within feastol=3.4e-09, reltol=8.9e-09, abstol=6.6e-09).\n",
      "Runtime: 0.000098 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "e1 = [1, zeros(10)...]\n",
    "dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, \"printlevel\" => 0))\n",
    "@variable(dual_model, y1 <= 0)\n",
    "@variable(dual_model, y2[1:11])\n",
    "@objective(dual_model, Max, q * y1 + dot(vcat(0, u0), y2))\n",
    "@constraint(dual_model, e1 - [0, p...] .* y1 - y2 .== 0)\n",
    "@constraint(dual_model, y2 in SecondOrderCone())\n",
    "optimize!(dual_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective_value(dual_model) = 0.7413476212168523\n"
     ]
    }
   ],
   "source": [
    "@show objective_value(dual_model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also have an equivalent formulation using a Rotated Second-Order Cone:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\min & t \\\\\n",
    "& \\;\\;\\text{s.t.} & p' \\cdot u = q \\\\\n",
    "& & (t, 1/2, u - u_{0})\\in Q_r^{n+2}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +0.000e+00  -3.294e-02  +9e+00  4e-01  2e-02  1e+00  5e+00    ---    ---    1  1  - |  -  - \n",
      " 1  +4.307e-03  +4.596e-03  +5e-01  3e-02  1e-03  5e-02  3e-01  0.9467  2e-04   1  1  1 |  0  0\n",
      " 2  +4.232e-01  +5.299e-01  +1e-01  2e-02  5e-04  1e-01  9e-02  0.9251  3e-01   2  2  2 |  0  0\n",
      " 3  +5.336e-01  +5.401e-01  +4e-03  8e-04  2e-05  8e-03  4e-03  0.9590  1e-04   2  1  2 |  0  0\n",
      " 4  +5.469e-01  +5.486e-01  +4e-04  1e-04  3e-06  2e-03  5e-04  0.9890  1e-01   2  2  2 |  0  0\n",
      " 5  +5.495e-01  +5.496e-01  +7e-06  3e-06  6e-08  4e-05  1e-05  0.9811  6e-04   2  2  2 |  0  0\n",
      " 6  +5.496e-01  +5.496e-01  +8e-08  3e-08  7e-10  4e-07  1e-07  0.9890  1e-04   2  1  1 |  0  0\n",
      " 7  +5.496e-01  +5.496e-01  +9e-10  3e-10  7e-12  5e-09  1e-09  0.9890  1e-04   2  1  1 |  0  0\n",
      "\n",
      "OPTIMAL (within feastol=3.2e-10, reltol=1.6e-09, abstol=8.7e-10).\n",
      "Runtime: 0.000112 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model = Model(optimizer_with_attributes(ECOS.Optimizer, \"printlevel\" => 0))\n",
    "@variable(model, u[1:10])\n",
    "@variable(model, t)\n",
    "@objective(model, Min, t)\n",
    "@constraint(model, [t, 0.5, (u - u0)...] in RotatedSecondOrderCone())\n",
    "@constraint(model, u' * p == q)\n",
    "optimize!(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "value.(u) = [0.3124638244248018, 0.7621182341177366, 0.5377408631710233, 0.049622730012400196, 0.5165848570203785, 0.8058883387998751, 0.08216859862006304, 0.019039831283705344, 0.22253665471226453, 0.2181345103716297]\n"
     ]
    }
   ],
   "source": [
    "@show value.(u);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference here is that the objective in the case of the Second-Order Cone is $||u - u_{0}||_2$,\n",
    "while in the case of a Rotated Second-Order Cone is $||u - u_{0}||_2^2$.\n",
    "However, the value of x is the same for both."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponential Cone"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An Exponential Cone is a set of the form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "K_{exp} = \\{ (x,y,z) \\in \\mathbb{R}^3 : y \\exp (x/y) \\le z, y > 0 \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is represented in JuMP using the MOI set `ExponentialCone`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Entropy Maximization\n",
    "As the name suggests, the entropy maximization problem consists of maximizing the entropy function,\n",
    "$H(x) = -x\\log{x}$ subject to linear inequality constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\max & - \\sum_{i=1}^n x_i \\log x_i \\\\\n",
    "& \\;\\;\\text{s.t.} & \\mathbf{1}' x = 1 \\\\\n",
    "& & Ax \\leq b\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can model this problem using an exponential cone by using the following transformation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "t\\leq -x\\log{x} \\iff t\\leq x\\log(1/x)  \\iff (t, x, 1) \\in K_{exp}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, our problem becomes,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\max & 1^Tt \\\\\n",
    "& \\;\\;\\text{s.t.} & Ax \\leq b \\\\\n",
    "& & 1^T x = 1 \\\\\n",
    "& & (t_i, x_i, 1) \\in K_{exp} && \\forall i = 1 \\ldots n \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n",
      "\n",
      "It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT\n",
      " 0  +0.000e+00  -1.380e+01  +6e+01  1e+00  6e-01  1e+00  1e+00    ---    ---    0  0  - |  -  - \n",
      " 1  -8.722e+00  -1.340e+01  +1e+01  3e-01  2e-01  4e-01  3e-01  0.7833  6e-02   1  1  1 |  2  1\n",
      " 2  -4.528e+00  -6.273e+00  +4e+00  1e-01  7e-02  2e-01  7e-02  0.9791  3e-01   1  1  1 |  5  0\n",
      " 3  -3.509e+00  -3.976e+00  +1e+00  4e-02  2e-02  5e-02  2e-02  0.7620  3e-02   1  1  1 |  1  1\n",
      " 4  -3.012e+00  -3.217e+00  +4e-01  2e-02  8e-03  2e-02  7e-03  0.6266  5e-02   1  1  1 |  2  2\n",
      " 5  -2.764e+00  -2.817e+00  +8e-02  5e-03  2e-03  5e-03  1e-03  0.9791  2e-01   1  1  1 |  4  0\n",
      " 6  -2.733e+00  -2.755e+00  +3e-02  2e-03  9e-04  2e-03  6e-04  0.7833  2e-01   1  1  1 |  4  1\n",
      " 7  -2.714e+00  -2.720e+00  +8e-03  6e-04  2e-04  6e-04  1e-04  0.7833  5e-02   1  1  1 |  2  1\n",
      " 8  -2.711e+00  -2.713e+00  +3e-03  2e-04  9e-05  2e-04  6e-05  0.6266  6e-02   1  1  1 |  2  2\n",
      " 9  -2.709e+00  -2.709e+00  +7e-04  5e-05  2e-05  5e-05  1e-05  0.7833  1e-02   2  1  1 |  1  1\n",
      "10  -2.708e+00  -2.708e+00  +3e-04  2e-05  9e-06  2e-05  5e-06  0.6266  5e-02   2  1  1 |  2  2\n",
      "11  -2.708e+00  -2.708e+00  +7e-05  5e-06  2e-06  5e-06  1e-06  0.7833  9e-03   1  1  0 |  1  1\n",
      "12  -2.708e+00  -2.708e+00  +3e-05  2e-06  8e-07  2e-06  5e-07  0.6266  5e-02   1  0  0 |  2  2\n",
      "13  -2.708e+00  -2.708e+00  +6e-06  4e-07  2e-07  4e-07  1e-07  0.7833  9e-03   1  0  0 |  1  1\n",
      "14  -2.708e+00  -2.708e+00  +2e-06  2e-07  7e-08  2e-07  4e-08  0.6266  5e-02   2  0  0 |  2  2\n",
      "15  -2.708e+00  -2.708e+00  +5e-07  4e-08  2e-08  4e-08  1e-08  0.7833  9e-03   1  0  0 |  1  1\n",
      "16  -2.708e+00  -2.708e+00  +2e-07  2e-08  6e-09  2e-08  4e-09  0.6266  5e-02   1  0  0 |  2  2\n",
      "17  -2.708e+00  -2.708e+00  +5e-08  4e-09  1e-09  3e-09  9e-10  0.7833  9e-03   1  0  0 |  1  1\n",
      "18  -2.708e+00  -2.708e+00  +2e-08  1e-09  6e-10  1e-09  4e-10  0.6266  5e-02   1  0  0 |  2  2\n",
      "\n",
      "OPTIMAL (within feastol=1.4e-09, reltol=7.2e-09, abstol=2.0e-08).\n",
      "Runtime: 0.000550 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "n = 15;\n",
    "m = 10;\n",
    "A = randn(m, n);\n",
    "b = rand(m, 1);\n",
    "\n",
    "model = Model(optimizer_with_attributes(ECOS.Optimizer, \"printlevel\" => 0))\n",
    "@variable(model, t[1:n])\n",
    "@variable(model, x[1:n])\n",
    "@objective(model, Max, sum(t))\n",
    "@constraint(model, sum(x) == 1)\n",
    "@constraint(model, A * x .<= b )\n",
    "# Cannot use the exponential cone directly in JuMP, hence we use MOI to specify the set.\n",
    "@constraint(model, con[i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone())\n",
    "\n",
    "optimize!(model);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective_value(model) = 2.708050216176813\n"
     ]
    }
   ],
   "source": [
    "@show objective_value(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Positive Semidefinite Cone\n",
    "The set of Positive Semidefinite Matrices of dimension $n$ form a cone in $\\mathbb{R}^n$.\n",
    "We write this set mathematically as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\mathcal{S}_{+}^n = \\{ X \\in \\mathcal{S}^n \\mid z^T X z \\geq 0, \\: \\forall z\\in \\mathbb{R}^n \\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A PSD cone is represented in JuMP using the MOI sets\n",
    "`PositiveSemidefiniteConeTriangle` (for upper triangle of a PSD matrix) and\n",
    "`PositiveSemidefiniteConeSquare` (for a complete PSD matrix).\n",
    "However, it is preferable to use the `PSDCone` shortcut as illustrated below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Largest Eigenvalue of a Symmetric Matrix\n",
    "Suppose $A$ has eigenvalues $\\lambda_{1} \\geq \\lambda_{2} \\ldots \\geq \\lambda_{n}$.\n",
    "Then the matrix $t I-A$ has eigenvalues $t-\\lambda_{1}, t-\\lambda_{2}, \\ldots, t-\\lambda_{n}$.\n",
    "Note that $t I-A$ is PSD exactly when all these eigenvalues are non-negative,\n",
    "and this happens for values $t \\geq \\lambda_{1} .$\n",
    "Thus, we can model the problem of finding the largest eigenvalue of a symmetric matrix as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\lambda_{1} = \\max t \\\\\n",
    "\\text { s.t. } t I-A \\succeq 0\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using SCS\n",
    "\n",
    "A = [3 2 4;\n",
    "     2 0 2;\n",
    "     4 2 3]\n",
    "\n",
    "model = Model(optimizer_with_attributes(SCS.Optimizer, \"verbose\" => 0))\n",
    "@variable(model, t)\n",
    "@objective(model, Min, t)\n",
    "@constraint(model, t .* Matrix{Float64}(I, 3, 3) - A in PSDCone())\n",
    "\n",
    "optimize!(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective_value(model) = 8.000000000000014\n"
     ]
    }
   ],
   "source": [
    "@show objective_value(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Cones and Functions\n",
    "For other cones supported by JuMP, check out the\n",
    "[MathOptInterface Manual](http://www.juliaopt.org/MathOptInterface.jl/dev/apimanual/#Standard-form-problem-1).\n",
    "A good resource for learning more about functions which can be modelled using cones is the\n",
    "MOSEK Modeling Cookbook[[1]](#c1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "<a id='c1'></a>\n",
    "1. MOSEK Modeling Cookbook — MOSEK Modeling Cookbook 3.1. Available at: https://docs.mosek.com/modeling-cookbook/index.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}