The Steiner tree problem

What is the shortest network of roads that joins a set of towns? The roads may meet at junctions of their own, Steiner points, and in the shortest network every Steiner point joins exactly three roads at angles of 120°. For the four corners of a square of side 1, the shortest network with no junctions but the corners, the minimum spanning tree, is 3 long, and the Steiner tree with two Steiner points is \(1 + \sqrt{3} \approx 2.732\) long. For many points, the problem is NP-hard, and the number of ways to connect them grows faster than exponentially. A soap film between two glass plates joined by pins forms a Steiner tree by itself, though not always the shortest one.

The example states the rule of 120° as constraints of a sketch, so that the solver places the Steiner points, and checks the trees against known lengths, for the Fermat point of a triangle, the two trees of a rectangle and the tree of a regular pentagon.

import math

import planex as px

A tree as a sketch

tree draws a line for every edge, from a rough position of each Steiner point. The towns are fixed. At each Steiner point, the ends of its three edges coincide, and two angles of 120° between its edges fix the third. The solver then moves the Steiner points to where the angles hold, and length adds up the lines:

def tree(towns, steiner, edges):
    """a sketch of the tree: towns {name: (x, y)} fixed, Steiner points {name: rough (x, y)}, edges [(a, b)]; and
    an end of a line at each Steiner point"""
    at = {**towns, **steiner}
    sk = px.Sketch()
    ends, lines = {s: [] for s in steiner}, {s: [] for s in steiner}
    for a, b in edges:
        if b in towns:
            a, b = b, a
        line = f"{a}_{b}"
        pen = sk.Move(at[a], px.Fix) if a in towns else sk.Move(at[a])
        pen.Line(at[b][0] - at[a][0], at[b][1] - at[a][1], id=line)
        for v, end in ((a, "start"), (b, "end")):
            if v in steiner:
                ends[v].append(f"{line}.{end}")
                lines[v].append(line)
    for s in steiner:
        first, *rest = ends[s]
        for other in rest:
            sk.Constrain(first, other, px.Coincident)
        e1, e2, e3 = lines[s]
        sk.Constrain(e1, e2, px.Angle, 120)
        sk.Constrain(e2, e3, px.Angle, 120)
    return sk, {s: ends[s][0] for s in steiner}


def length(sk):
    """the length of all lines of the sketch"""
    return sum(sk.Measure(line)["length"] for line in sk.elements["Line"])

Three towns

For a triangle whose angles are all below 120°, the one Steiner point is the Fermat point. Torricelli’s construction gives the length of the tree. With the equilateral triangle \(BCD\) erected on the side \(BC\), outside the triangle, the tree is as long as \(AD\), which the print compares (Figure 1):

A, B, C = (0, 0), (100, 0), (30, 70)
three, junctions = tree({"A": A, "B": B, "C": C}, {"S": (40, 25)}, [("A", "S"), ("B", "S"), ("C", "S")])
turn = math.radians(-60)
D = (B[0] + (C[0] - B[0]) * math.cos(turn) - (C[1] - B[1]) * math.sin(turn),
     B[1] + (C[0] - B[0]) * math.sin(turn) + (C[1] - B[1]) * math.cos(turn))  # C turned about B by -60°
print("status", three.status, "| tree", length(three), "| AD", math.dist(A, D))

status full | tree 157.87449335780033 | AD 157.87449335780033

Figure 1: Three towns, their Steiner point and Torricelli's construction, dashed

Four towns

The four corners of a rectangle can be joined in two ways with two Steiner points, by a bridge along the long side or by one across. The example is a rectangle of \(a = 150\) by \(b = 100\). Like a soap film, the solver finds the shortest tree of the topology it starts from, and in both topologies every angle of the tree is 120°. Their lengths are \(a + \sqrt{3}\,b\) and \(b + \sqrt{3}\,a\), and the bridge along the long side is shorter. Both are compared with these lengths (Figure 2):

rectangle = {"a": (0, 0), "b": (150, 0), "c": (150, 100), "d": (0, 100)}
along, _ = tree(rectangle, {"s": (40, 50), "t": (110, 50)}, [("a", "s"), ("d", "s"), ("s", "t"), ("b", "t"), ("c", "t")])
across, _ = tree(rectangle, {"s": (75, 30), "t": (75, 70)}, [("a", "s"), ("b", "s"), ("s", "t"), ("c", "t"), ("d", "t")])
print("along", length(along), "a + √3 b", 150 + math.sqrt(3) * 100)
print("across", length(across), "b + √3 a", 100 + math.sqrt(3) * 150)
print("minimum spanning tree", 150 + 2 * 100)

along 323.20508075688775 a + √3 b 323.2050807568877
across 359.8076211353316 b + √3 a 359.8076211353316
minimum spanning tree 350

The tree across the rectangle is even longer than the minimum spanning tree, three of the rectangle’s sides. A solver or a soap film that starts in the wrong topology stays in it.

(a) Along the long side, 323.2

(b) Across, 359.8

Figure 2: The two trees of a rectangle

Five towns

The corners of a regular pentagon need three Steiner points. Its Steiner tree is 3.8911568233 times its side long, a root of \(x^8 - 17x^6 + 29x^4 - 13x^2 + 1\), against 4 sides for the minimum spanning tree. For six towns at the corners of a regular hexagon, no Steiner point helps, and five of its sides are the shortest network, a result Du, Hwang and Weng proved for regular polygons in 1987. The pentagon’s tree for a side of 100 is compared with the closed form \(\sqrt{17 + 7\sqrt{5} + \sqrt{6\,(65 + 29\sqrt{5})}}\,/\,2\) (Figure 3):

R = 50 / math.sin(math.radians(36))  # the circumradius of a side of 100
corners = {f"t{k}": (R * math.cos(math.radians(90 + 72 * k)), R * math.sin(math.radians(90 + 72 * k))) for k in range(5)}


def toward(p, q, f):
    """the middle of p and q, pulled toward the center by f"""
    return (f * (p[0] + q[0]) / 2, f * (p[1] + q[1]) / 2)


starts = {"s1": toward(corners["t1"], corners["t2"], 0.7), "s2": (0, 0.1 * R), "s3": toward(corners["t3"], corners["t4"], 0.7)}
pentagon, junctions = tree(corners, starts, [("t1", "s1"), ("t2", "s1"), ("s1", "s2"), ("t0", "s2"), ("s2", "s3"),
                                             ("t3", "s3"), ("t4", "s3")])
exact = math.sqrt(17 + 7 * math.sqrt(5) + math.sqrt(6 * (65 + 29 * math.sqrt(5)))) / 2
print("status", pentagon.status, "| tree", length(pentagon), "| closed form", 100 * exact, "| spanning tree", 400)

status full | tree 389.1156823326853 | closed form 389.11568233268537 | spanning tree 400

Figure 3: The Steiner tree of a regular pentagon

See Mechanisms: the sketch as a solver for sketches that compute positions.