Abstract
The paper considers a packing optimization problem of different spheres and cuboids into a cuboid of the minimal height. Translations and continuous rotations of cuboids are allowed. In the paper, we offer a way of construction of special functions (-functions) describing how rotations can be dealt with. These functions permit us to construct the mathematical model of the problem as a classical mathematical programming problem. Basic characteristics of the mathematical model are investigated. When solving the problem, the characteristics allow us to apply a number of original and state-of-the-art efficient methods of local and global optimization. Numerical examples of packing from 20 to 300 geometric objects are given.
1. Introduction
At present, scientific research concerning mathematical modeling of optimization packing problems of 3D geometrical objects is intensively performed. The great interest to the given class of problems is motivated by the need of wide use of the problems both in scientific researches and applications in different branches of industry. Therefore, when tackling these classes of problems, development of fundamental bases and tools for mathematical and computer modeling is very important.
A state-of-the-art review of bin packing techniques is considered in [1]. It should be noted that there are many publications devoted to packing of cuboids which can be rotated only through 90° around of all coordinate axes.
Currently, 3D packing problems of cuboids, for which translations and continuous rotations are allowed, and spheres are poorly investigated.
The problem of packing nonoriented polyhedrons can be applied in CAD system for rapid prototyping which uses selective laser sintering process of a special powder [2]. Besides, the problem is applied in nanotechnologies for 3D modeling, visual and quality analysis of structural characteristics and mechanical properties of various composite, firm, liquid, glassy materials, granulated media, and biological systems [3–6]. Packing of spheres and cuboids is used to model heterogeneous and porous material morphologies such as concrete, sand, coal, porous explosives, and solid rocket propellants [7].
Also, the problem of packing nonoriented cuboids is applied to car design [8].
Various optimization search algorithms for solving 3D layout problems are considered in paper [9]. In the paper, the authors notice that simulated annealing algorithms and genetic algorithms are stochastic methods that are used in a wide variety of 3D problems.
A lot of authors traditionally use either spheres or cuboids or regular oriented polyhedrons to solve 3D packing problems. For solving problems, the Minkowski sum [10] is used. The sum can be successfully exploited for lattice packing of geometric objects [4, 11].
Torquato and Jiao in paper [4] try to find dense packing of tetrahedra. The authors formulate the problem of generating dense packing of polyhedra within an adaptive fundamental cell subject to periodic boundary conditions as an optimization problem which they call the adaptive shrinking cell scheme.
Paper [12] presents a heuristic approach based on mixed integer programming for solving nonstandard 3D problems with additional conditions. The described approach is addressed to standard MINLP solvers exploiting linear substructures of the mathematical model.
The packing of 3D tetris-like items with orthogonal rotations, within a convex domain (polyhedron) with additional conditions (separation planes, fixed position or orientation for some items, presence of “holes” within the domain, and balancing conditions), is considered in [13]. The author proposes approach which is based on mixed integer linear/nonlinear programming, from a global optimization point of view with an approximate starting solution.
Article [14] is devoted to the review of modern methods of modeling of inclusion, intersection, and contact relations of geometric objects. In the given work, it is noticed that at present one of the perspective approaches for the construction of adequate mathematical models of packing problems is the method of -functions. In work [15], it is shown that the method of -functions allows us to improve results of solution of packing problems of geometrical objects due to application of state-of-the-art nonlinear optimization methods.
Work [16] is devoted to the construction of -functions for oriented (rotations are not allowed) 3D objects whose frontiers are sphere, cuboid, cone, and cylinder. Article [17] considers the construction of the -function for two convex oriented polytopes. In paper [18], authors develop -functions that handle any 2D objects whose boundary is formed by linear segments and/or circular arcs. Constructed -functions support free translations and rotations of objects.
In this work, we consider packing problem of spheres and cuboids. Our mathematical model and approach support free translations and rotations of cuboids. In this respect, it differs from many other approaches that optimize the position of one object at a time. In addition, the solution of the packing problem reduces to minimization of an objective function on a multidimensional space, which can be done by mathematical programming, while most researchers use heuristics for solving packing problems.
The considered packing problem has the following statement. Let there be different cuboids , , different spheres , , and a container (cuboid) , where and are variables. We suppose that cuboids can be both translated by vectors and rotated through angles , ; that is, the location of is defined by a vector . Spheres can be translated by vectors , . Thus, is a motion vector of geometric object , where if and if , . Hence, a vector where determines the location of cuboids , and spheres , , in the Euclidean 3D space .
In what follows, moved by vector is denoted by designated as and of height is denoted by .
Basic Problem. Define a vector which ensures an arrangement of geometric objects , , within , without overlapping, so that the height of will attain the minimal value.
The remainder of the paper is organized as follows. In Section 2, we describe a construction of -functions that are used for construction of a mathematical model of the problem. In Section 3, the mathematical model of the problem is constructed and basic characteristics of the model are investigated. Section 4 defines how starting points are derived and Section 5 presents how a local extremum points are calculated. A smooth jump from one local maximum to one another is considered in Section 6. The construction of promising points is discussed in Section 7. In Section 8, we describe the solution algorithm. Section 9 reports the computational results. Some conclusions are drawn in Section 10.
2. -Functions
The method of -functions permits formulating adequate mathematical models of packing problems at which continuous rotations are allowed. The -functions allow describing the conditions of objects nonoverlapping as a collection of inequality systems whose left-hand sides are infinitely differentiable functions. This permits applying state-of-the-art effective methods of local and global optimization for solving the problem under consideration.
Definition 1. An everywhere defined and continuous function is called a -function for objects and if the function possesses the following characteristic properties: where is the interior of and is the frontier of .
Thus, for any location of two objects and in the Euclidean 3D space , the corresponding -function shows how far these objects are from each other, whether they touch each other, or whether they overlap (in the latter case, the one shows how large the overlap is).
Consequently, in order to formulate a mathematical model of the problem stated, it needs to construct -functions for pairs of cuboids, cuboid and sphere, and pairs of spheres.
2.1. -Function for Two Cuboids and
The article [19], describes a construction of -function for pairs of convex polytopes. As an example, the -function for two cuboids is given in the article.
2.2. -Function for Cuboid and Sphere
Let there be a cuboid and sphere . Vertex coordinates of are indexing as follows: , , , , , , , and (Figure 1).
The -function for a cuboid and a sphere is completely defined by an interaction of the following subsurfaces of cuboid and sphere:(1)a spherical surface of and facets of ;(2)a spherical surface of and edges of ;(3)a spherical surface of and vertices of .
If , then 0-level surface of has the form depicted on Figure 2(a). It is easily seen that the 0-level surface consists of subsurfaces of planes (Figure 2(b)), cylinders (Figure 2(c)), and spheres (Figure 2(d)).
(a)
(b)
(c)
(d)
The first tangency (Figure 2(b)) form generates functions , , which are specified as follows: where
It is easily verified that if at least one inequality of kind , then . Based on the reasoning and relations (2), we construct the function
Thus, if holds true, then . On the other hand, if , then there can be both and (see Figures 2(c) or 2(d)).
The second tangency form derives cylindrical parts of the 0-level surface of (see Figure 2(c)). The parts belong to cylinders whose axes pass through edges of and are specified by the following equations: where , , and , .
In order to cut necessary parts of the cylinders, we form the following equations of planes: which are parallel to edges of and pass through pairs of points lying at distance from vertices on the extending of edges of (see Figure 3).
Making use of the functions (5) and (6), , we form the following functions:
It is evident that if at least one of the inequalities , , holds true, then . The property allows us to construct the function
Thus, if , then . On the other hand, if , then there can be both or (see Figure 2(d)).
The third tangency form (see Figure 2(d)) generates spherical parts of the 0-level surface of . The parts belong to spheres that are formulated by equations where , ; that is, .
Note that centers of the spheres coincide with vertices , , of cuboid correspondently. We need to extract the necessary spherical sectors of the spheres (Figure 2(d)).
To this end, we construct for each of vertices , , of three elliptic cylinders oriented by a special way with respect to and a helper plane. Let us consider the case in detail.
Let equations of elliptic cylinders , , (Figure 4) have the form
We rotate the cylinders , , and through angles and round axes , , and , respectively. As a result, obtained cylinders , , (see Figure 5) are specified by the equations, respectively,
It is evident that the intersections of the coordinate plane and cylinders , , the coordinate plane and cylinders , , and the coordinate plane and cylinders , are circles of radius . We rotate cylinders , , through the angles , , , and as follows: cylinders and : we rotate round axes and, as a result, cylinders , , are obtained (see Figure 6); cylinders and : we rotate round axes and, as a result, cylinders , , are obtained; cylinders and : we rotate round axes and, as a result, cylinders , , are obtained.
Then, we translate the triple cylinders , , and , , by the vectors , . Thus, we obtain at each of the vertices , , of three elliptic cylinders , , and (Figure 7) which are described by the following equations: where , .
In order to cut necessary parts of the elliptic cylinders at each of the vertices of , we form the following equations of helper planes passing through triples of points lying at distance from vertices on the extending of edges of (Figure 8): where , .
Making use of the functions (9)–(13), we form the functions
It is easily verified that if at least one of the inequalities , , holds true, then . Based on (14), we construct the function
Thus, if , then .
Note that if at least one of inequalities either (4) or (8) or (15) is satisfied, then . It allows us to construct the -function for cuboid and sphere as follows:
2.3. -Function for and Object
It follows from paper [19] that the -function for and has the form where
2.4. -Function for and Object
On the ground of (5), the -function for and can be formulated as where
3. A Mathematical Model of the Basic Problem and Its Characteristics
In this Section, the mathematical model is constructed as a typical problem of nonsmooth mathematical programming. Characteristics of the mathematical model permit presenting the feasible region as a union of subregions which are specified by the inequality systems whose left-hand sides are infinitely differentiable functions. This allows us to apply the state-of-the-art nonlinear optimization methods to solve the problem.
Making use of -functions described in Section 2, a mathematical model of the stated problem can be formulated as where , , and ensure nonoverlapping and , and and and , respectively; and guarantee an arrangement of and within , respectively.
It follows from the form of -functions that the mathematical model possesses the following characteristics.(1)The feasible region is specified by an inequality system which includes operators “” and “” [20]. The operators contain logical constituents such as set-theoretical operations of conjunction and disjunction, respectively. Thus, in inequalities and , the operators “” generate a “union” or collection of inequality systems whereas the operators “” yield inequality systems. This means that the feasible region can be always represented by a finite union of subregions; that is, , where , 156, and 26 are quantity of inequality systems in and , respectively, , .(2)The feasible subregions are described by inequality systems where , , and , , are inequality systems of inequality collections which form inequalities and . Note that functions on the left-hand sides of the inequalities in (23) are nonlinear and infinitely differentiable. Thus, a point if at least one of inequality systems (23) is satisfied.(3)Since some subregions , , can be empty or can be contained in other subregions, then the number of the systems is much less than .(4)Some subregions of kind may have common points. This means that if a point is a local minimum with respect to , then one has to investigate all other subregions containing to prove that the point is a local minimum with respect to .(5)The matrix of inequality system (23) is strongly sparse.(6)The problem (21)-(22) is NP-hard [21].
The characteristics involve that a solution algorithm must include the following steps: a construction of starting points, a calculation of local minima, and a nonexhaustive search of local minima which ensures a good approximation to a global minimum.
4. Construction of Starting Points
In order to solve the problem (21)-(22), starting points belonging to the feasible region have to be constructed. When tackling cutting and packing problems, greedy and heuristic algorithms are usually exploited to derive starting points [22]. The algorithms do not permit us to obtain any starting points. This means that the algorithms significantly contract a set of local minima to be considered. Since cuboids are free rotated, then rules of their placement are very complex. This means that time expenditures when constructing starting points using greedy and heuristic algorithms essentially increase in 3D space [23]. The circumstances result in a need of development of new approaches to generate starting points.
For construction of starting points, we suggest a method which assumes that all metric characteristics of objects are variable and take values between the minimal and initial values. To find the starting point, a random generation of object placement parameters is executed. In addition, the sizes of the objects essentially diminish. After this, we solve the problem of nonlinear programming providing an increase of object sizes. If, as a result, the point of local maximum is obtained at which all sizes of objects reach the initial values, then the motion vector at the local maximum point and the given height are taken as a starting point for the basic problem.
Let us consider one of such approaches.
Let sizes of cuboids , , and radii of spheres , , be variables. Thus, -function for objects and depends on metric characteristics of ones as well; that is, the -function has the kind .
The variables , , form a vector , where . Thus, a vector of all variables is , where , .
We give and a point by a random way so that , , (Figure 9(a)). The coefficient 0.1 insures a correct construction of -functions.
(a)
(b)
Then, taking the starting point , we solve the problem where
Note, that as a result of local optimization of the problem, the metric characteristics increase. In addition, they are limited by their initial values due to the inequalities , , , .
As a solution result of the problem, a local maximum point is defined (Figure 9(b)).
The problem (24)-(25) possesses the characteristics that are similar to the ones of the problem (21)-(22). Furthermore, the problem (24)-(25) is master of additional properties.(1)If is a local maximum point of the problem (24)-(25) and , then ; that is, cuboids , , and spheres , , are packed into a cuboid . This means that, in this case, the point is a global maximum point of the problem (24)-(25).(2)If is such that at least one of components of is strictly less matched component of , then and the point is either a local maximum point or a global maximum point of the problem (24)-(25). In addition, if is a global maximum point, then cuboids , , and spheres , , cannot be packed into a cuboid .
Note that we always can choose so that an arrangement of geometric objects , , within is guaranteed.
5. Calculation of Local Minima
Since we deal with a classical problem of mathematical programming, then the nonlinear optimizations methods can be used to search for local minima. The problem characteristics allow reducing solving the problem to solving a sequence of the subproblems whose feasible regions are specified by inequality systems. In addition, the left-hand sides of inequalities are infinity differentiable functions; that is, for solving the subproblems, the optimization gradient methods can be used.
If (24), then (22). So, the point is taken as a starting point and the problem (21)-(22) is tackled. As a result, a local minimum point is computed (Figure 10).
Let be a local minimum point of the problem (21)-(22) obtained after iterations. We take
Assuming in (26), we derive a point . It is evident because of inequalities of kind , , , , which are broken at the point (Figure 11).
(a)
(b)
(c)
In order to obtain a point , primarily, we take a starting point , where , , and solve the helper problem where
It is easily seen that .
Let be a solution of the problem (27)-(28), where because of the inequality ; that is, is a global minimum point.
Then, taking the starting point , we tackle the problem (24)-(25) and obtain a local maximum point .
Two cases are possible: (Figure 11(c)) and (Figure 12).
Let ; that is, is a global maximum of the problem (24)-(25). Hence, . The point is not in the general case a local minimum point of the problem (21)-(22). So, taking the starting point , we solve the problem (21)-(22). As a result, a new local minimum point is computed. Then, we assume in (26), take a starting point (28), and sequentially solve the problems (27)-(28), (24)-(25), (21)-(22), and so on until is fulfilled.
If , we try to transit from the local maximum point to local maximum point so that (see Section 6). If we cannot execute such transition, then giving in (26) we take a starting point and sequentially solve the problems (27)-(28), (24)-(25), (21)-(22), and so on until , where is a solution accuracy, is attained.
To calculate local extrema of the problems (21)-(22), (24)-(25), and (27)-(28), the modification of the Zoutendijk method of feasible directions [24] along with the concept of -active inequalities [23, 25] is used.
A computation of local extrema of the problems (21)-(22), (24)-(25), and (27)-(28) is reduced to solving a sequence of optimization subproblems on the feasible subregions , , , , and , , respectively [26].
6. Transition from One Local Maximum to Another
Since exhaustive search of all local maxima of problem (24)-(25) is theoretically possible only, then heuristic algorithms, probabilistic approaches, or their combinations are exploited to find approximations to a global maximum.
We offer a special approach which provides a jump from a local maximum point to another one so that the objective function value may enlarge.
To this aim, we use a helper problem which maximizes sizes of geometric objects at a local maximum point of the problem (24)-(25) without restrictions on their values. It gives us a capability to derive a vector of the steepest ascent of the function (24). Since restrictions of the metric characteristics of objects are omitted in the helper problem, then some metric characteristics increase and other ones decrease at the point of local maximum point of the problem (24)-(25). In other words, the vector shows whether there exists “free space” which can be exploited to improve a value of the objective function in the problem (24)-(25) and consequently to increase sizes of objects if any.
Let be a local maximum point of the problem (24)-(25), , and . Since , then the point (22). In order to transit from the point to the point so that , we consider the helper problem where We note that the problem differs from the problem (21)-(22) by absence of restrictions This means that sizes of geometric objects can take any nonnegative values.
We compute the steepest ascent vector at the point for the problem (29)-(30) and construct the point
The positive elements of the vector corresponding to the direction of objects sizes change allow us to define objects whose sizes can be increased in the point of local maximum. Consequently, there is a “free space” in the neighborhoods of such objects which can be used to increase the local maximum value. Let us suppose that some coordinates of satisfy the inequalities , , , , , , and , . Let sets consist of , , elements and let consist of elements. Then, consists of elements and .
Remark 2. If , then there exists such that if , then .
We compute
Let , and , .
On the ground of points and (see Remark 2), we derive coordinates of the vector as follows.
(1) Let ( and are cuboids and ), or or and ; that is, at least one of metric characteristics of does not reach the initial value and one of metric characteristics of is strictly greater than its value at the point . For example, let the initial sizes of cuboids and be , , and and , , and and let sizes at the point be , , and and , , and , respectively. To the points and , there corresponds packing shown in Figure 13(a) and Figure 13(b), respectively.
(a)
(b)
(c)
(d)
Then if inequalities hold true then coordinates , , and are constructed as where is either or or (33). This means that cuboids and are “changed over” relatively their arrangement corresponding to the point (Figure 13(c)).
If for at least one of inequalities (34) is broken, then
(2) Let ( and are spheres and ), and , that is radius of does not reach the initial value and radius of is strictly greater than its value at the point . If inequalities are satisfied then coordinates , , and take the form where and .
Thus, spheres and are “changed over”.
If for at least one of the inequalities (37) is broken, then
(3) Let , ( and are cuboid and sphere ), , and ; that is, the radius of circle circumscribed round does not reach the initial value and radius of is strictly greater than its value at the point . If inequalities are satisfied, then Thus, cuboid and sphere are “changed over” relatively their placement corresponding to the point .
If for at least one of the inequalities (40) is broken, then
(4) Let , ( and are sphere and cuboid ), , and ; that is, radius of does not reach the initial value and the radius of circle circumscribed round is strictly greater than its value at the point . If inequalities are satisfied, then Thus, sphere and cuboid are “changed over” relatively their placement corresponding to the point .
If for at least one of the inequalities (43) is broken, then we take
Remark 3. It is easily seen that there can be ; that is, there exists such integer that, if , then .
Thus, if , then and are in attraction zones of different local maximum points; that is, to the point corresponds the packing of , , obtained from the packing of , , corresponding to the point in which some objects of , , “change over.”
Let be a local maximum point obtained from the starting point .
Theorem 4. If , then .
Proof. For the sake of simplicity, let and differ in values of variables , , , and only and let , and ; that is, , where , , , either or , either or , and either or .
Note that because of . Hence, the point is not in attraction zone of the local maximum point . It follows from Remark 2 that ; that is, the inequality system , at the point is satisfied. It follows from the relations (35) that inequalities , , , , and , , at the point , are satisfied as well. This means that . Bearing in mind the inequalities (34), we have ; that is, . It is evident that if is a local maximum point obtained from the starting point (see Figure 13(d)). If the number of coordinates of vectors and differs in the larger number of coordinates, then the theorem holds true all the more. For the cases 2, 3, and 4, the theorem is proved in the same way.
7. Construction of Promising Points
The theorem conditions above always ensure a jump from a local maximum point to another one at which the objective function has a more value. However, severe theorem constraints do not always allow us to fulfill such jump although the jump may execute. So, we extend a set of the starting points using points and (see Section 6). Since, in general case, local maximum points obtained from the starting points do not guarantee an increase of values of , then in what follows the points are called promising starting ones.
Let us consider a way of construction of promising starting points on the ground of the theorem and points (a local maximum of the problem (24)-(25)), (see (32) and Remark 2).
Firstly, making use of and (see the problem (29)-(30)), we derive coordinates of points , and as follows: