A power spot market simulation clearing method and system

Abstract

The application discloses a power spot market simulation clearing method, which combines electrical distance and dynamic boundary output to perform spatial division; adopts a dynamic programming algorithm integrated with a multi-target cost function and unit constraint to perform time division; establishes a mixed integer programming model of a time-space block and a master-slave iteration framework with the mixed integer programming model as a sub-problem; the specific steps are as follows: S1, data acquisition and pretreatment; S2, spatial division; S3, time division; S4, generation of a time-space block and a mixed integer programming model of the time-space block; S5, establishment of a master-slave iteration framework; and S6, iteration solving process. The application further discloses a power spot market simulation clearing system.

Description

Technical Field

[0001] This invention relates to the fields of electricity market operation and power system planning, and more specifically, to a simulated clearing method and system for the electricity spot market. Background Technology

[0002] In the fields of power market operation and power system planning, 8760-hour clearing simulation is one of the core technical means. Its purpose is to provide key decision-making basis for power market mechanism design, power source planning optimization, and grid operation and dispatch by simulating the hourly power supply and demand balance, unit operating status, and electricity price formation process throughout the year. With the continuous increase in the proportion of renewable energy installed capacity in the power system and the continuous access of flexible resources on the user side, the volatility and uncertainty of power system operation have increased significantly, which puts forward higher requirements for the computational efficiency and solution accuracy of 8760-hour clearing simulation. Currently, the mainstream 8760-hour clearing simulation methods in the industry generally adopt a single linear programming model and perform serial solutions hourly through time decoupling. The core idea of ​​this technical solution is to decompose the continuous clearing problem of 8760 hours throughout the year into 8760 independent hourly clearing subproblems. Each subproblem is based solely on the load data, renewable energy output forecast data, and unit parameters of the current hour. A linear programming algorithm is used to solve for the unit combination and power trading clearing results within that hour. The solution results of the previous hour are then used as the initial conditions for the next hour, and the simulation calculations for all hours are completed sequentially.

[0003] The drawbacks of existing technologies are: 1. Serial solving leads to a linear increase in computation time and hours. Since the 8760-hour clearing calculation needs to be performed hourly in chronological order, and each hour's linear programming solution consumes computational resources and time, the total computation time of the entire simulation process exhibits a strict linear positive correlation with the number of hours. Under the requirements of refined operation in the power market, if iterative verification of unit parameters and market rules is required across multiple scenarios, the cumulative computation time will increase significantly, severely reducing the efficiency of simulation analysis.

[0004] 2. The cross-time coupling characteristics of unit ramping constraints lead to reduced solution accuracy. In actual power system operation, unit ramping constraints (i.e., the limitation on the range of output change of units per unit time) are key technical constraints to ensure system frequency stability. This constraint has significant cross-time coupling characteristics. For example, the output adjustment range of a unit in a certain hour depends not only on the current hour's operating status but also on the direct limitation of the unit's actual output in the previous hour. However, the existing time-decoupled hourly solution method can only simplify the ramping constraint of the current hour based on the solution results of the previous hour. It cannot fully consider the continuity and constraint correlation of unit output changes between adjacent hours, resulting in deviations between the simulated unit output curves and actual operating requirements. This, in turn, affects the accuracy of electricity price calculation results and makes it difficult to truly reflect the operating rules of the electricity market.

[0005] In the prior art, the publication number CN121543934A, titled "A Simulated Clearing Method for the Energy Storage Electricity Spot Market Based on Hierarchical Model Control Prediction," discloses a simulated clearing method for the energy storage electricity spot market based on hierarchical model control prediction. This method includes the following steps: S1, establishing a globally simplified clearing model for the electricity spot market including energy storage as the upper-level predictive control model; S2, establishing a detailed clearing model for the electricity spot market including energy storage as the lower-level execution model; S3, dividing the long-term simulation into M solution cycles; in each solution cycle, first using the globally simplified clearing optimization to obtain the remaining energy storage final-state SOC constraint values ​​for each cycle; then, based on the energy storage final-state SOC constraint values ​​for this cycle, using the detailed clearing model to complete the market simulation clearing for this cycle; finally, updating the system state according to the current clearing results and rolling forward to the next cycle until all cycles are cleared. This prior art achieves dynamic global optimization of energy storage SOC constraints through a hierarchical model predictive control mechanism, balancing computational efficiency and market simulation clearing effectiveness. Summary of the Invention

[0006] In view of the above problems, the present invention provides a method and system for simulated clearing of the electricity spot market, the main purpose of which is to solve the computational efficiency of spatiotemporal dual decomposition and coordinate the boundaries of each spatiotemporal block.

[0007] To achieve the above-mentioned objectives, the technical solution adopted by this invention is: a simulated clearing method for the electricity spot market, the specific steps of which are as follows: Step S1: Data acquisition and preprocessing;

[0008] S11. Acquisition and preprocessing of annual load forecast data; The power system dispatch center periodically forecasts the subsequent electricity demand of the power grid within its dispatch area and stores it in the database. The load forecast data includes timestamps and corresponding load values, and is then processed. S12. Acquisition and standardization of unit parameters; Obtain parameters for various thermal power units and hydropower units, convert all unit parameters to per-unit values, and use the system baseline capacity as the benchmark: = / ,in Per unit value, This represents the actual active power of the generator unit; S13. Extraction and structuring of network constraint data; The network topology and operational constraint data of the power grid are obtained from the power grid energy management system (EMS), including node data, branch data, and transformer data. This data is then used to construct a node-branch correlation matrix and an admittance matrix. S14. Data consistency verification and time alignment; Finally, perform consistency checks across data sources: ensure all data uses the same timestamp system; verify the match between the total unit capacity and the system's maximum load; and check the correspondence between network constraints and unit locations. After completing the above steps, the three types of data are encapsulated into structured data objects in a unified format, providing standardized input for subsequent spatiotemporal decomposition steps.

[0009] Step S2: Spatial partitioning;

[0010] S21. Electrical distance calculation and similarity matrix construction; First, calculate the electrical distance between all nodes in the system as the basis for network partitioning; S22. Power Grid Diagram Model Construction and Community Discovery; The power grid is abstracted as a graph model G=(V,E,W), where: V is the set of nodes, containing all the main nodes; E is the set of edges, where each edge represents the connection between nodes; W is the edge weight matrix. ; Network partitioning is performed based on the power grid diagram model; S23. Zoning quality assessment and adjustment; S24. Boundary node identification and tie-line processing; After the space is divided, the boundary nodes and connecting lines between each area are identified. Boundary nodes are nodes that are electrically connected to multiple areas at the same time, and connecting lines are power transmission lines that connect different areas. S25. Final partitioning scheme output; The final output zoning scheme includes: (1) a regional division list: the nodes, units, and load information contained in each region; (2) a boundary node mapping table: recording the multiple regions to which each boundary node belongs; (3) tie line parameters: including impedance limit and thermal stability limit; (4) a regional characteristic summary: the power generation capacity, load level, and network density of each region.

[0011] Step S3: Time division;

[0012] S31. Load feature extraction; S32. Dynamic Programming: State Definition and State Transition; The goal of dynamic programming is to divide the time into several time blocks, so that the load fluctuation of each time block is small, the cost estimate is stable, and a balance is achieved between reducing computational complexity and minimizing accuracy loss. S33. Dynamic programming solution and post-optimization adjustment; The reverse dynamic programming algorithm is used to solve the problem recursively. After the solution is completed, a preliminary time block partitioning scheme is obtained by backtracking.

[0013] Step S4: Generate a spacetime block and a mixed-integer programming model for the spacetime block;

[0014] S41. Spatiotemporal block generation and data allocation; Based on the M region division results obtained in step S2 and the K time block division results obtained in step S3, M×K spatiotemporal blocks are generated, and each spatiotemporal block corresponds to a sub-problem of a region in a specific time period. Assign corresponding input data to each spatiotemporal block, including load data of all nodes in the region within the time block, extract technical parameters of all units in the region, network constraints and boundary connection constraints within the region; S42. Construction of a mixed-integer programming model; A mixed-integer programming model is constructed for each spatiotemporal block; the constraints include: power balance constraints, unit operation constraints, network security constraints, and time boundary consistency constraints.

[0015] Step S5. Establish a master-slave iteration framework;

[0016] S51. Master-slave architecture design; The master-slave iterative framework is a two-layer optimization framework, including a coordination layer for the master problem and a computation layer for the subproblems. The coordination layer is responsible for coordinating the boundary conditions between subproblems in various regions, maintaining global consistency constraints, ensuring the balance of power exchange between regions, and generating and updating coordination signals. The computation layer is used to solve the mixed integer programming model of each spatiotemporal block in parallel, submit local optimization results to the master problem, and receive and respond to the coordination signals issued by the master problem. S52. Modeling the main problem; The main problem is to use price signals as a coordinating tool to generate coordinating signals; S53. Sub-problem coordination mechanism; Each subproblem receives a coordination signal from the main problem during the solution process and incorporates it into its local optimization objective. After the subproblem is solved, it submits data to the main problem.

[0017] Step S6. Iterative solution process;

[0018] S61. Iterative solution process; The alternating direction multiplier method framework is used for master-slave collaborative solution. The specific process is as follows: Set the initial iteration counter k=0 and initialize the price signal. and At the same time, set the maximum number of iterations and the convergence tolerance; Then, each computing node solves its assigned spatiotemporal block mixed integer programming model in parallel, taking into account the coordination signal provided by the main problem during the solution process; after the management node collects the solution results of all subproblems, the main problem updates the price signal and calculates the degree of global constraint violation. If the degree of constraint violation is minimal and the convergence condition is met, the solution is considered complete; otherwise, an updated coordination signal is issued for the next round of calculation. After the solution is completed, a weighted average method is used to ensure consistency for boundary variables; S62. Clearing price calculation; Calculate the nodal marginal price (LMP) based on the final optimization results:

[0019] in, The system energy price component is equal to the dual variable of the global equilibrium constraint; The blocking price component of node i is equal to the sum of the dual variables of the corresponding network constraints; The network loss price component for node i is calculated using the network loss distribution factor. S63. Output report compilation; The final output report includes unit scheduling plan, clearing price curve, and auxiliary information.

[0020] This invention also discloses a simulated clearing system for the electricity spot market, characterized in that: The data acquisition and preprocessing module is used to perform the above step S1; The spatial partitioning module is used to perform step S2 described above; The time division module is used to execute step S3 above; The spatiotemporal block generation and mixed-integer programming module is used to perform step S4 above; The master-slave iteration framework establishment module is used to execute step S5 above; The iterative solution module is used to execute step S6 as described above.

[0021] The present invention also discloses a computer program product, including a computer program, characterized in that the computer program, when executed by a processor, implements the simulated clearing method for the electricity spot market of the invention. Beneficial effects

[0022] The present invention provides a simulated clearing method and system for the electricity spot market, which improves computational efficiency through spatiotemporal dual decomposition; and improves the cross-time period coupling constraint satisfaction rate by coordinating the boundaries of each spatiotemporal block through the master problem in the master-slave iterative framework. Attached Figure Description

[0023] Figure 1 This is a schematic diagram of the electricity spot market simulated clearing method according to an embodiment of the present invention. Detailed Implementation

[0024] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0025] This embodiment of the simulated clearing method for the electricity spot market combines electrical distance and dynamic boundary output for spatial partitioning; it employs a dynamic programming algorithm integrating multi-objective cost functions and unit constraints for time partitioning; it establishes a mixed-integer programming model for spatiotemporal blocks, and a master-slave iterative framework using this mixed-integer programming model as subproblems; and it achieves simulated clearing of the electricity market for 8760 hours (year-round); Figure 1 The flowchart shown below illustrates the specific steps as follows: S1. Data Acquisition and Preprocessing.

[0026] S11. Acquisition and preprocessing of annual load forecast data.

[0027] The power system dispatch center periodically forecasts the subsequent electricity demand of the power grid within its dispatch area and stores the forecasts in a database. In other words, the data processed by this invention pertains to the entire power grid.

[0028] Retrieve the latest 8760-hour load forecast data from the dispatch center's database, including timestamps and corresponding load values ​​(in MW). Then perform the following processing: First, a data quality check is performed to identify and handle outliers. A statistical method is used to calculate the Z-score (Z=(x-μ) / σ) of each hourly load forecast (hereinafter referred to as load value). When |Z|>3, it is identified as an outlier and replaced with the average load of the two hours before and after it.

[0029] Then, data normalization is performed, mapping the load values ​​to the [0,1] interval. The Min-Max normalization method is used: Where L is the original load value, and These are the minimum and maximum load values ​​for the entire year, respectively.

[0030] S12. Acquisition and standardization of unit parameters.

[0031] Detailed parameters of various generator sets are obtained from data submitted by power generation companies and from unit characteristic databases, including but not limited to: For thermal power units: collect rated capacity ( ), minimum technical output ( ), start-up and shutdown costs ( The parameters include operating cost coefficients (a, b, c, where the cost function is C = a + bP + cP²), and ramp rate. The cost coefficients need to be updated based on the latest fuel prices and recalibrated using a linear regression method.

[0032] For hydropower units: In addition to the parameters mentioned above, it is also necessary to collect constraint parameters such as reservoir regulation characteristics and minimum discharge flow.

[0033] Then, all unit parameters are converted to per-unit values ​​based on the system baseline capacity ( (Usually 100 MVA) is taken as the benchmark: = / ,in This is a per-unit value, without a unit, or the unit can be specified as pu; This represents the actual active power of the generator unit. This ensures that all unit parameters use the same baseline value, avoiding unit mismatch issues in subsequent parallel calculations.

[0034] S13. Extraction and structuring of network constraint data.

[0035] Obtain network topology and operational constraint data of the power grid from the power grid energy management system (EMS), including: Node data: Voltage limits for each node ( , ), reference voltage level; Branch data: Transmission line resistance (R), reactance (X), susceptance (B), thermal stability limit; Transformer data: turns ratio, impedance, tap range (range for fine-tuning the turns ratio).

[0036] Then, these data are used to construct a node-branch association matrix A (with dimensions N×B, where N is the number of nodes and B is the number of branches), and the matrix elements... This represents the connection relationship between node i and branch k (0 indicates no connection, 1 indicates connection in positive direction, and -1 indicates connection in negative direction). Simultaneously, the admittance matrix Y = G + jB is established, where G is the conductance matrix and B is the susceptance matrix.

[0037] S14. Data consistency verification and time alignment.

[0038] Finally, perform consistency checks across data sources: ensure all data uses the same timestamp system; verify the match between the total unit capacity and the maximum system load (requiring a reserve rate of no less than 15%); and check the correspondence between network constraints and unit locations.

[0039] After completing the above steps, the three types of data are encapsulated into structured data objects in a unified format, providing standardized input for subsequent spatiotemporal decomposition steps.

[0040] S2, Spatial division.

[0041] S21. Electrical distance calculation and similarity matrix construction.

[0042] First, the electrical distance between all nodes in the system is calculated as the basis for network partitioning. An electrical distance calculation method based on the impedance matrix is ​​used: Define the electrical distance matrix D, whose elements The electrical distance between node i and node j is represented by the following formula:

[0043] in: Let be the voltage phase angle (in radians) at node i. Let be the voltage amplitude (per unit) at node i. This is the system's rated voltage (per unit, usually 1.0). This is the transfer impedance (per unit) between nodes i and j. This is a weighting factor (taken as 0.2-0.5) used to balance the effects of voltage differences and impedance.

[0044] The α value is automatically adjusted according to the system's operating status, decreasing under heavy load and increasing under light load.

[0045] S22. Power Grid Graph Model Construction and Community Discovery.

[0046] The power grid is abstracted as a graph model G=(V,E,W), where: V is the set of nodes, containing all the main nodes; E is the set of edges, where each edge represents the connection between nodes; W is the edge weight matrix. .

[0047] Then, spectral clustering, the Kernighan-Lin algorithm, or the Louvain community detection algorithm can be used to partition the network based on the power grid graph model. When the power grid is small and the desired partitioning area is small, spectral clustering can be used for balanced and theoretically optimal partitioning; when the power grid is large and the desired partitioning area is large, the Louvain community detection algorithm can improve computational efficiency.

[0048] Here, we will take the Louvain community detection algorithm as an example to illustrate some points. The implementation process of the Louvain community detection algorithm is as follows.

[0049] First, assign an independent community label to each node (initial number of communities = total number of nodes), then calculate the total weight m of all edges, and the weighted degree of each node i. .

[0050] Iterative execution: Phase 1: Node Attribution Optimization Iterate through all nodes (in numerical order), for each node i: a. The modularity contribution of computing node i to its current community; b. Traverse all adjacent nodes of node i ( For nodes j with a value greater than 0, calculate the change in module degree contribution ΔQ when migrating i to the community where j resides; c. Filter migration schemes that have ΔQ > 0 and meet the constraints, and select the scheme with the largest ΔQ to perform the migration; d. If all migration schemes ΔQ≤0 or do not satisfy the constraints, maintain the original affiliation of node i; Repeatedly traverse the nodes until no node migrates after one round of traversal, or the ΔQ caused by migration is less than 10^{-4}, which is considered as local convergence.

[0051] Phase Two: Community Aggregation Treat each community as a "super node," calculate the edge weights between super nodes, which is the sum of the edge weights of all node pairs within two communities; then recalculate the weighted degree of the super nodes. 'Total weight m', proceed to the next iteration (repeated phase 3.2.1-3.2.2).

[0052] When the improvement in modularity Q after iteration is ΔQtotal=Qnew−Qold<10^{-5}, the iteration stops and the final community partitioning result is output.

[0053] In step a, the formula for calculating node i's modularity contribution to its current community C is:

[0054] Where m is the total weight. It is the total weighted degree of community C.

[0055] In step b, the change in module degree contribution of node i before and after migration is ΔQ = , The calculation formula is the same as above, except that community C adaptability is changed to D.

[0056] In the second stage, the modularity Q corresponds to the entire power grid diagram model and is used to measure the quality of community partitioning. Its value ranges from [-1, 1], with a better partitioning effect being closer to 1. The calculation formula is as follows:

[0057] Among them, indicator function This means that if node i and j are in the same community, the value is 1; otherwise, it is 0. These represent the communities where nodes i and j reside, respectively.

[0058] The verification process includes: 1. The final modularity must be greater than 0.3 (engineering statistics). 2. The following constraints are met: the size of all communities is greater than the preset minimum value; the key node combinations of the power grid are not split and belong to the same community; 3. The voltage deviation of nodes within each community is smaller than the voltage deviation between communities.

[0059] S23. Zoning quality assessment and adjustment.

[0060] Calculate the following metrics to evaluate the quality of the partition: (1) Coupling degree within the region: , The value must be greater than 0.85.

[0061] (2) Inter-regional coupling: Where A and B are indices for different regions. The value must be less than 0.15.

[0062] (3) Power balance: , The requirement is less than 0.1. Through iterative adjustments, ensure that the above partitioning quality requirements are met.

[0063] S24. Boundary node identification and tie line processing.

[0064] After the space is divided, the boundary nodes and connecting lines between each area are identified. Boundary nodes are nodes that are electrically connected to multiple areas at the same time, and connecting lines are power transmission lines that connect different areas.

[0065] For boundary nodes and tie lines, a dual attribution approach is adopted. This means that in subsequent calculations, these nodes simultaneously participate in the optimization calculations of adjacent regions and coordinate their power exchange through the master node. Additionally, virtual overlap zones are set at region boundaries. This involves creating additional virtual direct adjacent nodes of the boundary nodes in the optimization calculation model for each region. These virtual nodes participate in the model calculations, ensuring that the decomposed subproblems accurately reflect the coupling effects between regions. However, since they are essentially adjacent regions, they do not directly execute the scheduling results obtained from the model calculations corresponding to their "local region." Instead, they provide coupling references through a region coordination mechanism.

[0066] S25. Final partitioning scheme output.

[0067] The final output zoning scheme includes: (1) a regional division list: the nodes, units, and load information contained in each region; (2) a boundary node mapping table: recording the multiple regions to which each boundary node belongs; (3) tie line parameters: including impedance limits, thermal stability limits, etc.; (4) a regional characteristic summary: the power generation capacity, load level, network density, etc. of each region.

[0068] By introducing electrical distance weighting, dynamic boundary treatment, and a quality assessment system, the decomposed sub-problems are ensured to maintain the integrity of electrical characteristics while minimizing the degree of coupling between regions.

[0069] S3, Time Division.

[0070] S31. Load feature extraction.

[0071] For the 8760-hour load forecast curve obtained in step S1, calculate its hourly load change rate and daily load fluctuation coefficient.

[0072] Hourly load change rate : ,in, Let t be the load value for hour t.

[0073] Daily load fluctuation coefficient : ,in, These represent the maximum load, minimum load, and average load on day d, respectively.

[0074] S32. Dynamic Programming: State Definition and State Transition.

[0075] The goal of dynamic programming is to divide the 8760 hours into several time blocks, so that the load fluctuation of each time block is small, the cost estimate is stable, and a balance is achieved between reducing computational complexity and minimizing accuracy loss.

[0076] First, define the state variables for dynamic programming. Let represent the optimal partitioning scheme when the time block reaches the i-th hour. Each state contains the following information: the start time of the current time block. The end time of the current time block Maximum load fluctuation within the current time block Cost estimate for the current time block .

[0077] For example, This represents the optimal partitioning scheme from hour 1 to hour 10, with the current time block being the time block ending at hour 10.

[0078] Here is a brief explanation of the dynamic programming calculation logic used in this scheme: The problem of dividing 8760 hours is essentially about choosing some of the 8759 time intervals (1-2 hours, 2-3 hours, ..., 8759-8760 hours) as the dividing points. The total number of possible dividing schemes is 2^8759.

[0079] Dynamic programming can be viewed as sequentially applying... arrive The calculation, the calculation This requires enumerating all possible pre-split points i (where i ranges from 1 to j-1), calculating the total cost for each i, and then selecting the i with the minimum total cost as the optimal pre-split point. A specific example is as follows: right Only one calculation is required (no i can be selected); right : We need to enumerate i=1 (1 calculation); right : We need to enumerate i=1, 2 (2 calculations); ... right : We need to enumerate i=1 to 8759 (8759 calculations).

[0080] Thus, the total computational cost is 1 + 2 + 3 + ... + 8759, which is in the millions, reducing the computational cost by countless orders of magnitude compared to brute-force enumeration.

[0081] Returning to the solution, the state transition equation is:

[0082] in, Let represent the optimal partitioning cost from hour 1 to hour i, which is the sum of the estimated costs of the time blocks partitioned from hour 1 to hour i; Δ(i,j) represents the maximum load fluctuation rate within the time block from hour i to hour j. The length of the time block is represented by α, which is the load fluctuation weighting coefficient (0.6-0.8), and β is the time block length weighting coefficient (0.2-0.4).

[0083] Cost function for a single time block:

[0084] in, This represents the maximum load fluctuation rate within that time block. This represents the deviation between the time block length and the ideal length (72 hours). This is an estimate of the number of units that need to be started or stopped within this time block. The weighting coefficients must satisfy the condition that the sum of the three is 1.

[0085] The constraints include: (1) Minimum time block length constraint: Each time block must contain at least 24 hours (to ensure minimum unit operating time); (2) Maximum time block length constraint: Each time block contains a maximum of 168 hours (to avoid distortion of load characteristics); (3) Load fluctuation constraint: The maximum load fluctuation within a single time block shall not exceed 30%; (4) Unit operation constraints: The time block division needs to take into account the minimum start-up and shutdown time of the unit.

[0086] S33. Dynamic programming solution and post-optimization adjustment.

[0087] A reverse dynamic programming algorithm is used to recursively solve the problem from t=8760 hours to t=1 hour, storing the optimal partitioning scheme and corresponding cost value for each hour. The specific implementation is similar to the above idea, except that the order is reversed.

[0088] After solving the problem, a preliminary time block partitioning scheme is obtained through backtracking. The details are as follows: Starting from t=1, select the hour value k1 of the optimal partitioning scheme as the end hour of the first time block (denoted as k1), so the first time block is 1~k1 hours; jump to k1+1 hours, check the optimal hour value k2 corresponding to [k1+1], and get the second time block k1+1~k2 hours; repeat the above steps until all consecutive time blocks are obtained (such as 1~k1, k1+1~k2, ..., kn+1~8760).

[0089] Then, time blocks shorter than 24 hours are merged with adjacent time blocks; time blocks longer than 168 hours with large load fluctuations are split; and the boundary times are adjusted to ensure that the time block boundaries fall at times when load changes are relatively gentle.

[0090] By introducing a multi-objective cost function, unit constraint integration, and post-optimization adjustment mechanism, the partitioning results are ensured to conform to the load change pattern and meet the requirements of subsequent optimization calculations.

[0091] S4. Generate spacetime blocks and a mixed integer programming model for spacetime blocks.

[0092] S41. Spatiotemporal block generation and data allocation.

[0093] Based on the M region division results obtained in step S2 and the K time block division results obtained in step S3, M×K spatiotemporal blocks are generated, and each spatiotemporal block corresponds to a sub-problem of a region in a specific time period.

[0094] Assign corresponding input data to each spatiotemporal block, including load data of all nodes in the region within the time block, extract the technical parameters of all units in the region, and network constraints and boundary connection constraints within the region.

[0095] S42. Construction of mixed integer programming model.

[0096] For each spatiotemporal block, construct a mixed-integer programming model of the following form: The objective function is:

[0097] in, Let be the power generation cost function of unit j; , These represent the start-up and shutdown costs of unit j, respectively. , This represents the start / stop status of unit j, with a value of 1 or 0. The power generation of unit j at time t. The time range of the spatiotemporal block is given by t, where t is the hour index. It is a set of units within a spacetime block.

[0098] The constraints include: (1) Power balance constraint:

[0099] in, Let be the power transmitted by the tie line l at time t (positive indicates input, negative indicates output). Let be the set of connect lines within the spatiotemporal block region m. Let i be the index of the load node. Let i be the active power of load node i at hour t. Let be the estimated network loss of region m at time t.

[0100] (2) Unit operating constraints:

[0101]

[0102] Where I_j(t) is a 0-1 variable, representing the operating state of unit j at time t; Let be the rated gradeability of unit j.

[0103] (3) Network security constraints:

[0104]

[0105] in, Let be the voltage phase angle of node i at time t.

[0106] (4) Time boundary consistency constraints: ,in, , These represent the operating status and active power of unit j at the beginning hour of the time-space block.

[0107] S5. Establish a master-slave iteration framework.

[0108] S51. Master-slave architecture design.

[0109] The master-slave iterative framework of this scheme is a two-layer optimization framework, including a coordination layer for the main problem and a computation layer for the subproblems. The master-slave iterative framework is implemented on a distributed computing platform. The main problem is deployed on the management node and solved using CPLEX or Gurobi solvers; the subproblems are distributed across multiple computing nodes and solved in parallel.

[0110] The coordination layer is responsible for coordinating the boundary conditions between subproblems in each region, maintaining global consistency constraints, ensuring the balance of power exchange between regions, and generating and updating coordination signals (such as Lagrange multipliers or penalty terms).

[0111] The computational layer is used to solve the mixed-integer programming model of each spatiotemporal block in parallel, submit local optimization results to the main problem, and receive and respond to coordination signals issued by the main problem.

[0112] S52. Modeling the main problem.

[0113] The main problem is constructed as a linear programming problem of the following form: Its decision variables include the tie-line power price signal between regions m and n at time t. And the reserve capacity price signal between regions m and n at time t. .

[0114] The objective function is

[0115] Where t is the hour index, T equals 8760, m is the region index, M is the total number of regions, and N(m) represents the set of neighboring regions of region m. This indicates the power deviation of the tie line. This represents the deviation in reserve capacity. The deviation here is the difference between the planned value reported by the sub-problem and the ideal optimal value.

[0116] The constraints include power balance constraints: and system backup requirements constraints ,in, This represents the system's reserve capacity requirement at time t.

[0117] The main problem is to use price signals as a coordinating tool to generate coordinating signals:

[0118]

[0119] Where: ρ and σ are step size parameters that control the convergence speed, and k is the number of iterations.

[0120] S53. Sub-problem coordination mechanism.

[0121] Each subproblem receives a coordination signal from the main problem during the solution process and incorporates it into its local optimization objective:

[0122] The original objective is the objective function of the mixed integer programming model in step S4; , It is a coordination signal issued from the main issue.

[0123] After solving the subproblems, the data submitted to the main problem includes: planned power values ​​for the boundary tie lines. Boundary reserve capacity allocation value Locally optimize the objective function value and constrain violation information.

[0124] The convergence conditions are as follows:

[0125]

[0126]

[0127] in, , , To reduce convergence tolerance, Let $\frac{ ...

[0128] S6. Iterative solution process.

[0129] S61. Iterative solution process.

[0130] This invention employs the Alternating Direction Multiplier Method (ADMM) framework for master-slave collaborative solution, the specific process of which is as follows: Set the initial iteration counter k=0 and initialize the price signal. and Specifically, it can be set to a zero vector or a predicted value based on historical data. The maximum number of iterations and convergence tolerance can also be set.

[0131] Then, each computing node solves its assigned spatiotemporal block mixed-integer programming model in parallel, taking into account the coordination signal provided by the main problem during the solution process. After the management node collects the solution results of all subproblems, the main problem updates the price signal and calculates the degree of global constraint violation.

[0132] If the degree of constraint violation is minimal and the convergence condition is met, the solution is considered complete; otherwise, an updated coordination signal is issued for the next round of calculation.

[0133] After the solution is completed, a weighted average method is used for boundary variables to ensure consistency.

[0134] S62. Clearing price calculation.

[0135] Calculate the nodal marginal price (LMP) based on the final optimization results:

[0136] in, The system energy price component is equal to the dual variable of the global equilibrium constraint; The blocking price component of node i is equal to the sum of the dual variables of the corresponding network constraints; Let be the network loss price component of node i, calculated using the network loss distribution factor.

[0137] S63. Output report compilation.

[0138] The final output report includes: 1. Unit scheduling plan, which includes the following for each unit at each time: start / stop status (0 / 1 variable), power generation value (continuous variable), reserve capacity allocation, and operating cost information.

[0139] 2. The clearing price curve includes the following components for each node at each time point: energy price component, congestion price component, network loss price component, and final node marginal price.

[0140] 3. Auxiliary information: computational performance statistics (solution time, number of iterations, convergence status), constraint violation reports, system standby status, and network congestion status.

[0141] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit it. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features, and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A simulated clearing method for the electricity spot market, characterized in that, The specific steps are as follows: S1. Data acquisition and preprocessing; S2, Spatial Division; S3, Time Division; S4. Generate spacetime blocks and a mixed-integer programming model for spacetime blocks; S5. Establish a master-slave iteration framework; S6. Iterative solution process.

2. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S1, data acquisition and preprocessing, specifically includes: S11. Acquisition and preprocessing of annual load forecast data; The power system dispatch center periodically forecasts the subsequent electricity demand of the power grid within its dispatch area and stores it in a database. The load forecast data includes a timestamp and the corresponding load value, and is processed as follows: First, a data quality check is performed to identify and handle outliers. A statistical method is used to calculate the Z-score for each hourly load forecast: Z = (x - μ) / σ, where x is the original load value, μ is the historical average load for the same period, and σ is the standard deviation. Values ​​with |Z| > 3 are considered outliers and are replaced with the average load of the two preceding and following hours. Then, data normalization is performed, mapping the load values ​​to the [0,1] interval. The Min-Max normalization method is used: Where L is the original load value, and These are the minimum and maximum load values ​​for the entire year, respectively. S12. Acquisition and standardization of unit parameters; Obtain parameters for various thermal power units and hydropower units, convert all unit parameters to per-unit values, and use the system baseline capacity as the benchmark: = / ,in Per unit value, This represents the actual active power of the generator unit; S13. Extraction and structuring of network constraint data; The network topology and operational constraint data of the power grid are obtained from the power grid energy management system (EMS), including node data, branch data, and transformer data. This data is then used to construct a node-branch correlation matrix and an admittance matrix. S14. Data consistency verification and time alignment; Finally, perform consistency checks across data sources: ensure all data uses the same timestamp system; verify the match between the total unit capacity and the system's maximum load; and check the correspondence between network constraints and unit locations. After completing the above steps, the three types of data are encapsulated into structured data objects in a unified format, providing standardized input for subsequent spatiotemporal decomposition steps.

3. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S2, spatial partitioning, specifically includes: S21. Electrical distance calculation and similarity matrix construction; First, calculate the electrical distance between all nodes in the system as the basis for network partitioning; S22. Power Grid Diagram Model Construction and Community Discovery; The power grid is abstracted as a graph model G=(V,E,W), where: V is the set of nodes, containing all the main nodes; E is the set of edges, where each edge represents the connection between nodes; W is the edge weight matrix. ; Network partitioning is performed based on the power grid diagram model; S23. Zoning quality assessment and adjustment; S24. Boundary node identification and tie-line processing; After the space is divided, the boundary nodes and connecting lines between each area are identified. Boundary nodes are nodes that are electrically connected to multiple areas at the same time, and connecting lines are power transmission lines that connect different areas. S25. Final partitioning scheme output; The final output zoning scheme includes: (1) a regional division list: the nodes, units, and load information contained in each region; (2) a boundary node mapping table: recording the multiple regions to which each boundary node belongs; (3) tie line parameters: including impedance limit and thermal stability limit; (4) a regional characteristic summary: the power generation capacity, load level, and network density of each region.

4. The simulated clearing method for the electricity spot market according to claim 3, characterized in that, Step S21. Electrical distance calculation and similarity matrix construction, using an electrical distance calculation method based on the impedance matrix: Define the electrical distance matrix D, whose elements The electrical distance between node i and node j is represented by the following formula: ， in: Let be the voltage phase angle (in radians) at node i. Let be the voltage amplitude at node i. The system's rated voltage. The transfer impedance between nodes i and j These are weighting coefficients used to balance the effects of voltage differences and impedance. The α value is automatically adjusted according to the system's operating status, decreasing under heavy load and increasing under light load.

5. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S3, time division, specifically includes: S31. Load feature extraction; S32. Dynamic Programming: State Definition and State Transition; The goal of dynamic programming is to divide the time into several time blocks, so that the load fluctuation of each time block is small, the cost estimate is stable, and a balance is achieved between reducing computational complexity and minimizing accuracy loss. S33. Dynamic programming solution and post-optimization adjustment; The reverse dynamic programming algorithm is used to solve the problem recursively. After the solution is completed, a preliminary time block partitioning scheme is obtained by backtracking.

6. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S4 generates a spacetime block and a mixed-integer programming model for the spacetime block, specifically including: S41. Spatiotemporal block generation and data allocation; Based on the M region division results obtained in step S2 and the K time block division results obtained in step S3, M×K spatiotemporal blocks are generated, and each spatiotemporal block corresponds to a sub-problem of a region in a specific time period. Assign corresponding input data to each spatiotemporal block, including load data of all nodes in the region within the time block, extract technical parameters of all units in the region, network constraints and boundary connection constraints within the region; S42. Construction of a mixed-integer programming model; For each spatiotemporal block, construct a mixed-integer programming model of the following form: The objective function is: , in, Let be the power generation cost function of unit j; , These represent the start-up and shutdown costs of unit j, respectively. , This represents the start / stop status of unit j, with a value of 1 or 0. The power generation capacity of unit j at time t; The time range of the spatiotemporal block is given by t, where t is the hour index. A collection of units within a spacetime block; The constraints include: power balance constraints, unit operation constraints, network security constraints, and time boundary consistency constraints.

7. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S5 establishes the master-slave iteration framework, specifically including: S51. Master-slave architecture design; The master-slave iterative framework is a two-layer optimization framework, including a coordination layer for the master problem and a computation layer for the subproblems. The coordination layer is responsible for coordinating the boundary conditions between subproblems in various regions, maintaining global consistency constraints, ensuring the balance of power exchange between regions, and generating and updating coordination signals. The computation layer is used to solve the mixed integer programming model of each spatiotemporal block in parallel, submit local optimization results to the master problem, and receive and respond to the coordination signals issued by the master problem. S52. Modeling the main problem; The main problem is to use price signals as a coordinating tool to generate coordinating signals; S53. Sub-problem coordination mechanism; Each subproblem receives a coordination signal from the main problem during the solution process and incorporates it into its local optimization objective. After the subproblem is solved, it submits data to the main problem.

8. The simulated clearing method for the electricity spot market according to claim 1, characterized in that, Step S6, the iterative solution process, specifically includes: S61. Iterative solution process; The alternating direction multiplier method framework is used for master-slave collaborative solution. The specific process is as follows: Set the initial iteration counter k=0 and initialize the price signal. and At the same time, set the maximum number of iterations and the convergence tolerance; Then, each computing node solves its assigned spatiotemporal block mixed integer programming model in parallel, taking into account the coordination signal provided by the main problem during the solution process; after the management node collects the solution results of all subproblems, the main problem updates the price signal and calculates the degree of global constraint violation. If the degree of constraint violation is minimal and the convergence condition is met, the solution is considered complete; otherwise, an updated coordination signal is issued for the next round of calculation. After the solution is completed, a weighted average method is used to ensure consistency for boundary variables; S62. Clearing price calculation; Calculate the nodal marginal price (LMP) based on the final optimization results: ， in, The system energy price component is equal to the dual variable of the global equilibrium constraint; The blocking price component of node i is equal to the sum of the dual variables of the corresponding network constraints; The network loss price component for node i is calculated using the network loss distribution factor. S63. Output report compilation; The final output report includes the unit scheduling plan, clearing price curve, and ancillary information.

9. A simulated clearing system for the electricity spot market, characterized in that: The data acquisition and preprocessing module is used to perform the above step S1; The spatial partitioning module is used to perform step S2 described above; The time division module is used to execute step S3 above; The spatiotemporal block generation and mixed-integer programming module is used to perform step S4 above; The master-slave iteration framework establishment module is used to execute step S5 above; The iterative solution module is used to execute step S6 as described above.

10. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the invention's simulated clearing method for the electricity spot market.

Images