Train Localization During GNSS Outages: A
Minimalist Approach Using Track Geometry And
IMU Sensor Data

This work has been funded by Trafikverket (The Swedish Transport Administration),  through Excellence Area 6, TRV 2024/2731.

Wendi Löffler
Information Science and Engineering

KTH Royal Institute of Technology

Stockholm, Sweden

Mats Bengtsson
Information Science and Engineering

KTH Royal Institute of Technology

Stockholm, Sweden

Train localization during GNSS outages presents challenges for ensuring fail-safe and accurate positioning in railway networks. This paper proposes a minimalist approach exploiting track geometry and IMU sensor data. By integrating a discrete track map as a Look-Up Table (LUT) into a Particle Filter (PF) based solution, accurate train positioning is achieved with only an IMU sensor and track map data. The approach is tested on an open railway positioning data set, showing that accurate positioning (absolute errors below 10 m) can be maintained during GNSS outages up to 30 s in the given data. We simulate outages on different track segments and show that accurate positioning is reached during track curves and curvy railway lines. The approach can be used as a redundant complement to established positioning solutions to increase the position estimate's reliability and robustness.

train positioning, discrete track map, particle filter, statistical sensor fusion

ETCSEuropean Train Control System
GNSSGlobal Navigation Satellite Systems
IMUInertial Measurement Unit
SLAMSimultaneous localization and mapping
PFParticle Filter
EKFExtended Kalman Filter
LUTLook-Up Table
ENUEast North Up
DTMDiscrete Track Map
ZVUZero Velocity Update
KFKalman Filter
CDFCumulative Density Function
EKFMMExtended Kalman Filter with Map Matching
CTRAConstant Turn Rate and Acceleration

§ INTRODUCTION

Determination of train positions within the railway network must be fail-safe and of high accuracy. Inaccurate train position data can lead to delays in operation [1]. In systems operating today, train position is determined using trackside equipment such as axle counters and balises. These solutions are vulnerable to vandalism and weather influence; hence, to reduce maintenance costs on the trackside, train-borne solutions that are not dependent on trackside equipment have been a topic of interest in research for some years.

A widely used approach in tracking problems is fusion of Global Navigation Satellite Systems (GNSS) and Inertial Measurement Units (IMU) [2]. The challenges consist of adjusting the method to the high safety standards according to EN50126 [3] and the specifics of the railway environment.
Several KF based solutions have been proposed using GNSS, IMU, and odometer data [4, 5].

Considered set up in this paper. Our proposed solution exploits an IMU and a track map.

Unlike other transport modes such as aviation and automotive, train motion is restricted to the tracks; a constraint that can be used to improve the accuracy of the position estimate and different methods have been proposed fusing GNSS/IMU, and a track map [6, 7, 8, 9]. Geometric track maps can be used for train positioning in different ways: in [10, 11, 12], the authors use linear interpolation between discrete coordinates, whereas in [13], spline interpolation is used. On the other hand, [14] propose a formulation of the track map based on the geometrical elements in track construction, promising a more precise representation. However, [15] shows that linear interpolation between discrete coordinates with feasible distance offers sufficient accuracy in most applications.

Given the high safety requirements, GNSS cannot provide the required availability and accuracy in the complex propagation environments, which are encountered in railway settings [16, 17]. Tunnels, urban areas, heavily wooded districts, or deep valleys cause GNSS outages. Polar regions, as relevant in Sweden, cause a decrease in GNSS accuracy. In  [18] an extensive study is performed using positioning data from 80.000 km of rails and observe GNSS outages 18 % of the time, with the majority being below ten minutes. To overcome this problem, recent literature has considered GNSS-redundant solutions. [19, 20] find that railway tracks possess a unique ferromagnetic track signature that can be used for localization. However, changed track conditions or passing trains can disturb or change this signature. Other approaches have considered different kinds of signatures such as altitude [19]. In [21], the authors developed a map-aided train navigation based on IMU measurements for subway systems.

To meet the high safety requirements for train positioning solutions, there is a need for redundant positioning methods using complementing sensors. Therefore, in this paper, we develop a train-borne PF based positioning solution that can complement GNSS, odometer, or magnetometer-reliant methods. It can also be used as an extension to improve filter performance in areas where GNSS is sparsely available. <Ref> visualizes the examined situation. The key strength of the algorithm is minimalist sensor usage and an innovative formulation and usage of a DTM. Only measurements of tangential acceleration, centripetal acceleration, and yaw angle turn rate - all recorded by one IMU - and the railway track coordinates are required. We provide a flexible and adaptable formulation of the DTM, which is incorporated in a PF based positioning solution. The positioning problem is reduced to one dimension, leading to a two-dimensional state vector. Therefore, a low number of particles can be used, and computational costs are lower than in other PF based solutions. The method is tested on the open data set [22], which consists of a secondary line in a challenging environment.

§ ENHANCED GEOMETRIC TRACK MAP AND TRACK CONSTRAINT FORMULATION

The following section describes the design and usage of the *DTM used in the estimation. The dynamic motion model explained in <Ref> computes the train's position in 1D, the dimension being the traveled distance along the track. To acquire the corresponding map features and to map the position into 2D space, we require a geometric track map, which, in our proposed method, serves as a LUT. Consider a set of $M$ discrete points along the track. The $i$th point along the track consists of an identifier $d_i$ and six features $\mu_i$ recorded at this point. The track map $ \mathcal{M} $ is expressed as follows:

\begin{align}
\label{eqn:dtm}
\mathcal{M} &= \left\{ \left(
d_i, \boldsymbol{\mu}_i  \right) \right\} \\
\boldsymbol{\mu}_i &= \begin{bmatrix}
p_{x,i} & p_{y,i} & p_{z,i} & \kappa_i & \theta_{x,i} & \theta_{y,i} & \theta_{z,i}
\end{bmatrix}
\end{align}

$d_i$ is an identifier of a recorded point along the track corresponding to the distance along the track. Distance $d_i = 0$ is set arbitrarily, e.g., at the starting station or the beginning of the map. Note that $d_i-d_{i-1}$ must be the distance along the track, not the Euclidean distance, to avoid underestimating the traveled distance.
The map features included are the following:

* $p_{x,i}$, $p_{y,i}$ and $p_{z,i}$ correspond to the 3D geographical position of the point in a fixed reference frame.

* $\kappa_i$ is the computed curvature of the track, see <ref>.

* $\theta_{x,i}$, $\theta_{y,i}$ and $\theta_{z,i}$ are the track orientation expressed in roll, pitch and yaw angle, respectively.

<Ref> gives a visual overview of the variables stored in the track map.

Visualization of variables stored from the track map

The orientation information is used for IMU pre-processing in the filter, curvature is the main feature for tracking, and position is the filter output. Further features for filter extensions can be included in $\mathcal{M}$, such as details about stations, turnouts, or further track signatures.

§.§ Derivation of map curvature

To calculate the track curvature, a twice differentiable parametric representation of the 2D plane curve is required:

\begin{equation}
\gamma(d) = (x(d),y(d))
\end{equation}

In the DTM as presented in (<ref>), an arc-length parameterization on distance $d$ is given. A cubic spline interpolation between recorded map points ensures differentiability. The curvature is then estimated as

\begin{equation}
\kappa(d) = \frac{x' y'' - y' x''}{(x'^2+y'^2)^{3/2}}
\end{equation}

Spline interpolation suffers from some drawbacks, e.g., closely spaced data points can lead to oscillations in the interpolated curve. In [15], the author found the most accurate computation of railway track curvature by reducing the data points used for spline interpolation using the Ramer-Douglas-Peucker algorithm and interpolating to the original density after spline interpolation. Importantly, this strategy does not alter the 2D coordinates stored in the LUT.

§.§ Map Matching and interpolation of map data

To derive suitable parameters from the DTM, we find the largest $d_i^{(\mathcal{M})}$ that is smaller than the estimated distance $d_k$ at time step $k$:

\begin{align}
\label{eqn:mapping1}
q_1^{(\mathcal{M})} &= \mathrm{max}\{ d_i^{(\mathcal{M})} | d_i^{(\mathcal{M})} \in \mathcal{M} \, \mathrm{and} \, d_i^{(\mathcal{M})}<d_k \} \\
q_2^{(\mathcal{M})} &= q_1^{(\mathcal{M})}+1
\end{align}

The parameter set $\mu_k$ retrieved from the DTM is obtained by linear interpolation between $q_1^{(\mathcal{M})}$ and $q_2^{(\mathcal{M})}$:

\begin{align}
r_k &= \frac{d_k - q_1^{(\mathcal{M})}}{|q_2^{(\mathcal{M})}-q_1^{(\mathcal{M})}|} \\
\label{eqn:mapping2}
\mu_k &= (1-r_k)\cdot \mu(q_1^{(\mathcal{M})})+r_k \cdot \mu (q_2^{(\mathcal{M})})
\end{align}

For the simplicity of the formulation, we assumed a forward travel direction.
Note that

\begin{equation}
\left| \left| \begin{bmatrix}
p_x & p_y
\end{bmatrix}_{i+1}^{(\mathcal{M})} - \begin{bmatrix}
p_x & p_y
\end{bmatrix}_i^{(\mathcal{M})} \right| \right| \leq \left|d_{i+1}^{(\mathcal{M})} - d_i^{(\mathcal{M})}\right|
\end{equation}

We choose a linear interpolation between two recorded map points to compute $\mu_k$, which in a setup with 2D states would lead to underestimation of the traveled distance and velocity. Because of the format of the LUT, this problem is eliminated in our approach. Further, other interpolation methods can easily be implemented if required.

§ MAP-CONSTRAINED PARTICLE FILTER

The following section describes the algorithm developed in this paper. An overview is given in <Ref>.

Overview of complete algorithm

We use a PF based solution with a 2D linear state and a highly non-linear measurement update. The train's motion is modeled using a 1D constant acceleration model included in the particle diffusion step.

§.§ Dynamic Motion Model

The following section corresponds to the Prediction block in <Ref>.
Due to the constrained motion of the train along the tracks, we adopt a 1D motion model, as detailed below. Train motion consists of slow changes in acceleration. Hence, the constant acceleration model, as found in [23] is a suitable model of the train dynamics. The state $\mathbf{x}$ at iteration $k$ consists of two variables:

\begin{equation}
\mathbf{x}_k = \begin{bmatrix}
d_k & v_k
\end{bmatrix}^T
\end{equation}

with $v_k$ being velocity and $d_k$ traveled distance.
The constant acceleration model describes the motion dynamics as

\begin{align}
\label{eqn:motionmodel}
\begin{bmatrix}
d \\ v
\end{bmatrix}_{k+1} = \begin{bmatrix}
1 &T \\ 0& 1
\end{bmatrix} \cdot \begin{bmatrix}
d \\ v
\end{bmatrix}_k &+ \begin{bmatrix}
T^2/2 \\ T
\end{bmatrix} \cdot u_k + \nu_k, \\
\nu_k &\sim \mathcal{N}(0,\sigma_u^2) \notag
\end{align}

where $T$ is the sampling time. The system input $u$ is acceleration along the direction of movement measured by the IMU $a_x^{\mathrm{(IMU)}}$. The measurement is corrected using the estimated bias, as detailed in <Ref>.

§.§ Observation Model

This section corresponds to the Measurement Update block in <Ref>.
Consider the following measurements:

\begin{equation}
\mathbf{z}_k = [
\underbrace{ \begin{matrix} p_x^{(\text{GNSS})} & p_y^{(\text{GNSS})} & v^{(\text{GNSS})} \end{matrix} }_{\text{if available}} \quad \begin{matrix} \omega_z^{(\text{IMU})} &  a_y^{(\text{IMU})}
\end{matrix}]_k
\end{equation}

where $a_y^{(\text{IMU})}$ and $\omega_z^{(\text{IMU})}$ are measurements of the lateral acceleration and angular rate of the yaw angle. Both IMU measurements were corrected using the estimated bias, as explained in <Ref>.
GNSS measurements are generally not available except for a few irregular time intervals, which mostly reduces the measurement vector to

\begin{equation}
\mathbf{z}_k = \begin{bmatrix}
\omega_z^{(\text{IMU})} &  a_y^{(\text{IMU})}
\end{bmatrix}_k
\end{equation}

The predicted measurement $\hat{\mathbf{z}}_k$ is expressed as a function of the state $\hat{\mathbf{x}}_k$ and the DTM $ \mathcal{M} $:

\begin{equation}
\label{eqn:observationmodel}
\hat{\mathbf{z}}_k = h_k(\hat{\mathbf{x}}_k, \mathcal{M}(\hat{d}_k))
\end{equation}

The estimated distance $\hat{d}_k$ has two purposes: it is included in the state variable $\hat{\mathbf{x}}_k$, but also the input to the LUT $\mathcal{M}$.
Predicted measurements are calculated in three different ways:

* velocity $v$ is directly captured in the state

* position $p_x$, $p_y$ is a feature stored in the DTM

* angular rate $\omega_z$ and lateral acceleration $a_y$ are calculated from a combination of state variables and map features.

Visualization of observation model

$h_k$ is specified in the following way:

\begin{align}
\label{eqn:virtualmeas}
h_k ( \mathbf{x}_k,\mathcal{M} ) &= \begin{bmatrix}
d_k \rightarrow \mathcal{M} \rightarrow p_{x,k} \\
d_k \rightarrow \mathcal{M} \rightarrow p_{y,k} \\
v_k \\
v_k \cdot \kappa_k(d_k) \\
v_k^2 \cdot \kappa_k(d_k)
\end{bmatrix} + \mathbf{\eta}_k \\
\kappa_k &= d_k \rightarrow \mathcal{M} \rightarrow \kappa_k
\end{align}

Here, the right arrows indicate how $\mathcal{M}$ is used as a LUT.
The relation in (<ref>) is visualized in <Ref>.

Note that it is not necessary to express the relation between the variables stored in $\mathcal{M}$ and distance $d_k$ in closed form.

§.§ IMU Processing

This section describes the IMU processing, ZVU and Bias Correction block shown in <Ref>.

In the IMU processing block, the distortion caused by gravitational acceleration is removed. This is achieved by transforming the 3D accelerometer data from the sensor frame to the inertial frame and removing the gravitational acceleration from $a_z^{\text{inertial}}$. Sensor frame and inertial frame are related via the orientation matrix $\mathbf{R}(\boldsymbol{\theta})$:

\begin{equation}
\mathbf{a}^{\text{inertial}} = \mathbf{R}(\mathbf{\boldsymbol{\theta}}^{(\mathcal{M})}) \cdot \mathbf{a}^{\text{sensor}}
\end{equation}

where $\boldsymbol{\theta}$ is the orientation. In classic IMU applications, orientation is estimated from the gyroscope data, which is sensitive to integration errors and bias itself. Therefore, we choose to compute $\boldsymbol{\theta}$ from the estimated distance of the previous filter update $\hat{d}_{k-1}$ using the DTM, as described in <Ref>. The correction step is as follows:

\begin{equation}
\label{eqn:orientation}
\mathbf{a}^{\text{sensor}} =\Tilde{\mathbf{a}}^{\text{sensor}} - \mathbf{R}^{-1}(\boldsymbol{\theta}) \begin{bmatrix}
0 & 0 & g
\end{bmatrix}^T
\end{equation}

where $\Tilde{\mathbf{a}}^{\text{sensor}}$ is the raw IMU measurement. $\mathbf{a}^{\text{sensor}}$ is noted as $\mathbf{a}$ in the remaining paper.

Both accelerometer and gyroscope measurements contain several deterministic and stochastic error sources. The sensors are assumed to be calibrated in terms of scale factors, misalignment, and cross-coupling [24]. The scaling error is assumed to be small in comparison to bias; therefore, the simplified error model is assumed as follows:

\begin{equation}
\Tilde{\lambda} = \lambda + B + \delta B + \eta
\end{equation}

where $\lambda$ is the true value, $B$ is bias, $\delta B$ bias instability and $\eta$ the sensor noise. $\Tilde{\lambda}$ is the value measured by the accelerometer or gyroscope on one axis.
To estimate $B$ and account for $\delta B$, $\Tilde{\lambda}$ is measured during a stand-still in the filter initialization. At each further stand-still, the bias estimate is updated using a moving average:

\begin{equation}
\label{eqn:biasestimation}
B_k =
\begin{cases}
\frac{1}{n}\left( \sum_{i=k-n+1}^{k-1} B_i +\Tilde{\lambda}_k \right) &\quad \text{stand-still} \\
B_{k-1} &\quad \text{else}
\end{cases}
\end{equation}

Each measurement is corrected by

\begin{equation}
\label{eqn:biasupdate}
\lambda_k = \Tilde{\lambda}_k - B_k
\end{equation}

The paper refers to the estimated value after correction as sensor measurement.

To minimize integration errors and enable bias estimation updates, a Zero Velocity Update ZVU is implemented. Stand-still is detected using a combined acceleration and angular rate detector [25]. If a stand-still is detected, $a_{x,k}^{(\text{IMU})}$, $a_{y,k}^{(\text{IMU})}$, $\omega_{z,k}^{(\text{IMU})}$, and $v_k$ are set to zero.

§.§ Complete Particle Filter Setup

In cases of linear state-space models, Kalman filters are proven to be the optimal filter. In this project, the motion model is linear; however, the observation model is highly non-linear. Therefore, we propose using a PF to compute the position estimate using the motion model and observation model as described above. The complete algorithm is found in <Ref>.

Particle filter with map matching

[1]
Initialize particles $\hat{\mathbf{x}}_0^{(j)} \sim p(\hat{\mathbf{x}}_0)$ and particle weights $ w_0^{(j)} = \frac{1}{N}  $ for $j = 1, \dots, N$ particles
$k = 1, 2, \dots$
ZVU = 'true'
Bias estimation update (<ref>)
Set control input $u_k = 0$
Set measurements $a_{y,k} = \omega_{z,k} = 0$
Correct IMU measurements (<ref>), (<ref>)
$j = 1, \dots, N$
Particle diffusion according to (<ref>)
Map onto 2D plane using DTM: (<ref>)-(<ref>)
Compute predicted measurements $\hat{\mathbf{z}}_k^{(j)}$ using (<ref>)
Compute weights: $w_k^{(j)} \propto p(\mathbf{z}_k|\hat{\mathbf{z}}_k^{(j)})$
Normalize weights: $\tilde{w}_k^{(j)} = \frac{w_k^{(j)}}{\sum_n w_k^{(n)}}$

Resample particles: sample $M_p$ particles with replacement from the current set of particles, using the weights $\tilde{w}_k^{(j)}$

Compute $\hat{\mathbf{p}}_k$ (<ref>)

To avoid particle depletion, the particles are resampled when the effective sample size $N_{\mathrm{eff}}$ falls below a threshold $N_{\mathrm{th}} \in [1\, \, N]$. The effective sample size is commonly chosen as

\begin{equation}
N_{\mathrm{eff}} = \frac{1}{\sum_j\left( w_k^{(j)} \right)^2}
\end{equation}

Depending on the measurements, other resampling strategies may be sensible. If no absolute position measurements are available we suggest to resample when the train enters and travels on a curved track segment.
The final position estimate is calculated using

\begin{equation}
\label{eqn:finalestimate}
\hat{\mathbf{p}}_k = \arg \  \underset{\mathbf{p}_k}{\max} \ p(\mathbf{p}_k|\mathbf{z}_{1:k},\mathbf{u}_{1:k})
\end{equation}

where $\hat{\mathbf{p}}_k = \begin{bmatrix}
\hat{p}_{x,k} & \hat{p}_{y,k}
\end{bmatrix}^T$ is a sub-vector of $\hat{\mathbf{z}}_k$.

§ EXPERIMENTAL RESULTS

Comparison of different methods during a short GNSS outage

Complete mapped test track

Velocity profile of the test drive

For a proof-of-concept, we use the open Rail vehicle positioning data set, Lucy 2018 [22]. The data set contains 120 km of conventional and secondary lines, including tight curves, steep slopes, and forested embankments. A 3D geometrical track map with centimeter-level precision is available on 24.2 km of the track. The dataset contains 6 DOF IMU measurements and position data from two GNSS receivers. Since the positioning method presented in this paper heavily relies on the track map, we use the section of the dataset for which the 3D track map is available.

Since no ground truth is provided with the data set, the results in this paper are compared to the positions computed by the high-end GNSS receiver by Thales mounted on the test train. To verify its accuracy, the across-track error of each measurement is calculated, showing that 95% of all measurements lie within one meter from the track. Position updates are available every second with no outages during the test ride. To further enhance accuracy and establish what is referred to as ground truth throughout the remaining paper, the Thales GNSS measurements are corrected by projection onto the track map. The other receiver is then used to test our proposed method.

Absolute position error

Error CDF

All results are computed using our proposed method and compared to an EKF. The used data set includes an EKF based positioning solution, using tangential acceleration $a_x^{(\text{IMU})}$ and yaw angle turn rate $\omega_z^{(\text{IMU})}$ as input to a CTRA motion model. GNSS measurements of position and velocity are used in the measurement update. We combine this solution with the approach given in [11] to include map matching into the algorithm, which in the remaining paper will be called EKFMM. The parameter settings are adjusted to the sensor specification and are listed for both methods in <Ref>. The PF uncertainties are increased compared to the EKFMM to prevent particle depletion.

Straight passage from the track map

Curved passage of the track map

Curvature of straight track segment

Curvature of curved track segment

Euclidean position error at straight track segment

Euclidean position error at curved track segment

Performance evaluation during indefinite GNSS outage on a straight and on a curved path segment

Parameter values used in both filters

1c|     $\sigma_{a_x}^{(\text{IMU})}$     $\sigma_{a_y}^{(\text{IMU})}$     $\sigma_{\omega_z}^{(\text{IMU})}$     $\sigma_{\text{bias}}$
1c|PF     0.01 g     0.01 g     0.2 $^{\circ}$/s     5e-6 m/s$^2$

1c|EKFMM     0.005 g     -     0.05 $^{\circ}$/s     -

1c|     $\sigma_{p_x}^{(\text{GNSS})}$     $\sigma_{p_y}^{(\text{GNSS})}$     $\sigma_v^{(\text{GNSS})}$     $\sigma^{(\text{map})}$
1c|PF     2.45 m*     4.13 m*     0.4 m/s*     -

1c|EKFMM     2.04 m*     3.45 m*     0.35 m/s*     0.01 m

* Uncertainty updated dynamically according to GNSS receiver statistics. The displayed value is the mean uncertainty

Two scenarios are presented. First, the method is tested and evaluated on experimental data given in the data set. Second, the experimental data is altered by removing GNSS data after 50 s. In both scenarios, 1000 particles are used. The value was found empirically to balance computational cost and filter convergence.

§.§ Evaluation using experimental data

In this section, the proposed method is tested on experimental data. The map of the test track is shown in <Ref>.

The test train's starting position and travel direction are marked in <Ref>, and the velocity profile is shown in <Ref>. The recorded data includes stops, varying speeds, and several GNSS outages along the journey, the longest being 35 s. A representative example of the filter's performance can be seen in <Ref>, which shows position estimates and estimated $3\sigma$ confidence interval of the raw GNSS data, EKFMM solution and the proposed PF based solution, predicted according to the respective model. The confidence intervals of the EKFMM solution quickly increase during the outage, while the PF solution can maintain the position estimate with stable uncertainty.
The error statistics are shown in <Ref> and <Ref>. Both the EKFMM and the PF perform similarly when GNSS is availble, but during GNSS outages marked in grey, the PF solution outperforms the EKFMM. For the PF solution, the $3\sigma$ error lies at 11.3 m, whereas for the EKFMM solution, the $3\sigma$ error is 17.71 m.

§.§ Indefinite GNSS outage on different track characteristics

This section examines the performance of the PF solution during an indefinite GNSS outage. The experimental data as in <Ref> is used, but the GNSS signal is removed after 50 s, simulating an indefinite outage. The performance is evaluated on two sections of the test track with different track characteristics. <Ref> shows a segment consisting of several straight sections. A plot of the track curvature is shown in <Ref>. In contrast, <Ref> shows a segment consisting of several turns and short straight sections. The corresponding curvature is given in <Ref>. IMU measurements and track coordinates are as provided by the dataset, but GNSS reception stops after 50 s. We compare the results using the proposed method to the EKFMM by computing the absolute position error, see <Ref>. It can be seen that the PF solution performance is dependent on the track characteristic. The absolute position error increases on straight track sections. During curves, the accuracy reached by the proposed method is comparable to the accuracy when GNSS is available. During straight sections, both methods need to rely on the accelerometer measurements and diverge. However, it can be seen that the PF, in contrast to the EKFMM can recover to high accuracy during short curves. In <Ref> after 190 s, the absolute error does not increase further. This is due to four shorter stand-still phases and the resulting sensor bias update.

§ CONCLUSION

In this paper, we have proposed a minimalist approach for train localization during GNSS outages by utilizing track geometry and IMU sensor data. A discrete track map has been formulated as a look-up table and used in different ways: to map the 1D position estimate of the PF onto the 2D plane. Also, geometry features such as curvature have been used in the measurement update with gyroscope measurements. By this, the position estimation in the PF has been reduced to one dimension, resulting in lower computational complexity of the PF compared to 2D or 3D tracking filters. The formulation as a LUT ensures great flexibility and adaptability of the approach. Map coordinates do not need to be placed equidistantly. Furthermore, any interpolation method between map coordinates can be implemented and computed offline. Additional features and measurements can be added easily.

The proposed method was tested on experimental data from the open Lucy 2018 data set. The results show that the PF outperforms the EKFMM approach during shorter GNSS outages and providing similar accuracy when GNSS is available. We simulated an indefinite GNSS outage and tested the filter performance on two segments of the test track with different characteristics. It could be shown that the filter can provide an accurate positioning solution during track curves and curvy railway lines. However, the performance is dependent on the track characteristic and does, therefore, not work as a stand-alone positioning solution, when only using the IMU.

The key strength of the approach lies in its minimalist sensor usage and flexible track map incorporation. To meet the high safety standards for train positioning solutions, redundant solutions must be implemented. The presented approach does not use sensors incorporated in recent train positioning solutions, such as odometer data or magnetometers. If knowledge of the starting position is available, even the irregular GNSS measurements can be removed. Therefore, the proposed method is practical to include as redundant positioning method, increasing safety of the positioning solution. Another possible application is an implementation in a classical GNSS/IMU sensor fusion approach to increase positioning accuracy in challenging environments. Further research will include testing in different scenarios and an improved bias correction.

[1]
H. A. Hamid, G. L. Nicholson, and C. Roberts, “Impact of train positioning inaccuracies on railway traffic management systems: framework development and impacts on tms functions,” IET Intelligent Transport Systems, vol. 14, no. 6, pp. 534–544, 2020.

[2]
F. Gustafsson, Statistical Sensor Fusion.
Studentlitteratur, 2010.

[3]
“EN50126 - Railway applications - The specification and demonstration of Reliability, Availability, Maintainability, and Safety (RAMS),” Eur. Commitee for Electrotechnical Standardization Std.

[4]
C. Reimer, F. Müller, and E. Hinüber, “INS/GNSS/odometer data fusion in railway applications,” in 2016 DGON Intertial Sensors and Systems (ISS), pp. 1–14, IEEE, 2016.

[5]
J. Otegui, A. Bahillo, I. Lopetegi, and L. E. Díez, “Evaluation of experimental GNSS and 10-DOF MEMS IMU measurements for train positioning,” IEEE Transactions on Instrumentation and Measurement, vol. 68, no. 1, pp. 269–279, 2018.

[6]
O. Heirich, “Bayesian train localization with particle filter, loosely coupled GNSS, IMU, and a track map,” Journal of Sensors, vol. 2016, 2016.

[7]
J. Liu, B. Cai, and J. Wang, “Track-constrained GNSS/odometer-based train localization using a particle filter,” in 2016 IEEE Intelligent Vehicles Symposium (IV), pp. 877–882, IEEE, 2016.

[8]
D. Liu, W. Jiang, B. Cai, O. Heirich, J. Wang, and W. Shangguan, “Robust train localisation method based on advanced map matching measurement-augmented tightly-coupled GNSS/INS with error-state UKF,” The Journal of Navigation, pp. 1–24, 2023.

[9]
H. Winter, V. Willert, and J. Adamy, “Increasing accuracy in train localization exploiting track-geometry constraints,” in 2018 21st International Conference on Intelligent Transportation Systems (ITSC), pp. 1572–1579, IEEE, 2018.

[10]
O. Heirich, P. Robertson, and T. Strang, “RailSLAM-localization of rail vehicles and mapping of geometric railway tracks,” in 2013 IEEE International Conference on Robotics and Automation, pp. 5212–5219, IEEE, 2013.

[11]
W. Löffler and M. Bengtsson, “Using probabilistic geometrical map information for train localization,” in 2022 25th International Conference on Information Fusion (FUSION), pp. 01–08, IEEE, 2022.

[12]
W. Löffler and M. Bengtsson, “Evaluating the impact of map inaccuracies on path discrimination behind railway turnouts,” in 2022 IEEE 95th Vehicular Technology Conference: (VTC2022-Spring), pp. 1–5, IEEE, 2022.

[13]
C. Hasberg, S. Hensel, and C. Stiller, “Simultaneous localization and mapping for path-constrained motion,” IEEE Transactions on Intelligent Transportation Systems, vol. 13, no. 2, pp. 541–552, 2011.

[14]
H. Winter, S. Luthardt, V. Willert, and J. Adamy, “Generating compact geometric track-maps for train positioning applications,” in 2019 IEEE Intelligent Vehicles Symposium (IV), pp. 1027–1032, IEEE, 2019.

[15]
A. Andersson, “Compact digital track maps: Enhancing train traveller information at the crossing of accuracy and availability: A comparative analysis of algorithms for generating compact representations of railway tracks,” No. TRITA-EECS-EX-2023:712, (KTH Royal Institute of Technology), 2023.

[16]
S. Spinsante and C. Stallo, “Hybridized-GNSS approaches to train positioning: Challenges and open issues on uncertainty,” Sensors, vol. 20, no. 7, p. 1885, 2020.

[17]
N. Zhu, J. Marais, D. Bétaille, and M. Berbineau, “GNSS position integrity in urban environments: A review of literature,” IEEE Transactions on Intelligent Transportation Systems, vol. 19, no. 9, pp. 2762–2778, 2018.

[18]
P. d’Harcourt, J.-P. Michel, L. Poletti, and S. Glevarec, “Innovative and efficient inertial navigation system for train localization in GNSS-denied environments,” in 2023 DGON Inertial Sensors and Systems (ISS), pp. 1–18, IEEE, 2023.

[19]
O. Heirich, B. Siebler, and E. Hedberg, “Study of train-side passive magnetic measurements with applications to train localization,” Journal of Sensors, vol. 2017, 2017.

[20]
B. Kröper, M. Lauer, and M. Spindler, “Using the ferromagnetic fingerprint of rails for velocity estimation and absolute localization of railway vehicles,” in 2020 European Navigation Conference (ENC), pp. 1–10, IEEE, 2020.

[21]
M.-A. Lavoie and J. R. Forbes, “Map-aided train navigation with IMU measurements,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 3465–3470, IEEE, 2021.

[22]
H. Winter, “Rail vehicle positioning data set: Lucy, october 2018,” 2020.

[23]
X. R. Li and V. P. Jilkov, “Survey of maneuvering target tracking. Part I. Dynamic models,” IEEE Transactions on aerospace and electronic systems, vol. 39, no. 4, pp. 1333–1364, 2003.

[24]
H. Carlsson, I. Skog, and J. Jaldén, “Self-calibration of inertial sensor arrays,” IEEE Sensors Journal, vol. 21, no. 6, pp. 8451–8463, 2021.

[25]
J. Callmer, D. Törnqvist, and F. Gustafsson, “Probabilistic stand still detection using foot mounted IMU,” in 2010 13th International Conference on Information Fusion, pp. 1–7, IEEE, 2010.


# An Observer-Based Composite Identifier for Online Estimation of the Thévenin
Equivalent Parameters of a Power System

Daniele Zonetti, Romeo Ortega, Rafael Cisneros, Alexey Bobtsov, Fernando
Mancilla-David, Oriol Gomis-Bellmunt This work was supported by the Ministry
of Science and Higher Education of Russian Federation, passport of goszadanie
no. 2019-0898 and by FEDER/Ministerio de Ciencia, Innovación y Universidades-
Agencia Estatal de Investigación, Project RTI2018-095429-B-I00. The work of O.
Gomis-Bellmunt is supported by the ICREA Academia program.D. Zonetti and O.
Gomis-Bellmunt are with the Centre d’Innovació Tecnològica en Convertidors
Estàtics i Accionaments, Departament d’Enginyeria Elèctrica, Universitat
Politècnica de Catalunya, Barcelona 08028, Spain. (e-mail:
daniele.zonetti(oriol.gomis)@upc.edu). R. Ortega and R. Cisneros are with the
Departamento Académico de Sistemas Digitales, ITAM, Rio Hondo 1, Col. Progreso
Tizapan, 01080 Ciudad de México, Mexico (email:
romeo.ortega(rcisneros)@itam.mx).A. Bobtsov is with the Department of Control
Systems and Robotics, ITMO University, Kronverkskiy av. 49, Saint-Petersburg,
197101, Russia (e-mail: bobtsov@mail.ru).F. Mancilla–David is with the
Department of Electrical Engineering, University of Colorado Denver, Denver,
Colorado 80204, USA, (e-mail<EMAIL_ADDRESS>
###### Abstract

We consider a Thévenin equivalent circuit capturing the dynamics of a power
grid as seen from the point of common coupling with a power electronic
converter, and provide a solution to the problem of online identification of
the corresponding circuit parameters. For this purpose, we first derive a
linear regression model in the conventional $\tt abc$ coordinates and next
design a bounded observer-based composite identifier that requires local
measurements and knowledge of the grid frequency only. An extension that
guarantees exponential convergence of the estimates, under the additional
assumption of knowledge of the grid X/R ratio, is further provided. The
performance of the proposed identifier, which subsumes a conventional gradient
descent algorithm, is illustrated via detailed computer simulations.

##  I Introduction

With the widespread penetration of power electronic converters in the existing
power systems, the problem of establishing an accurate model that is
representative of the grid dynamics for the purpose of analysis and of control
design is timely and increasingly relevant [9]. The conventional approach is
to use a Thévenin equivalent (TE) circuit to represent the overall grid as
seen from the individual bus where the power converter is interfaced. The
theoretical basis standing behind this concept is provided by the TE theorem,
stating in its original version that any linear circuit constituted by current
and/or voltage sources and resistances can be described via an equivalent
circuit characterized by a voltage source combined in series with a
resistance—a formulation that has been later extended to single-phase and
three-phase AC circuits, with the notion of resistance being replaced by the
notion of impedance [8]. This representation is suited for the analysis of a
variety of key issues for a correct and safe operation of the power converter,
ranging from fault studies [13], voltage stability analysis [10], loadability
limit computation [11] to tuning of the related controllers [12]. In these
studies it is emphasized that, to make the problem tractable in practice, the
identification of the TE parameters shall be accomplished using only local
measurements of voltage and currents.
The problem of online estimation of the TE parameters has been widely studied
in the literature mainly using least-squares or extended Kalman filter
techniques. As is well-known, the excitation requirements of these estimators,
namely persistent excitation [7], is quite stringent and hard to satisfy
(without probing signals) in the current application. Another disadvantage of
these methods is that—if the excitation conditions are satisfied—the
adaptation gain converges to zero losing the alertness of the estimator. In
[1] a variation of the least-squares method proposed 32 years ago in [4], that
converges in finite time, was used for this problem. Although the excitation
conditions are weak the algorithm—that in its initial stage is akin to an off-
line estimator—involves a numerically sensitive matrix inversion and, as it
converges to a standard least-squares, loses its alertness.

The contributions of the paper are as follows. In Section II we discuss the
fundamental assumptions behind the formulation of the mathematical model of a
power converter interfaced to a TE circuit. Then, in Section III we derive a
LRE that is instrumental for the design of a bounded observer-based composite
identifier, a task that is carried-out in Section IV. By employing a further
assumption, we next generate in Section LABEL:sec5 an alternative, reduced
LRE, for which a conventional gradient descent algorithm can be applied to
ensure exponential convergence. The usefulness of the theoretical results are
illustrated via simulations in Section LABEL:sec6. We conclude the paper in
Section VII with some final remarks and guidelines for future investigation.

Notation. Given a vector $a\in\mathds{R}^{n}$, we denote the square of the
Euclidean norm as $|a|^{2}:=a^{\top}a$. The symbol $\mathsf{1}_{3}$ denotes a
three-dimensional vector of ones. Given a differentiable signal
$u(t)\in\mathds{R}^{r}$, we define the derivative operator ${du(t)\over
dt}=:{\mathfrak{p}}[u(t)]$ and denote the action of an LTI filter
$F({\mathfrak{p}})\in\mathds{R}({\mathfrak{p}})$ as $F({\mathfrak{p}})[u(t)]$.
For given $\omega>0$, $\varphi\in\mathbb{R}$, we denote the vector

$\mathcal{S}_{\varphi}(t):={\sqrt{\frac{2}{3}}}\begin{bmatrix}\sin(\omega
t+\varphi)\\\ \sin(\omega t+\varphi-\frac{2}{3}\pi)\\\ \sin(\omega
t+\varphi+\frac{2}{3}\pi)\end{bmatrix}.$

##  II Assumptions

The mathematical model, in $\tt abc$ reference frame, of the TE of a three-
phase, symmetrical configured system—see Fig. 1 for the corresponding circuit
schematic—is given by111A three-phase AC electrical system is said to be
symmetrically configured if a symmetrical feeding voltage yields a symmetrical
current and vice versa [8].

$L\frac{d}{dt}i=-Ri+v-e,$ (1)

where: $i(t)\in\mathbb{R}^{3}$, $e(t)\in\mathbb{R}^{3}$ denote the three–phase
current and equivalent voltage of the grid, respectively;
$v(t)\in\mathbb{R}^{3}$ denotes the three-phase voltage synthesized by the
converter at the point of common coupling (PCC); $L\in\mathbb{R}_{>0}$ and
$R\in\mathbb{R}_{>0}$ denote respectively the grid equivalent inductance and
resistance.

Figure 1: Single-line circuit schematic of a power converter connected to a TE
grid.

######  Assumption 1 (Grid voltage)

The TE voltage source is described by a three-phase balanced purely sinusoidal
signal:

$e=E\mathcal{S}_{0},$ (2)

with $E>0$.

######  Assumption 2 (Parameters)

The frequency $\omega>0$ is known and the positive parameters $L$, $R$, $E$
are unknown constants.

######  Assumption 3 (Measurements)

The three-phase signals $i(t)\in\mathds{R}^{3}$ and $v(t)\in\mathds{R}^{3}$
are measurable.

Assumption 1 is justified for AC grids that are characterized by a
sufficiently large number of synchronous rotating machines, providing large
inertia, and/or whenever a sufficiently tight regulation of the frequency is
implemented via grid-forming controllers. Assumption 2 can be instead enforced
via the design of a suitably defined phase-locked loop (PLL), which allows to
recover the exact value of the grid frequency at a faster time-scale.
Assumption 3 is quite natural and verified in all practical scenarios, as
those measurements are utilized for the control of the converter.

For an appropriate modeling of the three-phase voltage synthesized by the
converter at the PCC a caveat is necessary. In order to guarantee a correct
operation of the power system, the state and input variables $\\{i,v\\}$ are
required to converge to suitable balanced, three-phase AC signal with the
common grid frequency $\omega$, but different amplitudes and phase shifts.
These operating conditions are enforced by the converter control algorithms,
which are usually based on a $\tt dq$ reference frame ensuring that the
steady-states of interest, in new coordinates, correspond to constant
quantities. Instrumental for the definition of such reference frame is the
knowledge of the grid frequency $\omega$. Indeed, using this information the
overall system composed by the TE circuit and the power converter—transformed
in $\tt dq$ coordinates—is linear time-invariant (LTI) and adequate control
solutions can be established—see for example [6]. The design of an
exponentially stabilizing controller legitimizes then the following
assumption.

######  Assumption 4 (PCC voltage)

The voltage at the PCC is described by the three-phase signal

$v=V\mathcal{S}_{\phi}+{\boldsymbol{\epsilon}}_{t},$ (3)

with $V>0$, $\phi\in(-\frac{\pi}{2},\frac{\pi}{2})$, where
${\boldsymbol{\epsilon}}_{t}\in\mathds{R}^{3}$ is a signal exponentially
decaying to zero.

###### Remark 1

A different steady-state amplitude $V$, phase $\phi$ and exponential term
${\boldsymbol{\epsilon}}_{t}$ is triggered in (3) any time a change occurs in
the parameters of the TE or in the references provided to the power converter
controllers. The rapidity at which the signal ${\boldsymbol{\epsilon}}_{t}$
vanishes is determined by the tuning of the converter control algorithms.

##  III Derivation of a LRE

Following the standard procedure for parameter estimation [3, 7], we derive a
LRE for the system (1), as explained in the following lemma.

###### Lemma 1

Consider the system (1) verifying Assumptions A1-A4. There exist measurable
signals ${\bf Z}(t)\in\mathds{R}^{3}$ and
${\boldsymbol{\Psi}}_{f}(t)\in\mathds{R}^{3\times 3}$ such that the following
LRE is satisfied

${\bf Z}={\boldsymbol{\Psi}}_{f}\theta+{\boldsymbol{\epsilon}}_{t},$ (4)

where

$\theta:=\mbox{col}\left({R\over L},{1\over L},{E\over
L}\right)\in\mathds{R}^{3},$ (5)

and ${\boldsymbol{\epsilon}}_{t}\in\mathds{R}^{3}$ is a signal exponentially
decaying to zero.

###### Proof:

The equation (1) may be written as

$\frac{d}{dt}i={\boldsymbol{\Psi}}\theta,$ (6)

where we defined the three-dimensional square matrix

${\boldsymbol{\Psi}}:=\begin{bmatrix}-i&|&v&|&-\mathcal{S}_{0}\end{bmatrix}.$
(7)

The proof is completed applying to (6) the LTI, stable filter

$F({\mathfrak{p}})={\lambda\over{\mathfrak{p}}+\lambda},$

with $\lambda>0$ and defining

$\Box\Box\Box$

###### Remark 2

Since we assumed that the power system configuration is symmetric, the sum of
three phases of both the grid current and the voltage at the PCC is always
zero, that is:

$\mathsf{1}_{3}^{\top}i_{g}(t)=\mathsf{1}_{3}^{\top}v(t)=\mathsf{1}_{3}^{\top}\mathcal{S}_{0}(t)=0,\quad\forall
t\geq 0.$ (8)

A consequence of this fact is that the regression matrix ${\boldsymbol{\Psi}}$
satisfies ${\boldsymbol{\Psi}}\mathsf{1}_{3}=0$, hence it is singular and
cannot satisfy the persistent excitation requirement

$\int_{t}^{t+T}{\boldsymbol{\Psi}}(s){\boldsymbol{\Psi}}^{\top}(s)ds\geq\delta
I_{3},\;\forall t\geq 0,$ (9)

and some $T>0$ and $\delta>0$. Note also that in view of (8), only two of the
three $\tt abc$ phases are required for the construction of such matrices.

##  IV An Observer-Based Composite Identifier

Clearly, a classical gradient descent (or least squares) algorithm can be
applied to the LRE (4). However, in view of the lack of excitation indicated
in Remark 2, poor performances are observed in simulations, that cannot be
further improved by appropriate tuning of the estimator.

In view of this situation we propose instead the use of the composite
identifier of [5], which offers additional freedom in the design at the
expense of an higher-order dynamics.

###### Proposition 1

Consider the system (1) and the LRE (4). Define the observer-based composite
identifier:

$\displaystyle\dot{\hat{i}}$
$\displaystyle=-\alpha(\hat{i}-i)+{\boldsymbol{\Psi}}\hat{\theta}$ (10)
$\displaystyle\dot{\hat{\theta}}$
$\displaystyle=-\gamma_{P}{\boldsymbol{\Psi}}^{\top}(\hat{i}-i)+\gamma_{I}{\boldsymbol{\Psi}}_{f}^{\top}({\bf
Z}-{\boldsymbol{\Psi}}_{f}\hat{\theta}),$

with tuning gains $\alpha>0$, $\gamma_{P}>0$ and $\gamma_{I}>0$. Then the
system trajectories are globally bounded and

$\lim_{t\to\infty}\hat{i}(t)=i(t).$ (11)

###### Proof:

Let us define the incremental variables

$\tilde{i}:=\hat{i}-i,\quad\tilde{\theta}:=\hat{\theta}-\theta.$

Using (6) and (4) we obtain the error dynamics:

Moreover, if $\phi\neq 0$, the trajectories generated by the gradient descent
estimator

$\dot{\hat{\vartheta}}=\gamma{\boldsymbol{\Psi}}_{f,\tt ab}^{\top}({\bf
Z}_{\tt ab}-{\boldsymbol{\Psi}}_{f,\tt ab}\hat{\vartheta}),$ (14)

with tuning gain $\gamma>0$, verify

$\lim_{t\to\infty}\hat{\vartheta}(t)=\vartheta,\quad(exp).$

###### Proof:

Let us consider the phase $\tt a$ and $\tt b$ from equation (1), so that we
can write222Phases $a$ and $b$ are selected with no loss of generality, since
a similar proof applies for a different choice of the phases.

$\frac{d}{dt}i_{\tt ab}+\frac{\omega}{\rho}i_{\tt ab}={\boldsymbol{\Psi}}_{\tt
ab}\vartheta,$ (15)

where we defined the two-dimensional square matrix

${\boldsymbol{\Psi}}_{\tt ab}:=\begin{bmatrix}v_{\tt
ab}&|&-\mathcal{S}_{0,\tt{ab}}\end{bmatrix}.$ (16)

The first part of the proof is then completed, similar to Lemma 1, by applying
to (15) the LTI, stable filter

$F({\mathfrak{p}})={\lambda\over{\mathfrak{p}}+\lambda},$

with $\lambda>0$ and defining

Figure 2: Estimates of the grid resistance, inductance and voltage amplitude
(solid lines) obtained via the observer-based composite identifier (10), as
compared to the corresponding, actual values (dotted lines), and following: a)
a drop of the SCR from $3$ to $1.5$ at $T=1$ s (left); b) a drop of $10\%$ of
the TE voltage source at $2T$ s (right).

We evaluate the performance of the identifier following two types of
perturbations occurring respectively at time $T=1$ s and $2T$ s. The first
type of perturbation consists in a simultaneous change of the resistance and
inductance of the TE, following a drop in the SCR from $3$ to $1.5$. The
second type of perturbation takes the form of a drop of $10\%$ of the TE
voltage source. Extensive simulations are realized for these two relevant
scenarios, with the gains of the identifier (10) set to fixed
$\lambda=10^{3}$, $\gamma_{I}=10^{8}$ and variable $\alpha\in[0\;10^{6}]$,
$\gamma_{P}\in[0\;10^{8}]$. The obtained results for both scenarios and three
illustrative gain pairs $(\alpha,\gamma_{P})$ are reported in Fig. IV, where
it is shown that for all such gains the identifier allows to correctly recover
the actual values of the corresponding parameters. Unsurprisingly, we further
observe that acceptable performances strongly depend on an appropriate
selection of the gains of the identifier. A particular case is given whenever
we pick $\alpha=\gamma_{P}=0$, that is, the identifier coincides with a
gradient descent (GD) algorithm, see also Remark LABEL:rem2. We observe indeed
that with this design, both a change of the impedance and of the voltage of
the TE generate large peaks in the estimates provided by the identifier,
resulting in convergence times superior to $350$ ms—performances that can be
neatly improved with a different tuning of the gains. We thus conclude that
the additional degrees of freedom provided by the observer-based composite
identifier can be exploited to improve performances under all type of
perturbations. We complete this section briefly illustrating in Fig. IV the
performances of a gradient descent algorithm applied to the reduced LRE
(LABEL:lre-ab), with tuning gain $\gamma=10^{8}$. In addition to the
previously considered scenario, we evaluate the responses of the algorithm in
case that a change of the X/R ratio, initially set to $\rho(0)=5$, may occur
concurrently with the change of the SCR at time $T=1$ s, that is
$\rho(t)=\\{3,5,7\\}$, for $t\geq T$. It is shown that, as long as the X/R
ratio remains unchanged, the estimates fastly and asymptotically converge to
their actual values—a fact that stems from Fig. 4 where we plot

$\lambda_{\min}\Big{\\{}\int_{0}^{t}{\boldsymbol{\Psi}}_{f,\tt{ab}}(s){\boldsymbol{\Psi}}^{\top}_{f,\tt{ab}}(s)ds\Big{\\}},$

which grows to infinity. As is well-known [2, Proposition 4] the latter is a
necessary condition for global convergence of the gradient estimator.
Nevertheless, whenever the ratio $\rho$ changes at $T=1\;s$, Assumption A5 is
violated and large estimation errors are observed, suggesting that caution
must be taken in employing such solution.

Figure 3: Estimates of the grid resistance, inductance and voltage amplitude
(solid lines) obtained via a classical gradient descent algorithm applied to
the LRE (LABEL:lre-ab), as compared to the corresponding, actual values
(dotted lines), and following: a) a simultaneous change of the SCR from $3$ to
$1.5$ and of the ratio $\rho$ from $5$ to $\\{3,5,7\\}$ at $T=1$ s; b) a drop
of $10\%$ of the TE voltage source at $2T$ s.

Figure 4: Minimum eigenvalue of the matrix
$\int_{0}^{t}{\boldsymbol{\Psi}}_{f,\tt{ab}}(s){\boldsymbol{\Psi}}^{\top}_{f,\tt{ab}}(s)ds$.

##  VII Conclusions

In this paper we have addressed the problem of online identification of the
parameters of a Thévenin equivalent grid, where this is adopted to describe a
power system interfaced with a power electronic converter in closed-loop with
an exponentially stabilizing controller. Based on this setting, we have
derived a linear regression model for the system and next designed an
observer-based composite identifier ensuring that the resulting estimates are
bounded for all positive gains. An alternative, exponentially converging
design is further proposed, assuming either a dominantly inductive grid or the
a priori knowledge of the X/R ratio—conditions that allow to generate a
reduced linear regressor equation. The theoretical results are validated via
simulations, which suggest that convergence and suitable performances can be
enforced by an appropriate tuning of the gains. Future works will explore
solutions able to guarantee asymptotic convergence of the estimates to their
actual values by lifting Assumption A5. Validation of the obtained results on
a detailed benchmark and related experiments are under progress and will be
reported soon.

## References

  * [1] A. Arancibia, C. A. Soriano-Rangel, F. Mancilla-David, R. Ortega and K. Strunz, Finite-time Identification of the Thevenin Equivalent Parameters in Power Grids, International Journal of Electrical Power & Energy Systems, vol. 116, March 2020.
  * [2] N. E. Barabanov and R. Ortega, On global asymptotic stability of $\dot{x}=\phi(t)\phi^{\top}(t)x$ with $\phi(t)$ bounded and not persistently exciting, Systems and Control Letters, vol. 109, pp. 24-27, 2017.
  * [3] L. Ljung, System Identification: Theory for the User, Prentice Hall, New Jersey, 1987.
  * [4] R. Ortega, An online Least-squares Parameter Estimator with Finite Convergence Time, Proc. IEEE, Vol. 76, No. 7, July 1988.
  * [5] E. Panteley, R. Ortega and P. Moya, Overcoming the obstacle of detectability in certainty equivalent adaptive control, Automatica, vol. 38, no. 7, 2002.
  * [6] D. Zonetti, A., Bobtsov, R., Ortega, N., Nikolaev, N. and O. Gomis-Bellmunt. An Almost Globally Stable Adaptive Phase-Locked Loop for Synchronization of a Grid-Connected Voltage Source Converter. arXiv preprint arXiv:2201.05490.
  * [7] S. Sastry and M. Bodson, Adaptive Control: Stability, Convergence and Robustness, Prentice-Hall, New Jersey, 1989.
  * [8] J. Svoboda, and R. Dorf, Introduction to electric circuits John Wiley & Sons, 2013.
  * [9] O. Gomis-Bellmunt, J. Song, M. Cheah-Mane, E. Prieto-Araujo, Steady-state impedance mapping in grids with power electronics: What is grid strength in modern power systems?, International Journal of Electrical Power & Energy Systems, vol. 136, 2022.
  * [10] H. Sun et al., Review of Challenges and Research Opportunities for Voltage Control in Smart Grids, IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 2790-2801, July 2019.
  * [11] L. Peng, J. Zhao, Y. Tang, L. Mili, Z. Gu, Z. Zheng, Real-Time LCC-HVDC Maximum Emergency Power Capacity Estimation Based on Local PMUs, IEEE Transactions on Power Systems, vol. 36, no. 2, pp. 1049-1058, March 2021.
  * [12] A. Raghami, G. Ledwich, Y. Mishra, Improved Reactive Power Sharing Among Customers’ Inverters Using Online Thévenin Estimates, IEEE Transactions on Power Systems, vol. 34, no. 6, pp. 4168-4176, November 2019.
  * [13] S. Shen, D. Lin, H. Wang, P. Hu, K. Jiang, D. Lin, B. He, An Adaptive Protection Scheme for Distribution Systems With DGs Based on Optimized Thevenin Equivalent Parameters Estimation, IEEE Transactions on Power Delivery, vol. 32, no. 1, pp. 411-419, February 2017.
  * [14] K. S. Narendra and A. M. Annaswamy, Stable Adaptive Systems, Prentice-Hall, 1989.


# Environmental variability and network structure determine the optimal
plasticity mechanisms in embodied agents

Emmanouil Giannakakis1,2, Sina Khajehabdollahi 1 Anna Levina1,2,3

1Department of Computer Science, University of Tübingen, Tübingen, Germany
2Max Planck Institute for Biological Cybernetics, Tübingen, Germany
3Bernstein Center for Computational Neuroscience Tübingen, Tübingen, Germany

###### Abstract

The evolutionary balance between innate and learned behaviors is highly
intricate, and different organisms have found different solutions to this
problem. We hypothesize that the emergence and exact form of learning
behaviors is naturally connected with the statistics of environmental
fluctuations and tasks an organism needs to solve. Here, we study how
different aspects of simulated environments shape an evolved synaptic
plasticity rule in static and moving artificial agents. We demonstrate that
environmental fluctuation and uncertainty control the reliance of artificial
organisms on plasticity. Interestingly, the form of the emerging plasticity
rule is additionally determined by the details of the task the artificial
organisms are aiming to solve. Moreover, we show that co-evolution between
static connectivity and interacting plasticity mechanisms in distinct sub-
networks changes the function and form of the emerging plasticity rules in
embodied agents performing a foraging task.

## Introduction

One of the defining features of living organisms is their ability to adapt to
their environment and incorporate new information to modify their behavior. It
is unclear how the ability to learn first evolved Papini, (2012), but its
utility appears evident. Natural environments are too complex for all the
necessary information to be hardcoded genetically Snell-Rood, (2013) and more
importantly, they keep changing during an organism’s lifetime in ways that
cannot be anticipated Ellefsen, (2014); Dunlap and Stephens, (2016). The link
between learning and environmental uncertainty and fluctuation has been
extensively demonstrated in both natural Kerr and Feldman, (2003); Snell-Rood
and Steck, (2019), and artificial environments Nolfi and Parisi, (1996).

Nevertheless, the ability to learn does not come without costs. For the
capacity to learn to be beneficial in evolutionary terms, a costly nurturing
period is often required, a phenomenon observed in both biological Thornton
and Clutton-Brock, (2011), and artificial organisms Eskridge and Hougen,
(2012). Additionally, it has been shown that in some complex environments,
hardcoded behaviors may be superior to learned ones given limits in the
agent’s lifetime and environmental uncertainty Dunlap and Stephens, (2009);
Fawcett et al., (2012); Lange and Sprekeler, (2020).

The theoretical investigation of the optimal balance between learned and
innate behaviors in natural and artificial systems goes back several decades.
However, it has recently found also a wide range of applications in applied AI
systems Lee and Lee, (2020); Biesialska et al., (2020). Most AI systems are
trained for specific tasks, and have no need for modification after their
training has been completed. Still, technological advances and the necessity
to solve broad families of tasks make discussions about life-like AI systems
relevant to a wide range of potential application areas. Thus the idea of
open-ended AI agents Open Ended Learning Team et al., (2021) that can
continually interact with and adapt to changing environments has become
particularly appealing.

Many different approaches for introducing lifelong learning in artificial
agents have been proposed. Some of them draw direct inspiration from actual
biological systems Schmidhuber, (1987); Parisi et al., (2019). Among them, the
most biologically plausible solution is to equip artificial neural networks
with some local neural plasticity Thangarasa et al., (2020), similar to the
large variety of synaptic plasticity mechanisms Citri and Malenka, (2008);
Feldman, (2009); Caroni et al., (2012) that performs the bulk of the learning
in the brains of living organisms Magee and Grienberger, (2020). The
artificial plasticity mechanisms can be optimized to modify the connectivity
of the artificial neural networks toward solving a particular task. The
optimization can use a variety of approaches, most commonly evolutionary
computation.

The idea of meta-learning or optimizing synaptic plasticity rules to perform
specific functions has been recently established as an engineering tool that
can compete with state-of-the-art machine learning algorithms on various
complex tasks Burms et al., (2015); Najarro and Risi, (2020); Pedersen and
Risi, (2021); Yaman et al., (2021). Additionally, it can be used to reverse
engineer actual plasticity mechanisms found in biological neural networks and
uncover their functions Confavreux et al., (2020); Jordan et al., (2021).

Here, we study the effect that different factors (environmental fluctuation
and reliability, task complexity) have on the form of evolved functional
reward-modulated plasticity rules. We investigate the evolution of plasticity
rules in static, single-layer simple networks. Then we increase the complexity
by switching to moving agents performing a complex foraging task. In both
cases, we study the impact of different environmental parameters on the form
of the evolved plasticity mechanisms and the interaction of learned and static
network connectivity. Interestingly, we find that different environmental
conditions and different combinations of static and plastic connectivity have
a very large impact on the resulting plasticity rules.

## Methods

### Environment

We imagine an agent who must forage to survive in an environment presenting
various types of complex food particles. Each food particle is composed of
various amounts and combinations of $N$ ingredients that can have positive
(food) or negative (poison) values. The value of a food particle is a weighted
sum of its ingredients. To predict the reward value of a given resource, the
agent must learn the values of these ingredients by interacting with the
environment. The priors could be generated by genetic memory, but the exact
values are subject to change.

To introduce environmental variability, we stochastically change the values of
the ingredients. More precisely, we define two ingredient-value distributions
$E_{1}$ and $E_{2}$ Guttenberg, (2019) and switch between them, with
probability $p_{tr}$ for every time step. We control how (dis)similar the
environments are by parametrically setting $E_{2}=(1-2d_{e})E_{1}$, with
$d_{e}\in[0,1]$ serving as a distance proxy for the environments; when
$d_{e}=0$, the environment remains unchanged, and when $d_{e}=1$ the value of
each ingredient fully reverses when the environmental transition happens. For
simplicity, we take values of the ingredients in $E_{1}$ equally spaced
between -1 and 1 (for the visualization, see Fig. 3a, b).

### Static agent

The static agent receives passively presented food as a vector of ingredients
and can assess its compound value using the linear summation of its sensors
with the (learned or evolved) weights, see Fig. 1. The network consists of $N$
sensory neurons that are projecting to a single post-synaptic neuron. At each
time step, an input $X_{t}=(x_{1},\dots,x_{N})$ is presented, were the value
$x_{i},\ i\in\\{1,\dots,N\\}$ represents the quantity of the ingredient $i$.
We draw $x_{i}$ independently form a uniform distribution on the $[0,1]$
interval ($x_{i}\sim U(0,1)$). The value of each ingredient $w_{i}^{c}$ is
determined by the environment ($E_{1}$ or $E_{2}$).

The postsynaptic neuron outputs a prediction of the food $X_{t}$ value as
$y_{t}=g(WX_{t}^{T})$. Throughout the paper, $g$ will be either the identity
function, in which case the prediction neuron is linear, or a step-function;
however, it could be any other nonlinearity, such as a sigmoid or ReLU. After
outputting the prediction, the neuron receives feedback in the form of the
real value of the input $R_{t}$. The real value is computed as
$R_{t}=W^{c}X_{t}^{T}+\xi,$ where $W^{c}=(w_{1}^{c},\dots,w_{N}^{c})$ is the
actual value of the ingredients, and $\xi$ is a term summarizing the noise of
reward and sensing system $\xi\sim\mathcal{N}(0,\sigma)$.

Figure 1: An outline of the static agent’s network. The sensor layer receives
inputs representing the quantity of each ingredient of a given food at each
time step. The agent computes the prediction of the food’s value $y_{t}$ and
is then given the true value $R_{t}$; it finally uses this information in the
plasticity rule to update the weight matrix.

For the evolutionary adjustment of the agent’s parameters, the loss of the
static agent is the sum of the mean squared errors (MSE) between its
prediction $y_{t}$ and the reward $R_{t}$ over the lifetime of the agent. The
agent’s initial weights are set to the average of the two ingredient value
distributions, which is the optimal initial value for the case of symmetric
switching of environments that we consider here.

### Moving Agent

As a next step, we incorporate the sensory network of static agents into
embodied agents that can move around in an environment scattered with food. To
this end, we merge the static agent’s network with a second, non-plastic motor
network that is responsible for controlling the motion of the agent in the
environment. Specifically, the original plastic network now provides the agent
with information about the value of the nearest food. The embodied agent has
additional sensors for the distance from the nearest food, the angle between
the current velocity and the nearest food direction, its own velocity, and its
own energy level (sum of consumed food values). These inputs are processed by
two hidden layers (of 30 and 15 neurons) with $\tanh$ activation. The
network’s outputs are angular and linear acceleration, Fig. 2.

Figure 2: An outline of the network controlling the foraging agent. The sensor
layer receives inputs at each time step (the ingredients of the nearest food),
which are processed by the plastic layer in the same way as the static sensory
network, Fig. 1. The output of that network is given as input to the motor
network, along with the distance $d$ and angle $\alpha$ to the nearest food,
the current velocity $v$, and energy $E$ of the agent. These signals are
processed through two hidden layers to the final output of motor commands as
the linear and angular acceleration of the agent

The embodied agents spawn in a 2D space with periodic boundary conditions
along with a number of food particles that are selected such that the mean of
the food value distribution is $\sim 0$. An agent can eat food by approaching
it sufficiently closely, and each time a food particle is eaten, it is re-
spawned with the same value somewhere randomly on the grid (following the
setup of Khajehabdollahi et al., (2022)). After 5000 time steps, the
cumulative reward of the agent (the sum of the values of all the food it
consumed) is taken as its fitness. During the evolutionary optimization, the
parameters for both the motor network (connections) and plastic network
(learning rule parameters) are co-evolved, and so agents must simultaneously
learn to move and discriminate good/bad food.

Figure 3: a. Schematic representation of two-states Markov model with
transition probability $p_{tr}$ between two environments $E_{1}$ and $E_{2}$
defined by the ingredient value distributions. b. We vary the $E_{2}$
environment by changing the ingredient values linearly
$E_{2}=(1-2d_{e})E_{1}$, the $d_{e}$ is indicated by the color. c. The evolved
learning rate $\eta_{p}$ grows with the distance $d_{e}$ between the
environments and decreases with the reward variance $\sigma$. d. The
environment transition probability $p_{tr}$ (here for $d_{e}=1$ and
$\sigma=0.25$) has a non-monotonous relationship with the evolved learning
rate $\eta_{p}$. Up to a certain point, more rapid transitions lead to faster
learning, but too rapid environmental transition leads to a reduction of the
evolved learning rate. e. For slow environmental transition (top), the agent
fully adapts to the environment after each transition. If the transitions
happen fast (bottom), the agent maintains an intermediate position between the
two environments and never fully adapts to either of them.

### Plasticity rule parametrization

Reward-modulated plasticity is one of the most promising explanations for
biological credit assignment Legenstein et al., (2008). In our network, the
plasticity rule that updates the weights of the linear sensor network is a
reward-modulated rule which is parameterized as a linear combination of the
input, the output, and the reward at each time step:

$\Delta
W_{t}=\eta_{p}[R_{t}\cdot\overbrace{(\theta_{1}X_{t}y_{t}+\theta_{2}y_{t}+\theta_{3}X_{t}+\theta_{4})}^{\text{Reward
Modulated}}\\\
+\underbrace{(\theta_{5}X_{t}y_{t}+\theta_{6}y_{t}+\theta_{7}X_{t}+\theta_{8})}_{\text{Hebbian}}].$
(1)

Additionally, after each plasticity step, the weights are normalized by mean
subtraction, an important step for the stabilization of Hebbian-like
plasticity rules Zenke and Gerstner, (2017).

We use a genetic algorithm to optimize the learning rate $\eta_{p}$ and
amplitudes of different terms ${\theta}=(\theta_{1},\dots,\theta_{8})$. The
successful plasticity rule after many food presentations must converge to a
weight vector that predicts the correct food values (or allows the agent to
correctly decide whether to eat a food or avoid it).

To have comparable results, we divide ${\theta}=(\theta_{1},\dots,\theta_{8})$
by $\theta_{\mathrm{max}}=\max_{k}|\theta_{k}|$. So that
${\theta}/\theta_{\mathrm{max}}={\theta}^{\mathrm{norm}}\in[-1,1]^{8}$. We
then multiply the learning rate $\eta_{p}$ with $\theta_{\mathrm{max}}$ to
maintain the rule’s evolved form unchanged,
$\eta_{p}^{\mathrm{norm}}=\eta_{p}\cdot\theta_{\mathrm{max}}$. In the
following, we always use normalized $\eta_{p}$ and ${\theta}$, omitting norm.

### Evolutionary Algorithm

To evolve the plasticity rule and the moving agents’ motor networks, we use a
simple genetic algorithm with elitism Deb, (2011). The agents’ parameters are
initialized at random (drawn from a Gaussian distribution), then the sensory
network is trained by the plasticity rule and finally, the agents are
evaluated. After each generation, the best-performing agents (top 10 % of the
population size) are selected and copied into the next generation. The
remaining 90 % of the generation is repopulated with mutated copies of the
best-performing agents. We mutate agents by adding independent Gaussian noise
($\sigma=0.1$) to its parameters.

## Results

### Environmental and reward variability control the evolved learning rates of
the static agents

To start with, we consider a static agent whose goal is to identify the value
of presented food correctly. The static reward-prediction network quickly
evolves the parameters of the learning rule, successfully solving the
prediction task. We first look at the evolved learning rate $\eta_{p}$, which
determines how fast (if at all) the network’s weight vector is updated during
the lifetime of the agents. We identify three factors that control the
learning rate parameter the EA converges to: the distance between the
environments, the noisiness of the reward, and the rate of environmental
transition.

The first natural factor is the distance $d_{e}$ between the two environments,
with a larger distance requiring a higher learning rate, Fig. 3c. This is an
expected result since the convergence time to the “correct” weights is highly
dependent on the initial conditions. If an agent is born at a point very close
to optimality, which naturally happens if the environments are similar, the
distance it needs to traverse on the fitness landscape is small. Therefore it
can afford to have a small learning rate, which leads to a more stable
convergence and is not affected by noise.

A second parameter that impacts the learning rate is the variance of the
rewards. The reward an agent receives for the plasticity step contains a noise
term $\xi$ that is drawn from a zero mean Gaussian distribution with standard
deviation $\sigma$. This parameter controls the unreliability of the agent’s
sensory system, i.e., higher $\sigma$ means that the information the agent
gets about the value of the foods it consumes cannot be fully trusted to
reflect the actual value of the foods. As $\sigma$ increases, the learning
rate $\eta_{p}$ decreases, which means that the more unreliable an environment
becomes, the less an agent relies on plasticity to update its weights, Fig.
3c. Indeed for some combinations of relatively small distance $d_{e}$ and high
reward variance $\sigma$, the EA converges to a learning rate of
$\eta_{p}\approx 0$. This means that the agent opts to have no adaptation
during its lifetime and remain at the mean of the two environments. It is an
optimal solution when the expected loss due to ignoring the environmental
transitions is, on average, lower than the loss the plastic network will incur
by learning via the (often misleading because of the high $\sigma$)
environmental cues.

A final factor that affects the learning rate the EA will converge to is the
frequency of environmental change during an agent’s lifetime. Since the
environmental change is modeled as a simple, two-state Markov process (Fig.
3a), the control parameter is the transition probability $p_{tr}$.

Figure 4: The evolved parameters ${\theta}=(\theta_{1},\dots,\theta_{8})$ of
the plasticity rule for the reward prediction (a.) and the decision (b.)
tasks, for a variety of parameters ($p_{tr}=0.01$, $d_{e}\in{0,0.1,\dots,1}$,
and $\sigma\in{0,0.1,\dots,1}$ in all 100 combinations). Despite the
relatively small difference between the tasks, the evolved learning rules
differ considerably. For visual guidance, the lines connect $\theta$s from the
same run.

When keeping everything else the same, the learning rate rapidly rises as we
increase the transition probability from 0, and after reaching a peak, it
begins to decline slowly, eventually reaching zero (Fig. 3d). This means that
when environmental transition is very rare, agents opt for a very low learning
rate, allowing a slow and stable convergence to an environment-appropriate
weight vector that leads to very low losses while the agent remains in that
environment. As the rate of environmental transition increases, faster
learning is required to speed up convergence in order to exploit the
(comparatively shorter) stays in each environment. Finally, as the
environmental transition becomes too fast, the agents opt for slower or even
no learning, which keeps them near the middle of the two environments,
ensuring that the average loss of the two environments is minimal (Fig. 3d).

### The form of the evolved learning rule depends on the task: Decision vs.
Prediction

The plasticity parameters ${\theta}=(\theta_{1},\dots,\theta_{8})$ for the
reward-prediction task converge on approximately the same point, regardless of
the environmental parameters (Fig. 4a). In particular, $\theta_{3}\to 1$,
$\theta_{5}\to-1$, $\theta_{i}\to 0$ for all other $i$, and thus the learning
rule converges to:

$\Delta
W_{t}=\eta_{p}[\theta_{3}X_{t}R_{t}+\theta_{5}X_{t}y_{t}]\approx\eta_{p}X_{t}(R_{t}-y_{t}).$
(2)

Since by definition $y_{t}=g(W_{t}X_{t}^{T})=W_{t}X_{t}^{T}$ ($g(x)=x$ in this
experiment) and $R_{t}=W^{c}X_{t}^{T}+\xi$ we get:

$\Delta W_{t}=\eta_{p}X_{t}(W^{c}-W_{t})X_{t}^{T}+\eta_{p}\xi X_{t}^{T}.$ (3)

Thus the distribution of $\Delta W_{t}$ converges to a distribution with mean
0 and variance depending on $\eta_{p}$ and $\sigma$ and $W$ converges to
$W^{c}$. So this learning rule will match the agent’s weight vector with the
vector of ingredient values in the environment.

We examine the robustness of the learning rule the EA discovers by considering
a slight modification of our task. Instead of predicting the expected food
value, the agent now needs to decide whether to eat the presented food or not.
This is done by introducing a step-function nonlinearity ($g(x)=1$ if $x\geq
1$ and 0 otherwise). Then the output $y(t)$ is computed as:

$y_{t}=\begin{cases}1,&\text{if }\ W_{t}X_{t}^{T}\geq 0,\\\ 0,&\text{if }\
W_{t}X_{t}^{T}<0.\\\ \end{cases}$ (4)

Instead of the MSE loss between prediction and actual value, the fitness of
the agent is now defined as the sum of the food values it chose to consume (by
giving $y_{t}=1$). Besides these two changes, the setup of the experiments
remains exactly the same.

The qualitative relation between $\eta_{p}$ and parameters of environment
$d_{e},\sigma$ and $p_{tr}$ is preserved in the changed experiment. However,
the resulting learning rule is significantly different (Fig. 4). The evolution
converges to the following learning rule:

$\Delta W_{t}=\begin{cases}\eta_{p}X_{t}[\theta_{3}R_{t}+\theta_{7}],\
y_{t}=0,\\\
\eta_{p}X_{t}[(\theta_{1}+\theta_{3})R_{t}+(\theta_{5}+\theta_{7})],\
y_{t}=1.\\\ \end{cases}$ (5)

In both cases, the rule has the form $\Delta
W_{t}=\eta_{p}X_{t}[\alpha_{y}R_{t}+\beta_{y}]$. Thus, the $\Delta W_{t}$ is
positive or negative depending on whether the reward $R_{t}$ is above or below
a threshold ($\gamma=-\beta_{y}/\alpha_{y}$) that depends on the output
decision of the network ($y_{t}=0$ or $1$).

Both learning rules (for the reward-prediction and decision tasks) have a
clear Hebbian form (coordination of pre- and post-synaptic activity) and use
the incoming reward signal as a threshold. These similarities indicate some
common organizing principles of reward-modulated learning rules, but their
significant differences highlight the sensitivity of the optimization process
to task details.

### The learning rate of embodied agents depends on environmental variability

Figure 5: a. The trajectory of an agent (blue line) in the 2D environment. A
well-trained agent will approach and consume food with positive values (green
dots) and avoid negative food (red dots). b. The learning rate of the plastic
sensory network $eta_{p}$ grows with the distance between environments $d_{e}$
c. and decreases with the frequency of environmental change. d. The fitness of
an agent (measured as the total food consumed over its lifetime) increases
over generations of the EA for both the scalar and binary readouts in the
sensory network. e. The Pearson correlation coefficient of an evolved agent’s
weights with the ingredient value vector of the current environment ($E_{1}$
\- blue, $E_{2}$ \- red). In this example, the agent’s weights are anti-
correlated with its environment, which is not an issue for performance since
the motor network can interpret the inverted signs of food.

We now turn to the moving embodied agents in the 2D environment. To optimize
these agents, both the motor network’s connections and the sensory network’s
plasticity parameters evolve simultaneously. Since the motor network is
initially random and the agent has to move to find food, the number of
interactions an agent experiences in its lifetime can be small, slowing down
the learning. However, having the larger motor network also has benefits for
evolution because it allows the output of the plastic network to be read out
and transformed in different ways, resulting in a broad set of solutions.

The agents can solve the task effectively by evolving a functional motor
network and a plasticity rule that converges to interpretable weights (Fig.
5a). After $\sim$ 100 evolutionary steps (Fig. 5d), the agents can learn the
ingredient value distribution using the plastic network and reliably move
towards foods with positive values while avoiding the ones with negative
values.

We compare the dependence of the moving and the static agents on the
parameters of the environment: $d_{e}$ and the state transition probability
$p_{tr}$. At first, in order to simplify the experiment, we set the transition
probability to $0$, but fixed the initial weights to be the average of $E_{1}$
and $E_{2}$, while the real state is $E_{2}$. In this experiment, the distance
between states $d_{e}$ indicates twice the distance between the agent’s
initial weights and the optimal weights (the environment’s ingredient values)
since the agent is initialized at the mean of the two environment
distributions. Same as for the static agent, the learning rate increases with
the distance $d_{e}$ (Fig. 5b).

Then, we examine the effect of the environmental transition probability
$p_{tr}$ on the evolved learning rate $\eta_{p}$. In order for an agent to get
sufficient exposure to each environment, we scale down the probability
$p_{tr}$ from the equivalent experiment for the static agents. We find that as
the probability of transition increases, the evolved learning rate $\eta_{p}$
decreases (Fig. 5c). This fits with the larger trend for the static agent,
although there is a clear difference when it comes to the increase for very
small transition probabilities that were clearly identifiable in the static
but not the moving agents. This could be due to much sparser data and possibly
the insufficiently long lifetime of the moving agent (the necessity of scaling
makes direct comparisons difficult). Nevertheless, overall we see that the
associations observed in the static agents between environmental distance
$d_{e}$ and transition probability $p_{tr}$ and the evolved learning rate
$\eta_{p}$ are largely maintained in the moving agents. Still, more data would
be needed to make any conclusive assertions about the exact effect of these
environmental parameters on the emerging plasticity mechanisms.

### Rule redundancy in the embodied agents

Figure 6: The evolved parameters of moving agents’ plasticity rule for the
$g(s)=x$, identity (a.) and the step function (Eq. 4) (b.) sensory networks
(the environmental parameters here are $d_{e}\in[0,1],\ \sigma=0$ and
$p_{tr}=0.001$). The step function (binary output) network evolved a more
structured plasticity rule (e.g., $\theta_{3}>0$ for all realizations) than
the linear network. Moreover, the learned weights for the identity network
(c.) have higher variance and correlate significantly less with the
environment’s ingredient distribution compared to the learned weights for the
thresholded network (d.)

A crucial difference between the static and the moving agents is the function
the plasticity has to perform. While in the static agents, the plasticity has
to effectively identify the exact value distribution of the environment in
order to produce accurate predictions, in the embodied agents, the plasticity
has to merely produce a representation of the environment that the motor
network can evolve to interpret adequately enough to make decisions about
which food to consume.

To illustrate the difference, we plot the Pearson correlation coefficient
between an agent’s weights and the ingredient values of the environment it is
moving in (Fig. 5e). We use the correlation instead of the MSE loss (which we
used for the static agents in Fig. 3e) because the amplitude of the weight
vector varies a lot for different agents and meaningful conclusions cannot be
drawn from the MSE loss. For many agents, the learned weights are consistently
anti-correlated with the actual ingredient values (an example of such an agent
is shown in Fig. 5e). This means that the output of the sensory network will
have the opposite sign from the actual food value. While in the static
network, this would lead to very bad predictions and high loss, in the
foraging task, these agents perform exactly as well as the ones where the
weights and ingredients values are positively correlated, since the motor
network can simply learn to move towards food for which it gets a negative
instead of a positive sensory input.

This additional step of the output of the plastic network going through the
motor network before producing any behavior has a strong effect on the
plasticity rules that the embodied agents evolve. Specifically, if we look at
the emerging rules the top performing agents have evolved (Fig. 6a), it
becomes clear that, unlike the very well-structured rules of the static agents
(Fig. 4a), there is now virtually no discernible pattern or structure. The
difference becomes even clearer if we look at the learned weights (at the end
of a simulation) of the best-performing agents (Fig. 6c). While there is some
correlation with the environment’s ingredient value distribution, the variance
is very large, and they do not seem to converge on the “correct” values in any
way. This is to some extent expected since, unlike the static agents where the
network’s output has to be exactly correct, driving the evolution of rules
that converge to the precise environmental distribution, in the embodied
networks, the bulk of the processing is done by the motor network which can
evolve to interpret the scalar value of the sensory network’s output in a
variety of ways. Thus, as long as the sensory network’s plasticity rule co-
evolves with the motor network, any plasticity rule that learns to produce
consistent information about the value of encountered food can potentially be
selected.

To further test this assumption, we introduce a bottleneck of information
propagation between the sensory and motor networks by using a step-function
nonlinearity on the output of the sensory network (Eq. 4). Similarly to the
decision task of the static network, the output of the sensory network now
becomes binary. This effectively reduces the flow of information from the
sensory to the motor network, forcing the sensory network to consistently
decide whether food should be consumed (with the caveat that the motor network
can still interpret the binary sign in either of two ways, either consuming
food marked with $1$ or the ones marked with $0$ by the sensory network). The
agents perform equally well in this variation of the task as before (Fig. 5d),
but now, the evolved plasticity rules seem to be more structured (Fig. 6b).
Moreover, the variance of the learned weights in the best-performing agents is
significantly reduced (Fig. 6d), which indicates that the bottleneck in the
sensory network is increasing selection pressure for rules that learn the
environment’s food distribution accurately.

## Discussion

We find that different sources of variability have a strong impact on the
extent to which evolving agents will develop neuronal plasticity mechanisms
for adapting to their environment. A diverse environment, a reliable sensory
system, and a rate of environmental change that is neither too large nor too
small are necessary conditions for an agent to be able to effectively adapt
via synaptic plasticity. Additionally, we find that minor variations of the
task an agent has to solve or the parametrization of the network can give rise
to significantly different plasticity rules.

Our results partially extend to embodied artificial agents performing a
foraging task. We show that environmental variability also pushes the
development of plasticity in such agents. Still, in contrast to the static
agents, we find that the interaction of a static motor network with a plastic
sensory network gives rise to a much greater variety of well-functioning
learning rules. We propose a potential cause of this degeneracy; as the
relatively complex motor network is allowed to read out and process the
outputs from the plastic network, any consistent information coming out of
these outputs can be potentially interpreted in a behaviorally useful way.
Reducing the information the motor network can extract from the sensory system
significantly limits learning rule variability.

Our findings on the effect of environmental variability concur with the
findings of previous studies Lange and Sprekeler, (2020) that have identified
the constraints that environmental variability places on the evolutionary
viability of learning behaviors. We extend these findings in a mechanistic
model which uses a biologically plausible learning mechanism (synaptic
plasticity). We show how a simple evolutionary algorithm can optimize the
different parameters of a simple reward-modulated plasticity rule for solving
simple prediction and decision tasks. Reward-modulated plasticity has been
extensively studied as a plausible mechanism for credit assignment in the
brain Florian, (2007); Baras and Meir, (2007); Legenstein et al., (2008) and
has found several applications in artificial intelligence and robotics tasks
Burms et al., (2015); Bing et al., (2019). Here, we demonstrate how such rules
can be very well-tuned to take into account different environmental parameters
and produce optimal behavior in simple systems.

Additionally, we demonstrate how the co-evolution of plasticity and static
functional connectivity in different sub-networks fundamentally changes the
evolutionary pressures on the resulting plasticity rules, allowing for greater
diversity in the form of the learning rule and the resulting learned
connectivity. Several studies have demonstrated how, in biological networks,
synaptic plasticity heavily interacts with Butz et al., (2014); Stampanoni
Bassi et al., (2019); Bernáez Timón et al., (2022) and is driven by network
topology Giannakakis et al., (2023). Moreover, it has been recently
demonstrated that biological plasticity mechanisms are highly redundant in the
sense that any observed neural connectivity or recorded activity can be
achieved with a variety of distinct, unrelated learning rules Ramesh, (2023).
This observed redundancy of learning rules in biological settings complements
our results and suggests that the function of plasticity rules cannot be
studied independently of the connectivity and topology of the networks they
are acting on.

The optimization of functional plasticity in neural networks is a promising
research direction both as a means to understand biological learning processes
and as a tool for building more autonomous artificial systems. Our results
suggest that reward-modulated plasticity is highly adaptable to different
environments and can be incorporated into larger systems that solve complex
tasks.

## Future work

This work studies a simplified toy model of neural network learning in
stochastic environments. Future work could be built on this basic framework to
examine more complex reward distributions and sources of environmental
variability. Moreover, a greater degree of biological realism could be added
by studying more plausible network architectures (multiple plastic layers,
recurrent and feedback connections) and more sophisticated plasticity rule
parametrizations.

Additionally, our foraging simulations were constrained by limited
computational resources and were far from exhaustive. Further experiments can
investigate environments with different constraints, food distributions,
multiple seasons, more complex motor control systems and interactions of those
systems with different sensory networks as well as the inclusion of plasticity
on the motor parts of the artificial organisms.

## Acknowledgements

This work was supported by a Sofja Kovalevskaja Award from the Alexander von
Humboldt Foundation. EG and SK thank the International Max Planck Research
School for Intelligent Systems (IMPRS-IS) for their support. We acknowledge
the support from the BMBF through the Tübingen AI Center (FKZ: 01IS18039A). AL
is a member of the Machine Learning Cluster of Excellence, EXC number 2064/1 –
Project number 39072764.

## References

  * Baras and Meir, (2007) Baras, D. and Meir, R. (2007).  Reinforcement Learning, Spike-Time-Dependent Plasticity, and the BCM Rule.  Neural Computation, 19(8):2245–2279.
  * Bernáez Timón et al., (2022) Bernáez Timón, L., Ekelmans, P., Konrad, S., Nold, A., and Tchumatchenko, T. (2022).  Synaptic plasticity controls the emergence of population-wide invariant representations in balanced network models.  Phys. Rev. Res., 4:013162.
  * Biesialska et al., (2020) Biesialska, M., Biesialska, K., and Costa-jussà , M. R. (2020).  Continual lifelong learning in natural language processing: A survey.  In Proceedings of the 28th International Conference on Computational Linguistics. International Committee on Computational Linguistics.
  * Bing et al., (2019) Bing, Z., Baumann, I., Jiang, Z., Huang, K., Cai, C., and Knoll, A. (2019).  Supervised learning in snn via reward-modulated spike-timing-dependent plasticity for a target reaching vehicle.  Frontiers in Neurorobotics, 13.
  * Burms et al., (2015) Burms, J., Caluwaerts, K., and Dambre, J. (2015).  Reward-modulated hebbian plasticity as leverage for partially embodied control in compliant robotics.  Frontiers in Neurorobotics, 9.
  * Butz et al., (2014) Butz, M., Steenbuck, I., and van Ooyen, A. (2014).  Homeostatic structural plasticity increases the efficiency of small-world networks.  Frontiers in Synaptic Neuroscience, 6.
  * Caroni et al., (2012) Caroni, P., Donato, F., and Muller, D. (2012).  Structural plasticity upon learning: regulation and functions.  Nature reviews. Neuroscience, 13(7):478—490.
  * Citri and Malenka, (2008) Citri, A. and Malenka, R. (2008).  Synaptic plasticity: Multiple forms, functions, and mechanisms.  Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology, 33:18–41.
  * Confavreux et al., (2020) Confavreux, B., Zenke, F., Agnes, E., Lillicrap, T., and Vogels, T. (2020).  A meta-learning approach to (re)discover plasticity rules that carve a desired function into a neural network.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 16398–16408. Curran Associates, Inc.
  * Deb, (2011) Deb, K. (2011).  Multi-objective Optimisation Using Evolutionary Algorithms: An Introduction, pages 3–34.  Springer London, London.
  * Dunlap and Stephens, (2009) Dunlap, A. S. and Stephens, D. W. (2009).  Components of change in the evolution of learning and unlearned preference.  Proceedings of the Royal Society B: Biological Sciences, 276(1670):3201–3208.
  * Dunlap and Stephens, (2016) Dunlap, A. S. and Stephens, D. W. (2016).  Reliability, uncertainty, and costs in the evolution of animal learning.  Current Opinion in Behavioral Sciences, 12:73–79.  Behavioral ecology.
  * Ellefsen, (2014) Ellefsen, K. O. (2014).  The evolution of learning under environmental variability.  pages 649–656.
  * Eskridge and Hougen, (2012) Eskridge, B. E. and Hougen, D. F. (2012).  Nurturing promotes the evolution of learning in uncertain environments.  In 2012 IEEE International Conference on Development and Learning and Epigenetic Robotics (ICDL), pages 1–6.
  * Fawcett et al., (2012) Fawcett, T. W., Hamblin, S., and Giraldeau, L.-A. (2012).  Exposing the behavioral gambit: the evolution of learning and decision rules.  Behavioral Ecology, 24(1):2–11.
  * Feldman, (2009) Feldman, D. E. (2009).  Synaptic mechanisms for plasticity in neocortex.  Annual Review of Neuroscience, 32(1):33–55.  PMID: 19400721.
  * Florian, (2007) Florian, R. V. (2007).  Reinforcement Learning Through Modulation of Spike-Timing-Dependent Synaptic Plasticity.  Neural Computation, 19(6):1468–1502.
  * Giannakakis et al., (2023) Giannakakis, E., Vinogradov, O., Buendia, V., and Levina, A. (2023).  Recurrent connectivity structure controls the emergence of co-tuned excitation and inhibition.  bioRxiv.
  * Guttenberg, (2019) Guttenberg, N. (2019).  Evolutionary rates of information gain and decay in fluctuating environments.  ALIFE 2019: The 2019 Conference on Artificial Life:365–371.
  * Jordan et al., (2021) Jordan, J., Schmidt, M., Senn, W., and Petrovici, M. A. (2021).  Evolving interpretable plasticity for spiking networks.  eLife, 10:e66273.
  * Kerr and Feldman, (2003) Kerr, B. and Feldman, M. (2003).  Carving the cognitive niche: Optimal learning strategies in homogeneous and heterogeneous environments.  Journal of Theoretical Biology, 220(2):169–188.
  * Khajehabdollahi et al., (2022) Khajehabdollahi, S., Prosi, J., Giannakakis, E., Martius, G., and Levina, A. (2022).  When to Be Critical? Performance and Evolvability in Different Regimes of Neural Ising Agents.  Artificial Life, 28(4):458–478.
  * Lange and Sprekeler, (2020) Lange, R. T. and Sprekeler, H. (2020).  Learning not to learn: Nature versus nurture in silico.
  * Lee and Lee, (2020) Lee, C. and Lee, A. (2020).  Clinical applications of continual learning machine learning.  The Lancet Digital Health, 2:e279–e281.
  * Legenstein et al., (2008) Legenstein, R., Pecevski, D., and Maass, W. (2008).  A learning theory for reward-modulated spike-timing-dependent plasticity with application to biofeedback.  PLOS Computational Biology, 4(10):1–27.
  * Magee and Grienberger, (2020) Magee, J. C. and Grienberger, C. (2020).  Synaptic plasticity forms and functions.  Annual Review of Neuroscience, 43(1):95–117.  PMID: 32075520.
  * Najarro and Risi, (2020) Najarro, E. and Risi, S. (2020).  Meta-learning through hebbian plasticity in random networks.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 20719–20731. Curran Associates, Inc.
  * Nolfi and Parisi, (1996) Nolfi, S. and Parisi, D. (1996).  Learning to adapt to changing environments in evolving neural networks.  Adaptive Behavior, 5(1):75–98.
  * Open Ended Learning Team et al., (2021) Open Ended Learning Team, Stooke, A., Mahajan, A., Barros, C., Deck, C., Bauer, J., Sygnowski, J., Trebacz, M., Jaderberg, M., Mathieu, M., McAleese, N., Bradley-Schmieg, N., Wong, N., Porcel, N., Raileanu, R., Hughes-Fitt, S., Dalibard, V., and Czarnecki, W. M. (2021).  Open-ended learning leads to generally capable agents.
  * Papini, (2012) Papini, M. R. (2012).  Evolution of Learning, pages 1188–1192.  Springer US, Boston, MA.
  * Parisi et al., (2019) Parisi, G. I., Kemker, R., Part, J. L., Kanan, C., and Wermter, S. (2019).  Continual lifelong learning with neural networks: A review.
  * Pedersen and Risi, (2021) Pedersen, J. W. and Risi, S. (2021).  Evolving and merging hebbian learning rules: Increasing generalization by decreasing the number of rules.  In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’21, page 892–900, New York, NY, USA. Association for Computing Machinery.
  * Ramesh, (2023) Ramesh, P. (2023).  GANs schön kompliziert: Applications of Generative Adversarial Networks.  PhD thesis, University of Tübingen.
  * Schmidhuber, (1987) Schmidhuber, J. (1987).  Evolutionary principles in self-referential learning. on learning now to learn: The meta-meta-meta…-hook.  Diploma thesis, Technische Universitat Munchen, Germany.
  * Snell-Rood, (2013) Snell-Rood, E. C. (2013).  An overview of the evolutionary causes and consequences of behavioural plasticity.  Animal Behaviour, 85(5):1004–1011.  Special Issue: Behavioural Plasticity and Evolution.
  * Snell-Rood and Steck, (2019) Snell-Rood, E. C. and Steck, M. K. (2019).  Behaviour shapes environmental variation and selection on learning and plasticity: review of mechanisms and implications.  Animal Behaviour, 147:147–156.
  * Stampanoni Bassi et al., (2019) Stampanoni Bassi, M., Iezzi, E., Gilio, L., Centonze, D., and Buttari, F. (2019).  Synaptic plasticity shapes brain connectivity: Implications for network topology.  International Journal of Molecular Sciences, 20(24).
  * Thangarasa et al., (2020) Thangarasa, V., Miconi, T., and Taylor, G. W. (2020).  Enabling continual learning with differentiable hebbian plasticity.
  * Thornton and Clutton-Brock, (2011) Thornton, A. and Clutton-Brock, T. (2011).  Social learning and the development of individual and group behaviour in mammal societies.  Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 366:978–87.
  * Yaman et al., (2021) Yaman, A., Iacca, G., Mocanu, D. C., Coler, M., Fletcher, G., and Pechenizkiy, M. (2021).  Evolving Plasticity for Autonomous Learning under Changing Environmental Conditions.  Evolutionary Computation, 29(3):391–414.
  * Zenke and Gerstner, (2017) Zenke, F. and Gerstner, W. (2017).  Hebbian plasticity requires compensatory processes on multiple timescales.  Philosophical Transactions of the Royal Society B: Biological Sciences, 372(1715):20160259.