Quantification method of flexibility of island integrated energy system with flexible domain

Abstract

The application provides a flexible domain fusion island comprehensive energy system flexibility quantification method, and relates to the technical field of comprehensive energy system operation optimization. The application firstly constructs a multi-energy flow unified matrix operation model, constructs a Zono polyhedron according to source and load prediction error data, introduces the model, and determines the flexible domain which is feasible to all disturbances; the flexible domain is off-line described by using sparse polynomial chaos expansion and minimum volume circumscribed ellipsoid, and the ellipsoid volume and balance degree are obtained; a Gaussian process regression proxy model is trained to evaluate the candidate expansion scheme; and a multi-objective plant growth algorithm is used to solve the Pareto front. The application can solve the problems of insufficient source and load coupling description, difficult flexibility boundary quantification, and large online optimization calculation, and realizes fast evaluation of geometric characteristics and screening of expansion schemes.

Description

Technical Field

[0001] This invention relates to the field of integrated energy system operation optimization technology, and in particular to a method for quantifying the flexibility of an island integrated energy system that incorporates a flexible domain. Background Technology

[0002] Currently, the development and utilization of islands is gradually emerging globally, with an increasing number of islands being developed and utilized. Islands possess abundant, directly usable clean energy sources such as solar and wind power, while clean fossil fuels like natural gas are increasingly being transported to island households via undersea pipelines or LNG. Constructing an Integrated Energy System (IIES) as a comprehensive energy supply solution for islands can contribute to the rapid economic and social development of islands and improve the quality of life for island residents.

[0003] However, unlike the interconnected grids on the mainland, island systems are connected to the main grid via submarine cables with limited capacity. This weak interconnection characteristic determines the unique nature of IIES flexibility issues. First, external support capabilities are limited. The transmission capacity constraints of submarine cables make it difficult for the system to obtain timely and sufficient external power support when encountering fluctuations in wind and solar power output. The system must mainly rely on its own resources to cope with disturbances, and the flexibility of supply directly determines the safe operating boundary of the system, making it far more important than for mainland interconnected systems. Second, the superposition effect of uncertainties is significant. The penetration rate of renewable energy on islands is high, and wind and solar power are often affected by the ocean monsoon, showing a strong positive correlation in output. At the same time, the diverse loads on islands are significantly affected by seasonal climate, with summer cooling and peak electricity consumption highly coupled. Uncertainties on both the source and load sides are superimposed, further exacerbating the difficulty of system balancing. Third, the multi-energy coupling dimension is deep. Due to the special needs of islands in terms of freshwater supply and clean fuels, IIES often integrate equipment such as electrolytic hydrogen production, methanation, and seawater source heat pumps, forming a deeply coupled operating structure. The types of dispatchable equipment are numerous, the energy flow interaction relationships are complex, and the geometric shape of the system's feasible operating space is more complex. Under the combined effect of the above characteristics, the strength of the IIES's flexible adjustment capability directly determines the safe operating boundary of the system, and therefore, accurately quantifying the system's own flexibility is of great significance.

[0004] In the field of energy systems, flexibility is generally defined as the system's ability to adapt to random changes in source-load at a certain cost by optimizing the allocation of available resources. Currently, research on quantifying the flexibility of integrated energy systems includes two main approaches. Some studies have established flexibility supply-demand matching models to analyze the operational flexibility requirements within the net load variation range and the flexibility provided by various heterogeneous energy sources. These studies assess the system's regulatory capacity by constructing flexibility supply-demand balance constraints. Another approach integrates flexibility assessment into the scheduling optimization process. Under given uncertainties, it judges the system's flexibility by whether a feasible scheduling scheme can be solved, using expected values ​​of load shedding and energy curtailment as evaluation indicators. However, these methods typically output a scalar indicator or a specific scheduling scheme, rather than an explicit quantitative tool that intuitively reflects the overall regulatory boundary of the system. Decision-makers find it difficult to understand the complete feasible operating space of the system under multiple constraints, to directly compare different schemes in terms of flexibility, and to perceive whether there are flexibility shortcomings in various energy flow directions. For IIES, the above limitations are particularly prominent. When external regulation resources are not available, operators urgently need to grasp the regulation margin boundaries of the system in each energy flow dimension, rather than a general comprehensive score.

[0005] The feasible region method, as an explicit geometric quantification tool, can characterize the set of system operating points satisfying all security constraints as a high-dimensional geometric space. Its geometric form directly and comprehensively reflects the system's flexible adjustment capability under multiple constraints, and has received widespread attention in energy system security analysis and flexibility research in recent years. In the field of power systems, a relatively mature security region theory has been developed, covering various forms such as distribution network security region, inertia security region, and virtual power plant feasible region, providing a powerful tool for system security assessment and flexible resource potential analysis. In the field of integrated energy systems, related research has extended the feasible region method to multi-energy coupling scenarios such as electricity-gas and electricity-heat, and explored issues such as dimensionality reduction observation, geometric feature analysis, and optimization of the region. However, existing research still has significant limitations in applying the feasible region method to quantify flexibility: First, most existing studies focus on characterizing and visualizing the boundaries of feasible regions, lacking research on systematically utilizing the geometric characteristics of the region as an explicit indicator of flexibility; second, existing methods generally fail to provide detailed modeling for the uncertain structure of IIES with strong wind-solar correlation and deep multi-energy coupling, resulting in overly conservative or distorted characterization of the uncertainty set; third, there is a lack of systematic work on how to efficiently construct feasible regions for integrated energy systems under high-dimensional uncertainty and extract flexibility information from the geometric shape of the region. Summary of the Invention

[0006] The purpose of this invention is to provide a method for quantifying the flexibility of an island integrated energy system that integrates the flexible domain. This method can solve the problems of insufficient characterization of source-load coupling uncertainty, difficulty in explicitly quantifying the flexibility boundary, and large computational load in online optimization of candidate expansion schemes. It enables rapid evaluation of the geometric characteristics of the flexible domain and multi-objective screening of expansion schemes.

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

[0008] A method for quantifying the flexibility of an integrated island energy system that incorporates flexible domains includes:

[0009] Based on the dispatchable equipment output plan, energy storage status, and uncertain parameters formed by wind power output, photovoltaic power output, and diverse load demands of the island integrated energy system, a unified matrix operation model for multiple energy flows is constructed.

[0010] Based on the source-load prediction error data, the Qinuo polyhedron is constructed. The Qinuo polyhedron is used to characterize the strong correlation between wind and solar power and the coupling characteristics of multiple loads. The disturbances corresponding to the Qinuo polyhedron are introduced into the unified matrix operation model of multi-energy flow.

[0011] Based on the multi-energy flow unified matrix operation model and the Chino polyhedron, a set of operation schemes that remain feasible for all disturbances is determined, and the set of operation schemes is defined as the system's flexible domain.

[0012] The flexible domain of the system is characterized offline by sparse polynomial chaotic expansion and minimum volume circumscribed ellipsoid, and the ellipsoidal volume and equilibrium degree used to characterize the flexibility of the system are obtained.

[0013] Using the expansion capacity as input and the ellipsoidal volume and balance as training objectives, a Gaussian process regression surrogate model is trained. The flexible domain geometric features of candidate expansion schemes are then evaluated online using the Gaussian process regression surrogate model.

[0014] With the optimization objectives of maximizing the volume of the ellipsoid and minimizing the cost of expanding key equipment, a multi-objective plant growth algorithm is used to solve the Pareto front, generating a set of expansion optimization schemes.

[0015] Preferably, based on the dispatchable equipment output plan, energy storage status, and uncertain parameters formed by wind power output, photovoltaic power output, and diverse load demands of the island's integrated energy system, a unified matrix-based operation model for multiple energy flows is constructed, including:

[0016] A vector of operating points is generated based on the output plan of the schedulable equipment.

[0017] State variables are formed based on the energy storage state;

[0018] Uncertainty parameters are formed based on wind power output, solar power output, and diverse load demands;

[0019] Using the equipment in the island integrated energy system as nodes and the energy flow in the island integrated energy system as branches, an energy flow conversion matrix is ​​established to characterize the conversion relationship between electricity, heat, gas and hydrogen energy flow.

[0020] Based on the running point vector, state variables, uncertainty parameters, and energy flow transformation matrix, a unified matrix-based operating model for multiple energy flows, which includes energy flow balance relationships and operating constraints, is formed.

[0021] Preferably, a Chino polyhedron is constructed based on source-load prediction error data. This Chino polyhedron characterizes the strong correlation between wind and solar power and the coupling characteristics of multiple loads. The disturbances corresponding to the Chino polyhedron are then introduced into a unified matrix-based multi-energy flow operation model, including:

[0022] The prediction center value, the fluctuation half-width of each uncertainty, and the coupling amplitude between the source and the load are determined based on the source-load prediction error data.

[0023] An independent fluctuation generator is formed based on the fluctuation half-width of each uncertainty;

[0024] Based on the correlation between wind power output and photovoltaic power output and the corresponding coupling amplitude, a wind-solar coupling fluctuation generator is formed;

[0025] Based on the coupling relationship between electrical load and thermal load and the corresponding coupling amplitude, a load coupling fluctuation generator is formed;

[0026] The independent fluctuation generator, the wind-solar coupled fluctuation generator, and the load coupled fluctuation generator are merged into a generator matrix;

[0027] Constructing the Chino polyhedron based on the predicted center value and the generator matrix;

[0028] The perturbation corresponding to the Chino polyhedron is represented as the sum of the predicted center value and the perturbation term, which is determined by the generator matrix and the normalized perturbation coefficient vector.

[0029] Based on the generator matrix and the energy flow balance relationship in the multi-energy flow unified matrix operation model, an uncertainty transfer matrix is ​​formed.

[0030] By using the uncertainty transfer matrix, the source-load fluctuations represented by the Chino polyhedron are mapped to the energy flow balance relationship of the multi-energy flow unified matrix operation model.

[0031] Preferably, based on the multi-energy flow unified matrix operation model and the Chino polyhedron, a set of operation schemes that remain feasible for all disturbances is determined, and this set of operation schemes is defined as the system's flexible domain, including:

[0032] For any perturbation realization of the Chino polyhedron coverage, determine the set of operation schemes that satisfy the unified matrix operation model of multi-energy flow;

[0033] Take the intersection of the set of operational schemes corresponding to all disturbances covered by the Chino polyhedron;

[0034] The operational schemes in the intersection are determined as those that remain feasible for all disturbances.

[0035] The intersection is defined as the system's flexible domain.

[0036] Preferably, the flexible domain of the system is characterized offline using sparse polynomial chaotic expansion and a minimum volume circumscribed ellipsoid to obtain the ellipsoidal volume and equilibrium degree used to characterize the system's flexibility, including:

[0037] Determine feasible operating points within the flexible domain of the system and use these feasible operating points as search base points;

[0038] Set the search direction in the run point space, and form candidate run points along the search direction;

[0039] Construct a constraint violation function to characterize whether candidate running points satisfy the running feasibility under the perturbation corresponding to the Chino polyhedron;

[0040] Using the normalized running point and the perturbation corresponding to the Chino polyhedron as joint inputs, an approximate model of the constraint violation function is established by sparse polynomial chaotic expansion.

[0041] The critical distance between candidate running points and the search direction is determined based on the approximate model, and a set of flexible domain boundary sample points is formed based on the critical distance.

[0042] Construct a minimum volume circumscribed ellipsoid based on the sample point set of the flexible domain boundary, and determine the geometric features of the flexible domain based on the ellipsoid volume and balance of the minimum volume circumscribed ellipsoid.

[0043] Preferably, using the normalized running point and the perturbation corresponding to the Chino polyhedron as joint inputs, an approximate model of the constraint violation function is established using sparse polynomial chaotic expansion, including:

[0044] The normalized running point and the perturbation corresponding to the Chino polyhedron are combined into a joint input;

[0045] Polynomial chaotic expansion of constraint violation functions is performed using Legendre orthogonal polynomial basis;

[0046] The polynomial chaotic expansion is truncated to form a truncated model containing a finite number of expansion terms;

[0047] Construct an orthogonal polynomial information matrix based on the truth values ​​of the constraint violation functions corresponding to the training samples;

[0048] Sparse regression is used to screen expansion terms that contribute to the constraint violation function from the truncated model;

[0049] An approximate model of the constraint violation function is formed based on the expanded terms after filtering.

[0050] Preferably, a minimum volume circumscribed ellipsoid is constructed based on the set of sample points at the boundary of the flexible domain, and the geometric features of the flexible domain are determined based on the ellipsoidal volume and uniformity of the minimum volume circumscribed ellipsoid, including:

[0051] The minimum volume circumscribed ellipsoid is characterized by a center vector and a shape matrix;

[0052] Make the minimum volume circumscribed ellipsoid contain the boundary sample points in the set of boundary sample points of the flexible domain;

[0053] Given the boundary sample points, determine the smallest volume circumscribed ellipsoid.

[0054] Eigenvalue analysis is performed on the shape matrix to determine the principal axis directions of the circumscribed ellipsoid of the minimum volume and the corresponding semi-axis lengths for each principal axis direction;

[0055] The volume of the ellipsoid is determined based on the semi-axis lengths corresponding to each principal axis direction.

[0056] The degree of balance is determined based on the difference in the half-axis length corresponding to each main axis direction;

[0057] The volume and uniformity of the ellipsoid are defined as geometric features of the flexible domain.

[0058] Preferably, a Gaussian process regression surrogate model is trained using the expansion capacity as input and the ellipsoidal volume and balance as training objectives. The Gaussian process regression surrogate model is then used to perform online evaluation of the flexible domain geometric characteristics of candidate expansion schemes, including:

[0059] Generate expanded capacity samples within the feasible range of capacity expansion;

[0060] For each expanded capacity sample, sparse polynomial chaotic expansion and minimum volume circumscribed ellipsoid are used to offline characterize the corresponding system flexible domain, and the corresponding ellipsoid volume and balance are obtained.

[0061] Using expanded samples as training input and the corresponding ellipsoidal volume and equilibrium as training objectives, Gaussian process regression surrogate models for predicting ellipsoidal volume and equilibrium are established respectively.

[0062] For candidate expansion schemes, the predicted values ​​of ellipsoidal volume, equilibrium degree, and prediction variance are calculated using a Gaussian process regression surrogate model.

[0063] The flexible domain geometric characteristics of candidate expansion schemes are determined based on the predicted values ​​of ellipsoidal volume and equilibrium.

[0064] The online evaluation process is managed using a proxy model based on the predicted variance.

[0065] Preferably, the online evaluation process is managed by a proxy model based on the predicted variance, including:

[0066] Based on the predicted variance distribution of the initial training samples, determine the ellipsoid volume prediction variance threshold and the equilibrium prediction variance threshold.

[0067] When the ellipsoid volume prediction variance of a candidate expansion scheme exceeds the ellipsoid volume prediction variance threshold, or the equilibrium prediction variance of a candidate expansion scheme exceeds the equilibrium prediction variance threshold, the system flexible domain offline characterization is re-executed for the corresponding candidate expansion scheme to obtain the updated ellipsoid volume and equilibrium.

[0068] Add the corresponding candidate expansion schemes, as well as the updated ellipsoid volume and balance, to the training dataset.

[0069] Update the Gaussian process regression surrogate model based on the updated training dataset;

[0070] For candidate expansion schemes where the equilibrium prediction value satisfies the equilibrium constraint and the equilibrium prediction variance exceeds the equilibrium prediction variance threshold, the corresponding candidate expansion schemes are included in the proxy model update object.

[0071] Preferably, with the optimization objectives of maximizing the ellipsoidal volume and minimizing the expansion cost of key equipment, a multi-objective plant growth algorithm is used to solve the Pareto front, generating a set of expansion optimization schemes, including:

[0072] Generate a population of candidate expansion schemes based on equipment expansion capacity constraints;

[0073] The flexible domain geometric characteristics of each candidate expansion scheme in the candidate expansion scheme population were calculated using a Gaussian process regression surrogate model.

[0074] With the optimization objectives of maximizing the volume of the ellipsoid and minimizing the cost of expanding key equipment, the candidate expansion schemes are non-dominated and ranked.

[0075] Determine elite candidate expansion schemes based on non-dominated ranking results;

[0076] Based on the analytical gradient of the Gaussian process regression surrogate model, the negative gradient of the expansion cost of key equipment, and the random unit vector, the growth direction of elite candidate expansion schemes is determined.

[0077] The elite candidate expansion schemes are updated according to the growth direction, and the updated elite candidate expansion schemes are projected onto the feasible domain determined by the equipment expansion capacity constraints.

[0078] The Pareto front is formed through iteration, and a set of expansion optimization schemes is generated based on the Pareto front.

[0079] The present invention discloses the following beneficial effects:

[0080] This invention constructs a unified matrix-based operation model for multi-energy flow, incorporating the dispatchable equipment output plan, energy storage status, and uncertainties arising from wind power output, photovoltaic power output, and diverse load demands into a single operational expression for the island's integrated energy system. This addresses the problem of inconsistent operational feasibility foundations caused by fragmented modeling of multi-energy flow conversion, coupling, and external exchange relationships, ensuring that subsequent flexibility quantification is based on system operational constraints. Furthermore, by constructing a Chino polyhedron using source-load prediction error data and introducing the corresponding disturbances into the unified matrix-based operation model, this invention resolves the issue of distortion of uncertainty sets caused by independently processing the strong correlation between wind and solar power and the coupling characteristics of diverse loads, enabling the system's flexibility domain to reflect the feasible operational range under the superposition of source-load uncertainties. Finally, by defining the set of operational schemes that remain feasible for all disturbances as the system's flexibility domain, this invention addresses the problem that existing flexibility evaluations only form a single indicator or a single dispatch scheme. The problem of difficulty in representing a complete feasible operating space is addressed by limiting the flexible adjustment capability of the island integrated energy system to a set of operating schemes. The flexible domain of the system is characterized offline using sparse polynomial chaotic expansion and a minimum volume circumscribed ellipsoid, with ellipsoid volume and equilibrium degree representing system flexibility. This solves the problems of large computational cost at the boundary of high-dimensional flexible domains and difficulty in representing differences in adjustment margins across different energy flow directions, enabling flexibility evaluation to simultaneously reflect the size of the feasible operating space and the degree of equilibrium in the distribution of adjustment capability. Furthermore, the geometric characteristics of the flexible domain of candidate expansion schemes are evaluated online using a Gaussian process regression surrogate model. With maximizing ellipsoid volume and minimizing the expansion cost of key equipment as optimization objectives, a multi-objective plant growth algorithm is used to solve the Pareto front, addressing the problem of large computational cost in online optimization caused by individually characterizing the flexible domain of each candidate expansion scheme. This allows expansion optimization schemes to achieve a comparable trade-off between increased flexibility and expansion cost. Attached Figure Description

[0081] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0082] Figure 1 This is a flowchart of a method provided in an embodiment of the present invention.

[0083] Figure 2 The 24-time flexible domain ellipsoid time series distribution diagram provided in the embodiment of the present invention.

[0084] Figure 3 Scatter plot of sparse polynomial chaotic expansion prediction accuracy provided in an embodiment of the present invention.

[0085] Figure 4 Pareto front diagram provided for embodiments of the present invention.

[0086] Figure 5 A comparison diagram of the principal axes of the ellipsoid before and after optimization, provided for embodiments of the present invention. Detailed Implementation

[0087] 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 skilled in the art without creative effort are within the scope of protection of the present invention.

[0088] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0089] Figure 1 The method flowchart provided in the embodiments of the present invention is as follows: Figure 1 As shown, this embodiment provides a method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains, including:

[0090] Step 100: Based on the dispatchable equipment output plan, energy storage status, and uncertain parameters formed by wind power output, photovoltaic power output, and diverse load demands of the island integrated energy system, construct a unified matrix operation model for multiple energy flows;

[0091] Step 200: Construct the Chino polyhedron based on the source-load prediction error data, characterize the strong correlation between wind and solar power and the coupling characteristics of multiple loads through the Chino polyhedron, and introduce the disturbance corresponding to the Chino polyhedron into the unified matrix operation model of multi-energy flow.

[0092] Step 300: Based on the multi-energy flow unified matrix operation model and the Chino polyhedron, determine the set of operation schemes that remain feasible for all disturbances, and define the set of operation schemes as the system flexible domain;

[0093] Step 400: The flexible domain of the system is characterized offline by using sparse polynomial chaotic expansion and minimum volume circumscribed ellipsoid to obtain the ellipsoidal volume and equilibrium degree used to characterize the flexibility of the system.

[0094] Step 500: Using the expansion capacity as input and the ellipsoidal volume and balance as training objectives, train a Gaussian process regression surrogate model, and use the Gaussian process regression surrogate model to evaluate the flexible domain geometric features of candidate expansion schemes online;

[0095] Step 600: With the optimization objectives of maximizing the volume of the ellipsoid and minimizing the expansion cost of key equipment, the Pareto front is solved using a multi-objective plant growth algorithm to generate a set of expansion optimization schemes.

[0096] Specifically, in step 100 of this embodiment, the island integrated energy system achieves renewable energy consumption and system regulation capacity enhancement through multi-energy flow coupling and coordination. The island integrated energy system consists of a supply side, a conversion side, an energy storage side, and a load side. The supply side includes wind power generation equipment and photovoltaic power generation equipment; the conversion side includes combined heat and power units, proton exchange membrane electrolyzers, methanation units, hydrogen fuel cells, and seawater source heat pumps; the energy storage side includes batteries, thermal storage tanks, and hydrogen storage tanks; and the load side includes electrical load, thermal load, and gas load.

[0097] This embodiment defines the running point vector. Indicates time Output plans for all schedulable devices, defining state variables. Indicates time Energy storage state, defining uncertainty parameters This represents a vector composed of externally uncontrollable input parameters in the model, which takes the predicted value or a value for a specific scenario in a deterministic model. (Running point vector) State variables and uncertainty parameters Represented as:

[0098]

[0099]

[0100]

[0101] In the formula, , , For example, cogeneration units at Power generation, heat generation, and gas consumption at any given time; , For proton exchange membrane electrolyzers, respectively Power consumption and hydrogen production at any given time; , , For hydrogen fuel cells Hydrogen consumption power, electricity generation power, and heat generation power at any given moment; , The methanation unit is located at Hydrogen consumption power and gas production power at any given time; , Seawater source heat pumps are respectively Power consumption and heat generation at any given time; for External gas supply power at all times; , , The batteries are respectively The charging power, discharging power, and energy storage capacity at any given time; , , The heat storage tanks are respectively The charging power, dissipation power, and heat storage capacity at any given time; , , The hydrogen storage tanks are respectively The hydrogen charging power, hydrogen discharging power, and hydrogen storage capacity at any given time; , The outputs are wind power and solar power, respectively. , , These are respectively the electricity load demand, heat load demand, and gas load demand.

[0102] This embodiment defines a joint state vector. Used to fully characterize the island's integrated energy system at any time. Operating status:

[0103]

[0104] In the formula, For the running point vector Dimensions State variables Dimensions Joint state vector Dimensions.

[0105] Based on the concept of energy hub graph theory, this embodiment treats each device as a node and each energy flow as a branch. For nodes... The connection relationship between ports and all branches is represented by the port-branch association matrix. Description. Branch road Access Node When the input port is selected, the corresponding element is taken. branch road Outgoing Node When the output port is specified, the corresponding element is taken. In other cases, take , is represented as:

[0106]

[0107] The node energy conversion efficiency is determined by the conversion characteristic matrix. Description. Input port association conversion process. When, the corresponding element is the conversion efficiency. Output port association conversion process When, the corresponding element is In other cases, take , is represented as:

[0108]

[0109] In the formula, For the conversion process Energy conversion efficiency. Node Energy flow conversion matrix Port-Branch Association Matrix and transformation characteristic matrix get:

[0110]

[0111] Integrating the entire island's comprehensive energy system From the nodes, we obtain the system energy flow transformation matrix. :

[0112]

[0113] System energy flow conversion matrix This is used to describe the energy conversion relationships among the various devices within the system. The input correlation matrix is ​​further defined. and output correlation matrix Input the correlation matrix Associating external inputs with the system, the output is an association matrix. The load demand of the associated system. The system coupling and conservation relationships can be uniformly represented as:

[0114]

[0115] In the formula, The external input vector consists of wind power output, photovoltaic power output, and external gas supply. The load demand vector consists of electrical load demand, heat load demand, and gas load demand.

[0116] This embodiment integrates all operational constraints into a unified matrix form. The energy flow balance equation constitutes equality constraints, while the equipment capacity constraints, ramp rate constraints, and energy storage state evolution constraints constitute inequality constraints, expressed as follows:

[0117]

[0118]

[0119] In the formula, Energy flow conversion matrix Input correlation matrix and output correlation matrix The combination forms a system used to characterize the conversion, coupling, and external exchange relationships of electrical, thermal, gaseous, and hydrogen energy flows. Due to uncertainty parameters The predicted values ​​constitute; This is the inequality constraint coefficient matrix; and These are the constraint boundary vectors.

[0120] Optionally, in step 200 of this embodiment, the geographical scope of the island integrated energy system is limited, and both wind power output and photovoltaic power output are dominated by the marine climate, with output deviations exhibiting structural correlation characteristics; electrical load and heat load are affected by seasonal climate and have a synergistic coupling relationship. If the source-load correlation is ignored, the uncertainty set is prone to overconservatism. Therefore, this embodiment uses a Chino polyhedron to model the source-load uncertainty.

[0121] From the central vector and generator vector Defined Chino polyhedron Represented as:

[0122]

[0123] In the formula, Uncertainty parameter Predicted value ; For the first One normalized perturbation coefficient; For the first Generator vectors for each independent perturbation mode. The generator vectors characterize the fluctuation range of the corresponding variables and the correlation between variables through the position and value of their non-zero elements. (Cinno polyhedron) It can be simply remembered as:

[0124]

[0125] In the formula, This is the generator matrix.

[0126] The generator matrix consists of two parts: independent wave generators and coupled wave generators. The independent wave generators characterize the self-waves of each uncertainty. Let the first... The half-width of the fluctuation of an uncertain quantity is Independent fluctuation generator is represented as only in the first... A column vector with non-zero components:

[0127]

[0128] In the formula, For the first Standard basis vectors.

[0129] A coupled wave generator characterizes the correlation between source and load. Since wind and solar power outputs are governed by the same ocean climate system, a wind-solar coupled wave generator is introduced. The wind-solar coupled fluctuation generator takes values ​​from both the wind power output component and the photovoltaic power output component; the electrical load and thermal load are affected by seasonal climate, so a load coupled fluctuation generator is introduced. The load coupling fluctuation generator takes values ​​simultaneously on both the electrical load component and the thermal load component:

[0130]

[0131]

[0132] In the formula, The amplitude of wind power output coupling. For the photovoltaic output coupling amplitude, For electrical load coupling amplitude, This represents the thermal load coupling amplitude. The fluctuation half-width is... Both the coupling amplitude and the source-load prediction error data are determined by fitting the data.

[0133] By merging the independent wave generator and the coupled wave generator, we obtain the generator matrix. :

[0134]

[0135] The uncertainty parameter is decomposed into the predicted central value and the disturbance term:

[0136]

[0137] Substituting the above equation into the energy flow balance equation, the energy flow balance equation under uncertainty becomes:

[0138]

[0139] In the formula, The uncertainty transfer matrix is ​​used to transfer the generator matrix. The source-load fluctuations are mapped to the input and load components on the right-hand side of the energy flow equation. Through the above processing, the disturbances corresponding to the Chino polyhedron are introduced into the multi-energy flow unified matrix operation model, which can characterize the operational feasibility under source-load uncertainties.

[0140] Furthermore, in step 300 of this embodiment, for a given disturbance within the uncertain set... Define the feasible set of the system The set of all operating schemes that satisfy the system constraints:

[0141]

[0142] In the formula, and They are respectively , Chinese correspondence The coefficient matrix; The corrected right-hand side vector after removing the influence of state variables:

[0143]

[0144] In the formula, for Chinese correspondence The coefficient matrix.

[0145] Define time flexible domain of the system The intersection of the feasible sets under all disturbances covered by the uncertain set, i.e., the set of feasible operating schemes that are maintained for any possible source-load fluctuations:

[0146]

[0147] Equivalently, the flexible domain of the system is represented as:

[0148]

[0149] System Flexible Domain The shape and size are affected by the running point The influence of different operating points results in different flexible domain geometries. The system's flexible domain geometric volume can be used as a measure of the system's performance at different times. The overall flexible adjustment capability is a quantitative indicator. The difference in the extension of the system's flexible domain in different directions can characterize the difference in the distribution of flexibility in the directions of electric, heat, gas, and hydrogen energy flow.

[0150] Furthermore, in step 400 of this embodiment, the system flexible domain As a high-dimensional implicit set, the boundary is implicitly determined by the operational feasibility conditions. In scheduling optimization scenarios, it is necessary to evaluate the geometric features of the flexible domain in real time for a large number of candidate solutions. If the entire flexible domain characterization process is re-executed each time, the computational cost cannot meet the requirements for online use. This embodiment adopts an offline-online collaborative flexible domain characterization and optimization framework. In the offline stage, the system's flexible domain is characterized, geometric feature indicators are extracted, and a lightweight surrogate model is constructed. In the online stage, the surrogate model is called to drive multi-objective optimization, thereby reducing the computational time of geometric feature evaluation and optimization.

[0151] This embodiment employs a boundary sampling method based on direction search and binary search to discretize and approximate the boundary of the system's flexible domain. First, a deterministic optimization problem is solved to obtain a known feasible operating point within the system's flexible domain as the search base point. Then preset in the run point space. Search direction vectors For each search direction From the search base point Starting from the beginning, a one-dimensional search is performed along the ray, and the critical distance is determined through a binary search. The boundary points are obtained:

[0152]

[0153] In the formula, For the first Boundary points in each search direction, For the first Critical distance in each search direction.

[0154] The key to binary search lies in judging any candidate points on the path. Does it belong to the system's flexible domain? Define the constraint violation function. For disturbance The negative of the maximum constraint margin that the system can maintain is solved by the following linear programming problem:

[0155]

[0156]

[0157]

[0158] In the formula, Candidate running points The corresponding feasible adjustment amount, To constrain the margin, It is a vector of all 1s with the same dimension as the constraint boundary vector. At that time, the maximum margin is satisfied The system is under disturbance There exists a feasible adjustment that satisfies all constraints; when When this happens, the constraints will inevitably be violated.

[0159] For fixed candidate running points , The judgment can be transformed into taking the worst-case scenario among all possible perturbations. Define the worst-case constraint violation quantity. :

[0160]

[0161] Candidate running point Belongs to the flexible domain of the system The equivalent condition is:

[0162]

[0163] Thus, the binary search forms a two-layer nested structure: the outer layer along the search direction. Adjust candidate running points The inner layer has a fixed candidate running point. Solve for worst-case constraint violation. That is, for uncertain sets The process involves enumerating all vertices, solving for constraint violations, and finding the maximum value. Each step forward in the outer layer triggers a large number of linear programming problems to be solved in the inner layer. The total computational cost is high after combining multiple search directions and binary iteration.

[0164] To reduce computational cost, this embodiment employs sparse polynomial chaotic expansion techniques to construct constraint violation functions. Approximate model This replaces the repeated calls to linear programming in binary search. Worst-case constraint violation rate. The evaluation is transformed from iteratively solving linear programming problems into solving problems with uncertain sets. Vertex-based approximation model Evaluate by enumeration.

[0165] The decision variables to be solved in the constraint violation function are: running point and disturbance All of these appear in the right-hand side of the constraints. According to parametric linear programming theory, the optimal value of a linear programming problem is a continuous piecewise linear function with respect to the right-hand side parameters. Therefore... yes A continuous piecewise linear function; on a bounded domain, a continuous function is square-integrable. For uniformly distributed input variables... Any square-integrable function can be approximated by a linear combination of Legendre orthogonal polynomial bases.

[0166] Will After normalization and As a combined input:

[0167]

[0168] In the formula, the combined input Each component , , , To run the point vector dimension, For the perturbation dimension. Regarding joint input Perform a polynomial chaotic expansion:

[0169]

[0170]

[0171] In the formula, For joint input component index, For joint input The One portion, For the first The corresponding joint input components A single-variable Legendre polynomial of order 1. For multi-level subscripts, For the first 1 variable univariate Legendre polynomials; multivariate basis functions Composed of tensor products of all univariate basis functions, satisfying orthogonality and normalization:

[0172]

[0173] The expansion coefficients are to be determined.

[0174] Since the chaotic expansion of a polynomial is an infinite series, it cannot be directly calculated. This embodiment considers that the coefficients of higher-order terms in the expansion decay rapidly with the order, and calculates the series based on the total order. Truncate, retaining the set of indices :

[0175]

[0176] Corresponding number of expanded terms for:

[0177]

[0178] In one implementation, the total order Pick to This ensures calculation accuracy; in this embodiment... Truncation Model Represented as:

[0179]

[0180] After truncation One unknown coefficient. Extraction sample points For each sample point, apply the linear programming problem corresponding to the constraint violation function to obtain the true value. Construct sample output vector and orthogonal polynomial information matrix :

[0181]

[0182]

[0183] The expansion coefficients are obtained by minimizing the prediction residuals:

[0184]

[0185] Number of expanded terms With dimension sum order The growth is rapid, and a direct solution requires a large number of samples. Since the coefficients of most higher-order cross terms are close to zero, the expansion coefficients... It exhibits sparsity. This embodiment transforms the calculation of expansion coefficients into a sparse regression problem:

[0186]

[0187] In the formula, expansion factor The number of non-zero elements in neutron This represents the upper limit of sparsity. This embodiment uses the least angle regression algorithm to solve the problem, automatically identifying pairs from a finite training sample. A set of sparse terms that makes significant contributions, which avoids an explosion in the number of terms when expanding the model while ensuring model accuracy.

[0188] Obtain an approximate model Then, the worst-case constraint violation quantity The assessment is transformed into The polynomial is evaluated by enumerating it on the vertex set, ultimately yielding the set of boundary points. The set of boundary points constitutes the flexible domain for complex shape systems. A discrete sampling of the boundary.

[0189] The set of boundary points obtained based on sparse polynomial chaotic expansion and binary search can characterize the flexible domain of the original high-dimensional system. The system has complex geometries, but its non-convex, multi-cellular structure is not conducive to fast lookup and volume calculation. Therefore, this embodiment uses a minimum-volume circumscribed ellipsoid to perform a smooth, conservative approximation of the system's flexible domain. Unlike accurately describing complex polyhedra, which requires numerous inequality constraints, the ellipsoid is determined by a central vector. and shape matrix The only certainty:

[0190]

[0191] In the formula, The center of the ellipsoid represents the point where the system operates robustly. This is a shape matrix that describes the shape and orientation of the ellipsoid.

[0192] To ensure the conservatism of the ellipsoid approximation, the ellipsoid... Including all boundary points and minimizing the redundant space introduced by approximations, this embodiment constructs the following convex optimization problem:

[0193]

[0194]

[0195]

[0196] In the formula, As slack variables, For relaxation budget. Relaxed variables. and slack budget Used to enhance the model's robustness to sampling errors.

[0197] To reduce the difficulty of calculating the volume of a high-dimensional ellipsoid and to give the volume calculation a clear physical meaning, this embodiment introduces an adaptive dimensionality reduction method based on eigenvalue analysis on the basis of ellipsoid approximation. For the shape matrix... Perform eigenvalue decomposition:

[0198]

[0199]

[0200]

[0201] In the formula, An orthogonal matrix with column vectors Represents the ellipsoid of the first One main axis direction; For the first Eigenvalues ​​corresponding to each principal axis direction. The size reflects the system in the first The adjustment capability along each principal axis direction, where the eigenvalues ​​are arranged in descending order, i.e. Minimal eigenvalues ​​correspond to dimensions that contribute little to system flexibility. This can be addressed by setting a relative threshold. Preserving the principal dimension with significant eigenvalues ​​enables flexible automatic dimensionality reduction:

[0202]

[0203] In the formula, For a valid set of main dimensions, The relative threshold is used. The ellipsoid after dimensionality reduction is... Composed of several main dimensions, this approach reduces the complexity of volume calculation and scheduling optimization while preserving the main directions of flexibility.

[0204] The system's flexibility potential at a given moment is expressed through the system's flexibility domain. The volume is quantized. For the circumscribed ellipsoid of the minimum volume after dimensionality reduction, the semi-axis length is expressed as:

[0205]

[0206] The volume of an ellipsoid is expressed as:

[0207]

[0208] In the formula, For the volume of the ellipsoid, For the effective main dimension set The number of main dimensions included. This is the Gamma function. The ellipsoidal volume is used to quantify the overall flexible adjustment capability of the system.

[0209] The flexible domain volume reflects the total flexibility of the system, but it cannot reveal the distribution structure of flexibility in each energy flow direction. If the lengths of the principal axes differ significantly, even with a large total volume, the system may still have structural shortcomings in certain energy flow dimensions due to insufficient regulation capacity. Therefore, this embodiment defines an equilibrium index. The coefficient of variation of the lengths of each principal axis is used to measure the degree of equilibrium of the ellipsoid shape:

[0210]

[0211]

[0212]

[0213] In the formula, For the flexible domain ellipsoid in the th The length of the semi-axis in the direction of each principal axis This represents the average length of the semi-axis. This represents the standard deviation of the semi-shaft length. The closer The more consistent the lengths of each main axis are, the more balanced the system's flexible domain's adjustment capability is in the key energy flow direction; The larger the ellipsoid, the more skewed its shape, indicating a relatively insufficient adjustment margin in some dimensions.

[0214] This embodiment can be used for typical days. The system's flexibility domain is characterized at each moment, and the flexibility at each moment is explicitly quantified using ellipsoidal volume and isostatics. In one implementation, a typical day of the transitional season is used as the object, and the flexibility is analyzed over the entire day. The flexible domain ellipsoid is continuously characterized at each time step, and adaptive dimensionality reduction is performed through eigenvalue analysis, while retaining the principal axis directions that contribute significantly. , Visualize it.

[0215] Figure 2 This is a time-series distribution diagram of the flexible domain ellipsoid at 24 time points provided in an embodiment of the present invention. Figure 2 As shown, the horizontal time axis represents different times within a typical day, and the spatial coordinate axis... and These represent the two main flexibility directions obtained through eigenvalue analysis, and the ellipsoids at different times represent the minimum volume circumscribed ellipsoid approximation of the flexible domain of the system at the corresponding time. Figure 2 The color bar on the right indicates The color changes from a low value range to a high value range, corresponding to an increase in the volume of the ellipsoid. Figure 2 The red markers indicate the weakest points. The time corresponding to Blue marks indicate the most flexible moments. The time corresponding to . Figure 2 The data also includes several typical time points, reflecting the spatial location, shape, and volume variations of the flexible domain ellipsoid throughout the day. From early morning to dawn, the ellipsoid is relatively large, corresponding to a relatively stable wind power output, low load, and sufficient energy storage capacity. During the daytime, as electrical and thermal loads rise simultaneously, the ellipsoid gradually shrinks. During the nighttime peak transition period, wind and solar power output weaken while the load has not yet fully receded, increasing the uncertainty on both the source and load sides, causing the system's flexible domain to shrink to a lower level. Figure 2 Able to determine The nearby time is a time with low flexibility, and the corresponding time can be used as the object of expansion and optimization analysis.

[0216] Figure 3 This is a scatter plot showing the prediction accuracy of sparse polynomial chaotic expansion provided in an embodiment of the present invention. Figure 3 As shown, the horizontal axis represents the actual value of constraint violation. The vertical axis represents the predicted values ​​from the sparse polynomial chaotic expansion. ,dotted line As a reference line for predicting values ​​to equal true values, horizontal and vertical zero-value reference lines divide the scatter points into feasibility judgment areas. Figure 3 The blue scatter dots represent feasible samples, while the red scatter dots represent infeasible samples. The feasible and infeasible samples show a clear separation in the predicted values, indicating that the approximate model can be used to determine the feasibility of candidate running points. Figure 3 The coefficient of determination is marked in the bottom right corner. This indicates that the prediction results of sparse polynomial chaotic expansion have a high degree of consistency with the actual results of linear programming vertex enumeration, and can replace a large number of linear programming solutions in boundary search.

[0217] Table 1. Comparison of flexible domain characterization methods provided in the embodiments of the present invention.

[0218] Boundary sampling accuracy 100% 99.7% runtime 11792s 82s Coverage rate 85.5% 91.8% Number of parameters 200 54

[0219] As shown in Table 1, the boundary sampling accuracy of the flexible domain characterization method based on sparse polynomial chaotic expansion and minimum volume circumscribed ellipsoid reaches [percentage missing]. Approximate linear programming vertex enumeration and polyhedral approximation methods Runtime is determined by Reduce to The calculation time is approximately [time] times that of the original method. Accelerate approximately Times; Coverage rate from Upgraded to More realistic and feasible scheduling scenarios have been included in the flexible domain; the number of parameters has been increased from... Reduce to This demonstrates that the ellipsoidal approximation can characterize the geometric features of the flexible domain of a system with fewer parameters.

[0220] Furthermore, in step 500 of this embodiment, the online scheduling phase requires real-time evaluation of flexibility indicators for a large number of candidate expansion schemes. If the entire process of characterizing the flexibility domain is re-executed each time, the computational cost would be insufficient to meet real-time requirements. Therefore, this embodiment establishes a Gaussian process regression surrogate model from expansion capacity to ellipsoidal volume and equilibrium in the offline phase.

[0221] Sample generation within the feasible range of capacity expansion. Representative sample of capacity expansion:

[0222]

[0223] In the formula, For the first Each expanded sample uniformly covers the feasible decision space. A complete flexible domain characterization process is performed on each expanded sample to obtain the corresponding ellipsoidal volume and equilibrium.

[0224]

[0225] The training dataset consists of expanded samples and their corresponding ellipsoidal volumes and balance values. Considering the long computation time of a single calculation in the complete characterization process, the total computational cost in the offline stage is [missing information]. This independent computation can be completed through parallel computing, without needing to meet online real-time requirements. (Sample size) A trade-off is made between the accuracy of the surrogate model and the cost of offline computation.

[0226] This embodiment uses capacity expansion. Using ellipsoidal volume and uniformity as inputs, a surrogate model is trained using Gaussian process regression. Taking the logarithm of the ellipsoidal volume to model the model yields:

[0227]

[0228]

[0229] In the formula, Let be the mean function of the logarithmic prediction model for ellipsoidal volume. The kernel function for the logarithmic prediction model of ellipsoidal volume is... The mean function of the equilibrium prediction model. This is the kernel function for the equilibrium prediction model.

[0230] The kernel function can be a quadratic exponential kernel function:

[0231]

[0232] In the formula, For signal variance, The length scale parameter is used. The signal variance and length scale parameter are determined by maximizing the logarithmic marginal likelihood of the training data. Modeling the volume using the logarithm is because... Smoother changes in the operating space help improve the prediction accuracy of Gaussian process regression surrogate models.

[0233] The Gaussian process regression surrogate model provides both the predicted mean and the predicted variance. Candidate expansion schemes. The corresponding ellipsoid volume prediction variance and equilibrium degree prediction variance are expressed as follows:

[0234]

[0235]

[0236] Prediction variance is used to evaluate the prediction confidence of the surrogate model at the current scaling configuration during the online phase. As the optimization search progresses, candidate scaling solutions may enter sparse regions of the training samples, increasing prediction variance. Directly using the prediction mean in non-dominated ranking may lead to low-confidence optimal solutions appearing in the Pareto front. Therefore, this embodiment introduces a surrogate model management strategy based on Gaussian process regression prediction variance.

[0237] The two prediction variance thresholds are adaptively determined based on the variance distribution of the initial training samples:

[0238]

[0239]

[0240] In the formula, This represents the initial number of training samples; For the first One initial training sample; , This is a proportional coefficient used to control the update frequency.

[0241] For assessing the feasibility of the balance constraint, this embodiment adopts a conservative strategy. If the predicted balance value of a candidate expansion scheme meets the balance constraint, but the predicted variance of the balance exceeds the predicted variance threshold, then the expansion scheme will proceed as planned. If the corresponding candidate expansion scheme is considered to have uncertain feasibility due to constraints, it will be included in the subsequent proxy model update set first, rather than being directly accepted as a feasible solution.

[0242] Finally, in step 600 of this embodiment, two types of problems are identified by observing the geometric characteristics of the flexible domain of the system: one is insufficient overall volume of the flexible domain, indicating that the safe operating points available under source-load uncertainty are limited, and the overall flexible adjustment capability is restricted; the other is excessive equilibrium, indicating that the extension of the flexible domain in different energy flow directions varies significantly, and some energy flow dimensions have flexibility shortcomings. Based on this, this embodiment constructs a multi-objective optimization model based on the geometric characteristics of the flexible domain to balance the improvement of flexibility and economic costs.

[0243] The objective function for the volume of the ellipsoid is:

[0244]

[0245] The objective function for the cost of expanding key equipment is:

[0246]

[0247] In the formula, Let the objective function be the volume of the ellipsoid. The objective function for expanding the capacity of key equipment; For the first Cost of expanding the capacity of each device; For the first The capacity expansion of each device relative to the existing configuration; This represents the total number of expandable devices.

[0248] The optimization model also includes constraints on equipment expansion capacity and balance:

[0249]

[0250]

[0251] In the formula, For the first Maximum allowable expansion capacity for each device; The balance improvement rate represents the minimum improvement in balance relative to the initial state after optimization. This is the initial expansion configuration.

[0252] The aforementioned optimization model has two objectives: maximizing the ellipsoidal volume and minimizing the cost of expanding key equipment. The objective function is evaluated online using a Gaussian process regression surrogate model obtained through offline training. This embodiment introduces a prediction variance management strategy based on the multi-objective plant growth algorithm to improve the reliability of the Pareto front results.

[0253] The algorithm randomly generates the initial population within the flexible domain of the system ellipsoid:

[0254]

[0255] In the formula, For population size, each individual This represents a candidate expansion scheme that satisfies the equipment expansion capacity constraint. The Gaussian process regression surrogate model is invoked to calculate the corresponding ellipsoidal volume prediction mean. ellipsoid volume prediction variance Balance Prediction Mean and the variance of equilibrium prediction Complete the initial objective function and constraint evaluation.

[0256] In each generation In the iteration, feasible individuals are first selected, and the current population is sorted using fast non-dominated sorting. Divided into One non-dominated level:

[0257]

[0258] In the formula, To determine the optimal Pareto front, calculate the crowding distance for each candidate expansion scheme to maintain the Pareto front distribution:

[0259]

[0260] In the formula, The number of objective functions; For the first The solution is at the th solution. Function values ​​on each target; and The first The maximum and minimum values ​​of the objective function in the current undominated hierarchy.

[0261] Subsequently, the algorithm enters the growth phase. It selects the current generation of elite individuals based on non-dominated hierarchy and crowding distance. For each elite solution... Calculate the multi-objective guided growth direction The growth direction is determined by the analytical gradient of the Gaussian process regression surrogate model, the negative gradient of the cost of expanding key equipment, and the random unit vector.

[0262]

[0263] In the formula, The analytical gradient of the ellipsoidal volume target with respect to the expansion capacity can be analytically calculated by a Gaussian process regression surrogate model; The negative gradient direction for the cost target of expanding the capacity of key equipment; It is a random unit vector used to maintain population diversity; , , These are the weighting coefficients.

[0264] The weighting coefficients are dynamically adjusted as the iteration progresses:

[0265]

[0266]

[0267]

[0268]

[0269] In the formula, For the maximum number of iterations, This is for iterative algebra. The above weight settings cause the early search to favor expanding the ellipsoidal volume, while the later search gradually enhances cost control, while retaining random perturbations to maintain population diversity.

[0270] growth step Adaptive mechanism adopted:

[0271]

[0272] In the formula, The initial step size, The attenuation coefficient is used. The newly generated solution is mapped back to the feasible region using the projection operator:

[0273]

[0274] In the formula, These are candidate expansion schemes generated according to the growth direction and growth step size. These are the candidate expansion schemes after projection; The feasible region determined by the equipment capacity expansion constraint. Because... Given the constraints of the linear frame, the projection problem can be solved using a numerical method based on Lagrange multipliers.

[0275] After each iteration, for candidate expansion schemes whose ellipsoid volume prediction variance or equilibrium prediction variance in the current Pareto front exceeds the corresponding threshold, the true objective function and constraint evaluation are triggered. That is, the system flexible domain characterization process in step 400 is re-executed to obtain accurate ellipsoid volume and equilibrium:

[0276]

[0277] In the formula, Update the collection of objects for the proxy model. This represents the optimal non-dominated level corresponding to the current Pareto front. The precise calculation results corresponding to the candidate expansion schemes are added to the training dataset, and the Gaussian process regression surrogate model is updated. By iteratively executing the above process, the algorithm eventually converges to a well-distributed Pareto approximate front, providing a set of schemes that quantify the trade-off between the expansion cost of critical equipment and the ellipsoidal volume for expansion optimization.

[0278] Figure 4 The Pareto front diagram provided for embodiments of the present invention. For example... Figure 4 As shown, the horizontal axis The vertical axis represents the cost of expanding key equipment. The values ​​represent the logarithmic volume of the ellipsoid, the scatter plots represent candidate expansion schemes, and the color bar on the right represents the uniformity. The color range is approximately to . Figure 4 The rhombus marker represents the base point, and the logarithmic value of the ellipsoidal volume corresponding to the base point is approximately Star-shaped marking scheme ,plan Located on the Pareto front, the logarithmic value of the corresponding ellipsoidal volume is approximately . Figure 4 The plan was also marked in the text. ,plan Typical candidate expansion solutions are presented to demonstrate different trade-offs, such as low-cost moderate flexibility and high-cost high flexibility. As the cost of expanding critical equipment increases, The overall monotonic increase indicates that the candidate expansion schemes have expanded the system's flexible domain; the balance... The variation within the range corresponding to the color bar indicates that the candidate expansion scheme satisfies the equilibrium constraint while increasing the volume of the ellipsoid. Figure 4 The Pareto front distribution shown is relatively continuous, covering a trade-off range from low cost and low flexibility to high cost and high flexibility.

[0279] Figure 5 A comparison diagram of the principal axes of the ellipsoid before and after optimization, provided for embodiments of the present invention. Figure 5 As shown in the radar chart to These represent the seven principal axes obtained from eigenvalue analysis of the circumscribed ellipsoid of the minimum volume. The radar chart's radius values ​​represent the semi-axis lengths corresponding to each principal axis. The blue broken line represents the baseline scheme, and the isostatics of the baseline scheme are shown. The green line represents the optimized solution, and the corresponding equilibrium degree of the optimized solution. .from Figure 5 It can be seen that the semi-axis length distribution of the benchmark scheme is uneven in different principal axis directions. , The relatively shorter lengths in the corresponding directions indicate a weaker adjustment margin in the energy flow scheduling degrees of freedom; the optimization scheme expands outwards overall in multiple principal axis directions. , The orientation was reinforced, and the flexibility boundary of the system in multiple energy flow dimensions was broadened. The equilibrium of the optimized scheme is close to that of the baseline scheme, indicating that while the expansion optimization increases the ellipsoidal volume, it does not cause a significant deterioration in the distribution of the principal axis adjustment capability.

[0280] This embodiment establishes the foundation for system operational feasibility through a unified matrix-based multi-energy flow operation model; characterizes the strong correlation between wind and solar power and the coupling characteristics of multiple loads through the Chino polyhedron; defines a set of operational schemes that remain feasible for all disturbances through a flexible domain constraint; reduces the number of repeated linear programming solutions in boundary search through sparse polynomial chaotic expansion; transforms the flexible domain of the high-dimensional system into ellipsoidal volume and equilibrium through a minimum volume circumscribed ellipsoid; evaluates the online flexible domain geometric characteristics of candidate expansion schemes through a Gaussian process regression surrogate model; and forms a Pareto front through a multi-objective plant growth algorithm, obtaining a set of expansion optimization schemes.

[0281] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0282] Specific examples have been used to illustrate the principles and implementation methods of this invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of this invention. Furthermore, those skilled in the art will recognize that, based on the ideas of this invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this invention.

Claims

1. A method for quantifying the flexibility of an integrated island energy system that incorporates flexible domains, characterized in that, include: Based on the dispatchable equipment output plan, energy storage status, and uncertain parameters formed by wind power output, photovoltaic power output, and diverse load demands of the island integrated energy system, a unified matrix operation model for multiple energy flows is constructed. Based on the source-load prediction error data, a Qino polyhedron is constructed. The Qino polyhedron is used to characterize the strong correlation between wind and solar power and the coupling characteristics of multiple loads. The disturbances corresponding to the Qino polyhedron are introduced into the multi-energy flow unified matrix operation model. Based on the multi-energy flow unified matrix operation model and the Qino polyhedron, a set of operation schemes that remain feasible for all disturbances is determined, and the set of operation schemes is defined as the system flexible domain. The flexible domain of the system is characterized offline by using sparse polynomial chaotic expansion and minimum volume circumscribed ellipsoid to obtain the ellipsoidal volume and equilibrium degree used to characterize the flexibility of the system. Using the expansion capacity as input, and the ellipsoidal volume and the balance as training objectives, a Gaussian process regression surrogate model is trained, and the flexible domain geometric features of candidate expansion schemes are evaluated online through the Gaussian process regression surrogate model. With the optimization objectives of maximizing the volume of the ellipsoid and minimizing the cost of expanding key equipment, a multi-objective plant growth algorithm is used to solve the Pareto front, generating a set of expansion optimization schemes.

2. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 1, characterized in that, Based on the dispatchable equipment output plan, energy storage status, and uncertain parameters formed by wind power output, photovoltaic power output, and diverse load demands of the island integrated energy system, a unified matrix-based operation model for multiple energy flows is constructed, including: A running point vector is formed based on the schedulable equipment output plan; State variables are formed based on the energy storage state; Uncertainty parameters are formed based on the wind power output, the photovoltaic power output, and the diverse load demands. Using the equipment in the island integrated energy system as nodes and the energy flow in the island integrated energy system as branches, an energy flow conversion matrix is ​​established to characterize the conversion relationship between electricity, heat, gas and hydrogen energy flow. Based on the running point vector, the state variables, the uncertainty parameters, and the energy flow conversion matrix, a unified matrix-based multi-energy flow operating model is formed, which includes energy flow balance relationships and operating constraint relationships.

3. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 1, characterized in that, A Chino polyhedron is constructed based on source-load prediction error data. This polyhedron characterizes the strong correlation between wind and solar power and the coupling characteristics of multiple loads. The disturbances corresponding to the Chino polyhedron are then introduced into the unified matrix-based multi-energy flow operation model, including: The prediction center value, the fluctuation half-width of each uncertainty, and the coupling amplitude between the source and the load are determined based on the source-load prediction error data. An independent fluctuation generator is formed based on the fluctuation half-width of each uncertainty; Based on the correlation between wind power output and photovoltaic power output and the corresponding coupling amplitude, a wind-solar coupling fluctuation generator is formed; Based on the coupling relationship between electrical load and thermal load and the corresponding coupling amplitude, a load coupling fluctuation generator is formed; The independent fluctuation generator, the wind-solar coupled fluctuation generator, and the load coupled fluctuation generator are merged into a generator matrix; The Kino polyhedron is constructed based on the predicted center value and the generator matrix; The perturbation corresponding to the Kino polyhedron is represented as the sum of the predicted center value and the perturbation term, wherein the perturbation term is determined by the generator matrix and the normalized perturbation coefficient vector; An uncertainty transfer matrix is ​​formed based on the generator matrix and the energy flow balance relationship in the multi-energy flow unified matrix operation model. The source-load fluctuations represented by the Chino polyhedron are mapped to the energy flow balance relationship of the unified matrix-based multi-energy flow operation model through the uncertainty transfer matrix.

4. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 1, characterized in that, Based on the aforementioned multi-energy flow unified matrix operation model and the Chino polyhedron, a set of operation schemes that remain feasible for all perturbations is determined, and this set of operation schemes is defined as the system's flexible domain, including: For any disturbance realization covered by the Chino polyhedron, determine a set of operating schemes that satisfy the unified matrix operation model of the multi-energy flow; Take the intersection of the set of corresponding operation schemes for all disturbances covered by the Chino polyhedron; The operating schemes in the intersection are determined as those that remain feasible for all disturbances. The intersection is defined as the flexible domain of the system.

5. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 1, characterized in that, The flexible domain of the system is characterized offline using sparse polynomial chaotic expansion and a minimum volume circumscribed ellipsoid, yielding the ellipsoidal volume and equilibrium degree used to characterize the system's flexibility, including: Determine feasible operating points within the flexible domain of the system, and use these feasible operating points as search base points; A search direction is set in the run point space, and candidate run points are formed along the search direction; Construct a constraint violation function to characterize whether the candidate running point satisfies the running feasibility under the perturbation corresponding to the Kino polyhedron; Using the normalized running point and the perturbation action corresponding to the Kino polyhedron as joint inputs, an approximate model of the constraint violation function is established by sparse polynomial chaotic expansion. Based on the approximate model, the critical distance between the candidate running points and the search direction is determined, and a set of flexible domain boundary sample points is formed according to the critical distance. Construct a minimum volume circumscribed ellipsoid based on the set of sample points at the boundary of the flexible domain, and determine the geometric features of the flexible domain based on the ellipsoid volume and balance of the minimum volume circumscribed ellipsoid.

6. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 5, characterized in that, Using the normalized running point and the perturbation corresponding to the Chino polyhedron as joint inputs, an approximate model of the constraint violation function is established using sparse polynomial chaotic expansion, including: The normalized running point and the perturbation corresponding to the Kino polyhedron are combined into a joint input; The constraint violation function is expanded using a Legendre orthogonal polynomial basis; The polynomial chaotic expansion is truncated to form a truncated model containing a finite number of expansion terms; Construct an orthogonal polynomial information matrix based on the truth values ​​of the constraint violation functions corresponding to the training samples; Sparse regression is used to select expansion terms that contribute to the constraint violation function from the truncated model; An approximate model of the constraint violation function is formed based on the filtered expanded terms.

7. The method for quantifying the flexibility of an island integrated energy system with integrated flexible domains according to claim 5, characterized in that, A minimum volume circumscribed ellipsoid is constructed based on the set of boundary sample points of the flexible domain, and the geometric features of the flexible domain are determined based on the ellipsoidal volume and balance of the minimum volume circumscribed ellipsoid, including: The minimum volume circumscribed ellipsoid is characterized by a center vector and a shape matrix; Make the minimum volume circumscribed ellipsoid contain the boundary sample points in the set of boundary sample points of the flexible domain; Given the boundary sample points, determine the smallest circumscribed ellipsoid with the smallest volume; Eigenvalue analysis is performed on the shape matrix to determine the principal axis directions of the circumscribed ellipsoid of the minimum volume and the semi-axis lengths corresponding to each principal axis direction; The volume of the ellipsoid is determined based on the semi-axis lengths corresponding to each principal axis direction; The degree of balance is determined based on the difference in the semi-axis lengths corresponding to each main axis direction. The ellipsoidal volume and the uniformity are defined as the geometric features of the flexible domain.

8. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains according to claim 1, characterized in that, Using the expansion capacity as input, and the ellipsoidal volume and the equilibrium degree as training objectives, a Gaussian process regression surrogate model is trained. The Gaussian process regression surrogate model is then used to perform online evaluation of the flexible domain geometric features of candidate expansion schemes, including: Generate expanded capacity samples within the feasible range of capacity expansion; For each expanded capacity sample, the sparse polynomial chaotic expansion and the minimum volume circumscribed ellipsoid are used to offline characterize the corresponding system flexible domain to obtain the corresponding ellipsoid volume and balance. Using the expanded sample as training input and the corresponding ellipsoidal volume and balance as training objectives, Gaussian process regression surrogate models for predicting ellipsoidal volume and Gaussian process regression surrogate models for predicting balance are established respectively. For each candidate expansion scheme, the predicted ellipsoidal volume, the predicted equilibrium value, and the predicted variance are calculated using the Gaussian process regression surrogate model. The flexible domain geometric characteristics of the candidate expansion scheme are determined based on the predicted ellipsoidal volume and the predicted equilibrium degree. The online evaluation process is managed by a proxy model based on the predicted variance.

9. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains as described in claim 8, characterized in that, The online evaluation process is managed by a proxy model based on the predicted variance, including: Based on the predicted variance distribution of the initial training samples, determine the ellipsoid volume prediction variance threshold and the equilibrium prediction variance threshold. When the ellipsoid volume prediction variance of the candidate expansion scheme exceeds the ellipsoid volume prediction variance threshold, or the balance prediction variance of the candidate expansion scheme exceeds the balance prediction variance threshold, the system flexible domain offline characterization is re-executed for the corresponding candidate expansion scheme to obtain the updated ellipsoid volume and balance. Add the corresponding candidate expansion schemes, as well as the updated ellipsoid volume and balance, to the training dataset. The Gaussian process regression surrogate model is updated based on the updated training dataset; For candidate expansion schemes whose equilibrium prediction values ​​satisfy the equilibrium constraints and whose equilibrium prediction variance exceeds the equilibrium prediction variance threshold, the corresponding candidate expansion schemes are included in the proxy model update object.

10. The method for quantifying the flexibility of an island integrated energy system by incorporating flexible domains according to claim 1, characterized in that, With the optimization objectives of maximizing the ellipsoidal volume and minimizing the expansion cost of key equipment, a multi-objective plant growth algorithm is used to solve the Pareto front, generating a set of expansion optimization schemes, including: Generate a population of candidate expansion schemes based on equipment expansion capacity constraints; The flexible domain geometric characteristics of each candidate expansion scheme in the candidate expansion scheme population are calculated using the Gaussian process regression surrogate model. With the optimization objectives of maximizing the volume of the ellipsoid and minimizing the cost of expanding key equipment, the candidate expansion schemes are non-dominated and ranked. Determine elite candidate expansion schemes based on non-dominated ranking results; Based on the analytical gradient of the Gaussian process regression surrogate model, the negative gradient of the expansion cost of key equipment, and the random unit vector, the growth direction of the elite candidate expansion scheme is determined. The elite candidate expansion scheme is updated according to the growth direction, and the updated elite candidate expansion scheme is projected onto the feasible domain determined by the device expansion capacity constraint. The Pareto front is formed through iteration, and the set of scaling optimization schemes is generated based on the Pareto front.

Images