I am trying to solve a network design problem for gas pipeline networks (see PipeLayout). In the original problem statement, only the coordinates of the source and sink nodes are given. Pipes can be connected in any way, also using Steiner nodes at any location.

But I also consider a discretized variation of the problem, where the pipes are selected from a given ground structure. Initially, I thought of using a regular (square) grid as a ground structure, but this brings problems related to solution symmetry (among other things).

So, instead I will now try to find a (not so regular) triangulation of the area that starts with the given source & sink nodes, but may add additional points for refinement.

Mesh generation with triangles is a well-studied problems, with mature implementations (see Triangle), but is usually done in order to solve PDE systems or similar. So these tools define a good quality triangulation in terms of triangles with large angles or small areas.

For our application, the triangle cells themselves are not relevant, but rather the edges around them. These are used to connect the nodes, so they should provide good candidates for shortest paths and spanning trees. Therefore, we will evaluate several triangulations based on the minimum tree that spans all terminal nodes.

Implementation notes.¶

As usual, we will experiment with Julia code. The triangulation code is part of the PipeLayout package, and based on TriangleMesh. Visualizations are made with Makie and Cairo. Optimization models are formulated with JuMP and solved with SCIP.

In [1]:

# Activate custom environment
using Pkg
Pkg.activate("../../environments/refined-delaunay-for-flow-problems/");

# Julia standard library
using LinearAlgebra
using Random
using SparseArrays
using Statistics

# Third-party packages
using AbstractPlotting
using CairoMakie
using JuMP
using LightGraphs
using SCIP

# My own code (not registered)
using PipeLayout
using PipeLayout: delaunay_triangulation, make_pslg, TriangleSwitches, triangulate
using PipeLayout: antiparallel_digraph, euclidean_steiner_tree

Random Instance Data¶

We generate a single set of 7 nodes that is used throughout the experiment.

In [2]:

Random.seed!(0);
const N = 7
const WIDTH = 500
const HEIGHT = 300

x = round.(0.95 * WIDTH * rand(N))
y = round.(0.95 * HEIGHT * rand(N))
points = [x y];

Drawing with Makie¶

We define some drawing utitities to be used below.

In [3]:

function make_scene(width=WIDTH, height=HEIGHT)
    return Scene(resolution=(width, height), show_axis=false, scale_plot=false)
end

function draw_points!(points; markersize=4, color=:black)
    scatter!(points[:, 1], points[:, 2], markersize=markersize, color=color)
end

function draw_edges!(points, edges; color=:gray)
    linesegments!(points[edges'[:], :], color=color)
end

function draw_triangulation(points, edges)
    make_scene()
    draw_edges!(points, edges)
    draw_points!(points)   
end

draw_triangulation(trimesh) = draw_triangulation(trimesh.point', trimesh.edge');

Triangulation with TriangleMesh¶

TriangleMesh wraps the Triangle C library for Julia. They provide, in particular, an algorithm for Constrained Delaunay Triangulation. This takes as input a planar polygon and outputs a triangulation of the nodes, while ensuring that the existing edges are part of the triangulation (but may be subdivided).

In [4]:

# Create Delaunay Triangulation from points
deltri = delaunay_triangulation(points)

# Create Polygon from triangulation (as input for further refinements)
delpoly = make_pslg(deltri)

draw_triangulation(deltri)

Out[4]:

In [5]:

# Create triangulation of the raw points, specifying a minimum angle constraint
draw_triangulation(triangulate(points,
        TriangleSwitches(minimum_angle=20, conforming_delaunay=true)))

Out[5]:

In [6]:

# Create a CDT, keeping the edges of the original triangulation 
draw_triangulation(triangulate(delpoly,
        TriangleSwitches(minimum_angle=20, conforming_delaunay=true)))

Out[6]:

The simple Delaunay Triangulation makes sense visually and even has nice theoretical properties about approximating Euclidean distances with shortest paths!

The second triangle mesh introduces few additional nodes in order to get more regular triangles. But in turn, it also loses some of the "direct" connections between the terminal nodes.

Finally, the CDT keeps the direct edges of the Delaunay Triangulation, but refines the triangles, using many more additional nodes and edges.

Steiner Tree Model¶

We will now define an optimization model to determine the shortest subgraph that spans the given terminal nodes (Steiner tree in graphs). The model uses a multi-commodity flow formulation.

With the help of this model, we can see which edges of the triangulation are actually useful for our application, as well as quantitatively compare the triangulations.

In [7]:

function edge_lengths(points, edges)
    diff = points[edges[:, 1], :] - points[edges[:, 2], :]
    return dropdims(mapslices(norm, diff, dims=2), dims=2)
end

edge_lengths(trimesh) = edge_lengths(trimesh.point', trimesh.edge')

function edge_length_map(trimesh)
    edges = trimesh.edge'
    lengths = edge_lengths(trimesh)

    length_map = Dict{Tuple{Int64, Int64}, Float64}()
    for e = 1:size(edges, 1)
        s, t = edges[e, :]
        length_map[s, t] = lengths[e]
        length_map[t, s] = lengths[e]
    end
    return length_map
end

function steiner_tree(trimesh, terminals)
    # Compute length by Euclidean distance of nodes.
    lengths = edge_length_map(trimesh)
    
    # Build digraph with all antiparallel arcs, for nonnegative flow.
    graph = antiparallel_digraph(trimesh)
    nodes = collect(1:nv(graph))
    arcs = collect(keys(lengths))
    
    length(terminals) >= 2 || error("Need at least 2 terminals.")
    all(terminals .<= nv(graph)) || error("Terminals out of range.")
    root = terminals[1]
    sinks = terminals[2:end]
   
    demand(v, s) = 1.0*(v == s) - 1.0*(v == root)
    
    # Using arc length for fixed capacity cost and multi-commodity flow.
    model = JuMP.direct_model(SCIP.Optimizer(display_verblevel=0))
    @variable(model, select[a in arcs], Bin, container=SparseAxisArray)
    @variable(model, flow[a in arcs,s in sinks] ≥ 0, container=SparseAxisArray)
    @constraint(model, balance[v in nodes, s in sinks],
        sum(flow[(n, v), s] - flow[(v, n), s] for n in neighbors(graph, v))
        == demand(v, s))
    @constraint(model, capacity[a in arcs, s in sinks], flow[a, s] <= select[a])
    @objective(model, Min, sum(lengths[a] * select[a] for a in arcs))
    
    optimize!(model)
    
    # Return objective, as well as selected edges.
    return objective_value(model), value.(select)
end

Out[7]:

steiner_tree (generic function with 1 method)

Another utility will take a triangulation as input, compute the optimal subtree and visualize both.

In [8]:

function draw_tree(trimesh, terminals=1:N)
    points, edges = trimesh.point', trimesh.edge'
    obj, select = steiner_tree(trimesh, terminals)
    @show obj
    selected = [select[(s,t)] + select[(t,s)] for (s,t) in eachrow(edges)]
    active_edges = edges[selected .> 0.0, :]
    
    draw_triangulation(trimesh)
    linesegments!(points[active_edges'[:], :], color=:plum, linewidth=3)
    draw_points!(points[terminals, :], markersize=5, color=:teal)
end

Out[8]:

draw_tree (generic function with 2 methods)

In [9]:

draw_tree(deltri)

obj = 693.4867635250124

Out[9]:

In [10]:

tri = triangulate(points, TriangleSwitches(minimum_angle=20, conforming_delaunay=true))
draw_tree(tri)

obj = 700.5151721904267

Out[10]:

In [11]:

tri = triangulate(points, TriangleSwitches(maximum_area=sqrt(WIDTH*HEIGHT), conforming_delaunay=true))
draw_tree(tri)

obj = 697.572096498523

Out[11]:

In [12]:

tri = triangulate(delpoly, TriangleSwitches(maximum_area=sqrt(WIDTH*HEIGHT), conforming_delaunay=true))
draw_tree(tri)

obj = 693.4867635250124

Out[12]:

The Delaunay Triangulation allows for a solution of length 693.5. Other point cloud triangulations based on minimum angle or maximum area constraints yield poorer solutions, of length 700.5 and 687.5, while using a larger graph.

The CDT recovers the solution of the Delaunay Triangulation, but does not use any of the additional nodes or edges introduced in the subdivision.

Euclidean Steiner Tree¶

We can "cheat" and solve the EST problem directly (based on GeoSteiner). This will gives a lower bound on the length (what we can hope to achieve), and maybe a polygon to be triangulated.

In [13]:

est = euclidean_steiner_tree(points)

make_scene()
draw_edges!(est...)
draw_points!(est[1][1:N, :])

Out[13]:

We can see that 3 additional Steiner nodes are inserted.

Let's now build a triangulation of this tree and compute the length.

Then, we show a point cloud triangulation, based on the terminals and the Steiner nodes of the EST.

In [14]:

estpoly = make_pslg(est...)
draw_tree(triangulate(estpoly, TriangleSwitches(
            minimum_angle=20, maximum_area=sqrt(WIDTH*HEIGHT), convex_hull=true)))

obj = 664.018124937279

Out[14]:

In [15]:

draw_tree(triangulate(est[1], TriangleSwitches(minimum_angle=20, maximum_area=sqrt(WIDTH*HEIGHT))))

obj = 696.9337210443265

Out[15]:

The EST has a length of 664, which is significantly shorter than the 693.5 of the Delaunay Triangulation, so there is still some room for improvement.

Using only the additional Steiner nodes does not give a convincing solution.

Conclusion¶

Overall, the experiment is a little disappointing. Using any point cloud triangulation does not result in short trees, even if the triangulation becomes finer.

The Delaunay Triangulation is a strong candidate, given that it does not introduce any new nodes.

Improving on the Delaunay Triangulation with refinement based on CDT has not been fruitful on this example data, as the additional edges are never used for the tree.

Intuitively, to find short trees, one would like to connect the terminal nodes with direct, line-of-sight edges. But the irregular triangulation often runs orthogonal to that desired line, which would require frequent detours.

Finally, the actual solution of the EST can be used as the basis of a triangulation, which then contains the relevant edges and also looks regular and "natural". But from a methodology point of view, this means using the final result as input of the solution process.

Outlook¶

Luckily, the actual application has additional features, beyond the minimum length spanning tree. So the EST could still be used effectively in a preprocessing step. Maybe some pertubation of the node positions can be applied to yield triangles of more equal size, as in force-directed_graph_drawing.