diff --git a/.gitignore b/.gitignore index 9ff4c2a9acd..c769da3c251 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ *.sublime* *.opf *.cov +examples/Manifest.toml diff --git a/.travis.yml b/.travis.yml index 01c608b653f..1e1e6e92687 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,3 +29,9 @@ jobs: - julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.add(PackageSpec(path=pwd()))' - julia --project=docs --color=yes docs/make.jl after_success: skip + - stage: "Examples" + julia: 1.0 + os: linux + script: + - julia --project=examples -e 'using Pkg; Pkg.instantiate(); Pkg.add(PackageSpec(path=pwd()))' + - julia --project=examples --color=yes examples/run_examples.jl diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 00000000000..d622aa683fa --- /dev/null +++ b/examples/Project.toml @@ -0,0 +1,4 @@ +[deps] +GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" diff --git a/examples/basic.jl b/examples/basic.jl index a47f942349f..80ccb1c5acb 100644 --- a/examples/basic.jl +++ b/examples/basic.jl @@ -4,32 +4,51 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# basic.jl -# -# Solves a simple LP: -# max 5x + 3y -# st 1x + 5y <= 3 -# 0 <= x <= 2 -# 0 <= y <= 30 -############################################################################# -using JuMP, Clp +using JuMP, GLPK, Test + +""" + example_basic([; verbose = true]) + +Formulate and solve a simple LP: + max 5x + 3y + st 1x + 5y <= 3 + 0 <= x <= 2 + 0 <= y <= 30 + +If `verbose = true`, print the model and the solution. +""" +function example_basic(; verbose = true) + model = Model(with_optimizer(GLPK.Optimizer)) + + @variable(model, 0 <= x <= 2) + @variable(model, 0 <= y <= 30) + + @objective(model, Max, 5x + 3y) + @constraint(model, 1x + 5y <= 3.0) -m = Model(with_optimizer(Clp.Optimizer)) + if verbose + print(model) + end -@variable(m, 0 <= x <= 2) -@variable(m, 0 <= y <= 30) + JuMP.optimize!(model) -@objective(m, Max, 5x + 3y) -@constraint(m, 1x + 5y <= 3.0) + obj_value = JuMP.objective_value(model) + x_value = JuMP.value(x) + y_value = JuMP.value(y) -print(m) + if verbose + println("Objective value: ", obj_value) + println("x = ", x_value) + println("y = ", y_value) + end -JuMP.optimize!(m) + @test obj_value ≈ 10.6 + @test x_value ≈ 2 + @test y_value ≈ 0.2 +end -println("Objective value: ", JuMP.objective_value(m)) -println("x = ", JuMP.value(x)) -println("y = ", JuMP.value(y)) +example_basic(verbose = false) diff --git a/examples/cannery.jl b/examples/cannery.jl index 3091bdd1637..f48f16a783c 100644 --- a/examples/cannery.jl +++ b/examples/cannery.jl @@ -1,79 +1,75 @@ -############################################################################# # Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# JuMP implementation of the cannery problem -# Dantzig, Linear Programming and Extensions, -# Princeton University Press, Princeton, NJ, 1963. -# -# Author: Louis Luangkesorn -# Date January 30, 2015 +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -using JuMP, Clp +using JuMP, GLPK, Test const MOI = JuMP.MathOptInterface -solver = Clp.Optimizer - -function PrintSolution(is_optimal, plants, markets, ship) - println("RESULTS:") - if is_optimal - for i in 1:length(plants) - for j in 1:length(markets) - println(" $(plants[i]) $(markets[j]) = $(JuMP.value(ship[i,j]))") - end - end - else - println("The solver did not find an optimal solution.") - end -end +""" + example_cannery(; verbose = true) -function solveCannery(plants, markets, capacity, demand, distance, freight) - numplants = length(plants) - nummarkets = length(markets) - cannery = Model(with_optimizer(solver)) +JuMP implementation of the cannery problem from Dantzig, Linear Programming and +Extensions, Princeton University Press, Princeton, NJ, 1963. - @variable(cannery, ship[1:numplants, 1:nummarkets] >= 0) +Author: Louis Luangkesorn +Date: January 30, 2015 +""" +function example_cannery(; verbose = true) + plants = ["Seattle", "San-Diego"] + markets = ["New-York", "Chicago", "Topeka"] - # Ship no more than plant capacity - @constraint(cannery, capacity_con[i in 1:numplants], - sum(ship[i,j] for j in 1:nummarkets) <= capacity[i]) + # Capacity and demand in cases. + capacity = [350, 600] + demand = [300, 300, 300] - # Ship at least market demand - @constraint(cannery, demand_con[j in 1:nummarkets], - sum(ship[i,j] for i in 1:numplants) >= demand[j]) + # Distance in thousand miles. + distance = [2.5 1.7 1.8; 2.5 1.8 1.4] - # Minimize transporatation cost - @objective(cannery, Min, - sum(distance[i,j]* freight * ship[i,j] for i in 1:numplants, j in 1:nummarkets)) + # Cost per case per thousand miles. + freight = 90 - JuMP.optimize!(cannery) + num_plants = length(plants) + num_markets = length(markets) - term_status = JuMP.termination_status(cannery) - primal_status = JuMP.primal_status(cannery) - is_optimal = term_status == MOI.Optimal + cannery = Model(with_optimizer(GLPK.Optimizer)) - PrintSolution(is_optimal, plants, markets, ship) - return is_optimal -end + @variable(cannery, ship[1:num_plants, 1:num_markets] >= 0) + # Ship no more than plant capacity + @constraint(cannery, capacity_con[i in 1:num_plants], + sum(ship[i,j] for j in 1:num_markets) <= capacity[i] + ) -# PARAMETERS + # Ship at least market demand + @constraint(cannery, demand_con[j in 1:num_markets], + sum(ship[i,j] for i in 1:num_plants) >= demand[j] + ) -plants = ["Seattle", "San-Diego"] -markets = ["New-York", "Chicago", "Topeka"] + # Minimize transporatation cost + @objective(cannery, Min, sum(distance[i, j] * freight * ship[i, j] + for i in 1:num_plants, j in 1:num_markets) + ) -# capacity and demand in cases -capacitycases = [350, 600] -demandcases = [300, 300, 300] + JuMP.optimize!(cannery) -# distance in thousand miles -distanceKmiles = [2.5 1.7 1.8; - 2.5 1.8 1.4] + if verbose + println("RESULTS:") + for i in 1:num_plants + for j in 1:num_markets + println(" $(plants[i]) $(markets[j]) = $(JuMP.value(ship[i, j]))") + end + end + end -# cost per case per thousand miles -freightcost = 90 + @test JuMP.termination_status(cannery) == MOI.OPTIMAL + @test JuMP.primal_status(cannery) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(cannery) == 151200.0 +end -solveCannery(plants, markets, capacitycases, demandcases, distanceKmiles, freightcost) +example_cannery(verbose = false) diff --git a/examples/clnlbeam.jl b/examples/clnlbeam.jl new file mode 100644 index 00000000000..b0fb5b2e470 --- /dev/null +++ b/examples/clnlbeam.jl @@ -0,0 +1,55 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, Ipopt, Test +const MOI = JuMP.MathOptInterface + +""" +clnlbeam +Based on AMPL model +Copyright (C) 2001 Princeton University +All Rights Reserved + +Permission to use, copy, modify, and distribute this software and its +documentation for any purpose and without fee is hereby granted, provided that +the above copyright notice appear in all copies and that the copyright notice +and this permission notice appear in all supporting documentation. + +Source: H. Maurer and H.D. Mittelman, "The non-linear beam via optimal control +with bound state variables", Optimal Control Applications and Methods 12, pp. +19-31, 1991. +""" +function example_clnlbeam() + N = 1000 + h = 1/N + alpha = 350 + model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0)) + @variables(model, begin + -1 <= t[1:(N + 1)] <= 1 + -0.05 <= x[1:(N + 1)] <= 0.05 + u[1:(N + 1)] + end) + @NLobjective(model, Min, sum(0.5 * h * (u[i + 1]^2 + u[i]^2) + + 0.5 * alpha * h * (cos(t[i + 1]) + cos(t[i])) for i in 1:N) + ) + @NLconstraint(model, [i = 1:N], + x[i + 1] - x[i] - 0.5 * h * (sin(t[i + 1]) + sin(t[i])) == 0 + ) + @constraint(model, [i = 1:N], + t[i + 1] - t[i] - 0.5 * h * u[i + 1] - 0.5 * h * u[i] == 0 + ) + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 350.0 +end + +example_clnlbeam() diff --git a/examples/cluster.jl b/examples/cluster.jl index e60b9df64f3..17e01aa1b56 100644 --- a/examples/cluster.jl +++ b/examples/cluster.jl @@ -1,83 +1,82 @@ -############################################################################## +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# cluster.jl -# -# From "Approximating K-means-type clustering via semidefinite programming" -# By Jiming Peng and Yu Wei -# -# Given a set of points a_1, ..., a_m in R_n, allocate them to k clusters -############################################################################# -using JuMP -using SCS -using LinearAlgebra +using JuMP, SCS, LinearAlgebra, Test +""" + example_cluster() -solver = SCS.Optimizer +From "Approximating K-means-type clustering via semidefinite programming" By +Jiming Peng and Yu Wei. -# Data points -n = 2 -m = 6 -a = Any[[2.0, 2.0], [2.5, 2.1], [7.0, 7.0], - [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]] -k = 2 +Given a set of points a_1, ..., a_m in R_n, allocate them to k clusters. +""" +function example_cluster(; verbose = true) + # Data points + n = 2 + m = 6 + a = Any[[2.0, 2.0], [2.5, 2.1], [7.0, 7.0], + [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]] + k = 2 -# Weight matrix -W = zeros(m,m) -for i in 1:m - for j in i+1:m - dist = exp(-norm(a[i] - a[j])/1.0) - W[i,j] = dist - W[j,i] = dist + # Weight matrix + W = zeros(m, m) + for i in 1:m + for j in i + 1:m + W[i, j] = W[j, i] = exp(-norm(a[i] - a[j]) / 1.0) + end end -end - -mod = Model(with_optimizer(solver)) - -# Z >= 0, PSD -@variable(mod, Z[1:m,1:m], PSD) -@constraint(mod, Z .>= 0) - -# min Tr(W(I-Z)) -@objective(mod, Min, tr(W * (Matrix(1.0I,m,m) - Z))) -# Z e = e -@constraint(mod, Z*ones(m) .== ones(m)) + model = Model(with_optimizer(SCS.Optimizer, verbose = 0)) + # Z >= 0, PSD + @variable(model, Z[1:m, 1:m], PSD) + @constraint(model, Z .>= 0) + # min Tr(W(I-Z)) + @objective(model, Min, tr(W * (Matrix(1.0I, m, m) - Z))) + # Z e = e + @constraint(model, Z * ones(m) .== ones(m)) + # Tr(Z) = k + @constraint(model, tr(Z) == k) -# Tr(Z) = k -@constraint(mod, tr(Z) == k) + JuMP.optimize!(model) -JuMP.optimize!(mod) + Z_val = JuMP.value.(Z) -Z_val = JuMP.value.(Z) -println("Raw solution") -println(round.(Z_val, digits=4)) - -# A simple rounding scheme -which_cluster = zeros(Int,m) -num_clusters = 0 -for i in 1:m - Z_val[i,i] <= 1e-6 && continue - - if which_cluster[i] == 0 - global num_clusters += 1 - global which_cluster[i] = num_clusters - for j in i+1:m - if norm(Z_val[i,j] - Z_val[i,i]) <= 1e-6 - which_cluster[j] = num_clusters + # A simple rounding scheme + which_cluster = zeros(Int, m) + num_clusters = 0 + for i in 1:m + Z_val[i, i] <= 1e-6 && continue + if which_cluster[i] == 0 + num_clusters += 1 + which_cluster[i] = num_clusters + for j in i + 1:m + if norm(Z_val[i, j] - Z_val[i, i]) <= 1e-6 + which_cluster[j] = num_clusters + end end end end -end -# Print results -for cluster in 1:k - println("Cluster $cluster") - for i in 1:m - if which_cluster[i] == cluster - println(a[i]) + @test which_cluster == [1, 1, 2, 1, 2, 2] + + if verbose + # Print results + for cluster in 1:k + println("Cluster $cluster") + for i in 1:m + if which_cluster[i] == cluster + println(a[i]) + end + end end end end + +example_cluster(verbose = false) diff --git a/examples/corr_sdp.jl b/examples/corr_sdp.jl index 53d6548dbe4..df706ba4ffc 100644 --- a/examples/corr_sdp.jl +++ b/examples/corr_sdp.jl @@ -1,46 +1,53 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# corr_sdp.jl -# -# Given three random variables A,B,C and given bounds on two of the three -# correlation coefficients: -# -0.2 <= ρ_AB <= -0.1 -# 0.4 <= ρ_BC <= 0.5 -# We can use the following property of the correlations: -# [ 1 ρ_AB ρ_AC ] -# [ ρ_AB 1 ρ_BC ] ≽ 0 -# [ ρ_AC ρ_BC 1 ] -# To determine bounds on ρ_AC by solving a SDP -############################################################################# -using JuMP, CSDP +using JuMP, SCS, Test + +""" + example_corr_sdp() + +Given three random variables A,B,C and given bounds on two of the three +correlation coefficients: + -0.2 <= ρ_AB <= -0.1 + 0.4 <= ρ_BC <= 0.5 -m = Model(with_optimizer(CSDP.Optimizer)) +We can use the following property of the correlations to determine bounds on +ρ_AC by solving a SDP: + | 1 ρ_AB ρ_AC | + | ρ_AB 1 ρ_BC | ≽ 0 + | ρ_AC ρ_BC 1 | +""" +function example_corr_sdp() + model = Model(with_optimizer(SCS.Optimizer, verbose = 0)) + @variable(model, X[1:3, 1:3], PSD) -@variable(m, X[1:3,1:3], PSD) + # Diagonal is 1s + @constraint(model, X[1, 1] == 1) + @constraint(model, X[2, 2] == 1) + @constraint(model, X[3, 3] == 1) -# Diagonal is 1s -@constraint(m, X[1,1] == 1) -@constraint(m, X[2,2] == 1) -@constraint(m, X[3,3] == 1) + # Bounds on the known correlations + @constraint(model, X[1, 2] >= -0.2) + @constraint(model, X[1, 2] <= -0.1) + @constraint(model, X[2, 3] >= 0.4) + @constraint(model, X[2, 3] <= 0.5) -# Bounds on the known correlations -@constraint(m, X[1,2] >= -0.2) -@constraint(m, X[1,2] <= -0.1) -@constraint(m, X[2,3] >= 0.4) -@constraint(m, X[2,3] <= 0.5) + # Find upper bound + @objective(model, Max, X[1, 3]) + JuMP.optimize!(model) + @test JuMP.value(X[1, 3]) ≈ 0.87195 atol = 1e-4 -# Find upper bound -@objective(m, Max, X[1,3]) -JuMP.optimize!(m) -println("Maximum value is ", JuMP.value.(X)[1,3]) -@assert +0.8719 <= JuMP.value.(X)[1,3] <= +0.8720 + # Find lower bound + @objective(model, Min, X[1, 3]) + JuMP.optimize!(model) + @test JuMP.value(X[1, 3]) ≈ -0.978 atol = 1e-3 +end -# Find lower bound -@objective(m, Min, X[1,3]) -JuMP.optimize!(m) -println("Minimum value is ", JuMP.value.(X)[1,3]) -@assert -0.9779 >= JuMP.value.(X)[1,3] >= -0.9799 +example_corr_sdp() diff --git a/examples/sudoku.csv b/examples/data/sudoku.csv similarity index 100% rename from examples/sudoku.csv rename to examples/data/sudoku.csv diff --git a/examples/diet.jl b/examples/diet.jl index e35af475442..22c5f21e2be 100644 --- a/examples/diet.jl +++ b/examples/diet.jl @@ -4,90 +4,100 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# diet.jl -# -# Solve the classic "diet problem". -# Based on -# http://www.gurobi.com/documentation/5.6/example-tour/diet_cpp_cpp -############################################################################# -using JuMP, GLPK, LinearAlgebra +using JuMP, GLPK, Test const MOI = JuMP.MathOptInterface -solver = GLPK.Optimizer - -function PrintSolution(is_optimal, foods, buy) - println("RESULTS:") - if is_optimal - for i = 1:length(foods) - println(" $(foods[i]) = $(JuMP.value(buy[i]))") +""" + example_diet() + +Solve the classic "diet problem". Based on +http://www.gurobi.com/documentation/5.6/example-tour/diet_cpp_cpp +""" +function example_diet(; verbose = true) + function print_solution(is_optimal, foods, buy) + println("RESULTS:") + if is_optimal + for food in foods + println(" $(food) = $(JuMP.value(buy[food]))") + end + else + println("The solver did not find an optimal solution.") end - else - println("The solver did not find an optimal solution.") end -end - -function SolveDiet() # Nutrition guidelines - numCategories = 4 categories = ["calories", "protein", "fat", "sodium"] - minNutrition = [1800, 91, 0, 0] - maxNutrition = [2200, Inf, 65, 1779] + category_data = JuMP.Containers.DenseAxisArray([ + 1800 2200; + 91 Inf; + 0 65; + 0 1779 + ], categories, ["min", "max"] + ) + @test category_data["protein", "min"] == 91.0 + @test category_data["sodium", "max"] == 1779.0 # Foods - numFoods = 9 - foods = ["hamburger", "chicken", "hot dog", "fries", - "macaroni", "pizza", "salad", "milk", "ice cream"] - cost = [2.49, 2.89, 1.50, 1.89, 2.09, 1.99, 2.49, 0.89, 1.59] - nutritionValues = [410 24 26 730; - 420 32 10 1190; - 560 20 32 1800; - 380 4 19 270; - 320 12 10 930; - 320 15 12 820; - 320 31 12 1230; - 100 8 2.5 125; - 330 8 10 180] + foods = ["hamburger", "chicken", "hot dog", "fries", "macaroni", "pizza", + "salad", "milk", "ice cream"] + cost = JuMP.Containers.DenseAxisArray( + [2.49, 2.89, 1.50, 1.89, 2.09, 1.99, 2.49, 0.89, 1.59], + foods + ) + food_data = JuMP.Containers.DenseAxisArray( + [ + 410 24 26 730; + 420 32 10 1190; + 560 20 32 1800; + 380 4 19 270; + 320 12 10 930; + 320 15 12 820; + 320 31 12 1230; + 100 8 2.5 125; + 330 8 10 180 + ], foods, categories + ) + @test food_data["hamburger", "calories"] == 410.0 + @test food_data["milk", "fat"] == 2.5 # Build model - m = Model(with_optimizer(solver)) + model = Model(with_optimizer(GLPK.Optimizer)) - # Variables for nutrition info - @variable(m, minNutrition[i] <= nutrition[i in 1:numCategories] <= maxNutrition[i]) - # Variables for which foods to buy - @variable(m, buy[i in 1:numFoods] >= 0) + @variables(model, begin + # Variables for nutrition info + category_data[c, "min"] <= nutrition[c = categories] <= category_data[c, "max"] + # Variables for which foods to buy + buy[foods] >= 0 + end) # Objective - minimize cost - @objective(m, Min, dot(cost, buy)) + @objective(model, Min, sum(cost[f] * buy[f] for f in foods)) # Nutrition constraints - for j = 1:numCategories - @constraint(m, sum(nutritionValues[i,j]*buy[i] for i in 1:numFoods) == nutrition[j]) - end + @constraint(model, [c in categories], + sum(food_data[f, c] * buy[f] for f in foods) == nutrition[c] + ) # Solve - println("Solving original problem...") - JuMP.optimize!(m) - term_status = JuMP.termination_status(m) - primal_status = JuMP.primal_status(m) - is_optimal = term_status == MOI.Optimal - - PrintSolution(is_optimal, foods, buy) - - - # Limit dairy - @constraint(m, buy[8] + buy[9] <= 6) - println("Solving dairy-limited problem...") - JuMP.optimize!(m) - term_status = JuMP.termination_status(m) - primal_status = JuMP.primal_status(m) - is_optimal = term_status == MOI.Optimal - - PrintSolution(is_optimal, foods, buy) + verbose && println("Solving original problem...") + JuMP.optimize!(model) + term_status = JuMP.termination_status(model) + is_optimal = term_status == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 11.8288 atol = 1e-4 + verbose && print_solution(is_optimal, foods, buy) + + # Limit dairy (note that the problem will become infeasible). + @constraint(model, buy["milk"] + buy["ice cream"] <= 6) + verbose && println("Solving dairy-limited problem...") + JuMP.optimize!(model) + @test JuMP.termination_status(model) == MOI.INFEASIBLE_OR_UNBOUNDED + @test JuMP.primal_status(model) == MOI.NO_SOLUTION + verbose && print_solution(false, foods, buy) end -SolveDiet() +example_diet(verbose = false) diff --git a/examples/knapsack.jl b/examples/knapsack.jl index a2b30ab549d..9fc9e477caa 100644 --- a/examples/knapsack.jl +++ b/examples/knapsack.jl @@ -4,40 +4,44 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# knapsack.jl -# -# Solves a simple knapsack problem: -# max sum(p_j x_j) -# st sum(w_j x_j) <= C -# x binary -############################################################################# - -using JuMP, GLPK, LinearAlgebra - -# Maximization problem -m = Model(with_optimizer(GLPK.Optimizer)) - -@variable(m, x[1:5], Bin) - -profit = [ 5, 3, 2, 7, 4 ] -weight = [ 2, 8, 4, 2, 5 ] -capacity = 10 -# Objective: maximize profit -@objective(m, Max, dot(profit, x)) - -# Constraint: can carry all -@constraint(m, dot(weight, x) <= capacity) - -# Solve problem using MIP solver -JuMP.optimize!(m) - -println("Objective is: ", JuMP.objective_value(m)) -println("Solution is:") -for i = 1:5 - print("x[$i] = ", JuMP.value(x[i])) - println(", p[$i]/w[$i] = ", profit[i]/weight[i]) +using JuMP, GLPK, Test +const MOI = JuMP.MathOptInterface + +""" + example_knapsack(; verbose = true) + +Formulate and solve a simple knapsack problem: + max sum(p_j x_j) + st sum(w_j x_j) <= C + x binary +""" +function example_knapsack(; verbose = true) + profit = [5, 3, 2, 7, 4] + weight = [2, 8, 4, 2, 5] + capacity = 10 + model = Model(with_optimizer(GLPK.Optimizer)) + @variable(model, x[1:5], Bin) + # Objective: maximize profit + @objective(model, Max, profit' * x) + # Constraint: can carry all + @constraint(model, weight' * x <= capacity) + # Solve problem using MIP solver + JuMP.optimize!(model) + if verbose + println("Objective is: ", JuMP.objective_value(model)) + println("Solution is:") + for i in 1:5 + print("x[$i] = ", JuMP.value(x[i])) + println(", p[$i]/w[$i] = ", profit[i] / weight[i]) + end + end + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) == 16.0 end + +example_knapsack(verbose = false) diff --git a/examples/maxcut_sdp.jl b/examples/max_cut_sdp.jl similarity index 62% rename from examples/maxcut_sdp.jl rename to examples/max_cut_sdp.jl index 4d6682f4101..48389a5d60d 100644 --- a/examples/maxcut_sdp.jl +++ b/examples/max_cut_sdp.jl @@ -4,23 +4,15 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# maxcut_sdp.jl -# -# Demonstrates applying the SDP relaxation to the classic MAXCUT problem. -############################################################################# -using LinearAlgebra +using JuMP, SCS, Test, LinearAlgebra import Random -using JuMP -import SCS - - """ - solve_maxcut_sdp(num_vertex, weights) + solve_max_cut_sdp(num_vertex, weights) Solves a semidefinite programming relaxation of the MAXCUT graph problem: @@ -28,25 +20,23 @@ Solves a semidefinite programming relaxation of the MAXCUT graph problem: s.t. diag(X) == e X ≽ 0 -Where `L` is the weighted graph Laplacian. Uses this relaxation to generate -a solution to the original MAXCUT problem using the method from the paper: - - Goemans, M. X., & Williamson, D. P. (1995). - Improved approximation algorithms for maximum cut and satisfiability - problems using semidefinite programming. - Journal of the ACM (JACM), 42(6), 1115-1145. +Where `L` is the weighted graph Laplacian. Uses this relaxation to generate a +solution to the original MAXCUT problem using the method from the paper: +Goemans, M. X., & Williamson, D. P. (1995). Improved approximation algorithms +for maximum cut and satisfiability problems using semidefinite programming. +Journal of the ACM (JACM), 42(6), 1115-1145. """ -function solve_maxcut_sdp(num_vertex, weights) +function solve_max_cut_sdp(num_vertex, weights) # Calculate the (weighted) Lapacian of the graph: L = D - W. laplacian = diagm(0 => weights * ones(num_vertex)) - weights # Solve the SDP relaxation model = Model(with_optimizer(SCS.Optimizer)) @variable(model, X[1:num_vertex, 1:num_vertex], PSD) - @objective(model, Max, 1/4 * dot(laplacian, X)) + @objective(model, Max, 1 / 4 * dot(laplacian, X)) @constraint(model, diag(X) .== 1) JuMP.optimize!(model) - + # Compute the Cholesky factorization of X, i.e., X = V^T V. opt_X = Hermitian(JuMP.value.(X), :U) # Tell Julia its PSD. factorization = cholesky(opt_X, Val(true); check = false) @@ -61,7 +51,7 @@ function solve_maxcut_sdp(num_vertex, weights) Random.seed!(num_vertex) r = rand(num_vertex) r /= norm(r) - + # Iterate over vertices, and assign each vertex to a side of cut. cut = ones(num_vertex) for i in 1:num_vertex @@ -73,24 +63,14 @@ function solve_maxcut_sdp(num_vertex, weights) return cut, 0.25 * sum(laplacian .* (cut * cut')) end -function test1() +function example_max_cut_sdp() # [1] --- 5 --- [2] # # Solution: - # S = {1} - # S' = {2} - n = 2 - W = [0.0 5.0; - 5.0 0.0] - cut, cutval = solve_maxcut_sdp(n, W) + # (S, S′) = ({1}, {2}) + cut, cutval = solve_max_cut_sdp(2, [0.0 5.0; 5.0 0.0]) + @test cut[1] != cut[2] - @assert cut[1] != cut[2] - - println("Solution for Graph 1 = $cutval") - println(cut) -end - -function test2() # [1] --- 5 --- [2] # | \ | # | \ | @@ -100,23 +80,15 @@ function test2() # [3] --- 1 --- [4] # # Solution: - # S = {1} - # S' = {2,3,4} - n = 4 + # (S, S′) = ({1}, {2, 3, 4}) W = [0.0 5.0 7.0 6.0; 5.0 0.0 0.0 1.0; 7.0 0.0 0.0 1.0; 6.0 1.0 1.0 0.0] - cut, cutval = solve_maxcut_sdp(n, W) - @assert cut[2] != cut[1] - @assert cut[2] == cut[3] - @assert cut[2] == cut[4] + cut, cutval = solve_max_cut_sdp(4, W) + @test cut[1] != cut[2] + @test cut[2] == cut[3] == cut[4] - println("Solution for Graph 2 = $cutval") - println(cut) -end - -function test3() # [1] --- 1 --- [2] # | | # | | @@ -126,22 +98,15 @@ function test3() # [3] --- 2 --- [4] # # Solution: - # S = {2,3} - # S' = {1,4} - n = 4 + # (S, S′) = ({1, 4}, {2, 3}) W = [0.0 1.0 5.0 0.0; 1.0 0.0 0.0 9.0; 5.0 0.0 0.0 2.0; 0.0 9.0 2.0 0.0] - cut, cutval = solve_maxcut_sdp(n, W) - @assert cut[1] == cut[4] - @assert cut[2] == cut[3] - @assert cut[1] != cut[2] - - println("Solution for Graph 3 = $cutval") - println(cut) + cut, cutval = solve_max_cut_sdp(4, W) + @test cut[1] == cut[4] + @test cut[2] == cut[3] + @test cut[1] != cut[2] end -test1() -test2() -test3() +example_max_cut_sdp() diff --git a/examples/min_distortion.jl b/examples/min_distortion.jl new file mode 100644 index 00000000000..413b9db074e --- /dev/null +++ b/examples/min_distortion.jl @@ -0,0 +1,65 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, SCS, Test +const MOI = JuMP.MathOptInterface + +""" + example_min_distortion() + +This example arises from computational geometry, in particular the problem of +embedding a general finite metric space into a euclidean space. + +It is known that the 4-point metric space defined by the star graph: + x + \\ + x — x + / + x +where distances are computed by length of the shortest path between vertices, +cannot be exactly embedded into a euclidean space of any dimension. + +Here we will formulate and solve an SDP to compute the best possible embedding, +that is, the embedding f() that minimizes the distortion c such that + (1 / c) * D(a, b) ≤ ||f(a) - f(b)|| ≤ D(a, b) +for all points (a, b), where D(a, b) is the distance in the metric space. + +Any embedding can be characterized by its Gram matrix Q, which is PSD, and + ||f(a) - f(b)||^2 = Q[a, a] + Q[b, b] - 2 * Q[a, b] +We can therefore constrain + D[i, j]^2 ≤ Q[i, i] + Q[j, j] - 2 * Q[i, j] ≤ c^2 * D[i, j]^2 +and minimize c^2, which gives us the SDP formulation below. + +For more detail, see "Lectures on discrete geometry" by J. Matoušek, Springer, +2002, pp. 378-379. +""" +function example_min_distortion() + model = Model(with_optimizer(SCS.Optimizer, verbose = 0)) + D = [0.0 1.0 1.0 1.0; + 1.0 0.0 2.0 2.0; + 1.0 2.0 0.0 2.0; + 1.0 2.0 2.0 0.0] + @variable(model, c² >= 1.0) + @variable(model, Q[1:4, 1:4], PSD) + for i in 1:4 + for j in (i + 1):4 + @constraint(model, D[i, j]^2 <= Q[i, i] + Q[j, j] - 2 * Q[i, j]) + @constraint(model, Q[i, i] + Q[j, j] - 2 * Q[i, j] <= c² * D[i, j]^2) + end + end + @objective(model, Min, c²) + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 4/3 atol = 1e-4 +end + +example_min_distortion() diff --git a/examples/min_ellipse.jl b/examples/min_ellipse.jl new file mode 100644 index 00000000000..b52a09c242f --- /dev/null +++ b/examples/min_ellipse.jl @@ -0,0 +1,51 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, SCS, LinearAlgebra, Test +const MOI = JuMP.MathOptInterface + +""" + example_min_ellipse() + +This example is from the Boyd & Vandenberghe book "Convex Optimization". Given a +set of ellipses centered on the origin + E(A) = { u | u^T inv(A) u <= 1 } +find a "minimal" ellipse that contains the provided ellipses. + +We can formulate this as an SDP: + minimize trace(WX) + subject to X >= A_i, i = 1,...,m + X PSD +where W is a PD matrix of weights to choose between different solutions. +""" +function example_min_ellipse() + # We will use three ellipses: two "simple" ones, and a random one. + As = [ + [2.0 0.0; 0.0 1.0], + [1.0 0.0; 0.0 3.0], + [2.86715 1.60645; 1.60645 1.12639] + ] + # We change the weights to see different solutions, if they exist + weights = [1.0 0.0; 0.0 1.0] + model = Model(with_optimizer(SCS.Optimizer, verbose = 0)) + @variable(model, X[i=1:2, j=1:2], PSD, start = [2.0 0.0; 0.0 3.0][i, j]) + @objective(model, Min, tr(weights * X)) + for As_i in As + @SDconstraint(model, X >= As_i) + end + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 6.46233 atol = 1e-5 + @test JuMP.value.(X) ≈ [3.1651 0.8022; 0.8022 3.2972] atol = 1e-4 +end + +example_min_ellipse() diff --git a/examples/mindistortion.jl b/examples/mindistortion.jl deleted file mode 100644 index 19df154dcc1..00000000000 --- a/examples/mindistortion.jl +++ /dev/null @@ -1,60 +0,0 @@ -############################################################################# -# JuMP -# An algebraic modeling langauge for Julia -# See http://github.com/JuliaOpt/JuMP.jl -############################################################################# -# mindistortion.jl -# -# This example arises from computational geometry, in particular the problem -# of embedding a general finite metric space into a euclidean space. -# -# It is known that the 4-point metric space defined by the star graph: -# x -# \ -# x — x -# / -# x -# where distances are computed by length of the shortest path between -# vertices, cannot be exactly embedded into a euclidean space of any -# dimension. -# -# Here we will formulate and solve an SDP to compute the best possible -# embedding, that is, the embedding f() that minimizes the distortion c such -# that (1/c)*D(a,b) ≤ ||f(a)-f(b)|| ≤ D(a,b) for all points (a,b), where -# D(a,b) is the distance in the metric space. -# -# Any embedding can be characterized by its Gram matrix Q, which is PSD, -# and ||f(a)-f(b)||^2 = Q[a,a] + Q[b,b] - 2Q[a,b] -# We can therefore constrain -# D[i,j]^2 ≤ Q[i,i] + Q[j,j] - 2Q[i,j] ≤ c^2*D[i,j]^2 -# and minimize c^2, which gives us the SDP formulation below. -# -# For more detail, see -# "Lectures on discrete geometry" by J. Matoušek, Springer, 2002, pp. 378-379. - -using JuMP, CSDP - - -m = Model(with_optimizer(CSDP.Optimizer)) - -D = [0.0 1.0 1.0 1.0 - 1.0 0.0 2.0 2.0 - 1.0 2.0 0.0 2.0 - 1.0 2.0 2.0 0.0] - -@variable(m, cSq >= 1.0) - -@variable(m, Q[1:4,1:4], PSD) - -for i in 1:4 - for j in (i+1):4 - @constraint(m, D[i,j]^2 <= Q[i,i] + Q[j,j] - 2Q[i,j]) - @constraint(m, Q[i,i] + Q[j,j] - 2Q[i,j] <= D[i,j]^2*cSq ) - end -end - -@objective(m, Min, cSq) - -JuMP.optimize!(m) - -println(JuMP.value.(cSq)) diff --git a/examples/minellipse.jl b/examples/minellipse.jl deleted file mode 100644 index d0e780a9db2..00000000000 --- a/examples/minellipse.jl +++ /dev/null @@ -1,54 +0,0 @@ -############################################################################# -# JuMP -# An algebraic modeleling langauge for Julia -# See http://github.com/JuliaOpt/JuMP.jl -############################################################################# -# minellipse.jl -# -# This example is from the Boyd & Vandenberghe book "Convex Optimization". -# If PyPlot is installed, passing the argument `plot` will plot -# the solution, e.g. julia minellipse.jl plot -# -# Given a set of ellipses centered on the origin -# E(A) = { u | u^T inv(A) u <= 1 } -# find a "minimal" ellipse that contains the provided ellipses -# -# We can formulate this as an SDP: -# minimize trace(WX) -# subject to X >= A_i, i = 1,...,m -# X PSD -# where W is a PD matrix of weights to choose between different solutions -############################################################################# - -using JuMP -using SCS -using LinearAlgebra -#using PyPlot # Comment out if not installed - -# solver = SCSSolver(eps=1e-6) - -# We will use three ellipses -m = 3 -# Two "simple" ones -As = Any[ [2.0 0.0; - 0.0 1.0], - [1.0 0.0; - 0.0 3.0]] -# and a random one -randA = rand(2,2) -push!(As, (randA' * randA) * (rand()*2+1)) - -# We change the weights to see different solutions, if they exist -W = [1.0 0.0 - 0.0 1.0]; - -model = Model(with_optimizer(SCS.Optimizer)) -@variable(model, X[1:2, 1:2], PSD) -@objective(model, Min, tr(W*X)) -for i in 1:m - @SDconstraint(model, X >= As[i]) -end -JuMP.optimize!(model) - -X_val = JuMP.value.(X) -println(X_val) diff --git a/examples/mle.jl b/examples/mle.jl index 909689ee04e..f8d26ae83a5 100644 --- a/examples/mle.jl +++ b/examples/mle.jl @@ -2,36 +2,52 @@ # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -using JuMP, Ipopt, Statistics - -# Use nonlinear optimization to compute the maximum likelihood estimate (MLE) -# of the parameters of a normal distribution -# aka the sample mean and variance - -n = 1000 -data = randn(n) - -m = Model(with_optimizer(Ipopt.Optimizer, print_level=0)) - -@variable(m, μ, start = 0.0) -@variable(m, σ >= 0.0, start = 1.0) - -@NLobjective(m, Max, (n/2)*log(1/(2π*σ^2))-sum((data[i]-μ)^2 for i=1:n)/(2σ^2)) - -JuMP.optimize!(m) - -println("μ = ", JuMP.value(μ)) -println("mean(data) = ", mean(data)) -println("σ^2 = ", JuMP.value(σ)^2) -println("var(data) = ", var(data)) -println("MLE objective: ", JuMP.objective_value(m)) - -# constrained MLE? -@NLconstraint(m, μ == σ^2) - -JuMP.optimize!(m) -println("\nWith constraint μ == σ^2:") -println("μ = ", JuMP.value(μ)) -println("σ^2 = ", JuMP.value(σ)^2) - -println("Constrained MLE objective: ", JuMP.objective_value(m)) +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, Ipopt, Random, Statistics, Test + +""" + example_mle() + +Use nonlinear optimization to compute the maximum likelihood estimate (MLE) of +the parameters of a normal distribution aka the sample mean and variance +""" +function example_mle(; verbose = true) + n = 1_000 + Random.seed!(1234) + data = randn(n) + model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0)) + @variable(model, μ, start = 0.0) + @variable(model, σ >= 0.0, start = 1.0) + @NLobjective(model, Max, n / 2 * log(1 / (2 * π * σ^2)) - + sum((data[i] - μ)^2 for i in 1:n) / (2 * σ^2) + ) + JuMP.optimize!(model) + if verbose + println("μ = ", JuMP.value(μ)) + println("mean(data) = ", mean(data)) + println("σ^2 = ", JuMP.value(σ)^2) + println("var(data) = ", var(data)) + println("MLE objective: ", JuMP.objective_value(m)) + end + @test JuMP.value(μ) ≈ mean(data) atol = 1e-3 + @test JuMP.value(σ)^2 ≈ var(data) atol = 1e-2 + + # constrained MLE? + @NLconstraint(model, μ == σ^2) + + JuMP.optimize!(model) + if verbose + println("\nWith constraint μ == σ^2:") + println("μ = ", JuMP.value(μ)) + println("σ^2 = ", JuMP.value(σ)^2) + println("Constrained MLE objective: ", JuMP.objective_value(model)) + end + @test JuMP.value(μ) ≈ JuMP.value(σ)^2 +end + +example_mle(verbose = false) diff --git a/examples/multi.jl b/examples/multi.jl index 7de4fea9d03..6b6718eaf67 100644 --- a/examples/multi.jl +++ b/examples/multi.jl @@ -1,110 +1,103 @@ -############################################################################# # Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# Author: Louis Luangkesorn -# Date: February 26, 2015 -# -# JuMP implementation of the multicommodity transportation model -# AMPL: A Modeling Language for Mathematical Programming, 2nd ed -# by Robert Fourer, David Gay, and Brian W. Kernighan -# 4-1 multi.mod and multi.dat ############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, GLPK, Test +const MOI = JuMP.MathOptInterface -function PrintSolution(is_optimal, Trans, ORIG, DEST, PROD) - println("RESULTS:") - if is_optimal - for i in 1:length(ORIG) - for j in 1:length(DEST) - for p in 1:length(PROD) - print(" $(PROD[p]) $(ORIG[i]) $(DEST[j]) = $(JuMP.value(Trans[i,j, p]))\t") - end - println() +""" +Author: Louis Luangkesorn +Date: February 26, 2015 + +JuMP implementation of the multicommodity transportation model AMPL: A Modeling +Language for Mathematical Programming, 2nd ed by Robert Fourer, David Gay, and +Brian W. Kernighan 4-1. +""" +function example_multi(; verbose = true) + orig = ["GARY", "CLEV", "PITT"] + dest = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"] + prod = ["bands", "coils", "plate"] + numorig = length(orig) + numdest = length(dest) + numprod = length(prod) + + # supply(prod, orig) amounts available at origins + supply = [ + 400 700 800; + 800 1600 1800; + 200 300 300 + ] + + # demand(prod, dest) amounts required at destinations + demand = [ + 300 300 100 75 650 225 250; + 500 750 400 250 950 850 500; + 100 100 0 50 200 100 250 + ] + + # limit(orig, dest) of total units from any origin to destination + limit = [625.0 for j in 1:numorig, i in 1:numdest] + + # cost(dest, orig, prod) Shipment cost per unit + cost = reshape([[[ 30, 10, 8 , 10, 11 , 71, 6]; + [ 22, 7, 10, 7, 21, 82, 13]; + [ 19, 11, 12, 10, 25, 83, 15]]; + [[ 39, 14, 11, 14, 16, 82, 8]; + [ 27, 9, 12, 9, 26, 95, 17]; + [ 24, 14, 17, 13, 28, 99, 20]]; + [[ 41, 15, 12, 16, 17, 86, 8]; + [ 29, 9, 13, 9, 28, 99, 18]; + [ 26, 14, 17, 13, 31, 104, 20]]], + 7, 3, 3) + + # DECLARE MODEL + multi = Model(with_optimizer(GLPK.Optimizer)) + + # VARIABLES + @variable(multi, trans[1:numorig, 1:numdest, 1:numprod] >= 0) + + # OBJECTIVE + @objective(multi, Max, sum(cost[j, i, p] * trans[i, j, p] for + i in 1:numorig, j in 1:numdest, p in 1:numprod)) + + # CONSTRAINTS + + # Supply constraint + @constraint(multi, supply_con[i in 1:numorig, p in 1:numprod], + sum(trans[i, j, p] for j in 1:numdest) == supply[p, i]) + + # Demand constraint + @constraint(multi, demand_con[j in 1:numdest, p in 1:numprod], + sum(trans[i, j, p] for i in 1:numorig) == demand[p, j]) + + # Total shipment constraint + @constraint(multi, total_con[i in 1:numorig, j in 1:numdest], + sum(trans[i, j, p] for p in 1:numprod) - limit[i, j] <= 0) + + JuMP.optimize!(multi) + + @test JuMP.termination_status(multi) == MOI.OPTIMAL + @test JuMP.primal_status(multi) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(multi) == 225700.0 + + if verbose + println("RESULTS:") + for i in 1:length(ORIG) + for j in 1:length(DEST) + for p in 1:length(PROD) + print(" $(PROD[p]) $(ORIG[i]) $(DEST[j]) = $(JuMP.value(trans[i, j, p]))\t") + end + println() + end end - end - else - println(" No solution") end end -using JuMP, Clp -const MOI = JuMP.MathOptInterface - -solver = Clp.Optimizer - - -# PARAMETERS - -orig = ["GARY" "CLEV" "PITT"] -dest = ["FRA" "DET" "LAN" "WIN" "STL" "FRE" "LAF"] -prod = ["bands" "coils" "plate"] -numorig = length(orig) -numdest = length(dest) -numprod = length(prod) - -# supply(prod, orig) amounts available at origins -supply = [[ 400 700 800] - [ 800 1600 1800] - [ 200 300 300]] - -# demand(prod, dest) amounts required at destinations -demand = [[ 300 300 100 75 650 225 250] - [ 500 750 400 250 950 850 500] - [ 100 100 0 50 200 100 250]] - -# limit(orig, dest) of total units from any origin to destination -defaultlimit = float(625) - -limit = [[defaultlimit for j=1:numdest] for i=1:numorig] - -# cost(dest, orig, prod) Shipment cost per unit -cost = reshape([[[ 30, 10, 8 , 10, 11 , 71, 6]; - [ 22, 7, 10, 7, 21, 82, 13]; - [ 19, 11, 12, 10, 25, 83, 15]]; - [[ 39, 14, 11, 14, 16, 82, 8]; - [ 27, 9, 12, 9, 26, 95, 17]; - [ 24, 14, 17, 13, 28, 99, 20]]; - [[ 41, 15, 12, 16, 17, 86, 8]; - [ 29, 9, 13, 9, 28, 99, 18]; - [ 26, 14, 17, 13, 31, 104, 20]]], - 7, 3, 3) - -# DECLARE MODEL - -multi = Model(with_optimizer(solver)) - -# VARIABLES - -@variable(multi, Trans[1:numorig, 1:numdest, 1:numprod] >= 0) - -# OBJECTIVE - -length(cost) - -@objective(multi, Max, - sum(cost[j, i, p] * Trans[i,j, p] for - i in 1:numorig, j in 1:numdest, p in 1:numprod)) - -# CONSTRAINTS - -# Supply constraint -@constraint(multi, supply_con[i in 1:numorig, p in 1:numprod], - sum(Trans[i,j,p] for j in 1:numdest) == supply[p,i]) - -# Demand constraint -@constraint(multi, demand_con[j in 1:numdest, p in 1:numprod], - sum(Trans[i,j,p] for i in 1:numorig) == demand[p,j]) - -# Total shipment constraint -@constraint(multi, total_con[i in 1:numorig, j in 1:numdest], - sum(Trans[i,j,p] for p in 1:numprod) - limit[i][j] <= 0) - -JuMP.optimize!(multi) - -term_status = JuMP.termination_status(multi) -primal_status = JuMP.primal_status(multi) -is_optimal = term_status == MOI.Optimal -PrintSolution(is_optimal, Trans, orig, dest, prod) +example_multi(verbose = false) diff --git a/examples/optcontrol.jl b/examples/optcontrol.jl deleted file mode 100644 index 0a6b8f0cd0b..00000000000 --- a/examples/optcontrol.jl +++ /dev/null @@ -1,48 +0,0 @@ -# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. -using JuMP, Ipopt - -solver = Ipopt.Optimizer - -# clnlbeam -# Based on AMPL model -# Copyright (C) 2001 Princeton University -# All Rights Reserved -# -# Permission to use, copy, modify, and distribute this software and -# its documentation for any purpose and without fee is hereby -# granted, provided that the above copyright notice appear in all -# copies and that the copyright notice and this -# permission notice appear in all supporting documentation. - -# Source: -# H. Maurer and H.D. Mittelman, -# "The non-linear beam via optimal control with bound state variables", -# Optimal Control Applications and Methods 12, pp. 19-31, 1991. -let - N = 1000 - ni = N - h = 1/ni - alpha = 350 - - m = Model(with_optimizer(solver)) - - @variable(m, -1 <= t[1:(ni+1)] <= 1) - @variable(m, -0.05 <= x[1:(ni+1)] <= 0.05) - @variable(m, u[1:(ni+1)]) - - @NLobjective(m, Min, sum( 0.5*h*(u[i+1]^2 + u[i]^2) + 0.5*alpha*h*(cos(t[i+1]) + cos(t[i])) for i in 1:ni)) - - # cons1 - for i in 1:ni - @NLconstraint(m, x[i+1] - x[i] - (0.5h)*(sin(t[i+1])+sin(t[i])) == 0) - end - # cons2 - for i in 1:ni - @constraint(m, t[i+1] - t[i] - (0.5h)*u[i+1] - (0.5h)*u[i] == 0) - end - - JuMP.optimize!(m) -end diff --git a/examples/prod.jl b/examples/prod.jl index 31df1f44025..0d50ebc92ef 100644 --- a/examples/prod.jl +++ b/examples/prod.jl @@ -1,258 +1,257 @@ -############################################################################# # Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# Author: Louis Luangkesorn -# Date: February 26, 2015 -# -# This model determines a set of workforce levels that will most economically -# meet demands and inventory requirements over time. The formulation is -# motivated by the experiences of a large producer in the United States. The -# data are for three products and 13 numperiods -# -# Problem taken from the Appendix C of the expanded version of Fourer, Gay, -# and Kernighan, A Modeling Language for Mathematical Programming +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -function PrintSolution(is_optimal, CREWS, HIRE, LAYOFF) - println("RESULTS:") - if is_optimal - println("Crews") - for t = 0:length(CREWS.data)-1 - print(" $(JuMP.value(CREWS[t])) ") - end - println() - println("Hire") - for t = 1:length(HIRE.data) - print(" $(JuMP.value(HIRE[t])) ") - end - println() - println("Layoff") - for t = 1:length(LAYOFF.data) - print(" $(JuMP.value(LAYOFF[t])) ") - end - println() - else - println(" No solution") - end -end - - -using JuMP, Clp +using JuMP, GLPK, Test const MOI = JuMP.MathOptInterface -solver = Clp.Optimizer - -#### PRODUCTION SETS AND PARAMETERS ### - -prd = ["18REG" "24REG" "24PRO"] -# Members of the product group -numprd = length(prd) -pt= [1.194, 1.509, 1.509] -# Crew-hours to produce 1000 units -pc= [2304, 2920, 2910] -# Nominal production cost per 1000, used -# to compute inventory and shortage costs - - -### TIME PERIOD SETS AND PARAMETERS ### - -firstperiod = 1 -# Index of first production period to be modeled -lastperiod = 13 -# Index of last production period to be modeled -numperiods = firstperiod:lastperiod -# 'planning horizon' := first..last; - -### EMPLOYMENT PARAMETERS ### - -cs = 18 -# Workers per crew -sl = 8 -# Regular-time hours per shift -rtr = 16.00 -# Wage per hour for regular-time labor -otr = 43.85 -# Wage per hour for overtime labor -iw = 8 -# Crews employed at start of first period -dpp = [19.5, 19, 20, 19, 19.5, 19, 19, 20, 19, 20, 20, 18, 18] -# Regular working days in a production period -ol = [96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96] -# Maximum crew-hours of overtime in a period -cmin = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] -# Lower limit on average employment in a period -cmax = [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8] -# Upper limit on average employment in a period -hc = [7500, 7500, 7500, 7500, 15000, 15000, 15000, 15000, 15000, 15000, 7500, 7500, 7500] -# Penalty cost of hiring a crew -lc = [7500, 7500, 7500, 7500, 15000, 15000, 15000, 15000, 15000, 15000, 7500, 7500, 7500] -# Penalty cost of laying off a crew - -### DEMAND PARAMETERS ### - -d18REG = [63.8, 76, 88.4, 913.8, 115, 133.8, 79.6, 111, 121.6, 470, 78.4, 99.4, 140.4, 63.8] -d24REG = [1212, 306.2, 319, 208.4, 298, 328.2, 959.6, 257.6, 335.6, 118, 284.8, 970, 343.8, 1212] -d24PRO = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1102, 0, 0, 0, 0] -dem = Array[d18REG, d24REG, d24PRO] -# Requirements (in 1000s) -# to be met from current production and inventory - -pro = Array[ - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], - [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]] -# true if product will be the subject -# of a special promotion in the period - -### INVENTORY AND SHORTAGE PARAMETERS ### - -rir = 0.75 -# Proportion of non-promoted demand -# that must be in inventory the previous period -pir = 0.80 -# Proportion of promoted demand -# that must be in inventory the previous period -life = 2 -# Upper limit on number of periods that -# any product may sit in inventory -cri =[0.015, 0.015, 0.015] -# Inventory cost per 1000 units is -# cri times nominal production cost -crs =[1.1, 1.1, 1.1] -# Shortage cost per 1000 units is -# crs times nominal production cost -iinv =[82, 792.2, 0] -# Inventory at start of first period; age unknown -iil = [[max(0,iinv[p] - sum([dem[p][v] for v in firstperiod:t])) - for t in numperiods] - for p in 1:numprd] -# Initial inventory still available for allocation -# at end of period t -function checkpro(product, timeperiod, production, promotionalrate, regularrate) - if production[product][timeperiod + 1]==1 - value=promotionalrate - else - value =regularrate - end - value -end -minv = [[dem[p][t+1] * checkpro(p,t, pro, pir, rir) for t in numperiods] for p in 1:numprd] -# Lower limit on inventory at end of period t - -prod = Model(with_optimizer(solver)) - -### VARIABLES ### -@variable(prod, Crews[0:lastperiod] >= 0) -# Average number of crews employed in each period - -@variable(prod, Hire[numperiods] >= 0) -# Crews hired from previous to current period - -@variable(prod, Layoff[numperiods]>= 0) -# Crews laid off from previous to current period - -@variable(prod, Rprd[1:numprd, numperiods] >=0) -# Production using regular-time labor, in 1000s - -@variable(prod, Oprd[1:numprd, numperiods]>=0) -# Production using overtime labor, in 1000s - -@variable(prod, Inv[1:numprd, numperiods, 1:life] >=0) -# a numperiods old -- produced in period (t+1)-a -- -# and still in storage at the end of period t -@variable(prod, Short[1:numprd, numperiods]>=0) - # Accumulated unsatisfied demand at the end of period t - -### CONSTRAINTS ### - -@constraint(prod, [t=numperiods], - sum(pt[p] * Rprd[p,t] for p in 1:numprd) <= sl * dpp[t] * Crews[t]) -# Hours needed to accomplish all regular-time -# production in a period must not exceed -# hours available on all shifts - -@constraint(prod, [t=numperiods], - sum(pt[p] * Oprd[p,t] for p in 1:numprd) <= ol[t]) -# Hours needed to accomplish all overtime -# production in a period must not exceed -# the specified overtime limit - -@constraint(prod, Crews[firstperiod-1] == iw) -# Use given initial workforce - -@constraint(prod,[t in numperiods], - Crews[t] == Crews[t-1] + Hire[t] - Layoff[t]) -# Workforce changes by hiring or layoffs - -@constraint(prod, [t in numperiods], cmin[t] <= Crews[t]) -@constraint(prod, [t in numperiods], Crews[t] <= cmax[t]) -# Workforce must remain within specified bounds - -@constraint(prod, [p in 1:numprd], - Rprd[p, firstperiod] + Oprd[p,firstperiod] + Short[p, firstperiod] - Inv[p, firstperiod, 1] - == max(0, dem[p][firstperiod] - iinv[p])) -# 'first demand requirement - -# NOTE: JuMP xyconstr[] requires that indices be integer at compile time, -# so firstperiod +1 could not be an index within xycontr or triconstr -for t in (firstperiod+1:lastperiod) - @constraint(prod, [p in 1:numprd], - Rprd[p,t] + Oprd[p,t] + Short[p,t] - Short[p,t-1] + - sum(Inv[p, t-1, a] - Inv[p,t,a] for a in 1:life) == - max(0, dem[p][t] - iil[p][t-1])) +""" +Author: Louis Luangkesorn +Date: February 26, 2015 + +This model determines a set of workforce levels that will most economically meet +demands and inventory requirements over time. The formulation is motivated by +the experiences of a large producer in the United States. The data are for +three products and 13 numperiods. + +Problem taken from the Appendix C of the expanded version of Fourer, Gay, and +Kernighan, A Modeling Language for Mathematical Programming +""" +function example_prod(; verbose = true) + #### PRODUCTION SETS AND PARAMETERS ### + + prd = ["18REG" "24REG" "24PRO"] + # Members of the product group + numprd = length(prd) + pt= [1.194, 1.509, 1.509] + # Crew-hours to produce 1000 units + pc= [2304, 2920, 2910] + # Nominal production cost per 1000, used + # to compute inventory and shortage costs + + + ### TIME PERIOD SETS AND PARAMETERS ### + + firstperiod = 1 + # Index of first production period to be modeled + lastperiod = 13 + # Index of last production period to be modeled + numperiods = firstperiod:lastperiod + # 'planning horizon' := first..last; + + ### EMPLOYMENT PARAMETERS ### + + cs = 18 + # Workers per crew + sl = 8 + # Regular-time hours per shift + rtr = 16.00 + # Wage per hour for regular-time labor + otr = 43.85 + # Wage per hour for overtime labor + iw = 8 + # Crews employed at start of first period + dpp = [19.5, 19, 20, 19, 19.5, 19, 19, 20, 19, 20, 20, 18, 18] + # Regular working days in a production period + ol = [96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96] + # Maximum crew-hours of overtime in a period + cmin = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + # Lower limit on average employment in a period + cmax = [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8] + # Upper limit on average employment in a period + hc = [7500, 7500, 7500, 7500, 15000, 15000, 15000, 15000, 15000, 15000, 7500, 7500, 7500] + # Penalty cost of hiring a crew + lc = [7500, 7500, 7500, 7500, 15000, 15000, 15000, 15000, 15000, 15000, 7500, 7500, 7500] + # Penalty cost of laying off a crew + + ### DEMAND PARAMETERS ### + + d18REG = [63.8, 76, 88.4, 913.8, 115, 133.8, 79.6, 111, 121.6, 470, 78.4, 99.4, 140.4, 63.8] + d24REG = [1212, 306.2, 319, 208.4, 298, 328.2, 959.6, 257.6, 335.6, 118, 284.8, 970, 343.8, 1212] + d24PRO = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1102, 0, 0, 0, 0] + dem = Array[d18REG, d24REG, d24PRO] + # Requirements (in 1000s) + # to be met from current production and inventory + + pro = Array[ + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]] + # true if product will be the subject + # of a special promotion in the period + + ### INVENTORY AND SHORTAGE PARAMETERS ### + + rir = 0.75 + # Proportion of non-promoted demand + # that must be in inventory the previous period + pir = 0.80 + # Proportion of promoted demand + # that must be in inventory the previous period + life = 2 + # Upper limit on number of periods that + # any product may sit in inventory + cri =[0.015, 0.015, 0.015] + # Inventory cost per 1000 units is + # cri times nominal production cost + crs =[1.1, 1.1, 1.1] + # Shortage cost per 1000 units is + # crs times nominal production cost + iinv =[82, 792.2, 0] + # Inventory at start of first period; age unknown + iil = [[max(0,iinv[p] - sum([dem[p][v] for v in firstperiod:t])) + for t in numperiods] + for p in 1:numprd] + # Initial inventory still available for allocation + # at end of period t + function checkpro(product, timeperiod, production, promotionalrate, regularrate) + if production[product][timeperiod + 1]==1 + value=promotionalrate + else + value =regularrate + end + value + end + minv = [[dem[p][t+1] * checkpro(p,t, pro, pir, rir) for t in numperiods] for p in 1:numprd] + # Lower limit on inventory at end of period t + + prod = Model(with_optimizer(GLPK.Optimizer)) + + ### VARIABLES ### + @variable(prod, Crews[0:lastperiod] >= 0) + # Average number of crews employed in each period + + @variable(prod, Hire[numperiods] >= 0) + # Crews hired from previous to current period + + @variable(prod, Layoff[numperiods]>= 0) + # Crews laid off from previous to current period + + @variable(prod, Rprd[1:numprd, numperiods] >=0) + # Production using regular-time labor, in 1000s + + @variable(prod, Oprd[1:numprd, numperiods]>=0) + # Production using overtime labor, in 1000s + + @variable(prod, Inv[1:numprd, numperiods, 1:life] >=0) + # a numperiods old -- produced in period (t+1)-a -- + # and still in storage at the end of period t + @variable(prod, Short[1:numprd, numperiods]>=0) + # Accumulated unsatisfied demand at the end of period t + + ### CONSTRAINTS ### + + @constraint(prod, [t=numperiods], + sum(pt[p] * Rprd[p,t] for p in 1:numprd) <= sl * dpp[t] * Crews[t]) + # Hours needed to accomplish all regular-time + # production in a period must not exceed + # hours available on all shifts + + @constraint(prod, [t=numperiods], + sum(pt[p] * Oprd[p,t] for p in 1:numprd) <= ol[t]) + # Hours needed to accomplish all overtime + # production in a period must not exceed + # the specified overtime limit + + @constraint(prod, Crews[firstperiod-1] == iw) + # Use given initial workforce + + @constraint(prod,[t in numperiods], + Crews[t] == Crews[t-1] + Hire[t] - Layoff[t]) + # Workforce changes by hiring or layoffs + + @constraint(prod, [t in numperiods], cmin[t] <= Crews[t]) + @constraint(prod, [t in numperiods], Crews[t] <= cmax[t]) + # Workforce must remain within specified bounds + + @constraint(prod, [p in 1:numprd], + Rprd[p, firstperiod] + Oprd[p,firstperiod] + Short[p, firstperiod] - Inv[p, firstperiod, 1] + == max(0, dem[p][firstperiod] - iinv[p])) + # 'first demand requirement + + # NOTE: JuMP xyconstr[] requires that indices be integer at compile time, + # so firstperiod +1 could not be an index within xycontr or triconstr + for t in (firstperiod+1:lastperiod) + @constraint(prod, [p in 1:numprd], + Rprd[p,t] + Oprd[p,t] + Short[p,t] - Short[p,t-1] + + sum(Inv[p, t-1, a] - Inv[p,t,a] for a in 1:life) == + max(0, dem[p][t] - iil[p][t-1])) + end + # Production plus increase in shortage plus + # decrease in inventory must equal demand + + @constraint(prod, [p in 1:numprd, t in numperiods], + sum(Inv[p,t,a] + iil[p][t] for a in 1:life) >= minv[p][t]) + # Inventory in storage at end of period t + # must meet specified minimum + + @constraint(prod, [p in 1:numprd, v in 1:(life-1), a in v+1:life], Inv[p, firstperiod+v-1, a] ==0) + + # In the vth period (starting from first) + # no inventory may be more than v numperiods old + # (initial inventories are handled separately) + + @constraint(prod, [p in 1:numprd, t in numperiods], + Inv[p,t,1] <= Rprd[p,t] + Oprd[p,t]) + # New inventory cannot exceed + # production in the most recent period + + secondperiod = firstperiod + 1 + @constraint(prod, [p in 1:numprd, t in 2:lastperiod, a in 2:life], + Inv[p,t,a] <= Inv[p,t-1,a-1]) + # Inventory left from period (t+1)-p + # can only decrease as time goes on + + ### OBJECTIVE ### + + @objective(prod, Min, + sum(rtr * sl * dpp[t] * cs * Crews[t] for t in numperiods) + + sum(hc[t] * Hire[t] for t in numperiods) + + sum(lc[t] * Layoff[t] for t in numperiods) + + sum(otr * cs * pt[p] * Oprd[p,t] for t in numperiods, p in 1:numprd) + + sum(cri[p] * pc[p] * Inv[p,t,a] for t in numperiods, p in 1:numprd, a in 1:life) + + sum(crs[p] * pc[p] * Short[p,t] for t in numperiods, p in 1:numprd)) + + # Full regular wages for all crews employed, plus + # penalties for hiring and layoffs, plus + # wages for any overtime worked, plus + # inventory and shortage costs + + # (All other production costs are assumed + # to depend on initial inventory and on demands, + # and so are not included explicitly.) + + + JuMP.optimize!(prod) + + @test JuMP.termination_status(prod) == MOI.OPTIMAL + @test JuMP.primal_status(prod) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(prod) ≈ 4426822.89 atol = 1e-2 + + if verbose + println("RESULTS:") + println("Crews") + for t = 0:length(Crews.data) - 1 + print(" $(JuMP.value(Crews[t])) ") + end + println() + println("Hire") + for t = 1:length(Hire.data) + print(" $(JuMP.value(Hire[t])) ") + end + println() + println("Layoff") + for t = 1:length(Layoff.data) + print(" $(JuMP.value(Layoff[t])) ") + end + println() + end end -# Production plus increase in shortage plus -# decrease in inventory must equal demand - -@constraint(prod, [p in 1:numprd, t in numperiods], - sum(Inv[p,t,a] + iil[p][t] for a in 1:life) >= minv[p][t]) -# Inventory in storage at end of period t -# must meet specified minimum - -@constraint(prod, [p in 1:numprd, v in 1:(life-1), a in v+1:life], Inv[p, firstperiod+v-1, a] ==0) - -# In the vth period (starting from first) -# no inventory may be more than v numperiods old -# (initial inventories are handled separately) - -@constraint(prod, [p in 1:numprd, t in numperiods], - Inv[p,t,1] <= Rprd[p,t] + Oprd[p,t]) -# New inventory cannot exceed -# production in the most recent period - -secondperiod = firstperiod + 1 -@constraint(prod, [p in 1:numprd, t in 2:lastperiod, a in 2:life], - Inv[p,t,a] <= Inv[p,t-1,a-1]) -# Inventory left from period (t+1)-p -# can only decrease as time goes on - -### OBJECTIVE ### - -@objective(prod, Min, - sum(rtr * sl * dpp[t] * cs * Crews[t] for t in numperiods) + - sum(hc[t] * Hire[t] for t in numperiods) + - sum(lc[t] * Layoff[t] for t in numperiods) + - sum(otr * cs * pt[p] * Oprd[p,t] for t in numperiods, p in 1:numprd) + - sum(cri[p] * pc[p] * Inv[p,t,a] for t in numperiods, p in 1:numprd, a in 1:life) + - sum(crs[p] * pc[p] * Short[p,t] for t in numperiods, p in 1:numprd)) - -# Full regular wages for all crews employed, plus -# penalties for hiring and layoffs, plus -# wages for any overtime worked, plus -# inventory and shortage costs - -# (All other production costs are assumed -# to depend on initial inventory and on demands, -# and so are not included explicitly.) - - -JuMP.optimize!(prod) - -term_status = JuMP.termination_status(prod) -primal_status = JuMP.primal_status(prod) -is_optimal = status == MOI.Optimal -PrintSolution(is_optimal, Crews, Hire, Layoff) +example_prod(verbose = false) diff --git a/examples/qcp.jl b/examples/qcp.jl index 35e12305262..341393534ac 100644 --- a/examples/qcp.jl +++ b/examples/qcp.jl @@ -4,40 +4,40 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# qcp.jl -# -# A simple quadratically constrained probgram -# Based on http://www.gurobi.com/documentation/5.5/example-tour/node25 -############################################################################# - -using JuMP, Gurobi - -# Will require either Gurobi.jl, CPLEX.jl, or Mosek.jl to run -m = Model(with_optimizer(Gurobi.Optimizer)) - -# Need nonnegativity for (rotated) second-order cone -@variable(m, x) -@variable(m, y >= 0) -@variable(m, z >= 0) - -# Maximize x -@objective(m, Max, x) - -# Subject to 1 linear and 2 nonlinear constraints -@constraint(m, x + y + z == 1) -@constraint(m, x*x + y*y - z*z <= 0) -@constraint(m, x*x - y*z <= 0) - -# Print the model to check correctness -print(m) - -# Solve with Gurobi -JuMP.optimize!(m) -# Solution -println("Objective value: ", JuMP.objective_value(m)) -println("x = ", JuMP.value(x)) -println("y = ", JuMP.value(y)) +using JuMP, Ipopt, Test +const MOI = JuMP.MathOptInterface + +""" + example_qcp(; verbose = true) + +A simple quadratically constrained program based on +http://www.gurobi.com/documentation/5.5/example-tour/node25 +""" +function example_qcp(; verbose = true) + model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0)) + @variable(model, x) + @variable(model, y >= 0) + @variable(model, z >= 0) + @objective(model, Max, x) + @constraint(model, x + y + z == 1) + @constraint(model, x * x + y * y - z * z <= 0) + @constraint(model, x * x - y * z <= 0) + JuMP.optimize!(model) + if verbose + print(model) + println("Objective value: ", JuMP.objective_value(model)) + println("x = ", JuMP.value(x)) + println("y = ", JuMP.value(y)) + end + @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 0.32699 atol = 1e-5 + @test JuMP.value(x) ≈ 0.32699 atol = 1e-5 + @test JuMP.value(y) ≈ 0.25707 atol = 1e-5 +end + +example_qcp(verbose = false) diff --git a/examples/robust_uncertainty.jl b/examples/robust_uncertainty.jl index 196a17ac824..84eb8fe58c8 100644 --- a/examples/robust_uncertainty.jl +++ b/examples/robust_uncertainty.jl @@ -1,51 +1,53 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# robust_uncertainty.jl -# -# Computes the Value at Risk for a data-driven uncertainty set; see -# "Data-Driven Robust Optimization" (Bertsimas 2013), section 6.1 for -# details. Closed-form expressions for the optimal value are available. -############################################################################# - - -using JuMP, SCS, LinearAlgebra -R = 1 -d = 3 -𝛿 = 0.05 -ɛ = 0.05 -N = ceil((2+2log(2/𝛿))^2) + 1 +using JuMP, SCS, LinearAlgebra, Test -Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿))) -Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿))) +""" + example_robust_uncertainty() -μhat = rand(d) -M = rand(d,d) -Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat') +Computes the Value at Risk for a data-driven uncertainty set; see "Data-Driven +Robust Optimization" (Bertsimas 2013), section 6.1 for details. Closed-form +expressions for the optimal value are available. +""" +function example_robust_uncertainty() + R = 1 + d = 3 + 𝛿 = 0.05 + ɛ = 0.05 + N = ceil((2 + 2 * log(2 / 𝛿))^2) + 1 -m = Model(with_optimizer(SCS.Optimizer)) + c = randn(d) -@variable(m, Σ[1:d, 1:d], PSD) -@variable(m, u[1:d]) -@variable(m, μ[1:d]) + μhat = rand(d) + M = rand(d, d) + Σhat = 1 / (d - 1) * (M - ones(d) * μhat')' * (M - ones(d) * μhat') -@constraint(m, [Γ1(𝛿/2,N); μ-μhat] in SecondOrderCone()) -@constraint(m, [Γ2(𝛿/2,N); vec(Σ-Σhat)] in SecondOrderCone()) + Γ1(𝛿, N) = R / sqrt(N) * (2 + sqrt(2 * log(1 / 𝛿))) + Γ2(𝛿, N) = 2 * R^2 / sqrt(N) * (2 + sqrt(2 * log(2 / 𝛿))) -A = [(1-ɛ)/ɛ (u-μ)'; - (u-μ) Σ ] -@SDconstraint(m, A >= 0) + model = Model(with_optimizer(SCS.Optimizer, verbose = 0)) -c = randn(d) -@objective(m, Max, dot(c,u)) + @variable(model, Σ[1:d, 1:d], PSD) + @variable(model, u[1:d]) + @variable(model, μ[1:d]) + @constraint(model, [Γ1(𝛿 / 2, N); μ - μhat] in SecondOrderCone()) + @constraint(model, [Γ2(𝛿 / 2, N); vec(Σ - Σhat)] in SecondOrderCone()) + @SDconstraint(model, [((1 - ɛ) / ɛ) (u - μ)'; (u - μ) Σ] >= 0) + @objective(model, Max, dot(c, u)) -JuMP.optimize!(m) + JuMP.optimize!(model) -object = JuMP.objective_value(m) -exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*Matrix(1.0I,d,d))*c)) + exact = dot(μhat, c) + Γ1(𝛿 / 2, N) * norm(c) + sqrt((1 - ɛ) / ɛ) * + sqrt(dot(c, (Σhat + Γ2(𝛿 / 2, N) * Matrix(1.0I, d, d)) * c)) + @test JuMP.objective_value(model) ≈ exact atol = 1e-3 +end -println("objective value: $(object)") -println("error from exact: $(abs(exact-object))") +example_robust_uncertainty() diff --git a/examples/rosenbrock.jl b/examples/rosenbrock.jl index 1a6f7f8e3d1..91a3a70b659 100644 --- a/examples/rosenbrock.jl +++ b/examples/rosenbrock.jl @@ -2,20 +2,27 @@ # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -using JuMP, Ipopt -# rosenbrock +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# -let +using JuMP, Ipopt, Test +const MOI = JuMP.MathOptInterface - m = Model(with_optimizer(Ipopt.Optimizer, print_level=0)) - - @variable(m, x) - @variable(m, y) - - @NLobjective(m, Min, (1-x)^2 + 100(y-x^2)^2) - - JuMP.optimize!(m) - - println("x = ", JuMP.value(x), " y = ", JuMP.value(y)) +function example_rosenbrock() + model = Model(with_optimizer(Ipopt.Optimizer, print_level=0)) + @variable(model, x) + @variable(model, y) + @NLobjective(model, Min, (1 - x)^2 + 100 * (y - x^2)^2) + JuMP.optimize!(model) + @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 0.0 atol = 1e-10 + @test JuMP.value(x) ≈ 1.0 + @test JuMP.value(y) ≈ 1.0 end + +example_rosenbrock() diff --git a/examples/run_examples.jl b/examples/run_examples.jl new file mode 100644 index 00000000000..31195ea715e --- /dev/null +++ b/examples/run_examples.jl @@ -0,0 +1,23 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +# This script runs all the JuMP examples. +# Usage: julia --project=examples run_examples.jl + +using Test + +const EXAMPLES = filter(ex -> endswith(ex, ".jl") && ex != "run_examples.jl", + readdir(@__DIR__)) + +@testset "run_examples.jl" begin + @testset "$(example)" for example in EXAMPLES + include(example) + end +end diff --git a/examples/steelT3.jl b/examples/steelT3.jl index eb61c418226..9630e51ffad 100644 --- a/examples/steelT3.jl +++ b/examples/steelT3.jl @@ -1,128 +1,118 @@ -################################################################################ # Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# steelT3 model from -# AMPL: A Modeling Language for Mathematical Programming, 2nd ed -# by Robert Fourer, David Gay, and Brian W. Kernighan -# -# Author : Louis Luangkesorn -# Date: April 3, 2015 -# -################################################################################ +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# -using JuMP, Clp +using JuMP, GLPK, Test const MOI = JuMP.MathOptInterface -solver = Clp.Optimizer -# Data - -T = 4 -prod = ["bands", "coils"] -area = Dict([("bands", ("east", "north")), - ("coils", ("east", "west", "export"))]) - -avail = [40, 40, 32,40] - -rate = Dict("bands" => 200, - "coils" =>140) - -inv0 = Dict("bands" => 10, - "coils" => 0) - -prodcost = Dict("bands" => 10, - "coils" => 11) -invcost = Dict("bands" => 2.5, - "coils" => 3) - -revenue = Dict("bands" => Dict("east" => [25.0 26.0 27.0 27.0], - "north" => [26.5 27.5 28.0 28.5]), - "coils" => Dict("east" =>[30 35 37 39], - "west" => [29 32 33 35], - "export" => [25 25 25 28]) - ) - -market = Dict("bands" => Dict("east" => [2000 2000 1500 2000], - "north" => [4000 4000 2500 4500]), - "coils" => Dict("east" => [1000 800 1000 1100], - "west" => [2000 1200 2000 2300], - "export" => [1000 500 500 800]) - ) - -# Model - -Prod = Model(with_optimizer(solver)) -print(area["bands"]) -print(market["coils"]["west"][1]) - -# Decision Variables - -@variable(Prod, Make[p in prod, t in 1:T] >= 0) # tons produced -@variable(Prod, Inv[p in prod, t in 0:T] >= 0) # tons inventoried -@variable(Prod, market[p][a][t] >= Sell[p in prod, a in area[p], t in 1:T] >= 0) # tons sold - -@constraint(Prod, [p in prod, a in area[p], t in 1:T], - Sell[p, a, t] - market[p][a][t] <= 0) - - -@constraint(Prod, [t in 1:T], - sum((1/rate[p]) * Make[p,t] for p in prod) <= avail[t]) - -# Total of hours used by all products -# may not exceed hours available, in each week - -@constraint(Prod, [p in prod], - Inv[p,0] == inv0[p]) -# Initial inventory must equal given value - -@constraint(Prod, [p in prod, t in 1:T], - Make[p,t] + Inv[p, t-1] == sum(Sell[p,a,t] for a in area[p]) + Inv[p,t]) - -# Tons produced and taken from inventory -# must equal tons sold and put into inventory - - -@objective(Prod, Max, - sum( sum( - revenue[p][a][t] * Sell[p, a, t] - - prodcost[p] * Make[p,t] - - invcost[p]*Inv[p,t] for a in area[p]) for p in prod, t in 1:T)) -#maximize Total_Profit: -# Total revenue less costs for all products in all weeks - -function PrintSolution(is_optimal, area, Make, Inventory, Sell, product, Time) - println("RESULTS:") - if is_optimal - for p in product - println("Make $(p)") - for t in 1:T - print("$(JuMP.value(Make[p,t]))\t") - end - println() - println("Inventory $(p)") - for t in 1:T - print("$(JuMP.value(Inventory[p,t]))\t") +""" + example_steelT3(; verbose = true) + +steelT3 model from AMPL: A Modeling Language for Mathematical Programming, 2nd +ed by Robert Fourer, David Gay, and Brian W. Kernighan. + +Author: Louis Luangkesorn +Date: April 3, 2015 +""" +function example_steelT3(; verbose = true) + T = 4 + prod = ["bands", "coils"] + area = Dict( + "bands" => ("east", "north"), + "coils" => ("east", "west", "export") + ) + avail = [40, 40, 32, 40] + rate = Dict("bands" => 200, "coils" => 140) + inv0 = Dict("bands" => 10, "coils" => 0) + prodcost = Dict("bands" => 10, "coils" => 11) + invcost = Dict("bands" => 2.5, "coils" => 3) + revenue = Dict( + "bands" => Dict("east" => [25.0, 26.0, 27.0, 27.0], + "north" => [26.5, 27.5, 28.0, 28.5]), + "coils" => Dict("east" =>[30, 35, 37, 39], + "west" => [29, 32, 33, 35], + "export" => [25, 25, 25, 28]) + ) + market = Dict( + "bands" => Dict( + "east" => [2000, 2000, 1500, 2000], + "north" => [4000, 4000, 2500, 4500] + ), + "coils" => Dict( + "east" => [1000, 800, 1000, 1100], + "west" => [2000, 1200, 2000, 2300], + "export" => [1000, 500, 500, 800] + ) + ) + + # Model + model = Model(with_optimizer(GLPK.Optimizer)) + + # Decision Variables + @variables(model, begin + make[p in prod, t in 1:T] >= 0 + inventory[p in prod, t in 0:T] >= 0 + 0 <= sell[p in prod, a in area[p], t in 1:T] <= market[p][a][t] + end) + + @constraint(model, [p in prod, a in area[p], t in 1:T], + sell[p, a, t] - market[p][a][t] <= 0) + + # Total of hours used by all products may not exceed hours available, in + # each week + @constraint(model, [t in 1:T], + sum(1 / rate[p] * make[p, t] for p in prod) <= avail[t]) + + # Initial inventory must equal given value + @constraint(model, [p in prod], inventory[p, 0] == inv0[p]) + + # Tons produced and taken from inventory must equal tons sold and put into + # inventory. + @constraint(model, [p in prod, t in 1:T], + make[p, t] + inventory[p, t - 1] == sum(sell[p, a, t] for a in area[p]) + inventory[p, t]) + + # Maximize total profit: total revenue less costs for all products in all + # weeks. + @objective(model, Max, sum( + sum(revenue[p][a][t] * sell[p, a, t] - prodcost[p] * make[p, t] - + invcost[p] * inventory[p, t] for a in area[p]) + for p in prod, t in 1:T) + ) + + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) == 172850.0 + + if verbose + println("RESULTS:") + for p in product + println("make $(p)") + for t in 1:T + print(JuMP.value(make[p, t]), "\t") + end + println() + println("Inventory $(p)") + for t in 1:T + print(JuMP.value(Inventory[p, t]), "\t") + end + println() + for a in area[p] + println("sell $(p) $(a)") + for t in 1:T + print(JuMP.value(sell[p, a, t]), "\t") + end + println() + end end - println() - for a in area[p] - println("Sell $(p) $(a)") - for t in 1:T - print("$(JuMP.value(Sell[p,a,t])) \t") - end - println() - end - end - else - println(" No solution") end end -JuMP.optimize!(Prod) - -term_status = JuMP.termination_status(Prod) -primal_status = JuMP.primal_status(Prod) -is_optimal = status == MOI.Optimal - -PrintSolution(is_optimal, area, Make, Inv, Sell, prod, T) +example_steelT3(verbose = false) diff --git a/examples/sudoku.jl b/examples/sudoku.jl index 99a4dea755e..134ea3214fd 100644 --- a/examples/sudoku.jl +++ b/examples/sudoku.jl @@ -4,75 +4,70 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. ############################################################################# # JuMP -# An algebraic modeling langauge for Julia +# An algebraic modeling language for Julia # See http://github.com/JuliaOpt/JuMP.jl ############################################################################# -# sudoku.jl -# A sudoku solver that uses a MIP to find solutions -# -# We have binary variables x[i,j,k] which, if = 1, say that cell (i,j) -# contains the number k -# The constraints are: -# 1 - Each cell has one value only -# 2 - Each row contains each number exactly once -# 3 - Each column contains each number exactly once -# 4 - Each 3x3 subgrid contains each number exactly once -# We will take the initial grid as a CSV file, where 0s are "blanks -############################################################################# -using JuMP, GLPK +using JuMP, GLPK, Test const MOI = JuMP.MathOptInterface -solver = GLPK.Optimizer +""" + example_sudoku(filepath) -# Load data -function LoadData(filepath) - fp = open(filepath, "r") - initgrid = zeros(Int,9,9) - for row in 1:9 - line = readline(fp) - initgrid[row,:] = [parse(Int,s) for s in split(line,",")] +A sudoku solver that uses a MIP to find solutions. + +We have binary variables x[i, j, k] which, if = 1, say that cell (i, j) contains +the number k. The constraints are: + 1 - Each cell has one value only + 2 - Each row contains each number exactly once + 3 - Each column contains each number exactly once + 4 - Each 3x3 subgrid contains each number exactly once +We will take the initial grid as a CSV file at `filepath`, where 0s are blanks. +""" +function example_sudoku(filepath) + initial_grid = zeros(Int, 9, 9) + open(filepath, "r") do fp + for row in 1:9 + line = readline(fp) + initial_grid[row, :] .= parse.(Int, split(line, ",")) + end end - return initgrid -end -# Solve model -function SolveModel(initgrid) - m = Model(with_optimizer(solver)) + model = Model(with_optimizer(GLPK.Optimizer)) - @variable(m, x[1:9, 1:9, 1:9], Bin) + @variable(model, x[1:9, 1:9, 1:9], Bin) - @constraints m begin + @constraints(model, begin # Constraint 1 - Only one value appears in each cell + cell[i in 1:9, j in 1:9], sum(x[i, j, :]) == 1 # Constraint 2 - Each value appears in each row once only + row[i in 1:9, k in 1:9], sum(x[i, :, k]) == 1 # Constraint 3 - Each value appears in each column once only - cell[i in 1:9, j in 1:9], sum(x[i,j,:]) == 1 - row[i in 1:9, k in 1:9], sum(x[i,:,k]) == 1 - col[j in 1:9, k in 1:9], sum(x[:,j,k]) == 1 + col[j in 1:9, k in 1:9], sum(x[:, j, k]) == 1 # Constraint 4 - Each value appears in each 3x3 subgrid once only - subgrid[i=1:3:7,j=1:3:7,val=1:9], sum(x[i:i+2,j:j+2,val]) == 1 - end + subgrid[i=1:3:7, j=1:3:7, val=1:9], sum(x[i:i + 2, j:j + 2, val]) == 1 + end) # Initial solution for row in 1:9, col in 1:9 - if initgrid[row,col] != 0 - @constraint(m, x[row, col, initgrid[row, col]] == 1) + if initial_grid[row, col] != 0 + @constraint(model, x[row, col, initial_grid[row, col]] == 1) end end # Solve it - JuMP.optimize!(m) + JuMP.optimize!(model) - term_status = JuMP.termination_status(m) - primal_status = JuMP.primal_status(m) - is_optimal = term_status == MOI.Optimal + term_status = JuMP.termination_status(model) + primal_status = JuMP.primal_status(model) + is_optimal = term_status == MOI.OPTIMAL # Check solution if is_optimal - mipSol = JuMP.value.(x) - sol = zeros(Int,9,9) + mip_solution = JuMP.value.(x) + sol = zeros(Int, 9, 9) for row in 1:9, col in 1:9, val in 1:9 - if mipSol[row, col, val] >= 0.9 + if mip_solution[row, col, val] >= 0.9 sol[row, col] = val end end @@ -80,30 +75,39 @@ function SolveModel(initgrid) else error("The solver did not find an optimal solution.") end - -end - -# Initialization -if length(ARGS) != 1 - error("Expected one argument: the initial solution, e.g. julia sudoku.jl sudoku.csv") end -# Solve -sol = SolveModel(LoadData(ARGS[1])) - -# Display solution -println("Solution:") -println("[-----------------------]") -for row in 1:9 - print("[ ") - for col in 1:9 - print("$(sol[row,col]) ") - if col % 3 == 0 && col < 9 - print("| ") +function print_sudoku_solution(solution) + println("Solution:") + println("[-----------------------]") + for row in 1:9 + print("[ ") + for col in 1:9 + print(solution[row, col], " ") + if col % 3 == 0 && col < 9 + print("| ") + end + end + println("]") + if row % 3 == 0 + println("[-----------------------]") end end - println("]") - if row % 3 == 0 - println("[-----------------------]") - end +end + +if length(ARGS) == 1 + sol = example_sudoku(ARGS[1]) + print_sudoku_solution(sol) +else + solution = example_sudoku(joinpath(@__DIR__, "data/sudoku.csv")) + @test solution == [ + 3 1 7 9 5 8 2 6 4; + 4 6 9 3 2 7 8 1 5; + 8 2 5 1 6 4 7 9 3; + 2 7 1 6 4 5 3 8 9; + 9 3 8 7 1 2 5 4 6; + 5 4 6 8 3 9 1 7 2; + 1 8 4 2 9 3 6 5 7; + 6 9 2 5 7 1 4 3 8; + 7 5 3 4 8 6 9 2 1] end diff --git a/examples/sudokuhard.csv b/examples/sudokuhard.csv deleted file mode 100644 index edebecfc7e8..00000000000 --- a/examples/sudokuhard.csv +++ /dev/null @@ -1,9 +0,0 @@ -5, 4, 0, 0, 0, 3, 1, 0, 0 -0, 8, 0, 4, 0, 1, 0, 0, 6 -0, 0, 0, 5, 0, 0, 0, 2, 0 -0, 7, 0, 0, 0, 0, 6, 0, 0 -0, 0, 4, 0, 0, 0, 9, 0, 0 -0, 0, 6, 0, 0, 0, 0, 3, 0 -0, 5, 0, 0, 0, 4, 0, 0, 0 -1, 0, 0, 2, 0, 8, 0, 7, 0 -0, 0, 2, 7, 0, 0, 0, 1, 9 diff --git a/examples/transp.jl b/examples/transp.jl index 13fb3888015..a231a573ed2 100644 --- a/examples/transp.jl +++ b/examples/transp.jl @@ -1,70 +1,56 @@ -############################################################################# # Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. -# -# Author: Louis Luangkesorn -# Date: Jan 30, 2015 -# -# Allocation of passenger cars to trains to minimize cars required or car-miles run. -# Based on -# Fourer, D.M. Gay and Brian W. Kernighan, A Modeling Language for Mathematical -# Programming, http://www.ampl.com/REFS/amplmod.ps.gz -# Appendix D ############################################################################# -using JuMP, Clp, Printf -const MOI = JuMP.MathOptInterface - -solver = Clp.Optimizer - -ORIG = ["GARY", "CLEV", "PITT"] -DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"] +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# -supply = [1400 2600 2900] -demand = [ 900 1200 600 400 1700 1100 1000] +using JuMP, GLPK, Test +const MOI = JuMP.MathOptInterface -@assert sum(supply) == sum(demand) +""" + example_transp() -cost = [ -39 14 11 14 16 82 8; -27 9 12 9 26 95 17; -24 14 17 13 28 99 20 -] +Allocation of passenger cars to trains to minimize cars required or car-miles +run. Based on: Fourer, D.M. Gay and Brian W. Kernighan, A Modeling Language for +Mathematical Programming, http://www.ampl.com/REFS/amplmod.ps.gz Appendix D. -m = Model(with_optimizer(solver)) +Author: Louis Luangkesorn +Date: Jan 30, 2015 +""" +function example_transp() + ORIG = ["GARY", "CLEV", "PITT"] + DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"] -@variable(m, Trans[i in 1:length(ORIG), j in 1:length(DEST)] >= 0) + supply = [1_400, 2_600, 2_900] + demand = [900, 1_200, 600, 400, 1_700, 1_100, 1_000] + @assert sum(supply) == sum(demand) -@objective(m, Min, sum(cost[i,j] * Trans[i,j] for i in 1:length(ORIG), j in 1:length(DEST))) + cost = [ + 39 14 11 14 16 82 8; + 27 9 12 9 26 95 17; + 24 14 17 13 28 99 20 + ] -@constraint(m, [i in 1:length(ORIG)], sum(Trans[i,j] for j in 1:length(DEST)) == supply[i]) + model = Model(with_optimizer(GLPK.Optimizer)) -@constraint(m, [j in 1:length(DEST)], sum(Trans[i,j] for i in 1:length(ORIG)) == demand[j]) + @variable(model, trans[1:length(ORIG), 1:length(DEST)] >= 0) + @objective(model, Min, sum(cost[i, j] * trans[i, j] for i in 1:length(ORIG), j in 1:length(DEST))) -println("Solving original problem...") -JuMP.optimize!(m) + @constraint(model, [i in 1:length(ORIG)], + sum(trans[i, j] for j in 1:length(DEST)) == supply[i]) + @constraint(model, [j in 1:length(DEST)], + sum(trans[i, j] for i in 1:length(ORIG)) == demand[j]) -term_status = JuMP.termination_status(m) -primal_status = JuMP.primal_status(m) -is_optimal = term_status == MOI.Optimal + JuMP.optimize!(model) -if is_optimal - @printf("Optimal!\n") - @printf("Objective value: %d\n", JuMP.objective_value(m)) - @printf("Transpotation:\n") - for j in 1:length(DEST) - @printf("\t%s", DEST[j]) - end - @printf("\n") - for i in 1:length(ORIG) - @printf("%s", ORIG[i]) - for j in 1:length(DEST) - @printf("\t%d", JuMP.value(Trans[i,j])) - end - @printf("\n") - end -else - @printf("No solution\n") + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) == 196200.0 end + +example_transp() diff --git a/examples/urban_plan.jl b/examples/urban_plan.jl new file mode 100644 index 00000000000..ef480596679 --- /dev/null +++ b/examples/urban_plan.jl @@ -0,0 +1,74 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +############################################################################# +# JuMP +# An algebraic modeling language for Julia +# See http://github.com/JuliaOpt/JuMP.jl +############################################################################# + +using JuMP, GLPK, Test +const MOI = JuMP.MathOptInterface + +""" + example_urban_plan() + +An "urban planning" problem. Based on +http://www.puzzlor.com/2013-08_UrbanPlanning.html +""" +function example_urban_plan() + model = Model(with_optimizer(GLPK.Optimizer)) + + # x is indexed by row and column + @variable(model, 0 <= x[1:5, 1:5] <= 1, Int) + + # y is indexed by R or C, the points, and an index in 1:5. Note how JuMP + # allows indexing on arbitrary sets. + rowcol = ["R", "C"] + points = [5, 4, 3, -3, -4, -5] + @variable(model, 0 <= y[rowcol, points, 1:5] <= 1, Int) + + # Objective - combine the positive and negative parts + @objective(model, Max, sum( + 3 * (y["R", 3, i] + y["C", 3, i]) + + 1 * (y["R", 4, i] + y["C", 4, i]) + + 1 * (y["R", 5, i] + y["C", 5, i]) + - 3 * (y["R", -3, i] + y["C", -3, i]) + - 1 * (y["R", -4, i] + y["C", -4, i]) + - 1 * (y["R", -5, i] + y["C", -5, i]) + for i in 1:5) + ) + + # Constrain the number of residential lots + @constraint(model, sum(x) == 12) + + # Add the constraints that link the auxiliary y variables to the x variables + for i = 1:5 + @constraints(model, begin + # Rows + y["R", 5, i] <= 1 / 5 * sum(x[i, :]) # sum = 5 + y["R", 4, i] <= 1 / 4 * sum(x[i, :]) # sum = 4 + y["R", 3, i] <= 1 / 3 * sum(x[i, :]) # sum = 3 + y["R", -3, i] >= 1 - 1 / 3 * sum(x[i, :]) # sum = 2 + y["R", -4, i] >= 1 - 1 / 2 * sum(x[i, :]) # sum = 1 + y["R", -5, i] >= 1 - 1 / 1 * sum(x[i, :]) # sum = 0 + # Columns + y["C", 5, i] <= 1 / 5 * sum(x[:, i]) # sum = 5 + y["C", 4, i] <= 1 / 4 * sum(x[:, i]) # sum = 4 + y["C", 3, i] <= 1 / 3 * sum(x[:, i]) # sum = 3 + y["C", -3, i] >= 1 - 1 / 3 * sum(x[:, i]) # sum = 2 + y["C", -4, i] >= 1 - 1 / 2 * sum(x[:, i]) # sum = 1 + y["C", -5, i] >= 1 - 1 / 1 * sum(x[:, i]) # sum = 0 + end) + end + + # Solve it + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == MOI.OPTIMAL + @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT + @test JuMP.objective_value(model) ≈ 14.0 +end + +example_urban_plan() diff --git a/examples/urbanplan.jl b/examples/urbanplan.jl deleted file mode 100644 index 9058ba5733c..00000000000 --- a/examples/urbanplan.jl +++ /dev/null @@ -1,84 +0,0 @@ -# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. -############################################################################# -# JuMP -# An algebraic modeling langauge for Julia -# See http://github.com/JuliaOpt/JuMP.jl -############################################################################# -# urbanplan.jl -# -# An "urban planning" problem. -# Based on -# http://www.puzzlor.com/2013-08_UrbanPlanning.html -############################################################################# - -using JuMP, GLPK -const MOI = JuMP.MathOptInterface - -solver = GLPK.Optimizer - -function SolveUrban() - - m = Model(with_optimizer(solver)) - - # x is indexed by row and column - @variable(m, 0 <= x[1:5, 1:5] <= 1, Int) - - # y is indexed by R or C, and the points - # JuMP allows indexing on arbitrary sets - rowcol = ["R","C"] - points = [+5,+4,+3,-3,-4,-5] - @variable(m, 0 <= y[rowcol, points, 1:5] <= 1, Int) - - # Objective - combine the positive and negative parts - @objective(m, Max, sum( - 3*(y["R", 3, i] + y["C", 3, i]) - + 1*(y["R", 4, i] + y["C", 4, i]) - + 1*(y["R", 5, i] + y["C", 5, i]) - - 3*(y["R", -3, i] + y["C", -3, i]) - - 1*(y["R", -4, i] + y["C", -4, i]) - - 1*(y["R", -5, i] + y["C", -5, i]) for i in 1:5)) - - # Constrain the number of residential lots - @constraint(m, sum(x[i,j] for i in 1:5, j in 1:5) == 12) - - # Add the constraints that link the auxiliary y variables - # to the x variables - # Rows - for i = 1:5 - @constraint(m, y["R", 5,i] <= 1/5*sum(x[i,j] for j in 1:5)) # sum = 5 - @constraint(m, y["R", 4,i] <= 1/4*sum(x[i,j] for j in 1:5)) # sum = 4 - @constraint(m, y["R", 3,i] <= 1/3*sum(x[i,j] for j in 1:5)) # sum = 3 - @constraint(m, y["R",-3,i] >= 1-1/3*sum(x[i,j] for j in 1:5)) # sum = 2 - @constraint(m, y["R",-4,i] >= 1-1/2*sum(x[i,j] for j in 1:5)) # sum = 1 - @constraint(m, y["R",-5,i] >= 1-1/1*sum(x[i,j] for j in 1:5)) # sum = 0 - end - # Columns - for j = 1:5 - @constraint(m, y["C", 5,j] <= 1/5*sum(x[i,j] for i in 1:5)) # sum = 5 - @constraint(m, y["C", 4,j] <= 1/4*sum(x[i,j] for i in 1:5)) # sum = 4 - @constraint(m, y["C", 3,j] <= 1/3*sum(x[i,j] for i in 1:5)) # sum = 3 - @constraint(m, y["C",-3,j] >= 1-1/3*sum(x[i,j] for i in 1:5)) # sum = 2 - @constraint(m, y["C",-4,j] >= 1-1/2*sum(x[i,j] for i in 1:5)) # sum = 1 - @constraint(m, y["C",-5,j] >= 1-1/1*sum(x[i,j] for i in 1:5)) # sum = 0 - end - - # Solve it - JuMP.optimize!(m) - - term_status = JuMP.termination_status(m) - primal_status = JuMP.primal_status(m) - is_optimal = term_status == MOI.Optimal - - if ! is_optimal - error("The solver did not find an optimal solution.") - end - - # Print results - println("Best objective: $(round(JuMP.objective_value(m)))") - println(JuMP.value.(x)) -end - -SolveUrban()