A covert attack risk assessment method for power systems

Abstract

The present invention discloses a method for assessing the risk of covert attacks applied to power systems, comprising the following steps: inputting wind power and load data under typical scenarios, establishing the objective function and constraint conditions of an upper-level attack model that minimizes the attack cost; establishing the objective function and constraint conditions of a lower-level economic dispatch model that minimizes the power generation cost, and obtaining a two-level programming model of covert FDIA; processing the logical constraints in the upper-level attack model and converting them into mathematical forms; processing the lower-level economic dispatch model; converting the two-level programming model into a single-level mixed integer linear programming model, and using a solving software to calculate the single-level mixed integer linear programming model; and outputting the optimal attack vector. The present invention can effectively identify the weak links in the operation of the power system.

Description

Technical Field

[0001] The present invention relates to the technical field of information physical security of power systems, and in particular to a covert attack risk assessment method applied to power systems. Background Art

[0002] The safe and reliable operation of the power system requires monitoring the real-time data of many nodes in the network, summarizing and calculating the data at the dispatch center, and adjusting the power flow distribution in the network when necessary. With the high degree of informatization and intelligence of the power system and the rapid development of renewable energy, while bringing more uncertainties to the system, the security of communication data faces new challenges. The academic community has confirmed that if the attack vector is set properly, the false data injection attack (FDIA) can be successfully implemented and will not be identified by the existing bad data detection link.

[0003] However, existing research usually focuses on the attacker injecting false data to deceive the defender into making a user power outage strategy to ensure the safe operation of the line. Although this attack method also takes into account the operating characteristics of the power system as much as possible and limits the amount of power injection into the node, the zero power point, etc., in fact, this rare abnormal operating state is easy to be detected by the dispatcher, especially the sudden over-limit of the line, which is likely to be discovered by experienced dispatchers, and then doubt the authenticity of the data.

[0004] Therefore, there is a type of attack that focuses more on the secrecy of network attacks rather than maximizing load losses. In this type of attack, the attacker uses false data to deceive the dispatcher into believing that after issuing the dispatch instruction, the system has reached an economical and safe operating state, but in fact the system is already at a safety risk of exceeding the limit on some lines. Summary of the invention

[0005] In order to make up for the deficiencies of the above-mentioned background technology, the present invention proposes a covert attack risk assessment method applied to a power system to solve the problem that covert attacks are difficult to identify.

[0006] The technical problem of the present invention is solved by the following technical solutions:

[0007] The present invention discloses a method for evaluating the risk of stealth attacks applied to a power system, comprising the following steps: S6. Input the wind power and load data under typical scenarios, and establish the objective function and constraints of the upper-layer attack model for minimizing the attack cost; S7. Establish the objective function and constraints of the lower-layer economic dispatch model for minimizing the generation cost, and obtain the two-layer programming model of stealth FDIA; S8. Process the logical constraints in the upper-layer attack model and transform them into mathematical forms; S9. Process the lower-layer economic dispatch model; S10. Transform the two-layer programming model of stealth FDIA into a single-layer mixed-integer linear programming model, and use a solving software to calculate the single-layer mixed-integer linear programming model; S11. Output the optimal attack vector.

[0008] In some embodiments, before step S6, it further includes: S1. Input the power grid, conventional thermal power generation units, and deterministic load data; S2. Input the wind turbine generator units and uncertain load parameters, where the parameters include: the probability distribution of the wind turbine generator unit output and the probability distribution of the uncertain load; S3. Generate an initial scenario set N and perform scenario clustering; S4. Screen out k typical wind power and load scenarios; S5. Define ts as the serial number of the reduced typical scenario, and let ts = 1, that is, first take the first typical scenario generated by reduction as the calculation object; after step S11, it further includes: S12. After the calculation of the first typical scenario according to steps S6 to S11, determine whether ts < k holds, that is, whether all typical scenarios have been calculated; if ts < k, it means that there are still remaining typical scenarios not calculated, and it is necessary to make ts = ts + 1, and return to step S6, and repeat steps S6 to S11, and make a judgment once after the calculation of each typical scenario; when ts ≥ k, it means that all typical scenarios have been calculated, and the program ends running.

[0009] In some embodiments, the data in step S1 includes: branch reactance, upper and lower limits of the output of conventional thermal power generation units, generation cost, and deterministic load node data.

[0010] In some embodiments, the parameters in step S2 include: the probability distribution of the wind turbine generator unit output and the probability distribution of the uncertain load.

[0011] In some embodiments, step S3 specifically includes: using the Monte Carlo method to generate the initial scenario set N, and performing scenario clustering based on the maximum minimum distance algorithm to determine the clustering center.

[0012] In some embodiments, step S4 specifically includes: using the k-means scenario reduction method improved based on typicality to screen out k typical wind power and load scenarios; the definition of the typicality TD(s) of scenario s is shown in formula (1):

[0013]

[0014] Where, F={f1,f2,…,f i ,…,f n}, represents the allied scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes in i is the i-th element in F, n is the number of ally scenes; E = {e1, e2, …, e i ,…,e m}, represents the enemy scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes other than i is the i-th element in E, m is the number of adversary scenarios; IS(s) and ED(s) represent the internal similarity and external similarity calculated by the average Euclidean distance between scenario s and the ally scenario set F and the adversary scenario set E, respectively; Z ru , Z d is the normalization factor.

[0015] Furthermore, the normalization factor Z ru , Z d The specific calculation method is shown in formula (2):

[0016]

[0017] In the formula, s i and j are any two scenes in the initial scene set N; s a and b is a subset c u k is the preset number of scene reductions; u is the current scene clustering center.

[0018] In some embodiments, the calculation method of the objective function of the upper attack model in step S6 is shown in formula (3) to formula (5):

[0019]

[0020] Where N PD 、N PL are the node set and line set of the power system respectively; a i 、a j is a 0-1 integer variable, representing the attack strategy for node i and line j, and the false injection power ΔD of the node and line respectively. i , ΔPL j There is the following logical relationship:

[0021]

[0022]

[0023] Where, ΔD i , ΔPL j are the false power data injected by the attacker into node i and line j respectively.

[0024] Furthermore, the constraints of the upper-layer attack model in step S6 include: network security constraints, stealth attack constraints and attack target constraints.

[0025] Furthermore, the calculation method of the network security constraint is shown in formula (6) to formula (9):

[0026] PL T =H·PG–H·(PD-S) (6)

[0027] PL F =H·PG-H·(PD+ΔD-S) (7)

[0028] ΔPL=PL T -PL F =H·ΔD (8)

[0029] Where PL T The real data of the line flow, that is, the information that can be obtained from the attacker's perspective; PL F represents the false data of the line flow, that is, the information that can be obtained from both the attacker and the defender's perspective; ΔPL and ΔD represent the false data injected by the attacker into the line flow and node power; PG and PD represent the original generator output and load power of all nodes respectively; S represents the load power removed by the system; H is the conversion matrix, which represents the linear relationship between the line flow and the node power. The calculation formula is as follows:

[0030]

[0031] In the formula, B M and B S Respectively represent the mutual admittance matrix and the node self-admittance matrix; A A represents the node-branch association matrix;

[0032] The calculation method of the stealth attack constraint is shown in formula (10) to formula (13):

[0033]

[0034] -(1-η j )·PL j.lim ≤PL j.F ≤(1-η j )·PL j.lim(11)

[0035] -τ·PD i ≤ΔD i ≤τ·PD i (12)

[0036]

[0037] In the formula, S i represents the load power removed at node i; η j ∈[0,1), indicating the margin left by the dispatcher for the transmission capacity of the line (for example, 0.9); PL j.F and PL j.lim are the false power flow value and capacity limit value of line j respectively; PD i represents the original load power at node i; τ∈(0,+∞), represents the detection threshold of the defender for false data of the node, which depends on the accuracy of the detection equipment invested in the power grid;

[0038] The calculation method of the attack target constraint is shown in formula (14):

[0039]

[0040] Where PL j.T is the real power flow value of line j.

[0041] In some embodiments, the objective function of the lower-level economic dispatch model in step S7 is calculated as shown in formula (15):

[0042] min∑c g PG g (15)

[0043] In the formula, c g Represents the unit output cost of each generator set; PG g The output of each generator set.

[0044] Furthermore, the calculation method of the constraint conditions of the lower-level economic dispatch model in step S7 is shown in formulas (16) to (18):

[0045]

[0046]

[0047] 0≤S d ≤D d +ΔD d (18)

[0048] In the formula, represents the capacity limit of line l; Respectively represent the minimum output and maximum output limit of generator g; D d Represents the actual load data of node d; S d Represents the load shedding power of node d.

[0049] Further, step S8 specifically includes: for the logical constraint represented by formula (14), first define a 0-1 variable δ for each line j , as a sign of whether the real power flow of the line exceeds the line capacity limit, that is:

[0050]

[0051] To achieve the purpose of the attack, the attacker needs to inject false data, and the line overload constraint that should be satisfied is:

[0052]

[0053] Using the Big M method, the logical constraints of formula (14) are transformed into the following mixed integer form:

[0054] PL j.lim -|PL j.T |≤M·δ j -ε (21)

[0055] In the formula, M represents a sufficiently large positive number, and ε represents a sufficiently small positive number;

[0056] Introduce a 0-1 integer variable δ j.ori By removing the absolute value sign in formula (21), the calculation formula is:

[0057]

[0058] At this point, the logical constraints of formula (14) are transformed into the form of formula (20) and formula (22).

[0059] Furthermore, step S8 also includes: converting formula (4) and formula (5) into the form of formula (23) and formula (24), and the calculation formula is as follows:

[0060]

[0061]

[0062] In the formula, δ D+,d , δ D-,d , δ PL+,j and δ PL-,jis an additional 0-1 integer variable, which is related to the positive and negative direction of the injected power; M represents a sufficiently large positive number, ε represents a sufficiently small positive number; τ∈(0,+∞) represents the detection threshold of the defender for false node data, which depends on the accuracy of the detection equipment invested in the power grid.

[0063] In some embodiments, step S9 specifically includes: using the Kuhn-Tucker KKT method to transform the lower-level economic dispatch model into the form of formula (25) to formula (28):

[0064]

[0065]

[0066]

[0067]

[0068] Where λ, μ are Lagrange multipliers introduced to satisfy the equality constraints; H g T , H d T are the transpose of the g-th and d-th column vectors of matrix H respectively; cs d represents the load shedding loss at node d; PG g The output of each generator set; Respectively represent the minimum output and maximum output limit of generator g; S d represents the load shedding power of node d; α l , β g , γ d , Lagrange multiplier introduced to satisfy the upper and lower limit constraints of variables; ω α ,l , ω β ,g , ω γ ,d , They are 0-1 integer variables additionally introduced into the complementary relaxation conditions corresponding to the line power, generator output and node load shedding constraints.

[0069] In some embodiments, step S10 specifically includes: converting the two-layer programming model of the secret FDIA into a single-layer mixed integer linear programming model with formula (3) as the objective function and formulas (6) to (13), (20), (22) to (28) as constraints through steps S8 and S9, and using solving software to calculate the single-layer mixed integer linear programming model.

[0070] The present invention also discloses a computer-readable storage medium, wherein the computer-readable storage medium stores computer-executable instructions. When the computer-executable instructions are called and executed by a processor, the computer-executable instructions prompt the processor to implement the steps of any of the above methods.

[0071] The beneficial effects of the present invention compared with the prior art include:

[0072] The method for assessing the risk of covert attacks applied to power systems provided by the present invention first inputs wind power and load data under typical scenarios, establishes an upper-level attack model with the goal of minimizing the attack cost of the attacker and a lower-level economic dispatch model with the goal of minimizing the power generation cost of the power grid, then converts the constructed two-level programming model into a single-level mixed integer linear programming model, and calls the solving software for solving it, thereby effectively identifying the weak links in the operation of the power system. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a flowchart of a covert attack risk assessment method according to an embodiment of the present invention.

[0074] Figure 2 It is a modified IEEE-39 node system electrical wiring diagram of an embodiment of the present invention. DETAILED DESCRIPTION

[0075] The present invention is further described below with reference to the accompanying drawings and in combination with preferred embodiments. It should be noted that the embodiments and features in the embodiments of the present application can be combined with each other without conflict.

[0076] It should be noted that the directional terms such as left, right, up, down, top, and bottom in this embodiment are merely relative concepts or are based on the normal use status of the product, and should not be considered as restrictive.

[0077] The purpose of this invention is to analyze the minimum cost (i.e. the minimum number of state quantities that need to be tampered) required for an attacker to launch a covert FDIA to cause an arbitrary line overload, as well as the impact of wind power and load uncertainty on the success rate of such attacks, in order to effectively identify the weak links in the power system.

[0078] Figure 11 is a flow chart of a risk assessment method of a confidential FDIA considering wind power and load uncertainty provided by an embodiment of the present invention, the method comprising the following steps:

[0079] S1. Basic data of input power grid, conventional thermal power units and deterministic loads, including: branch reactance, upper and lower limits of conventional thermal power unit output, power generation cost and deterministic load node data.

[0080] S2. Input uncertainty variable information, including: probability distribution of wind turbine output (Weibull distribution), probability distribution of uncertain load (Gaussian distribution).

[0081] S3. Generate an initial scene set N using the Monte Carlo method, cluster the scenes based on the maximum and minimum distance algorithm, and determine the cluster center.

[0082] S4. Improved K-means Algorithm Based on Typicality Degree (IKBTD) is used to screen out k typical wind power and load scenarios. For scenario s, its typicality degree TD(s) is defined as shown in formula (1):

[0083]

[0084] Where, F={f1,f2,…,f i ,…,f n}, represents the allied scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes in i is the i-th element in F, n is the number of ally scenes; E = {e1, e2, …, e i ,…,e m}, represents the enemy scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes other than i is the i-th element in E, m is the number of adversary scenarios; IS(s) and ED(s) represent the internal similarity and external similarity calculated by the average Euclidean distance between scenario s and the ally scenario F and the adversary scenario E, respectively. ru , Z d is the normalization factor, and the specific calculation method is shown in formula (2):

[0085]

[0086] In the formula, s i and j are any two scenes in the initial scene set N; s aand b is a subset c u Any two scenes in ; k is the preset number of scene reductions; u is the current scene clustering center.

[0087] S5. Define ts as the serial number of the typical scene after reduction. Let ts = 1, that is, the first typical scene generated by reduction is taken as the research object.

[0088] S6. Input wind power and load data under the typical scenario ts, and establish the objective function (see formula (3) to formula (5)) and constraint conditions (see formula (6) to formula (14)) of the upper attack model that minimizes the attack cost.

[0089]

[0090] Where N PD , N PL are the node set and line set of the power system respectively; a i , a j is a 0-1 integer variable, representing the attack strategy for node i and line j, and the false injection power ΔD of the node and line respectively. i , ΔPL j There is the following logical relationship:

[0091]

[0092]

[0093] Where, ΔD i , ΔPL j are the false power data injected by the attacker into node i and line j respectively.

[0094] Constraints: network security constraints (see formula (6) to formula (9)), stealth attack constraints (see formula (10) to formula (13)), and attack target constraints (see formula (14)).

[0095] PL T =H·PG–H·(PD-S) (6)

[0096] PL F =H·PG-H·(PD+ΔD-S) (7)

[0097] ΔPL=PL T -PL F =H·ΔD (8)

[0098] Where PL T The real data of the line flow, that is, the information that can be obtained from the attacker's perspective; PL Frepresents the false data of the line flow, that is, the information that can be obtained from both the attacker and the defender's perspective; ΔPL and ΔD represent the false data injected by the attacker into the line flow and node power; PG and PD represent the original generator output and load power of all nodes respectively; S represents the load power removed by the system; H is the conversion matrix, which represents the linear relationship between the line flow and the node power. The calculation method is shown in formula (9):

[0099]

[0100] In the formula, B M and B S Respectively represent the mutual admittance matrix and the node self-admittance matrix; A A represents the node-branch association matrix.

[0101]

[0102] -(1-η j )·PL j.lim ≤PL j.F ≤(1-η j )·PL j.lim (11)

[0103] -τ·PD i ≤ΔD i ≤τ·PD i (12)

[0104]

[0105] In the formula, S i represents the load power removed at node i; η j ∈[0,1), indicating the margin left by the dispatcher for the transmission capacity of the line (for example, 0.9); PL j.F and PL j.lim are the false power flow value and capacity limit value of line j respectively; PD i represents the original load power at node i; τ∈(0,+∞) represents the detection threshold of the defender for false data of the node, which depends on the accuracy of the detection equipment invested in the power grid (for example, 0.5).

[0106]

[0107] Where PL j.T is the real power flow value of line j.

[0108] S7. Establish the objective function (see formula (15)) and constraints (see formula (16) to formula (18)) of the lower-level economic dispatch model that minimizes the power generation cost, and obtain the two-level planning model of the hidden FDIA.

[0109] min∑c g PG g (15)

[0110] In the formula, c g Represents the unit output cost of each generator set; PG g The output of each generator set.

[0111]

[0112]

[0113] 0≤S d ≤D d +ΔD d (18)

[0114] In the formula, represents the capacity limit of line l; Respectively represent the minimum output and maximum output limit of generator g; D d Represents the actual load data of node d; S d Represents the load shedding power of node d.

[0115] S8. Process the logical constraints in the upper-level attack model and convert them into mathematical form.

[0116] For the logical constraint expressed by formula (14), we first need to define a 0-1 variable δ for each line. j , as a sign of whether the real power flow of the line exceeds the line capacity limit, that is:

[0117]

[0118] To achieve the purpose of the attack, the attacker needs to inject false data, and the line overload constraint that should be satisfied is:

[0119]

[0120] Using the big M method, the logical constraints of formula (14) can be converted into the following mixed integer form:

[0121] PL j.lim -|PL j.T |≤M·δ j -ε (21)

[0122] In the formula, M represents a sufficiently large positive number (for example, 10 5 ), and ε represents a sufficiently small positive number (for example, 10 -5 ).

[0123] Considering the direction of the line flow, a 0-1 integer variable δ is introduced j.ori By removing the absolute value sign in formula (21), we can obtain:

[0124]

[0125] At this point, the logical constraints of formula (14) can be converted into the form of formula (20) and formula (22).

[0126] For formula (4) and formula (5), the same processing method is used to convert them into the form of formula (23) and formula (24):

[0127]

[0128]

[0129] In the formula, δ D+,d , δ D-,d , δ PL+,j and δ PL-,j is an additional 0-1 integer variable and is related to the positive and negative directions of the injected power; M represents a sufficiently large positive number (e.g., 10 5 ), and ε represents a sufficiently small positive number (for example, 10 -5 ), τ∈(0,+∞), represents the detection threshold of the defender for false node data, which depends on the accuracy of the detection equipment invested in the power grid (for example, 0.5).

[0130] S9. Use the Kuhn-Tucker (KKT) method to transform the lower-level economic dispatch model into the form of formula (25) to formula (28).

[0131]

[0132]

[0133]

[0134]

[0135] Where λ, μ are Lagrange multipliers introduced to satisfy the equality constraints; H g T , H d T are the transpose of the g-th and d-th column vectors of matrix H respectively; cs d represents the load shedding loss at node d; PG g The output of each generator set; Respectively represent the minimum output and maximum output limit of generator g; S dIndicates the load shedding power of node d; α l , β g , γ d , The Lagrange multipliers introduced to satisfy the upper and lower bound constraints of the variables; ω α ,l , ω β ,g , ω γ ,d , Are respectively the 0-1 integer variables additionally introduced in the complementary slackness conditions corresponding to the line power, generator output, and node load shedding constraints.

[0136] S10. Transform the two-layer programming model of the stealthy FDIA into a single-layer mixed-integer linear programming model with the formula (3) as the objective function and the formulas (6) to (13), formula (20), and formulas (22) to (28) as the constraint conditions. For this model, existing commercial solvers (such as Gurobi) can be used for calculation.

[0137] S11. Output the optimal solution, that is, the optimal attack vector of the attacker under the typical scenario ts.

[0138] S12. After the calculation of the first typical scenario according to steps S6 to S11, determine whether ts < k holds, that is, whether all typical scenarios have been calculated; if ts < k, it means there are still remaining typical scenarios to be calculated, and it is necessary to set ts = ts + 1, and return to step S6, and repeat steps S6 to S11, and make a judgment after each typical scenario is calculated; when ts ≥ k, it means that all typical scenarios have been calculated, and the program ends.

[0139] The embodiment of the present invention also provides a computer-readable storage medium. The computer-readable storage medium stores computer-executable instructions. When the computer-executable instructions are called and executed by a processor, the computer-executable instructions cause the processor to implement the steps of the method described above.

[0140] In a specific embodiment, taking the modified IEEE-39 bus system as an example, the electrical connection diagram of the system is as Figure 2 shown. The numbers 1 to 39 in the figure represent the bus numbers, the letter G represents the installation location of the conventional generating units, and two wind farms are connected to buses 31 and 35 respectively. The effectiveness of the established model and the solution method is verified through simulation analysis.

[0141] In this example system, the line capacity limit is set to 1.2 times the normal line flow, and the detection threshold of the node load data is set to 50%. The cut-in wind speed of the wind turbine is taken as 2m / s, the cut-out wind speed is taken as 20m / s, and the rated wind speed is taken as 12m / s. The wind speed obeys the Weibull distribution with shape parameters and scale parameters of 4 and 10.35 respectively, and the load power obeys the normal distribution with expectation and variance of 1 and 0.1 respectively.

[0142] Using the hidden attack risk assessment method applied to the power system proposed in the embodiment of the present invention, the Monte Carlo method is first used to generate an initial scenario set of the uncertainty of the wind power output of nodes 31 and 35 and all load nodes, which includes 5000 scenarios in total. Then, the IKBTD method is used to screen out a typical wind power output and load power scenario set, which includes 59 typical scenarios in total. Then, the data of each typical scenario is input to construct a two-level programming model of the hidden FDIA, and then converted into a single-level mixed integer linear programming model. Finally, the model is calculated using the Gurobi solver to obtain the minimum attack cost and corresponding attack strategy for implementing the hidden FDIA in each typical scenario.

[0143] The lowest attack costs and corresponding attack strategies in all typical scenarios are counted, and typical scenarios with the same lowest attack costs and corresponding attack strategies are defined as the same type of attack scenarios. The top five attack scenarios with the lowest attack costs are shown in Table 1:

[0144] Table 1 Top five attack scenarios with the lowest attack costs considering wind power and load uncertainty

[0145]

[0146] If the uncertainty of wind power and load is not considered, the minimum attack cost and corresponding attack strategy under the original scenario are calculated as shown in Table 2:

[0147] Table 2 Optimal attack strategy under the original scenario without considering wind power and load uncertainty

[0148]

[0149] By comparing Table 1 and Table 2, it can be seen that after considering the uncertainty of wind power and load, the attacker can only adopt the attack vector with the same minimum attack cost as the original scenario with a certain probability, and the minimum attack costs in other scenarios are higher than the original scenario without considering the uncertainty of wind power and load, indicating that the original evaluation result is more conservative. The evaluation method provided according to the embodiment of the present invention can allocate defense resources more efficiently, which illustrates the effectiveness of the embodiment of the present invention.

[0150] The embodiment of the present invention first obtains the uncertainty data of wind power and load based on the improved scenario analysis method, and then establishes an upper-level attack model with the goal of minimizing the attacker's attack cost and a lower-level economic dispatch model with the goal of minimizing the power generation cost of the power grid. Based on the KKT condition, the constructed two-level programming model is converted into a single-level mixed integer programming model, and a commercial solver is called to solve it, so as to effectively identify the weak links in the operation of the power system.

[0151] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art of the present invention, several equivalent substitutions or obvious variations can be made without departing from the concept of the present invention, and the performance or use is the same, which should be regarded as belonging to the protection scope of the present invention.

Claims

1. A method for assessing the risk of covert attacks applied to power systems, characterized in that: It includes the following steps: S6. Input the wind power and load data under typical scenarios, and establish the objective function and constraints of the upper-layer attack model for minimizing the attack cost; S7. Establish the objective function and constraints of the lower-layer economic dispatch model for minimizing the power generation cost, and obtain the two-layer programming model of the stealthy FDIA; S8. Process the logical constraints in the upper-layer attack model and transform them into mathematical forms; S9. Process the lower-layer economic dispatch model; S10. Transform the two-layer programming model of the stealthy FDIA into a single-layer mixed-integer linear programming model, and use a solving software to calculate the single-layer mixed-integer linear programming model; S11. Output the optimal attack vector; Before step S6, it also includes: S1. Input the power grid, conventional thermal power unit, and deterministic load data; S2. Input the wind turbine unit and uncertainty load parameters, where the parameters include: the probability distribution of the wind turbine unit output and the probability distribution of the uncertainty load; S3. Generate an initial scenario set N and perform scenario clustering; S4. Screen out k typical wind power and load scenarios; S5. Define ts as the serial number of the reduced typical scenario, and let ts = 1, that is, first take the first typical scenario generated by reduction as the calculation object; After step S11, it also includes: S12. After the calculation of the first typical scenario according to steps S6 to S11, judge whether ts < k holds, that is, whether all typical scenarios have been calculated; if ts < k, it means that there are still remaining typical scenarios not calculated, and it is necessary to make ts = ts + 1 and return to step S6, and repeat steps S6 to S11, and make a judgment once after the calculation of each typical scenario; when ts ≥ k, it means that all typical scenarios have been calculated, and the program ends running.

2. The method for assessing the risk of covert attacks applied to a power system according to claim 1, characterized in that: The data in step S1 includes: branch reactance, upper and lower limits of the output of conventional thermal power units, power generation cost, and deterministic load node data.

3. The method for assessing the risk of covert attacks applied to a power system according to claim 1, characterized in that: The parameters in step S2 include: the probability distribution of the wind turbine unit output and the probability distribution of the uncertainty load.

4. The method for assessing the risk of covert attacks applied to a power system according to claim 1, characterized in that: Step S3 specifically includes: generating the initial scenario set N by the Monte Carlo method and performing scenario clustering based on the maximum minimum distance algorithm to determine the clustering center.

5. The method for assessing the risk of covert attacks applied to a power system according to claim 1, characterized in that: Step S4 specifically includes: screening out k typical wind power and load scenarios by using the k-means scenario reduction method improved based on typicality; the definition of the typicality TD(s) of scenario s is shown in formula (1); Where, F={f1,f2,…,f i ,…,f n }, represents the allied scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes in i is the i-th element in F, n is the number of ally scenes; E = {e1, e2, …, e i ,…,e m }, represents the enemy scene set of scene s, that is, the subset c to which scene s belongs u The set of other scenes other than i is the i-th element in E, m is the number of adversary scenarios; IS(s) and ED(s) represent the internal similarity and external similarity calculated by the average Euclidean distance between scenario s and the ally scenario set F and the adversary scenario set E, respectively; Z ru , Z d is the normalization factor.

6. The method for assessing the risk of covert attacks applied to a power system according to claim 5, characterized in that: The normalization factor Z ru , Z d The specific calculation method is shown in formula (2): In the formula, s i and j are any two scenes in the initial scene set N; s a and b is a subset c u k is the preset number of scene reductions; u is the current scene clustering center.

7. The method for assessing covert attack risks applied to a power system according to claim 1, characterized in that: The calculation method of the objective function of the upper-layer attack model in step S6 is shown in formulas (3) to (5); Where N PD 、N PL are the node set and line set of the power system respectively; a i 、a j is a 0-1 integer variable, representing the attack strategy for node i and line j, and the false injection power ΔD of the node and line respectively. i , ΔPL j There is the following logical relationship: Where, ΔD i , ΔPL j are the false power data injected by the attacker into node i and line j respectively.

8. The method for assessing the risk of covert attacks applied to a power system according to claim 7, characterized in that: The constraints of the upper-layer attack model in step S6 include: network security constraints, stealthy attack constraints, and attack target constraints.

9. The method for assessing the risk of covert attacks applied to a power system according to claim 8, characterized in that: The calculation method of the network security constraints is shown in formulas (6) to (9); PL T =H·PG–H·(PD-S) (6) PL F =H·PG-H·(PD+ΔD-S) (7) ΔPL=PL T -PL F =H·ΔD (8) Where PL T The real data of the line flow, that is, the information that can be obtained from the attacker's perspective; PL F represents the false data of the line flow, that is, the information that can be obtained from both the attacker and the defender's perspective; ΔPL and ΔD represent the false data injected by the attacker into the line flow and node power; PG and PD represent the original generator output and load power of all nodes respectively; S represents the load power removed by the system; H is the conversion matrix, which represents the linear relationship between the line flow and the node power. The calculation formula is as follows: In the formula, B M and B S Respectively represent the mutual admittance matrix and the node self-admittance matrix; A A represents the node-branch association matrix; The calculation method of the stealthy attack constraints is shown in formulas (10) to (13); -(1-th j )·PL j.lim ≤PL j.F ≤(1-n) j )·PL j.lim (11) -τ·PD i ≤ΔD i ≤τ·PD i (12) In the formula, S i represents the load power removed at node i; η j ∈[0,1), indicating the margin left by the dispatcher for the transmission capacity of the line; PL j.F and PL j.lim are the false power flow value and capacity limit value of line j respectively; PD i represents the original load power at node i; τ∈(0,+∞), represents the detection threshold of the defender for false data of the node, which depends on the accuracy of the detection equipment invested in the power grid; The calculation method of the attack target constraints is shown in formula (14); Where PL j.T is the real power flow value of line j.

10. The method for assessing covert attack risks applied to a power system according to claim 7, characterized in that: The calculation method of the objective function of the lower-layer economic dispatch model in step S7 is shown in formula (15); min∑c g ·PG g (15) In the formula, c g Represents the unit output cost of each generator set; PG g The output of each generator set.

11. The method for assessing covert attack risks applied to a power system according to claim 10, characterized in that: The calculation method of the constraint conditions of the lower-level economic dispatch model in step S7 is shown in formulas (16) to (18): 0≤S d ≤D d +ΔD d (18) In the formula, represents the capacity limit of line l; Respectively represent the minimum output and maximum output limit of generator g; D d Represents the actual load data of node d; S d Represents the load shedding power of node d.

12. The method for assessing covert attack risks applied to a power system according to claim 9, characterized in that: Step S8 specifically includes: for the logical constraint represented by formula (14), first define a 0-1 variable δ for each line j , as a sign of whether the real power flow of the line exceeds the line capacity limit, that is: To achieve the purpose of the attack, the attacker needs to inject false data, and the line overload constraint that should be satisfied is: Using the Big M method, the logical constraints of formula (14) are transformed into the following mixed integer form: PL j.lim -|PL j.T |≤M·δ j -e (21) In the formula, M represents a sufficiently large positive number, and ε represents a sufficiently small positive number; Introduce a 0-1 integer variable δ j.ori By removing the absolute value sign in formula (21), the calculation formula is: At this point, the logical constraints of formula (14) are transformed into the form of formula (20) and formula (22).

13. The method for assessing covert attack risks applied to a power system according to claim 12, characterized in that: Step S8 also includes: converting formula (4) and formula (5) into the form of formula (23) and formula (24), and the calculation formula is as follows: In the formula, δ D+,d , δ D-,d , δ PL+,j and δ PL-,j is an additional 0-1 integer variable, which is related to the positive and negative direction of the injected power; M represents a sufficiently large positive number, ε represents a sufficiently small positive number; τ∈(0,+∞) represents the detection threshold of the defender for false node data, which depends on the accuracy of the detection equipment invested in the power grid.

14. The method for assessing covert attack risks applied to a power system according to claim 11, characterized in that: Step S9 specifically includes: using the Kuhn-Tucker KKT method to transform the lower-level economic dispatch model into the form of formula (25) to formula (28): Where λ, μ are Lagrange multipliers introduced to satisfy the equality constraints; H g T , H d T are the transpose of the g-th and d-th column vectors of matrix H respectively; cs d represents the load shedding loss at node d; PG g The output of each generator set; Respectively represent the minimum output and maximum output limit of generator g; S d represents the load shedding power of node d; α l , β g , γ d , Lagrange multiplier introduced to satisfy the upper and lower limit constraints of variables; ω α ,l , ω β ，g , ω γ ,d , They are 0-1 integer variables additionally introduced into the complementary relaxation conditions corresponding to the line power, generator output and node load shedding constraints.

15. The method for assessing covert attack risks applied to a power system according to claim 14, characterized in that: Step S10 specifically includes: through steps S8 and S9, the two-level programming model of the secret FDIA is converted into a single-level mixed integer linear programming model with formula (3) as the objective function and formulas (6) to (13), (20), (22) to (28) as constraints, and the single-level mixed integer linear programming model is calculated using a solving software.

16. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer-executable instructions, and when the computer-executable instructions are called and executed by a processor, the computer-executable instructions prompt the processor to implement the steps of the method according to any one of claims 1 to 15.

Images