0. Why do you want to optimize Julia?
This article is what I shared about optimizing modeling and solving with Julia at the 2018 User meeting of Julia Chinese Community.
JuMP. Jl, an open source AML(Algebraic Modeling Language) in Julia, is similar to AMPL, YALMIP, CVX, Poymo, GAMS, etc. JuMP’s optimized modeling + Solver memory transfer speed, compared to other commercial/open source AML, is shown in the table below. Jl is not involved, similar to CVX in Matlab. Boyd et al. Another type of AML developed in Julia.
We found that JuMP was as fast at optimizing modeling as commercial AML and significantly faster than other open source AML.
Another great benefit of JuMP is the seamless connection between various commercial/open source Solver: Gurobi, CPLEX, SCIP, Knitro, Mosek, etc……
1. Preliminary optimization modeling (most basic linear programming)
This part is mainly aimed at zero-based readers of optimization theory. Examples come from qin Hanzhang: How to understand complementary relaxation in operations research. How do we apply this property?
Imagine you are a carpenter, selling handmade wooden tables and wooden chairs. Now you want to make a production plan to maximize your profits for the month.
Note that it takes a certain amount of wood and time to produce a table, or a chair, and as a small workshop, you only have a limited amount of time to produce each month, and so does the amount of wood you can get in each month.
For simplicity, let’s assume a fixed profit of $10 for a table and $3 for a chair. It takes 5 units of wood and 3 units of time to produce a table, and 2 units of wood and 1 unit of time to produce a chair. And all of our tables and chairs are sellable. Let’s say we have 200 units of wood and 90 units of time for the month.
Obviously, this production planning problem can be modeled using linear programming:
According to the duality theory, we also know that the equivalent duality problem is:
Function CarpenterPrimal(c,A,b) Primal = Model(solver = GurobiSolver(OutputFlag = 0)) @variable(Primal, x[1:2]>=0) # Define inequality constr = @constraint(Primal, A*x.<=b) # define objective function @objective(Primal, Max, dot(c,x)) Return getobJectivevalue (Primal), getValue (x), getDual (constr) endCopy the codeCall CarpenterPrimal and the optimal solution is to produce 30 tables and no chairs:
We point out that JuMP can also be used for Sensitivity Analysis. Consider the following four variants of (P) :(of course, all can be solved manually and verified with CarpenterPrimal)
(1) Change the right end of (P1) in (P) to 250, that is, the total amount of wood increases. The rest is the same. What is the optimal solution? (2) Change the right end of (P2) in (P) to 120, that is, the production time increases. The rest is the same. What is the optimal solution? What is the relationship between increased profit (compared to original (P)) and shadow price of (P2)? (3) Assuming the unit profit of the table remains constant, to what extent does the unit profit of the chair increase, at which point does the optimal solution of (P) start producing chairs (instead of tables)? (4) Assuming that the total amount of time is still 90, when the total amount of wood is reduced to what extent, the shadow price of wood is strictly greater than 0 (the shadow price of time becomes 0 instead)?
If you are not familiar with these basic concepts, please refer to my answers mentioned at the beginning of this section.
2. The Column Generation instance in linear programming: large-scale additive concave regression
Given data
Notice how we add the regular term to the objective function
The above model is a linear program because it is equivalent to (n
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]) endCopy the codeFunction 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) endCopy the codeNext we introduce what Column Generation is and show you how to use JuMP and Gurobi for Column Generation. This method is based on simplex methopd. That is, all variables are 0 and not base variables by default at the beginning. Instead, at each iteration step, the reduced cost is calculated using the value of the dual variable of the previous iteration (only the matrix and vector multiplication step is needed), and the variable with the most negative reduced cost is added (becoming the base variable). It is also noted that because the constraint matrix of the original problem actually has the block property (according to
The Dantzig – Wolfe or column generation do not know much about the students, please see some classic information, such as: perso. Uclouvain. Be/Anthony. The pap…
Of course, the column(variable) generation constraint is similar to column(variable) generation.
Column generationization linear programming model function CG_ConvReg(Max_T) # build the LO: No δ included M_conv = Model(solver=GurobiSolver(TimeLimit = Max_T,OutputFlag = 0,Method = 0) 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 the column generation L [: 1, k] = ones (n, 1)] [:, 2, k = L - 'ones (n, 1)] [2, 3, k = L XX XX [1] [2] - [] 2, 4, k = L - XX XX [1] [2] + the 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 Constr_list = constrr [:,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 endCopy the codephi = 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])Copy the code
Next, we compared the solution speeds of these two formulations.
@benchmark CG_ConvReg(200)
Copy the code
@benchmark Full_ConvReg(600)
Copy the code
Surprisingly, the algorithms based on the old simplex method and Column Generation are much faster than Gurobi’s non-Deterministic Concurrent method for LP’s default solution! This proves that for a particular optimization problem, sometimes smarter modeling will give you a faster algorithm.
3. Robust Linear Programming
We consider a minimum cost flow in which the cost on each edge has uncertainty. That is, we consider a directed graph
# 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
Copy the codeLet’s consider the following three kinds of uncertainty sets:
In fact, the robust min cost flow problem easy to know (I) and (II) can still be written as a linear programming, (iii) The robust min cost flow problem can be written as a second-order cone optimization problem. But we show here that using the JuMPeR package saves the procedure of deriving the robust counterpart.
# RobustProblem 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) = = u + r [I, j] [I, j] * delta [I, j]) end end @ the 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) endCopy the codeHere we do some simulations to see the performance of robust solution. Assuming that G is completely integrated, and that all edges have a capacity of 0.5, we randomly generate
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 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, mu, the delta, Γ, Inf) sol_1 = RobustProblem (n, mu, the delta, Γ, 1) sol_2 = RobustProblem (n, mu, the delta, Γ, 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(); Output [8,:],xlabel = "n", label = "Uinf") plot = uncertainty set [10,25,50,75,100] (temp_p,25,50,75,100 [10], the output (9, :), label = "U1)" the plot! (temp_p,25,50,75,100 [10], the output (10, :), label = "U2")Copy the code
So we do see that box Constraint gives a more robust solution, while the robust optimization problem under L1/2 Uncertainty construction is more mild.
JuMPeR can also model multi-stage robust optimization problems, especially if it is using an affhydrocarbon policy (which is much looser than uncertainty set formulation). Here just mentioned, not detailed introduction, interested students can study by themselves, in fact, is also very useful.
4. Approximately solve Stable Number: semidefinite programming/integer programming
An undirected graph
using
Examples: Pena, Javier, Juan Vera, and Luis F. Zuluaga. “Computing the stability number of a graph via linear and semidefinite programming.” SIAM Journal Journal of 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
Copy the codefunction 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
Copy the code
Take this G, for example.
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; 1 0 0 0 1 1 1 0 0 0 0 1 0 0; 1 0 0 0 1 1 0 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 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 1 0 0 0 0; 1 0 0 0 0 0 0 0 1 1 1 0 0 0; 1 1 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 0 1; 0 0 0 0 1 1 0 0 0 0 0 0 0, 1]; # α(G14) = 5 StableBin(A14) # StableSDP(A14)Copy the code5. Lazy Constriant Generation: Use piecewise linear functions to approximate L2 spheres in integer programming
We pointed out that Gurobi does not yet support integer programming in combination with column Generation (branch and Price), while custom user cuts do support branch and cut! To put it simply, we can add only a few constraints at the beginning, and the constraints that violate will be added in the process of solving integer programming, which is also called lazy callback.
This example comes from: github.com/JuliaOpt/ju…
The problem we ask to solve is a second-order cone programming problem with integer constraints:
Function solve_ball(c, γ, ϵ=1e-6) n = length(c) m = Model(solver=GurobiSolver(OutputFlag=0)) # A box @variable(m, -γ ≤ x[1:n] ≤ γ, Int) @objective(m, Max, dot(c,x)) Define callback function, Num_callbacks = 0 function norm_callback(cb) num_callbacks += 1 N = getValue (x) # norm_callbacks = norm(N) # If L ≤ γ + ϵ return end # if not, add cut! Notice that this cut will make x infeasible, @lazyconstraint(cb, dot(N,x) ≤ γ *L) end Srand (1234) n = 2 c = rand(n) γ = 50.0 # Solve (m) return getValue (x), num_callbacks end # solve(m) return getValue (x), num_callbacks end Sol, num_callbacks = solve_ball(C, γ) println(norm(sol)) println(num_callbacks)Copy the codeAfter 11 times of Lazycallback, we get the optimal solution. Here’s another example of interaction, and I’m not going to animate it here, I’m just going to give code. Please experience the effect in Jupyter Notebook. In fact, it means that you can manually adjust the parameters on bar and see the optimal solution in different cases and the image of the feasible region in two dimensions.
Interact and Compose packages are shown here. 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 + SQL 300,0.5 + [1] / sol [2] / 300)]), Compose. The stroke (" black ")), Compose (context (), Circle ((0.5 + (100 / norm (sol)) * sol / 300)... , 0.02) and the fill (" red ")), compose (context (), circle (0.5, 0.5, 0.333), the fill (" lightblue "))) endCopy the code6. General nonlinear integer programming: Dense packing ball problem
Let’s consider the following question: I have k balls of radius R and a rectangle box of size D1 by D2. Is there a way to fit these balls into the rectangle? Obviously, a more general problem is to find the largest k for any rectangle. However, this general problem can be solved by solving the decision problem before the sequential solution (increasing k), so we will discuss the decision problem given k.
Define variables
Just like in the last video, we use piecewise linear functions and integer programming to approximate L2. Of course, we take a more general approach here and use the PiecewiseLinear package to facilitate the completion. In particular, we consider using integer linear programming to approximate non-convex constraints in nonlinear programming as follows.
Function rectangle (p1,p2) p = scatter(p2),(p1)) for I = 1:k plot! (p, (p2 [I]) + r * cos. (linspace (0, 2 * PI, 100)), (p1 [I]) + r * sin. (linspace (0, 2 * PI, 100)), color = "blue", legend = false, aspect_ratio = 1 ) end plot! (p, linspace (0, d2, 100), zeros (100, 1), color = "red") the plot! (p, zeros (100, 1), linspace (0, d1, 100), color = "red") the plot! (p, linspace (0, d2, 100), 'ones (100, 1) * d1, color = "red") the plot! (p,d2*ones(100,1),linspace(0,d1,100),color = "red"Copy the codeR = 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]) # approximate nonconvex constraint fun_dist = using the PiecewiseLinearOpt package Piecewiselinear (PackCircles, [I, j] s1 and s2 (I, j), - (j - I) * r * 2:5: (r/j - I) * r * 2, 0:5: (r/j - I) * r * 2, (u, v) - > (u v ^ ^ 2 + 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)Copy the code
This program can also be used as an exploration of this zhihu problem:
How many circles of diameter 1 can fit in a 2 by 1000 rectangle?
# Comparison of different bivariate piecewise linear methods 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") endCopy the code
We noticed that the calculation speeds of different linear piecewise approximation methods were quite different, ZigZag was the fastest, and formulation like CC/MC was generally not recommended. For detailed theoretical methods, see Huchette, Joey, and Juan Pablo Vielma. “Nonconvex piecewise Linear functions: Advanced formulations and simple modeling tools.”
7. More
GitHub links to this Tutorial, including other Julia related materials on Meetup:
JuliaCN/MeetUpMaterials
Other JuMP support optimization development package: www.juliaopt.org/packages/
Jl supports multi-objective optimization, Jumpchance. jl supports opportunity constraint, and Polyjum. jl supports polynomial optimization, as well as random dynamic programming, nonlinear control, and equilibrium constraint problems. Package of graph theory algorithms, meta-heuristic algorithms and so on.
I also wrote a Julia’s Amway before:
An efficient and easy-to-use numerical calculation/optimization programming language
Practice more, try more… Practice makes Perfect!
Any questions or exchanges, please feel free to contact me!