Computational design synthesis for Fabrication-Aware assembly problems using building objects with dimensional variations

State-of-the-art computational design approaches for designing the assembly of a structure often rely on high-precision building objects, which can be difficult and costly to fabricate. To consider a wider range of construction materials and fabrication processes, e.g. low-cost 3D printing and brick making, this research introduces a novel computational design method for fabrication-aware assembly problems using existing building components with dimensional variations. Given a finite set of building objects, this study develops a statistical preprocessing approach and a modified pattern search with an efficient integer linear programming method to generate the assembly of a target structure, such that it is statically stable and optimized for geometric smoothness and space utilization. In two case studies, the proposed algorithm exhibits significantly better boundary fitting and space utilization compared with a genetic algorithm and a standard pattern search method and lower computational cost compared with a genetic algorithm.


Introduction
In many scenarios, especially for developing countries and for small to medium scale civilian projects, a low-cost, rapid, and reliable design and fabrication process is often prioritized over accuracy. However, the inaccuracy in fabrication is usually compensated for by overproducing the building objects and only using the ones with satisfactory characteristics, which creates waste. As an example, building objects printed with low-cost 3D printers using sustainable and low-cost materials, e.g. industrial waste, can have low precision of dimensions and unpredictable variations in material properties. Further innovations in the fabrication process may set higher constraints on material properties and site conditions such as requiring skilled personnel, high-tech and highprecision equipment, and more resources, e.g. electricity. This research investigates a new computational design approach that is fabrication-aware and accommodates inaccuracies in the fabrication process during the design optimization to best use the building components that have been produced with low-precision fabrication processes.
The problem presented is similar computationally to packing and stacking problems. The main goal of a packing problem is to pack one or several empty domains with objects of different dimensions such that the summation of the objects' metrics (e.g. volume, prices, weights, profits, etc.) is optimized or the total number of occupied domains is minimized. Packing problems and their variants have been widely studied in many disciplines like scheduling [1], land planning [2], cutting stock and packing of industrial products [3,4], vehicle loading [5], memory allocation, and several other logistics and robotics-related problems [6]. While the models of packing problems address the full utilization of empty spaces or domains, and can be applied to various disciplines, stacking problems tend to focus on the structural behavior of a stack and thus emphasize the physical and structural interactions between two adjacent objects as well as the stability of the resulting wall assembly [7]. However, most existing research on stacking problems focus on the technical potential in simulation and modeling of irregular objects (e.g. natural stones) and stack them at the optimized contact points with robots [7][8][9]. This is usually inefficient and impractical in certain industries where the effort to simulate and model the irregular objects significantly exceeds the effort to physically preprocess and regulate the raw material. For example, it can be more demanding in terms of computational power and equipment to scan and model natural stones, and then simulate the optimal stacking point than simply cutting planar surfaces out of stones (e.g., cuboidal objects) and building simple models for positioning. Even if intensive computation and optimization are available, it is still difficult to stack more than five natural stones vertically without reinforcement [7], and the stability of the stack is sensitive to different loadings. Further, low-cost fabrication processes, e. g., stone-cutting and some low-cost additive manufacturing processes [19], can produce something in between irregular objects and highprecision objects. This paper focuses on low-precision objects that have planar surfaces (i.e. cuboidal objects). The detailed explanation is presented in Section 2.1.
In the Architectural, Engineering and Civil industry (AEC industry), a few designers and researchers have applied the techniques of packing and stacking problems to create assemblies of structures using discrete components, e.g. masonry structures such as vaults and brick walls [10][11][12][13][14] with predefined structural performing requirements and constraints. The advances in structural analysis and optimization methods enable designers to go beyond their limitations to explore alternative structural designs through automating design tasks in computational design synthesis [15]. Most existing assembly problems obtain the optimized results either by varying the sizes of the building objects and yielding the dimensions as output or by choosing from a predefined inventory of objects, while both are assuming the fabrication process produces bricks with actual dimensions similar or same as designed [16]. When the imprecise physical outcomes underperform the target functionality, most existing research investigates more accurate fabrication processes. As mentioned before, high-precision fabrication can not only cost unnecessary time and resources but may also set limitations on material choices to achieve target resolution and precision. Thus, there exists a research gap on how to use the dimensionally imperfect building objects produced with low-precision processes to assemble a stable stack with smooth boundaries and better space utilization.
Therefore, the research problem presented in this paper is formulated as follows: given a limited number of unique cuboidal building objects, whose measurable dimensions are known and fixed after fabrication, find the optimized subset of the available building objects, with their location and rotation angle to construct an assembly that best meets the defined objective. The design objective is to best fit the target geometry while satisfying predefined constraints, such as static stability [17,18], space utilization, and geometric smoothness. Fig. 1 illustrates the example case for a curved wall assembly with cuboidal building objects subject to dimensional variation. To distinguish this application from the packing and stacking problem, this paper refers to it as an assembly problem. This paper investigates the dimensional variation of 25 cuboid blocks printed by a powder bed 3D printer with geopolymers [19] and uses this information as the input to the computational design method. The details are described in Section 3.1.
This paper begins by specifying the building objects of interest and introducing the state-of-the-art methods in Section 2. Then in Section 3, a preprocessing procedure, including a data modeling method and a binning strategy, is introduced to accelerate the search process for certain building objects, followed by the proposed computational design  method using modified pattern search with an efficient integer linear programming sorting subproblem. The proposed method is applied to two example cases and benchmarked with a genetic algorithm and a standard pattern search without the sorting subproblem. The results are presented in Section 4 and discussed in Section 5. Section 6 includes the conclusion, the discussion of the limits of the proposed method, and future directions.

Background
This section introduces the basis for the method and related work, including the building objects of interest, uncertainty-based design optimization (UDO), existing discrete optimization methods for assembly problems, and the adopted representation model.

Building objects of interest
In this work, we aim to create a computational design synthesis method that can accommodate dimensional variations arising from fabrication processes to balance the cost of accuracy of fabrication with the cost of computation (Fig. 2). As mentioned in Section 1, assembly with in-situ, unprocessed building objects (e.g. natural stones) [7][8][9] is often computationally prohibitive. While the cost of fabrication is low, it takes a tremendous amount of computational resources to scan, model, and simulate the unprocessed building objects as well as find the right positions and angles to maintain the stability of the assembly. Highprecision components, on the contrary, demand accurate fabrication processes, which are costly and energy-intensive, and may limit the choices of usable materials. Hence, this research focuses on building objects manufactured with low to moderate precision. Such processes are typically lower cost in both computation and fabrication, compared with the other two scenarios, and have other potential benefits, e.g. using locally sourced material. Thus, uncertainties from the design, fabrication and assembly environment, that can result in fluctuations in the building object dimensions are accommodated without creating waste. It is noteworthy that the low-cost processes may also lead to poor material quality. However, the stability of a stack relies primarily on the geometry of the structure rather than on material strength [20]. Further, parts with poor quality should not be used if they fail the mechanical property test. Thus, material properties are not considered in this research.

Uncertainty-Based design optimization (UDO)
Uncertainty-Based Design Optimization (UDO) is an ongoing research topic to tackle uncertainties in design variables, parameters, and optimization systems, which can be traced back to the 1950 s [21,22]. Most existing literature on this topic focuses on the design and optimization algorithms under uncertainties that are quantified and propagated [23][24][25][26][27]. Successful applications have been observed in the civil engineering industry [28][29][30]. To generate designs that are robust to all the uncertain variables (e.g. dimensions, loads, material strength, etc), the conventional UDO scheme requires a large dataset for uncertainty quantification and a series of computationally expensive evaluations or analyses (e.g. finite element analysis) for uncertainty propagation of each uncertain variable. Since the computational cost grows exponentially with the dimensionality of the UDO, this research is motivated to approach the UDO problem in a different direction and take uncertain variables (i.e. dimensions) out of the propagation process by directly measuring the exact dimensions after fabrication and using them as input and constraints for the computational design method. It is more practical and efficient to measure the parameters that remain unchanged after the early phases (e.g. dimensional variations) than other uncertain parameters that occur after the design phase (e.g. loads) or parameters whose exact values cannot be measured without destroying the specific object (e.g. strength). Since the dimensions of a cuboidal block are easy to measure, low-tech measuring approaches, instead of high-tech methods such as laser scanning [31], can be adopted.

Existing discrete optimization methods for cuboid building objects
Since the building objects are treated as unique items, discrete optimization methods are considered. The literature on similar deterministic problems, such as the space packing problem and automated brick assembly construction problem [32], focus on packing or assembling using the limited types of given objects. The optimization methods include a split-and-merge-based greedy algorithm [33,34], cellular automata [35,36], shape annealing [37], simulated annealing [38,39], beam search [36], pattern/local search [35,38,40], and evolutionary algorithm [41,42]. However, these existing algorithms employ a limited set of building objects with replacements rather than dealing with unique objects. For example, the building objects are chosen from a predefined set of object types with identical copies, or duplicates, available for each type. The copies have identical dimensions and should not deviate from the design. Unique objects refer to building objects that have independent and random dimensional variations. Fig. 3 shows the differences between these two sets of building objects. Considering the building objects as unique objects drastically increases the complexity of the problem by introducing a large set of constraints that makes the design variables exclusive and strongly interdependent, i.e. if one object has been selected and placed into the current assembly, then it can not be used elsewhere in the same assembly, and thus the aforementioned algorithms become intractable. Moreover, stability is not addressed in these works as significantly as in the problem discussed in this study. This is because the structural analysis is not substantial for lego-like assemblies for their interlocking features can tightly connect two bricks. Therefore, this research fills in the gap by handling unique building objects and addressing structural stability, i.e. rigid body force equilibrium, of the whole assembly.

Representations
Typically, two types of representations are considered in the existing research: voxel and object. Voxelization, i.e. the discretization of the 3D space into a mesh of voxels, is suitable for space arrangement, especially for irregular containers or for space with irregular constraints, and performs well for unit-based objects with discrete dimensional variations and are identical within the same group, e.g. the LEGO® bricks. This representation method has been widely applied in the existing methods mentioned in Section 2.3. Applying voxelization in this paper would constrain the size of a voxel to be considerably small in order to better represent the continuously varying dimensions, leading to high demand for computational power. The second option, namely the object representation, represents an object using its geometric characteristics, which works well for components with simple geometries that could be explicitly defined by several metrics (e.g. length, width, radius), by the relative positions of vertices (e.g. polyhedron), or by formulas (e.g. sine-cosine surfaces).
As discussed in Section 2.1, this research focuses on objects produced with low-precision fabrication methods, which often have larger dimensional variation. The object representation adopted for this research is building objects represented geometrically with three dimensions (length, depth, height) and located with four degrees of freedom in the space (x,y,z coordinates and rotations along the z-axis). Considering the scale of the problem (e.g. a masonry wall consisting of hundreds of blocks), scanning the objects and representing them using point cloud [43] is computationally expensive. The cuboidal representation is sufficient as long as the surfaces are planar and pair-wise orthogonal/parallel, which can be achieved with a low-precision fabrication method.

Method
In this section, a novel fabrication-aware assembly method for building objects with unique dimensions is presented. To increase computational efficiency, a preprocessing procedure is introduced, including a data modeling method and a novel statistical binning strategy. Then, this paper employs computational design optimization to find designs that best fit the target geometry while satisfying certain predefined requirements for stability, space utilization, and geometric smoothness. For search methods, this research develops a modified pattern search with an efficient integer linear programming sorting algorithm, which is compared with two benchmarking algorithms, including a standard pattern search without sorting and a genetic algorithm. Fig. 4 shows the flowchart of the proposed method.

Preprocessing
This paper addresses the dimensional variations of the cuboidal building objects. However, searching through all of the building objects during each iteration in the process is inefficient and may lead to computational intractability. Before starting the optimization, a simple preprocessing procedure is introduced to 1) quantify the data with a statistical model and 2) reduce computational complexity by statistically clustering the dimensional variations. The example case of 25 binder-jet printed geopolymer-based bricks is presented in this section to demonstrate this preprocessing method.

Data quantification
For each brick, its length, depth, and height are measured according to ASTM C67/C67M and are shown in the appendix in Table A. 1. The mean and standard deviation, as well as the maximum and minimum values, are shown in Table 1.
The building objects of a large-scale structure, (e.g. bricks for a masonry wall), are usually produced in batches with low-precision fabrication that results in dimensional variations. Sorting the building objects in advance eases the search for an object with certain dimensions. One promising method that has been widely applied in the manufacturing industry is binning strategies, or data binning [44]. However, binning strategies generated directly from the raw data may lose generality and cause an overfitting issue. To obtain a robust binning strategy, the dimensional data measured in Section 3.1.1 is modeled with a distribution model. This research integrates expert knowledge into the quantification by fitting the height data to a statistical distribution model. According to the probabilistic model code edited by the Joint Committee on Structural Safety [45], a Gaussian distribution is recommended for dimensions and measurement errors. After statistical validation, this paper assumes that the height data follows a normal distribution with a mean of 45.402 and a standard deviation of 3.178 (N(45.402, 3.178 2 )). The same method applies to the other two dimensions.

Statistical binning strategies
Searching all the available building objects whenever a certain type of object is needed in the design synthesis is often unnecessary and can lead to computational intractability. To accelerate computation, this research takes in the distribution model generated from the raw data and establishes a binning strategy. The existing binning strategy is to slice the dimensional space into bins with equal spans, which works well for uniform distribution. However, most fabrication variations in dimensions are not uniform but rather follow a probabilistic model (e.g. a normal or log-normal distribution [45]), depending on the fabrication process. Another issue unaddressed in literature is how to choose the proper span of the bins, such that trade-offs between the reduction of computational cost and the representativeness of the building objects in the bins are balanced. Thus, this research introduces the following method to minimize the number of bins while preserving 1) the uniqueness of each building object in the bin as required, and 2) the practical amount of the building objects in each bin: where A is the matrix controlling the bin width and P is the matrix specifying the practical amount of building objects in each bin. Given a distribution model of a single parameter (e.g. a normal distribution model of the height of the building objects), the first step is to discretize the parameter space, from the lower bound h lb to the upper bound h ub , into n small equally spanned intervals, as shown in Fig. 5. Then the boundaries at each interval can be either activated or deactivated, serving as potential boundaries of the bins. A boolean variable x i is assigned to the right boundary of the i th interval, representing the activation of this specific boundary. For example, if all boundaries are switched on (i.e. x i = 1, ∀i ∈ {0, ⋯, n}), then this binning strategy generates n +1 bins, including two extreme bins, ( − ∞, h lb ) and [h ub , + ∞), that are not accounted for during the discretization. The objective function is to minimize the sum of the boolean variables, which is equivalent to minimizing the number of bins.
To preserve the representativeness of the building objects' uniqueness, the user is expected to set a maximum span for each bin diff r , defined by the maximum number of small intervals contained in each bin b max . Then the first constraint Ax ≥ 1 can be easily composed. For instance, if b max = 3, the corresponding formulation is shown in Eq. (2). This linear system indicates that for every three successive boundaries, at least one boundary is activated as a binning boundary so that the maximum variation allowed in one bin is three times the span of the small intervals. ⎡ ⎢ ⎢ ⎣ Since the wall is layer-based, this study encourages each bin to  contain more building objects required by this layer, such that mixed bins in one layer can be reduced. Therefore, the second constraint Px ≥ p r aims to ensure a sufficient number of building objects clustered into each bin, which is equivalent to constraining the probability of falling into each bin to be no smaller than p r . This constraint takes two input parameters from the user: the approximated number of building objects in each layer (n col ) and the approximated number of all building objects in this wall assembly (n all ). The minimum probability p r is defined as p r = n col /n all . Each row in Px ≥ p r shall be in the form of: where x i and x j are the left and right boundaries of a bin, meaning they are activated (x i = 1,x j = 1). p i and p j are the probability values at the heights corresponding to x i and x j , calculated by the cumulative distribution function of the distribution model. However, Eq.(3) has conditions in the boolean variables and cannot be integrated into an integer programming formulation. This paper thus reformulates Eq.(3) as follows.
Equation (4) is generalized for any x i and x j since it is equivalent to Eq. (3) when both x i and x j are activated (x i = 1,x j = 1), and it remains valid when any of x i and x j is not activated. One exception exists at the positive extreme but can be easily solved by setting up h ub properly such that the corresponding p n is no smaller than 1 − p r /2. To compose the matrix P, the linear system only needs to consider x i and x j within b min intervals since b min is given. Namely, for b min = 3, Eq. (4) can be expressed as.
For the distribution model in section 3.1.1 (N(45.402, 3.178 2 )), this study implements the above binning method and obtains nine bins shown in Fig. 6. The required minimum probability of each bin p r is set to 0.075, resulting in h lb = 40.83 and h ub = 49.98, which forms the 85% confidence interval. To initialize the design space, the range [40.83, 49.98] is discretized into 50 evenly spanned intervals with a userspecified maximum span limit for each bin diff r = 1.5. The 25 data points from Section 3.1.1 are plotted as blue dots along the x-axis. This binning strategy prepares the data for the fabrication-aware computational design synthesis method.
Since the proposed binning strategy does not alter the features of each object but rather pre-sort and bin them, the major variance caused by the proposed binning strategy is how much computational cost can be saved, which depends on the accuracy of the model. The accuracy of the model can be improved by sampling more data (i.e. dimensions of different objects with the same nominal design and from the same fabrication process) or acquiring a priori information from literature or experts.

Computational design synthesis using modified pattern search and a sorting method
This paper formulates an assembly problem given a limited number of cuboidal building objects, whose dimensions are known values after fabrication. Unlike most uncertainty-based optimization problems that integrate the dimensional uncertainty into the uncertainty propagation step, this research considers these building objects as unique items and clusters them into a limited number of bins, as demonstrated in Section 3.1.2. Next, a pattern search method and a sorting strategy are developed to solve the discrete optimization problem. In this section, three major modules of the proposed method are presented in detail: 1) Generation: generating the initial design randomly and updating the current design based on the guidance yielded from the previous iteration; 2) Evaluation: evaluating the performance of the current design on the predefined objective and constraints; 3) Guidance: generating and updating strategies in the search directions and step sizes according to the improvement in the performance of the current designs, passing search information to the next iteration. An efficient integer linear programming method is developed in Section 3.2.4, which significantly speeds up the approach and improves its performance, as demonstrated in Section 4.

Generation
This paper focuses on layer-based assembly problems, i.e. the base surfaces of the building objects in the same layer are placed at the same height. The height variations in one layer are smoothed to form the base for the next layer. In practice, this could be done by using infill materials such as earth, straw, etc., which do not necessarily have adhesive features like mortar. The void space is also minimized in the proposed method, as formulated in objective 2 in Section 3.2.2. As described in Section 3.1, each object in the material pool (Ω B ) is defined with three dimensions (length, depth, height) and is assigned to a bin during preprocessing. In this section, if a building object is chosen to be used in the structure, its location is defined by four degrees of freedom in the space (x,y,z coordinates and the rotational angle along the vertical/z-axis). Due to the characteristics of cuboids, rotations are simplified and constrained to two options: 0 and 90 degrees.
Similar to one of the common strategies used in combinatorial optimization, the generation of an assembled wall assembly is conducted with a certain number of alternating operations including in-bin swappings and out-of-bin swappings, which are specified by the modified pattern search in Section 3.2.3. Rotations are excluded from these options since rotations are integrated into an efficient integer programming method presented in Section 3.2.4.
In the design synthesis, the system first generates an initial wall assembly by randomly assigning one bin to each layer, which is called the default bin of this layer. Then, a sequenced set of building objects is sampled from the assigned bin and assembled sequentially to form the target geometry.
The current design is then sent to the evaluation module (Section 3.2.2) for performance evaluation, which serves as the basis for the guidance module to generate and update strategies in the search directions and step sizes accordingly, as described in detail in Section 3.2.3. With the new guidance strategy, the system loops back to this generation module, where it updates the current design accordingly. For instance, if the guidance requires s out out-of-bin swappings and s in in-bin swappings, the system will: 1) Select s out layers using the weighted random sampling method [46] based on the performance of each layer. Then replace their default bins with other bins. The building objects in these modified layers are updated accordingly by replacing these objects with unused building objects randomly sampled from their corresponding updated bins. If the number of unused building objects from a bin is not sufficient, switch to its neighboring bins without changing the default bin of this layer. 2) Randomly choose s in objects in the wall assembly and swap them with unused objects in the default bins of the layers they reside in.
An example operation with two out-of-bin swappings and one in-bin swapping is illustrated in Fig. 7. After applying the swap operations, the updated design of the wall assembly is ready for evaluation.

Evaluation
A newly generated design is given to the Evaluation module, where it will be evaluated and the objective function calculated according to the performance. This section introduces the requirements, constraints, and objective criteria of the assembly problem of interest and presents the overall formulation employed by the proposed method.

1) Constraint: Stability.
The stability of a wall assembly is essential since adhesives between each pair of adjacent objects are often unavailable, e.g. packing, or unfavorable, e.g. mortar between bricks is a carbon-intensive material. The stability is evaluated by modeling the building objects as rigid elements, computing force resultants on interfaces, and checking if tensile forces are required on these interfaces to maintain stability [17,18]. The algorithm from Whiting et al [17], with modifications, is applied and formulates the stability problem as a linear program, as shown in the equation below.
where B eq is an equilibrium matrix that records the affiliated vertices of the contacting interfaces (with directional information) of each building object; f i n is the normal force resultant at vertex i that is not restricted in sign when defined by a difference of two non-negative variables f i− n and f i+ n [47]; w is a vector composed of external forces (e.g. self-weight) of each building object; The first constraint in Eq. (6) is the equilibrium constraint that projects and sums up f i n for each axis as well as their introduced moments for each axis, such that both the force equilibrium Fig. 7. An example operation with two out-of-bin swappings and one in-bin swapping. and the moment equilibrium are maintained for each building object subject to external forces (e.g. self-weight) and moments, if any.
This formulation penalizes the tensile forces. If and only if the objective value C stab is zero, will this structure, without inter-object adhesives, be statically stable. To minimize the use of gap-filling material among adjacent building objects, the constraint is regarded as a hard constraint, provided that a feasible solution exists.

2) Objective 1: Target Smoothness.
Maintaining the target smoothness of the surface of an assembly is one of the key arguments for the proposed method to work with building objects that have dimensional variations. This paper considers the smoothness of both the top and side surfaces. The target smoothness T smoo is defined as the L-p norm of the differences between the actual surface and the target surface, as shown in Eq. (7).
where Ω b denotes the set of top and/or side surfaces. ∑ ‖x b − t xb ‖ p is the summation of the L-p norm of the vector of actual boundary lengths x b and its target value t xb . Based on user preferences, the properties of the smoothness, t xb can be defined as a combination or any of the following targets: 1) an upper bound, 2) a guiding boundary that does not restrict the signs of the difference, which are illustrated with the case studies in Section 4, or 3) the mean value of x b . The mean value representation does not push the design to fit the target but rather is an indicator for the variance and thus pushes the design to be as smooth as possible. In this paper, L-2 norm (i.e. the Euclidean distance) is selected to measure the variation.
To maximize the space utilization (e.g. for packing problems) or to minimize the usage of gap-filling material, the second objective minimizes void spaces formed by the neighboring objects within the whole assembly, which is represented as the difference between the bounding space and the total space taken by the building objects, as shown in Eq. (8): where Ω dim is a set of dimensions in the space (e.g. x, y, z),  The performance of the current design is calculated as a linear combination of all three objectives and a penalty function for stability when the stability indicator C stab is larger than zero, which is shown in Eq. (9): where w smoo , w void , and w stab are user-defined weights for target smoothness, void spaces, and stability, respectively. A large value for w stab is recommended since the stability is regarded as a hard constraint. If and only if the stability constraint is satisfied, max(C stab , 0) will be zero, imposing no penalty to the overall objective.
The objective function is formulated based on different practical needs. For example, for partition walls, smooth-line boundaries are needed to fit, e.g. door and window frames and the room ceiling. From straight walls in a conventional masonry building to curved walls in masonry towers and corbelled vaults, stability and smoothed boundaries are also needed so that the walls can meet closely with fewer gaps. The objective to minimize the void space relates to minimizing the infill material in practice such that the construction process can be simplified, lead time reduced and is more sustainable. The consideration of the dimensional variation is to lift the constraints on mid-to-high-precision building objects and to reduce material waste.

Guidance
If the current design evaluated in Section 3.2.2 is not converged, feedback on the search direction is given back to the generation modules in Section 3.2.1, where the current design is updated and re-optimized. This search direction is generated in the guidance module.
For the combinatorial problem addressed in this paper, it is difficult to define the "search direction" in the variable space, thus a gradientbased search method is not applicable. One of the widely-applied search methods for combinatorial optimization is local search by sequentially swapping the worst performed elements. Thus, instead of defining directions in the design space, this work manipulates the search space where the axes are the different alternating operations, e.g. the number of in-bin swappings and the number of out-of-bin swappings. Moreover, this work embeds these operations into a gradient-free search method modified from the Hooke and Jeeves pattern search method [48]. Pseudocode for the proposed algorithm, together with the aforementioned two modules, is shown in Fig. 8. Line 3, 8, 14, and 23 are the proposed sorting subproblems utilizing an efficient integer linear programming formulation, which will be introduced in detail in the  The search direction d, from an old design to a new design, is composed of several elements, as shown in Eq. (10): (10).where each p i represents the number of its corresponding operations i conducted to reach the new design from the old design. Ω OP is the set of available operations (e.g. out-of-bin swapping and in-bin swapping). The L-2 norm, or the Euclidean norm, of the search direction d is defined by the step size δ step : ‖d‖ 2 = δ step , which gets reduced each time after a successful move, i.e. a move that leads to improvement in the overall objective value.
For each temporary point (TP), at most n neigh neighbors at distance δ step are searched. If none of these neighbors are better than the previous base point (BP), reduce the step size by a factor r reduce , and search again.
This reduction iterates until the termination condition is satisfied (δ step < 1). With the number of consecutive step reductions n s being monitored, reheating can be triggered when n s reaches the predefined maximum limit n smax , which increases the step size by a factor r increase . However, since the reheating process can only be triggered once, the process will eventually terminate when a better BP is found or when δ step is reduced to smaller than one. It should be noted that δ step < 1 is the lower-bound for stopping criteria because δ step is the number of operations to be conducted, thus any δ step that is smaller than one has no practical meaning. The factor r reduce is in the range (0,1), which is also adopted in other classic optimization algorithms, such as simulated annealing [39]. The exact value is case-dependent and should be determined together with the initial step sizes and the computational cost for each iteration.

Efficient integer linear programming Sub-problem
Section 3.2.1 states that the rotating operation is omitted on purpose. This section explains the reason and integrates the rotating operation in an efficient integer linear programming subproblem at the end of the generation module. The corresponding objective is to minimize the gap between a target length and the actual length of this layer. The motivation of this formulation is to reduce redundant and computationally expensive operations by sorting the objects of each layer (i.e. switching their directions) rather than randomly rotating the objects from the worst-performing layer. The options for the rotation angle of an object, when it is placed on a layer, have been restricted to either 0 degrees or 90 degrees. So, either the edge in the length direction or the edge in the depth direction is contributing to the total length of the layer. For cuboidal building objects, this constraint actively reduces the void space since every two adjacent objects in the same layer are aligned.
Since each layer has a default bin that is not necessarily the same as the other layers, this sorting technique can be applied to each layer independently. Given a constant target t L bi for layer i, the sorting subproblem is formulated in Eq. (11).
where p ij , q ij are boolean variables representing the activation of placement in the length or depth direction of a building object j in layer i. So p i,j , q i,j are mutually exclusive, i.e. only one can be chosen, as shown by the second constraint. B i is the set of building objects in layer i, which belongs to the material pool Ω B . L i indicates the actual length calculated by accumulating the individual length, which is either the length or the depth of the building object based on its placement. This algorithm rearranges the directions of certain building objects in a layer and is implemented as a sorting algorithm for each row when a new design is generated by the method described in Section 3.2.1. It is computationally more efficient since it only takes each row of the structure as the subproblem rather than taking the whole structure as a mixed-integer linear programming problem, which is intractable.

Results
This section presents results for two case studies, including a curved wall assembly and a hollowed cylinder wall assembly. It also compares the proposed modified pattern search method with the integer linear programming subproblem (PS + S) with alternative algorithms including the genetic algorithm (GA) and a standard pattern search (PS) [48] without the sorting method. For a fair comparison and a feasible implementation, this work modifies the simple genetic algorithm [49], which allows for the implementation of the proposed operation scheme, including in-bin swapping, out-of-bin swapping, and rotation. The same objective function defined in Eq. (9) is used for all three algorithms. The hyperparameters of the modified genetic algorithm are set as follows: the population size is set to 60, the mutation rate is set to 0.25, the parent ratio is set to 0.3, the crossover rate is set to 0.9, and the maximum number of generations is set to 100. The selections of the population size, the mutation rate, and the parent ratio is based on the parameter study, with the result listed in Table A.3. The maximum number of generations is selected such that the convergence is guaranteed for both case studies, while the crossover rate is set to 0.9 considering the relatively small population size.
Both the curved wall and the hollowed cylinder wall are layer-based assemblies composed of a subset of objects from a material pool sampled  with the height distribution model developed in Section 3.1.2. The material pool contains 120 objects that are clustered into the bins generated in Section 3.1.2, as shown in Fig. 9. Since the geopolymerbased bricks in Section 3.1.1 have limited variations in length and depth, this paper demonstrates the efficiency and robustness of the proposed integer linear programming sorting method by introducing larger variations in length and depth (standard deviations increased to 10.937 and 10.984 for length and depth, respectively) while maintaining the same distribution model for the height of the objects. Specific dimensions of each generated object are attached in Table A The proposed method is implemented in Python 3.6 together with a CVXPY 1.0 solver for integer linear programming. Both of the optimization example cases are run on a standard laptop using a CPU with eight cores (i7-8550U) and 16 GB of RAM. All numerical results are average values obtained by running each algorithm 25 times.

Problem formulation for 9x9 curved wall assembly
The first case study is a curved wall assembly with 9x9 objects (9 layers with 9 objects in each layer). The target geometry, as defined by Eq. (12), is a sine-cosine surface that guides the placement of a given building object in the out-of-plane direction.
The goal of this case is to use a subset of objects in the pool to assemble a statically stable curved wall assembly, according to the target geometry, such that both the variations of the bonding surfaces and the void spaces within the wall are minimized. The general problem formulation follows that given in Section 3.2, which is not repeated in this section. The target smoothness objective sets a specific target value for the current boundary and a smoothness requirement for the current top surface, as shown in the following equation: where is the variation of the top surface x top , ‖x rb − t xrb ‖ 2 is the Euclidean distance between the current boundary of the surface on the right x rb and its specified target value t xrb . In this case study, t xrb denotes a vector of identical values defined by the expected number of objects in one layer times the average length of all objects. Without loss of generality, this example case smooths one of the side surfaces (i.e. the left side) by aligning the first object in every layer, so only x rb is involved in Eq. (13). The initial step size δ step is set to 15, with reduction rate r reduction set to 0.8. For reheating, the maximum limit for consecutive reductions n smax is set to 10 and the increase rate r increase is set to 3.

Results for 9x9 curved wall assembly
The performance of the proposed pattern search with sorting algorithm (PS + S) for the curved wall assembly, is illustrated in Fig. 10. From the first iteration to the last iteration, the proposed algorithm reduces the average objective value by 99.9% (from 3,235,704 to 3394), with void space V void reduced by 37.4% (from 5100 to 3194), Z-variation Var ( x top ) reduced by 49.2% (from 0.427 to 0.217), X-variation ‖x rb − t xrb ‖ 2 reduced by 99.8% (from 5.0174 to 0.0119). Fig. 11 shows the curved wall design from a sampled set of iterations, with the target right boundary represented by the grey panels in the 3D plots and the red vertical lines in the 2D plots below.
The performance comparison among the proposed approach (PS + S) and the alternative methods, including a genetic algorithm (GA) and a standard pattern search (PS), is shown in Table 2. For each algorithm,  the designs of the first iteration and the last iteration (optimized result) are illustrated in Fig. 12.
The data (mean values) given and illustrated above show that, compared with the other two algorithms, the proposed algorithm significantly smooths the right boundary (X-variation) by at least 99.9%, reduces the top variation (Z-variation) by half, and minimizes the void spaces by at least 64%. This algorithm also has a smaller variance for all of the objective values. In addition, the proposed algorithm requires about one-third of the computational cost compared with the GA.
In summary, when compared with a GA, the proposed PS + S method significantly reduces the computational cost and improves the performance, and compared with the conventional pattern search, it drastically improves the performance (i.e. lower objective values) with similar computational cost. It also exhibits more robust performance than the other two algorithms.

Problem formulation for 4x20 hollowed cylinder wall assembly
The second illustrative case study is a hollowed cylinder wall assembly with 4x20 building objects (4 layers with 20 building objects in each layer). The goal of this case study is to use a subset of a predefined library of building objects to assemble a statically stable cylinder wall assembly of radius R, defined by the users, such that both the variations of the bonding surfaces and the void spaces within the wall are minimized. This example sets the value of R with the following equation: where n B is the number of building objects expected for each layer, and μ l is the average length of the building objects. In this case, R is calculated using the data from Section 3.1 and is rounded to 176. Since this case study uses cuboidal building objects to configure a cylinder structure, the inner surface of the target cylinder wall is used as a reference surface for placing the building objects. When a new building object is to be placed, its interior surface (surface of the cuboidal objects facing the inner space of the cylinder) is attached to the previous building object on this layer. Unlike the curved wall assembly, the cylinder wall assembly requires its building objects to rotate along the center of the cylinder. The rotation angle of the newly-placed building object is calculated such that the perpendicular bisector of the interior surface goes through the center of the cylinder. The general problem formulation follows that shown in Section 3.2, which is not repeated in this section. The initial step size δ step is set to 25 while the other parameters remain the same from the curved wall assembly case. The sorting approach using the integer linear programming formulation can still be implemented, with minor modifications due to the existence of a constraint for the overlapping between the starting/ left boundary surface and the ending/right boundary surface. Namely, the specific target value for the surface on the right will be regarded as an upper bound limit, which requires, for each layer i its accumulated length L i must be smaller than the target length t L bi . Therefore, to remove the absolute operation in the objective function in Eq. (11) for the integer linear programming subproblem, Eq. (11) is modified to Eq. (15) in this case, where one additional constraint t x rb − L i ≥ 0 is added.

Results for 4x20 hollowed cylinder wall assembly
The performance of the proposed pattern search with sorting algorithm (PS + S) for the hollowed cylinder wall assembly is illustrated in Fig. 13. The plots exhibit a convergence towards the final result, except for some perturbation in the performance on void space. From the first iteration to the last, the proposed algorithm reduces the average objective value by more than 99.99% (from 560,315,394 to 5506), with void space V void reduced by 17.8% (from 6299 to 5176), Z-variation Var ( x top ) reduced by 47.9% (from 0.6340 to 0.3303), and X-variation ‖x rb − t xrb ‖ 2 reduced by almost 100 % (from 1.7852 to 2.77×10 -8 ). Fig. 14 shows the hollowed cylinder wall design and the corresponding interior surfaces from a sampled set of iterations.
The performance comparison among the proposed approach (PS + S) Fig. 15. Comparison of the optimized results of different methods for the cylinder wall assembly. and its alternative methodologies is shown in Table 3. For each algorithm, the designs of the first iteration and the last iteration (optimized result) are illustrated in Fig. 15. The data (mean values) shows that the proposed algorithm significantly smooths the right boundary by at least 99.9%, reduces the top variation (Z-variation) by half, and minimizes the void spaces by at least 35%. Besides, the proposed algorithm only needs approximately one fourth of the computational cost of the GA. This proposed algorithm yields a smaller variance for not only all the objective values but also the total number of objective calculations.

Discussion
This study develops a computational design synthesis for fabricationaware assembly problems in the context of predefined, unique components without replacement and subject to random, independent, and continuous dimensional variations. Two representative case studies are presented, including a curved wall assembly and a hollowed cylinder wall assembly, where the proposed algorithm is compared with a genetic algorithm, and a standard pattern search.
The optimized results obtained using the proposed method exhibit higher efficiency when compared with the initial design in the first iteration and when benchmarked with other alternative algorithms. From the first iteration to the last, both the curved wall assembly and the hollowed cylinder wall assembly have objective values reduced by at least 99.8%: void space V void reduced by at least 17.8%, Z-variation Var ( x top ) reduced by half, and X-variation ‖x rb − t xrb ‖ 2 reduced by more than 99.99% The relatively less significant improvement in several objectives (e.g. 17.8% improvement in minimizing the void space for the hollowed cylinder wall assembly) is due to its corresponding small weight on the void space, which can be adjusted based on user preferences. In the first two to three iterations of the proposed algorithm, increases in the average values of the void space, the Z-variation, and the X-variation can be observed in Fig. 11 and Fig. 13. This is because the first iteration of each run can start with an unstable assembly, leading to an extremely high objective value. To output a stable assembly, the proposed algorithm tends to compromise other objectives in the next few iterations. Variances of all of the objectives also decrease along with the iterations and convergence can be observed in Fig. 11 and Fig. 13.
When benchmarked with GA and PS, the proposed algorithm (PS + S) significantly smooths the right boundary by at least 99.9%, smooths the top boundary by half, and minimizes the void spaces by at least 35%. The drastic improvement in the X-variation is ascribable to the implementation of the efficient integer linear programming method for the sorting subproblem. As shown in Fig. 12 and Fig. 15, the initial design from the first iteration of each algorithm tends to be unstable, resulting in high objective values. Though all the algorithms can stabilize the assemblies at the final iteration, only the proposed PS + S yields better and optimized results without compromising all other objectives. Without the sorting algorithm, the GA and PS return designs with apparent gaps and roughness on the boundaries. In terms of computational cost, the proposed algorithm (PS + S) has approximately at most half of the computational cost of the GA in terms of objective function calculations and is similar to the computational cost of the PS. The mean time in seconds of the PS + S is higher than that of the PS, due to the sorting algorithm. However, the significant improvement in results merits the increase. Further, the proposed PS + S method also possesses remarkable robustness, i. e. significantly lower standard deviation in the performance metrics.
It is worth noting that although these two case studies consist of roughly the same numbers of building objects, the computational time varies significantly. The reasons are: 1) there is a minor difference in the objective formulation of case 2, as described in Section 4.2.1. Since exceeding the target bounds leads to overlapping in the cylinder wall, it is more challenging for the algorithm to find a feasible solution for this case; and 2) the optimization problem formulation is layer-based both in terms of the objective functions, i.e. the target smoothness, and in terms of the operations in the proposed method, i.e. each layer is assigned with a default bin and then in-bin or out-of-bin swappings are conducted. Therefore, with more variables in each row, the problem becomes larger, and thus needs more runtime even with a similar number of objective evaluations.
The method presented works best for cuboidal objects or cuboidal envelopes for certain irregularly configured objects. Thus, to preserve the efficiency of the proposed method, it is recommended that low-cost post-processing is conducted, as discussed in Section 2.1 and Section 2.4, such that the objects can be modeled as cuboidal. Future work considers extensions to other shapes, material quality variations, new low-cost fabrication processes, structural loads, and objective function formulations. It is worth mentioning that the novel statistical preprocessing approach proposed to cluster the building objects works not only for normal distributions but also for any distribution model derived from the a priori information related to the specific fabrication process considered. The only modification needed is to replace the cumulative distribution function with that of the distribution model considered.

Conclusion
This paper introduces a novel computational design synthesis method to efficiently solve an assembly problem with dimensional variations resulting from less precise fabrication processes or the construction environment. It addresses the unvisited gap between conventional methods for found building objects with no fabrication processes and precise building objects fabricated by high-tech processes that include excessive control on the tolerance. Given a set of fabricated building objects of given dimensions, this research develops a binning strategy and a modified pattern search with an efficient integer linear programming sorting algorithm to solve the assembly problem. The method tackles the combinatorial problem by embedding the operations into a gradient-free Hooke and Jeeves pattern search method. The two case studies show that the proposed method yields significantly improved results, which reduces the objective value by at least 99.8% compared with the initial design and by at least 88.2% compared with the final design generated by a genetic algorithm, while requiring less computation. The method remains compatible with a wide range of applications in different industries and disciplines, for example, the space packing problem, e.g. 3D knapsack problem, and stacking problem.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.    As shown in Table A.3, changing the parameters of the genetic algorithm influence the performance of the algorithm differently: 1) Population size (combination 1, 2, 3): reducing the population size (from 60 to 30) can drastically reduce the computation time as well as the number of objective calculations, which leads to an early convergence with large objective values. On the other hand, when increasing the population size from 60 to 120, the objective values and the number of objective calculations increase as well. However, the computation time is reduced by 20%, which is not in proportion to the increase in objective evaluations. This is because, with a larger population size, the GA converges with a smaller number of generations, thus less computation is needed. 2) Mutation rate (combination 1, 4, 5): increasing the mutation rate from 0.25 to 0.5 leads to large objective values while decreasing it from 0.25 to 0.125, though still increasing the objective values, has less influence on the performance. Both changes result in less computation time. 3) Parent ratio (combination 1, 6, 7): changing the parent ratio significantly influences the objective values, where the large objective values, resulting from reducing the parent ratio from 0.3 to 0.15, are due to early convergence. It can be concluded that, though changing the parameters does influence the final results and computational cost, the performance metrics are still within the same magnitude. Thus, combination 1, with the best objective values, is selected for the GA implemented for the case studies in this research.