# Robust Optimization of a Thermal Anti-Ice Protection System in Uncertain Cloud Conditions

Mariachiara Gallia, Bárbara Arizmendi Gutiérrez, Giulio Gori, Alberto Guardone, Pietro Marco Congedo

► To cite this version:

Mariachiara Gallia, Bárbara Arizmendi Gutiérrez, Giulio Gori, Alberto Guardone, Pietro Marco Congedo. Robust Optimization of a Thermal Anti-Ice Protection System in Uncertain Cloud Conditions. Journal of Aircraft, 2024, 61 (1), pp.1-15. 10.2514/1.C037223. hal-04346151

HAL Id: hal-04346151

https://inria.hal.science/hal-04346151v1

Submitted on 15 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Robust Optimization of a Thermal Anti-Ice Protection System in Uncertain Cloud Conditions

Mariachiara Gallia*, Bárbara Arizmendi Gutiérrez†, Giulio Gori ‡ and Alberto Guardone §

Department of Aerospace Science and Technology, Politecnico di Milano, Polytechnic Institute of Milan, Milano, 20156, Italy

Pietro M. Congedo¶

Platon Team, Centre de Mathématiques Appliquées, Ecole Polytechnique, Polytechnic Institute of Paris, IPP, Route de Saclay, 91128 Palaiseau, France

We present a framework for the robust optimization of the heat flux distribution for an anti-ice electro-thermal ice protection system (AI-ETIPS) under uncertain conditions. The considered uncertainty regards a lack of knowledge concerning the characteristics of the cloud i.e. the liquid water content and the median volume diameter of water droplets, and the accuracy of measuring devices i.e., the static temperature probe. Uncertain parameters are modelled as random variables and two sets of bounds are investigated. A forward uncertainty propagation analysis is carried out using a Monte Carlo approach exploiting a surrogate models. The optimization framework relies on a gradient-free algorithm (Mesh Adaptive Direct Search) and two different objective functions are considered. Namely, the 95 quantile of the freezing mass rate and the statistical frequency of the fully evaporative operating regime. The framework is applied to a reference test case, revealing a potential to improve the heat flux distribution of the baseline design. A new heat flux distribution is proposed that presents a more efficient use of the thermal power, increase flight safety even at non-nominal environmental conditions.

*

## Nomenclature

### Greek Letters

$ρ$ = Density [kg m⁻³]

### Parameters

<footnote>*PhD Student, mariachiara.gallia@polimi.it, AIAA University Student (Corresponding Author)</footnote>

<footnote>†Post-doctoral Researcher, barbara.arizmendi@polimi.it</footnote>

<footnote>‡Assistant Professor, giulio.gori@polimi.it</footnote>

<footnote>§Full Professor, alberto.guardone@polimi.it</footnote>

<footnote>¶Research Director, pietro.congedo@inria.fr</footnote>

<footnote>*The present work has been presented as an oral only presentation at the AIAA Aviation 2021 Forum, held Online August 2-6, 2021, https://doi.org/10.2514/MAVIAT21</footnote>
$c_p$ = Specific heat at constant pressure [J kg⁻¹ K⁻¹]

$\Delta l$ = Heater size [m]

$\dot{H}''$ = Vector of heat fluxes across the heaters [W m⁻²]

$k$ = Optimization algorithm iteration counter

$\text{LWC}$ = Liquid water content [g m⁻³]

$\dot{m}$ = Local mass flow rate [kg s⁻¹]

$\dot{M}$ = Total mass flow rate [kg s⁻¹]

$\dot{M}$ = Vector of total mass flow rate across the two airfoil sides [kg s⁻¹]

$||\dot{M}_{\text{ice}}||_2$ = Performance metric of runback ice formations [kg s⁻¹]

$\dot{M}_{\text{ice}}^{95}$ = Robust performance metric of runback ice formations [kg s⁻¹]

$\text{MVD}$ = Median volume diameter [m]

$N$ = Number of

$P$ = Probability

$P_b$ = Thermal power budget [W]

$\dot{Q}$ = Heat flow rate [W]

$s$ = Curvilinear abscissa [m]

$\text{SAT}$ = Static air temperature [K]

$V$ = Freestream air velocity [m s⁻¹]

$\mathbf{x}$ = Flight properties

$\xi$ = Uncertain cloud properties

## Subscripts

$\text{air}$ = Air

$g$ = Control elements

$h$ = Heater

$\text{ice}$ = Ice

$\text{IPS}$ = Ice protection system

$\text{lb}$ = Lower bound

$\text{MC}$ = Monte Carlo samples

$\text{train}$ = Samples for training

$\text{ub}$ = Upper bound

$\infty$ = Freestream

## Superscripts
95 = 95 quantile

best = Best design found so far

# I. Introduction

In-flight icing poses a severe threat to flight safety, unpredictably hindering the aircraft's performance. Electro-Thermal Ice Protection System (ETIPS) is a mature technology commonly preferred for mitigating the risk of in-flight icing. Thanks to their modularity, ETIPS are one of the most flexible technologies for the design and operation of Ice Protection Systems (IPS); e.g., ETIPS modules can be activated differently and independently, depending on the requirements. However, these systems require a tremendous amount of power (approx. 10 kW for a single-aisle commercial aircraft) to protect a large surface such as an aircraft wing. For instance, when operating under Anti-Ice (AI) requirements, the power absorbed by the protection system must be sufficient to ensure that all the impinging water is evaporated to avoid any runback liquid water freezing in unprotected parts, thus compromising the performance and safety of aircraft [1]. During a long icing event, the energy consumption to protect a large surface is enormous. In addition, real applications are unavoidably affected by uncertainty challenging the optimal operation of IPS. Indeed, the exact values of many characteristic cloud parameters, e.g., the droplets' Median Volume Diameter (MVD) or the Liquid Water Content (LWC), are generally unknown during flight because their measurement requires dedicated probes that are deployed only during test campaigns. Typically, the only parameter (concerning icing) known during actual flight operations is the static air temperature (SAT), measured through onboard devices. Nonetheless, as with any measurement, SAT readings are subject to uncertainty arising from the accuracy of the probes. It follows that, generally, icing events occur in scarcely known conditions; therefore, the design of protection systems requires a systematic treatment of the uncertainties. As a consequence, achieving the optimal design of IPS entails an inherent need to test a large number of different configurations. Experimental test campaigns would be just too costly and time-consuming, also considering the limited availability of test facilities and the cost of in-flight testing, which is approximately € 3-5M for Appendix C conditions and € 10-20M in Appendix O conditions, as defined in [2]. For this reason, numerous works in literature deal with the formulation of numerical models for the investigation and optimization of AI-ETIPS e.g., ANTICE [3], CANICE [4]. In particular, some works target the optimization of the layout of thermal sources. For instance, Pourbagian and Habashi [5] and Pourbagian et al. [6] performed a numerical optimization of an AI-ETIPS protecting an airfoil. The authors optimized the heat flux distribution across a predetermined design as well as the sizes of the heating bands employing the Mesh Adaptive Direct Search algorithm. The authors considered known and fixed cloud and flight properties in their study. They found potential efficiency improvements compared with an intuitive allocation of heat fluxes. On a similar note, Arizmendi et al. [7] found potential thermal power savings in the range of 10%, with no runback ice formations, thanks to an optimized layout of heat fluxes obtained through a genetic
algorithm. Moreover, there are additional optimization studies in the field of thermal-pneumatic IPS which investigate its geometrical parameters, such as the works of Saeed and Paraschivoiu [8] and Pellisier et al. [9].

In all the referenced works, it is assumed that cloud characteristics are deterministically known, while they are commonly uncertain (or even unknown) in real applications. This jeopardizes the consistency of the results and calls for systematically including the inherent uncertainty. Uncertainty Quantification (UQ) for in-flight icing applications is the subject of several works found in the literature. Existing studies focus on the probabilistic prediction of the severity category of icing encounters for uncertain environmental parameters in unprotected components. For instance, Zhang et al. [10] developed an analytical model to assess the ice formation rate accounting for operational environmental uncertain parameters, which were modeled as random variables. The range of values of each severity category was also considered uncertain, and a belonging class function was defined. The uncertainties were propagated employing Monte Carlo Sampling. Then, in an analogous study, Feng et al. [11] analytically propagated the cloud uncertainties to quantify the probability of belonging to a severity category defined in the work of Zhang et al. [10]. Later, Zhang et al. [12] performed a sensitivity analysis of ice severity on the environmental parameters through fuzzy-state assumption and the Profust model. Furthermore, De Gennaro et al. [13] conducted a UQ study of the effects entailed by the position and size of standard ice shapes on the performance degradation, using an approach based on Polynomial Chaos expansions. Gori et al. [14] conducted a study to understand and quantify the effects of parameter uncertainties on the numerical modeling of in-flight icing shapes, proposing a non-linear approach for the propagation of the uncertainties. To the authors' knowledge, no applications of Uncertainty Propagation techniques for the analysis and design of Thermal IPS are available in the open literature.

This work deals with the development of a robust optimization framework for the management of the onboard power required for operating an AI-ETIPS to protect the aircraft during a stratiform cloud encounter in uncertain icing conditions. Namely, a configuration is sought that minimizes the convective losses, i.e., the heat provided by the protection system but wasted because of the cooling effect of the external airflow, while maximizing the heat devoted to the evaporation of the liquid water film. The considered AI-ETIPS includes several modules capable of delivering arbitrary thermal power. The design parameters consist of the heat fluxes delivered by each module, and the optimal power distribution is sought under the constraint of maintaining an overall total power consumption. A computational Robust Design Optimization (RDO) framework is developed to execute the task. In RDO, the minimization of some statistics of the targeted Quantity of Interest (QoI) is sought. Therefore, a metric of robustness must be formulated. The pioneer of robustness in product design was Taguchi who, in the 1980s, envisioned the concept of robustness as a desired product feature even in their early stages [15]. He proposed quantifying deviations from target performance due to those noise factors, namely aleatory uncertainties. In this way, the consistency in performance can be quantified, seeking narrow performance probability distributions with predictable performance. Numerous RDO applications deal with the design of aeronautical systems such as an internal combustion engines [16], control systems [17], and
aerodynamic shapes, both clean [18] and iced ones [19]. In this work, we first perform the deterministic optimization of the IPS system to obtain a benchmark case. Afterward, we perform the robust optimization of the IPS w.r.t. the uncertain cloud parameters. The robust optimization is achieved considering two different cost functions, each entailing a specific statistic concerning the Quantity of interest. Namely, in the first case, we aim to reduce the risk of severe ice formation by minimizing a quantile of the ice accretion rate. In the second case, we also seek consistency in the fully evaporative operation of IPS outside the nominal design conditions. The performance of the robust designs is then compared to the baseline and the benchmark configuration. In all cases, during the optimization it is considered a fixed power budget that equals the baseline. This addresses an important constraint for aircraft design on the limited onboard power available. It is important to highlight that the proposed framework is highly flexible. Different objective functions, and constraints can be considered as well as different uncertain variables, uncertainty bounds and distributions. In the test case here presented, the QoI was selected to answer the question of whether the configuration presented in the baseline could be more efficient in evaporating impinging water by modifying the heat fluxes layout considering uncertain cloud conditions. Other questions could be easily addressed by choosing a different problem formulation.

The paper is structured as follows. In Section II, the case of study is introduced along with the cloud uncertainty characterization. After, in Section III, the mathematical formulation of the problem is presented. Next, the methodology for solving the problem is reported (Section IV), describing the numerical model, the optimization algorithm, the uncertainty propagation procedure, and the complete robust design optimization framework deployed. The results are discussed in Section V, and finally, the conclusions of the work are drawn in Section VI.

## II. Case of Study

### A. Geometric and Operational Parameters

The geometry considered is an extruded NACA0012 profile with a chord of 0.9144 m with seven heaters placed at the leading edge. The layout is taken from the experimental work of Al-Khalil et al. [20]. The IPS comprised a set of 7 multilayered heaters fitted at the leading edge and expanding span-wise. In the geometry, depicted in Fig. 1, due to a manufacturing issue, the heaters were shifted towards the airfoil suction side with an offset of 0.0145 m. Therefore, the configuration is asymmetric. Ultimately, a 2D profile is considered in this work, assuming that the three-dimensional effects are negligible on a straight geometry. For computational reasons, which will be explained in section IV, the airfoil is divided into two parts: the suction side (SS) and pressure side (PS), as shown in Fig. 1. The index to identify these two sides is $j = [SS, PS]$. From the experimental test cases reported in the work of Al-Khalil et al. [20], experiment 67A is selected as the baseline case for comparison with the optimal designs. The flight conditions, assumed to be fixed and known, are the freestream velocity (89.4 ms⁻¹), the pressure freestream static (90 000 Pa), and the angle of attack (0°). As for the uncertain cloud properties, they are presented in Section II.B. Concerning the AI-ETIPS, the baseline
<figure><img src="image_2.png" /><figcaption>Fig. 1 layout of the heaters of the ETIPS.</figcaption></figure>

operational heat flux distribution and size are reported in Table 1, while the total power consumption per unit span is $P_b = 4815 \text{ W m}^{-1}$. This configuration was designed to operate in fully evaporative regime for specific cloud conditions

<table><caption>**Table 1** Baseline operational heat flux distribution across the heaters and heaters size.</caption><thead><tr><th></th><th>heater G</th><th>heater E</th><th>heater C</th><th>heater A</th><th>heater B</th><th>heater D</th><th>heater F</th></tr></thead><tbody><tr><td>Heat flux [Wm<sup>-2</sup>]</td><td>18 600</td><td>18 600</td><td>26 350</td><td>43 400</td><td>32 550</td><td>21 700</td><td>20 150</td></tr><tr><td>Δl<sub>i</sub> [m]</td><td>0.038</td><td>0.025</td><td>0.025</td><td>0.019</td><td>0.025</td><td>0.025</td><td>0.038</td></tr></tbody></table>

having MVD = 20 µm and LWC = 0.55gm⁻³. The allocation of heat fluxes was performed intuitively, meaning that the largest heat fluxes were allocated on the heaters positioned where the largest amount of impinging water is predicted. This layout and operation conditions are selected as the baseline configuration since they are realistic and provide a reasonable benchmark for assessing the results obtained in this work.

## B. Cloud Uncertainty Characterization

Icing events associated with stratiform clouds are generally mild but can extend horizontally over many kilometers. The characteristic values of LWC and MVD are low compared to cumuliform clouds encounter, but the long exposure time produces a slow but potentially severe ice build-up. Two primary sources of uncertainty have been identified in the quantification of the cloud parameters. First, cloud properties cannot be controlled and vary over space and time. Unfortunately, the prediction of these parameters is complex, and the accuracy of remote measurement from weather stations limited [21]. Moreover, on-board ice detection systems only deliver qualitative information about the encounter. Consequently, the severity of the encounter is generally not quantifiable. Secondly, some icing-related quantities, e.g., the SAT, are measured using on-board probes. Nevertheless, these readings are subject to uncertainties due to the probe accuracy, which, for SAT, typically lies in the range of ±2°C.

Given the complexity of retrieving information about cloud properties, these are modeled as random and independent variables uniformly distributed within a min-max range. This choice limits the amount of information introduced by the
modeler regarding the random parameter. Indeed, only an upper and a lower limiting bound are introduced for each parameter, being any value included in the range equally probable. This approach is general and conservative because, we are assuming a larger envelope compared to, for example, Appendic C envelope. Indeed, in the envelope considered in this work conditions that may not occur or rarely occur in cloud encounters are considered like high LWC at low SAT. Our approach is overly-conservative since experimental flight campaigns revealed that lower LWC values are significantly more frequent [22, 23], especially at lower temperature. Concerning the SAT, we alternatively assume two different uncertain ranges B1 and B2. Namely, B1 emulates the limited knowledge about the time/space variation of stratiform icing clouds in the absence of direct real-time measurements. Therefore, the possible SAT temperature ranges between reference values specified in current regulation for, for in-flight icing, namely the values associated with stratiform clouds reported in Appendix C of the title 14 Code of the Federal Regulation part 25 for Continuous Maximum Events [2]. The range B2 emulates the limited accuracy of real-time measurements instead, assuming that a temperature probe is available on board the aircraft. Hence, the limiting bounds for B2 are established from the analysis of existing measurement systems, usually employed in experiments or on board. Note that the bounds selected for the LWC and MVD are equal in both sets because, in either case, their values are unknown. All the values considered are reported in Table 2.

<table><caption>**Table 2** Characteristic bounds of uncertain cloud parameters for lack of control of the variations (B1) and the limited precision on the measurement of the SAT (B2)</caption><thead><tr><th></th><th>LWC [gs<sup>-1</sup>]</th><th>MVD [µm]</th><th>SAT-B1 [K]</th><th>SAT-B2 [K]</th></tr></thead><tbody><tr><td>Lower Bound</td><td>0</td><td>10</td><td>238.15</td><td>249.27</td></tr><tr><td>Upper Bound</td><td>1</td><td>50</td><td>273.15</td><td>253.27</td></tr></tbody></table>

### III. Problem Formulation

The vector of design parameters $\dot{H}'' = \{\dot{H}''_i\} \in \mathbb{R}^{N_h}$, with $N_h = 7$, includes the heat fluxes $\dot{H}''_i$ allocated to each heater $i$. The layout of the heating system is depicted in Fig. 1. Each heater $i$ has a specific size $\Delta l_i$, corresponding to its spatial extension along the airfoil curvilinear abscissa reported in Table 1. The baseline heat flux profile, as well as the total power budget per unit span $P_b = 4815 \text{ W m}^{-1}$, are extracted from the test case 67A from the experimental campaign conducted by Al-Khalil et al. [20], that has been presented in Section II. Note that additional constraints apply to this problem. Namely, a maximum and a minimum heat flux value of, respectively, $H_{ub} = 46 000 \text{ W m}^{-1}$ and $H_{lb} = 4325 \text{ W m}^{-1}$, must be enforced to each heater.

The Quantity of Interest (QoI) targeted by the optimizer is the water freezing mass rate $\dot{M}_{ice,j}$, which is computed independently at the different locations (pressure and suction side of the airfoil). Due to the uncertain cloud characteristics
described in Sec. II.B, included here in vector $\xi$, the targeted QoI is a stochastic quantity, and it is defined as

$$ \dot{M}_{\text{ice},j}(\mathbf{x}, \xi, \dot{\mathbf{H}}'') = \sum_{n=1}^{N_g} \dot{m}_{\text{ice},j}(\mathbf{x}, \xi, \dot{\mathbf{H}}'', s_n), \quad (1) $$

where $\mathbf{x}$ are the flight conditions, $\xi$ are the uncertain environmental conditions, $s_k$ indicates the discretized curvilinear coordinate along the airfoil, and $N_g$ is the number of control elements on one side of the airfoil. This quantity gives a reasonable estimation of the ice accretion thickness and mass on the airfoil, which is commonly used for aircraft certification. Indeed, considering that $\dot{M}_{\text{ice},j}$ remains constant over time in anti-ice operations due to the long operational time that leads to steady-state conditions, it is linearly related to those two quantities. It is important to highlight that generally, for design purposes, an anti-ice IPS in fully evaporative condition should guarantee no ice formation in the complete App. C envelop. In order to guarantee this condition, for the baseline test case, a power increase would be required. Nonetheless, the aim of this work is to maintain the power consumption constant and only modify the heat flux distribution to reduce the probability and severity of ice formation on the surface. Therefore, run-back ice formations are sought to be minimized. As shown in Fig 1, the airfoil is divided in two sides $j = [SS, PS]$ and the QoI is computed separately for each side as in Eq. 1. Thus, among all possible configurations with the same overall rate of ice accretion, the ones with symmetric ice formations are selected during the optimization procedure. This is done to avoid zero accretion on one side and high on the other side, which would cause greater performance degradation. The optimization problem must be solved under the constraint that the overall power consumption is below the $P_b$ threshold; otherwise, the optimal solution is deemed infeasible. Moreover, the optimization process seeks a distribution that maximizes the amount of heat devoted to evaporation to avoid water running to unprotected regions while minimizing convective losses. In this work, we first focus on a deterministic optimization of the AI-ETIPS, i.e., considering the uncertain cloud parameters as exactly known. After, we restore the uncertainty and carry out robust optimization, providing a comparison of the two solutions.
## A. Deterministic formulation

In a deterministic setting, the formulation of the optimization problem reads:

$$ \underset{\dot{H}'' \in \mathbb{R}^{N_h}}{\text{minimize}} \quad \|\dot{\mathbf{M}}_{\text{ice}}\|_2, \qquad (2a) $$

$$ \text{subject to} \quad P_b \ge \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \qquad (2b) $$

$$ \dot{H}_{N_h}'' = P_b - \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \qquad (2c) $$

$$ \dot{H}_i'' \le H_{ub}'' \text{ W m}^{-2}, \qquad (2d) $$

$$ \dot{H}_i'' \ge H_{lb}'' \text{ W m}^{-2}, \qquad (2e) $$

where $P_b$ is the power budget available per unit span, $\dot{H}_i''$ is the heat flux provided by heater $i$, $\Delta l_i$ is the length of heater $i$, and $\|\dot{\mathbf{M}}_{\text{ice}}\|_2$, which is the quantity of interest, is formulated by applying Eq. (1) on both sides of the airfoil and computing the $L_2$ norm:

$$ \|\dot{\mathbf{M}}_{\text{ice}}\|_2 = \sqrt{\dot{M}_{\text{ice,PS}}^2 + \dot{M}_{\text{ice,SS}}^2}. \qquad (3) $$

Note that, in the following, the deterministic problem considers a worst-case scenario for impinging water. In this case, the selected LWC is the value associate to the continuous maximum icing conditions from the reference literature namely, LWC = 1 g m⁻³. The reference SAT is 251.35 K, whereas MVD = 20 µm. Under these conditions, the baseline AI-ETIPS design operates in running wet regime, and the optimal layout minimizes run-back ice formations.

## B. Robust formulation

In the robust optimization setting, the uncertainty associated to cloud parameters $\xi$ is handled by establishing a stochastic QoI associated to $\dot{M}_{\text{ice}}$. By optimizing a statistic of the QoI it is possible to enforce a reduction of the variability of the desired performance, thus promoting a consistent behavior robust to uncertain cloud parameters. The choice of the statistic is fundamental, as it determines the robustness characteristics of the optimal solution, and different criteria may be applied to foster the desired behavior. In the following, we consider two different statistics of $\dot{M}_{\text{ice}}$: the minimization of the 95 quantile of $\dot{M}_{\text{ice},j}$, $\dot{M}_{\text{ice}}^{95}$, and the maximization of the success probability of no-ice formation $P(\dot{M}_{\text{ice}} = 0)$ w.r.t. a binomial distribution modeling the ice/no-ice formation.

In the following, the first robust optimization problem will be labeled as A. Minimizing the 95 quantile corresponds to ensuring optimal performance for 95% of the random events. It indicates the magnitude of the 5% most unfavorable cases for a fixed heat flux distribution and flight conditions. The selection of this metric and its later minimization reduces the probability of experiencing severe ice formations. For a single ice formation location $j$, the robustness
measure is equal to:

$$ \dot{M}_{\text{ice},j}^{95}(\mathbf{x}, \xi, \dot{\mathbf{H}}'') = q_{\xi}^{95}[\dot{M}_{\text{ice},j}](\mathbf{x}, \dot{\mathbf{H}}'') \quad (4) $$

The quantities in each location $j = [PS, SS]$ are aggregated in the vector $\dot{M}_{\text{ice}}^{95}$. The objective function to be minimized $f(\mathbf{x}, \xi, \dot{\mathbf{H}}'')$ is the $L^2$ norm of this vector in the Euclidean space $\mathbb{R}^2$. In this way, it is sought to find a design configuration that ensures an even performance across each location compared to considering them as a single location. The formulation for the robust optimization task entailing the minimization of the 95 quantile then reads:

$$ \underset{\dot{\mathbf{H}}'' \in \mathbb{R}^{N_h}}{\text{minimize}} \quad ||\dot{M}_{\text{ice}}^{95}||_2, \quad (5a) $$

$$ \text{subject to} \quad P_b \ge \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \quad (5b) $$

$$ \dot{H}_{N_h}'' = P_b - \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \quad (5c) $$

$$ \dot{H}_i'' \le H_{ub}'' \text{ W m}^{-2}, \quad (5d) $$

$$ \dot{H}_i'' \ge H_{lb}'' \text{ W m}^{-2}, \quad (5e) $$

where $P_b$ is the power budget available per unit span, $\dot{H}_i''$ is the heat flux provided by heater $i$, $\Delta l_i$ is the length of heater $i$, and $||\dot{M}_{\text{ice}}^{95}||_2$ is the 95 quantile of the rate of ice accretion. In case the probability of no-ice formations $\mathcal{P}(\dot{M}_{\text{ice}} = 0)$ over the whole airfoil is maximized, we will be referring to robust optimization B. This metric directly ensures an increased consistency of complete evaporation of the impinging water regardless of the cloud conditions. The robust formulation problem then reads:

$$ \underset{\dot{\mathbf{H}}'' \in \mathbb{R}^{N_h}}{\text{maximize}} \quad P(\dot{M}_{\text{ice}} = 0) \quad (6a) $$

$$ \text{subject to} \quad P_b \ge \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \quad (6b) $$

$$ \dot{H}_{N_h}'' = P_b - \sum_{i=1}^{N_h-1} \dot{H}_i'' \Delta l_i \text{ Wm}^{-2}, \quad (6c) $$

$$ \dot{H}_i'' \le \dot{H}_{ub}'' \text{ W m}^{-2}, \quad (6d) $$

$$ \dot{H}_i'' \ge \dot{H}_{lb}'' \text{ W m}^{-2}, \quad (6e) $$

$$ ||\dot{M}_{\text{ice}}^{95}||_2(\mathbf{x}, \xi, \dot{\mathbf{H}}'') \le ||\dot{M}_{\text{ice}}^{95}||_2(\mathbf{x}, \xi, \dot{\mathbf{H}}_0''), \quad (6f) $$

where $P_b$ is the power budget available per unit span, $\dot{H}_i''$ is the heat flux provided by heater $i$, $\Delta l_i$ is the length of heater $i$.
heater *i*, and $||\dot{M}_{\text{ice}}^{95}(x, \xi, \dot{H}'')||_2$ is the 95 quantile of the rate of ice accretion for design $\dot{H}''$, while $||\dot{M}_{\text{ice}}^{95}(x, \xi, \dot{H}''_0)||_2$ refers to the quantile of the baseline design. Note that, to limit edge cases of ice formations that would be undetected using this metric, the additional constraint Eq. 6f is added. That is, the performance of the baseline configuration $\dot{H}''_0$ considering the $L^2$ norm of the quantile must be at least equaled. This second problem formulation follows the results obtained from RA. Indeed, as it will be shown in Sec. V, the quantile optimization does not guarantee a high probability of no ice formation.

# IV. Methodology

## A. Numerical Modeling

To assess the quantities of interest presented in Section III, the PoliMIce framework [24] for in-flight ice accretion has been employed. The modeling of in-flight icing is concerned with describing a complex multi-phase physical process. Regarding our purposes, the task presents an additional complexity arising from the necessity of including the effects of the anti-ice protection system.

The PoliMIce simulation framework [24] was initially devised to investigate the in-flight ice accretion problem on two- and three-dimensional geometries. The software implements state-of-the-art icing models [25], including features specifically tailored to the simulation and robust design optimization of thermal IPS [7, 26]. Moreover, the PoliMIce framework has also been employed to carry out UQ analysis in icing conditions [14]. In the context of the 1st AIAA Ice Prediction Workshop, numerical predictions from the PoliMIce framework have been assessed against high-quality experimental measurements taken from the NASA Glenn Icing Research Tunnel. In general, it was shown that predictions compare favorably with the experimental measurements, see [27]. In particular, key droplet impingement characteristics, including the collection efficiency peak and impingement limits, are well captured.

In short, the PoliMIce toolkit consists of a computational framework loosely coupling three different tools. First, the aerodynamic field developing around the geometry of interest is reconstructed using the open-source Computational Fluid Dynamics (CFD) solver SU2 [28]. In this work, the aerodynamic model relies on the Euler flow assumptions. An a-posteriori correction is employed at later stages, i.e., within the thermodynamic solver, to retrieve the viscous effects associated with the boundary layer developing at walls and model the temperature gradients at the surface of the airfoil. This correction is performed through an integral boundary layer model. Momentum and thermal boundary layer equations are solved both for laminar and turbulent regimes, considering smooth transition when moving from one regime to the other The numerical grid discretizing the domain is generated using the in-house software uhMesh [29].

Once the aerodynamic solution is available, the second step consists of calculating the distribution of water collected over the body's surface, the so-called collection efficiency $\beta$, which is conducted thanks to a Lagrangian Particle Tracking solver named PoliDrop [30]. PoliDrop uses the aerodynamic CFD solution to compute the trajectories of
water droplets, ultimately locating their impacting points over the aircraft's surface. Trajectories are computed using a Lagrangian approach.

Finally, the third tool is a thermodynamic solver, i.e., the actual ice accretion engine, which solves the energy and the mass conservation equations describing a multi-phase, multi-layered stratum. The governing equations are solved in a discrete domain divided into control elements. Within each cell, the ice accretion engine models the solid substratum corresponding to the anti-ice system, the ice and the liquid water layers, and the outer air flow in the close proximity of the airfoil (or ice) surface, evaluating the heat transfer between layers to ultimately compute the thickness of the ice and the liquid water fraction. Since the anti-ice system substratum is modeled, it is possible to simulate the action of the protection system under different operating conditions as presented in [26]. Due to the limited thickness, the IPS substratum is reduced to a unique layer with an equivalent thermal conductivity and a fixed thickness. It is noted that there are more sophisticated conduction models available in the open literature, especially for de-icing applications solving the heat equation in the multilayered substratum, such as in the work of Reid et al. [31]. For anti-icing applications, the equivalent thermal conductivity approach is considered standard practice. Examples of this are the works of Silva et al. [32] and Bu et al. [33]. In this work, a steady-state is assumed because Anti-Ice devices typically deal with long exposures to icing conditions. These could happen for example when encountering a stratiform icing cloud, which tends to extend horizontally. For instance, Ref. [34] reports measurements of icing encounters extending for over 100 km.

## B. Optimization Algorithm

The chosen algorithm is the Mesh Adaptive Direct Search (MADS), developed by Audet and Dennis [35], and deployed in the open-source package NOMAD [36]. MADS is a direct method, meaning that it is a derivative-free algorithm and belongs to the Pattern Search algorithm family. These rely on function evaluations and sorting and, consequently, are suitable for complex functions, including non-linearities or singularities. Besides, as pointed out by Rios and Sahinidis [37], these algorithms are particularly ideal for Simulation-Based Optimization whenever the mathematical formulation of the objective function is unavailable, computationally expensive, or subject to numerical noise. For these problems, the computation of local gradients can be unfeasible or unreliable.

A thorough description of the algorithm is out of the scope of this work. The interested reader could consult the reference sources for further information on the algorithm and its proof of convergence [35, 36, 38]. Since the objective functions considered here present a complex behavior, including non-linearities, the MADS algorithm may converge to different local minima, depending on the starting point. To assess the consistency of the optimization and promote the convergence to the global optimum, ten runs are produced for each test case presented. Each run assumes a different and random initial point.
### C. Uncertainty Propagation

A Monte Carlo approach was employed to forward propagate the uncertainty affecting cloud parameters through the full computational model to ultimately characterize the probability distribution of the targeted QoI i.e., $\dot{M}_{ice,j}$. One of the main advantages of the Monte Carlo approach is that it can deal with discontinuous and irregular models. On the other hand, a drawback is that it requires many model evaluations or realizations to obtain reliable statistical quantities estimates. In the present work, a single evaluation of the icing model entails the execution of an Euler CFD simulation, a particle tracking computation, and thermodynamic analysis of the multi-phase problem over the geometry of interest. Overall, the completion of these tasks requires about 1 h if executed on a single cluster node†. To achieve an adequate convergence of the estimated probability distribution, at least $N_{MC} = 1 \times 10^3$, or even more for low error, Monte Carlo samples input must be evaluated, i.e., a total of $1 \times 10^3$ h or about 42 days if run serially on a single cluster node. This would then be repeated for each design point of the optimization problem. Due to the limited computational resources available, such computational time would be prohibitive. Therefore, a Gaussian Process surrogate of the PoliMIce framework was built to perform the uncertainty propagation analysis.

### D. Surrogate Model

The main issue concerning building an accurate surrogate is capturing non-smooth features, e.g., discontinuities, associated with predictions from the original model. For instance, increasing the SAT causes $T_{rec}$ to reach values higher than the freezing point. Consequently, ice does not form regardless of the value of LWC and MVD. This results in a non-smooth behavior of the quantity of interest between different portions of the uncertainty input space. Depending on the particular input combination, $\dot{M}_{ice,j}$ can either assume a value equal or larger than zero. The large number of zero values in the training data set's output will lead to poor overall quality regression. For this reason, a standard Gaussian Process surrogate performs quite unsatisfactory. A combined classification plus regression Gaussian Process has been implemented to improve the accuracy of the surrogate. Namely, inputs resulting into $\dot{M}_{ice,j} = 0$ has been classified as class 0, whereas inputs resulting in a non-zero value of $\dot{M}_{ice,j}$ as class 1. The model takes as inputs the uncertain cloud properties $\xi(\omega)$ and returns the freezing mass rate on one side of the airfoil $\dot{M}_{ice,j}$ for a fixed heat flux distribution. The surrogate model is divided into two parts. First a binary classification is performed: ice formation $\dot{M}_{ice,j} \neq 0$ (class 1) or no ice formation $\dot{M}_{ice,j} = 0$ (class 0). Based on the posterior probability $P(1|x)$, the data is classified as:

$$ \text{Class} = \begin{cases} 1 & \text{if } P(1|x) > 0.5 \\ 0 & \text{if } P(1|x) \le 0.5 \end{cases} \qquad (7) $$

Then, if the data is classified in class 1, the freezing mass rate is determined through the regression surrogate model.

<footnote>†the single node includes two processors Intel Xeon X-5650 2.67 GHz, for a total of 12 cores</footnote>
The python library GPy [39] has been used to build a surrogate model based on Gaussian Processes. Several trials have been conducted to determine the most suitable kernel for classification and regression to reduce the prediction error. A composed kernel was used for classification by adding a white noise kernel and an Ornstein-Uhlenbeck exponential kernel. While for regression, an RBF kernel was chosen. The hyperparameters were then optimized by maximizing the likelihood of the observations. Before performing the uncertainty quantification through Monte Carlo sampling, the goodness of the surrogate models was tested against a random validation data set of 250 samples. The surrogate was trained with an increasing number of training samples drawn from a space-filling Halton sequence ranging from 10 to 250 samples. To assess the accuracy of the surrogate, the mean error has been computed by averaging among all test samples. Fig. 2 shows the mean surrogate error ($\frac{\sum_{i=1}^{N_{\text{samples}}} (\dot{M}_{\text{ice,surrogate},j}^i - \dot{M}_{\text{ice},j}^i)}{N_{\text{samples}}}$) and the percentage classification error ($\frac{N_{\text{misclassified}}}{N_{\text{samples}}} \times 100$) by varying the number of training samples. Moreover the precision, recall and F-Score for classification have been evaluated. Fig. 2 on the left shows the quality indexes for the model built for the data set with larger temperature bounds B1, while Fig. 2 on the right shows the quality indexes for the surrogate regarding the data set with smaller temperature bounds B2. The quality indexes for only regression and classification plus regression surrogate model are shown and compared in both cases.

One can notice how the mean surrogate error is higher when only regression is performed for both bounds. Furthermore, the classification error, i.e., the number of misclassified samples, is considerably larger when only the regression model is applied. The same consideration applies for the precision, recall, and F-Score where for the classification plus regression surrogate model the F-Score value is close to one, as it should be, while this value is not reached for the surrogate model built using only regression. Overall the classification plus regression surrogate is better performing. The surrogate model built for broader temperature bounds data set B1 is not as accurate as the B2 data set regarding classification. The narrower the uncertain temperature bounds considered, the more accurate the classification model. By analyzing the samples, most of the misclassified test data are close to the temperature upper bound, where there is no ice formation regardless of LWC and MVD values. This leads to poorer surrogate performance for data set B1 compared to data set B2. Finally, a reliable surrogate model with a limited computational cost is required to perform the uncertainty propagation. From the errors presented in Fig. 2 it is seen that for data set B1, the error remains quite constant after 100 training samples, while for data set B2 after 50 training samples. For the UQ the model has been trained with $N_{\text{Train}} = 100$ for B1 and $N_{\text{Train}} = 50$ for B2.

## E. Robust Design Optimization

Next, the integration of the optimization process is shown in Fig. 3 for the optimization of $||\dot{M}_{\text{ice}}^{95}||_2$. The optimization routine for $P(\dot{M}_{\text{ice}} = 0)$ is analogous. The optimization algorithm uses the model as a black-box function to compute the performance of each solution $H_j''$, and it does not interact with it. The user provides an initial incumbent point. A set of $i = 1:N$ trial points is then generated by search and poll as presented in Section IV.B. Then, for each trial point,
Bound set B1 - Pressure Side

<figure><img src="image_3.png" /></figure>

Bound set B2 - Pressure Side

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2 Errors, precision, recall and F-Score of the surrogate model varying the number of training samples for B1 and B2 uncertain cloud properties bounds</figcaption></figure>
$N_{\text{train}}$ samples are extracted from the uncertain variable space through Latin hypercube sampling. These are passed to the PoliMIce framework to obtain $\dot{M}_{\text{ice},j}$ and then train the surrogate model, recalling that the index $j$ represents the two side of the airfoil where the quantity of interest is computed independently. Once the surrogate model is built for that design point, the Monte Carlo uncertainty propagation is performed for evaluating the targeted statistic ($||\dot{M}_{\text{ice}}^{95}||_2$ or $P(\dot{M}_{\text{ice}} = 0)$). Then, the incumbent point is updated if the iteration $k$ is successful, and these actions are performed until a convergence threshold $\epsilon_f$ is reached.

# V. Results and Discussion

In this section, the numerical results are presented. First, the benchmark optimal design obtained from the deterministic procedure is presented along with the results of the uncertainty propagation for the baseline and this benchmark configuration. Finally, the robust optimization procedure results for each uncertainty bound are presented for both problem formulations.

## A. Deterministic optimization

The deterministic design process of the AI-ETIPS targets $||\dot{M}_{\text{ice}}||_2$, which is the aggregated rate of ice accretion on the two sides of the airfoil, and it relies on a multi-start approach entailing ten optimization runs. Each run is initialized using a different incumbent point to promote the discovery of the global optimum in case multiple local minima are present. The convergence history of the 5 best run is reported in Fig. 4, together with a dashed line indicating the performance of the baseline configuration $||\dot{M}_{\text{ice}}||_2 = 6.73 \times 10^{-4} \text{kgs}^{-1}$. Generally, it takes just less than a few steps to discover a design outperforming the baseline design configuration. The picture shows that, within about 200 iterations, each run converges to a design with the same improved performance, requiring about 500 function evaluations in run 10. Different runs with different starting points converged to the same optimum design, assuring that the optimizer is capable of identifying a minimum. Nonetheless, different runs converged to different minima, pointing out the presence of several local minima. Table 3 presents the details of the designs proposed from the best performing runs. The results in terms of heat flux distribution are very similar among the best-performing runs, showing that for these runs, an optimal configuration has been found. Following this, the values of the QoI of these five runs are also very similar. In Table 3, the column "% Sides" measures the existing performance deviations between pressure and suction sides, and % improvements are calculated with respect to the baseline design. One can see that, given that the AoA is 0° and the water distribution is symmetric, the evaporative performance is better across the suction side than the pressure side. As mentioned in Section II, the IPS has a manufacturing shift; therefore, a larger heat portion is supplied to the suction side. However, the performance deviations between the sides are generally reduced due to the integration of the objectives employing the L2 norm. Among these best-performing runs, the most promising design is the one obtained from run #5. This will be used for comparison with the robust optimal designs.
<figure><img src="image_7.png" /><figcaption>Fig. 3 Flowchart of the robust optimization procedure</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 4 Convergence history of the best five runs of the NOMAD algorithm, depicting the best solution found so far by the algorithm.</figcaption></figure>

<table><caption>Table 3 Results of the best performing runs of the deterministic optimization</caption><thead><tr><th>Run</th><th>M&#x0307;<sub>ice,SS</sub> &times; 10<sup>-1</sup><br>[gs<sup>-1</sup>]</th><th>M&#x0307;<sub>ice,PS</sub> &times; 10<sup>-1</sup><br>[gs<sup>-1</sup>]</th><th>||M&#x0307;<sub>ice</sub>||<sub>2</sub> &times; 10<sup>-1</sup><br>[gs<sup>-1</sup>]</th><th>% Sides</th><th>% Improvements</th></tr></thead><tbody><tr><td>Baseline</td><td>4.66</td><td>4.85</td><td>6.73</td><td>3.92</td><td>-</td></tr><tr><td>#5</td><td>3.94</td><td>3.91</td><td>5.547</td><td>0.78</td><td>17.58</td></tr><tr><td>#3</td><td>4.01</td><td>3.84</td><td>5.549</td><td>4.26</td><td>17.55</td></tr><tr><td>#2</td><td>3.95</td><td>3.90</td><td>5.550</td><td>1.14</td><td>17.54</td></tr><tr><td>#7</td><td>3.89</td><td>3.95</td><td>5.550</td><td>1.51</td><td>17.53</td></tr><tr><td>#10</td><td>3.95</td><td>4.00</td><td>5.552</td><td>1.07</td><td>17.50</td></tr></tbody></table>
The results of the optimized deterministic configuration will be compared in detail later in this section with the baseline, robust optimized A and B configurations.

The baseline and the deterministic optimal design are now characterized with respect to the input uncertainty, namely the surrogates of the baseline and optimal solutions are built. The surrogate considers fixed (baseline and optimal) heat flux configuration but random cloud properties as inputs. The Monte Carlo UQ analysis is performed by exploiting the surrogate model to forward propagate the input uncertainty and obtain a stochastic characterization of the QoI targeted by the optimizer. Recalling that the surrogate model is trained with 100 samples for data set B1 and 50 samples for data set B2. Fig. 5(a) and (b) report the frequency density of ||$\dot{M}_{ice}$||₂ resulting for the pressure side of the baseline configuration and the deterministic optimal configuration, respectively, considering the B1 and B2 uncertainty ranges for SAT. The pictures also report an enlargement of the histogram close to the axes' origin. Fig. 6 shows the counterpart obtained for the optimized configuration. From the comparison of Fig. 5 and Fig. 6 it is possible to identify some

<figure><img src="image_9.png" /><figcaption>Fig. 5 Histogram of $\dot{M}_{ice}$ under the uncertain cloud environment, zoomed around $\dot{M}_{ice} > 0$ for the baseline configuration on the pressure side of the airfoil, for both uncertain temperature bounds B1 and B2. The distribution of $\dot{M}_{ice} > 0$ has been fitted with the beta distribution of maximum likelihood.</figcaption></figure>

common features. First, the observation of the picture on a larger scale highlights the encoding of a boolean-valued outcome, namely, the ice/no-ice formation. Indeed, many realizations are associated with a null $||\dot{M}_{ice}||_2$, giving birth to a largely populated bin at $||\dot{M}_{ice}||_2 = 0$. In all four cases, the probability of no ice formation ranges between 0.4 and 0.5 for a bin size collecting in between 25k-30k realizations. Moving rightwards to values $||\dot{M}_{ice}||_2 > 0$, all the zoomed plots reveal a frequency distributions resembling a beta probability distribution. From the full frequency distribution it is then possible to retrieve the quantities necessary to the robust optimization namely, the 95 quantile, $\dot{M}_{ice,j}^{95}$, and the probability of no ice formation, $P(\dot{M}_{ice} = 0)$.

As mentioned, the UQ analysis of the considered application is very demanding, so taking advantage of a surrogate of the computational model is required. The surrogate is comparatively less computationally expensive to evaluate,
<figure><img src="image_10.png" /><figcaption>Fig. 6 Histogram of $\dot{M}_{ice}$ under the uncertain cloud environment, zoomed around $\dot{M}_{ice} > 0$ for the optimized configuration on the pressure side of the airfoil, for both uncertain temperature bounds B1 and B2. The distribution of $\dot{M}_{ice} > 0$ has been fitted with the maximum likelihood beta distribution.</figcaption></figure>

yet it requires several realizations to be adequately trained to reproduce the full CFD model accurately. In order to determine the best trade-off between computational time and accuracy, several sizes of the Monte Carlo analysis (from $N_{MC} = 1 \times 10^3$ to $N_{MC} = 1 \times 10^6$) have been considered. The baseline and the deterministic optimal configurations results are shown in Fig. 7, for temperature uncertainty bounds B1 and B2. Plots report the evaluation of the quantile and no-ice frequency estimation against the number of Monte Carlo samples.

The estimated quantile and probability of no-ice formation varies by varying the number of Monte Carlo samples. Nevertheless, for both bounds B1 and B2 the difference in the results between $N_{MC} = 1 \times 10^5$, $N_{MC} = 5 \times 10^5$ and $N_{MC} = 1 \times 10^6$ is very small. Therefore, for the optimization procedure, it was decided to use $N_{MC} = 1 \times 10^5$ samples to reduce the computational cost without losing accuracy.

## B. Robust optimization data set B1

This section presents the robust optimization targeting either the quantile or the probability of no-ice formation, considering the uncertain bound set B1. Results are compared against the baseline and the deterministic optimal configurations in terms of performance.

Figure 8 shows the convergence plots for both robust optimization problems. The x axis represents the number of function evaluations of the optimization algorithm, the y axis is the 95 quantile $||\dot{M}_{ice}^{95}||_2$ for the robust optimization problem A (RA) while is the probability of no ice formation for the robust optimization problem B (RB). For clarity, only the best five runs are shown. All runs converge within about 200 - 300 evaluations. The optimization procedure is capable to identify a minimum. The multi-start approach returns multiple occurrence of the same optimal design, assuring independence from the initialization point. In any case, the proposed designs outperform the baseline configuration in terms of targeted QoI. On the other hand, the deterministic optimal design is outperformed only by some designs
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 7 Estimated 95 quantile of the QoI and estimated probability of $\dot{M}_{ice} = 0$, for B1 and B2 uncertainty bounds for the pressure side of the airfoil.**</figcaption></figure>

achieved under the robust optimization framework, e.g., not all the robust designs are associated with a performance quantile lower than the deterministically optimized IPS. It is important to notice that all runs converge after a limited number of evaluations, considering, for example, similar optimizations performed with genetic algorithm would have required about $O(10^5)$ function evaluations, as in [7]. This dramatically reduces the computational time.

Table 4 compares the quantities of interest for the baseline, optimized deterministic, and optimized robust configuration. Table 5 and Fig. 10 report the design parameters, i.e. the heat flux distribution for both problem formulations presented in Section III. They are labeled as A and B, respectively. It is shown the probability of no ice formation $P(\dot{M}_{ice} = 0)$ and $||\dot{M}_{ice}^{95}||_2$ considering the complete airfoil. Besides, performance for suction and pressure sides is presented to evaluate possible asymmetries. As stated before, the objective of the robust optimization A was to minimize the 95 quantile of the mass rate of ice accretion, as it can be seen in Table 4, the goal was reached. Indeed, with the optimized robust configuration the QoI $||\dot{M}_{ice}^{95}||_2$ has been reduced by 7.72% with respect to the baseline configuration and by 0.2% with respect to the reference deterministic optimal configuration. Although, for the second robust optimization (B), this was not the primary objective. There is a reduction of $||\dot{M}_{ice}^{95}||_2$ by 5.51% relative to the baseline configuration. While compared to the deterministic optimized configuration $||\dot{M}_{ice}^{95}||_2$ is increased by 2.18%.

By having a look at the $\dot{M}_{ice,j}^{95}$ for the pressure and suction sides, it can be noticed that all configurations present more or less symmetric results between the two side of the airfoil. It means that on one side of the airfoil, the occurrence of severe ice formation is as likely as on the other side. It must be noticed that the layout of the heater is not perfectly
<figure><img src="image_13.png" /><figcaption>Fig. 8 Convergence plot for the MADS algorithm for the robust optimization A (top) and B (bottom) for uncertainty bounds B1</figcaption></figure>

<table><caption>Table 4 95 quantile of $\dot{M}_{ice}$ and probability of no ice formation for the baseline, optimized deterministic, optimized robust A and B configurations for uncertainty bounds B1.</caption><thead><tr><th></th><th>P(Μ̇<sub>ice</sub> = 0)<br>[%]</th><th>||Ṁ<sup>95</sup><sub>ice</sub>||<sub>2</sub><br>[g s<sup>-1</sup>]</th><th>P(Μ̇<sub>ice</sub> = 0), PS Ṁ<sup>95</sup><sub>ice, PS</sub><br>[%]</th><th>PS Ṁ<sup>95</sup><sub>ice, PS</sub><br>[g s<sup>-1</sup>]</th><th>P(Μ̇<sub>ice</sub> = 0), SS Ṁ<sup>95</sup><sub>ice, SS</sub><br>[%]</th><th>SS Ṁ<sup>95</sup><sub>ice, SS</sub><br>[g s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Base.</td><td>48.62</td><td>1.703</td><td>48.72</td><td>1.194</td><td>49.62</td><td>1.214</td></tr><tr><td>Opt. Det.</td><td>44.62</td><td>1.575</td><td>53.00</td><td>1.101</td><td>45.49</td><td>1.127</td></tr><tr><td>Opt. Rob. A</td><td>44.27</td><td>1.571</td><td>53.77</td><td>1.122</td><td>44.62</td><td>1.100</td></tr><tr><td>Opt. Rob. B</td><td>53.22</td><td>1.609</td><td>53.47</td><td>1.139</td><td>53.22</td><td>1.137</td></tr></tbody></table>
symmetric, but there is a shift towards the CEG heaters' side. This would generally lead to not symmetric results if a proper heat flux distribution was not considered. On the other hand, this symmetry is even more noticeable for the optimized robust B configuration both the 95 quantile and the probability of no ice formation have similar values for pressure and suction side; the difference is in the range of 0.2% between the two sides. This is to be addressed to the fact that the objective was to maximize the probability of no ice formation on the overall airfoil. Therefore, to have a high probability of no ice formation, ice should not form on the suction and pressure side under the same environmental conditions.

Table 4 also reports the probability of no ice formation on the whole airfoil for all considered configurations. In this case, $P(\dot{M}_{\text{ice}} = 0)$ is increased by 9.45% for the optimized robust B configuration compared to the baseline. While it is reduced by 8.95% for the optimized robust A configuration with respect to the baseline. The evident symmetry of the optimized robust B configuration results can be noticed even more by looking at the probability of no ice formation on the suction and pressure side in Table 4. On the suction and pressure sides, there is basically the same probability of no ice formation. Again this clear symmetry is to be addressed to the objective function formulation, presented in Sec. III. The two sides of the airfoil are considered separately and then aggregated through the L2 norm. Among all possible configurations with the same overall rate of ice accretion, the symmetric ones are preferred by the optimizer, as expected. It must be noticed that a quantile reduction does not imply an increase of $P(\dot{M}_{\text{ice}} = 0)$. This means that the two objectives are, in the end, non-cooperating. Indeed, the opposite happens: the quantile is reduced because the tail of the distribution is reduced. This means less probability of very severe ice formations. However, there is a larger probability of mild ice formations on the airfoil when considering the uncertain variables. It is important to remind that, as explained in Section II.B, the modeling of the uncertain variables is very conservative, meaning that severe ice formations are overestimated. Therefore, even the design with a longer tail in the probability distribution will have a lower frequency of severe ice formations at the end.

The probability distribution of $\dot{M}_{\text{ice},j} \neq 0$ is shown for the four configurations in the violin plots in Fig. 9. It can be noticed that, for the optimized deterministic and optimized robust A configuration where the quantile is lower, the tail is slightly shorter, and the peak of the distribution is higher and shifted towards lower values of rate of ice accretion with respect to the baseline and the optimized robust B configuration that have higher quantile values. Lastly, Fig. 10 and Table 5 report the heat flux corresponding to each heater of the IPS for the baseline, optimized deterministic, optimized robust A and B configuration. First of all, it can be noticed that for all the optimized configurations, there is a minimal heat flux applied to the heater at the limits of the protected region. This means that the heat fluxes provided by previous heaters are usually sufficient to evaporate the water, saving the overall heat consumption. The partially dry heater may lead to more considerable convective losses, as the heat would be directly dissipated to air. Therefore lower heat fluxes are provided in this region. Moreover, the similarity between the optimized deterministic and the optimized robust A configuration can be seen. This similar design brings comparable results in uncertainty propagation, returning
<figure><img src="image_14.png" /><figcaption>Fig. 9 Violin plot representing the probability distribution of $\dot{M}_{ice,PS} \neq 0$ on the pressure and suction side of the airfoil for the baseline, optimized deterministic, optimized robust A and B configuration for uncertainty bounds B1.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 10 Heat flux of each heater for the baseline, optimized deterministic, optimized robust A and B configuration for uncertainty bound B1.</figcaption></figure>
equivalent statistics, as stated previously.

It is remarkable, the optimized robust B configuration has a quite different distribution with higher heat fluxes in the whole protected region but at the limits. This could be due to the objective of increasing the probability of no ice formation requiring a more even distribution of the heat fluxes. Finally, it can be seen an asymmetry on the heat flux distribution in all cases with lower heat fluxes in the suction side, this is due to the fact that there is a slight manufacturing shift of the heaters towards that side. Therefore, asymmetric heat fluxes must be provided to guarantee performances in terms of ice formation.

<table><caption>**Table 5 Profile of the heat fluxes across the heaters for the baseline, benchmark and optimized robust configurations for data set B1**</caption><thead><tr><th>Heat flux [Wm<sup>-2</sup>]</th><th>heater G</th><th>heater E</th><th>heater C</th><th>heater A</th><th>heater B</th><th>heater D</th><th>heater F</th></tr></thead><tbody><tr><th>Base</th><td>18 600</td><td>18 600</td><td>26 350</td><td>43 400</td><td>32 550</td><td>21 700</td><td>20 150</td></tr><tr><th>Opt. Det.</th><td>2 346</td><td>44 539</td><td>45 682</td><td>2 392</td><td>45 032</td><td>45 812</td><td>2 306</td></tr><tr><th>Opt. Rob. A</th><td>2 333</td><td>41 237</td><td>46 215</td><td>2 534</td><td>46 407</td><td>46 392</td><td>2 818</td></tr><tr><th>Opt. Rob. B</th><td>2 880</td><td>43 000</td><td>21 462</td><td>46 000</td><td>42 958</td><td>34 729</td><td>5 458</td></tr></tbody></table>

By analyzing these results, it can be seen that a reduction of the 95 quantile ||$\dot{M}_{ice}^{95}$||₂ does not imply also a reduction of the probability of no ice formation. Actually, there is an important decrease in the latter quantity, leading to a broader set of uncertain environmental conditions where ice formation occurs, although they are mild ice accretions. In order to proceed with a proper robust ice protection system design, it must be analyzed whether it would be more dangerous to have milder but more frequent ice formations or the opposite, keeping in mind that from the modeling of the uncertain cloud parameters overestimates severe ice formations. Considering that usually, even a small ice formation can lead to performance degradation for a proper robust design ice encounters should be avoided as much as possible. Therefore the robust optimization B procedure should be preferred to A to have a safer and more reliable ice protection system. Although the optimized design might not be the global optimum, it points towards several characteristics of well-performing designs that go beyond the intuitive design suggested in the literature.

### C. Robust optimization data set B2

In this section, robust optimization is performed considering the uncertain bound set B2. All the results presented are analogous for bound set B1. First, the convergence plots of the MADS algorithm are presented in Fig. 11. Also in this case, it takes about 200 function evaluations for the best five runs to converge to a locally optimal solution. Indeed, as in bound set B1, each run convergences to a different value, suggesting that the optimization algorithm remains trapped in a locally optimal region. Table 6 reports the performances associated with the optimal designs and the reference values associated with the baseline and the deterministic optimized system. First for the A configuration there
<figure><img src="image_16.png" /><figcaption>Fig. 11 Convergence plot for the MADS algorithm for the robust optimization A (top) and B (bottom) for uncertainty bounds B2</figcaption></figure>

<table><caption>Table 6 95 quantile of $\dot{M}_{ice}$ and probability of no ice formation for the baseline, optimized deterministic, optimized robust A and B configurations for uncertainty bounds B2</caption><thead><tr><th></th><th>P(Μ̇<sub>ice</sub> = 0)<br>[%]</th><th>||Ṁ<sup>95</sup><sub>ice</sub>||<sub>2</sub><br>[g s<sup>-1</sup>]</th><th>P(Μ̇<sub>ice</sub> = 0), PS Ṁ<sup>95</sup><sub>ice,PS</sub><br>[%]</th><th>PS Ṁ<sup>95</sup><sub>ice,PS</sub><br>[g s<sup>-1</sup>]</th><th>P(Μ̇<sub>ice</sub> = 0), SS Ṁ<sup>95</sup><sub>ice,SS</sub><br>[%]</th><th>SS Ṁ<sup>95</sup><sub>ice,SS</sub><br>[g s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Base.</td><td>48.58</td><td>1.772</td><td>48.84</td><td>1.243</td><td>49.21</td><td>1.263</td></tr><tr><td>Opt. Det.</td><td>46.97</td><td>1.642</td><td>49.96</td><td>1.145</td><td>48.28</td><td>1.177</td></tr><tr><td>Opt. Rob. A</td><td>1.18</td><td>1.602</td><td>1.68</td><td>1.235</td><td>2.09</td><td>1.021</td></tr><tr><td>Opt. Rob. B</td><td>50.70</td><td>1.682</td><td>50.76</td><td>1.189</td><td>50.77</td><td>1.190</td></tr></tbody></table>
is a reduction of ||$\dot{M}_{\text{ice}}^{95}$||₂ of 9.59% with respect to the baseline configuration, and of 2.44% compared to the optimized deterministic one, so once again the objective of the optimization is achieved. On the other hand, Table 6 the probability of no ice formation is very low, close to 1%. This means that there is ice formation for most cloud encounters, although the severity is low. This is clear from Fig. 12, where the probability distribution presents a large peak for a low $\dot{M}_{\text{ice},j}^{95}$ (both on the pressure and suction side). On the other hand, the B optimal configuration reduces the QoI quantile by only

<figure><img src="image_17.png" /><figcaption>Fig. 12 Violin plot representing the probability distribution of $\dot{M}_{\text{ice}} \neq 0$ for the baseline, optimized deterministic and optimized robust configuration for uncertainty bounds B2.</figcaption></figure>

5.08% w.r.t. the baseline, whereas it is increased by 2.43% w.r.t. the deterministic optimal configuration. Nonetheless, the probability of no ice formation increases compared to all reference configurations. This means that there is a smaller probability of ice forming on the overall airfoil leading to an overall higher security and reliability of the ice protection system.

Figure 13 reports the heat flux distribution along the airfoil for the different configurations. In all cases, low heat
<figure><img src="image_18.png" /><figcaption>Fig. 13 Heat flux of each heater for the baseline, optimized deterministic, optimized robust A and B configuration for uncertainty bound B2.</figcaption></figure>

flux is assigned to heater A, which is at the stagnation point. Indeed, at the stagnation point, the convective heat transfer is maximum; therefore, a large amount of heat would dissipate through convection. Since the optimization limits the convective dispersion, the optimal heat flux associated with the A heater is lower than the remaining ones. It is interesting to compare the heat flux distribution of the different configurations. In the robust designs, the heat fluxes are low except for the very end of the protected region, where a large amount of heat is provided to fully evaporate the water and avoid the formation of runback ice. Similarly, the deterministic optima design is characterized by a low heat flux associated with the A heater. A relatively uniform distribution of thermal power results for the C, E, B, and D heaters, whereas heaters F and G are kept to have a minimum heat flux.

For the uncertainty bound set B2, it is even clearer that a reduction of the 95 quantile $||\dot{M}_{\text{ice}}^{95}||_2$ does not automatically imply an increase of the probability of fully evaporative IPS operation. As stated in the previous section, an evaluation of the most dangerous situation should be performed. Once again, to avoid even mild ice formation that could lead to significant performance degradation, the second robust optimization procedure should be preferred to have a safer and more reliable ice protection system. Again the optimized design might not be the global optimum, but it shows characteristics of well-performing designs that go beyond the intuitive design suggested in the baseline design.

<table><caption>Table 7 Profile of the heat fluxes across the heaters for the baseline, optimized deterministic, optimized robust A and B configurations for data set B2</caption><thead><tr><th>Heat flux<br>[Wm<sup>-2</sup>]</th><th>heater G</th><th>heater E</th><th>heater C</th><th>heater A</th><th>heater B</th><th>heater D</th><th>heater F</th></tr></thead><tbody><tr><th>Base</th><td>18 600</td><td>18 600</td><td>26 350</td><td>43 400</td><td>32 550</td><td>21 700</td><td>20 150</td></tr><tr><th>Opt. Det.</th><td>2 346</td><td>44 539</td><td>45 682</td><td>2 392</td><td>45 032</td><td>45 812</td><td>2 306</td></tr><tr><th>Opt. Rob. A</th><td>46 000</td><td>2 300</td><td>2 300</td><td>2 300</td><td>3 300</td><td>45 000</td><td>44 069</td></tr><tr><th>Opt. Rob. B</th><td>2 316</td><td>45 715</td><td>38 377</td><td>13 042</td><td>41 815</td><td>43 030</td><td>4 996</td></tr></tbody></table>
## VI. Conclusions

This work deals with the optimization of evaporative Electro-Thermal Ice Protection Systems commonly employed in aeronautical applications to protect aircraft against the risk of in-flight ice formation. The optimization framework is applied to a reference test case from a well-known experimental campaign in order to test the validity of the framework. The results presented show that there is room for a significant improvement in the efficiency of the protection system, particularly regarding the allocation of the onboard power among the available heaters.

Despite a simple formulation/implementation of the optimization problem, fair improvements have been achieved with respect to the baseline configuration. The improvements concern quantities of interest, which are representative of the in-flight accretion risk, e.g., the amount of liquid water flowing on the aircraft's surface. The optimization problem is also formulated by accounting for the uncertainty which unavoidably affects reality. In particular, the uncertain properties of clouds, which onboard devices can not measure, and the uncertain SAT, which is measured and therefore known only up to a certain degree of accuracy. In this case, a robust design is sought in the optimization process, targeting the minimization/maximization of some statistics of the QoI. Though improvements are achieved, it was found that the diverse statistics do not necessarily vary as desired when only one of them is optimized. For instance, the minimization of the 95 quantile led to a reduction in the frequency of severe ice formation. Nevertheless, the probability of ice formations, although associated with milder events, was increased. For this reason, it is concluded that minimizing the frequency of no ice formations while constraining the quantile value is more adequate to design a robust AI-ETIPS that can guarantee safer flight conditions. Also, due to the conservative modeling of the cloud uncertainties, the magnitudes of the most severe ice cases are magnified.

The computational cost of the entire framework is rather demanding, and the optimization would require high computational power in practical applications. For this reason, we exploit a surrogate model that reduces the cost of obtaining the statistics of the targeted QoI. Nonetheless, additional work is required to thoroughly evaluate surrogate models other than the Gaussian Processes employed here. Additionally, further investigations should be conducted to evaluate the performance of the optimization algorithm with respect to the parameter set. It is noted that the algorithm converged to local minima in numerous runs, and therefore, it cannot be ensured if the optimal values found are the global minima. The presented framework is highly flexible. Different uncertain input parameters, objective functions, geometry, and environmental parameters can be set. This means that, even if the framework is tested on a reference common case, it could be easily adapted to a real-world application and address different research questions. In particular, this framework could be of high interest when considering a closed-loop IPS where a feedback signal like surface temperature is used to determine the energy requirements. In this case, the uncertainty would also be linked to the measurement tools used in the control loop. Finally, it can be stated that a preliminary methodology for the robust optimization of AI-ETIPS has been developed that can help in the design process of a safe, robust, and more efficient ice protection system.
# Funding Resources

The work in this paper was supported in part by the H2020-MSCA-ITN-2016 UTOPIAE, grant agreement 722734

## References

[1] Whalen, E., Broeren, A., Bragg, M., and Lee, S., "Characteristics of Runback Ice Accretions on Airfoils and their Aerodynamics Effects," *43rd AIAA Aerospace Sciences Meeting and Exhibit*, 2005, p. 1065. https://doi.org/10.2514/6.2005-1065.

[2] Regulations, F. A., "Airworthiness Standard: Transport Category Airplanes, Part 25 Appendix C," *US Code of Federal Regulations*, Title 14, FAA, 2022, pp. 536-542.

[3] Al-Khalil, K. M., Keith Jr, T. G., and De Witt, K. J., "Development of an improved model for runback water on aircraft surfaces," *Journal of aircraft*, Vol. 31, No. 2, 1994, pp. 271-278. https://doi.org/10.2514/3.46484.

[4] Morency, F., Tezok, F., and Paraschivoiu, I., "Anti-Icing System Simulation using CANICE," *Journal of Aircraft*, Vol. 36, No. 6, 1999, pp. 999-1006. https://doi.org/10.2514/2.2541.

[5] Pourbagian, M., and Habashi, W. G., "Surrogate-based Optimization of Electrothermal Wing Anti-icing Systems," *Journal of Aircraft*, Vol. 50, No. 5, 2013, pp. 1555-1563. https://doi.org/10.2514/1.C032072.

[6] Pourbagian, M., Talgorn, B., Habashi, W. G., Kokkolaras, M., and Le Digabel, S., "Constrained Problem Formulations for Power Optimization of Aircraft Electro-Thermal Anti-icing Systems," *Optimization and Engineering*, Vol. 16, No. 4, 2015, pp. 663-693.

[7] Arizmendi Gutiérrez, B., Della Noce, A., Gallia, M., and Guardone, A., "Optimization of a Thermal Ice Protection System by Means of a Genetic Algorithm," *International Conference on Bioinspired Methods and Their Applications*, Springer, 2020, pp. 189-200. https://doi.org/10.1007/978-3-030-63710-1_15.

[8] Saeed, F., and Paraschivoiu, I., "Optimization of a Hot-Air Anti-Icing System," *41st Aerospace Sciences Meeting and Exhibit*, 2003, p. 733. https://doi.org/10.2514/6.2003-733.

[9] Pellissier, M., Habashi, W., and Pueyo, A., "Optimization via FENSAP-ICE of Aircraft Hot-Air Anti-Icing Systems," *Journal of Aircraft*, Vol. 48, No. 1, 2011, pp. 265-276. https://doi.org/10.2514/1.C031095.

[10] Zhang, F., Huang, Z., Yao, H., Zhai, W., and Gao, T., "Icing Severity Forecast Algorithm Under both Subjective and Objective Parameters Uncertainties," *Atmospheric Environment*, Vol. 128, 2016, pp. 263-267. https://doi.org/10.1016/j.atmosenv.2015.12.066.

[11] Feng, K., Lu, Z., and Yun, W., "Aircraft Icing Severity Analysis Considering Three Uncertainty Types," *AIAA Journal*, Vol. 57, No. 4, 2019, pp. 1514-1522. https://doi.org/10.2514/1.J057529.

[12] Zhang, X., Lu, Z., Feng, K., and Ling, C., "Reliability Sensitivity based on Profust model: an application to Aircraft Icing analysis," *AIAA Journal*, Vol. 57, No. 12, 2019, pp. 5390-5402. https://doi.org/10.2514/1.J058475.
[13] DeGennaro, A. M., Rowley, C. W., and Martinelli, L., "Uncertainty Quantification for Airfoil Icing using Polynomial Chaos Expansions," *Journal of Aircraft*, Vol. 52, No. 5, 2015, pp. 1404-1411. https://doi.org/10.2514/1.C032698.

[14] Gori, G., Congedo, P. M., Le Maître, O., Bellosta, T., and Guardone, A., "Modeling In-Flight Ice Accretion Under Uncertain Conditions," *Journal of Aircraft*, Vol. 59, No. 3, 2022, pp. 799-813. https://doi.org/10.2514/1.C036545.

[15] Roy, R. K., *A primer on the Taguchi Method*, Society of Manufacturing Engineers, 2010.

[16] McAllister, C. D., and Simpson, T. W., "Multidisciplinary Robust Design Optimization of an Internal Combustion Engine," *Journal of Mechanical Design*, Vol. 125, No. 1, 2003, pp. 124-130. https://doi.org/10.1115/1.1543978.

[17] Marrison, C. I., and Stengel, R. F., "Design of Robust Control Systems for a Hypersonic Aircraft," *Journal of Guidance, Control, and Dynamics*, Vol. 21, No. 1, 1998, pp. 58-63. https://doi.org/10.2514/2.4197.

[18] Quagliarella, D., and Iuliano, E., "Robust Design of a Supersonic Natural Laminar Flow Wing-Body," *IEEE Computational Intelligence Magazine*, Vol. 12, No. 4, 2017, pp. 14-27. https://doi.org/10.1109/MCI.2017.2742718.

[19] Ghisu, T., Jarrett, J. P., and Parks, G. T., "Robust Design Optimization of Airfoils with respect to Ice Accretion," *Journal of Aircraft*, Vol. 48, No. 1, 2011, pp. 287-304. https://doi.org/10.2514/1.C031100.

[20] Al-Khalil, K., Horvath, C., Miller, D., and Wright, W., "Validation of NASA Thermal Ice Protection Computer Codes. Part 3; The Validation of ANTICE," *35th Aerospace Sciences Meeting and Exhibit*, 1997, p. 51. https://doi.org/10.2514/6.1997-5.

[21] Bernstein, B. C., McDonough, F., Politovich, M. K., Brown, B. G., Ratvasky, T. P., Miller, D. R., Wolff, C. A., and Cunning, G., "Current Icing Potential: Algorithm Description and Comparison with Aircraft Observations," *Journal of Applied Meteorology*, Vol. 44, No. 7, 2005, pp. 969-986. https://doi.org/10.1175/JAM2246.1.

[22] Cober, S. G., Isaac, G. A., and Strapp, J., "Aircraft Icing Measurements in East Coast Winter Storms," *Journal of Applied Meteorology*, Vol. 34, No. 1, 1995, pp. 88-100. https://doi.org/10.1175/1520-0450-34.1.88.

[23] Politovich, M. K., and Bernstein, T. A., "Aircraft Icing Conditions in Northeast Colorado," *Journal of Applied Meteorology*, Vol. 41, No. 2, 2002, pp. 118-132. https://doi.org/10.1175/1520-0450(2002)041%3C0118:AICINC%3E2.0.CO;2.

[24] Gori, G., Zocca, M., Garabelli, M., Guardone, A., and Quaranta, G., "PoliMIce: A Simulation Framework for Three-Dimensional Ice Accretion," *Applied Mathematics and Computation*, Vol. 267, 2015, pp. 96-107. https://doi.org/10.1016/j.amc.2015.05.081, the Fourth European Seminar on Computing (ESCO 2014).

[25] Gori, G., Parma, G., Zocca, M., and Guardone, A., "Local Solution to the Unsteady Stefan Problem for In-Flight Ice Accretion Modeling," *Journal of Aircraft*, Vol. 55, No. 1, 2018, pp. 251-262. https://doi.org/10.2514/1.C034412.

[26] Arizmendi Gutiérrez, B., Noce, A. D., Gallia, M., Bellosta, T., and Guardone, A., "Numerical Simulation of a Thermal Ice Protection System including State-of-the-Art Liquid Film Model," *Journal of Computational and Applied Mathematics*, Vol. 391, 2021, p. 113454. https://doi.org/10.1016/j.cam.2021.113454.
[27] Morelli, M., Bellosta, T., Donizetti, A., and Guardone, A., "Assessment of the PoliMIce toolkit from the 1st AIAA Ice Prediction Workshop," *AIAA Aviation 2022 Forum*, 2022, p. 3307. https://doi.org/10.2514/6.2022-3307.

[28] Economon, T. D., Palacios, F., Copeland, S. R., Lukaczyk, T. W., and Alonso, J. J., "SU2: An Open-Source Suite for Multiphysics Simulation and Design," *AIAA Journal*, Vol. 54, No. 3, 2016, pp. 828–846. https://doi.org/10.2514/1.J053813.

[29] Dussin, D., Fossati, M., Guardone, A., and Vigevano, L., "Hybrid Grid Generation for Two-Dimensional High-Reynolds Flows," *Computers & fluids*, Vol. 38, No. 10, 2009, pp. 1863–1875. https://doi.org/10.1016/j.compfluid.2009.04.007.

[30] Bellosta, T., Parma, G., and Guardone, A., "A robust 3D particle tracking solver for in-flight ice accretion using arbitrary precision arithmetic," *8th International Conference on Computational Methods for Coupled Problems in Science and Engineering, COUPLED PROBLEMS 2019*, CIMNE, 2021, pp. 622–633.

[31] Reid, T., Baruzzi, G. S., and Habashi, W. G., "FENSAP-ICE: unsteady conjugate heat transfer simulation of electrothermal de-icing," *Journal of Aircraft*, Vol. 49, No. 4, 2012, pp. 1101–1109. https://doi.org/10.2514/1.C031607.

[32] Da Silva, G. A. L., de Mattos Silvares, O., de Jesus Zerbini, E. J. G., Hefazi, H., Chen, H.-H., and Kaups, K., "Differential Boundary-layer Analysis and Runback Water Flow Model Applied to Flow Around airfoils with Thermal Anti-ice," *1st AIAA Atmospheric and Space Environments Conference*, 2009, p. 3967. https://doi.org/10.2514/6.2009-3967.

[33] Bu, X., Lin, G., Yu, J., Yang, S., and Song, X., "Numerical simulation of an airfoil electrothermal anti-icing system," *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, Vol. 227, No. 10, (2013), pp. 1608–1622. https://doi.org/10.1177/0954410012463525.

[34] Sand, W. R., Cooper, W. A., Politovich, M. K., and Veal, D. L., "Icing Conditions encountered by a Research Aircraft," *Journal of Applied Meteorology and Climatology*, Vol. 23, No. 10, 1984, pp. 1427–1440. https://doi.org/10.1175/0733-3021-23.10.1427.

[35] Audet, C., and Dennis Jr, J. E., "Mesh Adaptive Direct Search Algorithms for Constrained Optimization," *SIAM Journal on Optimization*, Vol. 17, No. 1, 2006, pp. 188–217. https://doi.org/10.1137/040603371.

[36] Audet, C., Le Digabel, S., and Tribes, C., "NOMAD User Guide," Tech. Rep. G-2009-37, Les cahiers du GERAD, 2009.

[37] Rios, L. M., and Sahinidis, N. V., "Derivative-Free Optimization: a Review of Algorithms and Comparison of Software Implementations," *Journal of Global Optimization*, Vol. 56, No. 3, 2013, pp. 1247–1293.

[38] Abramson, M. A., Audet, C., Dennis Jr, J. E., and Digabel, S. L., "OrthoMADS: A Deterministic MADS instance with Orthogonal Directions," *SIAM Journal on Optimization*, Vol. 20, No. 2, 2009, pp. 948–966. https://doi.org/10.1137/080716980.

[39] GPy, "GPy: A Gaussian Process framework in Python," http://github.com/SheffieldML/GPy, since 2012.