The combined heat and power virtual power plant mainly includes a distributed electrical/thermal output module and a carbon recycling module. Among them, the distributed electrical/thermal output module includes distributed wind power, distributed photovoltaic, electric boiler, and controllable load. The carbon recycling module includes CHP unit, GPPCC, P2G, and carbon storage and hydrogen storage devices, which can recycle CO₂ generated by the CHP unit. Besides, GPPCC includes carbon capture and carbon storage. The VPP dispatching center will predict the available energy output in advance, obtain the operating status of each unit, and formulate an electric heating cooperative dispatching plan for VPP. In addition, VPPs can interact with power grids to fill power supply gaps or sell surplus power. The energy or material flow relationships between the components of VPP is shown in Figure 1.

In this paper, VPP mathematical model including CHP, P2G, GPPCC and other components is established. In addition, a multi-objective optimization model considering operational cost, carbon emission and operational risk is constructed. By solving the optimization model, the optimization objectives of VPP such as reducing carbon emission, reducing operation cost and optimizing power generation plan can be achieved, so as to make the operation of power system more efficient, stable and reliable.

This paper describes the uncertain factors of renewable energy output through the generation of uncertainty scenarios. In order to model the probability distribution of output power of fan and photovoltaic, it is necessary to mine the information of historical data to directly model the uncertainty of output power. Currently, Latin hypercube sampling is the most common method for scene generation (Zhang et al., 2023; Ju et al., 2024), but this method ignores the correlation between the renewable energy output at different times. Therefore, in order to take into account the randomness and correlation of renewable energy output at all times, this paper proposes a scenario generation method considering the temporal correlation of wind power and PV output. The steps of this method are as follows:

(1) First, the covariance matrix σ 24 × 24 of the full-cycle wind prediction error is constructed as shown in Equation 15:

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

 (2) Z 1 × 24 ∼ N 0 , σ 24 × 24 multivariate normal distribution of the full-cycle wind prediction error is constructed, and m v n r n d function in Matlab is called to generate N samples randomly.

In order to reduce the amount of computation, this paper uses k-means clustering to reduce scenes to n typical scenes.

The objective functions of VPP conventional dispatching optimization model include minimum operating cost and minimum carbon emission.

(1) Minimum operating cost

The operating cost includes the power generation cost of CHP unit C G , the maintenance and operation cost of various equipment C M , the controllable load cost C D R and the income from the buying and selling of electricity in the public grid I U G . The calculation formula is shown in Equation 16. min ⁡ F 1 = C G + C M + C D R − I U G (16)

The power generation cost of CHP units includes fuel costs and start-up and shutdown costs, which are calculated as Equation 17: C G = ∑ t = 1 T c C H 4 V C H 4 , t − V C H 4 , t m + ∑ t = 1 T c D T u c , t − u c , t − 1 (17)

Where, c C H 4 is the expense of natural gas; c D T is the start-stop cost; u c , t is the start-stop variables of CHP at time t .

Operation and maintenance costs include the operating costs of wind power, photovoltaic, GPPCC, P2G, and electric boilers, which are calculated as Equation 18: 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 H 2 , t + g C H 4 , t m + c 5 g e b , t (18)

Where, c 1 , c 2 , c 3 , c 4 and c 5 are the operating cost coefficients of wind power, photovoltaic, GPPCC, P2G, and electric boilers respectively.

Controllable load cost includes response output cost and standby output cost (Hao et al., 2023), which are calculated as Equation 19: 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 (19)

Where, the cost factor of c I , k u and c I , k d providing positive/negative response power to the k user; R k , t u and R k , t d represent the positive/negative spare capacity that can be provided by the k th user; c R , k u and c R , k d provide positive/negative spare power cost factors for the k the user.

The income from the buying and selling of electricity for the public grid is calculated as Equation 20: I U G = ∑ t = 1 T c U G , t g U G , t (20)

Where, c U G , t is the electricity cost of the public grid; g U G , t is the amount of electricity sold (purchased) via the VPP to public grids.

(2) Minimum carbon footprint

Considering that China is still dominated by thermal power generation, this paper will also include the equivalent carbon emissions of electricity purchased in the public grid into the carbon emissions of VPP. The calculation formula is shown in Equation 21. min ⁡ F 2 = ∑ t = 1 T Q c , t s − η U G ⁡ min g U G , t , 0 (21)

Where, η U G is the carbon emission coefficient ever unit of electricity.

VPPs routine dispatching optimization model includes the following constraints:

(1) Electrical/thermal power balance constraint

In order to achieve the electricity/heat supply and demand balance of the VPP in each period, Equation 22 is established: g W P P , t + g P V , t + g G e , t + Δ L I , t = L e , t + g G P P C C , t + g H 2 , t + g C H 4 , t m + g c o , t + g e b , t + g U G , t h G , t + h e b , t = L h , t (22)

(2) CHP unit output constraint

In order to limit the thermal and electrical output of CHP unit from exceeding its output range, Equation 23 is established: 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 ∑ w = t t + T S − 1 1 − u c , w ≥ T S u c , t − 1 − u c , t ∑ w = t t + T O − 1 u c , w ≥ T O u c , t − u c , t − 1 (23)

Where, h G , i , ⁡ max is the utmost value of thermal output; g G , i , ⁡ max and g G , i , ⁡ min are the greatest and least values of the total output; α i is the elastic coefficient of electric power and thermal power; β i is a constant; T S and T O are the minimum off/on time, respectively. Please refer to reference (Hao et al., 2023) for the specific climbing constraints of the unit.

(3) Controllable load constraint

In order to realize the reasonable transfer of controllable load, the transferable range is set. Therefore, Equation 24 is established: μ k , t u Δ L k , t u − μ k , t d Δ L k , t d + R k , t u ≤ Δ L k , ⁡ max μ k , t u Δ L k , t u − μ k , t d Δ L k , t d − R k , t d ≤ Δ L k , ⁡ max (24)

Where, Δ L k , ⁡ max and Δ L k , ⁡ min are the maximum positive/negative response forces that can be provided by the k th user.

(4) Equipment operation constraint

Similarly, GPPCC, P2G and other units also need to set their output ranges. Therefore, Equation 25 is established: g k , ⁡ min ≤ g k , t ≤ g k , ⁡ max − Δ g k , d ≤ g k , t − g k , t − 1 ≤ Δ g k , u (25)

Where, g k , ⁡ min and g k , ⁡ max are the least and greatest operating power of type k equipment; Δ g k , u and Δ g k , d represent up/down climbing capability.

(5) Gas storage device constraint

The decision maker needs to consider the capacity limit and output range of the gas storage device, so Equation 26 is established: 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 (26)

Where, E max is the utmost gas storage volume of the gas storage devices; V max i n and V max o u t are the maximum gas storage and venting rates; s t i n and s t o u t are the status of gas storage and venting of gas storage devices, respectively.

The conventional VPP scheduling model also includes system backup constraints, for details, please refer to Ref. (Zou et al., 2023).

On the basis of value at risk (VaR), CVaR considers the risk distribution beyond the confidence level, and can reflect the maximum possible loss within the full probability interval of the portfolio under the given confidence degree. Therefore, this paper adopts CVaR theory to quantify the loss of load risk in real-time VPP scheduling, and takes it as an optimization target to reflect the operational risk of VPP, so as to cope with the incertitude of wind power generations. The approximate calculation formula of CVaR is as shown in Equation 27 (Ju et al., 2022): F ^ β x , α = α + 1 N 1 − β ∑ n = 1 N f x , y n − α + (27)

Where, f x , y is the loss function; x is the portfolio vector, y n is the uncertainty scenario generated in Section 2.1; α and β represent VaR values and confidence levels; f x , y − α + is equivalent to max f x , y − α , 0 .

The measurement index of risk is often related to load loss and load loss duration (He et al., 2023), so this paper takes the loss penalty cost C e n s of VPP as the loss function, and the specific calculation is as shown in Equation 28: C e n s = ∑ t = 1 T c e n s , t Δ g W P P , t + Δ g P V , t − R t u (28)

Where, Δ g W P P , t and Δ g P V , t are the deviation of the actual power generation of the scenery; c e n s , t is the penalty cost coefficient of loss of load; R t u is the uplink standby capacity of VPP, which is mainly provided by the extraction steam unit, and the insufficient part is provided by the controllable load.

VPP multi-objective random dispatching optimization model is as shown in Equation 29: 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 − α (29)

s.t. Equation 22 – 26This section considers the uncertainty of new energy output, combined with the conventional VPP scheduling model, and then constructs the VPP multi-objective random scheduling optimization model. According to this idea, the decision maker can make the optimal VPP unit scheduling scheme.

The three objective functions in this paper have different orders of magnitude, so the reduced semi-gradient membership function is used to de-dimensionalize. For the specific method, please refer to reference (Xuan et al., 2021). Membership function is as shown in Equation 30: π F i = 0 , F i ≥ F i max F i max − F i F i max − F i min , F i min ≤ F i ≤ F i max 1 , F i ≤ F i min (30)

Where, F i is the value of the i TH objective function; F i min and F i max are the least and greatest values of the objective function.

There are two kinds of weighting methods: subjective weighting and objective weighting. The results of subjective weighting depend heavily on the subjective knowledge of experts. The results of objective weighting may not necessarily represent the actual importance of the indicators (Song et al., 2020). Neither subjective weighting method nor objective weighting method can perfectly reflect the importance of each objective function. Therefore, this paper chooses the analytic hierarchy process (AHP) as the subjective weighting method and the entropy weighting method as the objective weighting method. An integrated subjective and objective weighting method is proposed to assign weights to each optimization objective. Specific calculations are as shown in Equation 31: w i = u i r v i 1 − r ∑ i = 1 3 u i r v i 1 − r , i = 1 , 2 , 3 (31)

Where, w i is the weighting obtained by assigning weight to subjective and objective integration; v i and u i are the weights obtained by analytic hierarchy procedure and entropy weight method; r is the preference coefficient of decision makers for subjective and objective factors, with the value between 0 and 1.

After the objective function is de-dimensional and weighted by subjective and objective integration, Equation 29 can be converted into the form of Equation 32, and the result of VPP multi-objective scheduling optimization can be obtained by solving it. In addition, the first item of Equation 8 needs to be linearized. min ⁡ f = ∑ i = 1 3 w i π F i s . t Equation 22 – 26 . (32)

For the purpose of this paper, the VPP of a certain place in China is selected as the simulation object, VPP has two 0.8 MW CHP units with a combined capacity of 1 MW for wind and 0.4 MW for PV, a capacity of 0.15 MW for electric boilers, and a maximum response output of 0.03 MW for controllable loads. The penalty cost coefficient of loss of load is 800 yuan/MW, the confidence is 0.8, and the power upper limit of VPP interacting with the grid is set to 0.1 MW. Figure 2 shows the predicted next-day wind power and electric heating load. It can be seen that the scene generation method takes into account the randomness and correlation of the output of the scenery at every moment, and is more in line with the actual output of the scenery. In this paper, examples are simulated on Matlab R2016a. The time spent in generating and reducing uncertain scenes is about 5s, and the results can be obtained within 15s for model solving.

In order to analyze the carbon recycling capability of GPPCC and P2G and the effectiveness of the uncertainty coping method, the following four scenarios are set for simulation analysis.

Scenario 1: Basic scenario. GPPCC and P2G were not introduced, and the uncertainty coping method in this paper was not adopted.

Scenario 2: Carbon recycling scenario. GPPCC and P2G were introduced, but uncertainty coping methods were not adopted.

Scenario 3: Risk avoidance scenario. The uncertainty coping method was adopted, but GPPCC and P2G were not introduced.

Scenario 4: Integrated scenario. Both GPPCC and P2G are introduced, and uncertainty coping methods are adopted.

Firstly, the entropy weight method and analytic hierarchy procedure were used to calculate the weight of the objective function under each scenario. Then, the sensitivity analysis of decision-maker’s subjective and objective factor preference coefficient is carried out, and the weight of subjective and objective integration was calculated according to Equation 31, as shown in Table 2.

It can be seen that when r gradually increases, the weight of running cost F ₁ gradually decreases. In addition, the weight of carbon emission F ₂ and risk cost F ₃ are gradually increasing. This is because when r is small, the subjective weight has more influence, and the subjective weighting method pays more attention to the impact of operating costs on the system. It is worth pointing out that in Scenario 2-4, F ₁ is the most weighted optimization target regardless of how r changes. This also reflects the importance of running costs to the entire system.

Table 3 shows the optimization results of scenarios when r = 0.1, and Figure 3 shows the difference between the optimization results of scenarios when r ∈ [0.1,0.9] and when r = 0.1. In the legend of Figure 3C, m _carbon refers to the amount of carbon recycling. Combined with Table 3 and Figure 3, a comparative analysis of each scenario is carried out: In Scenario 2, compared with Scenario 1, the addition of carbon recycling devices increases the overall cost slightly, but the carbon emissions decrease significantly, by about 7%. In Scenario 3, compared with Scenario 1, the uncertainty coping method proposed in this paper reduces the overall cost by about 3.5%. In addition, carbon emissions are also reduced, by about 3%. In comparison with Scenario 2 and Scenario 3, the use of the comprehensive method makes the operating cost and carbon emission of Scenario 4 the lowest value among the four scenarios, and the carbon recovery amount of Scenario 4 is higher than that of Scenario 2. To sum up, adding carbon cycle device can greatly improve the environmental benefit; The use of CVaR method can greatly improve the economic benefit. Using a comprehensive approach is more effective than using a single approach.

In addition, the relationship between the value of the objective function and r can be analyzed from Figure 3. With the increase of r, the weight of F ₁ decreases, while the weight of F ₂ and F ₃ increases. Therefore, the operation cost of each scenario gradually increases, and the carbon emission and risk cost overall show a downward trend. In addition, for Scenario 2 and Scenario 4 where carbon recovery exists, the amount of carbon recovery increases significantly, with an increase ratio of 3.08% and 1.61% respectively. For Scenario 3 and 4 with risk cost, as r increases from 0.1 to 0.9, F ₃ decreases by 99.7% and 69.53% in Scenario 3 and 4, respectively.

According to the results of the sensitivity analysis in Figure 3, when the value of r is centered, the importance distribution of each objective in VPP is more balanced: VPP neither attaches too much importance to operational risks nor ignores carbon recycling. Therefore, in order to compare the results of different scenarios more conveniently, the results when the subjective and objective factor preference coefficient r is 0.5 are selected for further analysis.