Source-grid-load-storage collaborative optimization operation method and device based on digital twinning

Abstract

The application discloses a source-network-load-storage collaborative optimization operation method and device based on digital twinning, relates to the technical field of power system optimization scheduling, and comprises the following steps: constructing a digital twin of a source-network-load-storage physical system and a correlation history database; performing scene matching based on the history database, solving an optimization model under basic operation constraints with minimization of operation cost as a target, and generating a preliminary scheduling strategy; simulating the preliminary scheduling strategy in the digital twin, calculating system transient stability margin based on an extended equal-area criterion; if the margin is lower than a preset safety threshold, identifying a predicted period with insufficient margin, calculating sensitivity of the transient stability margin to device output to determine a key device set, generating additional stability constraints or penalty terms for the predicted period and the key device set, feeding back and re-solving the optimization model for iteration until the margin meets the requirement, and determining a final scheduling strategy; and converting the final scheduling strategy into control instructions and issuing the control instructions to physical devices.

Description

Technical Field

[0001] This invention relates to the field of power system optimization and dispatching technology, and more specifically, to a method and equipment for coordinated optimization of source-grid-load-storage operation based on digital twins. Background Technology

[0002] Coordinated and optimized operation of power generation, grid, load and storage has become a key technology to ensure the safe and economical operation of the power system. Dispatching agencies need to coordinate and dispatch multiple entities such as distributed power sources, energy storage devices and controllable loads. However, current coordinated and optimized dispatching methods still have many shortcomings. Traditional methods often adopt day-ahead dispatching or rolling optimization based on prediction curves, which lack the ability to perceive and respond dynamically to the actual operating status of the system. When the output of new energy sources fluctuates drastically or the load changes suddenly, the static dispatching strategy is difficult to adjust in time, which can easily lead to system frequency deviation, voltage over-limit, or even cascading failure risks.

[0003] In existing technologies, some scheduling systems still adopt a separation mode of optimization solution and safety verification. Their safety constraint verification and optimization solution are independent of each other. When the verification fails, manual intervention is required, and there is a lack of an automated closed-loop iterative mechanism. Digital twins are used to collect measurement data to update the twin state, but the focus is mainly on equipment monitoring. The simulation and prediction capabilities of the twin are not deeply integrated with the optimization scheduling, making it difficult to obtain future dynamic prediction sequences based on the fusion results of physical simulation and data prediction.

[0004] Furthermore, existing methods lack effective means of obtaining initial values ​​when generating scheduling strategies, resulting in low optimization efficiency. At the same time, stability analysis often relies on offline simulation verification, making it difficult to calculate the system's transient stability margin in real time and generate additional stability constraints to feed back to the optimization model for iterative solution. Especially in scenarios with a high proportion of renewable energy access, there is an urgent need for a source-grid-load-storage collaborative optimization operation technology solution that integrates real-time state synchronization of digital twins, historical scenario matching, operating cost optimization, and transient stability margin verification iterative coupling. Summary of the Invention

[0005] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a source-grid-load-storage collaborative optimization operation method and device based on digital twins, in order to solve the problem that optimization solution and safety verification are independent of each other and lack an automated closed-loop iterative mechanism in the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] The method and equipment for coordinated optimization of source-grid-load-storage operation based on digital twins include the following steps: constructing a digital twin of the source-grid-load-storage physical system and its associated historical database; constructing a digital twin of the source-grid-load-storage physical system and its associated historical database; performing scenario matching based on the historical database, solving the optimization model under basic operational constraints with the goal of minimizing operating costs, and generating a preliminary scheduling strategy; simulating the preliminary scheduling strategy in the digital twin, calculating the system's transient stability margin based on the extended equal area criterion; if the margin is lower than a preset safety threshold, identifying the prediction period with insufficient margin, calculating the sensitivity of the transient stability margin to equipment output to determine the set of key equipment, generating additional stability constraints or penalty terms for the prediction period and the set of key equipment, feeding back and resolving the optimization model iteratively until the margin meets the requirements, and determining the final scheduling strategy; and converting the final scheduling strategy into control commands and issuing them to the physical equipment.

[0008] In a preferred embodiment, the construction of the digital twin of the source-grid-load-storage physical system further includes: establishing measurement equations for the twin to describe the mapping relationship between measurement values ​​and system state variables; calculating the contribution factor of each measurement point based on the Jacobian matrix of the measurement equations; setting differentiated sampling periods based on the contribution factors and collecting data from the corresponding measurement points; and performing weighted least squares state estimation on the collected data using the contribution factors as weights to update the twin state.

[0009] In a preferred embodiment, the step of performing scene matching based on a historical database to obtain initial values ​​for the solution includes: extracting feature vectors of the current system state; matching the feature vectors with a pre-built historical scene database to obtain the scene with the highest similarity and its corresponding scheduling strategy; and using the weighted average of the scheduling strategies as the initial values ​​for the optimization model, wherein the weights are proportional to the similarity of the corresponding scenes.

[0010] In a preferred embodiment, generating the preliminary scheduling strategy includes: solving for the initial values ​​and using the minimum sum of generation cost, energy storage operation cost, and load regulation cost as the objective function; wherein the decision variables include the generator output, energy storage charging and discharging power, and controllable load regulation amount for each time period; solving the optimization model and outputting the preliminary scheduling strategy.

[0011] In a preferred embodiment, the basic operating constraints include: system power balance constraints, upper and lower limits of generator active and reactive power output constraints, state of charge constraints of energy storage devices, line transmission capacity constraints, and node voltage amplitude safety constraints.

[0012] In a preferred embodiment, the optimization model further includes an energy storage health constraint, which includes: calculating the current energy storage health index based on the historical operating data of the energy storage device; reducing the upper limit of the available power of the energy storage device when the health index is lower than a preset health threshold; and adding an energy storage aging cost term to the objective function of the optimization model, wherein the aging cost is positively correlated with the charging and discharging power amplitude and the number of switching times within the scheduling cycle.

[0013] In a preferred embodiment, simulating the preliminary scheduling strategy in the digital twin further includes obtaining a future dynamic prediction sequence based on the fusion result of simulation and data prediction. Specifically, this involves: simulating the preliminary scheduling strategy in the digital twin to calculate a physical baseline prediction sequence; predicting future power fluctuations based on historical state data using a time-series neural network, wherein the loss function of the time-series neural network embeds a node power imbalance penalty term and a branch power limit violation penalty term; and superimposing the power fluctuations with the corresponding physical baseline sequence to form the future dynamic prediction sequence.

[0014] In a preferred embodiment, the calculation process for the system's transient stability margin includes:

[0015] The time-series trajectories of frequency deviation, bus voltage amplitude, and cross-sectional active power are extracted from the dynamic prediction sequence; the critical energy value is determined based on the extended equal area criterion, and the transient energy function value corresponding to the time-series trajectory is calculated; the transient energy function value is compared with the critical energy value, and the normalized difference is used as the transient stability margin.

[0016] In a preferred embodiment, the additional stability constraints or penalty terms include: adding a limit on the rate of change of output, tightening the output range, or adding a penalty term for output deviation in the objective function.

[0017] In a preferred embodiment, a computer device includes a processor and a memory; the memory stores a computer program; when the processor executes the computer program, it implements a source-grid-load-storage collaborative optimization scheduling method based on a digital twin.

[0018] The technical effects and advantages of the source-grid-load-storage collaborative optimization operation method and equipment based on digital twins of this invention are as follows:

[0019] This invention integrates measurement data acquisition from a digital twin with historical data caching. Based on historical scene matching, it performs similarity retrieval on the current system state and obtains initial values ​​for the solution. An optimization model aiming to minimize operating costs transforms these initial values ​​into power scheduling values ​​for controllable devices at different time periods. Combined with basic operational constraints, a preliminary scheduling strategy is generated, ensuring rapid convergence and practicality of the optimization solution and effectively integrating historical experience with real-time conditions. Furthermore, by inputting the preliminary scheduling strategy into the digital twin simulation and obtaining future dynamic prediction sequences based on the fusion results of physical simulation and data prediction, the system's transient stability margin is calculated, improving the safety and reliability assessment capability of the scheduling strategy. By introducing an additional stability constraint feedback iteration mechanism, additional stability constraints are automatically generated and fed back to the optimization model for iterative solution when the margin is lower than a preset safety threshold. When the margin meets the requirements, the final scheduling strategy is determined and converted into control commands issued to physical devices. This helps to achieve closed-loop coupling between optimized scheduling and transient stability verification, effectively solving the problem in existing technologies where optimization and safety verification are independent and lack an automated closed-loop iteration mechanism. Attached Figure Description

[0020] Figure 1 A schematic diagram of the source-grid-load-storage collaborative optimization operation method and equipment process based on digital twin provided in an embodiment of the present invention;

[0021] Figure 2 A comparison diagram of the actual state of the physical system and the estimated state of the digital twin provided in this embodiment of the invention;

[0022] Figure 3 This is a similarity distribution map between the current scene and the historical scene library provided in this embodiment of the invention;

[0023] Figure 4 A comparison diagram of the physical simulation baseline sequence and the fusion prediction sequence provided in the embodiments of the present invention;

[0024] Figure 5 The graph shows the change in transient stability margin with the number of iterations provided in the embodiments of the present invention. Detailed Implementation

[0025] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0026] Example 1, Figure 1 The present invention provides a source-grid-load-storage collaborative optimization operation method based on digital twins, comprising the following steps:

[0027] S1, construct a digital twin of the source-grid-load-storage physical system and its associated historical database.

[0028] S101, the digital twin of the source-grid-load-storage physical system also includes:

[0029] (1) The node-branch description method is used to establish the topology of the power distribution system, and the spatial topology of the physical power distribution system is mapped into a digital graph model. The graph model includes: the set of nodes corresponding to the bus and equipment access points in the physical system and the set of branches corresponding to the physical line connection relationships. This embodiment is based on the IEEE 33-node power distribution system and is improved. The standard IEEE 33-node system is a 12.66kV radial power distribution network with a total of 33 nodes and 32 branches. This embodiment makes the following improvements on this basis: distributed photovoltaic power sources are connected at nodes 7, 15 and 25 respectively; energy storage systems are connected at nodes 10 and 31 respectively; node 1 is the main power source node and is connected to the upper 110kV power grid through the distribution transformer. The improved system includes 1 main power source node, 3 distributed photovoltaic access nodes, 2 energy storage access nodes and 27 load nodes. The specific configuration of the distributed power source and energy storage system is shown in Table 1.

[0030] Table 1

[0031]

[0032] (2) Establish a component parameter database and map the electrical characteristic parameters of the physical equipment to the digital twin model to ensure that the electrical behavior of the digital model is consistent with that of the physical equipment. The parameters are derived from the equipment nameplate and field test data; the main parameters include:

[0033] Line parameters: The resistance, reactance, and susceptance to ground of each branch are based on the IEEE 33-bus system standard parameters; the main branch parameters are shown in Table 2:

[0034] Table 2

[0035]

[0036] Transformer parameters: main transformer capacity 10MVA, turns ratio 110kV / 12.66kV, short-circuit impedance percentage 10.5%, connection group Ynd11;

[0037] Generator and inverter parameters: The photovoltaic inverter adopts constant power factor control mode with a power factor set to 1.0; the energy storage converter adopts PQ control mode and has four-quadrant operation capability.

[0038] Load parameters: A constant power load model is adopted, and the load data of each node is configured according to the IEEE 33-node system standard. The total active load of the system is 3715kW and the total reactive load is 2300kvar.

[0039] (3) Configure a measurement system to achieve real-time synchronization of the physical system's operating status with the digital twin; the measurement data serves as the input driver for the digital twin, enabling the digital model to dynamically track the actual operating conditions of the physical system. In this embodiment, 32 measurement points are configured to achieve real-time monitoring of the voltage, power, and energy storage status of key nodes; the main measurement point configuration is shown in Table 3:

[0040] Table 3

[0041]

[0042] (4) Configure a simulation solver to realize the dynamic simulation of the digital twin.

[0043] S102, Collect measurement data to update the twin's state, specifically including the following steps:

[0044] (1) Establish measurement equations based on the network topology and parameters of the twin. The measurement equations describe the mapping relationship between the measured values ​​of each measurement point and the system state variables. The formulas are as follows:

[0045] ,

[0046] ,

[0047] ,

[0048] in, Let be a measurement vector, containing measurement data collected from all measurement points at a given time, and its dimension is . ( (Total number of measurement points) This is the system state vector (containing the voltage magnitude and phase angle of each node); Let be the measurement function vector, and its first... Each component Description of the The nonlinear mapping relationship between the measured values ​​of each measurement point and the state variables; Let the measurement error vector represent the random error in the measurement process, and assume that it follows a Gaussian distribution with a mean of zero;

[0049] Measurement function vector The specific form depends on the type of measurement, as follows:

[0050] Formula for measuring node voltage amplitude Represented as:

[0051] ;

[0052] Formula for measuring active power injected at nodes Represented as:

[0053] ,

[0054] in, and They are nodes , The voltage amplitude; The elements are the real parts of the nodal admittance matrix; These are the imaginary elements of the nodal admittance matrix; For nodes and The phase angle difference between them; This represents the total number of system nodes.

[0055] Formula for measuring reactive power injected at nodes Represented as:

[0056] ;

[0057] Branch active power measurement formula Represented as:

[0058] ,

[0059] in, , These are the branch admittance parameters;

[0060] Specifically, the network topology refers to the node connection structure of the power system in the digital twin, including: the number and type of each node (generator node, load node, tie node, etc.), the connection relationship between the first and last nodes of each branch, the connection method and tap of the transformer, and the access location of distributed power sources and energy storage devices.

[0061] Specifically, the parameters refer to the electrical characteristic parameters of power system components, including: the resistance and reactance of each branch, the susceptance to ground of each branch, the turns ratio and impedance of the transformer, the parallel susceptance of each node, and the internal electromotive force and synchronous reactance of the generator.

[0062] Specifically, the measurement points include: voltage amplitude and phase angle measurements at each node, active and reactive power measurements at each branch, output power measurements at distributed power sources, and charging and discharging power and state of charge measurements at energy storage devices; the measurement data at the measurement points are acquired in real time through a data acquisition system and transmitted to the digital twin.

[0063] It should be noted that the system state vector contains complete operating state information of the power system at a certain moment, specifically including: the voltage amplitude of each node. ( =1, 2, ..., , (total number of nodes) and the voltage phase angle of each node (With the equilibrium node as a reference, its phase angle is set to 0); the state vector can be represented as: .

[0064] (2) Calculate the Jacobian matrix of the measurement equation. The Euclidean norm of the row vector corresponding to each measurement point is used as the contribution factor for the measurement point, specifically including:

[0065] The Jacobian matrix H is a measurement function. For the state vector The first-order partial derivative matrix has a dimension of . ; The first Jacobian matrix H The row corresponds to the first Each measurement point represents the sensitivity of the measurement value to each state variable, and its expression is:

[0066] ,

[0067] in, The Jacobian matrix is ​​the first... Row vector, representing the first row vector. The sensitivity of the measured values ​​at each measurement point to each state variable; For the first Measurement function vector for each measurement point; This is the system state vector;

[0068] Measurement points The contribution factor is calculated using the Jacobian matrix. Euclidean norm of row vectors The calculation formula is as follows:

[0069] ,

[0070] in, The Jacobian matrix is ​​the first... Row vectors, where n is the dimension of the state variables. Describes the Euclidean norm of a vector. For the Jacobian matrix, the first... Line 1 The element of the column represents the first element. The measurement value for the first Sensitivity of each state variable;

[0071] In this embodiment, the contribution factor reflects the degree of influence of the measurement point on the state estimation accuracy. The larger the contribution factor, the more significant the influence of the measurement point on the state estimation result, and therefore, it needs to be sampled at a higher frequency. Specifically, measurement points with contribution factors higher than the set threshold are sampled at a base period T0, and measurement points with contribution factors lower than the set threshold are sampled at intervals that are integer multiples of the base period to obtain the sampling data of each measurement point. As shown in Table 4, the table lists the contribution factor calculation results and corresponding sampling period configurations of typical measurement points. In this embodiment, the threshold is set to 0.6, and those below 0.4 are sampled using a 3T0 period.

[0072] Table 4

[0073]

[0074] The above sampling strategy can reduce communication and computational burden while ensuring the accuracy of state estimation.

[0075] (3) Perform state estimation, the goal of which is to find a set of state variables. This minimizes the weighted sum of squared residuals for all measurement points; that is, considering the contribution capacity and measurement accuracy of each measurement point, it seeks the system operating state that best matches all measurement data.

[0076] Before performing state estimation, the collected measurement data is preprocessed, including:

[0077] (a) Bad data detection: Calculate the standardized residuals of each measurement point. If the standardized residuals exceed the preset threshold (e.g., 3.0), the data is identified as bad data and removed.

[0078] (b) Data interpolation: For data loss caused by communication interruption, linear interpolation or historical data from the same period is used to supplement it;

[0079] (c) Time alignment: Align measurement data from different sampling periods according to a unified time reference;

[0080] Objective function of state estimation for:

[0081] ,

[0082] in, The total number of measurement points. Measurement point Contribution factor Let be the objective function, representing the weighted sum of squared residuals, which needs to be minimized through optimization. The system state vector to be estimated includes the voltage magnitude and phase angle of each node; Measurement point The actual measured value; Measurement point The measurement function value, that is, the theoretical measurement value calculated based on the current state estimate; Measurement residual; represents the deviation between the actual measured value and the theoretical calculated value; Measurement point The standard deviation of the measurement reflects the measurement accuracy at that measurement point. The smaller the value, the more accurate the measurement.

[0083] The objective function for state estimation is determined using the Gauss-Newton iteration method. The optimization solution uses the following iterative formula:

[0084] ,

[0085] ,

[0086] in, This is a weight matrix weighted by contribution factors; Measurement point Contribution factors; Measurement point The standard deviation of the measurement; This is the state estimate for the k-th iteration; H represents the state estimate in the (k+1)th iteration; H is the state estimate in... The Jacobian matrix calculated at point ; W is the weight matrix weighted by contribution factors; For measurement vectors; In order to be in The measurement function vector calculated at the location; For measuring the residual vector;

[0087] In this embodiment, the convergence criterion is that the duality gap is less than 10. -6 If the number of iterations exceeds 100, the converged state estimation results are updated to the digital twin to complete the state synchronization between the twin and the physical system.

[0088] It should be noted that the purpose of the weight matrix for weighted contribution factors is to assign different weights to different measurement points during the state estimation process. This allows measurement points with large contribution factors (significantly contributing to the state estimation) and high measurement accuracy (small measurement standard deviation) to occupy a larger weight in the estimation process, thereby improving the overall accuracy and reliability of the state estimation. The weight matrix combines contribution factors with measurement accuracy, enabling differentiated processing of measurement data.

[0089] like Figure 2As shown in the figure, the comparison results between the actual state of the physical system and the estimated state of the digital twin at a certain moment are presented. The horizontal axis represents the node number, and the vertical axis represents the node voltage amplitude. The solid line represents the actual voltage of the physical system, and the dashed line represents the estimated voltage of the digital twin. It can be seen from the figure that the twin state is highly consistent with the physical system state, and the voltage estimation error of each node is within ±0.5%, which verifies the accuracy of the state estimation by this method.

[0090] S103, after completing the state estimation, link the historical database, including:

[0091] The historical database uses a sliding window mechanism, with the cache window length set to N sampling periods. When new data is written, the oldest data is automatically deleted to keep the amount of historical database data constant.

[0092] The historical database stores the following data: active power output time-series data of each distributed power source (photovoltaic, wind power, etc.), active and reactive power time-series data of each load node, charging and discharging power and state of charge time-series data of each energy storage device, and corresponding timestamp information.

[0093] S2 performs scenario matching based on historical databases, aims to minimize operating costs, solves the optimization model under basic operating constraints, and generates a preliminary scheduling strategy.

[0094] S201, Scene matching based on a historical database further includes obtaining initial values ​​for the solution. The steps for obtaining the initial values ​​include:

[0095] In this embodiment, the system has a total of Each node Individual energy storage devices, For each renewable energy source, a feature vector of the current system state is extracted. This feature vector includes the active power deficit of each node, the state of charge of energy storage, and the proportion of renewable energy output. The feature vector is then defined as follows:

[0096] ,

[0097] ,

[0098] ,

[0099] ,

[0100] in, For nodes The vacancy due to merit. For energy storage devices The state of charge, For new energy power sources The percentage of output; For nodes The load power, node The local power generation capacity, where n is the maximum total number of nodes; For energy storage devices Current energy storage, For energy storage devices The rated capacity, M is the maximum total number of energy storage devices; For new energy power sources Actual output Let k be the rated capacity of the new energy power source; in this embodiment, =33, =2, =3, therefore the feature vector dimension is 38.

[0101] The feature vector is matched with a pre-built historical scene database to retrieve the most similar historical scenes and their corresponding scheduling strategies; the construction process of the historical scene database is as follows:

[0102] (1) Extract the operating data of the past W days from the historical database of S1. In this embodiment, W is set to 90 days. The data includes: the output of distributed power sources, load power, energy storage charge status and actual scheduling strategies implemented in each time period.

[0103] (2) Divided by scheduling period (1 hour), each scenario contains the feature vector of that time period and the corresponding actual scheduling strategy. The total number of scenarios is W×24=2160.

[0104] (3) Remove abnormal operating data (such as data during a fault) and retain normal operating scenarios. After filtering, retain about 2,000 valid scenarios.

[0105] (4) The historical scene library has a capacity of 2,000 scenes and is updated using a first-in-first-out strategy. When a new scene is added and the library capacity is full, the scene with the oldest time is deleted. The scene library data structure is shown in Table 5.

[0106] Table 5

[0107]

[0108] Specifically, the similarity is calculated using cosine similarity, and the formula is:

[0109] ,

[0110] in, The first in the historical scene library Feature vectors of a scene ;

[0111] The top k most similar historical scenes are retrieved, and the weighted average of the aforementioned historical scheduling strategies is used as the initial value for solving the optimization model. The weights are proportional to the similarity of the corresponding scenes.

[0112] Specifically, the retrieval strategy used to retrieve the top k most similar historical scenes is as follows:

[0113] (1) Set the minimum similarity threshold to 0.5;

[0114] (2) Retrieve all scenarios with a similarity ≥ 0.5;

[0115] (3) Select the top k=5 scenes with the highest similarity;

[0116] (4) If the number of scenarios that meet the conditions is less than k, then the flat start strategy (the output of all equipment is set to 50% of the rated value) is adopted as the initial value.

[0117] The initial value is obtained from the above process, and its calculation formula is as follows:

[0118] ,

[0119] in, For the first The historical scheduling strategies corresponding to each matching scenario The number of historical scenes to be matched.

[0120] like Figure 3 As shown in the figure, the similarity distribution between the current scene and each scene in the historical scene library is displayed. The horizontal axis represents the historical scene number, and the vertical axis represents the similarity value. The bar chart can be used to intuitively identify the set of historical scenes that best match the current scene. In the figure, green bars represent high-matching scenes with a similarity of ≥0.8, blue bars represent medium-matching scenes with a similarity between 0.5 and 0.8, and gray bars represent low-matching scenes with a similarity of <0.5.

[0121] S202, the specific steps for generating the preliminary scheduling strategy are as follows:

[0122] Based on the initial values, the objective function is to minimize the sum of generation cost, energy storage operation cost, and load regulation cost. The decision variables include the generator output, energy storage charging and discharging power, and controllable load regulation in each time period. The optimization model is solved to output a preliminary scheduling strategy.

[0123] Specifically, the optimization model takes minimizing the sum of power generation cost, energy storage operation cost, and load regulation cost as its objective function, and uses generator output, energy storage charging and discharging power, and controllable load regulation in each time period as decision variables; the objective function expression is:

[0124] ,

[0125] in, The number of scheduling cycles. Number of generator sets For the number of energy storage devices, Number of controllable loads; For the power generation cost function, For the energy storage operating cost function, The load adjustment cost function; For generator During the period contribution; For energy storage During the period The power; For controllable load During the period The adjustment amount, This is a time index. As shown in Table 6, the main parameter configurations of the optimization model in this embodiment are listed:

[0126] Table 6

[0127]

[0128] It should be noted that the operational constraints of the optimization model include system power balance constraints, upper and lower limits of active and reactive power output of generator sets, state of charge constraints of energy storage devices, transmission capacity constraints of lines, and safety constraints of node voltage amplitude, the specific expressions of which are as follows:

[0129] (1) System power balance constraints:

[0130] ,

[0131] in, For generator During the period Those who have made meritorious contributions; For energy storage During the period The power (positive discharge); For new energy power sources During the period contribution; For nodes During the period The load power; For time period Network loss;

[0132] (2) The formulas for the upper and lower limits of active and reactive power output of generator sets are as follows:

[0133] ,

[0134] ,

[0135] in, For generator The lower limit of meritorious contribution, For generator During the period Those who have made contributions For generator The upper limit of effective output; For generator The lower limit of no-efficiency output, Let i be the reactive power output of generator i during time period t. For generator The upper limit of reactive power output;

[0136] (3) The formula for the state of charge constraint of energy storage equipment is as follows:

[0137] ,

[0138] in, For energy storage The state of charge limit; For energy storage State of charge at time t; For energy storage The upper limit of the state of charge;

[0139] (4) The recursive formula for the SOC of energy storage is as follows:

[0140] ,

[0141] in, For energy storage The state of charge at the beginning of time period t+1 The duration of a single scheduling period. For energy storage Rated energy storage capacity, Efficiency factor;

[0142] (5) The formula for calculating the line transmission capacity constraint is as follows:

[0143] ,

[0144] in, For the line During the period The transmission power; For the line Maximum transmission capacity;

[0145] (6) Node voltage amplitude The formula for calculating safety constraints is as follows:

[0146] ,

[0147] in, This is the lower limit of the node voltage amplitude; For time period The voltage amplitude; This represents the upper limit of the node voltage amplitude.

[0148] It should be noted that the basic operational constraints also include energy storage health constraints;

[0149] Specifically, the current energy storage health index is calculated based on the historical charge / discharge cycles, depth of discharge, and operating temperature of the energy storage device. The formula is as follows:

[0150] ,

[0151] in, To accumulate the number of charge and discharge cycles, The average depth of discharge. For operating temperature deviation, , , The degradation coefficient is defined as follows; historical charge / discharge records are extracted from a historical database.

[0152] When the health index falls below a preset health threshold, the maximum available power of the energy storage device is reduced in the optimization model. , The calculation formula is as follows:

[0153] ,

[0154] in, To preset health thresholds, This is the upper limit of the rated power; in this embodiment, the value is taken as 0.8.

[0155] An energy storage aging cost term is added to the objective function of the optimization model, wherein the aging cost... The aging cost is positively correlated with the charging and discharging power amplitude and the number of switching cycles within this scheduling cycle. The calculation formula is as follows:

[0156] ,

[0157] in, The power amplitude aging factor, The aging factor is the number of switching cycles. For energy storage devices During the period The charging and discharging power; For energy storage devices during the dispatch cycle The number of charge / discharge switching cycles;

[0158] The number of charge / discharge switching times The calculation method is as follows: traverse each time period within the scheduling cycle, and count a handover when the power symbol changes between adjacent time periods. The calculation formula is as follows:

[0159] ,

[0160] in, For energy storage devices During the period The charging and discharging power, For the number of periods, For time period index, It is a sign function, and its value can be -1, 0, or 1.

[0161] This step obtains initial values ​​for the solution and generates a preliminary scheduling strategy based on historical scene matching. The weighted average of historical scheduling strategies is used as the initial value for optimization, making the starting point of optimization closer to the global optimum, significantly reducing the number of iterations, and shortening the solution time compared to traditional methods, thus solving the problem of low solution efficiency in traditional optimization methods.

[0162] S3 simulates the initial scheduling strategy in a digital twin and calculates the transient stability margin of the system based on the extended equal area criterion.

[0163] S301. Based on the fusion results of simulation and data prediction, obtain the future dynamic prediction sequence. The specific steps are as follows:

[0164] (1) Input the preliminary scheduling strategy into the digital twin, calculate the node voltage and branch power through the power flow equation, calculate the system frequency response through the generator rotor motion equation, and form a physical reference prediction sequence.

[0165] The power flow equation is:

[0166] ,

[0167] ,

[0168] in, For nodes Reactive power injection; For nodes Active power injection;

[0169] The equation of motion for the generator rotor is:

[0170] ,

[0171] in, Let be the generator's inertial constant. The rotor angular velocity, For synchronous angular velocity, For mechanical power, Electromagnetic power, The damping coefficient;

[0172] The above differential equations are solved using numerical integration methods (such as the improved Euler method or the Runge-Kutta method) to obtain the frequency deviation at each time step. Harmony and angle By combining the power flow calculation results, a physical baseline prediction sequence is formed;

[0173] Frequency deviation The calculation formula is as follows:

[0174] ;

[0175] Gongjiao The calculation formula is as follows:

[0176] ;

[0177] Physical benchmark prediction sequence as follows:

[0178] ,

[0179] in, These represent the voltage amplitudes at each node at time t. These represent the voltage phase angles at each node at time t.

[0180] (2) Extract the output of distributed power sources and load power of the previous cycle from the historical database of the twin as input, and predict the power fluctuation of the future step through a time-series neural network.

[0181] Specifically, the temporal neural network adopts an LSTM (Long Short-Term Memory) structure, including an input layer, hidden layers, and an output layer, specifically including:

[0182] Input layer: dimension is The calculation formula is as follows:

[0183] ,

[0184] in, The number of periods in the historical database;

[0185] Hidden layers: 2 LSTM layers, 128 neurons per layer, activation function is hyperbolic tangent activation function;

[0186] Output layer: dimension is The calculation formula is as follows:

[0187] ,

[0188] in, To predict the number of steps; in this embodiment, =24, The expression is 3 + 33 = 36, M = 12. It is 36;

[0189] Specifically, a node power imbalance penalty term and a branch power exceeding limit penalty term are embedded into the loss function of the temporal neural network. The formula for the loss function is as follows:

[0190] ,

[0191] ,

[0192] ,

[0193] ,

[0194] in, This is the mean square error loss; , The penalty coefficients are set to 0.1 and 0.2 in this embodiment, respectively. This is a penalty term for node power imbalance. This is a penalty for exceeding the branch power limit; To predict the number of steps, Let be the predicted power fluctuation at time t. The actual power fluctuation at time t; Let be the power generation of node i at time t. Let n be the load power at node n at time t; The total number of lines, For a moment line Transmission power, For the line Maximum transmission capacity.

[0195] (3) The power fluctuation amount With physical reference sequence Corresponding superpositions form the future dynamic prediction sequence. . The calculation formula is as follows:

[0196] ,

[0197] in, As a physical reference sequence; Here, M represents the fluctuation predicted by the time-series neural network, and M is the total number of time intervals in the sequence. For example... Figure 4 As shown in the figure, the comparison between the physical simulation baseline sequence, the data prediction fluctuation, and the fused future dynamic prediction sequence is presented. The horizontal axis represents the prediction time step, and the vertical axis represents the system frequency deviation (Hz). It can be seen that the fused prediction sequence can capture the dynamic characteristics of the system more accurately.

[0198] S302, Calculate the transient stability margin. The specific steps are as follows:

[0199] (1) Extract the time-series trajectories of frequency deviation, bus voltage amplitude and cross-sectional active power from the future dynamic prediction sequence.

[0200] Specifically, the extraction of the frequency deviation time trajectory is as follows: directly read the (2n+1)th component of the predicted sequence and construct the frequency deviation time trajectory vector;

[0201] Specifically, the extraction of the bus voltage amplitude time series trajectory is as follows: directly read the first to nth components of the predicted sequence, and construct a voltage amplitude time series trajectory matrix with a dimension of n×M;

[0202] Specifically, the time-series trajectory of the active power at the cross-section is calculated as follows: based on the voltage amplitude and phase angle in the predicted sequence, the active power at the cross-section is obtained using the branch power calculation formula. The formula is as follows:

[0203] ,

[0204] in, The set of branches that make up the cross section, and Branch roads The first and last nodes; and They are time points branch road Voltage amplitude at both ends of the node; and They are time points branch road The voltage phase angle at both ends of the node.

[0205] (2) Determine the critical energy value based on the extended equal area criterion, and calculate the transient energy function value corresponding to the time series trajectory.

[0206] The determination of the critical energy value based on the extended equal area criterion includes:

[0207] Based on the power angle trajectories of each generator after the fault, generators with leading power angles are assigned to the critical generator group, and the rest are assigned to the remaining generator group; the multi-generator system is equivalent to a single-generator infinite system, with equivalent inertia coefficients. The formula is as follows:

[0208] ,

[0209] in, This is the sum of the critical group inertia. This is the sum of the inertia of the remaining aircraft group;

[0210] Equivalent angle of work Defined as the difference between the critical group's center of inertia angle and the remaining group's center of inertia angle, the formula is as follows:

[0211] ,

[0212] in, The critical angle of inertia of the machine group. The angle of inertial center of the remaining aircraft group;

[0213] Based on equivalent mechanical power With equivalent electromagnetic power The equilibrium condition determines the work angle of the unstable equilibrium point. Calculate the critical energy The formula is as follows:

[0214] ,

[0215] in, The critical energy is the maximum disturbance energy that the system can withstand under its current operating conditions. The steady-state equilibrium point power angle is the equivalent power angle of the system during steady-state operation before the fault. The power angle at the unstable equilibrium point is determined by the equivalent power balance condition, and the equivalent power is represented by a sine model. Equivalent mechanical power; This is the equivalent electromagnetic power;

[0216] The transient energy function The formula is as follows:

[0217] ,

[0218] in, For synchronous angular velocity, Let g be the current equivalent power angle of generator g.

[0219] (3) The transient energy function value With critical energy value The normalized difference between the two values ​​is used as the transient stability margin. The transient stability margin The formula is as follows:

[0220] ,

[0221] It should be noted that, A value greater than 0 indicates that the system is in a stable state. A larger value indicates a more sufficient stability margin. A value less than 0 indicates that the system is at risk of instability.

[0222] This step integrates physical simulation and data prediction to obtain future dynamic prediction sequences and calculate transient stability margins. This solves the problem that pure physical simulation is difficult to accurately predict new energy fluctuations. By integrating the prediction results with the physical simulation baseline sequence, the prediction sequence can maintain the accuracy of the physical model while capturing the random fluctuation characteristics of new energy, thus improving the prediction accuracy compared to pure physical simulation.

[0223] S4. If the margin is lower than the preset safety threshold, identify the prediction period with insufficient margin, calculate the sensitivity of transient stability margin to equipment output to determine the set of key equipment, generate additional stability constraints or penalty terms for the prediction period and the set of key equipment, feed back and re-solve the optimization model for iteration until the margin meets the requirements, and determine the final scheduling strategy.

[0224] S401: Based on the simulation results, additional stability constraints or penalty terms are generated, fed back, and the optimization model is iterated again. The specific steps are as follows:

[0225] (1) Based on the temporal distribution of the transient stability margin, identify the prediction period when the margin is lower than the preset safety threshold; let the prediction period set be... for Then the set of unstable periods for: .

[0226] (2) Calculate the sensitivity of the transient stability margin to the output of each device, and determine the devices whose sensitivity exceeds the set value as the set of critical devices; sensitivity The calculation formula is as follows:

[0227] ,

[0228] in, The current output of device i;

[0229] Key equipment collection The definition is as follows:

[0230] ,

[0231] in, The sensitivity threshold is set to 0.05 in this embodiment.

[0232] (3) For the predicted time period and the set of key equipment, add a limit on the rate of change of output or tighten the output range in the optimization model, or add a penalty term for output deviation from the steady state point in the objective function of the optimization model:

[0233] The formula for limiting the rate of change of output is as follows:

[0234] ,

[0235] in, For the output of device i in time period t, Let i be the output of device i in time period t-1. Let be the rate of change of the original maximum output of device i. The coefficient for tightening the rate of change is 0.5 in this embodiment;

[0236] The formula for tightening the output range is as follows:

[0237] ,

[0238] in, This represents the lower limit of the original output power of device i. This is to reduce the amount of force applied. For the output of device i in time period t, This represents the initial output limit of device i;

[0239] Penalty for output deviation from steady state point The formula is as follows:

[0240] ,

[0241] in, Let be the penalty coefficient for device i. This is the steady-state output point of device i.

[0242] (4) After adding the above additional stability constraints or penalty terms to the optimization model, solve it again and input the transient stability margin into the digital twin simulation again to form a closed loop iteration;

[0243] like Figure 5 As shown in the figure, the transient stability margin changes with the number of iterations. The horizontal axis represents the number of iterations, and the vertical axis represents the minimum transient stability margin. It can be seen that after 3-5 iterations, the margin converges to above the safety threshold, that is, it moves from the unsafe region to the safe region.

[0244] (5) When the number of iterations reaches the preset maximum number of iterations and the transient stability margin is still lower than the preset safety threshold, the pre-stored safety benchmark strategy is adopted as the final scheduling strategy.

[0245] The security baseline strategy is a conservative scheduling scheme that is pre-calculated and stored offline to ensure that the system can maintain stable operation under extreme conditions.

[0246] If the margin requirement is met, that is, the transient stability margin of the system is greater than or equal to the preset safety threshold in all time periods, then the current scheduling strategy is determined as the final scheduling strategy.

[0247] This step determines the final scheduling strategy by generating additional stability constraints or penalty terms and solving the problem in a closed-loop iterative manner. This solves the problem that optimal scheduling and safety and stability analysis are independent of each other. This method establishes a closed-loop iterative mechanism of "optimization-simulation-constraint feedback-re-optimization". When the transient stability margin is lower than the preset safety threshold, additional stability constraints or penalty terms are automatically generated and fed back to the optimization model, realizing the automatic coordination of economic scheduling and safety constraints. A scheduling strategy that meets both practicality and safety requirements can be obtained without manual intervention.

[0248] S5 converts the final scheduling strategy into control commands and sends them to the physical devices.

[0249] In this embodiment, the power scheduling values ​​of each controllable device in each time period in the final scheduling strategy are converted into specific control instructions, including: active power setpoint and reactive power setpoint of generator set, charging power instruction or discharging power instruction of energy device and adjustment amount instruction of controllable load;

[0250] The control commands are sent to the local controllers of each physical device through the communication network. The local controllers adjust the operating status of the devices according to the commands to complete the closed-loop control of the coordinated optimization scheduling of source, grid, load and storage.

[0251] Example 2: A computer device, including a processor and a memory;

[0252] The memory stores computer programs;

[0253] When the processor executes the computer program, it implements a method for coordinated and optimized operation of source, grid, load, and storage based on digital twins.

[0254] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0255] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0256] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0257] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0258] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0259] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A source-grid-load-storage collaborative optimization operation method based on digital twins, characterized in that, include: Construct a digital twin of the source-grid-load-storage physical system and its associated historical database; Based on historical databases, scenario matching is performed, and an optimization model is solved with the goal of minimizing operating costs to generate a preliminary scheduling strategy. A preliminary scheduling strategy is simulated in a digital twin, and the power fluctuation is predicted by a time-series neural network. The loss function of the time-series neural network is embedded with node imbalance penalty terms and branch over-limit penalty terms. The predicted fluctuation is superimposed with the simulated physical baseline sequence to obtain a dynamic prediction sequence. The time-series trajectory is extracted from the dynamic prediction sequence, the transient energy function value is calculated based on the extended equal area criterion, and the transient stability margin is determined accordingly. If the transient stability margin is lower than a preset threshold, identify the period of insufficient margin, calculate the sensitivity of the transient stability margin to the equipment output to determine the set of key equipment, generate additional stability constraints or penalty terms for the key equipment and feed them back to the optimization model to solve again, iterate until the transient stability margin meets the requirements, determine the final scheduling strategy, and convert the final scheduling strategy into control commands and send them to the physical equipment.

2. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 1, characterized in that, The digital twin of the source-grid-load-storage physical system also includes: Establish measurement equations for the twin to describe the mapping relationship between measured values ​​and system state variables; Calculate the contribution factor of each measurement point based on the Jacobian matrix of the measurement equation. Differentiated sampling periods are set based on contribution factors, and data from corresponding measurement points are collected. Using contribution factors as weights, weighted least squares state estimation is performed on the collected data to update the twin state.

3. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 2, characterized in that, The scenario matching based on the historical database includes: Extract the feature vector of the current system state; The feature vector is matched with a pre-built historical scene library to obtain the scene with the highest similarity and its corresponding scheduling strategy; The weighted average of the scheduling strategy is used as the initial value for solving the optimization model, where the weights are proportional to the similarity of the corresponding scenarios.

4. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 3, characterized in that, The generation of the preliminary scheduling strategy includes: Based on the initial values, the objective function is to minimize the sum of power generation cost, energy storage operation cost, and load regulation cost. The decision variables include the power output of generator units, the charging and discharging power of energy storage, and the controllable load regulation amount in each time period. Solve the optimization model and output a preliminary scheduling strategy.

5. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 4, characterized in that, The basic operational constraints of the optimization model include: system power balance constraints, upper and lower limits of active and reactive power output of generator sets, state of charge constraints of energy storage devices, line transmission capacity constraints, and node voltage amplitude safety constraints.

6. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 1, characterized in that, The optimization model also includes energy storage health constraints, which include: Calculate the current energy storage health index based on historical operating data of energy storage devices; When the health index is lower than the preset health threshold, the upper limit of the available power of the energy storage device is reduced. An energy storage aging cost term is added to the objective function of the optimization model. The aging cost is positively correlated with the charging and discharging power amplitude and the number of switching times within the scheduling cycle.

7. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 1, characterized in that, The temporal neural network is trained based on historical state data.

8. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 1, characterized in that, The timing trajectory includes frequency deviation, bus voltage amplitude, and cross-sectional active power.

9. The source-grid-load-storage coordinated optimization operation method based on digital twins according to claim 1, characterized in that, The additional stability constraints or penalty terms include: adding a limit on the rate of change of output, tightening the output range, or adding a penalty term for output deviation in the objective function.

10. A computer device, characterized in that, Including processor and memory; The memory stores computer programs; When the processor executes the computer program, it implements the method as described in any one of claims 1 to 9.

Images