Minimizing form errors in additive manufacturing with part build orientation: An optimization method for continuous solution spaces


 For additive manufacturing (AM) to be successfully implemented in manufacturing systems, the geometric accuracy of components must be controlled in terms of form, fit, and function. Because the accuracy of AM products is greatly affected by the part build orientation, this factor dictates the achievable tolerances and thereby the ability to incorporate AM technologies in a large-scale production. This article describes a novel optimization method for minimizing form errors based on the geometric features of the part. The described method enables the combination of separate expressions for each feature to create a continuous solution space. Consequently, the optimal part build orientation can be precisely determined based on a mathematical description of the effect of build direction on each surface type. The proposed method is demonstrated in two case studies by step-by-step descriptions including discussions on viability and possible extensions. The results indicate good performance and enable flexible prioritization and trade-offs between tolerance characteristics.

Abstract: For additive manufacturing (AM) to be successfully implemented in manufacturing systems, the geometric accuracy of components must be controlled in terms of form, fit, and function. Because the accuracy of AM products is greatly affected by the part build orientation, this factor dictates the achievable tolerances and thereby the ability to incorporate AM technologies in a large-scale production. This article describes a novel optimization method for minimizing form errors based on the geometric features of the part. The described method enables the combination of separate expressions for each feature to create a continuous solution space. Consequently, the optimal part build orientation can be precisely determined based on a mathematical description of the effect of build direction on each surface type. The proposed method is demonstrated in two case studies by step-by-step descriptions including discussions on viability and possible extensions. The results indicate good performance and enable flexible prioritization and trade-offs between tolerance characteristics.
Keywords: additive manufacturing, part build orientation, quality, accuracy, build direction, optimization

Introduction
Additive manufacturing (AM) holds the potential to revolutionize the manufacturing industry through topologyoptimized lightweight structures and mass-customized designs. However, the full potential of the technology remains largely unexploited due to cost restrictions and quality issues. While AM is utilized in medical, aerospace, and automotive industries for small production volumes, this is only possible because every component is carefully engineered and validated through an iterative process before production is initiated. For mass customization to truly reach its potential in agile manufacturing, automated process planning for AM must be developed to ensure quality and consistency in the production. While AM offers design freedom to create innovative free form surfaces previously unattainable by conventional manufacturing methods, traditional shape features will still be present in novel designs, such as the interfaces between components of an assembly. The interconnection between AM and conventional manufacturing technologies in manufacturing systems warrants the control of geometric accuracy in terms of traditional tolerancing features.
To enable consistent production of unique components while meeting quality requirements, it is necessary to develop valid models and methods for predicting, mitigating, and adapting to variation in the build process. One of the major determining factors for geometric accuracy in AM is the part build orientation, i.e., the direction in which the material is added to the substrate to realize the geometry [1]. This article describes a method for the precise determination of optimal part build orientation to minimize the deviations from nominal to the actual geometry by considering the geometric features of the part.
The accuracy of an additively manufactured surface partly depends on its curvature (or lack thereof). Therefore, the proposed method enables separate mathematical models to be applied to each surface type. These models are then populated with data from the CAD model and combined into a single expression of quality as a function of part build orientation. The result is a continuous objective function for the entire solution space, which enables the identification of optimal part build orientations or feasible regions for secondary objectives.
The remainder of this article is structured as follows: First, the theoretic foundations and previous work are outlined in Section 2 before the method is described in Section 3. In Section 4, two case studies are presented to demonstrate the method step by step before a brief discussion on the viability and possible extensions is presented in Section 5. Finally, Section 6 summarizes this article and provides directions for the future work.

Theoretic background
In general, AM techniques successively add layers of material to create an object [1]. This layered manner of fabrication is a decisive factor in how accuracy errors occur in the AM process. Dantan et al. [2] identified a range of defect modes in fused filament fabrication, many of which can be extended to other AM technologies. Defect modes relating to the direction of material deposition, as well as errors in the machine/tool movement, influence the accuracy of the produced surface regardless of technology and material. Apart from the various surface defects and inaccuracies present in most AM processes, the products also exhibit anisotropic mechanical properties where the behavior depends on the build direction [3]. Consequently, the part build orientation is a decisive factor in the final part quality in terms of both accuracy and mechanical properties.
The staircase effect is perhaps the most illustrative example of how part build orientation is vital in AM. The layered manufacturing approach inevitably leaves a characteristic pattern on inclined surfaces. This pattern emerges when the contours of two subsequent layers cannot align perfectly. The result is a stepped surface as displayed in Figure 1, commonly referred to as the staircase effect.
The staircase effect can be modeled in two dimensions as the cusp height h cusp , i.e., the shortest distance from the inner corner of a step, to the hypotenuse of the right triangle of Figure 1 [4]. The cusp height may be calculated from the layer thickness h layer and surface angle θ, which is the angle between the build direction [ ] = ẑ 0 0 1 and the surface normal vector → n as described by Alexander et al. [5]: Another measure of the staircase effect is the volumetric error (VE) introduced by Masood et al. [6], which corresponds to the volumetric difference between the CAD model and the realized geometry. This solution soon becomes quite complex when extended from two to three dimensions. Figure 1 illustrates the deviations due to the staircase effect where the area below the dashed line is lost. The VE is the result of integrating this area along the perimeter of the layer. Complex surfaces drastically complicate the computation of VE, and simplifications are often necessary.
The obvious solution for mitigating the staircase effect is to make sure all surfaces are oriented either parallel or perpendicular to the build direction. However, when a part consists of several features in different directions (which is the case for most functional components), there will be no part build orientation in which all surfaces are either parallel or perpendicular to the build direction. Furthermore, curved surfaces do not benefit from perpendicular orientations as this will maximize the staircase effect. Therefore, a trade-off must be made to converge on a globally optimal solution, where the staircase effect is minimized. The orientation problem for accuracy is however composed of many more failure modes other than the staircase effect, all of which matters in modeling the final geometry.
While the part build orientation should be considered in the design stage, this is not always possible. Consequently, the information available when deciding the part build orientation may differ from full CAD model with tolerances, to the tessellated STereoLithography file format (STL) commonly used in the final stages of process planning for AM. Naturally, methods for dealing with the orientation problem with these varying knowledge levels have been proposed in the literaturesome relying on computational methods and others on a human expert. In addition, some may focus on finding a minimal solution fast, while others may prefer to spend more time to find the optimal orientation. The many possible combinations of methods, objectives, and limitations have resulted in a plethora of approaches ranging from the general to the highly specific. However, the crux of the orientation problem remains to achieve the best possible result with the information at hand, within a reasonable time.

Related work
Part build orientation is a key factor in AM as it is easily manipulated and has a clear influence on final properties. The effect of part build orientation has been the subject of many research efforts but remains an open issue [7]. Methods for determining the optimal part build orientation can be largely divided into two groups: (i) those evaluating a set of candidate solutions with respect to an objective function and (ii) those mathematically describing the solution space to explore this continuous space for the optimal solution. A selection of the most relevant related methods is showcased in Table 1. The interested reader may refer to refs [7,8] for recent reviews on the orientation problem.
The first group is the most common and starts by identifying candidate orientations. This process is either based on a set of rules (e.g., flat surfaces of the convex hull) or by discretization of the solution space, i.e., certain intervals of rotation about one or more predefined axes. In the next step, either an exhaustive search is performed where all candidates are considered or a guided search is conducted, e.g., by a genetic algorithm (GA). These discrete methods have certain advantages, but they inherently fail to consider the entire (continuous) solution space and therefore risk missing good orientations. In addition, any attempt to refine the search space by including more candidate orientations inevitably increases the computation time.
The second group, on the other hand, grants access to the entire solution space, where the complexity of this solution space generally correlates with the complexity of the geometry. This category is not as well explored, perhaps due to the simplicity of discrete solution spaces or the ability of discrete approaches to handle discontinuous functions. Nevertheless, a continuous solution space may facilitate more nuanced objective functions and complex solution spaces. Continuous solution spaces can be explored by evolutionary algorithms (EAs), or gradientbased methods where knowledge of the topology of the solution space is exploited in an iterative search for the global optimum.
The approach described herein constitutes a hybrid of the two methods outlined earlier: The solution space is described with a differentiable function, which is used to identify critical points. The critical points will then constitute the candidate orientations in the final evaluation, which identifies the global optimum. The novelty of this approach lies in the combination of simplicity from generalizing surface types, and the flexibility brought Table 1: Characteristics of related work in chronological order
forward by the general framework that can be populated with any objective function. Cheng et al. [9] considered all flat surfaces as candidates for determining the optimal part build orientation and evaluated every candidate orientation with respect to the accuracy, build time, and stability. A similar approach was proposed for the minimization of cost by Alexander et al. [5] who included the cost of postprocessing and build time in the cost calculations. Zhang et al. [16] generated a set of candidate orientations from every shape feature and evaluated them for several attributes including surface roughness and support volume. Similar methods have also been proposed with facet clustering, where groups of facets are considered collectively based on how they are affected by part build orientation [20][21]. All these methods benefit from the generalization of how different surfaces are affected by the build direction. While various methods are employed in these studies in the search for an optimal solution, none of them consider a continuous solution space generated from the geometry.
Zhang and Li [13] proposed a discretization of the solution space (i.e., a unit sphere), and let each facet of the STL file promote the two orientations parallel to the facet normal, as well as the great circle corresponding to the perpendicular of the normal vector. The authors describe methods to make the selection using an exhaustive search for minimizing VE, as well as a GA for combined optimization of VE and part height [13,14]. The authors argue for the use of GA over exhaustive search when the number of discretized points becomes large due to time concerns. Discretization allows for controlling the resolution of the solution space and therefore also the computational cost. However, the approach remains oblivious to any effects other than those of the predefined points in the solution space.
Das et al. [15,17] proposed a method for minimizing form errors and support structures by formulating a minimization problem to be solved in MATLAB with gradient-based optimization. The method is based on the onedimensional tolerance maps of Paul and Anand [23,24] for cylindricity and Arni and Gupta [25] for flatnessboth theoretically derived from the staircase effect. A similar approach was proposed by Chowdhury et al. [19] who continued to compensate the geometry for any expected deviations still present after optimization. However, both of these proposed methods involve a risk of convergence to local optima due to the gradient-based approach. The convergence to local optima is avoided in the study by Budinoff and McMains [18] who performs an exhaustive search of one-degree increments to identify feasible regions from which the final orientation may be selected.
The existing body of literature describes various approaches to the orientation problem with various benefits and drawbacks. Exploiting higher level information about local topology facilitates the generation of candidate orientations. This information can also be used to construct continuous solution spaces for achievable tolerances and other objectives. The method presented herein benefits from the generalization of part geometry to reduce the number of parameters in optimization, while also accessing the entire solution space for mathematical analysis.

Mathematical foundations
There is no shortage of mathematical formulations of the orientation problem in the academic literature. The plethora of formulations arises from the subtle differences in AM technologies, which have varying parameters with different effects. The number of formulations is further amplified by the deviating scope and objectives of previous works. For instance, the VE is a common measure of accuracy in AM. However, the calculation of VE is based on fundamental assumptions regarding the surface profile, typically assuming right-angled steps. However, the real surface will be filleted as various effects will round off the corners and hence throw off the theoretical models [26,27].
Nevertheless, the effect of part build orientation on surface accuracy is indisputable in the current AM systems, although of less concern when the layers are thinner [28]. Therefore, the modeling of quality as a function of part build orientation is warranted, and mathematical descriptions of orientations are necessary.
According to Euler's rotation theorem, any orientation of a rigid body in 3 can be described as a sequence of three basic rotations ( ) α β γ , , , where no two subsequent rotations are performed about the same axis. The basic rotations are performed about the x-, y-, and z-axis individually, and a rotation of θ degrees about the respective axes can be calculated using the × 3 3 rotation matrices R x , R y , and R z where The rotation matrices of equations (2)-(4) will rotate any column vector about the respective axis by the angle θ. The direction of the rotation is determined by the righthand rule, i.e., counterclockwise as seen from the positive end toward the origin. According to Euler, any orientation can be achieved by three successive rotations. However, because matrix multiplication is non-commutative, the sequence of rotations influences the final orientation of the body. Consequently, when three rotations are performed in succession following Euler's rules, there are still 12 possible combinations divided into two distinct groups: In the current work, the notation for rotation sequence is simply the axes in the order of which the rotations are performed with a subscript further emphasizing the sequence. For instance, the most common rotation following a proper Euler angle sequence is Z X Z 1 2 3 , which means first rotation about the z-axis, second rotation about the (new) x-axis, and finally another rotation about the (now current) z-axis. Table 2 tabulates all 12 rotation sequences following this notation.
The range of rotations is limited to the unit circle as rotations of more than π 2 radians make little sense. Hence, the range of α and γ can be restricted to [ ) π 0, 2 . Conversely, the range of β need not exceed π radians because larger rotations will only repeat previous spaces; hence, the range of β can be reduced to [ ) π 0, . In consequence, the orientation problem in AM is bound to the

Proposed method
The proposed optimization method is based on shape features to facilitate implementation in traditional tolerancing schemes. In this context, a shape feature is defined as either a geometric primitive, i.e., plane, cylinder, cone, sphere, or torus (Figure 2), or a free-form surface. These feature types are affected by the build direction in different ways and are also subject to different tolerance characteristics such as flatness, cylindricity, etc. Therefore, these features should be used in the construction of the objective function.
The current method requires information about the feature types and the relative orientation of the features. This information is readily available in many 3D-file formats such as STEP; however, when the geometry is converted into the popular intermediate STL format, the information about local topology is lost. For such tessellated file formats, the geometry may be deduced by feature recognition algorithms as described elsewhere [29][30][31].
Due to the ability to parameterize geometric primitives (e.g., height and diameter of a cylinder), these shape features may be generalized and are considered in the remainder of this article. The freeform surfaces on the other hand are more complex and therefore necessitate closer analysis and will be discussed only briefly. This prioritization is based on the assumption that surfaces that require tolerance are functional surfaces that predominantly are primitive in shape. The problem of freeform surfaces is left for future work.
Vectorial characterization of shape features enables the description of the relative location and orientation of shape features. Consider that each feature F is described by a location vector → p and an orientation vector 3 . This is demonstrated in Figure 3, where the PCS is located at the bottom left of the design and the black arrow represents the location vector for the highlighted horizontal through hole. The location vector points to the center of the hole where a feature coordinate system signifies the orientation of the feature, i.e., with the z-axis of the cylinder parallel to the y-axis of the PCS. When all features are described vectorially, the entire part P may be described as a set of all these features . The geometry of Figure 3 is recreated from the study by Cheng et al. [9] and will be used as a sample part for a simple geometry in this article.

Mathematical description of part build orientation
. These features may be of different sizes and feature types, i.e., planes, cylinders, spheres, cones, and tori. The quality of each feature type can be modeled as a function of some parameters . If a consistent measure of quality is applied (e.g., surface roughness Ra), the mean quality of the entire part Q part may be calculated as follows: where n is the number of features, Q i and A i are the quality and area of the ith feature, respectively, and A part is the total surface area of the part. While the area may be a suitable weight factor in many cases, it is also possible to introduce a separate weight factor that enables the prioritization of features. This weight factor can be applied on any level, e.g., certain feature types may be ignored, or individual features may attain a higher priority. Furthermore, the fundamental assumption is made that the quality of an additively manufactured surface can be described as a function of its orientation with respect to the build direction.
The orientation of a feature F with respect to the build direction ẑ may be described as two Euler angles , represents rotations about the xand y-axis respectively, both with reference to the original reference frame (i.e., Tait-Brian angles X Y Z 1 2 3 ). The angle θ between two arbitrary vectors → u and → v is found by: By utilizing unit vectors, | | | | → ⋅ → u v evaluates to 1 and reduces the expression to By inserting [ ] = ẑ 0 0 1 for → u , the aforementioned expression can be further reduced because In this context, z will be the z-component of the feature vector → v , denoted → v z . However, we want to express the angle θ as a function of the part's orientation in 3D space to enable optimization of orientation. For this purpose, the orientation vector can be expressed in terms of Tait-Bryan angles (Z Y X 1 2 3 ). These rotations are commonly referred to as yaw, pitch, and roll in engineering applications and describes the orientation of a rigid body with respect to the world coordinate system (WCS). The feature vector can then be derived from the rotation matrices in equations (2)-(4), and the sequential rotations about the x-, y-, and z-axis can then be performed in a single operation by multiplying the matrices as follows: This means that the orientation vector → v of a feature may be expressed as follows: which means When → v z from equation (11) is inserted in equation (7), the final expression simply becomes sin sin cos cos cos . (12) To enable the mathematical description of the entire geometry in a single expression, the orientation of each feature is described relative to a common coordinate frame, i.e., the WCS. The geometry may be regarded as a rigid body, which means that the relation between the surfaces remains constant and any transformation acts on all surfaces equally. In consequence, the orientation of all surfaces may be collectively calculated from the same values of α, β, and γ using equation (10), where x, y, and z are the components of the feature's initial orientation vector. Accordingly, equation (12) is valid for all surfaces with α and β as the only variables.

Finding the optimal part build orientation
Based on the aforementioned theory, it is possible to determine the optimal part build orientation mathematically by evaluating the critical points of the objective function. The critical points of a function ( ) f x y , are found, where Typically, there will be several solutions to equation (13), and each solution needs to be evaluated separately.
. These solutions will represent points, edges, and perhaps even entire areas in the 2D solution space. Provided the formalization in the previous section, the part build orientation can be described as a set of two rotations α and β. If a function is based on these two rotations, α and β will replace x and y in equation (13) in the search for critical points. The solution space will be the surface of a unit sphere, where [ ) ∈ α π 0, 2 and [ ) ∈ β π 0, . A point on this surface will correspond to a single unique orientation, while a line will correspond to a range of orientations. Note that the entire unit sphere is accessible already with only half a rotation of β, and still, a full rotation of β is used in the visualizations of this article to make their analysis more intuitive.
The fundamental assumption remains that the quality of a surface can be modeled as a function of its orientation with respect to the build orientation and that this function is differentiable. A conditional function (such as the one in equation (1)) introduces certain challenges to this method. However, such discontinuities would represent edges and areas in the 2D solution space that could be added to the list of candidate orientations.
Each feature of the geometry will add a term to the objective function, and each term will typically add one or more candidate orientations to the list. The exact number of additional candidates depends on the mathematical model as higher order functions will yield more candidates. It is therefore beneficial to limit or minimize the number of terms to avoid excessive computations. A simple way to minimize the number of terms is to join similar terms, e.g., two features of the same type and orientation can be combined into a single term of the objective function. Other measures include a manual selection of significant features and automatic filtering of features based on type, size, etc.
As surfaces are affected differently by build direction, separate models for each feature type are necessary to enable proper evaluation. This is easily implemented by inserting the appropriate expression for F Quality into equation (5) for each feature.

Case study
Two case studies are presented to demonstrate the proposed method: (1) A simple geometry to enable a step-by-step demonstration of the approach and all calculations. (2) A slightly more complex geometry to illustrate the applicability to more complex parts.
Simple mathematical models of quality are constructed for the illustrative purpose of this study. The implementation of empirical models is left for future work as this would obscure the central elements of this article. Before the case studies are presented, the construction of these mathematical models is detailed to provide the necessary foundations for the subsequent illustrations.

Constructing the mathematical models
The following models are based on the orientation rules described by Frank and Fadel [32], stating that cylinders should be oriented with the axis parallel to the build direction, while planes can be oriented either parallel or perpendicular to the build direction. Up-facing and down-facing surfaces are treated equally in the examples to keep the objective functions simple. This can be naturally incorporated in the objective function to account for any additional effects, e.g., overcure, support structures, and so on.
The central assumption in the current work is that the accuracy of a feature can be modeled as a function of its angle to the build direction (θ), which in turn is a function of the part's orientation as described in equation (12). Also, we are modeling deviations from nominal geometry, which can be considered a cost; hence, a low cost indicates a high fitness of a given orientation. However, the angle will not be sufficient to evaluate the fitness of a certain orientation. Consider, for instance, a cylinder oriented at a ∘ 45 angle from the build direction   Clearly, the objective function should be more sophisticated to incorporate this behavior.
Because the angle θ always will be in the interval [ ∘ ∘ 0 , 180 ], the sine of the angle will provide three desirable properties of an objective function: (i) the result is always a number between 0 and 1, (ii) the function is minimized at vertical orientations and maximized at horizontal orientations, and (iii) the function is periodic and symmetric. This study employs this function as an expression for the quality of cylinders: where θ is the angle between the feature normal vector and the build direction. This can be formulated as a function of ( ) α β , by inserting equation (12) Planes, on the other hand, require some additional configurations as we must consider both vertical and horizontal orientations as positive. To incorporate this new behavior, equation (14) is multiplied by θ cos 2 . This term ensures that planes are positively evaluated at both  vertical and horizontal orientations. In addition, the function is more sensitive to minor changes when a plane is horizontal than when the plane is vertical. This fit well together with the VE being large when the plane is close to horizontal, while not being as prominent in close-to-vertical orientations.
Finally, the function is normalized to facilitate comparison with other feature types, and so on. The normalization is achieved by introducing a divisor equal to the maximum of the function. The maximum is easily obtained by derivation and reveals a normalization factor of 2.598 when rounded to three decimal points. This yields the following expression for the quality of planes: where θ is the angle between the feature normal vector and the build direction. As with equation (14), the expression for planes may also be rewritten as a function of ( ) α β , by inserting equation (12) for θ. The result may be formulated as follows: The solution spaces of equations (15) and (17) are illustrated in Figure 4, where a 3D graph and a contour plot are presented for each of the equations.

Case 1: A simple geometry
The first case study is the simple geometry reconstructed from the study by Cheng et al. [9] and presented in Figure 3. This geometry provides a gentle introduction to the method by enabling a step-wise analysis of the geometry and the accuracy model. The first step of the method is to obtain a numeric description of the geometry in the appropriate format. Table 3 provides the positions and orientations of the features defined relative to the PCS as illustrated in Figure 3.
The data from Table 3 are inserted into equations (15) and (17) one row at a time as follows: (1) The feature type determines which equation to use (equation (15) for cylinders or equation (17) for planes) (2) The variables x, y, and z are substituted with E x , E y , and E z from Table 3 (3) The expression is multiplied with the relative surface area of the feature (final column of Table 3) Following the progression mentioned earlier, we start with the first row from Table 3 and perform the following steps: (1) The feature is identified as a planar type, and we, therefore, use equation (17). (2) The values of E x , E y , and E z Figure 7: Contour lines for the partial derivatives of equation (19). are inserted for x, y and z, respectively. For the first feature, this means that 0 is inserted for x and z, while 1 is inserted for y. (3) Finally, the weight factor is introduced by multiplying by the relative surface area as found in the final column of Table 3, namely, 0.093 (9.3%). These steps are displayed in the calculations of equation (18).
When all the data from Table 3 are inserted into equations (15) and (17), all the equations may be collected in a single expression for the entire geometry as follows: This yields the solution space illustrated in Figure 5. The solution space reflects the symmetry and regularity of the geometry as the orientations where feature vectors align with the build direction are clear.
In the next step, the partial derivatives are calculated and evaluated according to equation (13). Figure 6 shows the graphs of the partial derivatives where the dashed lines of Figure 6(a) and (b) correspond to the contour lines where the derivative evaluates to zero or are undefined.
The critical points are found where both derivatives either evaluate to zero or are undefined. This can be displayed graphically by plotting the dashed lines of Figure 6(a) and (b) in a single figure as shown in Figure 7. Finally, individual evaluation of these points and edges must be conducted to identify the global optimum as tabulated in Table 4.
The evaluation of [ ) ∈ α β π , 0,2 reveals four solutions with equal cost. However, these four solutions are pairwise identical, i.e., ( 270 , 0 ). Moreover, these two unique orientations are polar opposites, corresponding to the object lying on its left or right side as exemplified in Figure 8(a). This evaluation of the optimal orientation is also consistent with previous assessments of the same geometry [9,11,16].  Another orientation achieving a low cost is the upright position which is achieved for any value of α when β is ∘ 90 or ∘ 270 . This constitutes two edges in the solution space with identical solutions. These orientations correspond to the front or the back facing upwards, which aligns the largest cylinder with the build direction as displayed in Figure 8(b). Note that the up-facing and down-facing surfaces are not differentiated by the objective function. Incorporating this behavior would yield different results where repetition of the solution space is avoided.

Case 2: Example with more complex geometry
The second case is an original design created for the sole purpose of demonstrating the applicability of the proposed method on a slightly more complex geometry. The part displayed in Figure 9 is a joint with connectors in various directions. Forty-four features may be identified where 15 are cylindrical, and the remaining 29 are planes. A numeric description of the geometry is presented in Table 5.
Using the functions for planes and cylinders from equations (15) and (17) populated with the data of Table 5, the solution space of Figure 10 is obtained. Because many   features share orientation vectors, the overall objective function can be simplified to contain a minimal number of terms.
The solution space for the second case is also quite regular as demonstrated in Figure 11 where the partial derivatives of equation (20) are displayed. The regularity of these plots reflects the redundancy in the domain [ ] ∈ α β π , 0,2 . As demonstrated in case 1, the critical points are found where both partial derivatives are either zero or undefined, which corresponds to the intersections of red and black lines in Figure 12. All critical points are tabulated in Table 6 where the minima are highlighted in bold text.
For the second case study, the optimal orientation is achieved by ∘ 90 rotation about the x-axis as displayed in Figure 13(a). The same cost is also observed for three other combinations of rotations as displayed in Table 6. Due to redundancy in the solution space, these four combinations of rotations only correspond to two unique orientations in exactly opposite directions. With the hexagonal protrusions oriented parallel to the build direction, 25 out of the 29 planes are in a favorable orientation, while only two of the 15 cylinders are vertical.
Clearly, the large number of planes favors the orientation at ( ∘ ∘ 90 , 0 ). However, if cylinders are given a higher priority than planes, this orientation may no longer be as favorable. Figure 14 compares the effect of assigning a weight factor to cylindrical features for three orientations. A point of intersection is identified at a weight factor of 3.4, where the orientation ( ∘ ∘ 0 , 63 ) becomes more favorable than the previous best at ( ∘ ∘ 90 , 0 ) (see Figure 13(b)). This analysis indicates that for this part, cylinders should be at least 3.4 times as important relative to planes before the orientation at ( ∘ ∘ 0 , 63 ) is selected.
Both case studies exhibit repetitive symmetric patterns in the evaluation table. This symmetry arises from two sources: (i) as stated in Section 2.2, the domain of β can be constrained to [ ) π 0, as rotations exceeding this range will only repeat previous solutions, and (ii) the objective function of these case studies takes no regard of the difference between up-facing and down-facing surfaces. One may also question the alignment of many critical points with right-angle orientations. This is, however, a result of the initial orientation of the part, which originates from the design phase. If the initial orientation was different, the patterns observed in the plots would change due to the projection, but the solution space would remain unaltered.
A more complex geometry with more features would yield a solution space with more local extremes where the complexity of the solution space reflects the complexity of the geometry. Naturally, the objective functions for each surface type also contribute toward the topology of the solution space. Consequently, different objective functions, e.g., those obtained through experiments, would give different results to those reported earlier.
A mathematical description of the solution space provides a range of possibilities for finding the optimal orientation. What approach is best suited to solve the optimization problem depends on the purpose of the optimization as well as the complexity of the part. As the geometry becomes more complex, function evaluations are increasingly expensive.
Hence, complex geometries may benefit from intelligent methods where an effort is made to accelerate the computations. The simplicity of an exhaustive search may have certain advantages for simple geometries but may be ineffective for complex geometries.
The continuous solution space facilitates gradientbased methods as a qualified decision can be made concerning the next iteration of the search sequence. However,  unless all critical points are investigated, the risk of getting stuck in local optima is ever present. A stochastic component or a larger population of solutions may counter this problem, as is typically the case of EAs. When the solution space is clearly defined, it is possible to identify feasible regions in line with previous works [25,18]. These regions are defined from requirements and represent ranges of feasible orientations where tolerance requirements are met. Consequently, these regions can make up the boundaries for subsequent optimizations, e.g., with respect to mechanical properties. Adding steps to the optimization process will inevitably complicate and prolong the process, but interactive methods may be useful when the objectives are fuzzy or when flexibility is required.
Practical implementations of the proposed method would entail the formulation of objective functions for all relevant surface types. The relevant types and their definition may differ as long as they can be formulated as orientation vectors with accompanying objective functions. Furthermore, a solution for obtaining information on the orientation of the surfaces is required to automate the optimization process. An alternative implementation enables the identification of feasible regions for subsequent optimization. This would imply the integration of the method in a larger system with capabilities beyond what can be expressed by surfaces and their orientations.
The proposed optimization method has the benefit of being stable (i.e., no stochastic components), has no limitations with regards to search grid resolution, and provides the flexibility to incorporate separate expressions for each feature. Moreover, by considering the features of the part rather than every facet of a tessellated file, the effect of build direction may be generalized for each surface type. This drastically reduces the number of function evaluations. The feature-based approach also enables feature dimensions to be included in the objective function, which may affect the outcome.
The proposed method is not applicable to freeform surfaces due to the inability to formulate proper objective functions to handle the unknown. At present, this is not an issue for assembly features. However, as the potential of AM is unlocked, more complex surfaces may become widely used in the industrial design. The development of flexible formulations to handle this challenge is left for future work, along with the appropriate parametrization of such surfaces.

Conclusions and outlook
The approach described in this article provides the mathematical foundations for both deterministic and stochastic solutions to the orientation problem. By describing the solution space as a continuous function in a closed domain, the optimization can be performed mathematically. More importantly, this approach enables the determination of feasible orientation zones for the optimization of secondary objectives. Under the assumption that quality can be described as a function of build direction, the proposed method can be populated with any mathematical description of the relationship between quality and part build orientation.
The development of accurate mathematical models is crucial for optimization in process planning. AM technologies comprise many different processes that require separate models for predicting final part properties. Future research entails developing and validating prediction models that can be utilized in optimization processes. Practical implementation in a system with capabilities beyond the described method constitutes an interesting avenue for future research. Furthermore, the integration of all processing stages into a digital pipelinefrom design and process planning to quality assessment and verificationwill enable traceability throughout the manufacturing system and ultimately the entire product life cycle.

Acknowledgements:
The authors appreciate the support from colleagues at the Department of Manufacturing and Civil Engineering. Author contributions: TLL conceptualized the presented method with technical and theoretical support from OS. OS reviewed and edited the original draft put forward by TLL. Figure 14: Comparison of three orientations for case part 2 when an additional weight factor is introduced for cylinders.