# Scenario Based Cost Optimization of Water Distribution Networks Powered by
Grid-Connected Photovoltaic Systems

Mirhan Ürkmez  Aalborg University, Fredrik Bajers Vej 7c, DK-9220 Aalborg,
Denmark
(e-mail<EMAIL_ADDRESS>Carsten Kallesøe  Aalborg University,
Fredrik Bajers Vej 7c, DK-9220 Aalborg, Denmark
(e-mail<EMAIL_ADDRESS>Jan Dimon Bendtsen  Aalborg University,
Fredrik Bajers Vej 7c, DK-9220 Aalborg, Denmark
(e-mail<EMAIL_ADDRESS>John Leth  Aalborg University, Fredrik
Bajers Vej 7c, DK-9220 Aalborg, Denmark
(e-mail<EMAIL_ADDRESS>
###### Abstract

The paper presents a predictive control method for the water distribution
networks (WDNs) powered by photovoltaics (PVs) and the electrical grid. This
builds on the controller introduced in a previous study and is designed to
reduce the economic costs associated with operating the WDN. To account for
the uncertainty of the system, the problem is solved in a scenario
optimization framework, where multiple scenarios are sampled from the
uncertain variables related to PV power production. To accomplish this, a day-
ahead PV power prediction method with a stochastic model is employed. The
method is tested on a high-fidelity model of a WDN of a Danish town and the
results demonstrate a substantial reduction in electrical costs through the
integration of PVs, with PVs supplying $66.95\%$ of the required energy. The
study also compares the effectiveness of the stochastic optimization method
with a deterministic optimization approach.

##  I INTRODUCTION

Water distribution networks (WDNs) transport potable water from its source to
the end users. Approximately $7\%-8\%$ of the world’s energy is utilized for
water production and distribution [1]. Since many countries are trying to
increase installed renewable energy sources such as Photovoltaic (PV) and
wind, it is natural to ask whether WDNs can run efficiently with renewable
energy sources. Specifically, this paper focuses on controlling water pumps in
WDNs powered by grid-connected PV panels to increase the penetration of
renewable energy and reduce energy costs. Pump scheduling is a complex problem
because of the nonlinearities governing the network elements and the large
size of the networks. The problem is exacerbated by the introduction of PV
panels due to the uncertainties in PV power production. Another way to
approach the integration of photovoltaics (PVs) and WDNs is to consider it as
the design of an Energy Management System (EMS) for a microgrid. Microgrids
consist of energy sources, storage devices, and interconnected loads. In this
case, the grid can be thought of as having WDNs as its loads and water tanks
serving as its storage units.

There has been a significant amount of research focused on optimizing the
scheduling of pumps in WDNs to reduce costs. Due to the size of the network
and nonlinear pipe equations, it is common to approximate the network with a
simpler and sometimes linear model. In [2], an Economic Model Predictive
Control (EMPC) is employed with a Linear Parameter Varying (LPV) system model.
In [3], the network structure is simplified by dividing nodes into clusters
and representing each cluster with a single node. Then a system model is
developed from the simplified structure using a Deep Neural Network (DNN)
structure. In some works, pipe equations are replaced with linear equations or
inequalities to simplify the original problem. This approach is taken in [4],
where the pipe equations are linearized around an operating point and Model
Predictive Control (MPC) is applied, and in [5], where EMPC is constructed
using a set of linear inequalities derived from relaxing the nonlinear pipe
equations. In our study, we assume the presence of an elevated reservoir in
the network, which allows us to consider the node pressures to be within an
acceptable range. Consequently, we do not need to solve the entire network to
verify node pressure constraints. Instead, we adopt a simpler model to
represent changes in tank levels. Similarly, a network structure with an
elevated reservoir is considered and data is utilized to identify a reduced
system model in [6]. As there is an elevated reservoir, pressure constraints
are not included in problem formulation.

One of the main challenges in integrating photovoltaic (PV) panels with water
distribution networks (WDNs) is accurately predicting the power production of
the PV panels. To address this challenge, several different estimation models
have been developed that use weather measurements, historical power production
data, and numerical weather predictions (NWP). NWP has been demonstrated to be
effective for making predictions as it can take into account factors such as
rain and clouds that can significantly impact PV performance [7]. NWP data is
used to classify days into different categories such as sunny, cloudy, or
rainy, which allows for more accurate power production estimates to be made
[8, 9, 10]. Some models use physical system properties of PVs to formulate
power as a function of weather variables such as radiance, temperature, and
wind speed, and then use NWP data to make predictions [11, 12, 13]. Certain
studies have refrained from incorporating NWP due to its limitations. In [14],
the use of NWP is avoided due to concerns about its availability across all
locations. Similarly, in our PV forecasting method,we have chosen not to rely
on NWP data due to its limited effectiveness for short-term predictions within
a time frame of up to 6 hours [8]. Neural Networks (NN) with different
prediction horizons [15] and multiple inputs types (e.g. NWP, historical
production data) [8, 10] have also been used in the literature. Lately, long
short-term memory (LSTM) networks are being applied because of their success
in time series forecasting [16, 17, 9]. A limitation of most NN-based methods
is their inability to quantify prediction uncertainty rendering them
unsuitable for integration within a stochastic control framework.

The EMS of microgrids with PV panels connected to main grid has been studied
with different applications. In [18], an MPC formulated with mixed-integer
linear programming (MILP) is employed as EMS for multiple residential
microgrids powered with PVs. In [19], a distributed cooperative approach is
used for a network of microgrids with an aggregator determining the
distribution of power and individual microgrids able to deviate slightly from
the power profile determined by the aggregator. Since PV power production is
highly volatile, some works have used scenario-based optimization [20, 21].
Robust satisfaction of the constraints has also been considered alongside
scenario-optimization in [22].

The main contribution of this paper is a method for the control of the pumps
in the WDNs that are powered by grid-connected PV systems. PV power production
is modeled with a probabilistic model to account for the uncertainties in the
production. A stochastic predictive controller method based on a linear model
of the system is used to determine the pump flows. Scenario optimization is
used to solve the controller problem by utilizing the probabilistic PV model
to sample the uncertain variables. The method is applied to the EPANET, a
hydraulic simulation software, model of a medium-sized Danish town’s network
(Randers).

The outline of the rest of the paper is as follows. PV power estimation method
is given in Section II. The model of the network is derived in Section III.
The control method is explained in Section IV. The experimental results are
presented in Section V. The paper is concluded in Section VI.

##  II PV Power Forecasting

This section presents a method for predicting the PV power output, which was
first introduced in [23]. The method is designed to make day-ahead predictions
based on historical power data. At night, the prediction is done in two steps.
Firstly, the normalized shape of power output for the following day is
estimated, after which a multiplier value for this shape is determined. This
multiplier value is then applied to the normalized shape to give the final
prediction.

The method used for predicting the shape of the power production data is the
Exponentially Weighted Moving Average (EWMA). Let $\eta\in\mathbb{N}$ denote
the current day. At midnight between day $\eta-1$ and $\eta$, the production
shape of the current day given the historical data is estimated as

$\displaystyle X^{\prime}_{\eta-1}$
$\displaystyle=\frac{1}{\max(X_{\eta-1})}X_{\eta-1},$ (1a) $\displaystyle
Y_{\eta}$ $\displaystyle=\alpha X^{\prime}_{\eta-1}+(1-\alpha)Y_{\eta-1}$ (1b)

where $X_{\eta-1}\in\mathbb{R}^{N_{pv}}$ is the vector of production data for
the day $\eta-1$, $N_{pv}$ is the number of power data points in a day given a
sampling time of $\Delta_{pv}$, $X^{\prime}_{\eta-1}\in[0,1]^{N_{pv}}$ is the
normalized daily data, $\alpha\in[0,1]$, and $Y_{\eta}\in[0,1]^{N_{pv}}$ is
the moving average normalized daily production of day $\eta$. The parameter
$\alpha$ controls the relative weight given to more recent days versus older
days.

The next step is to determine the coefficient that will be used to multiply
the estimated normalized shape $Y_{\eta}$ to produce the final prediction. To
do this, the optimal multiplier value of the previous days is calculated using
the optimization problem

$\displaystyle
p_{\tau}=\operatorname*{arg\,min}_{p}\quad\sum_{i=1}^{N}{(pY_{\tau}^{i}-X_{\tau}^{i})^{2}}$
(2)

where $p\in\mathbb{R}$ is the multiplier, and $Y_{\tau}^{i}$ (resp.
$X_{\tau}^{i}$) are the $i^{th}$ coordinate of $Y_{\tau}$ (resp. $X_{\tau}$).
Since the optimal multiplier value $p_{\eta}$ for day $\eta$ can only be
calculated when all power production values $X_{\eta}$ of the day become
available, that is at midnight between day $\eta$ and $\eta+1$, the strategy
during the day $\eta$ is to estimate the optimal multiplier $p_{\eta}$.

To estimate $p_{\eta}$, the time series of optimal multipliers
$p_{\tau},\tau=1,\dots,\eta-1$ are modelled using an Autoregressive Moving
Average model ARMA(1,1) given by

$\displaystyle p_{\tau}=\mu+\phi
p_{\tau-1}+\theta\epsilon_{\tau-1}+\epsilon_{\tau}$ (3)

where $\epsilon_{\tau}$ and $\epsilon_{\tau-1}$ are the forecast errors at
step $\tau,\tau-1$ respectively and $\mu,\phi,\theta$ are the constants whose
values are set to minimize the sum of squared errors, $\epsilon_{\tau}$. The
error terms $\epsilon_{\tau}=p_{\tau}-\hat{p}_{\tau},~{}\tau=1,\dots,\eta$ are
assumed to be coming from independent zero-mean normal distributions. The
estimation of the next value $\hat{p}_{\eta}$ is made using the calculated
optimal value of the previous day $p_{\eta-1}$ as

$\displaystyle\hat{p}_{\eta}=\mu+\phi p_{\eta-1}+\theta\epsilon_{\eta-1}.$ (4)

Then, the prediction for the day $\eta$ is made as $\hat{p}_{\eta}Y_{\eta}$ at
midnight.

The problem of predicting during daytime is handled differently than
predicting during nighttime since the production values available during the
day can be used to predict future power values in the day. Firstly, sunrise is
determined when consecutive power values surpass a predetermined limit. Then,
an estimation of the optimal multiplier $p_{\eta}$ is made at each time step
using the power data after the sunrise by solving the problem given by

$\displaystyle\hat{p}^{m}_{\eta}=\arg\min_{p}\quad\sum_{i=s_{r}}^{s_{r}+m}{(pY_{\eta}^{i}-X_{\eta}^{i})^{2}}$
(5)

where $\hat{p}^{m}_{\eta}$ is the estimated optimal multiplier value $m$ time
steps after the sunrise time $s_{r}$ based on the available current day data.

We now have a prior estimation $\hat{p}_{\eta}$ and an observation
$\hat{p}^{m}_{\eta}$ of the optimal multiplier $p_{\eta}$. The prior is
assumed to be normally distributed as
$p_{\eta}=\hat{p}_{\eta}+\epsilon_{\eta}$ where $\epsilon_{\eta}\sim
N(0,\sigma^{2})$. Let $z^{m}_{\eta}$ denote the random variable associated
with the observations given by (5). We assume that $z^{m}_{\eta}$ and
$p_{\eta}$ are related as $z^{m}_{\eta}=p_{\eta}-\epsilon^{m}_{\eta}$ where
$\epsilon^{m}_{\eta}\sim N(0,\sigma_{m}^{2})$. Then, likelihood
$\mathcal{L}(p_{\eta}\mid\hat{p}^{m}_{\eta})=P(z^{m}_{\eta}=\hat{p}^{m}_{\eta}\mid
p_{\eta})$ can also be represented with a normal distribution. The variances
$\sigma^{2}$ and $\sigma_{m}^{2}$ are estimated from the sets of available
$\epsilon_{\tau}$ and $\epsilon^{m}_{\tau},\tau=1,\dots,\eta-1$ values from
the previous days. The observation $z^{m}_{\eta}=\hat{p}^{m}_{\eta}$ and the
prior estimation $\hat{p}_{\eta}$ are combined in a Bayesian setting for a
better estimation of $p_{\eta}$ as

$\displaystyle P(p_{\eta}=x\mid z^{m}_{\eta}=\hat{p}^{m}_{\eta})\propto
P(p_{\eta}=x)P(z^{m}_{\eta}=\hat{p}^{m}_{\eta}\mid p_{\eta}=x)$ (6a)
$\displaystyle\hat{p}_{\eta\mid m}=\arg\max_{x}P(p_{\eta}=x\mid
z^{m}_{\eta}=\hat{p}^{m}_{\eta})$ (6b)

The posterior probability distribution
$P(p_{\eta}|z^{m}_{\eta}=\hat{p}^{m}_{\eta})$ is also normally distributed
since both the likelihood and prior are normally distributed. The optimal
multiplier estimation $\hat{p}_{\eta|m}$ is selected as the value which
maximizes the posterior probability distribution, as expressed in (6b).

Differing from the original paper, we also model the error between the optimal
estimation $p_{\eta}Y_{\eta}$ and the actual values $X_{\eta}$ to incorporate
it into the calculation of the expected cost value of the controller
optimization problem outlined in Section IV. Let
$\delta^{i}_{\tau}=X_{\tau}^{i}-p_{\tau}Y_{\tau}^{i}$ denote the error at time
$i$ in day $\tau$. The daytime error values $\delta^{i}_{\tau}$ are assumed to
accept an ARMA(1,0) model given by

$\displaystyle\delta^{i}_{\tau}=\phi^{\delta^{i}_{\tau}}\delta^{i-1}_{\tau}+\epsilon^{\delta^{i}_{\tau}}$
(7)

where $\epsilon^{\delta^{i}_{\tau}}$ are independent and identically
distributed with $N(0,\sigma^{\delta^{i}_{\tau}})$. It is also assumed that
$\epsilon^{\delta^{i}_{\tau}}$ and $p_{\tau}$ are independent when the value
of $X_{\tau}^{i}$ is not known. Since there is no power production at night,
nighttime errors are not modeled. The model is fitted to the time series
$\delta^{i}_{\tau},\tau=1,\cdots\eta-1$ and the unknown parameters
$\sigma^{\delta^{i}_{\tau}},\phi^{\delta^{i}_{\tau}}$ have been determined.

##  III Network Model

A water distribution network typically consists of pipes, pumps, tanks,
junction nodes and reservoirs. Water networks are commonly divided into zones
sharing similar properties like altitude or water consumption. Water flows
through the network, driven by the difference in hydraulic head, which is a
measure of the fluid pressure and is equivalent to the height of a fluid held
in a static column at a given point. The only dynamic elements in a network
are the tanks whose water level change as

$A_{j}\dot{{h}}_{j}=\sum_{i\in\mathpzc{N}_{j}}q_{ij}$ (8)

where $A_{j}$ is the cross-sectional area, $h_{j}$ is the level of the tank,
$q_{ij}$ is the flow entering the tank $j$, $\mathpzc{N}_{j}$ denotes the set
of neighbor nodes of the tank $j$. Flow through pipes $q_{ij}$ connected to
the tanks are nonlinear functions $f_{i}$ of the demand at each node, tank
levels, and the amount of water coming from the pumps. To obtain explicit
forms of these functions, demand data for every node
$d=[d_{1}~{}d_{2}\cdots]^{T}$ must be available; however, this is often not
the case. We assume that the total demand of the zones supplied by the pumps
can be estimated through available data via time series analysis methods,
without needing to know the $d$ vector. Since the $f_{i}$ functions cannot be
found without the $d$ vector, we approximate them using linear models and
write out the equations for the tank level change as

$\dot{{h}}(t)=Ah(t)+B_{1}u(t)+B_{2}d_{a}(t)$ (9)

where $h(t)\in\mathbb{R}^{n}$ includes tank levels, $A\in\mathbb{R}^{n\times
n}$, $B_{1}\in\mathbb{R}^{n\times m}$, $B_{2}\in\mathbb{R}^{n\times 1}$ are
constant system matrices and $d_{a}(t)$ is the aggregated demand of controlled
zone at time $t$, $u(t)\in\mathbb{R}^{m}$ is the input containing pump flows.

##  IV Stochastic Predictive Controller

This section presents a predictive control method for WDNs that are powered by
both PVs and the electrical grid. The method builds on the controller
introduced in [24] and aims to minimize the economic costs associated with
operating the WDN. The problem at time $t$ is formulated as

$\displaystyle\min_{u_{0}^{t},u_{1}^{t}\cdots
u_{N(t)-1}^{t}}\sum_{j=0}^{N(t)-1}\operatorname*{\mathbb{E}}[J(h_{j}^{t},u_{j}^{t},t)]$
(10a) $\displaystyle
h_{j}^{t}=A_{d}h_{j-1}^{t}+B_{d1}u_{j-1}^{t}+B_{d2}d_{a}(j-1)$ (10b)
$\displaystyle h_{0}^{t}=h(t)$ (10c) $\displaystyle
u_{j}^{t}\in\mathcal{U}\subseteq\mathbb{R}^{m}$ (10d) $\displaystyle
h_{j}^{t}\in\mathcal{H}\subseteq\mathbb{R}^{n}$ (10e) $\displaystyle
h_{N(t)}^{t}\in\mathcal{H}_{tf}\subseteq\mathbb{R}^{n}$ (10f)

where $\operatorname*{\mathbb{E}}[J(h_{j}^{t},u_{j}^{t},t)]$ is the expected
value of economic cost function $J(h_{j}^{t},u_{j}^{t},t)$,
$h^{t}=[h_{1}^{t}\cdots h_{N(t)}^{t}]\in\mathbb{R}^{n\times N(t)}$ is the
predicted future states, $u_{j}^{t}\in\mathbb{R}^{m}$ is the input vector,
$N(t)$ is the prediction horizon, $\mathcal{U}\subseteq\mathbb{R}^{m}$ and
$\mathcal{H}\subseteq\mathbb{R}^{n}$ denotes the input and state constraints
respectively and $\mathcal{H}_{tf}\subseteq\mathbb{R}^{n}$ is the terminal
state set. The equation (10b) represents the discretized version of the
continuous system (9). At each time step, which is separated by a time
interval of $\Delta_{t}$, the optimization problem (10) is solved, and the
first term $u_{0}^{t}$ of the optimal input sequence
$\mathbf{u}^{t}=[u_{0}^{t}\cdots u_{N(t)-1}^{t}]\in\mathbb{R}^{m\times N(t)}$
is applied to the system.

The input constraints for the system are determined by the minimum and the
maximum flow rate capacity per unit of time that a pump can deliver. These
conditions are expressed as

$\mathcal{U}=\\{[u_{1}\cdots u_{m}]^{T}\in\mathbb{R}^{m}\mid\forall i:0\leq
u_{i}\leq\overline{u}_{i}\\}$ (11)

where $\overline{u}_{1}\cdots\overline{u}_{m}$ are upper flow limits. Tank
levels are also constrained so that there is always a sufficient reserve of
water on hand in case of emergencies, while also preventing overflow. The set
$\mathcal{H}$ can be defined as

$\mathcal{H}=\\{[h_{1}\cdots h_{n}]^{T}\in\mathbb{R}^{n}\mid\forall
i:\tilde{h}_{i}\leq h_{i}\leq\overline{h}_{i}\\}$ (12)

The linear model (9) does not completely capture the complexity of the entire
network, leading to a model-plant mismatch. That means the actual state
trajectory $h(t)$ might not remain within the set $\mathcal{H}$ even though
the predicted states $h^{t}$ satisfy state constraints. The problem is dealt
with barrier-like exponential functions. First, we rewrite state constraints
(12) as

$C_{i}(h)\leq 0,\quad i=0,1,\cdots 2(n-1)$ (13)

where $C_{0}(h)=\tilde{h}_{1}-h_{1}$ and the rest of the $C_{i}$ functions are
chosen in a similar manner. The cost function terms are then defined as

$J_{h_{i}}(h)=e^{a_{i}(C_{i}(h)+b_{i})}\quad i=0,1,\cdots,2(n-1)$ (14)

where $a_{i},b_{i}\in\mathbb{R}_{>0}$. The parameters $a_{i},b_{i}$ define an
area near the boundaries of the state constraints where the cost function
$J_{hi}$ becomes significantly high. To avoid these high cost values, the
predicted optimal state trajectories $h^{t}$ will steer clear of this region
if possible. This allows the actual states $h(t)$ to remain within the state
constraints (12) assuming the difference between the predicted state and the
actual state is small enough.

The cost function $J(h_{j}^{t},u_{j}^{t},t)$ includes the cost of electricity
purchased from the electrical grid to power the pumps. This electricity is
used to supplement the power generated by the PVs, and is only obtained when
the PV power is not sufficient to operate the pumps. The electricity purchased
from the grid is therefore represented as

$P_{grid}(t)=\max(0,P_{p}(t)-P_{pv}(t))$ (15)

where $P_{p}$ is the total power used by the pumps and $P_{pv}$ is the power
generated by the PV panels and $P_{grid}$ is the power purchased from the
grid. As this function is not differentiable at the origin we approximate it
with the softplus function which can be written as

$sp(x)=\frac{1}{\beta}\log(1+\exp{\beta x})$ (16)

where $\beta$ is a constant. The power provided to the network by the pump $i$
is equal to ${q}_{pi}(p^{out}_{i}-p^{in}_{i})$, where ${q}_{pi}$ is the pump
flow, $p_{out}^{i}$ and $p_{in}^{i}$ are the outlet and inlet pressures of the
pump $i$. The inlet pressures $p^{in}=[p^{in}_{1}~{}p^{in}_{2}]^{T}$ are the
pressures of the related reservoirs and are assumed to be constant. The outlet
pressures $p^{out}=[p^{out}_{1}~{}p^{out}_{2}]^{T}$ are given as the output of
the linear model

$p^{out}(t)=C_{p}h(t)+D_{p}u(t)$ (17)

where $C_{p}$ and $D_{p}$ are found using system identification on data
generated by the EPANET model. The power used by the pumps at time $t$ is then
equal to $P_{p}(t)=u(t)^{T}(p^{out}(t)-p^{in}(t))$. The power generated by the
PVs is equal to $P_{pv}(t)=X_{\tau}^{i}$ where $\tau=(t-t\bmod
T_{day})/T_{day}$ denotes the current day and $i=t/\Delta_{pv}\bmod N_{pv}$ is
the time step of the daily power profile $X_{\tau}$ corresponding to time $t$.
The overall cost function includes both the electricity expense term and the
constraint barrier functions can be expressed as

$J(h(t),u(t),t)=\sum_{i=0}^{\frac{\Delta_{t}}{\Delta_{pv}}-1}c(t)sp(P_{p}(t)-P_{pv}(t+i\Delta_{pv}))\\\
+\sum_{i=0}^{2(n-1)}J_{h_{i}}(h(t))$ (18)

where $c(t)$ is the electricity price. It is assumed that the sampling time
for PV data, $\Delta_{pv}$, is shorter than the system sampling time,
$\Delta_{t}$, resulting in multiple PV power values being available within the
time period from $t$ to $t+\Delta_{t}$. The electrical costs are then
evaluated individually for each distinct PV power value.

In order to address the control problem outlined in equation (10), we use a
scenario-based stochastic programming approach, in which we generate
$S\in\mathbb{N}$ number of scenarios. The only stochastic variable appearing
in the cost function $J(h(t),u(t),t)$ is the PV production data $X^{i}_{\tau}$
which is 0 during nighttime. The values of $X^{i}_{\tau}$ during daytime are
represented by the equation
$X^{i}_{\tau}=p_{\tau}Y^{i}_{\tau}+\delta^{i}_{\tau}$,
$\delta^{i}_{\tau}=\phi^{\delta^{i}_{\tau}}\delta^{i-1}_{\tau}+\epsilon^{\delta^{i}_{\tau}}$
where the values of $Y^{i}_{\tau}$ and
$\phi^{\delta^{i}_{\tau}}\delta^{i-1}_{\tau}$ have been previously calculated,
as explained in Section II. We use different sampling processes for problems
formulated during nighttime and daytime. When the optimization problem (10) is
formed at night, the values of $X^{i}_{\tau}$ are generated by sampling
$p_{\tau}$ using the equation $\hat{p}_{\tau}=p_{\tau}-\epsilon_{\tau}$ where
$\hat{p}_{\tau}$ is the estimated value and $\epsilon_{\tau}$ is normally
distributed with known mean and variance. Additionally, error values
$\delta^{i}_{\tau}$ are calculated using
$\delta^{i}_{\tau}=\phi^{\delta^{i}_{\tau}}\delta^{i-1}_{\tau}+\epsilon^{\delta^{i}_{\tau}}$,
where $\epsilon^{\delta^{i}_{\tau}}$ is a normally distributed value with a
mean of 0 and a known variance $\sigma^{\delta^{i}_{\tau}}$. After sunrise, we
start to receive non-zero power production data $X^{i}_{\tau}$, which provides
us with more information about the variables $\epsilon^{\delta^{i}_{\tau}}$
and $p_{\tau}$. Let the index corresponding to sunrise be denoted as $s$. At
time $i_{c}$ the $\epsilon^{\delta^{i}_{\tau}},i=s+1,s+2\cdots i_{c}$ values
become deterministic functions of $p_{\tau}$ given by the equations
$\epsilon^{\delta^{i}_{\tau}}=g_{i}(p_{\tau})=X^{i}_{\tau}-p_{\tau}Y^{i}_{\tau}-\phi\delta^{i-1}_{\tau}$
for $i=s+1,s+2\cdots i_{c}$ as the values of $X^{i}_{\tau},i=s+1,s+2\cdots
i_{c}$ are available at time $i_{c}$. Given this information, the probability
distribution of $p_{\tau}$ is updated using

$\displaystyle
P(p_{\tau}=x\mid\epsilon^{\delta^{s+1}_{\tau}}=X^{s+1}_{\tau}-p_{\tau}Y^{s+1}_{\tau}-\phi\delta^{s}_{\tau},$
(19a)
$\displaystyle\cdots,\epsilon^{\delta^{i_{c}}_{\tau}}=X^{i}_{\tau}-p_{\tau}Y^{i_{c}}_{\tau}-\phi\delta^{i_{c}-1}_{\tau})\propto
P(p_{\tau}=x,$
$\displaystyle\epsilon^{\delta^{s+1}_{\tau}}=g_{s+1}(x),\cdots,\epsilon^{\delta^{i_{c}}_{\tau}}=g_{i_{c}}(x))$
$\displaystyle
P(p_{\tau}=x,\epsilon^{\delta^{s+1}_{\tau}}=X^{s+1}_{\tau}-xY^{s+1}_{\tau}-\phi\delta^{s}_{\tau},$
(19b)
$\displaystyle\cdots\epsilon^{\delta^{i_{c}}_{\tau}}=X^{i}_{\tau}-xY^{i_{c}}_{\tau}-\phi\delta^{i_{c}-1}_{\tau})=P(p_{\tau}=x)$
$\displaystyle P(\epsilon^{\delta^{s+1}_{\tau}}=g_{s+1}(x))\cdots
P(\epsilon^{\delta^{i_{c}}_{\tau}}=g_{i_{c}}(x))$

The relation (19a) follows from the conditional probability $P(A\mid
B)=P(A,B)/P(B)$. In this case the probability of the given information
$P(B)=P(\epsilon^{\delta^{s+1}_{\tau}}=g_{s+1}(p_{\tau}),\cdots,\epsilon^{\delta^{i_{c}}_{\tau}}=g_{i_{c}}(p_{\tau}))$
is just a constant number, so (19a) is expressed as $P(A\mid B)\propto
P(A,B)$. The equation (19b) follows from the independency of $p_{\tau}$ and
$\epsilon^{\delta^{i}_{\tau}}$ terms. After sunrise, the $p_{\tau}$ is sampled
using the new distribution coming from (19b) and the error values
$\epsilon^{\delta^{i}_{\tau}},i=s+1,s+2\cdots i_{c}$ are computed directly as
they are functions of $p_{\tau}$. The remaining values of
$\epsilon^{\delta^{i}_{\tau}},i=i_{c}+1\cdots$ are sampled as they are sampled
at night.

The PV power production, electricity price $c(t)$ and total water demand
$d_{a}(t)$ signals can be seen as comprising a periodic signal with a period
of 1 day and a relatively minor deviation signal. This is leveraged to
increase the chance of finding a solution to the optimization problem (10).
The idea is to keep the tank levels at the start of each day within a certain
range. If a pumping schedule works for the first day, then it is probable that
a similar schedule could be applied to the other days since the initial
conditions will be similar. Therefore, the horizon $N(t)$ is chosen to target
the beginning of each day and expressed as

$N(t)=(T_{day}-t\bmod T_{day})/\Delta_{t}$ (20)

where $T_{day}$ is the duration of a whole day. The next step is to figure out
which tank levels the trajectories should return to at the start of each day.
We define the optimal periodic trajectory of the system as the solution of

$\displaystyle(\mathbf{u^{*}},\mathbf{h^{*}})=\operatorname*{arg\,min}_{{u_{i}},{h_{i}}}\sum_{i=0}^{(T_{day}/\Delta_{t})-1}J({h_{i}},{u_{i}})$
(21a)
$\displaystyle{h_{i}}=A_{d}{h_{i-1}}+B_{d1}{u_{i-1}}+B_{d2}d_{a}^{*}(i-1)$
(21b) $\displaystyle{u_{i}}\in\mathcal{U}\subseteq\mathbb{R}^{m}$ (21c)
$\displaystyle{h_{i}}\in\mathcal{H}\subseteq\mathbb{R}^{n}$ (21d)
$\displaystyle{h_{0}}=h_{T_{day}/\Delta_{t}}$ (21e)

where $d_{a}^{*}$ is the average daily demand profile obtained from the past
measurements. The resulting state trajectory $\mathbf{h^{*}}=[h^{*}_{0}\cdots
h^{*}_{T_{day}/\Delta_{t}}]\in\mathbb{R}^{n\times(T_{day}/\Delta_{t}+1)}$ is
the optimal periodic trajectory because of the periodicity constraint (21e).
Hence, it is sensible to drive tank levels close to
$h_{T_{day}/\Delta_{t}}^{*}$ at the end of each day. The terminal state
constraint set $\mathcal{H}_{tf}$ then can be written as

$\mathcal{{H}}_{tf}=\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$ (22)

where $\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$ is the open ball centered
at $h_{T_{day}/\Delta_{t}}^{*}$ with radius $r$. If the problem (10) becomes
infeasible at any time step $t$, we apply the second term of the input
sequence from the previous step $u_{1}^{t-\Delta_{t}}$.

##  V Application

The simulations were carried out utilizing the EPANET software, which is
widely used in the industry for creating realistic models and simulating the
behavior of networks with a specific pumping strategy. In this particular
case, the discussed method is applied to the EPANET implementation of the WDN
of Randers, a Danish city, as illustrated in Figure 1.

High ZoneLow Zone

Figure 1: Water Distribution Network of Randers. The pumping stations to be
controlled are shown in red and the remaining in blue. Tanks are shown with a
’T’ shaped symbol in yellow.

The network contains 4549 nodes and 4905 links connecting them and includes 8
pumping stations. Out of these, 6 are shown in the figure while the remaining
2 are located where tanks are placed. Our goal is to determine the schedules
for two pumping stations, while the others operate according to pre-determined
strategies. The stations that will be controlled are highlighted in red in the
figure, and they provide water mostly to the High Zone (HZ) and Low Zone (LZ).
Additionally, there are 3 tanks in the HZ, with two of them connected via
pipes and the third standing alone.

As the pipe connecting the two connected tanks is large enough, water levels
at the two connected tanks are almost equal $h_{1}\thickapprox h_{2}$ all the
time. That enables us to consider $h_{1},h_{2}$ together as one state variable
${h}_{1,2}$. We have utilized the EPANET model of the network to generate the
data needed to approximate the flow of the pipes connected to the tanks as
linear functions of the tank levels
$h(t)=[{h}_{1,2}(t)~{}{h}_{3}(t)]\in\mathbb{R}^{2}$, pump flows
$u(t)=[{q}_{p1}(t)~{}{q}_{p2}(t)]\in\mathbb{R}^{2}$ and the aggregated demand
$d_{a}(t)$. The total demand for the High and Low Zones is used as the
aggregated demand in the model since these areas are primarily supplied by the
controlled pumps. The model is simulated with various initial tank level
conditions and flow rates of the 2 pumping stations that are being controlled,
while the control laws for the remaining pumping stations have already been
defined in the EPANET model. After fitting the EPANET data to the linear model
of the pipe flows, we use the tank level equation (8) to obtain the dynamic
model of the system.

###  V-A Simulation Results

We test our proposed control method on the EPANET model of the Randers water
network. Using the Epanet-Matlab toolkit [25], we simulate the network by
setting the flow of the two pumps at each time step and controlling the
remaining pumps with previously defined rule-based control laws.

We set the parameters of exponential barrier functions $J_{hi}$ to $a_{i}=80,$
$b_{i}=0.3$ for all $i$, and the $\beta$ parameter of the softplus function to
$0.02$. It is assumed that the electricity prices are known in advance during
the test. The maximum tank levels are set to 3m for $h_{1}$ and $h_{2}$, and
2.8m for $h_{3}$, while the minimum tank level is set to half-full. The
maximum pump flow is set to 100 and the sampling time $\Delta_{t}$ to 1 hour,
so the control input is recalculated at each hour. The experiments were
conducted using a PV power production data set recorded in Albuquerque, US
[26]. The power values are adjusted by multiplying them with a specific
number, in order to align the average power generation throughout the entire
data set with the average power consumption of the controlled pumps on a
typical day. The average power consumption of the pumps is obtained by
simulating the EPANET model for one day. We assume that total demand
$d_{a}(t)$ of HZ and LZ can be estimated up to 1 day in advance from the
available data. We don’t have any actual historical data to measure the
demand, so we use a modified version of the actual demand found in the EPANET
simulation during MPC calculations. This modified version is based on a real
demand data set from a Danish facility, where the difference between the
average demand and the demand on a particular day is added to the EPANET
demand. This allows us to replicate an estimated demand. In each experiment a
different day from the data set is used, so the assumed estimated demand is
different each time.

The EPANET model of the network is run with the presented method and the
results are given in Figure 2. The data for PV production is obtained by
selecting a single day from the PV data set at random. To make the most
efficient use of electricity costs, a moderate amount of water is pumped
during the night when electricity prices are low, in order to maintain water
levels above the threshold until the PV power can be utilized to pump the
water. This strategy can be observed by looking at the water level of $h_{3}$
around hour 9, where despite the pumping that takes place at night, the water
level still approaches close to the lower threshold. To take advantage of the
PV energy available during the day, the majority of water is pumped during
this period. Once the sun sets and electricity prices decrease again, water
pumping continues steadily to ensure that the terminal water level constraints
are met. A total of $66.95\%$ of the energy needed for the operation of the
pumps is provided by the PV power.

(a)

(b)

(c)

(d)

(e)

Figure 2: Sample simulation. (a) evolution of tank levels through 1 day with
upper and lower level thresholds; (b) real total demand of HZ and LZ used in
EPANET simulation and the demand used in MPC calculations; (c) total flow
provided by the 2 pumps; (d) electricity price; (e) PV power production and
the power consumed by the pumps.

The proposed method is then compared with its deterministic optimization
counterpart. The key distinction between the two is that the cost function
used in the deterministic method is $J(h_{j}^{t},u_{j}^{t},t)$ rather than its
expected value $\operatorname*{\mathbb{E}}[J(h_{j}^{t},u_{j}^{t},t)]$. This
means that instead of using a scenario-based optimization, direct estimates of
PV power production are used for the calculation of the cost function. To
compare the performance of the two methods, the EPANET model is run for 10
days using both methods separately while keeping the demand patterns, PV
production profiles, and initial tank levels constant. The PV production
profiles used in the experiment are obtained from the dataset by selecting 10
consecutive days. To ensure the robustness and generality of the results, the
experiment is repeated 4 times, each with a different randomly selected PV
starting day, in order to cover different PV production patterns. The ratios
of the total energy costs and the energy used from the grid in the stochastic
optimization method with respect to those in the deterministic optimization
method are given in Table I for all 4 cases. In a total of 40 days, electrical
costs for the stochastic optimization were $95.92\%$ of the costs for the
deterministic optimization. While the costs for stochastic optimization were
found to be lower, it should be noted that more energy was drawn from the
grid. This is due to the strategy of pumping more water during the night when
electricity prices are lower, and less during the day when PV power is
available. The reasoning behind this is to avoid purchasing electricity at
high prices during times when PV production may fall short of predictions. By
using less PV energy and pumping more during the night, stochastic
optimization aims to minimize costs. This strategy is particularly beneficial
on cloudy days when PV power production is lower than expected. An
illustration of the energy consumption for both deterministic and stochastic
optimization scenarios on a cloudy day can be seen in Figure 3.

|  Ratio of Electrical Costs |  Ratio of Energy Used From The Grid
---|---|---
Case 1 |  0.9290 |  1.1138
Case 2 |  0.9759 |  1.0682
Case 3 |  0.9734 |  1.1323
Case 4 |  0.9729 |  1.0829
Total |  0.9592 |  1.1004

TABLE I: The ratios of the total energy costs and the energy used from the
grid between stochastic optimization (SO) method and deterministic
optimization (DO) method $SO/DO$. Figure 3: Comparisons of power consumed by
pumps with deterministic optimization (DO) and stochastic optimization (SO)
along with the power profile and its estimation before the sunset

##  VI Conclusion

A method for scheduling the pumps of WDNs powered by grid-connected PVs is
proposed to reduce economic costs while taking into account the uncertainties
of PV power production through stochastic optimization. Our findings showed
that $66.95\%$ of the energy needed for pump operation was supplied by PV
power. The integration of PVs resulted in a significant reduction in
electrical costs, with the use of stochastic optimization bringing an
additional $4\%$ cost reduction compared to deterministic optimization.
Stochastic optimization is most effective on days with high uncertainties.
However, it does draw more energy from the grid than deterministic
optimization, as more water is pumped from the grid when electricity prices
are low instead of pumping to the limit when PV power is present.

## References

  * [1] N. Sharif, H. Haider, A. Farahat, K. Hewage, and R. Sadiq, “Water energy nexus for water distribution systems: A literature review,” _Environmental Reviews_ , vol. 27, 03 2019.
  * [2] F. K. Pour, V. Puig, and G. Cembraño, “Economic mpc-lpv control for the operational management of water distribution networks,” _IFAC-PapersOnLine_ , 2019.
  * [3] F. Fiedler, A. Cominola, and S. Lucia, “Economic nonlinear predictive control of water distribution networks based on surrogate modeling and automatic clustering,” _IFAC-PapersOnLine_ , vol. 53, pp. 16 636–16 643, 2020.
  * [4] K. M. H. Baunsgaard, O. Ravn, C. S. Kallesøe, and N. K. Poulsen, “Mpc control of water supply networks,” _2016 European Control Conference (ECC)_ , pp. 1770–1775, 2016.
  * [5] Y. Wang, T. Alamo, V. Puig, and G. Cembraño, “Economic model predictive control with nonlinear constraint relaxation for the operational management of water distribution networks,” _Energies_ , vol. 11, p. 991, 2018.
  * [6] C. S. Kallesøe, T. N. Jensen, and J. D. Bendtsen, “Plug-and-play model predictive control for water supply networks with storage,” _IFAC-PapersOnLine_ , vol. 50, pp. 6582–6587, 2017.
  * [7] Y. Zhang, C. Jin, R. K. Sharma, and A. K. Srivastava, “Data-driven day-ahead pv estimation using hybrid deep learning,” _2019 IEEE Industry Applications Society Annual Meeting_ , pp. 1–6, 2019.
  * [8] S. Theocharides, G. Makrides, A. Livera, M. Theristis, P. Kaimakis, and G. E. Georghiou, “Day-ahead photovoltaic power production forecasting methodology based on machine learning and statistical post-processing,” _Applied Energy_ , vol. 268, p. 115023, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0306261920305353
  * [9] Y. Zhang, C. Jin, R. K. Sharma, and A. K. Srivastava, “Data-driven day-ahead pv estimation using hybrid deep learning,” in _2019 IEEE Industry Applications Society Annual Meeting_ , 2019, pp. 1–6.
  * [10] A. Mellit, A. M. Pavan, and V. Lughi, “Short-term forecasting of power production in a large-scale photovoltaic plant,” _Solar Energy_ , vol. 105, pp. 401–413, 2014.
  * [11] T. Huld, G. Friesen, A. Skoczek, R. P. Kenny, T. Sample, M. Field, and E. D. Dunlop, “A power-rating model for crystalline silicon pv modules,” _Solar Energy Materials and Solar Cells_ , vol. 95, no. 12, pp. 3359–3369, 2011. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0927024811004442
  * [12] F. Bizzarri, M. Bongiorno, A. Brambilla, G. Gruosso, and G. S. Gajani, “Model of photovoltaic power plants for performance analysis and production forecast,” _IEEE Transactions on Sustainable Energy_ , vol. 4, no. 2, pp. 278–285, 2013.
  * [13] F. E. Atencio Espejo, S. Grillo, and L. Luini, “Photovoltaic power production estimation based on numerical weather predictions,” in _2019 IEEE Milan PowerTech_ , 2019, pp. 1–6.
  * [14] J. Á. G. Ordiano, S. Waczowicz, M. Reischl, R. Mikut, and V. Hagenmeyer, “Photovoltaic power forecasting using simple data-driven models without weather data,” _Computer Science - Research and Development_ , vol. 32, pp. 237–246, 2017.
  * [15] M. G. De Giorgi, P. Congedo, and M. Malvoni, “Photovoltaic power forecasting using statistical methods: Impact of weather data,” _Science, Measurement & Technology, IET_, vol. 8, pp. 90–97, 05 2014.
  * [16] M. Park, J. Lee, W. Kang, J. Choi, and K. H. Lee, “Predictive model for pv power generation using rnn (lstm),” _Journal of Mechanical Science and Technology_ , vol. 35, pp. 1–9, 01 2021.
  * [17] H. Chen and X. Chang, “Photovoltaic power prediction of lstm model based on pearson feature selection,” _Energy Reports_ , vol. 7, pp. 1047–1054, 2021, 2021 International Conference on Energy Engineering and Power Systems. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2352484721009768
  * [18] A. Parisio, C. Wiezorek, T. Kyntaja, J. Elo, and K. H. Johansson, “An mpc-based energy management system for multiple residential microgrids,” _2015 IEEE International Conference on Automation Science and Engineering (CASE)_ , pp. 7–14, 2015.
  * [19] A. Parisio, C. Wiezorek, T. Kyntäjä, J. Elo, K. Strunz, and K. H. Johansson, “Cooperative mpc-based energy management for networked microgrids,” _IEEE Transactions on Smart Grid_ , vol. 8, pp. 3066–3074, 2017\.
  * [20] J. Shen, C. Jiang, Y. Liu, and X. Wang, “A microgrid energy management system and risk management under an electricity market environment,” _IEEE Access_ , vol. 4, pp. 2349–2356, 2016.
  * [21] Y. Xiang, J. Liu, and Y. Liu, “Robust energy management of microgrid with uncertain renewable generation and load,” _IEEE Transactions on Smart Grid_ , vol. 7, pp. 1034–1043, 2016.
  * [22] Z. Luo, W. Gu, Z. Wu, Z. Wang, and Y. Tang, “A robust optimization method for energy management of cchp microgrid,” _Journal of Modern Power Systems and Clean Energy_ , vol. 6, pp. 132–144, 2018.
  * [23] M. Ürkmez, C. Kallesøe, J. Dimon Bendtsen, and J. Leth, “Day-ahead pv power forecasting for control applications,” in _IECON 2022 – 48th Annual Conference of the IEEE Industrial Electronics Society_ , 2022, pp. 1–6.
  * [24] M. Ürkmez, C. Kallesøe, J. D. Bendtsen, and J. Leth, “Economic predictive control with periodic horizon for water distribution networks,” 2023. [Online]. Available: https://arxiv.org/abs/2301.13598
  * [25] D. G. Eliades, M. Kyriakou, S. Vrachimis, and M. M. Polycarpou, “Epanet-matlab toolkit: An open-source software for interfacing epanet with matlab,” in _Proc. 14th International Conference on Computing and Control for the Water Industry (CCWI)_ , The Netherlands, Nov 2016, p. 8.
  * [26] C. W. Hansen, W. F. Holmgren, A. Tuohy, J. Sharp, A. T. Lorenzo, L. J. Boeman, and A. Golnas, “The solar forecast arbiter: An open source evaluation framework for solar forecasting,” in _2019 IEEE 46th Photovoltaic Specialists Conference (PVSC)_ , 2019, pp. 2452–2457.

## ACKNOWLEDGMENT

This work is funded by Independent Research Fund Denmark (DFF). We acknowledge
Verdo company, Peter Nordahn, and Steffen Schmidt for providing us with the
EPANET model and the network information.


11institutetext: IE School of Sciences and Technology, P.º Castellana 259,
28029, Madrid, Spain 11email<EMAIL_ADDRESS>22institutetext:
Instituto de Ciencias Matematicas (CSIC-UAM-UC3M-UCM), Calle Nicolas Cabrera
13-15, 28049 Madrid, Spain 22email<EMAIL_ADDRESS>
33institutetext: LAGEPP, UCB Lyon 1 - CNRS UMR 5007, CPE Lyon - Bâtiment 308
G, Université Claude Bernard Lyon-1, 43, Bd. du 11 Novembre 1918, F-69622
Villeurbanne, France 33email<EMAIL_ADDRESS>
# Discrete gradient methods
for irreversible port-Hamiltonian systems

Alexandre Anahory Simoes 11 David Martín de Diego 2*2*
Bernhard Maschke 33

###### Abstract

In this paper we introduce discrete gradient methods to discretize
irreversible port-Hamiltonian systems showing that the main qualitative
properties of the continuous system are preserved using this kind
discretizations methods.

rreversible port-Hamiltonian system, discrete gradient system, geometric
integration, thermodynamics.

###### Keywords:

I

##  1 Introduction

In [1, 2, 3] irreversible port-Hamiltonian Systems have been suggested as a
class of port-Hamiltonian systems [4] adapted to represent open physical
systems in a way compatible with irreversible Thermodynamics [5, chap.14].
This formulation uses two thermodynamic potentials, the total energy and the
entropy functions, but is formulated as a Hamiltonian system defined on an
almost-Poisson bracket, contrary to other formulations which use a
metriplectic structure [6, 7]. Furthermore is encompasses open thermodynamic
systems by defining port variables and port maps associated with reversible
and irreversible phenomena at the interface of the system [8]. In this paper,
we shall define geometric integrators of these irreversible port-Hamiltonian
systems which preserve the energy and entropy balance equations as it has been
done for contact Hamiltonian systems [9] or metriplectic systems [10]. We
shall use discrete gradient methods that are typically designed for numerical
integration of ordinary differential equations (ODEs), while ensuring the
preservation of certain structures of the continuous system as for instance
energy preservation, dissipation, as it has also been used for Port
Hamiltonian Systems [11, 12]. We will show, that the suggested method
preserves both the invariance of energy and the irreversibility property of
original system (see [13, 14, 15] and references therein).

##  2 Irreversible Port-Hamiltonian systems

Irreversible Port Hamiltonian Systems have been introduced in [1, 2], as a
class of control systems for irreversible thermodynamic systems. Their
Hamiltonian function is the total energy of the system and they are defined
with respect to an almost-Poisson bracket which allows to embed the second
principle of Irreversible Thermodynamics, namely the irreversible entropy
creation term due to the irreversible phenomena. Furthermore they represent
also, open systems by expressing the interactions of the system with its
environment by using pairs of conjugated port variables, called port
variables. The definition given below, considers a more precise expression,
introduced in [8], than in [1, 2], of the input map associated with a
irreversible phenomena occurring at the interface of the system such as the
heat conduction.

###  2.1 Definition

###### Definition 1

An irreversible port-Hamiltonian Systems (IPHS) is the nonlinear control
system

$\displaystyle\left(\begin{array}[]{c}\frac{dx}{dt}\\\ -y\end{array}\right)=$
$\displaystyle\left[\sum_{i=1}^{k}\gamma_{i}\left(x,\tfrac{\partial
H}{\partial
x}\right)\left\\{S,H\right\\}_{J_{i}}\left(\begin{array}[]{cc}J_{i}&0\\\
0&0\end{array}\right)+\sum_{\alpha=1}^{l}\left(\begin{array}[]{cc}M_{\alpha}&0\\\
0&0\end{array}\right)\right.$ (2.7)
$\displaystyle+\sum_{j=1}^{\bar{k}}\gamma_{\textrm{port},j}\left(x,\tfrac{\partial
H}{\partial
x},u\right)\left\\{S_{\textrm{tot}},H_{\textrm{tot}}\right\\}_{J_{\textrm{port},g_{j}}}\left(\begin{array}[]{cc}0&g_{j}\\\
-g_{j}^{\top}&0\end{array}\right)$ (2.10)
$\displaystyle\left.+\sum_{\beta=1}^{\bar{l}}\left(\begin{array}[]{cc}0&g_{S\beta}\\\
-g_{S\beta}^{\top}&0\end{array}\right)\right]\left(\begin{array}[]{c}\frac{\partial
H}{\partial x}\\\ u\end{array}\right)$ (2.15)

where $x\left(t\right)\in\mathbb{R}^{n}$ is the _state vector_ ,
$u\left(t\right)\in\mathbb{R}^{m}$ is the _control input_ vector, and defined
by

  1. (i)

two (smooth) real functions called _Hamiltonian function_ $H\in
C^{\infty}(\mathbb{R}^{n})$ and _entropy function_ $S\in
C^{\infty}(\mathbb{R}^{n})$

  2. (i)

two sets of structure matrices $J_{i}(x)\in\mathbb{R}^{n\times n}$, $1\leq
i\leq k$ and $M_{\alpha}(x)\in\mathbb{R}^{n\times n}$, $1\leq\alpha\leq m$,
defining _almost-Poisson brackets_ $\\{\;,\;\\}_{J_{i}}$ and
$\\{\;,\;\\}_{M_{\alpha}}$, respectively, and the _real, strictly positive
functions_ $\gamma_{i}(x,\frac{\partial H}{\partial x})=\hat{\gamma}_{i}(x)\in
C^{\infty}(\mathbb{R}^{n},\,\mathbb{R}_{+}^{*})$ .

  3. (iii)

the irreversible port-maps are defined by the _input matrices_ $g_{j}(x)$ ,
$1\leq j\leq\bar{k}$, and the real strictly positive functions.
$\gamma_{\textrm{port},j}\left(x,\frac{\partial H}{\partial
x},u\right)=\widehat{\gamma_{\textrm{port},j}}(x,u)\in
C^{\infty}(\mathbb{R}^{n+m},\,\mathbb{R}_{+}^{*})$ . The _port bracket_ is
defined with respect to the additional functions
$H_{\textrm{tot}}\left(x,\xi\right)=H\left(x\right)+u^{\top}\xi\;,\xi\in\mathbb{R}^{m}$
and
$S_{\textrm{tot}}\left(x,\xi\right)=S\left(x\right)+\tau^{\top}\xi\;,\tau\in\mathbb{R}^{m}\;$
as

$\left\\{S_{\textrm{tot}},H_{\textrm{tot}}\right\\}_{J_{\textrm{port},g}}=\left[\left(g^{\top}\frac{\partial
S}{\partial x}\right)^{\top}u-\tau^{\top}\left(g^{\top}\frac{\partial
H}{\partial x}\right)\right]$

  4. (iii)

the reversible port-maps defined by the _input matrices_ $g_{S\beta}(x)$ in
$\mathbb{R}^{n\times m}$, $1\leq\beta\leq\bar{l}$ . We assume additionally
that
$\left\\{S_{\textrm{tot}},H_{\textrm{tot}}\right\\}_{J_{\textrm{port},g_{S\beta}}}=0$
for all $1\leq\beta\leq\bar{l}$

###### Remark 1

Note that the definition given here is more general than the one given in [8],
in the sense that the structure matrices of the port maps may be dependent on
the state.

From the definition of the Irreversible Port Hamiltonian Systems, one deduces
the following energy and entropy balance equations

$\displaystyle\frac{dH}{dt}-y^{\top}u=0$
$\displaystyle\frac{dS}{dt}-y^{\top}\tau=\sum_{i}\hat{\gamma_{i}}\left(x\right)\left\\{S,H\right\\}_{J_{i}}^{2}+\sum_{j}\widehat{\gamma_{\textrm{port},j}}\left(x,u\right)\left\\{S_{\textrm{tot}},H_{\textrm{tot}}\right\\}_{J_{\textrm{port},g_{j}}}^{2}\geq
0$

where the right-hand side of the second equation is equal to the irreversible
entropy creation due to the irreversible phenomena in the system and at its
boundary.

###  2.2 Example: gas-piston system

This system describes the thermodynamic properties of a gas contained in a
cylinder closed by a piston submitted to gravity ([2]) and an an external
force $u_{2}$, the gas is subject to an exchange of heat with a thermostat at
temperature $u_{1}$, obeying Fourier’s heat conduction law with coefficient
$\lambda_{e}>0$. Mathematically this system is described giving the mechanical
Hamiltonian $H_{mec}(q,p)=\frac{1}{2m}p^{2}+mgq\;.$ where $q$ describes the
altitude of the piston and $p$ its associated momentum. On the other hand, the
properties of the perfect gas may be defined by its internal energy $U(S,V)$
where $S$ is the entropy variable and $V$ the volume variable. Therefore, the
total energy if

$H(S,V,q,p)=U(S,V)+H_{mec}(q,p)$ (2.16)

Denote by $\frac{\partial H}{\partial S}=T,\quad\frac{\partial H}{\partial
V}=-P,\quad\frac{\partial H}{\partial q}=mg,\quad\frac{\partial H}{\partial
p}=v,\quad$ the corresponding co-energy variables, where $T$ represents the
temperature of the gas and $P$ its pressure, $mg$ is the gravity force and $v$
is the velocity of the piston. We assume that when the piston moves, a non-
adiabatic transformation due to mechanical friction and the viscosity of the
gas, induces a heat flow in the gas only, leading to the entropy balance law
$\frac{dS}{dt}=\frac{1}{T}\mu v^{2}\geq 0$ parametrized by $\mu\geq 0$.

The equations of motion of the gas-piston system can be written as:

$\left(\begin{array}[]{r}dS/dt\\\ dV/dt\\\ dq/dt\\\
dp/dt\end{array}\right)=\left(\begin{array}[]{cccc}0&0&0&\mu v/T\\\ 0&0&0&A\\\
0&0&0&1\\\ -\mu v/T&-A&-1&0\end{array}\right)\left(\begin{array}[]{r}T\\\
-P\\\ mg\\\
v\end{array}\right)+\left(\begin{array}[]{c}\lambda_{e}\left(\frac{1}{T}-\frac{1}{u_{1}}\right)\\\
0\\\ 0\\\ u_{2}\end{array}\right)$

where $A$ denotes the area of the piston. This system is an Irreversible Port-
Hamiltonian Systems generated by the Hamiltonian (2.16) and the entropy
function, denoted with an abuse of notation $S(S,V,q,p)=S$. The reversible
coupling between the kinetic and potential mechanical energy as well as with
the gas through the piston are given by the skew-symmetric matrix $M$ and the
coupling between the mechanical displacement and the heat flow due to the
dissipation is given by the skew-symmetric matrix $J$ defined by

$J=\left(\begin{array}[]{cccc}0&0&0&1\\\ 0&0&0&0\\\ 0&0&0&0\\\
-1&0&0&0\end{array}\right),\ M=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&0&0&A\\\ 0&0&0&1\\\ 0&-A&-1&0\end{array}\right)$

One may check that $\nabla S=(1,0,0,0)\in\ker M$, that is, the entropy
function is a Casimir function for $M$. Furthermore the Poisson backet is
$\\{S,H\\}_{J}=v$ which is the velocity of the piston and is the driving force
of the dissipation. The positive function giving the constitutive relation of
the dissipation is $\hat{\gamma}(S,V,q,p)=\frac{\mu}{T}\geq 0\;$ .

The input map are defined by the matrices
$J_{\textrm{port}}=\left(\begin{array}[]{cc}0&g\\\
-g^{\top}&0\end{array}\right)$ and
$M_{\textrm{port}}=\left(\begin{array}[]{cc}0&g_{S}\\\
-g_{S}^{\top}&0\end{array}\right)$ where
$g^{\top}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&0&0&0\end{array}\right)\;$
and $g_{S}^{\top}=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&0&0&1\end{array}\right).$ The extended Hamiltonian and entropy functions
are:

$\displaystyle H_{tot}(S,V,q,p,\xi_{1},\xi_{2})$ $\displaystyle=$
$\displaystyle H(S,V,q,p)+\xi_{1}u_{1}+\xi_{2}u_{2}$ $\displaystyle
S_{tot}(S,V,q,p,\xi_{1},\xi_{2})$ $\displaystyle=$ $\displaystyle S+\xi_{1}$

where we fixed the value $\tau_{1}=1$ associated with the port where the
irreversible phenomenon, the heat conduction with the thermostat at
temperature $u_{1}$, takes place with the positive function
$\gamma_{port}=\frac{\lambda_{e}}{Tu_{1}^{2}}$ and we fixed $\tau_{2}=0$ for
the port where the external force $u_{2}$ acting on the piston in a reversible
way. Observe that $\\{S_{tot},H_{tot}\\}_{J_{port}}=u_{1}-T$ and moreover
$\nabla S_{tot}$ is a Casimir of $M_{port}$.

The energy and the entropy balance equation become
$\frac{dH}{dt}=y_{1}u_{1}+y_{2}u_{2}$ and
$\frac{dS}{dt}-y_{1}=\frac{(u_{1}-T)^{2}}{Tu_{1}^{2}}+\frac{1}{T}\mu v^{2}\geq
0\;$.

###  2.3 Discrete gradient methods

For ODEs in skew-gradient form, i.e. $\dot{x}=J(x)\frac{\partial H}{\partial
x}(x)$ with $x\in\mathbb{R}^{n}$ and $J(x)$ a skew-symmetric matrix, we can
use discretizations of the gradient $\nabla H(x)\equiv\frac{\partial
H}{\partial x}$ to define a class of numerical integrators which exactly
preserve the first integral $H$.

###### Definition 2

[13] Let $H:\mathbb{R}^{n}\longrightarrow\mathbb{R}$ be a differentiable
function. Then $\bar{\nabla}H:\mathbb{R}^{2n}\longrightarrow\mathbb{R}^{n}$ is
a discrete gradient of $H$ if it is continuous and satisfies

$\displaystyle\bar{\nabla}H(x,x^{\prime})^{T}(x^{\prime}-x)$
$\displaystyle=H(x^{\prime})-H(x)\,,\quad\,\mbox{ for all
}x,x^{\prime}\in\mathbb{R}^{n}\,,$ (2.17a) $\displaystyle\bar{\nabla}H(x,x)$
$\displaystyle=\nabla H(x)\,,\quad\quad\quad\quad\mbox{ for all
}x\in\mathbb{R}^{n}\,.$ (2.17b)

Some examples of discrete gradient methods can be find in [13], [14] and [15].

Once a discrete gradient $\bar{\nabla}H$ has been chosen, it is
straightforward to define an energy-preserving integrator by
$\frac{x^{\prime}-x}{h}=\tilde{J}(x,x^{\prime},h)\bar{\nabla}H(x,x^{\prime})$
where $\tilde{J}$ is a differentiable skew-symmetric matrix approximating $J$.
In particular it satisfies $\tilde{J}(x,x,0)=J(x)$. As in the continuous case,
it is immediate to check that $H$ is exactly preserved since

$H(x^{\prime})-H(x)=\bar{\nabla}H(x,x^{\prime})^{\top}(x^{\prime}-x)=h\bar{\nabla}H(x,x^{\prime})^{\top}\tilde{J}(x,x^{\prime},h)\bar{\nabla}H(x,x^{\prime})=0\,.$

Some well-known examples of discrete gradients are:

  * •

The mean value (or averaged) discrete gradient introduced in [16] and given by

$\bar{\nabla}_{1}H(x,x^{\prime}):=\int_{0}^{1}\nabla
H((1-s)x+sx^{\prime})\,ds$

for $x^{\prime}\not=x$.

  * •

The midpoint (or Gonzalez) discrete gradient, introduced in [13] and given by

$\displaystyle\tiny\bar{\nabla}_{2}H(x,x^{\prime})$ $\displaystyle:=\nabla
H\left(\frac{1}{2}(x^{\prime}+x)\right)+\tfrac{H(x^{\prime})-H(x)-\nabla
H\left(\frac{1}{2}(x^{\prime}+x)\right)^{T}(x^{\prime}-x)}{|x^{\prime}-x|^{2}}(x^{\prime}-x)\,,$
$\displaystyle\mbox{ for }x^{\prime}\not=x\,.$

  * •

The coordinate increment discrete gradient, introduced in [14], with each
component given by

$\bar{\nabla}_{3}H(x,x^{\prime})_{i}:=\tfrac{H(x^{\prime}_{1},\ldots,x^{\prime}_{i},x_{i+1},\ldots,x_{n})-H(x^{\prime}_{1},\ldots,x^{\prime}_{i-1},x_{i},\ldots,x_{n})}{x^{\prime}_{i}-x_{i}}\,,\quad
1\leq i\leq N\,,$

when $x_{i}^{\prime}\not=x_{i}$, and
$\bar{\nabla}_{3}H(x,x^{\prime})_{i}=\frac{\partial H}{\partial
x_{i}}(x^{\prime}_{1},\ldots,x^{\prime}_{i-1},x^{\prime}_{i}=x_{i},x_{i+1},\ldots,x_{n})$
otherwise.

###  2.4 Irreversible port-Hamiltonian systems and discrete gradient methods

Assume for simplicity that $J_{i},J_{port,j}$, $M_{\alpha}$ and
$M_{port,\beta}$ are constant skew symmetric matrices. Choose a discrete
gradient $\bar{\nabla}$. Then the discretization of the IPHS is given by

$\displaystyle\frac{x_{k+1}-x_{k}}{h}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}\hat{\gamma_{i}}({x_{k+1/2}})\\{S,H\\}^{d}_{J_{i}}J_{i}\bar{\nabla}H(x_{k},x_{k+1})+\sum_{\alpha=1}^{l}M_{\alpha}\bar{\nabla}H(x_{k},x_{k+1})$
$\displaystyle+\sum_{j=1}^{\bar{k}}\widehat{\gamma_{Port,j}}({x_{k+1/2}},u_{k})\left[\\{S_{tot},H_{tot}\\}^{d}_{J_{port,g_{j}}}\right]g_{j}u_{k}+\sum_{\beta=1}^{\bar{l}}g_{S,\beta}u_{k}$
$\displaystyle y_{k}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{\bar{k}}\widehat{\gamma_{Port,j}}({x_{k+1/2}},u_{k})\left[\\{S_{tot},H_{tot}\\}^{d}_{J_{port,g_{j}}}\right]g_{j}^{\top}\bar{\nabla}H(x_{k},x_{k+1})$
$\displaystyle+\sum_{\beta=1}^{\bar{l}}g_{S,\beta}^{\top}\bar{\nabla}H(x_{k},x_{k+1})$

where we denote by
$\\{S,H\\}^{d}_{J_{i}}(x_{k},x_{k+1})=\bar{\nabla}S(x_{k},x_{k+1})^{\top}J_{i}\bar{\nabla}H(x_{k},x_{k+1})$
and
$\\{S_{tot},H_{tot}\\}^{d}_{J_{port,g_{j}}}(x_{k},x_{k+1})=(\bar{\nabla}S(x_{k},x_{k+1}),\tau_{k})^{\top}J_{port,g_{j}}(\bar{\nabla}H(x_{k},x_{k+1}),u_{k})$.
We also use the notation $x_{k+1/2}=(x_{k}+x_{k+1})/2$. In this case, $S$ is a
linear function, as a consequence, $\bar{\nabla}S(x_{k},x_{k+1})=\nabla
S(x_{k})$ is in the kernel of $M_{\alpha}$ and
$\left\\{S_{\textrm{tot}},H_{\textrm{tot}}\right\\}^{d}_{J_{\textrm{port},g_{S\beta}}}=0$,
assuming we are in the conditions of Definition 1. Therefore, we obtain the
discrete thermodynamic balance equations

$H(x_{k+1})-H(x_{k})-hy_{k}^{\top}u_{k}=0\;,\qquad
S(x_{k+1})-S(x_{k})-hy_{k}^{\top}\tau_{k}\geq 0\;.$

###  2.5 Simulating the gas-piston system using Discrete gradient methods

As a toy example to demonstrate the properties of the discrete gradient based
integrator, we have considered the gas-piston system with the expression of
the internal energy, for a closed system (constant number of moles $N_{0}$)

$\frac{U}{U_{\textrm{ref}}}=\exp\left(\frac{1}{c}\frac{\left(S-S_{\textrm{ref}}\right)}{N_{0}R}\right)\left(\frac{V_{\textrm{ref}}}{V}\right)^{\frac{1}{c}}$

where $R$ is the universal gas constant and $c$ is a positive number depending
on the nature of the gas (for a mono-atomic perfect gas $c=\frac{3}{2}$ for
instance) [5, page 68]. In this toy example we have chosen the reference
values $U_{\textrm{ref}}=1$, $V_{\textrm{ref}}=1$ and $S_{\textrm{ref}}=0$.
The heat diffusion coefficient is taken $\lambda_{e}=1$, the gravity constant
$g=10$ and the dissipation coefficient $\mu=0.5$. The initial state of the
system is chosen $S_{0}=2,V_{0}=4,q_{0}=1,p_{0}=0$. The controls are the
thermostat temperature $u_{1}=10$ and the external mechanical force
$u_{2}=-10$.

On the figure 1(c), one recovers that the numerical scheme preserves the power
balance equations and on the figure 1(b) that the internal energy production
term is indeed positive. In the figure 1(a) the evolution of the volume of the
gas and the displacement are shown and one sees that they are proportional and
can clearly distinguish two time interval, the first one, where the mechanical
oscillations is dominant, the gas acting as a nonlinear spring and the second
one, corresponding to the thermal relaxation to the equilibrium. In the figure
1(d), the time evolution of the temperature is shown and the two time
intervals appear clearly: the first one, corresponding to the damped
oscillatory behaviour where mechanical energy due to the external force is
converted to internal energy of the gas and the second one where the thermal
relaxation is dominant.

(a) Height and volume vs. time

(b) Positive Entropy balance in time.

(c) Energy balance vs time

(d) Piston temperature evolution

Figure 2.1: Simulation of the gas-piston system with controls.

##  3 Conclusions

In this paper, we have suggested a discretization scheme adapted for
irreversible port-Hamiltonian systems that allows to verify both the energy
balance equation and the positivity of entropy creation terms in the entropy
balance equation. In future work, we intend to extend our scheme to the
conservation of other geometric properties of this class of system to more
general classes of systems, e.g. systems with constraints or implicit systems
and apply this scheme to more realistic physical systems.

## Acknowledgements

A. Anahory Simoes and D. Martín de Diego acknowledges financial support from
the Spanish Ministry of Science and Innovation, under grants
PID2019-106715GB-C21 and CEX2019-000904-S.

## References

  * [1] Héctor Ramírez, Yann Le Gorrec, Bernhard Maschke, and Françoise Couenne.  Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the cstr.  Chemical Engineering Science, 89:223–234, 2013.
  * [2] Hector Ramirez, Bernhard Maschke, and Daniel Sbarbaro.  Modelling and control of multi-energy systems: an irreversible port-Hamiltonian approach.  Eur. J. Control, 19(6):513–520, 2013.
  * [3] Hector Ramírez.  Control of irreversible thermodynamical processes using port-Hamiltonian systems defined on pseudo-poisson and contact structures.  PhD thesis, Université Claude Bernhard (Lyon 1), 2012.
  * [4] V. Duindam, A. Macchelli, S. Stramigioli, and H. eds. Bruyninckx.  Modeling and Control of Complex Physical Systems - The Port-Hamiltonian Approach.  Springer, Sept. 2009.  ISBN 978-3-642-03195-3.
  * [5] H.B. Callen.  Introduction to Thermodynamics and Thermostatistics, 1985.  John Wiley and Sons, second edition.
  * [6] Philip J. Morrison.  A paradigm for joined Hamiltonian and dissipative systems.  volume 18, pages 410–419. 1986.  Solitons and coherent structures (Santa Barbara, Calif., 1985).
  * [7] Miroslav Grmela and Hans Christian Öttinger.  Dynamics and thermodynamics of complex fluids. I. Development of a general formalism.  Phys. Rev. E (3), 56(6):6620–6632, 1997.
  * [8] Bernhard Maschke and Jonas Kirchhoff.  Port maps of irreversible port hamiltonian systems.  https://arxiv.org/abs/2302.09026, 2023.
  * [9] Alexandre Anahory Simoes, Manuel de León, Manuel Lainz Valcázar, and David Martín de Diego.  Contact geometry for simple thermodynamical systems with friction.  Proc. A., 476(2241):20200244, 16, 2020.
  * [10] Hans Christian Öttinger.  GENERIC integrators: structure preserving time integration for thermodynamic systems.  J. of Non-equilibrium Thermodynamics, 43:89–100, 2018.
  * [11] A. Falaize and Th. Hélie.  Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach.  Applied Sciences, 6(10):273, 2016.
  * [12] Paul Kotyczka and Laurent Lefevre.  Discrete-time port-hamiltonian systems: A definition based on symplectic integration.  Systems & Control Letters, 133:104530, 2019.
  * [13] O. Gonzalez.  Time integration and discrete Hamiltonian systems.  J. Nonlinear Sci., 6(5):449–467, 1996.
  * [14] T. Itoh and K. Abe.  Hamiltonian-conserving discrete canonical equations based on variational difference quotients.  J. Comput. Phys., 76(1):85–102, 1988.
  * [15] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux.  Geometric integration using discrete gradients.  R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci., 357(1754):1021–1045, 1999.
  * [16] A. Harten, P. D. Lax, and B. van Leer.  On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.  SIAM Rev., 25(1):35–61, 1983.


# Long-range entanglement in the XXZ Heisenberg spin chain after a local
quench

Jie Ren1111E-mail<EMAIL_ADDRESS>Shiqun Zhu2222E-mail<EMAIL_ADDRESS>1Department of Physics and Jiangsu Laboratory of Advanced Functional
Materials, Changshu Institute of Technology, Changshu, Jiangsu 215500, China
2School of Physical Science and Technology, Suzhou University, Suzhou, Jiangsu
215006, China

###### Abstract

The long-range entanglement dynamics of an one-dimensional spin-1/2
anisotropic XXZ model are studied using the method of the adaptive time-
dependent density-matrix renormalization-group. The long-range entanglement
can be generated when a local quench on one of boundary bonds is performed in
the system. The anisotropic interaction has a strong influence both on the
maximal value of long-range entanglement and the time of reaching the maximum
long-range entanglement. The local coupling has a notable impact on the long-
range entanglement, but it can be neglected in the time of reaching the
maximal long-range entanglement.

PACS number: 03.67.Mn, 03.65.Ud, 75.10.Pq

Keywords: long-range entanglement, spin chain, adaptive time-dependent
density-matrix renormalization-group

##  I Introduction

Entanglement generation and distribution is one of the important problems in
performing quantum-information tasks, such as quantum computation and quantum
teleportation Nielsen ; Bennett ; Wootters . Many results showed that
entanglement existed naturally in the spin chain when the temperature is at
zero Osborne ; Osterloh . It is discouraging that the entanglement in many
spin systems is typically very short ranged. It exists only in the nearest
neighbors and the next nearest neighbors Amico . It is interesting that some
schemes for long-range entanglement, such as exploiting weak couplings of two
distant spins to a spin chain were proposed Venuti01 ; Venuti02 ; Giampaolo .
These methods have limited thermal stability or very long time scale of
entanglement generation. In recent years, many researches has shifted the
focus of the study in the dynamics of entanglement Bose0 ; Hartmann ; Li ; A ;
Bose01 . The dynamics of long range entanglement are obtained by a time quench
of magnetic field Galve ; Wang . The end-to-end entanglement can be shared
across a chain of arbitrary range Bose0 . It is excited that the most common
situation studies so far concerns a sudden quench of some coupling of the
model Hamiltonian Chiara ; Bose02 ; Bose03 . In Ref. Bose03 , it showed that
the long-range entanglement can be engineered by a non-perturbative quenching
of a single bond in a Kondo spin chain with impurity. This is the first
example that a minimal local action on a spin chain can generate long-range
entanglement dynamically. It would be interesting to investigate the
possibility of producing long-range entanglement in the XXZ Heisenberg spin
chain without impurity.

It is well-known that the Hamiltonian of an opened chain of N spin-$1/2$
system with nearest-neighbor XXZ interaction is given by

$H_{0}=\sum_{i=1}^{N-1}J[S^{x}_{i}S^{x}_{i+1}+S^{y}_{i}S^{y}_{i+1}+\Delta
S^{z}_{i}S^{z}_{i+1}],\\\ $ (1)

where $S^{\alpha}_{i}(\alpha=x,y,z)$ are spin operators on the $i$-th site,
$N$ is the length of the spin chain. The parameter $\Delta$ denotes the
couplings in the $z$-axis.This model can be realized in the Josephson-
junctionFazio and optical latticesDuan ; Trotzky .

In the paper, the long-range entanglement in the XXZ Heisenberg spin chain
after a local quench is investigated. An opened boundary condition (OBC) is
assumed because the antiferromagnetic Heisenberg spin chain with OBC can be
achieved artificially in the experiment hirjibehedin , and the coupling $J=1$
is considered for simplicity. In section II, the local quench is presented.
The concurrence is used as a measurement of the entanglement. In section III,
a more general situation and a single quench and its effect on the end-to-end
qubits entanglement are analyzed. The effects of the anisotropic interaction
and system size on the end-to-end qubits entanglement is also studied. The
robustness of the entanglement against an increase in temperature is
investigated. A discussion concludes the paper.

##  II Local Quench and Entanglement Measure

In the paper, the system is assumed to be in the ground state $Gs_{0}$ of
$H_{0}$ initially. A local quench change of the coupling between the first
qubit and the second qubit with the same anisotropy interaction $\Delta$ in
Eq. (1) is performed. The Hamiltonian of the system modifies to

$H_{1}=J_{1}[S^{x}_{1}S^{x}_{2}+S^{y}_{1}S^{y}_{2}+\Delta
S^{z}_{1}S^{z}_{2}]+\sum_{i=2}^{N-1}J[S^{x}_{i}S^{x}_{i+1}+S^{y}_{i}S^{y}_{i+1}+\Delta
S^{z}_{i}S^{z}_{i+1}].\\\ $ (2)

Since $[H_{0},H_{1}]\neq 0$, the the ground state $Gs_{0}$ of $H_{0}$ is not
one of the eigenstates of $H_{1}$. The state of the system will evolve as

$\psi(t)=\exp^{-iH_{1}t}Gs_{0}.$ (3)

In the paper, the concurrence is chosen as a measurement of the pairwise
entanglement Wootters . The concurrence $C$ is defined as

$C_{1,N}=\max\\{{\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4},0}\\},$ (4)

where the quantities $\lambda_{i}(i=1,2,3,4)$ are the square roots of the
eigenvalues of the operator
$\varrho=\rho_{12}(\sigma_{1}^{y}\otimes\sigma_{N}^{y})\rho_{1,N}^{\ast}(\sigma_{1}^{y}\otimes\sigma_{N}^{y})$.
They are in descending order. The case of $C_{1,N}=1$ corresponds to the
maximum entanglement between the two qubits, while $C_{1,N}=0$ means that
there is no entanglement between the two qubits.

##  III Long-Range Entanglement Dynamics

It is known that it is hard to calculate the dynamics of entanglement because
of the lack of knowledge of eigenvalues and eigenvectors of the Hamiltonian.
For models that are not exactly solvable, most of researchers resort to exact
diagonalization to obtain the ground state for small system size. The method
is difficult to be applied for large system size $N>20$. For a large system,
the adaptive time-dependent density-matrix renormalization-group can be
applied with a second order Trotter expansion of the Hamiltonian as described
in White ; Vidal02 . In order to check the accuracy of the results of the
adaptive time-dependent density-matrix renormalization-group, the results of
exact diagonalization can be considered as a benchmark for a small size
system.

The numerical error of the adaptive time-dependent density-matrix
renormalization-group comes from the discarded weight and the Trotter
decomposition. The error of discarded weight is dependent on the data
precision and truncated Hilbert space. The error of the Trotter decomposition
is relied on Trotter slicing. In our numerical simulations a Trotter slicing
$\delta t=0.05$ and Matlab codes of the adaptive time-dependent density-matrix
renormalization-group with double precision are performed with a truncated
Hilbert space of $m=100$. In turns out that a typically discarded weight of
$\delta\rho\leq 10^{-8}$. The error of Trotter decomposition
$\delta\propto(\delta t)^{3}$. These can keep the relative error $\delta C$ in
$C$ below $10^{-3}$ for a chain of $N=60$ sites with time $t\leq 50/J$.

In the adaptive time-dependent density-matrix renormalization-group, it is
hard to calculate the reduced density matrix $\rho_{1,N}$. By making use of
the relation between the correlation and the reduced density matrix, the
reduced density matrix can be expressed as

$\rho_{1,N}=\frac{1}{4}[I_{1,N}+\sum_{i,j=x,y,z}\langle\sigma_{1}^{i}\sigma_{N}^{j}\rangle\sigma_{1}^{i}\sigma_{N}^{j}],$
(5)

where $\sigma^{\alpha}_{k}(\alpha=x,y,z)$ are Pauli operators and $I_{1,N}$ is
identity matrix.

We obtain the reduced density matrix $\rho_{1,N}$ by Eq. (5), then calculate
the entanglement between two-ends qubits. The end-to-end entanglement
$C_{1,N}$ is plotted as a function of the anisotropic interaction $\Delta$ and
time $t$ with $J_{1}=-0.01$ in Fig. 1(a). It is seen that the long-range
entanglement $C_{1,N}$ can be generated after a short time. When the time $t$
increases, the entanglement reaches the maximal value and then disappears
quickly. That is, there is a peak in the entanglement. The peak of the end-to-
end entanglement is plotted as a function of the anisotropic interaction
$\Delta$ in Fig. 1(b). It is seen that the peak in long-range end-to-end
entanglement increases with the anisotropic interaction $\Delta$. It reaches
the maximal value when the anisotropic interaction $\Delta=1$ Gu ; Bose04 .
With the anisotropic interaction $\Delta$ increases further, the height of the
peak decreases.

The relation between the time of reaching the maximal end-to-end entanglement
labeled by $T_{max}$ and the anisotropic interaction $\Delta$ can be seen in
the inset of Fig. 1(b). The time of reaching the maximal long-range
entanglement decreases when the anisotropic interaction $\Delta$ increases,
while they are also not a linear relationship. Similar to Refs. Bose0 ; Bose03
, the system generates end-to-end entanglement periodically. In our
simulations, the error of the simulations increases when the time increases.
Since the coherent time of the system is not very long, the long-distance
entanglement is not shown when the time $t>50/J$. It is shown that a sharp
drop of long-distance entanglement occurs at $\Delta=-0.5$. Furthermore, the
time of the long-range entanglement reaching the maximal value goes up
drastically. This is similar to the anomalous behavior appeared in Ref. Bose04
.

The end-to-end entanglement $C_{1,N}$ is plotted as a function of the
interaction $J_{1}$ and the time $t$ with $\Delta=1$ in Fig. (2). The size of
the system $N=20$ is chosen. It does not include the case of $J_{1}=0$. The
long-range entanglement $C_{1,N}$ can be generated after a short time. There
is a peak in $C_{1,N}$ when $C_{1,N}$ is plotted as a function of time $t$. It
is seen that the influence of the coupling $J_{1}$ on the maximal long-range
entanglement is relatively large. The peak of the end-to-end entanglement
increases when the coupling interaction $J_{1}$ increases, and reaches the
maximal value when $J_{1}=-0.1$. It seems that the changing of local
interaction from antiferromagnetic to ferromagnetic may enhance the
entanglement creation Venuti01 . When $J_{1}$ increases further to
$J_{1}>0.32$, the long-range entanglement disappears. It is interesting that
the coupling interaction $J_{1}$ has a relatively small impact on the time
when the end-to-end entanglement generates and reaches the maximal value.

It is easy to obtain that
$[H_{0},\sum_{i=1}^{N}S^{z}_{i}]=[H_{1},\sum_{i=1}^{N}S^{z}_{i}]=0$. It means
that the ground states of $H_{0}$ is a total singlet $S_{tot}=0$. The
Hamiltonian $H_{1}$ has invariant on every excitation subspaces. During the
evolution, the boundary spin at $i=2$ will have a strong tendency to form a
singlet pair with its only nearest neighbor $i=3$ on the right-hand side. This
is similar for spin pairs $(4,5)$ and $(6,7)$, ……, etc. Then the two-end spin
qubits form a singlet. The local interaction $J_{1}$ decides the ability of
forming singlet for even bands. Thus, the long-distance entanglement creates.

The thermalization and relaxation during the period of generating long-
distance entanglement are neglected because the dynamical time scale is quite
short. When initial state is taken to be the relevant thermal state, the end-
to-end entanglement is plotted as a function of temperature in Fig. 3 after
bond quenching for different system sizes. It is shown that the entanglement
vanishes when $kT>1.16$ with the size of $N=8$, and $kT>1.06$ with the size of
$N=10$. It seems that our scheme is quite robust against temperature and is
similar to that in Ref. Bose03 . The initial thermal state does not have long-
distance entanglement in Eq. (1), but the energy gap between the ground state
and the first and other excited states are larger than the system with
impurity. This leads to our scheme is quite robust against temperature.

It is interesting to investigate the long-range entanglement creation even for
very long chain of large size. The maximal value of the end-to-end
entanglement $C_{1,N}$ is plotted as a function of the size of the system for
different anisotropic interaction $\Delta$ and different interaction $J_{1}$
in Fig. 4(a). It is found that the maximal value of the end-to-end
entanglement $C_{1,N}$ decreases when the size of system increases. When the
anisotropic interaction $\Delta=1$ and the interaction $J_{1}=-0.1$, the
maximal value is 0.5481 for $N=40$ and 0.4547 for $N=60$. It is shown that the
maximal value of the end-to-end entanglement label by $\xi(N)$ is given by
Bose0

$\xi(N)\simeq 1.35N^{-1/3}.$ (6)

In Ref. Bose0 , $\xi(40)\simeq 0.3947$ for $N=40$, and $\xi(40)\simeq 0.3448$
for $N=60$. While in the Kondo regime, $\xi(40)$ can be as large as $0.7$ for
$N=40$. It seems that our results are also quite good. The long distance
entanglement can be obtained by performing a local quench in either one of two
boundary bands. The time of reaching the maximal end-to-end entanglement
labeled by $T_{max}$ is plotted as a function of the size $N$ of the system
for different anisotropic interaction $\Delta$ and different interaction
$J_{1}$ in Fig. 4(b). The time of the end-to-end entanglement reaching the
maximal value $T_{max}$ is linear increase when the size $N$ increases Bose0 ;
Bose04 . It is seen that the slope of the line is almost dependent on the
anisotropic interaction $\Delta$.

In the paper, adiabatic quenches are studied and the decoherence effects in
the system are ignored. It is noted that the XXZ spin chains can be realized
in the Josephson-junction array Fazio and optical lattices Duan ; Trotzky . In
Josephson-junction array, the interactions and the anisotropy can be modulated
by varying voltages Allcock . The Josephson-junction array system does not
have significant decoherence over our time-scales ($T_{max}\simeq N/J$)
Lyakhov . The anisotropic interaction quenches can be achieved successfully
with long decoherence time in the optical lattices experiment Trotzky .

##  IV Discussion

By using the method of the adaptive time-dependent density-matrix
renormalization-group, the time evolution of the entanglement in a one-
dimensional spin-1/2 anisotropic XXZ model is investigated when a local quench
is performed in the system. The local quench is a abrupt change of interaction
between the first qubit and the second qubit. The dynamics of pairwise
entanglement between the two ends qubits in the spin chain is studied. The
entanglement of the two-ends spin qubits can be created after the local quench
is performed.The time when the long-range entanglement generates decreases in
a nonlinear manner with the anisotropic interaction increases. It reaches the
maximal value when the anisotropic interaction $\Delta=1$. The maximal value
of the end-to-end entanglement increases with the coupling interaction $J_{1}$
increases, and reaches the maximal when $J_{1}=-0.1$. It is interesting that
the coupling interaction has a relatively small impact on the time when the
long-range entanglement reaches the maximal value. This phenomenon may be used
to control the dynamics of the entanglement by varying the anisotropic
interaction and the local quench interaction of the Heisenberg spin chain.

Acknowledgments

It is a pleasure to thank Yinsheng Ling and Yinzhong Wu for their many helpful
discussions. The financial support from the National Natural Science
Foundation of China (Grant No. 10774108) is gratefully acknowledged.

## References

  * (1) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, England, 2000).
  * (2) C. H. Bennett, G. Brassard, C. Crepeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * (3) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * (4) T. J. Osborne and M. A. Nielsen, Phys. Rev. A 66, 032110 (2002).
  * (5) A. Osterloh, L. Amico, G. Falci, and R. Fazio, Nature 416, 608 (2002).
  * (6) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. 80, 517 (2008).
  * (7) L. C. Venuti, C. D. E. Boschi, and M. Roncaglia, Phys. Rev. Lett. 96, 247206 (2006).
  * (8) L. C. Venuti, S. M. Giampaolo, F. Illuminati, and P. Zanardi, Phys. Rev. A 76, 052328 (2007).
  * (9) S. M. Giampaolo and F. Illuminati, New J. Phys. 12, 025019 (2010).
  * (10) S. Bose, Phys. Rev. Lett. 91, 207901 (2003).
  * (11) M. J. Hartmann, M. E. Reuter, and M. B. Plenio, New J. Phys. 8, 94 (2006).
  * (12) Y. Li, T. Shi, B. Chen, Z. Song, and C.-P. Sun, Phys. Rev. A 71, 022301 (2005).
  * (13) A. Wójcik, T. Łuczak, P. Kurzyński, A. Grudka, T. Gdala, and M. Bednarska, Phys. Rev. A 72, 034303 (2005).
  * (14) J. Reslen and S. Bose, Phys. Rev. A 80, 012330 (2009).
  * (15) F. Galve, D. Zueco, S. Kohler, E. Lutz, and P. Hänggi, Phys. Rev. A 79, 032332 (2009).
  * (16) X. Wang, A. Bayat, S. G. Schirmer, S. Bose, Phys. Rev. A 81, 032312 (2010)
  * (17) G. D. Chiara, S. Montangero, P. Calabrese, and R. Fazio, J. Stat. Mech. P03001 (2006).
  * (18) H. Wichterich and S. Bose, Phys. Rev. A 79, 060302(R) (2009).
  * (19) P. Sodano, A. Bayat, and S. Bose, Phys. Rev. B 81, 100412(R) (2010).
  * (20) R. Fazio and H. van der Zant, Phys. Rep. 355, 235 (2001).
  * (21) L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
  * (22) S. Trotzky, P. Cheinet, S. Fölling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch, Science 319, 295 (2008).
  * (23) C. F. Hirjibedin, C. P. Lutz, and A. J. Heinrich, Science 312, 1021 (2006).
  * (24) S. R. White and A. E. Feiguin, Phys. Rev. Lett.  93, 076401 (2004).
  * (25) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * (26) S.-J. Gu, G.-S. Tian, and H.-Q. Lin, Phys. Rev. A 71, 052322 (2005)
  * (27) A. Bayat and S. Bose, Phys. Rev. A 81, 012304 (2010).
  * (28) J. Allcock and N. Linden, Phys. Rev. Lett. 102, 110501 (2009).
  * (29) A. Lyakhov and C. Bruder, New J. Phys. 7, 181 (2005).

Figure 1: (a). The long-range entanglement $C_{1,N}$ is plotted as a function
of anisotropic interaction $\Delta$ and time $t$ with $J_{1}=-0.1$. The size
of the system is $N=20$. (b). The maximal value of the end-to-end entanglement
is plotted as a function of the anisotropic interaction $\Delta$. The inset
shows the time of reaching the maximal end-to-end entanglement labeled by
$T_{max}$ as a function of anisotropic interaction $\Delta$.

Figure 2: The long-range entanglement $C_{1,N}$ is plotted as a function of
interaction $J_{1}$ and time $t$ when there is local quench of the interaction
$J_{1}$. The size of the system is $N=20$ and the anisotropic interaction
$\Delta=1$. It does not include the case of $J_{1}=0$.

Figure 3: The maximal value of end-to-end entanglement is plotted as a
function of temperature after band quenching for different sizes.

Figure 4:  (a). The maximal value of the long-range entanglement $C_{1,N}$ is
plotted as a function of the size of the system for different anisotropic
interaction $\Delta$ and different interaction $J_{1}$. (b). The time of
reaching the maximal long-range entanglement labeled by $T_{max}$ is plotted
as a function of the size of the system for different anisotropic interaction
$\Delta$ and different interaction $J_{1}$. The red lines are fixed lines.


# Smart Pixels: In-pixel AI for on-sensor data filtering

Benjamin Parpillon1,2, Chinar Syal1, Jieun Yoo2, Jennet Dickinson1, Morris
Swartz3, Giuseppe Di Guglielmo1,4,
Alice Bean5, Douglas Berry1, Manuel Blanco Valentin4, Karri DiPetrillo6,
Anthony Badea6, Lindsey Gray2,
Petar Maksimovic3, Corrinne Mills2, Mark S. Neubauer8, Gauri Pradhan1, Nhan
Tran1,4, Dahai Wen3, Farah Fahim1
1Fermi National Accelerator Laboratory (FNAL), 2University of Illinois
Chicago, 3 Johns Hopkins University, 4 Northwestern University,
5 University of Kansas, 6 The University of Chicago, 8 University of Illinois
Urbana-Champaign, Champaign

###### Abstract

We present a smart pixel prototype readout integrated circuit (ROIC) designed
in CMOS 28 nm bulk process, with in-pixel implementation of an artificial
intelligence (AI) / machine learning (ML) based data filtering algorithm
designed as proof-of-principle for a Phase III upgrade at the Large Hadron
Collider (LHC) pixel detector. The first version of the ROIC consists of two
matrices of 256 smart pixels, each 25$\times$25 µm2 in size. Each pixel
consists of a charge-sensitive preamplifier with leakage current compensation
and three auto-zero comparators for a 2-bit flash-type ADC. The frontend is
capable of synchronously digitizing the sensor charge within 25 ns.
Measurement results show an equivalent noise charge (ENC) of $\sim$30e- and a
total dispersion of $\sim$100e- The second version of the ROIC uses a fully
connected two-layer neural network (NN) to process information from a cluster
of 256 pixels to determine if the pattern corresponds to highly desirable
high-momentum particle tracks for selection and readout. The digital NN is
embedded in-between analog signal processing regions of the 256 pixels without
increasing the pixel size and is implemented as fully combinatorial digital
logic to minimize power consumption and eliminate clock distribution, and is
active only in the presence of an input signal. The total power consumption of
the neural network is $\sim$ 300 $\mu$W. The NN performs momentum
classification based on the generated cluster patterns and even with a modest
momentum threshold, it is capable of 54.4% – 75.4% total data rejection,
opening the possibility of using the pixel information at 40MHz for the
trigger. The total power consumption of analog and digital functions per pixel
is $\sim$ 6 $\mu$W per pixel, which corresponds to $\sim$ 1 W/cm2 staying
within the experimental constraints.

##  I Introduction

Detectors at forthcoming high-energy colliders will confront substantial
technical hurdles. Managing the unprecedented volume of particles anticipated
in each event demands highly detailed silicon pixel detectors with billions of
readout channels. With event frequencies reaching 40 MHz, these detectors will
produce immense amounts of data every second. To facilitate discoveries while
adhering to stringent bandwidth and latency constraints, future trackers must
possess the ability to swiftly, efficiently, and robustly reduce data at the
source while operating in high radiation environments. This prototype
development targets Phase III pixel detector upgrades projected for 2034, as
well as other future colliders.

##  II Analog Frontend for Hit Detection and Digitization with Test results

Reducing the pixel sizes by a factor of four while maintaining the same sensor
thickness leads to charge sharing which can be positively exploited to improve
the effective position resolution [1]. However, this requires reducing the
minimum detection threshold of the frontend and overcoming pileup errors.

Fig. 1: The pixel analog frontend for hit registration in a particle detector
containing a preamplifier with leakage current compensation and AC coupled
synchronous comparators.

###  II-A Architecture

The pixel architecture for in-pixel detection and digitization of collision
events is depicted in Fig. 1.

The charge collected at the sensor’s electrode is integrated, amplified, and
converted to voltage using a charge-sensitive preamplifier. An AC-coupled
2-bit flash-type ADC digitizes the signal. Due to the thermometric nature of
the flash ADC in our design, analog-to-digital conversion begins as soon as
the integrated charge output is above the first threshold, and continues until
the signal reaches its maximum value or the time for conversion runs out
(since all comparators are reset at the end of the bunch crossing cycle). The
first ROIC prototype was designed and fabricated in 2023. The topology,
implementation and simulation results are discussed in detail in [2].

Fig. 2: a. Picture of the first 1.5mm2 ROIC chip. The chip consists of two 16
× 16-pixel matrix variants with just the analog frontend. b. Layout of the
second prototype chip. The chip consists of two 8 × 32 - pixel matrix variants
with the analog frontend and the digital NN

###  II-B Test Results Of First ROIC Prototype

Fig. 2a shows the first ASIC chip, $1.5\text{\,}\mathrm{m}\mathrm{m}^{2}$ in
size. It is composed of a 512 pixels matrix divided into two halves with a
single-ended and differential frontend variants that were implemented. The
ROIC was characterized without connecting any sensor and all the comparators
in the pixel were set at 400e- equivalent threshold voltages. The ASIC
incorporates programmable front-end charge injection circuitry to generate
pixel cluster charges during characterization (bottom left in Fig. 1). The
programming didn’t work as reliably as expected which added a stochastic error
to our results. The total equivalent noise charge (ENC) of the preamplifier
and the comparator is presented in Fig. 3. The histograms were extracted by
scanning the entire input dynamic range in steps of 2 e- and injecting charges
200 times for each step at every pixel. After removing the incorrectly
programmed pixels, we find the ENC to be 32e- for both variants. The total
dispersion for each matrix variant is presented in Fig. 4. For both variants,
the effective threshold dispersion of the hit comparator across the entire
ASIC is 100 e-. This result includes gain dispersion in the preamplifier and
threshold dispersion of the comparators.

Fig. 3: Total pixel ENC for both variant with 400e- equivalent threshold
voltages and no sensor capacitance connected to the ROIC. a. ENC histogram of
the single ended variant. b. ENC histogram of the differential variant Fig. 4:
S-curves of total dispersion for all Hit Comparators in the Array for both
variants with 400e- equivalent threshold voltages and no sensor capacitance
connected to the ROIC

##  III On-chip Neural Network Implementation

A neural network classifier was designed and implemented on the second ROIC
prototype (Fig. 2b) to identify clusters associated with high $p_{T}$ charged
particles. The high $p_{T}$ signal class contains tracks with $p_{T}>200$ MeV.
Oppositely charged particle tracks curve in opposite directions in the
magnetic field, resulting in clusters of very different shape, therefore two
background classes are defined, corresponding to positively and negatively
charged particles with $p_{T}<200$ MeV. The choice of $p_{T}$ threshold in the
training class definition affects the physics performance and could be
adjusted depending on the physics goals. The details about the training
datasets and the algorithm development can be found in this reference [3].

We adopted hls4ml to translate the neural network classifier into optimal
hardware implementations [4]. hls4ml is an open-source Python framework that
facilitates the co-design of machine learning algorithms for hardware
deployment, supporting models from quantized models from QKeras and other
formats [5]. hls4ml allowed us to fine-tune the numerical precision and the
hardware parallelism to optimize area, performance, and power consumption
according to our system constraints. The conversion process started with the
quantized model of the classifier, which hls4ml translated into HLS-ready C++
code for Siemens Catapult HLS [6] that generates a hardware description at the
register-transfer level (RTL) suitable for the ASIC flow. We chose to fully
parallelize the hardware logic to minimize the latency of the neural network,
integrating the HLS-generated RTL design with system registers and data movers
for efficient operation shown in Fig. 5.

Fig. 5: Data flow through the digital implementation of the algorithm from the
summed ADC bits (on the left) through the neural network and the final
classification layer. At the top of the diagram we illustrate the
reconfigurability of the weights and biases in the algorithm stored in memory.

##  IV Conclusions and Future Works

We have presented the first prototype test results and showed the analog
frontend based on a synchronous ADC architecture is capable of processing and
digitizing input signals within a 25 ns bunch crossing period. The compact
area of the analog part of the pixel (<210µm2) results in $\times$4 improved
granularity while offering a $\times$2.5 power reduction ($\sim$4 µW/pixel)
and an $\times$2.5 improvement in the minimum threshold detection ($\approx$
475e-) compared to state of the art RD53B designs [7] such as CROCv1[8] and
ITKv1[9]. We also showed the implementation of machine learning-based
approaches for data processing in the pixelated region, such as analyzing
cluster shapes to determine track parameters and/or suppressing noise hits. By
leveraging the proposed architecture, only useful clusters will be transmitted
to the periphery for further processing and off-chip data transfer, thus
minimizing the necessary data bandwidth for data collection under much higher
collision rates for future upgrades. The second prototype chip is expected to
be received in May 2024 and test results will be presented at the conference.

## References

  * [1] E. A. A. Bergamaschi, M. Andrä1, “The mÖnch detector for soft x-ray, high-resolution, and energy resolved applications,” _Synchrotron Radiation News_ , vol. 36, no. 6, pp. 11–15, 2018.
  * [2] B. Parpillon, A. Trivedi, and F. Fahim, “Readout IC with 40 MSPS in-pixel ADC for future vertex detector upgrades of Large Hadron Collider,” _2023 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2023. [Online]. Available: https://ieeexplore.ieee.org/document/10182033
  * [3] J. Yoo, J. Dickinson, M. Swartz, G. Di Guglielmo, A. Bean, D. Berry, M. B. Valentin, K. DiPetrillo, F. Fahim, L. Gray _et al._ , “Smart pixel sensors: towards on-sensor filtering of pixel clusters with deep learning,” _arXiv preprint arXiv:2310.02474_ , 2023.
  * [4] F. Fahim, B. Hawks, C. Herwig, J. Hirschauer, S. Jindariani, N. Tran, L. P. Carloni, G. D. Guglielmo, P. Harris, J. Krupa, D. Rankin, M. B. Valentin, J. Hester, Y. Luo, J. Mamish, S. Orgrenci-Memik, T. Aarrestad, H. Javed, V. Loncar, M. Pierini, A. A. Pol, S. Summers, J. Duarte, S. Hauck, S.-C. Hsu, J. Ngadiuba, M. Liu, D. Hoang, E. Kreinar, and Z. Wu, “hls4ml: An open-source codesign workflow to empower scientific low-power machine learning devices,” 2021.
  * [5] C. N. Coelho, A. Kuusela, S. Li, H. Zhuang, J. Ngadiuba, T. K. Aarrestad, V. Loncar, M. Pierini, A. A. Pol, and S. Summers, “Automatic heterogeneous quantization of deep neural networks for low-latency inference on the edge for particle detectors,” _Nature Machine Intelligence_ , vol. 3, no. 8, pp. 675–686, 2021.
  * [6] Siemens, “Catapult HLS,” https://eda.sw.siemens.com/en-US/ic/ic-design/high-level-synthesis-and-verification-platfor.
  * [7] L. Gaioni, “Cms analog front-end: simulations and measurements,” https://cds.cern.ch/record/2746420/files/Linear%20AFE%20report.pdf.
  * [8] A. Papadopoulos and the CMS Tracker Group, “Analog performance of the crocv1 pixel readout chip for the cms phase-2 tracker upgrade,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 18, 2023.
  * [9] E. A. Stefano Terzo, Maurizio Boscardin, “Novel 3d pixel sensors for the upgrade of the atlas inner tracker,” _Frontiers in Physics_ , vol. 9, 2021.