A Revision of Empirical Models of Stirling Engine Performance Using Simple Artiﬁcial Neural Networks

: Stirling engines are currently of interest due to their adaptability to a wide range of energy sources. Since simple tools are needed to guide the sizing of prototypes in preliminary studies, this paper proposes two groups of simple models to estimate the maximum power in Stirling engines with a kinematic drive mechanism. The models are based on regression or ANN techniques, using data from 34 engines over a wide range of operating conditions. To facilitate the generalisation and interpretation of results, all models are expressed by dimensionless variables. The ﬁrst group models use three input variables and 23 data points for correlation construction or training purposes, while another 66 data points are used for testing. Models in the second group use eight inputs and 18 data points for correlation construction or training, while another 36 data points are used for testing. The three-input models provide estimations of the maximum brake power with an acceptable accuracy for feasibility studies. Using eight-input models, the predictions of the maximum indicated power are very accurate, while those of the maximum brake power are less accurate, but acceptable for the preliminary design stage. In general, the best results are achieved with ANN models, although they only employ one hidden layer.


Introduction
The oldest heat engines in use today are Stirling cycle engines.It is often emphasised that this type of cycle has the highest achievable efficiency under ideal conditions, but its greatest advantage might be the possibility of using a wide variety of energy sources, with low levels of chemical and noise pollution.Apart from refrigeration and cryogenic applications, where the Stirling cycle has achieved remarkable developments [1], technological and commercial vicissitudes over more than two hundred years have restricted the current applications of the Stirling engine practically, to cases where there is less competition.These range from anaerobic underwater equipment [2] to the combined production of heat, cooling, and electricity, especially in systems powered by waste heat [3] and renewable energies on a low power scale [4][5][6][7].Additionally, variants of the Stirling engine are among the alternatives being evaluated for power generation in space, both for space missions and for terrestrial use [8,9].
Comparisons with other technologies are based on technical indicators, such as efficiency, mean effective pressure and specific power, as well as economic indicators.Regarding the latter, units based on Stirling engines present, in general, higher investment costs than units based on internal combustion engines, since they are not currently massproduced, but they are economically competitive with fuel cell-based units.Their durability and relatively low maintenance costs [7] are additional strengths to the other advantages already mentioned.
The so-called Beale number is often used to estimate the specific power or mean effective pressure of a Stirling engine, not only in preliminary technical feasibility studies but also in academic analyses [10][11][12].William Beale observed that the maximum brake power in well-developed engines is roughly proportional to the mean cycle pressure, cycle volume amplitude, and engine rotational frequency.Walker called Beale's number, N B , the proportionality constant [13].Thus, the product of N B and the mean pressure would be approximately equal to the mean effective pressure which, multiplied by the rotational frequency, would be equal to the specific power of the engine.
The N B ≈ 0.15 value was obtained from the data of 24 engines, most of them operating with a heater temperature of about 650 • C. Nonetheless, using data from 14 engines with kinematic drive mechanisms, including the original 1843 Stirling engine, and 8 free piston Stirling engines (FPSE), West [14] deduced that the assumed proportionality constant varied in the range of 0.06 to 0.23, and could be expressed as a function of the ratio of the absolute temperatures of the heater and cooler.Reader and Hooper proposed that this function was linear [15].Other authors, using data from six operating points measured on five engines, two of them in the low temperature difference range, derived a correlation between the dimensionless maximum brake power and the main operating variables, i.e., pressure and temperatures, by a dimensional analysis [16].The proposal has been interpreted as a variant of West's correlation, with doubtful applicability to other engines due to the small number of parameters and influential variables considered in the model [17].However, it should be recognised that these authors pointed out the need to obtain an additional correlation to estimate the engine speed corresponding to the maximum power.
Previously, Organ [18] and, independently, Prieto et al. [19,20], had already deduced that the dimensionless maximum indicated power is a function of more than twenty dimensionless parameters and variables, laying the foundations for the similarity of Stirling engines and the application of scaling techniques for engines with a kinematic drive mechanism [21,22].Subsequently, Formosa and Fréchette [23] studied the similarity of FPSEs.Prieto and co-workers [24] deduced, from basic thermodynamic concepts and experimental data, that the engine speed corresponding to the maximum indicated power is not an independent variable in engines with a kinematic drive mechanism, but a function of the engine parameters and operating variables.They also analysed the leakage and mechanical losses of indicated power, and extended the model previously established for the gas circuit to the analysis of mechanical efficiency and brake power [25].
Correlations based on a few variables often lead to predictions of insufficient accuracy, even within the scope of their main objective, which is to guide the preliminary sizing of new prototypes.Therefore, correlations have been obtained based on a larger number of influential variables, using 21 operating points of eight engines, both for the dimensionless values of maximum indicated and brake power, as well as for the corresponding dimensionless speeds [26].Although the use of dimensionless variables facilitates the generalisation of the results, and the study has recently been extended to 54 operating points of 10 engines [27], the lack of data certainly limits the applicability of the correlations for engines not included in the databases.
Recently, soft computing methods based on genetic algorithms, particle swarm optimisation, fuzzy logic, and artificial neural networks (ANN), are an alternative to regressionbased correlations.Soft computing methods have, so far, mainly been applied to design or optimise particular Stirling engines, often FPSEs [28].A smaller number of soft-computing models are known to have been developed using the databases of a set of engines.The temperature of the heat source, rotation speed, fuel, and pressure have been selected as input variables in ANN-based approaches proposed to predict the torque and brake power of a Stirling engine, using experimental data, obtained mostly from the Philips M102C engine, to test, validate, and train the model [29,30].The temperature of the heat source, rotation speed, and type of engine configuration have given good results as input variables for other torque and brake power prediction models, which have used a larger amount of data from several engines for training, testing, and validation, but limited to the ranges of 453-1273 K, 46-1800 rpm, and 0.3-500 W [31].
The main novelties of this paper are the following: • Two groups of models are proposed to replace the Beale, West, and similar correlations for feasibility studies and the preliminary sizing of Stirling engines with a kinematic drive mechanism, as complementary tools to more advanced procedures for analysis, design, and optimisation.The aim is to cover the widest possible diversity of size, power, and ranges of operating variables, even though this may penalise the accuracy of the model predictions.It is rare to find experimental data on Stirling engines in the literature that are extensive enough to be used reliably for model building.In addition, the stages of design, manufacture, testing, optimisation, and commercial development of a prototype are a relatively long evolutionary process, so that not many engines have progressed beyond research and demonstration phases.It is, therefore, recognised that some objectives of the article are subject to the quality of the databases used and future improvements; • Group G1 includes models that, in combination, allow the maximum brake power to be estimated using three input variables.The G1-1R and G1-2R models are based on regression fitting techniques, while the G1-3ANN and G1-4ANN models are based on ANN methodologies.In G1 models, only the working gas, heater and cooler wall temperatures, mean pressure, and cycle volume amplitude are used as input variables.
The results are evaluated using data collected from the literature for 89 operating points of 34 engines with very different characteristics, which provides a wide overview of the ranges of parameters and operating variables; • Group G2 includes four correlations for estimating the maximum indicated and brake power, obtained by regression in a previous study [27], which are referred to in this paper as G2-5R, G2-6R, G2-7R, and G2-8R models.These correlations use eight input variables, often not available in the literature, and were, therefore, evaluated on data from a smaller number of engines.The results are compared in the article with those obtained using the new G2-9ANN, G2-10ANN, G2-11ANN, and G2-12ANN models, based on neural networks.
The use of dimensionless input and output variables allows the influence of a larger number of parameters and variables to be considered, and facilitates the generalisation of results.Comparisons between models are made in all cases using appropriate statistical indicators to highlight the advantages and limitations of each type of model.

Selection of Input Variables
The indicated power P ind developed by a Stirling engine with a kinematic drive mechanism depends on the following set of variables:

•
Heater wall temperature T wE and cooler wall temperature T wC ; • Mean pressure of the working gas, p m ; • Phase angle α and other geometrical parameters of the pistons and the drive mechanism, l 1 , . . ., l n , which determine the shape and volume amplitude V 0 of the thermodynamic cycle; Applying Buckingham's theorem, with p m , V 0 , T wC and n s as reference variables, leads to the following equivalent functional relationship between dimensionless variables [20]: where , and N TCR = ρ R c R T wC /p m .The characteristic Mach number N MA is a dimensionless velocity form representing the influence of the operating frequency n s .In order to make this influence explicit, the following equation has been proposed [24]: where ζ 0 is the dimensionless quasi-static work per cycle, i.e., a theoretical limit of the gas circuit performance that, irrespective of the configuration, operating point and working fluid, cannot be reached by real engines, while the coefficients Φ and Ψ are macroscopic representations of the indicated power losses associated with irreversibilities inherent to working gas friction and heat transfer.Depending on whether the drive mechanism is harmonic or non-harmonic, ζ 0 can be calculated with analytical equations or by numerical simulation using, in both cases, an isothermal model of the thermodynamic cycle.Functionally, ζ 0 depends on the temperature ratio, the parameters of the drive mechanism, and the dead volumes, but is independent of the working fluid and the mean pressure, so the following relationship can be written: The following conclusions can be drawn from Equation (2) [24]:

•
The dimensionless maximum indicated power, ζ ind,max , must be within the following range of variation: • The dimensionless speed corresponding to the maximum indicated power value, N MA,max , is not an independent design parameter, but a function which depends on the same parameters influencing the gas circuit performance;

•
The coefficients of indicated power losses are inversely proportional to N MA,max , as the following equations show: On the other hand, the mechanical losses of indicated power caused by leakage or friction depend on the same variables influencing the operation of the gas circuit, and on additional parameters characteristic of the seals and the drive mechanism, the influence of which can be expressed by mechanical efficiency.Therefore, the maximum brake power P B,max depends on more than thirty parameters and operating variables.

Input Variables for Group G1 Models
A simplified version of Stirling engine operation is to assume that the maximum brake power P B,max can be estimated by knowing only the main operating variables, i.e., where n * s,max is the engine speed corresponding to the peak power point.From this expression, Buckingham's theorem leads to the following equivalent functional relationship: where ζ B,max = P B,max / p m V 0 n * s,max is the dimensionless maximum brake power.Particular cases of this approach are the Beale, Reader and Hooper, and West correlations, which can be expressed, respectively, by the following equations: For the G1 models in this article, it is assumed that P B,max also depends on the properties of the working fluid and on the main variables of the thermodynamic cycle, i.e., from which Buckingham's theorem leads to the following functional relationship: On the other hand, since it is not possible to calculate P B,max from ζ B,max without knowing the corresponding speed, the following functional relationship is assumed analogously: which leads to: The Rayleigh method allows Equations ( 13) and (15) to be expressed as follows for the G1-1R and G1-2R models, respectively: For the G1-3ANN and G1-4ANN models, the output variables ζ B,max and N * MA,max will be estimated from the same three input variables, i.e., γ, τ and N p .

Input Variables for Group G2 Models
For the models of group G2, the maximum values of indicated and brake power, and their corresponding speeds, are estimated by the following functional relationships [27]: Note that a large number of influencing variables are implicit in the dimensionless variables, while others have been considered to be specific to a later stage of dimensioning, e.g., hydraulic radii and wetted areas of the heater and cooler, wetted area, porosity and material properties of the regenerator, as well as characteristics parameters of the seal rings and drive mechanism.
The data analysis procedure in models G2-5R, G2-6R, G2-7R, and G2-8R is based on regression fitting using functional relationships of the following type: In models G2-9ANN, G2-10ANN, G2-11ANN, and G2-12ANN, the eight independent variables in Equation (19) are also used as input variables to estimate the same four output variables ζ ind,max , N MA,max , ζ B,max and N * MA,max .

Characteristics of Neural Networks
A neuron is a single cell living in a network of cells that receives inputs, processes those inputs, and generates an output [32].Neurons can be classified into three groups or layers, as shown in Figure 1.The input layer receives the numerical data of the variables of interest x 1 , . . ., x n used to estimate the target values a 1 , . . ., a k of the network, and provides input values to the neurons that form the hidden layers of the structure, whose inputs and outputs are not accessible from outside the ANN.A neural network can have as many hidden layers as necessary, increasing the complexity of its structure, as only neurons from adjacent layers can be connected.
For the models of group G2, the maximum values of indicated and brake power, and their corresponding speeds, are estimated by the following functional relationships [27]: Note that a large number of influencing variables are implicit in the dimensionless variables, while others have been considered to be specific to a later stage of dimensioning, e.g., hydraulic radii and wetted areas of the heater and cooler, wetted area, porosity and material properties of the regenerator, as well as characteristics parameters of the seal rings and drive mechanism.

Characteristics of Neural Networks
A neuron is a single cell living in a network of cells that receives inputs, processes those inputs, and generates an output [32].Neurons can be classified into three groups or layers, as shown in Figure 1.The input layer receives the numerical data of the variables of interest  , … ,  used to estimate the target values  , … ,  of the network, and provides input values to the neurons that form the hidden layers of the structure, whose inputs and outputs are not accessible from outside the ANN.A neural network can have as many hidden layers as necessary, increasing the complexity of its structure, as only neurons from adjacent layers can be connected.The output of the neurons is calculated using Equation (20).The parameter  depicts the output value of the j-th neuron, which is calculated as the outcome of a certain activation function , given a linear combination of the output values of the previous layer.Thus,  is the input of the n-th neuron of the layer,  is the weight or gain factor of the input signals to the neuron, and  is the independent term or offset value of each neuron.According to the literature, the activation function could be, among other options, an arctangent, the absolute value, or a ramp function [33].The reason of this element is to introduce a non-linear component in the behaviour of the system.

𝑦
(20) The output of the neurons is calculated using Equation (20).The parameter y j depicts the output value of the j-th neuron, which is calculated as the outcome of a certain activation function f , given a linear combination of the output values of the previous layer.Thus, x n is the input of the n-th neuron of the layer, w nj is the weight or gain factor of the input signals to the neuron, and k j is the independent term or offset value of each neuron.According to the literature, the activation function could be, among other options, an arctangent, the absolute value, or a ramp function [33].The reason of this element is to introduce a non-linear component in the behaviour of the system.
Finally, the nodes of the output layer receive the outputs of the outermost hidden layer, and provide the final computation of the entire ANN.In this layer, the output of each node is calculated in the same way as in any of the previous layers.Given the structure of an ANN, even with a limited number of neurons and hidden layers, the number of variables involved in the network is large enough to make it difficult to find a direct solution to the best configuration of the network.Therefore, an optimisation or training process is required to reach the final configuration.
To carry out the training process, it is necessary to define a certain cost function and use learning algorithms, whose objective is to obtain for each neuron the weights and biases defined in Equation (20), that optimise the value of the cost function.The results are evaluated by comparing the ANN output with the target value using some statistical indicator.
There are several algorithms that can be used to optimise the network's functioning, such as genetic or gradient algorithms.In general, no one algorithm can be considered better than another, although, depending on the characteristics of the systems, some algorithms may achieve results quicker than others, or be more appropriate for complex network configurations.In this sense, genetic algorithms are more suitable for the latter, since they are more prone to avoid getting stuck in local solutions of the cost function.In any case, the aim of this work is to analyse the feasibility of the simplest possible solution, i.e., a network with a single hidden layer and a number of neurons to be determined.Thus, a gradient algorithm is used to optimise the parameters of the ANN during the training process.However, as the whole process is based on a stochastic process, the final optimisation of the ANN may be different, depending on the algorithm used, the execution of the training process, as well as its starting point.This means that the training must be performed repeatedly until an acceptable solution is reached.Additionally, aiming to avoid an overfitting of the ANN towards the training data, optimal neuron numbers are pending further adequacy analysis, based on the results obtained for the test data.
As mentioned above, this work considers two different approaches to characterise the performance of a Stirling engine, depending on the input variables used.Consequently, two different groups of ANNs are presented, depending on the number of variables.ANN models in the G1 and G2 groups have, respectively, three and eight different nodes in the input layer, one for each input variable.To facilitate the comparison with linear regressionbased models, all the networks in this paper have an output layer with a single neuron.

Experimental Data
Reasonably reliable performance data has been found for the 34 engines shown in Table 1, which comprise a representative variety of existing Stirling engines. 1 Free-piston engine.

Data for Analysis of Models in Group G1
Appendix A shows the 89 operating points corresponding to the engines listed in Table 1, which have been used as the basis for developing the G1-group models.Table A1 includes 23 operating points, corresponding to the first 11 engines listed in Table 1, operating with air, nitrogen, helium, and hydrogen, at mean pressures in the range of 1-150 bar, with heater wall temperatures ranging from 368 to 1173 K.These data have been used to build the correlations of the G1-1R and G1-2R models, and as training data for the development of the G1-3ANN and G1-4ANN models.Since the proportion of available data is not equal for all ranges of operating variables, this selection has been made to attempt to avoid bias towards any kind of engine, and to leave a significantly larger percentage of data than usual for testing purposes.Table A2 includes 66 operating points, used to validate the correlations of the G1-1R and G1-2R models, and as test data for the development of the G1-3ANN and G1-4ANN models.The data correspond to 29 engines, identified in Table 1 as No. 2-5, 7-9, and 12-33, being the last 11 operating points among those used as the basis for West's correlation.The ranges of variation of the variables used in the models of this group are probably the widest considered so far, as can be seen in Table 2.It should be noted that inconsistencies could be found in the literature regarding the volume used to express the dimensionless power developed by a Stirling engine.In general, both the symbols and the terms used may lead to confusion.For example, V sw is the most commonly used symbol for both the volume swept by the working piston and the amplitude of the total volume.Both values are the same for beta or gamma Stirling engines, but differ for alpha engines.To avoid ambiguity, the symbol V 0 is used in this article to designate the amplitude of the total volume.

Data for Analysis of Models in Group G2
Appendix B shows 54 operating points, corresponding to the 11 engines identified in Table 1 as No. 2-7, 10, 12, 13, and 34, which have been used as the basis for developing the G2 group models.Table A3 includes 18 operating points from nine engines, working with air, nitrogen, helium and hydrogen, at mean pressures in the range of 1-150 bar, with heater wall temperatures ranging from 403 to 1173 K.These data have been used to construct correlations for the G2-5R to G2-8R models, and as training data for the development of the G2-9ANN to G2-12ANN models.Table A4 includes 36 operating points from four engines, which have been used to validate the correlations of the G2-5R to G2-8R models, and as testing data for the development of the G2-9ANN to G2-12ANN models.Tables 3 and 4 show the ranges of the variables used in the models.

Statistical Indicators
Various criteria have been proposed to assess the fit between the experimental data and model estimates, but none of them is free of limitations [48].Since the percentage results facilitate interpretations, in this article, the performance of the models is assessed using dimensionless statistical indicators, namely the relative root mean square error RRMSE, the relative mean bias error RMBE, and the coefficient of determination R 2 .The RRMSE value is used as an optimisation criterion for both least squares regression-based models, and for training ANN-based models.To facilitate comparisons with results from other authors, the Nash-Sutcliffe coefficient of model efficiency, NSE, and the normalised values of the root mean square error and mean bias error, NRMSE and NMBE, respectively, have also been calculated, using the mean values of the experimental data as references for the latter two.R 2 values range from 0 to 1, and NSE values range from −∞ to 1. Negative NSE values suggest that the mean of the measurements is a better predictor than the model estimates themselves [49].Some authors warn against the risk of identifying the meaning of R 2 and NSE [50], while others try to avoid confusion by classifying the former as an indicator of dispersion and the latter as an indicator of overall performance [51].For example, while R 2 = 0.8 indicates that the model explains 80% of the variance in the observed data, the value NSE = 0.8 has a very different meaning, i.e., that the model mean squared error represents 20% of the observed variance.The cause of possible confusion may be that R 2 can be interpreted as a maximum potential value for NSE, as the following equation reduces to NSE = R 2 for the optimal case of __ Taylor diagrams have been recommended for visualising results in model analyses [52].In this article, this type of chart is used as an alternative to the usual scatterplots, because it provides a concise statistical summary of how well patterns match each other in terms of their correlation, their root-mean-square difference, and the ratio of their variances.The Taylor diagram is based on the definition of the centred pattern RMS difference E [52], which is related to key statistical indicators by means of the following equations: Equation ( 22) and the law of cosines are the basis for the graphical representation of the degree of closeness between a model and the reference data set, using dimensionless variables (Figure 2).Equation ( 22) and the law of cosines are the basis for the graphical representation of the degree of closeness between a model and the reference data set, using dimensionless variables (Figure 2).

Regression-Based Models of Group G1
Figure 3 provides a picture of the dispersion between the 89 dimensionless maximum brake power data points included in Tables A1 and A2, and the predictions based on Equations ( 9)- (11).

Regression-Based Models of Group G1
Figure 3 provides a picture of the dispersion between the 89 dimensionless maximum brake power data points included in Tables A1 and A2, and the predictions based on Equations ( 9)- (11).

Regression-Based Models of Group G1
Figure 3 provides a picture of the dispersion between the 89 dimensionless maximum brake power data points included in Tables A1 and A2, and the predictions based on Equations ( 9)- (11).Based on Equation ( 13), the following correlation has been obtained, which fits the 23 operating points included in Table A1   Based on Equation ( 13), the following correlation has been obtained, which fits the 23 operating points included in Table A1 with RRMSE = 22.2% and R 2 = 0.7847: This equation fits the 89 operating points included in both Tables A1 and A2 with RRMSE = 22.2% and R 2 = 0.7940.Figure 4 facilitates comparisons between data and model predictions.Based on Equation ( 15), the following correlation has been obtained, which fits the 23 operating points included in Table A1 with RRMSE = 30.1% and  0.9569: This equation fits, with RRMSE = 25.3% and  0.9598, the 82 total operating points that result after excluding the FPSE data, because of their different operating principle, and operating point No. 84, because its  , * value is significantly higher than those known for engines with a kinematic drive mechanism, from Tables A1 and A2.Based on Equation (15), the following correlation has been obtained, which fits the 23 operating points included in Table A1 with RRMSE = 30.1% and R 2 = 0.9569: This equation fits, with RRMSE = 25.3% and R 2 = 0.9598, the 82 total operating points that result after excluding the FPSE data, because of their different operating principle, and operating point No. 84, because its N * MA,max value is significantly higher than those known for engines with a kinematic drive mechanism, from Tables A1 and A2. Figure 5 facilitates comparisons between data and model predictions.Based on Equation ( 15), the following correlation has been obtained, which fits the 23 operating points included in Table A1 with RRMSE = 30.1% and  0.9569: This equation fits, with RRMSE = 25.3% and  0.9598, the 82 total operating points that result after excluding the FPSE data, because of their different operating principle, and operating point No. 84, because its  , * value is significantly higher than those known for engines with a kinematic drive mechanism, from Tables A1 and A2. Figure 5 facilitates comparisons between data and model predictions.

ANN-Based Models of Group G1
The optimal number of neurons in the hidden layer is usually determined statistically by some variant of the mean square error and the coefficient of determination [53].In this paper, the optimisation is based on RRMSE and R 2 , whose weighting is obtained by the normalised centred pattern RMSE, E n .Table 5 and Figure 6 show that the optimal number corresponds to 7 and 10 neurons for the G1-3ANN and G1-4ANN models, respectively.In addition, Figure 7 allows comparisons to be made between the data and the predictions of the models configured with the optimal number of neurons.

Regression-Based Models of Group G2
Model G2-5R is based on the following equation, which fits the 18 available experimental data points with RRMSE = 3.20% and RMBE = −0.17%: Figure 8 facilitates comparisons between data and model predictions.

Regression-Based Models of Group G2
Model G2-5R is based on the following equation, which fits the 18 available experimental data points with RRMSE = 3.20% and RMBE = −0.17%: Figure 8 facilitates comparisons between data and model predictions.Model G2-6R is based on the following equation, which fits the 18 experimental data points with RRMSE = 6.35% and RMBE = 0.43%: Figure 9 facilitates comparisons between data and model predictions.

Regression-Based Models of Group G2
Model G2-5R is based on the following equation, which fits the 18 available experimental data points with RRMSE = 3.20% and RMBE = −0.17%: (26) Figure 8 facilitates comparisons between data and model predictions.Model G2-6R is based on the following equation, which fits the 18 experimental data points with RRMSE = 6.35% and RMBE = 0.43%: Figure 9 facilitates comparisons between data and model predictions.Model G2-7R is based on the following equation, which fits the 18 experimental data points with RRMSE = 6.39% and RMBE = −0.44%: Figure 10 facilitates comparisons between data and model predictions.Model G2-7R is based on the following equation, which fits the 18 experimental data points with RRMSE = 6.39% and RMBE = −0.44%: Figure 10 facilitates comparisons between data and model predictions.
Model G2-7R is based on the following equation, which fits the 18 experimental data points with RRMSE = 6.39% and RMBE = −0.44%: (28) Figure 10 facilitates comparisons between data and model predictions.Model G2-8R is based on the following equation, which fits the 18 experimental data points with RRMSE = 10.93% and RMBE = −5.15%: Figure 11 facilitates comparisons between data and model predictions.
Figure 11 facilitates comparisons between data and model predictions.

ANN-Based Models of Group G2
Table 6 and Figure 12 allow the selection of optimal neuron numbers for the G2-9ANN to G2-12ANN models, while Figure 13 provides comparisons between the data and the predictions of the models with optimal neuron numbers.

ANN-Based Models of Group G2
Table 6 and Figure 12 allow the selection of optimal neuron numbers for the G2-9ANN to G2-12ANN models, while Figure 13 provides comparisons between the data and the predictions of the models with optimal neuron numbers.

Discussion and Future Work
The ability of a model to accurately predict the performance of a set of Stirling engines depends on the influence between the input and output variables of the model, the representativeness and quality of the data sample used for training and testing, and the mathematical procedure of the model.High accuracy cannot be expected from the models in group G1, because they only use three input variables, although they are undoubtedly influential on the outputs.Consistent with these observations, it follows from Figures 3, 4 and 7a that the G1-1R and G1-3ANN models are able to estimate  , with relative errors within ±25% for most of the Appendix A operating points, while the classical correlations are not adequate for a significant number of data.Table 7 and Figure 14 facilitate comparisons using statistical indicators.

Discussion and Future Work
The ability of a model to accurately predict the performance of a set of Stirling engines depends on the influence between the input and output variables of the model, the representativeness and quality of the data sample used for training and testing, and the mathematical procedure of the model.High accuracy cannot be expected from the models in group G1, because they only use three input variables, although they are undoubtedly influential on the outputs.Consistent with these observations, it follows from Figures 3, 4 and 7a that the G1-1R and G1-3ANN models are able to estimate ζ B,max with relative errors within ±25% for most of the Appendix A operating points, while the classical correlations are not adequate for a significant number of data.Table 7 and Figure 14 facilitate comparisons using statistical indicators.
With respect to particular operating points, for the G1-1R model, the largest percentage deviations are observed for the Ecoboy and Yamanokami-1 engines and, to a lesser extent, for the Karabulut-2 engine.It is interesting to note that this model predicts ζ B,max acceptably for the FPSE analysed, including the magnetic micro-actuator of only 10 mW and 0.05 cm 3 , values that are well outside the range of operating points used to build the model.The G1-3ANN model fits the training data somewhat better, but is less accurate for the testing data.The largest percentage deviations are observed for the Ecoboy engine, for several FPSEs and, especially, for the 2-cylinder Genoa engine.With respect to particular operating points, for the G1-1R model, the largest percentage deviations are observed for the Ecoboy and Yamanokami-1 engines and, to a lesser extent, for the Karabulut-2 engine.It is interesting to note that this model predicts  , acceptably for the FPSE analysed, including the magnetic micro-actuator of only 10 mW and 0.05 cm 3 , values that are well outside the range of operating points used to build the model.The G1-3ANN model fits the training data somewhat better, but is less accurate for the testing data.The largest percentage deviations are observed for the Ecoboy engine, for several FPSEs and, especially, for the 2-cylinder Genoa engine.
Regarding G1-model estimates of  , * , it is clear from Figure 5 that the G1-2R model gives results with relative errors within ±25% for most of the data, after excluding the FPSEs and Sunpower's Rice Husk engine, as mentioned before.The largest percentage deviations are observed for the Mitsubishi&Daihatsu and Karabulut-1 engines and, to a lesser extent, for several operating points of the Karabulut-2 engine, and for the point No. 35 of the MP102C engine in Table A2, previously associated with a measurement error [27].The estimates of the G1-4ANN model with 10 neurons in the hidden layer are acceptable for the training data, but worse for the testing ones, due to high percentage deviations for several operating points at low values of  , * , as shown in Figure 7b.However, it should be noted that such low values of  , * indicate that the degree of development of the gas circuit of the corresponding engines can be improved [24], i.e., the coefficients of indicated power losses are high, as can be seen from Equations ( 5) and (6).Table 8 shows the statistical indicators obtained for both G1-2R and G1-4ANN models.
In summary, in group G1, the proposed regression-based models are slightly better than the ANN-based models, with balanced relative errors for the total training and test data in Appendix A.  Regarding G1-model estimates of N * MA,max , it is clear from Figure 5 that the G1-2R model gives results with relative errors within ±25% for most of the data, after excluding the FPSEs and Sunpower's Rice Husk engine, as mentioned before.The largest percentage deviations are observed for the Mitsubishi&Daihatsu and Karabulut-1 engines and, to a lesser extent, for several operating points of the Karabulut-2 engine, and for the point No. 35 of the MP102C engine in Table A2, previously associated with a measurement error [27].The estimates of the G1-4ANN model with 10 neurons in the hidden layer are acceptable for the training data, but worse for the testing ones, due to high percentage deviations for several operating points at low values of N * MA,max , as shown in Figure 7b.However, it should be noted that such low values of N * MA,max indicate that the degree of development of the gas circuit of the corresponding engines can be improved [24], i.e., the coefficients of indicated power losses are high, as can be seen from Equations ( 5) and (6).Table 8 shows the statistical indicators obtained for both G1-2R and G1-4ANN models.
In summary, in group G1, the proposed regression-based models are slightly better than the ANN-based models, with balanced relative errors for the total training and test data in Appendix A.
With respect to G2-model estimates of ζ ind,max , both models G2-5R and G2-9ANN provide very accurate estimates for the training data, especially the G2-9ANN model.The results are also good for the testing data, with almost all relative deviations within ±10%, as shown in Figures 8 and 13a, and Table 9.Moreover, both models show the aforementioned anomaly of operating point No. 33 in Table A4.It is interpreted that the quality of these models is mainly due to the inclusion of the dimensionless quasi-static work per cycle ζ 0 among the input variables, since the ratio ζ ind,max /ζ 0 is bounded, as indicated by Equation ( 4).Accuracy is also enhanced by the inclusion of additional variables, reflecting the influence of dead volume and regenerator geometry.
Regarding G2-model estimates of N MA,max , Figure 9 shows that the predictions of the G2-6R model have a good degree of accuracy, with practically all relative errors within ±10%, again excluding point No. 33 in Table A4.
In contrast, Figure 13b shows that the predictions of the G2-10ANN model with two neurons in the hidden layer, although acceptable for the training data, show high relative deviations, not only as expected for the operating point No. 33 of Table A4, but also for the operating point No. 54.The statistical indicators listed in Table 10 show that the G2-6R model is somewhat better than the G2-10ANN model, especially for operating point No. 54.However, deviations for operation point No. 54 can be considered less relevant, as this prototype was never built [46], and the corresponding values of the input variables ζ 0 , r 2 hR L R /V 0 , ∑ µ dx and N p are outside the range of the training data variables and at the extremes of the test data range.In summary, it seems that a neural network with a more complex structure would be necessary to outperform the accuracy of the N MA,max predictions obtained by the regression-based model.
With respect to G2-model estimates of ζ B,max , Figure 10 and Table 11 show that the predictions of the G2-7R model fit most of the data used for its construction, but the deviations are relatively large for several of the test data.The highest percentage deviations are for some points of the MP102C engine operating at 873 K and the ST05G engine, although it should be noted that the latter data are from simulations.The relative deviations of points No. 33 and No. 54 are again outside the limits of ±10%.A4.
The results are similar using the G2-11ANN model, as shown in Table 11 and Figure 13c, apart from the 400 HP/cylinder prototype, for which this model predicts higher relative deviation.These observations show that both types of models can be acceptable for preliminary design, although their estimation capability is limited by not including input variables related to the mechanical losses of indicated power.Future improvements can be expected from the inclusion of input variables representative of mechanical efficiency, the availability of more operational data, and the use of neural networks with more complex structure.
In the case of the G2-model estimates of N * MA,max , Figure 11 shows that the predictions of the G2-8R model are accurate for most of the data used to construct the correlation, as they are virtually all included within the ±10% deviation limits.For the test data, in general, the model predictions are biased below the target data.Apart from operating point No. 33 in Table A4, the largest percentage deviations, in the order of −25%, correspond to points in the ST05G engine series.
Furthermore, Figure 13d shows that the predictions of the G2-12ANN model with two neurons in the hidden layer are acceptable for the training data, but show high relative deviations for operating point No. 54 of Table A4 and, to a lesser extent, for several points in the ST05G engine series.In order to improve the results, the G2-12ANN model was modified by arranging one neuron in the hidden layer, with the statistical effects shown in Table 12, depending on whether or not operating points No. 33 and 54 are included.Figure 15 shows the results with the modified model.As in the case of the  , models in group G2, future improvements of the G2-8R and G2-12ANN models could be based on more complexly structured neural networks and on seeking reliable training and test data, whose ranges of variation are similar once expressed in dimensionless variables.
All in all, the most pressing need for future work is to increase the amount of available experimental data.Uncertainty of existing data can be addressed by different techniques already successfully used in artificial intelligence and deep learning, such as clustering, backpropagation, Benders decomposition or option value [54][55][56].Furthermore, variants of the procedure used in this article for eight-input models could be interesting for other energy conversion systems, since Equation ( 2) is applicable to other thermodynamic cycles [27].

Conclusions
Regression-based models and simple ANNs have been developed and compared to estimate the peak power values for feasibility studies and the preliminary design of Stirling engines.The main findings are summarised below:  Models using three dimensionless input variables fit the 89 operating data points from 34 engines, with relative errors of about ±25%.The classical correlations are not adequate for a significant number of data, and do not provide criteria for estimating the speed corresponding to the peak power points;  Estimates of the dimensionless maximum indicated power from models using eight dimensionless input variables match practically all the operating data, with relative errors within ±10%, for both regression-based and ANN-based models.The regres-  As in the case of the ζ B,max models in group G2, future improvements of the G2-8R and G2-12ANN models could be based on more complexly structured neural networks and on seeking reliable training and test data, whose ranges of variation are similar once expressed in dimensionless variables.
All in all, the most pressing need for future work is to increase the amount of available experimental data.Uncertainty of existing data can be addressed by different techniques already successfully used in artificial intelligence and deep learning, such as clustering, backpropagation, Benders decomposition or option value [54][55][56].Furthermore, variants of the procedure used in this article for eight-input models could be interesting for other energy conversion systems, since Equation ( 2) is applicable to other thermodynamic cycles [27].

Conclusions
Regression-based models and simple ANNs have been developed and compared to estimate the peak power values for feasibility studies and the preliminary design of Stirling engines.The main findings are summarised below:

•
Models using three dimensionless input variables fit the 89 operating data points from 34 engines, with relative errors of about ±25%.The classical correlations are not adequate for a significant number of data, and do not provide criteria for estimating the speed corresponding to the peak power points; • Estimates of the dimensionless maximum indicated power from models using eight dimensionless input variables match practically all the operating data, with relative errors within ±10%, for both regression-based and ANN-based models.The regression-based model estimates are also quite accurate for the corresponding dimensionless engine speeds, while the ANN-based model predicts acceptable relative errors for the training data, but larger deviations for some test data; • Estimates of dimensionless maximum brake power from models using eight dimensionless input variables are also acceptable for data used in building regression models or training ANN models, but the deviations are somewhat less accurate for test data.
The predictions of the corresponding engine speeds show similar relative deviations for both regression-based and ANN models.Such deviations are related to the lack of input variables involved in the mechanical efficiency, whose inclusion is recommended for future works.

Figure 1 .
Figure 1.Architecture of an artificial neural network.

Figure 1 .
Figure 1.Architecture of an artificial neural network.

Figure 3 .
Figure 3.Comparison between classical correlations of dimensionless maximum brake power.

Figure 3 .
Figure 3.Comparison between classical correlations of dimensionless maximum brake power.
facilitates comparisons between data and model predictions.

Figure 6 .
Figure 6.Taylor diagrams to compare training results as a function of the number of neurons: (a) G1-3ANN model; (b) G1-4ANN model.

Figure 14 .
Figure 14.Comparison between the performances of the  , models of group G1 based on: (a) training data; (b) all data.

Figure 14 .
Figure 14.Comparison between the performances of the ζ B,max models of group G1 based on: (a) training data; (b) all data.

•
Dead volume V dx , wetted surface A wx , and hydraulic radius r hx of each gas circuit space x; • Physical properties of the working fluid, i.e., adiabatic coefficient γ, specific constant R, and viscosity µ; • Thermal diffusivity α R , volumetric specific heat capacity ρ R c R , and volumetric porosity ¶ V of the regenerator material; • Frequency of rotation, n s .

Table 1 .
Stirling engines with data analysed in this article.

Table 2 .
Ranges of variables in models of group G2.

Table 3 .
Ranges of ordinary variables in models of group G2.

Table 4 .
Ranges of dimensionless variables in models of group G2.

Table 5 .
Statistical results during training for ANN-based models in group G1 as a function of the number of neurons in the hidden layer.

Table 6 .
Statistical results as a function of the number of neurons during training for ANN-based models in group G2.

Table 6 .
Statistical results as a function of the number of neurons during training for ANN-based models in group G2.

Table 7 .
Statistical indicators of the performance of the  , models of group G1.

Table 7 .
Statistical indicators of the performance of the ζ B,max models of group G1.

Table 8 .
Statistical indicators of the performance of the  , * models of group G1.

Table 8 .
Statistical indicators of the performance of the N * MA,max models of group G1.

Table 9 .
Statistical indicators of the performance of the ζ ind,max models of group G2.
1Excluding the operating point No. 33 of TableA4.

Table 10 .
Statistical indicators of the performance of the N MA,max models of group G2.Excluding the operating point No. 33 of TableA4.2Excluding the operating points No. 33 and No. 54 of TableA4.

Table 11 .
Statistical indicators of the performance of the ζ B,max models of group G2.Excluding the operating point No. 33 of TableA4.2Excluding the operating points No. 33 and No. 54 of Table

Table 12 .
Statistical indicators of the performance of the N * MA,max models of group G2.Using 2 neurons in the hidden layer.2Using 1 neuron in the hidden layer.3Excluding the operating point No. 33 of TableA4.4Excluding the operating points No. 33 and No. 54 of TableA4.

Table A1 .
Dimensionless brake power= P B /(p m V 0 n s ) ζ ind Dimensionless indicated power = P ind /(p m V 0 n s ) ζ 0 Dimensionless quasi-static work per cycle = W 0 /(p m V 0 ) Cont.

Table A2 .
Database of operating points for testing models of group G1.