# Phy-Taylor: Physics-Model-Based Deep Neural Networks

Yanbing Mao  Engineering Technology Division, Wayne State University, Detroit,
MI 48201, USA  corresponding author: Yanbing Mao (e-mail<EMAIL_ADDRESS>Lui Sha  Department of Computer Science, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  Huajie Shao  Department of Computer Science,
College of William & Mary, Williamsburg, VA 23185, USA  Yuliang Gu  Department
of Mechanical Engineering, University of Illinois at Urbana–Champaign, Urbana,
IL 61801, USA  Qixin Wang  Department of Computing, Hong Kong Polytechnic
University, Hong Kong SAR, China  Tarek Abdelzaher  Department of Computer
Science, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

Purely data-driven deep neural networks (DNNs) applied to physical engineering
systems can infer relations that violate physics laws, thus leading to
unexpected consequences. To address this challenge, we propose a physics-
model-based DNN framework, called Phy-Taylor, that accelerates learning
compliant representations with physical knowledge. The Phy-Taylor framework
makes two key contributions; it introduces a new architectural physics-
compatible neural network (PhN), and features a novel compliance mechanism, we
call Physics-guided Neural Network Editing. The PhN aims to directly capture
nonlinearities inspired by physical quantities, such as kinetic energy,
potential energy, electrical power, and aerodynamic drag force. To do so, the
PhN augments neural network layers with two key components: (i) monomials of
Taylor series expansion of nonlinear functions capturing physical knowledge,
and (ii) a suppressor for mitigating the influence of noise. The neural-
network editing mechanism further modifies network links and activation
functions consistently with physical knowledge. As an extension, we also
propose a self-correcting Phy-Taylor framework that introduces two additional
capabilities: (i) physics-model-based safety relationship learning, and (ii)
automatic output correction when violations of safety occur. Through
experiments, we show that (by expressing hard-to-learn nonlinearities directly
and by constraining dependencies) Phy-Taylor features considerably fewer
parameters, and a remarkably accelerated training process, while offering
enhanced model robustness and accuracy.

##  1 Introduction

The paper proposes a novel physics-model-based deep neural network framework,
called Phy-Taylor, that addresses a critical flaw in purely data-driven neural
networks, when used to model aspects of physical engineering systems. Namely,
it addresses the potential lack of agreement between learned latent neural
network representations and prior physical knowledge – a flaw that sometimes
leads to catastrophic consequences [1]. As shown in Figure 1, the Phy-Taylor
framework introduces two contributions: the deep physics-compatible neural
networks and a physics-guided neural network editing mechanism, aiming at
ensuring compliance with prior physical knowledge.

Figure 1: Architectures of Phy-Taylor and physics-compatible neural network,
having neural network (NN) editing including link editing and activation
editing.

The work contributes to emerging research on physics-enhanced deep neural
networks. Current approaches include physics-informed neural networks [2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16], physics-guided neural-network
architectures [17, 18, 19, 20, 21, 22] and physics-inspired neural operators
[23, 24]. The physics-informed networks and physics-guided architectures use
compact partial differential equations (PDEs) for formulating loss functions
and/or architectural components. Physics-inspired neural operators, such as
the Koopman neural operator [23] and the Fourier neural operator [24], on the
other hand, map nonlinear functions into alternative domains, where it is
easier to train their parameters from observational data and reason about
convergence. These frameworks improve consistency with prior analytical
knowledge, but remain problematic in several respects. For example, (i) due to
incomplete knowledge, the compact or precise PDEs may not always be available,
and (ii) fully-connected neural networks can introduce spurious correlations
that deviate from strict compliance with available well-validated physical
knowledge. Instead, through the use of a Taylor-series expansion, the Phy-
Taylor is able to leverage partial knowledge. Moreover, thanks to the neural
editing mechanism, the framework removes links and reshapes activation
functions not consistent with physics-based representations.

The Phy-Taylor framework leverages the intuition that most physical relations
live in low-dimensional manifolds, shaped by applicable physical laws. It’s
just that the estimation of key physical variables from high-dimensional
system observations is often challenging. By expressing known knowledge as
relations between yet-to-be-computed latent variables, we force representation
learning to converge to a space, where these variables represent desired
physical quantities, shaped by the applicable (expressed) physical laws. In
effect, by shaping non-linear terms and relations in the latent space, we
arrive at a desired physics-compliant latent representation. More
specifically, Phy-Taylor offers the following two advantages:

  * •

Non-linear Physics Term Representation: Classical neural networks can learn
arbitrary non-linear relations by unfolding them into layers of linear
weighting functions and switch-like activations. This mechanism is akin to
constructing nonlinearities by stitching together piecewise linear behaviors.
Instead, by directly exploiting non-linear terms of the Taylor series
expansion, we offer a set of features that express physical nonlinearities
much more succinctly, thereby reducing the number of needed parameters and
improving accuracy of representation. Monomials of the Taylor series can
capture common nonliearities present in physics equations, such as kinetic
energy, potential energy, rolling resistance and aerodynamic drag force. The
(controllable) model error of the series drops significantly as the series
order increases [25]. The approach constructs input features that represent
monomials of the Taylor series and adds a compressor for mitigating influence
of noise on augmented inputs.

  * •

Removing Spurious Correlations: The general topology of neural networks allows
for models that capture spurious correlations in training samples
(overfitting) [26, 27]. In contrast, we develop a neural network (topology)
editing mechanism in the latent space that removes links among certain latent
variables, when these links contradict their intended physical behaviors,
thereby forcing the latent representation to converge to variables with the
desired semantic interpretation that obey the desired physical relations.

Through experiments with learning the dynamics of autonomous vehicles and
other non-linear physical systems, we show that Phy-Taylor exhibits a
considerable reduction in learning parameters, a remarkably accelerated
training process, and greatly enhanced model robustness and accuracy (viewed
from the perspective of long-horizon prediction of a trajectory). Experiments
with safe velocity regulation in autonomous vehicles further demonstrate that
the self-correcting Phy-Taylor successfully addresses the dilemma of
prediction horizon and computation time that nonlinear model-predictive
control and control barrier function are facing in safety-critical control.

##  2 Problem Formulation

Table 1: Table of Notation $\mathbb{R}^{n}$: set of n-dimensional real vectors |  $\mathbb{R}_{\geq 0}$: set of non-negative real numbers
---|---
$\mathbb{N}$: set of natural numbers |  $[\mathbf{x}]_{i}$: $i$-th entry of vector $\mathbf{x}$
$[\mathbf{x}]_{i:j}$: a sub-vector formed by the $i$-th to $j$-th entries of vector $\mathbf{x}$ |  $[\mathbf{W}]_{i,j}$: element at row $i$ and column $j$ of matrix $\mathbf{W}$
$[\mathbf{W}]_{i,:}$: $i$-th row of matrix $\mathbf{W}$ |  $\left[\mathbf{x}\leavevmode\nobreak\ ;\leavevmode\nobreak\ \mathbf{y}\right]$: stacked (tall column) vector of vectors $\mathbf{x}$ and $\mathbf{y}$
$\mathbf{0}_{n}$: $n$-dimensional vector of all zeros |  $\mathbf{1}_{n}$: $n$-dimensional vector of all ones
$\mathbf{O}_{m\times n}$: $m\times n$-dimensional zero matrix |  $\mathbf{I}_{n}$: $n\times n$-dimensional identity matrix
$||\cdot||$: Euclidean norm of a vector or absolute value of a number |  $\odot$: Hadamard product
$\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.8}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.8}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.8}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.8}{$\scriptscriptstyle\bullet$}}}}}$: multiplication operator |  $\mathrm{len}(\mathbf{x})$: length of vector $\mathbf{x}$
act: activation function |  sus: suppressor function
$\top$: matrix or vector transposition |  ina: a function that is inactive
$\boxplus$: a known model-substructure parameter |  $*$: an unknown model-substructure parameter

Consider the problem of computing some output vectors, $\mathbf{y}$, from a
set of observations, $\mathbf{x}$. The relation between $\mathbf{x}$ and
$\mathbf{y}$ is partially determined by physical models of known structure
(but possibly unknown parameter values) and partially unknown, thus calling
for representation learning of the missing substructures using neural network
observables. For example, $\mathbf{y}$ might denote the estimated momentum and
future position of a target as a function of a vector of $\mathbf{x}$, that
include its position, velocity, and type. In this formulation, position and
velocity might be directly related to output quantities via known physical
relations, but type is represented only indirectly by an image that requires
some representation learning in order to translate it into relevant parameters
(such as mass and maneuverability) from which the outputs can be computed. We
express the overall input/output relation by the function:

$\displaystyle\mathbf{y}=\underbrace{\mathbf{A}}_{\text{weight
matrix}}\cdot\underbrace{\mathfrak{m}(\mathbf{x},r)}_{\text{node-
representation vector}}+\underbrace{\mathbf{f}(\mathbf{x})}_{\text{model
mismatch}}\triangleq\underbrace{\mathbf{g}(\mathbf{x})}_{\text{ground truth
model}},$ (1)

where $\mathbf{y}$ and $\mathbf{x}$ are the output and input vectors of
overall system model, respectively, and the parameter $r\in\mathbb{N}$
controls model size. For convenience, Table 1 summarizes the remaining
notations used throughout the paper. Since Equation (1) combines known and
unknown model substructures, we distinguish them according to the definition
below.

######  Definition 1.

For all $i\in\\{1,2,\ldots,\mathrm{len}(\mathbf{y})\\}$,
$j\in\\{1,2,\ldots,\mathrm{len}(\mathfrak{m}(\mathbf{x},r))\\}$, element
$[\mathbf{A}]_{i,j}$ is said to be a known model-substructure parameter in
Equation (1) if and only if
$\frac{{\partial[\mathbf{f}(\mathbf{x})]_{i}}}{{\partial{[\mathfrak{m}}(\mathbf{x},r)]_{j}}}\equiv
0$. Otherwise, it is called an unknown model-substructure parameter.

Definition 1 indicates that the model-substructure knowledge includes:

  * •

(i) Known parameter values but completely unknown model formula (see e.g., the
Example 1).

  * •

(ii) Partially-known model formula. For example, in the (lateral) force
balance equation of autonomous vehicle [28]:

$\displaystyle
m\left({\ddot{y}+\dot{\psi}{V_{x}}}\right)={F_{\text{yf}}}+{F_{\text{fr}}}+{F_{\text{bank}}},$
(2)

the force formula due to road bank angle $\phi$, i.e.,
${F_{\text{bank}}}=mg\sin(\phi)$, is known while other force formulas are
unknown because of complex and unforeseen driving environments.

  * •

(iii) Known model formula but unknown parameter values.

######  Example 1 (Identify Known Model-Substructure Parameters).

DNN Design Goal: Use the current position $p(k)$, velocity $v(k)$ and mass $m$
to estimate a vehicle’s next velocity $v(k+1)$, sensed road friction
coefficient $r(k+1)$ and safety metric $s(k+1)$ of velocity regulation, in the
dynamic driving environments. With the knowledge of vehicle dynamics and
control [28], the problem can be mathematically described by

$\displaystyle v\left({k+1}\right)$
$\displaystyle={g_{1}}\left({p\left(k\right),\leavevmode\nobreak\
v\left(k\right),\leavevmode\nobreak\ m}\right),$ (3a) $\displaystyle
r\left({k+1}\right)$ $\displaystyle={g_{2}}\left({m,\leavevmode\nobreak\
v\left(k\right)}\right),$ (3b) $\displaystyle s\left({k+1}\right)$
$\displaystyle={g_{3}}\left({v\left(k\right)}\right),$ (3c)

We here assume the formulas of ${g_{1}}(\cdot)$–${g_{3}}(\cdot)$ are unknown.
While the Equations (3b) and (3c) indicate that given the inputs, the
available knowledge are (i) the road friction coefficient varies with a
vehicle’s mass and real-time velocity only, and (ii) the considered safety
metric depends on velocity only. We let $\mathbf{x}=[p(k);\leavevmode\nobreak\
v(k);\leavevmode\nobreak\ m]$, $\mathbf{y}=[v(k+1);\leavevmode\nobreak\
r(k+1);\leavevmode\nobreak\ s(k+1)]$,
$\mathbf{g}(\mathbf{x})=[{g_{1}}({p(k),v(k),m});\leavevmode\nobreak\
{g_{2}}({m,v(k)});\leavevmode\nobreak\ {g_{3}}({v(k)})]$, and
$\mathfrak{m}(\mathbf{x},r)=[1;\leavevmode\nobreak\ p(k);\leavevmode\nobreak\
v(k);\leavevmode\nobreak\ m;\leavevmode\nobreak\ p^{2}(k);\leavevmode\nobreak\
p(k)v(k);\leavevmode\nobreak\ mp(k);\leavevmode\nobreak\
v^{2}(k);\leavevmode\nobreak\ mv(k);\leavevmode\nobreak\ m^{2}]$. The ground
truth model (3) is then equivalently rewritten in the form of (1):

$\displaystyle\mathbf{y}=\underbrace{\left[{\begin{array}[]{*{20}{c}}*&*&*&*&*&*&*&*&*&*\\\
&0&*&*&0&0&0&*&*&*\\\
&0&*&0&0&0&0&*&0&0\end{array}}\right]}_{\mathbf{A}}\cdot\mathfrak{m}(\mathbf{x},r)+\underbrace{\mathbf{g}(\mathbf{x})-\mathbf{A}\cdot\mathfrak{m}(\mathbf{x},r)}_{\mathbf{f}(\mathbf{x})}=\mathbf{g}(\mathbf{x}),$
(7)

which thus encodes the available knowledge points (i) and (ii) to the known
model-substructure parameters (i.e., zeros) in system matrix $\mathbf{A}$.

Considering this definition, the problem addressed in this paper is formally
stated below.

######  Problem 1.

Given a time-series of inputs, $\mathbf{x}$, the corresponding outputs,
$\mathbf{y}$, and the known model substructures in Equation (1), it is desired
to develop an end-to-end neural network that directly estimates $\mathbf{y}$
(denoted by $\widehat{\mathbf{y}}$), given $\mathbf{x}$, consistently with all
known model substructures. In other words, the model must satisfy the property
that for each known model-substructure parameter, $[\mathbf{A}]_{i,j}$, the
end-to-end model must ensure that
$\frac{{\partial{{[\widehat{\mathbf{y}}}]_{i}}}}{{\partial{[\mathfrak{m}}\left({\mathbf{x},r}\right)]_{j}}}\equiv[\mathbf{A}]_{i,j}$
for any $\mathfrak{m}(\mathbf{x},r)$.

The above definition allows the system described by Equation (1) to have an
end-to-end model that intertwines well-known substructure properties with
high-order unmodeled correlations of unknown nonlinear structure. In this, our
problem differs from past seminal frameworks of physics-enhanced DNNs [2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 29,
30, 31, 32, 33, 34, 35], that use a compact partial differential equations
(PDEs) for formulating the PDEs-regulated loss function and/or DNN
architectures to count the degree mean of consistency with PDEs. The proposed
solution to Problem 1 is the Phy-Taylor framework, which will rely on two
building blocks: a deep physics-compatible neural network (PhN) and a physics-
guided neural network editing mechanism, presented in the next section.

##  3 Phy-Taylor Framework

The proposed Phy-Taylor for addressing Problem 1 is shown in Figure 1, which
is built on the conjunctive deep physics-compatible neural network (PhN) and
physics-guided neural network (NN) editing. In other words, implementing NN
editing according to Taylor’s theorem for embedding available physical
knowledge into deep PhN yields the Phy-Taylor. The PhN is a neural network
layer with a key component: a physics-inspired augmentation (called Phy-
Augmentation) for generating monomials in Equation (1) of Taylor series
expansion of nonlinear functions capturing physical knowledge. The physics-
guided NN editing – including link editing and activation editing – further
modifies network links and activation functions consistently with physical
knowledge. Specifically, the link editing performs removing and preserving
links according to the consistency with physical knowledge. Meanwhile, the
activation editing performs the physics-knowledge-preserving computing in
output channel of each PhN. Collaboratively through link and activation
editing, the input/output of Phy-Taylor strictly complies with the available
physical knowledge, which is a desired solution to Problem 1. Next, we detail
the two components.

###  3.1 The Physics-compatible Neural Network (PhN)

In order to capture non-linear features of physical functions, we introduce a
new type of network layer that is augmented with terms derived from Taylor
series expansion. The Taylor’s Theorem offers a series expansion of arbitrary
nonlinear functions, as shown below.  Taylor’s Theorem (Chapter 2.4 [25]): Let
$\mathbf{g}\\!:\leavevmode\nobreak\ \mathbb{R}^{n}\to\mathbb{R}$ be a
$r$-times continuously differentiable function at the point
$\mathbf{o}\in\mathbb{R}^{n}$. Then there exists
$\mathbf{h}_{\alpha}\\!:\leavevmode\nobreak\ \mathbb{R}^{n}\to\mathbb{R}$,
where $\left|\alpha\right|=r$, such that
$\displaystyle\mathbf{g}(\mathbf{x})=\sum\limits_{\left|\alpha\right|\leq
r}{\frac{{{\partial^{\alpha}}\mathbf{g}(\mathbf{o})}}{{\alpha!}}}{\left({\mathbf{x}-\mathbf{o}}\right)^{\alpha}}+\sum\limits_{\left|\alpha\right|=r}{{\mathbf{h}_{\alpha}}(\mathbf{x}){{({\mathbf{x}-\mathbf{o}})^{\alpha}}}},\hskip
5.69046pt\text{and}\leavevmode\nobreak\
\mathop{\lim}\limits_{\mathbf{x}\to\mathbf{o}}{\mathbf{h}_{\alpha}}\left(\mathbf{x}\right)=\mathbf{0},$
(8) where $\alpha=\left[{{\alpha_{1}};\leavevmode\nobreak\
{\alpha_{2}};\leavevmode\nobreak\ \ldots;\leavevmode\nobreak\
{\alpha_{n}}}\right]$,
$\left|\alpha\right|=\sum\limits_{i=1}^{n}{{\alpha_{i}}}$,
$\alpha!=\prod\limits_{i=1}^{n}{{\alpha_{i}}}!$,
${\mathbf{x}^{\alpha}}=\prod\limits_{i=1}^{n}{\mathbf{x}_{i}^{{\alpha_{i}}}}$,
and
${\partial^{\alpha}}\mathbf{g}=\frac{{{\partial^{\left|\alpha\right|}}\mathbf{g}}}{{\partial\mathbf{x}_{1}^{{\alpha_{1}}}\cdot\ldots\cdot\partial\mathbf{x}_{n}^{{\alpha_{n}}}}}$.
The Taylor’s theorem has several desirable properties:

  * •

Non-linear Physics Term Representation: The high-order monomials (i.e., the
ones included in $\left({\mathbf{x}-\mathbf{o}}\right)^{\alpha}$ with
$|\alpha|\geq 2$) of the Taylor series (i.e.,
$\sum\limits_{\left|\alpha\right|\leq
r}{\frac{{{D^{\alpha}}\mathbf{g}(\mathbf{o})}}{{\alpha!}}}{\left({\mathbf{x}-\mathbf{o}}\right)^{\alpha}}$)
capture core nonlinearities of physical quantities such as kinetic energy
($\triangleq\frac{1}{2}m{v^{2}}$), potential energy
($\triangleq\frac{1}{2}k{x^{2}}$), electrical power ($\triangleq V\cdot I$)
and aerodynamic drag force ($\triangleq\frac{1}{2}\rho{v^{2}}{C_{D}}A$), that
drive the state dynamics of physical systems.

  * •

Controllable Model Accuracy: Given ${\mathbf{h}_{\alpha}}(\mathbf{x})$ is
finite and $\left\|{\mathbf{x}-\mathbf{o}}\right\|<1$, the error
$\sum\limits_{\left|\alpha\right|=r}{{\mathbf{h}_{\alpha}}(\mathbf{x}){{({\mathbf{x}-\mathbf{o}})^{\alpha}}}}$
for approximating the ground truth $\mathbf{g}(\mathbf{x})$ will drop
significantly as the order $r=|\alpha|$ increases and
$\mathop{\lim}\limits_{|\alpha|=r\to\infty}{\mathbf{h}_{\alpha}}(\mathbf{x}){\left({\mathbf{x}-\mathbf{o}}\right)^{\alpha}}=\mathbf{0}$.
This allows for controllable model accuracy via controlling order $r$.

  * •

Knowledge Embedding: The Taylor series can directly project the known model
substructure parameters of the ground-truth model (1) into neural network
parameters including the weight matrix
(${\frac{{{D^{\alpha}}\mathbf{g}(\mathbf{o})}}{{\alpha!}}}$ with $|\alpha|>0$)
and bias (${\frac{{{D^{\alpha}}\mathbf{g}(\mathbf{o})}}{{\alpha!}}}$ with
$|\alpha|=0$), thus paving the way to embed the available physical knowledge
in the form of an appropriately weighted neural network layer.

Figure 2: Phy-Augmentation architecture.

We note the Taylor theorem relies on an assumption that the ground truth
$\mathbf{g}(\mathbf{x})$ is a $r$-times continuously differentiable function
at the point $\mathbf{o}$. If the assumption does not hold, what the Taylor
series will approximate is a proximity of ground truth that is $r$-times
continuously differentiable. For continuous functions, this is often a
sufficient approximation. Next, we describe how PhNs embed the Taylor series
expansion into neural network layers. The resulting architecture (of a single
PhN layer) is shown in Figure 1. Compared with a classical neural network
layer, we introduce the Phy-Augmentation, whose architecture is shown in
Figure 2. The Phy-Augmentation has two components: (i) augmented inputs that
represent monomials of a Taylor series expansion, and (ii) a suppressor for
mitigating the influence of noise on such augmented inputs (i.e., high-order
monomials). Next, we detail them.

####  3.1.1 Phy-Augmentation: Taylor Series Monomials

Input: augmentation order $r$, input $\mathbf{x}$, point $\mathbf{o}$,
suppressor mapping $\chi(\cdot)$.

1 Suppress input:
$[\mathbf{x}]_{i}\leftarrow\begin{cases}\chi([{\mathbf{x}}]_{i}),&\text{if
suppressor is active}\\\ [{\mathbf{x}}]_{i},&\text{otherwise}\end{cases}$,
$i\in\\{1,2,\ldots,\text{len}(\mathbf{x})\\}$;

2 Generate index vector of input entries:
$\mathbf{i}\leftarrow[1;\leavevmode\nobreak\ 2;\leavevmode\nobreak\
\ldots;\leavevmode\nobreak\ \mathrm{len}({\mathbf{x}})]$;

3 Generate augmentations:
${\mathfrak{m}}({\mathbf{x}},r)\leftarrow{\mathbf{x}}$;

4 for  _$\\_\ =2$ to $r$_ do

5 for  _$i=1$ to $\mathrm{len}({\mathbf{x}})$ _ do

6 Compute temporaries: ${\mathbf{{t}}}_{a}$
$\leftarrow[{\mathbf{x}}]_{i}\cdot[{\mathbf{x}}]_{\left[{[\mathbf{i}]_{i}\leavevmode\nobreak\
:\leavevmode\nobreak\ \mathrm{len}({\mathbf{x}})}\right]}$;

7 if  _$i==1$ _ then

8 Generate temporaries:
$\widetilde{\mathbf{{t}}}_{b}\leftarrow\widetilde{\mathbf{{t}}}_{a}$;

9 else

10 Generate temporaries:
$\widetilde{\mathbf{{t}}}_{b}\leftarrow\left[\widetilde{\mathbf{{t}}}_{b};\leavevmode\nobreak\
\widetilde{\mathbf{{t}}}_{a}\right]$;

11 end if

12 Update index entry: $[\mathbf{i}]_{i}\leftarrow\mathrm{len}({\mathbf{x}})$;

13 Update augmentations:
${\mathfrak{m}}({\mathbf{x}},r)\leftarrow\left[{\mathfrak{m}}({\mathbf{x}},r);\leavevmode\nobreak\
{{\mathbf{{t}}}_{b}}\right]$;

14

15 end for

16

17 end for

Output vector of augmented monomials:
${\mathfrak{m}}({\mathbf{x}},r)\leftarrow\left[1;\leavevmode\nobreak\
{\mathfrak{m}}({\mathbf{x}},r)\right]$.

Algorithm 1 Phy-Augmentation Procedure Figure 3:  An example of Algorithm 1 in
TensorFlow framework, where input $\tilde{\mathbf{x}}\in\mathbb{R}^{3}$ is
from Line 1 of Algorithm 1.

The function of physical-features augmentation in Figure 2 is to generate the
vector of physical features (i.e., node representations) in form of Taylor
series monomials, which is formally described by Algorithm 1. The Lines 1–1 of
Algorithm 1 guarantee that the generated node-representation vector embraces
all the non-missing and non-redundant monomials of Taylor series. The Line 1
shows that Algorithm 1 finally stacks vector with one. This operation means a
PhN node will be assigned to be one, and the bias (corresponding to
${\frac{{{D^{\alpha}}\mathbf{g}(\mathbf{o})}}{{\alpha!}}}$ with $|\alpha|=0$
in Taylor series) will be thus treated as link weights in PhN layers. As an
example shown in Figure 3, the Phy-Augmentation empowers PhN to well capture
core nonlinearities of physical quantities (see e.g., kinetic energy,
potential energy, electrical power and aerodynamic drag force) that drive the
state dynamics of physical systems, and then represent or approximate physical
knowledge in form of the Taylor series.

We note the Line 1 of Algorithm 1 means the noise suppressor is not applied to
all the input elements. The critical reason is the extra mapping induced by
suppressor on inputs can destroy the compliance with physical knowledge, when
the available model-substructure knowledge do not include the mapping of
suppressor. We next is going to present the function of suppressor.

####  3.1.2 Phy-Augmentation: Noise Suppressor

The suppressor in Figure 2 is to mitigate the influence of noise on the
augmented high-order monomials. Before proceeding on the working mechanism, we
present a metric pertaining to noise and true data.

######  Definition 2.

Consider the noisy data and define the data-to-noise ratio ($\mathrm{DNR}$):

$\displaystyle[\bar{\mathbf{x}}]_{i}=\underbrace{[\mathbf{h}]_{i}}_{\text{true
data}}+\underbrace{[{\mathbf{w}}]_{i}}_{\text{noise}}\in\mathbb{R},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{DNR}_{i}\triangleq\frac{[\mathbf{h}]_{i}}{[\mathbf{w}]_{i}}.$ (9)

The auxiliary Theorem 4 presented in Supplementary Information 9.1 implies
that the high-order monomials can shrink their DNRs due to nonlinear mapping.
This means the PhN can be vulnerable to the noisy inputs, owning to Phy-
Augmentation for generating high-order monomials. Hence, mitigating the
influence of noise is vital for enhancing the robustness of PhNs, consequently
the Phy-Taylor. As shown in Figure 2, we incorporate a suppressor into PhN to
process the raw input data, such that the high-order monomial from Phy-
Augmentation can enlarge their DNRs. Building on Definition 2, the proposed
noise suppressor mapping is

$\displaystyle\chi([\bar{\mathbf{x}}]_{i})=\chi([\mathbf{h}]_{i}+[\mathbf{w}]_{i})=\begin{cases}0,&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}<0\\\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i},&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}<0\\\
([\mathbf{h}]_{i}+[\mathbf{w}]_{i})\cdot\kappa_{i}+\rho_{i},&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}>0\end{cases},$ (10)

where the parameters $\rho_{i}$ and $\kappa_{i}$ satisfy

$\displaystyle|\rho_{i}|\geq|[\mathbf{h}]_{i}+[\mathbf{w}]_{i}|\cdot|\kappa_{i}|.$
(11)

We next present the suppressor properties in the following theorem, whose
proof appears in Supplementary Information 9.2.

######  Theorem 1.

Consider the noisy data $[\bar{\mathbf{x}}]_{i}$ and the suppressor described
in Equations (9) and (10), respectively. Under the condition (11), the
suppressor output, denoted by
$[\widehat{\mathbf{x}}]_{i}=\chi([\bar{\mathbf{x}}]_{i})$, has the properties:

The DNR magnitude of high-order monomial
$[\widehat{\mathbf{x}}]_{i}^{p}[\widehat{\mathbf{x}}]_{j}^{q}$ ($p+q\geq 2$)
is strictly increasing with respect to DNR magnitudes of
$[\widehat{\mathbf{x}}]_{i}$ and $[\widehat{\mathbf{x}}]_{j}$. (12) The true
data and the noise of suppressor output $[\widehat{\mathbf{x}}]_{i}$ are
$\displaystyle[\widetilde{\mathbf{h}}]_{i}=\begin{cases}\\![\mathbf{h}]_{i}\cdot\kappa_{i}+\rho_{i},\\!\\!&[\mathbf{h}]_{i}\\!+\\![\mathbf{w}]_{i}\\!\geq\\!0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}\\!>\\!0\\\
\\![\mathbf{h}]_{i},\\!\\!&\text{otherwise}\\\
\end{cases},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
[\widetilde{\mathbf{w}}]_{i}=\begin{cases}\\!-[\mathbf{h}]_{i},\\!\\!&[\mathbf{h}]_{i}\\!+\\![\mathbf{w}]_{i}\\!<\\!0\\\
\\![\mathbf{w}]_{i},\\!\\!&[\mathbf{h}]_{i}\\!+\\![\mathbf{w}]_{i}\\!\geq\\!0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}\\!<\\!0\\\
\\![\mathbf{w}]_{i}\cdot\kappa_{i},\\!\\!&[\mathbf{h}]_{i}\\!+\\![\mathbf{w}]_{i}\\!\geq\\!0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}\\!>\\!0\end{cases},$ (13)

such that
$[\widehat{\mathbf{x}}]_{i}=[\widetilde{\mathbf{h}}]_{i}+[\widetilde{\mathbf{w}}]_{i},\leavevmode\nobreak\
i\in\\{1,2,\dots,\mathrm{len}(\widehat{\mathbf{x}})\\}$.

The result (13) implies the parameters $\kappa$ and $\rho$ control the DNRs of
suppressed data, consequently the high-order monomials. Furthermore, the
result (12) suggests that through designing parameters $\kappa_{i}$,
$\rho_{i}$, $\kappa_{j}$ and $\rho_{j}$ for increasing the DNR magnitudes of
data $[\widehat{\mathbf{x}}]_{i}$ and $[\widehat{\mathbf{x}}]_{j}$, the DNR of
high-order monomial
$[\widehat{\mathbf{x}}]_{i}^{p}[\widehat{\mathbf{x}}]_{j}^{q}$ can be enlarged
consequently, such that the influence of noise is mitigated.

###  3.2 Physics-guided Neural Network Editing

Input: Available knowledge included in system matrix $\mathbf{A}$ of ground-
truth model (1), terminal output dimension $\text{len}(\mathbf{y})$, number
$p$ of PhNs, activation functions $\text{act}(\cdot)$,
$\mathbf{y}_{\left\langle 0\right\rangle}=\mathbf{x}$ and $r_{\left\langle
1\right\rangle}=r$.

1 for  _$t=1$ to $p$_ do

2 if  _$t==1$ _ then

3 Deactivate noise suppressor;

4 Generate node-representation vector $\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle})$ via Algorithm 1;

5 Generate knowledge matrix $\mathbf{K}_{\left\langle t\right\rangle}$:
$[\mathbf{K}_{\left\langle
t\right\rangle}]_{i,j}\leftarrow\begin{cases}[\mathbf{A}_{\left\langle
t\right\rangle}]_{i,j},&[\mathbf{A}_{\left\langle
t\right\rangle}]_{i,j}=\boxplus\\\ 0,&\text{otherwise}\end{cases}$;

6 Generate weight-masking matrix $\mathbf{M}_{\left\langle t\right\rangle}$:
$[\mathbf{M}_{\left\langle
t\right\rangle}]_{i,j}\leftarrow\begin{cases}0,&[\mathbf{A}_{\left\langle
t\right\rangle}]_{i,j}=\boxplus\\\ 1,&\text{otherwise}\end{cases}$;

7 Generate activation-masking vector $\mathbf{a}_{\left\langle
t\right\rangle}$:  $[\mathbf{a}_{\left\langle
t\right\rangle}]_{i}\\!\leftarrow\\!\begin{cases}0,\\!\\!\\!&[\mathbf{M}_{\left\langle
t\right\rangle}]_{i,j}\\!=\\!0,\forall
j\\!\in\\!\\{1,\ldots,\text{len}(\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle}))\\}\\\
1,\\!\\!\\!&\text{otherwise}\end{cases}$;

8

9 else

10 Generate node-representation vector $\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle})$ via Algorithm 1;

11 Generate knowledge matrix $\mathbf{K}_{\left\langle t\right\rangle}$:
$\displaystyle\mathbf{K}_{\left\langle
t\right\rangle}\leftarrow\left[\begin{gathered}\leavevmode\hbox
to270.7pt{\vbox to68.29pt{\pgfpicture\makeatletter\hbox{\hskip 0.2pt\lower
0.2pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{74.8947pt}{14.4014pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mbox{\small
O}_{(\text{len}(\mathbf{y}_{\left\langle
t\right\rangle})-\text{len}(\mathbf{y}))\times\text{len}(\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle}))}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{13.22885pt}{47.09738pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\mathbf{0}_{\text{len}(\mathbf{y})}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{73.00107pt}{47.24461pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mbox{\small
I}_{\text{len}(\mathbf{y})}$}} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{139.64648pt}{48.54462pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mbox{\small
O}_{\text{len}(\mathbf{y})\times(\text{len}(\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle
t\right\rangle}))-\text{len}(\mathbf{y})-1)}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{34.14322pt}\pgfsys@lineto{270.30121pt}{34.14322pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{51.21504pt}{34.14322pt}\pgfsys@lineto{51.21504pt}{68.28644pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setdash{3.0pt,3.0pt}{0.0pt}\pgfsys@invoke{
}{}\pgfsys@moveto{113.81104pt}{34.14322pt}\pgfsys@lineto{113.81104pt}{68.28644pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{gathered}\right];$
(15)

12 Generate weight-masking matrix $\mathbf{M}_{\left\langle t\right\rangle}$:
$[\mathbf{M}_{\left\langle
t\right\rangle}]_{i,j}\leftarrow\begin{cases}0,&\frac{{\partial[\mathfrak{m}(\mathbf{y}_{\left\langle
t\right\rangle},r_{\left\langle
t\right\rangle})]_{j}}}{{\partial[\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})]_{v}}}\neq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{M}_{\left\langle
1\right\rangle}]_{i,v}=0,\leavevmode\nobreak\ \leavevmode\nobreak\
v\in\\{1,2,\dots,\text{len}(\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle}))\\}\\\ 1,&\text{otherwise}\end{cases}$;

13 Generate activation-masking vector $\mathbf{a}_{\left\langle
t\right\rangle}\leftarrow\left[\mathbf{a}_{\left\langle
1\right\rangle};\leavevmode\nobreak\
\mathbf{1}_{\text{len}(\mathbf{y}_{\left\langle
t\right\rangle})-\text{len}(\mathbf{y})}\right]$;

14

15 end if

16 Generate weight matrix: ${\mathbf{W}_{\left\langle t\right\rangle}}$;

17 Generate uncertainty matrix $\mathbf{U}_{\left\langle
t\right\rangle}\leftarrow\mathbf{M}_{\left\langle
t\right\rangle}\odot{\mathbf{W}_{\left\langle t\right\rangle}}$;

18 Compute output: $\mathbf{y}_{\left\langle
t\right\rangle}\leftarrow\mathbf{K}_{\left\langle
t\right\rangle}\cdot\mathfrak{m}(\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle})+\mathbf{a}_{\left\langle
t\right\rangle}\odot\text{act}\left({\mathbf{U}_{\left\langle
t\right\rangle}\cdot\mathfrak{m}\left({\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle t\right\rangle}}\right)}\right)$ ;

19 end for

Output : terminal output:
$\widehat{\mathbf{y}}\leftarrow\mathbf{y}_{\left\langle p\right\rangle}$

.

Algorithm 2 Physics-guided NN Editing

Building on the deep PhN, this section presents the neural network (NN)
editing for embedding and preserving the available physical knowledge, through
the physics-guided link editing and activation editing. Specifically, the link
editing centers around removing and preserving the links according to the
consistency with physical knowledge. Meanwhile, the activation editing
performs the physical-knowledge-preserving computing in the output channels of
PhNs. Thanks to the concurrent link and activation editing, the input/output
of Phy-Taylor can strictly comply with the available physical knowledge. Using
the notation ‘$\boxplus$’ defined in the Table 1, the procedure of physics-
guided NN editing is described in Algorithm 2.

For the edited weight matrix $\mathbf{W}_{\left\langle t\right\rangle}$, if
its entries in the same row are all known model-substructure parameters, the
associated activation should be inactivate. Otherwise, the Phy-Taylor cannot
strictly preserve the available physical knowledge due to the extra nonlinear
mappings induced by the activation functions. This thus motivates the physics-
knowledge preserving computing, i.e., the Line 2 of Algorithm 2. Figure 4
summarizes the flowchart of NN editing in a single PhN layer:

  * •

(i) Given the node-representation vector from Algorithm 1, the original
(fully-connected) weight matrix is edited via link editing to embed assigned
physical knowledge, resulting in $\mathbf{W}_{\left\langle t\right\rangle}$.

  * •

(ii) The edited weight matrix ${\mathbf{W}_{\left\langle t\right\rangle}}$ is
separated into knowledge matrix $\mathbf{K}_{\left\langle t\right\rangle}$ and
uncertainty matrix $\mathbf{U}_{\left\langle t\right\rangle}$, such that
${\mathbf{W}_{\left\langle t\right\rangle}}=\mathbf{K}_{\left\langle
t\right\rangle}+\mathbf{U}_{\left\langle t\right\rangle}$. Specifically, the
$\mathbf{K}_{\left\langle t\right\rangle}$, generated in Lines 2 and 11,
includes all the known model-substructure parameters. While the
$\mathbf{M}_{\left\langle t\right\rangle}$, generated in Lines 2 and 2, is
used to generate uncertainty matrix $\mathbf{U}_{\left\langle t\right\rangle}$
(see Line 2) to include all the unknown model-substructure parameters, through
freezing the known model-substructure parameters of $\mathbf{W}_{\left\langle
t\right\rangle}$ to zeros.

  * •

(iii) The $\mathbf{K}_{\left\langle t\right\rangle}$,
$\mathbf{M}_{\left\langle t\right\rangle}$ and activation-masking vector
$\mathbf{a}_{\left\langle t\right\rangle}$ (generated in Lines 2 and 2) are
used by activation editing for the physical-knowledge-preserving computing of
output in each PhN layer. The function of $\mathbf{a}_{\left\langle
t\right\rangle}$ is to avoid the extra mapping (induced by activation) that
prior physical knowledge does not include.

Figure 4: Flowchart of NN editing in single PhN layer.

The flowchart of NN editing operating in cascade PhN is depicted in Figure 5.
The Lines 2-2 of Algorithm 2 means that $\mathbf{A}=\mathbf{K}_{\left\langle
1\right\rangle}+\mathbf{M}_{\left\langle 1\right\rangle}\odot\mathbf{A}$,
leveraging which and the setting $r_{\left\langle 1\right\rangle}=r$, the
ground-truth model (1) is rewritten as

$\displaystyle\mathbf{y}=(\mathbf{K}_{\left\langle
1\right\rangle}+\mathbf{M}_{\left\langle
1\right\rangle}\odot\mathbf{A})\cdot\mathfrak{m}(\mathbf{x},r)+\mathbf{f}(\mathbf{x})=\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+(\mathbf{M}_{\left\langle
1\right\rangle}\odot\mathbf{A})\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+\mathbf{f}(\mathbf{x}).$ (16)

Figure 5: Implementing NN editing, i.e., Algorithm 2, in the Example 1. (i)
Unknown substructures are formed by the grey links, the known substructures
are formed by the red and blue links. (ii) Cutting black links to avoid
spurious correlations, otherwise, the links introduce the dependence of
$s\left({k+1}\right)$ on mass m, thus contradicting physical knowledge.

We obtain from the Line 2 of Algorithm 2 that the output of the first PhN
layer is

$\displaystyle\mathbf{y}_{\left\langle
1\right\rangle}=\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+\mathbf{a}_{\left\langle
1\right\rangle}\odot\text{act}\left({\mathbf{U}_{\left\langle
1\right\rangle}\cdot{\mathfrak{m}}\left({\mathbf{x},r_{\left\langle
1\right\rangle}}\right)}\right).$ (17)

Recalling that $\mathbf{K}_{\left\langle 1\right\rangle}$ includes all the
known model-substructure parameters of $\mathbf{A}$ while the
$\mathbf{U}_{\left\langle 1\right\rangle}$ includes remainders, we conclude
from (16) and (17) that the available physical knowledge pertaining to the
ground-truth model (1) has been embedded to the first PhN layer. As Figure 5
shows the embedded knowledge shall be passed down to the remaining cascade
PhNs and preserved therein, such that the end-to-end Phy-Taylor model can
strictly with the prior physical knowledge. This knowledge passing is achieved
by the block matrix $\mathbf{K}_{\left\langle p\right\rangle}$ generated in
Line 11, due to which, the output of $t$-th PhN layer satisfies

$\displaystyle[\mathbf{y}_{\left\langle
t\right\rangle}]_{1:\text{len}(\mathbf{y})}=\underbrace{\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})}_{\text{knowledge
passing}}+\underbrace{[\mathbf{a}_{\left\langle
t\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
t\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
t-1\right\rangle},r_{\left\langle
t\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}}_{\text{knowledge
preserving}},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall t\in\\{2,3,\ldots,p\\}.$ (18)

Meanwhile, the $\mathbf{U}_{\left\langle
t\right\rangle}=\mathbf{M}_{\left\langle
t\right\rangle}\odot{\mathbf{W}_{\left\langle t\right\rangle}}$ means the
masking matrix $\mathbf{M}_{\left\langle t\right\rangle}$ generated in Line 2
is to remove the spurious correlations in the cascade PhN, which is depicted
by the cutting link operation in Figure 5.

###  3.3 The Solution to Problem 1: Phy-Taylor

Described in Figure 1, with the guidance of Taylor series, implementing the
physics-guided NN editing in the deep PhN yields the Phy-Taylor. The Phy-
Taylor embeds the available physical knowledge into each PhN, such that its
input/output strictly complies with the physical knowledge, which is formally
stated in the following theorem. The theorem proof appears in Supplementary
Information 9.3.

######  Theorem 2.

Consider the Phy-Taylor described by Figure 1. The input/output (i.e.,
${\mathbf{x}}$/$\widehat{\mathbf{y}}$) of Phy-Taylor strictly complies with
the available knowledge pertaining to the physics model (1) of ground truth,
i.e., if the $[\mathbf{A}]_{i,j}$ is a known model-substructure parameter,
then
$\frac{{\partial{{[\widehat{\mathbf{y}}}]_{i}}}}{{\partial{[\mathfrak{m}}\left({\mathbf{x},r}\right)]_{j}}}\equiv\frac{{\partial{[\mathbf{y}]_{i}}}}{{\partial{[\mathfrak{m}}\left({\mathbf{x},r}\right)]_{j}}}\equiv[\mathbf{A}]_{i,j}$
always holds.

##  4 Phy-Taylor Properties

Figure 6: (a): Two Fully-connected NN layers. (b): A single PhN layer.

Moving forward, this section focuses on the property analysis of Phy-Taylor.

###  4.1 Parameter Quantity Reduction

The Figures 2 and 3 show that the Phy-Augmentation, i.e., the Algorithm 1,
expands the input without involving weight-matrix multiplication. This trait
can be leveraged to significantly reduce the quantity of learning parameters
including weights and bias. For the demonstration, we consider the two network
models in Figure 6 (a) and (b), where the (a) describes a fully-connected two-
layer network, while the (b) describes a single PhN. Observing them, we obtain
that given the same dimensions of input and terminal output, the number of
learning parameters of Figure 6 (a) is $(m+1)n+(n+1)p$ (including $(m+p)n$
weights and $n+p$ bias), while the number of learning parameters of Figure 6
(b) is $(n+1)p$ (including $n\times p$ weights and $p$ bias). The number
difference of learning parameters is thus

$\displaystyle(m+1)n+(n+1)p-(n+1)p=(m+1)n.$ (19)

We note the quantity of reduced parameters (19) is the lower bound of PhN in
Phy-Taylor framework, since it is obtained without considering physics-guided
NN editing for removing and freezing links and bias according to the available
physical knowledge.

###  4.2 Single PhN v.s. Cascade PhN

We next investigate if the space complexity (i.e., the quantity of augmented
monomials) of Phy-Augmentation of a single PhN with a large augmentation order
can be reduced via cascade PhN with relatively small orders. To simplify the
presentation, a single PhN and cascade PhN are respectively represented in the
following equations.

$\displaystyle\widehat{\mathbf{y}}=\mathrm{PhN}({\left.{\mathbf{x}}\in\mathbb{R}^{n}\right|r})\in\mathbb{R}^{m},$
(20) $\displaystyle{\mathbf{x}}\in{{\mathbb{R}}^{n}}\leavevmode\nobreak\
\leavevmode\nobreak\ \longmapsto\leavevmode\nobreak\ \leavevmode\nobreak\
{\mathbf{y}_{\left\langle
1\right\rangle}}=\mathrm{PhN}({\left.{\mathbf{x}}\right|{r_{\left\langle
1\right\rangle}}})\in{{\mathbb{R}}^{{n_{\left\langle
1\right\rangle}}}}\leavevmode\nobreak\ \leavevmode\nobreak\
\longmapsto\leavevmode\nobreak\ \leavevmode\nobreak\
\ldots\leavevmode\nobreak\ \leavevmode\nobreak\
\longmapsto\leavevmode\nobreak\ \leavevmode\nobreak\
{\mathbf{y}_{{\left\langle
d-1\right\rangle}}}=\mathrm{PhN}({\left.{{\mathbf{y}_{{\left\langle
d-2\right\rangle}}}}\right|{r_{{\left\langle
d-1\right\rangle}}}})\in{{\mathbb{R}}^{{n_{{\left\langle
d-1\right\rangle}}}}}$ $\displaystyle\hskip
255.79042pt\longmapsto\leavevmode\nobreak\ \leavevmode\nobreak\
\widehat{\mathbf{y}}=\mathrm{PhN}({\left.{{\mathbf{y}_{{\left\langle
d-1\right\rangle}}}}\right|{r_{\left\langle
d\right\rangle}}})\in\mathbb{R}^{m},$ (21)

where the cascade architecture consists of $d$ PhNs. To guarantee the cascade
PhN (21) and the single PhN (20) have the same monomials, their augmentation
orders are required to satisfy

$\displaystyle\prod\limits_{v=1}^{d}{{r_{{\left\langle
v\right\rangle}}}}=r,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall r_{{\left\langle
v\right\rangle}},\leavevmode\nobreak\ \leavevmode\nobreak\ r\in\mathbb{N}.$
(22)

The space complexity difference of Phy-Augmentation is formally presented in
the following theorem, whose proof is presented in Supplementary Information
9.4.

######  Theorem 3.

Under the condition (22), the space complexity difference between single PhN
(20) and cascade PhN (21), due to Phy-Augmentation, is

$\displaystyle\mathrm{len}(\mathfrak{m}(\mathbf{x},r))-\sum\limits_{p=1}^{d}{\mathrm{len}(\mathfrak{m}(\mathbf{x},r_{{\left\langle
p\right\rangle}}))}=\sum\limits_{s={r_{\left\langle
1\right\rangle}}+1}^{r}{\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)!s!}}}-\sum\limits_{v=1}^{d-1}{\sum\limits_{s=1}^{{r_{{\left\langle
v+1\right\rangle}}}}{\frac{{\left({{n_{\left\langle
v\right\rangle}}+s-1}\right)!}}{{\left({{n_{\left\langle
v\right\rangle}}-1}\right)!s!}}}}+1-d.$ (23)

The Theorem 3 implies that the output dimensions and the augmentation orders
of intermediate PhNs are critical in the significant reduction of space
complexity via cascade PhN. However, an intuitive question arises: Does the
cascade PhN reduce the complexity at the cost of model accuracy? Without loss
of generality, we use the following example to answer the question.

######  Example 2.

For simplicity in explanation, we ignore the bias and consider the scenario
that both the activation and the suppressor are inactive. For the single PhN
(20), we let $\mathbf{x}\in\mathbb{R}^{2}$ and
$\widehat{\mathbf{y}}\in\mathbb{R}$ and $r=4$. Its output is then computed
according to

$\displaystyle\widehat{\mathbf{y}}={w_{1}}{[\mathbf{x}]_{1}}+{w_{2}}{[\mathbf{x}]_{2}}+{w_{3}}[\mathbf{x}]_{1}^{2}+{w_{4}}{[\mathbf{x}]_{1}}{[\mathbf{x}]_{2}}+{w_{5}}[\mathbf{x}]_{2}^{2}+{w_{6}}[\mathbf{x}]_{1}^{3}+{w_{7}}[\mathbf{x}]_{1}^{2}{[\mathbf{x}]_{2}}+{w_{8}}{[\mathbf{x}]_{1}}[\mathbf{x}]_{2}^{2}+{w_{9}}[\mathbf{x}]_{2}^{3}$
$\displaystyle\hskip
137.42685pt+{w_{10}}[\mathbf{x}]_{1}^{4}+{w_{11}}[\mathbf{x}]_{1}^{3}{[\mathbf{x}]_{2}}+{w_{12}}[\mathbf{x}]_{1}^{2}[\mathbf{x}]_{2}^{2}+{w_{13}}{[\mathbf{x}]_{1}}[\mathbf{x}]_{2}^{3}+{w_{14}}[\mathbf{x}]_{2}^{4},$
(24)

For the corresponding cascade PhN (21), we let $r_{{\left\langle
1\right\rangle}}=r_{{\left\langle 2\right\rangle}}=2$, the output dimension of
first PhN is 2 while dimension of terminal output is 1. We thus have

$\displaystyle\widehat{\mathbf{y}}$
$\displaystyle={{\hat{w}}_{6}}{\mathbf{y}_{\left\langle
1\right\rangle}}+{{\hat{w}}_{7}}\mathbf{y}^{2}_{\left\langle
1\right\rangle}={{\hat{w}}_{6}}({{{\tilde{w}}_{1}}{[\mathbf{x}]_{1}}+{{\tilde{w}}_{2}}{[\mathbf{x}]_{2}}+{{\tilde{w}}_{3}}[\mathbf{x}]_{1}^{2}+{{\tilde{w}}_{4}}{[\mathbf{x}]_{1}}{[\mathbf{x}]_{2}}+{{\tilde{w}}_{5}}[\mathbf{x}]_{2}^{2}})$
$\displaystyle\hskip
152.22241pt+{{\hat{w}}_{7}}{\left({{{\tilde{w}}_{1}}{[\mathbf{x}]_{1}}+{{\tilde{w}}_{2}}{[\mathbf{x}]_{2}}+{{\tilde{w}}_{3}}[\mathbf{x}]_{1}^{2}+{{\tilde{w}}_{4}}{[\mathbf{x}]_{1}}{[\mathbf{x}]_{2}}+{{\tilde{w}}_{5}}[\mathbf{x}]_{2}^{2}}\right)^{2}}.$
(25)

Observing (24) and (25), we discover that (i) both the single and cascade
architectures have the same monomials due to the satisfying condition (22),
and (ii) the single PhN layer has 14 weight parameters (i.e., the ${w_{1}}$,
${w_{2}}$, $\dots$, ${w_{14}}$ in Equation (24)), while the cascade layers
have only 7 weight parameters (i.e., the ${\hat{w}_{1}}$, $\ldots$, and
$\hat{w}_{7}$ in Equation (25)) in total. Intuitively, we can conclude that if
the reduced weights are associated with the links that contradict with
physical knowledge, the cascade PhN can further increase model accuracy,
otherwise, it can reduce the space complexity at the cost of model accuracy.

##  5 Extension: Self-Correcting Phy-Taylor

Figure 7: Self-correcting Phy-Taylor Architecture: $\mathbf{u}(k)$ denotes the
vector of real-time decisions, ${\bf{s}}(u(k))$ denotes the vector of real-
time safety metrics, $\bf{c}$ is the vector of safety bounds.

The recent incidents due to deployment of DNNs overshadow the revolutionizing
potential of AI in the physical engineering systems [36, 1], especially the
safety-critical systems, whose unintended behavior results in death or serious
injury to people, severe damage to equipment or environments [37, 38]. The
safe control and planning is a fundamental solution for enhancing safety
assurance of the AI-assisted physical systems often operating in environments
where time and safety are critical, such as the airplanes, medical drones and
autonomous vehicles. To comply with safety constraints in face of potential
conflicts from control objectives, the framework of control barrier function
(CBF) has been proposed for the computation of real-time safety-critical
control commands [39, 40]. The CBFs however use only current state information
without prediction, whose control policy is thus greedy and challenging for
proactive safe control. It is well known that model predictive control (MPC)
yields a less greedy safe control policy, since it takes future state
information into account [41, 42]. Motivated by the observations, MPC with
incorporation of CBF, i.e., MPC-CBF, was proposed [43]. Due to the nonlinear
dynamics, the MPC-CBF however faces a dilemma of prediction horizon and
computation time of safe control commands, which induces considerable feedback
delays and thus leads to failures in the time- and safety-critical operating
environments. To address the dilemma, we propose the self-correcting Phy-
Taylor, whose architecture is shown in Figure 7. Its one mission is learning
the safety relationship between the real-time decisions and the safety
metrics, with consideration of future information:

$\displaystyle\mathbf{s}(\mathbf{x}(k),\mathbf{u}(k),\tau)=\sum\limits_{t=k}^{k+\tau-1}{\tilde{\mathbf{f}}(\mathbf{x}(t),\mathbf{u}(t))},$
(26)

where $\tilde{\mathbf{f}}(\mathbf{x}(t),\mathbf{u}(t))$ is the predefined
vector of safety metrics at time $t$, and $\tau$ denotes the horizon of future
information of safety metrics.

Inside the self-correcting Phy-Taylor, the learned safety relationship for
approximating (26) will first be subject to the off-line verification of
available physical knowledge, based on which, the necessary revisions can be
needed. According to the off-line verified and revised (if needed) safety
relationship, the correcting of real-time decision $\mathbf{u}(k)$ will be
triggered if any safety metric
$[\mathbf{s}(\mathbf{x}(k),\mathbf{u}(k),\tau)]_{i}$,
$i\in\\{1,2,\ldots,h\\}$, exceeds (or leaves) the preset safety bounds (or
safety envelopes). However, the current learned formula corresponding to (26)
is not ready (if not impossible) for delivering the procedure, owning to the
complicated dependence of $[\mathbf{s}(\mathbf{x}(k),\mathbf{u}(k),\tau)]_{i}$
on both system state $\mathbf{x}(k)$ and decision $\mathbf{u}(k)$. To address
this problem, as shown in Figure 7, we decouple the real-time decisions from
the real-time system states. Specifically,

  * •

Given the real-time system state $\mathbf{x}(k)$ as the origin input, the
first Phy-Taylor outputs the real-time decision $\mathbf{u}(k)$, which is
motivated by the fact that the state-feedback control is used most commonly in
physical engineering systems [44]. In other words, the computation of raw
$\mathbf{u}(k)$ directly depends on real-time system state $\mathbf{x}(k)$.

  * •

Given the real-time decision $\mathbf{u}(k)$ (i.e., the output of the first
Phy-Taylor) as the input of the second Phy-Taylor, the terminal output is the
real-time safety metric $\mathbf{s}(\mathbf{u}(k))$, which is motivated by the
fact that the decision $\mathbf{u}(k)$ manipulates system state. In other
words, the safety metric $\mathbf{s}(\mathbf{u}(k))$ directly depends on
decision $\mathbf{u}(k)$ and indirectly depends on system state
$\mathbf{x}(k)$.

  * •

The two Phy-Taylors are trained simultaneously according to the training loss
function:

$\displaystyle\mathcal{L}=\alpha\left\|{{\bf{s}}(\mathbf{u}(k))-{\bf{s}}({\bf{x}}(k),{\bf{u}}(k),\tau)}\right\|+\beta\left\|{\breve{\bf{u}}(k)-\bf{u}(k)}\right\|,$
(27)

where the ${\bf{s}}({\bf{x}}(k),{\bf{u}}(k),\tau)$ given in Equation (26) is
ground truth of safety-metric vector, the $\breve{\bf{u}}(k)$ is ground truth
of decision vector, the $\alpha$ and $\beta$ are hyperparameters. The two
cascade Phy-Taylors thus depend on each other.

  * •

To render the learned safety relationship ${\bf{s}}(\mathbf{u}(k))$ tractable,
the activation and compressor inside the second Phy-Taylor are inactive, such
that the ${\bf{s}}(\mathbf{u}(k))$ is expressed in the form of Taylor series.

Given the verified and revised (if needed) relationship, the self-correcting
procedure will be triggered (if a safety metric exceeds the safety bound
$\mathbf{c}$) for correcting decision according to

$\displaystyle\mathbf{u}(k)\leftarrow\mathop{\arg\min}\limits_{\widetilde{\mathbf{u}}(k)\in{\mathbb{R}^{m}}}\left\\{{\left.{\left\|{\widetilde{\mathbf{u}}(k)-\mathbf{u}(k)}\right\|}\right|[\bf{s}(\widetilde{\mathbf{u}}(k)]_{i}<[\mathbf{c}]_{i},\leavevmode\nobreak\
\leavevmode\nobreak\
i\in\\{1,2,\ldots,\text{len}({\bf{s}}(\mathbf{u}(k)))\\}}\right\\}.$ (28)

We note the self-correcting mechanism and the safety revision of relationship
between ${\bf{s}}(\mathbf{u}(k))$ and $\mathbf{u}(k)$ for delivering (28) vary
with safety problems and physical systems. An example in this paper is the
safe control of autonomous vehicles: Algorithm 3 in the Experiments section.

##  6 Experiments

The demonstrations are performed on three different systems with different
degrees of available physical knowledge: autonomous vehicles, coupled
pendulums, and a subsystem of US Illinois climate.

###  6.1 Autonomous Vehicles

Figure 8: Vehicle’s driving environment.

The first experiment performs the demonstration of two functions: (i) the
learning of vehicle’s conjunctive lateral and longitudinal dynamics via Phy-
Taylor, and (ii) the safe velocity regulation via self-correcting Phy-Taylor.
The vehicle’s driving environment is shown in Figure 8, operating in the
AutoRally platform [45].

####  6.1.1 Vehicle Dynamics Learning

We first leverage our available physical knowledge to identify the known
model-substructure parameters for the physics-guided NN editing. According to
the Newton’s second law for motion along longitudinal and lateral axes [28],
we have the following governing equations:

$\displaystyle\bar{m}\ddot{\mathrm{p}}={F_{\mathrm{p}f}}+{F_{\mathrm{p}r}}-{F_{\text{aero}}}-{R_{\mathrm{p}f}}-{R_{\mathrm{p}r}},\quad\quad\bar{m}(\ddot{\mathrm{y}}+\dot{\psi}{v_{\mathrm{p}}})={F_{\mathrm{y}f}}+{F_{\mathrm{y}r}},$
(29)

where $\mathrm{p}$ is the longitudinal position, $\mathrm{y}$ is the lateral
position, $\psi$ is the vehicle yaw angle, $\bar{m}$ is the vehicle mass,
$v_{\mathrm{p}}\triangleq\dot{\mathrm{p}}$ is the longitudinal velocity,
${F_{\mathrm{p}f}}$ and ${F_{\mathrm{p}r}}$ denote the longitudinal tire force
at the front and rear tires, respectively, ${R_{\mathrm{p}f}}$ and
${R_{\mathrm{p}r}}$ denote the rolling resistance at the front and rear tires,
respectively, ${F_{\text{aero}}}$ represents the longitudinal aerodynamic drag
force, ${F_{\mathrm{y}f}}$ and ${F_{\mathrm{y}r}}$ are the lateral tire forces
of the front and rear wheels, respectively. With the notations of lateral
velocity $v_{\mathrm{y}}\triangleq\dot{\mathrm{y}}$ and yaw velocity
$v_{\psi}\triangleq\dot{\psi}$, the following state space model is derived
from the force balance equation (29) in the literature [28].

$\displaystyle\frac{\mathrm{d}}{{\mathrm{d}t}}\left[\begin{gathered}\mathrm{p}\hfill\\\
\mathrm{y}\hfill\\\ \psi\hfill\\\ {v_{\mathrm{p}}}\hfill\\\
{v_{\mathrm{y}}}\hfill\\\ {v_{\psi}}\hfill\\\
\end{gathered}\right]=\left[{\begin{array}[]{*{20}{c}}0&0&0&1&0&0\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ 0&0&0&*&0&0\\\ 0&0&0&0&*&*\\\
0&0&0&0&*&*\end{array}}\right]\underbrace{\left[\begin{gathered}\mathrm{p}\hfill\\\
\mathrm{y}\hfill\\\ \psi\hfill\\\ {v_{\mathrm{p}}}\hfill\\\
{v_{\mathrm{y}}}\hfill\\\ {v_{\psi}}\hfill\\\
\end{gathered}\right]}_{\triangleq{\mathbf{x}}}+\left[\begin{gathered}0\hfill\\\
0\hfill\\\ 0\hfill\\\ *\hfill\\\ 0\hfill\\\ 0\hfill\\\
\end{gathered}\right]\theta+\left[\begin{gathered}0\hfill\\\ 0\hfill\\\
0\hfill\\\ 0\hfill\\\ *\hfill\\\ *\hfill\\\ \end{gathered}\right]\delta,$ (64)

where ‘*’ can represent a state-dependent or time-dependent function or mixed
of them or just a scalar, but is unknown to us, and $\theta$ and $\delta$
denote throttle and steering, respectively. Given the practical physical
knowledge that the throttle computation depends on the longitudinal velocity
and position only, while the dependencies of steering are unknown, the state
space model (64) updates with

$\displaystyle\dot{{\mathbf{x}}}=\left[{\begin{array}[]{*{20}{c}}0&0&0&1&0&0\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ &0&0&*&0&0\\\ &*&*&*&*&*\\\
&*&*&*&*&*\end{array}}\right]{\mathbf{x}}.$ (71)

The sampling technique, with sampling period denoted by $T$, converts the
continuous-time state-space model above to the following discrete-time one:

$\displaystyle{\mathbf{x}}\left({k+1}\right)=\left[{\begin{array}[]{*{20}{c}}1&0&0&T&0&0\\\
0&1&0&0&T&0\\\ 0&0&1&0&0&T\\\ &0&0&*&0&0\\\ &*&*&*&*&*\\\
&*&*&*&*&*\end{array}}\right]{\mathbf{x}}\left(k\right).$ (78)

Figure 9: (a): Phy-Taylor${}_{\text{large order}}$ has a PhN with large order
$r=4$. (b): Phy-Taylor${}_{\text{small order}}$ has cascading PhNs with
relatively small orders satisfying $r_{\left\langle 1\right\rangle}\cdot
r_{\left\langle 2\right\rangle}=2\cdot 2=4=r$.

We first consider two Phy-Taylor models, named ‘Phy-Taylor${}_{\text{large
order}}$’ and ‘Phy-Taylor${}_{\text{small order}}$’, which can embed the
available knowledge (i.e., the known parameters included in system matrix of
model (78)). Their architectures are shown in Figure 9 (a) and (b). The Phy-
Taylor${}_{\text{large order}}$ has one PhN with a large augmentation order
while the Phy-Taylor${}_{\text{small order}}$ has two cascade PhN layers with
two relatively small augmentation orders. Meanwhile, the three orders satisfy
the condition (22) for having the same monomials of Taylor series. We also
consider the corresponding models without NN editing (i.e., without physical
knowledge embedding), which degrades the Phy-Taylor to the deep PhN (DPhN).
The two DPhN models are named ‘DPhN${}_{\text{large order}}$’ and
‘DPhN${}_{\text{small order}}$’. The final model we considered is the seminal
Deep Koopman [23], following the same model configurations therein. The
configurations of five models are summarized in Table 2.

Table 2: Model Configurations

| Layer 1 | Layer 2 | Layer 3 |  |
---|---|---|---|---|---
Model ID | #weights | #bias | #Weights | #bias | #weights | #bias | #parameter sum |  prediction error $e$
DPhN${}_{\text{large order}}$ | $3135$ | $15$ | $90$ | $6$ | $-$ | $-$ | $3246$ | nan
DPhN${}_{\text{small order}}$ | $270$ | $10$ | $520$ | $8$ | $48$ | $6$ | $862$ | $45.57277$
Phy-Taylor${}_{\text{large order}}$ | $2313$ | $12$ | $30$ | $3$ | $-$ | $-$ | $2358$ | $0.047758$
Phy-Taylor${}_{\text{small order}}$ | $167$ | $7$ | $265$ | $5$ | $18$ | $3$ | $465$ | $0.003605$

| Encoder | Decoder | Auxiliary Networks |  |
---|---|---|---|---|---
Model ID | #weights | #bias | #weights | #bias | #weights | #bias | #parameter sum | prediction error
Deep Koopman | $2040$ | $176$ | $2040$ | $176$ | $19920$ | $486$ | $24838$ | $0.232236$

Figure 10: Training and Testing. (a)-(c): The trajectories of averaged
training loss (5 random seeds) of different models described in Table 2.
(i)-(vi): Ground truth and predicted trajectories via trained models.

The trajectories of training loss are presented in Figure 10 (a)–(c). The
(training loss, validation loss) of trained DPhN${}_{\text{large order}}$,
DPhN${}_{\text{small order}}$, Phy-Taylor${}_{\text{large order}}$ and Phy-
Taylor${}_{\text{small order}}$ are (0.00389, 0.00375), (0.000344, 0.000351),
(0.000222, 0.000238) and (0.000915, 0.000916), respectively. To perform the
testing, we consider the long-horizon prediction of system trajectories, given
the same initial conditions. The prediction error is measured by the mean
squared error:
$e=\frac{1}{\kappa}\sum\limits_{t=k+1}^{k+\kappa}{\frac{1}{6}}\left\|{\underline{\bf{x}}\left(t\right)-{\bf{x}}\left(t\right)}\right\|\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\
\underline{\bf{x}}\left(k\right)={\bf{x}}\left(k\right)$, where
$\underline{\bf{x}}\left(t\right)$ is the prediction of ground truth
${\bf{x}}\left(t\right)$ at time $t$. The prediction errors over the horizon
$\kappa=300$ and initial time $k=950$ are summarized in Table 2. The ground-
truth trajectories and predicted ones from Deep Koopman and the Phy-Taylor
models are presented in Figure 10 (i)–(vi). Observing from Table 2 and Figure
10, we can conclude:

  * •

Figure 10 (i)–(vi) and Figure 10 (a)–(b): the physics-guided NN editing can
significantly accelerate model training, reduce validation and training loss
as well as improve model accuracy (viewed from long-horizon prediction of
trajectory).

  * •

Figure 10 (i)–(vi) and Figure 10 (c): with physics-guided NN editing, the
cascade PhN with small augmentation orders can further significantly reduce
training loss, and increase model accuracy. This can be due to the further
removed spurious correlations or NN links contradicting with physics law, via
the cascade architecture.

  * •

Figure 10 (i)–(vi) and Table 2: compared with the fully-connected DPhN models,
i.e., DPhN${}_{\text{large order}}$ and DPhN${}_{\text{small order}}$, the
seminal Deep Koopman strikingly increases model accuracy, viewed from the
perspective of long-horizon prediction of trajectory. Conversely, compared to
Deep Koopman, the Phy-Taylor models (both Phy-Taylor${}_{\text{large order}}$
and Phy-Taylor${}_{\text{small order}}$) notably reduce the model learning
parameters (weights and bias) and further remarkably increase model accuracy
simultaneously.

####  6.1.2 Self-Correcting Phy-Taylor

This experiment demonstrates the effectiveness of self-correcting Phy-Taylor
in guaranteeing vehicle’s safe driving in the environment shown in Figure 8.
The architecture of self-correcting Phy-Taylor is presented in Figure 11. Its
real-time input vector is ${\mathbf{x}}(k)=[w(k);\leavevmode\nobreak\
\mathrm{p}(k);\leavevmode\nobreak\ \mathrm{y}(k);\leavevmode\nobreak\
\psi(k);\leavevmode\nobreak\ v_{\mathrm{p}}(k);\leavevmode\nobreak\
v_{\mathrm{y}}(k);v_{\psi}(k)]$, where $w(k)$ is the average of four wheels’
velocities. The mediate output
$\mathbf{u}(k)=\left[\theta(k);\leavevmode\nobreak\ \gamma(k)\right]$ denotes
the vector of control commands, where
$\theta(k)\in[-0.156,\leavevmode\nobreak\ 0.156]$ is the throttle command and
$\gamma(k)\in[-0.6,\leavevmode\nobreak\ 0.6]$ is the steering command. The
considered safety-metric vector in Equation (26) is

$\displaystyle\mathbf{s}({\mathbf{x}}(k),\mathbf{u}(k),\tau)=\sum\limits_{t=k+1}^{k+\tau}{\left[{{{\left({{v_{\mathrm{p}}}(t)-\mathrm{v}}\right)}^{2}};\leavevmode\nobreak\
\leavevmode\nobreak\ {{\left({{v_{\mathrm{p}}}(t)-r\cdot
w(k)}\right)}^{2}}}\right]}\in{\mathbb{R}^{2}},$ (79)

Figure 11: Self-correcting Phy-Taylor for safe control of autonomous vehicle.

where $\mathrm{v}$ and $r$ denote the reference of longitudinal velocity and
the wheel radius, respectively. The safety metric (79) together with the
driving scenario in Figure 8 indicate the objective of safe control command is
to simultaneously steer vehicle’s longitudinal velocity to reference
$\mathrm{v}$, and constrain the slip (i.e., $(v_{\mathrm{x}}(t)-r\cdot
w(k))^{2}$) to prevent slipping and sliding. The hyperparameters in the
training loss function (27) are set to $\alpha=\beta=1$.

The testing results of trained model are presented in Figure 12 (a)–(d) (blue
curves), which in conjunction with the ground truth (orange curves)
demonstrate the trustfulness of trained model. We next output the learned
safety relationships for off-line verification and necessary revision:

$\displaystyle[{\bf{s}}({\bf{u}}(k))]_{1}$
$\displaystyle=0.00111007+{\left[\\!\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\\!\right]^{\top}}\left[{\begin{array}[]{*{20}{c}}{-0.04581441}&{0.00100625}\\\
{0.00100625}&{0.00342825}\end{array}}\right]\left[\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\right],$ (80g) $\displaystyle[{\bf{s}}({\bf{u}}(k))]_{2}$
$\displaystyle=0.14376973-{\left[\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\right]^{\top}}\left[{\begin{array}[]{*{20}{c}}{6.06750536}&{0.02701398}\\\
{0.02701398}&{0.00601609}\end{array}}\right]\left[\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\right].$ (80n)

The safety metrics (79) of ground truth are always non-negative. We thus need
to verify that given the ranges of control commands (i.e.,
$\theta(k)\in[-0.156,0.156]$, $\gamma(k)\in[-0.6,0.6]$, $\forall
k\in\mathbb{N}$), if both $[{\bf{s}}({\bf{u}}(k))]_{1}$ and
$[{\bf{s}}({\bf{u}}(k))]_{2}$ in Equation (80) can always be non-negative. If
a violation occurs, we will make revisions to the relationships. We can verify
from (80) that the non-negativity constraint does not always hold, such as
$\theta(k)=0.156$ and $\gamma(k)=0$. Therefore, the revision of safety
relationship is needed before working on the self-correcting procedure. The
regulated safety relationships (revisions are highlighted in red color) are
presented below.

$\displaystyle[{\bf{s}}({\bf{u}}(k))]_{1}$
$\displaystyle=\underbrace{0.00{\color[rgb]{1.00,0.00,0.00}\definecolor[named]{pgfstrokecolor}{rgb}{1.00,0.00,0.00}02}1007}_{{[\mathbf{b}]_{1}}}+{\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right]^{\top}}\underbrace{\left[{\begin{array}[]{*{20}{c}}{{\color[rgb]{1.00,0.00,0.00}\definecolor[named]{pgfstrokecolor}{rgb}{1.00,0.00,0.00}0.0018}1441}&{0.00100625}\\\
{0.00100625}&{0.00342825}\end{array}}\right]}_{{\triangleq{\bf{P}}_{1}}}\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right],$ (81g)
$\displaystyle[{\bf{s}}({\bf{u}}(k))]_{2}$
$\displaystyle=\underbrace{0.14376973}_{{[\mathbf{b}]_{2}}}-{\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right]^{\top}}\underbrace{\left[{\begin{array}[]{*{20}{c}}{{\color[rgb]{1.00,0.00,0.00}\definecolor[named]{pgfstrokecolor}{rgb}{1.00,0.00,0.00}5.90769724}}&{{\color[rgb]{1.00,0.00,0.00}\definecolor[named]{pgfstrokecolor}{rgb}{1.00,0.00,0.00}0.01201398}}\\\
{{\color[rgb]{1.00,0.00,0.00}\definecolor[named]{pgfstrokecolor}{rgb}{1.00,0.00,0.00}0.01201398}}&{0.00601609}\end{array}}\right]}_{\triangleq{\bf{P}}_{2}}\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right],$ (81n)

which satisfy $\left[{\bf{s}}({\bf{u}}(k))\right]_{1}\geq 0$ and
$\left[{\bf{s}}({\bf{u}}(k))\right]_{2}\geq 0$, for any
$\theta(k)\in[-0.156,0.156]$ and $\gamma(k)\in[-0.6,0.6]$, and can be
demonstrated by the green curves in Figure 12 (c) and (d).

Input: Real-time control-command vector
${\bf{u}}(k)=[\theta(k);\leavevmode\nobreak\ \gamma(k)]$, safety bounds
$[\mathbf{c}]_{1}$ and $[\mathbf{c}]_{2}$, and learned matrices
$\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ and bias $[\mathbf{b}]_{1}$ and
$[\mathbf{b}]_{2}$ defined in Equation (81).

1 Update original safety relationship with off-line verified and revised one:
${\bf{s}}(\bf{u}(k))\leftarrow\eqref{rchoo}$;

2 if  _$[{\bf{s}}({\bf{u}}(k))]_{1} >[\mathbf{c}]_{1}$ or
$[{\bf{s}}({\bf{u}}(k))]_{2}>[\mathbf{c}]_{2}$ _ then

3 if  _$[{\bf{s}}({\bf{u}}(k))]_{i}\geq[\mathbf{c}]_{i}$ , $i\in\\{1,2\\}$ _
then

4 Update safety metric:
$[\widehat{\mathbf{c}}]_{i}\leftarrow[\mathbf{c}]_{i},i\in\\{1,2\\}$;

5 else

6 Maintain safety metric:
$[\widehat{\mathbf{c}}]_{i}\leftarrow[{\bf{s}}({\bf{u}}(k))]_{i},i\in\\{1,2\\}$;

7 end if

8 Compute orthogonal matrix $\mathbf{P}_{1}$ and eigenvalues $\lambda_{1}$ and
$\lambda_{1}$ according to (88);

9 Compute matrix:
$\mathbf{S}\leftarrow\mathbf{Q}_{1}\cdot\mathbf{P}_{2}\cdot\mathbf{Q}_{1}$;

10 Compute $\widehat{\theta}(k)$ and $\widehat{\gamma}(k)$ according to (99);

11 Correct real-time control commands:
$\displaystyle\theta(k)\leftarrow\mathop{\arg\min}\limits_{\left\\{{\widehat{\theta}(k),-\widehat{\theta}(k)}\right\\}}\left\\{{|{\theta(k)-\widehat{\theta}(k)}|,|{\theta(k)+\widehat{\theta}(k)}|}\right\\},\leavevmode\nobreak\
\leavevmode\nobreak\
\gamma(k)\leftarrow\mathop{\arg\min}\limits_{\left\\{{\widehat{\gamma}(k),-\widehat{\gamma}(k)}\right\\}}\left\\{{|{\gamma(k)-\widehat{\gamma}(k)}|,|{\gamma(k)+\widehat{\gamma}(k)}|}\right\\}.$

12else

13 Maintain real-time control commands: $\theta(k)\leftarrow\theta(k)$ and
$\gamma(k)\leftarrow\gamma(k)$.

14 end if

Algorithm 3 Self-Correcting Procedure for Safe Control Commands

Figure 12: (a)–(b): Ground truth and inference of control commands. (c)–(d):
Ground truth, original inference and inference (according to the revised
safety relationships) of safety metrics (tracking error
${({{v_{\mathrm{p}}}(t)-\mathrm{v}})}^{2}$ and wheel slip
$(v_{\mathrm{p}}(t)-r\cdot w(k))^{2}$). (e): Safety metric: tracking error.
(f): Safety metric: wheel slip. (g): Average of wheel velocities. (h):
Vehicle’s driving curve (phase plot).

With the revised safety relationships (81) at hand, we are ready to develop
the self-correcting procedure. Considering the two matrices ${\bf{P}}_{1}$ and
${\bf{P}}_{2}$ defined in Equation (81) are symmetric, we have

$\displaystyle{\bf{P}}_{1}$
$\displaystyle=\underbrace{\left[{\begin{array}[]{*{20}{c}}{-0.934}&{-0.3572}\\\
{-0.3572}&{0.934}\end{array}}\right]}_{\triangleq{\bf{Q}}_{1}}\left[{\begin{array}[]{*{20}{c}}{\underbrace{0.0008}_{\triangleq{\lambda}_{1}}}&0\\\
0&{\underbrace{0.0038}_{\triangleq{\lambda}_{2}}}\end{array}}\right]\underbrace{\left[{\begin{array}[]{*{20}{c}}{-0.934}&{-0.3572}\\\
{-0.3572}&{0.934}\end{array}}\right]}_{={\bf{Q}}^{\top}_{1}},$ (88)
$\displaystyle{\bf{P}}_{2}$
$\displaystyle={\bf{Q}}_{1}\cdot{\bf{Q}}_{1}\cdot{\bf{P}}_{2}\cdot{\bf{Q}}_{1}\cdot{\bf{Q}}_{1},$
(89)

based on which, we further define:

$\displaystyle\widehat{\mathbf{u}}(k)\triangleq{\bf{Q}}_{1}\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\bf{S}\triangleq{\bf{Q}}_{1}\cdot{\bf{P}}_{2}\cdot{\bf{Q}}_{1}=\left[{\begin{array}[]{*{20}{c}}{{s_{11}}}&{{s_{12}}}\\\
{{s_{12}}}&{{s_{22}}}\end{array}}\right].$ (94)

We let $[\widehat{\mathbf{c}}]_{1}$ and $[\widehat{\mathbf{c}}]_{2}$ denote
the two assigned safety metrics. According to the derivations appearing in
Supplementary Information 9.5, the control commands included in the safety
formulas $[{\bf{s}}({\bf{u}}(k))]_{1}=[\widehat{\mathbf{c}}]_{1}$ and
$\left[{\bf{s}}({\bf{u}}(k))\right]_{2}=[\widehat{\mathbf{c}}]_{2}$ are
obtained as

$\displaystyle\left[\begin{array}[]{l}\pm\widehat{\theta}(k)\\\
\pm\widehat{\gamma}(k)\end{array}\right]\triangleq{\mathbf{Q}_{1}}\left[\begin{array}[]{l}\pm\sqrt{\frac{{{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}}{{{\lambda_{1}}}}-\frac{{{\lambda_{2}}}}{{{\lambda_{1}}}}\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}}\\\
\pm\sqrt{\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}}\end{array}\right],$
(99)

where

$\displaystyle{\varpi_{1}}$
$\displaystyle\triangleq{\left({\frac{{{\lambda_{2}}}}{{{\lambda_{1}}}}}\right)^{2}}s_{11}^{2}+s_{22}^{2}+\frac{{\left({4s_{12}^{2}-2{s_{11}}{s_{22}}}\right){\lambda_{2}}}}{{{\lambda_{1}}}},$
(100) $\displaystyle{\varpi_{2}}$
$\displaystyle\triangleq\frac{{2({{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}){s_{11}}{s_{22}}-4({{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}})s_{12}^{2}+2{\lambda_{2}}({{[\mathbf{b}]_{2}}-{[\widehat{\mathbf{c}}]_{2}}}){s_{11}}}}{{{\lambda_{1}}}}$
$\displaystyle\hskip
213.39566pt-\frac{{2({{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}){\lambda_{2}}s_{11}^{2}}}{{\lambda_{1}^{2}}}-2({{[\mathbf{b}]_{2}}-{[\widehat{\mathbf{c}}]_{2}}}){s_{22}},$
(101) $\displaystyle{\varpi_{3}}$
$\displaystyle\triangleq{\left({{[\mathbf{b}]_{2}}-{[\widehat{\mathbf{c}}]_{2}}}\right)^{2}}+{\left({\frac{{{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}}{{{\lambda_{1}}}}}\right)^{2}}s_{11}^{2}-\frac{{2\left({{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}\right)\left({{[\mathbf{b}]_{2}}-{[\widehat{\mathbf{c}}]_{2}}}\right){s_{11}}}}{{{\lambda_{1}}}}.$
(102)

The solution (99) has paved the way to delivering the self-correcting
procedure: Algorithm 3. The algorithm can be summarized as if the real-time
safety metric $[{\bf{s}}({\bf{u}}(k))]_{1}$ or $[{\bf{s}}({\bf{u}}(k))]_{2}$
is larger than the corresponding safety bound $[\mathbf{c}]_{1}$ or
$[\mathbf{c}]_{2}$, the real-time safety metric will be updated with the
corresponding safety bound (indicated by Line 3 of Algorithm 3). The corrected
control commands are then computed according to (99) (see Lines 3-3). The
solutions however are not unique. To address the problem, the Line 11 of
Algorithm 3 picks up the control commands that are most close to current ones.

Under the control of Phy-Taylor, with and without the self-correcting
procedure, the system performances are presented in Figure 12 (e)–(h). We can
observe from the figure that the self-correcting procedure can significantly
enhance the safe assurance of velocity regulation. The demonstration video of
implementing the self-correcting Phy-Taylor in AutoRally is available at
https://ymao578.github.io/pubs/taylorcontrol2.mp4.

###  6.2 Coupled Pendulums

Figure 13: Mechanical analog of three coupled pendulums.

The second experiment demonstrates the proposed approach to learn the dynamics
of three coupled pendulums, whose mechanical analog is shown in Figure 13.
Each pendulum is characterized by its phase angle $\theta_{i}$ and velocity
$\dot{\theta}_{i}$. The physical topology is described by
$a_{12}=a_{21}=a_{23}=a_{32}=1$ and $a_{13}=a_{31}=0$. We let
${\upsilon_{i}}\triangleq{\dot{\theta}_{i}}$. Referring to Figure 13, the
available physical knowledge for NN editing are summarized in Table 3.

Table 3: Models with Different Degrees of Embedded Physical Knowledge | Available Physical Knowledge |  |
---|---|---|---
Model ID |  Physics Law: $\upsilon=\dot{\theta}$ |  Sampling Period: $T$ |  Coupling Topology: $1\leftrightsquigarrow 2\leftrightsquigarrow 3$ |  Force Dependency | Training Loss |  Out-of-Distribution: Prediction Error $e$
Phy-Taylor-1 | $\surd$ | $\surd$ | $\surd$ | $\surd$ | $1.75146\\!\cdot\\!10^{-6}$ | $0.06486147$
Phy-Taylor-2 | $\surd$ | $\times$ | $\surd$ | $\times$ | $2.42426\\!\cdot\\!10^{-6}$ | $1.46647886$
FDNN | $\times$ | $\times$ | $\times$ | $\times$ | $6.63564\cdot 10^{-7}$ | $3.65772883$

Figure 14: (i): Phy-Taylor architecture. (ii): FDNN architecture. (a)–(f):
Ground truth and predicted trajectories.

We consider the dynamics learning via Phy-Taylor and fully-connected DNN
(FDNN), whose architectures are shown in Figure 14 (a) and (b). The inputs of
both networks are identical as $\mathbf{x}(k)$ $=$ [${\theta_{1}}({k})$;
${\theta_{2}}({k})$; ${\theta_{3}}({k})$; ${\upsilon_{1}}({k})$;
${\upsilon_{2}}({k})$; ${\upsilon_{3}}({k})$]. For a fair comparison, we let
each layer of FDNN has the same number of nodes as the Phy-Taylor (after
augmentation), which suggests the network configuration in Figure 14 (b).

We aim to demonstrate the robustness of Phy-Taylor owning to NN editing. To
achieve this, we consider the testing data that is out-of-distribution.
Specifically, the initial conditions of phase angles for generating training
data are ${\theta_{1}}(0),{\theta_{2}}(0),{\theta_{3}}(0)\in[-1,1]$, while the
initial conditions for generating testing data are outside the range:
${\theta_{1}}(0),{\theta_{2}}(0),{\theta_{3}}(0)\in[-1.5,-1)$. The training
loss (mean square error) of the considered models with different degrees of
embedded physical knowledge are summarized in Table 3. To review the
robustness of trained models, we consider the prediction of long-horizon
trajectory, given the same initial input. The prediction error is measured by
$e=\frac{1}{6}\sum\limits_{i=1}^{3}{\frac{1}{\tau}\left({\sum\limits_{t=k+1}^{k+\tau}{\left({{{\widehat{\theta}}_{i}}(t)-{\theta_{i}}(t)}\right)^{2}}+\sum\limits_{t=k+1}^{k+\tau}{\left({{{\widehat{\upsilon}}_{i}}(t)-{\upsilon_{i}}(t)}\right)^{2}}}\right)}$,
where ${{\widehat{\theta}}_{i}}(t)$ and ${{\widehat{\upsilon}}_{i}}(t)$ denote
the predicted angle and angular speed of $i$-th pendulum at time $t$,
respectively. The prediction errors are summarized in Table 3, which, in
conjunction with ground-truth and predicted trajectories in Figure 14 (c)–(h),
demonstrate that (i) the NN editing can significantly enhance the model
robustness, and (ii) the more embedded physical knowledge can lead to stronger
robustness, viewed from the perspective of long-horizon (out-of-distribution)
predictions of trajectories.

###  6.3 US Illinois Climate

In final example, we consider a climate system without any available physical
knowledge for NN editing, which degrades Phy-Taylor to deep PhN. The dataset
is the hourly climate normals in Illinois state, including five stations111
NOAA: https://www.ncdc.noaa.gov/cdo-
web/datasets/NORMAL_HLY/locations/FIPS:17/detail. The period of record is
01/01/2010–12/31/2010. The locations of considered stations are shown in
Figure 15, where $S_{1}$-$S_{4}$ denote stations GHCND:USW00094846,
GHCND:USW00094822, GHCND:USW00014842 and GHCND:USW00093822, respectively. The
input of deep PhN is

$\displaystyle\mathbf{x}(k)=\left[{{x_{1}}(k);\leavevmode\nobreak\
{x_{2}}(k);\leavevmode\nobreak\ {x_{3}}(k);\leavevmode\nobreak\
{x_{4}}(k);\leavevmode\nobreak\ {{\dot{x}}_{1}}(k);\leavevmode\nobreak\
{{\dot{x}}_{2}}(k);\leavevmode\nobreak\
}\right.{{\dot{x}}_{3}}(k);\leavevmode\nobreak\
{{\dot{x}}_{4}}(k);\leavevmode\nobreak\
{{\ddot{x}}_{1}}(k);\leavevmode\nobreak\
{{\ddot{x}}_{2}}(k);\leavevmode\nobreak\
{{\ddot{x}}_{3}}(k);\leavevmode\nobreak\ \left.{{{\ddot{x}}_{4}}(k)}\right],$

Figure 15: Station locations.

where
${{\dot{x}}_{i}}\left(k\right)={{x}_{i}}\left(k\right)-{{x}_{i}}\left(k-1\right)$
and
${{\ddot{x}}_{i}}\left(k\right)={\dot{x}_{i}}\left(k\right)-{\dot{x}_{i}}\left(k-1\right)=({{x}_{i}}\left(k\right)-{{x}_{i}}\left(k-1\right))-({{x}_{i}}\left(k-1\right)-{{x}_{i}}\left(k-2\right))$,
$i=1,2,3,4$, and the $x_{1}$, $x_{2}$, $x_{3}$ and $x_{4}$ denote the dew
point mean, the heat index mean, the wind chill mean and the average wind
speed, respectively. The training loss function is

$\displaystyle\mathcal{L}={\left\|{{[\widehat{\mathbf{y}}]_{1}}-{x_{2}}\left({k+1}\right)}\right\|}+{\left\|{{[\widehat{\mathbf{y}}]_{2}}-{x_{4}}\left({k+1}\right)}\right\|},$

which indicates the network models are to predict the heat index mean and the
average wind speed for next hour.

####  6.3.1 Phy-Augmentation

Figure 16: Architectures of deep PhN and classical DNN.

We use the data of station $S_{1}$ to study the influence of augmentation
order $r$ from the perspective of training loss. The considered network model
is shown in Figure 16 (a). The trajectories of training loss under two
different augmentation orders are presented in Figure 17 (a). It is
straightforward to observe from the figure that the input augmentation
significantly speeds up the training process and reduces the training loss.
The intuitive explanation is the input augmentation enlarges node number.

To validate the statement, we perform the comparisons with the classical DNN,
whose structure is given in Figure 16 (b). To guarantee the comparisons are
carried in the fair settings, we let the input dimension (i.e., $n$) of the
final layer of DNN in Figure 16 (b) equates to the output dimension of input
augmentation of PhN in Figure 16 (a). According to (112), we let
$n={\rm{len}}(\mathfrak{m}(\mathbf{x},r=5))=253$. The comparison results are
presented in Figure 17 (b), which indicates that given the same number of
nodes, the deep PhN still has much faster convergence speed and smaller
training loss of mean. The phenomenon can be explained by the fact that Phy-
Augmentation well captures physical features.

Figure 17: Deep PhN v.s. classical DNN: log of training loss (5 random seeds).

####  6.3.2 Large Noise

To test the robustness of deep PhN faced with large noise, we consider the
network structure in Figure 16 (c). To introduce noisy datasets for training,
we use the inputs $\mathbf{x}(k)$ from stations $S_{2}$–$S_{4}$, while the
dataset of outputs $[\widehat{\mathbf{y}}]_{1}$ and
$[\widehat{\mathbf{y}}]_{2}$ is from the station $S_{1}$. This setting means
we will use station $S_{1}$’s neighbors’ inputs to predict its heat index mean
and average wind speed. For the suppressor of deep PhN in (11), we let
$\beta=90$ and $\alpha=-1$. The trajectories of training loss in Figure 18
together with station locations in Figure 15 show that the training loss
decreases as the distance with station $S_{1}$ increases. Considering the fact
that noise of training (input) datasets can increase as distance with (truth)
station $S_{1}$ increases, the result demonstrates the superiority of
suppressor in mitigating the influence of large noise in augmented input
features.

Figure 18: Training loss with and without suppressor function (5 random
seeds).

##  7 Discussion

In this article, we have proposed a physics-model-based deep neural network
framework, called Phy-Taylor, that addresses the challenge the purely data-
driven deep DNNs are facing: the violation of inferred relations with physics
laws. The Phy-Taylor framework introduces two contributions: the deep physics-
compatible neural networks (PhNs) and a physics-guided neural network (NN)
editing mechanism, aiming at ensuring strict compliance with prior physical
knowledge. The PhN aims to directly capture nonlinearities inspired by
physical quantities, such as kinetic energy and aerodynamic drag force. The NN
editing mechanism further modifies network links and activation functions
consistently with physical knowledge to avoid spurious correlations. As an
extension, we have also proposed a self-correcting Phy-Taylor framework that
introduces a core capability of automatic output correction when violation of
safety occurs. The self-correcting Phy-Taylor successfully addresses the
dilemma of prediction horizon and computation time that nonlinear model-
predictive control and control barrier function are facing in safety- and
time-critical systems.

The experiments show that through physics-guided NN editing and direct capture
of hard-to-learn nonlinearities, the Phy-Taylor exhibits a considerable
reduction in learning parameters, a remarkably accelerated training process,
and greatly enhanced model robustness and accuracy. This suggests that
building on the Phy-Taylor framework, the concurrently energy-efficient,
robustness and high-accurate DNNs are promising for the energy-constraint
physical engineering systems (see e.g., internet of things and battery-powered
drones), which constitutes a part of our future research. The current Phy-
Taylor however suffers from curse of dimensionality, so it can hardly be
applied to high-dimensional data, such as image and text. Overcoming the
dimensionality curse problem will thus make the Phy-Taylor scalable. The
Tucker Decomposition is a potential to address this problem, since it can
decompose the higher-order derivatives in Taylor expansions parameterized by
DNNs into small core tensors and a set of factor matrices.

##  8 Methods

###  8.1 Datasets

All the datasets used in the experiments are publicly available at
https://waynestateprod-
my.sharepoint.com/:f:/g/personal/hm9062_wayne_edu/EoXO99WN8zJEidtGRj-
dISIBQPBWCnL_Ji6QOZ1uJACjug.

The datasets in Experiment 6.1 are collected from the AutoRally platform [45].
The vehicle therein is set to run 6 times, each running time is 30 minutes,
which can generate around 1100 samples of ellipse trajectories. The datasets
in Experiment 6.2 are generated by solving the differential equations of
coupled pendulum in MATLAB using the ode45 solver, given 1000 initial
conditions. The climate datasets in Experiment 6.3 are from the Climate Data
Online–National Oceanic and Atmospheric Administration222 NOAA:
https://www.ncdc.noaa.gov/cdo-
web/datasets/NORMAL_HLY/locations/FIPS:17/detail, whose period of record is
01/01/2010–12/31/2010. The datasets in Experiments 6.1 and 6.2 are evenly
distributed to 6 files: two files are preserved for testing and validating,
the remaining 4 files are used for training. For the Experiment 6.3, the data
over the final 100 hours is preserved for testing, the data over another 864
hours are used for validating, the remaining data are used for training.

###  8.2 Code

For the code, we use the Python API for the TensorFlow framework [46] and the
Adam optimizer [47] for training. The Python version is 2.7.12. The TensorFlow
version is 1.4.0. Our source code is publicly available at GitHub:
https://github.com/ymao578/Phy-Taylor. The source code of implementing self-
correcting Phy-Taylor in the AutoRally platform is publicly available at
https://waynestateprod-
my.sharepoint.com/:f:/g/personal/hm9062_wayne_edu/EnDmRmmbKlJFmlQfC3qLMHwBf4KG9FRGVVo3FFVk1TrZeg?e=TQZSgr.

###  8.3 Training

We set the batch-size to 200 for the Experiments 6.1 and 6.2, while the batch-
size of Experiment 6.3 is 150. The learning rates of Experiments 6.1–6.3 are
set to 0.0005, 0.00005 and 0.00005, respectively. In all the experiments, each
weight matrix is initialized randomly from a (truncated) normal distribution
with zero mean and standard deviation, discarding and re-drawing any samples
that are more than two standard deviations from the mean. We initialize each
bias according to the normal distribution with zero mean and standard
deviation.

## References

  * [1] Zachary, A. & Helen, T.  AI Accidents: An emerging threat.  _Center for Security and Emerging Technology_ DOI: https://doi.org/10.51593/20200072 (2021).
  * [2] Wang, R. & Yu, R.  Physics-guided deep learning for dynamical systems: A survey.  _arXiv preprint_ https://arxiv.org/abs/2107.01272.
  * [3] Willard, J., Jia, X., Xu, S., Steinbach, M. & Kumar, V.  Integrating scientific knowledge with machine learning for engineering and environmental systems.  _ACM Computing Surveys_ (2021).
  * [4] Jia, X. _et al._ Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles.  _ACM/IMS Transactions on Data Science_ 2, 1–26 (2021).
  * [5] Jia, X. _et al._ Physics-guided RNNs for modeling dynamical systems: A case study in simulating lake temperature profiles.  In _Proceedings of the 2019 SIAM International Conference on Data Mining_ , 558–566 (2019).
  * [6] Wang, S. & Perdikaris, P.  Deep learning of free boundary and stefan problems.  _Journal of Computational Physics_ 428, 109914 (2021).
  * [7] Lu, L. _et al._ Physics-informed neural networks with hard constraints for inverse design.  _SIAM Journal on Scientific Computing_ 43, B1105–B1132 (2021).
  * [8] Chen, Y. _et al._ Theory-guided hard constraint projection (HCP): A knowledge-based data-driven scientific machine learning method.  _Journal of Computational Physics_ 445, 110624 (2021).
  * [9] Wang, N., Zhang, D., Chang, H. & Li, H.  Deep learning of subsurface flow via theory-guided neural network.  _Journal of Hydrology_ 584, 124700 (2020).
  * [10] Xu, K. & Darve, E.  Physics constrained learning for data-driven inverse modeling from sparse observations.  _Journal of Computational Physics_ 110938 (2022).
  * [11] Karniadakis, G. E. _et al._ Physics-informed machine learning.  _Nature Reviews Physics_ 3, 422–440 (2021).
  * [12] Wang, R., Kashinath, K., Mustafa, M., Albert, A. & Yu, R.  Towards physics-informed deep learning for turbulent flow prediction.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 1457–1466 (2020).
  * [13] Daw, A., Karpatne, A., Watkins, W., Read, J. & Kumar, V.  Physics-guided neural networks (PGNN): An application in lake temperature modeling.  _arXiv preprint_ https://arxiv.org/abs/1710.11431.
  * [14] Cranmer, M. _et al._ Lagrangian neural networks.  In _ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations_ (2020).
  * [15] Finzi, M., Wang, K. A. & Wilson, A. G.  Simplifying Hamiltonian and Lagrangian neural networks via explicit constraints.  _Advances in Neural Information Processing Systems_ 33, 13880–13889 (2020).
  * [16] Greydanus, S., Dzamba, M. & Yosinski, J.  Hamiltonian neural networks.  _Advances in Neural Information Processing Systems_ 32 (2019).
  * [17] Muralidhar, N. _et al._ Phynet: Physics guided neural networks for particle drag force prediction in assembly.  In _Proceedings of the 2020 SIAM International Conference on Data Mining_ , 559–567 (2020).
  * [18] Masci, J., Boscaini, D., Bronstein, M. & Vandergheynst, P.  Geodesic convolutional neural networks on riemannian manifolds.  In _Proceedings of the IEEE international conference on computer vision workshops_ , 37–45 (2015).
  * [19] Monti, F. _et al._ Geometric deep learning on graphs and manifolds using mixture model cnns.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 5115–5124 (2017).
  * [20] Horie, M., Morita, N., Hishinuma, T., Ihara, Y. & Mitsume, N.  Isometric transformation invariant and equivariant graph convolutional networks.  _arXiv preprint_ https://arxiv.org/abs/2005.06316.
  * [21] Wang, R.  Incorporating symmetry into deep dynamics models for improved generalization.  In _International Conference on Learning Representations_ (2021).
  * [22] Li, Y., He, H., Wu, J., Katabi, D. & Torralba, A.  Learning compositional Koopman operators for model-based control.  _arXiv preprint_ https://arxiv.org/abs/1910.08264.
  * [23] Lusch, B., Kutz, J. N. & Brunton, S. L.  Deep learning for universal linear embeddings of nonlinear dynamics.  _Nature communications_ 9, 1–10 (2018).
  * [24] Li, Z. _et al._ Fourier neural operator for parametric partial differential equations.  _arXiv:2010.08895_ (2020).
  * [25] Königsberger, K.  _Analysis 2_ (Springer-Verlag, 2013).
  * [26] Yang, Y.-Y. & Chaudhuri, K.  Understanding rare spurious correlations in neural networks.  _arXiv preprint_ https://arxiv.org/abs/2202.05189.
  * [27] Sagawa, S., Raghunathan, A., Koh, P. W. & Liang, P.  An investigation of why overparameterization exacerbates spurious correlations.  In _International Conference on Machine Learning_ , 8346–8356 (2020).
  * [28] Rajamani, R.  _Vehicle dynamics and control_ (Springer Science & Business Media, 2011).
  * [29] Kani, J. N. & Elsheikh, A. H.  DR-RNN: A deep residual recurrent neural network for model reduction.  _arXiv preprint_ https://arxiv.org/abs/1709.00939.
  * [30] Belbute-Peres, F. D. A., Economon, T. & Kolter, Z.  Combining differentiable PDE solvers and graph neural networks for fluid flow prediction.  In _International Conference on Machine Learning_ , 2402–2411 (2020).
  * [31] Wu, D. _et al._ DeepGLEAM: A hybrid mechanistic and deep learning model for COVID-19 forecasting.  _arXiv preprint_ https://arxiv.org/abs/2102.06684 (2021).
  * [32] Guen, V. L. & Thome, N.  Disentangling physical dynamics from unknown factors for unsupervised video prediction.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 11474–11484 (2020).
  * [33] Garcia Satorras, V., Akata, Z. & Welling, M.  Combining generative and discriminative models for hybrid inference.  _Advances in Neural Information Processing Systems_ 32 (2019).
  * [34] Long, Y., She, X. & Mukhopadhyay, S.  HybridNet: Integrating model-based and data-driven learning to predict evolution of dynamical systems.  In _Conference on Robot Learning_ , 551–560 (2018).
  * [35] Yin, Y. _et al._ Augmenting physical models with deep networks for complex dynamics forecasting.  _Journal of Statistical Mechanics: Theory and Experiment_ 2021, 124012 (2021).
  * [36] AI incident database.  _AID_ DOI: https://incidentdatabase.ai/summaries/incidents.
  * [37] TOM, K. & STEFANIE, D.  Felony charges are 1st in a fatal crash involving autopilot.  _AP News_ DOI: https://apnews.com/article/tesla-autopilot-fatal-crash-charges-91b4a0341e07244f3f03051b5c2462ae (2022).
  * [38] Hawkins, A. J.  Tesla didn’t fix an autopilot problem for three years, and now another person is dead.  _The Verge_ DOI: https://www.theverge.com/2019/5/17/18629214/tesla-autopilot-crash-death-josh-brown-jeremy-banner (2019).
  * [39] Ames, A. D. _et al._ Control barrier functions: Theory and applications.  In _2019 18th European Control Conference_ , 3420–3431 (2019).
  * [40] Singletary, A., Swann, A., Chen, Y. & Ames, A.  Onboard safety guarantees for racing drones: High-speed geofencing with control barrier functions.  _IEEE Robotics and Automation Letters_ (2022).
  * [41] Williams, G., Drews, P., Goldfain, B., Rehg, J. M. & Theodorou, E. A.  Information-theoretic model predictive control: Theory and applications to autonomous driving.  _IEEE Transactions on Robotics_ 34, 1603–1622 (2018).
  * [42] Falcone, P., Borrelli, F., Asgari, J., Tseng, H. E. & Hrovat, D.  Predictive active steering control for autonomous vehicle systems.  _IEEE Transactions on Control Systems Technology_ 15, 566–580 (2007).
  * [43] Zeng, J., Zhang, B. & Sreenath, K.  Safety-critical model predictive control with discrete-time control barrier function.  In _2021 American Control Conference_ , 3882–3889 (2021).
  * [44] Doyle, J. C., Francis, B. A. & Tannenbaum, A. R.  _Feedback control theory_ (Courier Corporation, 2013).
  * [45] Goldfain, B. _et al._ AutoRally: An open platform for aggressive autonomous driving.  _IEEE Control Systems Magazine_ 39, 26–55 (2019).
  * [46] Kingma, D. P. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint_ https://arxiv.org/abs/1412.6980.
  * [47] Abadi, M. _et al._ Tensorflow: Large-scale machine learning on heterogeneous distributed systems.  _arXiv preprint_ https://arxiv.org/abs/1603.04467.
  * [48] Stanley, R. P.  What is enumerative combinatorics?  In _Enumerative combinatorics_ , 1–63 (Springer, 1986).

##  9 Supplementary Information

###  9.1 Auxiliary Theorems

######  Theorem 4.

The DNR magnitude of high-order monomial
$[\bar{\mathbf{x}}]_{i}^{p}[\bar{\mathbf{x}}]_{j}^{q}$, $p$, $q\in\mathbb{N}$,
is strictly increasing with respect to $|\mathrm{DNR}_{i}|$ and
$|\mathrm{DNR}_{j}|$, if

$\displaystyle\mathrm{DNR}_{i},\leavevmode\nobreak\
\mathrm{DNR}_{j}\in(-\infty,-1]\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{or}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{DNR}_{i},\leavevmode\nobreak\
\mathrm{DNR}_{j}\in[-\frac{1}{2},0)\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{or}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{DNR}_{i},\leavevmode\nobreak\ \mathrm{DNR}_{j}\in(0,\infty).$ (103)

###### Proof.

In view of definition 2, the true data can be equivalently expressed as
$[{\mathbf{h}}]_{i}=\mathrm{DNR}_{i}\cdot[{\mathbf{w}}]_{i}$, according to
which we have
$[\bar{\mathbf{x}}]_{i}=(1+\mathrm{DNR}_{i}){[{\mathbf{w}}]_{i}}$ such that

$\displaystyle[\bar{\mathbf{x}}]_{i}^{p}[\bar{\mathbf{x}}]_{j}^{q}={({1+\mathrm{DNR}_{i}})^{p}}{({1+\mathrm{DNR}_{j}})^{q}}[\mathbf{w}]_{i}^{p}[\mathbf{w}]_{j}^{q},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
[\mathbf{h}]_{i}^{p}[\mathbf{h}]_{j}^{q}=\mathrm{DNR}_{i}^{p}\cdot\mathrm{DNR}_{j}^{q}\cdot[\mathbf{w}]_{i}^{p}[\mathbf{w}]_{j}^{q}.$
(104)

We note the true data of high-order monomial
$[\bar{\mathbf{x}}]_{i}^{p}[\bar{\mathbf{x}}]_{j}^{q}$ is
$[{\mathbf{h}}]_{i}^{p}[{\mathbf{h}}]_{j}^{q}$, the corresponding noise can
thus be derived from the formula (104) as

$\displaystyle[\bar{\mathbf{x}}]_{i}^{p}[\bar{\mathbf{x}}]_{j}^{q}-[\mathbf{h}]_{i}^{p}[\mathbf{h}]_{j}^{q}=\left[{({1+\mathrm{DNR}_{i}})^{p}}{({1+\mathrm{DNR}_{j}})^{q}}-\mathrm{DNR}_{i}^{p}\cdot\mathrm{DNR}_{j}^{q}\right][\mathbf{w}]_{i}^{p}[\mathbf{w}]_{j}^{q},$

which, in conjunction with the second formula in Equation (104), lead to

$\displaystyle\left|\mathrm{DNR}^{p+q}_{ij}\right|\triangleq\left|\frac{{[\mathbf{h}]_{i}^{p}[\mathbf{h}]_{j}^{q}}}{{[\bar{\mathbf{x}}]_{i}^{p}[\bar{\mathbf{x}}]_{j}^{q}-[\mathbf{h}]_{i}^{p}[\mathbf{h}]_{j}^{q}}}\right|=\left|\frac{1}{{{{\left({1+\frac{1}{{\mathrm{DNR}_{i}}}}\right)}^{p}}{{\left({1+\frac{1}{{\mathrm{DNR}_{j}}}}\right)}^{q}}-1}}\right|,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ p,\leavevmode\nobreak\ q\in\mathbb{N}.$ (105)

We can straightforwardly verify from formula (105) that if
$\mathrm{DNR}_{i},\leavevmode\nobreak\ \mathrm{DNR}_{j}\in(0,\infty)$, we have

$\displaystyle\left|\mathrm{DNR}^{p+q}_{ij}\right|=\frac{1}{{{{\left({1+\frac{1}{{|\mathrm{DNR}_{i}|}}}\right)}^{p}}{{\left({1+\frac{1}{{|\mathrm{DNR}_{j}|}}}\right)}^{q}}-1}},$
(106)

which implies $|\mathrm{DNR}^{p+q}_{ij}|$ is strictly increasing with respect
to $|\mathrm{DNR}_{i}|$ and $|\mathrm{DNR}_{j}|$ under this condition.

The condition $\mathrm{DNR}_{i},\leavevmode\nobreak\
\mathrm{DNR}_{j}\in(-\infty,-1]$ means that

$\displaystyle\frac{{1}}{{\mathrm{DNR}_{j}}}\in[-1,0),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\frac{{1}}{{\mathrm{DNR}_{j}}}\in[-1,0),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
1+\frac{{1}}{{\mathrm{DNR}_{i}}}\in[0,1),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 1+\frac{{1}}{{\mathrm{DNR}_{j}}}\in[0,1),$ (107)

considering which, the formula (105) equivalently transforms to

$\displaystyle\left|\mathrm{DNR}^{p+q}_{ij}\right|=\frac{1}{{{{1-\left({1+\frac{1}{{\mathrm{DNR}_{i}}}}\right)}^{p}}{{\left({1+\frac{1}{{\mathrm{DNR}_{j}}}}\right)}^{q}}}}=\frac{1}{{{{1-\left({1-\frac{1}{{|\mathrm{DNR}_{i}|}}}\right)}^{p}}{{\left({1-\frac{1}{{|\mathrm{DNR}_{j}|}}}\right)}^{q}}}},$
(108)

which reveals that $|\mathrm{DNR}^{p+q}_{ij}|$ is strictly increasing with
respect to $|\mathrm{DNR}_{i}|$ and $|\mathrm{DNR}_{j}|$.

The condition $\mathrm{DNR}_{i},\leavevmode\nobreak\
\mathrm{DNR}_{j}\in[-\frac{1}{2},0)$ means

$\displaystyle\frac{{1}}{{\mathrm{DNR}_{j}}}\in(-\infty,-2],\leavevmode\nobreak\
\leavevmode\nobreak\
\frac{{1}}{{\mathrm{DNR}_{j}}}\in(-\infty,-2],\leavevmode\nobreak\
\leavevmode\nobreak\
1+\frac{{1}}{{\mathrm{DNR}_{i}}}\in(-\infty,-1],\leavevmode\nobreak\
\leavevmode\nobreak\ 1+\frac{{1}}{{\mathrm{DNR}_{j}}}\in(-\infty,-1],$ (109)

in light of which, the formula (105) can equivalently express as

  * •

if $p+q$ is even,

$\displaystyle\left|\mathrm{DNR}^{m+n}_{ij}\right|=\frac{1}{{{{\left|{\frac{1}{{|\mathrm{DNR}_{i}|}}-1}\right|}^{p}}{{\left|{\frac{1}{{|\mathrm{DNR}_{j}|}}-1}\right|}^{q}}}-1},$
(110)

  * •

if $p+q$ is odd,

$\displaystyle\left|\mathrm{DNR}^{p+q}_{ij}\right|=\frac{1}{{{{1+\left|{\frac{1}{{|\mathrm{DNR}_{i}|}}-1}\right|}^{p}}{{\left|{\frac{1}{{|\mathrm{DNR}_{j}|}}-1}\right|}^{q}}}}.$
(111)

We note both the functions (110) and (111) imply $|\mathrm{DNR}^{m+n}_{ij}|$
is strictly increasing with respect to $|\mathrm{DNR}_{i}|$ and
$|\mathrm{DNR}_{j}|$, which completes the proof. ∎

######  Theorem 5.

[48] For any pair of positive integers $n$ and $k$, the number of $n$-tuples
of non-negative integers whose sum is $r$ is equal to the number of multisets
of cardinality $n-1$ taken from a set of size $n+r-1$, i.e.,
$\left(\begin{array}[]{l}n+r-1\\\
n-1\end{array}\right)=\frac{{\left({n+r-1}\right)!}}{{\left({n-1}\right)!r!}}$.

######  Theorem 6.

The space complexity of Phy-Augmentation, i.e., the dimension of terminal
output generated by Algorithm 1, in PhN layer (20) is

$\displaystyle\mathrm{len}(\mathfrak{m}({\mathbf{x}},r))=\sum\limits_{s=1}^{r}{\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)!s!}}}+1.$
(112)

###### Proof.

We denote the output from Line 1 of Algorithm 1 by $\overline{\mathbf{x}}$.
Let us first consider the case $r=1$. In this case, the Algorithm 1 skips the
Lines 1–1 and arrives at $\mathfrak{m}(\mathbf{x},r)=[1;\leavevmode\nobreak\
\overline{\mathbf{x}}]$ in Line 1. Noticing from Line 1 that
$\overline{\mathbf{x}}\in\mathbb{R}^{n}$, we obtain
$\text{len}(\mathfrak{m}(\mathbf{x},r))=n+1$, which verifies the correctness
of (112) with $r=1$.

We next consider the case $r\geq 2$. Given the input dimension $n$ and an
order $s\in\\{2,\ldots,r-1,r\\}$, the Lines 1–1 of Algorithm 1 are to generate
all the non-missing and non-redundant monomials included in
${\left({\sum\limits_{i=1}^{n}{{[\overline{\mathbf{x}}]_{i}}}}\right)^{s}}$.
The problem of the number of generated monomials via Algorithm 1 is equivalent
to the problem that for any pair of positive integers $n$ and $s$, the number
of $n$-tuples of non-negative integers (whose sum is $s$) is equal to the
number of multisets of cardinality $n-1$ taken from a set of size $n+s-1$.
Additionally, we note that the vector generated in Line 1 of Algorithm 1,
denoted by $\widetilde{\mathfrak{m}}(\mathbf{x},s)$, stacks all the generated
monomials. According to the auxiliary Theorem 5 in Supplementary Information
9.1, we then have
$\text{len}(\widetilde{\mathfrak{m}}(\mathbf{x},s))=\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)!s!}}$.
Finally, we note that the Lines 1, and 1 of Algorithm 1 imply that the
generated vector $\mathfrak{m}(\mathbf{x},r)$ stack the $1$ with
$\widetilde{\mathfrak{m}}(\mathbf{x},s)$ over $s\in\\{1,\ldots,r-1,r\\}$,
respectively. We thus can obtain (112). ∎

###  9.2 Proof of Theorem 1

We note the $[\widetilde{\mathbf{h}}]_{i}$ given in Equation (13) can be
written as

$\displaystyle[\widetilde{\mathbf{h}}]_{i}=\begin{cases}[\mathbf{h}]_{i},&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}<0\\\
[\mathbf{h}]_{i},&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}<0\\\
[\mathbf{h}]_{i}\cdot\kappa_{i}+\rho_{i},&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}>0\end{cases},$ (113)

subtracting $[{\mathbf{h}}]_{i}$ from which yields

$\displaystyle\left|[\widetilde{\mathbf{h}}]_{i}-[{\mathbf{h}}]_{i}\right|$
$\displaystyle=\begin{cases}0,&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}<0\\\ 0,&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}<0\\\
\left|{[\mathbf{h}]_{i}\cdot\left({\kappa_{i}-1}\right)+\rho_{i}}\right|,&\text{if}\leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ [\mathbf{w}]_{i}>0\\\ \end{cases}.$ (114)

Referring to the output $\chi([\bar{\bf{x}}]_{i})$ of suppressor in Equation
(10), we can conclude the $[\widetilde{\mathbf{h}}]_{i}$ given in Equation
(113) is the true data of suppressor output. Subtracting the
$[\widetilde{\mathbf{h}}]_{i}$ from the $\chi([\bar{\bf{x}}]_{i})$ results in
the noise $[\widetilde{\mathbf{w}}]_{i}$ of suppressor output, given in
Equation (13).

To prove the property (12), we consider the following three cases:

  * •

If $[\mathbf{h}]_{i}+[\mathbf{w}]_{i}<0$, we obtain from the the first item of
$[\widetilde{\mathbf{h}}]_{i}$ in Equation (113) and
$[\widetilde{\mathbf{w}}]_{i}$ in Equation (13) that
$\mathrm{DNR}_{i}=\frac{[\widetilde{\mathbf{h}}]_{i}}{[\widetilde{\mathbf{w}}]_{i}}=-1$.

  * •

If $[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0$ and $[\mathbf{w}]_{i}<0$, we have
$[\mathbf{h}]_{i}>0$ and $[\mathbf{h}]_{i}>-[\mathbf{w}]_{i}|>0$. We then
obtain from the second item of $[\widetilde{\mathbf{h}}]_{i}$ in Equation
(113) and $[\widetilde{\mathbf{w}}]_{i}$ in Equation (13) that
$\mathrm{DNR}_{i}=\frac{[\widetilde{\mathbf{h}}]_{i}}{[\widetilde{\mathbf{w}}]_{i}}=\frac{[{\mathbf{h}}]_{i}}{[{\mathbf{w}}]_{i}}<-1$.

  * •

If $[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0$ and $[\mathbf{w}]_{i}>0$, we
obtain from the third item of $[\widetilde{\mathbf{h}}]_{i}$ in Equation (113)
and $[\widetilde{\mathbf{w}}]_{i}$ in Equation (13) that
$\mathrm{DNR}_{i}=\frac{[\widetilde{\mathbf{h}}]_{i}}{[\widetilde{\mathbf{w}}]_{i}}=\frac{{[\mathbf{h}]_{i}\cdot\kappa_{i}+\rho_{i}}}{{[\mathbf{w}]_{i}\cdot\kappa_{i}}}$.
Recalling $[\widetilde{\mathbf{w}}]_{i}>0$, if $\kappa_{i}>0$, the
$\frac{{[\mathbf{h}]_{i}\cdot\kappa_{i}+\rho_{i}}}{{[\mathbf{w}]_{i}\cdot\kappa_{i}}}\leq-1$
is equivalent to

$\displaystyle\rho_{i}\leq-([\mathbf{h}]_{i}+[\mathbf{w}]_{i})\kappa_{i}<0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{with}\leavevmode\nobreak\
\leavevmode\nobreak\ \kappa_{i}>0,\leavevmode\nobreak\ \leavevmode\nobreak\
[\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0.$ (115)

If $\kappa_{i}<0$, the
$\frac{{[\mathbf{h}]_{i}\cdot\kappa_{i}+\rho_{i}}}{{[\mathbf{w}]_{i}\cdot\kappa_{i}}}\leq-1$
is equivalent to

$\displaystyle\rho_{i}\geq-\left({[\mathbf{w}]_{i}+[\mathbf{h}]_{i}}\right)\kappa_{i}\geq
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \kappa_{i}<0,\leavevmode\nobreak\
\leavevmode\nobreak\ [\mathbf{h}]_{i}+[\mathbf{w}]_{i}\geq 0.$ (116)

We finally conclude from Equations (115) and (116) that
$\mathrm{DNR}_{i}=\frac{[\widetilde{\mathbf{h}}]_{i}}{[\widetilde{\mathbf{w}}]_{i}}\in(-\infty,-1]$
under the condition (11). Then, according to Theorem 4, we arrive in the
property (12), which completes the proof.

###  9.3 Proof of Theorem 2

Let us first consider the first PhN layer, i.e., the case $t=1$. The Line 2 of
Algorithm 2 means that the knowledge matrix $\mathbf{K}_{\left\langle
1\right\rangle}$ includes all the known model-substructure parameters, whose
corresponding entries in the masking matrix $\mathbf{M}_{\left\langle
1\right\rangle}$ (generated in the Line 2 of Algorithm 2) are frozen to be
zeros. Consequently, both $\mathbf{M}_{\left\langle
1\right\rangle}\odot\mathbf{A}$ and $\mathbf{U}_{\left\langle
1\right\rangle}=\mathbf{M}_{\left\langle
1\right\rangle}\odot\mathbf{W}_{\left\langle 1\right\rangle}$ excludes all the
known model-substructure parameters (included in $\mathbf{K}_{\left\langle
1\right\rangle}$). With the consideration of Definition 1, we thus conclude
that $\mathbf{M}_{\left\langle
1\right\rangle}\odot\mathbf{A}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+\mathbf{f}(\mathbf{x})$ in the ground-truth model (16) and
$\mathbf{a}_{\left\langle
1\right\rangle}\odot\text{act}\left(\mathbf{U}_{\left\langle
1\right\rangle}\cdot{\mathfrak{m}}\left(\mathbf{x},r_{\left\langle
1\right\rangle}\right)\right)$ in the output computation (17) are independent
of the term $\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle 1\right\rangle})$.
Moreover, the activation-masking vector (generated in Line 2 of Algorithm 2)
indicates that the activation function corresponding to the output’s $i$-th
entry is inactive, if the all the entries in the $i$-th row of masking matrix
are zeros (implying all the entries in the $i$-th row of weight matrix are
known model-substructure parameters). Finally, we arrive in the conclusion
that the input/output (i.e., $\mathbf{x}$/$\mathbf{y}_{\left\langle
1\right\rangle}$) of the first PhN layer strictly complies with the available
physical knowledge pertaining to the ground truth (16), i.e., if the
$[\mathbf{A}]_{i,j}$ is a known model-substructure parameter, the
$\frac{{\partial{[\mathbf{y}_{\left\langle
1\right\rangle}]_{i}}}}{{\partial[{\mathfrak{m}}\left({\mathbf{x},r}\right)]_{j}}}\equiv\frac{{\partial{[\mathbf{y}]_{i}}}}{{\partial{[\mathfrak{m}}({\mathbf{x},r})]_{j}}}\equiv{[\mathbf{A}]_{i,j}}$
always holds.

We next consider the remaining PhN layers. Considering the Line 2 Algorithm 2,
we have

$\displaystyle[\mathbf{y}_{\left\langle
p\right\rangle}]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=[\mathbf{K}_{\left\langle
p\right\rangle}\cdot\mathfrak{m}(\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle})]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=\mbox{\small
I}_{\text{len}(\mathbf{y})}\cdot[\mathfrak{m}(\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle})]_{2:(\text{len}(\mathbf{y})+1)}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$ (117a)
$\displaystyle=\mbox{\small
I}_{\text{len}(\mathbf{y})}\cdot[\mathbf{y}_{\left\langle
p-1\right\rangle}]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$ (117b)
$\displaystyle=[\mathbf{y}_{\left\langle
p-1\right\rangle}]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=[\mathbf{K}_{\left\langle
p-1\right\rangle}\cdot\mathfrak{m}(\mathbf{y}_{\left\langle
p-2\right\rangle},r_{\left\langle
p-1\right\rangle})]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=\mbox{\small
I}_{\text{len}(\mathbf{y})}\cdot[\mathfrak{m}(\mathbf{y}_{\left\langle
p-2\right\rangle},r_{\left\langle
p-1\right\rangle})]_{2:(\text{len}(\mathbf{y})+1)}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=\mbox{\small
I}_{\text{len}(\mathbf{y})}\cdot[\mathbf{y}_{\left\langle
p-2\right\rangle}]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=[\mathbf{y}_{\left\langle
p-2\right\rangle}]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot\breve{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$ $\displaystyle=\ldots$
$\displaystyle=[\mathbf{y}_{\left\langle
1\right\rangle}]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=[\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})]_{1:\text{len}(\mathbf{y})}+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})},$
$\displaystyle=\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})},$ (117c)

where (117a) and (117b) are obtained from their previous steps via considering
the structure of block matrix $\mathbf{K}_{\left\langle t\right\rangle}$
(generated in Line 11 of Algorithm 2) and the formula of augmented monomials:
$\mathfrak{m}(\mathbf{x},r)=\left[1;\leavevmode\nobreak\
\mathbf{x};\leavevmode\nobreak\
[\mathfrak{m}({\mathbf{x},r})]_{(\text{len}(\mathbf{x})+2):\text{len}(\mathfrak{m}(\mathbf{x},r))}\right]$
(generated via Algorithm 1). The remaining iterative steps follow the same
path.

The training loss function is to push the terminal output of Algorithm 2
(i.e., $\widehat{\mathbf{y}}=\mathbf{y}_{\left\langle p\right\rangle}$) to
approximate the real output $\mathbf{y}$, which in light of (117c) yields

$\displaystyle\widehat{\mathbf{y}}$ $\displaystyle=\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+[\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle
p\right\rangle}})}\right)]_{1:\text{len}(\mathbf{y})}$
$\displaystyle=\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle
1\right\rangle})+\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle p\right\rangle}})}\right),$ (118)

where (118) from its previous step is obtained via considering the fact
$\text{len}(\widehat{\mathbf{y}})=\text{len}(\mathbf{y})=\text{len}(\mathbf{y}_{\left\langle
p\right\rangle})$. Meanwhile, the condition of generating weight-masking
matrix in Line 2 of Algorithm 2 removes all the node-representations’
connections with the known model-substructure parameters included in
$\mathbf{K}_{\left\langle 1\right\rangle}$. Therefore, we can conclude that in
the terminal output computation (118), the term $\mathbf{a}_{\left\langle
p\right\rangle}\odot\text{act}\\!\left({\mathbf{U}_{\left\langle
p\right\rangle}\cdot{\mathfrak{m}}({\mathbf{y}_{\left\langle
p-1\right\rangle},r_{\left\langle p\right\rangle}})}\right)$ does not have
influence on the computing of knowledge term $\mathbf{K}_{\left\langle
1\right\rangle}\cdot\mathfrak{m}(\mathbf{x},r_{\left\langle 1\right\rangle})$.
Thus, the Algorithm 2 strictly embeds and preserves the available knowledge
pertaining to the physics model of ground truth (1), or equivalently the (16).

###  9.4 Proof of Theorem 3

Due to Theorem 6 in Supplementary Information 9.1, the number of augmented
monomials of $d$ cascading PhNs (21) is obtained as

$\displaystyle\sum\limits_{p=1}^{d}{\text{len}(\mathfrak{m}({\mathbf{x}},r_{{\left\langle
p\right\rangle}}))}$
$\displaystyle=\underbrace{\sum\limits_{s=1}^{{r_{\left\langle
1\right\rangle}}}{\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)s!}}}+1}_{\text{the
first
PhN}}+\underbrace{\sum\limits_{v=1}^{d-1}{\sum\limits_{s=1}^{{r_{{\left\langle
v+1\right\rangle}}}}{\frac{{\left({{n_{\left\langle
v\right\rangle}}+s-1}\right)!}}{{\left({{n_{\left\langle
v\right\rangle}}-1}\right)!s!}}}}+d-1}_{\text{the remaining PhNs}}.$ (119)

The condition (22) implies that $r>r_{{\left\langle 1\right\rangle}}$, which
in conjunction with (112), lead to

$\displaystyle\text{len}(\mathfrak{m}({\mathbf{x}},r))=\sum\limits_{s=1}^{{r_{\left\langle
1\right\rangle}}}{\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)!s!}}}+\sum\limits_{s={r_{\left\langle
1\right\rangle}}+1}^{r}{\frac{{\left({n+s-1}\right)!}}{{\left({n-1}\right)!s!}}}+1.$
(120)

Subtracting (119) from (120) yields (23).

###  9.5 Derivations of Solution (99)

With the consideration of (81)–(94), the safety formulas:
$\left[{\bf{s}}({\bf{u}}(k))\right]_{1}=[\widehat{\mathbf{c}}]_{1}$ and
$\left[{\bf{s}}({\bf{u}}(k))\right]_{2}=[\widehat{\mathbf{c}}]_{2}$ can be
rewritten as

$\displaystyle{\lambda_{1}}[\widehat{\mathbf{u}}(k)]_{1}^{2}+{\lambda_{2}}[\widehat{\mathbf{u}}(k)]_{2}^{2}$
$\displaystyle={[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}},$ (121)
$\displaystyle{\left[\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\right]^{\top}}{\mathbf{P}}_{2}\left[\begin{array}[]{l}\theta(k)\\\
\gamma(k)\end{array}\right]$
$\displaystyle={s_{11}}[\widehat{\mathbf{u}}(k)]_{1}^{2}+2{s_{12}}{[\widehat{\mathbf{u}}}(k)]_{1}{[\widehat{\mathbf{u}}}(k)]_{2}+{s_{22}}[\widehat{\mathbf{u}}(k)]_{2}^{2}=[\mathbf{b}]_{2}-{[\widehat{\mathbf{c}}]_{2}},$
(126)

We now define:

$\displaystyle\overline{\mu}_{1}\triangleq\frac{{{[\widehat{\mathbf{c}}]_{1}}-{[\mathbf{b}]_{1}}}}{{{\lambda_{1}}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\overline{\lambda}\triangleq\frac{{{\lambda_{2}}}}{{{\lambda_{1}}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\bar{b}\triangleq{[\mathbf{b}]_{2}}-{[\widehat{\mathbf{c}}]_{2}}.$ (127)

leveraging which, the (121) is rewritten as

$\displaystyle[\widehat{\mathbf{u}}(k)]_{1}^{2}=\overline{\mu}_{1}-\overline{\lambda}[\widehat{\mathbf{u}}(k)]_{2}^{2},$
(128)

and we can obtain from (126) that

$\displaystyle
4s_{12}^{2}[\widehat{\mathbf{u}}(k)]_{2}^{2}[\widehat{\mathbf{u}}(k)]_{1}^{2}={{\overline{b}}^{2}}+s_{11}^{2}[\widehat{\mathbf{u}}(k)]_{1}^{4}+s_{22}^{2}[\widehat{\mathbf{u}}(k)]_{2}^{4}-2\overline{b}{s_{11}}[\widehat{\mathbf{u}}(k)]_{1}^{2}-2\overline{b}{s_{22}}[\widehat{\mathbf{u}}(k)]_{2}^{2}+2{s_{11}}{s_{22}}[\widehat{\mathbf{u}}(k)]_{1}^{2}[\widehat{\mathbf{u}}(k)]_{2}^{2},$

substituting (128) into which yields

$\displaystyle{\varpi_{1}}[\widehat{\mathbf{u}}]_{2}^{4}\left(k\right)+{\varpi_{2}}[\widehat{\mathbf{u}}]_{2}^{2}\left(k\right)+{\varpi_{3}}\left(k\right)=0,$
(129)

where

$\displaystyle{\varpi_{1}}$ $\displaystyle\triangleq
s_{11}^{2}{\overline{\lambda}^{2}}+s_{22}^{2}-2{s_{11}}{s_{22}}\overline{\lambda}+4s_{12}^{2}\overline{\lambda},$
(130) $\displaystyle{\varpi_{2}}$ $\displaystyle\triangleq
2\bar{b}{s_{11}}\overline{\lambda}-2s_{11}^{2}{\overline{\mu}_{1}}\overline{\lambda}-2\bar{b}{s_{22}}+2{s_{11}}{s_{22}}{\overline{\mu}_{1}}-4s_{12}^{2}{\overline{\mu}_{1}},$
(131) $\displaystyle{\varpi_{3}}$
$\displaystyle\triangleq{\bar{b}^{2}}+s_{11}^{2}\overline{\mu}_{1}^{2}-2\bar{b}{s_{11}}{\overline{\mu}_{1}}.$
(132)

Considering $\widehat{u}_{2}^{2}(k)\geq 0$, the solution of (129) is

$\displaystyle[\widehat{\mathbf{u}}]_{2}^{2}(k)=\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}},$
(133)

substituting which into (128) yields

$\displaystyle[\widehat{\mathbf{u}}]_{1}^{2}(k)={{\overline{\mu}}_{1}}-\overline{\lambda}\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}.$
(134)

The $\widehat{\mathbf{u}}(k)$ is then straightforwardly obtained from (133)
and (134):

$\displaystyle\widehat{\mathbf{u}}(k)=\left[{\pm\sqrt{{{\overline{\mu}}_{1}}-\overline{\lambda}\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}};\pm\sqrt{\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}}}\right],$

which, in conjunction with (94) and $Q^{-1}=Q=Q^{\top}$, lead to

$\displaystyle\left[\begin{array}[]{l}\theta\left(k\right)\\\
\gamma\left(k\right)\end{array}\right]={\mathbf{Q}_{1}}\left[\begin{array}[]{l}\pm\sqrt{{{\overline{\mu}}_{1}}-\overline{\lambda}\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}}\\\
\pm\sqrt{\frac{{\sqrt{\varpi_{2}^{2}-4{\varpi_{1}}{\varpi_{3}}}-{\varpi_{2}}}}{{2{\varpi_{1}}}}}\end{array}\right].$
(139)

Substituting the notations defined in (127) into (139) and (130)–(132) results
in (99) and (100)–(102), respectively.

## Acknowledgements

This work was supported in part by National Science Foundation (award number:
CPS-1932529), Hong Kong Research Grants Council (RGC) Theme-based Research
Scheme T22-505/19-N (P0031331, RBCR, P0031259, RBCP), and RGC Germany/HK Joint
Research Scheme under Grant G-PolyU503/16.

## Author Contributions

Y.M., L.S., and H.S. designed research; Y.M. performed research; Y.M. led
experiments, H.S. conducted experiment of deep Koopman, and Y. L. collected
vehicle data and implemented self-correcting Phy-Taylor in AutoRally; Y.M.,
L.S., Q.W., and T.A. participated in research discussion and data analysis;
Y.M. and S.H. wrote this manuscript; T.A. edited this manuscript; L.S. and
T.A. led the research team.