# Data-driven aerodynamic shape design with distributionally robust
optimization approaches

Long Chena, Jan Rottmayera, Lisa Kuscha, Nicolas R. Gaugera and Yinyu Yeb L.
Chen. Email<EMAIL_ADDRESS>Rottmayer. Email:
<EMAIL_ADDRESS>Kusch. Email<EMAIL_ADDRESS>R.
Gauger. Email<EMAIL_ADDRESS>Ye. Email<EMAIL_ADDRESS>aChair for Scientific Computing, University of Kaiserslautern-Landau (RPTU),
Germany;
bDepartment of Management Science and Engineering and ICME, Stanford
University, USA

###### Abstract

We formulate and solve data-driven aerodynamic shape design problems with
distributionally robust optimization (DRO) approaches. Building on the
findings of the work [13], we study the connections between a class of DRO and
the Taguchi method in the context of robust design optimization. Our
preliminary computational experiments on aerodynamic shape optimization in
transonic turbulent flow show promising design results.

###### keywords:

Aerodynamic shape optimization; Distributionally robust optimization; Robust
design optimization and Taguchi method; Stochastic gradient methods

##  1 Introduction

Optimization under uncertainties remains one of the major challenges for
aerodynamic design optimization [25]. Uncertainty has the potential to render
an optimal shape design worthless, even if obtained using sophisticated
numerical approaches, as their conclusions are not realized in practice due to
inevitable variations in problem data [39]. On the other hand, the increasing
available data allows an unprecedented insight into uncertainties. For
example, it is nowadays not difficult to acquire data about various key flight
conditions (Mach number, altitude, temperature, etc.) [3], but the
incorporation of these data in design optimization is challenging. This paper
presents a work towards data-driven design optimization using distributionally
robust optimization (DRO) approaches [7] with applications in aerodynamic
shape design.

Mach
---
Lift
---
0.50
---
0.53
---
0.47
---
0.72
---
0.74
---
0.70
---
$w_{1}=0.2$
---
$w_{3}=0.4$
---
$w_{5}=0.2$
---
$w_{4}=0.1$
---
$w_{2}=0.1$
---
Figure 1: Multi-point aerodynamic shape optimization.

###  1.1 Motivation

To account for the variability in the flight conditions [8], a multi-point
optimization with the weighted sum formulation was introduced in [34], and is
considered a standard approach in academia and industry. Figure 1 illustrates
such a multi-point aerodynamic shape design optimization problem, in which the
objective function is defined as a composite function through a weighted sum
of five different flight conditions in Mach numbers and lift coefficients,

$f(x)=\sum_{i=1}^{5}\omega_{i}c_{d}\left(x;{Ma}_{i},{c_{l}}_{i}\right),$ (1)

where $c_{d}$ is the drag coefficient, $Ma$ is the mach number, $c_{l}$ is the
target lift coefficient, $\omega$ is the weight factor, and $x$ is the shape
design variable. The main shortcoming of (1) is the selection of appropriate
flight points and their associated weights [28].

Mach
---
0.66
---
0.68
---
0.70
---
0.72
---
0.74
---
0.76
---
0.78
---
Lift
---
0.56
---
0.54
---
0.52
---
0.50
---
0.48
---
0.46
---
0.44
---
Figure 2: Flight points during cruising of an A320 from Helsinki to Munich.

Figure 2 presents data points of a real flight of an A320 from Helsinki to
Munich111$Ma$ is computed using real-world data of altitude and ground speed
(see https://www.grc.nasa.gov/www/k-12/airplane/mach.html). As we do not have
the data for the lift force, $c_{l}$ is synthesized to be drawn from a normal
distribution $N(0.5,0.01)$ for demonstration purposes.. Following the multi-
point formulation, one way to incorporate these data into design optimization
is to consider the following weighted sum objective function,

$f(x)=\sum_{i=1}^{N}\frac{1}{N}c_{d}(x;{Ma}_{i},{c_{l}}_{i}),$ (2)

where $N$ is the total number of data points, and
$\omega_{i}=\frac{1}{N},i=1,2,...,N$. By assigning equal weight of probability
to all observations, (2) becomes an empirical mean optimization problem.

If the number of data points is large, (2) would be a more realistic objective
function than (1) almost surely. On the other hand, one should realize that
the empirical distribution is itself an approximation to the true
distribution. For example, it is subject to shift due to measurement errors
and varying conditions, such as different weather and routes. Furthermore, the
design process of an aircraft (wing) precedes its mass production, and
empirical data gathered (for example, from previous designs) cannot be the
true data of the to-be-optimized design. Since an engineering solution should
always be on the safe side, we want a design optimization framework that is
robust to distributional changes. Such a framework is called distributionally
robust optimization (DRO) [7]. Instead of minimizing the empirical expected
value, DRO seeks an optimal design against the worst-case expected value
evaluated on a set $\mathdutchcal{D}$ of candidate distributions
$Q=\left[q_{1},q_{2},...,q_{N}\right]$, assuming that $\mathdutchcal{D}$
contains the true distribution $P_{true}$,

$f(x)=\max_{Q\in\mathdutchcal{D}}\sum_{i=1}^{N}q_{i}c_{d}(x;{Ma}_{i},{c_{l}}_{i}).$
(3)

This work aims to formulate, study, and solve data-driven aerodynamic shape
design problems based on (3). We summarize our main contributions as follows.

###  1.2 Main contributions

  * 1.

We present data-driven design optimization frameworks for aerodynamic shape
design with DRO.

  * 2.

We study the connections and differences between DRO and the classical robust
design method of Taguchi.

  * 3.

We present stochastic algorithms for solving the proposed aerodynamic shape
design optimization problems that include additional inequality constraints.

  * 4.

We report preliminary results on transonic airfoil shape design by applying
our proposed optimization formulations and algorithms.

##  2 Related work

###  2.1 Aerodynamic shape optimization under uncertainties

Probabilistic airfoil design dates back to Huyse [15], who used a formulation
based on minimizing the expected value of the drag coefficient. The expected
drag coefficient was optimized under a uniformly distributed uncertain Mach
number in a transonic Euler flow. Different studies for probabilistic airfoil
design have been conducted afterwards. In [14], a design optimization is
performed using the simplified assumption of independent normally distributed
geometrical parameters of a 3D wing. In [30] geometrical and operational
uncertainties are considered separately for a 2D airfoil in Euler flow using a
worst-case formulation. Lee et al. [22] optimize a three-dimensional wing
under uncertainties in the Mach number. The sampled mean and the variance of
the drag coefficient are considered as objective functions for a multi-
objective evolutionary algorithm. In [38] different probabilistic design
formulations are compared for a lift-constrained problem in a 2D Euler flow
with uncertainties in the operating conditions. The Mach number and angle of
attack are distributed using a truncated normal distribution.

In the present work we consider stochastic optimization methods for
aerodynamic shape optimization under uncertainties. This approach has also
been pursued recently by other authors. In [11], the authors present a
stochastic optimization strategy for PDE-constrained optimization. In the
context of aerodynamic shape optimization, stochastic gradient descent has
been applied for the robust design optimization of the NACA-0012 airfoil in
[16], which considers uncertainties in operating conditions and model
uncertainty.

###  2.2 Distributionally robust optimization

Stochastic optimization is a classical framework that allows to model
uncertainty within a decision-making framework. Stochastic optimization
assumes that the decision maker has complete knowledge about the underlying
uncertainty through a known probability distribution and minimizes a
functional of the cost. Distributionally robust optimization (DRO)
acknowledges that we only have partial knowledge of the statistical properties
of the uncertain parameters. The core idea of DRO goes back to Scarf [37], who
proposed to define a set $\mathdutchcal{D}$ of probability distributions that
is assumed to include the true distribution $P_{true}$, and reformulated the
objective function with respect to the worst-case expected cost over the
choice of a distribution in this set. His idea was first generalized and
formally called DRO in [7]. This min-max DRO approach represents a systematic
modeling framework that bridges the gap between data and decision-making. In
the past decade, significant research has been done on DRO within the
operations research and statistical learning communities (see, e.g., [1, 10,
23, 44]), and we refer a comprehensive review to [33]. In engineering design,
DRO is a relatively new framework, with [6][17] being the two recent
literature. In this work, we formulate data-driven aerodynamic shape design
problems using DRO frameworks and present practical algorithms to solve them.

###  2.3 Robust design and robust design optimization

The Taguchi method lays the foundation of robust design and is considered the
pioneering work in bringing statistics into engineering design [2][42]. A
fundamental contribution of Taguchi was “his observation that one should
design a product in such a way as to make its performance insensitive to
variations in variables beyond the designer’s control” [40]. Prior to the era
of computational simulation and optimization, Taguchi method utilizes the
design of experiment for parameter design [5][26]. With the advances in high
performance computing, numerical analysis and optimization, a new field called
robust design optimization (RDO) has emerged. In RDO, the mean-variance
optimization (MV) is one of the most direct and fundamental mathematical
optimization formulations222While mean-variance optimization is often referred
to in the literature as the standard approach for robust design optimization,
we would like to keep this statement open, because the influence of the
Taguchi method is far-reaching [41]. Instead, we would interpret MV as a
representative mathematical optimization formulation of the Taguchi method.
following Taguchi’s idea of not only optimizing the mean performance but to
simultaneously reducing the variance [43],

$\min_{x}\mathbb{E}_{\xi\sim P}\left[f(x;\xi)\right]+\mu\mathbb{V}_{\xi\sim
P}\left[f(x;\xi)\right].$ (MV)

In this work, building on the findings of the work [13] showing that MV is
equivalent to $\phi$-divergence based DRO when considering a “small amount of
robustness”, we further study the connections and differences between MV and
DRO in the robust design optimization context. In particular, we show that
under certain conditions, the $\chi^{2}$-divergence based DRO can be viewed as
a dual optimization formulation of the Taguchi method.

##  3 Data-driven design optimization modeling frameworks

Aerodynamic shape optimization is a typical PDE-constrained optimization
problem, which is commonly described as

$\begin{split}&\min_{u,x}~{}~{}~{}F(x,u),~{}~{}~{}~{}~{}~{}~{}~{}F:X\times
U\rightarrow\mathbb{R}\\\ &~{}\text{s.t.}~{}~{}~{}~{}H(x,u)=0,~{}~{}H:X\times
U\rightarrow Z\end{split}$ (4)

where x is the optimization variable, $u$ is the state variable, $X,U,Z$ are
vector spaces, $F$ is the objective function, $H$ is the PDE constraint, and
$F,H$ are sufficiently smooth functions. By implicit function theorem, if the
Jacobian matrix of H about $u$ at a point $H(u_{p},x_{p})=0,(u_{p},x_{p})\in
X\times U$ is invertible, then there exists a unique continuously
differentiable vector-valued function $\psi$ such that in a neighborhood of
$(u_{p},x_{p})$ that satisfies

$\begin{split}\psi(x_{p})&=u_{p},\\\ H(x,\psi(x))&=0.\end{split}$ (5)

Suppose that $\psi$ exists everywhere in the domain of the function $H$, then
we can reformulate the PDE-constrained optimization (4) as an unconstrained
optimization of solely the optimization variables x,

$\min_{x}~{}~{}~{}f(x):=F(x,\psi(x)).$ (6)

By (6), in this work, we refer to a PDE-constrained optimization in the form
(4) as an unconstrained optimization problem. For practical aerodynamic shape
design, additional engineering constraints need to be considered. In a
deterministic setting, we therefore consider a constrained optimization
problem written as

$\begin{split}&\min_{x}~{}~{}~{}f(x)\\\
&~{}\text{s.t.}~{}~{}~{}~{}g_{i}(x)\leq 0,~{}i=1,2,...,m\end{split}$ ($DET$)

In the following, we present several data-driven optimization frameworks based
on ($DET$).

###  3.1 Empirical mean optimization for aerodynamic shape design

We start with the formulation of empirical mean optimization (EMO) for
aerodynamic shape design. Throughout this paper, we only consider the case
where the data in the objective function is uncertain. As the PDE constraint
is reformulated in the objective function in (6), the considered model is
general enough for many practical applications.

Suppose we want to minimize the objective function $f(x;\xi)$ where
$x\in\mathbb{R}^{d}$ is the (shape) design variable, and
$\xi\in\mathbb{R}^{r}$ is a random variable. Assume we have observed N
independent samples of $\xi_{i},i=1,2,...,N$, we want to minimize the
empirical mean subject to deterministic constraints,

$\begin{split}&\min_{x}~{}~{}~{}\mathbb{E}_{\xi\sim\hat{P}_{N}}\left[f(x;\xi)\right]=\sum_{i=1}^{N}\frac{1}{N}f(x;\xi_{i})\\\
&~{}\text{s.t.}~{}~{}~{}~{}g_{i}(x)\leq 0,~{}i=1,2,...,m,\end{split}$ ($EMO$)

where $\hat{P}_{N}$ denotes the empirical distribution. By the law of large
numbers, the empirical mean converges to the expected value as
$N\rightarrow\infty$. Compared to the multi-point optimization (1), we may
interpret ($EMO$) as a multi-point optimization with a large number of points.
While it is a heuristic choice of the weighing factors $\omega_{i}$ and data
points $\xi_{i}$ for multi-point objective functions, ($EMO$) minimizes the
empirical risk using the known observations of $\xi$ —a fundamental principle
in statistical learning. On the other hand, the statistically more meaningful
optimization modeling framework of ($EMO$) comes with the cost that the
objective function is much more computationally expensive to evaluate due to
the large $N$. To tackle the computational challenge, we rely on stochastic
optimization approaches (see section 5), which have proven to be efficient and
effective in countless machine learning applications.

###  3.2 Basics of distributionally robust optimization

In contrast to the empirical mean optimizations, for DRO it is desired that
the optimized design still has good performance under distribution shift.
Specifically, DRO proposes to minimize the worst-case expected performance
over a set $\mathdutchcal{D}$ of probability distributions $Q$, assuming that
$\mathdutchcal{D}$ contains the true distribution $P_{true}$,

$\begin{split}&\min_{x}\Psi(x):=\max_{Q\in\mathdutchcal{D}(P,\rho)}\mathbb{E}_{\xi\sim
Q}\left[f(x;\xi)\right],\\\ &~{}\textnormal{s.t.}~{}~{}g_{i}(x)\leq
0,~{}~{}i=1,2,...,m.\end{split}$ ($DRO_{con}$)

The set $\mathdutchcal{D}(P,\rho)$ is called the ambiguity set and is
typically defined as

$\mathdutchcal{D}(P,\rho):=\\{Q:d(Q,P)\leq\rho\\},$ (7)

where $d$ measures the distance between two probability distributions, $P$ is
the nominal distribution, and $\rho\geq 0$ corresponds to the size (radius) of
the ambiguity set.

Another well-established DRO formulation uses a soft penalty term instead of
imposing a hard constrained ambiguity set, resulting in the penalized problem:

$\begin{split}&\min_{x}\Psi(x):=\max_{Q}\\{\mathbb{E}_{\xi\sim
Q}\left[f(x;\xi)\right]-\delta d(Q,P)\\},\\\
&~{}\textnormal{s.t.}~{}~{}g_{i}(x)\leq 0,~{}~{}i=1,2,...,m,\end{split}$
($DRO_{pen}$)

where $\delta>0$ is the penalty coefficient. By duality, the larger $\delta$,
the less robustness is considered.

In the data-driven context, the nominal distribution is constructed from the
observed data $\xi_{i},i=1,2,...,N$, i.e., we set $P=\hat{P}_{N}$ in
($DRO_{con}$) and ($DRO_{pen}$).

###  3.3 DRO frameworks based on $\phi$-divergence

Obviously, different choice of the distance measure $d$ results in different
DRO solutions. In this work, we focus on one class of such a statistical
distance, the $\phi$-divergence [1], due to its deep connections to the
Taguchi method of robust design (see section 4).

We consider $\phi$-divergence between distributions $Q$ and $P$,

$d_{\phi}(Q,P):=\int\phi\left(\frac{dQ}{dP}\right)dP.$ (8)

Suppose that the distributions $P$ and $Q$ are represented with the
probability vectors, $\mathbf{p}=(p_{1},...,p_{n})^{T}$ and
$\mathbf{q}=(q_{1},...,q_{n})^{T}$. The $\phi$-divergence of $\mathbf{q}$ from
$\mathbf{p}$ is defined as

$d_{\phi}(\mathbf{q},\mathbf{p})=\sum_{i=1}^{n}p_{i}\phi(\frac{q_{i}}{p_{i}}).$
(9)

Different choices for $\phi$ have been proposed in literature (see, e.g.,
[31]). In this work, we focus on the $\chi^{2}$-divergence and Kullback-
Leibler (KL)/Burg entropy-based divergence and leave other useful divergences
for future investigations. The function $\phi(t)$ for (modified)
$\chi^{2}$-divergence reads

$\phi_{\chi^{2}}(t)=\frac{1}{2}(t-1)^{2}.$ (10)

For KL divergence,

$\phi_{kl}(t)=t\log t-t+1.$ (11)

The adjoint of KL divergence is Burg entropy-based divergence,

$\phi_{b}(t)=-\log t+t-1.$ (12)

Figure 3 illustrates different divergence balls with different radii. Note, KL
divergence and Burg entropy are “locally $\chi^{2}$” [32], i.e., they can be
well-approximated with $\chi^{2}$-divergence if the radius is small.

Figure 3: Ambiguity sets of different $\phi$-divergence in a barycentric
coordinate system.

###  3.4 Statistically principled determination of $\mathdutchcal{D}(P,\rho)$

DRO is a framework that naturally combines optimization with statistics. This
allows for statistically principled determinations of the ambiguity set
(distance measure and radius). We refer to [33, section 6] for an overview. In
the following, we give an example of the Burg entropy-based DRO. Recall,

$\phi_{b}(t)=-\log t+t-1.$

Consider the Burg entropy constrained ambiguity set with a uniform nominal
distribution $\hat{P}_{N}$,

$\mathdutchcal{D}(P,\rho):=\biggl{\\{}\mathbf{q}:\sum_{i=1}^{N}\frac{1}{N}\log(\frac{1}{Nq_{i}})\leq\rho,\mathbf{q}\in\Delta^{N}\biggl{\\}},$
(13)

Where $\Delta^{N}$ denotes the probability simplex. The inequality in the
above formula is equivalent to

$\sum_{i=1}^{N}\log(q_{i})\geq\sum_{i=1}^{N}\log(\frac{1}{N})-N\rho.$ (14)

(13) is closely related to likelihood robust optimization [44], which makes
assumptions that random variables taking values in set
$\Xi=\\{\xi_{1},\xi_{2},...,\xi_{\tilde{n}}\\}$ and the set $\Xi$ is finite
and known in advance. Assuming that N independent samples $\xi$ were observed,
with $N_{i}$ occurrences of $\xi_{i}$. Then (13) rewrites as

$\mathbb{D}(\gamma)=\biggl{\\{}\tilde{\mathbf{q}}=(\tilde{q}_{1},\tilde{q}_{2},...,\tilde{q}_{\tilde{n}}):\sum_{i=1}^{\tilde{n}}N_{i}\log\tilde{q}_{i}\geq\gamma,\tilde{\mathbf{q}}\in\Delta^{\tilde{n}}\biggr{\\}},$
(15)

where $\gamma=\sum_{i=1}^{\tilde{n}}N_{i}\log\frac{N_{i}}{N}-N\rho$, and
$\mathbb{D}(\gamma)$ is called the likelihood robust distribution set. Notice
that

$\sum_{i=1}^{\tilde{n}}N_{i}\log(\tilde{q}_{i})=\log(\Pi_{i=1}^{\tilde{n}}\tilde{q}_{i}^{N_{i}}),$
(16)

(15) implies that $\mathbb{D}(\gamma)$ contains all the distributions with
support in $\Xi$ such that the observed data achieve an empirical likelihood
of at least exp($\gamma$). With the set $\mathbb{D}(\gamma)$, the likelihood
robust optimization writes

$\min_{x}\Psi(x):=\max_{\tilde{\mathbf{q}}\in\mathbb{D}(\gamma)}\sum_{i=1}^{\tilde{n}}\tilde{q}_{i}f(x;\xi_{i}).$
(17)

With Bayesian statistics and empirical likelihood theory, it was shown in [44]
that one can approximately choose the threshold $\gamma$ as

$\gamma^{*}=\sum_{i=1}^{\tilde{n}}N_{i}\log\frac{N_{i}}{N}-\frac{1}{2}\chi^{2}_{{\tilde{n}}-1,1-\alpha},$
(18)

where $\chi^{2}_{{\tilde{n}}-1,1-\alpha}$ is the $1-\alpha$ quantile of a
$\chi^{2}$ distribution with ${\tilde{n}}-1$ degrees of freedom. We conclude
this section by remarking that there exist many other DRO frameworks that use
different statistical distance measures for the ambiguity set, and that their
applications to engineering design are worth investigating.

##  4 Connections and differences between DRO and Taguchi method

In this section, we first review recent results in the literature on the
relationship between DRO and mean-variance optimization. We then show that,
under certain conditions, DRO can be considered a dual optimization
formulation of the Taguchi method. We also discuss the differences between the
two methods when these conditions are not satisfied in order to give a
complete comparison within the scope of robust design optimization.

###  4.1 The equivalence of DRO and mean-variance optimization by small
ambiguity set

Let’s first consider the penalized $\phi$-divergence DRO objective with
nominal empirical distribution $\hat{P}_{N}$,

$\min_{x}\Psi(x;\hat{P}_{N}):=\max_{Q}\\{\mathbb{E}_{\xi\sim
Q}\left[f(x;\xi)\right]-\delta d_{\phi}(Q,\hat{P}_{N})\\}.$ (19)

It was shown in [13] that, if $\delta$ is sufficiently large and $\phi$
satisfies the conditions,

  * 1.

$\phi:\mathbb{R}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ is a closed, convex,
function,

  * 2.

$\phi(z)\geq\phi(1)=0$ for all $z$, is twice continuously differentiable
around $z=1$ with $\phi^{\prime\prime}(1)>0$,

then, the penalized DRO is asymptotically equivalent to the mean-variance
optimization up to $o(1/\delta)$,

$\min_{x}\Psi(x;\hat{P}_{N})=\min_{x}\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]+\frac{1}{2\delta}\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]+o(1/\delta).$
(20)

Notice that the $\chi^{2}$, KL, and Burg entropy-based divergence considered
in this work all satisfy these conditions. Hence, for large enough $\delta$,
the penalty-based $\phi$-divergence DRO (19) well-approximates the mean-
variance optimization, independent of the choice of the divergence measure.

Similar results to (20) are shown for the DRO objective with hard constraint
in the works in [21, 27],

$\min_{x}\max_{Q\in\mathdutchcal{D}(\hat{P}_{N},\rho)}\mathbb{E}_{\xi\sim
Q}\left[f(x;\xi)\right]=\min_{x}\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]+\sqrt{2\rho\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]}+o(\rho).$
(21)

With the DRO framework, we gain an important insight into the fundamental
robust design framework (MV) when $\mu$ is small: Adding variance/standard
deviation to the expectation in the objective function can be seen as a robust
optimization over certain distributions shifted from the nominal distribution.
In the following, we provide more insight to the Taguchi method beyond the
regime of “small amount of robustness” by studying $\chi^{2}$-divergence DRO
and its dual formulation.

###  4.2 $\chi^{2}$-divergence DRO: a dual view on Taguchi method

We show non-asymptotic result on the equivalence with a particular choice of
the divergence measure, the $\chi^{2}$-divergence. We give explicit conditions
for the parameter $\delta$ and show that under these conditions, the
$\chi^{2}$-divergence DRO can be viewed as a dual mathematical optimization
formulation of the Taguchi method.

######  Proposition 4.1.

Consider $\chi^{2}$-divergence distributionally robust optimization with

$\phi_{\chi^{2}}(t)=\frac{1}{2}(t-1)^{2}.$

Suppose that $\hat{P}_{N}$ is the empirical probability distribution on a set
of independently and identically distributed data $\\{\xi_{i}\\}_{i=1}^{N}$,
according to the true distribution $P_{true}$. Then, with any penalty
$\delta>0$ that satisfies

$\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]-f(x;\xi_{k})\leq\delta,~{}\forall
k,$ (22)

it holds

$\max_{Q}\\{\mathbb{E}_{Q}\left[f(x;\xi)\right]-\delta
d_{\phi_{\chi^{2}}}(Q,\hat{P}_{N})\\}=\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]+\frac{1}{2\delta}\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right].$
(23)

###### Proof.

We first derive an equivalent formulation of the DRO objective with
$\phi$-divergence in probability vector representation (9). The key is to use
duality theory on conjugate functions [35]. The conjugate of a function
$f:\mathbb{R}\rightarrow\mathbb{R}$ can be considered its dual object, and is
defined as a function $f^{*}:\mathbb{R}\rightarrow\mathbb{R}\cup\\{\infty\\}$,

$f^{*}(s)=\max_{t\geq 0}\\{st-f(t)\\}.$ (24)

Consider finite observations with the empirical distribution $\hat{P}_{N}$,
the DRO objective function can be formulated with probability vectors as (we
shorten the notation $q_{i}\geq 0,i=1,2,...,N$ as $q\geq 0$ )

$\Psi(x):=\max_{q\geq
0}\biggl{\\{}\sum_{i=1}^{N}q_{i}f(x;\xi_{i})-\delta\sum_{i=1}^{N}p_{i}\phi(\frac{q_{i}}{p_{i}})\biggr{\\}},$
(25)

with the constraint

$h(q)=\sum_{i=1}^{N}q_{i}-1=0.$ (26)

Introducing Lagrange Multiplier $\lambda$, we can rewrite the constrained
maximization problem (25)-(26) as

$\begin{split}\Psi(x)&=\min_{\lambda\in\mathbb{R}}\max_{q\geq
0}\biggl{\\{}\sum_{i=1}^{N}q_{i}f(x;\xi_{i})-\delta\sum_{i=1}^{N}p_{i}\phi(\frac{q_{i}}{p_{i}})\biggr{\\}}-\lambda\left(\sum_{i=1}^{N}q_{i}-1\right)\\\
&=\min_{\lambda\in\mathbb{R}}\max_{q\geq
0}\sum_{i=1}^{N}\left[q_{i}f(x;\xi_{i})-\delta
p_{i}\phi(\frac{q_{i}}{p_{i}})-\lambda q_{i}\right]+\lambda.\end{split}$ (27)

Since $\hat{P}_{N}$ is the empirical distribution, we have
$p_{i}=\frac{1}{N},i=1,...,N$. By the definition of conjugate function (24),
it then follows,

$\begin{split}\Psi(x)&=\min_{\lambda\in\mathbb{R}}\max_{q\geq
0}\sum_{i=1}^{N}\left[q_{i}f(x;\xi_{i})-\delta\frac{1}{N}\phi(Nq_{i})-\lambda
q_{i}\right]+\lambda\\\ &=\min_{\lambda\in\mathbb{R}}\max_{q\geq
0}\delta\frac{1}{N}\sum_{i=1}^{N}\left[\left(\frac{f(x;\xi_{i})-\lambda}{\delta}\right)Nq_{i}-\phi(Nq_{i})\right]+\lambda\\\
&=\min_{\lambda\in\mathbb{R}}\frac{1}{N}\sum_{i=1}^{N}\left[\delta\phi^{*}\left(\frac{f(x;\xi_{i})-\lambda}{\delta}\right)+\lambda\right].\end{split}$
(28)

Thus, the equivalent dual DRO objective with penalty formulation writes

$\Psi(x)=\min_{\lambda\in\mathbb{R}}\mathbb{E}_{\hat{P}_{N}}\left[\delta\phi^{*}\left(\frac{f(x;\xi)-\lambda}{\delta}\right)+\lambda\right].$
(29)

Recall the (modified) $\chi^{2}$-divergence,

$\phi_{\chi^{2}}(t)=\frac{1}{2}(t-1)^{2},$

whose conjugate function writes

$\phi_{\chi^{2}}^{*}(s)=\left\\{\begin{split}&-\frac{1}{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}s<-1,\\\
&\frac{1}{2}s^{2}+s,~{}~{}~{}~{}s\geq-1.\\\ \end{split}\right.$ (30)

Having $s_{k}=\frac{f(x;\xi_{k})-\lambda}{\delta},$ and suppose that
$s_{k}\geq-1$ holds for every $k=1,2,...,N$, then

$\begin{split}\Psi(x)&=\min_{\lambda\in\mathbb{R}}\mathbb{E}_{\hat{P}_{N}}\delta\left[\frac{1}{2}\left(\frac{f(x;\xi)-\lambda}{\delta}\right)^{2}+\left(\frac{f(x;\xi)-\lambda}{\delta}\right)\right]+\lambda\\\
&=\min_{\lambda\in\mathbb{R}}\mathbb{E}_{\hat{P}_{N}}\left[\frac{1}{2\delta}\left(f(x;\xi)-\lambda\right)^{2}+f(x;\xi)\right].\end{split}$
(31)

The minimizer $\lambda^{\star}$ for (31) is the arithmetic mean of
$f(x;\xi_{i}),i=1,...,N$, i.e.,
$\lambda^{\star}=\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]$. Thus,

$\Psi(x)=\frac{1}{2\delta}\mathbb{E}_{\hat{P}_{N}}\left[(f(x;\xi)-\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right])^{2}\right]+\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right],$
(32)

which is the mean-variance optimization (MV) objective function. Hence our
conclusion: $\chi^{2}$-divergence based DRO is equivalent to the mean-variance
optimization when
$\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]-f(x;\xi_{k})\leq\delta,k=1,2,...,n$.
∎

A similar result can be shown for DRO with hard constraint formulation using
the same procedure as in the proof of Proposition 4.1.

######  Lemma 4.2.

Consider $\chi^{2}$-divergence distributionally robust optimization with

$\phi_{\chi^{2}}(t)=\frac{1}{2}(t-1)^{2}.$

Suppose that $\hat{P}_{N}$ is the empirical probability distribution on a set
of independently and identically distributed data $\\{\xi_{i}\\}_{i=1}^{N}$,
according to the true distribution $P_{true}$. Then, with any parameter
$\rho>0$ that satisfies

$\left(\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]-f(x;\xi_{k})\right)\sqrt{2\rho}\leq\sqrt{\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]},~{}~{}\forall
k,$ (33)

it holds

$\max_{Q\in\mathdutchcal{D}_{\chi^{2}}(\hat{P}_{N},\rho)}\mathbb{E}_{\xi\sim
Q}\left[f(x;\xi)\right]=\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]+\sqrt{2\rho\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]}.$
(34)

###### Proof.

Introducing Lagrange Multiplier $\lambda,\nu$, we can rewrite the hard
constraint DRO objective as

$\begin{split}\Psi(x)&=\min_{\lambda\in\mathbb{R},\nu\geq 0}\max_{q\geq
0}\sum_{i=1}^{N}q_{i}f(x;\xi_{i})-\lambda\left(\sum_{i=1}^{N}q_{i}-1\right)-\nu\left(\sum_{i=1}^{N}\frac{1}{N}\phi(Nq_{i})-\rho\right)\\\
&=\min_{\lambda\in\mathbb{R},\nu\geq 0}\max_{q\geq
0}\frac{\nu}{N}\sum_{i=1}^{N}\left[\left(\frac{f(x;\xi_{i})-\lambda}{\nu}\right)Nq_{i}-\phi(Nq_{i})\right]+\lambda+\nu\rho\\\
&=\min_{\lambda\in\mathbb{R},\nu\geq
0}\frac{\nu}{N}\sum_{i=1}^{N}\phi^{*}\left(\frac{f(x;\xi_{i})-\lambda}{\nu}\right)+\lambda+\nu\rho.\end{split}$
(35)

Hence,

$\Psi(x)=\min_{\lambda\in\mathbb{R},\nu\geq
0}\mathbb{E}_{\hat{P}_{N}}\left[\nu\phi^{*}\left(\frac{f(x;\xi)-\lambda}{\nu}\right)+\lambda+\nu\rho\right].$
(36)

Suppose that $\left(\frac{f(x;\xi_{k})-\lambda}{\nu}\right)\geq-1$ holds for
every $k=1,2,...,N$, then with $\chi^{2}$-divergence we have

$\Psi(x)=\min_{\lambda\in\mathbb{R},\nu\geq
0}\mathbb{E}_{\hat{P}_{N}}\left[\frac{1}{2\nu}(f(x;\xi)-\lambda)^{2}+f(x;\xi)+\nu\rho\right].$
(37)

The right hand side of (37) attains the minimum at

$\begin{split}\lambda^{\star}&=\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right],\\\
\nu^{\star}&=\sqrt{\frac{\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]}{2\rho}}.\end{split}$
(38)

Finally,

$\Psi(x)=\mathbb{E}_{\hat{P}_{N}}\left[f(x;\xi)\right]+\sqrt{2\rho\mathbb{V}_{\hat{P}_{N}}\left[f(x;\xi)\right]}.$
(39)

∎

Comparing (25) -(26) with (29), the former maximization problem over
distributions in the DRO objective is the dual of the latter minimization
problem over the Lagrange multiplier $\lambda$. When using the
$\chi^{2}$-divergence, the analytical solution to (29) is the mean-variance
objective function. Therefore, under conditions (22), the
$\chi^{2}$-divergence DRO objective is the dual of the mean-variance
objective. The dual view on the Taguchi method with DRO allows new
possibilities in the modeling of robust engineering design problems:

  * 1.

If $f(x;\xi)$ is convex in $x$, then the DRO objectives are convex surrogates
for the variance/standard deviation regularized objectives, which are
generally non-convex even if $f$ is convex in $x$. Therefore, if the
underlying design optimization problem is convex (or can be equivalently
formulated as a convex problem), an equivalent DRO formulation is easier to
solve (see [27]).

  * 2.

DRO allows statistically principled and data-driven analysis and determination
of the robust parameter $\delta$ and $\rho$, which are good suggestions for
the factors in mean-variance/standard deviation formulations. Thus, it is an
immediate contribution to existing robust design frameworks and algorithms.

  * 3.

For sufficiently large robustness (i.e., $\delta$ small enough or $\rho$ large
enough), different $\phi$-divergence results in substantially different
ambiguity sets (see figure 3). In such situations, different $\phi$-divergence
DROs will generally give different design solutions. Intuitively, one may
think that different $\phi$-divergences lead to “mean”-“variance” (-“standard
deviation”) optimizations, in which the statistical moments are differently
measured.

###  4.3 Differences between DRO and mean-variance optimization

Proposition 4.1 states that if $\delta$ satisfies (22), $\chi^{2}$-divergence
DRO is the dual formulation of the mean-variance optimization. In the
following, we study the case when (22) cannot be satisfied and elaborate on
the differences between both optimization frameworks.

Suppose that $\delta$ is small enough such that
$s_{k}=\frac{f(x;\xi_{k})-\lambda^{\star}}{\delta}<-1$ for some sample
$\xi_{k}$, where $\lambda^{\star}$ is the minimizer of (29). Then,
$\phi_{\chi^{2}}^{*}(s_{k})=-\frac{1}{2}$ by the conjugate function (30). As a
result, the equivalence of DRO and MV shown in (23) no longer holds. When
computing the gradient of the DRO objective, the following differences appear
when compared to the equivalent MV objective (32):

  * 1.

the minimizer
$\lambda^{\star}\neq\mathbb{E}_{\hat{P}_{n}}\left[f(x;\xi)\right]$ in general;

  * 2.

the contribution of the gradient $\nabla f(x;\xi_{k})$ vanishes in
$\nabla\Psi(x)$.

To gain some intuition of these differences, we consider an illustrative
example with the objective function ([2, equation (28)])

$\ell(x;\xi)=\xi-(\xi-1)x^{2},$ (40)

where $0\leq x\leq 1$, and $\xi$ has the support
$\Xi=\\{-0.8,-0.4,0.0,0.4,0.8\\}$. Furthermore, we assume $\xi$ has a nominal
empirical distribution
$\hat{P}=\left[\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5}\right]$.
From figure 4, it is clear that,

$\begin{split}\mathop{\text{argmin}}_{0\leq x\leq
1}~{}\mathbb{E}_{\hat{P}}\left[\ell(x;\xi)\right]=0,\\\
\mathop{\text{argmin}}_{0\leq x\leq
1}~{}\mathbb{V}_{\hat{P}}\left[\ell(x;\xi)\right]=1.\\\ \end{split}$ (41)

Consider the mean-variance optimization for $\ell(x;\xi)$ with $\mu\geq 0$,

$\min_{0\leq x\leq
1}\mathbb{E}_{\hat{P}}\left[\ell(x;\xi)\right]+\mu\mathbb{V}_{\hat{P}}\left[\ell(x;\xi)\right].$
(42)

The solution of the above bi-objective optimization lies in the Pareto
frontier $\\{x\in\mathbb{R}:0\leq x\leq 1\\}$, depending on the choice of
$\mu$.

Figure 4: Function plots for (40) with different $\xi$.

Let’s now consider the $\phi$-divergence DRO formulation for (40),

$\min_{0\leq x\leq 1}\Psi(x)=\max_{q\geq
0}\biggl{\\{}\sum_{i=1}^{5}q_{i}\ell(x;\xi_{i})-\delta\sum_{i=1}^{5}p_{i}\phi(\frac{q_{i}}{p_{i}})\biggr{\\}}.$
(43)

Let
$q^{\star}(x)=\left[q_{1}^{\star}(x),q_{2}^{\star}(x),...,q_{5}^{\star}(x)\right]$
be the maximizer of the inner problem at $x$, we can write the gradient of the
DRO objective (43) as

$\nabla\Psi(x)=\sum_{i=1}^{5}q_{i}^{\star}(x)\nabla\ell(x;\xi_{i}),$ (44)

which is a convex combination of all sample gradients. By a gradient descent
procedure, we can see that

$\mathop{\text{argmin}}_{0\leq x\leq 1}~{}\Psi(x)=0,$ (45)

and the result hold for any $\delta\geq 0$. The result (45) makes sense,
because the solution of DRO is upper bounded by the solution of the related
robust optimization,

$\min_{0\leq x\leq 1}\max_{\xi\in\Xi}~{}~{}\xi-(\xi-1)x^{2},$ (46)

which seeks a robust solution over the supports instead of the distributions.
Obviously, $\xi=0.8$ maximizes the inner optimization problem (notice that
$0\leq x\leq 1$), and the minimizer for $\ell(x;\xi=0.8)$ is $x^{\star}=0$.

To summarize, while mean-variance optimization is a bi-objective optimization,
the min-max formulation of DRO results in solutions that are upper bounded by
the related robust optimization. As a consequence, DRO formulation provides a
degree of certainty to avoid generating overly conservative designs. On the
other hand, depending on applications, the whole mean-variance trade-off
Pareto front is worthwhile investigating.

##  5 Algorithms

The major challenges in applying the present data-driven optimization modeling
frameworks to aerodynamic shape design include 1) the linear scaling of the
number of function/gradient evaluations with the number of data points, and 2)
the high computational cost in solving even a single-point CFD computation. If
there is a surrogate model at hand that provides cheap (and accurate enough)
evaluations of $f(x;\xi)$ for any feasible $x$ and $\xi$, then one can resort
to the empirical mean formulation ($EMO$) or DRO formulations ($DRO_{con}$)
and ($DRO_{pen}$) and solve these data-driven optimization problems with well-
established deterministic optimization solvers. However, training such a
surrogate for practical aerodynamic shape design problems is yet an ongoing
endeavor, and it is especially demanding if the number of shape variable $x$
and the dimension of the uncertain parameter $\xi$ is large. In this work,
instead, we solve the present optimization formulations using stochastic
gradient approaches and make use of efficient gradient computations via the
adjoint method and algorithmic differentiation.

###  5.1 Stochastic gradient methods for EMO and DRO

Perhaps the most well-known method to solve the unconstrained empirical mean
optimization problem is the stochastic gradient descent (SGD). Stochastic
gradient methods solve EMOs with a number of gradient computations independent
of the total number of data points $\xi$ and the dimension of the design
variable $x$. If we only solve unconstrained EMOs, we can indeed rely on the
rich optimization algorithms for stochastic optimization with SGD and its
advanced variants.

Solving DRO with stochastic gradient methods is more challenging because the
inner maximization over candidate distributions makes cheap sampling-based
gradient estimates biased (i.e., it exists approximation error). The work [23]
shows that the bias upper bound is $O(\sqrt{1/n})$ and $O(1/n)$ for $\chi^{2}$
constraint and $\chi^{2}$ penalty formulation, respectively, where $n$ is the
batch size. Despite the somewhat pessimistic error bound, the same work
reported promising computational experiments that the error floor becomes
negligible even for batch sizes as small as 10. The work [12] considers the
same unconstrained DRO setting and suggests dynamically increasing the batch
sizes to balance the accuracy and computational efficiency. For practical
aerodynamic shape design, a large batch size is often prohibitively expensive.
Under a constrained computational budget, we choose a fixed batch size $n$ for
our algorithms.

###  5.2 Dealing with additional inequality constraints

One challenge left is to deal with the additional inequality constraints in
($EMO$), ($DRO_{con}$) and ($DRO_{pen}$). Recently, some of the authors have
proposed a gradient method for inequality constrained optimization (GDAM) [4].
In this work, we extend it for the setting of stochastic optimization.

####  5.2.1 GDAM in the deterministic setting

The method can be motivated with the interior-point method for constrained
optimization [29]. Consider the deterministic constrained optimization problem
($DET$), the logarithmic barrier function writes

$\Phi(x)=-\sum_{i=1}^{m}\log(-g_{i}(x)).$ (47)

The classical barrier interior-point method solves a sequence of
$\eta$-subproblems,

$\min_{x}~{}~{}f(x)+\eta\Phi(x),$ (48)

where $\eta>$ is the barrier parameter. As $\eta\rightarrow 0$, the solution
of (48) converges to the original problem ($DET$). Following the double-loop
structure of second-order IPM, one can design a gradient method: 1)
iteratively decrease the barrier parameter $\eta$ in the outer loop, and 2)
solve the subproblem (48) with the negative gradient direction

$d_{\Phi,\eta}=-\nabla f(x)-\eta\nabla\Phi(x).$ (49)

In contrast to the double-loop method, GDAM suggests a descent direction of
the objective $f(x)$ as

$\mathbf{s}_{\zeta}(x)=-\frac{\nabla f(x)}{\|\nabla
f(x)\|}-\zeta\frac{\nabla\Phi(x)}{\|\nabla\Phi(x)\|},~{}~{}\zeta\in[0,1),$
(GDAM)

where $\|\cdot\|$ denotes L2-norm. Scaling GDAM’s search direction $s_{\zeta}$
with $\|\nabla f(x_{k})\|$, we obtain

$\|\nabla f(x)\|s_{\zeta}=-\nabla f(x)-\zeta\frac{\|\nabla
f(x)\|}{\|\nabla\Phi(x)\|}\nabla\Phi(x).$ (50)

Comparing (50) with (49), it is then clear that the (scaled) GDAM suggests a
dynamic computation of the barrier parameter at each step $k$,

$\eta(x_{k})=\zeta\frac{\|\nabla f(x_{k})\|}{\|\nabla\Phi(x_{k})\|}.$ (51)

It was shown in [4] that GDAM is a computationally practical method for shape
optimizations, which motivated us to its extension for solving the present
data-driven shape design problems.

####  5.2.2 GDAM in the stochastic setting

For stochastic optimizations, solving $\eta$-barrier subproblems with
stochastic gradient methods (i.e., with a stochastic gradient estimation for
(49)) has been used in machine learning applications when additional
constraints are considered, see, e.g., [18][24]. The duality gap associated
with primal $x^{\star}$ of the $\eta$-barrier problem and implicit dual
feasible $\lambda^{\star}$ of the original inequality constrained problem
($DET$) is upper bounded by $m\eta$. However, a small $\eta$ leads to high
gradient oscillations near the boundary of the feasible set (in both the
deterministic and stochastic settings).

In this work, we extend GDAM to the stochastic setting. More precisely, the
gradient evaluation of the objective function is considered to be noisy. In
the deterministic setting, GDAM results in optimization trajectories that
travel alongside the central path within a $\Theta_{\zeta}$ neighborhood
defined by the hyperparameter $\zeta$ (a larger $\zeta$ indicates a smaller
cone neighborhood, see [4, Lemma 5.6]). In the stochastic setting, however,
the central path is no longer unique and varies with the gradient noise. To
ensure progress in reducing the objective function value, we relax the
parameter $\zeta=\tau\zeta,0<\tau<1$, when the optimization progress stalls.
Intuitively, we allow a larger central path neighborhood to bound a number of
central paths resulting from the imprecise gradient information. We show a
pseudocode of stochastic GDAM in Algorithm 1, which uses a fix-length step
size rule, i.e., $\|x_{k+1}-x_{k}\|=\beta.$ Furthermore, the algorithm can be
equipped with advanced stochastic optimization methods, such as ADAM [19], to
accelerate the optimization.

For aerodynamic shape design problems, our observation is that a design should
approach the boundary of the feasible set close enough so as to achieve good
performance in the objective function. In the next section, we report
computational experiments and show that stochastic GDAM achieves fairly good
design solutions at the constraint boundaries.

Algorithm 1 Stochastic GDAM

1:Initialization: $x_{0}$, $\zeta$, $\tau$, $\beta$, $n$, $k=0$,

2:while Stopping criteria is not met do

3: Sample $\xi_{1},...,\xi_{n}$ and compute $\nabla
f(x_{k};\xi_{1}),...,\nabla f(x_{k};\xi_{n})$.

4: Compute $g_{i}(x_{k}),i=1,...,m$.

5: Compute
$\tilde{s}_{\zeta}(x_{k})=\frac{1}{n}\sum_{i=1}^{n}\mathbf{s}_{\zeta}(x_{k};\xi_{i})$.

6: Update: $\alpha_{k}\leftarrow\frac{\beta}{\|\tilde{s}_{\zeta}(x_{k})\|}$,
$x_{k+1}\leftarrow x_{k}+\alpha_{k}\tilde{s}_{\zeta}(x_{k}),$ $k\leftarrow
k+1$

7: Update: $\zeta\leftarrow\tau\zeta$, if
$\mathbb{E}_{\hat{P}_{n}}\left[f(x_{k+1};\xi)\right]>\mathbb{E}_{\hat{P}_{n}}\left[f(x_{k};\xi)\right]$

8: (Line search to update $\beta$)

9:end while

##  6 Preliminary computational examples

We apply the above presented concepts frameworks and algorithms to aerodynamic
shape design optimization of the RAE2822 airfoil in transonic turbulent flow
with an uncertain free-stream Mach number that assumes to have a “true”
distribution as a normal distribution $N(\mu,\sigma)$, where $\mu=0.729$ and
$\sigma=0.01$. We draw samples from this distribution to construct an
empirical distributions $\hat{P}_{N}$ to facilitate the experiment. In the
following, we consider the deterministic aerodynamic shape optimization
problem, a probabilistic approach based on the empirical mean and a
distributionally robust design approach.

###  6.1 Deterministic aerodynamic shape optimization of RAE2822 in transonic
turbulent flow

First, we consider the deterministic shape optimization problem based on the
mean Mach number, i.e., $\bar{\xi}=0.729$ (which we refer to as the nominal
Mach number in the following). Furthermore, the lift coefficient is fixed to a
value of $0.724$. This is directly done in the CFD solver by adjusting the
angle of attack in an iterative fashion. The optimization goal is to minimize
the drag coefficient $c_{d}$ while constraining the the pitching moment
$c_{m}$ and the maximum thickness $t$ of the airfoil. The corresponding
optimization problem reads

$\begin{split}&\min_{x}~{}~{}c_{d}(x)\\\
&~{}~{}\text{s.t.}~{}~{}~{}c_{m}(x)\leq 0.093,t(x)\geq 0.12.\end{split}$ (52)

The underlying flow equations are modeled using the Reynolds-averaged Navier-
Stokes equations with a Spalart-Allmaras turbulence model and are solved with
the help of the open-source software SU2 [9]. The solution method is a finite
volume method based on a vertex-based median-dual grid using a Jameson-
Schmidt-Turkel scheme for the spatial discretization. The design variables $x$
are composed of $38$ Hicks-Henne parameters describing the shape. More
specifically, SU2 enables the provision of sensitivites of the objective
function with respect to the design variables with the help of the discrete
adjoint method based on a fixed-point formulation [36]. The software suite
also directly provides a Python framework for design optimization that also
enables the integration of optimizers like GDAM.

In a first experiment, the integrated Python optimizer SLSQP [20] and GDAM are
both applied to solve problem (52). We compare the results of GDAM to the
results of SLSQP since SLSQP is commonly applied in the context of constrained
aerodynamic shape optimization problems based on SU2.

GDAM converges to a design depicted in Figure 5 with a reduction of the drag
coefficient of the initial design $c_{d}=0.013462$ to a drag coefficient of
$c_{d}=0.010655$. Note that SLSQP converges to a similar design (see Figure 5)
with a comparable drag coefficient of $c_{d}=0.010620$.

Figure 5: Resulting shapes of deterministic shape optimization problem with
GDAM (bold line) and SLSQP (dashed line).

###  6.2 Probabilistic aerodynamic shape optimization of RAE2822 in transonic
turbulent flow

The empirical mean optimization problem with an uncertain Mach number as
introduced in ($EMO$) is given by

$\begin{split}&\min_{x}~{}~{}\mathbb{E}_{\xi\sim\mathcal{N}(0.729,0.01)}(c_{d}(x{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1};}\xi))\\\
&~{}~{}\text{s.t.}~{}~{}~{}c_{m}(x{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1};}\bar{\xi})\leq
0.093,t(x)\geq 0.12.\end{split}$ (53)

Based on the algorithm described in Section 5 we find a design by using GDAM
in a probabilistic setting. Here, we make use of $16$ samples from
$\hat{P}_{N}$ in each optimization step. In Figure 6, the resulting design is
compared to the design observed for the deterministic optimization. For the
optimal design based on the empirical mean, we observe that the trailing edge
is thicker. The mean value corresponding to this design, approximated using
$200$ samples, is $\mathbb{E}(c_{d})=0.011355$. The drag coefficient
corresponding to the nominal Mach number is given by $c_{d}=0.011148$. Note,
that the approximated mean value of the design found by the deterministic
optimization is $\mathbb{E}(c_{d})=0.011527$.

Figure 6: Resulting shapes of deterministic shape optimization problem with
GDAM (bold line) and empirical mean optimization problem with stochastic GDAM
(dashed line).

Details on the performance of the design for different Mach numbers in a given
range around the nominal Mach number can be seen in Figure 7. While the drag
coefficient of the optimal design (bold line) is sensitive around the Mach
number of 0.729, it is not very sensitive for the design found by the
probabilistic optimization. The strong increase of the drag coefficient is
also shifted to a higher Mach number. This already shows the advantages of
considering an optimization based on the empirical mean.

Figure 7: Behaviour of the drag coefficient in a range around the nominal
value of the Mach number for the deterministic (bold line) and the
probabilistic design (dashed line).

###  6.3 Distributionally robust aerodynamic design of RAE2822 in transonic
turbulent flow

We finally apply the two different strategies for distributionally robust
optimization using $\chi^{2}$-divergence as introduced in Section 3 for
aerodynamic design. We start with the constrained formulation according to
($DRO_{con}$) resulting in the optimization problem

$\begin{split}&\min_{x}\max_{Q\in\mathdutchcal{D}(\mathcal{N}(0.729,0.01),\rho)}\mathbb{E}_{\xi\sim
Q}\left(c_{d}(x;\xi)\right),\\\
&~{}\textnormal{s.t.}~{}~{}~{}c_{m}(x;\bar{\xi})\leq 0.093,t(x)\geq
0.12.\end{split}$ (54)

Given a value of
$\rho={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.5}$,
stochastic GDAM is applied to solve the optimization problem using $16$
samples in each step. The resulting design, with a drag coefficient of
$c_{d}=0.011272$ for the nominal value, is shown in Figure 8 and only small
differences to the design based found by optimizing the empirical mean.

Figure 8: Resulting shapes of deterministic shape optimization (bold line),
empirical mean optimization (dashed line) and constrained-based
distributionally robust optimization (blue line).

Even though the designs are very similar, Figure 9 shows why the result of the
distributionally robust optimization can be considered more robust regarding
changes in the Mach number. The behaviour for low Mach numbers is comparable
to the result of the probabilistic optimization, while for higher Mach number
the increase in the drag coefficient is further reduced.

Figure 9: Behaviour of the drag coefficient in a range around the nominal
value of the Mach number for the deterministic (bold line), the probabilistic
design (dashed line), and the distributionally robust design (blue line).

We also apply stochastic GDAM for solving the distributionally robust
optimization problem based on penalization $DRO_{pen}$, i.e.,

$\begin{split}&\min_{x}\max_{Q}\\{\mathbb{E}_{\xi\sim
Q}\left(c_{d}(x;\xi)\right)-\delta d(Q,\mathcal{N}(0.729,0.01))\\},\\\
&~{}\textnormal{s.t.}~{}~{}c_{m}(x;\bar{\xi})\leq 0.093,t(x)\geq
0.12\end{split}$ (55)

with a penalty factor of $\delta=0.001$. The optimization uses a number of
$16$ samples per iteration and converges in about $200$ iterations to a
similar design as for the constrained formulation, corresponding to a drag
coefficient of $c_{d}=0.011203$ in the nominal value. Figure 10 shows that the
behaviour for different Mach numbers is also comparable.

Figure 10: Behaviour of the drag coefficient in a range around the nominal
value of the Mach number for the distributionally robust design using the
constraint formulation (bold line) and the penalized formulation (dashed
line).

##  7 Discussions

We formulate, study, and solve data-driven aerodynamic shape design problems
with distributionally robust optimization approaches. Although DRO is a
relatively new framework for engineering design, we point out that there is a
strong connection between DRO and the robust design method of Taguchi, a
proven and successful method in many engineering industries. While the latter
is formulated based on the statistical moments of the cost (performance
measures), the former is based on observed data with uncertainties
(uncontrollable and unknown noise factors). Our preliminary computational
experiments on aerodynamic shape design in transonic turbulent flow show
promising design results.

## References

  * [1] A. Ben-Tal, D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen, _Robust solutions of optimization problems affected by uncertain probabilities_ , Management Science 59 (2013), pp. 341–357.
  * [2] H.G. Beyer and B. Sendhoff, _Robust optimization–a comprehensive survey_ , Computer methods in applied mechanics and engineering 196 (2007), pp. 3190–3218.
  * [3] S.L. Brunton, J. Nathan Kutz, K. Manohar, A.Y. Aravkin, K. Morgansen, J. Klemisch, N. Goebel, J. Buttrick, J. Poskin, A.W. Blom-Schieber, _et al._ , _Data-driven aerospace engineering: reframing the industry with machine learning_ , AIAA Journal 59 (2021), pp. 2820–2847.
  * [4] L. Chen, K.U. Bletzinger, N.R. Gauger, and Y. Ye, _A gradient descent akin method for constrained optimization: algorithms and applications_ , arXiv preprint arXiv:2302.11898 (2023).
  * [5] W. Chen, J.K. Allen, K.L. Tsui, and F. Mistree, _A procedure for robust design: minimizing variations caused by noise factors and control factors_ (1996).
  * [6] C. Dapogny, F. Iutzeler, A. Meda, and B. Thibert, _Entropy-regularized wasserstein distributionally robust shape and topology optimization_ , Structural and Multidisciplinary Optimization 66 (2023), p. 42.
  * [7] E. Delage and Y. Ye, _Distributionally robust optimization under moment uncertainty with application to data-driven problems_ , Operations research 58 (2010), pp. 595–612.
  * [8] M. Drela, _Pros & cons of airfoil optimization_, in _Frontiers of Computational Fluid Dynamics 1998_ , World Scientific, 1998, pp. 363–381.
  * [9] T.D. Economon, F. Palacios, S.R. Copeland, T.W. Lukaczyk, and J.J. Alonso, _SU2: An Open-Source Suite for Multiphysics Simulation and Design_ , AIAA Journal 54 (2016), pp. 828–846.
  * [10] P.M. Esfahani and D. Kuhn, _Data-driven distributionally robust optimization using the Wasserstein metric: Performance guarantees and tractable reformulations_ , arXiv preprint arXiv:1505.05116 (2015).
  * [11] C. Geiersbach and W. Wollner, _A Stochastic Gradient Method with Mesh Refinement for PDE Constrained Optimization under Uncertainty_ , SIAM Journal on Scientific Computing 42 (2020), pp. A2750–A2772.
  * [12] S. Ghosh, M. Squillante, and E. Wollega, _Efficient stochastic gradient descent for learning with distributionally robust optimization_ , arXiv preprint arXiv:1805.08728 (2018).
  * [13] J.y. Gotoh, M.J. Kim, and A.E. Lim, _Robust empirical optimization is almost the same as mean–variance optimization_ , Operations research letters 46 (2018), pp. 448–452.
  * [14] C.R. Gumbert, P.A. Newman, and G.J.W. Hou, _High-Fidelity Computational Optimization for 3-D Flexible Wings: Part II – Effect of Random Geometric Uncertainty on Design_ , Optimization and Engineering 6 (2005), pp. 139–156.
  * [15] L. Huyse, _Solving Problems of Optimization Under Uncertainty as Statistical Decision Problems_ , in _42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit_ , Seattle, WA, USA. American Institute of Aeronautics and Astronautics, AIAA 2001-1519, 2001.
  * [16] L. Jofre and A. Doostan, _Rapid aerodynamic shape optimization under uncertainty using a stochastic gradient approach_ , Structural and Multidisciplinary Optimization 65 (2022), p. 196.
  * [17] M.G. Kapteyn, K.E. Willcox, and A.B. Philpott, _Distributionally robust optimization for engineering design under uncertainty_ , International Journal for Numerical Methods in Engineering 120 (2019), pp. 835–859.
  * [18] H. Kervadec, J. Dolz, J. Yuan, C. Desrosiers, E. Granger, and I.B. Ayed, _Constrained deep networks: Lagrangian optimization via log-barrier extensions_ , in _2022 30th European Signal Processing Conference (EUSIPCO)_. IEEE, 2022, pp. 962–966.
  * [19] D.P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , arXiv preprint arXiv:1412.6980 (2014).
  * [20] D. Kraft, _A Software Package for Sequential Quadratic Programming_ , Tech. Rep., DFVLR, Oberpfaffenhofen, Germany, 1988.
  * [21] H. Lam, _Robust sensitivity analysis for stochastic systems_ , Mathematics of Operations Research 41 (2016), pp. 1248–1275.
  * [22] D.S. Lee, L.F. Gonzalez, J. Périaux, and K. Srinivas, _Robust Design Optimisation Using Multi-Objective Evolutionary Algorithms_ , Computers & Fluids 37 (2008), pp. 565–583.
  * [23] D. Levy, Y. Carmon, J.C. Duchi, and A. Sidford, _Large-scale methods for distributionally robust optimization_ , Advances in Neural Information Processing Systems 33 (2020), pp. 8847–8860.
  * [24] Y. Liu, J. Ding, and X. Liu, _IPO: Interior-point policy optimization under constraints_ , in _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 34. 2020, pp. 4940–4947.
  * [25] J.R. Martins, _Aerodynamic design optimization: Challenges and perspectives_ , Computers & Fluids 239 (2022), p. 105391.
  * [26] V.N. Nair, B. Abraham, J. MacKay, G. Box, R.N. Kacker, T.J. Lorenzen, J.M. Lucas, R.H. Myers, G.G. Vining, J.A. Nelder, _et al._ , _Taguchi’s parameter design: a panel discussion_ , Technometrics 34 (1992), pp. 127–161.
  * [27] H. Namkoong and J.C. Duchi, _Variance-based regularization with convex objectives_ , Advances in neural information processing systems 30 (2017).
  * [28] M. Nemec, D.W. Zingg, and T.H. Pulliam, _Multipoint and multi-objective aerodynamic shape optimization_ , AIAA journal 42 (2004), pp. 1057–1065.
  * [29] J. Nocedal and S.J. Wright, _Numerical optimization_ , Springer, 1999.
  * [30] Y.S. Ong, P.B. Nair, and K.Y. Lum, _Max-Min Surrogate-Assisted Evolutionary Algorithm for Robust Design_ , IEEE Transactions on Evolutionary Computation 10 (2006), pp. 392–404.
  * [31] L. Pardo, _Statistical inference based on divergence measures_ , CRC press, 2018.
  * [32] Y. Polyanskiy and Y. Wu, _Lecture notes on information theory_ , Lecture Notes for ECE563 (UIUC) and 6 (2014), p. 7.
  * [33] H. Rahimian and S. Mehrotra, _Distributionally robust optimization: A review_ , arXiv preprint arXiv:1908.05659 (2019).
  * [34] J. Reuther, A. Jameson, J. Alonso, M. Rimlinger, and D. Saunders, _Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers. AIAA paper 97-0103_ , in _35th Aerospace Sciences Meeting and Exhibit, Reno, Nevada_. 1997.
  * [35] R.T. Rockafellar, _Convex analysis_ , Vol. 11, Princeton university press, 1997.
  * [36] M. Sagebaum, T. Albring, and N.R. Gauger, _High-Performance Derivative Computations using CoDiPack_ , ACM Transactions on Mathematical Software 45 (2019).
  * [37] H.E. Scarf, K. Arrow, and S. Karlin, _A min-max solution of an inventory problem_ , Rand Corporation Santa Monica, 1957.
  * [38] V. Schulz and C. Schillings, _Problem Formulations and Treatment of Uncertainties in Aerodynamic Design_ , AIAA Journal 47 (2009), pp. 646–654.
  * [39] V. Schulz and C. Schillings, _Optimal aerodynamic design under uncertainty_ , in _Management and Minimisation of Uncertainties and Errors in Numerical Aerodynamics: Results of the German collaborative project MUNA_ , Springer, 2013, pp. 297–338.
  * [40] M.W. Trosset, _Taguchi and robust optimization_ , Tech. Rep., 1996.
  * [41] K.L. Tsui, _An overview of taguchi method and newly developed statistical methods for robust design_ , Iie Transactions 24 (1992), pp. 44–57.
  * [42] R. Unal and E.B. Dean, _Taguchi approach to design optimization for quality and cost: an overview_ , in _1991 Annual conference of the international society of parametric analysts_. 1990.
  * [43] G.G. Vining and R.H. Myers, _Combining taguchi and response surface philosophies: a dual response approach_ , Journal of quality technology 22 (1990), pp. 38–45.
  * [44] Z. Wang, P.W. Glynn, and Y. Ye, _Likelihood robust optimization for data-driven problems_ , Computational Management Science 13 (2016), pp. 241–261.


# Lepton masses and mixing in a two-Higgs-doublet model

Piotr Chaber  Institute of Physics, University of Silesia, 75 Pułku Piechoty
1, 41-500 Chorzów, Poland  Bartosz Dziewit  Institute of Physics, University
of Silesia, 75 Pułku Piechoty 1, 41-500 Chorzów, Poland  Jacek Holeczek
Institute of Physics, University of Silesia, 75 Pułku Piechoty 1, 41-500
Chorzów, Poland  Monika Richter  Institute of Physics, University of Silesia,
75 Pułku Piechoty 1, 41-500 Chorzów, Poland  Marek Zrałek  Institute of
Physics, University of Silesia, 75 Pułku Piechoty 1, 41-500 Chorzów, Poland
Sebastian Zajac  Faculty of Mathematics and Natural Studies, Cardinal Stefan
Wyszynski University, ul. Woycickiego 1/3, 01-938 Warsaw, Poland

###### Abstract

Within the framework of the two-Higgs Doublet Model (2HDM), we attempt to find
some discrete, non-abelian flavour symmetry which could provide an explanation
for the masses and mixing matrix elements of leptons. Unlike the Standard
Model, currently there is no need for the flavour symmetry to be broken. With
the GAP program we investigate all finite subgroups of the U3 group up to the
order of 1025. Up to such an order there is no group for which it is possible
to select free model parameters in order to match the masses of charged
leptons, masses of neutrinos, and the Pontecorvo-Maki-Nakagawa-Sakata mixing
matrix elements in a satisfactory manner.

##  I I. Introduction

Despite the success of the Standard Model (SM) in providing a description of
the current experimental data, there is a widespread belief that sooner or
later an increase in available energy or accuracy of measurements will lead to
the detection of a discrepancy between experimental results and theoretical
predictions. The SM cannot be considered a complete theory because it does not
provide answers to many pressing questions. One of the most important issues
to be resolved concerns the masses of the fundamental constituents of matter,
quarks and leptons. At the moment we still cannot theoretically predict their
masses; we are only able to obtain their values from experimental data.

The discovery of the Higgs particle offers a partial solution to the problem,
yet it does not resolve it completely. Particles acquire their mass by means
of interaction with the Higgs field, and rather than examine the numerical
values of the masses, we are currently more interested in the question of why
particles interact so differently with the Higgs field, or in other words, why
the Yukawa couplings cannot be theoretically predicted.

The solution to the problem of the elementary particles’ masses is important
in itself, because it would reduce the amount of unknown free parameters in
the SM. Another important reason for conducting such investigations is that
they create a great opportunity to understand the origin of the masses of
physical bodies. The main part of the mass of each physical body comes from
the interaction of the ingredients contained in it. But it is not the entire
mass; the remaining part (though small) being the masses of individual
fermions, which are still undetermined.

Several proposals to solve this problem, at least partially, can be found in
the literature (see e.g. Ishimori  _et al._ (2010); Blum  _et al._ (2008);
Frampton and Kephart (1995)). Although the problem concerns all matter
constituents, here we will concentrate on an attempt to explain the masses and
mixing angles of leptons. One of the most common approaches consists in the
imposition of a flavour symmetry on the leptonic part of Yukawa Lagrangian
(for review see e.g. King (2017); Altarelli, Guido and Feruglio, Ferruccio
(2010); King and Luhn (2013)). This approach was particularly popular and
successful before 2011 when it was discovered that the reactor-mixing angle
$\theta_{13}$ is non-vanishing Harrison  _et al._ (2002); Harrison and Scott
(2003). The models with an additional flavour symmetry are very popular, but
they are by no means the only ones – some papers have been published in which
the very existence of such a symmetry is denied de Gouvea and Murayama (2015);
Hall  _et al._ (2000).

Attempts at solving the problem of lepton masses by a horizontal symmetry are
dependent on the manner of the introduction of neutrino masses as well as on
the Higgs sector for spontaneous symmetry breaking. In the simplest case of
the conventional SM with one Higgs doublet in which neutrinos are massless,
only three additional right-handed neutrinos are introduced. In such an
extension of the SM, without introducing the Majorana term, neutrinos are
Dirac particles Giunti and Kim (2007). It is not necessary to introduce the
right-handed neutrinos to obtain their masses. Instead, it is possible to use
the existing left-handed neutrinos to form the Majorana masses Giunti and Kim
(2007). Both cases will be considered in this paper.

Within the framework of the Standard Model with one Higgs boson, a discrete
symmetry for Yukawa couplings provides the relations for the three-dimensional
mass matrices of charged leptons $(M_{l})$ and neutrinos $(M_{\nu})$ Lam
(2008a, b):

$\displaystyle A^{i{\dagger}}_{L}\left(M_{l}M_{l}^{\dagger}\right)A_{L}^{i}$
$\displaystyle=$ $\displaystyle\left(M_{l}M_{l}^{\dagger}\right),$ (1)
$\displaystyle
A^{i{\dagger}}_{L}\left(M_{\nu}M_{\nu}^{\dagger}\right)A_{L}^{i}$
$\displaystyle=$ $\displaystyle\left(M_{\nu}M_{\nu}^{\dagger}\right),$ (2)

where $A^{i}_{L}=A_{L}(g_{i}),\;i=1,2,\ldots,N$ are three dimensional
representation matrices for the left-handed lepton doublets for some N-order
flavour symmetry group $\mathcal{G}$.

In such a case, the I-st Schur’s lemma implies that $M_{l}M_{l}^{{\dagger}}$
and $M_{\nu}M_{\nu}^{{\dagger}}$ are proportional to the identity matrices,
which clearly entails the trivial lepton mixing matrix (known in the
literature as Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix Pontecorvo (1968);
Maki  _et al._ (1962); Dziewit  _et al._ (2013)).

In order to avoid the relations given in (2) and predict the non-trivial
lepton masses and their mixing in this case, the family symmetry has to be
broken. As a rule, the flavour symmetry $\mathcal{G}$ is spontaneously broken
by scalar singlet Higgs fields called flavons (see e.g. (King and Luhn, 2009;
Altarelli, Guido and Feruglio, Ferruccio, 2010; King  _et al._ , 2014)).
However, it can also be broken by introducing a bigger number of normal Higgs
multiplets (e.g. King  _et al._ (2014); Machado  _et al._ (2011)). The latter
way is indeed more economical, since the spontaneous gauge symmetry breaking
in this case gives rise to the particles’ masses, simultaneously leading to
the break of a family symmetry. Additional flavon scalar fields in this
framework are therefore redundant.

Such models were considered many times in the literature, but mostly in the
frame of supersymmetric models where Higgs doublets were singlets of a flavour
group Siyeon (2014); Nebot (2016), or, in a more general approach, where only
one selected flavour group was tested Kubo  _et al._ (2003); Nguyen (2014);
Morisi and Peinado (2009).

In the present work we attempt to explore how the Two Higgs Doublet Model
works in the context discussed above. As distinct from the previous works, we
do not consider a few selected discrete groups, but instead we try to find a
flavour symmetry in all groups up to 1025 order with one restriction, i.e.
each of our groups must have at least one faithful, three dimensional
irreducible representation.

In the next Section, we briefly introduce the flavour symmetry in the 2HDM
model and show how the symmetry transformation between two Higgs doublets
provides an opportunity to avoid the consequences of Schur’s lemma. We also
present all the formulas needed to conduct the computations in the case of
Dirac and Majorana neutrinos. In Sec. III the results of the final scan of the
Yukawa matrices, the lepton masses and the PMNS mixing matrix elements are
presented, and finally in Sec. IV we draw our conclusions.

##  II II. Two-Higgs-Doublet-Model with a flavour symmetry

###  II.1  A. Dirac neutrinos

To begin, the leptonic part of the Yukawa Lagrangian with Dirac neutrinos will
be considered. In contrast to the Standard Model, two Higgs doublets
$\Phi_{i}$ contribute to the lepton masses (the so-called Two-Higgs-Doublet-
Model of type III Hou (1992)) as follows Branco  _et al._ (2012)111Note that,
in comparison to our notation, in the paper Branco  _et al._ (2012):
$\Phi_{i}=\left(\phi_{i}^{+},\phi_{i}^{0}/\sqrt{2}\right)^{T},\;i=1,2.$

$\displaystyle\mathcal{L}_{Y}=$ $\displaystyle-$
$\displaystyle\sum\limits_{i=1,2}\sum\limits_{\alpha,\beta=e,\mu,\tau}\left((h_{i}^{(l)})_{\alpha,\beta}\left[\overline{L}_{\alpha
L}\tilde{\Phi}_{i}l_{\beta R}\right]\right.$ (3) $\displaystyle+$
$\displaystyle\left.(h^{(\nu)}_{i})_{\alpha\beta}\Big{[}\overline{L}_{\alpha
L}\Phi_{i}\nu_{\beta R}\Big{]}\right)+h.c.$

where:

$L_{\alpha L}=\left(\begin{array}[]{c}\nu_{\alpha L}\\\ \l_{\alpha
L}\end{array}\right),\quad\Phi_{i}=\left(\begin{array}[]{c}\phi_{i}^{0}\\\
\phi_{i}^{-}\end{array}\right),\;i=1,2$

are gauge doublets for the left-handed lepton and Higgs fields and the fields
$l_{\beta R},\;\nu_{\beta R}$ stand for the right-handed lepton and neutrino
fields, respectively. The couplings $h^{(l)}_{i}$ and $h^{(\nu)}_{i}$ create
the 3-dimensional Yukawa matrices.

The spontaneous gauge symmetry breaking gives non-zero vacuum expectation
values (VEVs) $v_{i}$ for the Higgs doublets:

$\left<\Phi_{i}\right>=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}v_{i}\\\
0\end{array}\right),$ (4)

and the mass matrices read as follows Branco  _et al._ (2012):

$\displaystyle M^{l}$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}\left(v_{1}^{*}h^{(l)}_{1}+v_{2}^{*}h^{(l)}_{2}\right),$
(5) $\displaystyle M^{\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left(v_{1}h^{(\nu)}_{1}+v_{2}h^{(\nu)}_{2}\right).$
(6)

In general the vacuum expectation values can be complex,
$v_{i}=|v_{i}|e^{i\varphi_{i}}$, but they are restricted by the Fermi coupling
constant:

$\sqrt{|v_{1}|^{2}+|v_{2}|^{2}}=\left(\sqrt{2}G_{F}\right)^{-1/2}\simeq 246\
GeV.$ (7)

Family symmetry of our theory implies, that after the transformation of fields
occurring in the 2HDM Lagrangian by the 3 dimensional
($A_{L},\;A^{R}_{l},\;A^{R}_{\nu}$) and 2 dimensional ($A_{\Phi}$)
representations of a flavour group $\mathcal{G}$:

$\displaystyle L_{\alpha L}\rightarrow L^{\prime}_{\alpha L}$ $\displaystyle=$
$\displaystyle\left(A_{L}\right)_{\alpha,\chi}L_{\chi L},\ l_{\beta
R}\rightarrow l^{\prime}_{\beta
R}=\left(A^{R}_{l}\right)_{\beta,\delta}l_{\delta R}$ $\displaystyle\nu_{\beta
R}\rightarrow\nu^{\prime}_{\beta R}$ $\displaystyle=$
$\displaystyle\left(A^{R}_{\nu}\right)_{\beta,\delta}\nu_{\delta R},\
\Phi_{i}\rightarrow\Phi^{\prime}_{i}=\left(A_{\Phi}\right)_{ik}\Phi_{k},$ (8)

the full 2HDM Lagrangian does not change:

$\mathcal{L}\left(L_{\alpha L},l_{\beta R},\nu_{\gamma
R},\Phi_{i}\right)=\mathcal{L}\left(L^{\prime}_{\alpha L},l^{\prime}_{\beta
R},\nu^{\prime}_{\gamma R},\Phi^{\prime}_{i}\right).$ (9)

Given that all the transformation matrices in Eq.(II.1) are unitary, the only
parts of the total 2HDM Lagrangian for which the aforementioned relations are
not automatically fulfilled are the Yukawa Lagrangian and the Higgs potential.
Imposition of symmetry on these terms of the model severely restricts their
forms.

The invariance of the Yukawa Lagrangian is expressed as follows:

$\displaystyle\mathcal{L}^{{}^{\prime}}_{Y}\equiv$
$\displaystyle-\sum_{i=1,2}\sum_{\alpha,\beta=e,\mu,\tau}\left(\left(h^{(l)}_{i}\right)_{\alpha,\beta}\left[\overline{L^{\prime}}_{\alpha
L}\tilde{\Phi^{\prime}}_{i}l^{\prime}_{\beta R}\right]\right.$ (10)
$\displaystyle\left.+\left(h^{(\nu)}_{i}\right)_{\alpha\beta}\Big{[}\overline{L^{\prime}}_{\alpha
L}\Phi^{\prime}_{i}\nu^{\prime}_{\beta R}\Big{]}\right)+h.c.=\mathcal{L}_{Y}.$

With regard to the Higgs potential, there appear to be two possibilities. The
first one assumes that:

$V\left(\Phi^{\prime}_{1},\Phi^{\prime}_{2}\right)=V\left(\Phi_{1},\Phi_{2}\right),$
(11)

which implies that before and after the transformation for Higgs fields
(Eq.(II.1)) the coefficients in the potential remain exactly the same and the
vacuum expectation values are equal, $v_{i}^{\prime}=v_{i}$.

There is also a second possibility, useful for phenomenological reasons, we
allow for the modification of VEVs, which transform in the same way as the
Higgs fields:

$v^{\prime}_{i}=\left(A_{\Phi}\right)_{ik}v_{k}.$ (12)

In this case, the form of the Higgs potential does not change, the terms in
the potential do not vary, while only the potential coefficients undergo
change. This kind of invariance is known in the literature as the form-
invariance (see Branco  _et al._ (2012); Davidson and Haber (2005)). After the
unitary transformation (Eq.(12)), the condition given in Eq.(7) is unchanged,
and hence also:

$\sqrt{|v_{1}^{{}^{\prime}}|^{2}+|v_{2}^{{}^{\prime}}|^{2}}\simeq 246\ GeV.$
(13)

The vacuum expectation values, whose sum of squares is constant, will need
adjustments to meet the experimental requirements. Thus, from the point of
view of flavour symmetry, the type of Higgs potential is irrelevant, so in our
approach, the issue of what symmetry for Higgs Lagrangian is chosen becomes
insignificant.

With reference to the Eq.(9), in order to find symmetric Yukawa matrices
$h^{(l)}_{i},\;h^{(\nu)}_{i}$, $i=1,2$, one can readily express the symmetry
conditions as the eigenequation for a direct product of unitary group
representations to the eigenvalue 1 (see e.g.Ludl (6 27)):

$\displaystyle\left((A_{\Phi})^{{\dagger}}\otimes(A_{L})^{{\dagger}}\otimes(A^{R}_{l})^{T}\right)_{k,\alpha,\delta;i,\beta,\gamma}(h^{l}_{i})_{\beta,\gamma}$
$\displaystyle=$ $\displaystyle(h^{l}_{k})_{\alpha,\delta},$
$\displaystyle\left((A_{\Phi})^{T}\otimes(A_{L})^{{\dagger}}\otimes(A^{R}_{\nu})^{T}\right)_{k,\alpha,\delta;i,\beta,\gamma}(h^{\nu}_{i})_{\beta,\gamma}$
$\displaystyle=$ $\displaystyle(h^{\nu}_{k})_{\alpha,\delta},$ (14)

for the charged leptons and neutrinos, respectively.

Both relations (Eq.(II.1)) need to be satisfied for any group’s element
$g\in\mathcal{G}$. It is however sufficient that they are fulfilled only for
the group generators Ludl (6 27), which considerably reduces the time of the
computation.

In such a model, the invariance equations for the mass matrices are not
trivial. For the symmetric Higgs potential (Eq.(11)) :

$A_{L}M^{l(\nu)}\left(A^{R}_{l(\nu)}\right)^{{\dagger}}=\frac{1}{\sqrt{2}}\sum_{i,k=1}^{2}h^{l(\nu)}_{i}\left(A_{\Phi}\right)_{i,k}v_{k}\neq
M^{l(\nu)}.$ (15)

then Eq.(1-2) are not satisfied and we avoid the consequences of the Schur’s
Lemma. The same can be shown for the form-invariant Higgs potential, where
Eq.(12) is satisfied. In this context we can obtain the non-trivial mass
matrices without the introduction of additional flavon fields.

###  II.2  B. Majorana neutrinos

For Majorana neutrinos the Yukawa term has to be changed. In 2HDM, the
simplest Yukawa Lagrangian can be taken as the non-renormalizable Weinberg
term in the form:

$\displaystyle\mathcal{L}_{Y}^{\nu}\equiv$
$\displaystyle-\frac{g}{M}\sum\limits_{i,k=1}^{2}\sum\limits_{\alpha,\beta=e,\mu,\tau}h^{(i,k)}_{\alpha,\beta}\left(\overline{L}_{\alpha
L}{\Phi}_{i}\right)\left({\Phi}_{k}{L}_{\beta R}^{c}\right)$ (16)
$\displaystyle+h.c,$

where ${L}_{\beta R}^{c}=C\overline{L}_{\beta L}^{T}$ is the charge conjugated
lepton doublet fields. After the spontaneous symmetry breaking, the neutrino
mass matrix is obtained:

$M^{\nu}_{\alpha,\beta}=\frac{g}{M}\sum_{i,k=1}^{2}v_{i}v_{k}h^{(i,k)}_{\alpha,\beta}.$
(17)

As in the preceding case, in compliance with the requirement of flavour
symmetry for the Yukawa Lagrangian (Eq.(16)), the neutrino Yukawa matrices
must satisfy the eigenvalue equation:

$\displaystyle\left((A_{\Phi})^{T}\otimes\left(A_{\Phi}\right)^{T}\otimes\left(A_{L}\right)^{{\dagger}}\otimes\left(A_{L}\right)^{{\dagger}}\right)_{k,m,\chi,\eta;\
i,j,\alpha,\beta}\left(h^{(i,j)}_{\alpha,\beta}\right)=\left(h^{(k,m)}_{\chi,\eta}\right).$
(18)

Such flavour symmetric Yukawa couplings restrict the neutrino mass matrix from
Eq.(17).

##  III III. Results of the investigation for finite groups up to 1025 order

###  III.1  A. The candidates for the flavour group $\mathcal{G}_{F}$

The flavour group $\mathcal{G}_{F}$, which is imposed on the 2HDM Lagrangian,
cannot be arbitrary. Due to the fields’ transformations defined in Eq.(II.1),
the group must possess at least one 2-dimensional (for $A_{\Phi}$) and at
least one 3-dimensional irreducible representation (for $A_{L}$, $A^{R}_{\nu}$
and $A^{R}_{l}$). The sole application of irreducible representations is
justified given the fact that were any of the representations
($A_{\Phi},A_{L}$ or $A^{R}_{\nu}(A^{R}_{l})$) to be reducible, invariant
Yukawa couplings would split up into independent sets for irreducible
representations (see e.g.Ludl (6 27)).

In the selection of flavour symmetry groups, we have limited ourselves to
finite-dimensional groups of the order of at most 1025, which are furthermore
subgroups of the $U(3)$ group (at least one of the 3-dimensional irreducible
or reducible representations must be faithful) Ludl (2010); Parattu and
Wingerter (2011). This additional condition is not necessary Grimus and Ludl
(2012) but it significantly reduces the number of groups to be processed.
Using the GAP (version 4.7.6) GAP  system for computational discrete algebra
with the included Small Groups Library Besche  _et al._ (2014) and REPSN
Dabbaghian (2011) packages, we have found in total 10862 groups with at least
one 2-dimensional and at least one 3-dimensional irreducible representation,
but only 413 of these groups are subgroups of the $U(3)$ group. They split
into two disjoint sets. Each group has either at least one faithful
3-dimensional irreducible representation (there are 173 such groups), or at
least one faithful 1+2 reducible representation (there are 240 such groups).
Some groups are also subgroups of the $U(2)$ group. They have at least one
faithful 2-dimensional irreducible representation (none of them have any
faithful 3-dimensional irreducible representation). None of the groups have
any faithful 1+1+1 reducible, faithful 1+1 reducible or faithful 1-dimensional
irreducible representation. All the groups which delivered any solutions
belong to polycyclic groups that use the polycyclic presentation for element
arithmetic (so called $PC$-groups).

Table 1:  Groups of the order of at most 100 subject to consideration (all the groups are listed in the ancillary file). Here: “$[\ o\ ,\ i\ ]$” the $i$-th group of the order $o$ in the Small Groups Library catalogue, “StructureDescription” a short string which provides some insight into the structure of the group under consideration, “2-D” the number of 2-dimensional irreducible representations, “3-D” the number of 3-dimensional irreducible representations, “$U(2)$” an indicator why the group is classified as a subgroup of the $U(2)$ group (at least one 2-dimensional irreducible is faithful), “$U(3)$” an indicator why the group is classified as a subgroup of the $U(3)$ group (either at least one 3-dimensional irreducible or one 1+2 reducible representation is faithful), “L” the number of different combinations of representations for charged leptons, “DN” the number of different combinations of representations for Dirac neutrinos, “MN” the number of different combinations of representations for Majorana neutrinos, “L+DN” the number of pairs of different combinations of representations for charged leptons and Dirac neutrinos, “L+MN” the number of pairs of different combinations of representations for charged leptons and Majorana neutrinos. Note that the “L” and the “DN” are always equal and that the “L+DN” is twice that number. All zero values are suppressed. $[\ o\ ,\ i\ ]$ | StructureDescription | 2-D | 3-D | $U(2)$ | $U(3)$ | L | DN | MN | L+DN | L+MN
---|---|---|---|---|---|---|---|---|---|---
$[\ 24\ ,\ 3\ ]$ | SL(2,3) | 3 | 1 | 2 | 1+2 |  |  | 3 |  |
$[\ 24\ ,\ 12\ ]$ | S4 | 1 | 2 |  | 3 | 4 | 4 | 2 | 8 | 4
$[\ 48\ ,\ 28\ ]$ | C2.S4=SL(2,3).C2 | 3 | 2 | 2 | 1+2 | 4 | 4 | 6 | 8 | 4
$[\ 48\ ,\ 29\ ]$ | GL(2,3) | 3 | 2 | 2 | 1+2 | 4 | 4 | 6 | 8 | 4
$[\ 48\ ,\ 30\ ]$ | A4:C4 | 2 | 4 |  | 3 | 16 | 16 | 8 | 32 | 16
$[\ 48\ ,\ 32\ ]$ | C2xSL(2,3) | 6 | 2 |  | 1+2 |  |  | 12 |  |
$[\ 48\ ,\ 33\ ]$ | SL(2,3):C2 | 6 | 2 | 2 | 1+2 |  |  |  |  |
$[\ 48\ ,\ 48\ ]$ | C2xS4 | 2 | 4 |  | 3 | 16 | 16 | 8 | 32 | 16
$[\ 54\ ,\ 8\ ]$ | ((C3xC3):C3):C2 | 4 | 4 |  | 3 | 32 | 32 |  | 64 |
$[\ 72\ ,\ 3\ ]$ | Q8:C9 | 9 | 3 | 2 | 1+2 |  |  | 9 |  |
$[\ 72\ ,\ 25\ ]$ | C3xSL(2,3) | 9 | 3 | 2 | 1+2 |  |  | 9 |  |
$[\ 72\ ,\ 42\ ]$ | C3xS4 | 3 | 6 |  | 3 | 36 | 36 | 6 | 72 | 12
$[\ 96\ ,\ 64\ ]$ | ((C4xC4):C3):C2 | 1 | 6 |  | 3 | 12 | 12 | 2 | 24 | 4
$[\ 96\ ,\ 65\ ]$ | A4:C8 | 4 | 8 |  | 3 | 64 | 64 | 16 | 128 | 32
$[\ 96\ ,\ 66\ ]$ | SL(2,3):C4 | 6 | 4 |  | 1+2 | 16 | 16 | 24 | 32 | 16
$[\ 96\ ,\ 67\ ]$ | SL(2,3):C4 | 6 | 4 | 2 | 1+2 | 16 | 16 | 8 | 32 | 16
$[\ 96\ ,\ 69\ ]$ | C4xSL(2,3) | 12 | 4 |  | 1+2 |  |  | 24 |  |
$[\ 96\ ,\ 74\ ]$ | ((C8xC2):C2):C3 | 12 | 4 | 2 | 1+2 |  |  |  |  |
$[\ 96\ ,\ 186\ ]$ | C4xS4 | 4 | 8 |  | 3 | 64 | 64 | 16 | 128 | 32
$[\ 96\ ,\ 188\ ]$ | C2x(C2.S4=SL(2,3).C2) | 6 | 4 |  | 1+2 | 16 | 16 | 24 | 32 | 16
$[\ 96\ ,\ 189\ ]$ | C2xGL(2,3) | 6 | 4 |  | 1+2 | 16 | 16 | 24 | 32 | 16
$[\ 96\ ,\ 192\ ]$ | (C2.S4=SL(2,3).C2):C2 | 6 | 4 | 2 | 1+2 | 16 | 16 | 8 | 32 | 16
$[\ 96\ ,\ 200\ ]$ | C2x(SL(2,3):C2) | 12 | 4 |  | 1+2 |  |  |  |  |

All the detected groups of the order of at most 1025 are listed in the
ancillary file, while Table 1 shows groups of the order of at most 100 only
which underwent investigation.

###  III.2 B. Yukawa couplings matrices in the model with Dirac neutrinos

In the case of all the groups subject to consideration, there exist 267 groups
that gave in total 748672 different combinations of 2- and 3-dimensional
irreducible representations that give 1-dimensional degeneration subspace for
all generators, which is the solution to the equations in Eq.(II.1). This
common vector gives Yukawa matrices for charged leptons ($h^{(l)}$) and for
neutrinos ($h^{(\nu)}$), which are interrelated. All the possible solutions
for Yukawa matrices for charged leptons and for Dirac neutrinos can be
expressed through 7 base forms ($\omega=e^{2\pi i/3}$):

$\displaystyle h_{1}^{(1)}=\left(\begin{array}[]{ccc}0&0&1\\\ 1&0&0\\\
0&1&0\\\ \end{array}\right),\ \quad
h_{2}^{(1)}=\left(\begin{array}[]{ccc}0&1&0\\\ 0&0&1\\\ 1&0&0\\\
\end{array}\right),$ (25)

$\displaystyle h_{1}^{(2)}=\left(\begin{array}[]{ccc}0&0&1\\\
\omega^{2}&0&0\\\ 0&\omega&0\\\ \end{array}\right),\quad
h_{2}^{(2)}=\left(\begin{array}[]{ccc}0&1&0\\\ 0&0&\omega\\\ \omega^{2}&0&0\\\
\end{array}\right),$ (32) $\displaystyle
h_{1}^{(3)}=\left(\begin{array}[]{ccc}0&0&1\\\ \omega&0&0\\\ 0&\omega^{2}&0\\\
\end{array}\right),\quad h_{2}^{(3)}=\left(\begin{array}[]{ccc}0&1&0\\\
0&0&\omega^{2}\\\ \omega&0&0\\\ \end{array}\right),$ (39)

the next three $h_{1}^{(i)}$ and $h_{2}^{(i)}$ for i = 4,5,6, are obtained
from those given in Eq.(25)-Eq.(39) by interchange, using the rule:

$h_{1}^{(3+i)}=h_{2}^{(i)},\ \ h_{2}^{(3+i)}=h_{1}^{(i)}$

for i= 1,2,3, and finally222$U(3)$ subgroups which have no faithful
3-dimensional irreducible representation give only these diagonal solutions.:

$\displaystyle h_{1}^{(7)}=\left(\begin{array}[]{ccc}1&0&0\\\
0&\omega^{2}&0\\\ 0&0&\omega\\\ \end{array}\right),\quad
h_{2}^{(7)}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&\omega&0\\\ 0&0&\omega^{2}\\\
\end{array}\right).$ (46)

Symmetric Yukawa matrices for any considered group and for any irreducible
representation within the group can be expressed now by seven (i=1,2,..,7)
basic matrix forms as follows:
a) for Dirac neutrinos:

$\displaystyle\\{h_{1}^{(\nu)},h_{2}^{(\nu)}\\}=\\{h_{1}^{(i)},e^{i\varphi}h_{2}^{(i)}\\},$
(47)

b) for charged leptons:

$\displaystyle\\{h_{1}^{(l)},h_{2}^{(l)}\\}=\\{h_{2}^{(i)},e^{-i(\delta_{l}+\varphi)}h_{1}^{(i)}\\},$
(48)

where $\varphi$ is a phase distinctive for a group and for irreducible
representations and $\delta_{l}=0,\pi.$
In order to find the lepton masses and mixing matrix we constructed the
hermitian matrices as in Eq.(1) $M^{l}M^{l\dagger}$ and
$M^{\nu}M^{\nu\dagger}$. For all the possible Yukawa matrices we have obtained
only three different forms $(x=l,\nu)$:

$\displaystyle
M_{x}M_{x}^{\dagger}=|c_{x}|^{2}\left(\begin{array}[]{ccc}1+\kappa^{2}&\kappa
e^{-i(\eta_{x}+2k\pi/3)}&\kappa e^{i(\eta_{x}-2k\pi/3)}\\\ \kappa
e^{i(\eta_{x}+2k\pi/3)}&1+\kappa^{2}&\kappa e^{-i\eta_{x}}\\\ \kappa
e^{-i(\eta_{x}-2k\pi/3)}&\kappa e^{i\eta_{x}}&1+\kappa^{2}\\\
\end{array}\right),$ (52)

with $k=-1,0,+1$ and $\kappa=|v_{2}|/|v_{1}|$, the same for neutrinos and for
charged leptons. The only difference lies in the phase $\eta_{x}$. For Dirac
neutrinos,

$\eta_{\nu}=\varphi+\varphi_{2}-\varphi_{1},$

and for charged leptons,

$\eta_{l}=\delta_{l}+\varphi+\varphi_{2}-\varphi_{1},$

where $\varphi_{i}(i=1,2)$ are phases of the VEVs $v_{i}$. After
diagonalization of Eq.(52)

$U^{\dagger}\left(M_{x}M_{x}^{\dagger}\right)U=diag\left(m_{x1}^{2},m_{x2}^{2},m_{x3}^{2}\right),$

we obtain:

$\displaystyle m_{x1}^{2}$
$\displaystyle=|c_{x}|^{2}\left(1+\kappa^{2}+2\kappa
cos\left(\eta_{x}\right)\right),$ (53) $\displaystyle m_{x2}^{2}$
$\displaystyle=|c_{x}|^{2}\left(1+\kappa^{2}+2\kappa
sin\left(\eta_{x}-\frac{\pi}{6}\right)\right),$ (54) $\displaystyle
m_{x3}^{2}$ $\displaystyle=|c_{x}|^{2}\left(1+\kappa^{2}-2\kappa
sin\left(\eta_{x}+\frac{\pi}{6}\right)\right),$ (55)

and the diagonalization matrix U:

$U=\frac{1}{\sqrt{3}}\left(\begin{array}[]{ccc}e^{-\frac{2}{3}\pi ik}&\omega
e^{-\frac{2}{3}\pi ik}&\omega^{2}e^{-\frac{2}{3}\pi ik}\\\
1&\omega^{2}&\omega\\\ 1&1&1\\\ \end{array}\right).$ (56)

This matrix does not depend on the phase $\eta_{x}$, so it is identical for
charged leptons and for the neutrino. Therefore, it is not possible to
reconstruct the correct mixing matrix. For groups and their irreducible
representations for which $\delta_{l}=0$, neutrinos and lepton masses are
proportional and the mixing matrix is 3x3 identity matrix. For groups and for
representations where $\delta_{l}=\pi$, masses of charged leptons and
neutrinos are not proportional. If the formulas in Eqs.(53-55) describe the
masses of neutrinos, then for the charged leptons there is:

$\displaystyle m_{l1}^{2}$
$\displaystyle=|c_{l}|^{2}\left(1+\kappa^{2}-2\kappa
cos\left(\eta_{\nu}\right)\right),$ (57) $\displaystyle m_{l2}^{2}$
$\displaystyle=|c_{l}|^{2}\left(1+\kappa^{2}-2\kappa
sin\left(\eta_{\nu}-\frac{\pi}{6}\right)\right),$ (58) $\displaystyle
m_{l3}^{2}$ $\displaystyle=|c_{l}|^{2}\left(1+\kappa^{2}+2\kappa
sin\left(\eta_{\nu}+\frac{\pi}{6}\right)\right).$ (59)

Similarly in this case we cannot reconstruct the PMNS mixing matrix for which
the 3x3 anti-diagonal identity matrix is obtained. Regardless of whether the
masses of charged leptons are described by formula Eqs.(53-55)) (groups with
$\delta_{l}=0$) or by Eq.(57-59) (groups with $\delta_{l}=\pi$), three
parameters $(c_{l},\kappa,\eta_{x})$ cannot be selected in such a way as to
obtain the physical masses of the electron, muon and tau (in any case, one
obtains: $0\leqslant m_{e}/m_{\mu}\leqslant 1$ and $1\leqslant
m_{\tau}/m_{\mu}\leqslant 2$).

The SM extended by one additional doublet of Higgs particles (2HDM) does not
posses a discrete family symmetry (in the groups under examination) that can
explain the masses of charged leptons, masses of neutrinos having the nature
of Dirac particles and the PMNS matrix.

###  III.3 C. Yukawa couplings matrices in the model with Majorana neutrinos

In the previous subsection, we did not receive masses of charged leptons
consistent with the experiment. This was the case when the neutrinos were
Dirac particles. In the case of Majorana neutrinos, we must once again look
for a possible symmetry, since a symmetric solution for the Weinberg component
may also deliver other groups of symmetry for charged leptons. We are
currently looking for symmetries satisfying the Eg.(II.1) for the charged
leptons and Eq.(18) for neutrinos. Out of all the groups under consideration,
there exist 195 groups that gave in total 20888 solutions. All the found
symmetries, as a solution of Eq.(18), give a two dimensional space which is
common to all generators of the groups in question. Each time we find two
36-dimensional vectors p and r and any linear combination of p and r, gives
symmetric Yukawa matrices for Majorana neutrinos. Thus these Yukawa matrices
are given by:

$h^{(i,k)}=xp_{i,k}+yr_{i,k},$ (60)

where x and y are two free complex numbers.
A detailed analysis of the solution for Majorana neutrinos is as follows:

$\displaystyle h^{(1,1)}=xh_{2}^{(7)},\quad$ $\displaystyle h^{(1,2)}=yI_{3},$
(61) $\displaystyle h^{(2,1)}=ye^{i\delta}I_{3},\quad$ $\displaystyle
h^{(2,2)}=xe^{i(\delta+2\varphi)}h_{1}^{(7)},$

and for the charged leptons Yukawa matrices:

$\displaystyle\\{h_{1}^{(l)},h_{2}^{(l)}\\}=\\{h_{2}^{(7)},e^{-i(\delta_{l}+\varphi)}h_{1}^{(7)}\\},$
(62)

where $h_{1}^{(7)}$, $h_{2}^{(7)}$ are given by Eq.(46) and $I_{3}$ is 3x3
identity matrix. As previously observed, the phases
$\delta=(0,\pi),\delta_{l}$ and $\varphi$ depend on the group and its
representations.The resulting neutrino mass matrix Eq.(17) has the form;

$\displaystyle M^{\nu}=$
$\displaystyle\frac{g}{2M}\left(x|v_{1}|^{2}e^{2i\varphi_{1}}h_{2}^{(7)}+y|v_{1}v_{2}|e^{i(\varphi_{1}+\varphi_{2})}I_{3}\left(1+e^{i\delta}\right)\right.$
(63)
$\displaystyle+\left.x|v_{2}|^{2}e^{i(\delta+2(\varphi_{2}+\varphi))}h_{1}^{(7)}\right),$

which gives the squares of neutrino masses:

$\displaystyle m_{1}^{2}=$
$\displaystyle|c_{\nu}|^{2}\left[1+\kappa^{4}+2\beta^{2}\kappa^{2}+2\beta^{2}\kappa^{2}\cos(\delta)+4\beta\kappa^{3}\cos\left(\frac{\delta}{2}\right)\cos\left(\frac{\delta}{2}+\tau+2\varphi\right)\right.$
(64)
$\displaystyle\left.+4\beta\kappa\cos\left(\frac{\delta}{2}\right)\cos\left(\frac{\delta}{2}+\tau\right)+2\kappa^{2}\cos(\delta+2\tau+2\varphi)\right],$
$\displaystyle m_{2}^{2}=$
$\displaystyle|c_{\nu}|^{2}\left[1+\kappa^{4}+2\beta^{2}\kappa^{2}+2\beta^{2}\kappa^{2}\cos(\delta)+4\beta\kappa^{3}\cos\left(\frac{\delta}{2}\right)\sin\left(\frac{\delta}{2}+\tau+2\varphi-\frac{\pi}{6}\right)\right.$
(65)
$\displaystyle\left.+4\beta\kappa\cos\left(\frac{\delta}{2}\right)\sin\left(\frac{\delta}{2}+\tau-\frac{\pi}{6}\right)-2\kappa^{2}\sin\left(\delta+2\tau+2\varphi+\frac{\pi}{6}\right)\right],$
$\displaystyle m_{3}^{2}=$
$\displaystyle|c_{\nu}|^{2}\left[1+\kappa^{4}+2\beta^{2}\kappa^{2}+2\beta^{2}\kappa^{2}\cos(\delta)-4\beta\kappa^{3}\cos\left(\frac{\delta}{2}\right)\sin\left(\frac{\delta}{2}+\tau+2\varphi+\frac{\pi}{6}\right)\right.$
(66)
$\displaystyle\left.-4\beta\kappa\cos\left(\frac{\delta}{2}\right)\sin\left(\frac{\delta}{2}+\tau+\frac{\pi}{6}\right)+2\kappa^{2}\sin\left(\delta+2\tau+2\varphi-\frac{\pi}{6}\right)\right],$

where $\beta=y/x$ , $\tau=\varphi_{2}-\varphi_{1}$ and
$c_{\nu}=gx|v_{1}|^{2}e^{2i\varphi_{1}}/(2M)$.

The masses of charged leptons are given by the same formulae as given above
Eqs.(57-59) and it is impossible to fit electron, muon and tau lepton masses.
Our mass matrices for neutrinos and charged leptons are diagonal, so we also
do not have the ability to fit the mixing matrix.

We also see that the symmetry condition for Majorana neutrinos (Eq.18) does
not give any new flavour symmetry group with new three dimensional
representation $A_{L}$ which would not be present in the symmetry equations
for charged leptons (Eq.II.1). The Yukawa matrices $h_{1}^{(l)}$ and
$h_{2}^{(l)}$ in Eq.62 for charged leptons are exactly the same as in the case
of Dirac neutrinos (Eq.48). This is not a good conclusion. In the set of
groups we are consider, regardless of the adopted neutrino sector, we will not
find a symmetry that gives acceptable solutions for mass of charged leptons.

##  IV IV. Conclusions

We have explored the possibility of using some discrete flavour symmetry to
explain the masses and mixing matrix elements of leptons in the Standard Model
with the Higgs particle sector extended by one additional Higgs doublet -
2HDM. In general, we have assumed that the total Lagrangian model has a full
flavour symmetry with one exception. We also admit that the Higgs potential in
our model is only invariant in form, which is often used in the model
description of experimental data. In such a model, we have avoided having to
break the family symmetry and introduce flavon fields.
We have investigated discrete groups that are subgroups of the continuous
$U(3)$ group up to the order of 1025. Models in which neutrinos have the
nature of Dirac particles and models with Majorana neutrinos have been
considered. Following a close analysis of these 413 groups and all their
possible combinations of 2- and 3-dimensional irreducible representations, it
is established that none of them can reproduce the current experimental data.
Thus, in the 2HDM model with the symmetrical or only form-invariant Higgs
potential, in the class of the groups under consideration, there is no
discrete family symmetry that would fully clarify the masses and parameters of
the mixing matrix for leptons. In addition, we have observed that the set of
symmetric Yukawa matrices for charged leptons is independent from the nature
of the neutrinos. This can serve as a guidance for further search for family
symmetry. In the models with two Higgs particles, regardless of the adopted
neutrino sector and the set of groups that we consider, we do not find a
symmetry that gives real masses of charged leptons, even approximately.

###### Acknowledgements.

This work has been supported by the Polish National Science Centre (NCN) under
grant no. UMO - 2013/09/B/ST2/03382

## References

  * Ishimori  _et al._ (2010) H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada, and M. Tanimoto, Prog. Theor. Phys. Suppl. 183, 1 (2010), arXiv:1003.3552 [hep-th] .
  * Blum  _et al._ (2008) A. Blum, C. Hagedorn, and M. Lindner, Phys. Rev. D77, 076004 (2008), arXiv:0709.3450 [hep-ph] .
  * Frampton and Kephart (1995) P. H. Frampton and T. W. Kephart, Int. J. Mod. Phys. A10, 4689 (1995), arXiv:hep-ph/9409330 [hep-ph] .
  * King (2017) S. F. King, Prog. Part. Nucl. Phys. 94, 217 (2017), arXiv:1701.04413 [hep-ph] .
  * Altarelli, Guido and Feruglio, Ferruccio (2010) Altarelli, Guido and Feruglio, Ferruccio, Rev. Mod. Phys. 82, 2701 (2010), hep-ph/1002.0211 .
  * King and Luhn (2013) S. F. King and C. Luhn, Rept. Prog. Phys. 76, 056201 (2013), arXiv:1301.1340 [hep-ph] .
  * Harrison  _et al._ (2002) P. F. Harrison, D. H. Perkins, and W. G. Scott, Phys. Lett. B530, 167 (2002), arXiv:hep-ph/0202074 [hep-ph] .
  * Harrison and Scott (2003) P. F. Harrison and W. G. Scott, Phys. Lett. B557, 76 (2003), arXiv:hep-ph/0302025 [hep-ph] .
  * de Gouvea and Murayama (2015) A. de Gouvea and H. Murayama, Phys. Lett. B747, 479 (2015), arXiv:1204.1249 [hep-ph] .
  * Hall  _et al._ (2000) L. J. Hall, H. Murayama, and N. Weiner, Phys. Rev. Lett. 84, 2572 (2000), arXiv:hep-ph/9911341 [hep-ph] .
  * Giunti and Kim (2007) C. Giunti and C. W. Kim,  _Fundamentals of Neutrino Physics and Astrophysics_ (2007).
  * Lam (2008a) C. S. Lam, Phys. Rev. D78, 073015 (2008a), arXiv:0809.1185 [hep-ph] .
  * Lam (2008b) C. S. Lam, Phys. Rev. Lett. 101, 121602 (2008b).
  * Pontecorvo (1968) B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968), [Zh. Eksp. Teor. Fiz.53,1717(1967)].
  * Maki  _et al._ (1962) Z. Maki, M. Nakagawa, and S. Sakata, Progress of Theoretical Physics 28, 870 (1962).
  * Dziewit  _et al._ (2013) B. Dziewit, S. Zajac, and M. Zralek, Acta Phys. Polon. B44, 2353 (2013).
  * King and Luhn (2009) S. F. King and C. Luhn, JHEP 10, 093 (2009), arXiv:0908.1897 [hep-ph] .
  * King  _et al._ (2014) S. F. King, A. Merle, S. Morisi, Y. Shimizu, and M. Tanimoto, New J. Phys. 16, 045018 (2014), arXiv:1402.4271 [hep-ph] .
  * Machado  _et al._ (2011) A. C. B. Machado, J. C. Montero, and V. Pleitez, Phys. Lett. B697, 318 (2011), arXiv:1011.5855 [hep-ph] .
  * Siyeon (2014) K. Siyeon, J. Korean Phys. Soc. 65, 1347 (2014), arXiv:1208.2645 [hep-ph] .
  * Nebot (2016) M. Nebot, Nucl. Part. Phys. Proc. 273-275, 1448 (2016), arXiv:1410.7759 [hep-ph] .
  * Kubo  _et al._ (2003) J. Kubo, A. Mondragon, M. Mondragon, and E. Rodriguez-Jauregui, Prog. Theor. Phys. 109, 795 (2003), [Erratum: Prog. Theor. Phys.114,287(2005)], arXiv:hep-ph/0302196 [hep-ph] .
  * Nguyen (2014) T. Nguyen,  _Communications in Physics_ ,  24, 113 (2014).
  * Morisi and Peinado (2009) S. Morisi and E. Peinado, Phys. Rev. D80, 113011 (2009), arXiv:0910.4389 [hep-ph] .
  * Hou (1992) W.-S. Hou, Phys. Lett. B296, 179 (1992).
  * Branco  _et al._ (2012) G. Branco, P. Ferreira, L. Lavoura, M. Rebelo, M. Sher, and J. P. Silva, Physics Reports 516, 1 (2012), theory and phenomenology of two-Higgs-doublet models.
  * Davidson and Haber (2005) S. Davidson and H. E. Haber, Phys. Rev. D 72, 035004 (2005).
  * Ludl (6 27) P. O. Ludl,  _Systematic analysis of finite family symmetry groups and their application to the lepton sector_ , Ph.D. thesis, Vienna U. (2009-06-27), arXiv:0907.5587 [hep-ph] .
  * Ludl (2010) P. O. Ludl, J. Phys. A43, 395204 (2010), [Erratum: J. Phys.A44,139501(2011)], arXiv:1006.1479 [math-ph] .
  * Parattu and Wingerter (2011) K. M. Parattu and A. i. e. i. f. Wingerter, Phys. Rev. D 84, 013011 (2011).
  * Grimus and Ludl (2012) W. Grimus and P. O. Ludl, J. Phys. A45, 233001 (2012), arXiv:1110.6376 [hep-ph] .
  * (32) GAP, “GAP – Groups, Algorithms, and Programming, Version 4.7.6,” http://www.gap-system.org (2014).
  * Besche  _et al._ (2014) H. U. Besche, B. Eick, and E. O’Brien, “small, the Small Groups library, Version 2.1,” http://www.icm.tu-bs.de/ag_algebra/software/small/ (2014), GAP 4.7.6 component.
  * Dabbaghian (2011) V. Dabbaghian, “Repsn, a GAP4 Package for constructing representations of finite groups, Version 3.0.2,” http://www.sfu.ca/~vdabbagh/gap/repsn.html (2011), refereed GAP package.