Sensitivity Analysis for Hybrid Systems and Systems with Memory

TL;DR

This paper introduces an adjoint sensitivity analysis method for hybrid systems with memory, enabling efficient computation of sensitivities in complex systems with mode switching and history dependence.

Contribution

It extends existing sensitivity methods to hybrid systems with memory and provides solvability conditions for mode transitions and state initialization.

Findings

Developed an adjoint sensitivity method for hybrid systems with memory.

Provided solvability conditions for mode switching and initialization.

Demonstrated the method with numerical examples including a mechanics problem.

Abstract

We present an adjoint sensitivity method for hybrid discrete -- continuous systems, extending previously published forward sensitivity methods. We treat ordinary differential equations and differential-algebraic equations of index up to two (Hessenberg) and provide sufficient solvability conditions for consistent initialization and state transfer at mode switching points, for both the sensitivity and adjoint systems. Furthermore, we extend the analysis to so-called hybrid systems with memory where the dynamics of any given mode depend explicitly on the states at the last mode transition point. We present and discuss several numerical examples, including a computational mechanics problem based on the so-called exponential model constitutive material law for steel reinforcement under cyclic loading.

Equations

$$F\left(\dot{x},x,p,t\right)=0$$

$$T^{(0)}(\dot{y},y,p,t)=0,$$

$$\left[\begin{array}{ccc}{T}_{\dot{y}}^{(0)}& T_y^{(0)}& T_z^{(0)}\\ F_{\dot{y}}^{(0)}& F_y^{(0)}& F_z^{(0)}\end{array}\right]\text{ is assumed to be nonsingular at }t=t_0,$$

$${\dot{y}}_i^{-}={\dot{y}}^{(i-1)}(p,t_i),y_i^{-}=y^{(i-1)}(p,t_i),z_i^{-}=z^{(i-1)}(p,t_i),$$

$$h^{(i-1)}\left({\dot{y}}_i^{-},y_i^{-},z_i^{-},p,t_i\right)=0,$$

$$T^{(i)}\left({\dot{y}}_i^{+},y_i^{+},z_i^{+},{\dot{y}}_i^{-},y_i^{-},z_i^{-},p,t_i\right)=0,$$

$$F^{(i)}\left({\dot{y}}_i^{+},y_i^{+},z_i^{+},p,t_i\right)=0,$$

$$\begin{array}{rl}& h_{\dot{y}}^{(i-1)}{\ddot{y}}_i^{-}+h_y^{(i-1)}{\dot{y}}_i^{-}+h_z^{(i-1)}{\dot{z}}_i^{-}+h_t^{(i-1)}\ne 0\\ \text{and}& \\ & \left[\begin{array}{ccc}{T}_{{\dot{y}}^{+}}^{(i)}& T_{y^{+}}^{(i)}& T_{z^{+}}^{(i)}\\ F_{\dot{y}}^{(i)}& F_y^{(i)}& F_z^{(i)}\end{array}\right]\text{ nonsingular at }t=t_i,\end{array}$$

$$F_{\dot{y}}\dot{s}+F_{y}s+F_{z}w+F_p=0,$$

$${\dot{s}}_i^{-}={\dot{s}}^{(i-1)}(p,t_i),s_i^{-}=s^{(i-1)}(p,t_i),w_i^{-}=w^{(i-1)}(p,t_i),$$

$$h_{\dot{y}}^{(i-1)}\left({\dot{s}}_i^{-}+{\ddot{y}}_i^{-}{\tau }_i\right)+h_y^{(i-1)}\left(s_i^{-}+{\dot{y}}_i^{-}{\tau }_i\right)+h_z^{(i-1)}\left(w_i^{-}+{\dot{z}}_i^{-}{\tau }_i\right)+h_t^{(i-1)}{\tau }_i+h_p^{(i-1)}=0,$$

$$\begin{array}{rl}& T_{{\dot{y}}^{+}}^{(i)}\left({\dot{s}}_i^{+}+{\ddot{y}}_i^{+}{\tau }_i\right)+T_{y^{+}}^{(i)}\left(s_i^{+}+{\dot{y}}_i^{+}{\tau }_i\right)+T_{z^{+}}^{(i)}\left(w_i^{+}+{\dot{z}}_i^{+}{\tau }_i\right)\\ & +T_{{\dot{y}}^{-}}^{(i)}\left({\dot{s}}_i^{-}+{\ddot{y}}_i^{-}{\tau }_i\right)+T_{y^{-}}^{(i)}\left(s_i^{-}+{\dot{y}}_i^{-}{\tau }_i\right)+T_{z^{-}}^{(i)}\left(w_i^{-}+{\dot{z}}_i^{-}{\tau }_i\right)\\ & +T_t^{(i)}{\tau }_i+T_p^{(i)}=0,\end{array}$$

$$F_{\dot{y}}^{(i)}\left({\dot{s}}_i^{+}+{\ddot{y}}_i^{+}{\tau }_i\right)+F_y^{(i)}\left(s_i^{+}+{\dot{y}}_i^{+}{\tau }_i\right)+F_z^{(i)}\left(w_i^{+}+{\dot{z}}_i^{+}{\tau }_i\right)+F_t^{(i)}{\tau }_i+F_p^{(i)}=0,$$

$$G(p)={\int }_{t_0}^{t_f}g(y,z,p,t)dt=\sum _{i=0}^N{\int }_{t_i}^{t_{i+1}}g(y,z,p,t)dt.$$

$$\frac{dG}{dp}=\sum _{i=0}^N\left[{\int }_{t_i}^{t_{i+1}}\left(g_{y}s+g_{z}w+g_p\right)dt+g\left(y_{i+1}^{-},z_{i+1}^{-},p,t_{i+1}\right){\tau }_{i+1}-g\left(y_i^{+},z_i^{+},p,t_i\right){\tau }_i\right],$$

$$\begin{array}{rl}\frac{dG}{dp}=\sum _{i=0}^N& [{\int }_{t_i}^{t_{i+1}}\left(g_{y}s+g_{z}w+g_p\right)dt+g_{i+1}^{-}{\tau }_{i+1}-g_i^{+}{\tau }_i\\ & +{\int }_{t_i}^{t_{i+1}}{\lambda }^T\left(F_{\dot{y}}\dot{s}+F_{y}s+F_{z}w+F_p\right)dt].\end{array}$$

$${\int }_{t_i}^{t_{i+1}}{\lambda }^T{F}_{\dot{y}}\dot{s}dt={\left({\lambda }^T{F}_{\dot{y}}s\right)|}_{t_i}^{t_{i+1}}-{\int }_{t_i}^{t_{i+1}}\frac{d}{dt}\left({\lambda }^T{F}_{\dot{y}}\right)sdt,$$

$$\begin{array}{rl}\frac{dG}{dp}=\sum _{i=0}^N& [{\int }_{t_i}^{t_{i+1}}\left(g_p+{\lambda }^T{F}_p\right)dt+g_{i+1}^{-}{\tau }_{i+1}-g_i^{+}{\tau }_i\\ & +{\left({\lambda }_{i+1}^{-}\right)}^T{F}_{\dot{y}}\left(t_{i+1}^{-}\right)s_{i+1}^{-}-{\left({\lambda }_i^{+}\right)}^T{F}_{\dot{y}}\left(t_i^{+}\right)s_i^{+}],\end{array}$$

$$\begin{array}{rl}& -\frac{d}{dt}\left({\lambda }^T{F}_{\dot{y}}\right)+{\lambda }^T{F}_y+g_y=0\\ & {\lambda }^T{F}_z+g_z=0.\end{array}$$

$$\begin{array}{rl}\frac{dG}{dp}& =\sum _{i=0}^N\left[{\int }_{t_i}^{t_{i+1}}\left(g_p+{\lambda }^T{F}_p\right)dt\right]\\ & +\sum _{i=0}^N\left[g_{i+1}^{-}{\tau }_{i+1}-g_i^{+}{\tau }_i+{\left({\lambda }_{i+1}^{-}\right)}^T{F}_{\dot{y}}\left(t_{i+1}^{-}\right)s_{i+1}^{-}-{\left({\lambda }_i^{+}\right)}^T{F}_{\dot{y}}\left(t_i^{+}\right)s_i^{+}\right]\\ & =\sum _{i=0}^N\left[{\int }_{t_i}^{t_{i+1}}\left(g_p+{\lambda }^T{F}_p\right)dt\right]\\ & +{\lambda }_f^T{F}_{\dot{y}}\left(t_f\right)s_f-{\lambda }_0^T{F}_{\dot{y}}\left(t_0\right)s_0\\ & +\sum _{i=1}^N\left[\left(g_i^{-}-g_i^{+}\right){\tau }_i+{\left({\lambda }_i^{-}\right)}^T{F}_{\dot{y}}\left(t_i^{-}\right)s_i^{-}-{\left({\lambda }_i^{+}\right)}^T{F}_{\dot{y}}\left(t_i^{+}\right)s_i^{+}\right],\end{array}$$

$${\lambda }_f^T{F}_{\dot{y}}=0.$$

$$\left[{\left({\lambda }_i^{-}\right)}^T{F}_{\dot{y}}+\left(g_i^{-}-g_i^{+}\right){\beta }_i\right]s_i^{-}={\left({\lambda }_i^{+}\right)}^T{F}_{\dot{y}}\left(t_i^{+}\right)s_i^{+},$$

$$\dot{y}$$

$$0$$

$$\left[\begin{array}{c}{T}_y^{(0)}\\ C^{(0)}\end{array}\right]\text{ nonsingular at }t=t_0.$$

$$y_i^{-}=y^{(i-1)}(p,t_i),z_i^{-}=z^{(i-1)}(p,t_i),$$

$$h^{(i-1)}\left(y_i^{-},z_i^{-},p,t_i\right)=0,$$

$$T^{(i)}\left(y_i^{+},y_i^{-},p,t_i\right)=0,$$

$${\dot{y}}_i^{+}=f^{(i)}\left(y_i^{+},z_i^{+},p,t_i\right),$$

$$0=k^{(i)}(y_i^{+},p,t_i),$$

$$0=C^{(i)}{\dot{y}}_i^{+}+k_t^{(i)}.$$

$$\begin{array}{rl}& h_y^{(i-1)}{\dot{y}}_i^{-}+h_z^{(i-1)}{\dot{z}}_i^{-}+h_t^{(i-1)}\ne 0\\ \text{and}& \\ & \left[\begin{array}{c}{T}_{y^{+}}^{(i)}\\ C^{(i)}\end{array}\right]\text{ nonsingular at }t=t_i,\end{array}$$

$$\dot{s}$$

$$0$$

$$\begin{array}{rl}\frac{dG}{dp}& =\sum _{i=0}^N\left[{\int }_{t_i}^{t_{i+1}}\left(g_p+{\lambda }^T{f}_p+{\mu }^T{k}_p\right)dt\right]\\ & -{\lambda }_f^T{s}_f+{\lambda }_0^T{s}_0\\ & +\sum _{i=1}^N\left[\left(g_i^{-}-g_i^{+}\right){\tau }_i-{\left({\lambda }_i^{-}\right)}^T{s}_i^{-}+{\left({\lambda }_i^{+}\right)}^T{s}_i^{+}\right],\end{array}$$

$${\dot{\lambda }}^T+{\lambda }^T{f}_y+{\mu }^T{k}_y+g_y$$

$${\lambda }^T{f}_z+g_z$$

Full text

Keywords: Forward Sensitivity Analysis (FSA), Adjoint Sensitivity Analysis (ASA), Hybrid discrete–continuous systems, Hybrid systems with memory

1. 1 Introduction

2. 2 Problem Formulation

3. 3 Sensitivity Analysis for Hybrid Systems

4. 3.1 Sensitivity Equations (FSA)

5. 3.2 Adjoint Sensitivity Analysis (ASA)

6. 3.3 Hessenberg index-2 problems

7. 4 Sensitivity Analysis for Hybrid Systems with Memory

8. 4.1 Sensitivity Equations (FSA)

9. 4.2 Adjoint Sensitivity Analysis (ASA)

10. 5 Numerical Examples

11. 5.1 A Simple Hybrid ODE

12. 5.2 A Computational Mechanics Problem

13. 6 Conclusions

1 Introduction

Sensitivity analysis, the study of how variation in a model output can be apportioned to different sources of variation, has a wide range of applications in science and engineering. These include model evaluation (finding the most and least influential parameters), generation of reduced-order models (reducing model complexity while preserving the input-output behavior), uncertainty quantification (characterizing and reducing uncertainty in model predictions), data assimilation (merging observations into a model to improve its accuracy), and dynamic optimization (finding model parameters for improved system response).

There is a significant body of work on methods and applications of sensitivity analysis, spanning at least the last three decades. For dynamic models described by systems of ordinary differential equations (ODEs) or differential-algebraic equations (DAEs), like those considered herein, two main approaches for continuous 1st order sensitivity analysis have been proposed and investigated. The forward sensitivity method (FSA) [1, 2, 3], which uses linearizations of the original model to calculate the state sensitivities, is relatively simple to formulate and implement, but can quickly become prohibitively expensive for large systems and large numbers of model parameters. Adjoint sensitivity analysis (ASA) provides an attractive alternative for problems with an arbitrary number of models parameters but relatively few output functionals for which gradients are sought. More complicated, both from a formulation and implementation point of view, ASA requires the formulation of adjoint models [4] and, in the case of time-dependent problems, backward in time integration of a system that may depend on the forward solution. The adjoint models and their properties have been studied for systems of ODEs, index-1 and index-2 DAEs [5, 6], as well as for the index-3 DAEs in multibody dynamics [7]. General-purpose software for FSA and ASA of systems described by ODEs and DAEs is available and widely used [8, 9, 10].

While many physical phenomena can be described with continuous differential models, there are applications that require models for interaction between continuous and discrete phenomena. Such hybrid discrete – continuous systems, with the continuous part described by systems of ODEs or DAEs and the discrete part modeled by finite automata, arise in control systems (e.g., safety interlock systems), biology (cell signaling models), chemical engineering (plant processing models), electrical engineering (heterogeneous electrical grid models), as well as mechanical engineering. In the area of computational mechanics, two applications stand out: collision and contact problems formulated in an event-driven framework; and rate-independent hysteretic phenomena in constitutive equations. Forward sensitivity analysis for hybrid systems of this type was treated in a comprehensive manner by Barton and co-workers [11, 12], who highlighted the problems associated with hidden discontinuities during sensitivity analysis and the importance of properly computing and incorporating potential sensitivity jumps. Adjoint methods for hybrid discrete – continuous problems have seen relatively less attention. Recent work by Corner et.al. [13] treats ASA for hybrid multibody dynamical systems in the context of non-smooth contact dynamics. In [14], the discrete adjoint method (adjoint of the discretized time-dependent problem) was used for computing derivative information for optimization of power system dynamics.

Here, we present an adjoint sensitivity method for hybrid discrete – continuous systems, extending the forward sensitivity methods presented in [11, 12]. We treat ordinary differential equations (ODEs) and DAEs of index up to two (Hessenberg) and provide sufficient solvability conditions for consistent initialization and state transfer at mode switching points, for both the sensitivity and adjoint systems. Furthermore, we extend the analysis to so-called hybrid systems with memory where the dynamics of any given mode depend explicitly on the states at the last mode transition point.

This paper is structured as follows. In §2 we provide the framework for the type of hybrid systems considered herein. In §3 we derive the sensitivity and adjoint systems for hybrid parameter-dependent problems whose dynamics are specified as index-0 (ODEs) or index-1 DAEs. The case of Hessenberg index-2 DAEs is treated separately in §3.3. In §4 we discuss the sensitivity formulation for systems with memory, also providing an algorithm for propagating backward in time the adjoint systems while properly accounting for jumps at transition points. Numerical examples are discussed in §5, including a computational mechanics problem based on the exponential model constitutive material law for steel structures under cyclic loading.

2 Problem Formulation

Let

[TABLE]

be a parameter-dependent DAE system, where $x=\{y,z\}\in\mathbb{R}^{N_{x}}=\mathbb{R}^{N_{y}+N_{z}}$, with $y$ and $\dot{y}$ being the differential states and their time derivatives, $z$ the algebraic variables, and $p\in\mathbb{R}^{N_{p}}$ a set of time-invariant problem parameters.

In this paper we consider hybrid discrete – continuous systems which can be described following the formalism in [11, 12]. Such systems can be in one of several modes $\mathbb{S}_{i}^{i+1}$, $i=0,\ldots,N$ where each mode and the transitions between modes are characterized by (see also Fig. 1):

Variables. The state of the system in mode $\mathbb{S}_{i}^{i+1}$ is described by a set of variables $\left\{\dot{y}^{(i)},y^{(i)},z^{(i)}\right\}$, functions of time $t$ and of the time-invariant parameters $p$;

Dynamics. In mode $\mathbb{S}_{i}^{i+1}$, the dynamics of the system are given by DAEs of the form of Eq. 1. In this paper, we consider problems for which $\text{rank}\left([F_{\dot{y}}|F_{z}]\right)=N_{y}+N_{z}$ which is true for index-0 (ODEs) and most index-1 DAEs. The case of different DAE structures (namely Hessenberg index-2) is treated separately in §3.3.

Transitions. The system transitions to a different mode based on the so-called transition conditions $h^{(i-1)}\left(\dot{y}^{(i-1)},y^{(i-1)},z^{(i-1)},p,t\right)$. These conditions implicitly define the transition time $t_{i},i=1,\ldots,N$ at which the system switches to mode $\mathbb{S}_{i}^{i+1}$.

The transition functions $T^{(i)}\left(\dot{y}^{(i)},y^{(i)},z^{(i)},\dot{y}^{(i-1)},y^{(i-1)},z^{(i-1)},p,t\right)$ define a system of equations that map the final states in mode $\mathbb{S}_{i-1}^{i}$ (at the transition time $t_{i}$) to the initial states in the next mode $\mathbb{S}_{i}^{i+1}$ (after the transition time).

Without loss of generality, we assume that the initial and final simulation times do not depend on the problem parameters $p$. Initial conditions are provided as:

[TABLE]

which are assumed to provide $N_{y}$ additional conditions consistent with the DAE of Eq. 1. In other words,

[TABLE]

where $F^{(0)}$ defines the system dynamics on the first mode $\mathbb{S}_{0}^{1}$.

With the above description, we assume that the system undergoes $N$ transitions at times $t_{i},i=1,\ldots,N$ and evolves through a sequence of $N+1$ modes $\mathbb{S}_{i}^{i+1}$, each over $t\in[t_{i},t_{i+1}]$. At each transition time $t_{i}$, we denote by $a^{-}_{i}$ a quantity $a$ evaluated before the transition; e.g., $y^{-}_{i}=y(t^{-}_{i})=y(t_{i}-0)$ and by $a^{+}_{i}$ the same quantity after the transition (possibly including a discontinuity); e.g, $y^{+}_{i}=y(t^{+}_{i})=y(t_{i}+0)$.

The transition times are defined implicitly by zero-crossings of the transition conditions $h^{(i-1)}$ as solution of the equations $h^{(i-1)}\left(\dot{y}^{-}_{i},y^{-}_{i},z^{-}_{i},p,t_{i}\right)=0$ and the initial values for the DAE problem over the next mode are dictated by the system of equations $T^{(i)}\left(\dot{y}^{+}_{i},y^{+}_{i},z^{+}_{i},\dot{y}^{-}_{i},y^{-}_{i},z^{-}_{i},p,t_{i}\right)=0$.

Solvability conditions of the resulting hybrid problem at each transition point stem from imposing, at each transition time $t_{i}$, the following nonlinear equations:

[TABLE]

where $F^{(i)}(\dot{y},y,z,p,t)$ defines the system dynamics over the next mode, $\mathbb{S}_{i}^{i+1}$. Sufficient conditions for the solvability of this problem (based on the implicit function theorem) are

[TABLE]

with the second simply being the conditions for consistent DAE initialization at $t_{i}$. Note that all partial derivatives above must be evaluated at $t_{i}$.

3 Sensitivity Analysis for Hybrid Systems

Forward sensitivity analysis and the general equations for parameteric sensitivity functions for hybrid systems like those described in §2 have been presented in  [11, 12]. We briefly review these results in §3.1 insofar as they are a required step in deriving the adjoint equations and performing adjoint sensitivity analysis for such systems, as presented in §3.2.

3.1 Sensitivity Equations (FSA)

Differentiating the system DAE of Eq. 1 with respect to the independent variables $p$, the sensitivity trajectories in a given mode $\mathbb{S}_{i}^{i+1}$ can be obtained by integrating the linear time varying sensitivity DAE simultaneously with the original DAE:

[TABLE]

where $s$ represents the $N_{y}\times N_{p}$ matrix of differential state sensitivities and $w$ is the $N_{z}\times N_{p}$ matrix of algebraic state sensitivities.

At each transition point $t_{i}$, in addition to computing consistent initial values for the DAE states in the new mode, the sensitivity system in the new mode must also be consistently initialized – often requiring calculating explicit jumps in the sensitivities. This procedure requires first calculating the sensitivity of the transition times themselves which can be obtained from a formal differentiation of the transition conditions of Eq. 3 to obtain:

[TABLE]

where $\tau_{i}={dt_{i}}/{dp}$ is the sensitivity of the transition time $t_{i}$ with respect to the problem parameters $p$, $s^{(i-1)}(p,t)$ and $w^{(i-1)}(p,t)$ are the state sensitivity trajectories over the previous mode, and all partial derivatives are evaluated at $t_{i}$. It can be easily seen that sufficient conditions for the solvability of the above system of equations are the same as those of Eq. 4.

Note that $\ddot{y}$ and $\dot{z}$ in Eq. 6 can be obtained by differentiating Eq. 1 with respect to the independent variable $t$: $F_{\dot{y}}\ddot{y}+F_{y}\dot{y}+F_{z}\dot{z}+F_{t}=0$, which is solvable under the assumption that $[F_{\dot{y}}|F_{z}]$ has full rank.

With the above, we can propagate forward in time the augmented DAE and sensitivity system through the sequence of modes $\mathbb{S}_{i}^{i+1},i=0,\ldots,N$, properly taking into account any required discrete updates (jumps) at the transition points $t_{i},i=1,N$.

3.2 Adjoint Sensitivity Analysis (ASA)

The augmented DAE and sensitvity system of Eqs. 1 and 5 can be quite large, especially for large problems. Even though there are efficient algorithms that exploit the structure of this augmented system [2, 3, 10], the cost of such computations can become prohibitive when there are many problem parameters. The adjoint sensitivity method has proven to be a viable, effective alternative, at least when interested in the parameter sensitivity of one or only a few functionals that depend on the system state [4, 5, 6, 10].

In this section, we present the derivation of the adjoints for hybrid systems and provide the necessary transition conditions for the adjoint variables, required for correct ASA of such problems. As before, we consider here index-0 (i.e., ODEs) and index-1 problems; the case of the Hessenberg index-2 DAE structure is treated in §3.3.

Consider a functional of the system state trajectory of the form:

[TABLE]

Its gradient with respect to the time-invariant problem parameters is then

[TABLE]

where we have used the Leibniz integral rule to account for the dependency of the transition times on $p$ (recall that $\tau_{i}={dt_{i}}/{dp}$). Introducing adjoint variables $\lambda\in\mathbb{R}^{N_{y}+N_{z}}$ and taking into consideration Eq. 5,

[TABLE]

Using integration by parts

[TABLE]

and rearranging terms, we obtain

[TABLE]

if the adjoint variables are selected to satisfy

[TABLE]

The computation of the gradient ${dG}/{dp}$ using the adjoint method is then completed if we impose appropriate jump conditions on $\lambda$ at the transition times $t_{i}$ so that the sensitivities $s_{i}$ cancel out in Eq. 11. To obtain these conditions, we start by rearranging terms and favorably changing the summation indices:

[TABLE]

where we used the fact that the initial ($t_{0}$) and final ($t_{f}=t_{N+1}$) times do not depend on the problem parameters $p$. For index-0 and index-1 DAEs, the initial conditions (at $t=t_{f}$) for the adjoint variables can be simply

[TABLE]

As shown in §3.3 this choice will not suffice for a Hessenberg index-2 system.

Finally, the jump conditions for the adjoint variables at each transition time are inferred from the components of the last term in Eq. 13 by imposing

[TABLE]

where we formally expressed $\tau_{i}=\alpha_{i}+\beta_{i}s^{-}_{i}$ from Eq. 6. These relations, together with the algebraic constraint in Eq. 12, allow us to calculate $\lambda^{-}_{i}$, the adjoint variables after crossing a transition point, in terms of $\lambda^{+}_{i}$, the adjoint variables at the end of (backward) integration over the previous mode. Note that the assumption $[F_{\dot{y}}|F_{z}]$ nonsingular again provides sufficient solvability conditions.

3.3 Hessenberg index-2 problems

The assumptions used previously for establishing the solvability conditions at transition points for the original DAE (Eq. 1), the linear sensitivity system (Eq. 5), and the adjoint problem (Eq. 12) will not suffice for a Hessenberg index-2 DAE system. We provide these derivations in this section. The main reason for including this type of problems in our analysis is that this is the form of DAEs obtained when modeling incompressible flow by discretized Navier-Stokes equations in CFD [15] or with a stabilized index reduction of the index-3 equations of motion in multibody system problems [16, 17].

Let the dynamics of a hybrid system such as those defined in §2 be described, in each mode $\mathbb{S}_{i}^{i+1}$, by DAEs in Hessenberg index-2 form:

[TABLE]

where $f^{(i)}$ is the right-hand side of the differential equations, $k^{(i)}$ denotes the algebraic equations, and $C^{(i)}B^{(i)}$ is nonsingular, with $A^{(i)}=f^{(i)}_{y}$, $B^{(i)}=f^{(i)}_{z}$, and $C^{(i)}=k^{(i)}_{y}$.

Obtaining consistent initial conditions for higher-index DAEs is more complicated and no general recipe exists. It is clear that not any choice of $N_{y}$ initial conditions for the differential variables will do, as $y_{0}=y(t_{0})$ must satisfy the algebraic equations 15b. Furthermore, higher-index DAEs have so-called hidden manifolds to which the solution is restricted; for the system above, additional constraints are given by the differentiated algebraic equations: $Cf+k_{t}=0$. See [18] for further details on consistent initialization of DAE and section §4.2 in [9] for a discussion of practical implementation considerations. Here, we assume that initial conditions $y_{0}$ consistent with Eq. 15 are provided. Furthermore, we assume these are given implicitly as $0=T^{(0)}(y,p,t)\in\mathbb{R}^{N_{y}-N_{z}}$. A sufficient condition (based on the implicit function theorem) for the solvability of the resulting nonlinear system is then

[TABLE]

Initial conditions for the algebraic variables are then obtained by solving the nonlinear system $C^{(0)}f^{(0)}+k^{(0)}_{t}=0$, evaluated at $\{t_{0},y_{0}\}$. Note that the regularity condition ($CB$ nonsingular) is critical in providing solvability conditions for this step. Finally, set $\dot{y}_{0}=f^{(0)}(y_{0},z_{0},p,t_{0})$.

As before, consider a hybrid discrete – continuous system whose dynamics in mode $\mathbb{S}_{i}^{i+1}$ are given by equations of the form 15, with transitions111For simplicity, we ignore here more complex transition conditions and functions which, in general, may also involve the time derivatives and/or algebraic variables. The case treated here is typical for problems in CFD or multibody dynamics. defined by the transition conditions $h^{(i)}\left(y^{(i)},z^{(i)},p,t\right)$ and the associated transition functions $T^{(i)}\left(y^{(i)},y^{(i-1)},p,t\right)$.

The state transition conditions of the resulting hybrid problem become

[TABLE]

Sufficient conditions for computing the state transitions are:

[TABLE]

in addition to the regularity condition on $\mathbb{S}_{i}^{i+1}$, namely $\exists\left(C^{(i)}B^{(i)}\right)^{-1}$. The transition time $t_{i}$ is solved from Eq. 16b; the differential states $y^{+}_{i}$ are obtained by simultaneously solving Eqs. 16c and 16e; finally, the algebraic states after the transition are obtained from the nonlinear system of Eq. 16f and the state derivatives $\dot{y}^{+}_{i}$ calculated with Eq. 16d.

The corresponding sensitivity equations are 222For clarity, we drop all subscripts and superscripts in the remainder of this section, except where necessary.

[TABLE]

The sensitivity transition conditions are obtained in a manner completely analogous to Eq. 6 and are not reproduced here. Similar to the case of index-0 and index-1 systems, the solvability conditions for computing the sensitivity of the transition times and the state sensitivity jump transitions are those of Eq. 17. Furthermore, we can formally write $\tau_{i}=\alpha_{i}+\beta_{i}s^{-}_{i}$ and $s^{+}_{i}=\Gamma_{i}+\Delta_{i}s^{-}_{i}$.

Considering the functional $G$ of Eq. 7 and following a derivation similar to that of §3.2, its gradient with respect to the problem parameters $p$ is:

[TABLE]

where the adjoint variables $\lambda$ and $\mu$ were chosen to satisfy

[TABLE]

in each interval $[t_{i},t_{i+1}]$. Because of the structure of the DAE, the choice $\lambda_{f}=0$ as initial conditions for the adjoint hybrid system will not suffice, as this could conflict with Eq. 20b if $g(y,z,p,t)$ depends explicitly on $z$. Instead, with all quantities evaluated at $t=t_{f}$ using functions corresponding to mode $\mathbb{S}_{N}^{N+1}$, we require that

[TABLE]

for some $\xi$ yet to be determined [6]. Inserting into Eq. 20b, we have

[TABLE]

Therefore,

[TABLE]

and since $Cs_{f}=-k_{p}$ (from Eq. 18b)

[TABLE]

Consistent initial conditions for $\mu_{f}$ are obtained from the hidden constraint resulting from differentiation of the algebraic equations 20b. Indeed,

[TABLE]

and inserting $\dot{\lambda}$ from Eq. 20a

[TABLE]

A similar procedure can be applied to obtain the transfer conditions (jumps) of the adjoint variables at each transition point. Indeed, we are seeking $\lambda^{-}_{i}$ such that each term in the last sum on the right-hand side of Eq. 19 can be evaluated without the need for computing the state sensitivities $s^{-}_{i}$ and $s^{+}_{i}$. Using $\tau_{i}=\alpha_{i}+\beta_{i}s^{-}_{i}$ and $s^{+}_{i}=\Gamma_{i}+\Delta_{i}s^{-}_{i}$, after simple algebraic manipulations, we have

[TABLE]

where

[TABLE]

Note that we cannot simply impose $\lambda^{-}_{i}=X_{i}$ as this could conflict with the algebraic equations in Eq. 20b. Instead, using quantities evaluated at $t^{-}_{i}$, we seek a vector $\zeta_{i}$ such that

[TABLE]

Inserting into Eq. 20b (evaluated at $t^{-}_{i}$ using the DAE for mode $\mathbb{S}_{i-1}^{i}$), we obtain

[TABLE]

Therefore,

[TABLE]

and since $Cs^{-}_{i}=-k_{p}$ (from Eq. 18b)

[TABLE]

The gradient of Eq. 19 then becomes

[TABLE]

4 Sensitivity Analysis for Hybrid Systems with Memory

In some situations, the system dynamics in mode $\mathbb{S}_{i}^{i+1}$ depends explicitly on the states at the transition into mode $\mathbb{S}_{i}^{i+1}$. Such a problem is the computational mechanics example discussed in §5.2 which uses the so-called exponential model [19] model for steel reinforcement bars under cyclic loading.

In such cases, this additional explicit dependency must be properly taken into account in deriving the sensitivity and adjoint systems and the corresponding transition rules. For brevity in exposition, we consider here only hybrid continuous-discrete systems whose dynamics between any two transition points are described by DAEs in Hessenberg index-1 form:

[TABLE]

over the interval $[t_{0},t_{f}]$, with $k_{z}$ nonsingular and $p\in\mathbb{R}^{N_{p}}$ a set of time-independent problem parameters. The differential variables are $y\in\mathbb{R}^{N_{y}}$, $z\in\mathbb{R}^{N_{z}}$ are the algebraic variables, and without loss of generality we assume that $f$ and $k$ do not explicitly depend on time. Furthermore, we assume that consistent initial conditions are provided and, for simplicity, we assume that $y_{0}\equiv y(t_{0})$ does not depend on $p$. Let the transition times $t_{i}$ between consecutive modes be implicitly defined by the transition conditions333Without loss of generality, we assume that the transition conditions $h^{(i-1)}$ do not depend explicitly on the algebraic variables, as those can always be formally expressed as $z=z(y)$ from the algebraic equations 22b.

[TABLE]

Let $t_{f}=t_{N+1}$ and let $\mathbb{S}_{i}^{i+1}$ be the mode over the interval $[t_{i},t_{i+1}]$. The transition into $\mathbb{S}_{i}^{i+1}$ at $t_{i}$ is defined by the transition functions

[TABLE]

which, in general, may also include explicit dependencies on time or the time derivatives ${\dot{y}}^{-}_{i}$ and ${\dot{y}}^{+}_{i}$, before and after the transition, respectively. Equation 24 implicitly defines $y^{+}_{i}$ and $z^{+}_{i}$, the system states after the transition. While in mode $\mathbb{S}_{i}^{i+1}$, the equations of motion are those of Eq. 22 with $y^{*}=y^{+}_{i}\equiv y(t^{+}_{i})$ and $z^{*}=z^{+}_{i}\equiv z(t^{+}_{i})$. Solvability of the hybrid problem requires that the transition conditions of Eq. 24 be consistent with the algebraic equations of Eq. 22b at time $t^{+}_{i}$. A sufficient condition is:

[TABLE]

In many physical phenomena, including the example problem described in §5.2, it is common for the system states to be $\mathcal{C}^{0}$ continuous at transition points; i.e., $y^{+}_{i}=y^{-}_{i}$ and $z^{+}_{i}=z^{-}_{i}$.

We call $y^{*}$ and $z^{*}$ memory states and the problem of Eq. 22 a system with memory. In this section, we extend the previous formulation to perform (forward and adjoint) sensitivity analysis of such systems. The new elements in this derivation are related to the implicit dependency of the memory states on problem parameters.

4.1 Sensitivity Equations (FSA)

Denoting by $s\equiv{dy}/{dp}$ and $w\equiv{dz}/{dp}$ the state sensitivities, a formal differentiation of Eq. 22 leads to the following forward sensitivity equations:

[TABLE]

where, in mode $\mathbb{S}_{i}^{i+1}$, we have $s^{*}=s^{+}_{i}\equiv s(t^{+}_{i})$ and $w^{*}=w^{+}_{i}\equiv w(t^{+}_{i})$.

Consistent initialization of the FSA problem after a transition can be carried out as follows. Consider the transition at $t_{i}$ and, for clarity, drop the superscripts $(i)$ and subscripts $i$. First, formal differentiation of the transition condition implies

[TABLE]

where we include the implicit dependency on $p$ through the time derivative $\dot{y}$. Recall that $\tau=dt_{i}/dp$. This determines the sensitivity of the transition time as:

[TABLE]

where all quantities are evaluated at $t^{-}_{i}$.

Next, a formal differentiation of the transition functions of Eq. 24 gives

[TABLE]

which, together with the consistency conditions of Eq. 25b, determines the sensitivity states after the transition, $s^{+}$ and $w^{+}$, given $s^{-}$ and $w^{-}$. Note that the time derivatives of the algebraic states, $\dot{z}$, can be obtained simultaneously with $\ddot{y}$ by taking the time derivative of Eq. 22 and solving, at time $t^{+}_{i}$, the resulting linear system with a nonsingular matrix

[TABLE]

(based on the Hessenberg index-1 condition). In other words, $s^{+}$, $w^{+}$, and $\dot{s}^{+}$ are solved for simultaneously considering Eqs. 28 and 25; this problem is well posed given the solvability conditions for the transition functions of the original problem.

Integration of the forward sensitivity equations can then be carried out simultaneously with the original problem of Eq. 22, including the necessary jumps in sensitivity variables at each mode transition time.

4.2 Adjoint Sensitivity Analysis (ASA)

Consider next a functional of the form

[TABLE]

and the problem of calculating the gradient of $G$ with respect to the problem parameters $p$.

Differentiating Eq. 29 and using the Leibniz integral rule to take into account the dependency of the transition times $t_{i}$ on the problem parameters $p$, we have

[TABLE]

where $g^{-}_{i+1}\equiv g(y^{-}_{i+1},z^{-}_{i+1},p,t_{i+1})$ and $g^{+}_{i}\equiv g(y^{+}_{i},z^{+}_{i},p,t_{i})$.

Introducing adjoint variables $\lambda$ and $\mu$ (of appropriate dimensions) and using Eq. 25, we have:

[TABLE]

where we have taken into account that, in mode $\mathbb{S}_{i}^{i+1}$, $s^{*}=s^{+}_{i}$ and $w^{*}=w^{+}_{i}$. Next, use integration by parts and rearrange terms to obtain

[TABLE]

where we define

[TABLE]

As before, the goal is to impose conditions on the adjoint variables $\lambda$ and $\mu$ so that the gradient of $G$ can be computed without the need to compute forward state sensitivities. We first impose that

[TABLE]

in each interval $[t_{i},t_{i+1}]$. Then, changing summation indices as appropriate and using the fact that the integration limits $t_{0}$ and $t_{f}$ are constant (i.e., $\tau_{0}=\tau_{N+1}=0$),

[TABLE]

The last term in the above equation dictates the choice of transfer conditions for the adjoint variable $\lambda$ such that all terms involving state sensitivities cancel out. To do this, we use the conditions of Eq. 28 in conjunction with the algebraic equations of Eq. 25b evaluated to the left and to the right of the transition at $t_{i}$:

[TABLE]

Also taking into account the dependency of $\tau_{i}$ on the sensitivities $s^{-}_{i}$ (Eq. 27), we can formally find a recursion formula for $w^{+}_{i}$ in terms of the sensitivities at all transition times to the left of $t_{i}$.

Setting $\lambda(t_{f})\equiv\lambda(t^{-}_{N+1})=0$ and computing $\mu(t_{f})$ to be consistent with Eq. 36b, we propagate backward in time the adjoint DAE of Eq. 36 until the transition time $t^{+}_{N}$. At this point, the adjoint variables $\lambda^{-}_{N}$, past the transition, are selected as to cancel all terms involving $s^{+}_{N}$, $s^{-}_{N}$, $w^{+}_{N}$, and $w^{-}_{N}$ in the last term of Eq. 37. This process is repeated in turn for all modes $\mathbb{S}_{i}^{i+1}$ with $i=N,\ldots,0$.

To illustrate this process, consider the following simplifying assumptions:

• •

State and state sensitivities are $\mathcal{C}^{0}$ continuous at all transitions:

[TABLE]

Note that this also implies that

[TABLE]

• •

Initial conditions do not depend on the problem parameters:

[TABLE]

• •

Without loss of generality (given that $k_{z}$ is nonsingular), the algebraic equations 25b can be written as:

[TABLE]

Under these assumptions, the gradient of $G$ becomes

[TABLE]

The transition conditions for adjoint variables, moving backward in time, are then

[TABLE]

where all partial derivatives are evaluated at the end of the time interval for the corresponding mode; e.g. $k^{(N)}_{y}\equiv k^{(N)}_{y}(t^{-}_{N})$. Note that, even though both the states and sensitivities are continuous across transitions, the adjoint variables are discontinuous with discrete jumps (see also the example in §5.2). The complete procedure, also including the required updates to $dG/dp$ at each transition point, is listed in Algorithm 1.

5 Numerical Examples

5.1 A Simple Hybrid ODE

The first example is reproduced from [11, 12] and is a single ODE which transitions between two modes:

[TABLE]

with initial condition $x(t_{0})=0$ and a nominal value $p=2.9$. The transition function for this example is simply $x^{+}=x^{-}$.

This example is used here to compare the solution and (forward) sensitivity results with those presented in [11, 12] and then further compare results of FSA and ASA in computing the derivative of $G=\int_{t_{0}}^{t_{f}}xdt$ with respect to the parameter $p$.

Figure 2 shows the state variable $x$ and its sensitivity $s=dx/dp$ over the interval $t\in[0,5]$. These results clearly indicate that explicitly accounting for the discrete jumps in the state sensitivity at the transition times is crucial in obtaining the correct result444As shown in [12], ignoring the hidden discontinuities results in $s\equiv 0$ at all times. The sensitivity discontinuities can be obtained from the differentiated transition function evaluated at a transition time $t^{*}$:

[TABLE]

where the sensitivity of the transition time is obtained from the transition condition $h\equiv x^{3}-5x^{2}+7x-p$:

[TABLE]

where all quantities are evaluated right before the transition at $t^{*}$.

With the state sensitivities such computed, we obtain $\frac{d}{dp}\int_{0}^{5}xdt=-2.31195$.

Using the adjoint method, the derivative of $G$ is obtained as:

[TABLE]

where the adjoint equations are

[TABLE]

with final conditions $\lambda(t_{f})=0$ and transition relationships

[TABLE]

The trajectory of the adjoint variable is shown in Fig. 3. Note that, for the purpose of evaluating $\frac{dG}{dp}$ backward integration for $\lambda$ could have been stopped once the first transition time $t_{1}$ is reached. The desired derivative is obtained as $\frac{d}{dp}\int_{0}^{5}xdt=-2.31195$.

5.2 A Computational Mechanics Problem

This subsection develops a one degree of freedom model using a general constitutive equation with application to computational mechanics. The internal forces are adapted from the exponential model (EM) of Vaiana et al. [19]. This material description is intended to fit the rate-independent hysteretic behavior of materials. The EM does not require the use of return-mapping algorithms and can describe a wide range of behaviors.

The dynamics of a single degree of freedom oscillator (see Fig. 4) can be cast as a Hessenberg index-1 hybrid system with memory:

[TABLE]

where $u$ is the deformation displacement, $z$ is the internal stress, $A$ is the cross section area, and $m$ is the mass of the system, concentrated in one node. The stress is determined by [19]

[TABLE]

and depends on the material parameters $p=[k_{a},k_{b},\alpha,\beta]$. Transitions are controlled by zero-crossings of the velocity $v$ (since $\xi=sign(v)$) and the dynamics in each mode depend on the state at the last reversal $i$ through the derived memory variable $u^{*}_{i}$:

[TABLE]

Here, the quantities $u_{0}$ and $\bar{f}$ are (constant) internal model parameters

[TABLE]

with $\delta=10^{-20}$.

At each transition point, the state is continuous; i.e., the transition functions are simply $u^{+}=u^{-}$ and $v^{+}=v^{-}$. In each mode, the solution depends on the state at the last transition. Therefore, consistency at each transition point implies that $z^{+}=z^{-}$ and $\dot{u}^{+}=\dot{u}^{-}$, $\dot{v}^{+}=\dot{v}^{-}$.

Starting with consistent initial conditions $u(t_{0})=v(t_{0})=z(t_{0})=0$, with a sinusoidal load $f(t)=0.5t\sin{2\pi t}$ applied to the free end of the bar, the system undergoes multiple reversals over a time interval $t\in[0,10]$ as illustrated in Figs. 5 and 6. The nominal values for the model parameters were $p=[k_{a},k_{b},\alpha,\beta]=[32\pi^{2},\pi^{2},205,0]$, with the mass of the node receiving the excitation $m=1$ and the cross section area $A=1$.

Following the derivation in §4, the sensitivity equations for this problem take the form

[TABLE]

with $s_{u}=\partial u/\partial p\in\mathbb{R}^{1\times 4}$, $s_{v}=\partial v/\partial p\in\mathbb{R}^{1\times 4}$, and $w=\partial z/\partial p\in\mathbb{R}^{1\times 4}$. The partial derivatives of the stress function $\sigma$ in the sensitivity algebraic equation were obtained analytically. It is easy to see that, at each transition time $t_{i}$, the sensitivities of the differential variables are $\mathcal{C}^{1}$ continuous, while the sensitivity of the algebraic variable is $\mathcal{C}^{0}$ continuous. See Fig. 7.

As output functional, we consider a measure of total deformation over the time interval considered defined as $G(p)=\int_{t_{0}}^{t_{f}}u^{2}dt$. During integration of the DAE and sensitivity equations, additional quadratures are carried out to compute $G$ and its gradient with respect to the problem parameters $p$. With the current choice of nominal model parameter values and a final simulation time $t_{f}=10$, the computed value is $G=0.04994$.

Finally, the procedure of §4.2 described in Algorithm 1 was employed to compute the gradient $dG/dp$ with the ASA approach. The adjoint systems for this problem assume the form:

[TABLE]

For this example problem, the recurrence relationship controlling the transfer of adjoint variables at transition points has the form:

[TABLE]

Starting with the consistent conditions $\lambda(t_{f})=[0\,,\,0]$ and $\mu(t_{f})=0$ at $t_{f}=10$, the adjoint solution is shown in Fig. 8. It is important to highlight that, even though both the states and sensitivity trajectories in this problem are $\mathcal{C}^{0}$ continuous across transitions, the adjoint variables are discontinuous with discrete jumps which must be properly accounted for and incorporated during the backward-in-time propagation of the associated hybrid adjoint system.

As an additional check on the output functional gradient computed using FSA and ASA approaches, we also provide below an estimate based on (forward) finite differences, using a relative perturbation factor $\bar{\epsilon}=10^{-4}$. The gradient of $dG/dp$ obtained with different techniques is:

[TABLE]

consistent with the integration tolerances and the relative finite difference perturbation factor. We note that integration of all DAEs (both forward and backward in time) was conducted with the ode15s Matlab solver, using a relative tolerance of $10^{-8}$ and an absolute tolerance of $10^{-12}$.

6 Conclusions

We have derived the sensitivity equations for forward sensitivity analysis and the adjoint equations for adjoint sensitivity analysis for hybrid discrete – continuous systems whose dynamics between transition points are described by index-0, index-1, and index-2 Hessenberg DAEs. Using a general framework for describing such systems, we provide sufficient solvability conditions and derive appropriate boundary (final) conditions, as well as the transfer relationships at transition points for the adjoint variables.

Furthermore, we have extended the derivation to include so-called hybrid systems with memory, where the dynamics depend explicitly on the state at the last transition and demonstrated this approach on a computational mechanics problem. Ongoing follow-up work focuses on extensions to more complex dependencies on past history, with application to different material constitutive laws, such as the Menegotto-Pinto model.