Julia JuMP example for cutting planes

problem.jl

using JuMP
using Gurobi

# define type to hold cutting planes info
# ---------------------------------------
type cuts
    RHS::Vector{Vector{Float64}}
    Coeff::Vector{Vector{Float64}}
    n::Int

    function cuts()
        new([ Float64[] for i in 1:0 ],
            [ Float64[] for i in 1:0 ],
            0)
    end
end


function main()

    time_sol = 0.0
    time_add = 0.0

    # repeats multiple time the same procedure
    for j in 1:3000

        # initialize type
        c = cuts()

        #initialize JuMP.Model
        m = createMod()

        # solve initial model with no cuts
        solve(m)

        # loop fo adding cutting planes
        # -----------------------------
        for i in 1:3

            # obtain new cut from oracle
            RHS, Coeff = oracle(m)

            # store cut into structure
            pushCut!(c, RHS, Coeff)

            # add cut to the model
            tic()
            addCut!(m,c)
            time_add += toq()

            # solve updated model
            tic()
            solve(m)
            time_sol += toq()
        end
    end

    @show time_sol
    @show time_add

    nothing
end


# creates model
function createMod()

    # model data difinitions
    # ----------------------
    NG = 10

    Gmax = 10.0*collect(1:NG)
    Cost = 1.0*collect(1:NG)
    Xmax = 200.0

    # difine LP
    # ----------
    m = Model( solver = GurobiSolver(OutputFlag = 0))

    @variable(m, y[i in 1:20] >= 0.0)

    @variable(m, 0.0 <= x[i in 1:20] <= Xmax)

    @variable(m, 0.0 <= g[j in 1:NG] <= Gmax[j])

    @constraint(m, sum{x[i], i in 1:20} + sum{ g[j], j in 1:NG} == 900 )

    @objective(m, Min, sum{y[i] , i in 1:20} +  sum{ Cost[j]*g[j] , j in 1:NG})

    return m
end

function oracle(m)

    Xstar = getvalue( getvariable(m,:x) )

    Ystar = 1.0 + ( Xstar - 30.0 ).^2

    Coeff = 2*( Xstar - 30.0 )

    RHS = Ystar -  Coeff*30.0

    return RHS, Coeff
end

function addCut!(m,cuts)

    # load variables into scope
    x = getvariable(m,:x)
    y = getvariable(m,:y)

    # add contraint to the model
    @constraint(m, cstr[i in 1:20], y[i] - cuts.Coeff[end][i]*x[i] >= cuts.RHS[end][i] )

    nothing
end

function pushCut!(cuts, RHS, Coeff)

    # adds ne cut to the cut pool
    push!(cuts.RHS, RHS)
    push!(cuts.Coeff, Coeff)
    cuts.n += 1

    nothing
end