*3.1. Analytic Hierarchy Process (AHP)*

The evaluation process of HRS site selection is a multi-index decision-making process, and different indices have different influences, so they should be accordingly weighed [13]. Methods to determine the weight value are usually AHP, the entropy weight method (EWM) and the Criteria Importance Though Intercrieria Correlation (CRITIC)method [37,42,57]. Researching the literature showed that AHP is usually used for criteria weighting [61]. It is rated, and the weight is calculated by using the relative size information of the value between two pairs of comparison. The entropy method calculates the weight of the amount of information of the data entropy value. The CRITIC method mainly uses data fluctuations or the correlation between data to calculate a weight [13,42,44,62].

Because the site selection of HRSs involves many indices, these indices lack accurate data descriptions. Through hierarchical analysis, the importance of the two indices in a multilevel system can be compared. Therefore, an expert questionnaire was designed in this study; 20 respondents from the hydrogen energy industry, scientific research institutions, and doctoral students in universities were selected to fill in the questionnaire and conduct detailed consultations to conduct a graded evaluation of our evaluation indicators, using AHP method to determine the most practical and reliable coefficients of each index.

The specific steps of AHP are as follows:

(1) Project weight calculation

An evaluation index for classification was set up for the calculation of the average values of experts for indices at all levels; it was used as the initial weight, and Table 4 Saaty 1–9 scale method was referenced to quantify [42]; the judgment matrix was constructed.

**Table 4.** Scales and meanings of the judgment matrix [42].


The product root method is used to calculate the geometric mean value of each row of judgment matrix *wi*:

$$\overline{w\_i} = \left(\prod\_{j=1}^n a\_{ij}\right)^{\frac{1}{n}} \tag{1}$$

Among them, *aij* is the relative importance of the index *i* and the index *j* compared in pairs, and *n* is the number of indices. Feature vectors *wi*, namely, the weight coefficients of each index, are obtained by normalizing the geometric mean values of each line [42,61].

#### (2) Consistency test of the judgment matrix

Calculate the maximal eigenvalue of the judgment matrix:

$$\lambda\_{\text{max}} = \frac{1}{n} \sum\_{i=1}^{n} \frac{\left(\sum\_{j=1}^{n} a\_{ij} w\_j\right)}{w\_i} \tag{2}$$

Calculating consistency indicators (*CI*). The specific formula is:

$$CI = \frac{\lambda\_{\text{max}} - 1}{n - 1} \tag{3}$$

Calculating consistency ratio (*CR*, where *RI* is the average random consistency index [47]:

$$CR = \frac{CI}{RI} \tag{4}$$

When *CR* < 0.10, the consistency of the judgment matrix is considered to be acceptable; otherwise, the judgment matrix should be appropriately modified [63,64]. The consistency of the judgment matrix in this paper was tested in MATLAB R2018a, and the consistency of the judgment matrix was less than 0.1, which was acceptable.

On the basis of expert judgments and statistical data, the constructed judgment matrices are shown in Tables 5–8.

**Table 5.** Judgment matrix U–U*i*.


Note: CI = 0.0368, CR = 0.0634, r = 3.0735.

**Table 6.** Judgment matrix U1–U1*j*.


Note: CI = 0.0231, CR = 0.0398, λ = 3.0005.

**Table 7.** Judgment matrix U2–U2*j*.


Note: CI = 0.0091, CR = 0.0158, λ = 3.0183.



Note: CI = 0.0024, CR = 0.0381, λ = 3.0000.

(3) Calculate the combined weight of each index

The combined weight of each index is calculated according to the weights of the criterion layer and the sub-criteria layer [63]. Taking the combined weight as an example, *A*<sup>11</sup> = *w*1*w*11, and other combined weights are analogous, so the combined weight A of the evaluation index system is:

*A* = (0.2255, 0.0883, 0.1152, 0.1280, 0.1491, 0.0977, 0.0393, 0.0589, 0.0982)

The combined weight not only takes into account the weight distribution of the secondary indicators in the indicators of this level, but also takes into account its overall evaluation. The weight distribution in the index system, the larger the combined weight, the more important the secondary index is in the whole evaluation index system [61].

#### *3.2. Principle of Fuzzy Comprehensive Evaluation(FCE) Method*

Fuzzy comprehensive evaluation (FCE) can reasonably combine these attributes and factors with items affected by multiple attributes or their overall strengths or weaknesses [43,48]. Therefore, the FCE method based on a fuzzy set evaluates subordinate-level status evaluation from multiple indicators [41]. The evaluation interval is, on the one hand, divided into changes, so that the evaluation criteria and the fuzziness of the influencing factors can reflect the level of the considered objects; on the other hand, the experience of people can be given full play in the evaluation, so that evaluation results are more objective and in line with the actual situation [46,63], FCE can achieve the combination of qualitative and quantitative aspects, expand the amount of information, thereby improving the amount of credible evaluation conclusions [43].

The steps of establishing the FCE model of HRSs site selection are as follows:

(1) A set of indices affecting the site selection of HRSs was established:

$$\mathcal{U} = \{\boldsymbol{u}\_1, \boldsymbol{u}\_2, \dots, \boldsymbol{u}\_n\} \tag{5}$$

There are 9 indices in *U*: In order to facilitate the following calculations, redefine *u*<sup>1</sup> construction scale; *u*<sup>2</sup> investment strength, *u*<sup>3</sup> operational costs, *u*<sup>4</sup> hydrogen production technology, *u*<sup>5</sup> hydrogen storage technology, *u*<sup>6</sup> transportation technology, *u*<sup>7</sup> the population density, *u*<sup>8</sup> environmental factors, and *u*<sup>9</sup> social identity.

(2) The evaluation set of influencing factors of hydrogenation station siting was established:

$$V = \{v\_1, v\_2, \dots, v\_m\}\tag{6}$$

Correctly determining the membership function is the basis of applying fuzzy theory to properly and quantitatively describe a fuzzy concept, and the key to solving various practical problems by using the fuzzy method. The corresponding membership function is selected according to different evaluation factors. Determine the weight of each evaluation factor:

$$A = \{a\_1, a\_2, \dots, a\_i\} \tag{7}$$

*ai* is the weight of the first factor, ∑*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *ai* = 1.

(3) Evaluate each factor *ui* to determine the degree to which it belongs to comment, and the evaluation result is written as *ri*:

$$r\_i = \{r\_{i1}, r\_{i2}, \dots, r\_{\prime}, r\_{im}\} \tag{8}$$

In this paper, 10 relatively authoritative experts were selected to evaluate all indices, and the expert evaluation level was converted into the specific quantitative value of each index of each HRS according to the assigned value. Due to the large number of experts, the average value was taken as the final evaluation result of a factor. Thus, fuzzy relation R from U to V could be determined by all single factor evaluation results [44,48].

R is the fuzzy comprehensive evaluation matrix.

$$R = \begin{bmatrix} r\_{11} & r\_{12} & \cdots & r\_{1m} \\ r\_{21} & r\_{22} & \cdots & r\_{2m} \\ \dots & \dots & \dots & \dots & \dots \\ r\_{n1} & r\_{n2} & \cdots & r\_{nm} \end{bmatrix} \tag{9}$$

(4) Calculate the value of *B*:

$$B = A \bullet R$$

This can be used as fuzzy comprehensive evaluation vector after the normalization of *B* obtained by the operation, where *A* is the index weight. The FCE result of HRSs is then obtained [44,48].

(5) Composition operator selection:

A basic step of the FCE of HRSs site selection is to determine the evaluation grade.

$$V = \{v\_1, v\_2, \dots, v\_m\} = \{\text{"very good"}, \text{"good"}, \text{"general"}, \text{"p}or \text{"}, \text{"very poor"}\}.$$

The FCE can synthesize and average the functions of various influencing factors according to the weight, which has obvious advantages in the comprehensive evaluation treatment of various influencing factors. In the calculation process, two types of fuzzy operators are often used, namely, principal factor prominence and weighted average [45]. The salient feature of the main factor operator is to reduce the interference of other data, while the salient feature of the weighted average operator is to avoid information loss. This can take into account the functions of various influencing factors as a whole. In view of the large number of factor sets in this study, in order to avoid information loss and make the collected data information play its role to the maximum, the weighted average fuzzy operator synthesis model is adopted [48,63].

#### *3.3. Artificial Neural Networks (ANN) Theory*

An artificial neural network (ANN) is an empirical model that mimics the function of a human neural network. There is often a highly nonlinear mapping relationship between its input and output, and it is generally difficult to write its expression, so it is called "black box" [28].The characteristic of ANN is that it can store information or knowledge distribution in a large number of neurons or the whole system. It has the potential of self-learning and self-organization [65]. In addition, it has strong fault tolerance and can deal with noisy or incomplete data [27,66]. Therefore, ANN provides a powerful tool in solving complex nonlinear multi-index comprehensive evaluation problems. The BP (Back propagation) neural network model is a multi-layer feedforward neural network trained according to the error back propagation algorithm, which is the most widely used neural network. So in this study, the back propagation algorithm (BP algorithm) is used, and the BP neural network model has the advantages of fast calculation speed and efficient solution. It can better simulate the process of comprehensive evaluation by evaluation experts [28,37]. In order to simulate the process of comprehensive evaluation by experts, considering the input, hidden and output layers, a BP neural network is designed. The comprehensive evaluation results of the 9 evaluation indices of 50 HRSs in coastal and major cities in China were taken as the input of the ANN, and the FCE results was obtained by fuzzy logic calculation are used as the output to train the neural network. The architecture is shown in Figure 4.

**Figure 4.** Architecture of artificial neural network for HRSs site selection.

ANN is an intelligence-based processing structure consisting of interconnected processing units called neurons. It includes input nodes, output nodes and one or more layers of hidden node. The input layer to the output layer needs to be processed by the hidden layer. The neurons of each hidden layer basically perform two tasks: (1) the weighted sum of all process inputs, and (2) the weighted sum of the nonlinear transformation of the neuron transfer function produces the output of each neuron. Common transfer functions include hyperbolic tangent S-type (Tansig), S-type (Logsig), positive linear (Poslin), and purelin linear (Purelin) transfer functions, as shown in Formulas (11)–(14) [28,67]. The output layer is used for use error estimates to predict outcomes. Initially in the backpropagation algorithm, the input is propagated to the hidden layer, which propagates the sensitivity back to reduce errors. At the end of the process, it updates the weights and biases. The performance index of the neural network is generally expressed by mean square error (*MSE*), mean relative error (*MRE*) and mean absolute error (*MAE*), which are used to judge the errors between network output values and target values, as shown in Formulas (15)–(17) [35,36].

$$\tan \text{sig}(n) = \frac{e^n - e^{-n}}{e^n + e^{-n}} \tag{11}$$

$$\log \text{sig}(n) = \frac{1}{1 + e^{-n}} \tag{12}$$

$$\text{poslim}(n) = \begin{cases} \text{ } n \text{ } n > 0\\ \text{ } 0 \text{ } n \le 0 \end{cases} \tag{13}$$

$$
overline{n}(n) = n\tag{14}$$

$$MSE = \frac{1}{n} \sum\_{t=1}^{n} \left( y - y' \right)^2 \tag{15}$$

$$MAE = \frac{1}{n} \sum\_{t=1}^{n} \left( y - y' \right) \tag{16}$$

$$MRE\% = \frac{100}{n} \sum\_{t=1}^{n} \left| \frac{y - y'}{y} \right| \tag{17}$$

where *n* represents the number of training sets, *y* represents the output result of the target, and *y* represents the output result of the neural network.

The basic idea of establishing an intelligent location decision model of HRSs was as follows.

Determine the evaluation index set, and the number of indices determines the input node of the neural network. If the number of specific evaluation indices is 9, the number of input nodes of neural network is 9.

Each index in the index set is the input factor of the neural network. According to the characteristics of each index, the input factor is graded, the input value of the neural network is determined by the expert evaluation method, and the expected output value of the artificial neural network is obtained by the FCE method. After grading each factor, each score value is normalized into the score value in the (0,1) domain to meet the requirements of the BP neural network range [28,66].

#### **4. Prediction of HRSs Site Selection Based on FCE-ANN**

In this study, 50 HRSs already in operation in coastal and major cities of China and 8 HRSs that are being planned in Shanghai were selected as analytical examples. According to the established evaluation index system for site selection of HRSs, the FCE-ANN model is verified by an example through the steps of model establishment, data analysis and calculation, and result comparison and analysis.

### *4.1. Overview of HRSs in Operation in China*

At present, coastal and major cities are the main developed areas for HRSs construction in China. The specific distribution is shown in Figure 5 [10], highlighted by the red dots. It can be seen that Guangdong is the province with the largest number of HRSs currently in operation. In order to better reflect the construction of HRSs in China, this paper selects 50 HRS from the provinces marked in the figure as research objects.

**Figure 5.** The Main Distribution of HRSs in China [10].

The collected information showed that the hydrogen storage, transportation, and hydrogenation of operational HRSs in China are mostly the same, and they were all constructed and are maintained at a relatively low cost. However, there are relatively few HRSs built with high costs and relatively advanced technologies such as hydrogen production and pipeline transportation. Among the 50 selected HRSs, only Shanghai YiLanShunGong has adopted hydrogen production and pipeline transportation in the industrial zone.

#### *4.2. Fuzzy Comprehensive Evaluation of Site Selection of HRSs*

In this study, 10 experts were invited to form an expert evaluation group, of which 3 were from the hydrogen energy industry with technical backgrounds such as hydrogen production and hydrogen storage for more than 5 years; 5 were from Chinese universities and other research institutions (including 2 doctoral students), specializing in HRSs research, and very concerned about the construction and planning of HRSs in China; the other two are government officials, familiar with the planning policies, population density and economic level of HRSs in various regions. With our organization, the expert evaluation group received as much research information as possible from 50 HRSs in China one week in advance. The formal evaluation is in a conference room. First, 10 experts interpret and discuss the information of 50 HRSs according to their respective areas of expertise. Secondly, experts conduct single-index evaluation of the 9 evaluation indices of each HRS. Finally, collect the evaluation grades of various evaluation indices of all HRSs by experts, and assign the quantitative results to each index of each HRS according to the membership function, and finally take the average value as the final evaluation result of a certain index, and the evaluation levels and corresponding assignments are shown in Table 9.

**Table 9.** Corresponding values of evaluation levels.


Among these indices, the evaluation requirements of each index are different. For example, when the cost of hydrogen is relatively low, the hydrogen cost level of the hydrogenation station is relatively high, both in construction cost and operation and maintenance cost. Among technical factors, the more advanced the technology of each index is, the higher the corresponding evaluation level is. For example, when the hydrogen production technology is more advanced, the hydrogen production technology of the HRS has a relatively high level, and the level of hydrogen storage technology is the same as that of transportation technology. Among the social influencing factors, the higher the population density, the more target users for the hydrogenation station, so the level of the population density index of the hydrogenation station is relatively high. However, for environmental factors, the less environmental damage caused, the higher the level. The higher the safety index is, the higher the corresponding evaluation level.

The values of 50 HRSs based on the evaluation levels of various indices of operating HRSs were collected and sorted here. The comprehensive evaluation results and neural network training data of HRSs are shown in Appendix A (Table A1) after the results had been arranged in order.

The above AHP calculation confirms that the combined weight *A* of each index in the sub-criteria layer to the target layer in the evaluation system is:

*A* = (*A*1, *A*2, *A*3)=(0.2255, 0.0883, 0.1152; 0.128, 0.1491, 0.0977; 0.0393, 0.0589, 0.0982)

The 50 HRSs were divided into 5 groups with 10 HRSs in each group. The 10 former HRSs were taken as examples.

The comprehensive results of economic factors are:

*B*1 = *A*1·*R*1 = (0.2255, 0.0883, 0.1152)· ⎡ ⎣ 0.98 0.98 0.98 0.93 0.98 0.93 0.91 0.91 0.94 0.93 0.96 0.96 0.98 0.98 0.97 0.98 0.93 0.96 0.80 0.98 0.98 0.98 0.88 0.91 0.97 0.93 0.93 0.96 0.96 0.93 ⎤ ⎦ = (0.4187, 0.4187, 0.4089, 0.4011, 0.4184, 0.4034, 0.3945, 0.4006, 0.3932, 0.4034)

The comprehensive results of technical factors are:


The comprehensive results of social factors are:

*B*3 = *A*3·*R*3 = (0.0393, 0.0589, 0.0982)· ⎡ ⎣ 0.71 0.93 0.88 0.98 0.83 0.94 0.76 0.94 0.89 0.96 0.91 0.91 0.91 0.88 0.87 0.88 0.88 0.91 0.86 0.88 0.91 0.91 0.88 0.88 0.90 0.88 0.88 0.88 0.90 0.88 ⎤ ⎦ = (0.1709, 0.1795, 0.1746, 0.1768, 0.1722, 0.1752, 0.1681, 0.1770, 0.1740, 0.1760)

So, the overall results of the top 10 are:

$$B = A \cdot R = [A1, A2, A3] \cdot \begin{bmatrix} R1 \\ R2 \\ R3 \end{bmatrix} = (0.9524, 0.9467, 0.9395, 0.939, 0.9313, 0.9223, 0.9037, 0.9221, 0.9172, 0.9205) \cdot R$$

Similarly, the fuzzy comprehensive evaluation results of the other 40 HRSs are calculated in the same way. Please refer to Appendix A (Table A1) for the specific results.

Comprehensive evaluation results of the above 50 HRSs were carried out in accordance with the comprehensive evaluation results reorder in Appendix A (Table A1). The last column sorting result shows that YiLanShunGong had the highest evaluation levels in the 50 HRSs, and the site of the construction scale, operational cost, hydrogen source, and the impact on the society were the most reasonable.

#### *4.3. ANN Training and Prediction*

#### 4.3.1. ANN Construction and Training

In this study, a BP neural network was modeled and simulated in the MATLAB software. The BP (Back propagation) neural network model is a multi-layer feedforward neural network trained according to the error back propagation algorithm, which is the most widely used neural network. The backpropagation algorithm is slow to converge and sometimes leads to overfitting. To address these issues, Levenberg–Marquardt (LM) backpropagation was developed to achieve fast convergence without overfitting [35]. By reduction in error response of the network became more smother and it also reduces the problem of overfitting. LM back propagation algorithm uses the conjugate gradient technique to reduce the sum of squares at each titration. The FCE results of the 9 evaluation indices of 50 HRSs in major cities in China were taken as the input of the ANN, and the ranking of FCE results was obtained by fuzzy logic calculation are used as expected output layer for constructing and training the neural network [27].

The number of hidden-layer neurons affects the performance and accuracy of a neural network. If the number of hidden layer neurons is too small, the data fitting of the neural network may not be carried out well. If the number of neurons is too large, the neural network may perform well in training, but poorly in the prediction results of unknown samples. Therefore, it is very important to determine the appropriate number of neurons. In order to ensure the best model performance, different neural networks with from 3 to 10 neurons using Tansig transfer function were tested, and the neural network model with the minimal *MAE*, *MRE*, *MSE* was selected [35,36]. Test results are shown in Table 10.


**Table 10.** Influence of number of hidden-layer neurons on neural network using Tansig transfer function.

As shown in Table 10, when the number of neurons is 7, the *MAE*, *MRE*, *MSE* reach the minimums, and the results are 0.0289, 0.1129 and 0.463×10−<sup>7</sup> espectively. Therefore, the number of neurons in this case model is 7. In addition to the number of neurons in the hidden layer affecting the performance of the neural network, different transfer functions in the hidden layer also affect the performance of the neural network. Figure 6 shows the *MAE*, *MRE*, *MSE* of different transfer functions in the hidden layer. It can be seen that the Purelin transfer function of this neural network structure has the best effect. The number of epochs is set 2000 so that the neural network can fully iterate. The three-layer BP neural network can approximate any nonlinear continuous function with arbitrary precision. For the small amount of data in this paper, the number of hidden layers is set to 1. The finally relevant parameters of this ANN model are shown in Table 11.

**Figure 6.** The *MAE*, *MRE*, *MSE* of different transfer functions in the hidden layer.


The neural network code was constructed in MATLAB R2018a, and the data in Appendix A (Table A1) were substituted into the neural network for training. The results of expert evaluation were taken as the output to obtain the structure of the neural

network. The correlation coefficient of neural network is shown in Figure 7. The samples number for training, validation, and testing data are 34 (about 70%), 8 (about 15%), and 8 (about 15%), respectively. The training set of a neural network learns the dataset of samples and constructs the required model by training the sample data. The function of the validation set is to stop the training of the neural network model in time when the performance of the model fails to continue to improve or change during training. Using ANN, the process of training is prone to overfitting and underfitting. The test set does not participate in the training and mainly tests the accuracy of the training model and can detect whether the model produces overfitting. If the results of the test set are too large, an overfitting may occur. The Figure 7 shows that the regression coefficients of the training, validation, and test sets of the model were all close to 1, indicating that the ANN model was well-trained and there is no phenomenon of underfitting and overfitting. So the ANN model could be used for the site selection prediction of HRSs.

**Figure 7.** Correlation coefficients between model objectives and ANN prediction based on training, validation, and test sets, and the whole dataset.

4.3.2. Location Prediction of HRSs

After training the ANN model, 8 HRSs that are being planned in Shanghai are selected for prediction, and the accuracy of the model is evaluated. The FCE results and ranking of the 8 HRSs are shown in Appendix A (Table A2). The scores of each index of {D1,D2,D3,D4,D5,D6,D7,D8} of the 8 HRSs were input into the ANN model for prediction and evaluation. The result of its operation is as follows:

*YP*= (0.9682, 0.9493, 0.9473, 0.94516, 0.9319, 0.8767, 0.8322, 0.8138)

A comparison of the predicted and evaluation results is shown in Table 12.

Table 12 shows that the expert evaluation results obtained by FCE were very close to the prediction result of the ANN model, the ANN prediction results are the same as the comprehensive evaluation. The comprehensive ranking of the eight HRSs was consistent. It shows that the ANN model we constructed can predict the location of the HRS well, and further shows that there is no underfitting and overfitting phenomenon. Among the two evaluation methods, the ShuiChan Road HRS was the best site selection result, which shows that HRSs site selection can be decided quickly and accurately through ANN model. According to the comprehensive evaluation results of experts, the values of the top four HRSs are relatively close, and the results were discussed and analyzed. According to the actual survey data, the planned daily hydrogen refueling capacity of these four HRSs is 500kg and the refueling pressure is 35Mpa. Except for ShangChai HRS, which is a fixed station, the other three HRSs are Skid-mounted stations. The skid-mounted station has a certain advantage in terms of construction scale and operation and maintenance cost, while the fixed station invests heavily. Therefore, the four HRSs have little difference in economic factors. Similarly, in terms of technical factors, experts generally believe that the ShuiChan Road HRS is located on the main road in BaoShan District, and the transportation is more convenient, which makes the cost of hydrogen relatively low. At the same time, the surrounding population density is relatively large, the demand for HRSs is relatively large, and the public acceptance is relatively high. The two HRSs planned in the lower-ranked JiaDing District are oil-hydrogen joint construction stations, which are far inferior to other HRSs in terms of major economic factors. Moreover, Jiading District has a low population density and is far from the core area [68]. The transportation cost is high, and the hydrogen storage technology in the station is far less than that of the independent HRS, so they rank low. After the above discussion and actual investigation, it is worthy of affirmation that the comprehensive score of the HRS of the ShuiChan road is the highest, which also proves the objectivity of the expert evaluation and the accuracy of the ANN.


**Table 12.** Comparison table of predicted sample results.

#### **5. Conclusions**

The evaluation index system of HRSs site selection is a complex nonlinear system, which needs to comprehensively consider economic, technical, and social factors. Therefore, it is particularly important to establish an efficient and accurate mathematical model to select the construction sites of hydrogen stations. Therefore, this paper provides a new site selection model for HRSs based on FCE-ANN, which evaluates the operational HRSs in China through the AHP and the FCE method, and uses the obtained data to train the ANN model. The input of the ANN model is the evaluation index results for site selection of HRSs, and the output is the FCE results of the site selection of HRSs. The trained ANN model is used to predict the site selection planning of Shanghai HRSs, and it is concluded that the ShuiChan Road HRS is the best site selection in Shanghai planning. The FCE-ANN model prediction results were compared with the results obtained by the FCE method. The results of the two methods were consistent, clearly showing that a FCE-ANN can objectively and accurately evaluate alternatives to reduce the influence of human factors, and that a FCE- ANN model can quickly and accurately solve HRSs site selection decision-making problems. It can also provide reference for scholars to study other site selection decision-making issues. According to the evaluation system of HRSs site selection

established in this paper, a comprehensive understanding of HRSs can be obtained, and on the basis of the evaluation results, it can provide reference for government decision-makers in formulating macro-control policies. In the evaluation system, economic factors have the largest weight, which also reflects that in the process of investment and construction of hydrogen refueling stations, it is very important to evaluate the economic strength of the investor and constructor. Currently, for China, the hydrogen energy industry is in its infancy, and investment in construction a HRS has a relatively long return on investment period, so the government should also give maximum policy and financial support, so as to promote the healthy and sustainable development of the hydrogen energy supply chain industry. This study uses the FCE-ANN model to solve the complex multi-factor evaluation problem, but requires sufficient historical data for model training. However, there are not many hydrogen refueling stations in China, and the historical data and evaluation indicators available for reference are limited. The problem can be better solved in the future. The HRSs evaluation index system established in this paper is based on the current development stage of the hydrogen energy industry. In the future, with the development of the hydrogen energy industry, the evaluation index system of the HRS should be continuously improved and supplemented. It is suggested to explore the deep influencing factors from the perspective of the public's recognition and satisfaction of HRSs. At the same time, the influence of uncertain factors on HRSs in the future should be considered, and fuzzy mathematics and mathematical logic should be combined with ANN to solve problems such as uncertainty and fuzziness, so as to consider the development strategy of HRSs in the long run.

**Author Contributions:** Conceptualization, X.Q. and Y.Z.; methodology, Y.Z.; software, C.L.; validation, Y.Z., C.L. and J.Z.; formal analysis, Y.Z.; investigation, J.Z.; resources, J.Z.; data curation, Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, X.Q.; visualization, Y.Z.; supervision, X.Q.; project administration, X.Q.; funding acquisition, X.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Major Project of Technological Innovation in Hubei Province, grant number 2020BED010 and 2019AAA075.

**Acknowledgments:** Gratitude is expressed to all the experts and scholars who participated in the establishment of the evaluation system of this study through interviews and questionnaires; and gratitude is expressed to 10 experts and scholars who participated in the comprehensive evaluation of hydrogen refueling stations in this study. Sincere thanks also to Orange Group Database for providing the information of HRSs.

**Conflicts of Interest:** The authors declare no conflict of interest.
