Optimal design for disc golf by computational fluid dynamics and machine learning

In this article, we introduce a computational methodology for golf disc shape optimization that employs a novel disc shape parameterization by cubic B-splines. Through application of batch Computational Fluid Dynamics simulations and Machine Learning, the disc parameterization yields functional relationships—so-called shape surrogate models—between the flying rotating disc shape and its flight characteristics. The shape surrogate models facilitate free and constrained optimization in both single- and multiobjective settings, such that both aerodynamic (drag and lift) and structural (mass and moment of inertia) features of the disc are addressed simultaneously. Further, the Professional Disc Golf Association rules for permissible golf discs can be cast as nonlinear constraints for the computational optimization problem. The proposed numerical optimization method yields disc drag coefficient values as low as 0.48 (unconstrained) and 0.52 (constrained) and lift coefficient values as high as 0.26 (unconstrained) and 0.19 (constrained). The presented numerical optimization results also describe the many design tradeoffs between the discs that target long flight range (so-called drivers) and the discs that target flight at low speeds (so-called putters). Moreover, novel optimal rule compliant designs are presented for driver-type and putter-type discs, as well as their compromise, the so-called mid-range discs.


List of symbols
Distance between P-k and X-axis Growth parameter used in parameterization g Symmetry curve for a disc , i.e., intersection of and XY−planê Any disc characteristic (y) that only depends on disc shape 1 Introduction

Background
Disc golf is a sport that is similar to traditional golf, but instead of a club and a ball, players use rotating flying discs thrown toward standing baskets in as few intermediate steps as possible. Although perhaps not yet as popular as traditional golf, the disc sport is spreading rapidly across the world. In 2020, the Professional Disc Golf Association (PDGA) reported 71016 active members worldwide, 33% growth from 2019. Further, the number of recreational disc golf players is considerably larger than this. For example, in Finland alone, single-disc golf courses have reportedly been visited by over 2000 players on a weekly basis in 2020. While the ball shape and appearance in traditional golf are rather fixed (clubs do vary), considerable shape variation is allowed in golf discs. Indeed, players attempt to choose the most suitable disc for each throw situation and course layout, by exploiting the discs' different flight characteristics. Golf discs are typically categorized as drivers, mid-rangers, or putters, and there are currently over 1300 different-shaped discs approved by PDGA. However, perhaps due to commercial interests or complexity of the physics of rotating flight, surprisingly little academic research has been reported on the disc design principles and methodologies, which are the focus of this article.
To describe their products, golf disc manufacturers typically list four integer parameters called speed, glide, turn, and fade. The first two relate to the disc's aerodynamic drag and lift characteristics in horizontal flight, whereas the latter two relate to its pitching moment in inclined flight conditions (Lorenz 2007;Menickelli and Pickens 2016). Unfortunately, the reported numerical values of these four parameters are not based on measurements, but opinions of testers nor are such values standardized across different disc manufacturers. Thus, in particular, the listed parameter values do not directly tell a player which disc he or she should use in any given throw situation. The problem is compounded by individual variation: An amateur player (with limited power and technique) might be better off with a different disc than a professional player, for (say) covering the most distance in a single throw (Menickelli and Pickens 2016). These observations make development of a systematic golf disc design methodology both challenging and an important research question, which is addressed herein.

Contribution and limitations
The purpose of the present article is to introduce a computational methodology for golf disc shape optimization, based on Computational Fluid Dynamics (CFD) and Machine Learning (ML). The proposed methodology employs a novel disc shape parameterization by cubic B-splines that can be harnessed for batch solution of CFD simulation cases. In contrast to the four vague disc shape parameters listed above (speed, glide, turn, fade), the proposed parameterization yields precise functional relationships between the disc shape and its flight characteristics. These relationships, or shape surrogate models (SSM), are first identified in the present article through extensive CFD simulations and ML methods and are then used in numerical optimization.
The SSMs facilitate free and constrained optimization in both single-and multiobjective settings, such that both aerodynamic (drag and lift) and structural (mass and moment of inertia) features of the shape are addressed simultaneously. The upshot is that the PDGA rules for permissible golf discs can be cast as nonlinear constraints for the computational optimization problem.
To demonstrate the efficacy of the proposed methodology, we carry out numerical solution of several optimization problems. They aim to -Establish general principles of rotating flying disc design for improving flight characteristics, without constraints on the shapes beyond a minimum thickness; -Describe tradeoffs between minimum-drag and maximum-lift shape solutions (see Subsect. 2.1); -Resolve such golf disc shapes that both meet the PDGA rules, cast as numerical optimization constraints, and yield the longest range (drivers), best flight at low speeds (putters), or their compromise (mid-rangers).
Throughout this analysis, we confine our attention to horizontal disc flight, with zero angle of attack (AoA). Moreover, we assume stationary air (i.e., windless) conditions. It is well known (see Subsect. 2.1) that, even in such ideal conditions, gyroscopic precession tends to change the disc's alignment during rotating flight. Indeed, disc golf players occasionally exploit this tendency by imposing an initial inclination in throw or by choosing a disc with significant curvature in flight trajectory. In the present work, one part of the disc shape optimization problem is to minimize this tendency for alignment change during flight by maximizing the disc's moment of inertia with respect to the axis of forced rotation. Thus, we optimize discs for linear motion, which, by the above, is not always preferred by players in practice. Moreover, we neglect all grip-related preferences in golf disc design and assume that the presented disc flight characteristics can be realized in actual throws. The approach proposed herein is computational, but we emphasize that the CFD-based golf disc aerodynamics modeling method we utilize is extensively validated by comparing our CFD predictions to measurement data and earlier CFD analyses reported in the literature. Although the optimal disc shape solutions presented in this article comply with first principles of fluid mechanics, they are not guaranteed to be globally optimal ones. This is because the proposed disc shape parameterization does not cover all possible disc shapes and because the optimization algorithms are not guaranteed to yield global solutions for such nonlinear problems. Finally, in this article, we do not address the manufacturability of the presented optimal shape discs beyond the requirement that the minimum surface thickness must always exceed a given threshold of 0.5 mm (see Subsect. 2.5.2).

Rotational flight
The dynamics of rotational flight is relevant for many different sports (Goff 2013), among others. Its relative complexity arises from the interplay of fluid mechanics-i.e., lift and drag-and structural mechanics-i.e., moments induced by rotation and translation. At the turn of the millennium, several researchers worldwide started working toward understanding these phenomena for rotating discs, using experimental and (then emerging) computational approaches. The history of flying discs up to that era is well documented in the work of Potts and Crowther (2002) and Potts (2005) who also reported extensive experimental results on the aerodynamics of flying discs. Their results were subsequently complemented by Kamaruddin (2011) who experimentally evaluated the effect of different disc shape parameters on the aerodynamics. One important conclusion of these experimental results, which is also at the core of the optimization approach proposed in this article, is that the aerodynamic lift and drag coefficients (cf. Subsect. 2.1) are independent of the spinning motion of the disc (Potts and Crowther 2002).
The experimental research cited above was carried out in wind tunnel laboratory conditions. Carrying out experimental work on rotating flying discs outside wind tunnels is not an easy challenge. We mention the recent work of Weizman et al. (2020) and that of Lee et al. (2018). Both proposed sensor architectures and algorithms for estimating the aerodynamic coefficients. Instrumentation was also considered by Mustafin et al. (2012). In the present article, we focus on computational optimization and validate our results by measurements and CFD simulation results reported in the academic literature.
CFD simulations facilitate determination of aerodynamic coefficients for different disc geometries, with good accuracy (Potts and Masters 2015;Lukes et al. 2014;Honeycutt 2020;Jin et al. 2020). Moreover, typically this computational approach requires less time and effort and less expensive research infrastructure than is required for experimental tests. Many other rotating structures immersed in air flow, such as Flettner rotors, have also been successfully studied by CFD methods (De Marco et al. 2016). We utilize these results in selecting and validating the proposed CFD simulation methodology.
If the aerodynamic and structural parameters of a flying rotating disc are known, then computer simulations can be used for predicting its flight trajectory (although validation of such trajectory predictions by experiments is, by the above, difficult). Since the pioneering work of Hummel (2003) and Hubbard and Hummel (2000), several simulation-based approaches have been proposed, including validation (Koyanagi et al. 2012). Crowther and Potts 12 Page 4 of 17 (2007) described a general six degrees-of-freedom simulation approach; Hummel and Hubbard (2002) proposed an algorithm for estimating the aerodynamic coefficients from measured flight trajectory data. Lissaman and Hubbard (2010) studied the maximum ranges of flying discs, and Kamaruddin (2011) addressed optimal initial launch conditions for maximal flight range. Optimal flight conditions for discus were considered in simulations by Hubbard and Cheng (2007). In the present article, we address disc trajectories indirectly by aerodynamic and structural coefficients, which are represented by parametric SSMs used for optimizing the disc shape for linear motion at different disc speeds.
Recently, Kamaruddin et al. (2018) investigated the aerodynamic properties of golf discs using experimental methods and carried out six degrees-of-freedom simulations for resolving disc trajectories. They suggested that, to improve the flight range of a disc, one should maximize its lift-todrag ratio c l ∕c d (see Subsect. 2.1 for definitions). Moreover, they suggested that, to retain a straight-line flight path for a disc, its coefficient of pitching moment should be minimized. This article complements their findings by providing a computational framework for optimizing disc shapes with respect to several similar objectives, while also addressing shape constraints imposed by PDGA.

CFD-based shape optimization
Several authors have utilized CFD methods for shape optimization. For an overview, we refer the reader to the work of Immonen (2017Immonen ( , 2018, Immonen et al. (2021), and the references therein. While in many practical shape optimization problems, fluid motion is complex [e.g., (Zong et al. 2018)] requiring specialist analysis methods [e.g., (Immonen 2019)], in disc golf the fluid is turbulent single-phase flow, which is elementary for modern CFD solvers. The complexity of the problem here arises from rotational motion and the coupling of structural and aerodynamic properties, which depend on the disc shape. Moreover, these dependencies are nonlinear and can have discontinuities, as demonstrated by Marchand et al. (2017).
The shape optimization method proposed in the present article falls into the category of parametric surrogate methods, which have been successfully used for CFD-based optimization for both internal and external flows. Such methods attempt to first establish a computationally simple approximate mathematical relationship, the surrogate model (referred to as SSM in this article), between the design inputs and the optimization objectives, and then carry out optimization on the surrogate. Among the many successful applications in this category, we mention aircraft engine nacelles (Song and Keane 2007), gas turbine pin-fin arrays (Ghosh et al. 2020), wing aerodynamics (Immonen 2017;Zhonghua et al. 2020;Liu et al. 2017), fluidic oscillators (Jeong and Kim 2018), and elbow draft tubes (Demirel et al. 2017). However, to the author's knowledge, shape optimization for rotating flying discs has not been considered (by any CFD-based method) before. Although potentially requiring extensive simulation data to identify, our motivation for carrying out shape optimization by the surrogate approach, instead of nonparametric alternatives, arises from flexibility: not only can complex nonlinear design constraints be easily included in the analysis (Immonen 2017), but aerodynamic and structural objectives can also be easily treated simultaneously in multiobjective optimization (see Subsect. 3.3.2).
Linear and quadratic polynomial surrogates [see e.g., (Immonen 2017;Mohammad et al. 2021)] are among the most widely used objective function model types in the literature because of their simplicity and because they require relatively little data for model identification. This is beneficial when data generation is carried out by computationally intensive CFD simulations, although methods for addressing efficient data generation have also been presented recently (Yan et al. 2019;Li et al. 2020). Quadratic polynomial surrogates are also able to represent simple nonlinear dependencies between the design inputs and design objectives. However, such simple surrogate models are apparently not sufficient for all CFD-based shape optimization applications. Indeed, Immonen (2018) utilized high-order polynomial surrogates for S-duct shape optimization and many authors are increasingly also utilizing more complex ML-based surrogates. Among others, Jahangirian and Shahrokhi (2011) utilized multilayer perceptron neural networks for describing the ratio c l ∕c d of lift and drag coefficients (cf. Subsect. 2.1) as a function of shape parameters in airfoil shape optimization. Further, Zhang et al. (2021) considered deep neural network surrogate models in the same application domai, and Bouhlel et al. (2020) utilized gradient-enhanced neural networks for airfoil shape optimization at different flight conditions. In the present article, we utilize Regression Tree Ensembles (RTE), with optimized hyperparameters, as ML-based surrogate models (cf. Subsect. 2.4). Ensemble methods have displayed robust performance across a spectrum of CFD-related optimization and prediction applications (Immonen et al. 2020;Molinaro et al. 2021;Umetani and Bickel 2018;Mondal et al. 2019;Wang et al. 2017), although to the author's knowledge they have not been utilized for describing parametric shapes in aerodynamic optimization before.

Organization of the article
This article is organized as follows. In Sect. 2, we introduce elements of the computational methodology: The kinematic disc characteristics (Subsect. 2.1), CFD simulation models for resolving them (Subsect. 2.2), disc shape parameterization (Subsect. 2.3), SSM representations (Subsect. 2.4), and optimization (Subsect. 2.5). We then apply this computational methodology to flying disc shape optimization in Sect. 3, with optimization stages described in Subsect. 3.1 and results from SSM identification and numerical optimization covered in Subsect. 3.2 and Subsect. 3.3, respectively. Golf disc design principles are discussed in Subsect. 3.4, along with a detailed analysis of the optimal disc shape solutions. Finally, the conclusions of this study are presented in Sect. 4.

Kinematic characteristics of flying discs
In this subsection, we mathematically describe those characteristics affecting the flight of rotating discs that will be later used in shape optimization. Consider, initially, discs denoted by and defined as homogeneous (density kg∕m 3 ) solid bodies that are rotationally symmetric about the X-axis, as in Fig. 1. In general, the disc motion is both translational and rotational: it is linear along the Z-axis at velocity v m∕s and rotational about the X-axis at angular velocity r rad∕s . These considerations imply that the disc's AoA is 0°, a standing assumption for the remainder of this article.
Air resistance induces a drag force F d [N] along the negative Z-axis that reduces the disc's linear velocity. It is defined by where a kg/m 3 is the density of air, A Z [m 2 ] is the disc's projection area in the Z-direction, and c d > 0 is the dimensionless drag coefficient which depends on the shape of the disc. Similarly, denoting the disc's projection area in the X-direction by A X [m 2 ], the disc motion is affected by a lift force F l [N] in the X-direction. It is defined by where c l ∈ ℝ is the dimensionless lift coefficient. It, also, depends on the shape of the disc and can be positive or negative. In the X-direction, the disc motion is also affected by the disc mass m [kg], which is a design parameter.
Several studies have shown that c d and c l are independent of r (Potts and Crowther 2002;Potts 2005;Kamaruddin 2011). Consequently, disc shape can be optimized with respect to c d and c l by assuming r = 0 . We make this standing assumption in the remainder of the article. The shape optimization problem thus becomes simpler both computationally and conceptually, as the effects of structural mechanics and fluid mechanics are less coupled in the analysis. In the remainder of this subsection, we discuss the implications of assuming r = 0 and how the present article addresses them.
By setting r = 0 , we neglect the deceleration of the rotation speed r due to air viscosity and the sideways force in the Y-direction. Potts and Crowther (2002) demonstrated that the sideways forces are negligible, compared to forces in the is the disc diameter. At higher rotation speeds, the Magnus effect becomes more visible, but it is not considered herein.
For nonzero v, the pressure distribution on the surfaces of the disc is not symmetric about the axis of rotation (Potts 2005;Kamaruddin 2011); see also Fig. 4b. Thus, the corresponding center of pressure (cop) does not coincide with the center of gravity (cog). Consequently, the lift force, acting through the cop in the X-direction, induces a nonzero pitching moment M [Nm] that attempts to change the disc's alignment at rate q rad∕s and, as a result, the flight trajectory. Whenever r > 0 , gyroscopic effects would dictate that this pitching moment M also results in a precessional rolling moment L [Nm] at rate p rad∕s . Potts (2005) and Kamaruddin (2011) derived the following simplified equations for p and q: where I X kg m 2 is the disc's moment of inertia about the X-axis: Here, M and L depend crucially on the flow conditions around the disc: They vary throughout the flight and are difficult to control. The moment of inertia I X , on the other hand, is a static design parameter that only depends on the disc's Flying rotating disc: Linear motion is along the positive Z-axis and rotation is around X-axis, which passes through the disc's center of gravity shape. It can, thus, be maximized by shape design in order to minimize the rates p and q in Eq. (3), which improves the disc's linear flight characteristics; this approach is taken in the present article. In summary, the disc characteristics relevant for shape optimization in the present article are as follows: c d and c l (aerodynamical characteristics) as well as I X and m (structural characteristics).

CFD analysis of golf discs
In this section, we describe the CFD simulation methodology utilized throughout this article for predicting the drag and lift coefficients, c d and c l , corresponding to a disc shape profile.

Virtual wind tunnel CFD model
For aerodynamics modeling, the flying disc setup described in Subsect. 2 is equivalent to a stationary disc with an incoming flow at speed v along the Z-axis. Thus, we consider the rectangular virtual wind tunnel arrangement as in Fig. 2. Similar to the approach of De Marco et al. (2016) and Immonen (2017Immonen ( , 2018; Immonen et al. (2021), here the disc is modeled as a (stationary) no-slip wall. The virtual wind tunnel contains a velocity inlet surface (velocity v m/s, normal to boundary) and a corresponding pressure outlet surface (0 Pa) at the opposing end, representing an opening. The four other planar surfaces of the virtual wind tunnel are symmetry zones, which imply, in particular, that there is no boundary layer development (that would possibly introduce a pressure disturbance on the disc) on these surfaces. Determination of c d and c l , as in Eqs. (1) and (2), respectively, for any given disc shape involves regressing the drag and lift forces F d and F l , obtained from CFD simulation results, against the corresponding flow velocities v.

Mesh generation
For shape optimization, it is important to establish the aerodynamic performance of the different shape candidates relative to each other. As described by Immonen (2017Immonen ( , 2018 and Immonen et al. (2021), this can be achieved by keeping the simulation conditions, including the computational mesh and iteration count, as constant (i.e., equivalent) as possible across the shape candidates. The upshot is that establishing relative performance is possible even if the CFD model accuracy, in absolute terms, is coarse. As shape optimization in the present context would involve thousands of simulation cases, in the mesh selection process, we favored a reasonable total computation time (i.e., smaller number of cells) over a high accuracy in absolute terms (i.e., larger number of cells). Another practical requirement was that a large number of equivalent meshes, for varying disc shapes, were to be generated without errors in an automated procedure. In the present context, this means that the local mesh dimensions, e.g., maximum edge length, must be compatible with the disc geometries obtained from the parameterization (described in Subsect. 2.3).
The resulting meshes used throughout this study consist of approximately 360000 cells each. There are polyhedral cells in the bulk flow domain and a prismatic boundary layer of 5 layers, at 1.3 growth rate, spread around the disc cavity. An example of such mesh structure is displayed in Fig. 3. Equivalent meshes with this structure were generated in ANSYS Fluent Meshing by an automated batch process utilizing CAD geometry files created in SpaceClaim. During the course of this process, the mass m and moment of inertia I X were recorded from SpaceClaim for each disc shape.  Mesh structure in the vicinity of the disc for a disc geometry candidate. The disc interior is not included in the CFD model

CFD simulations
In numerical CFD simulations, we used ANSYS Fluent 2019R3 with double precision arithmetic, SIMPLE pressure-velocity coupling and second-order discretization for resolving the Reynolds-Averaged Navier-Stokes (RANS) equations. Since the flow around a disc is turbulent (Lukes et al. 2014;Kamaruddin et al. 2018;Potts and Masters 2015), a suitable turbulence closure equation for the RANS system was required. Based on the earlier observations of Immonen (2017Immonen ( , 2018, Immonen et al. (2021) and Potts and Masters (2015) from similar simulation cases, the k-SST model was chosen for resolving turbulence in the present work.

Validation studies
We conclude this section by presenting validation results for the proposed aerodynamic CFD modeling methodology. Table 1 presents a comparison of c d and c l predictions from CFD simulations to values reported in previous studies, for three different disc shapes. Here, c d and c l were calculated based on the single given reference velocity only, i.e., without regression.
In all simulations, AoA = 0 • and a = 1.225 kg∕m 3 . Further, 1000 steady-state iterations were sufficient in all cases to ensure solution convergence by the residual curves and the forces were calculated as averages from the last 100 iterations. All simulations were carried out using 40 cores in parallel. The total simulation time for each CFD case was approximately 10 minutes. Figure 4 illustrates the flow field prediction obtained from CFD simulations for a commercial disc geometry used as part of the validation study (Table 1,  row 3).
Clearly, the CFD simulation results are in reasonably good agreement with the values reported in the literature. Besides insufficient mesh resolution, the differences can be attributed to lack of geometry details and potential accuracy variations in measurements. Indeed, Potts (2005) described the difficulty of measuring forces on flexible, deforming, and rotating discs in high-speed flows.

Remark
In the validation studies, the area A X was used in both drag force equation (1) and lift force equation (2). Consequently, the values reported for c d in Table 1 are lower than those reported later in this article for optimal disc shapes.

Computational shape parameterization
Let S denote the set of all closed hypersurfaces of ℝ 3 with the following two properties: (i) each ∈ S is symmetric with respect to all rotations about the X-axis and (ii) the curve on the intersection of and the XY-plane is twice continuously differentiable and non-selfintersecting. In this article, every ∈ S denotes a disc surface and we seek to find those discs with the best flight characteristics. To do that, we assume a parameterization g ∶ ℝ P → S ∶ → = g( ) , by which is represented by P real numbers.
For optimization, a computational representation of the function g is required. Here it is obtained from the coordinates of 19 points, labeled P-k ( k = 1, … , 19 ), on a plane. These so-called control points yield the disc's intersecting curve via a cubic B-spline approximation, as illustrated in Fig. 5a-c. The curve is revolved around the vertical axis to produce the surface (and hence the disc ) in 3D for CFD analyses. Figure 5d illustrates the admissible regions for the control points and the shown asymptotic average disc shape is obtained from uniformly sampling these regions. Here, the 5 control points P-8...P-12 vary in the (overlapping) rectangular regions shown, whereas the other points are constrained: the point P-1 is always at (100, 0), the points (P-k, P-20-k), k = 1, … , 7 always share the same (and fixed) horizontal coordinate point, and are at an equal vertical distance (thickness parameter t) from each other. Moreover, the vertical locations of the constrained points are assumed to follow a power law , where is the distance from vertical axis and is a growth parameter. No self-intersections are allowed for the spline curves. This  (Table 1, row 3). Flow ( v = 27.5 , m/s) is from left to right. The rendered disc surface is shown to highlight that the simulation is carried out in 3D computational parameterization yields, in total, P = 12 design parameters for the admissible disc shapes.

Shape surrogate models
Let y be a unique characteristic that only depends on the shape of the disc , such as any one of the 4 characteristics c d , c l , m, I X listed in Subsect. 2.1. An SSM is a func- Typically the functions f employ so-called hyperparameters, which are specified during SSM identification, as explained below.

Multivariate polynomial SSMs
Multivariate polynomials (MVP) have been utilized as SSMs by several authors, e.g., Immonen (2017Immonen ( , 2018. Their identification by least-squares optimization is straightforward and there are many systematic methods for adding and removing terms in MVP models based on their statistical significance in explaining the response variable (Desboulets 2018). Moreover, unlike many ML-based models, MVP model predictions can, especially for low-order polynomials, be interpreted by polynomial coefficients and term interactions.  P-1 P-2 P-3 P-4 P-5 P-6 P-7 P-8 P-9 P-10 P-11 P-12 P-13 P-14 P-15 P-16 P-17 P-18 P-19 (c) (d) Fig. 5 a-c Example disc shapes obtained from parameterization g by varying the location of points P-1 to P-19 (coordinates in millimeters). d The admissible regions for points P-1 to P-19 displayed in colored areas. The average disc shape from uniform sampling of these regions is displayed in black for reference MVPs thus establish a baseline performance which more complex SSMs should improve to merit attention. The general representation of a MVP p ∶ ℝ P → ℝ yields a prediction ŷ for y as: where the coefficients a i 1 ,…,i P are to be identified. The degree of the polynomial p is a design hyperparameter.

Regression tree ensemble SSMs
A regression tree model (RTM) is a decision tree, obtained through supervised recursive binary partitioning (branching), that attempts to represent the response y as a function of by a sequence of if-then conditions. The well-known CART algorithm of Breiman et al. (Breiman et al. 1984) can be used for RTM identification. As in many ML models, it has parameters learned from data (such as node split points) and hyperparameters (such as node split criterion), called , that are chosen by the user based on the literature or direct optimization. For CART, they were summarized by Mantovani et al. (2018).
A Regression Tree Ensemble (RTE) is a predictive model composed of a weighted combination of multiple ( N > 1 ) regression trees, which attempt to combine results from many weak learners into one higher-quality ensemble predictor. Among the algorithms for generating ensembles, each containing hyperparameters that influence the prediction outcome in the present article, we consider bootstrap aggregation [bagging, cf. (Breiman 1996)] and least-squares boosting with gradient descent [lsboost, cf. (Hastie et al. 2009)].
Bagging involves generating bootstrap replicas, i.e., randomly resampling the data set, and training regression trees on the replicas, as in the so-called Random Forest concept (Breiman 2001). The hyperparameters controlling the procedure are the minimum leaf size L , i.e., the minimum number of observations per leaf node in a regression tree and the number of trees N in the ensemble. The prediction of the bagged RTE is: where f k ( , ) are the individual regression tree predictions, utilizing the hyperparameters .
The lsboost algorithm, on the other hand, progresses sequentially such that the ensemble iteratively (during up to M learning cycles) fits a new learner (regression tree) to the difference between the observed response and the aggregated prediction of all learners grown previously. The ensemble fits to minimize mean-squared error through gradient descent controlled by a hyperparameter called learning rate . The prediction of the lsboost RTE is obtained recursively as: where y is the sample average of y across shapes in the training data and f k ( , k ) are the individual regression tree predictions. The values ( k , k ) are obtained sequentially by the lsboost algorithm.

SSM hyperparameter optimization
Hyperparameter values control the learning process in RTE identification and they can be optimized analogously to ordinary parameters, in order to obtain a better fit for the SSM representations. In the present work, we used Bayesian optimization (Snoek et al. 2012) in MATLAB 2020b for obtaining the hyperparameters that minimize log[1 + CV(5)] , where CV(5) is the fivefold cross-validation loss corresponding to the model's output with the hyperparameters. The optimal hyperparameters are, by definition, those of the model which achieves the best cross-validation score. To keep the optimization time manageable, in practical optimization runs, the CART algorithm was used with default values and the optimization was carried out for the ensemble method (bagging vs. lsboost), L , M , N , and .

Shape optimization
In this article, we consider both free and PDGA-constrained shape optimization problems. In the former, the only requirement is disc manufacturability, which is here imposed as a requirement for the minimum-disc surface thickness t m to always exceed 0.5 mm. For computational analysis, t m is defined as the minimum pairwise distance between nonneighboring points on the B-spline approximation for . The PDGA-constrained optimization problem refers to meeting the PDGA global approval conditions (PDGA 2021). It is a set of 15 design constraints for golf discs approved for all competitions, consisting of maximum disc weight (200 g) and several dimensional constraints, summarized concisely in Fig. 6 for reference.
Page 10 of 17 Remark There is some room for interpretation in the PDGA rules, which makes their strict numerical implementation difficult. For example, the linear measurement tolerances are up to 1 mm, which can be over 6% of the disc height in practice. Further, referring to Figure 6, PDGA does not define the "uniform" flight plate. The PDGA rules also do not require a straight vertical disc section to be present for inner rim depth measurement, even though such is displayed in Fig. 6.
The general optimization problem considered here is: where K is the number of shape surrogate models considered in the optimization. In free-shape optimization, the only constraint is (9b), which is the minimum thickness requirement.
For practical optimization, the PDGA dimensional constraints were programmed to the function in Equation (9d). The PDGA disc weight constraint in Equation (9c) was included via a SSM for disc mass obtained from the disc volumes by assuming a constant polyethylene disc material ( = 920 kg m 3 ). In numerical optimization, we used = 0.7 min i ( i ) and = 1.4 max i ( i ) , where the minimum and maximum are componentwise and taken across the designs i used for SSM identification (see Fig. 5d and Sect. 3). Numerical optimization was carried out in MAT-LAB 2020b using the built-in solver functions, as listed in Table 2. The choice of solvers was based on experimentation and recommendations given in the software documentation. We emphasize that although the proposed solvers are shown to yield reasonable optimal solutions, they are not claimed to be superior to other algorithms; establishing the best optimization method for the present application is left for future work.

Stages of optimization
The proposed disc shape optimization method progresses in two stages: (1) training data of N samples are generated for SSM identification and (2) the SSMs are used for parametric optimization. The steps within these two stages are described in detail below and are illustrated schematically in Fig. 7.
H e r e i = 1, … , 12 and the minima and maxima are componentwise, as in Fig. 5d and Subsect. 2.5.2.  After Stage 1, there are N samples of c d , c l , m and I X , corresponding to different disc shape profiles , for SSM identification such that the underlying parameters k and k = 1, … ,N are orthogonal (because D is orthogonal). In this study, we carried out batch training data generation such that N = 700 , which resulted in 2800 CFD cases to be solved. The number N was chosen as a compromise: high enough for identifying MVP SSMs up to degree 3 (cubic) as baseline models and low enough to result in manageable total computation time (about 1 week). The remainder of this section describes the results obtained from the parametric optimization stage, using that training data.

Results from SSM identification
For parametric flying disc shape optimization, separate SSMs are required for describing c d , c l , m, and I X . Table 3 summarizes results from their numerical identification. Clearly quadratic MVP surrogates provide good fits for m and I X , in terms of both coefficient of determination R 2 and mean absolute error (MAE). However, ML-based models with optimized hyperparameters perform significantly better at representing c l and c d . It is also clear that, for c d and c l , the MVP representation does not improve by increasing the model complexity beyond polynomial degree 2. This indicates that the nonlinear dependencies of c d and c l on the flying disc shape parameters cannot be described by simple polynomial expressions and the proposed ML-based SSMs indeed merit attention. The SSM model specifications displayed in boldface in Table 3 are chosen for parametric optimization and are denoted by , , , and in what follows.

Results from shape optimization
In this subsection, we describe the numerical results from shape optimization. The corresponding optimal shape profiles are analyzed and compared in detail in Subsect. 3.4.

Free optimization
Numerical results from free-shape optimization are presented in Table 4 and Fig. 8. These results describe the general minimum-drag and maximum-lift disc shapes, without reference to the PDGA rules. In Table 4, the design objectives are displayed on the first column; for multiobjective problems (last 4 rows), the objectives are listed in a row vector. The corresponding optimal shape solution figures are referenced in the last column. In the numerical results, the superscripts SSM, CFD, and CAD denote the origin of the result, i.e., SSM, CFD prediction, or SpaceClaim 3D CAD geometry, respectively. The second column describes the focus of the optimization problem: Drag and lift can be optimized separately or simultaneously in single-objective problems (first 3 rows), whereas multiobjective optimization facilitates deciding upon the tradeoff from the Pareto efficient frontier. For multiobjective problems, Table 4 describes minimumdrag and maximum-lift solutions only, for brevity, as they are the most relevant for golf disc design.

PDGA-constrained optimization
Numerical results from PDGA-constrained shape optimization are presented in Table 5 and Fig. 9, which also display CFD prediction of the static pressure near the discs. Optimization results for driver, putter, and mid-ranger type discs are included. In each of the presented 3 cases, the optimization objective is min , − , − , as in the last 2 rows of Table 4. Each disc profile presented in Table 5 satisfies the PDGA rules (Subsect. 2.5.2). The terminology used in Table 5 is analogous to that for Table 4.

Analysis of the optimal disc shapes
We now compare and analyze in detail the optimal shape solutions presented in Tables 4 and 5. In general, the numerical accuracies of the SSM predictions reported for m and I X are better than those for c d and c l . Moreover, the mean absolute error, between SSM and CFD predictions, is lower for c l than for c d . These observations are in agreement with the SSM model fitness values reported in Table 3. Although there are some numerical inaccuracies in the prediction data, the corresponding optimal disc shapes display distinct     characteristics that are logical and compliant with the principles of fluid mechanics, as explained in the next subsections.

Free optimization (Table 4)
Focusing only on minimizing drag (row 1) or maximizing lift (row 2) yields shapes that display the best performance (across the presented optimal solutions) with respect to these objectives, respectively, on the SSMs. The minimum-drag disc shape (Fig. 8a) has a clear and horizontally almost symmetrical wedge toward both the upstream and downstream flows, which is in line with the well-known fact that minimum-drag shapes resemble the American football (Immonen et al. 2021). The minimumdrag disc shape also hosts a relatively wide rim section (Fig. 6). Such are often seen in driver-type golf discs targeting long-distance throws. A wide rim affects the disc flight in two distinct ways: it increases gyroscopic stability by increasing the disc mass away from the axis of rotation and it also delays the onset of flow separation from the bottom surface. On golf balls, the latter effect is produced by surface dimples, which are designed to enhance the flight distance.
The maximum-lift disc shape (Fig. 8b), on the other hand, has virtually no rim and has a uniformly curved upper section. Such curvature, also present on the top side of airfoils, increases the length of the streamlines traversing on the disc's top surface, which results in higher velocities-and hence a lower static pressure-there. The absence of a rim results in immediate flow separation on the bottom of the disc, slowing down the flow there through turbulent circulation and thus inducing higher static pressure underneath the disc. The net result of these effects is a higher pressure difference across the disc and thus a higher upward force, which is expected by the optimization objective. It should also be noted that the upper surface profile of a maximumlift disc is strictly decreasing along the disc radius, whereas  Table 5. Optimal shape profiles from free optimization. The control points P-1 to P-19 for cubic B-splines are shown for reference, except for (8c) where they are omitted for clarity for the minimum-drag disc profile (Fig. 8a), the top side has a local peak (reaching above the flight plate) at approximately the inner diameter location.
By the above, it is clear that minimization of the drag coefficient and maximization of the lift coefficient are contradictory objectives for flying rotating discs-a well-known fact for airfoils (Immonen 2017). In analogy to the recommendation of Kamaruddin et al. (2018) to maximize the liftto-drag ratio c l ∕c d , we consider minimizing ac d − (1 − a)c l for a ∈ [0, 1] to prioritize low drag ( a ≈ 1 ) or high lift ( a ≈ 0 ). The disc shape shown in Fig. 8c (solid line) and row 3 represents the equal-weight solution a = 0.5 . The cases a = 0.15 and a = 0.85 are included in Fig. 8c for illustration. Obviously the optimal solution at a = 0.5 is not a simple average of the optimal solutions at a = 0.15 and a = 0.85 , indicating nonlinearity in the optimization problem. The extent to which one objective can be improved without making the others worse is addressed by Pareto optimality in multiobjective optimization.
The Pareto set resulting from multiobjective optimization is a collection of optimal solutions in which no one criterion can be improved without making the others worse. Choice of a point from the Pareto set then reflects preference for some design objectives relative to the others. The optimal solutions shown in Fig. 8d (row 4) and e (row 5) represent preference for drag and lift, respectively, on the same Pareto set. The difference between Fig. 8a and d is that the latter solution cannot be improved with respect to drag without reducing lift. Similarly, the difference between Fig. 8b and e is that the latter solution cannot be improved with respect to lift without increasing drag. The last two rows (6 and 7) continue by including the  effect of moment of inertia in the optimization problem. Figure 8f shows a minimum-drag profile which cannot be improved without reducing either lift or moment of inertia (or both). Similarly, Fig. 8g shows a maximum-lift profile which cannot be improved without either increasing drag or reducing moment of inertia (or both). It is essentially such an extension of the high-lift optimal disc shape shown in Fig. 8e that also maximizes gyroscopic stability by the added mass on the rim section.

PDGA-constrained optimization (Table 5)
All 3 optimal solutions have remarkably smooth corners throughout their cross-section profiles, which is in contrast to many existing golf discs in practice. Smooth corners are aerodynamically better, but sharp corners may be preferred for a better grip, which is not addressed in the present study.
The optimal driver disc (targeting long flight range) on row 1 is one which minimizes c d and this feature cannot be improved at the associated I X and c l levels. Similarly, the optimal putter disc on row 2 (targeting flight at slow velocities) is one which maximizes lift and cannot be improved without compromising drag or moment of inertia. The optimal mid-ranger on row 3, finally, is a compromise of low drag and high lift at the given moment of inertia.
In each of the 3 cases, the moment of inertia with respect to the axis of rotation ( I X ) is maximized as part of the optimization problem, to minimize the perturbation effects induced by moments on the flight trajectory as in Eq. (3). This is seen in the wide rim sections for all 3 disc types, including the putter. Many existing putter-type discs only have a thin rim section, especially in comparison to drivertype discs; here, the difference in rim widths is small. However, similar to existing putter-type discs, the curvature on the top side of the putter disc is higher than for the driver and the mid-ranger.
For the driver disc profile, the static pressure profile (Fig. 9b) is remarkably uniform both above and below the disc, indicating closeness of cop and cog. Consequently, the moments that would change the disc's alignment via Eq.
(3) during flight in practice, are small. This is a desirable characteristic when attempting a long throw with a linear trajectory, but may be difficult to attain in practice (see e.g., Fig. 4b). Here, flight stability is further enhanced by high moment of inertia relative to the axis of rotation. For the putter (Fig. 9d) and mid-ranger (Fig. 9f) discs, the net vertical pressure difference across the disc is higher than for the driver, indicating a higher upward force at the same disc velocity. This is reflected in the highest c l value for the putter disc: it is the most capable of maintaining flight in low speed throws.

Conclusion
In this article, we have introduced a novel computational methodology for golf disc shape optimization, addressing structural and fluid mechanical characteristics simultaneously. The flying disc shapes arise here from parametric cubic B-spline representations. Together with an extensive set of CFD simulations they yield ML-based parametric descriptions-shape surrogate models or SSMs-for the main shape characteristics ( c d , c l , m, I X ) that affect the disc flight. These SSMs are can then directly be used in parametric disc shape optimization.
In the present article, we have carried out both single objective and multiobjective disc shape optimization, with and without design constraints. The design constraints considered herein are the PDGA rules for admissible competitive golf discs. The numerical optimization results display both general principles for golf disc design, targeting either long throws (driver-type discs) via a small drag coefficient or flight at low velocities (putter-type discs) via a large lift coefficient. The presented numerical results on constrained optimization display three novel golf disc profiles that meet the PDGA rules and target these objectives or their compromise (mid-ranger-type discs).
There are several interesting directions for future work. Applying new and powerful optimization algorithms instead of the off-the-shelf ones employed herein could potentially improve the optimal disc designs further. Such an analysis should also carefully address convergence and parameter sensitivity of the optimal solutions, which is not considered herein beyond the observation that the optimal disc shapes display design nonlinearity. Since the pitching and rolling moments are available from numerical fluid flow simulations, besides merely maximizing the moment of inertia, one could attempt to directly minimize pitching and rolling rates by identifying surrogate models for the moments. The physics of flight could also be made more realistic. In this article, we have confined our attention to zero angle of attack conditions. It would, thus, be interesting to carry out disc shape optimization such that it is robust with respect to nonzero angle of attack. Due to the high number of CFD simulation cases required, this would also highlight the need to incorporate intelligent (and adaptive) methods for data generation for SSM identification. Moreover, future work should be directed toward addressing sensitivity of the results to nonstationary air conditions. Finally, future work should focus on ensuring manufacturability, e.g., by 3D printing, of the optimal disc shape profiles, as well as discussing potential design modifications for practical throw situations (e.g., grip). The latter is only possible via physical prototypes, not considered herein, that would also make experimental validation of the methods of this work possible.