A novel monotonic wind turbine power-speed characteristics model

Abbreviations: DE: Differential evolution; EA: Evolutionary algorithm; EP: Evolutionary programming; GA: Genetic algorithm; GCR: Generator control region; GWEC: Global wind energy council; LF: Logistic function; LSE: Least-Square error; MAE: Mean absolute error; MAPE: Mean absolute percent error; NRMSE: Normalized root mean square error; PDF: Power distribution function; PSO: Particle swarm optimization; RMSE: Root mean square error; STIR: Subspace trust-region; TI: Turbulence intensity; WT: Wind turbine; WTPSC: Wind turbine power-speed characteristics; SDS: Significant downwards shifting; SUS: Significant upwards shifting; MDS: Minor downwards shifting; MUS: Minor upwards shifting; CBS: Curve becomes steeper; CBF: Curve becomes flatter; NDS: Negligible downwards shifting; NUS: Negligible upwards shifting; DbS: Decrease (in $q(\vec{\theta}, { } v)$) by horizontal shifting; IoD: Increase or decrease (in $q(\vec{\theta}, { } v)$) by horizontal shifting

1.   Introduction

WT installations experienced an unprecedented annual growth of $53 \%$ in $2020$. GWEC reports that, with 93 GW of new installations in $2020$, the global cumulative capacity reached 743 GW ^[1].

Such growth is attributed to factors such as: 1. innovations in blades design and fabrication ^[2]; 2. improvement in nacelle components reliability (e.g., gears with fatigue and wear resistance and new softwares for system level modeling of the nacelle) ^[2]; 3. improved sensors and control algorithms ^[2]; 4. political support (e.g., low carbon energy goals, feed-in tariffs, or guaranteed access to the transmission grid) ^[3,4].

WTPSCs play significant roles in: 1. risk assessment ^[5]; 2. wind energy yield and WT selection ^[5]; 3. condition monitoring ^[6]; 4. sizing storage capacity for wind power integration ^[7]; 5. predictive control optimization ^[8]; 6. WTs troubleshooting ^[8]; 7. detection of degradation because of aging ^[6]; 8. optimal dispatching of wind farms ^[6].

Near cut-in wind speed, WTPSCs are difficult to model ^[9]. Although we theoretically expect WTs power to exhibit a cubic relationship with wind (when it is below rated speed), this is incorrect. The conversion efficiency of WTs varies with speed, which is most remarkable near cut-in and rated wind speeds. Just slightly above cut-in speed, virtually all WTs exhibit a steep growth in efficiency. Near rated wind speeds, WTs 'spill' wind energy ^[10], decreasing the efficiency. The result is complex and nonlinear relationship between wind speed and the WT power output ^[11].

The general shape of LFs (sometimes referred to as logistic distributions or activation functions) made them meaningfully suitable for scientific modeling of bacteria and plant growth since a growing population 'competes' for resources, placing an upper limit on the number of bacteria/plants ^[12]. LFs can also be useful for certain biological, chemical, linguistic, political science topics ^[13].

Since both of LFs and WTPSCs are 'S'-shaped, LFs are (at least in principle) candidates for accurate WTPSC mathematical models. LFs are parametric models (they are based on mathematical expressions with a fixed number of parameters), bringing potential to offer some analytical advantages. For instance, in wind energy assessment, it is possible to use them along with Weibull distribution to obtain explicit PDFs ^[14]. Lydia et al. ^[15] are among the earliest researchers to attempt modeling WTPSCs curves using LFs. They propose a $5$-parameter LF and estimate its parameters using GA, EP, PSO and DE. Many researchers have sought new methods to estimate the parameters of the LF proposed by Lydia et al. ^[15] (e.g., ^[6,16,17]). Sohoni et al. ^[9] reported that certain LFs (due to explicitly including an inflection point in their parameters) have the potential to increase accuracy and improve online-monitoring using LFs. Pei and Li ^[11] confirmed this by parameterizing various WTPSC models and comparing them using statistical metrics, such as MAE and RMSE.

Jing et al. ^[18] improved the model of Lydia et al. ^[15] by introducing 'quantile dependency' to its parameters. Jing et al. ^[18] parameterized the model proposed by Lydia et al. ^[15] using PSO for three wind farms and validating the accuracy of their work using MAPE and NRMSE.

Villanueva and Feijóo ^[19] proposed two 'generalized' $6$-parameter LFs and estimate their parameters using evolutionary optimization techniques. The proposed models show substantial accuracy improvement, but Villanueva and Feijóo ^[19] reported difficulties in optimizing their parameters.

Zou et al. ^[17] reviewed three, four, five, and six parameter LFs and their accuracy using MAE and RMSE. Results show that six parameter LF may sometimes provide lower accuracy than four and three parameter LF. Thus, adding new parameters to LFs does not necessarily contribute to increased accuracy. Villanueva and Feijóo ^[19] reported similar results.

In Section 2, we provide a mathematical background on WTPSC and LFs and highlight their limitations and our contributions. Section 3 proposes a new LF model, demonstrates its curve fitting merits, and develops upper and lower limits on its parameters. We validate our model in Section 4 by providing graphical and numerical accuracy results. To ensure a suitable comparison, we compare our work with other 6-parameter LFs. Section 5 draws final remarks, discusses the limitations of the study and puts forth suggestions for future investigations. A summary of the literature gap and the most significant contributions are then mentioned in Section 6.

2.   Background and contribution summary

2.1.   Mathematical representation of WTPSC

For stall or pitch-controlled WT, it is possible to characterize the WTPSC curves mathematically as ^[5]

where $v$ is the wind speed. $v_{ci}$, $v_{r}$, and $v_{co}$ are the cut-in, rated, and cut-out speeds of the WT, respectively. $P_{r}$ is the rated output power of the WT. $q_{e}(\vec{\theta}, { } v)$ is a mathematical expression which should (ideally) accurately estimate the manufacturer-provided data (i.e., empirical power-speed pairs in the GCR). Such estimation is typically done by fitting a vector of parameters $\vec{\theta}$ using the least-squares method ^[5].

2.2.   Preliminary concepts

Figure 1 depicts a WTPSC curve. When estimating $\vec{\theta}$ parameters to fit the manufacturer datasheet, some points of interest assist us in establishing its limits.

1. The curve is monotonic. This means it always increases (or at least remains the same) with $v$. Mathematically, this means the first derivative must remain greater than or equal to zero, as Eq 2:

2. For virtually all WTPSC curves, an inflection point exists. This is the point at which the curve switches from concave upwards to concave downwards. Mathematically, this means the second derivative must be zero, as Eq 3:

Where $v_{\mathrm{inf}}$ is the inflection point.

3. One major attribute of LFs is that they should approach an upper asymptote when their dependent variable approaches infinity. Mathematically, this means

Ideally, the upper asymptote of WTPSC should be close to $P_{r}$. This trend of approaching an upper asymptote is obvious in virtually all WTPSC curves when $v \to v_{r}$. In practice, $v$ is never $\infty$, but this mathematical formulation benefits in establishing upper and lower limits for $\vec{\theta}$ parameters, which, in turn, assists the optimization algorithm by searching in a suitable region instead of arbitrarily determining it by trial and error.

2.3.   Limitation of previous methods and summary of contribution

$6$-parameter LFs are discussed in Villanueva and Feijóo ^[19] and have one of the following forms

where $a$ and $d$ are the upper and lower asymptote, respectively. $b$ is the growth rate (sometimes called 'hill slope') around point $c$. $g$ is called the 'asymmetry factor' because it tunes the degree of asymmetry around point $c$ ^[11,15]. $c$ is the inflection point (the point at which $q(\vec{\theta}, { } v)$ turns from concave upward to concave downward). $\varepsilon$ has no special meaning, but has a value around one ^[8].

Disadvantages of fits in Eqs 4 and 5:

1. With some WTs, they inaccurately fit the WTPSC, especially near $v_{ci}$ and $v_{r}$. This is reported by as a disadvantage for piecewise models, which is the conventional way to account for the effect of turning off the WT when $v < v_{ci}$ and $v > v_{co}$ ^[9]. Yan et al. ^[20] reports such inaccuracy for LFs even with continuous WTPSC;

2. Contain a high number of parameters compared to other fits with fewer parameters, but better accuracy. For example, Zou et al. ^[17] review three, four, five, and six parameter LF and tested their accuracy using MAE and RMSE. Results showed that $6$-parameter LF may sometimes provide lower accuracy than three, four and five parameter LFs. Pei and Li ^[11] report similar results. Thus, adding new parameters does not contribute to increased accuracy. Intuitively, adding new parameters should provide certain advantages, such as accuracy improvement;

3. Their monotonicity is not guaranteed when the parameters are estimated using EAs. We overcome this disadvantage in Aldaoudeyeh et al. ^[8];

4. Sometimes LFs have the possibility of becoming ill-conditioned, making the estimation of $\vec{\theta}$ parameters difficult ^[9].

This paper contain the following contributions:

1. We propose a new $6$-parameter LF model for WTPSC curves; 2. We develop constraints on the parameters of our model to guarantee its monotonicity; 3. The model exhibits excellent fit that is at least as accurate as 6PL and 6PLE for some WTs, but also provides significant MAPE and NRMSE improvements. It also provides high accuracy near $v_{ci}$ and $v_{r}$; 4. The parameters of our model are easily optimized with the STIR algorithm.

3.   The proposed model

Our proposed model is

where $a$ is the upper asymptote. $d$ is a parameter to tune the lower asymptote near $v_{ci}$ (i.e., $d$ is not the lower asymptote, but it drastically influences it). $b$ is the growth rate around point $c$. $g$ is called the 'asymmetry factor' because it tunes the degree of asymmetry around point $c$. $c$ is the inflection point. $\zeta$ is the steepness tuning factor.

We call this model 6PLEZ. At a glance, 6PLEZ model resembles 6PL and 6PLE (Eqs 4 and 5), but, as we will see in Section 4, it is more accurate. To our best knowledge, the literature does not contain this model. Villanueva and Feijóo ^[21] and Zou et al. ^[17] review WTPSC models (including LFs models), but the 6PLEZ model does not appear on any of them.

We estimate $\vec{\theta}$ parameters with the same method mentioned in Aldaoudeyeh et al. ^[8] (LSEs objective function minimized using STIR algorithm).

3.1.   Characteristics of the 6PLEZ model

In this subsection, we illustrate the effect of $\vec{\theta}$ parameters variation on $q(\vec{\theta}, { } v)$ shape. Unless otherwise is specified, the parameters are $a = 1$, $b = 1.25$, $c = 6$, $d = 0$, $g = 1$, and $\zeta = 0$. The dependent parameter ($v$) is shown from $2.5$ (typical $v_{ci}$) to $12$ (typical $v_{r}$).

Figure 2 shows the effect of varying $a$ and $d$ on the shape of $q(\vec{\theta}, { } v)$. Clearly, increasing $a$ slightly shifts the curve upwards near $v_{ci}$. This effect becomes more pronounced the closer we get to $v_{r}$, where the upper asymptote varies linearly with $a$. The final value of $q(\vec{\theta}, { } v)$ is $a$. We thus call it the upper asymptote. We also note that increasing $a$ steepens the curve near $c$.

The effect of $d$ is better illustrated with $a = 10$ and when we extend the curve from $0$ to $v_{r}$. Increasing $d$ shifts the upwards curve near $v_{ci}$ and its effect diminishes at $v_{r}$. We also note that increasing $d$ flattens the curve near $c$.

Thus, most of $a$ and $d$ effect is to scale the maximum and minimum asymptotes of the LF curve, respectively. By contrast, $b$, $c$, $g$, and $\zeta$ (as we will discuss now) define the shape of the LF.

Figure 3 shows the effects of $b$, $c$, $g$, and $\zeta$ variation on $q(\vec{\theta}, { } v)$ shape. The increase in $b$ shifts the curve downwards near $v_{ci}$ and upwards near $v_{r}$. However, it is also obvious from the curve that the effect is more pronounced near $v_{ci}$. As $b$ increases, the curve becomes steeper near the inflection point ($v_{\mathrm{inf}} = c$).

The variation of $c$ only shifts the curve to the left and right. In fact, it is worth mentioning that all curves for varying $c$ are the same, but merely shifted horizontally. It is possible to conclude from Figure 3 that increasing $c$ decreases $q(\vec{\theta}, { } v)$ near $v_{ci}$, and increases or decreases $q(\vec{\theta}, { } v)$ near $v_{r}$.

The increase in $g$ drastically downshifts the curve near $v_{ci}$. It also substantially downshifts the curve near $v_{\mathrm{inf}}$ and makes it steeper, but has no effect on the power curve near $c$. As we increase $\zeta$, $q(\vec{\theta}, { } v)$ shifts slightly upward and becomes flatter near $c$, but has no effect on $q(\vec{\theta}, { } v)$ near $v_{r}$. In a sense, $\zeta$ adjusts the steepness around $c$ while having a negligible effect minor increase effect $v_{ci}$ and a negligible decrease near $v_{r}$.

Note that $b$ and $\zeta$ seem to provide the same effect. Remarkable differences, however, are: 1. $b$ tunes the steepness while having significant effect near $v_{ci}$ and a minor effect near $v_{r}$; 2. $\zeta$ tunes the steepness while having minor effect near $v_{ci}$ and a negligible effect near $v_{r}$.

Table 1 shows a summary of the previously mentioned effects. In Figure 4, we list some examples on why our model provides higher accuracy. The dashed lines show curves with low accuracy in some region and the high accuracy in others, while solid lines are improvements due to variations in $\vec{\theta}$. Circles are manufacturer-provided data.

In Figure 4a, we see low accuracy near $v_{ci}$ and $v_{\mathrm{inf}}$. The error is positive and tends to decrease the closer we get to $v_{r}$, but turns negative at $v_{r}$. In this case, increasing $b$ substantially increases the accuracy since it shifts the curve downward near $v_{ci}$, moderately downward near $v_{\mathrm{inf}}$, and slightly shifts the curve upward near $v_{r}$. In Figure 4b, we see a case where substantial errors occurring near $v_{ci}$ but they disappear as the speed increases and accuracy is very high near $v_{r}$. A suitable solution is to increase $g$, which shifts the curve downward with larger shifting occur closer to $v_{ci}$, but no effect at $v_{r}$. Figure 4c shows a curve providing high accuracy near $v_{ci}$ and $v_{r}$. The error, however, to the left/right of $v_{\mathrm{inf}}$ is negative/positive, with a slope greater than the slope of manufacturer-provided data. In this case, increasing $\zeta$ improves the accuracy by flattening the curve near $v_{\mathrm{inf}}$.

The previously mentioned examples demonstrate the flexibility of the 6PLEZ model and why it provides significant accuracy improvements. Our model adapts to improving the accuracy in some regions (near $v_{ci}$, $v_{\mathrm{inf}}$ or $v_{r}$) without sacrificing the accuracy in the others.

3.2.   Constraints on the proposed six-parameter LF (6PLEZ)

In this section, we devise some limits on the 6PLEZ model parameters (Eq 6), which help the STIR algorithm minimize the objective function (i.e., LSEs). We define the main limits on $\vec{\theta}$ as (the subscripts $\min$ and $\max$ denote the minimum and maximum limits of each parameter, respectively)

3.2.1.   Monotonicity conditions

To guarantee increasing monotonicity, the derivative of $q(\vec{\theta}, { } v)$ must be positive or zero. Thus, differentiating Eq 6 we get

The inequality of Eq 8 is satisfied when

These monotonicity conditions mean that Eqs 9–13 must be true

3.2.2.   Asymptotes conditions

Taking the limit of Eq 6 as $v \to \infty$ yields

Equation 14 is the exact value of the upper asymptote of Eq 6. However, in practice, $v < \infty$, but since we still want to help the STIR algorithm estimate $a$ such that $q(\vec{\theta}, { } v)$ still fits the manufacturer-provided data well, we allow the value of $a$ to vary within $\pm 10\%$ of $P_{r}$, resulting in the conditions

Equation 6 must be greater than $0$ for all $v > 0$. Thus, $a - d$ must always be positive. This requires that

Note that the conditions in Eqs 15 and 17 embody the condition in Eq 9.

3.2.3.   Other limits

$\left[\frac{v}{c} \right]^{\zeta} { } {\rm{e}}^{-b \left(v - c \right)}$ in Eq 6 must start to substantially decrease at some point between $v_{ci}$ and $v_{r}$. In other words, as $v$ increases, $\left[\frac{v}{c} \right]^{\zeta} { } {\rm{e}}^{-b \left(v - c \right)}$ should become much smaller than one at some point between $v_{ci}$ and $v_{r}$. Otherwise, the curve would never approach the upper asymptote (which is almost equal to $P_{r}$) when $v \to v_{r}$. Thus, the inflection point $c$ must be somewhere between $v_{ci}$ and $v_{r}$. We formulate the limits on $c$ as follows

Note that the condition in Eq 18 embodies the condition in Eq 11. $v_{\mathrm{inf}}$ (i.e., $c$) of an LF is the point at which a curve switches from concave upward to concave downward (i.e., the second derivative must be zero). Thus, taking the second derivative of Eq 6 with respect to $v$ and equating it with zero

Substituting $v = c$ in Eq 20 and simplifying

then solving Eq 21 for $g$ as Eq 22

by taking the limits for various extreme values of $\zeta$, we obtain Eq 23

this suggests that

We solve Eq 22 for $\zeta$

Equations 25a and 25b mean for every $b { } c$ (the product of the two parameters) and $g$ values, there are (in general) two $\zeta$ values to satisfy the condition $\left. \frac{{{d^2}{q_e}(\vec \theta ,v)}}{{d{v^2}}}\right|_{v = v_{\mathrm{inf}}} = 0$ (Eq 3). By plotting both equations in 3D for $g \sim \left[0, { } 0.95 \right]$ and $b { } c \sim [0, { } 100]$, we see that the maximum value of $\zeta$ is $180$. Thus, we say

Equations 13 and 26 dictate the upper limits for $\zeta$ to satisfy the monotonicity and the inflection point existence conditions, respectively. Since we must satisfy both conditions, it follows that

where '$\mathrm{minof}$' means the minimum of the two choices.

3.3.   Summary of the 6PLEZ LFs constraints

Equation 10, Eqs 12, 15–19 and 24 and 27 define nine constraints necessary for the STIR algorithm to quickly and reliably parameterize $\vec{\theta}$ in the 6PLEZ model (Eq 6). The three limits in which we did not derive analytically (or with reasoning) are: (1) $b_{\max}$ and (2) $d_{\min}$, (3) $\zeta_{\min}$.

We determined such limits by extensive simulations with the WT database of various ratings as follows: (1) $b_{\max} = 3$ and (2) $d_{\min} = -0.25 { } P_{r}$, (3) $\zeta_{\min} = -500$. Table 2 shows a summary of the limits for the 6PLEZ logistic model.

3.4.   Objective function definition

The main function to be minimized is the LSE between the fits of Eqs 4–6 and empirical power-speed data provided by the manufacturer. This is given as Eq 28:

where $q_{e}\big(\vec{\theta}_{k}, \thinspace v_{i}\big)$ are curves of Eqs 4–6; $v_{i}$, $q_{m}\left(v_{i}\right)$ are the $i^{th}$ wind speed and wind power output data point as provided by the manufacturer, respectively; $\vec{\theta}_{k}$ is the parameters of vector $\theta$ in Eqs 4–6 at the $k^{th}$ iteration; and $n$ is the number of data points in the GCR (region between $v_{ci}$ and $v_{r}$ in 1) as provided by the manufacturer. Aldaoudeyeh et al. ^[8] discussed further details on the application of STIR algorithm for WTPSC curve fitting.

Compared to other optimization algorithms, advantages of the STIR algorithm include: 1. an efficiency in solving large bound-constrained minimization problems ^[22,23]; 2. An appeal for solving non-linear problems ^[22]; 3. An ability to handle convex, nonconvex, and ill-conditioned objective functions with large number of variables ^[24].

4.   Results and discussion

We fit $\vec{\theta}$ parameters using STIR algorithm as described by Aldaoudeyeh et al. ^[8]. MAE, MAPE, RMSE and NRMSE as defined in ^[6,25] are the statistical metrics we use to demonstrate the accuracy of our proposed model.

4.1.   Accuracy near cut-in and rated wind speeds

To ensure as a fair comparison as possible, we compare our work to the closest one. Villanueva and Feijóo ^[19] provide fits for 6PL model (Eq 4) for multiple WTs. We choose four of them for comparison purposes: (1) Enercon E82 E2; (2) Repower MM82; (3) Siemens S82 SWT-2.3 82; (4) Vestas V164/8000.

Figure 5 shows curve fits for the 6PL model (as estimated by Villanueva and Feijóo ^[19]) and the 6PLEZ model (as estimated in this work). Although both models contain the same number of parameters, our model (despite being a parametric one) fits the WTPSC curve near $v_{ci}$ and $v_{r}$. Such inaccuracy is reported for LFs (see Yan et al. ^[20], for example), but it does not occur in our model.

4.2.   Benchmarking with WTs database

We test our model with $42$ on-shore and off-shore WTs from $29$ manufacturers. Ratings of the WTs range from 275 to 8000 kW. Tables 3 and 4 list the estimated $\vec{\theta}$ parameters, while Figures 6 and 7 show statistical metrics graphically. For convenience of presentation, we sorted WTs from with descending order of accuracy improvement. Clearly, the 6PLEZ model is either more accurate or at least as accurate as the 6PLE model.

Table 6 shows four statistical metrics (MAE, RMSE, MAPE, and NRMSE) with their mean, median, maximum and minimum values for the entire database of WTs. Villanueva and Feijóo ^[21] provided classifications of the accuracy of WTPSC models depending on MAPE values. Such categories are summarized in Table 5. Tables 5 and 6 show that: (1) Both mean and median MAPE values for the 6PLE and 6PLEZ models fall in the 'very high' category. However, the 6PLEZ model has a median MAPE of 0.002383, which is less than half the 'very high' category threshold and (2) both mean and median NRMSE for the 6PLEZ model are 32.3% and 38.5% less than for the 6PLE model. This is a substantial accuracy improvement while maintaining the same number of parameters.

5.   Final remarks, limitations and recommendations

One limitation of the study is the use of manufacturer-provided data. Because of several factors (e.g., wind power curtailment, accumulation of dirt and snow, sensor failures, and pitch angle control malfunctioning ^[26,27,28]), the real power-speed data points may not be exactly the same as the ones provided by the manufacturer. In short, the data was characterized by noise and outliers. According to the IEC 61400-12-1 standard, outliers are highly weighted by the ordinary LSE regression method ^[29]. This conclusion should also pertain to most WTPSC fitting methods, as virtually all of them employ the LSE as an objective function. Pei and Li ^[11] confirm that outliers compromise the accuracy machine-learning WTPSC fitting methods.

When fitting using the manufacturer datasheet, the data points were obtained under controlled conditions. For example, the air mass density was set to $1.225\ {\rm{kg}}\ {\rm{m}}^{-3}$, the wake effect was neglected (because only one WT was tested at a time), and the TI was assumed to be $10 \%$ ^[30]. The conditions in the actual WT location are usually not the same. The implementation of test condition control promoted fairness in comparing different WTPSC fitting models or estimation methods.

When fitting with real power-speed data, the air mass and TI may fluctuate for similar wind speed measurements, while the wake effect was greatly affected by the distance between WTs and their yaw angles. Additionally, several data cleaning techniques were devised in the literature (see ^[6,20,31]), which could be incorporated in future research. The amount of data we must remove, however, changes depends on the WT itself, making the entire data cleaning process susceptible to subjectivity ^[32], and sometimes outliers may still be present in the processed data ^[28].

In summary, although many papers addressed the effect of outliers on the accuracy of WTPSC models and the sensitivity of their parameters (e.g., ^[32,33]), the literature contained only few for LFs WTPSC models and most of them were in the last five years (e.g., ^[11,18,20]). Bilendo et al. ^[31] provide a recent review of WTPSC, but with no conclusive remarks on the effect of outliers on the performance of LFs WTPSC models. Thus, we suggest the following as future work:

1. Development of new robust LFs WTPSC curve fitting techniques (or objective functions) to guarantee model resiliency against outliers;

2. Incorporating additional parameters into the model to account for typical factors affecting outliers;

3. Examining the correlation between data filtering techniques and criteria and the accuracy of LFs WTPSC modeling;

4. Formulating probabilistic LFs WTPSC curves to accommodate fluctuations in instantaneous turbulence in wind speeds or enhanced LFs WTPSC to guarantee that the 10-minute averaged power closely approximates the empirical data.

6.   Conclusions

In this paper, we proposed the 6PLEZ model (Eq 6). We then derived upper and lower limits for its parameters (2) and estimated them using the STIR algorithm. The significance of this paper is:

1. Our model adapts to improving accuracy in some WTPSC regions without sacrificing the accuracy in others (3.1);

2. Unlike most parametric approaches, the 6PLEZ model provides accurate WTPSC modeling near $v_{ci}$ and $v_{r}$ (4.1);

3. The limits on the parameters of the 6PLEZ model are clearly specified, which:

(a) guarantee the monotonicity of the model (Eq 8);

(b) provide robust and reliable estimation of the parameters with no difficulties in 'guessing' their range (Table 2);

4. It provides substantial accuracy improvement over the 6PLE model despite having the same number of parameters as the 6PL and 6PLE models (Tables 5 and 6).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.