Optimal Allocation and Sizing of Decentralized Solar Photovoltaic Generators Using Unit Financial Impact Indicator

: A novel ﬁnancial metric denominated unit ﬁnancial impact indicator (UFII) is proposed to minimize the payback period for solar photovoltaic (PV) systems investments and quantify the ﬁnancial efﬁciency of allocation and sizing strategies. However, uncontrollable environmental conditions and operational uncertainties, such as variable power demands, component failures, and weather conditions, can threaten the robustness of the investment, and their effect needs to be accounted for. Therefore, a new probabilistic framework is proposed for the robust and optimal positioning and sizing of utility-scale PV systems in a transmission network. The probabilistic framework includes a new cloud intensity simulator to model solar photovoltaic power production based on historical data and quantiﬁed using an efﬁcient Monte Carlo method. The optimized solution obtained using weighted sums of expected UFII and its variance is compared against those obtained by using well-established economic metrics from literature. The efﬁciency and usefulness of the proposed approach are tested on the 14-bus IEEE power grid case study. The results prove the applicability and efﬁcacy of the new probabilistic metric to quantify the ﬁnancial effectiveness of solar photovoltaic investments on different nodes and geographical regions in a power grid, considering the unavoidable conditional and operational uncertainty.


Introduction
Photovoltaics (PVs) and decentralized generators play a fundamental role in today's developing world, where securing access to easily approachable, sustainable, cost-effective, and clean energy is one of the main objectives of all countries.According to International Energy Agency, renewable energy sources (RESs) will provide 30% of electrical power demand in 2023, and PVs are expected to face the highest growth [1].Unfortunately, higher penetration levels of RESs are hard to achieve in practice due to the low flexibility of thermal power plants and the stochastic nature of environmental conditions affecting the availability of renewable power.The inherent variability of renewable power production can introduce economic and technological risks, and scientific research must focus on viable ways to maintain the power quality and profitability of the grid despite the addition of PV systems.Substantial investments in renewable energy require the development of more effective planning and integration strategies, careful identification of the most suitable sizes and locations of RESs, and rigorous quantification of unavoidable sources of uncertainty.
It is important to understand the impact of renewable resources on power grids and their contribution to sustainability goals.Different types of generation technologies [2], alternative allocation and planning strategies in deregulated electricity markets [3], and technological, financial, and social aspects [4] are examples of some of the most important themes discussed in this context.Decentralized generators can be optimally allocated based on their financial, sustainability, and safety perspectives [5,6], whether adopting modern optimization techniques or via metaheuristic approaches [7].Computational frameworks for the optimal allocation of PV units often account for economic and environmental aspects [8], electrical engineering quantities such as power losses, voltage deviations [9], voltage stability indices [10], and probabilistic reliability and availability metrics [11,12].Examples of economic objective functions for optimally allocating the decentralized generators include the net profit of wind farm investors [13] or the costs of power losses, power supply, social welfare, or load interruptions, see for instance [14,15].Compared to economic objectives, the technical aspects of RES integration on existing power grids are more commonly studied in the literature.Technical objectives used for the optimal allocation of decentralized generators are classified and compared in [16,17], and the authors show some of the most common objectives focus on the minimization of power losses, improvement of voltage profile, current reduction in weak lines, the spinning reserve power, and network MVA capacity.Many research papers proposed static and deterministic analyses only, potentially leading to an incomplete understanding of relevant dynamics and risks due to a lack of execution of the system evaluations on a time horizon [18,19].Similarly, some limitations must be faced by deterministic optimization methods, which can lead to sub-optimal and non-robust investment solutions [20,21].A comprehensive discussion of PV technology and various mitigation methods to support decentralized PV generations occurring in low-voltage distribution networks can be found in [22].
Reliable and robust investment plans on RESs must include a clear quantification of relevant uncertainties and analysis of the grid operations over time due to their stochastic and time-varying nature [23].Robust solutions are obtained thanks to a combination of including dynamic simulations, Monte Carlo (MC) simulations [24,25], optimal power flow (OPF) models, sensitivity analysis methods [26], and population-based or metaheuristic optimization approaches [27].More advanced probabilistic frameworks have also been proposed for a risk-based optimization of renewables [28], for the simulation of extreme weather conditions [29], or for the minimization of expected costs and expected energy-notsupplied values of multi-energy networks [30].In these studies, the running and investment costs are often considered for financial analyses and optimization, where the minimization of the overall running cost of the system alone can lead to investment solutions characterized by a very high penetration of renewable generators which guarantees a reduction in the system running cost because of the negligible cost of renewable energy production.However, maximizing the PV units can lead to instability on the grid, and it is not always cost-effective because the system running cost is a non-probabilistic and nonlinear function of the total investment in PV capacity.A lack of financial efficiency becomes particularly clear when PV units are installed at cloudy locations or in the proximity of low-capacity cables behaving as bottlenecks for the network.Moreover, if PV systems are placed on buses where the power demands are low, this can prevent a localized consumption of the produced energy, thereby increasing the transmission/distribution losses.
This paper proposes a probabilistic framework for the optimal and robust allocation of utility-scale PV units in transmission networks, adopting a novel metric that gives a direct quantification of the payback time for the investment [31].The probabilistic framework proposed allows quantifying the effect of uncertainty, such as the random failures of components, extreme weather events, and non-stationary factors, such as load demand and solar irradiance (obtained using a novel probabilistic cloud covering model calibrated using real-world historical data).
The proposed metric denominated unit financial impact indicator (UFII) measures the financial efficiency of the investment.Minimizing the expectation value of the UFII allows for minimizing the payback time, while the variance of the UFII distribution represents the stability of the financial return (or the risk associated with the financial investment).This allows the identification of the optimal sizing and siting of PV units that minimize the time needed for the return on the investment to achieve efficient spending of resources under uncertainty.

Power Optimization
A widely used tool for energy system cost optimization is the AC-OPF method [32].It considers bus voltage angles, magnitudes, generators' active power, and reactive power flow.In this paper, an AC-OPF approach is used to evaluate the operational and maintenance costs of the grid, which satisfies the load demand while obeying the technical limitations and requirements.A compact formulation of the proposed AC-OPF problem is defined as follows [33]: subject to linearized nodal power balance constraints, g(x; δ) = 0, and inequality constraints, h(x; δ) ≤ 0, see Equations (A1)-(A9) in the Appendix A. The grid operational and maintenance cost, C(x; δ), is optimized by tuning the grid variables, x = (Θ, V, p c , q c ) T for the environmental-operational state of the grid δ, which comprises the load demand profile, the irradiance, the PV generation profile, and the topological structure of the network.The grid variables are the voltage angles and magnitudes, Θ and V respectively, that must satisfy the inequality constraints on the voltage profile (related to the stability of the grid) and the active and reactive power p c and q c , produced by the n c controllable generators.Hence, the optimal system cost C(x ; δ) represents the minimum network cost obtained by applying the optimal operational profile x to the uncertain scenario δ.If PV units are installed on the grid, the net generated power p = p c + p pv (G) is the sum of the controllable generators' production and uncontrollable active power produced by the on-grid PV systems, p pv (G).The latter depends on the uncontrollable solar irradiation, G. Similarly, the net reactive power is given by q = q c + q pv (G), however, for simplicity, in this work, it is assumed that the PV units produce active power only and q = q c (unity power factor).Thus, the integration of PV units will affect the power factor of the whole system, and conventional power plants need to adapt their active and reactive power outputs accordingly to secure grid reliability.Power curtailment involves a deliberate reduction of the PV power output to satisfy the grid's thermal and stability constraints [34].However, power curtailment is not considered within the scope of this work and the power output from PV units is regarded as uncontrollable.For a more detailed description of the AC-OPF problem, the reader is reminded of Appendix A.

Modeling of the Nominal PV Generation
The power generated by a solar panel, P PV [W], with respect to the total irradiance hitting a tilted PV module, G module , is defined as follows: where P max represents the power output of the   C] and finally T C [ • C] is the PV cell temperature.This model is adopted for finding the PV power generation normalized to a unit PV capacity of 1 MWP.The mathematical expressions for calculating the total irradiance hitting a tilted PV panel, i.e., G module , are represented in Appendix B. It is important to notice that the effect of temperature on the output production of PV units is not modeled directly.However, the cloud model detailed in Section 4.1 includes the effect of the temperature on the PV unit production since it is calibrated using effective power generated by the PV units that depend on different environmental factors but lumped in the cloud coverage model.

Modeling of the Nominal Power Demand
In warm-climate countries such as Cyprus, the industrial activity and the power consumption in governmental, private-sector, and domestic buildings remain almost the same throughout the year.However, the usage of air conditioning (A/C) units varies according to the season.The A/C unit operations directly influence the electrical power demand by using electricity for cooling in summer and heating in winter.Mild temperatures observed in spring and autumn result in negligible usage of A/C units with the lowest power demand values in the year, while the highest load demands occur in August due to the highest mean daily temperatures according to climatological data, comprising measurements between 1991 and 2005, by the Department of Meteorology of the Republic of Cyprus, August [35].
Because the case study is static and has single values for both active and reactive loads, the power demand values over 8760 h are normalized with respect to the case study's load values which are taken as the annual peak loads.In addition, it is assumed that all days of the year have similar load profiles with two daily peaks, one at noon and one in the evening.The effects of public holidays are neglected to ease the modeling efforts.

Unit Financial Impact Indicator (UFII)
The UFII quantifies the monetary efficiency of the transmission-connected PV plants added to the system compared to a baseline case.This defines the normalized financial influence of the total PV capacity added on the power grid computed as the percent change in the total running cost of the grid per one unit of PV capacity added.Mathematically the UFII is defined as: where C o represents the optimal running cost of the original system (€/h); C pv the optimal running cost of the system after the addition of PV generators (€/h) and PV T denotes the total PV capacity added to the system (MWp).
The optimal running cost of the power grid without PVs is obtained by solving the AC-OPF as follows: whilst the optimal running cost after the addition of PV units is: where f i P and f i Q are the marginal cost functions of the real and reactive power production of the i-th controllable generator and include both fuel and O&M costs.f PV (PV i ) represents the operations and maintenance cost of the i-th PV unit allocated on the grid, which is a constant as it only depends on the PV allocation vector PV = PV 1 , . . ., PV n b .PV will be optimized in a later stage of the proposed framework.C o (x; δ) and C pv (x; δ) both depend on design variables and the (uncertain) operational-environmental scenario δ.To ease the notation, these explicit dependencies are omitted in the remainder of the paper.The UFII metric in Equation (3) can also be extended to include other factors such as tax support determinants, tax deductions, inflation/deflation effects, and budget allowances.However, these extensions go beyond the scope of this work and will not be further considered.
To incorporate the unavoidable uncertainty affecting the environmental-operational scenarios δ, the expected value of the UFII, E[UFI I], and its standard deviation σ[UFII] are used.UFI I(δ) quantifies the economic gain due to the reduced generation cost, and therefore, it can be regarded as an instantaneous normalized reward.In contrast, the E[UFI I] and the σ[UFI I] quantify the average monetary gain and its variability as a percentage of per-unit PV capacity installed.Hence, a higher UFII means a lower instantaneous cost for the specific scenario, a higher E[UFI I] means a shorter payback period and indicates a more effective sizing and positioning of PV units that lead to better financial prospects in the long run.It is important to notice that the σ[UFI I] represents the total variability of the UFII function, including the variabilities due to the noise of the process and the "natural" variability of the signal δ, e.g., due to seasonal and daily variability of the weather conditions and power demand.
It is assumed that the UFII(δ) is an ergodic stochastic process.This implies that the random process UFII will not change its statistical properties with time and that its statistical properties can be obtained from single, sufficiently long realizations (samples).Hence, the expected efficiency of an investment can be estimated by averaging the UFII values identified for different scenarios n mc as follows: where f δ represents the probability density function of scenario δ.Similarly, the risk of the investment is quantified by the standard deviation of UFI I and calculated as follows: Different scenarios can be easily generated by means of the Monte Carlo method by generating random events from their probability density functions, as explained in Section 4.

Modeling Uncertainty
Environmental-operational uncertainties affect the power grid state and the goodness of PV investments.Therefore, uncertainties due to varying consumer power demand; seasonal and daily changes in the solar irradiance, and thus, PV power generation, for instance, due to stochastic cloud movements and intensities; and random component failures such as line outages and generator shut-down events are considered in this work.Several probabilistic models have been introduced in the literature to characterize power demand variability and component outages.
Power demand: The consumer power demand naturally varies due to different environmental conditions, special occasions, user behaviors, etc.In this work, the power demand profiles are generated according to the guidance in Reference [36], which studies the consumption trends with respect to consumer behavior.More specifically, the daily load profiles are generated from a realistic point of view, with two peaks per day at noon and evening times where relatively higher power demand levels are expected [37].Please note that the profiles of the 24-h demand curves (each profile representing a day of the year) resemble each other and only differ in peak values due to seasonal variations in consumption.In addition, the random floatations around the nominal hourly demand are modeled as Gaussian distribution with a coefficient of variation P f , which is used to represent the size of the random variability in the consumer power demand.The variability of the power demand is assumed to be limited to within three standard deviations from the nominal profile.The application of P f ensures that a load demand curve generated is always unique and allows studying various possible combinations of power demand.
Component outage: Malfunctions of the power grid components occur in real life due to equipment aging, weather-induced failures, and malicious attacks.Homogenous Poisson processes are often used to model the time-to-failure of different components and duration of component outages, see, e.g., [38].Here, a generic component i is considered to have an expected number of occurrences λ i per year while the duration τ i of failure events, τ i , are assumed to be normally distributed with a mean value of µ τ i and a standard deviation value of σ τ i .For each event type, the hour of the year associated with the component outage is represented by a vector F i = (F i1 , . . . ,F i n ) where n represents the total number of hours of the i-th component failure in one year.Although there are no limitations for the number and type of failures that can be added to the power grid model, to prove that the model supports their consideration, only a subset of the most probable failures is considered, namely, two branches and one generator are assumed to be the least reliable components and the only components affected by failures.
Cloud coverage: The power production of a PV system is directly proportional to solar irradiance [39,40].Cloud coverage is modeled as a non-stationary Gaussian process that is fully defined by its second-order statistics.The differences between the clear sky irradiation [41] and historical measurements measured in Kalkanli, Cyprus [35], are used to estimate the mean value and standard deviation for the irradiance for each hour of the year.This is obtained assuming PV modules have a tilt angle of 30 degrees facing south.The real measurements of irradiance data consist of hourly global horizontal irradiance and direct normal irradiance measurements covering a period of four years between 2013 and 2017.Clearly, irradiance cannot be lower than zero, and it is assumed that the higher irradiance value corresponds to the theoretical clear sky irradiance.Hence, the distribution is truncated to satisfy these two constraints.For more details on the stochastic irradiance model, see Appendix B.
Although the main contribution to solar irradiance is due to cloud coverage, the model takes indirectly into account other environmental factors such as temperature, humidity, wind velocity, altitude, and air pressure.Therefore, the cloud coverage model can be considered an effective model of all environmental factors affecting solar irradiance but not explicitly modeled for simplicity's sake.
Heat waves: Extreme weather conditions such as heat waves and freezing temperatures increase the demand for electric power due to air conditioning and electrical heating, respectively.The expected number of days in a month with heat waves, λ E , is sampled from a non-homogenous Poisson distribution with λ E > 0 during the warmest two months of the year, i.e., July and August, and negligibly small otherwise [35].During a heat wave, both active and reactive power demands are increased by 10%.This is used to simulate the extra electrical energy demanded for the increased use of A/C units.The heat wave model can also be extended to include the effect of climate change and the expected increase of extreme events.

Uncertainty Quantification
A realistic analysis of the integration of PV units into the power grid requires a proper quantification of the effect of uncertainty.The Monte Carlo method is a flexible and easily applicable tool to model complex problems under the effect of uncertainty.Monte Carlo simulators are especially useful for power grid analyses, which are often complex and high-dimensional, to simulate variable power demands, random component failures, and volatile renewable production and for uncertainty [42,43].The simulator generates power demand, solar irradiance, and component failures from the model of uncertainty defined in Section 4.1; then, these values are used to evaluate the AC-OPF and UFII; and statistical information of the output is computed.
The Monte Carlo simulation procedure can be summarized as follows: 1.
Initialization: The grid structure, generator and lines parameters, and the vector of PV sizes and capacities (i.e., the investment) are defined.The stochastic model is also defined and includes the parameters and family of the probability distributions.

2.
Sampling: Realizations of weather conditions and operational variables are generated for a year with an hourly resolution where: a.
Loads at each node of the power grid are sampled hourly from a truncated normal distribution; b.
Cloud covering for each individual region is sampled from normal distributions.The differences with the clear sky irradiance model are then used to calculate the net solar irradiance hitting the solar PV generators; c.
The presence of extreme heat waves and their durations are sampled from a non-homogeneous Poisson process; d.
Component failures are samples from a Poisson distribution and their duration from a truncated normal distribution.

3.
Scenario selection: A subset of hours is randomly selected from the year.The scenarios indicated with δ t with t = 1, 2, . . ., n mc comprises solar irradiance, the presence of clouds and heat waves, and power demand at n mc a random hour within a year.4.
Grid performance evaluation: From the scenarios {δ 1 . . . ,δ t , . . . ,δ n mc } defined at point 3, the active and reactive loads for each node of the power grid are calculated at each time t as follows: a.
Normalized PV power output: The net solar irradiance hitting a PV module is computed by subtracting from the clear-sky irradiance the effect of the cloud coverage.Finally, the PV generations per one unit of PV capacity (MWp) are computed for each node.b.
Power demand: The active and reactive power demand values are updated considering the occurrence of extreme weather conditions sampled at point 2c.c.
Component outages: From the component outages sampled at point 2d, the condition of each component is identified, defining the configuration of the power grid.d.
Grid evaluation: The UFI I t is calculated according to Equation (3), and the optimal costs C o (δ t ) and C pv (δ t ) are found by solving two AC-OPFs, one by evaluating the grid with no PV installed and one for the selected PV investment strategy.

5.
Post-processing: The statistical indicators of the network performance are computed, i.e., the expected E[UFI I] and standard deviation σ[UFI I].
Note that the number of randomly chosen scenarios n mc must be sufficiently large to capture the behavior of UFII over the year without sacrificing accuracy in the estimation of the UFII expectation for much less computational time.To check the accuracy of the simulation, the V(E[UFI I]) is quantified.Specifically, the variance of the estimator is obtained by resampling, where ten independent Monte Carlo runs obtained from the n mc scenarios are used to estimate ten realizations of E[UFI I] and associated variance.This allows us to estimate V(E[UFI I]) at no additional computational cost.

Optimization under Uncertainty
An optimization algorithm is introduced in this work to prescribe an optimal investment PV * defined as a vector of power capacities on each node of the system.To find a robust PV allocation vector that compromises between the minimization of the payback time and maximization of the revenue stream stability, an objective function is defined by a weighted sum of the negative expected UFII and the weighted standard deviation.This yields the following optimization problem: subjected to 0 ≤ PV j ≤ b u j where PV is the design vector defining the PV capacities on all nodes of the grid, b u j is the maximum capacity of the node j and PV * is the optimal investment that minimizes the output of interest.The Monte Carlo simulator procedure is embedded within the optimization algorithm, and it is necessary to evaluate the objective function − Therefore, E[UFI I] and σ[UFI I] are the outcomes of the Monte Carlo simulation procedure depending on different candidate solutions, i.e., the PV allocation vectors (provided in step 1 in Section 4.2).Moreover, note that the objective is to maximize E[UFI I] combined with a minimization of the uncertainty (second term).A well-established generic algorithm is used to solve the optimization problem numerically, which has often been used on distributed generators [44] because the gradients are not computable and the state space is usually quite large.
Four cases are considered to minimize the risk and to maximize the expected return of the investment, which combines the values of E[UFI I] and σ[UFI I] via a robustness parameter α.With a value of α = 0, the objective function is equivalent to the maximization of the expected UFII without any consideration for the variability of the instantaneous UFII returns.With α > 0, the objective function includes a robustness element by controlling the variance of the UFII metric by penalizing solutions having highly uncertain returns for the investments, i.e., unstable revenue streams.Note that the proposed approach may be computationally very burdensome, especially for large systems.In fact, for each candidate allocation vector in the population of the evolutionary algorithm (n pop ) and at each generation (n epochs ), a full Monte Carlo simulation is required for the evaluation of the performance of the power grid and the subsequent estimation of the UFII.Therefore, a total of 2 × n mc × n pop × n epochs AC-OPF runs are required.Hence, efficient sampling and caching procedures are proposed to reduce the computational cost of the optimization.The n mc samples are generated only once at the beginning of the optimization procedure, and the corresponding costs C o are computed and saved for future UFII evaluations.In fact, the costs of the initial grid C o are independent of the allocation vector and are only affected by the n mc random scenarios δ t .Hence, by fixing the randomized sample, the C o results can be conveniently re-used to compute the UFII of new PV investments.This procedure reduces half the computational cost.Further, a limited number of scenarios n mc are used to keep the computational cost of the analysis feasible.For more details and a discussion of a convergence analysis study, the reader is reminded of the case study section.

Computational Tool and Software
The proposed algorithm for optimal sizing and locating transmission-connected PV plants is summarized in Figure 1, where each PV plant to be installed in a power grid is a decentralized generator.The computational approach optimizes four objective functions for different levels of robustness (see Equation ( 8)) using an evolutionary algorithm, whereas the grid performance is simulated employing a combination of AC-OPF and Monte Carlo simulation.Matlab is the chosen working platform due to its easy applicability, high functionality, and flexibility, together with ready tools in its library.Genetic algorithm and nonlinear programming optimizations are executed via Matlab's ga and fminbnd builtin functions, respectively.Matpower [45] is used for modeling the power network and for running energy system functions such as OPF.OpenCossan software [46] is used for performing the uncertainty quantification part.The optimizations are conducted on a laptop computer with a 64-bit operating system, 1.8 GHz CPU, and 8 GB RAM.
Monte Carlo simulation.Matlab is the chosen working platform due to its easy applicability, high functionality, and flexibility, together with ready tools in its library.Genetic algorithm and nonlinear programming optimizations are executed via Matlab's ga and fminbnd built-in functions, respectively.Matpower [45] is used for modeling the power network and for running energy system functions such as OPF.OpenCossan software [46] is used for performing the uncertainty quantification part.The optimizations are conducted on a laptop computer with a 64-bit operating system, 1.8 GHz CPU, and 8 GB RAM.

14-bus IEEE Power Grid
The 14-bus IEEE power grid is a simplified version of a part of the American Electric Power System as of 1962 and is a commonly used power grid test case [45].It is a representation of the high-voltage transmission network consisting of five thermal power plants, and its peak total active and reactive power demands are 259 MW and 73.5 MVAr, respectively.Default  and  values are kept between 0.94 and 1.06 p.u., respectively, for  = 1, … ,14.Bus, branch, and generator data for the 14-bus IEEE power grid are listed in Appendix B. To demonstrate the applicability of the approach to a realistic system, the test case is modified by considering four geographical regions characterized by different climate conditions, i.e., cloud coverage.Figure 2 shows the test case and the four regions A, B, C, and D with increasing cloudiness.

Case Study 5.1. 14-bus IEEE Power Grid
The 14-bus IEEE power grid is a simplified version of a part of the American Electric Power System as of 1962 and is a commonly used power grid test case [45].It is a representation of the high-voltage transmission network consisting of five thermal power plants, and its peak total active and reactive power demands are 259 MW and 73.5 MVAr, respectively.Default v min i and v max i values are kept between 0.94 and 1.06 p.u., respectively, for i = 1, . . ., 14. Bus, branch, and generator data for the 14-bus IEEE power grid are listed in Appendix B. To demonstrate the applicability of the approach to a realistic system, the test case is modified by considering four geographical regions characterized by different climate conditions, i.e., cloud coverage.Figure 2 shows the test case and the four regions A, B, C, and D with increasing cloudiness.

Parameter Selection
Power demand: The active and reactive loads for each bus of the 14-bus IEEE power grid case study, as illustrated in Figure 3.These active and reactive load values are taken to be the annual peak values, and the 8760-h power demand curves are normalized ac-

Parameter Selection
Power demand: The active and reactive loads for each bus of the 14-bus IEEE power grid case study, as illustrated in Figure 3.These active and reactive load values are taken to be the annual peak values, and the 8760-h power demand curves are normalized according to these peaks.As mentioned in Section 4.1, the random variability of power demand around the hourly values is modeled as a Gaussian noise with a coefficient of variation of 5%.The variability of the consumer power demand due to the presence and effects of extreme weather conditions are also taken into consideration.

Parameter Selection
Power demand: The active and reactive loads for each bus of the 14-bus IEEE power grid case study, as illustrated in Figure 3.These active and reactive load values are taken to be the annual peak values, and the 8760-h power demand curves are normalized according to these peaks.As mentioned in Section 4.1, the random variability of power demand around the hourly values is modeled as a Gaussian noise with a coefficient of variation of 5%.The variability of the consumer power demand due to the presence and effects of extreme weather conditions are also taken into consideration.

Uncertainty Propagation and Convergence Analysis
Since it is not possible to produce and analyze an infinite number of scenarios, it is important to find the minimal number of Monte Carlo samples needed for an acceptable convergence of the result.A convergence analysis is conducted to analyze the behavior of the mean, standard deviation, and coefficient of variation of the UFII function as a function of the Monte Carlo samples.Twenty randomly generated  vectors are created where 0  100 for all  ∈ {1, … ,  } and 0 ≤  ≤ 150.The statistical indicators are computed using up to  = 10,000 samples.Each realization contains 45 input parameters, summarized in Table 2.  Grid failures and power demand: The expected number of failure events are modeled as a Poisson distribution with parameter λ i , while the duration of the event is modeled as a normal distribution with a mean µ i and standard deviation σ i .The numerical values of the parameters of the distributions are chosen to reflect realistic scenarios and are listed in Table 1.

Uncertainty Propagation and Convergence Analysis
Since it is not possible to produce and analyze an infinite number of scenarios, it is important to find the minimal number of Monte Carlo samples needed for an acceptable convergence of the result.A convergence analysis is conducted to analyze the behavior of the mean, standard deviation, and coefficient of variation of the UFII function as a function of the Monte Carlo samples.Twenty randomly generated PV vectors are created where 0 < PV i < 100 for all i ∈ {1, . . . ,n b } and 0 ≤ PV T ≤ 150.The statistical indicators are computed using up to n mc = 10, 000 samples.Each realization contains 45 input parameters, summarized in Table 2.The applicability of the optimization approach relies on the stability of the ranking of different investment vectors.Due to the stochastic nature of the fitness functions, the fitness hierarchy of the competitive investments could likely change in different epochs of the GA.This behavior is highly undesirable as it can lead to instability problems and suboptimal solutions.This work overcomes this issue by proposing an intelligent propagation scheme that maximizes solution ranking stability and simplifies the comparison of different candidate investments.The fitness of each investment PV is evaluated on an identical set of The samples of the uncertain factors are generated before running the optimization routine and used during each epoch to evaluate the fitness scores.In other words, all proposed solutions shared the same random scenarios.
The results of the convergence analysis are shown in Figure 6.The expected UFII and standard deviation converge around 500 samples, where the coefficient of variation is after 400 samples.However, the convergence graphs display parallel trends in the curves among different investment options (PVs), even at low sample amounts.This means that the ranking of the 20 investments is very stable, thanks to the proposed uncertainty propagation scheme.Therefore, 250 samples are used in the analysis, which represents a good compromise between practicality and accuracy for the estimation of the parameters of interest.

Optimal PV Allocation
The optimization approach based on the UFII is used to identify the optimal and robust allocation of decentralized PV generators.Optimal PV allocations are obtained using four values of the robustness parameter α, and the objective function is estimated using 250 samples.The optimization parameters, design variables, constraints, and linear inequalities are summarized in Table 3. b u is the upper limit for the PV unit capacity on a single node; and b T is the upper limit for the total capacity of all the PV units on all nodes.These parameters are chosen empirically with the aim of an accurate and efficient optimization, where b T is particularly set to prevent the uncontrollable PV generation from exceeding the total active power demand and assure OPF convergence.400 samples.However, the convergence graphs display parallel trends in the curves among investment options (s), even at low sample amounts.This means that the ranking of the 20 investments is very stable, thanks to the proposed uncertainty propagation scheme.Therefore, 250 samples are used in the analysis, which represents a good compromise between practicality and accuracy for the estimation of the parameters of interest.

Optimal PV Allocation
The optimization approach based on the UFII is used to identify the optimal and robust allocation of decentralized PV generators.Optimal PV allocations are obtained using four values of the robustness parameter , and the objective function is estimated using 250 samples.The optimization parameters, design variables, constraints, and linear inequalities are summarized in Table 3.  is the upper limit for the PV unit capacity on a  4. The optimization converged in 51 generations which took between h to complete.Faster results can be obtained using a computer with higher specifications or a computer cluster rather than a laptop computer.The results obtained give valuable information for identifying the critical buses for achieving a higher penetration of PV (and corresponding better performance of the UFII index) and the influence of the robustness parameter on the PV allocation strategy.For example, Buses 1 and 5 are the main driving force behind a high value of E[UFI I]; however, the total PV capacity needs to be distributed over the power grid, with lesser amounts on Buses 1 and 5, for achieving a smaller σ[UFI I].It is worth mentioning that relatively very small PV capacities are also suggested to be installed, e.g., 0.1 MW P PV capacity on Bus 7 according to the PV * obtained by the optimization procedure without consideration of robustness (α = 0).These very small PV installations at the transmission level are noises of the optimization algorithm and can be omitted since their impacts are negligible.
Figure 7 presents the optimal allocations of decentralized PV generators suggested for the four levels of robustness analyzed.The PV allocation strategy follows the characteristics of cloud coverage of the four geographical regions, with the most profitable investments obtained by installing the PVs in the areas with the least amount of cloud.Consider as an example region A which has the lowest amount of cloud coverage and includes Buses 1, 2, and 5.Although Buses 2 and 5 are the only load-carrying nodes in region A, the optimizer found that it is financially efficient to invest in Bus 1 and then forward its production to the remaining part of the network.
The results tend to be distributed over the entire power grid rather than focusing on a few buses to obtain higher robustness.In other words, a more uniformly distributed allocation of PV plants reduces the variance of UFII, thus the risk and uncertainty associated with the financial investment.In fact, as the robustness level, in order to minimize σ[UFI I] the optimization reduces the variability of the UFII values by increasing the lowest scores and reducing the highest scores of the UFII.For comparison purposes, the probability density functions (PDFs) and cumulative distribution functions (CDFs) of the resultant UFII values obtained from PV allocations suggested by α = 0 and α = 3 with 250 samples are plotted in Figure 8. Please note that there are also negative values of UFII obtained by the optimization procedure.This happens because the generator O&M costs (€/h) of PVs result in C pv > C 0 for scenarios with very low PV production, e.g., at night or for very cloudy weather conditions.Therefore, the extra O&M costs added to the system not compensated by the energy produced by the PVs, producing the corresponding negative UFII values.must be complied with.When planning and designing a system, one must pay attention to this trade-off as well as the needs and priorities of the power system operators and/or investors.must be complied with.When planning and designing a system, one must pay attention to this trade-off as well as the needs and priorities of the power system operators and/or investors.

Verification of the Solution
The identified optimal solutions estimated using a small number of Monte Carlo scenarios ( = 250) are verified against the results obtained by estimating the probabilistic performance of the PV investments on a larger set of scenarios.Ten years, corresponding to  = 876,000 samples are generated, and the first and second moments of the UFII are calculated.The estimators and errors (percentages) are reported in Table 5.Note that the estimators computed with 250 samples are reasonably close to the ones obtained using large sample sizes.Moreover, the ranking of the expectations, variance, and lowerand higher-risk solutions also remains unchanged.The greatest difference between the expectations is found to be only 2.4% for  = 1, which is very low considering the large computational gain attained by the low sample size.

Effect of Natural Variability and Randomness
The variability of the UFII scores is estimated via the Monte Carlo method and considers the natural variability of the parameters and randomness effects.The natural variability  [] is due to the variability in time of the nominal values of active/reactive loads and the clear sky irradiances.This can be estimated by generating Monte Carlo samples of time-varying factors but neglecting random effects such as failures, cloud cover-

Verification of the Solution
The identified optimal solutions estimated using a small number of Monte Carlo scenarios (n mc = 250) are verified against the results obtained by estimating the probabilistic performance of the PV investments on a larger set of scenarios.Ten years, corresponding to n mc = 876, 000 samples are generated, and the first and second moments of the UFII are calculated.The estimators and errors (percentages) are reported in Table 5.Note that the estimators computed with 250 samples are reasonably close to the ones obtained using large sample sizes.Moreover, the ranking of the expectations, variance, and lowerand higher-risk solutions also remains unchanged.The greatest difference between the expectations is found to be only 2.4% for α = 1, which is very low considering the large computational gain attained by the low sample size.

Effect of Natural Variability and Randomness
The variability of the UFII scores is estimated via the Monte Carlo method and considers the natural variability of the parameters and randomness effects.The natural variability σ NV [UFI I] is due to the variability in time of the nominal values of active/reactive loads and the clear sky irradiances.This can be estimated by generating Monte Carlo samples of time-varying factors but neglecting random effects such as failures, cloud coverages, and heat waves.The total variability σ[UFI I] is estimated to account for all the sources of variability of the parameters.
The relative contribution of randomness on σ[UFI I] is computed as: I] − σ NV [UFI I])/σ[UFI I] and showing a very high contribution to σ[UFII] from the natural variability of the signal.This contribution is also more significant for the solution obtained using larger values of α, i.e., the randomness contributes to only 0.45 percent of the total standard deviation for α = 3 and 10.4% for α = 0. Hence, the optimization with higher α gives more importance to the minimization of the σ[UFII] and effectively improves the stability (reduced uncertainty) of the instantaneous returns.However, even for high α this optimization procedure cannot reduce the σ NV [UFI I] much further.This is due to the natural and unavoidable variability in the reward process, e.g., higher/lower costs during nightly/daily hours.

Effect of Sampling Uncertainty
The epistemic uncertainty affecting the estimator of E[UFI I] due to the effect of random sampling can be estimated by computing the variance of the estimator: V(E[UFI I]).The variance of the estimator is shown in Table 2.As a figure of merit, an estimate for the coefficient of variation can be computed as CoV =

Comparison of UFII with Other Objective Functions Studied in the Literature
The UFII metric and the optimal decentralized PV generator allocation vector resulting from its maximization are compared with the results of well-established different economic objective functions used in the literature.The first objective maximizes the running cost of the power system whilst the second minimizes the levelized cost of energy (LCOE) [47].For this purpose, two more optimizations are executed with the same probabilistic framework, methodology, and constraints but with the following objective functions: where ∆cost = C o − C pv ; O&M PV is the total O&M cost of the PV systems [€/h]; P PV is the power output of the PV systems [MW].
The PV allocations obtained using the objective functions defined in Equations ( 9) and (10) are compared with the results obtained using the proposed UFII.The summary of the results is shown in Table 6.The results are obtained using n mc = 250 samples and the same optimization parameters listed in Table 6.The results show that the UFII metric provided a slightly reduced expected payback period of the investment and a smaller associated standard deviation (uncertainty).This demonstrates that the novel metric proposed in this work can be used as an effective indicator for economically analyzing a PV investment.The main reason why UFII performs better than LCOE is that LCOE directly measures the unit cost of the electrical energy generated by the PVs, but it does not take into account the changes in the power flows and respective possible power losses, which might create an increase in the running cost of the whole system.On the other hand, maximizing ∆cost is outperformed by maximizing UFII because the former suggests adding the maximum amount of PV allowed without paying attention to the investment's economic efficiency, and the findings of the comparison confirm this behavior.

Conclusions
This work has introduced a novel metric of the unit financial impact indicator, which allows for analyzing the efficacy of financial investments on PVs.A Monte Carlo approach is used to estimate the probability of the UFII score for different investment schemes, allowing unavoidable uncertainties affecting the system and the investment to be quantified.The probabilistic approach has been embedded within an evolutionary optimization algorithm to search for the optimal sizes and positions of PV units.The optimal solutions provide a compromise between expected profit and financial risks.The optimized (risk-profit) investments minimize a weighted sum of the expected UFII and its variance.The 14-bus IEEE power distribution grid tests the proposed method for finding an optimal PV allocation strategy.Based on the numerical results, we claim that the proposed framework can deliver investment solutions that are both robust (minimize the UFII uncertainty/variance) and profitable (maximize the expected long-term profit).The results confirmed that the sizing and siting of decentralized PV generators matter and play a fundamental role in achieving a successful financial investment.The unit financial impact indicator was compared with other commonly used economic metrics.The findings of the two case studies revealed that the UFII performs slightly better in achieving a PV unit positioning and sizing strategy that minimizes the payback period of the investment.The main advantage and novelty of the UFII are that it not only considers the costs and generations associated with PVs but also accounts for the generation, consumption, cost, and power flow parameters of the whole power system.Investors and system operators could consider utilizing the UFII, in combination with other financial performance scores, to support decision-making on how to allocate PV plants.
Because the optimization approach selects candidate solutions only based on the first and second-order moments, a small number of samples is needed to estimate the moment of the UFII distribution.This allows a significant reduction of the computational cost of the analysis yet allows a good treatment of the uncertainty.To further reduce the computational cost, an efficient sampling strategy was developed and used to improve the stability of the optimizer.These features make the proposed approach generally applicable to power grids of different sizes.The proposed method is generally applicable and extendable to include more technical parameters for the sake of achieving a resilient power grid.The results of the optimization prove that the proposed framework can yield investment solutions that are both robust (minimize the UFII uncertainty/variance) and profitable (maximize the expected long-term profit).

Future Recommendations
While the methodology initially aimed at evaluating the efficacy of large-scale PV investments considering unavoidable uncertainty and grid stability, it may also be utilized for evaluating individual rooftop photovoltaic systems.The innovative metric proposed for evaluating the impact of financial investments on photovoltaic systems can additionally be employed to analyze the optimal distribution of energy storage systems and assess the effectiveness of carbon taxes in stimulating the transition towards more sustainable energy sources.
2. Moreover, epistemic uncertainty may affect the sample-based approximation of E[UFI I] due to the limited number of samples used, and this effect is estimated by quantifying the variance of the estimator: V(E[UFI I]).The V(E[UFI I]) quantifies the variability of the estimator of E[UFI I] and differs from the standard deviation σ[UFI I], which is the second moment of the stochastic process.Therefore V(E[UFI I]) measures the quality of the estimator while σ[UFI I] relates to the inherent variability of instantaneous rewards, i.e., the variability of the σ[UFI I] under different scenarios.A low σ[UFI I] indicates that instantaneous gains are stable, independently of the operational and weather scenario, whilst a low V(E[UFI I]) indicates high confidence in the estimated payback period.

Figure 1 .
Figure 1.Flowchart representation of the proposed methodology.

Figure 1 .
Figure 1.Flowchart representation of the proposed methodology.

Figure 3 .
Figure 3. Annual peak active and reactive load values of 14 buses.

Figure 4
Figure 4 shows an example of the hourly active demand curves for Bus 3 resulting from the proposed load-generating method for four different days of the year.The original active load value of 94.2 MW for Bus 3 is taken as the annual peak value.Please note that none of these graphs include the effect of heat waves at this point.Additional graphs illustrating the irradiance, load, presence of heat waves, and component failures over a time horizon resulting from the Monte Carlo simulator are presented in Appendix C.

Figure 3 .
Figure 3. Annual peak active and reactive load values of 14 buses.

Figure 4
Figure 4 shows an example of the hourly active demand curves for Bus 3 resulting from the proposed load-generating method for four different days of the year.The original active load value of 94.2 MW for Bus 3 is taken as the annual peak value.Please note that none of these graphs include the effect of heat waves at this point.Additional graphs illustrating the irradiance, load, presence of heat waves, and component failures over a time horizon resulting from the Monte Carlo simulator are presented in Appendix C. Sustainability 2023, 15, 11715 11 of 25

Figure 4 .
Figure 4. Hourly active load curves on Bus 3 for the 14-bus IEEE power grid on the 55th, 144th, 216th, and 303rd days of the year and assuming an annual peak value of 94.2 MW.Heat waves: The number of days with a heat wave  are modeled with a Poisson distribution with an occurrence rate  = 5 over July and August, and  = 0 otherwise.Solar irradiance: The clear sky model prepared by the Solar Energy Research Institute for the U.S. Department of Energy is used to model solar irradiance over time [41].The input parameters are chosen for Cyprus, and the resultant dataset gives the direct normal irradiation (), zenith angle (θ), direct horizontal irradiation (), global horizontal irradiation (), and diffuse horizontal irradiation () with hourly resolution without considering the effect of clouds.Cloud coverage: Hourly irradiance samples, , are calculated for each of the geo-

Figure 4 .
Figure 4. Hourly active load curves on Bus 3 for the 14-bus IEEE power grid on the 55th, 144th, 216th, and 303rd days of the year and assuming an annual peak value of 94.2 MW.Heat waves: The number of days with a heat wave E e are modeled with a Poisson distribution with an occurrence rate λ E = 5 events month over July and August, and λ E = 0 otherwise.Solar irradiance: The clear sky model prepared by the Solar Energy Research Institute for the U.S. Department of Energy is used to model solar irradiance over time[41].The input parameters are chosen for Cyprus, and the resultant dataset gives the direct normal

Figure 5 25 Figure 5 .
Figure 5.The graph shows the mean value of the irradiance (in red) plus and minus one standard deviation (in gray) of two winter days with historically little clouds (leftmost curves) and two cloudy days in summer (two profiles on the right).The irradiance profile from the clear sky is shown in yellow.

Figure 5 .
Figure 5.The graph shows the mean value of the irradiance (in red) plus and minus one standard deviation (in gray) of two winter days with historically little clouds (leftmost curves) and two cloudy days in summer (two profiles on the right).The irradiance profile from the clear sky is shown in yellow.

Figure 6 .
Figure 6.Convergence analysis for the estimation of (a) the mean, (b) standard deviation, and (c) the coefficient of variation of UFII using 20 different investment options.

Figure 6 .
Figure 6.Convergence analysis for the estimation of (a) the mean, (b) standard deviation, and (c) the coefficient of variation of UFII using 20 different investment options.The sizes and locations of decentralized PV generators obtained for the four values of robustness selected and their corresponding E[UFI I] and σ[UFI I] values are listed in

Figure 8 .
Figure 8. Distribution of the UFII values obtained with  = 0 and  = 3: probability density functions (left panel) and the cumulative distribution functions (right panel).

Figure 9
Figure 9 shows the Pareto front of the obtained values of [] and [] corresponding to different values of robustness.It is obvious that there is compensation between [] and [] , i.e., to reduce [] , thus uncertainty, a lower []must be complied with.When planning and designing a system, one must pay attention to this trade-off as well as the needs and priorities of the power system operators and/or investors.

Figure 8 .
Figure 8. Distribution of the UFII values obtained with  = 0 and  = 3: probability density functions (left panel) and the cumulative distribution functions (right panel).

Figure 9
Figure 9 shows the Pareto front of the obtained values of [] and [] corresponding to different values of robustness.It is obvious that there is compensation between [] and [] , i.e., to reduce [] , thus uncertainty, a lower []must be complied with.When planning and designing a system, one must pay attention to this trade-off as well as the needs and priorities of the power system operators and/or investors.

Figure 8 .
Figure 8. Distribution of the UFII values obtained with α = 0 and α = 3: probability density functions (left panel) and the cumulative distribution functions (right panel).

Figure 9
Figure 9 shows the Pareto front of the obtained values of E[UFI I] and σ[UFI I] corresponding to different values of robustness.It is obvious that there is compensation between E[UFI I] and σ[UFI I], i.e., to reduce σ[UFI I], thus uncertainty, a lower E[UFI I] must be complied with.When planning and designing a system, one must pay attention to this trade-off as well as the needs and priorities of the power system operators and/or investors.

Figure 9 .
Figure 9. Obtained values of [] vs. for the different levels of robustness.

Figure 9 .
Figure 9. Obtained values of E[UFI I] vs. σ[UFI I] for the different levels of robustness.

Figure A2 .
Figure A2.Sampled irradiance on region A.

Figure A2 .
Figure A2.Sampled irradiance on region A.

Figure A2 .
Figure A2.Sampled irradiance on region A.
PV cell under standard test conditions [540 W]; G STC is the incident irradiation at standard test conditions [1000 W/m 2 ]; ∝ P is the temperature coefficient of P max [−0.341 %/ • C]; T C,STC is the PV cell temperature under standard test conditions [

Table 1 .
Line and thermal power plant outage model.Expected number of outages per year (. ), mean ( ) and standard deviation ( ) of the duration of the outage.

Table 2 .
Data set of each scenario analysis.

Table 1 .
Line and thermal power plant outage model.Expected number of outages per year (.λ i ), mean (τ i ) and standard deviation (σ i ) of the duration of the outage.

Table 2 .
Data set of each scenario analysis.

Table 3 .
Optimization parameter used for the robust PV allocation.

Table 4 .
PV allocation in MWp for different levels of robustness α.

Table 5 .
Verification of the optimal solution.

Table 5 .
Verification of the optimal solution.

Table 6 .
PV allocation in MWp for different levels of robustness α.