import spanning tree code from old repository

Author: BatyLeo

docs/src/benchmarks/two_stage_spanning_tree.md

@@ -0,0 +1,28 @@
+# Problem statement
+
+We consider a two-stage stochastic variant of the classic [minimum spanning tree problem](https://en.wikipedia.org/wiki/Minimum_spanning_tree).
+
+Rather than immediately constructing a spanning tree and incurring a cost ``c_e`` for each selected edge in the tree, we instead can build only a partial tree (forest) during the first stage and paying first stage costs ``c_e`` for the selected edges. Then, second stage costs ``d_e`` are revealed and replace first stage costs. The task then involves completing the first stage forest into a spanning tree.
+
+The objective is to minimize the total incurred cost in expectation.
+
+## Instance
+Let ``G = (V,E)`` be an undirected **graph**, and ``S`` be a finite set of **scenarios**.
+
+For each edge ``e`` in ``E``, we have a **first stage cost** ``c_e\in\mathbb{R}``.
+
+For each edge ``e`` in ``E`` and scenario ``s`` in ``S``, we have a **second stage cost** ``d_{es}\in\mathbb{R}``.
+
+# MIP formulation
+Unlike the regular minimum spanning tree problem, this two-stage variant is NP-hard.
+However, it can still be formulated as linear program with binary variables, and exponential number of constraints.
+```math
+\begin{array}{lll}
+\min\limits_{y, z}\, & \sum\limits_{e\in E}c_e y_e + \frac{1}{|S|}\sum\limits_{s \in S}d_{es}z_{es} & \\
+\mathrm{s.t.}\, & \sum\limits_{e\in E}y_e + z_{es} = |V| - 1, & \forall s \in S\\
+& \sum\limits_{e\in E(Y)} y_e + z_{es} \leq |Y| - 1,\quad & \forall \emptyset \subsetneq Y \subsetneq V,\, \forall s\in S\\
+& y_e\in \{0, 1\}, & \forall e\in E\\
+& z_{es}\in \{0, 1\}, & \forall e\in E, \forall s\in S
+\end{array}
+```
+where ``y_e`` is a binary variable indicating if ``e`` is in the first stage solution, and ``z_{es}`` is a binary variable indicating if ``e`` is in the second stage solution for scenario ``s``.

src/TwoStageSpanningTree/algorithms/anticipative.jl

@@ -0,0 +1,16 @@
+"""
+$TYPEDSIGNATURES
+
+Compute an anticipative solution for given scenario.
+"""
+function anticipative_solution(instance::TwoStageSpanningTreeInstance, scenario::Int=1)
+    (; graph, first_stage_costs, second_stage_costs) = instance
+    scenario_second_stage_costs = @view second_stage_costs[:, scenario]
+    weights = min.(first_stage_costs, scenario_second_stage_costs)
+    (; value, tree) = kruskal(graph, weights)
+
+    slice = first_stage_costs .<= scenario_second_stage_costs
+    y = tree[slice]
+    z = tree[.!slice]
+    return (; value, y, z)
+end

src/TwoStageSpanningTree/algorithms/benders_decomposition.jl

@@ -0,0 +1,164 @@
+function separate_benders_cut(instance::TwoStageSpanningTreeInstance, y, s; MILP_solver, tol=1e-5)
+	(; graph, second_stage_costs) = instance
+
+	E = ne(graph)
+
+	columns = BitVector[]
+
+	# Feasibility cut
+	model = Model(MILP_solver)
+
+	@variable(model, dummy, Bin)
+
+	@variable(model, νₛ <= 1)
+	@variable(model, 0 <= μₛ[e in 1:E] <= 1)
+
+	@objective(model, Max, νₛ + sum(y[e] * μₛ[e] for e in 1:E))
+
+	function feasibility_callback(cb_data)
+		μ_val = callback_value.(cb_data, μₛ)
+		ν_val = callback_value(cb_data, νₛ)
+
+		weights = -μ_val
+		val, tree = kruskal(graph, weights)
+
+		push!(columns, tree)
+		
+		if val + tol < ν_val
+			new_constraint = @build_constraint(
+				- sum(μₛ[e] for e in 1:E if tree[e]) - νₛ >= 0
+			)
+			MOI.submit(
+				model, MOI.LazyConstraint(cb_data), new_constraint
+			)
+		end
+	end
+
+	set_attribute(model, MOI.LazyConstraintCallback(), feasibility_callback)
+	optimize!(model)
+
+	if objective_value(model) > tol
+		return false, value.(νₛ), value.(μₛ), objective_value(model)
+	end
+	
+	# Else, optimality cut
+	optimality_model = Model(MILP_solver)
+
+	@variable(optimality_model, dummy, Bin)
+
+	@variable(optimality_model, νₛ)
+	@variable(optimality_model, μₛ[e in 1:E] >= 0)
+
+	@objective(
+		optimality_model, Max,
+		νₛ + sum(y[e] * μₛ[e] for e in 1:E) - sum(second_stage_costs[e, s] * y[e] for e in 1:E)
+	)
+
+	for tree in columns
+		@constraint(
+			optimality_model,
+			sum(second_stage_costs[e, s] - μₛ[e] for e in 1:E if tree[e]) >= νₛ
+		)
+	end
+
+	function my_callback_function(cb_data)
+		μ_val = callback_value.(cb_data, μₛ)
+		ν_val = callback_value(cb_data, νₛ)
+
+		weights = second_stage_costs[:, s] .- μ_val
+
+		val, tree = kruskal(graph, weights)
+
+		if val - ν_val + tol < 0
+			new_constraint = @build_constraint(
+				sum(second_stage_costs[e, s] - μₛ[e] for e in 1:E if tree[e]) >= νₛ
+			)
+			MOI.submit(
+				optimality_model, MOI.LazyConstraint(cb_data), new_constraint
+			)
+		end
+	end
+
+	set_attribute(optimality_model, MOI.LazyConstraintCallback(), my_callback_function)
+
+	optimize!(optimality_model)
+
+	# If primal feasible, add an optimality cut
+	@assert termination_status(optimality_model) != DUAL_INFEASIBLE
+	return true, value.(νₛ), value.(μₛ), objective_value(optimality_model)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Returns the optimal solution using a Benders decomposition algorithm.
+"""
+function benders_decomposition(
+    instance::TwoStageSpanningTreeInstance;
+    MILP_solver=GLPK.Optimizer,
+	tol=1e-6,
+	verbose=true
+)
+	(; graph, first_stage_costs, second_stage_costs) = instance
+	E = ne(graph)
+	S = nb_scenarios(instance)
+	
+    model = Model(MILP_solver)
+    @variable(model, y[e in 1:E], Bin)
+    @variable(
+        model,
+        θ[s in 1:S] >= sum(min(0, second_stage_costs[e, s]) for e in 1:E)
+    )
+    @objective(
+        model,
+        Min,
+        sum(first_stage_costs[e] * y[e] for e in 1:E) + sum(θ[s] for s in 1:S) / S
+    )
+
+    # current_scenario = 0
+    callback_counter = 0
+    function benders_callback(cb_data)
+        if callback_counter % 10 == 0
+            verbose && @info("Benders iteration: $(callback_counter)")
+        end
+        callback_counter += 1
+
+        y_val = callback_value.(cb_data, y)
+		θ_val = callback_value.(cb_data, θ)
+
+        for current_scenario in 1:S
+            optimality_cut, ν_val, μ_val =
+				separate_benders_cut(instance, y_val, current_scenario; MILP_solver)
+
+			# If feasibility cut
+            if !optimality_cut
+                new_feasibility_cut = @build_constraint(
+                    ν_val + sum(μ_val[e] * y[e] for e in 1:E) <= 0
+                )
+                MOI.submit(
+                    model,
+					MOI.LazyConstraint(cb_data),
+					new_feasibility_cut
+                )
+
+                return nothing
+            end
+
+			# Else, optimality cut
+			if θ_val[current_scenario] + tol < ν_val + sum(μ_val[e] * y_val[e] for e in 1:E) -
+				sum(second_stage_costs[e, current_scenario] * y_val[e] for e in 1:E)
+				con = @build_constraint(
+					θ[current_scenario] >=
+						ν_val + sum(μ_val[e] * y[e] for e in 1:E) - sum(second_stage_costs[e, current_scenario] * y[e] for e in 1:E)
+				)
+				MOI.submit(model, MOI.LazyConstraint(cb_data), con)
+				return nothing
+			end
+        end
+    end
+
+    set_attribute(model, MOI.LazyConstraintCallback(), benders_callback)
+    optimize!(model)
+
+	return solution_from_first_stage_forest(value.(y) .> 0.5, instance)
+end

src/TwoStageSpanningTree/algorithms/column_generation.jl

@@ -0,0 +1,111 @@
+"""
+$TYPEDSIGNATURES
+
+Solves the linear relaxation using a column generation algorithm.
+"""
+function column_generation(instance; MILP_solver=GLPK.Optimizer, tol=1e-6, verbose=true)
+	(; graph, first_stage_costs, second_stage_costs) = instance
+
+	V = nv(graph)
+	E = ne(graph)
+	S = nb_scenarios(instance)
+	
+	model = Model(MILP_solver)
+
+	@variable(model, dummy, Bin) # dummy binary variable to activate callbacks
+
+	@variable(model, ν[s in 1:S])
+	@variable(model, μ[e in 1:E, s in 1:S])
+
+	@objective(model, Max, sum(ν[s] for s in 1:S))
+
+	@constraint(
+		model, [e in 1:E, s in 1:S],
+		second_stage_costs[e, s] / S - μ[e, s] >= 0
+	)
+
+	@constraint(
+		model, [e in 1:E],
+		first_stage_costs[e] - sum(μ[e, s] for s in 1:S) >= 0
+	)
+
+	trees = BitVector[]
+
+	for s in 1:S
+		_, dummy_tree = kruskal(graph, min.(first_stage_costs, second_stage_costs[:, s]))
+		@constraint(
+        	model,
+        	ν[s] <= sum(μ[e, s] for e in 1:E if dummy_tree[e])
+    	)
+		push!(trees, dummy_tree)
+	end
+
+	callback_counter = 0
+	function my_callback_function(cb_data)
+		callback_counter += 1
+
+        ν_val = callback_value.(cb_data, ν)
+        μ_val = callback_value.(cb_data, μ)
+
+        for s in 1:S
+            val, T = kruskal(graph, @view μ_val[:, s])
+
+			if val + tol < ν_val[s]
+				push!(trees, T)
+                new_constraint = @build_constraint(
+                    ν[s] <= sum(μ[e, s] for e in 1:E if T[e])
+                )
+                MOI.submit(
+                    model, MOI.LazyConstraint(cb_data), new_constraint
+                )
+            end
+        end
+    end
+
+	set_attribute(model, MOI.LazyConstraintCallback(), my_callback_function)
+
+	optimize!(model)
+	verbose && @info "Optimal solution found after $callback_counter cuts"
+	return (; value=objective_value(model), ν=value.(ν), μ=value.(μ), columns=trees)
+end
+
+function column_heuristic(instance, columns; MILP_solver=GLPK.Optimizer)
+	(; graph, first_stage_costs, second_stage_costs) = instance
+	E = ne(graph)
+	S = nb_scenarios(instance)
+	T = length(columns)
+
+	model = Model(MILP_solver)
+
+	@variable(model, y[e in 1:E], Bin)
+	@variable(model, z[e in 1:E, s in 1:S], Bin)
+
+	@variable(model, λ[t in 1:T, s in 1:S], Bin)
+
+	@objective(
+		model, Min,
+		sum(first_stage_costs[e] * y[e] for e in 1:E) + sum(second_stage_costs[e, s] * z[e, s] for e in 1:E for s in 1:S) / S
+	)
+
+	@constraint(model, [s in 1:S], sum(λ[t, s] for t in 1:T) == 1)
+	@constraint(
+		model, [e in 1:E, s in 1:S],
+		y[e] + z[e, s] == sum(λ[t, s] for t in 1:T if columns[t][e])
+	)
+
+	optimize!(model)
+
+	return TwoStageSpanningTreeSolution(value.(y) .> 0.5, value.(z) .> 0.5)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Column generation heuristic, that solves the linear relaxation and then outputs the solution of the proble restricted to selected columns.
+Returns an heuristic solution.
+"""
+function column_heuristic(instance; MILP_solver=GLPK.Optimizer, verbose=true)
+	(; columns) = column_generation(instance; verbose=false, MILP_solver)
+	verbose && @info "Continuous relaxation solved with $(length(columns)) columns."
+	return column_heuristic(instance, columns)
+end