Shape-Informed Dimensional Reduction in Airfoil/Hydrofoil Modeling

1. Introduction

Free-form functional surfaces, such as wings, rudders, turbine blades, ship hulls, etc., play a vital role in the performance of a wide range of engineering products. Consequently, they have a significant impact on energy efficiency, the environmental footprint, along with structural integrity, safety, durability, and the financial viability of such endeavors. Such surfaces commonly possess complex geometries, which go far beyond simple primitive shapes, and therefore require advanced modeling techniques, e.g., parametric modeling, to efficiently represent and handle their free-form nature. Recent advancements in geometric modeling techniques, pioneered by mathematicians, computer scientists, and engineers, have significantly enhanced the capabilities of parametric modeling. This state-of-the-art technology effectively represents intricate shapes encountered in diverse fields such as naval architecture, ocean engineering, turbo-machinery, and automotive design [1,2].

To meet the need for cutting-edge designs, there is an increasing need for global optimization across increasingly wider design spaces with an ever-increasing number of design variables. However, such attempts are commonly hampered by the curse of dimensionality, which degrades the performance of optimization techniques as the dimension of the design space increases [3,4]. The first steps in design space dimensionality reduction involve methods which identify significant factors and variables and eliminate the remaining ones, by assigning appropriate constant values; see [5]. However, this approach is obviously simplistic and may fail to yield optimal results as it overlooks the potential importance of variables interdependency and their effect on optimization; see also [6]. Variance-based sensitivity analysis, employing Sobol indices [7], allows a more nuanced evaluation of the significance of the variables, but requires a statistically significant number of samples and becomes computationally intensive as dimensionality increases. To address such challenges dimensionality reduction techniques, ranging from principal component analysis to unsupervised learning and feature extraction techniques, that capture essential design space directions while preserving critical features, have been employed [8]. The last family of methods uncovers the underlying structures, enabling efficient optimization and knowledge acquisition.

Similarly, in the context of functional-surface design, high-dimensional simulation-driven design optimization problems can greatly benefit from the application of offline / upfront, design-space dimensionality reduction methods. Such contemporary methods assess the variability present within the design space and subsequently reduce its dimensionality accordingly. Hence, these techniques enable effective dimensionality reduction to be performed prior to initiating the optimization process. Diez et al. [9] have introduced a method based on the Karhunen-Loève expansion, also known as proper orthogonal decomposition, to assess shape modification variability and establish a latent space with reduced dimension for the shape modification vector; see also [10]. Despite the effectiveness of such approaches, they are not without drawbacks. One such drawback pertains to their failure to fully maintain the shape complexity and its underlying geometric structure. As a result, the resulting latent subspaces lack the ability to generate sufficiently diverse and/or valid shapes efficiently when used in shape optimization. This lack of representational capacity and compactness (i.e., shape validity) can hinder the success of the optimizer by spending the majority of the computational budget exploring infeasible or invalid shapes. In addition, these techniques often rely purely on geometric features and lack information related to performance criteria and the physics involved in the design assessment. Therefore, it is pivotal to employ more sophisticated representations that encompass high-level information regarding the shape’s structure and physics.

1.1. Proposed Approach

In this work we propose a methodology for tackling the aforementioned problems. The proposed approach augments the geometric shape description with integral properties, i.e., geometric moments, before constructing the design space of reduced dimension. Geometric moments allow us to capture both shape structure and performance information, as they are correlated to physical quantities (lift, drag, area, rigidity, etc.) which are assessed during foil optimization. Geometric moments provide important insights into the form, volume distribution, and potential performance characteristics of a design. Although the inclusion of the actual performance indices, as estimated via appropriate computational analysis, would be preferable, the cost of such simulations is, in most cases, significantly more expensive than moments’ calculation. Hence, constructing a latent space from a combined geometry/moment description results in a physics-informed, robust, and diverse design space, which can achieve high convergence rates when employed in design optimization.

Finally, as the employed dimensionality reduction approaches require discrete shape descriptions, we systematically study the effect of possible foil discretizations on dimensionality reduction outcomes and the quality of the resulting latent spaces.

1.2. Related Works

Locality Preserving Projection (LPP) [11] is a well-known unsupervised dimensionality reduction method based on feature extraction. It is a linear approximation of another common method known as Laplacian Eigenmaps (LE) [12]. The goal of LPP is to use a linear transformation to project the original data onto a lower-dimensional space while retaining links between nearest neighbors. In other words, LPP attempts to preserve local relationships and neighborhood structure during the dimensionality reduction process.

Several variants of LPP have been proposed to address specific limitations of the original method. One such variant is the ILLP-L1 [13], where the $L_2$ norm is replaced with the $L_1$ norm to enhance robustness. Despite LPP’s demonstrated effectiveness across various datasets, it is noteworthy that projection directions in LPP are highly sensitive to the coordinate system used for data points. Surprisingly, a mere rigid translation of data points, i.e., without altering their relative positions, can lead to significant changes in the projection directions determined by LPP. This indicates that LPP lacks the essential property of translation invariance, which is very significant in our case.

Generative adversarial networks (GANs) have been also used in pertinent literature to model forms of airfoils/hydrofoils. For example, in Chen et al. [14], authors demonstrated that their approach allowed them to achieve better optimization results with GAN mode shapes than they could with other parameterization techniques. At around the same time, Du et al. in [15] combined a B-spline representation approach with GANs, leading to the creation of BSplineGAN modes for foil design. However, the lack of prior knowledge regarding the proper number of modes to be used is a frequently voiced criticism against such dimension reduction strategies. Shape optimization with these reduced design spaces may fail to produce appropriate results when only a low number of modes is incorporated. On the other hand, using too many modes undermines the benefits gained by dimensionality reduction.

Using validity functions, which successfully eliminate aberrant foil designs, is another method to narrow the design space. This strategy depends on computationally efficient models to constrain the design space rather than reduce the number of design variables. Kedward et al. [16] presented a unique design constraint, based on curve derivatives/curvature, with the goal of ensuring smooth aerodynamic/hydrodynamic shapes during the shape design optimization process. This constraint significantly improved both the optimization’s convergence rate and the resulting optimized designs. Li et al. [17] developed a set of marginal functions to represent the link between dominant foil modes and higher-order modes derived from the UIUC [18] database airfoils. These functions effectively eliminated unwanted airfoils from the design space. However, defining acceptable boundaries for curvature-based constraints ahead of time is frequently difficult, and there was no guarantee that the margin-based validity functions did not also eliminate perfectly valid foil designs by accident. Hence, an ideal validity model should be an accurate discriminator of geometric irregularities without rejection of plausible foil designs.

Li et al. in [19] employed deep learning techniques, including a Deep Convolutional Generative Adversarial Network (DCGAN) for sampling realistic airfoils and a convolutional (CNN) discriminator to identify abnormal shapes, within a surrogate-based optimization framework for effective aerodynamic shape optimization. However, it is worth noting that the CNN-based validity function may have limitations when it comes to exploring geometric innovation beyond the scope of the UIUC airfoils. This is because the CNN model is trained using the airfoil GAN model, which itself is pre-trained exclusively on UIUC airfoils. Moreover, if the GAN model encounters modal collapse, the geometric filtering model could potentially exclude conventional airfoil shapes. Consequently, the deep-learning-based filtering model has faced criticism regarding its ability to hinder the discovery of optimal designs with innovative aerodynamic shapes. Jichao Li and Mengqi Zhang [20] introduced a GAN based on the Wasserstein metric (WGAN) to overcome the problem of model collapse observed in the airfoil GAN model. The WGAN-generated synthetic airfoils not only accurately represent the training airfoils but also extend into regions where there is a scarcity of training samples. As a result, the WGAN model demonstrated extrapolation capabilities, enabling it to capture additional geometric information beyond the ones present in the training set. Very recently, Shah et al. [21] presented an approach that constructs a generic parametric modeler on the basis of a deep convolutional generative adversarial network trained on a large dataset of ship-hull designs which had undergone several of the dimensionality reduction techniques discussed in this work.

2. Materials and Methods

For the application of the proposed dimensionality reduction method, we assume a rich space of foil-profile designs, $\left\{C\right\}$, which can be represented/modified by a design vector $\mathbf{v}\in \mathcal{V}\subseteq {\mathbb{R}}^n$. The design space $\mathcal{V}$ is then defined by an appropriate set of constraints that limit the design space to geometrically and physically valid foil profiles, i.e., $\mathcal{V}:=\{\mathbf{v}\in {\mathbb{R}}^n:g_i\left(\mathbf{v}\right)\le 0,i=1,\dots ,m\}$, where $g_i\left(\mathbf{v}\right),i=1,\dots ,m$ correspond to functional constraints capturing the design requirements. Furthermore, in design optimization, we assume at least one performance metric, e.g., $f:\mathcal{V}\to \mathbb{R}$, which evaluates the performance, $f\left(\mathbf{v}\right)$, of the corresponding design encoded by $\mathbf{v}$. Therefore, a simplified shape optimization problem targeting the maximization of the performance metric can be stated as:

$$\mathrm{Find}{\mathbf{v}}^{*}\in \mathcal{V}:f\left({\mathbf{v}}^{*}\right)=\underset{\mathbf{v}\in \mathcal{V}}{max}f\left(\mathbf{v}\right).$$

(1)

The computational cost of even the simple shape optimization in Eq. (1) increases exponentially with the dimension of $\mathcal{V}$, and increases even further when the evaluation of the performance metric, $f\left(\mathbf{v}\right)$, is complicated and time consuming. In this work, we begin by constructing the initial rich space of foil profile designs, $\left\{C\right\}$, with the help of publicly available foil design databases (UIUC foil database; see [18]), and its enrichment with the use of the parametric model described in Section 2.1.1. The details of its construction are discussed in Section 2.1 along with the applied discretization procedures required for the determination of $\mathbf{v}$ and its enrichment with appropriate integral shape properties (Geometric Moments) in Section 2.2. We conclude this section with a description of the approach followed for the proposed reduction of dimensionality in Section 2.3.

2.1. Design Data Set

Our first attempt to obtain a rich data set was solely based on random foil instances generated using the parametric model described in Section 2.1.1. Specifically, the generation process was initiated with a baseline design vector corresponding to NACA 2410 foil profile. Subsequently, 5000 random foil designs were generated by modifying the baseline design vector components in the range $\pm 50\%$ of their nominal values. However, the resulting design space proved to be insufficient as it lacked many profile families that are commonly found in pertinent literature. Instead of arbitrarily increasing the number of random designs in hope of reaching a sufficiently rich design space, we instead used, as a basis for our design space, the publicly available UIUC database of foil designs in [18]. This database comprises approximately 1500 foil profiles described by set of points with varying cardinality. This second approach poses two new problems: the number of foil designs is not sufficiently large, i.e., many gaps exist in the design space, and the profile representation is not only discrete, but it also lacks standardization, i.e., a common number of points or parameters describing each design. We tackle both problems with the same parametric model mentioned above. Specifically, the parametric model allows us to generate a high-accuracy continuous approximation of each foil profile as a NURBS curve, which can then be discretized, when needed, to any number of points. Additionally, we enrich the design space, by adding variations of each original profile shape. This is easily accomplished by modifying the design vector components in a range around the nominal values ($\pm 5\%$) for each design in the original database. Finally, a design space with 75781 unique profile shapes were generated and used in this work.

2.1.1. Employed Parametric Model

In this work, the airfoil/hydrofoil parametric model, firstly introduced in [22] and later extended to cover a wider range of designs in [23], has been extensively used, both for the generation of the original design space, and for providing benchmark shape optimization results; see Section 3. A wide variety of methods have been proposed in pertinent literature for the parametrization of airfoils and hydrofoils; see the works of Masters et al. [24,25] for a comprehensive review and comparison. However, the parametric model in [23] was specifically designed to accommodate the specific needs of design optimization, i.e, relatively small number of parameters, along with a high accuracy in representation of profile shapes.

The employed parametric model uses a range of a relatively small number of parameters that limits the dimensions of the design space. Specifically, the baseline parametric model with 11 parameters, presented in [23], can accurately describe a significant percentage (approximately 70%) of the profiles in UIUC database. The extended parametric models, in the same work, can represent the complete dataset, within Kulfan tolerance, with 15 or 28 parameters; see [23,26,27]. Furthermore, it employs parameters with physical meaning (widths, angles, lengths, etc.), which allow users to easily assign values to them, and generates shapes that do not exhibit undesirable oscillations, self-intersections, superfluous inflections points, or any other type of unwanted shape anomaly.

Each model-generated foil instance is represented by a NURBS curve (see [28]) of order k defined as follows:

$$\mathcal{C}(t;\mathbf{v})=(x(t;\mathbf{v}),y(t;\mathbf{v})):=\sum _{i=0}^n{\mathbf{b}}_i\left(\mathbf{v}\right)R_{i,k}\left(t\right),t\in [t_{k-1},t_{n+1}],$$

(2)

where, ${\left\{R_{i,k}\left(u\right)\right\}}_{i=0}^n$ is the set of rational B-spline basis of order k (degree $=k-1$) defined over the knot sequence, $\{t_0,t_1,t_2,.....,t_{n+k}\}$, and possessing nonnegative weights $w_i$, $i=0,1,2,\dots ,n$, while ${\mathbf{b}}_i\left(\mathbf{v}\right)$ are the corresponding control points determined by the parameter vector $\mathbf{v}$. In this work, we use a cubic representation ($k=4$) with 17 parameters and 13 control points ($n+1=13$); see Figure 1. Hence, we first find the parameters values for all designs in UIUC database and we then enrich the data set by generating a small number of variations, 5 for each foil profile in the database. Finally, we generate a rich superset of the UIUC database with 7578 foil profiles, represented by continuous cubic NURBS curves.

2.1.2. Design Discretization

The parametric modeler, as mentioned above, generates continuous curve representations that need to be discretized in the context of the dimensionality reduction techniques used in this work. The discretization in this case converts the continuous curve representation into a polygonal approximation, which can be stored as a vector of point coordinates. From a performance analysis perspective, we need to have a sufficient number of appropriately distributed points on the foil’s profile so that the performance evaluation, e.g., drag and lift estimation, produces an identical or insignificantly deviating result. In this work, XFOIL [29,30] has been used for foil performance evaluation, and 160 points, with a curvature-based distribution, have been proven in most cases sufficient to produce convergent performance results. However, we need to also account for the affect of the number of points and their distribution in the accurate representation of the data set that will be used for dimensionality reduction. Specifically, for minimizing the geometric deviation we pick a higher fixed number of points, i.e., 200, and examine the effect of points distribution by considering the following four cases:

• Uniform spacing over the parametric domain: We consider 200 parametric values, $\left\{t_i\right\}$ uniformly distributed on the parametric domain, $[t_{k-1},t_{n+1}]$; see Eq. (2). The corresponding 200 points on the profile are used in this case. Note that these 200 points are generally not uniformly distributed on the profile curve; see Figure 2a.
• Cosine spacing: We reparameterize the NURBS curve in Eq. (2) using the cosine function so that we achieve a concentration of points near the leading and trailing edge; see Figure 2b.
• Curvature-based spacing: This approach employs the profile’s curvature for calculating the distribution of parametric values. Specifically, parametric points are distributed in a way that guarantees an equal curvature integral for all parametric intervals. This approach results in a large concentration of curve points near the leading edge; see Figure 2c.
• Uniform spacing over the physical domain: This approach discretizes the profile by calculating segments of equal arc-length by dividing the curve length in 199 equal-length curve segments; see Figure 2d.

All discretization schemes mentioned above were applied in the representation of design space and their effect on the dimensionality reduction process and the quality of the latent space produced; see Section 3 for a detailed comparison.

2.2. Geometric Moments

The common design space dimensionality reduction techniques rely on the geometric representation of the design space (point sets, or polyhedral approximations in the general case) which quite often is insufficient for producing a latent space preserving fully shape complexity and the underlying geometric structure. This often leads to subspaces that lack diversity and/or capacity to generate fully valid shapes. Motivated by the work of Khan et al. [10], we enrich the design description with integral properties of the shape geometry. Specifically, the geometric information, residing in the discretized profiles, is coupled with shape moments leading to an enhanced design vector which more accurately describes the shape’s inherent characteristics and serves as its unique signature. Furthermore, the incorporation of geometric moments serves the additional goal of introducing quantities that are correlated to common performance criteria in a cost effective manner. More specifically, while most performance criteria (e.g., lift, drag, etc.) require generally costly simulations, some of these quantities are dependent on geometric moments and can be computed with a significantly lower computational cost. This approach has been also recently demonstrated in [31] for the case of airfoil and hydrofoil performance prediction and optimization where geometric moments have been included as input/output features when training the neural network models presented therein.

Assuming now that $\mathsf{\Omega }$ corresponds to the 2D domain bounded by the profile curve $\mathcal{C}=\partial \mathsf{\Omega }$, we may compute its $r^{th}$, $r=p+q$, order moments by:

$$\begin{array}{cc}\hfill M_{\left(r\right)}=M_{p,q}=& {\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{x}^p{y}^q\rho (x,y)dxdy,\hfill \\ \hfill & p,q\in \{0,1,2,\cdots \},\hfill \end{array}$$

(3)

where $\rho (x,y)=1$ when $(x,y)\in \mathsf{\Omega }$ and 0 otherwise. Obviously, the 2D domain area corresponds to the 0^th order moment, i.e., $\left|\mathsf{\Omega }\right|=M_{\left(0\right)}=M_{0,0}$, and its centroid can be calculated as $\mathbf{c}=(c_x,c_y)=\frac{M_{\left(1\right)}}{M_{\left(0\right)}}=\left(\frac{M_{1,0}}{M_{0,0}},\frac{M_{0,1}}{M_{0,0}}\right)$. A proper combination of geometry and its moments allows the construction of a vector that captures the shape’s inherent characteristics and serves as its unique identifier, the so-called signature vector in [32].

In our case, we may compute the Geometric Moments (GM) with the two-dimensional version of the divergence theorem, i.e., Green’s theorem. Specifically, using Green’s theorem and the representation of the foil profile in (2), we may equivalently write:

$$\begin{array}{ccc}\hfill & M_{0,0}=\frac{1}{2}{\int }_{t_{k-1}}^{t_{n+1}}x\left(t\right)\frac{dy\left(t\right)}{dt}-y\left(t\right)\frac{dy\left(t\right)}{dt}dt,\hfill & \hfill \\ \hfill & M_{1,0}=\frac{1}{2}{\int }_{t_{k-1}}^{t_{n+1}}x{\left(t\right)}^2\frac{dy\left(t\right)}{dt}dt,\hfill & \hfill M_{0,1}=-\frac{1}{2}{\int }_{t_{k-1}}^{t_{n+1}}y{\left(t\right)}^2\frac{dx\left(t\right)}{dt}dt,\\ \hfill & M_{2,0}=\frac{1}{3}{\int }_{t_{k-1}}^{t_{n+1}}x{\left(t\right)}^3\frac{dy\left(t\right)}{dt}dt,\hfill & \hfill M_{0,2}=-\frac{1}{3}{\int }_{t_{k-1}}^{t_{n+1}}y{\left(t\right)}^3\frac{dx\left(t\right)}{dt}dt,\\ \hfill & M_{1,1}=\frac{1}{2}{\int }_{t_{k-1}}^{t_{n+1}}x{\left(t\right)}^{2}y\left(t\right)\frac{dy\left(t\right)}{dt}dt,\hfill & \hfill M_{1,1}=-\frac{1}{2}{\int }_{t_{k-1}}^{t_{n+1}}y{\left(t\right)}^{2}x\left(t\right)\frac{dx\left(t\right)}{dt}dt,\\ \hfill \mathrm{etc}.& \hfill \end{array}$$

(4)

Finally, an approximate evaluation of the moments appearing in Eqs. (3) and (4) can be easily accomplished by employing a triangulation of $\mathsf{\Omega }$ or discretization of $\partial \mathsf{\Omega }\equiv \mathcal{C}\left(t\right)$, respectively.

The performance criteria considered in this work, i.e., lift and drag coefficient, are invariant to translation and uniform scaling, whereas the moments described above change value with translation (T) and scaling (S). Therefore, we need to employ TS-invariant Geometric Moments (TS-iGM). Xu and Li in [33] derive translation, scale, and rotation invariant moments (TSR-iGM). For our purposes, rotational invariance in unwanted and therefore, we employ central moments and normalization to achieve translational and scale invariance, respectively. Specifically, by using the centroid of $\mathsf{\Omega }$, i.e., $\mathbf{c}=(c_x,c_y)=\left(\frac{M_{1,0}}{M_{0,0}},\frac{M_{0,1}}{M_{0,0}}\right)$, we obtain the so-called central moments:

$${\overline{\mu }}_{p,q}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{(x-c_x)}^p{(y-c_y)}^q\rho (x,y)dxdy.$$

(5)

As for scaling, we can first check its effect on the value of central moments and normalize accordingly. Specifically, if we scale uniformly $\mathsf{\Omega }$ by a scale factor s, we get

$$\begin{array}{cc}\hfill {\widehat{\mu }}_{p,q}=& {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{(x-c_x)}^p{(y-c_y)}^q\rho (\frac{x}{s},\frac{y}{s})dxdy=\hfill \\ \hfill & {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{s}^p{(x-c_x)}^p{s}^q{(y-c_y)}^q\rho (x,s)s^{2}dxdy=s^{p+q+2}{\overline{\mu }}_{p,q}.\hfill \end{array}$$

We may now use any moment for the normalization, but we pick the 0^th order moment since it is easier and numerically more stable to compute, i.e., ${\widehat{\mu }}_{0,0}=s^2{\overline{\mu }}_{0,0}$. Hence, we can write:

$${\mu }_{p,q}=\frac{{\widehat{\mu }}_{p,q}}{{\left({\widehat{\mu }}_{0,0}\right)}^{\ell }}=\frac{s^{p+q+2}{\overline{\mu }}_{p,q}}{{\left(s^2{\overline{\mu }}_{0,0}\right)}^{\ell }}.$$

Finally, we set $\ell =\frac{p+q+2}{2}$ to eliminate the scale factor s and we get the $r^{th}$-order normalized central geometric moments

$${\mu }_{\left(r\right)}={\mu }_{p,q}=\frac{s^{p+q+2}{\overline{\mu }}_{p,q}}{{\left(s^2{\overline{\mu }}_{0,0}\right)}^{\frac{p+q+2}{2}}}=\frac{{\overline{\mu }}_{p,q}}{{\overline{\mu }}_{0,0}^{\frac{p+q+2}{2}}}.$$

(6)

As a side note, one may easily observe that for these TS-invariant moments, ${\mu }_{\left(0\right)}={\mu }_{0,0}=\frac{{\overline{\mu }}_{0,0}}{{\overline{\mu }}_{0,0}}=1$, and ${\mu }_{\left(1\right)}=(0,0)$.

2.3. Dimension Reduction

We assume here that our foil geometry corresponds to an appropriate vector of point coordinates $\overline{\mathit{\theta }}$. Additionally, a vector $\mathbf{v}\in \mathcal{V}$ (see Section 2 and Eq.(1)) allows to define a shape modification procedure: ${\overline{\mathit{\theta }}}^{†}=\overline{\mathit{\theta }}+{\mathbf{v}}_{\overline{\mathit{\theta }}}$, where ${\overline{\mathit{\theta }}}^{†}$ denotes the resulting shape when the modification vector, ${\mathbf{v}}_{\overline{\mathit{\theta }}}$, is applied to the original shape $\overline{\mathit{\theta }}$.

As previously discussed, the cost of shape optimization increases exponentially with the dimension of $\mathcal{V}$, especially when the performance criteria are many and computationally expensive. To overcome this issue, we propose to use feature extraction techniques to create a lower-dimensional latent-space. As already noted, we also intend to include additional features, i.e., geometric moments, to address some drawbacks of traditional Design Space Dimensionality Reduction (DSDR) techniques, and produce a rich latent-space which can be efficiently employed in shape optimization.

To create this latent-space, we combine ${\mathbf{v}}_{\overline{\mathit{\theta }}}$ with the corresponding geometric moment vector $\mathit{\mu }\left({\mathbf{v}}_{\overline{\mathit{\theta }}}\right)$. We can now define a combined geometry and moment vector $\mathit{p}\left({\mathbf{v}}_{\overline{\mathit{\theta }}},\mathit{\mu }\left({\mathbf{v}}_{\overline{\mathit{\theta }}}\right)\right)\in \mathcal{P}\subseteq {\mathbb{R}}^{n_p}$, where $n_p=n+n_{\mu }$, with n corresponding to the dimension of $\mathcal{V}$, and $n_{\mu }$ to the number of moments employed. For notational simplicity, we may now define $\mathit{\vartheta }=\left({\mathbf{v}}_{\overline{\mathit{\theta }}},\mathit{\mu }\left({\mathbf{v}}_{\overline{\mathit{\theta }}}\right)\right)$ and consequently the function $\mathit{p}\left(\mathit{\vartheta }\right)$, which encapsulates information about the foil’s shape and moments. This function belongs to a disjoint Hilbert space $L_f^2\left(\mathcal{P}\right)$ defined by a generalized inner product:

$${(\mathbf{a},\mathbf{b})}_f={\int }_{\mathcal{P}}f\left(\mathit{\vartheta }\right)\mathbf{a}\left(\mathit{\vartheta }\right)·\mathbf{b}\left(\mathit{\vartheta }\right)d\mathit{\vartheta }=f\left({\mathit{\mu }}_{\overline{\mathit{\theta }}}\right)\mathbf{a}\left({\mathit{\mu }}_{\overline{\mathit{\theta }}}\right)·\mathbf{b}\left({\mathit{\mu }}_{\overline{\mathit{\theta }}}\right)+{\int }_{\mathcal{V}}f\left(\overline{\mathit{\theta }}\right)\mathbf{a}\left(\overline{\mathit{\theta }}\right)·\mathbf{b}\left(\overline{\mathit{\theta }}\right)d\overline{\mathit{\theta }},$$

(7)

where the inner product is weighted by positive functions $f\left({\mathit{\mu }}_{\overline{\mathit{\theta }}}\right)$ and $f\left(\overline{\mathit{\theta }}\right)$ that correspond to moments and geometrical description, respectively.

When considering shape optimization in $\mathcal{V}$ one may employ engineers’ expert experience and identify regions that are more probable to produce the optimum design, i.e., the probability that $\mathbf{v}$ will correspond to the unknown optimum ${\mathbf{v}}^{*}$ varies in $\mathcal{V}$. However, as noted in [9,10], when we lack such prior knowledge, we may generally consider a uniform distribution for $\mathbf{v}$, i.e., we consider $\mathbf{v}$ to be an element of a stochastic space $\mathcal{V}$ with an associated probability density function $\rho \left(\mathbf{v}\right)$, which is uniform in our case. Under this assumption we can calculate the mean and variance of the signature vector function $\mathit{p}\left(\mathit{\vartheta }\right)$.

$$\langle \mathit{p}\rangle ={\int }_{\mathcal{V}}f\left(\mathit{\vartheta }\right)\mathit{p}\left(\mathit{\vartheta }\right)\rho \left(\mathbf{v}\right)d\mathbf{v},$$

(8)

$${\sigma }^2=\langle \parallel \overline{\mathit{p}}{\parallel }^2\rangle ={\int }_{\mathcal{V}}{\int }_{\mathcal{P}}f\left(\mathit{\vartheta }\right)\overline{\mathit{p}}\left(\mathit{\vartheta }\right)·\overline{\mathit{p}}\left(\mathit{\vartheta }\right)\rho \left(\mathbf{v}\right)d\mathit{\vartheta }d\mathbf{v},$$

(9)

where $\overline{\mathit{p}}$ is the deviation from the mean, i.e., $\overline{\mathit{p}}=\mathit{p}-\langle \mathit{p}\rangle$.

In dimensionality reduction, we are looking for a lower-dimensional vector $\mathbf{u}\in \mathcal{U}\subseteq {\mathbb{R}}^{\kappa }$, with $\kappa <<n$, which will substitute $\mathbf{v}\in \mathcal{V}\subseteq {\mathbb{R}}^n$ in the representation of $\overline{\mathit{p}}\left(\mathit{\vartheta }\right)$; recall that $\mathit{p}\left(\mathit{\vartheta }\right)\equiv \mathit{p}({\mathbf{v}}_{\overline{\mathit{\theta }}},{\mathit{\mu }}_{\overline{\mathit{\theta }}})$. We now consider this new $\kappa$-dimensional vector, the latent variable $\mathbf{u}$, which we want to construct with appropriate features from the original space. Consequently, if the latent space captures adequately the original space, this new latent variable will significantly reduce the computational burden in shape optimization since $\kappa <<n$.

2.3.1. Karhunen-Loève Expansion (KLE)

Once the initial shape signature vector (SSV), $\mathit{p}({\mathbf{v}}_{\overline{\mathit{\theta }}},{\mathit{\mu }}_{\overline{\mathit{\theta }}})$, is constructed to combine shape geometry with invariant geometric moments, we apply the so-called Karhunen-Loève Expansion (KLE) technique to find an appropriate basis of orthonormal functions to approximate our initial space, i.e.,

$$\overline{\mathit{p}}\left(\mathit{\vartheta }\right)\approx \sum _{i=1}^{\kappa }{u}_i{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right),$$

(10)

where ${\left\{{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right)\right\}}_{i=1}^{\kappa }$ is the set of the basis functions spanning $\mathcal{U}$, and $\mathbf{u}=(u_1,u_2,\dots ,u_{\kappa })$. We obviously want an appropriate latent space $\mathcal{U}$ that will preserve the variance, to the extent possible, exhibited in $\mathcal{P}$. Furthermore, for the components of $\mathbf{u}$, we get:

$$u_i={\left(\overline{\mathit{p}},{\mathit{\varphi }}_i\right)}_f={\int }_{\mathcal{P}}f\left(\mathit{\vartheta }\right)\overline{\mathit{p}}\left(\mathit{\vartheta }\right)·{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right)d\mathit{\vartheta }.$$

(11)

Additionally, as mentioned above, we aim in creating basis functions that preserve variance ${\sigma }^2$ in the original space; see also Eq. (9). If we now plug in Eqs. (10) and (11) in (9), it is easy to see that:

$${\sigma }^2=\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\langle u_i{u}_j\rangle {({\mathit{\varphi }}_i\left(\mathit{\vartheta }\right),{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right))}_f=\sum _{j=1}^{\infty }〈u_j^2〉=\sum _{j=1}^{\infty }〈{\left(\overline{\mathit{p}},{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right)\right)}_f^2〉.$$

(12)

As discussed in detail in [9,10], the basis exhibiting the required characteristics can be determined via the solution of the following constrained minimization problem:

$$\begin{array}{cc}& \underset{\mathit{\varphi }\in L_f^2\left(\mathcal{P}\right)}{min}〈{\left(\overline{\mathit{p}},\mathit{\varphi }\left(\mathit{\vartheta }\right)\right)}_f^2〉,\hfill \\ \hfill \mathrm{subject}\mathrm{to}:& {\left(\mathit{\varphi }\left(\mathit{\vartheta }\right),\mathit{\varphi }\left(\mathit{\vartheta }\right)\right)}_f^2=1.\hfill \end{array}$$

Finally, as demonstrated in [9,10], results in

$$\mathcal{L}\mathit{\varphi }\left(\mathit{\vartheta }\right)={\int }_{\mathcal{P}}f\left(\mathit{\theta }\right)〈\overline{\mathit{p}}\left(\mathit{\vartheta }\right)\otimes \overline{\mathit{p}}\left(\mathit{\theta }\right)〉\mathit{\varphi }\left(\mathit{\theta }\right)d\mathit{\theta }=\lambda \mathit{\varphi }\left(\mathit{\vartheta }\right),$$

(13)

where $\mathcal{L}$, ⊗ denote the self adjoint integral operator and the outer product, respectively. Both $\mathit{\theta }$ and $\mathit{\vartheta }$ belong to original space and $\left\{\lambda \right\}$ are the eigenvalues. Hence, the resulting eigenvectors ${\left\{{\mathit{\varphi }}_i\left(\mathit{\vartheta }\right)\right\}}_{i=1}^{\infty }$ correspond to the so-called KL-modes which are orthogonal and constitute the basis of $L_f^2\left(\mathcal{P}\right)$. Furthermore, the corresponding eigenvalues ${\left\{{\lambda }_i\right\}}_{i=1}^{\infty }$ represent the variance, i.e.,

$${\sigma }^2=\sum _{i=1}^{\infty }{\lambda }_i.$$

(14)

We now consider the first $\kappa$ eigenvectors form the basis for our approximation. Obviously, the selection of $\kappa$ depends on the desired level $\epsilon$, ($0<\epsilon \le 1$), of variance we want to capture. In other words, we can choose $\kappa$ such that the selected eigenvectors capture at least $(\epsilon \times 100)$ percent of the total variance in the data.

$$\sum _{i=1}^{\kappa }{\lambda }_i\ge \epsilon \sum _{i=1}^{\infty }{\lambda }_i=\epsilon {\sigma }^2,{\lambda }_i\ge {\lambda }_{i+1}.$$

(15)

As a final note, we need to add here that the method described in Diez et al. [9] and Khan et al. [10] is used to numerically solve Eq. (13).

2.4. Latent Space

Apart from the dimensionality reduction that leads to the construction of the latent space, as described in Section 3.2, we need to also consider the potential use of parameter bounds which will limit the exploration within a permissible space that lacks, or at least minimizes, the presence of infeasible or invalid designs. This step is quite critical since although in principle we want to be able to explore a large space, its size as well as the presence of invalid shapes at some regions constitute hampering factors in the subsequent design optimization. Traditionally, such bounds are defined via functional and/or side constraints that define the feasible design space. In the case of the original design space, the design vector $\mathbf{v}$ commonly comprises parameters with physical meaning and therefore introducing constraints can be a relatively straightforward task for designers.

However, there are two problems with this approach, the first one stems from the use of the latent space, while the second one pertains to the complexity of the resulting optimization problem. The latter problem is a common problem in constraint optimization where the inclusion of many and/or complicated functional constraints results in a complex and time-consuming optimization problem. Side constraints, i.e., bounds on parameter values, can be more easily handled by the optimization algorithms and are preferred when they can replace the more complicated, usually non-linear, functional constraints. However, the former problem, even when only considering side constraints, remains. Specifically, determining the side constraints, i.e., the parameter bounds $({\mathbf{u}}^{low},{\mathbf{u}}^{high})$, for the latent space poses additional challenges for designers since $\mathbf{u}$, unlike $\mathbf{v}$, is a latent design vector with no physical meaning, i.e, does not correspond to lengths, radii, angles, etc. Therefore, as no physical interpretation is generally associated with any component of $\mathbf{u}$, it is hard to determine the values of ${\mathbf{u}}^{low},{\mathbf{u}}^{high}$ which guarantee that designs within the feasible region of $\mathcal{U}$ are also valid designs in $\mathcal{V}$. In other words, it is hard to bound $\mathcal{U}$ in way that will guarantee the satisfaction of all design constraints and requirements in the original space.

A prevalent method for defining latent space bounds is to utilize the standard deviation from the mean shape, which is positioned at the center of the design space. This method achieves a balance between the number of infeasible shapes and the level of permissible diversity. We consequently follow this approach in our implementation with the bounds of the i^th component of the latent variable defines as

$$u_i\in \left[-\sqrt{\alpha {\lambda }_i},\sqrt{\alpha {\lambda }_i}\right],$$

(16)

with the $\alpha$ being an integer value between 1 and 3, with larger values allowing for larger spaces. In other words, each component of the latent design vector is bounded within a multiple of its deviation.

2.4.1. Latent Space Quality Analysis Measures

To evaluate the appropriateness and effectiveness of a latent space in generating a broad range of diverse and valid shapes, suitable metrics need to be introduced. A commonly used measure of diversity is based on Hausdorff distance [34], which is frequently employed in the quantification of distance between two subsets of a metric space. Similarly, it can also be employed to quantify the distance, i.e., similarity/dissimilarity, between two free-form shapes. In our case, we assume that both shapes can be appropriately discretized via a corresponding finite set of points, namely ${\widehat{C}}_1=\{{\mathbf{c}}_{1,i}:{\mathbf{c}}_1\in {\mathbb{R}}^2,i=1,2,\cdots ,n_{C_1}\}$ and ${\widehat{C}}_2=\{{\mathbf{c}}_{j,2}:{\mathbf{c}}_2\in {\mathbb{R}}^2,j=1,2,\cdots ,n_{C_2}\}$. We can then calculate the Hausdorff distance, H, between the finite sets ${\widehat{C}}_1$ and ${\widehat{C}}_2$ as follows:

$$\begin{array}{cc}\hfill H\left({\widehat{C}}_1,{\widehat{C}}_2\right)=& max\left(h({\widehat{C}}_1,{\widehat{C}}_2),h({\widehat{C}}_2,{\widehat{C}}_1\right),\mathrm{where}\hfill \\ \hfill h({\widehat{C}}_1,{\widehat{C}}_2)=& \underset{{\mathbf{c}}_1\in {\widehat{C}}_1}{max}\underset{{\mathbf{c}}_2\in {\widehat{C}}_2}{min}\parallel {\mathbf{c}}_1-{\mathbf{c}}_2\parallel ,\hfill \end{array}$$

(17)

where $h(A,B)$ is the directed Hausdorff distance from A to B, and $\parallel ·\parallel$ denotes the Euclidean norm. We further assume here that the Hausdorff distance between the two shapes is determined by the Hausdorff distance between the corresponding point sets as defined in Eq. (17). Finally, the diversity measure of a latent space is defined as the mean Hausdorff distance between a reference design and a set of densely sampled designs obtained from the latent space. Higher values for the diversity measure are desirable as they imply richer latent spaces; however this comes at a cost, i.e., higher chance of generating invalid geometries. Invalid shapes, in our case, include foils with self-intersections2, shapes with an oscillatory behavior, foils with upper or lower not corresponding to functions3, etc. Hence, the ideal latent space should not only have a high diversity value but a low or zero percentage of invalid geometries. We define here the validity measure of the latent space as the ratio of invalid designs to valid designs obtained from a dense sampling in the latent space. Obviously, a latent space with a validity measure close or equal to 0 is preferred.

2.5. Shape Optimization

After constructing the latent design space, $\mathcal{U}$, an optimization process is carried out to explore its capacity to substitute the original space when searching for optimal designs. We select the lift coefficient, $C_L$, as the quantity of interest (QoI) for the shape optimization assessment and the corresponding maximization problem can be formulated as follows:

$$\begin{array}{cc}\hfill \mathrm{Find}{\mathbf{u}}^{*}\in \mathcal{U}& \mathrm{such}\mathrm{that}\hfill \\ \hfill C_L\left({\mathbf{u}}^{*}\right)=& \underset{\mathbf{u}\in \mathcal{U}}{max}{C}_L\left(\mathbf{u}\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& |A\left(\mathbf{u}\right)-A_0|\le 0.001,\hfill \\ \hfill & {\mathbf{u}}^{low}\le \mathbf{u}\le {\mathbf{u}}^{high},\hfill \end{array}$$

(18)

where A corresponds to the airfoil’s area and the sub-index ${(}_0)$ denotes values for the reference design, i.e. NACA-2410 airfoil. For the above optimization problem, we employ Jaya Algorithm (JA) [35], which is a recently proposed gradient-free, yet effective, meta-heuristic optimization technique, whose excellent performance has been showcased in various engineering applications by different scholars [36,37,38,39]. Finally, the lift coefficient is evaluated with the help of the well-known XFOIL4 software package, which is a popular tool used in assessing the performance of subsonic airfoils and hydrofoils; see [29,30].

2.5.1. Jaya Algorithm (JA)

The Jaya algorithm is a widely used population-based optimization technique employed for solving optimization problems. It draws inspiration from social behavior, specifically the concept of individuals within a population enhancing their performance by learning from one another. Initially, a population of potential solutions is randomly generated, with each solution represented as a vector of design variables. Subsequently, the objective function is evaluated for each solution, and the best solution in the population is identified. During each iteration of the Jaya algorithm, the population is updated by considering the current best solution and improving the remaining solutions. This involves adjusting the design variables of each population member towards the current best solution.

The adjustment process in the Jaya algorithm encompasses two fundamental steps: exploration and exploitation. In the exploration phase, each solution within the population undergoes random modifications within a specified range. This enables the algorithm to explore the design space and discover promising new regions. In the exploitation phase, the design variables of each solution are updated based on the current best solution. This facilitates the exploitation of already identified regions that yield improved objective values. The adjustment process is guided by control parameters, such as learning rates and exploration factors, which determine the extent of exploration and exploitation.

The Jaya algorithm iterates until a predefined stopping criterion is met, such as reaching a maximum number of iterations or achieving a desired solution quality level. At the conclusion of the algorithm, the best solution obtained throughout the iterations is considered the optimal solution to the given optimization problem. The Jaya algorithm has proven successful in solving various optimization problems, including engineering design, scheduling, and parameter tuning. Its simplicity and efficacy have made it a popular choice for addressing real-world challenges; see again [36,37,38,39].

3. Results and Discussion

In this section we present the constructed latent design spaces for different SSVs and shape discretizations, and subsequently compare them with respect to the design space quality metrics introduced in Section 2.4.1. Finally, we demonstrate the use of these latent spaces in performance-based shape optimization, and assess their appropriateness and capacity to reconstruct and/or improve optimal designs.

3.1. Geometric Moment Invariants in SSVs

Although it is possible to calculate geometric moments invariants of any order using the formulation presented in Section 2.2, high-order geometric moments are susceptible to noise [40] and may be therefore easily affected by numerical inaccuracies [32] in their calculation. Moreover, a review of the literature across different fields, such as kinetic equations [41] and shape retrieval [33], indicates that moments of order higher than four are rarely beneficial. We have therefore limited our study to geometric moment invariants of up to fourth order. Table 1 includes all moment invariants, up to fourth order for the NACA 2410 airfoil profile, which will be considered as the baseline design in the shape optimization example presented herein.

One may easily observe from Table 1 that when considering moment invariants, only second order moments and above can be actually used in the combined geometry and moment vector, $\mathit{p}\left({\mathbf{v}}_{\overline{\mathit{\theta }}},\mathit{\mu }\left({\mathbf{v}}_{\overline{\mathit{\theta }}}\right)\right)$ or $\mathit{p}\left(\mathit{\vartheta }\right)$ which was introduced in Section 2.3, since ${\mu }_{\left(0\right)}=1$ and ${\mu }_{\left(1\right)}=(0,0)$ for all designs.

3.2. Construction of Latent Spaces

To evaluate and analyze the effectiveness of the proposed method, latent spaces with varying design/shape signature vectors, SSVs, and discretization type are considered. The combined geometry/moments SSVs considered in this study for each shape discretization along with the corresponding latent space notation are presented in Table 2:

The efficiency of the generated latent spaces is measured by the achieved dimensionality reduction whereas their quality is assessed using the criteria discussed in Section 2.4.1. Specifically, we begin by evaluating the achieved space-dimension reduction, depicted in Figure 3 and Figure 4. We then measure the diversity and robustness of these spaces (see Figure 5 and Figure 6, respectively), followed by their ability to reconstruct, or even improve on optimal designs acquired from the initial design space; see Table 3 and Figure 7 and Figure 8.

With regards to efficiency, we count the number of KL-modes needed to retain a given percentage of variance from the original space. Specifically, in Figure 3 the percentage of variance achieved for each discretization and SSV against the number of employed KL-modes is plotted. If we set the variance threshold to 95% and count the number of modes required to achieve this level of variance for each SSV and discretization type, we get the results depicted in Figure 4. As one may observe in Figure 4, we achieve a significant dimensional reduction of the original space, i.e., from 17 dimensions5 down to 4 or 5 dimensions while retaining 95% of variance exhibited in the original space. Furthermore, the inclusion of moments in the SSV has a slight beneficial effect in all discretizations studied, but the last one, i.e., uniform spacing at the physical domain. However, if we closely examine Figure 3d, we see that although the inclusion of moments for this discretization type negatively affects the achieved variance, the difference is almost negligible and therefore by slightly modifying the percentage threshold, we could easily claim that the dimension of the latent space for this case is equal to four regardless of the inclusion of moments or not. A similar picture can be drawn for the remaining discretization examined herein, with a slight positive impact to dimensionality reduction when second, third or fourth moments are included in the shape signature vector. However, since the effect is not very pronounced, we may in summary claim that the best performing discretization type is the uniform spacing at the physical domain, closely followed by the remaining three discretization types with a slight positive impact, for all cases, when moments are included in the SSV. However, this comparison does not reveal the full picture and we need to look into the quality of each latent space and its capacity to be used in design exploration and shape optimization. We subsequently examine these two aspects in Section 3.3 and Section 3.4, respectively.

3.3. Latent Space Quality Analysis

The analysis conducted here aims to determine whether the resulting subspaces are suitable for design exploration and optimization. To this end, we assess their ability to adequately capture the underlying shape structure using the latent parameter vector $\mathbf{u}$ and to generate valid and diverse geometries. The diversity metric is extremely important for design exploration and the capacity of the latent space to produce novel designs as it measures the difference between shapes included the latent space, as described in Section 2.4.1. Robustness, or shape validity, measures the percentage of invalid shapes residing in the latent space, and, as mentioned before, we ideally aim for latent spaces that eliminate, or at least minimize, the percentage of invalid designs. Needless to say that both quality metrics are important for shape optimization as both the existence of diverse designs as well as the absence of invalid designs have a positive impact on design optimization. For assessing these quality metrics, we perform five random Monte Carlo samplings with 10000 parameter vectors from each latent space. Subsequently, we estimate the Robustness with the average percentage of shapes with invalid geometries appearing in each sampling, whereas diversity is estimated by measuring the average distance between designs as described in Section 2.4.1. The acquired estimates for all discretization type and employed SSVs are presented in Figure 5 for diversity and in Figure 6 for robustness.

The estimated maximum diversity for each examined latent space, measured via Eq. (17), is shown in Figure 5. In all cases the estimated values correspond to rich latent spaces as they clearly surpass the mean value of $2.6$ corresponding to the original space, and approach or even surpass its maximum diversity value of $5.01$. Therefore, we can state that the examined latent spaces are sufficiently rich, with diverse and novel designs, regardless of the discretization type or employed SSV. However, some slight differences among discretization types can be observed. Specifically, the curvature-based discretization produces a slightly less diverse latent space and this result can be justified by the fact that this discretization exhibits a concentration of points in the area near the leading edge with a small number of points distributed along the remaining foil profile; see Figure 2c. This effectively limits the exhibited shape diversity in the area of the leading edge. The remaining discretization types have a more balanced distribution of points along the profile which benefits shape diversity as is demonstrated in Figure 5. Finally, as can be easily observed in the same figure, the inclusion of moments in the shape signature vector has a negligible effect on the estimated diversity.

Although a rather balanced performance was observed with respect to diversity in constructed latent spaces, the picture changes significantly when we examine robustness. As shown in Figure 6, both discretization type and employed SSVs play a significant role in the existence of invalid shapes in the constructed latent spaces. Specifically, the ${\mathcal{U}}_{(-1)}$ latent space, constructed using only geometry as in [9], includes a percentage of invalid shapes which ranges from 4.87% for cosine spacing, down to 2.56% when a uniform arc-length discretization is employed. For the first three discretizations, i.e., uniform (at the parametric domain), cosine and curvature-based spacing, robustness improves significantly when moments are included in the SSV, with invalid shapes being practically eliminated for the first two discretizations when fourth moments are included; see Figure 6. This behavior is however not observed for the last discretization type when second moments are included in the SSV. Specifically and in this case, the percentage of invalid shapes remains practically constant across varying SSVs with the exception of an exhibited increase when second moments, either separately or in combination, are included in the SSV. This is a rather unexpected result which requires further investigation. With the exception of the behavior exhibited for the last discretization type, our results for both diversity and robustness are in alignment with the results reported in [3] for the case of similar latent design spaces for ship hulls.

3.4. Latent Spaces in Shape Optimization

The last set of results pertains to the examination of latent spaces in performance-based shape optimization and assessment of their capacity to reconstruct, or even improve, optimal designs acquired with the initial design space which was produced by the employed parametric model presented in Section 2.1.1. The shape optimization problem considered here is the lift optimization problem subject to a foil area constraint formulated in Eq. (18). Apart from the area constraint appearing in Eq. (18), the side constraints in the same equation are implemented using Eq. (16). Finally, the employed optimization method is the Jaya algorithm discussed in Section 2.5.1.

As mentioned above, we perform a series of shape optimizations aiming at lift maximization under the area constraint that corresponds to NACA 2410 profile with the obtained results being collected in Table 3. The lift maximization problem, when the initial design space is used, results in a reference value of the lift coefficient which equals to 2.0518, as can be seen in the same table. All latent spaces lacking moments in their SSVs fail to attain this value and produce foil profiles with a lift coefficient value ranging between 1.12 and 1.42. Three optimized foil profiles for each discretization are depicted in Figure 7 which yield the average values shown in Table 3. The reference NACA 2410 foil profile is also drawn, with a dashed black line, on the same figure for reference purposes. The results for all discretizations show a rather large dispersion with no clear tendency and inferior performances when compared to the reference optimized design. However, this behavior changes significantly when moments are included in the SSVs. As we can observe in Table 3, the lift maximization problem results in improved values when moments are added in all cases, with the latent spaces based on uniform parametric spacing and curvature-based discretizations producing results that are superior to the outcome achieved with the initial design space; see rows 1 and 3 in the same table. The latent space based on cosine spacing exhibits improvement with the inclusion of moments and we also manage to reproduce the optimal design when fourth moments are included. However, for this case the results are generally inferior when compared to the previous two discretizations. The remaining discretization, i.e., uniform spacing at the physical domain, exhibits also some general improvement when moments are used, but does not manage to reconstruct or improve the initial optimum design. The SSVs including moments exhibit good or superior6 results across all discretizations and we have indicatively included here the optimized foils shapes for ${\mathcal{U}}_{\left(4\right)}$ in Figure 8. As one may easily observe in the same figure, the best performing discretizations, i.e., curvature-based, cosine, and uniform parametric spacing, tend to produce a very similar foil shape (see blue, red, and yellow-colored curves in Figure 8) which outperforms the reference design, while the fourth discretization fails to follow the trend exhibited by the remaining latent spaces. In summary, we can state that the curvature-based discretization, which adaptively distributes points along the foil profile, exhibits the best performance, closely followed by the uniform discretization in the parametric domain. The performance of the latter discretization is justified by the fact that shape information is indirectly incorporated in the discretization since the employed parametrization is coming from the parametric model. The cosine-spacing discretization shows rather inferior but acceptable performance, with the exception of the case where fourth moments are included in the SSV, in which optimimal behavior is also achieved. We should note here that the cosine and uniform spacing (in the physical domain) discretizations are shape-agnostic in the sense that they do not incorporate local shape features, but rely on a more general rule-set for point discretization, i.e., dense concentration of points at leading and trailing edge for cosine spacing and equal arc-lengths for the uniform parametrization at the physical domain. Therefore their inferior performance can be attributed to this feature. The worst performer is clearly the last discretization (uniform spacing in the physical domain) which fails to reproduce the optimal design across all examined spaces, although a slight improvement is observed when moments are included in its SSV.

As a final note regarding shape optimization using the various latent spaces, we need to briefly discuss the challenges faced which were mainly related to handling invalid shapes in some latent spaces; see also Figure 6. A series of additional filters and constraints were implemented, when needed, to address these problems and avoid the appearance of invalid shapes in the shape optimization procedure. Specifically, the specific problems and measures taken are collectively described below:

• The abscissas of the first and last foil profile point were not the same, for some designs, leading to performance inconsistencies when lift coefficient was assessed using XFOIL. In most cases, the problem was a minor difference in the coordinate value which could be easily fixed by equating the abscissas of the two points. In the rare occasion of a big difference, the design was penalized.
• A similar problem was faced with the ordinates of the same pair of points which was handled by penalizing design exhibiting blunt trailing edges with a gap exceeding the maximum gap found in the UIUC database.
• Shape irregularities including self intersections, abrupt shape changes, and lack of abscissa monotonicity along the upper or lower foil side. In all cases, these shape irregularities were monitored and penalized by the employed objective function.

4. Discussion and Conclusions

In this work, we have studied the effects of shape discretization and geometric moments’ inclusion on the performance and quality of latent design spaces produced by the generalized Karhunen-Loève expansion in the context of design space dimensionality reduction. Our results demonstrate a significant reduction in the dimensionality of the original design space while maintaining a high representational capacity and a large percentage of valid geometries that facilitate fast convergence to optimal solutions in performance-based optimization. We have also demonstrated that both the type of discretization and the composition of the shape signature vector significantly affect the resulting latent spaces and their capacity to be used in design exploration and shape optimization. In summary, we may claim that shape-informed discretizations augmented with geometric moments can produce rich latent spaces which exhibit robustness and can effectively substitute the original design space when performing shape optimization. This leads to equivalent, if not better, optimal designs with a significant reduction in the computational cost as the resulting latent spaces reduce the problem’s dimensionality by 60% or more.

The methodological approach and knowledge gained in this study can be expanded and employed in constructing discretizations and latent spaces that would address the problem of wings’ and/or blades’ design optimization, which obviously constitutes the next logical step in this line of work.