Support for compute_conflict!(model) #134

juliohm · 2023-03-08T17:46:14Z

Is it easy to add support for compute_conflict!(model) as described in JuMP's docs?

mtanneau · 2023-03-08T19:25:22Z

If you're working with linear programs (which is the only class supported by Tulip), then a conflict can be identified from a dual proof of infeasibility (which you can retrieve by getting the dual solution of an infeasible problem).

I don't know how the internals of compute_conflict! works, but I imagine the following should work:

solve the problem. If feasible, then there's no conflict.
If infeasible, get a dual proof of infeasibility
Every constraint that has a non-zero dual in that solution is in the conflict

Because Tulip is an interior-point solver, you get the following limitations:

the returned conflict is not irreducible (you might remove a constraint and the problem is still infeasible)
the dual solution might have close-to-zero-but-not-zero values, and you could end up with a very large conflict.

juliohm · 2023-03-08T19:55:26Z

I am trying to get access to the dual with the JuMP interface, but having problems still. My problem has linear objective and the following constraints:

A JuMP Model
Minimization problem with:
Variables: 416
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 103 constraints
`AffExpr`-in-`MathOptInterface.GreaterThan{Float64}`: 9 constraints
`AffExpr`-in-`MathOptInterface.Interval{Float64}`: 402 constraints

juliohm · 2023-03-08T19:56:00Z

The interval constraints are an issue in this case?

juliohm · 2023-03-08T20:00:22Z

I checked has_duals(model) and it returns true. I then tried dual.(LowerBoundRef.(model[:Vp])) but it returns an error. Is this because my constraints are interval, not one of LowerBoundRef, UpperBoundRef, or FixRef?

mtanneau · 2023-03-08T21:43:33Z

Can you give me a reproducible example?

juliohm · 2023-03-09T17:13:43Z

Here is a MWE where the infeasibility comes from two sources:

The consumer C1 is requesting 200 units of product A but the producer P1 only produces 100 units.
The consumer C2 is requesting 100 units of product B but the route R3 does not support more than 50 units.

Can we somehow highlight these two facts from the dual problem?

function mwe()
	model = Model(Tulip.Optimizer)

	### PRODUCERS
	
	pset = zip(["P1","P2"], ["A","B"])
	@variable(model, Vp[pset])
	@constraint(model, [0.,0.] .≤ Vp .≤ [100.,100.])

	### STATIONS

	sset = [("S1","A","IN"), ("S1","A","OUT"), ("S1","A-S","OUT"),
	        ("S2","B","IN"), ("S2","B","OUT")]
	@variable(model, Vs[sset])
	for s in ["S1","S2"]
		ssetin  = filter(x -> first(x) == s && last(x) == "IN",  sset)
		ssetout = filter(x -> first(x) == s && last(x) == "OUT", sset)
		@constraint(model, 0.0 ≤ sum(Vs[ssetin])  ≤ Inf)
		@constraint(model, 0.0 ≤ sum(Vs[ssetout]) ≤ Inf)
	end

	### BENEFICIATION
	
	bset = zip(["S1","S1","S2"], ["A","A","B"], ["A","A-S","B"])
	@variable(model, Vb[bset])
	@constraint(model, Vb .≥ 0)
	
	### CONSUMERS
	
	cset = zip(["C1","C2"], ["A","B"])
	@variable(model, Vc[cset])
	@constraint(model, [200.,100.] .≤ Vc .≤ [Inf,Inf])

	### ROUTES
	
	rset = [("R1","A","P1","S1"),("R2","A","S1","C1"),
	        ("R3","B","P2","S2"),("R4","B","S2","C2")]
	@variable(model, Vr[rset])
	@constraint(model, 0. ≤ sum(Vr[[("R1","A","P1","S1")]]) ≤ Inf)
	@constraint(model, 0. ≤ sum(Vr[[("R2","A","S1","C1")]]) ≤ Inf)
	@constraint(model, 0. ≤ sum(Vr[[("R3","B","P2","S2")]]) ≤ 50.)
	@constraint(model, 0. ≤ sum(Vr[[("R4","B","S2","C2")]]) ≤ Inf)

	### MASS BALANCE
	
	for (orig, prod) in pset
		rcon = filter(x -> x[3] == orig && x[2] == prod, rset)
		@constraint(model, Vp[(orig, prod)] == sum(Vr[rcon]))
	end

	for (station, prod, side) in sset
		if side == "IN"
			rcon = filter(x -> x[2] == prod && x[4] == station, rset)
			@constraint(model, Vs[(station, prod, side)] == sum(Vr[rcon]))
		elseif side == "OUT"
			rcon = filter(x -> x[3] == station && x[2] == prod, rset)
			@constraint(model, Vs[(station, prod, side)] == sum(Vr[rcon]))
		end
	end

	for (dest, prod) in cset
		rcon = filter(x -> x[2] == prod && x[4] == dest, rset)
		@constraint(model, Vc[(dest, prod)] == sum(Vr[rcon]))
	end

	@constraint(model,
		        Vs[("S1","A","IN")] == sum(Vb[[("S1","A","A"),("S1","A","A-S")]]))
	@constraint(model,
		        Vs[("S1","A","IN")] == sum(Vs[[("S1","A","OUT"),("S1","A-S","OUT")]]))
	@constraint(model,
		        Vs[("S2","B","IN")] == sum(Vb[[("S2","B","B")]]))
	@constraint(model,
		        Vs[("S2","B","IN")] == sum(Vs[[("S2","B","OUT")]]))

	### OBJECTIVE
	
	rcosts = sum(Vr[r] for r in rset)
	bcosts = sum(Vb[b] for b in bset)
	cvalue = sum(Vc[c] for c in cset)
	
	@objective(model, Min, rcosts + bcosts - cvalue)
	
	model
end

mtanneau · 2023-03-09T20:00:38Z

Here is how I (manually) print a conflict from an infeasibility proof

using JuMP
using Tulip
using Printf

# Add MWE code here
function mwe()
    ...
end

function conflict()
    model = mwe()
    optimize!(model)

    st = termination_status(model)
    dst = dual_status(model)

    if st == INFEASIBLE && dst == INFEASIBILITY_CERTIFICATE
        ctypes = list_of_constraint_types(model)
        for (F, S) in ctypes
            cons = all_constraints(model, F, S)
            for i in eachindex(cons)
                isassigned(cons, i) || continue
                con = cons[i]
                λ = dual(con)  # get dual of constraint
                # if dual is non-zero, the constraint is part of the conflict
                abs(λ) > 1e-4 && println("[λ = $λ] ", con) 
            end
        end
    else
        println("Termination status is $(st) (dual status $dst); no conflict to print")
    end
    return nothing
end

This code outputs the following:

[λ = 1.970914335410348] Vp[("P1", "A")] - Vr[("R1", "A", "P1", "S1")] = 0.0
[λ = -1.8082052821554653] Vs[("S1", "A", "IN")] - Vr[("R1", "A", "P1", "S1")] = 0.0
[λ = 1.310308006900506] Vs[("S1", "A", "OUT")] - Vr[("R2", "A", "S1", "C1")] = 0.0
[λ = 1.3103080069936996] Vs[("S1", "A-S", "OUT")] = 0.0
[λ = -1.3707587381583282] Vs[("S2", "B", "IN")] - Vr[("R3", "B", "P2", "S2")] = 0.0
[λ = 0.8477120338037647] Vs[("S2", "B", "OUT")] - Vr[("R4", "B", "S2", "C2")] = 0.0
[λ = -1.227954732915119] Vc[("C1", "A")] - Vr[("R2", "A", "S1", "C1")] = 0.0
[λ = -0.6537346401063114] Vc[("C2", "B")] - Vr[("R4", "B", "S2", "C2")] = 0.0
[λ = 0.27551734649770443] Vs[("S1", "A", "IN")] - Vb[("S1", "A", "A")] - Vb[("S1", "A", "A-S")] = 0.0
[λ = 1.4055739642095522] Vs[("S1", "A", "IN")] - Vs[("S1", "A", "OUT")] - Vs[("S1", "A-S", "OUT")] = 0.0
[λ = -1.0] Vs[("S2", "B", "IN")] - Vb[("S2", "B", "B")] = 0.0
[λ = 1.076829350357877] Vs[("S2", "B", "IN")] - Vs[("S2", "B", "OUT")] = 0.0
[λ = 0.2755173464754485] Vb[("S1", "A", "A")] ≥ 0.0
[λ = 0.2755173464754485] Vb[("S1", "A", "A-S")] ≥ 0.0
[λ = -1.9709143354282603] Vp[("P1", "A")] ∈ [0.0, 100.0]
[λ = 0.127113971372492] Vs[("S1", "A", "IN")] ∈ [0.0, Inf]
[λ = 0.09526595721585247] Vs[("S1", "A", "OUT")] + Vs[("S1", "A-S", "OUT")] ∈ [0.0, Inf]
[λ = 0.29392938776936095] Vs[("S2", "B", "IN")] ∈ [0.0, Inf]
[λ = 0.22911731651387665] Vs[("S2", "B", "OUT")] ∈ [0.0, Inf]
[λ = 1.2279547328428626] Vc[("C1", "A")] ∈ [200.0, Inf]
[λ = 0.6537346400692975] Vc[("C2", "B")] ∈ [100.0, Inf]
[λ = 0.16270905321115217] Vr[("R1", "A", "P1", "S1")] ∈ [0.0, Inf]
[λ = 0.08235327389692947] Vr[("R2", "A", "S1", "C1")] ∈ [0.0, Inf]
[λ = -1.370758738172468] Vr[("R3", "B", "P2", "S2")] ∈ [0.0, 50.0]
[λ = 0.1939773936510756] Vr[("R4", "B", "S2", "C2")] ∈ [0.0, Inf]

All these constraints have non-zero dual and therefore form part of the conflict. I just printed the dual proof of infeasibility, so it's not guaranteed to be irreducible (i.e., it's possible that removing one of these constraints doesn't render the model feasible).
For interval constraints, if the dual is positive (resp. negative), the lower (resp. upper) bound is the one that's active.

juliohm · 2023-03-09T20:07:16Z

Is this result similar to what we get with primal_feasibility_report in JuMP.jl? Em qui., 9 de mar. de 2023 17:00, mtanneau ***@***.***> escreveu:

…

Here is how I (manually) print a conflict from an infeasibility proof using JuMPusing Tulipusing Printf # Add MWE code herefunction mwe() ...end function conflict() model = mwe() optimize!(model) st = termination_status(model) dst = dual_status(model) if st == INFEASIBLE && dst == INFEASIBILITY_CERTIFICATE ctypes = list_of_constraint_types(model) for (F, S) in ctypes cons = all_constraints(model, F, S) for i in eachindex(cons) isassigned(cons, i) || continue con = cons[i] λ = dual(con) # get dual of constraint # if dual is non-zero, the constraint is part of the conflict abs(λ) > 1e-4 && println("[λ = $λ] ", con) end end else println("Termination status is $(st) (dual status $dst); no conflict to print") end return nothingend This code outputs the following: [λ = 1.970914335410348] Vp[("P1", "A")] - Vr[("R1", "A", "P1", "S1")] = 0.0 [λ = -1.8082052821554653] Vs[("S1", "A", "IN")] - Vr[("R1", "A", "P1", "S1")] = 0.0 [λ = 1.310308006900506] Vs[("S1", "A", "OUT")] - Vr[("R2", "A", "S1", "C1")] = 0.0 [λ = 1.3103080069936996] Vs[("S1", "A-S", "OUT")] = 0.0 [λ = -1.3707587381583282] Vs[("S2", "B", "IN")] - Vr[("R3", "B", "P2", "S2")] = 0.0 [λ = 0.8477120338037647] Vs[("S2", "B", "OUT")] - Vr[("R4", "B", "S2", "C2")] = 0.0 [λ = -1.227954732915119] Vc[("C1", "A")] - Vr[("R2", "A", "S1", "C1")] = 0.0 [λ = -0.6537346401063114] Vc[("C2", "B")] - Vr[("R4", "B", "S2", "C2")] = 0.0 [λ = 0.27551734649770443] Vs[("S1", "A", "IN")] - Vb[("S1", "A", "A")] - Vb[("S1", "A", "A-S")] = 0.0 [λ = 1.4055739642095522] Vs[("S1", "A", "IN")] - Vs[("S1", "A", "OUT")] - Vs[("S1", "A-S", "OUT")] = 0.0 [λ = -1.0] Vs[("S2", "B", "IN")] - Vb[("S2", "B", "B")] = 0.0 [λ = 1.076829350357877] Vs[("S2", "B", "IN")] - Vs[("S2", "B", "OUT")] = 0.0 [λ = 0.2755173464754485] Vb[("S1", "A", "A")] ≥ 0.0 [λ = 0.2755173464754485] Vb[("S1", "A", "A-S")] ≥ 0.0 [λ = -1.9709143354282603] Vp[("P1", "A")] ∈ [0.0, 100.0] [λ = 0.127113971372492] Vs[("S1", "A", "IN")] ∈ [0.0, Inf] [λ = 0.09526595721585247] Vs[("S1", "A", "OUT")] + Vs[("S1", "A-S", "OUT")] ∈ [0.0, Inf] [λ = 0.29392938776936095] Vs[("S2", "B", "IN")] ∈ [0.0, Inf] [λ = 0.22911731651387665] Vs[("S2", "B", "OUT")] ∈ [0.0, Inf] [λ = 1.2279547328428626] Vc[("C1", "A")] ∈ [200.0, Inf] [λ = 0.6537346400692975] Vc[("C2", "B")] ∈ [100.0, Inf] [λ = 0.16270905321115217] Vr[("R1", "A", "P1", "S1")] ∈ [0.0, Inf] [λ = 0.08235327389692947] Vr[("R2", "A", "S1", "C1")] ∈ [0.0, Inf] [λ = -1.370758738172468] Vr[("R3", "B", "P2", "S2")] ∈ [0.0, 50.0] [λ = 0.1939773936510756] Vr[("R4", "B", "S2", "C2")] ∈ [0.0, Inf] All these constraints have non-zero dual and therefore form part of the conflict. I just printed the dual proof of infeasibility, so it's not guaranteed to be irreducible (i.e., it's possible that removing one of these constraints doesn't render the model feasible). For interval constraints, if the dual is positive (resp. negative), the lower (resp. upper) bound is the one that's active. — Reply to this email directly, view it on GitHub <#134 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3KATTVZFISOX2TGIN3W3IZHDANCNFSM6AAAAAAVUCCQV4> . You are receiving this because you authored the thread.Message ID: ***@***.***>

mtanneau · 2023-03-10T01:33:15Z

Is this result similar to what we get with primal_feasibility_report in JuMP.jl?

No. JuMP's primal_feasibility_report takes as input a primal solution, and reports constraint violations.

The code I shared identifies a set of constraints that, together, render the model infeasible. There is no primal solution involved in that.

juliohm · 2023-03-10T12:21:16Z

Tank you @mtanneau for crafting the conflict function above, I will experiment with it here. Appreciate your help ❤️

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Support for compute_conflict!(model) #134

Support for compute_conflict!(model) #134

juliohm commented Mar 8, 2023

mtanneau commented Mar 8, 2023

juliohm commented Mar 8, 2023

juliohm commented Mar 8, 2023

juliohm commented Mar 8, 2023

mtanneau commented Mar 8, 2023

juliohm commented Mar 9, 2023

mtanneau commented Mar 9, 2023

juliohm commented Mar 9, 2023 via email

mtanneau commented Mar 10, 2023

juliohm commented Mar 10, 2023

Support for compute_conflict!(model) #134

Support for compute_conflict!(model) #134

Comments

juliohm commented Mar 8, 2023

mtanneau commented Mar 8, 2023

juliohm commented Mar 8, 2023

juliohm commented Mar 8, 2023

juliohm commented Mar 8, 2023

mtanneau commented Mar 8, 2023

juliohm commented Mar 9, 2023

mtanneau commented Mar 9, 2023

juliohm commented Mar 9, 2023 via email

mtanneau commented Mar 10, 2023

juliohm commented Mar 10, 2023