Off-grid wind and light storage system capacity configuration stochastic optimization method

Abstract

The invention relates to a stochastic optimization method for capacity configuration of an off-grid wind and light storage system. The method comprises the following steps: constructing a system mathematical model by considering equipment constraint and power balance; considering the uncertainty of wind speed and illumination intensity and performing probability distribution fitting by using a parameter estimation method; constructing a stochastic optimization algorithm by using an expected value model and a Monte Carlo method; a non-dominated sorting genetic algorithm is used for carrying out multi-objective optimization calculation on economic benefits and carbon emission; a new method for evaluating the permeability of renewable energy sources is provided by combining a Monte Carlo method definition, and investment decision suggestions are given by taking the new method as an index. The problem that the renewable energy source permeability of the off-grid energy system with uncertainty lacks specific quantification is solved. According to the method, the uncertainty of renewable energy sources is considered, compared with a deterministic method, the method has higher robustness, and the system uncertainty caused by volatility is effectively solved.

Description

Technical Field

[0001] The present invention belongs to the field of integrated energy, and specifically relates to a random optimization method for capacity configuration of an off-grid wind-solar-storage system taking into account the uncertainty of wind and solar. Background Art

[0002] In off-grid integrated energy systems including wind, solar and storage, capacity planning methods usually involve source-load input, system modeling and optimization algorithms. The precise quantitative characterization and prediction of "source-load" uncertainty is the basis for the optimization of integrated energy systems. Renewable energy sources such as wind and solar have significant intermittent and random fluctuations, which have a great impact on the stable operation of distributed energy systems. The uncertainty of renewable energy increases the difficulty of system planning, design and operation scheduling, especially the capacity configuration of energy storage equipment and the ultra-short-term energy regulation strategy of the system. In the uncertainty method, after analyzing the sources and characteristics of the uncertainty in the planning model, mathematical tools such as probability theory and sequence operations are used to model it, and the impact of uncertainty on the planning results is studied through sensitivity analysis, stochastic programming, robust programming and other methods. Energy planning is a multi-objective optimization problem. The Pareto frontier can be directly solved by the multi-objective optimization algorithm, or the multi-objective can be integrated into a single-objective optimization problem by normalization method. Economic feasibility and environmental sustainability are usually considered.

[0003] In off-grid systems, capacity planning and energy balance are affected by the volatility of renewable energy generation forecasts. Traditional deterministic optimization methods have input limitations, resulting in redundant equipment capacity and high costs. Multiple uncertainties are a difficult problem in planning. Among them, the uncertainty in renewable energy output is difficult to quantify, which increases the difficulty of modeling. Mainstream uncertainty methods usually consider random scenario generation for uncertainty optimization, lacking direct mathematical representation of volatility. At the same time, current technology still lacks comprehensive research on the renewable energy penetration rate of off-grid energy systems with uncertainty, and lacks specific quantification of the renewable energy penetration rate of the system. Summary of the invention

[0004] Aiming at the lack of specific quantification of renewable energy penetration rate in off-grid energy systems with uncertainty, a stochastic optimization method for capacity configuration of off-grid wind-solar-storage systems was proposed.

[0005] The technical solution of the present invention is:

[0006] A random optimization method for off-grid wind-solar-storage system capacity configuration comprises the following steps:

[0007] The mathematical model of the system is constructed by considering equipment constraints and power balance; the uncertainty of wind speed and light intensity is considered and the probability distribution is fitted using parameter estimation methods; a random optimization algorithm is constructed using the expected value model and the Monte Carlo method; a multi-objective optimization calculation of economic benefits and carbon emissions is performed using a non-dominated sorting genetic algorithm; a new method for evaluating the penetration rate of renewable energy is proposed in combination with the definition of the Monte Carlo method, and investment decision recommendations and reliability evaluation are given using it as an indicator;

[0008] Step 1: Construct a system mathematical model by considering equipment constraints and power balance:

[0009] Step 2: Establish a step-by-step model for carbon trading costs. A step-by-step model is used to calculate carbon trading costs. The model defines different emission ranges. If carbon emissions are lower than the specified quota, the cost is negative, indicating that energy suppliers can sell excess carbon emission rights in the carbon trading market and obtain technical subsidies. The transaction price increases as the emission range decreases. On the contrary, when emissions exceed the quota, the cost becomes positive, indicating that energy suppliers need to purchase additional carbon emission rights.

[0010] Step 3: Establish the objective function of the integrated energy system planning model:

[0011] The integrated energy system planning model aims to minimize the total investment cost and carbon trading cost; the total cost includes investment cost and operating cost;

[0012] Step 4: Establish constraints for the integrated energy system planning model:

[0013] The constraints mainly include system operation constraints, system installation constraints and power balance constraints:

[0014] Step 5: Determine the optimization variables of the integrated energy system planning model and the mathematical expression of the planning problem:

[0015] Step 6: Consider the uncertainty of wind speed and light intensity and use parameter estimation method to perform probability distribution fitting to study the renewable energy penetration rate REPR of the system;

[0016] How to deal with uncertainty:

[0017] Random distribution prediction can provide more uncertainty information for the operation of integrated energy and power systems; use parameter estimation methods; it presupposes the distribution form, performs statistical analysis on historical data, estimates the specific distribution form through parameter estimation based on historical data, and finally optimizes the objective function of the expected value model by generating multiple scenarios; in order to study the renewable energy penetration rate REPR of the system by quantifying uncertainty and characterizing economic feasibility;

[0018] Step 7: Use multiple new energy output scenarios and their corresponding probabilities to describe the uncertainty of output; establish a random expectation scheduling model, whose goal is to minimize the expected value of the multi-objective function;

[0019] Step 8: Monte Carlo simulation to solve the expected value model

[0020] Only the source-side scenario consisting of wind turbines and photovoltaics was considered, and Monte Carlo simulation was used to calculate to determine the optimal solution;

[0021] Step 9: Find the optimal algorithm for the planning problem

[0022] The proposed model considers both carbon emissions and economic costs, so a multi-objective solution algorithm is required; in order to ensure the practicality and accuracy of the model solution, NSGA-II is used to solve the model;

[0023] Step 10: Calculate confidence interval using bootstrap method

[0024] Choose an appropriate sample size to balance computational effort and precision; since the function is multivariate and complex, a nonparametric bootstrap method is used to calculate the confidence interval;

[0025] Step 11: Investment decision analysis and reliability assessment methods

[0026] The Monte Carlo sampling method is used to evaluate the system cost under different sample sizes. The Pareto optimal solution set obtained by the NSGA-II algorithm is used to obtain the carbon emissions and costs of different solutions. Reliability assessment is performed corresponding to different renewable energy penetration rates REPR and the use of diesel generators, and the impact of REPR on system costs is analyzed. The optimal investment plan under different situations can be analyzed numerically. The economic feasibility of the system is evaluated based on the payback period, net present value, internal rate of return, and levelized electricity cost analysis of different REPRs to find the economic optimum.

[0027] Further, the method specifically includes the following steps:

[0028] Step 1: Construct a system mathematical model by considering equipment constraints and power balance:

[0029] Step 1.1: Establish the mathematical model of electric boiler EB, a commonly used heating equipment. The mathematical model and operating constraints are as follows:

[0030]

[0031] Where: P EB (t) represents the power output of the electric boiler at time t (kW); H EB (t) represents the heating capacity of the electric boiler at time t (kW); d EBRepresents the capacity of the electric boiler (kW); Indicates the maximum output power limit of the electric boiler (kW); Indicates the minimum output power limit of the electric boiler (kW); Inv EB represents the investment cost of electric boiler; TC EB represents the total cost of the electric boiler; p ele Indicates the electricity price (kWh); H EB_limit Indicates the output power adjustment limit of the electric boiler (kW / h); η EB Indicates the thermal energy conversion efficiency of electric boilers;

[0032] Step 1.2: Establish the mathematical model of the electric refrigerator EC, a commonly used refrigeration equipment. The mathematical model and operating constraints are as follows:

[0033]

[0034] Where: P EC (t) represents the power output of the electric refrigerator at time t (kW); C EC (t) represents the cooling capacity of the electric refrigerator at time t (kW); d EC Represents the capacity of the electric refrigerator (kW); Indicates the maximum output power limit of the electric refrigerator (kW); Indicates the minimum output power limit of the electric refrigerator (kW); Inv EC represents the investment cost of the electric refrigerator; TC EC represents the total cost of the electric refrigerator; p ele Indicates the electricity price (kWh); C EC_limit Indicates the output power adjustment limit of the electric refrigerator (kW / h); η EC It represents the energy conversion efficiency of electric refrigerator;

[0035] Step 1.3: Establish the mathematical model of diesel generator DG. Diesel generator DG is a reliable source of electricity used to power various applications; it generates electricity by using diesel fuel. The mathematical model and operating constraints are described as follows:

[0036]

[0037] fc(t) represents the fuel consumption of the diesel generator at time t (L); P DGN Indicates the rated power of the diesel generator (kW); P DG (t) represents the output power of the diesel generator at time t (kW); A DG represents the fuel intercept (L / kW); B DG represents the fuel slope coefficient (L / kW); p fuelIndicates fuel price ( / L); Inv DG Represents the investment cost of diesel generator; TC DG represents the total cost of the diesel generator; Indicates the maximum output power limit of the diesel generator; Indicates the minimum output power limit of the diesel generator; CE DG represents the carbon emission of diesel generator (kg); η DG represents the carbon emission rate of diesel generators;

[0038] Step 1.4: Mathematical model of the electric energy storage system ES. The mathematical model of the electric energy storage system ES is represented based on the state of charge SOC, as described below:

[0039]

[0040] Where: P in (t) represents charging power (kW); P out (t) represents the discharge power (kW); E stored (t) represents the state of charge of the stored energy at time t (kWh); d ES Indicates the capacity of the energy storage system (kW); E max Indicates the maximum capacity of electrical energy storage (kWh); E min Indicates the minimum capacity of electrical energy storage (kWh); η charge Indicates charging efficiency; η discharge Indicates discharge efficiency; P in,max and P out,max Represents the maximum limit of charging and discharging power (kW); P in,min and P out,min Represents the minimum limit of charging and discharging power (kW); E initial Indicates the initial state of charge of the energy storage (kWh); Inv ES Represents the investment cost of the energy storage system; TC ES represents the total cost of the energy storage system;

[0041] Step 2: Establish a stepwise model for carbon trading costs. A stepwise model is used to calculate carbon trading costs. The model defines different emission ranges. If carbon emissions are lower than the specified quota, the cost is negative, which indicates that energy suppliers can sell excess carbon emission rights in the carbon trading market and obtain technical subsidies. The transaction price increases as the emission range decreases. On the contrary, when emissions exceed the quota, the cost becomes positive, indicating that energy suppliers need to purchase additional carbon emission rights. At this time, a larger emission range corresponds to a higher transaction price. The model divides actual carbon emissions into linear segments for analysis:

[0042]

[0043] CE p =CE DG

[0044] In this formula: represents the carbon trading cost in the integrated energy system; μ represents the reward coefficient; ctp represents the carbon trading price in the market; θ represents the growth rate of the carbon trading price at each step; CE p Represents the total carbon emissions of the system equipment; CE c represents the system's free carbon emission quota; ν represents the length of the carbon emission interval; CE DG represents the carbon emissions of diesel generators;

[0045] Step 3: Establish the objective function of the integrated energy system planning model:

[0046] The integrated energy system planning model aims to minimize the total investment cost and carbon trading cost; the total cost includes investment cost and operating cost; the multi-objective function is expressed as follows:

[0047]

[0048] Where N is the number of years of operation; r is the interest rate; F 1 represents the objective function of total investment cost; F 2 represents the objective function of carbon trading cost; Min means the optimization target is the minimum value;

[0049] Step 4: Establish constraints for the integrated energy system planning model:

[0050] The constraints mainly include system operation constraints, system installation constraints and power balance constraints:

[0051] System operation constraints: defined by the mathematical models of various devices in step 1;

[0052] Installation constraint: Ensure that the capacity of the running device does not exceed the capacity limit of the available devices, expressed as:

[0053]

[0054] Where: M is the set of all device types, M = [EB, EC, ES, DG]; i represents the device type; d i Indicates the capacity of the corresponding device; Indicates the maximum capacity limit of device type i;

[0055] Power balance constraint: includes the balance of electricity, heat energy and gas energy, expressed as follows:

[0056]

[0057] Where: P WD (t) represents the output power of the fan at time t; P PV (t) represents the output power of the photovoltaic at time t; H EB (t) represents the output heat energy of the electric boiler (EB) at time t; C EC (t) represents the output cooling power of the electric refrigerator (EC) at time t; L ELOAD (t) represents the power load at time t; L HLOAD (t) represents the heat load at time t; L CLOAD (t) represents the cooling load at time t; ξ represents the adjustment factor, which is usually equal to 0;

[0058] Step 5: Determine the optimization variables of the integrated energy system planning model and the mathematical expression of the planning problem:

[0059] The investment decision variable is denoted as x, and the operational optimization variable is denoted as y; they are expressed as follows:

[0060]

[0061] Where: d WT Indicates the capacity of the fan; d PV represents the capacity of photovoltaic; in summary, the general form of the planning problem is as follows:

[0062]

[0063] Among them: F(x) represents the objective function set; Min means that the optimization target is the minimum value; G i (x) represents the inequality constraint; J i (x) represents the equality constraint; m represents the number of inequality constraints; n represents the total number of constraints; i represents the constraint number;

[0064] Step 6: Consider the uncertainty of wind speed and light intensity and use parameter estimation method to perform probability distribution fitting to study the renewable energy penetration rate REPR of the system;

[0065] How to deal with uncertainty:

[0066] Random distribution prediction can provide more uncertainty information for the operation of integrated energy and power systems; use parameter estimation methods; it presupposes the distribution form, performs statistical analysis on historical data, estimates the specific distribution form through parameter estimation based on historical data, and finally optimizes the objective function of the expected value model by generating multiple scenarios; in order to study the renewable energy penetration rate REPR of the system by quantifying uncertainty and characterizing economic feasibility; the uncertainty processing method model is processed in the following steps:

[0067] Step 6.1: Use probability density function to represent the uncertain renewable energy generation; then, based on historical statistical data, propose parameter estimation to obtain the corresponding distribution characteristics;

[0068] Step 6.2: Construct a random expectation model to represent the uncertainty of the output in the form of a scenario set; use Monte Carlo simulation to generate an accurate output scenario set;

[0069] Step 6.3: Consider the constraints and optimization variables to obtain the planning problem of the off-grid integrated energy system; finally, use the multi-objective algorithm to solve the problem;

[0070] Historical statistical data are usually used to fit the distribution characteristics of uncertain renewable energy generation; this process aims to derive a probability density function or a cumulative distribution function; the probability models commonly used to simulate wind speed and solar irradiation include Weibull distribution and beta distribution; therefore, Weibull distribution and beta distribution are used to simulate wind speed and solar irradiation; the standardized solar irradiation G p =G / G max Following the beta distribution model, its probability density function f(G p |α,β) is expressed as follows:

[0071]

[0072] Where G represents the actual irradiance, the unit is kilowatt per square meter kW / m 2 ; G max is the maximum irradiance, also in kilowatts per square meter kW / m 2 α and β are the two shape parameters of the Beta distribution; B represents the Beta function; these parameters are affected by the mean value of the irradiance, μ p and standard deviation Influence of factors such as:

[0073]

[0074] The photovoltaic system output power is given by the following formula:

[0075]

[0076] Where: P pv Indicates the actual output power of the photovoltaic panel, in kilowatts (kW); P st It is the maximum output power under standard test conditions, in kilowatts (kW); G st represents standard irradiation conditions, in kilowatts per square meter (kW / m^2); η pv Represents the efficiency of the photovoltaic panel; η pvs Indicates the efficiency of the photovoltaic system;

[0077] The wind speed change v follows the Weibull distribution, and its probability density function f(v|c,λ) is expressed as follows:

[0078]

[0079] where c is the shape parameter and λ is the scale parameter; in the case of wind speed distribution, c determines the shape of the distribution: it indicates that higher or lower wind speeds are less common, with moderate wind speeds dominating in most cases; as c increases, the tail of the distribution becomes heavier, indicating that higher wind speeds are more likely and lower wind speeds are less likely; the scale parameter λ controls the scale or amplitude of the wind speeds; specifically, larger values ​​of λ indicate a wider range of wind speed amplitudes, which may have higher peaks and a wider range of fluctuations; conversely, smaller values ​​of λ indicate a narrower range of wind speed amplitudes, indicating that wind speeds do not vary as drastically; these parameters are affected by, for example, the average value of the irradiance, μ w and standard deviation σ w Influence of factors such as:

[0080]

[0081] Where: Γ represents the Gamma function; the fan system output power is given by the following formula:

[0082]

[0083] Where: P wt represents the actual output power of the wind turbine components in kilowatts (kW); v represents the wind speed in meters per second (m / s); v in is the starting wind speed, i.e. the minimum wind speed at which the wind turbine starts generating electricity, expressed in meters per second (m / s); v out is the cut-out wind speed, i.e. the maximum wind speed at which the wind turbine stops generating electricity to prevent mechanical damage, in meters per second (m / s); P n is the rated power under standard conditions, in kilowatts (kW); η wt Indicates the efficiency of the wind turbine;

[0084] Likelihood estimation is an important method for estimating the distribution parameters of a data set for subsequent analysis and modeling. The extracted irradiance and wind speed data will be analyzed using the likelihood estimation method. After preprocessing, these data are divided into multiple data points for each hour of the whole year, represented as PD i , where i represents the index of the illumination data point; in order to estimate the distribution parameters of the data, the negative log-likelihood function is chosen; the core goal of this function is to evaluate the degree of fit between the data and a specific probability distribution; specifically, the mathematical definition of the negative log-likelihood function is as follows:

[0085]

[0086] Where: α l ,β l ∈Ω;Ω={(α l ,β l )|α l >0,β l >0}; In the above equation, f(PD i |α l ,β l ) indicates that the parameter α l ,β l Data point PD under the condition i The probability density function of the data point is obtained by calculating the negative natural logarithm of the probability density function of each data point and summing it over all data points to obtain the negative log-likelihood value NL. In order to find the parameter value that best explains the data distribution, the fitting error between the data and the selected probability distribution must be minimized. By minimizing the negative log-likelihood function NL, we seek to find the distribution parameter α that best fits the data. l ,β l ;

[0087] Step 7: Use multiple new energy output scenarios and their corresponding probabilities to describe the uncertainty of output; establish a random expectation scheduling model, whose goal is to minimize the expected value of the multi-objective function:

[0088] F=E[F 1 ,F 2 ]

[0089] Where F represents the expected function of the random programming objective; E represents the mathematical expectation value;

[0090] The random expected value model constructed by this method characterizes the uncertainty of new energy output in the form of a scenario set; generating an accurate output scenario set is the key and prerequisite for efficiently solving the model; Monte Carlo simulation is used to generate the output scenario set; due to the uncertainty of wind speed and solar irradiation, the optimal cost of the expected value model is a function f(v,G p ), where G p Follows beta distribution, v follows weibull distribution; when f(v,G p ) contains two different random variables v and G p , respectively follow the probability density function f v (v) and When the Monte Carlo method is used to estimate the expected value of the function E(f(v,G p )); Generate N′ groups of independent samples where v i is from f v (v) the sample drawn, is from The sample is drawn, where i represents the sample number; the sample mean estimate can be calculated

[0091]

[0092] Next, calculate the expected value of this estimate Since v,G p is independent and has:

[0093]

[0094] Since each sample are drawn independently, and each expected value can be calculated separately:

[0095]

[0096] Sample mean estimation The expected value of is equal to the true expected value; in order to calculate the variance of the sample mean estimate

[0097]

[0098] Due to x i and i are independent, and taking into account the properties of the variance, the following relationship holds:

[0099]

[0100] Therefore, when the sample size N' tends to infinity, tends to zero; this means that as the sample size increases, the variance of the sample mean estimate decreases, resulting in a more stable estimate that is closer to the true value; therefore, when using the Monte Carlo method to estimate the sample mean, choosing a sufficiently large N′ is the key to ensuring calculation accuracy;

[0101] Step 8: Monte Carlo simulation to solve the expected value model

[0102] Only the source-side scenario consisting of wind turbines and photovoltaics is considered, and Monte Carlo simulation is used for calculation to determine the optimal solution; by using this method, the relationship between renewable energy penetration rate REPR and cost is revealed, thus providing suggestions for investment decisions; in addition, by quantifying REPR and uncertainty, the economic feasibility of the system is characterized; compared with directly calculating a certain proportion of renewable energy supply to the total load, the probabilistic description of renewable energy penetration rate is used; based on the original load data, it is assumed that renewable energy meets the load with a certain probability; in the sampling calculation process, it is convenient to screen out days without wind or sunshine and high load, and diesel power generation is used as a backup; the probabilistic form of renewable energy penetration rate is defined as:

[0103] REPR=Pr(P DG (t)=0,ξ≤0)

[0104] Where Pr represents probability; REPR can be understood as an event (P DG (t)=0,ξ≤0) occurs, and X is defined to follow a 0-1 distribution. If the event occurs, then X=1, and we can get: REPR=Pr(X=1); According to the expected value calculation formula of the 0-1 distribution, we can get: REPR=E(X); To calculate REPR, first, considering a fairly large sample size (N′), apply the law of large numbers, and for all

[0105]

[0106] Where: X 1 ,…,X N′ is a 0-1 random variable that is independent and identically distributed in X, ∈ is any positive real number;

[0107] Therefore, the sample mean Approximately estimate E(X), in the process of Monte Carlo simulation scenario generation, the scenario that successfully meets the load will be close to N′·REPR; in N′·(1-REPR) calculations, the constraints that do not meet the power balance condition are allowed to realize the probabilistic model; therefore, when using the Monte Carlo method for simulation, the choice of sample size greatly affects the accuracy of the probability calculation of the source output meeting the load;

[0108] Step 9: Find the optimal algorithm for the planning problem

[0109] The proposed model takes into account both carbon emissions and economic costs, so a multi-objective solution algorithm is required. To ensure the practicality and accuracy of the model solution, NSGA-II is used to solve the model. NSGA-II operates based on the non-dominated sorting principle and uses an elite strategy to divide the population into different non-dominated levels. Compared with other algorithms, NSGA-II has powerful search capabilities and high parallelization characteristics, making it particularly suitable for dealing with nonlinear multi-objective problems.

[0110] Step 10: Calculate confidence interval using bootstrap method

[0111] Therefore, it is necessary to calculate the confidence interval in order to choose the appropriate sample size and balance the computational effort and precision; since the function is multivariate and complex, a nonparametric bootstrap method is used to calculate the confidence interval, that is, multiple resampling and taking the corresponding percentiles to obtain a 95% confidence interval;

[0112] First, we define percentiles. For a given probability p X , the Monte Carlo sampling mean calculation result of discrete random variable X p X The percentile is defined as:

[0113]

[0114] in: represents the domain of real numbers; x′ is any real number;

[0115] When using the Bootstrap method to calculate the confidence interval of the mean, the sampling results are calculated repeatedly several times to obtain the initial sample of data;

[0116] Step 10.1, Bootstrap sampling: randomly draw samples from the data, allowing repeated draws, and each bootstrap sample has the same size; repeat this process s times, where s is the number of bootstrap samples;

[0117] Step 10.2, calculate the sample mean: calculate the sample mean of each Bootstrap sample;

[0118] Step 10.3. Calculate the confidence interval: Determine the confidence level based on the distribution of the Bootstrap sample mean. To calculate the 95% confidence interval for the population mean, use the 2.5% percentile and 97.5% percentile of the sample mean. The confidence interval can be expressed as:

[0119]

[0120] Here, p represents the sample mean XPercentile; the calculation result with a sample size of 1000 shows that the length of the 95% confidence interval is 5%; this indicates that in this planning scenario, the reliability deviation of the probability that the source side meets the load side demand is relatively low;

[0121] Step 11: Investment decision analysis and reliability assessment methods

[0122] The Monte Carlo sampling method is used to evaluate the system cost under different sample sizes. The Pareto optimal solution set obtained by the NSGA-II algorithm is used to obtain the carbon emissions and costs of different solutions. Reliability assessment is performed corresponding to different renewable energy penetration rates REPR and the use of diesel generators, and the impact of REPR on system costs is analyzed. The optimal investment plan under different situations can be analyzed numerically. The economic feasibility of the system is evaluated based on the payback period, net present value, internal rate of return, and levelized electricity cost analysis of different REPRs to find the economic optimum.

[0123] Furthermore, in the structure of the off-grid integrated energy system, wind turbines, photovoltaics and energy storage systems are used as main modules, electric refrigerators and electric boilers provide cooling and heating loads, and diesel generators serve as auxiliary backup power sources.

[0124] Furthermore, in step 1.4, P in (t) is the charging power, which is positive for charging and negative for discharging.

[0125] Furthermore, in step 8, only the source-side scenario consisting of wind turbines and photovoltaics is considered. Other scenarios require reconsidering the modeling method of the source side, but the random processing method proposed in steps 6 and 7 of this method can be applied.

[0126] Furthermore, in step 11, in actual calculations, it is usually necessary to determine a sample size of 1000-5000 as a model calculation benchmark to ensure the source side uncertainty processing effect and calculation efficiency.

[0127] The beneficial effects of the present invention are:

[0128] This method proposes a stochastic optimization algorithm for off-grid integrated energy system capacity planning. Taking into account the probability density function of wind speed and solar irradiation, the expected value model is used as the objective function to transform it into a deterministic planning problem. The Monte Carlo method is used to quantify uncertainty and renewable energy penetration, and reliability assessment and operation prediction are integrated. Equipment constraints and energy balance are considered in the model. A multi-objective approach is used to simultaneously consider carbon emissions and economic factors, and investment decisions are made by evaluating renewable energy penetration.

[0129] The advantages of this method are: 1. Considering the uncertainty of renewable energy, it has stronger robustness than deterministic methods and effectively solves the system uncertainty caused by volatility. 2. It proposes a calculation method for the renewable energy penetration index and integrates this index with the Monte Carlo method for reliability analysis, which can provide a more detailed and efficient analysis method than other methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0130] Figure 1 This is a structural diagram of the off-grid energy system of the present invention. DETAILED DESCRIPTION

[0131] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0132] The present invention discloses an uncertainty planning algorithm for an off-grid wind-solar-storage system, comprising the following steps: constructing a system mathematical model by considering equipment constraints and power balance; considering the uncertainty of wind speed and light intensity and using parameter estimation methods to perform probability distribution fitting; constructing a random optimization algorithm using an expected value model and a Monte Carlo method; using a non-dominated sorting genetic algorithm (NSGA-II) to perform multi-objective optimization calculations of economic benefits and carbon emissions; and proposing a new method for evaluating renewable energy penetration combined with the Monte Carlo method definition, and using it as an indicator to provide investment decision-making recommendations and reliability evaluation. . In addition, the present invention verifies the feasibility and effectiveness of the model by applying it to a simulation of a real project, showing the improvement of economic benefits and the effective reduction of pollutant emissions.

[0133] Figure 1 The structure of the off-grid integrated energy system proposed by the present invention is shown, in which wind turbines, photovoltaics and energy storage systems are the main modules. Electric refrigerators and electric boilers provide cooling and heating loads, while diesel generators serve as auxiliary backup power sources.

[0134] Step 1: Construct a system mathematical model by considering equipment constraints and power balance:

[0135] Step 1.1: Establish the mathematical model of electric boiler (EB), a commonly used heating equipment. The mathematical model and operating constraints are as follows:

[0136]

[0137] Where: P EB (t) represents the power output of the electric boiler at time t (kW); H EB (t) represents the heating capacity of the electric boiler at time t (kW); d EB Represents the capacity of the electric boiler (kW); Indicates the maximum output power limit of the electric boiler (kW); Indicates the minimum output power limit of the electric boiler (kW); Inv EB represents the investment cost of electric boiler; TC EB represents the total cost of the electric boiler; p ele Indicates the electricity price (kWh); H EB_limit Indicates the output power adjustment limit of the electric boiler (kW / h); η EB Indicates the thermal energy conversion efficiency of the electric boiler.

[0138] Step 1.2: Establish the mathematical model of the electric refrigerator (EC), a commonly used refrigeration equipment. The mathematical model and operating constraints are as follows:

[0139]

[0140] Where: P EC (t) represents the power output of the electric refrigerator at time t (kW); C EC (t) represents the cooling capacity of the electric refrigerator at time t (kW); d EC Represents the capacity of the electric refrigerator (kW); Indicates the maximum output power limit of the electric refrigerator (kW); Indicates the minimum output power limit of the electric refrigerator (kW); Inv EC represents the investment cost of the electric refrigerator; TC EC represents the total cost of the electric refrigerator; p ele Indicates the electricity price (kWh); C EC_limit Indicates the output power adjustment limit of the electric refrigerator (kW / h); η EC Represents the energy conversion efficiency of the electric refrigerator.

[0141] Step 1.3: Establish the mathematical model of diesel generator (DG). Diesel generator (DG) is a reliable source of electricity used to power various applications. It generates electricity by using diesel fuel. The mathematical model and operating constraints are described as follows:

[0142]

[0143] fc(t) represents the fuel consumption of the diesel generator at time t (L); P DGN Indicates the rated power of the diesel generator (kW); P DG (t) represents the output power of the diesel generator at time t (kW); A DG represents the fuel intercept (L / kW); B DG represents the fuel slope coefficient (L / kW); p fuel Indicates fuel price ( / L); InvDG Represents the investment cost of diesel generator; TC DG represents the total cost of the diesel generator; Indicates the maximum output power limit of the diesel generator; Indicates the minimum output power limit of the diesel generator; CE DG represents the carbon emission of diesel generator (kg); η DG Represents the carbon emission rate of diesel generators.

[0144] Step 1.4: Mathematical model of electric energy storage system (ES). The mathematical model of electric energy storage system (ES) is represented based on the state of charge (SOC), which is described as follows:

[0145]

[0146] Where: P in (t) represents the charging power (charging is a positive value, discharging is a negative value) (kW); P out (t) represents the discharge power (kW); E stored (t) represents the state of charge of the stored energy at time t (kWh); d ES Indicates the capacity of the energy storage system (kW); E max Indicates the maximum capacity of electrical energy storage (kWh); E min Indicates the minimum capacity of electrical energy storage (kWh); η charge Indicates charging efficiency; η discharge Indicates discharge efficiency; P in,max and P out,max Represents the maximum limit of charging and discharging power (kW); P in,min and P out,min Represents the minimum limit of charging and discharging power (kW); E initial Indicates the initial state of charge of the energy storage (kWh); Inv ES Represents the investment cost of the energy storage system; TC ES Represents the total cost of the energy storage system.

[0147] Step 2: Establish a step-by-step model of carbon trading costs. This study uses a step-by-step model to calculate carbon trading costs, which defines different emission ranges. If carbon emissions are lower than the specified quota, the cost is negative, which indicates that energy suppliers can sell excess carbon emission rights in the carbon trading market and obtain technical subsidies. The transaction price increases as the emission range decreases. On the contrary, when emissions exceed the quota, the cost becomes positive, indicating that energy suppliers need to purchase additional carbon emission rights. At this time, a larger emission range corresponds to a higher transaction price. The model divides actual carbon emissions into linear segments for analysis:

[0148]

[0149] CE p =CE DG

[0150] In this formula: represents the carbon trading cost in the integrated energy system; μ represents the reward coefficient; ctp represents the carbon trading price in the market; θ represents the growth rate of the carbon trading price at each step; CE p Represents the total carbon emissions of the system equipment; CE c represents the system's free carbon emission quota; v represents the length of the carbon emission interval (Note: Nu, capital Ν, lowercase v, is the thirteenth Greek letter); CE DG Represents the carbon emissions of diesel generators.

[0151] Step 3: Establish the objective function of the integrated energy system planning model:

[0152] The proposed integrated energy system planning model aims to minimize the total investment cost and carbon trading cost. The total cost includes investment cost and operating cost. The multi-objective function is expressed as follows:

[0153]

[0154] Where N is the number of years of operation. r is the interest rate; F 1 represents the objective function of total investment cost; F 2 Represents the objective function of carbon trading cost; Min means the optimization target is the minimum value.

[0155] Step 4: Establish constraints for the integrated energy system planning model:

[0156] The constraints mainly include system operation constraints, system installation constraints and power balance constraints:

[0157] System operating constraints: defined by the mathematical models of various devices in step 1.

[0158] Installation constraint: Ensure that the capacity of the running device does not exceed the capacity limit of the available devices, expressed as:

[0159]

[0160] Where: M is the set of all device types, M = [EB, EC, ES, DG]; i represents the device type (conflicting with the following symbol); d i Indicates the capacity of the corresponding device; Indicates the maximum capacity of device type i.

[0161] Power balance constraint: includes the balance of electricity, heat energy and gas energy, expressed as follows:

[0162]

[0163] Where: P WD (t) represents the output power of the fan at time t; P PV (t) represents the output power of the photovoltaic at time t; H EB (t) represents the output heat energy of the electric boiler (EB) at time t; C EC (t) represents the output cooling power of the electric refrigerator (EC) at time t; L ELOAD (t) represents the power load at time t; L HLOAD (t) represents the heat load at time t; L CLOAD (t) represents the cooling load at time t; ξ represents the adjustment factor, which is usually equal to 0.

[0164] Step 5: Determine the optimization variables of the integrated energy system planning model and the mathematical expression of the planning problem:

[0165] The investment decision variable is denoted as x, while the operational optimization variable is denoted as y. They are expressed as follows:

[0166]

[0167] Where: d WT Indicates the capacity of the fan; d PV In summary, the general form of the planning problem proposed by the present invention is as follows:

[0168]

[0169] Among them: F(x) represents the objective function set; Min means that the optimization target is the minimum value; G i (x) represents the inequality constraint; J i (x) represents an equality constraint; m represents the number of inequality constraints; n represents the total number of constraints; and i represents the constraint number.

[0170] Step 6: Consider the uncertainty of wind speed and sunlight intensity and use parameter estimation method to fit the probability distribution to study the renewable energy penetration rate (REPR) of the system;

[0171] How to deal with uncertainty:

[0172] Random distribution prediction can provide more uncertainty information for the operation of integrated energy and power systems. The present invention uses a parameter estimation method. It presupposes the distribution form, performs statistical analysis on historical data, estimates the specific distribution form through parameter estimation based on historical data, and finally optimizes the objective function of the expected value model by generating multiple scenarios. In order to study the renewable energy penetration rate (REPR) of the system by quantifying uncertainty and characterizing economic feasibility. The uncertainty processing method model is processed in the following steps:

[0173] Step 6.1: Use probability density function to represent the uncertain renewable energy generation. Then, based on historical statistical data, parameter estimation is proposed to obtain the corresponding distribution characteristics.

[0174] Step 6.2: Construct a random expectation model to represent the uncertainty of the output in the form of a set of scenarios. Use Monte Carlo simulation to generate an accurate set of output scenarios.

[0175] Step 6.3: Consider the constraints and optimization variables to obtain the planning problem of the off-grid integrated energy system. Finally, a multi-objective algorithm is used to solve the problem.

[0176] Historical statistical data is usually used to fit the distribution characteristics of uncertain renewable energy generation. This process aims to derive a probability density function or a cumulative distribution function. Commonly used probability models for simulating wind speed and solar irradiation include Weibull distribution and Beta distribution. Therefore, the present invention uses Weibull distribution and Beta distribution to simulate wind speed and solar irradiation. Standardized solar irradiance G p =G / G max Following the beta distribution model, its probability density function f(G p |α,β) is expressed as follows:

[0177]

[0178] Where G represents the actual irradiance, the unit is kilowatt per square meter kW / m 2 ; G max is the maximum irradiance, also in kilowatts per square meter kW / m 2 ; α and β are the two shape parameters of the Beta distribution; B represents the Beta function. These parameters are affected by the average value of the irradiance, μ p and standard deviation Influence of factors such as:

[0179]

[0180] The photovoltaic system output power is given by the following formula:

[0181]

[0182] Where: P pv Indicates the actual output power of the photovoltaic panel, in kilowatts (kW); P st It is the maximum output power under standard test conditions, in kilowatts (kW); G st represents standard irradiation conditions, in kilowatts per square meter (kW / m^2); η pv Represents the efficiency of the photovoltaic panel; η pvs Indicates the efficiency of the photovoltaic system.

[0183] The wind speed change v (note: English letter v) follows the Weibull distribution, and its probability density function f(v|c,λ) is expressed as follows:

[0184]

[0185] Where c is the shape parameter and λ is the scale parameter. In the case of wind speed distribution, c determines the shape of the distribution: it indicates that higher or lower wind speeds are less common, with moderate wind speeds dominating most of the time. As c increases, the tail of the distribution becomes heavier, indicating that higher wind speeds are more likely and lower wind speeds are less likely. The scale parameter λ controls the scale or amplitude of the wind speeds. Specifically, larger values ​​of λ indicate a wider range of wind speed amplitudes, which may have higher peaks and a wider range of fluctuations. Conversely, smaller values ​​of λ indicate a narrower range of wind speed amplitudes, indicating that the changes in wind speed are less dramatic. These parameters are affected by, for example, the average value of the irradiance, μ w and standard deviation σ w Influence of factors such as:

[0186]

[0187] Where: Γ represents the Gamma function. The fan system output power is given by the following formula:

[0188]

[0189] Where: P wt represents the actual output power of the wind turbine components in kilowatts (kW); v represents the wind speed in meters per second (m / s); v in is the starting wind speed, i.e. the minimum wind speed at which the wind turbine starts generating electricity, expressed in meters per second (m / s); v out is the cut-out wind speed, i.e. the maximum wind speed at which the wind turbine stops generating electricity to prevent mechanical damage, in meters per second (m / s); P n is the rated power under standard conditions, in kilowatts (kW); η wt Represents the efficiency of the wind turbine.

[0190] Likelihood estimation is an important method for estimating the distribution parameters of a data set for subsequent analysis and modeling. This paper will use the likelihood estimation method to analyze the extracted irradiance and wind speed data. After preprocessing, these data are divided into multiple data points per hour for the whole year, represented by PD i , where i represents the index of the illumination data point. In order to estimate the distribution parameters of the data, we choose to use the negative log-likelihood function. The core goal of this function is to evaluate the degree of fit between the data and a specific probability distribution. Specifically, the mathematical definition of the negative log-likelihood function is as follows:

[0191]

[0192] Where: α l ,β l ∈Ω;Ω={(ξα l ,β l )|α l >0,β l >0}. In the above equation, f(PD i |α l ,β l ) indicates that the parameter α l ,β l Data point PD under the condition i The probability density function of . By calculating the negative natural logarithm of the probability density function for each data point and summing over all data points, we obtain the negative log-likelihood value NL. In order to find the parameter value that best explains the data distribution, the fitting error between the data and the selected probability distribution must be minimized. By minimizing the negative log-likelihood function NL, we seek to find the distribution parameter α that best fits the data l ,β l .

[0193] Step 7: This method uses multiple new energy output scenarios and their corresponding probabilities to describe the uncertainty of output. A random expectation scheduling model is established, whose goal is to minimize the expected value of the multi-objective function:

[0194] F=E[F 1 ,F 2 ]

[0195] Where F represents the expected function of the random programming objective; E represents the mathematical expectation value.

[0196] The random expected value model constructed by this method characterizes the uncertainty of new energy output in the form of a scenario set. Generating an accurate output scenario set is the key and prerequisite for efficiently solving the model. The present invention uses Monte Carlo simulation to generate the output scenario set. Due to the uncertainty of wind speed and solar radiation, the optimal cost of the expected value model is a function f(v,G p ), where G pfollows the beta distribution, and v follows the weibull distribution. When f(v,G p ) contains two different random variables v and G p , respectively follow the probability density function f v (v) and When we use the Monte Carlo method to estimate the expected value of the function E(f(v,G p )) Generate N′ groups of independent samples where v i is from f v (v) the sample drawn, is from The sample is drawn, where i represents the sample number. We can calculate the sample mean estimate

[0197]

[0198] Next, we calculate the expected value of this estimate Since v,G p are independent, we have:

[0199]

[0200] Since each sample are drawn independently, we can calculate each expected value separately:

[0201]

[0202] Sample mean estimation The expected value of is equal to the true expected value. To calculate the variance of the sample mean estimate

[0203]

[0204] Due to x i and i are independent, and taking into account the properties of variance, we have the following relationship:

[0205]

[0206] Therefore, when the sample size N' tends to infinity, tends to zero. This means that as the sample size increases, the variance of the sample mean estimate decreases, resulting in a more stable estimate that is closer to the true value. Therefore, when using the Monte Carlo method to estimate the sample mean, choosing a sufficiently large N′ is the key to ensuring the accuracy of the calculation.

[0207] Step 8: Monte Carlo simulation to solve the expected value model

[0208] The present invention only considers the source-side scenario consisting of wind turbines and photovoltaics. Other scenarios need to reconsider the modeling method of the source side, but steps 6 and 7 proposed in the present invention can be applied. The present invention uses Monte Carlo simulation for calculation to determine the optimal solution. By using this method, the relationship between renewable energy penetration rate (REPR) and cost is revealed, thereby providing advice for investment decisions. In addition, by quantifying REPR and uncertainty, the economic feasibility of the system is characterized. Compared with directly calculating a certain proportion of renewable energy supply to the total load, the present invention proposes to use a probabilistic description of renewable energy penetration. Based on the original load data, it is assumed that renewable energy meets the load with a certain probability. During the sampling calculation process, days with no wind or no sunshine and high load are conveniently screened out, and diesel power generation is used as a backup. The probabilistic form of renewable energy penetration is defined as:

[0209] REPR=Pr(P DG (t)=0,ξ≤0)

[0210] Where Pr represents probability. REPR can be understood as an event (P DG (t)=0,ξ≤0) occurs, and we define X to follow a 0-1 distribution. If the event occurs, then X=1, and we can get: REPR=Pr(X=1). According to the expected value calculation formula of the 0-1 distribution, we can get: REPR=E(X). To calculate REPR, first, considering a fairly large sample size (N′), we apply the law of large numbers. For all

[0211]

[0212] Where: X 1 ,…,X N′ is a 0-1 random variable that is independent and identically distributed in X, and ∈ is an arbitrary positive real number.

[0213] Therefore, the present invention adopts the sample mean Approximately estimating E(X), during the Monte Carlo simulation scenario generation process, the scenario that successfully meets the load will be close to N′·REPR. In the N′·(1-REPR) calculations, the constraints that do not meet the power balance condition are allowed to achieve a probabilistic model. Therefore, when using the Monte Carlo method for simulation, the choice of sample size greatly affects the accuracy of the probability calculation of the source output meeting the load.

[0214] Step 9: Find the optimal algorithm for the planning problem

[0215] The proposed model takes into account both carbon emissions and economic costs, so a multi-objective solution algorithm is required. In order to ensure the practicality and accuracy of the model solution, NSGA-II is used to solve the model. NSGA-II operates based on the non-dominated sorting principle and uses an elite strategy to divide the population into different non-dominated levels. Compared with other algorithms, NSGA-II has powerful search capabilities and high parallelization characteristics, making it particularly suitable for dealing with nonlinear multi-objective problems.

[0216] Step 10: Calculate confidence interval using bootstrap method

[0217] Therefore, it is necessary to calculate the confidence interval to select the appropriate sample size and balance the computational effort and precision. Since the function is a multivariate and complex function, the present invention uses non-parametric bootstrap methods to calculate the confidence interval, i.e., multiple resampling and taking the corresponding percentiles to obtain a 95% confidence interval.

[0218] First, we define percentiles. For a given probability p X , the Monte Carlo sampling mean calculation result of discrete random variable X p X The percentile is defined as:

[0219]

[0220] in: represents the field of real numbers; x′ is any real number.

[0221] When the Bootstrap method is used to calculate the confidence interval of the mean, the sampling results are calculated repeatedly several times to obtain the initial sample of data.

[0222] 1. Bootstrap sampling: randomly draw samples from the data, allowing repeated draws, and each bootstrap sample has the same size. Repeat this process s times, where s is the number of bootstrap samples.

[0223] 2. Calculate the sample mean: Calculate the sample mean of each Bootstrap sample.

[0224] 3. Calculate the confidence interval: Determine the confidence level based on the distribution of the Bootstrap sample mean. To calculate the 95% confidence interval for the population mean, use the 2.5% percentile and 97.5% percentile of the sample mean. The confidence interval can be expressed as:

[0225]

[0226] Here, p represents the sample mean X Percentile. The calculation result with a sample size of 1000 shows that the length of the 95% confidence interval is 5%. This indicates that in this planning scenario, the reliability deviation of the probability that the source side will meet the load side demand is relatively low.

[0227] Step 11: Investment decision analysis and reliability assessment methods

[0228] The Monte Carlo sampling method is used to evaluate the system cost under different sample sizes. In actual calculations, it is usually necessary to determine a sample size of 1000-5000 as the model calculation benchmark to ensure the uncertainty processing effect and calculation efficiency on the source side. The Pareto optimal solution set obtained by the NSGA-II algorithm is used to obtain the carbon emissions and costs of different solutions. Reliability evaluation is performed corresponding to different renewable energy penetration rates (REPR) and the use of diesel generators, and the impact of REPR on system costs is analyzed. The optimal investment plan under different circumstances can also be analyzed numerically. The economic feasibility of the system is evaluated based on the payback period, net present value, internal rate of return, and levelized electricity cost analysis of different REPRs to find the economic optimum.

[0229] This invention was funded by the Shanghai Youth Science and Technology Talent Yangfan Program (22YF1414900).

[0230] The above-mentioned embodiment only expresses one implementation mode of the present invention, and its description is relatively specific and detailed, but it cannot be understood as limiting the scope of the invention patent. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be based on the attached claims.

Claims

1. A random optimization method for off-grid wind-solar-storage system capacity configuration, characterized in that: The following steps are involved: A mathematical model of the system was constructed by considering equipment constraints and power balance; uncertainty in wind speed and light intensity was considered and parameter estimation methods were used to fit probability distributions; a stochastic optimization algorithm was constructed using the expected value model and the Monte Carlo method; a non-dominated sorting genetic algorithm was used to perform multi-objective optimization calculations of economic benefits and carbon emissions; and a new method for evaluating renewable energy penetration was proposed, combining the Monte Carlo method definition, and using it as an indicator to provide investment decision recommendations and reliability assessments. Step 1: Construct a system mathematical model by considering equipment constraints and power balance: Step 2: Establish a step-by-step model for carbon trading costs. A step-by-step model is used to calculate carbon trading costs. The model defines different emission ranges. If carbon emissions are lower than the specified quota, the cost is negative, indicating that energy suppliers can sell excess carbon emission rights in the carbon trading market and obtain technology subsidies. The transaction price increases as the emission range decreases. Conversely, when emissions exceed the quota, the cost becomes positive, indicating that energy suppliers need to purchase additional carbon emission rights. Step 3: Establish the objective function of the integrated energy system planning model: The integrated energy system planning model aims to minimize the total investment cost and carbon trading cost; the total cost includes investment cost and operating cost; Step 4: Establish constraints for the integrated energy system planning model: Constraints mainly include system operation constraints, system installation constraints, and power balance constraints: Step 5: Determine the optimization variables of the integrated energy system planning model and the mathematical expression of the planning problem: Step 6: Consider the uncertainty of wind speed and sunlight intensity and use parameter estimation method to perform probability distribution fitting to study the renewable energy penetration rate REPR of the system; Uncertainty handling methods: Random distribution forecasting can provide more uncertainty information for the operation of integrated energy and power systems. A parameter estimation method is used. It presupposes the distribution form, performs statistical analysis on historical data, estimates the specific distribution form through parameter estimation based on historical data, and ultimately optimizes the objective function of the expected value model by generating multiple scenarios. This method aims to study the renewable energy penetration rate (REPR) of the system by quantifying uncertainty and characterizing economic feasibility. Step 7: Use multiple new energy output scenarios and their corresponding probabilities to describe output uncertainty; establish a stochastic expectation scheduling model whose goal is to minimize the expected value of the multi-objective function; Step 8: Monte Carlo simulation to solve the expected value model Only the source-side scenario consisting of wind turbines and photovoltaics was considered, and Monte Carlo simulation was used to determine the optimal solution; Step 9: Find the optimal algorithm for the planning problem The proposed model considers both carbon emissions and economic costs, so a multi-objective solution algorithm is required. To ensure the practicality and accuracy of the model solution, NSGA-II is used to solve the model. Step 10: Calculate confidence interval using the bootstrap method Select an appropriate sample size to balance computational effort and precision; since the function is multivariate and complex, a nonparametric bootstrap method is used to calculate the confidence interval; Step 11: Investment Decision Analysis and Reliability Assessment Methods The Monte Carlo sampling method is used to evaluate the system cost under different sample sizes. The Pareto optimal solution set obtained by the NSGA-II algorithm is used to obtain the carbon emissions and costs of different solutions. Reliability assessment is performed for different renewable energy penetration rates (REPR) and diesel generator usage, and the impact of REPR on system costs is analyzed. The optimal investment plan under different situations can be analyzed numerically. The economic feasibility of the system is evaluated based on the payback period, net present value, internal rate of return, and levelized cost of electricity analysis based on different REPRs, and the economic optimum is found.

2. The stochastic optimization method for off-grid wind-solar-storage system capacity configuration according to claim 1 is characterized in that: The specific steps include: Step 1: Construct a system mathematical model by considering equipment constraints and power balance: Step 1.1: Establish the mathematical model of the electric boiler EB, a commonly used heating equipment. The mathematical model and operating constraints are as follows: Where: P EB (t) represents the power output of the electric boiler at time t (kW); H EB (t) represents the heating capacity of the electric boiler at time t (kW); d EB Represents the capacity of the electric boiler (kW); Indicates the maximum output power limit of the electric boiler (kW); Indicates the minimum output power limit of the electric boiler (kW); Inv EB represents the investment cost of electric boiler; TC EB represents the total cost of the electric boiler; p ele Indicates electricity price (kWh); H EB_limit Indicates the output power adjustment limit of the electric boiler (kW / h); η EB Indicates the thermal energy conversion efficiency of electric boiler; Step 1.2: Build a mathematical model for an electric refrigerator (EC), a commonly used refrigeration device. The mathematical model and operating constraints are as follows: Where: P EC (t) represents the power output of the electric refrigerator at time t (kW); C EC (t) represents the cooling capacity of the electric refrigerator at time t (kW); d EC Represents the capacity of the electric refrigerator (kW); Indicates the maximum output power limit of the electric refrigerator (kW); Indicates the minimum output power limit of the electric refrigerator (kW); Inv EC represents the investment cost of the electric refrigerator; TC EC represents the total cost of the electric refrigerator; p ele Indicates electricity price (kWh); C EC_limit Indicates the output power adjustment limit of the electric refrigerator (kW / h); η EC Indicates the energy conversion efficiency of the electric refrigerator; Step 1.3: Establish a mathematical model for the diesel generator (DG). The diesel generator (DG) is a reliable source of electricity used to power various applications. It generates electricity by using diesel fuel. The mathematical model and operating constraints are as follows: fc(t) represents the fuel consumption of the diesel generator at time t (L); P DGN Indicates the rated power of the diesel generator (kW); P DG (t) represents the output power of the diesel generator at time t (kW); A DG represents fuel intercept (L / kW); B DG represents the fuel slope coefficient (L / kW); p fuel Indicates fuel price ( / L); Inv DG Indicates the investment cost of diesel generator; TC DG represents the total cost of the diesel generator; Indicates the maximum output power limit of the diesel generator; Indicates the minimum output power limit of the diesel generator; CE DG represents the carbon emissions of diesel generators (kg); η DG represents the carbon emission rate of diesel generators; Step 1.4: Mathematical model of the electric energy storage system ES. The mathematical model of the electric energy storage system ES is represented based on the state of charge (SOC) and is described as follows: Where: P in (t) represents the charging power (kW); P out (t) represents the discharge power (kW); E stored (t) represents the state of charge of the stored energy at time t (kWh); d ES Indicates the capacity of the energy storage system (kW); E max Indicates the maximum capacity of electrical energy storage (kWh); E min Indicates the minimum capacity of electrical energy storage (kWh); η charge Indicates charging efficiency; η discharge Indicates discharge efficiency; P in,max and P out,max Represents the maximum limit of charging and discharging power (kW); P in,min and P out,min Represents the minimum limit of charging and discharging power (kW); E initial Indicates the initial state of charge of the energy storage (kWh); Inv ES represents the investment cost of the energy storage system; TC ES represents the total cost of the energy storage system; Step 2: Establish a stepwise model for carbon trading costs. A stepwise model is used to calculate carbon trading costs. The model defines different emission ranges. If carbon emissions are lower than the specified quota, the cost is negative, indicating that energy suppliers can sell excess carbon emission rights in the carbon trading market and obtain technical subsidies. The transaction price increases as the emission range decreases. Conversely, when emissions exceed the quota, the cost becomes positive, indicating that energy suppliers need to purchase additional carbon emission rights. At this time, a larger emission range corresponds to a higher transaction price. The model divides actual carbon emissions into linear segments for analysis: WHAT p =EC DG In this formula: represents the carbon transaction cost in the integrated energy system; μ represents the reward coefficient; ctp represents the carbon transaction price in the market; θ represents the growth rate of the carbon transaction price at each step; CE p Represents the total carbon emissions of the system equipment; CE c represents the system's free carbon emission quota; ν represents the length of the carbon emission interval; CE DG represents the carbon emissions of diesel generators; Step 3: Establish the objective function of the integrated energy system planning model: The integrated energy system planning model aims to minimize the total investment cost and carbon trading cost; the total cost includes investment cost and operating cost; the multi-objective function is expressed as follows: Where N represents the number of years of operation; r represents the interest rate; F1 represents the objective function of total investment cost; F2 represents the objective function of carbon trading cost; Min represents that the optimization target is the minimum value; Step 4: Establish constraints for the integrated energy system planning model: Constraints mainly include system operation constraints, system installation constraints, and power balance constraints: System operation constraints: defined by the mathematical models of various devices in step 1; Installation constraints: Ensure that the capacity of the running device does not exceed the capacity limit of the available devices, expressed as: Where: M is the set of all device types, M = [EB, EC, ES, DG]; i represents the device type; d i Indicates the capacity of the corresponding device; Indicates the maximum capacity limit of device type i; Power balance constraint: includes the balance of electricity, heat energy, and gas energy, expressed as follows: Where: P WD (t) represents the output power of the fan at time t; P PV (t) represents the output power of photovoltaic at time t; H EB (t) represents the output heat energy of the electric boiler (EB) at time t; C EC (t) represents the output cooling power of the electric refrigerator (EC) at time t; L ELOAD (t) represents the power load at time t; L HLOAD (t) represents the heat energy load at time t; L CLOAD (t) represents the cooling load at time t; ξ represents the adjustment factor, which is usually equal to 0; Step 5: Determine the optimization variables of the integrated energy system planning model and the mathematical expression of the planning problem: The investment decision variable is denoted as x, and the operational optimization variable is denoted as y; they are expressed as follows: Where: d WT Indicates the capacity of the fan; d PV represents the capacity of photovoltaic; in summary, the general form of the planning problem is as follows: Among them: F(x) represents the objective function set; Min represents the optimization target is the minimum value; G i (x) represents the inequality constraint; J i (x) represents the equality constraint; m represents the number of inequality constraints; n represents the total number of constraints; i represents the constraint number; Step 6: Consider the uncertainty of wind speed and sunlight intensity and use parameter estimation method to perform probability distribution fitting to study the renewable energy penetration rate REPR of the system; Uncertainty handling methods: Random distribution prediction can provide more uncertainty information for the operation of integrated energy and power systems. Parameter estimation methods are used. They presuppose the distribution form, perform statistical analysis on historical data, estimate the specific distribution form through parameter estimation based on historical data, and ultimately optimize the objective function of the expected value model by generating multiple scenarios. In order to study the renewable energy penetration rate (REPR) of the system by quantifying uncertainty and characterizing economic feasibility, the uncertainty processing method model is processed in the following steps: Step 6.1: Use a probability density function to represent the uncertain renewable energy generation; then, based on historical statistical data, propose parameter estimates to obtain the corresponding distribution characteristics; Step 6.2: Construct a stochastic expectation model to represent the uncertainty of the output in the form of a set of scenarios; use Monte Carlo simulation to generate an accurate set of output scenarios; Step 6.3: Consider the constraints and optimization variables to obtain the planning problem of the off-grid integrated energy system; finally, use a multi-objective algorithm to solve the problem; Historical statistical data is often used to fit the distribution characteristics of uncertain renewable energy generation. This process aims to derive a probability density function or cumulative distribution function. Commonly used probability models for simulating wind speed and solar irradiation include the Weibull distribution and the Beta distribution. Therefore, the Weibull distribution and the Beta distribution are used to simulate wind speed and solar irradiation. The standardized solar irradiation G p =G / G max Following the beta distribution model, its probability density function f(G p |α,β) are expressed as follows: Where G represents the actual irradiance, the unit is kilowatt per square meter kW / m 2 ; G max is the maximum irradiance, also in kilowatts per square meter kW / m 2 α and β are the two shape parameters of the Beta distribution; B represents the Beta function; these parameters are affected by the mean value of the irradiance μ p and standard deviation Factors affecting: The output power of the photovoltaic system is given by the following formula: Where: P pv Indicates the actual output power of the photovoltaic panel, in kilowatts (kW); P st It is the maximum output power under standard test conditions, in kilowatts (kW); G st represents standard irradiation conditions, in kilowatts per square meter (kW / m^2); η pv Represents the efficiency of the photovoltaic panel; η pvs Indicates the efficiency of the photovoltaic system; The wind speed variation v follows the Weibull distribution, and its probability density function f(v|c,λ) is expressed as follows: where c is the shape parameter and λ is the scale parameter; in the case of wind speed distribution, c determines the shape of the distribution: it indicates that higher or lower wind speeds are less common, with moderate wind speeds dominating in most cases; as c increases, the tail of the distribution becomes heavier, indicating that higher wind speeds are more likely and lower wind speeds are less likely; the scale parameter λ controls the scale or amplitude of the wind speeds; specifically, larger values ​​of λ indicate a wider range of wind speed amplitudes, potentially with higher peaks and a wider range of fluctuations; conversely, smaller values ​​of λ indicate a narrower range of wind speed amplitudes, indicating that wind speed variations are less dramatic; these parameters are affected by, for example, the average value of the irradiance, μ w and standard deviation σ w Factors affecting: Where: Γ represents the Gamma function; the wind turbine system output power is given by the following formula: Where: P wt represents the actual output power of the wind turbine components in kilowatts (kW); v represents the wind speed in meters per second (m / s); v in is the starting wind speed, i.e. the minimum wind speed at which the wind turbine starts generating electricity, in meters per second (m / s); v out is the cut-out wind speed, i.e. the maximum wind speed at which the wind turbine stops generating electricity to prevent mechanical damage, in meters per second (m / s); P n is the rated power under standard conditions, in kilowatts (kW); η wt Indicates the efficiency of the wind turbine; Likelihood estimation is an important method for estimating the distribution parameters of a dataset for subsequent analysis and modeling. The extracted irradiance and wind speed data will be analyzed using the likelihood estimation method. After preprocessing, these data are divided into multiple data points for each hour of the entire year, represented as PD. i , where i represents the index of the illumination data point; in order to estimate the distribution parameters of the data, the negative log-likelihood function is used; the core goal of this function is to evaluate the degree of fit between the data and a specific probability distribution; specifically, the mathematical definition of the negative log-likelihood function is as follows: Where: α l ,β l ∈Ω;Ω={(α l ,β l )|α l >0,β l >0}; In the above equation, f(PD i |α l ,β l ) indicates that the parameter α l ,β l Conditional data point PD i The probability density function of the data point is obtained by calculating the negative natural logarithm of the probability density function of each data point and summing it over all data points to obtain the negative log-likelihood value NL. In order to find the parameter value that best explains the data distribution, the fitting error between the data and the selected probability distribution must be minimized. By minimizing the negative log-likelihood function NL, we seek to find the distribution parameter α that best fits the data. l ,β l ; Step 7: Use multiple new energy output scenarios and their corresponding probabilities to describe the uncertainty of output; establish a random expectation scheduling model whose goal is to minimize the expected value of the multi-objective function: F=E[F1,F2] Where F represents the expected function of the stochastic programming objective; E represents the mathematical expectation value; The random expected value model constructed by this method represents the uncertainty of renewable energy output in the form of a scenario set; generating an accurate output scenario set is the key and prerequisite for solving the model efficiently; Monte Carlo simulation is used to generate the output scenario set; due to the uncertainty of wind speed and solar irradiation, the optimal cost of the expected value model is a function f(v,G p ), where G p Follows beta distribution, v follows weibull distribution; when f(v,G p ) contains two different random variables v and G p , respectively follow the probability density function f v (v) and When the Monte Carlo method is used to estimate the expected value of the function E(f(v,G p )); Generate N′ groups of independent samples where v i It is from f v (v) the sample drawn, It is from The sample is drawn, where i represents the sample number; the sample mean estimate can be calculated Next, calculate the expected value of this estimate Since v,G p is independent and has: Since each sample are drawn independently, and each expected value can be calculated separately: Sample mean estimation The expected value of is equal to the true expected value; in order to calculate the variance of the sample mean estimate Since x i and y i are independent, and taking into account the properties of the variance, the following relationship holds: Therefore, when the sample size N' tends to infinity, tends to zero; this means that as the sample size increases, the variance of the sample mean estimate decreases, resulting in a more stable estimate that is closer to the true value; therefore, when using the Monte Carlo method to estimate the sample mean, choosing a sufficiently large N′ is key to ensuring calculation accuracy; Step 8: Monte Carlo simulation to solve the expected value model Considering only the source-side scenario consisting of wind turbines and photovoltaics, Monte Carlo simulation was used to determine the optimal solution. This method revealed the relationship between renewable energy penetration rate (REPR) and cost, providing advice for investment decisions. Furthermore, by quantifying REPR and uncertainty, the economic feasibility of the system was characterized. Compared with directly calculating a certain proportion of renewable energy supply to the total load, a probabilistic description of renewable energy penetration was used. Based on the raw load data, it was assumed that renewable energy could meet the load with a certain probability. During the sampling calculation process, days with no wind or sunshine and high load were easily screened out, and diesel power generation was used as a backup. The probabilistic form of renewable energy penetration is defined as: REPR=Pr(P DG (t)=0,ξ≤0) Where Pr represents probability; REPR can be understood as the event (P DG (t)=0,ξ≤0) occurs, we define X to follow the 0-1 distribution. If the event occurs, then X=1, we can get: REPR=Pr(X=1); According to the expected value calculation formula of the 0-1 distribution, we can get: REPR=E(X); To calculate REPR, first, considering a fairly large sample size (N′), apply the law of large numbers, for all Where: X1,…,X N′ is a 0-1 random variable that is independent and identically distributed in X, and ∈ is any positive real number; Therefore, the sample mean Approximately estimating E(X), during the Monte Carlo simulation scenario generation process, the scenario that successfully meets the load will be close to N′·REPR; in N′·(1-REPR) calculations, constraints that do not meet the power balance condition are allowed, thus realizing a probabilistic model; therefore, when using the Monte Carlo method for simulation, the choice of sample size greatly affects the accuracy of the probability calculation of the source output meeting the load; Step 9: Find the optimal algorithm for the planning problem The proposed model considers both carbon emissions and economic costs, necessitating the use of a multi-objective solution algorithm. To ensure the practicality and accuracy of the model solution, NSGA-II was employed to solve the model. NSGA-II operates on the non-dominated sorting principle, employing an elitist strategy to divide the population into different non-dominated hierarchies. Compared to other algorithms, NSGA-II has powerful search capabilities and a high degree of parallelization, making it particularly well-suited for handling nonlinear multi-objective problems. Step 10: Calculate confidence interval using the bootstrap method Therefore, it is necessary to calculate the confidence interval in order to choose an appropriate sample size and balance the computational effort and precision. Since the function is multivariate and complex, a nonparametric bootstrap method is used to calculate the confidence interval, that is, multiple resampling and taking the corresponding percentiles to obtain a 95% confidence interval. First, define the percentile, for a given probability p X , the Monte Carlo sampling mean calculation result of discrete random variable X p X Percentiles are defined as: in: represents the field of real numbers; x′ is any real number; When using the Bootstrap method to calculate the confidence interval of the mean, the sampling results are calculated repeatedly several times to obtain the initial sample of data; Step 10.

1. Bootstrap sampling: Draw random samples from the data, allowing for repeated draws, with each bootstrap sample of the same size. Repeat this process s times, where s is the number of bootstrap samples. Step 10.2, Calculate the sample mean: Calculate the sample mean of each bootstrap sample; Step 10.

3. Calculate the confidence interval: Determine the confidence level based on the distribution of the bootstrap sample mean. To calculate the 95% confidence interval for the population mean, use the 2.5% and 97.5% percentiles of the sample mean. The confidence interval can be expressed as: Here, p represents the sample mean X Percentile; the calculation result with a sample size of 1000 shows that the length of the 95% confidence interval is 5%; this indicates that in this planning scenario, the reliability deviation of the probability that the source side will meet the load side demand is relatively low; Step 11: Investment Decision Analysis and Reliability Assessment Methods The Monte Carlo sampling method is used to evaluate the system cost under different sample sizes. The Pareto optimal solution set obtained by the NSGA-II algorithm is used to obtain the carbon emissions and costs of different solutions. Reliability assessment is performed for different renewable energy penetration rates (REPR) and diesel generator usage, and the impact of REPR on system costs is analyzed. The optimal investment plan under different situations can be analyzed numerically. The economic feasibility of the system is evaluated based on the payback period, net present value, internal rate of return, and levelized cost of electricity analysis based on different REPRs, and the economic optimum is found.

3. The stochastic optimization method for off-grid wind-solar-storage system capacity configuration according to claim 2 is characterized in that: In the structure of the off-grid integrated energy system, wind turbines, photovoltaics and energy storage systems serve as the main modules, electric refrigerators and electric boilers provide cooling and heating loads, and diesel generators serve as auxiliary backup power sources.

4. The stochastic optimization method for off-grid wind-solar-storage system capacity configuration according to claim 2 is characterized in that: In step 1.4, P in (t) is the charging power, which is positive for charging and negative for discharging.

5. The stochastic optimization method for off-grid wind-solar-storage system capacity configuration according to claim 2 is characterized in that: In step 8, only the source-side scenario consisting of wind turbines and photovoltaics is considered. Other scenarios require reconsidering the source-side modeling method, but the random processing method proposed in steps 6 and 7 of this method can be applied.

6. The stochastic optimization method for off-grid wind-solar-storage system capacity configuration according to claim 2, characterized in that: In step 11, in actual calculations, it is usually necessary to determine a sample size of 1000-5000 as the model calculation benchmark to ensure the source side uncertainty processing effect and calculation efficiency.

Images