JuMP: Моделирование и решение оптимизации с Джулией

0. Зачем использовать Юлию для оптимизации?

Эта статья мояЮлия Китайское сообщество Встреча пользователей 2018 г.О совместном содержании оптимизационного моделирования и решения с Юлией.

В этом руководстве в основном представлен JuMP.jl, AML (язык алгебраического моделирования) с открытым исходным кодом на языке Julia, аналогичный AMPL, YALMIP, CVX, Poymo, GAMS и т. д. Оптимизированное моделирование + скорость передачи памяти решателя JuMP сравнивается с другими коммерческими / открытыми исходными кодами AML в таблице ниже. Convex.jl не задействован, как CVX в Matlab, еще один AML в Julia, разработанный Boyd et al.

Таблица взята из: Dunning, Iain, Joey Huchette и Miles Lubin. JuMP: язык моделирования для математической оптимизации. SIAM Review 59.2 (2017): 295-320. Дополнительные сведения о реализации JuMP и сравнения с другими AML см.

Мы обнаружили, что JuMP не уступает коммерческому AML в оптимизации скорости моделирования и значительно быстрее, чем другие AML с открытым исходным кодом.

Таблица взята с: https://www.juliaopt.org/

Еще одним большим преимуществом JuMP является бесшовная интеграция различных коммерческих/открытых решателей: Gurobi, CPLEX, SCIP, Knitro, Mosek и т. д.

1. Предварительное оптимизационное моделирование (самое основное линейное программирование)

(Эта часть в основном предназначена для читателей, у которых нет знаний в теории оптимизации. Примеры взяты из -Цинь Ханьчжан: Как понять дополнительный резерв в исследовании операций. Как следует использовать это свойство?)

Скажем, вы плотник, продающий деревянные столы и стулья ручной работы. Теперь вы хотите создать производственный план, чтобы максимизировать прибыль за месяц.

Обратите внимание, что для изготовления стола или стула требуется определенное количество дерева и времени, а в качестве небольшой ремесленной мастерской время, которое вы можете использовать для производства каждый месяц, ограничено, и вы можете ввести количество древесины также ограничено.

Для простоты предположим, что фиксированная прибыль составляет 10 долларов за стол и 3 доллара за стул. Для изготовления стола требуется 5 единиц дерева и 3 единицы времени, а для изготовления стула — 2 единицы дерева и 1 единица времени. И все столы и стулья, которые мы производим, можно продать. Предположим, у нас есть всего 200 единиц древесины и 90 единиц времени на месяц.

Очевидно, что эту задачу производственного планирования можно смоделировать с помощью линейного программирования:

$\begin{align*} \max~ & 10x_1+3x_2\tag{P}\\ \text{s.t. } & 5x_1+2x_2\leq 200 \tag{P1}\\ & 3x_1+x_2\leq 90 \tag{P2} \\ & x_1,x_2\geq 0. \end{align*}$

Согласно теории двойственности, мы также знаем, что эквивалентная двойственная задача:

$\begin{align*} \min~ & 200p_1+90p_2 \tag{D}\\ \text{s.t. } & 5p_1+3p_2\geq 10\tag{D1}\\ & 2p_1+p_2\geq 3\tag{D2} \\ & p_1,p_2\geq 0. \end{align*}$

# 利用JuMP求解原问题
function CarpenterPrimal(c,A,b)
    # 定义Model对象, OutputFlag = 0指不输出log
    Primal = Model(solver = GurobiSolver(OutputFlag = 0))
    # 定义变量，注意这里使用了宏（macro），宏的调用也是Julia&JuMP高效编译/元编程(metaprogramming)的重要技巧
    @variable(Primal, x[1:2]>=0)
    # 定义不等式约束
    constr = @constraint(Primal, A*x.<=b)
    # 定义目标函数
    @objective(Primal, Max, dot(c,x))
    # 求解
    solve(Primal)
    # 返回最优目标函数值，最优解（原问题），最优解（对偶问题）
    return getobjectivevalue(Primal), getvalue(x), getdual(constr)
end

При вызове функции CarpenterPrimal оптимальным решением будет изготовление 30 столов и ни одного стула: x_1^*=30,x_2^*=0. . Давайте посмотрим на другую информацию, стоящую за этим решением. В то же время теневые цены, соответствующие древесине и времени p_1^*,p_2^* равны 0 и 10/3.

Отметим, что анализ чувствительности также можно выполнить с помощью JuMP. Рассмотрим следующие четыре варианта (P) задач: (Конечно, все они могут быть решены вручную и проверены с помощью CarpenterPrimal)

(1) Измените правый конечный элемент (P1) в (P) на 250, то есть общее количество древесины увеличится. Остальные остаются без изменений. Какое оптимальное решение?
(2) Изменить правый конечный элемент (P2) в (P) на 120, то есть увеличить время производства. Остальные остаются без изменений. Какое оптимальное решение? Связь между увеличенной прибылью (по сравнению с первоначальной (P)) и теневой ценой (P2)?
(3) Если предположить, что удельная прибыль от стола постоянна, насколько увеличится удельная прибыль от стула, когда оптимальное решение (P) начнет производить стулья (без производства столов)?
(4) Предполагая, что общее количество времени по-прежнему равно 90, насколько уменьшится общее количество древесины, когда скрытая цена древесины строго больше 0 (вместо этого скрытая цена времени становится равной 0)?

Если вы не знакомы с этими основными понятиями, обратитесь к моему ответу, упомянутому в начале этого раздела.

2. Пример генерации столбца в линейном программировании: крупномасштабная аддитивная выпуклая регрессия

данные данные $Y\in \mathbb{R}^n, X\in \mathbb{R}^{n\times d},$ Рассмотрим следующую задачу линейного программирования:

$\begin{align*} \min~ & \sum_{i=1}^n \left|Y_i-\sum_{j=1}^d f_j(X_{ij})\right|\\ \text{s.t. } & f_j:\mathbb{R}\rightarrow\mathbb{R} \text{ is convex}, \forall~j=1,\ldots, d. \end{align*}$

Обратите внимание, что если мы добавим дополнительный регулярный член к целевой функции $\lambda\sum_{j=1}^d \left| \sum_{i=1}^n f_j(X_{ij}) \right|$ ( $\lambda>0$ ), Тогда мы можем получить разреженную аддитивную модель. См.: Равикумар, Прадип и др. «Разреженные аддитивные модели», Журнал Королевского статистического общества: Серия B (Статистическая методология) 71.5 (2009): 1009-1030.

Приведенная выше модель является линейной программой, поскольку она эквивалентна (обозначается [:]_j за $X_{:j}$ Последовательность, соответствующая исходному индексу после сортировки по возрастанию)

$\begin{align*} \min~ & \sum_{i=1}^n r_i^++r_i^-\\ \text{s.t. } & Y_i - \sum_{j=1}^d f_{ij} = r_i^+-r_i^-, \forall~i=1,\ldots,n\\ & (f_{[i]_j j}-f_{[i-1]_j j})(X_{[i+1]_j j}-X_{[i]_j j})\leq (f_{[i+1]_j j}-f_{[i]_j j})(X_{[i]_j j}-X_{[i-1]_j j}),\forall~i=2,\ldots,n-1,\forall~j=1,\ldots,d. \end{align*}$

# 生成样本
n = 5000
d = 2
function_list = [x->0.5*abs.(x).^2,x->2*abs.(x),x->5*x.^2,x->0.5*exp.(x),x->0.3*x,x->0.7*x,x->2*x.^2,x->0.1*exp.(x),x->x.^2,x->exp.(x)]
srand(123)
X = rand(n,d) - 0.5
Y = 0.02*randn(n,1)
for i = 1:d
    Y = Y + function_list[mod(i-1,10)+1](X[:,i])
end

# Full-sized的线性规划模型
function Full_ConvReg(Max_T)
    M_conv_full = Model(solver=GurobiSolver(TimeLimit=Max_T, OutputFlag=0))
    @variable(M_conv_full, f_hat[1:n,1:d])
    @variable(M_conv_full, Res_pos[1:n]>=0)
    @variable(M_conv_full, Res_neg[1:n]>=0)    
    for k = 1:d
        xx = X[:,k] 
        ids = sortperm(vec(xx),rev=false)
        xx = xx[ids]
        for i in 2:(n-1)
            d1 = (xx[i+1]-xx[i]); d2=(xx[i]-xx[i-1]);
            @constraint(M_conv_full, (f_hat[ids[i],k]-f_hat[ids[i-1],k])*d1<=(f_hat[ids[i+1],k]-f_hat[ids[i],k])*d2);
        end   
    end
    
        for i in 1:n
            @constraint(M_conv_full, Y[i] - sum(f_hat[i,k] for k=1:d) == Res_pos[i] - Res_neg[i]  )
        end  
        @objective(M_conv_full, :Min, sum(Res_pos[i] + Res_neg[i] for i = 1:n) )   
        solve(M_conv_full)   
    return getvalue(f_hat)
end

Далее мы расскажем, что такое генерация столбцов, и покажем, как использовать JuMP и Gurobi для выполнения генерации столбцов для решения вышеуказанных проблем. Этот метод основан на симплексном методе. То есть мы по умолчанию используем все переменные равные 0, а не базовые переменные в начале, и на каждом шаге итерации мы используем значение двойной переменной предыдущего раунда для расчета приведенной стоимости (нужно только вычислить умножение матрица и вектор за один шаг), и положить Переменная с наибольшей отрицательной приведенной стоимостью добавляется (становится базовой переменной). Также обратите внимание, что, поскольку матрица ограничений исходной задачи фактически имеет блочную характеристику (согласно j=1,…,d По сути, мы можем взять формулировку Данцига-Вульфа, то есть использовать метод циклического обновления для расчета приведенной стоимости отдельно для каждого блока в каждой итерации (быстрее, чем вычислять их вместе).раз). Когда все приведенные затраты неотрицательны, оптимальное решение находится: алгоритм завершает работу.

Для учащихся, не знакомых с Данцигом-Вульфом или генерацией столбцов, см. некоторые классические материалы, такие как:Пирс О. Калифорнийский университет в Лувене. Не голоден / Энтони. Боюсь…

Конечно, принцип генерации ограничений в пятом примере аналогичен принципу генерации столбца (переменной), и этот пример может быть легче читать учащимся с нулевыми базовыми знаниями.

# Column Generation化的线性规划模型
function CG_ConvReg(Max_T)
    # build the LO: no Δ included
    M_conv = Model(solver=GurobiSolver(TimeLimit = Max_T,OutputFlag = 0,Method = 0))
    @variable(M_conv, Res_pos[1:n]>=0)
    @variable(M_conv, Res_neg[1:n]>=0)
    @constraintref constr[1:n]
    
    MaxIter = 1000
    
    XX = zeros(n,d)
    L = zeros(n,n+2,d)
    ids = convert(Array{Int64,2},zeros(n,d))
    ids_rec = convert(Array{Int64,2},zeros(n,d))
    
    for k = 1:d
        # sort the sequences
        ids[:,k] = sortperm(vec(X[:,k]),rev=false)
        ids_rec[:,k] = sortperm(ids[:,k],rev=false)
        XX[:,k] = X[ids[:,k],k]  
        # needed for column generation
        L[:,1,k] = ones(n,1)
        L[:,2,k] = -ones(n,1)
        L[2,3,k] = XX[2] - XX[1]
        L[2,4,k] = - XX[2] + XX[1]
        for i = 3:n
            L[i,3,k] = XX[i,k] - XX[1,k]
            L[i,4,k] = - XX[i,k] + XX[1,k]
            L[i,5:i+2,k] = [XX[i,k]-XX[j+1,k] for j=1:i-2]'
        end
        L[:,:,k] = L[ids_rec[:,k],:,k] 
    end
    
    @constraint(M_conv, constr[i=1:n], Res_pos[i] - Res_neg[i] == Y[i])
    @objective(M_conv, :Min, sum(Res_pos[i] + Res_neg[i] for i = 1:n) )
    
    solve(M_conv)
    dual_var = getdual(constr)
                    
    NewColumns = [Variable[] for i=1:d]
    IdxList = [[] for i=1:d]                       
    
    # Implement column generation                 
    iter = 1   
    flag = vec(ones(d,1))
    while sum(flag) > 1e-3 && iter < MaxIter
    #for t = 1:30
    for k = 1:d    
        if flag[k] == 1 
            Δ_r = -L[:,:,k]'*dual_var 
        else
            continue
        end
        if findmin(Δ_r)[1] < -1e-3
            idx = findmin(Δ_r)[2]
            constr_list = constr[1:n]
            coeff_list = L[:,idx,k]
            @variable(M_conv, Δ_temp >= 0, objective = 0., inconstraints = constr_list, coefficients =  coeff_list ) 
            #print("CG iteration ", iter,": enter Δ ",idx," for covariate ",k, "\n")
            push!(NewColumns[k],Δ_temp)
            push!(IdxList[k],idx)
            solve(M_conv) 
            iter = iter + 1
            dual_var = getdual(constr);
        else
            #print("Skip covariate ", k,"\n")
            flag[k] = 0
        end
    end
    end
    #print("CG ends after ", iter-1, " itertaions.\n")
                
    # Retrieve the ϕ values
    delta = zeros(n+2,d)
    zz = zeros(n,d)
    phi = zeros(n,d)
    for k = 1:d
        delta[IdxList[k],k]=getvalue(NewColumns[k])
        if isempty(find(IdxList[k].==2))
                zz[1,k] = delta[1,k]
        else
                zz[1,k] = -delta[2,k]
        end
        if isempty(find(IdxList[k].==4))
                zz[2,k] = delta[3,k]
        else
                zz[2,k] = -delta[4,k]
        end
        phi[1,k] = zz[1,k]
        phi[2,k] = zz[1,k]+zz[2,k]*(XX[2,k]-XX[1,k])                           
        for i = 3:n
            zz[i,k] = zz[i-1,k] + delta[i+2,k]
            phi[i,k] = zz[i,k]*(XX[i,k]-XX[i-1,k])+phi[i-1,k]
        end
    end
    
    return phi
end

phi = CG_ConvReg(200);
# 画出第一维上的解
XX = zeros(n,d)
ids = convert(Array{Int64,2},zeros(n,d))
ids_rec = convert(Array{Int64,2},zeros(n,d))
    
for k = 1:d 
        # sort the sequences
        ids[:,k] = sortperm(vec(X[:,k]),rev=false)
        ids_rec[:,k] = sortperm(ids[:,k],rev=false)
        XX[:,k] = X[ids[:,k],k] 
end
plot(XX[:,1],phi[:,1])

Далее мы сравним скорость решения этих двух формулировок.

@benchmark CG_ConvReg(200)

@benchmark Full_ConvReg(600)

Удивительно, но алгоритм, основанный на древнем симплекс-методе и генерации столбцов, намного быстрее, чем недетерминированный параллельный метод Gurobi по умолчанию для решения LP! Это доказывает, что иногда более разумное моделирование дает более быстрые алгоритмы для конкретных задач оптимизации.

3. Надежное линейное программирование

Мы рассматриваем задачу потока минимальной стоимости, где стоимость на каждом ребре имеет неопределенность. То есть мы рассматриваем ориентированный граф G=(V,E)，V={1,…,n} , рассмотрим задачу оптимизации:

$\begin{align*} \min~ & \sum_{(i,j)\in E} c_{ij}x_{ij}\\ \text{s.t. } & \sum_{(1,i)\in E}x_{1i} = 1\\ & \sum_{(i,n)\in E}x_{in} = 1\\ & \sum_{(i,j)\in E}x_{ij} = \sum_{(j,k)\in E} x_{jk},\forall~j\in V \backslash \{1,n\}\\ & 0\leq x_{ij}\leq C_{ij}. \end{align*}$

# Nominal Problem
function NominalProblem(n,μ,δ)
    NetworkModel = Model(solver=GurobiSolver(OutputFlag=0))
    capacity = (ones(n,n)-eye(n,n))*0.5
    capacity[1,n] = 0
    capacity[n,1] = 0
    @variable(NetworkModel, flow[i=1:n,j=1:n]>=0)
    @constraint(NetworkModel, flow .<= capacity)
    @constraint(NetworkModel, sum(flow[1,i] for i=1:n)==1)
    @constraint(NetworkModel, sum(flow[i,n] for i=1:n)==1)
    for j = 2:n-1
        @constraint(NetworkModel, sum(flow[i,j] for i=1:n) == sum(flow[j,k] for k=1:n)  )
    end
    @objective(NetworkModel,Min, sum(flow[i,j]*μ[i,j] for i=1:n for j=1:n) );
    solve(NetworkModel)
    return getvalue(flow)
end

Мы рассматриваем следующие три набора неопределенностей:

$\begin{align*} & (i)~\text{Box Uncertainty}: \mathcal{U}_{L_\infty} = \{\textbf{c}:μ_{ij}-\delta_{ij}\gamma_{ij}\leq c_{ij}\leq \mu_{ij}+\delta_{ij}\gamma_{ij}, \|\gamma\|_{\infty}\leq \Gamma \}\\ & (ii)~\text{Polyhedral Uncertainty}: \mathcal{U}_{L_1} = \{\textbf{c}:μ_{ij}-\delta_{ij}\gamma_{ij}\leq c_{ij}\leq \mu_{ij}+\delta_{ij}\gamma_{ij}, \|\gamma\|_1\leq \Gamma \}\\ & (iii)~\text{Elliopsoidal Uncertainty}: \mathcal{U}_{L_2} = \{\textbf{c}:μ_{ij}-\delta_{ij}\gamma_{ij}\leq c_{ij}\leq \mu_{ij}+\delta_{ij}\gamma_{ij}, \|\gamma\|_2\leq \Gamma \}\\ \end{align*}$

На самом деле, легко понять, что робастная задача о потоке минимальных затрат (i)(ii) по-прежнему может быть записана в виде линейной программы, а робастная задача о потоке минимальных затрат в (iii) может быть записана как задача второго порядка. задача оптимизации конуса. Однако здесь мы показываем, что использование пакета JuMPeR может сэкономить этап получения надежного аналога.

# Robust Problem
function RobustProblem(n,μ,δ,Γ,norm_type)
    NetworkModel_robust = RobustModel(solver=GurobiSolver(OutputFlag=0))
    capacity = (ones(n,n)-eye(n,n))*0.5
    capacity[1,n] = 0
    capacity[n,1] = 0
    @variable(NetworkModel_robust, flow[i=1:n,j=1:n]>=0)
    @uncertain(NetworkModel_robust, cost[i=1:n,j=1:n])
    @uncertain(NetworkModel_robust, r[i=1:n,j=1:n])
    @variable(NetworkModel_robust, obj) 
    for i = 1:n
        for j = 1:n
        @constraint(NetworkModel_robust, cost[i,j] == μ[i,j]+r[i,j]*δ[i,j])
        end
    end  
    @constraint(NetworkModel_robust, norm(r,norm_type) <= Γ)
    @constraint(NetworkModel_robust, flow .<= capacity)
    @constraint(NetworkModel_robust, sum(flow[1,i] for i=1:n) == 1)
    @constraint(NetworkModel_robust, sum(flow[i,n] for i=1:n) == 1)
    for j = 2:n-1
        @constraint(NetworkModel_robust, sum(flow[i,j] for i=1:n) == sum(flow[j,k] for k=1:n))
    end
    @constraint(NetworkModel_robust, sum(cost[i,j]*flow[i,j] for i=1:n for j=1:n) <= obj )
    @objective(NetworkModel_robust, Min, obj)
    solve(NetworkModel_robust)
    return getvalue(flow)
end

Давайте проведем здесь несколько симуляций, чтобы увидеть производительность надежного решения. Предполагая, что G полностью ассоциативна, а пропускная способность всех ребер равна 0,5, мы случайным образом генерируем $μ_{ij}∼Uniform(0,10)$ $δ_{ij}∼Uniform(0,μ_{ij})$ , а затем случайным образом генерируем $c_{ij}$ Посмотрим на реальную стоимость. Мы берем решение номинальной задачи в качестве эталона и выводим разницу между стоимостью робастного решения и разностью.

function Evaluation()
    n_type = [10,25,50,75,100]
    Γ_type = [1e-4,1e-1,1e1,1e4]
    report_data = zeros(13,5)
    for n_idx = 1:5
        n = n_type[n_idx]
        count = 2
        μ = 10*rand(n,n)
        δ = zeros(n,n)
        for i = 1:n
            for j = 1:n
                δ[i,j] = μ[i,j]*rand()
            end
        end
        cost_sim = zeros(n,n)
                for i = 1:n
                    for j = 1:n
                        cost_sim[i,j] = μ[i,j] + δ[i,j] * (rand()-0.5) * 2
                    end
                end
        sol_0 = NominalProblem(n,μ,δ)
        report_data[1,n_idx] = sum(sum(cost_sim.*sol_0))
        for Γ in Γ_type
            print("n=",n," Γ=",Γ)
            sol_inf = RobustProblem(n,μ,δ,Γ,Inf)
            sol_1 = RobustProblem(n,μ,δ,Γ,1)
            sol_2 = RobustProblem(n,μ,δ,Γ,2)    
            report_data[count,n_idx] = sum(sum(cost_sim.*(sol_inf-sol_0)))
            report_data[count+1,n_idx] = sum(sum(cost_sim.*(sol_1-sol_0) ))
            report_data[count+2,n_idx] = sum(sum(cost_sim.*(sol_2-sol_0) ))
            count = count + 3
        end
    end
    return report_data
end
output = Evaluation();
# 输出Γ=10时三种uncertainty set给出的difference
temp_p = plot([10,25,50,75,100],output[8,:],xlabel = "n", label = "Uinf")
plot!(temp_p,[10,25,50,75,100],output[9,:], label = "U1")
plot!(temp_p,[10,25,50,75,100],output[10,:], label = "U2")

Тогда мы видим, что коробочное ограничение дает более надежное решение, а задача робастной оптимизации при построении неопределенности L1/2 является более мягкой.

JuMPeR также может моделировать многоэтапные задачи надежной оптимизации, особенно для политик с аффинной настройкой (которые намного свободнее, чем формулировка набора неопределенностей). Здесь он только упоминается, а не детализируется, желающие могут изучить его самостоятельно, что на самом деле очень полезно.

4. Приблизительное решение Стабильное число: положительное полуопределенное программирование/целочисленное программирование

неориентированный граф G=(V,E) Стабильное (независимое) множество V — это подмножество V, в котором все узлы не связаны друг с другом. Максимальная мощность стабильного множества определяется как α(G), также известный как число стабильности графа.Предполагая |V|=n.Естественно, его целочисленное программирование 0/1 можно записать как:

$\begin{align} \alpha(G) = & \max_x~ \sum_{i=1}^n x_i\\ &\text{s.t. } x_i+x_j\leq 1,~ \text{if } (i,j)\in E,\\ & x_i\in \{0,1\},\forall~i=1,\ldots,n. \end{align}$

использовать x_i+x_j≤1 и x_ix_j=0 Эквивалентно, мы можем получить релаксацию положительной полуопределенной программы:

$\begin{align} \alpha(G) = & \max_X~ \text{tr}(ee^TX) \\ &\text{s.t. } X_{ij} = 0,~ \text{if } (i,j)\in E,\\ & \text{tr}(X) = 1,\\ & X\succeq 0,\\ & X\geq 0. \end{align}$

Пример источника: Пена, Хавьер, Хуан Вера и Луис Ф. Сулуага, «Вычисление числа устойчивости графа с помощью линейного и полуопределенного программирования», SIAM Journal on Optimization 18.1 (2007): 87-105.

function StableBin(A)
    n = size(A,1)
    StableB = Model(solver = GurobiSolver(OutputFlag = 0))
    @variable(StableB, x[1:n], Bin)
    for i = 1:n
        for j in find(A[i,:].== 1)
            @constraint(StableB, x[i] + x[j] <= 1)
        end
    end
    @objective(StableB, Max, sum(x))
    solve(StableB)
    return getobjectivevalue(StableB)
end

function StableSDP(A)
    n = size(A,1)
    StableDD = Model(solver = MosekSolver(MSK_IPAR_LOG = 0))
    @variable(StableDD, X[1:n,1:n])
    @SDconstraint(StableDD, X>=0)
    @constraint(StableDD, X.>=0)
    @constraint(StableDD, sum(sum(A.*X)) == 0 )
    @constraint(StableDD, sum(X[i,i] for i=1:n) == 1)
    @objective(StableDD, Max, sum(sum(X)))
    solve(StableDD)
    return getobjectivevalue(StableDD)
end

Возьмите это G в качестве примера.

A14 = [0 1 1 0 0 0 0 0 0 0 0 0 0 0; 
       1 0 0 1 0 0 0 0 0 0 0 0 0 0;
       1 0 0 0 1 1 0 0 1 0 0 1 0 0;
       0 1 0 0 1 1 0 0 1 0 0 1 0 0;
       0 0 1 1 0 0 1 0 0 0 1 0 0 1;
       0 0 1 1 0 0 0 1 0 0 1 0 0 1;
       0 0 0 0 0 1 0 1 0 0 0 0 0 0;
       0 0 0 0 0 1 1 0 0 0 0 0 0 0;
       0 0 1 1 0 0 0 0 0 1 0 0 0 0;
       0 0 0 0 0 0 0 0 1 0 1 0 0 0;
       0 0 0 0 1 1 0 0 0 1 0 1 0 0;
       0 0 1 1 0 0 0 0 0 0 1 0 1 0;
       0 0 0 0 0 0 0 0 0 0 0 1 0 1;
       0 0 0 0 1 1 0 0 0 0 0 0 1 0];
# α(G14) = 5
StableBin(A14)
# 松弛之后的解：5.694362954824882
StableSDP(A14)

5. Генерация ленивых констант: аппроксимация сфер L2 с кусочно-линейными функциями в целочисленном программировании

Мы указали, что Gurobi еще не поддерживает комбинацию целочисленного программирования и генерации столбцов (т. е. ветвления и цены), но пользовательский разрез поддерживает это (т. е. ветвления и разрезания)! Короче говоря, мы можем добавить только небольшое количество ограничений в начале, и будут добавлены нарушенные ограничения в процессе решения целочисленного программирования, поэтому это также называется ленивым обратным вызовом.

Этот пример исходит из:GitHub.com/Julia opt/ Согласно…

Задача, которую мы хотим решить, — это задача программирования конусов второго порядка с целочисленными ограничениями:

$\begin{align*} \max ~ & c^T x\\ \text{s.t. } & \|x\|_2\leq \Gamma,\\ & x\in \mathbb{Z}^n. \end{align*}$

function solve_ball(c, Γ, ϵ=1e-6)
    n = length(c)
    m = Model(solver=GurobiSolver(OutputFlag=0))
    # 初始consrtaint: 一个box
    @variable(m, -Γ ≤ x[1:n] ≤ Γ, Int)
    # 定义目标函数
    @objective(m, Max, dot(c,x))
    #核心：定义callback function，记录加入cut的数量
    num_callbacks = 0
    function norm_callback(cb)
        num_callbacks += 1
        N = getvalue(x)
        # 求得当前x的L2 norm
        L = norm(N)
        # 如果足够小，说明已经得到一个可行解，即解最优
        if L ≤ Γ + ϵ
            return
        end
        # 不然的话，加入cut!注意这个cut将使得x变得不可行（infeasible），下步迭代必然会得到一个新的解
        @lazyconstraint(cb, dot(N,x) ≤ Γ*L)
    end
    # 将callback函数加入JuMP/Gurobi模型
    addlazycallback(m, norm_callback)
    
    #求解
    solve(m)
    return getvalue(x), num_callbacks
end

#产生一个随机样例
srand(1234)
n = 2
c = rand(n)
Γ = 50.0

# 求解，输出callback次数
sol, num_callbacks = solve_ball(c, Γ)
println(sol)
println(norm(sol))
println(num_callbacks)

После 11 ленивых обратных вызовов мы получаем оптимальное решение. Далее еще один пример взаимодействия, мне лень делать здесь анимацию, просто дайте код. Пожалуйста, испытайте эффект на себе в блокноте Jupyter. Фактически это означает, что вы можете вручную настроить параметры на панели, чтобы увидеть оптимальное решение в разных ситуациях и изображение допустимой области на двухмерной плоскости.

# 这里展示Interact和Compose包，可用来进行交互
# 画出2D情况下的解
set_default_graphic_size(8cm, 8cm)

@manipulate for c1 in -1:0.1:+1, c2 in -1:0.1:+1, logϵ in -4:2
    sol, _ = solve_ball([c1,c2], 100, 10.0^logϵ)

    compose(context(),
    compose(context(),
            line([(0.5,0.5),(0.5+sol[1]/300,0.5+sol[2]/300)]),
            Compose.stroke("black")),
    compose(context(),
            circle((0.5 + (100/norm(sol))*sol/300)...,0.02),
            fill("red")),
    compose(context(),circle(0.5,0.5,0.333),fill("lightblue"))
    )
end

6. Общее нелинейное целочисленное программирование: задача о плотно упакованной сфере

Давайте рассмотрим такую задачу: есть k шаров радиуса r и прямоугольная коробка размером d1 × d2, можно ли как-то вместить эти шары в прямоугольник? Очевидно, более общая задача состоит в том, чтобы найти наибольшее k для любого прямоугольника. Однако эту общую проблему можно решить, последовательно решая (увеличивая k) предыдущую проблему решения, поэтому мы обсуждаем проблему решения для данного k.

определить переменную $(p_{11},…,p_{1k}),(p_{21},…,p_{2k} )$ — координаты центра окружности по х и у. Тогда проблема принятия решения может быть решена с помощью следующей задачи невыпуклого нелинейного программирования, см. также: Биргин, Эрнесто Г., Дж. М. Мартинес и Дебора П. Ронкони. Оптимизация упаковка баллонов в прямоугольный контейнер: нелинейный подход». Европейский журнал операционных исследований 160.1 (2005): 19-33.

$\begin{align} \min~ & \sum_{i\neq j\in \{1,\ldots, k\} } \max\left(0,2r - \sqrt{(p_i^1-p_j^1)^2+(p_i^2-p_j^2)^2} \right)\\ \text{s.t. } & r\leq p_1^i\leq d_1-r,\forall~ i=1,\ldots,k, \\ & r\leq p_2^i\leq d_2-r,\forall~ i=1,\ldots,k. \\ \end{align}$

Подобно идее из предыдущего раздела, мы также используем здесь кусочно-линейные функции и целочисленное программирование для аппроксимации функции L2. Конечно, здесь мы используем более общий подход и используем пакет PiecewiseLinear для облегчения завершения. В частности, мы рассматриваем целочисленное линейное программирование для аппроксимации невыпуклых ограничений в нелинейном программировании следующим образом.

$\begin{align} \min~ & \sum_{i\neq j\in \{1,\ldots, k\} } z_{ij}\\ \text{s.t. } & z_{ij}\geq 2r - \sqrt{(p_i^1-p_j^1)^2+(p_i^2-p_j^2)^2}, \forall ~ i\neq j\in \{1,\ldots,k\},\\ &r\leq p_1^i\leq d_1-r,\forall~ i=1,\ldots,k, \\ & r\leq p_2^i\leq d_2-r,\forall~ i=1,\ldots,k. \\ & z_{ij}\geq 0, \forall ~ i\neq j\in \{1,\ldots,k\}. \end{align}$

# 画图
function PlotCirclesRectangle(p1,p2)
    p = scatter((p2),(p1))
    for i = 1:k
        plot!(p,(p2[i]).+r*cos.(linspace(0,2*π,100)), (p1[i]).+r*sin.(linspace(0,2*π,100)), color = "blue", legend=false, aspect_ratio = 1 )
    end
    plot!(p,linspace(0,d2,100),zeros(100,1),color = "red")
    plot!(p,zeros(100,1),linspace(0,d1,100),color = "red")
    plot!(p,linspace(0,d2,100),ones(100,1)*d1,color = "red")
    plot!(p,d2*ones(100,1),linspace(0,d1,100),color = "red")
    p
end

# 考虑这样一组问题
r = 25
d1 = 160
d2 = 190
k = 11;
# MIO
function PackCirclesGen(r,d1,d2,k,Method,MaxTime)
    tic()
    PackCircles = Model(solver = GurobiSolver(MIPGap = 1e-3, OutputFlag = 0, TimeLimit = MaxTime))
    @variable(PackCircles, p1[1:k]>=r)
    @variable(PackCircles, p2[1:k]>=r)
    @variable(PackCircles, s1[1:k,1:k])
    @variable(PackCircles, s2[1:k,1:k]>=0)
    @variable(PackCircles, z[1:k,1:k]>=0)

    @constraint(PackCircles, p1.<=d1-r)
    @constraint(PackCircles, p2.<=d2-r)

    for i = 1:k
        for j = i+1:k
        if j - i > 2 
            continue
        end
        @constraint(PackCircles,s1[i,j] == p1[j]-p1[i])
        @constraint(PackCircles,s2[i,j] == p2[j]-p2[i])
        #使用PiecewiseLinearOpt包近似nonconvex constraint
        fun_dist = piecewiselinear(PackCircles, s1[i,j] , s2[i,j], -(j-i)*r*2:r/5:(j-i)*r*2, 0:r/5:(j-i)*r*2, (u,v) -> (u^2+v^2)^0.5, method=Method)
        @constraint(PackCircles, z[i,j] >=  2*r - fun_dist )
        end
    end
    @objective(PackCircles, Min, sum(z[i,j] for i = 1:k for j = i+1:k))
    solve(PackCircles)
    return getvalue(p1),getvalue(p2),toc()
end
p1,p2,time = PackCirclesGen(r,d1,d2,k,:ZigZagInteger,Inf);
# 一个可行解
PlotCirclesRectangle(p1,p2)

Эту программу также можно использовать для исследования этого вопроса Чжиху:

Сколько кругов диаметром 1 может поместиться в прямоугольнике 2×1000?

# 不同bivariate分段线性方法的比较
r = 21
d1 = 120
d2 = 120
k = 6;
Method_List = [:CC,:MC,:DisaggLogarithmic,:SOS2,:Logarithmic,:ZigZag,:ZigZagInteger]
# The experiment, note the cutoff time set as 600 seconds
for Method in Method_List
    p1,p2,time = PackCirclesGen(r,d1,d2,k,Method,600)
    print("The method ", Method, " uses ", time, " seconds!\n")
end

Мы заметили, что скорость расчета различных методов линейной кусочной аппроксимации сильно различается, ZigZag — самый быстрый, а такие формулировки, как CC/MC, вообще не рекомендуются. Подробный теоретический подход см.: Huchette, Joey, and Juan Pablo Vielma, «Невыпуклые кусочно-линейные функции: расширенные формулировки и простые инструменты моделирования».

7. Еще

Ссылка на этот учебник на GitHub, включая другие материалы о встрече, связанные с Джулией:

JuliaCN/MeetUpMaterials

Другие оптимизированные пакеты, поддерживаемые JuMP, можно найти здесь:www.juliaopt.org/packages/

Например, MultiJuMP.jl, который поддерживает многокритериальную оптимизацию, JuMPChance.jl, который поддерживает случайное ограничение, PolyJuMP.jl, который поддерживает полиномиальную оптимизацию и поддерживает стохастическое динамическое программирование, нелинейное управление, задачи ограничения равновесия и алгоритмы теории графов, пакет для метаэвристика и др.

Я также написал текст Джулии Амвей раньше:

Цинь Ханьчжан: Julia: эффективный и простой в использовании язык программирования для численных вычислений/оптимизации

Практикуйтесь больше, пробуйте больше... Практика делает совершенным!

Любые вопросы или общение, добро пожаловать в контактя!