JuMP: 用Julia進行優化建模及求解

08-11

來自專欄 Optimizers Garden39 人贊了文章

0. 為什麼要用Julia做優化？

本文是我在Julia中文社區2018用戶見面會上關於用Julia進行優化建模及求解的分享內容。

本Tutorial主要介紹JuMP.jl,一種在Julia語言中的開源AML(Algebraic Modeling Language), 類似於AMPL, YALMIP, CVX, Poymo, GAMS等。JuMP的優化建模+solver內存傳遞速度，與其它商業/開源AML的比較見下表。不涉及Convex.jl，類似Matlab中的CVX，由Boyd et al. 開發的另一種Julia中的AML。

該表來自：Dunning, Iain, Joey Huchette, and Miles Lubin. JuMP: A modeling language for mathematical optimization. SIAM Review 59.2 (2017): 295-320. 更多JuMP的實現細節和與其它AML的比較可參

我們發現JuMP在優化建模速度上不輸商業AML，比其它的開源AML要明顯快出許多。

該表來自：https://www.juliaopt.org/

JuMP的另一大好處在於無縫銜接各種商業/開源solver: Gurobi, CPLEX, SCIP, Knitro, Mosek等......

1. 優化建模初步（最基本的線性規劃）

(這部分主要針對優化理論零基礎的讀者，例子來自於——覃含章：運籌學中應該如何理解互補鬆弛性。這條性質又該如何運用?)

假設你是一個木匠，出售手工製作的木頭桌子和木頭椅子。現在你想要制定一個生產計劃，讓這個月的利潤最大化。

注意，生產一張桌子，或者一把椅子，都需要消耗一定數量的木材和時間，而作為一個手工小作坊，你每個月能用來生產的時間是有限的，你每個月能進到的木材也是有限的。

簡單起見，我們假定桌子的利潤固定為一張10元，椅子為一把3元。生產一張桌子需要5單位木材和3單位時間，生產一把椅子需要2單位木材和1單位時間。且我們所有生產的桌子椅子都是能被賣掉的。先假設當月我們總共有200單位木材和90單位時間。

顯然，這個生產規劃問題可以用線性規劃來建模：

$egin{align*} max~ & 10x_1+3x_2 ag{P}\ ext{s.t. } & 5x_1+2x_2leq 200 ag{P1}\ & 3x_1+x_2leq 90 ag{P2} \ & x_1,x_2geq 0. end{align*}$

根據對偶理論，我們也知道等價的對偶問題為：

$egin{align*} min~ & 200p_1+90p_2 ag{D}\ ext{s.t. } & 5p_1+3p_2geq 10 ag{D1}\ & 2p_1+p_2geq 3 ag{D2} \ & p_1,p_2geq 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進行敏感度分析（Sensitivity Analysis）。考慮以下四個（P）的變種問題：（當然，其實都可以手解，用CarpenterPrimal進行驗證）

(1) 把(P)中的(P1)右端項改成250，即木材總量增加。其餘不變。最優解是什麼？

(2) 把(P)中的(P2)右端項改成120，即生產時間增加。其餘不變。最優解是什麼？增加的利潤（相比原來的(P)）和（P2）的影子價格的關係？

(3) 假設桌子的單位利潤不變，椅子的單位利潤增加到多少的時候(P)的最優解開始生產椅子（而不生產桌子）？

(4) 假設時間的總量仍為90，木材的總量減少到多少的時候木材的影子價格嚴格大於0（時間的影子價格反而變成0）？

如對這些基本概念不熟悉，可參閱本節一開始提到的我的回答。

2. 線性規劃中的Column Generation實例：大規模additive convex regression

給定數據 $Yin mathbb{R}^n, Xin mathbb{R}^{n imes d},$ 我們考慮以下線性規劃問題：

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

注意我們如果在目標函數中再額外加上正則項 $lambdasum_{j=1}^d left| sum_{i=1}^n f_j(X_{ij}) ight|$ （ $lambda>0$ ）, 則我們可以得到一個稀疏的(sparse) additive model。詳見：Ravikumar, Pradeep, et al. "Sparse additive models." Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71.5 (2009): 1009-1030.

上述模型是線性規劃因為等價於(記 $[:]_j$ 為 $X_{:j}$ 升序排序後的對應原index的序列)

$egin{align*} min~ & sum_{i=1}^n r_i^++r_i^-\ ext{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 = 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 2n-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

接下來我們介紹什麼是Column Generation，並展示如何對上述問題利用JuMP和Gurobi進行Column Generation。這種方法是基於單純形法（simplex methopd）的。即我們在一開始默認所有的變數都為0且不是基變數，而在每步iteration的時候利用前一輪的對偶變數的值計算reduced cost(只需要計算一步矩陣和向量的乘法)，並將reduced cost最負的變數加入（成為基變數）。同時注意到，因為原問題的約束矩陣實際上具有分塊特性（按照 $j=1,…,d$ 分塊），實際上我們可以採取Dantzig-Wolfe的formulation，即利用一個cyclic的update方式，在每個iteration對每個block單獨計算reduced cost（比一起算快了 $d$ 倍）。當所有reduced cost非負的時候，找到最優解：演算法終止。

對Dantzig-Wolfe或者column generation不太了解的同學，請見一些經典的資料，如：https://perso.uclouvain.be/anthony.papavasiliou/public_html/DantzigWolfe.pdf

當然，第五個例子里的constraint generation其實原理和column(variable) generation類似，那個例子對於零背景知識的同學可能更容易看。

# 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, " ") push!(NewColumns[k],Δ_temp) push!(IdxList[k],idx) solve(M_conv) iter = iter + 1 dual_var = getdual(constr); else #print("Skip covariate ", k," ") flag[k] = 0 end end end #print("CG ends after ", iter-1, " itertaions. ") # 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])

接下來我們比較這兩種formulation的求解速度。

@benchmark CG_ConvReg(200)

@benchmark Full_ConvReg(600)

令人驚訝的是，基於古老的單純形法和Column Generation的演算法竟然要遠快於Gurobi求解LP默認的non-deterministic concurrent方法！這證明針對特定的優化問題，有時候更聰明的建模會給你更快的演算法。

3. 魯棒線性規劃（Robust Linear Programming）

我們考慮最小費用流問題（minimum cost flow），其中每條邊上的cost具有不確定性（uncertainty）。即我們考慮一個有向圖 $G=(V,E)，V={1,…,n}$ ,我們考慮優化問題：

$egin{align*} min~ & sum_{(i,j)in E} c_{ij}x_{ij}\ ext{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~jin V ackslash {1,n}\ & 0leq x_{ij}leq C_{ij}. end{align*}$

# 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

我們考慮以下三種uncertainty set:

$egin{align*} & (i)~ ext{Box Uncertainty}: mathcal{U}_{L_infty} = { extbf{c}:μ_{ij}-delta_{ij}gamma_{ij}leq c_{ij}leq mu_{ij}+delta_{ij}gamma_{ij}, |gamma|_{infty}leq Gamma }\ & (ii)~ ext{Polyhedral Uncertainty}: mathcal{U}_{L_1} = { extbf{c}:μ_{ij}-delta_{ij}gamma_{ij}leq c_{ij}leq mu_{ij}+delta_{ij}gamma_{ij}, |gamma|_1leq Gamma }\ & (iii)~ ext{Elliopsoidal Uncertainty}: mathcal{U}_{L_2} = { extbf{c}:μ_{ij}-delta_{ij}gamma_{ij}leq c_{ij}leq mu_{ij}+delta_{ij}gamma_{ij}, |gamma|_2leq Gamma }\ end{align*}$

其實，容易知道(i)(ii)的robust min cost flow問題仍可寫成一個線性規劃，(iii)的robust min cost flow問題可寫成一個二階錐規劃問題(second-order cone optimization)。不過我們這裡展示利用JuMPeR包可以省去推導robust counterpart的步驟。

# 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

我們這裡做一些模擬來看一下robust solution的表現。假設G是完全連同的，所有邊的容量都是0.5，我們隨機產生 $μ_{ij}～Uniform(0,10)$ $δ_{ij}～Uniform(0,μ_{ij})$ ，然後我們再隨機產生 $c_{ij}$ 來看實際的cost的情況。我們以nominal problem的解為benchmark,輸出robust solution的cost與其的差。

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")

那麼我們確實看到box constraint給出了更robust的solution，而L1/2 uncertainty構建下的魯棒優化問題更mild一些。

JuMPeR也可以model多階段的魯棒優化問題，尤其比如是affinely-adjustable policy的(這比uncertainty set的formulation要松的多)。這裡只是提一下，不詳細介紹了，有興趣的同學可以自行研究，其實也是很有用的。

4. 近似解Stable Number: 半正定規劃/整數規劃

一個無向圖 $G=(V,E)$ 的stable(independent) set是V的一個子集，其中所有節點互相都不相連。 stable set的最大cardinality定義為α(G)，也叫做圖的stability number. 假設|V|=n. 自然，其0/1整數規劃可以寫成：

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

利用 $x_i+x_j≤1$ 和 $x_ix_j=0$ 等價，我們可以得到一個半正定規劃的relaxation：

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

例子來源：Pena, Javier, Juan Vera, and Luis F. Zuluaga. "Computing the stability number of a graph via linear and semidefinite programming." 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)

5. Lazy Constriant Generation: 在整數規劃中用分段線性函數近似L2球

我們指出Gurobi還暫不支持整數規劃和column generation結合（即branch and price）,然而自定義user cut確是支持的（即branch and cut）！簡單來說，即我們可以一開始只加入少量的約束，在一邊求解整數規劃的過程中將被violate的約束加入，因此也叫lazy callback。

這個例子來自於：https://github.com/JuliaOpt/juliaopt-notebooks/blob/master/notebooks/JuMP-LazyL2Ball.ipynb

我們要求解的問題是一個整數約束的二階錐規劃問題：

$egin{align*} max ~ & c^T x\ ext{s.t. } & |x|_2leq Gamma,\ & xin 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_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次lazycallback我們得到最優解。接下來還有個交互的例子，我這邊懶得做動畫了，只給code。效果請自行在jupyter notebook里體驗。實際上就是說可以手動在bar上調整參數看到不同情況下的最優解和可行域在二維平面的圖像。

# 這裡展示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} )$ 為圓心的x和y坐標. 那麼決策問題可以用如下非凸非線性規劃問題來求解, 也見：Birgin, Ernesto G., J. M. Mart?nez, and Débora P. Ronconi. "Optimizing the packing of cylinders into a rectangular container: a nonlinear approach." European Journal of Operational Research 160.1 (2005): 19-33.

$egin{align} min~ & sum_{i eq jin {1,ldots, k} } maxleft(0,2r - sqrt{(p_i^1-p_j^1)^2+(p_i^2-p_j^2)^2} ight)\ ext{s.t. } & rleq p_1^ileq d_1-r,forall~ i=1,ldots,k, \ & rleq p_2^ileq d_2-r,forall~ i=1,ldots,k. \ end{align}$

和上一節思路類似，我們這裡也利用分段線性函數和整數規劃來近似這裡的L2函數。當然，我們這裡採取更一般的做法，利用PiecewiseLinear包來方便完成。具體來說，我們考慮用整數線性規劃近似如下非線性規劃中的非凸約束。

$egin{align} min~ & sum_{i eq jin {1,ldots, k} } z_{ij}\ ext{s.t. } & z_{ij}geq 2r - sqrt{(p_i^1-p_j^1)^2+(p_i^2-p_j^2)^2}, forall ~ i eq jin {1,ldots,k},\ &rleq p_1^ileq d_1-r,forall~ i=1,ldots,k, \ & rleq p_2^ileq d_2-r,forall~ i=1,ldots,k. \ & z_{ij}geq 0, forall ~ i eq jin {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") 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/5j-i)*r*2, 0:r/5j-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)

這個程序也可以作為這個知乎問題的探索：

2×1000的矩形中能放多少個直徑為1的圓？?

www.zhihu.com

# 不同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! ")end

我們注意到不同的線性分段逼近方法的計算速度是大不一樣的，ZigZag的速度最快，CC/MC這樣的formulation一般不建議用。詳細的理論方法請見：Huchette, Joey, and Juan Pablo Vielma. "Nonconvex piecewise linear functions: Advanced formulations and simple modeling tools."

7. 更多

這篇Tutorial的GitHub鏈接，包括meetup上其它的julia相關資料：

JuliaCN/MeetUpMaterials?

github.com

其它JuMP支持的優化拓展包可見：https://www.juliaopt.org/packages/

比如，支持多目標優化的MultiJuMP.jl，機會約束（chance constraint）的JuMPChance.jl，多項式（polynomial）優化的 PolyJuMP.jl，還有支持隨機動態規劃，非線性控制，平衡約束問題，圖論演算法，元啟發式演算法等等的package。

之前我也寫過一個Julia的安利文：

覃含章：Julia：高效易用的數值計算/優化編程語言?

zhuanlan.zhihu.com

多練，多嘗試...Practice makes Perfect!

任何問題，或者交流，歡迎聯繫我！