1 Introduction

The layout design of a satellite module involves placing a certain number of objects, including various instruments and devices, in a particular satellite module, while satisfying various constraints with specific objectives. Stability and service life are essential requirements of a successful layout design, and other performances of the whole satellite module system should also be considered.

Due to engineering and mathematics complexities, satellite module layout design is known as an NP-hard problem. It stems from the bin packing problem (Crainic et al., 2011; Gonçalves and Resende, 2011; Khanafer et al., 2012), which can mainly be classified into two categories according to the study subjects: (1) 2D bin packing problem; (2) 3D bin packing problem. For the 2D bin packing problem, Gonçalves and Resende (2011) proposed a parallel multi-population genetic algorithm that hybridizes a novel placement procedure with a genetic algorithm based on random keys. Khanafer et al. (2012) proposed a tree decomposition method based heuristic strategy to decompose the original problem into several sub-problems. Li et al. (2012) proposed a heuristic particle swarm optimization approach with a quasi-human strategy. Zhang et al. (2005) proposed a hybrid heuristic algorithm based on divide-and-conquer and greedy strategies for the 2D rectangular packing problem. Liu et al. (2009) proposed an improved energy landscape paving (ELP) method by incorporating a new configuration update mechanism into the ELP method to solve the circular packing problem. For the 3D bin packing problem, Martello et al. (2000) proved that the asymptotical worst-case performance of the continuous lower bound is 1/8, and presented an exact branch-and-bound algorithm. Moon and Nguyen (2014) presented a mixed integer programming model which features upper and lower bounds for the 3D bin packing problem. To deal with the 3D parallelepiped bin packing problem, Allen et al. (2011) proposed a hybrid placement strategy, and the test results showed that it outperforms other methods from the literature.

In this paper, we study the layout design of a satellite module with performance constraints that can be considered a special case of the bin packing problem mentioned above. Various methods have been proposed to deal with the layout design of satellite modules. Sun and Teng (2003) proposed a two-stage layout method, the centripetal balancing method, for global layout design in the first stage, and an ant colony optimization algorithm for detailed layout design of a satellite module in the second stage. Zhang et al. (2008) presented a hybrid method by combining soft computing techniques, including the Hopfield neural network, genetic algorithm/particle swarm optimization, and quasi-principal component analysis. Teng et al. (2010) developed a dual-system framework based on the cooperative coevolutionary genetic algorithm with the merits of great population diversity and a decrease in premature convergence. He et al. (2013) proposed a quasi-physical algorithm based on coarse and fine adjustment. Tang and Teng (1999) presented a decimal coded adaptive genetic algorithm, which decreases combinatorial explosion and premature convergence. Wang and Teng (2009) proposed a knowledge fusion approach, which harnesses the potential of both humans and computers through evolutionary computation. Liu and Teng (2008) presented a human algorithm knowledge layout design (HAKD) method, which fuses human intelligence, computer intelligence (evolution algorithm), and prior knowledge (relevant layout diagrams) at the gene level of the evolution algorithm. Huo and Teng (2009) optimized a layout design using four steps. First, the whole layout problem was decomposed into several sub-layout problems. Second, a relaxation model was adopted to distribute all objects among subspaces. Third, a coevolutionary genetic algorithm was adopted to solve sub-layout problems. Finally, a heuristic combination-rotation method was adopted to adjust the constraints to obtain the final layout. Lei and Qiu (2006) presented a novel adaptive particle swarm optimizer based on multi-modified strategies, which can not only escape from the attraction of local optima, but also maintain the characteristic of a fast speed search in the early convergence phase. Zhou et al. (2005) presented a constraint handling strategy suitable for particle swarm optimization, and used direct search to intensify its local search ability. Jin and Teng (2007) proposed a case retrieval algorithm for reusing previously stored design solutions and optimizing layout design using prior knowledge and an evolutionary approach. Chen et al. (2008) presented an improved differential evolution algorithm, which makes random microperturbations to the location of components to obtain a better layout scheme. Zhang et al. (2013) proposed a parallel dual-system cooperative coevolutionary differential evolution algorithm with human-computer cooperation. Xu and Xiao (2008) presented an ant colony optimization algorithm based on a step-by-step positioning technique. Liu et al. (2011) proposed a heuristic algorithm based on tabu search, where tabu search was used to jump out of the local minima. Liu and Li (2010) proposed the basin filling algorithm by combining the improved energy landscape paving method, the gradient method, and the configuration update strategy. Liu et al. (2010) proposed a heuristic simulated annealing algorithm by incorporating the heuristic neighborhood search mechanism and the adaptive gradient method into the simulated annealing procedure.

Although the approaches mentioned above enable effective layout design of the satellite module, their efficiency still needs to be improved. In fact, a reasonable combination of the stochastic algorithm with global search ability, the local accurate search method, and some heuristic strategies may be an effective way to construct a high-performance algorithm for layout optimization of satellite modules. The Wang-Landau (WL) sampling algorithm (Wang and Landau, 2001; Landau et al., 2004) is an improved Monte Carlo algorithm. Unlike conventional Monte Carlo simulations that generate a probability distribution at a given temperature, the WL sampling method can estimate the density of states accurately via a random walk, which produces a flat histogram in the energy space. There have been many improvements on and applications of the WL sampling algorithm. For example, Zhou and Bhatt (2005) proved the convergence of the WL sampling algorithm, and found that the histogram increases uniformly with small fluctuations after a stage of initial accumulation. Seaton et al. (2010) used it to describe the thermodynamic behavior of a continuous homopolymer. Schulz et al. (2003) proposed a simple modification of the WL sampling algorithm. This modification removes the systematic error that occurs at the boundary of the range of energy over which the random walk takes place. The greatest challenge of the layout optimization of the satellite module is that the function to be optimized is characterized by a multitude of local minima separated by high-energy barriers. The WL sampling method can visit all the accessible states of the system, which means that it can jump out of these high-energy barriers. Therefore, the WL sampling method is an ideal global search algorithm for layout optimization of satellite modules. In this paper, we put forward a hybrid method called WL−LS by incorporating the WL sampling algorithm, the local search (LS) method, and the heuristic layout update strategies to deal with the layout optimization of a satellite module with performance constraints. Numerical results show that the proposed WL−LS algorithm is an effective method in designing the layout of satellite modules.

2 Problem statement

A simplified model of the international commercial communication satellite module (INTELSAT-III) (Teng et al., 2010) can be described as follows. The satellite module consists of several parts including a cylindrical container, a standing column, two bearing plates, and N objects which need to be packed. The column is located at the center of the module, perpendicular to the base surface of the satellite module (Fig. 1a). The objects can be fixed on either the top or bottom surfaces (P1−P4) of both plates (Fig. 1b). So, there are four surfaces available for locating these N objects, A i (i=1, 2, …, N). In this study, objects are considered to be homogeneous rigid bodies, and are simplified as cuboids and cylinders. The two circular bearing plates have the same thickness, mass, and radius R0. The objective of the layout optimization of the satellite module is to make the radius of the satellite module and the inertia moments of the whole system as small as possible. The final layout needs to meet the following constraints: (1) All the objects should be placed within the layout space and not overlap each other; (2) The error between the centroid of the whole system and the expected centroid should not exceed an allowable value, and be as small as possible; (3) The equilibrium degree error of the whole system should not exceed an allowable value, and be as small as possible.

Fig. 1
figure 1

Diagrammatic sketch of the simplified satellite layout: (a) 3D diagrammatic sketch of the layout; (b) 2D diagrammatic sketch of the layout; (c) layout of one object on a bearing plate

Full size image

In this paper, to describe clearly the mathematical model of the problem, we adopt four coordinate systems introduced by Sun and Teng (2003): Oxyz, O′x′y′z′, O″x″y″z″, and O‴x‴y‴z‴ (Fig. 2). In coordinate system Oxyz, the xOy plane coincides with the base surface of the satellite module, and the origin O is its geometric center. The upward direction of the z axis is from the satellite’s base to top. The x′O′y′ plane and the xOy plane are parallel, and the origin O′ of coordinate system O′x′y′z′ is the centroid of the whole satellite module. The z′ axis is parallel or coincides with the z axis, and they have the same direction. Coordinate system O″x″y″z″ is used to calculate the inertia moments of objects in relation to their own axes. The origin O″ is the centroid of the object, and the x″, y″, and z″ axes are the geometric symmetry axes of the object. The coordinate system O‴x‴y‴z‴ is used to calculate the inertia angles of the whole system. The origin O‴ coincides with O′, and the x‴, y‴, and z‴ axes have angles (θx′ θy′, θz′) with the x′, y′, and z′ axes. The counter-clockwise direction is positive.

Fig. 2
figure 2

Sketch map of coordinate systems: (a) coordinate system O″x″y″z″; (b) coordinate systems Oxyz, O′x′y′z′, and O‴x‴y‴z‴

Full size image

The mathematical model of the layout optimization problem discussed above can be described as follows:

$${\rm{min}}\;{f_1}\left( X \right),$$
((1))
$${\rm{min}}\;{f_2}\left( X \right) = {J_{x'}}\left( X \right) + {J_{y'}}\left( X \right) + {J_{z'}}\left( X \right)$$
((2))

subject to

  1. 1.

    non-overlapping constraints:

    $${g_1}\left( X \right) = \sum\limits_{i = 0}^{N - 1} {\sum\limits_{l = i + 1}^N {{\rm{int}}({A_i}) \cap {\rm{int}}({A_l}) = \emptyset ,} } $$
    ((3))
  2. 2.

    centroid error constraints:

    $${g_2}(X) = \left| {{x_{\rm{c}}} - {x_{\rm{e}}}} \right| \le 3.0,$$
    ((4))
    $${g_3}(X) = \left| {{y_{\rm{c}}} - {y_{\rm{e}}}} \right| \le 3.0,$$
    ((5))
    $${g_4}(X) = \left| {{z_{\rm{c}}} - {z_{\rm{e}}}} \right| \le 3.0,$$
    ((6))
  3. 3.

    equilibrium degree error constraints:

    $${g_5}(X) = \left| {{\theta _{x'}}(X)} \right| \le 0.03,$$
    ((7))
    $${g_6}(X) = |{\theta _{y'}}(X)| \le 0.03,$$
    ((8))
    $${g_7}(X) = |{\theta _{z'}}(X)| \le 0.03.$$
    ((9))

    Here X=(x1, y1, z1, α1, …, x i , y i , z i , α i , …, x N , y N , z N , α N ); α i are the orientation angles for the rectangles shown in Fig. 1c (for the circles, α i do not exist); N is the number of layout objects; (x i , y i , z i ) are the 3D coordinates of object A i ; f1(X) is the space utilization function, which is measured by radius R0 of the satellite module; f2(X) is the inertia moment of the whole system; Jx′(X), Jy′(X), and Jz′(X) are inertia moments of the whole system with respect to the coordinate system O′x′y′z′; int(A i ) denotes the internal part of object A i ; (xe, ye, ze) is the expected centroid of the whole system; (xc, yc, zc) is the real centroid position of the whole system; θx′(X), θy′(X), and θZ′(X) (Fig. 2) are the angles between the principal inertia axes of the module and the coordinate system O′x′y′z′ axes.

J x′ (X), J y′ (X), and J z′ (X) are calculated as follows:

$$\begin{array}{*{20}c} {{J_{x'}}(X) = } & {\sum\limits_{i = 0}^N {({J_{x''i}}{\rm{co}}{{\rm{s}}^2}{\alpha _i} + {J_{y''i}}{\rm{si}}{{\rm{n}}^2}{\alpha _i})} \quad \quad \quad \quad \quad \quad } \\ {} & { + \sum\limits_{i = 0}^N {{m_i}(y_i^2 + z_i^2) - (y_{\rm{c}}^2 + z_{\rm{c}}^2)\sum\limits_{i = 0}^N {{m_i}} ,} } \\ \end{array} $$
((10))
$$\begin{array}{*{20}c} {{J_{y'}}(X) = } & {\sum\limits_{i = 0}^N {({J_{y''i}}{\rm{co}}{{\rm{s}}^2}{\alpha _i} + {J_{x''i}}{\rm{si}}{{\rm{n}}^2}{\alpha _i})} \quad \quad \quad \quad \quad \quad } \\ {} & { + \sum\limits_{i = 0}^N {{m_i}(x_i^2 + z_i^2) - (x_c^2 + z_c^2)\sum\limits_{i = 0}^N {{m_i}} ,} } \\ \end{array} $$
((11))
$${J_{z'}}(X) = \sum\limits_{i = 0}^N {{J_{z''i}}\; + \sum\limits_{i = 0}^N {{m_i}(x_i^2 + y_i^2) - (x_{\rm{c}}^2 + y_{\rm{c}}^2)\sum\limits_{i = 0}^N {{m_i}} ,} } $$
((12))

where J x″i , J y″i , and J z″i are inertia moments of the ith object with respect to coordinate system O″x″y″z″, and m i is the mass of the ith object.

θx′(X), θy′(X), and θz′(X) are calculated as follows:

$${\theta _{x'}}(X) = {1 \over 2}{\rm{arctan}}\left( {{{2{J_{x'y'}}(X)} \over {{J_{x'}}(X) - {J_{y'}}(X)}}} \right),$$
((13))
$${\theta _{y'}}(X) = {1 \over 2}{\rm{arctan}}\left( {{{2{J_{x'z'}}(X)} \over {{J_{z'}}(X) - {J_{x'}}(X)}}} \right),$$
((14))
$${\theta _{z'}}(X) = {1 \over 2}{\rm{arctan}}\left( {{{2{J_{y'z'}}(X)} \over {{J_{z'}}(X) - {J_{y'}}(X)}}} \right),$$
((15))

where J x′y′ , J x′z′ , and J y′z′ are the products of inertia of the whole system with respect to coordinate system O′x′y′z′. Their calculations are as follows:

$$\begin{array}{*{20}c} {{J_{x'y'}}(X) = } & {\sum\limits_{i = 0}^N {\left[ {{m_i}{x_i}{y_i} + {{{J_{x''i}} + {m_i}(y_i^2 + z_i^2)} \over 2}{\rm{sin}}(2{\alpha _i})} \right]} \quad \quad } \\ {} & { - \sum\limits_{i = 0}^N {\left[ {{{{J_{y''i}} + {m_i}(x_i^2 + z_i^2)} \over 2}{\rm{sin}}(2{\alpha _i})} \right]} - {x_{\rm{c}}}{y_{\rm{c}}}\sum\limits_{i = 0}^N {{m_i}} ,} \\ \end{array} $$
((16))
$${J_{x'z'}}(X) = \sum\limits_{i = 0}^N {{m_i}{x_i}{z_i} - {x_{\rm{c}}}{z_{\rm{c}}}} \sum\limits_{i = 0}^N {{m_i}} ,$$
((17))
$${J_{y'z'}}(X) = \sum\limits_{i = 0}^N {{m_i}{y_i}{z_i} - {y_{\rm{c}}}{z_{\rm{c}}}} \sum\limits_{i = 0}^N {{m_i}} .$$
((18))

3 Framework for WL−LS

In this work, all tested objects were distributed onto the four bearing plate surfaces in advance according to the centripetal balancing method and human-computer cooperative method, as presented by Sun and Teng (2003) and Liu and Teng (2008). Therefore, in the process of layout optimization, every object always stays on the same surface; that is to say, the surface to which an object belongs is fixed. Suppose the layout X=(X1 X2, X3, X4)=(x11, y11, z11, α11, …, x1j, y1j, z1j, α1j; x21, y21, z21, α21, …, x2k, y2k, z2k, α2k; x31, y31, α31, …, x3m, y3m, z3m, α3m,; x41, y41, z41, α41, …, x4n, y4n, z4n, α4n), where X1, X2, X3, and X4 are the layouts of the four surfaces P1, P2, P3, and P4, respectively, and j, k, m, and n are the numbers of the objects distributed onto the four surfaces, respectively, j+k+m+n=N. Our goal is to obtain an optimal layout which minimizes the two objectives f1(X) and f2(X) and meanwhile satisfies the constraints given by Eqs. (3)–(9).

Suppose that all N objects, the satellite module, and its fixed components are smooth elastic solids. Since no object transfers between any two different surfaces are allowed, once one object is distributed onto a certain surface, it will never overlap another object on a different surface. So, we need to calculate simply the overlapping depth between each object and the fixed vessel, and the overlapping depth between any two different objects which are distributed onto the same surface. Using the quasi-physical strategy (Liu et al., 2010; 2011), the extrusive elastic potential energy of the whole system is as follows:

$$\begin{array}{*{20}c} {{f_3}(X) = \sum\limits_{i = 0}^{N - 1} {\sum\limits_{l = i + 1}^N {ud_{il}^2 = {f_3}({X_1}) + {f_3}({X_2}) + {f_3}({X_3}) + {f_3}({X_4})} } \quad \quad \quad \quad \quad \quad \quad } \\ { = \sum\limits_{i = 0}^{j - 1} {\sum\limits_{l = i + 1}^j {ud_{il}^2 + } \sum\limits_{i = 0}^{k - 1} {\sum\limits_{l = i + 1}^k {ud_{il}^2 + } } \sum\limits_{i = 0}^{m - 1} {\sum\limits_{l = i + 1}^m {ud_{il}^2 + } } \sum\limits_{i = 0}^{n - 1} {\sum\limits_{l = i + 1}^n {ud_{il}^2} } } ,} \\ \end{array} $$
((19))

where u is a physical coefficient. We set u=1 in this study. f3(X1), f3(X3), f3(X3), and f3(X4) denote the extrusive elastic potential energy which all objects allocated on four surfaces P1−P4 undergo, respectively. The calculation of d il is described in Section 4.

For a given radius R0 of the satellite module, we can convert the constrained optimization problem (2)–(9) into an unconstrained optimization problem by using the quasi-physical strategy and the penalty function method. The objective function is as follows:

$$\begin{array}{*{20}c} {{\rm{min}}\;E(X) = } & {{\omega _1}{f_2}(X) + {\omega _2}{f_3}(X)\quad \quad \quad \quad \quad } \\ {} & { + {\omega _3}[{g_2}(X) + {g_3}(X) + {g_4}(X)]\;\;} \\ {} & { + {\omega _4}[{g_5}(X) + {g_6}(X) + {g_7}(X)],} \\ \end{array} $$
((20))

where ω1, ω2, ω3, and ω4 are the penalty coefficients.

If a layout X satisfies constraints (3)–(9), then X is a feasible layout of the layout optimization problem for a satellite module. Therefore, if we find an effective algorithm which could solve the unconstrained optimization problem min(E(X)), we can obtain the optimal layout or approximate optimal layout of the original problem (1)–(9) by using effective search strategies, such as a dichotomous search, to obtain the smallest radii of the four surfaces. Suppose that the upper bound of the interval of the dichotomous search for a given surface is \(\overline {{R_i}} \) (i∈{1, 2, 3, 4}), which is initialized to the original radius R0 of the satellite module, and that the lower bound of the search interval is \(\underline {{R_i}} \) (i∈{1, 2, 3, 4}), which is initialized to the radius R of the standing column. The largest of the four radii obtained by the dichotomous search method is considered to be the smallest radius of the satellite module. The detailed procedure of the dichotomous search for a given surface P i (i∈{1, 2, 3, 4}) is as follows:

  1. (1)

    Set \(\overline {{R_i}} = {R_0},\;\;\underline {{R_i}} = R\).

  2. (2)

    Set \({r_i} = (\overline {{R_i}} + \underline {{R_i}} )/2\).

  3. (3)

    Run the WL−LS algorithm. If a feasible layout is obtained, set \(\overline {{R_i}} = {r_i}\); otherwise, set \(\underline {{R_i}} = {r_i}\).

  4. (4)

    If \(|\overline {{R_i}} - \underline {{R_i}} | > {10^{ - 4}}\), go to (2); otherwise, output \(\overline {{R_i}} \)as the smallest radius for this surface and exit.

4 Calculation of overlapping depth

To calculate the overlapping depth between two objects embedding each other, we can transfer the 3D model of the problem to a 2D one. Cylinders and cuboids are simplified as circles and rectangles, respectively. For a given bearing plate surface, a 2D Cartesian coordinate system Oxy is set up, which coincides with this surface, and its origin is located at the center of this surface. R0 is the radius of this surface. The set of rectangles to be located is RECT={Rect1, Rect2, …, Rect r }, where r is the number of rectangles distributed onto this surface; the set of circles to be located is CIR={Cir1, Cir2, …, Cir s }, where s is the number of circles distributed onto this surface. Cir indicates the fixed column. The ith rectangle is denoted by Rect i (p i , α i , a i , b i ), where p i =(x i , y i ) is the position of its centroid, α i ∈[0, π] is the orientation angle between the longer edge of the ith rectangle and the x axis in coordinate system Oxy, a, is the length of its longer edge, and b, is the length of its shorter edge. In this study, rectangles are placed orthogonally, i.e., α i =0 or π/2. The jth circle is denoted by Cir j (p j , r j ), where p j =(x j , y j ) is the position of its centroid and rj is its radius.

Various approaches have been proposed to compute the overlapping depth or the interference between objects, for example, the no-fit polygon (Bennell et al., 2001), the octree method (Wu et al., 1997), and the projection-separation approach (Li, 2010). In this study, we do not adhere to a single approach. Instead, according to the graphical characteristics of the two objects involved, we adopt different approaches that are most efficient. For example, if the two objects are rectangles, we adopt the projection approach; if the two objects are circles, we compare the distance between their centroids with the sum of the two radii; if the two objects are a rectangle and a circle, we adopt an approach based on a no-fit polygon. In the related previous work, discussions of approaches are quite general. Here, we give the specific calculation of overlapping depth which is suitable for orthogonally placed rectangles and circles.

4.1 Overlapping depth between a rectangle and a circle

Given a circle Cir j (p j , r j ) and a rectangle Rect i (p i , α i , a i , b i ), where α i =0 or π/2, and d ij is the overlapping depth between Rect i and Cir j , for every point p t on the circumference of Cir j , imagine there is a rectangle Rect t (p t , α i , a i , b i ), whose centroid is p t . Thus, there is a set of rectangles T={Rect t : t=1, 2, …}. Imagine every rectangle in T moves along the direction of vector V=(x t −x j , y t −y j ) until it just moves out of Cir j . After all rectangles in T move translationally in the plane, the new centroids of the rectangles in T form a continuous closed curve C (Fig. 3). C can be used to detect whether Rect i overlaps Cir j . If the centroid of Rect i falls inside the region enclosed by C, then Rect i overlaps Cir j , and vice versa.

Fig. 3
figure 3

A curve C for examining whether a rectangle overlaps a circle

TL (top line), LL (left line), RL (right line), and BL (bottom line) are lines parallel with the x or y axis. TLA (top-left arc), TRA (top-right arc), BLA (bottom-left arc), and BRA (bottom-right arc) are quadrants of a complete circle Cir j

Full size image

For simplicity, we assume p j =(0, 0), and α i =0. The equation for TL (top line in Fig. 3) is given by

$$y = {r_j} + {b_i}/2,\quad x \in [ - {a_i}/2,{a_i}/2].$$
((21))

The equation for TRA (top-right arc in Fig. 3) is given by

$$y = \sqrt {{r_j} - (x - {a_i}/2)} + {r_j},\quad x \in ({a_i}/2,{r_j} + {a_i}/2).$$
((22))

The equation for RL (right line in Fig. 3) is given by

$$x = {r_j} + {a_i}/2,\quad y \in [ - {b_i}/2,{b_i}/2].$$
((23))

Since C is symmetrical about the x and y axes, the equations for the other parts of C can be deduced easily. To calculate d ij , the region enclosed by C is divided into nine sub-regions (Fig. 4). If p i falls into S2, S4, S6, or S8, d ij equals the distance from p i to TL, LL, RL, or BL, respectively. For example, in Fig. 4, if p i falls into S6, the overlapping depth d ij equals the distance from p i to RL. If p i falls into S1, S 3 , S 7 , or S9, d ij equals r j minus the distance between p i and point a, b, c, or d, respectively. For example, in Fig. 4, if p i falls into S9, \({d_{ij}} = r_{j} - d_{dp_{i}}\), where \(d_{dp_{i}}\) is the distance between d and p i . If p i falls into S5, d ij is the shortest distance between p i and TL, LL, RL, and BL.

Fig. 4
figure 4

Sub-regions divided for calculating the overlapping depth between a circle and a rectangle

TL (top line), LL (left line), RL (right line), and BL (bottom line) are lines parallel with the x or y axis. TLA (top-left arc), TRA (top-right arc), BLA (bottom-left arc), and BRA (bottom-right arc) are quadrants of a complete circle Cir j . a, b, c, and d are the four vertices on the boundary of sub-region S5

Full size image

4.2 Overlapping depth between two rectangles

Given two rectangles Rect i (p i , α i , a i , b i ) and Rect j (p j , α j , a j , b j ), d ij denotes the overlapping depth between Recti and Rect j . l ix and l iy are the lengths of the edges of Rect i parallel with the x axis and the y axis, respectively, and l jx and l jy are the lengths of the edges of Rect j parallel with the x axis and y axis, respectively (Fig. 5). To examine whether two rectangles overlap, a straightforward method is to project rectangles Rec i and Rect j onto the x and y axes, and then to check whether the projecting lines overlap. The detection formulae are shown as follows:

$$|{x_i} - {x_j}| - ({l_{ix}} + {l_{jx}})/2 > 0,$$
((24))
$$|{y_i} - {y_j}| - ({l_{iy}} + {l_{jy}})/2 > 0.$$
((25))

If any two formulae hold, the two rectangles do not overlap, i.e., d ij =0. If neither of them holds, they overlap and the calculation of the overlapping depth is given by

$$\begin{array}{*{20}c} {{d_{ij}} = \left\{ {[({l_{ix}} + {l_{jx}})/2 - |{x_i} - {x_j}|} \right.{]^2}\quad \quad \quad \quad } \\ {{{\left. { + {{[({l_{iy}} + {l_{jy}})/2 - |{y_i} - {y_j}|]}^2}} \right\}}^{1/2}}.} \\ \end{array} $$
((26))
Fig. 5
figure 5

Two overlapping rectangles

Full size image

4.3 Overlapping depth between two circles

Given two circles Cir i (p i , r i ) and Cir j (p j , r j ), d ij is their overlapping depth (Fig. 6). The detection formula is given by

$${r_i} + {r_j} - \sqrt {{{({x_i} - {x_j})}^2} + {{({y_i} - {y_j})}^2}} > 0.$$
((27))

If expression (27) is true, circles Cir i and Cir j overlap, and the overlapping depth is calculated by

$${d_{ij}} = {r_i} + {r_j} - \sqrt {{{({x_i} - {x_j})}^2} + {{({y_i} - {y_j})}^2}} .$$
((28))
Fig. 6
figure 6

The overlapping depth of two circles

Full size image

4.4 Overlapping depth between a circle and a container

Given a circle Cir i (p i , r i ), if Cir i embeds into the standing column Cir, which is considered as a fixed circle, the calculation of the overlapping depth between Cir i and Cir is the same as in Section 4.3.

If circle Cir i (p i , r i ) overlaps the shell of the satellite module Cir0 (Fig. 7), i.e., \(\sqrt {x_i^2 + y_i^2} + {r_i} > {R_0}\), the calculation of overlapping depth d oi between Cir i and Cir0 is given by

$${d_{oi}} = \sqrt {x_i^2 + y_i^2} + {r_i} - {R_0}.$$
((29))
Fig. 7
figure 7

A circle embedding into the shell of the satellite module

Full size image

4.5 Overlapping depth between a rectangle and a container

Given a rectangle Rect j (p j , α j , a j , b j ), if Rect j embeds into the standing column Cir, which is considered as a fixed circle, the calculation of the overlapping depth between Rect j and Cir is the same as in Section 4.1.

Suppose α j =0. If rectangle Rect j (p j , α j , a j , b j ) overlaps the shell of satellite module Cir0 (Fig. 8), i.e., if \(\sqrt {{{(|{x_j}| + {a_j}/2)}^2} + {{(|{y_j}| + {b_j}/2)}^2}} > {R_0}\), the calculation of the overlapping depth d oj between Rect j and Cir0 is given by

$${d_{oj}} = \sqrt {{{(|{x_j}| + {a_j}/2)}^2} + {{(|{y_j}| + {b_j}/2)}^2}} - {R_0}.$$
((30))
Fig. 8
figure 8

A rectangle overlapping the satellite module shell

Full size image

5 Wang-Landau sampling method based on local search

5.1 Wang-Landau sampling method

The Wang-Landau (WL) sampling method is a novel Monte Carlo (MC) method introduced by Wang and Landau (2001). The energy levels of the models treated in the WL sampling method are discrete; however, the idea is very general and can be applied to any parameter. For a layout optimization problem where the energy is continuous, we must first bin the energy. Considering the value of energy in the range as a positive real number, we divide all possible energies in the energy landscape into finite intervals. For example, we divide [0, 5000] into 5000 individual intervals [0, 1), [1, 2), …, [4999, 5000), and numbers which are larger than 5000 are divided into a single energy interval [5000, 5000+). Thus, we obtain 5001 energy intervals. For simplicity, the energy interval [⎿E(X)⏌, ⎾E(X)⏋) is denoted by [E(X)], where ⎿E(X)⏌ rounds E(X) down to its nearest integer, and ⎾E(X)⏋ rounds E(X) up to its nearest integer. For example, E(X)=4.523 falls into energy interval [4, 5), which is denoted by [4.523]. Unlike conventional MC methods that directly generate a canonical distribution at a given temperature, this method is to estimate accurately the density of states g(E(X)) for the range of possible energies via a random walk which produces a flat histogram in the energy landscape, where E(X) is the energy of the whole system. The WL sampling method is based on the observation that if we perform a random walk in the energy landscape with a probability proportional to the reciprocal of the density of states 1/g(E(X), then a ‘flat’ histogram is generated for the energy distribution. By using a carefully controlled modification factor, the estimate for g(E(X)) is improved at each step of the random walk, which makes g(E(X)) converge to the correct value very quickly.

At the beginning of the WL sampling algorithm, all possible energies and the density of states function g(E(X)) are unknown. The density of states is set to be self-adaptive. If the random walk finds a new energy each time, we mark it as visited and set its density of states and the corresponding histogram to 1. In the simulations of the layout optimization for the satellite module, we begin the random walk in the energy landscape by heuristic layout update strategies (see Section 5.2), but the energy associated with each layout is accepted only with a probability proportional to the reciprocal of the density of states. Therefore, the acceptance probability from layout X1 to X2 is as follows: P(X1→X2)=min{exp(g(E(X1))−g(E(X2))), 1}. If X2 is accepted, then g(E(X2)) will be multiplied by a modification factor λ i , and its histogram H([E(X2)]) will be increased by one; that is to say, g(E(X2))=λ i *g(E(X2)), H([E(X2)])=H([E(X2)])+1 (i is initialized to 0). If X2 is not accepted, then g(E(X1))=λ i *g(E(X1)), H([E(X1)])=H([E(X1)])+1. If λ0 is too small, it will take a long time to find all the possible energies. On the contrary, if λ0 is too large, there will be statistical errors. In this study, we set λ0=1. The convergence of the WL sampling method is controlled by the flatness of the histogram. However, it is very difficult to obtain an absolutely flat histogram in practice. The so-called ‘flat histogram’ in the Wang-Landau sampling method means that all the entries of H([E(X)]) are not less than the histogram’s average <H([E(X)])> multiplied by k (0<k<1), where k is decided by the complexity of the system and the expected precision of g(E(X)). In this study, we set k=0.8, and check whether the histogram is flat every 103 MC sweep. When the histogram is flat, all the possible energies have been roughly visited an equal number of times, and the density of states converges to the true value with accuracy proportional to the modification factor lnλ i . Then we reduce the modification factor λ i to a finer one using a monotone decreasing function such as λi+1=0.5λ i , reset H([E(X)]) to 0 for all visited energy intervals [E(X)], and begin the next random walk. In this study, the modification factor and the corresponding decreasing function are selected by both experience and trial-and-error. From previous experience we set up candidate factors and their decreasing functions that might be suitable for the packing problem. Then from these candidates we select the best factor and its decreasing function through trial-and-error. When the modification factor λ i is less than a threshold λfinal, the algorithm is terminated, and g(E(X)) converges to its real value with precision. λfinal is the control parameter of g(E(X)) and determines the number of MC iterations in the whole simulation process. If λfinal is too small, the simulation will take a long time. On the contrary, if λfinal is too large, g(E(X)) will not converge to its real value. In this study, we set λfinal=0.00001.

Our goal is to find the layout with the lowest energy, so in the simulations we also keep the lowest energy Emin and the corresponding layout Xmin each time we find a new lower-energy layout.

5.2 Heuristic layout update strategies

An efficient layout update strategy is also impactful in the WL sampling simulations. According to the characteristics of the layout optimization of the satellite module, we propose the following heuristic layout update strategies:

Strategy 1 In each surface of current layout X, we pick out an object A j that has the largest relative extrusive elastic potential energy E j /S j to relocate, where \({E_j} = \sum\nolimits_{l = 0,l \ne j}^N {d_{ij}^2} \) is the extrusive elastic potential energy of the jth object exerted by other objects, and S j is the area of the jth object.

Strategy 2 Randomly generate 100 vacant points in each surface, where vacant points are the points that are inside this surface but not inside any object. The procedures are as follows:

  1. (1)

    Generate a random point inside the surface but outside the central column.

  2. (2)

    Judge whether this point falls into a certain object through computing the distance between this point and the centroid of the object. If this point does not fall into any object, save it as a vacant point and go to (3); otherwise, go to (1).

  3. (3)

    If 100 vacant points have been saved, exit; otherwise, go to (1).

Temporarily place the centroid of the chosen object at every vacant point. Here, if the chosen object is a rectangle, it is placed in two ways: one is with its long edge parallel with the x axis, and the other is with its long edge perpendicular to the x axis. Then, compute the extrusive elastic potential energy E j that associates to every vacant point and placing pattern. Finally, formally place the centroid of the chosen object A j at the vacant point, where the extrusive elastic potential energy of A j is the lowest. With the positions of other objects unchanged, we gain a new layout X′.

5.3 Local search

When a new layout X′ is obtained by heuristic layout update strategies, it may be very close to the global optimal layout, so any random layout update in the WL sampling algorithm may make search far from (and even farther and farther from) this global optimal layout. To avoid this, we adopt the local search (LS) method based on the gradient method (GM), which is a quasi-physical algorithm (Huang and Kang, 2004; He et al., 2013), to search for an optimal layout near X′.

GM is also known as the steepest descent method. The search direction is the negative gradient direction. We adopt an adaptive step size in this study. If the energy increases after one iteration, which indicates that the step size for this step is too large, we decrease h to 0.8h. If the energy decreases, we keep h. The procedures of the GM with an adaptive step size are described as follows:

  1. (1)

    Set h=1, hmin=10−4, ε=10−20.

  2. (2)

    Under the current layout X′, compute the gradient vector ∇(E(X′)) of E(X′) for each object A i (i=1, 2, …, N) in the x and y directions, and let X2=X′−h*∇(F(X′)).

  3. (3)

    If E(X2)>E(X′), set h=h*0.8.

  4. (4)

    Let X′=X2, X2=X′−h*∇(E(X′)).

  5. (5)

    If f3(X2)<ε or h<hmin, return X2; otherwise, go to (3).

5.4 Description of WL−LS

By combining the WL sampling method with heuristic layout update strategies and the LS procedure, we propose a hybrid WL−LS algorithm for the layout optimization of the satellite module. The calculation procedure of the WL−LS algorithm is outlined as follows:

  1. (1)

    Randomly produce an initial layout X1 based on four bearing plate surfaces. Set Xmin=X1, Emin=E(X1). Let the set of intervals containing visited energies be S={[E(X1)]}. Set the density of states function as g(E(X1))=1, and the histogram function as H([E(X1)])=1. Set i=0, l=0, λ0=1, k=0.8.

  2. (2)

    In each surface P i (i=1, 2, 3, 4) of the current layout X1, pick object A j (j=1, 2, 3, 4) with the largest E j /S j , where E j is the extrusive elastic potential energy of the jth object. Copy the current layout X1.

  3. (3)

    Relocate each picked object A j by using heuristic layout update strategy 2 in surface P i (i=1, 2, 3, 4), and gain a new layout X′.

  4. (4)

    Call the GM procedure. The outcome layout of GM is denoted by X2. Compute E(X2). Set g(E(X2))=1, H([E(X2)])=1, l=l+1.

  5. (5)

    If [E(X2)]∉S, let S=S◡{[E(X2)]}.

  6. (6)

    If random(0, 1)<min{exp[g(E(X1))−g(E(X2))], 1}, then accept X2 (i.e., let X1=X2, E(X1)=E(X2)); otherwise, do not accept X2.

  7. (7)

    Update the density of states g(E(X)) and the histogram H([E(X)]). That is, if X2 is accepted, let g(E(X2))=g(E(X2))*λ i , H([E(X2)])=H([E(X2)])+1, and go to (8); otherwise, let g(E(X1))=g(E(X1))*λ i , H([E(X1)])=H([E(X1)])+1, and go to (9).

  8. (8)

    If E(X2)<Emin, let Xmin=X2, and Emin=E(X2).

  9. (9)

    If l%1000=0, go to (10); otherwise, go to (2).

  10. (10)

    If H([E(X)])≥k<H([E(X)])> for all visited energy intervals [E(X)]∈S, then go to (11); otherwise, go to (2).

  11. (11)

    Set λi+1=λ i *0.5, i=i+1.

  12. (12)

    If λ i <0.00001, then output Emin and Xmin, and terminate the iteration; otherwise, reset H([E(X)])=0 and keep g(E(X)) for all visited energies E(X), and go to (2).

6 Experimental results and discussion

To test the computational performance of the WL−LS algorithm, we applied it in two instances. Both instances are based on the international commercial communication satellite module INTELSAT-III with different technological parameters, one with 51 objects and the other with 53 objects. We implemented the WL−LS algorithm in the Java programming language and ran it on a PC with 1.5 GHz CPU and 2.0 GB RAM. For each instance, the WL−LS algorithm was run 50 times independently to optimize the layout of objects on the four bearing plate surfaces. According to these 50 results, we can gain a pareto optimal set for each instance. In this study, we specify the solution with the smallest enveloping radius as the optimal solution for preference.

For the first instance, we suppose that 51 objects have been allocated on the four bearing surfaces in the satellite module. The objects’ dimensions and masses have been given by Liu and Teng (2008). The first 31 objects are cylinders, and the remaining 20 objects are cuboids. The parameters of the satellite module are as follows: the radius of the bearing plate is R0=500 mm; the radius of the column is R=100 mm; the heights from the base of the satellite to the bottoms of the lower plate and the upper plate are H1=500 mm and H2=1050 mm, respectively; the height of the column is H3=1400 mm; the thickness of each plate is Ht=20 mm; the mass of the satellite module (including the shell, two bearing plates, and the column) is 349.557 kg, and its centroid is (0, 0, 859) mm in the coordinate system Oxyz, with an inertia matrix

(31)

The technological requirements for the final layout of the whole system are given by Eqs. (3)–(9), with parameters xe=0, ye=0, ze=780 mm. The penalty coefficients are set as follows: ω1=10−1, ω2=106, ω3=104, and ω4=104.

For this instance, the pareto optimal set by WL−LS is P(Q, M)={(460.61, 690.22), (459.28, 691.70), (458.95, 692.65), (457.76, 693.07)}, where Q denotes the enveloping radius and M the inertia moment of the system. We choose the solution with Q=457.76, M=693.07 as the optimal solution for preference. The best and average computational results by WL−LS for instance 1 are shown in Tables 1 and 2, respectively, in comparison with those by the HAKD method (Liu and Teng, 2008). From Tables 1 and 2, we can see that the proposed WL−LS algorithm outperforms the HAKD method in every aspect. Compared with the best results by the HAKD method, the inertia moment of the best layout obtained by WL−LS decreases by (711.55−693.07)/711.55×100%=2.60%; the enveloping radius decreases by 0.35%; the centroid position error decreases by 29.12%; the inertia angle error reduces by 76.5%. Table 3 shows the data of the optimal layout obtained by WL−LS. Fig. 9 shows the diagram of the optimal layout obtained by WL−LS.

Fig. 9
figure 9

The 2D diagram of the optimal layout obtained by WL−LS for instance 1: (a) surface 1; (b) surface 2; (c) surface 3; (d) surface 4

Full size image
Table 1 Comparison of the best results by HAKD, WL, WL+GM, WL+HS, and WL−LS for instance 1
Full size table
Table 2 Comparison of the average results by HAKD, WL, WL+GM, WL+HS, and WL−LS for instance 1
Full size table
Table 3 Data for optimal layout by WL−LS for instance 1
Full size table

To further understand the effects of each component in WL−LS, we ran WL alone, WL with local search GM and without heuristic layout update strategies (denoted as WL+GM), and WL with heuristic layout update strategies (HS) and without local search GM (denoted as WL+HS) 50 times. The best and average results for 50 independent runs by WL, WL+GM, and WL+HS are listed in Tables 1 and 2, respectively, for comparison. The pareto optimal set by WL is P(Q, M)={(468.58, 696.78), (467.24, 697.40), (466.06, 699.08), (465.05, 710.96)}, and we choose the solution with Q=465.05, M=710.96 as the optimal solution for preference. The pareto optimal set by WL+GM is P(Q, M)={(462.51, 693.88), (462.03, 694.79), (460.63, 695.16), (459.17, 695.34)}, and we choose the solution with Q=459.17, M=695.34 as the optimal solution for preference. The pareto optimal set by WL+HS is P(Q, M)={(465.09, 705.86), (464.00, 706.16), (463.01, 707.47), (462.48, 708.31)}, and we choose the solution with Q=462.48, M=708.31 as the optimal solution for preference. Note that in each run of WL and WL+GM, we randomly update the layout. In WL−LS, the WL method is used mainly to execute a global search. The heuristic layout update strategies are used to generate new layouts, and the local search procedure based on the gradient method is used to search for lower-energy layouts near newly generated layouts. Through the comparison of computational results, one can see that WL+GM and WL+HS improve the performance of WL alone in different aspects. However, by combining the merits of GM and HS, WL−LS significantly outperforms the other three algorithms WL, WL+GM, and WL+HS in both the best and average results.

For the second instance, we suppose that 53 objects have also been allocated on the four bearing surfaces in the satellite module. The objects’ dimensions and masses have been given by Sun and Teng (2003). The first 24 objects are cuboids, and the remaining 29 objects are cylinders. The parameters of the satellite module are as follows: the radius of the bearing plate is R0=500 mm; the radius of the column is R=100 mm; the heights from the base of the satellite to the bottoms of the lower plate and the upper plate are H1=300 mm and H2=830 mm, respectively; the height of the column is H3=1150 mm; the thickness of both plates is Ht=20 mm; the mass of the satellite module (including the shell, two bearing plates, and the column) is 776.53 kg, and its centroid is (0, 0, 553.56) mm in coordinate system Oxyz, with an inertia matrix

(32)

The technological requirements for the final layout of the whole system are given by Eqs. (3)–(9) with parameters xe=0, ye=0, ze=523.26 mm. The penalty coefficients are set as follows: ω1=10−1, ω2=106, ω3=104, and ω4=104.

For this instance, the pareto optimal set by WL−LS is P(Q, M)={(479.31, 791.67), (478.14, 792.64), (476.07, 794.33), (475.11, 795.59)}, where Q denotes the enveloping radius and M the inertia moment of the system. We choose the solution with Q=475.11, M=795.59 as the optimal solution for preference. The best computational results by WL−LS for instance 2 are shown in Table 4, in comparison with those by the hybrid knowledge fusion (HKF) method (Wang and Teng, 2009). From Table 4, we can see that the proposed WL−LS algorithm outperforms the HKF method in every aspect. Compared with the best results by the HKF method, the inertia moment of the best layout obtained by WL−LS decreases by (796.15−795.59)/796.15×100%=0.07%; the centroid position error decreases by 99.70%; the inertia angle error decreases by 67.81%; the enveloping radius decreases by 4.98%. A diagram of the optimal layout obtained by WL−LS is shown in Fig. 10. Table 5 shows the data for the optimal layout obtained by WL−LS.

Fig. 10
figure 10

Optimal layout obtained by WL−LS for instance 2: (a) surface 1; (b) surface 2; (c) surface 3; (d) surface 4

Full size image
Table 4 Comparison of the best results by HKF, WL, WL+GM, WL+HS, and WL−LS for instance 2
Full size table
Table 5 Data for optimal layout by WL−LS for instance 2
Full size table

We also compare the best and average results for 50 independent runs by WL, WL+GM, and WL+HS in Tables 4 and 6, respectively. The pareto optimal set by WL is P(Q, M)={(491.40, 799.71), (489.43, 800.77), (488.30, 801.34), (487.05, 802.13)}, and we choose the solution with Q=487.05, M=802.13 as the optimal solution for preference. The pareto optimal set by WL+GM is P(Q, M)={(480.74, 795.27), (479.69, 796.67), (479.05, 796.96), (478.50, 797.01)}, and we choose the solution with Q=478.50, M=797.01 as the optimal solution for preference. The pareto optimal set by WL+HS is P(Q, M)={(486.41, 797.21), (485.58, 798.54), (484.32, 799.79), (483.49, 800.69)}, and we choose the solution with Q=483.49, M=800.69 as the optimal solution for preference. Through the comparison of computational results, we can draw the same conclusion as that for instance 1.

Table 6 Comparison of the average results by WL, WL+GM, WL+HS, and WL−LS for instance 2
Full size table

7 Conclusions and future work

When optimizing the layout design of a satellite module, it is easy for an algorithm to get trapped in local minima separated by high-energy barriers. To address this problem, we use a hybrid WL−LS method, which incorporates the LS procedure based on gradient descent and heuristic layout update strategies into the Wang-Landau sampling method. To improve the efficiency of WL−LS, we adopt an accurate and fast method for computing the overlapping depth between two objects (such as two rectangular objects, two circular objects, or a rectangular object and a circular object) embedding each other. Numerical results show that WL−LS outperforms methods in the literature. There are also several problems that need to be solved in the future: (1) In this study, the rectangles are placed orthogonally, which restricts the flexibility of the layout design. We need to find the methods for computing the overlapping depth with respect to arbitrarily placed rectangles. (2) The distribution of objects onto different surfaces is not discussed in this paper, and we simply adopt the allocation scheme proposed in the literature. In the future, we will focus research on these aspects.