COMPARING A STATIC EQUILIBRIUM BASED METHOD WITH THE SUPPORT FACTOR FOR HORIZONTAL CARGO STABILITY IN THE CONTAINER LOADING PROBLEM

ABSTRACT

This paper presents an approach to deal with horizontal cargo stability in container loading problems. Cargo stability has been explored mainly with support factors that constrain the minimum area of each box’s faces to be supported by other boxes. On the other hand, we propose an approach based on the static equilibrium of rigid bodies to check the static stability of a given packing. The approach is used as a cutting plane routine in a branch-and-cut framework to the single container loading problem. This framework considers the resolution of an integer linear programming model to obtain feasible packings next checked with the proposed approach to avoid unstable packings. The computational experiments consider 180 benchmark instances on which stable solutions of the proposed approach have better container fill rates than the support factor approach. In terms of lateral support, the proposed approach provides the minimum value inferior to 70% on average, which is satisfactorily smaller and less restrictive than the full support. Results also indicate that more unstable solutions emerge from refined grids and fewer types of boxes available.

Keywords:
container loading problem; horizontal cargo stability; equilibrium of material bodies

1 INTRODUCTION

This work is related to packing problems in which a set of items (e.g., boxes) must be packed in one or more containers to optimize an objective function (e.g., the container fill rate). We study the Single Container Loading Problem (SCLP). A subset of three-dimensional rectangular boxes is packed into a three-dimensional rectangular container to maximize the occupied container volume. We assume boxes are packed with each of their faces parallel to the container faces (orthogonal packing), and boxes cannot rotate (fixed orientation).

The literature has extensively addressed the SCLP. Exact methods were studied, for example, by ^{Hifi (2004}), who proposed a graph-search approach; ^{Junqueira et al. (2012b)}, who developed an integer programming model solved with a branch-and-cut algorithm; and, ^{Paquay et al. (2016}), who proposed a model that uses two reference positions to pack boxes. ^{Silva et al. (2019}) did a comparative study of the models and exact algorithms proposed by ^{Chen et al. (1995}), ^{Martello et al. (2000}), ^{Fasano (2004}), ^{Hifi et al. (2010)}, ^{Junqueira et al. (2012b)}, ^{Tsai et al. (2015}), and ^{Paquay et al. (2016)}. These algorithms solved all instances with few types of boxes (up to 7). For instances with more types (up to 15), the model of ^{Junqueira et al. (2012b)} outperformed the others. On the other hand, for instances with up to 250 types, the model of ^{Junqueira et al. (2012b)} required much memory, while the model of ^{Chen et al. (1995)}, followed by the algorithm of ^{Martello et al. (2000)}, performed better.

Recently, in ^{Nascimento et al. (2021}), a three-step exact algorithm was developed for the SCLP in the presence of twelve practical constraints. The authors combined integer linear and constraint programming models iteratively, with cuts to break symmetric solutions. In the computational experiments, the approach was compared with the model of ^{Junqueira et al. (2012b)} on instances proposed by the latter. The comparisons in the absence of practical constraints showed the model of ^{Junqueira et al. (2012b)} is a better choice for the problem. In this way, ^{Kurpel et al. (2020}) extended the model of ^{Junqueira et al. (2012b)} and developed new bounds to handle multiple container loading problems with practical constraints. The authors obtained new optimal or improved solutions for some problem variants.

Concerning heuristic methods, we may cite those where boxes are packed in layers/blocks with a filling method, as in ^{Bortfeldt & Gehring (2001}) and ^{Sheng et al. (2017}). Some authors have combined metaheuristics with filling methods, as in ^{Jamrus & Chien (2016}), with the use of genetic algorithm; ^{Tao & Wang (2015}), with the use of tabu search; and, ^{Egeblad & Pisinger (2009}) with the use of simulated annealing. Recently, ^{Silva et al. (2020}) proposed a matheuristic that combines a wall-building approach and two integer models to handle the SCLP with seven practical constraints.

One of the best heuristics for the SCLP is the beam search of ^{Araya & Riff (2014}) that packs blocks of boxes with higher utilization by using a tree search method. With this heuristic, the authors could outperform all previous results on benchmark instances in terms of occupied volume and runtime. This beam search was improved by ^{Araya et al. (2017)} with a new evaluation function for ranking and selecting blocks. The improved beam search results are significantly better compared to the previous ones on all sets of instances.

We consider the SCLP with the additional constraint of static vertical and horizontal cargo stability related to the capability of boxes to resist inertia, given external forces acting on them (see, e.g., ^{Queiroz et al. (2019})). The literature often considers a support factor to deal with this constraint, implying a minimum percentage area of each box’s faces to be supported by other boxes or by the container. To guarantee stability, some authors used the full support, as ^{Fanslau & Bortfeldt (2010}), ^{Junqueira et al. (2012a)}, ^{Araya et al. (2017}), ^{Silva et al. (2020}), and ^{Nascimento et al. (2021}).

Other authors instead have considered it as acceptable that only a percentage of the boxes’ faces should be supported to meet cargo stability, as ^{Christensen & Rousoe (2009}) and ^{Mack et al. (2004}). Some authors assumed a packing as vertically stable if adjacent boxes support the four bottom vertices of each box (see, e.g., ^{Egeblad et al. (2010})), or three of these vertices (see, e.g., ^{Paquay et al. (2016})). On the other hand, more realistic approaches to handle cargo stability are based on the static equilibrium of rigid bodies (^{Oliveira et al., 2021}).

Few authors used the mechanical equilibrium principle to handle cargo stability in packing problems. In this case, boxes are satisfying the static equilibrium equations for material bodies (^{Hibbeler, 2010}). In this line, the literature has focused on vertical cargo stability for the SCLP, assuming only the force that emerges from the gravitational field, with the contributions of ^{Silva et al. (2003}), ^{Ramos et al. (2016}), and more recently ^{Oliveira et al. (2021}).

Concerning horizontal cargo stability with the mechanical equilibrium principle, the number of contributions is indeed small. ^{Bracht et al. (2016}) used a physical simulation engine to heuristically check if the box’s center of gravity moves from its original position. The authors assumed two situations in which the gravitational, centripetal, and frictional forces acted on the packing. In the first situation, the container’s lateral faces were disregarded, while the lateral faces were taken into consideration in the second one. ^{Ramos et al. (2015}) commented that the literature has been considering two metrics. The first one is based on the average number of boxes supporting other boxes, excluding those placed directly on the floor. The second metric is based on the support of the boxes’ faces. They developed two new indicators for evaluating the horizontal stability of a cargo. These indicators are the number of fallen boxes and the number of boxes inside the damaged area of a damage curve limit.

^{Queiroz et al. (2019}) extended the mechanical equilibrium-based method in ^{Queiroz & Miyazawa (2014)} to horizontal stability, limited to two-dimensional packing problems. Recently, ^{Oliveira et al. (2021}) proposed two similar approaches for vertical stability. The first is based on the dynamic simulation of rigid bodies. The second is based on the equilibrium of buildings. The authors used the model of ^{Junqueira et al. (2012b)} solved with a branch-and-cut algorithm, where cutting planes to avoid unstable vertical packings were added on demand. The main difference between our work and the previous ones relies on the proposal of a mechanical equilibrium-based method for horizontal cargo stability in three-dimensional packing problems.

Differently from ^{Silva et al. (2003}), ^{Queiroz & Miyazawa (2014}), ^{Ramos et al. (2016}), and ^{Oliveira et al. (2021}) who only handled vertical cargo stability, we consider lateral forces, as centripetal and frictional ones, acting on the cargo. This is consistent with ^{Queiroz et al. (2019)}, but their method is limited to two-dimensional problems. Similarly to ^{Oliveira et al. (2021)}, we extend the method of ^{Baraff (1989}), which relies on the dynamic simulation of rigid bodies and investigates the impact of different grids of points in the presence of horizontal stability. It is essential to mention that we first check static vertical stability with the method of ^{Oliveira et al. (2021)} and next static horizontal stability for forces acting horizontally from left to right (along the x-axis) and from front to back (along the y-axis).

Following the literature (^{Junqueira et al., 2012b}; ^{Ramos et al., 2016}; ^{Oliveira et al., 2021}), one way to evaluate the proposed method is to compare with the support factor in terms of occupied volume. For that, we consider an existing solution method (i.e., the model of ^{Junqueira et al. (2012b)} solved with a branch-and-cut algorithm) and instances thought in the presence of cargo stability constraints. In the model of ^{Junqueira et al. (2012b)}, boxes are packed over a grid, which allows us to study the impact of the grid on the final solution.

This paper is organized as follows. In Section 2, we revisit the integer model of ^{Junqueira et al. (2012b)} and the support factor constraints to model cargo stability. Our solution method for the SCLP consists of solving this model with a branch-and-cut algorithm. For every feasible packing found, we first check it for vertical stability and next for horizontal stability. If the packing is unstable for vertical or horizontal stability, we add a valid inequality to avoid one of the boxes of being packed in its current position. The proposed method for horizontal stability is discussed in Section 3, taking advantage of the method for vertical stability. Section 4 follows the literature and presents computational experiments on 180 benchmark instances, where the proposed method is compared with the support factor on three different grids. Finally, Section 5 has some conclusions and suggestions for future works.

2 INTEGER FORMULATION

The SCLP considers a container with dimensions (L, W, H), where the length L, width W, and height H, are related to the x-, y-, and z-axes, respectively. We refer to the bottom-frontal-left vertex as the container’s origin point (0,0,0) in the Cartesian space. Each box of type $k=1,\dots ,n$, has dimensions (l _k , w _k , h _k ), volume v _k , and a maximum number of copies b _k . The dimensions of containers and boxes are assumed to be integers, and boxes are packed orthogonally, with fixed orientation (i.e., rotations are not allowed). In a feasible packing, boxes do not overlap each other and must respect the container dimensions. The objective is to obtain a feasible packing of maximum volume. The sets of available coordinates are X _c , Y _c , and Z _c , each one containing natural numbers from 0 to the respective container dimension along the x-, y-, and z-axes. From these sets, we define the feasible grid to pack a box of type k as $X_{ck}\times Y_{ck}\times Z_{ck}$:

$X_{ck}=\{0,1,\dots ,L-l_k\},$ (1)

$Y_{ck}=\{0,1,\dots ,W-w_k\},$ (2)

$Z_{ck}=\{0,1,\dots ,H-h_k\}.$ (3)

Let binary variables $x_{kpqr}\in \{0,1\}$ represent the decision of packing a box of type k with its bottom-frontal-left corner at the point $(p,q,r)\in X_{ck}\times Y_{ck}\times Z_{ck}$ The model (4)-(7) solves the SCLP (^{Junqueira et al., 2012b}):

$\text{max}\sum _{k=1}^n\sum _{p\in X_{ck}}\sum _{q\in Y_{ck}}\sum _{r\in Z_{ck}}{v}_k{x}_{kpqr}$ (4)

$$\begin{gathered} \text{s.t:}\sum _{k=1}^n\sum _{\{p\in X_{ck}|p\text{'}-l_k+1\le p\le p\text{'}\}}\sum _{\{q\in Y_{ck}|q\text{'}-w_k+1\le q\le q\text{'}\}}\sum _{\{r\in Z_{ck}|r\text{'}-h_k+1\le r\le r\text{'}\}}{x}_{kpqr}\le 1, \\ \forall p\text{'}\in X_c,q\text{'}\in Y_c,r\text{'}\in Z_c, \end{gathered}$$ (5)

$\sum _{p\in X_{ck}}\sum _{q\in Y_{ck}}\sum _{r\in Z_{ck}}{x}_{kpqr}\le b_k,k=1,\dots ,n,$ (6)

$x_{kpqr}\in \{0,1\},k=1,\dots ,n,p\in X_{ck},q\in Y_{ck},r\in Z_{ck}.$ (7)

The objective function (4) maximizes the container fill rate. Constraints (5) guarantee the nonoverlapping of boxes, while constraints (6) impose the limit on the number of boxes of each type in the solution. Constraints (7) define the variables domain.

2.1 Cargo Stability with Support Factor

The linear constraints to model cargo stability with the support factor are given next (^{Junqueira & Queiroz, 2022}). Such constraints are added to the integer model (4)-(7), resulting in the JMY approach. The support factor in vertical stability refers to the support of the bottom face of each box. In horizontal stability, the boxes’ lateral faces are supported by the lateral faces of other boxes.

For vertical stability, let $\alpha \in [0,1]$ indicate the percentage of support for the boxes’ bottom faces. Constraints 8 model the support of the boxes’ bottom faces:

$$\begin{gathered} \sum _{\Phi }\sum _{\Psi }\sum _{\Omega }{\overline{L}}_{kj}{\overline{W}}_{kj}{x}_{kpq(r\text{'}-h_k)}\ge \alpha l_j{w}_j{x}_{jp\text{'}q\text{'}r\text{'}}, \\ j=1,\dots ,n,p\text{'}\in X_{cj},q\text{'}\in Y_{cj},r\text{'}\in Z_{cj}\backslash \left\{0\right\}, \end{gathered}$$ (8)

with Φ as $\{k=1,\dots ,n|r\text{'}-h_k\ge 0\}$, Ψ as $\{p\in X_{ck}|p\text{'}-l_k+1\le p\le p\text{'}+l_j-1\}$, Ω as $\{q\in Y_{ck}|q\text{'}-w_k+1\le q\le q\text{'}+w_j-1\}$, ${\overline{L}}_{kj}=\min (p+l_k,p\text{'}+l_j)-\max (p,p\text{'})$, and ${\overline{W}}_{kj}=\min (q+w_k,q\text{'}+w_j)-\max (q,q\text{'})$.

For horizontal stability, let $\beta \in [0,1]$ and $\gamma \in [0,1]$ represent, respectively, the percentage of support of the left and frontal faces of each box. Constraints (9) and (10) model the horizontal stability for, respectively, the x and y-axes.

$$\begin{gathered} \sum _{\epsilon }\sum _{\kappa }\sum _{\xi }{\overline{W}}_{kj}{\overline{H}}_{kj}{x}_{k(p\text{'}-l_k)qr}\ge \beta w_j{h}_j{x}_{jp\text{'}q\text{'}r\text{'}}, \\ j=1,\dots ,n,\forall p\text{'}\in X_{cj}\backslash \left\{0\right\},\forall q\text{'}\in Y_{cj},\forall r\text{'}\in Z_{cj}, \end{gathered}$$ (9)

with $\epsilon =\{k=1,\dots ,n|p\text{'}-l_k\ge 0\}$, $\kappa =\{q\in Y_{ck}|q\text{'}-w_k+1\le q\le q\text{'}+w_j-1\}$, $\xi =\{r\in Z_{ck}|r\text{'}-h_k+1\le r\le r\text{'}+h_j-1\}$, ${\overline{W}}_{kj}=\min (q+w_k,q\text{'}+w_j)-\max (q,q\text{'})$, and ${\overline{H}}_{kj}=\min (r+h_k,r\text{'}+h_j)-\max (r,r\text{'})$.

$$\begin{gathered} \sum _{\rho }\sum _{\varrho }\sum _{\sigma }{\overline{L}}_{kj}{\overline{H}}_{kj}{x}_{kp(q\text{'}-w_k)r}\ge \gamma l_j{h}_j{x}_{jp\text{'}q\text{'}r\text{'}}, \\ j=1,\dots ,n,\forall p\text{'}\in X_{cj},\forall q\text{'}\in Y_{cj}\backslash \left\{0\right\},\forall r\text{'}\in Z_{cj}, \end{gathered}$$ (10)

$\rho =\{k=1,\dots ,n|q\text{'}-w_k\ge 0\}$, $\varrho =\{p\in X_{ck}|p\text{'}-l_k+1\le p\le p\text{'}+l_j-1\}$, $\sigma =\{r\in Z_{ck}|r\text{'}-h_k+1\le r\le r\text{'}+h_j-1\}$, ${\overline{L}}_{kj}=\min (p+l_k,p\text{'}+l_j)-\max (p,p\text{'})$, and ${\overline{H}}_{kj}=\min (r+h_k,r\text{'}+h_j)-\max (r,r\text{'})$.

2.2 Grids of Points

^{Herz (1972}) proposed an exact approach to generate two-dimensional cutting patterns based on a tree search. To limit the number of nodes explored, the author introduced the canonical dissections (CD). Following the CD, the positions along the dimensions of the container, where each box can have its bottom-frontal-left vertex, are given by the sets (11)-(13). These positions are obtained from non-negative integer linear combinations of boxes dimensions. The sets X _c , Y _c , and Z _c can be replaced by, respectively:

$X_d=\{p\in X_c|p=\sum _{k=1}^n{\delta }_k{l}_k,p\le L,{\delta }_k\le b_k,{\delta }_k\in {\mathrm{\mathbb{Z}}}_{+}\},$ (11)

$Y_d=\{q\in Y_c|q=\sum _{k=1}^n{\delta }_k{w}_k,q\le W,{\delta }_k\le b_k,{\delta }_k\in {\mathrm{\mathbb{Z}}}_{+}\},$ (12)

$Z_d=\{r\in Z_c|r=\sum _{k=1}^n{\delta }_k{h}_k,r\le H,{\delta }_k\le b_k,{\delta }_k\in {\mathrm{\mathbb{Z}}}_{+}\}.$ (13)

To obtain a smaller number of points, ^{Scheithauer & Terno (1996}) proposed the reduced raster points (RRP), which are represented by sets (14)-(16):

$X_r=\left\{\right(L-p{)}_x|p\in X_d\},$ (14)

$Y_r=\left\{\right(W-q{)}_y|q\in Y_d\},$ (15)

$Z_r=\left\{\right(H-r{)}_z|r\in Z_d\},$ (16)

with $(L-p{)}_x=\max \{s\in X_d|s\le L-p\}$, $(W-q{)}_y=\max \{t\in Y_d|t\le W-q\}$, and $(H-r{)}_z=\max \{u\in Z_d|u\le H-r\}$.

The CD and RRP can be used to reduce the number of variables of model (4)-(7). Both CD and RRP grids preserve an optimal solution for orthogonal packing problems (^{Almeida Cunha et al., 2020}). However, when cargo stability is enforced, grids of points may lead to the loss of stable packings.

3 CARGO STABILITY WITH THE STATIC EQUILIBRIUM OF RIGID BODIES

Based on Newton’s laws of motion, we consider a set of equations that a rigid material body must attend to stay in static equilibrium: the sum of all forces must be zero, and the sum of all torques must be zero. Details about the equilibrium of material bodies can be found in ^{Hibbeler (2010}). When these equations are satisfied for each box, the packing is assumed to be statically stable (^{Queiroz et al., 2019}). Our modeling assumes that each box is a rigid body, made of the same material, and has homogeneous density.

Considering the packing is being transported by road and horizontally stable, we need to check if it is stable. For that, we first check the packing for static vertical stability with the B approach described in ^{Oliveira et al. (2021}). Next, we extend such an approach to check static horizontal stability. In static stability, we model situations where the cargo, in particular, boxes, cannot translate nor rotate due to the action of external (vertical and horizontal) forces; otherwise, the packing is assumed as unstable (^{Queiroz et al., 2019}).

We consider four different kinds of forces to act on a box, although it is straightforward to assume the action of other forces. The four forces are: weight, due to the gravitational field (with a vertically downward gravitational acceleration of $g=9.8m/s^2$) permanently acting on boxes; normal and friction, due to contact that boxes have each other; and centripetal, since the cargo is following a curved path and then such force keeps the cargo on the path. Horizontal forces (e.g., centripetal) can make boxes slide or generate a moment that makes boxes tipping over (^{Hibbeler, 2010}).

Forces are assumed to be acting on the center of gravity of the boxes unless stated otherwise. As boxes are assumed homogeneous, the center of gravity coincides with the box’s geometric center. It is essential to mention that the normal and friction forces are reaction forces of the weight and centripetal ones, respectively.

Two adjacent boxes are touching themselves on infinite points, as shown in Figure 1. We assume that only the vertices, from now named contact points, of the contact surface are considered. On each contact point i acts a triad of forces ${\overrightarrow{f}}_i=(f_i,f_i^x,f_i^y)$, where f _i is the value (i.e., magnitude) of the normal force (on the direction of the z-axis), and $f_i^x$ and $f_i^y$ are the value of the friction forces on the direction of the x and y-axes, respectively.

Figure 1
Triad of forces acting on a contact point.

The approach to obtain the equilibrium equations for packing problems uses Newton’s laws of motion under the assumption that no box can interpenetrate each other (i.e., the forces must guarantee that different boxes will never share/occupy the same physical space) nor move away (^{Baraff, 1989}). In general, our approach for cargo stability consists of initially calculating the normal forces on each contact surface and verifying if they keep the initially adopted sign (i.e., check the static vertical stability). Next, we check if each box slides or tips over (i.e., check the static horizontal stability). In the following, we briefly describe the B approach for vertical stability (^{Oliveira et al., 2021}) since it is required to deal with horizontal stability.

3.1 Revisiting the B Approach

Each contact point i of a given contact surface $c\in {\mathcal{C}}_k$ (on the direction of the z-axis) of box k has acceleration $\ddot{p_i}$. This acceleration should be non-negative since boxes cannot accelerate so as to interpenetrate each other (^{Baraff, 1989}):

${\ddot{\chi }}_i=\ddot{p_i}\hat{n}=\left(\frac{{\overrightarrow{F}}_k}{m_k}+\frac{{\overrightarrow{\tau }}_k}{I_k}\times {\overrightarrow{r}}_i\right)\hat{n},$ (17)

where $\hat{n}=(0,0,1)$ represents the direction of the normal force, m _k is the mass of box k, I _k is the inertia tensor of box k, and ${\overrightarrow{r}}_i$ is the distance vector between the item’s center of gravity and point i. It follows that the resultant force ${\overrightarrow{F}}_k$ and resultant torque ${\overrightarrow{\tau }}_k$ are, respectively,

${\overrightarrow{F}}_k=m_k\overrightarrow{g}+f_1\hat{n}+f_2\hat{n}+f_3\hat{n}+f_4\hat{n},$ (18)

${\overrightarrow{\tau }}_k={\overrightarrow{r}}_1\times f_1\hat{n}+{\overrightarrow{r}}_2\times f_2\hat{n}+{\overrightarrow{r}}_3\times f_3\hat{n}+{\overrightarrow{r}}_4\times f_4\hat{n}.$ (19)

In general, a system with ε inequalities (i.e., ε contact points) can be written in the matrix form (20):

$A\overrightarrow{f}-\overrightarrow{e}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1\epsilon }\\ a_{21}& a_{22}& \cdots & a_{2\epsilon }\\ ⋮& ⋮& ⋮& ⋮\\ a_{\epsilon 1}& a_{\epsilon 2}& \cdots & a_{\epsilon \epsilon }\end{array}\right]\left[\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_{\epsilon }\end{array}\right]-\left[\begin{array}{c}g\\ g\\ ⋮\\ g\end{array}\right]=\left[\begin{array}{c}{\ddot{\chi }}_1\\ {\ddot{\chi }}_2\\ ⋮\\ {\ddot{\chi }}_{\epsilon }\end{array}\right]\ge \overrightarrow{0}.$ (20)

Matrix A has order $\epsilon \times \epsilon$. Vectors $\overrightarrow{f}$ and $\overrightarrow{e}$ have order $\epsilon \times 1$. Vector $\overrightarrow{e}$ has, in each component, the value of the gravitational acceleration and matrix A is composed of the coefficients of $\overrightarrow{f}$ extracted from (17) for each contact point. Other details are given in ^{Oliveira et al. (2021}).

In summary, the B approach for vertical cargo stability works as follows. For each box k, obtain the contact surfaces 𝒞_k of all boxes touching it. For each surface in 𝒞_k , there are four contact points, then resulting in $4\left|{\mathcal{C}}_k\right|$ normal forces acting on box k. Next, write $4\left|{\mathcal{C}}_k\right|$ equilibrium equations, where there is one equation per contact point. Each equation contains all normal forces. These steps are repeated for all boxes in the packing to obtain the resultant system of linear inequalities (20). Next, verify if there is a solution for the system satisfying $\overrightarrow{f}\ge \overrightarrow{0}$. If no solution satisfies such a condition, the packing is unstable. There is a normal force with a negative sign, meaning that such a negative normal force is pulling up the box (i.e., causing interpenetration) instead of pushing it down to guarantee the equilibrium.

It is important to mention that with up to two contact points between the bodies, the resultant system is easily solved. However, when a body has three or more contact points, some inequalities may be redundant, and the system may have more than one solution for the contact forces. In the last case, one solution is necessary for the ε inequalities in which no value of normal force is negative to guarantee the static equilibrium.

Even in a valid solution for the system of ine qualities (i.e., that guarantees the static equilibrium), the (normal) forces may not represent the real values the boxes are receiving. Therefore, the correct (i.e., the real) value of the forces in $\overrightarrow{f}$ is obtained from solving a quadratic programming problem. This problem has a quadratic objective function and linear constraints, that is to minimize $\sum _{i=1}^{\epsilon }{f}_i{\ddot{\chi }}_i={\overrightarrow{f}}^{T}A\overrightarrow{f}-{\overrightarrow{f}}^T\overrightarrow{e}$, subject to $A\overrightarrow{f}-\overrightarrow{e}\ge \overrightarrow{0}$ and $\overrightarrow{f}\ge \overrightarrow{0}$. To be correct, besides satisfying the system of inequalities (20) (problem constraints), forces on the contact points must be conservative (problem objective function), requiring f _i or ${\ddot{\chi }}_i$ to be null on each contact point i. More details are given in ^{Baraff (1989}).

3.2 Horizontal Stability

After confirming the packing is vertically stable, we check if it is horizontally stable. The latter is made under an inertial frame of reference (i.e., outside the cargo). We model a situation where the vehicle performs a curved path (i.e., centripetal acceleration), and it exists friction between boxes. The centripetal acceleration keeps the cargo moving on the curved path. It means the packing (and vehicle) follows the road because of the centripetal force. If a non-inertial frame of reference (i.e., from the boxes or driver point of view) was considered, a centrifugal force would arise (^{Hibbeler, 2010}).

From the above assumptions, the packing is checked for static horizontal stability (^{Queiroz et al., 2019}), in which different types of forces may emerge besides those related to vertical stability. As in the support factor, we assume that forces are handled from left to right, along the x-axis, and from front to back, along the y-axis.

The first type of horizontal force is the centripetal one. The value (magnitude) of the centripetal force $f_k^c$ acting on box k is, in which r is the radius of the curve and v is the velocity in which the cargo (vehicle) is moving:

$f_k^c=m_k\frac{v^2}{r}.$ (21)

The second type of horizontal force is the lateral one acting on each contact point i of each lateral contact surface of box k. Its value is given by $f_i^{lx}$ on the direction of the x-axis and $f_i^{ly}$ on the direction of the y-axis. They are transmitted from the adjacent boxes touching box k along each of these two directions. The container’s lateral faces are used when calculating the lateral forces. It is similar to when the container’s bottom face is used in vertical stability. Note that lateral contact surfaces appear due to the lateral faces of boxes touching each other.

We apply the B approach to obtain the lateral force acting on each contact point i for each direction. The steps are the same described in the previous section about the B approach, except for the resultant mass acting on boxes obtained with the formula for force. This formula considers the resultant normal force, which considers the influence of the mass above and the gravitational acceleration. Note that lateral forces consider the resultant mass on each box (instead of the original box’s mass) because boxes may receive a parcel of masses above acting (directly and indirectly) on it to guarantee the equilibrium. The value of the resultant horizontal force acting on k is expressed by $f_k^{hx}=f_k^c+f_k^{lx}$ and $f_k^{hy}=f_k^c+f_k^{ly}$, respectively, on the direction of the x-and y-axes.

On each contact point i of the contact surface related to the boxes’ bottom and top faces (on the direction of the z-axis) emerges two other forces. One force is perpendicular to the contact surface, which is the normal force of value f _i . The other forces are tangential to the surface, which are the friction forces of values $f_i^x$ and $f_i^y$ as calculated in (22) on the direction of the x-axis and (23) on the direction of the y-axis, respectively. Following ^{Queiroz et al. (2019}), the friction force that emerges from the lateral contact surface between boxes is not considered. The authors commented that these friction forces might allow solutions that did not comply with loading/unloading operations. Besides that, a stable packing without these friction forces remains stable in their presence.

$f_i^x=\mu f_i,$ (22)

$f_i^y=\mu f_i,$ (23)

in which µ is the coefficient of static friction, and the normal and friction forces are illustrated in Figure 1.

Figure 2 illustrates box k on a surface with a coefficient of static friction. For example, if a resultant horizontal force of value $f_k^h$ is acting on the left lateral face of k, a resultant friction force of value f ^x arises as a reaction. It follows that one of the following situations may happen:

• the horizontal force can exceed the friction force, so the box slides on the surface (sliding situation); or,

• the horizontal force cannot slide the box because of the friction force. However, the horizontal force has enough intensity to rotate the box at a contact point i (tip over situation).

Figure 2
Horizontal force acting on a box for the sliding situation on the direction of the x-axis.

The friction force depends on µ, which depends on the material of the contact surfaces. The coefficient of static friction is generally greater than that of dynamic friction (^{Hibbeler, 2010}). It follows that box k will not slide on the direction of x-axis if $\mu f_k\ge f_k^{hx}$, that is $f_k^x\ge f_k^{hx}$, where f _k and $f_k^x$ are the values of the resultant normal and friction forces acting on k, respectively. Similarly, we define the condition for the slide situation on the direction of the y-axis. Note that the lower is the friction coefficient, the easier it is for a box to slide.

When a box is about to slide, the normal force is redistributed to a rotation point to avoid the tip over. The normal force avoids the rotation of boxes when torque is generated. Figure 3 shows what happens when the resultant horizontal force cannot slide a box, but it is sufficient to generate torque (i.e., the tip over situation). The resultant normal force is shifted to the bottom-right corner of the contact surface to guarantee the equilibrium. The box k then tips over if $f_k^{hx}h>f_{k}x$ holds, where x and h are the distances between the box’s center of gravity and the box’s bottom and right lateral faces, respectively. Similarly, we define the condition for the tip over situation on the direction of the y-axis.

Figure 3
A horizontal force acting on a box for the tip over situation on the direction of the x-axis.

The steps for checking static vertical and horizontal stabilities are described in Algorithm 1. In lines 1 to 5, we define the variables with each box’s mass, the set of contact surfaces, the set of contact points of each contact surface, and the normal force of each contact point. Next, lines 6 to 9 are related to the vertical stability, following Section 3.1, to obtain the correct values of the normal force of each contact point on the direction of the z-axis. Suppose the value of one normal force f _i is negative. In that case, the packing is not vertically stable because this force is pulling up instead of pushing down, and the algorithm ends.

Algorithm 1:
Checking vertical and horizontal stabilities.

In the case the packing is vertically stable, the next lines (10 to 24) of Algorithm 1 are related to horizontal stability. We first calculate the centripetal, lateral, and friction forces to check the slide and tip over situations. In line 10, we obtain the value of the resultant centripetal force $f_k^c$ acting on each box k. In lines 11 and 12 (resp., 13 and 14), we calculate the lateral force $f_i^{lx}$ (resp., $f_i^{ly}$) that is transmitted at each lateral contact point i on the direction of the x-axis (resp., y-axis). These lateral forces are obtained with the B approach by solving a quadratic programming problem similarly to vertical stability. The main difference is that now the resultant mass of box k is obtained from the formula for force $f_k=m_{k}g$ because we consider the influence of the masses above/on the top of (i.e., acting directly and indirectly) k. We assume this because each box’s weight acts along the vertical direction (z-axis), perpendicular to the lateral/horizontal forces. Still, we also need to know such action of the box along the horizontal direction (x- and y-axes). Finally, in line 15, we calculate the resultant lateral forces acting on each box along each horizontal direction.

With the centripetal and resultant lateral forces, in line 16 of Algorithm 1 we calculate the resultant horizontal forces $f_k^{hx}$ and $f_k^{hy}$ acting on each box k along each horizontal direction. Lines 17 and 18 give the friction forces on the directions of the x- and y-axes. Then, we check the slide situation in line 19 (as illustrated in Figure 2). On the other hand, for the tip over situation, we need the resultant normal force f _k acting on each box k (line 21) and the vector with the distances between the box’s center of gravity and its back-right-bottom corner (as illustrated in Figure 3). Then, we check the tip over situation in line 23. Once the packing respects all vertical and horizontal stability conditions, we assume it is statically stable.

Differently from the support factor constraints described in (8)-(10) that are directly considered in the integer model (4)-(7) for the SCLP, Algorithm 1 is used as a cutting plane routine within a branch-and-cut algorithm. The proposed algorithm is a callback routine called whenever an integer solution (packing P) is found. If P is not stable, inequality (24) is added as a lazy constraint into the integer model for the SCLP to avoid one of the boxes in P of being packed in its current position.

$\sum _{(k,p,q,r)\in P}{x}_{kpqr}\le \left|P\right|-1,$ (24)

where packing P contains the tuples (k, p, q, r) associated with the box k and its position (p, q, r) inside the packing.

4 COMPUTATIONAL EXPERIMENTS

All algorithms were coded using the C++ programming language, and computational experiments were run on a computer with an Intel Core i7 − 4790K at 4.0 GHz with 32 GB of RAM, running under the Linux Ubuntu 14.04 LTS. In the experiments, we solve the integer model (4)-(7) for the SCLP with the branch-and-cut framework available in the Gurobi Optimizer 6.5.1 with its default parameters. The COIN-OR IPOPT 3.12.6 is used to solve the quadratic programming problems that appear in Algorithm 1. For the support factor (i.e., JMY approach), we solve directly the integer model with constraints (8)-(10). On the other hand, the proposed algorithm for horizontal cargo stability is used as a callback routine that returns lazy constraints (24) in unstable packings.

Due to the static horizontal stability assumptions, an optimal solution is attending a scenario where the vehicle moves on a curved path. We consider a curve of type “mountainous” from Class II for which the vehicle is limited to 40 km/h. The curve radius is 30 m (^{Departamento Nacional de Estradas de Rodagem, 1973}). Boxes are made of cardboard, for which we use the coefficient of static friction $\mu =0.25$. Preliminary experiments on other parameters have pointed to the same conclusions given next.

The experiments are conducted on the benchmark instances from ^{Junqueira et al. (2012b)}, which have $n=5$, 10, and 20 types of boxes. Containers have dimensions $L=W=H$, generating groups for dimensions equal to 10, 20 and 30. Each group has 10 different instances. Two groups of instances are used: A _n considers the dimensions of the boxes ranging from 25% to 75% of the container’s respective dimension; and B _n considers the dimensions ranging from 10% to 50% of the container’s respective dimension. Instances of the group A _n have predominantly large boxes compared to B _n , implying instances of the latter more challenging to solve. The maximum number of copies of each item k is randomly picked in the interval $\left[1,\lfloor \frac{L}{l_k}\rfloor ·\lfloor \frac{W}{w_k}\rfloor ·\lfloor \frac{H}{h_k}\rfloor \right]$, for $k=1,2,\dots ,n$.

There is a total of 180 instances, divided into 18 classes of 10 instances each. The time limit the branch-and-cut algorithm has to solve each instance is 3600 seconds. The following tables present the minimum, average, and maximum of the gap (in %) calculated by Gurobi, the runtime (in seconds), and the container fill rate (in %), for each group of 10 instances with the same container dimensions. For the proposed approach, we also present the number of cuts (unstable packings) generated. On the other hand, we do not show the total time is spent because this time is limited by the number of cuts and is less than 0.01 per call. It is essential to mention that these values are presented for instances with a feasible (optimal) solution. Columns with the entry “-” mean the time limit was exceeded, with no feasible solution. The discussion and comparison of results consider the average values if it is not stated otherwise. For results concerning only vertical stability, we refer to ^{Oliveira et al. (2021}).

In general, the proposed approach for horizontal stability has presented stable solutions with better container fill rates than the support factor (i.e., the JMY approach). This is demonstrated in the results for instances with predominantly large boxes, A _n in Table 1, and instances with mostly small boxes, B _n in Table 2. These results indicate how less restrictive the proposed approach is since it does not impose particular support for boxes’ faces. We noticed the JMY approach had problems handling instances with more types of boxes, presenting no feasible solution within the imposed time limit. The number of cuts added with the proposed approach is larger for groups with fewer types of boxes or the predominance of small boxes.

Table 1
Horizontal stability results for the groups A ₅, A ₁₀, and A ₂₀.

			Unitary Discretization							CD							RRP
Group	Dim.	Type	JMY approach			Proposed approach				JMY approach			Proposed approach				JMY approach			Proposed approach
			gap	Time	Vol.	#Cuts	gap	Time	Vol.	gap	Time	Vol.	#Cuts	gap	Time	Vol.	gap	Time	Vol.	#Cuts	gap	Time	Vol.
			(%)	(s)	(%)		(%)	(s)	(%)	(%)	(s)	(%)		(%)	(s)	(%)	(%)	(s)	(%)		(%)	(s)	(%)
A 5	10	Min.	0.00	0.10	76.00	0.00	0.00	0.20	77.20	0.00	0.10	76.00	0.00	0.00	0.07	77.20	0.00	0.05	76.00	0.00	0.00	0.06	77.20
		Avg.	0.00	6.15	88.75	50.30	1.11	391.22	90.04	0.00	2.74	88.75	58.90	0.00	15.14	90.74	0.00	2.35	88.75	58.90	0.00	17.99	90.74
		Max.	0.00	28.84	100.00	402.00	10.10	3600.08	100.00	0.00	10.58	100.00	402.00	0.00	92.60	100.00	0.00	12.82	100.00	402.00	0.00	124.72	100.00
	20	Min.	0.00	1.58	53.25	4.00	0.00	47.61	31.40	0.00	0.02	53.25	0.00	0.00	0.10	6.12	0.00	0.01	53.25	4.00	0.00	0.08	51.05
		Avg.	0.00	17.01	72.20	2125.40	14.30	2109.86	72.71	0.00	2.26	72.20	3130.80	20.60	2542.77	68.92	0.00	1.07	68.62	4641.40	12.97	2352.75	70.78
		Max.	0.00	76.84	87.00	10008.00	68.60	3601.37	90.10	0.00	15.40	87.00	17477.00	93.87	3603.93	90.10	0.00	1.07	68.62	28028.00	43.90	3632.73	88.25
	30	Min.	0.00	2.87	42.49	4.00	0.00	2568.87	9.48	0.00	0.00	42.49	4.00	0.00	0.13	29.67	0.00	0.00	42.49	4.00	0.00	0.09	34.20
		Avg.	0.00	7.64	58.58	15.00	34.60	3288.20	56.51	0.00	0.03	58.58	1211.30	12.91	1483.04	52.99	0.00	0.02	55.51	5054.90	8.45	1715.88	54.04
		Max.	0.00	16.12	76.18	28.00	90.50	3607.64	76.18	0.00	0.10	76.18	7769.00	41.90	3601.53	74.51	0.00	0.10	66.70	25885.00	35.00	3604.56	66.70
A 10	10	Min.	0.00	2.35	96.80	0.00	0.00	0.33	98.40	0.00	0.62	96.80	0.00	0.00	0.12	98.40	0.00	0.84	96.80	0.00	0.00	0.13	98.40
		Avg.	0.00	214.63	99.22	3.40	0.00	27.52	99.53	0.00	225.53	99.21	2.20	0.00	29.56	99.53	0.00	221.73	99.20	2.20	0.00	31.21	99.53
		Max.	0.00	1268.58	100.00	12.00	0.00	193.64	100.00	0.00	1232.70	100.00	4.00	0.00	235.2	100.00	0.00	1208.86	100.00	4.00	0.00	252.34	100.00
	20	Min.	0.00	13.26	70.76	0.00	0.00	325.74	4.90	0.00	1.93	70.76	0.00	0.00	5.99	49.88	0.00	1.95	71.76	0.00	0.00	1.73	18.50
		Avg.	2.14	789.42	85.14	4.20	19.51	1966.89	77.19	0.00	190.96	85.26	4.10	3.60	1532.41	81.71	0.00	125.05	82.28	7.00	18.85	1723.18	75.96
		Max.	15.50	3600.07	94.30	12.00	97.00	3601.22	94.29	0.00	1088.87	94.30	17.00	21.60	3603.41	94.30	0.00	582.21	94.30	26.00	81.50	3602.53	94.30
	30	Min.	0.00	68.20	63.83	4.00	0.00	3409.98	5.40	0.00	0.45	63.83	4.00	0.00	0.33	58.32	0.00	1.08	61.89	4.00	0.00	0.03	45.47
		Avg.	0.00	683.61	76.02	4.00	41.77	3658.79	57.54	0.00	53.13	76.02	4.00	3.27	1774.98	79.04	0.00	17.84	73.83	4.00	6.28	1291.41	77.50
		Max.	0.00	3601.33	93.56	4.00	94.60	3954.21	93.56	0.00	382.04	93.56	4.00	17.50	3600.92	93.56	0.00	136.95	93.56	4.00	40.80	3601.68	93.56
A 20	10	Min.	-	3600.00	-	0.00	0.00	1.84	99.60	0.00	16.86	98.40	0.00	0.00	1.70	99.60	0.00	16.97	98.40	0.00	0.00	1.70	99.60
		Avg.	-	3600.00	-	2.20	0.00	5.34	99.95	0.12	379.61	99.81	2.20	0.00	3.82	99.95	0.12	379.58	99.81	2.20	0.00	3.82	99.95
		Max.	-	3600.00	-	7.00	0.00	19.36	100.00	1.20	3601.01	100.00	7.00	0.00	11.30	100.00	1.20	3601.01	100.00	7.00	0.00	11.31	100.00
	20	Min.	-	3600.00	-	0.00	0.00	1588.51	19.25	0.00	37.24	100.00	0.00	0.00	11.33	90.64	0.00	183.03	81.93	0.00	0.00	79.78	90.64
		Avg.	-	3600.00	-	3.60	19.97	3123.07	79.14	0.47	1172.69	91.76	3.20	0.00	792.16	95.02	0.46	1021.18	91.53	3.20	1.36	730.02	95.01
		Max.	-	3600.00	-	4.00	80.70	3602.28	97.30	4.70	3601.03	97.75	4.00	0.00	3601.04	97.75	4.60	3601.12	97.75	4.00	0.00	3600.04	97.75
	30	Min.	-	3600.00	-	4.00	0.00	3257.34	7.70	0.00	51.63	77.96	4.00	0.00	37.24	86.24	0.00	27.09	77.96	4.00	0.00	21.18	86.42
		Avg.	-	3600.00	-	4.00	23.09	3643.17	75.58	0.28	806.01	88.51	4.00	1.27	1107.72	91.59	0.00	489.50	87.15	4.00	0.23	2366.33	91.27
		Max.	-	3600.00	-	4.00	92.30	4215.76	97.21	2.80	3601.03	92.89	4.00	2.70	2500.44	97.21	0.00	2663.68	92.89	4.00	1.20	14194.00	97.21

Table 2
Horizontal stability results for the groups B ₅, B ₁₀, and B ₂₀.

			Unitary Discretization							CD							RRP
Group	Dim.	Type	JMY approach			Proposed approach				JMY approach			Proposed approach				JMY approach			Proposed approach
			gap	Time	Vol.	#Cuts	gap	Time	Vol.	gap	Time	Vol.	#Cuts	gap	Time	Vol.	gap	Time	Vol.	#Cuts	gap	Time	Vol.
			(%)	(s)	(%)		(%)	(s)	(%)	(%)	(s)	(%)		(%)	(s)	(%)	(%)	(s)	(%)		(%)	(s)	(%)
B 5	10	Min.	0.00	4.69	95.20	0.00	0.00	0.55	95.20	0.00	1.35	95.20	0.00	0.00	0.63	95.20	0.00	1.36	95.20	0.00	0.00	0.38	95.20
		Avg.	0.14	751.97	98.50	6653.10	0.00	468.01	98.59	0.14	853.12	98.51	804.90	0.00	447.45	98.59	0.14	754.74	98.51	833.20	0.00	435.36	98.59
		Max.	0.80	3600.34	100.00	50063.00	0.00	3600.01	100.00	0.80	3600.06	100.00	7062.00	0.00	3601.67	100.00	0.80	3600.01	100.00	7062.00	0.00	3600.00	100.00
	20	Min.	0.00	204.57	58.50	45.00	0.00	108.52	77.60	0.00	32.74	69.65	61.00	0.00	11.53	77.60	0.00	13.60	69.05	281.00	0.00	2.87	48.60
		Avg.	14.56	3260.49	83.36	13203.80	3.70	2553.87	95.03	11.36	3154.16	85.87	73058.20	0.67	2358.07	94.47	11.49	3081.18	85.23	72782.10	13.03	2301.38	90.53
		Max.	32.70	3600.13	99.20	46127.00	2.20	3600.62	100.00	30.50	3600.12	99.20	162107.00	6.70	3609.72	100.00	30.90	3600.14	99.20	160379.00	92.00	3613.29	100.00
	30	Min.	2.90	3600.18	24.69	115.00	2.00	1031.63	84.18	3.10	3600.01	24.69	59.00	0.00	1128.52	47.50	0.00	971.14	49.01	59.00	0.00	946.94	11.61
		Avg.	18.51	3601.29	78.97	7026.20	10.40	3600.54	88.95	16.50	3601.32	79.36	64928.60	7.27	3182.20	92.72	18.63	3217.22	79.19	76743.40	26.22	3335.01	71.42
		Max.	75.30	3603.98	95.64	55356.00	27.00	4215.76	97.21	75.30	3610.04	95.38	171815.00	15.81	3600.5	99.00	82.70	3606.36	95.64	168289.00	88.30	3601.62	99.00
B 10	10	Min.	0.00	10.53	9 9.90	0.00	0.00	1.23	100.00	0.00	10.52	100.00	0.00	0.00	1.22	100.00	0.00	10.47	99.90	0.00	0.00	1.22	100.00
		Avg.	0.00	59.10	99.98	7.70	0.00	4.20	100.00	0.00	60.00	100.00	13.00	0.00	4.23	100.00	0.00	60.13	99.98	11.80	0.00	4.19	100.00
		Max.	0.00	273.47	100.00	22.00	0.00	7.07	100.00	0.00	276.49	100.00	56.00	0.00	7.06	100.00	0.00	278.84	100.00	56.00	0.00	7.06	100.00
	20	Min.	4.30	3600.12	60.00	42.00	0.00	3313.05	73.60	10.40	3600.13	60.00	32.00	0.00	1838.46	1.50	4.30	3600.92	60.00	45.00	0.00	2291.8	3.00
		Avg.	19.91	3602.11	81.73	12787.00	9.50	3571.42	90.07	24.48	3602.17	81.08	35940.30	18.28	3424.05	81.71	23.79	3603.25	81.01	33321.50	19.00	3469.40	81.95
		Max.	40.00	3605.10	95.63	54188.00	26.00	3600.64	100.00	40.30	3605.75	95.63	148589.00	98.50	3601.08	100.00	23.79	3603.25	81.01	159543.00	97.00	3600.63	100.00
	30	Min.	19.60	3660.61	68.47	13.00	9.30	3600.01	79.81	19.60	3622.96	68.47	46.00	9.30	3600.00	79.81	19.60	3606.69	68.47	7.00	9.00	3600.00	80.74
		Avg.	23.95	3724.07	76.00	270.10	15.17	3603.43	84.83	23.95	3646.85	76.00	42284.60	15.14	3600.01	84.83	25.58	3632.47	74.37	27184.60	15.02	3600.00	84.92
		Max.	31.50	3865.85	80.36	2295.00	20.19	3612.10	90.70	31.50	3671.19	80.36	136988.00	20.10	3600.03	90.70	31.50	3668.39	80.36	142789.00	19.20	3600.01	90.70
B 20	10	Min.	0.00	32.06	100.00	0.00	0.00	0.72	100.00	0.00	32.13	100.00	0.00	0.00	0.72	100.00	0.00	32.15	100.00	0.00	0.00	0.72	100.00
		Avg.	0.00	87.24	100.00	3.50	0.00	4.79	100.00	0.00	87.18	100.00	4.60	0.00	4.79	100.00	0.00	87.35	100.00	4.60	0.00	4.81	100.00
		Max.	0.00	156.05	100.00	9.00	0.00	8.07	100.00	0.00	156.26	100.00	16.00	0.00	8.06	100.00	0.00	157.20	100.00	16.00	0.00	8.08	100.00
	20	Min.	0.50	3628.63	82.99	24.00	0.30	3600.01	92.03	0.50	3612.72	82.99	6.00	0.30	3600.00	92.03	0.50	3612.30	82.99	14.00	0.30	3600.00	92.03
		Avg.	6.89	3630.96	93.08	296.40	3.78	3600.13	96.18	6.89	3629.82	93.08	3986.70	3.78	3600.02	96.18	6.89	3631.28	93.08	4301.00	3.98	3600.22	96.18
		Max.	17.00	3635.48	99.45	1825.00	7.90	3600.37	99.70	17.00	3672.34	99.45	21945.00	7.90	3600.44	99.70	17.00	3678.31	99.45	24174.00	7.90	3600.72	99.70
	30	Min.	-	-	-	5.00	11.60	3600.18	81.46	-	3600.00	-	7.00	11.60	3600.01	81.46	-	3600.00	-	6.00	11.60	3600.01	81.46
		Avg.	20.71	1689.04	59.29	13.70	14.63	3608.24	85.31	-	3600.00	-	3751.10	14.72	3600.02	85.27	-	3600.00	-	14117.00	14.62	3600.02	85.32
		Max.	33.50	5176.78	85.40	23.00	18.50	3667.17	88.32	-	3600.00	-	24952.00	18.53	3600.05	88.32	-	3600.00	-	83909.00	18.5	3600.06	88.32

In Table 1, for group A _n , the JMY approach has better container fill rates only when the gap of the proposed approach is worse (e.g., groups A ₅ and A ₁₀ with unitary discretization and containers dimensions equal to 30). In all groups whose gap is equal for both approaches, the JMY does not outperform the proposed approach, which in turn allows the container fill rates to increase about 2 percentage points in groups like A ₅ with RRP and containers dimensions equal to 10. With a more refined grid as the RRP, gaps are overall smaller, and better fill rates are achieved. On the other hand, such a grid may decrease the container fill rates, as in the results of the JMY approach for groups A ₅ and A ₁₀ with containers dimensions equal to 10 and 20. Refined grids are also associated with a high number of cuts to avoid unstable packings since the number of possibilities to pack boxes reduces, as in the results of the proposed approach for A ₅ with CD and RRP and containers dimensions equal to 20 and 30.

The results of Table 1 indicates the proposed approach has the best container fill rates for groups: A ₅, with CD or RRP and containers in size 10, and with unitary discretization and containers in size 20; A ₁₀, with all grids and containers in size 10, and with CD and containers in size 30; A ₂₀, with all grids and containers in size 10, and with CD and containers in sizes 20 and 30. In the others, the JMY approach is superior, although it could not present any solution for group A ₂₀ with unitary discretization. The JMY is superior in groups with fewer types of boxes and/or smaller containers in terms of the gap. The number of cuts added with the proposed approach is higher in groups with fewer types of boxes. On the other hand, the proposed approach obtains stable solutions with better fill rates from groups having more boxes.

Table 2 has results of group B _n whose conclusions are similar to those of A _n . In general, stable solutions tend to worsen, and more cuts are added in the presence of refined grids, although the gaps are better with such grids. Besides that, container fill rates tend to increase in groups with more types of boxes without worsening gaps in many cases. For B ₅, with unitary discretization and containers in sizes 10, 20, and 30, the proposed approach provides stable solutions with the best fill rates, with a difference of 0.10%, 12.28%, and 11.21% in comparison with the JMY, respectively. With CD and RRP, the best fill rates are again obtained with the proposed approach, except with RRP and containers in size 30 (but here, the proposed approach’s gap is worse). Similar results are obtained from B ₁₀ and B ₂₀, where the proposed approach has the best fill rates for all grids.

The results of groups A _n and B _n validate the ability of the proposed approach to allow stable solutions with larger container fill rates. We notice that group B _n is harder to solve due to the presence of small boxes, resulting in grids with more points. This is perceptible on the number of cuts and gaps in Table 2, which are overall larger for A _n in Table 1. Concerning the number of cuts, while the maximum is 28028 for group A ₅ with RRP, this value reaches 171815 for group B ₅ with CD. On the other hand, the container fill rates are superior in B _n compared to A _n . We notice that in B _n , the stable solutions with the proposed approach have overall larger container fill rates for all grids and containers dimensions.

The proposed approach has the advantage of not imposing support, which in turn is hard to determine. Thinking on that, we present Table 3 with the support provided with the proposed approach for horizontal stability, which may help decision-makers and practitioners in their decisions. We show the minimum, average, and maximum support for boxes’ faces of each group and container dimensions. This may be useful in practical situations when requiring solutions to be stable by satisfying a given support factor. The support comes from three directions: z-axis (bottom faces); x-axis (left faces); and, y-axis (frontal faces). The bottom faces resting on the floor, and the left and frontal faces in direct contact with the container’s lateral faces are not considered.

Table 3
Average support of the solutions obtained with the proposed approach.

		Horizontal Stability
Group	Dim.	On the direction of z-axis			On the direction of x-axis			On the direction of y-axis
		Min	Avg.	Max.	Min.	Avg.	Max.	Min.	Avg.	Max.
A 5	10	75.00	98.82	100.00	70.00	98.89	100.00	81.67	99.17	100.00
	20	56.10	95.47	100.00	53.38	94.94	100.00	53.80	96.13	100.00
	30	58.81	94.50	100.00	37.30	93.05	100.00	48.43	94.52	100.00
A 10	10	69.50	91.78	100.00	65.79	90.70	100.00	66.83	90.61	100.00
	20	59.76	78.30	93.98	41.78	70.28	90.36	38.69	67.05	91.82
	30	51.58	76.70	96.06	68.38	79.63	83.38	67.15	86.56	100.00
A 20	10	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00
	20	54.28	97.53	100.00	54.41	97.17	100.00	69.84	97.40	100.00
	30	36.78	93.14	100.00	21.01	91.58	100.00	25.83	89.50	100.00
B 5	10	94.67	99.41	100.00	83.29	98.38	100.00	94.17	99.69	100.00
	20	62.01	87.86	100.00	40.05	80.92	100.00	68.48	90.13	100.00
	30	58.77	81.67	96.46	43.94	72.56	90.47	52.68	74.79	87.82
B 10	10	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00
	20	34.50	97.17	100.00	43.96	97.04	100.00	46.25	97.41	100.00
	30	27.85	91.54	100.00	29.14	91.16	100.00	33.23	90.62	100.00
B 20	10	96.67	99.74	100.00	98.89	99.91	100.00	96.67	99.78	100.00
	20	64.98	91.81	100.00	63.38	89.56	97.78	69.53	91.72	100.00
	30	74.56	91.22	98.67	65.01	90.86	100.00	62.93	85.62	98.97

In Table 3, the average support for groups A _n and B _n are, respectively: 91.81% and 93.38%, for bottom faces; 90.69% and 91.15%, for left faces; and, 91.22% and 92.20%, for frontal faces. When comparing the groups, the support is smaller for A _n . Note that this group has larger boxes compared to B _n and has stable solutions with inferior container fill rates as shown in Tables 1 and 2. The difference is small in terms of average support compared to the JMY approach that uses full support. However, it is quite large when considering the column with the minimum support in each direction. For instance, the minimum values are below 40% in groups A ₂₀ and B ₁₀ with containers in size 30, and still, the container fill rates are better compared to JMY.

5 CONCLUDING REMARKS

In this paper, we handle the SCLP with cargo stability. We solve an integer model for the container loading problem while a routine generates cutting planes that prevent unstable packings. This routine verifies cargo stability by using an approach based on the static equilibrium of rigid bodies. The approach considers two situations in the presence of horizontal forces, sliding and tipping over.

Even though the proposed approach demands the correct value of the normal forces, which requires solving a quadratic programming model, it does not impose boxes on being side-by-side and touching themselves as demanded by the support factor approach. When considering all grids, the average results in terms of container fill rates for groups A ₅, A ₁₀, A ₂₀, B ₅, B ₁₀, and B ₂₀, are 72.44%, 86.24%, 62.06%, 87.50%, 85.57%, and 70.95%, with the support factor, and 71.95%, 83.06%, 91.94%, 92.10%, 89.81%, and 92.62%, with the proposed approach, respectively. Our approach allowed stable solutions with higher container fill rates, especially for groups in which small boxes predominate.

With the proposed approach for horizontal stability, the obtained stable solutions have container fill rates with the unitary discretization, CD, and RRP equal to 78.69%, 84.39%, and 83.87%, for group A _n , and 93.22%, 92.64%, and 88.67%, for B _n , respectively. In general, the more restricted grids allowed to find better solutions within a shorter computing time, but the results showed that the more restricted grid (RRP) does not necessarily lead to better results than a less restricted grid (CD). Besides that, instances with smaller and fewer types of boxes are susceptible to have a high number of cuts to avoid unstable solutions. In terms of the support obtained from the stable solutions with the proposed approach, the average value is superior to 90% in all three directions (x-, y-, and z-axes). On the other hand, the minimum support is inferior to 70% on average, allowing more flexibility than full support.

Future research directions are to consider other forces that may occur during transportation, as the influence of wind pressure, road curvature, and abrupt speed changes. The more we add details in the approach, the more it represents the real cargo stability. Another direction is to investigate the impact of other/multiple practical constraints cooperatively with cargo stability, for example, load balancing, load-bearing, and multi-drop requirements, besides coupling them in vehicle routing problems. In this way, the final cargo will comply with many aspects required for efficient transportation. It is also of interest the proposal of formulations that may have parts of the proposed approach as linear constraints similar to the support factor constraints.

Acknowledgements

The authors would like to thank the financial support of the National Council of Technological and Scientific Development (CNPq) [grants numbers 311499/2014-7, 308312/2016-3, 425340/2016-3, 314366/2018-0, 311185/2020-7], State of Goiás Research Foundation (FAPEG), and State of São Paulo Research Foundation (FAPESP) [grants numbers 2015/11937-9, 2016/01860-1, 2017/11831-1].

References

• ¹
  ALMEIDA CUNHA JG, DE LIMA VL & QUEIROZ TA. 2020. Grids for cutting and packing problems: a study in the 2D knapsack problem. 4OR A Quarterly Journal of Operations Research, 18: 293-339.
• ²
  ARAYA I, GUERRERO K & NEZ EN. 2017. VCS: A new heuristic function for selecting boxes in the single container loading problem. Computers & Operations Research, 82: 27-35.
• ³
  ARAYA I & RIFF MC. 2014. A beam search approach to the container loading problem. Computers & Operations Research , 43: 100-107.
• ⁴
  BARAFF D. 1989. Analytical methods for dynamic simulation of non-penetrating rigid bodies. Computer Graphics, 23(3): 223-232.
• ⁵
  BORTFELDT A & GEHRING H. 2001. A hybrid genetic algorithm for the container loading problem. European Journal of Operational Research, 131(1): 143-161.
• ⁶
  BRACHT EC, QUEIROZ TA, SCHOUERY RCS & MIYAZAWA FK. 2016. Dynamic Cargo Stability in Loading and Transportation of Containers. In: 2016 IEEE International Conference on Automation Science and Engineering (CASE). Fort Worth, TX, USA.
• ⁷
  CHEN C, LEE SM & SHEN Q. 1995. An analytical model for the container loading problem. European Journal of Operational Research , 80(1): 68-76.
• ⁸
  CHRISTENSEN SG & ROUSOE DM. 2009. Container loading with multi-drop constraints. International Transactions in Operational Research, 16(6): 727-723.
• ⁹
  DEPARTAMENTO NACIONAL DE ESTRADAS DE RODAGEM. 1973. Normas para o projeto das Estradas de Rodagem (Brasil). Tech. rep.. Ministério dos Transportes. Rio de Janeiro.
• ¹⁰
  EGEBLAD J, GARAVELLI C, LISI S & PISINGER D. 2010. Heuristics for container loading of furniture. European Journal of Operational Research , 200(3): 881-892.
• ¹¹
  EGEBLAD J & PISINGER D. 2009. Heuristic approaches for the two-and threedimensional knapsack packing problem. Computers & Operations Research , 36(4): 1026-1049.
• ¹²
  FANSLAU T & BORTFELDT A. 2010. Solving container loading problems by block arrangement. INFORMS Journal on Computing, 22(2): 222-235.
• ¹³
  FASANO G. 2004. A MIP approach for some practical packing problems: Balancing constraints and tetris-like items. Quarterly Journal of the Belgian, French and Italian Operations Research Societies, 2(2): 161-174.
• ¹⁴
  HERZ JC. 1972. Recursive computational procedure for two-dimensional stock cutting. Journal of Research and Development, 16(5): 462-469.
• ¹⁵
  HIBBELER RC. 2010. Statics & Mechanics of Materials. New Jersey, USA: Prentice Hall.
• ¹⁶
  HIFI M. 2004. Exact algorithms for unconstrained three-dimensional cutting problems: a comparative study. Computers & Operations Research , 31(5): 657-674.
• ¹⁷
  HIFI M, KACEM I, NÈGRE S & WU L. 2010. A linear programming approach for the three-dimensional bin-packing problem. Electronic Notes in Discrete Mathematics, 36: 993-1000.
• ¹⁸
  JAMRUS T & CHIEN CF. 2016. Extended priority-based hybrid genetic algorithm for the less-than-container loading problem. Computers & Industrial Engineering, 96: 227-236.
• ¹⁹
  JUNQUEIRA L, MORABITO R & YAMASHITA DS. 2012a. MIP-based approaches for the container loading problem with multi-drop constraints. Annals of Operations Research, 199(1): 51-75.
• ²⁰
  JUNQUEIRA L, MORABITO R & YAMASHITA DS. 2012b. Three-dimensional container loading models with cargo stability and load bearing constraints. Computers & Operations Research , 39(1): 74-85.
• ²¹
  JUNQUEIRA L & QUEIROZ TA. 2022. The static stability of support factor-based rectangular packings: an assessment by regression analysis. International Transactions in Operational Research , 29: 574-599.
• ²²
  KURPEL DV, SCARPIN CT, PÉCORA JE, SCHENEKEMBERG CM & COELHO LC. 2020. The exact solutions of several types of container loading problems. European Journal of Operational Research , 284(1): 87-107.
• ²³
  MACK D, BORTFELDT A & GEHRING H. 2004. A parallel hybrid local search algorithm for the container loading problem. International Transactions in Operational Research , 11(5): 511-533.
• ²⁴
  MARTELLO S, PISINGER D & VIGO D. 2000. The three-dimensional bin packing problem. Operations Research, 48(2): 256-267.
• ²⁵
  NASCIMENTO OX, QUEIROZ TA & JUNQUEIRA L. 2021. Practical constraints in the container loading problem: Comprehensive formulations and exact algorithm. Computers & Operations Research , 128: 105186.
• ²⁶
  OLIVEIRA LA, DE LIMA VL, QUEIROZ TA & MIYAZAWA FK. 2021. The container loading problem with cargo stability: a study about support factor, mechanical equilibrium, and grids. Engineering Optimization, 53: 1192-1211.
• ²⁷
  PAQUAY C, SCHYNS M & LIMBOURG S. 2016. A mixed integer programming formulation for the three-dimensional bin packing problem deriving from an air cargo application. International Transactions in Operational Research , 23(1-2): 187-213.
• ²⁸
  QUEIROZ TA, BRACHT EC, MIYAZAWA FK & BITTENCOURT ML. 2019. An extension of Queiroz and Miyazawa’s method for vertical stability in two-dimensional packing problems to deal with horizontal stability. Engineering Optimization , 56(6): 1049-1070.
• ²⁹
  QUEIROZ TA & MIYAZAWA FK. 2014. Order and static stability into the strip packing problem. Annals of Operations Research , 223(1): 137-154.
• ³⁰
  RAMOS AG, OLIVEIRA JF, GONÇALVES JF & LOPES MP. 2015. Dynamic stability metrics for the container loading problem. Transportation Research Part C: Emerging Technologies, 60: 480-497.
• ³¹
  RAMOS AG, OLIVEIRA JF & LOPES MP. 2016. A physical packing sequence algorithm for the container loading problem with static mechanical equilibrium conditions. International Transactions in Operational Research , 23(1-2): 215-238.
• ³²
  SCHEITHAUER G & TERNO J. 1996. The g4-heuristic for the pallet loading problem. Journal of the Operational Research Society, 47: 511-522.
• ³³
  SHENG L, XIUQIN S, CHANGJIAN C, HONGXIA Z, DAYONG S & FEIYUE W. 2017. Heuristic algorithm for the container loading problem with multiple constraints. Computers & Industrial Engineering , 108: 149-164.
• ³⁴
  SILVA E, AO AL, TOLEDO F & WAUTERS T. 2020. A matheuristic framework for the Three-dimensional Single Large Object Placement Problem with practical constraints. Computers & Operations Research , 124: 105058.
• ³⁵
  SILVA EF, TOFFOLO TAM & WAUTERS T. 2019. Exact methods for three-dimensional cutting and packing: A comparative study concerning single container problems. Computers & Operations Research , 109: 12-27.
• ³⁶
  SILVA JLC, SOMA NY & MACULAN N. 2003. A greedy search for the three-dimensional bin packing problem: the packing static stability case. International Transactions in Operational Research , 10(2): 141-153.
• ³⁷
  TAO Y & WANG F. 2015. An effective tabu search approach with improved loading algorithms for the 3L-CVRP. Computers & Operations Research , 55: 127-140.
• ³⁸
  TSAI JF, WANG PC & LIN MH. 2015. A global optimization approach for solving threedimensional open dimension rectangular packing problems. Optimization, 64(12): 2601- 2618.

Publication Dates

• Publication in this collection
  11 Oct 2021
• Date of issue
  2021

History

• Received
  21 Feb 2020
• Accepted
  21 Oct 2020