Originally Contributed by: Arpit Bhatia
This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP.
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 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:
The corresponding dual problem is:
where each $\mathcal{C}_i$ is a closed convex cone and $\mathcal{C}_i^*$ is its dual cone.
using JuMP
using ECOS
using LinearAlgebra
using Random
Random.seed!(1234);
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 in the docs.
The Second-Order Cone (or Lorentz Cone) of dimension $n$ is of the form:
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:
These cones are represented in JuMP using the MOI sets SecondOrderCone
and RotatedSecondOrderCone
.
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
u0 = rand(10)
p = rand(10)
q = rand();
We can model the above problem as the following conic program:
On comparing this with the primal form of a conic problem we saw above,
Thus, we can obtain the dual problem as:
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.000137 seconds.
@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]
e1 = [1, zeros(10)...]
dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(dual_model, y1 <= 0)
@variable(dual_model, y2[1:11])
@objective(dual_model, Max, q * y1 + dot(vcat(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 +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.000098 seconds.
@show objective_value(dual_model);
objective_value(dual_model) = 0.7413476212168523
We can also have an equivalent formulation using a Rotated Second-Order Cone:
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.000112 seconds.
@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.
An Exponential Cone is a set of the form:
It is represented in JuMP using the MOI set ExponentialCone
.
As the name suggests, the entropy maximization problem consists of maximizing the entropy function, $H(x) = -x\log{x}$ subject to linear inequality constraints.
We can model this problem using an exponential cone by using the following transformation:
Thus, our problem becomes,
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)
@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] 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.380e+01 +6e+01 1e+00 6e-01 1e+00 1e+00 --- --- 0 0 - | - - 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 OPTIMAL (within feastol=1.4e-09, reltol=7.2e-09, abstol=2.0e-08). Runtime: 0.000550 seconds.
@show objective_value(model);
objective_value(model) = 2.708050216176813
The set of Positive Semidefinite Matrices of dimension $n$ form a cone in $\mathbb{R}^n$. We write this set mathematically as
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.
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:
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)
@show objective_value(model);
objective_value(model) = 8.000000000000014
For other cones supported by JuMP, check out the MathOptInterface Manual. A good resource for learning more about functions which can be modelled using cones is the MOSEK Modeling Cookbook[1].