An Optimization Model of Carbon Sinks in CDM Forestry Projects Based on Interval Linear Programming

Abstract

This study describes the first general optimization model for complex systems with uncertain parameters and decision variables represented as intervals in CDM forestry projects. We work through a specific example of the optimization method developed for a Clean Development Mechanism (CDM) forestry project in Inner Mongolia, China. This model is designed to optimize the carbon sink capacity of the new forests, and can deal with uncertainties in the carbon sink capacity, average annual rainfall, ecological parameters, and biological characteristics of tree species. The uncertain inputs are presented in the form of intervals, as are several of the optimized output variables. Compared with the project’s originally recommended scheme, the optimized model will absorb and fix between 1,142 and 885,762 tonnes of extra carbon dioxide. Moreover, the ecological and environmental benefits of the project are also raised to various extents.

1. Introduction

The accumulation of greenhouse gases (GHGs) such as carbon dioxide in the atmosphere is one of the most serious environmental problems currently faced by human beings. Actions taken to reduce GHGs employ one of two methods: reducing emissions at the source, or increasing absorption from the atmosphere. A carbon sink is a process, activity, or mechanism that eliminates carbon dioxide from the atmosphere. Through reforestation, forest management, and protection, carbon dioxide can be absorbed from the atmosphere and fixed in vegetation and soil, thereby reducing the atmospheric concentration [1]. Such activities are referred to as forest carbon sinks. In the global terrestrial ecosystem, forests play an important and irreplaceable role in the maintenance of balance and the prevention of ecological crisis [2,3,4,5,6,7]. As the greenhouse effect and global warming intensify, forest carbon sinks have drawn increasing and extensive attention from the international community [8,9].

As forest carbon sinks have the ability to mitigate climate warming, afforestation and reforestation project activities have been added to the Clean Development Mechanism (CDM-AR). CDM-AR now certifies projects so that the carbon sequestered can be sold to interested buyers in Annex 1 countries [10]. This win-win arrangement not only helps developed countries implement their obligations at lower cost, but promotes the sustainable development of community economies and forestry in developing countries. More importantly, such carbon sink transactions transform the carbon dioxide fixed by forests into a commodity, so the value of forest carbon sinks is set by the market. Afforestation and reforestation projects have therefore become one of several human activities recognized by the Intergovernmental Panel on Climate Change (IPCC).

Many works over the past decades have explored forestation projects. For example, Van Elegem et al. [11] proposed a multiple criteria decision-making method to deal with the allocation of new urban forests. Espelta et al. [12] applied a discrete multi-criteria method to choose among alternative approaches in a post-fire reforestation program in Spain. Gilliams et al. [13] compared analytic hierarchy processes (AHP) with a discrete multi-criteria method to choose the best afforestation method in Belgium.

In real-world forestation problems, most of the data available for decision-making are highly uncertain, and their interrelationships can be extremely complicated [14]. Likewise, in CDM forestry projects, many uncertainties exist in the system parameters and their interrelationships. For example, spatial and temporal variations may exist in the number of carbon sinks, the average annual rainfall, various ecological parameters, and the biological characteristics of tree species. These uncertainties and complexities are further compounded not only by interactions among the parameters but also by additional economic and ecological implications.

In CDM forestry projects, uncertainty estimates are often limited to upper and lower bounds, so the parameters can be expressed as interval numbers. The quality of the information available in many practical problems is simply not good enough to present the system parameters as deterministic numbers.

The interval linear programming (ILP) method is a feasible approach for this type of planning problem. Many researchers have applied interval programming in other fields [15,16,17,18,19,20,21]. For example, Huang and Baetz provided a grey linear programming approach for municipal solid waste management planning under uncertainty [15].Huang developed an interval parameter water quality management model for water resources decision making under uncertainty [17]. However, to date this method has not been attempted in the context of the CDM forestry projects.

A CDM project investor desires the rapid increase of carbon sinks, so may prefer exotic tree species with fast growth rate and strong carbon sequestration abilities. However, carbon sequestration ability cannot be the sole criterion for tree species selection. The community implementing the project will want to promote sustainable development while improving the local economy and the ecological environment. As a result, they may prefer local tree species with important economic or ecological benefits. The effects of the selected tree species on the socio-economic objectives and environment must be comprehensively examined with an eye to sustainable development for the implementing party. Furthermore, according to the rules of CDM project, whether or not a CDM forestry project passes approval and certification mainly depends on the diversity of its benefits. That is, a successful CDM project should not only create carbon sinks to mitigate climate change, but also promote the development of the community, protect biodiversity, and so on. Thus, it is necessary to search for the optimal combination of parameters that brings the maximum carbon sink capacity of the project subject to suitable environmental and socio-economic standards in all categories. Unfortunately, no attempts have been made in this field to systematically determine the true optimum for a given project.

This paper develops an optimization model based on the ILP method for CDM forestry projects with uncertain parameters. This model of carbon sinks can deal with uncertainties (presented as intervals) in the total carbon sink capacity, average annual rainfall, the ecological parameters of the site, and the biological characteristics of proposed tree species. Its output is an afforestation plan which is flexible to possible variations in the system or deviations from the actual conditions caused by the input uncertainties. A case study will be provided to demonstrate the applicability of the optimization model. Our results indicate that it is indeed possible to find project solutions that maximize the carbon sinks while bringing additional benefits to the ecological environment and socio-economic objectives.

2. Methodology

In this section, we first provide the ILP formulation and its solution algorithm [15,16,17,18,19,20,21], and then provide the modelling formulation.

2.1. Interval Linear Programs [16,17]

Let denote a set of interval numbers. An interval linear programming (ILP) model can be defined as follows:

$$Max\hspace{0.17em}{f}^{\pm }=\hspace{0.17em}\hspace{0.17em}{C}^{\pm }{X}^{\pm }$$

(1a)

subject to:

$$A^{\pm }{X}^{\pm }\le B^{\pm }$$

(1b)

$$X^{\pm }\ge 0$$

(1c)

where $A^{\pm }\in {\left\{R^{\pm }\right\}}^{m\times n}$, $B^{\pm }\in {\left\{R^{\pm }\right\}}^{m\times 1}$, $C^{\pm }\in {\left\{R^{\pm }\right\}}^{1\times n}$, and $X^{\pm }\in {\left\{R^{\pm }\right\}}^{n\times 1}$.

2.2. Solution of the ILP Model

An interactive algorithm was developed by Huang et al. [15] to solve the ILP model by analyzing the interrelationships between parameters and variables. This type of solution may be favoured by decision-makers because of its flexibility. As our method does not lead to complicated intermediate models, it also has reasonable computational requirements. According to the algorithm, model (1) can be solved through a two-step process. The first step is to formulate a submodel corresponding to $f^{+}$, and solve it by maximizing the objective. The second step is to solve a submodel corresponding to $f^{-}$, based on the upper bound solution generated in the first step.

The submodel corresponding to $f^{+}$ is formulated as follows (assuming that $b_i^{\pm }\ge 0$):

$$\mathrm{Max}\hspace{0.17em}{f}^{+}=\hspace{0.17em}{\displaystyle \sum _{j=1}^{k_1}{c}_j^{+}}{x}_j^{+}\hspace{0.17em}+{\displaystyle \sum _{j=k_1+1}^n{c}_j^{+}{x}_j^{-}}$$

(2a)

subject to:

$${\displaystyle \sum _{j=1}^{k_1}{\left|a_{ij}\right|}^{-}\mathrm{sign}\left(a_{ij}^{-}\right)}{x}_j^{+}+{\displaystyle \sum _{j=k_1+1}^n{\left|a_{ij}\right|}^{+}\mathrm{sign}\left(a_{ij}^{+}\right)x_j^{-}}\le b_i^{+},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall i$$

(2b)

$$x_j^{+}\ge 0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall j$$

(2c)

Solutions of $x_{j\hspace{0.17em}\mathrm{opt}}^{+}$ ($j=1,\hspace{0.17em}2,\hspace{0.17em}...,k_1$), $x_{j\hspace{0.17em}\mathrm{opt}}^{-}$ ($j=k_1+1,\hspace{0.17em}{k}_1+2,\hspace{0.17em}...,n$), and $f_{\mathrm{opt}}^{+}$ can be obtained using submodel (2). The submodel corresponding to $f^{-}$ can then be formulated as follows (assuming that $b_i^{\pm }\ge 0$):

$$\mathrm{Max}\hspace{0.17em}{f}^{-}=\hspace{0.17em}{\displaystyle \sum _{j=1}^{k_1}{c}_j^{-}}{x}_j^{-}\hspace{0.17em}+{\displaystyle \sum _{j=k_1+1}^n{c}_j^{-}{x}_j^{+}}$$

(3a)

subject to:

$${\displaystyle \sum _{j=1}^{k_1}{\left|a_{ij}\right|}^{+}\mathrm{sign}\left(a_{ij}^{+}\right)}{x}_j^{-}+{\displaystyle \sum _{j=k_1+1}^n{\left|a_{ij}\right|}^{-}\mathrm{sign}\left(a_{ij}^{-}\right)x_j^{+}}\le b_i^{-},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall i$$

(3b)

$$x_j^{+}\ge 0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall j$$

(3c)

$$x_j^{-}\le x_{j\hspace{0.17em}\mathrm{opt}}^{+}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\hspace{0.17em}2,\hspace{0.17em}...,k_1$$

(3d)

$$x_j^{+}\ge x_{j\hspace{0.17em}\mathrm{opt}}^{-}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=k_1+1,\hspace{0.17em}{k}_1+2,\hspace{0.17em}...,n$$

(3e)

Solutions of $x_{j\hspace{0.17em}\mathrm{opt}}^{-}$ ($j=1,\hspace{0.17em}2,\hspace{0.17em}...,k_1$), $x_{j\hspace{0.17em}\mathrm{opt}}^{+}$ ($j=k_1+1,\hspace{0.17em}{k}_1+2,\hspace{0.17em}...,n$), and $f_{\mathrm{opt}}^{-}$ can be obtained using submodel (3). Thus, the general solutions can be obtained as follows:

$$f_{opt}^{\pm }=[f_{opt}^{-},f_{opt}^{+}]$$

(4)

$$x_{j\hspace{0.17em}opt}^{\pm }=[x_{j\hspace{0.17em}opt}^{-},x_{j\hspace{0.17em}opt}^{+}], \forall j$$

(5)

If the objective function is to be minimised, then the submodel corresponding to $f^{-}$ should be solved first.

2.3. The CDM Forestry Model

In CDM forestry projects, uncertainties exist in many system parameters and their interrelationships. For example, spatial and temporal variations may exist in the total carbon sink capacity, the average annual rainfall, various ecological parameters, and the biological characteristics of tree species. Given upper and lower bounds, all of these variables can be expressed as interval numbers. In many practical problems, the quality of the available information permits no more informative constraints.

We consider a carbon sink model with four possible planting patterns: pure forest, pure rotation forest, mixed forest, and mixed rotation forest. The objective is to maximize the total number of carbon sinks in CDM forestry projects subject to the decision variables and the relationships between the decision variables and the objective. A complete list of decision variables is provided after the model.

The optimization model can be formulated as follows:

$$\mathrm{Max} f^{\pm }={\displaystyle \sum _{i=1}^I{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{t=1}^T{D}_{imt}^{\pm }}}}\cdot E_{im}^{\pm }+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{t=1}^T{B}_{int}^{\pm }}}}\cdot A_{in}^{\pm }+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{t=1}^T{C}_{ijt}^{\pm }}}}\cdot S_{ij}^{\pm }+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{t=1}^T{F}_{iht}^{\pm }}}}\cdot G_{ih}^{\pm }$$

(6a)

where: $f^{\pm }$ = total carbon sink capacity (tonnes); $E_{im}^{\pm }$ = area of afforestation mode m in land type i when the planting pattern is pure forest (ha); $A_{in}^{\pm }$ = area of afforestation mode n in land type i when the planting pattern is pure rotation forest (ha); $S_{ij}^{\pm }$ = area of afforestation mode j in land type i when the planting pattern is mixed forest (ha); $G_{ih}^{\pm }$ = area of afforestation mode h in land type i when the planting pattern is mixed rotation forest (ha); t = computation period of the CDM forestry projects (yr); i = land type of the project (number of afforestation sites), i = 1,2,…,I; k = tree species selected for the project, k = 1,2,…,K; m = afforestation mode when the planting pattern is pure forest, m = 1,2,…,M; n = afforestation mode when the planting pattern is pure rotation forest, n = 1,2,…,N; j = afforestation mode when the planting pattern is mixed forest, $j$ = 1,2,…,J; h = afforestation mode when the planting pattern is mixed rotation forest, $h$ = 1,2,…,H; $D_{imt}^{\pm }$ = carbon sink capacity per unit area of afforestation mode m in land type i during the t^th year (tonnes/ha); $B_{int}^{\pm }$ = carbon sink capacity per unit area of afforestation mode n in land type i during the t^th year (tonnes/ha); $C_{ijt}^{\pm }$ = carbon sink capacity per unit area of afforestation mode j in land type i during the t^th year (tonnes/ha); $F_{iht}^{\pm }$ = carbon sink capacity per unit area of afforestation mode h in land type i during the t^th year (tonnes/ha).

subject to:

$${\displaystyle \sum _{i=1}^I{\displaystyle \sum _{k=1}^K\left\{\left(\begin{array}{l}{\displaystyle \sum _{m=1}^M{E}_{im}^{\pm }\cdot P{E}_{mk}}+{\displaystyle \sum _{n=1}^N{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }\left[\frac{T}{T_n}+1\right]\\ +{\displaystyle \sum _{j=1}^J{S}_{ij}^{\pm }\cdot P{S}_{jk}}+{\displaystyle \sum _{h=1}^H{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot \left[\frac{T}{T_h}+1\right]}\end{array}\right)\cdot {\rho }_k\cdot R_k^{\pm }\right\}}\le R^{\pm }}$$

(6b)

(Funding constraints)

where $P{E}_{mk}$ = presence of tree species k in afforestation mode m, takes the value 0 or 1 (pure forest); $P{A}_{nk}$ = presence of tree species k in afforestation mode n, takes the value of 0 or 1 (pure rotation forest); $P{S}_{jk}$ = proportion of tree species k in afforestation mode j (mixed forest); $P{G}_{hk}$ = proportion of tree species k in afforestation mode h (mixed rotation forest); $T_n$ = rotation of tree species k when the planting pattern is pure rotation forest (yr); $T_h$ = rotation of tree species h when the planting pattern is mixed rotation forest (yr); ${\rho }_k$ = planting density of tree species k (plants/ha); $R_k^{\pm }$ = seedling price of tree species k (million CNY/plant); $R^{\pm }$ = total investment in seedlings (million CNY).

$${\displaystyle \sum _{m=1}^M{E}_{im}^{\pm }+{\displaystyle \sum _{n=1}^N{A}_{in}^{\pm }}+{\displaystyle \sum _{j=1}^J{S}_{ij}^{\pm }+}}{\displaystyle \sum _{h=1}^H{G}_{ih}^{\pm }}\le A_i^{\pm }, \forall i$$

(6c)

$$A_i^{\pm }\le A_0^{\pm }$$

(6d)

(Total area constraints)

where $A_0^{\pm }$ = total afforestation area of the project (ha); $A_i^{\pm }$ = area of land type i (ha).

$${\displaystyle \sum _{m=1}^M{\displaystyle \sum _{i=1}^I{E}_{im}^{\pm }\cdot P{E}_{mk}}}+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{i=1}^I{A}_{in}^{\pm }\cdot P{A}_{nk}}}+{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{i=1}^I{S}_{ij}^{\pm }\cdot P{S}_{jk}}}+{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^I{G}_{ih}^{\pm }\cdot P{G}_{hk}}}\ge A{P}_k^{\pm },\hspace{0.17em} \forall k$$

(6e)

(Protected forest area constraints)

where $A{P}_k^{\pm }$ = area of tree species k for ecological forests in the project scheme (ha).

$${\displaystyle \sum _{k=1}^K\left\{\left({\displaystyle \sum _{m=1}^M{E}_{im}^{\pm }\cdot P{E}_{mk}}+{\displaystyle \sum _{n=1}^N{A}_{in}^{\pm }\cdot P{A}_{nk}}+{\displaystyle \sum _{j=1}^J{S}_{ij}^{\pm }\cdot P{S}_{jk}}+{\displaystyle \sum _{h=1}^H{G}_{ih}^{\pm }\cdot P{G}_{hk}}\right)\cdot L_k^{\pm }\right\}}\le L_i^{\pm },\hspace{0.17em} \forall i$$

(6f)

(Labour constraints)

where $L_k^{\pm }$ = labour coefficient per unit area of tree species k for meeting local production needs (man-day/ha); $L_i^{\pm }$ = gross labour force of land type i (man-day).

$$\left(\begin{array}{l}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{E}_{im}^{\pm }\cdot P{E}_{mk}}}}\cdot D{E}_{kt}^{\pm }+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }}}D{E}_{kt}^{\pm }\\ +{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{S}_{ij}^{\pm }\cdot P{S}_{jk}\cdot }}}D{E}_{kt}^{\pm }+{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot D{E}_{kt}^{\pm }}}}\end{array}\right)\cdot P_{DE}^{\pm }\ge E_{DEi}^{\pm },\hspace{0.17em} \forall i$$

(6g)

(Ecological and environmental constraints: benefits of windbreak and sand fixation)

where $D{E}_{kt}^{\pm }$ = amount of dust absorption by tree species k per unit area on the t^th year (tonnes/ha); $P_{DE}^{\pm }$ = dust removal costs (million CNY/tonne); $E_{DEi}^{\pm }$ = total benefit of windbreak and sand fixation in land type i during the project implementation period (million CNY).

$$\left(\begin{array}{l}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{E}_{im}^{\pm }\cdot P{E}_{mk}}}}\cdot {\beta }_{kt}^{\pm }+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }}}{\beta }_{kt}^{\pm }\\ +{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{S}_{ij}^{\pm }\cdot P{S}_{jk}\cdot }}}{\beta }_{kt}^{\pm }+{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot {\beta }_{kt}^{\pm }}}}\end{array}\right)\cdot Q_i^{\pm }\cdot P_w^{\pm }\cdot \left(1-{\alpha }_1^{\pm }\right)\cdot \left(1-{\alpha }_2^{\pm }\right)\ge E_{SWi}^{\pm },\hspace{0.17em} \forall i$$

(6h)

(Ecological and environmental constraints: benefits of water purification)

where $Q_i^{\pm }$ = average annual rainfall of land type i (mm); $P_w^{\pm }$ = treatment cost for urban domestic water (million CNY/m³); ${\beta }_{kt}^{\pm }$ = percentage of rainfall volume stored by tree species k during the t^th year (%); ${\alpha }_1^{\pm }$ = percentage of non-forest belt effect (%); ${\alpha }_2^{\pm }$ = physiological water consumption of forest (%); $E_{SWi}^{\pm }$ = total benefit of water purification in land type i during the project implementation period (million CNY).

$$\left(\begin{array}{l}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{E}_{im}^{\pm }\cdot P{E}_{mk}}}}\cdot \rho {\rho }_{kt}+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }}}\rho {\rho }_{kt}\\ +{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{S}_{ij}^{\pm }\cdot P{S}_{jk}\cdot }}}\rho {\rho }_{kt}+{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{t=1}^T{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot \rho {\rho }_{kt}}}}\end{array}\right)\cdot P_{RC}^{\pm }\cdot {\alpha }^{\pm }\ge E_{RCi}^{\pm },\hspace{0.17em} \forall i$$

(6i)

(Ecological and environmental constraints: benefits of microclimate regulation)

where $P_{RC}^{\pm }$ = savings in electricity consumption due to temperature regulation by new forest (million CNY/plant); $\rho {\rho }_{kt}$ = planting density of tree species k in a rotation forest on the t^th year (plant/ha); ${\alpha }^{\pm }$ = correction coefficient of temperature regulation by tree species (%); $E_{RCi}^{\pm }$ = total benefit of microclimate regulation in land type i during the project implementation period (million CNY).

$$\begin{array}{l}{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=k_w}^{K_w}{\displaystyle \sum _{t=1}^T{E}_{im}^{\pm }\cdot P{E}_{mk}}}}\cdot U_{kt}^{\pm }}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=k_w}^{K_w}{\displaystyle \sum _{t=1}^T{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }}}{U}_{kt}^{\pm }}\\ +{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=k_w}^{K_w}{\displaystyle \sum _{t=1}^T{S}_{ij}^{\pm }\cdot P{S}_{jk}\cdot }}}{U}_{kt}^{\pm }}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{k=k_w}^{K_w}{\displaystyle \sum _{t=1}^T{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot U_{kt}^{\pm }}}}}\ge U_0^{\pm }\end{array}$$

(6j)

(Benefits to the social economy)

where $k_w$ = the index of tree species for economic forest in the project scheme; $U_{kt}^{\pm }$ = economic benefit of tree species k per unit area in an economic forest the t^th year (million CNY/ha); $U_0^{\pm }$ = total economic benefit during the project implementation period (million CNY).

$$\left(\begin{array}{l}{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=k_v}^{K_v}{\displaystyle \sum _{t=1}^T{E}_{im}^{\pm }\cdot P{E}_{mk}}}}\cdot O_{kt}^{\pm }}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=k_v}^{K_v}{\displaystyle \sum _{t=1}^T{A}_{in}^{\pm }\cdot P{A}_{nk}\cdot }}}{O}_{kt}^{\pm }}+\\ {\displaystyle \sum _{i=1}^I{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=k_v}^{K_v}{\displaystyle \sum _{t=1}^T{S}_{ij}^{\pm }\cdot P{S}_{jk}\cdot }}}{O}_{kt}^{\pm }}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{k=k_v}^{K_v}{\displaystyle \sum _{t=1}^T{G}_{ih}^{\pm }\cdot P{G}_{hk}\cdot O_{kt}^{\pm }}}}}\end{array}\right)\cdot q^{\pm }/TI\ge O_0^{\pm }$$

(6k)

(Timber yield constraints)

where $k_r$ = the index of tree species for ecological forest in the project scheme; $O_{kt}^{\pm }$ = stand volume of tree species k per unit area in each year of a felling period (m³/ha); $TI$ = total timber yield (a); $q^{\pm }$ = timber production rate (%); $O_0^{\pm }$ = annual timber yield during a felling period (million CNY/yr).

3. Case Study

3.1. Overview of the Study System

A CDM forestry project in Inner Mongolia will be used to demonstrate the applicability of the ILP optimization model. In this new CDM project, carbon sequestration is to be increased by afforestation with the help of loans from the European Investment Bank. This project has taken both ecological benefits and farmers’ incomes into account. The project area includes four types of geographic regions. The map of the project area is shown in Figure 1. The total area of afforestation is 31,805.7 ha, all of which are suitable for woodland. The implementation period is 25 years, and the total investment in afforestation is 21,033 million CNY. According to the recommended scheme, the new forest will absorb and fix 3,862,195.07 tonnes of CO₂ over the implementation period. The main goods produced by the new forests are timber, Cistanche deserticola and Caragana korshinskii Kom products. The expected benefits of windbreak and sand fixation, water purification, microclimate regulation, and socio-economic objectives are 2,628, 2,354, 7,459 and 7,910 million CNY, respectively.

Figure 1. Map of the CDM project (color in green) in Inner Mongolia Autonomous Region, China.

Figure 1. Map of the CDM project (color in green) in Inner Mongolia Autonomous Region, China.

Afforestation sites and tree species will be selected for this project by applying different approaches to different land types and trees. To maintain the health of the forest and satisfy the project objectives, the objective is to select the suitable amount of mixed forest while maintaining a character at each site consistent with the selected species. Each tree species will be determined considering local resources. The selected tree species and the areas of the recommended scheme are shown in Table 1.

Table 1. Originally recommended scheme for the CDM forestry project.

Afforestation sites	Afforestation mode		Afforestation scale (ha)
	Tree species	Planting pattern
Land type 1	Pinus tabulaeformis Carr, Larix gmellini (Rupr.), A. sibirica (L) Lam	Mixed forest	737.1
	Pinus tabulaeformis Carr, Larix gmellini (Rupr.), Prunus tomentosa Thunb	Mixed forest	550
	Mongolian Scotch Pine, Salix matsudana Kodidz	Mixed forest	696
	Flos Caryophyllata, Rosa xanthina Lindl	Mixed forest	20.4
	Robinia pseudoacacia cv. idaho	Pure forest	14
	Salix matsudana Kodidz	Pure forest	14
	Poplar	Pure forest	57
	Mongolian Scotch Pine	Pure forest	100
	Sophora japonica Linn	Pure forest	26
	Catalpa ovata Don	Pure forest	20
Land type 2	Caragana Korsgunskii Kom, Elaegnus angustifolius	Mixed forest	6,000
	Caragana Korsgunskii Kom, Ulmus pumila L.	Mixed forest	6,000
	Poplar	Pure forest	666
	Caragana Korsgunskii Kom	Pure forest	3,275
Land type 3	Poplar, Ulmus pumila L.	Mixed forest	4,930.6
	Poplar, Mongolian Scotch Pine	Mixed forest	4,930.6
Land type 4	Elaegnus angustifolius, Tamarix chinensis Lour	Mixed forest	1,000
	Haloxylon ammodendron	Pure forest	2,769

Total carbon sink capacity (tonnes): 3,862,195.07.

Table 1. Originally recommended scheme for the CDM forestry project.

Afforestation sites	Afforestation mode		Afforestation scale (ha)
	Tree species	Planting pattern
Land type 1	Pinus tabulaeformis Carr, Larix gmellini (Rupr.), A. sibirica (L) Lam	Mixed forest	737.1
	Pinus tabulaeformis Carr, Larix gmellini (Rupr.), Prunus tomentosa Thunb	Mixed forest	550
	Mongolian Scotch Pine, Salix matsudana Kodidz	Mixed forest	696
	Flos Caryophyllata, Rosa xanthina Lindl	Mixed forest	20.4
	Robinia pseudoacacia cv. idaho	Pure forest	14
	Salix matsudana Kodidz	Pure forest	14
	Poplar	Pure forest	57
	Mongolian Scotch Pine	Pure forest	100
	Sophora japonica Linn	Pure forest	26
	Catalpa ovata Don	Pure forest	20
Land type 2	Caragana Korsgunskii Kom, Elaegnus angustifolius	Mixed forest	6,000
	Caragana Korsgunskii Kom, Ulmus pumila L.	Mixed forest	6,000
	Poplar	Pure forest	666
	Caragana Korsgunskii Kom	Pure forest	3,275
Land type 3	Poplar, Ulmus pumila L.	Mixed forest	4,930.6
	Poplar, Mongolian Scotch Pine	Mixed forest	4,930.6
Land type 4	Elaegnus angustifolius, Tamarix chinensis Lour	Mixed forest	1,000
	Haloxylon ammodendron	Pure forest	2,769

Total carbon sink capacity (tonnes): 3,862,195.07.

Carbon sequestration cannot be the sole criterion for tree species selection, because this ability also depends on the characteristics of the afforestation site. Consequently, the available tree species must be effectively allotted to maximize the carbon sink capacity of the CDM forestry projects subject to suitable environmental and socio-economic standards. Table 2 presents the carbon sink capacities per unit area of the tree species in different stand ages. The calculation method and main database were obtained from Fang et al. [22]. Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 describe the CDM forestry project, target demands, economic data and technical data.

Table 2. Carbon sink capacities per unit area of each tree species for different stand ages.

Tree species	1–10 years old	10–20 years old	Above 20 years old
Pinus tabulaeformis Carr	[1.674, 2.046]	[1.674, 2.046]	[1.881, 3.333]
Mongolian Scotch Pine	[3.42, 4.18]	[3.42, 4.18]	[3.609, 5.346]
Larix gmellini (Rupr.) Pupr	[4.383, 5.357]	[4.383, 5.357]	[4.410, 5.522]
Poplar	[6.39, 7.81]	[6.507, 7.93]	[6.507, 7.953]
Ulmus pumila L. (i)	[3.825, 4.675]	[3.825, 4.675]	[3.825, 4.675]
Ulmus pumila L. (ii)	[6.867, 8.393]	[6.867, 8.393]	[6.867, 8.393]
Salix matsudana Kodidz	[5.463, 6.677]	[5.463, 6.677]	[5.733, 7.876]
Sophora japonica Linn	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Catalpa ovata Don	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Robinia pseudoacacia cv. idaho	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Elaegnus angustifolius	[5.571, 6.389]	[5.571, 6.389]	[5.571, 6.389]
Haloxylon ammodendron	[6.489, 7.931]	[6.489, 7.931]	[6.489, 7.931]
Caragana Korsgunskii Kom	[15.048, 18.392]	[15.048, 18.392]	[15.048, 18.392]
Tamarix chinensis Lour	[3.366, 4.114]	[3.366, 4.114]	[3.366, 4.114]
Flos Caryophyllata	[3.366, 4.114]	[3.366, 4.114]	[3.366, 4.114]
Rosa xanthina Lindl	[3.366, 4.141]	[3.366, 4.114]	[3.366, 4.114]
A. sibirica (L.) Lam	[5.571, 6.809]	[5.571, 6.809]	[5.571, 6.809]
Prunus tomentosa Thunb	[3.336, 4.114]	[3.336, 4.114]	[3.336, 4.114]

Table 2. Carbon sink capacities per unit area of each tree species for different stand ages.

Tree species	1–10 years old	10–20 years old	Above 20 years old
Pinus tabulaeformis Carr	[1.674, 2.046]	[1.674, 2.046]	[1.881, 3.333]
Mongolian Scotch Pine	[3.42, 4.18]	[3.42, 4.18]	[3.609, 5.346]
Larix gmellini (Rupr.) Pupr	[4.383, 5.357]	[4.383, 5.357]	[4.410, 5.522]
Poplar	[6.39, 7.81]	[6.507, 7.93]	[6.507, 7.953]
Ulmus pumila L. (i)	[3.825, 4.675]	[3.825, 4.675]	[3.825, 4.675]
Ulmus pumila L. (ii)	[6.867, 8.393]	[6.867, 8.393]	[6.867, 8.393]
Salix matsudana Kodidz	[5.463, 6.677]	[5.463, 6.677]	[5.733, 7.876]
Sophora japonica Linn	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Catalpa ovata Don	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Robinia pseudoacacia cv. idaho	[5.867, 8.393]	[5.867, 8.393]	[5.867, 8.393]
Elaegnus angustifolius	[5.571, 6.389]	[5.571, 6.389]	[5.571, 6.389]
Haloxylon ammodendron	[6.489, 7.931]	[6.489, 7.931]	[6.489, 7.931]
Caragana Korsgunskii Kom	[15.048, 18.392]	[15.048, 18.392]	[15.048, 18.392]
Tamarix chinensis Lour	[3.366, 4.114]	[3.366, 4.114]	[3.366, 4.114]
Flos Caryophyllata	[3.366, 4.114]	[3.366, 4.114]	[3.366, 4.114]
Rosa xanthina Lindl	[3.366, 4.141]	[3.366, 4.114]	[3.366, 4.114]
A. sibirica (L.) Lam	[5.571, 6.809]	[5.571, 6.809]	[5.571, 6.809]
Prunus tomentosa Thunb	[3.336, 4.114]	[3.336, 4.114]	[3.336, 4.114]

Table 3. Technical parameters of the CDM forestry project.

Tree species	Planting density (plant/ha)	Labour coefficient (man-day/ha)	Amount of dust absorption (tonne/ha)
Pinus tabulaeformis Carr	1250	[84.5, 86]	[10.8, 11.7]
Mongolian Scotch Pine	1250	[84.5, 86]	[11.7, 12.6]
Larix gmellini (Rupr.) Pupr.	1250	[84.5, 86]	[10.8, 11.7]
Poplar	1250	[84.5, 86]	[10.8, 11.7]
Ulmus pumila L.	1666	[84.5, 86]	[11.7, 12.6]
Salix matsudana Kodidz	1666	[84.5, 86]	[11.7, 12.6]
Sophora japonica Linn	1250	[84.5, 86]	[10.8, 11.7]
Catalpa ovata Don	1250	[84.5, 86]	[10.8, 11.7]
Robinia pseudoacacia cv. idaho	1250	[84.5, 86]	[10.8, 11.7]
Elaegnus angustifolius	1666	[84.5, 86]	[10.8, 12.6]
Haloxylon ammodendron	1157	[41.5, 43]	[10.8, 12.6]
Caragana Korsgunskii Kom	2000	[41.5, 43]	[10.8, 12.6]
Tamarix chinensis Lour	1250	[41.5, 43]	[10.8, 12.6]
Flos Caryophyllata	1666	[41.5, 43]	[10.8, 11.7]
Rosa xanthina Lindl	1666	[41.5, 43]	[10.8, 11.7]
A. sibirica (L) Lam	1250	[41.5, 43]	[10.8, 12.6]
Prunus tomentosa Thunb	1250	[41.5, 43]	[10.8, 11.7]

Table 3. Technical parameters of the CDM forestry project.

Tree species	Planting density (plant/ha)	Labour coefficient (man-day/ha)	Amount of dust absorption (tonne/ha)
Pinus tabulaeformis Carr	1250	[84.5, 86]	[10.8, 11.7]
Mongolian Scotch Pine	1250	[84.5, 86]	[11.7, 12.6]
Larix gmellini (Rupr.) Pupr.	1250	[84.5, 86]	[10.8, 11.7]
Poplar	1250	[84.5, 86]	[10.8, 11.7]
Ulmus pumila L.	1666	[84.5, 86]	[11.7, 12.6]
Salix matsudana Kodidz	1666	[84.5, 86]	[11.7, 12.6]
Sophora japonica Linn	1250	[84.5, 86]	[10.8, 11.7]
Catalpa ovata Don	1250	[84.5, 86]	[10.8, 11.7]
Robinia pseudoacacia cv. idaho	1250	[84.5, 86]	[10.8, 11.7]
Elaegnus angustifolius	1666	[84.5, 86]	[10.8, 12.6]
Haloxylon ammodendron	1157	[41.5, 43]	[10.8, 12.6]
Caragana Korsgunskii Kom	2000	[41.5, 43]	[10.8, 12.6]
Tamarix chinensis Lour	1250	[41.5, 43]	[10.8, 12.6]
Flos Caryophyllata	1666	[41.5, 43]	[10.8, 11.7]
Rosa xanthina Lindl	1666	[41.5, 43]	[10.8, 11.7]
A. sibirica (L) Lam	1250	[41.5, 43]	[10.8, 12.6]
Prunus tomentosa Thunb	1250	[41.5, 43]	[10.8, 11.7]

Table 4. Technical parameters of the timber forest during a felling period.

Timber forest	Stand volume (m3/ha)	Timber-produced rate (%)
Poplar	[165, 181.5]	[62.5, 65]

Table 4. Technical parameters of the timber forest during a felling period.

Timber forest	Stand volume (m3/ha)	Timber-produced rate (%)
Poplar	[165, 181.5]	[62.5, 65]

Table 5. Other technical parameters of the CDM forestry project.

Parameters	Numerical value
Correction coefficient of temperature regulation (%)	[5.5, 6]
Percentage of rainfall stock volume (%)	[72.5, 75]
Percentage of non-forest belt effect (%)	[50, 50.1]
Physiological water consumption of forest (%)	[14.9, 15]

Table 5. Other technical parameters of the CDM forestry project.

Parameters	Numerical value
Correction coefficient of temperature regulation (%)	[5.5, 6]
Percentage of rainfall stock volume (%)	[72.5, 75]
Percentage of non-forest belt effect (%)	[50, 50.1]
Physiological water consumption of forest (%)	[14.9, 15]

Table 6. Mixed modes and ratios.

Mixed types	Mixed modes	Mixed ratios
Coniferous and broadleaved mixed forest	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr.	3:1
	Mongolian Scotch Pine, Larix gmellini (Rupr.) Pupr.	1:10
	Poplar, Pinus tabulaeformis Carr	8:2
	Ulmus pumila L., Mongolian Scotch Pine	8:2
	Mongolian Scotch Pine, Salix matsudana Kodidz	8:2
	Poplar, Mongolian Scotch Pine	8:2
Broadleaved mixed forest	Poplar, Ulmus pumila L.	8:2
	Poplar, Elaegnus angustifolius	8:2
	Ulmus pumila L., Elaegnus angustifolius	8:2
Arbers mixed forest	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr., A. sibirica (L) Lam	4:4:2
	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr., Prunus tomentosa Thunb	4:4:2
	Pinus tabulaeformis Carr, A.sibirica (L) Lam	2:3
	Pinus tabulaeformis Carr, Caragana Korsgunskii Kom	1:1
	Mongolian Scotch Pine, A.sibirica (L) Lam	2:3
	Mongolian Scotch Pine, Caragana Korsgunskii Kom	1:1
	Elaegnus angustifolius, Tamarix chinensis Lour	1:1
	Ulmus pumila L., Caragana Korsgunskii Kom	2:1
	Elaegnus angustifolius, Haloxylon ammodendron	1:1
Shrubs mixed forest	Flos Caryophyllata, Rosa xanthina Lindl	1:1
	Elaegnus angustifolius, Caragana Korsgunskii Kom	2:1

Table 6. Mixed modes and ratios.

Mixed types	Mixed modes	Mixed ratios
Coniferous and broadleaved mixed forest	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr.	3:1
	Mongolian Scotch Pine, Larix gmellini (Rupr.) Pupr.	1:10
	Poplar, Pinus tabulaeformis Carr	8:2
	Ulmus pumila L., Mongolian Scotch Pine	8:2
	Mongolian Scotch Pine, Salix matsudana Kodidz	8:2
	Poplar, Mongolian Scotch Pine	8:2
Broadleaved mixed forest	Poplar, Ulmus pumila L.	8:2
	Poplar, Elaegnus angustifolius	8:2
	Ulmus pumila L., Elaegnus angustifolius	8:2
Arbers mixed forest	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr., A. sibirica (L) Lam	4:4:2
	Pinus tabulaeformis Carr, Larix gmellini (Rupr.) Pupr., Prunus tomentosa Thunb	4:4:2
	Pinus tabulaeformis Carr, A.sibirica (L) Lam	2:3
	Pinus tabulaeformis Carr, Caragana Korsgunskii Kom	1:1
	Mongolian Scotch Pine, A.sibirica (L) Lam	2:3
	Mongolian Scotch Pine, Caragana Korsgunskii Kom	1:1
	Elaegnus angustifolius, Tamarix chinensis Lour	1:1
	Ulmus pumila L., Caragana Korsgunskii Kom	2:1
	Elaegnus angustifolius, Haloxylon ammodendron	1:1
Shrubs mixed forest	Flos Caryophyllata, Rosa xanthina Lindl	1:1
	Elaegnus angustifolius, Caragana Korsgunskii Kom	2:1

Table 7. Economic parameters of the CDM forestry project.

Economic parameters	Unit	Numerical value
Pinus tabulaeformis Carr seedling	CNY/plant	[1.0, 1.1]
Larix gmellini (Rupr.) Pupr. seedling	CNY/plant	[0.8, 0.88]
Mongolian Scotch Pine seedling	CNY/plant	[3.25, 3.5]
Poplar seedling	CNY/plant	[6.0, 6.5]
Ulmus pumila L. seedling	CNY/plant	[2.5, 2.75]
Salix matsudana Kodidz seedling	CNY/plant	[6.0, 6.6]
Caragana Korsgunskii Kom seedling	CNY/plant	[0.15, 0.18]
Other tree species seedling	CNY/plant	[0.4, 0.5]
Caragana Korshinskii Kom products	CNY/kg	[0.40, 0.56]
Cistanche Deserticola products	CNY/kg	[20, 22]
Wastewater treatment cost	CNY/m3	[0.825, 0.850]
Electricity price converted by temperature regulation	CNY/plant	[30.00, 30.01]
Dust removal costs	CNY/tonne	[82.84, 85]

Table 7. Economic parameters of the CDM forestry project.

Economic parameters	Unit	Numerical value
Pinus tabulaeformis Carr seedling	CNY/plant	[1.0, 1.1]
Larix gmellini (Rupr.) Pupr. seedling	CNY/plant	[0.8, 0.88]
Mongolian Scotch Pine seedling	CNY/plant	[3.25, 3.5]
Poplar seedling	CNY/plant	[6.0, 6.5]
Ulmus pumila L. seedling	CNY/plant	[2.5, 2.75]
Salix matsudana Kodidz seedling	CNY/plant	[6.0, 6.6]
Caragana Korsgunskii Kom seedling	CNY/plant	[0.15, 0.18]
Other tree species seedling	CNY/plant	[0.4, 0.5]
Caragana Korshinskii Kom products	CNY/kg	[0.40, 0.56]
Cistanche Deserticola products	CNY/kg	[20, 22]
Wastewater treatment cost	CNY/m3	[0.825, 0.850]
Electricity price converted by temperature regulation	CNY/plant	[30.00, 30.01]
Dust removal costs	CNY/tonne	[82.84, 85]

Table 8. Average annual rainfall at the CDM forestation sites and the gross labour force available in each area.

Afforestation sites	Average annual rainfall (mm)	Gross labour force (million man-days)
Land type 1	[434.9, 535.6]	[17.35, 17.69]
Land type 2	[300, 375]	[101.03, 103.42]
Land type 3	[400, 450]	[81.85, 83.33]
Land type 4	102.9	[17.23, 17.79]

Table 8. Average annual rainfall at the CDM forestation sites and the gross labour force available in each area.

Afforestation sites	Average annual rainfall (mm)	Gross labour force (million man-days)
Land type 1	[434.9, 535.6]	[17.35, 17.69]
Land type 2	[300, 375]	[101.03, 103.42]
Land type 3	[400, 450]	[81.85, 83.33]
Land type 4	102.9	[17.23, 17.79]

The problem under consideration is how to optimally allocate tree species in the presence of uncertainty. Since these uncertainties are expressed as intervals, the ILP method is considered a feasible approach for this type of planning problem.

3.2. Results and Discussion

Table 9 presents the solutions obtained by solving the model (6). Note that the solutions obtained for the objective function value and some decision variables are themselves intervals. Generally, solutions presented as intervals correspond to the upper and lower bounds of the objective function. Such results emphasise that the decisions are sensitive to uncertain inputs. In contrast, variables with deterministic solutions are not sensitive to the input uncertainties. Thus, alternative schemes can be generated by adjusting the uncertain solutions within their intervals according to the requirements.

Table 9. Optimal solutions obtained from the model.

Afforestation sites (i)	Afforestation mode		Afforestation scale (ha)
	Tree species	Planting pattern	Symbol	Value
Land type 1	Mongolian Scotch Pine, Salix matsudana Kodidz	Mixed forest	$S_{12}^{\pm }$	1,292.49
	Mongolian Scotch Pine, A. sibirica (L) Lam	Mixed forest	$S_{19}^{\pm }$	245.7
	Mongolian Scotch Pine, Poplar	Mixed forest	$G_{12}^{\pm }$	[0, 71.25]
	Mongolian Scotch Pine	Pure forest	$E_{12}^{\pm }$	494.66
	Haloxylon ammodendron	Pure forest	$E_{1\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	130.4
Land type 2	Caragana Korsgunskii Kom, Elaegnus angustifolius	Mixed forest	$G_{28}^{\pm }$	7,476.57
	Caragana Korsgunskii Kom, Ulmus pumila L.	Mixed forest	$G_{27}^{\pm }$	4,271.07
	Poplar, Ulmus pumila L.	Mixed forest	$G_{24}^{\pm }$	[480.88, 832.5]
	Mongolian Scotch Pine, Caragana Korsgunskii Kom	Mixed forest	$G_{26}^{\pm }$	3.48
	Caragana Korsgunskii Kom	Pure forest	$A_{22}^{\pm }$	3,357.38
Land type 3	Poplar, Ulmus pumila L.	Mixed forest	$G_{34}^{\pm }$	[9,669.82, 9,861.2]
Land type 4	Elaegnus angustifolius, Haloxylon ammodendron	Mixed forest	$S_{4\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	[973.5, 1,000]
	Tamarix chinensis Lour	Pure forest	$E_{4\hspace{0.17em}\hspace{0.17em}11}^{\pm }$	[400, 500]
	Haloxylon ammodendron	Pure forest	$E_{4\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	2,269

Total carbon sink capacity (tonnes):$f_{opt}^{\pm }=$ [3863337, 4747957].

Table 9. Optimal solutions obtained from the model.

Afforestation sites (i)	Afforestation mode		Afforestation scale (ha)
	Tree species	Planting pattern	Symbol	Value
Land type 1	Mongolian Scotch Pine, Salix matsudana Kodidz	Mixed forest	$S_{12}^{\pm }$	1,292.49
	Mongolian Scotch Pine, A. sibirica (L) Lam	Mixed forest	$S_{19}^{\pm }$	245.7
	Mongolian Scotch Pine, Poplar	Mixed forest	$G_{12}^{\pm }$	[0, 71.25]
	Mongolian Scotch Pine	Pure forest	$E_{12}^{\pm }$	494.66
	Haloxylon ammodendron	Pure forest	$E_{1\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	130.4
Land type 2	Caragana Korsgunskii Kom, Elaegnus angustifolius	Mixed forest	$G_{28}^{\pm }$	7,476.57
	Caragana Korsgunskii Kom, Ulmus pumila L.	Mixed forest	$G_{27}^{\pm }$	4,271.07
	Poplar, Ulmus pumila L.	Mixed forest	$G_{24}^{\pm }$	[480.88, 832.5]
	Mongolian Scotch Pine, Caragana Korsgunskii Kom	Mixed forest	$G_{26}^{\pm }$	3.48
	Caragana Korsgunskii Kom	Pure forest	$A_{22}^{\pm }$	3,357.38
Land type 3	Poplar, Ulmus pumila L.	Mixed forest	$G_{34}^{\pm }$	[9,669.82, 9,861.2]
Land type 4	Elaegnus angustifolius, Haloxylon ammodendron	Mixed forest	$S_{4\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	[973.5, 1,000]
	Tamarix chinensis Lour	Pure forest	$E_{4\hspace{0.17em}\hspace{0.17em}11}^{\pm }$	[400, 500]
	Haloxylon ammodendron	Pure forest	$E_{4\hspace{0.17em}\hspace{0.17em}10}^{\pm }$	2,269

Total carbon sink capacity (tonnes):$f_{opt}^{\pm }=$ [3863337, 4747957].

For example, the solutions of $G_{24}^{\pm }$ under the given targets reflect variations caused by the uncertain inputs. The upper bound of $G_{24}^{\pm }$ (i.e., $G_{24}^{+}$) corresponds to a greater benefit; the lower bound of $G_{24}^{\pm }$ (i.e., $G_{24}^{-}$) corresponds to a lesser benefit.

The solution of$f_{opt}^{\pm }=$ [3863337, 4747957] (tonnes) provides the range of carbon sink capacities possible under the optimal tree species allocation pattern. As the actual value of each variable or parameter may lie anywhere within its interval, the carbon sink capacities obtained from the CDM forestry project may change between$f_{opt}^{-}$ and$f_{opt}^{+}$ as the system variables change, either due to increased knowledge or over time. Planning for the lower-bound of the objective-function value will lead to a lower project carbon sinks. Conversely, planning with a higher system benefit will correspond to a higher project carbon sinks. Therefore, the optimal solutions obtained from the model are flexible in reflecting possible system condition variations caused by the existence of input uncertainties. Compared with the project’s original scheme, an extra [1,142, 885,762] tonnes of carbon dioxide will be absorbed and fixed by the optimal solutions obtained from the model. Figure 2 compares the total carbon sink capacities of the two models over the 25-year implementation period. In relative terms, the increased capacity lies in the range [0.03%, 22.93%]. The benefits of windbreak and sand fixation, water purification, microclimate regulation and socio-economic objectives obtained from the optimal solutions are [2,655, 3,067], [2,357, 3,064], [7,516, 8,382] and [8,515, 10,133] million CNY respectively. Figure 3 compares these benefits under the two schemes.

This optimization model of carbon sinks based on ILP has some advantages over other methods in terms of data availability and computational requirements. It is also particularly suitable to problems involving the carbon sink capacities of tree species, where the system information may be limited to estimates of upper and lower bounds.

Figure 2. Total carbon sink capacities of the CDM project’s original recommended scheme and the optimal solutions over a 25-year period.

Figure 2. Total carbon sink capacities of the CDM project’s original recommended scheme and the optimal solutions over a 25-year period.

Figure 3. Other benefits of the recommended scheme and the optimal solutions.

Figure 3. Other benefits of the recommended scheme and the optimal solutions.

The results of our case study show that decision-makers can not only identify the most desirable tree species and allocation plan, but also maximise the carbon sink capacity of the project while obtaining additional benefits for the ecological environment and socio-economic objectives. However, this optimisation model of carbon sinks may have limitations when the study system is very large and complicated. For example, the optimisation mode would be very complicated when the kinds of species are numerous. In such a situation, more complex models should be developed.

4. Conclusions

A mathematical programming method, interval linear programming, has been introduced for the first time in the context of CDM forestry projects. Based on ILP, we have developed a model for optimizing the capacity of carbon sinks in parallel with other decision criteria under uncertain input parameters. The results of the case study confirm that our method finds reasonable solutions for the problem of tree species allocation. Solution variables in interval form are useful for decision-makers, enabling them to form allocation plans that maximize carbon sinks while reaping diverse benefits in other criteria from the CDM forestry projects.

Decision alternatives can be generated by adjusting the solution variables within their intervals according to projects requirements. Compared with our case study’s originally recommended scheme, between 1,142 and 885,762 extra tonnes of carbon dioxide are absorbed and fixed by our optimized solutions over the 25-year implementation period. Furthermore, these solutions offer increased benefits in other criteria: windbreak and sand fixation, water purification, microclimate regulation and socio-economic objectives.

In general, although this study is the first attempt to plan a CDM forestry projects based on an ILP optimization model of carbon sinks, our results show that the method is applicable. More complex models would improve its range of applicability.