Tips and Tricks
|
Страница в процессе перевода. |
This tutorial was generated using Literate.jl. Download the source as a .jl file.
This tutorial was originally contributed by Arpit Bhatia.
This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP.
It uses the following packages:
using JuMP
import SCS
import LinearAlgebra|
This tutorial uses sets from MathOptInterface. By default, JuMP exports the |
```julia import MathOptInterface as MOI ```
|
A good resource for learning more about functions which can be modeled using cones is the MOSEK Modeling Cookbook. |
Background theory
A subset of a vector space is a cone if and positive scalars , the product .
A cone is a convex cone if , for any ], and any .
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 is a closed convex cone and is its dual cone.
Second-Order Cone
The SecondOrderCone (or Lorentz Cone) of dimension is a cone of the form:
It is most commonly used to represent the L2-norm of the vector :
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:3])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in SecondOrderCone())
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), value.(x)(0.5773503148525505, [0.3333333333638094, 0.33333333336381227, 0.33333333336381227])Rotated Second-Order Cone
A Second-Order Cone rotated by in the plane is called a RotatedSecondOrderCone. It is a cone of the form:
When u = 0.5, it represents the sum of squares of a vector :
data = [1.0, 2.0, 3.0, 4.0]
target = [0.45, 1.04, 1.51, 1.97]
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, θ)
@variable(model, t)
@expression(model, residuals, θ * data .- target)
@constraint(model, [t; 0.5; residuals] in RotatedSecondOrderCone())
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(θ), value(t)(0.4980000000000055, 0.004979422159620655)Exponential Cone
The MOI.ExponentialCone is a set of the form:
It can be used to model problems involving log and exp.
Exponential
To model , use (x, 1, z):
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x == 1.5)
@variable(model, z)
@objective(model, Min, z)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
value(z), exp(1.5)(4.481654619603649, 4.4816890703380645)Logarithm
To model , use (x, 1, z):
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x)
@variable(model, z == 1.5)
@objective(model, Max, x)
@constraint(model, [x, 1, z] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
value(x), log(1.5)(0.4055956586655662, 0.4054651081081644)Log-sum-exp
To model , use:
N = 3
x0 = rand(N)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[i = 1:N] == x0[i])
@variable(model, t)
@objective(model, Min, t)
@variable(model, u[1:N])
@constraint(model, sum(u) <= 1)
@constraint(model, [i = 1:N], [x[i] - t, 1, u[i]] in MOI.ExponentialCone())
optimize!(model)
value(t), log(sum(exp.(x0)))(1.4726849472138828, 1.472772274699489)Entropy
The entropy maximization problem consists of maximizing the entropy function, subject to linear inequality constraints.
We can model this problem using an exponential cone by using the following transformation:
Thus, our problem becomes,
m, n = 10, 15
A, b = randn(m, n), rand(m, 1)
model = Model(SCS.Optimizer)
set_silent(model)
@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)
@constraint(model, [i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone())
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)2.663206172203935The MOI.ExponentialCone has a dual, the MOI.DualExponentialCone, that offers an alternative formulation that can be more efficient for some formulations.
There is also the MOI.RelativeEntropyCone for explicitly encoding the relative entropy function
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x[1:n])
@objective(model, Max, -t)
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b)
@constraint(model, [t; ones(n); x] in MOI.RelativeEntropyCone(2n + 1))
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)2.6631736431749626PowerCone
The MOI.PowerCone is a three-dimensional set parameterized by a scalar value α. It has the form:
The power cone permits a number of reformulations. For example, when , we can model using the power cone with . Thus, to model with
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, x >= 1.5)
@constraint(model, [t, 1, x] in MOI.PowerCone(1 / 3))
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), value(x)(3.374754835774137, 1.49995522631532)The MOI.PowerCone has a dual, the MOI.DualPowerCone, that offers an alternative formulation that can be more efficient for some formulations.
P-Norm
The p-norm can be modeled using MOI.PowerCones. See the Mosek Modeling Cookbook for the derivation.
function p_norm(x::Vector, p)
N = length(x)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, r[1:N])
@variable(model, t)
@constraint(model, [i = 1:N], [r[i], t, x[i]] in MOI.PowerCone(1 / p))
@constraint(model, sum(r) == t)
@objective(model, Min, t)
optimize!(model)
@assert is_solved_and_feasible(model)
return value(t)
end
x = rand(5);
LinearAlgebra.norm(x, 4), p_norm(x, 4)(0.9316922467512209, 0.931687599495156)Positive Semidefinite Cone
The set of positive semidefinite matrices (PSD) of dimension form a cone in . 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.
Example: largest eigenvalue of a symmetric matrix
Suppose has eigenvalues . Then the matrix has eigenvalues . Note that is PSD exactly when all these eigenvalues are non-negative, and this happens for values . Thus, we can model the problem of finding the largest eigenvalue of a symmetric matrix as:
A = [3 2 4; 2 0 2; 4 2 3]
I = Matrix{Float64}(LinearAlgebra.I, 3, 3)
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* I - A in PSDCone())
optimize!(model)
@assert is_solved_and_feasible(model)
objective_value(model)8.000003377698668GeometricMeanCone
The MOI.GeometricMeanCone is a cone of the form:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, x[1:4])
@variable(model, t)
@constraint(model, sum(x) == 1)
@constraint(model, [t; x] in MOI.GeometricMeanCone(5))
optimize!(model)
value(t), value.(x)(0.0, [0.25000000110304854, 0.2500000011030519, 0.2500000011030504, 0.25000000110305026])RootDetCone
The MOI.RootDetConeSquare is a cone of the form:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; vec(X)] in MOI.RootDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))(2.236116364743881, 2.236067957402722)If X is symmetric, then you can use MOI.RootDetConeTriangle instead. This can be more efficient because the solver does not need to add additional constraints to ensure X is symmetric.
When forming the function, use triangle_vec to obtain the column-wise upper triangle of the matrix as a vector in the order that JuMP requires.
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; triangle_vec(X)] in MOI.RootDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
optimize!(model)
value(t), sqrt(LinearAlgebra.det(value.(X)))(2.2361792479821463, 2.2360679863317)LogDetCone
The MOI.LogDetConeSquare is a cone of the form:
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2])
@objective(model, Max, t)
@constraint(model, [t; u; vec(X)] in MOI.LogDetConeSquare(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))(1.497920445822091, 1.497866100640713)If X is symmetric, then you can use MOI.LogDetConeTriangle instead. This can be more efficient because the solver does not need to add additional constraints to ensure X is symmetric.
When forming the function, use triangle_vec to obtain the column-wise upper triangle of the matrix as a vector in the order that JuMP requires.
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, t)
@variable(model, u)
@variable(model, X[1:2, 1:2], Symmetric)
@objective(model, Max, t)
@constraint(model, [t; u; triangle_vec(X)] in MOI.LogDetConeTriangle(2))
@constraint(model, X .== [2 1; 1 3])
@constraint(model, u == 0.5)
optimize!(model)
@assert is_solved_and_feasible(model)
value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))(1.4979204481049089, 1.497866101098124)Other Cones and Functions
For other cones supported by JuMP, check out the MathOptInterface Manual.