Energy-based modelling of in-plane fragility curves for the 2D ultimate capacity of Italian masonry buildings

The high seismic hazard of the Italian territory and the vulnerability of its historic masonry heritage require the development of fragility curves that must be increasingly reliable and robustly correlated to exposure. To date, national-scale seismic risk analyses mainly use empirical curves derived from the statistical analysis of damage induced by past events. These curves have shown good reliability, but they correlate only with a few typological-structural characteristics of the building, such as the number of floors, the vertical structure typology or the construction period. The present research paper aims to overcome this limitation with a hybrid approach that provides a better exposure characterisation. Specifically, the proposed strategy integrates the SAVE and Piecewise Rigid Displacement (PRD) methods. SAVE is an empirical approach based on the damage assessment due to past seismic events used to identify a seismic behaviour of a structure, while the PRD method is a numerical approach that solves the boundary value problem for normal, rigid, no-tension material. It can model different structural typologies, and as a result, it also provides the value of the horizontal static multiplier that drives the masonry construction to collapse. An extended numerical campaign is carried out considering a sample of 750 masonry buildings distributed throughout the Italian territory and extracted from the PLINIVS typological database. Looking at each construction, first, a PRD analysis is conducted to define its seismic capacity, paying special attention to modelling construction details. After that, the SAVE method is used to classify the construction in a specific seismic vulnerability class, i.e., from A to C, with decreasing vulnerability. All the buildings belonging to the same class are then collected, and three fragility curves representative of the collapse state (one for each vulnerability class) are derived and validated against empirical and analytical ones commonly adopted in the Literature. The integrated methodology shows a good agreement between simulations and observations, confirming the viability of the proposed hybrid methodology for the large-scale assessment of masonry buildings, providing an effective strategy to plan mitigation and rehabilitation interventions.


Introduction
Italy is a country exposed to a high seismic hazard as a relevant part is estimated to be subjected to high seismic actions, i.e., with a 10% probability in 50 years of acceleration values exceeding 0.225 g (OPCM et al., 2006). The highest seismic values refer to Calabria, south-eastern Sicily, Friuli-Venezia Giulia, and the central-southern Apennines. Average seismic acceleration values refer to the Salento Peninsula, the Tyrrhenian coast, between Tuscany and Lazio, Liguria, most of the Po Valley, and the entire Alpine Arch. Sardinia is the least dangerous region, as only moderate shaking values are expected. Additionally, the corresponding seismic risk drastically increases as most areas show a vast presence of old civil buildings, either designed with outdated rules or suffering from deterioration due to age. According to ISTAT Census (2001), 60% of residential buildings were constructed before 1980, and 42.5% are over 50 years old.
Additionally, more than one-half of such buildings were built before 1970 without paying attention to any seismic rule, especially for masonry buildings in the historical centres. The combination of the elevated hazard and the vulnerability of civil buildings has already resulted in considerable high damage to building heritage and population, as evidenced by the catastrophic events of the last 50 years. Starting from the earthquake in the Belice area (1968), 15 seismic events with a magnitude greater than 5.5 occurred, and more than 120 billion euros were allocated to interventions due to damage. This condition has strongly fostered the development of structures for managing seismic risk to protect people and building heritage.
Starting with the 1976 Friuli earthquake, the Italian Civil Protection system has consolidated the DRM cycle (preparedness, response, recovery, mitigation) so much that it has become an international reference, boosting scientific research to develop models for the assessment of the seismic vulnerability of buildings. As a result of Civil Protection actions, such effort has been addressed to the development of shared analysis for the national seismic risk assessment (Dolce et al., 2021), according to European Guidelines (Poljansek et al., 2019), in response to the specific requirement of the Sendai Framework for Disaster Risk (UNDRR, 2015). The last National Risk Assessment for Italy dated back to the end of 2018 and was promoted by the Department of Civil Protection (Dolce et al., 2021) through empirical fragility curves based on available damage and vulnerability data for buildings inspected after past seismic events (Irpinia 1980;L'Aquila 2009). The methodology adopted has allowed the seismic risk assessment for the whole Italian territory thanks to a procedure based on the convolution of the hazard (OPCM et al., 2006), exposure (ISTAT Census 2001) and vulnerability. A seismic vulnerability model represents the expected damage due to a given level of Peak Ground Acceleration (PGA) for a building having known typological features. Different approaches can be adopted to construct a vulnerability model (Calvi et al., 2006): in particular, observation-based approaches require previous survey activities in areas that have experienced seismic events. Through specific post-earthquake seismic forms (i.e., AeDES (Baggio et al., 2007), Palazzi (Presidenza del Consiglio dei Ministri, 2006), and Churches (Presidenza del Consiglio dei Ministri, 2013) forms), devoted to rapid visual screening and structural safety checks through expert judgments, the information on the buildings' typological characteristics and level of damage are collected and reported in a damage database. This information is combined with shakemaps provided by geology and volcanology research centres, such as the National Institute of Geophysics and Volcanology-INGV in Italy, to associate a PGA to each surveyed building based on its position. Finally, it is possible to derive the correlation among collected typological characteristics, reached levels of damage and associated hazard. In the Literature, several observation-based approaches have been proposed in the context of the development of seismic fragility curves (the reader is referred to (Benedetti et al., 1988;Riuscetti et al., 1997;Rossetto et al., 2013;Zuccaro et al., 2021) and references therein).
Oppositely, mechanical approaches are based on analytical or numerical structural evaluations, which directly take into account type of materials, real geometries, the presence of reinforcements or any other structural feature. However, a high computational and modelling effort is required to assess several hundreds of buildings. In the scientific Literature, a few simplified approaches have been proposed to lower the computational effort for the large-scale assessment of the seismic vulnerability of masonry buildings. In the context of masonry structures, mechanical-based methods can be grouped into three main categories (Shabani et al., 2021): (i) collapse mechanism-based methods, in which kinematic chains are used to define and evaluate the collapse multipliers corresponding to possible failure mechanisms. Non-linear static or pushover analyses are then adopted to define capacity curves (Bernardini et al., 1990;Lagomarsino and Podesta, 2004;Zuccaro et al., 2017); (ii) capacity spectrum-based methods whose use allows to compute predetermined capacity curves for each building typology. The capacity curve is then intersected with the seismic demand to derive the performance points for different damage thresholds (D'Ayala et al., 2013;Ansal, 2012;Fajifar, 2000); (iii) fully displacement-based methods, where buildings are modelled through an equivalent single-degree-of-freedom system. The displacement capacity for each damage threshold is compared to the displacement demand in each corresponding period of vibration to derive the possibility of crossing the damage thresholds (Calvi, 1999;Crowley and Pinho eBommer, 2004;Borzi and Pinho, 2008;De Angelis et al., 2020).
Several methods have been proposed to provide collapse multipliers for in-and out-of-plane collapse mechanisms of ordinary masonry building stocks (Bernardini et al., 1990;Bernardini et al., 2008a;Bernardini et al., 2008b;Zuccaro et al., 2017;Donà et al., 2020), churches (Lagomarsino and Podesta, 2004;DPCM, 2011) and towers (Torelli et al., 2019) at the building scale. Sophisticated limit analysis-based methods have also been developed to define the most probable failure mechanism of 3D masonry towers at the building scale (Milani, 2019) considering predefined mechanisms such as rocking, Heyman's diagonal cracking, and base shear sliding; an optimisation algorithm finds the most probable mechanism by minimising the corresponding failure multipliers. The kinematic limit analysis theorem performed on NURBS-based elements has been developed to assess the buildings' aggregates and the 'domes' failure behaviour subjected to static horizontal loads (Grillanda and Milani, 2020a;Grillanda et al., 2020). Other strategies based on the kinematical approach include the Discrete Element Method (Cascini and Gagliardo ePortioli, 2018; Frontiers in Built Environment frontiersin.org Malena et al., 2019;Gobbin and Lemos, 2021) and Rigid Block Models (Gilbert and Melbourne, 1994;Baraldi eCecchi, 2017;Portioli and Cascini, 2017;Angelillo et al., 2018), ore, more in general, methods based on energy minimization (Gesualdo et al., 2019;Gesualdo et al., 2020;Fortunato et al., 2022). Other numerical approaches based on the Static theorem have been proposed. The collapse load is evaluated as the maximum load for which at least a statically admissible and equilibrated solution exists. Examples of applications of the Static theorem on masonry buildings are provided in (Monaco et al., 2018;Cusano et al., 2021a;Nodargi and Bisegna, 2021a;Cusano et al., 2021b;Cusano et al., 2021c;Montanino et al., 2021;Montanino et al., 2022). These mechanical-based approaches allow for numerically computing the collapse load of specific masonry structures and are independent of specific seismic events. Analytical methods for seismic fragility curves have also been proposed for different typologies of buildings, in particular, for Reinforced Concrete (RC). The most common approaches are based on non-linear dynamic analyses (Folić and Čokić, 2021;Manfredi et al., 2022;Ruggieri et al., 2022). The proposed methods often require a relevant modelling effort in defining the geometry and a high computational effort to solve the numerical problem. Some of them also require a detailed material characterisation to describe the mechanical response correctly (Monaco et al., 2021). A material description that all too often is impossible to achieve, particularly when looking at large-scale problems. This research aims to fill the gap between observation-and mechanical-based approaches by proposing a novel methodology to define fragility curves that, on the one hand, can directly take into account geometrical and material aspects and, on the other hand, provide a fast strategy for the large-scale assessment of masonry buildings. This approach founds upon a pipeline that combines the SAVE methodology (Zuccaro eCacace, 2015), the extensive PLINIVS database (Cacace et al., 2018) and the Piecewise Rigid Displacement (PRD) method (Iannuzzo et al., 2018;Iannuzzo and VanMele eBlock, 2020). The PRD approach is a fast energy-based method that allows framing and solving the Limit Analysis kinematic problem for notension materials through linear or second-order cone programming. Using a suitable objective function, PRD allows to take any boundary condition into account and to model several typological features. Based on the PRD results of the buildings sample extracted from the PLINIVS database, the SAVE method is then applied to define fragility curves. Specifically, the first step of the pipeline consists in defining a sample of structures from the PLINIVS database, which contains information on the structural-typological features of hundreds of thousands of Italian buildings. This sample comprises 750 buildings extracted through the most recurrent typologies in the PLINIVS databasea. The SAVE procedure is then applied to assign to each building a vulnerability class, which is defined as a function of the building's main structural-typological features and calibrated on damages observed after past seismic events. The analysis is conducted on the 2D main façade of the building, and a parametric description of their geometry is provided to generate digital models. The PRD is then applied to the corresponding digital models to find the pseudo-static collapse loads. Because of the fast numerical PRD solving, several hundred buildings are analysed, and the related results are used to construct fragility curves for each building vulnerability class.
The paper is organised as follows. Section 2 describes the theoretical framework, referring to the SAVE procedure and the mechanical framework on which the PRD method is based; Section 3 introduces the process employed to select the buildings sample and the numerical approach used to discretise the energy problem in a mathematical programming optimisation. Section 4 details geometries, material properties and typological features and how these are framed in a numerical optimisation referring to 2D buildings. The results are then presented and discussed in Section 5. A final section outlines the pro and cons of the present methodology along with future outcomes.

Methodology
The present Section briefly recalls the pipeline at the base of the proposed methodology, whose main structure is depicted in Figure 1. The first step looks at selecting a sample of 750 building typologies from the extensive PLINIVS database, which comprises about 240,000 buildings. All the sample buildings are assumed to be different from each other and have the same frequency of occurrence for each construction typology. Each building is then separately analysed using the SAVE procedure to assign a vulnerability class and the PRD method to define the maximum horizontal pseudo-static multiplier.
Specifically, as detailed in Section 2.1, the SAVE method reduces the uncertainty provided by the European Macroseismic Scale EMS98 (Grünthal, 1998) and is applied to link each building to a specific vulnerability class. Simultaneously, from an analytical standpoint, each building is numerically processed by the PRD method to find the maximum allowable horizontal load. The PRD method is based on a pure Heyman material model, i.e., explicitly considering the nosliding material failure. Indeed, whenever rigid block models are applied to model unilateral continua, such as old masonry buildings, the need for modelling diagonal cracks forces the introduction of a yield criterion to consider the shear-sliding behaviour. This criterion is usually represented by the classic Mohr-Coulomb relation. However, from a computational perspective, when a non-associative behaviour is considered, searching for a solution to the boundary value problem for which normal forces vanish in the presence of normal detachments requires more sophisticated computational and modelling strategies. Nonlinear programming methods (Nodargi and Bisegna, 2021b;Kao et al., 2022;Nodargi and Bisegna, 2022) have to be used for such a purpose, or heuristic, linear or convex algorithm procedures, to lower the computational burden, although the convergence is not mathematically proven.
The collapse multiplier is found as the maximum horizontal load that a building can sustain. It is later translated into a Peak Ground Acceleration (PGA) according to the Italian Technical Code for buildings (Ministero delle Infrastrutture e dei Trasporti, 2018) and used to construct fragility curves.
This Section is organised as follows. In section 2.1, the SAVE methodology and its procedure for the quick assignments of the vulnerability classes are briefly recalled. In section 2.2, the PRD's theoretical framework is illustrated.

SAVE methodology
The SAVE method provides an empirical procedure to quickly describe the overall seismic response of a building through the Frontiers in Built Environment frontiersin.org assignment of the seismic vulnerability classes. The 1998 European Macroseismic Scale, EMS'98 (Grünthal, 1998) considers five vulnerability classes (see Figure 2), denoted with letters ranging from A to E with decreasing levels of fragility, and assigns to each building the most likely vulnerability class as a function of the vertical structure typology, neglecting the influence of the other typological-structural characteristics (age, floor number, horizontal structure, etc). The SAVE procedure reduces the level of uncertainty on the attribution of vulnerability classes considering the effects on the building seismic response of additional typological features, defined modifiers, such as the horizontal structural typology, the number of floors, the presence of ties and other features. The method is empirical and built on the extensive PLINIVS database in which typological and damage information of about 250,000 masonry buildings are collected. In the first step, the correlation between the vertical typologies and the level of damage is estimated. The vertical structures defined in the SAVE method are three: generic masonry, denoted with V 0 , characterized by the absence of information on the quality of the wall structure; V 1 , weak and irregular masonry; V 2 , regular and good quality masonry. The levels of damage follow the EMS'98 scale and are: D0, that is, no damage; D1, light, non-structural damage; D2, light structural damage; D3, high structural damage; D4, partial collapse; and D5, global collapse. Referring to the vertical structure V i , SPD Vi is the Synthetic Parameter of Damage index, and it is computed as the barycentric abscissa of the damage distribution (Zuccaro eCacace,

FIGURE 1
Flowchart summarising the methodological pipeline at the base of the present research.

Frontiers in Built Environment
frontiersin.org 04 2015). As an outcome of this step, three SPD ranges representing the vertical structures are defined (Table 1): class A includes most of the buildings with weak and irregular masonry; class B includes most of the buildings with generic masonry; and class C includes most of the buildings characterized by good quality masonry.
In the second step, the influence of the parameters affecting the average seismic behaviour of the building is estimated through the Synthetic Parameter of Damage (SPD) variation. To this purpose, buildings with vertical structure V i , typology k of the j parameter are classified with the number SPD ViP jk , where j refers to the horizontal structure and k to the floor deformability. The influence of the modifier k on the parameter P j on the vertical structure V i is defined as the difference between the SPD ViPjk and SPD Vi . It is essential to point out that some typological parameters that characterise the seismic response of the building could be mutually dependent. For example, the presence of ties within a horizontal structure may be recurrent in the case of deformable floors. When calculating the influence of the presence of tie rods and deformable floors, it is, therefore, possible that the SPD ViP jk value is calculated on the same sample for both parameters. To avoid overestimation when determining the influence of a parameter, it is necessary to define its correlation with the remaining dependent parameters. For parameters independent of the others (e.g., the location of the building in the aggregate), it is not necessary to evaluate correlation factors. After that, once the typological characteristic of a building is known, the corresponding SPD can be calculated as follows: where q is the influence of each independent parameter, p is the influence of the dependent parameter, n is the total number of independent parameters, m is the number of dependent parameters, and c ij is the correlation coefficient between the classes of parameters p i and p j .

Boundary value problem and limit condition
This Section illustrates the boundary value problem (BVP) at the base of the PRD method referring to a no-tension Heyman material. Specific mathematical restrictions on the stress and latent strain are introduced to model the masonry response. A 2D masonry structure is modelled as a continuum occupying the region Ω of the Euclidean space E 2 . T denotes the stress, u the displacement of the material points x ∈ Ω, while E the infinitesimal strain field assuming small displacements. The Heyman model is enforced through the so-called Normal Rigid No-Tension (NRNT) material restrictions: with Sym − , Sym + the mutually polar cones of negative and positive semidefinite symmetric tensors. Equations 2 are equivalent to the wellknown normality rule: representing the necessary conditions to rigorously apply the classic theorem of limit analysis to unilateral masonry materials. For more information, the reader is referred to (Kooharian, 1952;Livesley, 1978;Giaquinta and Giusti, 1985;Fortunato et al., 2014;Fortunato et al., 2016;Chiozzi et al., 2017). Therefore the BVP on domains made up of NRNT material reads with n is the unit outward normal to the boundary zΩ partitioned into its constrained zΩ D and loaed zΩ N parts. A BVP solution is a triplet In the present formulation, only singular strains and stresses are considered, meaning that the displacement u and the stress vector s can admit jump discontinuities. The sets K and H of the admissible displacements and stresses are defined as follows: with S, S′ two suitable function spaces. For more mathematical details about the regularity of these spaces, the interested reader is referred to (Giaquinta and Giusti, 1985). Thereafter, only potential jumps whose supports are lines lying within the domain are considered. Therefore, stress and strain fields are everywhere zero except for those curves, where they are modelled with Dirac delta lines. The BVP can be solved using two possible variational strategies, i.e., either solving the equilibrium problem through the minimisation of the complementary energy or the kinematic problem by minimising the total potential energy (TPE). The PRD method uses the latter through a straightforward displacement approach, which results in the search for a displacement u ∈ K for which there exist a stress T ∈ H such that T · E(u) 0. The TPE for an NRNT material reads: with u ∈ K, i.e. the space of kinematically admissible displacements. This function depends exclusively on the displacement u, and for the case at hand, reduces to the potential energy of the external loads only. Notably, the minimiser u°, i.e., the element (if unique) or the elements of K on which the potential energy attains the minimum value (if bounded) implicitly guarantees the equilibrium of the loads imposed on the structure. From a mathematical standpoint, the minimiser u°can also be not unique, which means that multiple energy solutions to the same BVP can coexist. For detailed information, the interested reader can refer to (Giaquinta and Giusti, 1985) and (Anzellotti, 1985), where the existence of the minimum is proven for Elastic Normal No-Tension materials assuming as kinematically admissible space the set BD(Ω) and an additional restriction on the given loads, i.e., the safe load condition. For a discussion about the safe load condition and its numerical applications, the reader is referred to the CDF method (Iannuzzo et al., 2021). To the scope of the present paper, it must be noted that the TPE ℘(u) is always bounded from below as long as the load is compatible, i.e., the set H is not void. This concept provides a strategy to apply the TPE to find collapse multipliers. As it will be shown, it will also reconnect its use to classic lower-bound approaches.
Indeed, the main idea is to describe the external loads as an affine distribution so they remain proportional to the original self-weight through a scalar parameter. From a physical standpoint, many experimental tests are commonly performed through a tilting test to look at possible collapse mechanisms (DeJong, 2009;Ochsendorf, 2022). The numerical modelling of those tests is straightforward and can be performed through two different strategies. The first one is a pure geometrical strategy, and it consists of directly and fictitiously rotating the original angle geometry of an angle α until the structure collapses. An equivalent strategy to model horizontal actions is to consider an additional distribution of horizontal pseudo-static forces b ⊥ such that the total volume forces are: with (·) ⊥ the linear operator that clockwise rotates a vector of π/2, and λ is the scalar load multiplier. This strategy is equivalent to classic pushover analysis. In this regard, the value of the horizontal scalar multiplier represents the gravity loads' percentage applied on a horizontal plane that drives the structure to its last stable equilibrated configuration. This value λ c is called a pseudo-static horizontal collapse multiplier. The previous formulation refers only to body (mass) loads b. However, it must be extended whenever additional loads or tractions s depending on the mass are also involved. Because of (9), TPE reads: The last stable configuration, along with the collapse multiplier, can be defined through the following min-max problem: Because of the uniqueness of the multiplier due to the normality rule, this min-max problem can be solved following two strategies: the first one translates this in classic limit-analysis approaches, while the second one solves the problem through a sequence of minimisation problems until a solution can still be found. The present contribution follows this last procedure as it directly considers different typological features in the objective function as specified in Section 4.

Numerical methodology
The present Section illustrates the statistical procedure adopted to generate the buildings' sample from the data from the PLINIVS database to define the buildings sample (Section 3.1) along with the numerical strategy used to frame problem (11) in a linear programming problem (Section 3.2).

Sample of buildings
Following the SAVE method, the most relevant parameters in the vulnerability class assignment are the vertical and the horizontal structures, the presence of ties, the number of floors and the construction age. The first four parameters are used to generate the sample of virtual façades. Conversely, the last one is related to the exposure and cannot be directly considered to construct the analyticalfragility model. The construction age is a fundamental parameter when approaching the problem from an observational standpoint, as it often provides sufficient information about construction technology. When adopting an analytical approach, the construction technology is explicitly modelled, and the importance of including the construction age in the model vanishes. The four SAVE parameters used to scan the variability of the building stock under consideration: 1) vertical structure, which is defined in terms of hollow bricks, tuff, filled brick, irregular stonework, and regular stonework; 2) horizontal structure, which can vary among steel, reinforced concrete (RC), wooden floors as well as vaulted floor systems; 3) tie rods, whose presence is assumed as a binary variable, i.e., they can be considered or not considered at all. Moreover, when they are modelled, they are not considered at each floor level; 4) the number of floors, which ranges from 1 to 5.
From a technological perspective, not all horizontal structures are compatible with tie roads, and thus, parameters 2 and 3 are collected in a unique variable called horizontal technology. Therefore, these four parameters allow the definition of the following six horizontal technology classes: steel, RC, wooden with ties, wooden without tie rods, vaults with tie rods and vaults without tie rods. In addition to the SAVE parameters, other geometrical variables are considered, such as the façade baseline length and thickness, inter-storey height, and openings, whose size and numbers are assumed as variables. The inter-storey height has been fixed in all the analyses at 3.5 m. The façade baseline length ranges from 4.00 m to 7.00 m, which statistically includes most of the 2D masonry geometries catalogued in the PLINIVS database. The wall's thickness has been considered as a function of the number of floors, and its variation over the height goes from 0.60 m to 0.45 m, as detailed in Section 4.1. Table 2 collects all the typological parameters used to define the building stocks. As the reader can note, combining all of them, the number of façades under consideration is 750.  For each of these typologies, the following information has to be defined: the vulnerability class on the base of the SAVE parameters, the collapse multiplier, which will be found through numerical simulations, and the frequency of occurrence in a typological database.
Looking at the vulnerability class assignment, the exploitable parameters for the SAVE method are the vertical structure, the horizontal technology, and the number of floors. Based on these parameters, Table 3; Table 4; Table 5 report the vulnerability class assignments for buildings with 1-2, 3-4 and 5 floors, respectively. These values are used to develop the fragility curves as reported in Section 4.

Limit solution to the BVP by way of the PRD method
The approximate solution to the minimum problem (11) is obtained by restricting the search of the minimum in the class K pr of piecewise rigid displacements. This infinite-dimensional space is discretised by considering a finite partition of Ω (see Figure 3) into a number M of rigid pieces, such that P(Ω i ) being the perimeter of Ω i . In particular, restricting to convex polygonal elements, the boundary zΩ i of the n-polygon Ω i , is composed of n segments Γ, of length , whose unit normal and tangent vectors are denoted n, t, while the endpoints with 0 and 1. The edges shared by two adjacent elements or lying on the constrained boundary are called interfaces and are collected in the set Γ. K M pr represents then the finitedimensional approximation of K pr generated by this partition. The discretised version of the minimum problem reads: A generic piecewise rigid displacement u ∈ K M pr is in a one-to-one correspondence with the vector U of 3M components representing the 3M rigid body parameters. These parameters are restricted by the assumption that the strain must be positive semidefinite. For piecewise rigid displacements, the strain coincides with its singular part, namely: E s being concentrated along the interfaces among blocks, that is, within the present approximation, over the set Γ. From (4), it follows that sliding is not allowed on any interface and only detachment is possible: In Equation 16, [u] is the retive displacement between two element nodes of the discretisation. Notice that conditions (16), derived from the normality assumption, represent a condition of unilateral contact with no sliding among blocks. With N the number of the interfaces Γ, and v(0), v(1), w(0), w(1) the normal and tangential components of the relative displacements on the endpoints 0, 1 of any segments belonging to Γ, restrictions (16) are equivalent to the 2N inequalities and the 2N equalities Assuming homogeneous boundary conditions, restrictions (16) can be expressed in matrix forms as a function of the vector U collecting the unknown Lagrangian parameters and can be rewritten: Finally, the discretised version of the minimum problem (14), which approximates the minimum continuum problem (8), reads: where the objective function representing the potential energy of the external body loads b is minimised within the set: The limit analysis problem can be then solved through a sequence of LP problems as:t

FIGURE 2
Geometric features of three-story buildings with one (A) or two (B) openings per floor.

Frontiers in Built Environment
frontiersin.org 08 The minimisation problem (20) transforms the original minimisation problem (8) for a continuum, into a minimisation problem for a structure composed of rigid parts, acted on by given loads and given settlements and subject to unilateral contact conditions along the interfaces. Problem (20) is a standard linear finite-dimensional minimisation problem, as both objective function and constraints are linear. The existence of the solution of this approximate problem is trivially guaranteed if the original problem is bounded from below. For small problems, it can be solved exactly with the simplex method (Dantzig et al., 1955), while for large problems, interior-point algorithms represent efficient and fast alternatives (Mehrotra, 1992;Vanderbei, 2015;Dantzig, 2016).

Numerical implementation
The present Section illustrates the numerical pipeline adopted to define the fragility curves combining the SAVE and PRD methods. In particular, Section 4.1 provides a detailed overview of how masonry buildings coming from the PLINIVS database have been modelled, considering geometrical and mechanical information. Section 4.2 shows how these features have been modelled and, after that, framed in an LP problem to define the collapse load multiplier for each case. It is worth noting that the present analysis looks at the in-plane collapse load of the main 2D façade of a masonry building.

Geometrical and mechanical description of the buildings' sample
The present Section gives an overview of the procedure adopted to parametrically model geometries and mechanical features of the buildings' stock extracted from the PLINIVS database. The 2D analysis of a building is carried out by referring to its main façade. For example, looking at three-story buildings, Figures 4A, B describe the geometries and loading conditions of two facades having one and two openings per floor, respectively.
The inter-story height h is assumed to equal 3.5 m for all the façades independently of the floor number, while the openings' dimensions, as mentioned above, are fixed at 1.2 m x 2.5 m. The façade's depth d w varies over the building's height H b according to the following relation Table 6 reports geometric data needed to generate a masonry building of the sample referring to the symbols introduced in Figure 4. It is worth emphasising that these geometrical features as well as the openings number, do not depend on the building height but only on the building baseline length. Indeed, structures with baseline lengths equal to 4 m show only one opening per floor; the ones whose length is equal to 6 m or 6.6 m, show two openings per floor; and façades with a baseline length equal to 5 m can have one or two openings per floor as illustrated in Table 6.
Regarding the loads, mass density values depend on the masonry typology as reported in Table 7.
The presence of slabs and the corresponding transmitted live loads are modelled through additional linear loads placed as in Figure 4, and whose gross values are reported in Table 8. The liner loads adopted in the numerical analyses can be derived by multiplying these values by the slab's net orthogonal length, assumed to equal 5 m.
Lastly, the presence of tie roads depends on the number of floors and, according to the PLINIVS database, whenever they are present (i.e., for wooden and vaulted slabs only), one assumes: -One-, two-and three-story façade has only one tie-rod at the uppermost level; -Four-story facades have two tire rods at the second and fourth levels; and, -Five-story facades show two tire rods at the third and fifth levels.

FIGURE 3
The infinite-dimensional space K pr of piecewise rigid displacements having Ω as support (A) is discretised by partitioning the domain into convex polygonal elements (B). In (B), the partition of the domain Ω into M squares generates the discretised set K M pr .
TABLE 6 Summary of the geometric features needed to generate a masonry façade of the sample according to Figure 4. The symbol * refers to facades with two openings per floor ( Figure 4B).

FIGURE 4
The two strategies adopted to model horizontal loads in the presence of deformable (A) or rigid slabs (B).

FIGURE 5
Overview of five geometrical models and corresponding discretisations used to describe three-story masonry buildings. Yellow stripes denote linear loads due to and transmitted by the slabs, while red lines show the presence of tire rods. Buildings with reinforced concrete and steel slabs are assumed not to show any tire rod systems. The tie-rod strength (T r ), i.e., the maximum stress it can withstand depends on the strength of the cable (T c ) and of the anchor plate (T a ) as well as on the compressive resistance of the masonry (C a ), as from the following relation: It is assumed that the square anchor plate edge length is 0.30 m, the cable has a net diameter of 30 mm, and it is made up of steel with a yield tension of 240 MPa. Moreover, the calculation assumes a confidence factor of 1.35 according to the Italian national code (Ministero delle Infrastrutture e dei Trasporti, 2018) and the masonry compressive strength values as detailed in Table 9.
Based on these assumptions, once the masonry typology has been fixed, Eq. 23 yields the maximum tensile force that the tie-rod can sustain and that will be used in energy minimisation, as detailed in the next Section.

PRD numerical modelling of the buildings' sample
The minimum energy problem is formulated as an LP problem whose most general version reads:

TABLE 10
The number of elements adopted to discretise the buildings' geometries. The symbol * refers to buildings with a horizontal length of 5 m but with two openings per floor. In red are the numbers of digital models elements depicted in Figure 6. PRD-based fragility curves for the three vulnerability classes as from the SAVE procedure.

FIGURE 8
Approximation of the PRD-based fragility curves through lognormal functions.
Frontiers in Built Environment frontiersin.org with b g the vector modelling volume forces as lumped to the blocks' centroid, q s the vector representing the live loads transmitted by the slabs and translated to the blocks' centroid by the operator B q . It is worth pointing out that Q ⊥ s is the vector representing the seismic effects of the slabs' loads, as detailed below; B Ft F tr is the vector collecting the forces exerted by the tie rods and translated to the block's centroids closest to the tie-rod endpoints, while ΔU collects the unknown relative displacements between these blocks. The term B Ft F tr is implemented only on digital models with deformable floors, such as wooden or vaulted floor systems, when tired rods are explicitly considered. It should be mentioned that the horizontal rotated component of the loads Q ⊥ s is modelled following two strategies depending on the slab's typology Figure 5. For deformable floors (i.e., vaulted and wooden floor systems), it results: while for rigid slabs (i.e., steel and RC): where F ⊥ q s = L q ⊥ s q ⊥ s is the force vector that lumps the linear horizontal load q ⊥ s through the operator L q ⊥ s and is concentrated on the blocks' centroids placed on the left side of the building through the operator B F ⊥ q s . Looking at the constraints, the first enforces homogeneous boundary conditions, the second inequality the non-overlapping relation and the third one the no-sliding conditions.
Referring to three-story masonry buildings, Figure 6 shows five digital models and the corresponding discretisations for different façade baseline lengths. Particularly, the number of elements ranges from 6,194 for a baseline of 4 m, to 9,862 for a baseline of 6.6 m. Table 10 shows the number of elements used to discretise the digital models as a function of the building's baseline length and floor number. The robustness of the proposed discretisation against two different numerical approaches is validated here (Iannuzzo et al., 2021).
The computational time required to solve a single optimisation problem is 0.15 s using an agile laptop with an AMD Ryzen 7 5700U and using Mosek (Mosek, 2010) as a solver.

Results and discussions
PRD results, obtained in terms of load multipliers, are translated into peak ground acceleration (PGA) according to the procedure detailed in (Ministero delle Infrastrutture e dei Trasporti, 2018), as

FIGURE 9
Comparison between PRD-based fragility curves (dashed lines) and the ones (solid line) related to the first damage mechanism proposed in (Zuccaro et al., 2017).

FIGURE 10
Comparison among the PRD-based fragility curve for vulnerability class B with the ones proposed in (Cattari et al., 2014) for different damage levels. The PRD results match with the damage state D5.
Frontiers in Built Environment frontiersin.org with q the ductility factor, here assumed as q 2, S the subsoil factor, fixed to 1.25, and g 9.81 m/s 2 the gravitational acceleration. Results are then depicted in Figure 7 in terms of Damage State Probability (DSP) as a function of the three vulnerability classes assigned by the SAVE procedure. The DSP index represents the cumulative percentage of all structures whose collapse is activated by a given PGA value.
The PRD results obtained are then approximated through lognormal curves and the corresponding smooth curves are represented in Figure 8, again as a function of the three vulnerability classes. Equation 28 is the cumulative of a normal distribution with average value μ and standard deviation σ. In the present research, the three vulnerability classes are defined by the following values: In Figure 9, the PRD-based fragility curves are compared with the ones proposed in (Zuccaro et al., 2017), which evaluate the minimum acceleration triggering the first damage mechanism selected among predetermined mechanisms (including also out-of-plane ones) on a sample of buildings from the PLINIVS database. In this sense, it is worth noting that the PRD-based curves are associated with higher PGA levels since the curves evaluated in (Zuccaro et al., 2017) are strongly influenced by the out-of-plane mechanisms occuring at lower levels of PGA.
Restricting to vulnerability class B, shows the PRD-based fragility curve together with the ones proposed in (Cattari et al., 2014) for different damage levels. As expected, the numerical-based curve is in good agreement with the damage level D5, corresponding to the global collapse of the structure. Indeed, as applied, the PRD method can only detect the mechanism driving the structure to collapse without the possibility of accounting for lower damage levels.

Conclusion
The present paper proposes a novel analytic methodology to define masonry structures' fragility curves, employing a hybrid approach that combines the SAVE and PRD methods. SAVE is an observation-based approach to estimating the vulnerability class of buildings based on actual damages after past seismic events. It provides a quick procedure to identify the structure's seismic behaviour. The PRD method is a numerical approach to solving boundary value problems for normal, rigid, notension materials. These two methods were combined on an extensive sample of 750 masonry buildings extracted from the PLINIVS database, which collects typological-structural information on about 250,000 ordinary masonry buildings distributed throughout the Italian national territory. Each masonry structure of the sample was modelled referring to its 2D main façade. After that, the SAVE procedure was applied to assign a vulnerability class to each of them, while a PRD numerical campaign was conducted to compute the corresponding collapse loads. Specifically, PRD digital models were constructed to account for the relevant typological-structural and geometrical features. Lastly, the PRD and SAVE results were combined to produce fragility curves related to the three vulnerability classes as from the SAVE method.
The results of the proposed approach were benchmarked against other numerical strategies based on predetermined mechanisms. The comparison showed a good agreement between the PRD results and the fragility curves present in the Literature, confirming the robustness and effectiveness of the proposed approach, particularly when referring to the damage level D4-D5, as from EMS'98. The estimated curves can be used for the National Risk Assessment in Italy and can also be adapted for the seismic risk evaluation at different scales. Indeed, the main advantage of this methodology is its possibility to analytically evaluate fragility curves considering specific structural and geometric characteristics, which, combined with the ease of modelling and the solving time, allows it to be used on a large-scale assessment of masonry structures. Moreover, the proposed procedure can also be used to localise fragility curves modelling the seismic risk assessments of specific local areas (i.e., local, regional) accounting explicitly for the frequency of occurrence of the related typological-structural and geometrical features by also employing other databases (Zuccaro et al., 2023). This strategy will be pursued in forthcoming contributions while also validated against other empirical fragility curves in the Literature.
However, the current methodology accounts only for 2D in-plane collapse mechanisms. In contrast, out-of-plane ones, as well as information on lower levels of damage, such as D1/D3, cannot be provided. Further development will include extending the numerical approach to overcome these limitations while exploiting data provided by other databases.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Funding
This research was supported by the Italian Ministry of University and Research through the Programme "Rita Levi Montalcini for young researchers" (FFO 2020).