The Combination of Expert Judgment and GIS-MAIRCA Analysis for the Selection of Sites for Ammunition Depots

Ljubomir Gigović 1, Dragan Pamučar 2, Zoran Bajić 3,* and Milić Milićević 2 1 Department of Mathematics, University of Defence, Belgrade 11000, Serbia; ljubomir.gigovic@va.mod.gov.rs 2 Department of Logistics, University of Defence, Belgrade 11000, Serbia; dpamucar@gmail.com (D.P.); milic.milicevic@mod.gov.rs (M.M.) 3 Department of Military-Chemical Engineering, University of Defence, Belgrade 11000, Serbia * Correspondence: zoran.bajic@va.mod.gov.rs; Tel.: +381-11-3603259


Introduction
In recent years, much effort has been made to solve facility location problems (FLP).Furthermore, many additional features have been comprised in location models including multiple objectives and multiple facility types.Facility location problems continue to be an important part of strategic planning and decision making for logistic managers today.FLP should identify the best set of storage locations with the best mix of munitions' inventories that considers transportation and other logistical restrictions in satisfying a set of potential demands [1].The majority of the research in this area has been related to the problem of where to locate municipal facilities.However, several references are used to identify adequate locations for military related facilities that store munitions around the world [2][3][4][5][6][7][8][9].
Numerous Unplanned Explosions at Munitions Sites (UEMS) attract progressively more public attention due to the existing risk built from inadequate management of explosive materiel and its storage at the unsafe locations.UEMS leave behind numerous direct and indirect consequences.Along with the casualties and severely damaged buildings, there are additional expenses in medical care, reconstruction and rebuilding costs, unexploded ordnance removal costs, harmful impact on the environment and business profits loss.[10].
In 2006, in the Serbian army ammunition depot (AD) near Paraćin town, a huge explosion occurred, supposedly due to the inappropriate storage of ammunition.The explosion directly resulted in the total destruction of stored ammunition, the injury of 23 people, the pollution of the surrounding 30 square kilometers with nearly 90,000 unexploded ordnance items, the damage and demolition of 600 objects and 12 schools, and glass breakage in around 4740 buildings in the Paraćin town and around 2000 buildings in the nearby town of Ćuprija.Broader socioeconomic concerns were evident with the temporary highway and railroad shutting down and the evacuation of people from the nearby locations [11].These facts clearly point out the essential need to locate and build ammunition depots at safer locations, without losing their functional purposes.
The scope of this paper points to the wider public necessity for locating and building ammunition depots according to the legislations and law in order to establish safe surroundings and to preserve the environment.Therefore, finding suitable locations for ammunition depots is rather a complex issue requiring careful, joint analysis of numerous criteria.
Optimum selection of ammunition depot locations requires consolidation of mutually opposed economic and technical factors with ecological and safety restrictions.In other words, decision makers (DM) are confronted with a dual challenge.They have to satisfy the technological requirements of building AD while taking into account the economic aspects.Furthermore, it is needed to minimize risk to the environment and to meet location safety standards.In order to do this, it is necessary to set certain rules which will help DMs to evaluate different locations based on certain spatial, safety, environmental and economic criteria and restrictions.
Geographic information technologies (GIS) are ideal for this kind of studies because they can efficiently manage large numbers of spatial and attributive data collected from various sources.GIS have certain flaws when used for location selection, like the lack of multi-criteria decision making (MCDM) possibilities and optimum location recommendations.It is essential to combine tools for MCDM with GIS in order to obtain MCDM as well as spatial outputs.The combination of GIS and multi-criteria decision analysis (MCDA) is a widely accepted concept for spatial problem solving in many areas such as the environment, ecology, transportation, urban and regional planning, waste management, hydrology, agriculture, forestry, geology and site selection [12].
GIS-MCDA is currently the most advantageous concept used to assist DM in selecting the optimal choice from a number of feasible options/alternatives with multiple choice criteria and diverse criterion priorities.The multi-criterion choice can be attributed to spatial decision-making problems involving the search for the optimal AD site.Similar problems, often analyzed in GIS, include site selection for critical ecologically sensitive areas, hazardous waste disposal sites and emergency service locations [13][14][15][16][17][18][19].A large number of studies combine GIS with different MCDM: AHP [20,21], ANP [22,23], TOPSIS [24,25], PROMETHEE [26,27], VIKOR [28,29], etc.This paper presents the GIS-MAIRCA approach, thus investigating the new MCDM method and its potential use with GIS for the AD location problem.This paper comprises of the following parts: Definition of aim/problem and model structure; Identification of restrictions and evaluation criteria; Data gathering and building of the spatial criteria basis; GIS-MCDA evaluation; Sensitivity analysis and GIS visualization of the final solutions and recommendations.

Methodological Background
The methodological hierarchy model presented in this paper is based on spatial MCDA structure.Spatial multi-criteria decision analysis (SMCDA) consists of a number of procedures including the use of geographic data and DM preferences, and the manipulation of data and preferences according to specified decision rules [12].This approach uses the capabilities of GIS in spatial data management and the flexibility of MCDA to combine factual information (e.g., terrain slope, road and railways, land cover/use, etc.) with value-based information (e.g., expert assessments, standards, opinion polls, etc.) [30].The main advantage of GIS and MCDA integration is the ability to obtain unique capabilities which are complementary to each other.GIS has great possibilities for the manipulation, storage, management, analysis and visualization of geospatial data while MCDA offers an assortment of procedures, techniques and algorithms for decision making structure and design, and evaluating and prioritizing options [31].
From a methodological point of view, the suggested model for defining suitable locations for AD consists of the following steps (Figure 1): Sustainability 2016, 8, 372 3 of 30 the manipulation, storage, management, analysis and visualization of geospatial data while MCDA offers an assortment of procedures, techniques and algorithms for decision making structure and design, and evaluating and prioritizing options [31].
From a methodological point of view, the suggested model for defining suitable locations for AD consists of the following steps:

Definition of Aim/Problem and Model Architecture
Identification of suitable AD locations is the general aim of this paper.DM must justify their choice regarding AD location selection through a multisystem, transparent and documented process.The building of network hierarchy models involves an important creative phase in the problem solving process because the consideration of all relevant factors and their mutual interactions is the key to obtaining the correct solution.
The MCDA model is suggested as a relevant option for systematic analysis and rational decision making due to the complexity of the location selection problem.The MCDA model is based on the joint use of GIS and multi-criteria techniques such as DEMATEL, ANP and MAIRCA.
Authors have chosen the hybrid DEMATEL-ANP model because it eliminates subjectivity while generating model architecture, unlike the classic ANP method.The DEMATEL method creates a cause and effect diagram based on a mathematic formulation.Further on, this relationship is used as an initial data set for the ANP method.This paper also presents a novel MCDM method-MAIRCA.This method is proven to be more stable than other methods, like TOPSIS or ELECTRE.One reason is the different criteria normalization.In [32], it was shown that methods using a linear model of input data normalization have greater stability and rank consistently during sensitivity analysis.The MAIRCA method uses a linear normalization model.The second reason for choosing the MAIRCA method is its simple mathematic apparatus, solution stability and the possibility to combine this method with other ones [32].

Definition of Aim/Problem and Model Architecture
Identification of suitable AD locations is the general aim of this paper.DM must justify their choice regarding AD location selection through a multisystem, transparent and documented process.The building of network hierarchy models involves an important creative phase in the problem solving process because the consideration of all relevant factors and their mutual interactions is the key to obtaining the correct solution.
The MCDA model is suggested as a relevant option for systematic analysis and rational decision making due to the complexity of the location selection problem.The MCDA model is based on the joint use of GIS and multi-criteria techniques such as DEMATEL, ANP and MAIRCA.
Authors have chosen the hybrid DEMATEL-ANP model because it eliminates subjectivity while generating model architecture, unlike the classic ANP method.The DEMATEL method creates a cause and effect diagram based on a mathematic formulation.Further on, this relationship is used as an initial data set for the ANP method.This paper also presents a novel MCDM method-MAIRCA.This method is proven to be more stable than other methods, like TOPSIS or ELECTRE.One reason is the different criteria normalization.In [32], it was shown that methods using a linear model of input data normalization have greater stability and rank consistently during sensitivity analysis.The MAIRCA method uses a linear normalization model.The second reason for choosing the MAIRCA method is its simple mathematic apparatus, solution stability and the possibility to combine this method with other ones [32].

Identification of Restrictions and Evaluation Criteria
The location selection process is a MCDA problem requiring adequate consideration of the comprehensive criteria set in order to determine the suitability of certain areas for a particular land use.The MCDA process uses two criteria categories: restrictions and evaluation criteria.Restrictions are based on the criteria which limit (exclude) possible alternatives and that are based on Boolean relations (true/false).The evaluation criteria can be quantified according to the degree of suitability for all feasible alternatives [33,34].Criteria for the assessment of AD location suitability were adopted according to the established goals and based on UN and national regulations along with the leading military publications and literature [35][36][37][38][39][40][41].They comprise the following factors: technical, safety, environmental, spatial and socioeconomic.In order to determine restriction criteria, it is necessary to study the territory status carefully (settlements, roads, rivers, protected areas, etc.), along with the regulations (European, national, regional and local) and further determine the boundary values which will exclude some alternatives from the future analysis.

Data Gathering and Building of the Spatial Criteria Basis
It is necessary, during this step, to form individual map layers which correspond to each identified restriction and evaluation criteria.This can be achieved by different GIS operations used to enter, accept and transform spatial and thematic data into digital forms (map overlay, buffering, distance mapping, etc.).Initial data mainly come from digital and analog cartographic documents, terrain images, blueprints, statistic reports, etc. GIS data criteria are usually organized in vector or rasterized format.Successive arithmetic operations in GIS (reclassification, overlay functions, Boolean algebra, etc.) criteria maps are converted to the rasterized format, where every raster cell presents an alternative towards the goal.Calculations based on rasterized data formats are far less complicated than those on vector data [42].

GIS-MCDA Evaluation
This phase concerns standardization and evaluation of all criteria used in the analysis.

Formation and Aggregation of Restrictions Map
The implementation of restrictions in the thematic maps is accomplished using Boolean (logical) algebra, in the way that raster cells (alternatives) for the excluded area have zero value (0), while the other alternatives are given a value of one (1).The final constrained area is obtained through aggregation of all layers (criteria) with restrictions.Only the locations which are suitable on all restriction levels in the restrictions map can be considered as a potential location for the construction of AD and therefore are the subject of further analysis.

Individual Assessment of Evaluation Criteria and Standardization
Criteria are presented on maps in GIS layers in various ways and in various forms.The method of weighted linear combination (WLC) requires that all data criteria is standardized [33] or transformed (reclassified) into mutually comparable units.There are many useful approaches that can be used so the attributive layers are comparable [31].The decision as to which approach should be applied depends on the attribute nature and on DM.The most frequently used approach is linear transformation, or transformation on the interval scale, when every raster cell (attribute) gets a nondimensional value x i on the adopted scale in terms of the aim of the analysis.Hence, every attribute acquires an integer value x i according to the class i association.It is common in this method to determine the midpoint of the attribute interval because it represents the boundary between the desirable and undesirable solution.In real conditions, this method gives very good results even though it is quite arbitrary.

Comparison and Relative Assessment of Criteria Weight (DEMATEL-ANP)
While solving real problems, criteria do not have the same degree of significance so the DM has to define significance factors for some criteria using weighting coefficients (weights) or criteria evaluation.The WLC method demands normalized weights, meaning that the sum of all weights must be 1.
DEMATEL and ANP (DANP) method is used for the calculation of normalized criteria weights followed then by WLC method.The DEMATEL method was originally developed by the Science and Human Affairs Program of the Battelle Memorial Institute of Geneva, with the purpose of studying a complex and intertwined problem [43].It has been widely accepted as one of the best tools to solve the cause and effect relationship among evaluation criteria [43][44][45][46][47].The original DEMATEL method is modified by some researchers so as to make it comply with their problems.The modified DEMATEL method which is used in this paper is adapted from the study of [48] (Figure 2).While solving real problems, criteria do not have the same degree of significance so the DM has to define significance factors for some criteria using weighting coefficients (weights) or criteria evaluation.The WLC method demands normalized weights, meaning that the sum of all weights must be 1.
DEMATEL and ANP (DANP) method is used for the calculation of normalized criteria weights followed then by WLC method.The DEMATEL method was originally developed by the Science and Human Affairs Program of the Battelle Memorial Institute of Geneva, with the purpose of studying a complex and intertwined problem [43].It has been widely accepted as one of the best tools to solve the cause and effect relationship among evaluation criteria [43][44][45][46][47].The original DEMATEL method is modified by some researchers so as to make it comply with their problems.The modified DEMATEL method which is used in this paper is adapted from the study of [48] (Figure 2).The full procedure of the previously defined DEMATEL method is explained as follows: Step 1. Gather experts' opinions and calculate the average matrix Z.The modified DEMATEL method is built first creating the direct-relation matrix Z where Z is (n × n) matrix, n represents the number of criteria.The direct relation matrices are all obtained by holding a pair-wise comparison among the criteria themselves in which Zij indicates the degree to which the criterion i C affects criterion j C .Therefore, the relationship among the criteria could be held within a matrix.In this method, the effects of the criteria on each other are expressed in terms of linguistic equations.
A group of e experts and n factors are used in this step.Each expert is asked to view the degree of direct influence between two factors based on a pair-wise comparison.The degree to which the expert perceived factor i impacts factor j is denoted as ij z .For each expert, a n × n non-negative matrix is constructed (n represents the number of criteria) as e e ij Z z      , where e is the expert number of participating in evaluation process with 1 ≤ e ≤ k.Thus, 1 2  , , ..., k Z Z Z are the matrices The full procedure of the previously defined DEMATEL method is explained as follows: Step 1. Gather experts' opinions and calculate the average matrix Z.The modified DEMATEL method is built first creating the direct-relation matrix Z where Z is (n ˆn) matrix, n represents the number of criteria.The direct relation matrices are all obtained by holding a pair-wise comparison among the criteria themselves in which Z ij indicates the degree to which the criterion C i affects criterion C j .Therefore, the relationship among the criteria could be held within a matrix.In this method, the effects of the criteria on each other are expressed in terms of linguistic equations.
A group of e experts and n factors are used in this step.Each expert is asked to view the degree of direct influence between two factors based on a pair-wise comparison.The degree to which the expert perceived factor i impacts factor j is denoted as z ij .For each expert, a n ˆn non-negative matrix is constructed (n represents the number of criteria) as Z e " , where e is the expert number of participating in evaluation process with 1 ď e ď k.Thus, Z 1 , Z 2 , ..., Z k are the matrices from m experts.In this method, the effects of the criteria on each other are expressed in terms of linguistic equations.Aggregation of experts' opinions gives the final matrix Z " " where z e ij denotes preference of expert number e, k is total number of experts.where the elements of the initial direct-relation matrix is calculated from R " max where n denotes total number of criteria.
Step 3. Derive the total relation matrix T. The total-influence matrix T is obtained by utilizing Equation (7), in which, I is an n ˆn identity matrix.The element of t ij represents the indirect effects that factor i had on factor j, then the matrix T reflects the total relationship between each pair of system factors.Initial normalized direct-relation matrix can be separated into separate sub matrices i.e., (D 1 ,D 2 ,D 3 ).It was proven that lim kÑ9 pD s q w " 0 And lim kÑ9 ´I `Ds `... `D2 s `... `Dk s ¯" pI ´Ds q ´1 (6 where 0 is the null matrix and I is the identity matrix [49,50].Therefore, the total-relation matrix T can be acquired by calculating the following term [51]: As it is previously shown, the total relation matrix for the criteria (T) can be presented as: where t ij is the overall influence rating of decision maker for each alternative i against criterion j.
Step 4. Calculate the sums of rows and columns of matrix T. The sum of rows and sum of columns of the sub-matrices T 1 , T 2 and T 3 denoted by the numbers D i and R i , respectively, can be obtained through Equations ( 9) and (10), respectively: where n denotes number of criteria.
Step 5. Determination of the threshold value α.Threshold value α is used for creating a cause and effect relationship diagram.Threshold value α is calculated as the median of elements (t ij ) from the matrix T, according to Equation (11).This calculation aims eliminate some minor effects of elements in matrix T [52].
α " where N is the total number of elements in matrix T.
Step 6. Creating a cause and effect relationship diagram.The cause and effect diagram is constructed by mapping all coordinate sets of (D i + R j , D i ´Rj ) to visualize the complex interrelationship and provide information in order to determine the most important factors and how they influence the affected factors [53].The factors t ij that are greater than α are selected and shown in the cause and effect diagram [52].
After designing the cause and effect relationship diagram, weighted coefficients of criteria are calculated using the ANP method in the next phase.Afterwards, the implementation of DEMATEL method into the ANP method is elaborated.
The ANP method was developed to avoid the hierarchical restrictions that exist in the AHP method [54].ANP presents a generalized AHP method, designed so that the hierarchy is replaced by a network.The ANP method, unlike hierarchy structured problems, takes into account different forms of dependency and feedback.The structure of the feedback is not linear and is more close to the network where interdependent loops frequently appear.The calculation of relative weights of criteria using traditional ANP means the levels of interdependence of factors are treated as reciprocal values.On the contrary, using the DEMATEL method, the levels of interdependence of factors do not have reciprocal value, which is closer to real circumstances [55].
The traditional ANP method calculates weighted supermatrix by normalization of unweighted supermatrix.Each column of the unweighted supermatrix is divided by the number of criteria so that each column will add up to 1.It means that every criterion has the same weight.However, then, that is not a good assumption, because the mutual influence of two criteria can be different.Although it is easy to normalize using such a method, this neglects the fact that different groups should have different degrees of influence.Thus, we need to find another way of normalizing the unweighted supermatrix that relaxes this assumption of equal weight among criteria.Here, we turn to the total-influence matrix T in DEMATEL and its threshold value α for help [56].
In order to address the mentioned faults, the total relation matrix (matrix T) is used to calculate the weighted coefficient of criteria, according to the DEMATEL method (step 3, Equation ( 8)).Hence, the mentioned problem is easily resolved by merging DEMATEL with the ANP method (DANP) so the obtained results are more reflective of the real situation.The hybrid DANP method, or the integration of DEMATEL and ANP methods, is conducted through five steps.They are completed after finishing the effect relationship diagram (ERD).
Elements of the unweighted supermatrix are defined in the first step.Before unweighted supermatrix design is completed, it is necessary to define the network model for the ANP method.The network model is designed on the basis of the total relation matrix and ERD.
Unweighted supermatrix is created when each level with the total degree of influence from the total relation matrix T is normalized.Total sum of matrix elements by columns needs to be calculated before its normalization.
where T 11 c matrix contains factors from group K1 and influences in respect to the factors from group K1.Elements of the normalized total influence matrix for criteria T c α are calculated in the next step.After that, the normalization of the matrix T c is done.Normalized matrix T c α is shown in Equation ( 13) For an example, an explanation for the normalization of T c α11 on dimension K1 is shown with Equations ( 14) and (15).Factor sum c 11 , ..., c 1m1 inside the group K1 is obtained from Equation (14).
where t cj 11 denotes values of factor influences c 11 , ..., c 1m1 in relation to factors from group D1, and elements t c11 α11 and their normalized values.
The calculation procedure for other matrixes T c αnn inside the matrix T c α is identical and it will not be further explained.
The calculation of the unweighted supermatrix W is completed in the third step.Because the total influence matrix T c and fills the interdependence among dimensions and criteria, we can transpose the normalized total influence matrix T c α by the dimensions based on the basic concept of ANP resulting in the unweighted supermatrix W = (T c α )', as shown in Equation ( 16) The normalization of the total influence matrix T α D of the dimensions, Equation ( 17), gave a new normalized total influence matrix T α D , Equation (18).
The limit of the weighted supermatrix W α is calculated in the last step.Afterwards, the weight of each evaluation criteria is obtained.The weighted supermatrix can be raised to the limiting powers until the supermatrix has converged and becomes a long-term stable supermatrix.In this way, the global priority vectors, called DANP influence weights, were obtained such as lim kÑ8 " W k , where W denotes the limit supermatrix, while k denotes any integer.

Aggregation Results (Suitability Map)
With GIS-MCDA, during the process of criteria aggregation, Weighted Linear Combination (WLC) is generally used.The WLC aggregation method multiplies each standardized factor map (i.e., each raster cell within each map) by its factor weight and then sums up the results.Since the sum of the set of factor weights for an evaluation must be one, the resulting suitability map will have the same range of values as the standardized factor maps that were used.The cell (alternative) to which corresponds the maximum calculated value is "the best" solution.The following mathematical equation was used to combine the evaluation criteria (factors) according to the weighted linear combination method: where S is suitability, w i is normalized weight of factor i, and x i is the criterion score of factor i. In cases where Boolean restrictions also apply, the procedure can be modified by multiplying the suitability calculated from the factors by the product of the restrictions: where c j is the criterion score of the restriction j.
All GIS software systems provide the basic tools for evaluation of such models [57].WLC is relatively easy to comprehend and it can be applied in different situations.In addition, it is compensatory which means that low scores in one criterion can be compensated by high scores in another one, if that is desired for certain decision problems.For these reasons, WLC was selected as the method of aggregation.

Extraction and Ranking of Suitable Locations
The extraction of suitable locations presents the extraction of cells with the highest x i values from the final suitability map.They represent potentially the most suitable alternatives for the AD locations.The filtration of the most suitable alternatives (cells with the highest x i ) is conducted with the combined use of GIS arithmetic operations and inquiries.As a result, the suitable locations were obtained which are converted in the vector form meeting a certain level of suitability.
The final decision or recommendation should be based on the suitable locations' ranking.The ranking is completed in order to assess every location with the established goals of the analysis.This model proposes the use of MAIRCA (MultiAttributive Ideal-Real Comparative Analysis) method for the selection of the most suitable location.MAIRCA was developed in 2014 by the Center for Logistics Research at the University of Defence in Belgrade [58].
The main assumption of MAIRCA method is in the determination of the gap between ideal and empirical weights.The summation of gaps for each criterion gives the total gap for every observed alternative.At the end, it comes to the alternatives ranking, where the best ranked alternative is the one with the smallest value for the total gap.The alternative with the smallest value for the total gap is the alternative which had values closest to the ideal weights according to the largest number of criteria (ideal criteria values).The MAIRCA method is conducted in six steps: Step 1.The formation of the initial decision matrix (X).The criteria values (x ij , i = 1, 2, . . ..n; j = 1, 2, . . .m) are determined in the initial decision Equation (22) for each of the observed alternatives.Equation ( 22) are obtained according to the personal preferences of the DM or by aggregation of experts' decisions.

X "
Step 2. Preference determination according to the alternative selection P A i .During alternative selection, the DM is neutral to the process.In fact, he does not have preference to any of the proposed alternatives.The main assumption is that DM does not take into account the probabilities of any alternative selection.The DM perceives the alternatives as if any of them can have equal probability of appearance, so the preference to choose one of them from m possible alternatives is: where m denotes the total number of choices.
In the decision making analysis, with the given probabilities we assume that the DM is neutral towards the risk.In that case, all the preferences are equal according to the selection of certain alternatives where m denotes the total number of choices.
Step 3. The calculations of theoretical evaluation matrix elements (T p ). Theoretical evaluation matrix is created (T p ) with format n ˆm (n is the total criteria number, m is the total number of alternatives).Theoretical evaluation matrix elements (t pij ) are calculated as the multiplication of the preferences according to the alternatives P A i and criteria weights (w i , i " 1, 2, ..., n) T p " Because the DM is neutral to the initial selection of the alternatives, all preferences (P A i ) are equal for all alternatives.Then, the Equation ( 25) can be shown in the Equation ( 26) where n is the total number of criteria, t pi is theoretical evaluation.
where n is the total number of criteria, m is the total number of alternatives.
The calculation of real evaluation matrix elements (T r ) is conducted by multiplying the real evaluation matrix elements (T p ) and elements of the initial decision making matrix (X) according to the equation: (a) For the "benefit" type criteria (the bigger criterion value is desirable) For the "cost" type criteria (the smaller criterion value is desirable) where x ij , x ì i x í denote elements of the initial decision making matrix (X), x ì and x í are defined as: x ì " max px 1 , x 2 , ..., x m q are the maximum values of the marked criterion by its alternatives.
x í " min px 1 , x 2 , ..., x m q are the minimum values of the marked criterion by its alternatives.
Step 5.The calculation of the total gap matrix (G).The elements of the matrix G are calculated as the difference (gap) between theoretical (t pij ) and real evaluations (t rij ), or by actually subtracting the elements of the theoretical evaluation matrix (T p ) with the elements of the real evaluation matrix (T r ).
where n is the total number of criteria, m is the total number of alternatives.
Gap g ij takes values from the interval g ij P " 0, `tpij ´trij ˘‰, according to Equation ( 31) It is desirable that g ij goes to zero ( g ij Ñ 0 ) because the alternative with the smallest difference between the theoretical (t pij ) and real evaluation (t rij ) is chosen.If the alternative A i for the criterion C i has a theoretical evaluation value equal to the real evaluation value (t pij " t rij ), then the gap for the alternative A i for the criterion C i is g ij " 0. In fact, the alternative A i for the criterion C i is the best (ideal) alternative (A ì ).
If the alternative A i for the criterion.C i has a theoretical evaluation value.t pij and the real ponder value t rij " 0, than the gap for the alternative.A i for the criterion C i is g ij " t pij .In fact, the alternative A i is for the criterion C i the worst (anti-ideal) alternative (A í ).
Step 6.The calculation of the final values of the criteria functions (Q i ) for the alternatives.The values of the criteria functions are obtained from the sum of gaps (g ij ) for the alternatives, or just the sum of the matrix elements (G) in columns, using Equation (32)

Sensitivity Analysis
The following step is to examine the solution stability (for the best alternatives) in the change of certain initial data.The purpose of this is to identify the change effects of entry data (geographic data and DM preferences) to exit data (the best alternative).It is reasonable that the most interesting aspect is examining the solution stability of the changes of relative weights, which are representative of subjectivity in the MCDA.

GIS Visualization of the Final Solutions and Recommendations
GIS visualization of the final solutions for the location selection according to the DM and interest groups is generally conducted through cartographic material, graphs, tables, etc. [31].

Studied Area
Regarding the geomorphological structure of the mountain systems in Serbia, the most suitable locations for the AD are situation at the Carpathians-Balkan mountain ranges.The Carpathians mountain massif is especially suitable due to its location, significance and karst geological structure meaning that the terrain fits well with the AD operational criteria.This massif is the biggest karst area in Serbia.The existence of a military range in this region also supports the suitability of establishing AD in this area.
The Carpathian area is located in the eastern part of Serbia and it is geographically set between 43 ˝45 1 and 44 ˝15 1 north latitude and between 21 ˝26 1 and 22 ˝04 1 east longitude.The territory has 1643 km 2 .The whole area is defined by a plateau character with a complex tectonic and morphological structure along with distinctive small river valleys.This region is characterized by a greater number of forests.The dominant geological structure is karst.The additional benefit is the existence of local and regional roads and railroads which enables cheaper transport of heavy ammunition loads.(Figure 3).

Studied Area
Regarding the geomorphological structure of the mountain systems in Serbia, the most suitable locations for the AD are situation at the Carpathians-Balkan mountain ranges.The Carpathians mountain massif is especially suitable due to its location, significance and karst geological structure meaning that the terrain fits well with the AD operational criteria.This massif is the biggest karst area in Serbia.The existence of a military range in this region also supports the suitability of establishing AD in this area.
The Carpathian area is located in the eastern part of Serbia and it is geographically set between 43°45′ and 44°15′ north latitude and between 21°26′ and 22°04′ east longitude.The territory has 1643 km 2 .The whole area is defined by a plateau character with a complex tectonic and morphological structure along with distinctive small river valleys.This region is characterized by a greater number of forests.The dominant geological structure is karst.The additional benefit is the existence of local and regional roads and railroads which enables cheaper transport of heavy ammunition loads.(Figure 3).

Identification of Criteria that Influence the Selection of AD Location
According to the analysis of the environmental regulations, safety procedures when handling ammunition, public standards and a literature overview, nine restrictions are adopted in order to

Identification of Criteria that Influence the Selection of AD Location
According to the analysis of the environmental regulations, safety procedures when handling ammunition, public standards and a literature overview, nine restrictions are adopted in order to select a suitable AD location.These restrictions determine which areas should be ruled out from the analysis (Table 1).

R2
Urban areas.The social implications of the settlements that are near the AD locations are considered as a very important restriction due to safety issues.This influence should be assessed according to strict state regulations [39].<2000 m R3 Land use restrictions.Information about land use is a significant restriction.Among the numerous initiatives which provide land use information, the project CORINE is emphasized due to its methodology which allows providing sufficient and consistent information for Europe [60].

R4
Slope.Topography with steep slopes is generally considered as less suitable for AD location due to extra building and functioning costs [39].>15%

R5
Distance from roads and railways.Although AD were built according to the safety and construction standards and have solid construction, practice suggests ensuring a safe distance from roads and railways [39].<500 m R6 Electric power line.Electric power lines are forbidden above and across AD sites due to the high risk of accidents [39].<200 m

R7
The distance from important public objects.Industrial facilities must not be located near AD locations.Also, the presence of the AD can have a negative influence on tourism sites and areas [39].<1000 m

R8
Minerals.It is required that potential locations for the AD do not have metal ores (manganese, iron, etc.) which affect the frequency and the intensity of lightning, so the AD safety would be in jeopardy [39].

R9
Military ranges.Military facilities are under the special safety status of the Ministry of Defence.Due to army units' training and other military activities, it is necessary that AD are at a safe distance from such facilities [39].
<1000 m *: The first level of CORINE hierarchy nomenclature.
Unlike the restrictions, which can be suitable or unsuitable, evaluation criteria are generally assessed on a continuous scale.Their purpose is to express the suitability of certain alternatives (cell)) [61].Table 2 shows the adopted evaluation criteria.

Evaluation
Criteria Description

C1
Distance from roads and railways.Ammunition transport is specific because it requires roads to the AD location along with the use of local roads and railways.So, the new roads should be as short as possible [39].

C2
Distance from urban areas.An AD location near urban areas can have negative influences on the people and the environment.It is required that AD are as far as possible from settlements [39].

C4
Slope.AD should be sited at physically suitable locations.Areas with steep slopes are generally considered as less suitable or unsuitable due to extra expenses required for building and maintenance.It is recommended that chosen areas are as flat as possible [39].

C5
Distance from power lines.Due to the increased fire and explosion hazard, it is necessary that the AD is far from power lines [39].

C6
Distance from military ranges.The distance from the AD to military ranges should be short due to tactical, economic and time saving factors [39].
The presented criteria and restrictions for the AD location are not final.They are adjusted according to the characteristics of the region and the military and state regulations in Serbia [38,39,59,62].Other criteria and restrictions like distance to geological water courses or aquifers, areas which have a propensity to suffer earthquakes or high risk of floods are omitted from this study due to the compact karst nature of the Carpathian Mountains and their monolith pedological constituency so they would not have different spatial influences.

Data Gathering and Entry
Every criteria entered into GIS is presented in the form of spatial vector or rasterized maps, the layers of which present different alternatives with different intervals or scaled values.During data entry and the building of the spatial criteria and restrictions database, all GIS processes (digitalization, conversion, 3D analysis, etc.) are completed using integrated software tools in ArcGIS 10.2.The data sources and the design of GIS data sets are given in Table 3.

GIS-MCDA Evaluation
This phase consists of the aggregation of restrictions, standardization, weighting and sum analysis of all criteria considered in the decision making process.This is shown in Figure 4.
Map with slopes between 0°and 15°

GIS-MCDA Evaluation
This phase consists of the aggregation of restrictions, standardization, weighting and sum analysis of all criteria considered in the decision making process.This is shown in Figure 4.

The Design and Aggregation of the Restriction Map
Aggregation of restrictions given in Table 1 gives "free" territory in the Carpathian region suitable for the AD location.Restrictions presented in the shape of thematic layers are created using Boolean (logical) algebra.Every restriction map is converted into a rasterized form and subsequently reclassified.Raster cell (alternatives) for the excluded area gets zero value while the others get a value of one.(Table 1).The total restriction map is created when all restriction layers are aggregated.Afterwards, the newly created raster cell layer (score 1) represents suitable AD locations (Figure 4).
The analysis of the thematic layer displayed in Figure 4 shows that 325 km 2 or 19.8% of all Carpathian region territory has no AD location limits.The relatively low suitability of this region for

The Design and Aggregation of the Restriction Map
Aggregation of restrictions given in Table 1 gives "free" territory in the Carpathian region suitable for the AD location.Restrictions presented in the shape of thematic layers are created using Boolean (logical) algebra.Every restriction map is converted into a rasterized form and subsequently reclassified.Raster cell (alternatives) for the excluded area gets zero value while the others get a value of one.(Table 1).The total restriction map is created when all restriction layers are aggregated.Afterwards, the newly created raster cell layer (score 1) represents suitable AD locations (Figure 4).
The analysis of the thematic layer displayed in Figure 4 shows that 325 km 2 or 19.8% of all Carpathian region territory has no AD location limits.The relatively low suitability of this region for AD is the result of the huge diversity of the relief scarcity and rigid spatial and safety military required restrictions.This area consists of 519,245 cells (resolution 25 m ˆ25 m) denoting alternatives or the potential sites for the AD location which are the subject of the following suitability evaluation.

Individual Assessment of the Evaluation Criteria and Standardization
For further decision making analysis, detailed criteria maps are to be converted into joint intercomparable forms.Every criteria map is converted into raster and then reclassified according to the standardized class values from Table 4. Thus, every raster cell (alternative) is assigned a value on the shared scale from 1 to 5, where 1 denotes very low location suitability and 5 very high location suitability for the AD location.Obtained reclassified maps are thus converted into operational raster layers appropriate for the WLC.The alternative's suitability in the criteria maps is usually symbolized with different color nuances so the lighter nuances correspond to the least suitable areas while darker ones represent the most suitable alternatives for the AD locations.Experts used the linguistic scale in order to evaluate criteria in the first step of the DEMATEL method application: 0-no influence; 1-low influence; 2-medium influence; 3-high influence; 4-very high influence.The results of the expert's comparison in pairwise criteria are shown in Table 5.
Table 5. Expert comparison matrix in the paired criteria.
Aggregation of experts' decisions is completed according to Equation (1) and the average matrix Z is obtained as it is shown in Table 6.Equations ( 3) and ( 4) are used to normalize elements from matrix Z so the initial direct-relation matrix D is obtained, Equation (2), shown in Table 7. Equations ( 5)-( 7) are used to determine the elements of the total relation criteria matrix (T) (Table 8).Matrix T elements are used in the following step to calculate the cause and effect relationship diagram (CERD).Equations ( 9) and ( 10) sum the matrix T values in rows (D i ) and columns (R i ).The summed values of the matrix T in rows and columns are shown in Table 9.The threshold value was derived from the average of elements in matrix T, which was calculated by using Equation (11).According to the results shown in Tables 4 and 5 and Equation ( 11), CERD is designed (Figure 5) showing the criteria interrelationship.The values of t ij in Table 5, which are greater than α (1.351), are shown as t ij *, which present the interaction between perspectives.As an example, t 12 (1.6022)> α (1.351); hence, the arrow in the cause and effect diagram was drawn from C 1 to C 2 .
Equations ( 5)-( 7) are used to determine the elements of the total relation criteria matrix (T) (Table 8).Matrix T elements are used in the following step to calculate the cause and effect relationship diagram (CERD).9) and ( 10) sum the matrix T values in rows (Di) and columns (Ri).The summed values of the matrix T in rows and columns are shown in Table 9.The threshold value was derived from the average of elements in matrix T, which was calculated by using Equation (11).According to the results shown in Tables 4 and 5 and Equation (11), CERD is designed (Figure 5) showing the criteria interrelationship.The values of tij in Table 5, which are greater than α (1.351), are shown as t ij *, which present the interaction between perspectives.As an example, t 12 (1.6022)> α (1.351); hence, the arrow in the cause and effect diagram was drawn from C 1 to C 2 .The DEMATEL method gave CERD as the result of the data processing which serves as the basis for the ANP method utilization and better understanding relationship network between the criteria.The DEMATEL method gave CERD as the result of the data processing which serves as the basis for the ANP method utilization and better understanding relationship network between the criteria.Interrelationships between the criteria are defined according to the calculations in the matrix T, threshold value α and CERD.The understanding of the network considerably simplifies the understanding of the ANP method modeling used to calculate relative criteria weights.
The first step of the ANP method incorporates the total-influence matrix T into ANP model and, using Equations ( 13)- (18), the unweighted supermatrix is obtained (Table 10).In the following step, the elements of the weighted Equation ( 21) are calculated using Equations ( 19) and (20).Elements of the weighted supermatrix are shown in Table 11.The elements of the limited supermatrix are calculated in the last step of the ANP method, the vectors of which present the weighted coefficients of the model (Table 12).The WLC method is used to sum standardized criteria values, according to Equation ( 21).The criteria weights, obtained from the DANP method, are multiplied with the raster cell score of each criterion creating the final suitability map (Figure 6a).
The final suitability map is presented in the same value interval as for each criterion, from 1 to The higher cell values indicate areas more suitable for the AD location.

Extraction and Ranking of Suitable Locations
The final step of the analysis is the selection of cells with the highest values from the final suitability map.These are cells with pixels that have a value of 5.As the final result of the used model, 5803 potentially suitable locations (alternatives) for the AD are identified in the Carpathian region.Those mainly locations with very small areas of between 400 and 1 km 2 .The total area of the suitable locations is 44.5 km 2 (Figure 6a).
Taking into account that the optimum area for the AD location is a minimum of 0.5 km 2 , according to the explosive safety rules for the medium sized AD [39][40][41], GIS-SQL operation has extracted areas larger than 0.5 km 2 .Eight acceptable locations are identified in the Carpathian region after the analysis is conducted (Figure 6b).
The ranking of the suitable locations is completed after the identification, so the most suitable location is determined using the MAIRCA method.First of all, the initial decision matrix is designed (Table 13), and, afterwards, the criteria normalization is completed using Equations ( 26) and (27).Data shown in Table 13 represent average values of the location suitability obtained using GIS Zonal statistical operation from the reclassified evaluation criteria maps.Because all alternatives are equal and have a value of 1 for the criterion C5, it can be eliminated.It has no influence on the location ranking.Weight criteria is evenly distributed from C5 to other criteria by adding ΔwC5 = wC5/5 = 0.097/5 = 0.0194 to each of them.
After the formation of the initial matrix (Table 13), preferences according to the alternative is calculated using Equation (23).

Extraction and Ranking of Suitable Locations
The final step of the analysis is the selection of cells with the highest values from the final suitability map.These are cells with pixels that have a value of 5.As the final result of the used model, 5803 potentially suitable locations (alternatives) for the AD are identified in the Carpathian region.Those are mainly locations with very small areas of between 400 and 1 km 2 .The total area of the suitable locations is 44.5 km 2 (Figure 6a).
Taking into account that the optimum area for the AD location is a minimum of 0.5 km 2 , according to the explosive safety rules for the medium sized AD [39][40][41], GIS-SQL operation has extracted areas larger than 0.5 km 2 .Eight acceptable locations are identified in the Carpathian region after the analysis is conducted (Figure 6b).
The ranking of the suitable locations is completed after the identification, so the most suitable location is determined using the MAIRCA method.First of all, the initial decision matrix is designed (Table 13), and, afterwards, the criteria normalization is completed using Equations ( 26) and (27).Data shown in Table 13 represent average values of the location suitability obtained using GIS Zonal statistical operation from the reclassified evaluation criteria maps.Because all alternatives are equal and have a value of 1 for the criterion C5, it can be eliminated.It has no influence on the location ranking.Weight criteria is evenly distributed from C5 to other criteria by adding ∆w C5 = w C5 /5 = 0.097/5 = 0.0194 to each of them.
After the formation of the initial matrix (Table 13), preferences according to the alternative is calculated using Equation (23).
where m denotes the number of alternatives.In that case, all preferences are equal and are presented as follows, P A 1 " P A 2 " ... " P A 8 " 0.125 The calculation T p matrix elements (Table 14) is conducted according to Equation (25), t p32 " P A 3 ¨w2 " 0.125 ¨0.230 " 0.0288 After the calculation of T p , the next step calculates the matrix elements of the real evaluation T r .These elements (Table 15) are obtained by multiplying matrix T p with normalized elements of the initial matrix (Table 13).Elements' normalization is completed by using Equation ( 28) (for C2) and Equation ( 29) (for C1, C3, C4 and C6).The real evaluation matrix element at the position t r32 is obtained using Equation (28)  The total gap matrix elements are calculated using Equation (30) and they are shown in Table 16.The total gap matrix element at the position g 32 is obtained using Equations ( 30) and (31).g 22 " t p32 ´tr32 " 0.288 ´0.00 " 0.288 (37) It is desirable that g ij has values that are close to zero ( g ij Ñ 0 ), because the alternative with the smallest difference between the (t pij ) and real evaluation (t rij ) is to be chosen.The best ranking alternative, the one with the lowest gap value, is the location-L 1, as it is shown in Table 17.The best ranking alternative, the one with the lowest gap value, is the location-L 1.

Sensitivity Analysis
Sensitivity analysis is recommended in most of the MCDA processes as an instrument for testing the results' stability against the subjectivity of the DM [66].Sensitivity analysis is completed by changing the initial criteria weights [32], as it is shown in Table 18.According to the recommendations [32], the sensitivity analysis of the location rankings has been conducted through five scenarios.The first scenario (A) deals with the equalized criteria weights w i = 1/5 = 0.2000.The second scenario (B) favours criteria C 1 and C 2 , w 1 = w 2 = 0.3000, while other criteria weights were evenly distributed, w 3 = w 4 = w 6 = 0.1333.The remaining three scenarios (C, D and E) favours one of the criteria, C 3 , C 4 and C 6, respectively.The favoured criteria has a weight w = 0.3500, while the others have a value of 0.1625.
The resulting rankings are given in Table 19.Results show that assigning different weights (priorities) to the criteria leads to different rankings, i.e., the model is sensitive to these weights.Comparing the first ranked alternatives in all scenarios, it can be concluded that the alternative L1 remained first in three out of five scenarios.In scenario C, the alternative L1 is in third place and in scenario E, it is in sixth place.Thus, the alternative L1 is superior to the others and the ranking presented in Table 17 is acceptable.
MAIRCA presents a new MCDA method so it should be compared with other already proven methods.There are different studies associated with location selection decisions that have been commonly carried out by using multi-criteria decision-making techniques, location problems with MOORA and COPRAS method [67][68][69], selecting distribution centers for a firm and location choice of distribution centers with VIKOR [70], location selection of logistics centres based on Fuzzy AHP and TOPSIS [71], and selection of logistics centres' locations with fuzzy TOPSIS based on entropy weight [72].The methods VIKOR, TOPSIS, MOORA and COPRAS are proposed in the previously stated references as prominent for location ranking.The above mentioned methods are compared with the MAIRCA method in order to compare location ranking.
Table 20 shows the results of alternatives' ranking according to criteria from Table 13 using the methods VIKOR, TOPSIS, MOORA i COPRAS.Comparing results obtained using the MAIRCA method (Table 17) with the results obtained from the methods VIKOR, TOPSIS, MOORA and COPRAS (Table 20), it can be seen that the alternative ranking for MAIRCA and VIKOR is similar.They are favoring alternatives L 1 and L 7. On the other hand, other methods suggest quite different rankings favoring alternatives L 7 and L 5.
MAIRCA and VIKOR are similar due to the same method of data normalization, defined by the Equations ( 28) and (29).In contrast, MOORA and TOPSIS use vector normalization, given in the Equation (38).
where r ij denotes normalized elements of the matrix, x ij presents the elements of the initial matrix and m is the total number of alternatives.The COPRAS method uses simple normalization, benefit type, Equation (39), r ij " x ij x ìj (39) For the cost type criteria, Equation (40), x íj (40) where x ij presents the elements of the initial matrix, and x ´ij and x + ij are the minimum/maximum criteria values.
In order to compare alternatives with different methods in an adequate manner, linear normalization of initial data is used for all methods, keeping the original mathematic calculations.Obtained data are presented in Table 21.It can be concluded from Table 21 that the rankings for methods TOPSIS, COPRAS and MOORA are the same like for the MAIRCA, while it is slightly different for the VIKOR method.The alternative L1 is ranked first in all methods so it can be concluded that this alternative is superior and the ranking suggested in Table 17 is confirmed and found to be credible.

Conclusions
The selection of the optimum AD location using MCDA in GIS was for the first time conducted in this research for the Carpathian region in the Republic of Serbia.It represents a model that can be used in other areas with similar natural characteristics.
Multi-criteria techniques DEMATEL, ANP (DANP) and WLC inside GIS software are suggested in order to select suitable areas for AD locations.DANP is used for criteria weighting, and WLC is used for the weighting sums and final identification of suitable AD locations.The use of these methods is enabled by the ArcGis software solution ArcGIS Advanced 10.2, from ESRI.The implementation of the GIS-MCDA model was proven successful and justified, because, on the basis of nine restrictions and six evaluation criteria, suitable sites for the AD installation could be determined.
The model has identified a total 44.5 km 2 of the concerned region as the most suitable area for the AD location.Finally, eight suitable locations were found by restraining areas smaller than 0.5 km 2 .The MAIRCA method found location L1 to be the most appropriate.
This model expands the frame of knowledge for FLP.The existing problem is considered using a new methodology, thus creating the basis for further theoretical and practical advancement.Moreover, the methodology used in this study highlights new criteria never considered before, but that were found to be crucial for the given area of research.This procedure allows for the inclusion of other criteria, which are not considered here, such as distance to watercourses and streams, distance to geological water courses or aquifers, distance to protected areas, areas which have a propensity to suffer from earthquakes or at high risk of floods.These criteria and restrictions are not set out by military and state regulations.For example, slope is a crucial restriction here but elsewhere it could be treated differently.The addition or omission of certain restrictions and criteria can be designated in new regulations, rules and guidelines because, currently, Serbian law does not fully recognize criteria restrictions for hazardous material issues.All of these combined factors are important in the selection of suitable locations.
It should be mentioned that this study is not intended as an exhaustive overview but is intended to point future research in the right direction.Nevertheless, this study includes experts' opinions and gives recommendations which are certainly valuable for future research in this area.

Figure 1 .
Figure 1.Hierarchical model to identify the location of the ammunition depot.

Figure 1 .
Figure 1.Hierarchical model to identify the location of the ammunition depot.

Figure 2 .
Figure 2. The process of the DEMATEL method.

Figure 2 .
Figure 2. The process of the DEMATEL method.

Figure 3 .
Figure 3.The geographical position of the studied area.

Figure 3 .
Figure 3.The geographical position of the studied area.

Figure 5 .
Figure 5. Cause and effect relationship diagram.

Figure 5 .
Figure 5. Cause and effect relationship diagram.
t 11 t 12 ... t 1n t 21 t 22 ... t 2n ... ... ... ... t n1 t n2 ... t nn (18)calculation of weighted normalized supermatrix W α (Step 4) is conducted by multiplying elements of unweighted supermatrix W with corresponding elements of the normalized total influence matrix T α D .Elements of the normalized total influence matrix.T α D were obtained by normalization of the total influence matrix T D , according to the Equations (17) and(18).
The elements of new weighted supermatrix W α were calculated after determining the elements of the matrix T α D .The elements of matrix W α are obtained by multiplying normalized total influence matrix of the dimensions T α (19){d i .D with the unweighted supermatrix W, Equation(19) P A 1 w 1 P A 1 w 2 ... P A 1 w n P A 2 w 1 P A 2 w 2 P A 2 w n g 11 g 12 ... g 1n g 21 g 22 ...

Table 1 .
Description of restrictions criteria.

Table 2 .
Description of evaluation criteria.
The land use is considered as the environmental protection factor.Land use data were taken from the CORINE Land Cover 2006 (CLC2006) database.CORINE nomenclature has five main categories but three of them are excluded from the analysis.

Table 3 .
Data sources for input layers in land suitability assessment for the Carpathian region.

Table 4 .
Evaluation criteria description and standardization.
*: Level 3 of the CORINE nomenclature; **: Values are given for the mountain relief.3.4.3.The Comparison and Relative Assessment of the Criteria Weights (DEMATEL-ANP)

Table 6 .
The average matrix Z.

Table 7 .
The normalized initial direct-relation matrix D.

Table 8 .
The total relation matrix T.

Table 9 .
The sums of matrix T values for six criteria.
* Values greater than the threshold value.

Table 8 .
The total relation matrix T.

Table 9 .
The sums of matrix T values for six criteria.

Table 12 .
Weighted coefficients of the criteria.

Table 13 .
Initial matrix for the sustainable locations evaluation.

Table 13 .
Initial matrix for the sustainable locations evaluation.

Table 14 .
Matrix of the calculated theoretical evaluation T p . .

Table 15 .
Matrix of the calculated real evaluation T p .

Table 16 .
The total gap matrix.

Table 17 .
Alternative ranking according to MAIRCA.

Table 18 .
Combination of weights for the sensitivity analysis.

Table 19 .
Alternative rankings for different scenarios.