Exploring the Dynamic Landscape of Economic Dispatch: A Comprehensive Review of Stochastic Programming Techniques
D. Punsara Colambage1, W.D. Anura S. Wijayapala2, Tilak Siyambalapitiya3
1,2 Department of Electrical Engineering, University of Moratuwa, Sri Lanka.
3 Ceylon Electricity Board, Colombo 02, Sri Lanka.
ABSTRACT
Ensuring the provision of high-quality and continuous power at an affordable rate stands as the primary objective of power systems. To achieve this goal, utilities devise economic dispatch schedules for each quarter-hour time slot, leveraging stochastic programming methodologies. This becomes imperative due to the inherent uncertainty associated with renewable energy sources, heavily reliant on weather and climate conditions. Dynamic Programming, Stochastic Dynamic Programming, and Stochastic Dual Dynamic Programming techniques are deployed to formulate real-time economic dispatch schedules. Notably, while cost remains a pivotal factor, the emphasis shifts towards economic value, particularly water value, during the dispatch process. Various methodologies have been introduced to ascertain the optimal economic dispatch schedule, considering the dynamic landscape of economic dispatch. This landscape plays a crucial role in the design, analysis, planning, and modelling of smart energy systems, facilitating informed decisions for the green transformation of energy supply and demand systems towards future smart renewable energy solutions. This research paper offers a comprehensive literature review of stochastic programming techniques applied in economic dispatch, shedding light on their implications for smart energy systems.
Keywords: Economic Dispatch, Stochastic, Water Value, Dynamic Programming
INTRODUCTION
In the intricate dance of managing a power system, utility operators face the daunting task of orchestrating the ideal blend of power plants to meet fluctuating electricity demands. At the heart of their challenge lies the short-term optimisation puzzle: how to meticulously schedule generation to minimise total costs while adhering to a myriad of constraints. With an eye toward efficiency and sustainability, power system owners constantly seek to optimise the balance between constructing new, state-of-the-art power plants and retiring older, less efficient ones.
Responsibility for meeting electricity demand rests squarely on the shoulders of the utility company, underscoring the critical importance of ensuring that global generation matches forecasted demands on a quarter-hourly basis. Unit Commitment and Economic Dispatch emerge as pivotal factors influencing generation costs. Unit Commitment involves the strategic decision-making process of determining when and which generating units at each power station to activate or deactivate, while Economic Dispatch entails fine-tuning individual power outputs to align with scheduled generating unit capacities at each time interval [1][2].
At its core, the unit commitment process is driven by the overarching objective of minimising operating costs while satisfying all operational and systemic constraints. This intricate ballet of generation scheduling involves meticulous planning of startup, shutdown, and power levels for each unit over a defined scheduling horizon. The aim is to strike a delicate balance, minimising operational expenses, startup, and shutdown costs while ensuring seamless continuity of power supply [3].
Delving deeply into the concept, Unit Commitment intricately navigates the optimal scheduling and production levels of each generator unit within a power system for a given time frame, meticulously balancing operational constraints and equipment requirements. On one front, it meticulously delineates the commitment status of generators—whether they are poised to be activated or remain dormant—during the corresponding period. Subsequently, Economic Dispatch steps in to meticulously ascertain production levels within smaller time intervals [4][5].
Traditionally, Unit Commitment primarily addresses security-based constraints such as line outages and transmission line capacities, often overlooking network constraints. Conversely, from the perspective of market operators, Unit Commitment seamlessly integrates to minimise costs and maximise the benefits of generation targets [6].
Economic Dispatch unfolds into deterministic and stochastic realms, navigating uncertainties stemming from forecast deviations and equipment reliability issues. Fluctuations in forecasts can trigger load variations, compounded by the intermittent and volatile nature of renewable energy generation, leading to forecast errors with increasing renewable energy penetration. To minimise generation costs while meeting consumer electricity demands, it becomes imperative to strategically determine the output capacities and commitment statuses of each power plant, accounting for forecasts and constraints. Subsequently, a secondary dispatch round becomes indispensable to reconcile disparities between actual and forecasted demand [7][8].
In essence, addressing power system optimisation poses a significant challenge due to its capacity to encompass diverse technical, economic intricacies, and uncertainties. Attempting to incorporate all these elements into a single problem becomes computationally infeasible. Consequently, a typical approach involves considering a hierarchy of problems that tackle various time scales and perspectives. These include short-term dispatch (spanning a few days or weeks), mid-term operational planning (covering 1-2 years), and long-term operational planning (spanning 3-5 years). The outcomes derived from a long-term model can then be integrated into a model with a shorter horizon, but with more detailed considerations in other aspects of modelling [9].
Dynamic programming (DP) proves to be a potent method for addressing optimisation issues characterised by overlapping subproblems and optimal substructure. The enhanced forms of Dynamic Programming, namely Differential (or incremental) Dynamic Programming (DDP) and Discrete Differential Dynamic Programming (DDDP), contribute to achieving more precise outcomes in problems with overlapping subproblems and optimal substructure [10][11]. Various literature sources propose distinct methodologies; however, this discussion focuses on several successful and applicable approaches. One of these mathematical techniques were employed by a Sri Lankan Engineer, P.J. Perera, in the 1980s to establish a methodology for water resources planning [12].
P. J. Perera’s model of Dynamic Programming Techniques for Water Resource Planning
Dynamic Programming (DP) techniques which are based on Bellman’s Principle of Optimality might be conveniently used in water resources optimisation due to stage-by-stage applicability. The basic idea of Bellman’s Principle states that in a decision policy, whatever the initial states and the initial decisions are, the remaining decisions must constitute an optimal policy regarding the states resulting from the initial decision [12].
Basic derivation for the deterministic case
C=Optimize f[R_1 (S_1,U_0 ),…………………….,R_N (S_N,U_(N-1))]
Where,
C – Objective function
U_n∈ W_n ,S_n ∈ V_n
f – function representing the combined effect of the system transformations
succeeding from an initial state vector S_0 to a final state vector S_N resulting from a
series of decision vectors U_o,U_1,……………….,U_(N-1)
R_N (S_N,U_(N-1)) – the return at stage n with respect to the state vector S_n and the decision vector U_(N-1)
V_n,W_n – The administrative domains for decision vector U_n
S_n – State vector for stage n respectively
A cost minimisation, additive properties prevail in the function can be expressed as,
f[R_1 (S_1,U_0 ),…………………….,R_N (S_N,U_(N-1) )]= ∑_(n=1)^N▒〖R_n (S_n,U_(n-1) ) 〗
Where,
S_n – reservoir storage at stage n and
U_(N-1) – release of this storage at stage n
F_N (S_N ) – Minimum Cost for operating the system from an initial state S_0 to some final state S_N
F_N (S_N )=min∑_(n=1)^N▒〖R_n (S_n,U_(n-1) ) 〗
F_N (S_N )=min[R_N (S_N,U_(N-1) )]+min[∑_(n=1)^(N-1)▒〖R_n (S_n,U_(n-1) ) 〗]
By the definition,
F_(N-1) (S_(N-1) )= min[∑_(n=1)^(N-1)▒〖R_n (S_n,U_(n-1) ) 〗]
∴F_N (S_N )= min[R_N (S_N,U_(N-1) )+ F_(N-1) (S_(N-1) )]
Equation (1) represents the Dynamic Programming recursive algorithm, from an initial state vector S_0, to a final state S_N. Depending on whether the initial state or the final state is known, the algorithm could be applied in a forward or backward manner. If both states are known, it could be applied either way. Fig. 01 gives a graphical representation of the DP operation, when applied to a reservoir control problem.
Fig 01. Dynamic Programming Operation to Control the Reservoir Problem
Basic derivation for stochastic case
The stochastic element in reservoir operation is streamflow, and once this effect is introduced, the return function becomes dependent upon both storage 〖(S〗_n) and inflow (I_n), in addition to the decision variable U_(n-1).
Equation (1) can be revised as,
F_N (S_N )=min[∫_(-∝)^(+∝)▒〖P(I_N ) R_N (S_N,I_N,U_(N-1) )+ P_(N-1) (S_(N-1))〗]
Where,
P(I_N ) – The probability density function of I_N for a discrete probability distribution of I_N
F_N (S_N )=min[∑_(i=1)^K▒〖P(I_N^1 ).R_N (S_N,I_N^1,U_(N-1) )+ F_(N-1) (S_(N-1))〗]
Where,
I_N^1 – i^th discretisation of the streamflow at stage N
P(I_N^1 ) – Probability of occurrence
k – Total number of discretisation
Markov Chains and Serial Correlation
A Markov Chain of the first order is defined as a sequence of events x_0,x_1,x_2,………..,x_(n+1) of discrete random variables, with the property that the conditional distribution of x_(n+1), given x_0,x_1,x_2,………..,x_n , depends only on the value of x_n, but not further x_0,x_1,x_2,………..,x_(n-1). For any state values of h,j,…….,k belonging to the discrete state space,
Prob.[X_(n+1 )=k | X_0=h,………,X_n=j]=Prob.[X_(n+1)=k | X_n=j]
When a streamflow of a certain month is independently analysed using recorded previous data, a certain probability distribution could be derived for its representation. However, this distribution is usually conditional upon the actual streamflow of the previous month, as in the case of the occurrence in a first order Markov Chain. The basic derivation for the stochastic case with the introduction of this serial correlation is as follows.
F_N (S_N,I_(N+1)^1 )=min[∑_(i=1)^K▒〖P(I_N^1 I_(N+1)^1 〗).{R_N (S_N,I_N^1,U_(N-1) )+ F_(N-1) (S_(N-1),I_N^1)}]
Where,
P(I_N^1 I_(N+1)^1) – The conditional probability of streamflow at stage N falling into the
discretisation I, given I_(N+1)^1 inflow at stage N+1
It is evident that no unique optimal trajectory exists when serial correlation is considered.
Different Dynamic Programming Techniques
Differential Dynamic Programming (DDP)
Differential (or incremental) Dynamic Programming is a successive approximation technique, where the optimality of the objective function is considered locally, unlike in the conventional DP technique, where it is considered globally. DDP applies the principle of optimality in the neighbourhood of a nominal, possibly non-optimal trajectory. It is in effect, limiting the admissible domain of U_n and S_n to W_n and V_n respectively, which define the immediate neighborhood of U_n and S_n. These domains are successively updated in an iterative sequence to approach at an optimal solution.
2.2.2. Discrete Differential Dynamic Programming (DDDP)
Discrete Differential Dynamic Programming is an extension of DDP, in which the recursive equation of DP is used to search for an improved trajectory among discrete states within a neighbourhood (corridor) of a trial trajectory instead of the whole domain as in the case of conventional DP.
Let,
Δ_j S_n (k),j=1,2,…………….,I;k=1,2,…..,M;n=1,2,…….,N
be the j^th increment of the state variable S_n (k) at stage n in an Mdimentional vector space with a total of I possible increments. In the DDDP algorithm Δ_j S_n (k) can take any one value for j from 1 to I, from a set of assumed incremental values of the state domain.
Each of the M components of Δ_j S_n (k) can have I values, making the total number of combinations of Δ_j S_n (k) vectors at stage n equal to I. To have an accurate result, Δ_j S_n (k) must be sufficiently small. If the domains of the state variables span up to their entire extent, I will be very large, which is the main setback of conventional DP. In DDDP, the value of I is limited to 3 or 4 only. Defining a subdomain (D ̅_n) of the main domain such that,
D ̅_n= S ̅_n+ Δ ̅_j S ̅_n ,j=1,2,……,I
Where,
S ̅_n – The state vector at stage n represents a trial trajectory
All subdomains D ̅_n, n=1,2,…..,N together are called the corridor around the trial trajectory.
Procedure of DDDP
- Select a trial trajectory as close as possible to the optimal trajectory. This may be done with some insight into the actual operation of the system. It will ensure convergence on the global optimum rather than on a local optimum.
- Open corridors on state variables around this trial trajectory.
- Obtain the updated trajectory using the conventional DP technique within the corridor.
- Repeat steps (2) and (3) successively until the updated trajectory remains the same as the previous trajectory. They will be the optimum trajectory for that discretisation of the state variables within the corridor.
- If added accuracy is required, narrow down the corridor by having the same number of discretisation, but with a smaller increment of the state variable. Repeat steps (2), (3) and (4).
- If the accuracy is not sufficient, repeat step (5)
The return function in the problem under study is the thermal operating cost of the system. This is elevated as a function of the control variable which is equivalent to the hydro generations from the various plants. It is easily seen that the cost functions resulting from these generations are convex. The objective function, which is the sum of these convex functions is also convex, thus eliminating multiple optima. The optimal trajectory, therefore, becomes independent of the initial trial trajectory and this will be ensured by DDDP by selecting any arbitrary initial trajectory.
Dynamic Programming with successive approximations (DPSA)
DPSA can only be used where a unique optimal trajectory exists for the state variables. Since this is the basic requirement for DDDP too, both methods can be combined to give a fast-optimal result. In DPSA, the state variables are optimised one by one in a cyclic process, and the state variables which are not under optimisation at a particular time instant are frozen at their previous semi-optimal values. It can be proven as in DDDP, that if the return functions are convex, the result will converge to the global optimum.
Data Analysis
The main component of the stochasticity introduced to the problem is that due to inflows. A certain amount of statistical data processing is required to represent the inflows in an acceptable form to be input in the DP routine. It has been observed that streamflow patterns confirm to either normal or log-normal distributions with a high goodness- of- fit index. The log- normal distribution is useful in representing the skewness of the distribution, as is prominent in the case of daily inflow patterns whereas in the monthly patterns this skewness is damped to a certain degree. A chi-squared goodness of fit test will determine the suitability of the selected distribution.
When the proper distributions are obtained, discretisation is done with equal probability intervals. Five inflow conditions have been used in the study and hence each has a 20% probability of occurrence. Once discretisation is done, it is convenient to available the discrete conditional probabilities for the transition probabilities (or the transition probabilities) rather than the conditional distribution when treating the flow pattern of the year as events of a first order Markov Chain.
If I is the occurrence of the streamflow in month n within the class interval p,
Prob.(I_n^p│I_(n-1)^q ).m_p/m_q
Where,
m_q – number of occurrences of the streamflow in month n-1 within the class interval q
m_p – within the set of occurrences of m_q, the subset of occurrences of streamflow in month
n within the interval p
For k inflow intervals, the expectation of the inflow I_n for a certain month n, given the actual inflow interval of the previous month as q is,
I_n= ∑_(i=1)^k▒I_n^1 .prob(I_n^1 I_(n-1)^q)
Which may be quite different from the unconditional expectation,
I_n^0= 1/k.∑_(i=1)^k▒I_n^1
for equal probability intervals.
A transition probability matrix can be defined for the streamflow of any two adjoining months n and n+1 such that,
M= [P_ij ] ,i=1,2,…..,k ;j=1,2,…..,k
Where,
P_ij=Prob(I_(n+1)^j |I_n^i)
The Mid-Term schedule with deterministic inflows
This is implemented either on a monthly or weekly basis. The inflow regimes at various dam sites are found for various hydro-conditions using probability analysis and are classified under various dry-conditions. For example, 60%- dry would mean that the inflow will be greater or equal to the corresponding value for 60% of the time. It is illustrated by F in the inflow- duration curve given in Fig. 02.
Fig 02. Flow Duration Curve
The input of inflows to the model can be given as a combination of these dry conditions for 12 months. The model uses DDDPSA to arrive at the optimal trajectories for the reservoirs.
Evaluation of the return function
To obtain the thermal plant costs for each trajectory of the state variables, the load duration curve for that month is probabilistically simulated with the corresponding hydro energy included. The hydro components groups into the two basins of Mahaweli and Kelani rivers are fitted using Jacoby’s method with zero outages are assumed. For thermal plants, forced and scheduled outages are considered, and the plants are probabilistically stacked on the load duration curve in the order of their increasing operating cost.
Irrigation requirements
Irrigation only exists in the Mahaweli basin, at Polgolla and Minipe. The irrigation requirements off Polgolla are specified at three sites of Bowatenna, Elahera and Angammedilla. These requirements are given monthly at these sites. If an irrigation deficit is to be affected at Polgolla in an annual basis, there are numerous combinations possible in the adjusted diversions at these sites, and it should be initially investigated as to which distribution of the actual diversions are more realistic in actual practice. The DP routine can be made sensitive to irrigation diversions by the introduction of a penalty cost in the case of an irrigation deficit. The form of this penalty should incorporate not only the limitation of the total diversion deficits, but also how these deficits are distributed within the year.
For example, a total annual deficit of 100 MCM at Polgolla can give rise to a deficit at Bowatenna of 50 MCM in January and 50 MCM in December (leaving aside Elehera and Amgammedilla) or a 25 MCM deficit in January and a 75 MCM deficit in December. Since the diversion requirements at Bowatenna are 75 MCM for all months, the former case seems more appropriate in a practical sense as the optimising model does not consider the diversion storages as state variables. The introduction of a linear irrigation penalty can give rise to inappropriate monthly irrigation diversions because the results can exhibit insensitivity to the diversion pattern within the year. What is essentially required is the capability of the model to meet the requirements within the constraints imposed by the inflows, such that if there should be any deficits, these to be allocated monthly rather than in an annual basis. This is accomplished by introducing an exponential irrigational penalty rather than a linear penalty. For the example considered earlier, a linear penalty will give rise to identical irrigation penalty costs for the two cases, whereas the introduction of an exponential penalty of the form,
cost=EXP(deficit each month)
will minimise the penalty for the former case. The following expression illustrates this,
EXP(50)+EXP (50)=EXP(25)+EXP(75)
even though,
50 + 50 = 25 + 75
This penalty cost of course will have no physical meaning, but when added to the thermal costs for the various reservoir trajectories in the DP exercise, the resulting optimal trajectories will satisfy the irrigation in the most appropriate way monthly. This exponential irrigation penalty factor can also be used to enforce the trade-off between irrigation and power. A very low penalty will enforce negligible influence on irrigation deficits and in effect will optimise for the power only, whereas with a high penalty, the vice-versa will occur. Since quantification of actual irrigation losses are quite complex and inaccurate, under normal running the model is made to assume a very high penalty factor for irrigation deficits whereby the thermal cost is minimised only after meeting the irrigation requirements within the constraints of the inflows.
The mid-term schedule with stochastic inflows
A composite basin representation is used in the study due to the excessive number of state variables existent in the system. The composite representation agglomerates the reservoirs in each basin into one equivalent energy storage entity and an associated generating plant. The composite model is based on a single measure of potential energy which is indicative of the system’s generating capability. The one-dam representation of the multi-reservoir system in effect receives, stores and releases potential energy when operative in the DP routine. The inflows also should be represented by the potential energy input into the system. The dependence between the potential energy outflow and the electrical energy generated is given by the generation function of each composite basin [4].
The composite model is most applicable when the sequence of monthly decisions on the total hydro generation is of greater economic significance than the allocation of these among the various hydro plants. The decomposition of this total generation to those of the individual plants depends on the initial derivation of the generation function. Therefore, it is of basic importance to maintain as far as possible a consistent approach to this problem with that applied for achieving the main objective of the study.
Generation Functions
The power generation of a hydro plant as a function of the released water depends on the plant head and the overall plant efficiency. This is given by,
P=e.Q.H (2)
With,
P< P_max=e.Q_max.M_max
Where,
P – Power generated
Q – water flowrate through the turbine
H – Plant head
e -plant efficiency
Q_max is determined by the maximum turbine capacity, and H_max occurs when the reservoir level is at its maximum. If Q> Q_max, there will be no further contribution to the generation.
Equation (2) can be considered as the generation function for a single plant in the composite basin. The generation function of the composite model is an aggregate of the generation functions and spill characteristics of many similar plants. Its characteristics will depend on how the individual plants are operated integrally within the basin.
Operational curves for the generation function
The generation function for each basin and its inflow energy are obtained by a system simulation in a specific set of operational curves. In establishing these curves, it is necessary to consider that the optimisation routine within each basin must aim for similar goals as in the final scheduling programme, which is the system cost minimisation and serial correlation of the monthly inflows cannot be considered as a unique trajectory, which is essential for the simulation of the plants cannot be obtained that way.
With the above considerations in view, the most appropriate way of obtaining these curves would be a system cost minimisation using DDDPSA using the actual inflows of the available inflow regime. The system demands of the actual year of interest can be used as being incurring repetitively for every year for the available inflow period.
The short-term model with deterministic inflows
After the development of the mid-term models, a deterministic short-term model has been developed to operate daily. This is in fact an optimal dispatch model, to be used to recommend to the operators the optimal way to meet the demand on an hourly basis. Unlike in the mid-term case, the water transport time delays must be considered here as well as the true operational characteristics of the hydro and thermal machines. The operation of the Kotmale and Mahaweli complex ponds become significant in this case, and they must be introduced as state variable along with Ukuwela and Bowatenna. Victoria can be hydraulically decoupled from the system as the discharges from Polgolla would not make any significant changes in the Victoria storage when operating on an hourly basis.
The optimisation routine used here is DDDPSA. However, a major modification had to introduced into the conventional DDDPSA technique to fit it to the short-term problem. In normal DDDPSA technique to fit it to the short-term problem. In normal DDDPSA the state variables which are not under optimisation are frozen at their semi-optimal values. When relating this concept to the actual case under study the following problem arises,
Fig 03. Cascadic structure of Laxapana Hydro Complex in Sri Lanka
Consider the cascade of Castlereigh, Norton and Laxapana reservoirs (Fig. 03). The Wimalasurendra power plant has designed for peak operation whereas the Old Laxapana plant is running practically at baseload. When using conventional DDDPSA to these two plants, first the Castlereigh reservoir is optimised using the energy of the Wimalasurendra power plant while the Norton Pond trajectory is frozen. This means that al the water discharged through the Wimalasurendra power plant is forced through the Old Laxapana plant as well, thereby disabling the Norton Pond to act as a buffer to accommodate the highly irregular discharges of Wimalasurendra plant. Hence the feasible space of trajectories of the Norton Pond is severely curtailed due to the operation of the Old Laxapana plant. In the modified DDDPSA technique developed, the plant generations are frozen rather than the reservoir storages. In the example considered above when optimisation is done on the Castlereigh reservoir. The Old Laxapana generation is frozen, rather than the Norton Pond storage so that limitations on the Castlereigh discharge are not imposed by the operation of the Old Laxapana plant. After the optimisation of the Castlereigh reservoir, the Norton Pond trajectory is updated with the optimised discharges of the Wimalasurendra plant and used in the subsequent optimisation.
Analysis of METRO Model
METRO is applied for operations planning of integrated hydro, thermal, irrigation system. This uses the water value concept in deriving an operation policy for the power reservoirs of the system and uses optimal reservoir balancing factors to balance the reservoir drawdown (maintaining an optimal mix) on simulation of the operating policy. Water value is the incremental long-term replacement cost of hydro in storage by thermal and is a well proven concept used in medium term operational planning [13] [14] [15].
METRO can optimise the operation of hydro power resources in the system to minimise thermal operation costs and expenses associated with unserved energy in the long and medium term, considering varying hydrological conditions. It is used to identify thermal switching curves for the base load operation of thermal plants specifically designed for dry hydrological conditions. Moreover, it integrates the complete irrigation system in conjunction with the power station. Additionally, it excels in evaluating the opportunity cost of water in storage (water value) by assessing incremental thermal costs. Furthermore, it establishes reservoir operating policies (rule curves) that delineate reservoir drawdown philosophy based on reservoir storage and time, accounting for diverse hydrological conditions, and incorporates maintenance schedules for all hydro and thermal plants. Moving forward, METRO can optimise maintenance planning to minimise the impact of planned forced outages, treating them as truly random events rather than relying on capacity de-rating methods. It provides daily operational guidelines in the form of recommended generation schedules for all plants, projects expected system performance for one year ahead by simulating the past hydrological sequence and fuel cost budgeting, conducts quantitative and qualitative analyses of power system reliability, analyses irrigation reliability, examines power system and irrigation trade-offs, and incorporates all forthcoming expansions of the system through data file updates.
The precision of the short-term dispatch model is crucial, particularly since short-term economic dispatch can extend from 24 to 168 hours. During this timeframe, the variability in water inflow and fuel prices does not significantly impact the accuracy of the short-term model’s output. However, the accuracy of the model can be influenced by factors such as Addition and retirements of power plants, Varied shapes of load duration curves and Differing unit prices of thermal power plants.
Mathematical modelling of Economic Dispatch
The subsequent notations are employed in the formulation of the dispatch algorithm model.
t Time interval (hour) index
T Total number of time intervals (scheduling horizon)
D_t Total Demand (MW)
D_t^Th Total Thermal Power Supply
D_t^Hy Total Hydro Power Supply
I Number of available thermal plants
N_i Number of units in ith thermal plant
M_j Number of units in jth thermal plant
u_it Commitment state of ith thermal plant,
〖if u〗_it=1 committed; u_it=0 reserved
S_nit Status of ith thermal plant, nth unit at tth state
S_nit^D Discrete status of ith thermal plant, nth unit
if 0 <S_(ni,t-1)≤1 then S_nit^D=1
if -1 <S_(ni,t-1)≤0 then S_nit^D=0
〖Mt〗_nitMaintenance state of ith thermal plant, nth unit
if 〖Mt〗_nit=1 ,available; 〖Mt〗_nit=0,offline
〖Pt〗_nit Power output of ith thermal plant, nth unit (MW)
〖Ph〗_jt Power output for the reservoir j during the t hour
C_i (..) Cost function of ith thermal plant in LKR/hr; a function of P_it
〖Pt〗_ni^max Maximum power output of ith thermal plant, nth unit (MW)
〖Pt〗_ni^min Minimum power output ith thermal plant, nth unit (MW)
〖Ph〗_j^max Maximum power output of the hydro unit j
〖Ph〗_j^min Minimum power output of the hydro unit j
∆P_ni Ramp rate of ith thermal plant, nth unit (MW/hr)
T_ni^up Minimum up time of ith thermal plant, nth unit (hrs)
T_ni^dn Minimum down time ith thermal plant, nth unit (hrs)
T_ni^SD Reserved hours of ith thermal plant, nth unit (hrs)
St〗_nit Start-up cost of ith thermal plant, nth unit (LKR)
Sd〗_ni Shutdown cost of ith thermal plant, nth unit (LKR)
C_t (〖Pt〗_nit) Total Thermal Cost of ith thermal plant, tth hour
Dc〗_(1,t)^(l-1) Total demand when l-1 units are dispatched in tth hour
〖Rd〗_(l+1,t)^l Remaining demand when I-(l+1) units to be dispatched in tth hour
x_nit Consecutive time that ith thermal plant, nth unit has been down at time t
α Maximum capacity limit percentage
Objective function
The goal of resolving the economic dispatch problem in an electric power system is to ascertain the generation levels for all active units that both minimise the overall cost and adhere to a defined set of constraints.
Total cost=min∑_(t=0)^T▒〖{∑_(i=0)^(l-1)▒∑_(n=1)^Ni▒〖|S_nit | 〖Mt〗_nit (a_i 〖Pt〗_nit+c_i )+〖St〗_nit S_nit^D (1- S_(ni,t-1)^D ) + 〖 Sd〗_nit S_(ni,t-1)^D (1-S_nit^D ) 〗}〗
System Constrains
Several constraints of the system influence the unit commitment and economic dispatch.
Demand satisfaction for each hour t
Total Demand=Hydro generation+Thermal Generation+Losses
D_t= D_t^Hy+ D_t^Th+Losses
D_t= ∑_(j=1)^J▒〖〖Ph〗_jt+ ∑_(i=1)^I▒∑_(n=0)^Ni▒〖|S_nit | 〖Mt〗_nit 〖Pt〗_nit+ Losses〗〗
Technical operation limits of each unit
〖Pt〗_ni^min < 〖Pt〗_nit < 〖Pt〗_ni^max ∀t∈ [1,T],∀i∈[1,I],∀n∈[1,Ni]
〖Ph〗_j^min < 〖Ph〗_j < 〖Ph〗_j^max ∀t∈ [1,T],∀j∈[1,J]
Startup/ Shutdown Cost
〖St〗_nit ={█(〖St〗_ni^c if x_nit> t_cold@〖St〗_ni^w if t_hot<x_nit< t_cold@〖St〗_ni^h Otherwise )┤
startup cost= 〖St〗_nit S_nit^D (1- S_(ni,t-1)^D )
shutdown cost= 〖Sd〗_nit S_nit^D (1- S_(ni,t-1)^D)
Spinning reserve requirement
Spining reserve=Total amount of generation available from all synchrounous units-(present loads+losses)
∑_(i=1)^I▒∑_(n=0)^Ni▒〖|S_nit |(〖Pt〗_ni^max- 〖Pt〗_nit )+ ∑_(j=1)^J▒〖(〖Ph〗_j^max- 〖Ph〗_jt )≥ βt〗〗
For Sri Lankan context, β= 0
Every hydropower system exhibits distinctive hydraulic limitations, primarily dictated by geographical and hydrological factors. Occasionally, the release of water from one reservoir can impact the availability of water in another reservoir, referring to them as hydraulically coupled units.
Minimum uptime and downtime
Combined cycle power plants, coal power plants, and nuclear power plants exhibit variable minimum up times, contingent on the duration of the plant shutdown. In contrast, other thermal plants have distinct minimum up and down times.
Minimum startup time of a thermal power plant:
if (S_(ni,t-1)=0 and S_nit>0) then S_(ni,t+τ)=1 ; τ ≤ T_ni^up
T_ni^up= {█(T_ni^upC if x_nit> t_cold@T_ni^upW if t_warm< x_nit< t_cold@T_ni^upH otherwise)┤
Minimum downtime of a thermal power plant:
if (S_(ni,t-1)=1 and S_nit<0 )then S_(ni,t+τ)=0 ; τ ≤ T_ni^dn
〖Pt〗_nit= {█(mτ+K if (τ<T_ni^up)@〖Pt〗_nit otherwise)┤
Hydropower plants are characterised by rapid startup and shutdown capabilities. Therefore, the minimum startup and shutdown times are not constraints for hydropower plants.
Ramp rate
∆P_in=|〖Pt〗_nit- 〖Pt〗_(ni,t+1) |
It is assumed that rate of decrease in power is same as the rate of increase in power.
Status of the thermal plant
S_nit= {█(〖Pt〗_nit/(〖Pt〗_ni^min ) if (τ< T_ni^up)@1 if (τ≥ T_ni^up))┤
S_nit= {█(-〖Pt〗_nit/(〖Pt〗_ni^min ) if (τ< T_ni^dn)@ 0 if (τ≥ T_ni^dn))┤
Maximum Capacity Unit
αD_t ≥ 〖Pt〗_nit and αD_t ≥ 〖Ph〗_nit; α=0.2
Thermal Plant Cost Function
C_t (〖Pt〗_nit )= 〖aPt〗_nit+ 〖bPt〗_nit^2+c
Sri Lankan system does not contain with second order of cost function.
∴ C_t (〖Pt〗_nit )= 〖aPt〗_nit+ c
Algorithm of Metro for economic dispatch of Thermal Power plants
Input
- Half an hour demand curve of thermal and total demand
- Available plants and cost details
- Maintenance schedule
- Initial status of the plants (ON/OFF, Number of Units running, Number of hours stopped)
Step 1
- Generate initial schedule.
- Sort the power plants ascending order to unit cost.
- Combined cycle power plant consists of no-load cost per hour. Therefore, the average cost of the plant is taken as (a〖Pt〗_ni^max-C)/(〖Pt〗_ni^max )
Step 2
- Generate option for the next schedule (Number of Options = 2^I)
- Remove infeasible options by considering total minimum load and total maximum load (total minimum load < Load < total maximum load)3
∑_(i=1)^I▒∑_(n=1)^Ni▒〖u_it 〖Mt〗_nit 〖Pt〗_ni^min< D_t^Th< ∑_(i=1)^I▒∑_(n=1)^Ni▒〖u_it 〖Mt〗_nit 〖Pt〗_ni^max 〗〗
Step 3
- Options are generated for the plants but not for the plants’ units. Therefore, if a plant is committed that does not mean that all the available units are up running. But if a plant is reserved no units will be running.
- Even though it is a merit order list (MOL) all the committed power plant of an option should be dispatched regardless of the cost. Therefore, at least its operational minimum power will be dispatched.
- Plant Capacity ≤(Demand-Dispatched Capacity- Minimum capacity of total remaining power plants)
- Dispatched Capaticy= ∑_(i=1)^(l-1)▒∑_(n=1)^Ni▒〖|S_nit |〖Mt〗_nit 〖Pt〗_nit 〗
- Total mimimun load of remaining plants= ∑_(i=1+1)^I▒∑_(n=1)^Ni▒〖|S_nit |〖Mt〗_nit 〗 〖Pt〗_ni^min
- ∑_(n=1)^Ni▒〖|S_nit | 〖Mt〗_nit 〖Pt〗_nit≤(D_t^Th- ∑_(i=1)^(l-1)▒∑_(n=1)^Ni▒〖|S_nit | 〖Mt〗_nit 〖Pt〗_nit 〗- 〗 ∑_(i=1+1)^I▒∑_(n=1)^Ni▒〖|S_nit |〖Mt〗_nit 〗 〖Pt〗_ni^min)
- Four transitions can occur in a step change.
Transition 1
if u_(i,t-1)=0 and S_(ni,t-1)=0 then
if u_it=1 then
if (〖Pt〗_(nit )≤(D_t^Th-〖Dc〗_(0,t)^(l-1)-〖Rd〗_(l+1,t)^I )) then
if (τ<T_ni^up )then,0 < S_nit <1
if (τ ≥ T_ni^up )then,S_nit=1
if (〖Pt〗_nit> (D_t^Th-〖Dc〗_(0,t)^(l-1)-〖Rd〗_(l+1,t)^I ))then
S_nit=0
if u_it=0,then S_nit=0
Transition 2
if u_(i,t-1)=1 and S_(ni,t-1)=1 then
if u_it=1 then
if (〖Pt〗_(nit )≤(D_t^Th-〖Dc〗_(0,t)^(l-1)-〖Rd〗_(l+1,t)^I ))then,S_nit=1
if (〖Pt〗_(nit )>(D_t^Th-〖Dc〗_(0,t)^(l-1)-〖Rd〗_(l+1,t)^I )) then
if (τ<T_ni^dn )then,-1< S_nit <0
if (τ ≥ T_ni^dn )then,S_nit=0
if u_it=0,then
if (τ<T_ni^dn )then,-1< S_nit <0
if (τ ≥ T_ni^dn )then,S_nit=0
Transition 3
if u_(i,t-1)=1 and 0< S_(ni,t-1)<1 then u_nit ≠1
if u_nit=1 then,
if (τ< T_ni^up ) then,0<S_nit<1
if (τ ≥ T_ni^up ) then S_nit=1
Transition 4
if u_(i,t-1)=0 and-1<S_(ni,t-1)<0 then u_nit ≠1
if u_(nit )= 0 then
if (τ< T_ni^dn ) then-1<S_nit<0
if (τ ≥ T_ni^dn ) then,S_nit=0
Total cost of one option consists of energy cost and transition cost.
Total cost for hour=Energy Cost+Transition Cost
Transition Cost= 〖St〗_nit S_nit^D (1-S_(ni,t-1)^D )+ 〖Sd〗_nit S_(ni,t-1)^D (1- S_nit^D)
Energy Cost= ∑_(i=0)^(l-1)▒∑_(n=1)^Ni▒〖|S_nit |〖Mt〗_nit (a_i 〖Pt〗_nit+ c_i)〗
∴ Total cost for hour= ∑_(i=0)^(l-1)▒∑_(n=1)^Ni▒〖|S_nit |〖Mt〗_nit (a_i 〖Pt〗_nit+ c_i)〗+〖St〗_nit S_nit^D (1-S_(ni,t-1)^D )+ 〖Sd〗_nit S_(ni,t-1)^D (1- S_nit^D)
After completion of scheduling one option,
Full Total Cost=Previous Total Cost+Total Cost for hour
Short term power system planning
There are several methods of short-term forecasting including,
Classical time series and regression method
Time series models
Artificial intelligence and computational intelligence methods
Hybrid approaches
Absolute mean percentage error (MAPE) is used to measure the performance of these models.
MAPE= 100/T ∑_(t=1)^T▒〖(y_t-f_t)/y_t (%)〗
Where,
y_t= real value at time t
f_t= forecast at time t of period T
3.3.1 Water value models
Cost of fuel is the main component that is deciding the operation cost of a thermal power plant.
∴thermal plants are representing according to the operating cost∶ c_j ($/MWh)
Where, j=1,2,…….,n
Generation Limits,
g_t (j)≤G(j)
Where,
g_t (j)= energy production of plant j in stage t (MWh)
G(j) = maximum generation capacity of plant j
3.3.2 Determining the power plant combination to achieve the required load with minimising the fuel cost
z_t=min ∑_(j=1)^J▒〖c_j g_t (j)〗
∑_(j=1)^J▒g_t (t)= d_t
Where,
z(t) – the system operating cost in stage t
c – unit operation costs
d(t) – system load
g(t) – Power production
g – generation capacities
Immediate Cost Function (ICF) and Future Cost Function (FCF)
Hydropower plants utilise their stored water resources to avoid the fuel expenses of thermal power plant units. The availability of hydropower plants depends on the water storage of reservoirs. The main decision to be made at the scheduling stage is to release the required amount of water by evaluating immediate financial gain and the expected value of water in the future. The expected value of the water in the future is represented by a function of inflow, reservoir level, and time. The ICF (Inflow, Cost, Function) has a relationship with the cost of thermal energy at the t stage. When the final storage is gradually increasing, only a small amount of water is available for electricity generation at the same stage. As a result, more thermal energy is required, leading to an increase in cost. By the way, the FCF (Future Cash Flow) associated with predicted thermal expenses from (t+1) to the ending point of the period will decrease due to the availability of more water for future energy production, as shown in Fig. 04 [16] [17].
Fig. 04 Immediate and Future Costs Versus Final Storage
Optimisation Algorithm
p_k (k+1)= α_p (k) p_n (k)+ϵ(k)
Where,
p_k – average short price in the week
p_n (k) – normalised value of p_k
α_p (k) – auto correction coefficient determined from the future function
ϵ(k) – normal random variable
Principle of mass production Theory
When the utilised capacity is increasing, fixed cost is declined and increasing the variable.
tc= FC/x+vc
tc – total average cost
FC – total fixed cost
X – production quantity
vc – total variable cost
To obtain the variable cost of a thermal power plant, specific heat rate curve is used. Specific heat curve provides information on the necessary heat energy required for generating one unit of electrical energy H(P)
H(P)= α/p+ β+γP [kJ/kWh]
F(P)=P*H(P) [MJ/h]
C(P)=k*F(P)=a+bP+cP^2 [$/h]
P – Thermo power plant production
α,β,γ – coefficient
F(P) – input – output curve
C(P) – fuel cost function
k – fuel cost coefficient
Water value deterministic function
f_c (P){█(R_l (t)= R_pl (t),f_c (P)=0@R_l (t)< R_pl (t),f_c (P)=A@R_l (t)> R_pl (t),f_c (P)=B)┤
f_c (P) – function cost curve for water value
P – power needed in system
R_l – reservoir level
R_pl – planned reservoir level
System production optimisation [6]
min∑▒〖c(j) g_t (j)+ ∑▒〖ρ(k,v_(t+1) ) α_(t+1) (k) v_(t+1) (k)± E_i 〗〗
Where,
min∑▒〖c(j)g_t (j)〗 – immediate cost of thermal power plants Operating cost in stage t
∑▒〖ρ(k,v_(t+1) ) α_(t+1) (k) v_(t+1) (k) 〗 – future cost of hydro power plant
ρ(k,v_(t+1) ) – production coefficient of reservoir k
α_(t+1) (k) – water value function of the reservoir
v_(t+1) (k) – reservoir volume at the end stage t
E_i – power exchange
For the easiness of understanding, the function can be modified as:
min〖F= ∑_(i=1)^5▒〖c_i (TPP)+ ∑_(i=1)^6▒〖c_i (HPP)〗〗〗 ±c(PEX)
Where,
c_i (TPP) – cost of thermal power plant i
c_i (HPP) – cost of hydro power plant i
c(PEX) – energy sold or brought on power exchange
c_i (HPP)= ρ_i≈ V_(t+1)
Summation of power production and energy exchange has to be equal with the planned load for each period.
∑_(i=1)^5▒〖c_i (TPP)*P_i+ ∑_(i=1)^6▒〖C_i (HPP)* P_i ±c(PEX)*P=L〗〗
To achieve the maximum financial benefits,
C_i (HPP)+ C_i (TPP) ≤min〖[c(EPEX),c_i (TPP)]〗
CONCLUSION AND DISCUSSION
Cost accounting serves as a multifaceted method encompassing the accumulation, classification, summarisation, and interpretation of information crucial for operational planning, control, and strategic decision-making [18]. A primary aim of cost accounting lies in determining optimal selling prices. Beyond this, its objectives extend to enhancing control efficiency, facilitating financial statement preparations, and laying the groundwork for operational policies, such as delineating cost-volume-profit relationships and evaluating the feasibility of outsourcing from external suppliers. In the realm of electricity generation, these objectives align with tasks like tariff formulation, strategic decisions regarding the establishment of pump storage power plants, importing power from neighbouring countries via tie-lines, and electrification initiatives [19]. Broadly defined, cost signifies the expenditure associated with production, calculated either through actual or notional values. For power plants, unit cost considerations encompass operational and maintenance expenses, fuel costs, and depreciation. However, the valuation of energy units extends beyond these factors to include opportunity costs, reflecting the value in terms of the unit cost of the next best alternative [19] [20][21].
The primary objective of optimising the operational schedule of a hydrothermal system is to define efficient operational strategies for each planning stage and corresponding system states to achieve targeted power generation from each power plant. These strategies aim to minimise operational costs over time, including fuel expenses and penalties for load supply failures, while considering constraints such as stored water resources. The dynamic nature of water resource availability, influenced by present and future variations in value, underscores the importance of effectively managing reservoir levels [22][23]. This management involves linking operational decisions at each stage with their future implications, particularly regarding the depletion of hydroelectric energy and the potential need to resort to more costly thermal energy resources during periods of low inflow [24] [25]. Conversely, when reservoir levels reach capacity due to high inflow rates, it is economically advantageous to prioritise the use of inexpensive hydro sources over thermal sources. The inherent uncertainty of future inflow forecasts, owing to the stochastic nature of the problem, complicates decision-making. Moreover, the optimisation challenge is further compounded by the presence of multiple interconnected reservoirs and the need for multi-duration optimisation.
This research paper discusses two methodologies: one introduced by P.J. Perera and the other known as METRO. Both methodologies are employed to address practical problems, specifically in the realm of hydrothermal dispatching. In these systems, reservoirs of stored water are utilised. However, generators encounter an inventory problem. They seek to optimise the release of water from reservoirs to maximize profits using stochastic processes. Yet, this optimisation process is constrained by the unpredictability of future reservoir inflows.
While a substantial body of literature on economic dispatch focuses on enhancing mathematical computation methods to efficiently solve formulated objective functions, it has been noted that economic dispatch problems can be simplified by improving the quality of modelling.
CRediT authorship contribution statement.
D. Punsara Colambage: Formal analysis, Methodology, Writing – original draft.
W.D.Anura S. Wijayapala and Tilak Siyambalapitiya: Supervision, Writing – review & editing.
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Data availability
Data will be made available on request.
ACKNOWLEDGEMENT
The authors extend their sincere gratitude in memory of Prof. H.Y. Ranjit Perera, the pioneer of this research team, who has suddenly passed away.
During the preparation of this work the authors used Chatgpt in order to improve the language. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.
REFERENCES
- Dai H., Zhang N. and Su W. A, Literature Review of Stochastic Programming and Unit Commitment, Journal of Power and Energy Engineering, 2015, 3, 206-214, http://dx.doi.org/10.4236/jpee.2015.34029
- Colambage, D.P., Wijayapala, W.D.A.S., Siyambalapitiya, T. A Novel Approach of Exploring Water Value: A Study based on Sri Lanka. Engineer, Volume 57 issue 4, October 2024, pp 1-12. doi: 10.4038/engineer.v57i4.7660
- D.P. Colambage, W.D.A.S. Wijayapala, Tilak Siyambalapitiya. Water Value Concept in Generation Planning. Smart Grid Security and Protection, Chapter 41, 978-981-96-0823-2, Springer, May 2025. https://link.springer.com/book/9789819608232
- Colambage D. P., Wijayapala W. D. A. S and Siyambalapitiya T., Analysing the Impact of Diverse Factors on Electricity Generation Cost: Insights from Sri Lanka, Energy, Volume 311, 2024, 133451, ISSN 0360-5442, https://doi.org/10.1016/j.energy.2024.133451
- Pereira M.V.F., Pinto L.M.V.G. Multi-stage stochastic optimization applied to energy planning. Mathematical Programming 52, 359–375 (1991). https://doi.org/10.1007/BF01582895
- Dias B.H., Marcato A. L. M., Souza R. C., Soares M.P., Silva Jr I.C., Oliveira E. J. De, Brandi R.B.S. and Ramos T.P. , “Stochastic Dynamic Programming Applied to Hydrothermal Power Systems Operation Planning Based on the Convex Hull Algorithm”, Mathematical Problems in Engineering, vol. 2010, Article ID 390940, 20 pages, 2010. https://doi.org/10.1155/2010/390940
- Stüber M. and Odersky L., Uncertainty modeling with the open-source framework urbs, Energy Strategy Reviews, Volume 29, 2020, 100486, ISSN 2211-467X, https://doi.org/10.1016/j.esr.2020.100486.
- Masache S.P. and Barrera-Singaña C., “Short-Term Hydrothermal Economic Dispatch Applied on Hydraulic Coupled Power Plants Using Dynamic Programming,” 2020 IEEE ANDESCON, Quito, Ecuador, 2020, pp. 1-6, doi: 10.1109/ANDESCON50619.2020.9272182.
- Cerjan M., Marčić D. and Delimar M., “Short term power system planning with water value and energy trade optimisation,” 2011 8th International Conference on the European Energy Market (EEM), Zagreb, Croatia, 2011, pp. 269-274, doi: 10.1109/EEM.2011.5953022.
- Marín Cruz I., Badaoui M., Mota Palomino R., Medium-Term Hydrothermal Scheduling of the Infiernillo Reservoir Using Stochastic Dual Dynamic Programming (SDDP): A Case Study in Mexico. Energies 2023, 16, 6288. https://doi.org/10.3390/en16176288
- Golembiovsky, D., Pavlov, A. and Daniil, S. (2021) Experimental Study of Methods of Scenario Lattice Construction for Stochastic Dual Dynamic Programming. Open Journal of Optimization, 10, 47-60. doi: 10.4236/ojop.2021.102004.
- Perera, P.J. Dynamic Programming Techniques in Water Resource Planning. Hand-written document
- Andruszkiewicz J., Lorenc J. and Weychan A., Seasonal variability of price elasticity of demand of households using zonal tariffs and its impact on hourly load of the power system, Energy, Volume 196, 2020, 117175, ISSN 0360-5442, https://doi.org/10.1016/j.energy.2020.117175.
- Eid C., Koliou E., Valles M., Reneses J. and Hakvoort R., Time-based pricing and electricity demand response: Existing barriers and next steps, Utilities Policy, Volume 40, 2016, Pages 15-25, ISSN 0957-1787, https://doi.org/10.1016/j.jup.2016.04.001.
- Beltrán F., Finardi E.C., Oliveira W. de, Two-stage and multi-stage decompositions for the medium-term hydrothermal scheduling problem: A computational comparison of solution techniques, International Journal of Electrical Power & Energy Systems, Volume 127, 2021, 106659, ISSN 0142-0615, https://doi.org/10.1016/j.ijepes.2020.106659.
- Costa B.F.P. D. and Leclère V., Dual SDDP for risk-averse multistage stochastic programs, Operations Research Letters, Volume 51, Issue 3, 2023, Pages 332-337, ISSN 0167-6377, https://doi.org/10.1016/j.orl.2023.04.001.
- Chyong C.K. and Newbery D., A unit commitment and economic dispatch model of the GB electricity market – Formulation and application to hydro pumped storage, Energy Policy, Volume 170, 2022, 113213, ISSN 0301-4215, https://doi.org/10.1016/j.enpol.2022.113213.
- Horngren C. T., Datar S. M., Foster G., Cost Accounting: A Managerial Emphasis Charles T. Horngren series in accounting, ISSN 1559-2839 Prentice-Hall series in accounting, Pearson Prentice Hall, 2006
- Colambage D.P. and Perera H.Y.R., 2019. Assessment of cost of unserved energy for Sri Lankan industrial sector. IEEE International Conference on Industrial Electronics for Sustainable Energy Systems (IESES), Hamilton, New Zealand, 2018, 261-266, doi: 10.1109/IESES.2018.8349885.
- Colambage D.P. and Perera H.Y.R., Assessment of Cost of Unserved Energy for Sri Lankan Commercial Sector. Moratuwa Engineering Research Conference (MERCon), Moratuwa, Sri Lanka, 2018, 253-258, doi: 10.1109/MERCon.2018.8421892.
- Schledorn A., Guericke D., Andersen A. N. and Madsen H., Optimising block bids of district heating operators to the day-ahead electricity market using stochastic programming, Smart Energy, Volume 1, 2021, 100004, ISSN 2666-9552, https://doi.org/10.1016/j.segy.2021.100004.
- Lund H., Thellufsen J. Z., Østergaard P. A,, Sorknæs P., Skov I.R. and Mathiesen B. V., EnergyPLAN – Advanced analysis of smart energy systems, Smart Energy, Volume 1, 2021, 100007, ISSN 2666-9552, https://doi.org/10.1016/j.segy.2021.100007.
- Caldera U., Sadiqa A., Gulagi A. and Breyer C., Irrigation efficiency and renewable energy powered desalination as key components of Pakistan’s water management strategy, Smart Energy, Volume 4, 2021, 100052, ISSN 2666-9552, https://doi.org/10.1016/j.segy.2021.100052.
- Colambage, D.P. Economic Impact of Electricity Interruptions: A Case Study Based on Sri Lanka. American Journal of Electrical Power and Energy Systems, 11(1), 1-10. (2022). https://doi.org/10.11648/j.epes.20221101.11
- Colambage, D. P., Perera, H.Y. R. Variation of Cost of Unserved Energy Depending on Interruption Duration: A case study for Sri Lanka, IEEE 2nd International Conference on Power and Energy Applications (ICPEA), 27-30 April 2019, NTU- Singapore. doi: 10.1109/ICPEA.2019.8818530
