↑↑↑↑↑ нажмите выше синее слово Подписывайтесь на нас!
Оригинал "Стратегия ИЛИ стратегия"
Автор: Цинь Ханьчжан.
0 Зачем использовать Юлию для оптимизации?
Это статья, которую я поделился на собрании пользователей Julia Chinese Community 2018 об оптимизации моделирования и решения с помощью Julia.
В этом руководстве в основном представлен JuMP.jl, AML (язык алгебраического моделирования) с открытым исходным кодом на языке Julia, аналогичный AMPL, YALMIP, CVX, Poymo, GAMS и т. д. Оптимизированное моделирование + скорость передачи памяти решателя JuMP сравнивается с другими коммерческими / открытыми исходными кодами AML в таблице ниже. Convex.jl не задействован, как CVX в Matlab, еще один AML в Julia, разработанный Boyd et al.
Мы обнаружили, что JuMP не уступает коммерческому AML в оптимизации скорости моделирования и значительно быстрее, чем другие AML с открытым исходным кодом.
Еще одним большим преимуществом JuMP является бесшовная интеграция различных коммерческих/открытых решателей: Gurobi, CPLEX, SCIP, Knitro, Mosek и т. д.
1 Предварительное оптимизационное моделирование (самое основное линейное программирование)(Эта часть в основном предназначена для читателей, у которых нет базовых знаний в теории оптимизации. Пример взят из книги Цинь Ханьчжан: как понять дополнительный резерв в исследовании операций. Как следует использовать это свойство?)
Скажем, вы плотник, продающий деревянные столы и стулья ручной работы. Теперь вы хотите создать производственный план, чтобы максимизировать прибыль за месяц.
Обратите внимание, что для изготовления стола или стула требуется определенное количество дерева и времени, а в качестве небольшой ремесленной мастерской время, которое вы можете использовать для производства каждый месяц, ограничено, и вы можете ввести количество древесины также ограничено.
Для простоты предположим, что фиксированная прибыль составляет 10 долларов за стол и 3 доллара за стул. Для изготовления стола требуется 5 единиц дерева и 3 единицы времени, а для изготовления стула — 2 единицы дерева и 1 единица времени. И все столы и стулья, которые мы производим, можно продать. Предположим, у нас есть всего 200 единиц древесины и 90 единиц времени на месяц.
Очевидно, что эту задачу производственного планирования можно смоделировать с помощью линейного программирования:
Согласно двойственной теории мы также знаем, что эквивалентная двойственная задача:
# 利用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 столов и ни одного стула:. Давайте посмотрим на другую информацию, стоящую за этим решением. В то же время теневые цены, соответствующие древесине и времени
равны 0 и 10/3.
Отметим, что анализ чувствительности также можно выполнить с помощью JuMP. Рассмотрим следующие четыре варианта (P) задач: (Конечно, все они могут быть решены вручную и проверены с помощью CarpenterPrimal)
(1) Измените правый конечный элемент (P1) в (P) на 250, то есть общее количество древесины увеличится. Остальные остаются без изменений. Какое оптимальное решение? (2) Изменить правый конечный элемент (P2) в (P) на 120, то есть увеличить время производства. Остальные остаются без изменений. Какое оптимальное решение? Связь между увеличенной прибылью (по сравнению с первоначальной (P)) и теневой ценой (P2)? (3) Если предположить, что удельная прибыль от стола постоянна, насколько увеличится удельная прибыль от стула, когда оптимальное решение (P) начнет производить стулья (без производства столов)? (4) Предполагая, что общее количество времени по-прежнему равно 90, насколько уменьшится общее количество древесины, когда скрытая цена на древесину строго больше 0 (вместо этого скрытая цена времени становится равной 0)?
Если вы не знакомы с этими основными понятиями, обратитесь к моему ответу, упомянутому в начале этого раздела.
2 Экземпляр генерации столбцов в линейном программировании: Крупномасштабная аддитивная выпуклая регрессия
данные данныеРассмотрим следующую задачу линейного программирования:
Обратите внимание, что если мы добавим дополнительный регулярный член к целевой функции,
Тогда мы можем получить разреженную аддитивную модель. См.: Равикумар, Прадип и др. «Разреженные аддитивные модели», Журнал Королевского статистического общества: Серия B (Статистическая методология) 71.5 (2009): 1009-1030.
Вышеупомянутая модель является линейным программированием, потому что она эквивалентна (помните за
Последовательность, соответствующая исходному индексу после сортировки по возрастанию)
# 生成样本n = 5000d = 2function_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.5Y = 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 и не являются базовыми переменными в начале, и на каждом шаге итерации мы используем значение двойной переменной предыдущего раунда для расчета приведенной стоимости (нужно только вычислить умножение матрица и вектор за один шаг), и положить Переменная с наибольшей отрицательной приведенной стоимостью добавляется (становится базовой переменной). Также обратите внимание, что, поскольку матрица ограничений исходной задачи фактически имеет блочную характеристику (согласноНа самом деле можно взять формулировку Данцига-Вульфа, то есть использовать метод циклического обновления для расчета приведенного отдельно для каждого блока в каждой итерации
стоимость (быстрее, чем считать вместе)
раз). когда все уменьшилось
Когда стоимость неотрицательна, оптимальное решение найдено: алгоритм завершает работу.
Для студентов, которые не знакомы с Dantzig-Wolfe или генерацией столбцов, пожалуйста, обратитесь к некоторым классическим материалам, таким как: https://perso.uclouvain.be/anth ony.papavasiliou/public_html/DantzigWolfe.pdf
Конечно, принцип генерации ограничений в пятом примере аналогичен принципу генерации столбца (переменной), и этот пример может быть легче читать учащимся с нулевыми базовыми знаниями.
# 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 phiendскопировать код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] endplot(XX[:,1],phi[:,1])скопировать кодДалее мы сравним скорость решения этих двух формулировок.
@benchmark CG_ConvReg(200)скопировать код@benchmark Full_ConvReg(600)скопировать кодУдивительно, но алгоритм, основанный на древнем симплекс-методе и генерации столбцов, намного быстрее, чем недетерминированный параллельный метод Gurobi по умолчанию для решения LP! Это доказывает, что иногда более разумное моделирование дает более быстрые алгоритмы для конкретных задач оптимизации.
3 Надежное линейное программирование (Надежное линейное программирование)
Мы рассматриваем задачу потока минимальной стоимости, где стоимость на каждом ребре имеет неопределенность. То есть мы рассматриваем направленный, рассмотрим задачу оптимизации:
# Nominal Problemfunction 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скопировать кодМы рассматриваем следующие три набора неопределенностей:
На самом деле, легко понять, что робастная задача потока минимальных затрат (i)(ii) все еще может быть записана как линейная программа, а робастная задача потока минимальных затрат (iii) может быть записана как задача второго порядка. задача оптимизации конуса. Однако здесь мы показываем, что использование пакета JuMPeR может сэкономить этап получения надежного аналога.
# Robust Problemfunction 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, мы случайным образом генерируем
, а затем случайным образом генерируем
Посмотрим на реальную стоимость. Мы берем решение номинальной задачи в качестве эталона и выводим разницу между стоимостью робастного решения и разницей.
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_dataendoutput = Evaluation();# 输出Γ=10时三种uncertainty set给出的differencetemp_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 Приблизительное решение Стабильное число: Положительное полуопределенное программирование/целочисленное программирование
неориентированный графстабильный (независимый)
Множество — это подмножество V, в котором все узлы не связаны друг с другом. Максимальная мощность стабильного множества определяется как α(G), также известный как число стабильности графа.Предполагая |V|=n.Естественно, его целочисленное программирование 0/1 можно записать как:
использовать
Эквивалентно, мы можем получить релаксацию положительной полуопределенной программы:
Пример источника: Пена, Хавьер, Хуан Вера и Луис Ф. Сулуага, «Вычисление числа устойчивости графа с помощью линейного и полуопределенного программирования», 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) = 5StableBin(A14)# 松弛之后的解：5.694362954824882StableSDP(A14)скопировать кодМы указали, что Gurobi еще не поддерживает комбинацию целочисленного программирования и генерации столбцов (т. е. ветвления и цены), но пользовательский разрез поддерживает это (т. е. ветвления и разрезания)! Короче говоря, мы можем добавить только небольшое количество ограничений в начале, и будут добавлены нарушенные ограничения в процессе решения целочисленного программирования, поэтому это также называется ленивым обратным вызовом.
Этот пример исходит из:http://t.cn/ReWaORM
Задача, которую мы хотим решить, — это задача программирования конусов второго порядка с целочисленными ограничениями:
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_callbacksend#产生一个随机样例srand(1234)n = 2c = rand(n)Γ = 50.0# 求解，输出callback次数sol, num_callbacks = solve_ball(c, Γ)println(sol)println(norm(sol))println(num_callbacks)скопировать кодПосле 11 lazycallbacks получаем оптимальное решение. Далее еще один пример взаимодействия, мне тут лень делать анимацию, просто дайте код. Пожалуйста, испытайте эффект на себе в блокноте 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скопировать кодопределить переменную— координаты центра окружности по х и у. Тогда проблема принятия решений может быть решена с помощью следующей задачи невыпуклого нелинейного программирования, см. также: Биргин, Эрнесто Г., Дж. М. Мартинес и Дебора П. Ронкони. Оптимизация упаковка баллонов в прямоугольный контейнер: нелинейный подход». European Journal of Operational Research 160.1 (2005):
19-33.
Подобно идее из предыдущего раздела, мы также используем здесь кусочно-линейные функции и целочисленное программирование для аппроксимации функции L2. Конечно, здесь мы используем более общий подход и используем пакет PiecewiseLinear для облегчения завершения. В частности, мы рассматриваем целочисленное линейное программирование для аппроксимации невыпуклых ограничений в нелинейном программировании следующим образом.
# 画图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")
pendскопировать код# 考虑这样一组问题r = 25d1 = 160d2 = 190k = 11;# MIOfunction 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()endp1,p2,time = PackCirclesGen(r,d1,d2,k,:ZigZagInteger,Inf);# 一个可行解PlotCirclesRectangle(p1,p2)скопировать кодЭту программу также можно использовать для исследования этого вопроса Чжиху:
Сколько кругов диаметром 1 может поместиться в прямоугольнике 2×1000? (https://www.zhihu.com/question/276177865/answer/390216070)
# 不同bivariate分段线性方法的比较r = 21d1 = 120d2 = 120k = 6;Method_List = [:CC,:MC,:DisaggLogarithmic,:SOS2,:Logarithmic,:ZigZag,:ZigZagInteger]# The experiment, note the cutoff time set as 600 secondsfor 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, включая другие материалы о встрече, связанные с Джулией:
1. https://github.com/JuliaCN/MeetUpMaterials/tree/master/Beijing2018
2. https://github.com/JuliaOpt/juliaopt-notebooks/blob/master/notebooks/JuMP-LazyL2Ball.ipynb
Другие оптимизированные пакеты, поддерживаемые JuMP, можно найти здесь: https://www.juliaopt.org/packages/
Например, MultiJuMP.jl, который поддерживает многокритериальную оптимизацию, JuMPChance.jl, который поддерживает случайное ограничение, PolyJuMP.jl, который поддерживает полиномиальную оптимизацию и поддерживает стохастическое динамическое программирование, нелинейное управление, задачи ограничения равновесия и алгоритмы теории графов, пакет для метаэвристики. , и т.д.
Подбор персонала
【Присоединяйтесь к нам】Раздел [Событие] официального аккаунта «Стратегия ИЛИ Стратегия» набирает от одного до двух ответственных редакторов (вы также можете присоединиться в качестве добровольцев) для участия в планировании и написании статей о мероприятиях в этом разделе. Мы надеемся, что заинтересованные друзья присоединятся к нам и узнают друг друга больше друзей-единомышленников, работайте вместе, чтобы обеспечить хорошую платформу для коммуникации событий для аудитории 10W+, и будет определенная награда!
【Для выполнения условий】Я надеюсь, что вы являетесь аспирантом в области промышленной инженерии/исследования операций/управления цепочками поставок/статистики/прикладной математики или выдающимся студентом старших курсов бакалавриата с предпочтительными баллами CET-6 и TOEFL|IELTS! Опыт работы на министерском уровне или выше в студенческом союзе или обществе приветствуется! Приветствуется опыт организации общественных мероприятий!
【связаться с нами】Пожалуйста, присылайте свое резюме на juminaa@163.com, мы с нетерпением ждем вашего присоединения! !
написание эссеОт:Цинь Ханчжан
Ответственный редактор: Ван Юань
Редактор WeChat: Виноград
Источник статьи Отказ от ответственности: эта статья взята из статьи Zhihu: JuMP: Оптимизация моделирования и решения с Юлией ( https://zhuanlan.zhihu.com/p/40807662 )
Если вам нужно перепечатать, пожалуйста, получите инструкции по перепечатке в фоновом режиме общедоступной учетной записи.
Дружеское напоминание: Эта статья составлена и организована Operation Strategy OR Strategy. Она не используется в коммерческих целях. Если содержание нарушает авторские права, мы удалим его в любое время.
Рекомендуемые статьи в прошлом
※[Отчет] Каково работать в Google в области исследования операций?※
※[Интервью] Ху Вухуа из SF Technology: Развитие исследований операций в Китае станет началом весны※
※[Академия] Подробное обсуждение условий KKT※
◆
◆
◆
◆
◆
Если вы занимаетесь исследованием операций | магистром искусственного интеллекта или докторантом, пожалуйста, оставьте сообщение на фоне общедоступной учетной записи ниже: "Добавить группу WeChatСистема подскажет вам дальнейшие требования и шаги присоединения к группе и предложит вам присоединиться к группе исследований глобальных операций или группе исследователей ИИ (в группе собираются лидеры академических кругов и отрасли).
При этом имеем:[Операционные исследования | Оптимизация] [Цепочка поставок | Логистика] [Искусственный интеллект] [Наука о данных | Анализ]В группе QQ тысячи поклонников. Те, кто хочет присоединиться к группе, могут подписаться на общедоступную учетную запись ниже и нажать «Присоединяйтесь к сообществу", чтобы открыть групповой портал.
нажмитечитать оригинал, узнать больше