Cybersecurity Risk Assessment of Industrial Control Systems Based on Order-α Divergence Measures Under an Interval-Valued Intuitionistic Fuzzy Environment

With the increasing deployment of network technologies in industrial control systems (ICSs), cybersecurity has become a challenge in ICSs. Cybersecurity risk assessment (CRA) plays an important role in cybersecurity protection of ICSs. However, the weights of risk indices are constants in traditional CRA methods, and they do not fully consider the requirements of risk identification. In this paper, we define a novel order-<inline-formula> <tex-math notation="LaTeX">$\alpha $ </tex-math></inline-formula> divergence measure for interval-valued intuitionistic fuzzy numbers (IVIFNs) and further develop a novel CRA approach for ICSs based on the proposed divergence measure under an interval-valued intuitionistic fuzzy environment to contribute to the research gap. First, an order-<inline-formula> <tex-math notation="LaTeX">$\alpha $ </tex-math></inline-formula> divergence measure for IVIFNs is defined considering flexibility and robustness of divergence measures with the parameter. Next, a variable weight-based CRA approach for ICSs is developed. In this approach, IVIFNs are adopted to describe evaluation values of risk indices. The weights of risk indices are variable weight vectors and they are determined by the relative divergence closeness. Integration approaches of each node and each attack path in attack-defense trees (ADTs) are proposed based on the operations of IVIFNs, and risk scores of each attack path are calculated by using the score function. Finally, we apply the proposed method to the CRA of a civil aviation fuel supply automatic control system and verify its effectiveness and advantages by comparing it with other methods. This method can dynamically adjust the weights of risk indices considering the relationship between each risk index and the highest risk, and therefore, it can more effectively recognize the highest risk of ICSs than the traditional CRA method. In addition, it can also match the risk attitude of decision-makers by adjusting the parameter <inline-formula> <tex-math notation="LaTeX">$\alpha $ </tex-math></inline-formula>.


I. INTRODUCTION
Industrial control systems (ICSs) are widely used in electric power, petroleum and petrochemical, nuclear energy, aviation, railway, water treatment and other industries. They play an important role in today's industry [1]. In recent decades, the progress of computer and network technology has promoted the development of ICSs [2]. For example, the The associate editor coordinating the review of this manuscript and approving it for publication was Hualong Yu .
Internet of Things has entered ICSs to achieve connectivity among enterprises and savings of cost. The integration of advanced information technology and ICSs can realize the remote control and monitoring of field equipment [3]. However, advanced technology not only brings some advantages to ICSs but also makes ICSs more vulnerable and subject to various network attacks. For instance, hackers and criminal organizations use the loopholes of ICSs to destroy the normal operation of ICSs in various ways and cause great impact and loss to society and the economy. We know that the VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ ''Stuxnet'' virus, which swept the world's industries in 2010, used the loopholes of Windows system and SIMATIC WinCC of Siemens to attack ICSs so that the centrifuge in Iran ran out of control, covered up the failure, and sent back to the management department with ''normal operation'' records, which resulted in misjudgment of decision-making [4]. In 2015, the virus called blackenergy attacked an energy company in the Ivano frankisvik region of western Ukraine, which led to the power failure of 225,000 families in Ukraine. Since May 12, 2017, the variant Wannacry Blackmail virus has swept the world, and some governments, airports, hospitals, gas stations and other public institutions have been attacked.
In 2018, one chip manufacturing enterprise in Taiwan was also attacked by Wannacry Blackmail virus. All kinds of documents and databases were locked, which led to the shutdown of production lines. Although the number of attacks against ICSs is small compared with that against the internet, the consequences may be disastrous. Therefore, it is very important to ensure that ICSs are protected from cyber threats [5].
Aiming at cybersecurity problems, a series of standards have been formulated. ISO/IEC 27000 series standards on information security management solve the security problems of IT systems. The ISA/IEC 62443 standard provides a flexible framework to solve and mitigate the current and future security vulnerabilities in ICSs [6]. However, the cybersecurity problems of ICSs are difficult to eradicate. First, the standards of industrial control networks are mostly open to facilitate the application of users. As a result, it is not difficult for programmers familiar with the ICS to develop targeted malicious attack codes. In addition, the technology of hackers and criminal organizations is also improving, and it is difficult for ICSs to be free from network attacks. Risk assessment can help enterprises find the weakest links of ICSs and further take corresponding measures to optimize management, equipment and control [7]. Therefore, it is necessary to evaluate the cybersecurity risks of ICSs.
Quantitative risk assessment is an effective method of CRA. Some related works have been performed to assess cybersecurity risks of ICSs by exploiting different methods. For example, Li et al. [8] used a Petri net (PN) to establish an evaluation model and proposed a dynamic impact assessment method based on the full recognition of asset knowledge. Gemini and Sabu [9] presented a CRA method of the industrial internet of things (IoT) considering the largest loss stream based on an attack graph. Wang et al. [10] developed a CRA method combining the factor analysis of information risk (FAIR) model with Bayesian networks (BNs). From the accuracy-based perspective, they performed an in-depth analysis between FAIR and FAIR-BN under different situations based on the J-divergence measure. Qin et al. [11] designed an association network (AN) to infer the probabilities of cybersecurity incidents and built an association matrix with regard to the state variables and the key security variables to evaluate the cybersecurity risks of ICSs. According to the characteristics of risk assessment of power system, Sun et al. [12] proposed an incremental variable-based state enumeration method considering safety margins. Bolbot et al. [13] devised a risk assessment method for ship systems based on cyber preliminary hazard analysis. Chang et al. [14] applied a failure mode and effect analysis model to quantify the risk level combining evidential reasoning (ER) with a rule-based Bayesian network (RBN). Jha et al. [15] presented a risk assessment framework for smart grids, which applied hardware reliability and data reliability to evaluate risks.
Due to the limited prior knowledge of attacks, it is difficult for decision-makers to accurately evaluate the probabilities of attack events. Some researchers have implemented fuzzy theory into risk assessments of diverse fields, such as manufacturing corporations [16], construction project investment [17], traffic congestion [18], buildings [19], freight transportation systems [20], ship control systems [21], mines [22], and so on. The assessment data were expressed in their favor, such as fuzzy numbers [23]- [25], [30]- [32], intervals [33], [34], Z-numbers [35], triangular fuzzy numbers [16], [19]- [22], [26], [27], [36], trapezoidal fuzzy numbers [28], Pythagorean fuzzy numbers [37]- [40], linguistic term sets [17], [29], and double hierarchy hesitant fuzzy linguistic information [18]. In recent years, researchers have studied risk assessment from different technical models, different methods and different fuzzy evaluations. Wang et al. [18] combined double hierarchy hesitant fuzzy linguistic term sets with the ORESTE method and proposed a risk assessment method on the 5S traffic congestion model. Huang et al. [23] adopted entropy weights to calculate the relative importance of element layers in the fuzzy analytic hierarchy process (FAHP) model and improved the correlation between failure modes using the gray relation analysis (GRA) method. Considering the lack of sufficient historical data, Qi et al. [25] proposed a dynamic CRA method for ICSs by extending the traditional BN to a fuzzy BN. Gul et al. [26] combined a FAHP method with fuzzy VIKOR to construct a new risk assessment framework. Gul and Celik [27] proposed a risk assessment method by incorporating a fuzzy rule-based expert system with the Fine-Kinney method and applied it in rail transportation systems. Considering decision-makers' psychological behavior, interaction relationships, and uncertainty among risk indices, Wang et al. [28] presented a hybrid failure mode and effect analysis (FMEA) framework by combining the TODIM approach with the Choquet integral method. Li et al. [29] proposed a novel FMEA model taking linguistic term sets into account in fuzzy Petri nets (FPNs), which calculated the weights of decision-makers based on the TOPSIS method. Yu et al. [31] applied the cloud model to elaborate the risk indices, and the risk indices were integrated by the MAX-MIN operator. Ultimately, the method provided the risk levels under different situations, and a detailed and in-depth discussion was made. From the author's point of view, the randomness of the cloud model is very strong, which will lead to different calculation results even though the same input values are given. Tian et al. [33] advanced a risk assessment method considering intervals with self-confidence. In the approach, the weights of decisionmakers were calculated by combining the subjective weights with the objective weights, and the fuzzy inference laid the foundation for IF-THEN rules. Onari et al. [35] introduced the Z-number theory to risk assessment. In this approach, the risk indices including severity, occurrence, and detection, are expressed by Z-numbers. However, the operations for Z-numbers are difficult, and we need to convert Z-numbers to other fuzzy numbers for further processing. Wang et al. [36] implemented attack-defense tree models (ADTs) for CRA of an airport automatic fuel supply control system based on fuzzy theory. Akram et al. [38] provided a risk evaluation under a Pythagorean fuzzy environment using hybrid TOPSIS and the ELECTRE I method. Recently, some novel Pythagorean fuzzy interaction aggregation operators based on Archimedean t-conorm and t-norm (ATT) have been developed [40], and applying them to CRA of ICSs will be a good subject.
Although existing fuzzy risk assessment methods manage uncertainty to some extent, Atanassov and Gargov [41] introduced the concept of the interval-valued intuitionistic fuzzy set (IVIFS), which is more powerful in expressing uncertainty. The operations for the IVIFS have been defined [42]- [44]. The fundamental characteristic of the IVIFS is that the values of its membership degree and nonmembership degree are intervals rather than exact numbers. Therefore, the IVIFS is finer and smoother for representing the fuzzy evaluation information than the fuzzy set (FS) and intuitionistic fuzzy set (IFS). Recently, it has still attracted the focus of scholars [45]- [51]. Liu et al. [45] obtained variable weights of attributes by integrating the accuracy function and the subjective weights; on this basis, a novel multi-attribute group decision-making (MAGDM) approach was developed with IVIFS. Garg and Kumar [46] presented some exponential distance measures using the connection number (CN) of IVIFS. Kumar and Chen [47] proposed a novel score function of CN based on set pair analysis (SPA) theory under an interval-valued intuitionistic fuzzy information environment. Che et al. [48] constructed a new entropy measure in the context of IVIFS by introducing the proposed distance function. Zhou et al. [49] raised the ratio comparison rules of IVIFS. Bustince et al. [50] introduced some new similarity measures between interval-valued intuitionistic fuzzy numbers (IVIFNs) considering the width of intervals and admissible orders. Deveci et al. [51] proposed a new combinative distance-based ASsesment (CODAS) method based on Taxicab distance and the largest Euclidean distance between the IVIFNs. Information measure is an important topic in fuzzy decision-making problems. Bhandari and Pal [52] defined the fuzzy divergence measure and provided a new fuzzy measure method from the perspective of probability distribution. Since then, scholars have conducted extensive research on this topic [53]- [63]. In recent years, Wang and Wan [58] transformed IVIFS into IFS, defined a possibility degree of IVIFS and developed a divergence measure for IVIFNs based on the proposed possibility degree. Li et al. [59] presented a new crossentropy measure with parameters for IVIFNs based on the J-divergence measure. Mishra et al. [60] defined the Jensenlogarithmic divergence measure and the Jensen-exponential divergence measure for IVIFNs. Mishra et al. [61] defined some novel entropy and divergence measures, and applied them in the process of multicriteria service quality assessment combined with the TODIM method. Song et al. [62] introduced a divergence-based cross entropy measure with parameters for AIFSs, which is the intuitionistic fuzzy set defined by Atanassov and generally is also called IFS. Verma [63] advanced some new order-α divergence measures between two IFSs.
By combining the literature on risk assessment and IVIFNs, we are committed to developing a new CRA method for ICSs and a new fuzzy divergence measure under an interval-valued intuitionistic fuzzy environment for the following reasons: 1. The weights of risk indices are constants and do not vary with actual situations in the existing risk assessment methods. The weights of risk indices are determined by decision-makers considering their experience or determined by the entropy weight method. However, the constant-weighting approaches are unreasonable in the following situation. For example, it is easy to ignore the risk in this situation when the value of an index is very close to the highest risk and yet its weight is very small. However, the risks of ICSs easily cause serious consequences. Therefore, it is necessary to develop a new assignment method for the weights of risk indices.
2. FS only reflects the membership degree of belonging to one level. IFS embodies the membership degree of belonging to one level and the nonmembership degree of not belonging to one level. For FS and IFS, their membership degrees and nonmembership degrees are expressed by exact numbers. However, due to the limitations of the decision-maker's experience and uncertainty of risk indices, it is difficult for decision-makers to evaluate risk indices belonging to one level and not belonging to one level with exact numbers. IVIFS can better describe uncertainty because its membership degree and nonmembership degree are intervals.
3. The divergence measure with the parameter has the characteristics of flexibility and robustness, yet there is a lack of research on it for IVIFNs. Therefore, to address these issues, we define a novel orderα divergence measure for IVIFNs, and on this basis, we put forward a variable weight-based CRA method for ICSs under an interval-valued intuitionistic fuzzy environment. It arises from the ADT model. First, a novel divergence measure for IVIFNs is defined, and the weights of risk indices are calculated based on the proposed divergence measure. Second, a novel CRA method is established. Finally, the proposed method is applied to the CRA of a civil aviation fuel supply automatic control system, and its effectiveness is verified. VOLUME 10, 2022 The main contributions of this paper are as follows: 1. We define a divergence measure with the parameter for IVIFNs. It makes up for the gap that there is no divergence measure with parameters for IVIFNs. 2. We expand IVIFS to the CRA of ICSs considering the power of IVIFS and formulate integration approaches of all nodes and attack paths with IVIFNs in the ADT model. 3. We develop a novel CRA method for ICSs. In this method, the weights of risk indices can vary with actual situations in the risk assessment process and are calculated by using the proposed divergence measure.
The rest of this article is organized as follows: Section II reviews the concept of the ADT model and theory about IVIFS. Section III defines a novel order-α divergence measure for IVIFNs. Section IV describes the framework and implementation process of our proposed approach. In addition, it introduces the risk assessment scales with IVIFNs, the determination method for the weights of risk indices, and the integration expressions of all nodes and attack paths based on the operations of IVIFNs. Section V introduces a case involving the CRA of a civil aviation fuel supply automatic control system and makes a comparative analysis with other methods. Section VI summarizes our work, states the contributions and limitations of the proposed method, and expounds further work.

A. ATTACK-DEFENSE TREES
In 1999, Schneier [64] proposed attack trees (ATs) as a tool to evaluate the security of complex systems. Considering the limitation of ATs that they cannot reflect the interaction between attacks and defenses, Kordy et al. [65] extended ATs to attack-defense trees (ADTs), which constitute a graphical expression of the actions that attackers might take to attack one system and the defenses that defenders can adopt to protect the system. An ADT model is a tree-like graph of an attack scenario, as shown in Fig. 1.
It consists of one root node, attack leaf nodes, defense leaf nodes and combination nodes. The attack target is defined as the root node, the attack leaf nodes are the actions to be implemented by attackers, and defense leaf nodes are the measures to be taken by defenders. In Fig. 1., G0 represent the root node, L1, L2, and L3 represent the attack leaf nodes, and D1, D2, and D3 represent the defense nodes. In addition, combination nodes consist of AND nodes and OR nodes. AND nodes indicate different steps to attack the same goal, and OR nodes are alternatives. Any path from leaf nodes to the root node indicates a complete attack process to realize the attack target. All attack paths can be generated by traversing the whole attack tree.

B. THEORY ABOUT INTERVAL-VALUED INTUITIONISTIC FUZZY SET
Definition 1 [41]: Let X be a nonempty set, and then the interval-valued intuitionistic fuzzy set (IVIFS) is defined as whereµ L A (x) and µ H A (x) represent the upper and lower bounds of the membership degree, respectively, and represent the upper and lower bounds of the nonmembership degree, and π H A (x) represent the upper and lower bounds of the hesitation degree, respectively, π L for convenience. Definition 2 [42], [43]: be any two IVIFNs, then the operations of IVIFNs are defined as Definition 3 [66]: be an IVIFN, and then its score function S (A) is defined as

III. A NOVEL DIVERGENCE MEASURE FOR IVIFNs A. EXISTING DIVERGENCE MEASURES
Wang and Wan [58] gave a standard definition of divergence measures for IVIFSs below. Definition 4 [58]: Let X be a finite universe, and let IVIFSs (X ) be the set of all IVIFSs on X. A mapping D: IVIFSs (X ) × IVIFSs (X ) → [0, 1] is a standard divergence measure for IVIFSs if for every A, B ∈ IVIFSs (X ) which has the following properties: First, we review some recent divergence measures for IVIFNs.
Wan and Wang [58]: Li et al. [59]: . Mishra et al. [60]: Mishra et al. [61]: Now, we advance two examples to illustrate the weaknesses of the above developed divergence measures for IVIFNs. The calculation results are shown in Table 1.
It can be seen from Table 1 that the divergence value between A and B is 0 by using the divergence measure in [58]   it is not suitable to use them to calculate the divergence, and they do not meet the property DP1 of the divergence measure for IVIFNs. The divergence calculated by using the divergence measure D M 1 in [60] and the divergence measure D M in [61] is reasonable. However, we note that the divergence measure with parameters for IVIFNs is still a gap, which can provide better flexibility and robustness in real decision-making problems. Therefore, in the following, we are determined to develop a divergence measure with parameters for IVIFNs.

B. A NOVEL DIVERGENCE MEASURE FOR IVIFNs
First, we present three Lemmas in which we provide some support for developing the order-α divergence measure for IVIFNs.
Lemma 1: Let two finite discrete probability distributions be P = p 1 , p 2 , · · · , p t and Q = q 1 , q 2 , · · · , q t , 0 ≤ p k , q k ≤ 1 for k = 1, 2, · · · , t, and and the equality in (3) holds if and only ifp k = q k , ∀k. Proof of Lemma 1: Consider the function φ (x) = x 1−α for every x ∈ [0, ∞]. When the parameter α satisfies the condition 0 ≺ α ≺ 1, the function φ is concave. Then, according to the Jensen inequality, we can obtain: where Proof of Lemma 2: Based on the definition of IFS, we obtain the following: Therefore, we conclude that f I (A, B) ≤ 1 by Lemma 1, and the equality holds if and only if A = B, that is, where

Proof of Lemma 3:
Based on the definition of IVIFS, we obtain the following: be any two IVIFNs, and then we define the order-α divergence measure between two IVIFNs A and B given by

by Lemma 2, and the equality holds if and only if
To satisfy the symmetry of divergence measures, we define as an order-α divergence measure for two IVIFNs Their divergence values are also 1 in these cases by using (6) and (7). In other words, the maximum divergence is 1 between two IVIFNs A and B.
Property DP3 in fact holds when we define the order-α divergence measure for IVIFNs.
In summary, the order-α divergence measure in (7) is a standard divergence measure for IVIFNs.
Theorem 1: For all A, B, W ∈ IVIFS, the divergence measureD (A B ) satisfies the following properties: Their proofs are straightforward by the operations and the definition of the order-α divergence measure for IVIFNs. Hence, we omit their proofs from here. Now, we solve two examples in Table 1 again with our proposed measure, and the obtained divergence values are shown in Table 2. The results presented in Table 2 clearly show that the divergence values are all not 0 by using the proposed divergence measure when A is not equal to B. Hence, the proposed measure is a valid and flexible divergence measure for IVIFNs, which can be employed to handle some problems related to various application fields.

IV. THE PROPOSED CRA APPROACH FOR ICSs
The framework of our proposed variable weight-based CRA approach for ICSs is shown in Fig. 2. The architecture of our approach consists of three main aspects: establishment of a decision matrix, determination of VOLUME 10, 2022 variable weights, and risk assessment of ICSs. The establishment of a decision matrix includes the division of the risk level and the evaluation values of each attack accident given by experts. Evaluation values are expressed in the form of IVIFNs. The weights of risk indices are determined by the proposed divergence measure for IVIFNs. To ensure that the risks of ICSs can be better identified, the weights of risk indices are changed in the decision-making process, which is different from traditional CRA methods for ICSs. When the decision matrices and the weights of all attack accidents and defense strategies are obtained, all leaf nodes and each possible attack path can be integrated according to the operations for IVIFNs, and then, we can assess the highest risk with the score function of IVIFNs.
The implementation process of the CRA approach is shown in Table 3.

A. RISK ASSESSMENT SCALES
The risk assessment scales in this study are from [36]. We transform triangular fuzzy numbers (TFNs) into IVIFNs considering the same quantitative results. The quantization results of TFNs are calculated by the method in [67], in which the decision-maker's attitude is neutral. The scores of IVIFNs are calculated by using (2). The risk assessment scales in the form of IVIFNs are shown in Table 4.

B. DETERMINATION OF THE WEIGHTS
Since the maximum IVIFN is ([1, 1] , [0, 0]), the highest risk of leaf nodes A + in ICSs is defined as Next, we present the relative divergence closeness between the leaf nodes L ij and the highest risk A + given by (8) where i is the number of leaf nodes and j is the number of risk indices.
The greater the relative divergence closeness is, the farther the leaf node is from the highest risk, and the smaller its weight should be. Conversely, the weight should be larger. The normalized weight of the leaf nodes L ij is defined as where ω ij ∈ [0, 1] , n j=1 ω ij = 1.

C. THE CALCULATION OF THE RISK ASSESSMENT
We define the integration expressions of attack leaf nodes, defense leaf nodes, combination nodes and attack paths in the ADT model based on the operations of IVIFNs.

1) INTEGRATION OF THE NODES a: INTEGRATION OF ATTACK LEAF NODES
Suppose the evaluation values of the attack cost, attack difficulty and detection possibility for the ith attack leaf node are , respectively, and its comprehensive evaluation value Z A i is defined as where ω A i ,cos t , ω A i ,diff and ω A i ,det represent the weight of the attack cost, attack difficulty and detection possibility for the ith attack leaf node, respectively, and

b: INTEGRATION OF DEFENSE LEAF NODES
Suppose that the evaluation values of the defense cost, defense difficulty and defense time for the ith defense leaf node are µ L  evaluation value Z D i is defined as where ω D i ,cos t , ω D i ,diff and ω D i ,time represent the weight of the defense cost, defense difficulty and defense time for the ith defense leaf node, respectively, and sent the comprehensive value of n attack leaf nodes, and then the integration value Z i of the AND node is defined as The integration value of the OR node Z i is defined as

d: THE OPERATION BETWEEN THE ATTACK NODES AND DEFENSE NODES
Let represent the comprehensive value of one attack leaf node, and let represent the comprehensive value of one defense leaf node. Therefore, the operation between the attack nodes and defense nodes is defined . The comprehensive value of the attack path is expressed as

3) DETERMINATION OF RISK SCORES
Let p = µ L p , µ H p , ν L p , ν H p be the comprehensive value of any one attack path; and then the risk score of this attack VOLUME 10, 2022 path is obtained by (12).

V. APPLICATION AND ANALYSIS OF THE PROPOSED METHOD A. APPLICATION
To demonstrate the application of the proposed cybersecurity risk assessment method, we consider the case discussed in [36]. The civil aviation fuel supply automatic control system is mainly divided into three logical levels: the enterprise management layer, process monitoring layer and field control layer. The network structure is shown in Fig. 3. Fig. 4 is an ADT model against the civil aviation fuel supply automatic control system. The attack goal in the system is gaining access to the SCADA system. This is the root node in the ADT model. The possible risk events in the system include hardware failure, operator errors, Havex malware, Lnk file vulnerability, Printer Service Vulnerability, and MS08-067 vulnerability, which is to call the NetPath-Canonicalize function in the server service program through the MSRPC over the SMB channel, causing the stack buffer to overflow to obtain Remote Code Execution, U disk injected with virus, Wincc vulnerability, Denial of service attack, replay attack and eavesdropping attack are the attack leaf nodes in the ADT model. The defense measure adopted in the system is to set up a firewall, which is the defense leaf node in the ADT model.
The CRA process of the proposed method for the system is shown in Fig. 5.  The calculation process of the cybersecurity risk assessment for the civil aviation fuel supply automatic control system is as follows: Step 1: The evaluation values of attack leaf nodes and the defense leaf node are obtained by the method in Section 5.1, as shown in Table 5 and Table 6.
Step 2: The weights of all attack leaf nodes and the defense leaf node are obtained by the method in Section 4, respectively, as shown in Table 7 and Table 8 (α = 0.5).
Step 3: Comprehensive values of all leaf nodes are calculated by using (10) and (11). The calculation results are shown in Table 9.
Step 5: Scores of all attack paths are obtained by using (16). The results are shown in Table 10. Table 10 shows that the score of attack path AP2 is the highest, that is, the risk caused by operation error is the highest. The score of hardware failure is the second highest. The calculation results show that the key factor affecting the    securities of the civil aviation fuel supply automatic control system is internal security vulnerability. Measures should  be strengthened in the training of personnel operation and maintenance of aircraft.

B. SENSITIVITY ANALYSIS OF THE PARAMETER
In what follows, the influence of different values of the parameter α on the risk ranking results of the system is discussed. The risk scores of all attack paths and the ranking results under different parameter α in the interval [0.1, 0.9] are shown in Table 11.
From Table 11, we notice that the risk scores of all attack paths increase with the increase of the parameterα. If the attitude toward the risk level is conservative, set the parameter αto be small. If the attitude toward the risk levels is neutral, the parameterα is set to 0.5. If the attitude toward the risk levels is aggressive, the parameter α is set to be large. However, the risk ranking results varying with the parameter α are the same, and the attack paths of the highest risk are all AP2; that is, the changes in the parameter α do not have an effect on the risk ranking results, which indicates that the proposed order-α divergence measure for IVIFNs has good robustness.

C. COMPARISIONS AND DISCUSSION
To further prove the effectiveness of the proposed method in this paper, we still solve the CRA problem of the civil aviation fuel supply automatic control system by using these methods in [36], [45] and [61]. Triangular fuzzy numbers are adopted to describe evaluation values of risk indices and the weights of risk indices are determined by subjective weighting methods in [36]. IVIFNs are adopted to describe evaluation values of risk indices and the weights of risk indices are determined by the hybrid weighting method in [45]. IVIFNs are adopted to describe evaluation values of risk indices and the weights of risk indices are determined by the objective weighting method based on entropy in [61]. Among them, the weights of risk indices assigned by decision-makers are w=(0.500,0.333,0.167) in [36]. The range of the parameterγ is −12.27 ≤ γ ≤ 8.60 according to evaluation values of this example in [45]. The weights of risk indices obtained by using the entropy method are w=(0.287,0.340,0.373) in [61]. Ranking results under different methods are shown in Table 12. Table 12 shows that the highest risk is all AP2 by using these methods. The ranking results in this paper are fully consistent with those in [45] (γ = 0) and almost identical to those in [45] (γ = 12.87 and γ = 8.60) and [61], except that the rankings of AP3, AP4 and AP6 are slightly different. The method in [36] cannot distinguish the risks of AP5 and AP6. Therefore, our proposed method is effective, and it can be used to handle the cybersecurity risk assessment problems of ICSs.
To further prove the advantages of the developed method in this paper, we give a counterexample to show that the ranking results of risks based on the existing approaches are  unreasonable in some cases, while the proposed method can better identify risks.
Example 1: Suppose an expert evaluates the risk level of three attack leaf nodes (L1, L2, and L3) with regard to risk indices, including attack cost (G1), attack difficulty (G2) and detected probability (G3) in the form of IVIFNs, and their evaluation values are shown in Table 13.
We still utilize the above four methods to rank the risks, and the results are shown in Table 14. In this example, the range of the parameter γ is −13.54 ≤ γ ≤ 11.21 in [45].
From Table 14, we find that the highest risk is L1 by using the proposed method, and the ranking results of risks are the same as those in [61]. However, the highest risk is L2 by using the methods in [36] and [45]. Table 11 shows that for leaf node L1, the evaluation value of attack cost , which is very high and close to the highest risk. From the perspective of risk identification, the detected probability of leaf node L1 should be assigned a large weight, while the weight of attack cost should be assigned a small weight, so we allocate the weight of detected probability to 0.454 (α = 0.10), 0.535 (α = 0.50), and 0.573 (α = 0.90). At the same time, we allocate the weight of attack cost to 0.213 (α = 0.10), 0.132 (α = 0.50), and 0.094 (α = 0.90). However, the weight of attack cost assigned to all leaf nodes is always 0.500, and the weight of detected probability is always 0.167 in [36]. The weights of each index do not change with the evaluation information in the decision process; in that way, the identification of risks is easily confused when the evaluation of attack cost is relatively safe and the evaluation of detected probability is relatively dangerous. Although the weights of each index change with the evaluation information in the decision process in [45], the adjustments of weights depend on the discrimination between individual information and mean information with the parameter γ . As shown in Table 14, we find that the weight of attack cost is 0.220 (γ = −13.54), 0.500 (γ = 0), 0.732 (γ = 11.21), and the weight of detected probability is 0.073 (γ = −13.54), 0.167 (γ = 0), 0.244 (γ = 11.21) in the leaf node L1. The weights calculated by the method in [61] are 0.164, 0.557, and 0.279, respectively. Obviously, the assignment of weights is unreasonable and inconsistent with the actual situation in [36], [45] and [61]. Consequently, the proposed method fully considers the requirements of risk identification and can lead to a more reasonable risk assessment result.
Through the above discussion and analysis, we summarize the advantages of our methods compared with the other methods in detail below.
1. Different from traditional risk assessment methods, the weights of risk indices can be adjusted in the decision-making process by using our proposed method. They are determined by their closeness to the highest risk. Our proposed method can effectively recognize the highest risk because it fully considers the requirements of risk identification. 2. Our approach has good robustness and flexibility. Its ranking result of risks is stable and does not change with the change of the parameter α. In addition, decision-makers can also choose different values of the parameter α according to the different risk attitudes.

VI. CONCLUSION
In this paper, we define an order-α divergence measure for IVIFNs and develop a novel CRA method for ICSs under an interval-valued intuitionistic fuzzy environment. Finally, we apply the proposed method to the CRA of a civil aviation fuel supply automatic control system, verify its effectiveness and demonstrate its advantages. In summary, our research has three main contributions: 1. We define an order-α divergence measure for IVIFNs, which can make up for the gap that there is no divergence measure with the parameter for IVIFNs. 2. We expand IVIFS to the CRA of ICSs and formulate integration approaches of all nodes and attack paths with IVIFNs in the ADT model. 3. We propose a novel CRA method for ICSs. In our method, we regard the weights of risk indices as variable weight vectors and develop a new technology to determine the weights of risk indices based on the proposed divergence measure. The proposed method can effectively avoid irrationality in the results of risk assessment compared with traditional CRA methods.
However, the proposed method also has its limitations: 1. It is only applicable to CRA problems in which evaluation values are expressed in the form of IVIFS. 2. The risk indices are regarded as independent without considering their mutual interaction in the integration process of leaf nodes.
In future research, we are committed to the following two aspects: 1. We will extend other fuzzy sets to the CRA of ICSs, for instance, interval-valued q-rung orthopair fuzzy sets. 2. We will advance some novel CRA methods for ICSs considering the interactions among risk indices.