# Conic Programming

**Originally Contributed by**: Arpit Bhatia

This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP.

## What is a Cone?
A subset $C$ of a vector space $V$ is a cone if $\forall x \in C$ and positive scalars $\alpha$,
the product $\alpha x \in C$. A cone C is a convex cone if $\alpha x + \beta y \in C$,
for any positive scalars $\alpha, \beta$, and any $x, y \in C$.

## Conic Programming
Conic programming problems are convex optimization problems in which a convex function is minimized
over the intersection of an affine subspace and a convex cone.
An example of a conic-form minimization problems, in the primal form is:

$$
\begin{align}
& \min_{x \in \mathbb{R}^n} & a_0^T x + b_0 \\
& \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m
\end{align}
$$

The corresponding dual problem is:

$$
\begin{align}
& \max_{y_1, \ldots, y_m} & -\sum_{i=1}^m b_i^T y_i + b_0 \\
& \;\;\text{s.t.} & a_0 - \sum_{i=1}^m A_i^T y_i & = 0 \\
& & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m
\end{align}
$$

where each $\mathcal{C}_i$ is a closed convex cone and $\mathcal{C}_i^*$ is its dual cone.

## Some of the Types of Cones Supported by JuMP


```julia
using JuMP
using ECOS
using LinearAlgebra
using Random

Random.seed!(1234);
```

    ┌ Info: Precompiling JuMP [4076af6c-e467-56ae-b986-b466b2749572]
    └ @ Base loading.jl:1260
    ┌ Info: Precompiling ECOS [e2685f51-7e38-5353-a97d-a921fd2c8199]
    └ @ Base loading.jl:1260


By this point we have used quite a few different solvers.
To find out all the different solvers and their supported problem types, check out the
[solver table](http://www.juliaopt.org/JuMP.jl/v0.19.0/installation/#Getting-Solvers-1) in the docs.

### Second-Order Cone
The Second-Order Cone (or Lorentz Cone) of dimension $n$ is of the form:

$$
Q^n = \{ (t,x) \in \mathbb{R}^\mbox{n} : t \ge ||x||_2 \}
$$

A Second-Order Cone rotated by $\pi/4$ in the $(x_1,x_2)$ plane is called a Rotated Second-Order Cone.
It is of the form:

$$
Q_r^n = \{ (t,u,x) \in \mathbb{R}^\mbox{n} : 2tu \ge ||x||_2^2, t,u \ge 0 \}
$$

These cones are represented in JuMP using the MOI sets `SecondOrderCone` and `RotatedSecondOrderCone`.

#### Example: Euclidean Projection on a hyperplane
For a given point $u_{0}$ and a set $K$, we refer to any point $u \in K$
which is closest to $u_{0}$ as a projection of $u_{0}$ on $K$.
The projection of a point $u_{0}$ on a hyperplane $K = \{u | p' \cdot u = q\}$ is given by

$$
\begin{align*}
& \min & ||u - u_{0}|| \\
& \;\;\text{s.t.} & p' \cdot u = q  \\
\end{align*}
$$


```julia
u0 = rand(10)
p = rand(10)
q = rand();
```

We can model the above problem as the following conic program:

$$
\begin{align*}
& \min & t \\
& \;\;\text{s.t.} & p' \cdot u = q \\
& & (t, u - u_{0}) \in Q^{n+1}
\end{align*}
$$

On comparing this with the primal form of a conic problem we saw above,

$$
\begin{align*}
& x = (t , u) &\\
& a_0 = e_1 &\\
& b_0 = 0 &\\
& A_1 = (0, p) &\\
& b_1 = -q &\\
& C_1 = \mathbb{R}_- &\\
& A_2 = 1 &\\
& b_2 = -(0, u_0) &\\
& C_2 = Q^{n+1} &
\end{align*}
$$

Thus, we can obtain the dual problem as:

$$
\begin{align*}
& \max & y_1 + (0, u_0)^T y_2 \\
& \;\;\text{s.t.} & e_1 - (0,p)^T y_1 - y_2 = 0 \\
& & y_1 \in \mathbb{R}_- \\
& & y_2 \in Q^{n+1}
\end{align*}
$$


```julia
model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, (u - u0)...] in SecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)
```

    
    ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
    
    It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
     0  +0.000e+00  -0.000e+00  +1e+01  4e-01  2e-04  1e+00  6e+00    ---    ---    1  1  - |  -  - 
     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
     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
     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
     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
     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
     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
    
    OPTIMAL (within feastol=7.8e-11, reltol=1.4e-09, abstol=1.0e-09).
    Runtime: 0.000191 seconds.
    



```julia
@show objective_value(model);
@show value.(u);
```

    objective_value(model) = 0.7413476010760025
    value.(u) = [0.3124638243936739, 0.7621182340845496, 0.5377408631527775, 0.049622730071534914, 0.516584856932142, 0.8058883387736702, 0.08216859863776912, 0.019039831335664947, 0.22253665470728054, 0.21813451035590134]



```julia
e1 = [1.0, zeros(10)...]
dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(dual_model, y1 <= 0.0)
@variable(dual_model, y2[1:11])
@objective(dual_model, Max, q * y1 + dot(vcat(0.0, u0), y2))
@constraint(dual_model, e1 - [0.0, p...] .* y1 - y2 .== 0.0)
@constraint(dual_model, y2 in SecondOrderCone())
optimize!(dual_model)
```

    
    ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
    
    It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
     0  +0.000e+00  -0.000e+00  +3e+00  4e-01  4e-01  1e+00  1e+00    ---    ---    1  1  - |  -  - 
     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
     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
     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
     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
     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
     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
    
    OPTIMAL (within feastol=3.4e-09, reltol=8.9e-09, abstol=6.6e-09).
    Runtime: 0.000134 seconds.
    



```julia
@show objective_value(dual_model);
```

    objective_value(dual_model) = 0.7413476212168523


We can also have an equivalent formulation using a Rotated Second-Order Cone:

$$
\begin{align*}
& \min & t \\
& \;\;\text{s.t.} & p' \cdot u = q \\
& & (t, 1/2, u - u_{0})\in Q_r^{n+2}
\end{align*}
$$


```julia
model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, 0.5, (u - u0)...] in RotatedSecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)
```

    
    ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
    
    It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
     0  +0.000e+00  -3.294e-02  +9e+00  4e-01  2e-02  1e+00  5e+00    ---    ---    1  1  - |  -  - 
     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
     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
     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
     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
     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
     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
     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
    
    OPTIMAL (within feastol=3.2e-10, reltol=1.6e-09, abstol=8.7e-10).
    Runtime: 0.000134 seconds.
    



```julia
@show value.(u);
```

    value.(u) = [0.3124638244248018, 0.7621182341177366, 0.5377408631710233, 0.049622730012400196, 0.5165848570203785, 0.8058883387998751, 0.08216859862006304, 0.019039831283705344, 0.22253665471226453, 0.2181345103716297]


The difference here is that the objective in the case of the Second-Order Cone is $||u - u_{0}||_2$,
while in the case of a Rotated Second-Order Cone is $||u - u_{0}||_2^2$.
However, the value of x is the same for both.

### Exponential Cone

An Exponential Cone is a set of the form:

$$
K_{exp} = \{ (x,y,z) \in \mathbb{R}^3 : y \exp (x/y) \le z, y > 0 \}
$$

It is represented in JuMP using the MOI set `ExponentialCone`.

#### Example: Entropy Maximization
As the name suggests, the entropy maximization problem consists of maximizing the entropy function,
$H(x) = -x\log{x}$ subject to linear inequality constraints.

$$
\begin{align*}
& \max & - \sum_{i=1}^n x_i \log x_i \\
& \;\;\text{s.t.} & \mathbf{1}' x = 1 \\
& & Ax \leq b
\end{align*}
$$

We can model this problem using an exponential cone by using the following transformation:

$$
t\leq -x\log{x} \iff t\leq x\log(1/x)  \iff (t, x, 1) \in K_{exp}
$$

Thus, our problem becomes,

$$
\begin{align*}
& \max & 1^Tt \\
& \;\;\text{s.t.} & Ax \leq b \\
& & 1^T x = 1 \\
& & (t_i, x_i, 1) \in K_{exp} && \forall i = 1 \ldots n \\
\end{align*}
$$


```julia
n = 15;
m = 10;
A = randn(m, n);
b = rand(m, 1);

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1.0)
@constraint(model, A * x .<= b )
# Cannot use the exponential cone directly in JuMP, hence we use MOI to specify the set.
@constraint(model, con[i = 1:n], [t[i], x[i], 1.0] in MOI.ExponentialCone())

optimize!(model);
```

    
    ECOS 2.0.5 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
    
    It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
     0  +0.000e+00  -1.397e+01  +6e+01  1e+00  5e-01  1e+00  1e+00    ---    ---    0  0  - |  -  - 
     1  -8.340e+00  -1.331e+01  +2e+01  3e-01  2e-01  4e-01  3e-01  0.7833  8e-02   1  1  1 |  2  1
     2  -4.641e+00  -6.976e+00  +5e+00  2e-01  9e-02  2e-01  1e-01  0.9791  4e-01   1  1  1 |  6  0
     3  -3.703e+00  -4.235e+00  +1e+00  4e-02  2e-02  5e-02  2e-02  0.7833  9e-03   1  1  1 |  1  1
     4  -3.051e+00  -3.294e+00  +4e-01  2e-02  9e-03  3e-02  8e-03  0.6266  5e-02   1  1  1 |  2  2
     5  -2.789e+00  -2.868e+00  +1e-01  8e-03  3e-03  8e-03  2e-03  0.7833  1e-01   1  1  1 |  3  1
     6  -2.741e+00  -2.775e+00  +5e-02  3e-03  1e-03  3e-03  1e-03  0.9791  4e-01   1  1  1 |  6  0
     7  -2.709e+00  -2.717e+00  +1e-02  8e-04  3e-04  8e-04  2e-04  0.7833  9e-03   1  1  1 |  1  1
     8  -2.704e+00  -2.706e+00  +3e-03  2e-04  7e-05  2e-04  5e-05  0.7833  1e-02   2  1  1 |  1  1
     9  -2.703e+00  -2.703e+00  +6e-04  5e-05  2e-05  5e-05  1e-05  0.7833  5e-02   1  1  1 |  2  1
    10  -2.702e+00  -2.703e+00  +3e-04  2e-05  7e-06  2e-05  5e-06  0.6266  5e-02   2  1  1 |  2  2
    11  -2.702e+00  -2.702e+00  +6e-05  4e-06  2e-06  4e-06  1e-06  0.7833  1e-02   1  0  0 |  1  1
    12  -2.702e+00  -2.702e+00  +2e-05  1e-06  4e-07  1e-06  3e-07  0.7833  5e-02   2  0  0 |  2  1
    13  -2.702e+00  -2.702e+00  +6e-06  4e-07  2e-07  4e-07  1e-07  0.6266  5e-02   2  0  0 |  2  2
    14  -2.702e+00  -2.702e+00  +1e-06  1e-07  4e-08  1e-07  3e-08  0.7833  1e-02   1  0  0 |  1  1
    15  -2.702e+00  -2.702e+00  +3e-07  2e-08  1e-08  2e-08  6e-09  0.7833  5e-02   1  0  0 |  2  1
    16  -2.702e+00  -2.702e+00  +1e-07  9e-09  4e-09  1e-08  2e-09  0.6266  5e-02   1  0  1 |  2  2
    17  -2.702e+00  -2.702e+00  +3e-08  2e-09  9e-10  2e-09  6e-10  0.7833  9e-03   0  0  0 |  1  1
    18  -2.702e+00  -2.702e+00  +7e-09  5e-10  2e-10  5e-10  1e-10  0.7833  9e-03   0  0  0 |  1  1
    
    OPTIMAL (within feastol=5.1e-10, reltol=2.7e-09, abstol=7.4e-09).
    Runtime: 0.000841 seconds.
    



```julia
@show objective_value(model);
```

    objective_value(model) = 2.7021992720612427


### Positive Semidefinite Cone
The set of Positive Semidefinite Matrices of dimension $n$ form a cone in $\mathbb{R}^n$.
We write this set mathematically as

$$
\mathcal{S}_{+}^n = \{ X \in \mathcal{S}^n \mid z^T X z \geq 0, \: \forall z\in \mathbb{R}^n \}.
$$

A PSD cone is represented in JuMP using the MOI sets
`PositiveSemidefiniteConeTriangle` (for upper triangle of a PSD matrix) and
`PositiveSemidefiniteConeSquare` (for a complete PSD matrix).
However, it is preferable to use the `PSDCone` shortcut as illustrated below.

#### Example: Largest Eigenvalue of a Symmetric Matrix
Suppose $A$ has eigenvalues $\lambda_{1} \geq \lambda_{2} \ldots \geq \lambda_{n}$.
Then the matrix $t I-A$ has eigenvalues $t-\lambda_{1}, t-\lambda_{2}, \ldots, t-\lambda_{n}$.
Note that $t I-A$ is PSD exactly when all these eigenvalues are non-negative,
and this happens for values $t \geq \lambda_{1} .$
Thus, we can model the problem of finding the largest eigenvalue of a symmetric matrix as:

$$
\begin{align*}
\lambda_{1} = \max t \\
\text { s.t. } t I-A \succeq 0
\end{align*}
$$


```julia
using LinearAlgebra
using SCS

A = [3 2 4;
     2 0 2;
     4 2 3]

model = Model(optimizer_with_attributes(SCS.Optimizer, "verbose" => 0))
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* Matrix{Float64}(I, 3, 3) - A in PSDCone())

optimize!(model)
```

    ┌ Info: Precompiling SCS [c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13]
    └ @ Base loading.jl:1260



```julia
@show objective_value(model);
```

    objective_value(model) = 8.000000000000014


## Other Cones and Functions
For other cones supported by JuMP, check out the
[MathOptInterface Manual](http://www.juliaopt.org/MathOptInterface.jl/dev/apimanual/#Standard-form-problem-1).
A good resource for learning more about functions which can be modelled using cones is the
MOSEK Modeling Cookbook[[1]](#c1).

### References
<a id='c1'></a>
1. MOSEK Modeling Cookbook — MOSEK Modeling Cookbook 3.1. Available at: https://docs.mosek.com/modeling-cookbook/index.html.


```julia

```

7021c1118968be946fa061a7db9fbf01bbce7b87
ipynb
Jupyter Notebook
notebook/optimization_concepts/conic_programming.ipynb
mtanneau/JuMPTutorials.jl
f91839c6b127966c3cb17e8971cc95fe3adfa814
MIT
2020-06-04T22:17:54.000Z
2020-06-04T22:17:54.000Z
notebook/optimization_concepts/conic_programming.ipynb
mtanneau/JuMPTutorials.jl
f91839c6b127966c3cb17e8971cc95fe3adfa814
MIT
2020-08-04T18:36:58.000Z
2020-08-04T18:36:58.000Z
notebook/optimization_concepts/conic_programming.ipynb
mtanneau/JuMPTutorials.jl
f91839c6b127966c3cb17e8971cc95fe3adfa814
MIT
2020-05-18T01:29:51.000Z
2020-05-18T01:29:51.000Z
Qwen/Qwen-72B
1. YES
2. YES