Improved stability and flexibility

src/Oracle/oracle_mip.jl

@@ -41,6 +41,26 @@ mutable struct LindaOracleMIP <: AbstractLindaOracle
 
         oracle.index = index
 
+        (m, n) = size(A_sub)
+        m == length(row_lb) || throw(DimensionMismatch(
+            "A has $m rows but $(length(row_lb)) row lower bounds"
+        ))
+        m == length(row_ub) || throw(DimensionMismatch(
+            "A has $m rows but $(length(row_lb)) row upper bounds"
+        ))
+        n == length(col_lb) || throw(DimensionMismatch(
+            "A has $n cols but $(length(col_lb)) col lower bounds"
+        ))
+        n == length(col_ub) || throw(DimensionMismatch(
+            "A has $n cols but $(length(col_ub)) col upper bounds"
+        ))
+        n == length(vartypes) || throw(DimensionMismatch(
+            "A has $n cols but $(length(vartypes)) var types"
+        ))
+        n == length(costs) || throw(DimensionMismatch(
+            "A has $n cols but $(length(costs)) objective terms"
+        ))
+
         oracle.costs = costs
         oracle.A_link = A_link
         oracle.A_sub = A_sub
@@ -66,7 +86,8 @@ function query!(
     σ::Real;
     farkas::Bool=false,
     tol_reduced_cost::Float64=1.0e-6,
-    num_columns_max::Int=typemax(Int64)
+    num_columns_max::Int=typemax(Int64),
+    log=Dict()
 ) where{Tv<:Real}
 
     # Dimension checks

src/Oracle/oracle_pool.jl

@@ -46,7 +46,8 @@ function query!(
     farkas::Bool=false,
     tol_reduced_cost::Float64=1.0e-6,
     num_columns_max::Int=typemax(Int64),
-    prop_sp_priced_max::Float64=1.0
+    prop_sp_priced_max::Float64=1.0,
+    log::Dict=Dict()
 ) where{T1<:Real, T2<:Real}
 
     pool.new_columns = Set{Column}()
@@ -61,7 +62,6 @@ function query!(
 
     # price each sub-problem
     # go through sub-problems in random order
-    nsp_priced = 0
     perm = randperm(pool.n)
     for r in perm
 
@@ -72,13 +72,13 @@ function query!(
             farkas=farkas,
             tol_reduced_cost=tol_reduced_cost
         )
-        nsp_priced += 1
+        log[:nsp_priced] += 1
 
         # Check for infeasible sub-problem
         s = get_oracle_status(o)
         if s == PrimalInfeasible
             pool.status = PrimalInfeasible
-            @warn("Infeasible sub-problem")
+            @warn("Sub-problem $r is infeasible.")
             return pool.status
         end
 
@@ -108,10 +108,6 @@ function query!(
     # Pricing ended because all sub-problems were solved to optimality
     pool.status = Optimal
 
-    # println("\tNCols: ", length(pool.new_columns))
-    # println("\trc=: ", best_red_cost)
-    # println("\tPriced:", nsp_priced)
-
     return pool.status
 
 end

src/colgen.jl

@@ -6,37 +6,50 @@ Column-Generation algorithm.
 function solve_colgen!(
     env::LindaEnv,
     mp::LindaMaster{RMP},
-    oracle::Oracle.AbstractLindaOracle
+    oracle::Oracle.AbstractLindaOracle;
+    cg_log::Dict=Dict()
 ) where{RMP<:MPB.AbstractMathProgModel}
 
     # Pre-optimization stuff
     n_cg_iter::Int = 0
     time_start = time()
-    time_mp_total::Float64 = 0.0
-    time_sp_total::Float64 = 0.0
-    time_cg_total::Float64 = 0.0
+    cg_log[:time_mp_total] = 0.0
+    cg_log[:time_sp_total] = 0.0
+    cg_log[:time_cg_total] = 0.0
+
+    cg_log[:num_iter_bar] = 0  # Total number of barrier inner iterations
+    cg_log[:num_iter_spx] = 0  # Total number of simplex inner iterations
+
+    cg_log[:nsp_priced] = 0  # Number of calls to sub-problems
 
     if env[Val{:verbose}] == 1
-        println(" Itn    Primal Obj      Dual Obj        NCols    Time (s)")
+        println("Itn     Primal Obj        Dual Obj         NCols    MP(s)    SP(s)   Tot(s)  BarIter  SpxIter")
     end
 
     # Main CG loop
     while n_cg_iter < env[Val{:num_cgiter_max}]
+
+        # Check for time limit
+        if (time() - time_start) >= env[:time_limit]
+            env[Val{:verbose}] == 1 && println("Time limit reached.")
+            break
+        end
+
         # Solve RMP, update dual variables
         t0 = time()
         solve_rmp!(mp)
-        time_mp_total += time() - t0
+        cg_log[:time_mp_total] += time() - t0
 
         if mp.rmp_status == Optimal
             farkas=false
         elseif mp.rmp_status == PrimalInfeasible
             farkas=true
         elseif mp.rmp_status == PrimalUnbounded
-            println("Master Problem is unbounded.")
-            return mp.mp_status
+            env[Val{:verbose}] == 1 && println("Master Problem is unbounded.")
+            break
         else
-            error("RMP status $(mp.rmp_status) not handled. Exiting.")
-            return mp.mp_status
+            @warn("RMP status $(mp.rmp_status) not handled.")
+            break
         end
 
         # Price
@@ -45,9 +58,10 @@ function solve_colgen!(
             oracle, mp.π, mp.σ,
             farkas=farkas,
             tol_reduced_cost=env.tol_reduced_cost.val,
-            num_columns_max=env.num_columns_max.val
+            num_columns_max=env.num_columns_max.val,
+            log=cg_log
         )
-        time_sp_total += time() - t0
+        cg_log[:time_sp_total] += time() - t0
 
         cols = Oracle.get_new_columns(oracle)
         lagrange_lb = (
@@ -56,6 +70,9 @@ function solve_colgen!(
         )  # Compute Lagrange lower bound
         mp.dual_bound = lagrange_lb > mp.dual_bound ? lagrange_lb : mp.dual_bound
 
+        cg_log[:num_iter_bar] += MPB.getbarrieriter(mp.rmp)
+        cg_log[:num_iter_spx] += MPB.getsimplexiter(mp.rmp)
+
         # Log
         # Iteration count
         if env[Val{:verbose}] == 1
@@ -65,8 +82,13 @@ function solve_colgen!(
             @printf("%+16.7e", mp.dual_bound)
             # RMP stats
             @printf("%10.0f", mp.num_columns_rmp)  # number of columns in RMP
-            @printf("%9.2f", time_mp_total + time_sp_total)
+            @printf("%9.2f", cg_log[:time_mp_total])
+            @printf("%9.2f", cg_log[:time_sp_total])
+            @printf("%9.2f", time() - time_start)
+            @printf("%9d", cg_log[:num_iter_bar])
+            @printf("%9d", cg_log[:num_iter_spx])
             print("\n")
+            flush(Base.stdout)
         end
 
         # Check duality gap
@@ -77,28 +99,21 @@ function solve_colgen!(
         if mp_gap <= 10.0 ^-4
             mp.mp_status = Optimal
 
-            time_cg_total += time() - time_start
+            # time_cg_total += time() - time_start
             if env[Val{:verbose}] == 1
                 println("Root relaxation solved.")
-                println("Total time / MP: ", time_mp_total)
-                println("Total time / SP: ", time_sp_total)
-                println("Total time / CG: ", time_cg_total)
             end
 
-            return mp.mp_status
+            break
 
         elseif farkas && length(cols) == 0
 
             mp.mp_status = PrimalInfeasible
-            time_cg_total += time() - time_start
+            # time_cg_total += time() - time_start
             if env[Val{:verbose}] == 1
-                println("Master is infeasible.")
-                println("Total time / MP: ", time_mp_total)
-                println("Total time / SP: ", time_sp_total)
-                println("Total time / CG: ", time_cg_total)
+                println("Problem is infeasible.")
             end
-
-            return mp.mp_status
+            break
         else
             # add columns
             add_columns!(mp, cols)
@@ -107,11 +122,22 @@ function solve_colgen!(
         n_cg_iter += 1
     end
 
-    time_cg_total += time() - time_start
+    # Logs
+    cg_log[:time_cg_total] = time() - time_start
+    cg_log[:dual_bound] = mp.dual_bound
+    cg_log[:primal_bound] = mp.primal_lp_bound
+    cg_log[:n_cg_iter] = n_cg_iter
+    cg_log[:status] = mp.mp_status
+    cg_log[:num_cols_tot] = mp.num_columns_rmp
+
     if env[Val{:verbose}] == 1
-        println("Total time / MP: ", time_mp_total)
-        println("Total time / SP: ", time_sp_total)
-        println("Total time / CG: ", time_cg_total)
+        println()
+        @printf("Total time / MP: %.2fs\n", cg_log[:time_mp_total])
+        @printf("Total time / SP: %.2fs\n", cg_log[:time_sp_total])
+        @printf("Total time / CG: %.2fs\n", cg_log[:time_cg_total])
+        @printf("Inner barrier iterations: %d\n", cg_log[:num_iter_bar])
+        @printf("Inner simplex iterations: %d\n", cg_log[:num_iter_spx])
+        @printf("Pricing calls: %d\n", cg_log[:nsp_priced])
     end
 
     return mp.mp_status

src/master.jl

@@ -24,6 +24,7 @@ mutable struct LindaMaster{RMP<:MPB.AbstractMathProgModel}
             column is added to / removed from the RMP
 
     =#
+    num_var_link::Int       # Number of linking variables (including artificial ones)
     num_columns_rmp::Int  # Number of columns currently in RMP
     active_columns::Vector{Column}  # Active columns
 
@@ -87,18 +88,21 @@ mutable struct LindaMaster{RMP<:MPB.AbstractMathProgModel}
     #===========================================================================
         Constructor
     ===========================================================================#
-    
+
     function LindaMaster(
-        rmp::RMP,
         num_constr_cvxty::Int,
         num_constr_link::Int,
-        rhs_constr_link::AbstractVector{Tv}
-    ) where{RMP<:MPB.AbstractMathProgModel, Tv<:Real}
+        rhs_constr_link::AbstractVector{Tv},
+        num_var_link::Int,
+        initial_columns::Vector{Tc},
+        rmp::RMP,
+    ) where{RMP<:MPB.AbstractMathProgModel, Tv<:Real, Tc<:Column}
 
         # Dimension check
         n = MPB.numvar(rmp)
-        n == (2 * num_constr_link) || throw(ErrorException(
-            "RMP has $(n) variables instead of $(2*num_constr_link)"
+        ncols = length(initial_columns)
+        n == (num_var_link + ncols) || throw(ErrorException(
+            "RMP has $(n) variables instead of $(num_var_link + ncols)"
         ))
         m = MPB.numconstr(rmp)
         m == (num_constr_cvxty + num_constr_link) || throw(ErrorException(
@@ -111,18 +115,23 @@ mutable struct LindaMaster{RMP<:MPB.AbstractMathProgModel}
         # Instanciate model
         mp = new{RMP}()
 
-        mp.num_columns_rmp = 0
-        mp.active_columns = Vector{Column}(undef, 0)
+        # RMP columns
+        mp.num_var_link = num_var_link
+        mp.num_columns_rmp = ncols
+        mp.active_columns = deepcopy(initial_columns)
 
+        # RMP constraints
         mp.num_constr_cvxty = num_constr_cvxty
         mp.num_constr_link = num_constr_link
-        mp.rhs_constr_link = rhs_constr_link
+        mp.rhs_constr_link = copy(rhs_constr_link)
         mp.π = Vector{Float64}(undef, num_constr_link)
         mp.σ = Vector{Float64}(undef, num_constr_cvxty)
 
+        # Other RMP info
         mp.rmp = rmp
         mp.rmp_status = Unknown
 
+        # Column-Generation info
         mp.mp_status = Unknown
         mp.primal_lp_bound = Inf
         mp.primal_ip_bound = Inf
@@ -148,43 +157,49 @@ function LindaMaster(
     rmp = MPB.LinearQuadraticModel(lp_solver)
     # Add constraints
     for r in 1:num_constr_cvxty
-        MPB.addconstr!(rmp, Vector{Float64}(), Vector{Float64}(), 1.0, 1.0)
+        MPB.addconstr!(rmp, Float64[], Float64[], 1.0, 1.0)
     end
     for i in 1:num_constr_link
         MPB.addconstr!(
             rmp,
-            Vector{Float64}(),
-            Vector{Float64}(),
+            Float64[],
+            Float64[],
             rhs_constr_link[i],
             rhs_constr_link[i]
         )
 
     end
 
     # Add artificial variables
+    num_var_link = 0
     for i in 1:num_constr_link
         if senses[i] == '='
             # Artificial slack and surplus
             MPB.addvar!(rmp, [num_constr_cvxty+i], [1.0], 0.0, 0.0, 10^4) # slack
             MPB.addvar!(rmp, [num_constr_cvxty+i], [-1.0], 0.0, 0.0, 10^4) # surplus
+            num_var_link += 2
         elseif senses[i] == '<'
             # Regular slack ; artificial surplus
             MPB.addvar!(rmp, [num_constr_cvxty+i], [1.0], 0.0, Inf, 0.0) # slack
             MPB.addvar!(rmp, [num_constr_cvxty+i], [-1.0], 0.0, 0.0, 10^4) # surplus
+            num_var_link += 2
         elseif senses[i] == '>'
             # Artificial slack ; regular surplus
             MPB.addvar!(rmp, [num_constr_cvxty+i], [1.0], 0.0, 0.0, 10^4) # slack
             MPB.addvar!(rmp, [num_constr_cvxty+i], [-1.0], 0.0, Inf, 0.0) # surplus
+            num_var_link += 2
         else
-            error("Wrong input $(senses[i])")
+            error("Unknown sense for constraint $i: `$(senses[i])`")
         end
     end
 
     mp = LindaMaster(
-        rmp,
         num_constr_cvxty,
         num_constr_link,
-        rhs_constr_link
+        rhs_constr_link,
+        num_var_link,
+        Column[],
+        rmp
     )
 
     return mp
@@ -212,7 +227,7 @@ Solve Restricted Master Problem, and update current dual iterate.
 function solve_rmp!(master::LindaMaster)
 
     MPB.optimize!(master.rmp)
-    rmp_status = Status(Val{MPB.status(master.rmp)})
+    rmp_status = Status(Val(MPB.status(master.rmp)))
     master.rmp_status = rmp_status
 
     # Update dual iterate
@@ -232,6 +247,12 @@ function solve_rmp!(master::LindaMaster)
 
         # update dual variables
         y = MPB.getinfeasibilityray(master.rmp)
+        # Normalize dual ray, for stability purposes
+        nrm = norm(y[(master.num_constr_cvxty+1):(master.num_constr_cvxty+master.num_constr_link)], 1)
+        if nrm > 1e-6
+            y ./= nrm
+        end
+
         master.σ .= y[1:master.num_constr_cvxty]
         master.π .= y[(master.num_constr_cvxty+1):(master.num_constr_cvxty+master.num_constr_link)]
 

src/status.jl

@@ -21,8 +21,8 @@ These Status codes are based on those of MathOptInterface
 )
 
 Status(::T) where T = Unknown
-Status(s::Symbol) = Status(Val{s})
+Status(s::Symbol) = Status(Val(s))
 
-Status(::Type{Val{:Infeasible}}) = PrimalInfeasible
-Status(::Type{Val{:Unbounded}}) = PrimalUnbounded
-Status(::Type{Val{:Optimal}}) = Optimal
+Status(::Val{:Infeasible}) = PrimalInfeasible
+Status(::Val{:Unbounded}) = PrimalUnbounded
+Status(::Val{:Optimal}) = Optimal