feat: add objective offset field c0 to MILP

src/MOI_wrapper.jl

@@ -208,6 +208,7 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
     end
     if max_sense
         c .*= -one(T)
+        obj_constant = -obj_constant
     end
 
     dest.sets = cache.constraints.sets
@@ -216,7 +217,7 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
     lc = cache.constraints.constants.lower
     uc = cache.constraints.constants.upper
 
-    milp = MILP(; c, lv, uv, A, lc, uc)
+    milp = MILP(; c, lv, uv, A, lc, uc, c0 = obj_constant)
 
     algorithm = pop!(dest.options, :algorithm, PDLP)
 
@@ -246,8 +247,8 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
             + sum(safeprod_left.(lv, positive_part.(dest.z)))
             - sum(safeprod_left.(uv, negative_part.(dest.z)))
     )
-    dest.obj_value = (max_sense ? -raw_obj : raw_obj) + obj_constant
-    dest.dual_obj_value = (max_sense ? -raw_dual_obj : raw_dual_obj) + obj_constant
+    dest.obj_value = max_sense ? -raw_obj : raw_obj
+    dest.dual_obj_value = max_sense ? -(raw_dual_obj + milp.c0) : (raw_dual_obj + milp.c0)
     dest.solve_time = stats.time_elapsed
 
     cts = stats.termination_status

src/components/preconditioning.jl

@@ -60,7 +60,7 @@ end
 
 function precondition(milp::MILP, prec::Preconditioner)
     (;
-        c, lv, uv, A, At, lc, uc,
+        c, lv, uv, A, At, lc, uc, c0,
         int_var, var_names, dataset, name, path,
     ) = milp
     (; D1, D2) = prec
@@ -79,6 +79,7 @@ function precondition(milp::MILP, prec::Preconditioner)
         At = At_p,
         lc = lc_p,
         uc = uc_p,
+        c0,
         D1 = new_prec.D1,
         D2 = new_prec.D2,
         int_var,

src/problems/milp.jl

@@ -11,7 +11,7 @@ Represent a Mixed Integer Linear Program in "cuPDLPx form":
     MILP(;
         c, lv, uv, A, lc, uc,
         At=sametype_transpose(A),
-        [D1, D2, int_var, var_names, dataset, name, path]
+        [c0, D1, D2, int_var, var_names, dataset, name, path]
     )
 
 # Fields
@@ -39,6 +39,8 @@ struct MILP{
     lc::V
     "constraint upper bound"
     uc::V
+    "objective offset (constant added to cᵀx)"
+    c0::T
     "left preconditioner"
     D1::Diagonal{T, V}
     "right preconditioner"
@@ -62,6 +64,7 @@ struct MILP{
             At = sametype_transpose(A),
             lc,
             uc,
+            c0 = false,
             D1 = Diagonal(one!(similar(lc))),
             D2 = Diagonal(one!(similar(lv))),
             int_var = zero!(similar(c, Bool)),
@@ -93,6 +96,8 @@ struct MILP{
             throw(ArgumentError("Abstract type parameter"))
         end
 
+        c0 = T(c0)
+
         common_backend(c, lv, uv, A, At, lc, uc, D1, D2)
 
         if isempty(name) && !isempty(path)
@@ -107,6 +112,7 @@ struct MILP{
             At,
             lc,
             uc,
+            c0,
             D1,
             D2,
             int_var,
@@ -126,6 +132,7 @@ Construct a [`MILP`](@ref) from a `QPSData` object generated by [QPSReader.jl](h
 function MILP(qps::QPSData; kwargs...)
     return MILP(;
         c = qps.c,
+        c0 = qps.c0,
         lv = qps.lvar,
         uv = qps.uvar,
         A = sparse(qps.arows, qps.acols, qps.avals, length(qps.lcon), length(qps.lvar)),
@@ -211,6 +218,7 @@ function Base.isapprox(m1::MILP, m2::MILP; kwargs...)
             isapprox(m1.At, m2.At; kwargs...) &&
             isapprox(m1.lc, m2.lc; kwargs...) &&
             isapprox(m1.uc, m2.uc; kwargs...) &&
+            isapprox(m1.c0, m2.c0; kwargs...) &&
             isapprox(m1.D1, m2.D1; kwargs...) &&
             isapprox(m1.D2, m2.D2; kwargs...) &&
             m1.int_var == m2.int_var &&

src/problems/modify.jl

@@ -11,7 +11,7 @@ end
 
 function set_eltype(::Type{T}, milp::MILP) where {T}
     (;
-        c, lv, uv, A, At, lc, uc, D1, D2,
+        c, lv, uv, A, At, lc, uc, c0, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
     return MILP(;
@@ -22,6 +22,7 @@ function set_eltype(::Type{T}, milp::MILP) where {T}
         At = set_eltype(T, At),
         lc = set_eltype(T, lc),
         uc = set_eltype(T, uc),
+        c0 = T(c0),
         D1 = set_eltype(T, D1),
         D2 = set_eltype(T, D2),
         int_var,
@@ -51,7 +52,7 @@ end
 
 function set_indtype(::Type{Ti}, milp::MILP) where {Ti}
     (;
-        c, lv, uv, A, At, lc, uc, D1, D2,
+        c, lv, uv, A, At, lc, uc, c0, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
     return MILP(;
@@ -62,6 +63,7 @@ function set_indtype(::Type{Ti}, milp::MILP) where {Ti}
         At = set_indtype(Ti, At),
         lc,
         uc,
+        c0,
         D1,
         D2,
         int_var,
@@ -79,7 +81,7 @@ Convert the sparse matrices inside `milp` using constructor `M`.
 """
 function set_matrix_type(::Type{M}, milp::MILP) where {M}
     (;
-        c, lv, uv, A, At, lc, uc, D1, D2,
+        c, lv, uv, A, At, lc, uc, c0, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
     A_M = M(A)
@@ -94,6 +96,7 @@ function set_matrix_type(::Type{M}, milp::MILP) where {M}
         At = At_M,
         lc = adapt(backend, lc),
         uc = adapt(backend, uc),
+        c0,
         D1 = adapt(backend, D1),
         D2 = adapt(backend, D2),
         int_var = adapt(backend, int_var),
@@ -110,7 +113,7 @@ end
 
 function Adapt.adapt_structure(to, milp::MILP)
     (;
-        c, lv, uv, A, At, lc, uc, D1, D2,
+        c, lv, uv, A, At, lc, uc, c0, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
     return MILP(;
@@ -121,6 +124,7 @@ function Adapt.adapt_structure(to, milp::MILP)
         At = adapt(to, At),
         lc = adapt(to, lc),
         uc = adapt(to, uc),
+        c0,
         D1 = adapt(to, D1),
         D2 = adapt(to, D2),
         int_var = adapt(to, int_var),

src/problems/solution.jl

@@ -37,7 +37,7 @@ end
 
 Compute the value of the linear objective of `milp` at solution vector `x`.
 """
-objective_value(x, milp::MILP) = dot(x, milp.c)
+objective_value(x, milp::MILP) = dot(x, milp.c) + milp.c0
 
 """
     PrimalDualSolution

test/components/preconditioning.jl

@@ -68,6 +68,7 @@ end
     @test isapprox(sol, sol_unp)
     @test !isapprox(sol, sol_p)
 
+    @test milp.c0 == milp_p.c0
     @test objective_value(sol.x, milp) ≈ objective_value(sol_p.x, milp_p)
     @test dot(sol.y, milp.A, sol.x) ≈ dot(sol_p.y, milp_p.A, sol_p.x)
     @test CoolPDLP.proj_box.(sol.x, milp.lv, milp.uv) ≈ prec.D2 * CoolPDLP.proj_box.(sol_p.x, milp_p.lv, milp_p.uv)

test/problems/milp.jl

@@ -83,4 +83,7 @@ end
     qps, path = read_instance(Netlib, netlib[1])
     milp = MILP(qps; path, dataset = "Netlib")
     @test milp ≈ milp
+    milp_offset = MILP(qps; path, dataset = "Netlib", c0 = 5.0)
+    @test milp_offset ≈ milp_offset
+    @test !(milp ≈ milp_offset)
 end

test/problems/solution.jl

@@ -20,6 +20,10 @@ using Test
     @test @test_warn "Constraints not satisfied" !is_feasible([0.0, 0.0], milp)
     @test @test_warn "Variable bounds not satisfied" !is_feasible([2.0, -1.0], milp)
     @test objective_value([1.0, 1.0], milp) == 3
+
+    milp_offset = MILP(; c, lv, uv, A, lc, uc, int_var, c0 = 10.0)
+    @test objective_value([1.0, 1.0], milp_offset) == 13
+    @test objective_value([1.0, 1.0], milp_offset) - objective_value([1.0, 1.0], milp) == 10
 end
 
 @testset "Comparison with JuMP" begin