Economic/Environmental Optimal Power Flow Using a Multiobjective Convex Formulation

: This paper addresses the problem of economic/environmental optimal power ﬂow with a multiobjective formulation using a second-order conic programming (SOCP) optimization model. This problem formulation considers renewable energy sources (RES), fossil-fuel-based power generation units, and voltage control. The proposed SOCP model is a stochastic scenario-based approach to deal with RES and load behavior uncertainties. An ε -constrained algorithm is used to handle the following three objective functions: (1) the costs of power generation, (2) active power losses in the branches, and (3) the emission of pollutant gases produced by fossil-fuel-based power generation units. For comparative purposes, the SOCP model is also presented using a linearized formulation, and numerical results are presented using a 118-bus system. The results conﬁrm that changing the energy matrices directly affects the cost of objective functions. Additionally, using a linearized SOCP model signiﬁcantly reduces reactive power violation in the generation units when compared to the nonlinearized SOCP model, but also increases the computational time consumed.


Introduction
The optimal power flow (OPF) is a classical nonlinear (NL) optimization problem that aims to determine the operational point for a given electric power system (EPS) that would minimize objective function (OF), meet the electric demand of the system, and satisfy a set of physical and operating constraints.For this problem, the most common OF is to minimize the fuel cost of electric power generation; however, many alternative objectives can be studied to optimize the EPS operation, including, for example, modern electricity markets and renewable energy resources (RES) integration.In this regard, OPF tools must be adapted to solve more complex applications [1][2][3][4][5][6].
The EPS generation framework is undergoing continual development, with the intention to substitute conventional fossil-fuel-dependent generation facilities with sustainable and eco-friendly RES.However, the high RES penetration in the power system raises the complexity of the problem, due to the presence of uncertainties in the parameters that define the functioning of diverse types of RES [7].Thus, solving the economic/environmental OPF (EEOPF) problem is still a timely topic in power systems research, and one that could guide us toward environmental preservation by decreasing the emission of those greenhouse gases (GHG) which result from producing electric power with fossil fuels [8][9][10].In this study, probabilistic and stochastic approaches have been frequently used to include uncertainties in the EEOPF formulation, and thus obtain a more realistic (if complex) model [11][12][13].
Multiobjective approaches are being developed to solve the EEOPF problem by minimizing both the generation costs and the GHG emissions produced in energy power generation [14][15][16][17][18].For example, in [14], a multiobjective evolutionary algorithm was applied to a deterministic economic and environmental dispatch (EED) problem, minimizing not only the cost of generation, but also the number of active power losses of the system, in addition to the emission of pollutants from fossil-fuel-based thermoelectric generations (TG).
Metaheuristic algorithms are an alternative approach to providing solutions to the OPF problem when classical approaches are not able to deal with the complexity of the model.However, the quality of their solutions cannot be considered the optimal global solutions to the problem.On the other hand, the OPF problem can be expressed using a well-defined mathematical formulation and solved using classical optimization techniques.Moreover, various convexification techniques have been developed to determine the optimal global solution of the OPF problem, employing algebraic manipulations of the objective function and constraints [1,19].Among them, the second-order conic programming (SOCP) models are widely used, striking a balance between computational effort and the precision of the solutions [20][21][22].In addition, a convex formulation of the OPF problem enables it to be adaptable to the evolving challenges of modern EPS planning/operation [19,23,24].
Decentralized approaches are also an alternative, suited to addressing large and complex problems in power systems [25].For example, a decentralized model can be used to solve the OPF problem by considering electricity market environments, as is presented in [26][27][28], due to the improvement of security and privacy of each area of the EPS; this concept can also be extended to distribution systems [29].

Related Literature
The multiobjective probabilistic EEOPF model, outlined in [15], aimed to minimize the pairs of OFs: power generation costs vs. GHG emissions, power generation costs vs. active power losses in the transmission lines, and GHG emissions vs. active power losses in the transmission resulting from high RES penetration.The Nondominated Sorting Genetic Algorithm II (NSGA-II) was used to deal with the nonlinearities of the problem.The incorporation of uncertainties related to both electrical demands and RES was facilitated through a fast and effective point-estimate probabilistic approach.The operational status of the system was ascertained utilizing the conventional Newton AC power flow technique.A similar probabilistic multiobjective EEOPF problem has been formulated in [16], with the objective of minimizing the cost of generation and GHG emissions.Their version of the problem incorporates TG and wind generator (WG) resources, while also accounting for uncertainties by employing random variables to represent variations in electrical and environmental behaviors.A biogeography-based optimization algorithm was utilized to solve the problem, and the point estimate method was employed to analyze the uncertainty parameters.In [30], by contrast, the authors formulated an active and reactive EED, adopting a dynamic wind-thermal generation scheduling approach to minimize the generation costs and pollutant emissions in the EPS.The probability density functions of WGs were estimated based on the gaussian mixture model.A population -heuristic technique was proposed to solve the EED problem and improve accuracy during the optimization.In [11], a hybrid algorithm (composted using phasor particle swarm optimization and a gravitational search algorithm) was applied to solve the EEOPF problem, with particular consideration given to WG and photovoltaic (PV) generation in the system.National Renewable Energy Laboratory data were used to estimate the output of WG and PV, while Monte Carlo simulation was used to evaluate the statistical features of the EEOPF problem.In [12], a stochastic OPF-based active-reactive power dispatch problem was proposed, with consideration to the renewable generators.The problem was solved through an efficient differential evolution algorithm (DE) that aimed to minimize the generation costs of all units installed in the EPS.In [13], a modified quantum-behavior lightning search algorithm was applied to solve the EED problem of minimizing the emission cost of a given power system with TG, small hydroelectric generation (HG), and WG generations.A particle swarm optimization algorithm with an artificial neural network was used to capture the probability factor of wind speed.The proposed approach was applied in a small system (six-bus).In [26], the authors presented an iterative algorithm for the decentralized DC single-objective OPF problem, under various contingencies.Reference [28] proposed an MINLP formulation of the decentralized AC single-objective OPF problem, with consideration to prohibited operational zones within the generation units.Both [26,28] dismissed the presence of RES, minimizing just the power generation costs, and consumed more computational time than the centralized approach presented in the literature.

Research Gap and Novel Contributions
Based on the literature review presented below, it can be concluded that metaheuristic algorithms are frequently used to tackle the EEOPF problem.Nevertheless, the conventional approach needs to be further investigated and developed to adequately handle the intricate challenges associated with the problem, including nonconvexity, nonlinearity, and other related issues.Table 1 summarizes the main works from literature related to the EEOPF problem which consider an EPS approach (centralized/decentralized), operational uncertainty, voltage control, fossil-fuel-based generation, RES, OFs of costs of power generation, active power losses, and GHG emissions in the power system.This paper provides a multiobjective stochastic SOCP optimization model for the EEOPF problem, giving consideration to voltage control, the presence of fuel-based generation units, and RES.The voltage control was considered by optimally adjusting capacitor banks and transformer taps.This study analyzed three distinct OFs: (1) power generation costs, (2) active power losses, and (3) GHG emissions.The ε-constraint method was employed to deal with these objectives and present a representative subset of the Pareto set.Historical data of uncertain variables were used to generate a set of representative operational scenarios for RES and load behavior.Finally, a linearized version of the proposed SOCP model was used to validate the robustness of the proposed approach to solve the EEOPF multiobjective problem.The main contributions of this paper are presented as follows: 1.
Proposing a stochastic-scenario-based SOCP model to solve the EEOPF problem, including RES, fossil-fuel-based generation units, and voltage control.Moreover, this paper also presents a linearized version of the SOCP model for the problem.

2.
Developing a multiobjective method for the EEOPF problem, considering three objective functions.This approach used the ε-constraint method to obtain Pareto optimal solutions.3.
Analyzing and comparing the SOCP models (traditional and linearized) results to verify their efficacy.The results highlighted that the linearized SOCP model can significantly improve the precision of the reactive power subproblem, in regard to the traditional SOCP model.The reactive power dispatches were verified using the Newton-Raphson method of MATPOWER.
The remainder of this paper is organized as follows: Section 2 presents the modeling of the uncertainties and the linear and SOCP models for the EEOPF problem.Then, in Section 3, the results and discussion are presented.Finally, relevant conclusions are drawn in Section 4.

Mathematical Models
This section first describes the scenario generation strategy to define both RES and load behavior over a given planning horizon.Then, the RES power generation models were presented for use in PV-based and wind-generation units.Finally, the SOCP and linearized models used to solve the EEOPF multiobjective problem are defined.

Modeling Uncertainty
In power systems operations, the level of electrical demand and climatic conditions are among the main factors that contribute to the introduction of uncertainties in the OPF problem [1,2,31].The uncertain behaviors of electrical demand, wind speed, and solar irradiation are considered through historical measurements, and reduced into a treatable set of representative scenarios, labelled Ω [31].This methodology has been described in the following steps: 1.
Historical data of electrical demand (D), wind speed (W), and solar irradiation (I) were divided by their peak values, respectively, to obtain the normalized data.

2.
The data were sorted into descending order, according to the electric demand data, so as to trace curve D, maintaining the hourly correlation with W and I.

3.
Electric demand levels have a significant impact on the system's operation; in this regard, planning horizons were divided into time blocks, (t ∈ Ω t ) according to curve D.Then, the information of W and I in each time block was sorted in descending order, as presented in Figure 1.Note that this strategy loses the chronology of the data; however, this was not a significant issue for medium-term operation planning.

4.
Using the information of D, W, and I, normalized, in each time block, it was then possible to determine the cumulative distribution functions (CDFs) of each variable according to their respective time blocks.5.
We then divided the CDFs into segments.The mean value of the information in each segment was used to determine the value of the scenario (π).Therefore, the probability of the scenario (ρ) was determined according to the size of the segments.For illustrative purposes, Figure 2 illustrates this process for the CDF curve of W at t 3 , where horizontal and vertical dotted lines represent the ρ W t 3 and π W t 3 , respectively.In this case, we used three segments of 30%, 70%, and 100% for light, nominal, and heavy levels, respectively.The probabilities of the scenarios were 30%, 40%, and 30% for light, nominal, and heavy levels, respectively.6.
The subsets of stochastic scenarios (ss) in each time block were obtained according to Ω ss t = (π ss t , ρ ss t ), where ss ∈ (D, W, I), ∀(t ∈ Ω b ).Then, a one-by-one combi- nation of scenarios contained in the same time block was carried out as follows:  4. Using the information of D, W, and I, normalized, in each time block, it was then possible to determine the cumulative distribution functions (CDFs) of each variable according to their respective time blocks.5. We then divided the CDFs into segments.The mean value of the information in each segment was used to determine the value of the scenario (π).Therefore, the probability of the scenario (ρ) was determined according to the size of the segments.For illustrative purposes, Figure 2 illustrates this process for the CDF curve of W at  , where horizontal and vertical dotted lines represent the  and  , respectively.In this case, we used three segments of 30%, 70%, and 100% for light, nominal, and heavy levels, respectively.The probabilities of the scenarios were 30%, 40%, and 30% for light, nominal, and heavy levels, respectively.4. Using the information of D, W, and I, normalized, in each time block, it was then possible to determine the cumulative distribution functions (CDFs) of each variable according to their respective time blocks.5. We then divided the CDFs into segments.The mean value of the information in each segment was used to determine the value of the scenario (π).Therefore, the probability of the scenario (ρ) was determined according to the size of the segments.For illustrative purposes, Figure 2 illustrates this process for the CDF curve of W at  , where horizontal and vertical dotted lines represent the  and  , respectively.In this case, we used three segments of 30%, 70%, and 100% for light, nominal, and heavy levels, respectively.The probabilities of the scenarios were 30%, 40%, and 30% for light, nominal, and heavy levels, respectively.Then, a one-by-one combination of scenarios contained in the same time block was carried out as follows:  = ( ,  ,  ) ∀(,  ∈  ) . Figure 3 shows the  combination considering three different levels to D, W, and I, respectively.The probability of each scenario was determined using the product of the individual probabilities of each variable presented in the same time block, such as  =    , ∀(,  ∈  ).

RES Power Output
The active power injected by the PV and WG units was calculated according to (1) and (2).The probability of each scenario was determined using the product of the individual probabilities of each variable presented in the same time block, such as

RES Power Output
The active power injected by the PV and WG units was calculated according to (1) and (2).
Constraints ( 1) and ( 2) are linear models of the power output of PV and WG units, respectively.For each stochastic scenario, the power output of each renewable unit depends on its available primary energy font λ i and λ w .
It is worth mentioning that modern EPS systems can include energy storage devices to mitigate the effects of RES' uncertainties in the EPS operation [32,33].However, this study focuses on verifying the efficiency/performance of the SOCP model to solve the OPF problem in EPS that have high levels of RES penetration.Due to this, energy storage devices were not considered in the problem formulation.

Multiobjective EEOPF Problem: SOCP Formulation
The EEOPF problem was formulated through a multiobjective stochastic scenariobased SOCP formulation ( 7)- (25).The OF Ψ, as presented in (3), considered the following three objectives: the energy generation cost $ pg of all generation technologies available in the system (Γ G ), presented in (4); the creation of an Equation (5) to determine the active power losses cost (ploss) in the EPS.Finally, the equation outlined in (6), to determine the GHG cost (GHGC) in nature, calculated in TGs units [14].
In each t and stochastic scenario ss, the physical and operational constraints of the EEOPF problem were modeled.Equations ( 7) and ( 8) are the active/reactive power balances, respectively.
In this formulation, the square of the voltage and current magnitudes are represented with the variables v SQ k = v 2 k and i SQ km = i 2 km , respectively.Constraints ( 9) and ( 10) determined the reactive power injected by shunt compensators through a disjunctive formulation.Thus, if the binary variable σ sh k = 0, then q sh k = 0, otherwise q sh k = v SQ k b sh k .Constraints (11) and (12) were the magnitude and angular voltage limits in the bus.
Constraints ( 13)-( 15) described Kirchhoff's voltage law, considering transformers with on-load tap changers (OLTC).Equation ( 14) served as an approximation for the voltage angle difference, where θ km = (θ k − θ m ) and vk was the estimated voltage magnitude at the bus k.Conic constraint (15) defined the voltage drop calculation through use of the active/reactive power flow and the square variables of both voltage and current in the bus.Constraint (16)  This study gave particular consideration to a continuous formulation of the OLTC tap, wherein an ideal transformer's tap ratio (τ) and impedance were adopted in the modeling.A fictitious node was adopted between the ideal transformer and the transformer's impedance.Equation (17) describes the square voltage magnitude in the fictitious node (fn), while Equation ( 18) determines the voltage magnitude drop at buses k and m.
∀km ∈ Γ T The transformer tap ratio was defined based on the regulation percentage of the OLTC (∆ km ) installed in the branch τ km = ∆ km + 1, and could rewrite (18), as shown in (19).Equation ( 20) defined τ km in function of δ km and v k .Constraint (21) represents the operational limits of δ km , adopting the relaxing of ( 18) and the maximum regulation percentage in OLTCs .

Solution Technique: ε-Constraint Method
This study used the ε-constraint method to determine a Pareto set for the proposed multiobjective EEOPF model [35,36].The ε-constraint method aims to optimize one OF, while the others become constraints of the optimization model, as shown in (34).Where n is the quantity of OFs analyzed, − ε j is the upper bounds for j th ε-constraint, and X is the feasible region of decision variables, x.

s.t EEOPF constraints
OF j (x) ≤ − ε j ; ∀j ∈ {1, 2, . . . ,n}, j = i; x ∈ X Determining appropriate upper bounds for the ε-constraint formulation was a crucial task to avoid this approach missing a potentially viable solution close to the feasible region's boundary [36].The proposed approach determined some values for − ε j through the difference between the upper and lower bounds in each ε-constraint, as shown in (35).The lower bound of the j th ε-constraint (OF j ) was defined by solving the EEOPF problem with consideration given only to the OF j .On the other hand, the upper bound (OF j ) was estimated as being the highest value of the OF j after solving the EEOPF problem with consideration given to the other OFs, individually.The optimization model (35) was solved by adopting several values of the continuous parameter ε j and restricted among zero and one (0 ≤ ε j ≤ 1). Figure 4 describes the proposed methodology, where n is the number of objective functions to be analyzed.
OF j (x) ≤ OF j − ε j OF j − OF j ∀j ∈ {1, 2, . . . ,n}, j = i; 0 ≤ ε j ≤ 1 The EEOPF problem which considered the ε-constraint method and both SOCP models is presented in (36).14), ( 16)- (33) SOCP linearized model (35) ε-constraint Table 2 shows a brief comparison between the number of variables and constraints present in the SOCP models proposed in this work, where |.| is the cardinality of the set.
The first SOCP model presented three continuous (v SQ k , θ k , and q sh k ) and one binary (σ sh ) variables per bus, six variables (i SQ km , p km , q km , p mk , q mk , and δ km ) per branch, and two variables (p g , q g ) per generation unit.In addition, the SOCP model must respect six constraints ( 7)-( 12) per bus (two equality and four inequality); six constraints ( 13)-( 16) and (21) per branch (two equality, one quadratic, and two inequality); three constraints ( 22)-( 24) for TG, HG, and WG (two inequality and one quadratic); and two constraints (25) per PV unit (two inequality).On the other hand, the SOCP linearized model presented two more variables (∆ p km and ∆q km ) per bus, and eight constraints ( 26)-( 33) (five equality and three inequality) per branch, in addition to L blocks of piecewise linearization adopted in the modeling.However, when the stochastic approach was applied to the EEEOPF problem, each scenario analyzed in the model represented a single problem; the number of variables/constraints were then multiplied by all scenarios adopted.The EEOPF problem which considered the ε-constraint method and both SOCP models is presented in (36).

Tests and Results
The proposed models were implemented in AMPL [37] and solved using the commercial optimization solver CPLEX 20.1 with default settings.A 118-bus system from PGLib-OPF v19.05 [38] was used to validate the proposed models and the ε-constraint method.Numerical experiments were carried out on a computer with an Intel i7 processor of 3.2 GHz and 16 GB of RAM.
The generation costs can be found in [14,15,24,38,39], and the linear coefficients of WG and PV units were selected through the energy auctions [40].The minimum and maximum voltage limits were 0.95 and 1.05 p.u, respectively.The voltage regulation of all OLTCs of the system was ∆ km = 10%.
Two different energy matrices case in the EPS were analyzed: (i) Case E1: only TG present.
In this study, the WG VESTAS v80 2 MW and the PV panel GCL-P6/72 330 were used to provide a maximum power of 1 GW.The c ploss and c GHG were equal to 120.00 US$/MWh [41] and 45.00 US$/ton [42], respectively.The stochastic nature of the electrical demand was modeled on the weekly demand curve of the Midwestern region of the United States in 2021 [43].It was assumed, for the sake of simplicity, that all load buses exhibit identical consumption patterns.
The historical data of wind speed and solar radiation during 2021 in Illinois (Midwestern United States) were obtained from [44], with a time resolution of 8760 h.RES data were divided into four blocks, while power demand, wind speed, and solar irradiation were divided into three load levels each (light, nominal, and heavy).
Table 3 displays the environmental conditions that were adopted in the simulations.In addition, uniform meteorological and seasonal conditions were adopted for locations with RES units; based on previous experiments, we concluded that this would demonstrate acceptable system behavior for one year.The ε-constraint method employs ε values ranging from 0 to 0.9 in increments of 0.1, to construct Pareto curves.Step 1.0 was dismissed in the ε-constraint method because it cannot find OF values lower than the minimal OF.In Case E2, HG units were located at buses 4, 24, 26, 31, 40, 42, and 69; WG units were located at buses 1, 15, 19, and 56; PV units were located at buses 73, 91, 99, and 107; and TG units were in the remaining buses.
Table 4 presents the average values and the probabilities of demand scenarios adopted.Complete results are available in [45].

Results: SOCP Model
Table 5 presents the bound values obtained using individual OF j minimizations, which were then used in the ε-constraint method for the SOCP method.By analyzing the OFs bounds results, it is possible to verify that the RES introduction in the EPS reduced almost all bounds except in OF ploss , which rose 12.00% in Case E1.This increase was because the installation of RES units in the EPS did not adopt an optimal allocation strategy.   .This is expected, since the search space is affected by modifying the OFs boundaries, according to the ε-constraint method.
By analyzing the Pareto curves presented in Figure 5a for Case E1, it was possible to verify that  -or  -constraints significantly impact  $ when  = 0.7, increasing the  $ value more that 5.73%.By contrast, the value of  = 0.7 rose more than 4.50% of  $ value.The largest increase observed for the value of  in  $ was 16.55%, while for  , it was 8.63%.The Pareto curves presented in Figure 5b highlighted that the  -constraint did not modify the  values, with the highest increase in values equal to 0.6%.The  $constraint increased in 5.06% the  values when  $ = 0.4, and presented largest increase, equal to 24.57%.
The Pareto curves shown in Figure 5c revealed that  -constraint does have a significative impact on  values, where the largest increase obtained was 0.66%.On the other hand, the  $ -constraint increased the  values in 4.89% of cases when The horizontal axis represents the upper limits of each OF according to the values of ε.The solution of the SOCP model cannot determine a feasible solution for some values of ε, for OF ploss (b) and OF GHGC .This is expected, since the search space is affected by modifying the OFs boundaries, according to the ε-constraint method.
By analyzing the Pareto curves presented in Figure 5a for Case E1, it was possible to verify that ε ploss -or ε GHGC -constraints significantly impact OF $ pg when ε ploss = 0.7, increasing the OF $ pg value more that 5.73%.By contrast, the value of ε EGHG = 0.7 rose more than 4.50% of OF $ pg value.The largest increase observed for the value of ε ploss in OF $ pg was 16.55%, while for ε EGHG , it was 8.63%.
The Pareto curves presented in Figure 5b highlighted that the ε GHGC -constraint did not modify the OF ploss values, with the highest increase in values equal to 0.6%.The ε $ pg -constraint increased in 5.06% the OF ploss values when ε $ pg = 0.4, and presented largest increase, equal to 24.57%.
The Pareto curves shown in Figure 5c revealed that ε ploss -constraint does have a significative impact on OF GHGC values, where the largest increase obtained was 0.66%.On the other hand, the ε $ pg -constraint increased the OF GHGC values in 4.89% of cases when ε $ pg = 0.5, presenting the largest increase, equal to 35.64%.
The introduction of RES in the EPS (Case E2), shown in Figure 6a-c, highlighted that the OFs values were better distributed compared to curves presented in Figure 5.In addition, the RES generations reduced the fuel costs by 17.69%, 12.46%, and 18.54% in higher ε-constraint steps for minimizing the power generation costs, active power losses, and GHG emissions, respectively.
By analyzing the Pareto curves of Figure 6a, the OF $ pg value increased in 5.25% of cases when ε ploss = 0.5, with its largest increment equal to 26.11%.The OF $ pg value presented a rise of 3.33% when ε GHGC -constraint was equal to 0.7, whereas its maximum increment in OF $ pg value was 7.92%.The Pareto curves presented in Figure 6b indicated that ε GHGC -constraint rose in 4.57% the OF ploss values wherein ε GHGC = 0.5, with the highest increase in values equal to 23.96%.The ε $ pg -constraint resulted in an increase in 7.13% of the OF ploss values when ε $ pg = 0.4, whereas the largest increase was 40.86%.The Pareto curves shown in Figure 6c highlighted that the ε ploss -constraint rose in 4.40% the OF GHGC values when ε ploss = 0.5, with the largest increase obtained as 52.08%.On the other hand, the ε $ pg -constraint increased the OF GHGC values in 4.65% when ε $ pg = 0.7, and presented their largest increase as being 33.79%.
Table 6 shows the modification in the OF values regarding the lower bound of each OF analyzed.ε GHGC 4.65% (0.7) 33.79% (0.9) 4.40% (0.5) 52.08% (0.9) − − −: result did not present an increment because the OF analyzed was the same as the ε-constraint.X: the increment observed did not present a significative impact in the OF value.
The CPU times consumed by optimizing OFs (power generation costs, active power losses, and GHG emissions) were about 32, 16, and 27 s, respectively, for Case E1.On the other hand, the CPU times were less than 44, 17, and 35 s, for the same OFs, for Case E2.

Results: SOCP Linearized Model
The SOCP linearized model ( 7)-( 14), ( 16)- (33) was solved using L = 10 blocks for piecewise linearization.The number of blocks was chosen through empirical means to achieve a favorable balance between the computational effort required and the quality of the solution obtained.The voltage magnitude estimated and the ∆S km used in the simulations were based on the operational points present in PGLib-OPF v19.05 [38].
Table 7 presents the boundary values used in the ε-constraint method for the SOCP linearized method.Similar to the SOCP model, the upper bound for OF ploss rose 10.46% when RES generations are factored into the EPS.The other bounds for OFs likewise reduced their values when the RES are allocated in the EPS.The linearized model presented similar behaviors to those results obtained using the SOCP model by optimizing the  $ (a),  (b), and  (c).Similarly, this model could not determine a feasible solution for some ε values due to the reduction in the ε-constraints boundaries.The Pareto curves shown in Figures 7 and  8 were equivalent to those in Figures 5 and 6, respectively.The RES penetration reduced the TG costs by 17.72%, 12.59%, and 18.57% in the worst ε-constraint steps for the minimization of the  $ ,  , and  , respectively.The CPU times consumed for the SOCP linearized model for  $ ,  , and  were about 190, 78, and 144 s, respectively, for Case E1.On the other hand, for Case E2, the CPU times were about 207, 64, and 221 s, respectively.These CPU times are sensitive to the number of linearization blocks in the piecewise linearization.Note that each block increased the number of constraints of the problem.As a result, fewer linearization blocks could compromise the quality of the solutions obtained.
The developed conic models showed practically identical results, but significantly differed in regard to computational times.Therefore, the feasibility of reactive powers from the generating units was investigated.The MATPOWER Software [46] was utilized Similarly, this model could not determine a feasible solution for some ε values due to the reduction in the ε-constraints boundaries.The Pareto curves shown in Figures 7 and 8 were equivalent to those in Figures 5 and 6, respectively.The RES penetration reduced the TG costs by 17.72%, 12.59%, and 18.57% in the worst ε-constraint steps for the minimization of the OF $ pg , OF ploss , and OF GHGC , respectively.
The CPU times consumed for the SOCP linearized model for OF $ pg , OF ploss , and OF GHGC were about 190, 78, and 144 s, respectively, for Case E1.On the other hand, for Case E2, the CPU times were about 207, 64, and 221 s, respectively.These CPU times are sensitive to the number of linearization blocks in the piecewise linearization.Note that each block increased the number of constraints of the problem.As a result, fewer linearization blocks could compromise the quality of the solutions obtained.
The developed conic models showed practically identical results, but significantly differed in regard to computational times.Therefore, the feasibility of reactive powers from the generating units was investigated.The MATPOWER Software [46] was utilized to apply nodal voltage levels, active power from generating units, and control variable values, such as transformer taps and compensating shunts [24].The results indicated that the SOCP linearized model significantly reduced the infeasibility of reactive power supply for cases E1 and E2 by 76.07% and 64.48%, respectively, compared to the SOCP model.

Conclusions
This work presented a multiobjective second-order conic programming (SOCP) model to resolve the economic/environmental optimal power flow (EEOPF) problem in the EPS.Both of the SOCP models we developed were capable of efficiently solving the stochastic EEOPF problem while considering the objectives of minimizing OF $ pg , OF ploss , and OF GHGC , and could do so with satisfactory computational time.It should be noted that the SOCP model outperformed the SOCP linearized approach in all analyzed cases in terms of speed.Furthermore, the ε-constraint method efficiently determined Pareto fronts for each OF.Finally, for comparative purposes, a linearized model of constraint (15) was also proposed utilizing the reference, which used first-order Taylor expansion [20].
The results indicated that modifying a system's energy matrix by installing renewable energy sources (RES) reduced pollutant emissions and system generation costs.However, the values of active power losses in the network increased due to RES allocation issues and dependence on climatic conditions.This outcome was due to the lack of a prior study of the optimal allocation of RES in the system.Conducting a prior analysis to determine optimal RES allocation in the EPS would result in only a few instances of increased active power loss values.In addition, the results confirmed that conventional and linearized SOCP models can determine feasible Pareto sets of solutions for the EEOPF problem, whether or not they considered RES in the EPS with tight OF values, and they can do so with respect to all operational and control constraints.However, the linearized model presented higher computational times than the SOCP model due to the piecewise linearization.This drawback can be overcome by reducing the number of linearization blocks, but it will impact the accuracy and the computational time consumed by the computational implementation of the model.

Figure 3
shows the Ω ss t combination considering three different levels to D,W, and I, respectively.

Figure 2 .
Figure 2. CDF curve for the third t of wind speed factor.

Figure 2 .
Figure 2. CDF curve for the third t of wind speed factor.Figure 2. CDF curve for the third t of wind speed factor.

Figure 2 .
Figure 2. CDF curve for the third t of wind speed factor.Figure 2. CDF curve for the third t of wind speed factor.

Figure 4 .
Figure 4.An ε-constrained strategy for multiobjective models applied to the EEOPF problem.

Figures 5 and 6 Figure 5 .
Figures 5 and 6 show the Pareto curves obtained by optimizing the OFs for Cases E1 and E2, respectively.Energies 2023, 16, x FOR PEER REVIEW 13 of 21

Figures 7 and 8
Figures 7 and 8 show the Pareto curves by optimizing the OF $ pg (a), OF ploss (b), and OF GHGC (c).Horizontal lines represent the OFs values considering the ε-constraint method for both Cases.

Figure 8 .
Figure 8. Curves for the SOCP linearized model Case E2: (a) OF $ pg ; (b) OF ploss ; and (c) OF GHGC .The linearized model presented similar behaviors to those results obtained using the SOCP model by optimizing the OF $ pg (a), OF ploss (b), and OF GHGC (c).Similarly, this model could not determine a feasible solution for some ε values due to the reduction in the ε-constraints boundaries.The Pareto curves shown in Figures7 and 8were equivalent to those in Figures5 and 6, respectively.The RES penetration reduced the TG costs by 17.72%, 12.59%, and 18.57% in the worst ε-constraint steps for the minimization of the OF $ pg , OF ploss , and OF GHGC , respectively.The CPU times consumed for the SOCP linearized model for OF $ pg , OF ploss , and OF GHGC were about 190, 78, and 144 s, respectively, for Case E1.On the other hand, for Case E2, the CPU times were about 207, 64, and 221 s, respectively.These CPU times are sensitive to the number of linearization blocks in the piecewise linearization.Note that each

Table 1 .
Brief comparative analysis of the some works from literature.
ss)∈Ω ss t OF $ pg + OF ploss + OF GHGC hr t ρ ss g t,ss β g + α g c GHGC ∀t ∈ Ω t ; ∀ss ∈ Ω ss t defined the voltage angular opening limits in the branch.,ss + δ km,t,ss = 2(r km p km,t,ss + x km q km,t,ss ) + i SQ km z (14)m (θ km,t,ss ) = x km p km,t,ss − r km q km,t,ss(14)B ; ∀t ∈ Ω t ; ∀ss ∈ Ω ss t defined the active/reactive power flow limits in each l th block.

Table 2 .
Dimension of the proposed SOCP formulations.

Table 3 .
Illinois wind speed and solar irradiation levels.

Table 5 .
Bound values for the SOCP model.

Table 6 .
Summary of increments of SOCP results in relation to OF J .

Table 7 .
Bound values for the SOCP linearized model.
Shunt susceptance/conductance in the bus k [Ω] b km , g km Series susceptance/conductance in the branch km [Ω] b sh km , g sh km Shunt susceptance/conductance in the branch km [Ω] c 2 , c 1 , c 0 Quadratic/linear/constant cost coefficients of power generation US$/W 2 , US$/W, US$ c ploss , c GHGC Active power loss and greenhouse gas emission cost [US$/W, US$/ton] h Duration of each time block in hours [h] Slope of linearization power flow in l th block P D , Q D Active/reactive power demand in the buses [W, VAr] P, Q Lower active/reactive power injected limit [W, VAr] P,Q Upper active/reactive power injected limit [W, VAr] P W RW , P I RI Rated active power of the WG and PV [W] S km Upper apparent power flow limit in the branch [VA] S g Apparent power limit of a generation unit [VA] tg θ CAP , Capacitive/Inductive reactive power injection factors tg θ I ND V, V, Lower/Upper voltage magnitude limit [V] vk Voltage magnitude estimated in the node k [V] γ, β, α, Emission coefficients for thermal power generation ton/W 2 , ton/W, ton Nominal solar irradiation level W/m 2 λ W , λ W , λ W R Lower, upper, and nominal wind speed limits [m/s] ∆ km Regulation percentage of OLTC in the branch km Square of the current magnitude in the branch km [A] Active/reactive power injected by the generation technology g installed in bus k [W, VAr] p km , q km Active/reactive power flow in the branch km [W, VAr] Reactive power injection by the shunt compensator installed in bus k [VAr] v k , θ k Voltage magnitude and angle in the bus k [V, rad] v Square of the voltage magnitude in the bus k [V] δ km Auxiliary variable for the OLTC model ∆p km,l , Value of active/reactive power flow in l th block [W, VAr] ∆q km,l Binaries Variables: σ sh Status of the shunt compensator