IMPLEMENTING MINIMUM AREA HARVESTING BLOCKS IN AN OPTIMIZED FOREST PLANNING MODEL¹

IMPLEMENTANDO BLOCOS DE COLHEITA FLORESTAL COM ÁREA MÍNIMA EM UM MODELO DE PLANEJAMENTO FLORESTAL OTIMIZADO

ABSTRACT

The spatial relationships between forest stands scheduled for harvesting are crucial for forest planning, once they affect directly environmental impacts and logistic aspects related to this operation. In this context, the objective of this study was to perform the optimized forest planning of a Pinus taeda, Pinus elliottii and Eucalyptus dunnii forest, composed by 236 stands, including minimum area constraints. Two approaches were applied for clustering harvesting activities. The first approach considered the inclusion of a set of constraints so-called ring inequalities, proposed by Carvajal et al. (2013). The second approach was based on the formulation proposed by Rebain and McDill (2003) for the creation of old-growth forest areas. Both formulations were capable to generate harvesting blocks with minimum area requirements, causing a reduction up to 5.1% in the objective function for the most restrictive scenario. We conclude that the formulation based on Carvajal et al. (2013) is the best alternative when only minimum area constraints are considered, due to the inferior number of constraints and superior performance in terms of the NPV generated.

Keywords:
Minimum area constraints; Forest planning; Forest harvest

RESUMO

As relações espaciais entre os talhões designados para a colheita são cruciais para o planejamento florestal, uma vez que estas relações afetam diretamente os impactos abientais decorrentes da colheita, bem como os aspectos logísticos destas operações. Nexte contexto, o objetivo deste trabalho foi realizar o planejamento florestal otimizado de uma floresta formada por 236 talhões de Pinus taeda, Pinus elliottii e Eucalyptus dunnii incluindo restrições para formar blocos de colheita respeitando um limite mínimo de área. Duas abordagens foram utilizadas para a formação dos blocos de colheita. A primeira abordagem considerou a aplicação de restrições denominadas ring inequalities, propostas por Carvajal et al. (2013). A segunda abordagem foi baseada na modelagem de Rebain e McDill (2003), para a formação de áreas de floresta madura. A segunda abordagem utilizada. Ambas as abordagens foram capazes de gerar blocos de colheita respeitando as áreas mínimas determinadas, causando uma redução de 5.1% no cenário mais restritivo. A formulação baseda em Carvajal et al. (2013) se mostrou a melhor alternativa quando apenas restrições de área mínima são consideradas, devido ao número inferior de restrições e performance superior em termos do VPL gerado.

Palavras-chave:
Restrições de área mínima; Planejamento florestal; Colheita florestal

1 INTRODUCTION

The inclusion of spatial relationships has become an essential aspect of forest planning during the past decades, especially the issues related to landscape patterns resulting from forest harvesting activities (Constantino et al., 2008). Optimization models aiming to include these spatial relationships into forest planning have been proposed in the literature, notably through adjacency constraints. The focus of spatial planning problems has been the inclusion of maximum clear-cut area constraints and developing efficient solution methods for these problems (Öhman and Lämas, 2003).

Problems considering the connectivity of forest areas, however, have been gaining importance in forest planning, such as the creation of ecological corridors and continuous core areas, aiming to ensure the maintenance of biodiversity and other ecosystem services (e.g. water protection and scenic beauty). Numerous studies addressed problems related to the creation of ecological corridors and core areas for biodiversity ( Rebain and Mcdill, 2003; Tóth et al., 2006; Wei and Hoganson, 2007;Carvajal et al., 2013).

While the creation of core areas aims to set contiguous stands aside from production, other important class of connectivity problems targets at clustering harvesting areas. These models also aim at protecting biodiversity, once spreading harvesting activities may reduce the habitat availability for interior forest species and enhance negative edge effects (Gustafson, 1996). Different approaches were proposed in the literature for dealing with clustering of harvesting activities, such as dynamic zoning of harvesting areas and multi-objective models for aggregating stands scheduled for harvesting (Gustafson, 1998; Öhman and Lämas, 2003; Öhman and Eriksson, 2010; Smaltschinski et al., 2012).

Problems involving clustering of harvesting areas and creation of core areas are closely related, allowing the use of similar approaches for its solution. Nevertheless, in the case of clustering harvesting activities the problem can become more complex, as different contiguous harvesting blocks must be formed at each period. On the other hand, core areas are usually maintained during the whole planning horizon (PH). Therewith, the clustering of harvesting areas may demand more processing time.

The objective of this study was to include adjacency constraints into an optimized forest planning model, aiming to maintain interior forest habitat and reduce edge effects. To this end, we applied two approaches, proposed initially for the creation of core areas, generating Integer Linear Programming models. The first approach was based on the dynamic connectivity model proposed by Carvajal et al. (2013) and the second was based in the old-growth forest formulation from Rebain and McDill (2003). We evaluated the efficiency of both formulations for creating harvesting blocks and their impacts on the objective function.

2. MATERIAL AND METHODS

2.1. Study area

Our research area was a planted forest with Pinus taeda, Pinus elliottii and Eucalyptus dunnii, located in the municipalities of Bituruna and General Carneiro, Brazil, between the coordinates 26º13'58.31" and 26º22'5.634" S and 51º34'14.6 and 51º30'26.14" W. The climate according to KÖppen classification is Cfb: mesotermic humid subtropical, with mild summer and occurrence of frequent frosts during winter.

The area was composed by 236 stands, with 18 Eucalyptus dunnii stands, 21 Pinus elliottii stands and 197 Pinus taeda stands, with a total area of 2365.8 ha, average stand area of 10.3 ha and the average number of adjacencies per stand was 3 (Figure 1).

Figure 1
Study area. The figure shows the map of the study area, with the productive stands highlighted.
Figura 1
Área de estudo. A figura mostra o mapa da área de estudo com os talhões produtivos em destaque.

2.2. Optimization model

The first formulation applied for ensuring the minimum contiguous harvesting area was the approach proposed by Carvajal et al. (2013), so-called ring inequalities, described in Eq. (1):

(1) $\sum _{\mathit{u}\in \partial \mathit{C}}{\mathit{x}}_{\mathit{uj}}\ge \sum _{\mathit{i}∊\mathit{C}}{\mathit{x}}_{\mathit{ij}}-{\mathit{n}}_{\mathit{c}}\forall \mathit{C}∊\mathit{C}\left(\mathit{A}\min \right),\forall \mathit{j}$

Where: x_ij: binary values that take value 1 case stand i is harvested in year j or 0 otherwise; n_c: number of stands in set C minus one; : block of stands that violate the minimum area requirement; C(A_min): all sets of stand blocks with area inferior to the established limit; ∂C: set of stands in the neighborhood of C.

Eq. (1) requires that in case a block of stands is selected for harvesting in a defined year of the planning horizon, and the sum of their areas is inferior to the established minimum area, at least one stand adjacent to the considered block and not part of the block is selected for harvesting in the same year. This constraint is applied for all sets of stands that violate the minimum area limit. Therewith, the optimized forest planning model considering minimum harvesting block areas can be defined as follows:

(2) $\mathit{Max}=\sum _{\mathit{i}=1}^{\mathit{N}}\sum _{\mathit{j}=1}^{\mathit{HP}}{\mathit{c}}_{\mathit{ij}}{\mathit{x}}_{\mathit{ij}}$

(3) $\sum _{\mathit{j}=1}^{\mathit{HP}}{\mathit{x}}_{\mathit{ij}}=1\forall \mathit{i}$

(4) $\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{v}}_{\mathit{ij}}{\mathit{x}}_{\mathit{ij}}\ge 0.9\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{v}}_{\mathit{i}1}{\mathit{x}}_{\mathit{i}1}\mathit{j}=2\mathrm{,...}\mathit{HP}$

(5) $\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{v}}_{\mathit{ij}}{\mathit{x}}_{\mathit{ij}}\le 1.1\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{v}}_{\mathit{i}1}{\mathit{x}}_{\mathit{i}1}\mathit{j}=2\mathrm{,...}\mathit{HP}$

(6) $\sum _{\mathit{u}∊\partial \mathit{C}}{\mathit{x}}_{\mathit{uj}}\ge \sum _{\mathit{i}∊\mathit{C}}{\mathit{x}}_{\mathit{ij}}-{\mathit{n}}_{\mathit{c}}\forall \mathit{C}∊\mathit{C}\left(\mathit{A}\min \right),\mathit{j}=1\mathrm{,...}5$

(7) ${\mathit{x}}_{\mathit{ij}}∊\left\{0,1\right\},\forall \mathit{i}\forall \mathit{j}$

Where: N: number of stands; NP number of years in the planning horizon; c_ij: NPV generated by stand i in case it is harvested in year j; x_ij: binary variable that takes value 1 case stand is harvested in year j or 0 otherwise; v_ij: volume generated by stand i in year j; A_min: minimum area limit; ; C: block of stands that violate the minimum area requirement; C(A_min): all sets of stand blocks with area inferior to the established limit; ∂C: set of stands in the neighborhood of C.

The objective function of the optimization model is represented in Eq. (2), maximizing the total forest total NPV, from the sum of the NPVs generated by each stand under a certain management regime. Constraint Eq. (3) imposes that stands are harvested only once during the planning horizon. Constraints Eq. (4) and Eq. (5) are the wood flow constraints, ensuring that the harvested volume during each period of the planning horizon respects the established limits. In our study, the allowed fluctuation for the harvested volume at each period of the PH was 10% in comparison to the volume harvested in the first year. Constraint (6) is the minimum area set of constraints, ensuring that the area of stands scheduled for harvesting are greater or equal the established limit. We tested different minimum harvesting area limits, creating scenarios with minimum area constraints of 30, 40, 50 and 60 ha. Constraint Eq. (7) requires that the decision variables x_ij take binary values.

The second approach, based on the model from Rebain and McDill (2003), creates harvesting blocks a priori applying a minimum and maximum block area threshold and ensures that all stands selected for harvesting in a defined year of the PH are part of at least one of these blocks:

(8) $\sum _{\mathit{i}∊{\mathit{C}}_{\mathit{k}}}{\mathit{x}}_{\mathit{ij}}\ge {\mathit{n}}_{\mathit{k}}{\mathit{B}}_{\mathit{kj}}\forall \mathit{j},\forall \mathit{k}∊\mathit{S}$

(9) $\sum _{\mathit{k}∊{\mathit{P}}_{\mathit{i}}}{\mathit{B}}_{\mathit{kj}}\ge {\mathit{x}}_{\mathit{ij}}\forall \mathit{i},\forall \mathit{j}$

Where x_ij: bynary variable that takes value 1 case stand is harvested in year or 0 otherwise; n_k: number of stands in block ; : binary variable that takes value 1 case block k is selected for harvesting in year j; C_k: set of stands in block k; S: set of all blocks of stands formed applying the minimum and maximum area limits; P_i: set of harvesting blocks containing stand i.

Constraint Eq. (8) guarantees that in case block k is selected for harvesting in a defined year of the planning horizon, all stands in this block are selected for harvesting at the same year. Constraint (9) aims to ensure that case stand i is selected for harvesting in a defined year of the planning horizon, a harvesting block containing stand i is also selected for harvesting. For the formulation based on Rebain and McDill (2003) the set of minimum area constraints Eq. (6) was replaced by the set of constraints Eq. (8) and Eq. (9), with j varying from 1 to 5, the objective function and other constraints remained unchanged. We highlight that this set of constraints allows harvesting blocks to overlap, thus generating contiguous harvesting areas greater than the limits applied for creating the blocks.

The adopted planning horizon was 16 years, applying management regimes without thinning. The wood flow constraints allowed a fluctuation of 10% for more or less, compared to the wood production in the first year of the PH. For producing harvesting blocks according to the formulation from Rebain and McDill (2003), we applied ranges of 10 and 20 ha for the harvesting block area. Hence, we obtained 13 scenarios (Table 1). The planning scenarios were generated through the software Optimber-LP and the optimization models were solved through Gurobi 5.5 in a PC with Intel® Core™ Duo CPU 2.93 GHz processor and 4Mb of memory. The optimization time limit was established as 90 minutes.

Table 1
Tested scenarios. The table shows all scenarios tested in the optimization, with its respective constraints.
Tabela 1
Cenários testados. A tabela mostra todos os cenários testados na otimização, com as respectivas restrições.

Scenario	Name	Constraints
1	baseline	2,3,4,5,7
2	Ring 30 ha and wood flow	2,3,4,5,6,7
3	Ring 40 ha and wood flow	2,3,4,5,6,7
4	Ring 50 ha and wood flow	2,3,4,5,6,7
5	Ring 60 ha and wood flow	2,3,4,5,6,7
6	McDill 30 – 40 ha and wood flow	2,3,4,5,7,8,9
7	McDill 40 – 50 ha and wood flow	2,3,4,5,7,8,9
8	McDill 50 – 60 ha and wood flow	2,3,4,5,7,8,9
9	McDill 60 – 70 ha and wood flow	2,3,4,5,7,8,9
10	McDill 30 – 50 ha and wood flow	2,3,4,5,7,8,9
11	McDill 40 – 60 ha and wood flow	2,3,4,5,7,8,9
12	McDill 50 – 70 ha and wood flow	2,3,4,5,7,8,9
13	McDill 60 – 80 ha and wood flow	2,3,4,5,7,8,9

3. RESULTS

Figure 2 shows the forest harvest scheduling for selected scenarios (1, 5 and 13). The stands selected for harvesting in years in which the minimum area constraints were applied appear highlighted.

Figure 2
Spatial harvest scheduling. The stands highlighted show the harvest scheduling during the first 5 years of the planning horizon, in which the adjacency constraints were included, for scenarios 1,5 and 13.
Figura 2
Agendamento espacial da colheita. Em destaque aparecem os talhões agendados para a colheita nos primeiros 5 anos do horizonte de planejamento, em que foram incluídas as restrições de adjacência, para os cenários 1, 5 e 13.

Considering scenario 1 in Figure 2, without inclusion of minimum area constraints, we perceived a dispersion of harvesting areas, especially in year 5 of the planning horizon. In the remaining years isolated harvesting stands occurred, but harvesting blocks were generated as well. The inclusion of the minimum area constraints leads to a reorganization in the harvest scheduling for all scenarios. Notably, it is possible to see the creation of harvesting blocks with area greater than 60 ha in scenarios 5 and 13, especially in year 5, in the southern part of the study area.

In all scenarios, the minimum area limit was respected and the largest changes occurred in 5^th year. The stands scheduled for harvesting during the 1^st and 2^nd years showed alteration in the block configuration, but most of harvesting areas remained unchanged.

Regarding the harvesting block formation, all models showed satisfactory results, being capable of creating harvesting blocks respecting the established thresholds, as exemplified by scenarios 5 and 13 in Figure 2. Therewith, there was no major differences between the two formulations regarding the compliance with a minimum area limit. Therefore, aiming at evaluating the formulation performance, additional criteria was assessed and the results are displayed in Figure 3.

Figure 3
Formulations performance. The figure shows for each scenario, including the formulations based on Carvajal et al. (2013) and Rebain and McDill (2003), in terms of the constraint number and best bound, considering the relaxed optimum.
Figura 3
Performance das formulações. A figura mostra para cada cenário testado, incluindo as formulações baseadas em Carvajal et al. (2013) e Rebain and McDill (2003), em termos do número de restrições e melhor bound, considerando o ótimo relaxado.

The results shown in Figure 3 explicit the effect of the increase in the area limit on the optimization model. With the increase in the area limit, there was an exponential increase in the number of constraints, due to the growing number of possible combination of stands to form blocks respecting the established threshold, resulting in a greater complexity for solving the optimization model. The formulation based on Rebain and McDill (2003) with a 10 ha range and the formulation based on Carvajal et al. (2013) showed similar constraint number for all minimum area limits, with the formulation based on Rebain and McDill (2003) presenting constraint number slightly smaller for minimum area limits of 30 and 40ha. For all tested area limits with a 20ha block area range the number of constraints was superior.

We observed that the bounds of scenarios based on Rebain e McDill (2003) formulation with a 10ha block area range were inferior to the bounds obtained in scenarios with a 20ha block range. With smaller range, the block formation may not include all possible stand combinations that respect the minimum area limit. Additionally, stands may not be included in any block. Therefore, the solution space is reduced, generating a limitation in the model, reflected in the scenarios' bounds. The formulation based on Carvajal et al. (2013) presented the best results, with smaller number of constraints and higher bounds. This behavior can also be recognized in the optimization results, as presented in Table 2.

Table 2
Optimization results. The table shows the values obtained for the total volume, Net Present Value, the difference related to the free scenario and the number of constraints for each scenario after the optimization process.
Tabela 2
Resultados da otimização. A tabela mostra os valores obtidos para o volume total, Valor Presente Líquido, a redução em relação ao cenário livre e o número de restrições para cada cenário após a otimização.

Scenario	Volume (m3)	NPV (BRL)	Variation (%)	Constraint number
1	1,609,073.3	62,107,951.21	-	82,742
2	1,607,172.4	61,725,398.78	-0.6%	88,212
3	1,601,475.4	61,389,986.58	-1.2%	95,172
4	1,597,700.7	61,227,921.00	-1.4%	121,382
5	1,589,151.9	60,976,371.14	-1.8%	193,357
6	1,547,129.6	58,938,743.84	-5.1%	88,187
7	1,561,452.5	59,707,283.27	-3.9%	94,952
8	1,577,441.1	60,357,385.42	-2.8%	121,767
9	1,574,370.8	60,166,499.03	-3.1%	203,867
10	1,598,461.7	61,133,267.85	-1.6%	93,452
11	1,602,819.9	61,321,530.09	-1.3%	106,967
12	1,589,450.7	60,882,169.36	-2.0%	164,417
13	1,588,114.4	60,497,667.87	-2.6%	334,527

Based on the results displayed in Table 2, we found that the greatest impact on the NPV occurred in scenario 6, with a 5.1% reduction. In general, the minimum area constraints caused reasonable impacts on the NPV. The best outcomes were achieved applying the ring inequalities (Carvejal et al., 2013), with reductions ranging from 0.6 to 1.8%. The constraints based on Rebain and McDill (2003) with a 20 ha block range showed reductions ranging from 1.3 to 2.6% and the scenarios with a 10ha block range displayed the worst outcomes in terms of NPV.

With the increase in the area limit, the number of constraints in the model increased substantially (Table 2), at an exponential rate, and consequently the solution process demanded a higher processing time. Despite the binary variables and the high number of constraints generated for the harvesting blocks formation, it was possible to achieve suitable responses in reasonable processing times for all scenarios. We obtained gaps between the relaxed optimum and the current solution below 2% after one hour of processing time.

4. DISCUSSION

All tested formulations could create harvesting blocks respecting the proposed minimum area limits. The ring inequalities formulation, proposed by Carvajal et al. (2013), presented the best results, with inferior number of constraints and suitable responses in terms of the NPV. Nevertheless, in problems with simultaneous consideration of minimum and maximum clear-cut limits, e.g. resulting from certification policies, the formulation based on Rebain and McDill (2003) may be suitable, sufficing to ensure that each stand is part of only one harvesting block. This formulation can be particularly superior when the difference between the minimum and maximum harvesting area limits is small.

The range of the block area created through the Rebain and McDill (2003) formulation had a significant impact in the number of constraints and obtained solution. In this sense, with small area range, the number of constraints is reduced. However, with small area range, the probability of stands not being included in any harvesting block increases and the number of possible solutions reduces, which may lead to inferior NPV.

The increase in the minimum area limit led to a reduction in forest NPV. Considering that stands are not uniformly distributed in a forest, in terms of the standing stock, the occurrence of neighbor stands with different areas and standing volume becomes challenging when the wood production flow is required. Thus, it is expected that the impacts on the objective function are enhanced when minimum area constraints are applied in combination with wood flow constraints. For this reason, the inclusion of minimum area constraints might lead to infeasibilities in the optimization model. In this sense, it may be necessary a relaxation of limits of wood flow constraints and management regimes, meaning harvesting stands in ages outside the usual management range, so as to increase the fluctuation limits of wood production.

The impacts of the minimum area constraints were acceptable for all scenarios tested and are compatible with values reported in the literature. For example, Öhman and Eriksson (2010) obtained a 2.6% NPV reduction for achieving a 51.8% reduction in the perimeter of harvesting areas. Our results are also compatible with related adjacency problems, such as maximum harvesting area constraints (Murray and Weintraub, 2002; Zhu and Bettinger, 2007; Gomide et al., 2010; Binoti et al., 2012; Augustynczik et al., 2015). Moreover, with the aggregation of harvesting areas and increase in operational efficiency related to road maintenance and forest harvesting, a cost reduction resulting from blocking harvesting activities is expected.

The inclusion of minimum area constraints increased substantially the complexity of the optimization model, with the addition of a considerable number of constraints with binary variables. In this sense, the application of heuristic solutions may be suitable for problems involving forests with large number of stands. A series of heuristics were proposed for solving harvest scheduling problems (Pukkala and Kurttila, 2005; Bettinger and Zhu, 2006; Gomide et al., 2009; Bachmatiuk et al., 2015), displaying suitable outcomes.

5. CONCLUSION

The formulations based on Carvajal et al. (2013) and Rebain and McDill (2003) were capable to successfully incorporate minimum area requirements, aiming to reduce negative edge effects and promote biodiversity maintenance, with acceptable impacts on the economic performance of the forest. Moreover, the formulation based on Carvajal et al. (2013) is superior when the unique manager objective is to incorporate minimum area constraints, due to the reduced number of constraints and a higher efficiency in the solution process.

6. REFERENCES

• Augustynczik ALD, Arce JE, Silva AL. Planejamento operacional da colheita considerando áreas máximas operacionais. Cerne, 2015;21(4):649-56.
• Bachmatiuk J, Garcia-Gonzalo J, Borges JG. Analysis of the performance of different implementations of a heuristic method to optimize forest harvest scheduling. Silva Fennica. 2015;49(4)e1326
• Bettinger P, Zhu J. A new heuristic method for solving spatially constrained forest planning problems based on mitigation of infeasibilities radiating outward from a forced choice. Silva Fennica. 2006;40(2):315.
• Binoti DHB, Silva Binoti MLM, Leite HG, Gleriani JM, Campos JCC. Regulation of even-aged forest with adjacency constraints. Forest Policy and Economics. 2012;20:49-57.
• Carvajal R, Constantino M, Goycoolea M, Vielma JP, Weintraub A. Imposing connectivity constraints in forest planning models. Operations Research. 2013;61(4):824-36.
• Constantino M, Martins I, Borges JG. A new mixed-integer programming model for harvest scheduling subject to maximum area restrictions. Operations Research. 2008;56(3):542-51.
• Gomide LR, Arce JE, Silva AL. Efeito das restrições espaciais de adjacência no planejamento florestal otimizado. Floresta. 2010;40(3):573-84.
• Gomide LR, Arce JE, Silva ACL. Uso do algoritmo genético no planejamento florestal considerando seus operadores de seleção. Cerne. 2009;15(4):460-7.
• Gustafson EJ. Clustering timber harvests and the effect of dynamic forest management policy on forest fragmentation. Ecosystems. 1998;1(5):484-92.
• Gustafson EJ. Expanding the scale of forest management: allocating timber harvests in time and space. Forest Ecology and Management. 1996;87(1):27-39.
• Murray AT, Goycoolea M, Weintraub A. Incorporating average and maximum area restrictions in harvest scheduling models. Canadian Journal of Forest Research. 2004;34:456-64.
• Öhman K, Eriksoon L. Aggregating Harvest Activities in Long Term Forest Planning by Minimizing Harvest Area Perimeters. Silva Fennica. 2010;44:77-89.
• Öhman K, Lämas T. Clustering of harvest activities in multi-objective long-term forest planning. Forest Ecology and Management. 2003;176:161-71.
• Pukkala T, Kurttila M. Examining the performance of six heuristic optimisation techniques in different forest planning problems. Silva Fennica. 2005;39(1):67-80.
• Rebain S, Mcdill M. A Mixed-Integer formulation of the minimum patch size problem. Forest Science. 2003;49(4):608-18.
• Tóth SF, Mcdill M, Rebain S. Finding the efficient frontier of a bi-criteria, spatially explicit, harvest scheduling problem. Forest Science. 2006;52(1):93-107.
• Wei Y, Hoganson HM. Scheduling forest core area production using mixed integer programming. Canadian Journal of Forest Research. 2007;37:1924-32.
• Zhu J, Bettinger P, Li R. Additional insight into the performance of a new heuristic for solving spatially constrained forest planning problems. Silva Fennica. 2007;41:687-98.

Publication Dates

• Publication in this collection
  2017

History

• Received
  16 Feb 2016
• Accepted
  16 Mar 2017