\section{Introduction}

We consider the problem of learning input-to-state stable (ISS) dynamics from observed input and output trajectories. Stability constraints provide an inductive bias that can help a dynamics learning algorithm choose one model over another and may improve the generalization accuracy of the dynamical system model when it is simulated in novel interconnection configurations. In some cases, imposing a stability constraint can guide the learning algorithm and improve performance.

This work was motivated by a need for fast-to-simulate behavioral models of electronic circuits.
A system-on-a-chip (SoC), or larger microelectronic system, is composed of multiple functional blocks, often referred to as \emph{IP blocks}, where IP denotes intellectual property. Prior to manufacturing, simulation is used to verify system performance. Each IP block may be represented by its transistor-level \emph{netlist}, which specifies the interconnection of transistors that comprise its nonlinear dynamics.
Transient simulation of the complete model of the non-linear system can be prohibitively slow, so it is preferable to replace the transistor-level description of each IP block with a cheaper-to-simulate behavioral model.


Transient circuit simulation presents a challenge for dynamics learning. At test time, an adaptive timestep solver is used to simulate the learned dynamics model, which is usually interconnected (in feedback) with external circuits and other IP blocks \citep{Hajj2016}. We refer to these external elements collectively as the \emph{load}. In order to be useful, the model must be accurate when interconnected with a variety of loads, and those loads are generally not known a priori. Due to uncertainty over the load, a model that performs well in training and validation may fail when embedded in a circuit simulation at test time. 


In this work, we assume that the circuit of interest is well-modeled by a system of ODEs. 
A natural approach is to directly learn a parametric system of controlled ODEs
\begin{align}
    \dot x &= f(x,u), \label{eq:node1} \\
    y &= h(x), \label{eq:node2}
\end{align}
which has state $x\in\mathbb{R}^n$, input $u\in\mathbb{R}^m$, and output $y\in\mathbb{R}^p$. $u$ and $y$ consist of node voltages and currents, and possibly their time derivatives. 

In this work, we directly learn a \emph{neural ODE} model of the form \eqref{eq:node1} -- \eqref{eq:node2}, which may be trained by either directly backpropagating through an ODE solver or implicitly differentiating through the solution to the ODEs using an adjoint method \citep{ChenRBD2018}. Our models are trained using interpolated trajectories of $u(t)$ and $y(t)$ obtained from a circuit simulator; this approach is similar to prior works on neural ODE models of physical systems with continuous inputs \citep{KidgerMFL2020, ZhongDC2020a}. 

We focus on the case where $f$ takes the form of a continuous-time recurrent neural network (CTRNN) and $h$ is an affine function of the state. Models of this form are universal approximators on finite time intervals \citep{FunahashiN1993} and, in some cases, infinite time intervals \citep{HansonR2020}. In addition, a CTRNN may be implemented as a generic circuit block using the Verilog-A behavioral modeling language and subsequently simulated by commercial circuit simulators \citep{ChenRR2017}.

In Section \ref{sec:parametrization}, we propose a CTRNN parametrization that is guaranteed to be ISS with respect to a Lur'e-Postnikov type (quadratic plus integral) ISS-Lyapunov function $V(x)$, which has parameters that are learned jointly with the CTRNN model parameters. ISS is a natural \textit{a priori} assumption for many circuits; the state is guaranteed to be bounded given bounded inputs, and the state converges to a unique equilibrium if the input is set to zero.

There has been recent interest in learning neural ODEs jointly with a Lyapunov-like function $V(x)$. For autonomous systems, a sufficient condition for global asymptotic stability (GAS) is to ensure that $V$ is strictly decreasing along any system trajectory, i.e.,
\begin{equation}\label{eq:dissipation_gas}
	\dot V(x)<0 \quad \forall x\ne 0.
\end{equation}
\citet{RichardsBK2018} encourage \eqref{eq:dissipation_gas} via regularization, but do not guarantee that the dissipation inequality holds everywhere. \citet{ManekK2019} and \citet{MassaroliPBPYA2020} define the model dynamics as a function of $V(x)$ such that \eqref{eq:dissipation_gas} holds for all $x$. \citet{CranmerGHBSH2020} and \citet{ZhongDC2020a} considered the related problem of learning dynamics with Lagrangian and Hamiltonian structure, respectively. 

Our approach is similar to that of \citet{ManekK2019} in the sense that we guarantee that a dissipation inequality on $V$ holds everywhere, although we consider ISS, which can be seen as a generalization of GAS to systems with inputs. We use a stability condition that generalizes the ISS condition derived by \citet{Ahn2011}, which is based on a quadratic ISS-Lyapunov function. In another related work, \citet{cao2006state} first learn an unconstrained circuit model and subsequently stabilize the model using nonlinear constrained optimization. In contrast, we build the stability constraint directly into the model parametrization.

Besides providing stability guarantees, we observed that our proposed model parametrization can accelerate training convergence. In this sense, it is related to prior works on regularization methods for accelerating neural ODE training. For example, \citet{FinlayJNO2020} penalized the complexity of the model dynamics while \citet{KellyBJD2020} penalized the forward ODE solution time. Unlike those methods, our stability constraint does not introduce additional penalty terms, which can be difficult to tune.

We also show that our ISS parametrization is directly compatible with \emph{aging-aware} circuit modeling. The dynamics of a circuit drift over time due to semiconductor degradation. Aged dynamics, estimated using physics-based approaches \citep{tu1993berkeley}, can be used to verify lifetime specifications and identify aging-induced failures. \citet{rosenbaum2020machine} directly learn an aging-aware circuit model by choosing the dynamics $f$ in \eqref{eq:node1} and output map $h$ in \eqref{eq:node2} to themselves be learned functions of a periodic \emph{stress waveform} $u_{\text{stress}}$, which is assumed to have been applied to the circuit continuously for an operating time $T_{\text{op}}$ on the order of years. Aging analysis can greatly benefit from fast-to-simulate surrogate models since separate aging simulations are needed to characterize different possible use condition profiles, each of which is specified by a pair ($u_{\text{stress}}$, $T_{\text{op}}$).




Section \ref{sec:parametrization} presents our ISS-constrained model and describes how it can be used for transient circuit simulation, with and without aging effects. Section \ref{sec:experiments} evaluates the proposed methods on a variety of circuit modeling tasks.






\section{Input-to-State Stable Continuous-Time Recurrent Neural Networks}\label{sec:parametrization}

\subsection{Continuous-Time Recurrent Neural Networks}\label{subsec:ctrnns}

We consider controlled neural ODEs of the form
\begin{align}
	\dot x &= -\frac{1}{\tau} x + W\sigma_{\ell}(Ax + Bu + \mu) + \nu, \label{eq:1layer_ctrnn} \\
	     y &= Hx + b, \label{eq:ctrnn_output}
\end{align}
where $x\in\mathbb{R}^{n}$ is the state, $u\in\mathbb{R}^m$ is the input, and $y\in\mathbb{R}^p$ is the output. $\tau>0$ is a positive scalar time constant, and $W,A^\top\in\mathbb{R}^{n\times \ell}$, $B\in\mathbb{R}^{\ell\times m}$, $\mu\in\mathbb{R}^\ell$, and $\nu\in\mathbb{R}^n$ are parameters. The element-wise function $\sigma_{\ell}:\mathbb{R}^{\ell}\to\mathbb{R}^{\ell}$ has the form $\sigma_{\ell}(w) = \begin{bmatrix} \sigma(w_1), \dots, \sigma(w_\ell) \end{bmatrix}^\top$, where $\sigma$ is a strictly increasing, continuous, and subdifferentiable scalar-valued nonlinearity that satisfies $\sigma(0)=0$ and the slope condition 
\begin{equation}\label{eq:slope_condition}
	0 \le \frac{\sigma(r)-\sigma(r')}{r-r'} \le 1, \quad \forall r,r'\in\mathbb{R}, r\ne r'.
\end{equation}
Geometrically, \eqref{eq:slope_condition} means that the graph of $\sigma$ lies within a sector in the first and third quadrants, between the horizontal axis and the line with slope one. 
For example, the conditions on $\sigma(\cdot)$ are satisfied by the rectified linear unit $\text{ReLU}(\cdot)=\max\{0,\cdot\}$ and the hyperbolic tangent $\tanh(\cdot)$.

The dynamics \eqref{eq:1layer_ctrnn} may be interpreted as a feedforward neural network with a single hidden layer of dimension $\ell$ and a stabilizing term $-x/\tau$, which is similar to ``skip-connections'' in residual networks \citep{HeZRS2016}. We also assume that $\ell\ge n$. Universal approximation results guarantee that a dynamical system with state dimension $n$ can be approximated arbitrarily well by a CTRNN of the form \eqref{eq:1layer_ctrnn} -- \eqref{eq:ctrnn_output}, as long as $\ell$ is sufficiently large \citep{FunahashiN1993,HansonR2020}. 


\subsection{Input-to-State Stability}

The notion of input-to-state stability (ISS) was developed as a state-space approach to analyzing the stability of systems with inputs \citep{Sontag2008}. Suppose that \eqref{eq:1layer_ctrnn} has an equilibrium point and, without loss of generality, that the equilibrium is at the origin.
\begin{definition}\label{def:iss}
    The system \eqref{eq:node1} is input-to-state stable (ISS) if there exist a class $\mathcal{KL}$ function\footnote{
    	See \citep{Khalil2002} for definitions of class $\mathcal{K}$, $\mathcal{K}_\infty$, and $\mathcal{KL}$ functions.
    } $\beta$ and class $\mathcal{K}_\infty$ function $\gamma$ such that
	\begin{equation}
		\|x(t)\| \le \beta(\|x_0\|, t) + \gamma(\|u\|_\infty)
	\end{equation}
	for all $t\ge 0$, given any bounded input $u:[0,\infty)\to\mathbb{R}^m$ and initial condition $x(0)=x_0$.
\end{definition}
The ISS property captures the idea that bounded inputs result in bounded state. Moreover, the effect of the initial condition on the trajectory (the transient response) should diminish to zero as $t\to\infty$, with rate bounded by the function $\beta$. 
A sufficient condition for ISS can be found by identifying an appropriate ISS-Lyapunov function $V:\mathbb{R}^n\to\mathbb{R}_+$.

\begin{theorem}{\citep{Khalil2002}}\label{thm:iss_lyapunov}
The system \eqref{eq:node1} is ISS if it admits an ISS-Lyapunov function, i.e., a smooth, positive definite, and radially unbounded
function $V$
for which there exist a positive definite function $\alpha$ and class $\mathcal{K}$ function $g$ such that, for bounded inputs $u$,
    \begin{equation}\label{eq:dissipation_iss}
        \dot V(x,u) := \nabla V(x)^\top f(x,u) \le -\alpha(x)\quad\text{if}\,\, \|x\|\ge g(\|u\|).
    \end{equation}
\end{theorem}
The dissipation inequality \eqref{eq:dissipation_iss} ensures that $V$, and therefore $\|x\|$, cannot grow too large relative to the magnitude of the input. Note that in the absence of inputs, $g(0)=0$, and Theorem \ref{thm:iss_lyapunov} reduces to a sufficient condition for global asymptotic stability. In that case, we refer to the associated function $V$ simply as a Lyapunov function.



\subsection{Lyapunov Diagonal Stability Condition}\label{subsec:lds}

\citet{Forti1995} derived a sufficient condition for which \eqref{eq:1layer_ctrnn} is GAS for constant input $u(t)\equiv u_0$. 

\begin{proposition}\label{prop:forti_tesi}
If the matrix $A$ is full rank, i.e. $\emph{rank}(A)=n$, and there exists a positive diagonal matrix $\Omega=\mathrm{diag}(\omega_1,\dots,\omega_\ell)$ with $\omega_i>0$ for each $i=1,\dots,\ell$ such that
\begin{equation}\label{eq:lds_condition}
	\Omega\Big(AW - \frac{1}{\tau}I\Big) + \Big(W^\top A^\top - \frac{1}{\tau}I\Big)\Omega \prec 0,
\end{equation}
then \eqref{eq:1layer_ctrnn} is GAS for constant input $u(t)\equiv u_0$.
\end{proposition}

If the conditions of Proposition \ref{prop:forti_tesi} hold, we say that the matrix $AW-(1/\tau)I$ is \emph{Lyapunov Diagonally Stable} (LDS). The rank condition on $A$ is not restrictive, since the set of rank-deficient $A$ has measure zero. The proof of Proposition \ref{prop:forti_tesi} makes use of a Lyapunov function of the form
\begin{equation}\label{eq:v_lure}
	V(x) = x^\top P x + 2\sum_{i=1}^\ell \omega_i\int_0^{A_i x} (\sigma_\ell)_i (r)\, \mathrm{d}r,
\end{equation}
where $P\succ 0$ is a positive definite matrix, $A_i$ denotes the $i^{\text{th}}$ row of $A$ in \eqref{eq:1layer_ctrnn} for each $i$, and $\omega_i\ge 0$. If $V$ of the form \eqref{eq:v_lure} can be used to prove 0-GAS for \eqref{eq:1layer_ctrnn}, then it it can also serve as an ISS-Lyapunov function.

\begin{proposition}\label{prop:iss_0gas}
	If the conditions in Proposition \ref{prop:forti_tesi} are satisfied, then \eqref{eq:1layer_ctrnn} is also ISS.
\end{proposition}

The direct extension of Proposition \ref{prop:forti_tesi} to ISS is a consequence of the fact that the LDS condition guarantees that \eqref{eq:1layer_ctrnn} is globally exponentially stable when $u\equiv 0$. In general, 0-GAS is a necessary, but not sufficient, condition for ISS. 


\subsection{An Input-to-State Stable Model Parametrization}\label{subsec:iss_new}

Observe that as the matrix $AW$ approaches the zero matrix, the matrix on the left hand side of \eqref{eq:lds_condition} approaches $-\frac{2}{\tau}\Omega$, which is negative definite. Therefore, we may stabilize a given CTRNN by scaling $AW$. Here, we consider a parametrized matrix $A_\theta\in\mathbb{R}^{\ell\times n}$, and define $A$ to be
\begin{equation}\label{eq:parametrization_A}
	A = \frac{1}{\rho(\tau,A_\theta,W,\Omega)+1}A_\theta,
\end{equation}

In the following, let $\lambda_{\max}(M)$ denote the largest eigenvalue of a symmetric matrix $M$. 
\begin{theorem}\label{thm:stable_ctrnn}
For any $\delta>0$, let
\begin{equation}\label{eq:parametrization_rho}
	\rho(\tau,A_\theta,W,\Omega) = \emph{\text{ReLU}}\bigg(\frac{\tau}{2}\lambda_{\max}\big(\Omega^{1/2} A_\theta W\Omega^{-1/2} + \Omega^{-1/2}W^\top A_\theta^\top\Omega^{1/2}\big) - 1 + \delta\bigg),
\end{equation}
Then, \eqref{eq:1layer_ctrnn} with $A$ given by \eqref{eq:parametrization_A} is ISS. 
\end{theorem}

Using \eqref{eq:parametrization_rho}, the scalar $\rho$ is made large enough to ensure that the LDS condition holds, if the condition does not already hold with $A=A_\theta$. The hyperparameter $\delta$ controls the minimum dissipation rate of $V$; smaller values of $\delta$ allow for longer transients. The model may be directly trained using a gradient descent method, since the stability constraint is built into the definition of $A$. $\Omega$ parametrizes $V$ in \eqref{eq:v_lure}, and may be either learned with the rest of the model parameters or fixed, e.g. to $\Omega=I$.






\subsection{Training CTRNN Neural ODEs}\label{subsec:training}

At the start of a transient simulation, a circuit simulator sets the initial condition of its state variables to an equilibrium point given the initial input. For the model \eqref{eq:node1} -- \eqref{eq:node2}, the circuit simulator sets $x(0)=x_0$, where $x_0$ satisfies
\begin{equation}\label{eq:rootfind}
	0 = f(x_0,u(0)).
\end{equation}
This is done using a numerical root-finding scheme, such as the Newton-Raphson method. The equilibrium condition \eqref{eq:rootfind} does not necessarily uniquely define $x_0$, and a user-specified initial condition may need to be provided. Fortunately, if $f$ is given by a CTRNN \eqref{eq:1layer_ctrnn} satisfying the LDS condition, then Proposition \ref{prop:forti_tesi} guarantees the existence of a unique $x_0$ that satisfies \eqref{eq:rootfind} for any $u(0)$. During model training, we set the initial condition by numerically solving \eqref{eq:rootfind}. In order to obtain the necessary derivatives for optimization, we implicitly differentiate through the root-finding operation \citep{BaiKK2019}.

\begin{figure}[htbp]
\centering
\floatconts
  {fig:interconnection}
  {\caption{Simulation diagram for an interconnection between a circuit of interest and a load.}}
  {%
\resizebox{0.45\textwidth}{!}{%
\begin{circuitikz}[line width=0.75pt, american]
    \node[draw, minimum width=2cm, minimum height=2cm, anchor=south west] at (1.5,-.5) {\begin{tabular}{c} Circuit \\ of Interest \end{tabular}};
    \node[draw, minimum width=2cm, minimum height=2cm, anchor=south west] at (5.5,-.5) {\begin{tabular}{c} Load \end{tabular}};
    \draw (8.5,0.9) to[short, i_=$u_{\text{ext}}$] (7.5,0.9);
    \draw (3.75,0.9) to[short, i^<=$u$] (5.5,0.9);
    \draw (5.5,0.1) to[short, i^<=$\tilde y$] (3.75,0.1);
\end{circuitikz}
    }
}
\end{figure}

In our experiments, we consider the simulation setup illustrated in Figure \ref{fig:interconnection}. In a simulation, the circuit of interest is interconnected with a load, which itself may be driven by an external signal $u_{\text{ext}}$. We assume that both the load and $u_{\text{ext}}$ are random and have known distributions. We train our models using $N$ input and output trajectories $u^{(i)}$ and $\tilde y^{(i)}$ for $i=1,\dots,N$. Each pair ($u^{(i)},\tilde y^{(i)})$ is obtained by simulating the system in Figure \ref{fig:interconnection} on a time interval $[0,T]$ with initial condition defined by \eqref{eq:rootfind}, using a random instantiation of the load and $u_{\text{ext}}$. Like \citet{KidgerMFL2020}, we obtain continuous trajectories $u^{(i)}$ and $\tilde y^{(i)}$ by interpolating the solution points provided by the circuit simulator.

The model parameter learning problem is given by the optimization problem
\begin{equation}\label{eq:loss}
\text{minimize}\,\, \frac{1}{N}\sum_{i=1}^N \frac{1}{T}\int_0^T \big(\tilde y^{(i)}(t) - y^{(i)}(t)\big)^2 \, \mathrm{d}t,
\end{equation}
where $y^{(i)}$ is the predicted output. We estimate the integral in \eqref{eq:loss} by a Monte Carlo estimate as follows. Let $S$ be a random variable uniformly distributed on $[0,T]$. Then, we have $\frac{1}{T}\int_0^T \big(\tilde y^{(i)}(t) - y^{(i)}(t)\big)^2 \, dt = \mathbb{E}[\big(\tilde y^{(i)}(S) - y^{(i)}(S)\big)^2]$ for each $i$, and so we may estimate the expectation using $\mathbb{E}[\big(\tilde y^{(i)}(S) - y^{(i)}(S)\big)^2] \approx \frac{1}{K}\sum_{j=1}^K \big(\tilde y^{(i)}(S_j) - y^{(i)}(S_j)\big)^2$, where $S_1,\ldots,S_K$ are i.i.d.\ copies of $S$.


\subsection{Aging-Aware Neural ODEs}\label{subsec:aging}

Conventionally, circuit aging simulation involves two transient simulations of the complete transistor-level netlist \citep{tu1993berkeley}. In the first step, the fresh circuit ($T_{\text{op}}=0$) is simulated subject to  $u_{\text{stress}}$ on a short time horizon $T_{\text{stress}}\llT_{\text{op}}$\footnote{Modern circuits have nanosecond-scale signal periods; $T_{\text{stress}}$ on the order of tens of nanoseconds is usually sufficient.} 
to estimate the per-transistor stress profile. Each transistor's dynamics is subsequently age-adjusted assuming that the $T_{\text{stress}}$-periodic input $u_{\text{stress}}$ is applied for time $T_{\text{op}}$, which is typically on the order of years. In the second step, the circuit is re-simulated using the age-adjusted transistor models, subject to a possibly new input $u$.

\citet{rosenbaum2020machine} proposed a two-step learning approach to learn aging-aware models; we extend that approach to include the ISS constraint. In the first step, we learn a ``fresh'' CTRNN 
\begin{align*}
    \dot x &= -\frac{1}{\tau_0} x + W_0\sigma_\ell(A_0x + B_0u + \mu_0) + \nu_0, \\
    y &= H_0x + b_0,
\end{align*}
corresponding to $T_{\text{op}}=0$ using the approach in Subsection \ref{subsec:training}. In the second step, we fix $\tau_0,W_0,\dots,b_0$ and form an aging-aware CTRNN \eqref{eq:1layer_ctrnn} -- \eqref{eq:ctrnn_output} whose parameters are given by 
\begin{align*}
    \tau &= \tau_0 + \Delta_\tau(u_{\text{stress}},T_{\text{op}}), \\
    W &= W_0 + \Delta_W(u_{\text{stress}},T_{\text{op}}), \\
    &\dots \\
    b &= b_0 + \Delta_b(u_{\text{stress}},T_{\text{op}}),
\end{align*}
where $\Delta_\tau,\Delta_W,\dots,\Delta_b$ are learned parameter perturbation functions. The model structure is suitable because, for realistic use conditions, the stress-induced drift in the dynamics will be relatively small, as illustrated by the example in Figure \hyperref[subfig:predictedinv]{2\textit{(d)}}. An ISS aging-aware model may be obtained by setting $A_\theta = A_0 + \Delta_A(u_{\text{stress}},T_{\text{op}})$ in \eqref{eq:parametrization_A}. The perturbation functions are learned using randomly-generated $u_{\text{stress}}$, $T_{\text{op}}$, and corresponding output trajectories $\tilde y$ obtained from the circuit simulator.




\section{Experiments}\label{sec:experiments}

\subsection{Test Cases}

\paragraph{Common Source Amplifier.} This circuit is a one-transistor, two-port voltage amplifier connected to resistor-capacitor (RC) loads with randomly generated values. The system input $u_{\text{ext}}$ is driven by a random piecewise linear voltage source. The goal is to predict the currents at the input and output ports, given the port voltages. Our CTRNN models had dimensions $n=6$, $\ell=14$, and $m=p=2$.


\paragraph{Continuous-Time Linear Equalizer (CTLE).} The CTLE is a five-transistor differential amplifier with two input ports and two output ports. It is designed to compensate for signal distortion that occurs when digital data are transmitted between two chips. The system input $u_{\text{ext}}$ is given by the output of a pseudorandom bit sequence generator passed through a USB serial link, and the ports are connected to randomly-generated RC loads. The goal is to predict the input port currents and output port voltages, given the input port voltages, their time derivatives, and the output port currents. Our CTRNN models had dimensions $n=20$, $\ell=30$, $m=6$, and $p=4$.


\paragraph{Large IP Block.} 
This test circuit contains between 1000 and 2000 transistors; the exact number is unknown because the circuit is described by an encrypted netlist. Unlike the previous two test cases, this circuit is used with a known, fixed load. The goal is to predict two output voltages given seven input voltages. The inputs are driven by $u_{\text{ext}}$ given by the outputs of seven pseudorandom bit sequence generators. Our CTRNN models had dimensions $n=20$, $\ell=30$, $m=7$, and $p=2$. 



\paragraph{Inverter Chain with Aging Effects.} The last test circuit is a chain of nine cascaded digital inverters; this circuit is often used to benchmark aging. For this test case, aging-induced degradation slows down the dynamics; Figure \hyperref[subfig:predictedinv]{2\textit{(d)}} illustrates the delay between the outputs of a fresh and aged circuit. The model outputs two port currents given two port voltages, given the stress profile $(u_{\text{stress}},T_{\text{op}})$. The system input $u_{\text{ext}}$ was driven by a random piece-wise linear voltage source, and random capacitive loads were connected to the output port. Aging analysis was performed using random piece-wise linear $u_{\text{stress}}$ and random $T_{\text{op}}$ sampled from a log uniform distribution from 0.001 to 10 years. We considered a special case of \eqref{eq:1layer_ctrnn} with $W=I$ and $\nu=0$ with dimensions $n=\ell=20$ and $m=p=2$, and learned parameter perturbations only for $A$, $B$, and $\mu$. $\Delta_A$, $\Delta_B$, and $\Delta_\mu$ were implemented by single-layer gated recurrent unit (GRU) network \citep{cho2014learning} with hidden state dimension 20.


\subsection{Results}

We trained the CTRNNs by directly backpropagating through the order three Bogacki-Shampine ODE solver
with the ADAM optimizer \citep{KingmaB2015}. We took $\sigma_\ell$ to be $\text{ReLU}$ with a bias term, and we used $\delta=10^{-3}$ in the stability constraint \eqref{eq:parametrization_rho}. Model parameters were initialized randomly, with the constraint that \eqref{eq:lds_condition} held with $\Omega=I$.  Each dimension of the inputs and outputs in the training data was separately normalized to $[-1,1]$ prior to training, and the time horizon $T$ was scaled up to be on the order of seconds (instead of nanoseconds). The models used in the first three test cases were trained using the Julia package DiffEqFlux \citep{RackauckasMMCZSSRE2020}; the aging-aware models were trained using the Python package torchdiffeq \citep{ChenRBD2018}. Learned models were implemented in Verilog-A and simulated using the Spectre circuit simulator \citep{Cadence2020}.

Table \ref{table:performance} compares three different training methods: CTRNN with no constraints (Baseline), the stability constraint \eqref{eq:parametrization_A} (Proposed), and the stability constraint with $\Omega$ fixed to be the identity matrix (Proposed, $\Omega=I$). The table shows the mean squared error (MSE) of the predicted model outputs measured on a held-out validation set of input and output waveforms (``Valid.'') and measured when the model is simulated by the circuit simulator as a Verilog-A model (``Test''). The ``Test'' MSE values are averaged across 100 simulations, each with random instantiations of load, $u_{\text{ext}}$, and stress profile ($u_{\text{stress}},T_{\text{op}}$), in the aging-aware inverter chain test case. The MSE for both ``Valid.'' and ``Test'' are computed after applying the aforementioned normalization to $[-1,1]$. Figure \ref{fig:predicted_waveforms} shows example simulations of the proposed ISS CTRNN models, carried out by Spectre.



\begin{table}[htbp]
\centering
\tableconts
{table:performance}
{\caption{MSE metrics. The reported MSE have been multiplied by a factor of 1000.}}
{%
\begin{tabular}{|l||cc||cc||cc||cc|}
\hline
 & \multicolumn{2}{c||}{Amplifier}     & \multicolumn{2}{c||}{CTLE}          & \multicolumn{2}{c||}{IP Block}     & \multicolumn{2}{c|}{Inverter Chain}          \\ 
 \hline
Model Type & \multicolumn{1}{c|}{Valid.} & Test & \multicolumn{1}{c|}{Valid.} & Test & \multicolumn{1}{c|}{Valid.} & Test & \multicolumn{1}{c|}{Valid.} & Test  \\ \hline\hline
Proposed                \rp{0.234}{0.263}   \rp{0.805}{\bf 0.939}    \rp{\bf 0.031}{\bf 0.250}   \rp{0.280}{\bf 0.918}   \\ \hline
Proposed ($\Omega=I$)   \rp{0.260}{0.314}   \rp{1.03}{4.44}    \rp{0.142}{0.2943}   \rp{\bf 0.232}{1.93}    \\ \hline
Baseline                \rp{0.243}{0.279}   \rp{0.714}{3.1447}    \rp{2.93}{15.9}     \rp{0.238}{3.50}        \\ \hline
\end{tabular}
}
\end{table}


Across all test cases, the errors on the held out test set were lower than the errors accumulated when tested in the circuit simulator. This is expected, since the model is simulated in open-loop with the validation set data rather than in closed-loop with the circuit simulator. The proposed constraints uniformly outperformed the baseline learning method, and the stability constraint with learned $\Omega$ almost always outperformed the stability constraint with $\Omega$ fixed to the identity.

For the IP block test case, the ISS constraint was necessary for learning an accurate model. Without the stability constraint, the training struggled to escape a local minimum. This can be seen in the validation losses over the course of training, shown for the IP block and CTLE test cases in Figure \ref{fig:valid_losses}. Although the effect was less pronounced, the ISS constraints were able to stabilize training for the CTLE test case as well.


\section{Conclusion}

In this work, we proposed provably ISS CTRNN neural ODE models for transient circuit simulation. The CTRNN models are constructed in terms of an ISS-Lyapunov function $V$ such that a dissipation inequality on $V$ is satisfied, and can yield better models than baseline unconstrained training approaches. In principle, this type of approach may be extended to impose other dissipation or invariance conditions on $V$, for example conditions based on passivity or energy conservation; that is a suitable subject for future work. 



\begin{figure}[hbp]
\centering
\floatconts
    {fig:predicted_waveforms}
    {\caption{Examples of the proposed ISS CTRNN predictions when simulated by a circuit simulator.}}
    {%
      \subfigure[Amplifier]{%
        \label{subfig:predictedamp}%
        \includegraphics[width=0.48\textwidth]{cs_amp_i2.pdf}
      }
      \subfigure[CTLE]{%
        \label{subfig:predictedctle}%
        \includegraphics[width=0.48\textwidth]{ctle_v4.pdf}
      }
      \subfigure[IP block]{%
        \label{subfig:predictedadi}%
        \includegraphics[width=0.48\textwidth]{adi_2.pdf}
      }
      \subfigure[Inverter Chain]{%
        \label{subfig:predictedinv}%
        \includegraphics[width=0.48\textwidth]{inv9_aging.pdf}
      }
    }
\end{figure}



\begin{figure}[hbp]
\centering
\floatconts
  {fig:valid_losses}
  {\caption{Comparison of validation losses over the course of training.}}
  {%
      \subfigure[IP Block]{%
        \label{subfig:predicted_amp}
        \includegraphics[width=0.45\textwidth]{baseline_vs_constr_fixed_free_Adaptive_ADI_1.pdf}
      }
      \subfigure[CTLE]{%
        \label{subfig:predicted_ctle}
        \includegraphics[width=0.45\textwidth]{baseline_vs_constr_fixed_free_Adaptive_CTLE_1.pdf}
      }
  }
\end{figure}


\acks{This work was funded in part by the NSF under CNS 16-24811 and the industry members of the CAEML I/UCRC, and in part by the Illinois Institute for Data Science and Dynamical Systems (iDS${}^2$), an NSF HDR TRIPODS institute, under award CCF-1934986.}


