Urban power grid chain risk quantification method, system, equipment and medium

Abstract

The invention discloses an urban power grid chain risk quantification method, system and device and a medium, and the method comprises the steps: obtaining the fault probability of each device in each scene based on the meteorological data and network security data of each link at the current moment in combination with a pre-constructed multi-drive-fault probability mapping relation, and determining a tripping device; when the equipment trips, acquiring meteorological data and network security data of each link at the current moment and a power grid state vector at the previous moment; based on the state vector of the power grid at the previous moment, combining with a pre-constructed full-link chain transmission model to obtain a state vector at the current moment; calculating the total risk of the system by adopting an analytic hierarchy process based on the state vector at the current moment; wherein the full-link chain transmission model is constructed based on a state association relationship in historical data of a power grid state vector, a coupling relationship between links and influences of meteorological and network security factors on states of the links. According to the method, the cross-link and cross-time scale dynamic coupling effect is quantified, and the risk assessment precision and timeliness are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of power system safety evaluation, in particular to a city power grid chain risk quantification method, system, device and medium. BACKGROUND

[0002] With the evolution of urban new power systems to high proportion of renewable energy and high proportion of power electronic devices ("double high"), such systems have chain-type unsafe state risks. The volatility of high proportion of wind and light power sources and extreme weather are coupled, which easily leads to multi-link chain failure of "power grid-load-energy storage" (such as new energy off-grid→voltage collapse→load interruption). At the same time, the "source network load storage" collaborative network relied on digital dispatch may become a risk amplifier under network security attacks, leading to cascading effects of supply chain nodes (such as energy storage, carbon trading center) failure. The existing technology lacks quantitative means for such cross-link and cross-time scale chain risks, and only evaluates single links (such as power grid transient stability) or static scenarios (such as device redundancy), which cannot reflect the dynamic coupling characteristics of multiple types of power source complementary failure, multi-agent interaction vulnerability, and low-carbon technology maturity difference in new power systems, resulting in delayed city power safety risk warning and inefficient emergency response.

[0003] In recent years, major power accidents caused by chain-type unsafe states have occurred in many places. For example, in 2021, during a typhoon in a coastal city, photovoltaic power stations were off-grid due to sudden wind speed increase, causing voltage fluctuations in the distribution network, and superimposing energy storage system communication interruption (caused by network attacks), ultimately causing commercial area load interruption for 4 hours. In 2022, during the cold wave in the north, wind power output dropped, triggering the overload of thermal power units, while the carbon trading center could not adjust the carbon emission quota in real time due to data transmission failure, causing the delay of emergency power supply start, which amplified the load loss. Such accidents highlight the limitations of the existing evaluation system - it cannot capture the chain logic driven by "natural disturbance-network attack", and it is difficult to quantify the dynamic conflict between "low-carbon target-safety redundancy".

[0004] The existing power system safety evaluation technology has many defects in dealing with the risk assessment of urban new power systems with high proportion of renewable energy access under extreme weather conditions, mainly in the following three aspects:

[0005] Single-link assessment limitation: The traditional N-1 criterion is a commonly used method for power system security assessment, which checks the security and stability of the system after the failure of a single component (such as a generator, transmission line, transformer, etc.). However, in urban new-type power systems with a high proportion of renewable energy access, there is a close coupling relationship between the generation and consumption links. The failure of a certain link may spread to other links through a complex network, triggering a chain reaction. For example, when the output of distributed power generation suddenly decreases due to extreme weather, it will not only affect local power supply, but also cause changes in power distribution of transmission lines, thereby affecting power consumption in other areas. The traditional N-1 criterion only focuses on the failure of a single component and ignores the coupling effect of the generation-consumption chain, making it impossible to comprehensively assess the overall risk of the system. In some extreme cases, a single component failure assessment may show that the system is safe, but a large-scale power outage may occur due to a chain reaction.

[0006] Coarse weather modeling: Weather conditions are an important factor affecting the safe operation of power systems, especially extreme weather, which has a significant impact on transmission lines, towers, and other equipment. However, the spatial resolution of meteorological data used in existing technologies is usually greater than 10 km, which cannot accurately reflect the local wind load impact of transmission lines caused by micro-topography. Micro-topography (such as valleys, high-rise buildings, etc.) can cause sharp changes in local wind speed and direction, forming strong turbulence and causing large local wind loads on transmission lines. For example, in valley areas, the wind speed may be much higher than in surrounding areas due to the channeling effect of the terrain, causing transmission lines in this area to bear more tension and increasing the risk of line breakage and tower collapse. Coarse weather modeling cannot capture these local weather characteristics, resulting in large errors in the assessment of equipment failure probability and affecting the accuracy of risk assessment.

[0007] Single risk indicator: Existing risk assessment methods mostly use static indicators such as expected energy not supplied (EENS) to measure the risk level of power systems. EENS mainly reflects the total amount of expected power shortage in a certain period of time, which can reflect the power supply reliability of the system to some extent. However, it is a static indicator and cannot reflect the dynamic changes of the system during the failure process. The failure evolution of a power system is a dynamic process that may result in voltage collapse, frequency instability, and other dynamic phenomena. For example, when a major failure occurs in the system, voltage and frequency may fluctuate dramatically, and if not controlled in time, it may lead to the overall collapse of the system. Static indicators such as EENS cannot quantify these dynamic processes, making it difficult to comprehensively and accurately assess the risk of the system and provide timely and effective information for emergency decision-making.

[0008] Therefore, there is an urgent need for a method to quantify the probability and impact of the unsafe state chain of the "source-grid-load-storage-carbon" whole chain, which can support the resilience optimization and low-carbon transformation safety of the new urban power system. SUMMARY

[0009] In order to solve the problem that the prior art lacks quantitative means for chain risk across links and time scales, only evaluates a single link (such as power grid transient stability) or a static scenario (such as device redundancy), and cannot reflect the dynamic coupling characteristics in the new power system, resulting in delayed urban power safety risk warning and inefficient emergency response, the present application provides a city grid chain risk quantification method, comprising:

[0010] Based on the meteorological data and network security data of each link at the current time and the pre-constructed multi-driving-fault probability mapping relationship, the fault probability of each device under each scenario is obtained, and the tripping device is determined based on the fault probability of each device under each scenario.

[0011] When the device is tripped, the meteorological data, network security data and the power grid state vector of the previous time of each link at the current time are obtained.

[0012] Based on the power grid state vector of the previous time and the pre-constructed whole-link chain model, the state vector of each link at the current time is obtained.

[0013] Based on the state vector of each link at the current time, the total risk of the system is calculated by using the analytic hierarchy process.

[0014] The whole-link chain model is constructed based on the state correlation relationship in the historical data of the power grid state vector, the coupling relationship between each link, and the influence of meteorological factors and network security factors on the state of each link.

[0015] Optionally, the construction process of the whole-link chain model comprises:

[0016] Based on the historical data of the state vector of each link, the weight matrix of each link is obtained by training with a deep neural network.

[0017] Based on the correlation relationship between the historical data of the state vector of each link, the initial matrix composed of the historical data of the state vector of each link is modified to obtain the coupling matrix of each link to the next link.

[0018] Based on the historical data of the state vector, a joint disturbance function is constructed based on the functional relationship between the meteorological factors and the state vector in each link and the functional relationship between the network security attack factors and the state vector in each link.

[0019] The full-link chain model is constructed based on the weight matrix of each link, the coupling matrix of each link, the joint disturbance function and the state vector historical data of each link.

[0020] Optionally, the state vector of each link is calculated based on the grid operation parameters.

[0021] Optionally, the joint disturbance function is constructed based on the state vector historical data, the functional relationship between the meteorological factor and the state vector in each link and the functional relationship between the network security attack factor and the state vector in each link, comprising:

[0022] The functional relationship between the meteorological factor and the state vector in each link and the functional relationship between the network security attack factor and the state vector in each link are combined based on the state vector historical data to calculate the meteorological factor influence sub-function and the network security attack factor influence sub-function.

[0023] The weights of the meteorological factor influence sub-function and the network security attack factor influence sub-function are determined according to the contribution degrees of the meteorological factor and the network security attack in the historical faults.

[0024] The meteorological factor influence sub-function and the network security attack factor influence sub-function and the corresponding weights are weighted and fused to obtain the joint disturbance function.

[0025] Optionally, the calculation formula of the joint disturbance function is as follows:

[0026] F(M(t),A(t))=α×F m (M(t))+(1-α)×F a (A(t))

[0027] In the formula, F(M(t),A(t)) is the joint disturbance function, α is the weight coefficient, F m (M(t)) is the functional relationship between the meteorological factor and the state variable, F a (A(t)) is the functional relationship between the network security attack factor and the state variable.

[0028] Optionally, the full-link chain model is as follows:

[0029] X i (t+Δt)=ReLU(W i ×X i (t)+Σ(J ij ×X j (t))+F(M(t),A(t))×Δt)

[0030] In the formula, X i (t+Δt) is the state vector of link i at time t+Δt, X i(t) is the state vector of the link j at time t, X j (t) is the state vector of the link j at time t, W i is the weight matrix of the link i, J ij is the coupling matrix of the link j to the link i, F(M(t), A(t)) is the joint disturbance function of weather and network security attacks, and Δt is the time step.

[0031] Optionally, the historical data of the state vectors of each link are combined with a deep neural network to train and obtain the weight matrix of each link, including:

[0032] Obtaining historical data of state vectors of each link, the historical data of state vectors of each link being sampled at intervals of Δt, including operation data under different seasons, different weather conditions, and different network security states;

[0033] Performing outlier rejection or correction and missing value supplement processing on the historical data of state vectors of each link to obtain processed data;

[0034] Constructing a sample set from the processed data, each sample including state vectors of the same link at times t and t+Δt, t being a time;

[0035] Dividing the sample set into a training set, a validation set, and a test set according to a set proportion;

[0036] Defining the number of input layer nodes and the number of output layer nodes of the deep neural network according to the dimension of the state vectors of each link; taking the state vector at time t as the input of the deep neural network and taking the state vector at time t+Δt as the output;

[0037] Training, validating, and testing the deep neural network based on the training set, the validation set, and the test set to obtain the weight matrix of each link.

[0038] Optionally, the training, validating, and testing of the deep neural network based on the training set, the validation set, and the test set to obtain the weight matrix of each link include:

[0039] Using the training set to perform forward propagation training on the deep neural network, taking the mean square error as a loss function to calculate the training loss;

[0040] Calculating the gradient by back propagation and updating the weight matrix W i of the deep neural network using an Adam optimizer.

[0041] After completing K training rounds, evaluating the validation loss value of the current deep neural network model using the validation set, K being a preset period.

[0042] If the variation of the validation loss value of the N continuous periods is less than a set threshold, the training is stopped, a deep neural network model passing the validation is obtained, N is a preset positive integer;

[0043] The generalization performance of the deep neural network model passing the validation is evaluated by using the test set, and the weight matrix W output by the deep neural network model meeting the set threshold requirement of generalization performance is obtained i As the weight matrix of each link.

[0044] Optionally, the total risk of the system is calculated by using the analytic hierarchy process based on the state vector of the current moment of each link, comprising:

[0045] Based on the power loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission exceeding loss, network security risk degree, and dynamic risk calculation formula, the dynamic risk value of each node at each time in each scene is calculated;

[0046] Based on the dynamic risk value of each node at each time in each scene, the accumulated risk contribution degree of each node in the chain process of each scene is calculated by using the accumulated risk contribution degree calculation formula;

[0047] Based on the accumulated risk contribution degree, the total risk of the system is calculated by using the total risk calculation formula;

[0048] The state vector includes: power loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission exceeding loss, and network security risk degree.

[0049] Optionally, the dynamic risk calculation formula is as follows:

[0050] DRI i,j (t)=ω1×L i,j (t)+ω2×U i,j (t)+ω3×F i,j (t)+ω4×E i,j (t)+ω5×C i,j (t)+ω6×S i,j (t)

[0051] In the formula, DRI i,j (t) is the dynamic risk value of node i in scene j and time t, L i,j (t) is the relative power supply of node i in scene j and time t, U i,j (t) is the voltage out-of-limit degree, F i,j (t) is the frequency out-of-limit degree, E i,j (t) is the energy storage reserve capacity shortage degree, C i,j (t) is the carbon emission exceeding loss degree, S i,j(t) is a network security risk degree, ω1 is a weight coefficient of the relative power supply shortage of node i at scene j and time t, ω2 is a weight coefficient of the voltage overrun degree, ω3 is a weight coefficient of the frequency overrun degree, ω4 is a weight coefficient of the energy storage reserve capacity shortage degree, ω5 is a weight coefficient of the carbon emission overrun loss degree, ω6 is a weight coefficient of the network security risk degree, i is a node serial number, j is a scene serial number, and t is a time value.

[0052] Optionally, the accumulated risk contribution degree calculation formula is as follows:

[0053] CVI i,j =∫(t0 to t1)DRI i,j (t)×e^(-λ×(t-t0))dt

[0054] In the formula, CVI i,j is the accumulated risk contribution degree of node i in the chain process of scene j, t0 is the initial time of failure, t1 is the first time of system overrun, λ is a decay coefficient dynamically adjusted with the time scale, to represents to or to, and e is a natural constant.

[0055] Optionally, the total risk calculation formula is as follows:

[0056] SR j =Σ(CVI i,j ×α i ×β i )

[0057] In the formula, SR j is the total risk of the system, α i is the sensitive load weight supplied by node i, and β i is the load level coefficient.

[0058] Optionally, the construction of the multi-drive-fault probability mapping relationship comprises:

[0059] Based on the current time meteorological data, each extreme weather scenario is generated by using a WRF-CFD coupling model;

[0060] Based on the current time network security data, each network security scenario is generated by using an MITREATT&CK framework;

[0061] Each extreme weather scenario and each network security scenario are aligned according to a time axis to generate an extreme weather-network security joint scenario set;

[0062] An improved K-medoids algorithm is used to reduce the extreme weather-network security joint scenario set to obtain key scenarios;

[0063] Based on the device physical characteristics and data driving in the key scene, a function relationship of extreme weather, network security attack and power equipment failure probability is established;

[0064] The function relationship of extreme weather, network security attack and power equipment failure probability includes a relationship formula of transmission line failure probability, a relationship formula of photovoltaic inverter failure probability and a relationship formula of energy storage battery failure probability.

[0065] Optionally, the relationship formula of transmission line failure probability is as follows:

[0066] P t Failure = P ase _ base × [1 + β1 × (v(t) / v nitical -1) + β2 × (i(t)i nitical + β3 × p(t)]

[0067] In the formula, P t Failure is the transmission line failure probability, P ase _base is the basic failure rate, β1 is the first coefficient, v(t) is the 1km resolution real-time wind speed output by CFD, v nitical is the line design critical wind speed, β2 is the second coefficient, i(t) is the ice thickness, i nitical is the line design critical ice thickness, β3 is the third coefficient, and p(t) is the data tampering probability.

[0068] The relationship formula of the photovoltaic inverter failure probability is as follows:

[0069] P +＝ PV Failure = P ase _PV × [1 + β4 × (t(t)t max -1) + β5 × s(t)]

[0070] In the formula, P +＝ PV Failure is the photovoltaic inverter failure probability, P ase _PV is the basic failure rate of the photovoltaic inverter, β4 is the fourth coefficient, t(t) is the environmental temperature, t max is the maximum allowable temperature of the inverter, β5 is the fifth coefficient, and s(t) is the infection rate of the inverter.

[0071] The relationship formula of the energy storage battery failure probability is as follows:

[0072] P e Failure = P 8ase _Storage × [1 + β6 × (1-SOC) + β7 × r(t)]

[0073] In the formula, P e Failure is the energy storage battery failure probability, P aseThe energy storage is a basic failure rate of the energy storage battery, β6 is a sixth coefficient, SOC is a state of charge, β7 is a seventh coefficient, and r(t) is a malicious code propagation rate.

[0074] Optionally, the method further comprises:

[0075] According to the total system risk from large to small, the scenarios corresponding to the first set percentage of the total system risk are high-risk scenarios.

[0076] According to the cumulative risk contribution degree from large to small, the nodes corresponding to the second set percentage of the cumulative risk contribution degree are high-risk nodes.

[0077] In still another aspect, the present application also provides a city power grid chain risk quantification system, comprising:

[0078] A fault probability evaluation module is configured to obtain device fault probabilities under each scenario based on meteorological data and network security data of each link at the current time and a multi-driving-fault probability mapping relationship, and determine tripped devices based on the device fault probabilities under each scenario.

[0079] A parameter acquisition module is configured to acquire meteorological data, network security data of each link at the current time, and a power grid state vector at the previous time when a device is tripped.

[0080] A chain module is configured to obtain a state vector of each link at the current time based on the power grid state vector at the previous time and a pre-constructed full-link chain model.

[0081] A risk quantification module is configured to calculate a total system risk based on the state vector of each link at the current time using an analytic hierarchy process.

[0082] The full-link chain model is constructed based on state correlation relationships in historical data of power grid operation parameters, coupling relationships between links, and influences of meteorological factors and network security factors on link states.

[0083] The power grid operation parameters are used to calculate link state vectors.

[0084] The historical data of link state vectors are used to train a deep neural network to obtain a weight matrix of each link.

[0085] The initial matrix composed of link state vector historical data is modified based on correlation relationships between the link state vector historical data to obtain a coupling matrix of each link to the next link.

[0086] The joint disturbance function is constructed based on historical data of the state vectors, function relationships between meteorological factors and the state vectors in each link, and function relationships between network security attack factors and the state vectors in each link.

[0087] The full-link chain model is constructed based on the weight matrix of each link, the coupling matrix of each link, the joint disturbance function, and historical data of the state vectors of each link.

[0088] Optionally, the full-link chain model is as follows:

[0089] X i (t+Δt)=ReLU(W i ×X i (t)+Σ(J ij ×X j (t))+F(M(t),A(t))×Δt)

[0090] In the formula, X i (t+Δt) is the state vector of link i at time t+Δt, X i (t) is the state vector of link i at time t, X j (t) is the state vector of link j at time t, W i is the weight matrix of link i, J ij is the coupling matrix from link j to link i, F(M(t), A(t)) is the joint disturbance function of meteorological factors and network security attack factors, and Δt is the time step.

[0091] Optionally, the specific implementation steps of constructing the joint disturbance function based on historical data of the state vectors, function relationships between meteorological factors and the state vectors in each link, and function relationships between network security attack factors and the state vectors in each link in the model construction module include:

[0092] Based on historical data of the state vectors, function relationships between meteorological factors and the state vectors in each link, and function relationships between network security attack factors and the state vectors in each link, meteorological factor influence sub-functions and network security attack factor influence sub-functions are calculated.

[0093] The weights of the meteorological factor influence sub-functions and the network security attack factor influence sub-functions are determined according to the contribution degrees of meteorological factors and network security attacks in historical faults.

[0094] The joint disturbance function is obtained by weighted fusion based on the meteorological factor influence sub-functions and the network security attack factor influence sub-functions and the corresponding weights.

[0095] Optionally, the calculation formula of the joint disturbance function is as follows:

[0096] F(M(t), A(t)) = a x F m (M(t)) + (1-a) x F a (A(t))

[0097] In the formula, F(M(t), A(t)) is a joint disturbance function, a is a weight coefficient, F m (M(t)) is a function relationship between weather factors and state variables, F a (A(t)) is a function relationship between network security attack factors and state variables.

[0098] Optionally, the model construction module trains based on historical data of state vectors of each link in combination with a deep neural network to obtain a weight matrix of each link, and the specific implementation steps include:

[0099] Obtain historical data of state vectors of each link, the historical data of state vectors of each link being sampled at an interval of Δt and including operation data under different seasons, different weather conditions and different network security states;

[0100] Perform outlier rejection or correction and missing value supplement processing on the historical data of state vectors of each link to obtain processed data;

[0101] Construct a sample set from the processed data, each sample including state vectors of the same link at times t and t+Δt, t being a time;

[0102] Divide the sample set into a training set, a validation set and a test set according to a set proportion;

[0103] Define the number of input layer nodes and the number of output layer nodes of the deep neural network according to the dimension of the state vectors of each link; take the state vector at time t as the input of the deep neural network and take the state vector at time t+Δt as the output;

[0104] Train, validate and test the deep neural network based on the training set, the validation set and the test set to obtain the weight matrix of each link.

[0105] Optionally, the risk quantification module is specifically configured to:

[0106] Calculate the dynamic risk value of each node at each time in each scenario based on power shortage loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission exceeding loss and network security risk degree in combination with a dynamic risk calculation formula;

[0107] Calculate the accumulated risk contribution degree of each node in the chain process of each scenario based on the dynamic risk value of each node at each time in each scenario in combination with an accumulated risk contribution calculation formula;

[0108] The accumulated risk contribution degree is combined with a total risk calculation formula to calculate the system total risk.

[0109] The state vectors of the links include power loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission exceeding loss, and network security risk degree.

[0110] Optionally, the dynamic risk calculation formula is as follows:

[0111] DRI i,j (t) = ω1*L i,j (t) + ω2*U i,j (t) + ω3*F i,j (t) + ω4*E i,j (t) + ω5*C i,j (t) + ω6*S i,j (t)

[0112] In the formula, DRI i,j (t) is a dynamic risk value of node i in scenario j and at time t, L i,j (t) is a relative power supply shortage of node i in scenario j and at time t, U i,j (t) is a voltage out-of-limit degree, F i,j (t) is a frequency out-of-limit degree, E i,j (t) is an energy storage reserve capacity deficiency degree, C i,j (t) is a carbon emission exceeding loss degree, S i,j (t) is a network security risk degree, ω1 is a weight coefficient of the relative power supply shortage of node i in scenario j and at time t, ω2 is a weight coefficient of the voltage out-of-limit degree, ω3 is a weight coefficient of the frequency out-of-limit degree, ω4 is a weight coefficient of the energy storage reserve capacity deficiency degree, ω5 is a weight coefficient of the carbon emission exceeding loss degree, ω6 is a weight coefficient of the network security risk degree, i is a node number, j is a scenario number, and t is a time value.

[0113] Optionally, the accumulated risk contribution degree calculation formula is as follows:

[0114] CVI i,j = ∫(t0 to t1) DRI i,j (t) * e^(-λ*(t-t0)) dt

[0115] In the formula, CVI i,j is an accumulated risk contribution degree of node i in a chain process of scenario j, t0 is a fault initial time, t1 is a system first out-of-limit time, λ is a decay coefficient dynamically adjusted with a time scale, to represents to or to, and e is a natural constant.

[0116] Optionally, the total risk calculation formula is as follows:

[0117] SR j =∑(CVI i,j ×α i ×β i )

[0118] wherein, SR j is a total risk of the system, α i is a weight of sensitive load supplied by the node i, and β i is a load level coefficient.

[0119] Optionally, the method further comprises a mapping relationship construction module, configured to:

[0120] generating each extreme weather scenario based on weather data at a current time using a WRF-CFD coupling model;

[0121] generating each network security scenario based on network security data at the current time using an MITRE ATT&CK framework;

[0122] aligning each extreme weather scenario and each network security scenario according to a time axis to generate a set of extreme weather-network security joint scenarios;

[0123] performing scenario reduction on the set of extreme weather-network security joint scenarios using an improved K-medoids algorithm to obtain key scenarios;

[0124] establishing a functional relationship between extreme weather, network security attacks and power equipment failure probability based on physical characteristics of devices and data driving in the key scenarios;

[0125] The functional relationship between extreme weather, network security attacks and power equipment failure probability includes a relationship formula of transmission line failure probability, a relationship formula of photovoltaic inverter failure probability and a relationship formula of energy storage battery failure probability.

[0126] Optionally, the calculation formula of the transmission line failure probability is as follows:

[0127] P t failure=P ase ×[1+β1×(v(t) / v nitical -1)+β2×(i(t)i nitical )+β3×p(t)]

[0128] wherein, P t failure is the transmission line failure probability, P ase is a basic failure rate, β1 is a first coefficient, v(t) is a 1km resolution real-time wind speed output by CFD, v nitical is a line design critical wind speed, β2 is a second coefficient, i(t) is an icing thickness, i nitical is a basic icing thickness, and β3 is a third coefficient.The critical icing thickness is designed for the line, beta3 is the third coefficient, and p(t) is the data tampering probability.

[0129] The calculation formula of the photovoltaic inverter failure probability is as follows:

[0130] P +＝ Photovoltaic failure = P ase _ photovoltaic x [1 + beta4(t(t) t max -1) + beta5 x s(t)]

[0131] In the formula, P +＝ failure is the photovoltaic inverter failure probability, P ase _ photovoltaic is the basic failure rate of the photovoltaic inverter, beta4 is the fourth coefficient, t(t) is the environmental temperature, t max is the maximum allowable temperature of the inverter, beta5 is the fifth coefficient, and s(t) is the inverter infection rate.

[0132] The calculation formula of the energy storage battery failure probability is as follows:

[0133] P e failure = P ase _ energy storage x [1 + beta6 x (1-SOC) + beta7 x r(t)]

[0134] In the formula, P e failure is the energy storage battery failure probability, P ase _ energy storage is the basic failure rate of the energy storage battery, beta6 is the sixth coefficient, SOC is the state of charge, beta7 is the seventh coefficient, and r(t) is the malicious code propagation rate.

[0135] In another aspect, the application further provides an electronic device, comprising: at least one processor and a memory; the memory and the processor are connected through a bus;

[0136] The memory is used for storing one or more programs.

[0137] When the one or more programs are executed by the at least one processor, the urban power grid chain risk quantification method is realized.

[0138] In another aspect, the application further provides a readable storage medium, which has an execution program stored thereon, and the execution program, when executed, realizes the urban power grid chain risk quantification method.

[0139] Compared with the prior art, the application has the following beneficial effects:

[0140] A city power grid chain failure risk quantification method, comprising: obtaining the failure probability of each device under each scenario based on the meteorological data and network security data of each link at the current time combined with the pre-constructed multi-driving-failure probability mapping relationship, and determining the tripping device based on the failure probability of each device under each scenario; obtaining the meteorological data and network security data of each link at the current time, and the power grid state vector at the previous time when the device trips; obtaining the state vector of each link at the current time based on the power grid state vector at the previous time combined with the pre-constructed full-link chain failure model; calculating the total risk of the system based on the state vector of each link at the current time using the analytic hierarchy process; wherein the full-link chain failure model is constructed based on the state correlation relationship in the historical data of power grid operation parameters, the coupling relationship between each link, and the influence of meteorological factors and network security factors on the state of each link. This method uses a full-link chain failure model, considers the complex energy flow and information flow exchange relationship between the power generation link (including centralized power stations and distributed power sources), the power transmission link, the power distribution link and the power consumption link, and clearly defines the relationship between the state variables of each link and the failure propagation law. In addition, the analytic hierarchy process is used to quantify the dynamic risk contribution of each node in the chain failure process. BRIEF DESCRIPTION OF DRAWINGS

[0141] Figure 1 A city power grid chain failure risk quantification method flow chart of the present application;

[0142] Figure 2 An electronic device structure schematic diagram of the present application. DETAILED DESCRIPTION

[0143] The existing power system safety evaluation technology has many defects in dealing with the risk assessment of the city new-type power system with high proportion of renewable energy access under extreme weather conditions, mainly in the following three aspects:

[0144] Single-link evaluation limitation: the traditional N-1 criterion is a commonly used method in power system safety evaluation, and its core idea is to check the safety and stability of the system after the failure of a single component (such as a generator, a transmission line, a transformer, etc.). However, in the city new-type power system with high proportion of renewable energy access, there is a close coupling relationship between the power generation link and the power consumption link, and the failure of a certain link may propagate to other links through complex network, causing a chain reaction. For example, when the output of distributed power sources suddenly decreases due to extreme weather, it will not only affect local power supply, but also may cause changes in power distribution of transmission lines, and further affect power consumption in other areas.

[0145] The traditional N-1 criterion is a commonly used method in power system safety evaluation, and its core idea is to check the safety and stability of the system after the failure of a single component (such as a generator, a transmission line, a transformer, etc.). However, in the city new-type power system with high proportion of renewable energy access, there is a close coupling relationship between the power generation link and the power consumption link, and the failure of a certain link may propagate to other links through complex network, causing a chain reaction. For example, when the output of distributed power sources suddenly decreases due to extreme weather, it will not only affect local power supply, but also may cause changes in power distribution of transmission lines, and further affect power consumption in other areas.

[0146] The N-1 criterion only focuses on the failure of a single element and ignores the coupling effect of the entire power generation and power consumption chain, and thus cannot comprehensively evaluate the overall risk of the system. In some extreme cases, a single element failure evaluation result may be safe, but a large-scale power outage of the system may occur due to cascading effects.

[0147] Coarse weather modeling: Weather conditions are an important factor affecting the safe operation of a power system, especially extreme weather, which has a significant impact on equipment such as transmission lines and towers. However, the spatial resolution of meteorological data used in existing technologies is usually greater than 10 km, which cannot accurately reflect the local wind load impact of micro-topography on transmission lines. Micro-topography (such as valleys, high-rise buildings, etc.) can cause sharp changes in local wind speed and direction, forming strong turbulence and causing large local wind loads on transmission lines. For example, in a valley area, due to the channeling effect of the terrain, the wind speed can be much higher than in the surrounding areas, causing the transmission lines in this area to bear more tension and increasing the risk of line breakage and tower collapse. Coarse weather modeling cannot capture these local weather characteristics, resulting in large errors in the evaluation of equipment failure probability and affecting the accuracy of risk assessment.

[0148] Risk indicator single:

[0149] One: Existing risk assessment methods mostly use static indicators such as expected energy not supplied (EENS) to measure the risk level of a power system. EENS mainly reflects the total amount of expected power shortage in the system within a certain period of time, and although it can reflect the power supply reliability of the system to some extent, it is a static indicator and cannot reflect the dynamic changes of the system during the failure process. The failure evolution of a power system is a dynamic process, which can result in a series of dynamic phenomena such as voltage collapse and frequency instability. For example, when a major failure occurs in the system, the voltage and frequency will fluctuate dramatically, and if not controlled in time, it can lead to the overall collapse of the system. Static indicators such as EENS cannot quantify these dynamic processes, making it difficult to comprehensively and accurately assess the risk of the system and provide timely and effective information for emergency decision-making.

[0150]

[0151] In order to better understand the present application, the contents of the present application will be further described below in conjunction with the drawings and examples of the specification. The present application provides a method for quantifying the cascading risk of a city power grid, as shown in FIG. 1, which includes the following steps:

[0152] In order to better understand the present application, the contents of the present application will be further described below in conjunction with the drawings and examples of the specification. The present application provides a method for quantifying the cascading risk of a city power grid, as shown in FIG. 1, which includes the following steps: Figure 1

[0153] ​​Step 1: Based on the current time of each link meteorological data and network security data combined with the pre-constructed multi-drive-fault probability mapping relationship, the fault probability of each device under each scene is obtained, and the tripping device is determined based on the fault probability of each device under each scene;

[0154] Step 2: When the device trips, the current time of each link meteorological data, network security data, and the last time of the power grid state vector are obtained;

[0155] Step 3: Based on the last time of the power grid state vector combined with the pre-constructed full-link chain model, the state vector of the current time of each link is obtained;

[0156] Step 4: Based on the state vector of the current time of each link, the total risk of the system is calculated by using the analytic hierarchy process;

[0157] Wherein, the full-link chain model is constructed based on the state correlation relationship in the historical data of power grid operation parameters, the coupling relationship between each link, and the influence of meteorological factors and network security factors on the state of each link.

[0158] The present application focuses on the risk quantification of chain failure type unsafe state under extreme weather for the new urban power system, and needs to solve three core problems: high-precision meteorological-fault mapping, full-link chain model construction and dynamic risk contribution quantification.

[0159] Before step 1, it also includes constructing a multi-drive-fault probability mapping relationship, which will be described in detail below.

[0160] The construction of multi-drive-fault probability mapping relationship includes:

[0161] Based on the current time of meteorological data, WRF-CFD coupling model is used to generate each extreme weather scenario;

[0162] Based on the current time of network security data, MITREATT&CK framework is used to generate each network security scenario;

[0163] Align each extreme weather scenario and each network security scenario according to the time axis to generate an extreme weather-network security joint scene set;

[0164] The improved K-medoids algorithm is used to reduce the scene of the extreme weather-network security joint scene set to obtain the key scene;

[0165] Based on the physical characteristics of the device in the key scene and data-driven, the functional relationship between extreme weather, network security attack and power device fault probability is established;

[0166] The function relationship of the extreme weather, the network security attack and the power equipment failure probability includes a relationship formula of a transmission line failure probability, a relationship formula of a photovoltaic inverter failure probability and a relationship formula of a storage battery failure probability.

[0167] Further, the calculation formula of the transmission line failure probability is as follows:

[0168] P t Failure=P ase _0×[1+β1×(v(t) / v nitical -1)+β2×(i(t)i nitical )+β3×p(t)]

[0169] In the formula, P t Failure is the transmission line failure probability, P ase _0 is a basic failure rate, β1 is a first coefficient, v(t) is a 1km resolution real-time wind speed output by a CFD, v nitical is a line design critical wind speed, β2 is a second coefficient, i(t) is an icing thickness, i nitical is a line design critical icing thickness, β3 is a third coefficient, and p(t) is a data tampering probability.

[0170] The calculation formula of the photovoltaic inverter failure probability is as follows:

[0171] P +＝ PV Failure=P ase _PV×[1+β4×(t(t)t max -1)+β5×s(t)]

[0172] In the formula, P +＝ PV Failure is the photovoltaic inverter failure probability, P ase _PV is a basic failure rate of the photovoltaic inverter, β4 is a fourth coefficient, t(t) is an environmental temperature, t max is a maximum allowable temperature of the inverter, β5 is a fifth coefficient, and s(t) is an inverter infection rate.

[0173] The calculation formula of the storage battery failure probability is as follows:

[0174] P e Failure=P ase _Storage×[1+β6×(1-SOC)+β7×r(t)]

[0175] In the formula, P e Failure is the storage battery failure probability, P ase _Storage is a basic failure rate of the storage battery, β6 is a sixth coefficient, SOC is a state of charge, β7 is a seventh coefficient, and r(t) is a malicious code propagation rate.

[0176] The MITRE ATT&CK framework mentioned in this embodiment is a general knowledge framework proposed by MITRE Corporation in 2013, and its Chinese name is "Countermeasures, Techniques and Common Sense". The MITRE ATT&CK framework is based on real cyber attack and defense cases and data, and adopts the TTPs (Tactics, Techniques & Procedures) methodology in military war, that is, the combination of tactics, techniques and procedures, and rearranges the network security knowledge system. The purpose is to establish a general language of network security.

[0177] High-resolution multi-drive-fault probability mapping relationship generation:

[0178] High-resolution extreme weather scenario set generation:

[0179] Meteorological modeling: adopt WRF-CFD coupling model, first generate 10km resolution regional meteorological field through weather research and forecast model WRF (Weather Research and Forecasting) model, time resolution is 1 hour. The parameter settings of WRF model include: microphysical process adopts WSM6 scheme, cumulus convection parameterization adopts Kain-Fritsch scheme, and boundary layer scheme adopts YSU scheme. Then, the output results of WRF model are used as the boundary conditions of computational fluid dynamics model CFD (Computational Fluid Dynamics) model, and the micro-terrain wind speed, wind direction and other parameters of 1km resolution are calculated through CFD model nesting, and the time resolution is 10 minutes. CFD model adopts Reynolds average Navier-Stokes equation, and k-ε model is used as turbulence model to accurately simulate the micro-terrain influence such as building flow around and mountain valley pipe effect.

[0180] Network security scenario generation: based on MITRE ATT&CK framework, construct network security attack scenario library, including attack type (such as phishing attack, malicious code injection, DDoS attack), attack path (such as terminal device→edge computing node→dispatch center), time sequence of attack intensity parameters (r(t), p(t), s(t)). Through Monte Carlo simulation, 100,000+ network security scenarios with different attack intensity and attack path are generated, and each scenario has a time span of 24 hours and a time resolution of 1 minute.

[0181] Joint scenario set construction: align the high-resolution meteorological scenario and the network security scenario according to the time axis to generate "meteorological-network security" joint scenario set, and the total number is meteorological scenario number x network security scenario number = 500 x 100,000 = 50,000,000+.

[0182] Failure probability mapping: based on the physical characteristics of the equipment and data-driven, the functional relationship between extreme weather, network security attacks and power equipment failure probability is established, and the formula of transmission line failure probability (P t Failure) is as follows:

[0183] P t Failure = P ase _ base × [1 + β1 × (v(t) / v nitical -1) + β2 × (i(t) i nitical ) + β3 × p(t)]

[0184] Where, P ase base is the basic failure rate (times / year, such as 0.02 times / year, based on historical failure statistics); v nitical is the critical wind speed of line design (such as 25 m / s, determined by tower strength); i nitical is the design critical ice thickness of line (such as 10 mm); β1, β2, β3 are coefficients, which are obtained by historical failure data regression, β1=0.05, β2=0.03, β3=0.02; v(t) is the real-time wind speed of 1 km resolution output by CFD (m / s); i(t) is the ice thickness (mm); p(t) is the data tampering probability (%) reflecting the influence of network security attacks on line state monitoring data errors, delayed maintenance or protection action. Photovoltaic inverter failure probability (P +＝ PV failure) is:

[0185] P +＝ PV failure = P ase _ PV × [1 + β4 × (t(t) t max -1) + β5 × s(t)]

[0186] Where, P ase _ PV is the basic failure rate of photovoltaic inverter (such as 0.01 times / year); t max is the maximum allowable operating temperature of the inverter (such as 60℃); β4=0.04, β5=0.03, obtained by historical data regression; t(t) is the ambient temperature (℃); s(t) is the infection rate of inverter (%) reflecting the influence of network security attacks on inverter control function. Energy storage battery failure probability (P e Failure) is:

[0187] P e Failure = P ase _ storage × [1 + β6 × (1-SOC) + β7 × r(t)]

[0188] Where, P aseThe energy storage is the basic failure rate of the energy storage battery (such as 0.015 times / year); β6=0.06, β7=0.025; SOC is the state of charge, the battery cycle life is shortened when the SOC is too low, and the failure probability is increased; r(t) is the malicious code propagation rate (times / hour), which reflects the influence of the diffusion speed of network security attacks in the energy storage system on the failure probability.

[0189] Scene reduction algorithm: the number of generated initial scene sets is huge (50 million +), and the calculation amount is huge, so the improved K-medoids algorithm is used to retain key scenes, and the specific steps are as follows:

[0190] Scene feature extraction: for each joint scene, a feature vector is extracted, including the maximum value, average value, and change rate of meteorological parameters (v(t), i(t), t(t)), the maximum value, average value, and duration of network security parameters (r(t), p(t), s(t)), and the peak value of the failure probability of each device. Initial clustering: the K-medoids algorithm is used to divide the scenes into K clusters (K=1000), and the center scene (medoid) of each cluster is the scene with the smallest distance to other scenes in the cluster. The distance measurement uses the Euclidean distance, that is, the square root of the sum of the square differences of the elements of the feature vectors of the two scenes. Objective function optimization:

[0191] minΣΣd(S i ,C j )×P(S i )

[0192] Where S i is the i-th original scene; C j is the center scene of the j-th scene cluster; d(S i ,C j ) is the distance between scene S i and center scene C j ; P(S i ) is the occurrence probability of scene S i (based on meteorological data and network security attack frequency statistics, the meteorological scene probability is calculated by the frequency of historical meteorological data, the network security scene probability is calculated by the frequency of attack type, and the joint scene probability is the product of the two); constraint condition: retain extreme scenes with risk exceeding threshold value (such as risk value>0.8), to ensure that the tail risk is not ignored; the risk value is obtained by preliminary calculation of load loss, economic loss, etc. Comprehensive evaluation; each cluster contains at least one scene with different network security attack intensity and one scene with different meteorological conditions to ensure the diversity of the scenes; after scene reduction, 1000 key scenes are retained, which cover more than 95% of the risk information in the original scene set.

[0193] Step 1: based on the meteorological data and network security data of each link at the current time and combined with the pre-constructed multi-drive-fault probability mapping relationship, the fault probability of each device under each scene is obtained, and the tripping device is determined based on the fault probability of each device under each scene, including:

[0194] The meteorological data and network security data of each link at the current time are substituted into the relationship formula of the transmission line fault probability, the relationship formula of the photovoltaic inverter fault probability and the relationship formula of the energy storage battery fault probability to obtain the transmission line fault probability, the photovoltaic inverter fault probability and the energy storage battery fault probability;

[0195] According to the fault probability, the device with the highest fault probability in the same scene can be regarded as the tripping device, so that the device is tripped and exits operation, or according to the task requirement, the device with a fault probability greater than a set threshold in different scenes can be selected as the tripping device.

[0196] Step 1 is to convert meteorological-network parameters into device fault probability, and based on the device fault probability, it is determined whether to "roll the dice" to trip the device at this time. Once tripped, the full-link chain model in step 2 is triggered to generate a new state vector X(t+Δt) at the t+Δt time. That is, at each time step, the device fault probability is calculated, and based on the device fault probability, it is determined which device to exit. Once exited, the new state vector is immediately updated through the full-link chain model.

[0197] Step 2: when the device is tripped, the meteorological data, network security data of each link at the current time, and the power grid state vector at the last time are obtained, including:

[0198] After the device is tripped, the actual active power, reference power, rated power, actual voltage, rated voltage, reliability margin, available capacity, peak load, actual charge and discharge power, rated charge and discharge power, dispatching instruction power, current price, average price, etc. of each link at the current time are obtained.

[0199] Before step 3, it also includes constructing a full-link chain model, which will be further introduced below:

[0200] The construction process of the full-link chain model includes:

[0201] Based on the historical data of the state vector of each link and combined with the deep neural network, the weight matrix of each link is obtained;

[0202] Based on the association relationship between the historical data of the state vector of each link, the initial matrix composed of the historical data of the state vector of each link is modified to obtain the coupling matrix of each link across to the next link;

[0203] A joint disturbance function is constructed based on historical data of the state vector, a function relationship between meteorological factors and the state vector in each link, and a function relationship between network security attack factors and the state vector in each link.

[0204] A full-link chain model is constructed based on the weight matrix of each link, the coupling matrix of each link, the joint disturbance function, and historical data of the state vector of each link.

[0205] The state vector of each link is calculated based on power grid operation parameters.

[0206] The construction process of the full-link chain model specifically includes:

[0207] The power generation (including centralized power plants and distributed power sources), power transmission, power distribution, and power consumption links are deeply coupled through energy flow and information flow, and a single link failure can trigger a chain reaction (such as distributed power source off-grid leading to power distribution network voltage fluctuation, and then triggering load shedding). A mathematical model needs to be established to describe the dynamic transmission law of the state variables of each link, and to quantify the chain effect of "failure-propagation-amplification". The specific steps are as follows:

[0208] Define the state variables of each link, and describe the chain process through differential equations:

[0209]

[0210] Wherein: s k (t) is the state vector of link k∈{G,T,D,E} at time t; P k (t) is the active power deviation rate, P k (t) = (P act (t) - P ref ) / P rated , P act (t) is the actual active power, P ref is the reference power, P rated is the rated power; Q k (t) is the voltage deviation (unit value), Q k (t) = V k (t) - 1.0, V k (t) is the actual voltage; R k (t) is the reliability margin, R k (t) = C avail (t) / L peak , C avail (t) is the available capacity, L peak is the peak load.

[0211] Power generation link state variable (G):

[0212] Active power deviation rate (ΔP+＝ ):

[0213]

[0214] where P +＝ actual is the actual active power of the power generation link, P +＝ reference is the planned active power, P +＝ rated is the rated active power. This variable reflects the deviation of the power generation link from the plan, the greater the deviation, the more unstable the power generation link, the more likely to cause chain failure. For example, the rated power of a certain photovoltaic power station is 100 MW, the planned active power is 80 MW, and the actual active power is 50 MW due to the sudden drop in light intensity, then ΔP +＝ = (50-80) / 100 = -0.3.

[0215] voltage deviation (ΔU +＝ ) : ΔU +＝ = U +＝ actual - U +＝ rated (unit value), where U +＝ actual is the actual voltage of the power generation node, U +＝ rated is the rated voltage (unit value 1.0). This variable reflects the voltage stability of the power generation node, and voltage deviation exceeding a certain threshold (such as ±0.05 unit value) may cause grid-connected equipment protection action.

[0216] new energy penetration rate (λ +＝ ) : where P +＝ new energy is the actual output of new energy (wind power, photovoltaic, etc.), P +＝ total is the total output of the power generation link. This variable reflects the proportion of new energy in the power generation link, the higher the penetration rate, the greater the impact of weather, and the higher the chain failure risk.

[0217] network security state (S +＝ ) : S +＝ = (1 - infection rate x data integrity), where the infection rate is the proportion of power generation link control terminals (such as wind turbine controllers, photovoltaic inverters) infected by malicious software, and the data integrity is the probability that the power generation data is not tampered with during transmission. This variable reflects the network security status of the power generation link, the lower the value, the higher the network security risk.

[0218] Transmission link state variables (T) :

[0219] active power deviation rate (ΔP t ) : ΔP t = (P t actual - P t flow) / P t rated, where Pt Actual active power of transmission line, P t Actual planned power flow, P t Rated power of transmission line. This variable reflects the deviation of power transmission of transmission section from the plan, and the deviation is too large, which will lead to overload of the line.

[0220] Voltage deviation amount (ΔU t ): ΔU t = U t Actual - U t Rated (per unit), which reflects the voltage stability of the transmission node, U t Actual voltage of transmission line, U t Rated voltage of transmission line.

[0221] Line temperature (θ t ): Actual measured temperature of transmission line, which exceeds the allowable temperature (such as 70℃) will lead to increased line resistance, reduced transmission capacity, and even cause the line to be blown.

[0222] Reliability margin (M t ): M t = (P t Rated - P t Actual) / P t Peak, where P t Peak is the historical maximum transmission power. This variable reflects the backup capacity of the transmission section, and the smaller the margin, the weaker the ability to resist faults.

[0223] Distribution section state variable (D):

[0224] Voltage deviation amount The voltage stability of the distribution section has a greater impact on user equipment, and the deviation is too large, which will lead to damage to the user equipment, Actual voltage of distribution section, Rated voltage of distribution section.

[0225] Load transfer rate Where the transferred load is the load transferred to other lines through the distribution automation system. This variable reflects the flexibility of the distribution section to deal with faults, and the higher the transfer rate, the stronger the resilience of the distribution system.

[0226] Network security state Where the smart switch infection rate is the proportion of smart switches of the distribution section infected by malicious software, and the communication link availability is the probability of normal operation of the communication link of the distribution automation system.

[0227] Load segment state variable (L):

[0228] Load deviation ratio (ΔL): ΔL = (Lactual - Lpredicted) / Lpredicted, where Lactual is the actual load and Lpredicted is the predicted load. This variable reflects the uncertainty of the load, and a large deviation can lead to a supply-demand imbalance.

[0229] Sensitive load proportion (γ l ): γ l = sensitive load / total load, where the sensitive load refers to the load sensitive to voltage and frequency fluctuations (such as hospital ICU equipment and precision manufacturing equipment). This variable reflects the vulnerability of the load, and a higher proportion means greater loss caused by chain failure.

[0230] Demand response capability (R l ): R l = reducible load / total load, where the reducible load is the load that can be temporarily reduced through demand response measures. This variable reflects the active adjustment capability of the electricity consumption segment in response to chain failure.

[0231] Energy storage segment state variable (E):

[0232] State of charge (SOC): SOC = actual storage capacity / rated storage capacity, reflecting the available capacity of the energy storage segment, and a too low or too high SOC will affect its regulation capability.

[0233] ΔP e = (P e actual - P e command) / P e rated, where P e actual is the actual charging and discharging power, P e command is the dispatch command power, P e rated is the rated charging and discharging power, and ΔP e is the power deviation. This variable reflects the accuracy of the energy storage segment in executing the dispatch command, and a large deviation will affect its support for the power grid.

[0234] Network security state (S e ): S e = (1 - energy storage controller infection rate) x data transmission delay -1 , where the energy storage controller infection rate is the proportion of energy storage system controllers infected by malicious software, and the data transmission delay is the delay time (seconds) of data transmission between the energy storage system and the dispatch center. The smaller the delay, the higher the value of S e .

[0235] Carbon segment state variable (C):

[0236] Carbon intensity deviation (ΔC): ΔC = (C actual - C quota) / C quota, where C actual is the actual carbon intensity (tCO2 / MWh), and C quota is the carbon emission quota per unit of electricity. This variable reflects the deviation of carbon emissions from the quota, with positive values indicating overage and negative values indicating surplus.

[0237] Carbon trading price fluctuation (σc): σc = (current price - average price) / average price, reflecting the price stability of the carbon trading market. The greater the fluctuation, the higher the economic risk of the carbon link.

[0238] Carbon data integrity (Ic): Ic is the probability that carbon emission data has not been tampered with. The lower the data integrity, the higher the risk of the carbon link, which may result in errors in quota calculation and economic losses.

[0239] The state vector of each link can be represented as:

[0240] X +＝ (t) = [ΔP +＝ , ΔU +＝ , λ +＝ , S +＝ ] T ;

[0241] X t (t) = [ΔP t , ΔU t , θ t , M t ] T

[0242]

[0243] X l (t) = [ΔL, γ l , R l ] T

[0244] X e (t) = [SOC, ΔP e , S e ] T

[0245] Xc(t) = [ΔC, σc, Ic] T

[0246] X +＝ (t), X t (t), X1(t), X e (t), Xc(t) are the state vectors of power generation, power distribution, power transmission, power consumption, energy storage and carbon trading, respectively. The chain propagation equation is (describing the evolution of state over time):

[0247] X i (t+Δt) = ReLU(W i X i (t) + Σ(J ij X i (t)) + F(M(t), A(t)) x Δt)

[0248] Wherein: X i (t) is the state vector of the link i at time t; W i is the weight matrix of the link i (obtained by training the deep neural network DNN, reflecting the state association in the historical data); for example, for the power generation link, the element W i in W +＝11 represents the influence weight of ΔP +＝ on the next time state of itself, which can be obtained by training the association relationship between the value of ΔP +＝ at t time and the value at t+Δt time, W +＝11 is the first row and first column element of the power generation link weight matrix W +＝ , ΔP +＝ is the active power deviation rate of the power generation link at a certain time; J ij is the coupling matrix from link j to link i (quantifying the influence strength across links); for example, J t+＝ represents the coupling matrix of the power generation link to the power transmission link, wherein the element J t+＝11 represents the influence strength of the ΔP +＝ of the power generation link on the ΔP t of the power transmission link, the larger the value, the more significant the influence of the power generation power deviation on the power transmission power deviation, ΔP t is the power deviation of the power transmission link; the determination of the coupling matrix needs to be combined with the physical mechanism and the historical data, such as the voltage support of the power transmission link to the power distribution link, which can be calculated based on the circuit theory to obtain the basic value, and then combined with the historical voltage data for correction; M(t) is a meteorological disturbance vector, M(t) = [v(t), h(t), t(t), i(t)] T , wherein v(t) is the wind speed (m / s), h(t) is the relative humidity (%), t(t) is the temperature (°C), and i(t) is the ice thickness (mm); A(t) is a network security attack vector, A(t) = [r(t), p(t), s(t)] t , wherein r(t) is the malicious code propagation rate (times / hour), p(t) is the data tampering probability (%), and s(t) is the terminal device infection rate (%); F(M(t), A(t)) is a joint disturbance function of meteorology and network security attack, which is obtained by fusing the influence of meteorological factors and network security factors on the state of each link; for example, for the power transmission link θ tThe temperature component in F(M(t), A(t)) directly affects θ t , while cyber-attacks can cause temperature monitoring data tampering, indirectly affecting the calculation of θ t ; the ReLU activation function (ReLU(x) = max(0, x)) introduces nonlinearity, simulating threshold characteristics such as device saturation, protective actions, etc., where x is the independent variable (intermediate calculation result after linear combination); for example, when the power deviation of the transmission line exceeds a certain threshold, the protective device acts, at which point the ReLU function can truncate the part that exceeds the threshold, simulating the state change after the protective action; Δt is the time step (default 1 minute, which can be adjusted according to the dynamic characteristics of different links, such as voltage and frequency-related state variables, which can use a time step of 1 second).

[0249] wherein, the training method of the weight matrix (W i ) specifically includes:

[0250] Data collection: Collect historical data of state variables of each link in the past 5 years, with a sampling interval consistent with the time step Δt (e.g. 1 minute), the data should cover different seasons, different weather conditions, and different network security states, with a total data volume of not less than 1 million, data preprocessing: missing data is supplemented by interpolation method, abnormal data (such as values obviously exceeding the reasonable range) is removed or corrected, the data is divided into training set (70%), validation set (20%) and test set (10%), deep neural network (DNN) structure: DNN with 3 layers of hidden layers, the number of input layer nodes is the dimension of the link state vector (e.g. 4 nodes for the power generation link), the number of hidden layer nodes is 64, 32, and 16 respectively, the number of output layer nodes is the same as the input layer, the activation function uses ReLU, training process: the state vector X i (t) at time t is input, the state vector X i (t+Δt) at time t+Δt is output, the Adam optimizer is used, the loss function is mean square error (MSE), and the training is iterated until the loss function of the validation set no longer decreases, obtaining the weight matrix W i .

[0251] The determination method of the coupling matrix (J ij ) specifically includes:

[0252] Physical mechanism analysis: According to the energy flow, information flow and carbon flow transmission law of the power system, the influence relationship and direction between each link are determined. For example, the power generation link has a direct power impact on the power transmission link, the power transmission link has a voltage support impact on the power distribution link, and the load change of the power consumption link will affect the power generation and energy storage link, etc.; Initial matrix construction: Based on physical formulas and empirical data, the initial value of the coupling matrix is constructed. For example, the power coupling coefficient from the power generation link to the power transmission link can be initially determined according to the power flow calculation formula, that is, J t+＝ The initial value of the influence coefficient of ΔP +＝ on ΔP t is the impedance parameter related value of the transmission line; Historical data correction: The initial matrix is corrected by using the correlation of the state variables of each link in the historical data. The least square method is adopted, and the historical data of the state variables of each link are used as samples, so that the deviation between the cross-link influence calculated according to the coupling matrix and the actual historical data is minimized, thereby obtaining the final coupling matrix J ij .

[0253] The construction method of the joint disturbance function (F(M(t), A(t))) includes:

[0254] The meteorological factor influence sub-function F m (M(t)): For different links, the functional relationship between meteorological factors and state variables is established. For example, for the line temperature θ t of the power transmission link, the temperature component in F m (M(t)) is t(t), the wind speed v(t) affects the line heat dissipation efficiency to affect θ t , which can be expressed as the change amount of θ t is inversely proportional to v(t); the ice thickness i(t) increases the weight of the line, affects the mechanical stress of the line, and then indirectly affects ΔP t , which can be expressed as the change amount of ΔP t is proportional to i(t); the network security attack influence sub-function F a (A(t)): For different links, the functional relationship between network security attack factors and state variables is established. For example, for S e (network security state) of the energy storage link, the higher r(t) (malicious code propagation rate) in F a (A(t)) is, the faster the energy storage controller infection rate rises, and the more obvious S e decreases; the higher p(t) (data tampering probability) is, the larger the data transmission delay is, and the lower S e value is; Joint function fusion: F(M(t), A(t)) = α × F m (M(t)) + (1-α) × F a(A(t)) where a is a weight coefficient determined according to the contribution degree of meteorological factors and network security attacks in historical failures. For example, a has a larger value (such as 0.7) in a failure case mainly caused by extreme weather, and a has a smaller value (such as 0.3) in a failure case mainly caused by network security attacks.

[0255] Step 3: Based on the power grid state vector of the last time and the pre-constructed full-link chain model, the state vector of each link at the current time is obtained, including:

[0256] The power grid state vector of each link at the last time is substituted into the pre-constructed full-link chain model, and the state vector at the current time is updated by the full-link chain model based on the power grid state vector of each link at the last time, that is, the power grid state vector X(t) at the last time and the external disturbance M(t), A(t) are read in every time interval Δt, and the current power grid state vector X(t+Δt) is output through the coupling matrix and the ReLU activation function.

[0257] Step 4: Based on the state vector of each link at the current time, the total risk of the system is calculated by the analytic hierarchy process, including:

[0258] Based on the power shortage loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission loss, network security risk degree, and dynamic risk calculation formula, the dynamic risk value of each node at each time in each scene is calculated;

[0259] Based on the dynamic risk value of each node at each time in each scene, the accumulated risk contribution degree of each node in the chain process of each scene is calculated by the accumulated risk contribution degree calculation formula;

[0260] Based on the accumulated risk contribution degree and the total risk calculation formula, the total risk of the system is calculated;

[0261] The state vector of each link includes: power shortage loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission loss, and network security risk degree.

[0262] Step 4 specifically includes:

[0263] Node dynamic risk quantification:

[0264] Dynamic risk index (DRI i,j (t)): Comprehensive power shortage loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission loss, and network security risk, quantifying the dynamic risk value (0-1) of node i at time t in scene j:

[0265] DRI i,j (t) = ω1×L i,j (t) + ω2×U i,j(t) + ω3×F i,j (t) + ω4×E i,j (t) + ω5×C i,j (t) + ω6×S i,j (t)

[0266] where, L i,j (t) is the relative lack of power supply of node i in scenario j at time t, L i,j (t) = lack of power supply i,j (t) / total system load j , lack of power supply i,j (t) is the total lack of power of node i at the scene and time, total system load j is the total system load under scenario j; U i,j (t) is the voltage out-of-limit degree, U i,j (t) = max(0, |ΔU i,j (t) | - 0.05) / 0.2, ΔU i,j (t) is the voltage deviation of node i (unit value), when the deviation is within ±0.05, U i,j (t) = 0, when it exceeds 0.2, U i,j (t) = 1; F i,j (t) is the frequency out-of-limit degree, F i,j (t) = max(0, |Δf i,j (t) | - 0.2) / 0.5, Δf i,j (t) is the frequency deviation (Hz), the rated frequency is 50 Hz, when the deviation is within ±0.2 Hz, F i,j (t) = 0, when it exceeds 0.5 Hz, F i,j (t) = 1; E i,j (t) is the lack of energy storage reserve capacity degree, E i,j (t) = max(0, 0.2 - SOC i,j (t)) / 0.2, when SOC i,j (t) ≥ 0.2, E i,j (t) = 0, when SOC i,j (t) = 0, E i,j (t) = 1, SOC i,j (t) is the state of charge of energy storage of node i in scenario j at time t; C i,j (t) is the loss degree of carbon emission exceeding the standard, C i,j (t) = max(0, ΔC i,j (t)) / 0.5, ΔC i,j (t) is the carbon emission intensity deviation, when ΔC i,j (t) ≤ 0, C i,j (t) = 0, ΔC i,jC i,j (t) = 0.5; S i,j (t) = 1; S i,j (t) = 1; S i,j (t) = 1; S i,j (t) = 1; S

[0267] AHP method: invite 10 experts in the field of power system, network security, carbon emission, construct judgment matrix, determine the subjective weight of each index: ω1_AHP = 0.3, ω2_AHP = 0.2, ω3_AHP = 0.15, ω4_AHP = 0.1, ω5_AHP = 0.1, ω6_AHP = 0.15; Entropy weight method: based on the historical data of 1000 key scenes, calculate the information entropy of each index, the lower the entropy value, the higher the index differentiation, the greater the weight. The objective weight is calculated: ω1_entropy = 0.25, ω2_entropy = 0.22, ω3_entropy = 0.18, ω4_entropy = 0.12, ω5_entropy = 0.1, ω6_entropy = 0.13; Joint weight: ω k = 0.5 × ω k _AHP + 0.5 × ω k _entropy, finally ω1 = 0.275, ω2 = 0.21, ω3 = 0.165, ω4 = 0.11, ω5 = 0.1, ω6 = 0.14.

[0268] Chain vulnerability index (CVI i,j ):

[0269] The formula for quantifying the cumulative risk contribution of node i in the chain launch process in scene j is:

[0270] CVI i,j = ∫(t0 to t1)DRI i,j (t) × e^(-λ × (t-t0))dt

[0271] Where t0 is the initial time of failure; t1 is the first time the system exceeds the limit (such as the time when voltage collapses or frequency instability occurs); λ is the decay coefficient, which is dynamically adjusted with time scale: λ = 0.3 for seconds (recent risk weight is very high), λ = 0.1 for minutes, λ = 0.05 for hours; e^(-λ × (t-t0)) is the decay factor, which makes the recent risk contribution weight greater. If the initial failure is a network attack, CVI i,j is multiplied by a factor of 1.2 (because the network attack chain launches faster).

[0272] System-level risk aggregation and emergency decision output:

[0273] System total risk (SR j ): the system total risk under scenario j is calculated by integrating the chain vulnerability index CVI of each node, the importance of load, and the impact of carbon emission:

[0274] SR j =∑(CVI i,j ×α i ×β i )

[0275] wherein, α i is the weight of sensitive load supplied by node i, α i = sensitive load i / system total sensitive load, sensitive load i is the power of sensitive load (such as hospital, data center) within the power supply range of node i; β i is the load level coefficient (1-5 levels, 5 levels are the most important), 5 level load (hospital ICU) β i = 5, 4 level (data center) β i = 4, 3 level (commercial building) β i = 3, 2 level (residential life) β i = 2, 1 level (interruptible industrial load) β i = 1; the increment of carbon emission is added to the correction: SR j ' = SR j × (1+γ×ΔC total,j ), wherein, SR j ' is the system total risk under scenario j after adding the increment of carbon emission correction, ΔC total,j is the total carbon emission over-standard amount (tCO2) under scenario j, and γ = 0.001 is the carbon emission impact coefficient.

[0276] The city power grid chain risk quantification method disclosed in the application further comprises:

[0277] According to the system total risk from large to small, the scenario corresponding to the first set percentage of system total risk is a high-risk scenario;

[0278] According to the cumulative risk contribution degree from large to small, the node corresponding to the second set percentage of cumulative risk contribution degree is a high-risk node.

[0279] The key risk node identification specifically comprises:

[0280] For each scenario j, the CVI of each node is calculated according to the following formula: i,jSort from large to small, take the top percentage (take 20% as an example) of nodes as the key risk nodes in this scenario. By comparing the key risk nodes in different scenarios, the "core key nodes" that are key risk nodes in multiple scenarios are summarized, such as power transmission line hub nodes, energy storage scheduling center nodes, carbon trading center nodes, etc.

[0281] Risk ranking and early warning, specifically including:

[0282] Sort the SR of 1000 key scenarios j Take the top 10% of scenarios as high-risk scenarios to generate a risk warning list, including scenario occurrence probability, main impact nodes, expected loss, etc.

[0283] Emergency decision output:

[0284] Risk heat map: Real-time display of each node's chain vulnerability index CVI value, marked with red (CVI>0.8), yellow (0.5<CVI≤0.8), and green (CVI≤0.5) three colors, highlighting high contribution nodes.

[0285] Targeted control strategy specifically includes:

[0286] For high-CVI new energy nodes, it is recommended to start energy storage backup 1 hour in advance, and increase rotating backup capacity by 5%-10%; for high-CVI power transmission line nodes, according to the extreme weather predicted by WRF-CFD, arrange line inspection 24 hours in advance, reinforce the tower, and remove the icing hazard; for high-CVI network security nodes (such as dispatching center, smart meter), strengthen firewall protection, improve data encryption level, scan terminal equipment for malicious software, and start isolation measures when the infection rate exceeds 5%; for high-CVI carbon link nodes, when ΔC approaches the quota threshold, adjust the power generation plan in advance, increase low-carbon power output, reduce high-carbon power use, or purchase carbon quota reserves, and ΔC is the deviation of node carbon intensity.

[0287] Emergency resource allocation suggestion: According to the CVI value and geographical location of each node, optimize the deployment location and quantity of emergency repair teams, backup power, and energy storage devices to ensure that resources can reach key risk nodes within 30 minutes when high-risk scenarios occur.

[0288] Through the above method, the present application can comprehensively quantify the chain risk of the "source-grid-load-storage-carbon" whole chain in the urban new power system under extreme weather and network security attacks, provide accurate quantitative support for system resilience optimization, emergency dispatch and low-carbon transformation, and effectively improve the safety guarantee ability of urban power system.

[0289] The present application has the following advantages:

[0290] 1) Break through the limitations of traditional single-link static assessment:

[0291] Incorporate the "source-network-load-storage-carbon" whole chain into chain risk assessment, quantify the dynamic coupling effect across links and time scales (such as new energy off-grid → carbon quota imbalance → storage scheduling disorder), and solve the problem that existing technologies cannot capture multi-link chain failure.

[0292] 2) Improve risk assessment accuracy and timeliness:

[0293] Through high-resolution weather modeling (1km resolution), network security attack coupling analysis and dynamic risk indicators, the device failure probability prediction error is reduced by more than 30%, the risk warning response speed is improved by 50%, and key risk nodes (such as power transmission lines under extreme weather and smart meters under network attacks) can be identified in advance.

[0294] 3) Support safe-low-carbon collaborative decision-making:

[0295] Integrate carbon emission loss, carbon trading fluctuations and other low-carbon factors to quantify the "safety-low-carbon" coupling risk and provide targeted strategies for emergency dispatch (such as preferentially starting low-carbon backup power and dynamically adjusting carbon quotas), taking into account power system safety and "double carbon" goals.

[0296] Based on the above, the present application mainly includes the following three parts:

[0297] (1) Establish a mathematical model of whole-link chain propagation of power generation and power consumption: There is a complex exchange of energy flow and information flow between the power generation link (including centralized power plants and distributed power sources), the power transmission link, the power distribution link and the power consumption link. When a fault occurs in a link, the fault will propagate in various forms between links. A mathematical model is needed to accurately describe this whole-link chain propagation process, and to clearly define the relationship between the state variables of each link and the rules of fault propagation.

[0298] (2) Generate high-precision extreme weather-device failure probability mapping relationship: There is a close relationship between extreme weather (such as different wind speeds, ice thickness, humidity, etc.) and power equipment (such as power transmission lines, transformers, generators, etc.) failure probability. Due to the complexity of extreme weather and the diversity of equipment characteristics, a method is needed to generate a high-precision extreme weather-equipment failure probability mapping relationship to improve the accuracy of risk assessment. This requires high-precision meteorological data to be obtained, and a model that accurately reflects the internal relationship between meteorological factors and equipment failure probability to be established.

[0299] (3) Quantify the dynamic risk contribution of each node in the chain process: In the chain process of power system fault, different nodes (such as power generation node, power transmission node, power distribution node, power consumption node, etc.) have different contribution to the overall risk of the system. An effective method is needed to quantify the dynamic risk contribution of each node in the chain process, and to determine which nodes are the key risk points in the system, so as to take targeted measures in emergency decision-making.

[0300] The relationship between the three is: "rolling by time step, embedding each other" microcirculation closed loop:

[0301] (1) is a "state engine", which reads in the power grid state vector X(t) and external disturbance M(t), A(t) at the last moment in each moment Δt, and outputs X(t+Δt) at the current moment through the coupling matrix and ReLU evolution equation; (2) is an "external disturbance translator", which immediately maps high-resolution meteorological elements (wind speed, icing, humidity, etc.) and network security indicators (infection rate, data tampering probability, etc.) into the conditional failure probability of each device at the same moment, providing the real-time vulnerability of "how easy is the device to break" for risk calculation; (3) is a "risk accountant", which simultaneously substitutes the power grid state given by (1) and the vulnerability given by (2) into the dynamic risk index DRI, weights the six dimensions of power shortage, voltage, frequency, energy storage, carbon emission and network security, integrates to get node-level CVI and system-level SR, and sorts in the massive "meteorological-network" joint scene, directly filters out the top first set percentage (such as 10%) of high-risk scenarios, and returns the SR result to correct the weight matrix of (1), realizing the closed loop of "pushing state-translating vulnerability-calculating risk-revising".

[0302] The failure probability of each device in each scenario output by (2) above is the premise, that is, in each time step, the failure probability of each device in each scenario is randomly determined to determine which device to exit; once the device exits, the chain equation of (1) immediately updates the new power flow, voltage, frequency, energy storage SOC, carbon emission and network security state; these updated quantities are sent into the DRI formula to get the L, U, F values at this step. Therefore, the higher the failure probability of each device in each scenario, the easier the device is to break, and the more likely the DRI dimensions are to jump, and the CVI and SR after integration are amplified. The failure probability is not directly multiplied into the DRI, but indirectly but decisively shapes the total risk of the system through "controlling the scene trajectory".

[0303] Whole-link chain propagation model construction technology: define the multi-link state variables (including energy flow: active power / voltage deviation, reliability margin; information flow: network security state, data integrity; carbon flow: carbon emission intensity deviation, carbon quota balance) of "source-network-load-storage-carbon", describe the chain propagation law through differential equation, and quantify the "fault-propagation-amplification" chain effect. Innovatively introduce weight matrix (DNN training) and coupling matrix (physical mechanism + historical data correction) to dynamically depict the cross-link influence strength (such as voltage support of transmission to distribution, influence of carbon trading center on energy storage scheduling).

[0304] High-precision multi-drive-fault probability mapping technology: generate 1km resolution micro-topography meteorological data using WRF-CFD coupling model, combine network security attack scenarios (malicious code propagation rate, data tampering probability), and construct a three-dimensional mapping relationship of "meteorology-network security-equipment failure". Propose a multi-factor failure probability formula (such as the fusion of wind speed, icing, and data tampering probability for transmission line failure probability; the association of temperature and terminal infection rate for photovoltaic inverter failure probability), breaking through the limitations of traditional single-factor mapping.

[0305] Dynamic risk quantification and key node identification technology: construct a dynamic risk index (DRI) that integrates power loss, voltage / frequency stability, energy storage reserve, carbon emission exceeding, and network security risk, determines the weight through AHP-entropy weight method, and realizes risk value quantification in the range of 0-1. Propose a chain propagation vulnerability index (CVI) to quantify the cumulative risk contribution of nodes in the chain propagation process through a time decay coefficient (second-level / minute-level / hour-level dynamic adjustment), and accurately identify core key nodes (top 20% high contribution nodes).

[0306] System-level risk aggregation and emergency decision-making technology: integrate node risk contribution, sensitive load level, and carbon emission increment to calculate the system total risk (SR), and generate a risk heat map and a hierarchical warning list. Output targeted emergency strategies (such as starting energy storage in advance for high-risk new energy nodes and dynamically adjusting quotas for carbon link nodes), ensure that resources reach key nodes within 30 minutes, and improve emergency response efficiency.

[0307] Embodiment 2:

[0308] The application also provides a city power grid chain propagation risk quantification system, comprising:

[0309] A fault probability evaluation module is configured to obtain the failure probability of each device under each scenario based on the meteorological data and the network security data of each link at the current time and the multi-drive-fault probability mapping relationship, and determine the tripped device based on the failure probability of each device under each scenario.

[0310] The parameter acquisition module is configured to acquire meteorological data, network security data, and a power grid state vector of a previous time point when the device trips;

[0311] The chain generation module is configured to obtain a state vector of a current time point of each link based on the power grid state vector of the previous time point and a pre-constructed full-link chain generation model;

[0312] The risk quantification module is configured to calculate a total risk of the system based on the state vector of the current time point of each link by using an analytic hierarchy process;

[0313] The full-link chain generation model is constructed based on a state correlation relationship in historical data of power grid operation parameters, a coupling relationship between links, and influences of meteorological factors and network security factors on states of the links, and the model construction module is configured to:

[0314] calculate a state vector of each link based on the power grid operation parameters;

[0315] train each link weight matrix based on historical data of the state vector of each link and a deep neural network;

[0316] correct an initial matrix composed of the state vector historical data of each link based on a correlation relationship between the state vector historical data of each link, to obtain a coupling matrix of each link across to a next link;

[0317] construct a joint disturbance function based on a function relationship between the state vector historical data and the meteorological factors in each link and a function relationship between the state vector historical data and the network security attack factors in each link;

[0318] construct the full-link chain generation model based on the link weight matrix, the coupling matrix of each link, the joint disturbance function, and the state vector historical data of each link.

[0319] Optionally, the full-link chain generation model is as follows:

[0320] X i (t+Δt)=ReLU(W i ×X i (t)+Σ(J ij ×X j (t))+F(M(t),A(t))×Δt)

[0321] In the formula, X i (t+Δt) is a state vector of link i at time t+Δt, X i (t) is a state vector of link i at time t, X j (t) is a state vector of link j at time t, W iis the weight matrix of link i, J ij is the coupling matrix from link j to link i, F(M(t),A(t)) is a joint disturbance function of weather and cyber-attack, and Δt is a time step.

[0322] Optionally, the step of constructing the joint disturbance function based on the historical data of the state vectors in the model construction module in combination with the function relationship between the weather factors and the state vectors in each link and the function relationship between the cyber-attack factors and the state vectors in each link comprises:

[0323] calculating a weather factor influence sub-function and a cyber-attack factor influence sub-function based on the historical data of the state vectors in combination with the function relationship between the weather factors and the state vectors in each link and the function relationship between the cyber-attack factors and the state vectors in each link;

[0324] determining the weights of the weather factor influence sub-function and the cyber-attack factor influence sub-function according to the contribution degrees of the weather factors and the cyber-attack in the historical faults;

[0325] performing weighted fusion based on the weather factor influence sub-function and the cyber-attack factor influence sub-function and the corresponding weights to obtain the joint disturbance function.

[0326] Optionally, the calculation formula of the joint disturbance function is as follows:

[0327] F(M(t),A(t))=α×F m (M(t))+(1-α)×F a (A(t))

[0328] In the formula, F(M(t),A(t)) is the joint disturbance function, α is a weight coefficient, F m (M(t)) is the function relationship between the weather factors and the state variables, and F a (A(t)) is the function relationship between the cyber-attack factors and the state variables.

[0329] Optionally, the step of training based on the historical data of the state vectors of each link in combination with the deep neural network in the model construction module to obtain the weight matrix of each link comprises:

[0330] obtaining the historical data of the state vectors of each link, the sampling interval of the historical data of the state vectors of each link being Δt, and the historical data of the state vectors of each link including operating data under different seasons, different weather conditions and different cyber-attack states;

[0331] performing outlier rejection or correction and missing value supplement processing on the historical data of the state vectors of each link to obtain processed data;

[0332] The processed data constitutes a sample set, each sample including state vectors of the same link at time t and t+Δt, t being a time point;

[0333] The sample set is divided into a training set, a verification set and a test set according to a set proportion;

[0334] The dimensions of the state vectors of the links are used to define the number of input layer nodes and the number of output layer nodes of the deep neural network; the state vector at time t is used as the input of the deep neural network, and the state vector at time t+Δt is used as the output;

[0335] The deep neural network is trained, verified and tested based on the training set, the verification set and the test set, and the weight matrix of each link is obtained.

[0336] Optionally, the risk quantification module is specifically configured to:

[0337] Based on the power loss, voltage stability, frequency stability, energy storage reserve capacity, carbon emission exceeding loss and network security risk degree, a dynamic risk value of each node at each time point in each scenario is calculated by using a dynamic risk calculation formula;

[0338] Based on the dynamic risk value of each node at each time point in each scenario, an accumulated risk contribution degree of each node in the chain process of each scenario is calculated by using an accumulated risk contribution calculation formula;

[0339] Based on the accumulated risk contribution degree, a total risk of the system is calculated by using a total risk calculation formula;

[0340] The state vector of each link includes the power loss, the voltage stability, the frequency stability, the energy storage reserve capacity, the carbon emission exceeding loss and the network security risk degree.

[0341] Optionally, the dynamic risk calculation formula is as follows:

[0342] DRI i,j (t)=ω1×L i,j (t)+ω2×U i,j (t)+ω3×F i,j (t)+ω4×E i,j (t)+ω5×C i,j (t)+ω6×S i,j (t)

[0343] In the formula, DRI i,j (t) is the dynamic risk value of node i at time t in scenario j, L i,j (t) is the relative power supply of node i at time t in scenario j, U i,j (t) is the voltage overrun degree, F i,j (t) is the frequency overrun degree, E i,j(t) is the degree of insufficient energy storage backup capability, C i,j (t) is the degree of carbon emission exceeding the standard, S i,j (t) is the degree of network security risk, ω1 is the weight coefficient of the relative power supply shortage of node i in scenario j and at time t, ω2 is the weight coefficient of the voltage out-of-limit degree, ω3 is the weight coefficient of the frequency out-of-limit degree, ω4 is the weight coefficient of the insufficient energy storage backup capability, ω5 is the weight coefficient of the carbon emission exceeding the standard, ω6 is the weight coefficient of the network security risk, i is the node number, j is the scenario number, and t is the time value.

[0344] Optionally, the accumulated risk contribution degree calculation formula is as follows:

[0345] CVI i,j =∫(t0 to t1)DRI i,j (t)×e^(-λ×(t-t0))dt

[0346] In the formula, CVI i,j is the accumulated risk contribution degree of node i in the chain process of scenario j, t0 is the initial time of failure, t1 is the first time of system out-of-limit, and λ is a decay coefficient dynamically adjusted with the time scale.

[0347] Optionally, the total risk calculation formula is as follows:

[0348] SR j =Σ(CVI i,j ×α i ×β i )

[0349] In the formula, SR j is the total risk of the system, α i is the sensitive load weight supplied by node i, and β i is the load level coefficient.

[0350] Optionally, it further includes a mapping relationship construction module, configured to:

[0351] generate each extreme weather scenario based on the current time weather data using the WRF-CFD coupling model;

[0352] generate each network security scenario based on the current time network security data using the MITRE ATT&CK framework;

[0353] align each extreme weather scenario and each network security scenario according to the time axis to generate an extreme weather-network security joint scenario set;

[0354] reduce the extreme weather-network security joint scenario set using an improved K-medoids algorithm to obtain key scenarios;

[0355] establish a function relationship of extreme weather, cyber-attack and power equipment failure probability based on the physical characteristics and data driven of the equipment in the key scene;

[0356] The function relationship of extreme weather, cyber-attack and power equipment failure probability includes a relationship formula of transmission line failure probability, a relationship formula of photovoltaic inverter failure probability and a relationship formula of energy storage battery failure probability.

[0357] Optionally, the calculation formula of the transmission line failure probability is as follows:

[0358] P t Failure = P ase _ base × [1 + β1 × (v(t) / v nitical -1) + β2 × (i(t) / i nitical ) + β3 × p(t)]

[0359] In the formula, P t Failure is the transmission line failure probability, P ase _base is the basic failure rate, β1 is the first coefficient, v(t) is the 1km resolution real-time wind speed output by CFD, v nitical is the line design critical wind speed, β2 is the second coefficient, i(t) is the ice thickness, i nitical is the line design critical ice thickness, β3 is the third coefficient, and p(t) is the data tampering probability.

[0360] The calculation formula of the photovoltaic inverter failure probability is as follows:

[0361] P +＝ PV Failure = P ase _basePV × [1 + β4 × (t(t) / t max -1) + β5 × s(t)]

[0362] In the formula, P +＝ PV Failure is the photovoltaic inverter failure probability, P ase _basePV is the basic failure rate of the photovoltaic inverter, β4 is the fourth coefficient, t(t) is the environmental temperature, t max is the maximum allowable temperature of the inverter, β5 is the fifth coefficient, and s(t) is the infection rate of the inverter.

[0363] The calculation formula of the energy storage battery failure probability is as follows:

[0364] P e Failure = P ase _baseEnergy × [1 + β6 × (1-SOC) + β7 × r(t)]

[0365] In the formula, P e Failure is the energy storage battery failure probability, Pase Energy storage is the energy storage battery basic failure rate, β6 is the sixth coefficient, SOC is the state of charge, β7 is the seventh coefficient, r(t) is the malicious code propagation rate.

[0366] Embodiment 3

[0367] As Figure 2 shown, the present application also provides an electronic device, which can be a computer device, a single-chip microcomputer device, a smart mobile device, etc. The electronic device in this embodiment can include a processor, a memory, a transceiver component, etc. The memory, the processor and the transceiver component are connected through a bus; the memory can be used to store an execution program, and the exemplary execution program can include instructions; the processor is used to execute the instructions stored in the memory. The memory can also be used to store data, which can be called and / or modified when the instructions are executed.

[0368] The processor can be a central processing unit (CPU), and can also be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc., which are the computing core and control core of the terminal, and are suitable for implementing one or more instructions, and are specifically suitable for loading and executing one or more instructions in the storage medium to implement a corresponding method flow or a corresponding function, to implement the steps of the city power grid chain risk quantification method in the above embodiment.

[0369] Embodiment 4

[0370] Based on the same inventive concept, the application further provides a readable storage medium, specifically, an electronic device readable storage medium (Memory). The electronic device readable storage medium is a memory device in the electronic device, and is used to store programs and data. It can be understood that the storage medium herein can include a built-in storage medium in the electronic device, and of course can also include an extended storage medium supported by the electronic device. The storage medium provides a storage space, and the storage space stores an operating system of the terminal. In addition, one or more instructions suitable for being loaded and executed by the processor are also stored in the storage space, and the instructions can be one or more execution programs (including program codes). It should be noted that the storage medium herein can be a high-speed RAM memory or a non-volatile memory such as at least one disk memory. Loading and executing one or more instructions stored in the storage medium by the processor can realize the steps of the urban power grid chain risk quantification method in the above embodiment.

[0371] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage, etc.) containing computer-usable program code.

[0372] The present application is described with reference to flowcharts and / or block diagrams according to the methods, devices (systems), and computer program products of the embodiments of the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, and the combination of flows and / or blocks in the flowcharts and / or block diagrams can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing apparatus to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing apparatus produce a device that implements the functions specified in the flowcharts and / or block diagrams. Figure 1 one or more flows and / or blocks Figure 1 means for performing the functions specified in the flowchart

[0373] These computer program instructions can also be stored in a computer-readable memory capable of guiding a computer or other programmable data processing apparatus to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured product including instruction means, which implements the functions specified in the flowcharts and / or block diagrams. Figure 1 one or more flows and / or blocks Figure 1 means for performing the functions specified in the flowchart

[0374] These computer program instructions can also be loaded into computer or other programmable data processing devices, so that a series of operation steps are performed on the computer or other programmable data processing devices to generate computer-implemented processes, thus the instructions executed on the computer or other programmable data processing devices provide a process for implementing the functions specified in the flowcharts Figure 1 one flow or multiple flows and / or one block or multiple blocks. Figure 1 the steps of the functions specified in the flowcharts

[0375] The above only is the embodiment of the present application, and is not used to limit the present application, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application, are included in the claim scope of the application to be approved of the present application.

Claims

1. A method for quantifying the risk of power grid chain-like events, characterized in that, include: Based on the meteorological data and network security data at the current moment of each link, combined with the pre-constructed multi-drive-fault probability mapping relationship, the failure probability of each device in each scenario is obtained, and the tripping device is determined based on the failure probability of each device in each scenario. When the equipment trips, it acquires the current meteorological data, network security data, and the power grid state vector of the previous moment for each component. Based on the power grid state vector of the previous moment and combined with the pre-constructed full-link chain-generation model, the state vector of each link at the current moment is obtained. The total risk of the system is calculated using the analytic hierarchy process based on the current state vector of each component. The full-chain generation model is constructed based on the state correlations in the historical data of the power grid state vector, the coupling relationships between each link, and the impact of meteorological and cybersecurity factors on the state of each link.

2. The method according to claim 1, characterized in that, The construction process of the full-chain transmission model includes: The weight matrix of each stage is obtained by training a deep neural network based on historical data of the state vectors of each stage. Based on the correlation between the historical state vector data of each stage, the initial matrix composed of the historical state vector data of each stage is corrected to obtain the coupling matrix from each stage to the next stage. Based on historical state vector data, a joint perturbation function is constructed by combining the functional relationships between meteorological factors and state vectors in each stage and the functional relationships between network security attack factors and state vectors in each stage. A full-chain transmission model is constructed based on the weight matrix, coupling matrix, joint perturbation function, and historical state vector data of each stage.

3. The method according to claim 2, characterized in that, The joint perturbation function, constructed based on historical state vector data and combining the functional relationships between meteorological factors and state vectors at each stage, as well as the functional relationships between cybersecurity attack factors and state vectors at each stage, includes: Based on historical state vector data, and combining the functional relationships between meteorological factors and state vectors in each stage, as well as the functional relationships between cybersecurity attack factors and state vectors in each stage, the influence sub-functions of meteorological factors and cybersecurity attack factors are calculated. The weights of the meteorological factor influence sub-function and the network security attack factor influence sub-function are determined based on the contribution of meteorological factors and network security attack factors in historical failures. The joint perturbation function is obtained by weighting and fusing the sub-functions of meteorological factors and network security attack factors, along with their corresponding weights.

4. The method according to claim 2, characterized in that, The formula for calculating the joint perturbation function is as follows: F(M(t),A(t))=α×F m (M(t))+(1-α)×F a (A(t)) In the formula, F(M(t),A(t)) is the joint perturbation function, α is the weighting coefficient, and F m (M(t)) represents the functional relationship between meteorological factors and state variables, F a (A(t)) represents the functional relationship between network security attack factors and state variables.

5. The method according to claim 2, characterized in that, The full-chain transmission model is shown in the following formula: X i (t+Δt)=ReLU(W i ×X i (t)+Σ(J ij ×X i (t))+F(M(t),A(t))×Δt) In the formula, X i (t+Δt) is the state vector of element i at time t+Δt, X i (t) is the state vector of element i at time t, X j (t) is the state vector of element j at time t, W i Let J be the weight matrix of link i. ij Let F(M(t),A(t)) be the coupling matrix from stage j to stage i, and let Δt be the time step.

6. The method as described in claim 2, characterized in that, The weight matrix for each stage is obtained by training a deep neural network based on historical data of the state vectors of each stage, including: The historical state vector data of each stage is obtained, and the sampling interval of the historical state vector data of each stage is Δt, including: operation data under different seasons, different weather conditions, and different network security states; The historical state vector data of each stage are processed by removing or correcting outliers and supplementing missing values ​​to obtain the processed data. The processed data are used to form a sample set. Each sample includes the state vector of the same link at time t and t+Δt, where t is the time. The sample set is divided into a training set, a validation set, and a test set according to a set ratio; The number of input layer nodes and output layer nodes of the deep neural network are defined according to the dimension of the state vector of each stage; the state vector at time t is used as the input of the deep neural network, and the state vector at time t+Δt is used as the output; The deep neural network is trained, validated, and tested based on the training set, validation set, and test set to obtain the weight matrix of each stage.

7. The method as described in claim 6, characterized in that, The process of training, validating, and testing the deep neural network based on the training set, validation set, and test set to obtain the weight matrix for each stage includes: The deep neural network is trained by forward propagation using the training set, and the training loss is calculated using the mean squared error as the loss function. The gradient is calculated through backpropagation, and the weight matrix W of the deep neural network is updated using the Adam optimizer. i ; After every K training rounds, the validation set is used to evaluate the validation loss value of the current deep neural network model, where K is a preset period. If the change in the validation loss value over N consecutive periods is less than a set threshold, then training is stopped, and a validated deep neural network model is obtained, where N is a preset positive integer. The test set is used to evaluate the generalization performance of the validated deep neural network model, and the weight matrix W output by the deep neural network model whose generalization performance meets the set threshold requirement is used. i This serves as the weight matrix for each stage.

8. The method as described in claim 1, characterized in that, The calculation of the total system risk based on the current state vectors of each stage using the analytic hierarchy process includes: Based on the dynamic risk calculation formula, the dynamic risk value of each node at each time in each scenario is calculated, taking into account factors such as power outage loss, voltage stability, frequency stability, energy storage backup capacity, carbon emission exceedance loss, and network security risk level. The cumulative risk contribution of each node in the chain-sending process in each scenario is calculated based on the dynamic risk value of each node at each time in each scenario and the cumulative risk contribution calculation formula. The total risk of the system is calculated based on the accumulated risk contribution rate combined with the total risk calculation formula. The state vector includes: power outage loss, voltage stability, frequency stability, energy storage backup capacity, carbon emission exceedance loss, and network security risk level.

9. The method as described in claim 8, characterized in that, The dynamic risk calculation formula is shown below: DRI i,j (t)=ω1×L i,j (t)+ω2×U i,j (t)+ω3×F i,j (t)+ω4×E i,j (t)+ω5×C i,j (t)+ω6×S i,j (t) In the formula, DRI i,j (t) represents the dynamic risk value of node i in scenario j at time t, L i,j (t) represents the relative power shortage of node i in scenario j and time t, U i,j (t) represents the degree of voltage exceedance, F i,j (t) represents the degree of frequency violation, E i,j (t) represents the degree of inadequacy in energy storage backup capacity, C i,j (t) represents the degree of loss due to excessive carbon emissions, S i,j (t) represents the degree of network security risk, ω1 is the weighting coefficient of the relative power shortage of node i in scenario j and time t, ω2 is the weighting coefficient of the degree of voltage exceeding the limit, ω3 is the weighting coefficient of the degree of frequency exceeding the limit, ω4 is the weighting coefficient of the degree of insufficient energy storage backup capacity, ω5 is the weighting coefficient of the degree of carbon emission exceeding the standard loss, ω6 is the weighting coefficient of the degree of network security risk, i is the node number, j is the scenario number, and t is the time value.

10. The method as described in claim 9, characterized in that, The formula for calculating the accumulated risk contribution is as follows: CVI i,j =∫(t0to t1)DRI i,j (t)×e^(-λ×(t-t0))dt In the formula, CVI i,j Let t0 be the cumulative risk contribution of node i in the chain transmission process of scenario j, t1 be the initial time of the fault, λ be the first time the system exceeds the limit, e be the decay coefficient that is dynamically adjusted with time scale, and e be the natural constant.

11. The method as described in claim 10, characterized in that, The formula for calculating the total risk is as follows: SR j =Σ(CVI i,j ×α i ×β i ) In the formula, SR j For the total system risk, α i The sensitive load weight supplied to node i, β i This is the load level coefficient.

12. The method according to claim 1, characterized in that, The construction of the multi-drive-failure probability mapping relationship includes: Based on the meteorological data at the current moment, various extreme weather scenarios are generated using a WRF-CFD coupled model. Based on the current cybersecurity data, various cybersecurity scenarios are generated using the MITREATT&CK framework; By aligning various extreme weather scenarios and various cybersecurity scenarios along the timeline, a joint extreme weather-cybersecurity scenario set is generated. An improved K-medoids algorithm was used to reduce the number of scenarios in the extreme weather-cybersecurity joint scenario set to obtain key scenarios; Based on the physical characteristics of equipment and data-driven approaches in the aforementioned key scenarios, a functional relationship is established between extreme weather, cybersecurity attacks, and the probability of power equipment failure. The functional relationships between extreme weather, cybersecurity attacks, and power equipment failure probabilities include: the relationship between transmission line failure probabilities, the relationship between photovoltaic inverter failure probabilities, and the relationship between energy storage battery failure probabilities.

13. The method according to claim 12, characterized in that, The relationship between the fault probabilities of the transmission lines is shown in the following formula: P t Fault = P ase ×[1+β1×(v(t) / v nitical -1)+β2×(i(t)i nitical )+β3×p(t)] In the formula, P t The fault is the probability of a transmission line fault, P. ase The base failure rate, β1 is the first coefficient, v(t) is the real-time wind speed at 1km resolution output by the CFD, and v nitical The critical wind speed for line design, β2 is the second coefficient, i(t) is the icing thickness, i nitical The critical icing thickness for the line design is given by β3, which is the third coefficient, and p(t) is the probability of data tampering. The relationship between the failure probabilities of the photovoltaic inverter is shown below: P +＝ Photovoltaic fault = P ase Photovoltaics × [1+β4(t(t)t] max -1)+β5×s(t)] In the formula, P +＝ Photovoltaic fault is the probability of photovoltaic inverter failure, P. ase _PV refers to the basic failure rate of the photovoltaic inverter, β4 is the fourth coefficient, t(t) is the ambient temperature, t max β5 is the maximum allowable temperature of the inverter, s(t) is the fifth coefficient, and s(t) is the inverter infection rate. The relationship between the failure probability of the energy storage battery is shown below: P e Fault = P ase Energy storage × [1 + β6 × (1 - SOC) + β7 × r(t)] In the formula, P e Fault is the probability of failure in an energy storage battery, P ase _Energy storage refers to the basic failure rate of the energy storage battery, β6 is the sixth coefficient, SOC is the state of charge, β7 is the seventh coefficient, and r(t) is the malicious code propagation rate.

14. The method as described in claim 1, characterized in that, Also includes: Sorting the total system risk from highest to lowest, the scenarios corresponding to the first set percentage of total system risk are designated as high-risk scenarios. Based on the cumulative risk contribution, nodes are sorted from largest to smallest, and the nodes corresponding to the second-highest cumulative risk contribution percentage are designated as high-risk nodes.

15. A system for quantifying the risk of power grid chain-like events, characterized in that, include: The fault probability assessment module is used to obtain the fault probability of each device in each scenario based on the meteorological data and network security data at the current moment of each link, combined with the pre-constructed multi-drive-fault probability mapping relationship, and to determine the tripping device based on the fault probability of each device in each scenario. The parameter acquisition module is used to acquire the current meteorological data, network security data, and the power grid state vector of each link when the equipment trips. The chain-generation module is used to obtain the current state vector of each link based on the power grid state vector of the previous moment and the pre-built full-link chain-generation model. The risk quantification module is used to calculate the total system risk based on the current state vector of each link using the analytic hierarchy process (AHP). The full-chain generation model is constructed based on the state correlations in the historical data of the power grid state vector, the coupling relationships between each link, and the impact of meteorological and cybersecurity factors on the state of each link.

16. The system as described in claim 15, characterized in that, Also includes: a model building module, used for: The weight matrix of each stage is obtained by training a deep neural network based on historical data of the state vectors of each stage. Based on the correlation between the historical state vector data of each stage, the initial matrix composed of the historical state vector data of each stage is corrected to obtain the coupling matrix from each stage to the next stage. Based on historical state vector data, a joint perturbation function is constructed by combining the functional relationships between meteorological factors and state vectors in each stage and the functional relationships between network security attack factors and state vectors in each stage. A full-chain transmission model is constructed based on the weight matrix, coupling matrix, joint perturbation function, and historical state vector data of each stage.

17. The system according to claim 15, characterized in that, The full-chain transmission model is shown in the following formula: X i (t+Δt)=ReLU(W i ×X i (t)+Σ(J ij ×X i (t))+F(M(t),A(t))×Δt) In the formula, X i (t+Δt) is the state vector of element i at time t+Δt, X i (t) is the state vector of element i at time t, X i (t) is the state vector of element j at time t, W i Let J be the weight matrix of link i. ij Let F(M(t),A(t)) be the coupling matrix from stage j to stage i, and let Δt be the time step.

18. The system as described in claim 15, characterized in that, The risk quantification module is specifically used for: Based on the dynamic risk calculation formula, the dynamic risk value of each node at each time in each scenario is calculated, taking into account factors such as power outage loss, voltage stability, frequency stability, energy storage backup capacity, carbon emission exceedance loss, and network security risk level. The cumulative risk contribution of each node in the chain-sending process in each scenario is calculated based on the dynamic risk value of each node at each time in each scenario and the cumulative risk contribution calculation formula. The total risk of the system is calculated based on the accumulated risk contribution rate combined with the total risk calculation formula. The state vector includes: power outage loss, voltage stability, frequency stability, energy storage backup capacity, carbon emission exceedance loss, and network security risk level.

19. An electronic device, characterized in that, include: At least one processor and memory; The memory and processor are connected via a bus; The memory is used to store one or more programs; When the one or more programs are executed by the at least one processor, a method for quantifying urban power grid chain-related risks as described in any one of claims 1 to 14 is implemented.

20. A readable storage medium, characterized in that, It contains an execution program, which, when executed, implements a method for quantifying urban power grid chain-related risks as described in any one of claims 1 to 14.

Images