Computing tight bounds via piecewise linear functions through the example of circle cutting problems

Abstract

This paper discusses approximations of continuous and mixed-integer non-linear optimization problems via piecewise linear functions. Various variants of circle cutting problems are considered, where the non-overlap of circles impose a non-convex feasible region. While the paper is written in an “easy-to-understand” and “hands-on” style which should be accessible to graduate students, also new ideas are presented. Specifically, piecewise linear functions are employed to yield mixed-integer linear programming problems which provide lower and upper bounds on the original problem, the circle cutting problem. The piecewise linear functions are modeled by five different formulations, containing the incremental and logarithmic formulations. Another variant of the cutting problem involves the assignment of circles to pre-defined rectangles. We introduce a new global optimization algorithm, based on piecewise linear function approximations, which converges in finitely many iterations to a globally optimal solution. The discussed formulations are implemented in GAMS. All GAMS-files are available for download in the Electronic supplementary material. Extensive computational results are presented with various illustrations.

Appendices

Appendix 1: A general properties of approximator systems

In this appendix, we provide some definitions for terminology and we generalize some of the concepts used in Sects. 3–7.

We are interested in the following non-linear optimization problem:

$$\begin{aligned} \mathrm {(NOP)}: \quad \min&\,f(x) \end{aligned}$$

(54)

$$\begin{aligned} \text {s.t.}&\, g(x)=0 \end{aligned}$$

(55)

$$\begin{aligned}&\, h(x)\le 0 \end{aligned}$$

(56)

$$\begin{aligned}&\,x \in \mathbb {D} := [X_-, X_+]^n \subset \mathbb {R}^n \end{aligned}$$

(57)

with continuous functions $f\,{:}\,\mathbb {D}\rightarrow \mathbb {R}, g\,{:}\,\mathbb {D} \rightarrow \mathbb {R}^{m_{1}}$ and $h\,{:}\,\mathbb {D}\rightarrow \mathbb {R}^{m_{2}}$. We also allow for special cases when the number of equality constraints, $m_1$, is zero, or the number of inequality constraints, $m_2$, is zero. We are particularly interested in those cases where problem (54)–(57) is a non-convex optimization problem, for example, f(x) is a non-convex function in x, and/or constraint set (55)–(56) defines a non-convex feasible region.

Notice that integrality on the decision variable x can be enforced via special non-convex constructs (Horst et al. 2000). Thus, (NOP) models the general class of mixed-integer non-linear programming problems where lower and upper bounds on all decision variables are present. With that, the continuous NLP cutting problems $(\mathcal {C}^\mathrm{length})$ and $(\mathcal {C}^\mathrm{area})$ as well as the MINLP assignment–cutting problem $(\mathcal {A})$ are special cases of (NOP).

We are not interested in solving problem (54)–(57) to global optimality. In fact, we assume that there are global solvers available, which can solve this optimization problem (or variations of it) to global optimality. Rather, we are interested in finding linear $\varepsilon$-approximations of the optimization problem (54)–(57) in the following sense:

Definition 1

($\varepsilon$-approximated problem) Let $x^*$ be a globally optimal solution to (54)–(57) and $| \cdot |$ denote the absolute value function. We call an optimization problem an $\varepsilon$ -approximated problem, iff for any optimal solution $y^*$ of the $\varepsilon$-approximated problem, $y^*$ satisfies the following properties:

$$\begin{aligned}&P_1: \quad \left| g_i(y^*) \right| \le \varepsilon \quad i \in \{ 1, \ldots , m_1 \} \\&P_2: \quad h_j(y^*) \le \varepsilon \quad j \in \{ 1, \ldots , m_2 \} \\&P_3: \quad \left| f(x^*) - f(y^*) \right| \le \varepsilon . \end{aligned}$$

Linear $\varepsilon$-approximations satisfying properties $P_{1}$–$P_{3}$ often require non-convex and non-concave piecewise linear functions, which can be expressed via linear functions and breakpoints. Thus, one obtains an MILP, approximating the non-linear optimization problem (54)–(57). The size of the resulting MILP depends crucially on the number of breakpoints and one can expect that the MILP has significantly more variables than the original (MI)NLP (54)–(57).

The definition of the $\varepsilon$-approximated problem assumes that both the original non-linear optimization problem as well as the approximated problem are feasible. However, to safely detect infeasibility of the non-linear optimization problem (54)–(57) is an interesting problem by itself. Particularly, for an $\varepsilon$-approximated problem, one might demand the additional property

$$\begin{aligned} P_4:&\text {if the }\, \varepsilon \text {-approximated problem is infeasible, then the global optimization problem}\\&(54){-}(57)\,\text {is also infeasible}. \end{aligned}$$

If an optimization problem is infeasible, then we use the convention for minimization problems that the “optimal” objective function value is infinity, e.g., $f(x^{*})=\infty$. We note that feasible solutions to the $\varepsilon$-approximated problem might not be feasible to the (NOP). This relaxation enables us to obtain valid lower bounds on the optimal objective function value of (NOP).

Let us recall the definition of a $\delta$-approximator (in n-dimensions)

Definition 2

($\delta$-approximator, Rebennack and Kallrath 2015a) Let $f{:}\, \mathbb {D} \rightarrow \mathbb {R}$ be a function on the compactum $\mathbb {D}$ and let scalar $\delta > 0$. A piecewise linear, continuous function $\ell \,{:}\,\mathbb {D} \rightarrow \mathbb {R}$ is called a $\delta$ -approximator for f, iff the following property holds

$$\begin{aligned} \max _{x \in \mathbb {D}} \left| \ell (x)-f(x)\right| \le \delta . \end{aligned}$$

(58)

as well as a $\delta$-underestimator and a $\delta$-overestimator

Definition 3

($\delta$-underestimator/$\delta$-overestimator, Rebennack and Kallrath 2015a) Let scalar $\delta > 0$. We call function $\ell \,{:}\, \mathbb {D} = D_1 \rightarrow \mathbb {R}$ a $\delta$ -underestimator of function $f\,{:}\,\mathbb {D} \rightarrow \mathbb {R}$, iff condition (58) is satisfied along with

$$\begin{aligned} \ell (x) \le f(x) \quad \forall x \in \mathbb {D}. \end{aligned}$$

(59)

We call function $\ell$ a $\delta$ -overestimator of function f, iff $-\ell$ is a $\delta$-underestimator of $-f$.

The approximation quality of the $\delta$-approximators, $\delta$-underestimators and $\delta$-overestimators are defined through the maximal absolute difference between the function and the approximators over the given compactum $\mathbb {D}$. Later in this appendix, we see how that is useful when approximating (NOP). However, other quality measures for approximators exist, for instance, the area between the approximator and the function (Kallrath and Rebennack 2014a).

We call a $\delta$-approximator, $\delta$-underestimator and $\delta$-overestimator tight, if there is no approximator, under- or overestimator with the same number of breakpoints having a maximal absolute function deviation smaller than $\delta$.

Next, we apply $\delta$-approximator functions towards NLP (54)–(57). We obtain

Theorem 1

Given is problem (54)–(57). For a $\delta$-approximator $\ell$ over $\mathbb {D}$ for objective function f with $\delta := \frac{\varepsilon }{2}$, any optimal solution of

$$\begin{aligned} \min&\, \ell (y) \end{aligned}$$

(60)

$$\begin{aligned} \mathrm{s.t.}&\, g(y) = 0 \end{aligned}$$

(61)

$$\begin{aligned}&\, h(y) \le 0 \end{aligned}$$

(62)

$$\begin{aligned}&\, y \in \mathbb {D} \end{aligned}$$

(63)

satisfies properties $P_1$–$P_3$. Furthermore, (60)–(63) satisfies $P_4$. Replacing $\ell$ with an $\varepsilon$-underestimator for f yields the same properties.

Proof

Assume that (60)–(63) is infeasible. As the feasible region of (61)–(63) is identical to the feasible region of (55)–(57), one obtains $P_4$.

Now, assume that (60)–(63) is feasible. For $P_3$, let $y^*$ be an optimal solution to (60)–(63) and $x^*$ be an optimal solution to (54)–(57). Assume that

$$\begin{aligned} f(x^*) \notin \overline{B}_{\varepsilon }\big (f(y^*)\big ), \end{aligned}$$

where the closed ball is defined as

$$\begin{aligned} \overline{B}_{\mu }(z):= \big \{ x \, : \left| x - z \right| \le \mu \big \}. \end{aligned}$$

By (58), we obtain

$$\begin{aligned}&\ell (y^*) \in \overline{B}_\delta \big ( f(y^*) \big )\quad \text {and} \quad \ell (x^*) \in \overline{B}_\delta \big ( f(x^*) \big ). \end{aligned}$$

By assumption, $\overline{B}_\delta \big ( f(y^*) \big ) \cap \overline{B}_\delta \big ( f(x^*) \big ) = \emptyset$ and because $f(x^*) \le f(y^*)$ due to the optimality of $x^*$, we have that $\ell (x^*) < \ell (y^*)$ which contradicts optimality of $y^*$.

We have that $f(y^*) - \ell (y^*) \le \varepsilon$ which implies $f(y^*) - f(x^*) \le \epsilon$. $\square$

One might expect that $x^* = y^*$ in most cases. However, note that $x^*$ does not even have to be “close” to $y^*$ in general.

Theorem 1 shows that if the “difficult” non-linear relations only occur in the objective function, then linear $\varepsilon$-approximations work fine to approximate the original non-linear optimization problem. However, if non-linear terms occur in difficult to handle constraints, such as equality constraints, we rather resort to under- and overestimators.

By relaxing the feasible region (outer approximation), lower bounds on the optimal objective function of the original problem can be calculated. This is formalized in the next

Theorem 2

(Outer approximation) Given is problem (54)–(57). With functions $\ell _-, g_{i-}, h_{j-}$ being $\varepsilon$-underestimators for $f, g_i$ and $h_i$, respectively, and $g_{i+}$ being a $\varepsilon$-overestimator for $g_{i}$, any optimal solution $y^*$ of

$$\begin{aligned} \min&\, \ell _-(y) \end{aligned}$$

(64)

$$\begin{aligned} \text {s.t.}&\, g_{i-}(y) \le 0, \quad \quad i = 1, \ldots , m_1 \end{aligned}$$

(65)

$$\begin{aligned}&\, g_{i+}(y) \ge 0, \quad \quad i = 1, \ldots , m_1 \end{aligned}$$

(66)

$$\begin{aligned}&\, h_{j-}(y) \le 0, \quad \quad j = 1, \ldots , m_2 \end{aligned}$$

(67)

$$\begin{aligned}&\, y \in \mathbb {D} \end{aligned}$$

(68)

satisfies properties $P_1, P_2$ and $\ell _-(y^*) \le f(x^*)$. Furthermore, (64)–(68) satisfies $P_4$.

Proof

$P_1$ and $P_2$ hold by construction. The condition $\ell _-(y^*) \le f(x^*)$ and $P_4$ hold true because any feasible solution x to (54)–(57) is also feasible for (64)–(68). $\square$

By Theorem 2, any optimal solution of (64)–(68) defines a (valid) lower bound to the original problem (54)–(57). An upper bound can be obtained, e.g., by running a local solver (using the obtained solution $y^*$).

Note that Theorem 2 remains true when any of the exact functions f, g or h are used. With that, our MILP cutting problem $(\mathcal {P}^\mathrm{length})$ is a special case for (64)–(68) when the quadratic functions in the non-overlap constraints are overestimated. The non-overlap constraints are parts of the $-h(x)$ function in (54)–(57). Similarly for $(\mathcal {P}^\mathrm{area})$ when in addition underestimators are used in the objective.

One may wonder, if it is possible to prove a stronger version of Theorem 2 along the lines of Theorem 1. To do so, consider the following one-dimensional example for any positive scalars M and $\vartheta \le 1$:

$$\begin{aligned} \ \min&M \cdot x_1 \\ \text {s.t.}&-x_1^2 + \vartheta \le 0 \\&x_1 \in \{0,1\}, \end{aligned}$$

with the (unique) optimal solution $x_1^* = 1$ and optimal objective function value $f(x_1^*) = M$. Let $\ell _-$ be an $\varepsilon$ -underestimator of $h(x):= - x_1^2 + \vartheta$ with $\varepsilon \ge \vartheta$. Thus, the optimization problem of the form

$$\begin{aligned} \min&\, M \cdot y_1 \end{aligned}$$

(69)

$$\begin{aligned} \text {s.t.}&\, \ell _-(y_1) \le 0 \end{aligned}$$

(70)

$$\begin{aligned}&\, y_1 \in \{0,1\}, \end{aligned}$$

(71)

may allow $y_1^* = 0$ as an optimal solution of (69)–(71) with objective function value 0. (The interested reader may construct an example over a continuous domain, by replacing the domain by [0, 1] and adding the term $M \cdot x_1 \cdot (1-x_1)$ to the objective, and working with an $\varepsilon$-underestimator for the objective). For this example, the phenomenon of arbitrary large differences in the objective function is avoided by using an $\varepsilon$-underestimator for h(x) with $\varepsilon < \vartheta$.

When shrinking the feasible region (inner approximation), one may obtain an upper bound on the objective function of the original problem if the approximated problem is feasible.

Theorem 3

(Inner approximation) Given is problem (54)–(57). With functions $\ell _{+}$ and $h_{j+}$ being $\varepsilon$-overestimators for f and $h_j$, any optimal solution $y^*$ of

$$\begin{aligned} \min&\, \ell _+(y) \end{aligned}$$

(72)

$$\begin{aligned} \text {s.t.}&\, g(y) = 0 \end{aligned}$$

(73)

$$\begin{aligned}&\, h_{j+}(y) \le 0, \quad \quad j = 1, \ldots , m_2 \end{aligned}$$

(74)

$$\begin{aligned}&\, y \in \mathbb {D} \end{aligned}$$

(75)

satisfies properties $P_1, P_2$ and $\ell _+(y^*) \ge f(x^*)$.

Proof

$P_1$ and $P_2$ hold by construction. Because any feasible point to (72)–(75) is also feasible for (54)–(57), and function $\ell _+$ is an overestimator for objective function $f, \ell _+(y^*) \ge f(x^*)$ holds. $\square$

Note that $P_4$ does not hold for the inner approximation (72)–(75).

Similar to the discussion on outer approximations, $(\mathcal {P}^\mathrm{length})$ and $(\mathcal {P}^\mathrm{area})$ define an inner approximation if the quadratic functions in the non-overlap constraints are underestimated.

Appendix 2: Discretization of rectangles

Instead of approximating the non-overlap constraints (1) of the circles via piecewise linear constructs, we restrict the placement of the circle centers on predefined grid points. Therefore, let $[0,S_1^+]$ be discretized at points $\chi _l$ with

$$\begin{aligned} 0< \chi _l< \chi _{l+1} < S_1^+ \quad \forall l \in \left\{ 1, \ldots , K^\mathrm{R}_1-1 \right\} , \end{aligned}$$

where set $\mathbb {B}^\mathrm{rect}_1 = \{1, \ldots , K^\mathrm{R}_1\}$ collects all indices of the discretization points. Similarly, the width is discretized as

$$\begin{aligned} 0< \chi _w< \chi _{w+1} < S_2^+ \quad \forall w \in \left\{ 1, \ldots , K^\mathrm{R}_2-1 \right\} , \end{aligned}$$

with set $\mathbb {B}^\mathrm{rect}_2$. We label all grid points inside the rectangle accordingly by lw.

Binary variable $x^\mathrm{b}_{ilw}$ then has value 1, if the circle center $i \in \mathbb {I}$ is placed on the grid point lw and 0 otherwise. The non-overlap constraints (1) then simplify to the linear constraints

$$\begin{aligned} x^\mathrm{b}_{ilw} + x^\mathrm{b}_{jl'w'} \le 1 \quad \forall i \in \mathbb {I}, \; j \in \mathbb {I}, \; l \in \mathbb {B}^\mathrm{rect}_1, \; w \in \mathbb {B}^\mathrm{rect}_2 \quad \text {with}\;\; i<j \quad \text {and}\quad l'w'\in O_{ijlw}.\nonumber \\ \end{aligned}$$

(76)

If circle center i is placed on grid point lw, then set $O_{ijlw}$ contains all such grid points $l'w'$ at which circle center j cannot be placed on. For i, j and lw, set $O_{ijlw}$ contains grid point $l'w'$, if

$$\begin{aligned} (\chi _{l} - \chi _{l'})^2 + (\chi _{w} - \chi _{w'})^2 < (R_i + R_j)^2. \end{aligned}$$

(77)

Each circle needs to be placed exactly once

$$\begin{aligned} \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} x^\mathrm{b}_{ilw} = 1 \quad \forall i \in \mathbb {I}. \end{aligned}$$

(78)

To ensure that the circle fits inside the rectangle, we formulate conditions (2) and (3) in terms of the binary variable $x^\mathrm{b}_{ilw}$. Therefore, we approximate $x_{ id}$ through $x^\mathrm{b}_{ilw}$. For the length, or dimension $d=1$, we identify

$$\begin{aligned} x_{i1} \approx \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _l \cdot x^\mathrm{b}_{ilw}\quad \forall i \in \mathbb {I} \end{aligned}$$

and for the width, or $d=2$,

$$\begin{aligned} x_{i2} \approx \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _w \cdot x^\mathrm{b}_{ilw}\quad \forall i \in \mathbb {I}. \end{aligned}$$

Circle fitting constraints (2) are then enforced by fixing all decision variables $x^\mathrm{b}_{ilw}$ to value 0 for which

$$\begin{aligned} \chi _l \le R_i \quad \text {or} \quad \chi _w \le R_i. \end{aligned}$$

Instead of adding constraints of the type “$x^\mathrm{b}_{ilw}=0$”, we can eliminate these decision variables from the model formulation in the pre-processing step.

Constraints (3) translate then to

$$\begin{aligned} \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _l \cdot x^\mathrm{b}_{ilw} + R_i \le x_1^\mathrm{P} \quad \text {and} \quad \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _w \cdot x^\mathrm{b}_{ilw} + R_i \le x_2^\mathrm{P} \quad \forall i \in \mathbb {I}.\nonumber \\ \end{aligned}$$

(79)

We obtain the following MILP problem, which minimizes the design rectangle’s length over the discretized rectangle, as

In our second cutting problem $(\mathcal {C}^\mathrm{area})$, we minimize the area of the design rectangle. This leads to the bilinear term

$$\begin{aligned} x_{1}^{\mathrm {P}} \cdot x_{2}^{\mathrm {P}}. \end{aligned}$$

(80)

Because both $x_{1}^{\mathrm {P}}$ and $x_{2}^{\mathrm {P}}$ are continuous decision variables, we cannot find an exact linear representation for (80). Thus, we use an approximation. The idea here is again to use a discretization of the length and the width, just like we did for the circle center coordinates in $(\mathcal {D}^\mathrm{length})$. Let $\chi _m$ discretize the rectangle’s length in (80) ($\mathbb {B}^\mathrm{obj}_1 = \{ 1 \ldots , K^\mathrm{O}_1 \}$)

$$\begin{aligned} S_1^-< \chi _m< \chi _{m+1} < S_1^+ \quad \forall m \in \left\{ 1, \ldots , K^\mathrm{O}_1-1 \right\} \end{aligned}$$

and let $\chi _o$ discretize its width ($\mathbb {B}^\mathrm{obj}_2 = \{ 1 \ldots , K^\mathrm{O}_2 \}$)

$$\begin{aligned} S_2^-< \chi _n< \chi _{n+1} < S_2^+ \quad \forall n \in \left\{ 1, \ldots , K^\mathrm{O}_2-1\right\} . \end{aligned}$$

We obtain the approximations

$$\begin{aligned} x_{1}^{\mathrm {P}} \approx \sum _{m \in \mathbb {B}^\mathrm{obj}_1} \chi _m \cdot x_{1m}^\mathrm{P,b} \quad \text {and} \quad x_{2}^{\mathrm {P}} \approx \sum _{n \in \mathbb {B}^\mathrm{obj}_2} \chi _n \cdot x_{2n}^\mathrm{P,b} \end{aligned}$$

with the additional binary variables $x_{1m}^\mathrm{P,b}$ and $x_{2n}^\mathrm{P,b}$ respecting the constraints

$$\begin{aligned} \sum _{m \in \mathbb {B}^\mathrm{obj}_1} x_{1m}^\mathrm{P,b} = 1 \quad \text {and} \quad \sum _{n \in \mathbb {B}^\mathrm{obj}_2} x_{1n}^\mathrm{P,b} = 1. \end{aligned}$$

(81)

Alternatively, we can use “$\le$” constraints in (81)—we expect slightly better computational performance when using “$\le$” instead of “$=$”.

The design rectangle’s area (80) is then approximated by

$$\begin{aligned} \sum _{m \in \mathbb {B}^\mathrm{obj}_1} \sum _{n \in \mathbb {B}^\mathrm{obj}_2} \chi _m \cdot \chi _n \cdot x_{1m}^\mathrm{P,b} \cdot x_{2n}^\mathrm{P,b}. \end{aligned}$$

(82)

In contrast to (80), the bilinear terms in (82) are the product of two binary variables. This can be linearized exactly. For that we define for $m \in \mathbb {B}^\mathrm{obj}_1$ and $n \in \mathbb {B}^\mathrm{obj}_2$

$$\begin{aligned} x^\mathrm{a}_{mn}\le & {} x_{1m}^\mathrm{P,b} \end{aligned}$$

(83)

$$\begin{aligned} x^\mathrm{a}_{mn}\le & {} x_{2n}^\mathrm{P,b} \end{aligned}$$

(84)

$$\begin{aligned} x^\mathrm{a}_{mn}\ge & {} x_{1m}^\mathrm{P,b} + x_{2n}^\mathrm{P,b} - 1 \end{aligned}$$

(85)

$$\begin{aligned} x^\mathrm{a}_{mn}&\,&\text {binary} \end{aligned}$$

(86)

and observe that $x^\mathrm{a}_{mn} = x_{1m}^\mathrm{P,b} \cdot x_{2n}^\mathrm{P,b}$.

It remains to substitute the approximated length and with into the fitting constraint (79) to yield

$$\begin{aligned} \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _l \cdot x^\mathrm{b}_{ilw} + R_i\le & {} \sum _{m \in \mathbb {B}^\mathrm{obj}_1} \chi _m \cdot x_{1m}^\mathrm{P,b} \quad \text {and}\nonumber \\ \sum _{l \in \mathbb {B}^\mathrm{rect}_1} \sum _{w \in \mathbb {B}^\mathrm{rect}_2} \chi _w \cdot x^\mathrm{b}_{ilw}+\,R_i\le & {} \sum _{n \in \mathbb {B}^\mathrm{obj}_2} \chi _n \cdot x_{2n}^\mathrm{P,b}. \end{aligned}$$

(87)

This allows us to formulate the following pure integer programming (IP) problem, which approximates $(\mathcal {C}^\mathrm{area})$,

$$\begin{aligned} (\mathcal {D}^\mathrm{area}) \quad \tilde{A}^* = \min&\tilde{A} = \sum _{m \in \mathbb {B}^\mathrm{obj}_1} \sum _{n \in \mathbb {B}^\mathrm{obj}_2} \chi _m \cdot \chi _n \cdot x^\mathrm{a}_{mn} \\ \text {s.t.}&(76), (78), (83){-}(87) \\&x^\mathrm{b}_{ilw} \quad \text { binary} \quad \forall i \in \mathbb {I}, \; l \in \mathbb {B}^\mathrm{rect}_1, \; w \in \mathbb {B}^\mathrm{rect}_2 \\&x^\mathrm{a}_{mn}, \; x_{1m}^\mathrm{P,b}, \; x_{2n}^\mathrm{P,b} \quad \text { binary} \quad \forall m \in \mathbb {B}^\mathrm{obj}_1, \; n \in \mathbb {B}^\mathrm{obj}_2. \end{aligned}$$

By construction of the non-overlap constraints (76) via the set $O_{ijlw}$, we are able to recover a feasible point for $(\mathcal {C}^\mathrm{length})$ and $(\mathcal {C}^\mathrm{area})$ from an optimal solution of $(\mathcal {D}^\mathrm{length})$ and $(\mathcal {D}^\mathrm{area})$, respectively, assuming that there is at least one. Thus, by solving $(\mathcal {D}^\mathrm{length})$ and $(\mathcal {D}^\mathrm{area})$, we obtain an upper bound on the minimal rectangle’s length and minimal area, i.e.,

$$\begin{aligned} \tilde{L}^* \ge L^* \quad \text {and} \quad \tilde{A}^* \ge A^* \end{aligned}$$

for any discretization (assuming feasibility).

Because we are discretizing the bilinear objective function term via the discretization of the rectangle $[S_-^1,S_+^1] \times [S_-^2, S_+^2]$ and the binary decision variables $x_{1m}^\mathrm{P,b}$ and $x_{2n}^\mathrm{P,b}$, the area of the rectangle associated with the computed (feasible) placement of the circles is actually given by

$$\begin{aligned} \breve{A}= & {} \max _{i \in \mathbb {I}} \max _{l \in \mathbb {B}^\mathrm{rect}_1} \max _{w \in \mathbb {B}^\mathrm{rect}_2} \left\{ \chi _l \cdot x^\mathrm{b}_{ilw} + R_i \right\} \cdot \max _{i \in \mathbb {I}} \max _{l \in \mathbb {B}^\mathrm{rect}_1} \max _{w \in \mathbb {B}^\mathrm{rect}_2} \left\{ \chi _w \cdot x^\mathrm{b}_{ilw} + R_i \right\} \\= & {} l^* \cdot w^* \\\le & {} \tilde{A}^*, \end{aligned}$$

where $l^*$ and $w^*$ is the optimal length and width of the design rectangle computed by $(\mathcal {D}^\mathrm{area})$.

Next, we briefly discuss symmetry. In the NLP formulation, we used constraints (5) to break the symmetry caused by axial reflection (both horizontal and vertical) of any placement of the circles. However, this symmetry is not likely to be present in the derived formulation $(\mathcal {D}^\mathrm{area})$ because the grid is likely to not be symmetric with respect to the constructed design rectangle. The same holds true for $(\mathcal {D}^\mathrm{length})$ for vertical reflection. Horizontal reflection, however, can be broken by eliminating decision variables $x^\mathrm{b}_{1lw}$ with $\chi _w > S_2^+/2$ in the preprocessing step from the problem formulation. Thus, the number of variables is reduced.

The quality of the upper bounds computed depends on the discretization chosen. Actually, if we refine the grid, i.e., we add additional grid points to an existing grid (and do not change the grid points on the existing grid), we obtain a sequence of monotone decreasing length and area. Intuitively, this sequence converges to $L^*$ and $A^*$, respectively. However, we cannot readily provide a formula which bounds the maximal deviation, dependent on the discretization.

Lower bounds on $L^*$ and $A^*$ can also be obtained via the solution of $(\mathcal {D}^\mathrm{length})$ and $(\mathcal {D}^\mathrm{area})$ by allowing “enough” overlap of the circles. For instance, we can define the members of the set $O_{ijlw}$, for i, j and lw, if

$$\begin{aligned} (\chi _{l} - \chi _{l'})^2 + (\chi _{w} - \chi _{w'})^2 < ( \max \{R_i + R_j - 2\cdot \varrho , 0\} )^2, \end{aligned}$$

(88)

where $\varrho$ is the maximal diagonal of all rectangles imposed by the chosen discretization. In order to obtain valid bounds, the circles must respect the boundary conditions, i.e., they must fit inside the design triangle.

Table 11 Results for discretization of rectangle

Fig. 15

Graphical representation of the solutions to the discretization computed by CPLEX. a TC04 21-11. b TC04 21-11. c TC05 31-16. d TC05 31-16. e TC05 31-16. f TC05 31-16. g TC06 41-21. h TC06 41-21

Further, we obtain the following bounds on the length and area

$$\begin{aligned} \tilde{L} \le L^* \le \tilde{L} + (I-1) \cdot 2 \cdot \varrho \quad \text {and} \quad \breve{A} \le A^* \le (l^* + (I-1) \cdot 2 \cdot \varrho )\cdot (w^* + (I-1) \cdot 2 \cdot \varrho ), \end{aligned}$$

To illustrate the two formulations, computational results⁸ are reported in Table 11 for the instances of Table 1. Consider now Table 11. Columns two to four report on the discretization choice and the resulting maximal diagonal value $\varrho$. We have placed the grid equidistantly and the distance between two discretization points is equal in both dimensions. For formulation $(\mathcal {D}^\mathrm{area})$, we have chosen $K^\mathrm{O}_1 = K^\mathrm{R}_1$ and $K^\mathrm{O}_2 = K^\mathrm{R}_2$. Columns five to eight belong to formulation $(\mathcal {D}^\mathrm{length})$. $\underline{L}$ is the computed lower bound on $L^*$ while $\bar{L}$ is an upper bound. Similarly, columns nine to twelve provide the lower and upper bounds on the area $A^*$ as computed by $(\mathcal {D}^\mathrm{area})$. Computational times are given in seconds.

We observe that all bounds reported in Table 11 are correct but quite weak. The computational times are considerable, especially compared to global solvers, see Table 1. The reason for the weak bounds is the quite large discretization error, as measured in $\varrho$, even for a grid with as many as 52 and 27 segments. This is particularly visible when tracking the lower bound on the rectangle’s length, $\underline{L}$. We also notice that it is computationally more expensive to compute upper bounds (both for the length and the area) than lower bounds, though the formulations are identical. The reason is that problems computing the lower bounds contain much fewer constraints due to the less restrictive overlap conditions.

A few example plots in Fig. 15 illustrate the computed circle cuttings (upper bounds) as well as the resulting arrangements when allowing overlap of circles (lower bounds).