A multiobjective GIS-based land use planning algorithm

Computers, Environment and Urban Systems xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Computers, Environment and Urban Systems journal homepage: www.elsevier.com/locate/compenvurbsys

A multiobjective GIS-based land use planning algorithm Theodor J. Stewart a,b,⇑, Ron Janssen c a

Department of Statistical Sciences, University of Cape Town, Rondebosch 7701, South Africa Manchester Business School, The University of Manchester, Manchester M15 6PB, UK c Institute for Environmental Studies (IVM), Vrije Universiteit, De Boelelaan 1087, 1081 HV Amsterdam, The Netherlands b

a r t i c l e

i n f o

Article history: Received 15 November 2013 Received in revised form 3 April 2014 Accepted 3 April 2014 Available online xxxx Keywords: Multiple criteria decision making Reference point approach Genetic algorithm Land use planning Spatial optimization

a b s t r a c t This paper purposes an enhanced land use optimization model for land-use planning with a new spatial component. This component uses a simple representation of the proximity of related land uses to each other as a function of distances between parcel centroids. A special purpose genetic algorithm is developed for solving the resultant optimization problems for both the direct (additive) objectives and the indirect (spatial) objective. The context relates to interactive decision support for land use planning in which the data are stored in a vector-based GIS, and the requirement was to integrate the multiobjective optimization with the GIS structure. The present work thus extends earlier work by the authors which used a grid (raster) structure. The model is based on a reference point approach in which both additive and spatial goals can be specified. Numerical testing of the algorithm, and experimentation with possible user inputs, are described in the context of a real case study from a region of The Netherlands. It is shown that the simplified spatial proximity measure and the associated algorithm produce consistent results in which the spatial distribution of activities are essentially the same as with more complex modeling of spatial goals, achievable in the particular case study with little loss in terms of the additive objectives. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction During the last decade a large number of spatial decision support systems have been developed to assist decision makers in the field of resource allocation and, in particular, spatial planning issues (Geertman & Stillwell, 2009). When alternatives or objectives are spatial, data are needed on the geographical locations of alternatives, spatial formulations of objectives and data on the spatial pattern of criterion values. This requires a combination of multicriteria methods with a geographical information system (GIS) (Arciniegas, Janssen, & Omtzigt, 2011; Malczewski, 2010). This combination is referred to as a spatial decision support system (SDSS) (Carver, 1991). GIS is used to produce thematic maps and to perform spatial operations. Multicriteria methods are used to translate these maps into value maps, optimal or compromise maps and to rank spatial alternatives (Arciniegas, Janssen, & Rietveld, 2013; Alexander et al., 2012; Janssen, Arciniegas, & Verhoeven, 2013). Spatial multicriteria analysis (MCA) typically starts with a set of land-use alternatives that is defined beforehand. This set of ⇑ Corresponding author at: Department of Statistical Sciences, University of Cape Town, Rondebosch 7701, South Africa. Tel.: +27 21 6503224. E-mail addresses: [email protected] (T.J. Stewart), [email protected] (R. Janssen).

alternatives can be defined using the inputs from experts such as spatial planners or landscape architects (Arciniegas & Janssen, 2012), land-use models, (e.g. Lau & Kam, 2005; Oxley, McIntosh, Mulligan, Winder, & Engelen, 2004), or by applying design methods based on multiobjective optimization techniques (MOOT). Such design techniques generate an ‘optimal’ solution for a specific preference structure from a large or possibly infinite set of alternatives, where the set of alternatives to choose from is implicitly defined through constraints and exogenous influences. In other words, the optimal solution is created or ‘designed’ by the SDSS using techniques based on tools such as multiobjective linear programming (Aerts, Eisinger, Heuvelink, & Stewart, 2003; Cova & Church, 2000). As a special case of design methods, interactive optimization offers solutions to the planner in a number of steps where, after each step, the planner can change the conditions for optimization (Janssen, van Herwijnen, Stewart, & Aerts, 2008; Stewart, Janssen, & van Herwijnen, 2004). Using optimization to generate the optimal solution requires that all objectives can be described in mathematical terms and incorporated in the optimization model. Geographical information systems make use of two types of data: spatial data and attribute data (Longley, Goodchild, Maguire, & Rhind, 2005). The spatial data describe location and shape of spatial entities; the attributes are the properties of these

http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002 0198-9715/Ó 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

2

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

spatial entities. In using GIS as input to an optimization model, the decision units are related to spatial entities, in the sense that choices have to be made for land uses at every spatial location. An attribute table includes both the decision variables, e.g. types of land use, together with any attribute values which need to be included in the objective function(s). Spatial data in a GIS are typically arranged using one of two models: vector or raster. Entities in vector format are represented by strings of coordinates (points). Two points can be connected to form a line segment, while sequences of lines can be connected to each other, terminating back at the starting point to form a polygon (parcel or area). Attribute data are stored for each polygon (which can be of varying sizes). Data in a raster model are stored in a two-dimensional matrix of uniform cells on a regular grid. Depending on the model used, each polygon or grid cell is assumed to have homogeneous properties. By their nature raster data are substantially easier to include in mathematical representations of the world for purposes of optimization. As a result most GIS-based applications of MOOT use gridbased data as their input. Using a grid based representation of a planning region, Stewart et al. (2004) and Janssen et al. (2008) showed that it was possible to formulate a spatial planning problem in mathematical terms and apply MOOT to generate optimal solutions interactively. Unfortunately, they also had to conclude that using a grid size which would realistically describe the planning region leads to unrealistically long computation times, as a result of the large number of decision variables. This prevented use in a fully interactive setting where short response times are essential. Other examples of grid-based applications can be found in Cao, Huang, Wang, and Lin (2011), Ligmann-Zielinska, Church, and Jankowski (2008), and Santé-Riveira, Boullon-Magan, Crecente-Maseda, and Miranda-Barriós (2008). Looking at Fig. 1 it is immediately clear that moving to a vector based representation will lead to a more efficient representation of the problem. The vector presentation (Fig. 1(b)) shows the border of the spatial decision unit (e.g. a parcel) with type of land use as its decision variable. The raster representation of the same spatial unit requires a large number of cells with an equally large number of decision variables bound together by the constraint that the land use for all grid cells must be the same. Although it is clear that a vector representation leads to a more efficient representation of the problem, it is also clear that the switch from grid to raster creates new complications. In a raster each grid cell has the same shape and size, borders on exactly four other grid cells and has four borders of equal length. In a vector format each polygon can have any shape and size, and have any number of borders of various lengths. Closely related land use optimization problems based on a GIS link include Matthews, Sibbald, and Craw (1999), who do not however consider the spatial objectives which we shall later discuss, Demetriou, See, and Stillwell (2013), who examine

partitioning of a region into distinct land uses, and Porta et al. (2013) which is perhaps most closely allied to our problem, but with some differences we shall highlight later. Janssen et al. (2008) differentiate between two types of objectives, namely simple additive objectives, which associate costs and/or benefits with the allocation of any particular land use to a specific cell, which are then cumulated additively across all cells; and spatial objectives which indicate the extent to which the different land uses are contiguous (i.e., the extent to which activities are or are not fragmented across the region). The shift from a grid to vector based representation do not create great difficulties for the additive objectives. The differences in area size of the decision units can easily be accommodated using area size as a weighting factor in calculating overall performance. The shift from raster to vector is not so easily implemented for the spatial objectives, as we shall discuss in the next section. The present paper describes a genetic algorithm that can be used to generate land use plans that maximize both additive and spatial objectives in a vector based GIS environment. The algorithm will be applied as part of a series of collaborative planning workshops to support a land use allocation problem in a peat-meadow polder in The Netherlands. Results from the algorithm will be used to generate a number of land use plans that can serve as reference plans at the start of these workshops. This article is organized as follows. In the next section, alternative formulations of the land use planning problem are discussed, giving rise in general to a non-linear combinatorial optimization problem. In Section 3 we describe a practical case study providing background for the type of situation in which our algorithms would be applied. Algorithmic methods for solving the non-linear problem are described in Section 4, including reference to a simplified linear formulation excluding spatial objectives. The specially designed genetic algorithm (GA) for the full optimization step is then defined in some detail, and numerical testing of this GA is reported in Section 5. Finally, in Section 6, some experimental runs for the case study are described. 2. Mathematical formulation Suppose that the total region is represented by I parcels of land, labeled i ¼ 1; . . . ; I. In a raster-based GIS, these would be rectangles of uniform size, but for a vector-based GIS each parcel would be represented by ‘‘polygons’’ of unique size and shape. The underlying assumption is that parcels are defined to a level of size and resolution such that a single land use would need to be assigned to each parcel. Parcels for which land uses are fixed (e.g. streams, roads, etc.) can either be omitted from the data base, or the land use assignment can be pre-specified as a hard constraint. (Our model implementation allows the user both options.) This may be a restriction in other application contexts than our own, where division of a polygon into more than one land use may be desired.

Fig. 1. Raster (a) and vector (b) data representation.

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

Let the number of possible land uses be M, labeled m ¼ 1; . . . ; M. Each parcel is assigned to one and only one land use. A complete solution is thus defined by I  M binary variables, say dim taking on the value 1 if parcel i is assigned to land use m (and 0 otherwise). We note that this problem context is closely allied to that discussed in Matthews et al. (1999) and Porta et al. (2013), but quite different to the very interesting alternative problem of land partitioning (Demetriou et al., 2013). For notational convenience we shall denote the full two dimensional array of these decision variables, i.e. a full allocation, by d. The requirement that all parcels be assigned one and only one land use imposes the constraints: M X

dim ¼ 1 for i ¼ 1; . . . ; I:

ð1Þ

m¼1

Typically, the allocation would be subject to constraints on the total area allocated to each land use, expressed in general form by:

amL 6

I X ai dim 6 amU

for m ¼ 1; . . . ; M:

ð2Þ

i¼1

where ai is the area of parcel i, and amL and amU are respectively lower and upper bounds on total area to be allocated to land use m. We recognize both direct and spatial management objectives. Direct objectives relate to (for example) economic, social and environmental concerns, such as value from agriculture, from recreation and from nature. Similar structures of objectives appear in Matthews et al. (1999). As illustrated by the case study of Section 3, in most cases such values would be assessed by stakeholders or experts as in standard value measurement theory. For purposes of this paper, we shall assume that such values to be additive across parcels (see discussion in Belton & Stewart (2002) Chapter 4 for a theoretical background to this assumption), and we thus term these as additive objectives. Let bimk be a standardized measure of the value in terms of objective k (k ¼ 1; . . . ; K) realized by allocating land use m to parcel i. Typically we might define such measures so that bimk ¼ 1 for the ideal land-use/parcel combination in terms of objective k, falling to bimk ¼ 0 for a land-use/parcel combination that makes no contribution to objective k. Overall performance on objective k implied by the allocation d; V k ðdÞ is then defined by an area-weighted average across all polygons, namely:

PM PI V k ðdÞ ¼

m¼1

i¼1 ai bimk dim

PI

i¼1 ai

ð3Þ

In the simpler raster-based formulation each parcel is a grid cell, where all cells have the same area, and the above expression can be simplified. But for vector-based representation, the sizes of each polygon can vary enormously. As an aside, we comment that in practice it may sometimes be useful to structure one or more of the additive objectives in terms of a set of sub-objectives. Suppose that objective k can be disaggregated into Lk subobjectives, and let bimk‘ be the value in terms of subobjective ‘ of k from allocation of land use m to parcel i. Under PLk mild conditions we then be able to define bimk ¼ ‘¼1 wk‘ bimk‘ , where wk‘ is the relative importance of subobjective ‘ to objective k. In the software described below, we did allow for such disaggregation into subobjectives, but this option plays no role in the further description of the optimizing algorithms. The spatial objective relates to the need to ensure that each land use is not unduly fragmented, but has adequate contiguity. Spatial objectives have been recognized by a number of authors (e.g. Aerts et al., 2003; Chen, Li, Liu, & Liu, 2010; Santé-Riveira et al., 2008; Vanegas, Cattrysse, & Van Orshoven, 2011). Recent work relating to such issues include Porta et al. (2013), who however deal with a single aggregated compactness objective over all uses, and Demetriou et al. (2013) who address the difficult, but

3

rather different, problem of the shape of parcels of land generated by partitioning. In previous work (Stewart et al., 2004) we introduced a number of performance measures related to the spatial objectives, derived from initially quite qualitative expressions of goals. These measures were based on a clustering approach: one cluster is a set of parcels sharing the same land use, such that all parcels in the set are interconnected (i.e. each is contiguous to another parcel in the set). Three such measures were then defined as follows (modified slightly from the original reference to take account of the possibly varying sizes and shapes of the polygonal parcels):  The square root of the number of disjoint clusters for each land use (to be minimized); the use of the square root operator is a variation from Stewart et al. (2004), introduced to recognize that marginal increases from a low base are relatively much more influential than from a larger base.  Ratio of the area of largest contiguous cluster to total area allocated to the land use (to be maximized for each land use, to reflect a preference that as much as possible of the land use should be in a single cluster).  A shape index defined as the area weighted average of a compactness measure for clusters involving this land use, where for a single cluster this measure is defined by perimeter divided by the square root of area. This shape index should also be minimized. Implementation of a clustering algorithm for a grid based structure is relatively straightforward. Two cells are contiguous if they follow one another on the same row, column or diagonal, while the perimeter of a cluster is directly calculated from the number of edges of cells shared with cells of other land uses. For a vector-based, structure, however, a number of practical problems emerge, for example:  The search for contiguous polygons involves checking all the border lines of the polygons.  The perimeter needs to calculated by the sum of all perimeters of polygons in the cluster less the lengths of common border lines.  (Perhaps the most critical problem for practical implementation), the standard GIS structure contains non-usable parcels such as ditches and streams between parcels, but these parcels should still be viewed as otherwise ‘‘contiguous’’. In spite of these difficulties, a version of the genetic algorithm from Stewart et al. (2004) was implemented for the vector based structure. This was primarily for purposes of validity checks on the final algorithm developed in Section 4.2, and discussed in the context of the numerical results of Section 5. The complications of interfacing with the GIS data base to address the above difficulties appeared not to make this approach suitable for the envisaged decision support system. (The mode of use of the algorithm is to extract relevant GIS data into separate file to be read by the optimizer as a stand-alone executable file; special GIS skills were needed to design a region-specific macro which would supply the necessary data to the original optimizer.) A simpler representation of the spatial objectives was thus needed, where the necessary data can be extracted from the GIS with no special pre-processing. One approach may be to ignore the spatial objectives entirely, in which case the problem is a simple allocation problem which can in principle be solved by integer linear programming, although in Section 4.1 a further simplification is described, allowing use of continuous linear programming. The primary contribution of this paper, however, is to retain consideration of the spatial objectives, but in a simpler form. The

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

4

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

coordinates of the centroids (i.e. the ‘‘centre of gravity’’) for each polygon are provided directly by the GIS software (without need for any intermediate processing). Let us denote the centroid for polygon i by ðxi ; yi Þ. A Euclidean distance between any pair of polygons, say i and j is then given by:

Dij ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi  xj Þ2 þ ðyi  yj Þ2 :

ð4Þ

We then define two polygons i and j to be proximate (in the proximity of each other) if Dij 6 q for a specified run-time parameter q. This is illustrated in Fig. 2, in which are shown four illustrative polygons and their centroids. (Of course, the entire region is in fact filled by polygons.) In the figure, polygons A and B are contiguous (having a shared border) as well as proximate; polygons A and C are not contiguous, but are close enough (although close to the limit) to be proximate for the specified value of q; polygons A and D are not proximate. We have observed that polygons may be ‘‘proximate’’ even if not sharing a border, which may seem to conflict with our basic spatial objectives. However, we shall be demonstrating in the numerical results reported in Section 5 that maximizing overall proximity of each land use does tend to cluster polygons assigned to this land use in a largely contiguous manner. In our implementation q was defined by a given quantile (specified as a model parameter) of the distribution of all Dij values in the data set (although the distribution was estimated by sampling rather than by a full enumeration). Note that definition of the proximate relationship does not use precise contiguity directly, and is thus not influenced by intervening non-use areas such as canals or roads. Then for a given solution (represented by d), we can identify for each polygon i which of the polygons proximate to it are assigned the same land use as for i. We then define pi as the area weighted proportion of all polygons proximate to i which are assigned the same land use. In other words, if we let uðiÞ be the land use allocated to polygon i in the solution given by d, then pi is given by:

P

pi ¼

fj:Dij 6q;uðjÞ¼uðiÞg aj

P

j:Dij 6q aj

:

ð5Þ

Clearly, pi ¼ 1 if all polygons proximate to i share the same land use assignment as i, while pi ¼ 0 if i is an isolated polygon with land use different to all others in its proximity. Finally, for each land use m, we can define the degree of contiguity by the areaweighted average of the pi for those polygons assigned to land use m, i.e. by:

P

i:d ¼1 P m ¼ P im

ai pi

i:dim ¼1 ai

:

ð6Þ

Each P m represents a further objective. As users may have different priorities for the degree of contiguity required for different land uses, we shall retain the P m as separate objectives in the multiobjective formulation below. For each land use m, a maximum value of P m ¼ 1 is achieved if all polygons assigned this land use are proximate to each other, while P m ! 0 if this land use is widely dispersed across the region. In summary, therefore, we have a multiobjective optimization, in which the I  M binary variables dim are to be selected in order to maximize a total of K þ M objectives, namely the K additive objectives which we have denoted by V k ðdÞ and the M objectives denoted by P m , subject to the constraints (1) and (2). Additional constraints can be added, either to fix land use in certain polygons, or to prohibit certain land uses in particular polygons. Such constraints are easily incorporated into any of the optimization models described below, and for this reason we will not make further explicit reference to them. One approach to multiobjective optimization that has gained much prominence is Evolutionary Multiobjective Optimization (EMO) (Deb, 2001). This approach strives to provide the user with a visual representation of the efficient frontier, i.e. of the set of Pareto optimal solutions, not dominated by any others, from which a direct choice may be made. Such a visual approach is less useful for much more than 2 or 3 objectives. Our formulation can easily result in 10 or more objectives (e.g. if there were 4 additive objectives and 6 land uses). We thus selected to use an interactive method in which objectives can be aggregated in some way, and in which the user can experiment with different parameter values in the aggregation. The use of simple weighted sums of objectives in non-linear multiobjective optimization presents many pitfalls, as explained for example in the last paragraph of Section 3.1.4 of Miettinen (1999). As in Stewart et al. (2004), we made use of a form of reference point approach method as an alternative (see Miettinen, Ruiz, & Wierzbicki (2008, chapter 2) for a discussion of reference point approaches). One important feature of reference point methods is that the solutions obtained are always Pareto optimal (efficient). In fact, by experimenting with different reference points as defined below, all points on the efficient frontier can in principle be reached, in contrast to the use of weighted sums in which for non-linear problems, some Pareto optimal points may not be reachable. In our implementation, based on the experience reported in Stewart et al. (2004), we adopted a slightly modified (smoother) form of scalarizing function from the conventional form described in Miettinen et al. (2008, chapter 2). With this modification, the aggregation procedure may then be described as follows: let z1 ; z2 ; . . . ; zR (where R ¼ K þ M) be the collection of all performance measures, representing both the additive and spatial objectives. Let zr be the ideal, i.e. the maximum value for zr , and let g r be an aspiration level or performance goal for objective r. Our scalarizing function is then defined by:

SðdÞ ¼

R X r¼1

Fig. 2. Illustration of proximity relationships between four polygons.

 Wr

zr  zr zr  g r

c ð7Þ

for some c > 1, and where W r is a relative measure of the importance of objective r. In previous experience, we have found c ¼ 4 to produce useful results, which retains the strong emphasis on minimizing the largest deviation from the conventional scalarising function (which is based on an augmented Chebychev norm), while not allowing a single large deviation to completely dominate the process. For purposes of the present paper, which was to assess the numerical behavior of our approach, fairly arbitrary values were chose for the W r and g r . In follow-up research, we will experiment with modes of guiding users through experimentation with W r and

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

g r , although experience with the earlier model suggests that results are more sensitive to choice of g r than to W r . For any given set of aspiration levels, the problem becomes that of minimizing SðdÞ subject to (1) and (2) (as well as any constraints fixing or prohibiting specific land uses in some polygons). The intention in the use of the algorithm is to allow the aspiration levels to be adjusted in response to user preferences expressed in an interactive manner.

3. Illustrative case study Before proceeding to details of the algorithm which was implemented, we illustrate the problem structure we have discussed by means of a real case study. In the last few years we have used spatial tools in collaborative planning workshops to support a land use allocation problem in a peat-meadow polder in The Netherlands. Specific tools were integrated into the design of each workshop and implemented in an interactive mapping device called the touch table. Within these workshops it is accepted that the problem cannot be formulated completely in mathematical terms. The end product of the workshops was a negotiated land use plan for the polder (Arciniegas & Janssen, 2012). The tools used focused on evaluation of plans generated by the participants within the workshops. In the later numerical studies, we will assess how the algorithms described in the next section can be used to generate land use plans for this polder. It is our intention to integrate the genetic algorithm in our workshop design to be able to generate reference plans as starting points for negotiation and as reference for the performance of the negotiated plans, as well as to generate re-optimized plans in the light of discussions. The Bodegraven polder is a multifunctional low-lying peat meadow area of 4072 hectares, located in the province of South-Holland, The Netherlands. Historically, water tables have been artificially controlled in the polder to enable several functions, primarily commercial dairy farming. Despite still being used widely for intensive agriculture, Bodegraven is also important for their natural, cultural and historical values. Main driving forces in the area are ground water level and land use. A high ground water level will reduce soil subsidence and is generally good for nature. For profitable agriculture a low water level is required. Land use in the polder will need to be changed to accommodate the various uses in combination with the predicted ground water levels. Multiple stakeholders such as the local water board, the city of Bodegraven, the province of South Holland, farmers organizations, and nature conservation organizations, as well as individual farmers, residents, and recreational visitors are involved. The provincial authorities have started a planning process to change both water management and land use in the area. Water management is the driving force within this process and land use has to adapt to changing water conditions (Strategiegroep Gouwe Wiericke, 2007). This implies that agreement is needed on both the watermanagement strategy and the consequent land use. At the start of the process the terms of reference were defined in a framing workshop. This involved an inventory of the objectives of the various stakeholders and the evaluation criteria linked to these objectives (Cornelisse, Janssen, Arciniegas, & Omtzigt, 2007; Strategiegroep Gouwe Wiericke, 2010). Three key objectives were identified: (1) profitability of intensive agriculture, (2) maximization of the visual quality of the landscape and (3) maximization of the natural value of the area. The objectives visual quality of landscape and natural value included several sub-objectives. Visual quality was defined as the weighted sum of landscape perception, cultural historic value and recreational value. Natural value was defined as the weighted sum of the qualities of meadow birds, species-rich grasslands and marsh

5

birds (Janssen et al., 2013). Each combination of land use and water level is valued for all sub-objective by experts according to their field of expertise (see Arciniegas & Janssen (2012) for details on this valuation approach). We distinguish three types of land use: 1. Intensive agriculture; 2. Extensive agriculture; and 3. Nature. Current land use of the polder is shown in Fig. 3. Each polygon on the map represents a parcel in the field. This information was used to generate the attribute table used as input for the genetic algorithm. Input data for the simple optimizer model as well as for the genetic algorithm model the genetic algorithm are stored in an attribute table. This table includes values for all objectives and criteria for each parcel and for each land use (i.e. for the bimk parameters). As part of the reallocation of land use in the polder, the provincial authorities decided to allocate 860 ha of nature and 1600 ha of land for extensive agriculture. This means that agricultural land must be purchased for conversion to nature and subsidies must be made available to support the transition to extensive agriculture. This information is used to set constraints for the optimization process. 4. Solution methods 4.1. Simplified models As mentioned above, the primary goal of the research reported here was to develop a comprehensive optimization model to address both additive and spatial objectives in an efficient manner, so that the software could be integrated into interactive workshops. Nevertheless, it was worth investigating whether there was any advantage in providing a simpler, less rigorous optimization that could very rapidly generate tentative solutions for discussion and provide performance bounds on the additive objectives, by ignoring the spatial objectives completely. This reduces the problem to maximization of the K linear objectives given by (3) subject to (1) and (2). Although we have argued above that a linear weighted sum of the V k ðdÞ would not be justifiable for a complete solution, for the purposes of a simple module with the above-stated purposes such a weighted sum may be useful, especially if linked to sensitivity analysis on the weights. In other words, the simplified model would be based on maximization of:

UðdÞ ¼

K K M X I X X X ai bimk W k V k ðdÞ ¼ Wk dim A m¼1 i¼1 k¼1 k¼1

which can be written as:

UðdÞ ¼

M X I X kim dim

ð8Þ

m¼1 i¼1

where we define:

kim ¼

K X W k ai bimk k¼1

A

:

ð9Þ

The decision variables remain binary, and in a typical problem we could well expect 2000 or more polygons and perhaps 5 or more additive objectives, resulting in 10,000 or more binary variables, which would be a large number for achieving a rapid ‘‘quick-anddirty’’ solution so that further simplifications were sought. 4.1.1. Formulation as an assignment problem Suppose we replace the constraints (2) by the following: I X dim ¼ mm

for m ¼ 1; . . . ; M

ð10Þ

i¼1

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

6

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

Fig. 3. Current land use of the Bodegraven polder.

for specified integers mm . Maximization of (8) subject to (1) and (10) is a classical assignment problem. Provided that the coefficients kim are integer (which can be ensured to any desired level of approximation by appropriate re-scaling), the problem is totally unimodular, and the solution will automatically be binary without the need to specify this requirement as a constraint. This simple LP problem is rapidly solved. The only practical difficulty is in selecting the mm in order to achieve satisfaction of the area constraints (2). For implementation of the simple scheme, it is convenient to assume that the constraints (2) are nearly equalities within a small tolerance, i.e. that amL ¼ ð1  sÞam and amU ¼ ð1 þ sÞam , where the am are desired areas for each land use, and s is a small positive tolerance parameter. With this assumption, suitable values for the mm are obtained in a simple iterative scheme which is in outline as follows:  Start with an initial guess at values for mm , e.g. by dividing the total number of polygons in proportion to the am .  Solve the assignment problem. Let Am be the resultant areas allocated to each land use. If amL 6 Am 6 amU for all m, then stop.  Otherwise, adjust the mm up or down, dependent on the ratio of Am to am , and repeat from the previous step. We have found that no more than about 6 or 7 iterations are typically required to reach feasibility. The overall computational time to perform the optimization was around 1 s on a Dell Vostro Notebook computer. This is effectively instantaneous in any interactive session. However, as expected, and as will be illustrated further in Section 6, the solutions are spatially more fragmented than would be desirable. Nevertheless, the solution provides a good benchmark for achievements possible for each additive objective. 4.1.2. Ignoring integer constraints The above is the simplified version we have currently been using to obtain results reported in Section 6. However, there does

seem to be some promise in using an alternative simplification. This is based on maximizing (8) subject to (1) and (2), but ignoring the binary constraints on the dim . Non-negativity together with (1) ensures dim 6 1, but some values may turn out to be fractional. Preliminary experimentation has suggested that no more than one or two of the many thousands of dim variables are non-integer. We have not been able to prove this result rigorously, but if this turns out to be a general principle, then this would be a useful simplification. Bounds on objective function achievement are likely to be quite tight, while simple rounding off of the solutions (taking into consideration the need to have an assignment for every polygon) is likely to be an excellent starting point for further interactions. However, there is little practical advantageous in terms of computational time. 4.2. Comprehensive optimization model The primary challenge remains that of generating good solutions in terms of both additive and spatial objectives in an efficient and effective manner. As indicated earlier, the approach is to minimize (7), subject to (1) and (2), other constraints on land uses that can be allocated to each polygon and to all dim being binary. The scalarizing function (7) is strongly dependent on choice of the aspiration levels g r , values of which are determined by interaction with users. We shall first discuss the optimization algorithm for a fixed set of the g r , after which we shall return to discussion of the manner in which these aspiration levels are set. The optimization problem is similar in structure to that addressed in our earlier work (Stewart et al., 2004) for the simpler context of a grid structure. There a genetic algorithm (GA) approach had worked well, and for this reason the same approach was adopted here. As discussed in Stewart et al. (2004), for algorithmic efficiency (especially for interactive decision support), there is considerable advantage in specially designing the GA to exploit the special features of the problem. A GA is broadly defined

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

by a sequence of three processes applied iteratively, namely (i) generation of an initial population of solutions; (ii) the crossover step to generate child solutions; and (iii) gene mutation of the child solutions. 4.2.1. Initial generation of solutions This is perhaps the most critical point at which special features need to be exploited, as simple random selection of the binary variables will result in a considerable waste of computational time in processing infeasible solutions and highly fragmented solutions (i.e. very poor performance in terms of the spatial objectives). The essential aspects of the process adopted may briefly be described as follows. (Technical details of precise numerical values for parameters of the algorithm are suppressed in interests of brevity, but are available from the first author on request.). 1. Allocate any land use that is fixed by the constraints (e.g. where only one land use has been designated as suitable for a particular polygon). 2. If all polygons have been assigned a land use, then the process terminates. Otherwise, select one as-yet unallocated polygon at random. Assign an allowable land use to this polygon randomly, with probabilities proportional to the PK additive objective performance (defined by k¼1 W k bimk ) for the chosen polygon i. These probabilities are further adjusted as the algorithm progresses by a factor related to how much area has already been assigned to the this land use (relative to the required area). 3. The selected polygon becomes the seed for a new cluster of parcels with common land use. An ‘‘adjacency list’’ defined by those polygons which are as-yet unalloacted, proximate to the selected polygon and suitable for the chosen land use. If this list is empty, return to step 2; otherwise continue. 4. Randomly select a polygon from the adjacency list and allocate the same land use to it. 5. Terminate growth of the current cluster if maximum area for the chosen land use has been exceeded, or by random choice with a non-zero probability that increases as the maximum area is approached. If growth is terminated in this way, then return to step 2; otherwise continue to the next step. 6. Extend the adjacency lists by appending further suitable polygons not in the list, but proximate to the last polygon added to the cluster. Return to step 4. 4.2.2. Crossover step For each crossover step, a random selection of j solutions (where j is a model parameter) is made from the parent population, to form a parent tournament. These are sorted by values of the scalarizing function, but with a model specified probability of erroneously reversing the order for any pair of solutions compared. The top two ranking solutions are then identified as the parents. For polygons in which both parents have the same land use, this use is retained in the child solution. All other polygons are partitioned into sets defined by pairs of land uses: polygon i belongs to set Mmn if land use m is allocated in one parent and land use n in the other (where m < n, so that each pair occurs only once). Let Amn be the total area of the polygons in Mmn ; the aim will then be to assign an area of approximately Amn =2 to each of land uses m and n, in such a way that these uses are as contiguous as possible. To this end, for each ðm; nÞ, the following assignments are made. 1. If Mmn contains only one polygon, this polygon is assigned land use m or n with equal probabilities, and the process is complete.

7

2. Otherwise find the two polygons in Mmn which are maximally distant from each other. Randomly assign one of these to land use m. Order the remaining polygons in order of increasing distance from the selected polygon, start assigning land use m to these until the desired area Amn =2 is achieved. 3. Remaining polygons are assigned land use n. 4.2.3. Mutation step Two polygons are selected at random, and a list of proximate polygons to each is constructed, so that the total area of polygons in each list is as close as possible to a given fraction (a model parameter) of the total area of the region. Then the land uses of in each list are swopped. 5. Numerical testing of the genetic algorithm The genetic algorithm was based on a population size of 200, and on 200 children per generation. After generation of the 200 child solutions, the population for the next iteration was taken to be the best 200 out of the set of all 400 parent and child solutions (a simple ‘‘elitist’’ structure). These properties were chosen as they had worked well for the original cluster-based approach described in (Stewart et al., 2004). Numerical studies were then undertaken in order primarily to investigate the quality of solutions obtained using the simpler proximity objectives in the multiobjective optimization. Pure numerical convergence of the algorithm was less of an issue, and are not reported in detail here. The stopping rule for the GA used in the tests was based on the earlier experience from Stewart et al. (2004): GA iterations were terminated when the relative improvement per generation in the scalarizing function was not more than 0.5% in each of 20 consecutive generations. Note that the scalarizing function value is bounded below by 0, so that an absolute bound on possible improvements is available. This is not to say that further experiments on convergence may not have value, but as we shall see, there is little scope for further improvements in results obtained for practical planning purposes. The initial set of studies examined the effects of two key parameters defining the proximity measures and their implementation in the scalarizing function, namely:  The population quantile used to define q.  The goal level used in the scalarizing function for the proximity objectives – this is an important tuning parameter as numerical values for the proximity measure may not be naturally interpreted by users. Parameters related to the genetic algorithm itself were held fixed at this stage (informed by the earlier experiences), namely (i) the proportion of the total area to which the mutation step is applied (set at 0.03); (ii) the number of solutions included in the parent tournament (set at 5); and (iii) the probability of erroneous reversal of pairwise ranking in the tournament (set at 0.1). Table 1 gives computational results for each case considered. For each case, the full optimization was executed 5 times, and the averages are reported in the Table. The column headed ‘‘aggregate value’’ displays the (average) value for the aggregated additive objectives, while the three columns headed ‘‘spatial objectives’’ display the (average) values, also averaged across land uses, for the three original spatial objectives (numbers of clusters, relative size of largest cluster, and the shape index). These are recorded both as indication of the quality of solutions obtained, and as comparison (see below) with results of explicitly optimizing all three sets of spatial objectives. Recall that the number of clusters and the shape index need to be minimized and the relative size of the largest cluster

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

8

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

Table 1 Computational results for variations in proximity measures. Proximity quantile

Proximity goal

Time (s)

Aggregate value

0.04 0.04 0.04 0.02 0.02 0.02 0.06 0.06 0.06 0.03 0.03

0.5 0.6 0.7 0.5 0.6 0.7 0.5 0.6 0.7 0.6 0.65

93.4 102.0 107.8 87.2 98.0 83.0 107.8 125.2 89.8 80.6 92.6

0.524 0.522 0.520 0.525 0.525 0.524 0.522 0.521 0.521 0.522 0.522

maximized. The reported computational times were for implementation on a Dell Vostro computer, running Windows XP on a 2.27 GHz Intel processor. All but the last two rows in Table 1 represent a systematic study of the effects of the two parameters. The effects are generally small, especially in the light of variations between individual optimization runs. There is (perhaps not surprisingly) no real effect on the aggregate additive value. There is some indication that the higher value for the proximity goal reduced running times, but at the expense of slightly poorer spatial performance, while larger proximity quantiles increased running times. The last two rows in Table 1 report on some further finer tuning, leading to the decision to maintain a proximity quantile of 0.03 and a proximity goal of 0.65. It is interesting to record here that running a full cluster-based optimization on the same problem, placing an equal weighting on the three cluster-based components of the spatial objectives led to averages of 0.524 for the aggregative additive objectives, 36.5 for number of clusters per land use, 0.613 for the relative size of the largest cluster, and 12.2 for the shape index. The proximity based approach thus gives very similar results for performance in terms of the objectives, and in fact, visual inspection of the resulting maps reveals very little obvious difference. Interestingly, the proximity based approach does substantially better for the relative size of the largest cluster, although somewhat worse on the shape index. As the proximity-based model focusses on keeping similar activities close together, but not on the shape of clusters, this observed shift in tradeoffs between the relative size and shape objectives should not be unexpected. We conclude, therefore, that the proximity based measure is a valid and convenient tool for the practical purposes at hand. Finally we briefly examined the effects of the genetic algorithm tuning parameters, i.e. the proportion of the total area to which the mutation step is applied, the number of solutions included in the parent tournament and the probability of erroneous reversal of pairwise ranking in the tournament. Numerical results are displayed in Table 2 (in the same format as that of Table 1). To a large extent, performance is robust to the tuning parameters. The main effect is that of mutation. Recall that the purpose of

Spatial objectives No. clusters

Rel. size

Shape

44.1 35.9 29.5 42.2 40.9 36.0 41.1 34.7 36.9 42.2 35.8

0.824 0.848 0.902 0.784 0.806 0.845 0.811 0.824 0.892 0.840 0.851

16.0 13.7 12.9 15.4 15.0 12.9 15.2 14.2 13.6 14.9 13.7

mutation is to maintain diversity in the population to avoid premature convergence. If the mutation step is omitted, then running time is (not surprisingly) drastically reduced. Although lower values of the scalarizing function are achieved, it is clear from the Table that the deterioration the spatial objectives is small (although in subsequent model runs for other case studies, omission of the mutation step has led on occasions to greater deterioration in solution quality). It was decided to leave the option of implementing the mutation as a choice to be made by the user at run time. 6. Validity tests on the case study The previous section focussed on numerical behavior of the algorithm, and especially the effectiveness of the proximity measures in achieving spatial compactness. As a final step, we undertook a further set of practically oriented tests, using the case study of Section 3. In effect, we simulated how decision makers might use such a system, based on our experiences with the ‘touch table’ workshops described earlier, in order to evaluate whether performance of the algorithm was consistent with intuitive expectations in this real-life setting. We do emphasize that the results of this final section are hypothetical, and do not represent actual inputs from or a complete study with real decision makers. The simulations are, nevertheless, based on our experiences in working with real decision makers using the touch table. 6.1. Additive objectives As a start we focus on the additive objectives maximizing the value for agriculture, landscape and nature, in order to understand ranges of outcomes that are achievable when spatial objectives are ignored. These results are obtained using the simplified model based on the assignment algorithm. For purposes of comparison, four land use plans were generated, one in which the three additive objectives were assigned equal weights, and each of the others placing 98% of the weight on one of the three objectives (effectively maximizing one objective alone). The latter three alternatives

Table 2 Effects of genetic algorithm parameters on performance. Mutation proportion

No. in tournament

Prob. of reversal

Time (s)

Aggregate value

0 0.015 0.015 0.015 0.015 0.015

5 5 4 7 5 5

0.1 0.1 0.1 0.1 0.2 0.02

37.8 88.4 67.0 88.6 76.6 73.8

0.525 0.524 0.524 0.523 0.524 0.524

Spatial objective No. clusters

Rel. size

Shape

44.5 33.5 41.3 36.7 41.3 34.5

0.863 0.826 0.852 0.867 0.808 0.817

14.4 13.6 14.3 14.0 13.7 13.1

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx Table 3 Objective values of land use plans using different weights. Situation

Agriculture

Landscape

Nature

Current situation Equal weights Best for agriculture Best for landscape Best for nature

563 580 589 558 537

661 677 675 681 672

199 267 247 255 286

9

maximize additive values given by (8), subject to the constraints (1) and (2), and to the dim being binary. The same four sets of weights were used. The resulting ranges of values for the three objectives were extremely similar (547–581 for agriculture, 672– 681 for landscape, and 255–286 for nature). The similarity in ranges with and without the spatial objective foreshadows the results of the next paragraphs.

6.2. Inclusion of spatial objectives define the extremes of the possible solutions, which is useful as a starting set for subsequent negotiations. In optimizing the scalarizing function defined in (7), the reference levels g r may potentially have as much or more effect than the weights. Thus although a full investigation of interactive assessment of the reference points will form part of later implementation research, it was appropriate to allow some shift of the reference levels consistent with changed weights. For this purpose we defined the reference levels by f  ðv alue in the current solutionÞ þ ð1  f Þ  ð the estimated ideal v alueÞ. In the case of the equally-weighted scenario, we used f ¼ 0:7 (i.e. putting most emphasis on the current solution values) for each objective. For the other three cases, we set f ¼ 0:05 for the objective receiving maximum weight (i.e. the reference point close to ideal), and f ¼ 1:5 for the other objectives (i.e. shifting further away from ideal than the current solution). Recall that in the reference point approach (Miettinen et al., 2008, chapter 2), there is no requirement that the reference point be feasible. Table 3 shows the values for the three objectives for the four alternative land use plans, as well as for the existing (current) situation in the area. For purposes of interpretation, note that the subjectively assessed value scores were scaled so that values could in theory range from 0 (if all parcels have no value for this objective) to 1000 (if all parcels represent optimal conditions for this objective). Clearly the weights have the expected effects in the last three extreme alternatives, although the influences are small (the largest effect being a trade-off between the agriculture and nature objectives). The reason for the relatively moderate influences is probably that the areas where objectives linked to different land uses can be traded off is only part of the whole polder. Other parts of the polder are either exclusively suitable for agriculture or for nature. In order to put the range of variation across the four plans into perspective, we also used binary integer linear programming to

Having established the range of outcomes achievable for the additive objectives, we finally examine the influence of the spatial objectives on results for essentially the same problem. We have seen that the results are relatively insensitive to the weight selection for the additive objectives and thus the effect of introducing the spatial objective was undertaken for a moderate set of weights for the three additive objectives (in the ratios 1:2:1), representing the preferences observed in the earlier workshops. Fig. 4(a) shows the land use plan when only the additive objectives are maximized (using binary integer programming for the additive objectives as described above). The map shown in Fig. 4(b) results from inclusion of the spatial objective, using the proximity-based genetic algorithm, the same weights on the additive objectives, and with high priority on the proximity goals (a reference level of 0.825 for each of the contiguity measures P m ). The most obvious difference is the nearly total separation of the three land uses, and especially the large connected area of extensive agriculture in the center of the polder. The dramatic reduction in the numbers of small clusters of parcels of the same use is visually evident. In a practical implementation, the resulting maps would be used as a starting point for negotiations within the relevant groups using the map table. Users might modify the maps directly, or further optimizations could be run using different weights and/or reference points for the additive and spatial objectives. In follow-up research with real groups we plan to experiment with various forms of user interaction, but this would represent a substantial study in its own right. It is interesting to note that in this case study, the performance on the additive objectives are very similar when ignoring spatial goals, and when performing the full optimization including the spatial objective, a feature we had already observed above in

Fig. 4. Land use plan without (a) and with (b) spatial objectives.

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002

10

T.J. Stewart, R. Janssen / Computers, Environment and Urban Systems xxx (2014) xxx–xxx

Table 4 Objective values of land use plans with and without spatial goals. Objective (maximization of benefits to)

Agriculture

Landscape

Nature

Without spatial goals With spatial goals

581 555

678 672

281 266

comparing ranges with and without the spatial objective. For the three key additive objectives, the values (relative to a maximum of 1000) obtained for the two cases were as shown in Table 4. In this case study, therefore, it turns out that the spatial objectives can largely be realized without penalty to the additive objectives. This result is specific to the Bodegraven polder, and would not be expected to occur generally.

7. Conclusions The primary focus of this article is to describe a simplified definition of spatial objectives and an associated genetic algorithm that can be used to generate land use plans maximizing both additive and spatial objectives in a vector based GIS environment. The problem defined in this article is a non-linear combinatorial optimization problem. An important contribution to the non-linear characteristic is the set of indirect (spatial contiguity) objectives. Earlier grid based approaches showed that it was not possible to combine sufficient detail in the presentation of the problem with acceptable calculation times for use in an interactive setting. Using a vector representation it proved possible to represent the problem realistically with a response times of around one minute or less (if the mutation step is omitted), which from previous experiments had been shown to be an upper limit for effective use in an interactive setting. The GA approach is however supplemented by much faster algorithms which provide problem bounds although not satisfaction of spatial objective goals. Comparison of the land use maps generated with the simple and spatial algorithm demonstrated that it is possible to substantially improve the values of the spatial objectives for each land use without substantial loss in terms of the additive objectives. Most importantly the number of clusters and therefore fragmentation is dramatically reduced. In particular, the increased compactness for intensive agriculture is important for efficient management of the land. A useful extension to the algorithm would be to include the possibility to generate corridors. This would require specification of a entry and exit zones of the corridor and additional objectives such as length, width and quality. More importantly, perhaps, attention needs to be given to improving the integration of the algorithm into the GIS environment, so that changes made by users on the touch table can lead to immediate re-optimization. An important feature of such integration would the operationally meaningful guidance to users in adjusting importance weights and/or reference points. Such follow-up research has been commenced.

Acknowledgments The first author acknowledges with thanks the financial support from the South African National Research Foundation (NRF). This research was part of one of the projects of the consortium on the development of decision support tools funded by the ‘Knowledge for Climate’ research program of the Netherlands.

References Aerts, J. C. J. H., Eisinger, E., Heuvelink, G. B. M., & Stewart, T. J. (2003). Using linear integer programming for multi site land use allocation. Geographical Analysis, 35, 148–169.

Alexander, K. A., Janssen, R., Arciniegas, G. A., O’Higgins, T. G., Eikelboom, T., & Wilding, T. A. (2012). Interactive marine spatial planning: Siting tidal energy arrays around the Mull of Kintyre. PLoS ONE, 7, 1–9. Arciniegas, G. A., & Janssen, R. (2012). Spatial decision support for collaborative land use planning workshops. Landscape and Urban Planning, 107, 332–342. Arciniegas, G. A., Janssen, R., & Omtzigt, N. (2011). Map-based multicriteria analysis to support interactive land use allocation. International Journal of Geographical Information Science, 1–17. Arciniegas, G. A., Janssen, R., & Rietveld, P. (2013). Effectiveness of collaborative map-based decision support tools: Results of an experiment. Environmental Modelling & Software, 39, 159–175. Belton, V., & Stewart, T. J. (2002). Multiple criteria decision analysis: An integrated approach. Boston: Kluwer Academic Publishers. Cao, K., Huang, B., Wang, S., & Lin, H. (2011). Sustainable land use optimization using boundary-based fast genetic algorithm. Computers, Environment and Urban Systems Available on line. Carver, S. J. (1991). Integrating multi-criteria evaluation with geographical information systems. International Journal of Geographical Information Science, 5, 321–339. Chen, Y., Li, X., Liu, X., & Liu, Y. L. (2010). An agent-based model for optimal land allocation (AgentLA) with a contiguity constraint. International Journal of Geographical Information Science, 24, 1269–1288. Cornelisse, A. C., Janssen, R., Arciniegas, G. A., & Omtzigt, A. Q. A. (2007). Report on the framing workshop – Bodegraven stakeholders. Technical Report. Institute for Environmental Studies, VU University, Amsterdam. Provinciehuis Zuid Holland 26 oktober 2007 (In Dutch). Cova, T. J., & Church, R. L. (2000). Exploratory spatial optimization in site search: A neighborhood operator approach. Computers, Environment and Urban Systems, 24, 401–419. Deb, K. (2001). Multi-objective optimization using evolutionary algorithms. Chichester: John Wiley & Sons. Demetriou, D., See, L., & Stillwell, J. (2013). A spatial genetic algorithm for automating land partitioning. International Journal of Geographical Information Science, 27, 2391–2409. Geertman, S., & Stillwell, J. (Eds.). (2009). Planning support systems: Best practice and new methods. New York: Springer. Janssen, R., Arciniegas, G. A., & Verhoeven, J. T. A. (2013). Spatial evaluation of ecological qualities to support interactive design of land use plans. Environment and planning B: Planning and design, 40, 427–446. Janssen, R., van Herwijnen, M., Stewart, T. J., & Aerts, J. C. J. H. (2008). Multiobjective decision support for land use planning. Environment and planning B: Planning and design, 35, 740–756. Lau, K. H., & Kam, B. H. (2005). A cellular automata model for urban land-use simulation. Environment and planning B: Planning and design, 32, 247–263. Ligmann-Zielinska, A., Church, R. L., & Jankowski, P. (2008). Spatial optimization as a generative technique for sustainable multiobjective land-use allocation. International Journal of Geographical Information Science, 22, 601–622. Longley, P., Goodchild, M., Maguire, D., & Rhind, D. (2005). Geographic information systems and science. New York: Wiley. Malczewski, J. (2010). Multiple criteria decision analysis and geographic information systems. In M. Ehrgott, J. Figueira, & S. Greco (Eds.), Trends in multiple criteria decision analysis (pp. 369–395). New York: Springer. Matthews, K. B., Sibbald, A. R., & Craw, S. (1999). Implementation of a spatial decision support system for rural land use planning: Integrating geographic information system and environmental models with search and optimization algorithms. Computers and Electronics in Agriculture, 23, 9–26. Miettinen, K. (1999). Nonlinear multiobjective optimization. Norwell, Massachusetts: Kluwer Academic Publishers. Miettinen, K., Ruiz, F., & Wierzbicki, A. P. (2008). Introduction to multiobjective optimization: Interactive approaches. In J. Branke, K. Deb, K. Miettinen, & R. Słowin´ski (Eds.), Multiobjective optimization – Interactive and evolutionary approaches. Lecture notes in computer science (pp. 27–57). Berlin, Heidelberg: Springer-Verlag GmbH. Oxley, T., McIntosh, B., Mulligan, M., Winder, N., & Engelen, G. (2004). Integrated modelling and decision support tools: A mediterranean example. Environmental Modelling and Software, 19, 999–1010. Porta, J., Parapar, J., Doallo, R., Rivera, F. F., Santé & Crecente, R. (2013). High performance genetic algorithm for land use planning. Computers, Environment and Urban Systems, 37, 45–58. Santé-Riveira, I., Boullon-Magan, M., Crecente-Maseda, R., & Miranda-Barriós, D. (2008). Algorithm based on simulated annealing for land-use allocation. Computers and Geosciences, 34, 259–268. Stewart, T. J., Janssen, R., & van Herwijnen, M. (2004). A genetic algorithm approach to multiobjective land use planning. Computers and Operations Research, 32, 2293–2313. Strategiegroep Gouwe Wiericke (2007). Strategic Agenda Gouwe Wiericke. Technical Report. Royal Haskoning, Rotterdam. In Dutch. Strategiegroep Gouwe Wiericke (2010). Agreement on fen meadow Gouwe Wiericke fen meadow area (Veenweideconvenant), Provincie Zuid Holland, City Council of Bodegraven, Reeuwijk, Gouda, Vlist and Waddinxveen, Waterboard of Rijnland. Technical Report. Waterboard de Stichtse Rijnlanden. In Dutch. Vanegas, P., Cattrysse, D., & Van Orshoven, J. (2011). A multiple criteria heuristic solution method for locating near to optimal contiguous and compact sites in raster maps. In B. Murgante et al. (Eds.), Geocomputation, sustainability & environmental planning, SCI 348 (pp. 35–56). Berlin, Heidelberg: SpringerVerlag.

Please cite this article in press as: Stewart, T. J., & Janssen, R. A multiobjective GIS-based land use planning algorithm. Computers, Environment and Urban Systems (2014), http://dx.doi.org/10.1016/j.compenvurbsys.2014.04.002