Quantitative Assessment Method for Wide Area Control System Cyber ​​Attacks Based on Reachability Analysis

Abstract

The present invention discloses a quantitative assessment method for network attacks on wide-area control systems based on reachability analysis, comprising: establishing a nonlinear differential equation model of a wide-area control system with wind farm access based on a simplified admittance matrix; using Zeno polyhedrons to model the network attack amplitude range and the wind farm output variation range, respectively, and adding the established models to the nonlinear differential equation model of the wide-area control system; based on reachability analysis theory, improving the reachable set calculation algorithm using a Krylov subspace iteration method based on a stable biconjugate gradient method and zonotope bundle technology; applying the improved reachable set calculation algorithm to calculate the reachable set of the wide-area control system under the influence of different types of network attacks and wind power disturbances. The present invention quantitatively evaluates the damping performance of the wide-area control system under different network attack scenarios and operating conditions based on the reachable set results, providing important guidance for improving the damping performance of the wide-area controller and ensuring the safe and stable operation of the power system.

Description

Technical Field

[0001] The present invention relates to the technical field of power system security assessment, and in particular to a quantitative assessment method for wide-area control system network attacks based on reachability analysis. Background Art

[0002] As the proportion of renewable energy continues to increase, the uncertainty and volatility brought by renewable energy sources have exacerbated the complexity of inter-regional power system operations, further exacerbating inter-regional low-frequency oscillations. To effectively suppress inter-regional low-frequency oscillations, wide-area damping controllers (WACS) based on wide-area measurement systems have been widely used in wide-area control systems. These controllers use widely available wide-area signals, such as synchronous generator angular frequency deviation, relative power angle, and tie-line power, as inputs to detect oscillations in the regional system in real time. By adjusting the active power output of the synchronous generators, they provide additional damping to suppress the magnitude and duration of inter-regional low-frequency oscillations. However, during the information exchange between the wide-area control system through data measurement, network communication, and monitoring and control systems, hackers can launch cyberattacks against physical devices such as phasor measurement units and communication devices in the WACS. These attacks can intrude and tamper with data collected and transmitted, thereby compromising the integrity, availability, and authenticity of measurement or control signals. Such cyberattacks can not only weaken the damping performance of the WACS but also negatively impact regional oscillations and pose a serious threat to the safe and stable operation of the power system. Therefore, in order to ensure the effectiveness of the wide-area damping controller and maintain the safety and reliability of power system operation, it is crucial to study the quantitative analysis method of the impact of cyber attacks on the damping performance of the wide-area control system in the context of new energy access.

[0003] In existing research, domestic and international scholars have conducted extensive research on the implementation, detection, and mitigation of cyberattacks to mitigate the impact of cyberattacks on wide-area control systems. However, current research rarely involves quantitative assessment of the impact of cyberattacks, making it difficult to directly determine the reliability of wide-area damping controllers when encountering different types of cyberattacks with variable amplitudes, which can easily affect the effectiveness of detecting and preventing cyberattacks. Furthermore, existing research primarily analyzes linearized models through limited time-domain simulations, making it difficult to capture regional oscillations under multiple uncertainties and disturbances, which can easily lead to overly optimistic analysis results. Therefore, there is a lack of a quantitative assessment method for the impact of cyberattacks on wide-area control systems that consider wind power fluctuations. Summary of the Invention

[0004] The purpose of the present invention is to overcome the shortcomings and deficiencies of the existing technology and provide a quantitative assessment method for network attacks on wide-area control systems based on reachability analysis. First, a simplified admittance matrix is ​​used to establish a wide-area control system model. While improving computational efficiency, the error accumulation problem caused by the interaction between differential and algebraic equations in the calculation of reachable sets is effectively solved. Secondly, the Krylov subspace iteration method based on the stable biconjugate gradient method and the zonotope bundle technology are used to improve the reachable set calculation algorithm, thereby improving the calculation accuracy and stability of the reachable set. Finally, the improved reachable set calculation algorithm is used to effectively predict all variation ranges of the wide-area control system state under various operating conditions, providing an important reference for determining the impact of various types of attack injections and the operating status of the wide-area control system under different working conditions.

[0005] To achieve the above objectives, the present invention provides a technical solution: a method for quantitatively assessing network attacks on a wide area control system based on reachability analysis, comprising the following steps:

[0006] Step 1: Based on the simplified admittance matrix, a nonlinear differential equation model of the wide-area control system with wind farm integration is established;

[0007] Step 2: Considering the variability of the cyber attack amplitude and the uncertainty of wind power disturbances, the cyber attack amplitude range and wind farm output variation range are modeled using Zeno polyhedra. These models are then incorporated into the nonlinear differential equation model of the wide-area control system.

[0008] Step 3: To effectively analyze the nonlinear differential equation model of the wide-area control system, an improved reachable set calculation algorithm is designed based on reachability analysis theory. The nonlinear differential equation model is converted into a linear differential equation model. The reachable set of the nonlinear differential equation model is obtained by combining the solution results of the linear differential equation model and the linearization error estimation results. The specific improvements to the reachable set calculation algorithm include: using the Krylov subspace iteration method based on the stable biconjugate gradient method to calculate the homogeneous solution of the linear differential equation model, and using the zonotope bundle technique to combine the Zeno polyhedra in the reachable set calculation process to improve the calculation accuracy of the reachable set.

[0009] Step 4: Apply the improved reachable set calculation algorithm to calculate the reachable sets of the wide-area control system under the influence of different types of network attacks and wind power disturbances. Based on the reachable set results, quantitatively analyze the impact of network attacks and wind power disturbances on the damping performance of the wide-area control system, thereby completing the quantitative evaluation of network attacks on the wide-area control system.

[0010] Furthermore, in step 1, the wide-area control system includes a synchronous generator, a wind farm, a wide-area damping controller, a transformer, and a power load; assuming that the wide-area control system has N nodes, g Synchronous generators, N l load nodes and N w There are N wind farms. p A wide-area damping controller is deployed in the synchronous generator of the wide-area control system. The wide-area damping controller uses the collectible wide-area signal to perceive the regional frequency oscillation in the wide-area control system in real time, and provides additional damping for the wide-area control system by adjusting the active power output of the synchronous generator, thereby suppressing the regional frequency oscillation.

[0011] Furthermore, in step 1, first, a nonlinear differential-algebraic equation model of the wide-area control system with wind farm access is constructed considering the impact of cyber attacks; secondly, based on the simplified admittance matrix, the node admittance matrix of the wide-area control system is reduced to an admittance matrix containing only the nodes of each generator; finally, the nonlinear differential-algebraic equation model of the wide-area control system is converted into a nonlinear differential equation model.

[0012] Furthermore, the nonlinear differential-algebraic equation model includes: a wind farm model, a synchronous generator model, a wide-area damping controller model, a network attack model, and a power system flow model;

[0013] The wind farm model is represented by a synchronous generator simulation model of a doubly-fed wind turbine generator. The synchronous generator simulation model of the i-th doubly-fed wind turbine generator is specifically represented as follows:

[0014]

[0015] Where i = 1, 2, ..., N w ; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the internal electromotive force of the doubly fed wind turbine generator; θ ei 、ω ri 、E si are the virtual rotor angle, rotor speed and internal electromotive force of the doubly fed wind turbine generator respectively; is the reference value of the rotor speed of the doubly fed wind turbine generator; I ri0 、i qri They represent the initial value of the rotor current of the doubly fed wind turbine generator and the q-axis component of the rotor current respectively; K I1i With K p1i are the integral gain and proportional gain of the PI controller in the rotor-side converter of the doubly fed wind turbine generator; Hi With L mi are the inertia constant and mutual inductance parameters of the doubly fed wind turbine generator respectively; the reactance of the doubly fed wind turbine generator is x si =ω s L si , L si is the stator inductance of the doubly fed wind turbine generator, ω s is the rated value of the rotor angle of the doubly fed wind turbine generator; V pi represents the voltage of the common coupling node connected to the doubly fed wind turbine generator, θ pi is the common coupling node voltage V pi Phase angle; P wmi is the mechanical power of the doubly fed wind turbine generator;

[0016] The synchronous generator model adopts the classic second-order model, which is specifically expressed as follows:

[0017]

[0018] Where j = 1, 2, ..., N g , is the derivative value of the power angle of the jth synchronous generator, is the derivative value of the angular velocity of the jth synchronous generator, ω j represents the angular velocity of the jth synchronous generator; D j 、M j are the damping coefficient and moment of inertia of the j-th synchronous generator respectively; P mj is the mechanical power of the jth synchronous generator, P ej is the electromagnetic power output by the jth synchronous generator, specifically expressed as Among them E j is the internal potential of the synchronous generator, V j is the node voltage connected to the jth synchronous generator, δ j represents the power angle of the jth synchronous generator, θ j is the node voltage V j The phase angle of G gj With B gj are the conductance and susceptance values ​​of the line between the synchronous generator and the connected node, respectively;

[0019] The wide-area damping controller model uses the difference in angular velocity of synchronous generators located in two different areas as a feedback signal and controls the output of the wide-area damping controller model through filtering, phase compensation, gain, and saturation adjustment. The specific expression of the wide-area damping controller model is:

[0020]

[0021] In the formula, the coefficient matrix Coefficient matrix B wm =[0;1], coefficient matrix m=1,2,...,N p ; is the derivative value of the state variable vector of the wide-area damping controller model, P WADCm represents the output value of the wide-area damping controller model, x wm is the state variable vector of the wide-area damping controller model, K m 、T m and T wm are the gain parameter of the wide-area damping controller model, the time constant of the phase compensation component, and the time constant of the filter; Δω ab is the angular velocity difference between the a-th and b-th synchronous generators, which is set as the control signal of the wide-area damping controller model. The output of the wide-area damping controller model can be used as part of the active power output command of the synchronous generator, enabling the synchronous generator to respond to the oscillation of the regional frequency and enhance the damping effect of the wide-area control system.

[0022] The network attack model considers the network attacks that the wide-area control system may suffer during the data acquisition and transmission process, namely, false data injection attacks and delay attacks, wherein the false data injection attacks include disturbance-type false data injection attacks and proportional-type false data injection attacks; when the control signal Δω of the wide-area damping controller model is ab After being affected by the injected network attack, Δω ab The control signal of the wide-area damping controller model after being affected by the attack is converted into

[0023] The specific representation of the perturbation-type false data injection attack is as follows:

[0024]

[0025] Where t is the running time of the wide area control system, d is the injection component of the disturbance type false data injection attack, and t c Indicates the injection time of the perturbation-type false data injection attack;

[0026] The specific representation of the proportional false data injection attack is as follows:

[0027]

[0028] Where, α sc is the proportional gain of proportional false data injection attack, t sc Indicates the injection time of proportional false data injection attack;

[0029] For delayed attacks, the second-order modified Padé approximation is used to model the time delay, and the time delay model is obtained. The specific state space equation is expressed as: in, is the derivative value of the state variable vector of the time delay model, the coefficient matrix Coefficient matrix Coefficient matrix C d =[0,1],x d Represents the state variable vector of the time delay model, T d Indicates the delay time. represents the angular velocity difference of the synchronous generator after being affected by the time delay; the control signal of the wide-area damping controller model after being affected by the time delay attack is expressed as:

[0030]

[0031] Where, t d Indicates the injection time of the delayed attack;

[0032] When the network attack component directly contaminates the control signal of the wide-area damping controller model, the wide-area damping controller model will become:

[0033]

[0034] Where, is the state variable vector of the wide-area damping controller model after being affected by the network attack, for The derivative value of is the abnormal output value of the wide-area damping controller model after being affected by the cyber attack;

[0035] The power system flow model is specifically expressed as:

[0036]

[0037] Where V k is the voltage of the kth node of the wide area control system, V h is the voltage of the hth node of the wide-area control system, P k is the active power flowing into the kth node, P ek is the active power output by the synchronous generator at the kth node, P lk is the load power at the kth node, P wk is the output power of the wind farm at the kth node, θ k is the phase angle of the voltage at the kth node, θ h is the phase angle of the voltage at the hth node, B kh With G khare the conductance and susceptance of the line between the kth node and the hth node respectively;

[0038] The above-mentioned wind farm model, synchronous generator model, wide-area damping controller model, network attack model, and power system flow model are transformed into a nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system. At the same time, in order to record the time changes during the calculation process, an extended variable t representing the operating time of the wide-area damping control system is introduced. Its state equation is expressed as follows: When no network attack is added, define the state variable vector matrix Algebraic variables vector matrix Among them, δ p 、ω p ,θ ep 、ω rp 、E sp 、x wp denote the power angle of the synchronous generator at the p-th node, the angular velocity of the synchronous generator, the virtual rotor angle of the doubly fed wind turbine, the rotor speed of the doubly fed wind turbine, the internal electromotive force of the doubly fed wind turbine, and the state variable vector of the wide-area damping controller model, respectively. is the field of real numbers, E p 、V p ,θ p are the internal potential of the synchronous generator at the pth node, the node voltage amplitude and the node voltage phase angle respectively; the wide area control system can be modeled as:

[0039]

[0040] Where x(t), y(t) and u(t) are the state variable vector, algebraic variable vector and input variable set of the wide area control system under normal operating conditions, respectively. is the derivative value of the state variable vector of the wide-area control system under normal operating conditions, f(·) and g(·) represent the set of nonlinear differential equations and algebraic equations describing the system dynamic characteristics of the wide-area control system, respectively;

[0041] After adding the network attack, the nonlinear differential-algebraic equation model of the wide-area control system is transformed into:

[0042]

[0043] Where x c (t), y c (t) and u c (t) are the state variable vector, algebraic variable vector and input variable set of the wide area control system after being affected by the network attack, is the derivative value of the state variable vector of the wide-area control system after being affected by the network attack;

[0044] In order to solve the error accumulation problem caused by the interaction between differential and algebraic equations in the calculation of reachable sets, a simplified admittance matrix construction method is adopted to achieve the reduction of the node admittance matrix of the wide-area control system through matrix operation. Based on the obtained reduced admittance matrix, the electromagnetic power output by the j-th synchronous generator is expressed as:

[0045]

[0046] Where G jj is the self-conductance of the j-th synchronous generator, E j With E q are the internal potentials of the jth and qth synchronous generators, δ j and δ q are the power angles of the jth and qth synchronous generators respectively, B jq With G jq are the conductance and susceptance between the jth and qth synchronous generators, respectively. According to the above formula (11), it can be seen that the power distribution in the wide-area control system can be obtained without calculating the algebraic equations;

[0047] In addition, considering the access of wind farms, it is necessary to retain the internal nodes of the wind farm and the common coupling nodes connected to the wind farm. Assuming that the common coupling node of the wind farm is the rth node in the wide-area control system, according to Kirchhoff's current law, we can obtain:

[0048]

[0049] Where G wp With B wp are the conductance and susceptance between the internal nodes of the wind farm and the common coupling node, G p0 With B p0 are the ground conductance and susceptance values ​​of the common coupling node, G gr With B gr are the conductance and susceptance between the common coupling node and the external g-th node, V g and δ g are the voltage and phase angle of the g-th node respectively; According to the above formula (12), the internal potential E of the wind farm can be si The voltage of the common coupling node is obtained, which is specifically expressed as:

[0050]

[0051] Where B qq is the self-susceptance of the common coupling node of the wind farm; the electromagnetic power output by the wind farm obtained by the above formula (13) can be used to establish the coupling relationship between the wind farm and other synchronous generators;

[0052] Therefore, the nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system can be transformed into a nonlinear differential equation model in the form of a differential equation system, which can be specifically expressed as:

[0053]

[0054] According to the above formula (14), without explicitly processing the algebraic equations, the dynamics of the wide-area control system can be obtained directly by solving the constructed nonlinear differential model.

[0055] Furthermore, in step 2, considering the variability of the network attack amplitude and the impact of uncertain wind power disturbances caused by random changes in wind speed, Zeno polyhedrons are used to model the network attack amplitude range and the wind farm output variation range respectively, and the constructed models are added as uncertain input sets to the nonlinear differential equation model of the wide-area control system, thereby reflecting the impact of network attack amplitude changes and wind power disturbances on the damping performance of the wide-area control system.

[0056] Furthermore, based on the construction method of Zeno polyhedron, the model of the network attack amplitude range is expressed as:

[0057]

[0058] Where,

[0059] Indicates the range of changes in the magnitude of network attacks; represents the amplitude center of the network attack amplitude, l is a variable used to record the change in the number of generators of the Zeno polyhedron during the modeling process, is the range of network attack amplitude variation represented by the generator of Zeno polyhedron, where α l is the proportionality coefficient, g l is the generator vector; w1 is the number of generators of the Zeno polyhedron used to model the range of network attack amplitude variations; by setting the Zeno polyhedron, a model of network attack amplitude range of any size can be generated, and this model is added to the network attack model to quantitatively analyze the impact of variable attacks;

[0060] Considering the impact of wind power fluctuations on wide-area control systems, the Zeno polyhedron is used to model the wind farm output variation range and integrated into the wind farm model. The model of the wind farm output variation range is specifically expressed as:

[0061]

[0062] Where, The output variation range of the wind farm; is the center of change of wind farm output, set to 0; is the fluctuation range of wind farm output represented by the generator of Zeno polyhedron, where β l is the proportional coefficient, h l is the generator vector; w2 is the number of generators of the Zeno polyhedron used to model the output variation range of the wind farm.

[0063] Furthermore, in step 3, the reachability analysis theory obtains the reachable set of the wide-area control system by conservatively predicting the dynamic operation trajectory of the wide-area control system under different operating conditions. The reachable set can effectively analyze the operating characteristics of the wide-area control system under random and uncertain working conditions.

[0064] When performing reachability analysis on a wide-area control system with nonlinear characteristics, it is necessary to first approximate the nonlinear differential equation model into a linear differential equation model through the method of linear abstraction of the state space, and then use the reachable set calculation algorithm applicable to the linear differential equation model to perform calculations; for the reachable set calculation algorithm of the linear differential equation model, the Krylov subspace iteration method based on the stable biconjugate gradient method is used to improve the reachable set calculation algorithm of the linear differential equation model. Specifically, the Krylov subspace is defined as in, For vector The vector space generated by the linear combination of , span(·) represents the linear span algorithm that returns a set of vectors, n represents the dimension of the Krylov subspace, The variable vector represents the linear equation system in the linear differential equation model, and A is the coefficient matrix of the linear equation system. In order to obtain the homogeneous solution of the linear differential equation system, the problem of solving the linear differential equation system is transformed into the problem of solving the best approximation in the Krylov subspace. In order to improve the stability of the calculation, the stable biconjugate gradient method is used instead of the Arnoldi method to solve the problem, thereby obtaining the solution of the linear differential equation model. In addition, in order to compensate for the linearization error caused by the linear abstraction method of the state space, the over-approximation method is used to determine the linearization calculation error, and the linearization calculation error is added to the reachable set calculation result of the linear differential equation model in the form of an uncertain input set, and the reachable set calculation result of the nonlinear system is obtained within the error tolerance range.

[0065] In order to improve the processing efficiency and calculation accuracy of Zeno polyhedra during the reachable set calculation process, the zonotope bundle technology is introduced to combine Zeno polyhedra during the reachable set calculation process. The zonotope bundle technology is specifically expressed as follows: Where o represents the number of Zeno polyhedra, is a reachable set described in terms of a Zeno polyhedron, is the reachable set after combining Zeno polyhedra using zonotope bundle technology, M is the total number of Zeno polyhedra; after calculating the reachable set results at each moment, the sth calculation step Δt s As a unit, apply the convex hull calculation CH(·) to obtain the reachable set within the specified calculation step size The calculation process is expressed as: in t s time and t s+1 The reachable set of time; Finally, in order to obtain the specified time range [t0,t f ] reachable set within The reachable sets computed at each individual time step need to be combined: where t0, t f The start and end times are set respectively.

[0066] Furthermore, in order to consider the reachable sets under various operating conditions of the wide-area control system, a block calculation method is adopted to calculate the reachable sets of different operating states in sequence, and the reachable set result of the previous operating state is used as the initial set of the next operating state. Finally, the union of the reachable sets of the four operating states consisting of normal operating state, attack state, fault state, and fault removal state is taken as the final calculation result, which is specifically expressed as: in, and Represent the reachable sets under four operating states respectively, is the final reachable set result, t0, t1, t2, t3, and t4 are the start and end times of the four operating states, respectively. The obtained reachable set result can reflect the frequency oscillation between regions of the wide-area control system under the influence of the set network attack and wind power disturbance.

[0067] Furthermore, in step 4, based on the established nonlinear differential equation model of the wide-area control system, the improved reachable set calculation algorithm in step 3 is used to calculate the reachable sets of the wide-area control system affected by different types of network attacks and wind power disturbances. The obtained reachable set results are used to quantitatively analyze the impact of the set network attacks and wind power disturbances on the damping performance of the wide-area control system. At the same time, in order to further quantify the impact of uncertainties, two performance indicators are calculated to respectively reflect the maximum deviation degree and the maximum fluctuation range of the inter-regional frequency difference of the wide-area control system. The specific conditions of these two performance indicators are as follows:

[0068] Indicator-JD : Where: e(t) is Δω ab The difference between the instantaneous value and the steady-state value, t1 and t2 are the starting time and the ending time of the indicator calculation respectively; this indicator is used to capture Δω ab The maximum deviation, J D The higher the value of , the worse the damping performance of the wide-area control system;

[0069] Indicator 2J F : Where: Δω abmax and Δω abmin Represent Δω respectively ab By calculating the maximum and minimum values ​​of this index, we can effectively reflect the degree of influence of various uncertain factors on the wide area control system. F The larger the value, the greater the impact of fluctuations caused by uncertain factors.

[0070] Furthermore, by combining the reachable set results with the two performance index results, the inter-regional frequency oscillations of the wide-area control system after the injection of different types of network attacks are compared, so that the type of network attack that has the greatest impact on the damping performance of the wide-area control system can be obtained; in addition, by adding the network attack amplitude range model to the network attack model, and adding the wind farm output variation range model to the wind farm model, the inter-regional frequency oscillations of the wide-area control system under the dual influence of network attack and wind power disturbance are calculated, reflecting the additional impact of the uncertainty of large-scale renewable energy on the effect of network attack; based on the reachable set results and the two performance index results, the vulnerability of the wide-area control system under different attack scenarios and operating conditions can be evaluated, providing important guidance for improving the damping performance of the wide-area control system and ensuring the safe and stable operation of the power system.

[0071] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0072] 1. Based on a simplified admittance matrix, the present invention establishes a nonlinear differential model of a wide-area control system with wind farm access. This model can calculate the active power distribution of the wide-area control system without constructing algebraic equations. While improving computational efficiency, it effectively solves the problem of error accumulation in reachable set calculations caused by the interaction between differentials and algebraic equations.

[0073] 2. The present invention takes into account the variability of the attack amplitude and the uncertainty of wind power disturbance, and uses Zeno polyhedron to model the attack amplitude and wind power disturbance variation range. The established model of network attack amplitude range and wind farm output variation range is added as uncertain input variables to the nonlinear differential model of the wide-area control system, providing a convenient and effective modeling method for quantitatively considering the impact of the dual effects of attack amplitudes and wind power disturbances with different variation ranges on the inter-regional oscillation of the wide-area control system.

[0074] 3. The present invention improves the reachable set calculation algorithm by using the Krylov subspace iteration method based on the stable biconjugate gradient method and the zonotopebundle technology, thereby improving the calculation accuracy of the reachable set and the stability of the reachable set calculation algorithm.

[0075] 4. This paper utilizes a reachability analysis method to effectively predict the full range of operating state variations in the wide-area control system under various operating conditions, exploring the maximum fluctuation range of regional frequency oscillations in the wide-area control system under different types of cyberattacks and wind power disturbances. Furthermore, two performance indicators are introduced to record the maximum deviation of frequency differences between different regions of the wide-area control system and the cumulative fluctuation range, respectively. This provides an effective tool for determining the operating state of the wide-area control system under various attack injection types and different operating conditions, and offers important guidance for improving the performance of the wide-area damping controller and the safety and reliability of system operation.

[0076] 5. The model studied in this invention retains its original nonlinear characteristics. The reachability analysis method adopted does not require manual simplification of the nonlinear part of the model and can effectively handle the nonlinearity, randomness and uncertainty in the actual system. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 Flowchart of the method of the present invention.

[0078] Figure 2 This is a structural diagram of a wide area control system according to an embodiment of the present invention.

[0079] Figure 3 This is a flow chart of reachability analysis of a wide area control system according to an embodiment of the present invention. DETAILED DESCRIPTION

[0080] The present invention will be described in further detail below with reference to the embodiments and drawings, but the embodiments of the present invention are not limited thereto.

[0081] like Figure 1 As shown, this embodiment discloses a method for quantitatively assessing network attacks on a wide area control system based on reachability analysis, comprising the following steps:

[0082] Step 1: Based on the simplified admittance matrix, a nonlinear differential equation model of the wide-area control system with wind farm integration is established as follows:

[0083] First, a nonlinear differential-algebraic equation model of a wide-area control system with wind farm integration is constructed, taking into account the impact of cyber attacks. The specific model is as follows:

[0084] like Figure 2 As shown in Figure 1, the wide-area control system consists of wind farms, synchronous generators, wide-area damping controllers, cyber attacks, and power loads. The nonlinear differential-algebraic equation model includes the wind farm model, synchronous generator model, wide-area damping controller model, cyber attack model, and power system flow model. The specific models are:

[0085] The wind farm model is represented by the synchronous generator simulation model of the doubly fed wind turbine generator. The synchronous generator simulation model of the i-th doubly fed wind turbine generator is specifically expressed as follows:

[0086]

[0087] Where i = 1, 2, ..., N w ; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the internal electromotive force of the doubly fed wind turbine generator; θ ei 、ω ri 、E si are the virtual rotor angle, rotor speed and internal electromotive force of the doubly fed wind turbine generator respectively; is the reference value of the rotor speed of the doubly fed wind turbine generator; I ri0 、i qri They represent the initial value of the rotor current of the doubly fed wind turbine generator and the q-axis component of the rotor current respectively; K I1i With K p1i are the integral gain and proportional gain of the PI controller in the rotor-side converter of the doubly fed wind turbine generator; H i With L mi are the inertia constant and mutual inductance parameters of the doubly fed wind turbine generator respectively; the reactance of the doubly fed wind turbine generator is x si =ω s L si , L si is the stator inductance of the doubly fed wind turbine generator, ω s is the rated value of the rotor angle of the doubly fed wind turbine generator; V pi represents the voltage of the common coupling node connected to the doubly fed wind turbine generator, θ pi is the common coupling node voltage V piPhase angle; P wmi is the mechanical power of the doubly fed wind turbine generator;

[0088] The synchronous generator model adopts the classic second-order model, which is specifically expressed as:

[0089]

[0090] Where j = 1, 2, ..., N g , is the derivative value of the power angle of the jth synchronous generator, is the derivative value of the angular velocity of the jth synchronous generator, ω j represents the angular velocity of the jth synchronous generator; D j 、M j are the damping coefficient and moment of inertia of the j-th synchronous generator respectively; P mj is the mechanical power of the jth synchronous generator, P ej is the electromagnetic power output by the jth synchronous generator, specifically expressed as Among them E j is the internal potential of the synchronous generator, V j is the node voltage connected to the jth synchronous generator, δ j represents the power angle of the jth synchronous generator, θ j is the node voltage V j The phase angle of G gj With B gj are the conductance and susceptance values ​​of the line between the synchronous generator and the connected node, respectively;

[0091] The wide-area damping controller model uses the difference in angular velocity of synchronous generators located in two different areas as the feedback signal. It controls the output of the wide-area damping controller model through filtering, phase compensation, gain, and saturation adjustment. The specific expression of the wide-area damping controller model is:

[0092]

[0093] In the formula, the coefficient matrix Coefficient matrix B wm =[0;1], coefficient matrix m=1,2,...,N p ; is the derivative value of the state variable vector of the wide-area damping controller model, P WADCm represents the output value of the wide-area damping controller model, x wm is the state variable vector of the wide-area damping controller model, K m 、T m and T wmare the gain parameter of the wide-area damping controller model, the time constant of the phase compensation component, and the time constant of the filter; Δω ab is the angular velocity difference between the a-th and b-th synchronous generators, which is set as the control signal of the wide-area damping controller model; the output of the wide-area damping controller model can be used as a part of the active power output instruction of the synchronous generator, so that the synchronous generator can respond to the oscillation of the regional frequency and enhance the damping effect of the wide-area control system;.

[0094] The network attack model mainly considers the network attacks that the wide-area control system may suffer during the data acquisition and transmission process, namely, false data injection attacks and delay attacks. Among them, false data injection attacks include disturbance-type false data injection attacks and proportional-type false data injection attacks. When the control signal Δω of the wide-area damping controller model is ab After being affected by the injected network attack, Δω ab The control signal of the wide-area damping controller model after being affected by the attack is converted into

[0095] The specific representation of the perturbation-type false data injection attack is as follows:

[0096]

[0097] Where t is the running time of the wide area control system, d is the injection component of the disturbance type false data injection attack, and t c Indicates the injection time of the perturbation-type false data injection attack;

[0098] The specific representation of the proportional false data injection attack is as follows:

[0099]

[0100] Where, α sc is the proportional gain of proportional false data injection attack, t sc Indicates the injection time of proportional false data injection attack;

[0101] For delayed attacks, the second-order modified Padé approximation is used to model the time delay, and the time delay model is obtained. The specific state space equation is expressed as: in, is the derivative value of the state variable vector of the time delay model, the coefficient matrix Coefficient matrix Coefficient matrix C d =[0,1]; Represents the state variable vector of the time delay model, T d Indicates the delay time. represents the angular velocity difference of the synchronous generator after being affected by the time delay; the control signal of the wide-area damping controller model after being affected by the time delay attack is expressed as:

[0102]

[0103] Where, t d Indicates the injection time of the delayed attack;

[0104] When the network attack component directly contaminates the control signal of the wide-area damping controller model, the wide-area damping controller model will become:

[0105]

[0106] Where, is the state variable vector of the wide-area damping controller model after being affected by the network attack, for The derivative value of is the abnormal output value of the wide-area damping controller model after being affected by the cyber attack;

[0107] The power system flow model is specifically expressed as:

[0108]

[0109] Where V k is the voltage of the kth node of the wide area control system, V h is the voltage of the hth node of the wide-area control system, P k is the active power flowing into the kth node, P ek is the active power output by the synchronous generator at the kth node, P lk is the load power at the kth node, P wk is the output power of the wind farm at the kth node, θ k is the phase angle of the voltage at the kth node, θ h is the phase angle of the voltage at the hth node, B kh With G kh are the conductance and susceptance of the line between the kth node and the hth node respectively;

[0110] The above-mentioned wind farm model, synchronous generator model, wide-area damping controller model, network attack model, and power system flow model are transformed into a nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system. At the same time, in order to record the time changes during the calculation process, an extended variable t representing the operating time of the wide-area damping control system is introduced. Its state equation is expressed as follows: When no network attack is added, define the state variable vector matrix Algebraic variables vector matrix Among them, δ p 、ω p ,θ ep 、ω rp 、E sp 、x wp denote the power angle of the synchronous generator at the p-th node, the angular velocity of the synchronous generator, the virtual rotor angle of the doubly fed wind turbine, the rotor speed of the doubly fed wind turbine, the internal electromotive force of the doubly fed wind turbine, and the state variable vector of the wide-area damping controller model, respectively. is the field of real numbers, E p 、V p ,θ p are the internal potential of the synchronous generator at the pth node, the node voltage amplitude and the node voltage phase angle respectively; the wide area control system can be modeled as:

[0111]

[0112] Where x(t), y(t) and u(t) are the state variable vector, algebraic variable vector and input variable set of the wide area control system under normal operating conditions, respectively. is the derivative value of the state variable vector of the wide-area control system under normal operating conditions, f(·) and g(·) represent the set of nonlinear differential equations and algebraic equations describing the system dynamic characteristics of the wide-area control system, respectively;

[0113] After adding the network attack, the nonlinear differential-algebraic equation model of the wide-area control system is transformed into:

[0114]

[0115] Where x c (t), y c (t) and u c (t) are the system state variable vector, algebraic variable vector and input variable set after being affected by the network attack, respectively, x c (t) is the derivative value of the state variable vector of the wide area control system after being affected by the network attack;

[0116] Furthermore, based on the simplified admittance matrix construction theory, the node admittance matrix of the wide-area control system is reduced to a matrix containing only the nodes of each synchronous generator, thereby transforming the differential-algebraic model of the wide-area control system into a nonlinear differential model, specifically:

[0117] First, all loads in the wide-area control system are converted into equivalent admittances. The internal reactance of the synchronous generator is added to the impedance of the transformer. The internal nodes of the synchronous generator are isolated to obtain the N-order node admittance matrix Y of the wide-area control system and divide it into two parts. Specifically, it is expressed as: Where n represents the node connected to the synchronous generator, represents the rest of the nodes in the network, Y nn is the self-admittance of the synchronous generator at node n, and is the mutual admittance between the synchronous generator node n and the remaining nodes n, For the remaining nodes Furthermore, the admittance matrix is ​​reduced by matrix operation, and the reduced node admittance matrix is ​​expressed as Finally, based on the obtained reduced admittance matrix, the electromagnetic power output of the j-th synchronous generator can be expressed as

[0118]

[0119] Among them, G jj is the self-conductance of the j-th synchronous generator, E j With E q are the internal potentials of the jth and qth synchronous generators, δ j and δ q are the power angles of the jth and qth synchronous generators respectively, B jq With G jq are the conductance and susceptance between the jth and qth synchronous generators, respectively. According to the above formula (11), it can be seen that the power distribution in the wide-area control system can be obtained without calculating the algebraic equations;

[0120] Considering the integration of wind farms, in order to retain the internal nodes of the wind farm and the common coupling nodes connected to the wind farm, assuming that the common coupling node is the rth node in the wide-area control system, according to Kirchhoff's current law, we can obtain:

[0121]

[0122] Where G wp With B wp are the conductance and susceptance between the internal nodes of the wind farm and the common coupling node, G p0 With B p0 are the ground conductance and susceptance values ​​of the common coupling node, G gr With B gr are the conductance and susceptance between the common coupling node and the external g-th node, V g and δ g are the voltage and phase angle of the g-th node respectively; According to the above formula (12), the internal potential E of the wind farm can be si The voltage of the common coupling node is obtained, which is specifically expressed as:

[0123]

[0124] Where B qq is the self-susceptance of the common coupling node of the wind farm; the electromagnetic power output by the wind farm obtained by the above formula (13) can be used to establish the coupling relationship between the wind farm and other synchronous generators;

[0125] Therefore, the nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system can be transformed into a nonlinear differential equation model in the form of a differential equation system, which can be specifically expressed as:

[0126]

[0127] According to the above formula (14), without explicitly processing the algebraic equations, the dynamics of the wide-area control system can be obtained directly by solving the constructed nonlinear differential model.

[0128] Step 2: Considering the variability of the cyber attack amplitude and the impact of uncertain wind power disturbances caused by random changes in wind speed, the attack amplitude range and the wind farm output variation range are modeled using the regional isotope-based Zeno polyhedron construction method. The constructed models are then added to the wide-area control system model as uncertain input sets to reflect the impact of cyber attack amplitude changes and wind power disturbances on the damping performance of the wide-area control system. The specific implementation method is as follows:

[0129] Based on the construction method of Zeno polyhedron, the model of network attack amplitude range is expressed as:

[0130]

[0131] Where, Indicates the range of changes in the magnitude of network attacks; represents the amplitude center of the network attack amplitude, l is a variable used to record the change in the number of generators of the Zeno polyhedron during the modeling process, is the range of network attack amplitude variation represented by the generator of Zeno polyhedron, where α l is the proportionality coefficient, g l is the generator vector; w1 is the number of generators of the Zeno polyhedron used to model the variation range of network attack amplitude; by setting the Zeno polyhedron, a model of network attack amplitude range of any size can be generated, and this model is added to the network attack model to quantitatively analyze the impact of variable attacks.

[0132] Considering the impact of wind power fluctuations on wide-area control systems, the Zeno polyhedron is used to model the wind farm output variation range and integrated into the wind farm model. The model of the wind farm output variation range is specifically expressed as:

[0133]

[0134] Where, The output variation range of the wind farm; is the center of change of wind farm output, set to 0; is the fluctuation range of wind farm output represented by the generator of Zeno polyhedron, where β l is the proportional coefficient, h l is the generator vector; w2 is the number of generators of the Zeno polyhedron used to model the output variation range of the wind farm.

[0135] Step 3: Based on the reachability analysis theory, an improved reachable set calculation algorithm is designed to automatically transform the nonlinear differential equation model of the wide-area control system into a linear differential equation model. The reachable set of the nonlinear differential equation model is obtained by combining the solution results of the linear differential equation model and the linearization error estimation results. The details are as follows:

[0136] An improved reachable set calculation algorithm is designed based on reachability analysis theory. The specific improvements are as follows: based on the stable biconjugate gradient method, the Krylov subspace iteration method is used to calculate the homogeneous solution of the linear equation to improve the calculation efficiency and stability of the reachable set; at the same time, the zonotope bundle technology is used to combine a series of Zeno polyhedra in the reachable set calculation process to capture the running trajectory of the state variables of the wide-area control system in the dynamic process, thereby improving the calculation accuracy of the reachable set. Figure 3 As shown in Figure 2, the improved reachable set calculation process is as follows:

[0137] The reachability analysis of nonlinear differential equation models is mainly based on the state space linear abstraction method to calculate the reachable set, and the over-approximation method is used to determine the linearization calculation error. The error is added in the form of uncertain input, and the calculation result is compensated for the error, so as to obtain the reachable set of the nonlinear system within the error tolerance range. s is the starting time, t s+1 The specified time interval Δt for the end time s =[t s ,t s+1 ], calculate the reachable set of the nonlinear differential equation model The specific steps are as follows:

[0138] First, the Taylor expansion method is used to linearize the nonlinear differential equation model at the set linearization point. It is abstracted into a linear differential equation model: in: and They represent the state variable vector and the transposed value of the input variable set of the wide-area control system at the linearization point, respectively. * Indicated by and The variable vector set of the wide-area control system at the linearization point composed of x c (t) and u c (t) is the variable vector set of the wide area control system, represents the derivative value of the state variable vector of the wide-area control system, represents the Minkowski addition operator; is the Lagrange remainder, ξ is z * The value at any point in the domain. Through the approximate transformation, the nonlinear differential equation model is converted into the standard form of the linear differential equation model: in is x c (t) relative to The amount of change, for The derivative value of is the approximate coefficient matrix, Represents the disturbance input term containing approximate error and uncertain input. Through the above transformation, the reachable set calculation method of the linear differential equation model can be used to calculate the homogeneous solution of the linear differential equation model. In the traditional reachable set calculation algorithm, the linear differential equation model is s+1 Homogeneous solution of time The calculation method is: in, t s The reachable set at the moment, Δt is the specified calculation step, Represented by the center point u of the input set c Determine the reachable set. This calculation method directly substitutes into the matrix Exponential calculations can easily amplify matrix calculation errors and cannot handle situations where the matrix values ​​are large, which can easily lead to large errors in the calculation results and even non-convergence of the calculation. In order to deal with this problem, the Krylov subspace iteration method based on the stable biconjugate gradient method is used. Part. In the Krylov subspace, the problem of solving the linear equations is transformed into the problem of solving the best approximation problem. At the same time, the stable biconjugate gradient method is used instead of the Arnoldi method to solve the problem, thereby obtaining the solution of the linear model.

[0139] Next, calculate the error caused by linearization The calculated error result is used as part of the uncertain input to expand the reachable set to obtain t s+1 The reachable set at a moment The specific calculation formula is: in, The reachable set is determined by the uncertain input. To control the errors introduced by linearization and set calculation, the actual error needs to be calculated based on the obtained reachable set. If the actual error is within the predicted range, it means that the error meets the accuracy requirements. However, if the error exceeds the predicted value, the reachable set must be further modified by decomposing the reachable set or reducing the iteration step size.

[0140] Finally, the convex hull is used to calculate CH(·) to obtain the reachable set within the specified calculation step size. Specifically expressed as: In order to obtain the specified time range [t0,t f ] within the reachable set Where t0 is the calculation start time, t f To calculate the termination moment, the reachable sets computed at each individual time step need to be combined: Where s is the number of calculation intervals, ∪ represents the union calculation;

[0141] In order to improve the processing efficiency and calculation accuracy of Zeno polyhedra during the reachable set calculation process, the zonotope bundle technology is introduced to combine the Zeno polyhedra within each calculation step during the reachable set calculation process. The specific expression is: Where o represents the number of Zeno polyhedra, is a reachable set described in terms of a Zeno polyhedron, is the reachable set of Zeno polyhedra combined using zonotope bundling technology, and M is the total number of Zeno polyhedra.

[0142] To account for reachable sets under various system operating conditions, a block-by-block calculation method is used to sequentially calculate reachable sets for different operating states, using the reachable set result for the previous operating state as the initial set for the next operating state. Finally, the union of the reachable sets for the four operating states—normal operating state, attack state, fault state, and post-fault removal state—is taken as the final calculation result, specifically expressed as: in, and Represent the reachable sets under four operating states respectively, is the final reachable set result, t0, t1, t2, t3, and t4 are the starting and ending times of the four operating states respectively.

[0143] Step 4: Apply the improved reachable set calculation algorithm to calculate the reachable sets of the wide-area control system under the influence of different types of network attacks and wind power disturbances. Based on the reachable set results, quantitatively analyze the impact of network attacks and wind power disturbances on the damping performance of the wide-area control system. Specifically:

[0144] Based on the established nonlinear differential equation model of the wide-area control system, the improved reachable set calculation algorithm in step 3 is used to calculate the reachable set of the wide-area control system affected by different types of network attacks and wind power disturbances. That is, all possible operating trajectories of the wide-area control system state under certain disturbance conditions are solved, and the obtained reachable set results are used to quantitatively analyze the impact of the set network attacks and wind power disturbances on the damping performance of the wide-area control system. At the same time, in order to further quantify the impact of uncertain factors, two performance indicators are calculated to reflect the maximum offset degree and maximum fluctuation range of the frequency difference between regions of the wide-area control system. The indicators are specifically expressed as:

[0145] Indicator-J D : Where: e(t) is Δω ab The deviation value from its steady-state value, t1 and t2 are the start and end time of indicator calculation respectively. This indicator is used to capture Δω ab The maximum deviation, J D The higher the value of , the worse the damping performance of the wide-area control system.

[0146] Indicator 2J F : Where: Δω abmax and Δω abmin Represent Δω respectively ab By calculating this index, we can effectively reflect the impact of various uncertain factors on the wide area control system. F The larger the value, the greater the impact of fluctuations caused by uncertain factors.

[0147] By combining the reachable set results obtained from the reachability analysis with the two calculated performance indicators, the inter-regional frequency oscillations of the wide-area control system after the injection of different network attacks are compared, and the type of network attack that has the greatest impact on the damping performance of the wide-area control system can be obtained; in addition, by adding a network attack amplitude range model to the network attack model and a wind farm output variation range model to the wind farm model, the oscillation of the regional frequency of the wide-area control system under the dual influence of network attacks and wind power disturbances is calculated, reflecting the additional impact of the uncertainty of large-scale renewable energy on the effectiveness of network attacks.

[0148] Next, we take a 3-machine 9-node wide-area control system as an example to introduce the application method of the present invention in an actual power system.

[0149] First, the synchronous generators in area 3 of a traditional 3-machine, 9-bus power system are replaced with a wind farm consisting of multiple doubly-fed wind turbines. A wide-area damping controller is added to synchronous generator 1 to control the frequency oscillation between areas 1 and 2, forming a 3-machine, 9-bus wide-area control system. The 3-machine, 9-bus wide-area control system is modeled based on the models of its components, including the synchronous generators, wind farm, and wide-area damping controller. Considering the error accumulation caused by the interaction between differential and algebraic equations in the calculation of reachable sets, the node admittance matrix of the 3-machine, 9-bus wide-area control system is reduced to a third-order matrix containing only the generator nodes based on a simplified admittance matrix. Due to the inclusion of the wind farm, the coupling nodes of the wind farm must be retained, resulting in a fourth-order node admittance matrix, thus transforming the differential-algebraic equation model into a nonlinear differential equation model.

[0150] Considering the variability of cyber attack amplitudes and the impact of uncertain wind power disturbances, Zeno polyhedra are used to model the cyber attack amplitude range and wind farm output variation range, respectively. These are then added as uncertain input sets to the wide-area damping controller model and wind farm model, respectively. Based on the established nonlinear differential equation model of the wide-area control system with cyber attacks and wind power disturbances, the system operating state is specifically divided into four stages: normal operation, cyber attack injection, short-circuit fault occurrence, and short-circuit fault removal. A block calculation method is used to sequentially calculate the reachable set of the system state for each operating stage. The fluctuation domain of the wide-area control system regional frequency obtained from the reachable set results is used to quantitatively analyze the impact of the dual disturbance on the wide-area control system's damping performance. Furthermore, to further quantify the impact of uncertain factors, two performance indicators are calculated, reflecting the maximum deviation degree and maximum fluctuation range of the regional frequency.

[0151] The above description is merely a specific embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of this application. Therefore, the scope of protection of this application should be based on the scope of protection of the claims.

Claims

1. A quantitative assessment method for wide area control system network attacks based on reachability analysis, characterized by: The following steps are involved: Step 1: Based on the simplified admittance matrix, a nonlinear differential equation model of the wide-area control system with wind farm integration is established; Step 2: Considering the variability of the cyber attack amplitude and the uncertainty of wind power disturbances, the cyber attack amplitude range and wind farm output variation range are modeled using Zeno polyhedra. These models are then incorporated into the nonlinear differential equation model of the wide-area control system. Step 3: To effectively analyze the nonlinear differential equation model of the wide-area control system, an improved reachable set calculation algorithm is designed based on reachability analysis theory. The nonlinear differential equation model is converted into a linear differential equation model. The reachable set of the nonlinear differential equation model is obtained by combining the solution results of the linear differential equation model and the linearization error estimation results. The specific improvements to the reachable set calculation algorithm include: using the Krylov subspace iteration method based on the stable biconjugate gradient method to calculate the homogeneous solution of the linear differential equation model, and using the zonotope bundle technique to combine the Zeno polyhedra in the reachable set calculation process to improve the calculation accuracy of the reachable set. Step 4: Apply the improved reachable set calculation algorithm to calculate the reachable sets of the wide-area control system under the influence of different types of network attacks and wind power disturbances. Based on the reachable set results, quantitatively analyze the impact of network attacks and wind power disturbances on the damping performance of the wide-area control system, thereby completing the quantitative evaluation of network attacks on the wide-area control system.

2. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 1 is characterized in that: In step 1, the wide-area control system includes synchronous generators, wind farms, wide-area damping controllers, transformers, and power loads. Assume that the wide-area control system has N nodes. g Synchronous generators, N l load nodes and N w There are N wind farms. p A wide-area damping controller is deployed in the synchronous generator of the wide-area control system. The wide-area damping controller uses the collectible wide-area signal to perceive the regional frequency oscillation in the wide-area control system in real time, and provides additional damping for the wide-area control system by adjusting the active power output of the synchronous generator, thereby suppressing the regional frequency oscillation.

3. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 2 is characterized in that: In step 1, first, a nonlinear differential-algebraic equation model of the wide-area control system with wind farm access is constructed, considering the impact of cyber attacks. Secondly, based on the simplified admittance matrix, the node admittance matrix of the wide-area control system is reduced to an admittance matrix containing only the nodes of each generator. Finally, the nonlinear differential-algebraic equation model of the wide-area control system is converted into a nonlinear differential equation model.

4. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 3 is characterized by: The nonlinear differential-algebraic equation model includes: a wind farm model, a synchronous generator model, a wide-area damping controller model, a network attack model, and a power system flow model; The wind farm model is represented by a synchronous generator simulation model of a doubly-fed wind turbine generator. The synchronous generator simulation model of the i-th doubly-fed wind turbine generator is specifically represented as follows: Where i = 1, 2, ..., N w ; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the virtual rotor angle of the doubly fed wind turbine; is the derivative value of the internal electromotive force of the doubly fed wind turbine generator; θ ei 、ω ri 、E si are the virtual rotor angle, rotor speed and internal electromotive force of the doubly fed wind turbine generator respectively; is the reference value of the rotor speed of the doubly fed wind turbine generator; I ri0 、i qri They represent the initial value of the rotor current of the doubly fed wind turbine generator and the q-axis component of the rotor current respectively; K I1i With K p1i are the integral gain and proportional gain of the PI controller in the rotor-side converter of the doubly fed wind turbine generator; H i With L mi are the inertia constant and mutual inductance parameters of the doubly fed wind turbine generator respectively; the reactance of the doubly fed wind turbine generator is x si =ω s L si , L si is the stator inductance of the doubly fed wind turbine generator, ω s is the rated value of the rotor angle of the doubly fed wind turbine generator; V pi represents the voltage of the common coupling node connected to the doubly fed wind turbine generator, θ pi is the common coupling node voltage V pi Phase angle; P wmi is the mechanical power of the doubly fed wind turbine generator; The synchronous generator model adopts the classic second-order model, which is specifically expressed as follows: Where j = 1, 2, ..., N g , is the derivative value of the power angle of the jth synchronous generator, is the derivative value of the angular velocity of the jth synchronous generator, ω j represents the angular velocity of the jth synchronous generator; D j 、M j are the damping coefficient and moment of inertia of the j-th synchronous generator respectively; P mj is the mechanical power of the jth synchronous generator, P ej is the electromagnetic power output by the jth synchronous generator, specifically expressed as Among them E j is the internal potential of the synchronous generator, V j is the node voltage connected to the jth synchronous generator, δ j represents the power angle of the jth synchronous generator, θ j is the node voltage V j The phase angle of G gj With B gj are the conductance and susceptance values ​​of the line between the synchronous generator and the connected node, respectively; The wide-area damping controller model uses the difference in angular velocity of synchronous generators located in two different areas as a feedback signal and controls the output of the wide-area damping controller model through filtering, phase compensation, gain, and saturation adjustment. The specific expression of the wide-area damping controller model is: In the formula, the coefficient matrix Coefficient matrix B wm =[0;1], coefficient matrix is the derivative value of the state variable vector of the wide-area damping controller model, P WADCm represents the output value of the wide-area damping controller model, x wm is the state variable vector of the wide-area damping controller model, K m 、T m and T wm are the gain parameter of the wide-area damping controller model, the time constant of the phase compensation component, and the time constant of the filter; Δω ab is the angular velocity difference between the a-th and b-th synchronous generators, which is set as the control signal of the wide-area damping controller model; The output of the wide-area damping controller model can be used as part of the active power output command of the synchronous generator, enabling the synchronous generator to respond to regional frequency oscillations and enhance the damping effect of the wide-area control system. The network attack model considers the network attacks that the wide-area control system may suffer during the data acquisition and transmission process, namely, false data injection attacks and delay attacks, wherein the false data injection attacks include disturbance-type false data injection attacks and proportional-type false data injection attacks; when the control signal Δω of the wide-area damping controller model is ab After being affected by the injected network attack, Δω ab The control signal of the wide-area damping controller model after being affected by the attack is converted into The specific representation of the perturbation-type false data injection attack is as follows: Where t is the running time of the wide area control system, d is the injection component of the disturbance type false data injection attack, and t c Indicates the injection time of the perturbation-type false data injection attack; The specific representation of the proportional false data injection attack is as follows: Where, α sc is the proportional gain of proportional false data injection attack, t sc Indicates the injection time of proportional false data injection attack; For delayed attacks, the second-order modified Padé approximation is used to model the time delay, and the time delay model is obtained. The specific state space equation is expressed as: in, is the derivative value of the state variable vector of the time delay model, the coefficient matrix Coefficient matrix Coefficient matrix C d =[0,1],x d Represents the state variable vector of the time delay model, T d Indicates the delay time. represents the angular velocity difference of the synchronous generator after being affected by the time delay; the control signal of the wide-area damping controller model after being affected by the time delay attack is expressed as: Where, t d Indicates the injection time of the delayed attack; When the network attack component directly contaminates the control signal of the wide-area damping controller model, the wide-area damping controller model will become: Where, is the state variable vector of the wide-area damping controller model after being affected by the network attack, for The derivative value of is the abnormal output value of the wide-area damping controller model after being affected by the cyber attack; The power system flow model is specifically expressed as: Where V k is the voltage of the kth node of the wide area control system, V h is the voltage of the hth node of the wide-area control system, P k is the active power flowing into the kth node, P ek is the active power output by the synchronous generator at the kth node, P lk is the load power at the kth node, P wk is the output power of the wind farm at the kth node, θ k is the phase angle of the voltage at the kth node, θ h is the phase angle of the voltage at the hth node, B kh With G kh are the conductance and susceptance of the line between the kth node and the hth node respectively; The above-mentioned wind farm model, synchronous generator model, wide-area damping controller model, network attack model, and power system flow model are transformed into a nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system. At the same time, in order to record the time changes during the calculation process, an extended variable t representing the operating time of the wide-area damping control system is introduced. Its state equation is expressed as follows: When no network attack is added, define the state variable vector matrix Algebraic variables vector matrix Among them, δ p 、ω p ,θ ep 、ω rp 、E sp 、x wp denote the power angle of the synchronous generator at the p-th node, the angular velocity of the synchronous generator, the virtual rotor angle of the doubly fed wind turbine, the rotor speed of the doubly fed wind turbine, the internal electromotive force of the doubly fed wind turbine, and the state variable vector of the wide-area damping controller model, respectively. is the field of real numbers, E p 、V p ,θ p are the internal potential of the synchronous generator at the pth node, the node voltage amplitude and the node voltage phase angle respectively; the wide area control system can be modeled as: Where x(t), y(t) and u(t) are the state variable vector, algebraic variable vector and input variable set of the wide area control system under normal operating conditions, respectively. is the derivative value of the state variable vector of the wide-area control system under normal operating conditions, f(·) and g(·) represent the set of nonlinear differential equations and algebraic equations describing the system dynamic characteristics of the wide-area control system, respectively; After adding the network attack, the nonlinear differential-algebraic equation model of the wide-area control system is transformed into: Where x c (t), y c (t) and u c (t) are the state variable vector, algebraic variable vector and input variable set of the wide area control system after being affected by the network attack, is the derivative value of the state variable vector of the wide-area control system after being affected by the network attack; In order to solve the error accumulation problem caused by the interaction between differential and algebraic equations in the calculation of reachable sets, a simplified admittance matrix construction method is adopted to achieve the reduction of the node admittance matrix of the wide-area control system through matrix operation. Based on the obtained reduced admittance matrix, the electromagnetic power output by the j-th synchronous generator is expressed as: Where G jj is the self-conductance of the j-th synchronous generator, E j With E q are the internal potentials of the jth and qth synchronous generators, δ j and δ q are the power angles of the jth and qth synchronous generators respectively, B jq With G jq are the conductance and susceptance between the jth and qth synchronous generators, respectively. According to the above formula (11), it can be seen that the power distribution in the wide-area control system can be obtained without calculating the algebraic equations; In addition, considering the access of wind farms, it is necessary to retain the internal nodes of the wind farm and the common coupling nodes connected to the wind farm. Assuming that the common coupling node of the wind farm is the rth node in the wide-area control system, according to Kirchhoff's current law, we can obtain: Where G wp With B wp are the conductance and susceptance between the internal nodes of the wind farm and the common coupling node, G p0 With B p0 are the ground conductance and susceptance values ​​of the common coupling node, G gr With B gr are the conductance and susceptance between the common coupling node and the external g-th node, V g and δ g are the voltage and phase angle of the g-th node respectively; According to the above formula (12), the internal potential E of the wind farm can be si The voltage of the common coupling node is obtained, which is specifically expressed as: Where B qq is the self-susceptance of the common coupling node of the wind farm; the electromagnetic power output by the wind farm obtained by the above formula (13) can be used to establish the coupling relationship between the wind farm and other synchronous generators; Therefore, the nonlinear differential-algebraic equation model of the wide-area control system in the form of a differential-algebraic equation system can be transformed into a nonlinear differential equation model in the form of a differential equation system, which can be specifically expressed as: According to the above formula (14), without explicitly processing the algebraic equations, the dynamics of the wide-area control system can be obtained directly by solving the constructed nonlinear differential model.

5. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 4 is characterized in that: In step 2, considering the variability of the network attack amplitude and the impact of uncertain wind power disturbances caused by random changes in wind speed, Zeno polyhedrons are used to model the network attack amplitude range and wind farm output variation range respectively. The constructed models are added as uncertain input sets to the nonlinear differential equation model of the wide-area control system, thereby reflecting the impact of network attack amplitude changes and wind power disturbances on the damping performance of the wide-area control system.

6. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 5 is characterized by: Based on the construction method of Zeno polyhedron, the model of network attack amplitude range is expressed as: Where, Indicates the range of changes in the magnitude of network attacks; represents the amplitude center of the network attack amplitude, l is a variable used to record the change in the number of generators of the Zeno polyhedron during the modeling process, is the range of network attack amplitude variation represented by the generator of Zeno polyhedron, where α l is the proportionality coefficient, g l is the generator vector; w1 is the number of generators of the Zeno polyhedron used to model the range of network attack amplitude variations; by setting the Zeno polyhedron, a model of network attack amplitude range of any size can be generated, and this model is added to the network attack model to quantitatively analyze the impact of variable attacks; Considering the impact of wind power fluctuations on wide-area control systems, the Zeno polyhedron is used to model the wind farm output variation range and integrated into the wind farm model. The model of the wind farm output variation range is specifically expressed as: Where, The output variation range of the wind farm; is the center of change of wind farm output, set to 0; is the fluctuation range of wind farm output represented by the generator of Zeno polyhedron, where β l is the proportional coefficient, h l is the generator vector; w2 is the number of generators of the Zeno polyhedron used to model the output variation range of the wind farm.

7. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 6 is characterized in that: In step 3, the reachability analysis theory conservatively predicts the dynamic operating trajectory of the wide-area control system under different operating conditions to obtain the reachable set of the wide-area control system. The reachable set can effectively analyze the operating characteristics of the wide-area control system under random and uncertain working conditions. When performing reachability analysis on a wide-area control system with nonlinear characteristics, it is necessary to first approximate the nonlinear differential equation model into a linear differential equation model through the method of linear abstraction of the state space, and then use the reachable set calculation algorithm applicable to the linear differential equation model to perform calculations; for the reachable set calculation algorithm of the linear differential equation model, the Krylov subspace iteration method based on the stable biconjugate gradient method is used to improve the reachable set calculation algorithm of the linear differential equation model. Specifically, the Krylov subspace is defined as in, For vector The vector space generated by the linear combination of , span(·) represents the linear span algorithm that returns a set of vectors, n represents the dimension of the Krylov subspace, The variable vector represents the linear equation system in the linear differential equation model, and A is the coefficient matrix of the linear equation system. In order to obtain the homogeneous solution of the linear differential equation system, the problem of solving the linear differential equation system is transformed into the problem of solving the best approximation in the Krylov subspace. In order to improve the stability of the calculation, the stable biconjugate gradient method is used instead of the Arnoldi method to solve the problem, thereby obtaining the solution of the linear differential equation model. In addition, in order to compensate for the linearization error caused by the linear abstraction method of the state space, the over-approximation method is used to determine the linearization calculation error, and the linearization calculation error is added to the reachable set calculation result of the linear differential equation model in the form of an uncertain input set, and the reachable set calculation result of the nonlinear system is obtained within the error tolerance range. In order to improve the processing efficiency and calculation accuracy of Zeno polyhedra during the reachable set calculation process, the zonotopebundle technology is introduced to combine Zeno polyhedra during the reachable set calculation process. The zonotope bundle technology is specifically expressed as follows: Where o represents the number of Zeno polyhedra, is a reachable set described in terms of a Zeno polyhedron, is the reachable set after combining Zeno polyhedra using zonotope bundle technology, M is the total number of Zeno polyhedra; after calculating the reachable set results at each moment, the sth calculation step Δt s As a unit, apply the convex hull calculation CH(·) to obtain the reachable set within the specified calculation step size The calculation process is expressed as: in t s time and t s+1 The reachable set of time; Finally, in order to obtain the specified time range [t0,t f ] reachable set within The reachable sets computed at each individual time step need to be combined: where t0, t f The start and end times are set respectively.

8. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 7 is characterized in that: In order to realize the reachable set under various operating conditions of the wide-area control system, a block calculation method is adopted to calculate the reachable sets of different operating states in sequence, and the reachable set result of the previous operating state is used as the initial set of the next operating state. Finally, the union of the reachable sets of the four operating states consisting of normal operating state, attack state, fault state, and fault removal state is taken as the final calculation result, which is specifically expressed as: in, and Represent the reachable sets under four operating states respectively, is the final reachable set result, t0, t1, t2, t3, and t4 are the start and end times of the four operating states, respectively. The obtained reachable set result can reflect the frequency oscillation between regions of the wide-area control system under the influence of the set network attack and wind power disturbance.

9. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 8 is characterized in that: In step 4, based on the established nonlinear differential equation model of the wide-area control system, the improved reachable set calculation algorithm in step 3 is used to calculate the reachable sets of the wide-area control system affected by different types of network attacks and wind power disturbances. The obtained reachable set results are used to quantitatively analyze the impact of the set network attacks and wind power disturbances on the damping performance of the wide-area control system. At the same time, in order to further quantify the impact of uncertainties, two performance indicators are calculated to reflect the maximum deviation degree and maximum fluctuation range of the inter-regional frequency difference of the wide-area control system. The specific conditions of these two performance indicators are as follows: Indicator-J D : Where: e(t) is Δω ab The difference between the instantaneous value and the steady-state value, t1 and t2 are the starting time and the ending time of the indicator calculation respectively; this indicator is used to capture Δω ab The maximum deviation, J D The higher the value of , the worse the damping performance of the wide-area control system; Indicator 2J F : Where: Δω abmax and Δω abmin Represent Δω respectively ab By calculating the maximum and minimum values ​​of this index, we can effectively reflect the degree of influence of various uncertain factors on the wide area control system. F The larger the value, the greater the impact of fluctuations caused by uncertain factors.

10. The method for quantitatively assessing network attacks on a wide area control system based on reachability analysis according to claim 9 is characterized in that: By combining the reachable set results with the two performance index results, the inter-regional frequency oscillations of the wide-area control system after the injection of different types of network attacks are compared, so that the type of network attack that has the greatest impact on the damping performance of the wide-area control system can be obtained; in addition, by adding a network attack amplitude range model to the network attack model and a wind farm output variation range model to the wind farm model, the inter-regional frequency oscillations of the wide-area control system under the dual influence of network attacks and wind power disturbances are calculated, reflecting the additional impact of the uncertainty of large-scale renewable energy on the effect of network attacks; based on the reachable set results and the two performance index results, the vulnerability of the wide-area control system under different attack scenarios and operating conditions can be evaluated, providing important guidance for improving the damping performance of the wide-area control system and ensuring the safe and stable operation of the power system.

Images