Packing Oblique 3D Objects

: 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.


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][6][7][8][9].
Various modelling and solution approaches, exact and approximated, are known for the regular packing (see, e.g., [20][21][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][24][25][26][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][29][30][31][32][33]. An alternative methodology is using genetic algorithms, directly or in combination with the first two approaches [34][35][36][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 1.
New tools of mathematical modelling are presented to describe analytically non-overlapping and containment constraints for packing irregular 3D objects composed by the union of oblique and right basic shapes (circular cylinders, cones, truncated cones and spheres). The objects can be freely translated and rotated.

2.
An exact mathematical model for the irregular packing problem is formulated in the form of nonlinear continuous programming problem.

3.
A solution algorithm for the irregular packing problem is developed.

4.
New benchmark instances are provided to illustrate the efficiency of the approach.
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.

Problem Formulation
The packing problem is considered in the following setting. Denote a cuboid of variable length l, width w and height h by Ω (see Figure 1). Let a set of objects T q ⊂ R 3 , q ∈ J N = {1, 2, . . . , N} be given. The location and orientation of each 3D object T q is defined by a vector u q = (v q , θ q ) of its (variable) placement parameters in the fixed coordinate system OXYZ.
is a rotation matrix of the form: i denotes a basic object from a family of oblique and right circular cylinders, cones, truncated cones denoted by and spheres (Figure 2a). An object T q for n q ≥ 2 is referred to the composed object ( Figure 2b).  1  3  1  2  3  1  3  1  2  3  1  2   1  3  1  2  3  1  3  1  2  3  1  2   2  3  2  3  2   cos  cos  sin  cos  sin  cos  sin  sin  cos  cos  sin  sin   ( ) = sin  cos  cos  cos  sin  sin  sin  cos  cos  cos  cos  sin   sin  sin  sin cos cos , where q i T denotes a basic object from a family of oblique and right circular cylinders, cones, truncated cones denoted by ℑ and spheres ( Figure 2a). An object q T for 2 q n ≥ is referred to the composed object ( Figure 2b). Let a sphere be defined by its centre ( , , ) and a radius q i r . Each basic object from the family ℑ is defined by three vectors T is used, where v q is the translation vector and Θ(θ q ) is the rotation matrix of the object T q . The placement constraints can be stated in the following form: Conditions (1) describe non-overlapping for all pairs of objects T q (u q ) and T g (u g ) for q > g ∈ J N (further, non-overlapping constraints), while conditions (2) assure containing T q (u q ) in the container Ω for q ∈ J N (further, containment constraints).
Irregular packing problem. Pack the set of objects T q , q ∈ J N , within a cuboidal container Ω of minimal volume κ = l · w · h, taking into account the placement constraints (1)-(2).

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.

Modeling Non-Overlapping Constraints
Consider a pair of the composed objects T q (u q ) = n q ∪ i=1 T q i (u q ) and T g (u g ) = n g ∪ j=1 T g j (u g ). To describe non-overlapping constraints (1) a phi-function for two composed objects T q (u q ) and T g (u g ) is introduced. It can be written in the form x ≤ 0} be a half space. Denote the half-space P translated along OX on value µ qg ij and rotated by angles Mathematics 2020, 8, 1130 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 (u q , τ qg ij ) and Φ * g and having its radius r q i . The phi-function for the sphere T q i (u q ) and a half space P qg ij has the form while the phi-function for the sphere T g j (u g ) and the half space P * qg ij can be defined as follows and has its radius r q i2 (the cone corresponds to r q i2 = 0). The phi-function for the object T q i (u q ) ∈ and a half space P qg ij has the form 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 i = ( n x i , n y i , n z i ) stands for a unit vector of the external normal to the object T q i (u q ) (see Figure 3). The phi-function for the object T g j (u g ) and a half space P * qg ij has the form  and a half space * qg ij P  has the form

Modeling the Containment Constraints
Let us express the containment constraint ( ) The phi-function for the objects and * Ω can be defined in the form: is the phi-function for the basic convex object ( ) q q i T u and the object * Ω .
The inequality

Containment of a sphere
and the object * Ω can be defined by the equation The inequality Φ qg (u q , u g , τ qg ) ≥ 0 assures the non-overlapping condition int T q (u q ) ∩ int T g (u g ) = ∅ in (1).

Modeling the Containment Constraints
Let us express the containment constraint T q (u q ) ⊂ Ω in the equivalent form as intT q (u q ) ∩ Ω * = ∅, where Ω * = R 3 \intΩ.
The phi-function for the objects T q (u q ) = n q ∪ i=1 T q i (u q ) and Ω * can be defined in the form: where Φ q i (u q ) is the phi-function for the basic convex object T q i (u q ) and the object Ω * . The inequality Φ q (u q ) ≥ 0 implies fulfilling the condition intT q (u q ) ∩ Ω * = ∅ that describes the containment condition (2).
Containment of a sphere T q i (u q ) in a cuboid Ω. The phi-function for the sphere T q i (u q ) and the object Ω * can be defined by the equation Containment of the object T q i (u q ) from the family in a cuboid Ω. The phi-function for the objects T q i (u q ) and Ω * has the form Φ where

Mathematical Model
Now the irregular packing problem can be formulated as the following nonlinear optimization problem: where τ qg = (τ qg ij , i = 1, . . . , n q , j = 1, . . . , n q ) is a vector of the auxiliary variables for the quasi-phi-function Φ qg (u q , u g , τ qg ) (3)-(5) for the objects T q (u q ) and T g (u g ), Φ q (u q ) is the phi-function (6)- (7) for the objects T q (u q ) and Ω * , 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 + 2N + m), where m = q>g∈I N n q n g . The model 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.

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).

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 S q , q ∈ J N around the objects T q (u q ), q ∈ J N .
Step 2. Construct a container Ω 0 with sufficiently large starting length l 0 , width w 0 and height h 0 allowing placement of all the spheres S q , q ∈ J N .
Step 3. Generate a set of N randomly chosen centre points (x 0 q , y 0 q , z 0 q ) for the spheres S q , q ∈ J N inside the container Ω 0 .
Step 4. Grow up the spheres S q of the radii λr q , q ∈ J N , starting from λ = 0 to the full size (λ = 1). Here the decision variables are the centres of S q and the homothetic coefficient λ, where 0 ≤ λ ≤ 1.

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.

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. The corresponding packing is shown in Figure 4a.   The corresponding packings are shown in Figure 5.  The best local minimum obtained for 3478.23 sec. is κ * = l * · w * · h * = 14.889393 · 24.925430 · 17.165791 = 6370.6459961746.
The corresponding packings are shown in Figure 5.
(2) for N = 3 the best local minimum κ * = l * · w * · h * = 12.682284 · 11.050980 · 6.000000 = 840.910031 was obtained for 26.146 sec. starting from the feasible solution  The corresponding packings are shown in Figure 5. The corresponding packings are shown in Figure 6. The corresponding packings are shown in Figure 7. The corresponding packings are shown in Figure 8. The corresponding packings are shown in Figure 8.  The corresponding packings are shown in Figure 9.  The corresponding packings are shown in Figure 9.  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 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. As can be seen from Table 1, different starting points lead to different local minima.

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).