Mixed integer linear programming (MILP) approach to find an optimal investment order for cycles

[tsmbland]

Description

This PR introduces a mixed integer linear programming (MILP) approach to find an optimal investment order for cycles.

Just to give a bit of background/reminder on the current approach for finding the investment order in complete systems. The first step is to build a graph where the nodes are markets (i.e. a specific combination of commodity and region, such as the ELCTRI market in GBR), and the edges are processes. An edge exists from node A to node B if there exists a process that can convert A to B (currently these edges are all literal processes, but in the future we may introduce trade edges, e.g. COAL in GBR to COAL in FRA). In a graph with no cycles, we can find an investment order with no "conflicts" by doing a reverse topological sort. By "lack of conflicts" I mean that markets further along the supply chain will always be solved earlier.

If a system contains a cycle, then it's not possible to create an investment order with no conflicts - whatever order we come up with, there will have to be at least one pairwise violation, i.e. where market X is solved before market Y despite Y being a consumer of X. Such a conflict may then require us to redo investments for X after we've solved Y if the demand profile for X is significantly influenced by Y. Since this is extra work we'd like to minimise the number of times we potentially have to do this.

So far, I've written code to identify SCCs and compress them into "supernodes". The resulting condensed graph is guaranteed to have no cycles, so we can perform a topological sort as normal, where these "supernodes" as a whole will appear somewhere in the order. But the question remains as to what we do when a supernode is reached in the investment algorithm - currently we just exit the simulation.

An optimal approach is going to have two components:

• order the commodities in such a way that minimises the number of pairwise conflicts
• an approach that can deal with conflicts when they inevitably arise

I've tried a few approaches for the ordering problem, and have had the most success with a mixed integer linear programming (MILP) approach using highs, which is what I'm proposing in this PR (big credit to ChatGPT for suggesting this, but there seems to be a lot of literature on this approach as well)

The MILP approach builds a set of binary decision variables x_i,j over all combinations of i and j (nodes in the SCC), which is 1 if i comes before j in the ordering, or 0 otherwise. Constraints are build in such a way that these decision variables must be self-consistent:

• if x_i,j + x_j,i = 1 (i.e. if i comes before j then j cannot come before i)
• x_i,j + x_j,k + x_k,i <= 2 (i.e. if i comes before j and j comes before k, then k cannot come before i)

We apply costs to x_i,j based on the connections within the SCC:

• the cost of x_i,j = 1 if theres an outgoing edge from i to j (i.e. we penalise putting i before j in the ordering if j is a consumer of i)

Additionally, we want to encourage, nodes with connections going out of the SCC to come earlier in the ordering (the SCC needs to "see" any external demands early on, otherwise there will be no demand for investment):

• I've added a small bias to all x_i,j if i has outgoing edges to a node(s) outside the SCC (i.e. encouraging i to come before all other nodes in the SCC order)
• I've hardcoded the bias value to 0.1, but there are some implications of this: a small value means it will try to put exporter markets first, but it will never introduce conflicts for the sake of doing so, whereas a large (>1 value) means that it will introduce conflicts (if necessary) to put exporter markets first
• Not sure what's better really, something to keep in mind

The solver then optimises x_i,j to minimise penalty_i,j * x_i,j, with the constraints to ensure that x_i,j is self-consistent. We can then work out the optimal position of each i in the order by summing over all x_i,j (i.e. the number of times it appears before other nodes in all pairwise orderings). We then save this order within the InvestmentSet::Cycle, so that when the investment algorithm reaches the supernode, it will solve investments in this order (this will be another PR).

Some considerations:

• I haven't taken any consideration of flow coefficients. Intuitively, you'd might happier to break ranks for an X<-Y with a small coefficient (i.e. production of X requires a small amount of Y) than one with a large coefficient. But in practice that's not really feasible because of units (e.g. how do you compare consumption of wood in kg with consumption of electricity in PJ), and because each conversion may have multiple processes with different coefficients and you don't know which one will be picked
• Similarly, I don't take into account the amplitude of external demands. Again, intuitively you'd want to put markets with high external demand earlier, but there's the same problem around units, and we're also pre-computing the order at the start of the simulation when we won't know what the external demands will be
• The number of constraints is n^3. Supposedly the solver has a worst-case complexity of O(2^n), which could be very bad for large SCCs. I haven't tried this yet on anything bigger than ~5 nodes, but I think we'll just have to see what happens with more complex systems. I've also experimented a bit with a greedy algorithm, which may be much quicker for large SCCs, but it's been a bit fiddly and doesn't seem to behave as well. It may be that we end up providing a choice of algorithms

See also #999

Fixes # (issue)

Type of change

• Bug fix (non-breaking change to fix an issue)
• New feature (non-breaking change to add functionality)
• Refactoring (non-breaking, non-functional change to improve maintainability)
• Optimization (non-breaking change to speed up the code)
• Breaking change (whatever its nature)
• Documentation (improve or add documentation)

Key checklist

• All tests pass: $ cargo test
• The documentation builds and looks OK: $ cargo doc

Further checks

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works

---

Note

Introduce HiGHS-backed MILP to order nodes within SCCs, integrate into cycle compression, and add tests; make InvestmentSet hashable.

• Graph/Investment ordering:
  • Replace naive SCC handling with MILP-based internal ordering using highs to minimise forward edges within cycles, with a small external-export bias.
  • Integrate ordering into compress_cycles(&InvestmentGraph) and subsequent topological layering.
  • Adjust function signatures to pass graphs by reference and add logging for solver issues.
• Types:
  • Derive Eq and Hash for InvestmentSet to support lookup during SCC ordering.
• Tests:
  • Add test_order_sccs_simple_cycle validating MILP ordering on a toy cycle.

^{Written by Cursor Bugbot for commit a19a151. This will update automatically on new commits. Configure here.}

[codecov]

Codecov Report

❌ Patch coverage is 90.97744% with 12 lines in your changes missing coverage. Please review.
✅ Project coverage is 84.26%. Comparing base (c0e05fb) to head (6aaf758).
⚠️ Report is 15 commits behind head on main.

Files with missing lines	Patch %	Lines
src/graph/investment.rs	90.97%	6 Missing and 6 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1004      +/-   ##
==========================================
+ Coverage   84.12%   84.26%   +0.14%     
==========================================
  Files          52       52              
  Lines        5844     5974     +130     
  Branches     5844     5974     +130     
==========================================
+ Hits         4916     5034     +118     
- Misses        689      695       +6     
- Partials      239      245       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[cursor]

✅ Bugbot reviewed your changes and found no bugs!

---

[dalonsoa]

It clearly works - at least for the example given - and I think I follow the logic, but it takes quite a bit of effort to understand what's going on. I've left a couple of comments.

Maybe it will be useful to include explicitly in the docstring the matrix being solved, with all its constrains (maybe reduce the system to 4 nodes instead of 5, so you have less elements to worry about)

[dalonsoa]

I would include here a link to somewhere where this method is explained. Otherwise, the following discussion on antisymmetry, transitivity, etc. it is pretty hard to follow.

[dalonsoa]

This might need some inline comments, for clarity. Something like:

Suggested change

	for edge in original_graph.edges_directed(idx, Direction::Outgoing) {
	if let Some(&j) = index_position.get(&edge.target()) {
	penalties[i][j] = 1.0;
	} else {
	has_external_outgoing[i] = true;
	}
	}
	// We loop over the edges going out of this node
	for edge in original_graph.edges_directed(idx, Direction::Outgoing) {
	if let Some(&j) = index_position.get(&edge.target()) {
	// If they lead to another node within the SCC, we apply a penalty
	penalties[i][j] = 1.0;
	} else {
	// Otherwise they lead to an external node so we take note of it
	has_external_outgoing[i] = true;
	}
	}

[dalonsoa]

If I get this right, this enforces that 1 = x_ij + x_ij. However,

1. What enforces that x_ij and x_ji are not float numbers? x_ij = x_ji = 0.5 would be a valid solution, right?
2. What enforces that one of them is not negative, eg. x_ij = -5 and x_ji = 6`

In summary, it seems that the x are meant to be 0 or 1, but I do not see where that is enforced. You mention binary variables, but don't know where that happens.

[dalonsoa]

Oh, wait, is it here where you enforce that by using an integer column with value between 0 and 1... and so, only 0 and 1 are valid options?

[dalonsoa]

Very nice! I think this is way clearer now.