Identiﬁcation of airfoil polars from uncertain experimental measurements

. A new method is described to identify the aerodynamic characteristics of blade airfoils directly from operational data of the turbine. Improving on a previously published approach, the present method is based on a new maximum likelihood formulation that includes both errors in the outputs and in the inputs, generalizing the classical error-in-the-outputs only formulation. Since many parameters are necessary to meaningfully represent the behavior of airfoil polars as functions of angle of attack and Reynolds number, the approach uses a singular value decomposition to solve for a reduced set of observable 5 parameters. The new approach is demonstrated by identifying high quality polars for small-scale wind turbines used in wind tunnel experiments for wake and wind farm control research. The new proposed formulation was applied to the estimation of the aerodynamic characteristics of the blades of small-scale wind turbine models. This is a particularly difﬁcult problem, because an extended set of parameters is necessary in order to give a meaningful description of the polars, taking into account their variability with blade span, angle of attack and Reynolds number; invariably, this results in an ill-deﬁned problem because of the many unknown parameters and their possible collinearity. In addition, measurement errors affect both the outputs and the inputs, the latter being particularly relevant and representing the operating conditions of the turbines. On the other hand, accurate estimates of the polars are of crucial importance for the accuracy of simulation models based on lifting lines. Results indicate that a higher quality of the estimates is achieved by the proposed method, compared to an error-in-the-outputs only approach. Indeed, the estimated polars were for the ﬁrst time able to correctly model also derated operating conditions, which were not included in the parameter estimation process. All prior attempts at modeling these conditions failed to a various extent, when using polars estimated by the standard maximum likelihood formulation. In addition, results indicate that the present approach was successfully able to cope with the ill-posedness of the problem caused by the low observability 390 of the many unknown parameters, which is an aspect of importance for the practical applicability of the method to complex problems as the one considered here.

control (Campagnolo et al., 2016Frederik et al., 2019). In addition to the collection of valuable data sets in the known, repeatable, and controllable environment of the wind tunnel, the development and validation of digital copies of these experiments has been one of the main ambitions of this research effort. Both aeroelastic BEM  and LES-ALM (Wang et al., 2019) models of the experiments have been developed, in the latter case including not only the wind turbines, but also the wind tunnel and the passive generation of a sheared and turbulent flow. Results collected to date demonstrate an 45 excellent ability of the simulation models in reproducing the experiments, including multiple wake interactions and conditions relevant to wind farm control (Wang et al., 2019;Wang, Muñoz-Simón , 2020;Wang, Sharma, et al., 2020;. One crucial component of the simulation chain has been a method for estimating the polars directly from operational data of the turbines . In fact, the blades of scaled wind turbine models operate in low Reynolds regimes, where even relatively small changes in the operating conditions can cause significant changes in the aerodynamic characteristics of the 50 blade sections. In addition, given the small size of these models, even modest manufacturing imperfections and normal wear of the blades can lead to deviations from their nominal shape. Using the method of , the nominal airfoil polars are augmented with parametric correction terms, which are identified using a maximum likelihood (ML) criterion based on operational power and thrust measurements. These data points are collected on the turbine at various operating conditions, selected in order to span a desired range of angles of attack and Reynolds numbers. Since a large number of free parameters 55 are necessary to represent the correction terms, the resulting problem is ill-posed and the parameters are collinear. To address this issue, the original parameters are transformed into an new orthogonal set by using the singular value decomposition (SVD). Because the new parameters are uncorrelated with each other, one can select an observability threshold, discard the wind speed, the rotor angular velocity and the blade pitch setting. Errors in such quantities have a non-negligible effect on the outputs, and should be taken into account in a rigorous statistical sense.
To address this issue, the present paper proposes a new general formulation of ML identification that includes errors both in the outputs and in the inputs. This generalized formulation leads to an optimization problem in the model parameters and the unknown model inputs, which can now differ from their measured values. The proposed method is again cast within the SVD-70 based reformulation of the unknowns, to deal with the ill-posedness and redundancy of the parameters. The new formulation is applied to the identification of the polars of small-scale controlled wind turbines, developed to support wind farm control and wake research (Wang et al., 2019;Campagnolo et al., 2020;Frederik et al., 2019;Bottasso et al., 2020). Results indicate that the new formulation delivers polars of superior quality with respect to the original error-in-the-outputs only formulation.
Specifically, the new polars were able for the first time to correctly predict the turbine power outputs in derated conditions, 75 which had always defied previous efforts.
The paper is organized according to the following plan. Section 2 describes first the classical ML approach in §2.1 and its reformulation in terms on uncorrelated parameters in §2.2; §2.3 presents the novel ML method with errors in both outputs and inputs, while §2.4 discusses a way to take into account a priori information on the errors. Section 3 specializes the general formulation of §2.3 to the identification of the polars of scaled wind turbines. Finally, Sect. 4 presents the results, and conclusions 80 are drawn in Sect. 5.

Formulation
2.1 Classical maximum likelihood estimation with errors in the outputs

85
where u ∈ R l are the inputs (or, in the present context, the operating conditions), p ∈ R n the model parameters and y ∈ R m the outputs. In correspondence to the N inputs U = {u * 1 , u * 2 , . . . , u * N }, N experimental measurements of the outputs are available and noted Y = {y * 1 , y * 2 , . . . , y * N }. Because of modeling and measurement errors, the experimental measurements are in general not identical to the outputs predicted by model (1), a difference that can be quantified by the residual r = y * − y. The goal of the estimation problem is to find the model parameters p that minimize the residuals r.
A classical approach to this parameter estimation problem is the ML method (Klein, 2006). The idea of maximum likelihood estimation is to find the parameters p that maximize the probability J of obtaining the measurement sample Y, where J is written as R being the residual covariance and w i a weight assigned to the i-th residual. In this work, weights are introduced to account 95 for the fact that not all operating conditions appearing in the sample U might have the same importance. For example, it might happen that some u i 's represent frequent typical operating conditions of the system, whereas others are less frequent or relevant conditions. It might then be desirable to better match these more frequent conditions than the less frequent ones. One way to achieve this behavior from the ML estimator is to assign weights to the residuals. The weights could be proportional to the relative frequency of each operating condition in the lifetime of the system, or be inversely proportional to the distance of that 100 operating condition to some nominal behavior, a concrete example of this latter case being explained later in the results section.
A robust implementation of this optimization problem is obtained by the following iteration (Klein, 2006): 1. Assuming temporarily frozen parameters equal to p, minimize J with respect to R, which yields the following expression for the covariance matrix (Jategaonkar, 2015) 2. Assuming a temporarily frozen error covariance R, solve the minimization problem 3. Return to step 1, and repeat until convergence.
In the following, alternating between steps 1 and 2 is termed a "major" iteration. The internal iterations necessary for the 110 solution of the optimization problem at step 1 are termed in the following "minor" iterations.

Maximum likelihood estimation in terms of uncorrelated parameters
The estimation problem expressed by Eqs. (3,4) can be ill-posed, because of low observability and collinearity of the unknowns.
This is a classical difficulty in parameter estimation: on the one hand one would typically prefer a rich set of parameters that give ample freedom to adjust the behavior of a model in order to accurately match the measurements; on the other hand,

115
it might be difficult -if not altogether impossible-to always guarantee that there is enough informational content in the measurements to correctly identify and distinguish the effects of each one of the unknown parameters.
Indeed, the well posedness of the identification problem is associated with the curvature of the likelihood function with respect to changes in the parameters. Around a flat maximum, different values of the parameters yield similar values of the likelihood. A measure of the curvature of the solution space is provided by the Fisher information matrix (Jategaonkar, 2015).

120
The inverse of this matrix is also useful because it bounds the variance of the estimates (Cramér-Rao bound) (Jategaonkar, 2015). Unfortunately, the Fisher information by itself does not offer a constructive way of reformulating a given ill-posed problem.
To overcome this difficulty,  proposed to transform the original physical parameters of the model into an orthogonal parameter space. This mapping is obtained by diagonalizing the Fisher matrix using the SVD. As the new set 125 of variables are now statistically independent, one can readily select and retain in the analysis only the parameters that are associated with a sufficiently high level of confidence. Once the problem is solved, the uncorrelated parameters are mapped back onto the original physical space.
This approach enables one to solve an identification problem with many free parameters, some of which might be interdependent or not observable in a given data set. Furthermore, the SVD diagonalization reduces the problem size, retaining only 130 the orthogonal parameters that are indeed observable. Finally, this approach reveals, through the singular vectors generated by the SVD, the inter-dependencies that may exist among some parameters of the model, which may provide useful insight into the problem itself.
A detailed description of the SVD-based version of ML identification is given in . The same formulation is used also in the present paper.

Maximum likelihood estimation with errors in the inputs and outputs
The standard formulation of the ML identification presented in §2.1 considers the presence of noise in the outputs y. Indeed, outputs are affected by measurement errors but also, being computed through a model, by the deficiencies of the model itself.
Although errors in the outputs are typically the primary source of uncertainty in a parameter estimation problem, there are situations where significant errors may also be associated with the inputs u, which is the case of the present application. A 140 formulation of ML that accounts for errors both in the outputs and inputs is presented next.
The parametric model (1) is expanded as Because of modeling and measurement errors, the experimental output measurements y * are in general not identical to the model-predicted outputs y. Similarly, because of measurement errors and an imperfect realization of the operating conditions, 145 the experimental inputs u * are in general not identical to the nominal ones u. These differences can be synthetically quantified by the residual r = y * − y, where now y * is an expanded vector that contains measurements of both outputs and inputs: The goal of the estimation problem is to find the model parameters p and system inputs u i that maximize the probability of obtaining the measurements y * and u * . According to the maximum likelihood criterion, Eq. (4) becomes and Eq. (3) is now Instead of solving the problem in a monolithic fashion, the following iteration can be conveniently used: 3. Assuming temporarily frozen inputs u i , solve This is formally identically to the classical (error-in-the-output only) ML formulation, which can be solved by the SVDbased re-formulation in terms of uncorrelated parameters . 4. Assuming temporarily frozen parameters p, solve These are N decoupled small size problems, which return the values of the model inputs.
5. Return to step 2, and repeat until convergence.
This way the solution of the identification problem with input and output errors is obtained by using the classical error-in-the-165 output only ML implementation (using Eq. (9)), followed by a sequence of inexpensive optimizations to compute the model inputs (using Eq. (10)). Notice that, as long as it converges, this iteration returns the same result as the monolithic solution of problem (7,8).

Filtering of measurements based on a priori uncertainties
Often, a priori information on the expected uncertainties may be available. In such cases, the unknown true inputs u i can be 170 bounded as where ∆u are the expected uncertainty bounds. This a priori information can be used to retain in the cost function J only those measurements for which the corresponding residual cannot be simply explained by the uncertainties (11), but must be due to the model parameters p.

175
To this end, notice first that the residual r i is a function of p and u i , i.e.
Indicating the j-th component of residual r i as r ij , its maximum and minimum values for a given p are computed as If the maximum r M ij and minimum r m ij have different signs, then r ij = 0 lies somewhere within this range and hence this residual component can be fully explained by input uncertainties. Therefore, it cannot drive meaningful changes in the parameters, and should be neglected. Otherwise, this residual carries valuable information and should be retained. To account for this, a filtered residualr ij is defined as The a priori estimates are used to initialize the parameters p at step 1 of the iterative algorithm formulated in §2.3. A standard ML method is used for the initialization, considering only errors in the outputs and using Eqs. (3,4) where the residual components r ij are replaced by the filtered onesr ij . Filtering accelerates the optimization, because it avoids meaningless tuning of parameters caused by measurement noise. Once this initial estimate of the parameters is obtained, it is further refined 190 by considering the a posteriori effects of noise in inputs and outputs by stepping through points 2-5 of the algorithm. Residual filtering is not used further, because it is based on a priori assumptions relying on knowledge of the measurement chain, which can only estimate bounds and might not reflect the actual noise effectively experienced for any given measurement.
In practice, a naive implementation of filtering can be very expensive. In fact, as the residual r i depends on p, one would have to recompute the optimization problems (13) each time the parameters are updated, which becomes prohibitively expensive.

195
The cost of filtering can be drastically reduced with a simple approximation, as graphically illustrated in Fig. 1. The figure shows with a blue dotted line the residual component r ij as a function of the input u i for a given value of the model parameters p (0) . The counter (·) (0) refers to the values that the parameters assume at the beginning of each major iteration used to solve problem (4). The minimum and maximum of this curve, corresponding to r m ij and r M ij , are respectively indicated with blue lower and upper pointing triangles. These stationary points are computed at the beginning of each major iteration, by 200 solving problems (13). For simplicity, this is obtained by a simple evaluation of the residuals over a regular subdivision of the unknowns.
At the k-th minor iteration of the solution of problem (4), the model parameters have been updated and they now assume the value p (k) . The corresponding function r ij is depicted in the figure with a red solid line, together with its new minimum ij , i.e. the difference in the residual evaluated at the nominal inputs u * i for the two parameter values p (k) and p (0) . The shifted function is shown by the black dashed curve in Fig. 1. This is an inexpensive operation since it does not require any optimization. This nominal difference is then used for shifting the minimum and maximum residuals from their 210 initial value at p (0) to the new value at p (k) . By this approximation, the maximum and minimum residuals are readily and inexpensively updated at each iteration as Based on these updated values, the residual filtering condition expressed by Eq. (13) can be readily updated.

215
This approximation works very well in practice since the interval [u * i − ∆u, u * i + ∆u] is small. In addition, by a standard Taylor series analysis, one can show that this approximation entails neglecting terms that are quadratic in the changes of the parameters within a major iteration, which are typically small. Finally, the approximation does not affect the quality of the results, as the true stationary points are recomputed at each new major iteration of the ML algorithm. In this sense, the approximation only speeds up the calculations of the minor iterations, but the results -at convergence of the major and minor 220 loops-are the same that would have been obtained by a straightforward (but more expensive) solution of problem (13).
The parameter identification problem setting described in the previous pages is completely general, and could be used for a wide range of applications. However, for the specific problem at hand and with reference to Eq. (1), the outputs are defined as y = (C P , C T ) T , where C P = 2P/(ρAV 3 ) and C T = 2T /(ρAV 2 ) are respectively the rotor power and thrust coefficients, 225 and P is power, T thrust, ρ air density, A = πR 2 the rotor swept area, R the rotor radius and V the wind speed. The inputs describe the rotor operating conditions and are defined as u = (ρ, V, Ω, β) T , where Ω is the rotor angular velocity and β the blade collective pitch angle. To obtain the power and trust coefficients, nominal values of the inputs are used for both the measured and predicted cases.
The airfoil lift and drag coefficients, respectively noted C L and C D , are now assumed to be in error, and the goal of the 230 estimation problem it to calibrate them in order to match a given set of measurements. This is achieved by defining changes ∆C L and ∆C D with respect to nominal values C L0 and C D0 , i.e.
where η is the spanwise location along the blade (because different airfoils are typically used at different stations along a 235 rotor blade), α is the local angle of attack and Re = uc/ν the local Reynolds number, u being the relative flow speed, c the chord length and ν the kinematic viscosity of air. The dependency of these functions on spanwise location, angle of attack and Reynolds number is approximated using assumed shape functions and their associated nodal parameters p C L and p C D , which therefore represent the tunable algebraic parameters of the model, i.e.
Following , instead of working directly with p = (p C L ; p C D ), which might not be all observable, these variables are first transformed by the SVD into an uncorrelated set of parameters, which are then truncated with a variance threshold, calibrated according to the ML criterion and are finally projected back to the original functional space ∆C L and ∆C D .

245
The dependency of y on p and u is expressed through model (1) using Blade Element Momentum (BEM) theory (Manwell et al., 2009), as implemented in FAST (OpenFast, 2020).
The typical Reynolds number distribution along a wind turbine blade is almost constant for the majority of its span, but assumes smaller values close to the blade tip and root. The implementation of this paper, improving on the work of , specifically considers that the airfoil polars depend on Re. The expected range of Reynolds numbers is discretized by 250 linear shape functions and associated nodal values, and the local Reynolds number is computed at each spanwise station based on local geometry and flow conditions. The results presented later on consider scaled wind turbine models for wind tunnel testing. For these rotors, the chord-based Reynolds number is much lower than in typical full-scale applications, and ad hoc low-Reynolds airfoils (Lyon and Selig, 1998) are used. Because of the special flow regime of these airfoils, the formulation is complemented by the conditions ∂C L /∂Re > 0 and ∂C D /∂Re < 0. The first of these conditions accounts for the earlier 255 reattachment of the laminar separation bubble on the suction side of the airfoil for increasing Re, and the second for the shorter chord extent of that same bubble (Selig and McGranahan, 2004). They are enforced as soft penalty constraints in problem (4) by modifying the cost function as J = J + J p , where

Experimental setup
A scaled wind turbine model of the G1 type (Campagnolo et al., 2016) was operated in the boundary layer wind tunnel of the Politecnico di Milano in low turbulence (1%) conditions. The rotor blade design is based on one single low-Reynolds airfoil of the RG14 type (Lyon and Selig, 1998). Measurements of the rotor thrust and power were obtained for 158 different operational  Table 1 reports a priori estimates of the uncertainties associated with the various measured quantities. Given the uncertainties on the measurements, worst-case uncertainties on the power and thrust coefficients can be readily computed as The wind speed V was measured by a Mensor CPT-6100 pitot transducer (Mensor, 2016), which is affected by pressure and alignment errors. The pitot tube measures the dynamic pressure, i.e. the difference ∆p = 1/2 ρV 2 between the total and the 275 static pressures. Since the wind speed is computed by inverting the dynamic pressure expression, errors in ∆p and ρ directly pollute V . Additionally, a yaw and tilt misalignment may exist between the pitot axis and the incoming wind vector, increasing the error in V . The uncertainty of the air density was estimated from the hygrometer and barometer installed in the wind tunnel.
After considering all relevant factors, the uncertainty of the wind speed was determined using the guidelines described in ISO (2008). The uncertainty in the blade pitch angle β was estimated by calibrating the actuator encoder with a Wyler Clinotronic

280
Plus inclinometer (Campagnolo, 2013). Power was computed as P = QΩ, where Q is the torque, which was measured by strain gages at the rotor shaft. These sensors were calibrated by applying a known torque to the locked rotor. The rotor speed Ω was measured by an optical incremental encoder with a count per revolution N e = 10, 000 and an observation window t ow = 4 ms, which results in an error ∆Ω = 1/N e t ow ≈ 1.5 rpm. The thrust T was obtained by measuring with a strain gage bridge the fore-aft bending moment at tower base; here again, the strain gages were calibrated by applying a known load to the turbine 285 by a pulley and weight system. The contribution to the bending moment due to the drag of nacelle and rotor was obtained by a dedicated experiment in the wind tunnel without the blades. Additional details on sensors and error quantification are discussed in Campagnolo (2013) and .
For each wind speed V , a turbine should operate at a specific TSR λ and blade pitch β, which are computed in region II to maximize power capture and in region III to limit power output to the rated value. On the other hand, for the task of identifying was assigned based on its distance to the nominal conditions, computed as where (·) * indicates a nominal value and 1/2/3 are scaling factors. All data points were divided into four groups according to their distance. Data points within each group were assigned the same weight, longer mean distances corresponding to lower weights.

300
Nominal values of the blade polars are defined as the ones previously computed with the method of .
Although of a good quality, these polars are however not always able to correctly represent the behavior of the turbine, for example in derated conditions. To improve on this situation, the method proposed here was used to further correct the polars and provide improved estimates.
The lift and drag coefficients were parameterized in terms of bi-linear shape functions, using seven nodal values for Reynolds 305 and 21 for angle of attack for each one of the two coefficients. Since the G1 blades use one single airfoil type along their entire span, it was not necessary to introduce the dependency on η appearing in the general expressions of Eqs. (16).
For the nominal polars, Fig. 2 plots the variance σ 2 (which is the inverse of the singular values produced by the SVD analysis) for the seven considered Reynolds numbers and the lowest 25 modes. The identification first used nominal model inputs u * and the residual filtering technique of §2.4 to identify an initial guess to the system parameters p, a process that converged after nine major iterations of Eqs. (3,4). For the converged solution, Fig. 3 315 shows the nominal model inputs (two upper plots) and the output residuals (two lower plots), including the nominal residual r * , the maximal residual r M , the minimal residual r m and the filtered residualr (see Eqs. (13,14)) for each one of the measured data points. The filtered residualsr are zero for most conditions, indicating that the information carried by these data points cannot be distinguished further from input measurement noise. In addition, all non-zero filtered residuals are small, indicating an almost singular R, which is in fact used as termination criterion.

320
An a priori estimate of the maximal uncertainties on the power and thrust coefficients can be computed based on Eqs. (19) and Table 1, which yields On the other hand, an a posteriori estimate of the uncertainties evaluated with nominal inputs u * is available by the covariance 325 matrix R of Eq. (3) that, using unfiltered residuals, gives As expected, the a posteriori estimates are smaller than the a priori ones, since the latter represent a worst case scenario.
The process was then continued, using the previously converged parameters as an initial guess, and now adding to the 330 identification also the model inputs u to include the effects of their uncertainties. After three iterations, a converged solution was obtained. The final identified inputs are denoted in the following as u I . For all operational conditions, Fig. 4 shows the differences ∆u = u I − u * between identified and nominal values. In all subplots, two horizontal dashed lines indicate the a priori uncertainties reported in Table 1. It is interesting to observe that most estimated inputs are within the a priori bounds, indicating a good coherence between a priori and a posteriori statistics. The right part of the same figure reports the distributions 335 of the errors that, except for wind speed, are close to normal. On the other hand, density appears to have a small bias, which violates one of the assumptions of ML estimation.    Table 2 reports the correlation coefficients, computed from the extended covariance matrixR at convergence as ij = R ij /(σ i σ j ), where σ k = R kk . Because of symmetry, only the upper triangle is shown.
The correlation coefficient between the two outputs, ∆C P and ∆C T , is negative. This means that, on average, at the end 345 of the identification process the power and thrust residuals have opposite signs. This is expected, since it allows to minimize the cost function of problem (7). Additionally, each input induces variations of the same sign in the two outputs; for example, a larger wind speed or density imply higher power and thrust coefficients, whereas a larger blade pitch implies lower power and thrust coefficients. Given that ∆C P and ∆C T have a negative correlation, the input-output correlation coefficients always have different signs for both outputs, e.g. (∆C P , ∆β) and (∆C T , ∆β) have opposite signs. The signs of the input-input 350 correlations can be explained in similar terms. For example, the correlation between density and blade pitch is negative because these two inputs have correlations of opposite sign to the outputs, whereas the correlation between blade pitch and wind speed is positive because these two inputs have correlations of the same sign to the outputs. From the extended covariance matrix at convergence, the mean absolute a posteriori uncertainties of the inputs |u I − u * | were found to be 0.06 m/s for speed V , 0.09 deg for blade pitch angle β, 0.5 rpm for rotor speed Ω, and 0.005 kg/m 3 for density 355 ρ. By comparison with Table 1, all a posteriori uncertainties are smaller than the a priori ones, as expected.

Power derating cases
To verify the quality of the identified polars, derated operational conditions were considered. It should be stressed that these conditions were not included in the identification data set, and therefore provide for a verification of the generality of the results. These additional conditions are listed in Table 3 and correspond to values equal to 100%, 97.5%, 95% and 92.5% of 360 rated power. Figure 6 shows the results in terms of power (on the left) and thrust (on the right) coefficients, as functions of derating percentage. In all plots, the experimental results are shown using a solid blue line with * symbols; whiskers indicate the uncertainties according to Eq. (19) and Table 1. Simulation results are computed with nominal measured inputs u * , both for the nominal polars p * according to  and for the newly identified polars p I , and they are marked with  .

Conclusions
This paper has presented a new maximum likelihood identification method that, departing from the classical formulation, accounts for errors both in the outputs and in the inputs. The new method is a generalization of the classical approach, where 370 the system parameters are estimated together with the system inputs, which this way can differ from their actually measured quantities because of noise. The new expanded formulation is solved using a partitioned approach, resulting in an iteration between the standard parameter estimation and a series of decoupled and inexpensive steps to compute the inputs. To cope with the ill-posedness of the problem caused by low observability of the parameters, the formulation uses an SVD-based transformation into a new set of uncorrelated unknowns, which, after truncation to discard unobservable modes, are mapped 375 back onto the original physical space. The formulation is further improved by an initialization step that accounts for a priori information on the errors affecting the measurements, discarding all data points whose residuals can be simply explained by uncertainties.
The new proposed formulation was applied to the estimation of the aerodynamic characteristics of the blades of smallscale wind turbine models. This is a particularly difficult problem, because an extended set of parameters is necessary in 380 order to give a meaningful description of the polars, taking into account their variability with blade span, angle of attack and Reynolds number; invariably, this results in an ill-defined problem because of the many unknown parameters and their possible collinearity. In addition, measurement errors affect both the outputs and the inputs, the latter being particularly relevant and representing the operating conditions of the turbines. On the other hand, accurate estimates of the polars are of crucial importance for the accuracy of simulation models based on lifting lines.

385
Results indicate that a higher quality of the estimates is achieved by the proposed method, compared to an error-in-theoutputs only approach. Indeed, the estimated polars were for the first time able to correctly model also derated operating conditions, which were not included in the parameter estimation process. All prior attempts at modeling these conditions failed to a various extent, when using polars estimated by the standard maximum likelihood formulation. In addition, results indicate that the present approach was successfully able to cope with the ill-posedness of the problem caused by the low observability 390 of the many unknown parameters, which is an aspect of importance for the practical applicability of the method to complex problems as the one considered here.
Code and data availability. An implementation of the polar identification method and the data used for the present analysis can be obtained by contacting the authors. Author contributions. CW developed the a priori residual filtering method, wrote the software, performed the simulations and analyzed 435 the results; FC was responsible for the wind tunnel experiments and the analysis of the measurements; CLB devised the original idea of estimating polars from operational turbine data, developed the ML formulation with errors in inputs and outputs and supervised the work; CW and CLB wrote the manuscript. All authors provided important input to this research work through discussions, feedback and by improving the manuscript.
Competing interests. The authors declare that there is no conflict of interest.