Integer programming for urban design

Integer programming for urban design

Accepted Manuscript Integer Programming for Urban Design Hao Hua, Ludger Hovestadt, Peng Tang, Biao Li PII: DOI: Reference: S0377-2217(18)30923-8 ht...

7MB Sizes 0 Downloads 15 Views

Accepted Manuscript

Integer Programming for Urban Design Hao Hua, Ludger Hovestadt, Peng Tang, Biao Li PII: DOI: Reference:

S0377-2217(18)30923-8 https://doi.org/10.1016/j.ejor.2018.10.055 EOR 15447

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

22 March 2018 16 October 2018 29 October 2018

Please cite this article as: Hao Hua, Ludger Hovestadt, Peng Tang, Biao Li, Integer Programming for Urban Design, European Journal of Operational Research (2018), doi: https://doi.org/10.1016/j.ejor.2018.10.055

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • The integer program maximizes the floor area in an urban site under the shadowfree rule or the n-hour sunlight rule. • The integer program fills a building volume with 3D room templates.

CR IP T

• The Special Ordered Sets are employed to encode translational symmetry in urban layout.

AC

CE

PT

ED

M

AN US

• Both the Steiner tree model and the coverage network are extended to create reasonable route network.

1

ACCEPTED MANUSCRIPT

Integer Programming for Urban Design Hao Huaa,b,∗, Ludger Hovestadtc , Peng Tanga,b , Biao Lia,b Laboratory of Urban and Architectural Heritage Conservation (Southeast University), Ministry of Education,China. b School of Architecture, Southeast University, 2 Sipailou, Nanjing 210096, China. c ETH, Z¨ urich

CR IP T

a Key

Abstract

ED

M

AN US

We present an integer program (IP) for urban design that 1) maximizes the floor area; 2) fills building volume with room templates; 3) encodes translational symmetry in urban layout; and 4) constructs economical urban routes. Regardless that integer programming is intensively studied in operational research (OR), its role in solving geometrical and topological problems in urban design was overlooked. Based on a regular grid, our 01 IP formulates the sunlight-gain rules, which give urban sites their shapes, especially for residential projects. With predefined plot templates, the gross floor area (volume) within a given site can be maximized under various sunlight requirements. Subsequently, the IP fills each building volume with 2D/3D room templates. Finally, an IP-based algorithm constructs routes that connect all plots and the site’s entrances to public transportation. Both the classical Steiner tree model and the latest coverage network model are extended to create reasonable routes. In addition, this work extends the concept of special ordered sets (SOS) to encode translational symmetry in urban layouts. Encoding layout symmetry can benefit from the solvers’ SOS2 mechanism in the Branchand-Bound search algorithm. The results indicate that traditional decision making for cities could be partially automated by IP and an abundance of valid solutions are available for designers.

CE

PT

Keywords: linear programming, urban design, special ordered set, translational symmetry, Steiner tree

1. Introduction

AC

Making cities efficient, sustainable, and livable has become a critical social issue. Governments have been encouraging and funding research projects to design better urban environments with state-of-art technologies. For example, the Space Syntax Network in the UK, the computer graphics team at Purdue University, and the Future City Laboratory as a Singapore-ETH Z¨ urich corporation strove to apply computational intelligence to design problems and decision making arising in cities. Sophisticated evaluation models (e.g. foot traffic simulation and evaluation), interactive procedural modeling (Smelik ∗ corresponding

author Email address: [email protected] (Hao Hua) Preprint submitted to European Journal of Operational Research

November 1, 2018

ACCEPTED MANUSCRIPT

et al., 2014; Talton et al., 2011), inverse modeling and reconstruction (Aliaga et al., 2008; Bokeloh et al., 2010), and layout synthesis (Merrell et al., 2010; Rodrigues et al., 2013) have lately inspired many designers and architects. 1.1. Urban design

PT

ED

M

AN US

CR IP T

Urban design arranges building groups, street networks and public spaces. Urban designers focus on the layout of urban sites, as a composition of buildings, terraces, routes, green spaces, etc. Urban design is at the interface between architectural design and urban planning but is distinct form both disciplines (Moughtin, 2007). Urban planning emphasizes on large-scale urban structures, such as regional land use and transportation, and is less concerned with individual site layouts. Architectural design deals with individual buildings, usually after the building’s envelope is specified by urban design. Rowley (1994) pointed out that the results of urban design tend to be more relativistic and less deterministic than architecture, but more definite than planner’s prescription. Today’s urban design inevitably involves a multitude of rules, stipulations, and design principles. There are both aesthetic/idealism principles and performance-based principles. Schenk (2013) identified the qualitative traits (symmetry, gestalt cognitive laws, proportions and so on) and various measurable performances (functions of buildings, accesses to natural resources, circulation efficiency, etc.). This work deals with the following aspects of urban design: 1) gross floor area; 2) the percentages of distinct building/room types; 3) sunlight gain; 4) routes and circulation; and 5) symmetry in layout. The first four aspects are functional while the fifth aspect emphasizes geometric patterns. Architects have been enthusiastic in developing computational models for urban design, especially after the 2000s. For instance, the Architectural Association School of Architecture (AA School) in London pioneered parametric urbanism (Leach, 2009; Schumacher, 2009); the Berlage Institute at Rotterdam investigated associative design and synthetic vernacular (Trummer, 2008); and Computer-Aided Architectural Design (CAAD) at ETH Z¨ urich advanced self-organization methods in urban design (Hovestadt, 2009). However, to date, few works have been made for applying mathematical programming to urban design. 1.2. Integer programming

AC

CE

Linear programming (LP) is a well-developed technique for optimizing linear objective functions, subject to linear equality/inequality constraints (Dantzig, 2016). Given a problem, LP completely represents both the goal (objective function) and the constraints only with linear formula. In integer programming (IP) (Schrijver, 1998; J¨ unger et al., 2009), the decision variables are restricted to integers. This paper mostly involves 01 binary variables. The 0-1 IP is a thriving area of optimization because it is closely associated with combinatorial optimization (Wolsey & Rinaldi, 1993). The 0-1 IP is classified as NP-hard and belongs to Karp’s 21 NP-complete problems. Researchers have developed efficient algorithms such as Branch-and-Bound (B&B) algorithm (Conforti et al., 2014) and the cutting plane method for solving IP. Nowadays, cross-platform commercial solvers (Meindl & Templ, 2012) such as Gurobi and IBM CPLEX often employ both methods within a Mixed-Integer Programming (MIP) framework. Urban layout bears some similarity with the facility layout (Amaral, 2006; Drira et al., 2007; Friedrich et al., 2018), block layout for enterprises (Li & Smith, 2018), 3

ACCEPTED MANUSCRIPT

AN US

CR IP T

urban facility location (Farahani et al., 2010), the Very Large Scale Integration (VLSI) layout (Kim & Kim, 2003; Srinivasan et al., 2006). The facility layout problem (FLP) involves the arrangement of given facilities so that the total cost is minimized (Heragu & Kusiak, 1991). Either continuous variables (e.g., representing the relative position of each facility) or boolean variables (e.g, representing whether a facility occupies a specific cell in a regular grid) were employed to model FLP (Kothari & Ghosh, 2012; Amaral, 2006). In urban facility location problems, a major factor is the desirability of the facility (Coutinho-Rodrigues et al., 2012). Hammad et al. (2017) formulated a mixed integer non-linear programming model to minimize noise pollution and network congestion. Wu et al. (2018) developed an MIP-based model to make layout design for building interiors. The MIP has been a common solution approach to layout problems (Singh & Sharma, 2006; Amaral, 2008a,b; Hathhorn et al., 2013). Our model employs boolean variables and is solved by common MIP solvers. A latest breakthrough in the field of computer graphics is integrating IP with urban street networks and building floorplans based on a deformable grid (Peng et al., 2014, 2016).

ED

M

1.3. Special ordered set In 0-1 IP, the special ordered set (SOS) is an ordered set of variables satisfying particular constraints. This paper works with the special ordered sets of type 2 (SOS2): at most, two variables in the set can be non-zero, and if the two are non-zero they must be consecutive in the set. A historical motivation for developing SOS is to model piecewise linear functions (Dantzig, 1960; Beale & Forrest, 1976). However, we discovered a connection between SOS2 and translational symmetry in layout. We regard SOS2 as both a mathematical concept and a data structure. As an abstract model, it represents translational symmetry within the IP framework. As a data structure, SOS2 gives the B&B algorithm a cleverer way to handle the search procedure (Ernee Filho, 2014). Most commercial solvers have implemented the SOS2 constraints. 1.4. Contributions This paper’s contributions include:

AC

CE

PT

1. Using an IP to create urban layout that maximizes gross floor area, subject to desirable percentages of distinct building/room types and sunlight gain (shadowfree rule and n-hour sunlight rule). The IP also packs 2D/3D room templates into the building volumes. Our shadow-free model can be viewed as a variation of facility layout programs (Singh & Sharma, 2006; Amaral, 2008b) if the noon shadows are appropriately transformed into the template formulations. While our n-hour model is novel as it develops an OR logical treatment of shadows at different times during a day. 2. Using the IP to construct economical routes that connect all buildings and the site’s entrances to public transportation. Both the classical Steiner tree model and the latest coverage network model (Peng et al., 2016) are extended to create reasonable routes. The straightness of routes is also considered. Compared with the coverage network of Peng et al. (2016), our formulation involves fewer types of variables/constraints to achieve coverage and avoid islands. Peng et al. (2016) treated Zig-zag and T-junction respectively, while our model addresses the straightness of routes instead. 4

CR IP T

ACCEPTED MANUSCRIPT

M

AN US

Figure 1: Left: The T-shape template consists of four pixels. The bottom pixel is selected as the reference point. Right: The site P is in white and the void pixels (boundary) Q are in light grey. Two templates overlap at (4, 3). One template overlaps with the boundary at (1, 7).

ED

Figure 2: A plot template characterizes the extent of the building (red), the south rooms (green) that require sunlight, the shadows (dark grey), and ground terrace (cyan). The shadow pixels are determined via ray trace given the sun position at noon.

PT

3. Employing SOS2 to encode the translational symmetry in the layout. Cyclic SOS2 and composite SOS2 are defined to facilitate the pattern formulation.

AC

CE

Our method has a three-level scheme: 1) the peripheral coding such as specifying the plot templates or choosing a regular lattice (for example, square or triangular). Using different templates or grids will not change the IP formulation ; 2) the IP formulation; and 3) solving the IP with common solvers. Table 1 compares the results using Gurobi and IBM CPLEX. This paper mostly focuses on the IP and its references to the elements of urban design. 1.5. Paper structure Section 2 introduces the IP for shadow-free urban layouts. Various sunlight-gain rules are modeled with linear expressions. In Section 3, we investigate the symmetry in the layout where plot overlapping is allowed. We view such layouts as SOS2 layouts. In Section 4, we construct economical routes based on the Steiner tree and the coverage network. 5

AN US

CR IP T

ACCEPTED MANUSCRIPT

2. Shadow-free Layouts

M

Figure 3: Fee-shadow layouts constructed by the IP. The grey areas in the plans depict the noon shadows of the buildings. The envelopes of the buildings in the 3D view are predefined for each plot template. The square cell is 6 × 6 m (a-c). The side length of the triangular cell is 6 m (d-f).

PT

ED

A plot typically presents one building and its surrounding area (ground terrace, route, green space, etc.) (Meeda et al., 2007). A site usually has a network of roads or paths that subdivide the site into a set of plots (Lynch et al., 1984). So an essential step in urban design is to place the desirable plots inside a given site subject to legal constraints and restrictions. Some researchers regarded the road network as a major factor in shaping the urban structure (Samaniego & Moses, 2008; Peng et al., 2016), which will be elaborated in the section 3. 2.1. 2D/3D packing

AC

CE

An elementary principle of arranging multiple plots in a site is non-overlapping. In other words, the given plots should be compactly packed within the given site without having them overlap. Let P denote the site, namely the pixels of buildable area on a regular grid. A rectangular site could be represented by P = {(i, j)|1 ≤ i ≤ I, 1 ≤ j ≤ J} where i and j refer to the ith row and jth column of the site lattice, respectively. The set Q identifies the pixels that should be kept void. P ∩ Q = ∅ holds. The pixels in Q should enclose that in P , see Fig.1. A plot template located at (i, j) is represented by a set of pixels A(i, j), namely the plot’s footprint. One pixel is selected as a reference point, so A(i, j) = {(i, j), (i − 1, j − 1), (i − 1, j), (i − 1, j + 1)} describes the four-pixel template in Fig.1 (left). Each position in A refers to the reference point’s relative position to every pixel of the template. For instance, (i − 1, j) refers to the pixel above the reference point in the T-shape template(Fig.1, left), since the reference point’s relative position to this pixel is (−1, 0). 6

ACCEPTED MANUSCRIPT

Let binary variables xij , ∀ij ∈ P represent whether position (i, j) is occupied by the reference point of a template. Then, the layout problem is to maximize X xij (1) ij∈P

X

mn∈A(i,j)∩P

X

mn∈A(i,j)∩P

xmn = 0, ∀(i, j) ∈ Q xmn ≤ 1, ∀(i, j) ∈ P

CR IP T

subject to (2)

(3)

AN US

The first constraint keeps the void pixels empty while the second serves as the nonoverlapping principle. The parameter (i, j) in (2-3) could be a pixel from either P or Q, however, the subscript of xmn is only defined on P . Thus, the constraint (2) must address mn ∈ A(i, j) ∩ P . There are |P | + |Q| linear constraints. The boundary constraint (2) is violated at position (1,7) in Fig.1 because X xmn = x17 + x06 + x07 + x08 = x06 = 1 6= 0 mn∈A(1,7)∩P

mn∈A(4,3)∩P

M

The non-overlapping constraint (3) is violated at position (4,3) in Fig.1 because X xmn = x43 + x32 + x33 + x34 = x32 + x34 = 2

ED

Filling a building volume with 3D room templates can also be performed if the IP above replaces the pixels with voxels.

AC

CE

PT

2.2. Maximize floor area with shadow-free rule An urban site often comes across several types of plot templates. Inside a plot template, there are distinct parts such as buildings and (ground) terraces. Suppose there are L types of plot templates, we can characterize a plot template by Al (i, j) : the building’s footprint of the lth type plot, 1 ≤ l ≤ L. Bl (i, j) : the south rooms which require sunlight (for example, the bedroom). Cl (i, j) describes the building’s shadow at noon. The shadow’s dimension is conventionally reduced to the gap/height ratio stipulated by the local government. Dl (i, j) identifies the building’s additional extents, such as a ground terrace, courtyard, or garden. B ⊆ A, A ∩ C = ∅ (the building’s footprint does not intersect with the building’s shadow), A ∩ D = ∅ hold. The height (number of stories) of the building is specified so the shadows can be precalculated. The shadows of low buildings (no more than six stories) are often officially reduced to a gap/height ratio (Fig.2). For high buildings, one needs to calculate the exact hour of sunlight duration within a day (midwinter). The shape of shadows is a function of the sun’s position within the day. A building’s gross floor area equals the area of a single floor (the same as the footprint’s area in our cases) multiplying the number of 7

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 4: N -hour sunlight layout on a rectangular grid. The square cell is 7 × 7m. (a)(c): pre-calculate the building’s shadows at each hour for each type of plot template. Footprint A is marked in orange, south rooms B green, ground terrace D cyan. (b): The 3D room templates for the building in the plot template(a). (d): The 3D room templates for the building in the plot template(c). (e): STPSP constructs the walking routes after the layout IP optimizes the buildings’ locations. (f): The shadow conditions in the optimized layout. (g): 3D view of the layout.

ED

M

stories. Typically, wl in the objective function equals (or be proportional to) the number of stories of the lth type building so that the objective function (4) characterizes the gross floor area of the entire site. A common rule in arranging plots is that the south rooms which require sunlight should not be covered by the shadows of any other buildings. Thus, the IP is to maximize L X X wl xijl (4) l=1 ij∈P

AC

CE

PT

over binary variables xijl ∈ P . The xijl denotes whether the position (i, j) is occupied by the reference point of the lth template. subject to X xmnl = 0, ∀l, ∀(i, j) ∈ Q (5) mn∈Al (i,j)∩P

L X

X

l=1 mn∈Al (i,j)∩P

L X

X

l=1 mn∈[Bl (i,j)∪Cl (i,j)]∩P

xmnl ≤ 1, ∀(i, j) ∈ P

(6)

xmnl ≤ 1, ∀(i, j) ∈ P

(7)

ξl ≤

X

ij∈P

xijl ≤ τl , ∀l

(8)

There are L|Q| + 2|P | + L constraints. The constraints in (7) prevent the south rooms from overlapping with shadows. The constraint (8) specifies the minimum (and 8

ACCEPTED MANUSCRIPT

the maximum) number of each type of plot template. As a result, the number of the lth type plot will between ξl and τl in the generated layout. An alternative of (8) is to control the relative percentages. An additional constraint is xmnl +

mn∈Dl (i,j)∩P

L X

X

l=1 mn∈Al (i,j)∩P

xmnl ≤ 1, ∀l, ∀(i, j) ∈ P

(9)

CR IP T

X

which excludes (ground)terrace-building overlapping. A ground terrace is illustrated in Fig.2, sometimes a designer wants to prevent such terrace from occupying by other facilities. An alternative of (9) is to exclude terrace-terrace overlapping (if the template’s terrace encloses its building): L X

X

xmnl ≤ 1, ∀(i, j) ∈ P

(10)

AN US

l=1 mn∈Dl (i,j)∩P

M

The program in Fig.4 employs this constraint to prevent the ground terraces of adjacent plots from overlapping each other. Fig.3 illustrates two examples of the shadow-free layouts constructed by the IP. The grey areas in the plans (Fig.3(b)(d)) depict the noon shadows of the allocated buildings. Thus, the amount of white pixels indicates how much space is left after packing the plot templates (buildings). The site is 100% used if there is no white pixel in the generated layout. There is one white pixel in Fig(b), which indicates that the result is near optimum. The 3D envelopes of the buildings are predefined for each plot template. The choices on the lattices (rectangular, triangular, and others), the templates, and the sites can be separated from the IP.

AC

CE

PT

ED

2.3. N -hour sunlight rule A more precise way of evaluating the amount of sunlight gain is to trace the sun’s position within a day. Many countries have stipulated special requirements for sunlight availability in dense construction and for residential projects. For instance, Z¨ urich employs the so-called two-hour shadow rule, which only permits a building to cast shadows on residential buildings for (at most) two hours in a day (Hovestadt, 2009). A variation around the east coast of China is that the residual building must receive at least two-hour sunlight on the ‘worst’ day, that is, during the winter solstice. The general principle of such n-hour rules is to bind the shadowed hours below a threshold. The program first pre-calculates the building’s shadows at each hour (or at a smaller time interval) for each plot template, then constructs Clt (i, j) which describes the shadowed pixels at a specific time t, 1 ≤ t ≤ T . We used ray trace to determine whether a pixel is shadowed by the building at time t. First, the sun’s position in the sky at a particular time during a year follows a parametric model. Second, each face of the building’s 3D envelope is projected to the ground according to the sun position. Third, A pixel is counted as the lth building’s shadow (at time t) if it is inside any projected faces (at time t) from the building. Finally, all resultant shadowed pixels are saved as the constant arrays Clt (i, j) before the IP starts. The n-hour sunlight IP contains a set of binary variables ytij that tells whether position (i, j) is shadowed at time t. If the shadow at time t from any plot template 9

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 5: n-hour sunlight layout on a triangular grid. The side length of the triangular cell is 7 m. (a)(c): pre-calculate the building’s shadows at each hour for each type of plot template. Footprint A is marked in orange, south rooms B green, ground terrace D cyan. (e): STPSP constructs the walking paths after the layout IP optimizes the buildings’ locations. (f): The shadow conditions in the optimized layout. (g): 3D view of the layout.

covers the pixel (i, j), we assign 1 to ytij . This logical condition can be written as ytij = OR{xmn |mn ∈ C t (i, j) ∩ P }. Therefore, the IP maximizes the objective function in (4) subject to

M

X

PT

ED

mn∈A(i,j)∩P

X

mn∈A(i,j)∩P

X

mn∈D(i,j)∩P

xmn = 0, ∀(i, j) ∈ Q

(11)

xmn ≤ 1, ∀(i, j) ∈ P

(12)

xmn ≤ 1, ∀(i, j) ∈ P

(13)

CE

ytij = OR{xmn |mn ∈ C t (i, j) ∩ P } T X t=1

ytij + (T − α)

X

mn∈B(i,j)∩P

∀t, ∀(i, j) ∈ P

(14)

xmn ≤ T, ∀(i, j) ∈ P

(15)

AC

The logical OR could be transformed into linear inequalities X 0 ≤ M ytij − xmn ≤ M − 1, ∀t, ∀(i, j) ∈ P mn∈C t (i,j)∩P

where M = |C t |. Constraints in (14,15) mean that the shadowed hours must be below the threshold α if the position is occupied by a sunlight-demanding room. Imposing a slightly stronger constraint, one can replace (14,15) with 10

ACCEPTED MANUSCRIPT

mn∈C t (i,j)∩P T X

X

t=1 mn∈C t (i,j)∩P

xmn + (M T − α)

xmn ≤ M, ∀t, ∀(i, j) ∈ Q

(16)

xmn ≤ M T, ∀(i, j) ∈ P

(17)

X

mn∈B(i,j)∩P

CR IP T

X

P where 1 ≤ M ≤ |C t | (usually set M to an estimated value, e.g., 2). mn∈B(i,j) xmn is limited to 0 or 1 by extents-extents exclusion. Notice that the constraints in (16,17) no longer involves the auxiliary variable ytji . Therefore, the n-hour sunlight IP for multiple plot templates is to maximize L X X

wl xij

(18)

AN US

l=1 ij∈P

where the weight wl typically equals Hl |Al |. Hl denotes the number of stories in the lth plot template. subject to (5)(6)(10) and L X

X

L T X X

X

xmnl + (M T − α)

ED

t=1 l=1 mn∈Clt (i,j)∩P

M

l=1 mn∈Clt (i,j)∩P L X

X

xmnl ≤ M, ∀t, ∀(i, j) ∈ Q

(19)

xmnl ≤ M T, ∀(i, j) ∈ P

(20)

l=1 mn∈Bl (i,j)∩P

CE

PT

Constraints (19-20) are multi-type extension of constraints (16-17). Fig.4 and Fig.5 illustrate two examples of layout construction under the n-hour sunlight rule, on a rectangular grid and on a triangular grid, respectively. The illustrated results involve three IPs: the layout IP (n-hour sunlight rule), the 3D room templates packing IP, and the circulation IP (Section 4). The 3D packing and circulation are performed after the layout IP optimizes the buildings’ locations. 3. Symmetry and SOS2

AC

Symmetry is a constant topic of urban and architectural design in both theory and in practice. Schenk (2013) regarded the symmetry as an essential quality when evaluating or designing an urban area. The urban symmetry has a powerful psychology effect and was thought invincible in reality (Lynch, 1984). The symmetry may reinforce the organization of government, the disposition of social ranks, and the behavior of citizens. In terms of modeling techniques, shape grammar, L-system, and procedural models have dominated the landscape of symmetry modeling (Aliaga et al., 2008; Bokeloh et al., 2010; Hua, 2017). However, we discovered a “hidden” mechanism to encode translational symmetry in SOS2-based IP. What we need to do is to manipulate the order in which the elements (pixels) of a template are presented to the IP. 11

CR IP T

ACCEPTED MANUSCRIPT

M

AN US

Figure 6: Left: the ordered elements of the Knight8 template. Within its SOS2 layout, any overlapped pair of templates must be associated with a horse-jump translation. Right: the ordered elements of the Lemniscate16 template. Within its CS2 layout, any overlapped pair of templates must be associated with a two-step diagonal or vertical translations.

ED

Figure 7: Urban layouts constructed with Knight8 (left) and Hexa9 (right). The vehicle route network is constructed via CST before the arrangement of plot templates.

AC

CE

PT

Modeling urban layouts with SOS2 allows two plots to overlap. As a result, a position (pixel) in the site can be 1) unoccupied, 2) occupied once, and 3) occupied twice by the plot templates. How the three situations correspond to the presence of building, (ground) terrace, or other spaces is left to who designs the plot templates. Two examples of such mapping are given in Fig.7. Since the SOS2 constraint only allows two consecutive binary variables to be 1, the rules for overlapping templates are implicitly encoded by the order in which the template’s pixels are listed in A. In contrast to the IP in Section 2 which is invariant to the predefined order of elements of a template, such order defined in SOS2 IP informs the displacements of plot templates in the layout. As a result, bringing a specific symmetry into the layout can be achieved at no additional expense. 3.1. Generalizing SOS2 In IP formalism, one role of SOS2 is addressing constraints: no more than two variables in the set may be non-zero, with the further condition that if there are two they must be adjacent (Beale & Forrest, 1976; Ernee Filho, 2014). Another role is offering a special data structure that facilitates the B&B search algorithm. This work focuses on 12

ACCEPTED MANUSCRIPT

AN US

CR IP T

the confluence of the first role with translational symmetry, as most MIP solvers have implemented SOS2 algorithms. The SOS2 model was developed for modeling piecewise linear functions. One way to formulate nonlinear functions is by adding binary variables and new inequalities to the model, yielding a mixed-integer program (MIP) (Keha et al., 2004). The well-known MIP formulations of piecewise linear functions include the incremental model and the convex combination model (Croxton et al., 2003; Vielma et al., 2010). The SOS2 formulation is similar to convex combination, however, the nonlinearity is implemented in the B&B algorithm of solver. Given a set of binary variables {xi }, the notation SOS2{xi } means that the set {xi } satisfies the SOS2 constraints. In terms of conventional linear expressions, SOS2{xi } means (Dantzig, 1960)  x1 ≤ y1      x2 ≤ y1 + y2    x ≤ y i i−1 + yi  · · ·     xI ≤ yI−1   PI−1 yi = 1 i

M

where yi are auxiliary binary variables. It is natural to introduce a cyclic version CS2{xi } ( xi ≤ yi + y(i+1)%I , i = 1, 2, · · · , I PI i yi = 1

ED

where % denotes the modulo operation. I is the set’s cardinality. One can further extend the SOS2 concept to a set of sets. CS2{{x11 , x12 , · · · }, {x21 , x22 , · · · }, · · · }, namely composite cyclic SOS2, imposes the following constraints ( xij ≤ yij + yi,(j+1)%J , ∀i, j P P i j yij = 1

PT

where yij are binary variables. Here we assume that the cardinalities the of subsets are the same (equal J). SOS2{{}} has the same structure as CS2{{}} except for the cyclic feature. The classical SOS2{xi } is a special case of composite SOS2{{xi }}.

AC

CE

3.2. Encoding translational symmetry The SOS2 mechanism and the user-defined order of elements in a template lead to an implicit formulation of translational symmetry. Take the eight-element template in Fig. 6 (top left) for example. We arrange the elements in the sequence of the knight move (as in chess). We named this plot template with ‘Knight8’: A(i, j) = {(i, j − 1), (i − 1, j + 1), (i + 1, j), (i − 1, j − 1), (i, j + 1), (i + 1, j − 1), (i − 1, j), (i + 1, j + 1)}

The SOS2 layout program maximizes X

xij

ij∈P

13

(21)

ACCEPTED MANUSCRIPT

subject to X

xmn = 0,

∀(i, j) ∈ Q

(22)

SOS2{xmn |mn ∈ A(i, j)},

∀(i, j) ∈ P

(23)

mn∈A(i,j)∩P

CR IP T

If (m, n) ∈ / P occurs for some (m, n) ∈ A(i, j), one can place a free binary variable at the corresponding position in the ordered set. This program yields the layout in Fig.7 (left) with the MIP solver’s built-in SOS2 function. We can see that the relative displacement between every pair of overlapping templates complies with the horse-jump pattern, because 1. Only two consecutive positions can be simultaneously activated (which causes overlapping) by SOS2. 2. All the consecutive positions are predefined in a horse-jump sequence.

AN US

The following example explains the horse-jump mechanism. Let (i, j) = (0, 0) in (23), then the SOS2 constraint says that only the consecutive pair in {x0,−1 , x−1,1 , x1,0 , x−1,−1 , x0,1 , x1,−1 , x−1,0 , x1,1 } can equal 1. All the possible pairs are

M

x0,−1 = x−1,1 = 1. x−1,1 = x1,0 = 1. x1,0 = x−1,−1 = 1. x−1,−1 = x0,1 = 1. x0,1 = x1,−1 = 1. x1,−1 = x−1,0 = 1. x−1,0 = x1,1 = 1.

ED

1. 2. 3. 4. 5. 6. 7.

AC

CE

PT

The relative displacement in each pair follows the horse-jump pattern, e.g. (0, −1) against (−1, 1), or (−1, 1) against (1, 0). One can see that the resultant displacement between overlapping templates is a function of the predefined order of the elements in the template. If the horse-jump sequence is a closed tour, then we should use CS2 (cyclic) constraints instead of the standard SOS2. How to find such a knight-move sequence that exactly fills an arbitrary pixelated template? This is the classic Knight’s tour problem. The Lemniscate16 template (Fig.6 top right) illustrates the application of the composite SOS2. The template consists of 16 pixels organized into four groups. Lemniscate16 is defined by Ak (i, j) = {(i − p, j − q), (i − 2 − p, j − q),

(i − p, j − 2 − q), (i − 2 − p, j − 2 − q)}, p = (k − 1)%2, q = (k − 1)/2

Within each subset Ak (i, j) of four pixels, the sequence follows a lemniscate-shape pattern. Then the layout program employs CS2{{xmn |mn ∈ A1 (i, j)}}, {xmn |mn ∈ A2 (i, j)}, · · · }, 14

∀(i, j) ∈ P

ACCEPTED MANUSCRIPT

CR IP T

which allows two templates to overlap by a four-pixel square (Fig.6 bottom right) with either a diagonal or a vertical translation. The SOS2 encoded symmetry allows for a certain degree of freedom in regulating the pattern of covering. Firstly, a plot’s position is not SOS2 constrained if it does not overlap others. Secondly, a SOS2 structure may prescribe multiple transitional rules for overlapped plots (for example, the Lemniscate16 permits relative translations at 45◦ , 90◦ , and 135◦ ). Finally, what the overlapped parts mean is an essential part of the template design. The layout automation program still leaves much of the design decisions up to designers. 4. Route network

ED

M

AN US

The circulation between the plots (buildings) is planned before or after the layout composition is computed by the IP in Section 2. The grid P , excluding the parts occupied by buildings, presents a graph from which the routes will be identified to connect all buildings (terminals). This task raises classical topics such as the minimum spanning tree (MST) and the network flow problems. Another significant development, namely network design from functional specifications (or coverage network), was recently proposed by Peng et al. (2016). We treat both the Steiner tree and the coverage network as a flow formulation of network. The major difference between the Steiner tree and the coverage network is that the former connects a set of predefined destinations (nodes) while the latter makes any node within a certain range of the routes. They are suited to different scenarios in urban design. An urban designer is usually concerned with two practical criteria. First, the routes should allow the inhabitants to conveniently walk from one building to other buildings (interior pairwise connections). Second, the routes should provide easy access to urban traffic by reaching certain boundaries of the site (interior-exterior connections). Differing from classical graph-based problems, a desirable circulation demands both topological and geometric features:

CE

PT

1. The walking path network must connect all buildings and certain predefined boundaries of the site. 2. The routes should be as straight as possible. This is equivalent to minimizing the number of the routes’ turns on a regular grid. This geometric property is beyond the scope of graph theory.

AC

First, a graph (V, E) is constructed from the grid P (see examples in Fig.4(e) and Fig.5(e)). People can either regard the cells (pixels) as the nodes (like a chess board) or view the intersection points as the nodes (like the Go game). Given an undirected graph (V, E) and a subset of K nodes as the terminals, a flow formulation of the route network typically includes: |E| positive weights wij , denotes the Euclidean length of the edge (i, j). Here we let wij = 1, ∀(i, j) ∈ E. |E| binary variables zij , undirected, denotes whether the edge (i, j) is in the route network. 2|E| binary variables yij , directed, denotes whether the directed edge (i, j) is in the route network. 15

ACCEPTED MANUSCRIPT

yij + yji = zij , ∀(i, j) ∈ E which

CR IP T

2|E| integer variables dij (for coverage network), directed, denotes the distance from the edge (i, j) to a terminal. The lower bound of dij is 0, the upper bound is the longest among the shortest distances between any pair of nodes. k 2K|E| binary variables fij (for Steiner tree), denotes the flow of commodity k from ith node to jth node. Notice that zij corresponds to the original graph (undirected) while yij renders all edges as bi-directional. Our flow formalism imposes (24)

AN US

1. prevents two flow goes through an edge with opposite directions, because yij = yji = 1 leads to yij + yji = 2 6= zij (notice that zij ≤ 1 ). 2. imposes zij = 1 when either yij or yji equals 1, i.e., makes an undirected edge active if there is any flow goes through it 4.1. Straightness

ED

M

Urban design prefers straight paths because people do not like detours - whether a segment of a route is straight on a rectangular grid can be identified by znI + znIII and znII + znIV (see Fig 8, bottom left). zij refers to an edge by the indices of the two end nodes; while znI , znII , . . . refers to an edge by a node and the directions. The IP program does not create variables for znI , znII , . . . , as they refer to the same set of variables {zij }. If znI + znIII equals 1, there is a right-angle turn or a dead end; otherwise (0 or 2), the path is straight or there is no path. It is the same for the term znII + znIV . To count the number of turns in the tree, the IP introduces the binary variables ( hn = znI XOR znIII ∀n ∈ V vn = znII XOR znIV

PT

If the solver does not support the XOR logical constraint, one can use the following equalities instead ( znI + znIII = hn + 2mn ∀n ∈ V (25) znII + znIV = vn + 2on

AC

CE

where mn , on are auxiliary binary variables. The triangular grid pertains to a similar formulation with a slightly different index method (see Fig.8, top left). P To impose straightness, one can add a term n∈V (hn + vn ) in the objective function. This technique is employed in our extensions of the Steiner tree and the coverage network. 4.2. Steiner tree and extensions Given an undirected graph and a subset of vertices called terminals, the Steiner tree is a minimum weight subgraph that contains all terminals. We choose an arbitrary terminal as the source and other terminals as sinks. The terminals are arranged as the first K + 1 nodes of the graph. Index 0 denotes the source, then the sinks are indexed with k, 1 ≤ k ≤ K. In an urban layout, one building is the source, while other buildings 16

CR IP T

ACCEPTED MANUSCRIPT

Figure 8: Both Steiner trees have the same length while the one on the right prefers straight routes. The source is marked in green and sinks in blue.

AN US

(and the site’s entrances) serve as sinks. One can imagine there are a group of people walking from the source towards each sink. The total foot travelling distance should be minimized. The flow formulation renders the edges bi-directional. If E ∗ denotes the set of directed edges, |E ∗ | = 2|E| holds. The multicommodity flow IP (Chopra & Tsai, 2001; Polzin & Daneshmand, 2001) minimizes X wij zij (26) ij∈E

subject to (24) and

ED

ij∈E

 −1 if j = 0   k k fij − fji = 1 if j = k   0 otherwise

M

X

k fij ≤ yij

∀k, ∀ij ∈ E ∗

∀k, ∀j ∈ V

(27)

(28)

CE

PT

(27) presents K|V | constraints. (28) presents 2K|E| constraints (edges are bi-directed). (24) presents |E| constraints. The straight Steiner tree minimizes X X zij + s (hn + vn ) (29) ij∈E

n∈V

AC

subject to (25). s is a constant weight for imposing straightness. Source-sink distances. The Steiner tree characterizes the length of the entire tree but ignores the lengths of pairwise paths. As the distance from the source to the kth k k sink is simply fij + fji , one strategy for controlling source-sink distances is 1) choosing one terminal around the site’s center as the source, and 2) imposing X  k k fij + fji ≤ Mk , ∀k (30) ij∈E

where Mk is a positive threshold for each sink. If the threshold is too small, the IP is infeasible; if the threshold is too big, this IP is reduced to the standard Steiner tree. The 17

ACCEPTED MANUSCRIPT

AN US

CR IP T

consequence of regulating source-sink distances are illustrated in Fig.9. Combining the constraints on straightness and on source-sink distances leads to more reasonable routes in urban layouts. Steiner tree with pairwise shortest paths (STPSP). An essential property of the Steiner tree is that it excludes cycles. Our program transforms certain dead ends (the leaves of the Steiner tree) into cycles by adding pairwise shortest paths after the Steiner tree is constructed. The program identifies a pair of terminals whose Euclidean distance is small but the travelling distance through the tree is large. The added paths should coincide with the Steiner tree as much as possible, so the number of added segments is minimized. Given a Steiner tree T as a set of edges, the pseudo-codes of adding pairwise paths are 1. Construct a set of node pairs {(i, j)} where dij > lij β. lij denotes the Euclidean distance between the ith and jth nodes, while dij denotes the length of i − j path through the tree. β is a constant threshold. 2. Sort the pairs by dij in descending order. 3. For each pair in {(i, j)}, Do a. run a shortest path IP to find a short circuit between the ith and jth nodes. b. add the new path into T (the resulting edge network). The shortest path IP is a special case of the Steiner tree IP with one source and one sink. Given a Steiner tree T , the pairwise IP minimizes X X X zij − t zij + s (hn + vn ) (31) ij∈T

M

ij∈E

n∈V

AC

CE

PT

ED

subject to (24,25,27,28,30). where t is the constant weight for the coincidence between the new path and the given tree T . A typical value of t is 0.1. The second term tends to overlap the added path with the existing Steiner tree. An illustration of adding pairwise paths are shown in Fig.9 (d). The routes in Fig.4 and Fig.5 are also constructed by the STPSP method. Coupled Steiner trees (CST). If there are two flows travelling through distinct paths from the source to a sink, then there is a cycle through the source and the sink. Based on this insight, we can construct a couple of Steiner trees to create cycles instead of dead-end branches. We termed this model with Coupled Steiner trees (CST). To distinguish the two trees in the formulation, we denote the variables of the male tree † with zij and the female with zij . The CST program minimizes X X † (zij + zij − ρcij ) + s (hn + vn ) (32) ij∈E

n∈V

subject to (24,25,27,28,30) and the auxiliary variables cij satisfy † cij = zij ∧ zij ,

∀ij ∈ E

(33)

The first term in the objective function governs the coincidence of the two trees. A negative ρ enforces coincidences between the two trees, while a positive ρ will reduce coincidences. Two examples are shown in Fig 10. One tree is marked in black, the other in red. If ρ =0, the program results in two independent Steiner trees. 18

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

ED

M

Figure 9: The Steiner trees in (a) and (b) have the same length while (b) keeps the source-sink distances below a threshold so that the distance from the source to A and B are significantly reduced. The tree in (c) controls straightness and source-sink distances. Adding pairwise shortest paths (C to C’, D to D’, E to E’) leads to the network in (d).

AC

Figure 10: Two series of Coupled Steiner trees (CST) tests. One tree is marked in black, the other in red. Negative ρ enforces coincidences between the two trees, while positive ρ prevents coincidences.

4.3. Coverage network The coverage network takes a distinct approach to networking planning. It constructs a connected subgraph which reach all nodes within a predefined range (Peng et al., 2016). So the network can evenly cover the entire site even with a single terminal. The objective prefers a network with 1) smaller total length; and 2) short distances to the terminals (similar to the constraints in (30) of our Steiner tree model). Given an undirected graph (V, E) and the terminals indexed as the first K nodes, the 19

CR IP T

ACCEPTED MANUSCRIPT

coverage network is to minimize X

ij∈E

AN US

Figure 11: Coverage networks with distinct coverage ranges. Top row: without straightness objective. Bottom row: with straightness objective.

zij + ρ

X

dij

M

where E ∗ denotes the set of directed edges. subject to (24) and

IF yij = 1, dij − djk ≥ 1, ∀ij ∈ E ∗ ∧ j > K X yjk ≥ yij , ∀ij ∈ E ∗ ∧ j > K X yij ≥ 1, ∀l ∈ V

ED

(34)

ij∈E ∗

(35) (36) (37)

ij∈N (l)

AC

CE

PT

where (35) and (36) present roughly 2|E| constraints, receptively. While (37) presents V constraints. The IF condition in (35) can be directly implemented by the indicator function of MIP solvers, or transformed into inequality by big-M method (Peng et al., 2016). If the edges are not of uniform length, one should P replace the inequality dij − djk ≥ 1 with d − d ≥ w and replace the term ij jk ij ij∈E zij in the objective function with P w z . ij ij ij∈E The constraints in (36) mean that at least one of outward edges of the jth node (except for the terminals) should be active if there is an active inward edge (i, j). In other words, a non-terminal node should allow a flow travels through it. The set N (l) denotes the edges which should be covered by the lth node. So (37) requires that the topological distance of any node to active edges should be within a threshold. This is why we named this model with coverage network. The value of the coverage range (distance threshold) is critical (see Fig.11). Compared with original formulation of Peng et al. (2016), our program involves fewer variables (the original involves five types of variables while our program has three: zij , yij 20

ACCEPTED MANUSCRIPT

and dij ) and fewer constraints (the original has six kinds of constraints while our program has four) to achieve 1) the no-island goal, and 2) the coverage goal. To make the routes as straight as possible, one can use constraints in (25) and minimize X X X dij + s (hn + vn ) (38) zij + ρ ij∈E ∗

n∈V

CR IP T

ij∈E

4.4. Comparison

5. Results and Discussions

M

AN US

Both the coverage network and the Steiner tree have been treated as a flow problem defined in a graph. The coverage network uses a distance indicator dij (integer) for each directed edge to achieve global connectivity (i.e., no islands) while Steiner tree k employs a binary variable fij for each directed edge (and for each sink) to achieve global connectivity . The coverage network is more suitable for the design tasks where the goal is to roughly cover the entire site and few specific terminals are known. In this situation, the designer assumes that the buildings can reach the route network via ground terraces within a distance threshold. Comparing STPSP with CST, the former is more flexible since the designer may choose particular pairs to add the short circuits, while the latter systematically create cycles. In terms of the urban/street circulations, STPSP which probably still contain dead ends is more suitable for pedestrians, while CST which facilitates round-trips is more appropriate for automobiles.

PT

ED

We proposed an IP treatment of urban design, including shadow-free layout synthesis, room templates allocation, and circulation planning. Various design criteria and restrictions are completely modeled with linear expressions. Our method consists of three levels: 1) specification of the site and the plot templates; 2) the IP formalism; and 3) IP solving by professional solvers. Hence, a particular urban design can be implemented without altering the IP. 5.1. Applications

AC

CE

The increasing conviction of AI now makes the automation of urban/architectural design an eager expectation among architects. Our IP solves a number of common problems in urban design, such as maximizing floor area under sunlight demand, endowing the urban layout with symmetric patterns, filling the building volumes with room templates, and planning interior routes. Manual arrangement of these elements requires expertise and is time-consuming. The parameterization and automation of these tasks can quickly report the statistics of a given site and visualize the near-optimal layouts. The interaction between investors, governments, inhabitants, and designers could be reshaped by the automation programs. The IP of n-hour sunlight rule fits residential projects that emphasize on the sunlight gain of individual rooms. We can also define proper ground terrace for each plot template (working with constraints in (9) or (10)) to regulate the distances between generaltype buildings in high-density regions. The symmetry-encoding mechanism based on 21

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 12: Left: comparison of the density of empirical data from Shanghai (the maximum and minimum curves in the figure) and that of the generated layouts. Right: comparison of the FAR of empirical data from Shanghai and that of the generated layouts. Curve a-d: shadow-free layouts. Curve e-k: 2-hour sunlight layouts.

SOS2 allows designers to manipulate translational rules when allocating plots/buildings. Although shape grammar and procedural models have been immensely employed for modeling geometric symmetry (Bokeloh et al., 2010; Hua, 2017), this work presents an alternative approach using pure 0-1 IP. Cheng et al. (2006) investigated the complex relationship between high density, urban form, and solar potential. The randomness was elaborated:

ED

M

1. Horizontal randomness (against uniform distribution of buildings in the layout) is preferable. Our generated layouts (Fig.3,4,5,7,13) are non-uniform. 2. Vertical randomness (buildings with different heights) is preferable. This is implemented in our IPs by employing plot templates with buildings of different heights (Fig.3,4,5,7,13).

AC

CE

PT

We compared the statistics of our generated layouts with the empirical data from a survey of residential buildings in Shanghai (YANG & CHEN, 2005). Fig.12 (left) compares the density and Fig.12 (right) compares the floor area ratio (FAR), namely the ratio of total floor area to the site area. The shadow-free IP is tested for the low buildings (2-story to 6-story) in a 45288 m2 site in Shanghai, using uniform building with footprint of 114 m2 (curve a), 216 m2 (b), 288 m2 (c), and 324 m2 (d), respectively. We also tested the 2-hour IP for the high buildings (6-story to 18-story) in a 49980 m2 site in Shanghai, using 294 m2 (e), 392 m2 (f), 441 m2 (g), 588 m2 (h), and 784 m2 (k) footprint, respectively. Fig.12 (left) shows that the densities of the most generated layouts are between the maximum and the minimum values from the real world. Fig.12 (right) indicates that 1) the FAR of the generated layouts are around the empirical maximum, probably because the IPs maximize the FAR while the real projects may choose a medium value of FAR to reconcile other criteria; and 2) The larger the building’s footprint is, the higher the FAR is. The designer may manipulate the settings of the plot template to enhance the design space for a given project. For instance in the shadow-free program (section 2.2), one can design a template as the void region (as a sunken courtyard) on the basement level. And 22

ACCEPTED MANUSCRIPT

nodes 1 1 164 247 17568 3411 3311 2137 43 286 607554 3361 3538702 254 11 158 1 8 1 1

gap(%) 0 0 0.2 0 2.5 5.9 27.6 29.9 48.5 49.3 199 0 130 2 0.8 0 1.9 1.1 0.1 0

value 198 198 2473.5 2473.5 1557.9 1497.9 111 111 57882 57882 101 105 97 105 142.8 142.8 151.3 150.3 179.1 179

solver CPLEX ∗ Gurobi CPLEX ∗ Gurobi CPLEX ∗ Gurobi ∗ CPLEX Gurobi ∗ CPLEX Gurobi ∗ CPLEX SOS2 CPLEX linear Gurobi SOS2 Gurobi linear CPLEX ∗ Gurobi CPLEX ∗ Gurobi CPLEX ∗ Gurobi

CR IP T

time(s) <1 <1 5 39 180 180 30 120 30 60 60 60 60 60 180 74 600 215 230 98

AN US

Program packing room templates (Fig. 4) packing room templates (Fig. 4) shadow-free layout (Fig. 3 left) shadow-free layout (Fig. 3 left) shadow-free layout (Fig. 3 right) shadow-free layout (Fig. 3 right) n-hour sunlight layout (Fig. 4) n-hour sunlight layout (Fig. 4) n-hour sunlight layout (Fig. 5) n-hour sunlight layout (Fig. 5) SOS2 layout (Fig. 7 left) SOS2 layout (Fig. 7 left) SOS2 layout (Fig. 7 left) SOS2 layout (Fig. 7 left) Steiner tree (Fig. 4 (e)) Steiner tree (Fig. 4 (e)) Steiner tree (Fig. 5 (e)) Steiner tree (Fig. 5 (e)) CST (Fig. 7 left) CST (Fig. 7 left)

M

Table 1: The timing, number of B&B nodes, gap, and value of solving IPs with commercial solvers (Gurobi & IBM CPLEX). The ∗ symbol means that the solution is illustrated in the image.

ED

subsequently these sunken courtyards can interact with the buildings above the ground. Fig.13 shows such an application where the sunken courtyards do not overlap with the buildings’ footprints. 5.2. Implementation and timings

AC

CE

PT

One appeal of modeling urban design problems with IP lies in the clear separation of formulation and implementation. The latest professional solvers usually support a multitude of computer platforms and programming languages, which makes the implementation of IP easy. We implemented our programs using Gurobi and IBM CPLEX through their Java APIs and often adopted the default configuration of their MIP functionality. The solvers take seconds or minutes to find the optimal solution for many programs presented in this paper. For larger-scale IPs, we opt for deriving a relatively good solution within minutes. Table 1 lists the timing of some cases, on a 64-bit Windows computer with 3.6GHz i7-4790 CPU and 16GB RAM. The performances of Gurobi and IBM CPLEX seem similar. Both solvers support SOS2 constraints. The performances of using SOS2 and linear inequalities are similar. 5.3. Further work All models proposed in this paper are based on regular grid, which is a fundamental limitation of our IP approaches to urban design. Although many current projects of urban design still commence with a regular grid, our future work will integrate deformed grid (Peng et al., 2016). It would be interesting to consider the inverse engineering 23

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

M

Figure 13: Layout with sunken courtyards. Four terminals (marked with green dots) are pre-defined on the site’s boundary. The sunken courtyards (i.e., void space on the basement level) marked in green do not overlap with buildings’ footprints. The buildings’ shadows are marked in grey in the plan. The coverage network program (coverage range = 2) performs after the buildings and the courtyards are optimized.

AC

CE

PT

of exiting cities based on IP. Successful works on inverse modeling (Pauly et al., 2008; ˇ St’ava et al., 2010; Kalogerakis & S. Chaudhuri, 2012) train the parameter values of a parametric model from given samples. The inverse modeling of street networks (Aliaga et al., 2008; Fiser et al., 2016) is probably most close to this purpose. However, encoding the hidden rules in terms of IP (or other models of mathematical programming) is still very much a work in progress. Another future work is integrating GPU parallel computing with algorithms for solving 1-0 IP. Classical algorithms such as B&B were not tailored for special hardware, such as GPUs or quantum computers. The parallelization of B&B algorithms (Borisenko et al., 2017; Shen et al., 2017) and simplex type algorithms (Ploskas & Samaras, 2015) has recently become advanced. It would be rewarding to study the parallelization of SOS2 constraints. The SOS-based operations that are native to the hardware are expected. 6. Conclusions This paper demonstrates that a series of tasks in urban design can be fully modeled by IP and subsequently solved by state-of-art solvers. The multitude of design criteria, stipulations, and physical rules have long delayed the mathematical programming of urban design. However, our work indicates that both functional objectives (e.g. n-hour sunlight 24

ACCEPTED MANUSCRIPT

CR IP T

rule) and aesthetic principles (e.g. translational symmetry) could be transformed into linear expressions. Unlike many existing approaches that mix mathematical programming with other optimization methods (Singh & Sharma, 2006), our approach mostly employs 0-1 IP to handle various layout problems. Eventually, current MIP solvers will allow us to obtain near-optimal layouts for middle-scale sites (about 30,000 m2 ) within seconds or minutes. We employ the plot template as the essential element in layout synthesis. A plot can represent a building and its neighborhood; the overlap between two plots may represent a building or a terrace. Such flexibility constitutes a design space outside of the IP formalism. Therefore, the designers still play an essential role in devising the plot template. The interplay between a designer and the IP may yield novel results without modifying the IP algorithm. Acknowledgments

AN US

Thank Zifeng Guo at ETH Z¨ urich for the introduction to Gurobi. This work was supported by the National Natural Science Foundation of China [51778118, 51478116, 51578123]; the Ministry of Housing and Urban-Rural Development of China [UDC2017020212]. References

AC

CE

PT

ED

M

Aliaga, D. G., Vanegas, C. A., & Benes, B. (2008). Interactive example-based urban layout synthesis. In ACM transactions on graphics (TOG) (p. 160). ACM volume 27. Amaral, A. R. S. (2006). On the exact solution of a facility layout problem. European Journal of Operational Research, 173 , 508–518. Amaral, A. R. S. (2008a). An exact approach to the one-dimensional facility layout problem. Operations Research, 56 , 1026–1033. Amaral, A. R. S. (2008b). A MIP model for the layout problems with departments of unequal area. In XL SBPO Simp´ osio Brasileiro de Pesquisa Operacional (pp. 2017–2028). Beale, E., & Forrest, J. J. (1976). Global optimization using special ordered sets. Mathematical Programming, 10 , 52–69. Bokeloh, M., Wand, M., & Seidel, H.-P. (2010). A connection between partial symmetry and inverse procedural modeling. In ACM transactions on graphics (TOG) (p. 104). ACM volume 29. Borisenko, A., Haidl, M., & Gorlatch, S. (2017). A GPU parallelization of branch-and-bound for multiproduct batch plants optimization. The Journal of Supercomputing, 73 , 639–651. Cheng, V., Steemers, K., Montavon, M., & Compagnon, R. (2006). Urban form, density and solar potential. In PLEA 2006 LESO-PB-CONF-2006-008. Chopra, S., & Tsai, C.-Y. (2001). Polyhedral approaches for the steiner tree problem on graphs. In Steiner Trees in Industry (pp. 175–201). Springer. Conforti, M., Cornuejols, G., & Zambelli, G. (2014). Integer programming. Springer. Coutinho-Rodrigues, J., Tralh˜ ao, L., & Al¸cada-Almeida, L. (2012). A bi-objective modeling approach applied to an urban semi-desirable facility location problem. European Journal of Operational Research, 223 , 203–213. Croxton, K. L., Gendron, B., & Magnanti, T. L. (2003). A comparison of mixed-integer programming models for nonconvex piecewise linear cost minimization problems. Management Science, 49 , 1268– 1273. Dantzig, G. (2016). Linear programming and extensions. Princeton university press. Dantzig, G. B. (1960). On the significance of solving linear programming problems with some integer variables. Econometrica, Journal of the Econometric Society, (pp. 30–44). Drira, A., Pierreval, H., & Hajri-Gabouj, S. (2007). Facility layout problems: A survey. Annual reviews in control, 31 , 255–267. Ernee Filho, K. (2014). Valid inequalities and computational results for SOS1-, SOS2-, and cardinalityconstrained linear programs. Ph.D. thesis Texas Tech University.

25

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Farahani, R. Z., SteadieSeifi, M., & Asgari, N. (2010). Multiple criteria facility location problems: A survey. Applied Mathematical Modelling, 34 , 1689–1709. Fiser, M., Benes, B., Galicia, J. G., Abdul-Massih, M., Aliaga, D. G., & Krs, V. (2016). Learning geometric graph grammars. In Proceedings of the 32nd Spring Conference on Computer Graphics (pp. 7–15). ACM. Friedrich, C., Klausnitzer, A., & Lasch, R. (2018). Integrated slicing tree approach for solving the facility layout problem with input and output locations based on contour distance. European Journal of Operational Research, 270 , 837–851. Hammad, A. W., Akbarnezhad, A., & Rey, D. (2017). Sustainable urban facility location: Minimising noise pollution and network congestion. Transportation Research Part E: Logistics and Transportation Review , 107 , 38–59. Hathhorn, J., Sisikoglu, E., & Sir, M. Y. (2013). A multi-objective mixed-integer programming model for a multi-floor facility layout. International Journal of Production Research, 51 , 4223–4239. Heragu, S. S., & Kusiak, A. (1991). Efficient models for the facility layout problem. European Journal of Operational Research, 53 , 1–13. Hovestadt, L. (2009). Beyond the Grid–Architecture and Information Technology: Applications of a Digital Architectonic. Walter de Gruyter. Hua, H. (2017). A bi-directional procedural model for architectural design. In Computer Graphics Forum (pp. 219–231). Wiley Online Library volume 36. J¨ unger, M., Liebling, T. M., Naddef, D., Nemhauser, G. L., Pulleyblank, W. R., Reinelt, G., Rinaldi, G., & Wolsey, L. A. (2009). 50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-art. Springer Science & Business Media. Kalogerakis, E., & S. Chaudhuri, V. K., D. Koller (2012). A probabilistic model for component-based shape synthesis. ACM transactions on graphics (TOG), 31 . Keha, A. B., de Farias Jr, I. R., & Nemhauser, G. L. (2004). Models for representing piecewise linear cost functions. Operations Research Letters, 32 , 44–48. Kim, J.-G., & Kim, Y.-D. (2003). A linear programming-based algorithm for floorplanning in VLSI design. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 22 , 584– 592. Kothari, R., & Ghosh, D. (2012). The single row facility layout problem: state of the art. Opsearch, 49 , 442–462. Leach, N. (2009). Digital cities. Architectural Design, 79 , 6–13. Li, J., & Smith, A. E. (2018). Block layout for attraction-based enterprises. European Journal of Operational Research, 266 , 1100–1112. Lynch, K. (1984). Good city form. MIT press. Lynch, K., Lynch, K. R., & Hack, G. (1984). Site planning. MIT press. Meeda, B., Parkyn, N., & Walton, D. S. (2007). Graphics for urban design. Thomas Telford. Meindl, B., & Templ, M. (2012). Analysis of commercial and free and open source solvers for linear optimization problems. Eurostat and Statistics Netherlands within the project ESSnet on common tools and harmonised methodology for SDC in the ESS , (p. 20). Merrell, P., Schkufza, E., & Koltun, V. (2010). Computer-generated residential building layouts. ACM Transactions on Graphics, 29 , 485–501. Moughtin, C. (2007). Urban design: street and square. Routledge. Pauly, M., Mitra, N. J., Wallner, J., Pottmann, H., & Guibas, L. J. (2008). Discovering structural regularity in 3D geometry. ACM transactions on graphics (TOG), 27 , 43. Peng, C.-H., Yang, Y.-L., Bao, F., Fink, D., Yan, D.-M., Wonka, P., & Mitra, N. J. (2016). Computational network design from functional specifications. ACM Transactions on Graphics (TOG), 35 , 131. Peng, C.-H., Yang, Y.-L., & Wonka, P. (2014). Computing layouts with deformable templates. ACM Transactions on Graphics (TOG), 33 , 99. Ploskas, N., & Samaras, N. (2015). Efficient GPU-based implementations of simplex type algorithms. Applied Mathematics and Computation, 250 , 552–570. Polzin, T., & Daneshmand, S. V. (2001). A comparison of steiner tree relaxations. Discrete Applied Mathematics, 112 , 241–261. ´ (2013). An evolutionary strategy enhanced with a local Rodrigues, E., Gaspar, A. R., & Gomes, A. search technique for the space allocation problem in architecture, part 1: Methodology. ComputerAided Design, 45 , 887–897. Rowley, A. (1994). Definitions of urban design: The nature and concerns of urban design. Planning Practice and Research, 9 , 179–197.

26

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Samaniego, H., & Moses, M. E. (2008). Cities as organisms: Allometric scaling of urban road networks. Journal of transport and land use, 1 , 21–39. Schenk, L. (2013). Designing cites: basic - principles - projects. Birkh¨ auser. Schrijver, A. (1998). Theory of linear and integer programming. John Wiley & Sons. Schumacher, P. (2009). Parametricism: A new global style for architecture and urban design. Architectural Design, 79 , 14–23. Shen, J., Shigeoka, K., Ino, F., & Hagihara, K. (2017). An out-of-core branch and bound method for solving the 0-1 knapsack problem on a GPU. In International Conference on Algorithms and Architectures for Parallel Processing (pp. 254–267). Springer. Singh, S. P., & Sharma, R. R. (2006). A review of different approaches to the facility layout problems. The International Journal of Advanced Manufacturing Technology, 30 , 425–433. Smelik, R. M., Tutenel, T., Bidarra, R., & Benes, B. (2014). A survey on procedural modelling for virtual worlds. Computer Graphics Forum, 33 , 31–50. Srinivasan, K., Chatha, K. S., & Konjevod, G. (2006). Linear-programming-based techniques for synthesis of network-on-chip architectures. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 14 , 407–420. ˇ St’ava, O., Beneˇs, B., Mˇ ech, R., Aliaga, D. G., & Kriˇstof, P. (2010). Inverse procedural modeling by automatic generation of L-systems. Computer Graphics Forum, 29 , 665–674. Talton, J. O., Lou, Y., Lesser, S., Duke, J., Mˇ ech, R., & Koltun, V. (2011). Metropolis procedural modeling. ACM Transactions on Graphics (TOG), 30 , 11. Trummer, P. (2008). Engineering ecologies. Architectural Design, 78 , 96–101. Vielma, J. P., Ahmed, S., & Nemhauser, G. (2010). Mixed-integer models for nonseparable piecewiselinear optimization: Unifying framework and extensions. Operations research, 58 , 303–315. Wolsey, L., & Rinaldi, G. (1993). Integer programming and combinatorial optimization. CIACO. Wu, W., Fan, L., Liu, L., & Wonka, P. (2018). MIQP-based layout design for building interiors. Computer Graphics Forum, 37 , 511–521. YANG, S., & CHEN, W. (2005). Study on reasonable density of residence [j]. City Planning Review , 3 , 014.

27