University of Southern California Dissertations and Theses

Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms

(USC Thesis Other)

Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Novel Multi-Stage and CTLS-Based Model Updating Methods and Real-Time Neural Network-Based Semiactive Model Predictive Control Algorithms by Tianhao Yu A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Civil Engineering) May 2022 Copyright 2022 Tianhao Yu Acknowledgments The author gratefully acknowledges the financial support of the Provost’s fellowship, the Gammel Scholarship and the assistantship appointments from the University of Southern California (USC) as well as the National Science Foundation (NSF) through grants: 13-44937/13-44622 and 14- 36018; the computing resources from the USC’s Center for Advanced Research Computing (CARC); the assistance of Dr. Wael M. Elhaddad (former group member); and the partnerships of E-Defense senior researchers Dr. Eiji Sato and Dr. Tomohiro Sasaki (now at Obayashi R&D). I would like to thank my advisor Professor Erik A. Johnson for his guidance, patience and trust, without which I could never make this far. I also would like to express my gratitude to my doctoral committee members for spending their valuable time to help me go through this important step in my study, career and life, and for providing helpful and invaluable advice. Finally, I would like to thank my family for their continued and unwavering support so that I can focus on my course work and research. ii Table of Contents Acknowledgments ii List of Tables vi List of Figures viii Abstract xi Chapter 1: Introduction 1 1.1 Introduction to Modeling and Model Updating of a Full-Scale Experimental Base- Isolated Building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Introduction to Cross-Model Cross-Mode Model Updating Method Based on Con- strained Total Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Introduction to Real-Time Neural Network based Semiactive Model Predictive Control of Structural Vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Chapter 2: Modeling and Model Updating of a Full-Scale Experimental Base-Isolated Building 17 2.1 Experimental Set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 System Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.1 Test 010 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2 Test 013 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3 Validation of Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 Introduction to the Finite Element Model . . . . . . . . . . . . . . . . . . . . . . 27 2.5 Linear Updating and Mode Matching . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.5.1 Mode Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.5.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6 Results of Linear Model Updating . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.7 Formulation of the Partially Nonlinear FEM . . . . . . . . . . . . . . . . . . . . . 43 2.7.1 Modeling the Isolation-Layer Devices . . . . . . . . . . . . . . . . . . . . 44 2.7.2 Time-History Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Chapter 3: Cross-Model Cross-Mode Model Updating Method Based on Constrained Total Least Squares 52 3.1 Brief introduction to the CMCM Method . . . . . . . . . . . . . . . . . . . . . . . 52 iii 3.1.1 CMCM Method Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.1.2 Model Expansion for Spatially Incomplete Experimental Mode Shapes . . 54 3.2 Comparison of TLS, CTLS and WTLS Approaches . . . . . . . . . . . . . . . . . 57 3.2.1 TLS Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2.2 WTLS Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2.3 CTLS Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2.4 Line-Fitting Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.2.5 Algorithm for Model Updating . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 Numerical Example: 3D Frame Structure . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Introduction to FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.2 Examples and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.3.2.1 Comparison of the TLS and CTLS Approaches . . . . . . . . . . 67 3.3.2.2 Performance Evaluation of the CTLS Approach . . . . . . . . . 68 3.4 Experimental Study: Base-isolated Reinforced-Concrete Frame Building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.4.1 Model Updating Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.4.2 Model Updating Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.4.3 Time Histories from Experimental and Updated FEM . . . . . . . . . . . . 75 Chapter 4: Real-Time Neural Network Based Semiactive Model Predictive Control of Structural Vibrations 82 4.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Offline neural network training . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.2.1 Shrinking input parameter space — homogeneity . . . . . . . . . . . . . . 86 4.2.2 Sampling in the input parameter space . . . . . . . . . . . . . . . . . . . . 87 4.2.3 Breaking a large NN into smaller NNs . . . . . . . . . . . . . . . . . . . . 88 4.3 Real-time implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.3.1 Adaptingκ at each time step . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.3.2 Bounding the sub-optimality . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.3.3 Warm-start technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3.4 Formulating the fMPC problem at each time step . . . . . . . . . . . . . . 94 4.3.5 Steps for real-time implementation of the proposed algorithm . . . . . . . 95 4.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.4.1 Single degree of freedom system . . . . . . . . . . . . . . . . . . . . . . . 99 4.4.2 Three degree of freedom system . . . . . . . . . . . . . . . . . . . . . . . 101 4.4.3 Five degree of freedom system with one controllable damper . . . . . . . . 105 4.4.3.1 Control design and performance evaluation . . . . . . . . . . . . 105 4.4.3.2 Robustness of Control Performance . . . . . . . . . . . . . . . . 108 4.4.4 Five degree of freedom system with four controllable dampers . . . . . . . 112 Chapter 5: Conclusions and Directions for Future Research 119 5.1 Conclusions — Modeling and Model Updating of a Full-Scale Experimental Base- Isolated Building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 Conclusions — Cross-Model Cross-Mode Model Updating Method Based on Con- strained Total Least Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 iv 5.3 Conclusions — Real-Time Neural Network based Semiactive Model Predictive Control of Structural Vibrations . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.4 Directions for Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Bibliography 126 Appendices 134 A Two-Level Strategy for Optimizing Isolation-Layer Device Parameters . . . . . . . 134 B Data and Code Availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 C Derivation for Eq. (3.16) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 D Derivation for Eq. (3.22) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 E Strategy Switch When Relative Velocity Switch Sign . . . . . . . . . . . . . . . . 140 F Matrices and Vectors in Eq. (4.10) . . . . . . . . . . . . . . . . . . . . . . . . . . 141 v List of Tables 2.1 Descriptions of excitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Identified natural frequencies and damping ratios (Test 010) (adapted from Brewick et al., 2018) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3 Identified effective natural frequencies and damping ratios (Test 013) . . . . . . . . 23 2.4 Baseline 3DOF model natural frequencies and damping ratios, and the natural fre- quencies of 3DOF models with 20% softer isolation or 20% softer superstructure (∆ denotes the frequency change relative to the baseline; the greater-magnitude light gray entries are those expected to change and are not relevant in this analysis) 25 2.5 Nominal compressive strength f ′ ck and Young’s modulus E c , as well as the corre- sponding post-update Young’s modulusE c , for each superstructure concrete material 29 2.6 Relationship between characteristic compressive strength f ′ ck and Young’s modulus E c in the Japanese concrete design code (Ueda, 2007) . . . . . . . . . . . . . . . . 30 2.7 Nominal and updated stiffnesses and damping coefficients of isolation-layer devices 32 2.8 Frequency comparisons of the experimental building, and the original and updated FEMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.9 Effective mass and moment of inertia ratios of the FEM modes matched to identi- fied modes 1–8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.10 Percentage errors, using metric (2.8), in acceleration response and isolation-layer device force RMSs using the linear FEM before (after) updating . . . . . . . . . . 42 2.11 Parameters in the RB, ESB and SDP models . . . . . . . . . . . . . . . . . . . . . 47 2.12 Percentage errors, using metric (2.8), in acceleration response and isolation-layer device force RMSs using the partially nonlinear system model . . . . . . . . . . . 49 3.1 Slope and intercept statistics for fitting a line to 100 points . . . . . . . . . . . . . 62 3.2 Comparing modal characteristics of the unupdated and “damaged” FEMs . . . . . 67 vi 3.3 Average and standard deviation of magnitude change ˆ x for comparison of the TLS and the CTLS approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.4 Average and standard deviation of magnitude change ˆ x estimated by the CTLS approach for different noise levels . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.5 Nominal and updated stiffnesses and damping coefficients of isolation-layer devices 74 3.6 Frequency comparisons of the experimental building, and the unupdated and up- dated FEMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.7 Effective mass ratio of FEM modes 1–7 and 9 . . . . . . . . . . . . . . . . . . . . 77 3.8 Percentage errors in RMS acceleration responses and isolation-layer device forces before (after) updating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.1 Computation time (in milliseconds) for running one step of the hMPC and the NN+fMPC algorithms (one-controllable-damper examples) . . . . . . . . . . . . . 94 4.2 Information for Earthquake Records . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.3 Average, largest and smallest performance improvements of the NN+fMPC rela- tive to the CLQR in terms of the cost and the peak and RMS metrics of inter-story drifts, absolute accelerations, and controllable damper forces when the structures are subjected to 20 historical earthquake records (one-controllable-damper examples)101 4.4 Average, largest and smallest performance improvements of the NN+fMPC rela- tive to the CLQR in terms of the cost and the peak and RMS metrics of inter-story drifts, absolute accelerations, and controllable damper forces when the structures are subjected to 20 historical earthquake records (four-controllable-damper example)116 vii List of Figures 1.1 Demonstration of dissipative and nondissipative forces . . . . . . . . . . . . . . . 11 1.2 Comparison of control algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 An example of neural network classifiers . . . . . . . . . . . . . . . . . . . . . . . 12 2.1 Experimental specimen (photo E.A. Johnson) . . . . . . . . . . . . . . . . . . . . 18 2.2 Accelerometer configuration on floors 0–3 . . . . . . . . . . . . . . . . . . . . . 18 2.3 Isolation-layer layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4 Isolation-layer devices (photos E.A. Johnson & T. Sasaki) . . . . . . . . . . . . . . 19 2.5 Shake table center accelerations, in the x and y directions, in Test 010 (random) and Test 016 (scaled March 2011 Tohoku-Oki earthquake) . . . . . . . . . . . . . 20 2.6 Representative linear and nonlinear force-displacement relationships in Tests 010 and 013, respectively (note that the scales differ) . . . . . . . . . . . . . . . . . . . 21 2.7 Identified mode shape vector elements in the complex plane . . . . . . . . . . . . 22 2.8 MAC values of the identified mode shapes against themselves (each natural fre- quency corresponds to two mode shapes that are complex conjugates of each other) 24 2.9 Baseline numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.10 Original and softer model mode shape comparisons: (a) baseline model and one with a softer isolation layer; (b) baseline model with one with a softer superstructure 26 2.11 Finite element model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.12 Mode shapes of the nominal FEM . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.13 Possible mode matchings for the isolation-layer device parameter updating . . . . . 34 2.14 Flowchart for computing superstructure mode-matching cost function values (the damping ratio could be included but is omitted because λ ζ i is set to zero herein; elements drawn assumingλ MAC i is unity) . . . . . . . . . . . . . . . . . . . . . . 35 viii 2.15 Comparisons of frequencies and mode shapes (MAC values) between experimen- tally identified modes and those from the (a) unupdated and (b) updated FEM models 39 2.16 Measured and calibrated-FEM-predicted accelerations, both downsampled to 100 Hz, from the accelerometer near grid 1A at floor 0 (base) in random excitation Test 010: (a) full time duration; (b) narrower time window to show detailed comparison . . . 42 2.17 Measured force-displacement loops, and those predicted by the entire updated lin- ear FEM, of the SDP1 in thex direction in random excitation Test 010 . . . . . . . 43 2.18 RB1x-direction force-displacement loops in Test 016 . . . . . . . . . . . . . . . . 44 2.19 Rubber bearing bilinear elastic stiffness model . . . . . . . . . . . . . . . . . . . . 44 2.20 Measured and calibrated-FEM-predicted accelerations from the accelerometer near grid 1A on floor 0 in Test 013 (scaled Tohoku-Oki earthquake): the left graphs show the full time duration; the right graphs show a narrower time window to show detailed comparison (note: the right z-direction graph uses a narrower time window to show the detail for the higher-frequency vertical motion) . . . . . . . . 50 2.21 Measured and calibrated-FEM-predicted force-displacement loops during the scaled Tohoku-Oki Test 013 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.1 Exact line and noisy coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.2 FEM of 3D frame structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.3 Comparison of FEM and experimental mode shapes as measured by the MAC value 76 3.4 Experimental and updated-FEM-predicted absolute acceleration time histories at the accelerometer near grid 1A on floor 0: full duration and a corresponding zoomed-in version to show some details (the zoomed-in version for z-direction uses a narrower time window for better visibility of the high-frequency motions) . . 80 3.5 Experimental and updated-FEM-predicted force-displacement relationships . . . . 81 4.1 Multi-DOF shear building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.2 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (SDOF system) . . . . . . . . . . . . . . . . . . . . 102 4.2 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (SDOF system) (cont.) . . . . . . . . . . . . . . . . . 103 4.3 Comparisons of the inter-story drifts, the absolute accelerations and the control forces computed by the NN+fMPC with and without force saturation as well as by the CLQR (SDOF system subjected to earthquake No. 15) . . . . . . . . . . . . . 104 ix 4.4 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (3DOF system) . . . . . . . . . . . . . . . . . . . . 106 4.4 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (3DOF system) (cont.) . . . . . . . . . . . . . . . . . 107 4.5 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with one controllable damper) . . . . 109 4.5 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with one controllable damper) (cont.) . 110 4.6 Comparisons of the first story inter-story drifts, the top story absolute accelera- tions and the control forces computed by the hMPC, the NN+fMPC and the CLQR (5DOF building subjected to earthquake No. 5) . . . . . . . . . . . . . . . . . . . 111 4.7 Generation of 3 7 = 2,187 parameter combinations for robustness study (note: the side-branches growing from the side nodes, which are omitted with “··· ”, all grow to the bottom level in the same pattern as the central branch growing from the central node on the same level as the side nodes) . . . . . . . . . . . . . . . . . . . 112 4.8 Cost when the system is controlled by the NN+fMPC and the hMPC algorithms normalized by the corresponding cost when the system with the same parameteri- zation is controlled by the CLQR algorithm . . . . . . . . . . . . . . . . . . . . . 113 4.9 Multi-DOF shear building with multiple controllable dampers . . . . . . . . . . . 114 4.10 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with four controllable dampers) . . . . 117 4.10 Comparisons of the costs and the RMS and peak response and control force met- rics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with four controllable dampers) (cont.) 118 x Abstract Large-scale seismic structural tests are crucial to validating both structural design methodologies and the effectiveness of seismic isolation devices. However, considering the significant costs of such tests, it is essential to leverage data from completed tests by taking advantage of numerical models of the tested structures, updated using data collected from the experiments, to complete additional studies that may be difficult, unsafe or impossible to physically test. Such studies may include responses to additional earthquake records, validations of damage detection techniques, studies of novel seismic control strategies, simulating the effects of retrofits, verification of new design methodologies, etc. However, updating complex numerical models with high-dimensional parameter spaces poses its own challenges. Following the introduction in Chapter 1, two chapters are dedicated to develop model updating approaches suitable for high-order models of a base-isolated structure, motivated and evaluated through modeling and model updating of a full-scale four-story base-isolated reinforced-concrete frame building that was tested in 2013 at the NIED E-Defense laboratory in Japan. Chapter 2 focuses on the modeling of the structural components, the development of a multi-stage model updating algorithm and a mode-matching algorithm. Chapter 3 features the constrained total least squares (CTLS) based cross-model cross-mode (CMCM) model updating method, which is a more accurate and robust substitute for the current state-of-the-art total least squares (TLS) based CMCM model updating method. In most studies involving model updating, all to-be-updated parameters are typically updated simultaneously; however, given the observation that the superstructure in this study predominantly moves as a rigid body in low-frequency modes and the isolation layer plays a minor role in all xi other modes, Chapter 2 proposes updating parameters in stages: first, the linear superstructure pa- rameters are updated so that its natural frequencies and mode shapes match those identified via a subspace system identification of the experimental building responses to low-level random excita- tions; then, the isolation-layer device linear parameters are updated so that the natural frequencies, damping ratios and mode shapes of the three isolation modes match. This two-step process breaks a large-scale linear model updating problem into two smaller problems, thereby reducing the search space for the to-be-updated parameters, which generally reduces computational costs regardless of what optimization algorithm is adopted. Due to the limited instrumentation, the identified modes constitute only a subset of all modes; to match each identified mode with a FEM mode, a pro- cedure is proposed to compare each identified mode with a candidate set of FEM modes and to select the best match. Finally, the nonlinear isolation-layer device models are updated based on their force-displacement relationships under both random excitations and strong ground motions as well as the force time histories of the resulting combined nonlinear system model. Combining the isolation-layer devices’ nonlinear models with the updated superstructure linear FEM, the final result is a data-calibrated nonlinear numerical model suitable for, and that will be used for, further studies of controllable damping and validation of new design methodologies, and is being made available for use by the research community. Continuing the investigation of model updating algorithms, Chapter 3 focuses on a novel solu- tion to the CMCM model updating method, which is a physics-based method that avoids the diffi- culty of matching experimental modal information computed from measured responses with that of a corresponding FEM mode. Due to measurement noise in, and the spatial incompleteness of, the experimental mode shapes, this method involves solving a set of ill-conditioned linear equations. The state-of-the-art approach to solve these equations is based on the total least squares (TLS) ap- proach. Instead, this study proposes replacing the TLS solution with one based on the constrained total least squares (CTLS), and shows that CTLS can provide more accurate solutions given the same level of noise and spatial incompleteness. To demonstrate this observation, the TLS and the CTLS approaches are first used to solve a simple line-fitting problem with noisy measurement xii data; the pros and cons of the approaches are highlighted through this rudimentary example. Then, they are both applied to the model updating of a three-dimensional (3D) frame structure using a structure model from previous CMCM studies. Results indicate that the CTLS approach is a better alternative to the TLS approach when dealing with the set of ill-conditioned linear equations in the CMCM method. Finally, the CTLS-based CMCM method is applied to update a large-scale FEM of a full-scale four-story base-isolated reinforced-concrete frame building. Results affirm that the proposed method is applicable to real engineering practice. Starting from Chapter 4, the focus of this dissertation is shifted to the development of a real- time neural network based semiactive model predictive control strategy for reducing structural vibrations. Semiactive model predictive control (MPC) algorithms perform well in a range of sce- narios. However, because their implementation is computationally time-consuming due to the in- herent mixed integer quadratic programming (MIQP) optimization problem, they have never been adopted in real-time structural vibration control. Inspired by neural network (NN) approaches for MIQP solutions and a fast quadratic programming (QP) solver, this paper proposes that NN clas- sifiers be trained to predict in real-time only the integer variable values (each unique combination of these is called a strategy) for a given current state variable vector; during real-time operation, the NN predicts the integer variables, reducing MPC control force optimization to a much simpler QP problem, which is solved using a fast QP solver adapted to this type of problem. The number of possible strategies grows exponentially with the number of integer variables, which itself grows with the numbers of prediction horizon steps and damping devices; accurately matching the current state vector (the input to the MIQP) with one of these strategies typically requires a massive NN, which is time-consuming to train offline and, more importantly, may be too computationally inten- sive to utilize online. To reduce the number of possible strategies/predictions, the control algorithm proposed herein first exploits the homogeneity-of-order-one nature of this structural control prob- lem and the statistics of the system states to shrink the input parameter space, leading to a smaller set of corresponding output strategies; additionally, this algorithm uses a few smaller trained NNs running simultaneously, each linking an input to a unique group of strategies, instead of one large xiii NN that takes much longer to run, to jointly predict the correct strategy for a given state vector input to the NN. Further, occasional inaccurate NN predictions can lead to sub-optimal solutions; to enforce a bound on the sub-optimality, the control algorithm will switch to a simpler clipped linear quadratic regulator (CLQR) algorithm, if the CLQR cost function is smaller than that cal- culated based on the predicted strategy. Numerical examples using shear buildings of one, three and five degrees-of-freedom (DOFs) demonstrate that the proposed real-time algorithm speeds up the online computation by approximately one order of magnitude and offers a control performance close to the state-of-the-art offline semiactive MPC algorithm. Future planned work is listed in Chapter 5, together with the importance and feasibility, before conclusions are drawn in Chapter 6. The future work includes: (1) exploring new techniques so that the real-time semiactive control strategy can be applied to structure with a significantly greater number of DOFs; (2) combining the established control strategy with algorithms that estimate the changing structural parameters, so that the control strategy can offer control signals adaptive to the changing structural parameters, maintaining the robustness of performance; (3) applying the online semi-active model predictive control strategy in the context of a decentralized control paradigm in which each controller is only responsible for controlling a local subset of DOFs; (4) combining the established control strategies with dual-control methods to determine control forces that can reduce the vibration of structures while helping with the estimation of structural parameters. xiv Chapter 1 Introduction This chapter focuses on the motivations of my work and reviews previous research investigations in the corresponding fields. Since my work is three-folded, this chapter is divided into three parts, each focusing on one topic of my research. 1.1 Introduction to Modeling and Model Updating of a Full- Scale Experimental Base-Isolated Building In many countries around the world, including the United States, earthquakes have caused severe damage to buildings and infrastructure and are one of the major concerns in the field of civil en- gineering. Many protective measures have been taken to minimize the seismic responses of struc- tures; among these measures, base isolation has been widely adopted in recent decades (Constanti- nou et al., 1990; Mokha et al., 1993; Sato et al., 2013). Although the behaviors of base-isolated buildings have been extensively analyzed in previous research, full-scale seismic tests of these buildings are limited in quantity due to the high costs of such experiments (Brewick et al., 2018). However, such tests, like the series of tests investigated in this study, are vital to understanding the true performance of structures and for validating and innovating new base isolation techniques. Although some conclusions drawn from the analysis of a particular base-isolated building may 1 be unique to that building, other conclusions can deepen our knowledge of general base-isolated building behavior. Given the high costs of full-size seismic tests, it is essential to leverage data from completed tests by utilizing a test structure’s numerical models, that have been subsequently calibrated from the experimental data, to pursue additional studies that could not be tested in the laboratory or in the field. Such studies may include responses to additional earthquake records — including those too severe to safely test or even beyond laboratory capability — as well as damage detection tech- nique validations, novel seismic control strategies, retrofit effects, and new design methodology verification. In these scenarios, finite element models (FEMs) of the experimental structure can be formulated and become powerful tools for researchers to evaluate the structural performance in subsequent studies. However, before ensuring that the model’s key features (natural frequencies, mode shapes, energy dissipation pattern,etc.) are sufficiently close to those of the tested building, the trustworthiness of model results is indeterminate. Thus, however challenging it may be, model updating is a vital step between the experiments and subsequent model use in future studies. Model updating methods may be characterized into four categories of methods: modal (Moaveni et al., 2009; Bakir et al., 2007, 2008; Friswell and Mottershead, 2013), frequency response function (FRF) (Lin and Ewins, 1994; Imregun et al., 1995; Lin and Zhu, 2006; Lin, 2011; Shadan et al., 2016), modal strain energy (MSE) (Wang and Xu, 2019; Stubbs et al., 1995; Shi et al., 1998; Yang et al., 2020; Cornwell et al., 1999), and probabilistic (Bayesian) (Beck and Arnold, 1977; Beck and Katafygiotis, 1998a,b; Vanik et al., 2000; Beck et al., 2001; Jang and Smyth, 2017a). In modal methods, the model parameters are updated so that the FEM’s modal properties (natural frequen- cies, damping ratios and mode shapes) match those of the experimental structure. FRF methods — proposed in Lin and Ewins (1994) for analytical models and in Imregun et al. (1995) for numer- ical models, and utilized in other studies (Lin and Zhu, 2006; Lin, 2011; Shadan et al., 2016) — update the model parameters so that the FEM transfer functions reproduce those computed from the experiment’s input and output measurements; thus, modal property estimation, through system 2 identification from the experimental data is not required (as it is for modal methods). MSE meth- ods explore the differences in modal strain energy between the FEM and the experimental structure to update the model parameters; since proposed in 1995 (Stubbs et al., 1995), MSE methods have been successfully applied to model updating of beam structures (Yang et al., 2020) and plate struc- tures (Cornwell et al., 1999), and a variant called the cross-MSE method (Hu et al., 2007, 2006; Wang et al., 2015) was successfully applied to three-dimensional truss structures. It is believed that the modal strain energy is more sensitive to model parameter changes, though it is also sensitive to noise, which limits its effectiveness in tests with noisy experimental data. Bayesian methods, the most common stochastic model updating methods, utilizes Bayes’ inference theorem to pro- vide the probability distribution of (instead of just the best) model parameters, making it robust to uncertainties. However, the significant computational cost to update the model of a complicated structure can be a challenge; for example, updating the 22 parameters of a bridge structure FEM, constructed with 19,632 beam elements and 1,464 truss elements, took three months to generate enough posterior samples (Jang and Smyth, 2017a). Although these applications of the aforemen- tioned methods used linear model updating, some studies (e.g., Song et al., 2009a,b, 2012) have also incorporated nonlinearities (e.g., nonlinear material models), but are often applied only to laboratory-scale data. In this study, several approaches are proposed for updating a partially nonlinear, high-order FEM of a base-isolated structure, motivated by and applied to updating a detailed model of a full- scale base-isolated structure using response data from the specimen’s 2013 shake table tests at Japan’s NIED “E-Defense” facility (Nakashima et al., 2018). The FEM involves approximately 85,000 degrees of freedom (DOFs); this level of model complexity effectively precludes capitaliz- ing on Bayesian probabilistic methods because of computational cost. The MSE-based total least square cross-model cross-mode linear updating (Li et al., 2020) was attempted but resulted in a poor match of the FEM and experimentally measured modal properties, likely due to MSE’s sen- sitivity to the the sensor noise in this series of experiments. This leaves FRF and modal methods to consider. 3 One aspect that is particular to this structure is its base-isolated, not fixed-base, design that results in the superstructure moving approximately as a rigid body in the first three modes (iso- lation modes) and the deformations in all other modes (superstructure modes) mostly a result of superstructure behavior. For a base-isolated building like this, the superstructure model parameters can be updated first, and the isolation-layer parameters can be updated subsequently. As a direct result, one larger optimization problem is broken into two smaller problems (fewer modes to match and a smaller parameter space) that are easier and less time-consuming to solve (a previous study trying to update all parameters together took much longer and resulted in a sub-optimal match to the experimental results). Because modal methods treat each mode separately, and naturally sepa- rate the isolation modes from the superstructure modes, they are better choices than the FRF-based methods for the model updating of the base-isolated building in this study. Nevertheless, there are several challenges in applying this proposed two-stage modal-based model updating method; these hurdles and the approaches proposed herein to address them are: 1. The rationale for separating the update into two stages must be established. To do so, this study introduces a simple example with a stiffness distribution similar to that of the ex- perimental building, and empirically demonstrates that: (a) isolation-layer linear stiffness changes have minor, and typically negligible, effects on the superstructure modes’ natural frequencies and mode shapes; (b) superstructure linear stiffness changes have only marginal effect on the isolation modes’ natural frequencies and mode shapes; and (c) superstructure damping changes have negligible influence on the isolation modes’ damping ratios. 2. While the isolation modes are likely non-classically damped, due to the significant differ- ences in stiffness and damping in the isolation layer relative to the superstructure, the super- structure modes are assumed to be proportionally damped. An analysis of the experimentally- identified mode shapes demonstrates that this assumption is reasonable. 4 3. Successful modal-based model updating reduces the error between each experimentally- identified mode and its corresponding FEM mode, which is often not the same mode num- ber (when modes are sorted by increasing frequency) because: the unupdated FEM model’s modes may be in a different order due to uncalibrated model parameters; identification from a limited sensor set may not be able to detect some modes; and measurement noise can result in spurious non-physical “identified” modes. While previous studies have success- fully matched experimental and model modes for relatively simple structures based on fre- quency differences, different definitions of modal assurance criteria (Allemang and Brown, 1982b; Lieven and Ewins, 1988), mass matrix compatibilities (Avitabile and O’Callahan, 1988; O’Callahan, 1995), comparisons of modal energies (Bugeat and Lallement, 1976) and comparison of FRFs (Ibrahim, 1993), or for ones in which the mode matching is straight- forward (e.g., Jang and Smyth, 2017b), this specimen’s limited sensor coverage and the high-order FEM poses a significant mode-matching challenge. To address this difficulty, a new efficient and effective mode-matching procedure is proposed. 4. Modal approaches are inherently based on linear system characteristics, but this structure’s isolation layer is highly nonlinear, as is typical of most isolated buildings; thus, it must be determined whether superstructure nonlinearity is a factor. Herein, it is demonstrated for this structure that the identified natural frequencies of the corresponding superstructure modes, based on the random and the earthquake excitation test data, are nearly identical, which is consistent with base-isolation design intent for superstructure deformation; thus, this su- perstructure’s nonlinearity is negligible during significant earthquake excitation, even with peak horizontal accelerations as much as 0.5g. Finally, to simulate the building’s nonlinear responses, the calibrated superstructure model is combined with isolation-layer device mod- els — some bi-directional hysteretic parameter-calibrated (Brewick et al., 2020) Bouc-Wen models (Park et al., 1986; Wen, 1980) and some bilinear stiffness models. The result is a locally nonlinear FEM that can faithfully reproduce the building’s responses to earthquake ground motions. 5 In summary, the model updating is conducted through the following steps. (a) Perform sys- tem identification using experimental input-output (accelerometer) measurements to estimate the modal properties of the building. (b) Update the superstructure linear stiffnesses so that the natural frequencies and mode shapes of the FEM superstructure modes match those identified. (c) Update the linear stiffnesses and damping coefficients of the isolation-layer devices so that the natural frequencies, mode shapes and damping ratios of the FEM isolation modes match those identified. (d) Perform time-history analyses and compare the simulated responses with the experimental data to evaluate and demonstrate the success of model updating. (e) Calibrate nonlinear models of the isolation-layer devices, combine them with the superstructure model, and evaluate the result- ing partially nonlinear FEM through its time-history responses to demonstrate model calibration success. The remainder of this paper is organized as follows. Section 2.1 provides a brief review of the experiments with a focus on the experimental set-up used in this study. Section 2.2 summarizes the relevant system identifications that were partially reported in Brewick et al. (2018) as well as those new in this study. Section 2.3 explains and validates the aforementioned assumptions and challenges. Sections 2.4–2.6 detail the linear model updating, including an introduction to the FEM, the formulation of the optimization algorithm, a summary of the updated FEM and a time-history analysis using the updated model. Section 2.7 focuses on the construction of the partially nonlinear FEM, including the isolation-layer devices’ modeling and comparisons of the model-predicted and measured earthquake responses. 1.2 Introduction to Cross-Model Cross-Mode Model Updating Method Based on Constrained Total Least Squares A dictionary of methods have been developed in the past few decades to update a finite element model (FEMs) using experimental data (Jang and Smyth, 2017b; Moaveni et al., 2009; Bakir et al., 2007, 2008; Song et al., 2009a, 2012, 2009b; Hu et al., 2007, 2006). Among these, the cross-model 6 cross-mode method (CMCM) shows some distinct advantages; the basics of CMCM are provided by Hu et al. (2007, 2006), and it has been applied to simultaneous updating of stiffness, mass and damping matrices (Hu and Li, 2007) and a practical engineering application (Wang et al., 2015; Li et al., 2020). The first advantage of CMCM is that it is a physics-based method, which can reveal the nature and magnitude of the change of a physical parameter, as opposed to data-driven methods that have nonphysical parameters. Second, CMCM avoids the challenge in some model updating approaches of determining which FEM mode to match with which experimental mode, which is often a problem with no exact solution but multiple imperfect and competing solutions. In particular, the set of experimental modes may be incomplete because sensor placement may preclude detecting some modes; further, those modes discovered are spatially sparse relative to the FEM. With undiscovered modes, the experimental and FEM modes cannot be matched in a sequential fashion when modes are numbered by increasing frequency. (For example, a FEM may have modes that are dominated by vertical vibration of floor slabs, which may be impossible to de- tect if the experiment included no sensors on the slabs.) Third, CMCM is far more computationally efficient than methods involving the genetic algorithm or the Bayes’ inference theory — typically used when the cost function has no gradient — making CMCM a powerful alternative for updating models of full-scale civil structures with high-dimensional parameter spaces. Despite these appealing characteristics, CMCM still poses one technical difficulty. CMCM explores the equations for the modal strain energy balance of the experimental structure, in which the mass, damping and stiffness matrices’ deviations from those of a nominal FEM are parameter- ized by a set of magnitude coefficients of (generally low-rank) matrix perturbations, and aims to determine the magnitude coefficients that maintain the energy balance. When the unknown mag- nitudes are assembled into a column vector x, the set of modal strain balance equations can be cast as Ax= b. The formulation of this set of linear equations requires the experimental mode shape data at all degrees of freedom (DOFs)as shown in Section 3.1, which can only be obtained by expanding the available experimental mode shapes at limited sensor locations. Because of the necessarily imperfect modal expansion procedure and the noise in experimental modal information 7 identified from (noisy) measured structural responses, both A and b are inaccurate, making Ax= b ill-conditioned and challenging to solve. The simplest solution to CMCM’s resulting set of linear equations of the form Ax= b is a least squares (LS) solution for x that essentially involves the Moore-Penrose pseudo-inverse of A; LS is generally effective when the matrix A is known exactly and only b is uncertain, i.e., A(x+∆x)= b+∆b, where∆b represents some error. The more robust and accurate approach, used in CMCM (Li et al., 2020) prior to this study, is to solve using the total least squares (TLS) approach that assumes there may be errors in both matrix A and vector b —i.e.,(A+∆A)(x+∆x)= b+∆b — which is certainly the case in CMCM as both A and b depend on uncertain modal information, as demonstrated in Section 3.1. However, TLS inappropriately assumes (Golub and Van Loan, 1980; Van Huffel and Vandewalle, 1985; Silvia and Tacker, 1982) that the errors∆A in matrix A and the errors∆b in vector b are uncorrelated; this assumption is clearly violated (Section 3.1), and still introduces errors to the model updating solution. The major contribution of this paper is to account for the correlated errors in ∆A and ∆b by introducing the constrained total least squares (CTLS) approach to solve CMCM’s ill-conditioned linear equations for accurately (Markovsky and Van Huffel, 2007) updating a FEM. The CTLS was first developed in the field of signal processing (sound and image restoration) (Abatzoglou and Mendel, 1987; Abatzoglou et al., 1991; Yang et al., 2009; Ng et al., 2002; Mesarovic et al., 1995), and its solution was obtained by solving an optimization problem using either the Newton’s approach or an iterative approach involving quadratic programming. The Newton’s approach may lead to a more accurate solution but without guaranteed convergence. Thus, the iterative approach is used with crucial additional steps proposed herein to adapt the CTLS to CMCM problems. Since the inception of the CTLS, beyond its applications in electrical and computer engineering, it has also been successfully implemented in geodesy to compute coordinates of certain locations (Chen et al., 2011; Tsai and Kao, 2006), and in biochemistry to identify genetic networks (Guner et al., 2015). Surprisingly, despite its widespread applications in various fields of studies, to the authors’ 8 knowledge, the CTLS has not yet been used in the model updating and structural health monitoring communities, in which researchers frequently deal with uncertainties in parameters. A second contribution of this paper is to expand the experimental mode shapes based on the updated FEM mass, damping and stiffness matrices instead of those before the current round of model updating. This formulation may be mathematically complex to derive, but is a more math- ematically sound formulation for modal expansions, especially when the near-true mass, damping and stiffness matrices of the experimental building are obtained after just one round of model updating (in general, model updating may be performed a couple of times until the perturbation magnitude coefficients are small but, in certain cases like the experimental example in Section 3.4.1 when updating the isolation-layer device parameters, only one iteration is necessary). A direct consequence of these contributions is that the CTLS-based CMCM method produces more accurate results than its TLS-based counterpart (Li et al., 2020), and works well with reason- able level of measurement noise even without the strong regulation procedure (discarding solutions that do not satisfy some preconceived expectations) introduced in Li et al. (2020); in fact, eliminat- ing this regulation procedure is necessary for model updating applications in which the magnitude coefficients can legitimately be positive or negative. Chapter 3 is organized as follows. Section 3.1 provides a brief introduction to CMCM to famil- iarize the reader with the necessary background. Section 3.2 explains the possible algorithms — including the TLS, the CTLS and the weighted total least square (WTLS) — to solve the set of uncertain linear equations (A+∆A)(x+∆x)= b+∆b (in some mathematical literature (Carroll et al., 2006), this equation is called the error-in-variable, or EIV , problem). The state-of-the-art for CMCM is the TLS solution; thus, this section compares the effectiveness of the TLS, the CTLS and the WTLS (which has been posed as a possible alternative to the CTLS) algorithms for a line-fitting problem with noisy data and analyzes their respective advantages and disadvantages. Section 3.3 compares the efficacy of the TLS-based and the CTLS-based CMCM methods for solving an analytical problem, a 3D frame structure used in previous CMCM studies (Hu et al., 9 2006; Li et al., 2020), to demonstrate that the CTLS-based method provides more accurate solu- tions than the TLS-based method at the same noise intensity and level of spatial incompleteness. Section 3.4 is dedicated to a practical model updating example utilizing the CTLS-based CMCM method to update a large-scale FEM based on experimental responses from a full-scale civil engi- neering structure, demonstrating the effectiveness of the CTLS-based CMCM method in real-world applications. 1.3 Introduction to Real-Time Neural Network based Semiactive Model Predictive Control of Structural Vibration Civil engineering structures are vulnerable to earthquakes and strong winds. The traditional way of earthquake prevention design is to increase the ductility of the structures, so that they deform but do not collapse in case of extreme events. However, such an approach can be cost prohibitive in areas subjected to severe earthquake ground motions. Also, these structures usually suffer damage in the course of extreme events, and the structural retrofits can be challenging and expensive. A more cost effective approach to reduce the vibration of the structures is to implement vibration control devices dictated by state-of-the-art control algorithms. In the control engineering community, various control strategies have been developed, which can generally be sorted into three categories: passive, active, and semiactive. In the field of vibration control for civil structures subjected to earthquakes or strong winds, semiactive controllable damping strategies are preferred because they can often provide far better control performance compared with passive control algorithms, consume less energy compared with active control algorithms, and are inherently stable because the control forces are naturally dissipative (the control forces always act to oppose motion across the device, as shown in Fig- ure 1.1); in certain cases, controllable damping can perform just as well as active control algo- rithms (Dyke et al., 1996). Some practical applications of semiactive clipped optimal control in civil engineering first determine control forces by ignoring the dissipativity constraint through the 10 use of linear control algorithms such as a linear quadratic regulator (LQR) (Figure 1.2a), and then clip (e.g., set to zero, or as close to zero as is possible) the non-dissipative forces with a secondary control device. Because the primary control device does not explicitly consider the constraint, the separate clipping mechanism takes a toll on the original (pre-clipping) control performance, especially when frequent clipping is required. One control algorithm works around this issue by using a dissipative force guaranteed to be, at least instantaneously, superior to an optimal passive control (Scruggs et al., 2007). Another control algorithm that avoids clipping underperformance is the semiactive model predictive control (MPC) or hybrid MPC (hMPC), in which integer variables are introduced to explicitly enforce the dissipative constraint (Figure 1.2b). However, finding the hMPC solution involves solving mixed-integer quadratic programming (MIQP) problems, making it too slow for real-time implementation in structural vibration control problems. In Elhaddad and Johnson (2013), to move most calculations offline, a pre-calculated look-up table was introduced to find the control forces for systems with up to 4 states; while the control performance was effec- tive, it would be infeasible to store and efficiently access the look-up table when the state-space is larger: e.g., the look-up table for a ten-state system would occupy computer memory on the order of petabytes. damper velocity control force Dissipative Dissipative Nondissipative Nondissipative Figure 1.1: Demonstration of dissipative and nondissipative forces One way of accelerating the hMPC while not overloading the memory is to apply the algorithm in Bertsimas and Stellato (2019, 2021): create a list of candidate combinations of integer variable values and inequality constraints that are active at optimality (active constraints) — each unique combination of integer variable values and active constraints is termed a “strategy” and denoted 11 LQR Clipping device Measurements Optimal control force Dissipative control force (a) Clipped LQR hMPC Measurements Optimal dissipative control force (b) Hybrid MPC Figure 1.2: Comparison of control algorithms x NN 0,1 x NN 0,2n P(S 0 ) P(S M− 1 ) . . . NN input that is uniquely determined from state vector x 0 . . . Hidden layer(s) Hidden layer(s) Input layer Output layer Figure 1.3: An example of neural network classifiers herein byS s II ,s= 0,1,...,M− 1 — and train offline a feedforward neural network (NN), as in Fig- ure 1.3, to predict online the probabilitiesP(S s II ) that a strategy S s II leads to the optimal solution when given the current state vector x 0 . Adopting the strategy S s ∗ II , where s ∗ = argmax s P(S s II ), that is predicted most likely to be correct, the online computation of controllable damping forces simplifies to solving a set of linear equations. However, a direct implementation of this approach is difficult due to two challenges. First, because the NN’s input vector (the current state vector) is generally high-dimensional and the number of possible strategies is huge, many samples are re- quired to effectively train the NN, which is computationally expensive; for example, a preliminary study that trained a NN to predict the probabilities of the S II strategies for the ten-state numerical example in Section 4.4.3 required about 19 CPU-years (over two weeks wall-clock time) on a mas- sively parallel cluster in the Center for Advanced Research Computing (CARC) at the University 12 of Southern California (USC); for someone with only a standard personal computer, such a control design task is infeasible. Second, because there are too many strategies —e.g., in the same prelim- inary study, the training samples constituted millions of unique S II strategies — any NN classifier with a fair performance must be of gigantic scale, with many hidden layers and many nodes, again leading to both long training time and prohibitive prediction time. To address these challenges, the algorithm proposed herein both uses multiple smaller NNs in parallel via the “codeword” approach (Zhang et al., 2018) and predicts only the values of the integer variables, not the active constraints, which reduces the online optimization problem from MIQP to quadratic programming (QP), which is more complex than just a set of linear equations but can be efficiently solved by the fast MPC (fMPC) algorithm (Wang and Boyd, 2010) adapted with certain variations. The idea of predicting the integer variables was mentioned previously in the field of robotic control (Cauligi et al., 2020), but was also independently developed and presented by the first author (Yu, 2020) a few months before Cauligi et al. (2020) was made public. Further, the problems addressed in Cauligi et al. (2020) are of far smaller scale compared with the MIQP problem in structural vibration control; thus, there are challenges unique to large-scale problems that Cauligi et al. (2020) never attempted to handle: 1. Straightforward NN training would use a vector of the system states as the input and a vec- tor of the probabilities that each candidate strategy can lead to the optimal solution of the MIQP problem as the output. However, by exploiting the homogeneity-of-order-one nature of the control algorithm (i.e., control forces scale with the magnitude of the state vector), the prediction does not change when the input vector is proportionally scaled. Thus, all input state vectors for NN training are normalized so that their two-norms are unity, leading to a moderate shrinkage of the sampling space and less effort for sampling. 2. To train the NNs on the resulting unit hypersphere in a manner that approximately uniformly represents the relative magnitudes of the states, a purely offline hMPC algorithm is run for the simulated structure subject to a Gaussian random excitation to approximate the state 13 response covariance matrix, which will be used to generate state vector samples for NN training. 3. Although the applications in Cauligi et al. (2020) permit a certain group of integer variables to be determined by a corresponding group of input parameters, thereby allowing some un- correlated small NNs to be trained offline and applied online separately and simultaneously, this simplification does not apply to the semiactive controllable damping strategies in this study; thus, the NN remains very large and computationally prohibitive. To cut both training time and prediction time, a “codewords” (Zhang et al., 2018) paradigm is introduced to re- place a single large NN with multiple smaller NNs, each predicting the likelihoods of each strategy within a unique set of strategies, chosen so that the NNs together uniquely predict individual strategies: strategy number s is represented as a set of digits of some base, and each smaller NN predicts one such digit to uniquely determines. 4. A small-scale problem can be solved efficiently by a commercial optimization solver after the integer variables are available; however, available commercial solvers are too slow for the real-time implementation required herein. To guarantee sufficient computing speed, a C-language implementation of the fMPC in Wang and Boyd (2010) is tailored to find the control forces after the integer variables are predicted by the NN(s). In addition to these considerations, three techniques are proposed to promote accuracy and/or the computational efficiency of the proposed algorithm. 5. Instead of using a fixed penalty coefficient when solving the QP problem ( i.e., when imple- menting the fMPC), this study proposes a way to adapt the penalty coefficient at every time step, promoting the accuracy of the QP solution. 6. An easily computed, guaranteed dissipative, but not necessarily optimal semiactive control strategy (run in parallel with the proposed algorithm) is taken as the initial solution for the QP solver for the fMPC algorithm, speeding up the computation in the same fashion as 14 what some may call a “warm-start” technique (Wright, 2000; Yildirim and Wright, 2002) in solving optimization problems. 7. Previous studies never considered the effects of inaccurate NN predictions, which are in- evitable in real practice; further, the fMPC optimization can, occassionally, result in a non- converged sub-optimal solution. To address both of these, this study proposes using several strategies predicted most likely by the NNs, not just the one most probable, comparing the those strategies’ predicted costs with each other and against that from an alternate easily computed semiactive control strategy, and choosing at each time step the one that leads to the smaller predicted cost. This ensures that the proposed algorithm will always perform at least comparably to, and usually better than, the suboptimal-but-easy strategy. Herein, the alternate strategy used in #6 and #7 is the clipped linear quadratic regulator (CLQR) solution; however, other strategies could be used, such as a passive strategy or a semiactive strategy such as that proposed by Scruggs et al. (2007). These proposed solutions to the challenges enable an effective real-time semiactive MPC for structural vibration control. Simulations herein using single-degree-of-freedom (SDOF), 3DOF and 5DOF structure models show that the proposed algorithm can achieve control performance nearly identical to that of the hMPC but at a much faster speed that is sufficient for real-time implementation. For surveys of structural control, the interested reader is directed to Housner et al. (1997), Soong and Spencer (2002), Symans et al. (2008), and particularly Spencer and Nagarajaiah (2003) that has a special focus on semiactive control. For applications of MPC in structural vibration control and in a range of other industries, see Mei et al. (2001, 2002, 2004) and Qin and Badgwell (2003), respectively. Modeling systems with “on-off” logic (e.g., logic enforcing the dissipativity constraints) through MIQP was covered in Sontag (1981), Bemporad and Morari (1999), Heemels et al. (2001) and Bemporad (2002). Transforming a semiactive structural vibration control algo- rithm into a MIQP problem was detailed in Elhaddad and Johnson (2013). 15 The idea of training NNs offline to accelerate solution of online optimization problems was first proposed in Bertsimas and Stellato (2019, 2021). A similar study (Cauligi et al., 2020) was conducted independently and simultaneously with this study, but is incapable of handling large- scale problems arising from structural vibration control. Separating one huge NN into several smaller NNs was studied in Zhang et al. (2018), which is a generalization of the error-correcting output codes (ECOC) method (Dietterich and Bakiri, 1994; Allwein et al., 2000; Passerini et al., 2004) that has various applications (Wright, 1997; Ghani, 2001; Escalera et al., 2009; Yang et al., 2015; Qin et al., 2017). A fast solver for the active MPC (with no semiactive constraint) was introduced in Wang and Boyd (2010), which exploited the block-diagonal structure of the matrices (Wright, 1997; Boyd and Vandenberghe, 2004), applied a warm-start technique (Wright, 2000; Yildirim and Wright, 2002), and made approximations to accelerate the computation. 16 Chapter 2 Modeling and Model Updating of a Full-Scale Experimental Base-Isolated Building 2.1 Experimental Set-up A 686-ton full-scale four-story base-isolated reinforced-concrete frame building (Figure 2.1) was tested in 2013 (Sato et al., 2013) on Japan’s NIED (National Research Institute for Earth Science and Disaster Resilience) 6DOF shake table at the Hyogo Earthquake Engineering Research Center, best known as “E-Defense” (Ogawa et al., 2001). The geometric asymmetry and the corner stairway-core walls create significant lateral-torsional coupling in the superstructure responses. For the 8 August 2013 experiments concerned in this study, the isolation layer consisted of two rubber bearings (RBs), two elastic sliding bearings (ESBs) and two passive U-shaped steel yielding damper pairs (SDPs) (Figures 2.3 and 2.4). A series of experiments were conducted in August 2013 with excitations including low-intensity filtered white noise as well as scaled versions of both historical and synthetic earthquake records. Various types of sensors — including displacement sensors, accelerometers and force trans- ducers — were mounted on the superstructure and the shake table to monitor the time-history responses of the entire system. Specifically, tri-directional accelerometers were placed at three cor- ners on floors 0–3 (Figure 2.2) — where floor 0 is the superstructure base — and at two corners on the roof (floor 4), for a total of 14 locations and 42 (14 × 3) superstructure acceleration channels (x, 17 Figure 2.1: Experimental speci- men (photo E.A. Johnson) Stairs 1 2 3 A B C Sensor Sensor Sensor x y Figure 2.2: Accelerometer configuration on floors 0–3 y and z directions). Further, the shake table was instrumented with four tri-directional accelerom- eters, for a total of 12 (4× 3) shake table acceleration channels, from which the 6× 1 table-center acceleration vector ¨ x g (t) — composed of three translational and three rotational accelerations — is computed in a least squares sense, assuming that the table is rigid, using a transformation ma- trix determined from the sensor location geometry (Brewick et al., 2018). Force transducers were mounted under each isolation-layer device to measure the restoring forces. Unless otherwise noted, all signals were collected at a 1 kHz sampling rate and then low-pass filtered with a 35 Hz cutoff frequency (below the 60 Hz electrical frequency used at the testing facility). This study focuses on the series of tests conducted on 8 August 2013 — including Tests 010– 012 (different realizations of a random excitation) and Tests 013–016 (scaled versions of the March 2011 Mw9.0 Tohoku-Oki earthquake record from the KNET Furukawa station). Figure 2.5 shows the horizontal table-center accelerations in random excitation Test 010 and earthquake excitation Test 016 (the most intensive, as ranked by root-mean-square ground acceleration, of the Tohoku- Oki earthquake tests). The excitation records’ peak accelerations and durations are provided in Table 2.1. When subjected to the random excitation, both the superstructure and the isolation-layer de- vices remain linear so that the isolation-layer devices can be reasonably modeled with a linear 18 0 1 2 3 4 5 6 7 8 x [m] 0 1 2 3 4 5 6 7 8 9 10 y [m] x y TA01 TA02 TA03 shake table TA04 base overhang structure base ESB2 RB2 RB1 ESB1 SDP1 SDP2 D02 D01 D03 D04 tridirectional accelerometer laser displacement sensor Figure 2.3: Isolation-layer layout RB ESB SDP Figure 2.4: Isolation-layer devices (photos E.A. Johnson & T. Sasaki) stiffness and a constant damping coefficient in each horizontal direction. The earthquake excita- tions, however, as discussed in Section 2.2.2, induce mild nonlinearity in the rubber bearings and strong nonlinearity in the elastic sliding bearings and the U-shaped steel damper pairs (Figure 2.6), though the superstructure remains essentially linear. 19 Figure 2.5: Shake table center accelerations, in the x and y directions, in Test 010 (random) and Test 016 (scaled March 2011 Tohoku-Oki earthquake) Table 2.1: Descriptions of excitations Peak Acceleration (cm/s 2 ) Test No. Excitation x y z Duration (s) 010 6DOF Random 52.08 57.98 68.41 172.06 013 243.89 300.32 170.10 314.90 014 415.81 432.98 363.74 314.91 015 Tohoku-Oki EQ 347.66 480.95 224.81 314.93 016 356.21 398.14 222.99 314.95 2.2 System Identification To estimate the modal properties from the experimental data — so that the FEM can be updated in Section 2.5 to match those properties — identification is performed based on the experimental data from Test 010, in which the building was subjected to low-level random excitations in all six shake table DOFs and remained primarily linear in both superstructure and isolation layer. 20 Figure 2.6: Representative linear and nonlinear force-displacement relationships in Tests 010 and 013, respectively (note that the scales differ) Further, given that calibrated energy dissipation is crucial for reproducing time-history responses and that the damping ratios in low-level motion might be different from the effective damping in large earthquake response, system identification is also performed based on data from Test 013, in which the input is a scaled March 2011 Tohoku-Oki earthquake record (Table 2.1), to determine effective natural frequencies and effective damping ratios. 2.2.1 Test 010 The primarily linear responses to the low-level random excitations were used, in a previous study (Brewick et al., 2018), to estimate modal properties (natural frequencies, damping ratios and mode shapes) through the N4SID (Subspace State Space System Identification) system identification method (Van Overschee and De Moor, 1994). The 12 table acceleration responses were used as inputs to the building model and the 42 base and superstructure acceleration responses were the 21 outputs. The table and structure accelerations were detrended, additionally low-pass filtered at 28.6 Hz (8th-order Chebyshev Type I filtered, first forward and then backward) and downsampled from 1 kHz to 71.4 Hz. When implementing the N4SID method, selecting the order of the model is always a challenge: too low and some physically important modes may fail to be identified; too high and spurious or redundant modes may also be “identified.” Thus, the stabilization diagram strategy (Peeters, 2000) was used to choose the order and identify “stable” modes, which are de- fined to be modes identified at the optimal model order that remain, at the next higher model order, within 1% frequency deviation and 5% damping ratio deviation, and have modal assurance crite- rion (MAC) values larger than 0.95 (stabilization diagrams for Test 010 are shown in Brewick et al., 2018). The MAC(φ φ φ i ,φ φ φ j )=|φ φ φ H i φ φ φ j |/(φ φ φ H i φ φ φ i φ φ φ H j φ φ φ j ) 1/2 (superscript H denotes conjugate transpose) ranges from 0 to 1, where 0 indicates no similarity (orthogonal vectors) and 1 indicates a perfect match (parallel vectors) (Peeters, 2000; Vacher et al., 2010). The identified frequencies and damping ratios with model order 66 are shown in Table 2.2. The first three modes have comparatively lower natural frequencies as intended, and are labeled “isolation modes.” The remaining modes are labeled “superstructure modes.” Because there are Re Im M o d e 1 M o d e 2 M o d e 3 M o d e 4 M o d e 5 M o d e 6 M o d e 7 M o d e 8 Figure 2.7: Identified mode shape vector elements in the complex plane 22 42 output channels, each identified mode shape is a vector with 42 complex-valued entries; Fig- ure 2.7 shows the directions of the mode shape’s elements in the complex plane for each of the first eight identified modes. The mode shapes for identified modes 4–6, which move primarily in the horizontal directions, all have approximately one dominant direction, indicating that they are near-monophase and can be well approximated by their real parts; modes 7–8, which move primarily vertically, are less monophase but still much more so than isolation modes 1–3. Thus, it is fair to assume that the superstructure modes are proportionally damped; as a result, the damp- ing matrix need not be incorporated when updating the superstructure FEM to achieve a mode Table 2.2: Identified natural frequencies and damping ratios (Test 010) (adapted from Brewick et al., 2018) Mode Frequency [Hz] Damping Ratio [%] Isolation        1 0.6853 7.63 2 0.6975 8.62 3 0.7095 7.92 Superstructure                                    4 4.7812 3.21 5 5.1749 3.41 6.1199 76.13 6 7.2931 3.17 10.0301 86.06 7 10.8364 3.50 11.5684 13.19 14.5622 63.46 8 15.3463 3.26 bold indicates fully stable identified mode Table 2.3: Identified effective natural fre- quencies and damping ratios (Test 013) Mode Frequency [Hz] Damping Ratio [%] Isolation        1 0.6079 19.95 2 0.6182 16.65 3 0.7682 25.68 Superstructure                                                                    4 4.6415 4.60 5 5.0268 3.59 5.4095 98.95 6 7.1137 3.40 7.4211 98.99 7.7266 68.17 8.3521 81.80 9.0055 91.56 9.5116 57.01 7 10.6261 4.14 11.3949 32.37 8 11.4646 4.57 13.6537 41.00 14.5471 29.69 9 14.9861 3.54 bold indicates fully stable identified mode 23 Figure 2.8: MAC values of the identified mode shapes against themselves (each natural frequency corresponds to two mode shapes that are complex conjugates of each other) shape match. On the contrary, matching the nonclassically damped isolation modes must update the damping and stiffness simultaneously to achieve a mode shape match. This finding can also be demonstrated by evaluating the MAC values of the identified mode shapes against themselves (each mode has two complex-valued mode shapes that are complex conjugates of each other), as shown in Figure 2.8 (a larger square indicates a larger MAC value and a better match of two mode shapes). When a mode shape is near-monophase (modes 4–8), its conjugate is close to the original mode shape, and their MAC value will be near unity; this is not the case for nonmonophase mode shapes (modes 2–3). 2.2.2 Test 013 A similar methodology for system identification is applied to the experimental data from Test 013; in this case, the input (table) and output (structure) accelerations were additionally low-pass filtered with a 30 Hz low-pass filter (FIR filter with no more than 0.1 dB variation below 30 Hz, and at least two orders of magnitude reduction above 35 Hz) and then downsampled (from 1kHz to 333.3Hz). With model order set to 66, the identified effective frequencies and damping ratios are shown in 24 Table 2.4: Baseline 3DOF model natural frequencies and damping ratios, and the natural fre- quencies of 3DOF models with 20% softer isolation or 20% softer superstructure (∆ denotes the frequency change relative to the baseline; the greater-magnitude light gray entries are those ex- pected to change and are not relevant in this analysis) Baseline Softer Isolation Softer Superstructure Mode # Freq. [Hz] Damping Ratio [%] Freq. [Hz] ∆ ∆ ∆ [%] Freq. [Hz] ∆ ∆ ∆ [%] 1 0.7043 8.13 0.6313 –10.37 0.7024 –0.27 2 4.7232 3.19 4.7136 0.20 4.2355 –10.33 3 9.1976 3.71 9.1968 0.01 8.2274 –10.55 Table 2.3. Comparing Tables 2.2 and 2.3, the natural frequencies of identified modes 4–8 do not change much, indicating at most a very mild nonlinearity in the superstructure during Test 013; indeed, most of the small changes are likely due to the strongly nonlinear isolation-layer devices. 2.3 Validation of Assumptions To validate the assumptions noted in the introduction, this section uses a simple numerical example, with characteristics resembling those of the the experimental building, to show that: (1) the super- structure modes’ natural frequencies and mode shapes are insensitive to isolation-layer stiffness changes; (2) the isolation modes’ natural frequencies and mode shapes are only slightly affected by the superstructure stiffness changes; (3) isolation modes’ damping ratios are mostly unaffected by the superstructure damping changes. The baseline simple model is an isolated building with a two-story superstructure, modeled by a 3DOF shear model with masses, stiffnesses and damping coefficients as shown in Figure 2.9. To simulate the isolation layer, k 0 = 40kN/m is chosen much smaller than k 1 = k 2 = 600kN/m. Table 2.4 shows the natural frequencies and damping ratios of this baseline structure; the first modal frequency 0.7043 Hz is similar to those of the isolation modes of the experimental building; the second modal frequency 4.7232 Hz is similar to that of the first superstructure mode of the 25 m 2 = 0.5 Mg x 2 m 1 = 0.5 Mg x 1 m 0 = 1.0 Mg x 0 Foundation k 0 = 40kN/m c 0 = 1.5Mg/s k 1 = 600kN/m c 1 = 1.0Mg/s k 2 = 600kN/m c 2 = 0.6Mg/s Figure 2.9: Baseline numerical model 0 0.5 1 0 1 2 mode 1 mode 2 mode 3 baseline softer isolation 0 0.5 1 0 1 2 mode 1 mode 2 mode 3 baseline softer superstr. Figure 2.10: Original and softer model mode shape comparisons: (a) baseline model and one with a softer isolation layer; (b) baseline model with one with a softer superstructure experimental building. The damping ratios are similar as well. The mode shapes are demonstrated with the solid lines in Figure 2.10. As shown in Table 2.4, after reducing the isolation-layer stiffness by 20% (the variation bound in the subsequent linear model updating in Section 2.5), the natural frequencies of the superstruc- ture modes (modes 2–3) change by 0.2% or less. Further, the changes in mode shapes, as shown in Figure 2.10a, are almost imperceptible. Thus, the assumption is valid that the superstructure 26 modes’ natural frequencies and mode shapes are essentially unaffected by isolation-layer linear stiffness changes. Similarly, after reducing the superstructure stiffnesses by 20%, the natural fre- quency of the isolation mode (mode 1) changes by less than 0.3%; furthermore, the original and modified mode shapes, as shown in Figure 2.10b, are essentially identical. Thus, it is valid to assume that the isolation modes’ natural frequencies and mode shapes are insensitive to super- structure stiffness changes. Finally, when reducing the superstructure damping coefficents c 2 and c 3 by 40% each, the damping ratio of the isolation mode (mode 1) only reduces by 0.12%, thereby validating the assumption that superstructure damping has negligible effect on the isolation mode damping ratios. 2.4 Introduction to the Finite Element Model The nominal FEM, shown in Figure 2.11, was developed based on the structural design drawings. The beams, columns, and shear walls were modeled by solid concrete elements and embedded reinforcing steel bars were modeled with truss elements; the floor slabs and the nonstructural walls (autoclaved lightweight concrete [ALC] panels) were modeled with shell elements; and the isolation-layer devices were modeled with spring elements. Approximations typically must be made when constructing FEMs, either because some details have little influence on the behavior floor 4 (roof) floor 3 floor 2 floor 1 floor 0 (base) Figure 2.11: Finite element model 27 of the building or because of the vast effort required to include every minute detail; for exam- ple, it is assumed here that each structural member is internally homogenous (constant density, Young’s modulus, etc.), floor slab loading from nonstructural components are ignored, and each nonstructural wall with window(s) is modeled as one homogeneous shell element. To conveniently update the FEM, the global mass matrix M, the nominal global stiffness matrix K 0 , and the global element stiffness matrices K i induced by unit changes to the i th to-be-updated parameterθ i ,i= 1, ...,n θ , are extracted for the model updating analysis. Then, the stiffness matrix K given a specific parameter vector θ θ θ can be written as K= K 0 +∑ n θ i=1 K i (θ i − θ nominal i ), where θ θ θ nominal is the vector of nominal initial parameter values that correspond to nominal initial stiffness K 0 . The first six mode shapes of this nominal FEM are shown in Figure 2.12. (a) Mode 1 (b) Mode 2 (c) Mode 3 (d) Mode 4 (e) Mode 5 (f) Mode 6 Figure 2.12: Mode shapes of the nominal FEM 28 A preliminary study explored two choices for the set of FEM parameters to be updated. In Case I, the to-be-updated parameters θ θ θ included the Young’s moduli of the x- and y-direction beams under floors 0–4 (10 parameters), the columns in each story (4 parameters), the shear walls in each story (4 parameters), the floor slabs on each floor (5 parameters), the stairs in each story 1–3 (3 parameters; the stairs end at floor 3), and the stiffnesses and damping coefficients of the isolation- layer devices in both x and y directions (24 parameters). A simpler Case II groups superstructure elements composed of concrete with the same characteristic compressive strength (Table 2.5), each with a to-be-updated Young’s modulus (3 parameters), and the same isolation-layer coefficients as in Case I (24 parameters). Case I was found to lead to an easier match of modal properties, as expected, because it has a greater number of free parameters, but promotes excessive parameter changes that violate some clear underlying physical properties of the building. Thus, only the updating of the Case II parameters is reported herein. Table 2.5: Nominal compressive strength f ′ ck and Young’s modulusE c , as well as the corresponding post-update Young’s modulusE c , for each superstructure concrete material Nominal Compressive Strength f ′ ck [GPa] Nominal Young’s Modulus E c [GPa] Updated Young’s Modulus E c [GPa] Superstructure Components 36 32.78 34.42 elements in and below floor 0 (base) 27 29.15 30.61 elements above floor 0 and no higher than floor 3 21 25.85 27.14 elements above floor 3 The nominal superstructure concrete Young’s moduli are chosen through interpolation of the relationship between concrete compressive strength f ′ ck and the corresponding Young’s modulus E c provided by the Japanese concrete design code (Ueda, 2007) summarized in Table 2.6. This design code notes that, “in the case of repeated loading or when the applied stress level is low, the values in [Table 2.6] may be preferably increased by 10%” – thus, the concrete Young’s moduli in Table 2.5 are interpolated from the compressive strength in the design drawings using the design code Table 2.6 and then multiplied by 110%. These moduli will be allowed to change by± 5% from their nominal values in the updating optimization. 29 Table 2.6: Relationship between characteristic compressive strength f ′ ck and Young’s modulus E c in the Japanese concrete design code (Ueda, 2007) f ′ ck [GPa] 18 24 30 40 E c [GPa] 22 25 28 31 Although the ALC panels have a Young’s modulus of about 2 GPa, already much smaller than the other concrete elements, the installation allows some in-plane sliding between adjacent panels as well as out-of-plane rotation where connected to the main structure (beams and columns), so the effective stiffnesses of ALC panels are negligible compared with other concrete elements (in fact, a preliminary model updating that included the ALC panels’ effective Young’s modulus resulted in an updated ALC Young’s modulus that was more than three orders of magnitude smaller than those of the concrete elements — again, indicating the ALC plate effective moduli are negligible). Thus, while the panels are included in the modeling as placeholders in the event that future studies on ALC plates or other exterior wall panels — both deemed beyond the scope of this study — are needed, their stiffnesses herein are set to zero and their masses are applied to the adjacent main structural components. In each horizontal direction, each isolation-layer device is modeled in low-level linear motion with a spring with a constant stiffnessk and a viscous damper with a constant damping coefficient c (vertical flexibility is ignored). Thus, the restoring force at time t l is approximated by f(t l )=kx(t l )+c ˙ x(t l ) (2.1) 30 where x(t l ) and ˙ x(t l ) are the displacement and velocity at time t l , respectively, across the device. Combining (2.1) evaluated atN time points,l= 1,2,...,N, into a matrix equation results in                  f(1∆t) f(2∆t) . . . f(N∆t)                  | {z } f =          x(1∆t) ˙ x(1∆t) x(2∆t) ˙ x(2∆t) . . . . . . x(N∆t) ˙ x(N∆t)          | {z } A=[a k a c ]      k c      (2.2) Then, the least squares fit ( ˆ k, ˆ c) for device parameters (k, c) can be calculated using the Moore- Penrose pseudoinverse of matrix A; the resulting values [ ˆ k ˆ c] T = A † f, shown in Table 2.7, are used as the nominal device stiffness and damping coefficient. The restoring force prediction error, or residue, is r= A[ ˆ k ˆ c] T − f. The average relative prediction error err kc =∥r∥ 2 /∥f∥ 2 , where ∥·∥ 2 denotes the vector two-norm, is also shown in Table 2.7. To determine how much the sub- sequent model updating will allow k and c to vary from their nominal initial values, define the residual r ′ = A[ ˆ k+∆k ˆ c+∆c] T − f= r+ a k ∆k+ a c ∆c that results from stiffness∆k and damping ∆c changes from nominal. The bounds on k are chosen to be the tightest range for which∥r ′ ∥ 2 = 1.5∥r∥ when∆c is zero, resulting in∆k={− r T a k ± [(r T a k ) 2 +(1.5 2 − 1)∥a k ∥ 2 2 ∥r∥ 2 2 ] 1/2 }/∥a k ∥ 2 2 ; if these two solutions are denoted∆k + and∆k − , then the lower bound (LB) and upper bound (UB) are max{0, ˆ k+∆k − } (constrained so that the stiffness is never negative) and ˆ k+∆k + , respectively. The corresponding expression for∆c is identical except replacing “k” with “c” throughout. As will be seen in Section 2.6, constraining the isolation-layer coefficients within the resulting bounds, reported in Table 2.7, provides a high fidelity match. 31 Table 2.7: Nominal and updated stiffnesses and damping coefficients of isolation-layer devices RB1 RB2 ESB1 ESB2 SDP1 SDP2 x y x y x y x y x y x y ˆ k [kN/m] 1128 1104 1077 1115 1539 1413 1495 1453 3936 4070 4004 3772 LB 1062 1026 1016 1044 1367 1225 1022 1096 3634 3689 3634 3387 UB 1195 1182 1138 1187 1711 1602 1969 1810 4239 4452 4373 4157 updatedk 1173 1044 1130 1165 1587 1571 1595 1349 3793 3709 3881 3995 ˆ c [Mg/s] 16.74 17.35 15.55 17.84 26.23 29.84 60.06 55.55 101.75 106.65 108.67 119.25 LB 15.88 0.00 1.36 1.17 0.00 0.00 0.00 0.00 30.84 19.15 22.08 30.06 UB 31.90 35.25 29.73 34.51 65.35 73.53 170.20 137.53 172.66 194.15 195.26 208.45 updatedc 24.83 20.47 26.90 24.20 45.42 45.64 82.06 60.85 159.20 134.68 182.22 180.45 err kc [%] 5.26 6.29 5.05 5.74 9.92 11.79 26.87 21.19 6.81 8.30 8.16 9.00 2.5 Linear Updating and Mode Matching Most model updating procedures result in optimization problems that minimize some cost (or objective) functions. Here, the error metric J(θ θ θ) is defined to be the sum of the error of each identified mode with the “matched” corresponding FEM mode: J(θ θ θ)= ∑ i∈I ID J i ′ i i (θ θ θ) (2.3a) whereI ID is a set of experimentally identified mode numbers used in the optimization, i ′ i is the FEM mode number matched to experimentally identified mode i, andJ i ′ i i (θ θ θ) is a weighted sum of the relative differences between FEM modei ′ i and identified mode i: J i ′ i i (θ θ θ)= ˆ f i ′ i (θ θ θ)− f i /f i +λ MAC i 1− MAC ˆ φ φ φ i ′ i (θ θ θ),φ φ φ i (2.3b) +λ ζ i ˆ ζ i ′ i (θ θ θ)− ζ i /ζ i where f i ,φ φ φ i andζ i are thei th experimentally identified natural frequency, mode shape and damping ratio, whereas ˆ f i ′ i (θ θ θ), ˆ φ φ φ i ′ i (θ θ θ) and ˆ ζ i ′ i (θ θ θ) are the corresponding values for the FEM’si ′ i th mode. The FEM mode numberi ′ i may be the same as the identified mode number i, but this is often not the case 32 because of FEM modes that remain unidentified from the experimental data (further discussed in the following section). This cost function uses the relative frequency and damping ratio differences so that all modes are evaluated on the same scale. 1− MAC( ˆ φ φ φ i ′ i (θ θ θ),φ φ φ i ) quantifies the differences between the identified mode shape φ φ φ i and the FEM mode shape ˆ φ φ φ i ′ i (θ θ θ). λ MAC i and λ ζ i are weights that adjust the contribution to cost function (2.3) by errors from mode shapes and damping ratios, relative to those from the frequencies. Reducing relative fre- quency error and mode shape error are considered approximately of equal importance to ensure reproduction of structure time-history responses, so the mode shape weightλ MAC i is chosen to be unity for each mode i. The superstructure modes’ damping ratios, which are sufficiently small to generate primarily monophase mode shapes, are of relatively little utility for matching those modes’ frequencies and mode shapes; thus, the damping relative error weightλ ζ i is set to zero for superstructure modes i≥ 4. In contrast, because the isolation mode damping is much larger, with nonproportional mode shapes, accurate damping modeling in the isolation modes is essential for accurate response prediction, soλ ζ i = 1 for modesi= 1,2,3. Finally, the superstructure modes are first updated using I ID str ={4,5,...,8}, and then the isolation modes usingI ID iso ={1,2,3}. 2.5.1 Mode Matching With the cost functions defined, the next challenge is to find the FEM mode i ′ i that corresponds to identified mode i. In other words, although there are 8 identified stable modes in Table 2.2, one must determine, for each identified mode, the FEM mode that corresponds: i.e., similar frequency and mode shape. For example, stable identified mode 8 is found in Section 2.6 to not correspond to FEM mode 8, but rather FEM mode 9. This mismatch is expected: due to the limited number of measurement directions and/or locations, it is possible (perhaps even likely) that some modes of the experimental building cannot be identified; further, sensor noise can result in spurious non- physical modes appearing in the identification results. When updating the parameters of the isolation-layer devices, mode matching is typically not a challenge because the number of identified modes and the number of FEM modes that can possibly 33 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 i= 1 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 i= 1 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 i= 1 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 i= 1 i= 1 i= 2 i= 3 i ′ = 1 i ′ = 2 i ′ = 3 Figure 2.13: Possible mode matchings for the isolation-layer device parameter updating be matched to them are both very small (three each), so there are only six possible match permu- tations, as shown in Figure 2.13. A cost function can be computed for each of the six scenarios, the smallest of which not only provides the updated parametersθ θ θ, but also indicates the best mode match. In contrast, in the general case of n ID identified modes, each to be matched to one of n FEM modes (which likely does not constitute all FEM modes but rather those within a frequency range reasonably close to those of the identified modes), the number of permutations is n FEM !/(n FEM − n ID )!, which is typically far too large for a tractable exhaustive search. For a particular choice of parameter vector θ θ θ, a modal analysis will provide a set of n FEM candidate modal frequencies ˆ f j (θ θ θ) and corresponding mode shapes ˆ φ φ φ j (θ θ θ) for j∈I FEM ; for identifying the superstructure modes of the example herein,I FEM ={j| ˆ f j (θ θ θ)∈[2,17]Hz} because the modes below 2 Hz are the isolation modes, and the experimentally-identified modes only range up to about 17 Hz. Let U ID ⊂I ID andU FEM ⊂I FEM be the sets of yet unmatched experimentally-identified and FEM modes, respectively; when the mode-matching algorithm begins, these are identically I ID and 34 I FEM , respectively. Assume that the modes are numbered by increasing frequency; i.e., f i < f i+1 and ˆ f j (θ θ θ)< ˆ f j+1 (θ θ θ) for any positivei or j. FEM θ θ θ Modal Analysis i= arg min m∈U ID f m l = arg min m∈U FEM ˆ f m (θ θ θ) i ′ i = 0, J 0i =∞ ˆ f l (θ θ θ) ? ∈ f i (1± r f ) and 1− MAC( ˆ φ φ φ l (θ θ θ),φ φ φ i ) ? <r MAC J li (θ θ θ) = |( ˆ f l (θ θ θ)/f i )− 1|+ 1− MAC( ˆ φ φ φ l (θ θ θ),φ φ φ i ) J li (θ θ θ) ? <J i ′ i i (θ θ θ) i ′ i =l l ? = arg max m∈U FEM ˆ f m (θ θ θ) l = arg min m∈U FEM ˆ f m (θ θ θ)> ˆ f l (θ θ θ) ˆ f m (θ θ θ) Remove i fromU ID . Remove i ′ i (or all{ j| j∈U FEM , j max i∈I ID ,i/ ∈U ID ˆ f i ′ i (θ θ θ)}. 2.5.2 Optimization To minimize the error metric J(θ θ θ), two optimization algorithms were evaluated: the Genetic Al- gorithm (GA) in the MATLAB ® Global Optimization Toolbox and a Nelder Mead Simplex method via thefminsearch function in MATLAB ® (a traditional hill climbing optimizer). For this opti- mization, the GA uses a population of 400 for 40 generations, with defaults for other parameters (5% elite rate, 80% crossover fraction, 1% mutation rate). When updating the superstructure pa- rameters, no further tuning of the algorithm was needed. However, updating the isolation-layer parameters by directly minimizingJ(θ θ θ) in (2.3) results in no value ofθ θ θ that is within the permis- sible range[θ θ θ LB ,θ θ θ UB ] and also provides a “close” match, defined as satisfying all of the following ˆ f i ′ i (θ θ θ)− f i /f i < 5% (2.6a) ˆ ζ i ′ i (θ θ θ)− ζ i /ζ i < 5% (2.6b) 1− MAC ˆ φ φ φ i ′ i (θ θ θ),φ φ φ i < 5% (2.6c) for every i∈I ID iso (if λ ζ i or λ MAC i were zero, which they are not herein for the isolation modes, the corresponding criterion would be ignored for a “close” match). To overcome this difficulty, 37 the first 20 GA generations are devoted to finding a permissible θ θ θ population that gives a “close” match; this is accomplished by instead minimizing a somewhat relaxed cost function ˜ J(θ θ θ)= ∑ i∈I ID iso max n 0, 0.05− ˆ f i ′ i (θ θ θ)− f i /f i o + ∑ i∈I ID iso λ MAC i max n 0, 0.05− h 1− MAC ˆ φ φ φ i ′ i (θ θ θ),φ φ φ i io + ∑ i∈I ID iso λ ζ i max n 0, 0.05− ˆ ζ i ′ i (θ θ θ)− ζ i /ζ i o (2.7) that eventually converged to ˜ J(θ θ θ)= 0, indicating that the “close” criteria (2.6) were met. The remaining 20 generations, starting from the 20 th generation’s population, minimized J(θ θ θ) to seek an even closer match. The hill-climber failed to outperform the GA in the following respects. (a) As is typical, this hill-climber tended to get stuck at local minima and failed to search some key parts of the pa- rameter space. (b) Hill-climbing optimizers, in general, offer only one solution to the problem; however, the GA usually produces a population of solutions with extremely similar cost function values. In preliminary studies, one (or more) of the competing GA solutions outperformed the only solution obtained through the hill-climbing optimizer. (c) The GA is gradient-free and can mini- mize discontinuous cost functions. In this study, although cost function (2.3) seems to be smooth, it is indeed piecewise discontinuous due to the incorporation of the mode-matching decision pro- cess (Figure 2.13 or 2.14). Thus, J(θ θ θ) does not have finite gradients everywhere, which poses a problem for most hill-climbers. Thus, while preliminary studies evaluated the fminsearch hill-climber, it is subsequently omitted due to its inferior performance relative to the GA. 2.6 Results of Linear Model Updating Figure 2.15 shows comparisons of the natural frequencies and mode shapes between the experi- mental building and the FEM before and after linear model updating; the natural frequencies of the 38 identified and FEM modes are shown on the horizontal and vertical axes, respectively, as well as listed in Table 2.8. For most modes, the natural frequency errors have been reduced, especially for those with initially relatively large errors: the relative frequency differences for modes 4–6 reduce from 18%–29% to less than 1.5%. Further, after model updating, the identified mode shapes of modes 1–3 and mode 8 achieve much better matches with the corresponding FEM mode shapes. Since modes 1–3 (either identified or FEM) are the dominant modes for the horizontal motions, it should be expected that the dynamic behavior of the updated FEM in the horizontal directions matches much more closely with those of the experimental building. The effective masses (and moments of inertia) of the updated FEM’s modes, computed by means of modal analysis, are listed in Table 2.9, whereM x i ,M y i andM z i are the effective masses of the i th mode in the x, y, and z directions, respectively, and I x i , I y i and I z i are the effective moments of inertia of thei th mode about thex,y, andz axes, respectively. The total massM total and the total moments of inertia of the FEM are readily available after constructing the FEM. For the x- and y-direction translations, as well as the rotation around the z axis, the sum of the effective masses (or moments of inertia) from the 8 listed modes are almost 100% of the total. For the z-direction translation, the FEM mode corresponding to the identified mode 7 has an effective mass as large 0.676 0.688 0.703 5.626 6.666 9.149 10.728 15.504 15.763 0.685 0.698 0.710 4.781 5.175 7.293 10.836 15.346 0 0.2 0.4 0.6 0.8 1 MAC (a) Before linear model updating 0.688 0.698 0.706 4.767 5.239 7.297 10.724 12.291 15.247 0.685 0.698 0.710 4.781 5.175 7.293 10.836 15.346 0 0.2 0.4 0.6 0.8 1 MAC (b) After linear model updating Figure 2.15: Comparisons of frequencies and mode shapes (MAC values) between experimentally identified modes and those from the (a) unupdated and (b) updated FEM models 39 Table 2.8: Frequency comparisons of the experimental building, and the original and updated FEMs N4SID Identified Original FEM Updated FEM Mode # Freq. [Hz] Mode # Freq. [Hz] (Err. [%]) Mode # Freq. [Hz] (Err. [%]) 1 0.685 1 0.676 (–1.352) 1 0.686 ( 0.174) 2 0.698 2 0.688 (–1.355) 2 0.700 (–0.376) 3 0.710 3 0.703 (–0.912) 3 0.702 (–1.051) 4 4.781 4 5.626 (17.668) 4 4.767 (–0.292) 5 5.175 5 6.666 (28.813) 5 5.240 ( 1.248) 6 7.293 6 9.149 (25.447) 6 7.296 ( 0.041) 7 10.836 7 10.728 (–1.000) 7 10.724 (–1.041) 8 15.504 n/a 8 12.292 n/a 8 15.346 9 15.763 ( 2.540) 9 15.247 (–0.645) as 87% of the total mass, making it the dominant mode for vertical motion. Further, the moments of inertia around the x and y axes, together, exceed 0.85 and 0.90, respectively. In sum, the FEM modes that are matched with the identified modes can capture the dominant dynamic behaviors of the experimental building. The (updated) isolation-layer device damping coefficients provide some damping in both iso- lation and superstructure modes. The three isolation-mode damping ratios are are now 0.0761, 0.0842 and 0.0788, which are only 0.29%, 2.34% and 0.45%, respectively, below those identi- fied and reported in Table 2.2; since these errors are quite small, no additional modal damping is needed for isolation modes 1–3. However, the superstructure-mode damping ratios induced by the isolation-layer device linear viscous damping coefficients are quite small: 0.0079, 0.0066 and 0.0037 in modes 4–6 and effectively zero in all other superstructure modes. Thus, additional modal damping in modes 4–6 is added so that their total damping ratios match those reported in Ta- ble 2.2;i.e., 0.0321− 0.0079= 0.0242, 0.0341− 0.0066= 0.0275 and 0.0317− 0.0037= 0.0280, respectively. For mode 7, the updated model uses damping ratio 0.0414, taken from the Test 013 vertical mode identification in Table 2.3, which was found to well reproduce (and better than using Test 010’s 0.0350 from Table 2.2) vertical motion in both the low-level random excitation and the 40 Table 2.9: Effective mass and moment of inertia ratios of the FEM modes matched to identified modes 1–8 Identified Mode # (i) Matched FEM Mode # (i ′ i ) M M M x x x iii M M M total M M M y y y iii M M M total M M M z z z iii M M M total I I I x x x iii I I I x x x I I I y y y iii I I I y y y I I I z z z iii I I I z z z 1 1 0.097 0.744 0.000 0.067 0.039 0.123 2 2 0.902 0.082 0.000 0.007 0.315 0.725 3 3 0.000 0.173 0.000 0.013 0.000 0.152 4 4 0.000 0.000 0.000 0.008 0.337 0.001 5 5 0.000 0.000 0.011 0.152 0.076 0.000 6 6 0.000 0.000 0.002 0.010 0.015 0.000 7 7 0.000 0.000 0.871 0.562 0.135 0.000 8 9 0.000 0.000 0.004 0.033 0.011 0.000 ∑ 8 i=1 : 1.000 1.000 0.889 0.854 0.928 1.000 stronger earthquake excitations. For modes 8 and up, the damping ratios are set to 0.0331, which is the average identified Test 010 superstructure mode damping ratio ( i.e., average of the damping ratios of modes 4–8 in Table 2.2). To evaluate the effect of model updating on the time-history responses, the updated FEM is subjected to the random excitation in Test 010. While various implicit or explicit time integration algorithms could be used, a modal approach is used here, simulating the response of each of the first 36 modes (those less than 35 Hz) to a piecewise linear form of the 1 kHz-sampled table- center acceleration ¨ x g (t); the modal responses are then combined with the mode shapes to provide the global response. For a measured experimental response w(t) and the corresponding FEM prediction ˆ w(t;θ θ θ), one metric of time-history relative prediction error is Err RMS w (θ θ θ)= |RMS[ ˆ w(t;θ θ θ)]− RMS[w(t)]| RMS[w(t)] × 100% (2.8) where RMS[w(t)]=[t − 1 f R t f 0 w 2 (t)dt] 1/2 , or discrete time[n − 1 t ∑ n t − 1 k=0 w 2 (k∆t)] 1/2 , is the root-mean- square response (RMS) over time[0,t f ]. Table 2.10 shows that the relative errors — using error metric (2.8) — in RMS acceleration responses in the x and y directions drop significantly, by about two thirds, after updating: from 41 Table 2.10: Percentage errors, using metric (2.8), in acceleration response and isolation-layer de- vice force RMSs using the linear FEM before (after) updating Direc. Status Errors* [%] in Acceleration RMS Errors [%] in Isolation-Layer Device Force RMS RB1 RB2 ESB1 ESB2 SDP1 SDP2 x before 7.3 – 34.5 13.8 31.0 13.3 26.7 16.9 16.7 x (after) (0.2 – 12.9) (2.9) (6.1) (6.4) (1.4) (10.5) (9.8) y before 9.2 – 26.9 9.9 14.8 14.1 7.8 9.8 14.5 y (after) (0.3 – 8.0) (4.4) (2.3) (8.2) (7.8) (7.7) (4.0) z before 0.0 – 8.5 n/a n/a n/a n/a n/a n/a z (after) (1.0 – 7.5) (n/a) (n/a) (n/a) (n/a) (n/a) (n/a) * minimum and maximum errors across the 14 structure accelerometers as much as 34.5% to at most 12.9% in the x direction, and from 26.9% to at most 8.0% in the y direction; the z-direction change is not as significant, but the vertical accelerations from the unupdated FEM were already more accurate than those in the x and y directions. Further, most of the relative errors in RMS forces in the x and y directions drop significantly after updating, with the largest among them also dropping by more than two thirds. FEM Experiment (a)x-direction acceleration FEM Experiment (b) Zoomed-in version Figure 2.16: Measured and calibrated-FEM-predicted accelerations, both downsampled to 100 Hz, from the accelerometer near grid 1A at floor 0 (base) in random excitation Test 010: (a) full time duration; (b) narrower time window to show detailed comparison The acceleration response with the largest relative error in RMS, thex-direction acceleration at the base-level sensor near grid 1A on Figure 2.2, is shown in Figure 2.16, both as measured (and downsampled to 100 Hz) and as predicted with the updated FEM (similar graphs for the responses 42 FEM Experiment Figure 2.17: Measured force-displacement loops, and those predicted by the entire updated linear FEM, of the SDP1 in thex direction in random excitation Test 010 with smaller prediction errors, omitted for brevity, exhibit FEM predicted accelerations that look almost exactly the same as those measured). Even for this worst case, the acceleration time history from the experiment is still fairly well reproduced. Figure 2.17 shows a comparison of the simulated and the measured force-displacement loops for the SDP1 in the x direction; even though this force exhibits the largest prediction error among all restoring forces, visual inspection shows that the equivalent linear stiffness and damping of the measured restoring force are sufficiently represented by the linear simulation. 2.7 Formulation of the Partially Nonlinear FEM The rubber bearings exhibited slight stiffness reductions during large earthquake-induced displace- ments, so linear models are close but not perfect for modeling their behavior; a directionally un- coupled bilinear model is utilized herein. Models of the two elastic sliding bearings and the two metallic-yielding steel damper pairs were developed in a previous study (Brewick et al., 2020) based on the responses during Test 016 (Tohoku) and are mildly suboptimal during random exci- tations, primarily due to insufficient damping. Thus, slightly different ESB and SDP models are developed in this section. Then, the methodology for coupling the nonlinear isolator models with the updated linear superstructure model is provided, followed by an evaluation of the time-history predictions of the combined building model. 43 Figure 2.18: RB1 x-direction force- displacement loops in Test 016 ¯ u RBi − ¯ u RBi force drift 1 k RBi 1 ¯ k RBi Figure 2.19: Rubber bearing bilinear elastic stiffness model 2.7.1 Modeling the Isolation-Layer Devices As exemplified in Figure 2.6, the RBs behaved linearly in each of the x andy directions under ran- dom excitations (Tests 010–012); linear regression of the force-displacement data estimates stiff- nesses of approximately 1100 kN/m. During the earthquake excitations (Tests 013–016), the RBs also appear quite linear, but a linear regression of Test 016 force-displacement data (Figure 2.18) estimates stiffnesses of approximately 950 kN/m, indicating some stiffness reduction when the dis- placement is large. Thus, RBi (i∈{1,2}) is modeled herein with two uncoupled identical bilinear elastic stiffness elements — as depicted in Figure 2.19 with initial stiffnesses k RBi =k RBi xx =k RBi yy , postyield stiffnesses ¯ k RBi = ¯ k RBi xx = ¯ k RBi yy , and yield displacement ¯ u RBi = ¯ u RBi x = ¯ u RBi y — in paral- lel with the linear viscous damping coefficients c RBi xx and c RBi yy found in the previous section and reported in Table 2.7. Thus, RBi’s force in thex direction is given by ˆ f RBi x =k RBi u RBi x +c RBi xx ˙ u RBi x +(k RBi − ¯ k RBi )(u RBi x + ¯ u RBi )[1− H(u RBi x + ¯ u RBi )] +(k RBi − ¯ k RBi )(u RBi x − ¯ u RBi )H(u RBi x − ¯ u RBi ) (2.9) whereH(·) is the Heaviside unit step function; the formula for they-direction force is identical by replacing everyx withy. 44 The ESBs and SDPs, which exhibit clear hysteretic behavior in all of the tests, are each modeled using a bidirectionally coupled Bouc-Wen hysteretic behavior of the form (Park et al., 1986) ˙ Z x =A ˙ u x − β| ˙ u x Z x |Z x − γ ˙ u x Z 2 x − β ˙ u y Z y Z x − γ ˙ u y Z x Z y (2.10a) ˙ Z y =A ˙ u y − β ˙ u y Z y Z y − γ ˙ u y Z 2 y − β| ˙ u x Z x |Z y − γ ˙ u x Z x Z y (2.10b) where u x and u y are the relative displacements across the device in the x and y directions, respec- tively, and A, β and γ control the shape of the hysteresis. The two SDPs have nearly identical behavior but the two ESBs act somewhat differently from each other (Brewick et al., 2020); thus, the SDPs are assumed to have identical parameters, but the ESB1 and ESB2 parameters are allowed to differ. The resulting ESB and SDP model forces are then provided by      ˆ f ESBi x ˆ f ESBi y      =µ ESBi (t)W ESBi (t)      Z ESBi x Z ESBi y      +      c ESBi xx ˙ u ESBi x c ESBi yy ˙ u ESBi y      (2.11a)      ˆ f SDPi x ˆ f SDPi y      = K SDP    α      u SDPi x u SDPi y      +(1− α)      Z SDPi x Z SDPi y         +      c SDPi xx ˙ u SDPi x c SDPi yy ˙ u SDPi y      (2.11b) µ ESBi (t)=µ 0 [S − 1 W ESBi (t)· (mm 2 /N)] µ 1 [V ESBi (t)· (s/mm)] µ 2 is the coefficient of friction, where S= 0.1018m 2 is the slider’s effective surface area,W ESBi (t) is the weight carried by ESBi at time t, V ESBi (t)=∥[ ˙ u ESBi x (t) ˙ u ESBi y (t)]∥ 2 is the sliding velocity, and the friction formula coefficients areµ 0 = 0.0606,µ 1 =− 0.5883 andµ 2 = 0.1148 as determined in Brewick et al. (2020); K SDP = [k SDP xx k SDP xy ] T [k SDP xy k SDP yy ] T is a symmetric SDP stiffness matrix;α defines the SDP postyield-to- preyield stiffness ratio; andc (·) xx andc (·) yy are damping coefficients in x andy directions, respectively (which are not fixed at the linear model’s values in Table 2.7 but rather subsequently optimized). Several metrics M (·) can be used to evaluate the fidelity of a single isolation-layer device model or of a set of device models. MetricM ˆ f(t,u(t)) : The model error for a single device can be quantified by the mean square error between the measured device forces’ time histories and 45 those predicted using the device model in (2.9) or (2.10)–(2.11) driven by the measured cross- device drifts and drift velocities; this metric is defined mathematically in (A.1) in Appendix A. MetricM ˆ f(t;θ θ θ,¨ x g ) : The force response errors for a candidate set of models for the devices are com- puted by simulating the response of a partially nonlinear FEM (linear FEM coupled with candidate models of the isolation-layer devices) to the measured shake table accelerations, extracting the resulting device force time histories, and quantifying the relative errors in the forces’ RMS using (2.8). Metric M ˆ ζ : Modal damping evaluation is determined for a candidate set of models for the devices by comparing the measured responses’ identified modal damping in Section 2.2 with the effective model damping estimated, using the same system identification strategy, from the 42 acceleration time-history responses of a partially nonlinear FEM at the 14 accelerometer locations. Preliminary studies showed that omitting the linear viscous terms in (2.11) leads to isolation- mode damping ratios computed from time-history responses to low-level random excitations such as Test 010 (M ˆ ζ ) that are up to about a quarter smaller than those identified and reported in Table 2.2. Further, straightforward optimization of the ESB and SDP parameters in (2.10)–(2.11) to match the experimental force-displacement relationships in Tests 010 and 016 (i.e., minimizing M ˆ f(t,u(t)) ) still does not lead to sufficient levels of damping ( M ˆ ζ ). Finally, the RB parameters in (2.9) that minimize the mean square RB force prediction error (M ˆ f(t,u(t)) ) do not consistently lead to smaller errors in simulated RB force RMS (M ˆ f(t;θ θ θ,¨ x g ) ); rather, slightly different RB parameters that are very mildly suboptimal in reducing mean square RB force error (M ˆ f(t,u(t)) ) can provide smaller simulated RB force RMS (M ˆ f(t;θ θ θ,¨ x g ) ). To address these trade offs, a two-level hierarchical optimization strategy — detailed in Appendix A — is used to introduce sufficient damping in low-level response and balance trade offs between local device-specific model fidelity and that of the overall partially nonlinear FEM. The resulting isolation-layer device parameters are listed in Table 2.11. It may be noted that the optimal parameters for the two ESBs are quite different from each other, and the damping coefficients for ESB1 in the two horizontal directions differ significantly as well. This is not surprising given that ESB1 was seen to exhibit a strong stick-slip behavior that 46 Table 2.11: Parameters in the RB, ESB and SDP models Parameter RB1 RB2 ESB1 ESB2 SDP k xx 1.075 MN/m 1.025 MN/m 3.975 MN/m k yy 1.075 MN/m 1.025 MN/m 3.779 MN/m k xy 0 0 0 ¯ k 0.786 MN/m 0.772 MN/m ¯ u 0.0258 m 0.0309 m c xx 24.83 Mg/s 26.90 Mg/s 3.79 Mg/s 30.00 Mg/s 90.00 Mg/s c yy 20.47 Mg/s 24.20 Mg/s 50.04 Mg/s 20.00 Mg/s 90.00 Mg/s A 47.13m –1 54.10m –1 1 β 10.41 m –1 25.41 m –1 635 m –2 γ 36.72 m –1 28.69 m –1 1102 m –2 α 0.1225 Bold indicates optimized parameters; italic indicates parameters computed from opti- mized parameters; other parameters are fixed. was quite different from ESB2 Brewick et al. (2020). Further, the SDP Bouc-Wen parameters here differ mildly from those reported in Brewick et al. (2020), which optimized solely for local device behavior in Test 016; the parameters here generalize better to the random excitation response in Test 010 as well as the lower-level earthquake response in Test 013 (which exhibited larger postyield SDP stiffnesses, possibly because of smaller motion than in Test 016 or because the SDP stiffnesses degraded in Tests 013–015). 2.7.2 Time-History Analysis The linear superstructure is combined with the nonlinear isolation-layer device models to simulate the time-history responses of the building; the resulting equation of motion is: M¨ x+ C˙ x+ Kx=− MP¨ x g + K nom iso x iso − B iso f iso (x iso , ˙ x iso ) (2.12) where M, C and K are the n× n (n is the number of DOFs) global mass, damping and stiffness matrices of the updated linear FEM, in which M remains unchanged throughout this study, C represents the damping in the superstructure (computed from the added modal damping discussed 47 in Section 2.6) and the linear viscous damping in the isolation layer, and K is the updated linear stiffness matrix; P is ann× 6 influence matrix for the 6 × 1 ground (shake table center) excitation vector ¨ x g ; x is a vector of nodal displacements relative to the shake table; x iso = G iso x is a vector of isolation-layer device drifts relative to the shake table (G iso is a transformation matrix extracting these drifts from the vector of all displacements); K nom iso is a stiffness matrix that cancels out the nominal linear isolation-layer restoring forces (from the updated linear FEM) so that the actual nonlinear restoring forces can be included as− B iso f iso (x iso , ˙ x iso ), where f iso (x iso , ˙ x iso ) is a vector of isolation-layer device restoring forces that depend on x iso and ˙ x iso and B iso is an influence matrix that reintroduces f iso (x iso , ˙ x iso ) back into the global force vector. The advantages of this formulation are three-folded. First, the left-hand side of (2.12) re- mains linear, so the modal superposition method remains applicable. Since the experimental data is low-pass filtered at 35 Hz, modes with frequencies above 35 Hz can be neglected, reducing the 84,669 DOFs down to 36 modal DOFs. Second, regardless of the excitation, the eigenvalues and eigenvectors of the dynamic system remain unchanged, so modal decomposition need be per- formed only once. Third, since the restoring forces of the rubber bearings, elastic sliding bearings and the U-shape steel damper pairs are introduced as external forces in this formulation, differ- ent models of the isolation-layer devices can be easily incorporated into this formulation only by changing the form of f iso (x iso , ˙ x iso ) on the right-hand side, instead of making extensive changes to the structural system (the left-hand side), which is a benefit in the subsequent studies. Utilizing (2.12), the partially nonlinear model’s time-history responses to the earthquake ex- citation from Test 013 are computed and compared with the experimentally measured responses (note that Test 013 is intentionally chosen here as it was not used in calibrating the isolation-layer devices); relative errors in RMS are shown in Table 2.12. The z-direction accelerations are very well reproduced, indicating that the superstructure FEM is well constructed and updated; the small discrepancies are likely due to approximations in constructing the FEM. The slightly lower fidelity of the x- and y-direction predicted accelerations could result from the limitations of the bidirec- tional Bouc-Wen model of the isolation layer that cannot model the stick-slip behavior of the ESB1 48 Table 2.12: Percentage errors, using metric (2.8), in acceleration response and isolation-layer de- vice force RMSs using the partially nonlinear system model Direction Errors* [%] in Acceleration RMS Errors [%] in Isolation-Layer Device Force RMS RB1 RB2 ESB1 ESB2 SDP1 SDP2 x 1.4 – 11.4 11.2 5.4 23.8 9.6 2.1 0.2 y 3.6 – 9.7 9.9 8.0 7.3 11.7 2.1 5.5 z 2.5 – 6.8 n/a n/a n/a n/a n/a n/a * minimum and maximum errors across the 14 structure accelerometers (especially in thex direction) observed in this study Brewick et al. (2020), as exhibited by the larger error in RMS x-direction ESB1 restoring forces. The unmodeled stick-slip behavior also affects the y-direction restoring forces due to the rotation of the building. The rubber bearing restoring forces are better reproduced than those of the ESBs due to accurate estimation of the stiffnesses; the SDP restoring forces are reproduced very accurately, with errors in RMS at most a few percent. Figure 2.20 shows comparisons of the accelerations on floor 0 at the sensor near grid 1A on Figure 2.2: thex-direction response measured by this sensor had the largest relative errors in RMS prediction (11.4% as shown in Table 2.12). In the x and y directions, the large responses (e.g., around 50 and 100 sec) are reproduced well, with a very good match in terms of magnitudes and frequencies. Visual inspection shows that the y-direction response is reproduced slightly better than the x-direction response. The z-direction response is simulated very well, which is typical for all responses in this direction, and the slight differences may be due to small differences in damping for modes with vertical motion. Further, the good matches of accelerations in all three directions indicate that a linear superstructure model is sufficient, which is consistent with the intent of typical base-isolation design. As observed in Figure 2.21, the stiffnesses of the RB bilinear model match those from the experimental tests in both directions. For the elastic sliding bearings and the steel damper pairs, the restoring forces generated by the modified Bouc-Wen models match well the experimental data in terms of slopes and outer loop excursions (perfect matches cannot be obtained due to the 49 FEM Experiment FEM Experiment FEM Experiment FEM Experiment FEM Experiment FEM Experiment Figure 2.20: Measured and calibrated-FEM-predicted accelerations from the accelerometer near grid 1A on floor 0 in Test 013 (scaled Tohoku-Oki earthquake): the left graphs show the full time duration; the right graphs show a narrower time window to show detailed comparison (note: the right z-direction graph uses a narrower time window to show the detail for the higher-frequency vertical motion) irregular stick-slip behavior of ESB1, which poses greater modeling challenges beyond the scope of this study). 50 FEM Experiment FEM Experiment FEM Experiment FEM Experiment FEM Experiment FEM Experiment Figure 2.21: Measured and calibrated-FEM-predicted force-displacement loops during the scaled Tohoku-Oki Test 013 51 Chapter 3 Cross-Model Cross-Mode Model Updating Method Based on Constrained Total Least Squares 3.1 Brief introduction to the CMCM Method 3.1.1 CMCM Method Basics The equations of motion of an unforced linear structural dynamical system can be written as M¨ q+ C˙ q+ Kq= 0 (3.1) where M, C and K are (unknown) n× n mass, damping and stiffness matrices; and ¨ q, ˙ q and q are n-element acceleration, velocity and displacement vectors. An eigen-analysis of (3.1) leads to (λ 2 i M+λ i C+ K)φ φ φ i = 0 (3.2) where λ i =− ω i (ζ i ± j q 1− ζ 2 i ) is the eigenvalue and φ φ φ i is the eigenvector of the i th mode of the damped system, where j= √ − 1, and ω i and ζ i are the natural frequency and damping ratio, respectively, of thei th mode (in practice,λ i andφ φ φ i are obtained by performing system identification based on the inputs and outputs of the structural system, and are uncertain because the inputs and outputs are inevitably polluted by noise). In practice, the matrices M, C and K of the experimental 52 structure are unknown, so they are often expanded as the sum of the corresponding unupdated FEMn× n matrices ˆ M, ˆ C and ˆ K and a series of (typically low-rank) changes of to-be-determined magnitudes: M= ˆ M+ ∑ h µ h ˆ M h , C= ˆ C+ ∑ h χ h ˆ C h , K= ˆ K+ ∑ h κ h ˆ K h (3.3) with magnitudes quantified by the coefficients µ h , χ h andκ h (note that the upper limits on h need not be the same for mass, damping and stiffness). Let ˆ λ i and ˆ φ φ φ i denote thei th eigenvalue and eigenvector of the unupdated FEM given by ˆ λ 2 i ˆ M ˆ φ φ φ i + ˆ λ i ˆ C ˆ φ φ φ i + ˆ K ˆ φ φ φ i = 0. (3.4) Then, premultiplying (3.2) by ˆ φ φ φ H j (the Hermitian transpose * of ˆ φ φ φ j ) yields ˆ φ φ φ H j λ 2 i Mφ φ φ i + ˆ φ φ φ H j λ i Cφ φ φ i + ˆ φ φ φ H j Kφ φ φ i = 0. (3.5) Expanding (3.5) with the matrix series in (3.3) and rearranging terms gives a set of equations, linear in the unknown magnitudes, of the form ∑ h µ h λ 2 i m h,ij + ∑ h χ h λ i c h,ij + ∑ h κ h k h,ij =− λ 2 i m ij − λ i c ij − k ij , i= 1,...,n exp , j= 1,...,n FEM (3.6) where n exp and n FEM are the number of modes used from the experiment and FEM, respectively, and m ij = ˆ φ φ φ H j ˆ Mφ φ φ i , c ij = ˆ φ φ φ H j ˆ Cφ φ φ i , k ij = ˆ φ φ φ H j ˆ Kφ φ φ i , m h,ij = ˆ φ φ φ H j ˆ M h φ φ φ i , c h,ij = ˆ φ φ φ H j ˆ C h φ φ φ i , k h,ij = ˆ φ φ φ H j ˆ K h φ φ φ i . (3.7) In some cases, it may be that mass is known exactly or with sufficient accuracy that the µ h are small enough that the corresponding terms can be neglected; (3.6) simplifies to ∑ h χ h λ i c h,ij + ∑ h κ h k h,ij =− λ 2 i m ij − λ i c ij − k ij . (3.8) * Herein,(·) H denotes the conjugate transpose,(·) T denotes the transpose and(·) ∗ denotes the conjugate. 53 If only the stiffness matrix is to be updated (which is assumed for the 3D frame numerical example herein), then (3.8) can be further reduced to ∑ h γ h k h,ij =− λ 2 i m ij − λ i c ij − k ij . (3.9) Eqs.(3.6), (3.8) and (3.9) can each be transformed into a matrix equation Aˆ x= b (3.10) where ˆ x is [µ T χ T κ T ] T or [χ T κ T ] T or [κ], and both A and b have (n exp n FEM ) rows. It should be noted that matrix A and vector b are complex-valued; for the TLS approach, (3.10) can be converted into a corresponding real-valued equation    ℜ(A) ℑ(A)    ˆ x=      ℜ(b) ℑ(b)      (3.11) whereℜ andℑ are the real and imaginary operands, respectively. The TLS solution for (3.11) is well studied; for the CTLS approach, however, matrix manipulations detailed in Section 3.2.3 are required so that the solution can be found efficiently through quadratic programming. 3.1.2 Model Expansion for Spatially Incomplete Experimental Mode Shapes The standard CMCM uses full mode shapes, including all DOFs, for both the experimental and the FEM mode shapes; i.e., the dimension of φ φ φ i is assumed the same as of ˆ φ φ φ i . However, given the limited number of sensors in a practical experiment, it is rarely possible, especially outside the laboratory, to estimate the experimental mode shapes at all DOFs. Thus, a model expansion technique is required. Several such techniques have been adopted in previous studies, including the optimal filtering method (Wang et al., 2015) and the Guyan method (Li et al., 2020). However, the dynamic model expansion technique is chosen herein because, compared with these other methods, 54 it is known to best preserve the modal information (Flanigan, 1998), and the ultimate goal is to expand the mode shapes. Let the complete i th experimental mode shape φ φ φ i =[φ φ φ uT i φ φ φ kT i ] T , estimated from measured re- sponses, be partitioned into unknown φ φ φ u i at the DOFs that are not instrumented and the known φ φ φ k i (n k -element column vector) at the DOFs with sensors. Then, the first step of the dynamic ex- pansion method is to rearrange and partition the structural matrices and dynamic stiffness matrix S i =λ 2 i M+λ i C+ K=λ 2 i ( ˆ M+∑ h µ h ˆ M h )+λ i ( ˆ C+∑ h χ h ˆ C h )+( ˆ K+∑ h κ h ˆ K h ) as M=    M uu M uk M ku M kk    , C=    C uu C uk C ku C kk    , K=    K uu K uk K ku K kk    , S i =    S iuu S iuk S iku S ikk    , (3.12) with dimensions compatible with the partition of the mode shape. Substituting these partitioned matrices into (3.2) and solving forφ φ φ u i yields φ φ φ i =      φ φ φ u i φ φ φ k i      =    − S i − 1 uu S iuk I    φ φ φ k i (3.13) where I is an identity matrix of appropriate dimensions. In practice, M, C and K are unknown, and can only be approximated with (3.3); thus, they, and theφ φ φ i , are functions ofµ,χ andκ and so areφ φ φ i in (3.2). While some previous studies (Wang et al., 2015; Li et al., 2020) have used ˆ M, ˆ C and ˆ K instead of M, C and K in the model expansion procedure because of computational efficiency, this approximation is inaccurate compared with the 55 following approach. Let FEM structural matrices and their low-rank changes be rearranged and partitioned as ˆ M=    ˆ M uu ˆ M uk ˆ M ku ˆ M kk    , ˆ C=    ˆ C uu ˆ C uk ˆ C ku ˆ C kk    , ˆ K=    ˆ K uu ˆ K uk ˆ K ku ˆ K kk    , ˆ M h =    ˆ M huu ˆ M huk ˆ M hku ˆ M hkk    , ˆ C h =    ˆ C huu ˆ C huk ˆ C hku ˆ C hkk    , ˆ K h =    ˆ K huu ˆ K huk ˆ K hku ˆ K hkk    , (3.14) with dimensions compatible with the partition of the mode shape, then, in (3.13), M uu = ˆ M uu + ∑ h µ h ˆ M huu , C uu = ˆ C uu +∑ h χ h ˆ C huu , K uu = ˆ K uu +∑ h κ h ˆ K huu , M uk = ˆ M uk +∑ h µ h ˆ M huk , C uk = ˆ C uk +∑ h χ h ˆ C huk and K uk = ˆ K uk +∑ h κ h ˆ K huk . Thus, (3.5) can be rewritten as ˆ φ φ φ H j    S iuu S iuk S iku S ikk       − S i − 1 uu S iuk I    φ φ φ k i = ˆ φ φ φ H j    0 S ikk − S iku S i − 1 uu S iuk    φ φ φ k i = ˆ φ φ φ k j H (S ikk − S iku S i − 1 uu S iuk )φ φ φ k i = 0 (3.15) where ˆ φ φ φ k j is obtained by partitioning ˆ φ φ φ j =[ ˆ φ φ φ uT j ˆ φ φ φ kT j ] T in the same way. Linearizing S ikk − S iku S i − 1 uu S iuk with respect to (w.r.t.) µ,χ andκ leads to S ikk − S iku S i − 1 uu S iuk ≈ ˆ S ikk − ˆ S iku ˆ S i − 1 uu ˆ S iuk + h ˆ S iku ˆ S i − 1 uu − I i ∑ h µ h ˆ M h + ∑ h χ h ˆ C h + ∑ h κ h ˆ K h !    ˆ S i − 1 uu ˆ S iuk − I    (3.16) 56 where ˆ S i =λ 2 i ˆ M+λ i ˆ M+ ˆ K is partitioned like matrix S i into ˆ S i =[[ ˆ S i T uu ˆ S i T ku ] T [ ˆ S i T uk ˆ S i T kk ] T ] T (see Appendix C for derivations). Thus, for the case that requires matrix expansion, the A matrix in (3.10) is redefined as A= h A M 1 A M 2 ··· A C 1 A C 2 ··· A K 1 A K 2 ··· i = h F M 1 w F M 2 w ··· F C 1 w F C 2 w ··· F K 1 w F K 2 w ··· i (3.17) where w= h φ φ φ k 1 T φ φ φ k 2 T ··· φ φ φ k T n exp λ 1 λ 2 ··· λ n exp i T . The((i− 1)n FEM + j) th row of F M h is h F M h i (i− 1)n FEM +j = ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    P i , i= 1,...,n exp , j= 1,...,n FEM (3.18) where P i = h 0 n k × (i− 1)n k I n k × n k 0 n k × (n exp − i)n k 0 n k × n exp i andh ranges from 1 to the number of mass perturbations ˆ M h . The formulas for F C h and F K h are the same as that for F M h except replacing every ˆ M h with ˆ C h and ˆ K h , respectively, and replacing λ 2 i with λ i and 1, respectively. The b vector in (3.10) is redefined as b= F b w where the((i− 1)n FEM + j) th row of F b is h F b i (i− 1)n FEM +j =− ˆ φ φ φ k j H ( ˆ S ikk − ˆ S iku ˆ S i − 1 uu ˆ S iuk )P i . (3.19) 3.2 Comparison of TLS, CTLS and WTLS Approaches 3.2.1 TLS Approach Before the development of the TLS, the least squares (LS) approach was widely used to solve the overdetermined EIV problem Aˆ x≈ b under the assumption that the error lies entirely (or mostly) in b. In this case, ˆ x LS =(A H A) − 1 A H b uses the Moore-Penrose pseudo-inverse; this solution may be far from accurate when A is not exact and the condition number of A H A is large because a small 57 change (uncertainty) in A will produce a large change in solution ˆ x LS . The TLS was introduced by Golub and Van Loan (1980) and Golub (1973) to address this problem. TLS assumes that the following equality (A+∆A)ˆ x TLS = b+∆b (3.20) and the solution ˆ x TLS = min ˆ x ∥[∆A ∆b]∥ F =− V nq v q /∥v q ∥ 2 2 minimizes the Frobenius norm of the matrix of uncertainties, where then× q matrix V nq andq× 1 vector v q — along withn× p matrix V np and p× 1 vector v p where p=n− q+ 1 — are partitions of the right singular matrix of the singular value decomposition of[A b]: [A b]= U p u q    Σ Σ Σ p 0 0 σ σ σ q            V np V nq v T p v T q         T (3.21) whereΣ Σ Σ p andσ σ σ q are diagonal matrices containing the p largest andq smallest singular values, re- spectively; and the unity matrices are partitioned accordingly. It should be noted that an alternative TLS,i.e., mixed TLS, is also available to deal with cases in which certain columns of[∆A ∆b] are known to be zeros (Van Huffel and Vandewalle, 1989); however, this is not used herein because no columns of [∆A ∆b] herein are known to be zeros. For a thorough review of the TLS approach, the interested reader is directed to the abundance of literature on this topic (Golub and Van Loan, 1980; Golub, 1973; Li et al., 2020). Although the TLS approach is a major advancement over the LS approach, it fails to address a key issue: the underlying interdependence of∆A and∆b. Also, in practical applications, determining p against q is also a challenge, especially when there is no obviously large gap between any two adjacent singular values after the singular values are arranged in the descending order — an overly large (small) p may lead to underestimation (overestimation) of the uncertainties and inaccurate ˆ x TLS . A new generation of TLS approaches that account for the correlation between∆A and∆b are generally classified into the structured total least squares and the weighted total least squares (WTLS); the CTLS approach used herein belongs to the former 58 class. Both the CTLS and the WTLS approaches, summarized in the next two subsections, seek to solve the same problem, differing only in how to account for the relationship between∆A and∆b. 3.2.2 WTLS Approach This approach considers the relationship between∆A and∆b by computing the correlation and the cross-correlation between each and every entry of [∆A ∆b], making it obviously better than the TLS, which does not consider the relationship at all. However, it is not as accurate as the CTLS approach when it is feasible to directly mathematically express the relationships between∆A and ∆b. After all, the WTLS approach only uses up to the second-order statistics of the elements in ∆A and∆b. For a more comprehensive introduction to the WTLS approach, the interested reader is directed to the relevant references (Schaffrin and Wieser, 2008; Fang, 2013, 2015). 3.2.3 CTLS Approach The CTLS approach first traces the sources of the errors so as to express ∆A and∆b as functions of these source errors. When combining the CTLS with the CMCM method, the sources of the errors are assumed to be in the experimental mode shapes (after model expansion) and eigenvalues. Thus, the errors include both (functions of) the measurement noise and the inaccuracies introduced during the mode expansion procedure. Assuming model updating leads to the correct structural matrices, then, the inaccuracies introduced by the model expansion will eventually be very small; how- ever, the uncertainties introduced by measurement noise are persistent throughout model updating. Thus, the following derivations assume that the errors are from noisy experimental mode shapes and eigenvalues. Assume thei th exact specimen’s mode shape (at measured DOFs) isφ φ φ k,exact i ; then, the error in φ φ φ k i is∆φ φ φ k i =φ φ φ k i − φ φ φ k,exact i , and the error in λ i is∆λ i =λ exact i − λ i where λ exact i is the exact specimen’s eigenvalue. By combining all source errors into one column vector, the source error vector∆w is defined as ∆w= h ∆φ φ φ k 1 T ∆φ φ φ k 2 T ··· ∆φ φ φ k T n exp ∆λ 1 ∆λ 2 ··· ∆λ n exp i T . Based on these assumptions and notations, the expressions for∆A and∆b can be derived as follows. 59 From (3.17), the errors in matrix A is ∆A= ∆A M 1 ∆A M 2 ··· ∆A C 1 ∆A C 2 ··· ∆A K 1 ∆A K 2 ··· , or a subset of these if the assump- tions of (3.8)–(3.9) are used, where ∆A M h = ¯ F M h ∆w, ∆A C h = ¯ F C h ∆w and ∆A K h = ¯ F K h ∆w and the ((i− 1)n FEM + j) th row of ¯ F M h is obtained by linearizing ∆A M h w.r.t. ∆w (see Appendix D for derivations): h ¯ F M h i (i− 1)n FEM +j = h F M h i (i− 1)n FEM +j + ˆ φ φ φ k j H      h ˆ S ′ iku ˆ S i − 1 uu − ˆ S iku ˆ S i − 1 uu ˆ S i ′ uu ˆ S i − 1 uu 0 i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    + h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S ′ iuk − ˆ S i − 1 uu ˆ S i ′ uu ˆ S i − 1 uu ˆ S iuk 0    + h ˆ S iku ˆ S i − 1 uu − I i (2λ i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I         φ φ φ k i q T i (3.22) where ˆ S ′ i = ∂ ˆ S i /∂λ i = 2λ i ˆ M+ ˆ C is rearranged and partitioned, like ˆ M and ˆ C, into ˆ S ′ i = [ ˆ S ′ iuu ˆ S ′ iuk ; ˆ S ′ iku ˆ S ′ ikk ]; and q T i = h 0 1× (n exp n k +i− 1) 1 0 1× (n exp − i) i . The formulas for ¯ F C h (i− 1)n FEM +j and ¯ F K h (i− 1)n FEM +j are the same as that for ¯ F M h (i− 1)n FEM +j except replacing M with C and K, respectively, replacingλ 2 i ˆ M h withλ i ˆ C h and ˆ K h , respectively, and replacing 2λ i ˆ M h withλ i ˆ C h and 0, respectively. Similarly, the error in vector b is∆b= ¯ F b ∆w, after linearizing∆b w.r.t.∆w, where the((i− 1)n FEM + j) th row of ¯ F b is h ¯ F b i (i− 1)n FEM +j = h F b i (i− 1)n FEM +j + ˆ φ φ φ k j H ˆ S i ′ kk − ˆ S i ′ ku ˆ S i − 1 uu ˆ S iuk − ˆ S iku ˆ S i − 1 uu ˆ S i ′ uk + ˆ S iku ˆ S i − 1 uu ˆ S i ′ uu ˆ S i − 1 uu ˆ S iuk φ φ φ k i q T i . (3.23) Define n exp n FEM × n k n exp matrix G(x) as follows: G(x)≡ G([µ T χ T κ T ] T )= ∑ h µ h ¯ F M h ! + ∑ h χ h ¯ F C h ! + ∑ h κ h ¯ F K h ! − ¯ F b . (3.24) 60 Then the CTLS solution to (3.10) is (Abatzoglou et al., 1991) ˆ x= x CTLS = min x∈R n x [x T − 1][A b] H G(x)G H (x) − 1 [A b][x T − 1] T (3.25) wheren x is the length of the vector[µ T χ T κ T ] T (bounds may be imposed on elements of x if they are known to change within certain ranges). Note that a prerequisite for this solution to hold is that the number of rows in∆w is smaller than the product of number of rows in A andn x + 1; thus, this condition will be checked throughout. This nonlinear optimization problem can be solved either using a Newton’s approach (Abatzoglou et al., 1991) or the iterative method in Bresler and Macov- ski (1986). Because the Newton’s approach was found hard to convergence, the iterative method is chosen with some modifications. The steps for finding the solution are summarized as follows. (1) Initialize by settingr= 0 and initial guess of x r to 0 (i.e., x r = 0). (2) Compute G(x r ) using (3.24). (3) Compute the complex-valued matrix Q(x r ) as Q(x r )=[A b] H G(x r )G H (x r ) − 1 [A b]. (4) BecauseℑQ(x r ) is a real skew-symmetric matrix and x∈R n x ,[x T − 1]ℑQ(x r )[x T − 1] T ≡ 0, so the minimization problem in (3.25) becomes x r+1 = min x∈R n x[x T − 1]ℜQ(x r )[x T − 1] T (lower and upper bounds on x may be enforced here), which can be solved using thequadprog.m (quadratic programming) function running the trust-region-reflective algorithm in M ATLAB ® (MATLAB, 2020; Coleman and Li, 1996). (5) Convergence check:∥x r+1 − x r ∥>tolerance ? If yes, setr=r+1 and go to step (2); if no, let ˆ x= x CTLS = x r+1 and exit. (Note: without using only the real part of Q(x r ) in step (4), it is possible to solve the optimization, but it was found to be much slower to converge, possibly due to some numerical difficulties, than using the real part only; thus, step (4) is crucial for computational efficiency.) 61 Table 3.1: Slope and intercept statistics for fitting a line to 100 points Mean [Relative Standard Deviation/Relative Root-Mean-Square Error] Exact Value LS TLS CTLS WTLS ˆ α 20 20.7591 [2.85%/4.75%] 20.2870 [2.78%/3.13%] 20.0137 [0.24%/0.25%] 20.0012 [3.86%/3.86%] ˆ β 2 1.9854 [0.84%/1.11%] 1.9948 [0.84%/0.88%] 1.9960 [0.42%/0.46%] 2.0011 [3.24%/3.24%] 3.2.4 Line-Fitting Example A simple line-fitting problem is used to establish a comparison between the TLS, the CTLS and the WTLS approaches. In short, the goal of this problem is to fit a straight line described by v= αu+β given noisy data, in which u and v are coordinates of the to-be-fitted line. In this example, the data (before corrupting it with noise) falls along the exact line v=αu+β where α = 2 andβ = 20. One hundred random points(u i ,v i ) are generated as follows: (1) Lett i ∼ U(0,100),i= 1,...,100, be a set of values uniformly distributed between 0 and 100. (2) Let u i =t i + 0.05t i ξ i and v i =αt i +β+ 0.01(αt i +β)ξ i where ξ i ∼ U(− 1,1) is a noise term. (Note that the noise in thex and y coordinates are intentionally chosen as fully correlated because this is the easiest to demonstrate the differences between the TLS and CTLS solutions.) Given the noisy data represented by u=[u 1 u 2 ... u 100 ] T and v=[v 1 v 2 ... v 100 ] T , the line fitting problem is solved using the LS, the TLS, the CTLS and the WTLS approaches. In this example, the A matrix is A=[u 1], where 1 is a column vector of ones, and vector b= v. The points from one realization of the random points and the exact line are shown in Figure 3.1. For this problem,∆A and∆b are defined as ∆A=[∆u 0] and∆b=∆v. In this example, the number of elements in the noise vectorξ ξ ξ is 100, which is smaller than the product of number of rows in A (which is 100) and n x + 1 (which is 3), thus, the prerequisite for (3.25) to hold is satisfied. To develop statistics (mean and standard deviation) of the estimates ˆ α and ˆ β by the four methods, 1000 independent trials are performed (each with 100 points); these statistics are provided in Table 3.1. According to Table 3.1, the LS approach leads to the least accurate results because, as men- tioned previously, it incorrectly assumes∆A= 0, which the other three methods do not; the TLS 62 0 50 100 0 20 100 220 u v Exact line Noisy coordinates Figure 3.1: Exact line and noisy coordinates approach performs better than the LS approach in terms of the accuracy of mean values, but worse than the CTLS and the WTLS approaches in general, because it does not account for the relation- ship between∆A and∆b; the WTLS approach performs similarly to the CTLS approach in terms of the accuracy of mean values, but much worse than the CTLS approach in terms of standard de- viations, because it only takes advantage of up to the second-order statistics of∆A and∆b instead of using their relationships directly. Thus, the CTLS approach is preferred in this study. 3.2.5 Algorithm for Model Updating For a generic model updating problem, in which the parameters (in a candidate set of to-be-updated parameters) that need updating is unknown (as opposed to the two known parameters,α andβ, in the line-fitting example), determining which parameter(s) need(s) to be changed requires additional mathematical/programming strategies. In a previous study (Li et al., 2020), an iterative strategy for damage detection was proposed. First, it assumes that each ˆ K h is the stiffness matrix of an individual element. After solving (3.10) for the (stiffness) perturbation magnitudes ˆ x, if any ˆ x h is smaller than –1 or larger than 0, then the element is assumed to be undamaged, and ˆ x h is set to zero in the next iteration and the cor- responding column of A is removed, resulting in the A matrix getting thinner and, hopefully, the 63 accuracy of the solution will improve. While this iterative strategy may be appropriate for a dam- age detection problem, it is poorly suited for model updating because the range of ˆ x h for updating is not as clear (can no longer assume− 1≤ ˆ x h ≤ 0). Further, a positive ˆ x h in one iteration, if not removed from A in (3.10), may turn out to be negative in the next iteration. Thus, the rationale for deleting elementh when going from one iteration to the next may be faulty, since ˆ x h changing sign could simply be because the model updating has not converged (the previous strategy only has one iteration loop, which roughly corresponds to the inner-loop iteration of the strategy to be proposed next). In this study, a different iterative strategy used with the CTLS approach for CMCM model updating is proposed as follows: (1) Assume, via engineering judgment, an appropriate range for how much the stiffness of each member can change. Denote the vectors of lower and upper bounds of ˆ x as LB and UB; i.e., the constraint LB≤ ˆ x≤ UB will be enforced throughout the model updating process (for the examples in Section 3.3.2, entries of LB and UB are –0.25 and 0.25, respectively). (2) Initialization: setr= 0 and ˆ x= 0; determinen exp andn FEM (for the examples in Section 3.3.2, n exp = 6 andn FEM = 40). (3) Letr=r+1,s= 0, LB r = LB− ˆ x and UB r = UB− ˆ x. Formulate (3.10) assuming the stiffness of every member is subject to change (n x equals to the total number of members, which is 40 for the structure considered herein). (4) Let s= s+ 1. Estimate the value of stiffness change magnitudes vector ˆ x, denoted ˆ x r,s for the s th inner-loop iteration of the r th outer-loop iteration, by solving (3.10) under the condition LB r ≤ ˆ x r,s ≤ UB r . If any entry of|ˆ x r,s | * has magnitude smaller than a preset value (e.g., 0.01), called the “resolution”, of the model updating algorithm, the value of that entry is assumed to be an artifact of noise and is set to zero. Then, if ˆ x r,s has no zero-entry, or the number of non-zero entries is too small to satisfy the prerequisite for (3.25) to hold (fewer non-zero entries corresponds to fewer columns in the matrix A in the coming inner-loop iteration, which makes satisfying the * |·| denotes the element-wise absolute value operator. 64 prerequisite more difficult), let ˆ x r = ˆ x r,s and jump out of the inner-loop by going to step (6); other- wise, go to step (5). (5) Remove from the current matrix A the columns corresponding to the zero-entries of ˆ x r,s , refor- mulate LB r and UB r by removing the entries corresponding to the zero-entries of ˆ x r,s , and update n x to reflect the number of columns of the current matrix A. Then, go to step (4). (6) Let ˆ x= ˆ x+ ˆ x r . If max h | ˆ x r h |>“resolution”, go to step (3); otherwise, ˆ x is the final solution. This strategy is different from the one in Li et al. (2020) in four main aspects: (a) the outer- loop updates the unupdated FEM based on the current ˆ x, reducing the errors introduced by model expansion as the unupdated FEM gets closer to the true FEM; (b) the outer-loop allows the stiffness of each member to change as opposed to only a subset of members; (c) the loop suppresses the influence of noise by rejecting small stiffness changes ( i.e., small changes to ˆ x) that are below the resolution of the model updating algorithm given the specified level of noise (applicable to both model updating and damage detection); (d) when solving for ˆ x r , the bounds are enforced directly, which generally leads to more accurate results and/or fewer iterations than the previous study (Li et al., 2020) approach that ignores the bounds while solving for the perturbation magnitude and then saturate to within the bounds afterwards (this is tailored for use with the CTLS approach, not with the LS and TLS approaches). 3.3 Numerical Example: 3D Frame Structure 3.3.1 Introduction to FEM To compare the TLS and the CTLS approaches in the context of solving (3.10), a numerical ex- ample of the dynamics of a 3D frame structure, which was utilized in previous CMCM studies (Hu et al., 2006; Li et al., 2020), is used here. The base of the structure is rigidly connected to the ground; the dimensions of the model and the node numbers are specified in Figure 3.2 (not in- cluding the base, 40 members in total). Each member is modeled as a beam element with Young’s modulus 210 GPa, cross-section area 28.25 cm 2 , and moment of inertia 289 cm 4 (these parameters 65 are the same as in previous studies of this structure (Hu et al., 2006; Li et al., 2020), which specify neither density nor Poisson’s ratio, but using ρ = 7800 kg/m 3 and ν = 0 gives frequencies that match those references, so these values are used in this example). This FEM is considered the unupdated FEM. The data will be generated from a “damaged” FEM with the Young’s moduli of elements 6–8 (dashed in Figure 3.2) and 9–10 (dotted) both reduced by 20% (it is assumed for this example that only stiffness is unknown: mass is assumed known and damping is ignored). Where parameter constraints are needed in applying CMCM to this structure, the stiffness of each mem- ber is allowed to change by± 25% (a range that includes the magnitudes of the assumed stiffness reductions). Table 3.2 shows a comparison of the unupdated FEM and the “true” FEM natural frequencies and mode shapes (measured by the modal assurance criteria (Allemang and Brown, 1982a), or MAC, defined as MAC (φ φ φ i ,φ φ φ j )=|φ φ φ H i φ φ φ j | 2 /(φ φ φ H i φ φ φ i φ φ φ H j φ φ φ j )) of the first six modes. The goal is to update the stiffness matrix of the unupdated FEM based on the (noise-free or noise- polluted) modal information of the “true” FEM, so the model matches, as best as possible, the “true” FEM. 0 1 2 y [m] 3 z [m] 0 4 5 x [m] 0 1.5 0.5 Figure 3.2: FEM of 3D frame structure 66 Table 3.2: Comparing modal characteristics of the unupdated and “damaged” FEMs Mode # Mode Type Natural Frequencies [Hz] MAC True FEM Unupdated FEM (∆ ∗ ) 1 1st bending iny direction 6.797 6.884 (1.28%) 0.9999 2 1st bending inx direction 9.300 9.364 (0.69%) 0.9877 3 1st torsion aroundz axis 9.831 9.885 (0.56%) 0.9929 4 2nd bending iny direction 22.310 22.317 (0.03%) 1.0000 5 2nd bending inx direction 29.389 29.560 (0.58%) 0.9837 6 2nd torsion aroundz axis 30.611 30.730 (0.39%) 0.9887 * change of the unupdated FEM natural frequency relative to the true FEM natural frequency 3.3.2 Examples and Results The examples analyzed in this section have two main objectives: first, demonstrate that the CTLS approach is superior to the TLS approach (Section 3.3.2.1); second, evaluate the performance of the CTLS approach under different levels of noise and different selections of “resolution” (Sec- tion 3.3.2.2). 3.3.2.1 Comparison of the TLS and CTLS Approaches When the noise-free modal information of the true FEM is known at all DOFs, both the TLS and the CTLS approaches lead to the exact solution; the results of this trivial case are omitted. The examples used herein apply both methods to perform model updating (not damage detection) together with the proposed strategy from Section 3.2.5. To resemble the “spatially incomplete and noise-polluted” example in Li et al. (2020), the true FEM’s natural frequencies are assumed noise- free but the mode shapes’ elements have independent zero-mean Gaussian noise with standard deviations 0.1% of the elements” magnitudes. The number of elements in the noise vector∆w is n exp × n k = 6× 80= 480, which is smaller than the product of number of rows in A (n exp × n FEM = 67 6× 40= 240) and n x + 1= 2+ 1= 3; thus, the prerequisite for (3.25) to hold is satisfied as long as n x remains no smaller than 2 in step (4) in Section 3.2.5. The “resolution” is first set to be 1%; and the effect of this value on the resulting ˆ x will be discussed afterwards in Section 3.3.2.2. The model updating problem is solved 100 times with independent realizations of the mode shape noise by both the TLS and the CTLS approaches; the statistics of the stiffness changes for both the damaged and the undamaged elements are shown in Table 3.3. Note that the statis- tics of the undamaged elements are obtained by finding the largest absolute stiffness change (a positive value) among the undamaged elements in each (TLS or CTLS) run and computing the statistics (mean and standard deviation) of the 100 values from the 100 runs. Table 3.3: Average and standard deviation of magnitude change ˆ x for comparison of the TLS and the CTLS approaches Average [Relative Standard Deviation] Exact TLS CTLS Element 6–8 –0.2000 –0.2078 [6.224%] –0.2004 [0.037%] Element 9–10 –0.2000 –0.1996 [3.294%] –0.2023 [0.042%] Undamaged Elements 0.0000 0.0157 [0.0116 † ] 0.0000 [0.0000 † ] † An undamaged element’s exact stiffness change is 0, so this standard deviation is not normalized. Table 3.3 demonstrates that, under the given noise level, the TLS and the CTLS have similar accuracy in term of mean values when estimating the stiffness changes of the damaged elements; however, the CTLS leads to much smaller corresponding error standard deviations. Also, the CTLS exactly estimates the (zero) stiffness changes of the undamaged elements, yet the TLS leads to some (mildly) incorrect estimations. Thus, the CTLS is clearly superior to the TLS. 3.3.2.2 Performance Evaluation of the CTLS Approach Section 3.3.2.1 demonstrates that the CTLS performs better than the TLS; however, two questions arise in applying the CTLS approach: (a) how does the “resolution” affect the estimation of ˆ x? (b) how much noise can the CTLS approach handle? Note that noise considered in the previous 68 section (also in Li et al. (2020)) has two main issues: first, the possible noise in frequency mea- surements is not considered; second, the noise level for the mode shape measurements, which is 0.1% of the absolute value of the exact value, is rather lower than is typical in real applications. Thus, when addressing the second question, the frequency measurements of the true FEM will be polluted by noise and the mode shape measurements will be polluted by a higher level of noise than before. To address question (a), two cases are considered, with the “resolution” set as 1.0% and 1.5% for cases A and B, respectively. In both cases, the error added to each frequency measurement is zero-mean with a standard deviation 0.015% of the exact value, and the error added to each element of the mode shape vector is zero-mean with a standard deviation 1.5% of (the absolute value of) the exact value (errors are all independently distributed). Note that the noise level for frequencies is smaller because the frequency differences between the true and unupdated FEMs are very small according to Table 3.2. These noise levels are much larger and more realistic than those in Li et al. (2020), and solving this problem using the TLS approach leads to obviously erroneous ˆ x — the TLS approach fails to handle such levels of noise that are realistic in the practical applications; again, indicating the CTLS approach is better than the TLS approach. Further, since the noise in identified frequencies are considered, the dimension of ∆w is increased by n exp = 6 to 486; it can be easily verified, like in the previous subsection, that the prerequisite for (3.25) to hold is still satisfied. Model updating for both cases is conducted 100 times with independent realizations of the frequency and mode shape noise, and the statistics of the results are summarized in Table 3.4. Table 3.4: Average and standard deviation of magnitude change ˆ x estimated by the CTLS approach for different noise levels Average [Relative Standard Deviation] Exact Case C Res.=1.0% Case D Res.=1.0% Case A Res.=1.0% Case B Res.=1.5% Element 6–8 –0.2000 –0.2000 [0.098%] –0.2000 [0.222%] –0.1995 [2.177%] –0.1999 [0.542%] Element 9–10 –0.2000 –0.2000 [0.080%] –0.2000 [0.715%] –0.2012 [5.735%] –0.2021 [0.864%] Undamaged Elements 0.0000 0.0000 [0.0000 † ] 0.0013 [0.0093 † ] 0.0058 [0.0223 † ] 0.0018 [0.0093 † ] † An undamaged element’s exact stiffness change is 0, so this standard deviation is not normalized. Res. is short for “Resolution”. 69 Table 3.3 demonstrates that, given the assumed realistic noise level, cases A and B both lead to reasonably good results. Case B (resolution=1.5%) leads to better ˆ x estimations compared with case A (resolution=1.0%) in terms of standard deviations for all ˆ x entries as wells as mean estimations for ˆ x entries corresponding to the undamaged elements, because the 1.0% resolution in case A is not sufficient to reject the noise satisfactorily. However, it should be noted that selecting too large a resolution may lead to small changes in stiffness being ignored; thus, in reality, the resolution should initially be chosen as a very small value and increased gradually until the results do not improve. In this example, for each structural member, a resolution of 1.0% of the unupdated FEM stiffness can be arguably considered small enough, and because such a resolution leads to reasonably good ˆ x estimations under realistic noise levels, the effectiveness of the proposed CTLS approach has been proven. To address question (b), two cases (cases C and D) are considered in addition to case A, with noise levels in both frequencies and mode shapes one-third and two-thirds of those in case A for cases C and D, respectively (other parameters are the same as those in case A). Model updating for cases C and D is conducted 100 times with independent realizations of the frequency and mode shape noise, and the statistics of the results are summarized in Table 3.4. According to Table 3.4, the algorithm performs better under lower noise levels, and vice versa; the magnitude estimation is almost exact for case C, not as good in case D in terms of the estimates for the undamaged elements as wells as all the standard deviations, and worst in Case A though the estimates may still be deemed to be good in practice. 3.4 Experimental Study: Base-isolated Reinforced-Concrete Frame Building This section focuses on validating the proposed model updating algorithm by implementing it on the realistic base-isolated building and its nominal FEM introduced in Sections 2.1–2.4. 70 3.4.1 Model Updating Procedure Although most of the theoretical aspects have been covered up until this point, there are still certain details crucial to effectively and accurately updating the FEM of this particular structure. (1) Given that the damping in the superstructure is small, it is fair to assume (as in many studies of lightly damped systems) that the damping in the superstructure has infinitesimal influence on the frequencies and the mode shapes. Also, given that the linear damping provided by the isolation- layer devices is mild, it is assumed that every isolation-layer device can be modeled with a constant stiffness plus a linear viscous damping with constant damping coefficient. (2) The second detail concerns which modes to select when updating certain parameters. When updating the stiffnesses of members in the superstructure, the first three modes are ignored be- cause they are very insensitive to superstructure parameters and preliminary studies found that including the isolation modes injected greater uncertainty than helpful information. When updat- ing the stiffnesses and the damping coefficients of the isolation-layer devices, ideally, only the first three modes should be used; however, because there are 24 parameters that need updating (twelve stiffnesses and twelve damping coefficients) and using the first three modes resulting in matrix A having only 3× 3= 9 < 24 rows, more modes need to be incorporated. Thus, the eight experi- mental modes, the three FEM isolation modes and the five FEM superstructure modes that have already been matched to their experimental counterparts, are used to update the isolation-layer parameters. It should be noted that other FEM modes can also be used in theory, however, using these modes may unnecessarily introduce uncertainties into the system; thus, using these modes is not suggested when better options are available. (3) It should be noted that, when non-proportional damping is considered, each mode has two state- space complex eigenvalues (namely,λ i andλ ∗ i for thei th mode) and two corresponding state-space complex mode shapes. Since the damping is, strictly speaking, nonproportional in this system, the modal information from the system identification always follows this rule. However, for the FEM, this is the case only after the damping of the isolation-layer devices is considered. In other words, 71 when the damping of the isolation-layer devices is not considered (superstructure damping is neg- ligible, so C= 0), which is the case when updating superstructure parameters, each FEM mode has one real eigenvalue and one real mode shape, and the following two equations are identical (M and K are real and symmetric, so ˆ φ φ φ j is real). Thus, only one of (3.26a) and (3.26b) should be used when formulating the A matrix and the b vector in (3.10). ˆ φ φ φ H j λ 2 i Mφ φ φ i + ˆ φ φ φ H j λ i Cφ φ φ i + ˆ φ φ φ H j Kφ φ φ i = 0 (3.26a) ˆ φ φ φ H j λ ∗ 2 i Mφ φ φ ∗ i + ˆ φ φ φ H j λ ∗ i Cφ φ φ ∗ i + ˆ φ φ φ H j Kφ φ φ ∗ i = 0 (3.26b) where the mass matrix M of the experimental building is assumed the same as that of the FEM, i.e., M= ˆ M; the stiffness matrix K is constructed following (3.3), where ˆ K h is either the stiffness matrix corresponding to each of the 24 different combinations of superstructure components — including the beams under floors 0–4 (5 matrices), the columns in stories 1–4 (4 matrices), the shear walls in each story (4 matrices), the floor slabs on each floor (5 matrices), the stairs in stories 1–3 (3 matrices) and the nonstructural walls in stories 2–4 (3 matrices), when updating the superstructure parameters, or the stiffness matrix of an isolation-layer device in thex ory direction when updating the isolation-layer parameters (twelve parameters and twelve corresponding ˆ K h matrices); the damping matrix C is ignored when updating the superstructure parameters (only superstructure stiffnesses are updated); when updating the isolation-layer parameters (stiffnessnes and damping coefficients in x and y direction for all isolation-layer devices), ˆ C is the damping matrix constructed based on the twelve nominal damping coefficients ( ˆc) listed in Table 3.5 and ˆ C h is the damping matrix corresponding to each isolation-layer device damping coefficient in the x or y direction (twelve ˆ C h matrices in total). (4) When updating the isolation-layer parameters, ˆ φ φ φ j is no longer real and (3.26a)–(3.26b) are not equivalent, so it is desirable to use both equations to make use of all the available information. However, in such a case, the errors in ˆ φ φ φ j and ˆ φ φ φ ∗ j must be conjugate of each other — a condition 72 very challenging to enforce directly. Thus, to avoid using ˆ φ φ φ ∗ j , (3.27) is used instead of (3.26b), which is obtained by taking the conjugate of left- and right-hand-side sides of (3.26b). ˆ φ φ φ T j λ 2 i Mφ φ φ i + ˆ φ φ φ T j λ i Cφ φ φ i + ˆ φ φ φ T j Kφ φ φ i = 0 (3.27) (5) In all model updating strategies, a certain amount of expert judgment is required to decide how much the to-be-updated parameters can deviate from their nominal values. The upper and lower bounds of each isolation-layer device parameter (either linear stiffness or viscous damping coeffi- cient), summarized in Table 3.5, were determined in a previous study (Yu et al., 2022). Although the Young’s modulus of nonstructural walls (autoclaved aerated concrete panels, also ALC panels) is around 2.0 GPa, because the installation allows some in-plane sliding between the adjacent pan- els and out-of-plane rotation where connected to the main structure, the lower and upper bounds of the effective Young’s modulus are set to be zero and 2.5 GPa, respectively. The Young’s moduli of the concrete are allowed to change 5% from their nominal values, which is arguably a reasonable assumption in the real engineering practice. Based on the aforementioned details, the following procedures are applied. (a) The unupdated FEM modes above the highest-frequency experimental mode are ignored; thus, a total ofn FEM = 9 modes of of the unupdated FEM are used. (b) FEM modes 4–9 and experimental modes 4–8 (i.e., the superstructure modes) are used to update the stiffnesses of superstructure components (24 parameters in total, as mentioned in the details), repeating iteratively until their changes are negligible. Since no damping has been added to the FEM up until this point, ˆ φ φ φ j remains real. The A matrix in (3.10) is 5(n FEM − 3)× 24 or 30× 24. The number of elements in the noise vector∆w is n exp × n k +n exp = 5× 42+ 5= 215, which is smaller than 750 (the product of 30, the number of rows in A, andn x + 1= 24+ 1= 25); thus, the prerequisite for (3.25) to hold is satisfied. (c) Experimental modes 1–3, the three FEM isolation modes and the five FEM superstructure modes that have been matched in (b) are used to update the isolation-layer device parameters. In 73 this case, there are three φ φ φ i (i= 1,2,3), eight ˆ φ φ φ j , and eight ˆ φ φ φ ∗ j (j = 1,2,3,4,5,6,7,9). Strictly speaking, ˆ φ φ φ j is (close to) a real vector, i.e., ˆ φ φ φ j = ˆ φ φ φ ∗ j , for j = 4,5,6,7,9, thus, there are three ˆ φ φ φ ∗ j (j= 1,2,3). Then, the 33× 24 matrix A and the 33× 1 vector b are formulated using (3.26a) and (3.27). Similar to (b), perform the updating iteratively until the parameter changes are negligible. The number of elements in the noise vector∆w is n exp × n k +n exp = 3× 42+ 3= 129, which is smaller than 825 (the product of 33, the number of rows in A, andn x + 1= 24+ 1= 25), thus, the prerequisite for (3.25) to hold is satisfied. 3.4.2 Model Updating Results The updated Young’s moduli for the nonstructural walls are around 30 MPa (around two orders of magnitude smaller than the nominal value and three orders of magnitude smaller than those of the concrete in this study), which is expected, given how the ALC panels are installed. For the updated Young’s moduli of the concrete and the updated stiffnesses and damping coefficients of isolation-layer devices, see Table 2.5 and Table 3.5, respectively. Table 3.5: Nominal and updated stiffnesses and damping coefficients of isolation-layer devices RB1 RB2 ESB1 ESB2 SDP1 SDP2 x y x y x y x y x y x y ˆ k [kN/m] 1128 1104 1077 1115 1539 1413 1495 1453 3936 4070 4004 3772 LB 1062 1026 1016 1044 1367 1225 1022 1096 3634 3689 3634 3387 UB 1195 1182 1138 1187 1711 1602 1969 1810 4239 4452 4373 4157 updatedk 1184 1051 1138 1184 1611 1602 1616 1382 3852 3819 3963 4157 ˆ c [Mg/s] 16.74 17.35 15.55 17.84 26.23 29.84 60.06 55.55 101.75 106.65 108.67 119.25 LB 15.88 0.00 1.36 1.17 0.00 0.00 0.00 0.00 30.84 19.15 22.08 30.06 UB 31.90 35.25 29.73 34.51 65.35 73.53 170.20 137.53 172.66 194.15 195.26 208.45 updatedc 26.51 26.48 27.37 25.20 45.87 45.87 92.10 72.19 162.67 129.42 183.16 182.16 As seen in the frequency comparisons in Table 3.6, the errors in the natural frequencies of the superstructure modes dropped significantly after model updating, especially for modes 4–6 dropping from 25%–29% to less than 1%. The mode shapes comparisons in Figure 3.3 show that the mode shapes of the isolation modes match much better after model updating, with MAC values rising from as low as 0.68 (for mode 2 before updating) to larger than 0.95 for all three modes; 74 the mode shape of mode 9 also matched much better. Further, with the updated isolation-layer device damping ratios listed in Table 3.5, the damping ratios of the first three modes (and their differences relative to the identified damping ratios in Table 2.2) are 7.64 (0.13%), 8.47 (–2.18%) and 8.11 (2.13%), indicating very good matches. Damping ratios of superstructure FEM modes that are matched to the identified superstructure modes will be taken from Table 2.2, whereas damping ratios of all other superstructure modes will be taken as the average of the damping ratios of identified modes 4–8 in Table 2.2, for the time-history analyses in the next section. One metric for evaluating whether a subset of modes captures the salient dynamics of the structure is to quantify the fraction of effective mass ratio provided by the mode subset. Let ˆ M x i , ˆ M y i and ˆ M z i denote the effective modal masses of mode i in the x, y and z directions, and let ˆ M (·) =∑ n FEM i=1 ˆ M (·) i denote the effective modal mass sum (theoretically, ˆ M x , ˆ M y and ˆ M z should be the same; in practice, they may not exactly do so, and are replaced by an aggregate total M total , which is assumed to be the sum of all modal effective masses of modes with natural frequencies smaller than 35 Hz in the corresponding direction). Similarly, let ˆ I x i , ˆ I y i and ˆ I z i be the effective modal moments of inertia, and ˆ I x , ˆ I y and ˆ I z be their respective sums. Given these definitions, Table 3.7 shows the effective mass ratios and moment-of-inertia ratios provided by FEM modes 1–7 and 9, demonstrating that this subset of modes, contributing almost all of the effective mass/moment of inertia in the x and y directions and rotation around the vertical z axis, and most of them for the other directions. 3.4.3 Time Histories from Experimental and Updated FEM Time-history analysis was performed on the updated FEM, in which a modal approach was used to simulate the response of each mode with a natural frequency smaller than 35 Hz (first 36 FEM modes) to the 1 kHz table-center excitations (six excitations in thex,y,z, Rx, Ry and Rz directions) 75 Table 3.6: Frequency comparisons of the experimental building, and the unupdated and updated FEMs Experimental (N4SID) Unupdated FEM Updated FEM Mode # Freq. [Hz] Mode # Freq. [Hz] (Err. [%]) Mode # Freq. [Hz] (Err. [%]) 1 0.685 1 0.676 (–1.352) 1 0.688 ( 0.351) 2 0.698 2 0.688 (–1.355) 2 0.696 (–0.228) 3 0.710 3 0.703 (–0.912) 3 0.704 (–0.815) 4 4.781 4 5.626 (17.668) 4 4.759 (–0.469) 5 5.175 5 6.666 (28.813) 5 5.226 ( 0.993) 6 7.293 6 9.149 (25.447) 6 7.281 (–0.160) 7 10.836 7 10.728 (–1.000) 7 10.718 (–1.092) 8 15.504 N/A 8 12.223 N/A 8 15.346 9 15.763 ( 2.540) 9 15.241 (–0.683) 0.676 0.688 0.703 5.626 6.666 9.149 10.728 15.504 15.763 0.685 0.698 0.710 4.781 5.175 7.293 10.836 15.346 0 0.2 0.4 0.6 0.8 1 MAC (a) Before updating 0.688 0.698 0.706 4.767 5.239 7.297 10.724 12.291 15.247 0.685 0.698 0.710 4.781 5.175 7.293 10.836 15.346 0 0.2 0.4 0.6 0.8 1 MAC (b) After updating Figure 3.3: Comparison of FEM and experimental mode shapes as measured by the MAC value as piecewise linear; global responses are then modal coordinates times the corresponding mode shapes. A root mean square (RMS) error metric is introduced to quantify the difference between an experimentaly exp (t) and updated FEMy FEM (t) time histories. Err RMS = |RMS[y exp (t)]− RMS[y FEM (t)]| RMS[y exp (t)] × 100% (3.28) 76 Table 3.7: Effective mass ratio of FEM modes 1–7 and 9 FEM Mode# (i) ˆ M x i /M x ˆ M y i /M y ˆ M z i /M z ˆ I x i /I x ˆ I y i /I y ˆ I z i /I z 1 0.009 0.607 0.000 0.055 0.005 0.004 2 0.980 0.020 0.000 0.002 0.337 0.757 3 0.012 0.373 0.000 0.029 0.002 0.237 4 0.000 0.000 0.000 0.030 0.250 0.001 5 0.000 0.000 0.012 0.130 0.157 0.000 6 0.000 0.000 0.006 0.021 0.028 0.000 7 0.000 0.000 0.882 0.557 0.135 0.000 8 0.000 0.000 0.000 0.000 0.000 0.000 9 0.000 0.000 0.008 0.056 0.015 0.000 Summation 1.000 1.000 0.908 0.880 0.929 0.999 Lighter color indicates FEM modes that are not matched to the experi- mental modes. The summation only includes matched FEM modes. where RMS[y(t)]=[t − 1 f R t f 0 y 2 (t)dt] 1/2 is the root-mean-square response over time interval[0,t f ] of a continuous-time signal, and, for a discrete-time signal,[n − 1 t ∑ n t − 1 k=0 y 2 (k∆t)] 1/2 . Signals compared herein include the acceleration time histories on the superstructure and the time histories of the restoring forces in the isolation layer. Table 3.8 shows the ranges of RMS errors of absolute acceleration time histories at various floor levels in the x, y and z directions, as well as the RMS errors of the isolation-layer device restoring force time histories. After model updating, the maximum percentage errors in x- and y-direction acceleration responses drop by more than half and two-thirds, respectively; the z-direction errors did not change much, but they were already on the same level as the post-updating x- and y-direction errors. The errors in most of the isolation-layer device force predictions drop significantly, except for the little increase and minor decrease of errors in they-direction restoring forces of ESB2 and SDP1, respectively. Figure 3.4 compares the experimental and updated-FEM-predicted accelerations that exhibit the largest errors according to Table 3.8 in the horizontal (Figures 3.4b and 3.4d) and vertical (Fig- ures 3.4f and 3.4h) directions (both cases occur at the accelerometer near grid 1A on floor 0). Visual inspections show that the frequencies and magnitudes of the absolute accelerations are well predicted, except that some high-frequency content in the x-direction is predicted somewhat less 77 accurately. Figure 3.5 compares the experimental and updated-FEM-predicted force-displacement relationships of isolation-layer devices exhibiting the largest errors (according to Table 3.8) among devices of the same kind. The loops for the y-direction RB1 and x-direction SDP1 are well simu- lated, yet the loop fory-direction ESB2 is not simulated as well, mainly because of the mild sliding behavior occurred in the experiment that is not captured in the linear FEM; nevertheless, given the aforementioned good matches of time histories, it is reasonable to conclude that the energy dissi- pation features are reasonably simulated. Table 3.8: Percentage errors in RMS acceleration responses and isolation-layer device forces be- fore (after) updating Direc. Status Accel.* Restoring Forces of Isolation-layer Devices RB1 RB2 ESB1 ESB2 SDP1 SDP2 x before 7.3 – 34.5 13.8 31.0 13.3 26.7 16.9 16.7 x (after) (0.1 – 13.5) (4.8) (3.3) (5.3) (3.1) (11.9) (10.8) y before 9.2 – 26.9 9.9 14.8 14.1 7.8 9.8 14.5 y (after) (1.6 – 9.1) (7.8) (0.4) (5.7) (9.1) (9.2) (3.7) z before 0.0 – 8.5 N/A N/A N/A N/A N/A N/A z (after) (0.8 – 8.3) (N/A) (N/A) (N/A) (N/A) (N/A) (N/A) * the RMS error minima and maxima across the 14 structure accelerometers 78 FEM Experiment (a) (b)x-direction FEM Experiment (c) (d)x-direction, zoomed-in 79 FEM Experiment (e) (f)z-direction FEM Experiment (g) (h)z-direction, zoomed-in Figure 3.4: Experimental and updated-FEM-predicted absolute acceleration time histories at the accelerometer near grid 1A on floor 0: full duration and a corresponding zoomed-in version to show some details (the zoomed-in version for z-direction uses a narrower time window for better visibility of the high-frequency motions) 80 FEM Experiment (a)y-direction RB1 FEM Experiment (b)y-direction ESB2 FEM Experiment (c)x-direction SDP1 Figure 3.5: Experimental and updated-FEM-predicted force-displacement relationships 81 Chapter 4 Real-Time Neural Network Based Semiactive Model Predictive Control of Structural Vibrations 4.1 Problem formulation For annDOF building model subjected to an external excitation and controlled bym smart damping devices, the equation of motion is M ¨ e q(t)+ C ˙ e q(t)+ Ke q(t)= ˆ B u e u(t)+ ˆ B w e w(t) (4.1) where M, C and K are the n× n symmetric mass, damping and stiffness matrices of the building; e q is an n× 1 vector of displacements;e u is an m× 1 vector of control forces with n× m influence matrix ˆ B u that depends on the control device locations; and e w is an n w × 1 vector of external zero-mean excitations, with influence matrix ¯ B w . Equation (4.1) can then be written in state-space form: ˙ e x(t)= ¯ Ae x(t)+ ¯ B u e u(t)+ ¯ B w e w(t) (4.2) with statese x=[e q T ˙ e q T ] T , state matrix ¯ A= [0 n KM − 1 ] T [I n CM − 1 ] T where I n and 0 n are ann× n identity and zero matrices, respectively, and influence matrices ¯ B u =[0 m× n ˆ B T u M − 1 ] T and ¯ B w = [0 n w × n ˆ B T w M − 1 ] T , where 0 (·)× (·) is an all-zero matrix of the given size. With a zero-order-hold time 82 discretization with sampling time interval∆t, the discrete-time state-space equation corresponding to (4.1) is e x l+1 = Ae x l + B u e u l + B w e w l (4.3) wheree x l ,e u l and e w l are the state, control force and excitation vectors, respectively, at time l∆t, A=e ¯ A∆t , B u = ¯ A − 1 (A− I 2n ) ¯ B u , and B w = ¯ A − 1 (A− I 2n ) ¯ B w . The control objective is to minimize a quadratic cost e J=∆t n t ∑ l=0 e x T l Qe x l + 2e x T l Ne u l +e u T l Re u l (4.4) where R is a symmetric positive-definite matrix, [Q N] T [N T R] T is symmetric and positive semidefinite, and n t is the number of time steps (n t ∆t is the final time). At time step l, the goal is to determine the optimal choice of control forcee u l when the state ise x l . In the absence of constraints, the linear quadratic regulator (LQR) minimizes the expected costE[ e J]; for an infinite horizon ( i.e., large n t , tending to infinity), the discrete-time LQR becomes e u l =− K LQR e x l , with control gain K LQR =(B T u Q p B u + R) − 1 (B T u Q p A+ N T ), where Q p is the solution to the algebraic Riccati equation A T Q p A− Q p − (A T Q p B u + N)(B T u Q p B u + R) − 1 (B T u Q p A+ N T )+ Q= 0. Instead, MPC minimizes, subject to the constraints (dissipativity) and state equation, the pre- dicted cost over the next p time steps. Let x k be the model-predicted value of future statee x k+l when the future control force sequence is{u 0 ,...,u p− 1 }, where u k is a possible value of future control forcee u k+l . Herein, the cost is further augmented with a term x T p Q p x p that approximates the residual cost beyond p time steps by emulating a subsequent infinite-horizon discrete LQR cost: J(X,U;x 0 )=∆t x T p Q p x p +∆t p− 1 ∑ k=0    x k u k    T    Q N N T R       x k u k    (4.5) where X=[x T 1 x T 2 ... x T p ] T and U=[u T 0 u T 1 ... u T p− 1 ] T collect the future state predictions and candidate control forces, respectively. The dissipativity constraints for an ideal controllable damper 83 can be expressed as nonlinear equalities sgnv rel k,j = sgnu k,j or nonlinear inequalities v rel k,j u k,j ≥ 0 at time step k for each damper k= 1,2,...,m, where v rel k =− ˆ B T u ˙ q k =[0 m× n − ˆ B T u ]x k = Vx k is the vector of velocities across the damping devices. The constrained minimization of the expected cost E[J] can be recast as an hMPC (MIQP) problem. Define ∆ ∆ ∆≡ [δ δ δ T 1 δ δ δ T 2 ... δ δ δ T p− 1 ] T , where δ δ δ k is an m× 1 vector of binary numbers δ k,j (the “integer” parts of the MIQP) that will, upon solution, beH(v rel k,j ), whereH(·) is the Heaviside unit step function (i.e.,δ k,j = 1 implies v rel k,j > 0);δ δ δ 0 is known a priori because x 0 is known. The MIQP problem in search spaceξ ξ ξ≡ [X T U T ∆ ∆ ∆ T ] T is, then, ξ ξ ξ ∗ = argmin ξ ξ ξ J(ξ ξ ξ ;x 0 ) (4.6a) s.t.: x k+1 = Ax k + B u u k (4.6b) − (v min − ϵ )◦ δ δ δ k ≤ v rel k − v min , k> 0 (4.6c) − (v max +ϵ )◦ δ δ δ k ≤− v rel k − ϵ , k> 0 (4.6d) − (u min − ϵ )◦ δ δ δ k ≤ u k − u min (4.6e) − (u max +ϵ )◦ δ δ δ k ≤− u k − ϵ (4.6f) x min ≤ x k+1 ≤ x max (4.6g) where: k = 0, ..., p− 1 throughout the remainder of the paper unless otherwise noted; ϵ is a tiny positive scalar (required for solvers unable to handle strict inequality constraints); “◦ ” is the Hadamard (element-by-element) product operator of two matrices/vectors; and the lower x min = − x max < 0 and upper x max > 0 bounds on states x k , lower v min =− v max < 0 and upper v max > 0 bounds on cross-damper velocities v rel k , and lower u min =− u max < 0 and upper u max > 0 bounds on control forces u k are introduced as actual bounds — if needed — or as artificial bounds (of sufficiently large magnitude that the corresponding inequalities are never active) that improve the numerical conditioning of the solver. (For the numerical examples herein, the lower and upper bounds for each quantity are of equal magnitude but opposite sign, and chosen large enough that they are inactive.) Constraints (4.6c)–(4.6d) ensure that v min ≤ v rel k ≤ v max (element-wise) and that 84 δ k,j =H(v rel k,j ); constraints (4.6e)–(4.6f) do likewise for u k . Each constraint in (4.6c)–(4.6g) can be cast in the form f T i ξ ξ ξ ≤ h i ; these can be combined to form a vector inequality Fξ ξ ξ ≤ h, where F=[f 1 f 2 ... f 4(m+n)p− 2m ] T and h=[h 1 h 2 ... h 4(m+n)p− 2m ] T . It should be noted that the external excitation w k , if zero-mean, vanishes from (4.6b) as proven in Elhaddad and Johnson (2013). Finally, the components ofξ ξ ξ can be extracted as        X U ∆ ∆ ∆        ≡     I 2np 0 2np× mp 0 2np× m(p− 1) 0 mp× 2np I mp 0 mp× m(p− 1) 0 m(p− 1)× 2np 0 m(p− 1)× mp I m(p− 1)     ξ ξ ξ, (4.7a) where: x k+1 ≡ [0 2n× 2nk I 2n 0 2n× 2n(p− k− 1) ]X, u k ≡ [0 m× mk I m 0 m× m(p− k− 1) ]U and δ δ δ k ≡ [0 m× m(k− 1) I m 0 m× m(p− k− 1) ]∆ ∆ ∆, k> 0; (4.7b) 4.2 Offline neural network training This section focuses on training NNs to predict binary variable vector∆ ∆ ∆ ∗ =[δ δ δ T 1 ... δ δ δ T p− 1 ] T for a given x 0 . Two classes of strategies can be defined. This study focuses on a strategy class, denoted S I herein, that is defined solely by a particular value of binary vector ∆ ∆ ∆. The second strategy class, denotedS II and used in previous studies, is additionally defined by a vector θ θ θ=[θ θ θ T v θ θ θ T u ] T of integer values indicating which constraints (4.6c)–(4.6g) are active; he(mk+ j) th element (1≤ j≤ m) of mp× 1 vectorθ θ θ v is− 1 if the corresopnding element of the lower bound (4.6c) is active,+1 if the upper bound (4.6d) is active, and 0 if both are inactive;i.e., θ v,mk+j =            − 1, (1− δ k,j )v min j +δ k,j ϵ =v rel k,j 0, (1− δ k,j )v min j +δ k,j ϵ <v rel k,j <δ k,j v max j − (1− δ k,j )ϵ +1, v rel k,j =δ k,j v max j − (1− δ k,j )ϵ ; (4.8) θ θ θ u is defined similarly but for constraints (4.6e) and (4.6f). In this study, only fully-connected feedforward NNs are used, in which the numbers of nodes in the input and output layers match the dimensions of the input and the output, respectively, but 85 the numbers of hidden layers and nodes in each hidden layer are chosen case-by-case so that the NNs can be trained to have prediction accuracy (fraction of NN most-likely-strategy predictions that are correct) of around 95% for both the training data (randomly selected from the depositories of all samples ) and the testing data (the remaining data not used for NN training). The activation functions involved are theTanh function for the input and hidden layers and theSoftmax func- tion for the output layer. It should be noted that the NNs are fixed after training — they do not change with regard to either the excitations or the control forces. 4.2.1 Shrinking input parameter space — homogeneity The input to the optimization (4.6) is a vector of system states x 0 ∈R 2n ; the effort for NN training grows exponentially with 2n. Thus, it is ideal to shrink this space so that fewer samples are needed to train the NN(s) and so that the NN(s) is(are) less complex,i.e., of smaller scale. One way to reduce the optimization input space dimension by one is to exploit the homogeneity-of-order-one nature of the problem. For this homogeneity to hold, two assumptions are made: (a) the building can be safely controlled so that its states remain bounded and inequal- ity (4.6g) will never be active; (b) the control device is capable of producing the required control forces so that force constraints u min ≤ u k ≤ u max will never be active (otherwise, a force saturation can be introduced, but is never required in the numerical examples herein). Because homogeneity- of-order-one holds with these assumptions, the sampling space for x 0 may be contracted fromR 2n to only those on the unit hypersphere because the control force for any other state vector can be scaled based on the state vector’s distance from the origin. Further, without loss of generality, the entire unit hypersphere need not be sampled, but only half of it (as the reflection across the state space origin will give control forces of identical magnitudes but opposite signs); thus, only v rel 0,1 = ¯ v T x 0 ≥ 0 is sampled, with matrix subdivision V≡ [¯ v ¯ V T ] T to extract the first row of V. As proven in Appendix E, the optimal strategy for the opposite sign is ∆ ∆ ∆| v rel 0 = 1 m(p− 1)× 1 − ∆ ∆ ∆| − v rel 0 and θ θ θ| v rel 0 =− θ θ θ| − v rel 0 , where 1 (·)× (·) is a matrix of ones of the given size. Thus, the number of possible S I strategies is 2 mp− 1 , where m is the dimension of the force vector and p is the MPC 86 prediction horizon, though many of these strategies are unlikely to occur (the “–1” in the power is because δ 0,1 is always 1 because the first damper’s relative velocity v rel 0,1 is chosen positive for training). Note that 2 mp− 1 is far smaller than the corresponding number of S II strategies, which is 2 mp− 1 3 2mp . 4.2.2 Sampling in the input parameter space In the previous section, the sampling space for x 0 is limited to x 0 |x 0 ∈R 2n , ∥x 0 ∥ 2 = 1, ¯ v T x 0 ≥ 0 . One can quite easily sample (uniformly) on this half unit hypersphere: choose ψ ψ ψ to be a 2n× 1 vector of independent standard unit normal random variables, and then use sample x 0 =ψ ψ ψ∥ψ ψ ψ∥ − 1 2 sgnψ n+1 (discarding, of course, an all zero vector if it were to occur). However, the state vector during simulations will not contain independent and identically-distributed states, but rather having relative magnitudes roughly similar to the closed-loop state vector covariance matrix (the response is, of course, not exactly Gaussian, but the second order statistics are sufficient to improve the distribution of samples). Thus, because the proposed algorithm seeks to mimic the hMPC, except for faster online execution, the offline NN training can use samples scaled according to an approximation of the covariance matrix of hMPC-controlled state response to a Gaussian white noise excitation. That said, the sampling need not follow this covariance exactly (indeed, one could even use a response controlled by CLQR or another suitable strategy), so one may use an approximate state covariance matrix from the response to some earthquake ground motion, and the time duration of the hMPC simulation need not be infinitely long (but also should not be so short that the covariance is poorly approximated). Thus, samples for NN training are generated through the following steps: 1. Run the hMPC algorithm on the to-be-controlled structure while it is subjected to a real- ization of Gaussian white noise excitation for 200 seconds (the numerical examples herein use a zero-order-hold with the sampling time 0.02 s; thus, the hMPC runs for 10,000 time steps, which was found to be sufficient for the numerical examples herein), and estimate the covariance matrixX of the state vector response. 87 2. For each of the N samples for NN training, generate system state sample x 0 = ˆ ψ ψ ψ∥ ˆ ψ ψ ψ∥ − 1 2 sgn ˆ ψ n+1 , where ˆ ψ ψ ψ = L X ψ ψ ψ is a Gaussian random vector with covariance matrix X= L X L T X , and run the hMPC algorithm with that initial state to find the corresponding strategyS (eitherS I orS II ), providing (x 0 ,S) as an input/output pair for NN training. 4.2.3 Breaking a large NN into smaller NNs In this study, the number of strategies M can be huge, growing exponentially with the number of MPC horizon steps p. For the 5DOF shear building numerical example in Section 4.4,M≃ 13,400 different strategies are found among seven million samples generated for NN training. Because the number of neurons in the output layer must matchM, the number of strategies, and the number of neurons in the last hidden layer, just before the output layer, must be at least M for good NN prediction performance (Zhang et al., 2018), large M means that the scale of the NN is huge, rendering training the NN generally infeasible. The scale of a NN classifier can be reduced by using the “codeword” approach (Zhang et al., 2018), using the intersection of predictions from P smaller-scale NNs — which can be trained offline in parallel and run online to make predictions in parallel — to uniquely predict the correct strategy. Instead of labeling each unique strategy with an integers∈[0,M− 1] and one-hot encode (Harris and Harris, 2010) it, strategy S s can be labeled with a set of smaller integer values, using one of two approaches. If no ordering of these integers is considered, choose a set{λ 1 ,λ 2 ,...,λ P } of relatively prime integersλ j > 1 (i.e., fori̸= j,λ i andλ j have no common positive factors aside from 1) such that∏ P j=1 λ j ≥ M; because theλ j values are relatively prime, strategyS s is uniquely defined by the set {s 1 ,s 2 ,...,s P } where s j =s modλ j . An alternate approach (not used herein) represents s as the number s P ...s 2 s 1 , with digits s j =⌊s/∏ j− 1 i=0 λ j ⌋ modλ j of base λ j > 1 (not necessarily relatively prime), where⌊·⌋ is the “floor” operator (largest integer that is no more than its argument; i.e., ⌊a/b⌋ represents integer division), and λ 0 ≡ 1; if all λ i =λ, i > 0, then this representation is a normal base λ numbering system. With either approach, the j th smaller-scale NN classifier is trained to predict the probabilities of each of the λ j possibly values{0,1,...,λ j − 88 1} of s j so that the P individual predictions of most likely s j together uniquely correspond to the optimal strategy S s . Because λ j ≪ M, the scale of each NN classifier is much smaller than the original single large NN classifier; further, both offline training and online prediction using these smaller NNs can be performed in parallel. If the relatively prime numbers λ j are chosen so that∏ P− 1 j=1 λ j <M≤ ∏ P j=1 λ j , then a correct strategy prediction requires that all P NNs simultaneously predict correctly, which may be chal- lenging because, in practice, no trained NN is perfect. One way to relax this concern is to use γ j ≥ 1 most likely predictions from the j th NN, instead of using only the single most likely predic- tion, leading toΓ=∏ P j=1 γ j possible strategy predictions, and computing the corresponding cost for each; then the strategy leading to the smallest cost will be taken as the correct/best strategy. A second approach would be using additional relatively prime numbers so that∏ P− L− 1 j=1 λ j <M≤ ∏ P− L j=1 λ j for some positive integerL, leading to as many asΓ ′ = P C P− L =P!/[(P− L)!L!] possible strategy predictions, all of which would need to be evaluated as in the previous approach. Although both approaches could be used together to promote the accuracy of the NN strategy prediction, the second is not explored in this study. 4.3 Real-time implementation In Wang and Boyd (2010), a few techniques were adopted to computationally accelerate the QP so- lution inherent in the fMPC, including a primal barrier method, an infeasible start Newton method, fast computation of the Newton step, a warm-start technique, using fixed κ > 0 — the coefficient of the log barrier terms associated with the inequalities in (4.10b), also called the barrier param- eter in Wang and Boyd (2010) — and using a fixed iteration limit for the primal barrier method. Most of these techniques or their variants are applicable herein except for using fixed κ; note that although Wang and Boyd (2010) chose to fix κ to accelerate the algorithm, they also note that κ typically varies in similar algorithms. 89 The applications in Wang and Boyd (2010) used a fixed κ because the state and the control force ranges were small and because their choice of a warm-start technique requires it — neither of which is the case herein. The warm-start herein (Section 4.3.3) does not require a constant κ. Further, ifκ is fixed and too large, the control performance is good when the cost in (4.6a) is large, but becomes overly sub-optimal when the cost is small; in contrast, algorithm convergence within a fixed iteration limit is problematic when the cost is large and κ is fixed and too small, leading to poor control performance. Thus, avoiding a fixed κ is both reasonable and promotes greater accuracy. Another feature of the proposed algorithm, unlike that in Wang and Boyd (2010), is to dis- card the occasional non-converged fMPC solution, which cannot be trusted to satisfy the equality constraints and/or optimality conditions, and to instead fall back to an easily computed alterna- tive control strategy (the CLQR herein) for that time step. To implement this alternative, the C-language fMPC implementation is augmented with an additional output indicating whether a solution converged. This remainder of this section focuses on (a) adaptingκ at each time step, (b) applying a variant of the warm-start technique, (c) bounding the sub-optimality of the final solution, (d) formulating the fMPC algorithm, and (e) implementing the proposed algorithm in real-time. 4.3.1 Adaptingκ at each time step According to Wang and Boyd (2010) and Cauligi et al. (2020), the cost of the optimal solution derived from the primal barrier method is suboptimal by up to κ times the number of effective inequality constraints, which herein is(4mp− 2m)/2=m(2p− 1): ignoring the state bounds (4.6g) because they will not be active in practice given that they are chosen to be of large magnitude; minus 2m because cross-damper velocity bounds (4.6c)–(4.6d) are automatically satisfied at the current time k = 0; and divided by two because, given a value of δ δ δ k and assuming the cross- damper velocity bounds are chosen of sufficiently large magnitude, either the lower bound (4.6c) or the upper (4.6d) will never be active — and similar for inequality constraints (4.6e)–(4.6f) on 90 u k ). In other words, if the NN predicts binary vector ˆ ∆ ∆ ∆ 1 with corresponding theoretical minimal cost J ∗ 1 = min X,U J([X T U T ˆ ∆ ∆ ∆ T 1 ] T ;x 0 ), but the fMPC solver instead finds the suboptimal cost ˆ J 1 = min X,U J([X T U T ˆ ∆ ∆ ∆ T 1 ] T ;x 0 )≥ J ∗ 1 , then the level of sub-optimality∆J 1 = ˆ J 1 − J ∗ 1 ≤ m(2p− 1)κ. Thus, to maintain a constant level of relative sub-optimality — i.e.,∆J 1 /J ∗ 1 ≈ constant or, nearly equivalently assuming∆J 1 is small,∆J 1 / ˆ J 1 ≈ constant — it is reasonable to let barrier weight κ be proportional to J ∗ 1 or ˆ J 1 . However, this means that the fMPC’s input κ value depends on the fMPC result, which would require a costly algorithm that iteratively solves fMPC solutions. One approach is to choose an initial value of κ, find the optimal cost ˆ J 1 through the fMPC solution, and then run the fMPC just one more time; this would be faster than a full iterative solution but still requires two fMPC solutions. Instead, this study proposes approximating J ∗ 1 using the corresponding cost J CLQR = J(ξ ξ ξ CLQR ;x 0 ) that is computed assuming a CLQR algorithm is used to compute the control forces over the next p time steps. The CLQR algorithm uses the previously-defined LQR control gain K LQR but clips the control force to zero if nondissipative;i.e., the p-step ahead prediction of the CLQR-controlled system is u d k =− K LQR x CLQR k u CLQR k,j =u d k,j H(u d k,j v CLQR k,j ), j= 1,2,...,m x CLQR k+1 = Ax CLQR k + B u u CLQR k        , k= 0,1,...,p− 1 ξ ξ ξ CLQR =[X CLQRT U CLQRT ∆ ∆ ∆ CLQRT ] T X CLQR =[x CLQR 1 T x CLQR 2 T ... x CLQR p T ] T U CLQR =[u CLQR 0 T u CLQR 1 T ... u CLQR p− 1 T ] T (4.9) where u d k =[u d k,1 u d k,2 ... u d k,m ] T is the “desired” force vector, and ∆ ∆ ∆ CLQR is ignored as it is not needed for computing the resulting cost (4.5). Preliminary numerical simulations demonstrated that the hMPC algorithm can result in a cost reduction of up to 50% relative to that using the CLQR algorithm; thus,J CLQR is on the same order of magnitude asJ ∗ 1 and ˆ J 1 . Given thatκ can take on a wide range of values of the correct magnitude 91 and still lead to sufficiently accurate solutions (Wang and Boyd, 2010), the choice of κ need not be perfect and using aκ proportional toJ CLQR should provide solutions that are sufficiently accurate. 4.3.2 Bounding the sub-optimality By adjusting barrier coefficient κ, the sub-optimality is properly bounded when ˆ ∆ ∆ ∆ 1 ≃∆ ∆ ∆ ∗ . However, even well-trained NNs will occasionally and inevitably provide incorrect predictions (prediction accuracy can be close, but not strictly equal, to 100%). To maintain consistently good control performance, the sub-optimality resulting from these incorrect predictions must be bounded. To reduce the effect of incorrect NN prediction, this study proposes using not just the NNs’ prediction of the single most likely strategy ˆ ∆ ∆ ∆ 1 , but rather to evaluate the costs of theΓ> 1 most likely strategies — each of which requires a fMPC solution ˆ J r = min X,U J([X T U T ˆ ∆ ∆ ∆ T r ] T ;x 0 ), r = 1,2,...,Γ — as well as CLQR cost J CLQR , and choose the control force at each time step corresponding to the smallest cost. This approach ensures that the proposed algorithm performs at least as well as the CLQR algorithm regardless of the performance of the NNs. When the NN predictions are extremely accurate, the proposed real-time control algorithm will outperform the CLQR approximately just as much as the hMPC; in the worst case, when the NN predictions are undesirable, the proposed algorithm will perform very closely to the CLQR algorithm. Thus, the overall sub-optimality of the proposed algorithm is bounded by the CLQR algorithm. (While the CLQR algorithm is used herein, any other easily computed control algorithm, preferably with an analytical solution, could be used.) For the numerical examples demonstrated in Section 4.4, given the excitations in Table 4.2, the percentage of time the NN does not output the optimal strategy ranges from less than 1% to 17%, and the percentage of time the CLQR control force is adopted (instead of the fMPC control force) ranges from less than 1% to 19% — note that, with the optimal strategy, the (approximate) fMPC performs better than the CLQR most of the time; when it does not, the fallback CLQR control force generally does not deviate much from the optimal semiactive control force, so the 92 proposed algorithm performs almost identical to the hMPC, as shown in the numerical examples (Section 4.4). 4.3.3 Warm-start technique For MPC, a “warm start” refers to initiating the optimization problem at the current time step using most of the optimal solution from the previous time step. Let X ∗ =[x T 1 x T 2 ... x T p ] T and U ∗ =[u T 0 u T 1 ... u T p− 1 ] T be the solutions from the previous time step; one warm-start is to use [x T 2 ... x T p 0 1× 2n ] T and [u T 1 ... u T p− 1 0 1× m ] T as the initial guesses for the current time step’s X and U, respectively (Wright, 2000; Yildirim and Wright, 2002). Three drawbacks arise herein for this warm-start: first, u 1 might fail to satisfy the dissipativity constraints; second, letting the last n and m entries of the initial guesses of X and U take on 0 n× 1 and 0 m× 1 , respectively, is purely heuristic with no mathematical justification (arguably one could instead use x p again or Ax p + B u u p− 1 , and u p− 1 again, but that last term has little effect); third, the sub-optimality of this initial guess cannot be estimated, so this warm-start technique cannot be used to estimate the initial κ. Thus, the typically used warm-start technique is not applicable herein. However, as noted in Wang and Boyd (2010), a warm-start technique can take other forms. In this study, at each time step, x CLQR k+1 and u CLQR k ,k= 0,1,...,p− 1, used to generate the costJ CLQR are taken as the initial guess for minimizingJ([X T U T ˆ ∆ ∆ ∆ T r ] T ;x 0 ). Note that although x CLQR k+1 and u CLQR k ,k= 0,1,...,p− 1, satisfy (4.6b)–(4.6g) for some∆ ∆ ∆, they may not satisfy these inequalities for the NN predicted ˆ ∆ ∆ ∆ 1 , in which case the elements of v rel k CLQR = Vx CLQR k and u CLQR k causing the (4.6b)–(4.6g) violation(s) will be set to zero(s); note: when making these adjustments, state equation (4.6b) is not required to initially hold because the infeasible start Newton method adopted herein to solve the fMPC problem does not require the initial solution to satisfy equality constraints. To evaluate the improvement provided to the proposed algorithm (NN+fMPC) by using this warm-start technique instead of the “cold start” initial guesses of X= 0 and U= 0, both are used on the first three numerical examples introduced in Section 4.4 when the structure is subjected 93 to a realization of Gaussian white noise excitation. The fMPC parameters (e.g., maximum itera- tion limit) are tuned so that the costs computed with the proposed algorithm differ by less than 5% from those with the hMPC algorithm. Table 4.1 summarizes the average computation time for running the hMPC as well as the proposed algorithm with both starting approaches, and the speedups achieved by the proposed algorithm. The warm-start technique proposed herein provides no significant improvement for the hMPC computation time (because of the branch-and-bound logic inherent in the baseline hMPC algorithm solver); in contrast, this warm-start speeds up the proposed NN+fMPC algorithm by 1.3–2.2 times compared with the zero cold start. Further, from Table 4.1, it is reasonable to expect a speedup of around one order of magnitude to be achieved by using the proposed algorithm (compared against the hMPC algorithm). For the SDOF system, the Table 4.1: Computation time (in milliseconds) for running one step of the hMPC and the NN+fMPC algorithms (one-controllable-damper examples) Cold start (X= 0, U= 0) Warm-start using CLQR solution Comp. Time (ms) Time Ratio (Speedup) Comp. Time (ms) Time Ratio (Speedup) hMPC NN+fMPC hMPC NN+fMPC SDOF 37.7 1.9 19.8 34.3 1.4 24.5 3DOF 64.6 15.3 4.2 56.5 6.9 8.2 5DOF 102.0 14.5 7.0 98.0 10.7 9.2 speedup provided by the NN+fMPC seems to be much greater than for the 3DOF and 5DOF exam- ples; however, this is an artifact of the SDOF exploring only the first Γ= 2 most likely strategies, thus requiring two fMPC solutions per time step, whereas the 3DOF and 5DOF examples solve Γ= 4 fMPC solutions; if this difference were eliminated, the SDOF NN+fMPC speedup would be of the same order of magnitude as the other two examples, particularly with the proposed warm start. Additionally, because multiple fMPC solutions are evaluated, there is additional room for a 2–4-fold increase in speedup if the fMPC solutions are pursued in parallel. 4.3.4 Formulating the fMPC problem at each time step The C implementation of the fMPC introduced in Wang and Boyd (2010) handles a separable purely quadratic cost function with box constraints; this C code is generalized in two ways for this 94 study. First, as mentioned previously, a convergence flag is output. Second, to accommodate the off-diagonal N matrix in the cost function, a change of variable is used: let u ′ k = u k + R − 1 N T x k and U ′ =[u ′ 0 T u ′ 1 T ... u ′ T p− 1 ] T , A ′ = A− R − 1 N T and Q ′ = Q− NR − 1 N T . Then, the optimization problem can be written min X,U ′ J(ξ ξ ξ ′ ;x 0 )= x T 0 Q ′ x 0 + min ξ ξ ξ ′ ξ ξ ξ ′ T Hξ ξ ξ ′ (4.10a) subject to ξ ξ ξ ′ =[X T U ′ T ˆ ∆ ∆ ∆ T r ] T and Gξ ξ ξ ′ = b, ¯ Fξ ξ ξ ′ ≤ ¯ h (4.10b) where the final two conditions impose the variable-changed state equation (4.6b) and the con- straints (4.6c)–(4.6g). The matrices in (4.10) are detailed in Appendix F. 4.3.5 Steps for real-time implementation of the proposed algorithm Given current state vector x 0 , the following steps should be taken to find the control forces u 0 for the current step using the proposed real-time semiactive control algorithm. 1. Based on x 0 , run the CLQR in the background for p time steps and store J CLQR , X CLQR and U CLQR . 2. For then× 1 input vector x 0 ∥x 0 ∥ − 1 2 sgnx 0,n+1 , let theP trained NNs make predictions (while not implemented for the numerical examples herein, these NNs can predict in parallel). 3. If the j th (1≤ j≤ P) NN providesγ j possible predictions, thenΓ=∏ P j=1 γ j codewords are available; for each r = 1,2,...,Γ, check a predefined look-up table for the corresponding strategy ˆ ∆ ∆ ∆ r . 4. If v rel 0,1 = ¯ v T x 0 < 0, let the value of each ˆ ∆ ∆ ∆ r be replaced by 1 m(p− 1)× 1 − ˆ ∆ ∆ ∆ r (i.e., its comple- ment),r= 1,2,...,Γ. 5. For each strategy ˆ ∆ ∆ ∆ r , formulate the fMPC optimization problem as in Section 4.3.4. 6. Determineκ based onJ CLQR . 95 7. For each strategy ˆ ∆ ∆ ∆ r , run the fMPC (while not implemented for the numerical examples herein, multiple fMPC algorithms can run in parallel) with initial solution guess provided by the warm-start technique: useξ ξ ξ CLQR directly if it satisfies (4.10b) or adjust it to so that it is in the space defined by (4.6b)–(4.6g). For the fMPC runs that converge, make a record of both ˆ J r ≡ J([X T r U T r ˆ ∆ ∆ ∆ T r ] T ;x 0 ) and the optimal zeroth step control force u 0 . 8. Find the smallest value among { ˆ J 1 , ˆ J 2 ,..., ˆ J Γ ,J CLQR }, selecting the corresponding control force vector to apply to the structure. Note that, if the NNs are used to predict the active constraints as well as the binary variables ˆ ∆ ∆ ∆ r —i.e., strategies of theS II type — then instead of running fMPC in step 7, an optimization problem with only equality constraints would be solved through simple matrix manipulations at very high speed (Bertsimas and Stellato, 2019) (though the effort for generating samples and training the NN is excessive and not implemented for the examples herein). 4.4 Numerical examples In this section, four shear-building numerical examples are used to demonstrate the NN training and the tuning of fMPC parameters, and to evaluate the performance of the proposed algorithm. In the first three examples, one controllable damper is installed between the ground and the first building level, as shown in Figure 4.1, for one-, three-, and five-story buildings, respectively. In the last example, four controllable dampers are installed in stories 1–4, one in each story, of a five-story building (Figure 4.9). In these examples: the displacements are ground-relative; the excitation w l =− ¨ q g l M1 n× 1 where ¨ q g l is the earthquake ground acceleration at timel; the j th column of ˆ B u is [− 1 0 1× (n− 1) ] T if i= 1 and[0 1× (i− 2) 1 − 1 0 1× (n− i) ] T for i> 1; and the state and control force bounds are chosen large enough that they are inactive, and the cross-damper velocity bounds are chosen to be v max i =− v min i =x max n+i − x min n+i− 1 where x min 0 ≡ 0. These examples use a serviceability goal, minimizing the absolute accelerations ¨ e q abs l = ¨ e q l + ¨ q g l 1 n× 1 (plus an optimalρ-weighted penalty 96 on control force) — i.e., the cost function is e J=∆t∑ n t l=0 || ¨ e q abs l || 2 2 +∆t∑ n t l=0 ρ||u l || 2 2 — resulting in weighting matrices Q=[M − 1 K M − 1 C] T [M − 1 K M − 1 C] N=[M − 1 K M − 1 C] T ¯ B u R= ¯ B T u ¯ B u +ρI m . (4.11) In the first three examples, the control forces have reasonable magnitudes and need not be explicitly penalized, soρ is set to zero. The hMPC MIQP problems are solved using the Gurobi ® solver (Gurobi, 2020) — with multi- core parallelism allowed by default) and programmed with its MATLAB ® application programming interface (API). The fMPC QP problems are solved with the C implementation introduced in Wang and Boyd (2010) but adapted with the generalizations explained in Section 4.3. The computing platform used herein is a self-assembled workstation, with an Intel ® i7-10700 eight-core processor clocked at 4.8 GHz (when boosted) and 96 GB of RAM (DDR4 2666 MHz) running MATLAB ® R2020b (MATLAB, 2020) on a 64-bit Windows ® 10 operating system. m n m 2 m 1 Foundation q n q 2 q 1 q g Damper Ground motion k n ,c n k 2 ,c 2 k 1 ,c 1 Figure 4.1: Multi-DOF shear building 97 The earthquake ground motion records used in the simulations, selected from the NGA-West 2 ground motion database (Ancheta et al., 2013) provided by the Pacific Earthquake Engineering Research (PEER) center, are summarized in Table 4.2; earthquakes Nos. 1–10 are far-field earth- quakes, and earthquakes Nos. 11–15/Nos. 16–20 are near-field earthquakes without/with impulses. Table 4.2: Information for Earthquake Records No. Name Year Station Mag. PGA [g] 1 Landers 1992 Yermo Fire Station 7.28 0.152 2 Hector Mine 1999 Hector 7.13 0.264 3 Kocaeli Turkey 1999 Duzce 7.51 0.310 4 Chi-Chi Taiwan 1999 CHY101 7.62 0.313 5 Imperial Valley-02 1940 El Centro Array #9 6.95 0.349 6 Imperial Valley-06 1979 El Centro Array #11 6.53 0.367 7 Kobe Japan 1995 Nishi-Akashi 6.90 0.483 8 Manjil Iran 1990 Abbar 7.37 0.515 9 Northridge-01 1994 Beverly Hills-14145 Mulhol 6.69 0.601 10 Duzce Turkey 1999 Bolu 7.14 0.739 11 Chi-Chi Taiwan 1999 TCU054 7.62 0.140 12 Christchurch New Zealand 2011 Christchurch Botanical Gardens 6.20 0.141 13 Gazli USSR 1976 Karakyr 6.80 0.599 14 Nahanni Canada 1985 Site 1 6.76 0.945 15 Cape Mendocino 1992 Cape Mendocino 7.01 1.494 16 Irpinia Italy-01 1980 Bagnoli Irpinio 6.90 0.190 17 Superstition Hills-02 1987 Parachute Test Site 6.54 0.368 18 Loma Prieta 1989 Saratoga-Aloha Ave 6.93 0.512 19 Niigata Japan 2004 NIGH11 6.63 0.598 20 Chi-Chi Taiwan 1999 TCU065 7.62 0.761 Note: “Mag.” = Moment Magnitude;g= 9.81m/s 2 is the gravitational acceleration. Time-history analyses are performed using MATLAB’s ode45 function with the excitations taken from Table 4.2. The control forces are computed at each time step by the hMPC, the NN+fMPC, and the CLQR algorithms, respectively. A relative tolerance of 10 − 7 and an absolute tolerance of 10 − 9 are used inode45 to ensure the convergence of the peak absolute acceleration metric. At each time step, barrier coefficient κ is determined so that ∆J/J CLQR ≤ 5%, i.e., κ = (5%/(2mp− m))J CLQR , which isJ CLQR /780 for the first three examples and J CLQR /1520 for the fourth 98 example (in which four controllable dampers are used). The iteration limit for the infeasible-start Newton method, which is used to solve the optimization problem formulated by the primal-barrier method (Wang and Boyd, 2010), is set to be 13 for the four examples considered herein, so that pro- posed algorithm (NN+fMPC) performs almost the same as the hMPC algorithm (further increasing this iteration limit does not lead to noticeable control performance improvements). The simulation results include the cost e J, defined in (4.4), and the peak and root-mean- square (RMS) metrics of inter-story drifts, absolute accelerations, and (normalized by the build- ing weight W =∑ n i=1 m i g) controllable damper forces. The i th inter-story drift at time l is de- fined as d l,i = q l,i − q l,i− 1 , i = 1,2,...,n, where q l,0 ≡ 0. The RMS metrics of an n× 1 vec- tor response y l =[y l,1 y l,2 ... y l,n ] T over times l = 0,1,...,n t are y rms i =[n − 1 t ∑ n t l=0 y 2 l,i ] 1/2 and y rms =[n − 1 n − 1 t ∑ n i=1 ∑ n t l=0 y 2 l,i ] 1/2 . The corresponding peak metrics arey peak i = max l |y l,i | andy peak = max i,l |y l,i |. The RMS control force is computed slightly differently, only counting those time steps with nonzero control forces: u rms = h ∑ m j=1 ∑ n t l=0 u 2 l,j I >0 (|u l,j |) i 1/2 h ∑ m j=1 ∑ n t l=0 I >0 (|u l,j |) i − 1/2 for j= 1,2,...,m, where I >0 (·) is an indicator function that is unity when the argument is nonzero (in practice here, a threshold of 10 − 16 N is used) and zero otherwise; thus, u rms is the RMS of the control forces when activated (nonzero). 4.4.1 Single degree of freedom system Following the notations in Figure 4.1, the system parameters of the SDOF system arem 1 = 100Mg, k 1 = 3.948MN/m andc 1 = 62.833kN· s/m, resulting in a natural frequency of 1Hz and a damping ratio of 5%. When constructing the control weight matrix R in (4.5), letρ = 0 so that the control force is not explicitly penalized. This system is the same as what was used in a previous study on semiactive control through hMPC (Elhaddad and Johnson, 2013). The sampling frequency is set to be 50Hz, to which all the earthquake records used in the numerical simulations are downsampled (keep a sample at every 0.02 s and eliminate the others). The hMPC horizon is set to p= 20 time steps, which provides much better performance than the CLQR while not running unnecessarily slow. 99 To train the NN (only one NN is needed for the SDOF system) offline to predict S I =∆ ∆ ∆ ∗ , around 4.8 million system state vectors (x 0 ) are generated so that, by Theorem 9 in McAllester and Schapire (2000), the probability of finding a new strategy with additional samples is less than 0.5%. Solving the MIQP for each x 0 leads to 34 unique strategies. Because the number of possible strategies is small, one moderate scale NN is sufficient to make accurate predictions for the SDOF example. A typical NN training uses the majority (70%–90%) of the samples for training, and the remaining samples for evaluating the NN prediction accuracy. Here, however, the number of samples, chosen to ensure a small likelihood of finding additional strategies with additional samples, is found to be much larger than necessary for high NN prediction accuracy. As a result, for the four examples, 10% of the samples corresponding to each strategy are used for NN training and another 20% used as testing data; this results in a NN prediction accuracy of about 95% when evaluated with the training data and the testing data — which is considered sufficient given the excellent performance detailed in the following sections. The trained NN has one input layer with two neurons, three hidden layers each with 50 neurons, and one output layer with 34 neurons. When making predictions, the bestΓ= 2 predictions will be taken; thus, the fMPC algorithm must be run twice at each time step. Table 4.3 summarizes the average, largest and smallest performance improvements of the NN+fMPC relative to the CLQR, when the structure is subjected to 20 historical earthquake records, in terms of the cost and the peak and RMS metrics of inter-story drifts, absolute acceler- ations, and controllable damper forces; these values are demonstrated in detail in Figure 4.2. For the SDOF RMS responses (including the cost), as detailed in Figures 4.2a–4.2d, the NN+fMPC well outperforms the CLQR without exception, and performs almost as well as the hMPC, which generally uses control force levels smaller than those of the CLQR while delivering better control performance. The peak structure responses (detailed in Figures 4.2e and 4.2f) also decrease with- out exception (easily seen in Table 4.3). The 13% increase in peak control force u peak (see last row of Table 4.3) when the structure is subjected to earthquake ground motion No. 15 (shown in 100 Figure 4.2g) can be traced back to the performance of the hMPC (the target of the NN+fMPC), which also shows an increase of the peak controllable damper force levelu peak by 11%. If it is desired to keep the hMPC (and NN+fMPC) peak control force no higher than that with CLQR, the commanded damper force can be saturated,e.g., at the same level as the CLQR: instead of commanding forceu hMPC l,j , instead command its saturated counterpart min(|u hMPC l,j |,u CLQR peak )sgnu hMPC l,j . A time-history analysis of this saturated system, and the unsaturated one, both subjected to earth- quake No. 15 (Figure 4.3; only the first 15 seconds are shown because the motion is small after that), shows that the hMPC performance improvement in both the cost and response metrics are not compromised, indicating that the peaks in control force have little impact on the overall sys- tem performance and the occasional large values can be clipped without concern, and the superior NN+fMPC performance can be achieved with no greater force demand than that of CLQR. Table 4.3: Average, largest and smallest performance improvements of the NN+fMPC relative to the CLQR in terms of the cost and the peak and RMS metrics of inter-story drifts, absolute accelerations, and controllable damper forces when the structures are subjected to 20 historical earthquake records (one-controllable-damper examples) Average / Largest / Smallest Improvement [%] SDOF 3DOF 5DOF (single device) d rms –30 / –39 / –25 –25 / –36 / –16 –25 / –36 / –14 ¨ q abs rms –12 / –16 / –1 –15 / –25 / –5 –18 / –30 / –4 u rms –42 / –49 / –32 –25 / –38 / –5 –22 / –32 / –10 cost –22 / –30 / –2 –27 / –44 / –10 –31 / –50 / –8 d peak –25 / –33 / –9 –19 / –34 / –2 –18 / –33 / 0 ¨ q abs peak –24 / –32 / –7 –19 / –46 / 18 –22 / –54 / 8 u peak –24 / –33 / 13 –15 / –42 / 12 –15 / –39 / 10 Note: This table shows the percent changes of the values computed from the NN+fMPC relative to those computed from the CLQR; negative numbers denote improvements. 4.4.2 Three degree of freedom system The masses and stiffnesses of a three-story shear building, modeled as a 3DOF system, are de- signed asm 1 =m 2 =m 3 = 100Mg,k 1 = 3.323MN/m,k 2 = 18.988MN/m andk 3 = 9.494MN/m, 101 1 5 10 15 20 0 2 4 hMPC NN+fMPC CLQR (a) RMS drift 1 5 10 15 20 0 0.5 1 hMPC NN+fMPC CLQR (b) RMS abs. accel. 1 5 10 15 20 0 8 16 hMPC NN+fMPC CLQR (c) RMS control force 1 5 10 15 20 0 35 70 hMPC NN+fMPC CLQR (d) Cost Figure 4.2: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (SDOF system) resulting in natural frequencies 0.50Hz, 1.83Hz and 3.42Hz. Rayleigh damping is applied, re- sulting in damping ratios 2.42%, 2.00% and 3.00%, respectively. ρ = 10 − 16 is used in the control weight to avoid numerical issues due to finite precision computation. All other parameters are the same as those of the SDOF system. To train the NNs offline to predict S I =∆ ∆ ∆ ∗ , around 5.5 million system state vectors (x 0 ) are generated (as per Section 4.2.2) so that the probability of finding a 102 1 5 10 15 20 0 9 18 hMPC NN+fMPC CLQR (e) Peak drift 1 5 10 15 20 0 3.5 7 hMPC NN+fMPC CLQR (f) Peak abs. accel. 1 5 10 15 20 0 35 70 hMPC NN+fMPC CLQR (g) Peak control force Figure 4.2: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (SDOF system) (cont.) new strategy with additional samples is again less than 0.5%; solving the MIQP for each x 0 leads to 5,224 unique strategies. Samples for training are chosen in the same way as in the SDOF example. Each of the two trained NNs has one input layer with six neurons, three hidden layers with 200 neurons, and one output layer with 73 and 79 neurons, respectively. When making predictions, the best two predictions from each NN will be taken; thus, the fMPC algorithm must be runΓ= 2· 2= 4 times (or, occasionally, fewer if some predictions do not correspond to valid strategies) at each time step. In terms of RMS responses (including the cost), as summarized in Table 4.3 and detailed in Figures 4.4a–4.4d, the NN+fMPC outperforms the CLQR without exception. The peak inter- story drifts (shown in Figure 4.4e) always decrease, just like in the SDOF example. The peak 103 0 5 10 15 0 8 16 NN+fMPC without saturation NN+fMPC with saturation CLQR (a) Inter-story drifts 0 5 10 15 0 3.5 7 NN+fMPC without saturation NN+fMPC with saturation CLQR (b) Absolute accelerations 0 5 10 15 0 35 70 NN+fMPC without saturation NN+fMPC with saturation CLQR (c) Normalized control force Figure 4.3: Comparisons of the inter-story drifts, the absolute accelerations and the control forces computed by the NN+fMPC with and without force saturation as well as by the CLQR (SDOF system subjected to earthquake No. 15) absolute accelerations ¨ q abs peak (shown in Figure 4.4f) decrease except when the structure is subjected to earthquake No. 12, in which case ¨ q abs peak is up by 18%. Here, the NN+fMPC approach is simply following the hMPC, which sometimes results in a similar or modestly larger ¨ q abs peak (also 18% increase when the structure is subjected to earthquake No. 12) compared with that of the CLQR; thus, this is not a problem with the proposed NN+fMPC algorithm. Further, because the objective 104 is to minimize the RMS responses, some increases in the peak responses are possible. As shown in Figure 4.4g, small increases inu peak (12% when the structure is subjected to earthquake No. 18) occur when the hMPC uses peak forces similar to, or slightly larger than, CLQR. When the larger control forces commanded by the NN+fMPC are saturated to the same level as CLQR’su peak when the controlled structure is subjected to the same excitation, the other performance metrics are not visibly compromised, which was also seen in the SDOF example. 4.4.3 Five degree of freedom system with one controllable damper This system is designed to model a realistic five-story base-isolated building with masses m 1 = 225Mg, m 2 = m 3 = m 4 = m 5 = 150Mg, and stiffnesses k 1 = 6.546MN/m, k 2 = k 3 = k 4 = 37.730MN/m and k 5 = 16.365MN/m, resulting in natural frequencies 0.50Hz, 1.67Hz, 2.79Hz, 4.15Hz and 5.32Hz. Rayleigh damping is applied, resulting in damping ratios 1.58%, 2.00%, 3.00%, 4.31% and 5.45%, respectively. All other parameters are the same as those of the 3DOF system. 4.4.3.1 Control design and performance evaluation To train the NNs offline to predict S I =∆ ∆ ∆ ∗ , around 7.1 million system state vectors (x 0 ) are gener- ated (as per Section 4.2.2) so that the probability of finding a new strategy with additional samples is again less than 0.5%. Solving the MIQP for each leads to 13,403 unique strategies. Samples for training are chosen in the same way as in the previous two examples. Each of the two trained NNs has one input layer with ten neurons, five hidden layers each with 200 neurons, and one output layer with 113 and 127 neurons, respectively. When making predictions, the best two predictions of each NN will be taken; thus, the fMPC algorithm must be run Γ= 2· 2= 4 times at each time step (or, occasionally, fewer if some predictions do not correspond to valid strategies). 105 1 5 10 15 20 0 3 6 hMPC NN+fMPC CLQR (a) RMS drift 1 5 10 15 20 0 2.5 5 hMPC NN+fMPC CLQR (b) RMS abs. accel. 1 5 10 15 20 0 6 12 hMPC NN+fMPC CLQR (c) RMS control force 1 5 10 15 20 0 550 1100 hMPC NN+fMPC CLQR (d) Cost Figure 4.4: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (3DOF system) In terms of RMS responses (including the cost), as summarized in Table 4.3 and detailed in Figures 4.5a–4.5d, the NN+fMPC outperforms the CLQR without exception. In terms of peak re- sponses, the NN+fMPC leads to decreases in peak drifts (shown in Figure 4.5e) without exception. The NN+fMPC generally leads to the decreases of peak accelerations ¨ q abs peak (shown in Figure 4.5f), except for the small increase when the structure is subjected to earthquake No. 7 (< 1% increase), 106 1 5 10 15 20 0 10 20 hMPC NN+fMPC CLQR (e) Peak drift 1 5 10 15 20 0 10 20 hMPC NN+fMPC CLQR (f) Peak abs. accel. 1 5 10 15 20 0 40 80 hMPC NN+fMPC CLQR (g) Peak control force Figure 4.4: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (3DOF system) (cont.) No. 9 (4% increase), No. 11 (1.5% increase), No. 12 (8.3% increase) or No. 15 (4.7% increase) — all because the hMPC either performs similarly to [earthquake No. 11 (10.5% decrease) or No. 15 (8.3% decrease)] or slightly worse than [earthquake No. 7 (2.7% increase), No. 9 (< 1% increase) or No. 12 (12.0% increase)] the CLQR in these cases. Again, because the objective is to minimize the RMS responses, some increases in the peak responses are possible. Further, the increases of ¨ q abs peak when the structure is controlled by either the hMPC or the NN+fMPC usually occur on the level above the isolation layer because of the swift changes in control force level when going from one time step to another. In reality, when the control force is commanded to change gradually within the∆t= 0.02second time step, these acceleration increases will be alleviated. 107 As shown in Figure 4.5g, small increases in peak force u peak occur when the hMPC uses peak forces slightly larger than [earthquake No. 1(9.1% increase) or No. 7 (9.3% increase)] or similar to [earthquake No. 11 (< 1% decrease)] CLQR. When the larger control forces commanded by the NN+fMPC are saturated to the same level as CLQR’s u peak , the other metrics’ performance improvements are, again, not visibly compromised. Further, the peak control forces are generally less than 20% of the total weight of the building except when the structure is subjected to earth- quake No. 20, in which case a large-amplitude impulse leads to a huge control demand — very challenging for all control algorithms to handle. To demonstrate sample time histories, the isolation-layer drifts, the roof absolute accelerations and the control forces (normalized by the building weight), when the structure is subjected to the 1940 El Centro earthquake ground motion record (earthquake No. 5), are shown in Figure 4.6, in which it can be easily observed that the NN+fMPC resembles the hMPC and well outperforms the CLQR. Aside from the earthquake ground motions listed in Table 4.2, the 5DOF building model was also subjected to a Gaussian white noise excitation as well as Kanai-Tajimi stationary models of earthquake-induced ground accelerations (Soong and Grigoriu, 1993; Johnson et al., 1995, 1997; Ramallo et al., 2002) with the fundamental frequency of the ground model ranging from 1Hz to 5Hz. In all of these cases, the performance improvements of the NN+fMPC relative to the CLQR are consistent with the historical earthquakes, so the performance metrics of these stochastic responses are omitted for brevity. 4.4.3.2 Robustness of Control Performance Controllable damper forces are calculated by the control algorithms based on the best knowledge of the structural parameters. However, the nominal structural parameters adopted in the control design may differ somewhat from the real/exact structural parameters — large or small, a gap between the nominal structural model and the real structure is inevitable. Thus, it is essential to 108 1 5 10 15 20 0 5 10 hMPC NN+fMPC CLQR (a) RMS drift 1 5 10 15 20 0 1.5 3 hMPC NN+fMPC CLQR (b) RMS abs. accel. 1 5 10 15 20 0 4 8 hMPC NN+fMPC CLQR (c) RMS control force 1 5 10 15 20 0 250 500 hMPC NN+fMPC CLQR (d) Cost Figure 4.5: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with one controllable damper) study the deterioration of control performance when small discrepancies between the real and the nominal structural parameters exist. To study the robustness of control performance, 3 n+2 = 2,187 system parameter combinations are generated, in each of which the actual stiffness parameters k i (i= 1,2,...,n; n= 5 herein) 109 1 5 10 15 20 0 20 40 hMPC NN+fMPC CLQR (e) Peak drift 1 5 10 15 20 0 6 12 hMPC NN+fMPC CLQR (f) Peak abs. accel. 1 5 10 15 20 0 20 40 hMPC NN+fMPC CLQR (g) Peak control force Figure 4.5: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with one controllable damper) (cont.) takes on a value from{0.9k nom i , 1.0k nom i , 1.1k nom i }, where k nom i is the i th nominal stiffness pa- rameter used to design the control, and the Rayleigh damping C=α R M+β R K coefficients α R and β R of the actual structure take on values chosen from {0.8α nom R , 1.0α nom R , 1.2α nom R } and {0.8β nom R , 1.0β nom R , 1.2β nom R }, where α nom R and β nom R are the corresponding nominal Rayleigh damping parameters used for control design. In Figure 4.7, each case denotes a unique combina- tion of stiffness and damping parameters and is assigned a number from 1 to 2,187;e.g., case No. 1 corresponds to the parameter combination with all actual stiffnesses 10% below and damping pa- rameters 20% below the nominal values used in control design; case No. 1094 corresponds to the parameter combination with all nominal values; and case No. 2187 corresponds to the parameter 110 0 10 20 30 0 10 20 hMPC NN+fMPC CLQR No control (a) First-story inter-story drifts 0 10 20 30 0 2 4 hMPC NN+fMPC CLQR No control (b) Top-story absolute accelerations 0 10 20 30 0 6 12 hMPC NN+fMPC CLQR (c) Normalized control force Figure 4.6: Comparisons of the first story inter-story drifts, the top story absolute accelerations and the control forces computed by the hMPC, the NN+fMPC and the CLQR (5DOF building subjected to earthquake No. 5) combination with stiffnesses 10% above nominal and Rayleigh parameters 20% above nominal used for control design. Time-history analyses are performed for systems parameterized with each parameter combina- tion, subjected to earthquake No. 5 (1940 El Centro) and controlled by the CLQR, the NN+fMPC or the hMPC algorithm (2,187· 3= 6,561 time-history analyses in total). The costs computed 111 k 1 : 0.9k nom 1 1.0k nom 1 1.1k nom 1 ··· ··· ··· ··· ··· ··· k 2 : 0.9k nom 2 1.0k nom 2 1.1k nom 2 ··· ··· ··· ··· ··· ··· k 3 : 0.9k nom 3 1.0k nom 3 1.1k nom 3 ··· ··· ··· ··· ··· ··· k 4 : 0.9k nom 4 1.0k nom 4 1.1k nom 4 ··· ··· ··· ··· ··· ··· k 5 : 0.9k nom 5 1.0k nom 5 1.1k nom 5 ··· ··· ··· ··· ··· ··· α R : 0.8α nom R 1.0α nom R 1.2α nom R ··· ··· ··· ··· ··· ··· β R : 0.8β nom R 1.0β nom R 1.2β nom R Case No.: ··· ··· ··· ··· ··· 1 1093 1094 1095 2187 ··· ··· Figure 4.7: Generation of 3 7 = 2,187 parameter combinations for robustness study (note: the side- branches growing from the side nodes, which are omitted with “··· ”, all grow to the bottom level in the same pattern as the central branch growing from the central node on the same level as the side nodes) when the structure is controlled by the NN+fMPC and the hMPC algorithms are normalized by that computed when the same structure is controlled by the CLQR, and displayed in Figure 4.8. It can be observed that the robustness in terms of control performance for the NN+fMPC or the hMPC algorithms are essentially the same. Further, the NN+fMPC and the hMPC algorithms consistently outperform the CLQR algorithm. 4.4.4 Five degree of freedom system with four controllable dampers The aforementioned examples are controlled by one device installed between the ground and the first level of the building. However, depending on the dynamic characteristics and the performance 112 1 729 1458 2187 0.1 0.5 1 NN+fMPC relative to CLQR hMPC relative to CLQR Figure 4.8: Cost when the system is controlled by the NN+fMPC and the hMPC algorithms nor- malized by the corresponding cost when the system with the same parameterization is controlled by the CLQR algorithm requirements of the building, a multiple-controllable-damper design may be necessary. In this example, the building model is a fixed-base 5DOF structure, with one controllable damper in each of the first four stories, as shown in Figure 4.9, with all floor masses 100Mg and all story stiffnesses 48.727MN/m, respectively, resulting in natural frequencies 1.00Hz, 2.92Hz, 4.60Hz, 5.91Hz and 6.74Hz. Rayleigh damping is applied, resulting in damping ratios 2.00%, 2.52%, 3.56%, 4.43% and 5.00%. ρ = 10 − 11 so that the maximum CLQR control force when the building is subjected to 1940 El Centro earthquake (earthquake No. 5) is 11.84% of the building weight, as shown in Figure 4.10g. The number of hMPC horizon steps is chosen to be ten, i.e., p= 10, leading to much better control performance than that of the CLQR. Note that while p= 20 was used in the previous examples, for the example herein, preliminary studies show that using p> 10 does not lead to significant hMPC control performance improvement. Each of the two trained NNs has one input layer with ten neurons, five hidden layers each with 500 neurons, and one output layer with 347 and 349 neurons, respectively. When making predictions, the best two predictions of each NN will be taken; thus, the fMPC algorithm must be run at mostΓ= 2· 2= 4 times during each time step. To train these NNs offline, to predict S I =∆ ∆ ∆ ∗ , 113 m 5 m 4 m 3 m 2 m 1 Foundation q 5 q 4 q 3 q 2 q 1 q 0 Damper 4 Damper 3 Damper 2 Damper 1 Ground motion k 5 ,c 5 k 4 ,c 4 k 3 ,c 3 k 2 ,c 2 k 1 ,c 1 Figure 4.9: Multi-DOF shear building with multiple controllable dampers 114 using enough samples to make the probability of finding a new strategy less than 0.5%, as in the previous examples, would require 40.6 million system state vectors resulting in 248,755 unique strategies, posing a challenge for computer memory during training if only two large NNs are used and (because more NNs may require a larger number Γ of predictions to ensure a low level of incorrect-NN-prediction-induced suboptimality) longer online computation time for the multiple (Γ) serial fMPC runs (a future parallel implementation would permit a larger number of smaller NNs). Instead, the likelihood of missing a strategy is raised to 1%, requiring around 10.1 million system state vectors (x 0 ) that result in 117,295 unique strategies. The performance improvements of the NN+fMPC (standard and with damper forces saturated to the same level as the CLQR with the same excitation) relative to the CLQR are summarized in Table 4.4 and demonstrated in detail in Figure 4.10. For the case without force saturation, all structure performance metrics are improved; the peak control forceu peak on average is comparable to that of CLQR but is modestly larger for some earthquakes. In Table 4.4, for the case with- out force saturation, the significant improvements in ¨q abs rms and ¨ q abs peak indicate that the NN+fMPC is much better at attenuating the accelerations (the main goal of the control design) than the CLQR; the inter-story drifts (see the rows for d rms and d peak ) are additionally modestly reduced. These improvements are not compromised after imposing force saturation, indicating that the NN+fMPC as well as the hMPC achieve better control performance by using the controllable dampers effi- ciently, not by using higher control force levels, which is echoed by the improvements in terms of u rms in both cases (cases without/with force saturation). Figure 4.10 again demonstrates that the NN+fMPC performs similarly to the hMPC as in the previous one-controllable-damper examples. It may be noted that earthquake No. 15 induces unusually large control force levels because it has a PGA of 1.494g, the largest among not only just all the near-field earthquake records without impulse but also all the earthquake records considered in this study. Finally, the computational cost of the NN+fMPC is around 7 ms per time step, 13 times faster than that of the hMPC (taking 95 ms per time step to run), and also faster than what was achieved in the one-controllable-damper 5DOF example, because there are fewer parameters for the fMPC to 115 optimize —(n+m)p=(10+4)10= 140 parameters herein against(n+m)p=(10+1)20= 220 parameters in the previous one-controllable-damper 5DOF example — due to half the number of prediction steps. Table 4.4: Average, largest and smallest performance improvements of the NN+fMPC relative to the CLQR in terms of the cost and the peak and RMS metrics of inter-story drifts, absolute accelerations, and controllable damper forces when the structures are subjected to 20 historical earthquake records (four-controllable-damper example) Average / Largest / Smallest Improvement [%] Without force saturation With force saturation d rms –4 / –7 / –1 –4 / –7 / –2 ¨ q abs rms –18 / –23 / –8 –18 / –23 / –9 u rms –8 / –11 / –4 –8 / –11 / –4 cost –25 / –34 / –12 –25 / –34 / –12 d peak 0 / –6 / 7 0 / –6 / 7 ¨ q abs peak –40 / –68 / –2 –38 / –65 / 3 u peak 0 / –11 / 29 –3 / –11 / 0 Note: This table shows the percent changes of the values computed from the NN+fMPC controlled case relative to those computed from the CLQR controlled case; negative numbers denote improvements. 116 hMPC NN+fMPC CLQR (a) RMS drift hMPC NN+fMPC CLQR (b) RMS abs. accel. hMPC NN+fMPC CLQR (c) RMS control force hMPC NN+fMPC CLQR (d) Cost Figure 4.10: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with four controllable dampers) 117 hMPC NN+fMPC CLQR (e) Peak drift hMPC NN+fMPC CLQR (f) Peak abs. accel. hMPC NN+fMPC CLQR (g) Peak control force Figure 4.10: Comparisons of the costs and the RMS and peak response and control force metrics calculated by the hMPC, the NN+fMPC and the CLQR algorithms using 20 historical earthquake records (5DOF system with four controllable dampers) (cont.) 118 Chapter 5 Conclusions and Directions for Future Research 5.1 Conclusions — Modeling and Model Updating of a Full- Scale Experimental Base-Isolated Building To leverage the data from the expensive full-scale seismic tests conducted on 8 August 2013 at Japan’s E-Defense facility to be used for future investigations on new seismic protective designs, controllable damping strategies, retrofit procedures, and so forth, FEMs were assembled based on the structural design drawings, updated to have approximately the same modal properties as the experimental building, and combined with optimal nonlinear models of the isolation-layer devices to reproduce the responses of the building under earthquake excitations. Insights gained through this study are listed as follows. 1. For base-isolated structures resembling the one in this study, the modes can be separated into isolation modes and superstructure modes. Superstructure linear stiffness changes have minor effects on the frequencies and mode shapes of the isolation modes, and vice versa. Further, the superstructure damping has a negligible effect on the damping ratios of the iso- lation modes. Finally, the superstructure modes are approximately proportionally damped, while the isolation modes are not. 119 2. A mode-matching algorithm was proposed to efficiently match the experimental modes with the corresponding FEM modes, by searching for FEM modes with modal properties resem- bling those of the to-be-matched experimental mode and choosing the one leading to the smallest deviation in modal characteristics. 3. To reduce the computational burden of base-isolated building model updating, a two-step model updating approach was proposed, in which the superstructure linear stiffnesses are first updated so that the superstructure modes’ natural frequencies and mode shapes of the ex- perimental building match those of the FEM; then, the isolation-layer stiffnesses and damp- ing coefficients are updated so that the isolation modes’ natural frequencies, damping ratios and mode shapes match. The modal properties of the unupdated and updated FEMs were compared with those from the experiment; better agreement after model updating demon- strated the proposed approach successful. The matched modes’ effective masses (moments of inertia), relative to the total mass (total moments of inertia), indicated that the matched modes capture the major dynamic properties of the building in all six DOFs. This conclu- sion was echoed by the good match between the simulated time-history responses using the updated FEM and those from the experimental data. 4. A final partially nonlinear FEM is determined by combining the updated superstructure FEM, rubber bearings modeled with directionally uncoupled bilinear stiffness models with a constant damping coefficient in each horizontal direction, and elastic sliding bearings and steel damper pairs modeled with bidirectionally coupled Bouc-Wen models with constant damping coefficients in each horizontal direction. This combined model satisfactorily repro- duced the nonlinear time-history responses of the experimental building and can be used in future studies to predict the responses of the base-isolated building under earthquakes with peak accelerations in a similar range. 120 The updated partially nonlinear FEM will be made open source and, together with an accompa- nying documentation, will be posted on the Natural Hazards Engineering Research Infrastructure (NHERI) repository (see B). 5.2 Conclusions — Cross-Model Cross-Mode Model Updating Method Based on Constrained Total Least Squares The CTLS approach was utilized to solve a set of ill-conditioned linear equations formulated by the CMCM model updating method to provide the magnitudes of change for the to-be-updated structural parameters. The main theoretical contributions include deriving the formulas for com- bining the CTLS approach with the CMCM method, deriving new formulas for model expansion, and adapting the CTLS solution to better fit the CMCM model updating method. A line-fitting example and a numerical 3D frame structure example demonstrated that the CTLS approach per- forms much better than the state-of-the-art TLS approach and somewhat better than the WTLS approach, in terms of robustness to noise. Thus, the CTLS approach is the best option to be used together with the CMCM model updating method. Using the CTLS-based CMCM method, a FEM of a full-scale four-story structure was updated based on the experimentally identified natural fre- quencies, damping ratios and mode shapes, so that its corresponding modal properties resemble those of the experimental building. Practical skills and details necessary to implement the pro- posed method on a base-isolated building (or base-isolated structures of this type) were carefully discussed. The time-history responses of the building to measured shake table accelerations were simulated and compared with those measured in the experiment. The results show that the model updating was successful, which validates the applicability of the CTLS-based CMCM method to real engineering practice. 121 5.3 Conclusions — Real-Time Neural Network based Semiactive Model Predictive Control of Structural Vibrations Proposed herein is a semiactive structural control algorithm that performs almost as well as the hMPC and runs fast enough for real-time implementation. The success of this algorithm is con- tingent upon the skillful selection of the sampling space for the system state vector x 0 , a novel technique to tune the fMPC penalty coefficient κ, a novel application of the warm-start technique, the proposed way of bounding the sub-optimalities (either from inaccurate NN predictions or the QP solver) and advanced programming (efficient manipulations of matrices and exploitations of possible parallelisms). Based on the numerical simulations conducted, the following conclusions can be drawn: 1. By skillfully selecting the sampling space instead of sampling naively, the NN(s) can be trained to accurately predict the binary variables in the hMPC algorithm. 2. By letting NN(s) predict only the binary variables and using multiple NNs to predict the most possible codewords/strategies, the offline NN training and online NN prediction efforts are reduced dramatically, making the proposed algorithm capable of handling large-scale problem arising from structural vibration control. 3. When solving the fMPC QP problem, a fixed penalty coefficient κ is not an ideal option with large variations of the control algorithm input. The proposed approach for adaptingκ leads to solutions with a consistent level of accuracy regardless of the magnitudes of the input. 4. The proposed warm-start technique, which is compatible with the control problem formu- lated herein, can further accelerate the fMPC algorithm. 122 5. To limit the influence of inaccurate predictions from the NN(s), a fast control algo- rithm (e.g., CLQR) that minimizes a similar cost function is required to offer a lower bound on the overall control performance. 6. In general, the proposed algorithm can achieve control performance close to that of the hMPC, which is typically much better than that of the CLQR. 7. The proposed algorithm is fast enough to be applied in real-time to control structural vibra- tions. Even without exploiting possible parallelisms in programming, this has been proven to be true for two realistic 5DOF shear buildings; this enables application of the proposed algorithm to structures with global dynamics that can be approximated by up to about ten DOFs, or to large-scale structures with decentralized control strategies. When parallelisms are fully exploited (making this algorithm around 3–4 times faster) in the future theoretical and experimental studies, the proposed algorithm will be able to accommodate even more complex structure models. Future research that can further reduce the sample generation computational effort (currently, the sample generation process, i.e., solving the MIQP problems given the randomly generated initial state vectors, preferably requires high performance computing facilities) and that of the fMPC are welcomed. 5.4 Directions for Future Research Based on the successful numerical implementation real-time NN based semiactive control, the following works are proposed to be done in the future. (1) As discussed early, the number of samples collected is based on Good-Turing estima- tion theory (McAllester and Schapire, 2000), which says that “The probability of encountering a parameter θ θ θ N+1 corresponding to an unseen strategy S(θ θ θ N+1 ) satisfies with confidence at least 1− β: P(θ θ θ N+1 ∈ R p |S(θ) / ∈ S(Θ Θ Θ N ))≤ G+c p (1/N)ln(3/β); where G= N 1 /N (N 1 is number 123 of distinct strategies that appeared exactly once; N is total number of samples); c= 2 √ 2+ √ 3; Θ Θ Θ N ={θ θ θ 1 , θ θ θ 2 ,...,θ θ θ N } are the samples; S(Θ Θ Θ N )={S 1 , S 2 ,...,S N } are the unique strategies found.”. According to this theory, for a 5DOF building, around seven million samples need to be generated to make sure enough unique strategies are found. This forms the bottleneck for implementing the proposed methods in terms of designing time. One way to alleviate this burden is to start training a smaller neural network after getting a certain amount of samples. Then for each state vector sample thereafter, the optimization problem will be solved both by using the neural network together with the fMPC and by using Gurobi ® to solve it in the brute force way. If the later only offers negligible amount of improvements, then the sample will be discarded; if not, the sample will be collected. Every time when enough new samples are collected, an expanded neural network will be trained based on all current samples. By checking and discarding samples accordingly, the value of N 1 will be reduced for a given number of samples, resulting in the decrease of G for a given number of samples — the probability criteria can be met with fewer samples generated. (2) In general, the parameters of the building (masses, stiffnesses, damping ratios, et al.) can change with the elapse of time. Thus, a control design based on the parameters of the building at a given time point might become sub-optimal later. Thus, when collecting samples for NN training, it is reasonable to make these changing parameter values as part of the inputs so that their changes within certain ranges can be incorporated into the NNs. During the implementation of the proposed method, parameter identification algorithms will be used to find the structural parameters (given the current computational speed, there is time to incorporate the parameter identification process), then the neural network predictions as well as the subsequent fMPC runs will be based on the identified parameters. (3) Under certain circumstances, practical constraints (e.g., the scale of the structure to be con- trolled and the complexity of the control algorithm) will prevent the use of centralized control approach Bakule (2008). Then the control problem will need to be divided into smaller problems with little interdependence. This is a great scenario for implementing the proposed method con- sidering that it has been proved to work well with a 5DOF building. Thus, as part of the future 124 research investigations, the proposed method will be implemented on a benchmark structure in a decentralized control setting. (4) In recent years, “dual control” has become an important topic, in which the control efforts are not only used to achieve desirable control performances but also used to help facilitate the identification of system parameters. According to Filatov and Unbehauen (2004), “a dual control system should have two properties: (i) the system output cautiously tracks the desired reference value; (ii) the control efforts excite the plant sufficiently for accelerating the parameter estima- tion process so that the control quality becomes better in future time intervals.”. Dual control ap- proaches have been developed in conjunction with many popular control methods (LQG controller, pole-placement controller, generalized minimum-variance controller, et al.), by making modifica- tions to the original control strategies. However, it is hard to extend semiactive control strategies that involve clipping (clipped-LQR) to form dual control algorithms because it will stay dormant for a significant amount of time when the control forces are highly nondissipative. However, when using the proposed methods as well as the original hMPC approach, this problem no longer exists, because the control forces do not stay dormant (are not clipped) as often. Thus, formulating a dual control algorithm based on the proposed methods is an interesting future research topic. 125 Bibliography Abatzoglou, T.J., Mendel, J.M., 1987. Constrained total least squares, in: IEEE International Conference on Acoustics, Speech, and Signal Processing, pp. 1485–1488. Abatzoglou, T.J., Mendel, J.M., Harada, G.A., 1991. The constrained total least squares technique and its applications to harmonic superresolution. IEEE Transactions on Signal Processing 39, 1070–1087. Allemang, R.J., Brown, D.L., 1982a. A correlation coefficient for modal vector analysis, in: Pro- ceedings of the 1st International Modal Analysis Conference & Exhibit, Orlando, Florida, USA, 1982. Allemang, R.J., Brown, L.D., 1982b. A correlation coefficient for modal vector analysis, in: Pro- ceedings of the 1st International Modal Analysis Conference, pp. 110–116. Allwein, E.L., Schapire, R.E., Singer, Y ., 2000. Reducing multiclass to binary: A unifying ap- proach for margin classifiers. Journal of Machine Learning Research 1, 113–141. Ancheta, T.D., Darragh, R.B., Stewart, J.P., Seyhan, E., Silva, W.J., Chiou, B.S., Wooddell, K.E., Graves, R.W., Kottke, A.R., Boore, D.M., Kishida, T., Donahue, J.L., 2013. PEER NGA-West2 Database. Technical Report. Pacific Earthquake Engineering Research Center. Avitabile, P., O’Callahan, J.C., 1988. Model correlation and orthogonality criteria, in: Proceedings of the 6th International Modal Analysis Conference, pp. 1039–1047. Bakir, P.G., Reynders, E., De Roeck, G., 2007. Sensitivity-based finite element model updating using constrained optimization with a trust region algorithm. Journal of Sound and Vibration 305, 211–225. Bakir, P.G., Reynders, E., De Roeck, G., 2008. An improved finite element model updating method by the global optimization technique ‘Coupled Local Minimizers’. Computers & Structures 86, 1339–1352. Bakule, L., 2008. Decentralized control: An overview. Annual reviews in control 32, 87–98. Beck, J.L., Au, S.K., Vanik, M.W., 2001. Monitoring structural health using a probabilistic mea- sure. Computer-Aided Civil and Infrastructure Engineering 16, 1–11. Beck, J.L., Katafygiotis, L.S., 1998a. Updating models and their uncertainties. I: Bayesian statis- tical framework. Journal of Engineering Mechanics 124, 455–461. 126 Beck, J.L., Katafygiotis, L.S., 1998b. Updating models and their uncertainties. II: Model identifi- ability. Journal of Engineering Mechanics 124, 463–467. Beck, J.V ., Arnold, K.J., 1977. Parameter estimation in engineering and science. John Wiley & Sons. Bemporad, A., 2002. An efficient technique for translating mixed logical dynamical systems into piecewise affine systems, in: Proceedings of the IEEE Conference on Decision and Control, Las Vegas, NV , USA. pp. 1970–1975. Bemporad, A., Morari, M., 1999. Control of systems integrating logic, dynamics, and constraints. Automatica 35, 407–427. Bertsimas, D., Stellato, B., 2019. Online mixed-integer optimization in milliseconds. arXiv preprint arXiv:1907.02206 . Bertsimas, D., Stellato, B., 2021. The voice of optimization. Machine Learning 110, 249–277. Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press, Cambridge, U.K. Bresler, Y ., Macovski, A., 1986. Exact maximum likelihood parameter estimation of superimposed exponential signals in noise. IEEE Transactions on Acoustics, Speech, and Signal Processing 34, 1081–1089. Brewick, P.T., Johnson, E.A., Sato, E., Sasaki, T., 2018. Constructing and evaluating generalized models for a base-isolated structure. Structural Control and Health Monitoring 25, 2243. Brewick, P.T., Johnson, E.A., Sato, E., Sasaki, T., 2020. Modeling the dynamic behavior of iso- lation devices in a hybrid base-isolation layer of a full-scale building. Journal of Engineering Mechanics 146, 040220127. Bugeat, L.P., Lallement, G., 1976. Methods of matching calculated and identified eigensolutions. Storjnicky Casopis 32, 162–167. Carroll, R.J., Ruppert, D., Stefanski, L.A., Crainiceanu, C.M., 2006. Measurement error in non- linear models: a modern perspective. Chapman and Hall/CRC. Cauligi, A., Culbertson, P., Stellato, B., Bertsimas, D., Schwager, M., Pavane, M., 2020. Learning mixed-integer convex optimization strategies for robot planning and control, in: Proceedings of the IEEE Conference on Decision and Control, Jeju Island, Republic of Korea. pp. 1698–1705. Chen, W., Yuan, Q., Chen, Y ., 2011. Application of constrained total least-squares to cloud point registration. Journal of Geodesy and Geodynamics 31, 137–141. Coleman, T.F., Li, Y ., 1996. A reflective Newton method for minimizing a quadratic function subject to bounds on some of the variables. SIAM Journal on Optimization 6, 1040–1058. Constantinou, M.C., Mokha, A.S., Reinhorn, A.M., 1990. Teflon bearings in base isolation. II: modeling. ASCE Journal of Structural Engineering 116, 455–474. 127 Cornwell, P., Doebling, S.W., Farrar, C.R., 1999. Application of the strain energy damage detection method to plate-like structures. Journal of Sound and Vibration 224, 359–374. Dietterich, T.G., Bakiri, G., 1994. Solving multiclass learning problems via error-correcting output codes. Journal of Artificial Intelligence Research 2, 263–286. Dyke, S.J., Spencer, Jr., B.F., Sain, M.K., Carlson, J.D., 1996. Modeling and control of magne- torheological dampers for seismic response reduction. Smart Materials and Structures 5, 565. Elhaddad, W.M., Johnson, E.A., 2013. Hybrid MPC: An application to semiactive control of structures, in: Conference Proceedings of the Society for Experimental Mechanics Series, New York, NY , USA. pp. 27–36. Escalera, S., Pujol, O., Radeva, P., 2009. Separability of ternary codes for sparse designs of error- correcting output codes. Pattern Recognition Letters 30, 285–297. Fang, X., 2013. Weighted total least squares: necessary and sufficient conditions, fixed and random parameters. Journal of geodesy 87, 733–749. Fang, X., 2015. Weighted total least-squares with constraints: a universal formula for geodetic symmetrical transformations. Journal of geodesy 89, 459–469. Filatov, N.M., Unbehauen, H., 2004. Adaptive dual control: Theory and applications. volume 302. Springer Science & Business Media. Flanigan, C.C., 1998. Model reduction using Guyan, IRS, and dynamic methods, in: spie the inter- national society for optical engineering, SPIE INTERNATIONAL SOCIETY FOR OPTICAL. pp. 172–176. Friswell, M., Mottershead, J.E., 2013. Finite element model updating in structural dynamics. Springer Science & Business Media. Ghani, R., 2001. Using error-correcting codes for efficient text classification with a large number of categories. Technical Report. KDD Lab Project Proposal. URL: http://citeseerx. ist.psu.edu/viewdoc/summary?doi=10.1.1.18.7018. Golub, G.H., 1973. Some modified matrix eigenvalue problems. Siam Review 15, 318–334. Golub, G.H., Van Loan, C.F., 1980. An analysis of the total least squares problem. SIAM journal on numerical analysis 17, 883–893. Guner, U., Jang, H., Realff, M.J., Lee, J.H., 2015. An extended constrained total least-squares method for the identification of genetic networks from noisy measurements. Industrial & Engi- neering Chemistry Research 54, 10583–10592. Gurobi, 2020. Gurobi optimizer reference manual (Version 9.0). Gurobi Optimization, LLC, Beaverton, OR, USA. Harris, D.M., Harris, S.L., 2010. Digital design and computer architecture. Second edition ed., Morgan Kaufmann, Waltham, MA, USA. 128 Heemels, W.P.M.H., De Schutter, B., Bemporad, A., 2001. On the equivalence of classes of hybrid dynamical models, in: Proceedings of the IEEE Conference on Decision and Control, Orlando, FL, USA. pp. 364–369. Housner, G.W., Bergman, L.A., Caughey, T.K., Chassiakos, A.G., Claus, R.O., Masri, S.F., Skel- ton, R.E., Soong, T.T., Spencer, Jr., B.F., Yao, J.T.P., 1997. Structural control: past, present and future. Journal of Engineering Mechanics 123, 897–971. Hu, S.L.J., Li, H., 2007. Simultaneous mass, damping, and stiffness updating for dynamic systems. AIAA journal 45, 2529–2537. Hu, S.L.J., Li, H., Wang, S., 2007. Cross-model cross-mode method for model updating. Mechan- ical Systems and Signal Processing 21, 1690–1703. Hu, S.L.J., Wang, S., Li, H., 2006. Cross-modal strain energy method for estimating damage severity. Journal of Engineering Mechanics 132, 429–437. Ibrahim, S.R., 1993. Correlation and updating methods: finite element dynamic model and vibra- tion test data, in: International Conference on Structural Dynamics Modelling, Test, Analysis and Correlation, pp. 323–347. Imregun, M., Visser, W.J., Ewins, D.J., 1995. Finite element model updating using frequency response function data: I. theory and initial investigation. Mechanical Systems and Signal Pro- cessing 9, 187–202. Jang, J., Smyth, A.W., 2017a. Bayesian model updating of a full-scale finite elementmodel with sensitivity-based clustering. Struct Control and Health Monitoring 24, e2004. Jang, J., Smyth, A.W., 2017b. Model updating of a full-scale FE model with nonlinear constraint equations and sensitivity-based cluster analysis for updating parameters. Mechanical Systems and Signal Processing 83, 337–355. Johnson, E.A., Wojtkiewicz, S.F., Bergman, L.A., 1995. Some experiments with massively parallel computation for Monte Carlo simulation of stochastic dynamical systems. Balkema, Rotterdam, Netherlands. Johnson, E.A., Wojtkiewicz, S.F., Bergman, L.A., Spencer, Jr., B.F., 1997. Observations with regard to massively parallel computation for Monte Carlo simulation of stochastic dynamical systems. International Journal of Non-Linear Mechanics 32, 721–734. Third International Stochastic Structural Dynamics Conference. Li, Y ., Wang, S., Xia, Z., An, W., 2020. An iterative total least squares-based estimation method for structural damage identification of 3D frame structures. Structural Control and Health Mon- itoring . Lieven, N.A.J., Ewins, D.J., 1988. Spatial correlation of mode shapes, the coordinate modal assur- ance criterion (COMAC), in: Proceedings of the 6th International Modal Analysis Conference, pp. 690–695. 129 Lin, R.M., 2011. Identification of modal parameters of unmeasured modes using multiple FRF modal analysis method. Mechanical Systems and Signal Processing 25, 151–162. Lin, R.M., Ewins, D.J., 1994. Analytical model improvement using frequency response functions. Mechanical Systems and Signal Processing 8, 437–458. Lin, R.M., Zhu, J., 2006. Model updating of damped structures using FRF data. Mechanical Systems and Signal Processing 20, 2200–2218. Markovsky, I., Van Huffel, S., 2007. Overview of total least-squares methods. Signal processing 87, 2283–2302. McAllester, D.A., Schapire, R.E., 2000. On the convergence rate of Good-Turing estimators, in: Proceedings of the 13th Annual Conference on Computational Learning Theory, Palo Alto, CA, USA. pp. 1–6. Mei, G., Kareem, A., Kantor, J.C., 2001. Real-time model predictive control of structures under earthquakes. Earthquake Engineering & Structural Dynamics 30, 995–1019. Mei, G., Kareem, A., Kantor, J.C., 2002. Model predictive control of structures under earthquakes using acceleration feedback. Journal of Engineering Mechanics 128, 574–585. Mei, G., Kareem, A., Kantor, J.C., 2004. Model predictive control of wind-excited building: Benchmark study. Journal of Engineering Mechanics 130, 459–465. Mesarovic, V .Z., Galatsanos, N.P., Katsaggelos, A.K., 1995. Regularized constrained total least squares image restoration. IEEE Transactions on Image Processing 4, 1096–1108. Moaveni, B., Conte, J.P., Hemez, F.M., 2009. Uncertainty and sensitivity analysis of damage identification results obtained using finite element model updating. Computer-Aided Civil and Infrastructure Engineering 24, 320–334. Mokha, A.S., Constantinou, M.C., Reinhorn, A.M., 1993. Verification of friction model of Teflon bearings under triaxial load. ASCE Journal of Structural Engineering 119, 240–261. MATLAB, 2020. User’s Guide (R2020a). The MathWorks Inc., Natick, Massachusetts, USA. Nakashima, M., Nagae, T., Enokida, R., Kajiwara, K., 2018. Experiences, accomplishments, lessons, and challenges of E-defense—Tests using world’s largest shaking table. Japan Archi- tectural Review . Ng, M.K., Koo, J., Bose, N.K., 2002. Constrained total least-squares computations for high- resolution image reconstruction with multisensors. International Journal of Imaging Systems and Technology 12, 35–42. O’Callahan, J.C., 1995. Development of a general pseudo orthogonality correlation procedure, in: Proceedings of the 13th International Modal Analysis Conference, pp. 1013–1021. 130 Ogawa, N., Ohtani, K., Katayama, T., Shibata, H., 2001. Construction of a three-dimensional, large-scale shaking table and development of core technology. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 359, 1725–1751. Park, Y .J., Wen, Y .K., Ang, A.H.S., 1986. Random vibration of hysteretic systems under bi- directional ground motions. Earthquake Engineering and Structural Dynamics 14, 543–557. Passerini, A., Pontil, M., Frasconi, P., 2004. New results on error correcting output codes of kernel machines. IEEE Transactions on Neural Networks 15, 45–54. Peeters, B., 2000. System identification and damage detection in civil engeneering. Ph.D. thesis. Katholieke Universiteit te Leuven. Qin, J., Liu, L., Shao, L., Shen, F., Ni, B., Chen, J., Wang, Y ., 2017. Zero-shot action recognition with error-correcting output codes, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Honolulu, HI, USA. pp. 2833–2842. Qin, S.J., Badgwell, T.A., 2003. A survey of industrial model predictive control technology. Con- trol Engineering Practice 11, 733–764. Ramallo, J.C., Johnson, E.A., Spencer, Jr., B.F., 2002. “Smart” base isolation systems. Journal of Engineering Mechanics 128, 1088–1099. Sato, E., Sasaki, T., Fukuyama, K., Tahara, K., Kajiwara, K., 2013. Development of innovative base isolation system based on E-Defense full-scale shake table experiments part I: Outline of project research, in: AIJ Annual Meeting. Hokkaido, Japan, pp. 751–752. (in Japanese). Schaffrin, B., Wieser, A., 2008. On weighted total least-squares adjustment for linear regression. Journal of geodesy 82, 415–421. Scruggs, J.T., Taflanidis, A.A., Iwan, W.D., 2007. Non-linear stochastic controllers for semiac- tive and regenerative systems with guaranteed quadratic performance bounds-Part 2: Output feedback control. Structural Control and Health Monitoring 14, 1121–1137. Shadan, F., Khoshnoudian, F., Esfandiari, A., 2016. A frequency response-based structural damage identification using model updating method. Structural Control and Health Monitoring 23, 286– 302. Shi, Z., Law, S.S., Zhang, L.M., 1998. Structural damage localization from modal strain energy change. Journal of Sound and Vibration 218, 825–844. Silvia, M.T., Tacker, E.C., 1982. Regularization of Marchenko’s integral equation by total least squares. Journal of the Acoustical Society of America 72, 1202–1207. Song, W., Dyke, S.J., Harmon, T., 2012. Application of nonlinear model updating for a reinforced concrete shear wall. Journal of Engineering Mechanics 139, 635–649. 131 Song, W., Dyke, S.J., Harmon, T., So, M., 2009a. Nonlinear model updating in concrete struc- tures based on ambient response data, in: Proceedings of the 27th international modal analysis conference, In: Proceedings of the 27th international modal analysis conference, Orlando, USA. Song, W., Dyke, S.J., Yun, G., Harmon, T., 2009b. Improved damage localization and quantifica- tion using subset selection. Journal of Engineering Mechanics 135, 548–560. Sontag, E., 1981. Nonlinear regulation: The piecewise linear approach. IEEE Transactions on Automatic Control 26, 346–358. Soong, T.T., Grigoriu, M., 1993. Random vibrations of mechanical and structural systems. Prentice Hall, Englewood Cliffs, New Jersey, USA. Soong, T.T., Spencer, Jr., B.F., 2002. Supplemental energy dissipation: state-of-the-art and state- of-the-practice. Engineering Structures 24, 243–259. Spencer, Jr., B.F., Nagarajaiah, S., 2003. State of the art of structural control. Journal of Structural Engineering 129, 845–856. Stubbs, N., Kim, J.T., Farrar, C.R., 1995. Field verification of a nondestructive damage localization and severity estimation algorithm, in: Proceedings of the 13th International Modal Analysis Conference, 13–16 February 1995, Nashville, TN, Society for Experimental Mechanics, Inc., Bethel, CT. pp. 210–218. Symans, M.D., Charney, F.A., Whittaker, A.S., Constantinou, M.C., Kircher, C.A., Johnson, M.W., McNamara, R.J., 2008. Energy dissipation systems for seismic applications: current practice and recent developments. Journal of Structural Engineering 134, 3–21. Tsai, C.L., Kao, W.W., 2006. Constrained total least-square solution for GPS compass attitude determination. Applied Mathematics and Computation 183, 106–118. Ueda, T. (Ed.), 2007. Standard specifications for concrete structures – 2007 ‘Design’. JSCE Guidelines for Concrete, JSCE, No. 15. Vacher, P., Jacquier, B., Bucharles, A., 2010. Extensions of the MAC criterion to complex modes, in: Proceedings of the international conference on noise and vibration engineering, Leuven, Belgium. pp. 2713–2726. Van Huffel, S., Vandewalle, J., 1985. The use of total linear least squares techniques for identifi- cation and parameter estimation, in: IFAC Proceedings V olumes, pp. 1167–1172. Van Huffel, S., Vandewalle, J., 1989. Analysis and properties of the generalized total least squares problem AX≈ B when some or all columns in A are subject to error. SIAM Journal on Matrix Analysis and Applications 10, 294. Van Overschee, P., De Moor, B., 1994. N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica 30, 75–93. Vanik, M.W., Beck, J.L., Au, S., 2000. Bayesian probabilistic approach to structural health moni- toring. Journal of Engineering Mechanics 126, 738–745. 132 Wang, S., Li, Y ., Li, H., 2015. Structural model updating of an offshore platform using the cross model cross mode method: An experimental study. Ocean Engineering 97, 57–64. Wang, S., Xu, M., 2019. Modal strain energy-based structural damage identification: a review and comparative study. Structural Engineering International 29, 234–248. Wang, Y ., Boyd, S., 2010. Fast model predictive control using online optimization. IEEE Trans- actions on Control Systems Technology 18, 267–278. Wen, Y .K., 1980. Equivalent linearization for hysteretic systems under random excitation. ASME Journal of Applied Mechanics 47, 150–154. Wright, S.J., 1997. Applying new optimization algorithms to model predictive control. Chemical Process Control-V 93, 147–155. Wright, S.J., 2000. Interior-point methods. Journal of Computational and Applied Mathematics 124, 281–302. Yang, K., An, J., Bu, X., Sun, G., 2009. Constrained total least-squares location algorithm us- ing time-difference-of-arrival measurements. IEEE Transactions on Vehicular Technology 59, 1558–1562. Yang, S., Luo, P., Loy, C.C., Shum, K.W., Tang, X., 2015. Deep representation learning with target coding, in: Proceedings of the AAAI Conference on Artificial Intelligence, Austin, TX, USA. pp. 3848–3854. Yang, X., Ouyang, H., Guo, X., Cao, S., 2020. Modal strain energy-based model updating method for damage identification on beam-like structures. Journal of Structural Engineering 146, 04020246. Yildirim, E.A., Wright, S.J., 2002. Warm-start strategies in interior-point methods for linear pro- gramming. SIAM Journal on Optimization 12, 782–810. Yu, T., 2020. Novel model updating approaches and machine learning based semiactive model predictive control algorithms. Technical Report. Ph.D Dissertation Proposal 28 August 2020, University of Southern California. Ph.D Dissertation Proposal, University of Southern Califor- nia. Yu, T., Johnson, E.A., Brewick, P.T., Christenson, R.E., 2022. Modeling and model updating of a full-scale experimental base-isolated building. Engineering Structures , in review. Zhang, Q., Lee, K.C., Bao, H., You, Y ., Li, W., Guo, D., 2018. Large scale classification in deep neural network with label mapping, in: Proceedings of the IEEE International Conference on Data Mining Workshops, Singapore, Singapore. pp. 1134–1143. 133 Appendices A Two-Level Strategy for Optimizing Isolation-Layer Device Parameters To optimize the isolation-layer device parameters for the locally nonlinear FEM, a two-level hier- archical strategy is used. For each isolation-layer device D, the first level of the strategy selects a set of candidate values for a subset of the device parameters (or their bounds) and then, for each candidate subset, optimizes the other device parameters using a cost metricE D : E D (θ θ θ D )=E (Test010) D (θ θ θ D )+E (Test016) D (θ θ θ D ) (A.1a) E (T) D (θ θ θ D )= " ∑ nt− 1 k=0 [ ˆ f D x (k∆t;θ θ θ D )− f D x (k∆t)] 2 +[ ˆ f D y (k∆t;θ θ θ D )− f D y (k∆t)] 2 ∑ nt− 1 k=0 [f D x (k∆t)] 2 + f D y (k∆t) 2 # T (A.1b) whereθ θ θ D is a vector of the parameters of deviceD, and the device force predictions ˆ f D x (t;θ θ θ D ) and ˆ f D y (t;θ θ θ D ) are computed from the device model in (2.9) or (2.10)–(2.11) driven by the drifts and corresponding drift velocities across device D measured during test T . The steps of the first level are: 1a. RBs: For a candidate RBi preyield stiffness k RBi = k RBi xx = k RBi yy , taken from a set of n RBi (which is 7 herein) candidate values evenly spanning the range [1025,1175] kN/m (approxi- mately the LB and UB of RB stiffnesses in Table 2.7), determine the postyield stiffness ¯ k RBi and yield displacement ¯ u RBi that minimizeE RBi ([k RBi ¯ k RBi ¯ u RBi ] T ) using (A.1) for RBi. 134 1b. ESB1: For a pair of lower bounds on ESB1 damping coefficients c ESB1 xx and c ESB1 yy — each taken from a set of n ESB1 (which is 6 herein) bounds evenly spanning [0,50]Mg/s, thus forming a grid of n 2 ESB1 pairs — find the ESB1 parameter vector θ θ θ ESB1 = [c ESB1 xx c ESB1 yy β ESB1 γ ESB1 ] T that minimizeE ESB1 (θ θ θ ESB1 ). (Note: A ESBi ≡ β ESBi +γ ESBi .) 1c. ESB2: Identical to ESB1 except the damping coefficient lower bounds span [0,30]Mg/s (and n ESB2 is taken to be 4 herein). 1d. SDPs: For a candidate lower bound on c SDP xx and c SDP yy (the same lower bound for both di- rections and both SDPs) — taken from a set of n SDP (which is 11 herein) SDP damping coefficient lower bounds evenly spanning [0,100]Mg/s — determine the set of SDP pa- rameters θ θ θ SDP = [k SDP xx k SDP yy c SDP xx c SDP yy α SDP β SDP γ SDP ] T that minimizes the sum E SDP (θ θ θ SDP )=E SDP1 (θ θ θ SDP )+E SDP2 (θ θ θ SDP ). Lower bounds on the ESB and SDP damping coefficients are needed to ensure that the effective modal damping ratios identified from the predicted partially nonlinear FEM responses to the low- level random excitations rise to the level reported in Table 2.2; the upper end of the ranges for the sets of lower bounds are chosen to limit the increases in corresponding optimal device force cost metric — i.e., E ESB1 (θ θ θ ESB1 ), E ESB2 (θ θ θ ESB2 ) or E SDP (θ θ θ SDP ) — to no more than 15% above the optimal metric when the lower bound is zero. The second level of the optimization strategy then chooses from among the first level’s set of n combos =n RB1 · n RB2 · n 2 ESB1 · n 2 ESB2 · n SDP combinations of RB1, RB2, ESB1, ESB2 and SDP optimal parameters, which may be denotedθ θ θ 1 ,...,θ θ θ n combos . For each of then combos combinations, the partially nonlinear FEM time-history responses to low-level random excitation in Test 010 are computed; the relative errors in the predicted RMS of each of the twelve (six devices, each with two directions) isolation-layer device horizontal force time histories are computed. The second 135 level then chooses the combination of parameters that provides the smallest of the largest device force errors. In other words, choose the combination arg min k∈{1,2,...,n combos } max i∈{1,2} D∈{RB,ESB,SDP} max w∈{f Di x ,f Di y } Err RMS w (θ θ θ k ) (A.2) The resulting combination of optimal parameters is then listed in Table 2.11; this optimal combi- nation indeed results in relative errors in Test 010 RMS forces only marginally larger than those reported in Table 2.10 for the purely linear FEM (to force them the same would mean less faithful reproduction of the earthquake responses such as in Test 016). Note that, in some of the first level’s combinations, some optimal damping coefficients do converge to their corresponding lower bounds; indeed, the optimal parameters chosen in the second level is such a case, resulting in the damping coefficients for ESB2 and the SDPs to be at their corresponding discrete lower bound values, as shown in Table 2.11. B Data and Code Availability To facilitate future simulation studies on this full-scale base-isolated building, the following data and MATLAB ® codes will be made available in the DesignSafe-CI “Data Depot”: (a) table-center acceleration records computed from the shake table acceleration measurements during random excitation and earthquake ground motion tests; (b) system matrices M and K of the unupdated and updated superstructure linear FEMs; (c) codes for finding the node closest to any given 3D coordinates; (d) codes for time-history analyses using different isolation-layer device models; and (e) codes to simulate the response of the partially nonlinear FEM. 136 C Derivation for Eq. (3.16) This appendix shows the derivation for (3.16). Let f S (x)= f S (µ,χ,κ)= S ikk − S iku S i − 1 uu S iuk and linearize it w.r.t. x (orµ,χ andκ). Then expanding f S (x) as a first-order Maclaurin series f S (x)≈ f S | x=0 + ∑ h ∂ f S ∂x h x=0 x h = f S | x=0 + ∑ h ∂S ikk ∂x h x=0 − ∂S iku ∂x h S i − 1 uu S iuk x=0 + S iku S i − 1 uu ∂S iuu ∂x h S i − 1 uu S iuk x=0 − S iku S i − 1 uu ∂S iuk ∂x h x=0 x h (C.1) where f S | x=0 =λ 2 i ˆ M kk +λ i ˆ C kk + ˆ K kk − (λ 2 i ˆ M ku +λ i ˆ C ku + ˆ K ku )(λ 2 i ˆ M uu +λ i ˆ C uu + ˆ K uu ) − 1 (λ 2 i ˆ M uk +λ i ˆ C uk + ˆ K uk ) (C.2a) ∑ h ∂S ikk ∂x h x=0 x h = ∑ h µ h λ 2 i ˆ M hkk + ∑ h χ h λ i ˆ C hkk + ∑ h κ h ˆ K hkk (C.2b) ∑ h ∂S iku ∂x h S i − 1 uu S iuk x=0 x h = ∑ h µ h λ 2 i ˆ M hku + ∑ h χ h λ i ˆ C hku + ∑ h κ h ˆ K hku ! (λ 2 i ˆ M uu +λ i ˆ C uu + ˆ K uu ) − 1 (λ 2 i ˆ M uk +λ i ˆ C uk + ˆ K uk ) (C.2c) ∑ h S iku S i − 1 uu ∂S iuu ∂x h S i − 1 uu S iuk x=0 x h =(λ 2 i ˆ M ku +λ i ˆ C ku + ˆ K ku )(λ 2 i ˆ M uu +λ i ˆ C uu + ˆ K uu ) − 1 ∑ h µ h λ 2 i ˆ M huu + ∑ h χ h λ i ˆ C huu + ∑ h κ h ˆ K huu ! (λ 2 i ˆ M uu +λ i ˆ C uu + ˆ K uu ) − 1 (λ 2 i ˆ M uk +λ i ˆ C uk + ˆ K uk ) (C.2d) ∑ h S iku S i − 1 uu ∂S iuk ∂x h x=0 x h =(λ 2 i ˆ M ku +λ i ˆ C ku + ˆ K ku )(λ 2 i ˆ M uu +λ i ˆ C uu + ˆ K uu ) − 1 ∑ h µ h λ 2 i ˆ M hku + ∑ h χ h λ i ˆ C hku + ∑ h κ h ˆ K hku ! (C.2e) 137 Note that matrix ˆ S i was defined as ˆ S i = λ 2 i ˆ M + λ i ˆ C + ˆ K and partitioned as ˆ S i = [[ ˆ S i T uu ˆ S i T ku ] T [ ˆ S i T uk ˆ S i T kk ] T ] T , then f S (µ,χ,κ) can be simplified to f S (µ,χ,κ)≈ ˆ S ikk − ˆ S iku ˆ S i − 1 uu ˆ S iuk + ∑ h µ h λ 2 i ˆ M h + ∑ h χ h λ i ˆ C h + ∑ h κ h ˆ K h ! kk − ∑ h µ h λ 2 i ˆ M h + ∑ h χ h λ i ˆ C h + ∑ h κ h ˆ K h ! ku ˆ S i − 1 uu ˆ S iuk + ˆ S iku ˆ S i − 1 uu ∑ h µ h λ 2 i ˆ M h + ∑ h χ h λ i ˆ C h + ∑ h κ h ˆ K h ! uu ˆ S i − 1 uu ˆ S iuk − ˆ S iku ˆ S i − 1 uu ∑ h µ h λ 2 i ˆ M h + ∑ h χ h λ i ˆ C h + ∑ h κ h ˆ K h ! uk = ˆ S ikk − ˆ S iku ˆ S i − 1 uu ˆ S iuk + h ˆ S iku ˆ S i − 1 uu − I i ∑ h µ h λ 2 i ˆ M h + ∑ h χ h λ i ˆ C h + ∑ h κ h ˆ K h !    ˆ S i − 1 uu ˆ S iuk − I    (C.3) which is linear w.r.t. x (orµ,χ andκ). This concludes the derivation. D Derivation for Eq. (3.22) This appendix shows the derivation for (3.22). According to (3.17)–(3.18) and noting thatφ φ φ k i = P i w, the((i− 1)n FEM +j) th row of matrix A M h is h A M h i (i− 1)n FEM +j = ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    φ φ φ k i . (D.1) Thus, the((i− 1)n FEM + j) th row of changes∆A M h to matrix A M h due to∆λ i and∆φ φ φ k i is h ∆A M h i (i− 1)n FEM +j = ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i ({λ i +∆λ i } 2 ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    (φ φ φ k i +∆φ φ φ k i )− h A M h i (i− 1)n FEM +j . (D.2) 138 Ignoring terms with∆λ r i ∆φ φ φ k i s andr+s> 1 leads to h ∆A M h i (i− 1)n FEM +j = ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    ∆φ φ φ k i + ∂ ∂λ i      ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    φ φ φ k i      ∆λ i . (D.3) Note that ∂ ∂λ i      ˆ φ φ φ k j H h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    φ φ φ k i      = ˆ φ φ φ k j H      ∂ h ˆ S iku ˆ S i − 1 uu − I i ∂λ i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    + h ˆ S iku ˆ S i − 1 uu − I i ∂λ 2 i ˆ M h ∂λ i    ˆ S i − 1 uu ˆ S iuk − I    + h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h ) ∂ ∂λ i    ˆ S i − 1 uu ˆ S iuk − I         φ φ φ k i , (D.4) and that ˆ S i − 1 uu ∂λ i =− ˆ S i − 1 uu ∂ ˆ S iuu ∂λ i ˆ S i − 1 uu , ∂ ˆ S iuu ∂λ i = ˆ S ′ iuu , ∂ ˆ S iuk ∂λ i = ˆ S ′ iuk , ∂ ˆ S iku ∂λ i = ˆ S ′ iku , (D.5) 139 then h ∆A M h i (i− 1)n FEM +j = h F M h i (i− 1)n FEM +j ∆w+ ˆ φ φ φ k j H      h ˆ S ′ iku ˆ S i − 1 uu − ˆ S iku ˆ S i − 1 uu ˆ S i ′ uu ˆ S i − 1 uu 0 i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    + h ˆ S iku ˆ S i − 1 uu − I i (2λ i ˆ M h )    ˆ S i − 1 uu ˆ S iuk − I    + h ˆ S iku ˆ S i − 1 uu − I i (λ 2 i ˆ M h )    ˆ S i − 1 uu ˆ S ′ iuk − ˆ S i − 1 uu ˆ S i ′ uu ˆ S i − 1 uu ˆ S iuk 0         φ φ φ k i q T i ∆w (D.6) This concludes the derivation. E Strategy Switch When Relative Velocity Switch Sign This Appendix proves that, when the state bounds (4.6g) are inactive, if ξ ξ ξ ∗ = h X ∗ T U ∗ T ∆ ∆ ∆ ∗ T i T minimizes (4.6) given initial state vector x 0 with strategyS I ≡ ∆ ∆ ∆ ∗ or strategyS II ≡ [∆ ∆ ∆ ∗ T θ θ θ T ] T , then ˆ ξ ξ ξ ∗ = h − X ∗ T − U ∗ T ˆ ∆ ∆ ∆ ∗ T i T minimizes (4.6) given initial state vector− x 0 for strategy ˆ S I ≡ ˆ ∆ ∆ ∆ ∗ = 1 m(p− 1)× 1 − ∆ ∆ ∆ ∗ or strategy ˆ S II ≡ [ ˆ ∆ ∆ ∆ ∗ T ˆ θ θ θ T ] T =[ ˆ ∆ ∆ ∆ ∗ T − θ θ θ T ] T . To prove this, consider each element of (4.6): (4.6a) J is quadratic in the vector h x T 0 X ∗ T U ∗ T i T , soJ is invariant when changing the signs of the states and control forces. (4.6b) If the state equation holds for h x T 0 X ∗ T U ∗ T i T , then it also holds for h − x T 0 − X ∗ T − U ∗ T i T , because multiplying both the left- and right-hand-side by a constant coefficient − 1 does not change the equality. (4.6c) If δ ∗ k,j and v rel k,j ∗ (determined from X ∗ ) satisfy (4.6d), and v min =− v max , then 1− δ ∗ k,j and − v rel k,j ∗ (determined from− X ∗ ) satisfy (4.6c). 140 (4.6d) Similarly, ifδ ∗ k,j and v rel k,j ∗ satisfy (4.6c), and v min =− v max , then 1− δ ∗ k,j and− v rel k,j ∗ satisfy (4.6d). Further, if (4.6c) is active for δ ∗ k,j and v rel k,j , then (4.6d) is active for 1− δ ∗ k,j and− v rel k,j ∗ , and vice versa — this switch of active constraints is reflected by ˆ θ θ θ v =− θ θ θ v . (4.6e) Same as (4.6c) but replace v,v and X with u,u, and U respectively. (4.6f) Same as (4.6d) but replace v,v and X with u,u, and U respectively. The result is that ˆ θ θ θ u =− θ θ θ u . (4.6g) If x k satisfies the state bounds, and x min =− x max , then− x k also satisfies the bounds; al- ternately, if both state bounds are chosen large enough that (4.6g) is always inactive, then (4.6g) may be ignored regardless of the states’ signs. Thus, the proof is concluded. F Matrices and Vectors in Eq. (4.10) Using the notation that blkdiag(Y 1 ,Y 2 ,...) denotes a block diagonal matrix of size(∑ j n rows Y j )× (∑ j n cols Y j ) where matrix Y j is n rows Y j × n cols Y j , and diag(y) is a square matrix with diagonal elements taken from the vector y, the matrices and vectors in (4.10) are: H=blkdiag(Q ′ ,...,Q ′ | {z } p− 1 times ,Q p ,R,R,...,R | {z } p times ,0 m(p− 1) ) G=[G 1 G 2 0 2np× m(p− 1) ] G 1 = I 2np −    0 2n× 2n(p− 1) 0 2n blkdiag( ¯ A ′ ,..., ¯ A ′ | {z } p− 1 times ) 0 2n(p− 1)× 2n    141 G 2 =− blkdiag(B u ,B u ,...,B u | {z } p times ) b=[x T 0 ¯ A ′ T 0 1× 2n(p− 1) ] T ¯ F=                 − ¯ F 1 0 m(p− 1)× mp − ¯ F 2 ¯ F 1 0 m(p− 1)× mp − ¯ F 3 ¯ F 4 − I mp − ¯ F 5 − ¯ F 4 I mp − ¯ F 6 I 2np 0 2np× mp 0 2np× m(p− 1) − I 2np 0 2np× mp 0 2np× m(p− 1)                 ¯ F 1 =blkdiag(V,...,V | {z } p− 1 times ) ¯ F 2 =blkdiag(diag(v min − ϵ ),...,diag(v min − ϵ ) | {z } p− 1 times ) ¯ F 3 =blkdiag(diag(v max +ϵ ),...,diag(v max +ϵ ) | {z } p− 1 times ) ¯ F 4 =    0 m× 2n(p− 1) 0 m× 2n , blkdiag(R − 1 N T ,...,R − 1 N T | {z } p− 1 times ) 0 m(p− 1)× 2n    ¯ F 5 = 0 m(p− 1)× m blkdiag(diag(u min − ϵ ),...,diag(u min − ϵ ) | {z } p− 1 times ) T T ¯ F 6 = 0 m(p− 1)× m blkdiag(diag(u max +ϵ ),...,diag(u max +ϵ ) | {z } p− 1 times ) T T ¯ h=[ ¯ h T 1 ¯ h T 2 ¯ h T 3 ¯ h T 4 ¯ h T 5 ¯ h T 6 ] T ¯ h 1 =− [v min T ... v min T | {z } p− 1 times ] T ¯ h 2 =− ϵ 1 m(p− 1)× 1 ¯ h 3 =− h {R − 1 N T x 0 − (u min − ϵ )δ δ δ 0 + u min } T u min T ... u min T | {z } p− 1 times i T 142 ¯ h 4 = {R − 1 N T x 0 +(u max +ϵ )δ δ δ 0 − ϵ } T − ϵ 1 1× m(p− 1) T ¯ h 5 =+[x maxT x maxT ... x maxT | {z } p times ] T ¯ h 6 =− [x min T x min T ... x min T | {z } p times ] T 143

Abstract (if available)

Abstract This dissertation seeks to further knowledge of two important topics in structural dynamics: obtaining high-fidelity numerical models for physical structures by updating nominal finite element models (FEM) and developing control algorithms to mitigate the structural vibration response to natural hazards. ❧ The first contribution of this dissertation is to develop a multi-stage model updating method for base-isolated structures together with a mode-matching algorithm. Given the observation that the superstructure predominantly moves as a rigid body in low-frequency modes and the isolation layer plays a minor role in all other modes, this study proposes updating parameters in stages (instead of altogether): first, the linear superstructure parameters are updated so that its natural frequencies and mode shapes match those identified via a subspace system identification of the experimental building responses to low-level random excitations; then, the isolation-layer device linear parameters are updated so that the natural frequencies, damping ratios and mode shapes of the three isolation modes match. This two-step process breaks a large-scale linear model updating problem into two smaller problems, thereby reducing the search space for the to-be-updated parameters, which generally reduces computational costs regardless of what optimization algorithm is adopted. Due to the limited instrumentation, the identified modes constitute only a subset of all the modes; to match each identified mode with a finite element model mode, a procedure is proposed to compare each identified mode with a candidate set of FEM modes and to select the best match. Finally, nonlinear models of the isolation-layer devices are calibrated so that the model well-predicts the physical specimen’s response to earthquake excitation. ❧ The second contribution is to introduce a novel solution algorithm into the cross-model cross-mode (CMCM) model updating method, which is a physics-based method that avoids mode matching. Due to measurement noise in, and the spatial incompleteness of, the experimental mode shapes, this method involves solving a set of ill-conditional linear equations. This study shows that the constrained total least squares (CTLS) algorithm can provide more accurate solutions, given the same level of noise and spatial incompleteness, compared with the current state-of-the-art total least squares (TLS). To demonstrate this observation, both methods are first used to solve a simple line-fitting problem with noisy measurement data, with pros and cons of the methods highlighted. Then, they are applied to the model updating problem of a three-dimensional frame structure using a structure model that closely resembles one used in previous CMCM studies. Results indicate that the CTLS method is a better alternative to the TLS method when dealing with the set of ill-conditional linear equations in the CMCM method. ❧ Both model updating methods are validated with the experimental data of a full-scale four-story base-isolated reinforced-concrete frame building that was tested in 2013 at the NIED “E-Defense” laboratory in Japan, proving their applicability to FEMs of full-scale structures. ❧ The final contribution is to develop a real-time neural network based semiactive model predictive control (MPC) strategy for reducing structural vibrations. Semiactive MPC algorithms perform well over a range of scenarios; however, their computational costs due to the inherent mixed integer quadratic programming (MIQP) optimization problem prevent them from being adopted in real-time. This dissertation proposes that neural network (NN) classifiers be trained to predict in real-time the integer variable values (each unique combination of these is called a strategy) for a given current state variable vector, reducing MPC control force optimization to a quadratic programming (QP) problem, which is solved using a fast QP solver adapted to this type of problem. The number of possible strategies is typically huge, as it grows exponentially with the number of integer variables (determined by the numbers of prediction horizon steps and damping devices); accurately matching the current state vector (the input to the MIQP) with one of these strategies typically requires an infeasibly massive NN. To reduce the number of possible strategies, the control algorithm proposed herein first exploits the homogeneity-of-order-one nature of this structural control problem and the statistics of the system states to shrink the input parameter space, leading to a smaller set of corresponding output strategies; additionally, it is proposed to use a few smaller trained NNs running simultaneously, each linking an input to a unique group of strategies, instead of one large NN that takes much longer to run, to jointly predict the correct strategy for a given state vector input to the NN. Further, occasional inaccurate NN predictions can lead to sub-optimal solutions; to enforce a bound on the sub-optimality, the control algorithm will switch to a simpler clipped linear quadratic regulator (CLQR) algorithm if the CLQR cost function is smaller than that calculated based on the predicted strategy. Numerical examples using shear buildings of one, three and five degrees-of-freedom (DOFs) demonstrate that the proposed real-time algorithm speeds up the online computation by approximately one order of magnitude and offers a control performance close to the state-of-the-art offline semiactive MPC algorithm.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Structural nonlinear control strategies to provide life safety and serviceability

PDF

Optimal clipped linear strategies for controllable damping

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Performance monitoring and disturbance adaptation for model predictive control

PDF

Provable reinforcement learning for constrained and multi-agent control systems

PDF

Economic model predictive control for building energy systems

PDF

Non‐steady state Kalman filter for subspace identification and predictive control

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Robustness of gradient methods for data-driven decision making

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Efficient control optimization in subsurface flow systems with machine learning surrogate models

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

Shrinkage methods for big and complex data analysis

PDF

Iterative path integral stochastic optimal control: theory and applications to motor control

PDF

Algorithms and frameworks for generating neural network models addressing energy-efficiency, robustness, and privacy

PDF

Closing the reality gap via simulation-based inference and control

PDF

Learning and control for wireless networks via graph signal processing

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

Asset Metadata

Creator Yu, Tianhao (author)

Core Title Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Degree Conferral Date 2022-05

Publication Date 02/14/2022

Defense Date 01/20/2022

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag base-isolated building,constrained total least squares,cross-model cross-mode method,frame structure,model predictive control,model updating,neural network,OAI-PMH Harvest,real-time control,semiactive control,total least squares

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Johnson, Erik A. (committee chair), Jovanovic, Mihailo R. (committee member), Masri, Sami F. (committee member)

Creator Email chasing1990yth@outlook.com,yutianha@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC110702284

Unique identifier UC110702284

Legacy Identifier etd-YuTianhao-10390

Document Type Dissertation

Format application/pdf (imt)

Rights Yu, Tianhao

Type texts

Source 20220214-usctheses-batch-912 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu