A virtual power plant optimization scheduling method based on physical-electric decoupling

Abstract

The application discloses a kind of virtual power plant optimization scheduling methods based on matter electricity decoupling, it is related to power system optimization scheduling technical field, including the following steps: constructing each time period node injection power feasible region and cross time period timing constraint;Using Lagrange relaxation to convert timing constraint into dual function, introduce the regularized penalty term based on resource physical response time constant to carry out smooth approximation, calculate gradient correction term;The correction term is updated to power feasible region boundary, and the corrected feasible region sequence is obtained;In the sequence, the power instruction mapping relationship about prediction error is established, and then the reference output plan and power output correction instruction are output.The present application is used to solve the problem of virtual power plant scheduling model dimension disaster, insufficient robustness and poor decoupling ability of heterogeneous resource matter electricity characteristics in the prior art.

Description

Technical Field

[0001] This invention relates to the field of power system dispatch automation technology, and more specifically, to a virtual power plant optimization dispatch method based on physical-electrical decoupling. Background Technology

[0002] In existing virtual power plant dispatch architectures, the spatial constraints of distribution network nodes and the temporal constraints of internal heterogeneous resources (such as ramp-up and energy storage state transitions) are usually directly coupled into a high-dimensional model for centralized solution. This approach faces two major computational and control bottlenecks in practical engineering: First, the extremely high spatiotemporal coupling dimension leads to a significant increase in computational load, making it difficult to meet the computational requirements for real-time command correction; second, centralized solutions often ignore the significant differences in the physical response time constants of various regulating resources. Due to the lack of an effective damping feedback mechanism, control command oscillations are easily triggered at the constraint boundaries, thereby exacerbating the fatigue and wear of underlying equipment.

[0003] Furthermore, traditional dispatching relies heavily on deterministic baseline plans and lacks the dynamic adaptability to random fluctuations in renewable energy output and forecast errors. Once real-time operating conditions deviate from the day-ahead forecast, the system struggles to smoothly distribute and adjust errors among multiple heterogeneous resources, making it easy for real-time correction commands to exceed the physical response boundaries of equipment, thereby posing a risk to the safe operation of the distribution network. Summary of the Invention

[0004] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a virtual power plant optimization scheduling method based on physical-electrical decoupling, in order to solve the problems of dimensionality curse, insufficient robustness, and poor decoupling capability of heterogeneous resource physical-electrical characteristics in the virtual power plant scheduling model in the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: A virtual power plant optimization scheduling method based on physical-electrical decoupling includes the following steps: A feasible region for node injection power in each time period is constructed, and cross-time period ramping and state transition constraints are established to form temporal constraints. The temporal constraints are transformed into a dual function using Lagrange relaxation, and a regularization penalty term based on the physical response time constant of adjustable resources is introduced. A smooth approximation is performed on the dual function, and a gradient correction term is calculated. The gradient correction term is updated to the boundary of the node injection power feasible region to obtain a corrected feasible region sequence. The physical response time constant is defined as the shortest time required for an adjustable resource to ramp up from zero power to full power at its maximum ramp rate, measured in units of the number of scheduling periods. Within the corrected feasible region sequence, a power command mapping relationship regarding prediction error is established, and a baseline output plan and power output correction command are output based on this mapping relationship.

[0006] The technical effects and advantages of the virtual power plant optimization scheduling method based on physical-electrical decoupling of the present invention are as follows: First, it reduces the solution dimensionality and improves real-time computation efficiency. By constructing a feasible region for node-injected power and Lagrange relaxation, complex cross-time-period physical constraints are decoupled in the time dimension. This "physical-electrical decoupling" mechanism effectively overcomes the surge in computational load caused by high-dimensional strong coupling and meets the computational power requirements for real-time system correction.

[0007] Second, eliminate control oscillations and extend equipment lifespan. The physical response time constant of the equipment is incorporated into the smooth approximation penalty term. By dynamically updating the boundary through gradient correction terms, frequent oscillations of control commands near the physical boundary are eliminated, achieving smooth control of heterogeneous resources and reducing equipment wear.

[0008] Third, enhance system robustness and ensure the physical safety of the power grid. Within the corrected feasible domain sequence, establish a power command mapping relationship for prediction errors. When real-time renewable energy output deviates, a safe output correction command can be directly output without requiring global re-optimization, strictly preventing equipment from exceeding limits. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of a virtual power plant optimization scheduling method based on physical-electrical decoupling, provided for an embodiment of the present invention.

[0010] Figure 2 This is a comparison chart of the convergence trajectories of the dual variables before and after correction. Detailed Implementation

[0011] 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.

[0012] Example 1, Figure 1 This invention presents a virtual power plant optimization scheduling method based on physical-electrical decoupling, comprising the following steps: S1, construct the feasible domain of injected power at each time period node, and establish cross-time period ramping and state transition constraints to form time-series constraints; S2, using Lagrange relaxation to transform the temporal constraints into a dual function, introducing a regularization penalty term based on the physical response time constant of the adjustable resource, performing a smooth approximation on the dual function, and calculating a gradient correction term; updating the gradient correction term to the feasible region boundary of the node injection power to obtain the corrected feasible region sequence; wherein, the physical response time constant is defined as the shortest time required for the adjustable resource to rise from zero power to full power at its maximum ramp rate, with the number of scheduling periods as the unit; S3. Within the corrected feasible domain sequence, establish a power command mapping relationship for the prediction error, and output a baseline output plan and a power output correction command based on the power command mapping relationship.

[0013] It should be noted that the "physical-electrical decoupling" described in this invention refers to the process of using Lagrange relaxation and function smoothing approximation techniques to extract the characteristics of the time-series constraints that originally belong to global nonlinear coupling from the complex physical state evolution equations, and mapping them to dynamic boundary correction formulas for the feasible domain of injected power at each time period. This ensures the efficiency and robustness of control while achieving decoupling and coordination between the dynamic constraints of the physical layer (physical) and the electrical optimization of the scheduling layer (electrical).

[0014] The specific implementation process of S1 includes: S101, Construct the feasible domain of node injection power for each time period; S1011, Obtain the distribution network security operation domain. Specifically: Performing second-order cone relaxation on the non-convex AC power flow equations yields a set of convex inequalities. Specifically, based on the original AC power flow equations (DistFlow branch power equations), performing second-order cone relaxation results in a set of convex inequalities composed of second-order cones (Lorentz cones) plus the linear constraints in the network: (1) In the formula: For scheduling periods; Number the nodes; For the active and reactive power flow of the branch; The square of the voltage magnitude at the node; The square of the current amplitude of the branch is given.

[0015] It should be noted that this application is applicable to radial or weak-loop distribution networks; in radial networks, the second-order cone relaxation is exactly equivalent (single-shot), and weak-loop networks can be supplemented with loop interruption or loop constraint processing to ensure that the security domain modeling is not distorted.

[0016] The convex inequality set is approximated by a polyhedron using a piecewise linearization method to obtain the linearized constraint set. The piecewise density is dynamically adjusted according to the rate of change of the load prediction curve. The tuning calculation of the piecewise density (number of segments) is shown in Equation (2). (2) In formula (2): Indicates time period Segment density (number of segments); Indicates the preset minimum number of segments (preferably) ); This represents the preset gradient sensitivity gain constant, and its unit is MW. - ¹·Time period, making The increment for the dimensionless segment number is preferably in the range of [0.5, 2.0]. and These represent the day-ahead load forecast curves in the time period. and The predicted active power load value, where the load forecast curve is obtained using a conventional load forecasting algorithm, is represented as follows: ; Indicates the time interval step size; This represents a preset basic segmentation constant (preferably 3); This represents the function for rounding up.

[0017] The linearized constraint set can be represented in matrix form as follows: (3) In the formula: the total node injected power vector can be divided into the grid-connected part and the non-grid-connected part. ,in Inject a power vector into the grid connection point (the first component is active power, and the second component is reactive power). Inject power into non-grid-connected nodes. This represents the injected power vector (containing the active power component and reactive power component of the node) of each node in the distribution network. This represents a vector of internal state variables containing the power flow, voltage, and current squares of the entire network. Represents the node injection power mapping matrix; Represents the power flow-voltage linear coupling matrix; This represents the boundary limit vector for safe operation of the distribution network (including line thermal stability capacity limits and node voltage deviation limits). It should be noted that the node-injected power mapping matrix is ​​the correlation matrix between nodes and branches, used to characterize the linearized mapping of Kirchhoff's current law in the matrix domain; the power flow-voltage linear coupling matrix is ​​obtained by applying the branch voltage drop equation... The linear coefficient matrix obtained by rearranging (i.e., the linear DistFlow approximation); the safe operation boundary limit vector is a fixed value vector preset offline, whose dimension is determined by the total number of distribution network branches and the total number of nodes, defining the feasible solution space of the distribution network under the power grid safety constraints (N-0 criterion).

[0018] Power injection vector at grid connection point To preserve variables, Fourier-Motzkin elimination is performed on the linearized constraint set (Equation (3)) to eliminate internal variables. The safe operating domain of the distribution network is obtained, as shown in equation (4): (4) In equation (4): Indicates time period The projected safe operating domain of the power distribution network; and These represent the projected safe operating boundary coefficient matrix (dimensions). ) and the dimension of the limit vector It consists of the set of explicit polyhedral normal vectors and intercepts retained after eliminating internal variables using Fourier-Motzkin elimination. Since the retained variables are only two-dimensional (active and reactive), the elimination does not produce a dimensionality explosion.

[0019] To further reduce computational load, implicit projection algorithms based on Benders decomposition or key vertex search methods can be used in engineering implementation to directly construct approximate values. and This reduces the computational dimensionality while ensuring the safety margin of the power grid.

[0020] S1012, Obtain the feasible output domain of each adjustable resource; Extract the output boundary parameters of each feasible output region and construct a linear inequality-constrained polyhedron. Specifically, for the M adjustable resources within the virtual power plant, extract the output boundary parameters of each adjustable resource. For the ... Each resource has a corresponding feasible output domain constructed from a set of linear inequalities as shown in equation (5): (5) In formula (5): and Representing resources The effective and ineffective contributions; and These represent the upper and lower physical limits of active power output, respectively. and These represent the upper and lower physical limits of reactive power output, respectively.

[0021] Rewrite equation (5) in matrix form: , in In this embodiment, the following are provided: .

[0022] Each of the polyhedra undergoes a Minkowski Sum aggregation operation to generate an initial running domain, as shown in equation (6): (6) In formula (6): This represents the initial operating domain generated by the aggregation (defined in the resource output space); Minkowski sum operations are represented. This represents the total number of adjustable resources participating in aggregation within the virtual power plant; Indicates the number index of the adjustable resource; Indicates the first An adjustable resource during the time period The corresponding initial operating domain.

[0023] The injected power at the grid connection point of the virtual power plant is the sum of the outputs of all resources, and a resource output stacking vector is defined. , With aggregation mapping matrix : (7) The initial running domain Through aggregation mapping matrix Projecting onto the grid connection point space yields the feasible region for grid connection point aggregation. : (8) In the formula, These are the boundary coefficient matrix and limit vector of the aggregation domain at the grid connection point, respectively, given by equation (9).

[0024] Since the output constraints of each resource are independent boxes, equation (8) degenerates into a rectangle (ultra-cubic prism) in the grid connection point space: (9) The secure operating domain The intersection of the feasible region of the grid connection point and the tangent plane algorithm is used for iterative optimization to generate an inscribed polyhedron as the feasible region for node injection power. As shown in equation (10): (10) In formula (10): Indicates time period The feasible region of node injection power (i.e., the feasible region of grid connection point injection power). and Let represent the boundary coefficient matrix and limit vector of the inscribed polyhedron constraint, respectively; This represents the total number of valid constraint hyperplanes determined after iterative optimization. It should be noted that each row of the boundary coefficient matrix of the inscribed polyhedron constraint... For the first The normal vectors of the effective tangent planes ( ), which is from the projected safety boundary Aggregation boundary with grid connection point The normal vectors are normalized and concentrated, and obtained through iterative selection using the tangent plane algorithm, satisfying the following conditions: : (11) Limit vector Each component is obtained by taking the support function along the normal direction in the safe operation domain and the grid connection point aggregation domain respectively (i.e., along...). Taking the smaller of the upper bounds of the directions, we get: (12) (13) In equations (12) and (13): For the first Limits (scalars) of the hyperplane; minimizing the support functions of the two domains. Just the intersection along An upper bound for the directional support function, minus the nonnegative convergence margin. This ensures that the resulting half-space does not overflow the intersection, thus achieving inscribedness; For the tangent plane algorithm in the first The convergence margin on each hyperplane is updated by equation (14): (14) In the formula, Let the relaxation step size (convergence coefficient) be taken as... Control the convergence margin Convergence rate; This step, while ensuring that the obtained polyhedra are completely inscribed within the original safe intersection (ensuring the absolute physical safety of the power grid), significantly reduces the constraint dimension and effectively avoids the real-time optimization delay caused by a massive number of linear inequalities.

[0025] This invention effectively solves the problem that traditional high-dimensional constraint models cannot match in real time when dealing with complex time-series climbing by using a dimensionality reduction mapping strategy, thus ensuring the physical feasibility of scheduling instructions throughout the entire time period.

[0026] S102, establish cross-time-period ramp-up and state transition constraints to form temporal constraints; No. The cross-time ramping constraint for each resource is shown in equation (15): (15) In equation (15): and These represent the maximum upward and downward ramp limits of the equipment (unit: MW / time period).

[0027] When the When a resource is an energy storage type, the energy storage state transition constraints, such as the decoupled distribution relationship of active power output in adjacent time periods, the energy state transition equation across time periods, the capacity boundary constraints, and the energy balance constraints at the beginning and end of the daily cycle, are shown in equations (16) to (19): (16) (17) (18) (19) In equations (16) to (19): and These represent the active power of the energy storage unit during charging and discharging, respectively. This represents the charged energy of the energy storage unit (unit: MWh). and These represent the charge-discharge conversion efficiency, respectively. and These represent the upper and lower limits of energy storage operating capacity (unit: MWh); This indicates the total number of time periods in the day-ahead scheduling cycle; This represents the initial charged energy.

[0028] To eliminate simultaneous charging and discharging, complementary relaxation constraints are further introduced. Alternatively, a binary variable can be introduced to ensure that at most one of the charging and discharging powers is positive.

[0029] The timing constraints formed by equations (15) and (16) to (19) can be uniformly represented as a compact matrix form: (20) In the formula: This is a stacked vector of decision variables for the entire time period, where the output variable affects all... All resources exist, and the variable of charged energy. Only for energy storage resource sets It exists, and its dimension is denoted as . ; The temporal constraint correlation matrix (block-like, with each row coupling variables from at most two adjacent time periods); The limit vector (composed of the physical ramp limits of each resource, the physical extreme values ​​of energy storage capacity, and the energy consistency requirements at the beginning and end of the daily cycle); the total number of time-series constraints is denoted as... , equal to the correlation matrix The number of rows (number of coupling constraints) is the dimension of the Lagrange multiplier.

[0030] The specific implementation process of S2 includes: S201, the temporal constraints are transformed into dual functions using Lagrange relaxation, and a penalty term is introduced to perform smooth approximation, and the gradient correction term is calculated. S2011, The time series constraints are transformed into dual functions using Lagrange relaxation; Establish the matrix form of timing constraints Introducing global Lagrange multiplier vectors ( The scheduling objective function is to minimize the total operating cost of all adjustable resources, and is expressed as: (twenty one) In the formula, This indicates the total number of time periods in the day-ahead scheduling cycle.

[0031] The corresponding Lagrange function is: (twenty two) In the formula: The global temporal constraint Lagrange multiplier vector (dimension equal to) (The row numbers are based on the coupling constraints, not on a single time period). They represent the first The economic cost coefficients of the quadratic, linear, and constant terms of a resource (determined by equipment depreciation, fuel costs, or operation and maintenance costs).

[0032] because Given a block-band matrix (each row couples at most two adjacent time periods), Inner layer minimization It can be decomposed into independent subproblems of each resource in the time dimension, and its dual function is defined as: (twenty three) This step transforms the strong physical coupling across time periods into an economic penalty through Lagrange relaxation, achieving mathematical decoupling in the time dimension.

[0033] S2012, A regularization-based penalty term is introduced to perform a smooth approximation of the dual function, generating a smooth dual function; wherein the regularization coefficient is determined by the physical response time constant of the adjustable resource; Due to dual function Given a piecewise linear concave function with non-differentiable points, a quadratic regularization term is introduced. Performing a smooth approximation on it (equivalent to applying Moreau / proximal type regularization to the multiplier space, which can also be regarded as a direct implementation of Moreau-Yosida type smoothness) generates a smooth dual function. : (twenty four) Wherein, the regularization coefficient matrix It is a diagonal matrix, consisting of the physical response time constants of the resources corresponding to each timing constraint. Sure: (25) (26) In the formula: The smoothness reference scaling factor (preferably 1.0); For the first The index of the resource bound by each time-series constraint. ; For the first The physical response time constant (unit: number of scheduling periods) of the resource bound to each timing constraint. For resources The upper limit of effective output, , Set its ramp limit. Response time constant. The smaller (the faster the resources). The larger the value, the stronger the smoothness, thus exerting stronger inertial damping on fast resources.

[0034] It should be noted that the physical response time constant mentioned in equation (26) Defined as a response time characteristic parameter of adjustable resources on a scheduling timescale, it is measured in units of scheduling periods and characterizes the shortest time required for the resource to ramp up from zero power to full power at the maximum ramp rate, distinct from the electromechanical time constant in the electrical transient sense. This definition allows for the quantification and comparison of response speed differences among heterogeneous resources in virtual power plant scheduling. A smaller value indicates a faster resource response time, and the corresponding regularization coefficient is... The larger the value, the stronger the smoothing effect, thus applying stronger inertial damping to fast resources.

[0035] By introducing this regularization term, the non-smooth dual space is transformed into a continuously differentiable strongly concave function space, effectively eliminating the "start-stop control" generated near the optimal boundary by the traditional subgradient optimization method, making the update of the dual variable more stable and possessing physical inertial damping characteristics.

[0036] like Figure 2 The figure shown is a comparison of the convergence trajectories of the dual variables before and after gradient correction according to the present invention: the horizontal axis represents the number of iterations, and the vertical axis represents the dual variables. The representative component is represented by the dashed line, which is the convergence trajectory of the traditional subgradient optimization method without introducing a physical response time constant penalty term. It exhibits constant amplitude oscillations (start-stop control) near the optimal dual boundary. The solid line is the convergence trajectory after smooth approximation, which shows smooth and asymptotic convergence, thus demonstrating the technical effect of the present invention in suppressing command jitter.

[0037] S2013, Solve for the gradient vector of the smooth dual function, as the gradient correction term; According to the envelope theorem, the gradient of the smooth dual function is the Lagrange function with respect to the Lagrange multipliers. The partial derivatives. Taking the derivative of equation (24) yields the gradient vector (gradient correction term). : (27) In the formula: For Lagrange multipliers The current optimal decision vector is obtained by solving the subproblem (Equation (24)); Given the current physical resource response, the real-time violation residual of the timing constraints (i.e., the degree of timing constraint violation). This is the corrected bias generated by the regularization term.

[0038] During the iterative solution process, this gradient correction term serves as a virtual feedback signal, mapped to the underlying resource control system. When a component of the gradient vector... When this occurs, it indicates a severe violation of the corresponding temporal constraints, and the gradient correction term will cause the dual variable in the next iteration to... This increases the penalty for violating timing constraints, forcing the scheduling plan to move in a direction that satisfies the laws of physical ramping and energy transfer.

[0039] S202, update it to the feasible region boundary of node injection power to obtain the corrected feasible region sequence; S2021, establish a sliding buffer for the sequence of differences between dual variables and gradient vector differences, wherein the sliding window step size is dynamically set according to the ratio of the fastest and slowest physical response time constants of the adjustable resources. Specifically, establish the difference between dual variables. Difference with gradient vector sliding buffer : (28) in, Indicates the first Global Lagrange multiplier vector at each iteration step; sliding window step size Based on the fastest response time constant of the resource With the slowest response time constant The ratio is dynamically set: (29) In the formula: The historical information retention coefficient for the sliding window is preferably 0.25. , These are the minimum and maximum values ​​of the physical response time constants of all resources, respectively.

[0040] S2022, recursively calculate the sliding buffer according to the rank-2 update rule in the quasi-Newton method (BFGS) to construct an approximate matrix of the Hessian inverse matrix. ; (30) In the formula: It is the identity matrix; For curvature correction scalar, Let be the curvature gradient difference vector. The dual problem of this invention is the maximization of a concave function, ensuring that the maximization of the concave function... To maintain Positive definiteness, damped BFGS prevents pathological conditions. If If so, the failure protection for this rank-2 update will be skipped.

[0041] S2023, Perform matrix multiplication between the approximate matrix and the current gradient vector to determine the update direction of the dual variable. : (31) S2024, Calculate the obstacle correction vector of the dual boundary in the timing constraint; Specifically, the obstacle correction vector The logarithmic barrier gradient (dimensions / MW) defined in the dual space characterizes the compactness of the physical constraints approximating the dual variables (i.e., The repulsive potential energy at that time is used to suppress Crossing this threshold leads to numerical overflow, which, together with the BFGS direction, forms a stable update direction: Let the real-time relaxation margin of the m-th timing constraint be... (The units of measurement are constrained units, such as MW), define the logarithmic barrier function of the dual space: (32) Among them, barrier parameters Dimensions (i.e., the multiplier) The square of the dimension plays a role in scaling and dimensionless adjustment. The obstacle correction vector is taken with respect to its... The total derivative (gradient in dual space): (33) In the formula: and Global timing constraints No. row (i.e., the first row) The coefficient vector of each time-series constraint and the limit constant (constraint is...) ); For the current dual variable The optimal solution to the original problem (Equation (24)); For the first The real-time relaxation margin of each time-series constraint (positive when the Lagrange relaxation solution is in the feasible neighborhood of the constraint; the original solution of the Lagrange relaxation may have slight violations, so this barrier term is used as a stability correction), when When the constraint boundary is increased to make the solution of the original problem approximate the boundary condition... The repulsive potential energy tends to infinity, thus repelling... Exceed the critical threshold and avoid numerical overflow.

[0042] Sensitivity It can be obtained from the dual Hessian approximated by BFGS through the KKT adjoint equation or estimated by finite difference, without the need to build a separate solver for the obstacle term.

[0043] It should be noted that the equation (33) gives yes Regarding dual variables The dual space gradient, rather than the old one. The primitive sensitivity, which contains only the inverse of the relaxation margin, is also not a smooth dual function. right The analytical total derivative; its dimensions, after being bridged by the chain rule, are related to the BFGS direction. (Equation (31), the dimension is also / MW) are completely consistent, so the two can be directly weighted and merged according to Equation (37) without any contradiction in the dimension.

[0044] S2025, Calculate the margin from the current dual variable to the temporal constraint boundary, whereby the margin is defined as the ratio of the remaining settling time of the ramp constraint to the total settling time; (34) in The total adjustable time corresponding to the climbing constraint. This represents the adjustable time that has been consumed up to now.

[0045] S2026, Calculate the mapping weight coefficients based on the rate of change of the margin; Based on the rate of change of margin : (35) Calculate the mapping weight coefficients : (36) In the formula is the weighted smoothing constant.

[0046] S2027, according to the weight coefficient, perform linear weighted fusion on the update direction and the obstacle correction vector, update the dual variable, and substitute the updated dual variable into the boundary correction vector to obtain the corrected feasible domain sequence; According to the weight Regarding the direction of updates With obstacle correction vector Perform linear weighted fusion and update the dual variable. : (37) In the formula The iteration step size is set to satisfy the Wolfe condition; after the update, the projection is forced to the non-negative quadrant: .

[0047] Updated dual variables Substitute the boundary correction vector .remember This belongs to a time period The set of temporal constraint subscripts (including Climbing constraints and energy storage (Government time period energy / capacity constraints), then: (38) (39) (40) (41) In the formula: For the boundary correction vector, its components correspond to The translation distance of each constrained plane; For time period The translation direction vector of the aggregated virtual boundary; For the first The sensitivity vector of a time-series constraint in the grid connection point space is specifically formed by the aggregation mapping of the normal vector of the constraint. The resources obtained by mapping to the grid connection point space Active ramp constraint correspondence (Active power axis direction at grid connection point), corresponding reactive power constraints (Reactive axis direction at grid connection point), energy storage energy / capacity constraints, and its grid connection point spatial sensitivity. ; The normalized weighting coefficient for resource response speed ( (The sensitivity index is preferably 1.0). For the first The dual variables after the time-series constraints are updated.

[0048] Boundary matrix of the feasible region for injecting power into nodes and limit vector The updated feasible region sequence is obtained by performing an update. : (42) The corrected feasible region sequence By the boundary intercept Adaptive translation enables real-time tracking and dynamic response of the scheduling space to temporal physical constraints.

[0049] The specific implementation process of S3 includes: S301, centered on the empirical distribution of historical prediction errors, constructs a fuzzy set of the probability distribution of the true error based on the Wasserstein distance, where the order of the Wasserstein distance is... Its radius is dynamically and adaptively set based on the time autocorrelation decay coefficient of the prediction error. The fuzzy set of the probability distribution is: (43) in, The distance is the second-order Wasserstein distance; Based on historical prediction errors and empirical observations Empirical distribution ( Sourced from the historical database of virtual power plants SCADA actual records for one calendar day, preferred ); Distance radius Based on the time autocorrelation decay characteristics of the prediction error sequence, the settings are dynamically and adaptively adjusted according to the look-ahead time interval: (44) In the formula: The base radius (calibrated offline based on the dispersion of the historical error distribution); For the current decision-making moment, For time period The timeframe relative to the current decision-making moment; The autocorrelation coefficient of the prediction error sequence AR(1) Derived uncertainty growth rate (autocorrelation at time intervals) The uncertainty radius decreases with time, therefore the uncertainty radius decreases with time. The greater the time interval, the greater the prediction uncertainty. The upper limit of the radius (calibrated offline by the maximum tolerable error within the cycle) is set to prevent excessive expansion of the remote time distance radius during the scheduling period. hour, Pick To maintain dimensional consistency.

[0050] S302, Perform a strong dual transformation on the expected cost objective function containing the fuzzy set, and transform it into a set of linear matrix inequalities (LMI) constraints; Construct the expected cost objective function containing fuzzy sets : (45) In the formula: For the first Each resource during the time period The baseline power output plan (in MW) is a deterministic decision variable for the day-ahead scheduling phase; For the first Each resource during the time period The real-time adjustment allocation coefficient (dimensionless) characterizes the proportion of the resource's response to fluctuations in new energy sources. For the first The real-time prediction error value (in MW, scalar) of renewable energy output during a given period follows a probability distribution as a random variable. ; The penalty factor for the physical lifespan depreciation of equipment due to high-frequency regulation (data sourced from the nameplate parameters of the physical equipment). For the first The basic generation / operation cost function for each resource (fitted and calibrated by historical fuel consumption or grid purchase and sale contract electricity price); For resources During the period The high-frequency regulating power (i.e., the real-time deviation of the actual output power relative to the reference), its square term Characterize the second-order penalty that increases adjustment costs as deviation increases.

[0051] Considering that the virtual power plant's baseline output plan must be balanced with the projected net load, a power balance constraint is established: (46) In the formula For time period The projected net load; its Lagrange multiplier is denoted as (Scalar, due to a single balance equation for each time period). The power deviation in the real-time phase is constrained by the allocation coefficient of equation (49). Resource-by-resource allocation is used to increase the actual grid-connected power. Real-time tracking of the actual net load, including errors.

[0052] By utilizing Kantorovich-Rubinstein strong duality and semidefinite programming reconstruction (Schur's complement lemma), the complex probabilistic analytical problem is transformed into a set of finite-dimensional linear matrix inequalities (LMI) constraints. For convex quadratic loss under second-order Wasserstein fuzzy sets, the semidefinite programming reconstruction of its worst expectation is based on the conventional dual theory of Wasserstein partial Bruker optimization. A second-order expansion of the cost function is performed near the historical error samples, retaining the convex quadratic penalty term. And record its details. The first-order fitting coefficient is The constant term is The quadratic coefficient is For each historical sample LMI takes the form of: (47) In the formula: For Wasserstein scalar dual variables, represent the weight of the influence of the fuzzy set radius on the objective function; The Lagrange multiplier (scalar) for power balance constraint is a deterministic balance constraint multiplier, independent of this LMI, so it does not enter equation (46); The auxiliary dual variable introduced (the first) (each sample corresponds to) This indicates that The matrix is ​​a positive semi-definite matrix (because) It is a scalar, and the block dimension is ); , , The scheduling objective function is respectively in the prediction error The first-order coefficients, constant term, and second-order coefficients of the nearby expansion; For the first One historical prediction error empirical observation.

[0053] It is worth noting that retaining the quadratic penalty term avoids the degradation of the optimal solution caused by excessive linearization (i.e., start-stop control of control commands); this quadratic structure is the mathematical basis for applying Schur's complement lemma to project it into a set of linear matrix inequalities (LMI) constraints, thereby eliminating the obstacles of nonlinear calculation while retaining the system's second-order damped response characteristics to meteorological error fluctuations and the ability to smoothly allocate adjustment costs.

[0054] S303, the set of linear matrix inequality constraints is used as a deterministic penalty term and included in the scheduling objective function for equivalent solution; Introducing scalar auxiliary variables The upper bound of the "maximum potential operational risk cost of the system" under the worst weather conditions is represented, which satisfies the following linear inequality constraint (the penalty term under the second-order Wasserstein fuzzy set is proportional to the square of the radius): (48) The scheduling objective function, taking into account deterministic penalties, is: (49) in, The robustness adjustment factor (preferably 1.0) is set by the power grid dispatcher based on the extreme weather warning level of the day.

[0055] In solving equation (49), the decision variables , The following conditions must be met simultaneously: ① The feasible region of node injection power defined by equation (10) after being modified by equation (41) ; ② Cross-period ramping and energy storage constraints of equations (16) to (19); ③ Power balance of equation (45); ④ LMI constraint of equation (46) and risk upper bound constraint of equation (47); ⑤ Allocation coefficient constraints of S304 and S305 below. This ensures that the baseline output plan is physically feasible throughout the time period and has the ability to defend against partial brooms.

[0056] S304 sets the output command for each time period as a power command mapping relationship between the baseline output plan and the real-time adjustment allocation coefficient, and applies a balance constraint that the sum of each resource allocation coefficient is one. Output command Set as the baseline output plan With real-time adjustment of allocation coefficient Power command mapping relationship: (50) S305, substitute the power command mapping relationship into the feasible region of each time period to obtain the baseline output plan. and allocation coefficient Solving for joint decision variables, where the response capability boundary of the allocation coefficients is determined based on the magnitude ratio of the gradient correction term; Specifically, the power command mapping relationship is substituted into the corrected feasible domain sequence. In China, with Solving for joint decision variables. Based on the gradient correction term magnitude ratio, defining the target share of each resource. Boundary of controlled deviation: (51) (52) In the formula: For gradient correction term vector According to resources Subvectors formed by corresponding components (i.e., constraint multipliers) Binding resources The sum of the gradient components corresponding to several components), its magnitude ratio Reflecting resources The proportion of physical capacity to perform regulatory tasks, and ; For controlled deviations around the target share; because and ,constraint Automatic fulfillment; ( The bias bandwidth coefficient (preferably 0.5) is determined by the gradient magnitude ratio. The defined "response capability boundary" ensures that resources with faster response times undertake a larger share of regulation, while avoiding the occurrence of regulation dead zones for some resources due to excessively small gradient correction term magnitudes. The minimum proportion coefficient is preferably 0.05.

[0057] S306, set the depreciation cost coefficient for resource response under different frequency fluctuations as a frequency domain penalty; perform spectral analysis on the prediction error sequence and separate high and low frequency components, wherein the cutoff frequency is determined by the boundary point of the response time constant of each resource; Specifically, spectral analysis is performed on the historical prediction error sequence, using the threshold of the boundary time constant. Separate the error into high-frequency components. with low-frequency components ; Set the depreciation cost factor as the frequency domain penalty Its differential tuning is as follows: (53) In the formula For resource response time constant, The cutoff frequency is the threshold separating the response time constants of power-type and energy-type devices (i.e., the time constant boundary point corresponding to the "cutoff" in the frequency domain, determined offline by the device's response characteristics, with an empirical range of values). (Within a scheduling period, those shorter than this interval are power-type devices, and those longer than this interval are energy-type devices). The curvature coefficient of the frequency domain penalty exponent (preferably 1.0). , These are the lower and upper limits of the low-frequency / high-frequency depreciation penalty, respectively (determined by the depreciation curve on the equipment nameplate).

[0058] S307, based on the frequency domain penalty, allocate high-frequency components to power devices and low-frequency components to energy devices; Specifically, based on frequency domain penalty, high-frequency components are forcibly allocated to power-type devices ( ), distributing low-frequency components to energy-type devices ( To be implemented during the daytime phase, the frequency division allocation coefficient is defined as follows: (54) satisfy , This enables frequency division and coordinated scheduling with optimal lifetime loss.

[0059] The power command mapping relationship is used to determine the baseline power output plan during the day-ahead scheduling phase and to generate power output correction commands during the intraday real-time phase.

[0060] During the daytime operation phase, when real-time monitoring reveals a deviation between the actual output of new energy sources and the predicted value... In this case, there is no need to re-invoke the large-scale nonlinear optimization model; the real-time correction logic is triggered directly. Real-time controller acquires error Using a bandpass filter , Extract the error components of each resource allocation frequency band: (55) In the formula, This represents the in-window error sequence ending at the current time t; The bandpass filter , It is implemented using a causal recursive digital filter (IIR) with phase compensation, with parameters preset to approximate zero-phase characteristics; the filtering object is the error sequence within a real-time sliding time window.

[0061] Based on the allocation coefficients of the high and low frequency components stored recently and Calculate the real-time output correction commands for each resource. : (56) Final instructions The grid connection point injection power is obtained through aggregation mapping. The power tracking is achieved by transmitting the data to the underlying distributed terminal through the power dispatch communication protocol. To prevent real-time correction commands from causing physical limits to be exceeded by the equipment, the controller introduces real-time safety monitoring, projecting the grid connection point aggregation command onto the corrected feasible domain. : (57) in, To the modified feasible region The projection. If the real-time correction command deviates from the safety boundary, the controller will force it to be projected onto the target area. The most feasible approach is to ensure that the power correction remains within the physical limits of safe operation of the distribution network; the aggregated instructions after projection are replaced back to the individual machine instructions according to the original allocation ratio of each resource.

[0062] In the above formula, each parameter is obtained from the nameplate parameters of the actual physical system, historical operating data, or offline experimental calibration, and its specific value is determined by those skilled in the art based on the actual equipment configuration and operating conditions of the virtual power plant.

[0063] 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.

[0064] 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.

[0065] 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.

[0066] 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.

[0067] 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.

[0068] 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 virtual power plant optimization scheduling method based on physical-electrical decoupling, characterized in that, Includes the following steps: Construct the feasible domain of injected power at each time period node, and establish cross-time period ramping and state transition constraints to form time-series constraints; The temporal constraints are transformed into a dual function using Lagrange relaxation. A regularization penalty term based on the physical response time constant of the adjustable resource is introduced. The dual function is then subjected to smooth approximation, and a gradient correction term is calculated. The physical response time constant is defined as the shortest time required for the adjustable resource to ramp up from zero power to full power at its maximum ramp rate, and is expressed in terms of the number of scheduling periods. The gradient correction term is updated to the feasible region boundary of the node injection power to obtain the corrected feasible region sequence; Within the corrected feasible domain sequence, a power command mapping relationship for prediction error is established, and a baseline power output plan and a power output correction command are output based on the power command mapping relationship.

2. The method according to claim 1, characterized in that, The steps for constructing the feasible region of node injection power include: Obtain the safe operation domain of the distribution network and the output feasible domain of each adjustable resource; Extract the output boundary parameters of each of the said feasible output regions, and construct a linear inequality constrained polyhedron; Perform Minkowski and aggregation operations on each of the polyhedra to generate an initial running domain; The intersection of the initial operating domain and the safe operating domain is taken, and the tangent plane algorithm is used for iterative optimization to generate an inscribed polyhedron as the feasible domain for node injection power.

3. The method according to claim 2, characterized in that, The steps for obtaining the safe operation domain of the distribution network include: Performing second-order cone relaxation on the AC power flow equations yields a set of convex inequalities; The set of convex inequalities is approximated by a polyhedron using a piecewise linearization method to obtain a linearized set of constraints, wherein the piecewise density is dynamically adjusted according to the rate of change of the load forecast curve. Using the power injection vector at the grid connection point as the retained variable, Fourier-Motzkin elimination is performed on the linearized constraint set to eliminate internal state variables and obtain the safe operation domain of the distribution network.

4. The method according to claim 1, characterized in that, The calculation steps for the gradient correction term include: Construct a regularization coefficient matrix, wherein each diagonal element of the regularization coefficient matrix is ​​determined by the physical response time constant of the adjustable resource constrained by each time-series constraint; By applying a quadratic regularization term to the dual function using the regularization coefficient matrix, a smooth dual function is generated. According to the envelope theorem, the gradient vector of the smooth dual function with respect to the Lagrange multipliers is solved, and used as the gradient correction term.

5. The method according to claim 4, characterized in that, The step of updating the gradient correction term to the feasible region boundary of the node injection power to obtain the corrected feasible region sequence includes: Establish a sliding buffer between the difference of dual variables and the difference sequence of gradient vectors, wherein the step size of the sliding window is dynamically set according to the ratio of the fastest and slowest physical response time constants of the adjustable resource; The sliding buffer is recursively calculated according to the rank-2 update rule in the quasi-Newton method to construct an approximate matrix of the Hessian inverse matrix; Perform matrix multiplication between the approximate matrix and the current gradient vector to determine the update direction of the dual variable.

6. The method according to claim 5, characterized in that, The obtained modified feasible region sequence further includes: Calculate the obstacle correction vector of the dual boundary in the temporal constraint; Calculate the margin from the current dual variable to the temporal constraint boundary, whereby the margin is defined as the ratio of the remaining settling time of the ramp constraint to the total settling time; Calculate the mapping weight coefficients based on the rate of change of the margin; The update direction and the obstacle correction vector are linearly weighted and fused according to the weight coefficients to update the dual variables.

7. The method according to claim 1, characterized in that, The establishment of the power command mapping relationship regarding the prediction error includes: Centered on the empirical distribution of historical prediction errors, a fuzzy set of the probability distribution of the true error is constructed based on the Wasserstein distance, wherein the radius of the Wasserstein distance is dynamically and adaptively set according to the time autocorrelation decay coefficient of the prediction error. Perform a strong dual transformation on the expected cost objective function containing the fuzzy set, transforming it into a set of linear matrix inequality constraints; The set of linear matrix inequality constraints is used as a deterministic penalty term and included in the scheduling objective function for equivalent solution.

8. The method according to claim 7, characterized in that, The steps for constructing the power command mapping relationship include: The output commands for each time period are set as a power command mapping relationship between the baseline output plan and the real-time adjustment allocation coefficient, and a balance constraint is applied to ensure that the sum of each resource allocation coefficient is 1. The power command mapping relationship is substituted into the modified node injection power feasible region, and the solution is performed using the baseline output plan and allocation coefficient as joint decision variables. The response capability boundary of each resource adjustment allocation coefficient is determined based on the magnitude ratio of the gradient correction term.

9. The method according to claim 8, characterized in that, The steps for solving the joint decision variables further include: Set a depreciation cost coefficient for resource response under different frequency fluctuations as a frequency domain penalty; Perform spectral analysis on the prediction error sequence and separate the high-frequency and low-frequency components; Based on the frequency domain penalty, high-frequency components are allocated to power-type devices, and low-frequency components are allocated to energy-type devices.

Images