Achieving a Sustainable Urban Form through Land Use Optimisation : Insights from Bekasi City ’ s Land-Use Plan ( 2010 – 2030 )

Cities worldwide have been trying to achieve a sustainable urban form to handle their rapid urban growth. Many sustainable urban forms have been studied and two of them, the compact city and the eco city, were chosen in this study as urban form foundations. Based on these forms, four sustainable city criteria (compactness, compatibility, dependency, and suitability) were considered as necessary functions for land use optimisation. This study presents a land use optimisation as a method for achieving a sustainable urban form. Three optimisation methods (particle swarm optimisation, genetic algorithms, and a local search method) were combined into a single hybrid optimisation method for land use in Bekasi city, Indonesia. It was also used for examining Bekasi city’s land-use-plan (2010–2030) after optimising current (2015) and future land use (2030). After current land use optimisation, the score of sustainable city criteria increased significantly. Three important centres of land use (commercial, industrial, and residential) were also created through clustering the results. These centres were slightly different from centres of the city plan zones. Additional land uses in 2030 were predicted using a nonlinear autoregressive neural network with external input. Three scenarios were used for allocating these additional land uses including sustainable development, government policy, and business-as-usual. Future land use allocation in 2030 found that the sustainable development scenario showed better performance compared to government policy and business-as-usual scenarios.


Introduction
Unsustainable city forms have been perceived as a source of environmental problems [1,2].They directly affects habitats, ecosystems, endangered species, water quality, habitat fragmentation, car use, open spaces, pollution, noise, global climate, replace natural cover with impervious surfaces, and other effects [2][3][4][5].Therefore, a city with sustainable development goals (SDGs) has to achieve a sustainable urban form.Such a city adheres to the sustainable development concept, i.e., development which caters both to present and future needs [6,7].
Many sustainable urban forms have been proposed.These include (i) the compact city; (ii) the eco city; (iii) neo-traditional development; and (iv) urban containment [2].The compact city and eco city were chosen in this study as a foundation to achieve a sustainable urban form within the study area, Bekasi city, Indonesia.The eco city was also included as a base model in this study since the sustainability of the compact city is still being debated [8].Based on the compact and eco city models, four criteria were chosen, compactness, compatibility, dependency and suitability, to determine the sustainability score for the city when optimising urban land use (LU).These criteria are similar to those of a previous study on current LU allocation [9].Although difficult to change an existing city form with newly optimised form, some activities dealing with urban structure and land use pattern can be implemented to achieve sustainable urban form as in the current study.Some studies implemented sustainable urban development and urban intensification concepts for those activities [3,[10][11][12][13].In the current study, the optimisation process accounted also for future LU.
Simple LU zoning has been implemented worldwide.These processes allocate each LU to appropriate locations based on their classes.In contrast, Alexander [14] criticised the LU zoning because of its tree formation instead of semi-lattice that makes the city unnatural and difficult to grow.Semi-lattice with inter-connections among LUs makes a city more natural, easy to grow and sustainable.Thus, Alexander suggested the planner calculating all LU relations instead of just proposing the LU zones based on the group of LUs.However, the semi-lattice concept of natural city is still under debate and difficult to manually calculate.The current study proposed a tool and method for calculating the sustainability score of all LU relationships while automatically generating a sustainable urban form using LU optimisation as an alternative to traditional methods.Additionally, the current study can be used for examining policy regarding sustainability since the study area has already implemented LU zoning in its city plan (2010-2030).A zoning system is affected by the street-block structure that sometimes creates problems.Since the city government has already established the street-block structure, the proposed model excluded it and only optimised other specific LUs.
The current study optimised both current LU and future LU through a transformation and an allocation model respectively.Whereas the transformation model created a suggested LU zone, the allocation model allocated newly predicted LUs in the future.In allocation model three exploratory scenarios were used: (i) a sustainable development scenario (SD); (ii) a government policy scenario (GP); and (iii) a business-as-usual scenario (BAU).
The main part of the current study is a LU optimisation model.Many techniques have been used and developed as candidates for LU optimisation methods such as fuzzy logic [15], genetic algorithms (GA) [15][16][17][18][19][20], particle swam optimization (PSO) [9,[21][22][23], and simulated annealing (SA) [24].These heuristic methods have advantages in handling discontinuities, non-linear or stochastic problems [25].Nowadays, there is a trend to use hybrid methods that combine several optimisation methods into a single and more powerful one [26][27][28].In the current study, the PSO, GA and local search optimisation methods were combined into a hybrid method for LU optimisation.
The specific objectives of the current study were: (1) to do an exploratory analysis of LU optimisation models and their use to achieve a sustainable urban form; (2) analysis of current LU optimisation with comparison to the Bekasi city plan (2010-2030); and (3) to develop future allocation of newly predicted LU through the devised modelling approach to examine different allocation scenarios in the study area.By achieving these objectives, the current study intended to fill the research gap about LU optimisation usage for achieving a sustainable urban form.

Study Area
The study area, Bekasi city, is located to the east of Jakarta province, with Bekasi regency on the east and north, and Bogor regency on the south (Figure 1).This city is part of the Jakarta Metropolitan Region (JMR) with an area of 210.49km 2 .More than 90 per cent of the land is residential and the remaining is for commercial, industrial, education, agricultural, and other uses [29].Agricultural production is limited due to the lack of cropland.Therefore, agricultural products must be obtained from other regions.Commercial and market places supply the basic needs of people [29].production is limited due to the lack of cropland.Therefore, agricultural products must be obtained from other regions.Commercial and market places supply the basic needs of people [29].Obtaining basic needs is a problem that is worsened by rapid population growth (Figure 1B).The population of more than 2.7 million in 2015 (population density 12,827/km 2 ), increases slightly every year due to the birth rate and migration [30].Serving people with proper housing, a liveable environment, and good roads are the main tasks of city government.Fortunately, Bekasi city is relatively flat (less than 2 per cent grade) and geologically suitable for housing, large buildings, and roads [31].Toll roads across this city ease access to other regions but create harmful air and sound Obtaining basic needs is a problem that is worsened by rapid population growth (Figure 1B).The population of more than 2.7 million in 2015 (population density 12,827/km 2 ), increases slightly every year due to the birth rate and migration [30].Serving people with proper housing, a liveable environment, and good roads are the main tasks of city government.Fortunately, Bekasi city is relatively flat (less than 2 per cent grade) and geologically suitable for housing, large buildings, and roads [31].Toll roads across this city ease access to other regions but create harmful air and sound pollution nearby.Other types of pollutions (water and soil pollution) also happen in Bantargebang sub-district as it is the central waste disposal of Bekasi and the capital city, Jakarta.
Since the area is relatively flat, low lying regions are prone to flooding in the rainy season from November to February [32].This problem is further compounded by its location next to Bogor regency which has a higher elevation and more rainfall [31].The local government has been struggling to minimize flooding by managing rivers, wetlands, and LUs.

Dataset
In addition to direct surveying, a variety of datasets were used, namely, Bing basemap on ArcGIS, Google Earth imagery (30 m spatial resolution), commercial sites, official government sites, thematic maps from the Geospatial Information Bureau of Indonesia (BIG), the yellow pages, and other spatial and non-spatial data sources.In this study, ten LU classes was chosen following the United States Geological Survey (USGS) standards [33] (10) high density residential areas.The residential area class was divided into high and low-density residential classes because of their different characteristics.As an example, regarding compatibility aspect between two LUs, a high-density residential LU usually conflict with another high-density residential LU but compatible with the low-density residential LU.
The land use and land cover (LULC) map from BIG is difficult to use for classifying built-up areas such as commercial, residential, and other LU classes.Hence, we classified each LU class manually using GIS tools based on available data sources.However, remote sensing data from USGS (August 2015) was used to validate current built-up LU data.For the residential LU class, an approximation model was used.Instead of one point for one building, we used one point for a 100 m radius with 150 to 200 buildings for high density residential areas and 10 to 150 buildings for low density residential areas.This is similar to lowering the resolution in raster dataset and such an approximation in modelling is common to enable more efficient computation [34].
Other thematic maps-roads, rivers, recreation areas, flood zones, open spaces, land price, and polluted areas (air, noise and soil) were used for suitability analysis.The locations of polluted areas were adapted from earlier works which showed regions near toll roads and central business as pollution centres [35].
For optimising future land use in 2030, historical imagery from Google Earth was used to see LU changes over periods of three consecutive years (2003, 2006, 2009, 2012, and 2015).This historical data was used for predicting future LU (2030).

Optimisation Module
This study intended to examine the feasibility of LU optimisation for achieving a sustainable urban form.Two robust evolutionary algorithms were used, particle swarm optimisation (PSO) and a genetic algorithm (GA).They were combined into a hybrid algorithm.Another method, a local search, was also employed to refine the final results.This hybrid algorithm was used to optimise the fitness function, which was created from four sustainability city criteria using an aggregating function method as hybrid multi-criteria evolutionary algorithm (HMCEA).This optimisation module adopted a constrained optimisation problem, a kind of optimisation with several restricted ranges or regions as constraints [36].
Figure 2A shows the HMCEA procedure for optimising both current (Figure 2B) and future LUs (Figure 2C).While the current LU optimisation was used for all LUs within the study area, future LU optimisation was only used for the predicted new LUs.As in other studies [9,23], the optimisation used the current LUs as the initial PSO particles.Future LU optimisation used random initial LUs before optimisation process.The first stage was PSO, which is a fast optimisation method that mimics a flock of birds or a school of fish (called particles) searching for food [37].Earlier studies showed that PSO is faster than GA but less accurate [9,21,23].A small initial weight might be used to increase accuracy, but it lowers the searching distance of particles and weakened their ability to jump over a wide constraint.One way to overcome this is excluding the particles with lower fitness scores [38].This cannot be implemented in the current study since every particle represented a particular land use.However, instead of deleting the particles with lower fitness scores, half of them were optimised in the next stage using GA [39].Without the moving characteristic in PSO, GA uses crossover and, especially, mutation, to jump over wide constraint regions.Additionally, computation is faster because it only optimises half of LUs from the PSO stage.At the end of this stage, the result is merged again with the particles with higher fitness scores from the PSO stage before continuing to the next stage.
The last stage was a local search method that refined the results by searching some locations around every optimised LU for the possibility of finding the more optimal result.Previous studies showed the benefit of integrating GA with a local search [40,41].In this study, a triangle pattern search method was used.It checked three locations around every result location (250 m coverage).
After the local search stage, the optimisation process continued for other LU classes before reaching a maximum number of iterations as a stop condition.Increasing the fitness score of every LU class also increased the total fitness score.PSO, GA, local search, and the GIS conversion module in Matlab software (Version 7.7 Release 2008b, Mathworks, Natick, MA, USA) were used and integrated.
HMCEA was used for optimising both current and future LUs. Figure 2B shows the current LU optimisation for creating cluster-based centres of LU classes that can also be used to propose LU zones.This optimisation used empty and allowable LUs within the study area as constraints.Constraints for future LU optimisation were based on a scenario that is discussed below.

Criteria Functions
LU optimisation is a kind of LU transformation.It needs some conversions between LU classes [42].Other studies included resistance-to-change or conversion cost factors in their criteria to manage the conversion problem [17,23].Relocating a building is difficult not only in Bekasi city but also in most of the regions in Indonesia.The current study did not use conversion cost factor, so every LU freely searched for the optimum location.LU composition suggested from LU optimisation was The first stage was PSO, which is a fast optimisation method that mimics a flock of birds or a school of fish (called particles) searching for food [37].Earlier studies showed that PSO is faster than GA but less accurate [9,21,23].A small initial weight might be used to increase accuracy, but it lowers the searching distance of particles and weakened their ability to jump over a wide constraint.One way to overcome this is excluding the particles with lower fitness scores [38].This cannot be implemented in the current study since every particle represented a particular land use.However, instead of deleting the particles with lower fitness scores, half of them were optimised in the next stage using GA [39].Without the moving characteristic in PSO, GA uses crossover and, especially, mutation, to jump over wide constraint regions.Additionally, computation is faster because it only optimises half of LUs from the PSO stage.At the end of this stage, the result is merged again with the particles with higher fitness scores from the PSO stage before continuing to the next stage.
The last stage was a local search method that refined the results by searching some locations around every optimised LU for the possibility of finding the more optimal result.Previous studies showed the benefit of integrating GA with a local search [40,41].In this study, a triangle pattern search method was used.It checked three locations around every result location (250 m coverage).
After the local search stage, the optimisation process continued for other LU classes before reaching a maximum number of iterations as a stop condition.Increasing the fitness score of every LU class also increased the total fitness score.PSO, GA, local search, and the GIS conversion module in Matlab software (Version 7.7 Release 2008b, Mathworks, Natick, MA, USA) were used and integrated.
HMCEA was used for optimising both current and future LUs. Figure 2B shows the current LU optimisation for creating cluster-based centres of LU classes that can also be used to propose LU zones.This optimisation used empty and allowable LUs within the study area as constraints.Constraints for future LU optimisation were based on a scenario that is discussed below.

Criteria Functions
LU optimisation is a kind of LU transformation.It needs some conversions between LU classes [42].Other studies included resistance-to-change or conversion cost factors in their criteria to manage the conversion problem [17,23].Relocating a building is difficult not only in Bekasi city but also in most of the regions in Indonesia.The current study did not use conversion cost factor, so every LU freely searched for the optimum location.LU composition suggested from LU optimisation was difficult for implementation.However, some important centres created through clustering the result can be used as LU zoning suggestion and an instrument for analysing the sustainability of the city LU plan (2010-2030).
In the current study, four criteria functions necessary for sustainable urban form were considered when modelling the problem.The first criterion is compactness based on the compact city form.Many studies have shown that this criterion will minimise travel distances and fuel consumption [4,[43][44][45].However, other studies suggested addition of mixed use areas to avoid social inequality, environment problems, and low quality of life [8].The second criterion was compatibility, a criterion for avoiding conflicting between two LU classes and minimising the negative effects.In the current study, compatibility would work together with the third criterion, dependency.Whereas compatibility described the conflicting of two land uses, dependency described the attracting of two land uses.They were chosen as additional criteria to ensure mixed use areas (will be discussed in detail in the next section).The last criterion was suitability.This was based on the eco city form to increase the quality of life by implementing suitability analysis for every land use class [2,7,46].HMCEA was used for maximizing compactness, compatibility, dependency, and suitability of four criteria functions (F 1 , F 2 , F 3 , and F 4 ) as defined in Equations ( 1)-( 4).
where n is the number of LUs in the study area, i and j are the current LU and its neighbour respectively.Compactness, Comp ij , Dep ij , and Suitability Score are criteria values discussed in the following sections.

Compactness
There are three aspects of compact city: density, mix of uses, and intensification [8].In the current study, compactness intended to gain the density and to avoid urban sprawl.Other criteria (compatibility and dependency) intended to enhance another aspect, mix of uses.
The scattered nature of urban sprawl can be both costly and unsustainable [47].Therefore, compactness criterion was used to reduce the scattered of sprawl.It based on the fact that some LU classes such as residential, industrial, and commercial are preferred to be in the vicinity of similar LUs.The optimisation module used compactness to attract other similar LUs, especially in scattered locations.Equation ( 5) was used to calculate compactness [9,47].
where I equals one if a land use is the same as a neighbour; n is the number of LUs in class i; h is the number of all LU inside the effectual region as shown in Figure 3.If all LUs inside region R 1 have the same class, the compactness score is maximum (one), but the other compact city criteria (compatibility and dependency) force to enhance the mix of uses (in optimisation process) through mixing with other appropriated LU classes in this region.This situation will decrease the compactness score but increase overall optimum score after aggregation with other criteria.mixing with other appropriated LU classes in this region.This situation will decrease the compactness score but increase overall optimum score after aggregation with other criteria.The choice of R 1 and R 2 were based on previous studies [45] in which it was found that at a distance of 5 to 15 km, the influence (e.g., economic, health, socio-cultural) of its neighbour was not significant.In this study, 5 km distance was chosen for R 1 and 6 km distance was chosen for R 2 .If a neighbour was located between 5 km and 6 km, a fuzzy set factor was used (Equation ( 6)) and assumed β = 1 (linear) for more efficient calculation [9,48].
where  is a membership function value (from 0 to 1) scored from the distance to its neighbour, ; i and j are the LU under consideration and its neighbour respectively; is minimum distance of the effectual region; is maximum distance of the effectual region.The effectual region was used not only for determining I variable but also the h variable in Equation ( 5).Compatibility and dependency criteria (which are highly correlated to compactness criterion in supporting the mixed use of compact city form) also used Figure 3 and Equation (6) for scoring.

Compatibility and Dependency
A previous study [8] showed that although compact city has been accepted as a sustainable urban form, without mixed use areas, the effect of minimal travel distance will not be significant.People who live in a city interact with more than one LU class to work, learn, and for other social activities.Therefore, to achieve a sustainable urban form, one LU class needs other classes.However, the influence of one class upon another must be considered since it is possible to create negative effects.In the current study, compatibility, i.e., the degree to which there is a coexistence among LU classes without a significant negative impacts, was required as criteria functions [9,49].Dependency, which is the amount of satisfaction with positive effects among land use classes, was also used [9].
To describe the compatibility and dependency between LUs, the Delphi method was used as a framework for creating compatibility and dependency matrices.This method is an iterative process to achieve consensus among a group of experts on a particular topic and is especially useful in situations where no standard criteria exist for evaluation [50].Two-round Delphi survey was chosen and the questionnaire was sent to 20 research participants (potential users, local government staff, lecturers, and experts within the study area).Before scoring compatibility and dependency between The choice of R 1 and R 2 were based on previous studies [45] in which it was found that at a distance of 5 to 15 km, the influence (e.g., economic, health, socio-cultural) of its neighbour was not significant.In this study, 5 km distance was chosen for R 1 and 6 km distance was chosen for R 2 .If a neighbour was located between 5 km and 6 km, a fuzzy set factor was used (Equation ( 6)) and assumed β = 1 (linear) for more efficient calculation [9,48].
where u A is a membership function value (from 0 to 1) scored from the distance to its neighbour, d ij ; i and j are the LU under consideration and its neighbour respectively; d min ij is minimum distance of the effectual region; d max ij is maximum distance of the effectual region.The effectual region was used not only for determining I variable but also the h variable in Equation (5).Compatibility and dependency criteria (which are highly correlated to compactness criterion in supporting the mixed use of compact city form) also used Figure 3 and Equation (6) for scoring.

Compatibility and Dependency
A previous study [8] showed that although compact city has been accepted as a sustainable urban form, without mixed use areas, the effect of minimal travel distance will not be significant.People who live in a city interact with more than one LU class to work, learn, and for other social activities.Therefore, to achieve a sustainable urban form, one LU class needs other classes.However, the influence of one class upon another must be considered since it is possible to create negative effects.In the current study, compatibility, i.e., the degree to which there is a coexistence among LU classes without a significant negative impacts, was required as criteria functions [9,49].Dependency, which is the amount of satisfaction with positive effects among land use classes, was also used [9].
To describe the compatibility and dependency between LUs, the Delphi method was used as a framework for creating compatibility and dependency matrices.This method is an iterative process to achieve consensus among a group of experts on a particular topic and is especially useful in situations where no standard criteria exist for evaluation [50].Two-round Delphi survey was chosen and the questionnaire was sent to 20 research participants (potential users, local government staff, lecturers, and experts within the study area).Before scoring compatibility and dependency between two LU classes, the questionnaire first explained to the research participants that two attracting LUs should have high dependency score and vice versa.However, between two conflicting LUs, they should have the low compatibility score and vice versa.The first-round responses were discussed, shared, and used as a basis for second-round.The second-round questionnaire was released to the research participants and when completed, returned for analysis.The participants were given the opportunity to change or expand their round-one response based on previous discussion and other participants' answer.Various participants with their various scores might cause there were exist some LU classes having the same scores both for compatibility and dependency after averaging.However, the results reflected the characteristic of Bekasi city as the study area (Table 1).The levels considered for compatibility and dependency of LUs are vhgh, high, med, low, and vlow which represented values of very high, high, medium, low, and very low, respectively.In the implementation, these variables were converted to numerical values of 0.428, 0.275, 0.176, 0.081, and 0.041, respectively, by following the Analytic Hierarchy Process (AHP) conversion standard before putting them in matrix form for calculating fitness scores [51].Equations ( 7) and ( 8) were used for calculating the compatibility (Comp ij ) and dependency (Dep ij ) of a LU.
where µ A d ij is membership function based on distance of a LU to its neighbor using a fuzzy set (Equation ( 6) and Figure 3).

Suitability
A sustainable urban form is complex and considers travel distance, fuel consumption, and other compact city characteristics [3].This type of form should be managed to avoid environmental and ecological degradation [2].Hence, suitability analysis that considers other important aspects (environment, access, and disaster risk, among others) must be used in LU optimisation to achieve another kind of sustainable urban form, the eco-city.Figure 4 shows the process to create a suitable region.another kind of sustainable urban form, the eco-city.Figure 4 shows the process to create a suitable region.Six thematic maps were used: distance from a pollution source, distance to a road, distance to floodplain, distance to an open space, distance to a lake/recreational area, and land price.Weighted ratios from each thematic map were gathered from a survey and used in an overlay process.The pairwise comparison with AHP method was used in the survey to generate weighted overlay ratio.Nine experts (planners, lecturers, and local government staff) were chosen as research participants to fill the questionnaire.The expert choice application was used to find the weighted overlay ratios of Figure 4A,B.Both inconsistency scores were below 0.1 that proved the weighted overlay ratios validity for implementation (based on AHP references) [51,52].The weighted overlay ratio of Figure 4A for distance to pollution, flood zone, recreation, open space, main road, and land price were 0.37, 0.275, 0.112, 0.107, 0.073, and 0.064, respectively.The weighted overlay ratio of Figure 4B for distance to roads and flood risk were 0.77 and 0.23, respectively.Three regions were created (low suitability, medium suitability, and high suitability) for residential areas (Figure 4A) and for commercial, industrial, health, and educational areas (Figure 4B).The suitability analysis in Figure 4 only considered the static factors and used other criteria, i.e., compatibility and dependency, to handle the dynamic factors (the location of residential, commercial, industrial, medical, etc.) that would probably move when the optimisation process run.In the current study, the scores for high, medium, and low for suitability regions were 0.428, 0.176, and 0.041, respectively.
To calculate a suitability score, every LU was checked to determine whether it was inside a low, medium, or high suitability region.An in-polygon algorithm in Matlab software was chosen to do these calculations [53,54].For example, if a residential LU is inside a medium region after an in-polygon algorithm process, the suitability score is 0.176.

Multi-Criteria and Constraint Handling
For handling multiple criteria, the current study used the aggregating function method [28].This method is appropriate when used for comparing many scenarios.An aggregating function, F, used a weight ratio for calculating total fitness in optimisation problem (Equations ( 9) and ( 10)),

= max
= Inside Allowable Location (10) Six thematic maps were used: distance from a pollution source, distance to a road, distance to floodplain, distance to an open space, distance to a lake/recreational area, and land price.Weighted ratios from each thematic map were gathered from a survey and used in an overlay process.The pairwise comparison with AHP method was used in the survey to generate weighted overlay ratio.Nine experts (planners, lecturers, and local government staff) were chosen as research participants to fill the questionnaire.The expert choice application was used to find the weighted overlay ratios of Figure 4A,B.Both inconsistency scores were below 0.1 that proved the weighted overlay ratios validity for implementation (based on AHP references) [51,52].The weighted overlay ratio of Figure 4A for distance to pollution, flood zone, recreation, open space, main road, and land price were 0.37, 0.275, 0.112, 0.107, 0.073, and 0.064, respectively.The weighted overlay ratio of Figure 4B for distance to roads and flood risk were 0.77 and 0.23, respectively.Three regions were created (low suitability, medium suitability, and high suitability) for residential areas (Figure 4A) and for commercial, industrial, health, and educational areas (Figure 4B).The suitability analysis in Figure 4 only considered the static factors and used other criteria, i.e., compatibility and dependency, to handle the dynamic factors (the location of residential, commercial, industrial, medical, etc.) that would probably move when the optimisation process run.In the current study, the scores for high, medium, and low for suitability regions were 0.428, 0.176, and 0.041, respectively.
To calculate a suitability score, every LU was checked to determine whether it was inside a low, medium, or high suitability region.An in-polygon algorithm in Matlab software was chosen to do these calculations [53,54].For example, if a residential LU is inside a medium region after an in-polygon algorithm process, the suitability score is 0.176.

Multi-Criteria and Constraint Handling
For handling multiple criteria, the current study used the aggregating function method [28].This method is appropriate when used for comparing many scenarios.An aggregating function, F, used a weight ratio for calculating total fitness in optimisation problem (Equations ( 9) and ( 10)), H = Inside Allowable Location (10) where w i ≥ 0 are the weighting coefficients representing the relative importance of the k criteria functions; F i is the criteria function of criterion i from Equations ( 1) to ( 4) and H is constraint handling based on LU class and scenario.The weighted ratio was normalised using Equation (11).
As in previous research, a survey of experts was conducted to determine the appropriate weighted ratios of the criteria function [17,23].In the current study, the weighted ratio from the survey was checked to ensure that no criterion dominated another.If there is no dominant criterion in a weighted ratio, it is feasible for use.In the current study, a Pareto-front optimal set was used to test each weighted ratio.
To create a Pareto-front optimal set, HMCEA runs several times with different weighted ratios.Every weighted ratio had one total score and four criteria scores.After normalising criteria scores, the Pareto-front function in Matlab software was used for sorting the non-dominant weighted ratios from the dominant ones.The resulting set became a Pareto-front optimal set as a base for checking each new weighted ratio to determine if it was dominant or non-dominant.
Constraint handling is also important since this study adopted constraint-based optimisation.Previous studies [42] faced difficulties handling multiple constraints in GA.Research by [9] used mapping of continuous particles to discrete forms to avoid some constraints in PSO.In the current study, a death penalty was used that forced every violated candidate to its previous location.It was chosen since it is adequate, fast and simple [36].The penalty function was adopted an in-polygon function in Matlab.This is an algorithm to check whether a tested point location is inside or outside a polygon [53,54].Constraints used in this study were based on LU class and scenario.

Scenario
Three scenarios were used for allocating predicted new LUs in 2030.These were government policy (GP), sustainable development (SD) and business-as-usual (BAU).The scenarios were different in regard to criteria function, constraints and weighted ratios (Table 2).First, the GP scenario was based on the Bekasi city LU plan (2030) which adopted a LU zoning system [32].This LU zoning describes the control of specific LU classes including commercial, industrial, and residential areas [32,55].From these zones, constraints were created for commercial, industrial and residential classes (Figure 5A).This scenario used criteria functions similar to the SD scenario, but only locating LUs according to appropriate zones.Secondly, the SD scenario used all criteria functions that were sustainable development factors with two regions kept from new LU as constraints: the open space areas and industrial zones as constraints.In this study, open space (20 per cent of total city area) was chosen from low suitability regions.For health issue considerations, the SD scenario also adopted the industrial zones of GP constraint, which was only used for industrial classes (Figure 5B).Finally, BAU scenario used empty LUs as constraints.This scenario purely used the compact city criteria (compactness, compatibility, and dependency) without suitability criterion (eco city criteria).Without a suitability criterion, this scenario depicts real conditions in this city where environmental factors are usually ignored (pollution, floods, etc.) when choosing LU locations (Figure 5C).
Using parameters in Table 2, HMCEA was run for allocating new LUs in 2030.These three scenarios were used for analysing the effect of sustainable city criteria in optimisation and ensuring that optimisation would occur in the future.environmental factors are usually ignored (pollution, floods, etc.) when choosing LU locations (Figure 5C).Using parameters in Table 2, HMCEA was run for allocating new LUs in 2030.These three scenarios were used for analysing the effect of sustainable city criteria in optimisation and ensuring that optimisation would occur in the future.

Prediction of New LU (2030)
This study not only optimised the current LU but also optimised new LUs in the future (2030) using an allocation model.This is a model that allocates every new LU to an empty and allowable location [42].Additional LUs in 2030 were predicted based on historical data using a nonlinear autoregressive neural network with external input (NARXNET) module in Matlab software.This module was based on Multilayer Perceptron Neural Network (MLPNN) that had better performance in predicting non-linear data than other conventional methods [5].In the current study, the NARXNET module used one hidden layer with 10 neurons and 2 delays.Delay in NARXNET is a parameter that uses prior trend to be trained for prediction.Previous LU data in time-series format and population projection as external input were trained and validated in Matlab before used to predict additional LUs in 2030.The selection of prediction year (2030) was based on the fact that the local development authorities have prepared a land-use plan (2010-2030).

Pareto-Front Optimal Set
The Pareto-front optimal set was used in the current study to ensure that all criteria participated in LU optimisation process.Figure 6 shows a Pareto-front optimal set (contained more than sixty weight ratios) that also contains weighted ratio from the survey.Since the weighted ratio from survey was the member of Pareto optimal set, there was no dominant criterion in this ratio and sufficient for simulation.Figure 6 also shows that there were conflicts between criteria.In these cases, one criterion score was increasing while the others were decreasing.Compatibility and dependency criteria were not conflict with each other where both criteria scores were either increasing or decreasing.Since there is no conflict between compatibility and dependency criteria, one of these criteria is allowable to be excluded without affecting the result and create a three-dimensional chart as shown in Figure 6

Prediction of New LU (2030)
This study not only optimised the current LU but also optimised new LUs in the future (2030) using an allocation model.This is a model that allocates every new LU to an empty and allowable location [42].Additional LUs in 2030 were predicted based on historical data using a nonlinear autoregressive neural network with external input (NARXNET) module in Matlab software.This module was based on Multilayer Perceptron Neural Network (MLPNN) that had better performance in predicting non-linear data than other conventional methods [5].In the current study, the NARXNET module used one hidden layer with 10 neurons and 2 delays.Delay in NARXNET is a parameter that uses prior trend to be trained for prediction.Previous LU data in time-series format and population projection as external input were trained and validated in Matlab before used to predict additional LUs in 2030.The selection of prediction year (2030) was based on the fact that the local development authorities have prepared a land-use plan (2010-2030).

Pareto-Front Optimal Set
The Pareto-front optimal set was used in the current study to ensure that all criteria participated in LU optimisation process.Figure 6 shows a Pareto-front optimal set (contained more than sixty weight ratios) that also contains weighted ratio from the survey.Since the weighted ratio from survey was the member of Pareto optimal set, there was no dominant criterion in this ratio and sufficient for simulation.Figure 6 also shows that there were conflicts between criteria.In these cases, one criterion score was increasing while the others were decreasing.Compatibility and dependency criteria were not conflict with each other where both criteria scores were either increasing or decreasing.Since there is no conflict between compatibility and dependency criteria, one of these criteria is allowable to be excluded without affecting the result and create a three-dimensional chart as shown in Figure 6 (3D plot of three criteria).Excluding dependency is recommended since its weight is smaller than compatibility.(3D plot of three criteria).Excluding dependency is recommended since its weight is smaller than compatibility.
Figure 6.Pareto-front optimal set showing weighted ratio from survey data with ratios for compatibility, dependency, compactness, and suitability of 0.402, 0.237, 0.115, and 0.247 respectively.

Current LU Optimisation
After checking its Pareto-optimality, a weighted ratio from survey data was used for optimising the current LU.The constraint of current LU optimisation was every empty and allowable location within Bekasi city.It restricted HMCEA to search optimum location outside Bekasi city and protected areas for new LUs (roads, lakes, rivers, etc.).To run the HMCEA, the existing LUs (Figure 1C) and suitability regions (Figure 6) were imported from ArcGIS.For ecological conservation reasons, the park class was unchanged but was still used in scoring.More than six thousand LU points were optimised until the total fitness score reached its saturation.
Figure 7 shows an increase in the total fitness score from 2.68 to 2.785 (3.92 per cent).This current LU optimisation was ideal since it ignored conversion cost or resistance to change factors.Despite this, it might be difficult or impossible to follow.However, the optimised LUs can be used to depict the spread of each of the LU classes.It also can be used to show the LU class centres for examining the Bekasi city LU plan (2010-2030).
Figure 8A shows the optimised LUs and that they are more compact than they initially were (Figure 1C).The simple fuzzy c-means clustering (FCM) method in Matlab software was used for clustering three important LU classes in the city plan, i.e., the residential, commercial, and industrial classes.Figure 8B shows three high and low density centres of the residential class.This residential class became denser and was located inside the high density residential zone called for by government policy.The residential centres might be appropriate for vertical buildings, which has been stated in the city plan (vertical building should be done in the high density region) [32].Figure 8C shows six commercial centres with one centre that is a more densely commercial LU.The denser LUs can be used as a central business district and the others as sub-central business districts.Though

Current LU Optimisation
After checking its Pareto-optimality, a weighted ratio from survey data was used for optimising the current LU.The constraint of current LU optimisation was every empty and allowable location within Bekasi city.It restricted HMCEA to search optimum location outside Bekasi city and protected areas for new LUs (roads, lakes, rivers, etc.).To run the HMCEA, the existing LUs (Figure 1C) and suitability regions (Figure 6) were imported from ArcGIS.For ecological conservation reasons, the park class was unchanged but was still used in scoring.More than six thousand LU points were optimised until the total fitness score reached its saturation.
Figure 7 shows an increase in the total fitness score from 2.68 to 2.785 (3.92 per cent).This current LU optimisation was ideal since it ignored conversion cost or resistance to change factors.Despite this, it might be difficult or impossible to follow.However, the optimised LUs can be used to depict the spread of each of the LU classes.It also can be used to show the LU class centres for examining the Bekasi city LU plan (2010-2030).
Figure 8A shows the optimised LUs and that they are more compact than they initially were (Figure 1C).The simple fuzzy c-means clustering (FCM) method in Matlab software was used for clustering three important LU classes in the city plan, i.e., the residential, commercial, and industrial classes.Figure 8B shows three high and low density centres of the residential class.This residential class became denser and was located inside the high density residential zone called for by government policy.The residential centres might be appropriate for vertical buildings, which has been stated in the city plan (vertical building should be done in the high density region) [32].Figure 8C shows six commercial centres with one centre that is a more densely commercial LU.The denser LUs can be used as a central business district and the others as sub-central business districts.Though these plans might be difficult or impossible to implement, the centres can be used for guiding the LU zones.these plans might be difficult or impossible to implement, the centres can be used for guiding the LU zones.these plans might be difficult or impossible to implement, the centres can be used for guiding the LU zones.

LU Allocation in 2030
Table 3 shows the prediction results using the NARXNET module for allocation in 2030.The number of LUs tended to gradually increase except for parks and the low density residential class.Additional LU predictions in 2030 (as shown in Table 3 in italic) will be optimised based on three scenarios using HMCEA.Figure 9 shows fitness scores of the three scenarios after 38 optimisation runs.In 2015 all fitness scores were similar (2.68) and decreased when additional LUs were randomly added.HMCEA optimised LUs until it reached the saturation value with the fitness scores for GP, SD, and BAU scenarios of 2.676, 2.683, and 2.643, respectively.Figure 9 also shows that additional LUs tend to lower the fitness scores, although LU optimisation has been implemented as a signal for city government to control LU growth, except the SD scenario (0.11 per cent above the initial fitness score in 2015).

LU Allocation in 2030
Table 3 shows the prediction results using the NARXNET module for allocation in 2030.The number of LUs tended to gradually increase except for parks and the low density residential class.Additional LU predictions in 2030 (as shown in Table 3 in italic) will be optimised based on three scenarios using HMCEA.Figure 9 shows fitness scores of the three scenarios after 38 optimisation runs.In 2015 all fitness scores were similar (2.68) and decreased when additional LUs were randomly added.HMCEA optimised LUs until it reached the saturation value with the fitness scores for GP, SD, and BAU scenarios of 2.676, 2.683, and 2.643, respectively.Figure 9 also shows that additional LUs tend to lower the fitness scores, although LU optimisation has been implemented as a signal for city government to control LU growth, except the SD scenario (0.11 per cent above the initial fitness score in 2015).The BAU scenario only used purely compact city criteria (excluding suitability criterion) in LU optimisation and had smallest fitness score since the final process suitability was also used in scoring.GP and SD scenarios had the same use of criteria functions, but the GP scenario had limited locations for LU candidates because it must follow LU zoning.This might be the source of its lower fitness The BAU scenario only used purely compact city criteria (excluding suitability criterion) in LU optimisation and had smallest fitness score since the final process suitability was also used in scoring.GP and SD scenarios had the same use of criteria functions, but the GP scenario had limited locations for LU candidates because it must follow LU zoning.This might be the source of its lower fitness score.After LU optimisation, the result was exported from Matlab into a shapefile for analysis of the optimisation results using a GIS tool.
As shown in Figure 10, the SD scenario had not only a higher fitness score but also wider open spaces.New residential LUs spread in south-west region can be considered as a high density residential zone in next city plan (2030-2050).
score.After LU optimisation, the result was exported from Matlab into a shapefile for analysis of the optimisation results using a GIS tool.
As shown in Figure 10, the SD scenario had not only a higher fitness score but also wider open spaces.New residential LUs spread in south-west region can be considered as a high density residential zone in next city plan (2030-2050).

The Influence of Compact City, Eco City and LU Zoning in Optimisation
The current study only used two kinds of sustainable urban form as a foundation: compact city and eco city.However, these forms were successfully used for LU optimisation with their four indicators as criteria function.From the compact city form, compactness, compatibility, and dependency criteria were used not only for increasing the density of the LUs but also ensuring the diversity of LUs in a region.This situation would reduce the fuel consumption and make people walk or use other green transportations.
Eco-city form needs a suitability analysis to check the suitability of an LU regarding the environment factors such as pollution, water, vegetation, disaster risk, etc.With the ratio of 0.247, the suitability criterion has a significant impact on the sustainability score in the study area.BAU scenario had lower sustainable score than the other scenarios because it did not use suitability criterion in its criteria function (Figure 9).SD scenario only had slight advantage compared to GP scenario which adopted LU zoning.Zones acted as constraints that limit the searching area in optimisation.It caused the lower sustainable score compared to SD scenario.However, since the LU zones within the study area was created after suitability analysis by city planners [32], the effect of limited searching area was not too significant.

Limitations
The current study has some limitations.Firstly, despite this research only used four criteria for LU optimisation, they have represented the compact city and eco city indicators.However, additional criteria might be useful to add.Secondly, the HMCEA has successfully handled ten LU classes.Some more specific classes might be involved in optimisation for illustrating the optimised LUs in detail.For example, commercial areas might be divided into modern market, traditional market, services,

The Influence of Compact City, Eco City and LU Zoning in Optimisation
The current study only used two kinds of sustainable urban form as a foundation: compact city and eco city.However, these forms were successfully used for LU optimisation with their four indicators as criteria function.From the compact city form, compactness, compatibility, and dependency criteria were used not only for increasing the density of the LUs but also ensuring the diversity of LUs in a region.This situation would reduce the fuel consumption and make people walk or use other green transportations.
Eco-city form needs a suitability analysis to check the suitability of an LU regarding the environment factors such as pollution, water, vegetation, disaster risk, etc.With the ratio of 0.247, the suitability criterion has a significant impact on the sustainability score in the study area.BAU scenario had lower sustainable score than the other scenarios because it did not use suitability criterion in its criteria function (Figure 9).SD scenario only had slight advantage compared to GP scenario which adopted LU zoning.Zones acted as constraints that limit the searching area in optimisation.It caused the lower sustainable score compared to SD scenario.However, since the LU zones within the study area was created after suitability analysis by city planners [32], the effect of limited searching area was not too significant.

Limitations
The current study has some limitations.Firstly, despite this research only used four criteria for LU optimisation, they have represented the compact city and eco city indicators.However, additional criteria might be useful to add.Secondly, the HMCEA has successfully handled ten LU classes.Some more specific classes might be involved in optimisation for illustrating the optimised LUs in detail.For example, commercial areas might be divided into modern market, traditional market, services, etc. Thirdly, the current research did not consider the street-block structure as LU class since they have already been established around Bekasi city.Fourthly, residential LU approximation in this study might reduce the accuracy.Therefore, additional LU optimisation for higher resolution scale (below 1:20,000 scales) or sub-district region might be useful.Finally, optimising and analysing the other regions around Bekasi city (Bogor and Bekasi regency) should be conducted since they might affect Bekasi city, especially for residential LU class.

Conclusions
This study has shown that LU optimisation can be used to achieve a sustainable urban structure as a path to a sustainable form.Despite difficult to follow the optimisation result, current LU optimisation found centres of important LU classes within the study area that can be used for analysing the current city plan (2010-2030).Four sustainable city criteria were successfully optimised using HMCEA for the present (2015) and future LU (2030).Future LU optimisation allocated new LUs in 2030 and can be used as a tool for revising and preparing the next city plan (2030-2050).The SD scenario, which only used sustainable city criteria, showed a better result compared to GP (combined with LU zoning) and BAU scenarios (without suitability criterion).Although the sustainability of the compact city is still being debated, its combination with the other sustainable urban forms might be more reliable.The results of proposed LU optimisation could help planners in considering the relations among LUs since they are difficult to calculate manually.Additionally, the results can be used to quantitatively analyse other urban forms (centrism or decentrism) or other kinds of regions such as peri-urban, rural, desert, and forests.

Figure 1 .
Figure 1.Study area: (A) map of Bekasi city; (B) population growth in the study area from 2003 to 2015 (Bureau of Statistics 2015); and (C) existing LU.

Figure 1 .
Figure 1.Study area: (A) map of Bekasi city; (B) population growth in the study area from 2003 to 2015 (Bureau of Statistics 2015); and (C) existing LU.

Figure 2 .
Figure 2. Diagrams showing (A) HMCEA as optimisation module; (B) optimising current land use and clustering; (C) optimising new land uses.

Figure 2 .
Figure 2. Diagrams showing (A) HMCEA as optimisation module; (B) optimising current land use and clustering; (C) optimising new land uses.

Figure 3 .
Figure 3.The influence region of a LU to the others.(A) Effectual Region; and (B) Fuzzy area between R 1 and R 2 .

Figure 3 .
Figure 3.The influence region of a LU to the others.(A) Effectual Region; and (B) Fuzzy area between R 1 and R 2 .

Figure 4 .
Figure 4. Spatial analysis using a weighted overlay technique to create suitable regions for scoring (A) residential areas; and (B) commercial, industrial, health, and educational areas.

Figure 4 .
Figure 4. Spatial analysis using a weighted overlay technique to create suitable regions for scoring (A) residential areas; and (B) commercial, industrial, health, and educational areas.

Figure 5 .
Figure 5. HMCEA constraint for three scenarios, (A) Bekasi city LU plan (2010-2030) as GP constraint; (B) open space and industrial zone as SD constraint; and (C), empty LUs for BAU constraint.

Figure 5 .
Figure 5. HMCEA constraint for three scenarios, (A) Bekasi city LU plan (2010-2030) as GP constraint; (B) open space and industrial zone as SD constraint; and (C), empty LUs for BAU constraint.

Figure 6 .
Figure 6.Pareto-front optimal set showing weighted ratio from survey data with ratios for compatibility, dependency, compactness, and suitability of 0.402, 0.237, 0.115, and 0.247 respectively.

Figure 7 .
Figure 7. Result graph showing current LU optimisation performance.

Figure 8 .
Figure 8. (A) Current optimisation results; (B) cluster of residential areas; (C) clusters of commercial areas; and (D) clusters industrial areas.

Figure 7 .
Figure 7. Result graph showing current LU optimisation performance.

Figure 7 .
Figure 7. Result graph showing current LU optimisation performance.

Figure 8 .
Figure 8. (A) Current optimisation results; (B) cluster of residential areas; (C) clusters of commercial areas; and (D) clusters industrial areas.

Figure 8 .
Figure 8. (A) Current optimisation results; (B) cluster of residential areas; (C) clusters of commercial areas; and (D) clusters industrial areas.

Figure 9 .
Figure 9. New LU allocation results that show fitness scores of each scenario.

Figure 9 .
Figure 9. New LU allocation results that show fitness scores of each scenario.
Figure 10 suggested the allocation of future new LUs because of LU growth based on three scenarios that useful for planning.Maintaining the open space and park, locating the vertical building, new high-density residential areas, and new industrial location are the main task of city planners for the next city plan (2030-2050).

10 .
Allocation-based optimisation results of the (A) SD scenario; (B) GP scenario; and (C) BAU scenario.
Figure 10 suggested the allocation of future new LUs because of LU growth based on three scenarios that useful for planning.Maintaining the open space and park, locating the vertical building, new high-density residential areas, and new industrial location are the main task of city planners for the next city plan (2030-2050).

Table 1 .
Comp cicj and Dep cicj are compatibility and dependency scores, respectively, based on Table 1.For example, if an industrial class LU has an elementary-school-LU neighbour with distance lower than R 1 , it has compatibility score of vlow (0.041) based on Table 1 (shown in bold).Compatibility and Dependency Scores from Survey. A.

Table 2 .
Scenarios for LU allocation in 2030.

Table 3 .
LU growth projection with population as an external input.

Table 3 .
LU growth projection with population as an external input.