diff --git a/Project.toml b/Project.toml index 5c52909..d31c5db 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RingStarProblems" uuid = "bed4ad48-c7de-4cc5-bbc0-d98473e6c0f4" authors = ["jkhamphousone "] -version = "0.1.10" +version = "0.2.0" [deps] Cairo = "159f3aea-2a34-519c-b102-8c37f9878175" diff --git a/README.md b/README.md index 30d30d8..9988a10 100755 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ The references: Introduce the Resilient Ring Star Problem (called 1-R-RSP). The package can solve 1-R-RSP thanks to: - - A Branch-and-Benders-cut algorithm (refered to as B&BC) + - A Branch-and-Benders-cut algorithm (referred to as B&BC) - An Integer Linear Programming model (ILP) ## Ring Star Problem @@ -26,7 +26,7 @@ When setting `backup_factor=0` or `tildeV=0`, 1-R-RSP reduces to RSP # Requirements -[JuMP.jl](https://github.com/jump-dev/JuMP.jl) must be installed. You can use CPLEX, GLPK, Gurobi, Xpress or any [solver that MOI.supports Lazy Constraints callback in JuMP](https://jump.dev/JuMP.jl/stable/manual/callbacks/#Available-solvers). +[JuMP.jl](https://github.com/jump-dev/JuMP.jl) must be installed. You can use CPLEX, GLPK, Gurobi, Xpress or any [solver that supports Lazy Constraints callback in JuMP](https://jump.dev/JuMP.jl/stable/manual/callbacks/#Available-solvers). # Installation ```julia @@ -45,7 +45,7 @@ julia> pars = RSP.SolverParameters( s_ij = RSP.Euclidian(), # star costs r_ij = RSP.Euclidian(), # ring costs backup_factor = 0.01, # backup_factor c'=0.01c and d'=0.01c - alpha = 3, # See [`Labbé et al., 2004`](ttps://doi.org/10.1002/net.10114) + alpha = 3, # See [`Labbé et al., 2004`](ttps://doi.org/10.1002/net.10114) tildeV = 100, # uncertain nodes set writeresults = RSP.WHTML(), # output results locally, WLocal(), WHTML() or no output false plotting = false, # plot results to subfolder src/plots/results/ @@ -74,6 +74,23 @@ julia> symbolinstance = :berlin52 julia> RSP.rspoptimize(pars, symbolinstance, Gurobi.Optimizer) ``` +## Solving with nodes coordinates + +Either: +```julia +julia> x = 1:10 +julia> y = rand(1:10, 10) +julia> RSP.rspoptimize(pars, x, y, Gurobi.Optimizer) +``` + +Or: +```julia +julia> xycoors = tuple.(1:10, rand(1:10, 10)) +julia> RSP.rspoptimize(pars, xycoors, Gurobi.Optimizer) +``` + + + ## Plotting To plot the solutions in the folder ext/results diff --git "a/ext/results/TinyInstance_12_2____\316\261-3.0_benders_Poly()/TinyInstance_12_2____\316\261-3.0_benders_Poly()__F=0.0_UB=903.0).pdf" "b/ext/results/TinyInstance_12_2____\316\261-3.0_benders_Poly()/TinyInstance_12_2____\316\261-3.0_benders_Poly()__F=0.0_UB=903.0).pdf" index 3fcd96e..70ab9f2 100644 Binary files "a/ext/results/TinyInstance_12_2____\316\261-3.0_benders_Poly()/TinyInstance_12_2____\316\261-3.0_benders_Poly()__F=0.0_UB=903.0).pdf" and "b/ext/results/TinyInstance_12_2____\316\261-3.0_benders_Poly()/TinyInstance_12_2____\316\261-3.0_benders_Poly()__F=0.0_UB=903.0).pdf" differ diff --git "a/ext/results/eil51________________\316\261-3.0_benders_LP()/eil51________________\316\261-3.0_benders_LP()__F=0.0_UB=1278.0).pdf" "b/ext/results/eil51________________\316\261-3.0_benders_LP()/eil51________________\316\261-3.0_benders_LP()__F=0.0_UB=1278.0).pdf" index f05b283..bbbf956 100644 Binary files "a/ext/results/eil51________________\316\261-3.0_benders_LP()/eil51________________\316\261-3.0_benders_LP()__F=0.0_UB=1278.0).pdf" and "b/ext/results/eil51________________\316\261-3.0_benders_LP()/eil51________________\316\261-3.0_benders_LP()__F=0.0_UB=1278.0).pdf" differ diff --git a/src/benders_rrsp.jl b/src/benders_rrsp.jl index 423f722..9a1861c 100755 --- a/src/benders_rrsp.jl +++ b/src/benders_rrsp.jl @@ -123,10 +123,10 @@ function rrspcreatebenders_modellazy(filename, inst, pars; optimizer) function f(x, y) - sum(sum(c[i, j] * x[i, j] for j in V if i < j; init = 0) for i in V) + - sum(c[1, i] * x[i, n+1] for i ∈ 2:n) + - sum(sum(d[i, j] * y[i, j] for j in V if i != j) for i in V) + - sum(o[i] * y[i, i] for i in V) + sum(sum(c[i, j] * x[i, j] for j in V if i < j; init = 0) for i in V; init = 0) + + sum(c[1, i] * x[i, n+1] for i ∈ 2:n ; init = 0) + + sum(sum(d[i, j] * y[i, j] for j in V if i != j; init = 0) for i in V; init = 0) + + sum(o[i] * y[i, i] for i in V; init = 0) end bestsol, bestobjval = three(inst.n, inst.o, inst.c, inst.d, tildeV) diff --git a/src/instance.jl b/src/instance.jl index ff5ff83..2b6767f 100755 --- a/src/instance.jl +++ b/src/instance.jl @@ -127,8 +127,9 @@ function printinst(inst::RRSPInstance, pars) end function createinstance_rrsp(filename, α, pars) + random_filepath = - eval(@__DIR__) * "/instances/Instances_journal_article/RAND/$(filename[1]).dat" + eval(@__DIR__) * "/instances/Instances_journal_article/RAND/$(filename[2]).dat" if pars.nrand > 0 && !isfile(random_filepath) n = pars.nrand data = [1:1000 rand(1:1000, 1000) rand(1:1000, 1000)] @@ -142,47 +143,21 @@ function createinstance_rrsp(filename, α, pars) if α == 0 error("Please define α") end - x_coors = data[1:n, 2] - y_coors = data[1:n, 3] - @assert pars.r_ij isa Euclidian - c = Dict( - (e[1], e[2]) => ceil_labbe( - distance([x_coors[e[1]], y_coors[e[1]]], [x_coors[e[2]], y_coors[e[2]]]), - ) for e in E - ) - - - d = Dict( - (a[1], a[2]) => - a[1] == a[2] ? 0 : round(1 / rand(Uniform(n / 2, 3n / 2)), digits = 5) - for a in A - ) - if pars.s_ij isa Euclidian - d = Dict( - (a[1], a[2]) => ceil_labbe( - distance( - [x_coors[a[1]], y_coors[a[1]]], - [x_coors[a[2]], y_coors[a[2]]], - ), - ) for a in A - ) - end - for i ∈ 2:n - c[i, n+1] = c[1, i] - end - c[1, n+1] = 0 + xs = data[1:n, 2] + ys = data[1:n, 3] + c, d = createweights(xs, ys, A, pars.r_ij, pars.s_ij) o = zeros(Float64, n) if pars.o_i == RandomInterval o = rand(RandomInterval.a:RandomInterval.b, n) - elseif pars.o_i == 1 - o = ones(Float64, n) + elseif typeof(pars.o_i) == UInt + o = pars.o_i*ones(Float64, n) end open(random_filepath, "w") do f println(f, "$n 0.0 $α") for i ∈ 1:n - println(f, "$i $(x_coors[i]) $(y_coors[i])") + println(f, "$i $(xs[i]) $(ys[i])") end println(f, "opening costs") for i ∈ 1:n @@ -205,14 +180,12 @@ function createinstance_rrsp(filename, α, pars) for kv in d d′[kv[1]] = d[kv[1]] * pars.backup_factor end - return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, x_coors, y_coors) + return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, xs, ys) elseif pars.nrand == 0 data = readdlm(filename[1]) n = filename[2] - @assert α != 0 - V = 1:n tildeV = 2:Int(ceil(n * pars.tildeV / 100)) E = [(i, j) for i in V, j in V if i < j] @@ -224,8 +197,8 @@ function createinstance_rrsp(filename, α, pars) o = zeros(Float64, n) if pars.o_i == RandomInterval o = rand(RandomInterval.a:RandomInterval.b, n) - elseif pars.o_i == 1 - o = ones(Float64, n) + elseif typeof(pars.o_i) == UInt + o = pars.o_i*ones(Float64, n) end if filename[1][end-12:end] == "brazil58.tsp2" @@ -270,15 +243,14 @@ function createinstance_rrsp(filename, α, pars) if filename[1][end-7:end] == "120.tsp2" shift_n = 414 - 7 end - succesfullread = true - x_coors, y_coors = readcoordinates(data) + xs, ys = readcoordinates(data) c = Dict( (e[1], e[2]) => ceil_labbe( distance( - [x_coors[e[1]], y_coors[e[1]]], - [x_coors[e[2]], y_coors[e[2]]], + [xs[e[1]], ys[e[1]]], + [xs[e[2]], ys[e[2]]], ) * (α), ) for e in E ) @@ -293,8 +265,8 @@ function createinstance_rrsp(filename, α, pars) d = Dict( (a[1], a[2]) => ceil_labbe( distance( - [x_coors[a[1]], y_coors[a[1]]], - [x_coors[a[2]], y_coors[a[2]]], + [xs[a[1]], ys[a[1]]], + [xs[a[2]], ys[a[2]]], ) * (10 - α), ) for a in A ) @@ -307,7 +279,7 @@ function createinstance_rrsp(filename, α, pars) for kv in d d′[kv[1]] = d[kv[1]] * pars.backup_factor end - return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, x_coors, y_coors) + return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, xs, ys) end else @@ -318,19 +290,11 @@ function createinstance_rrsp(filename, α, pars) E = [(i, j) for i in V, j in V if i < j] A = [(i, j) for i in V, j in V] Ac = [(i, j) for i in V, j in V if i != j] - x_coors, y_coors = readcoordinates(data) + xs, ys = readcoordinates(data) o = data[n+3, 1:n] - @assert pars.r_ij isa Euclidian - c = Dict( - (e[1], e[2]) => ceil_labbe( - distance([x_coors[e[1]], y_coors[e[1]]], [x_coors[e[2]], y_coors[e[2]]]), - ) for e in E - ) - for i ∈ 2:n - c[i, n+1] = c[1, i] - end - c[1, n+1] = 0 + c, d = createweights(xs, ys, A, pars.r_ij, pars.s_ij) + d_data = data[n+5:2n+4, 1:n] d = Dict((a[1], a[2]) => d_data[a[1], a[2]] for a in A) c′ = Dict{Tuple{Int,Int},Float64}() @@ -342,42 +306,152 @@ function createinstance_rrsp(filename, α, pars) d′[kv[1]] = d[kv[1]] * pars.backup_factor end - return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, x_coors, y_coors) + return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, xs, ys) + end +end + +""" + createinstance_rrsp(xycoordinates, α, pars) + `x` and `y` are coordinates of points + + Return an RRSPInstance +""" +function createinstance_rrsp(instdataname::Tuple{Vector{Tuple{Int, Int}},Int}, α, pars) + n = length(instdataname[2]) + + xs, ys = first.(instdataname[1]), last.(instdataname[1]) + + V = 1:n + tildeV = 2:Int(ceil(n * pars.tildeV / 100)) + E = [(i, j) for i in V, j in V if i < j] + A = [(i, j) for i in V, j in V] + Ac = [(i, j) for i in V, j in V if i != j] + + o = zeros(Float64, n) + if pars.o_i == RandomInterval + o = rand(RandomInterval.a:RandomInterval.b, n) + elseif typeof(pars.o_i) == UInt + o = pars.o_i*ones(Float64, n) + end + + c = Dict( + (e[1], e[2]) => ceil_labbe( + distance( + [xs[e[1]], ys[e[1]]], + [xs[e[2]], ys[e[2]]], + ) * (α), + ) for e in E + ) + + + for i ∈ 2:n + c[i, n+1] = c[1, i] end + c[1, n+1] = 0 + + + d = Dict( + (a[1], a[2]) => ceil_labbe( + distance( + [xs[a[1]], ys[a[1]]], + [xs[a[2]], ys[a[2]]], + ) * (10 - α), + ) for a in A + ) + + c′ = Dict{Tuple{Int,Int},Float64}() + d′ = Dict{Tuple{Int,Int},Float64}() + for kv in c + c′[kv[1]] = c[kv[1]] * pars.backup_factor + end + for kv in d + d′[kv[1]] = d[kv[1]] * pars.backup_factor + end + return RRSPInstance(n, V, tildeV, pars.F, o, α, c, c′, d, d′, xs, ys) end + + + + + """ readcoordinates(data) + first line and column of data contains the number of nodes and is ignored Return two vectors of coordinates read from data """ function readcoordinates(data) n_start = 2 - x_coors = Float64[] - y_coors = Float64[] - while length(x_coors) == 0 + xs = Float64[] + ys = Float64[] + while length(xs) == 0 try if data[end, 2] == "" - x_coors = Float64.(data[n_start:end-1, 2]) + xs = Float64.(data[n_start:end-1, 2]) else - x_coors = Float64.(data[n_start:end, 2]) + xs = Float64.(data[n_start:end, 2]) end catch n_start += 1 end end n_start = 2 - while length(y_coors) == 0 + while length(ys) == 0 try if data[end, 2] == "" - y_coors = Float64.(data[n_start:end-1, 3]) + ys = Float64.(data[n_start:end-1, 3]) else - y_coors = Float64.(data[n_start:end, 3]) + ys = Float64.(data[n_start:end, 3]) end catch n_start += 1 end end - return x_coors, y_coors + return xs, ys end + + +function createweights(xs, ys, A, r_ij, s_ij) + if r_ij isa RandomInterval + c = Dict( + (a[1], a[2]) => + a[1] == a[2] ? 0 : rand(pars.r_ij[1]:pars.r_ij[2]) + for a in A + ) + end + if r_ij isa Euclidian + c = Dict( + (e[1], e[2]) => ceil_labbe( + distance([xs[e[1]], ys[e[1]]], [xs[e[2]], ys[e[2]]]), + ) for e in E + ) + end + + if s_ij isa RandomInterval + d = Dict( + (a[1], a[2]) => + a[1] == a[2] ? 0 : rand(pars.r_ij[1]:pars.r_ij[2]) + for a in A + ) + end + if s_ij isa Euclidian + d = Dict( + (a[1], a[2]) => ceil_labbe( + distance( + [xs[a[1]], ys[a[1]]], + [xs[a[2]], ys[a[2]]], + ), + ) for a in A + ) + end + + for i ∈ 2:n + c[i, n+1] = c[1, i] + end + c[1, n+1] = 0 + + + return c, d +end \ No newline at end of file diff --git a/src/main.jl b/src/main.jl index 072cbf1..aa3534b 100755 --- a/src/main.jl +++ b/src/main.jl @@ -4,7 +4,7 @@ """ rspoptimize(pars::SolverParameters, symbolinstance::Symbol, optimizer, solutionchecker = false) -Return exit code 0 +Return exit succesful code 0 """ function rspoptimize( pars::SolverParameters, @@ -39,6 +39,43 @@ function rspoptimize( end + + +""" + rspoptimize(pars::SolverParameters, xycoordinates::Vector{Tuple{Int,Int}}, optimizer, solutionchecker = false) + + Return exit succesful code 0 +""" +function rspoptimize( + pars::SolverParameters, + xycoordinates::Vector{Tuple{Int,Int}}, + optimizer, + solutionchecker = false, +) + + instdataname = "xycoordinates_$(length(xycoordinates))", (xycoordinates, length(xycoordinates)) + + if pars.redirect_stdio + # redirect terminal outputs/stdio txycoordinateso file + now_file = Dates.format(Dates.now(), "yyyy-mm-dd_HHhMM") + now_folder = Dates.format(Dates.now(), "yyyy-mm-dd") + output_path = joinpath(@__DIR__, "debug", "stdio", "$now_folder") + mkpath(output_path) + redirect_stdio( + stdout = "$output_path/stdout_$(symbolinstance)_$now_file.txt", + stderr = "$output_path/stderr_$(symbolinstance)_$now_file.txt", + ) do + main(pars, instdataname, optimizer, solutionchecker) + GC.gc() + end + end + main(pars, instdataname, optimizer, solutionchecker) + GC.gc() + return 0 +end + + + function main(pars::SolverParameters, instdataname, optimizer, solutionchecker = false) @@ -58,13 +95,13 @@ function main(pars::SolverParameters, instdataname, optimizer, solutionchecker = nstr_random = pars.nrand > 0 ? "_" * "$(pars.nrand)" : "" if instdataname[1] == :RandomInstance - instdataname[2] = "$(instdataname[1])$(nstr_random)_o=$(replace(pars.o_i,":"=>"-"))_$(pars.r_ij)_s($(pars.s_ij))_ID$(rand_inst_id)" + instdataname[2][1] = "$(instdataname[1])$(nstr_random)_o=$(replace(pars.o_i,":"=>"-"))_$(pars.r_ij)_s($(pars.s_ij))_ID$(rand_inst_id)" if pars.r_ij == pars.s_ij - instdataname[2] = "$(instdataname[1])$(nstr_random)_o$(replace(pars.o_i,":"=>"-"))_rs=l_ij_ID$(rand_inst_id)" + instdataname[2][1] = "$(instdataname[1])$(nstr_random)_o$(replace(pars.o_i,":"=>"-"))_rs=l_ij_ID$(rand_inst_id)" end end - + inst = createinstance_rrsp(instdataname[2], α, pars) println("\nInstance: ", "$(instdataname[1])") @@ -118,12 +155,20 @@ function main(pars::SolverParameters, instdataname, optimizer, solutionchecker = ilp_table = read_ilp_table(input_filepath, pars.plot_id) end - if pars.plotting && pars.writeresults != "" + if pars.plotting && pars.writeresults if pars.solve_mod == ILP() - perform_plot(pars, inst, instdataname[1], ilp_table, true) + if ilp_table.sol.n != 0 + perform_plot(pars, inst, instdataname[1], ilp_table, true) + else + @info "No feasible solution has been found by ILP within the solving time, can not perform plotting" + end end if pars.solve_mod == BranchBendersCut() - perform_plot(pars, inst, instdataname[1], benders_table, false) + if benders_table.sol.n != 0 + perform_plot(pars, inst, instdataname[1], benders_table, false) + else + @info "No feasible solution has been found by B&BC within the solving time, can not perform plotting" + end end end if pars.writeresults != "" @@ -143,3 +188,20 @@ function main(pars::SolverParameters, instdataname, optimizer, solutionchecker = GC.gc() return 0 end + +""" + rspoptimize(pars::SolverParameters, x::Vector{Int}, y::Vector{Int}, optimizer, solutionchecker = false) + + Return exit succesful code 0 +""" +rspoptimize( + pars::SolverParameters, + x, y, + optimizer, + solutionchecker, +) = rspoptimize( + pars, + tuple.(x, y), + optimizer, + solutionchecker, +) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 8d4c6ec..da9f1d2 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,10 +16,10 @@ include("aqua.jl") pars = SolverParameters( - solve_mod = BranchBendersCut(), # ILP, or BranchBendersCut + solve_mod = BranchBendersCut(), # ILP(), or BranchBendersCut() F = 7, # total failing time F in days per year, see [`PhD manuscript`](https://theses.hal.science/tel-04319443) sp_solve = Poly(), - o_i = 0, # opening costs + o_i = 1, # opening costs s_ij = Euclidian(), # star costs r_ij = Euclidian(), # ring costs alpha = 3, # See [Labbé et al., 2004](ttps://doi.org/10.1002/net.10114) @@ -36,18 +36,41 @@ include("aqua.jl") include("solutionchecker.jl") + xycoordinates = tuple.(1:10, rand(1:10, 10)) + + + @test rspoptimize( + pars, + xycoordinates, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), + true, + ) == 0 + + + x = 1:10 + y = rand(1:10, 10) + @test rspoptimize( + pars, + x, + y, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), + true, + ) == 0 + @test rspoptimize( pars, :TinyInstance_10_3, optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), true, ) == 0 + @test rspoptimize( pars, :Instance_15_1_0_3_1, - optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), - true, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), + true, ) == 0 + @test rspoptimize( pars, :eil51, @@ -60,16 +83,16 @@ include("aqua.jl") @test rspoptimize( pars, - :TinyInstance_12_2, - optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), + :TinyInstance_12_2, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2), true, ) == 0 pars.sp_solve = LP() pars.redirect_stdio = true @test rspoptimize( pars, - :eil51, - optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2, "tm_lim" => 15_000), + :eil51, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2, "tm_lim" => 15_000), true, ) == 0 @@ -82,6 +105,15 @@ include("aqua.jl") ) == 0 pars.tildeV = 0 + pars.o_i = 0 + @test rspoptimize( + pars, + :TinyInstance_10_3, + optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2, "tm_lim" => 1_000), + true, + ) == 0 + + pars.solve_mod = ILP() @test rspoptimize( pars, :TinyInstance_10_3,