The structural index of sensitivity equation systems

ABSTRACT

This work presents a new methodology for computing parameter sensitivities of differential algebraic system of equations with higher differential index. This methodology is particularly adequate for performing sensitivity analysis of object-oriented models described by modern universal modelling languages. By employing the same concepts and tools adopted by these languages for structural analysis of systems of equations, it is shown that the computational graphs of a differential algebraic system of equations and its corresponding sensitivity equation are structurally isomorphic. As a consequence, the structural index of both systems of equations are proven to be equal. Based on this result, an efficient strategy for index reduction of sensitivity equations is designed.

KEYWORDS:

• Modelica
• differential algebraic equations
• structural index

1. Motivation

1.1. Object-oriented modelling and simulation

The recent decades encountered a tremendous progress in the field of multidisciplinary simulation tools for continuous time discrete variable systems. A new generation of universal tools and languages for modelling and simulation multi-physical domain applications as for example Modelica [1], VHDL-AMS [2] and others emerged and became widely accepted. Although this development is still in progress, commercial and academic high-quality implementations are available and broadly used, cf. www.modelica.org/tools.

Two major drivers can be named for this success story: First, from the modelling perspective, the principles of modular domain-independent acausal modelling have been fully developed and powerful modelling languages have been built upon these concepts [3–5]. Second, from the compiler perspective, advanced algorithms for transforming high-level object-oriented models into efficient simulation code via symbolic manipulation and automatic code generation became mature in the recent decades. These algorithms allow to generate large-scale differential algebraic systems of equations (DAEs), automatically simplify and finally transform them into highly efficient simulation code. The combination of such modelling languages and the mentioned algorithms yielded universal object-oriented software tools that are already competitive with domain-specific modelling and simulation tools.

In the context of automatic DAEs generation of object-oriented models, a well-known obstacle is the occurrence of higher index DAEs [5–10]. Briefly spoken, higher index DAEs suffer from a rank deficiency in the subset of algebraic equations which hinders a direct numerical integration using common ODE solvers, solvability of DAE systems is comprehensively presented in [11] (see Section 2.1 for more details). Although some numerical methods for higher index DAE problems exist [12–14], a more general and frequently implemented approach is to perform automatic index reduction. This can be achieved via regularization [15] sometimes using domain-specific knowledge of the model. A universal approach is to combine Pantelides index reduction algorithm [16] with the method of dummy derivatives [17,18]. Using this approach, it is well understood how to reduce higher index DAEs into systems of index one (or zero) by which the well-known drift-problem is avoided. Such DAEs can then be solved by standard solvers like IDA [19] and DASSL [6]. Because the process of index reduction fundamentally relies on symbolically computed derivatives of equations, algorithmic differentiation (AD) techniques for program codes [20,21] are usually involved [22].

A highly related concept to index reduction is the efficient initialization of DAE systems [17,23]. In object-oriented equation-based models, initial values of state variables of components are usually feed as an initial guess of the set of typically over-determined equations composed of

• The hidden algebraic equations implied by the DAE index

• The further algebraic equations resulting by index reduction

A typical approach is to view consistent initialization as an optimization problem attempting to find solution(s) minimizing the residual of the above equations, for example as in [24].

1.2. Index reduction and parameter sensitivity analysis

For the first time, the present contribution links the computational task of DAE index reduction and DAE sensitivity analysis. While sensitivity analysis of DAEs of index one and even index two have been intensively investigated in literature [25,26], the more general case regarding DAEs of arbitrary index has not been considered yet. This work is restricted to DAE systems explicitly given in symbolic form, namely those which has been automatically generated from equation-based object-oriented models. Formally, for a DAE system with an arbitrary index

(1) $F(x^{(k)},x^{(k-1)},\dots ,x^{(1)},x,p,t)=0$(1)

together with a set of consistent initial conditions

$x^{(k-1)}(t_0)=x_{k-1,0}(p),\dots ,x^{(1)}(t_0)=x_{1,0}(p),x(t_0)=x_{0,0}(p)$

where $x^{(i)}(t)=(d^{i}x/d{t}^i)(t)$, $x(t)\in R^n$ are state variables, $p\in R^m$ are model parameters and $F:R^{(k+1)n+m+1}\to R^n$, Parameter Sensitivities (PS) $(\mathrm{\partial }x/\mathrm{\partial }p)(t)\in R^{n\times m}$ are required. These are described by the sensitivity equation composed of the system (1) augmented with the sensitivity subsystem

(2) $F_{x^{(k)}}{x}_p^{(k)}+\dots +F_{x^{(1)}}{x}_p^{(1)}+F_x{x}_p+F_p=0$(2)

together with the initial conditions

$x_p^{(k-1)}(t_0)=\frac{\mathrm{\partial }{x}_{k-1,0}(p)}{\mathrm{\partial }p},\dots ,x_p^{(1)}(t_0)=\frac{\mathrm{\partial }{x}_{1,0}(p)}{\mathrm{\partial }p},x_p(t_0)=\frac{\mathrm{\partial }{x}_0(p)}{\mathrm{\partial }p}$

where $u_v$ stands for the partial derivative of $u$ w.r.t. $v$. General properties of sensitivity equations can be found in [27]. Equation (2) can be automatically derived with the help of AD techniques as done in [28,29]. This sensitivity equation like the original DAE can suffer from a higher differential index. Two strategies can now be imagined for computing PS [28]

1. IR-AD strategy: Index reduction of the original DAE (1) is done first. Then, the resulting index one (or index zero) DAE is algorithmically differentiated w.r.t. its parameters

2. AD-IR strategy: Vice versa, the original DAE is first algorithmically differentiated w.r.t. its parameters. Afterwards, index reduction is applied to the resulting sensitivity Equations (1) and (2).

An example following the IR-AD strategy is the Modelica-based optimization prototyping environment JModelica [30] (www.jmodelica.org). DAEs possibly represented in object-oriented form and combined with an objective optimization function is prototyped with JModelica. The JModelica compiler employs AD techniques to compute the required partial derivatives [31] for index reduction and derivative-based optimization. On the other hand, an example following the AD-IR strategy is the Modelica library ADGenKinetics [32] (www.modelica.org/libraries/adgenkinetics). The library is composed of two subpackages:

• a subpackage containing all components for equation-based modelling of biochemical reaction networks [33] and

• another extended package of identical structure implementing sensitivity equations using equation-based AD technique [28]

Using a Modelica simulation environment, the modelled sensitivity equation is subject to index reduction and numerical integration.

Despite the apparent simplicity and practicability of the IR-AD approach, we show in this contribution that the AD-IR strategy is straightforward and results in a new reliable algorithm for computing PS. In fact, the proper applications of both strategies to a DAE give identical results. The reasons are given by the following facts that will be proven

1. There is a structural isomorphism between the original system (1) and the sensitivity subsystem (2)

2. This isomorphism can be exploited to utilize the Pantelides algorithm [16] in such a way that the index of sensitivity equation can be reduced without having to analyse its structure

3. The structural indices of a DAE and its associated sensitivity equations are always the same even in the case of a higher index

The proposed strategy can be exploited for enabling direct applications of AD techniques to descriptive DAE-based object-oriented models [28,29,34].

The rest of the work is structured as follows. The fundamental tools for dealing with DAEs on a structural level – computational graphs, structural index and Pantelides algorithm – are revised in Section 2. Having these tools available, the previously given statements are exactly proven in Section 3. The presented concepts as well as the AD-IR strategy will be illustrated by demonstrative examples in Section 4. Section 5 ends with a discussion of practical consequences.

2. The structural index

2.1. The index notion

Typically, a given (first-order) ordinary differential system of equations (ODE) is solved by integrating the continuous-time state variables from their derivatives (e.g. using Euler formula). On the other side, due to the potential presence of hidden algebraic constraints among state variables within a given first-order DAE, applying integration alone on the structurally singular DAE generally leads to numerical errors except for special cases, for example [35–38]. A theoretical though inaccurate way (cf. Subsection 2.4) for simulating a DAE system is to transform it to an ODE and then to apply an ODE solver.

For instance and w.l.o.g., given a DAE system of the form:

(3) $F(\dot{x},x,p,t)=0,x(t_0)=x_0(p)$(3)

with state variables $x(t)\in R^n$, model parameters $p\in R^m$, $F:R^{2n+m+1}\to R^n$, consistent initial conditions $x_0(p)\in R^n$. By differentiating (3) w.r.t. time, one obtains:

(4) $F_{\dot{x}}\ddot{x}+F_x\dot{x}+\dot{F}=0$(4)

An ODE could be reformulated if $F_x$ is non-singular. Otherwise, the transformation to an ODE system is done by an iterative combination of differentiations and algebraic substitutions for selected subsets of equations. The minimum possible number of the required differentiation steps to selected equation subsets with which the DAE system can be solved is referred to as the differential index cf. [6,11] for variants of informal definitions. The higher the differential index of a given DAE is, the more difficulties its solution is supposed to impose, although it is not always the case [39]. ODEs as special forms of DAEs correspond to a differential index of zero.

A variant of a definition of differential index is given as follows [6]:

Definition 2.1.

(The Differential Index of a DAE): For a general DAE (3), the differential index along a solution $x(t)$ is the minimum number of differentiations of the system which would be required together with the help of algebraic substitutions to derive $\dot{x}$ uniquely as a continuous function in terms of $x,p$ and $t$. Thus, the index is defined in terms of the over-determined system

(5) $\begin{array}{c}F(\dot{x},x,p,t)=0\\ \frac{dF}{dt}(\dot{x},x,p,t)=0\\ ⋮\\ \frac{d^{q}F}{d{t}^q}(\dot{x},x,p,t)=0\end{array}$(5)

to be the smallest $q$ so that $\dot{x}$ can be solved in terms of $x,p$ and $t$.

The phrase ‘can be solved’ in the last definition does not precisely specify how exactly to derive a solvable system of equations. According to the specific meaning given to this phrase, different index concepts arise as, for example structural index [40], perturbation index [39], tractability index [41], strangeness index [42] and many others [43,44]. ODEs as special forms of DAEs usually correspond to index zero. The higher an index of a given DAE is, the more difficulties its’ solution imposes classified according to the used index. The different indices are in general not equal, although comparative relationship studies have been done as in [43–46]. By targeting the reduction of an index for a given DAE, it is hoped that the resulting system of equations can be solved by common ODE solvers [47].

In the context of automatically generated DAEs from object-oriented models [17], there is no reliable algorithm that correctly calculates the differential index and performs index reduction upon. A computationally efficient variant of the differential index is the structural index practically adequate for performing index reduction. This index attempts to estimate the differential index via maximum matchings of computational graph representation of DAEs (cf. Subsection 2.2 for more details). Nevertheless, while this index purely relies on the topological graph representation of DAEs, the differential index imposes numerical aspects that cannot be captured by the topology alone, e.g. Jacobian singularities. Consequently, there are cases where the structural index does not correctly estimate the differential index [46]. Nevertheless, the structural index concept is a fundamentally established concept to most computational tools. Except for singular cases, it can successfully estimate the differential index in most of the cases for which the statements made in this work are valid.

2.2. Definition of the structural index

The structural index definition is based on computational bipartite graph representations of DAEs defined as follows [3,48]

Definition 2.2.

(Computational graph of a DAE: The computational bipartite Graph $G=(V,E,w)$ representation of a DAE of the form (1) consists of a set of vertices

$V=V_e\cup V_x=\left\{e_1,e_2,\dots ,e_n\right\}\cup \left\{x_1,x_2,\dots ,x_n\right\}$

where $e_i$ and $x_i$ correspond to the $i$-th equation and the $i$-th variable in the given DAE, respectively. The set of edges $E$ is defined as follows

$E=\left\{(e_i,x_j):x_{j}arisesinequation{e}_i\right\}i,j\in \left\{1,2,\dots ,n\right\}$

combined with a weight function $w:E\to N$

$w((e_i,x_j))=\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}{x}_j\in V_x\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}{e}_i\in V_e.$

In the previous definition $G$ is clearly a bipartite graph and $\left|V_e\left|=\right|V_x\right|=n$ for well-determined systems of equations. The structural index is based on maximum matchings [49] of the corresponding computational bipartite graph defined as follows

Definition 2.3.

(Matching, perfect matching, weight of a matching, maximum matching): Given a graph $G=(V,E,w)$, where $V$ and $E$ are sets of vertices and edges, respectively, and given $w:E\to N$ as a weight function, a matching $M^s(G)$ is a set of $s$ edges

$M^s(G)=\left\{d_k:d_k=(e_i,x_j)\in E,k=1,2,..,s\right\}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}|M^s(G)|=s$

s.t. for any $(e_{i_1},x_{j_1})\in M^s(G)$ and $(e_{i_2},x_{j_2})\in M^s(G)\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}{i}_1\ne i_2$ and $j_1\ne j_2$, i.e. no two edges are incident in a vertex. A perfect matching is a matching where all vertices $v\in V$ are covered. The weight of a matching $M^s(G)$ is given by

$w(M^s(G))=\sum _{k=1}^{s}w(d_k),\hspace{0.17em}{d}_k\in M^s(G)$

A maximum matching $M_{max}^s(G)$ of size $s$ is a matching where $w(M_{max}^s(G))\ge w(M_k^s(G))$ $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{s}{M}_k^s(G)$.

A maximum perfect matching of the computational graph from Definition 2.2 is necessarily of size $n$ and can be identified with polynomial-runtime algorithms [50–52]. The structural index is now defined as follows [40]

Definition 2.4

(Structural Index of a DAE-System): Given the DAE (1) with $n\ge 2$ and let $G=(V,E,w)$ be the corresponding computational bipartite graph according to Definition 2.2. The structural index is given by

$I_{str}(G)=w(M_{max}^{n-1}(G))-w(M_{max}^n(G))+1.$

2.3. The pantelides algorithm

Based on the structural index, index reduction can be performed via a derivation of Pantelides algorithm [16] originally motivated by the problem of consistent initialization of higher index DAEs. Algorithm 1, verbally reproduced from [53] (cf. Section 18.2.4), shows a simplified version of Pantelides algorithm for selecting the equations subject to differentiation targeted at reducing the structural index. It accepts a DAE system of arbitrary order as a set of input equations, and it returns a corresponding DAE system of structural index one. If the structural index of the input DAE system is equal to one or zero, nothing needs to be done.

Algorithm 1 Pantelides algorithm for index reduction

1: function $Equations=$ Pantelides($Equations$)

2: $[G=(V,E,w)]\leftarrow$ SetUpComputationalGraph($Equations$) ▹Def. 2.2

3: $[M_{max}^n(G),M_{max}^{n-1}(G),I_{str}(G)]\leftarrow$StructuralIndex($G$)

4: while$I_{str}(G)>1$do

5: $e\leftarrow$ the unique $e\in V_e$ with a $u\in V_x$ but no $v\in V_x$ s.t.

$(e,u)\in M_{max}^n(G)\mathrm{a}\mathrm{n}\mathrm{d}(e,v)\in M_{max}^{n-1}(G)$

6: $Equations(e)$ $\leftarrow d/dt$ $Equations(e)$

▹replace Equation e with the derivative of e w.r.t. t

7: $w((e,v))\leftarrow w((e,v))+1\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}v\in V_{x}s.t.(e,v)\in E$

8: $[M_{max}^n(G),M_{max}^{n-1}(G),I_{str}(G)]\leftarrow$StructuralIndex($G$)

…

12: end while

13:return $Equations$

14: end function

15:

16: function$[M_{max}^n(G),M_{max}^{n-1}(G),I_{str}(G)]$= StructuralIndex$G$

17: $n\leftarrow \left|V\right|/2$

18: $M_{max}^n(G)\leftarrow M_i^n(G)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}w(M_i^n(G))\ge w(M_j^n(G))\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}j$

19: $M_{max}^{n-1}(G)\leftarrow M_i^{n-1}(G)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}w(M_i^{n-1}(G))\ge w(M_j^{n-1}(G))\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}j$

20: $I_{str}(G)\leftarrow w(M_{max}^{n-1}(G))-w(M_{max}^n(G))+1$

21: Return $[M_{max}^n(G),M_{max}^{n-1}(G),I_{str}(G)]$

22: end function

Based on two maximum matchings of size $n$ and $n-1$, i.e. $M_{max}^n(G)$ and $M_{max}^{n-1}(G)$, of the computational graph $G$ of a given DAE, the only equation $e\in V_e$ not present in $M^{n-1}(G)$ is selected for differentiation towards index reduction. This process is repeated until the structural index becomes equal to one. An illustrative example is given in section 4.

2.4. The method of dummy derivatives

Applying numerical solvers on the reduced resulting systems alone via Pantelides algorithm may not lead to accurate numerical results. Numerical results tend to drift away from the original algebraic constraints that have been differentiated by the process of index reduction [53]. The numerical integration of systems of differential equations at each time step resolves to problem of nonlinear system equation where hidden constraints are not explicitly included. In practice, the Pantelides algorithm is combined with the method of dummy derivatives [18]. This method rather augments a higher-index DAE with the differentiated equations selected by the Pantalides algorithm. However for each new generated equation, a new dependent algebraic variable is introduced replacing the derivative of a selected variable (i.e. a dummy derivative). In this way, the generated reduced system of equations is well-determined in terms of the number of unknowns. Overall, the index gets reduced while the original algebraic constraints are kept in the system of equations and hence they are still subject to satisfaction by the numerical solvers. Algorithm 2 lists a modification of Algorithm 1 that combines the method of dummy derivatives with Pantelides algorithm. Here, line 6 from Algorithm 1 is modified and lines 9–11 from Algorithm 2 are inserted in Algorithm 1. Section 4 provides two examples demonstrating the index reduction process.

Algorithm 2 Pantelides algorithm combined with the method of dummy derivatives

1: function $Equations$ = Pantelides($Equations$)

…

4: while$I_{str}(G)>1$ do

…

6: $Equations\leftarrow Equations{\cup }^{E}quations(e)$

$Equations(e)\leftarrow \frac{d}{dt}Equations(e)$

▹augment $Equations$ with the derivative of $e$

…

9: $u\leftarrow$ select one element from $\left\{v:(e,v)\in E\right\}$

10: $m=w((e,u))$

11: $Equations\leftarrow Equations(u_d\leftarrow \frac{d^{m}u}{d{t}^m})$

▹replace occurrences of $d^{m}u/d{t}^m$ with a new variable $u_d$

12: end while

13: return $Equations$

14: end function

3. Structural characteristics of a sensitivity equation

3.1. Isomorphism of computational graphs

In this section, it is shown that the structural index of sensitivity Equations (1) and (2) is equal to that of the original system (1). This is done by utilizing common structural characteristics between the systems of Equations (1) and (2). Namely, their corresponding computational graphs are shown to be isomorphic [54]. The following definition is needed to prove the stated arguments.

Definition 3.1

(The induced subgraph): Given a graph $G=(V,E,w)$ where $V$ is a set of vertices, $E$ is a set of edges and $w:E\to N$ is a weight function and suppose that $U\subset V$, then the induced subgraph of $G$ with the vertex set $U$ is defined as

$G\left[U\right]=(U,E_U,w_U)$

where $E_U=\left\{(u_i,u_j)\in E:u_i\in U\mathrm{a}\mathrm{n}\mathrm{d}{u}_j\in U\right\}$ and $w_U:E_U\to N$ s.t. $w_U(u_i,u_j)=w(u_i,u_j)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}(u_i,u_j)\in E_U$.

Theorem 3.2:

Let $G=(V,E,w)$ be the computational bipartite graph representation (cf. Definition 2.2) of the DAE (1) and let $G^{\ast }=(V^{\ast },E^{\ast },w^{\ast })$ be the bipartite graph representation of the corresponding sensitivity Equation (1) and (2), where

$V^{\ast }=V\cup V^{\prime }=V_x\cup V_e\cup V_x^{\prime }\cup V_e^{\prime }$

with $V_x,V_e$ as in Definition 2.2 and

$V_x^{\prime }=\left\{{x^{\prime }}_1,{x^{\prime }}_2,\dots ,{x^{\prime }}_n\right\},\hspace{0.17em}{V}_e^{\prime }=\left\{{e^{\prime }}_1,{e^{\prime }}_2,\dots ,{e^{\prime }}_n\right\}$

where $e_i^{\prime }$ ($x_i^{\prime }$) is the derivative of equation $e_i$ (variable $x_i$) w.r.t. $p$ combining $p$ equations (variables) in one vertex. The weight function $w^{\ast }:E^{\ast }\to N$ is defined as follows:

$w^{\ast }((e_i^{\ast },x_j^{\ast }))=Thehighesttimederivativeorderof{x}_j^{\ast }\in V_x^{\ast }inequation{e}_i^{\ast }\in V_e^{\ast }$

Consider the graph

$G^{\prime }=G^{\ast }\left[V^{\prime }\right]=(V^{\prime },E^{\prime },w^{\prime })$

then $G^{\prime }\cong G$, that is $G$ and $G^{\prime }$ are isomorphic.

Proof. Clearly, $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{y}{x}_i\in V_x,\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}{x}_i^{\prime }\in V_x^{\prime }$ and similarly, $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{y}{e}_i\in V_e,\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}{e}_i^{\prime }\in V_e^{\prime }$ for $i=1,2,...,n$. Additionally, $\mathrm{a}\mathrm{n}\mathrm{y}(e_i,x_j)\in E\mathrm{i}\mathrm{f}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{y}\mathrm{i}\mathrm{f}(e_i^{\prime },x_j^{\prime })\in E^{\prime }\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{y}(e_i,x_j)\in E\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}w(e_i,x_j)=w^{\prime }(e_i^{\prime },x_j^{\prime })$. Hence, the graph $G^{\prime }=(V^{\prime },E^{\prime },w^{\prime })=G^{\ast }[V^{\prime }]\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{s}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}Gi.e.G^{\prime }\cong G$ by definition of isomorphism □

3.2. The structural index of sensitivity equation

The previous theorem is in general a way to conclude that the structural characteristics remain preserved in the derived sensitivity Equation (2). Since the definition of the structural index is based on maximum matchings in the bipartite graph representation, it can be concluded that the image of a maximum matching $M_{max}^n(G)$ in $G^{\prime }$ is also a maximum matching. Consequently, it can be shown that the structural index of the DAE (1) is equal to the structural index of the sensitivity subsystem (2) as shown in the followings using the same notations in Theorem 3.2. Note that a perfect matching in $G^{\ast }$ is necessarily of size $2n$.

Corollary 3.3:

(1) $M^s(G)$ is a matching

$implies{M^{\prime }}^s(G^{\prime })=\left\{(e^{\prime },{x^{\prime }}_i):(e,x_i)\in M^s(G)\right\}isamatching$

(2) $M^n(G)$ is a maximum matching $implies{M^{\prime }}^n(G^{\prime })$ is a maximum matching

(3) The weight of a maximal matching of $G$ is equal to the weight of a maximum matching of $G^{\prime }$, i.e. $w(M_{max}^n(G))=w^{\prime }(M_{max}^n(G^{\prime }))$

Proof. Straightforward from Theorem 3.2 □

Corollary 3.4:

A maximum perfect matching of $G^{\ast }$ is composed of a maximum perfect matching in $G$ and a maximum matching in $G^{\prime }$, i.e.

$M_{max}^{\ast 2n}(G^{\ast })=M_{max}^n(G)\cup M_{max}^{\prime n}(G^{\prime })$

Proof. $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{y}\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}d\in M_{max}^{\ast 2n}(G^{\ast }),d=(e^{\ast },x^{\ast })$, $e^{\ast }\in V_e^{\ast }$ and $x^{\ast }\in V_x^{\ast }$, there are three cases concerning $d$:

First case: Either $d\in E$ or $d\in E^{\prime }$.

In this case, the proof is straightforward using Corollary 3.3 and Theorem 3.2.

Second case: $e^{\ast }\in V_e$ and $x^{\ast }\in V_x^{\prime }$.

Then the variable $x^{\ast }$ representing a $\frac{\mathrm{\partial }{x}_i}{\mathrm{\partial }p},i\in \left\{1,2,...,n\right\}$ is present in an equation $e^{\ast }\in V_e$ in Equation (2) $impliescontradiction$

Third case: $e^{\ast }\in V_e^{\prime }$ and $x^{\ast }\in V_x$.

Since $M_{max}^{\ast 2n}(G^{\ast })$ is a perfect matching, $\left|{V^{\prime }}_x\left|=\right|{V^{\prime }}_e\right|$ and any $x^{\prime }\in V_x^{\prime }$ must be present exactly in one edge $u\in M_{max}^{\ast 2n}(G^{\ast })$ $impliesthereisu\in M_{max}^{\ast 2n}(G^{\ast })$ s.t. $u=(e,x^{\prime })$, where $e\in V_e,x^{\prime }\in V_x^{\prime }$ $impliescontradiction(accordingtosecondcase)$ □

Similarly, a matching of size $2n-1$ in $G^{\ast }$ necessarily contains one maximum perfect matching either in $G$ or $G^{\prime }$.

Corollary 3.5:

A maximum match of size $2n-1$ in $G^{\ast }$ is either composed of:

• a maximum perfect matching in $G$ and a maximum matching of size $n-1$ in $G^{\prime }$ or

• a maximum matching of size $n-1$ in $G$ and a maximum perfect matching in $G^{\prime }$

i.e. either

$M_{max}^{\ast 2n-1}(G^{\ast })=M_{max}^n(G)\cup M_{max}^{\prime n-1}(G^{\prime })$

or

$M_{max}^{\ast 2n-1}(G^{\ast })=M_{max}^{n-1}(G)\cup M_{max}^{\prime n}(G^{\prime })$

is true.

Proof. Suppose that $e^{\ast }\in V_e{\cup }^{V_e^{\prime }}$ is the only element in $V$ s.t. there is no $x^{\ast }\in V_x{\cup }^{V_x^{\prime }}$ with $(e^{\ast },x^{\ast })\in M_{max}^{\ast 2n-1}(G^{\ast })$. Then there are two cases

First case: Suppose that $e^{\ast }\in V_e$, then since $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}e\in V_e$ and $e\ne e^{\ast }\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{i}\mathrm{s}x\in V_x$ s.t. $(e,x)\in M_{max}^{\ast 2n-1}(G^{\ast })\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}{M}_{max}^{n-1}(G)\subset M_{max}^{\ast 2n-1}(G^{\ast })$. Moreover, since $w^{\ast }((e^{\prime },x))=w^{\prime }((e^{\prime },x^{\prime }))\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}{e}^{\prime }\in V_e^{\prime },x\in V_x$ s.t. $(e^{\prime },x)\in E^{\ast }\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}{M}_{max}^{\ast 2n-1}(G^{\ast })=M_{max}^{n-1}(G)\cup M_{max}^{\prime n}(G^{\prime })$.

Second case: $e^{\ast }\in V_e^{\prime }$ analogous to first case using Theorem 3.2 $\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{s}{M}_{max}^{\ast 2n-1}(G^{\ast })=M_{max}^n(G)\cup M_{max}^{\prime n-1}(G^{\prime })$ □

Using the presented corollaries, the following corollary is introduced.

Corollary 3.6:

(1) The weight of a maximal perfect matching in $G^{\ast }$ is twice the weight of a maximum perfect matching in $G$, i.e.

$w^{\ast }(M_{max}^{\ast 2n}(G^{\ast }))=2\cdot w(M_{max}^n(G))$

(2) Similarly, the weight of a matching of size $2n-1$ in $G^{\ast }$ is equal to the sum of the weight of a maximum perfect matching and a matching of size $n-1$ in $G$, i.e.

$w^{\ast }(M_{max}^{\ast 2n-1}(G^{\ast }))=w(M_{max}^n(G))+w(M_{max}^{n-1}(G))$

Proof. Straightforward using Theorem 3.2 and corollaries 3.3 3.4 3.5. □

Now, the main result of this work can be shown: .

Theorem 3.7:

The structural index of $G^{\ast }$ is equal to the structural index of $G$, i.e.

$I_{str}(G^{\ast })=I_{str}(G)$

Proof. $I_{str}(G^{\ast })=w(M_{max}^{\ast 2n-1}(G^{\ast }))-w(M_{max}^{\ast 2n}(G^{\ast }))+1$ (Definition 2.4)

$=w(M_{max}^n(G))+w(M_{max}^{n-1}(G))-2\cdot w(M_{max}^n(G))+1$ (Corollary 3.6)

$=w(M_{max}^{n-1}(G))-w(M_{max}^n(G))+1$

$=I_{str}(G)$ □

3.3. Index reduction of sensitivity equation

With this theorem, it can be concluded that when applying AD-IR strategy for computing PS of a DAE of the form (1), there is no need to apply Pantelides algorithm on the large sensitivity subsystem (2). For index reduction of sensitivity equation, the lines 5, 6, 10 and 11 in Algorithm 2 are modified to the corresponding statements in Algorithm 3. Algorithm 3 expects both of the DAE system (1) (Equations) and its corresponding sensitivity subsystem (2) (Equations’) as inputs and returns the reduced systems of equations of both systems of equations as outputs. Comprehensive algorithms for computing the sensitivity equation subsystems (Equations’) are given in [28]. When an equation $e$ in (1) is selected for index reduction, it is enough to consider the equation $e^{\prime }$ in (2) for the index reduction as well. Similarly, when a variable $u$ is selected as a dummy derivative, the corresponding variable $u^{\prime }$ is also selected as a dummy derivative as well.

Algorithm 3 Extended version of Pantelides algorithm for index reduction of sensitivity equation

1: function $[Equations,Equation{s}^{\prime }]$ = Pantelides($Equations,Equation{s}^{\prime }$)

…

4: while $I_{str}(G)>1$ do

5: $e\leftarrow$

$e^{\prime }\leftarrow$ the corresponding equation to $e$ in $Equation{s}^{\prime }$

6: $Equations\leftarrow \dots$

$Equation{s}^{\prime }\leftarrow Equation{s}^{\prime }{\cup }^{E}quation{s}^{\prime }(e)$

$Equation{s}^{\prime }(e)\leftarrow \frac{d}{dt}Equation{s}^{\prime }(e)$

▹augment $Equation{s}^{\prime }$ with the derivative of $e^{\prime }$

…

9: $u\leftarrow$ select one element from $\left\{v:(e,v)\in E\right\}$

$u^{\prime }\leftarrow$ the derivative of $u$ w.r.t. $p$

10: $m=w((e,u))$

11: $Equations$ $\leftarrow$ …

$Equation{s}^{\prime }\leftarrow Equation{s}^{\prime }(u_d^{\prime }\leftarrow \frac{d^m{u}^{\prime }}{d{t}^m},u_d\leftarrow \frac{d^{m}u}{d{t}^m})$

▹replace occurrences of $d^m{u}^{\prime }/d{t}^m$ with a new variable $u_d^{\prime }$

12: end while

13: return [$Equations,Equation{s}^{\prime }$]

14: end function

4. Illustrative examples

4.1. The electrical circuit model example

Figure 1(a) shows an electrical circuit model (cf [53]. Section 18.2.3.2). The model can be assembled to the following equation system:

$e_1:C_1{\dot{V}}_1+V_1/R_1-I_1=0$

$e_2:C_2{\dot{V}}_2+V_2/R_2-I_2=0$

$e_3:V_1-V_2=0$

$e_4:I_1+I_2-I=0$

$I_1(0)=I_{1,0},V_1(0)=V_{1,0}$

Figure 1. An electrical circuit model and its computational graphs representations before and after one iteration of Algorithm 2.

To simplify the illustration, the parameter values of capacitances, i.e. $C_1$ and $C_2$, and resistances, i.e. $R_1$ and $R_2$, are fixed to the value one. The only parameter remaining is the input current $I$. The voltage and the current through an electrical circuit with $I$ as a fixed input current is described through the following DAE

$e_1:{\dot{V}}_1+V_1-I_1=0$

$e_2:{\dot{V}}_2+V_2-I_2=0$

$e_3:V_1-V_2=0$

$e_4:I_1+I_2-I=0$

$I_1(0)=I_{1,0},V_1(0)=V_{1,0}$

where $e_i$s are equation labels. Additionally, its computational graphs before and after applying Algorithm 2 are shown in Figure 1(b). It can be shown that the electrical circuit DAE turns out to have a differential index of two.

4.1.1. Index reduction with the structural index

To accurately solve this system of equations using, say, a DAE solver, it needs to get transformed into a DAE of index one including the algebraic constraints. Figure 1 shows the computational graphs of the system of equations before and after one iteration of Algorithm 2. Using Definition 2.4 and denoting the computational graph $G$, the structural index of the system of equations of the electrical circuit is computed as follows

$M_{max}^3(G)=\left\{(e_1,V_1),(e_2,V_2),(e_4,I_2)\right\}$

$M_{max}^4(G)=\left\{(e_1,V_1),(e_2,I_2),(e_3,V_2),(e_4,I_1)\right\}$

$\Rightarrow I_{str}(G)=2-1+1=2$

By applying Algorithm 1 on the system of equations one obtains

$e_1:{\dot{V}}_1+V_1-I_1=0$

$e_2:{\dot{V}}_2+V_2-I_2=0$

$e_3:{\dot{V}}_1-{\dot{V}}_2=0$

$e_4:I_1+I_2-I=0$

The above system of equations is overdetermined if the algebraic constraint $V_1-V_2=0$ is included (i.e. five equations in four unknowns). By rather applying Algorithm 2, a consistent system of equations is constructed via the replacement of the time derivative ${\dot{V}}_2$ with an algebraic variable $V_{d2}$

$e_1:{\dot{V}}_1+V_1-I_1=0$

$e_2:V_{d2}+V_2-I_2=0$

$e_3:{\dot{V}}_1-V_{d2}=0$

$V_1-V_2=0$

$e_4:I_1+I_2-I=0$

Note that the new generated equation is not included in the underlying computational graph. The algorithm influences the computational graph only by incrementing the weight of edges incident to the selected equation. Now it can be shown that the structural index of the underlying computational graph $G$ is equal to one

$M_{max}^3(G)=\left\{(e_1,V_1),(e_2,I_2),(e_3,V_2)\right\}$

$M_{max}^4(G)=\left\{(e_1,V_1),(e_2,I_2),(e_3,V_2),(e_4,I_1)\right\}$

$\Rightarrow I_{str}(G)=2-2+1=1$

The previous system of equations has structural index one and it can be solved via a DAE solver.

4.1.2. Structural index of sensitivity equation

Using Algorithm 3, the structural index of sensitivity equation

$e_1:{\dot{V}}_1+V_1-I_1=0e_1^{\prime }:{{\dot{V}}^{\mathrm{\prime }}}_1+V_1^{\prime }-I_1^{\prime }=0$

$e_2:{\dot{V}}_2+V_2-I_2=0e_2^{\prime }:{{\dot{V}}^{\mathrm{\prime }}}_2+V_2^{\prime }-I_2^{\prime }=0$

$e_3:V_1-V_2=0e_3^{\prime }:V_1^{\prime }-V_2^{\prime }=0$

$e_4:I_1+I_2-I=0e_4^{\prime }:I_1^{\prime }+I_2^{\prime }-1=0$

with $X^{\prime }=\mathrm{\partial }X/\mathrm{\partial }I$ can be reduced in a straightforward way as follows

$e_1:{\dot{V}}_1+V_1-I_1=0e_1^{\prime }:{{\dot{V}}^{\mathrm{\prime }}}_1+V_1^{\prime }-I_1^{\prime }=0$

$e_2:V_{d2}+V_2-I_2=0e_2^{\prime }:V_{d2}^{\prime }+V_2^{\prime }-I_2^{\prime }=0$

$e_3:{\dot{V}}_1-V_{d2}=0e_3^{\prime }:{{\dot{V}}^{\mathrm{\prime }}}_1-V_{d2}^{\prime }=0$

$V_1-V_2=0V_1^{\prime }-V_2^{\prime }=0$

$e_4:I_1+I_2-I=0e_4^{\prime }:I_1^{\prime }+I_2^{\prime }-1=0$

It is easy to see that the AD-IR strategy results in the same system of equations as the IR-AD strategy. There is no need to apply Pantelides algorithm on the whole computational graph of the sensitivity equation.

4.2. The pendulum model example

Figure 2(a) demonstrates the motion of a free pendulum in the Cartesian space (cf [53]. Section 12.7.2). This is described by a second-order system of equations

$e_1:m\ddot{x}+\frac{x}{L}F=0$

$e_2:m\ddot{y}+\frac{y}{L}F+mg=0$

$e_3:x^2+y^2-L^2=0$

Figure 2. The pendulum model and its computational graphs representation for two iterations of Algorithm 2.

Additionally, its computational graphs after applying Algorithm 2 for two iterations are shown in Figure 2(b).

4.2.1. The structural index of the pendulum system of equations

It can be shown that the Pendulum system of equations is of differential index three (i.e. it requires three differentiation steps to get transformed into an explicit ODE system). The structural index can be shown to be equal to three as follows

$M_{max}^2(G)=\left\{(e_1,x),(e_2,y)\right\},M_{max}^3(G)=\left\{(e_1,x),(e_2,F),(e_3,y)\right\}$

$\Rightarrow I_{str}(G)=4-2+1=3$

4.2.2. Structural index of sensitivity equation

The sensitivity equation of the pendulum system of equations w.r.t. the parameters $p=(m,L{)}^T$ is as follows

$e_1:m\ddot{x}+\frac{x}{L}F=0e_1^{\prime }:m^{\prime }\ddot{x}+m{\ddot{x}}^{\prime }+\frac{x^{\prime }L-L^{\prime }x}{L^2}F+\frac{x}{L}{F}^{\prime }=0$

$e_2:m\ddot{y}+\frac{y}{L}F+mg=0e_2^{\prime }:m^{\prime }\ddot{y}+m{\ddot{y}}^{\prime }+\frac{y^{\prime }L-L^{\prime }y}{L^2}F+\frac{y}{L}{F}^{\prime }+m^{\prime }g=0$

$e_3:x^2+y^2-L^2=0e_3^{\prime }:x{x}^{\prime }+y{y}^{\prime }-L{L}^{\prime }=0$

where $X^{\prime }=\mathrm{\partial }X/\mathrm{\partial }p$. The structural index of the sensitivity equation is reduced to one by applying Algorithm 3. Two iterations are performed in which equations $e_3$ and $e_3^{\prime }$ are differentiated two times. After one iteration, the sensitivity equation becomes

$e_1:m\ddot{x}+\frac{x}{L}F=0e_1^{\prime }:m^{\prime }\ddot{x}+m{\ddot{x}}^{\prime }+\frac{x^{\prime }L-L^{\prime }x}{L^2}F+\frac{x}{L}{F}^{\prime }=0$

$e_2:m\ddot{y}+\frac{y}{L}F+mg=0e_2^{\prime }:m^{\prime }\ddot{y}+m{\ddot{y}}^{\prime }+\frac{y^{\prime }L-L^{\prime }y}{L^2}F+\frac{y}{L}{F}^{\prime }+m^{\prime }g=0$

$e_3:x\dot{x}+y{y}_d=0e_3^{\prime }:\dot{x}{x}^{\prime }+x{\dot{x}}^{\prime }+y_d{y}^{\prime }+y{y}_d^{\prime }=0$

$x^2+y^2-L^2=0x{x}^{\prime }+y{y}^{\prime }-L{L}^{\prime }=0$

In the previous system of equations $\dot{y}$ and ${\dot{y}}^{\prime }$ were replaced by $y_d$ and $y_d^{\prime }$, respectively, according to the method of dummy derivatives. Alternatively, a dummy variable $x_d$ and its derivatives $x_d^{\prime }$ could have been introduced mainly satisfying line 5 in Algorithm 1 and line 9 in Algorithm 2. This will result in a rather different equation system, however, similarly leading to a reduced-index system after executing the main loop in Algorithm 3. After the second iteration, the sensitivity equation becomes as follows

$e_1:m{x}_d+\frac{x}{L}F=0e_1^{\prime }:m^{\prime }{x}_d+m{x}_d^{\prime }+\frac{x^{\prime }L-L^{\prime }x}{L^2}F+\frac{x}{L}{F}^{\prime }=0$

$e_2:m\ddot{y}+\frac{y}{L}F+mg=0e_2^{\prime }:m^{\prime }\ddot{y}+m{\ddot{y}}^{\prime }+\frac{y^{\prime }L-L^{\prime }y}{L^2}F+\frac{y}{L}{F}^{\prime }+m^{\prime }g=0$

$e_3:x{x}_d+{\dot{x}}^2{e}_3^{\prime }:x_d{x}^{\prime }+2\dot{x}{\dot{x}}^{\prime }+x{x}_d^{\prime }$

$+\dot{y}{y}_d+y{\dot{y}}_d=0+{\dot{y}}_d{y}^{\prime }+y_d{\dot{y}}^{\prime }+\dot{y}{y}_d^{\prime }+y\dot{y}{ }_d^{\prime }=0$

$x\dot{x}+y{y}_d=0\dot{x}{x}^{\prime }+x{\dot{x}}^{\prime }+y_d{y}^{\prime }+y{y}_d^{\prime }=0$

$x^2+y^2-L^2=0x{x}^{\prime }+y{y}^{\prime }-L{L}^{\prime }=0$

Similarly, $\ddot{x}$ and ${\ddot{x}}^{\prime }$ were replaced by $x_d$ and $x_d^{\prime }$, respectively. The resulting system can be solved by a DAE solver if it is transformed into a first-order system (by introducing a dummy variable $v_y=\dot{y}$).

5. Conclusion

In this work, a novel result concerning the structural index of sensitivity equation is presented. The relationship between the differential index of a DAE and that of the corresponding sensitivity equation is emphasized. It is concluded that both DAEs have an identical structural index. This implies that the differential index of both systems of equations is equal for a considerable class of systems, namely those DAEs with structural index identical to the differential index.

Due to the common structural characteristics of a DAE (1) and its corresponding sensitivity Equations (1) and (2), the index reduction process is computationally equal in terms of efforts in the context of sensitivity equation. As concluded by Theorem (3.2), if an equation $e_i\in V_e$ is selected for differentiation towards structural index reduction, the corresponding equations ${e_i}^{\prime }\in {V_e}^{\prime }$, obtained by differentiating $e_i$ w.r.t. $p$, should be also chosen for index reduction of the whole sensitivity subsystem of Equation (2) (cf. Theorem 3.2 for the used notions).

As Definition 2.2 does not exclude DAEs of higher time derivatives, regardless of the fact that they can be rewritten as DAEs of order one, the conclusions made in this work are also applicable to DAEs of higher orders. This can be observed from the second-order system of equations in the presented example. The same statements also apply to sensitivity equation w.r.t. some model start values by which derivatives can be also present.

These theoretical results encourage the usage of the AD-IR strategy for performing sensitivity analysis of high-index differential algebraic equations. From one side, the IR-AD strategy is typically employed as the case with JModelica [30]. The compiler capabilities of JModelica enhanced with AD techniques [31] are exploited for linking generated index-reduced DAEs with highly specialized solvers within the SUNDIALS suite [19]. The underlying state-of-the-art solver is based on a modern method for performing sensitivity analysis of DAEs. It exploits the structure of the sensitivity equation and the linearity of its extended Jacobian for efficient numerical integration. Moreover, specialized interfaces are provided for automatically differentiated routines (i.e. computer programs that have been automatically generated by AD tools, cf. www.autodiff.org) computing the required partial derivatives in Equation (2).

On the other side, the non-conventional AD-IR is an attractive choice for sensitivity analysis of high-level abstraction models. An example is given by the Modelica library ADGenKinetics [32], which includes algorithmically differentiated components for modelling biochemical reaction networks. By simulating base models additionally embedding sensitivity equation, PSs are evaluated using arbitrary DAE solvers provided by the chosen Modelica simulation environment. From one side, direct numerical integration of large-scale systems becomes inefficient. From the other side, the advantages of employing modern DAE solvers for sensitivity analysis vanish with increasing nonlinearity of large-scale sensitivity systems [55]. Nevertheless, an easily applicable modification of the parallelizable method from Dickinson and Gelinas [56] can be used for retaining the efficiency of the AD-IR strategy. In this method, several smaller copies of sensitivity equation w.r.t. few parameters are simulated in parallel as shown in [57] Section 16.3. Based on this work, similar tools and libraries based on the AD-IR strategy for AD of object-oriented models can be designed in future, for example the ADMSL library [58].

A significant practical contribution of this work is to assess a potential enhancement to the Modelica language w.r.t. to the operator der, expressing time derivatives of state variables, e.g. der(x) $\equiv dx/dt$. A way to express PS at Modelica language level is to allow a second parameter argument for the operator der to take place, i.e. der(x,p) $\equiv \mathrm{\partial }x/\mathrm{\partial }p$. This work shows that such an enhancement would not increase the compilation algorithmic complexity at least w.r.t. the process of index reduction. Another interesting point that needs to be investigated is whether the structure of sensitivity equations can be exploited to handle its robust initialization efficiently.

Acknowledgments

The first author acknowledges Austrian Institute of Technology as this article has been partially written during my employment there.

Disclosure statement

No potential conflict of interest was reported by the authors.