A stochastic programming model for an energy planning problem: formulation, solution method and application

The paper investigates national/regional power generation expansion planning for medium/long-term analysis in the presence of electricity demand uncertainty. A two-stage stochastic programming is designed to determine the optimal mix of energy supply sources with the aim to minimise the expected total cost of electricity generation considering the total carbon dioxide emissions produced by the power plants. Compared to models available in the extant literature, the proposed stochastic generation expansion model is constructed based on sets of feasible slots (schedules) of existing and potential power plants. To reduce the total emissions produced, two approaches are applied where the first one is performed by introducing emission costs to penalise the total emissions produced. The second approach transforms the stochastic model into a multi-objective problem using the ϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document}-constraint method for producing the Pareto optimal solutions. As the proposed stochastic energy problem is challenging to solve, a technique that decomposes the problem into a set of smaller problems is designed to obtain good solutions within an acceptable computational time. The practical use of the proposed model has been assessed through application to the regional power system in Indonesia. The computational experiments show that the proposed methodology runs well and the results of the model may also be used to provide directions/guidance for Indonesian government on which power plants/technologies are most feasible to be built in the future.


Introduction
Global electricity production has grown continuously year by year since 1974. According to International Energy Agency (www.iea.org), world gross electricity production has increased from 6299 to 26,730 TWh between 1974 and 2018 with an average annual growth rate of 3.3%. Electricity generation can be divided into three types, namely fossil fuel-based, nuclear, and renewable energy power generation. The fossil fuel-based plants include coal, petroleum, and natural gas which produce large quantity of cheap energy with greenhouse gas (GHG) B Chandra Ade Irawan chandra.irawan@nottingham.edu.cn emissions. The renewable energy consists of hydro, solar, wind, geothermal and biomass which generate clean energy (without GHG emissions or potentially carbon neutral in the case of biomass). However, the cost and feasibility of renewable energy are highly dependent on the potential of regions and its weather conditions (Thangavelu et al. 2015). The strategic energy planning includes the medium-to long-term expansion of the electricity generating capacity. Cost effective and low emission energy planning is very important to fulfil high demand growth for electricity in a specific country or region. An optimal energy planning will help reducing GHG emissions at the lowest cost while satisfying the electricity demand. Carbon dioxide (CO 2 ) makes up the vast majority of GHG emissions, which accounts about 72% out of GHG. CO 2 is considered to be the principle gas responsible for global warming and climate change.
This paper investigates the power generation expansion planning problem (GEP) with the objective to obtain the timing and technology choices for power plant investment over a long planning horizon in order to minimise the total cost while meeting the electricity demand and considering the environmental impact. We consider the presence of electricity demand uncertainty during the planning horizon. As the demand is predicted to increase in the coming years, optimal strategic and operational decisions are required. First, each existing power plant needs to be assessed as to whether its current generator is still being used or requires to be retrofitted with an other type of technology. Another possible decision is to shut down an existing power plant due to economic problems or environmental issues. Moreover, the optimal period for retrofitting or shutting down the power plant must also be determined. Second, to fulfil the increased electricity demand, it is important to determine the best period to construct new power plants together with their energy source type (technology) and their capacity. Finally, the optimal amount of electricity generated by each power plant need to be addressed. The deterministic multi-period mixed-integer non-linear programming (MINLP) formulation for the power generation planning problem is first built in order to minimise the total cost. Different with the models that have been proposed in the literature, the model proposed in this paper is constructed based on sets of feasible slots (schedules) of existing and potential power plants that have been generated in priori. This methodology is proposed to ensure that some important factors related to the energy planning problem are considered including the decommissioning and depreciation costs. Moreover, this also enforces that there is no idle period for an operated power plant to produce the electricity over planning horizon. As the electricity demand uncertainty is taken into account, a two-stage stochastic programming is then proposed based on the deterministic model. The strategic decisions of retrofitting and terminating the existing power plants together with constructing new plants are planned for the year ahead, then the amount of electricity generated by each power plant considered as operational decisions are determined after uncertainty is revealed.
To reduce the total emissions produced by the power plants, two approaches are put forward where the first one is performed by penalising the emissions produced by the power plants.
Here, CO 2 emissions price ($/ton CO 2 ) is introduced to calculate the total emission cost. In the second approach, the implementation of a multi-objective mathematical programming is considered to deal with multiple objectives, namely minimising total cost and minimising total emissions problems. We apply -constraint method for producing the Pareto optimal solutions. Based on preliminary experiment, the proposed energy problem is challenging to solve by an exact method using the commercial solver such as IBM CPLEX 12.8. Therefore, we also design an effective solution method that iteratively decomposes the problem into a set of smaller sub-problems to obtain good solutions within an acceptable computational time. The proposed model and its solution methods are evaluated on a realistic instance based on a power system data in the Indonesian region of Java and Bali.
The contributions of the study are as follows: (i) The introduction of a deterministic multi-period MINLP formulation for the power generation planning problem based on sets of feasible slots for existing and potential power plants. (ii) The development of a two-stage stochastic programming model for the energy planning problem under uncertainty. (iii) The construction of an effective hybrid decomposition method to solve the stochastic energy planning problem. (iv) The analysis on the effect of emission allowance price and the implementation of a multi-objective model using -constraint method to reduce total emissions.
The paper is organized into 6 main sections. In Sect. 2, an overview of the literature regarding stochastic optimisation models for energy planning is given along with the description of our contribution. Section 3 presents the technique and formulation of the proposed stochastic energy planning model. The description of the proposed hybrid decomposition method for solving the problem is provided in Sect. 4. In Sect. 5, the computational study is presented followed by conclusions and suggestions for future work in Sect. 6.

Literature review
Studies on power generation expansion planning problem (GEP) have received extensive interest due to environmental concerns and huge investment in construction of power plants. Optimisation models are developed to obtain the optimal power generation configuration in order to minimise the total cost while satisfying electricity demand and emission constraints (Hashim et al. 2005;Elkamel et al. 2009;Mirzaesmaeeli et al. 2010;Ba-Shammakh 2011;Koltsaklis et al. 2014; Koltsaklis and Georgiadis 2015b;Ahmadi et al. 2015;Elsholkami and Elkamel 2017). A number of previous works on the GEP has considered renewable energy penetration including Muis et al. (2011), Bakirtzis et al. (2012, Zhang et al. (2013) and Lara et al. (2018). Earlier studies also investigate the implementation of emission allowance prices to penalise the emissions produced by power plants (Rentizelas et al. 2012;Tolis and Rentizelas 2011;Chen et al. 2016).
The GEP have successfully been solved by a set of solution methods. Mathematical programming approach (exact method) is commonly used to solve relatively small problems. The mathematical models are usually solved by commercial solvers including CPLEX, Xpress and GAMS. The GEP can be solved by meta-heuristic techniques including, for examples, genetic algorithms (Ozcan et al. 2014;Kannan et al. 2005; Pereira and Saraiva 2011), particle swarm optimization (Kannan et al. 2005;Neshat and Amin-Naseri 2015;Moghddas-Tafreshi et al. 2011), tabu search (Kannan et al. 2005;Sadegheih and Drake 2008;Yoza et al. 2014) and simulated annealing (Kannan et al. 2005(Kannan et al. , 2007. The decomposition method is also investigated to deal with the problem such as in Sirikum et al. (2007), Flores-Quiroz et al. (2016 and Lara et al. (2018).
In this section, we review optimisation models for the GEP where uncertain factors are taken into account. Murphy et al. (1982) and Gorenstin et al. (1993) are among the first to investigate the GEP in the presence of uncertainty. They study the formulation of a stochastic program to minimise the expected total cost. Malcolm and Zenios (1994) and Mulvey et al. (1995) propose robust optimization models to obtain capacity expansion plans in the presence of uncertain electricity demand. Their models are a generalization of the model proposed by Murphy et al. (1982). Costa et al. (2017) use a robust optimization to tackle the uncertainty in investment cost for constructing a power plant. The problem is modelled using second order cone programming and semidefinite programming. Moret et al. (2019) propose a robust optimization framework considering multiple uncertain parameters in objective function and constraints. The model and methods are applied to the case study of a national energy system.
A two-stage stochastic programming model is used to deal with the stochastic GEP. Krukanont and Tezuka (2007) propose a two-stage stochastic programming model for the GEP under various uncertainties and the value of information. The worst-case scenario is used to analyse the role of various energy technologies. Feng and Ryan (2013) investigate a new scenario reduction heuristic in order to reduce the complexity of the problem. Park and Baldick (2015) propose a decomposed two-stage stochastic integer program to control CO 2 emissions where uncertain load and wind are approximated using the Gaussian copula method. Nie et al. (2017) investigate risk management of an energy system where the interval-stochastic risk management approach is proposed to address the variability of the recourse cost in stochastic programming.
The multi-stage stochastic programming can be considered as the most popular method to tackle the stochastic GEP. Min and Chung (2013) consider conditional value-at-risk in the unexpected events but expensive if occurs. Li et al. (2010) and Li and Huang (2012) analyse tradeoffs among costs, energy utilization and GHG emission control. Thangavelu et al. (2015) built a model to minimise the variations in the desired performance criteria such as energy security and costs. Betancourt-Torcat and Almansoori (2015) take into account the natural gas price uncertainty in their model. Cano et al. (2016) present an optimisation model for energy systems planning and risk management at the building level. Ioannou et al. (2019) put forward a hybrid method to handle uncertain inputs where a scenario tree configuration and Monte Carlo simulation are used. Hu et al. (2014) and Li et al. (2014) use fuzzy stochastic programming for the GEP under uncertainty. Hu et al. (2014) implement inexact fuzzy chance-constraint programming to generate solutions in order to minimise cost. Li et al. (2014) combine fuzzy, single/dual interval and multi-stage programming approaches to tackle interactions between GHG mitigation and energy-related activities.
Uncertain parameters in the stochastic GEP are usually approached by the scenario tree configuration or Monte Carlo simulation. The scenario tree is used by, among others, Thangavelu et al. (2015) Feng and Ryan (2013) and Ioannou et al. (2019). This method transforms continuous distribution into discrete scenarios where optimization at each realization of stochastic parameter is weighted with the corresponding discrete probability (Betancourt-Torcat and Almansoori 2015). Koltsaklis and Nazos (2017), Vithayasrichareon and MacGill (2012), Tekiner et al. (2010) and Min and Chung (2013) apply the Monte Carlo simulation to address the uncertain parameters. This approach generates random scenarios based on continuous distributions that can be based on historical data. In this paper, we use the scenario tree configuration to deal with the uncertain parameters. Table 1 summarises relevant previous research and highlights the position of the current study. Some important aspects that are not considered in cited papers include the possibility of plants to be retrofitted, the lead time to construct/retrofit power plants, and the depreciation and decommissioning costs. This is mainly because most available models are not constructed based on individual power plant. In other words, the data of power plants are aggregated based on the type of power generator. Therefore, it is difficult to perform a comparative analysis between our proposed stochastic model and the models available in literature. However, some deterministic models have considered retrofitted power plants including Hashim et al. (2005), Mirzaesmaeeli et al. (2010) and Bakirtzis et al. (2012). Decommissioning cost of nuclear power plant is also taken into account in the stochastic model proposed by Betancourt-Torcat and Almansoori (2015).
In this paper we propose a new methodology to build the stochastic power generation expansion planning model where a set of feasible slots (schedules) of each existing and potential power plant is initially generated in priori. A slot consists of relevant factors considered in the model including the status of a power plant (in operation or not), the decommissioning and depreciation costs, along with the construction time. The annual depreciation cost is used instead of the capital cost to accommodate the life time difference among new power plants usually determined based on the type of energy source used. A two-stage stochastic programming is designed to select the best slot/schedule for each existing and new power plant in order to achieve the total emissions target at the lowest cost considering uncertain electricity demand. Moreover, we also develop an effective hybrid decomposition method to solve the proposed stochastic model, as an exact method using a commercial solver (CPLEX) is not able to solve large problems.

The proposed optimisation model
In this section, five subsections are presented with the first discussing the proposed energy planning problem. The second subsection describes the procedures for generating feasible slots/schedule for existing or new power plants. A deterministic mathematical model based on MINLP is given in the third subsection. Then, the stochastic version of the model is presented in the fourth subsection. Finally, the description of the total emissions reduction approaches is given in the last subsection.

Problem statement
A set of existing and potential generators (power plants) is given where each power plant is classified based on the energy source namely fossil fuelled and non-fossil fuelled power plants. The fossil fuelled power plant includes coal, natural gas and petroleum whereas the non-fossil fuelled one consists of wind, solar and hydro. Each power plant (existing or potential) has known information as follows: the energy source used, the maximum installed capacity, the capacity factor, fixed and variable operating cost, fuel cost, heat rate, CO 2 emissions produced, and lifetime. As there is a possibility for an existing power plant to be retrofitted, a set of potential retrofits for each existing power plant is given. In the proposed model, the construction lead time and the capital cost required for retrofitting a power plant are needed. The most common retrofit is to transform a coal power plant into a natural gas plant to reduce the amount of CO 2 emissions produced. As a new power plant may be required to satisfy the electricity demand, similarly, each potential power plant also needs information on the construction lead time and the capital cost to build the new power plant. The construction time and the cost required to built a power plant will depend on the type of energy source used and the maximum capacity installed.
The main objective of the proposed model is to determine the optimal mix of energy supply sources in order to minimise the total cost while satisfying system constraints especially electricity demand. Other strategic decisions related to existing and new power plants are also determined as follows: The following notations are used to describe the sets and parameters of the proposed energy planning model.

Sets
P: set of power plants indexed by p ∈ P P F : set of fossil fuelled power plants (P F ⊂ P) P N F : set of non-fossil fuelled power plants (P N F ⊂ P) P E : set of existing power plants (P E ⊂ P) P R : set of existing power plants that can be retrofitted (P R ⊂ P E ) P N : set of potential power plants (P N ⊂ P) R p : set of possible alternative technologies that can be used for retrofitting power plant p ∈ P R with r as its index J : set of generation types (e.g. coal, gas, petroleum, wind, hydro, etc.) indexed by j ∈ J T : set of time periods (annually) indexed by t ∈ T Parameters k E p , k R pr , and k N p : installed capacity (MW) of the existing power plant p ∈ P E , the retrofitted power plant p ∈ P R with r ∈ R p , and the potential power plant p ∈ P N respectively. μ E p , μ R pr and μ N p : capacity factor (%) of the power plant representing the ratio of its actual output over a period of time to its potential output which is calculated based on the installed capacity. π E p , π R pr and π N p : minimum electricity generated (% of the installed capacity) by the power plant. f E pt , f R pr t and f N pt : fixed operating cost ($/MW) of the power plant in period t ∈ T . v E pt , v R pr t and v N pt : variable operating cost ($/MWh) of the power plant in period t ∈ T . u E pt , u R pr t and u N pt : fuel cost ($/MMBtu) of the power plant in period t ∈ T . h E p , h R pr and h N p : heat rate (MMBtu/MWh) of the power plant representing the efficiency of electrical power plant that concert a fuel into electricity. In other word, this parameter is the amount of energy used by a power plant to generate one megawatthour of electricity. γ : total operating hours of power plant per period. β R pr : lead time for retrofitting power plant p ∈ P R using alternative r ∈ R p (years). β N p : construction lead time for new power plant p ∈ P N (years). d t : annual electricity demand during period t ∈ T (MWh). ρ t : peak load electricity demand in period t ∈ T (MW).
The time period is based on an annual (yearly) period where the annual electricity demand (d t ) needs to be satisfied. The total installed generation capacity must be able to handle with the demand in peak load hour (ρ t ) by considering the capacity factor for each power plant (μ p ). To ensure the production continuity over the planning horizon, a minimum electricity generated by a power plant (π p ) requires to be fulfilled. If a power plant is decided to be retrofitted, the plant will not able to produce the electricity during the construction period.
In the proposed methodology, during the planning horizon (e.g. 10 or 25 years), a set of possible slots (schedule) for each existing and potential power plant is generated where a slot provides information whether a power plant is in operation or not. The slot may also give information whether a power plant is under construction as the power plant needs to be built or retrofitted. Moreover, the slot may provide related cost information including decommissioning and depreciation costs. The idea of our methodology is to explore all feasible slots for all existing and potential power plants. Then, the best slot for each power plant need to be selected in order to minimise the total cost. The main steps to generate feasible slots are presented in the next subsection. Figure 1 illustrates a simple example on how to generate feasible slots for a power plant. The slots of an existing power plant that indicate whether the power plant is in operation or not are presented Fig. 1a whereas Fig. 1b shows an example of feasible slots generated for an existing power plant to be retrofitted by a certain technology. An example of possible schedule to build a new power plant is given in Fig. 1c. Figure 1a shows the slots for an existing power plant to produce electricity within |T | periods (years). There are |T | + 1 feasible slots for each power plant where there is no intermittent in producing electricity in each slot. The status of a power plant is considered in operation at a certain period if the plant produces electricity more than the minimum that has been set (π p ). Once an existing power plant is decided to retire, it will not produce electricity again unless retrofitting this power plant is conducted. The set of feasible slots (Q E p indexed by q) of each existing power plant p ∈ P E consists of several parameters as follows:

The procedures for generating feasible slots
0 otherwise m E pq = the decommissioning cost of existing power plant p in slot q In Fig. 1b, the construction time to retrofit the power plant (β R pr ) is set to 2 years. There are (|T | − β R pr ) feasible slots for each existing power plant to be retrofitted. Retrofitting a power plant can be very expensive as new equipment needs to be installed. The capital expenditure costs for this project are depreciated over the life of the asset. Therefore, the depreciation cost for each slot needs to be determined. The set of feasible slots (Q R p indexed by q) of existing power plant p to be retrofitted consists of several parameters as follows: 0 otherwise a R prqt = the depreciation cost of power plant p for period t in slot q due to retrofitting the plant using alternative technology r In Fig. 1c, the construction time to build a new power plant (β N p ) is set to 3 years. There will be (|T | − β N p ) feasible slots to build a new power plant. The capital cost for the power plant is also depreciated over the plant life. The set of feasible slots (Q N p indexed by q) of potential power plant p ∈ P N consists of several parameters as follows: 1 if in slot q potential power plant p is under construction in period t, 0 otherwise 1 if in slot q potential power plant p is in operation in period t, 0 otherwise a N pqt = the depreciation cost for period t in slot q due to building power plant p The main steps of generating all feasible slots for each power plant is given in Algorithm 1 which consists of three stages. The first stage generates all feasible slots (schedule) for each existing power plant where the slot indicates whether the power plant is in operation or not. The second stage aims to generate all possible slots for an existing power plant p ∈ P R to be retrofitted using alternative technology r ∈ R p . The third stage deals with potential power plants that will be built if necessary. The objective of this algorithm is to explore all possible slots for potential power plants to be built. It is noted that other related information (parameters) can be easily added into each slot generated making the mathematical model closer to the real energy planning problem.

Algorithm 1
The algorithm for generating all feasible slots for each power plant 1: Stage 1: Generating operating slots for each existing power plant 2: for p = 1 to |P E | do 3: else determine m E pq for slot q 10: end for 11: end for 12: Stage 2: Generating retrofitting slots for each existing power plant 13: for p = 1 to |P R | do 14: for r = 1 to |R p | do 15: for q = 1 to (|T | − β R pr ) do 16: for t = 1 to |T | do 17: if (t < (q + β R pr )) then set c R prqt = 0 and the depreciation cost (a R prqt ) to 0. 18: else set c R prqt = 1 and determine the depreciation cost (a R prqt ) 19: if (t ≥ q and t < (q + β R pr )) then set o R prqt = 1 20: else set o R prqt = 0 21: end for 22: end for 23: end for 24: end for 25: Stage 3: Generating construction slots for each potential power plant 26: for p = 1 to |P N | do 27:

A deterministic energy planning model
In this subsection, we present a deterministic mathematical model based on an MINLP for the proposed energy planning problem. Six decision variables are constructed which consists of three binary variables and three continuous variables. The binary variables (X ) is designed to select the optimal slot that will be used by each power plant whereas the continuous variables (E) determine the amount of electricity produced (MWh) by each power plant per year. The description of the decision variables is given as follow: Decision Variables pr t = annual electricity generated by retrofitted power plant p ∈ P R using alternative technology r ∈ R p in period t ∈ T E N pt = annual electricity generated by potential power plant p ∈ P N in period t ∈ T The MINLP model is expressed as follows: where subject to: The objective function (1) consists of two components which are related to the strategic decision variables X (Z x ) and the operational decision variables E (Z e ). The first objective function (Z x ) includes the fixed cost, the decommissioning cost if the existing power plants are decided to retire, the retrofitting cost of existing power plants and the constructing cost of building new power plants. The second component cost (Z e ) includes the variable O&M and fuel costs together with the heat rate, which is calculated based on the power plant capacity and amount of electricity produced. The forth term of the objective function Z x in Equation (2) makes the model non-linear as there is a product of two decision variables (X E pq and X R prq ). However, in the implementation, the model can be linearized by adding new variables and constraints.
Constraints (4) ensure that one slot (schedule) of an existing power plant must be selected. Constraints (5) guarantee that if an existing power plant is retrofitted, only one type of technology can be used over the planning horizon. Constraints (6) aim to decide whether a new power plant should be built or not. These constraints also determine the best period (year) when to build it. Constraints (7-8) impose that an existing power plant cannot be in operation when the power plant is under construction for retrofitting or after it has been retrofitted. In other words, the time for an existing power plant to produce electricity cannot overlap with time for retrofitting the power plant. Constraints (9) enforce that the total installed generator capacity is able to satisfy the peak load electricity demand for each period by taking into account the capacity factor for each power plant. Constraints (10-12) ensure that the amount of electricity produced by a power plant is less than its capacity factor and more than the minimum amount of electricity that should be generated. Here, the annual production is calculated based on the power plant capacity (MW) and the total operating hours per year (γ = 8760 hours). Constraints (15) ensure that the annual electricity generated by the power plants must satisfy the annual electricity demand. Constraints (14)(15)(16)(17)(18)(19) indicate that X variables are binary whereas E variables are treated as continuous.

A stochastic programming model
In this section, a stochastic version of the proposed energy planning problem is introduced. In particular, it is assumed that annual electricity and peak demands over the planning horizon are not known in advance and are considered as uncertain parameters. However, the demands can be captured by a probability distribution. Here, the decisions to retrofit and terminate the existing power plants together with to build new ones have a long-lasting effect. Therefore, these are strategic decisions which are typically determined before having precise information about the demands in future. On the other hand, the decisions about the amount of electricity generated by each power plant are often considered as operational decisions that can be made when perfect information is available, after uncertainty is revealed.
To include the random nature of the electricity demands in the problem, we consider a twostage stochastic programming with multi-period model (see Birge and Louveaux 2011) which reflects the way in which the uncertainty in electricity demands is expressed and revealed during the planning horizon. In the two-stage stochastic program, we assume that strategic decision variables (X ) for all periods must be determined before the actual realization of the electricity demands as the first stage. The second-stage concerns the amount of electricity generated which is related to operational decision variables E. This stage is considered as recourse decisions because they are determined in such a way that the best response on the occurring scenario is given to the setting defined by the first-stage decisions.
In this approach, for every t ∈ T , the annual demand (d t ) and peak load (ρ t ) per period are assumed to be random variables. Each realization of these random variables is a scenario and it is assumed that it is possible to compute or estimate accurately the probability associated with each scenario. Here, all possible scenarios of future conditions are explored. Stochastic optimization aims to minimise the expected cost of the system. Let S = {1, . . . , n s } be the set of scenarios indexed by s ∈ S. The stochastic parameters can be denoted as d ts and ρ ts for the annual and peak demands respectively in period t ∈ T under scenario s ∈ S. Accordingly, decision variables E E pts , E R pr ts and E N pts are the proportion of the electricity generated in period t ∈ T by existing, retrofitted and new power plants respectively under scenario s ∈ S. Let w s be the probability associated with scenario s ∈ S. The main objective is to obtain the optimal mix of energy supply sources and to determine the electricity generated by each power plant in a way that performs well in every possible situation (scenario). The stochastic energy planning model with scenario-based approach can be formulated as follows: subject to: (2), (4)- (8), (14)

The total emissions reduction approaches
To reduce the total emissions produced by the power plants, two approaches are applied where the first is performed by penalising the CO 2 emissions produced by the power plants.
In the second approach, a multi-objective mathematical programming based on -constraint method is implemented to deal with multiple objectives, namely minimising total cost and minimising total emissions.

The implementation of the emission allowance price
One way to reduce the total CO 2 emissions is to penalise for each tonne of CO 2 produced by power plants. Here, a CO 2 emissions allowance price ($/ton CO 2 ) is introduced to determine the total emissions cost for producing CO 2 emissions. To accommodate this approach, additional parameters are introduced as follows: ε E p , ε R pr and ε N p : CO 2 emissions (tonne of CO 2 /MWh) caused by power plant p to generate electricity. e t : CO 2 emissions price ($/ton CO 2 ) in period t ∈ T .
The total emissions cost is expressed by the following formulation: The objective function (20) of the stochastic model is then enhanced by considering the total emissions cost (29). The stochastic model that considers the emission allowance price is then reformulated as follow: s.t.
(2), (4)- (8), (14)- (16) and (21)-(28) By adding the total emissions cost in the objective function, it is assumed that the solution of the model will try to reduce the total CO 2 emissions produced by the power plants. This can be done by switching the use of fossil-fuelled power plants to the renewable energy-based generators.

A multi-objective mathematical programming based on the -constraint method
A multi-objective mathematical programming is used to deal with the problems that have more than one objective function. This technique has been widely implemented for addressing real-life applications (Song et al. 2019;Yakavenka et al. 2019;Kaveh et al. 2019). The generic multiple objective programming model (Steuer 1986) is formulated as follows: where x is the vector of decision variables, ( f 1 (x), f 2 (x), . . . , f n (x)) are the n objective functions. F is the feasible region and is defined by a number of constraints given as functions of x. Several methods can be used to address multi-objective problems including Pareto efficient set generation, compromise programming and goal programming. The nature of the achievement function is affected by the multi objective method technique used (Jones 2011). In the efficient set generation method, it is considered as an optimisation of a vector to produce Pareto efficient solutions (Steuer 1986;Maddah et al. 2019;Torabi Yeganeh and Zegordi 2020). In compromise programming, a distance function is applied to measure the closeness between the ideal and efficient solutions. In the goal programming method, unwanted deviations from goal levels are minimised.
In this study, -constraint method (Haimes et al. 1971) is used to solve the multi-objective energy planning problem. -constraint method is one of Pareto efficient set generation methods that works by pre-defining a virtual grid in the objective space and solving different single-objective problems constrained to each grid cell (Laumanns et al. 2006). The problem is transformed to a single-objective problem with limiting the others by some allowable values i , i ∈ {1, . . . , n}. In the -constraint method, one of the objective functions is optimized by incorporating the other objective functions as constraints of the model which is expressed as follows: By varying the values i , the Pareto efficient solutions of the problem are obtained. In this approach, the proposed energy planning problem is treated as a bi-objective problem where two objective functions are considered, namely minimising the total cost ( f 1 ) and minimising the total emissions produced ( f 2 ). The former function is based on Eq. (20) where f 1 = Z c . The latter is determined based on the total electricity produced by power plants together with their energy sources. Similar to the previous subsection, parameters ε E p , ε R pr and ε N p are used indicating the CO 2 emissions (tonne of CO 2 /MWh) caused by power plant p to generate electricity. The bi-objective problem is formulated as follows: s.t.
(2), (4)- (8), (14)- (16) and (21)-(28) Algorithm 2 presents the main procedure of the proposed -constraints method where the first step is to transform the bi-objective problem into a single-objective problem. Here, f 1 (minimising the total cost) is treated as the main objective function whereas f 2 (total emissions) is considered as a constraint of the problem. The single-objective problem can be defined as follows: (2), (4)- (8), (14)- (16) and (21)-(28) Algorithm 2 The proposed -constraints method 1: Transform the bi-objective problem into a single-objective problem by considering the total cost ( f 1 ) as the main objective function. 2: Determine the upper bound of 2 (¯ 2 ) by solving the single-objective problem without Constraint 36. 3: Define the lower bound of 2 ( 2 ) where 2 = λ ·¯ 2 4: Define parameter ϕ for varying the value of 2 . 5: Set 2 =¯ 2 6: while 2 > 2 do 7: Solve the bi-objective mathematical model by considering the 2 -constraint 8: Store the solutions and the objective function values ( f 1 and f 2 ) 9: 2 = 2 − ϕ ·¯ 10: end while 11: Generate a set of Pareto optimal solutions based on the solutions generated in the previous step. Decision makers will select the preferred solution from the Pareto optimal solutions.
The upper bound of 2 (¯ 2 ) is obtained by solving the single-objective problem that minimises the total cost only. Here the stochastic model represented by Eqs. (20), (2), (4)-(8), (14)-(16) and (21)-(28) are optimally solved. In other word, Constraint 36 is not incorporated in the model. The solution obtained is used to calculated the total emissions produced ( f 2 ) which is considered as the upper bound of 2 (¯ 2 ). The lower bound of 2 ( 2 ) is then determined based on¯ 2 . This can be performed by introducing parameter λ representing the realistic maximum total emission reduction (%) that can be achieved ( 2 = λ ·¯ 2 ). In the next step, parameter ϕ is defined to vary the value of 2 where a set of bi-objective problems are constructed. Each bi-objective mathematical model is solved and the solutions are then used to generate Pareto efficient solutions. Decision makers will then select the preferred solution from the Pareto optimal solutions.

A hybrid decomposition method
An exact method can be used to solve the proposed MINLP model for obtaining the optimal solution for the medium-and long-term energy planning. Here, we use a commercial optimizer called IBM ILOG CPLEX ver. 12.8. However, the exact method was not able to solve the problem optimally especially when dealing with a relatively large problem as CPLEX terminates due to being out of memory before attaining the optimal solution. Therefore, we propose a hybrid decomposition method to overcome this limitation and to obtain good solutions in an acceptable time. The main steps of the hybrid decomposition method are depicted in Algorithm 3 which consists of two stages. In the first stage, the main problem is decomposed into several sub-problems based on the period where |T | sub-problems are constructed. This stage aims to obtain a good initial solution where the slot for each power plant (arraysX E pq ,X R prq andX N pq for the existing, retrofitted and new power plants respectively) is generated. This relatively good initial solution is then fed to the second stage where the main problem is iteratively reduced into several restricted problems based on the type of power plants.

Algorithm 3
The proposed decomposition method 1: Stage 1: period-based decomposition method 2: Generate an initial solution using the procedure given in Algorithm 4. 3: LetX E pq andX N pq be the obtained solution for the existing and new power plants respectively. 4: Determine the objective function value (Z ) by solving the original stochastic model using obtained solution (X E pq andX N pq ). 5: Stage 2: restricted problem-based iterative method 6: UseX E pq andX N pq generated from the first stage as an initial solution. 7: Implement the proposed Algorithm 5 withZ ,X E pq ,X R prq andX N pq as inputs and outputs. The amount of electricity produced by each power plant (E E pts , E R prts and E N pts ) is also determined.

The period-based decomposition method
The description of the proposed period-based decomposition for solving the stochastic energy planning problem is presented in this subsection. The main objective of this method is to obtain a relatively good solution to be used as the initial solution for the next stage. Algorithm 4 describes the procedure of the period-based decomposition approach which consists of two phases.
In the first phase, the main problem is decomposed into |T | sub-problems where each sub-problem is solved optimally using the exact method. The sub-problem that represents the problem for period |T | is first constructed and solved. The solution of this sub-problem is then used to construct the sub-problem for period |T | − 1. This process is repeated until the sub-problem for the first period is designed and solved. In the second phase, the initial solution for the main problem is then constructed based on the solutions obtained from solving the sub-problems.
First, array G pt is constructed to store the solutions obtained from solving the subproblems. An MINLP model is designed for each sub-problem that represents the energy planning problem for each period (year) starting from the end of period. The MINLP model for each sub-problem is the reduction of the proposed stochastic model where we deal with one-period problem instead of multi-period problem. Here, decision variables X E pt , X N pt and X R pr t are replaced byX E pq ,X R prq andX N pq respectively. The problems for period t where t < |T | are even have little computational effort as some of decision variables are fixed due to the solution that has been obtained for the following period. For instance, if an existing power plant is in operation in period (t + 1) then this power plant needs to operate in period t.
In Lines 14-25 of Algorithm 4, the MINLP model of the sub-problem is optimally solved by an exact method where the obtained solution is stored in array G pt . This array has 3 types of value where 0, 1 and 2 indicate that the power plant is not in operation, under construction and in operation respectively. The value in this array for each plant in period t is determined based on the obtained solution from solving the sub-problem in period t (X E p andX N p ) and the value in the following period (t + 1). Finally, once all sub-problems have been solved, Algorithm 4 The period-based decomposition approach 1: Construct arrays G pt p ∈ P, t ∈ T 2: Phase 1: 3: for t = |T | to 1 do 4: Consider the stochastic model for period t only (sub-problem) 5: if t < |T | then 6: In the sub-problem, fix the value of decision variables as follows: 7: for each existing power plant p ∈ P E do the following step: 8: if G p(t+1) = 2 then the existing plant p is in operation (i.e.X E p = 1). 9: end for 10: for each potential power plant p ∈ P N do the following steps: 11: if G p(|T |) = 0 or G p(t+1) ≤ 1 or t < β N p then the potential plant p is not in operation (i.e.X N p = 0) 12: end for 13: end if 14: Solve optimally the sub-problem using an exact method and the solution is used to update arrays G pt with the following procedure: 15: for each existing power plant p ∈ P E do the following step: 16: ifX E p = 1 then G pt = 2 17: end for 18: for each potential power plant p ∈ P N do the following step: 19: ifX N p = 1 then 20: the initial solution for the main problem is constructed based on arrays G pt andG pr t . The procedure to generate this solution is presented in Lines 28-40 of Algorithm 4. This solution is then fed to the second stage of Algorithm 3 as an initial solution. In the next sub-section, the description of the problem restriction-based iterative approach is presented.

The restricted problem-based iterative method
In the restricted problem-based iterative approach, a set of sub-problems are constructed based on the type of power plants. The main steps of the proposed approach is presented in Algorithm 5 which aims to seek the best configuration of existing and potential power plants for each generation type. To construct a sub-problem, the main stochastic model presented in Sect. 3.4 is reduced by considering the power plants with a certain generation type only. This can be done by fixing the binary variables (X E pq ,X R prq andX N pq ) for the other power plants with different generation types based on the incumbent solution. The restricted problem is then solved using the exact method within τ 1 seconds. If an improvement is found then we update the incumbent solution (X E pq ,X R prq ,X N pq , E E pts , E R pr ts and E N pts ) along with its objective function (Z ). This process is repeated iter max times and the best solution (X E pq , X R prq andX N pq ) is taken. The main problem is then optimally solved using the exact method where all the binary variables are fixed based on the values of (Z ,X E pq ,X R prq andX N pq ) which are obtained in previous steps. Here, the problem is transformed into a linear programming (LP) model which is relatively easy to solve. The amount of electricity produced by each power plant and each scenario (E E pts , E R pr ts and E N pts ) together with its objective function value (Z ) are then obtained.

Computational study
This section presents the experimental results that have been carried out to evaluate the performance of the proposed model and solution method. The experiments are conducted based on a real case study of a regional power system in Indonesia (Java and Bali). To assess the

Algorithm 5 The restricted problem-based iterative method
Require:Z together with arraysX E pq ,X R prq andX N pq . 1: forî = 1 to iter max do 2: for each generation type j ∈ J do 3: Reduce the main problem by considering existing power plants with generation type j only. In other words, we fixed the binary variables X E pq , X R prq and X N pq for the other power plants based on the incumbent solution (X E pq ,X R prq andX N pq ).

4:
Solve the restricted problem using an exact method within τ 1 s. LetZ be its objective function value and store the obtained binary variables into arraysX E pq ,X R prq andX N pq .

5:
ifZ <Z then 6: UpdateZ =Z ,X E pq ←X E pq ,X R prq ←X R prq andX N pq ←X N pq . 7: end if 8: end for 9: Do the same as Lines 2-8 with a minor revision where in Line 4, the best configuration of potential power plants with generation type j is determined instead of existing power plants. 10: end for 11: Solve the main model optimally using an exact method where all the binary variables are fixed based on the values of (Z ,X E pq ,X R prq andX N pq ). In other words, the problem is transformed into a linear programming (LP) model. The amount of electricity produced by each power plant and each scenario (E E pts , E R prts and E N pts ) together with its objective function value (Z ) are then obtained. 12: ReturnZ ,X E pq ,X R prq andX N pq along with E E pts , E R prts and E N pts .
performance of the proposed hybrid decomposition method (HDM), the experiments are first performed for the medium term planning horizon (10 years). Then, the analysis on experimental results for long term planning horizon (25 years) is presented. The implementation of the proposed hybrid decomposition method was written in C++ .Net 2017 and used the IBM ILOG CPLEX version 12.8 Concert Library. The tests were run on a PC with an Intel Xeon W-2133 CPU @3.60 GHz processor, 64.00 GB of RAM. Indonesia is the largest region in South-East Asia with excess availability of natural gas, geothermal, and hydro resources. In this country, the electricity demand is concentrated in Java and Bali islands wherein more than half the population live. The detailed information on existing and potential power plants along with the electricity demand is taken from the report provided by PLN (2019). In summary, in 2019 there are 135 existing power plants in this region with different capacity and generation type. The generation sources used for the existing power plants include coal, petroleum, natural gas, hydro and geothermal. There are 23 coal power plants that account for almost 47% of the total installed capacity (approximately 44 GW). Indonesia has committed to reducing CO 2 emissions under the climate agreement. Current commitment is 29% of greenhouse gas reduction in 2030 relative to 2010 under a business as usual baseline (Wijaya et al. 2017). Table 2 presents the cost and technical details of power generation types considered in this case study. The depreciation and decommissioning costs are determined based on capital cost and life time. The straight-line depreciation method is used to calculate those costs where salvage value is used to determine the decommissioning cost. The construction lead time data is mainly taken from Chen et al. (2016) with minor revision. The forecasted cost of fossil fuels for the electric power industry is taken from U.S. Energy Information Administration (EIA-www.eia.gov). In this study, we consider one type of method for retrofitting an existing power plant, which is transforming a coal/petroleum power plant to natural power plant. The detailed data of retrofitted power plants is taken from EIA.

Experimental results on medium-term planning horizon
The report provided by PLN (2019) reveals the government plan to build new power plants for the next 10 year (2019-2028) to satisfy the future electricity demand. In this paper, these power plants are considered as potential power plants. Based on their plan, there are 116 potential power plants that can be built in this region with different capacity and generation type including the ones with small capacity (< 10 MW). The generation types used in the potential power plants include hydro, wind, biomass, natural gas, geothermal, solar and pump storage.
PLN (2019) also provides the estimated annual electricity demand for Java and Bali regions from 2019 to 2028 where the demand is expected to increase from 180,806 GWh to 301,085 GWh. Moreover, the report also presents the estimated peak load demand which is projected to rise from 28 GW in 2019 to 45 GW in 2028. As it is expected that the annual demand and the peak load have the same trend, it is assumed that the growth of these demands is the same. According to Thangavelu et al. (2015), the annual demand growth for Indonesia can be classified into three categories, namely low, medium and high with the occurrence probability of 0.3, 0.55 and 0.15 respectively. This growth rate is dependent on population growth, urbanization and international actions. It is assumed that the low, medium and high demand growth are 2%, 4% and 6% respectively which is based on the forecasted demand given by PLN (2019). In this case study, the optimal energy mix for the next 10 years is determined and we assume that the growth rate may change every two years. In other words, Table 2 Cost and technical details of power generation types . Source: IEA (www.eia.gov) and Thangavelu et al. (2015) Generation To evaluate the performance of the proposed HDM, we compare the solutions obtained by the HDM with the ones found using the exact method (EM). In the EM, IBM ILOG CPLEX version 12.8 is used to solve the stochastic energy model presented in Sect. 3.4. In solving the problem (for a single-objective or a multi-objective), the computational time (CPU) using the exact method (CPLEX) is limited to 3 h (10,800 s). By limiting the computational time in CPLEX, the lower bound (LB) and upper bound (UB) can be obtained. The performance of the proposed HDM will be measured by Gap (%) between the Z value attained by the proposed method and the lower bound (LB) obtained from the exact method. Gap (%) is calculated as follows: where Z m refers to the feasible solution cost obtained by either the exact method (UB) or the proposed decomposition method.
In the proposed decomposition method, the value of τ 1 is set to 50 s which is determined based on our preliminary study. There are two computational experiments that are conducted to deal with total costs and total emissions. First, computational experiences are carried out to analyse the effect of the emission allowance price to the total emission reduction. Second, a multi-objective energy planning problem where two conflicting objectives namely minimising total cost and total emission (CO 2 ) is addressed using the -constraint method.

The effect of the emission allowance price
In this study, four emission allowance price scenarios are considered based on the work of Rentizelas et al. (2012). The emission allowance price (e t ) is classified into four categories namely zero, low, medium and high prices which are given in Table 3. Here, the proposed energy problem presented in Sect. 3.5.1 is solved to minimise the total cost considering the emission allowance price.
In Table 4, results of the energy planning problem with several emission price scenarios are given which are obtained using the exact method (EM) and the proposed hybrid decomposition method (HDM). Here, Z best is the best objective function found by either EM and HDM. During the experiments, we noticed that the EM addressed the problems well as the Gaps (%) between UB and LB obtained are very small. We noticed that the EM produced the optimal solutions when low, medium and high CO 2 price are used. Based on the table, the proposed HDM also runs well as it produced good solutions (0.3937% Gap on average) within an acceptable computational time (373 s on average). The solution generated is then used to calculate the amount of CO 2 produced by power plants. It can be noted that the total emissions (ton CO 2 ) produced using the solutions generated by the EM are relatively similar with the one obtained by the proposed HDM.
According to Table 4, penalising the emissions produced by power plants generates the energy mix that reduces the total CO 2 emissions produced. The use of low, medium and high CO 2 emission allowance prices reduces the total emissions by 22.59%, 23.66% and 26.12% respectively compared to the one produced using zero CO 2 emission allowance price. In this case study, the implementation of CO 2 emission allowance price significantly changes the structure of energy mix. When the high price is applied, the capacity generation mix is arranged to reduce the emission cost. In other words, new non fossil-fuelled power plants are preferred to build rather than the fossil ones to satisfy the demand. Table 5 presents the cost using several emission price scenario breakdown by existing plant costs (EPC), retrofitted plant costs (RPC), new plant costs (NPC), depreciation costs (DPC), decommissioning costs (DCC) and emission costs (EMC). The cost components of EPC, RPC and NPC consist of fixed, variable and fuel costs. The cost of existing plants contributes the highest portion to the total cost. As expected, the emission costs increase with the CO 2 emission allowance price. Although some existing plants are terminated, no decommissioning cost occurs as the plants have reached their life time. CPU is measured in seconds  Table 5 also reveals the cost breakdown when the emission cost is not included to form the total cost. According to the table, the existing plants costs decrease with the emission price which is in contrast with the new plants cost. New non-fossil fuelled plants require to be built in order to reduce the emission cost. However, this leads to higher depreciation cost.

The application of a multi-objective problem using the -constraint method
This subsection presents the computational experiments when the proposed stochastic energy planning problem is treated as a multi-objective problem which is addressed by theconstraint method. According to the previous experiments, solving the minimising total cost problem will have the energy mix configuration that produces CO 2 emission of 1,589,803,226 tonnes for the next 10 years. Therefore, the upper bound of total emissions (¯ 2 ) is set to 1,589,803,226. It is assumed that a realistic maximum total emissions that can be reduced over the planning horizon is set to 35% of the upper bound of total emissions. Parameter λ is then set to 0.65 meaning that lower bound of total emissions ( 2 ) is set to 1,033,372,097. As we aim to generate Pareto efficient solutions, we set parameter ϕ to 2.5% to vary the value of 2 . Here, we vary the emissions reduction target from 20 to 35% out of the total CO 2 emissions produced at the lowest cost (¯ 2 ).
Tables 6 shows the computational results for a set of the energy planning problems for medium-term planning horizon where each problem has different CO 2 emission reduction. The problems are solved by the exact method (EM) and the proposed hybrid decomposition method (HDM). The first two columns of Table 6 indicate the emission target ( 2 ) that needs to be achieved. The same as previous experiments, the CPU time in solving each problem using the EM (CPLEX) is limited to 3 h (10,800 s). Based on the table, CPLEX was able to obtain upper and lower bounds (U B and L B) for three problems only. CPLEX failed to attain the U B for the problems with 25% and 32.5% total emissions reduction. However, relatively good lower bounds (L B) are obtained for those problems.  Table 6 also reveals that the proposed HDM runs much faster than the EM while producing good Gap(%). Based on the average result, the proposed HDM yields a gap of 0.6468% whereas the EM produces 0.0696% (calculated based on the obtained gaps). In general, it can be concluded that the proposed HDM performs well in term of the quality of solutions obtained and the computational time. In these experiments, the EM produced a very good lower bound (LB) for all problems. Table 7 presents the cost breakdown on each CO 2 emission reduction target. It can be noted that a higher total emission reduction target yields a higher total cost. As expected, more new non-fossil fuelled plants require to be built if higher emission target is implemented. Moreover, there is a significant increment in the retrofitted cost for 35% emission target as many coal plants need to be transformed into gas plants to reduce the emissions produced. Figures 3 and 4 show the generation capacity mix and the energy mix for the mediumterm planning horizon respectively. The figures compare the results between 0% (without emission constraints) and 35% emission reduction. The figure shows that the electricity demand is mainly satisfied from fossil-fuelled power plants (coal and natural gas). However, the use of coal power plants is gradually reduced over the planning horizon and replaced by natural gas power plants as the main energy source as natural gas produces less CO 2 emissions than coal.
Without total emission constraints, around 5.8 GW new non-fossil fuelled plants and 11.9 GW gas plants are constructed during the planning horizon. Moreover, almost all petroleum plants with the capacity of 2.1 GW in total require to be transformed (retrofitted) into gas power plants at the beginning of period. When 35% emission reduction target is implemented, the power plant configuration is changed where almost 7.5 GW new non-fossil fuelled plants together with 9.6 GW new gas plants are built. Over the planning horizon, all petroleum plants with capacity of 2.3 GW together with 14.7 GW coal plants need to be converted into natural gas. In summary, more non-fossil fuelled power plants are added when 35% emission reduction is applied. Figure 5 presents the comparison of total emissions produced between 0 and 35% emission reduction. The total emissions produced decrease at the beginning of period, however it increases over the year with electricity demand.
The solutions generated by the HDM are used to generate Pareto efficient solutions presented in Fig. 6. Decision makers will then choose the solution from the Pareto optimal Table 7 The cost breakdown on each CO 2 emission reduction target for medium-term plan % Emissions reduction  solutions based on their preference. It is worthwhile noting that the Pareto solutions cannot be constructed using the EM as not all problems can be solved by this method within 10,800 s. It indicates that the stochastic bi-objective problem is hard to solve by the EM. The proposed HDM is then designed to overcome this issue.

Experimental results on long-term planning horizon
In this subsection, experimental results on long term analysis are conducted where the problem considers 25 years planning horizon from 2019 to 2043. In addition to 116 potential power plants that can built during 2019-2028, PLN (2019) also provides other 83 potential non-fossil fuelled generations that can be constructed for the long-term planning. To satisfy the long-term electricity demand, we generate 94 other potential power plants so the problem consists of 293 potential power plants in total that can built during the next 25 years period. The same as previous experiments, the growth for annual demand and peak load can be classified into three categories (low, medium and high). In this long-term analysis, the growth demand rate for each scenario may change every five years. Therefore, there are 3 5 = 243 possible demand scenarios presented in Fig. 7. In summary, the problem comprises 283,475 binary and 2,849,175 continuous variables.
In these experiments the proposed stochastic energy planning problem is first solved to minimise the total cost without considering the maximum total emissions produced by using the zero CO 2 emission allowance price. Based on the solution generated, the total emissions is calculated. Then, the model is solved to meet 35% emissions reduction target at the lowest cost. The long-term problem was solved by the proposed Hybrid Decomposition Method (HDM) only as the exact method (EM) was not able to generate feasible solution for the problem. CPLEX used in the EM unfortunately terminates when solving the problem due to being out of memory. Tables 8 presents the experimental results for long-term energy planning problems obtained by the proposed HDM. The table reveals the breakdown of the costs which includes the total cost ($) together with the percentage of each cost type to the total cost.
As expected, the total cost increases with the emission reduction target which confirm the previous experimental results. According to Table 8, new power plants contributes the highest portion (41.47%) to the total cost when 35% emission reduction is applied. Implementing 35% emission reduction target increases new non-fossil fuelled plants to be installed and the number of existing coal power plants to be transformed into natural gas. The solutions generated indicate that several power plants require to be shut down, however these plants have reached their life time so no decommissioning costs exists. Figures 8 and 9 presents the comparison of capacity and energy mix for 2019 to 2043 between 0 (without emission constraints) and 35% emission reduction. The figures shows that the electricity demand is mainly satisfied from fossil-fuelled power plants (coal and natural gas). However, the use of coal power plants is gradually reduced over the planning horizon and replaced by natural gas power plants as the main energy source as natural gas produces less CO 2 emissions than coal. For the 35% emission reduction target, more non-fossil plants are used to generate electricity to reduce the total emissions at the expense of a higher total cost.
Minimising total cost without total emission constraints leads to the following results: 54.4 GW new non-fossil fuelled plants and 65.9 GW natural gas plants are built; 2,254 MW petroleum plants are transformed into natural gas at the beginning of period. When 35% emission reduction target is implemented, the power plant configuration is changed where 60.4 new non-fossil fuelled power plants together with 59.9 GW new fossil plants are constructed over the horizon planning. In addition to retrofitting the petroleum plants, almost 2 GW coal plants also converted into natural gas. Figure 10 presents the comparison of total emissions produced between 0 and 35% emission reduction.

Conclusions
In this paper, a stochastic programming model is proposed to obtain the optimal mix of energy supply sources for medium to long-term planning horizon in order to minimise total electricity generation cost. The deterministic MINLP model is first constructed based on feasible slots (schedule) for each power plant. Based on the deterministic model, a twostage stochastic model is then built to accommodate uncertain annual demand and peak load. Two approaches are put forward to reduce the total emissions produced where the first one is conducted by introducing carbon emission prices to penalise the total carbon emissions. The second approach is performed by implementing a multi-objective problem using the -constraint method for producing the Pareto optimal solutions. A hybrid decomposition method is then designed to generate good solutions within an acceptable computational time.
A dataset based on a real case study of the regional power system in Indonesia (Java and Bali) is used to evaluate the performance of the solution method. In this study, the energy mix for the next 10 and 25 years for this region is obtained. Based on the results of computational experiments, the introduction of CO 2 emission allowance price for penalising the total emissions reduces the total CO 2 produced. The implementation of high emission price decreases the total emissions by 26.12% for the medium-term planning horizon. The implementation of a multi-objective problem using the -constraint method is also found to be an effective approach to reduce the total CO 2 emissions. The Pareto optimal solutions can efficiently be generated by the proposed decomposition method. Here, decision makers will then choose the preferred solution from Pareto solutions. The application of 35% emission reduction target significantly changes the structure of energy mix for medium-and long-term planning horizon.
The following research directions may be worthy of investigation in future. In the proposed model, considering every year in the scenario tree exponentially increases the number of scenarios. This also increases the complexity of the model. Scenario reduction technique can be designed to address this issue. This study could be extended by integrating the proposed long-term generation expansion planning (GEP) methodology with the well-known shortterm unit commitment problem (daily energy planning). This can be performed to increase the decision accuracy power networks' stability (Koltsaklis and Georgiadis 2015a).