Integrated Optimization for Biofuel Management Associated with a Biomass-Penetrated Heating System under Multiple and Compound Uncertainties

Abstract

The biofuel management of a biofuel-penetrated district heating system is complicated due to its association with multiple and polymorphic uncertainties. To handle uncertainties and system dynamic complexities, an inexact two-stage compound-stochastic mixed-integer programming technique is proposed, innovatively based on the integration of different uncertain optimization approaches. The proposed technique can not only address the inexact recourse problems sourced from multiple and compound uncertainties existing in the pre-regulated biofuel supply–demand match mode, but can also quantitatively analyze the conflicts between the economic target that minimizes the system cost and the risk preference that maximizes the heating service satisfaction. The developed model is applied to a real-world biofuel management case study of a district heating system to obtain the optimal biofuel management schemes subject to supply–demand, policy requirement constraints, and the financial minimization objective. The results indicate that biofuel allocation and expansion schemes are sensitive to the multiple and compound uncertainty inputs, and the corresponding biofuel-deficit change trends of three heat sources are obviously distinct with the system’s condition, varying due to the complicated interactions of the system’s components. Beyond that, a potential trade-off relationship between the heating cost and the constraint-violation risk can be obtained by observing system responses with thermalization coefficient varying.

1. Introduction

Against the background of worldwide carbon neutrality, biomass energy has received growing attention and been gradually penetrating into heating systems in the recent years in Northern China [1,2]. A typical biofuel-penetrated district heating system (BDHS) involves a main coal-fired heat source (MCHS) and several biomass-based peak-shaving heat sources (BPHSs) to meet the overall heating demand. The biofuel management (BM) of a BDHS is complicated and associated with different activities, including energy resource supply and distribution, market enthusiasm, demand variation, policy guidance, and environmental impact. As a result, controversial and conflict-laden issues related to the biofuel/heat supply–demand, pollutant emission-limitation, and cost expectation may exist during BM [3]. Therefore, the effective BM of a BDHS is important, because it can not only improve the heating comfort level and system economy, but can also reduce and sequester carbon dioxide and offset its emissions from fossil fuel sources effectively. Beyond the complicated interaction, the BM complexity of a BDHS can be further intensified by various uncertainties (e.g., supply intermittency, volatility, and demand fluctuation) and their various interaction [4]. If the uncertainties are ignored or simplified, the modelling condition may not reflect the reality, and the generated decisions may lead to mismatching or even serious accidents.

During the past decades, there was little research focusing on the uncertainty within the fuel management of a heating system, while a number of inexact optimization techniques were proposed to tackle uncertainties and complexities in the field of energy management and analysis (EMA) [5,6]. Among them, the interval linear programming (ILP) proposed by Huang et al. has been considered to be an effective approach to handling the inexact system parameters and model inputs, which cannot be expressed with precision but can only be described as discrete intervals (i.e., interval numbers with deterministic lower and upper bounds) [7]. For example, Huang et al. applied the ILP model to a hypothetical problem of solid waste management for the first time, indicating that reasonable solutions can be generated for both the upper and lower bounds of the objective function cases [8]. Guo et al. developed an ILP-based method for a regional energy-management system, and the solutions expressed with interval numbers can be obtained to provide energy decision alternatives [9]. Cai et al. addressed the uncertainties in the energy planning model through the ILP approach and improved the robustness of an interactive decision support system for the Region of Waterloo in Canada [10]. From previous application studies, it can be seen that, in addition to the advantage of a lower data requirement, the effective two-step interactive solution algorithm of ILP can generate the interval solutions with a low uncertainty degree, leading to the widespread acceptance of ILP, not only in EMA [11,12], but also in other resource management fields [13,14]. Thus, the ILP technique is suitable for coping with the economic and biomass-pellet-fuel quality parameters during BM. However, the ILP is incapable of dealing with the probabilistic variable (e.g., biofuel-availability scenarios) that are associated with different levels of economic penalties when the planned biofuel-allocation targets are violated.

Considering a real BDHS with several BPHSs, the local available biofuels cannot be quantified precisely because of market enthusiasm, biomass-pellet-fuel production capacity and policy orientation, and the biofuel availability can be represented with an interval variable corresponding to different stochastic scenarios such as abundant, medium, and scarce available levels (i.e., a type of compound uncertainty) [15]. Beyond the compound uncertainty existing on the supply side, another type of compound uncertainty lying on the heat-demand side could also aggravate the decision-making difficulty. The heat provisions undertaken by heat sources may be simultaneously affected by the continuous-type random residents’ behaviors (e.g., the random ventilation and the selective use of district heating service) and discrete-type random meteorological condition during heating periods. In this situation, the compound-random parameters (i.e., a birandom variable), which is a measurable mapping from a probability space to a collection of random variables, would be generated for reflecting reality accurately and make it difficult for decision makers to provide a satisfactory heating service to the heat consumers [16]. For coping with the complex decision problem induced by the multi-stochastic uncertainties, a class of effective recourse optimization techniques known as two-stage stochastic programming (TSP), and its derivatives, were proposed previously and have been successfully applied to the energy management studies [17,18]. Through merging different randomness into both decision stages and running the recourse mechanism, the TSP-based techniques can effectively lower the risk of decision-making misplay when using the traditional deterministic optimization methods.

Specifically, Lin and Huang constructed an inexact stochastic two-stage energy-planning model for Beijing Municipality to manage the energy systems and greenhouse-gas emissions under uncertainty, suggesting that the proposed recourse model was applicable in reflecting the complexities of multi-uncertainty, dynamic, and interactive municipal energy-management systems [19]. Zhou et al. applied the TSP technique in the optimal design of a distributed energy system in a hotel. The proposed model presented the advantage of the recourse mechanism during the decision process, compared with a deterministic optimization model [20]. Ji et al. integrated the TSP into the day-ahead dispatch model for an electricity system management with wind power under uncertainty. The integrated model was capable of providing different optimal dispatch strategies corresponding to different scenarios for thermal power units and wind turbines [21]. Fu et al. improved the general TSP by introducing the discrete random variable into the first decision stage and applied the improved model into a real-world district heating system, obtaining hybrid fuel management schemes under different scenario combinations [22]. Nevertheless, the previous studies reveal that the TSP method and its derivatives can deal with discrete-type stochastic uncertainties in both decision stages effectively while they can hardly handle the compound-random uncertainty embedded within the constrained right-hand-side parameters or address the quantitative risk of violating uncertain constraints.

Fortunately, chance-constrained programming (CCP) is an alternative capable of tackling continuous-type random parameters in the constraints. CCP requires that all of the constraints should be satisfied in a proportion of cases under given probability levels, which could be linked to a system risk or the constraint satisfaction degree in the resource and energy management [23]. Previously, a plethora of CCP methods was proposed to deal with the risk-oriented power planning and resource management issues, whereas the applications of this method in the field of fuel management of a BDHS were scarce [24,25]. Through combining the recourse model (i.e., the TSP-based technique) with CCP, the compound-random uncertainties can be effectively decoupled and directly communicated into the solving optimization process, such that a wealth of pertinent information could be effectively merged into the decision process and multiple decision alternatives could be generated through the risk-based interpretation for the solutions. In addition to the complexity of the mentioned uncertainties with different formats, the dynamic complexities, such as timing, sizing, and siting decisions in terms of the heating-capacity expansion schemes of heat sources, should also be considered for the high-quality fuel management. From a long-term planning point of view, the heating capacities of existing heat sources will have cumulative or time-sharing limits, while the residents’ heat demands keep growing, owing to the population increase and society development. This tendency may bring about the insufficient heating capacities of the existing heat sources to cover the overall increasing heat demand, indirectly impacting fuel utilization and adjustment. Therefore, the heating capacity expansion is also a crucial issue in BM for a BDHS, where the decision should contain whether a particular heat source development or expansion option needs to be undertaken. Mixed-integer linear programming (MILP) is a beneficial tool for this purpose by using integer variables to indicate whether an expansion action is conducted [26]. MILP has been widely used in the research fields of power, industry, logistics, and transportation, which could be taken as a helpful reference for BM [27,28].

Consequently, the main objective of this study is to propose an inexact two-stage compound-stochastic mixed-integer programming technique for facilitating the real BM of a BDHS under multiple and compound uncertainties. The novelty and contributions of this research paper are outlined in the following aspects. (1) The developed technique incorporates CCP, TSP-derivative, and MILP optimization techniques into a general ILP model framework innovatively, and will be applied to a practical BDHS case in northeastern China for the first time to support BM. (2) The generated scheme results will be beneficial for (a) coping with multiple and compound uncertainties in the formats of interval numbers, probability distributions, and interval-stochastic and compound stochastic variables lying in the BM of a BDHS; (b) identifying the biofuel allocation patterns of different heat sources under different system conditions and various scenario combinations; (c) facilitating the dynamic analysis of heating-capacity-expansion decisions; and (d) addressing the conflicts between economic objectives and system risk levels during the BM of a BSHS. This paper is structured as follows: Section 1 provides an overview of existing work conducted in the area of the uncertain optimization techniques for energy and resource management, the object of our study and the main contributions; Section 2 introduced the devised methodology, including the quantitative technique of the compound-stochastic heat provisions and integrated biofuel-management model technique under multiple and compound uncertainties; Section 3 and Section 4, respectively, describe the real case study utilized for demonstrating the practicality and validity of the proposed model and states the model results and discussion, which are followed by Section 5, which outlines the conclusion and potential challenges in the future work.

2. Methodology

Consider a problem wherein a decision maker of the BM in a BDHS is in charge of optimizing the biofuel allocation among several BPHSs and identifying their capacity-expansion schemes, with the objective of minimizing the total system cost. Various uncertain factors lying in the BDHS, such as biofuel availability, residential heat demands, and the associated economic implications, are linked to the BM process and presents multiple and even compound features. To tackle these uncertainties and system dynamic complexities, a novel optimization approach for BM of a BDHS is proposed in this paper. The traditional energy and resource-management recourse model under uncertainty is initially introduced as the basis of the proposed model; and then the quantitative technique of the compound stochastic heat provisions is provided to compute the model-constraint boundaries; last, an inexact two-stage compound-stochastic mixed integer programming and its solution method are developed via the integration of other different uncertain optimization models into the traditional recourse model.

2.1. Inexact Two-Stage Dual-Stochastic Programming

The traditional interval two-stage stochastic programming (ITSP) proposed by Huang and Loucks is capable of tackling the interval and stochastic uncertainties via the combination of ILP and TSP, and the recourse–decision can be obtained in terms of interval values after the random event have occurred in the second stage [29,30]. Nevertheless, such decision complexity of the BM in a BDHS could be further exacerbated by the other uncertainty existing in the first stage, which cannot be handled by ITSP. For instance, the annual fluctuation of the local meteorological condition imposes stochastic characteristics on the users’ heat demands, further influencing the predetermined biofuel allocation plan in the first decision stage. According to Fu et al., one potential approach for better addressing additional randomness is to predefine the various first-stage decisions corresponding to the foreseeable random events, and the corrective decisions associated with economic penalties can be obtained subsequently under different stochastic scenario combinations after another random event has taken place in the second stage [22]. Based on this technical route, an inexact two-stage dual-stochastic programming (ITDSP) method can be developed, which is capable not only of handling the uncertainties presented as discrete intervals (e.g., economic and technical parameters), but also of reflecting a novel type of complex running mechanism of “recourse” caused by dual stochastic uncertainties in both decision stages. For the BM, a typical ITDSP model can be formulated as follows.

$$\min f^{\pm }={\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{p}_l^{hd}CB{F}_{nl}^{\pm }TB{F}_{nl}^{\pm }}}+{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{m=1}^M{p}_l^{hd}{p}_m^{ba}PB{D}_n^{\pm }DB{F}_{nml}^{\pm }}}}$$

(1)

subject to:

$${\displaystyle \sum _{n=1}^N(TB{F}_{nl}^{\pm }-DB{F}_{nml}^{\pm })}\le AVB{F}_m^{\pm },\forall m,l$$

(2)

$$H{D}_{nl}^{\pm }\le \sigma \cdot TB{F}_{nl}^{\pm },\forall n,l$$

(3)

$$TB{F}_{nl}{}_{,\max }\ge TB{F}_{nl}^{\pm }\ge DB{F}_{nml}^{\pm }\ge 0,\forall n,m,l,$$

(4)

where the mathematical symbol “±” in superscript denotes that the corresponding parameters are interval parameters/variables; the subscripts n, m and l stand for different heat sources, biofuel-availability levels, and different residents’ heat-demand levels, respectively, n = 1, 2, …, N, m = 1, 2, …, M and l = 1, 2, …, L; f means the expectation of the biofuel management cost; TBF means the projected biofuel amount supplied to the heating system (i.e., the first-stage decision variable); CBF stands for the unit price of the purchased biomass-pellet fuel; DBF is the biofuel deficit amount owing to the mismatch between the projected biofuel supply and the real requirement (i.e., the second-stage decision variable); PBD is the economic penalty (including the additional transport and management cost) because of the biofuel deficit; AVBF represents the available biofuel amount from the local supply channels; HD means the space-heating demand of the residents; σ is the conversion coefficient; p^hd represents the probability of a certain heat-demand level; p^ba denotes the probability corresponding to a certain biofuel-availability level, which is independent of p^hd, and scenario combinations can thus be generated via sampling the underlying probability distributions p^hd and p^ba.

2.2. Heat Provisions Undertaken by Heat Sources

To conveniently map the heat-load profile and compute the design heat provisions, the non-dimensional comprehensive equations (NCEs) method can be used. The NCEs method has been widely adopted in Northern China owing to the demand of the relatively small-scale meteorological and building data when compared with those of the recent data-driven methods and the traditional simulation techniques [5]. The design heat provisions undertaken by BPHSs can be obtained and presented as follows:

$$Q_{nb}=24Q_n^{\prime }\left[(1-\beta )N_{\beta }-\frac{{\beta }_0{(N_{\beta }-5)}^{1+b}}{(1+b){(N_{zh}-5)}^b}\right],$$

(5)

where Q_nb denotes the design heat provision undertaken by BPHSs, $Q_n^{\prime }$ represents the total design space-heating load β denotes the thermal coefficient, β₀ denotes the temperature correction coefficient, N_β is the cumulative up-time of a BPHS, b stands for the exponential value of non-dimensional heating duration, and N_zh is the total heating duration.

In the practical BM for a BDHS, beyond the objective meteorological change during the heating season as mentioned before, the residents’ subjective random actions, including their ventilation actions and willingness to use heating services, will also impact biofuel consumption by changing the real heat provision indirectly. Thus, the compound-stochastic heat provisions can be formed via merging two layer of randomness and should be embedded within the modelling parameter (i.e., heat provisions or residents’ heating demands). Such compound-stochastic uncertainty can be formulated as shown in Figure 1.

To be more specific, in the first layer, the traditionally designed heat provisions under discrete-random meteorological conditions can be achieved by setting the temperature parameters of NCEs lower/higher than those under the “normal” condition to a certain extent and adjusting the heat durations to be longer/shorter than those of the “normal”. In the second layer, the randomness of residents’ actions is further taken into consideration, and, given the expectation and variance of heat provisions under different meteorological conditions (e.g., different “freezing-degree” levels), the continuous-random heat provisions (e.g., Gaussian distribution) corresponding to a discrete-probabilistic “freezing degree” level, can be generated. Therefore, the proposed compound-stochastic heat provisions can also be explained as a “random variable” (i.e., continuous-random heat provisions), taking random variable values (i.e., discrete-probabilistic meteorological conditions).

2.3. Inexact Two-Stage Compound-Stochastic Mixed-Integer Programming and Its Solution Method

In addition to the compound-stochastic heat provisions on the constraint right-hand sides mentioned in Section 2.2, considering the difference between the existing heating capacity and the future heating demand, the dynamic complexity of the heating-capacity expansion of a BDHS needs to be merged into the ITDSP model framework as well. As a consequence, an inexact two-stage compound-stochastic mixed-integer programming (ITCS-MIP) method can be developed innovatively by integrating the CCP and MIP techniques into the ITDSP model framework to tackle the BM problem related to the complicated component interaction and dynamic evolution, as well as multiple and compound uncertainties. The developed ITCS-MIP model can be formulated as follows.

$$\min f^{\pm }={\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{p}_l^{hd}CB{F}_{nl}^{\pm }TB{F}_{nl}^{\pm }}}+{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{m=1}^M{p}_l^{hd}{p}_m^{ba}PB{D}_n^{\pm }DB{F}_{nml}^{\pm }}}}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{n=1}^{N}X{D}_{nk}^{\pm }\cdot C{D}_k^{\pm }}}$$

(6)

subject to:

$${\displaystyle \sum _{n=1}^N(TB{F}_{nl}^{\pm }-DB{F}_{nml}^{\pm })}\le AVB{F}_m^{\pm },\forall m,l$$

(7)

$$\mathrm{Pr}\{H{D}_{nl}^t\le \sigma \cdot TB{F}_{nl}^{\pm }\}\ge 1-p,\forall n,l$$

(8)

$$TB{F}_{nl}{}_{,\max }\ge TB{F}_{nl}^{\pm }\ge DB{F}_{nml}^{\pm }\ge 0,\forall n,m,l$$

(9)

$$X{D}_{nk}^{\pm }=\left\{\begin{array}{l}1;ifheatingcapacityexpansionisundertaken\\ 0;ifotherwise\end{array}\right\}$$

(10)

$$D{E}_n^{\pm }\le {\displaystyle \sum _{k=1}^{K}X{D}_{nk}^{}}\cdot D{H}_k^{}+E_n^{\pm },\forall n,$$

(11)

where the subscript k stands for the different heating capacity expansion options, k = 1, 2, …, K; XD is a binary decision variable for identifying whether or not a heating capacity expansion action needs to be executed; CD represents the cost of a heating-capacity expansion choice; Pr{}denotes the probability of the events in {}; p is a set of predetermined constraint-violation (i.e., system risks) probability levels, and thus HD is a type of compound-stochastic parameter; DE denotes the heat load undertaken by a BPHS; DH means a heating capacity expansion alternative; E stands for the existing heating capacity of a BPHS.

Among the multiple and polymorphic uncertainties (e.g., compound-stochastic variables, interval-stochastic variables, and discrete intervals) existing in the developed ITCS-MIP model, the nonlinear chance constraints ought to be initially converted into the “crisp constraint” [31,32]. Thus, according to the equivalent transformation method [33,34], the constraint (i.e., constraint (8)) can be reformulated equivalently as follows.

$$H{D}_{nl}^{t,1-p}\le \sigma \cdot T_{nl}^{\pm },$$

(12)

where $H{D}_{nl}^{t,1-p}=F^{-1}(1-p)$, given the cumulative distribution function of F(p) and the probability (i.e., p) of violating constraint.

In such a manner, the compound-stochastic uncertainty existing in the model can be converted into the traditional discrete probability distribution and the ITCS-MIP model is transformed to a general ITDS-MIP model under a constraint-violation probability (i.e., system risk) level. After that, the interactive two-step solution algorithm can be used to cope with the interval parameters in the reformulated model [8,29]. In general, via the crisp conversion of chance constraints of CCP in conjunction with the interactive two-step solution algorithm of ILP, the developed ITCS-MIP model is able to cope with the discrete probabilistic, compound-stochastic, and interval uncertainties lying in the model parameters/coefficients. The detailed modelling process and solution procedure is summarized in the Figure 2.

3. Case Study

3.1. Biofuel Management Problem Statement of the Investigated BDHS

A biofuel-penetrated district heating system (BDHS) located in Jinpu New District, Dalian City of China, is used for demonstrating the practicality and validity of the ITCS-MIP model on BM. The main BM problem to be solved is how to allocate limited local biomass-pellet fuel during three consecutive heating seasons to different biomass-based peak-shaving heat sources (BPHSs) in this BDHS by minimizing the total system cost with considering heating capacity expansion and the biofuel-shortage penalty.

The BDHS has three BPHSs (fueled with solid pellet biofuels) marked as BPHS_1, BPHS_2, and BPHS_3 with the heating efficiencies of [0.845, 0.855], [0.835, 0.84], and [0.825, 0.83], respectively. They provide the corresponding peak-shaving heating service to Heping, Hongqi, and Shengli sub-districts. In the normal meteorological condition, the heating period is 130 days. However, owing to both the existing climatic change and the residents’ random actions, the compound-stochastic uncertainties would be embedded within the heat provisions undertaken by different BPHSs. In this context, the traditional NECs method applicable to the “normal” heating season will lead to serious deviation from the practical heat provisions in the “severe” or “mild” heating seasons.

To quantify the uncertain heat provisions under different “freezing-degree” levels, on one hand, the appearance probabilities (i.e., the first layer of the birandom uncertainty) of the “severe”, “normal”, and “mild” heating seasons are, respectively, assumed to be 0.25, 0.45, and 0.30, and the corresponding meteorological parameters for NECs modeling are shown in

Table 1. Modelling parameters for NCEs method [5].

“Freezing Degree” Level	Average Outdoor Temperature, °C	Design Outdoor Temperature, °C	Design Indoor Temperature, °C	Space-Heating Durations, day
“Severe” (phd = 0.125)	−2.4	−11.5	18	140
“Normal” (phd = 0.55)	−1.9	−11	18	130
“Mild” (phd = 0.325)	−1.4	−10.5	18	120

. On the other hand, considering residents’ irregular behaviors in the heating season, such as the random actions for ventilating and employing a district heating service, this type of uncertainty (i.e., the second layer of the birandom uncertainty) can be expressed as a continuous random variable (e.g., Gaussian distribution). In detail, the design heat provisions corresponding to the “normal” heating season can be used as the mean value of the practical heat provision distribution. After that, according to Figure 1, the practical heat provisions under different “freezing degree” levels and predetermined risk probability levels can be achieved with the aid of the given variance (i.e., 1000 GJ used herein) and the NECs method mentioned in Section 2.2. The obtained heat provisions under the compound-stochastic uncertainty are shown

Table A1. Heat provisions undertaken by different peak-shaving heating sources under α = 0.5, GJ.

Percentile	Heat Source	“Freezing Degree”	Period 1	Period 2	Period 3
1%	BPHS_1	Severe	153,898.6	169,521.1	186,705.9
		Normal	143,440.4	158,017.1	174,051.4
		Mild	132,948.9	146,476.4	161,356.6
	BPHS_2	Severe	64,627.19	71,322.54	78,687.43
		Normal	60,145.1	66,392.24	73,264.1
		Mild	55,648.72	61,446.22	67,823.48
	BPHS_3	Severe	19,991.5	22,223.28	24,678.24
		Normal	18,497.47	20,579.85	22,870.47
		Mild	16,998.67	18,931.17	21,056.93
5%	BPHS_1	Severe	154,580	170,202.5	187,387.3
		Normal	144,121.8	158,698.5	174,732.8
		Mild	133,630.3	147,157.8	162,038
	BPHS_2	Severe	65,308.69	72,004.04	79,368.93
		Normal	60,826.6	67,073.74	73,945.6
		Mild	56,330.22	62,127.72	68,504.98
	BPHS_3	Severe	20,673	22,904.78	25,359.74
		Normal	19,178.97	21,261.35	23,551.97
		Mild	17,680.17	19,612.67	21,738.43
10%	BPHS_1	Severe	154,943.3	170,565.8	187,750.6
		Normal	144,485.1	159,061.8	175,096.1
		Mild	133,993.6	147,521.1	162,401.3
	BPHS_2	Severe	65,671.99	72,367.34	79,732.23
		Normal	61,189.9	67,437.04	74,308.9
		Mild	56,693.52	62,491.02	68,868.28
	BPHS_3	Severe	21,036.3	23,268.08	25,723.04
		Normal	19,542.27	21,624.65	23,915.27
		Mild	18,043.47	19,975.97	22,101.73
90%	BPHS_1	Severe	157,506.5	173,129	190,313.8
		Normal	147,048.3	161,625	177,659.3
		Mild	136,556.8	150,084.3	164,964.5
	BPHS_2	Severe	68,235.09	74,930.44	82,295.33
		Normal	63,753	70,000.14	76,872
		Mild	59,256.62	65,054.12	71,431.38
	BPHS_3	Severe	23,599.4	25,831.18	28,286.14
		Normal	22,105.37	24,187.75	26,478.37
		Mild	20,606.57	22,539.07	24,664.83
95%	BPHS_1	Severe	157,869.8	173,492.3	190,677.1
		Normal	147,411.6	161,988.3	178,022.6
		Mild	136,920.1	150,447.6	165,327.8
	BPHS_2	Severe	68,598.39	75,293.74	82,658.63
		Normal	64,116.3	70,363.44	77,235.3
		Mild	59,619.92	65,417.42	71,794.68
	BPHS_3	Severe	23,962.7	26,194.48	28,649.44
		Normal	22,468.67	24,551.05	26,841.67
		Mild	20,969.87	22,902.37	25,028.13
99%	BPHS_1	Severe	158,551.2	174,173.7	191,358.5
		Normal	148,093	162,669.7	178,704
		Mild	137,601.5	151129	166,009.2
	BPHS_2	Severe	69,279.89	75,975.24	83,340.13
		Normal	64,797.8	71,044.94	77,916.8
		Mild	60,301.42	66,098.92	72,476.18
	BPHS_3	Severe	24,644.2	26,875.98	29,330.94
		Normal	23,150.17	25,232.55	27,523.17
		Mild	21,651.37	23,583.87	25,709.63

,

Table A2. Heat provisions undertaken by different peak-shaving heating sources under α = 0.55, GJ.

Percentile	Heat Source	“Freezing Degree”	Period 1	Period 2	Period 3
1%	BPHS_1	Severe	119,844.1	132,061.2	145,499.9
		Normal	111,664.1	123,063.2	135,602.1
		Mild	103,468.1	114,047.5	125,684.9
	BPHS_2	Severe	50,032.4	55,268.27	61,027.74
		Normal	46,526.69	51,411.99	56,785.83
		Mild	43,014.09	47,548.13	52,535.58
	BPHS_3	Severe	15,126.57	16,871.86	18,791.68
		Normal	13,958	15,586.43	17,377.71
		Mild	12,787.13	14,298.48	15,960.96
5%	BPHS_1	Severe	120,525.5	132,742.6	146,181.3
		Normal	112,345.5	123,744.6	136,283.5
		Mild	104,149.5	114,728.9	126,366.3
	BPHS_2	Severe	50,713.9	55,949.77	61,709.24
		Normal	47,208.19	52,093.49	57,467.33
		Mild	43,695.59	48,229.63	53,217.08
	BPHS_3	Severe	15,808.07	17,553.36	19,473.18
		Normal	14,639.5	16,267.93	18,059.21
		Mild	13,468.63	14,979.98	16,642.46
10%	BPHS_1	Severe	120,888.8	133,105.9	146,544.6
		Normal	112,708.8	124,107.9	136,646.8
		Mild	104,512.8	115,092.2	126,729.6
	BPHS_2	Severe	51,077.2	56,313.07	62,072.54
		Normal	47,571.49	52,456.79	57,830.63
		Mild	44,058.89	48,592.93	53,580.38
	BPHS_3	Severe	16,171.37	17,916.66	19,836.48
		Normal	15,002.8	16,631.23	18,422.51
		Mild	13,831.93	15,343.28	17,005.76
90%	BPHS_1	Severe	123,452	135,669.1	149,107.8
		Normal	115,272	126,671.1	139,210
		Mild	107,076	117,655.4	129,292.8
	BPHS_2	Severe	53,640.3	58,876.17	64,635.64
		Normal	50,134.59	55,019.89	60,393.73
		Mild	46,621.99	51,156.03	56,143.48
	BPHS_3	Severe	18,734.47	20,479.76	22,399.58
		Normal	17,565.9	19,194.33	20,985.61
		Mild	16,395.03	17,906.38	19,568.86
95%	BPHS_1	Severe	123,815.3	136,032.4	149,471.1
		Normal	115,635.3	127,034.4	139,573.3
		Mild	107,439.3	118,018.7	129,656.1
	BPHS_2	Severe	54,003.6	59,239.47	64,998.94
		Normal	50,497.89	55,383.19	60,757.03
		Mild	46,985.29	51,519.33	56,506.78
	BPHS_3	Severe	19,097.77	20,843.06	22,762.88
		Normal	17,929.2	19,557.63	21,348.91
		Mild	16,758.33	18,269.68	19,932.16
99%	BPHS_1	Severe	124,496.7	136,713.8	150,152.5
		Normal	116,316.7	127,715.8	140,254.7
		Mild	108,120.7	118,700.1	130,337.5
	BPHS_2	Severe	54,685.1	59,920.97	65,680.44
		Normal	51,179.39	56,064.69	61,438.53
		Mild	47,666.79	52,200.83	57,188.28
	BPHS_3	Severe	19,779.27	21,524.56	23,444.38
		Normal	18,610.7	20,239.13	22,030.41
		Mild	17,439.83	18,951.18	20,613.66

and

Table A3. Heat provisions undertaken by different peak-shaving heating sources under α = 0.6, GJ.

Percentile	Heat Source	“Freezing Degree”	Period 1	Period 2	Period 3
1%	BPHS_1	Severe	90,304.58	99,567.65	109,757.1
		Normal	84,107.47	92,750.85	102,258.6
		Mild	77,907.64	85,931.04	94,756.77
	BPHS_2	Severe	37,372.62	41,342.52	45,709.4
		Normal	34,716.72	38,421.02	42,495.76
		Mild	32,059.64	35,498.24	39,280.7
	BPHS_3	Severe	10,906.64	12,229.94	13,685.57
		Normal	10,021.34	11,256.11	12,614.35
		Mild	9,135.652	10,281.85	11,542.67
5%	BPHS_1	Severe	90,986.08	100,249.1	110,438.5
		Normal	84,788.97	93,432.35	102,940
		Mild	78,589.14	86,612.54	95,438.27
	BPHS_2	Severe	38,054.12	42,024.02	46,390.9
		Normal	35,398.22	39,102.52	43,177.26
		Mild	32,741.14	36,179.74	39,962.2
	BPHS_3	Severe	11,588.14	12,911.44	14,367.07
		Normal	10,702.84	11,937.61	13,295.85
		Mild	9817.146	10,963.35	12,224.17
10%	BPHS_1	Severe	91,349.38	100,612.4	110,801.8
		Normal	85,152.27	93,795.65	103,303.3
		Mild	78,952.44	86,975.84	95,801.57
	BPHS_2	Severe	38,417.42	42,387.32	46,754.2
		Normal	35,761.52	39,465.82	43,540.56
		Mild	33,104.44	36,543.04	40,325.5
	BPHS_3	Severe	11,951.44	13,274.74	14,730.37
		Normal	11,066.14	12,300.91	13,659.15
		Mild	10,180.45	11,326.65	12,587.47
90%	BPHS_1	Severe	93,912.48	103,175.6	113,365
		Normal	87,715.37	96,358.75	105,866.5
		Mild	81,515.54	89,538.94	98,364.67
	BPHS_2	Severe	40,980.52	44,950.42	49,317.3
		Normal	38,324.62	42,028.92	46,103.66
		Mild	35,667.54	39,106.14	42,888.6
	BPHS_3	Severe	14,514.54	15,837.84	17,293.47
		Normal	13,629.24	14,864.01	16,222.25
		Mild	12,743.55	13,889.75	15,150.57
95%	BPHS_1	Severe	94,275.78	103,538.9	113,728.3
		Normal	88,078.67	96,722.05	106,229.8
		Mild	81,878.84	89,902.24	98,727.97
	BPHS_2	Severe	41,343.82	45,313.72	49,680.6
		Normal	38,687.92	42,392.22	46,466.96
		Mild	36,030.84	39,469.44	43,251.9
	BPHS_3	Severe	14,877.84	16,201.14	17,656.77
		Normal	13,992.54	15,227.31	16,585.55
		Mild	13,106.85	14,253.05	15,513.87
99%	BPHS_1	Severe	94,957.28	104,220.3	114,409.7
		Normal	88,760.17	97,403.55	106,911.2
		Mild	82,560.34	90,583.74	99,409.47
	BPHS_2	Severe	42,025.32	45,995.22	50,362.1
		Normal	39,369.42	43,073.72	47,148.46
		Mild	36,712.34	40,150.94	43,933.4
	BPHS_3	Severe	15,559.34	16,882.64	18,338.27
		Normal	14,674.04	15,908.81	17,267.05
		Mild	13,788.35	14,934.55	16,195.37

. in Appendix A.

In reality, the biofuel availability for a BDHS is affected by a range of events, such as the local biomass-harvest level, the forest residue accessibility, biomass-pellet-fuel production capacity, and the within-year meteorological condition. Thus, a kind of stochastic-interval uncertainty is utilized herein to represent the practical biofuel amount corresponding to a certain available level (i.e., a scarce, medium, or abundant level). The specific biofuel amounts under different available levels are provided in

Table 2. Biofuel available amount, 10³ t [5].

Biofuel Available Level	Period 1	Period 2	Period 3
Abundant (pba = 0.3)	[9.60, 9.97]	[11.53, 11.94]	[13.04, 13.52]
Medium (pba = 0.45)	[11.92, 12.36]	[13.85, 14.29]	[15.37, 15.91]
Scarce (pba = 0.25)	[13.27, 13.65]	[15.14, 15.56]	[17.28, 17.81]

. When the supplied biomass-pellet fuel amount cannot meet the demands of the heat sources, the biofuel deficit would occur, and the BM manager has to turn to the other markets, causing the biofuel deficit cost obviously higher than normal. The biomass-pellet fuel prices in normal conditions are [310, 350], [332, 361], and [341, 372] CNY·t⁻¹, corresponding to planning periods 1, 2, and 3, respectively. Related economic parameters are presented in

Table A4. Related economic parameters.

Economic Parameter	Heat Source	Period 1	Period 2	Period 3
Biofuel deficit price for a BPHS, CNY·t−1	BPHS_1	[410, 435]	[415, 440]	[420, 445]
	BPHS_2	[400, 417]	[405, 423]	[410, 437]
	BPHS_3	[385, 405]	[395, 412]	[405, 424]
Heat supply price of a PHS, CNY·GJ−1	BPHS_1	[0.81, 0.92]	[0.83, 0.94]	[0.85, 1.01]
	BPHS_2	[0.82, 0.91]	[0.85, 0.96]	[0.87, 1.08]
	BPHS_3	[0.94, 1.07]	[0.99, 1.13]	[1.04, 1.16]
Pollutant removal price in a PHS, CNY·GJ −1	BPHS_1	[0.95, 1.11]	[0.97, 1.13]	[0.99, 1.15]
	BPHS_2	[0.99, 1.15]	[1.01, 1.17]	[1.03, 1.20]
	BPHS_3	[1.02, 1.19]	[1.04, 1.22]	[1.07, 1.25]

in Appendix A. The heating value of the biomass-pellet fuel is [17.02, 17.83] GJ·t⁻¹.

Furthermore, the existing heating capacities of the BPHSs may not sufficiently meet the increasing heat-load demand owing to the incremental heating area year by year; and beyond that, the potential thermalization coefficient variation could also result in the undertaken heat-load variation of BPHSs, causing the heat-load mismatch between the supply and demand sides and indirectly impacting BM. Thus, the heating capacity expansion issue is incorporated into the modeling framework. The existing heat capacities of BPHSs_1, _2, and _3 are, respectively, 29MW, 14MW, and 7MW, and the corresponding heat loads to be undertaken are provided in

Table A5. Heating load undertaken by different peak-shaving heating sources, MW.

Thermalization Coefficient	Heat Source	“Freezing Degree”	Period 1	Period 2	Period 3
α = 0.5	BPHS_1	Severe	37.77948	38.92424	40.10368
		Normal	34.34498	35.38567	36.45789
		Mild	30.91048	31.8471	32.8121
	BPHS_2	Severe	16.19121	16.68182	17.18729
		Normal	14.71928	15.16529	15.62481
		Mild	13.24735	13.64876	14.06233
	BPHS_3	Severe	5.397069	5.560605	5.729097
		Normal	4.906426	5.055096	5.20827
		Mild	4.415783	4.549586	4.687443
α = 0.55	BPHS_1	Severe	34.00153	35.03181	36.09331
		Normal	30.91048	31.8471	32.8121
		Mild	27.81944	28.66239	29.53089
	BPHS_2	Severe	14.57209	15.01363	15.46856
		Normal	13.24735	13.64876	14.06233
		Mild	11.92262	12.28388	12.6561
	BPHS_3	Severe	4.857362	5.004545	5.156187
		Normal	4.415783	4.549586	4.687443
		Mild	3.974205	4.094627	4.218699
α = 0.6	BPHS_1	Severe	30.22358	31.13939	32.08294
		Normal	27.47599	28.30854	29.16631
		Mild	24.72839	25.47768	26.24968
	BPHS_2	Severe	12.95296	13.34545	13.74983
		Normal	11.77542	12.13223	12.49985
		Mild	10.59788	10.91901	11.24986
	BPHS_3	Severe	4.317655	4.448484	4.583278
		Normal	3.925141	4.044076	4.166616
		Mild	3.532627	3.639669	3.749954

in Appendix A. The available expansion choices are, respectively, 7MW, 14 MW, and 28 MW for three BPHSs, and the corresponding costs are [600, 720] × 10³, [1100, 1250] × 10³, and [2000, 2300] × 10³ CNY.

3.2. Modeling Formulation

Based on the overview of the investigated BDHS, different uncertain inputs can be identified, quantified, and expressed, and the developed ITCS-MIP model is suitable for coping with the practical BM problem mentioned above. The modeling process based on the ITCS-MIP framework is presented in Figure 3, and the specific model objective and constraints are given hereunder:

$$\mathit{Min}{f}_{cost}^{\pm }=f_1^{\pm }+f_2^{\pm }+f_3^{\pm }+f_4^{\pm }$$

(13)

$$f_1^{\pm }={\displaystyle \sum _{i=1}^I{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{p}_k^{hd}\cdot CC{B}_t^{\pm }\cdot XC{B}_{ikt}^{\pm }}}}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{p}_j^{ba}\cdot CD{B}_{it}^{\pm }\cdot p_k^{hd}\cdot }XD{B}_{ijkt}^{\pm }}}}$$

(14)

$$f_2^{\pm }={\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^{T}Q{b}^{\pm }\cdot l{b}_i^{\pm }\cdot CSP{H}_{it}^{\pm }\cdot (p_j^{ba}\cdot p_k^{hd}\cdot XD{B}_{ijkt}^{\pm }+p_k^{hd}\cdot XC{B}_{ikt}^{\pm })}}}}$$

(15)

$$f_3^{\pm }={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{i=1}^{I}Y{D}_{im}^{\pm }\cdot CG{E}_m^{\pm }}}$$

(16)

$$f_3^{\pm }={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{i=1}^{I}Y{D}_{im}^{\pm }\cdot CG{E}_m^{\pm }}}$$

(17)

subject to:

Constraints for the biomass-pellet fuel allocation to BPHSs:

$${\displaystyle \sum _{i=1}^I(\frac{TSH{T}_{ikt}^{a,1-pr}}{Q{b}^{\pm }\cdot l{b}_i^{\pm }}-XD{B}_{ijkt}^{\pm }})\le AVB{F}_{jt}^{\pm },\forall j,k,t$$

(18)

$$XD{B}_{ijkt}^{\pm }\cdot Q{b}^{\pm }\cdot l{b}_i^{\pm }\le {\eta }^{\pm }\cdot TSH{T}_{ikt}^{a,pr},\forall i,j,k,t$$

(19)

$$TSH{T}_{ikt}^{a,1-pr}\le XC{B}_{ikt}^{\pm }\cdot Q{b}^{\pm }\cdot l{b}_i^{\pm },\forall i,k,t$$

(20)

$$0\le XD{B}_{ijkt}^{\pm },\forall i,j,k,t$$

(21)

$$0\le XC{B}_{ikt}^{\pm },\forall i,k,t$$

(22)

Constraints for the heating capacity expansion:

$$MH{L}_{ikt}^{\alpha }\le {\displaystyle \sum _{m=1}^{M}Y{D}_{im}^{\pm }}\cdot GE{H}_m^{}+EH{C}_i,\forall i,k,t$$

(23)

$${\displaystyle \sum _{m=1}^{M}Y{D}_{im}^{\pm }\le 1,}\forall i$$

(24)

$$Y{D}_{im}^{\pm }=\left\{\begin{array}{l}1;ifheatingcapacityexpansionisundertaken\\ 0;ifotherwise\end{array}\right\}$$

(25)

Nonnegative constraints:

$$0\le XC{B}_{ikt}^{\pm },XD{B}_{ijkt}^{\pm },Y{D}_{im}^{\pm },\forall i,j,k,t,m$$

(26)

The detailed nomenclature for decision variables and modeling parameters is provided in Nomenclature.

4. Result Analysis and Discussion

The ITCS-MIP model developed in this research work integrates a variety of uncertain optimization techniques, such as the TDSP, IPP, CCP, and MIP methods, to handle the multiple and polymorphic uncertainties. Consequently, the results of the BM scheme (i.e., biofuel deficit assignment) and its derivative (i.e., heating capacity expansion) can be obtained in terms of the solutions under different scenario combinations (i.e., the biofuel availability and “freezing degree” level of a heating season). Moreover, the detailed sensitivity of various system conditions, including the violation-risk probability level, planning period, and thermalization coefficient, can also be examined on the obtained results for the in-depth analysis. The representative thermalization coefficient values (α) are set to 0.5, 0.55, and 0.6 based on the latest “China Energy Saving Law”, and the control coefficient for biofuel deficit (i.e., η) is [0.45, 0.47] in the constraints.

4.1. Result Analysis

Figure 4, Figure 5 and Figure 6 jointly present the biofuel-deficit allotment pattern among three BPHSs under various scenario combinations with different system conditions, indicating that the biofuel deficits of different BPHSs would vary markedly due to the uncertain modeling inputs and the temporal and spatial variations of the system conditions. In detail, the biofuel deficit fluctuates significantly in BPHS_1 under different scenario combinations, while the fluctuations of BPHSs_2 and _3 are relatively insignificant. More specifically, on one hand, with the “freezing degree” of the heating season changing from the mild level to the severe, the biofuel deficits of the three BPHSs would increase with different magnitudes. Among them, BPHS_1 has the highest deficit with the largest amplitude of variation, followed by BPHS_2, and the deficit of BPHS_3 is the lowest. For instance, under the system conditions of α = 0.5, Pr = 0.01 and the scarce biofuel availability, the biofuel deficits of BPHS_1 would be, respectively, [3870.34, 4337.59], [2879.69, 3291.18], and [1885.89, 2241.48] tons, corresponding to severe, normal, and mild “freezing-degree” in Period 2; the corresponding deficits of BPHS_2 would be [2129.79, 2360.12], [1982.59, 2196.97], and [1834.87, 2033.30] tons; while the deficits of BPHS_3 would be [671.61, 744.30], [621.95, 689.28] and [572.12, 634.04] tons, respectively. That is mainly because, with the outdoor temperature declining (when season “freezing-degree” varying from the “mild” to the “normal” and then to the “severe”), more biofuels should be fed into BPHSs for satisfying the increasing heat demand when the other system conditions are unchanged. Moreover, since the penalty cost of BPHS_3 is the lowest for the biofuel deficit, and the cost of BPHS_2 ranks the second, the deficits of the BPHSs_3 and _2 should be assigned initially for economic optimality when the biofuel availability is at a fixed level. However, due to the deficit control constraints [i.e., in Equation (19)] existing in the ITCS-MIP model, the deficit of BPHSs_2 and _3 would reach their control-constraint boundaries when the “severe” heating season appears, and thus more biofuel deficits should be assigned to the BPHS_1; by contrast, during some “normal” or “mild” heating seasons, the deficit in BPHSs_2 and _3 would still be within their permissible constraint ranges, causing a lower or even no biofuel deficit in BPHS_1.

On the other hand, with the biofuel availability changing from the scarce level to the abundant, the deficit result of BPHS_1 shows a significant downward trend, while the deficit variations of BPHSs_2 and _3 are insignificant. Specifically, as shown in Figure 4, the biofuel deficit of BPHS_1 would be [3554.45, 3935.28], [1234.45, 1545.28], and [0, 255.28] tons corresponding to the scarce, medium, and abundant biofuel availability at the normal “freezing degree” level in Period 1; the corresponding deficits of BPHS_2 would be [1796.02, 1990.29], [1796.02, 1990.29], and [1680.47, 1990.25] tons, respectively; while the deficits of BPHS_3 would remain unchanged, which is [559.02, 619.52] tons. It is mainly due to the fact that when the biofuel demand and its influencing factors are at the fixed levels, the increasing biofuel supply will lead to a decrease in the total biofuel deficit in the BDHS. Subsequently, considering the deficit penalty difference among the BPHSs mentioned above, a large decline of the biofuel deficit in BPHS_1 would thus occur with the biofuel availability changing from the scarce level to the abundant, and the deficit in BPHS_1 may be even less than those in BPHSs_2 and _3 at the medium and abundant levels, as Figure 4, Figure 5 and Figure 6 show.

Figure 4 also examines the biofuel deficit variation trend from Period 1 to Period 3 under various scenario combinations, indicating that the biofuel deficit in BPHS_1 would decline (or decline to even no deficit under some scenarios) observably while the corresponding biofuels in BPHSs_2 and _3 would rise slightly. For example, from Period 1 to Period 3, the deficit in BPHS#1 would be [1234.45, 1545.28], [559.69, 941.18], and [420.46, 779.67] tons under the scenario combination of the medium biofuel availability and at the normal “freezing degree” level, respectively; correspondingly, the deficit in BPHS#2 would be [1796.02, 1990.29], [1982.57, 2196.97], and [2187.77, 2424.37] tons and the deficit in BPHS#3 would be [559.02, 619.51], [621.95, 689.26], and [691.17, 765.98] tons. In fact, the biofuel demands could grow over time with the heating area increasing, but the biofuel supply growth could be stronger than the demands due to the subsidy policy and market effect. The difference between the supply and demand would thus lead to the total biofuel deficit decreasing, and the downward variation trend of the deficit in BPHS_1 would be consequently dramatic when the biofuel deficits in BPHSs_2 and _3 reach their corresponding constraint bounds.

Figure 5 describes the biofuel deficit variation trend of different BPHSs under various scenario combinations with different Pr levels at α = 0.5 in Period 1. It can be found that, with the constraint–violation risk level increases (i.e., the Pr level increases), the biofuel deficit would be lowered in BPHS_1, and the deficit in BPHSs_2 and _3 would mount up slightly under most scenario combinations. For instance, when the Pr level varying from 0.01 to 0.1, the deficit in BPHS#1 would be [1234.45, 1545.28], [1058.09, 1355.96], and [964.67, 1255.02] tons under the scenario combination of the medium biofuel availability and at the normal “freezing degree” level, respectively. In comparison, the deficits in BPHSs_2 and _3 would be correspondingly [1796.02, 1990.25] and [559.02, 619.52] tons, [1816.37, 2012.8] and [579.61, 642.34], [1827.22, 2024.82] and [590.59, 654.50] tons. This can be mainly explained by the fact that the increased constraint–violation risk level (i.e., Pr) would lower the totaling heat provisions and the associated biofuel deficit under each scenario combination. Thereafter, when the penalty difference of the BPHSs makes the biofuel deficit preferentially allocated to BPHSs_2 and _3, the biofuel deficit in BPHS_1 would decrease dramatically when the Pr level raised, while the corresponding deficits in BPHSs_3 and _2 would increase slightly owing to the controlling effect of the deficit-associated constraint (i.e., in Equation (19)). It is worth noting that, although a lower Pr level would result in a lower risk of biofuel shortage and a higher heating satisfaction degree, there would be a potential waste of biofuel and heat supply when the “freezing degree” of heating season is mild, and the biofuel availability level is abundant.

Figure 6 shows the biofuel deficit variation trend of different BPHSs in Period 2 under various scenario combinations with different α values at Pr = 0.95. It can be seen that, with the α value increasing from 0.5 to 0.6, the biofuel deficits in three BPHSs would all decline significantly. For instance, at α = 0.5, the biofuel deficits in BPHSs_1, _2, and _3 would be [1709.52, 2052.11], [1855.22, 2055.86], and [592.72, 656.87] tons under the scenario combination of the scarce biofuel availability and the mild “freezing degree” level, respectively; at α = 0.55, the corresponding deficits would be 0, [357.18, 698.34], and [452.71, 501.71] tons under the same scenario combination; at α = 0.6, there would be no biofuel deficits in three BPHSs. Meanwhile, the results also show that the biofuel deficits in different BPHSs would become 0 tons under nearly half of the scenario combinations, especially the deficits in BPHS_1 under more than two-thirds of the scenario combinations. In reality, the thermalization coefficient (i.e., α) determines the heat-provision assignment between the MCHS and the BPHSs. The increase in the thermalization coefficient means that the heat provisions undertaken by the BPHSs would be lowered, which indirectly reduces the overall biofuel demand and the corresponding biofuel deficit in each BPHS. Beyond that, the penalty-cost difference results in the falling range of the biofuel deficit in BPHS_1 to be the largest among different BPHSs when the α value is adjusted to a higher level. Thus, a 0 biofuel deficit would appear frequently in BPHS_1 under different scenario combinations.

Figure 7 presents the heating-capacity expansion result for each BPHS. The results indicate that the expansion schemes of the three BPHSs would be different from each other under varied thermalization coefficients (i.e., α). In detail, BPHS_1 would reach the expansion amounts of 14 MW, 14 MW, and 7 MW under α = 0.5, 0.55, and 0.6, which would have the highest capacity–expansion amount in comparison with the other two BPHSs; whereas BPHS_3 would not expand over the entire planning horizon, along with the α value varying from 0.5 to 0.6; at the same time, the capacity–expansion amount of BPHS_2 is intermediate among that of three BPHSs. The various capacity–expansion schemes for three BPHSs are because the gap between the actual peak heat-load demand of each community and the existing heating capacity of the corresponding BPHS is different. Beyond that, the arising α value can reduce the actual peak heat-load demand of a community, and a relatively low capacity–expansion is able to compensate for the heat-load shortage for each BPHS. Consequently, as mentioned above, the expansion amounts of BPHS_1 and _2 are gradually lowered with the α value increasing.

Figure 8 shows the solution of the objective function value representing the total heating-cost expectation under different thermalization coefficients (i.e., α) and constraint–violation risk levels (i.e., Pr) over the whole planning period. It can be found that all the solutions under different system conditions can be obtained in terms of interval numbers, demonstrating that the developed model is valid and sensitive to the uncertain modelling inputs. With the α or Pr values increasing, the heating cost would fall to a different extent. On one hand, under Pr = 0.05, the heating cost corresponding to α = 0.5, 0.55, and 0.6 would be [24.44, 29.00] × 10⁶, [16.97, 20.01] × 10⁶, and [11.94, 14.01] × 10⁶ RMB, respectively. Due to the α value representing the heat provision assignment between the MCHS and the BPHSs, a high α value makes the heat provisions undertaken by BPHSs lowered, and the biofuel demand and the corresponding deficit would also decline, causing the final heating cost decreased. On the other hand, the heating cost would be lowered, along with the Pr value raising (i.e., with the heating satisfaction decreasing), for example, under α = 0.55, the heating cost would be [17.19, 20.27] × 10⁶, [16.97, 20.01] × 10⁶, and [16.86, 19.88] × 10⁶ RMB, respectively, corresponding to Pr = 0.01, 0.05, and 0.1, reflecting that a low Pr level would lower the constraint-violation risk by utilizing a relatively high system cost; conversely, a high Pr level would sacrifice the quality of the heating service (i.e., supply insufficient heat to communities) in order to reduce the heating cost. Therefore, in general, the adjustment of the Pr level could also reveal the decision-maker’s preference regarding the tradeoff between the heating cost and the constraint–violation risk.

4.2. Discussion

The practicality and validity of the established ITCS-MIP model can be verified by the real case we used, and the model can be applied widely to the similar district heating systems partly fueled by biomass-pellet fuels. However, in the context of the “reaching carbon neutral by 2060” strategy of China, this development direction encouraged by the government has driven remarkable technological progress in terms of heat storage and other renewable-energy heating technologies (e.g., wind, solar energy, or ground-source heat pumps) in China. Different renewable energy and heat storage technologies can be incorporated simultaneously into the traditional district heating systems [35,36]. The integrated heating systems tend to be more complex due to higher amounts of system components or modules (e.g., solar, wind energy heating equipment, and thermal storage devices) and their interactions lying in the system. Beyond that, the concomitant multiple and polymorphic uncertainties will also be introduced into the modelling system and could indirectly lead to changes in fuel–energy management patterns by affecting heat loads (or heat provisions). In fact, although the ITCS-MIP model has merits, including its high efficiency in obtaining the optimum solution under uncertainties and decreasing the calculation complexity, there are potential limitations that exist in the proposed model for the future of fuel–energy management and should be addressed in future research. For example, the developed model cannot deal with the fuzzy uncertainty embedded within the heating duration or the compound uncertainties originated from other renewable-energy heating technologies. Meanwhile, its solution may be one-sided, since only the single objective is considered instead of multiple objectives, and the impacts among different targets are ignored. The corresponding interactive complexity and uncertainties need to be identified, quantified, and handled by developing advanced inexact optimization techniques. Thus, other optimization techniques, such as multi-objective optimization, mixed-integer programming, dynamic programming, fuzzy optimization, and intelligent optimization algorithms, should be merged into the ITCS-MIP model framework to cope with the complicated energy-optimization problem within a multi-energy, multi-module, multi-interaction, and multi-uncertainty context. Moreover, considering that other single and compound uncertainties could exist in the future district heating system, factorial analysis, causal analysis, or other advanced techniques need to be applied to obtain the factor-interaction impacts on the model response.

5. Conclusions

In this paper, an inexact biofuel-management model is developed by integrating different uncertain programming techniques for real biofuel-based heating sources. The proposed model can address multiple and compound uncertainties lying in the system and generate the optimal biofuel management schemes, in terms of biofuel allocation planning and heating capacity expansion subject to supply–demand, policy requirement constraints, and the financial minimization objective. Beyond that, the model can also quantitatively analyze the conflict between economic targets that minimize the system cost and risk preference that maximize heating-service satisfaction.

Due to the penalty difference of the BPHSs and the complicated interaction model constraints, the results indicate that (1) with the biofuel availability changing from the scarce level to the abundant, the deficit result of BPHS_1 shows a significant downward trend while the deficit variations of BPHSs_2 and _3 are insignificant. (2) When the “freezing degree” of the heating season changes from the mild level to the severe, BPHS_1 has the highest deficit with the largest amplitude of variation, followed by BPHS_2, and the deficit of BPHS_3 is the lowest. (3) With the α value increasing from 0.5 to 0.6, the biofuel deficits in three BPHSs would all decline significantly. Meanwhile, BPHS_1 would reach the expansion amounts of 14 MW, 14 MW, and 7 MW, which would have the highest capacity–expansion amount in comparison with the other two BPHSs, whereas BPHS_3 would not be expanded over the entire planning horizon. (4) It can be found that, with the constraint–violation risk level increasing (i.e., the Pr level increases), the biofuel deficit would be lowered in BPHS_1, while the deficit in BPHSs_2 and _3 would mount up slightly under most scenario combinations. Moreover, a low constraint–violation risk (i.e., a high heating satisfaction level) would potentially lead to a high heating cost.