Crossing Minimal Edge‐Constrained Layout Planning using Benders Decomposition

Abstract

We present a new crossing number problem, which we refer to as the edge‐constrained weighted two‐layer crossing number problem (ECW2CN). The ECW2CN arises in layout planning of hose coupling stations at BASF, where the challenge is to find a crossing minimal assignment of tube‐connected units to given positions on two opposing layers. This allows the use of robots in an effort to reduce the probability of operational disruptions and to increase human safety. Physical limitations imply maximal length and maximal curvature conditions on the tubes as well as spatial constraints imposed by the surrounding walls. This is the major difference of ECW2CN to all known variants of the crossing number problem. Such as many variants of the crossing number problem, ECW2CN is $\mathcal{NP}$ ‐hard. Because the optimization model grows fast with respect to the input data, we face out‐of‐memory errors for the monolithic model. Therefore, we develop two solution methods. In the first method, we tailor Benders decomposition toward the problem. The Benders subproblems are solved analytically and the Benders master problem is strengthened by additional cuts. Furthermore, we combine this Benders decomposition with ideas borrowed from fix‐and‐relax heuristics to design the Dynamic Fix‐and‐Relax Pump (DFRP). Based on an initial solution, DFRP improves successively feasible points by solving dynamically sampled smaller problems with Benders decomposition. Because the optimization model is a surrogate model for its time‐dependent formulation, we evaluate the obtained solutions for different choices of the objective function via a simulation model. All algorithms are implemented efficiently using advanced features of the GuRoBi‐Python API, such as callback functions and lazy constraints. We present a case study for BASF using real data and make the real‐world data openly available.

Introduction

BASF is the world's leading chemical company (Tullo 2018). It has customers in over 170 countries and supplies about 8.000 products to almost all industries. BASF operates more than 350 sites worldwide, six of them being major, so called “Verbund” sites. The largest site is in Ludwigshafen, Germany, where BASF has its headquarters. The Verbund system is one of the core competencies of BASF; it creates efficient value chains that extend from basic chemicals to high‐value‐added products such as coatings and crop protection agents. In addition, the by‐products of one plant can be used as the starting materials of another.

As part of its Verbund system, BASF operates more than 2.800 kilometers of pipelines in Ludwigshafen. While most of the pipes are dedicated to a certain product or utility, some of them are used for different products alternatingly. To maximize flexibility, BASF operates switching hubs, where dozens of pipes end in the same location, and any two of them might be connected on demand by flexible tubes. The work of connecting and disconnecting pipes was a manual and very unpopular task, due to the weight of the steel pipes and the need for full body safety equipment in the switching hub. However, necessary sub‐tasks, such as leaking test, flushing, and drying made it difficult to obtain an automated solution. Currently, the engineering department of BASF is designing a dual‐arm robotic system to accomplish that task. As the robot is mounted on the ceiling, with its two hands, it cannot remove a pipe crossed by another one. For that, the robot must either remove and rebuild all blocking connections, or release the blocked pipe to the ground, where it has to be removed manually later. Both possibilities impede the efficiency of the system. Thus, the Advanced Business Analytics department within BASF was asked for a mathematical solution of this problem, which translates into a crossing minimal design where the number of blocking connections (so‐called FIFO crossings) is to be reduced to a minimum.

To that end, we are given a two‐layer hose coupling station where units are positioned on two opposing sides surrounded by walls. These units are to be connected by tubes. Figure 1 depicts such a hose coupling station. The connecting tubes (edges) have to obey three physical restrictions (constraints). First, tubes must lie inside the box drawn by the units because the units abut walls. A close distance to the surrounding walls is needed since the fluid flows through units and the tubes emerge from taps integrated in the walls. Due to security reasons, these taps have to be close to the units without any disturbing elements in between, that is, moving the units away from the wall to gain more freedom for tube placement is not an option. Second, the tubes’ length has an upper bound which is imposed by the size of the tube storage. Third, the tubes can only be bent up to a maximal curvature; otherwise, they break. Moreover, the connections between the units change over time, that is, the tubes have to be reconfigured. Possible failure during the disconnection process may cause severe health risks due to leaking chemicals. Therefore, BASF is taking very high safety precautions, like full body protection. In the future, this is supposed to be executed by a dual‐arm robot. As mentioned above, for the robot, crossing tubes pose a challenge and may lead to additional efforts, as all tubes stacked above need to be disconnected first. Therefore, a layout plan for the units, which is crossing minimal with respect to the connecting tubes, is highly valuable.

Figure 1

An Illustration of a Possible Assignment of 9 Tube‐connected Units to 10 Positions. Due to the Maximal Curvature Condition, no Direct Neighbors may be Connected and a Maximal Tube Length has to be Respected. This is Illustrated for the Vertically Hatched Position (I) where all Reachable Positions are Horizontally Hatched (B, C, E, G, H). The Thickness of the Edges is Proportional to its Historical Connection Frequency [Color figure can be viewed at wileyonlinelibrary.com]

In the real‐world problem at BASF, the locations of 75 units have to be determined. The connections between the units change over time, and we recorded historical connection data over a time span of more than 3 years. Next, we performed a thorough data cleansing in close cooperation with the operating engineering department at BASF and dropped all historical connections which did not generalize well for future considerations. Eventually, we agreed on 375 representative connections whose weights were determined by their historical connection frequency. This leads to a graph with 75 nodes and 375 weighted edges, having an edge‐to‐node ratio of 5.

The available data do not allow a meaningful probabilistic evaluation of possible future connections, because the number of observations is very low compared to the huge number of possible pairings. Consequently, the data are too sparse to extract distributional information, let alone the fact that some sets of pairings exclude each other, while other sets are likely to co‐exist. Therefore, we do not set up a stochastic optimization model. Instead, we separate the problem into a deterministic design optimization problem, which serves as surrogate model, and a simulation problem, that is used for quality evaluation.

The deterministic design optimization problem leads to what we call the edge‐constrained weighted two‐layer crossing number problem (ECW2CN). We formulate this problem as a mixed integer linear program (MILP). The MILP finds an optimal assignment using aggregated data, that is, time independent data, considering both the number of connections between each pair of units and the aggregated connection frequencies. A natural objective function is the minimization of the number of reconfigurations over time which require the removal of FIFO connections. Because our proposed MILP is static, such an objective function cannot be modeled directly. To overcome this challenge, we set up multiple MILPs that differ in their treatments of the weighted edges by employing four different objective functions. Using these four different objective functions, we generate a series of feasible points for the MILP models. These points are then evaluated in a simulation model. This simulation model captures the complexity of the design problem by taking time dynamics and historical connections into account. In a second simulation model, only so‐called FIFO crossings are counted, that is, reconfigurations which require the removal of blocking connections which cause operational disruptions. The best performing point in the second simulation model is then chosen as the optimal design. Computing the FIFO crossings for all feasible points is not possible since the second simulation model is very time consuming.

The ECW2CN is a non‐standard edge‐weighted crossing minimization problem (Schaefer 2020). Because this problem has not been discussed in the literature before, we develop an optimization model and tailor solution techniques to solve this real‐world problem.

Our unique contributions are:

We introduce the edge‐constrained weighted two‐layer crossing number problem (ECW2CN). This is a new variant of the crossing number problem, motivated by a real‐world problem at BASF.

We analyze the computational complexity of ECW2CN and present MILP formulations, as well as two tailored decomposition algorithms: a Benders decomposition with direct cut calculation (BDC) and a Dynamic Fix‐and‐Relax Pump (DFRP). The master problem of BDC is strengthened by incorporating some information from relevant parts of the Benders subproblems. These additional cuts are provably stronger than the Benders optimality cuts. All cuts are computed analytically yielding a very efficient algorithm. The DFRP employs the BDC algorithm in a fix‐and‐relax fashion where the problem sizes are dynamically adjusted. This way, DFRP computes high quality solutions while being able to theoretically provide an optimality certificate.

We provide a case study using real data of BASF. The real‐world problem is solved by state‐of‐the‐art implementations of all algorithms. Furthermore, final solutions of the real‐world problem for different choices of the objective function are evaluated and compared in a time‐dependent simulation model.

The remainder is organized as follows. We place the ECW2CN problem in the literature in section 2. In section 3, we present a mathematical programming model for the ECW2CN problem and show that ECW2CN is

$\mathcal{NP}$

‐hard. This is followed by our two solution algorithms in section 4. The case study of the real‐world problem is presented in section 5, and we conclude with section 6.

Crossing Number Literature

The crossing number problem for a given graph G consists of finding a drawing of G, such that this drawing possesses the lowest number of edge crossings (Schaefer et al. 2008, Székely 2004). Many variants of graph crossing number problems have been studied during the last decades. The review paper by Schaefer (2020) lists about 86 different variants. Therefore, we focus our review on important aspects of crossing number problems that are relevant in our context.

The crossing number problem appeared in the 1940s in the context of railway crossing minimization (Turán 1977) and in the 1960s in circuit design (Sinden 1966). However, it took until 2005 for the first exact integer linear programming (ILP) formulations to appear (Buchheim et al. 2005), although, for the related maximum planar subgraph problem, ILP formulations were published a decade before (Jünger and Mutzel 1996a,b. For a general graph G(V, E), the ILP formulation by Buchheim et al. (2005, 2008) models edge crossing via binary decision variables x _e,f for all edges e ∈ E and f ∈ E. The difficulty lies then in ensuring that a graph exists for given values x _e,f . This is known as the realizability problem, which is itself an

$\mathcal{NP}$

‐complete problem. Buchheim et al. (2005, 2008) ensure this realizability by a family of (potentially exponential many) constraints. In 2008, Chimani et al. (2008) presented an improved ILP formulation, compared to the one presented in Buchheim et al. (2005, 2008), requiring one order of magnitude fewer binary variables by finding a compacter way to model the graph realizability. The three mentioned ILP formulations are solved by tailored branch‐and‐cut algorithms, which combine advanced prepossessing (Buchheim et al. 2008) with combinatorial column generation (Chimani et al. 2008). To avoid the difficult realizability problem, we exploit the special structure of two‐layer graphs present in ECW2CN (cf. section 3). This allows us to directly model the placement of the “nodes” (that is, units

$u\in \mathcal{U}$

on possible places

$p\in \mathcal{P}$

) via binary decision variable y _u,p . The crossing is then modeled by considering all quadruples

$p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }\in \mathcal{P}$

. This model allows us to consider various different ways of counting the edge crossings and to model the edge constraints regarding maximal length and maximal curvature explicitly.

In the literature, a variety of multi‐layer and two‐layer crossing problems are presented. Carpano (1980) discusses two types of hierarchical drawings (also called level drawings): horizontal drawings (nodes are placed on parallel lines) and radial drawings (nodes are placed on concentric circles) in a multi‐layer context. All edges have to be drawn as straight lines, in contrast to general drawing problems where edges can contain bends. The level specification of each node has to be respected. A semidefinite programming approach for hierarchical drawing is presented by Chimani et al. (2012). Hierarchical drawings are extended to allow for intra‐layer connections by Bachmaier et al. (2010); the intra‐layer edges are allowed to bend. With that, Bachmaier et al. (2010) comes close to the ECW2CN problem, when restricting to two layers, which is illustrated in Figure 2. In the ECW2CN problem, we have also a two‐layer set‐up with inter‐layer connections, but the units can be placed freely among the two sides in contrast to Bachmaier et al. (2010). Another similarity of the ECW2CN problem to Bachmaier et al. (2010) is that all edges need to stay inside the rectangle spanned by the location positions. Different to Bachmaier et al. (2010), the ECW2CN problem has restrictions on the inter‐ and intra‐layer connections.

Figure 2

A Comparison of ECW2CN with Related Settings in the Literature, Yielding to Different Number of Minimal Crossings

A special case of two‐layer hierarchical drawings is obtained for bipartite graphs; in the context of multigraphs, Garey and Johnson (1983) call this problem the bipartite crossing number problem. Because only straight lines are allowed for the edges, Jünger and Mutzel (1996a, b) refer to this problem as the two‐layer straight line crossing problem. This is also the usual setting for the k‐layer case where only straight lines are allowed. In the two‐layer straight line crossing problem, the bipartition of the node set has to be respected (hence, there are no intra‐layer edges). An ordering of the nodes, for each of the two layers, determines the graph and its crossings. This allows for efficient permutation approaches (Kobayashi et al. 2014, Laguna and Marti 1999, Palubeckis et al. 2019). As such, the bipartite crossing number problem is highly related to the (quadratic) linear ordering problem (Buchheim et al. 2010, Shahrokhi et al. 2001). Edge‐weighted variants of the bipartite crossing number problem are presented by Çakiro

$\overline{\mathrm{g}}$

lu et al. (2009) and Kobayashi et al. (2014), where the weights are handled in a multiplicative manner. The ECW2CN problem also possess edge weights.

In sum, ECW2CN has the following key specifications:

Maximal length and maximal curvature conditions for the tube (that is, edges) have to be met.

Though we have a two‐layer set‐up, intra‐layer connections are allowed; like the formulation in Bachmaier et al. (2010).

All connections have to stay inside a rectangle spanned by the position locations; similar to the bipartite crossing number problem (Garey and Johnson 1983).

Each node (that is, unit u ∈ U) can be placed freely among the two layers; thus, this problem differs from bipartite or two‐layer crossing minimization problems and yields some additional overlapping types, see Figure 3.

The edges are weighted due to the historical connection frequencies. However, in contrast to the set up of Kobayashi et al. (2014), weights can be handled in a very flexible framework. In this work, we examine four different manners.

The corresponding graph does not have to be connected. Connectivity is a standard assumption as the crossing number problems decompose with their connected components. This is not true for ECW2CN because of the maximal length and maximal curvature conditions.

The number of units and places do not have to be the same, that is, there might exist isolated nodes in the corresponding graph. This is a special case of a disconnected graph (see previous point). The “placement” of the isolated nodes is non‐trivial due to the maximal length and maximal curvature conditions.

We are facing an edge‐to‐node ratio of 5. Most instances found in the literature have a ratio between 1 and 2.

Figure 3

An Overview Over all Different Crossing Types in the ECW2CN

ECW2CN, its Complexity and MILP Formulations

Based on the definition of the ECW2CN in section 3.1, we prove its

$\mathcal{NP}$

‐hardness (section 3.2) and present our MILP models (section 3.3).

The Edge‐Constrained Weighted two‐Layer Crossing Number Problem (ECW2CN)

We begin by defining the notation used in the definition of the ECW2CN as well as in the MILP.

Indices, sets and parameters

$u\in \mathcal{U}$

: Units

$p\in \mathcal{P}$

: Ordered set of positions; each position p has an (X _p , Y _p )‐coordinate in the plane and all positions are aligned on two layers, that is, there are only two different values for X _p

D > 0: Maximal Euclidean distance between two units due to a maximal tube length

E > 0: Minimal distance between two units on the same layer due to maximal curvature

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$

: Four tuples of positions, such that connections of units on positions (p, p′) and

$(\widehat{p},{\widehat{p}}^{\prime })$

respect the maximal tube length as well as the maximal curvature condition parameterized by D and E and induce one of the crossing types illustrated in Figure 3.

$F_{u,u^{\prime }}\ge 0$

: Aggregated historical connection frequency of two units

$u,u^{\prime }\in \mathcal{U}$

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\ge 0$

: A penalty term that is used to penalize crossings weighted by their frequencies; therefore,

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=0$

if and only if

$F_{u,u^{\prime }}=0$

and/or

$F_{\widehat{u},{\widehat{u}}^{\prime }}=0$

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}$

: Four tuples of pairwise different connected units, that is,

$F_{u,u^{\prime }}>0$

and

$F_{\widehat{u},{\widehat{u}}^{\prime }}>0$

$(p,p^{\prime })\in \mathcal{N}$

: Pairs (p, p′) that violate the maximal length or maximal curvature condition

We are now ready to formally define the ECW2CN.

Definition 1 Edge‐constrained weighted two‐layer crossing number problem

Given is a set of units

$\mathcal{U}$

, positions

$\mathcal{P}$

on two layers with coordinates (X _p , Y _p ) for each position

$p\in \mathcal{P}$

, maximal Euclidean distance D between two units, minimal Euclidean distance E > 0 between two units on the same layer, edge weights

$F_{u,u^{\prime }}\ge 0$

for pairs of units

$u,u^{\prime }\in \mathcal{U}$

and a function

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\ge 0$

of edge weights for

$u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }\in \mathcal{U}$

. Then, the edge‐constrained weighted two‐layer crossing number problem (ECW2CN) assigns each unit

$u\in \mathcal{U}$

to at most one position

$p\in \mathcal{P}$

, such that the Euclidean distance of any pair of connected units is not greater than D, at least E if both units are on the same layer and all edges stay within the rectangle spanned by the convex hull of the position candidates, while minimizing the sum over

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}$

for all units

$u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }\in \mathcal{U}$

whose placement defines a crossing.

ECW2CN is

$\mathcal{NP}$

‐Hard

To show that ECW2CN is

$\mathcal{NP}$

‐hard, we first define a decision version of ECW2CN, denoted by ECW2CN‐D.

Definition 2 ECW2CN decision version

Instance: A set of units

$\mathcal{U}$

, positions

$\mathcal{P}$

on two layers with coordinates (X _p ,Y _p ) for each position

$p\in \mathcal{P}$

, maximal Euclidean distance D between two units, minimal Euclidean distance E > 0 between two units on the same layer, edge weights

$F_{u,u^{\prime }}\ge 0$

for pairs of units

$u,u^{\prime }\in \mathcal{U}$

, a function

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\ge 0$

of edge weights for

$u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }\in \mathcal{U}$

and bound C < ∞.

Question: Is there an assignment of each unit

$u\in \mathcal{U}$

to at most one position

$p\in \mathcal{P}$

, such that the Euclidean distance of any pair of connected units is not greater than D, at least E if both units are on the same layer, all edges stay within the rectangle drawn by the position locations, and such that the sum over

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}$

for all units

$u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }\in \mathcal{U}$

whose placement defines a crossing does not exceed C?

In our reduction, we use the decision version of the bipartite crossing number (BCN) problem, denoted by BCN‐D.

Definition 3 BCN decision version, (BCN decision version, Garey and Johnson (1983)

Instance: A connected bipartite multigraph G(V ₁,V ₂,E) and an integer K< ∞.

Question: Can G be embedded in a unit square so that all vertices of V ₁ are on the northern boundary, all vertices in V ₂ are on the southern boundary, all edges are within the square and there are at most K crossings?

We need the following complexity result.

Lemma 1 Garey and Johnson (1983)

The BCN‐D is

$\mathcal{NP}$

‐complete.

For their reduction, Garey and Johnson (1983) utilize the Optimal Linear Arrangement problem which has itself been reduced from the MAX‐Cut problem in (Garey et al. 1974). Therefore, BCN‐D is

$\mathcal{NP}$

‐complete in the strong sense.

This allows us to state and prove our main result.

Theorem 1

The decision version ECW2CN is

$\mathcal{NP}$

‐complete in the strong sense.

Proof 1

Observe that ECW2CN‐D is in

$\mathcal{NP}$

.

Given an instance of the BCN‐D, that is, a connected bipartite multigraph G(V ₁,V ₂,E), we construct an instance of the ECW2CN‐D as follows:

$\mathcal{U}=V_1\cup V_2$

$\mathcal{P}=\mathcal{U}$

with

$(X_i,Y_i)=\left(0,\frac{i-1}{|V_1|-1}\right)$

for

$i=1,\dots ,|V_1|$

and

$(X_i,Y_i)=\left(1,\frac{i-1}{|V_2|-1}\right)$

for

$i=1,\dots ,|V_2|$

D = E = 2

$F_{u,u^{\prime }}=0$

for all

$u,u^{\prime }\in \mathcal{U}$

;

$F_{u,u^{\prime }}\leftarrow +1$

for all

$(u,u^{\prime })\in E$

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=F_{u,u^{\prime }}·F_{\widehat{u},{\widehat{u}}^{\prime }}$

for all

$F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}>0$

; o/w 0.

The idea is to use the restriction on the minimal distance between two units on the same layer to ensure that only nodes from V ₁ (and V ₂) are assigned to the same layer. This can be achieved by assigning any number E > 1 (see below). The transformation is shown in Figure 4.

Figure 4

A Transformation of a “YES” Instance for K ≥ 5 of the BCN‐D to a “YES” Instance for C ≥ 5 of the ECW2CN‐D Using the Technique Described in the Proof of Theorem 1

The maximal distance is no restriction as long as

$\bar{D}\ge \sqrt{2}$

. This construction is (strongly) polynomial in the size of G.

Next, we show that an instance for BCN‐D with bound K is a “YES” instance, if and only if the transformed graph is a “YES” instance for ECW2CN‐D with bound C = K. For that, we observe that both problems count the number of crossings. Thus, if the solution of one problem maps to a solution of the other, then the bounds are the same. It remains to show the existence of such a mapping.

“⇒” Any feasible point for ECW2CN‐D yields a partition of

$\mathcal{U}$

into two disjoint sets

${\bar{V}}_1$

and

${\bar{V}}_2$

where

${\bar{V}}_1=V_1$

or

${\bar{V}}_1=V_2$

. To see this, assume that

$v_1\in V_1$

and

$v_2\in V_2$

are both assigned to

${\bar{V}}_1$

. If

$(v_1,v_2)\in E$

, then v ₁ and v ₂ do not respect E. If (v ₁,v ₂) ∉ E, then there exits a path from v ₁ to v ₂ in G, because G is connected. Choose any such path p. Then, there exists a node

$v_3\in V_1$

on p with (v ₂, v ₃) ∉ E. If

$v_3\in {\bar{V}}_1$

, then v ₂ and v ₃ do not respect E; otherwise,

$v_3\in {\bar{V}}_2$

and (v ₁, v ₃) ∉ E. Continuing this argument along path p from v ₃ to v ₁ shows the contradiction.

“⇐” Any feasible point for BCN‐D respects D and E. With that, it is a feasible point for ECW2CN‐D.

The Mixed Integer Linear Programming Models

We use the indices, sets, and parameters as introduced in section 3.1. Furthermore, we require the following decision variables

$y_{u,p}\in \{0,1\}$

: Binary decision variable which is 1 if unit u is located at position p, 0 otherwise

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge 0$

: Continuous decision variable that represents the maximal crossing frequency of arbitrary units located on the positions

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })$

.

The ECW2CN can be modeled as the following MILP formulation, which we call the monolithic model.

$$\begin{array}{cc}& (\mathcal{M}):z^{*}:=\underset{y,o}{\text{min}}\sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}{o}_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\hfill \\ & \text{s.t.}(y_{u,p}+y_{u^{\prime },p^{\prime }}+y_{\widehat{u},\widehat{p}}+y_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\le o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\hfill \end{array}$$

1

$$\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}\forall (u,u^{\prime },\widehat{u},\widehat{u})\in \mathcal{Q}$$

2

$$\sum _{p\in \mathcal{P}}{y}_{u,p}=1\forall u\in \mathcal{U}$$

3

$$\sum _{u\in \mathcal{U}}{y}_{u,p}\le 1\forall p\in \mathcal{P}$$

4

$$y_{u,p}+\sum _{u^{\prime }:F_{u,u^{\prime }}>0}{y}_{u^{\prime },p^{\prime }}\le 1\forall (p,p^{\prime })\in \mathcal{N}\forall u\in \mathcal{U}$$

5

$$y_{u,p}\in \{0,1\}\forall u\in \mathcal{U}\forall p\in \mathcal{P}$$

6

$$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge 0\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$$

7

The objective function (1) minimizes the weighted number of crossings. This weighted crossing is modeled through decision variables

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}$

. Constraints (2) together with the nonnegativity of

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}$

ensures that whenever a crossing is induced by locating two pairs of connected units

$u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }$

on positions

$p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }$

, we have

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}=P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}$

, that is, the crossing is weighted with

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}>0$

. Constraint blocks (3)–(4) are assignment restrictions. Each unit

$u\in \mathcal{U}$

has to be put at exactly one place

$p\in \mathcal{P}$

as ensured by constraints (3). That at most one unit

$u\in \mathcal{U}$

is located at each place

$p\in \mathcal{P}$

is given by the SOS‐1 constraints (4). Constraints (5) ensure the maximal length and maximal curvature restrictions. We discuss these constraints now in detail.

Constraints (5) can be more naturally written as

$$y_{u,p}+y_{u^{\prime },p^{\prime }}\le 1\forall (p,p^{\prime })\in \mathcal{N}\forall u,u^{\prime }\in \mathcal{U}:F_{u,u^{\prime }}>0.$$

8

These constraints (8) define a so‐called conflict graph, where y _u,p and

$y_{u^{\prime },p^{\prime }}$

cannot both equal 1. Finding the maximal number of possible y _u,p  = 1 values leads to the so‐called stable set problem which is also known as the independent set problem. The stable set problem is related to the maximum clique problem through the inverse graphs. Both problems are classical and well‐studied combinatorial optimization problems (Bomze 1997, Bomze et al. 1999, Rebennack et al. 2011). They are both

$\mathcal{NP}$

‐hard (Garey and Johnson 1979). Constraints (8), together with nonnegativity of y _u,p , define the so‐called stable set polytope. Therefore, we can make use of the knowledge about the stable set polytope in our context. Specifically, constraints (8) are facet defining for the stable set polytope, if and only if the edge in the conflict graph defines as maximal clique (Rebennack et al. 2012). This inspired constraints (5), which are a strengthened version of (8) as can be observed by adding the corresponding

$F_u=|\{u^{\prime }:F_{u,u^{\prime }}>0\}|$

constraints to obtain

$$F_u{y}_{u,p}+\sum _{u^{\prime }:F_{u,u^{\prime }}>0}{y}_{u^{\prime },p^{\prime }}\le F_u\forall (p,p^{\prime })\in \mathcal{N}\forall u\in \mathcal{U}.$$

Also note that there are

$|\mathcal{N}|{\sum }_{u\in \mathcal{U}}{F}_u$

constraints with two non‐zeros each, while there are only

$|\mathcal{N}||\mathcal{U}|$

constraints (5) with F _u  + 1 non‐zeros each. Therefore, (5) yield to fewer non‐zeros, to less constraints and to stronger constraints compared to (8).

Nevertheless, constraints (5) are not facet defining for polytope (3)‐(5) as they are dominated by

$$\sum _{u^{\prime }\in \mathcal{U}}{y}_{u^{\prime },p}+\sum _{u^{\prime }:F_{u,u^{\prime }}>0}{y}_{u^{\prime },p^{\prime }}\le 1\forall (p,p^{\prime })\in \mathcal{N}\forall u\in \mathcal{U}.$$

9

Because (9) are much denser than (5) and the maximal violation of constraints (9) is at most 1 for any point of the LP relaxation of (3)‐(5), we consider (5) instead.

From the above discussions, we summarize

Proposition 1

Formulation

$(\mathcal{M})$

models ECW2CN correctly.

We close this section with a few remarks.

Remark 1

The subsequent theoretical developments do not the depend on a particular choice of the crossing weights, that is, they hold true for all

$P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}:{\mathbb{R}}^{>0}\times {\mathbb{R}}^{>0}\to {\mathbb{R}}^{>0}$

. However, for the problem at hand we consider the following four candidates:

1. $P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=F_{u,u^{\prime }}+F_{\widehat{u},{\widehat{u}}^{\prime }}$
   (“Sum”)

2. $P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=\min (F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }})$
   (“Minimum”)

3. $P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=F_{u,u^{\prime }}·F_{\widehat{u},{\widehat{u}}^{\prime }}$
   (“Product”)

4. $P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}=1$
   (“One”)

For illustrative purposes, suppose we have a crossing induced by the positions of (u, u ^′) and

$(\widehat{u},{\widehat{u}}^{\prime })$

. The durations of these connections are given by

$F_{u,u^{\prime }}$

and

$F_{\widehat{u},{\widehat{u}}^{\prime }}$

. Then it is intuitive to see that the probability of having an actual crossing at a point in time increases for larger values of

$F_{u,u^{\prime }}$

and

$F_{\widehat{u},{\widehat{u}}^{\prime }}$

. That motivates the “Sum” and the “Product” objective function for the time static MILP formulation. The second choice “Minimum” provides an upper bound on the expected number of frequencies since a crossing cannot last longer than

$\text{min}(F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }})$

. The last candidate counts the unweighted number of crossings for benchmark purposes. These different variations are illustrated in Figure 5 in Example 1. In section 5, we show that there is indeed a high correlation between good feasible points with respect to the objective functions above and actual crossings which are computed by evaluating the time‐dependent simulation model.

Figure 5

Optimal Configurations of the Graph in Figure 1 with Respect to Four Different Objective Functions

Example 1

The graph depicted in Figure 1 possesses nine nodes and nine edges whose weights are given by: (A, F): 10, (A, G): 8, (B, D): 6, (B, F): 2, (C, E): 3, (C, F): 5, (D, H): 1, (E, I): 15, (H, I): 12. Depending on the choice of the objective function, we obtain the following optimal assignments as illustrated in Figure 5.

Remark 2

In order to break some symmetry, we extend formulation

$(\mathcal{M})$

by the constraint that the unit with the most connections must be located in the upper left half of the rectangle spanned by the positions.

Remark 3

Observe that it is possible to reformulate the monolithic model

$(\mathcal{M})$

as a minimax model of the form

$$\begin{array}{cc}& \underset{y}{\text{min}}\sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}\underset{(u,u^{\prime },\widehat{u},\widehat{u})\in \mathcal{Q}}{\text{max}}(y_{u,p}+y_{u^{\prime },p^{\prime }}+y_{\widehat{u},\widehat{p}}+y_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\hfill \\ & \text{s.t.}(3)-(6)\hfill \end{array}$$

which reveals the intrinsic piecewise linear structure of its objective function.

Remark 4

For our real‐world application, we have

$|\mathcal{C}|=532,401$

and

$|\mathcal{Q}|=129,922$

, that is, Equation (2) translates into 69,170,602,722 constraints. Even to build the mathematical model

$(\mathcal{M})$

on a personal computer with default settings is a challenge without running into memory issues. Due to its size, it is currently not possible to solve the monolithic model

$(\mathcal{M})$

directly using state‐of‐the‐art MILP solvers. Therefore, in section 4, we develop tailored solution approaches to tackle this problem.

Solution Methods

The monolithic formulation

$(\mathcal{M})$

is a MILP problem of large scale, due to constraints (2) and variables

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}$

. We propose a tailored Benders decomposition to exploit its special structure (section 4.1). This Benders decomposition is then used within a fix‐and‐relax‐inspired algorithm (section 4.2) in an effort to compute good feasible points for the real‐world problem at hand.

A Tailored Benders Decomposition with Direct Cut Calculation (BDC)

Benders Decomposition

The idea of Benders (1962) was to separate a, what he called, mixed variables problem into a master problem and a subproblem. The mixed variables in our case are the y _u,p and the

$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}$

variables. The master problem is a relaxation of the original problem and does not contain the so‐called complicating constraints (in our case, constraints (2)), which connect both types of variables.

These constraints are passed onto the subproblem which computes, for a given solution of the master problem (in our case, trial values of variables y _u,p ), the objective function value associated with an optimal “response” to y _u,p . Exploiting dual information, a feasibility or optimality cut is generated and passed on to the master problem. This way, the master problem is extended by information from the subproblem, in terms of dual extreme directions or extreme points. For a minimization problem, the master problem yields a lower bound and the trial solution together with the subproblem yields an upper bound. By iteratively solving master and subproblems until the upper and lower bounds are sufficiently close to each other, the original problem is solved. This so‐called Benders algorithm always terminates after finitely many iterations because the dual subproblem contains only a finite number of extreme directions and extreme points. For further details, we refer to the literature (Rahmaniani et al. 2017, Rebennack 2016).

For a given trial solution

${\widehat{y}}_{u,p}$

of the master problem, the subproblem

$(\mathcal{S})$

is given by

$$\begin{array}{cc}& z_S^{*}(\widehat{y}):=\text{min}\sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}{o}_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\hfill \\ & \text{s.t.}{o}_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge ({\widehat{y}}_{u,p}+{\widehat{y}}_{u^{\prime },p^{\prime }}+{\widehat{y}}_{\widehat{u},\widehat{p}}+{\widehat{y}}_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\hfill \end{array}$$

10

$$\forall (u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q},(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$$

11

$$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge 0\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$$

12

Notice that

$(\mathcal{S})$

is a linear programming (LP) problem which is feasible for any trial values

${\widehat{y}}_{u,p}\in [0,1]$

. Consequently, we do not need feasibility cuts.

Remark 5

Observe that for given

$\widehat{y}$

the unique solution of

$(\mathcal{S})$

is given by

$$o_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}^{*}=\underset{(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}}{\text{max}}\left\{0,({\widehat{y}}_{u,p}+{\widehat{y}}_{u^{\prime },p^{\prime }}+{\widehat{y}}_{\widehat{u},\widehat{p}}+{\widehat{y}}_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\right\}\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$$

which, in particular, implies its non‐degeneracy. This rules out the possibilities to add further cuts of different strengths from degenerate solutions as elaborated by Magnanti and Wong (1981), Papadakos (2008) and Sherali and Lunday (2013).

Let

${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}$

be an optimal dual (basic) solution associated with constraints (11). Because of strong duality, we obtain

$$z_S^{*}(\widehat{y})=z_D^{*}(\widehat{y}):=\sum _{(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}}\sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}\left(({\widehat{y}}_{u,p}+{\widehat{y}}_{u^{\prime },p^{\prime }}+{\widehat{y}}_{\widehat{u},\widehat{p}}+{\widehat{y}}_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\right){\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}.$$

For some (cut) index c, we define the (variable) cut coefficient

$${\pi }_{u,p,c}^{\mathrm{var}}:=\sum _{u^{\prime },\widehat{u},{\widehat{u}}^{\prime }:(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}}\sum _{p^{\prime },\widehat{p},{\widehat{p}}^{\prime }:(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}{\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}\forall u\in \mathcal{U},p\in \mathcal{P}$$

13

and cut constant

$${\pi }_c^{\mathrm{con}}:=-3\sum _{(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}}\sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}}{P}_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}{\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}.$$

Because the dual subproblem is a maximization problem, we obtain

$$z_S^{*}(y)\ge \sum _{u\in \mathcal{U},p\in \mathcal{P}}{\pi }_{u,p,c}^{\mathrm{var}}{y}_{u,p}+{\pi }_c^{\mathrm{con}}\forall y\in {\{0,1\}}^{|\mathcal{U}|\times |\mathcal{P}|}.$$

This yields the optimality cut

$$\eta \ge \sum _{u\in \mathcal{U},p\in \mathcal{P}}{\pi }_{u,p,c}^{\mathrm{var}}{y}_{u,p}+{\pi }_c^{\mathrm{con}}$$

for the master problem with free variable η.

The master problem

$(\mathcal{MP})$

can then be summarized as

$$z_M^{*}:=\text{min}\eta$$

14

$$\text{s.t.}(3)-(6)$$

15

$$\eta \ge \sum _{u\in \mathcal{U},p\in \mathcal{P}}{\pi }_{u,p,c}^{\mathrm{var}}{y}_{u,p}+{\pi }_c^{\mathrm{con}}\forall c\in \mathcal{O}$$

16

$$\eta \ge 0$$

17

where set

$\mathcal{O}$

is the cut set of Benders optimality cuts. The master problem is a MILP.

For a given cut set

$\mathcal{O}$

, the master problem yields a lower bound on the optimal solution value, that is,

$z_M^{*}\le z^{*}$

. In our case, the subproblem yields an upper bound for any trial

$\widehat{y}$

, that is,

$z_S^{*}(\widehat{y})\ge z^{*}$

. If

$$z_S^{*}(\widehat{y})-z_M^{*}\le \epsilon ,$$

then an ɛ‐optimal solution to

$(\mathcal{M})$

has been computed; for some ɛ > 0. Otherwise, an optimality cut is generated for trial

$\widehat{y}$

and the master problem is re‐solved.

Direct Cut Calculation

The idea of Benders decomposition is to exploit the special structure of the original problem

$(\mathcal{M})$

. Sometimes, the subproblem also exhibits a special structure which may allow a more efficient solution than solving it as an LP. For example, the stochastic minimum s − t cut problem features a maximum flow subproblem (Rebennack et al. 2020). Similarly, the nonlinear (but convex) subproblem of the power system expansion problem in Lohmann and Rebennack (2017) is a market clearing mechanism which can be solved efficiently by a tailored sorting algorithm. Another example arises in a stochastic production‐inventory planning problem where the dual subproblems are solved using some ordering property (Golari et al. 2017). An entire class of Benders algorithms exploiting the special structure of the subproblems is logic‐based Benders decomposition (Naderi et al. 2021). In our case, the subproblem

$(\mathcal{S})$

, and its dual, can also be solved very efficiently.

For subproblem

$(\mathcal{S})$

, it suffices to obtain an optimal dual (basic) solution

${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}$

. With such a dual point, we can calculate the cut coefficient

${\pi }_{u,p,c}^{\mathrm{var}}$

, the cut constant

${\pi }_c^{\mathrm{con}}$

and the optimal objective function value

$z_S^{*}(\widehat{y})$

. Specifically, we do not require the solution of the primal subproblem

$(\mathcal{S})$

.

The dual subproblem is a continuous knapsack problem of special form. For trial

$\widehat{y}$

, observe that an optimal dual (basic) solution is given by

$${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}=\left\{\begin{array}{cc}1,\hfill & \text{if}(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })={\text{argmax}}_{(\widehat{u},{\widehat{u}}^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}}\left\{({\widehat{y}}_{\widehat{u},p}+{\widehat{y}}_{{\widehat{u}}^{\prime },p^{\prime }}+{\widehat{y}}_{\widehat{u},\widehat{p}}+{\widehat{y}}_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{\widehat{u},{\widehat{u}}^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}\right\}>0\hfill \\ 0,\hfill & \text{o/w}\hfill \end{array}\right.$$

18

for all

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}$

and

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$

. The maximum is positive, if and only if there is a tuple

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}$

of units which are placed on positions

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$

.

Assume that

$\widehat{y}$

is binary. Because at most one unit u can be placed at any position p, this tuple

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

must be unique, if it exists. Consequently, for binary values of

$\widehat{y}$

, argmax is unique (implying that the solution is basic) and can be computed by inspecting all tuples

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

. Thus,

${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}=1$

if and only if units

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}$

cause an overlap for positions

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$

. These overlaps can be computed efficiently by inspecting whether the tuples

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

cause an overlap. This yields a runtime of

$O(|\mathcal{Q}|)$

to identify the 1 entries for

${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },·,·,·,·}^{*}$

. The cut coefficients

${\pi }_{u,p,c}^{\mathrm{var}}$

can also be calculated in

$O(|\mathcal{Q}|)$

through (13). This is extremely efficient when observing that subproblem

$(\mathcal{S})$

contains

$|\mathcal{Q}|·|\mathcal{C}|$

functional constraints.

In case

$\widehat{y}$

is fractional, the argmax in (18) may contain more than a single value  > 0. Again, it suffices to run through all tuples

$(u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}$

, retrieve the corresponding position tuples

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })$

, check if they induce a crossing (that is, check if

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}$

) and save them if necessary for later computation of

$$(y_{u,p}+y_{u^{\prime },p^{\prime }}+y_{\widehat{u},\widehat{p}}+y_{{\widehat{u}}^{\prime },{\widehat{p}}^{\prime }}-3)·P_{F_{u,u^{\prime }},F_{\widehat{u},{\widehat{u}}^{\prime }}}$$

and to choose the maximum value among them.

Benders Based Branch‐and‐Cut

Next to the efficient calculation of the cut coefficients, the particular implementation of the Benders decomposition can have a significant effect on the algorithm's computational performance. We chose a state‐of‐the‐art implementation based on lazy constraints using callback functions. For that, observe that a valid Benders optimality cut is associated with each trial

$\widehat{y}$

, not necessarily being binary. This allows us to compute and add Benders optimality cut(s) each time a trial has been computed in the branch‐and‐bound tree. The resulting algorithm is also known in the literature as “Benders based branch‐and‐cut” (Naoum‐Sawaya and Elhedhli 2013).

We follow a strategy similar to Adulyasak et al. (2015) and add Benders cuts in two cases: (1) at the root node of the branch‐and‐bound tree, these are global cuts, and (2) for feasible (incumbent) solutions, these are local cuts. This strategy avoids computing an excessively large number of Benders cuts. Specifically, we add Benders optimality cuts in the root node until

$$z_S^{*}(\widehat{y})-z_{M,\mathrm{LP}}^{*}\le \epsilon$$

for

$\widehat{y}$

of the master problem, where

$z_{M,\mathrm{LP}}^{*}$

is the optimal objective function value of the LP relaxation of the master problem

$(\mathcal{MP})$

. This ensures that

$$z_{M,\mathrm{LP}}^{*}-z_{\mathrm{LP}}^{*}\le \epsilon ,$$

19

with optimal objective function value

$z_{\mathrm{LP}}^{*}$

of the LP relaxation of the monolith

$(\mathcal{M})$

. Condition (19) certifies that we have solved the LP relaxation of the monolith with the Benders decomposition approach before the branching starts. The solution of the LP relaxation of

$(\mathcal{M})$

via Benders decomposition is efficient because the master problem works only with the y _u,p variables, avoiding the excessively large number of variables and constraints in (2). Additional efficiency is gained by using the direct cut calculation via (18) instead of solving LPs in order to compute the optimality cuts.

Strengthening the Master Problem

Inspired by Rahmaniani et al. (2017), we strengthen the master problem (

$\mathcal{MP}$

) guided by the LP relaxation as follows. We first solve the LP relaxation of the monolith with the Benders decomposition approach as described in section 4.1.3 using only the Benders optimality cuts (16) and the continuous y _u,p and η variables. Once the LP relaxation has been solved to optimality, we inspect all computed Benders cuts. Specifically, we collect all occurring tuples

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

in a set

${\mathcal{T}}^{\mathit{pu}}\subset \mathcal{C}\times \mathcal{Q}$

, that is, we collect those tuples with a positive dual variable

${\pi }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime }}^{*}>0$

. Next, we project

${\mathcal{T}}^{\mathit{pu}}$

onto the “position‐space,” that is, we define

$${\mathcal{T}}^p:=\{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in \mathcal{C}|\exists (u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in \mathcal{Q}:(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in {\mathcal{T}}^{\mathit{pu}}\}\subset \mathcal{C}.$$

Note that for most instances,

${\mathcal{T}}^p$

contains much fewer elements than

$\mathcal{C}$

.

This allows us to define the strengthened master problem (

${\mathcal{MP}}^{\star }$

), where besides some initial Benders optimality cuts (16), we have also integrated the additional cuts (22).

$$z_{M^{*}}^{*}:=\text{min}\eta$$

20

$$\text{s.t.}(15)-(17)$$

21

$$\eta \ge \sum _{(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in {\mathcal{T}}^p}{\tau }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}$$

22

$${\tau }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge {\pi }_{u,p,c}^{\mathrm{var}}{y}_{u,p}+{\pi }_c^{\mathrm{con}}\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })\in {\mathcal{T}}^{\mathit{pu}},c\in \mathcal{O}$$

23

$${\tau }_{p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime }}\ge 0\forall (p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in {\mathcal{T}}^p.$$

24

Remark 6

Note that the cuts in (22) are valid cuts and, at the same time, stronger than the related Benders cuts. To see that, observe that a Benders cut is constituted by the sum over all tuples of units and positions that cause crossings. Cuts in (22) separate these sums onto multiple cuts chosen by the tuple of overlapping positions, implying that these cuts are valid and stronger. This comes at the cost that we have to introduce an additional variable for each cut of type (22) whereas Benders cuts are not enlarging the space of decision variables.

While (

${\mathcal{MP}}^{\star }$

) might look similar to the monolithic model (

$\mathcal{M}$

), the key aspect in the formulation of

$({\mathcal{MP}}^{\star })$

is that considerably fewer variables and constraints of type (2) are needed. Thus,

$({\mathcal{MP}}^{\star })$

is a hybrid version that combines the strength of the monolithic formulation with the nimbleness of its Benders decomposition. This way, the Benders algorithm resembles the delayed row generation approach in the restricted set

${\mathcal{T}}^p$

.

Starting from a feasible solution

$\bar{y}$

of

$({\mathcal{MP}}^{\star })$

we suggest the following procedure:

Compute the tuples

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

induced by all crossings of

$\bar{y}$

and check whether we have

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in {\mathcal{T}}^p$

.

If

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\notin {\mathcal{T}}^p$

, add a Benders cut of type (16) over all

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

with

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\notin {\mathcal{T}}^p$

.

If

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in {\mathcal{T}}^p$

, add cuts of type (22) for all

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime },u,u^{\prime },\widehat{u},{\widehat{u}}^{\prime })$

with

$(p,p^{\prime },\widehat{p},{\widehat{p}}^{\prime })\in {\mathcal{T}}^p$

.

The Tailored Benders Based Branch‐and‐Cut Algorithm

The resulting algorithm is summarized in pseudo‐code in Algorithm 1. The algorithm is initialized in step 1. Because the master problem

$(\mathcal{MP})$

is unbounded for empty cut set

$\mathcal{O}$

, we set an artificial lower bound on η. This lower bound can be any value ≤ 0 because 0 is a trivial lower bound of the LP relaxation of

$(\mathcal{M})$

.

In the first “repeat” loop (steps 2–6), the LP relaxation of the monolith

$(\mathcal{M})$

is solved through Benders decomposition (McDaniel and Devine 1977). Every trial solution from the master problem (step 3) is evaluated to obtain its “true” objective function value

$z_S^{*}(\widehat{y})$

(step 4). If this objective function value matches the objective function value of the master problem, then the LP relaxation has been solved and the “repeat” loop is left (step 5). Otherwise, a Benders cut for the fractional‐valued trial

$\widehat{y}$

is added (step 6). Step 7 checks integrality of the computed incumbent

$\widehat{y}$

; integrality implies that

$\widehat{y}$

is an optimal solution of the monolith

$(\mathcal{M})$

. Step 8 prepares the sets required for the strengthened Benders master problem cuts (22).

In the second “repeat” loop (steps 9–20), the master problem is solved to optimality by adding Benders cuts until all integral trial solutions are evaluated with the correct objective function value. Therefore, the master problem is solved (step 10) either until optimality (then the algorithm terminates in step 11) or until an integral trial solution

$\widehat{y}$

has been computed. If the objective function value of this trial solution is evaluated correctly in the master problem, then the solution is checked for possible update (step 14) and the corresponding node in the branch‐and‐bound tree is fathomed (step 15). Otherwise, Benders optimality cuts (16) or (22) are added (step 18 and 19, respectively).

The Dynamic Fix‐and‐Relax Pump (DFRP)

To improve the best feasible point found by BDC, we introduce the so‐called Dynamic Fix‐and‐Relax Pump (DFRP). DFRP is inspired by other numerical methods in large‐scale optimization such as block coordinate descent (Hildreth 1957, Wright 2015) and stochastic gradient descent (Bottou 2010, Robbins and Monro 1951) where, due to the model size, not all the information available is used in each iteration. We also borrow ideas from the fix‐and‐relax heuristic, which successively improves a feasible solution by fixing parts of the variables and optimizing the resulting problem only in the free remaining variables (Toledo et al. 2015). Because ECW2CN does not exhibit a sequential decision‐making process which can be exploited algorithmically, the proposed method differs from rolling horizon approaches in that any decision variables can be freed or fixed throughout the algorithm (Thevenin et al. 2021).

Based on a feasible solution

$\bar{y}$

, we select randomly a batch of units

$B\subseteq \mathcal{U}$

. For these units, an optimal positioning is computed using BDC. The batch size of units to be optimized is controlled in a dynamic way: If we see the same optimal value for some iterations, the batch size is increased. Once we reach a batch size, such that the corresponding problem cannot be solved to optimality within the time limit, we reset the sample size to an initial value. DFRP is stated explicitly in Algorithm 2.

Stopping Criteria and Deterministic Variations of DFRP

We want to discuss some stopping criteria and alternative versions of DFRP:

A maximal number of iterations, a time limit or a maximal number of iterations without changes in the objective value are classical stopping criteria which, of course, may also be used for DFRP.

Note that DFRP computes a global optimal point of

$(\mathcal{M})$

, if we have z _best = 0 or if the problem in step 4 of Algorithm 2 is solved to optimality for

$|B|=|\mathcal{U}|$

.

One deterministic variation of the DFRP is to drop step 5 if Algorithm 2 and to increase the batch sizes until

$|B|=|\mathcal{U}|$

. This variation can also be seen as warm‐started BDC, preceded by a “trajectory” of feasible points generated by DFRP.

Solution of the Real‐World Problem at BASF

Problem Specifications

In this use case, we have 75 units with 375 connections that have to be located on 76 positions. The maximal tube length permits connecting two units on the same side with a distance of 20 positions and two units on opposing sides whose distance does not exceed 18 positions. The maximal curvature condition translates into the constraint that connected units may not be located on neighbored positions.

Implementation Details

All computations are run on an Intel Xeon Processor with 4 cores at 3.7 GHz and 128 GB RAM. The code was implemented using Python 3.7 and GuRoBi 9.1 using the following GuRoBi methods/attributes:

GRB.Callback.MIPSOL: This attribute indicates the occurrence of a new incumbent. If a new feasible point was computed, we add a new optimality cut with Model.cbLazy().

Model.cbLazy() which adds the cuts in a lazy manner.

The attribute GRB.Callback.MIPNODE can be used to get an entry points for actions that are triggered by the computation of new feasible points of the continuous relaxation, for example, for adding fractional cuts in the strengthening step of the master problem.

GRB.Callback.MIPNODE_NODCNT and master.terminate() are used to disrupt the solve routine after all fractional cuts are collected in the root node.

For Algorithm 2, we set b _init  = 6, b ₊ = 1, ɛ = 10⁻⁸ and s _max  = 5 for all computations. The time limit for each problem within DFRP is set to 2 hours and we set the tolerable relative MIP gap to 1%.

DFRP Statistics

To strengthen the master problem as described in section 4.1.4, we collect 47,265 fractional cuts while solving the root node, that is, we have

$|{\mathcal{T}}^p|=47,265$

. Based on the same initial solution, we run DFRP for each of the four objective functions for 48 hours. In each iteration of DFRP, an ECW2CN instance is solved for partially fixed units. The solution statistics are summarized for the different objective functions in Tables 1–4. The column “Sample Size” refers to the number of non‐fixed units whose position is to determine optimally. “#Entries” is the number of instances for each row, for example, we have 25 iterations of DFRP for the objective function “min” with sample size 6. Column “Av. # Opt. Cuts” describes the average number of optimality cuts introduced for each problem instance in that row, whereas its average relative optimality is given in column “Av. Gap.” Finally, “Av. CPU Time [s]” and “Av. # BB Nodes” describe the average CPU time in seconds and the number of searched branch‐and‐bound nodes, respectively.

Table 1

Numerical Results for Objective Function “minimum”

Sample size	# Entries	Av. # Opt. Cuts	Av. Gap	Av. CPU Time [s]	Av. # BB Nodes
6	25	428,843	0.0041	100.3	11.08
7	27	426,632	0.0034	61.6	13.37
8	29	515,654	0.0045	128.0	17.76
9	64	837,583	0.0061	138.6	44.59
10	41	1,974,287	0.0073	570.7	161.68
11	57	4,505,610	0.0079	873.3	439.19
12	30	6,621,153	0.0094	1820.4	815.70
13	8	11,584,127	0.0287	3955.1	2027.88

Table 2

Numerical Results for Objective Function “product”

Sample size	# Entries	Av. # Opt. Cuts	Av. Gap	Av. CPU Time [s]	Av. # BB Nodes
6	24	306,481	0.0033	148.4	4.92
7	27	521,725	0.0038	307.1	24.67
8	29	600,600	0.0045	543.2	25.76
9	43	683,431	0.0035	525.7	34.35
10	16	2,086,097	0.0149	1572.5	184.62
11	10	2,577,913	0.0078	1724.9	226.30
12	38	3,032,248	0.0079	1825.7	313.95
13	4	3,582,362	0.0290	2542.0	781.50

Table 3

Numerical Results for Objective Function “one”

Sample size	# Entries	Av. # Opt. Cuts	Av. Gap	Av. CPU Time [s]	Av. # BB Nodes
6	36	422,135	0.0034	361.1	15.00
7	24	835,595	0.0046	758.0	47.46
8	27	756,887	0.0047	877.9	37.63
9	43	1,290,207	0.0061	1396.5	77.67
10	22	1,837,311	0.0062	1603.9	131.55
11	9	3,165,046	0.0081	2572.3	252.33

Table 4

Numerical Results for Objective Function “sum”

Sample size	# Entries	Av. # Opt. Cuts	Av. Gap	Av. CPU Time [s]	Av. # BB Nodes
6	45	373,946	0.0035	93.6	11.04
7	38	378,917	0.0022	91.1	10.97
8	48	452,494	0.0039	135.1	14.04
9	44	830,338	0.0047	221.3	44.32
10	75	1,442,428	0.0062	245.2	101.83
11	39	3,294,341	0.0081	990.6	324.33
12	35	2,753,871	0.0071	732.3	270.54
13	19	9,086,700	0.0183	3413.9	1453.05

In Tables 1–4 we observe that the sample sizes vary between 6 and 13 for most objective functions and we see that in general the number of optimality cuts increases with the sample size. The optimality gaps reported by GuRoBi are close to zero and exceed the given threshold of 1% only if the model cannot be solved within the time limit of 2 hours. As expected, the average computing time as well as the number of searched branch‐and‐bound nodes increases for larger sample sizes.

The “one” objective function seems to generate rather hard problem instances since we always faced a problem instance of sample size 11 that could not be solved within the time limit of 2 hours before increasing the sample size (cf. Table 3).

Results

Each objective functions leads between 60 and 83 feasible points. These points are then evaluated by a simulation model that counts the total number of crossings in a time dependent framework, based on the historic connections. The result is depicted in Figure 6.

Figure 6

For Each Feasible Point, the Objective Function Value on the Left Axis and the Number of Crossings as Output of the Simulation Model are Shown. Each Graph Corresponds to a Different Objective Function and a Different Run of Algorithm Dynamic Fix‐and‐Relax Pump [Color figure can be viewed at wileyonlinelibrary.com]

Each of the four graphs in Figure 6 is generated as follows. For a specific objective function, all computed feasible points by DRFP are sorted according to decreasing objective function value. This is shown by the blue dots. For each feasible point, the resulting crossings, as evaluated by the simulation model, are shown by the red dots. Therefore, in a perfect (that is, dynamic) model, the red dots would lie on the blue line. We observe that the two y‐axis values live on different scales.

Next, we want to evaluate how good the obtained solutions (“blue dots”) mimic the simulation results (“red dots”); therefore, consider Table 5. This table lists both the Pearson correlation and the Spearman's rank correlation of the objective function value to the simulation result. Such as the correlation, the rank correlation lies between −1 and +1 and a value close to +1 indicates strong dependence. While the correlation measures linear dependence, the rank correlation captures monotonic relationships, that is, only the ranks of the data points are compared and not their values. Therefore, the rank correlation is a more relevant measure for our data than the correlation because we are just interested which solutions perform best.

Table 5

Comparison of Relationship between Objective Function Value and Simulation Result

Objective function	Correlation	Rank correlation
Minimum	0.91	0.67
Product	0.90	0.96
One	0.91	0.90
Sum	0.84	0.73

For the results in Figure 6, the “product” objective function value performs very strong for both the correlation and the rank correlation. The “sum” and “minimum” objectives perform well with respect to the correlation of objective values and simulation results but show weaknesses with respect to rank correlation. That can also be seen in Figure 6 since for both objectives there exist points with very good objective values but poor simulation performance. It is striking to see that also the “one” objective performs very well for both performance measures. However, while the results for “sum,” “minimum,” and “product” have been very similar in previous versions of the numerical results, “one” did behave much worse in these runs. So, the good performance of “one” might be a ”lucky shot” caused by the stochastic nature of DFRP.

Finally, for each objective function, the best feasible points with respect to their objective value and the simulation result are evaluated with the second, much more detailed simulation model. Because for the “one” objective these points coincide, we computed the number of FIFO crossings for seven feasible points and ended up with the results in Table 6.

Table 6

Computation of Critical FIFO Crossings for the Best Seven Solutions Found by Dynamic Fix‐and‐Relax Pump (DFRP) Evaluated on a Time Period of Three Years

Instance	Number of FIFO crossings
Minimum ‐ best_obj	17
Minimum ‐ best_sim	8
Product ‐ best_obj	10
Product ‐ best_sim	0
One ‐ best_obj and best_sim	26
Sum ‐ best_obj	26
Sum ‐ best_sim	18

The best solution found without any FIFO crossings is illustrated in Figure 7.

Figure 7

Best Solution with Regard to the Number of Critical FIFO Crossings Computed by Dynamic Fix‐and‐Relax Pump

The final solution was presented to the engineering department. Although we computed a fully feasible layout in the sense that no FIFO crossings occurred on historical connection data, we were aware of the fact that that might not hold true for futures connections. However, the best layout generated by our algorithm reduces the expected number of such events to a few per year, which was deemed to be acceptable by the engineers.

As a side effect, our models provided valuable insights for the engineers. In particular, the question of how many pipes of the different lengths will be needed was an important parameters for the dimensioning of the tube handling system. Also, more general questions, such as the effect of the width of the switching hub (the distance between the two lines of connectors) on other parameters, could be answered easily.

Our results removed the last concerns toward the applicability of the robotic invention, and a decision was made to build it. The robotic system will automatically prevent connections between wrong units, decrease the risk of spillage and increase the safety and satisfaction of the workforce.

Conclusions

In an effort to increase human safety, BASF is developing a dual‐arm robotic system to reconnect pipes in a switching hub. The robot can undertake almost all reconnections, except such reconnections which require the removal of multiple pipes at ones, which is left to the humans. This raised the question on how to design the layout of a new switching hub in order to minimize the number of reconnections involving humans. We have modeled this real‐world problem as a special minimal crossing number problem in conjunction with two simulation models. The resulting edge‐constrained weighted two‐layer crossing number problem (ECW2CN) leads to instances of gigantic size (approximately 70 billion constraints), when tackled via MILP approaches. Therefore, we have developed a tailored Benders decomposition which is employed in a dynamic fix‐and‐relax pump. For different objective functions, we have identified a total of 454 feasible layouts. These layouts are then evaluated by two different simulation models. The obtained solution not only increases human safety but also provides valuable feedback to engineers regarding the design of the switching hub.

We see three possible directions for future academic research. The first direction is to increase the computational efficiency of the tailored solution algorithms or to propose better performing models and methods. For Benders decomposition, there might be possibilities to further strengthen the master problem through tailored cuts. Another possibility might arise when combining column generation (on the τ‐variables) with Benders decomposition. A second direction is to enhance the model by dynamic and stochastic aspects, capturing the entire complexity of the real‐world problem. This comes at the cost that: (1) no real‐world data are available, and (2) the model becomes even more challenging to solve. A third direction is to obtain better theoretic lower bounds through specific crossing number lemmas, taking the special structure of ECW2CN into account.