Adjoint-based optimal design of adaptive negative stiffness devices for inelastic frame structures

Abstract

This paper presents a time-domain optimal design procedure for adaptive negative stiffness devices (ANSDs) installed in inelastic building frames. The structure is modeled as a nonlinear shear-type system with hysteretic restoring forces, classical viscous damping with approximately uniform modal damping ratios, and story-wise ANSDs composed of a negative stiffness device and a nonlinear fluid viscous damper. A state-space formulation suitable for gradient-based optimization is derived, and an objective functional is defined that combines a shear-energy measure based on squared story-shear histories over an ensemble of ground motions with penalties on changes in ANSD stiffness and damping coefficients. Sensitivities of this functional with respect to the device parameters are obtained via an adjoint (Lagrange multiplier) formulation, evaluated in discrete time, which enables an efficient iterative tuning scheme. The methodology is demonstrated using an eight-story reinforced-concrete frame, idealized as a planar shear building, and subjected to 15 recorded earthquakes. The tuned devices reduce the shear-energy measure for most records, yield modest decreases in peak base shear and floor accelerations for critical motions, and prevent collapse under one previously critical record with only moderate adjustments of the ANSD parameters.

Keywords

seismic vibration control adaptive negative stiffness devices inelastic frame structures adjoint-based optimal design nonlinear fluid viscous dampers

Introduction

Modern seismic design must address conflicting objectives: ensuring structural safety, limiting damage and economic loss, and maintaining serviceability under strong ground motions. Achieving these goals is challenging when structures are driven into the nonlinear or inelastic range. Optimal control provides a framework in which seismic performance targets—such as reduced displacements, accelerations, or energy dissipation—are encoded in a cost functional subject to the structural dynamics, guiding the tuning of passive and semi-active devices in a systematic way rather than by trial and error (Takewaki and Akehashi, 2021; Zakian and Kaveh, 2023).

The integration of optimal and robust control into structural engineering has advanced the design of passive and semi-active systems that strategically adjust stiffness and damping to enhance seismic performance. Early work by Kelly et al. (1987), Gluck et al. (1996), Takewaki (1997), and Takewaki and Yoshitomi (1998) established key concepts for damper placement and parameter tuning using transfer-function norms and LQR-based formulations. Subsequent studies extended these ideas to more complex devices and scenarios, for example, analysis-redesign procedures for drift-constrained frames (Shmerling et al., 2018), robust semi-active base isolation (Zhang et al., 2020), short-horizon acceleration-predictive control with fluid viscous dampers (Shmerling and Gerdts, 2023), and frequency-response-based optimization of viscous dampers (Gao and Lu, 2024). In this context, Wang et al. (2021) proposed a fixed-point optimum design and performance evaluation of a tuned inerter–negative-stiffness damper for seismic protection of single-degree-of-freedom structures, and Gao et al. (2023) extended this line of work to performance improvement and demand-oriented optimum design of tuned negative-stiffness inerter dampers for base-isolated structures. Many of these approaches, however, rely on linear-elastic or equivalent-linear representations of earthquake response (Cimellaro, 2009; Etedali and Tavakoli, 2017; Tehrani and Harvey, 2019; Özuygur and Gündüz, 2018).

Nonlinearity and interaction effects have been considered in several contributions, including soil–structure interaction (Özuygur and Gündüz, 2018), dual-mode systems under extreme loading (Tehrani and Harvey, 2019), magnetorheological devices under nonlinear excitation (Pohoryles and Duffour, 2015), and broader reviews of optimal-control techniques for building frames (Ho, 2024). Nevertheless, linear or equivalent-linear models cannot fully capture the degradation of stiffness, pinching, or residual deformations in reinforced concrete and steel frames under strong shaking (Oh et al., 2023). This motivates methods that can operate directly on inelastic time-domain response while remaining computationally tractable.

In parallel, negative stiffness devices (NSDs) have attracted growing interest as tools for vibration and seismic control. Reviews by Li et al. (2020) and Wang et al. (2019) summarize how negative stiffness can be exploited to enhance isolation and damping. Analytical and experimental studies have demonstrated the potential of NSDs for seismic protection, including analytical and shake-table investigations of negative-stiffness devices for isolated structures (Sarlis et al., 2011, 2016), nonlinear viscous-damping-based negative stiffness isolation systems (Hu et al., 2024), and H₂-optimized negative-stiffness and inerter-based dampers (Islam and Jangid, 2023). From a performance-based perspective, Idels and Lavan (2021) proposed an optimization-based seismic retrofitting methodology for frame structures using negative stiffness devices and fluid viscous dampers, explicitly targeting drift and acceleration criteria. These works demonstrate that negative stiffness mechanisms can be powerful components in seismic control strategies, while also highlighting that device performance is sensitive to modeling assumptions and the choice of objective functions.

Within this context, there remains a need for methodologies that operate directly on inelastic time-domain response, use objective functionals with clear mechanical interpretation, and are compatible with gradient-based tools for realistic frame models and ensembles of ground motions. Many control-based formulations require practitioners to assign abstract weighting coefficients that balance state and control terms, often tuned heuristically, which obscures the link between optimization results and structural performance. The present study addresses these issues by focusing on the optimal tuning of adaptive negative stiffness devices (ANSDs) installed in an inelastic frame model. The frame is idealized as a planar shear building with hysteretic restoring forces and classical viscous damping, characterized by approximately uniform modal damping ratios. Story-wise ANSDs, each composed of a negative stiffness device and a nonlinear fluid viscous damper, are incorporated into the system. A state-space formulation suitable for adjoint-based sensitivity analysis is developed, and a control-inspired objective functional is introduced that aggregates squared story-shear trajectories over a set of ground motions, together with regularization terms on ANSD stiffness and damping coefficients. This gives the weighting matrices a direct mechanical interpretation in terms of shear-force energy and device effort, and the resulting adjoint equations yield gradients that enable an iterative tuning procedure based on the most critical records in the ensemble.

Inelastic system state-space formulation

The inelastic state-space formulation is necessary for control applications. Its formulation is derived in the following section. We begin with the differential equation that governs the inelastic dynamic response, specifically the equilibrium between the horizontal forces applied at the ceiling levels of the structure. The equilibrium is given by

f^{R} (δ, \dot{δ}) + f^{inh} (\dot{δ}) + f^{I} (\ddot{δ}) = p (t) + u (t)

(1)

So that for a structure with N-stories (ceiling levels), $f^{R}$ is the N-dimensional resisting force vector, $f^{inh}$ is the N-dimensional inherent damping force vector, $f^{I}$ is the N-dimensional inertia force vector, $u$ is the N-dimensional control force vector applied by the ANSD design variables, $p$ is the N-dimensional vector representing the applied dynamic load (i.e., earthquake load), and $δ$ is the N-dimensional interstory drift vector whose components are depicted in Figure 1. For the ideal state-space formulation, the inelastic resisting force $f^{R}$ is divided into its hysteretic portion (relating to the elasto-plastic and plastic behavior) and the linear-elastic portion, as depicted in Figure 1 This separation is expressed as

f^{R} (δ, \dot{δ}) = f^{H} (δ, \dot{δ}) + k^{E} δ (t)

(2)

where

f^{H}

is the N-dimensional hysteretic force vector and

k^{E}

is the 2 N × 2 N linear-elastic stiffness matrix. Furthermore, the force vectors

f^{inh}

and

f^{I}

are defined as

f^{inh} (\dot{δ}) = c^{inh} \dot{δ} (t)

(3)

f^{I} (\ddot{δ}) = m \ddot{δ} (t)

(4)

So that

m

is the diagonal N × N mass matrix (composed of the ceiling masses), and

c^{inh}

is the N × N inherent viscous damping (defined according to Caughey classical damping).

Figure 1.

Frame structure under horizontal deflections.

When referting to the horizontal force equilibrium about the drift coordinates (i.e., referring to $δ$ and its derivatives), the mass matrix is of the form:

m = T_{f \to v} [\begin{array}{l} m_{1} \\ ⋱ \\ m_{N} \end{array}] T_{f \to v}^{'}

(5)

So that $m_{n}$ is the nth story mass, and $T_{f \to v}$ is the transformation matrix from applied lateral forces into applied shear forces, which is given by

T_{f \to v} = [\begin{array}{l} 1 & \dots & 1 \\ ⋱ & ⋮ \\ 0 & 1 \end{array}]

(6)

Also, in drift coordinates, the applied horizontal earthquake load is expressed as

p (t) = - T_{f \to v} m 1 a_{g} (t)

(7)

where

1

is an N-dimensional vector of ones, and

a_{g} (t)

is the ground acceleration function.

The inelastic state-space formulation takes after the closed-loop paradigm of Figure 2, which enhances the nonlinear control model suggested by Shmerling and Gerdts (2022) and represents equation (1) as

\dot{z} (t) = A z (t) + B (- f^{H} (z) + u (t) + p (t))

(8)

where

A

is the 2 N × 2 N state matrix,

z

is the 2 N-dimensional state vector,

B

is the 2 N × N state matrix,

u

is the N-dimensional ANSD control force vector,

p (t)

is the N-dimensional effective earthquake load vector, and

\dot{z} (t)

is the 2N-dimensional state vector.

Figure 2.

Inelastic system state-space paradigm.

The state vector $z$ is composed of the interstory drift deflection and its derivatives (i.e., drift velocities):

z (t) = {[\begin{array}{l} δ (t) & \dot{δ} (t) \end{array}]}^{'}

(9a)

\dot{z} (t) = {[\begin{array}{l} \dot{δ} (t) & \ddot{δ} (t) \end{array}]}^{'}

(9b)

The state and input-to-state matrices are defined as

A = [\begin{array}{c} 0 & I \\ - m^{- 1} k^{E} & - m^{- 1} c^{inh} \end{array}]

(10)

B = [\begin{array}{c} 0 \\ - m^{- 1} \end{array}]

(11)

So that,

0

the N × N zero matrix,

I

the N × N identity matrix.

The control force generated by ANSD is composed of two components: the force from the NSD and the force from the NFVD, expressed as

u (t) = - f^{NSD} (t) - f^{NFVD} (t)

(12)

where

\begin{array}{l} f^{NSD} (t) = [\begin{array}{c} f_{1}^{NSD} (t) \\ ⋮ \\ f_{N}^{NSD} (t) \end{array}] & ; & f^{FVD} (t) = [\begin{array}{c} f_{1}^{NFVD} (t) \\ ⋮ \\ f_{N}^{NFVD} (t) \end{array}] \end{array}

(13)

To support the optimal control formulation developed later in this study, the traditional NFVD force model is reformulated. The classical model defines the damper force at each story level n as

\begin{array}{l} f_{n}^{NFVD} (t) = Δ c_{n} {| {\dot{δ}}_{n} (t) |}^{α} sgn ({\dot{δ}}_{n}) & \leftrightarrow & 0 < α \leq 1.0 & \forall & n = 1, \dots N \end{array}

(14)

where coefficients

Δ c_{n}

are the n^th story damping coefficients, and

α

is the parameter controlling the nonlinearity of the NFVD force. When α = 1.0, the NFVD force is linear to the velocity through

Δ c_{n}

. When α < 1.0, nonlinearity is introduced. The closer α is to 0, the more the NFVD force behavior becomes similar to that of a friction force. Equation (14) formulation, however, presents challenges for gradient-based optimization due to its non-smooth behavior near zero velocity. To address this, the NFVD force is equivalently expressed using a formulation that leverages the Dirac delta function:

\begin{array}{l} f_{n}^{NFVD} (t) = Δ c_{n} {| {\dot{δ}}_{n} (t) - D ({\dot{δ}}_{n}) |}^{α - 1} {\dot{δ}}_{n} (t) & \forall & n = 1, \dots N \end{array}

(15)

where

D ({\dot{δ}}_{n})

is a smoothing function designed to regularize the nonlinearity near zero velocity. It should be emphasized that the term

D ({\dot{δ}}_{n})

is introduced solely as a numerical regularization of the NFVD power-law near zero velocity. In particular, it does not represent an additional physical force, but serves only to avoid evaluating expressions of the form

{| {\dot{δ}}_{n} (t) |}^{α - 1}

{\dot{δ}}_{n} (t) = 0

when

α < 1

, for which the limit is undefined. For

| {\dot{δ}}_{n} (t) |

away from zero, the regularizing term vanishes, and the damper force reduces exactly to the standard power-law formulation, so that the force–velocity relation and its derivative are modified only in an arbitrarily small neighborhood of

{\dot{δ}}_{n} (t) = 0

. In this way, the smoothing has no role beyond ensuring a well-posed and differentiable NFVD model for time integration and adjoint-based gradient computations, and it does not affect the damper behavior in the velocity range where most seismic energy is dissipated. Additionally, when expressed in terms of the state-space representation, the NFVD force vector can be written compactly as

f^{NFVD} (t) = d i a g {Δ c} [[\begin{array}{l} 0 & I \end{array}] z (t) ⊙ {| [\begin{array}{l} 0 & I \end{array}] z (t) - D ([\begin{array}{l} 0 & I \end{array}] z (t)) |}^{α - 1}]

(16)

where

Δ c = {[\begin{array}{l} Δ c_{1} & \dots & Δ c_{N} \end{array}]}^{'}

, and

⊙

denotes the Hadamard (element-wise) product. This compact form is particularly useful for embedding the NFVD model into the overall control system formulation and supports efficient numerical implementation.

Figure 3 illustrates the geometric configuration of the Negative Stiffness Device (NSD), while Figure 4 presents its mechanical behavior under lateral deflection. The control force generated by the NSD at the n^th story is governed by the following nonlinear force law:

\begin{array}{l} f_{n}^{NSD} (t) = - [\frac{P_{n}^{in} + Δ k_{n} l_{n}^{p}}{l_{n}^{s} (t)} - Δ k_{n}] δ_{n} (t) \frac{l_{n}^{α}}{l_{n}^{β}} [2 + \frac{l_{n}^{β}}{l_{n}^{α}} + \frac{l_{n}^{p} + l_{n}^{α}}{\sqrt{| {(l_{n}^{β})}^{2} - {(δ_{n} (t))}^{2} |}}] + F_{n}^{g} (t) & \forall & n = 1, \dots N \end{array}

(17)

Here,

P_{n}^{i n}

is the pre-compressed force applied to the vertical spring,

Δ k_{n}

is the stiffness of the vertical spring,

l_{n}^{p}

is the initial vertical spring length,

l_{n}^{α}

and

l_{n}^{β}

denote the lengths of the pivot arms below and above the rotation point, respectively,

l_{n}^{s}

is the vertical spring length, and

F_{n}^{g}

is the additional spring force introduced by the NSD bottom spring assembly. The vertical spring length

l_{n}^{s}

is a nonlinear function of the lateral deflection, defined by

\begin{array}{l} l_{n}^{s} (t) = \sqrt{{[l_{n}^{p} + l_{n}^{α} - l_{n}^{α} \sqrt{| 1 - {(\frac{δ_{n} (t)}{l_{n}^{β}})}^{2} |}]}^{2} + {(δ_{n} (t))}^{2} {(1 + \frac{l_{n}^{α}}{l_{n}^{β}})}^{2}} & \forall & n = 1, \dots N \end{array}

(18)

Figure 3.

Negative stiffness damper geometrical parameters.

Figure 4.

Negative stiffness damper under lateral deflection.

The bottom spring force $F_{n}^{g} (t)$ is modeled as a two-stage linear spring system:

\begin{array}{l} F_{n}^{g} (t) = {\begin{cases} k_{n}^{s 1} δ_{n} (t) & , & 0 \leq | δ_{n} (t) | \leq δ_{n}^{y^{'}} \\ k_{n}^{s 1} δ_{n}^{y^{'}} + \frac{k_{n}^{s 2} k_{n}^{s 1}}{k_{n}^{s 2} + k_{n}^{s 1}} (δ_{n} (t) - δ_{n}^{y^{'}}) & , & | δ_{n} (t) | \geq δ_{n}^{y^{'}} \end{cases} & \leftrightarrow & δ_{n}^{y^{'}} = \frac{p_{n}^{is 2}}{k_{n}^{s 1}} & \forall & n = 1, \dots N \end{array}

(19)

where

k_{n}^{s 1}

and

k_{n}^{s 2}

are the stiffnesses of the first and second sub-springs, respectively,

p_{n}^{i s 2}

is the preload force applied to the second sub-spring, and

δ_{n}^{y^{'}}

is the deflection threshold at which the spring behavior transitions from the first to the second stage.

Equation (17) reveals that the NSD force comprises two components: a primary term governed by the interstory drift $δ_{n}$ , and a secondary term, $F_{n}^{g}$ , induced by the bottom springs. In this study, the design process emphasizes optimization of the first component by calibrating the vertical spring stiffness coefficients $Δ k_{n} \forall n = 1, \dots N$ . To facilitate this process, the following relationship is introduced:

\begin{array}{l} P_{n}^{in} = Δ k_{n} Δ_{n}^{p} & \forall & n = 1, \dots N \end{array}

(20)

where

Δ_{n}^{p}

represents the target pre-compression length associated with the vertical spring of the NSD at story n. Substituting equation (20) into equation (17), we obtain the reformulated NSD force law:

\begin{array}{l} f_{n}^{NSD} (t) = - Δ k_{n} [\frac{Δ_{n}^{p} + l_{n}^{p}}{l_{n}^{s} (t)} - 1] δ_{n} (t) \frac{l_{n}^{α}}{l_{n}^{β}} [2 + \frac{l_{n}^{β}}{l_{n}^{α}} + \frac{l_{n}^{p} + l_{n}^{α}}{\sqrt{| {(l_{n}^{β})}^{2} - {(δ_{n} (t))}^{2} |}}] + F_{n}^{g} (t) & \forall & n = 1, \dots N \end{array}

(21)

In state-space vector form, with respect to the state vector

z (t)

, the complete NSD force vector can be compactly written as

f^{NSD} (t) = - diag {Δ k} [\frac{Δ^{p} + l^{p}}{l^{s} (t)} - 1] ⊙ [\begin{array}{l} I & 0 \end{array}] z (t) ⊙ \frac{l^{α}}{l^{β}} ⊙ [2 + \frac{l^{β}}{l^{α}} + \frac{l^{p} + l^{α}}{\sqrt{| {(l^{β})}^{2} - {([\begin{array}{l} I & 0 \end{array}] z (t))}^{2} |}}] + F^{g} (t)

(22)

Where all vector terms are defined component-wise as

\begin{array}{l} l^{α} = [\begin{array}{c} l_{1}^{α} \\ ⋮ \\ l_{N}^{α} \end{array}]; & l^{β} = [\begin{array}{c} l_{1}^{β} \\ ⋮ \\ l_{N}^{β} \end{array}]; & l^{s} (t) = [\begin{array}{c} l_{1}^{s} (t) \\ ⋮ \\ l_{N}^{s} (t) \end{array}]; & l^{p} = [\begin{array}{c} l_{1}^{p} \\ ⋮ \\ l_{N}^{p} \end{array}]; & F^{g} (t) = [\begin{array}{c} F_{1}^{g} (t) \\ ⋮ \\ F_{N}^{g} (t) \end{array}]; & Δ^{p} = [\begin{array}{c} Δ_{1}^{p} (t) \\ ⋮ \\ Δ_{N}^{p} (t) \end{array}] \end{array}

(23)

Δ k = [\begin{array}{c} {Δ k}_{1} \\ ⋮ \\ {Δ k}_{N} \end{array}]

(24)

This vectorized formulation enables the use of gradient-based optimization algorithms by explicitly expressing the NSD force vector in terms of design variables $Δ k$ , and by ensuring compatibility with the state-space dynamics and numerical integration schemes employed in subsequent sections.

Optimization approach

The proposed methodology aims to optimize the entries of the diagonal matrices $Δ k$ and $Δ c$ , which represent the vertical spring stiffnesses of the NSDs and the damping coefficients of the NFVDs, respectively. The development process is structured into three main stages: (i) formulation of the objective function, (ii) construction of the optimization problem, and (iii) derivation of the solution strategy and implementation algorithm.

Objective function

Before defining the objective function, the full expression of the ANSD control force is established. To this end, the NFVD and NSD force models given in equations (16) and (22) are substituted into equation (12), resulting in the explicit control-law formulation:

\begin{array}{l} u (t) = - diag {Δ k} [\frac{Δ^{p} + l^{p}}{l^{s} (t)} - 1] ⊙ \frac{l^{α}}{l^{β}} ⊙ [2 + \frac{l^{β}}{l^{α}} + \frac{l^{p} + l^{α}}{\sqrt{| {(l^{β})}^{2} - {([\begin{array}{l} I & 0 \end{array}] z (t))}^{2} |}}] ⊙ [\begin{array}{l} I & 0 \end{array}] z (t) + F^{g} (t) + \dots \\ - diag {Δ c} [[\begin{array}{l} 0 & I \end{array}] z (t) ⊙ {| [\begin{array}{l} 0 & I \end{array}] z (t) - D ([\begin{array}{l} 0 & I \end{array}] z (t)) |}^{α - 1}] \end{array}

(25)

Since the state vector $z$ is composed of drift-related coordinates, the control force $u$ directly represents the story-level shear forces applied by the ANSD. This force expression serves as the core of the optimal control formulation and is used to define the following objective function:

J = \int_{0}^{t_{f}} u^{'} (t) Q_{1} u (t) d t + \int_{0}^{t_{f}} z^{'} (t) Q_{2} z (t) d t + \int_{0}^{t_{f}} {f^{H}}^{'} (z) Q_{3} f^{H} (z) d t + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c

(26)

Here, the objective function minimizes the energy of the control forces

u

, the structural state trajectories

z

, and the hysteretic force vector

f^{hys}

, each weighted by matrices

Q_{1}

Q_{2}

, and

Q_{3}

, respectively. Additionally, the cost of implementing the design variables

Δ k

and

Δ c

is penalized through the weighting matrices

R_{k}

and

R_{c}

The total story shear force $V$ comprises three main contributions: the control force $u$ , the inherent damping force, and the resisting force composed of both elastic and hysteretic components. Accordingly, it can be expressed in either of the following equivalent forms:

V (t) = [\begin{array}{l} k^{E} & c^{inh} \end{array}] z (t) + f^{H} (z) - u (t)

(27a)

V (t) = m [\begin{array}{l} 0 & I \end{array}] (A z (t) + B f^{H} (z) - Bu (t))

(27b)

where

V

is the N-dimensional vector of the total story shear force acting on the stories. In the proposed formulation, specific weighting matrices are introduced to relate the objective function directly to the shear force trajectories. These matrices are defined as

\begin{array}{l} Q_{1} = [\begin{array}{l} I & 0 \\ 0 & I \end{array}] & ; & Q_{2} = [\begin{array}{l} Δ k & 0 \\ 0 & Δ c \end{array}] [\begin{array}{l} Δ k & 0 \\ 0 & Δ c \end{array}] & ; & Q_{3} = [\begin{array}{l} I & 0 \\ 0 & I \end{array}] \end{array}

(28a)

\begin{array}{l} R_{k} = [\begin{array}{l} R_{k, 1} & 0 \\ ⋱ \\ 0 & R_{k, N} \end{array}] & ; & R_{c} = [\begin{array}{l} R_{c, 1} & 0 \\ ⋱ \\ 0 & R_{c, N} \end{array}] \end{array}

(28b)

Substituting these definitions into the original formulation, the objective function reduces to

J = \int_{0}^{t_{f}} V^{'} (t) V (t) d t + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c

(29)

This revised form highlights that the optimization aims to minimize the total story shear force trajectories (in an energy norm sense), while simultaneously penalizing excessive or impractical adjustments to the stiffness ( $Δ k$ ) and damping ( $Δ c$ ) parameters. The free parameters $R_{k, 1} \dots, R_{k, N}$ and $R_{c, 1} \dots, R_{c, N}$ , allow tailored tuning of each ANSD component contribution to the overall control effort.

Classical optimal control strategies such as LQR, H₂, and H∞ are often sensitive to the choice of weighting matrices in the cost function, which are typically tuned heuristically to balance structural response and control effort. In this work, the weighting matrices associated with stiffness weakening and energy dissipation ( $R_{k}$ and $R_{c}$ ) are treated as given design preferences with clear mechanical meanings, and the optimization computes the corresponding updates in the ANSD stiffness and damping parameters. This provides a more transparent link between the selected weights and the resulting device configuration, reducing reliance on extensive trial-and-error calibration of abstract weighting coefficients.

Problem formulation

Building upon the objective function defined in equation (29), the complete optimization problem is formulated by introducing the following constraints:

\begin{array}{l} \begin{array}{l} \min_{{Δ k}_{1}, \dots, {Δ k}_{N}, {Δ c}_{1}, \dots, {Δ c}_{N}} {J = \int_{0}^{t_{f}} V^{'} (t) V (t) d t + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c} \\ s . t . \\ \dot{z} (t) = A z (t) + B (- f^{H} (z) + u (t) + p_{r} (t)) \\ p_{r} (t) = - T_{f \to v} m 1 a_{g, r} (t) \\ z (0) = 0 \\ \begin{array}{l} {Δ k}_{n}^{\min} \leq {Δ k}_{n} \leq {Δ k}_{n}^{\max} & \forall & n = 1, \dots, N \end{array} \\ \begin{array}{l} {Δ c}_{n}^{\min} \leq {Δ c}_{n} \leq {Δ c}_{n}^{\max} & \forall & n = 1, \dots, N \end{array} \end{array}} & \forall & r = 1, \dots, R \end{array}

(30)

The constraints comprise the dynamic equilibrium of equation (8), the initial conditions, and the bounds on the design variables.

In the present formulation, the weighting matrices $R_{k}$ and $R_{c}$ are symmetric positive-definite matrices that penalize changes in the vertical spring stiffnesses and damping coefficients, respectively, in the objective function. When these matrices are chosen to be diagonal, each diagonal entry acts as a scalar weight on the corresponding design variable, so that larger values discourage changes in that parameter, while smaller values allow more aggressive updates. In contrast to classical LQR-type formulations, where weighting matrices are tuned indirectly until a satisfactory response is obtained, here $R_{k}$ and $R_{c}$ are interpreted as given design preferences reflecting acceptable parameter variability and implicit device cost, and the optimization computes the corresponding modifications of ${Δ k}_{n}$ and ${Δ c}_{n}$ . In the numerical example $R_{k}$ and $R_{c}$ are taken as diagonal matrices with uniform entries, so that all stories are penalized equally, and these terms act primarily as regularization of the parameter update. Increasing these weights shifts the trade-off toward smaller parameter changes and more conservative device configurations, whereas decreasing them permits larger adjustments and potentially greater reductions of the shear-energy term in the objective.

The dynamic equilibrium is formulated using the first-order state-space representation presented in equation (8), incorporating the horizontal seismic excitation defined in equation (7). Zero initial conditions are imposed, as they are necessary for consistently solving the optimization problem. Additionally, design constraints are introduced to define the permissible lower and upper bounds of the design variables, namely, the vertical spring stiffness coefficients ${Δ k}_{n}$ and the damping coefficients ${Δ c}_{n}$ for each story.

Given that infinitely many combinations of $R_{k}$ and $R_{c}$ can be specified in the objective function, there exists an equally infinite number of corresponding optimal solutions for $Δ k$ and $Δ c$ . To address this challenge, the proposed optimization strategy begins with an initial set of design variables assumed to be optimal. From these, the implied weighting components of $R_{k}$ and $R_{c}$ are inferred. The relative magnitudes of these components are then used to adjust the design variables iteratively. The resulting update rule for refining the design variables based on the derived weighting components is given as follows:

\begin{array}{l} {Δ k}_{n}^{(j + 1)} = {Δ k}_{n}^{(j)} {(I - \frac{R_{k, n}^{(j)}}{\sum | R_{k, 1} | + . . . + | R_{k, N} |})}^{γ_{k}} & \forall & n = 1, \dots, N \end{array}

(31a)

\begin{array}{l} {Δ c}_{n}^{(j + 1)} = {Δ c}_{n}^{(j)} {(I - \frac{R_{c, n}^{(j)}}{\sum | R_{c, 1} | + . . . + | R_{c, N} |})}^{γ_{c}} & \forall & n = 1, \dots, N \end{array}

(31b)

Here,

j

denotes the iteration index, while

γ_{k}

and

γ_{c}

are the convergence parameters governing the update of the stiffness and damping variables, respectively. To compute the components of

R_{k}

and

R_{c}

the Lagrange multiplier method is employed, enabling the integration of design constraints within the optimization framework.

Lagrange multiplier method

To effectively optimize the stiffness and damping parameters for our ANSD, we rigorously apply the principles of variational calculus. The core of this optimization lies in the definition of the Lagrangian, which is meticulously constructed to capture the system’s dynamics and constraints. The Lagrangian, associated with equation (30), is expressed as

L = \int_{0}^{t_{f}} V^{'} (t) V (t) + λ^{'} (t) [A z (t) + B (f^{H} (z) - u (t) - p_{r} (t)) - \dot{z} (t)] d t + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c

(32)

Here,

λ (t)

represents the critical time-varying 2N-dimensional Lagrange multiplier vector, serving to enforce the system’s dynamic constraints.

A crucial step in simplifying this complex formulation involves judicious application of integration by parts. Specifically, by applying this technique to the sub-term $\int_{0}^{t_{f}} λ {(t)}^{T} \dot{z} (t) d t$ we obtain:

\begin{array}{l} L = \int_{0}^{t_{f}} V^{'} (t) V (t) + λ^{'} (t) [A z (t) + B (f^{H} (z) - u (t) - p_{r} (t))] + {\dot{λ}}^{'} (t) z (t) d t + \dots \\ - {λ (t_{f})}^{T} z (t_{f}) + λ {(0)}^{T} z (0) + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c \end{array}

(33)

This strategic manipulation, combined with the practical boundary conditions where $z (0) = 0$ and the imposition $λ (t_{f}) = 0$ (reflecting the terminal and initial states, respectively), profoundly simplifies the Lagrangian to

L = \int_{0}^{t_{f}} V^{'} (t) V (t) + λ^{'} (t) [A z (t) + B (f^{H} (z) - u (t) - p_{r} (t))] + {\dot{λ}}^{'} (t) z (t) d t + {Δ k}^{'} R_{k} Δ k + {Δ c}^{'} R_{c} Δ c

(34)

The minimization of this refined Lagrange function is paramount to our optimization strategy. This is achieved by systematically setting the partial derivatives with respect to the state vector $z$ , the control gains $Δ k$ and $Δ c$ to zero, yielding the following set of governing optimality conditions:

L_{Δ k} = \int_{0}^{t_{f}} 2 V_{Δ k} (t) V (t) - u_{Δ k} (t) B^{'} λ (t) d t + 2 Δ k R_{k} = 0

(35a)

L_{Δ c} = \int_{0}^{t_{f}} 2 V_{Δ c} (t) V (t) - u_{Δ c} (t) B^{'} λ (t) d t + 2 Δ c R_{c} = 0

(35b)

L_{z} = \int_{0}^{t_{f}} 2 V_{z}^{'} (t) V (t) + A^{'} λ (t) + \dot{λ} (t) d t = 0

(35c)

We can derive explicit and highly valuable expressions from these fundamental optimality conditions. Leveraging the relationships $V_{Δ k} (t) = - u_{Δ k} (t)$ and $V_{Δ c} (t) = - u_{Δ c} (t)$ , we elegantly express $\dot{λ}$ and the optimal matrices $R_{k}$ and $R_{c}$ as follows:

\begin{array}{l} \dot{λ} (t) = - 2 V_{z}^{'} (t) V (t) - A^{'} λ (t) & \leftrightarrow & λ (t_{f}) = 0 \end{array}

(36)

R_{k} = 0.5 {Δ k}^{- 1} \int_{0}^{t_{f}} u_{Δ k} (t) (2 V (t) + B^{'} λ (t)) d t

(37a)

R_{c} = 0.5 {Δ c}^{- 1} \int_{0}^{t_{f}} u_{Δ c} (t) (2 V (t) + B^{'} λ (t)) d t

(37b)

Furthermore, the crucial state vector $z$ is determined by solving equation (8). The subsequent sections will detail the advanced numerical evaluation schemes developed for accurately calculating the terms within equations (36), (37a), and (37b), demonstrating the practical implementability of our theoretical framework.

Numerical evaluation of (t)

The precise calculation of $Δ k$ and $Δ c$ presents a formidable challenge. Performing this analysis analytically is simply unfeasible, particularly when confronted with the inherent complexities of unoriented earthquake loading and the formidable uncertainties associated with inelastic structural response. These real-world scenarios introduce layers of complexity that conventional analytical methods cannot robustly address.

To surmount this critical hurdle, we have developed pioneering numerical differentiation equations specifically tailored for the accurate and reliable evaluation of $Δ k$ and $Δ c$ . Our approach masterfully integrates reversed time calculation for the formulation of $λ$ from equation (36) in discrete time steps, leading to

\begin{array}{l} {\dot{λ}}_{i} = - 2 V_{z_{i}}^{'} V_{i} - A^{'} λ_{i} & \forall & i = \frac{t_{f}}{Δ t}, \dots, 1, 0 \end{array}

(38)

Here,

i

denotes the time-step index. Furthermore, the calculation of the partial derivative

V_{z}

is performed with unprecedented stability using the extended mean-value theorem:

V_{z_{i}} = {\begin{cases} (V_{i} - V_{i - 1}) {(z_{i} - z_{i - 1})}^{'}^{- 1} & \forall & i = \frac{t_{f}}{Δ t} \\ 2 (V_{i + 1} - V_{i}) {(z_{i + 1} - z_{1, i})}^{'}^{- 1} - V_{z_{i + 1}} & \forall & i = \frac{t_{f}}{Δ t} - 1, \dots, 1, 0 \end{cases}

(39)

It is crucial to emphasize that the strategic application of the extended mean-value theorem is a cornerstone of our method’s robustness, guaranteeing superior numerical stability compared to other explicit methods, such as the standard finite difference method, which are often prone to instability in such complex scenarios.

Our optimization framework critically requires that $λ (t_{f}) = 0$ . To meet this stringent condition and ensure a complete numerical solution, the Lagrange multiplier vector at the final time step is meticulously calculated using the following scheme:

λ_{i} = {\begin{array}{c} 0 & \forall & i = \frac{t_{f}}{Δ t} \\ λ_{i + 1} - \frac{1}{2} Δ t ({\dot{λ}}_{i} + {\dot{λ}}_{i + 1}) & \forall & i = \frac{t_{f}}{Δ t} - 1, \dots, 1, 0 \end{array}

(40)

Equation (40) is not merely a calculation; it represents a key enabler for the practical implementation of our sophisticated optimization. It provides the definitive numerical scheme for accurately evaluating $λ_{i}$ , completing the necessary conditions for our analytical framework. With these powerful numerical schemes firmly established, we now proceed to outline the comprehensive solution algorithm that fully leverages these advancements.

Optimization procedure

A highly efficient and robust optimization algorithm is developed to translate our sophisticated theoretical framework into practical, implementable solutions. This procedure is meticulously designed to navigate the complex interplay of seismic forces and device parameters, yielding optimal stiffness and damping for Adaptive Negative Stiffness Devices. The algorithm is strategically divided into three interconnected stages, each playing a vital role in ensuring a comprehensive and convergent solution:

(i) Algorithm Preparations

(ii) Earthquake Response Calculation

(iii) Iterative optimization of $Δ k$ and $Δ c$

This structured, iterative approach is not merely a sequence of steps; it represents a powerful, self-correcting mechanism that converges on the most effective stiffness and damping profiles. The detailed process is fully described herein, demonstrating the sophistication and efficacy of our proposed methodology.

1. Algorithm preparations

1.1 Choose the load cases of R ground acceleration records $a_{g, 1} (t), \dots a_{g, R} (t)$ .

1.2 Define the initial NSD vertical spring stiffness coefficients and their design limitations:

\begin{array}{l} {Δ k}^{(0)} = [\begin{array}{c} {Δ k}_{1}^{(0)} \\ ⋮ \\ {Δ k}_{N}^{(0)} \end{array}] & ; & {Δ k}^{\min} = [\begin{array}{c} {Δ k}_{1}^{\min} \\ ⋮ \\ {Δ k}_{N}^{\min} \end{array}] & ; & {Δ k}^{\max} = [\begin{array}{c} {Δ k}_{1}^{\max} \\ ⋮ \\ {Δ k}_{N}^{\max} \end{array}] \end{array}

1.3 Define the initial NFVD damping coefficients and their design limitations:

\begin{array}{l} {Δ c}^{(0)} = [\begin{array}{c} {Δ c}_{1}^{(0)} \\ ⋮ \\ {Δ c}_{N}^{(0)} \end{array}] & ; & {Δ c}^{\min} = [\begin{array}{c} {Δ c}_{1}^{\min} \\ ⋮ \\ {Δ c}_{N}^{\min} \end{array}] & ; & {Δ c}^{\max} = [\begin{array}{c} {Δ c}_{1}^{\max} \\ ⋮ \\ {Δ c}_{N}^{\max} \end{array}] \end{array}

1.4 Define the convergence coefficients $γ_{k}$ and $γ_{c}$ , convergence tolerance $T O L$ , and the maximum number of iterations $j^{\max}$ .

1.5 Set the iteration index to $j = 0$ and the record number of most severe earthquake response to $ρ^{o l d} = 0$ .

2. Earthquake response calculation for all ground acceleration records

2.1 Perform numerical evaluation for $z_{r}^{(j)} (t)$ and $V_{r}^{(j)} (t)$ :

\begin{array}{l} {f^{H}}_{r}^{(j)} (δ, \dot{δ}) + k^{E} δ_{r}^{(j)} (t) + c^{inh} {\dot{δ}}_{r}^{(j)} (t) + m {\ddot{δ}}_{r}^{(j)} (t) = - T_{f \to v} m 1 a_{g, r} (t) + u_{r}^{(j)} (t) \\ ↓ \\ z_{r}^{(j)} (t) = {[\begin{array}{l} δ_{r}^{(j)} (t) & {\dot{δ}}_{r}^{(j)} (t) \end{array}]}^{'} \\ ↓ \\ V_{r}^{(j)} (t) = [\begin{array}{l} k^{E} & c^{inh} \end{array}] z_{r}^{(j)} (t) + {f^{H}}_{r}^{(j)} (z) - u_{r}^{(j)} (t) \end{array}

2.2 Evaluate $λ_{r}^{(j)} (t)$ using equation (40).

2.3 Derive ${(R_{k}^{(j)})}_{r}$ and ${(R_{c}^{(j)})}_{r}$ according to equations (37a) and (37b).

2.3 Calculate the objective function:

J_{r}^{(j)} = \int_{0}^{t_{f}} {V_{r}^{(j)}}^{'} (t) V_{r}^{(j)} (t) d t + {Δ k}^{(j)}^{'} {(R_{k}^{(j)})}_{r} {Δ k}^{(j)} + {Δ c}^{(j)}^{'} {(R_{c}^{(j)})}_{r} {Δ c}^{(j)}

2.5 Take the largest objective function by

\begin{array}{l} J_{r}^{(j)} > \max_{i = 1, \dots, r - 1} {J_{i}^{(j)}} & \to & ρ = r & & & J_{\max}^{(j)} = J_{r}^{(j)} \end{array}

2.6 If $ρ = ρ^{o l d}$ then Finish. Else, proceed to stage 3.

3. Iterative optimization of $Δ k$ and $Δ c$

3.1 Update ${Δ k}^{(j)}$ and ${Δ c}^{(j)}$ :

\begin{array}{l} \begin{array}{l} {Δ k}_{n}^{(j *)} = {Δ k}_{n}^{(j)} {(1 - \frac{R_{k, n}^{(j)}}{\sum | R_{k, 1}^{(j)} | + . . . + | R_{k, N}^{(j)} |})}^{γ_{k}} & \forall & n = 1, \dots, N \end{array} & \to & {Δ k}_{n}^{(j *)} = {\begin{cases} {Δ k}_{n}^{(j *)} \leq {Δ k}_{n}^{\min} : & {Δ k}_{n}^{\min} \\ {Δ k}_{n}^{(j *)} \geq {Δ k}_{n}^{\max} : & {Δ k}_{n}^{\max} \end{cases} \end{array}

\begin{array}{l} \begin{array}{l} {Δ c}_{n}^{(j *)} = {Δ c}_{n}^{(j)} {(1 - \frac{R_{c, n}^{(j)}}{\sum | R_{c, 1}^{(j)} | + . . . + | R_{c, N}^{(j)} |})}^{γ_{c}} & \forall & n = 1, \dots, N \end{array} & \to & {Δ c}_{n}^{(j *)} = {\begin{cases} {Δ c}_{n}^{(j *)} \leq {Δ c}_{n}^{\min} : & {Δ c}_{n}^{\min} \\ {Δ c}_{n}^{(j *)} \geq {Δ c}_{n}^{\max} : & {Δ c}_{n}^{\max} \end{cases} \end{array}

3.2 Perform numerical evaluation for $z^{(j *)} (t)$ and $V^{(j *)} (t)$ :

\begin{array}{l} f_{ρ}^{H (j *)} (δ, \dot{δ}) + k^{E} δ_{ρ}^{(j *)} (t) + c^{inh} {\dot{δ}}_{ρ}^{(j *)} (t) + m {\ddot{δ}}_{ρ}^{(j *)} (t) = - T_{f \to v} m 1 a_{g, ρ} (t) + u_{ρ}^{(j *)} (t) \\ ↓ \\ z_{ρ}^{(j *)} (t) = {[\begin{array}{l} δ_{ρ}^{(j *)} (t) & {\dot{δ}}_{ρ}^{(j *)} (t) \end{array}]}^{'} \\ ↓ \\ V_{ρ}^{(j *)} (t) = [\begin{array}{l} k^{E} & c^{inh} \end{array}] z_{ρ}^{(j *)} (t) + {f^{H}}_{ρ}^{(j *)} (z) - u_{ρ}^{(j *)} (t) \end{array}

3.3 Evaluate $λ_{ρ}^{(j)} (t)$ using equation (40).

3.4 Derive ${(R_{k}^{(j)})}_{ρ}$ and ${(R_{c}^{(j)})}_{ρ}$ according to equations (37a) and (37b).

3.5 Calculate the objective function:

J_{ρ}^{(j *)} = \int_{0}^{t_{f}} {V_{ρ}^{(j *)}}^{'} (t) V_{ρ}^{(j *)} (t) d t + {Δ k}^{(j *)}^{'} R_{k}^{(j *)} {Δ k}^{(j *)} + {Δ c}^{(j *)}^{'} R_{c}^{(j *)} {Δ c}^{(j *)}

3.6 Verify that the shear force trajectory has been improved:

\int_{0}^{t_{f}} {V_{ρ}^{(j *)}}^{'} (t) V_{ρ}^{(j *)} (t) d t < \int_{0}^{t_{f}} {V_{ρ}^{(j)}}^{'} (t) V_{ρ}^{(j)} (t) d t

If not, update the iteration index “ $j + 1 \to j$ ,” and go back to 2.0 with the a priori design variables ${Δ k}^{(j)}$ and ${Δ c}^{(j)}$ .

If yes, check the convergence and update the matrices:

\begin{array}{l} 1 - 0.5 (\frac{{Δ k}^{(j *)}^{'} {Δ k}^{(j *)}}{{Δ k}^{(j)}^{'} {Δ k}^{(j)}} + \frac{{Δ c}^{(j *)}^{'} {Δ c}^{(j *)}}{{Δ c}^{(j)}^{'} {Δ c}^{(j)}}) < TOL & \to & \begin{array}{l} \begin{array}{l} {Δ k}_{n}^{(j *)} \to {Δ k}_{n}^{(j)} \\ {Δ c}_{n}^{(j *)} \to {Δ c}_{n}^{(j)} \end{array} & \forall & n = 1, \dots, N \end{array} \end{array}

If satisfied, update the most severe earthquake record number “ $ρ^{o l d} = ρ$ ,” increase iteration index “ $j + 1 \to j$ ,” and go to 2.0.

If not, go to 3.1.

Case study

A numerical case study is conducted on a reinforced concrete moment-resisting frame that forms part of the lateral force-resisting system of an eight-story residential building in Ramat-Gan, Israel. The frame is modeled as a planar eight-story system. Figure 5(a) shows the elevation scheme, and Figure 5(b) presents a photograph of the actual building.

Figure 5.

Studied moment-resisting frame: (a) elevation scheme and (b) photo of the building containing the planar frame.

The structural model uses the actual story heights, bay spans, and lumped story masses derived from the building data. Beams and columns are represented by 30 cm × 60 cm rectangular cross-sections, and their longitudinal and transverse reinforcement layouts follow the drawings in Figure 6. Axial force–bending moment interaction diagrams for both beams and columns, as shown in Figure 7, are used in post-processing to assess demand–capacity ratios and define a collapse indicator; however, they do not directly modify the equations of motion in the dynamic analysis.

Figure 6.

Reinforcement design: (a) beam and (b) column.

Figure 7.

Cross-section eccentric compression interaction curve: (a) RC beam and (b) RC column.

The dynamic properties of the bare frame are characterized by the first three vibration modes, whose periods and mode shapes are reported in Figure 8. A 5% inherent damping ratio is assigned to all modes using a classical viscous damping model, which produces approximately uniform modal damping and represents typical energy dissipation in reinforced concrete buildings.

Figure 8.

Moment-resisting frame: the three most dominant modal system periods and shapes.

Seismic demand is represented by an ensemble of 15 recorded ground acceleration histories, summarized in Table 1. The records are selected to be compatible with the seismic hazard of Ramat-Gan (2% probability of exceedance in 50 years) for relatively soft soil conditions and were previously scaled to match the target elastic response spectrum at the site. For each record, Figure 9 illustrates representative dynamic responses of the bare frame. Table 2 reports the maximum roof displacement, maximum base shear, peak absolute acceleration, and collapse flag. In particular, the El Mayor–Cucapah record at station Chihuahua (record r = 8) produces a collapse mechanism characterized by bending failure at both ends of all columns in the third story. In contrast, the remaining records drive the structure into a severe but non-collapsing inelastic response.

Table 1.

Ground acceleration records user for the case study.

Rec. r	Earthquake name	Recording station	Mag.	PGA (g)	PGV (cm/s)	PGD (cm)
1	Landers	Morongo Valley Fire Station	7.28	0.164	22.22	10.19
2	Kocaeli	Arcelik	7.51	0.134	40.05	37.08
3	Kocaeli	Gebze	7.51	0.144	32.63	29.76
4	Kocaeli	Lamont1059	7.14	0.152	12.86	10.19
5	Kocaeli	Lamont1059	7.14	0.137	10.32	8.47
6	Kocaeli	Lamont1061	7.14	0.131	12.13	8.79
7	Kocaeli	Lamont531	7.14	0.160	10.84	7.65
8	El Mayor-Cucapah	Chihuahua	7.2	0.197	34.02	31.21
9	El Mayor-Cucapah	Ejido Saltillo	7.2	0.152	45.60	38.96
10	El Mayor-Cucapah	Ejido Saltillo	7.2	0.148	36.72	38.67
11	Darfield	Canterbury Aero club	7	0.198	47.46	45.73
12	Darfield	Canterbury Aero club	7	0.186	34.04	32.49
13	Darfield	Riccarton High school	7	0.190	30.53	19.83
14	Darfield	RKAC	7	0.191	17.89	11.72
15	Darfield	RKAC	7	0.167	10.55	3.62

Figure 9.

Stories n = 1, 3, 8 response plots under selected ground acceleration records r = 8, 10, 13 (under initial design variables configuration).

Table 2.

Earthquake response of the bare frame under the 15 ground acceleration records.

Rec. r	Max. roof displacement (m)	Max. base shear (kN)	Peak absolute acceleration (g)	Collapsed?
1	0.225	505	0.311	No
2	0.016	308	0.209	No
3	0.025	388	0.324	No
4	0.017	353	0.296	No
5	0.011	338	0.328	No
6	0.012	265	0.248	No
7	0.024	447	0.361	No
8	0.238	522	0.419	Yes
9	0.023	440	0.275	No
10	0.200	481	0.303	No
11	0.021	450	0.265	No
12	0.017	376	0.327	No
13	0.032	475	0.406	No
14	0.025	389	0.346	No
15	0.017	398	0.380	No

The performance measures examined in the case study include maximum inter-story drift at each story, roof displacement, base shear, peak absolute floor accelerations, and an energy-type objective functional based on the story-shear force trajectories. Inter-story drifts are primarily used to assess damage and potential collapse, while roof displacement and base shear describe the global response, and floor accelerations relate to non-structural demands.

The ANSD configuration adopted in the case study follows the prototype geometry tested by Sarlis et al. (2016). The lever-arm lengths, vertical spring initial length and pre-compression, and the two-stage horizontal spring assembly at the device base are defined as in Section 3. The viscous damper nonlinearity parameter α is kept constant, while the diagonal matrices of vertical spring stiffnesses and damping coefficients form the design variables. Initial values and bounds for these parameters, listed in Tables 3 and 4, are selected such that the initial ANSD effectiveness is negligible and the supplemental damping contribution is approximately 20% in the first modal system.

Table 3.

NFVD coefficient design variables in each iteration.

Story N	j = 0 initial design (kNs/m)	j = 1 optimized under record ρ = 13 (kNs/m)	j = 2 optimized under record ρ = 10 (kNs/m)	j = 3 optimized under record ρ = 13 (kNs/m)
8	400	392	396	396
7	400	400	400	400
6	400	396	396	396
5	400	408	416	416
4	400	404	404	404
3	400	404	412	412
2	400	404	428	428
1	400	412	412	412

Table 4.

NSD vertical spring stiffness design variables in each iteration.

Story n	j = 0 initial design (kN/m)	j = 1 optimized under EQ ρ = 13 (kN/m)	j = 2 optimized under EQ ρ = 10 (kN/m)	j = 3 optimized under EQ ρ = 13 (kN/m)
8	20,000	19,244	19,344	19,344
7	20,000	19,948	19,936	19,936
6	20,000	19,832	19,768	19,768
5	20,000	20,552	20,756	20,756
4	20,000	20,496	20,396	20,396
3	20,000	20,184	20,556	20,556
2	20,000	20,044	21,332	21,332
1	20,000	20,232	20,236	20,236

The optimization algorithm described in Section 4 is then applied to this frame. Following initialization, the bare-frame response to all 15 records is computed, and the objective-function value is evaluated for each case. Table 5 compiles the corresponding objective-function values, along with peak response measures, and identifies record r = 13 as the most severe excitation according to the shear-force-based objective, which considers the entire story-shear time history rather than just its maximum value.

Table 5.

Maximum responses of the bare frame under the 15 ground acceleration records at iteration j = 0 (before introducing the design variables).

Rec. r	Max. roof displacement (m)	Max. base shear (kN)	Peak absolute acceleration (g)	Objective function (kN²)
1	0.017	434	0.209	1.95 E + 10
2	0.007	234	0.209	3.68 E + 08
3	0.013	377	0.209	8.15 E + 09
4	0.005	183	0.266	1.00 E + 06
5	0.005	181	0.265	1.10 E + 10
6	0.008	245	0.265	5.43 E + 09
7	0.009	262	0.266	1.62 E + 10
8	0.021	430	0.267	6.47 E + 08
9	0.010	290	0.209	2.95 E + 10
10	0.012	347	0.209	3.99 E + 10
11	0.011	333	0.209	1.06 E + 11
12	0.010	305	0.209	1.18 E + 11
13	0.016	378	0.240	3.82 E + 11
14	0.012	335	0.223	9.53 E + 09
15	0.008	238	0.209	1.84 E + 09

Subsequent iterations update the design variables and re-analyze the frame until no further reduction in the shear-force trajectory is obtained or the maximum number of iterations is reached. For this case study, the first optimization attempt targets the severe record ρ = 13, the second targets the record ρ = 10, and the third again targets the record ρ = 13; no additional critical records emerge thereafter. Table 6 summarizes the evolution of the shear-force trajectory component of the objective function across these attempts (see the ρ = 10 and ρ = 13 iterations in bold fonts). For most records, this quantity decreases monotonically, confirming the effectiveness of the iterative procedure. Records r = 4 and r = 5 exhibit essentially unchanged, low shear-force trajectories, whereas record r = 8 shows a higher value due to the later onset of the collapse mechanism in the optimized configuration.

Table 6.

Objective function values during the iterative process.

Rec. r	j = 0 initial design	j = 1 optimized under severe record ρ = 13		j = 2 optimized under severe record ρ = 10		j = 3 optimized under severe record ρ = 13
1	1.19 E + 06	1.17 E + 06	(Lower)	1.14 E + 06	(Lower)	1.14 E + 06	(Unchanged)
2	2.19 E + 05	2.17 E + 05	(Lower)	2.16 E + 05	(Lower)	2.16 E + 05	(Unchanged)
3	3.58 E + 05	3.54 E + 05	(Lower)	3.48 E + 05	(Lower)	3.48 E + 05	(Unchanged)
4	1.27 E + 05	1.27 E + 05	(Unchanged)	1.27 E + 05	(Unchanged)	1.27 E + 05	(Unchanged)
5	1.08 E + 05	1.08 E + 05	(Unchanged)	1.08 E + 05	(Unchanged)	1.08 E + 05	(Unchanged)
6	1.37 E + 05	1.36 E + 05	(Lower)	1.35 E + 05	(Lower)	1.35 E + 05	(Unchanged)
7	2.83 E + 05	2.79 E + 05	(Lower)	2.74 E + 05	(Lower)	2.74 E + 05	(Unchanged)
8	2.10 E + 06	2.09 E + 06	(Lower)	2.19 E + 06	(Higher)	2.19 E + 06	(Unchanged)
9	1.72 E + 06	1.70 E + 06	(Lower)	1.69 E + 06	(Lower)	1.69 E + 06	(Unchanged)
10	1.67 E + 06	1.65 E + 06	(Lower)	1.62 E + 06	(Lower)	1.62 E + 06	(Unchanged)
11	1.02 E + 06	1.01 E + 06	(Lower)	1.00 E + 06	(Lower)	1.00 E + 06	(Unchanged)
12	6.95 E + 05	6.88 E + 05	(Lower)	6.79 E + 05	(Lower)	6.79 E + 05	(Unchanged)
13	9.65 E + 05	9.46 E + 05	(Lower)	9.22 E + 05	(Lower)	9.22 E + 05	(Unchanged)
14	3.25 E + 05	3.21 E + 05	(Lower)	3.16 E + 05	(Lower)	3.16 E + 05	(Unchanged)
15	3.14 E + 05	3.09 E + 05	(Lower)	3.02 E + 05	(Lower)	3.02 E + 05	(Unchanged)

Tables 3 and 4 report the story-wise evolution of the nonlinear fluid viscous damper coefficients and negative stiffness device vertical spring stiffnesses, respectively, and show that the optimized ANSD design remains within practical parameter ranges with moderate variations across stories. Table 7 provides the final maximum roof displacements, base shears, peak absolute accelerations, maximum inter-story drifts, and objective-function values for all records after convergence. Compared with the bare-frame response in Table 2, the optimized ANSD configuration generally reduces the shear-energy measure and moderates peak responses for the most demanding records, while preserving the global collapse pattern associated with record r = 8. Figure 10 presents representative hysteretic loops for stories 1, 3, and 8, illustrating the improved energy dissipation and the redistribution of inelastic demand achieved by the proposed optimization strategy.

Table 7.

Maximum responses of the bare frame under the 15 ground acceleration records at iteration j = 3 (following optimization).

Rec. r	Max. roof displacement (m)	Max. base shear (kN)	Peak absolute acceleration (g)	Objective function (kN²)
1	0.0208	404	0.216	2.54 E + 10
2	0.0083	237	0.218	1.23 E + 09
3	0.0145	356	0.222	6.61 E + 09
4	0.0054	171	0.299	4.06 E + 08
5	0.0055	186	0.299	4.44 E + 09
6	0.0071	215	0.299	3.3 E + 09
7	0.0087	250	0.299	4.01 E + 08
8	0.1574	367	0.215	8.34 E + 08
9	0.0113	293	0.216	7.69 E + 09
10	0.0124	320	0.216	6.44 E + 09
11	0.0132	329	0.215	4.28 E + 10
12	0.0124	313	0.216	9.71 E + 10
13	0.0154	329	0.215	3.88 E + 11
14	0.0113	294	0.216	5.03 E + 09
15	0.0087	240	0.215	2.82 E + 08

Figure 10.

Story n = 1, 3, 8 responses under ground acceleration records r = 8, 10, 13 under final design variables configuration.

In this case study, only moderate adjustments to the ANSD parameters are introduced. The optimized vertical spring stiffnesses and damping coefficients remain close to their initial design values, with changes on the order of approximately 5–10%. Thus, the optimization primarily acts as a fine-tuning procedure for an already feasible ANSD configuration rather than prescribing a radically different device layout. Within the chosen parameter bounds, this means that the final design is primarily driven by the gradients of the shear-energy objective and is expected to be only weakly sensitive to modest perturbations of the initial ANSD values. That is, different starting configurations with negligible device effects would lead to similar tuned patterns. From a practical standpoint, such changes can be implemented by adjusting pre-compression levels or spring and damper characteristics within typical manufacturing tolerances, without requiring a redesign of the frame or replacement of the devices. At the same time, the example highlights a limitation: if the performance objectives were to demand changes beyond these bounds, the current framework would still indicate the required trends in stiffness weakening and energy dissipation, but the resulting design might be challenging to realize without more substantial modifications to the ANSD hardware or the underlying structural system, and a dedicated sensitivity study with respect to the initial design would then be required.

For record r = 8, the bare frame develops a concentration of inelastic demand and large interstory drift in the third story, leading to bending moments that drive the demand points outside the interaction diagrams for all columns at that level and triggering the collapse indicator. In the tuned ANSD configuration, the additional negative-stiffness and viscous-damping contributions reduce and redistribute the story shears and drifts at this story, so that the peak bending moments and axial–bending interaction ratios remain within the capacity envelope. Consequently, the collapse indicator is not activated for record r = 8, and the tuned ANSD design is interpreted as preventing the collapse mechanism observed in the bare frame.

Summary and conclusions

This paper has presented a time-domain optimal design framework for tuning adaptive negative stiffness devices (ANSDs) in inelastic frame structures subjected to earthquake excitation. The structure is modeled as a nonlinear shear building with hysteretic restoring forces, with ANSDs—each composed of a negative stiffness device and a nonlinear fluid viscous damper—placed at every story. A state-space representation suitable for control and optimization is derived, and a control-inspired objective functional is defined that aggregates squared story-shear trajectories over an ensemble of ground motions, together with terms that penalize changes in ANSD stiffness and damping coefficients. An adjoint (Lagrange multiplier) formulation provides the gradients of this objective with respect to the device parameters, which are evaluated numerically in discrete time and embedded in an iterative parameter-update algorithm.

The framework is applied to an eight-story reinforced concrete frame representing an existing building in Ramat-Gan, Israel, idealized as a planar shear system and subjected to fifteen recorded ground motions scaled to the site-specific hazard. For each record, the analysis evaluates roof displacement, base shear, peak floor acceleration, interstory drift at selected stories and a collapse indicator derived from section interaction diagrams. Starting from an initial configuration in which the ANSDs have negligible effect, the procedure updates the vertical spring stiffnesses and damping coefficients using the earthquakes that yield the largest values of the shear-energy objective.

The numerical results indicate that, for most ground motions, the tuning procedure reduces the shear-energy measure and yields moderate decreases in peak base shear and floor accelerations for the critical records, while keeping changes in device parameters relatively small. In the considered example, the tuned design also prevents a collapse that occurs for one record in the absence of ANSDs. At the same time, the focus on global shear-energy minimization rather than explicit code-based limits allows interstory drifts and response distributions to vary non-uniformly; therefore, the formulation should be viewed primarily as a control-inspired tuning tool rather than a complete performance-based design procedure.

Several modeling simplifications suggest directions for further work. The structural model is planar and shear-type, without soil–structure interaction, torsional effects, or higher-mode localization, and the nonlinear fluid viscous dampers are represented as pure dashpots. The optimization objective does not enforce explicit constraints on drift, acceleration or member demand–capacity ratios. Future research will extend the framework to multi-objective and constrained formulations, incorporate more detailed three-dimensional and soil–structure interaction models, and investigate the experimental or hybrid-simulation validation of tuned ANSD designs.

In this study, the structure is idealized as a two-dimensional moment-resisting frame, and torsional effects, soil–structure interaction, and higher-mode localization in the transverse direction are neglected. In asymmetric or irregular buildings, these phenomena may alter the quantitative seismic response; therefore, the present results should be interpreted as representative of the dominant in-plane translational behavior of the examined frame. It is noted, however, that such planar frame idealizations, combined with codified provisions for torsional amplification (for example, the accidental eccentricity factors in EN 1998-1:2004, 4.3.3.2.4), are consistent with common seismic design practice, and future work will extend the proposed ANSD optimization to more detailed three-dimensional models including torsion, soil–structure interaction and coupled flexural or flexure–shear frame behavior.

Footnotes

ORCID iD

A. Shmerling

Ethical considerations

This article does not contain any studies with human participants or animals performed by any of the authors. Ethical approval was therefore not required.

Funding

The author received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data Availability Statement

The numerical models and simulation codes used during the current study are available from the corresponding author on reasonable request. Ground motion records employed in the analyses are available from public strong-motion databases as cited in the manuscript.*

References

Cimellaro

(2009) Optimal weakening and damping using polynomial control for seismically excited nonlinear structures. Earthquake Engineering and Engineering Vibration 8(4): 607–616. https://doi.org/10.1007/s11803-009-9124-2

Etedali

Tavakoli

(2017) PD/PID controller design for seismic control of high-rise buildings using multi-objective optimization: a comparative study with LQR controller. Journal of Earthquake and Tsunami 11(3): 1750009. https://doi.org/10.1142/S1793431117500099

Gao

(2024) Supplemental damping and energy dissipation of structural systems installed with viscous dampers considering connection stiffness. Structures 70: 107629. https://doi.org/10.1016/j.istruc.2024.107629

Gao

Xing

Wang

, et al. (2023) Performance improvement and demand-oriented optimum design of the tuned negative stiffness inerter damper for base-isolated structures. Journal of Building Engineering 63: 105488. https://doi.org/10.1016/j.jobe.2022.105488

Gluck

Reinhorn

Gluck

, et al. (1996) Design of supplemental dampers for control of structures. Journal of Structural Engineering 122(12): 1394–1399. https://doi.org/10.1061/(ASCE)0733-9445(1996)122:12(1394)

GW-M

(2024) Some techniques in optimal control of building frames. In: Ho

GW-M

(ed). Structural Dynamic Systems: Computational Techniques and Optimization. Routledge, 99–162.

Zhao

Barredo

, et al. (2024) Nonlinear viscous damping-based negative stiffness isolation system for over-track complex structures. Journal of Earthquake Engineering 28(16): 3146–3174. https://doi.org/10.1080/13632469.2024.2328643

Idels

Lavan

(2021) Performance-based seismic retrofitting of frame structures using negative stiffness devices and fluid viscous dampers via optimization. Earthquake Engineering & Structural Dynamics 50(12): 3116–3137. https://doi.org/10.1002/eqe.3502

Islam

Jangid

(2023) Closed form expressions for H2 optimal control of negative stiffness and inerter-based dampers for damped structures. Structures 50: 791–809. https://doi.org/10.1016/j.istruc.2023.02.065

10.

Kelly

Leitmann

Soldatos

(1987) Robust control of base-isolated structures under earthquake excitation. Journal of Optimization Theory and Applications 53(2): 159–180. https://doi.org/10.1007/BF00939213

11.

(2020) Negative stiffness devices for vibration isolation applications: a review. Advances in Structural Engineering 23(8): 1739–1755. https://doi.org/10.1177/1369433219900311

12.

Kim

Song

(2023) Bouc–Wen class models considering hysteresis mechanism of RC columns in nonlinear dynamic analysis. International Journal of Non-Linear Mechanics 148: 104263. https://doi.org/10.1016/j.ijnonlinmec.2022.104263

13.

Özuygur

Gündüz

(2018) Optimal control of structures under earthquakes including soil–structure interaction. Journal of Earthquake Engineering 22(8): 1317–1335. https://doi.org/10.1080/13632469.2016.1277567

14.

Pohoryles

Duffour

(2015) Adaptive control of structures under dynamic excitation using magnetorheological dampers: an improved clipped-optimal control algorithm. Journal of Vibration and Control 21(13): 2569–2582. https://doi.org/10.1177/1077546313510543

15.

Sarlis

Pasala

DTR

Constantinou

, et al. (2011) Negative stiffness device for seismic protection of structures: an analytical and experimental study. In Proceedings of the 3rd ECCOMAS thematic conference on computational methods in structural dynamics and earthquake engineering (COMPDYN 2011), 25-28 May, 2011, Corfu, Greece.

16.

Sarlis

Pasala

DTR

Constantinou

, et al. (2016) Negative stiffness device for seismic protection of structures: shake table testing of a seismically isolated structure. Journal of Structural Engineering 142(5): 04016005. https://doi.org/10.1061/(ASCE)ST.1943-541X.0001455

17.

Shmerling

Gerdts

(2022) A design methodology for the seismic retrofitting of existing frame structures post-earthquake incident using nonlinear control systems. Buildings 12(11): 1886. https://doi.org/10.3390/buildings12111886

18.

Shmerling

Gerdts

(2023) Short-horizon acceleration-predictive control for reducing lateral seismic inertia forces of inelastic frame structures using semi-active fluid viscous dampers. Computers & Structures 281: 107032. https://doi.org/10.1016/j.compstruc.2023.107032

19.

Shmerling

Levy

Reinhorn

(2018) Seismic retrofit of frame structures using passive systems based on optimal control. Structural Control and Health Monitoring 25(3): e2038. https://doi.org/10.1002/stc.2038

20.

Takewaki

(1997) Optimal damper placement for minimum transfer functions. Earthquake Engineering & Structural Dynamics 26(11): 1113–1124. https://doi.org/10.1002/(SICI)1096-9845(199711)26:11

21.

Takewaki

Akehashi

(2021) Comprehensive review of optimal and smart design of nonlinear building structures with and without passive dampers subjected to earthquake loading. Frontiers in Built Environment 7: 631114. https://doi.org/10.3389/fbuil.2021.631114

22.

Takewaki

Yoshitomi

(1998) Effects of support stiffnesses on optimal damper placement for a planar building frame. The Structural Design of Tall Buildings 7(4): 323–336. https://doi.org/10.1002/(SICI)1099-1794(199812)7:4

23.

Tehrani

Harvey

(2019) Enhanced passive control of dual-mode systems under extreme seismic loading: an optimal control approach. Structural Control and Health Monitoring 26(3): e2367. https://doi.org/10.1002/stc.2367

24.

Wang

Sun

Nagarajaiah

(2019) Simplified optimal design of MDOF structures with negative stiffness amplifying dampers based on effective damping. The Structural Design of Tall and Special Buildings 28(1): e1664. https://doi.org/10.1002/tal.1664

25.

Wang

Gao

, et al. (2021) Optimum design and performance evaluation of the tuned inerter-negative-stiffness damper for seismic protection of single-degree-of-freedom structures. International Journal of Mechanical Sciences 212: 106805. https://doi.org/10.1016/j.ijmecsci.2021.106805

26.

Zakian

Kaveh

(2023) Seismic design optimization of engineering structures: a comprehensive review. Acta Mechanica 234(4): 1305–1330. https://doi.org/10.1007/s00707-022-03470-6

27.

Zhang

Pan

Zeng

, et al. (2020) A novel robust optimum control algorithm and its application to semi-active controlled base-isolated structures. Bulletin of Earthquake Engineering 18(6): 2431–2460. https://doi.org/10.1007/s10518-019-00761-7

Story N	j = 0 initial design (kNs/m)	j = 1 optimized under record ρ = 13 (kNs/m)	j = 2 optimized under record ρ = 10 (kNs/m)	j = 3 optimized under record ρ = 13 (kNs/m)
8	400	392	396	396
7	400	400	400	400
6	400	396	396	396
5	400	408	416	416
4	400	404	404	404
3	400	404	412	412
2	400	404	428	428
1	400	412	412	412

Story N	j = 0 initial design (kNs/m)	j = 1 optimized under record ρ = 13 (kNs/m)	j = 2 optimized under record ρ = 10 (kNs/m)	j = 3 optimized under record ρ = 13 (kNs/m)
8	400	392	396	396
7	400	400	400	400
6	400	396	396	396
5	400	408	416	416
4	400	404	404	404
3	400	404	412	412
2	400	404	428	428
1	400	412	412	412

Story N	j = 0 initial design (kNs/m)	j = 1 optimized under record ρ = 13 (kNs/m)	j = 2 optimized under record ρ = 10 (kNs/m)	j = 3 optimized under record ρ = 13 (kNs/m)
8	400	392	396	396
7	400	400	400	400
6	400	396	396	396
5	400	408	416	416
4	400	404	404	404
3	400	404	412	412
2	400	404	428	428
1	400	412	412	412