Physical motion models offer interpretable predictions for the motion of vehicles.
However, some model parameters, such as those related to aero- and hydrodynamics, are expensive to measure and are often only roughly approximated reducing prediction accuracy.
Recurrent neural networks achieve high prediction accuracy at low cost, as they can use cheap measurements collected during routine operation of the vehicle, but their results are hard to interpret.
To precisely predict vehicle states without expensive measurements of physical parameters,
we propose a hybrid approach combining deep learning and physical motion models including a novel two-phase training procedure.
We achieve interpretability by restricting the output range of the deep neural network as part of the hybrid model, which limits the uncertainty introduced by the neural network to a known quantity.
We have evaluated our approach for the use case of ship and quadcopter motion.
The results show that our hybrid model can improve model interpretability with no decrease in accuracy compared to existing deep learning approaches.

§ INTRODUCTION

Models for multi-step prediction yield a sequence of future system states given the initial system state and a sequence of control inputs. Control methods, such as model predictive control <cit.>, rely on the predictions of future states for computing optimal control inputs given a set of objectives like heading, speed or minimal fuel consumption. Multi-step prediction models should be accurate and, in addition, they should be interpretable such that a domain expert can understand their predictions and identify their limitations, e.g. their failure states.
Interpretability is considered a crucial safety concern for ML systems by governmental entities, e.g. the European Commission <cit.>, or technical inspection and product certification services, e.g. the Technical Inspection Association (TÜV) <cit.>.
This paper develops a model for multi-step prediction, which is accurate and interpretable, to allow for future certification and usage in vehicle control.

Traditional multi-step prediction models use physical (motion) models <cit.>, which are commonly based on first-principles or linear regression. The parameters of these models correspond to physical parameters of the system, which makes these approaches interpretable. However, physical models suffer from the high cost and difficulty of measuring
parameters concerning, e.g., aero- or hydrodynamics. Although these costs can be reduced by applying simplifications, such as linearization of the model around a task-dependent operating point (e.g. a constant speed for manoeuvring <cit.>), these simplifications introduce prediction errors, which aggregate over multiple prediction steps causing the predicted states to diverge from the actual states <cit.>.

The lack of accuracy can be reduced by deep neural networks (DNN), which have shown to be capable at simulating dynamical systems over large time spans for various domains, such as unmanned aerial vehicles <cit.>, ships <cit.>, engines <cit.>, aerodynamics <cit.>, and lake temperatures <cit.>.

DNN-based models require little domain knowledge in comparison to physical models.
In addition, they are less expensive to build, as the parameters are automatically learned from the data.
However, DNNs suffer from low interpretability, as their parameters do not directly correspond to the parameters of the system and are therefore not physically meaningful (Challenge C1: Interpretability). As a consequence, they may not be used without further developments in practice, as certification bodies will not licence such technologies for use in autonomous vehicles <cit.>, at least in the current state.

In order to benefit from both approaches, hybrid models combine physical models and DNN, while minimizing the respective disadvantages. Improved interpretability compared to deep learning is achieved by partially retaining the interpretability of the physical model via certain structures in the model architecture or by constraining the dynamics of the DNN.
Recent approaches have employed a residual hybrid architecture, in which the DNN is trained to compensate for the error residual of the physical model <cit.>.
<cit.> and <cit.> achieve a higher interpretability compared to DNN due to the combination of physical model and DNN output for the overall prediction. <cit.> feeds the output of an additional DNN as input to the physical model. This has a negative impact on the interpretability compared to <cit.> and <cit.>, since the input-output relation of the physical model is obscured by the DNN.
However, the approach allows the inclusion of unknown parameters, which are identified during training of the DNN, which makes it more flexible than <cit.> and <cit.>.
However, all these three approaches only consider purely linear models as physical models, which is a considerable limitation, since most vehicle dynamics are non-linear (Challenge C2: Non-linear physical models).

The training of hybrid models that use non-linear or inaccurate physical models for multi-step prediction is difficult (Challenge C3: Trainability of hybrid models).
Due to feedback of the predicted state as input to the physical model,
errors in the prediction aggregate over time, which causes divergence in the weights of the DNN.

Another limitation of the existing work is its evaluation in use cases with low environmental disturbances, such as quadcopters operated indoors or unmanned surface vehicles operated on lakes <cit.>.
Environmental disturbances have a large effect on vehicle dynamics due to strong coupling between states (Challenge C4: Environmental disturbances).

To address these challenges, we propose a novel hybrid method for multi-step prediction of vehicle states consisting of a residual hybrid architecture with a respective training procedure. Our proposed hybrid architecture combines long short-term memory networks (LSTM) and physical models similar to <cit.>, but with the capacity to model a large variety of physical models composed of first-principles and regression models (C2).
Using our two-phase training procedure, we address the existing problem of divergence during training of hybrid models in multi-step prediction, thus allowing for inaccurate and non-linear physical models as part of the hybrid model (C3).
Our procedure ensures convergence of the DNN by training it to correct the single-step error of the physical model and the aggregated error over multiple prediction steps.
In addition, we guarantee interpretability of our approach by applying a constraint on the output range of the LSTM (C1). This enables a trade-off between interpretability of the hybrid model and accuracy by varying the strictness of the constraint. Finally, we evaluate our approach on the use-case of ship motion under environmental disturbances (C4) and on the use-case of quadcopter motion without disturbances.

§ RELATED WORK

We have identified three families of related approaches that merge knowledge of a physical system with deep learning.

Network architectures.
Deep feedforward neural networks (FNN) are frequently used to approximate complex non-linear dynamics <cit.>.
These FNN are known as non-linear auto-regressive exogenous (NARX) networks.
They are given a fixed window of prior states and control inputs and predict subsequent states.
Because the size of their input window is fixed, NARX networks do not cope well with long-term dynamics.

Recurrent neural networks (RNN) are a superior alternative to FNNs for simulating dynamical systems, as their recurrent structure permits them to model complex non-linear dynamics over large time spans <cit.>.

The hidden state of RNN can learn latent temporal dependencies, e.g. environmental effects such as waves, which are not directly measurable with sensors.
Thus, RNNs may provide more accurate state predictions than FNNs.
A network architecture that models the physical processes of the system can also improve prediction accuracy.
Schürholz et al. <cit.> employs a recurrent architecture with additional forward connections between  recurrent units, which correspond to the physical information flow between components further improving the prediction accuracy.

Physics-based loss functions.
Loss functions can encode known dynamics and physical parameters as regularization terms.
In <cit.> and <cit.> physics-based restrictions motivated by the law of energy conservation are used on the loss function to simulate lake temperatures with FNNs and RNNs.
Their evaluations show that this method allows the network to converge faster than a network with unmodified loss.
Approaches dealing with object tracking  <cit.> employ loss functions for semi-supervised learning by encoding dynamics in the loss function and punishing trivial solutions.
Including knowledge about the physical system in loss functions reduces the amount of data and time required for training the DNN and results in physically plausible predictions.
However, physics-based loss functions do not improve interpretability, as they only affect the training of the neural network.

Hybrid models.
Hybrid models combine different types of models, such as physical models and DNN.
In a serial architecture, approaches like <cit.> and <cit.> either give the output of a DNN as input to a physical model or vice versa.
In either variant, the output of the second model is returned as overall prediction.
The rationale underlying these architectures is that certain inputs to a physical model, such as hydrodynamic forces, are difficult to model with first-principles, but can be predicted by a DNN from control inputs and the system state.
The serial connection limits the interpretability of the hybrid model, as the input-output relation is obscured by the neural network.
On the other hand, a parallel architecture adds the output of a physical model and the output layer of a DNN into an overall prediction <cit.>.
This approach retains the interpretability of the physical model, if the DNN contribution is sufficiently small in comparison.
Recent work uses residual architectures, which employ serial and parallel connections in joint <cit.>.
The serial connection, which gives the output of the physical model as input to the DNN, allows the network to rectify errors of the physical model and improve accuracy over multiple prediction steps. The parallel connection sums the error residual predicted by the DNN and the physical model output for the overall prediction thereby ensuring the interpretability of the prediction.
Lastly, current approaches utilize physical models, which are either relatively accurate but expensive <cit.> or use simplifications such as linear models <cit.>.
The former type of model is expensive to develop, while the latter limits the possible parameters present in the physical model.

§ HYBRID MODEL AND TRAINING FOR MULTISTEP PREDICTION

Our proposed approach consists of a hybrid architecture and a two-phase training procedure.
We combine a physical model with an LSTM, which predicts the error residual of the physical model. The physical model and LSTM are described in Sections <ref> and <ref>, respectively.
The goal of the two-phase training process is to achieve convergence.
In the first phase, we reduce the one-step prediction error by teacher-forcing <cit.>, i.e. at each time step, the hybrid model receives the true prior state as input.
In the second phase, we train with multistep prediction feeding back each predicted state as input for the next prediction.
The variables used for the model definition are summarized in Table <ref>.

Summary of variables and physical parameters used in the model.
Matrices are denoted by bold uppercase letters.
Vectors are denoted by lowercase letters with upper arrow.
Scalars are denoted by lowercase letters.

Variable      3lDescription

$\vect{\eta}$     3lposition/attitude vector in inertial frame with

$x$     x-coordinate
$y$     y-coordinate

$\phi$     roll angle
$\psi$     yaw angle

$\vect{v}$     3lvelocity vector in body frame with

$u$     surge velocity (x-axis)
$w$     sway velocity (y-axis)

$p$     roll rate (x-axis)
$r$     yaw rate (z-axis)

$\matr{M}$     3lmass and inertia matrix

$\matr{D}(\vect{v})$     3lhydrodynamic damping matrix

$\matr{C}_\rb(\vect{v})$     3lrigid-body Coriolis and centripetal matrix

$\vect{g}(\vect{\eta})$     3lrestoring forces

$\force{control}$     3lcontrol forces (propulsion, steering)

$\force{env}$     3lenvironmental forces

$\matr{J}(\vect{\eta})$     3lrotation from body to inertial frame

$\vect{z}_t$     3lmeasured system state at time step $t$

3l$= [u, w, p, r, \phi]$

$\vect{c}_t$     3lcontrol input at time step $t$

§.§ Overall Architecture

As shown in Figure <ref>, the hybrid model consists of an LSTM and a physical model. The physical model can be a regression or first-principles component, as well as a combination of both. A first-principles component is a set of differential equations derived directly from Newton's second law of motion. A regression component allows learning of an interpretable model from data, which is useful estimating hydrodynamical coefficients <cit.>. It is common practice to combine both approaches to yield more accurate physical models.

The physical model $\vect{z}_{t+1}^{\;\text{phy}}$ computes its prediction from the control input $\vect{c}_t$ and the previous state $\vect{z}_{t}$.
The LSTM $\vect{z}^{\;\text{lstm}}_{t+1}$ predicts the error residual of the physical model  given the previous hidden state $\vect{h}_{t}$, the control input $\vect{c}_t$, and the current state prediction of the physical model $\vect{z}_{t+1}^{\;\text{phy}}$.
The previous state $\vect{z}_{t}$ is not needed as input, since it is composed of the prior physical model output $\vect{z}_{t}^{\;\text{phy}}$, which was provided as input, and its own output $\vect{z}^{\;\text{lstm}}_{t}$.
Accordingly, this information is already encoded in the hidden state and its inclusion yielded no performance improvements in our initial experiments.
The state prediction of the hybrid model $\hat{\vect{z}}_{t+1}$ is the sum of $\vect{z}_{t+1}^{\;\text{phy}}$ and $\vect{z}^{\;\text{lstm}}_{t+1}$.
During multi-step prediction, the auto-regressive loop, represented by the dashed arrow, feeds back this state prediction $\hat{\vect{z}}_{t+1}$ as input to the physical model.
Because of the feedback loop, the multi-step prediction error may aggregate over time.
The forward computation for our hybrid model for the time step $t+1$ is defined as follows:

\begin{align}
\begin{split}
\vect{z}^{\;\text{phy}}_{t+1} &= \text{Regression}(\vect{c}_{t}, \vect{z}_{t}) + \text{1st-Principles}(\vect{c}_{t}, \vect{z}_{t})
\end{split} \label{eq:physical_true}\\
\vect{h}_{t+1} &=  \text{LSTM}(\vect{c}_t, \vect{z}^{\;\text{phy}}_{t+1}, \vect{h}_{t}) \\
\vect{z}^{\;\text{lstm}}_{t+1} &= \matr{W}^{\text{hx}} \, \vect{h}_{t+1} \\
\begin{split}
\hat{\vect{z}}_{t+1} &= \vect{z}^{\;\text{phy}}_{t+1} + \vect{z}^{\;\text{lstm}}_{t+1} \\
\end{split}
\end{align}

[>=stealth,->,line width=0.75pt]

[matrix of nodes, column sep=0.75em, align=center, nodes=rectangle, anchor=center]
[value] (state) $\vect{z}_{t}$;
[concat] (cat_input) ;
[weight] (physical) 1st-Principles;
[operator] (sum_whitebox) $+$;

[operator] (sum_all) $+$;
[value] (state_est)   $\hat{\vect{z}}_{t+1}$;

[value] (control) $\vect{c}_t$;

[weight] (semiphysical) Regression;

[concat] (cat_blackbox) ;
[weight] (blackbox) LSTM;

(control) to [out=0,in=270,loop,looseness=0.75] (cat_input);
(state) to (cat_input);
[to] (cat_input) – (physical);
(cat_input) to [out=0,in=180,loop,looseness=0.75]  (semiphysical);
(control) to [out=270,in=270,loop,looseness=0.5] (cat_blackbox);
[to] (cat_blackbox) – (blackbox);
(semiphysical) to [out=0,in=270,loop,looseness=0.75] (sum_whitebox);
[to] (physical) –  (sum_whitebox);
(sum_whitebox) to [out=0,in=90,loop,looseness=0.75]  (cat_blackbox);
[to] (sum_whitebox) – (sum_all);
(blackbox) to [out=0,in=270,loop,looseness=0.75]  (sum_all);
[to] (sum_all) – (state_est);

(blackbox) to [out=315,in=225,loop,looseness=2] (blackbox) node[xshift=0.3em,yshift=-2.6em] $\vect{h}_{t}$;
[dashed] (state_est) to [out=90, in=90, loop, looseness=0.2] (state) node[fill=white,midway,yshift=4em] $t := t + 1$;

The proposed hybrid architecture consists of physical model (regression and first-principles) and LSTM.
Concatenation of vectors is represented by a black circle. The $+$-operators corresponds to an element-wise addition.

Training. The training of the hybrid approach is divided into training of the physical model and training of the LSTM.
The regression component is fit on the residual of the first-principles model for the single-step objective:

\begin{align}
\begin{split}
\mathcal{L}^{\text{reg}}(\vect{\theta}) =& \text{MSE}(\vect{z}_{t+1} - \text{1st-Principles}(\vect{c}_{t}, \vect{z}_{t}),\\
&\phantom{\text{MSE}(}\text{Regression}(\vect{c}_t, \vect{z}_{t}; \vect{\theta}) )
\end{split}
\end{align}

If no first-principles component is used, the regression fits on the task of predicting the system state instead of the first-principles residual.

The training of the LSTM on the multi-step prediction objective, i.e. predicting the state over multiple time steps given a sequence of control inputs, is difficult, because each next-step prediction of the physical model depends on the inaccurate previous prediction. Thus, the error aggregates over several time steps causing divergence between predicted state and true state.
LSTMs are trained with backpropagation-through-time (BPTT), which is sensitive to large errors.
A combination of increasingly large errors and BPTT causes exploding gradients and prevents successful training of the network.

We solve this problem with our two-phase training procedure, which enables training hybrid models with divergent and inaccurate physical models.

First phase.
In the first phase, we employ teacher forcing, where the physical model receives the true state $\vect{z}_{t}$ as input at each time step $t+1$.
This mechanism enables the LSTM to learn the one-step error of the physical model, which occurs due to missing parameters or simplifications. The loss function for the first phase is defined as:

\begin{align}
\hspace*{-0.8em}\mathcal{L}^\text{lstm}(\vect{\theta}) = \sum_{t=1}^{T} \text{MSE}(\vect{z}_{t+1} - \vect{z}^{\;\text{phy}}_{t+1}, \text{LSTM}(\vect{c}_t, \vect{z}^{\;\text{phy}}_{t+1}, \vect{h}_{t}; \vect{\theta}))
\end{align}

where $T$ is the number of time steps and MSE is the mean squared error.

The first training phase is only performed to reduce the estimation error until feedback of the state prediction does not cause divergence during multi-step prediction. Therefore the first-phase loss function is only employed for a limited number of epochs.

Second phase.
In this phase, the model learns to predict states over multiple time steps without receiving true states as input. This is achieved by using the state prediction from the previous time step as input in the current prediction step.
The training remains effective, because the first-phase ensured that the hybrid model achieves sufficiently accurate predictions to prevent the strong divergence of DNN gradients.
The loss function for the second phase only differs from the first-phase in the physical model prediction $\vect{z}^{\;\text{phy}}_{t+1}$, as the prediction now depends on the predicted prior state $\hat{\vect{z}}_t$ and not on the true prior state:

\begin{align}
\vect{z}^{\;\text{phy}}_{t+1} &= \text{1st-Principles}(\vect{c}_t, \hat{\vect{z}}_{t}) + \text{Regression}(\vect{c}_t, \hat{\vect{z}}_{t})
\end{align}

§.§ Physical Models

We define (cf. <cit.>) physical models for ships combining a first-principles model with a regression model.
The simple dynamics are modelled precisely using the corresponding physical equations. The complex system dynamics are approximated using least squares regression.

First-principles. The models of motion are represented by differential equations, which define the equilibrium of forces resulting in the acceleration of a vehicle.
We use the motion model for surface vessels developed by Fossen <cit.>:

\begin{align}
\deriv{\eta} &= \matr{J}(\vect{\eta}) \vect{v} \label{eq:kinematics} \\
\deriv{v} &= \matr{M}^{{-1}} (\force{control} + \force{env} - \matr{D}(\vect{v})\vect{v}  - \matr{C}_\rb(\vect{v})\vect{v} - \vect{g}(\vect{\eta}))
\label{eq:kinetics}
\end{align}

Equation <ref> defines the kinematics of the vessel, as it describes the mapping of velocities in the ship's body frame to a motion in the inertial reference frame, i.e. earth.
The ship's kinetics represented by Equation <ref> define the relation of forces affecting the vessel and the resulting acceleration.

Developing first-principles models is challenging due to the high cost and difficulty associated with measuring their physical parameters.
Additionally, environmental variables, such as the sea state, cannot be measured during the operation of a vehicle, which makes their prediction with first-principles models infeasible.
Due to these limitations, we consider two examples of incomplete first-principles models, which entirely omit immeasurable parameters.
First, the minimal model () uses only parameters that are easy to measure or approximate, such as those related to the vessel's mass, inertia and geometry. The restoring forces are linearized as $\matr{G}\vect{\eta}$ around a roll angle of $0$ yielding a good approximation with fewer parameters <cit.>:

\begin{align}
\deriv{v} &= \matr{M}^{{-1}} (-\matr{C}_\rb(\vect{v})\vect{v} - \matr{G} \vect{\eta}))
\end{align}

All the other parameters that are expensive to measure are set to $0$, e.g. ones concerning hydrodynamic forces.

Second, a propulsion-based model (), which extends the minimal model by parameters for propulsion and steering  forces.
This model includes $\force{control}$, which represents the control forces generated by propellers and rudders:

\begin{align}
\deriv{v} &= \matr{M}^{{-1}} (\force{control}(\vect{c}, \vect{v})-\matr{C}_\rb(\vect{v})\vect{v} - \matr{G} \vect{\eta}))
\end{align}

The model for control forces $\force{control}$ is highly dependent on the propulsion system of the modeled vessel, as it depends on the thruster configuration and the type of steering mechanism. The implementation details of a standard propeller and rudder may be found in the MSS toolbox <cit.>.

Linear regression. We employ linear regression in two models components as a coarsely approximating physical model <cit.>.
The linear regression model () receives  prior system state $\vect{z}_{t}$ and control input $\vect{c}_{t}$:

\begin{align}
\hat{\vect{z}}_{t+1} &= \matr{Q} \vect{c}_{t} + \matr{P} \vect{z}_{t} + \vect{b}\label{eq:linear}
\end{align}

with input matrix $\matr{Q}$, system matrix $\matr{P}$, and bias vector $\vect{b}$.

Second, we extend the linear model with a non-linear input mapping $f$ to better approximate hydrodynamic forces:

\begin{align}
\hat{\vect{z}}_{t+1} &= \matr{Q} \vect{c}_{t} +\matr{P} f(\vect{z}_{t}) + \vect{b}
\end{align}

We use the set of non-linear features developed by <cit.> consisting of polynomial and absolute terms, which are well-suited for modeling the damping and drag effects of a ship hull with the surrounding fluid.

We adapt their features for our hydrodynamic regression model () as follows:

\begin{align}
u_{t+1} &= \vect{a}_{u}^\intercal \, [ u_{t} \hspace{0.4em} \abs{u_{t}}u_{t} \hspace{0.4em} w_{t}r_{t}]^\intercal \\
v_{t+1} &= \vect{a}_{w}^\intercal \, [
w_{t} \hspace{0.4em}
\abs{u_{t}}w_{t} \hspace{0.4em}
u_{t}r_{t} \hspace{0.4em}
\abs{w_{t}}w_{t} \hspace{0.4em}
\abs{r_{t}}w_{t} \hspace{0.4em}
\abs{w_{t}}r_{t} \hspace{0.4em}
\abs{u_{t}w_{t}}\phi \hspace{0.4em}
u_{t}^2\phi ]^\intercal\\
p_{t+1} &= 0\\
r_{t+1} &= \vect{a}_{r}^\intercal \, [ r_{t} \ \abs{u_{t}}w_{t} \hspace{0.4em} \abs{w_{t}}r_{t} \hspace{0.4em} \abs{u_{t}}r_{t} \hspace{0.4em} \abs{r_{t}}r_{t}\hspace{0.4em} \abs{u_{t}\phi_{t}}\phi_{t} \hspace{0.4em} \abs{r_{t}}u_{t}\phi_{t} \hspace{0.4em} \abs{u_{t}}u_{t-1}\phi_{t} ]^\intercal\\
\phi_{t+1} &= \vect{a}_{\phi}^\intercal \, [ p_{t} \hspace{0.4em} \phi_{t} ]^\intercal
\end{align}

where weights for each state are encoded by the vectors $\vect{a}_u$, $\vect{a}_w$, $\vect{a}_p$, $\vect{a}_r$, and $\vect{a}_\phi$. The vertical bars, e.g. $\abs{u_t}$, denote the absolute value of a scalar.

§.§ Recurrent Neural Network

We employ a stacked LSTM including an initialization mechanism as DNN component in our hybrid model, as it is capable of memorizing long-term dynamics and of modeling dynamics at different time scales <cit.>. The latter capability is important for modeling ships, since waves and control inputs introduce dynamics with different frequencies.

The models for multi-step prediction, such as a LSTMs, rely on an initial state to yield accurate predictions.
Thus, the initial state should encode all available information about the system and the surrounding environment, e.g. the velocity and position of a ship but also the sea state.
The state of an LSTM is its cell and hidden state. An LSTM can be initialized at $t=0$ by providing a cell and hidden state.
Mohajerin et al. <cit.> propose the initialization with an additional LSTM acting as initializer.
The initializer receives as input a window of system states and control inputs, which precede the predicted time frame.
The recurrence of the initializer enables it to extract latent dynamics from the data, which cannot be directly measured such as the sea state.
The prediction model is then initialized with the final hidden state of the initializer.

§ EVALUATION

The proposed method is evaluated extensively for modeling dynamics of a patrol ship under environmental disturbances (Section <ref> and <ref>). To prove generalizability Section <ref> summarizes results for quadcopter motion.

§.§ Experiment Setup

Dataset. The dataset is generated with a 4 degrees-of-freedom (DOF) maneuvering model of a patrol ship, which simulates the horizontal plane motion and rolling due to waves.
The original model and physical parameters are provided by Perez et al. <cit.> and can be found in the MSS toolbox <cit.>.
The model is extended to include wind forces according to Isherwood <cit.>, ocean waves generated with the JONSWAP spectrum <cit.>, and corresponding wave forces computed with force response amplitude operators.
Two symmetrically-placed fixed-pitch rudder propellers are implemented for control of heading and speed of the vessel.
Control inputs to the simulation model are generated via an open-control loop, which mimic a human operator.
The measurements of the vessel are sampled with a rate of 1.
The dataset consists of $96$ hours separated into $\SI{1}{\hour}$ simulations, where each simulation is initialized with a random sea state and control inputs.
The dataset is split into training, validation, and test set with a 60-10-30 split.

Metrics. The prediction performance is measured with regards to state and trajectory prediction over a time span of $\SI{15}{\minute}$ or $900$ time steps.
A window containing $\SI{1}{\minute}$ of prior control inputs and state variables is provided as initial state.
The five predicted states are the velocities for surge $u$, sway $w$, roll $p$, and yaw $r$, as well as the roll angle $\phi$.
The root mean squared error (RMSE) is averaged over all time steps for each state to measure the state prediction performance.
The trajectory RMSE measures the average distance between the predicted and true trajectory for each time step.
The trajectory error emphasizes surge, sway, and yaw, as they have the largest impact on the trajectory.
Interpretability is measured using an output range constraint enforced on the LSTM. An unconstrained LSTM results in no interpretability, while an increasingly smaller output range increases interpretability. We employ the relative threshold applied to the LSTM as measure of interpretability. The relative threshold is the ratio between the sizes of the LSTM output range and the expected output range of the physical model. The expected output range of the physical model is the $95\%$ range of its outputs computed on the training and validation set. For example, a physical model has an expected range of $[-5, 5]$ and the LSTM is constrained to a range of $[-1, 1]$, then the relative threshold is $20\%$. A lower relative threshold implies a higher degree of interpretability.

Model selection. Hyperparameter search is performed via grid search on the validation set.
The best model per class is chosen based on the RMSE summed over each state.
Quadratic control lag () is a regression model developed by <cit.> as best-performing physical baseline in their work.
It uses a window of prior control inputs and states, and includes control inputs as quadratic terms to model the relation between propulsion and velocity.
Stacked LSTM with an LSTM initializer (, <cit.>) is used as DNN baseline.
Hybrid models are trained with each potential physical model configuration using the default one-phase training and the proposed two-phase training ({, }+{, }-{, }).
This allows us to evaluate whether our two-phase training can successfully train hybrid models with incomplete and non-linear physical models.
The models are denoted by abbreviations as defined in the previous text.
is similar to the hybrid model by <cit.>.

§.§ Results

State prediction.
Table <ref> shows the state prediction performance of the baselines  and , as well as each variation of the hybrid model. Error values across different state variables are not directly comparable, as they are not normalized and each variable has a different magnitude.
is the only hybrid model without two-phase training that outperforms its counterpart with two-phase training. However,  performs worse overall compared to all other two-phase hybrid models.
For all other configurations, the hybrid model with one-phase training performs worse by a large factor.
For , no model was identified to be capable of simulating the vessel without diverging completely.
These results indicate that the proposed two-step training procedure allows integration of a large range of physical models with DNN compared to the default training procedure.
achieves the lowest error for three state variables. Additionally, it achieves a trajectory error comparable to  as their confidence intervals overlap. We therefore consider  as the best performing model in our evaluation. Overall, most hybrid models achieve performance comparable to or better than the DNN baseline . Hybrid models that employ  show better results on average. We hypothesize that the resulting loss function for the LSTM is less complex due to the linearity of  compared to more complex non-linear physical models and therefore easier to learn.

Figure <ref> visualizes a single test sample of predicted and true states over a time span of $\SI{15}{\minute}$.
Performance between  and  are similar and outperform  and .
Prediction of surge speed $u$ is very accurate over 15 for , , . It is the easiest state to predict, since the effect of propulsion outweighs all other dynamics. However, sway $w$ and yaw rate $r$ are non-linearly coupled and strongly influenced by external forces, which makes them to difficult to predict. The impact of waves on the yaw rate $r$ is observable in the oscillations of the predicted and true state.
and  successfully predict the mean yaw rate, but fail to predict the wave-induced oscillations due to its randomness.
has the best performance predicting roll rate $p$. However, the time series plot shows that the prediction is not accurate and tends towards predicting a value of $0$.
Accordingly, none of the evaluated methods are suited for use-cases with strong waves, which can be attributed to the randomness of ocean waves.

Trajectory prediction.
The results for trajectory prediction are summarized in the last of column of Table <ref>.
The errors mirror the results for state prediction.
Hybrid models with two-phase training and  achieve similar trajectory errors.
The hybrid models ,  and  perform significantly better than all other models.
Two-phase training provides better trajectory prediction compared to one-phase training for all cases but .
Figure <ref> shows boxplots for the trajectory error per minute for the baseline model  and the best-performing hybrid model .
Our hybrid approach achieves a lower mean and standard deviation in the trajectory error than the . Accordingly, the approach is more reliable for trajectory prediction over long time frames than a pure deep learning method.

An example of a trajectory prediction is visualized in Figure <ref>.
The trajectories are computed from the state predictions in Figure <ref>.
and  predict similar trajectories, which are very close to the true trajectory.
The divergence from the true trajectory in each model occurs due to the difficulty of predicting yaw rate $r$, i.e. rotation around the z-axis, since the divergence occurs as a result of changing the heading of the vessel.
This difficulty can be observed specifically for , as it underpredicts the yaw rate $r$, as shown in Figure <ref> starting at approx. 300, and consequently diverges from the true trajectory.

We apply a threshold to the output range of the LSTM to create an interpretable hybrid model.
Figure <ref> shows the RMSE of our best-performing hybrid model  for surge $u$, sway $v$ and yaw rate $r$ for various relative thresholds in orange. The blue horizontal line indicates the error of the model without constraints.
The figure shows that a relative threshold of ${\sim} 15\%$ is sufficient to achieve prediction accuracy similar to the corresponding unconstrained model for all three state variables. Therefore, the loss in interpretability due to the LSTM can be limited effectively. The LSTM applies small adjustments to correct errors introduced by the physical model, which matches the training objective.

In summary,  outperforms all other models including  and the other hybrid model variants w.r.t. ship state and trajectory prediction.
We show that  can be constrained effectively to achieve model interpretability with no loss in accuracy. Accordingly, our hybrid approach enables interpretable and accurate multistep predictions in the use of case of predicting ship motion.

§.§ Generalizability

To show generalizability of our approach, we summarize our results for the task of quadcopter motion with the dataset used in <cit.>.

Dataset. The dataset by <cit.> consists of multiple acrobatic maneuvers performed by an expert pilot in an indoor environment. The sampling rate is $\delta = \SI{100}{\hertz}$. Compared to the ship motion use case, measurements are more precise and frequent due to lack of environmental disturbances and use of a motion capture system.
Unlike a ship, motion of the quadcopter is mainly determined by its propulsion, as its mass and surface area are low and therefore inertial and aerodynamic forces are small. Due to the performed aerobatics, the control inputs are irregular unlike the ship use case, where control inputs were generated to simulate routine operation.
The state vector is defined as $\vect{z} = \begin{bmatrix}\dot{x} & \dot{y} & \dot{z} & p & q & r\end{bmatrix}$, which are velocities along each axis in the inertial frame and the angular rates for roll, pitch, and yaw.

Quadcopter Models. We perform a train-validation-test split of 60-10-30 and use grid-search to optimize hyperparameters.
As suggested by <cit.>, the evaluated prediction horizon is 1 corresponding to $100$ time steps. Similar to , we develop a minimal physical model  for quadcopters, which models rigid-body forces and gravity. We reuse the linear regression model  (Equation <ref>). We introduce a quadratic regression model , which includes the square of each control input and each state to more effectively represent forces, such as air resistance. <cit.> and <cit.> are again used as baseline models.

Results. Table <ref> shows state and trajectory prediction performance of the baseline models and each hybrid model using two-phase training for the quadcopter dataset.  outperforms all other models in state and trajectory prediction. All hybrid models achieve similar results and outperform . The best performing hybrid model is . Figure <ref> shows that  constrained to a relative threshold of ${\sim} 8\%$ achieves comparable results to the unconstrained model.
The results confirm that our hybrid models are preferable w.r.t. accuracy over physical models (). The hybrid models perform slightly worse than pure deep learning for the quadcopter motion, but allows more interpretable predictions.

§ CONCLUSIONS

In this paper, we developed a residual hybrid model including a two-phase training procedure. We define relative threshold as a metric for interpretability and employ constraints on the LSTM to achieve interpretable predictions.

Our approach is more flexible than existing work as it is capable of employing a large variety of physical models, such as models with missing parameters and non-linearities (C2: Non-linear physical models). Our results indicate that linear models are preferable as physical models w.r.t. accuracy, as they achieve the best results for ship motion prediction.
The proposed two-phase training procedure enables accurate predictions over large time spans for various physical models, while the default training approach either fails to train the network or results in worse prediction accuracy (C3: Trainability of hybrid models).
The developed hybrid models outperform pure deep learning for ship motion with environmental disturbances (C4: Environmental disturbances) and performs only slightly worse for quadcopter motion.

The licensing of our method as part of a control system requires that the DNN contribution has to be low to ensure interpretability. Additionally, failure conditions of the DNN have to be understood to ensure safety. While we achieve the former requirement with a constraint on the DNN output range without loss of accuracy (C1: Interpretability), the latter one requires future research that answers what dynamics are learned by the DNN.

Additionally, our current approach only uses a simple constraint that applies to the output. Mathematical properties desirable in models for control are, for example, Lipschitz continuity and bounded-input, bounded-output stability, as they allow integration within robust model-based control frameworks. These properties require that a relation between input and output of the model is enforced. Future work will extend our hybrid approach with such model properties to integrate with model-based control frameworks.

Acknowledgments. Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2075 - 390740016. We acknowledge the support by the Stuttgart Center for Simulation Science (SimTech).


# Dissipative soliton interaction in Kerr resonators with high-order
dispersion

A. G. Vladimirova, M. Tlidib, and M. Takic aWeierstrass Institute,
Mohrenstrasse 39, 10117 Berlin, Germany  bDépartement de Physique, Faculté des
Sciences, Université Libre de Bruxelles (U.L.B.), CP 231, Campus Plaine,
B-1050 Bruxelles, Belgium  cUniversité de Lille, CNRS, UMR 8523 - PhLAM -
Physique des Lasers Atomes et Molécules, F-59000 Lille, France

###### Abstract

We consider an optical resonator containing a photonic crystal fiber and
driven coherently by an injected beam. This device is described by a
generalized Lugiato-Lefever equation with fourth order dispersion. We use an
asymptotic approach to derive interaction equations governing the slow time
evolution of the coordinates of two interacting dissipative solitons. We show
that Cherenkov radiation induced by positive fourth-order dispersion leads to
a strong increase of the interaction force between the solitons. As a
consequence, large number of equidistant soliton bound states in the phase
space of the interaction equations can be stabilized. We show that the
presence of even small spectral filtering not only dampens the Cherenkov
radiation at the soliton tails and reduces the interaction strength, but can
also affect the the bound state stability.

##  I Introduction

. Optical frequency combs generated by micro-cavity resonators have
revolutionized many fields of science and technology, such as high-precision
spectroscopy, metrology, and photonic analog-to-digital conversion (Fortier
and Baumann, 2019). A particular interest is paid to the soliton frequency
combs associated with the formation in the time domain of the so-called
temporal cavity solitons – nonlinear localized structures of light, which
preserve their shape in the course of propagation. Temporal dissipative
solitons often called cavity solitons were reported experimentally in mode-
locked lasers, micro-cavity resonators (Kippenberg et al., 2018; Herr et al.,
2014), and in coherently driven fiber cavities (Leo et al., 20103).

In this work we consider a photonic crystal fiber cavity driven by a coherent
injected beam. When operating close to the zero dispersion wavelength, high-
order chromatic dispersion effects could play an important role in the
dynamics of the system. Taking into account these effects together with
spectral filtering the dimensionless model equation in the mean-field limit
reads

$\partial_{t}U=S-(1+i\theta)U+iU|U|^{2}+\left(\delta+i\right)\partial^{2}_{\tau}U+\beta_{3}\partial^{3}_{\tau}U+i\beta_{4}\partial^{4}_{\tau}U,$
(1)

where $U(\tau,t)$ is the complex electric field envelope, $\tau$ is time and
$t$ is the slow time variable describing the number of round trips in the
cavity. The parameter $S$ measures the injection rate, $\theta$ describes
frequency detuning, second order dispersion and Kerr nonlinearity coefficients
are normalized to unity, $\beta_{3}$ and $\beta_{4}$ are the third and fourth-
order dispersion coefficients, respectively, and $0<\delta\ll 1$ is the small
spectral filtering coefficient. The optical losses are determined by the
mirror transmission and the intrinsic material absorption. This losses are
normalized to unity.

In the absence of high-order dispersion and spectral filtering, we recover
from Eq. (1) the Lugiato-Lefever equation (Lugiato and Lefever, 1987) which is
a paradigmatic model to study temporal cavity solitons (see overview (Chembo
et al., 2017; Lugiato et al., 2018)). It is widely applied to describe two
important physical systems: passive ring fiber cavity with coherent optical
injection and driven optical microcavity used for frequency comb generation
(Haelterman et al., 1992; Maleki et al., 2010; Matsko et al., 2011; Chembo and
Menyuk, 2013). The inclusion of the fourth-order dispersion allows the
modulational instability to have a finite domain of existence delimited by two
pump power values (Tlidi et al., 2007). As a consequence, upper homogeneous
steady state solution becomes modulationally stable and dark dissipative
solitons sitting this solution can appear (Tlidi and Gelens, 2010). In the
presence of third order dispersion bright and dark dissipative solitons become
asymmetric and acquire an additional group velocity shift associated with this
asymmetry (Akhmediev et al., 1990; Tlidi et al., 2013; Milián and Skryabin,
2014; Vladimirov et al., 2018).

Being well separated from one another dissipative solitons can interact via
their exponentially decaying tails and form bound states characterized by
fixed distances between the solitons. This weak interaction can be strongly
affected by different perturbations, such as periodic modulation (Soto-Crespo
et al., 2003; Turaev et al., 2012) and high-order dispersions (Olivier et al.,
2006), which lead to the appearance of the so-called soliton Cherenkov
radiation at the soliton tails (Akhmediev and Karlsson, 1995). Soliton
interaction in the presence of high-order dispersions was theoretically
studied in several works (Olivier et al., 2006; Milián and Skryabin, 2014;
Parra-Rivas et al., 2016, 2017; Vladimirov et al., 2018). However, they were
either focused on the asymmetric soliton interaction in the presence of third-
order dispersion or based mainly on the numerical calculation of the soliton
interaction potential. Unlike these works, here we present an analytical
theory of the interaction of two dissipative solitons of the Lugiato-Lefever
equation with fourth-order dispersion term based on the asymptotic approach
developed in (Gorshkov and Ostrovsky, 1981; Karpman and Solov’ev, 1981).
Furthermore, we show that similarly to the case of the interacting oscillatory
solitons (Turaev et al., 2012), a small spectral filtering effect can strongly
affect the interaction force and the stability properties of the bound soliton
states.

##  II Single peak dissipative soliton

Without high-order dispersion and spectral filtering terms,
$\beta_{3}=\beta_{4}=\delta=0$, Eq. (1) supports a single or multipeak
dissipative solitons characterized by damped oscillatory tails (Scroggie et
al., 1994). Stable dissipative solitons have been found in a strongly
nonlinear regime, where the modulational instability is subcritical, i.e., for
$\theta>41/30$. More precisely, they have been found in the pinning region,
where the lower stationary homogeneous solution coexist with a periodic one.
The number of dissipative solitons and their distribution in the cavity are
determined by the initial conditions while their maximum peak power remains
constant for fixed values of the system parameters (Scroggie et al., 1994).
For $\theta>41/30$ Eq. (1) supports a single peak dissipative soliton solution
in the form $U(t,\tau)=U_{0}+u_{0}(\tau)$, where $I_{0}=|U_{0}|^{2}=const$ is
the intensity of the stationary homogeneous solution of Eq. (1) and
$u_{0}(\tau)$ decays exponentially at $\tau\to\pm\infty$. This solution
persists also at sufficiently small $\beta_{3}$, $\beta_{4}$, and $\delta$. It
remains motionless for $\beta_{3}=0$ and becomes uniformly moving otherwise,
$U(t,\tau)=U_{0}+u_{0}(\tau-vt)$. Asymptotic analytic theory of the asymmetric
dissipative soliton interaction via Cherenkov radiation induced by the third-
order dispersion coefficient $\beta_{3}$ was developed in (Vladimirov et al.,
2018). Below we consider the case when only small fourth-order dispersion is
present, $\beta_{3}=0$ and $|\beta_{4}|\ll 1$. In this case due to the
symmetry property of Eq. (1), $\tau\to-\tau$, the soliton velocity is always
zero, $v=0$.

The dispersion relation for the small amplitude waves is determined by
substituting $U(t,\tau)=U_{0}+A_{0}e^{ik\tau-i\Lambda t}$ into Eq. (1) and
linearizing the resulting equation at $U=U_{0}$. This yields

$\Lambda=-2I_{0}^{2}+i\sqrt{\left(1+\delta
k^{2}\right)^{2}-I_{0}^{2}}+k^{2}-\beta_{4}k^{4}.$

The phase velocity of the dispersive waves $V=\Re(\Lambda)/k$ is shown in Fig.
1(a) for positive (1(a)) and negative (1(b)) $\beta_{4}$, as a function of the
wave number $k$. Cherenkov radiation appears when the phase velocity $V$
coincides with zero soliton velocity as shown in Fig.1(a). It is seen from
this figure that the Cherenkov radiation emitted from the soliton tail occurs
only when $\beta_{4}$ is positive. Therefore, below we consider only the case
of positive fourth-order dispersion coefficient $0<\beta_{4}\ll 1$ when the
Cherenkov radiation is present. For negative $\beta_{4}$ the soliton
interaction is only weakly affected by the small fourth-order dispersion term.

Linear stability of the dissipative soliton solution $u_{0}\left(\tau\right)$
is determined by calculating the eigenvalue spectrum $\lambda$ of the operator

$\hat{L}\left(\mathbf{u}_{0}\right)=\hat{L}_{0}+\hat{L}_{1}\left(\mathbf{u}_{0}\right),$
(2)

obtained by linearization of Eq. (1) around the soliton solution. Here
$\mathbf{u}_{0}=\left(\begin{array}[]{c}u_{0}\\\ u_{0}^{*}\end{array}\right)$
, $\hat{L}_{0}=\hat{L}(0)$ is the linear differential operator evaluated at
the stationary homogeneous solution $\mathbf{U}=\mathbf{U}_{0}$:

$\hat{L}_{0}=\left(\begin{array}[]{cc}-1-i\theta+2iI_{0}^{2}+\left(i+\delta\right)\partial^{2}_{\tau}+i\beta_{4}\partial^{4}_{\tau}&iU_{0}^{2}\\\
-iU_{0}^{*2}&-1+i\theta-2iI_{0}^{2}-\left(i-\delta\right)\partial^{4}_{\tau}-i\beta_{4}\partial^{4}_{\tau}\end{array}\right),$

and

$\hat{L}_{1}\left(\mathbf{u}_{0}\right)=\left(\begin{array}[]{cc}2iU_{0}^{*}u_{0}+2iU_{0}u_{0}^{*}+2i|u_{0}|^{2}&2iU_{0}u_{0}+iu_{0}^{2}\\\
-2iU_{0}^{*}u_{0}^{*}-iu_{0}^{*2}&-2iU_{0}^{*}u_{0}+2iU_{0}u_{0}^{*}-2i|u_{0}|^{2}\end{array}\right).$

We have calculated numerically the soliton solution and the eigenvalue
spectrum $\lambda$ of the operator $\hat{L}\left(\mathbf{u}_{0}\right)$ by
discretizing Eq. (1) on an uniform grid of $2000$ points on the interval
$\tau\in[0,80]$ with periodic boundary conditions. The result is shown in Fig.
2 for $\beta_{3}=\beta_{4}=\delta=0$. The continuous spectrum lies on the line
$\Re(\lambda)=-1$, while the discrete spectrum of the soliton is symmetric
with respect to this line (Barashenkov and Smirnov, 1996). For the parameter
values of Fig. 2 apart from two real eigenvalues: zero eigenvalue,
$\lambda=0$, associated with the translational symmetry of the Lugiato-Lefever
equation and symmetric one, $\lambda=-2$, soliton has two symmetric pairs of
complex conjugated eigenvalues. The right pair of these complex eigenvalues is
responsible for an Andronov-Hopf bifurcation taking place with the increase of
the injection parameter $S$.

Figure 1: Phase velocity $V$ of small dispersive waves with positive (a) and
negative (b) fourth-order dispersion coefficient $\beta_{4}$, and
$\beta_{3}=0$. Solid line corresponds to $\beta_{4}=0.025$ (a) and
$\beta_{4}=-0.025$ (b). Dashed line corresponds to $\beta_{4}=0$. The
parameter values are $S=1.8$, $\theta=3.5$, and $\delta=0.02$.

The decay rates of the soliton tails depend on the eigenvalues $\mu$
satisfying the characteristic equation

$\beta_{4}^{2}\mu^{8}+2\beta_{4}\mu^{6}+\left[1+\delta^{2}+2\beta_{4}(2I_{0}-\theta)\right]\mu^{4}+2(2I-\theta-\delta)\mu^{2}+\left[1-I_{0}^{2}+(2I_{0}-\theta)^{2}\right]=0.$
(3)

obtained by lineariazation of the Eq. (1) with $\partial_{t}U=0$ at the
homogeneous steady state solution $U=U_{0}$.

In the case when the high-order dispersion and spectral filtering are absent
$\beta_{3}=\beta_{4}=\delta=0$, Eq. (3) gives two pairs of complex conjugate
eigenvalues:

$\mu_{1,2}^{(0)}=\pm\sqrt{\theta-2I_{0}+i\sqrt{1-I_{0}^{2}}}$ (4)

and $\mu_{1,2}^{(0)*}$, which determine the decay and oscillation rates of the
soliton tails. For example, for $S=2.0$ and $\theta=3.5$ we have
$\mu_{1,2}^{(0)}=\pm\left(1.6837+0.275817i\right)$, which means that in the
absence of high-order dispersions the soliton tail oscillations are strongly
damped. This might explain the fact that without soliton Cherenkov radiation
it is hardly possible to observe soliton bound state formation experimentally
(Leo et al., 20103).

Figure 2: Soliton solution of the Lugiato-Lefever equation (1) with
$\beta_{3}=\beta_{4}=\delta=0$ (a) and eigenvalue spectrum obtained by
numerical linear stability analysis of this solution (b). Other parameters are
the same as in Fig. 1.

For nonzero but sufficiently small fourth-order dispersion coefficient,
$0<\beta_{4}\ll 1$, the eigenvalues (4) of Eq. (3) are only slightly
perturbed. However, in addition to (4) two more pairs of complex conjugate
eigenvalues, $\mu_{3,4}$ and $\mu_{3,4}^{*}$, appear. For zero spectral
filtering coefficient, $\delta=0$, they are given by

$\mu_{3,4}=\mp
i\sqrt{\frac{1}{2\beta_{4}}\left[1+\sqrt{1-4\beta_{4}\left(2I_{0}-\theta+i\sqrt{1-I_{0}^{2}}\right)}\right].}$
(5)

It is seen that real (imaginary) parts of $\mu_{3,4}$ in Eq. (5) vanish
(diverge) in the limit $\beta_{4}\to 0$. When the spectral filtering
coefficient is nonzero, $\delta>0$, analytical expressions for the eigenvalues
$\mu_{3,4}$ become very cumbersome. However, in the limit $\beta_{4}={\cal
O}(\delta)\ll 1$ we get:

$\displaystyle\mu_{3,4}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\beta_{4}}\left[\frac{\sqrt{\left(1+\delta/\beta_{4}\right)^{2}-I_{0}^{2}}}{2}\right.$
(6) $\displaystyle+$ $\displaystyle
i\left.\left(\frac{1}{\beta_{4}}+\frac{\theta-2I_{0}}{2}\right)+{\cal
O}\left(\delta\right)\right].$

Due to the presence of the eigenvalues $\mu_{3,4}$ and $\mu_{3,4}^{*}$ the
tails of the soliton of Eq. (1) with $\beta_{3}=0$ and $0<\beta_{4}\ll 1$
become weakly decaying and fast oscillating, which favors the formation of
soliton bound states, and can be referred to as the soliton Cherenkov
radiation (Akhmediev and Karlsson, 1995). Note, that when $\beta_{4}$ is
sufficiently small, the term $\delta/\beta_{4}$ describing in Eq. (6) the
contribution of spectral filtering into the real part of $\mu_{3,4}$ can lead
to a considerable increase of the decay rate of the soliton tails without
significant change of their oscillation frequency. For example, for $S=2.0$,
$\theta=3.5$, $\beta_{4}=0.025$ , and $\delta=0.02$ we get
$\mu_{3}=-0.123-6.529i$, while for the same parameter set and $\delta=0$ one
obtains $\mu_{3}=0.063-6.528i$. Numerically calculated intensity profile of
the soliton solution of Eq. (1) with small fourth-order dispersion coefficient
$\beta_{4}=0.025$ is depicted in Fig. 3 together with the corresponding
eigenvalue spectrum of the operator $\hat{L}\left(\mathbf{u}_{0}\right)$
defined by Eq. (2).

Note, that the proof of the reflectional symmetry property of the discrete
soliton spectrum with respect to the $\Re\lambda=-1$ line given in
(Barashenkov and Smirnov, 1996) is trivially generalized to the case when even
high order dispersions are present. Nevertheless, the soliton spectrum shown
in Fig. 3 does not possess this symmetry property due to the presence of
nonzero spectral filtering coefficient $\delta=0.02$. Furthermore, as it is
seen from Fig. 3, for $\delta=0.02$ real parts of the complex conjugate
eigenvalues, responsible for the Andronov-Hopf bifurcation of the soliton, are
shifted to the left from the imaginary axis as compared to those shown in Fig.
2 obtained for $\delta=0$.

Figure 3: Soliton solution of the Lugiato-Lefever equation (1) with
$\beta_{4}=0.025$ and $\delta=0.02$ (a); eigenvalue spectrum obtained by
numerical linear stability analysis of this solution (b). Other parameters are
the same as in Fig. 1.

Sufficiently far away from the soliton core its trailing tail can be
represented in the form

$u_{0}(\tau)\sim
a_{3}e^{\mu_{3}\tau}+a_{4}e^{\mu_{3}^{*}\tau},{\text{when}}\quad\tau\to+\infty,$
(7)

where the Cherenkov radiation amplitude $a_{3}$ is exponentially small in the
limit $\beta_{4}\to 0$ (Akhmediev and Karlsson, 1995; Karpman, 1993),
$a_{4}=p_{a}a_{3}^{*}$, and for $\beta_{4}={\cal O}(\delta)\ll 1$ we get

$p_{a}=i\frac{1-\sqrt{1-I_{0}^{2}}}{A_{0}^{*2}\left(\frac{\delta}{\beta_{4}\sqrt{1-I_{0}^{2}}}+1\right)}+{\cal
O}\left(\delta\right),$ (8)

where $p_{a}$ is independent of $\beta_{4}$ at $\delta=0$. Numerically for
$S=2.0$, $\theta=3.5$, $\delta=0.02$, and $\beta_{4}=0.025$ we obtain
$p_{a}\approx 0.0571+0.0833i$.

##  III Interaction between dissipative solitons

. Two or more solitons will interact through their overlapping oscillatory
tails when they are sufficiently close to one another. In what follows, we
investigate the interaction between two dissipative solitons. We consider the
limit of weak overlap when the solitons are well separated from each other and
derive the interaction equations describing the slow time evolution of the
soliton coordinates denoted by $\tau_{1,2}$. To this end, let us first rewrite
Eq. (1) in a general form:

$\partial_{t}\mathbf{U}=\hat{F}\mathbf{U},$ (9)

where $\mathbf{U}=\left(\begin{array}[]{c}U\\\ U^{*}\end{array}\right)$,
$\hat{F}\mathbf{U}=\left(\begin{array}[]{c}\hat{f}U\\\
\hat{f}^{*}U^{*}\end{array}\right)$, and $\hat{f}$ is the differential
operator defined by the RHS of Eq. (1). We look for the solution of Eq. (9) in
the form

$U(\tau,t)=U_{0}+u_{1}+u_{2}+\Delta u(\tau,t).$ (10)

Here $u_{1,2}=u_{0}\left(\tau-\tau_{1,2}\right)$ are two unperturbed soliton
solutions, with slowly evolving in time coordinates
$\tau_{1,2}\left(\varepsilon t\right)$, $\Delta u(\tau,t)={\cal
O}(\varepsilon)$ is a small correction to the superposition of two solitons,
and small parameter $\varepsilon$ describes the weakness of the overlap of the
two solitons. Substituting Eq. (10) into the model equation (9) and collecting
first order terms in $\varepsilon$ we obtain the following linear
inhomogeneous equation for
$\mathbf{\mathbf{\Delta}u}=\left(\begin{array}[]{c}\Delta u\\\ \Delta
u^{*}\end{array}\right)$:

$\hat{L}(\mathbf{u}_{1}+\mathbf{u}_{2})\mathbf{\Delta
u}=-\partial_{x}\mathbf{u}_{1}\partial_{t}\tau_{1}-\mathbf{\partial_{\mathbf{\mathrm{\mathit{x}}}}u_{\mathsf{\mathrm{2}}}}\partial_{t}\tau_{2}-\hat{F}(\mathbf{u}_{1}+\mathbf{u}_{2}),$
(11)

where the linear operator $\hat{L}\left(\mathbf{u}\right)$ is defined by Eq.
(2). Due to the transnational invariance of Eq. (1) this linear operator
evaluated at the soliton solution $\mathbf{u}_{0}$ has zero eigenvalue
corresponding to the so-called transnational neutral (or Goldstone) mode
$\mathbf{v}_{0}=\left(\begin{array}[]{c}v_{0}\\\ v_{0}^{*}\end{array}\right)$
with $v_{0}=du_{0}/d\tau$,
$\hat{L}\left(\mathbf{u}_{0}\right)\mathbf{v}_{0}=0$. The adjoint linear
operator $\hat{L}^{\dagger}\left(\mathbf{u}\right)$ obtained from
$\hat{L}\left(\mathbf{u}\right)$ by transposition and complex conjugation also
has zero eigenvalue with the eigenfunction
$\mathbf{w}_{0}=\left(\begin{array}[]{c}w_{0}\\\ w_{0}^{*}\end{array}\right)$,
which is referred below as the “adjoint neutral mode”,
$\hat{L}^{\dagger}\left(\mathbf{u}_{0}\right)\mathbf{w}_{0}=0$. Below we will
assume that $\mathbf{w}_{0}$ satisfies the normalization condition
$\left\langle\mathbf{w}_{0}\cdot\mathbf{u}_{0}\right\rangle=\int_{-\infty}^{\infty}\left(\mathbf{w}_{0}\cdot\mathbf{u}_{0}\right)d\tau=2\int_{-\infty}^{\infty}\Re\left(w_{0}^{*}u_{0}\right)d\tau=1$.

Figure 4: RHS of Eq. (16) as a function of the soliton separation
$\tau_{2}-\tau_{1}$. Black (red) dots indicate the separations of the two
solitons in stable (unstable) bound states calculated numerically. Parameter
values: $S=2.0$, $\theta=3.5$, $\delta=0.02$ and $\beta_{4}=0.025$.

Far away from the soliton core the asymptotic behavior of adjoint neutral mode
is given by:

$w_{0}(\tau)\sim
b_{3}e^{\mu_{{}^{3}}^{*}\tau}+b_{4}e^{\mu_{3}\tau},\quad\tau\to+\infty,$ (12)

with $b_{4}=p_{b}b_{3}^{*}$, where asymptotic expression for $p_{b}$ coincides
with that of $p_{a}$ given by Eq. (8).

Figure 5: Stable bound state of two dissipative solitons, $\delta=0.02$. (a) –
intensity distribution, (b) – eigenvalue spectrum. Other parameter values are
the same as for Fig. 4.

Figure 6: Unstable bound state of two dissipative solitons, $\delta=0.02$. (a)
– intensity distribution, (b) – eigenvalue spectrum. Other parameter values
are the same as for Fig. 4.

.

Figure 7: The same bound bound state as shown in Fig. 5, but calculated for
$\delta=0$. (a) – intensity distribution, (b) – eigenvalue spectrum. Bound
state is unstable with respect to an Andronov-Hopf bifurcation.

When the two interacting solitons are located sufficiently far away from one
another the solvability conditions of Eq. (11) can be written as

$\partial_{t}\tau_{1,2}\approx G_{1,2},\quad
G_{1,2}=\left\langle\mathbf{w}_{1,2}\cdot\hat{F}(\mathbf{u}_{1}+\mathbf{u}_{2})\right\rangle,$
(13)

where we approximated the adjoint neutral modes of the operator
$\hat{L}^{\dagger}(\mathbf{u}_{1}+\mathbf{u}_{2})$ by the adjoint neutral
modes $\mathbf{w}_{1,2}=\mathbf{w}_{0}(\tau-\tau_{1,2})$ of the operators
$\hat{L^{\dagger}}(\mathbf{u}_{1,2})$.

In order to derive the soliton interaction equations we need to calculate
$G_{1,2}$ in Eq. (13). To this end we split the integral in Eq. (13) into two
parts and using the relations
$\hat{L}^{\dagger}(\mathbf{u}_{1,2})\mathbf{w}_{1,2}=0$ together with the fact
that $\mathbf{u}_{1}$ and $\mathbf{w}_{1}$ ($\mathbf{u}_{2}$ and
$\mathbf{w}_{2}$) are small for $\tau\in[0,+\infty)$ ($\tau\in(-\infty,0]$),
where the origin of coordinates $\tau=0$ corresponds to the central point
between two solitons, $\left(\tau_{2}+\tau_{1}\right)/2=0$, we get

$\displaystyle G_{1,2}$ $\displaystyle=$
$\displaystyle\left\langle\mathbf{w}_{1,2}\cdot\hat{F}(\mathbf{u}_{1}+\mathbf{u}_{2})\right\rangle_{1,2}+\left\langle\mathbf{w}_{1,2}\cdot\hat{F}(\mathbf{u}_{1}+\mathbf{u}_{2})\right\rangle_{2,1}$
(14) $\displaystyle\approx$
$\displaystyle\left\langle\mathbf{w}_{1,2}\cdot\hat{F}(\mathbf{u}_{1}+\mathbf{u}_{2})\right\rangle_{1,2}$
$\displaystyle\approx$
$\displaystyle\left\langle\mathbf{w}_{1,2}\cdot\hat{L}(\mathbf{u}_{1,2})\mathbf{u}_{2,1}\right\rangle_{1,2}-\left\langle\hat{L}^{\dagger}(\mathbf{u}_{1,2})\mathbf{w}_{1,2}\cdot\mathbf{u}_{2,1}\right\rangle_{1,2}$
$\displaystyle=$ $\displaystyle\left(\delta+i\right)\left[\left\langle
w_{1,2}\partial^{2}_{\tau}u_{2,1}\right\rangle_{1,2}-\left\langle
u_{2,1}\partial^{4}_{\tau}w_{1,2}\right\rangle_{1,2}\right]$ $\displaystyle+$
$\displaystyle i\beta_{4}\left[\left\langle
w_{1,2}\partial^{4}_{\tau}u_{2,1}\right\rangle_{1,2}-\left\langle
u_{2,1}\partial^{4}_{\tau}w_{1,2}\right\rangle_{1,2}\right]+c.c.$

with
$\left\langle\mathbf{w}\cdot\mathbf{u}\right\rangle_{1}=\int_{-\infty}^{0}\left(\mathbf{w}\cdot\mathbf{u}\right)d\tau$,
$\left\langle\mathbf{w}\cdot\mathbf{u}\right\rangle_{2}=\int_{0}^{\infty}\left(\mathbf{w}\cdot\mathbf{u}\right)d\tau$
and $\hat{L}^{\dagger}(\mathbf{u}_{1,2})\mathbf{w}_{1,2}=0$.

Next, performing integration by parts and using the symmetry properties of the
soliton and its neutral modes, $u_{0}(\tau)=u_{0}(-\tau)$,
$\partial_{\tau}u_{0}(\tau)=-\partial_{\tau}u_{0}(-\tau)$,
$w_{0}(\tau)=-w_{0}(-\tau)$, and
$\partial_{\tau}w_{0}(\tau)=\partial_{\tau}w_{0}(\tau)$ we get:

$\displaystyle G_{1,2}$
$\displaystyle\approx\pm\left[\left(\delta+i\right)\left(w_{1,2}^{*}\partial_{\tau}u_{2,1}-u_{2,1}\partial_{\tau}w_{1,2}^{*}\right)\right.$
(15)
$\displaystyle+i\beta_{4}\left(w_{1,2}^{*}\partial^{3}_{\tau}u_{2,1}-u_{2,1}\partial^{3}_{\tau}w_{1,2}^{*}-\partial_{\tau}w_{1,2}^{*}\partial^{2}_{\tau}u_{2,1}\right.$
$\displaystyle+\left.\left.\partial^{2}_{\tau}w_{1,2}^{*}\partial_{\tau}u_{2,1}\right)\right]_{\tau=0}+c.c.=$
$\displaystyle\pm\left[\left(\delta+i\right)\partial_{\tau}(w_{0}^{*}u_{0})-i\beta_{4}\left(w_{0}^{*}\partial^{3}_{\tau}u_{0}+u_{0}\partial^{3}_{\tau}w_{0}^{*}\right.\right.$
$\displaystyle+\left.\left.\partial_{\tau}\left(\partial_{\tau}w_{0}^{*}\partial_{\tau}u_{0}\right)\right)\right]_{\tau=(\tau_{2}-\tau_{1})/2}+c.c.$

Finally, substituting into Eq. (15) the asymptotic relations (7) and (12) we
obtain:

$\frac{d\left(\tau_{2}-\tau_{1}\right)}{dt}\approx-\frac{12}{\sqrt{\beta_{4}}}e^{-\gamma(\tau_{2}-\tau_{1})}\Re\left[\left(1-i\frac{\delta}{3}\right)\left(a_{3}b_{3}^{*}e^{-i\Omega(\tau_{2}-\tau_{1})}-p_{a}p_{b}^{*}a_{3}^{*}b_{3}e^{i\Omega(\tau_{2}-\tau_{1})}\right)\right],$
(16) $\frac{d\left(\tau_{2}+\tau_{1}\right)}{dt}=0,$ (17)

where
$\gamma=\Re(\mu_{3})\approx\left(\sqrt{\beta_{4}}/2\right)\sqrt{\left(1+\delta/\beta_{4}\right)-I_{0}^{2}}$
, $\Omega=-\Im(\mu_{3})\approx
1/\sqrt{\beta_{4}}+\sqrt{\beta_{4}}\left(\theta-2I_{0}\right)$, and the
Cherenkov radiation coefficients $a_{3}$ and $b_{3}$ are exponentially small
in the limit $\beta_{4}\to 0$. For $S=2.0$, $\theta=3.5$, $d=0.02$, and
$\beta_{4}=0.025$ numerically we get $a_{3}\approx-0.158+0.149i$ and
$b_{3}\approx 0.017+0.136i$. Finally, neglecting ${\cal O}\left(\delta\right)$
terms and taking into account that in the leading order in $\delta$ we have
$p_{a}=p_{b}\equiv p$, Eq. (16) can be rewritten in the form:

$\frac{d\left(\tau_{2}-\tau_{1}\right)}{dt}\approx\frac{12}{\sqrt{\beta_{4}}}e^{-\gamma(\tau_{2}-\tau_{1})}|a_{3}b_{3}|\left(|p|^{2}-1\right)\cos\left[\Omega\left(\tau_{2}-\tau_{1}\right)+\arg\left(b_{3}/a_{3}\right)\right].$
(18)

The RHS of Eq. (18) is plotted in Fig. 4, where the intersections of the black
solid line with axis of abscissas correspond to the soliton bound states.
Examples of stable and unstable soliton bound states calculated numerically
are shown in Figs. 5 and 6, respectively, together with the most unstable
eigenvalues of the operator $\hat{L}$ evaluated on the bound state solutions.

Finally, in Fig. 7 we present the same soliton bound state as the one shown in
Fig. 5, but calculated for $\delta=0$. It is seen that the eigenvalue spectrum
of this state contains many discrete eigenvalues, which split from the
continuous spectrum, and that it is oscillatory unstable due to the presence
of two complex conjugate eigenvalues with positive real parts. Therefore, we
can conclude that in the absence of spectral filtering the one-dimensional
asymptotic equations (16)-(18) can be insufficient to describe the soliton
interaction. The derivation of the interaction equations taking into account
an Andronov-Hopf bifurcation of the soliton bound states in the presence of
fourth-order dispersion is beyond the scope of this study. A related problem
concerning the effect of oscillatory instability on the soliton interaction
was studied in (Turaev et al., 2012).

##  IV Conclusions

. We have considered an all fiber photonic crystal cavity coherently driven by
an injected field. The intracavity field inside the fiber experiences self-
phase modulation, dispersion, optical injection, and optical losses. Its
space-time evolution can be described by the Lugiato-Lefever equation with
high order dispersion, where, in addition, we have taken into account small
spectral filtering term. We have first discussed the properties of a single
dissipative soliton and derived asymptotic expressions for the soliton
Cherenkov radiation amplitudes. We have focused our analysis on the regime,
where the forth order dispersion and the spectral filtering coefficients are
small, $0<\beta_{4},\delta\ll 1$. Second, we have investigated the interaction
between two dissipative solitons in the case when they are well separated from
each other. Assuming a weak overlap of soliton tails, we have established
analytically the interaction law (Eqs. 18) governing the slow time evolution
of the coordinates of two interacting solitons. We have shown that although
the Cherenkov radiation due to the small fourth-order dispersion can strongly
enhance the soliton interaction and thus lead to the formation of a large
number of soliton bound states, in the absence of spectral filtering these
states can be unstable with respect to an oscillatory instability even when
single soliton is well below the Andronov-Hopf bifurcation threshold. This
means that one-dimensional equation (18) can be insufficient to describe the
interaction of solitons in the generalized Lugiato-Lefever model (1) with zero
spectral filtering coefficient, $\delta=0$. On the other hand, the inclusion
of small but sufficiently large spectral filtering, $0<\delta\ll 1$, allows to
stabilize oscillatory unstable bound states and validate the one-dimensional
interaction equation (18).

###### Acknowledgements.

We also acknowledge the support from the Deutsche Forschungsgemeinschaft (DFG-
RSF project No. 445430311), French National Research Agency (LABEX CEMPI,
Grant No. ANR-11- LABX-0007) as well as the French Ministry of Higher
Education and Research, Hauts de France council and European Regional
Development Fund (ERDF) through the Contrat de Projets Etat-Region (CPER
Photonics forSociety P4S). M. Tlidi is a Research Director at the Fonds
National de la Recherche Scientifique (Belgium). A. Vladimirov and M. Taki
acknowledge the support from Invited Research Speaker Programme of the Lille
University.

## References

  * Fortier and Baumann (2019) T. Fortier and E. Baumann, Communications Physics 2, 1 (2019).
  * Kippenberg et al. (2018) T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, Science 361 (2018).
  * Herr et al. (2014) T. Herr, V. Brasch, J. D. Jost, C. Y. Wang, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, Nature Photonics 8, 145 (2014).
  * Leo et al. (20103) F. Leo, S. Coen, P. Kockaert, S.-P. Gorza, P. Emplit, and M. Haelterman, Nature Photonics 4, 471 (20103).
  * Lugiato and Lefever (1987) L. A. Lugiato and R. Lefever, Phys. Rev. Lett. 58, 2209 (1987).
  * Chembo et al. (2017) Y. K. Chembo, D. Gomila, M. Tlidi, and C. R. Menyuk, The European Physical Journal D 71, 198 (2017).
  * Lugiato et al. (2018) L. Lugiato, F. Prati, M. Gorodetsky, and T. Kippenberg, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20180113 (2018).
  * Haelterman et al. (1992) M. Haelterman, S. Trillo, and S. Wabnitz, Optics communications 91, 401 (1992).
  * Maleki et al. (2010) L. Maleki, V. Ilchenko, A. Savchenkov, W. Liang, D. Seidel, and A. Matsko, in _2010 IEEE International Frequency Control Symposium_ (IEEE, 2010), pp. 558–563.
  * Matsko et al. (2011) A. Matsko, A. Savchenkov, W. Liang, V. Ilchenko, D. Seidel, and L. Maleki, Optics letters 36, 2845 (2011).
  * Chembo and Menyuk (2013) Y. K. Chembo and C. R. Menyuk, Phys. Rev. A 87, 053852 (2013).
  * Tlidi et al. (2007) M. Tlidi, A. Mussot, E. Louvergneaux, G. Kozyreff, A. G. Vladimirov, and M. Taki, Optics letters 32, 662 (2007).
  * Tlidi and Gelens (2010) M. Tlidi and L. Gelens, Optics letters 35, 306 (2010).
  * Akhmediev et al. (1990) N. Akhmediev, V. Korneev, and N. Mitskevich, Radiophysics and quantum electronics 33, 95 (1990).
  * Tlidi et al. (2013) M. Tlidi, L. Bahloul, L. Cherbi, A. Hariz, and S. Coulibaly, Physical Review A 88, 035802 (2013).
  * Milián and Skryabin (2014) C. Milián and D. Skryabin, Opt. Express 22, 3732 (2014).
  * Vladimirov et al. (2018) A. G. Vladimirov, S. V. Gurevich, and M. Tlidi, Physical Review A 97, 013816 (2018).
  * Soto-Crespo et al. (2003) J. M. Soto-Crespo, N. Akhmediev, P. Grelu, and F. Belhache, Optics Letters 28, 1757 (2003).
  * Turaev et al. (2012) D. Turaev, A. G. Vladimirov, and S. Zelik, Phys. Rev. Lett. 108, 263906 (2012).
  * Olivier et al. (2006) M. Olivier, V. Roy, and M. Piché, Optics Letters 31, 580 (2006).
  * Akhmediev and Karlsson (1995) N. Akhmediev and M. Karlsson, Phys. Rev. A 51, 2602 (1995).
  * Parra-Rivas et al. (2016) P. Parra-Rivas, E. Knobloch, D. Gomila, and L. Gelens, Physical Review A 93, 063839 (2016).
  * Parra-Rivas et al. (2017) P. Parra-Rivas, D. Gomila, P. Colet, and L. Gelens, The European Physical Journal D 71, 198 (2017).
  * Gorshkov and Ostrovsky (1981) K. Gorshkov and L. Ostrovsky, Physica D: Nonlinear Phenomena 3, 428 (1981).
  * Karpman and Solov’ev (1981) V. Karpman and V. Solov’ev, Physica D: Nonlinear Phenomena 3, 487 (1981).
  * Scroggie et al. (1994) A. Scroggie, W. Firth, G. McDonald, M. Tlidi, R. Lefever, and L. A. Lugiato, Chaos, Solitons & Fractals 4, 1323 (1994).
  * Barashenkov and Smirnov (1996) I. Barashenkov and Y. S. Smirnov, Physical Review E 54, 5707 (1996).
  * Karpman (1993) V. I. Karpman, Phys. Rev. E 47, 2073 (1993).