The ongoing load-shedding and energy crises due to mismanagement of energy produced by different sources in Pakistan and increasing dependency on those sources which produce energy using expensive fuels have contributed to rise in load shedding and price of energy per kilo watt hour. In this paper, we have presented the linear programming model of 95 energy production systems in Pakistan. An improved multiverse optimizer is implemented to generate a dataset of 100000 different solutions, which are suggesting to fulfill the overall demand of energy in the country ranging from 9587 MW to 27208 MW. We found that, if some of the power-generating systems are down due to some technical problems, still we can get our demand by following another solution from the dataset, which is partially utilizing the particular faulty power system. According to different case studies, taken in the present study, based on the reports about the electricity short falls been published in news from time to time, we have presented our solutions, respectively, for each case. It is interesting to note that it is easy to reduce the load shedding in the country, by following the solutions presented in our dataset. Graphical analysis is presented to further elaborate our findings. By comparing our results with state-of-the-art algorithms, it is interesting to note that an improved multiverse optimizer is better in getting solutions with lower power generation costs.

1. Introduction

Proper management and finding out optimal solutions to utilize the available power production sources is an important optimization problem in Pakistan, which is the main theme of this paper [111]. According to the reports in [12, 13], the demand of the power is increasing exponentially and which is causing an increase in load shedding each year [14]. The cause of this shortage is mismanagement of power production sources.

According to a report published by National Electric Power Regulatory Authority (NEPRA) in 2016 [15], the total capacity of 27240 MW is installed in the country, whereas the demand rises to 22000 MW in summer [16]. These production units can produce enough energy if all are managed optimally. In [17], these production units were normally producing up to 13240 MW with 4760 MW of shortfall. In nomenclature, we have indicated all the 95 installed capacities of electricity production with majority of them depending on thermal sources [15]. In [14], it is reported that the gap between demand and production is expected to increase in the coming years. It is further reported in the literature [14, 18] that urban areas are facing a load shading in summer season as of 10–12 hours, and in rural areas, it is even worse as 16–18 hours a day. Urbanization and rapid increase in population have caused an increase in the industrial zones [19], and it is worth noting that the demand of electricity was 16000 (MW) during summer of 2012 while the short fall was 8500 (MW) [20].

In this study, we have explored a scientific decision-making approach in terms of mathematical programming for providing a strategy to utilize the available production sources of Pakistan to meet the demand of energy in different seasons. However, in [21], a general mechanism was presented in terms of 5 major energy production types as to . They have fixed the per unit price to estimate the overall mathematical model. We have extended the model and included 95 production sources, which generates electricity using different means (oil, gas, wind, hydropower, and coal). We have implemented per unit price according to the actual price of fuel used to run these sources. The production capacities of each source are given in Table 1 [22]. Our model will help to reduce the shortfall by utilizing the power production sources properly.

In the last few years, evolutionary algorithms (EAs) are implemented to solve different real-world optimization problems. EAs simulate social behavior or natural procedures in order to handle optimization problems with best solution. Among the popular EAs are the Bat algorithm [23, 24], grasshopper optimization algorithm (GOA) [25], and firefly algorithm (FA) [26, 27]. The efficiency of metaheuristics is better as compared to classical optimization techniques in solving optimization problems with iterations and random search behavior. The main idea behind all EAs is the survival of the fittest, which in return increases the fitness of individuals in population. EAs are implemented with different frameworks, which includes population structure and search equations [28]. EAs with different frameworks generate random populations. Then, all solutions are evaluated according to a given objective function. The search space is explored and exploited in two phases. EAs get stuck in local optima due to their random search. Several papers are published in the literature to overcome this problem in EAs [29, 30].

The multiverse optimization algorithm (MVO) is a nature-inspired optimization technique given in [31]. The key idea of MVO is inspired from the theory of multiverses. Since its invention, it is applied to solve several real-world problems. In [31], it is shown that MVO obtained very competitive results compared with other metaheuristic optimization algorithms. Also, Faris et al. [32] employed the MVO for training the multilayer perceptions in neural network. The efficiency of their techniques was evaluated on different medical datasets from UCI repository. The outcome of their experiments is then compared to different metaheuristics, namely, Particle Swarm Optimization (PSO), Genetic Algorithm (GA), Differential Evolution (DE) Algorithm, and Cuckoo Search (CS). It was observed that MVO was improved in terms of convergence and avoidance of local optima. MVO still lacks behind in accuracy of results and convergence speed. Fewer studies, like Levy-flight-based MVO [33] and a quantum version of MVO [34], are recently proposed to further enhance the capabilities of MVO. Experimental results show that they have improved the solutions to some extent.

Currently, as can be seen in the former review, there are limited works on improving the performance of MVO. This work presents a new method of initialization with the MVO algorithm called improved multiverse optimization algorithm (IMVO). The experimental results on mathematical model of economic dispatch problem of electricity generation system in Pakistan and its three case studies show that IMVO is very competitive over FA, Bat algorithm, and GOA.

Furthermore, we have presented solutions to five case studies for proper management and production by all the 95 energy producers in the country. Our model can be implemented to reduce the load shedding and distribute the burden of production on respective power production units.

This work is organized as follows: Section 2 describes the concept of multiverse optimizer and improvements in it. The mathematical model for the problem is described in Section 3. Experimental settings, results, and discussions based on numerical simulations are given in Section 4. Conclusion of our findings is presented in Section 5.

2. An Improved Multiverse Optimizer

The MVO algorithm simulates the theory of multiverse, where three verses play main part in the whole algorithmic procedure: these are white holes, black holes, and wormholes. This is a population-based algorithm, where the population is searched in two phases: exploration in which the search space is visited very well and exploitation performs the search in neighborhood of a solution extensively. White hole and black hole, in the theory of multiverses, are used as exploration agents, and wormhole is acting as exploitation agent in this algorithm. In Figure 1, black, yellow, green, and pink are the objects which are moving through the white/black hole, and white points represent objects which move through the wormhole. There are 5 assumptions that are applied in MVO to the population/universe [31]:(a)The probability of having the white holes is directly proportional to the inflation rate.(b)Individuals (universes) having higher inflation rate send objects through white holes with higher probability.(c)Individuals (universes) having higher inflation rate send objects through white holes with higher probability.(d)Individuals (universes) having lower rate of inflation send objects through black holes with higher probability.(e)The objects present in all universes are moved randomly towards the current best individual (universe) through wormholes. This random process occurs without taking into account the inflation rate. A sketch of this algorithm is shown in Figure 1.

The initialization phase of MVO can be written as follows:where represents the matrix of candidate solutions, is the dimensions of the problem, and is the number of universes (population of solutions) which can be represented as follows:where is the parameter of universe, denotes the universe, NI () is the scaled inflation rate of the universe, , and is the parameter of universe selected by a roulette wheel selection mechanism. The objects are changed among the universes without perturbation. To maintain the balance between exploration and exploitation during the searching process in MVO, each universe is randomly treated to have wormholes to transport its objects through space. In order to ensure the exploitation around the current best solutions, particular wormhole tunnels are always established between a universe and the best universe formed so far, which can be represented mathematically as follows:where is the parameter of best universe formed so far, (travelling distance rate) and (wormhole existence probability) are coefficients, and are the lower bound and upper bound of variable, respectively, is the parameter of universe and r2, r3, and r4 are random numbers between 0 and 1. and are treated as adaptive formula as follows:where is the minimum (in this paper, it is set to 0.2), is the maximum (in this paper, it is set to 1), is the current iteration, is the maximum iterations, and is the exploitation accuracy over the iterations (in this paper, it is set to 6). Higher the value of , the sooner and more accurate exploitation/local search is performed [31]. It is worth to highlight that the MVO algorithm depends on number of iterations, number of universes, roulette wheel mechanism, and universe sorting mechanism. The quick sort algorithm is used to sort universe at each iteration, and roulette wheel selection is run for each variable in every universe in all iterations. Detailed description of MVO can be seen in [31].

3. Our Proposed Initialization Strategy

According to the trend in stochastic techniques, random initialization in each generation is the frequently used method to initialize algorithms without knowing any priory information about the problem in hand. The idea of random initialization is often pushing the algorithm to explore the given domain very well, while the exploitation around best solutions is not carried out by this type of initialization. We have used an initialization strategy in which IMVO is balanced in terms of exploration and exploitation. We have designed a two-stage initialization strategy: in the first stage, the algorithm is initialized randomly for a certain number of scattered solutions in the whole search domain. In stage two, after the algorithm completes a certain number of iterations and collects best points in the search domain, we initialize the IMVO with population of these best solutions to concentrate the algorithm on exploitation around the best solutions in the space. Hence, in current studies, we have assigned the function evaluations, fixed to solve a problem, in two parts. The first part initializes with random populations and collects the best solutions of each iteration by using the following equation:

In the second phase, the algorithm is directed to search in the neighborhood of the best solutions obtained in the first part and exploits the search space around the population of these best solutions.

4. Modeling of Electricity Production Cost as a Single Objective Linear Programming Problem

In this paper, we have presented an optimal mixture of energy production system as a linear programming (LP) model consisting of different decision variables, representing different electricity production sources of Pakistan, see Nomenclature [21]. Each variable is bounded and generates electricity through different fuel types like hydro, thermal, nuclear, wind, and solar energies. The suggested LP model includes different production sources as and represent different costs per for each source. The objective function represents the total cost () incurred to meet the demand of energy in different situations. Mathematically,where represents the cost per and are the different sources given in nomenclature. The objective of this study is to produce a data set of solutions so that we can select the required optimal solution fulfilling the given demand with lowest cost per as compared to other solutions in the data set. For more details, about the data set generated, the reader may refer to [35]. The problem is to search for best solutions among several candidate solutions and to meet the required power demand and satisfy constraints on each variable. These constraints are given in Equation (7) and Table 1 as follows.

Total demand constraint:

5. Experimental Settings, Results, and Discussion

We have implemented IMVO which is coded in Matlab R2016 [31]. The main features of our computer included 8 GB RAM with Intel® Core™ i7-7500U CPU @ 2.7 GHz 2.9 GHz with a 64-bit operating system Windows 10. For a fair comparison, the population size was taken as 100 and the number of iterations was 100. We have repeated our simulations 10 times. To elaborate the efficiency of IMVO, we have considered three special cases. The experimental outcome is compared with state-of-the-art algorithms: firefly algorithm, bat algorithm, and grasshopper optimization algorithm, as in Tables 2, 3, and 4. Also, we have picked 5 solutions from overall 100000 solutions and presented these solutions as different case studies with their respective costs per MWh (Table 5), which are based on the power demand in summers as provided in [20, 3638]. The power demands and decision variables, which represent the 95 production sources in Pakistan and cost incurred along with the total power produced, are given in Table 5. The individual variations of 20 best variables required to produce a particular solution in all case studies are furnished in Table 6, and a weighted comparison of the 95 energy sources is given in Figure 2. Looking at these graphs, the reader can deduce the importance of each production source of electricity in national grid of the country. One can deduce which solution to choose if a particular source is down due to some fault or which sources are very impotent and acting as a backbone in generating the required electricity. We have depicted five case studies graphically as shown in Figure 3. These graphs are collective representation of 95 production sources acting to produce a certain demand. It is interesting to note that the production of electricity from available sources is manageable, if proper techniques, like IMVO algorithm, are implemented to solve the problem. It is worth to note that we have pointed out the top 20 production sources of electricity in Pakistan according to the five case studies listed as in Table 5. These sources play vital part in the energy production as shown in Table 6. It is obvious that certain sources like are producing higher demand which are overburdened and could be relaxed by installing or expanding alternate sources. In Table 7, prices depending on different fuel types are listed for a megawatt of electricity. The convergence plots for these three case studies are depicted in Figures 46. It is interesting to note that, in all cases with demands of 15000, 19000, and 23000 megawatts, IMVO produced better solutions with lowest costs in millions (Table 24).

6. Conclusions

In this paper, an improved multiverse optimization algorithm is implemented for solving a linear programming model for proper management of electricity production sources in Pakistan. To address the proper management of power production sources, several bound constraints, such as constraint for total demand, and bounds on all variables are considered. The objective function is linear with 95 decision variables and different coefficients. The IMVO algorithm is easy to implement with few parameters. We have found solutions for five different case studies and analyzed all solutions graphically. A link to the data set of 100000 different solutions is provided in this paper. In our results, we have obtained better solutions to meet the power demand in summer. It is recommended that current resources can produce the required power demand but some of the resources are required to be expanded in capacity along with installation of new sources. In future, we intend to include other factors, like power loss, in our model.


TD:Tarbela Dam ()
MD:Mangla Dam ()
GBHPP:Ghazi Barota Hydropower Project ()
WD:Warsak Dam ()
CB:Chashma Barrage ()
DKD:Duber Khwar Dam ()
AKPH:Allai Khwar Hydropower Plant ()
KKPL:Khan Khwar Hydropower Plant ()
JHP:Jagran Hydropower Plant ()
JabHP:Jabban Hydropower Plant ()
RB:Rasal Barrage ()
DHP:Dargai Power Plant ()
GZD:Gomal Zam Dam ()
NHP:Nandipur Hydropower Plant ()
SHP:Shadiwal Hydropower Plant ()
KHP:Khuram Hydropower Plant ()
RKHP:Renela Khuard Hydropower Plant ()
CHPP:Chatral Hydropower Plant ()
GTPP:Guduu Thermal Power Plant ()
MTPP:Muzaffargarh Thermal Power Plant ()
JTPP:Jamshoro Thermal Power Plant ()
FGTPP:Faisalabad Gas Turbine Power Plant ()
MGTPP:Multan Gas Turbine Power Plant ()
KGTPP:Kotri Gas Turbine Power Plant ()
LTPP:Larkana Thermal Power Plant ()
FSPP:Faisalabad Steam Power Plant ()
SGTPP:Shahdra Gas Turbine Power Plant ()
PGTPP:Panjgur Gas Turbine Power Plant ()
QTPP:Quetta Thermal Power Plant ()
PTPP:Pasni Thermal Power Plant ()
KPC:Korangi Power Complex ()
KGTPS:Korang Gas Turbine Power Plant Station ()
GTPS:Gas Turbine Power Station ()
KPS:Thermal Power Plant ()
CCPP1:Combined Cycle Power Plant 1 ()
CCPP2:Combined Cycle Power Plant 2 ()
AESLL:AES Lalpir Limited ()
AESPG:AES Pak Gen ()
AEL:Altern Energy Limited ()
AP:Altes Power ()
AGL:Attock Gen Limited ()
CMECP:CMEC Power Pvt. ()
DHACL:DHA Cogen Limited ()
EPC:Eastern Power Company ()
EPQL:Engro Powergen Qaidpur Limited ()
FKPP:Fauji Kabirwala Power Plant ()
FTSM:First-Star Modaraba ()
FPCDL:Foundation power company Daharki Limited ()
GHLLP:Grang Holding Limited Power Plant ()
GEPtvl:Green Electric Pvt. Limited ()
GEL:Gujranwala Energy Limited ()
DAEL:Gul Ahmad Energy Limited ()
HCPC:Habibullah Coastal Power Company ()
HPGC:Halmore Power Generation Company ()
HUBCONPP:HUBCO Narowal Power Plant ()
HUBCOHPP:HUBCO Hub Power Plant ()
IPtvL:Intergen PTV Limited ()
JPG:Japan Power Plant ()
KEL:Kohinor Energy Limited ()
KAPCL:Kot Addu Power Company Limited ()
LPL:Liberty Power Limited ()
LPTL:Liberty Power Tech. Limited ()
LEPCL:Lucky Electric Power Company Limited ()
NCP:Nishat Chunian Power ()
NPL:Nishat Power Limited ()
OPCPrvL:Orient Power Company (Pvt.) Limited ()
REPGC:Radian Energy Power Generation Company ()
RPK:Rousch Power, Khanewal ()
SPC:Saba Power Company ()
SPPQ:Saif Power Plant, Qadirabad ()
SECL:Sapphire Electric Company Limited ()
SSEL:Siddiqsons Energy Limited ()
SNPCL:Sindh Nooriabad Power Company (Pvt.) Limited ()
SE:Sitara Energy ()
SEPC:Southern Electric Power Company Limited ()
SPGL:Star Power Generation Limited ()
TEL:Tapal Energy Limited ()
UchPL:Uch (Uch-I and Uch-II) Power Limited ()
NPP:Nandipur Power Project ()
QeaSP:Quaid-e-Azam Solar Park ()
YEL:Yunus Energy Limited ()
MWPCL:Metro Wind Power Co Limited ()
TGL:Tenaga Generai Limited ()
GAWPL:Gul Ahmed Wind Power Limited ()
MWEL:Master Wind Energy Limited ()
FFCEL:FFC Energy Limited, Jhimpir ()
ZEPJ:Zorlu Enerji Pakistan, Jhimpir ()
TWEL:Tapal Wind Energy Limited ()
HCDPL:HydroChina Dawood Power Limited ()
FEW-1L:Foundation Wind Energy-I Limited ()
FEW-2PT:Foundation Wind Energy-II Private Limited ().

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Authors’ Contributions

All authors contributed equally to the writing of this paper. And all authors read and approved the final manuscript.


This paper was supported by the Department of Mathematics and the Office of Research, Innovation, and Commercialization of Abdul Wali Khan University Mardan, Pakistan.