# Efficient Industrial Refrigeration Scheduling with Peak Pricing

Rohit Konda, Jordan Prescott, Vikas Chandan, Jesse Crossno,
Blake Pollard, Dan Walsh, Rick Bohonek, and Jason R. Marden R. Konda
<EMAIL_ADDRESS>and J. R. Marden (jrmarden@ece.ucsb.edu)are with the
Department of Electrical and Computer Engineering at the University of
California, Santa Barbara, CA. Vikas Chandan<EMAIL_ADDRESS>Jesse
Crossno<EMAIL_ADDRESS>Blake Pollard<EMAIL_ADDRESS>and Dan
Walsh<EMAIL_ADDRESS>are with CrossnoKaye®. Rick Bohonek
<EMAIL_ADDRESS>is with Butterball LLC®. This work is supported by
funding from CrossnoKaye.

###### Abstract

The widespread use of industrial refrigeration systems across various sectors
contribute significantly to global energy consumption, highlighting
substantial opportunities for energy conservation through intelligent control
design. As such, this work focuses on control algorithm design in industrial
refrigeration that minimize operational costs and provide efficient heat
extraction. By adopting tools from inventory control, we characterize the
structure of these optimal control policies, exploring the impact of different
energy cost-rate structures such as time-of-use (TOU) pricing and peak
pricing. While classical threshold policies are optimal under TOU costs,
introducing peak pricing challenges their optimality, emphasizing the need for
carefully designed control strategies in the presence of significant peak
costs. We provide theoretical findings and simulation studies on this
phenomenon, offering insights for more efficient industrial refrigeration
management.

##  I Introduction

Industrial refrigeration systems are widely utilized across diverse sectors,
not limited to plastics manufacturing, chemical processing, food storage, and
electronics production [1, 2, 3, 4]. Collectively, industrial refrigeration
contributes to approximately $8.4\%$ of total energy consumption in the United
States [5]. Consequently, there exist significant opportunities for energy
conservation within industrial refrigeration. These opportunities extend
beyond hardware upgrades to include enhancements in the control algorithms
implemented in these systems. Improvement in the algorithm design can
potentially be more appealing, as they can yield substantial energy savings
with minimal capital investment required for system retrofitting.

The intention behind improvements in algorithm design is to strategically
adjust the control strategies of the components within the industrial
refrigeration process to achieve the necessary heat extraction while
minimizing operation costs. These costs can be measured through total power,
electric costs, carbon emissions, or other relevant metrics. For example, the
set points or steady state configurations of the components can be optimized
to raise the energy efficiency [6, 7, 8, 9]. Furthermore, model predictive
control or trajectory optimization can be utilized to dynamically optimize the
energy efficiency of the refrigeration cycle over a time horizon [10, 11, 12].
Another approach that has garnered significant attention is _thermal load
shifting_ , where refrigeration loads are dynamically managed to leverage
variable energy cost-rate structures [13, 14, 15, 16, 17]. While these results
are quite encouraging in conventional cost structures, there remains a need to
characterize the qualitative behavior of optimal control policies when
considering more varied cost structures.

There exist many different energy cost-rate structures that may depend on
various factors. _Fixed rate_ pricing, the most traditional rate structure,
involves charging a flat rate per unit of energy consumed, irrespective of
external conditions. In contrast, _time-of-use_ (TOU) pricing varies the per-
unit rate based on the time or season. Typically, this results in higher
energy costs during peak demand periods and lower energy costs during off-peak
periods. In this setup, there may be significant economic benefit for the
operator to shift their energy usage to off-peak hours. Deciding on how to
shift is exactly the focus of the thermal load shifting literature. To further
regulate energy usage, _peak pricing_ can also be introduced, where the cost
is dependent on the maximum power usage of the system over a time period and
spikes in energy consumption are highly disincentived. In many industrial
refrigeration systems, peak pricing may comprise of a large portion of the
energy costs. In this work, we focus on when both peak pricing and TOU costs
are present and quantify the impact on optimizing the scheduling of
refrigeration loads.

Figure 1: We depict a prototypical power consumption profile over different
pricing regions.

Case Study. We introduce an example of a rate structure that a refrigeration
facility may be charged with for its energy consumption.111This data was
obtained through a direct collaboration with CrossnoKaye (see
crossnokaye.com), a company that focuses in integrating control systems for
industrial refrigeration systems within the cold food and beverage sector. The
time-of-use costs are dependent on time of day, where the cost rates vary
between three time windows: on-peak hours (4:00PM - 9:00PM), mid-peak hours
(9:00PM - 8:00AM), and off-peak hours (8:00AM - 4:00PM). Typically, energy-
costs are higher in on-peak hours and lower in off-peak hours. Additionally,
the facility can incur additional peak costs: a peak charge over the maximum
energy consumption over a month and a peak charge over the month for each of
the specific time windows (on, mid, off). We display a possible demand profile
under this cost structure in Figure 1.

Designing control policies using optimization methods that account for peak
pricing, though not as common as TOU pricing, has been studied previously [18,
19, 20]. However, the primary objective of this work is to theoretically
evaluate how peak pricing influences the qualitative structure of the optimal
control policies. We approach this through uniquely adopting the perspective
of _inventory control_. Inventory control [21, 22] is a classical branch of
multi-stage decision problems that explores purchasing policies of inventory
to ensure optimal warehouse stock levels. When energy rates solely consist of
time-of-use (TOU) costs, classical findings from inventory control suggest
that the optimal control algorithms should adopt a threshold approach: when
the facility’s temperature exceeds a specified threshold, the refrigeration
load is increased to return the facility to a desired buffer temperature.
Threshold policies are commonly implemented in practice for industrial
refrigeration; we verify its optimality under TOU costs in Proposition 1.
However, the introduction of peak pricing disrupts the optimality of threshold
policies. In fact, in Proposition 3, we characterize the structure of the
optimal control policy under peak pricing; we see that the optimality of
simple threshold policies is lost in these settings. We also verify our
findings through simulation in Section IV. Our work suggests that designing
control policies should be done carefully if significant peak costs are
present.

##  II Peak Pricing Model

In this paper, we aim to devise scheduling strategies for refrigeration
systems that strike a balance between minimizing overall costs and adequately
meeting the necessary cooling demands of the facility. Additionally, we
investigate the impact of _peak pricing_ , a common energy pricing mechanism
involving charging based on the maximum energy consumption over a specific
period. While such pricing mechanisms are commonly employed, strategies that
directly account for these costs are not as well explored. Therefore we
leverage concepts from _inventory control_ to analyze how peak pricing
influences the optimal scheduling policies.

To focus on this, we simplify the operational process of the industrial
refrigeration, and solely focus on the relationship between the total cooling
and energy expenditures of the refrigeration process. The primary system state
is the facility temperature, denoted as $x_{t}\in\mathbb{R}$ for a given time
$t\in\\{1,\dots,T\\}$. Here, $T$ represents the length of the horizon under
consideration. Using a first-order model of specific heat222According to the
specific heat equation, $x_{t+1}-x_{t}=C_{f}Q_{\mathrm{net}}$, where
$Q_{\mathrm{net}}$ is the net heat transfer and $C_{f}$ is the heat capacity.
For simplicity of notation, we assume that $C_{f}=1$ for the manuscript., we
can succinctly describe the dynamics of the facility’s temperature as follows:

$x_{t+1}=x_{t}-u_{t}+q_{t},$ (1)

where $u_{t}\geq 0$ represents the heat removed from the facility via the
refrigeration system and $q_{t}$ represents the heat influx from the
surrounding environment. We assume that over the horizon, the incoming heat
$q_{t}\geq 0$ is a non-negative random variable that is drawn from a known
distribution $\mathbf{Q}_{t}$ that may be time varying. As an illustration, we
present a scatter plot of possible heat demands over a day333Over June 2023,
power and coefficient of power (COP) estimates were collected from a
refrigeration site of Butterball LLC ® to estimate the incoming refrigeration
loads.  in Figure 2.

Figure 2: We present distributions and respective averages of $\mathbf{Q}_{t}$
over the horizon of a day for a particular refrigeration facility.  Figure 3:
We display the power-heat curves of different compressor types as taken from
[23]. We see that for compressors without variable frequency drives (VFDs),
the power draw and the respective thermal capacity share an affine
relationship.

The control task of the refrigeration system is to generate a sequence of
refrigeration loads $\mathbf{u}\equiv u_{1},\dots,u_{T}$ that minimizes the
overall system cost, dependent on power consumption, while providing the
necessary heat extraction. For a food storage facility, the required
refrigeration comes in the form of a temperature constraint $x_{t}\leq 0$,
where we would like to maintain temperatures to be below freezing. We capture
both the energy costs and temperature constraint violations in the following
stage cost,

$\mathrm{c}_{t}(x_{t},u_{t})=\mathrm{o}_{t}(u_{t})+\mathbb{E}_{q_{t}}\left[\mathrm{h}_{t}(x_{t+1})\right].$
(2)

Here, $\mathrm{o}_{t}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ describes the
cost associated with running the refrigeration at a load of $u_{t}$ and
$\mathrm{h}_{t}:\mathbb{R}\to\mathbb{R}_{\geq 0}$ represents the penalty costs
for set-point deviations from the desired temperature $x_{t}=0$. For each $t$,
we assume that $\mathrm{h}_{t}$ is continuous, convex, and has a minimum at
$0$, i.e. our desired temperature set point. If $x>0$ is positive,
$\mathrm{h}_{t}(x)$ represents the penalty cost for temperature constraint
violations. If $x<0$ is negative, then $\mathrm{h}_{t}(x)$ represents the cost
of excessive cooling which may lead to system inefficiencies. An example of a
reasonable penalty function is

$\mathrm{h}_{t}(x)=\begin{cases}b_{t}\cdot x^{2}&\text{ if }x\geq 0\\\
d_{t}\cdot x^{2}&\text{ if }x<0.\end{cases}$ (3)

with $b_{t}\gg d_{t}$ for each $t$. Since the energy costs are proportional to
the total power consumption of the refrigeration system, we assume an affine
structure for the energy costs with respect to a given refrigeration load,
written as

$\mathrm{o}_{t}(u_{t})=\begin{cases}K+a_{t}\cdot u_{t}&\text{ if }u_{t}>0\\\
0&\text{ if }u_{t}=0,\end{cases}$ (4)

where $K$ represents the setup cost of having the refrigeration system in
operation and $a_{t}$ determines the per-unit cost for refrigeration capacity
at time $t$. The per-unit cost $a_{t}$ may depend heavily on $t$, representing
a potential time-of-use cost structure. This cost model is reflected in the
power-heat curves of the compressors, which represent the majority ($\sim
90\%$) of the energy consumption of the refrigeration process, as shown in
Figure 3.

The cost for peak pricing is reflected in the maximum refrigeration load over
the horizon. More formally, the peak price can be written as
$P\cdot\max\\{u_{t}\text{ for }t\leq T\\}$, where $P\geq 0$ is the scaling
factor associated with the peak costs. To represent peak cost as a terminal
cost, we can introduce an auxiliary state variable $y\in\mathbb{R}_{\geq 0}$
with dynamics

$y_{t+1}=\max\\{y_{t},u_{t}\\}.$ (5)

We note that while the structure of problem without peak pricing is classical
in the inventory control literature (for this, see [24, 25, 26] for similar
results with regards to convex, piecewise affine $\mathrm{o}_{t}$), the
addition of peak pricing is a novel consideration. While not as studied, peak
pricing is extremely significant to the cost structures for industrial
processes.

Consolidating the stage and terminal costs, the total cost of refrigeration
process under the sequence of loads $\mathbf{u}$ with a given peak $y$ is

$J(x,y,\mathbf{u})=\mathbb{E}\left[\sum_{t=1}^{T}\mathrm{c}_{t}(x_{t},u_{t})\right]+P\max\Big{\\{}y,\\{u\\}_{t\leq
T}\Big{\\}},$ (6)

where $x_{t}$ follows from the respective state transition probabilities from
the initial state $x_{1}=x$ and the expectation is taken over the possible
incoming heat $q_{t}\sim\mathbf{Q}_{t}$. The optimal total cost can be written
recursively via the Bellman equation,

$\displaystyle V_{t}(x,y)=\min_{u\geq
0}\Big{\\{}c_{t}(x,u)+\mathbb{E}_{q_{t}}\Big{[}V_{t+1}(x^{+}+q_{t},y^{+})\Big{]}\Big{\\}},$
(7) $\displaystyle V_{T+1}(x,y)=P\cdot y.$ (8)

$V_{T+1}$ represents the terminal cost of the dynamic program. We use
$x^{+}=x-u$ and $y^{+}=\max\\{y,u\\}$ to denote the successor states for
simplicity of notation. The optimal loads can be written in feedback form with
a policy function $\pi^{*}_{t}:\mathbb{R}\times\mathbb{R}_{\geq
0}\to\mathbb{R}_{\geq 0}$ as the argument to the previous optimization
formulation:

$\pi^{*}_{t}(x,y)\in\arg\min_{u\geq
0}\left\\{c_{t}(x,u)+\mathbb{E}_{q_{t}}\left[V_{t+1}(x^{+},y^{+})\right]\right\\}.$
(9)

The main concern of this work is on characterizing the structure of these
optimal policies with respect to peak pricing. We do this through analytical
characterizations in Section III and through simulations in Section IV.

##  III Optimal Policy Characterizations

We first characterize the structure of optimal refrigeration policies with no
peak pricing costs, i.e. when $P=0$ in Eq. (6), where there are only TOU costs
present. In this case, our models align with the standard ones present in
inventory control, and we can directly invoke classical results to get the
structure of the optimal policies. Interestingly, the optimal policies
simplify to a threshold strategy.

######  Proposition 1 ([27]).

Consider the industrial refrigeration problem with a total cost in Eq. (6)
with no peak cost ($P=0$). The optimal policy $\pi^{*}_{t}$ is a threshold
policy of the form

$\pi^{*}_{t}(x,y)=\begin{cases}x-S_{t}&\text{ if }x>s_{t},\\\ 0&\text{ if
}x\leq s_{t},\end{cases}$ (10)

for some $s_{t}\geq S_{t}\in\mathbb{R}$ for every $t\leq T$.

###### Proof.

We can directly use classical results in inventory control with sign changes
(for e.g., see [27], [22, Chapter 2.6], or [28, Chapter 4]) to show the claim.
∎

From the above proposition, we see that the optimal structure comes in the
form of a simple threshold policy, where the refrigeration system is turned
off if the facility is below a certain buffer temperature $s_{t}$ but if
turned on, sets the facility to a lower buffer temperature $S_{t}$. We remark
that if there is no setup costs, or that $K=0$ in Eq. (6), the two buffer
temperatures align, with $s_{t}=S_{t}$.

Threshold policies are commonly implemented in industrial refrigeration to
schedule refrigeration loads. While we have shown that these threshold
policies are optimal without peak pricing, computing closed form solutions for
$s_{t}$ and $S_{t}$ is not possible in general. However, these buffer
temperatures can be derived through data-driven strategies to produce well-
performing threshold parameters. We do this in Section IV.

When peak pricing is introduced to the total cost, however, we will see that
the optimal policies stop corresponding to threshold policies. In fact, for
simple models of refrigeration, we see that the optimal policies can become
quite complicated in the next example.

######  Example 1 (Peak Pricing).

Consider a toy refrigeration scenario, where there is external heat input
$q_{t}\equiv 0$ for all $t$. Let the temperature penalty function be
$\mathrm{h}_{t}(x)=b|x|$ for all $t\leq T$. We set the parameters $P\geq b\geq
a=1$ for the peak, penalty and refrigeration costs, and fix the setup cost to
$K=0$ to frame the problem naturally. Even within this basic model, the
optimal policy can be surprisingly intricate.

To see this, we solve for the value function explicitly in Eq. (7) via the
backwards recursion. Under the model assumptions, the value function $V_{T}$
at time $T$ can be simply expressed as

$V_{T}(x,y)=\min_{u\geq 0}\Big{\\{}u+b|x-u|+P\max\\{y,u\\}\Big{\\}}.$

From the value function, there are three regimes to describe the optimal
policy. These three regimes are depicted in Figure 4, where the optimal policy
is explicitly characterized below.

$\pi^{*}_{T}(x,y)=\begin{cases}0&\text{ if }x\leq 0,\\\ x&\text{ if }0\leq
x\leq y,\\\ y&\text{ if }x\geq y.\end{cases}$

Figure 4: We depict the possible optimal inputs for each of the three regimes,
dependent on $x$. Note that since $P\geq b$, the optimal $u^{*}$ must live
between $[0,y]$, as delineated by the orange markers.

The optimal policy across one time step is a threshold policy with a cap at
$u=y$ due to the peak cost. However, adding another step in the backwards
recursion complicates the optimal policy $\pi^{*}_{T-1}$. For the next step in
the value recursion, we have that

$V_{T-1}(x,y)=\min_{u\geq 0}\Big{\\{}u+b|x-u|+V_{T}(x^{+},y^{+})\Big{\\}},$

where the value function $V_{T}$ can be described in closed form via the
characterization of the optimal policy $\pi^{*}_{T}$ as

$V_{T}(x,y)=\begin{cases}Py-bx&\text{ if }x\leq 0,\\\ Py+x&\text{ if }0\leq
x\leq y,\\\ (P+1-b)y+bx&\text{ if }x\geq y.\end{cases}$

With the expression of the value function $V_{T}$, we can describe the optimal
policy for a horizon of $2$. Since $V_{T}$ is piecewise-affine, the optimal
inputs must occur at the boundary conditions, where $u=0$ or $u=x/2$ or $u=y$
or $u=x$. According to this, we can solve for the optimal policy algebraically
to be

$\pi^{*}_{T-1}(x,y)=\begin{cases}0&\text{ if }x\leq 0,\\\ x&\text{ if }0\leq
x\leq y,\\\ y&\text{ if }y\leq x\leq 2y\text{ or}\\\ &\text{ if }2y\leq
x\text{ and }\\\ &(P+2-3b)y\leq(\frac{P}{2}-\frac{3b}{2}+1)x,\\\
\frac{x}{2}&\text{ otherwise}.\\\ \end{cases}$

As we can see, even with a horizon of two steps, the optimal policy differs
greatly from the original threshold policy in Eq. (10).

While generating closed form expressions for optimal policies are hard to do
in general (as seen in Example 1), if we limit to a horizon of $1$, we have
that a modified threshold policy depicted in Figure 5 is optimal. We
characterize the structure in the following Proposition.

Figure 5: We depict the modified policy.

######  Proposition 2.

Consider the industrial refrigeration problem with a total cost in Eq. (6)
with a horizon $T=1$ with no set up cost ($K=0$ in Eq. (4)). The optimal
policy $\pi^{*}_{1}$ is characterized by

$\pi^{*}_{1}(x,y)=\begin{cases}x-\hat{S}&\text{ if }x\geq\hat{S}+y,\\\
y&\text{ if }x\in[S+y,\hat{S}+y],\\\ x-S&\text{ if }x\in[S,S+y],\\\ 0&\text{
if }x\leq S.\end{cases}$ (11)

for some $s\geq S\in\mathbb{R}_{\geq 0}$ and $\hat{S}\geq S\in\mathbb{R}_{\geq
0}$.

###### Proof.

First, we express the value function $V_{1}$ as below

$V_{1}(x,y)=\min_{u>0}\Big{\\{}f(x^{+})+Py^{+}\Big{\\}}+ax,$

with $f(z)\triangleq\mathbb{E}[h(z+q_{t})]-az$ being a convex scalar function.
Let $S$ be the minimum of $f$. Note that if $x<S$, then by convexity of
$f(z)$, the optimal input is $u=0$. Additionally, if $S\leq x\leq S+y$, the
optimal input is $u=x-S$. Now let $\hat{S}$ be such that the subderivative
$\partial f/\partial z|_{z=\hat{S}}\ni P$. Since $f$ is convex, the minimum is
defined by the first order condition on the subderivative $-\partial
f/\partial z|_{x-u}+P\ni 0$. This gives the first two conditions of the
optimal policy. ∎

While we can compute the optimal one-step policy in closed form, this is not
true in general, as shown in Example 1. However, as shown in the next
Proposition, we can still derive a qualitative threshold-like characterization
of the optimal policy in general. We depict the optimal policy pictorially in
Figure 6.

######  Proposition 3.

Consider the industrial refrigeration problem with a total cost in Eq. (6)
with no set up cost ($K=0$ in Eq. (4)). The optimal policy $\pi^{*}_{t}$ is
characterized by

$\pi^{*}_{t}(x,y)\in\begin{cases}[y,z^{*}_{t}]&\text{ if }g_{t}(y)+y\leq x,\\\
\\{x-g_{t}(y)\\}&\text{ if }x-y\leq g_{t}(y)\leq x,\\\ \\{0\\}&\text{ if
}x\leq g_{t}(y),\end{cases}$ (12)

where $g_{t}(y):\mathbb{R}_{\geq 0}\to\mathbb{R}$ and $z^{*}$ is such that
$x=g_{t}(z^{*}_{t})+z^{*}_{t}$.

###### Proof.

By assumption, we have that $K=0$ in Eq. (6). First we show that for every
$t$, $V_{t}$ is non-decreasing in $y$ via backwards induction. For the base
case, $V_{T+1}(x,y)=Py$ is non-decreasing in $y$, as $P>0$. From Eq. (7), we
have that $V_{t+1}(x^{+},y^{+})$ is a composition of non-decreasing functions
in $y$ by the induction assumption. Since monotonicity is preserved under
expectations and infimum projections, the iterate $V_{t}$ is thus non-
decreasing in $y$ as well, and we have the claim.

Now we show that $V_{t}$ is convex for every $t$. Similarly as before, we show
this via backwards induction. For the base case, $V_{T+1}(x,y)=Py$ is affine
in $y$ and thus convex. For the inductive case, we have that $\mathrm{c}_{t}$
is convex in $x$ and $u$. Additionally, $V_{t+1}$ is convex by the induction
assumption. Moreover, since $x^{+}$ is affine in $x$ and $u$; $y^{+}$ is
convex in $y$ and $u$; and $V_{t+1}$ is non-decreasing in $y$,
$V_{t+1}(x^{+},y^{+})$ is a convex function of $x$, $y$, and $u$. Since
convexity is preserved under expectation, $V_{t}(x,y)$ can be concisely
expressed as

$V_{t}(x,y)=\min_{u\geq 0}f_{t}(x^{+},y^{+})+ax$ (13)

where
$f_{t}(x,y)=\mathbb{E}_{q_{t}}[\mathrm{h}_{t}(x+q_{t})+V_{t+1}(x+q_{t},y)]-ax$
is a convex function. Since $\\{u\in\mathbb{R}:u\geq 0\\}$ is convex,
convexity is preserved under the minimization, and thus $V_{t}(x,y)$ is a
convex function. Thus the claim is shown.

Now we describe the optimal policy $\pi_{t}(x,y)$. Note that we can describe
the optimal policy as

$\pi_{t}(x,y)=\arg\min_{u\geq 0}\big{\\{}f_{t}(x^{+},y^{+})\big{\\}}+ax.$

where $f_{t}(x,y)$ is defined as in Eq. (13). For each $y$, define the
function $g_{t}(y)$ to be

$g_{t}(y)=\max\Big{\\{}\arg\min_{x}f_{t}(x,y)\Big{\\}}.$ (14)

We describe the three cases for the optimal policy as shown in Figure 6 and
Eq. (12). If $x\leq g_{t}(y)$, then observe that $f_{t}(x^{+},y^{+})\geq
f_{t}(x^{+},y)$, since $V_{t+1}$ is non-decreasing in $y$. Additionally, as
$f_{t}$ is convex, and $x\leq g_{t}(y)$, we have that $f_{t}(x^{+},y)\geq
f_{t}(x,y)$ for all $u\geq 0$, and thus $\pi^{*}_{t}(x,y)=0$ when $x\leq
g_{t}(y)$. Likewise, if $x-y\leq g_{t}(y)\leq x$ we have that
$f_{t}(x^{+},y^{+})\geq f_{t}(x^{+},y)\geq f_{t}(g_{t}(y),y)$ for all $u\geq
0$, and thus $\pi_{t}(x,y)=x-g_{t}(y)$.

For the third condition $g_{t}(y)+y\leq x$, we first note that for $u=y$
generates a cost of $f_{t}(x-y,y)\leq f_{t}(x-u,y)$ for any $u\leq y$. Thus
the optimal policy must satisfy $\pi^{*}(x,y)\geq y$. Moreover, let $u=z^{*}$
be the input in which $(x^{+},y^{+})$ intersects the curve $x=g_{t}(y)$. Note
that $f(x^{+},y^{+})\geq f(x^{+},z^{*})\geq f(g_{t}(z^{*}),z^{*})$ for any
$u\geq z^{*}$, since $f$ is non-decreasing in $y$ and $f$ is convex. Thus the
optimal policy must satisfy $\pi^{*}(x,y)\leq z^{*}$. ∎

Figure 6: We depict the possible optimal inputs for each of the three regimes
corresponding to Eq. (12), given the curve $g_{t}(y)$.

######  Remark.

From Proposition 3, we note that if $g_{t}(y)+y\leq x$, the optimal control
action $u$ may lie anywhere in the interval $[y,z^{*}_{t}]$. Where on this
interval depends non-trivially on the relative peak cost $P$ and the specific
holding cost $\mathrm{h}_{t}$ and the current time $t$. However, we maintain a
threshold-like policy structure, where below a buffer temperature $g_{t}(y)$
(dependent on the current peak), the optimal decision is to turn off the
refrigeration system.

######  Example 2.

The optimal policy $\pi^{*}(x,y)$ in Eq. (12) may incur non-intuitive
dependencies on the current peak value $y$; we outline a simplified scenario
that outlines this phenomena. Consider a horizon of $T=2$, where the starting
temperature is $x_{1}=-2$ and the incoming heat is $q_{t}=2$ for all $t$. Let
there be no per-unit cost ($a=0$) and setup cost ($K=0$). Furthermore, we set
the peak costs to $P\gg 1$ and the temperature penalty function to be

$\mathrm{h}_{t}(x)=\begin{cases}bx&\text{ if }x\geq 0\\\ -x&\text{ if
}x<0,\end{cases}$

where $b\gg 1$. We compare two scenarios: where the current peak is $y=1$ and
the current peak is $y=2$. In the first case ($y=1$), observe that the optimal
control sequence is $u_{1}^{*}=1$ and $u_{2}^{*}=1$ to not incur any
temperature violation costs. However, when $y=2$, the optimal control sequence
is $u_{1}^{*}=0$ and $u_{2}^{*}=2$. As such, the optimal policy $\pi_{1}(x,y)$
actually increases with the current peak $y$, counter to intuition that
increasing that the peak should introduce more conservatism to the optimal
policy.

##  IV Simulation Study on Policy Designs

In this section, we evaluate the performance of different policy designs
through a case scenario based on the thermal data (also represented in Figure
2) acquired from a refrigeration facility owned by Butterball LLC®. We depict
the possible incoming heat distributions during the on and off peak hours in
Figure 7.

Figure 7: Incoming heat distributions of on and off peak time periods.

The specifics used in the case study are outlined as follows. We evaluate the
total cost as according to Eq. (6), where each time step $t$ approximately
corresponds to a $12$-hour time window and the trajectory cost is evaluated
over horizon of a month. If the time $t$ corresponds to the on-peak, then the
incoming heat distribution $\mathbf{Q}_{t}$ is exactly the one depicted in
Figure 7 (and respectively for the off-peak time periods). We additionally
normalize the incoming heat values to Megawatts to prevent numerical issues.
We assume that the starting temperature and starting peak value are both $0$.
We use the temperature violation cost of the form in Eq. (3) with $b_{t}=20$
and $d_{t}=1$. We additionally set the respective coefficients for peak costs
to $P=30$, setup costs to $K=0$, and per-unit costs to $a=1$ to reasonably
model a possible refrigeration scenario.

Methodology | Total Cost
---|---
Static Threshold Policy | 36.1
Dynamic Threshold Policy | 26.7
Modified Threshold Policy | 30.7
Dynamic Modified Threshold Policy | 25.5
TABLE I: Cost of Algorithms

Under this example setup, we evaluate the performance of different types of
threshold policies according to the total cost in Eq. (6). We first examine
the performance of the threshold policy in Eq. (10), where $s_{t}=S_{t}$ is
not time varying. Additionally, we examine the dynamic version, where $s_{t}$
is allowed to depend on if $t$ corresponds to an on-peak or off-peak period.
We also examine the performance of the policies stated in Proposition 1 for
time-dependent and time-independent values as well. For these policies, we
optimize for the coefficients through a Monte Carlo search. The performance of
these algorithms are depicted in Table I.

##  V Conclusion

Improving the control algorithms for industrial refrigeration systems can have
significant impacts on energy usage, given their widespread presence and
substantial scale. This work has thus focused on control algorithm design,
specifically investigating the implications of various energy cost-rate
structures, such as time-of-use and peak pricing. To do this, we leverage
tools from inventory control to characterize the structure of the optimal
policies. While simple threshold policies are optimal under solely time-of-use
cost structures, adding significant peak costs can markedly change the policy
structure of the optimal control algorithms. We provide characterizations of
the one-step optimal policies in Proposition 1 and generalize to arbitrary
horizons in Proposition 3, highlighting significant differences from the
classical threshold policies. Through theoretical analysis and simulations, we
underscore the necessity for mindfulness in constructing control policies in
the face of significant peak costs - thus providing valuable insights for more
effective refrigeration management.

## References

  * [1] I. Dincer, Refrigeration systems and applications.  John Wiley & Sons, 2017.
  * [2] W. F. Stoecker, Industrial refrigeration handbook.  McGraw-Hill Education, 1998.
  * [3] F. Fabrega, J. Rossi, and J. d’Angelo, “Exergetic analysis of the refrigeration system in ethylene and propylene production process,” Energy, vol. 35, no. 3, pp. 1224–1231, 2010.
  * [4] S. Tassou, J. S. Lewis, Y. Ge, A. Hadawey, and I. Chaer, “A review of emerging technologies for food refrigeration applications,” Applied Thermal Engineering, vol. 30, no. 4, pp. 263–276, 2010.
  * [5] “Manufacturing energy consumption survey,” tech. rep., US Energy Information Administration, 2018.
  * [6] L. Zhao, W. Cai, X. Ding, and W. Chang, “Model-based optimization for vapor compression refrigeration cycle,” Energy, vol. 55, pp. 392–402, 2013.
  * [7] L. S. Larsen, C. Thybo, J. Stoustrup, and H. Rasmussen, “Control methods utilizing energy optimizing schemes in refrigeration systems,” in 2003 European Control Conference (ECC), pp. 1973–1977, IEEE, 2003.
  * [8] L. F. S. Larsen, Model based control of refrigeration systems.  Department of Control Engineering, Aalborg University, 2006.
  * [9] K. A. Manske, “Performance optimization of industrial refrigeration systems,” 2000.
  * [10] T. G. Hovgaard, S. Boyd, L. F. Larsen, and J. B. Jorgensen, “Nonconvex model predictive control for commercial refrigeration,” International Journal of Control, vol. 86, no. 8, pp. 1349–1366, 2013.
  * [11] X.-H. Yin and S.-Y. Li, “Model predictive control for vapor compression cycle of refrigeration process,” International Journal of Automation and Computing, vol. 15, no. 6, pp. 707–715, 2018.
  * [12] S. E. Shafiei, J. Stoustrup, and H. Rasmussen, “Model predictive control for flexible power consumption of large-scale refrigeration systems,” in 2014 American Control Conference, pp. 412–417, IEEE, 2014.
  * [13] Y. Sun, S. Wang, F. Xiao, and D. Gao, “Peak load shifting control using different cold thermal energy storage facilities in commercial buildings: A review,” Energy conversion and management, vol. 71, pp. 101–114, 2013.
  * [14] Y. Yao and D. K. Shekhar, “State of the art review on model predictive control (mpc) in heating ventilation and air-conditioning (hvac) field,” Building and Environment, vol. 200, p. 107952, 2021.
  * [15] R. Pattison, C. Touretzky, T. Johansson, I. Harjunkoski, and M. Baldea, “Optimal process operations in fast-changing electricity markets: framework for scheduling with low-order dynamic models and an air separation application,” Industrial & Engineering Chemistry Research, vol. 55, pp. 4562–4584, 2016.
  * [16] R. Pattison, C. Touretzky, I. Harjunkoski, and M. Baldea, “Moving horizon closed‐loop production scheduling using dynamic process models,” AIChE Journal, vol. 63, pp. 639–651, 2017.
  * [17] A. Vishwanath, V. Chandan, and K. Saurav, “An iot-based data driven precooling solution for electricity cost savings in commercial buildings,” IEEE Internet of Things Journal, vol. 6, no. 5, pp. 7337–7347, 2019.
  * [18] M. J. Risbeck and J. B. Rawlings, “Economic model predictive control for time-varying cost and peak demand charge optimization,” IEEE Transactions on Automatic Control, vol. 65, no. 7, pp. 2957–2968, 2019.
  * [19] Y. Mo, Q. Lin, M. Chen, and S.-Z. J. Qin, “Optimal online algorithms for peak-demand reduction maximization with energy storage,” in Proceedings of the twelfth ACM international conference on future energy systems, pp. 73–83, 2021.
  * [20] F. Oldewurtel, A. Ulbig, A. Parisio, G. Andersson, and M. Morari, “Reducing peak electricity demand in building climate control using real-time pricing and model predictive control,” in 49th IEEE conference on decision and control (CDC), pp. 1927–1932, IEEE, 2010.
  * [21] S. Axsäter, Inventory control, vol. 225.  Springer, 2015.
  * [22] B. Liu and A. O. Esogbue, Decision criteria and optimal inventory processes, vol. 20.  Springer Science & Business Media, 2012.
  * [23] “Industrial refrigeration facts and figures.” https://bptraining.ornl.gov/wp-content/uploads/2021/05/Ind.-Refrigeration-Facts-Figures.pdf.  Accessed: 2024-03.
  * [24] E. L. Porteus, “Stochastic inventory theory,” Handbooks in operations research and management science, vol. 2, pp. 605–652, 1990.
  * [25] Y. Lu and M. Song, “Inventory control with a fixed cost and a piecewise linear convex cost,” Production and Operations Management, vol. 23, no. 11, pp. 1966–1984, 2014.
  * [26] M. J. Sobel, “Making short-run changes in production when the employment level is fixed,” Operations Research, vol. 18, no. 1, pp. 35–51, 1970.
  * [27] H. Scarf, K. Arrow, S. Karlin, and P. Suppes, “The optimality of (s, s) policies in the dynamic inventory problem,” Optimal pricing, inflation, and the cost of price adjustment, pp. 49–56, 1960.
  * [28] J.-S. J. Song, Research handbook on inventory management.  Edward Elgar Publishing, 2023.


#  LabelGit: A Dataset for Software Repositories Classification using
Attributed Dependency Graphs

Cezar Sas  Bernoulli Institute
University of Groningen
Groningen, Netherlands
<EMAIL_ADDRESS>Andrea Capiluppi  Bernoulli Institute
University of Groningen
Groningen, Netherlands
<EMAIL_ADDRESS>
###### Abstract

Software repository hosting services contain large amounts of open-source
software, with GitHub hosting more than 100 million repositories, from new to
established ones. Given this vast amount of projects, there is a pressing need
for a search based on the software’s content and features. However, even
though GitHub offers various solutions to aid software discovery, most
repositories do not have any labels, reducing the utility of search and topic-
based analysis. Moreover, classifying software modules is also getting more
importance given the increase in Component-Based Software Development.
However, previous work focused on software classification using keyword-based
approaches or proxies for the project (e.g., README), which is not always
available. In this work, we create a new annotated dataset of GitHub Java
projects called LabelGit. Our dataset uses direct information from the source
code, like the dependency graph and source code neural representations from
the identifiers. Using this dataset, we hope to aid the development of
solutions that do not rely on proxies but use the entire source code to
perform classification.

###### Index Terms:

github, software classification, dependency graph, graph neural networks

##  I Introduction

Given the overwhelming amount of repositories hosted on services like GitHub,
there is a need for a functionality that allows for a search and retrieval
based on semantics. As a solution, GitHub proposed a service called
Topics111https://github.com/topics, which allows user to annotate projects
manually and others to search using these topics; and
Collections222https://github.com/collections a curated list of topics where
repositories are displayed. However, these solutions are not perfect and have
various issues. For example, Topics are optional, and GitHub does not suggest
or restrict the user in any way. As a result, there are plenty of similar (or
identical) topics written with various surface forms, making the search less
effective. On the other hand, the Collections list is manually curated;
therefore, it is not scalable to all topics, reducing the effectiveness of
finding repositories, especially those annotated with non-popular topics.

Recently, there has been an increase in the interest of classification of
GitHub repositories.

[1] uses an approach based on bytecode, and the external dependency of the
project, with information from stack overflow to generate a tag cloud. Their
dataset is unavailable.

Sharma et al. [2] uses an combined solution of topic modeling and genetic
algorithms called LDA-GA[3]. They apply topic modeling on the README files,
and optimize the hyperparameters using genetic algorithms. While LDA is an
unsupervised solution, humans are needed to label the topics from the
identified keywords. They release a list of 10,000 examples annotated by their
model (half belonging to ’Others’ category), which was evaluated using 400
manually labeled projects.

ClassifyHub [4] uses an ensemble of 8 naïve classfiers, each using different
features (e.g. file extensions, README, GitHub metadata and more) to perform
the classification task. They use the InformatiCup
2017333https://github.com/informatiCup/informatiCup2017 dataset, which
contains 7 categories.

HiGitClass [5] uses approach for modeling the co-occurence of multimodal
signal in a repository (e.g. user, name of repository, tags, README and more).
They perform the annotation according to a hierarchical taxonomy that is given
as an input with keyword for each leaf node. The release a datset with
taxonomies for two domains: a machine learning taxonomy with 1600 examples,
and a bioinformatics one with 876 projects.

In [6], they use the content of the README files and source code, represented
using TFIDF, as input to a probabilistic model called Multinomial Naïve
Bayesian Network to recommend possible topics. They used 120 popular topics
from GitHub, and released a dataset of around 10,000 distantly annotated
projects in different programming languages.

[7] uses names, descriptions, READMEs, wiki pages, and file names concatenated
together as input to BERT[8], a neural language model, that creates a dense
representation of the input text. Then, to this representation, they apply a
fully connected neural network to predict multiple labels. Their dataset is no
currently available and contains 152K GitHub repositories in various languages
classified using 228 topics.

However, all these solutions perform software classification by classifying
proxies of the projects (e.g., README), which might not always be available.
This problem gets more accentuated if we want to perform classification on
smaller parts of code like components, which require models trained on source
code.

This work’s contributions are a new dataset, called LabelGit[9], containing
annotated GitHub projects for software classification. The classification task
is performed given an attributed graph with continuous features representing
the source code files. We hope this will act as a catalyst for further
research on machine learning for software engineering. Moreover, this can also
be used to develop and evaluate new deep learning methods for real-world
graphs.

Figure 1: Our dataset creation pipeline. We start with an annotated list from
which the annotations are reduced to a smaller set using a manually created
mapping. On the other side, we extract the git project history, and for each
commit in the master branch (green nodes), we run the dependency graph and
feature extraction methods. The final example is made by a graph $G$, a
feature matrix $X$, and a label $y$.

##  II Preliminaries

In the following section, we introduce the concept of creating document
embeddings, a dense vector representation of documents used to capture the
semantic and syntactic information carried by the specific data source.

###  II-A Source Code Representations

The learning of a document’s embedding is usually performed using learned
representations of words first, then combining those to create a document’s
embedding. We used two different types of features to represent documents into
a continuous space (i.e., its embeddings): the resulting embeddings of the
source code files reflect two different interpretations. The first one is
semantic, and based on natural text information. It is created using a neural
language model architecture called fastText. The second solution, called
code2vec, is trained on source code and learns representations by taking
advantage of the structural information of source code.

#### fastText [10]

is a neural language model that creates word embeddings by learning subword
information [10]. fastText words are split into n-grams, and for each it
learns an embedding. The final word embedding is an aggregation of the
embeddings of all n-grams that the word is made of. An advantage of using
subword information like n-grams, is that the smaller information level means
that there is no risk of out-of-vocabulary words when applying to more
technical domains.

#### code2vec[11]

is a neural embedding model specifically designed for source code. The
embeddings are learned using abstract syntax trees paths between variables in
a method, and the task is self supervised by predicting the method name. The
available pretrained model is learned from Java projects, therefore, the
variables representations contain domain specific information.

##  III Datasources

We created our dataset using two sources, the first is GitHub, while the
second is a list of curated GitHub Java projects called Awesome-Java.

#### GitHub

is a repository hosting service that provides an ecosystem for storing,
development, and sharing of software. It is home to many open source
communities and commercial software developers worldwide and contains more
than 100 million repositories.

#### Awesome Lists

are curated repositories containing resources that are useful for a particular
domain. In our case we use Awesome-Java444https://github.com/akullpp/awesome-
java, a GiHub project that aggregates several hundreds of curated Java
frameworks, libraries and software. It contains 69 categories, and an overall
700 repositories. We removed tutorials, URLs to websites, and projects that
failed to be analyzed, and obtained a grand total of 495 projects.

##  IV LabelGit

LabelGit[9] is our new dataset for classification of Java repositories that
uses source code semantics and dependency graphs. It is based on a curated and
annotated list of Java projects. The pipeline used to create the LabelGit
dataset is represented in Figure 1, and can be divided in the following steps:
1. extraction of the dependency graph;  2. extraction of the features;  3.
label remapping;  4. data augmentation.

We made our source code555https://github.com/SasCezar/LabelGit and
dataset666https://zenodo.org/record/4459080#.YBPENOhKiF4 publicly available.
It contains the annotations in CSV format, the zip file for the graphs, and
the two other zip files for the features. The dataset size is around 15 GB
compressed, and 52 GB when uncompressed.

###  IV-A Dependency Graph

The first step of our approach is the extraction of the dependency graph for
each project in our sample. Using the Arcan [12] tool, we obtained the
complete set of nodes and the edges describing the dependencies between
classes. Edge weights were also extracted to describe the number of uses [13]
of one class by the other. For example, the antlr4 project contains a
dependency between ParserATNFactory.java (from the org.antlr.v4.automata
package) and LexerGrammar.java (from the org.antlr.v4.tool package), as the
former imports the latter.

The distribution of the number of vertices and edges for the graphs in our
dataset are shown in Figure 2.

Figure 2: Distribution of the number of vertices (top) and edges (bottom) in
the augmented dataset’s graphs.

The dependency graph is stored using GraphML (.graphml), a XML based graph
storage file format. The nodes in the graph contain an attribute with the
relative path of the file in the project folder, this is used to identify the
source file in the feature files and perform alignment.

The .graphml file name allows to easily identify the project and the git
commit, and it is structured as follows:

[project name]-[commit number]-[commit sha].graphml

where [commit number] is the index of the commit in the version history tree,
and [commit sha] is the id of the commit.

###  IV-B Semantic Features

The extraction of features representing the documents is performed with two
different language models: fastText and code2vec. The first one was trained
with natural text, the second with source code in various programming
languages.

We apply the two embedding methods on the identifiers extracted from the
source code file. For the extraction of the identifiers contained in the
source code, we used the tree-sitter777https://github.com/tree-sitter/tree-
sitter parser generator tool. It makes easy to get the identifiers, without
keywords, from the annotated concrete syntax tree created using a grammar for
Java code. We clean the identifiers by separating the camel case strings into
words, lower case every token, and remove common Java terms that do not add
much semantically (e.g., ‘main’, ‘println’, etc).

The features are stored in a textual file (.vec), with each file containing
the embeddings for a particular project at a specific commit. The project’s
source code files are described with a row containing the relative path of the
file, and the vector of features with each feature separated by a space.

The filename follows the same structure as before, with the only difference in
the extension:

[project name]-[commit number]-[commit sha].vec

###  IV-C Label Mapping

The Awesome-Java list contains 69 annotations (or categories), and on average,
each category contains some 8 projects, making any classification task very
challenging. Moreover, those categories have different granularity levels and
can represent very general concepts or very detailed keywords.

For these reasons, we decided to manually reduce the original categories to a
smaller set of 13 (see Label column in Table I). This mapping was performed
manually, in a hierarchical fashion, by the authors.

The initial and final annotated labels are stored as a CSV file
(annotations.csv). The file has the followings schema:

  * •

project.name: name of the project;

  * •

project.desc: short description of the project from Awesome-Java;

  * •

project.link: URL to the GitHub repository;

  * •

category: Awesome-Java annotation;

  * •

category.desc: short description of the category from Awsome-Java;

  * •

label: mapping of the Awsome-Java category into one of the labels of the
smaller set.

###  IV-D Data Augmentation

Given the small amount of projects contained in the Awesome-Java list, and the
complexity of the classification task, we decided to perform a ‘data
augmentation’[14] which has been shown to aid generalization in deep learning
models in various domains[15, 16], including software engineering[17] step.
Given our data’s nature, we can synthetically generate a large amount of data,
while maintaining its original integrity. This is done by making use of the
history of a project’s repository: instead of taking a single snapshot of it,
we take all the snapshots from its master branch. Therefore, for each
annotated project, we are able to extract multiple different examples that can
be used as added ‘synthetic’ projects in the training phase.

This augmentation approach allowed us to increase the initial 495 samples to
11,502, with an average increase of 23 times the amount of samples per
project. In Table I we can see the number of examples before and after the
augmentation separated per class.

TABLE I: Distribution of the number of examples before (Projects), and after (Samples) the data augmentation process. Label | Projects | Samples | Increase (x)
---|---|---|---
Introspection | 32 | 744 | 23
CLI | 8 | 142 | 17
Data | 49 | 1,088 | 22
Development | 100 | 2,306 | 23
Graphical | 11 | 226 | 22
Miscellaneous | 59 | 1,729 | 20
Networking | 25 | 503 | 20
Parser | 41 | 935 | 22
Scientific/Engineering | 39 | 915 | 23
Security | 14 | 249 | 18
Server | 37 | 727 | 20
Testing | 42 | 974 | 23
Web | 38 | 964 | 25
Total | 495 | 11,502 | 23

To evaluate the overlap resulting from the augmentation approach, we measured
the graph difference between two consecutive commits of a project. In order to
do so, we take the absolute value of the difference of size (sum of the number
of vertices and the number of edges) of the two graphs, and normalize by the
maximum size of the considered graphs888We chose to evaluate this proxy
because computing an exact measure would be extremely demanding, as graph edit
distance on large graph requires substantial amounts of memory and time..

Formally, given two graphs representing consecutive analyzed versions of a
specific project, $\textbf{G}_{i}=(V_{i},E_{i})$,
$\textbf{G}_{j}=(V_{j},E_{j})$, their difference is:

$diff(\textbf{G}_{i},\textbf{G}_{j})=\frac{abs((|V_{i}|+|E_{i}|)-(|V_{j}|+|E_{j}|))}{max(|V_{i}|+|E_{i}|,|V_{j}|+|E_{j}|)}$

where $|\cdot|$ is the cardinality of the set.

The average difference between the graphs of two subsequent commits is around
12%, with a mean interval ($abs(i-j)$) between the analyzed consecutive
versions of around 550 commits.

##  V Uses and Future Work

#### Uses

the dataset that we have proposed will aid the development and evaluation of
new deep learning models, especially those designed to work on non-euclidean
data[18], such as Graph Neural Networks (GNN)[19]. These models will allow to
solve tasks such as classification of software projects, and more in general
of real-world graph data. Moreover, by performing the classification, these
models can also learn better graph representation: this, in turn, can be used
for other downstream tasks, such as the identification of architectural
smells[20], and software similarity[21]. Furthermore, the learned
representation and the longitudinal nature of our dataset allow us to analyse
how software semantics changes with the development and increase in the
project’s size. Lastly, we can transfer the model’s knowledge learned from the
coarse project to finer parts of software, like components, to classify them
and aid their retrieval for reuse purposes during development.

#### Future Work

as future improvements to the dataset, we plan to extract a bottom-up taxonomy
from the data and provide it with multi-label classification capabilities.
This is because a single project might encompass different characteristics and
contain components that belong to different categories. We foresee that this
process will be based on GitHub’s topics and that those will be reduced using
hierarchical clustering, as shown in our process above. Moreover, we plan to
expand the sample by increasing the number of projects and examples, allowing
for a better generalization of models.

#### Threats to validity

the main threat to our approach is the use of the curated list for both the
gathering of the projects, and as the starting point for creating our
taxonomy. This choice can be considered arbitrary; however, we use this as a
starting point and something to improve on; moreover, the development of
various approaches will not be affected by the taxonomy.

Another threat is the data augmentation technique adopted as it might not
allow for good generalization as they are still only 495 projects; however,
this is reduced by the large difference between the graphs and the fact that
various studies showed how data augmentation improves generalization in low
data domains.

##  VI Conclusions

In this work, we proposed LabelGit, a new dataset for software classification
of Java repositories. Compared to previous solutions, ours is designed to be a
direct representation of the source code by creating an attributed dependency
graph, where each node is represented with semantic features. We expect that
LabelGit will increase the interest in developing better solutions for
software classification, and more in general, for real-world attributed graph
classification.

## Acknowledgements

We would like to thank the Center for Information Technology of the University
of Groningen for their support and for providing access to the Peregrine high
performance computing cluster.

## References

  * [1] S. Vargas-Baldrich, M. Linares-Vásquez, and D. Poshyvanyk, “Automated tagging of software projects using bytecode and dependencies (n),” in _2015 30th IEEE/ACM International Conference on Automated Software Engineering (ASE)_ , 2015, pp. 289–294.
  * [2] A. Sharma, F. Thung, P. S. Kochhar, A. Sulistya, and D. Lo, “Cataloging github repositories,” in _Proceedings of the 21st International Conference on Evaluation and Assessment in Software Engineering_ , ser. EASE’17. New York, NY, USA: Association for Computing Machinery, 2017, p. 314–319.
  * [3] A. Panichella, B. Dit, R. Oliveto, M. Di Penta, D. Poshynanyk, and A. De Lucia, “How to effectively use topic models for software engineering tasks? an approach based on genetic algorithms,” in _2013 35th International Conference on Software Engineering (ICSE)_ , 2013, pp. 522–531.
  * [4] M. Soll and M. Vosgerau, “Classifyhub: An algorithm to classify github repositories,” in _KI 2017: Advances in Artificial Intelligence_ , G. Kern-Isberner, J. Fürnkranz, and M. Thimm, Eds. Cham: Springer International Publishing, 2017, pp. 373–379.
  * [5] Y. Zhang, F. F. Xu, S. Li, Y. Meng, X. Wang, Q. Li, and J. Han, “Higitclass: Keyword-driven hierarchical classification of github repositories,” in _2019 IEEE International Conference on Data Mining (ICDM)_ , Nov 2019, pp. 876–885.
  * [6] C. Di Sipio, R. Rubei, D. Di Ruscio, and P. T. Nguyen, “A multinomial naïve bayesian (MNB) network to automatically recommend topics for github repositories,” in _Proceedings of the Evaluation and Assessment in Software Engineering_ , ser. EASE ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 71–80. [Online]. Available: https://doi.org/10.1145/3383219.3383227
  * [7] M. Izadi, S. Ganji, and A. Heydarnoori, “Topic recommendation for software repositories using multi-label classification algorithms,” _ArXiv_ , vol. abs/2010.09116, 2020.
  * [8] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of deep bidirectional transformers for language understanding,” in _Proceedings of the 2019 Conference of the North American Chapter of the ACL: HLT, Volume 1_. Minneapolis, Minnesota: Association for Computational Linguistics, Jun. 2019, pp. 4171–4186.
  * [9] C. Sas and A. Capiluppi, “LabelGit: A dataset for software repositories classification using attributed dependency graphs,” Feb. 2021. [Online]. Available: https://doi.org/10.5281/zenodo.4459080
  * [10] P. Bojanowski, E. Grave, A. Joulin, and T. Mikolov, “Enriching word vectors with subword information,” _Transactions of the Association for Computational Linguistics_ , vol. 5, pp. 135–146, 2017. [Online]. Available: https://www.aclweb.org/anthology/Q17-1010
  * [11] U. Alon, M. Zilberstein, O. Levy, and E. Yahav, “Code2vec: Learning distributed representations of code,” _Proc. ACM Program. Lang._ , vol. 3, no. POPL, Jan. 2019. [Online]. Available: https://doi.org/10.1145/3290353
  * [12] F. A. Fontana, I. Pigazzini, R. Roveda, D. Tamburri, M. Zanoni, and E. Di Nitto, “Arcan: A tool for architectural smells detection,” in _2017 IEEE International Conference on Software Architecture Workshops (ICSAW)_. IEEE, 2017, pp. 282–285.
  * [13] L. Pruijt, C. Köppe, J. M. van der Werf, and S. Brinkkemper, “The accuracy of dependency analysis in static architecture compliance checking,” _Software: practice and Experience_ , vol. 47, no. 2, pp. 273–309, 2017.
  * [14] C. Shorten and T. M. Khoshgoftaar, “A survey on image data augmentation for deep learning,” _Journal of Big Data_ , vol. 6, no. 1, pp. 1–48, 2019.
  * [15] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [16] L. Perez and J. Wang, “The effectiveness of data augmentation in image classification using deep learning,” _arXiv preprint arXiv:1712.04621_ , 2017\.
  * [17] Q. Mi, Y. Xiao, Z. Cai, and X. Jia, “The effectiveness of data augmentation in code readability classification,” _Information and Software Technology_ , vol. 129, p. 106378, 2021. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0950584920301464
  * [18] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, “Geometric deep learning: Going beyond euclidean data,” _IEEE Signal Processing Magazine_ , vol. 34, no. 4, pp. 18–42, 2017.
  * [19] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 1, pp. 61–80, 2009.
  * [20] I. Pigazzini, “Automatic detection of architectural bad smells through semantic representation of code,” in _Proceedings of the 13th European Conference on Software Architecture - Volume 2_ , ser. ECSA ’19. New York, NY, USA: Association for Computing Machinery, 2019, p. 59–62.
  * [21] A. Capiluppi, D. Di Ruscio, J. Di Rocco, P. T. Nguyen, and N. Ajienka, “Detecting java software similarities by using different clustering techniques,” _Information and Software Technology_ , vol. 122, p. 106279, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S095058492030029X


11institutetext: Physics Department, Tsinghua University, Beijing 100084,
China 22institutetext: Helmholtz Research Academy Hesse for FAIR (HFHF), GSI
Helmholtz Center for Heavy Ion Physics, Campus Frankfurt, 60438 Frankfurt,
Germany 33institutetext: Institut für Theoretische Physik, Johann Wolfgang
Goethe-Universität,Max-von-Laue-StraÃe 1, D-60438 Frankfurt am Main, Germany

# Gluon decay into heavy quark pair under a strong magnetic field

Shile Chen 11<EMAIL_ADDRESS>Jiaxing Zhao 2233 Pengfei Zhuang 11

###### Abstract

Due to the extreme large magnetic field produced in the initial stage of non-
central heavy-ion collision, the dynamical process of gluon decay into heavy
quark pair will take place under an external field rather than in vacuum.
Unlike in the vacuum case, where the process is forbidden by energy momentum
conservation, under the external field, a process emerges considering the
background energy which recovers the conservation. We calculate the gluon
decay rate at leading order under a uniform magnetic field.

It is widely accepted that the strongest electromagnetic field in nature is
generated in high energy nuclear collisions Kharzeev:2007jp ; Voronyuk:2011jd
; Deng:2012pc ; Tuchin:2013ie . The maximum of the magnetic field can reach
$5m_{\pi}^{2}\sim 0.1\ \rm GeV^{2}$ in Au+Au collisions at top RHIC energy and
almost $70m_{\pi}^{2}\sim 1\ \rm GeV^{2}$ in Pb+Pb collisions at LHC energy
Deng:2012pc ; Tuchin:2013ie , where $m_{\pi}$ is the pion mass. The strong
magnetic field will cause many novel phenomena like chiral magnetic
effect(CME) Kharzeev:2007jp , photoproduction of dileptons and quarkonia in
peripheral and ultra-peripheral collisions ALICE:2015mzu , splitting of
$D^{0}$ and $\bar{D}^{0}$ directed flows Das:2016cwd , and spin-polarized
difference between $\Lambda$ and $\bar{\Lambda}$ Guo:2019joy .

In our previous study, we calculate the heavy quark production with the
elementary process $gg\to Q\bar{Q}$ under the extreme strong magnetic field
and show the analytical results in the strong field limit LLL(Lowest Landau
Level). Because of the dimension reduction in phase space, the cross section
of this process has a divergence around the incoming energy threshold. Then
the heavy quark pair production in nucleus-nucleus collisions will be
extremely enhanced at low transverse momentum and suppressed at high
transverse momentumChen:2024lmp . The existence of external magnetic field
which breaks the rotational invariance and only the momentum along the
magnetic field is conserved. The production process becomes anisotropic and
strongly dependent on the direction of motion.

To verify the validity of this LLL approximation, we put our calculation a
step forward, extending to the next Landau level. From the result, we could
conclude that if the partons mainly distributed in the small x region and the
heavy quark produced at low $p_{T}$ region with external magnetic field $\sim
10GeV^{2}$ then the process could be well described at LLL. In the
calculation, we found that the cross section of such an elementary process
will encounter a divergence for $u-$ and $t-$ channel if we set the final
quarks up to $n=1$ state. The inducement of this divergence is not from the
threshold of the initial energy but the $on-shell$ of the internal quark,
which means the Feynman diagram could be cut, and the sub-process is just the
decay of the on shell gluon to quark anti-quark pair. In vacuum, the
elementary process like Fig. 1 is forbidden due to the energy-momentum
conservation because the massless gluon has no rest mass. However, with the
background field, this process will no longer be forbidden any more because
the strong magnetic field could compensate the energy. This process has been
widely studied in the astro physics when they consider the QED processes on a
magnetars. In heavy-ion collision, when considering the pre equilibrium stage
with magnetic field, the gluon decay has barely contributed to the quark
statistics. In this proceeding paper, we will calculate and discuss this
process under LLL.

Figure 1: Elementary process of the gluon decay into quark anti-quark pair.

The Feynman rules associated with quarks are controlled by Dirac equation

$\left[i\gamma^{\mu}\left(\partial_{\mu}+iqA_{\mu}\right)-m\right]\psi=0,$ (1)

where $m$ is the quark mass, $q$ the quark electric charge, and $A_{\mu}$ the
electromagnetic potential. We consider an external magnetic field $B$ in the
direction of $z$-axis and choose the Landau gauge with $A_{0}=0$ and
${\bm{A}}=Bx{\bm{e}}_{y}$. In this case, the momentum along the $x$-axis is
not conserved. Taking into account the Landau energy levels for a fermion
moving in an external magnetic field, the stationary solution of the Dirac
spinor can be written as

$\psi_{n,\sigma}^{-}(x,p)=e^{-ip\cdot x}u_{n,\sigma}({\bm{x}},p),\ \ \ \
\psi_{n,\sigma}^{+}(x,p)=e^{ip\cdot x}v_{n,\sigma}({\bm{x}},p)$ (2)

where the quantum number $n$ is Landau level, and the four-momentum is defined
as $p_{\mu}=(\epsilon,0,p_{y},p_{z})$ with $p_{y}=aqB$ controlled by the
center of gyration $a$.

With the stationary solution of fermion field above, we could obtain the
S-matrix of the decay process $g\to q\bar{q}$ as

$S_{fi}=-ig\int
d^{4}x\bar{\psi}_{-}^{(\sigma_{-})}(x)_{p_{z-},0,a_{-}}\gamma_{\mu}t^{c}A^{\mu}_{c}(x)\psi_{+}^{(\sigma_{+})}(x)_{p_{z+},0,a_{+}}$
(3)

when we proceed our calculation under LLL, we could obtain

$\displaystyle S_{fi}$ $\displaystyle=$ $\displaystyle
t^{c}\frac{g}{2\sqrt{2\omega}L^{7/2}}e^{-\frac{i}{4}\omega\sin\theta[4a_{-}\cos\phi+\lambda_{B}^{2}\omega\sin\theta(-i+\sin(2\phi))]}$
(4) $\displaystyle\times$
$\displaystyle\epsilon_{z}\frac{p_{z-}(-p_{z-}+\omega\cos\theta)+(m+E_{-})(m+E_{+})}{\sqrt{m(m+E_{-})}\sqrt{m(m+E_{+})}}(2\pi)^{3}\delta_{fi}^{3}(E,p_{y},p_{z})$

where
$\omega\hat{k}=\omega(1,\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ is
the momentum of the initial gluon and $\epsilon_{z}$ is the z-component of its
polarization vector.

The two independent polarization modes of gluon is $\mathcal{O}$ mode and
$\mathcal{E}$ mode, with
$\epsilon_{\mathcal{O}}=(-\cos\theta\cos\phi,-\cos\theta\sin\phi,\sin\theta)$
and $\epsilon_{\mathcal{E}}=(\sin\phi,-\cos\phi,0)$. Since the S-matrix is
proportional to $\epsilon_{z}$, only $\mathcal{O}$-mode remains. If here we
set $\theta=\pi/2$, then $\epsilon_{z}=1$ and $p_{z+}=-p_{z-},\
E_{+}=E_{-}=E$, we could obtain a simplified formulation

$S_{fi}=t^{c}\frac{gm}{E\sqrt{2\omega}L^{7/2}}e^{-\frac{1}{4}\lambda_{B}^{2}\omega^{2}+i\Phi}(2\pi)^{3}\delta_{fi}^{3}(E,p_{y},p_{z})$
(5)

where $\Phi$ is a phase term which will not contribute to the decay rate. And
this S-matrix has the exact the same formula with Kostenko:2018cgv ;
Kostenko:2019was  except the color factor. The decay rate of the process $g\to
Q\bar{Q}$

$\displaystyle\Gamma$ $\displaystyle=$ $\displaystyle
g^{2}\sum_{c}\frac{1}{T}\int L\frac{dp_{z+}}{2\pi}\int
L\frac{dp_{z-}}{2\pi}\int L\frac{da_{+}}{2\pi\lambda_{B}^{2}}\int
L\frac{da_{-}}{2\pi\lambda_{B}^{2}}|S_{fi}|^{2}$ (6) $\displaystyle=$
$\displaystyle g^{2}\frac{4q^{2}m^{2}}{12\pi
E^{2}\omega\lambda_{B}^{2}}e^{-\lambda_{B}^{2}\omega^{2}/2}\int\frac{da_{-}}{L}\int
dp_{z-}\delta(\omega-2E)$ $\displaystyle=$ $\displaystyle
g^{2}\frac{4qBm^{2}}{3\pi\omega^{2}\sqrt{\omega^{2}-4m^{2}}}e^{-\lambda_{B}^{2}\omega^{2}/2}$

Figure 2: The decay rate of $g\to Q\bar{Q}$ with final quark mass $m=1.5$GeV
changes with initial gluon energy(left) and polarization angle(right).

From the numerical result, we could see that the decay rate decreasing with
incoming energy has the same property with the scattering cross section $gg\to
Q\bar{Q}$, which has a threshold depending on the rest mass of final fermions
and increasing with magnetic field. The difference appears when considering
the polarization angle dependence. The decay rate is zero if $\theta=0,\pi$,
which means the gluon is moving along the magnetic field. In that case, the
elementary process will not be affected by the external magnetic field and is
absolutely forbidden by the energy momentum conservation. Also since the
exponential decay term $e^{-\lambda_{B}^{2}\omega^{2}\sin^{2}\theta/2}$
exists, the decay rate will not reach its maximum at $\theta=\pi/2$.

In summary, we derived the gluon decay rate under the external magnetic field
with strong field limit. This process is barely considered in the initial
stage quark production since the magnetic field is believed to decay at a very
fast speed when the system enters the QGP phase. From the result, we could
conclude that the decay rate of gluon increases with the magnitude of magnetic
field which means a stronger magnetic field can induce a faster decay of the
gluon even if the gluon has no interaction with magnetic field. As our
previous work mentioned, the LLL is more suitable for a lighter quark, and we
believed that the light quark production will also be affected by this
process.

Acknowledgement: The work is supported by the NSFC grants No. 12075129, the
Guangdong Major Project of Basic and Applied Basic Research
No.2020B0301030008, the funding from the European Unionâs Horizon 2020
research, and the innovation program under grant agreement No. 824093
(STRONG-2020).

## References

  * (1) D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, The Effects of topological charge change in heavy ion collisions: ’Event by event P and CP violation’, Nucl. Phys. A 803 (2008) 227–253, [arXiv:0711.0950].
  * (2) V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, and S. A. Voloshin, (Electro-)Magnetic field evolution in relativistic heavy-ion collisions, Phys. Rev. C 83 (2011) 054911, [arXiv:1103.4239].
  * (3) W.-T. Deng and X.-G. Huang, Event-by-event generation of electromagnetic fields in heavy-ion collisions, Phys. Rev. C 85 (2012) 044907, [arXiv:1201.5108].
  * (4) K. Tuchin, Particle production in strong electromagnetic fields in relativistic heavy-ion collisions, Adv. High Energy Phys. 2013 (2013) 490495, [arXiv:1301.0099].
  * (5) ALICE Collaboration, J. Adam et al., Measurement of an excess in the yield of $J/\psi$ at very low $p_{\rm T}$ in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV, Phys. Rev. Lett. 116 (2016) 222301, [arXiv:1509.08802].
  * (6) S. K. Das, S. Plumari, S. Chatterjee, J. Alam, F. Scardina, and V. Greco, Directed Flow of Charm Quarks as a Witness of the Initial Strong Magnetic Field in Ultra-Relativistic Heavy Ion Collisions, Phys. Lett. B 768 (2017) 260–264, [arXiv:1608.02231].
  * (7) Y. Guo, S. Shi, S. Feng, and J. Liao, Magnetic Field Induced Polarization Difference between Hyperons and Anti-hyperons, Phys. Lett. B 798 (2019) 134929, [arXiv:1905.12613].
  * (8) S. Chen, J. Zhao and P. Zhuang, [arXiv:2401.17559 [hep-ph]].
  * (9) Z. Wang, J. Zhao, C. Greiner, Z. Xu, and P. Zhuang, Incomplete electromagnetic response of hot QCD matter, Phys. Rev. C 105 (2022) L041901, [arXiv:2110.14302].
  * (10) L. Yan and X.-G. Huang, Dynamical evolution of magnetic field in the pre-equilibrium quark-gluon plasma, arXiv:2104.00831.
  * (11) Y. Chen, X.-L. Sheng, and G.-L. Ma, Electromagnetic fields from the extended Kharzeev-McLerran-Warringa model in relativistic heavy-ion collisions, Nucl. Phys. A 1011 (2021) 122199, [arXiv:2101.09845].
  * (12) A. Kostenko and C. Thompson, QED Phenomena in an Ultrastrong Magnetic Field. I. Electron–Photon Scattering, Pair Creation, and Annihilation, Astrophys. J. 869 (2018) 44, [arXiv:1904.03324].
  * (13) A. Kostenko and C. Thompson, QED Phenomena in an Ultrastrong Magnetic Field. II. Electron-Positron Scattering, $e^{\pm}$-Ion Scattering, and Relativistic Bremsstrahlung, Astrophys. J. 875 (2019) 23, [arXiv:1904.03325].