Genetic Algorithms-Based Optimum PV Site Selection Minimizing Visual Disturbance

: In this paper, an integrated methodology is developed to determine optimum areas for Photovoltaic (PV) installations that minimize the relevant visual disturbance and satisfy spatial constraints associated with land use, as well as environmental and techno-economic siting factors. The visual disturbance due to PV installations is quantiﬁed by introducing and calculating the “Social Disturbance” (SDIS) indicator, whereas optimum locations are determined for predeﬁned values of two siting preferences (maximum allowable PV locations—grid station distance and minimum allowable total coverage area of PV installations). Thematic maps of appropriate selected exclusion criteria are produced, followed by a cumulative weighted viewshed analysis, where the SDIS indicator is calculated. Optimum solutions are then determined by developing and employing a Genetic Algorithms (GAs) optimization process. The methodology is applied for the municipality of La Palma Del Condado in Spain for 100 different combinations of the two siting preferences. The optimization results are also employed to create a ﬂexible and easy-to-use web-GIS application, facilitating policy-makers to choose the set of solutions that better fulﬁls their preferences. The GAs algorithm offers the ability to determine distinguishable, but compact, regions of optimum locations in the region, whereas the results indicate the strong dependence of the optimum areas upon the two siting preferences.


Introduction
The solar Photovoltaic (PV) market presents one of the most dynamic renewable energies markets, and its growth is expected to rise [1]. Due to their environmental and economic benefits, PV plants provide energy for numerous applications, in addition to the electrical grid supply, such as solar irrigation [2], seawater desalination [3], groundwater pumping [4], as well as the combined production of agricultural crops and power [5]. The technology of PVs is continuously being developed to improve solar cells efficiency (e.g., [6,7]). Regarding the deployment of PVs in residential areas, many researchers have developed and applied various relevant tools and methods, including, for example, tools that provide access to PV datasets [8] and advanced integrated approaches for rooftop solar energy assessments (e.g., [9][10][11]). Furthermore, there is a number of studies focusing on PV plants' allocation in rural areas (e.g., [12,13]).
When considering the site selection criteria for installing PV plants, it is necessary to take into consideration the associated potential risks (however small they might be) [14,15], while at the same time setting the regulations that will not prevent the PV plants' further expansion. Authorities often regulate the PV plants' deployment based on environmental and social criteria [16]. As an inaccurate sitting of PV plants might lead to public oppositions, it is necessary, when proposing regulations, to take into account all those parameters that will make citizens feel safe and, thus, will increase public acceptance [17][18][19]. Specifically, bance, while satisfying, at the same time, spatial constraints associated with land use, and environmental and techno-economic siting factors. The visual disturbance of the public due to PV installations is quantified by introducing and calculating an indicator called herein the "Social Disturbance" (SDIS) indicator. The largest the SDIS indicator for an area is, the higher the visibility of this area will be, advocating higher visual disturbance of the public (i.e., social disturbance) due to PV installations in this area. The proposed methodology is applied for the municipality of La Palma Del Condado located in in the province of Huelva, Andalucia, Spain. Initially, thematic maps of exclusion siting criteria with their incompatibility zones and a Digital Surface Model (DSM) of the region are developed using GIS. A cumulative weighted viewshed analysis follows to calculate the SDIS indicator. Successively, the optimization problem is solved by developing a GAs-driven optimization process, where optimum locations are determined for predefined values of the following two siting preferences: (a) the maximum allowable PV locations-grid station distance; and (b) the minimum allowable total coverage area of PV installations. By modifying the values of the aforementioned preferences, different sets of optimum solutions are obtained. These solutions are, finally, employed as a basis to create a flexible and easy-to-use web-GIS application, facilitating policy-makers to choose the set of solutions that better fulfils their preferences/strategies. The mitigation of social oppositions resulting from visual disturbances of PV installations presents an important factor that can contribute to the sustainable site selection of these renewable energy systems. The present paper tackles this problem and fills existing research gaps by: (a) explicitly formulating the PV site selection problem as an optimization problem, where areas that minimize the PVs' visual disturbance and satisfy spatial constraints are being sought; and (b) proposing an integrated methodology to solve the aforementioned optimization problem, where GIS databases, viewshed analysis, and GAs are efficiently combined.
The remainder of the paper is organized as follows: in Section 2, the study area is presented. Section 3 includes a detailed description of the components of the developed methodology. The results of the present paper are presented and discussed in Section 4, and, finally, in Section 5, the concluding remarks and key findings of this investigation are cited.

Study Area
The region under investigation corresponds to the municipality of La Palma Del Condado located in Andalucia, Spain, in the province of Huelva (Figure 1a). The municipality has (2021) a population of approximately 10,700 people [45] and covers an area of 60.5 km 2 . The study area (Figure 1b) consists of agricultural lands, including vineyards that are critical for the local economy, residences, and urban areas, main roads, and railways. Furthermore, it is surrounded by Natura 2000 areas. Specifically, the Rio Tinto River (Figure 1b), as well as the wider Natura 2000 zone in the northern part of the area, are famous attractions to tourists and scientists because of the special characteristics of their aquatic environment.
High temperatures and solar radiation values make the study area favorable in terms of PV productivity [46]. It should be mentioned, however, that, currently, a specific spatial planning framework for deploying PVs in the region does not exist and the municipality adopts general siting criteria applicable to any type of structure (e.g., buildings) for new PV installations. Along these lines, the municipality targets to develop a relevant framework that will make it possible for investors to install PVs, while, at the same time, respecting the physical and the socio-economic environment.

Methodology
The developed methodology in this paper aims at identifying optimum areas for PV installations on municipality scale that minimize the relevant visual disturbance, quantified by the SDIS indicator, and satisfy spatial constraints related to land use, environmental, as well as to techno-economic siting factors. By minimizing the SDIS indicator, mitigation of potential social oppositions and negative impacts on land activities (e.g., tourism, residency construction) can be achieved.
The whole methodology is shown schematically in Figure 2. Initially, a GIS database is created, including thematic maps of exclusion criteria with their incompatibility zones (spatial constraints), and a DSM of the region. The latter model is successively deployed in order to perform a cumulative weighted viewshed analysis, where viewshed maps are created and the SDIS indicator is calculated and mapped. The spatial constraints and the SDIS maps are, then, used to perform a GAs-driven optimization process, where optimum locations are determined for predefined values of two siting preferences corresponding to: (a) the maximum allowable PV locations-grid station distance; and (b) the minimum allowable total coverage area of PV installations. By modifying the values of the aforementioned siting preferences, different sets of optimum solutions are obtained, which are, finally, used as a basis to create a web-GIS application. The latter tool enables the visualization of PV plants' optimum locations in the study area for different bounds of the PV locations-grid station in-between distance and of the PV locations' total coverage area, facilitating policy-makers to choose a set of solutions that fulfils better their preferences on the above factors. In the following sections, the components of the present methodology are described in detail.

Methodology
The developed methodology in this paper aims at identifying optimum areas for PV installations on municipality scale that minimize the relevant visual disturbance, quantified by the SDIS indicator, and satisfy spatial constraints related to land use, environmental, as well as to techno-economic siting factors. By minimizing the SDIS indicator, mitigation of potential social oppositions and negative impacts on land activities (e.g., tourism, residency construction) can be achieved.
The whole methodology is shown schematically in Figure 2. Initially, a GIS database is created, including thematic maps of exclusion criteria with their incompatibility zones (spatial constraints), and a DSM of the region. The latter model is successively deployed in order to perform a cumulative weighted viewshed analysis, where viewshed maps are created and the SDIS indicator is calculated and mapped. The spatial constraints and the SDIS maps are, then, used to perform a GAs-driven optimization process, where optimum locations are determined for predefined values of two siting preferences corresponding to: (a) the maximum allowable PV locations-grid station distance; and (b) the minimum allowable total coverage area of PV installations. By modifying the values of the aforementioned siting preferences, different sets of optimum solutions are obtained, which are, finally, used as a basis to create a web-GIS application. The latter tool enables the visualization of PV plants' optimum locations in the study area for different bounds of the PV locations-grid station in-between distance and of the PV locations' total coverage area, facilitating policy-makers to choose a set of solutions that fulfils better their preferences on the above factors. In the following sections, the components of the present methodology are described in detail.

Methodology
The developed methodology in this paper aims at identifying optimum areas for PV installations on municipality scale that minimize the relevant visual disturbance, quantified by the SDIS indicator, and satisfy spatial constraints related to land use, environmental, as well as to techno-economic siting factors. By minimizing the SDIS indicator, mitigation of potential social oppositions and negative impacts on land activities (e.g., tourism, residency construction) can be achieved.
The whole methodology is shown schematically in Figure 2. Initially, a GIS database is created, including thematic maps of exclusion criteria with their incompatibility zones (spatial constraints), and a DSM of the region. The latter model is successively deployed in order to perform a cumulative weighted viewshed analysis, where viewshed maps are created and the SDIS indicator is calculated and mapped. The spatial constraints and the SDIS maps are, then, used to perform a GAs-driven optimization process, where optimum locations are determined for predefined values of two siting preferences corresponding to: (a) the maximum allowable PV locations-grid station distance; and (b) the minimum allowable total coverage area of PV installations. By modifying the values of the aforementioned siting preferences, different sets of optimum solutions are obtained, which are, finally, used as a basis to create a web-GIS application. The latter tool enables the visualization of PV plants' optimum locations in the study area for different bounds of the PV locations-grid station in-between distance and of the PV locations' total coverage area, facilitating policy-makers to choose a set of solutions that fulfils better their preferences on the above factors. In the following sections, the components of the present methodology are described in detail.

Development of the GIS Database
In order to solve the examined siting optimization problem, the development of an adequate GIS database is initially required. This database includes: (a) thematic maps of exclusion criteria (with their incompatibility zones) that impose spatial limitations for the PV installations in the study area based on utilization restrictions, and environmental and techno-economic considerations; and (b) a DSM that is deployed in the preceding viewshed analysis. In the present analysis, five (5) exclusion criteria have been taken into account (Table 1). EC1 contributes to mitigate the potential environmental risk associated with PV installations. By applying this criterion, all Natura 2000 zones of recognized natural and ecological value (e.g., [15]) are not considered suitable for the deployment of PVs. The same holds true for the areas within a distance of 150 m from the aforementioned zones. This distance has been selected based on relevant limits available in the literature (e.g., 500 m in [49,50], 100 m in [51]), taking into account the size of the examined study area.
EC2 and EC3 are considered for safety and technical reasons. EC2 excludes the existing major roads and all areas within a distance of 100 m from them, which has been defined according to the relevant limits found in the literature (e.g., [49,52]). In a similar manner, EC3 excludes the existing railway network and a 50 m buffer safety zone from it. The latter limit has been selected according to [53].
Regarding EC4, all non-agricultural areas (e.g., urban areas, residences etc.) are excluded from the optimization process. The same holds true for vineyards, since, according to the municipality authorities, they are considered extremely important for the local economy.
As for EC5, the deployment of PVs in areas with large terrain slopes leads to high investment costs. Accordingly, an upper bound of the terrain slope is taken into account for technical-economic reasons. A wide range of the aforementioned bound can be found in the literature: 2% in [46], 15% in [15], and 25% in [51,54]. In the present paper, the value of 5% has been selected.
All spatial data related to the exclusion criteria were collected from various sources (Table 1), and they were appropriately processed to develop the GIS database and the relevant thematic maps using QGIS. In the case of EC5, the gradient field consideration requires an adequate Digital Terrain Model (DTM). For this purpose, Lidar data have been used and processed using R package, "lidR" [55]. The developed DTM model was inserted in QGIS, and the thematic map of EC5 (raster format) was created using the "slope" function. The high accuracy of the produced map (the Lidar data accuracy was 1.5 points/m 2 ) can result in standalone points with very large slope terrain values (representing, for example, a tree). To reduce the effect of these outliers, the average slope of the polygons examined was calculated. If this average value was smaller than 5%, the polygon was considered as a potential area for PV installations in terms of EC5. It is noted that Lidar data were also used to create the DSM model required in the viewshed analysis, including both the ground elevation and the objects' elevations (buildings, trees etc.). Finally, it should be mentioned that the exclusion criteria of Table 1 have been selected by taking into account the special features of the examined region, although numerous other criteria can be found in the relevant literature. For example, historical heritages or archeological sites have not been considered, since they do not exist in the study area. The altitude has also been omitted, since the maximum altitude of the examined region (<200 m) is far away from the upper bound of 1500-2000 m found in the literature [15]. In a similar manner, solar radiation (e.g., [46]) and temperature (e.g., [56]) have not been considered as exclusion criteria, since the limited area of the municipality examined leads to negligible variations of the two aforementioned quantities.

Vieswhed Analysis
In the present investigation, a cumulative weighted viewshed analysis is performed by aggregating the visibility from multiple observers, while taking into account the effect of the distance between the observers and the locations of the PV installations. The overall aim of this analysis, which is realized using the "QGIS visibility analysis" plugin [57], is the quantification and the mapping of the SDIS indicator in the study area.
Starting with the observers, two categories of observers have been taken into account corresponding to residents and drivers. Thus, the spatial placement of the observers has been performed according to the residences' (buildings) location and the roads' traffic, respectively. The residences have been represented as spatial polygons (based on the relevant data obtained from the municipality), and 10% of them were randomly chosen by applying a uniform distribution function to the polygons' IDs. Successively, each of the selected polygons has been replaced by a point in QGIS representing one observer. Regarding the roads, these spatial entities have been divided into major and minor ones based on their traffic characteristics. In the case of the major roads, 10% of the daily traffic was chosen for reasons of conciseness. This percentage was then divided over the total length of each road, leading to 1 observer per 50 m. Due to the absence of data for the minor roads, a proportion of 1/5 was chosen between the major and the minor roads, resulting in 1 observer per 250 m in the latter roads' category. The railroad in the area was treated similarly to the major roads (i.e., 1 observer per 50 m). Furthermore, some extra observers not belonging to the aforementioned categories were placed in the Natura 2000 region, since this zone and especially the Rio Tinto River attract many visitors due their great environmental importance. The location of the observers in a region around the La Palma Del Condado city is shown in Figure 3.
Regarding the distance between the observers and the locations of the PV installations, the closer the distance is, the higher the disturbance perceived by the observers will be. The concept presented is similar to "utility functions" that are often used in economics [58]. Let us consider M different classes of viewshed maps, with an jth, j = 1, . . . , M, class corresponding to a predefined maximum visibility distance, r j . By weighting and aggregating the aforementioned viewshed maps, the disturbance, DIS i,pix , of an ith, i = 1, . . . , N, observer, relevant to a point/pixel, pix, can be quantified according to Equation (1):  Regarding the distance between the observers and the locations of the PV installations, the closer the distance is, the higher the disturbance perceived by the observers will be. The concept presented is similar to "utility functions" that are often used in economics [58]. Let us consider different classes of viewshed maps, with an th, = 1, … , , class corresponding to a predefined maximum visibility distance, . By weighting and aggregating the aforementioned viewshed maps, the disturbance, , , of an th, = 1, … , , observer, relevant to a point/pixel, , can be quantified according to Equation (1): In the above equation, , = 1, … , , corresponds to the weight of the th viewshed class defined below, whereas , , = 1, … , , = 1, … , , is equal to 1 if the examined point/pixel for the th viewshed class can be seen by the th observer, or equal to 0 if the opposite holds true.
By summing up the values of , for all the observers, the SDIS indicator for each point/pixel, , is obtained (see Equation (2) below). for all points/pixels in the study area forms the SDIS map.
The question that arises is how the , = 1, … , , values are determined. To answer this question, it is necessary to take into consideration two factors. The first one is that the disturbance of an observer decreases in a non-linear way with respect to the distance. The non-linear relationship between an object and an observer has been already indicated in the literature [31,59,60]. The second factor is that the disturbances of all the observers will be aggregated; therefore, the principle of additivity needs to be fulfilled, in a way that the By summing up the values of DIS i,pix for all the N observers, the SDIS indicator for each point/pixel, SDIS pix , is obtained (see Equation (2) below). SDIS pix for all points/pixels in the study area forms the SDIS map.
The question that arises is how the c j , j = 1, . . . ,M, values are determined. To answer this question, it is necessary to take into consideration two factors. The first one is that the disturbance of an observer decreases in a non-linear way with respect to the distance. The non-linear relationship between an object and an observer has been already indicated in the literature [31,59,60]. The second factor is that the disturbances of all the observers will be aggregated; therefore, the principle of additivity needs to be fulfilled, in a way that the final SDIS map will offer trustworthy results. To fulfil these two requirements, Equation (3) is introduced to calculate c j , j = 1, . . . , M: where r M is the maximum distance that an observer can see in the case of the M class viewshed (largest maximum distance among all classes), whereas the value of "+1" is used in order to consider the M class as the reference class. For applying Equation (3), the parameter y needs to be defined, having in mind that the smaller the y value is, the more emphasis will be given on objects cited close to an observer. In the present paper, y has been taken as equal to e and, thus, c j , j = 1, . . . , M, is finally calculated suing Equation (4) as follows: In total, 10 classes of viewshed maps have been considered with r j and c j , j = 1, . . . , 10, values shown in Table 2. It can be seen that the increase of r j leads to a nonlinear decrease of c j . The maximum distance of the 10th class, r 10 , is of high importance, since an area located at distance from an observer larger than r 10 does not disturb the observer. The intermediate classes support a better discretization of the distances, and, furthermore, they introduce some extra importance to the smaller distances. The reader should have in mind that as the weighted viewshed maps are aggregated (Equation (1)), for a PV placed between 0 m and 100 m from an observer, the DIS i,pix will be equal to 28.859 and not equal to 4.912. Figure 4 includes two cumulative viewshed maps for r 1 = 100 m and r 6 = 1000 m, assuming c 1 = c 6 =1 in a region around the La Palma Del Condado city. Cumulative weighted viewshed maps, along with the final DSIS map, are cited and discussed in the Results section. viewshed (largest maximum distance among all classes), whereas the value of "+1" is used in order to consider the class as the reference class. For applying Equation (3), the parameter needs to be defined, having in mind that the smaller the value is, the more emphasis will be given on objects cited close to an observer. In the present paper, has been taken as equal to and, thus, , = 1, … , , is finally calculated suing Equation (4) as follows: In total, 10 classes of viewshed maps have been considered with and , = 1, … , 10, values shown in Table 2. It can be seen that the increase of leads to a nonlinear decrease of . The maximum distance of the 10th class, , is of high importance, since an area located at distance from an observer larger than does not disturb the observer. The intermediate classes support a better discretization of the distances, and, furthermore, they introduce some extra importance to the smaller distances. The reader should have in mind that as the weighted viewshed maps are aggregated (Equation (1)), for a PV placed between 0 m and 100 m from an observer, the , will be equal to 28.859 and not equal to 4.912. Figure 4

Optimization Process
The present GAs-driven optimization process aims at identifying optimum areas for PV installations on municipality scale that minimize the SDIS indicator and satisfy spatial constraints introduced with the exclusion criteria of Table 1. Optimum locations are determined for predefined values of two siting preferences: (a) the maximum allowable PV locations-grid station distance; and (b) the minimum allowable total coverage area of PV installations.
In order to form the examined optimization problem, SP is defined as a spatial entity corresponding to a spatially projected vector (polygon or spatial polygon data frame). The overall region, OR, examined (municipality) is divided into K spatial entities; namely, OR = {SP 1 , . . . , SP k , . . . , SP K }, where |K| is the cardinality of the OR superset. Each SP k , k = 1, . . . , K, has a set of spatial and non-spatial attributes, AT k , defined as: In Equation (5), AT k,use includes the land use information of SP k related to EC4 (Table 1), AT k,s is the average slope terrain gradient (%) of SP k related to EC5 (Table 1), AT k,SDIS corresponds to the SDIS of SP k , AT k,area is the coverage area (m 2 ) of SP k , whereas AT k,DG is the distance (m) between the electrical grid station and the most distant point of SP k . Furthermore, AT k,set represents the set (category) where an SP k belongs, and can take the labelled values of "A", "B" or "C" according to Table 3. Set A includes the optimum solutions; namely, all SPs that minimize the DSIS indicator, satisfy none of the exclusion criteria, and are consistent with the two siting preferences. Set C includes SPs that satisfy any of the EC1, EC2, or EC3 exclusion criteria, while, at the same time, correspond to pure non-agricultural areas (EC4). Finally, in set B, SPs that do not belong to either A or B sets are included. Specifically, set B contains SPs that satisfy only EC5 and/or EC4 (i.e., they correspond to vineyards). Accordingly, in this set, areas that are not eligible for PV installations only due to economic reasons are classified. Furthermore, it contains SPs that although do not satisfy any of the exclusion criteria, they do not correspond to solutions that lead to minimum DSIS values and/or are consistent with the two siting preferences. Based on all the above, the GAs optimization process developed herein seeks to determine the SPs that belong to set A by matching and classifying each SP k , k = 1, . . . , K, to one of the A, B, or C sets. This problem has many similarities with graph theory and matching problems [61], as the SPs will be matched to a set according to their attributes. The aforementioned matching and classification are implemented in two successive stages. In the first stage, SPs that satisfy EC1-EC4 are determined by deploying the relevant data and the thematic maps of the GIS database, and by performing the required spatial intersections using a relevant R package. If, for example, an SP k intersects with the 150 m buffer of an environmentally protected area (EC1), then AT k,set = "C". Similarly, for EC4, if AT k,use = "agricultural", then AT k,set = "C". The intersections are performed in each GAs iteration following an approach that will be presented later in this section. The aforementioned SPs are classified to the C set and are subtracted from the OR superset. Each SP k / ∈ C is taken into account in the second stage, where the GAs algorithm determines the elements (SPs) of the A set according to Table 3. Thus, the objective function of the present optimization problem is defined as follows: where Q = A ∪ B. Equation (6) implies that the GAs algorithm seeks for optimum SPs in terms of minimizing the sum of their SDIS indicators. Furthermore, the solution is subjected to the following constraints: AT k,use = "agricultural" & AT k,use = "vineyard", for any SP k ∈ Q AT k,s ≤ 5%, for any SP k ∈ Q (8) In Equation (9), DG max denotes the maximum allowable distance from the grid station for installing PVs, whereas in Equation (10), Area min corresponds to the minimum allowable total coverage area of PV installations. By setting different values for DG min and Area min based on relevant preferences, different optimum solutions can be obtained. Having solved the optimization problem described by Equations (6)-(10) and, thus, having determined the SPs that belong to the A set, the remaining SPs are classified to set B.
In order to label the AT k,use attribute (Equation (7)), the land use data included in the GIS database are used. Similarly, by deploying the slope terrain raster map of the GIS database, the AT k,s attribute (Equation (8)) is quantified and is taken as equal to the average value of the raster pixels that fall inside an SP. As for the AT k,SDIS attribute (Equation (6)), the final SDIS raster produced from the viewshed analysis is "clipped" using the spatial features of an SP as a "mask". The masked SDIS pix (Equation (2)) values are then summed up to quantify AT k,SDIS .
In order to solve the present optimization problem, GAs are deployed, which are numerically realized using the R package, "GA" [62]. GAs are based on the iterative generation of populations of chromosomes, representing possible optimum solutions based on the "survival of the fittest" rule [38][39][40]. Accordingly, the generation of populations of possible SPs with AT k,set = "A" could have been considered. This approach, however, would result to extensively sparse optimum solutions within the examined region, which, in turn, would be difficult to be realized in terms of regional planning. For this reason, another approach is deployed in the present paper, facilitating the formation of compact regions with SP k ∈ A. Specifically, a chromosome is represented by two distinctive triangles, which intersect with and include a number of SPs ( Figure 5). At each iteration, the GAs algorithm initially finds the SPs of the triangles that belong to the C set and excludes them from further analysis. For the remaining SPs of the two triangles, Equations (6)-(10) are then applied, and the SPs satisfying those equations (i.e., SPs∈A) correspond to a potential optimum solution.
A X − Y grid of 10 m is deployed in the area that has the dimensions, 11.227 km × 14.593 km. Considering that 2 10 = 1024 and 2 11 = 2048, any of the X, Y coordinates require 11 bits to be represented in a binary form. Accordingly, 22 bits are required to create a point, 66 to create one triangle, and 132 to create two triangles. The latter 132 bits correspond to a chromosome (i.e., a possible solution to the optimization problem containing SPs that belong to the A set). Three-hundred generations of a population of 125 chromosomes each have been used. The mutation and the crossover probability have been set to 35% and 85%, respectively, determining which of the chromosomes will survive to the next generation. Elitism has been also deployed, letting the best five chromosomes survive always to the next generation. Penalties with successively decreasing values are assigned in the algorithm in a sequential manner as follows: (i) for points or triangle areas that fall outside of the examined region, as well as for triangles that intersect with each other (penalties with the largest values); (ii) for limitations resulting from the exclusion criteria of Table 1 related to set C (penalties with intermediate values); and, finally, (iii) for limitations resulting from the constraints described by Equations (7)-(10) (penalties with the smaller values). In this way, the algorithm efficiently promotes the chromosomes that fulfill as many requirements as possible.

regions with
∈ . Specifically, a chromosome is represented by two distinctive triangles, which intersect with and include a number of ( Figure 5). At each iteration, the GAs algorithm initially finds the of the triangles that belong to the set and excludes them from further analysis. For the remaining of the two triangles, Equations (6)-(10) are then applied, and the satisfying those equations (i.e., ∈ ) correspond to a potential optimum solution. Figure 5. One of the two triangles of a chromosome, intersecting with and including a number of (red blocks). , , = 1, 2, 3, are the spatial coordinates of the triangle's vertices in the deployed − grid. Red lines correspond to minor roads, the pink line represents the railway, and the background is the SDIS raster.
A − grid of 10 m is deployed in the area that has the dimensions, 11.227 km x 14.593 km. Considering that 2 10 = 1024 and 2 11 = 2048, any of the , coordinates require 11 bits to be represented in a binary form. Accordingly, 22 bits are required to create a point, 66 to create one triangle, and 132 to create two triangles. The latter 132 bits correspond to a chromosome (i.e., a possible solution to the optimization problem containing that belong to the set). Three-hundred generations of a population of 125 chromosomes each have been used. The mutation and the crossover probability have been set to 35% and 85%, respectively, determining which of the chromosomes will survive to the next generation. Elitism has been also deployed, letting the best five chromosomes survive always to the next generation. Penalties with successively decreasing values are assigned in the algorithm in a sequential manner as follows: (i) for points or triangle areas that fall outside of the examined region, as well as for triangles that intersect with each other (penalties with the largest values); (ii) for limitations resulting from the exclusion criteria of Table 1 related to set  (penalties with intermediate values); and, finally, (iii) for limitations resulting from the constraints described by Equations (7)-(10) (penalties with the smaller values). In this way, the algorithm efficiently promotes the chromosomes that fulfill as many requirements as possible.
In the present paper, optimum solutions have been found for 100 different combinations of and values. More specifically, 10 different values have been taken into account, varying from 2.5 km to 7.0 km, with a step of 0.5 km. Regarding , 10 different values of this quantity have been also considered, from 0.5 km 2 (≈0.8% of the overall region) to 5.0 km 2 (8% of the overall region), with a step of 0.5 km 2 . To perform all the required computations, the computer cluster, Aristotelis of Aristotle Figure 5. One of the two triangles of a chromosome, intersecting with and including a number of SPs (red blocks). X i , Y i , i = 1,2,3, are the spatial coordinates of the triangle's vertices in the deployed X − Y grid. Red lines correspond to minor roads, the pink line represents the railway, and the background is the SDIS raster.
In the present paper, optimum solutions have been found for 100 different combinations of DG max and Area min values. More specifically, 10 different DG max values have been taken into account, varying from 2.5 km to 7.0 km, with a step of 0.5 km. Regarding Area min , 10 different values of this quantity have been also considered, from 0.5 km 2 (≈0.8% of the overall region) to 5.0 km 2 (8% of the overall region), with a step of 0.5 km 2 . To perform all the required computations, the computer cluster, Aristotelis of Aristotle University of Thessaloniki, has been used, deploying 125 CPUs simultaneously, one for every chromosome of each generation. The corresponding results are, finally, used as a basis to create a relevant web-GIS application, which is presented in the next section. This application is realized by using the R packages, "shiny" [63], "leaflet" [64], "tmap" [65], and "plotly" [66].

Results and Discussion
In this section, the results of the present analysis are presented and discussed. Starting with the thematic maps of the exclusion criteria, Figure 6a shows the incompatibility zones of EC1-EC4 (Table 1), whereas in Figure 6b, the slope terrain thematic map (EC5, Table 1) is cited. Areas that fall within the incompatibility zones of Figure 6a cover 26.9 km 2 and belong to set C. Most of these areas correspond to environmentally protected areas, followed by minor roads and non-agricultural areas.
With regard to the viewshed analysis, Figure 7a shows a cumulative weighted viewshed map indicatively for a maximum observation distance of 5000 m (10th viewshed class of Table 2, with r 10 = 5000 m and c 10 = 1), whereas in Figure 7b, the finally produced SDIS map, where all viewshed classes of Table 2 are taken into account, is presented. In Figure 7a, the symbol SDIS 10 has been used in the legend to denote results that have been obtained for all N observers (Equation (1)), but for DIS i,pix = c 10 V i 10,pix , i = 1, . . . , N in Equation (2). The results of Figure 7a indicate large SDIS 10 values in the high elevation areas, since those areas correspond to the most visible ones from a large distance. However, it can be seen that the SDIS map (Figure 7b) offers better results in terms of disturbance, since it accounts for the largest effect of the nearest objects on the visibility in an efficient manner. Accordingly, large SDIS values are bounded in areas located close to the roads and the residencies. and "plotly" [66].

Results and Discussion
In this section, the results of the present analysis are presented and discussed. Starting with the thematic maps of the exclusion criteria, Figure 6a shows the incompatibility zones of EC1-EC4 (Table 1), whereas in Figure 6b, the slope terrain thematic map (EC5, Table 1) is cited. Areas that fall within the incompatibility zones of Figure 6a cover 26.9 km 2 and belong to set . Most of these areas correspond to environmentally protected areas, followed by minor roads and non-agricultural areas. With regard to the viewshed analysis, Figure 7a shows a cumulative weighted viewshed map indicatively for a maximum observation distance of 5000 m (10th viewshed class of Table 2, with = 5000 m and = 1), whereas in Figure 7b, the finally produced SDIS map, where all viewshed classes of Table 2 are taken into account, is presented. In Figure 7a, the symbol SDIS10 has been used in the legend to denote results that have been obtained for all observers (Equation (1)), but for , = , , = 1, … , in Equation (2). The results of Figure 7a indicate large SDIS10 values in the high elevation areas, since those areas correspond to the most visible ones from a large distance. However, it can be seen that the SDIS map (Figure 7b) offers better results in terms of disturbance, since it accounts for the largest effect of the nearest objects on the visibility in an efficient manner. Accordingly, large SDIS values are bounded in areas located close to the roads and the residencies. As for the optimization results, the minimum values of the objective function described by Equation (6) (i.e., sum of the SDIS values over all optimum areas) obtained for the 100 different and combinations are summarized in Table 4. In this table, values ≥ 2.5 correspond to penalty values and denote and combinations without any possible solutions (i.e., = {}). The results of Table 4 indicate that for a given value, the increase of the minimum allowable total coverage area for PV installations ( ) generally leads to larger SDIS values. This is attributed to the fact that by increasing , larger areas suitable for installing PVs can be allocated, leading to a larger social disturbance. The opposite holds true for , since for a given , the increase of reduces the space suitable for PV installations, and, thus, As for the optimization results, the minimum values of the objective function described by Equation (6) (i.e., sum of the SDIS values over all optimum areas) obtained for the 100 different DG max and Area min combinations are summarized in Table 4. In this table, values ≥ 2.5 correspond to penalty values and denote DG max and Area min combinations without any possible solutions (i.e., A = {}). The results of Table 4 indicate that for a given DG max value, the increase of the minimum allowable total coverage area for PV installations (Area min ) generally leads to larger SDIS values. This is attributed to the fact that by increasing Area min , larger areas suitable for installing PVs can be allocated, leading to a larger social disturbance. The opposite holds true for DG max , since for a given Area min , the increase of DG max reduces the space suitable for PV installations, and, thus, it leads, in general, to smaller SDIS values. The effect of the DG max and Area min siting preferences on the optimization results is also shown schematically in Figure 8, where the minimum SDIS values of Table 4 for the 67 non-empty solution sets are plotted. It is also interesting to note that empty solution sets are mainly obtained for numerous Area min values, when DG max ≤ 3.5 km. For example, for DG max = 2.5 km, only one optimum solution has been obtained in the case of the smallest examined Area min (equal to 0.5 km 2 ). This, in turn, illustrates that for small DG max values, extensive areas for PV installations cannot be found in the region.     Figure 9 shows the entities of the A set (i.e., optimum solutions) for DG max = 7.0 km and Area min = 0.5 km 2 (smallest examined low bound of PV locations' total coverage area), whereas in Figure 10, the A set entities corresponding to the aforementioned DG max , but to Area min = 5.0 km 2 (largest examined low bound of PV locations' total coverage area), are presented.  It can be easily concluded that for both siting preferences combinations, the optimum areas for PV installations correspond to agricultural parcels that have low visibility, while  It can be easily concluded that for both siting preferences combinations, the optimum areas for PV installations correspond to agricultural parcels that have low visibility, while It can be easily concluded that for both siting preferences combinations, the optimum areas for PV installations correspond to agricultural parcels that have low visibility, while not intersecting with the incompatibility zones of any of the exclusion criteria considered herein. Furthermore, distinctive compact regions of closely located agricultural parcels are formed. In the case of Area min = 0.5 km 2 (Figure 9), optimum locations are bounded in the north-eastern region of the study area. Some of these locations also correspond to optimum solutions when Area min = 5.0 km 2 ( Figure 10). However, for this Area min value, a larger number of optimum areas are obtained, which are also distributed in the form of compact regions within the whole examined study area. Although only two triangles are deployed in the GAs algorithm to identify optimum spatial entities (see Section 3.3), the results of Figure 10 clearly indicate that distinguishable compact regions of optimum locations for PV installations can be realized. This is attributed to the ability of the algorithm to create large-size triangles, which, at the end of the iterative optimization process, will include only spatial entities belonging to set A.
Finally, Figure 11 shows a snapshot of the web-GIS application developed by taking into account the results of the GAs-driven optimization process for all the examined DG max and Area min combinations. The user can: (a) set, via bars, his/her preferences regarding the values of DG max and Area min ; (b) visualize all the results (optimum areas, as well as spatial entities belonging to B and C sets); and (c) download a matrix that contains all spatial information (e.g., official IDs) related to the areas of set A. not intersecting with the incompatibility zones of any of the exclusion criteria considered herein. Furthermore, distinctive compact regions of closely located agricultural parcels are formed. In the case of = 0.5 km 2 (Figure 9), optimum locations are bounded in the north-eastern region of the study area. Some of these locations also correspond to optimum solutions when = 5.0 km 2 ( Figure 10). However, for this value, a larger number of optimum areas are obtained, which are also distributed in the form of compact regions within the whole examined study area. Although only two triangles are deployed in the GAs algorithm to identify optimum spatial entities (see Section 3.3), the results of Figure 10 clearly indicate that distinguishable compact regions of optimum locations for PV installations can be realized. This is attributed to the ability of the algorithm to create large-size triangles, which, at the end of the iterative optimization process, will include only spatial entities belonging to set . Finally, Figure 11 shows a snapshot of the web-GIS application developed by taking into account the results of the GAs-driven optimization process for all the examined and combinations. The user can: (a) set, via bars, his/her preferences regarding the values of and ; (b) visualize all the results (optimum areas, as well as spatial entities belonging to and sets); and (c) download a matrix that contains all spatial information (e.g., official IDs) related to the areas of set . Figure 11. A view of the web-GIS application for selecting optimum locations for PV installation in La Palma del Condado municipality.

Conclusions
In this paper, an integrated methodology has been developed to identify optimum areas for PV installations on municipality scale that minimize the relevant visual disturbance and satisfy spatial constraints related to land use, environmental, and techno-economic siting factors. The visual disturbance of the public due to PV installations is quantified by introducing and calculating the SDIS indicator, whereas optimum locations are determined for the predefined values of two siting preferences (maximum allowable PV locations-grid station distance; and minimum allowable total coverage area of PV installations). In addition to the development of standard thematic maps of exclusion criteria, a cumulative weighted viewshed analysis is realized to quantify the SDIS indicator, and

Conclusions
In this paper, an integrated methodology has been developed to identify optimum areas for PV installations on municipality scale that minimize the relevant visual disturbance and satisfy spatial constraints related to land use, environmental, and techno-economic siting factors. The visual disturbance of the public due to PV installations is quantified by introducing and calculating the SDIS indicator, whereas optimum locations are determined for the predefined values of two siting preferences (maximum allowable PV locationsgrid station distance; and minimum allowable total coverage area of PV installations). In addition to the development of standard thematic maps of exclusion criteria, a cumulative weighted viewshed analysis is realized to quantify the SDIS indicator, and then optimum solutions are determined by developing a GAs-driven optimization process. The proposed methodology is applied for the municipality of La Palma Del Condado located in the province of Huelva, Andalucia, Spain, for 100 different combinations of DG max (maximum allowable PV locations-grid station distance) and Area min (minimum allowable total coverage area of PV installations). The main conclusions of the present investigation can be summarized as follows: • The map of the proposed SDIS indicator can be easily created by both researchers and practitioners with low computational cost, and it accounts for the larger effect of the nearest objects on the visibility in an efficient manner. Accordingly, it can offer more realistic results than traditional viewsheds for assessing the visual effect of PV installations to the public. • For a given DG max value, the increase of Area min facilitates the allocation of larger optimally suitable areas for installing PVs; thus, the aforementioned increase leads generally to larger SDIS values. The opposite holds true for DG max , since, from a physical point of view, the increase of DG max for a given Area min value reduces the space suitable for PV installations. For the examined study area, the GAs-driven optimization process has led to empty optimum solution sets for numerous Area min values, especially when DG max ≤ 3.5 km, demonstrating that for small DG max values, extensive areas for PV installations cannot be found in the region.

•
The developed GAs-driven optimization process offers the ability to determine distinguishable, but compact, regions of optimum locations for PV installations within the examined region, facilitating relevant regional planning decisions. The consideration of the SDIS indicator in the objective function can contribute to the mitigation of potential social oppositions and negative impacts on land activities, since optimum areas correspond to those that will have the minimum visual impact to the community.

•
The developed web-GIS application presents a flexible and easy-to-use tool that enables the visualization of PV plants' optimum locations in the study area for different bounds of the PV locations-grid station in-between distance and of the PV locations' total coverage area. Accordingly, it facilitates policy-makers to choose the set of solutions that better fulfils their preferences/strategies related to the above factors. The flexibility of the tool can also contribute to the reduction of bureaucracy, as well as to the further boost of the local solar energy market in an environmentally friendly and socially accepted manner.
The proposed methodology is general and can be easily applied to other regions in or outside Spain in order to identify optimum areas for PV installations by appropriately modifying, if necessary, the set of the exclusion criteria according to the spatial conditions, legislation limitations, and/or any special features of the examined region. Furthermore, the calculation of the SDIS indicator by considering additional factors that can affect visibility (e.g., temporal variations), as well as the deployment of another evolutionary optimization method to solve the examined optimization problem and the comparison of results, including computational cost aspects, with the ones of the present study, could represent items for future research.