The CHP-VPP in this paper mainly includes distributed power/heat output module and carbon cycle module. The distributed power/heat output module includes distributed wind power and photovoltaic, electric boiler, controllable load, and the power storage device. Carbon cycle module mainly includes the gas CHP unit, GPPCC, P2G, and gas storage device. The VPP realizes the recycling of CO₂ through GPPCC and P2G. GPPCC captures CO₂ generated by the CHP unit, and P2G converts CO₂ into CH₄. The carbon storage and hydrogen storage devices can be used to store excess CO₂ and H₂ at a certain time, so as to decouple carbon capture and electric conversion process. The electric boiler can use the surplus renewable energy to generate electricity to supply heat for the system, reduces the dependence of VPP on the heat output of CHP unit, and increases the flexibility of CHP unit operation. The controllable load and power storage device can cut peak and fill valley, and provide spare output for VPP. Figure 1 shows the energy flow diagram of CHP-VPP.

VPP coordination control center conducts information interaction with each unit in the VPP through communication technology, so that it can sense the operation status of each device and issue dispatching instructions to each unit. On this basis, the VPP forecasts the WPP and PV output of the next day, and then consider the operating status of each unit, the demand for VPP internal electrical load and thermal load, and formulate the next day’s operation plan of each unit and form the next day’s electricity purchase and sale strategy in the public grid.

Wind and photovoltaic power generation often have strong uncertainty, which will bring risks to the real-time operation of VPPs. Since the uncertainty of new energy output mainly comes from the prediction error, this paper constructs the joint probability distribution function according to the correlation of wind power and photoelectric output error. Then, the inverse transformation method is adopted to generate typical scenarios of wind-photoelectric output, and the random model is transformed into a deterministic model through the generation of uncertainty scenarios while retaining the wind-photoelectric output correlation. In order to take into account the randomness and correlation of the scene output at each moment, the scene output scene is generated.

(1) Constructing the covariance matrix σ 24 × 24 of the full cycle wind and solar forecast error, as follows:

where σ i j represents the covariance period i and period j of time t and ε is the covariance key parameter, which is used to control the correlation strength.

(2) The multivariate normal distribution Z 1 × 24 ∼ N 0 , σ 24 × 24 of the prediction error of full cycle scenery is constructed, and each random variable follows the standard normal distribution. Then, the m v n r n d function is called in MATLAB to randomly generate N samples.

Then, in order to reduce the amount of computation, k-means clustering is used to reduce scenes N to typical scenes n. The specific steps are as follows:

(1) Initial cluster centers D i 0 i = 1,2 , ⋯ , n are randomly generated within the value range of the n above N scenarios.

Finally, the typical output scenarios for wind and PV are combined to obtain the final typical output scenario n for wind and PV. The flow chart of n 2 uncertainty scenario generation in this paper is shown in Figure 3.

To improve the economy, promote the low-carbon development of VPP, and respond to the national call for “double carbon”, operating costs and carbon emissions are used as the optimization objectives of the VPP in this paper.

(1) Operating cost

The operation cost of the VPP includes the generation cost C G of CHP units, the operation and maintenance cost C M of various equipment, the cost C D R of controllable load, and the revenue I U G from the electricity trading on the public grid. min ⁡ F 1 = C G + C M + C D R − I U G . (18)

The generation cost of the CHP unit includes fuel cost and startup and shutdown cost, which are calculated as follows: C G = ∑ t = 1 T c C H 4 V C H 4 , t − V C H 4 , t t m + C D T , (19)

where c C H 4 represents the price of natural gas and C D T represents the start-up/shut-down costs.

The operation and maintenance costs include operation costs of wind power, photovoltaic, GPPCC, P2G, electric boilers, and power storage devices, which are calculated as follows: C M = ∑ t = 1 T c 1 g W P P , t + c 2 g P V , t + c 3 g G P P C C , t + c 4 g P 2 G , t + c 5 g e h , t + c 6 g e s , t d i s + g e s , t c h ) , (20)

where c 1 、 c 2 、 c 3 、 c 4 、 c 5 and c 6 represent the operating cost coefficients of wind power, photovoltaic, GPPCC, P2G, electric boiler, and power storage device, respectively.

The controllable load cost includes the response output cost and the standby output cost. The specific calculation is as follows: C D R = ∑ t = 1 24 ∑ k = 1 N I c I , k u Δ L k , t u + c I , k d Δ L k , t d + c R , k u R k , t u + c R , k d R k , t d , (21)

where c I , k u and c I , k d denote the cost coefficient of providing positive/negative response output for the user k, respectively; R k , t u and R k , t d are the positive/negative spare capacity that can be provided by the user k, respectively; and c R , k u and c R , k d denote the cost coefficient of providing positive/negative standby output for user k, respectively.

The revenue from electricity purchase and sale of public grid is calculated as follows: I U G = ∑ t T c U G , t g U G , t , (22)

where c U G , t denotes the electricity price of public power grid and g U G , t indicates electricity sold (purchased) to the public grid for VPP.

(2) Carbon emissions

Considering that China is still dominated by thermal power generation, the equivalent carbon emissions of purchased public grid electricity are also reckoned in the carbon emissions of CHP-VPP. The expression is written in the following form: min ⁡ F 2 = ∑ t = 1 T Q c , t s − η U G ⁡ min g U G , t , 0 , (23)

where η U G represents the carbon emission coefficient per unit of electricity.

The constraints of VPP conventional dispatching model mainly include electric/thermal power balance constraints, CHP unit output constraints, controllable load constraints, equipment operation constraints, and gas storage device constraints.

(1) Electric/thermal power balance constraints

The VPP proposed in this paper includes two kinds of energy flows, electric and thermal, and needs to meet both power/thermal balance constraints. g W P P , t + g P V , t + g G e , t + g e s , t d i s + Δ L I , t = L e , t + g e s , t c h + g G P P C C , t + g P 2 G , t + g e b , t + g U G , t h G , t + h e b , t = L h , t , (24)

(2) CHP unit output constraints

The CHP unit output constraints primarily include the upper and lower limit constraints of the unit thermal output, electrical output, and total output: 0 ≤ h G , i , t ≤ h G , i , ⁡ max s G , i , t ⁡ max g G , i , ⁡ min − η e h , i h G , i , t , α i h G , i , t + β i ≤ g G e , i , t ≤ s G , i , t g G , i , ⁡ max − η e h , i h G , i , t s G , i , t g G , i , ⁡ min ≤ g G , i , t ≤ s G , i , t g G , i , ⁡ max , (25)

where h G , i , ⁡ max represents the maximum value of thermal output; g G , i , ⁡ max and g G , i , ⁡ min are the max/min total output, respectively; α i represents the elastic coefficient of electric power and thermal power and can be considered as a constant; and β i represents a constant.

(3) Controllable load constraints

Controllable load constraints mainly include upper limit constraints u k , t u L k , t u − u k , t d Δ L k , t d + R k , t u ≤ Δ L k , ⁡ max u k , t u L k , t u − u k , t d Δ L k , t d − R k , t u ≥ Δ L k , ⁡ min , (26)

where Δ L k , ⁡ max and Δ L k , ⁡ min represent the maximum positive/negative response output that can be provided by user k, respectively.

(4) Equipment operation constraints

The equipment operating constraints consist primarily of upper and lower limit constraints and climb constraints for the GPPCC, P2G, and electric boilers. g k , ⁡ min ≤ g k , t ≤ g k , ⁡ max − Δ g k , d ≤ g k , t − g k , t − 1 ≤ Δ g k , u , (27)

where g k , ⁡ min and g k , ⁡ max are the min/max operating power of type equipment, respectively. Δ g k , u and Δ g k , d represent uphill/downhill climbing ability, respectively.

(5) Energy storage/gas device constraints

Constraints on energy or gas storage devices mainly include energy storage/gas capacity constraints, upper limit of charging and discharging rate constraints, charging and discharging state constraints, and equal energy storage/gas capacity limitations at the beginning and end of the cycle. Taking the gas storage devices as an example: 0 ≤ E t ≤ E max 0 ≤ V t i n   ≤ s t i n   V max i n   0 ≤ V t o u t   ≤ s t o u t   V max o u t   0 ≤ s t i n   + s t o u t   ≤ 1 E 0 = E 24 (28)

where E max represents the maximum storage capacity of the gas storage unit; s t i n and s t o u t represent the storage and venting states, respectively, and are 0–1 variables; and V max i n and V max o u t represent the maximum rates of gas storage and venting, respectively.

(6) System backup constraints

Because of the uncertainty of variable renewable energy, the conventional dispatching model of the VPP also requires consideration of system reserve constraints. This paper emphasizes the effect of load loss on the system when the actual generation power of wind power and PV is lower than the predicted power. The upper rotation reserve constraint is considered. R t u ≥ r W P P g W P P , t + r P V g P V , t R t u = R k , t u + R e s , t u , (29)

where r W P P and r P V represent upper rotational reserve coefficients of WPP and PV, respectively, and R t u is an upper rotation backup available for VPP. The reserve capacity R e s , t u provided for the power storage device. The operation mode of controllable load and power storage device is flexible, which can provide a certain reserve capacity for the VPP. However, the CHP unit has poor flexibility, so this paper does not consider it as a standby power supply. R e s , t u = ⁡ min g e s , ⁡ max − g e s , t , E e s , t − g e s , t Δ t Δ t , (30)

where g e s , ⁡ max represents the maximum input or output power of the power storage device and g e s , t represents the operating power of the ESD, which is equal to g e s , t d i s when it is positive and equal to g e s , t c h when it is negative.

Based on value at risk (VaR), CVaR takes into account the distribution of risk outside the confidence level, and can reflect the maximum possible loss in the full probability interval of the portfolio under a given level of confidence. Therefore, in this paper, the CVaR theory is utilized to quantify the risk of load loss in real-time dispatching of VPPs and is used as an optimization objective reflecting the operational risk of VPPs to cope with the uncertainty of variable renewable energy. The approximate formula of CVaR is as follows: F β x , α = α + 1 1 − β ∫ f x , y − α + ρ y d y , (31)

where x and y represent portfolio vectors and random vectors, respectively; f x , y represents the loss function; β represents the confidence; α represents the VaR value; ρ y is the joint probability density function of the random vector y ; and f x , y − α + represents max f x , y − α , 0 .

When the analytic formula ρ y is difficult to obtain, the integral term of Eq. 31 can be estimated by historical data or sample data obtained by Monte Carlo simulation. In this paper, the scenarios generated in Section 2.1 are used as samples, which are expressed as follows: F ^ β x , α = α + 1 N 1 − β ∑ n = 1 N f x , y n − α + , (32)

where y 1 , y 1 , ⋯ , y n represent N samples of y . The loss function values f n of each sample is arranged from large to small, and the β N first is the value of α .

Risk metrics are often related to the amount and duration of load loss, so by taking the penalty cost of VPP load loss as a loss function, and the specific calculation is as follows: C e n s = ∑ t = 1 T c e n s , t Δ g W P P , t + Δ g P V , t − R t u , (33)

where Δ g WPP , t and Δ g P V , t indicate deviations from actual wind and PV generation, respectively, and c e n s , t represents the penalty cost coefficient of load loss.

A multi-objective random dispatching optimization model for the VPP is as follows: min ⁡ F 1 = C G + C M + C D R − I U G min ⁡ F 2 = ∑ t = 1 T Q c , t s − η U G ⁡ min g U G , t , 0 min ⁡ F 3 , β = α + 1 N 1 − β ∑ k = 1 N C r i s k G , g k − α + . s . t . E q u a t i o n 22 − 26 (34)

It can be seen from Eq. 9 that the calculation process of Q c , t needs to be linearized by multiplying λ c , t and Q G , c , t . First, λ c , t will be discretized into 100 linear combinations of 0–1 variables. Since the value of λ c , t is between 0 and 1, this operation is equivalent to limiting the precision of λ c , t to 0.01. The details are as follows: λ c , t = 0.01 ∑ i = 1 100 λ i , t , (35)

where λ i , t represents the 0–1 variable. The results showed that Q c , t = 0.01 ∑ i = 1 100 λ i , t Q G , c , t . (36)

Then, by making Q i , t = λ i , t Q G , c , t , and adding the appropriate constraints, the goal of linearizing Q c , t calculation process is achieved. The details are as follows: Q c , t = 0.01 ∑ i = 1 100 Q i , t 0 ≤ Q i , t ≤ Q G , c , t Q i , t ≤ M λ i , t Q i , t ≥ Q G , c , t − M 1 − λ i , t , (37)

where M represents a large enough number. Similarly, the formula for multiplying other binary variable and continuous variable can be linearized.

Since the three objective functions in this paper have different orders of magnitude, the method based on fuzzy satisfaction is used for dimensioning the objective function (Gong et al., 2011). The fuzzy satisfaction theory can reflect the satisfaction degree of the objective function compared with the single-objective optimization, and its principle is to use the membership function of the fuzzy theory to quantify the solution of the objective function. First, each objective function is taken as the optimization object, the single-objective model is solved, and the values of other objective functions are calculated. See Table 1 for details. * denotes that the objective function is used as the optimization object.

The optimal values of each objective function can be obtained from Table 1, namely, F 1 min , F 2 min , and F 3 min . Then, the maximum value F 1 max , F 2 max , and F 3 max , is determined and can be scaled appropriately according to the preferences of the decision maker and the situation on the ground. F 1 min ≤ F 1 max ≤ F 1 2 , F 1 3 F 2 min ≤ F 2 max ≤ F 2 1 , F 2 3 . F 3 min ≤ F 3 max ≤ F 3 1 , F 3 2 (38)

Finally, the objective functions are all optimized in the direction of minimization, and each objective function uses ascending semi-linear membership functions as membership functions. The details are as follows: π i F i = 0 , F i ≤ F i min F i − F i min F i max − F i min , F i min < F i < F i max , 1 , F i ≥ F i max (39)

where π F i represents the membership function of objective function F i .

The entropy weight method is the most widely used method for solving VPP multi-objective problems. However, the entropy weight method mainly empowers through the degree of dispersion of each objective, ignoring the horizontal influence generated by the correlation between the objectives. CRITIC is an objective weighting method that considers the impact of index correlation. The principle is to determine the weight according to the contrast strength of the evaluation index and the correlation between the indexes, which can reduce the influence of the correlation between the indexes on the final weight and make the results more objective and reasonable. The general process of the CRITIC method is as follows:

(1) First, suppose there are m plans and n goals, respectively. Taking the solutions of F 1 , F 2 , and F 3 as objectives are taken as three CRITIC weighted schemes, and the following evaluation matrix is obtained. X = x 11 x 12 ⋯ x 1 m x 21 x 22 ⋯ x 2 m ⋮ ⋮ ⋮ x n 1 x n 2 ⋯ x n m , (40)

where x i j denotes the dimensioned value of the first j target of the first i scheme.

(2) Then, the standard deviation and correlation coefficient were calculated for each target, as follows: σ i = 1 m ∑ j = 1 m x i j − x i ¯ 2 ρ i k = cov X i , X k / σ i σ k , (41)

where σ i is the standard deviation of the target i ; ρ i k indicates the correlation coefficient between target i and target k; and cov X i , X k is the covariance of lines i and k .

(3) Calculating the amount of information contained in each goal and acquiring the weight of each goal, as follows: G i = σ i ∑ k = 1 n 1 − ρ i k u i = G i ∑ k = 1 n G k , (42)

where G i represents the information amount of the target and ∑ k = 1 n 1 − ρ i k represents the quantitative indicator of the conflict between the first goal i and other goals.

Finally, the combined objective function is as follows: F = ∑ i = 1 3 u i π i F i . (43)

For the purpose of this study, an industrial park in Lankao County, Henan Province is selected as the research object. The VPP of the park has two 0.8 MW CHP units, the total capacity of wind and PV is 1.2 MW and 0.4 MW, and the energy storage capacity is 0.1 MW. The maximum response outputs for the electric boiler capacity and controllable loads are 0.15 MW and 0.03 MW, respectively. The maximum operating power of carbon capture device is 0.1 MW, and the maximum operating power of electrolytic cell and methane reactor is 0.3 MW and 0.15 MW, respectively. In the conventional dispatching model, the spinning reserve coefficients of WPP, PV, and load are 0.25, 0.15, and 0.1, respectively. In the uncertain dispatching model, the penalty cost coefficient of load loss is 800 yuan/MW, and the confidence level of the CVaR value is 0.8. Figure 4 shows the wind power, photovoltaic output and electrothermal load predicted by the VPP dispatching center in day ahead. Figures 5, 6 show the actual output scenarios and the reduced typical scenarios of wind charge photovoltaic generated in this paper, respectively. In a typical output scenario for wind and photovoltaic power generation, there is a certain correlation between the output values, while the output values at each time also retain a certain degree of randomness, which is more in line with the actual output of wind power and photovoltaic.

This paper proposes a carbon recycling module considering the carbon capture device and power-to-gas device, and creatively decouples the generation and utilization process of CO₂ through carbon storage and hydrogen storage devices, while realizing the time shift of surplus renewable energy power. In addition, a risk quantification method based on CVaR theory is proposed. For the sake of verifying, the conclusiveness of the method propounded in this study, the following four scenarios are set up for simulation and analysis.

Scenario 1: Basic scenario. This scenario does not include carbon recycling module and the risk quantification method, but the conventional system backup constraint is applied to deal with the uncertainty of new energy.

Scenario 2: Carbon recycling scenario. This scenario introduces the carbon recycling module and does not adopt the risk quantification method in this paper.

Scenario 3: Risk quantification scenario. This scenario adopts the risk quantification method in this paper, without introducing the carbon recycling module.

Scenario 4: Comprehensive scenario. This scenario introduces the carbon recycling module and adopts the risk quantification method.

According to the multi-objective weighting method in Section 3, the weights of the objective functions of minimum operation cost, minimum carbon emissions, and minimum operation risk in Scenario 3 and Scenario 4 are 0.26, 0.3, and 0.44, respectively. Since Scenario 1 and Scenario 2 do not use the risk quantification method, and only include the minimum operating cost and the minimum carbon emissions, using the entropy weight method to calculate the weight, which are 0.59 and 0.44, respectively. Table 2 shows the optimization results of each scenario.

According to Table 2, the operation cost, carbon emissions, and operation risk of Scenario 1 are 10,606.46¥, 8,594.14 kg, and 7.6¥, respectively. Compared with Scenario 1, Scenario 2 utilize the surplus wind power generation to achieve the recycling of CO₂ owing to the introduction of carbon recycling module, reduce the fuel cost of CHP units, and reduce the operating cost and carbon emissions by 23.95¥ and 280.6 kg, respectively. Scenario 3 measures the risk level in the real-time operation of the VPP by adopting the risk quantification method, and develops a dispatching scheme with risk and economy, which reduces the operation cost and carbon emissions by 456.38¥ and 153.75 kg, respectively, while the operation risk only increases by 81.02¥. Based on Scenario 2 and Scenario 3, the operating cost and carbon emissions of Scenario 4 are further reduced by 689.95¥, 257.52¥, 245.47 kg, and 372.32 kg. Figure 7 shows the operating power of each unit in the VPP under each scenario.

According to Figure 7, the CHP unit is limited by the thermoelectric ratio and the minimum output, and maintain high output level all the time. The electric boiler uses wind power to supply heat for the system in periods 1–8 and 22–24, and conducts thermoelectric decoupling. The controllable load and power storage device mainly maintain the power balance of the VPP, providing access space for wind power and photovoltaic, and reserve capacity for the VPP. During periods 1–8, 11–16, and 23–24, the output of WPP and PV is high, and VPP sells surplus renewable energy power on the main network. On this basis, this section will further analyze the carbon emission reduction capability of the proposed carbon recycling module and the uncertainty response capability of the risk quantification method. Compared with Scenario 1, the operating power of the CHP unit in Scenario 2 increases slightly, the operating power of electric boilers is higher, and more electric energy is sold in the electricity market. Scenario 2 introduces the carbon recycling module, which requires more power consumption. The consumption of wind power and photovoltaic is greatly increased, increasing of downlink calls of controllable load, to improve the uplink spare space.

Compared with Scenario 1, the operating power of the CHP unit in Scenario 3 is slightly lower, and more electric energy is sold in the power market because Scenario 3 adopts the risk quantification method, and chooses to absorb more scenic calls to improve the economy of VPP, while taking certain risks. Therefore, the number of calls of controllable loads in Scenario 3 is less, to save the backup cost of VPP.