Giant Trevally Optimization Approach for Probabilistic Optimal Power Flow of Power Systems Including Renewable Energy Systems Uncertainty

Abstract

In this study, the Giant Trevally Optimizer (GTO) is employed to solve the probabilistic optimum power flow (P-OPF) issue, considering Renewable Energy Source (RES) uncertainties, achieving notable cost reduction. The objective function is established to minimize the overall generation cost, including the RES cost, which significantly surpassing existing solutions. The uncertain nature of the RES is represented through the employment of a Monte Carlo Simulation (MCS), strengthened by the K-means Clustering approach and the Elbow technique. Various cases are investigated, including various combinations of PV systems, WE systems, and both fixed and fluctuating loads. The study demonstrates that while considering the costs of solar, wind, or both might slightly increase the total generation cost, the cumulative generation cost remains significantly less than the scenario that does not consider the cost of RESs. The superior outcomes presented in this research underline the importance of considering RES costs, providing a more accurate representation of real-world system dynamics and enabling more effective decision making.

1. Introduction

In recent years, the power grid has seen a significant surge in the incorporation of RESs, notably wind and photovoltaic systems [1,2]. These sources offer significant advantages, such as reduced pollution and the conservation of resources. However, their integration into power systems has introduced new challenges due to growing uncertainties in load demands and variations in power production. Consequently, effectively managing these sources is pivotal for their successful integration into the network, ensuring a safe and financially viable energy market [1,2,3]. To achieve this, it is crucial to account for the probabilistic nature of WE and PV units while integrating such intermittent power sources into the grid. By taking these factors into consideration, power system operators can successfully handle the fluctuations and uncertainty of these renewable sources and ensure a reliable and cost-effective power supply.

Optimal Power Flow has been instrumental in the power system community for almost 50 years [4]. It is a multidimensional issue that seeks to optimize a particular objective function while adhering to operational and physical limitations set by the equipment and power system demand. The objective function, whether single or multi-objective, encompasses parameters like fuel costs, grid losses, emissions, and the voltage security index. The OPF incorporates inequality and equality limitations that concern the equilibrium equations of power and restrictions on control and state variables. Control variables encompass bus voltage, real power, and transformer tap ratios. State variables include imaginary power outputs, network line flow, and bus load voltages. That is why the OPF is considered a crucial tool for electrical utilities, as it outlines reliable and economically beneficial operating conditions for power systems. However, solving the P-OPF issue has become increasingly complex due to the integration of renewables and the network’s size and complexity, making it computationally demanding.

The deterministic OPF (D-OPF) was initially addressed using conventional mathematical programming methods, which proved to be effective in demonstrating their feasibility [5]. Despite being commonly used to effectively solve the D-OPF issue, conventional optimization techniques, like Newton’s technique [6], gradient projection technique [7], linear programming approach, and interior point approach [8], are reported to be associated with various challenges [9]. As optimization problems continue to evolve, new techniques encompassing artificial intelligence, as well as the metaheuristic search-based optimization approaches were designed to tackle the D-OPF problem. Recent efforts focused on search-based optimization approaches, which include the genetic algorithm (GA) optimization method [10], particle swarm optimizer (PSO) method [11,12], differential evolution optimization method [13,14], enhanced genetic algorithms optimization method [15], gravitational searching algorithm (GSA) method [16,17], multi-phase searching optimization algorithm [18,19], improving colliding bodies method [20], improved PSO method [21], biogeography-based optimizing approach [22], fuzzy-based hybrid PSO method [23], black-hole optimization approach [24], imperialist competitive optimization algorithm [25], harmony search optimization algorithm [26], PSO hybrid with GSA method [27], grey wolf optimization technique [28], and bee colony optimization approach [29]. Additionally, many multi-objective functions have been introduced for the D-OPF in [30,31].

The unpredictable nature of integrating intermittent RESs into the system has made solving the P-OPF problem considerably challenging where the metaheuristic approaches can provide best optimal solutions [32]. Moreover, solving the P-OPF problem requires computationally intensive and impractical simulations to consider a vast range of possible operational conditions. Although both the traditional and the intelligence-based approaches have been used as deterministic approaches for addressing the D-OPF approach in the previous paragraph, probabilistic methods have been shown to offer more satisfactory solutions and higher precision, especially in the presence of uncertainties [3,33]. Many of these techniques have been employed for tackling the P-OPF issue within power systems containing a significant presence of PV and WE installations. In [34], the two-point estimation method (2PEM), derived from the moment’s method, was employed to calculate the PDF of the photovoltaic output power. Nevertheless, it is important to note that the method of moments can sometimes produce results that fall outside the parameter range, making them unreliable as sufficient statistics. In [35], the Cornish–Fisher expansion method was introduced to tackle uncertainties associated with PV systems. Yet this technique may not offer accurate estimations for problems with an intricate structure and continuous return functions [36]. In [37], a P-OPF method for wind power plants was presented, and its PDF was evaluated using a heuristic technique. However, it is worth noting that obtaining accurate results for the PDF necessitates the availability of real data. A P-OPF issue, which explores the effects of wind speed’s high dimensional dependencies on P-OPF, was introduced in [38]. In this particular method, the Kernal density estimation technique was employed for estimating the wind speed PDF. However, this method can be impractical due to its dependence on the initial position of the bins and the exponential increase in the number of bins as dimensions increase. In [39], a comprehensive analysis was conducted on a model with multi-objectives within stochastic and non-convex OPF incorporating wind power. The study utilized the Latin hypercube sampling (LHS) technique to examine various aspects, including mean, variance, and skewness. Nevertheless, it is important to note that the Latin hypercube sampling (LHS) method may encounter a lack of statistical independence among sample points and may not demonstrate significant superiority over a random sampling of practical computational analysis. In [40], the unscented transformation method is introduced to tackle uncertainties in Renewable Energy Sources and load demand forecasts. In [41], a model predictive control approach, complemented with the autoregressive integrated moving average (ARIMA) prediction method, is utilized to forecast environmental and load variations in standalone PV/battery systems, addressing the intermittent nature of Renewable Energy Sources. As described in [36,37], the MCS and its variants were employed to determine the WE system’s power output PDF. A P-OPF issue was designed in [42,43] according to the MCS method of the power system that includes PV and WE units.

As stated earlier, integrating wind and solar energy in power systems presents critical challenges because of their variable characteristics. Typically, solar PV farms or wind facilities are privately owned, and the independent system operator (ISO) or the grid engages in a contractual agreement to buy the power scheduled from these operators. However, due to the unpredictable nature of these RESs, their power output might sometimes exceed the scheduled power, resulting in an underestimation of the accessible supply. In such cases, the ISO incurs additional costs, often termed ‘penalty costs’. This is because any overproduction that is not effectively utilized can go to waste [44]. Conversely, overestimation occurs when the power produced falls short of the levels scheduled previously. To balance demand of power, it is imperative for the ISO to maintain a spinning reserve, which subsequently drives up the operating expenses of the system [45]. This paper’s objective function integrates the uncertainty costs of RESs in conjunction with the generation expenses of thermal power units.

This study provides a comprehensive approach to the P-OPF issue, taking into account the uncertain nature of the RES and the various costs associated with this uncertainty. The key contributions made in this paper are as follows:

• Implemented and solved the classical optimal power flow problem using the Giant Trevally Optimizer (GTO) approach. The obtained results were then compared with other established algorithms, such as GA, SFO, and HHO, as mentioned in the paper, to confirm the efficiency and robustness of the proposed approach.
• Integrated an RES into the OPF problem to create the P-OPF problem, taking into account the uncertainty inherent in these sources. The P-OPF problem was tackled using the GTO approach, which was instrumental in handling the complexity and uncertainty of the problem.
• The uncertainty of the RES was addressed using the Monte Carlo simulation method, enhanced by the Elbow and K-means clustering methods. This innovative approach reduced the computational time and improved the precision of the simulation.
• Introduced the concept of three types of costs associated with RES uncertainty: direct cost, reserve cost, and penalty cost. These costs were incorporated into the P-OPF problem to provide a more accurate representation of real-world energy systems.
• Collected and analyzed actual historical data for the summers of four years (2018, 2019, 2020, and 2021). The data were used to simulate a more accurate 24 h summer day, with curve fitting performed for each hour of data for PV and WE using the Beta and Weibull Probability Density Functions (PDFs), respectively.
• The P-OPF problem, including the cost of RES uncertainty, was solved for two system scales: a 30-bus system and a 118-bus system. This was performed for both fixed and variable loads, providing a comprehensive view of the impact of RES costs on different system sizes and load conditions.

This article is organized as follows: Section 2 delves into the mathematical modelling and problem definition. The modeling of stochastic wind/solar power, along with the associated uncertainties, are detailed in Section 3. The proposed Giant Trevally Optimizer algorithm is introduced in Section 4. Section 5 presents the simulation and subsequent results. To conclude, Section 6 wraps up the paper by offering final insights and reflections.

2. Mathematical Modelling and Problem Formulation

Wind, solar, and thermal generators play a significant role as primary generation elements in power systems and are crucial in the P-OPF problem. The inherent uncertainties in these Renewable Energy Sources (RESs), due to their fluctuating nature, result in considerable variations in generation power. This underlines the necessity for a precise simulation of PV and Wind Turbine generators, discussed subsequently. Another crucial consideration, addressed in this section is the integration of the associated costs of these RESs into the mathematical model. Thus, the section at hand introduces a mathematical model that encapsulates these cost considerations with a goal to minimize total generation cost. Meanwhile, a thorough discussion on addressing inherent uncertainties is reserved for Section 3.

2.1. Modelling of Thermal Generators’ Costs

Conventional power generators, which rely on fossil fuels for operation, exhibit a quadratic relationship between the fuel expenses ($/h) and the power output (MW) at hour $h$ as follows [18]:

$$C_{Tj,h}\left(P_{TG,h}\right)={\displaystyle \sum }_{j=1}^{N_{TG}}{a}_j{P}_{TGj,h}^2+b_j{P}_{TGj,h}+c_j$$

(1)

In the equation, $a_j$, $b_j$, and $c_j$ denote the cost coefficients for the $j^{\mathrm{th}}$ conventional generator, which generates a power output of $P_{TG,h}$ at hour $h$, and $N_{TG}$ signifies the total count of the thermal generators.

2.2. Modelling the Direct Cost Associated with Wind/Solar Power

Unlike conventional thermal power units, wind and solar PV units operate without the need for fuel. If the ISO is the owner of the wind and solar PV power plants, the cost function might not be necessary because these energy sources do not demand fuel. This is unless the ISO chooses to assign a repayment cost to the initial investment for the PV or wind power plants or considers this as maintenance and renewal expenses [46]. In contrast, when solar PV or wind facilities are under private ownership, the ISO offers compensation corresponding to the contracted scheduled power.

The direct cost related to the $j^{\mathrm{th}}$ wind power plant at hour $h$ is represented through a function of the scheduled output power:

$$C_{d,w,i,h}=D_{w,i}{P}_{ws,i,h}$$

(2)

where $D_{w,i}$ denotes the coefficient of the direct cost tied to $i^{\mathrm{th}}$ wind generator, while $P_{ws,i,h}$ denotes the power scheduled at hour $h$ from the same generator.

Likewise, the direct cost linked to the $m^{\mathrm{th}}$ solar plant at hour $h$ is calculated as follows:

$$C_{d,s,m,h}=D_{s,m}{P}_{ss,m,h}$$

(3)

$D_{s,m}$ denotes the coefficient of the direct cost tied to the $m^{\mathrm{th}}$ solar plant, and $P_{ss,m,h}$ represents the scheduled output power at hour $h$ from the same solar power plant.

2.3. Modelling Costs of the Uncertainties Associated with Wind Power

In certain cases, the actual power output from the wind farm might fall short of the estimated value, leading to an overestimation of power from this uncertain source. The ISO must maintain spinning reserve capacity to ensure uninterrupted supply to consumers. The expense associated with allocating reserve power units to compensate for the overestimated amount is referred to as the reserve cost [47]. This reverse cost at hour $h$ is formulated as follows:

$$\begin{array}{l}{C}_{Rw,j,h}\left(P_{ws,j,h}-P_{wav,j,h}\right)=G_{Rw,j}\left(P_{ws,j,h}-P_{wav,j,h}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =G_{Rw,j}{\int }_0^{P_{ws,j,h}}\left(P_{ws,j,h}-p_{w,j,h}\right)f_w\left(p_{w,j,h}\right)d{p}_{w,j,h}\end{array}$$

(4)

where $G_{Rw,j}$ signifies the coefficient of the reserve cost linked to the $j^{\mathrm{th}}$ wind plant, while $P_{wav,j,h}$ represents real power available from the same plant at hour $h$. The term $f_w\left(p_{w,j,h}\right)$ denotes the PDF for wind power specific to the $j^{\mathrm{th}}$ wind power plant and hour $h$. Section 3.2 elaborately covers the procedure for calculating the probability of power output at different wind speeds.

In contrast to a situation of overestimation, there can be instances within the system where the wind farm’s actual power surpasses estimates. In these cases, the output of the RES is underestimated. In such instances where power output from renewable sources exceeds the scheduled amount, it can create a surplus. If this surplus cannot be offset by reducing power output from conventional generators, it may not be fully utilized. In such cases, the ISO may need to bear additional costs, often referred to as ‘penalty costs’, corresponding to the surplus amount [44]. This penalty cost can be formulated as follows:

$$\begin{array}{l}{C}_{Pw,j,h}\left(P_{wav,j,h}-P_{ws,j,h}\right)=G_{Pw,j}\left(P_{wav,j,h}-P_{ws,j,h}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =G_{Pw,j}{\int }_{ws,j,h}^{P_{wr,j}}\left(p_{w,j,h}-P_{ws,j,h}\right)f_w\left(p_{w,j,h}\right)d{p}_{w,j,h}\end{array}$$

(5)

where $G_{Pw,j}$ signifies the cost coefficient for penalties linked to the $j^{\mathrm{th}}$ wind facility, and $P_{wr,j}$ denotes rated power output derived from the same plant.

2.4. Modelling the Uncertainty Costs Associated with Solar Power

Similar to wind plants, solar plants exhibit uncertain and intermittent output. The strategy for handling the overestimation and underestimation of solar power should fundamentally align with that of wind power. For ease of calculation, the models for penalty and reserve costs are constructed according to the concept outlined in reference [48]. Section 3.3 further elaborates on the detailed approach used to compute the probability of power output at different solar irradiance.

The reserve cost related to the $k^{\mathrm{th}}$ solar power facility at hour $h$ is presented as follows:

$$\begin{array}{l}{C}_{Rs,k,h}\left(P_{ss,k,h}-P_{sav,j,h}\right)=G_{Rs,k}\left(P_{ss,k,h}-P_{sav,k,h}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=G_{Rs,k}{\int }_0^{P_{ss,k,h}}\left(P_{ss,k,h}-p_{s,k,h}\right)f_s\left(p_{s,k,h}\right)d{p}_{s,k,h}\end{array}$$

(6)

In the above equation, $G_{Rs,k}$ denotes the coefficient of the reserve cost linked to the $k^{\mathrm{th}}$ solar plant, while $P_{sav,k,h}$ stands for the actual power available from the same plant at hour $h$.

$P_{ss,k,h}$ refers to the scheduled power production from the solar power facility at hour $h$.

The penalty cost of the $k^{\mathrm{th}}$ solar plant at hour $h$ is presented as follows:

$$\begin{array}{l}{C}_{Ps,k,h}\left(P_{sav,k,h}-P_{ss,k,h}\right)=G_{Ps,k}\left(P_{sav,k,h}-P_{ss,k,h}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=G_{Ps,k}{\int }_{ss,k,h}^{P_{sr,k}}\left(p_{s,k,h}-P_{ss,k,h}\right)f_s\left(p_{s,k,h}\right)d{p}_{s,k,h}\end{array}$$

(7)

In this equation, $G_{Ps,k}$ refers to the coefficient of the penalty cost linked to the $k^{\mathrm{th}}$ solar plant.

The coefficients of the uncertainty costs associated with PV and wind power can be found in

Table 1. The coefficients of uncertainty costs associated with PV and wind power.

Cost Coefficient	Value
$D_w$	1.75
$D_s$	1.6
$G_{Rw}$	3
$G_{Pw}$	1.5
$G_{Rs}$	3
$G_{Ps}$	1.5

[44].

2.5. The Objective Function

The P-OPF problem’s target function can be determined by integrating all of the various cost functions previously covered.

$$\begin{array}{l}{F}_T=C_{T,h}\left(P_{TG,h}\right)+{\displaystyle {\displaystyle \sum }_{j=1}^{N_{WG}}}\left[C_{d,w,j,h}\left(P_{ws,j,h}\right)+C_{Rw,j,h}\left(P_{ws,j,h}-P_{wav,j,h}\right)\right.\\ \hspace{1em} \left.+C_{Pw,j,h}\left(P_{wav,j,h}-P_{ws,j,h}\right)\right]+{\displaystyle {\displaystyle \sum }_{k=1}^{N_{SG}}}\left[C_{d,s,k,h}\left(P_{ss,k,h}\right)+C_{Rs,k,h}\left(P_{ss,k,h}-P_{sav,k,h}\right)\right.\\ \hspace{1em} \left.+C_{Ps,k,h}\left(P_{sav,k,h}-P_{ss,k,h}\right)\right]\end{array}$$

(8)

The objective function is recalculated on an hourly basis to determine the optimal generation cost for each hour. This target function is bounded by the system’s equality and inequality constraints.

2.6. Equality Constraints

The power balance equations in power system optimization define equality restrictions. These stipulate that the generated real and imaginary power within the system must match the network’s demand and power losses.

$$P_{Gi}-P_{Di}=V_i{\displaystyle \sum }_{j=1}^{NB}{V}_j\left[G_{ij}cos\left({\delta }_{ij}\right)+B_{ij}sin\left({\delta }_{ij}\right)\right] \forall i \in NB$$

(9)

$$Q_{Gi}-Q_{Di}=V_i{\displaystyle \sum }_{j=1}^{NB}{V}_j\left[G_{ij}sin\left({\delta }_{ij}\right)-B_{ij}cos\left({\delta }_{ij}\right)\right] \forall i \in NB$$

(10)

In these equations, $P_{Gi}$ denote the real generated power at bus i, while $Q_{Gi}$ denotes the imaginary generated power at the same bus. $P_{Di}$ corresponds to the real demand power, while $Q_{Di}$ corresponds to the imaginary demand power. NB signifies the overall count of buses. $G_{ij}$ represents transfer conductance among buses $i$ and $j$, while $B_{ij}$ refers to the susceptance between the same buses. The term ${\delta }_{ij}$ $=$ ${\delta }_i$ $-$ ${\delta }_j$ expresses the difference in voltage angle among buses $i$ and $j$. Finally, $V_i$ and $V_j$ designate the respective voltages at buses $i$ and $j$, respectively.

2.7. Inequality Constraints

The inequality limits pertain to operating restrictions of equipment and components within the electric grid, as well as the security constraints placed on the load buses and transmission lines.

The generator constraints:

$$P_{TG,i}^{min}\le P_{TG,i,h}\le P_{TG,i}^{max} ,i=1,\dots ,{\mathrm{N}}_{TG}$$

(11)

$$P_{ws,j}^{min}\le P_{ws,j,h}\le P_{ws,j}^{max} ,j=1,\dots ,{\mathrm{N}}_{WG}$$

(12)

$$P_{ss,k}^{min}\le P_{ss,k,h}\le P_{ss,k}^{max} ,k=1,\dots ,{\mathrm{N}}_{SG}$$

(13)

$$Q_{TG,i}^{min}\le Q_{TG,i,h}\le P_{TG,i}^{max} ,i=1,\dots ,{\mathrm{N}}_{TG}$$

(14)

$$Q_{ws,j}^{min}\le Q_{ws,j,h}\le P_{ws,j}^{max} ,j=1,\dots ,{\mathrm{N}}_{WG}$$

(15)

$$Q_{ss,k}^{min}\le Q_{ss,k,h}\le P_{ss,k}^{max} ,k=1,\dots ,{\mathrm{N}}_{SG}$$

(16)

$$V_{Gi}^{min}\le V_{Gi}\le V_{Gi}^{max} ,i=1,\dots ,\mathrm{NG}$$

(17)

The security constraints:

$$V_{L_p}^{min}\le V_{L_p}\le V_{L_p}^{max} ,\mathrm{p}= 1,\dots ,\mathrm{NPQ}$$

(18)

$$S_{l_q}\le S_{l_q}^{max}, \mathrm{q}=1,\dots ,\mathrm{NTL}$$

(19)

Equations (11)–(13) define the limits of real power generated by thermal, wind, and solar power generators, respectively. Equations (14)–(16) specify the imaginary power capabilities of these generators. NG represents the overall count of generators (known as generator buses). Equation (17) sets the restrictions of the voltage of generator buses. Equation (18) outlines limits on the voltage of the load buses (also known as PQ buses). $\mathrm{NPQ}$ denotes the number of PQ buses. Equation (19) provides the line capacity constraints, which apply to the overall number of transmission lines (NTL) in the network.

It is critical to recognize that the convergence of power flow towards an optimal solution inherently ensures compliance with equality restrictions established by power equilibrium equations. Within inequality restrictions, both the generator’s real power (excluding the slack generator) and the generator bus voltage are considered as constrained control variables. The proposed algorithm picks a viable value within the variable’s range for each of the control variables. However, inequality constraints concerning the slack generator’s power, voltage constraints at the load bus, the imaginary output power of remaining power units, and line capabilities require particular attention, where their values should not exceed limits.

3. Modeling Stochastic Wind/Solar Power along with the Associated Uncertainties

Given that solar power and wind generation hinge upon solar irradiance and wind speed, respectively, fluctuations in the wind speed may be represented through a PDF, like the Weibull probability distribution [49]. Conversely, solar radiation presents a complex layer of uncertainties. The Beta distribution function proves to be an effective tool for modeling the PDF of these solar radiation variations [50]. The following equations represent the Weibull PDF of wind speed $f\left(V_w\right)$ and the cumulative distribution function (CDF) $F\left(V_w\right)$, respectively.

$$f\left(v\right)=\left(\frac{K}{C}\right){\left(\frac{v}{C}\right)}^{K-1}{e}^{-{\left(\frac{v}{C}\right)}^{K }},v>0$$

(20)

$$F\left(v\right)=1-e^{-{\left(\frac{v}{C}\right)}^{K }},v>0$$

(21)

where, $v$ stands for the current wind speed. While K and C denote the shape and scale parameters of the Weibull PDF, respectively.

The PDF of solar irradiance, denoted as $f\left(R\right)$, is represented using the Beta distribution function as follows:

$$f\left(R\right)=\frac{\mathsf{\Gamma }\left(\alpha +\beta \right)}{\mathsf{\Gamma }\left(\alpha \right)\mathsf{\Gamma }\left(\beta \right)}{R}^{\alpha -1}{(1-R)}^{\beta },$$

(22)

In this context, $R$ represents the current solar irradiance. The parameters β and α represent shape components of the Beta distribution. The symbol Γ(.) indicates the incomplete gamma function, while R stands for solar irradiance, which is quantified in W/m².

In this study, the historical data was leveraged for wind speed and solar irradiance during the summer seasons across four years (from 2018 to 2021), as reported in [51]. The data was utilized to emulate a precise 24 h summer day, with curve fitting performed for each hour across that day across the four years using either the Weibull for wind speed or the Beta probability distribution functions for solar irradiance and obtain the best-fit values of scale and shape parameters for each hour that represent the historical data [3]. The wind speed data was gathered at the height of 10 m, but as per product datasheet of Enercon E82-E4, the wind turbine, characterized by a 3 MW rating, usually has its hub height at 84 m. Therefore, to acquire the speed of wind at this elevation, the Weibull distribution extrapolation method is used as expressed below [3]:

$$C_2=C_1{\left(\frac{h2}{h1}\right)}^n$$

(23)

$$k_2=k_1\frac{1-0.0881ln\left(\frac{h_1}{hr}\right)}{1-0.0881ln\left(\frac{h_2}{hr}\right)}$$

(24)

$$n=\frac{0.37-0.0881ln\left(C_1\right)}{1-0.0881ln\left(\frac{h_1}{hr}\right)}$$

(25)

where $C_1$ and $K_1$ represent the Weibull old scale and shape parameters, while $C_2$ and $K_2$ represent the new scale and shape parameters, and $hr$ represents a 10 m reference height.

3.1. Modeling Wind and Solar Power

After fitting the data using Weibull and Beta PDF for each hour, a Monte Carlo Simulation (MCS) is run with 8000 scenarios to obtain the fitting curves of each hour and the scenarios obtained are used to determine the PDF of output power for solar and wind generators using the following equations.

The output power from wind generators is computed using the subsequent equation:

$$P_{wav}\left(v\right)=\left\{\begin{array}{cc}0& v<v_{\mathrm{in}}\mathrm{or}v\ge v_{\mathrm{out}}\\ P_{wr}\frac{v-v_{\mathrm{in}}}{v_r-v_{\mathrm{in}}}& v_{\mathrm{in}}\le v<v_r\\ P_{wr}& v_{\mathrm{in}}\le v<v_r\end{array}\right.$$

(26)

Here, $v_{\mathrm{r}}$, $v_{\mathrm{in}}$and $v_{\mathrm{out}}$represent rated, cut-in, and cut-out speeds, respectively; and $P_{wr}$ represents rated power output from the wind plant.

The power output from solar generators is computed using the subsequent equation:

$$P_{sav}\left(R\right)=\left\{\begin{array}{c}{P}_{sr}\left(\frac{R^2}{R_{std}{R}_c}\right)\mathrm{for}0<R<R_c\\ P_{sr}\left(\frac{R}{R_{std}}\right)\mathrm{for}R⩾R_c\end{array}\right.$$

(27)

where $R_C$ and $R_{STD}$ denote a specific radiation point and the standard condition solar radiation, respectively; and $P_{sr}$ represents the rated power output from the solar power plant. This study uses parameters from the specification details of the Enercon E82-E4 product datasheet, with the wind turbine parameters established as follows: $v_{\mathrm{in}}$ = 3 m/s, $v_{\mathrm{out}}$ = 25 m/s, and $v_{\mathrm{r}}$ = 16 m/s. At the same time, the rated power output from each turbine $P_{wr}$ is set at 3 MW. For the photovoltaic generator, $R_C$ and $R_{STD}$ are set at 120 W/m² and 1000 W/m², respectively.

3.2. Calculating the Probability of Wind Power

To calculate the reserve and penalty costs in Equations (4) and (5), the probability of wind power $f_w\left(p_w\right)$ needs to be calculated. As observed from Equation (26), wind power as a variable is discrete within certain wind speed ranges. Specifically, the output power is equal to 0 when $v<v_{\mathrm{in}}$or $v\ge v_{\mathrm{out}}$ and equal to $P_{wr}$ when $v_{\mathrm{in}}\le v<v_r$. The probabilities of such distinct zones can be presented as follows, according to [52]:

$$\begin{gathered} f_w\left(p_w\right)\left\{p_w=0\right\}=f_w\left(v<v_{\mathrm{in}}\right)+f_w\left(v>v_{\mathrm{out}}\right)=F\left(v_{\mathrm{in}}\right)+\left(1-F\left(v_{\mathrm{out}}\right)\right) \\ \hspace{1em} =1-\exp \left[-{\left(\frac{v_{\mathrm{in}}}{c}\right)}^k\right]+\exp \left[-{\left(\frac{v_{\mathrm{out}}}{c}\right)}^k\right] \end{gathered}$$

(28)

$$\begin{gathered} f_w\left(p_w\right)\left\{p_w=p_{wr}\right\}=f_w\left(v_r<v<v_{\mathrm{out}}\right)=F\left(v_{\mathrm{out}}\right)-F\left(v_{\mathrm{r}}\right) \\ =\exp \left[-{\left(\frac{v_r}{c}\right)}^k\right]-\exp \left[-{\left(\frac{v_{\mathrm{out}}}{c}\right)}^k\right] \end{gathered}$$

(29)

The power output of a wind unit is linear within the range from cut-in to rated speeds. Although the probability of wind speed $f\left(v\right)$ is readily available, it is also necessary to compute the probability of power at every specific wind speed $f_w\left(p_w\right)$. In order to achieve this, the change of variables rule is applied. This rule states that if the PDF of one random variable $f\left(v\right)$ is known, then the PDF of another random variable $f_w\left(p_w\right)$ can be calculated by considering the relationship between these two variables as follows:

$$f_w\left(p_w\right)=f\left(v\right)\times \frac{dv}{dp}$$

(30)

$$f_w\left(p_w\right)=\frac{k\left(v_r-v_{in}\right)}{c^k\ast p_{wr}}{\left[v_{in}+\frac{p_w}{p_{wr}}\left(v_r-v_{in}\right)\right]}^{k-1}\exp \left[-{\left(\frac{v_{in}+\frac{p_w}{p_{wr}}\left(v_r-v_{in}\right)}{c}\right)}^k\right]$$

(31)

3.3. Calculating the Probability of Solar Power

To calculate the reserve and penalty costs in Equations (6) and (7), the probability of solar power $f_s\left(p_s\right)$ needs to be calculated. As observed from Equation (27), the solar power output is continuous, but it is divided into two regions: from 0 to $R_C$ and greater then $R_C$. Here, the probability of solar irradiance $f\left(R\right)$ is known, so the change of variables rule is applied to calculate the probability of solar output power $f_s\left(p_s\right)$ in each region, such that

$$f_s\left(p_s\right)=f\left(R\right)\times \frac{dR}{dp}$$

(32)

Region 1: $0<R<R_c$

$$f_s\left(p_s\right)=f\left(R\right)\times 0.5\times \sqrt{\frac{R_c\times R_{std}}{p_{sav}\times P_{sr}}}$$

(33)

Region 2: $R⩾R_c$

$$f_s\left(p_s\right)=f\left(R\right)\times \frac{R_{std}}{P_{sr}}$$

(34)

It is important to note that the scheduled power for PV or wind can take any value between zero and the rated power according to the contract between the ISO and the private PV or wind power plant owner. When the power available falls below or exceeds the power scheduled, the uncertainty costs, whether reserve or penalty costs, must be paid. So, the scheduled power must be chosen properly. Here, the scheduled power is chosen to be the power that has the most probability of occurrence in each hour, as illustrated in Figure 1.

The rated power of the wind plant is selected to be 30 MW and 500 MW for the 30- and 118-bus systems, respectively. The rated power of the solar plant is chosen to be 20 MW and 250 MW for 30- and 118-bus systems, respectively [3]. Figure 1 offers a representation of the PDF for the output power of the wind plant integrated into the IEEE 30-bus system at hour 2. In the figure, curve (1) represents the relationship between the power and the speed of the wind, curve (2) shows the Weibull PDF of wind speed, and curve (3) shows the probability distribution of output wind power corresponding to different wind speeds. Figure 2 provides an example of the probability distribution of output power for the solar plant in the 30-bus system at hour 12. In the figure, curve (1) represents the relationship between solar power and solar irradiance, curve (2) shows the Beta PDF of solar irradiance, and curve (3) shows the probability distribution of output solar power corresponding to different solar irradiance levels. From curve (3) in Figure 1 and Figure 2, the scheduled power is the power that has the highest probability of occurrence in each hour.

3.4. Monte Carlo Simulation

In the previous sections, it was highlighted that the data of summer for the past 4 years was depicted by a corresponding 24 h day, applying curve fitting techniques for every hourly segment to obtain the PDF that represents each hour such that the Weibull PDF represents each hour in the wind data and the Beta PDF represents each hour in the solar data. To effectively integrate wind (WE) and photovoltaic (PV) energy sources into the power system, it is necessary for the system to handle every possible scenario due to the uncertain nature of RESs. Therefore, the use of a stochastic technique becomes indispensable, superseding the deterministic approach. In this research, the chosen technique was the Monte Carlo Simulation (MCS). The MCS technique was selected as it is deemed the most appropriate approach for managing large and intricate systems that contain a substantial quantity of uncertainties [53]. This simulation ran 8000 scenarios for each hour, which were subsequently included in the analysis. Figure 3 provides an example of the Weibull PDF for wind speed with MCS running for 8000 scenarios at hour 13. Figure 4 shows an example of the Beta PDF for solar irradiance with MCS running for 8000 scenarios at hour 11.

3.5. Proposed Clustering Technique

The MCS is performed with 8000 scenarios for each hour for wind and solar data. Therefore, a huge number of scenarios need to be considered in the P-OPF problem, which is a complex and time-consuming process. To decrease the overall count of scenarios and iterations, the MCS is integrated with the K-means clustering approach [3,54]. The K-means clustering technique condenses the 8000 scenarios per hour into a smaller number of equivalent clusters that accurately represent these scenarios [3]. However, one limitation of the K-means Clustering method is its initial random assignment of the number of clusters, necessitating the use of a cluster optimization technique to ascertain the optimal number of clusters [55]. This is where the Elbow method is employed. The method determines the optimal number of equivalent clusters for the dataset by performing K-means clustering for assumed cluster quantities ranging from 2 to 10. For each of these quantities, a corresponding sum of squared errors (SSE) is computed. The SSE values are then plotted against the assumed number of clusters. In this plot, the ‘elbow’ point, representing the most noticeable change, indicates the ideal number of clusters [56]. This is because beyond this point, the error increases only marginally despite an increase in the number of clusters. From our analysis, the most distinct Elbow point occurred when the number of clusters was 3, leading us to conclude that the optimal number of clusters for our study is 3. More comprehensive understanding and graphical representations of the K-means clustering and Elbow method can be found in [3].

4. Proposed Optimization Algorithm

The Giant Trevally Optimizer (GTO) approach was chosen to address the P-OPF issue across various scenarios. The GTO technique is a metaheuristic algorithm that derives inspiration from the hunting behaviors of the giant trevally [57]. To facilitate understanding, the GTO’s flowchart is provided in Figure 5, with comprehensive descriptions provided for each step of the algorithm. The giant trevally uses techniques that involve patterned foraging movements, selection of an optimal hunting region, and leaping out of the water to capture prey. The GTO algorithm replicates these strategies into a three-step process: extensive search, selection of area, and attack.

4.1. Extensive Search

The GTO technique simulates the long distances giant trevallies travel to locate food using a mathematical model based on Levy flights, a form of random walk. This phase improves the algorithm’s exploration capability and aids in avoiding local optima. The equation used in this phase can be illustrated as depicted below:

$$\begin{array}{c}X\left(t+1\right)={\mathrm{Best}}_P\times R+\left(\mathrm{Maximum}-\mathrm{Minimum}\right)\\ \times R+\mathrm{Minimum}\times Levy\left(\mathrm{Dim}\right)\end{array}$$

(35)

where $X\left(t+1\right)$ is the giant trevally’s position vector in the next iteration,${\mathrm{Best}}_P$ represents the best position obtained, $R$ represents a random number ranging from 0 to 1, and $\mathrm{Levy}\left(\mathrm{Dim}\right)$ represents the Levy flight.

4.2. Choosing Area

In this phase, the algorithm finds the optimum hunting region based on food availability inside the search space. The following equation is used to mathematically replicate this behavior:

$$X\left(t+1\right)={\mathrm{Best}}_P\times \mathcal{A}\times R+\mathrm{Mean}\_\mathrm{Info}-Xi\left(t\right)\times R$$

(36)

where $\mathcal{A}$ is a parameter that controls position change, $Xi\left(t\right)$ denotes the current position, and $R$ is a random integer. The effective utilization of all data obtained from prior locations is implied by the term Mean_Info for these giant trevallies.

4.3. Attacking

The algorithm’s last stage simulates the trevally’s attack on the prey. Light refraction, which affects the trevally’s eyesight, influences its behavior. To imitate this behavior, the algorithm calculates the visual distortion $\mathcal{V}$ using Snell’s equation, then simulates the trevally’s assault using:

$$X\left(t+1\right)=\mathcal{L}+\mathcal{V}+\mathcal{H}$$

(37)

where $X\left(t+1\right)$ denotes the next position, $\mathcal{L}$ is the launch speed, $\mathcal{V}$ is the visual distortion, and $\mathcal{H}$ is the leaping slope function, thereby allowing the algorithm to move from the phase of exploration to the phase of exploitation.

5. Simulation and Results

This research proposes an innovative method for addressing OPF and P-OPF issues by utilizing the capabilities of the GTO algorithm. The proposed GTO-based solution was evaluated for efficacy and viability on 30 and 118 bus systems, the details of which are displayed in

Table 2. The characteristics of the used standard test systems.

Features	Value	Detail	Value	Detail
Test system	30-bus	[63]	118-bus	[64]
Branches	41		186
Generators	6	[3]	54	[3]
Transformers	4	Branches: 11, 12, 15, and 36	9	Branches: 8, 32, 36, 51, 93, 95, 102, 107, and 127
Number of loads	21		99
Loads (MVA)	283.4 + j126.2		4242 + j1438
losses (MVA)	5.28 + j23.14		132.86 + j783.79
Bus voltage restrictions	[0.9–1.1]		[0.9–1.1]

[53]. To demonstrate the effectiveness of the suggested technique, it was compared to several well-established algorithms for classical OPF, including Sunflower Optimization (SFO) [58], Genetic Algorithm (GA) [59] and Novel Harris Hawk Optimization (HHO) [60]. Upon affirming the efficacy of the suggested technique, RESs were incorporated, and the P-OPF problem was addressed through different scenarios.

Table 3. Details of the case studies.

Case Study	Details
1	Address the classical OPF issue (without RES integration) for fixed loads, and assess the efficacy of the suggested approach in comparison to the selected alternative methods.
2	Integration of PV at the designated optimal bus location. Initially, the P-OPF is solved without considering the cost associated with PV, and then it is repeated with the cost included. Perform this case under both fixed and variable load conditions, then compare the results.
3	Integration of wind into the optimal bus location. Solve the P-OPF first without considering the cost associated with PV, then repeat with the cost included. Perform this case under both fixed and variable load conditions, then compare the results.
4	Integration of wind and PV into optimal bus location. Solve the P-OPF first without considering the cost associated with PV, then repeat with the cost included. Perform this case under both fixed and variable load conditions, then compare the results.

summarizes the specific details of these instances. This section presents the results of these scenarios, with all of them employing the GTO algorithm, and offers accompanying explanations. Each scenario illustrates the influence of integrating RES into the grid, following two strategies: one disregards the cost associated with RES uncertainties, while the other takes this cost into account. A comparative analysis of the results derived from these two strategies is also presented. The RESs are integrated into the system at optimal locations [61], as outlined in

Table 4. The optimal bus locations of RES in different systems.

Optimal Bus	30-Bus	118-Bus
PV	4	114
Wind	21	15

. The case studies are executed under both fixed load conditions and variable load conditions. The term fixed load refers to the standardized load corresponding to the abovementioned systems. On the other hand, the variable load, defined as a percentage of the standard test system’s load on a regular summer weekday, is outlined in reference [62]. Corresponding visualizations of these loads can be found in reference [3].

The computational framework of all proposed algorithms, including the OPF and P-OPF procedures, was realized within the MATLAB environment. Performance assessments of the simulations were executed on a system equipped with an Intel(R) Core (TM) i7-5500U CPU @ 2.40 GHz 2.40 GHz, 8 GB RAM, operating on a 64-bit operating system. The parameters for the GTO method, similar to any other metaheuristic optimization technique, were chosen through a systematic process of trial and error. This involved running several independent trials, then evaluating the algorithm’s performance to determine the optimal parameter values.

5.1. Case 1: Evaluating GTO Efficacy in Classical OPF

The proposed method was used in this case study to solve the classical OPF issue under fixed load conditions. The effectiveness of this approach was compared to well-known algorithms, such as SFO, GA, and HHO. With no RES fed into the system, the major objective was to optimize the thermal generators’ operating costs. Each algorithm was tested with specific parameters. The population size for the 30-bus system was set at 15 with a limit of 2000 iterations, whereas the population size for the 118-bus system was increased to 30 with the iteration limit remaining at 2000. Figure 6 and Figure 7 show a visual comparison of the objective function progress across the various methods for the 30- and 118-bus systems, respectively. As visibly substantiated by these figures, the GTO technique shows its superiority by obtaining the lowest operational cost and attaining this objective with fewer iterations than competing solutions.

5.2. Case 2: Adding the PV Power Plant to the Optimal Bus Location

For the IEEE 30-bus system, a 20 MW PV system is integrated at the optimal location, as indicated in Table 4. The population size is limited to 15 with a maximum of 200 iterations. The P-OPF is solved under fixed and variable load conditions, initially without accounting for the cost associated with the PV system and subsequently with these costs included. For both constant and fluctuating loads, the integration of PV units into the grid resulted in a reduction in total costs between hour 7 and hour 20, coinciding with the period of solar irradiance availability. Figure 8 illustrates the impact of integrating a solar power plant into the IEEE 30-bus system on the overall generation cost under two scenarios: one that overlooks the cost associated with PV power and another that takes it into account, both for a fixed load scenario. Similarly, Figure 9 presents the corresponding results for the scenario with variable loads.

The PV system with a rating of 250 MW is integrated at the optimal bus location in the IEEE 118-bus system. The population size is limited to 30, with a maximum of 400 iterations. Following a similar methodology as the previous case, the P-OPF is solved under the same scenarios. These calculations are carried out twice: once considering the cost associated with PV and once disregarding it. Figure 10 visually demonstrates the influence of incorporating a solar power plant into the 118 bus system on the overall power generation cost under two scenarios: one that overlooks the cost associated with PV power and another that takes it into account, both for a fixed load scenario. Similarly, Figure 11 presents the corresponding results for the scenario with variable loads.

The outcomes of this case clearly indicate that introducing the cost of PV does increase the total generation cost. However, the integration of PV into the grid still results in an overall decrease in generation costs. Incorporating the cost of PV power provides a more accurate representation and better reflects the realities of the system.

5.3. Case 3: Adding Wind Power Plant to the Optimal Location

For this case study, a wind system with a rating of 30 MW is integrated at the optimal bus location of the 30-bus system, as indicated in Table 4. The population size is limited to 15 with a maximum of 200 iterations. The P-OPF is executed under the same scenarios as those utilized in a prior case study. Figure 12 illustrates the impact of integrating a wind power plant into the IEEE 30-bus system on overall generation cost under two scenarios: one that overlooks the cost associated with wind power and another that takes it into account, both for a fixed load scenario. Similarly, Figure 13 presents the corresponding results for the scenario with variable loads.

The wind system with a rating of 500 MW is integrated at the optimal bus location in IEEE 118-bus system. The population size is limited to 30 with a maximum of 400 iterations. The P-OPF is executed under the same scenarios as those utilized in prior case studies. Figure 14 visually demonstrates the influence of incorporating a wind power plant into a 118-bus system on the overall power generation cost under two scenarios: one that overlooks the cost associated with wind power and another that takes it into account, both for a fixed load scenario. Similarly, Figure 15 presents the corresponding results for the scenario with variable loads.

The outcomes from this case clearly demonstrate that incorporating the cost associated with wind power does elevate total generation cost. Nevertheless, the incorporation of wind power into the grid continues to yield an overall reduction in generation costs. Incorporating the cost of wind power ensures a more accurate representation and better reflects the realities of the system.

5.4. Case 4: Adding Both PV and Wind Power Plant to the Optimal Bus Location

In this study, the capacity for PV and wind systems were specifically chosen to maintain an equivalent penetration percentage of RESs in both the 30-bus and the 118-bus systems. For each system, the RES penetration was set at 17.6%. This was achieved by incorporating a PV capacity of 20 MW and a wind capacity of 30 MW for the 30-bus system, which has a total load of 283.4 MW. Similarly, for the 118-bus system with a total load of 4242 MW, the RES penetration was maintained by integrating a PV capacity of 250 MW and a wind capacity of 500 MW. This consistent RES penetration percentage across different system sizes allows for more reliable comparisons and insights into the impacts of RES integration and cost considerations on overall power system operations.

For the IEEE 30-bus system, the wind and PV systems are integrated at the optimal bus location, as indicated in Table 4. The population size is limited to 15 with a maximum of 200 iterations. The P-OPF is solved in line with the same scenarios described in the prior case studies. Figure 16 demonstrates the impact of the insertion of wind and solar power plants into the IEEE 30-bus system, reflecting on total generation cost under two distinct scenarios: one excluding the costs associated with wind and PV power plants and another considering these costs, both in a fixed load situation. Concurrently, Figure 17 represents equivalent results but in the context of variable load scenarios. Complementing these figures,

Table 5. Percentage cost reduction comparison with and without considering the cost of RES for 30-bus (fixed load).

Hour	1	2	3	4	5	6	7	8	9	10	11	12
RES without cost	2.97%	2.91%	3.05%	3.04%	3.10%	2.87%	2.67%	2.74%	3.64%	5.06%	6.55%	7.51%
RES with cost	0.68%	0.64%	0.69%	0.66%	0.70%	0.61%	0.64%	0.87%	1.24%	1.64%	1.88%	2.02%
Hour	13	14	15	16	17	18	19	20	21	22	23	24
RES without cost	7.56%	7.66%	7.39%	6.61%	5.75%	4.77%	3.64%	2.95%	2.63%	2.69%	2.74%	2.91%
RES with cost	2.02%	2.03%	1.52%	1.61%	1.81%	1.57%	0.78%	0.89%	0.62%	0.61%	0.63%	0.64%

and

Table 6. Comparison of cost reduction percentage with and without considering the cost of RES for 30 Bus (variable load).

Hour	1	2	3	4	5	6	7	8	9	10	11	12
RES without cost	3.20%	3.24%	3.46%	3.57%	3.55%	3.30%	2.89%	2.78%	3.55%	5.12%	6.58%	7.46%
RES with cost	0.28%	0.29%	0.33%	0.35%	0.33%	0.28%	0.28%	0.37%	0.79%	1.51%	1.86%	1.98%
Hour	13	14	15	16	17	18	19	20	21	22	23	24
RES without cost	7.64%	7.49%	7.40%	6.56%	5.75%	4.91%	3.70%	2.88%	2.68%	2.61%	2.58%	3.05%
RES with cost	2.04%	1.86%	1.52%	1.39%	1.63%	1.57%	0.62%	0.63%	0.50%	0.37%	0.16%	0.01%

provide a detailed breakdown of the cost reduction percentages under the same scenarios for a more comprehensive comparison.

For the IEEE 118 bus system, the PV and wind systems are integrated at the optimal bus locations. The population size is limited to 30 with a maximum of 400 iterations. The P-OPF is executed following the same scenarios as described in the preceding case studies. Figure 18 visually demonstrates the influence of incorporating PV and wind power plants into the IEEE 118-bus system on overall power generation cost under two scenarios: one that overlooks the cost related to both wind and PV power and another that takes it into account, both for a fixed load scenario. Similarly, Figure 19 presents the corresponding results for the scenario with variable loads. Complementing these figures,

Table 7. Comparison of cost reduction percentage with and without considering the cost of RES for 118 Bus (fixed load).

Hour	1	2	3	4	5	6	7	8	9	10	11	12
RES without cost	4.09%	4.11%	4.26%	4.29%	4.31%	4.08%	3.67%	3.59%	4.36%	5.84%	7.56%	8.72%
RES with cost	3.82%	3.84%	3.99%	4.01%	4.03%	3.81%	3.43%	3.39%	4.12%	5.50%	7.10%	8.18%
Hour	13	14	15	16	17	18	19	20	21	22	23	24
RES without cost	8.75%	8.82%	8.77%	8.01%	7.04%	6.13%	4.94%	3.99%	3.63%	3.74%	3.76%	4.07%
RES with cost	8.20%	8.26%	8.19%	7.49%	6.63%	5.78%	4.61%	3.75%	3.40%	3.49%	3.52%	3.80%

and

Table 8. Comparison of cost reduction percentage with and without considering the cost of RES for 118 Bus (variable load).

Hour	1	2	3	4	5	6	7	8	9	10	11	12
RES without cost	5.86%	6.49%	6.69%	7.03%	6.82%	6.69%	5.48%	4.57%	5.04%	6.17%	7.64%	8.74%
RES with cost	5.38%	5.98%	6.13%	6.44%	6.22%	6.15%	5.05%	4.28%	4.75%	5.82%	7.18%	8.20%
Hour	13	14	15	16	17	18	19	20	21	22	23	24
RES without cost	8.88%	8.82%	8.74%	8.29%	7.41%	6.45%	5.35%	4.38%	3.98%	4.02%	4.40%	5.46%
RES with cost	8.33%	8.26%	8.16%	7.76%	6.97%	6.07%	5.00%	4.11%	3.72%	3.75%	4.10%	5.06%

provide a detailed breakdown of the cost reduction percentages under the same scenarios for a more comprehensive comparison.

Incorporating both wind and solar power plants contributes to a more significant cost reduction compared to individual use. Even when the costs of these Renewable Energy Sources (RES) are included, there are still notable overall savings, providing a realistic reflection of system operation. Considering the cost of RES thus offers a more accurate portrayal of the true system dynamics.

6. Conclusions

This paper proficiently formulates and solves the P-OPF problem for a hybrid power system integrating photovoltaic and wind energy sources under both fixed and variable load conditions, achieving substantial reductions in the operational cost of RES-integrated power systems. This feat has been achieved by taking into account the inherent uncertainties of Renewable Energy Sources and utilizing curve fitting for historical data, along with the most suitable Probability Density Functions. This approach of including costs of uncertainties associated with RES, encompassing direct, reserve, and penalty costs, has provided a superior solution in terms of cost reduction compared to existing solutions. These results clearly show that although considering the costs of solar, wind, or both might increase the total generation cost, the cumulative generation cost remains substantially less than the scenario that does not consider the cost of RESs. For instance, with a RES penetration of 17.6% in the 30-bus system (20MW PV and 30MW wind) and the 118-bus system (250MW PV and 500MW wind), substantial cost savings were observed even after accounting for RES costs.

This outcome strongly underscores the superiority of our approach in the planning and functioning of modern power systems. As a future direction, it is suggested to consider the integration of energy storage systems and electric vehicles into the power system. Additionally, expanding the optimization framework to incorporate more objectives, such as system losses, emission reduction, voltage stability, and resilience, could provide a more comprehensive solution to power system planning and operation.