Seismic Vulnerability Assessment from Earthquake Damages Historical Data Using Constrained Optimization Technique

Affiliation(s)

^{1}
BRGM (French Geological Survey), Orléans, France.

^{2}
Materials Sciences and Environment Laboratory, Hassiba Ben Bouali University of Chlef, Chlef, Algeria.

Abstract

This present work falls within the context of efforts that have been made over the past many years, aimed in improving the seismic vulnerability modelling of structures when using historical data. The historical data describe the intensity and the damages, but do not give information about the vulnerability, since only in the ’90 the concept of vulnerability classes was introduced through the EMS92 and EMS98 scales. Considering EMS98 definitions, RISK-UE project derived a method for physical damage estimation. It introduced an analytical equation as a function of an only one parameter (Vulnerability Index), which correlates the seismic input, in term of Macroseismic Intensity, with the physical damage. In this study, we propose a methodology that uses optimization algorithms allowing a combination of theoretical-based with expert opinion-based assessment data. The objective of this combination is to estimate the optimal Vulnerability Index that fits the historical data, and hence, to give the minimum error in a seismic risk scenario. We apply the proposed methodology to the El Asnam earthquake (1980), but this approach remains general and can be extrapolated to any other region, and more, it can be applied to predictive studies (before each earthquake scenarios). The mathematical formulation gives choice for regarding, to the optic of minimizing the error, either for the: 1) very little damaged building (D0-D2 degree) or 2) highly damaged building (D4-D5 degree). These two different kinds of optics are adapted for the people who make organizational decisions as for mitigation measures and urban planning in the first case and civil protection and urgent action after a seismic event in the second case. The insight is used in the framework of seismic scenarios and offers advancing of damage estimation for the area in which no recent data, or either no data regarding vulnerability, are available.

This present work falls within the context of efforts that have been made over the past many years, aimed in improving the seismic vulnerability modelling of structures when using historical data. The historical data describe the intensity and the damages, but do not give information about the vulnerability, since only in the ’90 the concept of vulnerability classes was introduced through the EMS92 and EMS98 scales. Considering EMS98 definitions, RISK-UE project derived a method for physical damage estimation. It introduced an analytical equation as a function of an only one parameter (Vulnerability Index), which correlates the seismic input, in term of Macroseismic Intensity, with the physical damage. In this study, we propose a methodology that uses optimization algorithms allowing a combination of theoretical-based with expert opinion-based assessment data. The objective of this combination is to estimate the optimal Vulnerability Index that fits the historical data, and hence, to give the minimum error in a seismic risk scenario. We apply the proposed methodology to the El Asnam earthquake (1980), but this approach remains general and can be extrapolated to any other region, and more, it can be applied to predictive studies (before each earthquake scenarios). The mathematical formulation gives choice for regarding, to the optic of minimizing the error, either for the: 1) very little damaged building (D0-D2 degree) or 2) highly damaged building (D4-D5 degree). These two different kinds of optics are adapted for the people who make organizational decisions as for mitigation measures and urban planning in the first case and civil protection and urgent action after a seismic event in the second case. The insight is used in the framework of seismic scenarios and offers advancing of damage estimation for the area in which no recent data, or either no data regarding vulnerability, are available.

Keywords

Seismic Damage Scenario, Damage Database, Vulnerability Assessment, Inverse Optimization, Model Calibration, El Asnam Earthquake, RISK-UE

Seismic Damage Scenario, Damage Database, Vulnerability Assessment, Inverse Optimization, Model Calibration, El Asnam Earthquake, RISK-UE

1. Introduction

Optimization algorithms have started to be widely used in the field of earthquake engineering because of their capability to calibrate the model when no data are available for important parameters associated to behavior of the system (e.g. a system can be a building, a bridge, a network…). The first utilization of the optimization algorithms has been at the scale of one structure generally for fitting the model parameters (e.g. as the Young modulus, reinforcement bar diameter, compression resistance…). Several studies have also focused on effective computational methods of cost optimization for bridges [1] [2] [3] and structures in general [4] . Recently, [5] used optimization process for the design of structures. The application of the framework proposed by them is illustrated in a bridge design optimization problem where the column reinforcement bar diameter and concrete cover are considered as design parameters. [6] presents a hybrid optimization methodology for the probabilistic finite element model updating of structural systems. In paper [6] , the model updating process is formulated as an inverse problem, analyzed by Bayesian inference, and solved using a hybrid optimization algorithm.

Even if the optimization process is widely used at the scale of one structure, it is less common to apply it at the scale of a group of buildings, networks or the scale of town analysis. More recently, the optimization algorithms were also used for the natural hazard studies at large scales, and more particularly, for manage the decision after the catastrophic event. [7] proposes a simulation model to find out optimum evacuation routes, during a tsunami using Ant Colony Optimization (ACO) algorithms. [8] proposes a framework that clarifies the interrelationships between notions of coping capacity, preparedness, robustness, flexibility, recovery capacity, and resilience, previously espoused as independent measures, and provides a single mathematical decision problem for quantifying these measures congruously and maximizing their values. [9] proposes a fuzzy multi-criteria model to deal with “qualitative” (unquantifiable or linguistic) or incomplete information and illustrates it within the post-earthquake reconstruction problem in Central Taiwan, including the restoration concerning the safe and serviceable operation of “lifeline” systems, such as electricity, water, and transportation networks, immediately after a severe earthquake.

The seismic risk scenario is now one of the most powerful tools for assessing damages either for prevention and mitigation aims or for crises management situation. Many studies [10] [11] [12] [13] , more or less detailed, are published on this subject, which demonstrate the demand of stakeholders and first responders (e.g. civil protection) to have a better understanding and evaluation of damages and the various outcomes from these kind of studies. Different methodologies are applied for the assessment of the earthquake damages, which are generally based on two mandatory terms: the level of hazard and the level of vulnerability of the exposure (i.e. element at the risk). A very complete state of art of the development of seismic vulnerability assessment methodologies for variable geographical scales over the past 30 years is presented in [14] . Different methodologies are facing over the years: the widest used methodologies based on a vulnerability index [15] [16] [17] [18] ; pushover-based vulnerability analysis [19] ; displacement-based vulnerability assessment procedure [20] and buildings’ vulnerability assessment using the parameter less scale of seismic intensity [21] . [22] [23] provide a detailed description (Summary of Software, Methodology, IT Details, Exposure Module, Hazard Module, Vulnerability Module and Output) of the existing software for seismic risk assessment that is either open source or has been made available to the GEM Risk Team. Here we are dealing with empirical methodologies based on the vulnerability index obtained from observational data after earthquakes events. These methods are very useful for the representation of the vulnerability at large scale (see [13] for more details of the used methodology).

Concerning empirical methods, starting from the first “Earthquake damage probability matrices” derived by [24] after the saint Fernando earthquake and presented to the “Fifth World conference on earthquake engineering”, all the empirical methods for assessment of the seismic damage of structure are based on the relation between the observed intensity and the observed damage. One key moment in the collecting and the utilisation of the Earthquake damage data is the development of the European Macroseismic Scale―EMS98 [25] which is the first intensity scale that clearly defines the concept of vulnerability. Table 1 gives the classification of damage in EMS98 according to masonry and reinforced concrete buildings. Before this scale, the previous scales Mercalli-Cancani Sieberg Scale―MCS [26] , Modified Mercalli scale (MM-31 and MM-56) [27] and Medvedev-Sponheuer-Karnik scale (MSK-64) [28] , do not clearly correlate the vulnerability of structure with the intensity and the damage scale. One of the main objects of the EMS98 scale was the fact to be consistent with the pervious scale; hence, that data collected using previous scales can be adapted to the definitions of the EMS98 scale. The data sets are a very laborious and time-consuming task, but is the key of the all these methods, and their representatively depend on the accuracy of the collected data. In Italy, the concerted effort to collect earthquake damage data over the past 30 years has led to the development of an extensive database from which vulnerability predictions for the Italian building stock can be derived [29] . In France, a database which repertory the historical earthquake [30] [31] is filled, but the information related to the damages are quite thin. Some other databases exist, developed by other countries, or within

Table 1. Classification of damage to masonry and reinforced concrete buildings (according to EMS98).

framework of some European Commission projects, but their accessibility, especially to original “raw” data, is fastidious.

[32] systematically compared statistical modelling techniques for different empirical datasets and explored many of the issues raised regarding the treatment of uncertainty.

Database typologies and their typical issues depend on the manner on which the observations were obtained they were been obtained: for Detailed “Engineering” Surveys and Surveys by Reconnaissance Teams the main issue is the possibility of unrepresentative samples; for the Rapid Surveys methods, generally the evaluation concern the habitability or the safety not the evaluation of the damage; while for remotely sensed survey method, which is quite new and not already able to efficiently evaluate other damage degrees than the collapse or very heavy damage [33] . For all these cases it is quite rare to have the three essential terms at the scale of the building or census block: intensity measure, building typology, and hence it vulnerability, and damage degree. Moreover, once this “original/unchanged” dataset is defined, usually data manipulation and combination were done in order to associated to some related parameters and develop new method for vulnerability assessment.

Hence, in our paper we are using a literature available “original/unchanged” dataset that represent the large majority of the available historical dataset, namely dataset that have only the intensity and the damages degrees by census block. Our paper proposes a new mathematical manner of treating the existing datasets, in which vulnerability information are not available, and of estimating the best vulnerability index and its characteristics to be use in the seismic risk scenarios using optimization procedure.

2. Historical Case Study and Available Observed Data on Damages

A retro-scenario of the October 10, 1980 El Asnam earthquake was performed by the authors [34] [35] using Armagedom tools [13] developed by the French Geological Survey (BRGM). In this paper, we are just using the intensity and damage observations in order to find the vulnerability index that will give the smaller error between the observed damage and the damages calculated numerically with Armagedom software.

Chlef city (formerly named El Asnam) situated 200 km west of the capital of Algeria Algiers, is an area that suffers of important seismic activities due to the interaction of the two Eurasian and African plates. Major earthquakes rocked the city during the last century, the last event that caused major damages being the one of October 10, 1980; known as El Asnam earthquake and killed around 3000 persons. This event was qualified as the largest earthquake ever known in the west-Mediterranean region. Different international studies [36] [37] [38] were performed based, or in collaboration, with the Algerian Technical Inspection of Construction (CTC) that conducted a large field investigation and damage inventory in El Asnam city on 5131 buildings. Figure 1 presents the ten sectors in which the city was divided after the CTC’s survey (adapted from [37] ). Based on this data, [36] established a buildings damage classification, in which the damage levels are: Green―very little damage. Can be occupied immediately; Orange: needs further study before it can be either occupied or condemned; Red:

Figure 1. Identification of the ten El Asnam sectors for CTC survey (adapted from [37] ): collapse―black solid fill, very heavy damage―dark checkerboard pattern, heavy damage―dark diagonal line pattern, moderate damage―lighter diagonal line pattern and slight or no damage―lighter checkerboard pattern.

Table 2. Damage classification of El-Asnam buildings by [36] . ID: district Id. I: EMS98 intensity. Ntot: total number of buildings in percentage. Gr: green damage grade (D0 to D2). Or: orange damage grade (D3). Re: red damage grade (D4 to D5). Un: undefined damage grade.

condemned and should be demolished. Table 2 presents the number of buildings and the damage classification (number of red, orange and green) performed by [36] in each of the ten sectors of the city. The seismic vulnerability assessment of the existing buildings has never been done either.

The data given in Table 2 represents the collected observational dataset. This representation of information is the very common for post-earthquake damages assessment after seismic event. This concern ancient’s earthquake as well as the most recent earthquake (Aquila 2009), for which intuitively we can imagine to have more detailed information’s, but in reality, this is not the case.

3. Model Calibration

3.1. Model Description

Figure 2 schematize the general methodology for seismic damage calculation. It presents each of the most important phases necessary for creating an earthquake scenario: regional seismic hazard, local seismic hazard, etc. In the case of El Asnam, the hazard modules (module 1, 2 and 3) are not performed, since we got the observed macroseismic intensity. Here all performed simulations (9800 runs) were executed on the Armagedom tools [13] developed by the French Geological Survey (BRGM) and used for simulation of the damage scenarios.

For damage assessment, in Armagedom software the procedure developed within RISK-UE is used [39] . Firstly, the vulnerability to earthquake shaking of exposed elements at risk (for this study, buildings) are characterized by vulnerability indices (V_{i}), which range from zero (not vulnerable) to one (building is

Figure 2. General methodology and modulus in Armagedom software for seismic damage calculation.

highly vulnerable). Depending on the buildings proprieties (materials, age of construction, height…), these last are regrouped in categories; commonly known as typologies (see [40] for more details). Therefore, a value of V_{i} is assigned to each building typology. For each typology, the mean damage degree (µ_{D}: between zero and five) is estimated based on a vulnerability function:

${\mu}_{D}=2.5\cdot \left[1+\mathrm{tan}h\left(\frac{I+6.25\cdot {V}_{i}-13.1}{\varphi}\right)\right]$ (1)

I represent the seismic hazard described in terms of macroseismic intensity (EMS98 scale), V_{i} the vulnerability index, and
$\varphi $ the ductility index, which is evaluated taking into account the building typology and its constructive features ( [39] ); it controls the slope of the curves and assumes different values to fit the data obtained through damage surveys. For residential buildings, it takes a value of 2.3.

Finally, the damage distribution is derived using beta probability density function (Equation (2)) and the beta cumulative density function (Equation (3)).

${P}_{\beta}\left(x\right)=\frac{\text{\Gamma}\left(t\right)}{\text{\Gamma}(r)\cdot \text{\Gamma}\left(t-r\right)}\cdot \frac{{\left(x-a\right)}^{r-1}\cdot {\left(b-x\right)}^{t-r-1}}{{\left(b-a\right)}^{t-1}}$ (2)

${P}_{\beta}\left(x\right)={\displaystyle {\int}_{x}^{a}{P}_{\beta}\left(\epsilon \right)\cdot d\epsilon}$ (3)

With 𝑎 < 𝑥 < 𝑏 and $r=t\left(0.007{\mu}_{D}^{3}-0.052{\mu}_{D}^{2}+0.287{\mu}_{D}\right)$ . The parameters a, b, t, and r are the parameters of the distribution and Γ is the gamma function. The parameter r is the variance and controls the shape of the distribution.

The parameter t (named after T_{beta}) resulting from the probability calculation by the law β determines the dispersion of the values; it therefore represents the propagation of uncertainties that can come from different sources as the input data, or the variable behaviour of the same typology structure to the seismic hazard. In the Risk-EU method, this one is fixed to t = 8 which represent the best fits for the European buildings [40] . However, the study posits to leave it variable in a range of values, hence the optimal value for each typology can be found numerically.

In order to use the beta distribution, it is necessary to refer to the damage grades D_{k} (k = 0 to 5) defined by EMS scale; for this purpose, it is advisable to assign value 0 to the parameter a and value 6 to the parameter b [18] .

The outcome is the distribution in terms of the six levels defined in EMS98: D0 (undamaged), D1 (slight damage), D2 (moderate damage), D3 (heavy damage), D4 (partial collapse) and D5 (total collapse) for each location that is separately considered, given by the equation (Equation (4)).

$P\left(D\ge {D}_{k}\right)=1-{P}_{\beta}\left(k\right)$ (4)

To perform seismic damage scenarios, very often, many typologies of buildings exist in the same studying area. Therefore, to run the simulation model (Armagedom software) we need to describe:

・ the total number of typologies;

・ the vulnerability indices (V_{i}) and T_{beta} parameters for each typology;

・ the spatial distribution of the typologies in each district (polygon) of the affected area vulnerability (nbBat).

In Armagedom software, these parameters are stored in file (*.txt) as follows:

・ the V_{i} and T_{beta} are regrouped in the same file (Figure 3) called *.tvi_t. The file is structured in four columns;

・ type of building as character, vulnerability index as float, T_{beta} parameters as float and period as float. Each line represents a specific typology;

・ the nbBat are stocked in another file (Figure 3) called *.nbbat. This file is given as table in which, each line represents the area, the number of columns represent the number of typology and the value the number of buildings.

Examples of *.nbbat and *.tvi_t files for the case with three typologies (A1, A2, A3) are illustrated in Figure 3.

3.2. Vulnerability Assessment Using Inverse Optimization

Here, the proposed methodology that answer to the question: “What is (are) the value(s) of vulnerability parameter(s) that optimize the fit of the model’s prediction

Figure 3. Required file for Armagedom software for the case with three typologies (A1, A2, A3): top―the *.tvi_t file containing the vulnerability index V_{i} and bottom―the *.nbbat file containing the distribution of the typologies in each polygon.

damages to the observed ones?” This question is the main issue for all the seismic risk analysis either produced for prediction without a situation of seismic event or produced in quasi-real time situation during a seismic crisis.

The observed data in Table 2 give two major information:

i) the macroseismic EMS98 intensity, which is uniform (IX) for all the city;

ii) the observed damages, regrouped on three classes: green (regroups D0, D1 and D2), orange (represents D3) and red (regroups D4 and D5).

Our problem is clearly considered underdetermined (for number of typologies higher than 1) because there is fewer equations than unknowns. Here, the total number of parameters needed to run the model depends on the total number of typology. Therefore, for 𝑁_{PQRS} the number of parameters to determine is
$\text{2}\times {N}_{\text{PQRS}}+\left({N}_{\text{PQRS}}-\text{1}\right)\times {N}_{\text{V}<\text{WC}}$ . Therefore, for this ill-posed problem, the optimization algorithm leads to local solutions. This limitation, often problematic for solving optimization problem is very interesting for us. Where, not only one final solution is obtained but all solutions minimizing the objective function are retained. The final decision comes to the expert to choose/validate the most appropriate one.

The expert can also reduce the solution space by introducing additional constraints to get more “realistic solutions”. These constraints can be related to the operational needs, better knowledge on the studying area, etc. Follow are examples of questions than can be answered by experts:

・ What error needs to be reduced?

1) Uniform error on all damage degree D0 to D5;

2) Error on D4 to D5 damage, useful for civil protection in case of seismic catastrophe;

3) Error on D0 to D2 damage, useful for mitigation and planning.

・ Which parameters can be fixed?

1) Number of typology and its corresponding V_{i} and T_{beta};

2) The nbBat parameters; repartition of number of buildings per typology in each area.

・ Which parameter can be constrained?

1) Number of typology between 1 and 6;

2) Values of V_{i}, not varying from 0 to 1 but constrained by the fuzzy function (Risk-UE) method (Figure 4);

3) T_{beta} between 4 and 16.

In our study, we have explored a large number of solution in order to evaluate the impact of expert’s constraints on the final solutions. Finally, validate the most appropriate one. We underline that, the range of variation of the vulnerability index is governed by the fuzzy function method (Figure 4), which adapt with the number of typologies [18] [40] [41] .

For our simulations, we used Armagedom for the forward modelling and the augmented Lagrange multiplier method described in [42] for the inverse one. Figure 5 illustrate the conceptual scheme followed for solving the problem.

The proposed methodology was implemented in R programming language and environment [43] using the SOLNP algorithm available in “Rsolnp” package [44] . This process was repeated 100 times to generate 100 values of V_{i}, T_{beta} and nbBat values. These solutions represent potential solution, which will be analyzed by the expert in order to select the most appropriate one. We will show in the next section, that all the 100 solutions present a very little difference. The algorithmic parameters used are as follows: population size = 50, generations = 200, crossover fraction = 0.8, and stall generations = 100. In fitting convex function

Figure 4. Fuzzy function for vulnerability classes (from A to F following EMS98) that constrain the vulnerability index for the optimization method.

Figure 5. Conceptual scheme for optimization process.

models for the data, we use ellipsoidal, super-ellipsoidal, and Nataf transformation based methods.

For each observed damage grade, we compute the misfit error between observed and estimated damage (Equation (5)). In order to have a uniform distributed damage grade error for all polygons we define for each damage grade the entity defined in Equation (6). Finally, the objective function to minimize can be defined as a weighted sum of the theses entities (Equation (7)).

$\{\begin{array}{l}{E}_{{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}}=\left|{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}{}^{Est}-{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}{}^{Obs}\right|/{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}{}^{Obs}\\ {E}_{{\text{D}}_{3}}=\left|{\text{D}}_{3}{}^{Est}-{\text{D}}_{3}{}^{Obs}\right|/{\text{D}}_{3}{}^{Obs}\\ {E}_{{\text{D}}_{4}{\text{D}}_{5}}=\left|{\text{D}}_{4}{\text{D}}_{5}{}^{Est}-{\text{D}}_{4}{\text{D}}_{5}{}^{Obs}\right|/{\text{D}}_{4}{\text{D}}_{5}{}^{Obs}\end{array}$ (5)

$\{\begin{array}{l}{J}_{1}=\mathrm{max}\left({E}_{{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}}\right)-\mathrm{min}\left({E}_{{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}}\right)+{\stackrel{^}{E}}_{{\text{D}}_{0}{\text{D}}_{1}{\text{D}}_{2}}\\ {J}_{2}=\mathrm{max}\left({E}_{{\text{D}}_{3}}\right)-\mathrm{min}\left({E}_{{\text{D}}_{3}}\right)+{\stackrel{^}{E}}_{{\text{D}}_{3}}\\ {J}_{3}=\mathrm{max}\left({E}_{{\text{D}}_{4}{\text{D}}_{5}}\right)-\mathrm{min}\left({E}_{{\text{D}}_{4}{\text{D}}_{5}}\right)+{\stackrel{^}{E}}_{{\text{D}}_{4}{\text{D}}_{5}}\end{array}$ (6)

$J={\omega}_{1}\cdot {J}_{1}+{\omega}_{2}\cdot {J}_{21}+{\omega}_{3}\cdot {J}_{3}$ (7)

The choice of the weights
${\omega}_{1},{\omega}_{2},{\omega}_{3}$ depends on the application and should be fixed by the expert. For example, for people who make organizational decisions J_{1} must be important for them, therefore
${\omega}_{1}$ should be greater than
${\omega}_{2}$ _{ }and
${\omega}_{3}$ . In the other hand for civil protection and urgent action after a seismic event,
${\omega}_{3}$ should be greater than
${\omega}_{2}$ and
${\omega}_{1}$ .

In this study, we have tested 4 hypotheses:

1) Case 1: error on D4 to D5 useful for civil protection in case of seismic catastrophe; in this case, we take $\left({\omega}_{1},{\omega}_{2},{\omega}_{3}\right)=\left(0.0\text{1},0.0\text{1},\text{1}\right).$

2) Case 2: error on the damage degree D0 to D2, useful for mitigation and planning; in this case, we take $\left({\omega}_{1},{\omega}_{2},{\omega}_{3}\right)=\left(\text{1},0.0\text{1},0.0\text{1}\right).$

3) Case 3: error on damage degree D3 the most complicated to evaluate; in this case, we take $\left({\omega}_{1},{\omega}_{2},{\omega}_{3}\right)=\left(0.0\text{1},\text{1},0.0\text{1}\right).$

4) Case 4: uniform error on all damage degree; in this case, we take $\left({\omega}_{1},{\omega}_{2},{\omega}_{3}\right)=\left(\text{1},\text{1},\text{1}\right).$

Figure 6 illustrates the organigram of explored solutions.

The number of simulations runs needed was determined in order that providing stable predictions for numerical results. The choice of 100 was guided by two raisons:

-The variability of 100 estimated V_{i} indices values is very low. This is shown in section 4 (Figure 8(a1) and Figure 8(a2)).

-We first execute 10, 20, 50 then 100 model runs during the experimentation phase. We clearly observe that first (mean) and second (variance) moment-order for 50 and 100 models run were similar. As optimization process is not time-consuming process, we fix to number of model runs to 100.

4. Results and Comments

In this section, we present a synthetic analysis of obtained results. We remember that our main objective here, is the determination of optimal vulnerabilities parameters (V_{i}, T_{beta} and nbBat) fitting at best the historical observed damages. The originality of the method relies on practical and operational aspects. Indeed, the method gives a panel of optimal solutions (ill-posed problem), and final choice return experts. For clarity of exposition, we will first present the typical results obtained for a specific hypothesis (here case 4; uniform distributed error for all damages) and for a fixed number of typologies (here 4 typologies). Then, synthetic results for two cases (from D0 to D2 in Figure 9 and, from D4 to D5 in Figure 8) of all executed simulations (2400 for each case).

Figure 6. Number of runs for all the study by each step.

Figure 7 represents an example of obtained results for a number of typologies fixed to 4 for the case 4. The optimization process provides 100 possible solutions to experts for the vulnerability index (Figure 7(a)), the T_{beta} parameters (Figure 7(b)) and nbBat to choose/validate the most appropriate one. In Figure 7(a) (respectively Figure 7(b)), different colors (blue, green, red and black color) of the points represent vulnerability index (respectively T_{beta}) for the 4 typologies. Its show that:

・ V_{i}’s values are very close even the variation range are very high (fuzzy function in Figure 4). For example, blue points (first typology) in Figure 7(a) gives values of V_{i} around 0.8 and can take values between 0.785 and 0.965.

・ Even if, the values of V_{i} present a small variability, the value of T_{beta} parameter (Figure 7(b)) have a high dispersion (between 4 and 16).

Figure 7(c) and Figure 7(d) present the errors on the number of structures given in the number and in percentage for two arbitrary solutions (samples 17 and 61). The barplots presents the error on the damages degree D0 to D2 (green color), D3 (orange color) and D4 to D5 (red color) for the 10 districts. It shows that the distribution is very different for the two case, in the other hand the error

Figure 7. Mathematical solutions for the vulnerability parameters (V_{i} and T_{beta}) for case 1
$\left({\omega}_{1}=1,{\omega}_{2}=1,{\omega}_{3}=1\right)$ with 4 typologies: (a) 100 possible solutions for V_{i}; (b) 100 possible solutions for T_{beta}; (c) Damages errors on the number of buildings represented in number and percentage for sample 17; (d) Damages errors on the number of buildings represented in number and percentage for sample 61. Each typology is represented by one color (blue, green, red and black).

remains distributed uniformly. This can be explained by the fact that the problem is ill posed and infinity of local minima are possible, and very different solution can satisfy the same optimization criteria. However, this limitation is very interesting in our study, were experts must validate the optimal solution. These results emphasise the idea of global consideration of results, by looking to the error at the scale of the city and in all the districts of the studied area. Looking separately at a zoomed area the interpretation of results can be inexact, since, for example for the solution 17, for the district five the error is very reduced, and vice versa, for example for the district one of the solution 61.

In follow, we will present the obtained results for two practical cases.

・ The first one, very useful for crisis management. Where, the authority is interested by the number of highly damaged buildings; here named as D4 to D5. All other categories, remains less important in this case.

・ The second one, very useful for mitigation and planning. Where, concerned authority need information about habitability.

For each case presented below, we have executed 2400 simulations in order to test some hypothesis and it allows to answer to two important questions:

1) Does higher number of typologies implies lower misfit error?

2) In the Risk-UE method, T_{beta} parameter was calibrated according to the European buildings and fixed to 8. So, changing its value improves the accuracy of the model predictions?

4.1. Case1: Crisis Management

In this section, we will summarize the obtained results for the 2400 executed simulations. In Figure 8, the graphics in the left column show results for fixed T_{beta} parameters (equal to 8) according to the European buildings, and the right one for optimal T_{beta}. We present the results for the three sets of observational damages grade described below (D0 to D2, D3, D4 to D5).

In Figure 8(a1) and Figure 8(a2), boxplots represent results for possible value for vulnerability indices (y-axis) for each number of typologies (x-axis). Each boxplot contains 100 values. We observe that the variability of V_{i} indices values is very low for the two cases; fixed and optimal T_{beta}. In Figure 8(a1) for one typology, the V_{i} value is always the same because the optimization process converges for the global minimum. In this specific case, the well-posed problem has a unique solution, which is V_{i} equal to 0.70; T_{beta} is fixed to 8 and nbBat is equal to total number of buildings.

In Figure 8(b1) and Figure 8(b2), boxplots represent results for possible value for T_{beta} parameters (y-axis) for each number of typologies (x-axis). Each boxplot contains 100 values. In Figure 8(b1), all T_{beta} values are equal to 8 because it is fixed input parameter; the graphic plotted only to show for comparison with Figure 8(b2). This last, show clearly that for all number of typologies, the T_{beta} parameters varies from 4 to 16 (expert’s constraint) with a high dispersion.

In Figure 8(c1) and Figure 8(c2), boxplots represent results for mean

Figure 8. Results for crisis management case: (a1) the value of V_{i} indices for different typologies with fixed T_{beta}, (a2) the value of V_{i} indices for different typologies with optimal T_{beta}, (b1) the value of T_{beta} parameters for different typologies with fixed T_{beta}, (b2) the value of T_{beta} parameters for different typologies with optimal T_{beta}, (c1) the mean absolute error in percentage for different typologies with fixed T_{beta}, (c2) the mean absolute error in percentage for different typologies with optimal T_{beta}, (d) Comparison between minimal mean absolute error in percentage for fixed and optimal T_{beta}.

absolute error (y-axis) for each number of typologies (x-axis). Each boxplot contains 100 values. Here, we have compare the results for three indicators: mean absolute error on D0 to D2 (green), on D3 (orange) and on D4 to D5 (red). We remember that here, we are minimizing according to high weight on D4 to D5. Therefore, we expect a best performance on D4 to D5 (red) and less on the other one. The two graphics Figure 8(c1) and Figure 8(c2) confirm that error on D4 to D5 is very small (less than 10%) for the 2 cases; with fixed and optimal T_{beta} comparing to the error on D3 and D4 to D5.

Figure 8(d) reports the lowers band in Figure 8(c1) (continues lines) et Figure 8(c2) (dashed lines) in the same graphics. Each plot represents the minimum error (among the 100) for the considered indicator for each number of typologies. Bottom and top of the boxplot in Figure 8 are the 25^{th} and 75^{th} percentile, the ban near the middle of the box is the median and the ends of the whiskers are the minimum and maximum. This graphics allows the comparison between the two cases (fixed and optimal T_{beta}).

We observe that for:

・ D4 to D5 indicator that for fixed and optimal T_{beta} the misfit error on D4 to D5 is less than 2%, except for the case with one typology. Moreover, the increase of the number of typologies does not improve the model’s results. This can be explained by the fact that increasing a number of typologies, the V_{i} index have more constraints for variation range given by the fuzzy function. And the second one that T_{beta} fixed to 8 was accurately calibrated in the Risk-UE project [15] and it is representative of European buildings.

・ D0 to D2 and D3 indicators, we show that considering T_{beta} variable and increasing the number of typologies reduce the misfit error (up to 20%). However, it’s meaningless because, the optimization criterion is highly weighted on D4 to D5.

For El Asnam area, for the crisis management situation, after the analysis of the mathematical solutions, we suggest to use for the predictive damage scenario five or six typologies, with the V_{i} equal to the values presented in the Figure 4 and the T_{beta} fix equal to 8. The variation of the T_{beta} do not ameliorate the results and the use of less typologies give less information for the same error (less than 2%). The choice of the expert is constantly balanced between the acceptable error and the needed information. In the case of El Asnam, if we are looking only to the D4 to D5 damages, as is generally the case for crisis management, we prefer to choose more typologies (that will represent better the geographical distribution of the damage) for the same error. However, if for one reason, we are also interested by D3 damages, in this case the error increase at around 30% for five typologies and hence, we recommend using two typologies for which the error is less than 10%.

4.2. Case 2: Mitigation and Planning

Figure 9 summarize results of the obtained 2400 simulations. It’s organized in

Figure 9. Results for mitigation and planning case: (a1) the value of V_{i} indices for different typologies with fixed T_{beta}, (a2) the value of V_{i }indices for different typologies with optimal T_{beta}, (b1) the value of T_{beta} parameters for different typologies with fixed T_{beta}, (b2) the value of T_{beta} parameters for different typologies with optimal T_{beta}, (c1) the mean absolute error in percentage for different typologies with fixed T_{beta}, (c2) the mean absolute error in percentage for different typologies with optimal T_{beta} and (d) comparison between minimal mean absolute error in percentage for fixed and optimal T_{beta}.

same manner as Figure 8. In Figure 9(a1) and Figure 9(a2), boxplots represent results for possible value for vulnerability index (y-axis) for each number of typologies (x-axis). Each boxplot contains 100 values. We observe that minor differences in the V_{i} indices value variability exist for the two cases; fixed and optimal T_{beta}. In Figure 9(a1) for 1, 3 or 5 typologies, the V_{i} value always converge to a global minimum. This is very specific to study case and cannot be generalized. The same analysis can be reported in Figure 9(a2) for 5 typologies.

In Figure 9(b1) and Figure 9(b2), boxplots represent results for possible value for T_{beta} parameters (y-axis) for each number of typologies (x-axis). The same analysis can be done here as in Figure 9(b1) et Figure 9(b2); the values of T_{beta} parameters are highly dispersive.

In Figure 9(c1) and Figure 9(c2), boxplots represent results for mean absolute error (y-axis) for each number of typologies (x-axis). Each boxplot contains 100 values. Here, we have compare the results for three indicators: mean absolute error on D0 to D2 (green), on D3 (orange) and on D4 to D5 (red). We remember that we are minimizing with high weight on D0 to D2. Therefore, we expect a best performance on D0 to D2 (green) and less on the other one. The two graphics seems giving equivalent results in order of magnitude for D0 to D2 errors. In addition, we observe that this error grow with the increase of the number of typologies.

Figure 9(d) reports the lowers band in Figure 9(c1) (continues lines) and Figure 9(c2) (dashed lines). Each plot represents the minimum error (among the 100) for the considered indicator for each number of typologies. This graphics allows the comparison between the two cases (fixed and optimal T_{beta}). We observe that:

・ For D0 to D2 indicator that for fixed and optimal T_{beta} the misfit error on D0 to D2 is less than 20%. We observe that the increase of the number of typologies increase the misfit error. This can be explained by the fact that increasing a number of typologies, the V_{i}_{ }index have more constraints for variation range given by the fuzzy function. And the second one that T_{beta} fixed to 8 was accurately calibrated in the Risk-UE project (Mouroux et al., 2004) and it is representative of European buildings.

・ For D0 to D2 and D3 indicators, we show that considering T_{beta} variable and increasing the number of typologies reduce the misfit error (up to 2%) for D4 to D5 indicator. We also notice, that, for 2 typologies the result for D4 to D5 indicators with fixed T_{beta} is better than with optimal one. These results depend on the study case and it’s meaningless to extrapolate the results except for D0 to D2 indicator in this case.

For El Asnam area, for the mitigation and planning situation, after the analysis of the mathematical solutions, we suggest to use for the predictive damage scenario two typologies, with the V_{i} equal to (0.6 and 0.87) and the T_{beta} fix equal to 8. If for one reason, the used decide to use a different number of typologies, in this case we recommend using variable values of T_{beta}, since in this case this ameliorate the results (from Figure 9(d). we can see that the error is quite stable for 3, 4 and 5 typologies with optimal T_{beta}, but increase for the same number of typologies with T_{beta} fixed to 8).

In Figure 9, the bottom and top of the boxplot are the 25^{th} and 75^{th} percentile, the ban near the middle of the box is the median and the ends of the whiskers are the minimum and maximum. From the results of the two cases, we can deduce that the methodology is calibrated for table the important damage degrees, that can be managed by the optimization process as can be noticed in the Figure 8(d), where the error is less than 2% for any number of typology. This is more hardly to manage for the error relative to the case 2: mitigation and planning, for which, the error, although stable for 1, 2 and 3 typologies, grow up for the situations with 4, 5 and 6 typologies.

The fact that we use different typologies introduces more error, but that also gives supplementary information relative to the location of the damages. If the question of the user is just to know the total error without knowing the distribution by districts, in this case the use of only one typology can be useful. However, generally, even in crises management situation, in addition to this question, the second question is the identification of the most affected districts in order to send the rescue. In this case, the geographical description of the error is required.

5. Conclusions

We have aimed in this work to develop unified computational method for assessment of vulnerability of structures, to be used for seismic risk scenarios at the town scale, when the underlying uncertainties are modelled using random variables, interval analysis, and (or) fuzzy variables. The developed approach is based on the minimization of the error between observed and estimated damages and it allows the determination of the number of typologies classes, the estimation of their vulnerability index and associated T_{beta} value as well as the spatial distribution in the studied area (nbBat). This method remains general and can be applied for any observed dataset, for which generally the information related to the typologies and their vulnerability are not available. The main insights after the analysis of the results are:

・ The Risk-EU methodology and the model that it proposed fit very well the damages D4 to D5. When the attention is played to the damages D0 to D2 the estimation is correct but damage D3 remains the most difficult to be assessing because the diagnosis experts on this category remain very subjective and the description remains very vague.

・ The results of our work corroborate the affirmation on the Risk-EU methodology, where the T_{beta} was fixed equal to 8, based on the statistic treatments of all the European data. Our analysis (Figure 9(b2) and Figure 8(b2)) shows that by modifying it between 4 and 16, the best results are improved by 10%. Thus, our approach confirms that the latter does not bring too much variability on the results and the choice of T_{beta} = 8 must be respected.

・ Contrary to the expectance, increasing the number of typologies does not reduce the error for the results. Indeed, the fact of constraining by the expert the V_{i} by the fuzzy function restricts the search space. Therefore, the total error is increased but, by increasing the number of typologies, the geographic description of the vulnerability in the studied area is improved and brings supplementary needed information (e.g. the hierarchizing of the most damaged areas).

The two cases (crisis management and planning and mitigation) have two different trends response: i) for mitigation case, increasing the number of typology increases the misfit error; ii) on the other hand in the case of crisis management increasing the number of typology improves the model prediction. In some cases, the optimization algorithm can detect the number of classes that minimize the error, that means on one hand that, a smaller number of classes are insufficient (small number of degree of freedom), and on the other hand that, spending time for dividing and refining the number of classes is useless. In the case of El Asnam situation, the optimum number of typologies was not clearly identified. The mathematical formulation of the objective function gives the opportunity to have a perfect fit of the vulnerability parameters (error that tends to zero) by minimizing the error for a specific damage grade (e.g. a minimization related to D4 to D5 damage gives the most accurate vulnerability parameters for civil protection for urgent action after a seismic event). Hence, the approach can answer to different actors that work at the scale of the town in different situations (e.g. planning, mitigation action, retrofitting, crises management). Of course, the total error (total number of building as a sum of building in each area) is much less big than the error related to each area. This observation should be considered in the studies in which only the total error is computed, which have certainly big utilities in certain almost “real time crises situation” when the first question is the rough estimate of number of victims; but quickly after, the spatial localization of damages is asked, in order to know where to send the first aids, and hence the error at the area level become very important. The method that we propose presents a better integration of the vulnerability and loss results could allow city councils or regional authorities to plan interventions based on a global view of the site under analysis, leading to more accurate and comprehensive risk mitigation strategies that support the requirements of safety and emergency planning.

Cite this paper

Benaïchouche, A. , Negulescu, C. , Sedan, O. and Boutaraa, Z. (2018) Seismic Vulnerability Assessment from Earthquake Damages Historical Data Using Constrained Optimization Technique.*Journal of Geoscience and Environment Protection*, **6**, 89-111. doi: 10.4236/gep.2018.62007.

Benaïchouche, A. , Negulescu, C. , Sedan, O. and Boutaraa, Z. (2018) Seismic Vulnerability Assessment from Earthquake Damages Historical Data Using Constrained Optimization Technique.

References

[1] Adeli, H. and Ge, Y. (1989) A Dynamic Programming Method for Analysis of Bridges under Multiple Moving Loads. International Journal for Numerical Methods in Engineering, 28, 1265-1282.

https://doi.org/10.1002/nme.1620280604

[2] Adeli, H. and Mak, K.Y. (1990) Interactive Optimization of Plate Girder Bridges Subjected to Moving Loads. Computer-Aided Design, 22, 368-376.

https://doi.org/10.1016/0010-4485(90)90087-S

[3] Sirca Jr., G.F. and Adeli, H. (2005) Cost Optimization of Prestressed Concrete Bridges. Journal of Structural Engineering, 131, 380-388.

https://doi.org/10.1061/(ASCE)0733-9445(2005)131:3(380)

[4] Adeli, H. and Balasubramanyam, K.V. (1988) Expert Systems for Structural Design: A New Generation. Prentice Hall, Upper Saddle River.

[5] Bisadi, V. and Padgett, J.E. (2015) Explicit Time-Dependent Multi-Hazard Cost Analysis Based on Parameterized Demand Models for the Optimum Design of Bridge Structures. Computer-Aided Civil and Infrastructure Engineering, 30, 541-554.

https://doi.org/10.1111/mice.12131

[6] Sun, H. and Betti, R. (2015) A Hybrid Optimization Algorithm with Bayesian Inference for Probabilistic Model Updating. Computer-Aided Civil and Infrastructure Engineering, 30, 602-619.

https://doi.org/10.1111/mice.12142

[7] Forcael, E., González, V., Orozco, F., Vargas, S., Pantoja, A. and Moscoso, P. (2014) Ant Colony Optimization Model for Tsunamis Evacuation Routes. Computer-Aided Civil and Infrastructure Engineering, 29, 723-737.

https://doi.org/10.1111/mice.12113

[8] Faturechi, R. and Miller-Hooks, E. (2014) A Mathematical Framework for Quantifying and Optimizing Protective Actions for Civil Infrastructure Systems. Computer-Aided Civil and Infrastructure Engineering, 29, 572-589.

[9] Opricovic, S. and Tzeng, G.H. (2002) Multicriteria Planning of Post-Earthquake Sustainable Reconstruction. Computer-Aided Civil and Infrastructure Engineering, 17, 211-220.

https://doi.org/10.1111/1467-8667.00269

[10] Strasser, F.O., Bommer, J.J., Sesetyan, K., Erdik, M., Cagnan, Z., Irizarry, J., Crowley, H., et al. (2008) A Comparative Study of European Earthquake Loss Estimation Tools for a Scenario in Istanbul. Journal of Earthquake Engineering, 12, 246-256.

https://doi.org/10.1080/13632460802014188

[11] Molina, S., Lang, D.H. and Lindholm, C.D. (2010) SELENA—An Open-Source Tool for Seismic Risk and Loss Assessment using a Logic Tree Computation Procedure. Computers & Geosciences, 36, 257-269.

https://doi.org/10.1016/j.cageo.2009.07.006

[12] Hancilar, U., Tuzun, C., Yenidogan, C. and Erdik, M. (2010) ELER Software—A New Tool for Urban Earthquake Loss Assessment. Natural Hazards and Earth System Sciences, 10, 2677-2696.

https://doi.org/10.5194/nhess-10-2677-2010

[13] Sedan, O., Negulescu, C., Terrier, M., Roullé, A., Winter, T. and Bertil, D. (2013) Armagedom—A Tool for Seismic Risk Assessment Illustrated with Applications. Journal of Earthquake Engineering, 17, 253-281.

https://doi.org/10.1080/13632469.2012.726604

[14] Calvi, G.M., Pinho, R., Magenes, G., Bommer, J.J., Restrepo-Vélez, L.F. and Crowley, H. (2006) Development of Seismic Vulnerability Assessment Methodologies over the Past 30 Years. ISET Journal of Earthquake Technology, 43, 75-104.

[15] Mouroux, P., Bertrand, E., Bour, M., Le Brun, B., Depinois, S. and Masure, P. (2004) The European RISK-UE Project: An Advanced Approach to Earthquake Risk Scenarios. The 13th World Conference on Earthquake Engineering, Vancouver, 1-6 August 2004, Paper No. 3329.

[16] Petrini, V. (1993) Rischio sismico di edifici pubblici, Parte I Aspetti metodologici. GNDT-CNR Report, Tipografia Moderna.

[17] HAZUS99, F.E.M.A. (1999) Earthquake Loss Estimation Methodology: User’s Manual. Federal Emergency Management Agency, Washington DC.

[18] Giovinazzi, S. (2005) The Vulnerability Assessment and the Damage Scenario in Seismic Risk Analysis. PhD Thesis.

[19] Borzi, B., Pinho, R. and Crowley, H. (2008) Simplified Pushover-Based Vulnerability Analysis for Large-Scale Assessment of RC Buildings. Engineering Structures, 30, 804-820.

https://doi.org/10.1016/j.engstruct.2007.05.021

[20] Crowley, H., Pinho, R. and Bommer, J.J. (2004) A Probabilistic Displacement-Based Vulnerability Assessment Procedure for Earthquake Loss Estimation. Bulletin of Earthquake Engineering, 2, 173-219.

https://doi.org/10.1007/s10518-004-2290-8

[21] Orsini, G. (1999) A Model for Buildings’ Vulnerability Assessment using the Parameterless Scale of Seismic Intensity (PSI). Earthquake Spectra, 15, 463-483.

https://doi.org/10.1193/1.1586053

[22] Crowley, H., Colombi, M., Crempien, J., Erduran, E., Lopez, M., Liu, H., Milanesi, M., et al. (2010) GEM1 Seismic Risk Report: Part 1. Technical Report, GEM Foundation, Pavia, 5.

[23] Crowley, H., Colombi, M., Crempien, J., Erduran, E., Lopez, M., Liu, H., Milanesi, M., et al. (2010) GEM1 Seismic Risk Report: Part 2. Technical Report, GEM Foundation, Pavia, 5.

[24] Whitman, R.V., Reed, J.W. and Hong, S.T. (1973) Earthquake Damage Probability Matrices. In: Proceedings of the 5th World Conference on Earthquake Engineering, Vol. 2, Palazzo DEI CONGRESSI (EUR), Rome, 2531-2540.

[25] Grunthal, G. (1998) European Macroseismic Scale 1998. Cahiers du centre Européen de Géodynamique et de Séismologie, 15, 1-99.

[26] Sieberg, A. (1930) Geologie der Erdbeben. Handbuch der Geophysik, 2, 552-555.

[27] Wood, H.O. and Neumann, F. (1931) Modified Mercalli Intensity Scale of 1931. Bulletin of the Seismological Society of America, 21, 277-283.

[28] Medvedev, S.V., Spohneuer, W. and Karnik, V. (1965) Seismic Intensity Scale Version MSK 1964. Akad. Nauk SSSR, Geofiz. Kom., 10. (In Moscow)

[29] Colombi, M., Borzi, B., Crowley, H., Onida, M., Meroni, F. and Pinho, R. (2008) Deriving Vulnerability Curves Using Italian Earthquake Damage Data. Bulletin of Earthquake Engineering, 6, 485-504.

https://doi.org/10.1007/s10518-008-9073-6

[30] Lambert, J., Levret-Albaret, A., Cushing, M. and Durouchoux, C. (1996) Mille ans de séismes en France. Ouest Editions, Nantes, 79.

[31] SisFrance (2005) Catalogue des séismes francais métropolitains, BRGM, EDF, IRSN.

http://www.sisfrance.net

[32] Rossetto, T., Ioannou, I. and Grant, D.N. (2013) Existing Empirical Fragility and Vulnerability Functions: Compendium and Guide for Selection. GEM Technical Report, GEM Foundation, Pavia, 65.

[33] Rossetto, T., Ioannou, I., Grant, D.N. and Maqsood, T. (2014) Guidelines for the Empirical Vulnerability Assessment. GEM Technical Report, GEM Foundation, Pavia, 65.

[34] Boutaraa, Z., Negulescu, C., Arab, A. and Sedan, O. (2018) Buildings Vulnerability Assessment and Damage Seismic Scenarios at Urban Scale: Application to Chlef City (Algeria). KSCE Journal of Civil Engineering.

[35] Boutaraa, Z., Negulescu, C., Arab, A. and Sedan, O. (2015) Scénarios de risque sismique pour la ville de Chlef (Ex El Asnam) Algérie. 9eme Colloque AFPS, Marne la vallée. (In French)

[36] Talaganov, K., Aleksovski, D., Milutinovic, Z., Ameur, B., Arsovski, M., Jancevski, J. and Andreevski, V. (1982) Studies for Elaboration of the Code for Repair and Strengthening of Damaged Buildings in the Region of El Asnam: Engineering Geology, Geotechnical and Geophysical Characteristics of the Town of El Asnam and Other Sites. Report IZIIS, 82-55.

[37] Bertero, V.V. and Shah, H.C. (1983) El-Asnam Algeria Earthquake Oktober 10, 1980: A Reconnaissance and Engineering Report. EERI.

[38] EERI, Database of Concrete Buildings Damaged in Earthquakes.

https://www.eeri.org/2013/11/onlinedatabase-of-concrete-buildings-damaged-in-earthquakes/

[39] Lagomarsino, S. and Giovinazzi, S. (2006) Macroseismic and Mechanical Models for the Vulnerability and Damage Assessment of Current Buildings. Bulletin of Earthquake Engineering, 4, 415-443.

https://doi.org/10.1007/s10518-006-9024-z

[40] Milutinovic, Z.V. and Trendafiloski, G.S. (2003) Risk-UE An Advanced Approach to Earthquake Risk Scenarios with Applications to Different European Towns. Contract: EVK4-CT-2000-00014, WP4: Vulnerability of Current Buildings.

[41] Demartinos, K. and Dritsos, S. (2006) First-Level Pre-Earthquake Assessment of Buildings using Fuzzy Logic. Earthquake Spectra, 22, 865-885.

https://doi.org/10.1193/1.2358176

[42] Ye, Y. (1987) Interior Algorithms for Linear, Quadratic, and Linearly Constrained Non-Linear Programming. PhD Thesis, Department of ESS, Stanford University, Stanford.

[43] Team, R.C. (2015) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna.

[44] Ghalanos, A. and Theussl, S. (2012) Rsolnp: General Non-Linear Optimization Using Augmented Lagrange Multiplier Method. R Package Version, 1.