运筹系列47：JuMP数学建模语言

本文参考：https://jump.dev/
首先，JuMP是一种AML（数学建模语言），而不是solver。可用的sovler列表如下：

JuMP是MOI的一个封装。

1. 快速开始

加载一个solver（这里使用GLPK），也可以先定义Model，再加载solver：

然后依次定义变量、目标函数、约束，调用optimize!(model)进行求解。

1.1 变量

使用base_name可以定义变量别名，也可以使用set_name方法，也可以使用 JuMP.variable_by_name。

has_lower_bound、lower_bound、lowerBoundRef、delete_lower_bound三者分别用于判断是否有lb，lb的值，lb的引用，删除lb。

使用下述方法定义变量数组：

如果使用表达式的话，会生成DenseAxisArrays：

还可以手动指定：

如果下标不构成矩形的话，会生成sparseAxisArrays：


一样可以手动指定：

两种定义0-1变量的方法：

julia> @variable(model, x, Bin)
julia> @variable(model, x, binary=true)

两种定义整型变量的方法

julia> @variable(model, x, Int)
julia> @variable(model, x, integer=true)

四种定义半正定变量的方法：

julia> @variable(model, x[1:2, 1:2], PSD)
julia> @variable(model, x[1:2, 1:2] in PSDCone())
julia> x = @variable(model, [1:2, 1:2], set = PSDCone())
julia> @variable(model, x[1:2, 1:2], Symmetric)

下面是定义初始值的几种方法：

最后，可以用@container宏来指定生成下标的方法：

julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(Base.OneTo(3), Base.OneTo(3)))
3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(1:3, 1:3))
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 1:3
    Dimension 2, 1:3
And data, a 3×3 Array{Int64,2}:
 2  3  4
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.vectorized_product(2:3, Base.OneTo(3)))
2-dimensional DenseAxisArray{Int64,2,...} with index sets:
    Dimension 1, 2:3
    Dimension 2, Base.OneTo(3)
And data, a 2×3 Array{Int64,2}:
 3  4  5
 4  5  6

julia> Containers.container((i, j) -> i + j, Containers.nested(() -> 1:3, i -> i:3, condition = (i, j) -> isodd(i) || isodd(j)))
SparseAxisArray{Int64,2,Tuple{Int64,Int64}} with 5 entries:
  [1, 2]  =  3
  [2, 3]  =  5
  [3, 3]  =  6
  [1, 1]  =  2
  [1, 3]  =  4

1.2 表达式

表达式可以用于目标函数和约束

model = Model()
@variable(model, x)
@variable(model, y)
ex = @expression(model, 2x + y - 1)
@objective(model, Min, 2 * ex - 1)

下面是二次的例子

二次问题还可以用标准的QP方法进行定义：

也可以使用下面的方法进行创建：

model = Model()
@variable(model, x)
@variable(model, y)
ex = AffExpr(-1.0)
add_to_expression!(ex, 2.0, x)
add_to_expression!(ex, 1.0, y)

# output

2 x + y - 1

使用drop_zeros删除0变量

非线性的表达式只能用 @NLexpression进行构建，并只能用在@NLobjective和@NLconstraint之中。

1.3 目标函数

使用@objective，下面是两种定义方式

1.4 约束

使用@constraint宏

注意，定义完约束后，生成了一个JuMP变量和model索引，两者同名。修改JuMP变量不影响索引。

如果要批量创建约束，同样有三种方法：
Arrays方法如下：

DenseAxisArray方法如下：

SparseAxisArrays方法如下：

另外，也可以使用向量化的方法：

最后，也可以对单变量使用in约束：




indicator constraints示例如下，当a为1时满足x+y<=1

使用set_normalized_rhs设置rhs，或者使用add_to_function_constant给表达式添加一定的数值。

使用set_normalized_coefficient设置系数，如下：

使用delete方法删除约束：

使用all_constraints获得约束

使用⟂定义互补约束，定义为：

布尔变量约束转换如下：

半整变量如下：

SOS1变量：最多一个不为0，如下

SOS2变量：最多两个不为0，如下

1.5 求解器

有两种指定参数的方法：

julia model = Model(Gurobi.Optimizer) 
setoptimizerattribute(model, "Presolve", 0) 
setoptimizerattribute(model, "OutputFlag", 1) 

下面是设置参数的示例：

使用write_to_file/read_from_file方法将模型读、写入文件。
调用TerminationStatusCode获得求解器的终止状态：


调用ResultStatusCode获得解的状态：

如果解的状态不是MOI.NO_SOLUTION，那么最优解可以通过value方法获得。下面是一个推荐的流程：

2. 回调函数

回调函数在混合整数规划中经常使用到，比如：分支定界法中修改braching decision，添加user-cut等等，仅限于CPLEX, GLPK, Gurobi。

2.1 Lasy constraint

下面是一个简单的例子，在x<=2不满足时，添加x<=2的约束。
回调函数传入的参数用 callback_value方法获取值
其中MOI.submit(optimizer::AbstractOptimizer,sub::AbstractSubmittable,values)在这里指将约束条件提交到model的LazyConstraint中。

model = Model(GLPK.Optimizer)
@variable(model, x <= 10, Int)
@objective(model, Max, x)

function my_callback_function(cb_data)
    x_val = callback_value(cb_data, x)
    if x_val > 2 + 1e-6
        con = @build_constraint(x <= 2)
        MOI.submit(model, MOI.LazyConstraint(cb_data), con)
    end
end

MOI.set(model, MOI.LazyConstraintCallback(), my_callback_function)

下面看一个实战案例，使用lazy constraint求解tsp问题。
首先是定义tsp问题的松弛模型和添加lazy constraint的subtour条件：

function build_tsp_model(d, n)
    model = Model(GLPK.Optimizer)
    @variable(model, x[1:n, 1:n], Bin, Symmetric)
    @objective(model, Min, sum(d .* x) / 2)
    @constraint(model, [i in 1:n], sum(x[i, :]) == 2)
    @constraint(model, [i in 1:n], x[i, i] == 0)
    return model
end

function subtour(edges::Vector{Tuple{Int,Int}}, n)
    ......
    return shortest_subtour
end

接下来我们定义subtour_elimination_callback方法，一旦违反subtour约束则进行调用：

lazy_model = build_tsp_model(d, n)
function subtour_elimination_callback(cb_data)
    cycle = subtour(callback_value.(cb_data, lazy_model[:x]))
    if !(1 < length(cycle) < n)
        return  # Only add a constraint if there is a cycle
    end
    S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
    con = @build_constraint(
        sum(lazy_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
    )
    MOI.submit(lazy_model, MOI.LazyConstraint(cb_data), con)
    return
end
MOI.set(lazy_model, MOI.LazyConstraintCallback(), subtour_elimination_callback)
optimize!(lazy_model)
objective_value(lazy_model)

2.2 User cut

User cut不改变整数解空间，它是用来剔除非整数解的。
操作方法和上面如出一辙：

model = Model(GLPK.Optimizer)
@variable(model, x <= 10.5, Int)
@objective(model, Max, x)
function my_callback_function(cb_data)
    x_val = callback_value(cb_data, x)
    con = @build_constraint(x <= floor(x_val))
    MOI.submit(model, MOI.UserCut(cb_data), con)
end
MOI.set(model, MOI.UserCutCallback(), my_callback_function)

2.3 Heuristic solutions

可以使用启发式算法添加cut，使用方法也类似。其中第三项是JuMP变量，第四项是每个变量的值。

model = Model(GLPK.Optimizer)
@variable(model, x <= 10.5, Int)
@objective(model, Max, x)
function my_callback_function(cb_data)
    x_val = callback_value(cb_data, x)
    status = MOI.submit(
        model, MOI.HeuristicSolution(cb_data), [x], [floor(Int, x_val)]
    )
    println("I submitted a heuristic solution, and the status was: ", status)
end
MOI.set(model, MOI.HeuristicCallback(), my_callback_function)