Introduction

Fossil fuel-fired power plants continue to be the primary method of generating electric power. The need to investigate alternative energy sources has increased due to the rapid rise in global electricity usage, the continuous depletion of fossil fuel reserves, and the growing environmental impact caused by the burning of fossil fuels in power plants1,2. Society's attention has been directed towards sustainable energy solutions due to the urgent need to reduce the negative effects of electricity generation on climate change3. Solar and wind power have become noticeable alternatives in this situation, acknowledged for their economic feasibility and ability to meet energy needs without causing harmful emissions4,5. However, the incorporation of these environmentally aware energy sources, such as wind and solar technologies, has brought about a level of intricacy and uncertainty in the energy sector. The emerging transition to renewable energy requires a detailed comprehension of the challenges associated to the inherent irregularity and fluctuation of solar and wind resources6. This requires a thorough examination of the dynamic properties that arise from integrating these renewable sources into the power grid. As the discussion about sustainable energy progresses, it is crucial to understand the complexities of utilizing solar and wind power to achieve their best possible integration into the overall energy system7. These insights are essential for progressing the discussion on sustainable energy usage and developing effective strategies to align the shift towards green energy with the needs of a reliable and robust power grid8. The load variation is unaffected by the unpredictability of solar irradiation and wind speed. These resources' unpredictability and sporadic nature present serious obstacles to solving the generation scheduling issue. The inherent variability and irregular characteristics of renewable energy sources, such as wind and solar power, present a potential risk to the stability and dependability of the power grid. This oscillating behavior, commonly known as "blinking," can have detrimental effects on the grid as a whole9,10. In order to address these challenges and improve the ability of the power grid to withstand disruptions, the incorporation of pumped hydroelectric energy storage is seen as a feasible solution11. Pumped-storage hydropower (PSH) units are widely recognized globally for their ability to effectively manage fluctuations in generation and supply. The growing popularity of PSH units arises from their inherent capacity to efficiently store electrical energy. Pumped-storage hydroelectric (PSH) units play a crucial role in the electric power systems by storing excess electrical energy, which is usually available and cost-effective during low-demand periods, as hydraulic potential energy12. This complex procedure entails the movement of water from the lower reservoir of the PSH unit to its upper reservoir. During times of high demand, the stored hydraulic potential energy is used to meet the increased load requirements, thereby assisting to maintain stability in the power grid. PSH units, operating on a daily or weekly basis, provide an efficient solution to mitigate the effects of renewable energy intermittency on the power grid13. Implementing Pumped Storage Hydro (PSH) units results in a gradual decrease in the overall fuel expenditure in a power system. The cost-effectiveness of this approach is due to the strategic placement of PSH units, which helps to stabilize fluctuations in energy supply and demand and optimize the operation of the power system14. Overall, the integration of pumped hydroelectric energy storage, demonstrated by PSH units, is an effective approach to mitigate the intermittent nature of renewable energy sources. By utilizing the storage capabilities of PSH, the power grid can attain heightened stability, decreased operational expenses, and enhanced flexibility to accommodate the ever-changing landscape of renewable energy generation15,16.

A modest sovereign system's optimal generation scheduling using renewable energy sources has been covered in17. Although clean and pollution-free, renewable energy sources have a limited ability to provide electricity. The optimal approach to address the economic dispatch quandary lies in dynamic economic dispatch (DED). This approach efficiently distributes the time-varying load demand across all active generating units, while taking into account the limitations presented by thermal generator ramp rates18. In the realm of Dynamic Economic Dispatch (DED), decisions made at one time significantly influence subsequent decisions. Addressing this, a novel Enhanced Non-Dominated Sorting Crisscross Optimization (ENSCSO) algorithm was introduced to solve the multi-objective Dynamic Economic Emission Dispatch problem19. This algorithm was tested via simulations on a ten-unit generation system that integrates wind power and a time-of-use demand response program. Ameliorated dragonfly algorithm (ADFA) was applied to solve static economic load dispatch and dynamic economic load dispatch problem in20. Static economic dispatch was carried out on three different test systems and dynamic economic dispatch was implemented on two different test systems. In21, a Levy Interior Search Algorithm was crafted with a focus on resolving the multi-objective economic load dispatch issue, integrating the incorporation of wind power. The objective functions considered were operation cost and system risk. A simulation was conducted using a modified IEEE 30-bus system, incorporating the integration of wind power. A distributed structure and stochastic linear programming game were presented, allowing for the scheduling of appliances and storage units as well for the energy payments in22. A distributed primal–dual continuous time consensus algorithm was implemented for solving dynamic economic dispatch problem23. Simulation was carried out on three different test systems. An improved version of Circle Search Algorithm was introduced in ref.24 to resolve the economic emission dispatch problem by incorporating demand response integration. Improved circle search algorithm was investigated on IEEE 6-bus and IEEE 30-bus system to implemented the multi-objective economic emission dispatch25. In26, multi-objective particle swarm optimization was proposed to solve the dynamic economic emission dispatch problem. Within the Demand-Side Management (DSM) process, a strategy utilizing day-ahead load shifting techniques was implemented to manage residential loads. The primary objective involved minimizing the utility's energy bill. The application of the Interior Search Algorithm was utilized to address the economic load dispatch problem within a microgrid setting, as referenced in27. Multi-objective dynamic optimal power flow problem was implemented using harmony search algorithm. In27, the day-ahead load shifting DSM technique was enacted using a day-ahead pricing strategy combined with an energy consumption game. Additionally, in28, the successful implementation of the normal boundary intersection method effectively addressed the centralized multi-objective dynamic economic dispatch incorporating demand side management for individual residential loads and electric vehicles. Generation costs, emissions, and energy loss are considered as objective functions. A suite of innovative optimization algorithms was developed to tackle various complex challenges within energy management systems. In29, the Improved Mayfly Optimization Algorithm was devised to solve the combined economic emission dispatch problem within a microgrid setting. Ref.30 introduced the Chaotic Fast Convergence Evolutionary Programming (CFCEP) aimed at resolving the combined heat and power dynamic economic dispatch problem. This solution incorporated demand side management, renewable energy sources, and pumped hydro energy storage. The Social Group Entropy Optimization (SGEO) technique, highlighted in reference31, was proposed to address the fuel-constrained dynamic economic dispatch problem. This strategy combined demand-side management, renewable energy sources, and a pumped hydro storage plant. It implemented a Multi-Objective Dynamic Economic Emission Dispatch by incorporating game theory-based demand response techniques32. Lastly, in ref.33, a Multi-Objective Dynamic Economic Emission Dispatch approach was applied within a microgrid context. This implementation incorporated demand response strategies along with a zero-balance approach.

As outlined in the International Energy Agency's strategic plan, DSM stands as the major choice for energy policy decisions. DSM programs offer various advantages, such as cost reduction and heightened security within power systems. Here's an overview of the ongoing research contributions in this domain:

  1. 1.

    In our research paper, we introduce the Enhanced Cheetah Optimizer Algorithm (ECOA) to address dynamic economic dispatch while integrating renewable energy sources and demand side management. We've integrated chaotic sine map and levy flight mechanism and into this algorithm to improve solution quality and convergence speed. This learning method involves simultaneously considering an estimate and its opposite counterpart, aiming to refine the current candidate solution more effectively.

  2. 2.

    The inherent variability of wind and solar power generators is depicted through the utilization of the most reliable probability density functions (PDFs).

  3. 3.

    The alteration in the generation costs of wind and solar power in relation to the respective scheduled power adjustments is thoroughly investigated.

  4. 4.

    We subjected our proposed algorithm to a comprehensive evaluation to assess its effectiveness in addressing dynamic economic dispatch challenges associated with pumped-storage hydroelectric units and demand-side management. The ECOA algorithm we introduced plays a vital role in determining optimal times for both pumping water to the upper reservoir and releasing it for power generation, taking into consideration factors such as electricity prices, demand patterns, and the availability of renewable energy.

  5. 5.

    The algorithm we proposed was thoroughly examined for its efficacy in resolving dynamic economic dispatch problems involving unconventional energy sources and demand side management. We compared the optimization results of our proposed algorithm with those obtained using COA and GWO for comprehensive analysis.

Mathematical formulation of dynamic economic dispatch

The primary aim of integrating renewable energy sources into the Dynamic Economic Dispatch (DED) system is to achieve a dual objective of minimizing two factors simultaneously34. The primary objective is to minimize the overall expenses linked to thermal power plants by enhancing their operating efficiency. Furthermore, the integration aims to reduce costs associated with the functioning of wind-power producing units and solar Photovoltaic (PV) facilities. This extensive framework of DED expands its scope to incorporate the integration of pumped hydroelectric energy storage, acknowledging its crucial role in mitigating the intermittent nature of renewable energy sources35,36. The study attempts to achieve an efficient and cost-effective balance between traditional and renewable energy sources within the dynamic economic dispatch framework using this integrated method37.

The formulation of the DED problem with DSM encompasses defining the resultant objective function along with its associated constraints. The fuel cost function for the ith thermal generator at time \(t\), accounting for the valve-point effect38,39, is expressed as:

$$F\left( {P_{G} } \right) = \mathop \sum \limits_{i = 1}^{{N_{TH} }} (a_{i} + b_{i} P_{Gi} + c_{i} P_{Gi}^{2} ) + |e_{i} {\text{sin}}(f_{i} *\left( {P_{Gimin} - P_{i} } \right)|$$
(1)

where \(a_{i}\), \(b_{i}\) and \(c_{i}\) are fuel cost coefficients of \(i^{th}\) generator, k is the total number of generating units, \(P_{Gi}\) is the output power of the \(i^{th}\) generator in megawatts. Here \(e_{i}\) and \(f_{i}\) represents the generating cost coefficients of the \(i^{th}\) unit are utilized to model valve point loading effect.

Modelling the costs of renewable energy sources

Assessment of direct costs for wind and solar photovoltaic power

The functioning of energy generation from RESs doesn't require any fuel. Hence, in cases where Independent System Operators (ISO) own Renewable Energy Sources (RESs), only maintenance costs are incurred without any associated cost function40. Yet, if private organizations manage RESs, the ISO compensates them as per the mutually agreed-upon contract for the scheduled electricity generation40.

The assessment of direct costs for wind turbines and solar photovoltaic (PV) power involves a detailed examination of the expenses associated with the design, construction, installation, operation, and maintenance of these renewable energy systems41,42. Direct costs are those directly attributable to the development and operation of the specific technology.

The literature offers the direct cost function for the \(ith\) wind farm concerning the planned power40.

$$C_{W} \left( {P_{W} } \right) = K_{W} P_{W}$$
(2)

Here, \(P_{W}\) represents the generated power and \(K_{W}\) represents the direct cost coefficient related to the wind turbine. In connection with is, the direct cost involved in solar PV with scheduled power \(P_{PV}\) and cost coefficient, \(K_{PV}\) is represented by the following equation

$$C_{PV} \left( {P_{PV} } \right) = K_{PV} P_{PV}$$
(3)

In this context, \(P_{PV}\) denotes the generated power, while \(K_{PV}\) represents the direct cost coefficient associated with solar photovoltaic generation.

Assessment of reserve cost and penalty cost associated with wind power

As wind energy is inherently unpredictable, the power generated by wind turbines fluctuates over time, potentially surpassing or falling short of the scheduled power43. Therefore, the ISO needs to have backup generating capacity to meet demand. The assessment of reserve cost and penalty cost associated with wind power involves evaluating the expenses and penalties incurred due to the intermittent and variable nature of wind energy. Reserve costs and penalty costs are critical aspects in the economic evaluation and operational planning of power systems that include wind power44.

The reserve cost for the wind unit is presented based on the literature40.

$$\begin{aligned} C_{{R_{W,i} }} \left( {P_{W sh,i} - P_{W ac,i} } \right) & = k_{rw,i}^{\prime } \left( {P_{W sh,i} - P_{W ac,i} } \right) \\ & = k_{rw,i}^{\prime } \mathop \smallint \limits_{0}^{{P_{W sh,i} }} \left( {P_{W sh,i} - p_{w,i} } \right)f_{w} \left( {p_{w,i} } \right)dp_{w,i} \\ \end{aligned}$$
(4)

When wind generators produce more output power than is scheduled, the ISOs must pay the fine by reducing the power of thermal generators when they do not consume the extra power.

$$\begin{aligned} C_{{P_{W,i} }} \left( {P_{W\; ac,i} - P_{W\; sh,i} } \right) & = k_{rw,i}^{\prime } \left( {P_{W\; ac,i} - P_{W\; sh,i} } \right) \\ & = k_{pw,i}^{\prime } \mathop \smallint \limits_{{P_{W\; sh,i} }}^{{P_{W r,i} }} \left( {p_{w,i} - P_{W\; sh,i} } \right)f_{w} \left( {p_{w,i} } \right)dp_{w,i} \\ \end{aligned}$$
(5)

Here \(P_{W \;sh,i}\) represents the scheduled wind power and \(P_{W\; ac,i}\) represents the actual power generated by the wind turbine. The rated power is represented by \(P_{W\; r,i}\) while the Probability Density Function (PDF) of the wind power is represented as \(f_{w} \left( {p_{w,i} } \right)\). Accordingly, it is possible to compute the reserve and penalty costs for solar-only and solar-and-hydro combined generators45. The required input data for modeling the cost of renewable energy sources is obtained from existing literature30.

Assessment of the penalty and reserve cost associated with PV power

Assessing the cost of generating wind power aligns closely with formulating the stochastic generation cost of solar PV electricity40. Moreover, the lognormal Probability Density Function (PDF) proves useful in representing solar radiation40. Additionally, reserve and penalty cost models for solar PV-powered plants are devised based on the methodology outlined in Reference46. Section 2.3 does the computation for the solar photovoltaic unit's generated power output. The assessment of reserve cost and penalty cost associated with photovoltaic (PV) power involves evaluating the expenses and penalties incurred due to the intermittent and variable nature of solar energy. Reserve costs and penalty costs are critical aspects in the economic evaluation and operational planning of power systems that include solar PV.

The reserve cost for overestimating solar PV power is characterized as40:

$$\begin{aligned} C_{{R_{PV,i} }} \left( {P_{PV\; sh,i} - P_{S ac,i} } \right) & = k_{rpv,i}^{\prime } \left( {P_{PV\; sh,i} - P_{PV\; ac,i} } \right) \\ & = k_{rpv,i}^{\prime } f_{pv} \left( {P_{PV\; ac,i} < P_{PV\; sh,i} } \right)\left[ { P_{PV\; sh,i} - E\left( {P_{PV\; ac,i} < P_{PV\; sh,i} } \right) } \right] \\ \end{aligned}$$
(6)

where \(P_{PV\; ac,i}\) represents the actual power generated by the solar PV plant, and \(k_{rpv,i}^{\prime }\) denotes the reserve cost coefficient related to the solar PV plant. The expectation of solar PV power below P_ is represented by \(E\left( {P_{PV \;ac,i} < P_{PV\; sh,i} } \right)\), and the likelihood of a solar power shortage from the scheduled solar PV power is given by \(f_{pv} \left( {P_{PV\; ac,i} < P_{PV\; sh,i} } \right).\) The cost of the penalty for underestimating solar PV power is characterised as40:

$$\begin{aligned} C_{{P_{PV,i} }} \left( {P_{PV\; ac,i} - P_{PV\; sh,i} } \right) & = k_{ppv,i}^{\prime } \left( {P_{PV \;ac,i} - P_{PV\; sh,i} } \right) \\ & = k_{ppv,i}^{\prime } f_{pv} \left( {P_{PV\; ac,i} > P_{PV\; sh,i} } \right)\left[ { E\left( {P_{PV\; ac,i} > P_{PV\; sh,i} } \right) - P_{PV \;sh,i} } \right] \\ \end{aligned}$$
(7)

where \(k_{ppv,i}^{\prime }\) represents the penalty cost coefficient for the solar PV plant and \(f_{pv} \left( {P_{PV\; ac,i} > P_{PV \;sh,i} } \right)\) characterizes the likelihood that solar power will be above the scheduled power (\(P_{PV \;sh,i}\)), and \(E\left( {P_{PV\; ac,i} > P_{PV\; sh,i} } \right)\) indicates the expectation that solar PV power will be above \(P_{PV \;sh,i}\).

Formulation of overall generation cost with the integration of renewable energy sources

The overall operation cost is a critical metric that reflects the economic efficiency of the power system operation. The overall operation cost considers the intermittent nature of renewable energy sources, accounting for periods of high and low generation, and the associated economic implications.

The overall operation cost within the DED problem is structured as follows30,40:

Minimize

$$\begin{aligned} F_{C} & = F\left( {P_{G} } \right) + \mathop \sum \limits_{i = 1}^{{N_{WG} }} \left[ { C_{W} \left( {P_{W} } \right) + C_{{R_{W,i} }} \left( {P_{W\; sh,i} - P_{W \;ac,i} } \right) + C_{{P_{W,i} }} \left( {P_{W \;ac,i} - P_{W \;sh,i} } \right)} \right] \\ & \quad + \mathop \sum \limits_{i = 1}^{{N_{PV} }} \left[ {\left( {C_{PV} (P_{PV} } \right) + C_{{R_{PV,i} }} \left( {P_{PV\; sh,i} - P_{PV \;ac,i} } \right) + C_{{P_{PV,i} }} \left( {P_{PV\; ac,i} - P_{PV\; sh,i} } \right)} \right] \\ \end{aligned}$$
(8)

Equality and inequality constraints

Equality constraints

Generator power output constraint

The total power generation, when combined with demand-side management, can be expressed through the following Eq30:

$$\mathop \sum \limits_{i = 1}^{{N_{TH} }} (P_{Git} ) + \mathop \sum \limits_{i = 1}^{{N_{W} }} (P_{Wit} ) + \mathop \sum \limits_{i = 1}^{{N_{PV} }} (P_{PVit} ) + \mathop \sum \limits_{i = 1}^{{N_{pump} }} (P_{GHit} ) = \left( {1 - DR_{t} } \right) \times L_{Base,t} + L_{{S_{t} }} + P_{loss}$$
(9)
$$\mathop \sum \limits_{i = 1}^{{N_{TH} }} (P_{Git} ) + \mathop \sum \limits_{i = 1}^{{N_{W} }} (P_{Wit} ) + \mathop \sum \limits_{i = 1}^{{N_{PV} }} (P_{PVit} ) - \mathop \sum \limits_{i = 1}^{{N_{pump} }} (P_{PHit} ) = \left( {1 - DR_{t} } \right) \times L_{Base,t} + L_{{S_{t} }} + P_{loss}$$
(10)

where \(N_{TH}\), \(N_{W}\), \(N_{PV}\) and \(N_{pump}\) denotes the quantity of thermal power units, wind power units, solar photovoltaic units, and pumped storage units, respectively. The power generated by the ith thermal, wind, solar photovoltaic, and pumped storage units is represented as \(P_{Gi}\), \(P_{Wi}\), \(P_{PVi}\) and \(P_{GHi}\) respectively.

In order to achieve optimal economic load dispatch, one must include transmission line losses. The transmission line losses are calculated using Newton–Raphson methods and B-coefficient methods. In order to calculate the active power loss \(P_{loss} .\) Newton–Raphson method is used in conjunction with the power flow solution. The subsequent equation defines the actual power loss while adhering to equality prerequisites40.

$$P_{Gj} - P_{Dj} - V_{j} \mathop \sum \limits_{j = 1}^{NB} V_{k} \left[ {G_{jk} \cos \left( {\delta_{j} - \delta_{k} } \right) + B_{jk} sin\left( {\delta_{j} - \delta_{k} } \right)} \right] = 0$$
(11)
$$Q_{Gj} - Q_{Dj} - V_{j} \mathop \sum \limits_{j = 1}^{NB} V_{k} \left[ {G_{jk} \sin \left( {\delta_{j} - \delta_{k} } \right) + B_{jk} cos\left( {\delta_{j} - \delta_{k} } \right)} \right] = 0$$
(12)

With \(j = 1, 2, \ldots NB\); in this case, \(NB\) represents the total number of buses. \(V_{j}\) and \(V_{k}\) represents the \(jth\) bus and \(kth\) bus voltage respectively. \(Q_{gj}\) denotes the \(jth\) bus reactive power output and \(\delta_{j}\) and \(\delta_{k}\) characterizes the voltage angle at bus \(j\) and bus \(k\) respectively.\(B_{jk}\) and \(G_{jk}\) represents the transfer susceptance and conductance between buses j and \(k\) respectively. \(P_{Dj}\) and \(Q_{Dj}\) represents the \(jth\) bus active and reactive power load respectively. In order to determine the equality constraints, the Newton–Raphson load flow technique solution is used. Bus voltage magnitudes and angles can be determined using the power flow solution.

Inequality constraints

Limits on the lowest and highest generation capacities

Each generator's active power generation output needs to stay within specific minimum and maximum limits40. Power generation constraints refer to the limitations and restrictions imposed on the operation of power generation units over time. These constraints are crucial for ensuring the secure and reliable operation of the power system.

$$P_{Gimin} \le P_{Gi} \le P_{Gimax} \forall i \in N_{TH}$$
(13)
$$P_{w}^{min} \le P_{w} \le P_{w}^{max}$$
(14)
$$P_{PV}^{min} \le P_{PV} \le P_{PV}^{max}$$
(15)
Pumped-storage constraints

The integration of pumped-storage hydro units adds a dynamic and flexible component to the system, enabling better balancing of supply and demand.

The net water usage of the pumped-storage hydropower (PSH) unit should balance out to zero as the final and initial water volumes in the upper reservoir are considered equal within this scenario30.

$$V_{{res,j\left( {t + 1} \right)}} = V_{res,jt} + Q_{phjt} \left( {P_{phjt} } \right), j \in N_{pump} , t \in T_{pump}$$
(16)
$$V_{{res,j\left( {t + 1} \right)}} = V_{res,jt} - Q_{ghjt} \left( {P_{ghjt} } \right), j \in N_{pump} , t \in T_{gen}$$
(17)
$$P_{ghj}^{min} \le P_{ghj} \le P_{ghj}^{max} , j \in N_{pump} , t \in T_{gen}$$
(18)
$$P_{phj}^{min} \le P_{phj} \le P_{phj}^{max} , j \in N_{pump} , t \in T_{pump}$$
(19)
$$V_{res,j}^{min} \le V_{res,jt} \le V_{res,j}^{max} , j \in N_{pump} , t \in T$$
(20)

Given the equality between the initial and final water volumes of the upper reservoir in the pumped-storage hydroelectric (PSH) unit for this scenario, the total net water used by the PSH unit should equate to zero30.

$$V_{res,j0} = V_{res,jT} = V_{res,j}^{start} = V_{res,j}^{end}$$
(21)
$$\begin{aligned} Q_{net,spent,j} & = Q_{spent,TOT,j} - Q_{pump,TOT,j} \\ & = \mathop \sum \limits_{{t \in T_{gen} }} Q_{ghjt} \left( {P_{ghjt} } \right) - \mathop \sum \limits_{{t \in T_{pump} }} Q_{phjt} \left( {P_{phjt} } \right) = 0 \\ \end{aligned}$$
(22)
Ramp rate limits of thermal generator

The ramp rate limits of thermal generators are crucial parameters in power system operation and control. The ramp rate refers to the maximum rate at which the power output of a generator can change over a specified time interval. Rapid and large changes in power output from generators can lead to instability in the power grid. By imposing ramp rate limits, the system operators ensure that the changes in power output are gradual, helping to maintain grid stability.

$$P_{Git} - P_{{Gi\left( {t - 1} \right)}} \le UR_{i} , i \in N_{t} , t \in T$$
(23)
$$P_{{Gi\left( {t - 1} \right)}} - P_{Git} \le DR_{i} , i \in N_{t} , t \in T$$
(24)

Wind, solar and hydro uncertainty models

To represent the unpredictable output power from Renewable Energy Sources (RESs), a range of Probability Density Functions (PDFs) are utilized.

The wind speed determines how much power the wind turbines can produce. According to past research investigations40,46, the likelihood of wind speed follows Weibull PDF.

The Weibull distribution is commonly used in the field of wind energy because it is well-suited for modeling the variability of wind speeds at a particular location.

$$f_{wv} \left( v \right) = \left( {\frac{\alpha }{\lambda }} \right)\left( {\frac{v}{\lambda }} \right)^{{\left( {\alpha - 1} \right)}} exp[^{{ - \left( {\frac{v}{\lambda }} \right)}} ]^{\alpha } \;for\; 0 < v < \infty$$
(25)

where \(\alpha\) represents the scale of the Weibull PDF and stands for the shape parameter of the Weibull PDF. These variables' values were collected from30. Weibull PDF's median is provided by:

$$M_{w} = \lambda *\Gamma \left( {1 + \alpha^{ - 1} } \right)$$
(26)

The gamma (\(\Gamma )\) function is crucial in the context of the Weibull probability density function (PDF) for wind distribution because it is used to normalize the Weibull distribution and ensure that it integrates to 1 over its entire range.

\(\Gamma\) function can be represented as:

$$\Gamma \left( {x^{\prime } } \right) = \mathop \smallint \limits_{0}^{\infty } e^{ - t} t^{{x^{\prime } - 1}} dt$$
(27)

As shown in Fig. 1, the frequency distribution is derived from Weibull fitting using wind speed results obtained through simulating 8000 Monte Carlo scenarios. The values for the scale and shape parameters are sourced from30. Consistent with the literature30, the PDF parameter values have been selected. Achieving a cumulative rated output of 175 MW involves the collective output from 35 wind generators, each possessing a capacity of 5 MW. The subsequent equation delineates the power generated by the wind turbines, contingent upon the wind speed.

$$P_{WG} = \left\{ {\begin{array}{*{20}l} 0 \hfill & {for\; v \le v_{in} } \hfill \\ {P_{{W_{r} }} \left( {\frac{{v - v_{in} }}{{v - v_{out} }}} \right)} \hfill & {for \;v_{in} \le v \le v_{r} } \hfill \\ {P_{{W_{r} }} } \hfill & {for \;v_{r} \le v \le v_{out} } \hfill \\ \end{array} } \right.$$
(28)

where \(P_{{W_{r} }}\) denotes the rated power of a single turbine. \(v_{in}\) signifies the cut-in speed, \(v_{out}\) denotes the cut-out speed whereas \(v_{r}\) is the rated speed. The study investigated different Weibull parameters that dictate the distribution of wind speeds, in line with the selections made in previous studies30. Equation (40) emphasizes the discrete nature of the wind generator's output power, notably in specific regions. Specifically, wind farm output remains at zero when wind speed falls below the cut-in speed or exceeds the cut-out speed. The wind generators operate at their rated power within the range delineated between the cut-out and cut-in regions. Previous studies30,40 detail the probability associated with these discrete zones.

$$f_{{P_{WG} }} = 1 - {\text{exp}}\left[ { - \left( {\frac{{v_{in} }}{\lambda })^{\alpha } } \right] + {\text{exp}}} \right[ - \left( {\frac{{v_{out} }}{\lambda })^{\alpha } } \right]\;for\; \left( {P_{WG} = 0} \right)$$
(29)
$$f_{{P_{WG} }} = {\text{exp}}\left[ { - \left( {\frac{{v_{r} }}{\lambda })^{\alpha } } \right] - {\text{exp}}} \right[ - \left( {\frac{{v_{out} }}{\lambda })^{\alpha } } \right]\;for\; \left( {P_{WG} = P_{WR } } \right)$$
(30)
Figure 1
figure 1

Wind speed variation in wind power generation unit.

In the continuous domain, the probability distribution for wind power is expressed as follows40,46:

$$f_{{P_{WG} }} = \frac{{\alpha \left( {v_{r} - v_{in} } \right)}}{{\lambda^{\alpha } *P_{wr} }}\left[ {v_{in} + \frac{{P_{WG } }}{{P_{wr} }}\left( {v_{r} - v_{in} } \right)} \right]^{\alpha - 1} exp - \left( {\frac{{(v_{in} + \frac{{P_{WG } }}{{P_{wr} }}\left( {v_{r} - v_{in} } \right)}}{\lambda }} \right)^{\alpha }$$
(31)

This Weibull PDF is utilized to characterize and model the probability distribution of wind speeds, which is crucial for assessing the potential power output of wind turbines.

Furthermore, the solar photovoltaic (PV) output power is solely contingent on solar irradiance (G), conforming to the parameters of the lognormal Probability Density Function (PDF)40,46. A previous study40 outlined the probability distribution of solar irradiance, specifying its mean and standard deviation. The lognormal distribution is often used in PV modeling because it provides a good fit for the skewed and positive-valued nature of solar irradiance and power output data. Many natural processes, including solar irradiance, exhibit lognormal characteristics, making the lognormal distribution a suitable choice for modeling. The lognormal distribution is well-suited for data with a positively skewed distribution, capturing the asymmetric behavior often observed in solar irradiance data. The parameters in the lognormal distribution have physical interpretations, such as the mean and standard deviation, which can provide insights into the characteristics of the solar resource.

$$f_{PV} \left( G \right) = \frac{1}{{G\sigma \sqrt {2\pi } }}\exp \left[ { - \frac{{\left( {\ln x - \mu } \right)^{2} }}{{2\sigma^{2} }}} \right] \;for\; G > 0$$
(32)

The subsequent equation represents the mean of the lognormal distribution (\(M_{Lgn}\))

$$M_{Lgn} = \exp \left( {\mu + \frac{{\sigma^{2} }}{2}} \right)$$
(33)

After running 8,000 Monte Carlo simulations, a frequency distribution for solar irradiance is derived, and Fig. 2 illustrates the lognormal fitting, demonstrating the solar PV output power.

$$P_{PV} \left( G \right) = \left\{ {\begin{array}{*{20}l} {P_{PVr} \left( {\frac{{G^{2} }}{{G_{std} R_{C} }}} \right)} \hfill & {for\; 0 \le G \le R_{C} } \hfill \\ {P_{PVr} \left( {\frac{G}{{G_{std} }}} \right)} \hfill & { for\; G \ge R_{C} } \hfill \\ \end{array} } \right.$$
(34)
Figure 2
figure 2

Distribution of solar irradiance for solar PV.

The critical value (\(R_{C} )\) introduces a threshold beyond which the model transitions to a simpler form. This threshold may represent a point where the PV system behavior changes, possibly due to system constraints, saturation effects, or other factors.

In the standard environmental conditions, standard deviation of solar irradiance is represented by \(G_{std}\) and certain irradiance is characterized by \(R_{C}\). The assumed value for \(G_{std}\) stands at 1000 W/m2, whereas for \(R_{C}\), it amounts to 150 W/m2. Regarding the PV module, the rated output power \(P_{PVr}\) is specified as 175 MW.

Demand-side management

DSM initiatives bring forth numerous benefits such as cost efficiency and improved power system security47. These programs encompass various categories, prominently featuring demand response. Among these, the time-of-use (TOU) program48 stands out—it redistributes a segment of the load demand from peak hours to off-peak periods or times of lower cost, while maintaining the overall load demand. This TOU program served as the foundational inspiration for the demand response program applied in this study. This flattens the load curve and lowers the expected operation cost. The numerical model for the TOU program is created in line with Eq. (35) and is constrained by Eqs. (36) - (39).

$$L_{t} = \left( {1 - DR_{t} } \right) \times L_{Base} + L_{st}$$
(35)
$$\mathop \sum \limits_{t = 1}^{T} L_{st} = \mathop \sum \limits_{t = 1}^{T} DR_{t} \times L_{Base,t}$$
(36)
$$L_{{Inc_{t} }} = Inc_{t} \times L_{Base,t}$$
(37)
$$DR_{t} \le DR^{max} , t \in T$$
(38)
$$Inc_{t} \le Inc^{max} , t \in T$$
(39)

Enhanced Cheetah optimizer algorithm

Akbari et al.49 introduced the COA algorithm, drawing inspiration from the hunting techniques of cheetahs. This method integrates three primary strategies: prey search, ambush tactics, and active attacks. Significantly, it implements a mechanism to navigate away from a prey location and return to a home position, effectively avoiding entrapment in local optimal points. Each cheetah's potential hunting patterns correspond to potential solutions for the problem at hand. The algorithm operates on the premise that the population's best position determines the optimal solution, akin to identifying the prey. Cheetahs adapt their hunting patterns to enhance their performance over the hunting period. By mimicking these strategies, the COA algorithm49 effectively seeks optimal solutions for intricate problems.

When a cheetah scans its surroundings, it can detect potential prey, giving it the option to either wait for the prey to approach or to initiate an immediate attack upon spotting it. The attack itself involves two distinct phases: a rapid approach followed by capture. However, several factors might prompt the cheetah to abandon the hunt, such as low energy reserves or if the prey is too agile. In such scenarios, the cheetah might retreat to its resting spot, preparing for a fresh hunting opportunity. The cheetah carefully assesses the prey's condition, the environment, and the distance involved before choosing between these strategies. The COA algorithm encapsulates this entire hunting process, relying on the strategic utilization of these tactics across multiple hunting cycles or iterations49. Essentially, the COA algorithm leverages these intelligent hunting strategies iteratively throughout the hunting process.

  1. i.

    Searching: Cheetahs engage in scanning or active search within their territories or the surrounding area to locate prey within the search space.

  2. ii.

    Sitting-and-waiting: Upon detecting prey but under unfavorable conditions, cheetahs may opt to sit and wait, allowing the prey to approach or for a better opportunity to arise.

  3. iii.

    Attacking: This strategy involves two crucial phases:

    1. a.

      Rushing: Once committed to an attack, cheetahs sprint toward the prey at maximum speed.

    2. b.

      Capturing: Leveraging speed and agility, cheetahs capture the prey by closing in swiftly.

  4. iv.

    Returning home and leaving prey: This strategy comes into play under two circumstances. Firstly, if the cheetah fails to catch its prey, it may choose to relocate or return to its territory. Secondly, when a certain time lapses without successful hunting, the cheetah may reposition itself to the last known prey location and conduct further searches in that area49. Detailed mathematical models for these hunting strategies are expounded upon in subsequent sections.

The CO algorithm has demonstrated strong capabilities in tackling expansive problems. However, as the upcoming experimental results will demonstrate, there remains an opportunity for improvement in terms of convergence speed and computational time, particularly when fine-tuning the parameters of photovoltaic models. To overcome these limitations, we present an upgraded iteration of the COA algorithm tailored explicitly to tackle these drawbacks.

Searching strategy

In the exploration phase of the COA algorithm, each cheetah adjusts its position by referencing its prior location. Cheetahs commonly follow the lead of the leader within their group. Expanding upon this notion, the search approach detailed in Eq. (16) is adapted based on the position of the group's second-best cheetah, designated as \(X_{L,j}^{t}\), influencing the modification process. This adjustment is detailed as follows50:

$$X_{i,j}^{t + 1} = X_{L,j}^{t} + \hat{r}^{t} \cdot \alpha_{i,j}^{t}$$
(40)

where the randomization parameter (\(\hat{r}^{t}\)) and the random step length (\(\alpha_{i,j}^{t}\)) undergo modifications as follows:

The value of the randomization parameter (\(\hat{r}^{t}\)) in Eq. (40) can be ascertained through the implementation of a sine map, where the initial values for \(C_{t}\) and \(a\) are specifically set at 0.36 and 2.8 as indicated in reference51.

$$C_{t + 1} = \frac{a}{4}\sin \left( {\pi C_{t} } \right), 0 < a < 4$$
(41)

where \(t\) represents the current iteration number.

The random step length (\(\alpha_{i,j}^{t}\)) can be represented as

$$\alpha_{i,j}^{t} = X_{{k^{\prime } ,j}}^{t} - X_{{i^{\prime } ,j}}^{t}$$
(42)

Here \(X_{{k^{\prime},j}}^{t}\) and \(X_{{i^{\prime},j}}^{t}\) are the positions of \(kth\) and \(ith\) cheetahs in the sorted population, respectively.

Emphasizing the alignment of every cheetah's position around the group leader plays a crucial role in the local search phase. Furthermore, the second term in Eq. (40) enhances solution diversity, actively aiding in the global search or exploitation phase. In addition, introducing substantial strides during the hunting phase via the random parameter extends solutions beyond variable ranges. Subsequently, these are substituted by fresh random solutions within the population. This dual purpose not only broadens the spectrum of solutions but also shields the algorithm from being stuck in local optimum points.

Attacking strategy

To bolster the optimization capabilities of COA algorithm, the researcher crafted the Enhanced Cheetah Optimizer (ECOA) algorithm. This new approach combines principles inspired by Levy flights, mirroring the flight patterns observed in birds. Adopting a Levy flight-based approach for system identification offers expedited convergence without relying on derivative information40. This method employs stochastic random searches based on Levy flight concepts52. Integrating the Levy flight approach bolsters local search capabilities, mitigating the risk of local entrapment for the optimal solution52.

Furthermore, the attacking strategy within the ECOA algorithm undergoes reformulation as follows40:

$$X_{i,j}^{t + 1} = X_{B,j}^{t} + Levy\left( \lambda \right) \cdot \beta_{i,j}^{t}$$
(43)
$$Levy\left( \lambda \right) = 0.01 \frac{{r_{1} \sigma }}{{\left| {r_{2} } \right|^{{\frac{1}{\beta }}} }}$$
(44)

where σ can be calculated as40:

$$\sigma = [\Gamma \left( {1 + \lambda } \right){\text{sin}}\left( {\pi \frac{\lambda }{2}} \right)/\left( {\Gamma \left( {\frac{1 + \lambda }{2}} \right)\lambda \left[ {2^{{\frac{{\left( {\lambda - 1} \right)}}{2}}} } \right]} \right)]^{1/\lambda }$$
(45)

The function \(\Gamma \left( x \right)\) represents the factorial of (x-1), while \(r_{1}\) and \(r_{2}\) denote indiscriminate numbers within the range of \(\left[ {0,1} \right].\) For \(1 < \beta \le 2\), a constant value (e.g., 1.5) for \(\beta\) is specifically applied in this research49. The symbol \(Levy\left( \lambda \right)\) signifies step length, employing the Levy distribution characterized by infinite variance and a mean of \(1 < \lambda < 3\). \(\lambda\) serves as the distribution factor, with \(\Gamma \left( . \right)\) representing the gamma distribution function.

Within the COA algorithm, the interaction factor considers the position of neighbouring cheetahs. Ordinarily, cheetahs hunt individually, adapting their positions in response to their prey's whereabouts. Therefore, in this newly suggested attack strategy, each cheetah adjusts its position relative to the prey during the attack phase, advancing toward it following this formula50:

$$\beta_{i,j}^{t} = X_{B,j}^{t} - X_{i,j}^{t}$$
(46)

This refined attack strategy significantly accelerates the COA algorithm's ability to approach near-optimal solutions swiftly. It bolsters the algorithm's local search prowess (exploitation phase), thus amplifying its convergence speed. Figure 3 showcases the schematic of the enhanced Cheetah Optimizer Algorithm as proposed.

Figure 3
figure 3

Flowchart of the proposed enhanced cheetah optimizer algorithm.

Results and discussion

Test system: I

The proposed approach has been deployed to address the dynamic economic dispatch problem, both with and without DSM. To gauge its effectiveness, optimization outcomes were compared against COA, GWO, CFCEP30, FCEP30, CCDE30, and HSPSO30. The MATLAB 9.12 software was utilized to implement the ECOA, COA, and GWO53 models on a Laptop with an AMD Athlon processor, 1 TB storage, and 3.0 GHz processing speed. The test system encompasses 10 thermal power plants, one equivalent wind turbine, a solar photovoltaic plant, and a pumped-storage hydroelectric plant. The scheduling spans 24 intervals, considering the valve-point loading effect on thermal generators. The input data, including bus data, PDF parameters, and cost coefficients, were gathered from a preceding study30. Notably, during intervals 11, 12, and 13, peak loads are identified, prompting DSM to redistribute 10% of the load from these hours to the 2nd, 3rd, and 4th intervals. It's important to note that the pumped-storage hydroelectric (PSH) plant operates in generating mode specifically when both the power generated and discharge rate are positive. Conversely, it functions in pumping mode when pumping power and pumping rate are negative30.

The Weibull PDF parameters in this case are chosen from Ref.30. The direct cost coefficients, penalty cost coefficients, and reserve cost coefficients for wind power are sourced from literature30. Notably, the direct cost of renewable power is lower than the average cost of thermal power. Additionally, the penalty incurred for underutilizing available wind power is less than the direct cost40. Examining the scheduled power range from 0 to the wind farm's rated power, Fig. 4 illustrate the variations in reserve, penalty, direct, and total costs for the two wind farms. The total cost comprises the combined direct, reserve, and penalty costs corresponding to the scheduled power. The direct cost shows a linear relationship with scheduled power; as scheduled power rises, a larger spinning reserve becomes necessary, leading to increased reserve costs and consequently driving up the overall generation cost. The penalty cost decreases, albeit at a slower rate, as scheduled power increases. Similarly, the cost variations for solar power over/under-estimation against scheduled power are portrayed in Fig. 5. The yearly operating and maintenance costs for solar PV power plants align within a comparable range to those of onshore wind power plants30. Lognormal PDF parameters for solar irradiance are adopted from Ref.30 as well. Furthermore, the direct cost coefficients, penalty cost coefficient, and reserve cost coefficient for solar power are also referenced from literature30. Yet, using the chosen PDF parameters for solar irradiance, the overall cost of solar power doesn't follow a strictly upward trajectory.

Figure 4
figure 4

Variation in the cost of wind power relative to scheduled power for wind generators.

Figure 5
figure 5

Fluctuation in the cost of solar power versus scheduled power for solar PV units.

The solar PV plant's stochastic power output is shown as a histogram in Fig. 6. The solar PV system's scheduled electricity delivery to the grid is shown by the red dotted line. As previously said, the schedule power can be any amount of electricity that ISO and the owner of the solar PV firm mutually agreed upon. Figure 7 represents the stochastic power generated by the wind farm. The red dotted line represents the scheduled electricity delivery to the grid by the wind farm. Tables 1 and 2 presents the optimal scheduling of the ten-unit system with and without DSM respectively. The best, average and worst cost and average CPU time among 100 runs of solutions acquired from the proposed ECOA, COA and GWO with and without DSM are summarized in Table 3. It is observed from Table 3 that execution time for ECOA algorithm is lesser compared to COA, GWO, HSPSO, CCDE, FCEP and CFCEP. Reduced computation time enables more effective implementation of demand-side management strategies since quick response times are essential for implementing demand response programs, load shedding, or load shifting, contributing to improved demand-side management and grid reliability. Furthermore, faster computation facilitates better integration of variable renewable energy sources by adapting quickly to their inherent variability.

Figure 6
figure 6

Distribution of real power (MW) from solar PV.

Figure 7
figure 7

Distribution of real power (MW) from wind farm.

Table 1 Optimal Scheduling of the 10-Unit System to Minimize Operational Costs without DSM.
Table 2 Optimal Scheduling of the 10-Unit System to Minimize Operational Costs with DSM.
Table 3 Statistical analysis of optimization results for test system – I.

Figures 8 and 9 illustrate the cost convergence patterns obtained from the proposed ECOA, COA, and GWO algorithms, both with and without DSM.

Figure 8
figure 8

Characteristics of convergence in a 10-unit system without DSM.

Figure 9
figure 9

Characteristics of convergence in a 10-unit system with DSM.

It is apparent from Fig. 8’s convergence characteristics that the proposed ECOA algorithm achieves convergence after 163 iterations in the context of dynamic economic dispatch with demand-side management. In comparison, the conventional COA and GWO algorithms converge at the end of 167 and 172 iterations, respectively. It is evident from Fig. 8 that the convergence behavior indicates the proposed ECOA algorithm reaches convergence after 170 iterations, while the conventional COA and GWO algorithms converge at the conclusion of 173 and 180 iterations, respectively. The findings suggest that the convergence of the proposed ECOA was not only swift but also exhibited a smoother trajectory compared to COA and GWO. Table 3 reveals that the operational cost is minimized when dynamic economic dispatch incorporates Demand-Side Management (DSM), as opposed to dynamic economic dispatch without DSM. Furthermore, the cost derived from the proposed ECOA remained the most economical among all methods. The achieved mean cost value closely approached the minimum, showcasing ECOA's competence in reaching global optimal solutions. Moreover, Table 3 reveals that the proposed ECOA algorithm exhibits a lower standard deviation in comparison to COA and GWO. This reduced standard deviation suggests greater stability and consistency in the performance of the ECOA algorithm, highlighting its potential for reliable and predictable outcomes. Due to space limitations, results acquired from COA and GWO53 cannot be given here. A sensitivity analysis was performed based on 100 trial test runs. Table 4 displays the results of the sensitivity analysis conducted for the proposed ECOA algorithm applied to Test System I and II. The results lead to the conclusion that a population size of 30 for the provided test system yields the global optimum for the test system—I. Consequently, the simulation outcomes firmly support the conclusion that the ECOA algorithm, as introduced in this study, holds significant potential for delivering high-quality solutions when contrasted with alternative algorithms.

Table 4 Sensitivity Analysis for the test systems I and II.

Test system: II

This system comprises twenty thermal power plants, two similar wind power generation units, two equivalent solar photovoltaic (PV) facilities, and two pumped-storage hydroelectric plants. The data for this test system are derived by mirroring the information from test system 1. Notably, the power demand in this configuration is twice that of test system 1. Specifically, hours 11, 12, and 13 represent peak load periods. During Demand-Side Management (DSM), 10% of the load during the 11th, 12th, and 13th hours is shifted to the 2nd, 3rd, and 4th hours. The optimal scheduling of the 20-unit system with and without DSM respectively for the 20-unit system is presented in Tables 5 and 6. Tables 5 and 6 provides an analysis of the best, average, and worst costs and average CPU time for 100 runs of solutions obtained from the proposed ECOA, COA and GWO with and without DSM. Table 7 reveals that the computational time for the ECOA algorithm is notably shorter than that of COA, GWO, HSPSO, CCDE, FCEP, and CFCEP algorithms. This accelerated computational speed enables swift decision-making in the face of dynamic system conditions, including abrupt shifts in demand or renewable energy generation. Additionally, the proposed ECOA algorithm adeptly harnesses available renewable energy while safeguarding system stability, thereby optimizing the equilibrium between conventional and renewable generation. Furthermore, faster algorithms may require fewer computational resources, making them more efficient and cost-effective for implementation on various hardware platforms, including embedded systems or edge devices.

Table 5 Optimal scheduling of the 20-unit system to minimize operational costs without DSM.
Table 6 Optimal scheduling of the 20-unit system to minimize operational cost with DSM.
Table 7 Statistical analysis of optimization results for test system—II.

Figures 10 and 11 show the cost convergence characteristics obtained from planned ECOA, COA, and GWO54 with and without DSM respectively. It is evident from Fig. 10’s convergence characteristics that the proposed ECOA algorithm achieves convergence after 221 iterations, while the conventional COA and GWO algorithms converge at the end of 226 and 230 iterations, respectively. It is noted from Fig. 11 that the convergence characteristics indicate the proposed ECOA algorithm achieves convergence after 193 iterations in the scenario of dynamic economic dispatch with demand-side management. In contrast, the conventional COA and GWO algorithms converge at the conclusion of 215 and 217 iterations, respectively. According to the findings, the proposed ECOA's convergence characteristic was faster and smoother than those of COA and GWO. Table 6 showcases that the inclusion of DSM results in lower costs compared to scenarios without DSM. Furthermore, among all the approaches, the proposed ECOA exhibits the most economical cost. The achieved mean cost value was close to the lowest value. Table 7 illustrates that, in comparison to COA and GWO, the proposed ECOA algorithm displays a diminished standard deviation. This decrease in standard deviation implies enhanced stability and consistency in the performance of the ECOA algorithm, underscoring its potential for delivering reliable and predictable outcomes. This proves that ECOA has the efficacy to create a global optimal solution. The findings from COA and GWO cannot be presented here due to space restrictions. The outcomes of the sensitivity analysis for the proposed ECOA algorithm on Test System I and II are presented in Table 4. These results indicate that, for Test System II, a population size of 45 results in the global optimum. Based on the simulation outcomes, it is evident that the ECOA algorithm proposed in this study possesses a greater probability of generating superior-quality results compared to alternative algorithms.

Figure 10
figure 10

Characteristics of convergence in a 20-unit system without DSM.

Figure 11
figure 11

Characteristics of convergence in a 20-unit system with DSM.

Conclusion and future research directions

The current study introduces an enhanced Cheetah Optimizer Algorithm that addresses the unpredictability of renewable energy sources and the involvement of pumped-storage hydroelectric units. This enhancement serves as a practical solution for real-life Distributed Energy Dispatching (DED) scenarios, both with and without Demand-Side Management (DSM). The proposed ECOA, COA and GWO are used to resolve two test systems. Optimization results indicate that the operational expenses associated with Demand-Side Management (DSM) are lower compared to those incurred without its implementation. Furthermore, research indicates that the introduced ECOA algorithm surpasses COA and GWO in performance metrics. The proposed ECOA approach exhibits adaptability and reliability, making it a viable solution for tackling multi-objective energy management challenges within a microgrid, especially when integrating demand response mechanisms. Future endeavors will involve exploring the capabilities of the enhanced cheetah optimization algorithm in addressing multi-objective optimization problems that encompass constraints. This investigation will specifically focus on navigating the trade-offs between conflicting objectives and constraints. Additionally, there is an opportunity to delve into hybridization with other optimization techniques, aiming to enhance convergence speed and improve solution accuracy. The suggested ECOA can be analyzed for its application in realizing the multi-objective optimal operation of bipolar DC microgrids. Furthermore, the suggested ECOA can be applied to elucidate the multi-objective dynamic optimal power flow problem in multi-microgrid systems which involve the integration of electric vehicles and renewable energy sources.