**Alexander Pankratov 1,2, Tatiana Romanova 1,2,\* and Igor Litvinchev <sup>3</sup>**


Received: 23 May 2020; Accepted: 3 July 2020; Published: 10 July 2020

**Abstract:** Packing irregular 3D objects in a cuboid of minimum volume is considered. Each object is composed of a number of convex shapes, such as oblique and right circular cylinders, cones and truncated cones. New analytical tools are introduced to state placement constraints for oblique shapes. Using the phi-function technique, optimized packing is reduced to a nonlinear programming problem. Novel solution approach is provided and illustrated by numerical examples.

**Keywords:** packing; irregular 3D objects; quasi-phi-function s; nonlinear optimization

#### **1. Introduction**

Packing problems aim to allocate a set of objects in a container subject to placement constraints. The latter typically stipulates non-overlapping between the objects and the boundary of the container. Additional placement constraints may include weight distribution and stacking, cargo stability, balance constraints, loading and unloading preferences, etc. [1]. In optimized packing certain criteria have to be optimized, e.g., maximizing the number of the packed objects, minimizing the waste or optimizing characteristics of the container, its volume or shape [2]. Packing problems are proved to be NP-hard [3].

Packing issues traditionally are important in logistics, e.g., in maritime transportation and container loading, cutting industrial materials in furniture and glass manufacturing [4]. Packing problems also arise in modelling liquid and glass structures, in analyses of powder and granular materials in mineral industry, in molecular and nanotechnologies [5–9].

Different classifications for packing problems are proposed (see, e.g., [2,4,10] and the references therein). Focusing on the shapes involved, packing problems can be divided into two large groups: regular and irregular. While regular 3D packing deals with relatively simple shapes (spheres, ellipsoids, convex polyhedrons), irregular packing focuses basically on nonconvex figures (see, e.g., [2,11–19]).

Various modelling and solution approaches, exact and approximated, are known for the regular packing (see, e.g., [20–22] and the references therein). For irregular 3D packing, heuristics are widely used [2]. A large group of heuristics is based on representing complex irregular shapes by corresponding collections of simpler (regular) figures thus reducing the problem approximately to a regular case [23–27]. Techniques in the other group combine a local search with simple decision rules, such as the deepest bottom-left approach or random allocation [28–33]. An alternative methodology is using genetic algorithms, directly or in combination with the first two approaches [34–37].

In this paper, packing irregular 3D objects in a cuboid of minimum volume is considered. Each complex object is composed of a number of convex shapes, such as oblique and right circular

cylinders, cones and truncated cones. Studying composed objects requires new modelling tools, different from those used previously for simpler objects (see, e.g., [20–22] and the references therein). In this paper, the phi-function approach is applied to represent analytically containment and non-overlapping conditions. Using the concept of quasi phi-functions [38], an exact mathematical model is formulated and a corresponding nonlinear programming problem is stated. A solution algorithm is proposed and computational results are presented to illustrate the approach. To the best of our knowledge, exact mathematical models for packing complex objects composed of a mixture of oblique and right convex shapes nether were considered before.

Nonspherical particles presented by the superquadric equations are widely used in different industrial production, and significantly affect the macro- and microcharacteristics of granular materials, see, e.g., [39] and the references therein. However, the particle shapes constructed by the superquadric equations are geometrically symmetrical and strictly convex, which significantly limits their further engineering applications [40]. In recent years, the composed element method has been successfully used. In this approach a complex nonconvex object is composed of basic convex elements, e.g., spheres, cylinders, super-quadratic elements and other convex (irregular) shapes [41,42].

Another source of packing complex composed shapes is additive manufacturing (AM), also known as 3D printing. AM refers to technologies for producing complex parts in a layer-by-layer material deposition process. The process takes place inside the machine in an enclosed build container or a "build volume". AM does not use any conventional physical tooling such as moulds, cutting implements or dies. Using AM, products previously designed and manufactured as assemblies of multiple components can now be manufactured as single items [43]. As a parallel manufacturing process, AM permits producing various complex parts in a single build volume. This gives rise to a build volume packing problem arising during the machine setup process. Thus, packing complex objects composed of different shapes plays an important role 3D printing.

Our interest in studying oblique objects is motivated by modelling particulate systems of nonspherical shapes [44,45] and by build volume packing problems in 3D printing [10].

The main contributions of the paper are as follows:


The paper is organized as follows. Section 2 provides the general formulation for the irregular packing problem. Quasi-phi-function s for the composed 3D objects are defined to describe analytically placement constraints in Section 3. The nonlinear programming model for the irregular packing problem is presented and solved by an algorithm described in Section 4. Computational results are given in Section 5, while Section 6 concludes. Definitions of the phi-functions and quasi-phi-function s are provided in Appendix A.

#### **2. Problem Formulation**

The packing problem is considered in the following setting. Denote a cuboid of variable length *l*, width *<sup>w</sup>* and height *<sup>h</sup>* by <sup>Ω</sup> (see Figure 1). Let a set of objects *Tq* <sup>⊂</sup> <sup>R</sup>3, *<sup>q</sup>* <sup>∈</sup> *JN* <sup>=</sup> {1, 2, ... , *<sup>N</sup>*} be given.

**Figure 1.** Container Ω.

The location and orientation of each 3D object *Tq* is defined by a vector *uq* = (*vq*, θ*q*) of its (variable) placement parameters in the fixed coordinate system *OXYZ*. Here *vq* = (*xq*, *yq*, *zq*) is a translation vector and θ*<sup>q</sup>* = (θ<sup>1</sup> *<sup>q</sup>*, θ<sup>2</sup> *<sup>q</sup>*, θ<sup>3</sup> *<sup>q</sup>* ) is a vector of rotation parameters, where θ<sup>1</sup> *<sup>q</sup>*, θ<sup>2</sup> *<sup>q</sup>*, θ<sup>3</sup> *<sup>q</sup>* are Euler angles.

The notation *Tq*(*uq*) = {*<sup>p</sup>* <sup>∈</sup> <sup>R</sup><sup>3</sup> : *p* = *vq* + Θ(θ*q*) · (*p*) *<sup>T</sup>*, <sup>∀</sup>*<sup>p</sup>* <sup>∈</sup> *Tq*} is used for translated and rotated object *Tq*, where Θ(θ*q*) = Θ(θ<sup>1</sup> *<sup>q</sup>*, θ<sup>2</sup> *<sup>q</sup>*, θ<sup>3</sup> *<sup>q</sup>* ) is a rotation matrix of the form:

$$\Theta(\boldsymbol{\theta}\_{q}) = \begin{pmatrix} \cos\theta\_{q}^{1}\cos\theta\_{q}^{3} - \sin\theta\_{q}^{1}\cos\theta\_{q}^{2}\sin\theta\_{q}^{3} & -\cos\theta\_{q}^{1}\sin\theta\_{q}^{3} - \sin\theta\_{q}^{1}\cos\theta\_{q}^{2}\cos\theta\_{q}^{3} & \sin\theta\_{q}^{1}\sin\theta\_{q}^{2} \\ \sin\theta\_{q}^{1}\cos\theta\_{q}^{3} + \cos\theta\_{q}^{1}\cos\theta\_{q}^{2}\sin\theta\_{q}^{3} & -\sin\theta\_{q}^{1}\sin\theta\_{q}^{3} + \cos\theta\_{q}^{1}\cos\theta\_{q}^{2}\cos\theta\_{q}^{3} & -\cos\theta\_{q}^{1}\sin\theta\_{q}^{2} \\ \sin\theta\_{q}^{2}\sin\theta\_{q}^{3} & \sin\theta\_{q}^{2}\cos\theta\_{q}^{3} & \cos\theta\_{q}^{2} \end{pmatrix}$$

.

Assume that *Tq*(*uq*) = *nq* ∪ *i*=1 *Tq <sup>i</sup>* (*uq*), where *<sup>T</sup><sup>q</sup> <sup>i</sup>* denotes a basic object from a family of oblique and right circular cylinders, cones, truncated cones denoted by and spheres (Figure 2a). An object *Tq* for *nq* ≥ 2 is referred to the composed object (Figure 2b).

**Figure 2.** Placement objects: (**a**) basic; (**b**) composed.

Let a sphere be defined by its centre *p q <sup>i</sup>* = (*x q i* , *y q i* , *z q i* ) and a radius *r q i* . Each basic object from the family is defined by three vectors *p q <sup>i</sup>*<sup>1</sup> = (*x q <sup>i</sup>*1, *y q <sup>i</sup>*1, *z q <sup>i</sup>*1), *p q <sup>i</sup>*<sup>2</sup> = (*x q <sup>i</sup>*2, *y q <sup>i</sup>*2, *z q <sup>i</sup>*2) and **<sup>n</sup>***<sup>q</sup> <sup>i</sup>* = (n*<sup>x</sup> <sup>i</sup>* , <sup>n</sup>*<sup>y</sup> <sup>i</sup>* , <sup>n</sup>*<sup>z</sup> i* ), as well as a pair of parameters *r q <sup>i</sup>*<sup>1</sup> and *r q <sup>i</sup>*2. Here *p q <sup>i</sup>*1, *p q <sup>i</sup>*<sup>2</sup> are the centres and *r q <sup>i</sup>*1, *r q <sup>i</sup>*<sup>2</sup> are the radii of the bottom and top bases of *T<sup>q</sup> <sup>i</sup>* , **<sup>n</sup>***<sup>q</sup> <sup>i</sup>* is the unit normal vector to the bottom (top) base of *<sup>T</sup><sup>q</sup> i* .

Note that if *r q <sup>i</sup>*<sup>1</sup> = *r q <sup>i</sup>*<sup>2</sup> <sup>&</sup>gt; 0 then *<sup>T</sup><sup>q</sup> <sup>i</sup>* is a circular cylinder; if *r q <sup>i</sup>*<sup>1</sup> *r q <sup>i</sup>*<sup>2</sup> and *r q <sup>i</sup>*<sup>1</sup> > 0,*r q <sup>i</sup>*<sup>2</sup> <sup>&</sup>gt; 0 then *<sup>T</sup><sup>q</sup> <sup>i</sup>* is a circular truncated cone; if *r q <sup>i</sup>*<sup>1</sup> > 0,*r q <sup>i</sup>*<sup>2</sup> = 0 or *r q <sup>i</sup>*<sup>1</sup> = 0,*r q <sup>i</sup>*<sup>2</sup> <sup>&</sup>gt; 0 then *<sup>T</sup><sup>q</sup> <sup>i</sup>* is a circular cone. The height of each object *T<sup>q</sup> <sup>i</sup>* ∈ is denoted by *h q i* .

Let *p q <sup>i</sup>* = (*x q i* , *y q i* , *z q i* ) be a reference point of the basic object *T<sup>q</sup> <sup>i</sup>* ⊂ *Tq*: the centre point of a sphere or the central point for a circular base of the objects from the family . In what follows, notation *p q <sup>i</sup>* = (*x q i* , *y q i* ,*z q i* ) = *vq* + Θ(θ*q*)·(*p q i* ) *<sup>T</sup>* is used, where *vq* is the translation vector and <sup>Θ</sup>(θ*q*) is the rotation matrix of the object *Tq*.

The *placement constraints* can be stated in the following form:

$$\text{dist}T\_{\mathfrak{q}}(\mathfrak{u}\_{\mathfrak{q}}) \cap \text{int}T\_{\mathfrak{F}}(\mathfrak{u}\_{\mathfrak{F}}) = O \text{ for } \mathfrak{q} > \mathfrak{g} \in \mathfrak{f}\_{\mathbf{N}\_{\prime}} \tag{1}$$

$$T\_{\emptyset}(u\_{\emptyset}) \subset \Omega \text{ for } q \in \mathfrak{h}. \tag{2}$$

Conditions (1) describe non-overlapping for all pairs of objects *Tq*(*uq*) and *Tg*(*ug*) for *q* > *g* ∈ *JN* (further, *non-overlapping constraints*), while conditions (2) assure containing *Tq*(*uq*) in the container Ω for *q* ∈ *JN* (further, containment constraints).

*Irregular packing problem*. Pack the set of objects *Tq*, *q* ∈ *JN*, within a cuboidal container Ω of minimal volume κ = *l* · *w* · *h*, taking into account the placement constraints (1)–(2).

#### **3. Analytical Tools**

In this section geometric tools for the mathematical modelling of the placement constraints (1) and (2) are presented. To describe analytically the relations between a pair of objects considered in the placement constraints, the phi-functions [38] and quasi-phi-functions [18] are used (see Appendix A for the details. These functions for irregular objects composed of oblique and right shapes are introduced in this paper for the first time.

#### *3.1. Modeling Non-Overlapping Constraints*

Consider a pair of the composed objects *Tq*(*uq*) = *nq* ∪ *i*=1 *Tq <sup>i</sup>* (*uq*) and *Tg*(*ug*) = *ng* ∪ *j*=1 *Tg <sup>j</sup>* (*ug*).

To describe non-overlapping constraints (1) a phi-function for two composed objects *Tq*(*uq*) and *Tg*(*ug*) is introduced. It can be written in the form

$$\boldsymbol{\Phi}'\_{\boldsymbol{q}\mathcal{K}}(\boldsymbol{u}\_{\boldsymbol{q}\prime},\boldsymbol{u}\_{\boldsymbol{\mathcal{K}}\prime}\boldsymbol{\pi}\_{\boldsymbol{q}\mathcal{K}}) = \min\_{\boldsymbol{i}=1,\ldots,n\_{\boldsymbol{q}\prime}\in\mathcal{I},\ldots,n\_{\boldsymbol{\mathcal{K}}}} \boldsymbol{\Phi}'^{\boldsymbol{q}\mathcal{K}}\_{\boldsymbol{i}\mathcal{j}}(\boldsymbol{u}\_{\boldsymbol{q}\prime},\boldsymbol{u}\_{\boldsymbol{\mathcal{K}}\prime}\boldsymbol{\pi}^{\boldsymbol{q}\mathcal{K}}\_{\boldsymbol{i}\mathcal{j}}) \tag{3}$$

where τ*qg* = (τ *qg ij* , *<sup>i</sup>* <sup>=</sup> 1, ... , *nq*, *<sup>j</sup>* <sup>=</sup> 1, ... , *nq*), <sup>Φ</sup> *qg ij* (*uq*, *ug*, τ *qg ij* ) is an adjusted quasi-phi-function for convex objects *T<sup>q</sup> <sup>i</sup>* (*uq*) and *<sup>T</sup><sup>g</sup> <sup>j</sup>* (*ug*).

Let *Pqg ij* = {(*x*, *y*, *z*) : *x* ≤ 0} be a half space. Denote the half-space *P* translated along OX on value μ*qg ij* and rotated by angles <sup>θ</sup>1, <sup>θ</sup>2, by *<sup>P</sup> qg ij* <sup>=</sup> {*<sup>p</sup>* : <sup>Ψ</sup>*qg ij* (*p*, τ *qg ij* ) ≤ 0}, where

$$\Psi\_{ij}^{q\overline{g}}(p,\pi\_{ij}^{q\overline{g}}) = \cos\theta\_{ij}^{1q\overline{g}} \cdot \cos\theta\_{ij}^{2q\overline{g}} \cdot \mathbf{x} - \sin\theta\_{ij}^{2q\overline{g}} \cdot y + \sin\theta\_{ij}^{1q\overline{g}} \cdot \cos\theta\_{ij}^{2q\overline{g}} \cdot z + \mu\_{ij}^{q\overline{g}} \tag{4}$$

*p* = (*x*, *y*, *z*), τ *qg ij* = (θ 1*qg ij* , θ 2*qg ij* , <sup>μ</sup>*qg ij* ), θ 1*qg ij* and θ 2*qg ij* are rotation angles of the half space *Pqg ij* around the axes *OY* and *OZ* in the fixed coordinate system. Define the plane-*L qg ij* <sup>=</sup> {(*x*, *<sup>y</sup>*, *<sup>z</sup>*) : <sup>Ψ</sup>*qg ij* (*p*, τ) = 0}.

*A quasi-phi-function for convex basic objects T<sup>q</sup> <sup>i</sup>* (*uq*) *and T<sup>g</sup> <sup>j</sup>* (*ug*) can be defined in the form

$$\Phi^{r\emptyset\emptyset}\_{\vec{i}\vec{j}}(\boldsymbol{u}\_{\boldsymbol{q}\prime}\boldsymbol{u}\_{\boldsymbol{\mathcal{K}}\prime}\boldsymbol{\tau}^{\otimes\boldsymbol{\mathcal{K}}}\_{\vec{i}\vec{j}}) = \min\{\Phi^{q}\_{\vec{i}}(\boldsymbol{u}\_{\boldsymbol{q}\prime}\boldsymbol{\tau}^{\otimes\boldsymbol{\mathcal{K}}}\_{\vec{i}\vec{j}}), \Phi^{\ast\boldsymbol{\mathcal{K}}}\_{\vec{j}}(\boldsymbol{u}\_{\boldsymbol{\mathcal{K}}\prime}\boldsymbol{\tau}^{\otimes\boldsymbol{\mathcal{K}}}\_{\vec{i}\vec{j}})\}\tag{5}$$

where Φ*<sup>q</sup> i* (*uq*, τ *qg ij* ) is a phi-function for the object *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) and the half space *P qg ij* , <sup>Φ</sup>∗*<sup>g</sup> <sup>j</sup>* (*ug*, τ *qg ij* ) is a phi-function for the object *T<sup>g</sup> <sup>j</sup>* (*ug*) and the half space *P* -∗*qg ij* <sup>=</sup> <sup>R</sup>3\int*<sup>P</sup> qg ij* . Here τ *qg ij* = (θ 1*qg ij* , θ 2*qg ij* , <sup>μ</sup>*qg ij* ) is the vector of auxiliary variables of the quasi-phi-function Φ *qg ij* .

Therefore to define the *quasi-phi-functions* (5) *for each pair of convex basic objects* the phi-functions Φ*q i* (*uq*, τ *qg ij* ) and <sup>Φ</sup>∗*<sup>g</sup> <sup>j</sup>* (*ug*, τ *qg ij* ) have to be derived.

First, define the phi-function Φ*<sup>q</sup> i* (*uq*, τ *qg ij* ) for the basic object *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) and the half space *P qg ij* . Let *T<sup>q</sup> <sup>i</sup>* (*uq*) be a sphere centred at the point *p q <sup>i</sup>*<sup>1</sup> = (*x q i* , *y q i* ,*z q i* ) and having its radius *r q i* . The phi-function for the sphere *T<sup>q</sup> <sup>i</sup>* (*uq*) and a half space *P qg ij* has the form

$$\Phi\_i^q(\boldsymbol{u}\_{\boldsymbol{\eta}}, \boldsymbol{\tau}\_{ij}^{q\boldsymbol{\eta}}) = \cos\theta\_{ij}^{1q\boldsymbol{\eta}} \cdot \cos\theta\_{ij}^{2q\boldsymbol{\eta}} \cdot \overrightarrow{\mathbf{x}}\_{i1}^q - \sin\theta\_{ij}^{2q\boldsymbol{\eta}} \cdot \overrightarrow{\boldsymbol{y}}\_{i1}^q + \sin\theta\_{ij}^{1q\boldsymbol{\eta}} \cdot \cos\theta\_{ij}^{2q\boldsymbol{\eta}} \cdot \overrightarrow{\boldsymbol{z}}\_{i1}^q + \boldsymbol{\mu}\_{ij}^{q\boldsymbol{\eta}} - \boldsymbol{r}\_i^q = 0$$

while the phi-function for the sphere *T<sup>g</sup> <sup>j</sup>* (*ug*) and the half space *P* -∗*qg ij* can be defined as follows

$$\Phi\_{\vec{j}}^{\ast \mathsf{S}}(\boldsymbol{u}\_{\mathcal{S}},\boldsymbol{\mathsf{T}}\_{\vec{i}\boldsymbol{j}}^{\otimes \mathsf{S}}) = -\cos\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{1g\underline{\mathsf{g}}}} \cdot \cos\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{2g\underline{\mathsf{g}}}} \cdot \overrightarrow{\mathbf{x}}\_{\vec{i}1}^{\mathsf{q}} + \sin\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{2g\underline{\mathsf{g}}}} \cdot \overrightarrow{\boldsymbol{y}}\_{\vec{i}1}^{\mathsf{q}} - \sin\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{1g\underline{\mathsf{g}}}} \cdot \cos\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{2g\underline{\mathsf{g}}}} \cdot \overrightarrow{\boldsymbol{z}}\_{\vec{i}1}^{\mathsf{q}} - \boldsymbol{\mu}\_{\vec{i}\boldsymbol{j}}^{\mathrm{q\underline{\mathsf{g}}}} - r\_{\vec{i}}^{\mathsf{q}} = -\sin\theta\_{\vec{i}\boldsymbol{j}}^{\mathrm{1g\underline{\mathsf{g}}}}.$$

For *T<sup>q</sup> <sup>i</sup>* (*uq*) from the family , let the bottom base of *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) is a circle centred at the point *p q <sup>i</sup>*<sup>1</sup> = (*x q i* , *y q i* ,*z q i* ) and having its radius *r q <sup>i</sup>*1, while the top circular base of *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) is centred at *p q <sup>i</sup>*<sup>2</sup> = (*x q i* , *y q i* ,*z q <sup>i</sup>* + *h q i* ) and has its radius *r q <sup>i</sup>*<sup>2</sup> (the cone corresponds to *r q <sup>i</sup>*<sup>2</sup> = 0).

The phi-function for the object *T<sup>q</sup> <sup>i</sup>* (*uq*) ∈ and a half space *P qg ij* has the form

$$\begin{aligned} \Phi\_i^{\mathfrak{q}}(\boldsymbol{\upmu\_{\boldsymbol{q}}}, \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}}) &= \min \{ f\_1(\boldsymbol{\upmu\_{\boldsymbol{q}}}, \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}}), f\_2(\boldsymbol{\upmu\_{\boldsymbol{q}}}, \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}}) \}, \\ f\_1(\boldsymbol{\upmu\_{\boldsymbol{q}}}, \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}}) &= \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}} \cdot \boldsymbol{\uppi\_{i1}^{\mathfrak{q}}} + \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}} - r\_{i1}^{\mathfrak{q}} \sqrt{1 - \left( \overleftarrow{\mathbf{n}}\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}} \cdot \overrightarrow{\mathbf{n}}\_{i}^{\mathfrak{q}} \right)^{2}}, \\ f\_2(\boldsymbol{\upmu\_{\boldsymbol{q}}}, \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}}) &= \overleftarrow{\mathbf{n}\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}} \cdot \overrightarrow{p}\_{i2}^{\mathfrak{q}} + \boldsymbol{\uppi\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}}} - r\_{i2}^{\mathfrak{q}} \sqrt{1 - \left( \overleftarrow{\mathbf{n}}\_{ij}^{\mathfrak{q}\boldsymbol{\mathcal{S}}} \cdot \overrightarrow{\mathbf{n}}\_{i}^{\mathfrak{q}} \right)^{2}}. \end{aligned}$$

where **n***qg ij* = (**n***x ij*,**n***y ij*,**n***z ij*) denotes a unit vector of the external normal to the half space *P qg ij* and **n***q <sup>i</sup>* = (**n***x <sup>i</sup>* ,**n***y <sup>i</sup>* ,**n***z <sup>i</sup>*) stands for a unit vector of the external normal to the object *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) (see Figure 3).

The phi-function for the object *T<sup>g</sup> <sup>j</sup>* (*ug*) and a half space *P* -∗*qg ij* has the form

$$\begin{split} \Phi\_{j}^{\ast \mathcal{S}}(\boldsymbol{u}\_{\mathcal{S}'}\boldsymbol{\tau}\_{ij}^{q\mathcal{S}}) &= \min \{ f\_1(\boldsymbol{u}\_{\mathcal{S}'}\boldsymbol{\tau}\_{ij}^{q\mathcal{S}}), f\_2(\boldsymbol{u}\_{\mathcal{S}'}\boldsymbol{\tau}\_{ij}^{q\mathcal{S}}) \}, \\ \{ f\_1(\boldsymbol{u}\_{\mathcal{S}'}\boldsymbol{\tau}\_{ij}^{q\mathcal{S}}) &= -\overrightarrow{\mathbf{n}}\_{ij}^{q\mathcal{S}} \cdot \overrightarrow{p}\_{j1}^{\mathcal{S}} - \mu\_{ij}^{q\mathcal{S}} - r\_{j1}^{\mathcal{S}} \sqrt{1 - \left(\overrightarrow{\mathbf{n}}\_{ij}^{q\mathcal{S}} \cdot \overrightarrow{\mathbf{n}}\_{i}^{q}\right)^2}, \\ \{ f\_2(\boldsymbol{u}\_{\mathcal{S}'}\boldsymbol{\tau}\_{ij}^{q\mathcal{S}}) &= -\overrightarrow{\mathbf{n}}\_{ij}^{q\mathcal{S}} \cdot \overrightarrow{p}\_{j2}^{\mathcal{S}} - \mu\_{ij}^{q\mathcal{S}} - r\_{j2}^{\mathcal{S}} \sqrt{1 - \left(\overrightarrow{\mathbf{n}}\_{ij}^{q\mathcal{S}} \cdot \overrightarrow{\mathbf{n}}\_{i}^{q}\right)^2}. \end{split}$$

**Figure 3.** Non-overlapping: (**a**) convex basic objects *T<sup>q</sup> <sup>i</sup>* and the half space *P qg ij* ; (**b**) convex basic objects *Tq <sup>i</sup>* and the half space *P* -∗*qg ij* .

The inequality Φ *qg*(*uq*, *ug*, τ*qg*) ≥ 0 assures the non-overlapping condition int *Tq*(*uq*) ∩ int *Tg*(*ug*) = ∅ in (1).

#### *3.2. Modeling the Containment Constraints*

Let us express the containment constraint *Tq*(*uq*) <sup>⊂</sup> <sup>Ω</sup> in the equivalent form as int*Tq*(*uq*) <sup>∩</sup> <sup>Ω</sup><sup>∗</sup> <sup>=</sup> <sup>∅</sup>, where <sup>Ω</sup><sup>∗</sup> = *<sup>R</sup>*3\intΩ.

The phi-function for the objects *Tq*(*uq*) = *nq* ∪ *i*=1 *Tq <sup>i</sup>* (*uq*) and Ω<sup>∗</sup> can be defined in the form:

$$\Phi\_{\emptyset}(\mu\_{\emptyset}) = \min\_{i=1,\ldots,n\_{\emptyset}} \Phi\_{i}^{q}(\mu\_{\emptyset}),\tag{6}$$

where Φ*<sup>q</sup> i* (*uq*) is the phi-function for the basic convex object *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) and the object Ω<sup>∗</sup> .

The inequality <sup>Φ</sup>*q*(*uq*) <sup>≥</sup> 0 implies fulfilling the condition int*Tq*(*uq*) <sup>∩</sup> <sup>Ω</sup><sup>∗</sup> <sup>=</sup> <sup>∅</sup> that describes the containment condition (2).

*Containment of a sphere T<sup>q</sup> <sup>i</sup>* (*uq*) *in a cuboid* <sup>Ω</sup>. The phi-function for the sphere *<sup>T</sup><sup>q</sup> <sup>i</sup>* (*uq*) and the object Ω∗ can be defined by the equation

$$\Phi\_i^q(u\_q) = \min \{ l - \overrightarrow{p}\_{i1}^{qx} - r\_{i1'}^q, w - \overrightarrow{p}\_{i1}^{qy} - r\_{i1'}^q, h - \overrightarrow{p}\_{i1}^{qz} - r\_{i1'}^q, \overrightarrow{p}\_{i1}^{qx} - r\_{i1'}^q, \overrightarrow{p}\_{i1}^{qy} - r\_{i1'}^q, \overrightarrow{p}\_{i1}^{qz} - r\_{i1'}^q \}.$$

*Containment of the object T<sup>q</sup> <sup>i</sup>* (*uq*) *from the family in a cuboid* Ω. The phi-function for the objects *Tq <sup>i</sup>* (*uq*) and Ω<sup>∗</sup> has the form

$$\Phi\_i^q(u\_q) = \min\{\varphi\_k(u\_q), k = 1, \dots, 6\} \tag{7}$$

where

$$\begin{split} \varphi\_{1}(u\_{\emptyset}) &= \min(\overline{p}\_{i1}^{qx} - r\_{i1}^{q}\sqrt{1 - \left(\overline{n}\right)^{qx}})^{2}, \overline{p}\_{i2}^{qx} - r\_{i2}^{q}\sqrt{1 - \left(\overline{n}\right)^{qx}})^{2}), \\ \varphi\_{2}(u\_{\emptyset}) &= \min(\overline{p}\_{i1}^{qy} - r\_{i1}^{q}\sqrt{1 - \left(\overline{n}\right)^{qy}})^{2}, \overline{p}\_{i2}^{qy} - r\_{i2}^{q}\sqrt{1 - \left(\overline{n}\right)^{qy}})^{2}), \\ \varphi\_{3}(u\_{\emptyset}) &= \min(\overline{p}\_{i1}^{qz} - r\_{i1}^{q}\sqrt{1 - \left(\overline{n}\right)^{qz}})^{2}, \overline{p}\_{i2}^{qz} - r\_{i2}^{q}\sqrt{1 - \left(\overline{n}\right)^{qz}})^{2}), \\ \varphi\_{4}(u\_{\emptyset}) &= \min(l - \overline{p}\_{i1}^{qx} - r\_{i1}^{q}\sqrt{1 - \left(\overline{n}\right)^{qx}})^{2}, l - \overline{p}\_{i2}^{qx} - r\_{i2}^{q}\sqrt{1 - \left(\overline{n}\right)^{qx}})^{2}), \end{split}$$

$$\begin{split} \varphi\_{5}(u\_{q}) &= \min \{ w - \overrightarrow{p}\_{i1}^{qy} - r\_{i1}^{q} \sqrt{1 - \left(\overrightarrow{\overline{n}\_{i}^{qy}}\right)^{2}}, w - \overrightarrow{p}\_{i2}^{qy} - r\_{i2}^{q} \sqrt{1 - \left(\overrightarrow{\overline{n}\_{i}^{qy}}\right)^{2}} \}, \\ \varphi\_{6}(u\_{q}) &= \min \{ h - \overrightarrow{p}\_{i1}^{qz} - r\_{i1}^{q} \sqrt{1 - \left(\overrightarrow{\overline{n}\_{i}^{qz}}\right)^{2}}, h - \overrightarrow{p}\_{i2}^{qz} - r\_{i2}^{q} \sqrt{1 - \left(\overrightarrow{\overline{n}\_{i}^{qz}}\right)^{2}} \}. \end{split}$$

The inequality Φ*<sup>q</sup> i* (*uq*) <sup>≥</sup> 0 implies the condition int*T<sup>q</sup> <sup>i</sup>* (*uq*) <sup>∩</sup> <sup>Ω</sup><sup>∗</sup> <sup>=</sup> <sup>∅</sup>.

#### **4. Mathematical Model and Solution Algorithm**

#### *4.1. Mathematical Model*

Now the irregular packing problem can be formulated as the following nonlinear optimization problem:

$$\min \kappa \text{ s.t. } (u, \pi) \in \mathcal{W}, \tag{8}$$

$$\mathcal{W} = \left\{ (\mathfrak{u}, \mathfrak{r}) : \Phi'\_{\mathfrak{g}\mathbb{K}}(\mathfrak{u}\_{\mathfrak{q}}, \mathfrak{u}\_{\mathfrak{g}\mathbb{K}}, \mathfrak{r}\_{\mathfrak{g}\mathbb{K}}) \ge 0, \, q > \mathfrak{g} \in \mathbf{I}\_{\mathbb{N}}, \, \Phi\_{\mathfrak{q}}(\mathfrak{u}\_{\mathfrak{q}}) \ge 0, \, q \in \mathbf{I}\_{\mathbb{N}} \right\}, \tag{9}$$

where τ*qg* = (τ *qg ij* , *i* = 1, ... , *nq*, *j* = 1, ... , *nq*) is a vector of the auxiliary variables for the quasi-phi-function Φ *qg*(*uq*, *ug*, τ*qg*) (3)–(5) for the objects *Tq*(*uq*) and *Tg*(*ug*), Φ*q*(*uq*) is the phi-function (6)–(7) for the objects *Tq*(*uq*) and Ω<sup>∗</sup> , *q* ∈ **I***N*.

The feasible region *W* in (9) is defined by a system of nonsmooth inequalities that can be reduced to a system of inequalities with differentiable functions.

The model (8)–(9) is a nonconvex and continuous nonlinear programming problem. This is an exact formulation in the sense that it gives all optimal solutions to the irregular packing problem.

The number of the problem variables is σ = 3(1 + 2*N* + *m*), where *m* = *q*>*g*∈**I***<sup>N</sup> nqng*. The model

(8)–(9) involves *O*(*N*2) nonlinear inequalities and *O*(*N*2) variables.

The model (8)–(9) represents all globally optimal solutions to the original irregular packing problem. It can be solved by any available global solver, e.g., BARON or LGO included in AMPL [46]. However, due to large number of variables and constraints, a direct solution to this problem may be time consuming and complicated. In the next section a solution approach is proposed to search for a local minimum of the problem (8)–(9). This solution can be used either as a reasonable approximation to the original global solution, or as a starting point for a global solver or heuristics.

#### *4.2. Solution Algorithm*

The following multistart strategy is used to solve the problem (8)–(9). A number of feasible starting points is generated. Then a local maximum of the problem (8)–(9) is obtained starting from each feasible point generated at the first stage. Finally, the best solution is selected from those obtained at the second stage. This result is considered as the solution of the problem (8)–(9).

#### 4.2.1. Feasible Starting Points

To find feasible starting points for the problem (8)–(9) an algorithm based on the homothetic (scaling) transformation of objects is applied. The basic steps of the algorithm are as follows.

**Step 1.** Circumscribe spheres *Sq*, *q* ∈ *JN* around the objects *Tq*(*uq*), *q* ∈ *JN*.

**Step 2**. Construct a container Ω<sup>0</sup> with sufficiently large starting length *l* 0, width *w*<sup>0</sup> and height *h*<sup>0</sup> allowing placement of all the spheres *Sq*, *q* ∈ *JN*.

**Step 3**. Generate a set of *N* randomly chosen centre points (*x*<sup>0</sup> *<sup>q</sup>*, *y*<sup>0</sup> *<sup>q</sup>*, *z*<sup>0</sup> *<sup>q</sup>* ) for the spheres *Sq*, *q* ∈ *JN* inside the container Ω0.

**Step 4**. Grow up the spheres *Sq* of the radii λ*rq*, *q* ∈ *JN*, starting from λ = 0 to the full size (λ = 1). Here the decision variables are the centres of *Sq* and the homothetic coefficient λ, where 0 ≤ λ ≤ 1.

**Step 5**. Form a vector of feasible parameters of our objects *Tq*(*uq*) ⊂ *Sq*(*vq*).

Now proceed with a more detailed description of the algorithm.

First, fix *l* = *l* 0, *w* = *w*0, *h* = *h*<sup>0</sup> and then starting from the point ω<sup>0</sup> = (*x*<sup>0</sup> <sup>1</sup>, *<sup>y</sup>*<sup>0</sup> <sup>1</sup>, *<sup>z</sup>*<sup>0</sup> <sup>1</sup>, ... , *<sup>x</sup>*<sup>0</sup> *<sup>N</sup>*, *<sup>y</sup>*<sup>0</sup> *<sup>N</sup>*, *<sup>z</sup>*<sup>0</sup> *<sup>N</sup>*, <sup>λ</sup><sup>0</sup> <sup>=</sup> <sup>0</sup>) solve the following nonlinear programming subproblem:

$$\max \lambda \text{ s.t. } \omega \in W\_{\lambda} \subset \mathbb{R}^{3N+1}\_{\text{-}}$$

$$\mathcal{W}\_{\lambda} = \langle \omega : \Phi^{\mathbb{S}\_{\mathbb{P}} \mathbb{S}\_{\mathbb{S}}}(\omega) \ge 0, \Phi^{\mathbb{S}\_{\mathbb{P}} \Omega^\*}(\omega) \ge 0, q < \emptyset \in \mathbf{I}\_{\mathbb{N}}, 1 - \lambda \ge 0, \lambda \ge 0 \rangle.$$

Here ω = (*x*1, *y*1, *z*1, ... , *xN*, *yN*, *zN*, λ),

$$\Phi^{\mathbb{S}\_{\mathbb{S}}\mathbb{S}\_{\mathbb{S}}}(\omega) = (\mathbf{x}\_{\mathfrak{q}} - \mathbf{x}\_{\mathfrak{J}})^2 + (y\_{\mathfrak{q}} - y\_{\mathfrak{J}})^2 + (z\_{\mathfrak{q}} - z\_{\mathfrak{J}})^2 - (\lambda r\_{\mathfrak{q}} + \lambda r\_{\mathfrak{J}})^2$$

is the phi-function for the sphere *Sq* of the radius λ*rq* and the sphere *Sg* of the radius λ*rg*,

$$\Phi^{S\_{\emptyset}\Omega^\*}(\omega) = \min \{ \varphi(\omega)\_{k\emptyset}, k = 1, \dots, 6 \}$$

is the phi-function for the sphere *Sq* of the radius λ*rq* and the object Ω<sup>∗</sup> , where

$$\begin{aligned} \varphi\_{1q}(\omega) &= l^0 - \chi\_{\overline{q}} - \lambda r\_{\overline{q}\prime} \, \, \varphi\_{2q}(\omega) \, = \mathfrak{x}\_{\overline{q}} - \lambda r\_{\overline{q}\prime}, \\\\ \varphi\_{3q}(\omega) &= \mathfrak{w}^0 - \mathcal{y}\_{\overline{q}} - \lambda r\_{\overline{q}\prime} \, \, \varphi\_{4q}(\omega) \, = \mathcal{y}\_{\overline{q}} - \lambda r\_{\overline{q}\prime}, \\ \varphi\_{5q}(\omega) &= l^0 - \mathcal{z}\_{\overline{q}} - \lambda r\_{\overline{q}\prime} \, \, \varphi\_{6q}(\omega) \, = \mathcal{z}\_{\overline{q}} - \lambda r\_{\overline{q}\prime}. \end{aligned}$$

Denote the global maximum point of the above subproblem by (*x*∗ <sup>1</sup>, *y*<sup>∗</sup> <sup>1</sup>, *z*<sup>∗</sup> <sup>1</sup>, ... , *x*<sup>∗</sup> *<sup>N</sup>*, *y*<sup>∗</sup> *<sup>N</sup>*, *z*<sup>∗</sup> *<sup>N</sup>*, λ<sup>∗</sup> = 1) and form a vector of feasible parameters ς0=(*l* 0, *w*0, *h*0, *u*<sup>0</sup> <sup>1</sup>, ... , *<sup>u</sup>*<sup>0</sup> *<sup>q</sup>*, ... , *u*<sup>0</sup> *<sup>N</sup>*). Here *<sup>u</sup>*<sup>0</sup> *<sup>q</sup>* = (*x*<sup>∗</sup> *<sup>q</sup>*, *y*<sup>∗</sup> *<sup>q</sup>*, *z*<sup>∗</sup> *<sup>q</sup>*, θ<sup>0</sup> *<sup>q</sup>* ) and θ0 *<sup>q</sup>* is a vector of randomly generated rotation parameters of the objects *Tq*, *q* ∈ *JN*.

To generate a starting point *u*<sup>0</sup> = (ς0, τ0) for a subsequent search for a local minimum of the problem (8)–(9), define a vector τ<sup>0</sup> = (τ *qg ij* , *i* = 1, ... , *nq*, *j* = 1, ... , *nq*, *q* > *g* ∈ **I***N*) by solving the following optimization subproblems: max τ*qg ij* Φ *qg ij* (*uq*, *ug*, τ *qg ij* ) for *i* = 1, ... , *nq*, *j* = 1, ... , *ng*, *q* > *g* ∈ **I***N*.

#### 4.2.2. Local Optimization

The definition of the feasible set *W* in (9) involves a large number *O*(*N*2) of inequalities and variables. To cope with this large-scale problem the decomposition algorithm [47] is used that reduces the problem (8)–(9) to a sequence of nonlinear programming subproblems with a smaller number *O*(*N*) of inequalities and variables. The key idea of the algorithm is as follows. For each vector of feasible placement parameters of the objects, fixed individual cubic containers are constructed containing spheres that circumscribe the appropriate convex basic object. Each sphere is allowed to move within the appropriate individual container. The motion of each sphere is described by a system of six linear inequalities. Then a subregion of the feasible region *W* is formed as follows. For all spheres, *O*(*N*) inequalities are added to the system (9) and *O*(*N*2) phi-inequalities corresponding to the pairs of basic objects with individual containers non-overlapping each other are deleted. Moreover, some redundant containment constraints are also deleted. This auxiliary local minimization subproblem has *O*(*N*) variables and nonlinear constraints. The solution to this problem is used as a starting feasible point for the next iteration. On the last iteration of the algorithm a local minimum to the problem (8)–(9) is obtained.

#### **5. Computational Results**

In this section, five new benchmark instances are provided to demonstrate the efficiency of the proposed methodology. All experiments were running on an AMD FX(tm)-6100, 3.30 GHz computer (Ultra A0313). Programming Language C++, Windows 7. For the local optimisation, the IPOPT code (https://projects.coin-or.org/Ipopt) reported in [48] was used under default options. The multistart approach was used for the problem (8)–(9) as follows. For each problem instance, 10 starting points were generated using the algorithm of SubSection 4.2.1. Then 10 corresponding local minima were obtained by the algorithm described in SubSection 4.2.2 and the best local mimimum was selected as an approximate solution to the problem (8)–(9). The CPU time indicated for each problem instance is the total time for all 10 runs.

**Example 1.** *The optimized packings of N* = 25 *irregular objects. Each object Tq*(*uq*), *q* = 1, ... , 25 *is composed by the union of two cones T<sup>q</sup> <sup>i</sup>* (*uq*), *i* = 1, 2 *given by the following parameters: p q* <sup>11</sup> = (0,0,0), *p q* <sup>12</sup> <sup>=</sup> (9,0,0), **<sup>n</sup>***<sup>q</sup>* <sup>1</sup> = (1,0,0), *r q* <sup>11</sup> = 3, *r q* <sup>12</sup> <sup>=</sup> <sup>0</sup> *and p<sup>q</sup>* <sup>21</sup> = (7,0,0), *p q* <sup>22</sup> <sup>=</sup> (−2,0,0), **<sup>n</sup>***<sup>q</sup>* <sup>2</sup> = (1,0,0), *r q* <sup>21</sup> = 3, *r q* <sup>22</sup> = 0 *respectively.*

The best local minimum obtained for 2857.99 sec is κ∗ = *l* <sup>∗</sup> · *w*<sup>∗</sup> · *h*<sup>∗</sup> = 17.673918 · 20.065788 · 24.972631 = 8856.3211208954.

The corresponding packing is shown in Figure 4a.

**Figure 4.** Optimized packings of composed objects: (**a**) Example 1; (**b**) Example 2.

**Example 2.** *The optimized packings of N* = 25 *basic and irregular objects, including:*


The best local minimum obtained for 3478.23 sec. is κ∗ = *l* <sup>∗</sup> · *w*<sup>∗</sup> · *h*<sup>∗</sup> = 14.889393 · 24.925430 · 17.165791 = 6370.6459961746.

The corresponding packing is shown in Figure 4b.

**Example 3.** *The optimized packings of N* = 2, 3, 4, 5 *objects from Example 1:*

(1) *for N* = 2 *the best local minimum*

$$\kappa^\* = l^\* \cdot w^\* \cdot h^\* = 8.085071 \cdot 10.392305 \cdot 6.000000 = 504.135155$$

*was obtained for 11.638 sec. starting from the feasible solution*

$$\kappa^0 = l^0 \cdot w^0 \cdot h^0 = 26.147623 \cdot 14.547834 \cdot 14.550066 = 5534.7182137988.1$$

*The corresponding packings are shown in Figure 5.*

(2) *for N* = 3 *the best local minimum*

$$\kappa^\* = l^\* \cdot w^\* \cdot h^\* = 12.682284 \cdot 11.050980 \cdot 6.000000 = 840.910031$$

*was obtained for 26.146 sec. starting from the feasible solution*

$$\kappa^0 = l^0 \cdot w^0 \cdot h^0 = 27.617932 \cdot 26.392891 \cdot 13.953370 = 10170.849561976.1$$

**Figure 5.** Packings of two composed objects: (**a**) the feasible starting point; (**b**) the corresponding local optimal packing.

*The corresponding packings are shown in Figure 6.*

(3) *for N* = 4 *the best local minimum*

$$\kappa^\* = l^\* \cdot w^\* \cdot h^\* = 8.091432 \cdot 13.454112 \cdot 10.099594 = 1099.4724285295$$

*was obtained for 53.602 sec. starting from the feasible solution*

κ<sup>0</sup> = *l* <sup>0</sup> · *<sup>w</sup>*<sup>0</sup> · *<sup>h</sup>*<sup>0</sup> = 14.771808 · 27.984482 · 24.620742 = 10177.75667595.

**Figure 6.** Packings of three composed objects: (**a**) the feasible starting point; (**b**) the corresponding local optimal packing.

*The corresponding packings are shown in Figure 7.*

(4) *for N* = 5 *the best local minimum*

κ∗ = *l* <sup>∗</sup> · *w*<sup>∗</sup> · *h*<sup>∗</sup> = 11.492738 · 11.267728 · 10.652737 = 1379.4979707504

*was obtained for 79.186 sec. starting from the feasible solution*

κ<sup>0</sup> = *l* <sup>0</sup> · *<sup>w</sup>*<sup>0</sup> · *<sup>h</sup>*<sup>0</sup> = 26.126900 · 26.524042 · 23.378807 = 16201.302676444.

**Figure 7.** Packings of four composed objects: (**a**) the feasible starting point; (**b**) the corresponding local optimal packing.

*The corresponding packings are shown in Figure 8.*

**Figure 8.** Packings of five composed objects: (**a**) the feasible starting point; (**b**) the corresponding local optimal packing.

**Example 4.** *The optimized packings of N* = 12 *objects (three objects of each type from Example 2). The best local minimum*

> κ∗ = *l* <sup>∗</sup> · *w*<sup>∗</sup> · *h*<sup>∗</sup> = 15.084109 · 15.073729 · 16.129947 = 3667.5268798149

*was obtained for 934.103 sec. starting from the feasible solution*

κ<sup>0</sup> = *l* <sup>0</sup> · *<sup>w</sup>*<sup>0</sup> · *<sup>h</sup>*<sup>0</sup> = 20.944159 · 24.547251 · 32.847595 = 16887.65573111.

*The corresponding packings are shown in Figure 9.*

**Figure 9.** Packings of twelve composed objects: (**a**) the feasible starting point; (**b**) the corresponding local optimal packing.

Table 1 below provides objective function values for the feasible starting point (κ0) and the corresponding local optimal solution (κ∗ ) for all 10 starting points. These values are indicated for N = 2, 3, 4, 5 (Example 3) and N = 12 (Example 4). The best local optimal solutions are highlighted in bold.


**Table 1.** Objective values for feasible starting points and corresponding local optima.

As can be seen from Table 1, different starting points lead to different local minima.

#### **6. Conclusions**

Packing irregular 3D objects in a cuboid of minimum volume is considered. Each irregular object is composed by convex shapes from the family of oblique and right circular cylinders, cones, truncated cones and spheres. Continuous translations and rotations for all objects are allowed. The optimized packing problem is formulated for 3D regular and irregular objects. New analytical tools (quasi-phi-functions and phi-functions) are defined for the first time to describe non-overlapping and containment constraints for irregular objects composed by oblique shapes.

The phi-function technique is used to state the irregular packing in the form of nonlinear programming problem. The solution approach is proposed and illustrated by the numerical examples. The problem instances were selected to demonstrate the ability of the proposed modelling techniques to work with complex objects composed by different convex shapes used in applications.

The multistart algorithm used in this paper consists of two stages: constructing a number of initial feasible solutions and using local minimization (compaction) procedure to improve starting points. A simple and fast heuristic was implemented to get a starting solution. Using more advanced heuristics to construct better starting points may result in improving the overall optimization scheme. Some results in this direction are on the way.

The model (8)–(9) provides all global solutions to the original irregular packing problem. It can be solved by global solvers, e.g., BARON or LGO included in AMPL [46]. However, due to a large number of variables and constraints, the direct solution to this problem is time consuming and this approach was not used in the paper. Instead, the decomposition technique was used for the large-scale problem (8)–(9) and combined with IPOPT for solving NLP subproblems. An interesting direction for the future research is using alternative decomposition techniques [49] or constructing Lagrangian relaxations with respect to binding constraints (see, e.g., [50] and the references therein).

**Author Contributions:** Investigation, A.P., T.R., I.L.; Methodology, A.P., T.R., I.L.; Programming Algorithms, A.P.; Project administration, T.R., I.L.; Writing—original draft, T.R., I.L.; Writing—review and editing, T.R., I.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** A. Pankratov and T. Romanova were partially supported by the "Program for the State Priority Scientific Research and Technological (Experimental) Development of the Department of Physical and Technical Problems of Energy of the National Academy of Sciences of Ukraine" (#6541230).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

To place feasibly two objects within a container, an analytical description of the relations between the pair of objects is required. In this paper the phi-function technique is used to express this relation. Let *A* and *B* be three-dimensional objects. The position of the object *A* is defined by a vector of placement parameters *uA* = (*vA*, θ*A*), where: *vA* = (*xA*, *yA*, *zA*) is a translation vector and θ*<sup>A</sup>* is a vector of rotation parametrs. The object *A*, rotated by θ*<sup>A</sup>* and translated by *vA* is denoted by *A*(*uA*). Phi-functions allow us to distinguish the following three cases: *A* and *B* are intersecting so that *A* and *B* have common interior points; *A* and *B* do not intersect, i.e., *A* and *B* do not have common points; *A* and *B* are in contact, i.e., *A* and *B* have only common frontier points.

A continuous function Φ*AB*(*uA*, *uB*) is called a phi-function of the objects *A*(*uA*) and *B*(*uB*) if the following conditions are fulfilled [38]: <sup>Φ</sup>*AB*(*uA*, *uB*) <sup>&</sup>gt; 0, for *<sup>A</sup>*(*uA*) <sup>∩</sup> *<sup>B</sup>*(*uB*) <sup>=</sup> <sup>∅</sup> (see Figure A1a); <sup>Φ</sup>*AB*(*uA*, *uB*) = 0, for int*A*(*uA*) <sup>∩</sup> int*B*(*uB*) <sup>=</sup> <sup>∅</sup> and *f rA*(*uA*) <sup>∩</sup> *f rB*(*uB*) - ∅ (see Figure A1b); <sup>Φ</sup>*AB*(*uA*, *uB*) <sup>&</sup>lt; 0, for int*A*(*uA*) <sup>∩</sup> int*B*(*uB*) -∅ (see Figure A1c).

**Figure A1.** Attributes of a phi-function: (**a**) non-overlapping, Φ*AB*(*uA*, *uB*) > 0; (**b**) touching, Φ*AB*(*uA*, *uB*) = 0; (**c**) interior overlapping, Φ*AB*(*uA*, *uB*) < 0.

Here *f rA* denotes the boundary (frontier) of the object *A*, while int*A* stands for its interior. Thus, <sup>Φ</sup>*AB*(*uA*, *uB*) <sup>≥</sup> <sup>0</sup> <sup>⇔</sup> int *<sup>A</sup>*(*uA*) <sup>∩</sup> int *<sup>B</sup>*(*uB*) <sup>=</sup> <sup>∅</sup>.

A function Φ *AB*(*uA*, *uB*, *u* ) is called a *quasi-phi-function* for two objects *A*(*uA*) and *B*(*uB*) if max *u* ∈*U*<sup>Φ</sup> *AB*(*uA*, *uB*, *<sup>u</sup>* ) is a phi-function for the objects [12]. Here *u* denotes a vector of auxiliary variables depending on the object shapes. This function is defined for all values of *uA*, *uB* and has to be continuous in all its variables.

The main property of the quasi-phi-function for two objects *A*(*uA*) and *B*(*uB*) is as follows: if Φ *AB*(*uA*, *uB*, *u* ) ≥ 0 for some *u* , then int *<sup>A</sup>*(*uA*) <sup>∩</sup> int *<sup>B</sup>*(*uB*) <sup>=</sup> <sup>∅</sup>.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
