Skip to content
This repository has been archived by the owner on Apr 11, 2020. It is now read-only.

Commit

Permalink
update notebooks to use GLPK and Ipopt
Browse files Browse the repository at this point in the history
  • Loading branch information
@joaquimg
joaquimg committed Jun 29, 2018
1 parent e84d479 commit 9ba8b16
Show file tree
Hide file tree
Showing 10 changed files with 5,313 additions and 1,010 deletions.
135 changes: 74 additions & 61 deletions notebooks/Chiwei Yan - Cutting Stock.ipynb
View file
Expand Up @@ -106,8 +106,16 @@
],
"source": [
"#import necessary packages and define model\n",
"# Load JuMP\n",
"using JuMP\n",
"master = Model()\n",
"using MathOptInterface\n",
"# Load solver package\n",
"using MathOptInterfaceXpress\n",
"# shortcuts\n",
"const MOI = MathOptInterface\n",
"const MOIU = MathOptInterface.Utilities\n",
"\n",
"master = Model(optimizer = XpressOptimizer())\n",
"\n",
"#If Gurobi is installed (requires license), you may uncomment the code below to switch solvers\n",
"#using Gurobi\n",
Expand Down Expand Up @@ -179,7 +187,7 @@
],
"source": [
"#define initial variables\n",
"@defVar(master, x[1:2] >= 0)\n",
"@variable(master, x[1:2] >= 0)\n",
"\n",
"#width\n",
"w=[14 31 36 45]\n",
Expand All @@ -190,15 +198,16 @@
"\n",
"\n",
"#define constraint references\n",
"@defConstrRef myCons[1:4]\n",
"myCons = Any[]\n",
"\n",
"#define constraints\n",
"for i=1:4\n",
" myCons[i] = @addConstraint(master, dot(x, vec(A[i,:]))>=b[i])\n",
" myCon = @constraint(master, dot(x, vec(A[i,:]))>=b[i])\n",
" push!(myCons, myCon)\n",
"end\n",
"\n",
"#define objective\n",
"@setObjective(master, Min, sum(x))\n",
"@objective(master, Min, sum(x))\n",
"\n",
"master"
]
Expand All @@ -224,8 +233,9 @@
}
],
"source": [
"status=solve(master)\n",
"getValue(x)\n",
"JuMP.optimize(master)\n",
"status = JuMP.terminationstatus(master)\n",
"JuMP.resultvalue.(x)\n",
"\n",
"#get the optimal solution\n",
"println(\"\\nOptimal Solution is:\\n\")\n",
Expand All @@ -236,8 +246,8 @@
"\n",
"for i=1:size(A,2)\n",
" \n",
" if getValue(x[i])>epsilon \n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", getValue(x[i]))\n",
" if JuMP.resultvalue(x[i])>epsilon \n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(x[i]))\n",
" end\n",
"end"
]
Expand Down Expand Up @@ -286,23 +296,24 @@
"source": [
"r=[getDual(myCons)[1:4]]\n",
"\n",
"sub = Model() \n",
"sub = Model(optimizer = XpressOptimizer()) \n",
"\n",
"#width\n",
"w=[14,31,36,45]\n",
"\n",
"#define cutting pattern variables\n",
"@defVar(sub, a[1:4]>=0, Int)\n",
"@variable(sub, a[1:4]>=0, Int)\n",
"\n",
"#define feasible cutting constraint\n",
"@addConstraint(sub, dot(w,a)<=100)\n",
"@constraint(sub, dot(w,a)<=100)\n",
"\n",
"#define objective\n",
"@setObjective(sub, Max, dot(r,a))\n",
"@objective(sub, Max, dot(r,a))\n",
"\n",
"sub\n",
"\n",
"status=solve(sub)\n",
"JuMP.optimize(master)\n",
"status = JuMP.terminationstatus(master)\n",
"\n",
"#print new cutting pattern\n",
"pattern=[getValue(a)[1:4]]\n",
Expand Down Expand Up @@ -401,7 +412,8 @@
],
"source": [
"#column-wise adding new variable z\n",
"@defVar(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)\n",
"# TODO column addtion\n",
"@variable(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)\n",
"\n",
"#look at the master problem again\n",
"master"
Expand Down Expand Up @@ -430,7 +442,8 @@
],
"source": [
"#solve the master problem again\n",
"status=solve(master)\n",
"JuMP.optimize(master)\n",
"status = JuMP.terminationstatus(master)\n",
"\n",
"#get the optimal solution\n",
"println(\"\\nOptimal Solution is:\\n\")\n",
Expand All @@ -439,14 +452,14 @@
"\n",
"for i=1:length(x)\n",
" \n",
" if getValue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", getValue(x[i]))\n",
" if JuMP.resultvalue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(x[i]))\n",
" end\n",
"end\n",
"\n",
"\n",
"if getValue(z)>epsilon\n",
" println(\"Cutting Pattern: \", int(pattern), \", Number of Paper Rolls Cut Using this Pattern: \", getValue(z))\n",
"if JuMP.resultvalue(z)>epsilon\n",
" println(\"Cutting Pattern: \", int(pattern), \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(z))\n",
"end\n"
]
},
Expand Down Expand Up @@ -535,51 +548,51 @@
"#master = Model(solver=GurobiSolver(Method=0)) # Switch LP algorithm to Primal Simplex, in order to enjoy warm start\n",
"\n",
"#define initial variables\n",
"@defVar(master, x[1:2] >= 0)\n",
"@variable(master, x[1:2] >= 0)\n",
"\n",
"#constraint coefficient for initial variables\n",
"A=[1 0; 1 0; 0 2; 1 0]\n",
"b=[211; 395; 610; 97]\n",
"\n",
"#define constraint references (why?)\n",
"@defConstrRef myCons[1:4]\n",
"@constraint myCons[1:4]\n",
"\n",
"#define constraints\n",
"for i=1:4\n",
" myCons[i] = @addConstraint(master, dot(x, vec(A[i,:]))>=b[i])\n",
" myCons[i] = @constraint(master, dot(x, vec(A[i,:]))>=b[i])\n",
"end\n",
"\n",
"#define objective\n",
"@setObjective(master, Min, sum(x))\n",
"@objective(master, Min, sum(x))\n",
"\n",
"#solve master problem\n",
"solve(master)\n",
"JuMP.optimize(master)\n",
"\n",
"println(\"Iteration 1, Master Problem Objective Value:\", getObjectiveValue(master))\n",
"\n",
"#subproblem to iteratively generate new columns\n",
"\n",
"#get optimal dual prices from the master problem\n",
"r=[getDual(myCons)[1:4]]\n",
"r=[JuMP.resultdual(myCons)[1:4]]\n",
"\n",
"sub=Model() \n",
"sub=Model(optimizer = XpressOptimizer()) \n",
"\n",
"#width\n",
"w=[14,31,36,45]\n",
"\n",
"#define cutting pattern variables\n",
"@defVar(sub, a[1:4]>=0, Int)\n",
"@variable(sub, a[1:4]>=0, Int)\n",
"\n",
"#define feasible cutting constraint\n",
"@addConstraint(sub, dot(w,a)<=100)\n",
"@constraint(sub, dot(w,a)<=100)\n",
"\n",
"#define objective\n",
"@setObjective(sub, Max, dot(r,a))\n",
"@objective(sub, Max, dot(r,a))\n",
"\n",
"#solve the subproblem\n",
"solve(sub)\n",
"JuMP.optimize(sub)\n",
"\n",
"sub_obj=getObjectiveValue(sub);\n",
"sub_obj=JuMP.objectivevalue(sub);\n",
"\n",
"epsilon=1e-6; \n",
"\n",
Expand All @@ -593,10 +606,10 @@
"while sub_obj>1+epsilon #why?\n",
"\n",
" #cutting pattern (constraint coefficients) for the new variable\n",
" pattern=getValue(a)[1:4]\n",
" pattern=JuMP.resultvalue(a)[1:4]\n",
" \n",
" #column-wise adding new variable z\n",
" @defVar(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)\n",
" @variable(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)\n",
" \n",
" println(\"\\tAdd a new variable with cutting pattern: \", pattern, \", reduced cost: \", (1-sub_obj))\n",
" \n",
Expand All @@ -605,19 +618,19 @@
" #add new cutting pattern to pattern list\n",
" append!(A_new, pattern)\n",
" \n",
" solve(master)\n",
" JuMP.optimize(master)\n",
" \n",
" println(\"\\nIteration \",iter, \", Master Problem Objective Value:\", getObjectiveValue(master))\n",
" println(\"\\nIteration \",iter, \", Master Problem Objective Value:\", JuMP.objectivevalue(master))\n",
" \n",
" #get new optimal dual prices\n",
" r=[getDual(myCons)[1:4]]\n",
" r=[JuMP.resultdual(myCons)[1:4]]\n",
" \n",
" #modify the objective of the subproblem based on new dual prices\n",
" @setObjective(sub, Max, dot(r,a))\n",
" @objective(sub, Max, dot(r,a))\n",
" \n",
" solve(sub)\n",
" JuMP.optimize(sub)\n",
" \n",
" sub_obj=getObjectiveValue(sub)\n",
" sub_obj=JuMP.objectivevalue(sub)\n",
" \n",
" iter=iter+1\n",
" \n",
Expand All @@ -632,15 +645,15 @@
"\n",
"for i=1:length(x)\n",
" \n",
" if getValue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", getValue(x[i]))\n",
" if JuMP.resultvalue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(x[i]))\n",
" end\n",
"end\n",
"\n",
"for i=1:length(newColumns)\n",
" \n",
" if getValue(newColumns[i])>epsilon\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", getValue(newColumns[i]))\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(newColumns[i]))\n",
" end\n",
"end"
]
Expand All @@ -659,7 +672,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"###How to obtain INTEGER solution? ###\n",
"### How to obtain INTEGER solution? ###\n",
"\n",
"We solve the LP relaxation successfully using column generation, however the original cutting stock problem is an integer program. Can we apply column generation to obtain optimal integer solution? \n",
"\n",
Expand All @@ -672,7 +685,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"####Method 1: Rounding####\n",
"#### Method 1: Rounding####\n",
"\n",
"Rounding a fractional solution to its _nearest_ and _feasible_ is a common heuristic for solving integer program. It's pretty problem specific. In cutting stock problem, we observe that if we _round up_ all the fractional solutions, feasibility will maintain. Thus we get our first integer solution:"
]
Expand Down Expand Up @@ -708,17 +721,17 @@
"\n",
"for i=1:length(x)\n",
" \n",
" if getValue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", ceil(getValue(x[i])))\n",
" summation=summation+ceil(getValue(x[i]))\n",
" if JuMP.resultvalue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", ceil(JuMP.resultvalue(x[i])))\n",
" summation=summation+ceil(JuMP.resultvalue(x[i]))\n",
" end\n",
"end\n",
"\n",
"for i=1:length(newColumns)\n",
" \n",
" if getValue(newColumns[i])>epsilon\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", ceil(getValue(newColumns[i])))\n",
" summation=summation+ceil(getValue(newColumns[i]))\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", ceil(JuMP.resultvalue(newColumns[i])))\n",
" summation=summation+ceil(JuMP.resultvalue(newColumns[i]))\n",
" end\n",
"end\n",
"\n",
Expand All @@ -736,7 +749,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"####Method 2: Branch-and-Bound on Root Node ####\n",
"#### Method 2: Branch-and-Bound on Root Node ####\n",
"\n",
"It is troublesome to implement column generation on every node of the branch and bound tree. A common industry / research practice is to directly branch-and-bound the model only with columns generated from solving the LP relaxation. This is a heuristic because optimal set of cutting patterns for the IP might not be the same as the LP relaxation, i.e. we might lose some \"good columns\" to reach optimal integer solution. The upside is, it is very easy to implement with commercial solvers."
]
Expand Down Expand Up @@ -787,7 +800,7 @@
" setCategory(newColumns[i],:Int)\n",
"end\n",
"\n",
"solve(master)\n",
"JuMP.optimize(master)\n",
"\n",
"\n",
"#print optimal solution\n",
Expand All @@ -802,17 +815,17 @@
"\n",
"for i=1:length(x)\n",
" \n",
" if getValue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", getValue(x[i]))\n",
" summation=summation+getValue(x[i])\n",
" if JuMP.resultvalue(x[i])>epsilon\n",
" println(\"Cutting Pattern: \", A[:,i], \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(x[i]))\n",
" summation=summation+JuMP.resultvalue(x[i])\n",
" end\n",
"end\n",
"\n",
"for i=1:length(newColumns)\n",
" \n",
" if getValue(newColumns[i])>epsilon\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", getValue(newColumns[i]))\n",
" summation=summation+getValue(newColumns[i])\n",
" if JuMP.resultvalue(newColumns[i])>epsilon\n",
" println(\"Cutting Pattern: \", int(A_new[:,i]), \", Number of Paper Rolls Cut Using this Pattern: \", JuMP.resultvalue(newColumns[i]))\n",
" summation=summation+JuMP.resultvalue(newColumns[i])\n",
" end\n",
"end\n",
"\n",
Expand All @@ -838,15 +851,15 @@
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.11",
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.3"
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.3.11"
"version": "0.6.0"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 9ba8b16

Please sign in to comment.