12 Reducing concrete consumption is important as part of the global effort of fighting the climate 13 change, and specifically in concrete flat slabs as these are among the largest cement consumers. In 14 this study we formulate an efficient gradient-based optimization of column locations, that minimizes 15 the slabs’ thickness with constraints on the deflections, bending moments and shear stresses while 16 accounting for architectural considerations. The results show that the columns’ optimal locations 17 are not trivial and that the slab thickness is very sensitive to the columns’ exact locations. Thus, 18 concrete savings in slabs of up to 20% are possible with minor modification to traditional layouts of 19 columns, and up to 50% with more pronounced updates, which emphasizes the importance of early 20 collaboration between architects and engineers. The results indicate the critical trade-off between 21 structural efficiency and architectural freedom and demonstrate the potential of formal optimization 22 in structural design. 23 Keywords: Concrete floors, Structural optimization, Columns layout, Structural Engineering, 24 Architectural constraints 25 1. Introduction 26 Concreteisoneofthemosthighlyconsumedmaterialsintheworld,beingthethirdlargestsourceof 27 carbondioxideemissions[1]. Consideringstructuralelementsinbuildings,alargeportionofconcrete 28 is used in slabs. In fact, several recent studies investigated the usage of cement in different structural 29 components in buildings and infrastructure, and it was shown that slabs hold the highest share of 30 cement [2, 3, 4]. Therefore, reducing the volume of concrete in slabs has high potential for reducing 31 Email address: yzelick1@jhu.edu (Yakov Zelickman) Preprint submitted to Architecture, Structures and Construction March 2, 2024 the environmental burden caused by cement production [4]. Moreover, slabs in buildings contribute 32 significantly to the mass of the structure, and consequently to the gravitational and earthquake loads 33 thatthebuildingmustwithstand. Thus,reducingtheslabs’masswillleadtofurtherconcretesavings 34 in other structural elements, such as columns and foundations. 35 Structural optimization is a design approach where a structural design problem is formulated as a 36 constrained minimization problem and solved with mathematical programming tools [5]. It has been 37 shownasaneffectivedesigntoolinmanybranchesofengineeringthatoftenleadstosignificantsavings 38 in material and improvements in performance [6, 7]. Thus, structural optimization is a promising 39 design approach to reduce the environmental impact of concrete structures [4, 8, 9]. 40 Optimizationofconcretefloorsystemswherethecolumnlocationsarefixed,isthesubjectofmany 41 studies,aimingtominimizeobjectivefunctionssuchasmaterialconsumption,cost,andenvironmental 42 footprint. To name a few, Varaee and Ahmadi-Nedushan [10] minimized the cost of uni-directional 43 flat plates with a single span, whereas cost minimization of flat plates with arbitrary shapes can be 44 foundin[11]. CostoptimizationofawaffleslabwaspresentedbyOlawaleetal.[12],whoformulated 45 a compact geometrical parametrization and therefore used a Genetic Algorithm (GA) for solving the 46 optimization. Richerparametrization,thatallowstheshapeoftheribstovarywasrecentlypresented 47 by Ismail and Mueller [13]. Some papers proposed optimization methods that consider multiple 48 options for the floor structural system, for example [14]. 49 Thelayoutofcolumns,andmoregenerallythelayoutofsupports,significantlyaffectsthestructural 50 responseofplates. Therefore,optimizingthelocationsofthesupportingelementsisaneffectivewayof 51 reducingtheenvironmentalfootprintofconcreteslabs[15,16]. Inanearlypaper,theauthorsminimize 52 the cost of a rectangular flat plate by optimizing a comprehensive set of parameters, including the 53 span lengths [17]. Therein, a two-step framework is presented where the floor is optimized using 54 GA for given span lengths, which are then updated following a heuristic scheme. Optimization of 55 a rectilinear floor was presented in Shaw et al. [18], where the authors used GA to optimize the 56 2 layoutofprefabricatedslabelementsandthesupportingcolumns. Inamorerecentstudy,theauthors 57 use Ant Colony Optimization to optimize the layout of an orthogonal-supported rectilinear building 58 [19]. Additionally, the floor plan is optimized with a constraint on the total floor area. The objective 59 function includes the cost of the frames and the slabs, and the eccentricity between the mass- and the 60 rigidity- centers. Recently, Building Information Modelling (BIM) was coupled with Finite Element 61 (FE) analysis and GA to create a framework for preliminary design of concrete structures, including 62 spacingbetweencolumngrid-lines[15]. Inanotherrecentstudy,theauthorsuseMonteCarlomethod 63 tofindtheoptimallocationsofsupportsofconcreteplates,minimizingthestrainenergy,reinforcement 64 steel and maximal deflection [20]. 65 Allstudiesthatwerementionedsofar,andmostoftheavailableliteraturethatdiscussesoptimiza- 66 tion of concrete floor systems, adopts meta-heuristic and zero-order optimization algorithms, which 67 allowtocopewiththenon-differentiableanddiscontinuousconstraints,butalsobecomesveryexpen- 68 sive computationally in high dimensional optimization [21]. Therefore, the design space includes a 69 small number of design variables, restricting the optimization to regular layouts of columns or to a 70 limited number of columns. 71 Gradient-based optimization algorithms are more likely to converge to local minima than meta- 72 heuristic algorithms, but offer superior numerical efficiency and therefore were also considered in 73 many studies. In a straightforward approach for optimization of supports’ locations, the coordinates 74 of the supported nodes are being optimized [22, 23]. This approach requires constant remeshing and 75 control over the FE mesh, and therefore is numerically expensive and may encounter stability issues. 76 Another approach, that uses a SIMP-like parametrization, was proposed by Buhl [24]. Mathematical 77 continuity is obtained by adding springs to all nodes and assigning penalized topological design 78 variablestoeachoneofthesprings. Thus,byaddingaconstraintonthesumofthetopologicaldesign 79 variables, the most effective springs remain, designating the optimal locations. This approach was 80 adopted in several studies, for example Jihong and Weihong [25], Denli and Sun [26], and recently 81 3 usedbyMengetal.[27]tominimizethecomplianceofplateroofstructures. Inarecentpaperbythe 82 authors we introduce the stiffness projection method for support location optimization, which is both 83 numerically efficient and mesh-independent, and therefore much less prone to convergence to local 84 minimathantheothergradient-basedapproaches[28]. SimilarlytoMengetal.[27],theoptimization 85 formulationthereispurelyacademic,whereunconstrainedcomplianceminimizationwasadoptedfor 86 the purpose of establishing the stiffness projection method. 87 From the discussion above it is apparent that existing studies on column layout optimization 88 of concrete floors were using meta-heuristic algorithms, mainly GA. As a result, the design space is 89 limitedtoasmallnumberofdesignvariables,andthereforetheexistingmethodsfocusonaregulargrid 90 of columns and simple floor plans. In most examples, architectural constraints were not considered, 91 or they were limited to restricting the span length of the regular grid of columns. On the other hand, 92 studiesthatusedefficientgradient-basedalgorithms,whichresultinarichdesignspace,consideronly 93 global structural performance and lack the necessary structural and architectural constraints for the 94 design of concrete floors. 95 In this study we aim to fill this gap by proposing a gradient-based optimization for the layout of 96 the columns in floors with arbitrary shapes, considering the major structural design requirements of 97 concrete plates as well as imposing general architectural-geometrical constraints. In this regard we 98 notethatthefocusofthepaperisonresidentialandofficebuildings,howeverothertypesofbuildings 99 andfloorsystemscanbeconsideredaswellwithappropriateloadconditions,andstructuralmodelsfor 100 theslab. Additionally,weinvestigatethesensitivityoftheplatethicknesstothelocationsofcolumns. 101 Specifically,webuildoffthestiffnessprojectionmethodthatwaspresentedintheauthors’previous 102 work [28], and extend it significantly by adding deflection, punching shear, bending moment, and 103 architectural constraints, as well as adding the plate thickness to the design space and explicitly 104 minimizing the concrete volume. As a result, the columns are not restricted to any location or 105 pattern(otherthanthearchitecturalconstraints),whichgivesrisetonon-triviallayoutsandsignificant 106 4 reduction in concrete mass. As expected, there is a clear trade-off between concrete volume and the 107 architectural design freedom. Less expected is how significant are the concrete savings when only 108 mildchangesaremadetothecolumnlocations,imposingonlyminorcompromiseonthearchitectural 109 design. This study demonstrates the importance of close collaboration between structural engineers 110 and architects from the preliminary design stages [18], when the column locations are determined. 111 The remainder of this paper is arranged as follows. In the next section we briefly present the 112 mathematical model, thereafter in Section 3 we discuss the optimization formulation extensively. In 113 Section 4 we present three numerical examples that are followed by a brief discussion and some 114 concluding remarks in Section 5. The paper has two appendices: The first presents the analytical 115 sensitivityanalysisandthesecondprovidessomedetailsabouttheimplementationoftheoptimization 116 method. 117 2. Mathematical model 118 Instructuraloptimization,themathematicalmodelisastructuralmodelthatpredictsthestructural 119 response to a given set of loads for a given set of parameters, including the design parameters. In the 120 contextofthecurrentstudy,thestructuralmodelisaplatemodelwherethesupportslocationsandthe 121 thicknessmayvarythroughouttheoptimization. Wenotethatthemathematicalmodelthatweusein 122 the current study was already presented and discussed extensively in our previous work [28]. Herein 123 the mathematical model is described briefly for completeness. 124 The slabs are modeled with plate finite elements using Mindlin plate theory [29, 30]. Following 125 commonpracticeintheanalysisofconcreteslabs,weassumesmalldisplacementsandstrainsaswell 126 aslinearelasticbehavioroftheconcrete. Thus,thefloorismodeledwith4-nodedplateelementswith 127 mixed interpolation, that are known to be accurate and insensitive to shear locking [31]. 128 Since we optimize the locations of the columns, the boundary conditions of the slab change 129 throughout the optimization. Generally, this class of problems suffers from several difficulties: 1) 130 Possible discontinuity of the design space; 2) High computational cost if remeshing is used; and 3) 131 5 Tendency to converge to poor local optima. Therefore, in this study we use the stiffness projection 132 method that was presented in our previous work to overcome these challenges [28]. 133 As the name suggests, the basic idea is to project the stiffness of the columns upon the plate’s FE 134 mesh instead of modeling the columns explicitly. Thus, all nodes within a circular projection area 135 Ω𝑖defined by a projection radius of 𝜂𝑖, have added nodal stiffness. This added stiffness equals to the 136 column’s stiffness matrix multiplied by a weight factor 𝑤𝑖𝑗that relates the 𝑖thcolumn with the 𝑗th137 node. Thereafter, the added nodal stiffness matrices are assembled into a global equivalent stiffness 138 matrix of the 𝑖thcolumn 139 K𝑐𝑝,𝑖=∑︁ Ω𝑖𝑤𝑖𝑗K𝑐,𝑖withΩ𝑖= 𝑗|𝑟𝑖𝑗≤𝜂 . (1) In the above expression, K𝑐𝑝,𝑖andK𝑐,𝑖are the𝑖thcolumn equivalent and nominal stiffness matrices; 140 𝑟𝑖𝑗is the distance between the 𝑖thcolumn and the 𝑗thnode; and the sum operator represents assembly 141 according the degrees of freedom of the model. Because we use gradient based optimization in this 142 study, all functions have to be differentiable and therefore we use a smooth radial super-Gaussian 143 function for the projection weight factors 144 ˜𝑤𝑖𝑗=˜𝑤𝑟𝑖𝑗=exp −0.5𝑟𝑖𝑗 𝜂2𝛽! , (2) where𝛽isaparameterthatcontrolsthesharpnessofthetransitionacrosstheboundaryoftheprojection 145 area. ThismeansthatmathematicallythestiffnessofanycolumnisprojectedontoallnodesoftheFE 146 mesh, with practically zero projection weight outside the desired projection area. To ensure that no 147 excess stiffness is generated by the projection, we normalize the projection weights 148 𝑤𝑖𝑗=˜𝑤𝑖𝑗Í 𝑘˜𝑤𝑖𝑘with𝑘=[1...𝑁𝑛], (3) where𝑁𝑛is the total number of nodes. After the equivalent stiffness matrices of all columns have 149 been computed, they are added to the plate’s stiffness matrix K𝑝, which results in the stiffness matrix 150 of the supported plate, K. Finally, walls are modeled as pinned supports that are added to all nodes 151 6 within the projection area of the walls on the FE mesh. To ensure that all walls are modeled, in all 152 examples we use regular FE meshes with element size that is less than the wall thickness. 153 3. Optimization problem formulation 154 In this section we present the proposed design-oriented problem formulation. Thus, we minimize theconcretevolume,considerthemajorserviceanddesignlimitstaterequirements,aswellasadding architectural constraints (AC). Arranging the optimization problem into standard form, we obtain minimize X𝑓0=𝑉 s.t.𝑓1=˜𝛿 ˜𝛿∗−1≤0 𝑓2=˜𝜎𝑡𝑠 ˜𝜎∗ 𝑡𝑠−1≤0 (4) 𝑓3=˜𝜇 ˜𝜇∗−1≤0 ˜X∈Ω𝑎𝑟𝑐 with: Kus=fs Kud=fd. In the formulation above: 𝑓1,𝑓2and𝑓3are respectively the deflection, shear stress, and bending moment constraints. The next set of constraints are the architectural constraints, which are imposed through the design space definition. The last two expressions represent the equilibrium equations that are considered in a nested configuration. Although the formulation is general and any number of load cases can be accommodated, in this study all examples have only two different load cases with uniformly distributed loads that correspond to service and design limit states. The service limit state and the design limit state load vectors, f𝑠andf𝑑, are given by   f𝑠=g+Δg+q f𝑑=1.4(g+Δg)+1.6q, 7 where g,Δgandqarethenodalselfweight,deadload,andliveloadvectors,respectively. Wenotethat 155 including pattern loading would represent the expected loads on the floor more accurately. However, 156 the optimized design will be influenced only marginally by a pattern loading, because different load 157 patters balance each other in the context of the optimal location of the columns. Furthermore, in 158 residentialbuildingsaswellasinmostofficebuildings,theliveloadsaremuchsmallerthanthedead 159 loads, diminishing the effect of pattern loading. Finally, as will be apparent from the results, the 160 optimizeddesignstendtohaveevencolumndistribution,andthereforetheeffectofpatternloadingis 161 evensmaller. Infact,undersuchcircumstancessomebuildingcodesdonotrequiretoconsiderpattern 162 loading [32, 33]. 163 The plate forces and moments are obtained in design limit state by computing 164 ˆS=DBu𝑑, (5) where ˆSis a vector with the plate forces and moments evaluated at the Gauss points, Dis the plate’s 165 constitutive matrix, and Bis a differentiation matrix. The nodal forces and moments are computed 166 using the SPR technique [34] 167 S= M𝑇 𝑥𝑥M𝑇 𝑦𝑦M𝑇 𝑥𝑦𝝈𝑇 𝑥𝑧𝝈𝑇 𝑦𝑧 𝑇=W𝑇ˆS, (6) where Sis a vector with the nodal bending moments and transverse shear forces, and Wis a constant 168 transformation matrix. Finally, Xis the normalized mathematical design vector, whereas ˜Xis the 169 physicaldesign vector that holds the actual design parameters that we wish to optimize. 170 3.1. Architectural design space 171 Duetotheexplicitparametrization,inwhichthecoordinatesofthecolumnlocationsareconsidered 172 asdesignvariables,architectural-geometricallimitationsmaybeintroducedbyappropriatedefinition 173 ofthedesignspace. Thisinherentlyensuresthatthearchitecturalrequirementsaremetwhileexploring 174 the design space, hence the title Architectural Design Space. However, for brevity we will refer to 175 the architectural designs space, simply as the design space. Another advantage of considering the 176 8 architecturalrequirementsthroughthedesignspaceisthattheydonothavetobedifferentiable,which 177 increases the range of architectural constrains that can be imposed. 178 Thus,thearchitecturalconstraints(AC)inEq.(4)areimposedthroughthedesignspacedefinition, 179 ˜X∈Ω𝑎𝑟𝑐, whereΩ𝑎𝑟𝑐is a set of all architecturally admissible designs. This set encodes the specific 180 architectural requirements for a design problem and defines the allowable range for each design 181 variable. Considering a floor that is supported on 𝑁𝑐𝑜𝑙columns, the design variables vector is 182 ˜X𝑇= x𝑇 𝑐,y𝑇 𝑐,ℎ . (7) where, x𝑐andy𝑐arevectorswiththe 𝑥and𝑦coordinatesofthecolumnsand ℎisthethicknessofthe 183 slab. Accordingly,thedesignspacehas 𝑁𝑑𝑣=2𝑁𝑐𝑜𝑙+1designvariables,andtheadmissibleset Ω𝑎𝑟𝑐184 has the same dimension. 185 The limits of the thickness design variable are straightforward: ℎ𝑚𝑖𝑛≤ℎ≤ℎ𝑚𝑎𝑥, whereℎ𝑚𝑖𝑛186 arisesfromstructuralbuildingcodesandregulations,and ℎ𝑚𝑎𝑥isanarchitecturalparameter,defining 187 the available space for the structural slab. However, for general floors with non-convex contours 188 and arbitrary geometrical-architectural constraints, the limits of the column locations variables are 189 design-dependent and therefore more complex. In this regard, we distinguish between the trivial AC 190 that ensure that the columns remain within the floor area, and the characteristic AC that represent all 191 theadditionalrequirements,suchasrestrictingacolumnmovementtoaspecificregion. Asingeneral 192 constrainedoptimization,notallconstraintsmustbeactiveandacharacteristicACusuallywillmake 193 the corresponding trivial AC inactive. 194 Considering the trivial AC, we require that at each design iteration, the updated location of a columnwillremaininthecircledefinedbythecurrentlocationofthecolumnandtheshortestdistance to the boundary, which includes both the contour of the floor and the openings. Therefore, for a columnwithshortestdistanceof 𝑑𝑚𝑖𝑛,thedesignlimitsineachdirectionareconservativelysetto𝑑𝑚𝑖𝑛 2, 9 which results in the following limits  x𝑐,𝑚𝑎𝑥,x𝑐,𝑚𝑖𝑛 =x𝑐±1 2d𝑚𝑖𝑛(x𝑐,y𝑐) (8)  y𝑐,𝑚𝑎𝑥,y𝑐,𝑚𝑖𝑛 =y𝑐±1 2d𝑚𝑖𝑛(x𝑐,y𝑐), (9) where d𝑚𝑖𝑛is a vector with the shortest distance from all columns to the boundary of floor. 195 To compute 𝑑𝑚𝑖𝑛we discretize the floor boundaries to sampling points (SP) with spacing of 196 roughly 0.1[𝑚]between adjacent SP. Thereafter, we compute the distance from the column to all SP, 197 select the two closest SP, and approximate 𝑑𝑚𝑖𝑛with shortest distance to the line connecting both SP. 198 Additionally,weusethederivativesoftheshortestdistancewithrespecttothecolumncoordinatesto 199 identifythedirectiontothenearestboundary. Forexample,consideringthe 𝑥coordinateofacolumn, 200 a positive derivative indicates that the closest SP is somewhere to the left of the column location. 201 Therefore,thedesignlimittotherightmaybelargerandisdefinedbythemaximalmovelimitvalue, 202 which is discussed in Appendix B.2. Similar logic applies also to a negative sign of the derivative 203 and when considering the 𝑦coordinate. We note that in a case of close vicinity of a column to an 204 ear vertex of the boundary polygon, the proposed strategy may allow the column to exit the domain. 205 However, since the columns naturally prefer to remain strictly within the floor area, setting small 206 enough distance between the SP resolves any related issues. 207 The characteristic AC represent additional limitations on the column locations. In engineering 208 practice it is often that the multidisciplinary design is performed sequentially. Thus, first an architect 209 designsthefunctionalityandthelayoutofthepartitioningwalls,andthenthelocationofthecolumns 210 aredefinedbythestructuralengineeralongthosewalls. Therefore,oneclassofcharacteristicACare 211 pathconstraints,whereforagivencolumnwedefineapathalongwhichitisallowedtomove. Froma 212 mathematicalperspective,eachcolumnisassignedwithascalarparameter 𝑡𝑖,andbothcoordinatesof 213 the column are defined by explicit path functions, which encode the allowed path: 𝑥𝑐,𝑖=𝑥𝑐,𝑖(𝑡𝑖)and 214 𝑦𝑐,𝑖=𝑦𝑐,𝑖(𝑡𝑖). In section 4.2 we will present example of path constraints. A natural extension of the 215 pathconstraintsarearealconstraints,wherethecolumnsareallowedtomoveonlywithinapredefined 216 10 2D subdomain. We will use such areal constraints in Section 4.3 to investigate the relation between 217 the design freedom and the concrete savings. 218 Asmentioned,wedistinguishbetweenthe physicaldesignvariables ,whichreferstotheactualpa- 219 rametersthatwewanttofind,andthe mathematicaldesignvariables ,thatwesolveintheoptimization 220 problem 221 X𝑇= r𝑇 𝑐,s𝑇 𝑐,𝜔 . (10) Themathematicaldesignvariablesarenormalizedandlinearlyrelatedtothephysicaldesignvariables 222 X=N˜X, (11) where Nis the diagonal normalization matrix. The entries on the diagonal of Nare1/𝐵𝑥,1/𝐵𝑦, or 223 1/(ℎ𝑚𝑎𝑥−ℎ𝑚𝑖𝑛)for the column locations in 𝑥and𝑦directions, and the slab thickness, respectively. 224 This normalization generally results in more stable optimization and conveniently separates the opti- 225 mizationprocedurefromthespecificgeometricalparametersoftheproblembeingsolved. Thelimits 226 on the mathematical design variables are obtained by normalization of the physical design limits. 227 3.2. Volume Objective function 228 As stated, we wish to minimize the concrete consumption, and therefore minimize the concrete 229 volume. We measure the concrete volume by summing the volumes of the individual finite elements 230 𝑉=𝑁ℓ∑︁ ℓ=1ℎ𝐴ℓ, (12) where𝐴ℓis the area of the ℓthfinite element, and 𝑁ℓis the total number of elements in the FE mesh. 231 3.3. Deflection Constraint 232 Many standards define the allowed deflection in concrete elements as a fraction of their span. In 233 general, floors have multiple spans, each possibly with different length, and therefore different areas 234 of a floor might have different allowed deflection. To successfully impose deflection constraints we 235 11 definetherelativedeflectionateachnode 𝛿𝑗astheratiobetweentheactualelasticdownwarddeflection 236 in service limit state and the allowed deflection at this node 237 𝛿𝑗=𝑤𝑗 𝑤𝐴,𝑗, (13) where𝑤𝑗and𝑤𝐴,𝑗are the actual and allowed deflections in 𝑧direction at the 𝑗thnode, respectively. 238 Theconstraintaggregatesallnodalrelativedeflectionsbyconsideringthemaximalrelativedeflection, 239 which is approximated using a 𝑝-norm function 240 ˜𝛿=©­ «𝑁𝑛∑︁ 𝑗=1𝛿𝑝 𝑗ª® ¬1 𝑝 . (14) In the equation above, ˜𝛿is the approximate maximal relative deflection, 𝑁𝑛is total number of nodes 241 intheFEmesh,and 𝑝isanevennumberallowingtoaccountforbothpositive(upward)andnegative 242 (downward) deflections. Moreover, since the deflections surface is smooth, we use high power value 243 of𝑝=30. This approximation overestimates the real maximum, ˜𝛿 > max(𝜹), which may lead to 244 undesired conservativeness. Therefore, the threshold value of the constraint is dynamically updated 245 as follows 246 ˜𝛿∗=˜𝛿 max(𝜹)𝛿∗, (15) every𝑁𝐼𝑐=5iterations, where the nominal required relative deflection is 𝛿∗=1.0. 247 The definition of the allowed deflection follows the recommendations in Eurocode 2 (EC2) [32], 248 wherethelongtermdeflectionshouldbelessthan1 250ofthespanlength. Thus,assumingalongterm 249 deflection coefficient of 3.0, the allowed deflection at node 𝑗is 250 𝑤𝐴,𝑗=𝐿𝑒𝑞,𝑗 750, (16) where𝐿𝑒𝑞,𝑗is the equivalent span length at this node. 251 Inatraditionalapproachforobtainingthespan lengthsinirregularfloors,thespansareidentified 252 bymanualinspectionofthedeflectionsurface,andthenthespanlengthsaresetasthediameterofthe 253 maximalinscribedcircleineachspan. However,sincethecolumnschangetheirlocationsthroughout 254 12 the optimization, it is not clear how to identify the irregular spans in a systematic and differentiable 255 manner. Therefore, we take a slightly different approach and directly approximate this diameter as 256 𝑑=√ 2𝑟𝑚𝑖𝑛,where𝑟𝑚𝑖𝑛isthedistancetotheclosestcolumnorwall. Thisapproximationisequivalent 257 to the traditional approach for regular and rectangular spans, and conservative otherwise. Thus, we 258 define the equivalent span length at any node 𝑗as follows, 259 𝐿𝑒𝑞,𝑗=𝑟0+√ 2𝑟𝑚𝑖𝑛,𝑗with𝑟𝑚𝑖𝑛,𝑗=min 𝑖𝑟𝑖𝑗, 𝑖∈[1,...,𝑁𝑐𝑜𝑙+𝑁𝑤𝑎𝑙𝑙],(17) where𝑟𝑚𝑖𝑛,𝑗is the distance from the 𝑗thnode to the closest column or wall and 𝑟0is a constant value 260 that we add to allow some minimal deflectionat the supports. This allowed deflection at the columns 261 is necessary to accommodate for the inevitable deflection at the supports, as the supports have finite 262 stiffness,asdiscussedinSection2. Wechosethevalue 𝑟0=0.7[𝑚],whichallowsanelasticdeflection 263 atthesupportsofapproximately 1×10−3[𝑚]. Finally,weapproximatethenon-differentiabledistance 264 to the closest column in Eq. (17) with a 𝑝-norm function 265 𝑟𝑚𝑖𝑛,𝑗≈ 𝑁𝑐𝑜𝑙+𝑁𝑤𝑎𝑙𝑙∑︁ 𝑖𝑟−𝑝 𝑖𝑗!−1 𝑝 . (18) 3.4. Shear Constraint 266 Shearinslabs,orpunchingshear,isakeyconsiderationinthedesignofconcreteslabsandhenceis 267 addedtotheformulation. Wedefineasufficientthicknessoftheslabsuchthatthepunchingresistance 268 at each point can be provided by steel details only, without further thickening. Thus, following the 269 recommendations in EC2, we will require for each node 𝑗that 270   𝜎𝑥𝑧,𝑗≤𝜈𝑅𝑑,𝑚𝑎𝑥 𝜎𝑦𝑧,𝑗≤𝜈𝑅𝑑,𝑚𝑎𝑥with𝜈𝑅𝑑,𝑚𝑎𝑥=0.4·0.6 1−𝑓𝑐𝑘 250 𝑓𝑐𝑑≈0.2𝑓𝑐𝑑. (19) In the expression above, 𝜎𝑥𝑧,𝑗and𝜎𝑦𝑧,𝑗are the plate transverse shear stresses acting at node 𝑗in 271 design limit state, 𝜈𝑅𝑑,𝑚𝑎𝑥is the maximal allowed shear stress, 𝑓𝑐𝑘and𝑓𝑐𝑑are the characteristic and 272 compression design strengths of the concrete (in [𝑀𝑃𝑎]). We note that we omit the eccentricity 273 13 parameter𝛽suggested by the EC2 [32], because the shear stresses are computed directly and thus 274 the actual structural response is already taken into account. Similarly to the deflection constraint, we 275 constrainthemaximalshearstressratherthanhavingseparatenodalconstraints. Thus,theapproximate 276 maximal shear stress is 277 ˜𝜎𝑡𝑠=©­ «2𝑁𝑛∑︁ 𝑗=1𝜎𝑝 𝑡𝑠,𝑗ª® ¬1 𝑝 with 𝝈𝑡𝑠=𝝈𝑥𝑧 𝝈𝑦𝑧. (20) Wenotethattheshearmaybebothpositiveandnegativeandthereforethevalueofthepower 𝑝should 278 be even. The threshold is updated in the same way as in the deflection constraint, 279 ˜𝜎∗ 𝑡𝑠=˜𝜎𝑡𝑠 max(𝝈𝑡𝑠)𝜈𝑅𝑑,𝑚𝑎𝑥. (21) For convenient presentation of the results, we define the relative shear stress as the ratio between the 280 nodal shear stress and the maximal allowed shear stress 281 𝝉𝑥𝑧=𝝈𝑥𝑧 𝜈𝑅𝑑,𝑚𝑎𝑥,𝝉𝑦𝑧=𝝈𝑦𝑧 𝜈𝑅𝑑,𝑚𝑎𝑥. (22) 3.5. Bending Moment Constraint 282 Another important design consideration in concrete elements is the bending moment capacity. In 283 slabs,itiscommonthatnocompressivesteelisneeded. Thus,inthisstudyweaimforstructuraldepth 284 that will subsequently allow a design with tensile steel only. Following recommendations in many 285 design codes including EC2 [32] and ACI[33], we assume a simplified rectangular stress block with 286 maximalheightof 0.4𝑑,where𝑑=ℎ−𝑑𝑠istheeffectivestructuraldepthand 𝑑𝑠istheconcretecover 287 ofthesteelrebars. Thus,themaximalbendingcapacityperunitwidth,withoutcompressionrebarsis 288 given by 289 𝑀𝑐=0.32(ℎ−𝑑𝑠)2𝑓𝑐𝑑, (23) which provides good approximation for 𝑓𝑐𝑑≤28[𝑀𝑝𝑎], especially as the moments approach 𝑀𝑐. 290 Followingcommonpractice,wetakeintoaccountthetorsionmomentsintheslabbyconsidering the Wood and Armer (W&A) moments [35]. Thus, we combine the pure bending moments with the 14 torsional moments to create the design moments 𝑀𝑟𝑥,𝑚𝑎𝑥=𝑀𝑥𝑥+ 𝑀𝑥𝑦 𝑀𝑟𝑥,𝑚𝑖𝑛=𝑀𝑥𝑥− 𝑀𝑥𝑦 𝑀𝑟𝑦,𝑚𝑎𝑥=𝑀𝑦𝑦+ 𝑀𝑥𝑦 𝑀𝑟𝑦,𝑚𝑖𝑛=𝑀𝑦𝑦− 𝑀𝑥𝑦 , where𝑀𝑥𝑥,𝑀𝑦𝑦,𝑀𝑥𝑦are the plate moments in design limit state. For convenient presentation of the 291 bending of the plate, we define the relative moment as the ratio between the nodal moments and the 292 moments capacities. Thus, the relative 𝑀𝑟𝑥,𝑚𝑎𝑥moment at any node 𝑗is 293 𝜇𝑟𝑥,𝑚𝑎𝑥,𝑗=𝑀𝑟𝑥,𝑚𝑎𝑥,𝑗 𝑀𝑐, (24) and similarly for the other moments. In order to constrain all moments at all nodes, we constrain the 294 approximate maximum relative moment 295 ˜𝜇=©­ «4𝑁𝑛∑︁ 𝑗=1𝜇𝑝 𝑗ª® ¬1 𝑝 with 𝝁=𝝁𝑟𝑥,𝑚𝑎𝑥 𝝁𝑟𝑥,𝑚𝑖𝑛 𝝁𝑟𝑦,𝑚𝑎𝑥 𝝁𝑟𝑦,𝑚𝑖𝑛. (25) Finally, the threshold value of the moment constraint is updated similarly to the shear and deflection 296 constraints, with normalized desired threshold value 𝜇∗=1 297 ˜𝜇∗=˜𝜇 max(𝝁). (26) 3.6. Optimization Sequence 298 It was observed during our numerical experiments that often only the displacement constraint is 299 active. Thus, in many cases the bending moment constraint and the shear constraint may be omitted. 300 This results in much faster optimization because it spares computing u𝑑as well as the corresponding 301 adjoint vectors, each requires solving a set of equilibrium equations which is the most expensive 302 15 computationaltask. Obviously,onecannotknowinadvancewhetherthedesignlimitstateconstraints 303 willbeactive. Therefore,inthisstudyweimplementedahierarchicaloptimizationsequence. Initially, 304 weoptimizeonlywiththedisplacementconstraintandcheckuponconvergencetheresultantmoment 305 and shear distribution. In a case that both the moment and shear values are within the desired limits, 306 theoptimizeddesignisconsideredasthesolutionoftheoptimizationproblem. Otherwise,weupdate 307 the optimized design by another optimization. This time, all constraints are included and the initial 308 design is the optimized design from the previous optimization. 309 3.7. Sensitivity analysis 310 In this study we use gradient-based optimization that allows to effectively cope with multidimen- 311 sionaloptimization,andspecificallyweadopttheMMAalgorithm[36]. Therefore,thederivativesof 312 allfunctionalsinEq.(4)withrespecttoalldesignvariablesshouldbederived,aprocessthatisoften 313 referred to as Sensitivity Analysis (SA). The SA of a functional 𝑓𝑖with respect to the mathematical 314 design variables is given by 315 𝜕𝑓𝑖 𝜕X=𝜕𝑓𝑖 𝜕˜XN−1(27) where𝜕𝑓𝑖 𝜕˜Xare the derivatives with respect to the physicaldesign variables, and presented in detail in 316 Appendix A. 317 4. Numerical examples 318 In this section we demonstrate the ability of the proposed method to reduce concrete volume 319 in slabs by optimizing the column locations. Additionally, the results presented here demonstrate 320 the critical trade-off between the structural efficiency and the architectural cost and emphasize the 321 importance of collaboration between architects and engineers at early stages of the project. The first 322 example is a simple design problem that validates the proposed optimization method, illustrates the 323 sensitivity of the slab thickness to the exact column location, and shows that the optimal location 324 of columns might be non-intuitive. The other two examples are inspired by real projects with more 325 16 complicatedgeometriesanddemonstratetheabilityoftheproposedmethodtocontributetoconcrete 326 savings in complex, real-life projects while satisfying the architectural requirements. Furthermore, 327 we take advantage of the more realistic geometries to investigate the effects of the design freedom on 328 the concrete savings. 329 Table 1.Material properties and other design parameters used in the current study symbol value units description 𝑑𝑠0.025 [m] concrete cover 𝑓𝑐𝑑17.40[𝑀𝑃𝑎]concrete compression design strength 𝐸30000[𝑀𝑃𝑎]concrete modulus of elasticity 𝜈0.3 [-] concrete Poisson’s ratio 𝛾 5/6[-] shear strain correction factor Δ𝑔4[𝑘𝑁/𝑚2] dead load 𝑞1.5[𝑘𝑁/𝑚2] live load 𝛾𝑐25[𝑘𝑁/𝑚3] concrete weight density In all designs presented in this section, the deflection constraint is active. In fact, we see that 330 the deflections often reach the allowed values in many regions, indicating better utilization of the 331 feasible space and hence better design. Therefore, we introduce the mean relative deflection ¯𝛿as 332 an additional measure of the optimality, providing another basis for comparison between different 333 designs, where higher ¯𝛿values indicate better designs. Furthermore, we use ¯𝛿also as an absolute 334 gauge of the optimality, where we first need to establish the maximal possible value that corresponds 335 to the best design. In this regard, we note that since the deflection surface is smooth, whereas the 336 alloweddeflectionsurfaceisanintersectionofconesandplanes, ¯𝛿=1.0cannotbereached. Toassess 337 theactualbestdesignintermsofthedeflectionwerecallthattheslopeofthedeflectedsurfaceshould 338 be zero above the optimal location of the supports, as has been shown in other studies dealing with 339 17 supports optimization, for example [37, 28]. Thus, we can estimate the maximal theoretical mean 340 relative deflection, ¯𝛿𝑚𝑎𝑥, by considering representative cases of fixed-fixed and cantilever beams, 341 which have known analytical deflection curves. Following the reasoning in Section 3.3, we set linear 342 alloweddeflectionsfunctionswithzerovalueatthesupports,integratetherelativedeflectionsalongthe 343 representative beams, and divide by their lengths, which yields ¯𝛿0.77and ¯𝛿0.64, respectively. 344 Therefore, we can assume that ¯𝛿𝑚𝑎𝑥∈[0.64,0.77]. 345 4.1. Example 1: Single Column Optimization 346 In the first example we wish to validate the optimization method, demonstrate that the optimal 347 location of columns might not be trivial, and illustrate the sensitivity of the slab thickness to the 348 column location. For this purpose we consider a minimalist problem and assume no architectural 349 constraints, which together allow us to explore the entire design space. Thus we optimize a single 350 column in an apartment located in a typical residential tower. The floor plan and all dimensions are 351 plotted in Figure 1a, where the columns have square cross section with 0.35[𝑚]side lengths. In 352 addition to the column being optimized, the boundary conditions of the floor include seven columns 353 along the contour of the floor, symmetry boundary conditions along the inner edges, and a portion of 354 aninternalcorewithwallthicknessof 0.25[𝑚]. Thefloorisdiscretizedwith 0.24×0.24[𝑚]elements 355 and is subjected to additional external loads as listen in Table 1. The design space includes the 356 coordinatesofthecolumn,whichcanbeanywherewithintheapartment,asmarkedbythehexagonal 357 pattern in Figure 1a, and the thickness of the slab ℎ, which may vary between ℎ𝑚𝑖𝑛=0.05[𝑚]and 358 ℎ𝑚𝑎𝑥=0.5[𝑚]. 359 Theoptimizedcolumnlocationis (𝑥𝑐,𝑦𝑐)≈(10.3,5.6)[𝑚]andthecorrespondingslabthicknessis 360 ℎ=0.228[𝑚],withconcretevolumeof 𝑉=49.9 𝑚3 . Thisoptimizeddesign,aswellastheresultant 361 deflections𝑤, the allowed deflections 𝑤𝐴, and the relative deflection 𝛿, are presented in Figures 1b, 362 1c and 1d, respectively. Inspecting these figures, it is evident that the deflection constraint is active, 363 with the maximal deflection at the resultant three spans reaching the maximal allowed deflection, 364 18 (a)Floor plan, all dimensions in [𝑐𝑚] (b)Deflections𝑤[𝑚𝑚] (c)Allowed deflections 𝑤𝐴[𝑚𝑚] (d)relative deflections 𝛿 Fig. 1.Problem setup for the single column optimization, optimized column location and the resultant deflections maps. The deflection constraint is active with deflections in all three spans reaching the allowed deflection values. hence utilizing effectively the feasible space with ¯𝛿=0.516. Figure 2a presents the distribution of 365 the minimal relative moment in 𝑦direction, where−0.833≤𝜇𝑟𝑦,𝑚𝑖𝑛≤0.278indicating that 𝜇𝑟𝑦,𝑚𝑖𝑛366 is within the allowed range. In fact, this is true for all normalized moments and shear stresses, 367 which have similar distributions and yields extremum values of 𝜇𝑚𝑖𝑛=0.8706,𝜇𝑚𝑎𝑥=0.3645, and 368 𝜏𝑡𝑠,𝑚𝑎𝑥=0.985. Thus, following the discussion in Section 3.6, re-optimization with all constraints 369 was not necessary. 370 Once the optimized column location has been found, we verify it by manually investigating the 371 entire design space. Thus, we generate an optimal surface by finding the minimum required slab 372 thickness for every column location. We discretize the design space of the column location to a grid 373 of0.3×0.3[𝑚]points and minimize the slab thickness for a column fixed at each one of those grid 374 19 (a)𝜇𝑟𝑦,𝑚𝑖𝑛 (b)Optimal thickness surface, ℎ(𝑥𝑐,𝑦𝑐)[𝑚] Fig. 2.(a) Normalized minimal moment in 𝑦direction, showing that 𝜇𝑟𝑦,𝑚𝑖𝑛is feasible. (b) Investigation of the entire design space of(𝑥𝑐,𝑦𝑐). Each pixel is 0.3×0.3[𝑚]and represents the required slab thickness ℎ[𝑚]for each column location, with minimal value of ℎ=0.2322[𝑚]. The optimized location of the column is at the region with the minimal requiredthickness,implyingsuccessfulconvergencetotheoptimum. Theblackcrossesarethesuggestedlocationsforthe column, as obtained by 26 practicing structural engineers, and show that the optimal column location is not trivial. points. Figure 2b presents the optimal surface, where the color of each pixel represents the required 375 slab thickness for a column located at the centre of the pixel. The filled gray square is the optimized 376 columnlocationasobtainedfromthestraightforwardoptimizationanditisevidentthatitislocatedat 377 the optimum. In fact, the minimal required thickness obtained by the optimization is slightly smaller 378 than the minimal value obtained by the design space exploration, probably due to the continuous 379 nature of the design variables [𝑥𝑐,𝑦𝑐]. 380 Theoptimallocationofthecolumnisnotintuitive,asamoretraditionalsolutionwouldbetolocate 381 the column closer to the center of the large span. To test this hypothesis we performed a poll among 382 26 practicing structural engineers and asked them where would they locate the column, ignoring any 383 architectural considerations. The black crosses in Figure 2b represent the answers received from the 384 participants of the poll, where the numbers indicate multiple answers for a certain location. It can be 385 seenthatmostparticipantslocatedthecolumnapproximatelyatthecenterofthelargebay,markedby 386 theintersectionofthedottedgridlines. Theseresultsconfirmthattheoptimalcolumnlocationmight 387 notbeobvious,eveninasimplefloorgeometryandwithoutconsideringanyarchitecturalconstraints. 388 20 Foreachcolumnlocationsuggestedbytheparticipantsofthepoll,wecomputedtherequiredslab 389 thickness, which vary between ℎ𝑟𝑒𝑓∈[0.2688,0.2987][𝑚]with mean value of ¯ℎ𝑟𝑒𝑓=0.2823[𝑚]. 390 Thus, we can say that the concrete savings vary between 0%in the unlikely case that the column 391 was originally located at the most efficient location, and 23.8%when considering ¯ℎ𝑟𝑒𝑓. Moreover, 392 inspectingtheoptimalthicknesssurfaceinFigure2b,itseemsthatithassteepgradientvaluesaround 393 the optimal point, such that moving the column by one pixel, or 0.3[𝑚](in the correct direction), 394 results in thickness reduction of almost 10[𝑚𝑚]. This implies a high sensitivity of the slab thickness 395 tothecolumnlocation,atleastinregionswithlargegradientsoftheoptimalthickness. Inthecurrent 396 example,thisregionreachesthetraditionallocationofthecolumnandincludesmostofthesuggested 397 locations in the poll. Thus, indicating that significant savings could be achieved even with small 398 changes to a traditional design. 399 4.2. Example 2: Rounded Triangular Floor 400 The second example that we present is inspired by a floor plan of an actual building that was 401 presented in [38] in the context of post-tensioning optimization. Therefore, this example shows the 402 ability of the proposed method to deal with real-life problems characterized with many columns as 403 well as non-convex shapes of floors. Furthermore, we add path constraints to ensure that the design 404 comply with an architectural design of this floor plan. 405 The floor has a triangular shape with rounded corners, has three rectangular openings, and sup- 406 ported on 19 square columns as well as a central concrete core. The thickness of the core walls is 407 0.35[𝑚]andthereforethefloorismodeledwithslightlysmallersquareFE,with 0.333[𝑚]sidelength. 408 All other parameters are the same as in the previous example. Figure 3a depicts the floor plan and 409 some measures, whereas all geometrical data can be found in the Supplementary Material file. 410 Thereferencelayoutofthecolumnsfollowsthegenerallayoutin[38],andshowninFigure3a. To 411 computetheslabthicknessforthereferencedesignwefixthecolumnlocationsandoptimizeonlythe 412 thickness, which equals to ℎ𝑟𝑒𝑓=0.331[𝑚], with resultant concrete volume of 𝑉𝑟𝑒𝑓=285.81 𝑚3 . 413 21 In Figure 3b we present the relative deflections, where the deflections reach the allowed values only 414 atthecantileverspanatthebottomoftheplan,whileinbroadareasofthefloorthedeflectionsarefar 415 form the allowed values, yielding ¯𝛿𝑟𝑒𝑓=0.331. Re-optimization is not necessary because the design 416 limit state requirements are met, as can be seen in Table 2, that summarise all the results. 417 After establishing the reference design, we add architectural path constraints to the optimization, 418 where columns are allowed to move only along the partitioning walls. The partitioning walls are 419 plotted in Figure 3a with thin lines, and are inspired by the layout of the partitioning walls in an 420 architectural design of this floor, as can be found in San Francisco Planning Department’s website 421 [39], whereas the corresponding allowed paths are plotted in Figures 3b and 3c with white lines. 422 The optimized design is presented in Figure 3c, where the slab thickness is ℎ𝑎𝑟𝑐=0.271[𝑚]and the 423 concretevolumeis 𝑉𝑎𝑟𝑐=234.17 𝑚3 ,whichrepresentssavingsof 18.1%,relativelytothereference 424 design. Thus, significant savings were possible while satisfying all architectural constraints. The 425 second row of Table 2 summarises the results of the optimization with the path constraints. 426 Figure3calsopresentstherelativedeflections,wheremoreregionsofthefloorreachtheirallowed 427 deflection,andconsequentlytheaveragerelativedeflectionincreasesto ¯𝛿𝑎𝑟𝑐=0.427. Keepinginmind 428 thatthebottomcantileverwasthegoverningspanatthereferencedesign,weseethattheoptimization 429 increases the length of the adjacent span, which reduced the deflections at the cantilever and enabled 430 some thickness reduction. Another notable and somewhat less expected change in the design is the 431 movementofthecolumnattheverticalcenter-lineofthefloor(topoftheplan)totheedgeofthefloor, 432 increasingthelengthoftheendspan. Thisreducedthedeflectionsattheneighboringspansfromboth 433 sides and allowed for further reduction in the slab thickness. 434 So far we showed that significant savings are possible even when the design space is restricted to 435 theportioningwalls,asitwouldbeinasequentialarchitectural-engineeringdesignprocess. However, 436 early collaboration between architects and engineers may increase the design freedom and facilitate 437 further savings. To demonstrate the potential of a concurrent design process, we performed another 438 22 (a)Floor plan, Dimensions are in [𝑐𝑚] (b)Reference design with ℎ𝑟𝑒𝑓=0.331[𝑚]. (c)Optimized design with AC, ℎ𝑎𝑟𝑐=0.271[𝑚](−18.1%) (d)Free optimized design. ℎ=0.1836[𝑚](−44.5%) Fig. 3.Reference and optimization results for the rounded triangular floor. At the reference design the deflections reach the allowed value only in one region of the floor, indicating sub-optimal design with ¯𝛿𝑟𝑒𝑓=0.331. When optimized with architectural constraints, to ensure that the columns remain along the partitioning walls (marked with white lines), the deflection constraint becomes active in more spans with ¯𝛿𝑎𝑟𝑐=0.427. At the free optimized design the deflections reachtheallowedvalueinmostspanofthefloor,thusindicatinggoodutilizationofthefeasiblespacewith ¯𝛿𝑎𝑟𝑐=0.659. Surprisingly, one column merged with the core walls, thus practically eliminating this column. optimization, this time without any limitation on the column layout. The optimized slab thickness 439 isℎ=0.1836[𝑚]and the resultant concrete volume is 𝑉=158.607 𝑚3 , which represents volume 440 saving of 44.5%. Figure 3d presents the relative deflection map and the optimized column layout, 441 23 ℎ[𝑚]𝑉 𝑚3 𝛿𝑚𝑎𝑥¯𝛿 𝜇𝑟𝑥,𝑚𝑎𝑥𝜇𝑟𝑥,𝑚𝑖𝑛𝜇𝑟𝑦,𝑚𝑎𝑥𝜇𝑟𝑦,𝑚𝑖𝑛𝜏𝑥𝑧,𝑚𝑎𝑥𝜏𝑦𝑧,𝑚𝑎𝑥 ref 0.331 285.81 1.0 0.331 0.252 -0.451 0.220 -0.320 0.492 0.528 opt. arc. 0.271234.17 (-18.1%) 1.0 0.427 0.296 -0.574 0.260 -0.557 0.570 0.538 opt. free 0.184158.61 (-44.5%) 1.0 0.659 0.415 -0.822 0.436 –0.793 0.681 0.672 Table 2. Results of the rounded triangular floor. top row: reference design, middle row: optimized design with architectural constraints, bottom row: free optimization. It can be see that in all cases the deflection constraint is the only active constraint, where the rate at which the feasible space is utilized increases with the design freedom. whichdiffersnotablyfromthereferencelayoutinFigure3b. Surprisingly,oneofthecolumnsthatwas 442 originallylocatedinsidethecorehasmergedwiththecorewall. Thus,thecolumnisnotactiveandthe 443 optimization effectively converged to a solution with fewer columns. Looking at Figure 3d it is clear 444 thatthedeflectionsacrossmostofthefloorapproachthealloweddeflection,andaccordingly ¯𝛿=0.659. 445 All other constraints are inactive, as can be seen at the third row of Table 2, and re-optimization was 446 not needed. 447 In this example we included path constraints in the optimization and showed that even when the 448 column are allowed to move only along the partitioning walls, thus having no ’architectural cost’, 449 significant savings are possible. Furthermore, we showed that when the optimization is granted 450 with complete freedom, the concrete savings increases dramatically, but also the architectural cost 451 increases. In the next example we explore further this trade off between the structural efficiency and 452 the architectural cost. 453 4.3. Example 3: Irregular Residential Floor 454 In this section, we investigate the relation between the design freedom and the potential concrete savings. We approach this issue by examining the effect of the maximal allowed modification to the locations of columns, compared to a reference configuration – namely, a given architectural plan. 24 Fig. 4.Floor plan for example 3, including the reference layout of columns, dimensions are in [𝑐𝑚]. Thus, we introduce an areal architectural constraint through an additional layer of parametrization: 𝑥𝑖(𝑡𝑖,𝑠𝑖)=Δ𝑚𝑎𝑥𝑡𝑖+𝑥𝑖0,and𝑦𝑖(𝑡𝑖,𝑠𝑖)=Δ𝑚𝑎𝑥𝑠𝑖+𝑦𝑖0,with−1≤𝑡𝑖,𝑠𝑖≤1, whereΔ𝑚𝑎𝑥is the maximal allowed change in a column location, and (𝑥𝑖0,𝑦𝑖0)are the column’s 455 coordinates in the reference design. Thus, each column is allowed to move only within a local box 456 that is centered at the reference location of this column and has side lengths of 2Δ𝑚𝑎𝑥. Other shapes 457 can be as easily parameterized to account for a variety of areal architectural constraints. 458 The selected floor plan in this example is inspired by a floor geometry that was presented by He 459 etal.[40]inthecontextofyieldlineidentification. Thisisanirregularfloorinaresidentialbuilding, 460 supported by 19 square columns with side length of 0.35[𝑚]and several walls, as can be seen in the 461 floor plan in Figure 4. The thickness of the walls is 0.25[𝑚], and accordingly the FE mesh consists 462 of square elements with 0.2[𝑚]side length. The applied forces and other parameters are the same as 463 in the previous examples. All geometrical data, including the column locations, can be found in the 464 Supplementary Material file. 465 In the following we will perform series of optimizations with increasing values of Δ𝑚𝑎𝑥, repre- 466 sentingincreasinglevelsofarchitecturalfreedom,andtheninvestigatethetrade-offswiththeconcrete 467 savings. For convenience, all numerical results discussed herein are summarized at the end of this 468 25 section in Table 3. 469 We begin with Δ𝑚𝑎𝑥=0, which is essentially optimizing the thickness of the slab while keeping 470 thecolumnsattheirreferencelocation. Theresultingreferencethicknessis ℎ𝑟𝑒𝑓=0.2226[𝑚]andthe 471 corresponding concrete volume is 𝑉𝑟𝑒𝑓=82.15 𝑚3 . The first row in Table 3 summarise the results 472 ofthereferencedesign,whereitisevidentthatthedeflectionsconstraintistheonlyactiveconstraint, 473 and re-optimization was not necessary. 474 Next, we set Δ𝑚𝑎𝑥=0.1[𝑚], allowing for very minor adjustment of the column locations, which 475 probably has very little architectural cost. The optimized slab thickness is ℎ=0.2109[𝑚]which 476 reflects a reduction of 5.3%in concrete volume with respect to the reference design. Again, the 477 deflection constraint is the only active constraint and re-optimization is not necessary. 478 IncreasingΔ𝑚𝑎𝑥further leads to greater concrete savings, as can be seen in Figure 5 that depicts 479 the concrete volumes for different values of Δ𝑚𝑎𝑥. The right most point on the curve corresponds to 480 Δ𝑚𝑎𝑥=∞andrepresentstheextremecasewherethecolumnsarefreetomove,wheremathematically 481 weomittedthearchitecturalconstraintsforthisoptimization. Theoptimizedslabthicknessinthiscase 482 isℎ∞=0.1126, with concrete volume of 𝑉𝑜𝑝𝑡=41.531 𝑚3 , and not insignificant concrete savings 483 of 49.4%. 484 Increasingthedesignfreedomincreasestheextentatwhichtheoptimizationutilizestheavailable 485 feasible space. This is evident from Figure 5 that displays also the relative deflections for Δ𝑚𝑎𝑥= 486 {0,0.5,3,∞}[𝑚]where ¯𝛿increaseswiththedesignfreedom. Likewise,beginningat Δ𝑚𝑎𝑥=1.1[𝑚] 487 andon,themomentconstraintbecomesalsoactive. Thus,themomentandshearvaluesinitiallyexceed 488 thedesiredthresholdvalues,andre-optimizationwasnecessary,afterwhichallconstraintsaremetand 489 the objective function value increases (i.e., get worse). However, for Δ𝑚𝑎𝑥=3.0[𝑚]andΔ𝑚𝑎𝑥=∞ 490 the re-optimized design results in slightly thinner slab than the corresponding optimized thicknesses 491 whenonlydeflectionconstraintwasconsidered. Thus,indicatingconvergencetoalocalminimaatthe 492 first optimization attempts, which is not unlikely in non-convex optimization. Nevertheless, in most 493 26 Fig. 5.Optimized concrete volume for different values of Δ𝑚𝑎𝑥. The steep slope at low values of Δ𝑚𝑎𝑥indicates that even a small update of the column layout can significantly affect the concrete volume. The red asterisks represent infeasibleoptimizationresultsthatwereobtainedwiththedeflectionconstraintonly. Increasingthedesignfreedomresults in larger concrete savings that reach 49.4%. The color maps present the distribution of the relative deflection, 𝛿, for Δ𝑚𝑎𝑥={0,0.5,3.0,∞}[𝑚]. It is evident that increasing Δ𝑚𝑎𝑥results in more efficient design with higher ¯𝛿. cases adding more constraints leads to higher (worse) objective function values, as apparent from the 494 results in Table 3 and in Figure 5. Therein, the optimization trials that were re-optimized are marked 495 with an asterisk. 496 AnotherinterestingobservationfromthetherelativedeflectionplotsinFigure5isthatasthedesign 497 freedom increases, the columns distribution tend be more uniform, with small differences between 498 bay lengths. The reason for this is that large differences in adjacent bay lengths result in non-zero 499 slope of the deflection surface, and therefore are generally not optimal [37, 28]. Thus, we expect 500 optimized column layouts to be characterised with relatively uniform distribution, which can be used 501 tosetagoodinitialdesign. Interestingly,sincetheeffectofpatternloadingreduceswiththedifference 502 between bay lengths, including pattern loading in the formulation could result in larger savings in 503 27 concrete volume. Therefore, the obtained savings are possibly somewhat on the conservative side. 504 Additionally,wenotethattheoptimizationswith Δ𝑚𝑎𝑥=4.0[𝑚]andΔ𝑚𝑎𝑥=5.0[𝑚]convergedto 505 the same optimum. A possible explanation for this is the non-convexity of the optimization problem. 506 Thus,theoptimalsolutionmighthaveadiscretedependenceonthedesignspacefreedom. Thiscould 507 also explain why the optimal concrete volume that corresponds to Δ𝑚𝑎𝑥=∞is lower than one would 508 expect based on the graph in Figure 5. 509 Finally,since Δ𝑚𝑎𝑥canberegardedasameasureofthearchitecturalcost,thecurveinFigure5can 510 be interpreted as the trade-off between the architectural cost and the concrete volume. This trade-off 511 curve is convex and therefore small increase in the architectural cost with respect to a traditionally 512 obtainedreferencedesign,mayleadtosignificantreductioninconcretevolume. Forexample,allowing 513 Δ𝑚𝑎𝑥=0.9[𝑚]results in almost 30% reduction. 514 5. Discussion and conclusions 515 We presented a method to minimize the concrete consumption in slabs by optimizing the column 516 locations, and then use it to investigate the sensitivity of the thickness to the column locations. We 517 formulateadeflectionconstraintthatisdifferentiablewithrespecttothechangingcolumnslocations, 518 as well as moment and shear constraints, and include architectural constraints through the explicit 519 design parametrization. For any given floor plan, the method generates an optimized layout of 520 columnsandthecorrespondingminimalrequiredslabthickness. Weusegradient-basedoptimization 521 with analytically derived sensitivities, which results in a very effective numerical method that can 522 be used for problems with a large number of design variables, that would have not been practical 523 with zero-order optimization methods. For example, simultaneous optimization of a large number 524 of columns within an extended framework that includes also the column dimensions and a slab with 525 varying thickness. 526 Through three different examples we show the ability of the proposed method to cope with any 527 floor geometry, any number of columns and a range of architectural-geometrical constraints. The 528 28 Table 3.Optimization of the irregular slab with increasing level of design freedom.∗Indicates infeasible result obtained with deflection constraint only. Δ𝑚𝑎𝑥ℎ[𝑚]𝑉 𝑚3 𝛿𝑚𝑎𝑥¯𝛿 𝜇𝑟𝑥,𝑚𝑎𝑥𝜇𝑟𝑥,𝑚𝑖𝑛𝜇𝑟𝑦,𝑚𝑎𝑥𝜇𝑟𝑦,𝑚𝑖𝑛𝜏𝑥𝑧,𝑚𝑎𝑥𝜏𝑦𝑧,𝑚𝑎𝑥 0(ref) 0.223 82.146 - 1.000 0.376 0.380 -0.680 0.225 -0.703 0.530 0.753 0.1 0.211 77.829 5.3% 1.000 0.397 0.388 -0.694 0.238 -0.757 0.534 0.721 0.3 0.189 69.856 15.0% 1.000 0.453 0.408 -0.784 0.279 -0.853 0.524 0.680 0.5 0.175 64.546 21.4% 1.000 0.486 0.394 -0.834 0.338 -0.892 0.494 0.694 0.7 0.165 60.958 25.8% 1.000 0.516 0.382 -0.838 0.374 -0.919 0.473 0.805 0.9 0.156 57.693 29.8% 1.000 0.511 0.400 -0.863 0.391 -0.949 0.536 0.946 1.1* 0.148 54.765 33.3% 1.000 0.530 0.412 -0.903 0.402 -1.070 0.675 1.092 1.1 0.149 54.808 33.3% 1.000 0.524 0.412 -0.901 0.399 -0.999 0.679 0.997 1.4* 0.138 51.007 37.9% 1.000 0.548 0.468 -1.046 0.424 -1.170 0.874 1.116 1.4 0.142 52.219 36.4% 0.998 0.542 0.556 -0.937 0.420 -1.000 0.637 0.850 2.0* 0.127 46.976 42.8% 0.999 0.575 0.497 -1.215 0.466 -1.145 1.134 0.931 2.0 0.129 47.524 42.1% 0.998 0.552 0.502 -0.952 0.463 -1.000 0.776 0.781 3.0* 0.123 45.269 44.9% 1.000 0.601 0.516 -1.019 0.547 -1.095 0.556 0.794 3.0 0.121 44.820 45.4% 1.001 0.610 0.592 -0.949 0.544 -1.000 0.554 0.826 4.0* 0.119 43.745 46.7% 0.996 0.622 0.591 -1.040 0.510 -1.114 0.582 0.744 4.0 0.120 44.174 46.2% 0.996 0.596 0.608 -0.926 0.481 -0.999 0.572 0.676 5.0* 0.119 43.745 46.7% 0.996 0.622 0.591 -1.040 0.510 -1.114 0.582 0.744 5.0 0.120 44.174 46.2% 0.996 0.596 0.608 -0.926 0.481 -0.999 0.572 0.676 6.0* 0.119 43.745 46.7% 0.996 0.622 0.591 -1.040 0.510 -1.114 0.582 0.744 6.0 0.120 44.174 46.2% 0.996 0.596 0.608 -0.926 0.481 -0.999 0.572 0.676 ∞* 0.113 41.793 49.1% 1.000 0.638 0.587 -1.149 0.528 -1.342 0.694 1.029 ∞0.113 41.531 49.4% 1.000 0.613 0.521 -0.947 0.501 -1.001 0.576 0.774 29 results show that even when the columns are allowed to move only along the partitioning walls, 529 notablesavingsofupto 18%arepossible. Furthermore,weshowthatincreasingthegrantedfreedom 530 totheoptimizationdramaticallyincreasesthepotentialforconcretesavings,whichmayreachashigh 531 as50%. Interestingly,evenminorchangesincolumnlocations(intheorderof 0.1[𝑚]−0.9[𝑚])with 532 respect to a traditional design, may result in substantial savings of 5%−30%. 533 The results of this study indicate that the optimal column layout is not trivial and that traditional 534 designs are often sub-optimal. Thus, given the high sensitivity of the slab thickness to the columns 535 locations,theproposedmethodcanbeusedtofine-tunecolumnlayoutstoreducetheslabthicknessin 536 many buildings. Furthermore, collaborative architectural and structural design from the preliminary 537 stages when the architectural layout and consequently the column layout are determined, is key to 538 achieve significant concrete savings. 539 Another interesting observation that can be made is the direct relation between the optimality of 540 the column layout and the rate at which the deflection constraint is satisfied. Thus, the mean relative 541 deflection can be used as an indicator for the effectiveness of a design, with an estimated theoretical 542 maximum value in the range ¯𝛿𝑚𝑎𝑥∈[0.64,0.77]. 543 Throughout this research the architectural requirements were imposed as hard constraints based 544 onexistingarchitecturaldesigns. However,explicitandquantitativeconsiderationofthearchitectural 545 cost is an interesting direction for future research, which will allow to construct more meaningful 546 Paretofrontsthatillustratethetrade-offsbetweenstructuralefficiencyandarchitecturalperformance. 547 Additionally,althoughasubstantialreductioninconcretevolumeispossiblebytheproposedmethod, 548 itispossiblethatmoresteelreinforcementwillbeneeded. Thus,theoptimalbalancebetweenconcrete 549 and steel in terms of cost and environmental impact still remains open for future work. Furthermore, 550 including also non-linear material response and plasticity might reveal interesting failure modes. 551 30 6. Data Availability Statement 552 Some or all data, models, or code that support the findings of this study are available from the 553 corresponding author upon reasonable request. 554 7. Acknowledgments 555 This research was funded by the Ministry of Construction and Housing, Israel. The financial 556 support is gratefully acknowledged. 557 8. Conflict of Interest 558 The authors declare that they have no conflict of interest. 559 Appendix A. Sensitivity analysis 560 Since we implement gradient based optimization, the first-order derivatives should be provided. 561 In this section we present in detail all calculations involved in the computation of these derivatives. 562 Wenotethattheanalyticalsensitivitieswereverifiedbycomparingtonumericalderivativesobtained 563 with finite differences method and were found to be accurate. 564 The derivatives with respect to the mathematical design variables are obtained by the chain rule 565 𝜕𝑓𝛼 𝜕X=𝜕𝑓𝛼 𝜕˜X𝜕˜X 𝜕X,with𝛼∈{0,1,2,3}, (A.1) where𝜕˜X 𝜕X=N−1is the Jacobian matrix, and𝜕𝑓𝛼 𝜕˜xare the derivatives of the 𝛼functional with respect 566 the physical design variables and discussed in following sub-sections. 567 Appendix A.1. Volume Objective function 568 The sensitivities of the volume can be obtained explicitly because it does not depend on the 569 structural response. Thus, we differentiate Eq. (12) 570 𝜕𝑉 𝜕˜𝑋=𝑁ℓ∑︁ ℓ=1𝜕ℎℓ 𝜕˜𝑋𝐴ℓ, (A.2) 31 where ˜𝑋is any of the physical design variables. The derivative of the elemental thickness ℎℓwith 571 respect the slab thickness is simply𝜕ℎℓ 𝜕ℎ=1.0and zero with respect the columns locations. 572 Appendix A.2. Deflection Constraint 573 IntheperspectiveoftheindividualMMAiteration,thethresholdvalueoftheconstraintisconstant. 574 Therefore the derivative of the deflection constraint equals to the derivative of the maximal relative 575 deflection, scaled by 1/ˆ˜𝛿∗576 𝜕𝑓1 𝜕˜𝑋=1 ˜𝛿∗𝜕˜𝛿 𝜕˜𝑋. (A.3) Thus, we focus on the derivative of ˜𝛿. 577 Thedeflectionconstraintisanimplicitfunctionofthedesignvariablesandthereforeweadoptthe 578 adjointapproach. Thebasicideaistoaugmentthefunctionalwiththeequilibriumresidualmultiplied 579 byanadjointvectorthatwillbeselectedsuchthattheimplicittermswillvanish. Thus,theaugmented 580 functional is 581 ˜𝛿𝑎=˜𝛿−𝝀𝑇 𝛿(Ku𝑠−f𝑠). (A.4) Sincetheequilibriumresidualequalstozero,theaugmentedfunctionalequalstotheoriginalfunctional 582 and so are the derivatives. 583 Thus, we differentiate the augmented constraint with respect to the design variables. Keeping in 584 mindthatthedeflectionconstraintalsodependsexplicitlyonthedesignvariablesthroughtheallowed 585 deflection, we get 586 𝜕˜𝛿𝑎 𝜕˜𝑋=𝜕˜𝛿 𝜕u𝑠𝜕u𝑠 𝜕˜𝑋+𝜕˜𝛿 𝜕w𝐴𝜕w𝐴 𝜕˜𝑋−𝝀𝑇 𝛿𝜕K 𝜕˜𝑋u𝑠+K𝜕u𝑠 𝜕˜𝑋−𝜕f𝑠 𝜕˜𝑋 . (A.5) Since the derivatives of the augmented and original functionals are the same, we switch back to the 587 original functional. As mentioned, the adjoint vector is computed such that the terms𝜕u 𝜕˜𝑋will cancel 588 each other. Thus, the derivative of the deflection is 589 𝜕˜𝛿 𝜕˜𝑋=𝜕˜𝛿 𝜕w𝐴𝜕w𝐴 𝜕˜𝑋−𝝀𝑇 𝛿𝜕K 𝜕˜𝑋u𝑠−𝜕f𝑠 𝜕˜𝑋 with K𝑇𝝀𝛿=𝜕˜𝛿 𝜕u𝑠𝑇 (A.6) 32 The adjoint vector 𝝀𝛿and𝜕˜𝛿 𝜕ˆu𝑠are the same as presented in [41]. The derivative of the maximal 590 approximated relative deflection with respect to the allowed deflection is obtained by substituting 591 Eq. (13) into Eq. (14) and differentiating 592 𝜕˜𝛿 𝜕w𝐴=−˜𝛿"∑︁ 𝑗𝛿𝑝 𝑗#−1∑︁ 𝑗𝛿𝑝−1 𝑗𝑤𝑗(𝑤𝐴)−2 𝑗. (A.7) Thederivativeofthealloweddeflectionswithrespecttothedesignvariablesisobtainedbyreplacing 593 𝑟𝑚𝑖𝑛,𝑗in Eq. (17) with its derivative, and multiplying by1 750according to Eq. (16). Thus, by 594 differentiating Eq. (18) we obtain the derivative of the distance form the 𝑗thnode to the closest 595 column, 596 𝜕𝑟𝑚𝑖𝑛,𝑗 𝜕˜𝑋=𝑟𝑚𝑖𝑛,𝑗∑︁ 𝑖𝑟𝑝−1 𝑖𝑗𝜕𝑟𝑖𝑗 𝜕˜𝑋. (A.8) Thederivative𝜕𝑟𝑖𝑗 𝜕˜𝑋iscomputedbydifferentiatingthedistancebetweenthe 𝑖thcolumnandthe 𝑗thnode, 597 where for all design variables other than the 𝑖thcolumn location, the derivative is equal to zero. 598 The next term in Eq. (A.6) is the derivative of the stiffness matrix with respect to the design 599 variables, which were discussed in [28] and are brought here for completeness. 600 As mentioned, the stiffness matrix of the supported plate is simply summation of the plate’s 601 stiffness matrix and the equivalent matrices of the columns 602 K=K𝑝+𝑁𝑐𝑜𝑙∑︁ 𝑖=1K𝑐𝑝,𝑖. (A.9) Thus, the derivatives with respect the column locations affect only the added equivalent column 603 stiffnessmatrices. Thus,bydifferentiatingEq.A.9withrespectthe 𝑥coordinateofthe 𝑖thcolumnwe 604 get 605 𝜕K 𝜕𝑥𝑐,𝑖=𝜕K𝑐𝑝,𝑖 𝜕𝑥𝑐,𝑖. (A.10) The derivative of equivalent stiffness matrix of the 𝑖thcolumn with respect 𝑥𝑐,𝑖is obtained by differ- 606 entiating Eq. (1) 607 𝜕K𝑐𝑝,𝑖 𝜕𝑥𝑐,𝑖=𝑁𝑛∑︁ 𝑗𝜕K𝑐𝑝,𝑖 𝜕𝑥𝑐,𝑖 𝑗=𝑁𝑛∑︁ 𝑗𝜕𝑤𝑖𝑗 𝜕𝑥𝑐,𝑖K𝑐,𝑖. (A.11) 33 The summation sign stands for assembly according the nodal DOF. The derivative of the projection 608 weight is obtained by differentiating Eq. (2) and substituting into Eq. (3), 609 𝜕𝑤𝑖𝑗 𝜕𝑥𝑐,𝑖=𝜕˜𝑤𝑖𝑗 𝜕𝑥𝑐,𝑖Í 𝑘˜𝑤𝑖𝑘−˜𝑤𝑖𝑗Í 𝑘𝜕˜𝑤𝑖𝑘 𝜕𝑥𝑐,𝑖 (Í 𝑘˜𝑤𝑖𝑘)2(A.12) with 610 𝜕˜𝑤𝑖𝑗 𝜕𝑥𝑐,𝑖=−𝛽 𝜂𝑟𝑖𝑗 𝜂2𝛽−1𝜕𝑟𝑖𝑗 𝜕𝑥𝑐,𝑖˜𝑤𝑖𝑗. (A.13) The derivatives with respect to 𝑦𝑐,𝑖are computed in the same way. 611 The derivative of the stiffness matrix with respect the thickness design variable affect the plate’s 612 stiffnessmatrix, K𝑝,andareobtainedbydifferentiatingtheelementalstiffnessmatricesandthereafter 613 assembling in a regular manner. The plate’s stiffness matrix is assembled in a standard manner, for a 614 mesh with identical elements 615 K𝑝=∑︁ ℓKℓ=∑︁ ℓB𝑇 ℓDℓBℓ=∑︁ ℓB𝑇DB, (A.14) where BandDare the elemental generalized differentiation and constitutive matrices. Thus, after 616 differentiating we get 617 𝜕K𝑝 𝜕ℎ=∑︁ ℓB𝑇𝜕D 𝜕ℎB, (A.15) where the derivative of the constitutive matrix is computed by explicit differentiation. 618 ThefinalterminEq.(A.6)isthederivativeoftheexternalforcesvectorwithrespecttothedesign 619 variables. The external forces depend on the design through the thickness and the concrete mass 620 density 621 𝜕f𝑠 𝜕ℎ=∑︁ ℓ𝛾𝑐𝐴ℓ 4, (A.16) where𝛾𝑐is the mass density of the concrete, 𝐴ℓis the area of the elements, and the summation sign 622 stands for assembly according the elemental DOF. 623 Appendix A.3. Shear Constraint 624 Similarly to the deflection constraint, the derivative of the shear constraint is, 625 𝜕𝑓2 𝜕˜𝑋=1 ˜𝜎∗ 𝑡𝑠𝜕˜𝜎𝑡𝑠 𝜕˜𝑋. (A.17) 34 Since ˜𝜎𝑡𝑠is an implicit function of the design variables, we use the adjoint approach again. The 626 augmented functional is 627 (˜𝜎𝑡𝑠)𝑎=˜𝜎𝑡𝑠−𝝀𝑇 𝜏(Ku𝑑−f𝑑). (A.18) This time, there is no explicit dependence and therefore after differentiating and replacing the aug- 628 mented functional with the original one, we get 629 𝜕˜𝜎𝑡𝑠 𝜕˜𝑋=𝜕˜𝜎𝑡𝑠 𝜕u𝑑𝜕u𝑑 𝜕˜𝑋−𝝀𝑇 𝜏𝜕K 𝜕˜𝑋u𝑑+K𝜕u𝑑 𝜕˜𝑋−𝜕f𝑑 𝜕˜𝑋 . (A.19) Selecting the adjoint vector such that the terms involving𝜕u𝑑 𝜕˜𝑋will vanish, we get 630 𝜕˜𝜎𝑡𝑠 𝜕˜𝑋=−𝝀𝑇 𝜏𝜕K 𝜕˜𝑋u𝑑−𝜕f𝑑 𝜕˜𝑋 with K𝑇𝝀𝜏=𝜕˜𝜎𝑡𝑠 𝜕u𝑑𝑇 . (A.20) The only term that is unknown is𝜕˜𝜎𝑡𝑠 𝜕ˆu𝑑which is obtained by differentiation of Eq. (20) 631 𝜕˜𝜎𝑡𝑠 𝜕u𝑑=ˆ˜𝜎𝑡𝑠©­ «2𝑁𝑛𝑜𝑑𝑒𝑠∑︁ 𝑗=1𝜎𝑝 𝑡𝑠,𝑗ª® ¬−1 𝝈◦(𝑝−1) 𝑡𝑠𝑇𝜕𝝈𝑡𝑠 𝜕u𝑑, (A.21) where𝜕𝝈𝑡𝑠 𝜕u𝑑isobtainedbydifferentiating(5),multiplyingwith W,andselectingtheappropriateterms 632 𝜕S 𝜕u𝑑=W𝑇DB (A.22) Appendix A.4. Moment Constraint 633 The derivative of the moment constraint is 634 𝜕𝑓3 𝜕˜𝑋=1 ˜𝜇∗𝜕˜𝜇 𝜕˜𝑋(A.23) The augmented moment functional is 635 (˜𝜇)𝑎=˜𝜇−𝝀𝑇 𝜇(Ku𝑑−f𝑑). (A.24) Therelativemomentisrelatedtothedesignvariablesbothimplicitlyandexplicitly, ˜𝜇=˜𝜇˜𝑋,u𝑑˜𝑋. 636 Therefore,wedistinguishbetweenthetotalderivativeandthepartialderivativeoftherelativemoment 637 35 by using different operators notations of 𝑑and𝜕, respectively. Thus, after differentiating the above 638 equation and getting back to the original moment functional we get 639 𝑑˜𝜇 𝑑˜𝑋=𝜕˜𝜇 𝜕u𝑑𝜕u𝑑 𝜕˜𝑋+𝜕˜𝜇 𝜕˜𝑋−𝝀𝑇 𝜇𝜕K 𝜕˜𝑋u𝑑+K𝜕u𝑑 𝜕˜𝑋−𝜕f𝑑 𝜕˜𝑋 . (A.25) After eliminating the derivatives𝜕u 𝜕˜𝑋by finding a proper adjoint vector, the derivative of the moment 640 constraint is 641 𝑑˜𝜇 𝑑˜𝑋=𝜕˜𝜇 𝜕˜𝑋−𝝀𝑇 𝜇𝜕K 𝜕˜𝑋u𝑑−𝜕f𝑑 𝜕˜𝑋 with K𝑇𝝀𝜇=𝜕˜𝜇 𝜕u𝑑𝑇 . (A.26) The explicit derivative can be written in the following form 642 𝜕˜𝜇 𝜕˜𝑋=𝜕˜𝜇 𝜕ℎ𝜕ℎ 𝜕˜𝑋,with𝜕˜𝜇 𝜕ℎ=˜𝜇∑︁ 𝝁𝑝−1𝜕𝝁 𝜕ℎ𝑇 𝝁◦(𝑝−1), (A.27) where◦indicates elementwise operation. The derivative of the relative W&A moments with respect 643 the slab thickness is given by 644 𝜕𝝁 𝜕ℎ=("𝜕M 𝜕ℎ𝑇 M𝑐−M𝑇𝜕M𝑐 𝜕ℎ# ◦M◦−2 𝑐) . (A.28) In the equation above, Mis a vector with all W&A moments at all nodes and M𝑐is a vector with 645 the moment capacities. All W&A moments have similar structure, thus for example the derivative of 646 𝑀𝑟𝑥,𝑚𝑎𝑥is given by 647 𝜕𝑀𝑟𝑥,𝑚𝑎𝑥 𝜕ℎ=𝜕𝑀𝑥𝑥 𝜕ℎ+sign𝑀𝑥𝑦𝜕𝑀𝑥𝑦 𝜕ℎ. (A.29) The derivatives of the plate moments are obtained by differentiating Eq. (5), multiplying with W, 648 and selecting the moments components 649 𝜕S 𝜕ℎ=W𝑇𝜕D 𝜕ℎBˆu𝑑. (A.30) The derivative of the moment capacities is obtained by differentiating Eq. (23), where the only 650 derivative with non zero value is the derivative with respect the slab thickness 651 𝜕M𝑐 𝜕ℎ=1(ℎ−𝑑𝑠)0.64𝑓𝑐𝑑. (A.31) 36 The last component is the derivative of the approximate maximum relative moment with respect to 652 the displacements which is obtained by differentiating Eq. (25) 653 𝜕˜𝜇 𝜕u𝑑=˜𝜇∑︁ 𝝁𝑝−1𝜕M 𝜕u𝑑𝑇 M◦−1 𝑐 𝝁◦(𝑝−1)(A.32) where<·>is a diagonal operator and the derivatives of the nodal moments were computed in 654 Eq.(A.22). AllothercomponentsaregiveninpreviousderivationsoftheSAoftheotherfunctionals. 655 Appendix B. Implementation 656 We solve the optimization problem using a gradient based algorithm due to its efficiency in 657 dealing with large number of design variables. Specifically, the MMA algorithm [36] which is 658 common algorithm in structural optimization. However, a successful optimization requires also 659 several implementational techniques which are described in the following sub-sections together with 660 some related considerations. Thereafter, we summarise all the geometrical data that is used in the 661 examples that are presented in this study. 662 Appendix B.1. Convergence Criteria 663 The basic convergence criterion is related to the change in the objective function. Because the 664 objectivefunctionmighthavenoisybehavior,weconsidertheaveragechangeintheobjectivefunction 665 over the previous 𝑁𝑓0iteration. We define a cumulative convergence parameter 𝑓0𝑐that is promoted 666 each iteration that the change in average objective function is less than 𝑓∗ 0𝑐and demoted otherwise. 667 The objective function is converged when the cumulative convergence parameter is equal to 𝑓0𝑐𝑖. 668 Additionally, we require that at convergence the solution is feasible, such that the maximum of all 669 constraints is less than 𝑓∗=0.01 670 Appendix B.2. Dynamic Move Limits 671 It was observed that the optimization may have oscillatory behavior of the design variables, and 672 as a result the objective function, do not converge. Therefore we implement a dynamic move limit 673 37 mechanismsuchthatthemovelimitofanoscillatingdesignvariableistightenedandthemovelimitof 674 monotonically behaving design variables gets wider. Thus, each design variable has a stability index 675 𝑆𝐼that is promoted each time that the change in design variable value is the same as in the previous 676 iterationanddemotedotherwise. Thestabilityindexofthe 𝑚thdesignvariableatthethe 𝑛thiteration 677 is given by 678 𝑆𝐼𝑛 𝑚=𝑆𝐼𝑛−1 𝑚+signh 𝑋𝑛 𝑚−𝑋𝑛−1 𝑚  𝑋𝑛−1 𝑚−𝑋𝑛−2 𝑚i . (B.1) Oncethestabilityindexofadesignvariablereachesthepositiveornegativethresholdvalues, 𝑆𝐼+and 679 𝑆𝐼−, the move limit is updated accordingly as follows 680 𝑀𝐿𝑛 𝑚=  𝑀𝐿𝑛−1 𝑚𝛼 𝑆𝐼𝑛 𝑚=𝑆𝐼+ 𝑀𝐿𝑛−1 𝑚𝛼 −𝑆𝐼+ 𝑆𝐼− 𝑆𝐼𝑛 𝑚=𝑆𝐼− 𝑀𝐿𝑛−1 𝑚 otherwise,with𝛼>1. (B.2) Additionally, it was observed that the oscillations may occur on a larger scale, where the design 681 variables behave monotonically with respect the neighboring iterations but the optimization fail to 682 converge. In order to deal with this problem we monitor the number of times that the objection 683 function crosses the average objective function at a predefined sampling widow of iterations. Thus, 684 we define a threshold value for the number of intersections between the average and non-average 685 objective functions, beyond which all move limits of all design variables are narrowed down. Herein 686 we consider two sampling windows, representing two different scales of iterations, of 10 and 100 687 iterations and set the threshold value of intersections to 3 and 10 respectively. Thus, each time that 688 anyofthethresholdvaluesisreached,allmovelimitsnarroweddownbyfactorof0.9. Finally,weset 689 minimum and maximum values for the move limits of 1×10−2and1×10−4, respectively. 690 Appendix B.3. Numerical Damping And Continuation Of The Projection Radius 691 It was shown in [28] that the numerical performance of optimization of supports location can be 692 significantlyimprovedbyimplementingthreetechniquespresentertherein. Namely: Controlofinitial 693 38 design,continuationoftheprojectionradiiandnumericaldampingofthederivatives. Inthisstudywe 694 implemented the numerical damping and the three stage continuation scheme of the projection radii 695 as presented in [28]. The initial design control has not been implemented directly, since the initial 696 designs herein are obtained manually and comply with the conditions of the initial control as defined 697 in [28]. 698