# 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
```

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 Lorenz 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(with_optimizer(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  +1.534e-01  +1.559e-01  +9e-01  5e-02  2e-05  9e-02  5e-01  0.9177  5e-04   2  1  1 |  0  0
     2  +9.409e-01  +9.665e-01  +7e-02  6e-03  2e-06  3e-02  5e-02  0.9533  6e-02   2  2  2 |  0  0
     3  +9.743e-01  +9.746e-01  +8e-04  6e-05  2e-08  4e-04  6e-04  0.9890  1e-04   1  1  1 |  0  0
     4  +9.747e-01  +9.747e-01  +9e-06  7e-07  2e-10  4e-06  7e-06  0.9890  1e-04   1  1  1 |  0  0
     5  +9.747e-01  +9.747e-01  +1e-07  8e-09  2e-12  5e-08  7e-08  0.9890  1e-04   1  1  1 |  0  0
     6  +9.747e-01  +9.747e-01  +1e-09  8e-11  2e-14  5e-10  8e-10  0.9890  1e-04   2  1  1 |  0  0
    
    OPTIMAL (within feastol=8.4e-11, reltol=1.1e-09, abstol=1.1e-09).
    Runtime: 0.000177 seconds.
    



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

    objective_value(model) = 0.9747011153011627
    value.(u) = [0.659213, 0.309677, 0.511505, 0.285757, 0.146255, 0.241544, -0.352419, -0.0583324, -0.107894, 0.511019]



```julia
e1 = [1, zeros(10)...]
dual_model = Model(with_optimizer(ECOS.Optimizer, printlevel = 0))
@variable(dual_model, y1 <= 0)
@variable(dual_model, y2[1:11])
@objective(dual_model, Max, q * y1 + [0, u0...]' * y2)
@constraint(dual_model, e1 - [0, p...] .* y1 - y2 .== 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  +4e+00  4e-01  4e-01  1e+00  2e+00    ---    ---    1  1  - |  -  - 
     1  -5.671e+00  -1.407e+00  +4e-01  3e-01  3e-01  5e+00  2e-01  0.8790  2e-03   1  1  1 |  0  0
     2  -1.487e+00  -1.569e+00  +1e-01  9e-02  4e-02  3e-02  5e-02  0.8344  1e-01   2  1  1 |  0  0
     3  -9.975e-01  -9.795e-01  +5e-03  3e-03  2e-03  2e-02  3e-03  0.9786  4e-02   2  2  2 |  0  0
     4  -9.750e-01  -9.748e-01  +6e-05  4e-05  2e-05  2e-04  3e-05  0.9890  1e-04   2  1  1 |  0  0
     5  -9.747e-01  -9.747e-01  +6e-07  4e-07  2e-07  3e-06  4e-07  0.9890  1e-04   2  1  1 |  0  0
     6  -9.747e-01  -9.747e-01  +7e-09  4e-09  2e-09  3e-08  4e-09  0.9890  1e-04   1  1  1 |  0  0
    
    OPTIMAL (within feastol=4.4e-09, reltol=7.1e-09, abstol=6.9e-09).
    Runtime: 0.000099 seconds.
    



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

    objective_value(dual_model) = 0.9747011467702701


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(with_optimizer(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  -7.846e-17  -3.254e-02  +9e+00  4e-01  2e-02  1e+00  5e+00    ---    ---    1  1  - |  -  - 
     1  +2.476e-02  +2.543e-02  +6e-01  4e-02  1e-03  7e-02  4e-01  0.9304  2e-04   1  1  1 |  0  0
     2  +6.892e-01  +9.125e-01  +3e-01  5e-02  1e-03  3e-01  2e-01  0.8037  4e-01   1  2  2 |  0  0
     3  +9.190e-01  +9.297e-01  +7e-03  2e-03  3e-05  1e-02  5e-03  0.9731  1e-04   1  1  1 |  0  0
     4  +9.371e-01  +9.412e-01  +2e-03  4e-04  9e-06  5e-03  1e-03  0.7864  7e-02   1  1  1 |  0  0
     5  +9.450e-01  +9.479e-01  +7e-04  2e-04  5e-06  3e-03  6e-04  0.7406  3e-01   1  1  1 |  0  0
     6  +9.483e-01  +9.490e-01  +2e-04  6e-05  1e-06  8e-04  2e-04  0.7802  5e-02   1  1  1 |  0  0
     7  +9.500e-01  +9.500e-01  +1e-05  3e-06  6e-08  3e-05  8e-06  0.9715  3e-02   1  1  1 |  0  0
     8  +9.500e-01  +9.500e-01  +1e-06  3e-07  6e-09  3e-06  8e-07  0.9062  3e-03   1  1  1 |  0  0
     9  +9.500e-01  +9.500e-01  +4e-08  1e-08  3e-10  2e-07  4e-08  0.9890  4e-02   1  1  1 |  0  0
    10  +9.500e-01  +9.500e-01  +1e-09  4e-10  9e-12  6e-09  1e-09  0.9702  1e-04   2  1  1 |  0  0
    
    OPTIMAL (within feastol=4.3e-10, reltol=1.2e-09, abstol=1.2e-09).
    Runtime: 0.000125 seconds.
    



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

    value.(u) = [0.659213, 0.309677, 0.511505, 0.285757, 0.146255, 0.241544, -0.352419, -0.0583324, -0.107894, 0.511019]


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 (1, x, t) \in K_{exp}
$$

Thus, our problem becomes,

$$
\begin{align*}
& \max & 1^Tt \\
& \;\;\text{s.t.} & Ax \leq b \\
& & 1^T x = 1 \\
& & (1, x_i, t_i) \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(with_optimizer(ECOS.Optimizer, printlevel = 0))
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1)
@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], [1, x[i], t[i]] 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  +6.662e+00  +6e+01  1e+00  8e-01  1e+00  1e+00    ---    ---    0  0  - |  -  - 
     1  +5.698e+01  +2.834e+02  +2e+01  2e+01  2e+01  1e+02  4e-01  0.6305  7e-02   1  1  1 |  0  0
     2  -2.863e+01  +1.302e+01  +1e+01  7e-01  6e-01  4e+01  3e-01  0.7833  6e-01   1  1  1 |  0  1
     3  -8.836e+01  +3.194e+01  +5e+00  9e-01  7e-01  1e+02  1e-01  0.5949  9e-02   1  2  2 |  2  2
     4  -5.666e+02  +1.039e+02  +2e+00  1e+00  1e+00  7e+02  6e-02  0.7833  2e-01   1  1  1 |  3  1
     5  -4.486e+03  +3.220e+02  +7e-01  3e+00  2e+00  5e+03  1e-02  0.9791  2e-01   1  1  1 |  4  0
     6  -1.438e+04  +3.825e+02  +1e-01  3e+00  2e+00  1e+04  3e-03  0.7815  1e-02   1  1  1 |  1  1
     7  -5.045e+04  -5.628e+02  +3e-02  2e+00  2e+00  5e+04  6e-04  0.7826  1e-02   2  1  1 |  1  1
     8  -1.894e+05  -7.687e+03  +6e-03  2e+00  1e+00  2e+05  1e-04  0.7833  9e-03   2  1  1 |  1  1
     9  -7.396e+05  -4.634e+04  +1e-03  2e+00  1e+00  7e+05  3e-05  0.7833  9e-03   1  1  1 |  1  1
    10  -2.955e+06  -2.334e+05  +3e-04  2e+00  1e+00  3e+06  5e-06  0.7833  9e-03   1  1  1 |  1  1
    11  -1.197e+07  -1.090e+06  +5e-05  2e+00  1e+00  1e+07  1e-06  0.7833  9e-03   2  1  1 |  1  1
    12  -1.110e+08  -1.198e+07  +7e-06  3e+00  2e+00  1e+08  1e-07  0.9791  1e-01   2  0  0 |  3  0
    13  -3.626e+08  -4.717e+07  +1e-06  3e+00  2e+00  3e+08  3e-08  0.7833  1e-02   2  0  0 |  1  1
    14  -1.728e+09  -2.822e+08  +4e-07  6e+00  3e+00  1e+09  7e-09  0.9791  3e-01   1  0  0 |  5  0
    15  -6.699e+09  -1.226e+09  +6e-08  6e+00  2e+00  5e+09  1e-09  0.7833  1e-02   1  0  0 |  1  1
    16  -6.766e+10  -1.311e+10  +2e-08  2e+01  6e+00  5e+10  3e-10  0.9791  3e-01   1  0  0 |  5  0
    
    UNBOUNDED (within feastol=3.5e-09).
    Runtime: 0.000772 seconds.
    



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

    objective_value(model) = 6.765923931549602e10


### 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 prefferable 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 CSDP

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

model = Model(with_optimizer(CSDP.Optimizer, printlevel = 0))
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* Matrix{Float64}(I, 3, 3) - A in PSDCone())

optimize!(model)
```


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

    objective_value(model) = 8.000000000000036


## 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.

8efe0d62ae61a957ed193b978ae4c84de3eceff3
ipynb
Jupyter Notebook
notebook/optimization_concepts/conic_programming.ipynb
carlosal1015/JuMPTutorials.jl
4d9a86ea310ecc7a22de7f14b783dbd218e4b612
MIT
notebook/optimization_concepts/conic_programming.ipynb
carlosal1015/JuMPTutorials.jl
4d9a86ea310ecc7a22de7f14b783dbd218e4b612
MIT
notebook/optimization_concepts/conic_programming.ipynb
carlosal1015/JuMPTutorials.jl
4d9a86ea310ecc7a22de7f14b783dbd218e4b612
MIT
Qwen/Qwen-72B
1. YES
2. YES