# Transverse momentum weighted Sivers asymmetries in SIDIS and Drell-Yan
processes at COMPASS

Shi-Chen Xue  School of Physics and Microelectronics, Zhengzhou University,
Zhengzhou, Henan 450001, China  Shuailiang Yang  School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China  Xiaoyu
Wang<EMAIL_ADDRESS>School of Physics and Microelectronics, Zhengzhou
University, Zhengzhou, Henan 450001, China  De-Min Li<EMAIL_ADDRESS>School
of Physics and Microelectronics, Zhengzhou University, Zhengzhou, Henan
450001, China  Zhun Lu<EMAIL_ADDRESS>School of Physics, Southeast
University, Nanjing, Jiangsu 211189, China

###### Abstract

We investigate the transverse momentum weighted Sivers asymmetries in the
processes with transversely polarized proton target, including the
$\frac{P_{hT}}{zM_{p}}\sin(\phi_{h}-\phi_{S})$ weighted asymmetry in charged
hadron production in semi-inclusive deeply inelastic scattering (SIDIS) and
the $\frac{q_{T}}{M_{p}}\sin(\phi_{S})$ weighted asymmetry in $\pi^{-}p$
Drell-Yan process. Due to the integration over the transverse momentum, the
weighted asymmetries can be expressed as the product of the transverse-moments
of the transverse momentum dependent (TMD) parton distribution functions
(PDFs) and fragmentation functions (FFs). Using the parametrization for the
Sivers function of the proton, the unpolarized PDFs of proton and pion, and
the unpolarized FF of charged hadron, we present the numerical calculation for
weighted Sivers asymmetries in SIDIS and Drell-Yan processes at the kinematics
of COMPASS, and compare them with experimental data. We find that our
prediction on the weighted Sivers asymmetry in SIDIS process is in agreement
with the recent COMPASS measurement. Due to the relatively large uncertainties
of the preliminary data from COMPASS Drell-Yan program, high precision
experimental data are needed to test the sign change property of the Sivers
function between SIDIS and Drell-Yan processes and to constrain the sea quark
Sivers function.

##  I INTRODUCTION

Since the European Muon Collaboration performed the measurement of the
fraction carried by the internal quark of the proton spin EuropeanMuon:1987isl
; EuropeanMuon:1989yki , the significant deviation from the theoretical
prediction based on the quark model has inspired a large number of
experimental and theoretical research for the proton spin structure.
Understanding the spin structure of the nucleon has become one of the main
goals in QCD and hadronic physics. Normally, the 3-dimensional partonic
structure of the nucleon is described by TMD PDFs. Sivers function
$f_{1T}^{\perp}$ Sivers:1989cc  is one of the eight TMD PDFs at leading twist,
which is a time-reversal odd (T-odd) distribution function that denotes the
asymmetric distribution of unpolarized quarks inside a transversely polarized
nucleon. Arising from the correlation between the internal quark transverse
momentum and the nucleon transverse spin, the Sivers function manifests novel
spin structure of hadrons within the twist-2 approximation of QCD parton
model. Because of its T-odd property, the Sivers function and its chiral-odd
partner the Boer-Mulders function are forbidden by the naive time-reversal
invariance of QCD, thus the very existence of the T-odd distribution functions
was not so obvious. However, the situation was changed after the calculations
in Refs. Brodsky:2002cx ; Brodsky:2002rv ; Boer:2002ju , which show that the
T-odd distributions can actually survive using spectator model calculations
incorporating gluon exchange between the struck quark and the spectator. The
time-reversal-invariance argument was reexamined in Ref. Collins:2002kn ,
which shows that the gauge-link in the operator definition of the correlator
guarantees the T-odd distribution functions to be nonzero. Particularly, the
presence of the gauge-link indicates that the Sivers function and the Boer-
Mulders function have opposite sign between semi-inclusive deeply inelastic
scattering (SIDIS) and Drell-Yan processes Brodsky:2002rv ; Brodsky:2002cx ;
Collins:2002kn , a significant prediction by QCD. The verification of this
sign change is one of the most fundamental tests of our understanding of the
QCD dynamics and the factorization scheme, and it is also the main pursue of
the existing and future Drell-Yan facilities.

The transverse single-spin asymmetries (TSSAs) related to the Sivers function
in high energy scattering processes turn out to be important experimental
tools to investigate the information of Sivers function. During the past two
decades, there are plenty of experimental measurements in SIDIS and Drell-Yan
processes. Measurements in SIDIS process have been made by the HERMES
Collaboration HERMES:2004mhh ; HERMES:2009lmz ; HERMES:2020ifk , COMPASS
Collaboration COMPASS:2008isr ; COMPASS:2010hbb ; COMPASS:2012dmt ;
COMPASS:2016led , and JLab Collaboration JeffersonLabHallA:2011ayy ;
JeffersonLabHallA:2014yxb . Meanwhile, the first measurement of TSSAs from
Sivers function in the pion-induced Drell-Yan process is reported by COMPASS
COMPASS:2017jbv . Measurement of TSSAs of weak boson produced in transversely
polarized proton-proton collisions has also been performed by STAR experiment
at RHIC STAR:2015vmv . The TSSAs data were adopted by several groups
Anselmino:2005ea ; Efremov:2004tp ; Collins:2005ie ; Vogelsang:2005cs ;
Anselmino:2008sga ; Anselmino:2012aa ; Bacchetta:2011gx ; Echevarria:2014xaa ;
Anselmino:2016uie ; Martin:2017yms ; Boglione:2018dqd ; Bury:2021sue  to
extract the Sivers function from parametrization. From the theoretical side,
the quark Sivers function has been intensively studied by various models, such
as the spectator model Brodsky:2002cx ; Boer:2002ju ; Bacchetta:2003rz , the
light-cone quark model Lu:2004hu ; Pasquini:2010af , the non-relativistic
constituent quark model Courtoy:2008vi , and the MIT bag model Yuan:2003wk ;
Courtoy:2008dn . Very recently, the sea quark Sivers function has been
estimated from the light-cone wave function in Ref. Luan:2022fjc . Although
there are extensive studies about Sivers function by theoretical studies,
experimental measurements, and phenomenological analyses, what we know about
Sivers function is still limited. One of the reasons is that the transverse
momentum dependence of the Sivers function is tricky theoretically. The
transverse spectrum in the large transverse momentum (the transverse momentum
of the final-state hadron in SIDIS process and the transverse momentum of the
final-state dilepton in Drell-Yan process) region of the processes is expected
to be described by the fixed order perturbative calculation of QCD, while in
the small transverse momentum region, the transverse spectrum is described by
the resummation of the soft-gluon. Thus, it is relatively difficult to
describe the physical observables because of the complicated transverse
momentum dependent effects.

Choosing proper transverse momentum dependent weighting function $W$ with the
transverse momentum integrated out, the authors in Refs. Kotzinian:1995cz ;
Boer:1997nt  proposed the idea of transverse momentum weighted asymmetry,
which has been widely applied to SIDIS process Pasquini:2011tk ;
Anselmino:2011ay ; xue:2021svd ; Bacchetta:2010si  and Drell-Yan process
Sissakian:2005yp ; Lu:2011qp ; Wang:2017onm ; Bacchetta:2010si ; Liu:2021boj .
Weighted asymmetry can be expressed as the product of collinear functions
instead of the complicated convolution of the TMD functions considering the
TMD effects. COMPASS experiment has the unique advantage to explore the sign
change of the Sivers function since it has almost the same setup for SIDIS and
Drell-Yan process, which may reduce the uncertainties in the extraction of the
Sivers function from the two kinds of measurements. Recently, COMPASS
Collaboration has reported the measurements of weighted Sivers asymmetry in
charged hadron produced SIDIS process in Ref. COMPASS:2018ofp . The
preliminary measurement at COMPASS Drell-Yan program was also presented in
Ref. Longo:2019bih , in which $\pi^{-}$ beam was scattered off the
transversely polarized NH3 target. These data provide an ideal opportunity to
verify the sign change of T-odd PDFs through simultaneous analysis of the
weighted Sivers asymmetries in both SIDIS and Drell-Yan processes.

In this work, we first numerically calculate the weighted
$\frac{P_{hT}}{zM_{p}}\sin(\phi_{h}-\phi_{S})$ asymmetry in charged hadron
produced SIDIS process with a lepton beam scattering off the transversely
polarized proton target at COMPASS and compare with the corresponding
experimental data. Due to the agreement between the theoretical calculation
and the experimental data in SIDIS process, then we change the sign of the
Sivers function applied in the SIDIS process to estimate the weighted
$\frac{q_{T}}{M_{p}}\sin(\phi_{S})$ asymmetry in Drell-Yan process with
$\pi^{-}$ beam scattering off the transversely polarized proton at COMPASS.

The rest of the paper is organized as follows. In Sec. II, we provide the
theoretical expressions of the weighted Sivers asymmetries in charged hadron
produced SIDIS process and $\pi^{-}p$ Drell-Yan process. In Sec. III, we
numerically estimate the weighted
$\frac{P_{hT}}{zM_{p}}\sin(\phi_{h}-\phi_{S})$ asymmetry in SIDIS process and
the weighted $\frac{q_{T}}{M_{p}}\sin(\phi_{S})$ asymmetry in $\pi^{-}p$
Drell-Yan process at COMPASS. We summarize the paper and discuss the results
in Sec. IV.

##  II FORMALISM OF THE WEIGHTED SIVERS ASYMMETRIES

###  II.1  Weighted Sivers asymmetry in SIDIS process

The process under study is the SIDIS process, particularly an unpolarized
lepton beam scatters off the transversely polarized proton target to produce
an unpolarized charged hadron in the final-state:

$l(\ell)+p^{\uparrow}(P_{p})\longrightarrow
l(\ell^{\prime})+h^{\pm}(P_{h})+X(P_{X}),$ (1)

where $\ell$ and $\ell^{\prime}$ denote the four-momenta of the incoming and
outgoing leptons, respectively; $P_{p}$ and $P_{h}$ represent the four-momenta
of the proton target $p$ and the final state hadron $h^{\pm}$, respectively;
$\uparrow$ denotes that the proton target is transversely polarized.
$q=\ell-\ell^{\prime}$ denotes the four-momentum of the virtual photon with
invariant mass squared $Q^{2}=-q^{2}$. The reference frame of the process is
depicted in Fig. 1, in which the momentum direction of the virtual photon
defines the $z$ axis in accordance with the Trento conventions
Bacchetta:2004jz . $P_{hT}$ and $S_{T}$ are the component of $P_{h}$ and the
proton spin vector $S$ that is transverse to the virtual photon momentum
direction ($z$ axis). $\phi_{h}$ denotes the azimuthal angle between the
lepton and hadron planes which has the definition in the target rest frame as

$\displaystyle\cos\phi_{h}=\frac{(\bm{\hat{q}}\times\bm{\ell})}{|\bm{\hat{q}}\times\bm{\ell}|}\cdot\frac{(\bm{\hat{q}}\times\bm{P_{h}})}{|\bm{\hat{q}}\times\bm{P_{h}}|}.$
(2)

The azimuthal angle in any frames can be obtained by a boost along
$\bm{\hat{q}}$, therefore, one can write the Lorentz invariant form of Eq. (2)
into

$\displaystyle\cos\phi_{h}=-\frac{g^{\mu\nu}_{\perp}l_{\mu}P_{h\nu}}{|l_{T}||P_{hT}|},\quad\sin\phi_{h}=-\frac{\epsilon^{\mu\nu}_{\perp}l_{\mu}P_{h\nu}}{|l_{T}||P_{hT}|},$
(3)

where $|l_{T}|=\sqrt{g^{\mu\nu}_{\perp}l_{\mu}l_{\nu}}$,
$|P_{hT}|=\sqrt{g^{\mu\nu}_{\perp}P_{h\mu}P_{h\nu}}$. $g^{\mu\nu}_{\perp}$ and
${\epsilon^{\mu\nu}_{\perp}}$ are the perpendicular projection tensors

$\displaystyle
g^{\mu\nu}_{\perp}=g^{\mu\nu}-\frac{q^{\mu}P_{p}^{\nu}+P_{p}^{\mu}q^{\nu}}{P_{p}\cdot
q(1+\gamma^{2})}+\frac{\gamma^{2}}{1+\gamma^{2}}(\frac{q^{\mu}q^{\nu}}{Q^{2}}-\frac{P_{p}^{\mu}P_{p}^{\nu}}{M_{p}^{2}}),$
(4)
$\displaystyle\epsilon^{\mu\nu}_{\perp}=\epsilon^{\mu\nu\rho\sigma}\frac{P_{p\rho}q_{\sigma}}{P_{p}\cdot
q\sqrt{(1+\gamma^{2})}},$ (5)

with $\gamma=\frac{2M_{p}x}{Q}$, $M_{p}$ being the proton mass and $x$ being
the Bjorken variable $x=\frac{Q^{2}}{2P_{p}\cdot q}$. $\phi_{S}$ stands for
the azimuthal angle between the direction of $S_{T}$ and the lepton scattering
plane which has the similar definition in Eqs. (2) and (3) with $P_{h}$
replaced by the covariant spin vector $S$ of the target proton.

Figure 1: Reference frame of the transversely polarized SIDIS process.

Besides the Bjorken variable $x$ and $\gamma$, the following Lorentz
invariants are introduced to express the differential cross section of the
SIDIS process as well as the experimental observables

$\displaystyle y=\frac{P_{p}\cdot q}{P_{p}\cdot\ell},\quad z=\frac{P_{p}\cdot
P_{h}}{P_{p}\cdot q}\,,\quad s=(P_{p}+\ell)^{2}\,,$ (6)

where $s$ is the center of mass energy squared of the $e$-$p$ system, $y$ is
the inelasticity, $z$ is the longitudinal momentum fraction of the final-state
hadron. The variables $x$ and $y$ are related to $Q^{2}$ through
$xy=\frac{Q^{2}}{s}$ (actually, $xy=\frac{Q^{2}}{s-m_{\ell}^{2}-M_{p}^{2}}$,
while the lepton mass and proton mass are neglected in general).

Assuming single-photon exchange, the SIDIS cross section can be expressed in
terms of 18 structure functions Bacchetta:2006tn  in a model-independent way.
Here, we only consider the terms relevant to the Sivers function, with other
terms absorbed into the ellipsis Bacchetta:2006tn ; Bacchetta:2010si

$\displaystyle\frac{d\sigma}{dxdydzd\phi_{S}d\phi_{h}d\bm{P}^{2}_{hT}}=\frac{\alpha^{2}}{xyQ^{2}}\,\frac{y^{2}}{2(1-\epsilon)}\left(1+\frac{\gamma^{2}}{2x}\right)\,\biggl{[}F_{UU}+|\bm{S}_{T}|\,\sin(\phi_{h}-\phi_{S})\,F_{UT}^{\sin(\phi_{h}-\phi_{S})}+\ldots\biggr{]},$
(7)

where $\alpha=e^{2}/(4\pi)$ is the fine structure constant, and $\epsilon$ is
the ratio of the longitudinal and transverse photon flux
$\epsilon=\frac{1-y-\frac{1}{4}\gamma^{2}y^{2}}{1-y+\frac{1}{2}y^{2}+\frac{1}{4}\gamma^{2}y^{2}}$
such that the depolarization factor as
$\frac{y^{2}}{2(1-\epsilon)}\approx(1-y+\frac{1}{2}y^{2})$. $F_{UU}$ stands
for the unpolarized structure function and
$F_{UT}^{\sin\left(\phi_{h}-\phi_{S}\right)}$ represent the transverse spin-
dependent structure function, with the subscripts $U$ (unpolarized) or $T$
(transverse polarized) denoting the polarization states of the beam (first
subscript) and the target (second subscript). The structure functions can be
expressed as the following convolutions of the TMD PDFs and FFs
Bacchetta:2006tn :

$\displaystyle F_{UU}$ $\displaystyle=$ $\displaystyle{\cal
C}\biggl{[}f_{1}D_{1}\biggr{]}\;,$ (8) $\displaystyle
F_{UT}^{\sin(\phi_{h}-\phi_{S})}$ $\displaystyle=$ $\displaystyle{\cal
C}\biggl{[}-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M_{p}}\,f_{1T}^{\perp}\,D_{1}\biggr{]}\;$
(9)

with the notation $\mathcal{C}$ representing the convolution of the transverse
momentum:

$\mathcal{C}\bigl{[}\omega fD\bigr{]}=x\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\bigl{(}z\bm{p}_{T}+\bm{k}_{T}-\bm{P}_{hT}\bigr{)}\omega(\bm{p}_{T},-\bm{k}_{T}/z)f^{q}(x,\bm{p}_{T}^{2})\,D^{q}(z,\bm{k}_{T}^{2}),$
(10)

where the sum runs over all the quark and antiquark flavors,
$\hat{\bm{h}}=\bm{P}_{hT}/|\bm{P}_{hT}|$, $\bm{p}_{T}$ and $\bm{k}_{T}$ are
the transverse momentum of quarks in the target proton and that of the final-
state hadron relative to the fragmentation quark, respectively; $\omega$ is a
function of $\bm{p}_{T}$ and $\bm{k}_{T}$, $f_{1}(x,\bm{p}_{T}^{2})$ and
$f_{1T}^{\perp}(x,\bm{p}_{T}^{2})$ are the unpolarized TMD PDF and the Sivers
function, respectively; $D_{1}(z,\bm{k}_{T}^{2})$ is the unpolarized FF.

Since the evolution effects of the TMD PDF and FF encoded by the Collins-Soper
equation are complicated and can not be analytically solved without the
parametrization of the nonperturbative evolution kernel Collins:2011zzd , for
simplicity, we assume Gaussian form for TMD PDF and FF Anselmino:2008jk ;
Anselmino:2007fs ; Anselmino:2013vqa ; Bradamante:2017yia ; Lefky:2014eia :

$\displaystyle f(x,\bm{p}_{T}^{2})=f(x)\,\frac{1}{\pi\langle
p_{T}^{2}\rangle}\,\exp\left(-\frac{\bm{p}_{T}^{2}}{\langle
p_{T}^{2}\rangle}\right)\;,$ $\displaystyle
D(z,\bm{k}_{T}^{2})=D(z)\,\frac{1}{\pi\langle
k_{T}^{2}\rangle}\,\exp\left(-\frac{\bm{k}_{T}^{2}}{\langle
k_{T}^{2}\rangle}\right)\;.$ (11)

Here, $f(x)$ and $D(z)$ are the collinear PDF and FF that depend on $Q^{2}$,
that is, they follow the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP)
evolution. $\langle p_{T}^{2}\rangle$ and $\langle k_{T}^{2}\rangle$ are the
average values of $\bm{p}_{T}^{2}$ and $\bm{k}_{T}^{2}$, which will be
eliminated in the integral calculation over the intrinsic transverse momenta.
This Gaussian form of TMD PDF and FF is at tree level and is suitable to
describe nonperturbative effects at small transverse momenta.

In order to simplify the complicated TMD evolution effects as well as the
convolution in the transverse momentum space, we will consider the transverse
momentum weighted asymmetry, which can be expressed as the simple product of
the collinear PDFs and FFs. Following Refs. Boer:1997nt ; Bacchetta:2010si ,
weighted asymmetry with proper weighting function $W$ is defined as:

$\displaystyle A_{UT}^{W}=\frac{\int d^{2}\bm{P}_{hT}WF_{UT}}{\int
d^{2}\bm{P}_{hT}F_{UU}}.$ (12)

After choosing the weighting function $W=\frac{P_{hT}}{zM_{p}}$, the weighted
Sivers asymmetry can be written as

$\displaystyle
A_{UT}^{\sin(\phi_{h}-\phi_{S})\,\frac{P_{hT}}{zM_{p}}}=\frac{\int
d^{2}\bm{P}_{hT}\frac{P_{hT}}{zM_{p}}F_{UT}^{\sin\left(\phi_{h}-\phi_{S}\right)}}{\int
d^{2}\bm{P}_{hT}F_{UU}}.$ (13)

with the denominator $F_{UU}$ being

$\displaystyle\int d^{2}\bm{P}_{hT}F_{UU}=$ $\displaystyle\int
d^{2}\bm{P}_{hT}\,{\cal C}\bigg{[}f_{1}\;D_{1}\bigg{]}$ $\displaystyle=$
$\displaystyle\ x\sum_{q}e_{q}^{2}\int
d^{2}\bm{P}_{hT}d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\bigl{(}z\bm{p}_{T}+\bm{k}_{T}-\bm{P}_{hT}\bigr{)}f_{1}(x,\bm{p}_{T}^{2})D_{1}(z,\bm{k}_{T}^{2})$
$\displaystyle=$
$\displaystyle\sum_{q}e_{q}^{2}\left[xf_{1}(x)\,D_{1}(z)\right].$ (14)

Similarly, the $P_{hT}$-weighted spin-dependent structure function
$\frac{P_{hT}}{zM_{p}}F_{UT}^{\sin\left(\phi_{h}-\phi_{S}\right)}$ can be
written as

$\displaystyle\int
d^{2}\bm{P}_{hT}\frac{P_{hT}}{zM_{p}}F_{UT}^{\sin\left(\phi_{h}-\phi_{S}\right)}=$
$\displaystyle\int d^{2}\bm{P}_{hT}\frac{P_{hT}}{zM_{p}}{\cal
C}\biggl{[}-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M_{p}}\,f_{1T}^{\perp}\,D_{1}\biggr{]}$
$\displaystyle=$ $\displaystyle~{}x\sum_{q}e_{q}^{2}\int
d^{2}\bm{P}_{hT}\frac{P_{hT}}{zM_{p}}d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\bigl{(}z\bm{p}_{T}+\bm{k}_{T}-\bm{P}_{hT}\bigr{)}\frac{-\hat{\bm{h}}\cdot\bm{p}_{T}}{M_{p}}f_{1T}^{\perp}(x,\bm{p}_{T}^{2})D_{1}(z,\bm{k}_{T}^{2})$
$\displaystyle=$
$\displaystyle~{}2\sum_{q}e_{q}^{2}\bigl{[}xf_{1T}^{\perp(1)}(x)\,D_{1}(z)\bigr{]},$
(15)

where $f_{1T}^{\perp(1)}(x)$ is the first $\bm{p}_{T}$-moment of the Sivers
function:

$\displaystyle
f_{1T}^{\perp(1)}(x)=\int{d}^{2}\bm{p}_{T}\frac{\bm{p}_{T}^{2}}{2M_{p}^{2}}f_{1T}^{\perp}(x,\bm{p}_{T}^{2}).$
(16)

In SIDIS process, $P_{hT}$-weighted Sivers asymmetry with weighting function
$W=\frac{P_{hT}}{zM_{p}}$ can be rewritten as

$\displaystyle
A_{UT}^{\sin(\phi_{h}-\phi_{S})\,\frac{P_{hT}}{zM_{p}}}=2\,\frac{\sum_{q}e_{q}^{2}\bigl{[}xf_{1T}^{\perp(1)}(x)\,D_{1}(z)\bigr{]}}{\sum_{q}e_{q}^{2}\bigl{[}xf_{1}(x)\,D_{1}(z)\bigr{]}}\;.$
(17)

###  II.2  Weighted Sivers asymmetry in Drell-Yan process

In this subsection, we will set up the theoretical framework of the transverse
momentum weighted Sivers asymmetry in Drell-Yan process in which a $\pi^{-}$
beam scatters off a transversely polarized proton target:

$\displaystyle\pi^{-}(P_{\pi})+p^{\uparrow}(P_{p})\longrightarrow\gamma^{*}(q)+X(P_{X})\longrightarrow
l^{+}(\ell)+l^{-}(\ell^{\prime})+X(P_{X}),$ (18)

where $P_{\pi}$, $P_{p}$, and $q$ represent the four-momenta of the pion beam
$\pi^{-}$, the proton target $p$, and the virtual photon $\gamma^{*}$,
respectively. One should note that $q$ is a time-like vector different from
that in the SIDIS process, namely, $Q^{2}=q^{2}>0$. $Q^{2}$ is the invariant
mass square of the final-state dilepton. The target rest frame of the process
under study is shown in Fig. 2, in which the target proton is at rest and the
momentum direction of the incident pion defines $z$ axis. $q_{T}$ stands for
the transverse component of $q$. $S_{T}$ is the transverse spin of the target
proton, and $\phi_{S}$ denotes the angle between the direction of $S_{T}$ and
the hadron plane.

Figure 2: Target rest frame of the $\pi^{-}p$ Drell-Yan process.

The following Lorentz invariants are introduced to express the differential
cross section of the Drell-Yan process as well as the experimental
observables:

$\displaystyle x_{\pi}=\frac{Q^{2}}{2P_{\pi}\cdot q},\quad
x_{p}=\frac{Q^{2}}{2P_{p}\cdot q},\quad x_{F}=2q_{L}/s=x_{\pi}-x_{p},$
$\displaystyle s=(P_{\pi}+P_{p})^{2},\quad\tau=Q^{2}/s=x_{\pi}x_{p},\quad
y=\frac{1}{2}\mathrm{ln}\frac{q^{+}}{q^{-}}=\frac{1}{2}\mathrm{ln}\frac{x_{\pi}}{x_{p}},$
(19)

where $s$ represents the center of mass energy squared of the $\pi^{-}p$
system. $x_{\pi}$ and $x_{p}$ are the Bjorken variables of the pion and
proton, respectively; $x_{F}$ is the Feynman $x$ variable with $q_{L}$ being
the longitudinal momentum of the virtual photon; $y$ is the rapidity of the
dilepton.

At leading order, the cross section of the transversely polarized $\pi^{-}p$
Drell-Yan process can be written as Arnold:2008kf ; COMPASS:2010shj :

$\displaystyle\frac{d\sigma}{dq^{4}d\Omega}=\frac{\alpha^{2}_{em}}{Fq^{2}}\,$
$\displaystyle\biggl{[}(1+\cos^{2}\theta)F_{UU}+|\bm{S}_{T}|(1+\cos^{2}\theta)\sin\phi_{S}\,F_{UT}^{\sin\phi_{S}}+\cdots\biggr{]},$
(20)

where the solid angle $\Omega$ specifies the orientation of the dilepton. $F$
represents the flux of incoming hadrons. $F_{UU}$ and $F_{UT}^{\sin\phi_{S}}$
are the relevant structure functions. Using the notation $\cal C$,

$\displaystyle\mathcal{C}\bigl{[}\omega\,f_{q/\pi}\,f_{\bar{q}/p}\bigr{]}=\frac{1}{N_{c}}\;\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T\pi}\,d^{2}\bm{p}_{Tp}\,\delta^{(2)}(\bm{q}_{T}-\bm{p}_{T\pi}-\bm{p}_{Tp})\,\omega\,f_{q/\pi}(x_{\pi},\bm{p}_{T\pi}^{2})f_{\bar{q}/p}(x_{p},\bm{p}_{Tp}^{2})\,,$
(21)

the two structure functions can be expressed as the following convolutions
Arnold:2008kf :

$\displaystyle F_{UU}={\cal C}\bigg{[}f_{1,q/\pi}\;f_{1,\bar{q}/p}\bigg{]},$
(22) $\displaystyle F_{UT}^{\sin\phi_{S}}={\cal
C}\bigg{[}\frac{\hat{\bm{h}}\cdot\bm{p}_{Tp}}{M_{p}}\;f_{1,q/\pi}\;f_{1T,\bar{q}/p}^{\perp}\;\bigg{]},$
(23)

where $\hat{\bm{h}}=\bm{q}_{T}/|\bm{q}_{T}|$. $M_{p}$ is the proton mass.
$f_{1,\bar{q}/p}(x_{p},\bm{p}_{Tp}^{2})$ and
$f_{1T,\bar{q}/p}^{\perp}(x_{p},\bm{p}_{Tp}^{2})$ are the unpolarized TMD PDF
and the Sivers function of proton, which depend on the Bjorken variable
$x_{p}$ and the transverse momentum $\bm{p}_{Tp}$ of the quark inside proton.
$f_{1,q/\pi}(x_{\pi},\bm{p}_{T\pi}^{2})$ is the unpolarized TMD PDF of pion
that depend on $x_{\pi}$ and $\bm{p}_{T\pi}$.

After choosing the weighting function $W=\frac{q_{T}}{M_{p}}$, the weighted
Sivers asymmetry in Drell-Yan process can be written as:

$\displaystyle A_{UT}^{\sin(\phi_{S})\frac{q_{T}}{M_{p}}}=\frac{\int
d^{2}\bm{q}_{T}\frac{q_{T}}{M_{p}}F_{UT}^{\sin\phi_{S}}}{\int
d^{2}\bm{q}_{T}F_{UU}}.$ (24)

with the denominator being

$\displaystyle\int d^{2}\bm{q}_{T}F_{UU}=$ $\displaystyle\int
d^{2}\bm{q}_{T}\,{\cal C}\bigg{[}f_{1,q/\pi}\;f_{1,\bar{q}/p}\bigg{]}$
$\displaystyle=$ $\displaystyle\frac{1}{N_{c}}\sum_{q}e_{q}^{2}\int
d^{2}\bm{q}_{T}\,d^{2}\bm{p}_{T\pi}\,d^{2}\bm{p}_{Tp}\,\delta^{(2)}(\bm{q}_{T}-\bm{p}_{T\pi}-\bm{p}_{Tp})\,f_{1,q/\pi}(x_{\pi},\bm{p}_{T\pi}^{2})f_{1,\bar{q}/p}(x_{p},\bm{p}_{Tp}^{2})$
$\displaystyle=$
$\displaystyle\frac{1}{N_{c}}\sum_{q}e_{q}^{2}\bigl{[}f_{1,q/\pi}(x_{\pi})\,f_{1,\bar{q}/p}(x_{p})\bigr{]}.$
(25)

For the spin dependent structure function
$F_{UT}^{\sin\left(\phi_{S}\right)}$, we have

$\displaystyle\int d^{2}\bm{q}_{T}\frac{q_{T}}{M_{p}}F_{UT}^{\sin\phi_{S}}=$
$\displaystyle\int d^{2}\bm{q}_{T}\frac{q_{T}}{M_{p}}\,{\cal
C}\bigg{[}\frac{\hat{\bm{h}}\cdot\bm{p}_{Tp}}{M_{p}}\;f_{1,q/\pi}\;f_{1T,\bar{q}/p}^{\perp}\;\bigg{]}$
$\displaystyle=$ $\displaystyle\frac{1}{N_{c}}\sum_{q}e_{q}^{2}\int
d^{2}\bm{q}_{T}\,\frac{q_{T}}{M_{p}}\,d^{2}\bm{p}_{T\pi}\,d^{2}\bm{p}_{Tp}\,\delta^{(2)}(\bm{q}_{T}-\bm{p}_{T\pi}-\bm{p}_{Tp})\,\frac{\hat{\bm{h}}\cdot\bm{p}_{Tp}}{M_{p}}\;f_{1,q/\pi}(x_{\pi},\bm{p}_{T\pi}^{2})f^{\perp}_{1T,\bar{q}/p}(x_{p},\bm{p}_{Tp}^{2})$
$\displaystyle=$ $\displaystyle
2\,\frac{1}{N_{c}}\sum_{q}e_{q}^{2}\bigl{[}f_{1,q/\pi}(x_{\pi})\,f_{1T,\bar{q}/p}^{\perp(1)}(x_{p})\bigr{]}.$
(26)

In Drell-Yan process, the $q_{T}$-weighted Sivers asymmetry with weighting
function $W=\frac{q_{T}}{M_{p}}$ can be rewritten as:

$\displaystyle A_{UT}^{\sin(\phi_{S})\frac{q_{T}}{M_{p}}}$
$\displaystyle=2\frac{\sum_{q}e_{q}^{2}\bigl{[}f_{1,q/\pi}(x_{\pi})\,f_{1T,\bar{q}/p}^{\perp(1)}(x_{p})\bigr{]}}{\sum_{q}e_{q}^{2}\bigl{[}f_{1,q/\pi}(x_{\pi})\,f_{1,\bar{q}/p}(x_{p})\bigr{]}}.$
(27)

##  III NUMERICAL CALCULATION

Based on the above formalism, we present the numerical calculation for the
weighted Sivers asymmetries in charged hadron produced SIDIS process and the
$\pi^{-}p$ Drell-Yan process at the kinematical regions of COMPASS, and
compare the numerical results with the SIDIS data COMPASS:2018ofp  and Drell-
Yan preliminary data Longo:2019bih .

To do this, the collinear unpolarized PDF of the proton $f_{1,p}(x_{p})$ and
unpolarized FF of the final-state hadron $D_{1}^{q}(z)$ are adopted from the
CT10 parametrization Lai:2010vv  and the DSS parametrization deFlorian:2007aj
, respectively. For the unpolarized PDF of the pion $f_{1,\pi}(x_{\pi})$, we
use the SMRS parametrization Sutton:1991ay .

As mentioned above, one needs the information of the first $\bm{p}_{T}$-moment
of proton Sivers function $f_{1T}^{\perp(1)}(x)$ in both SIDIS and Drell-Yan
processes. $f_{1T}^{\perp(1)}(x)$ and the usual Qiu-Sterman function
$T_{q,F}(x,x)$ have the following relation Sun:2013hua ; Kang:2011mr :

$\displaystyle
T_{q,F}(x,x)=\int{d}^{2}\bm{p}_{T}\frac{p_{T}^{2}}{M_{p}}f_{1T}^{\perp}(x,\bm{p}_{T}^{2})=2M_{p}f_{1T}^{\perp(1)}(x).$
(28)

In Ref. Echevarria:2014xaa , the authors have extracted Qiu-Sterman function
by using the preliminary SIDIS data from HERMES HERMES:2009lmz , COMPASS
COMPASS:2008isr ; COMPASS:2012dmt , and JLab JeffersonLabHallA:2011ayy  on
Sivers asymmetry. Assuming the Qiu-Sterman function is proportional to the
unpolarized PDF $f_{1,p}(x_{p})$, we adopt the parametrization
Echevarria:2014xaa  for $T_{q,F}(x,x)$ at the initial scale
$Q^{2}_{0}=2.4~{}\rm{GeV}^{2}$:

$\displaystyle
T_{q,F}(x,x)=N_{q}\frac{(\alpha_{q}+\beta_{q})^{(\alpha_{q}+\beta_{q})}}{\alpha_{q}^{\alpha_{q}}\beta_{q}^{\beta^{q}}}x^{\alpha_{q}}(1-x)^{\beta_{q}}f_{1,p}(x_{p}),$
(29)

and the values of the free parameters in global fit are Echevarria:2014xaa :

$\displaystyle\alpha_{u}=1.051_{-0.180}^{+0.192},\quad\alpha_{d}=1.552_{-0.275}^{+0.303},\quad\alpha_{sea}=0.851_{-0.305}^{+0.307},\quad\beta=4.857_{-1.395}^{+1.534},$
$\displaystyle N_{u}=0.106_{-0.009}^{+0.011},\quad
N_{d}=-0.163_{-0.046}^{+0.039},\quad N_{s}=0.103_{-0.604}^{+0.548},$
$\displaystyle N_{\bar{u}}=-0.012_{-0.020}^{+0.018},\quad
N_{\bar{d}}=-0.105_{-0.060}^{+0.043},\quad N_{\bar{s}}=-1.000\pm 1.757.$ (30)

The extraction can give a rather good description for all the corresponding
data of the differential cross sections. We assume that Qiu-Sterman function
follows the proportional relation with $f_{1,p}(x)$ in all the energy region.
Note that the Sivers function applied in Drell-Yan process has a opposite sign
with respect to the one extracted from SIDIS data.

The kinematical ranges in SIDIS process covered by COMPASS are as follows
COMPASS:2018ofp ; Bradamante:2017yia :

$\displaystyle 0.004<x<0.7,\quad 0.1<y<0.9,\quad z>0.2,$ $\displaystyle
Q^{2}>1\ \mathrm{GeV}^{2}\;,\quad 5\ \mathrm{GeV}<W<18\ \mathrm{GeV},\quad
s=301~{}\mathrm{GeV}^{2}.$ (31)

where $W^{2}=(P_{p}+q)^{2}\approx\frac{1-x}{x}Q^{2}$ is the invariant mass
squared of the virtual photon-proton system. As for the $\pi^{-}p$ Drell-Yan
process at COMPASS, we adopt the following kinematical cuts:

$\displaystyle 0.05<x_{p}<0.4,\quad 0.05<x_{\pi}<0.9,\quad-0.3<x_{F}<1,$
$\displaystyle 4.3\ \mathrm{GeV}<Q<8.5\ \mathrm{GeV},\quad
s=357~{}\mathrm{GeV}^{2}.$ (32)

We apply Eq. (17) to calculate the $\frac{P_{hT}}{zM_{p}}$-weighted Sivers
asymmetry of charged hadron produced in SIDIS process and plot the results in
Fig. 3. The upper two panels show the asymmetry for $h^{+}$ production; the
lower panels give the results for $h^{-}$ production. The left and right
panels show the asymmetry as functions of $x$ and $z$, respectively. In each
case the other kinematical variables are integrated out. In Fig. 3, the solid
lines correspond to the results from the central values of the parameters in
the parametrization of $T_{q,F}(x,x)$, and the shaded areas represent the
uncertainty bands due to the uncertainties of the parameters in Eq. (30). As a
comparison, we also show the experimental data measured by the COMPASS
Collaboration COMPASS:2018ofp . As shown in Fig. 3, our results are in
agreement with the COMPASS data within errors.

Figure 3: The weighted Sivers asymmetry of charged hadron produced in SIDIS
process. The solid squares with error bars represent the COMPASS data for
comparison COMPASS:2018ofp .

Figure 4: The weighted Sivers asymmetry in $\pi^{-}p$ Drell-Yan process. The
solid squares with error bars are COMPASS preliminary data Longo:2019bih .

The agreement between our results and the COMPASS SIDIS data encourages us to
further estimate the weighted Sivers asymmetry
$A_{UT}^{\frac{q_{T}}{M_{p}}}\sin(\phi_{S})$ in $\pi^{-}p$ Drell-Yan process
by utilizing the relation
$f_{1T}^{\perp\mathrm{(DY)}}=-f_{1T}^{\perp\mathrm{(SIDIS)}}$, which can be
used to verify the sign-change of Sivers function. The corresponding results
are plotted in Fig. 4, in which the four panels show the asymmetry as
functions of $x_{p}$ (upper left), $x_{\pi}$ (upper right), $x_{F}$ (lower
left) and the asymmetry integrated over the entire kinematical range (lower
right). In Fig. 4, the solid lines denote the central results, and the
uncertainty bands are determined by the uncertainties of the parameters in Eq.
(30). The corresponding preliminary COMPASS data with error bars Longo:2019bih
are also shown in Fig. 4. One can find the theoretical results of
$x_{\pi}$-dependent asymmetry agree with the preliminary data. Although there
is difference between the preliminary data and the theoretical calculation of
the tendency for the $x_{p}$-dependent asymmetry as well as the
$x_{F}$-dependent asymmetry, the asymmetry integrated over the entire
kinematical range still agrees with the preliminary COMPASS data within
errors. Compared to the shaded areas in Fig. 3, the uncertainty bands of
theoretical calculation for asymmetry in $\pi^{-}p$ Drell-Yan process are
relatively narrow. The reason may be that in Eq. (27) the contribution of sea
quarks is small, thus the asymmetry is dominated by
$2f_{1T,p}^{\perp(1)u(\mathrm{DY})}(x_{p})/f_{1,p}^{u}(x_{p})$. It indicates
that only the uncertainty of the parameters for $u$ quark in Eq. (30) is the
main source of the uncertainty bands, which means that the Sivers function of
sea quark has not been well constrained.

At present, we are still unable to conclude that the Sivers function has
opposite sign in SIDIS process and in Drell-Yan process due to the limited
Drell-Yan data with large errors. More and high precision experimental data
are needed to verify the opposite sign of the Sivers function predicted by QCD
and to constrain the sea quark Sivers function.

##  IV CONCLUSION

In this work, we have studied the weighted Sivers asymmetries in charged
hadron produced SIDIS process and $\pi^{-}p$ Drell-Yan process. The
asymmetries are contributed by the first transverse-moment of the Sivers
function, for which we have utilized the parameterization from literature. The
weighted $\frac{P_{hT}}{zM_{p}}\sin(\phi_{h}-\phi_{S})$ asymmetry in SIDIS
process and the weighted $\frac{q_{T}}{M_{p}}\sin(\phi_{S})$ asymmetry in
Drell-Yan process at the kinematics configurations of COMPASS SIDIS and Drell-
Yan programs have been estimated, respectively. Our numerical estimate of the
asymmetry in charged hadron produced SIDIS process is consistent with the
COMPASS measurement. However, due to the preliminary Drell-Yan data with large
errors, our numerical results are still unable to demonstrate clearly the
sign-change property of the Sivers function, and high precision experimental
data are needed to clarify this point and to constrain the sea quark Sivers
function.

##  V ACKNOWLEDGMENTS

This work is partially supported by the NSFC (China) grants 11905187,11847217,
12150013. X. Wang is supported by the China Postdoctoral Science Foundation
under Grant No. 2018M640680.

## References

  * (1) J. Ashman et al. [European Muon], Phys. Lett. B 206, 364 (1988)
  * (2) J. Ashman et al. [European Muon], Nucl. Phys. B 328, 1 (1989)
  * (3) D. W. Sivers, Phys. Rev. D 41, 83 (1990).
  * (4) S. J. Brodsky, D. S. Hwang and I. Schmidt, Phys. Lett. B 530, 99 (2002).
  * (5) S. J. Brodsky, D. S. Hwang and I. Schmidt, Nucl. Phys. B642, 344 (2002).
  * (6) D. Boer, S. J. Brodsky and D. S. Hwang, Phys. Rev. D 67, 054003 (2003).
  * (7) J. C. Collins, Phys. Lett. B 536, 43-48 (2002).
  * (8) A. Airapetian et al. [HERMES], Phys. Rev. Lett. 94, 012002 (2005).
  * (9) A. Airapetian et al. [HERMES], Phys. Rev. Lett. 103, 152002 (2009).
  * (10) A. Airapetian et al. [HERMES], J. High Energy Phys. 12 (2020) 010.
  * (11) M. Alekseev et al. [COMPASS], Phys. Lett. B 673, 127-135 (2009).
  * (12) M. G. Alekseev et al. [COMPASS], Phys. Lett. B 692, 240-246 (2010).
  * (13) C. Adolph et al. [COMPASS], Phys. Lett. B 717, 383-389 (2012).
  * (14) C. Adolph et al. [COMPASS], Phys. Lett. B 770, 138-145 (2017).
  * (15) X. Qian et al. [Jefferson Lab Hall A], Phys. Rev. Lett. 107, 072003 (2011).
  * (16) Y. X. Zhao et al. [Jefferson Lab Hall A], Phys. Rev. C 90, 055201 (2014).
  * (17) M. Aghasyan et al. [COMPASS], Phys. Rev. Lett. 119, 112002 (2017).
  * (18) L. Adamczyk et al. [STAR], Phys. Rev. Lett. 116, 132301 (2016).
  * (19) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia and A. Prokudin, Phys. Rev. D 72, 094007 (2005) [erratum: Phys. Rev. D 72, 099903 (2005)].
  * (20) A. V. Efremov, K. Goeke, S. Menzel, A. Metz and P. Schweitzer, Phys. Lett. B 612, 233-244 (2005).
  * (21) J. C. Collins, A. V. Efremov, K. Goeke, S. Menzel, A. Metz and P. Schweitzer, Phys. Rev. D 73, 014021 (2006).
  * (22) W. Vogelsang and F. Yuan, Phys. Rev. D 72, 054028 (2005).
  * (23) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, S. Melis, F. Murgia, A. Prokudin and C. Turk, Eur. Phys. J. A 39, 89 (2009).
  * (24) M. Anselmino, M. Boglione and S. Melis, Phys. Rev. D 86, 014028 (2012).
  * (25) A. Bacchetta and M. Radici, Phys. Rev. Lett. 107, 212001 (2011).
  * (26) M. G. Echevarria, A. Idilbi, Z. B. Kang and I. Vitev, Phys. Rev. D 89, 074013 (2014).
  * (27) M. Anselmino, M. Boglione, U. D’Alesio, F. Murgia and A. Prokudin, J. High Energy Phys. 04 (2017) 046.
  * (28) A. Martin, F. Bradamante and V. Barone, Phys. Rev. D 95, 094024 (2017).
  * (29) M. Boglione, U. D’Alesio, C. Flore and J. O. Gonzalez-Hernandez, J. High Energy Phys. 07 (2018) 148 .
  * (30) M. Bury, A. Prokudin and A. Vladimirov, J. High Energy Phys. 05 (2021) 151 .
  * (31) A. Bacchetta, A. Schaefer and J. J. Yang, Phys. Lett. B 578, 109-118 (2004).
  * (32) Z. Lu and B. Q. Ma, Phys. Rev. D 70, 094044 (2004).
  * (33) B. Pasquini and F. Yuan, Phys. Rev. D 81, 114013 (2010).
  * (34) A. Courtoy, F. Fratini, S. Scopetta and V. Vento, Phys. Rev. D 78, 034002 (2008).
  * (35) F. Yuan, Phys. Lett. B 575, 45-54 (2003).
  * (36) A. Courtoy, S. Scopetta and V. Vento, Phys. Rev. D 79, 074001 (2009).
  * (37) X. Luan and Z. Lu, arXiv:2204.06854 [hep-ph].
  * (38) A. M. Kotzinian and P. J. Mulders, Phys. Rev. D 54, 1229-1232 (1996).
  * (39) D. Boer and P. J. Mulders, Phys. Rev. D 57, 5780-5786 (1998).
  * (40) B. Pasquini and P. Schweitzer, Phys. Rev. D 83, 114044 (2011).
  * (41) M. Anselmino, H. Avakian, D. Boer, F. Bradamante, M. Burkardt, J. P. Chen, E. Cisbani, M. Contalbrigo, D. Crabb and D. Dutta, et al. Eur. Phys. J. A 47, 35 (2011).
  * (42) S. C. xue, X. Wang, D. M. Li and Z. Lu, Phys. Lett. B 820, 136598 (2021).
  * (43) A. Bacchetta, M. Radici, F. Conti and M. Guagnelli, Eur. Phys. J. A 45, 373-388 (2010).
  * (44) A. Sissakian, O. Shevchenko, A. Nagaytsev, O. Denisov and O. Ivanov, Eur. Phys. J. C 46, 147-150 (2006).
  * (45) Z. Lu, B. Q. Ma and J. She, Phys. Lett. B 696, 513-517 (2011).
  * (46) Z. Wang, X. Wang and Z. Lu, Phys. Rev. D 95, 094004 (2017).
  * (47) X. Liu and B. Q. Ma, Eur. Phys. J. C 81, 635 (2021).
  * (48) M. G. Alexeev et al. [COMPASS], Nucl. Phys. B940, 34-53 (2019).
  * (49) R. Longo [COMPASS], PoS DIS2019, 186 (2019).
  * (50) A. Bacchetta, U. D’Alesio, M. Diehl and C. A. Miller, Phys. Rev. D 70, 117504 (2004).
  * (51) A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders and M. Schlegel, J. High Energy Phys. 02 (2007) 093.
  * (52) J. Collins, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 32, 1-624 (2011).
  * (53) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia, A. Prokudin and S. Melis, Nucl. Phys. Proc. Suppl. 191, 98 (2009).
  * (54) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia, A. Prokudin and C. Turk, Phys. Rev. D 75, 054032 (2007).
  * (55) M. Anselmino, M. Boglione, U. D’Alesio, S. Melis, F. Murgia and A. Prokudin, Phys. Rev. D 87, 094019 (2013).
  * (56) F. Bradamante [COMPASS], arXiv:1702.00621 [hep-ex].
  * (57) C. Lefky and A. Prokudin, Phys. Rev. D 91, 034010 (2015).
  * (58) S. Arnold, A. Metz and M. Schlegel, Phys. Rev. D 79, 034005 (2009).
  * (59) F. Gautheron et al. [COMPASS], SPSC-P-340.
  * (60) H. L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin and C. P. Yuan, Phys. Rev. D 82, 074024 (2010).
  * (61) D. deFlorian, R. Sassot and M. Stratmann, Phys. Rev. D 75, 114010 (2007).
  * (62) P. J. Sutton, A. D. Martin, R. G. Roberts and W. J. Stirling, Phys. Rev. D 45, 2349-2359 (1992).
  * (63) P. Sun and F. Yuan, Phys. Rev. D 88, 114012 (2013).
  * (64) Z. B. Kang, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 107, 152002 (2011).


# Warmstarting column generation for unit commitment

Nagisa Sugishita, Andreas Grothey††footnotemark: , Ken McKinnon††footnotemark:
School of Mathematics, University of Edinburgh, James Clerk Maxwell Building,
Edinburgh (UK), EH9 3FD, Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This work focuses on a solution method for unit commitment problems which are
to be solved repeatedly with different data but with the same problem
structures. Dantzig-Wolfe decomposition with a column generation procedure has
been proved successful for the unit commitment problem. We demonstrate use of
a neural network to generate initial dual values for the column generation
procedure, and show how the neural network can be trained efficiently using
dual decomposition. Our numerical experiments compare our methods with
baselines and reveal that our approach solves test instances to high precision
in shorter computational time and scales to handle larger-scale instances.

##  1 Introduction

The unit commitment (UC) problem is an important optimization problem in the
energy industry to compute the optimal operating schedules of power plants
given demand over a fixed time period. This problem is solved by electricity
generating companies on a daily basis to determine which generators are to be
used. Given a demand forecast over a planning horizon, the timing of switching
and the amount of power dispatch are to be optimized simultaneously. This
gives rise to a large-scale combinatorial problem and due to its practical
importance it has been extensively studied over the last few decades. For a
recent survey, see [1].

This work focuses on UC problems which are to be solved repeatedly with
different data but with the same problem structures. In practice, when a UC
problem is solved as a day-ahead planning problem, the characterizations of
generators such as generation costs and ramping rates remain the same and
problems are solved with different demand forecasts.

One popular solution method for UC problems is to use Dantzig-Wolfe
decomposition to decompose the problem by generators. The reformulated problem
is then solved with a column generation procedure. This procedure can be seen
as the dual of a cutting plane approach to Lagrangian relaxation, as discussed
by [5]. The DW reformulation is a relaxation of the original problem, however
[3] and [2] reported that the optimality gap introduced by the reformulation
is typically small especially if the problem size is large. In such cases,
Dantzig-Wolfe decomposition with a suitable primal heuristic is likely to
solve the UC problem to high precision.

In the aforementioned decomposition-based approaches, the dual variable plays
a key role. It is expected that using an appropriate dual value as an initial
point may speed up the solution method. One approach to generate such a dual
value is to solve an approximation of the original problem. [11] solved the
linear programming relaxation (LPR) of the original problem and obtained its
optimal dual value. It was then used to initialize Dantzig-Wolfe
decomposition. [12] constructed an approximation by dropping some constraints,
such as integrality constraints, minimum up/down time constraints and minimum
power output constraints. The optimal dual value of the approximation was used
to warmstart Lagrangian relaxation with the subgradient algorithm. However, as
[11] reported, solving such an approximation may require a non-trivial amount
of time.

To warmstart the column generation procedure, we propose to use a pre-trained
model to generate initial dual values. Our work is inspired by the method for
generating dual policies proposed by [10]. In their work, neural networks were
trained to compute dual values that yielded strong lower bounds on two-stage
integer stochastic programming. In this work, we propose to combine the dual
policy with the aforementioned column generation procedure to solve UC.
Namely, a model is used with problem data as input to generate dual values
that yield tight lower bounds, as demonstrated by [10]. In our work, the
generated dual values are then used to warmstart the column generation
procedure. With this approach, we can exploit the strength of the dual policy
while maintaining the desirable property of the solution method, such as
exactness. The above approach were tested with numerical experiments on large-
scale UC problems, where our method outperformed baselines.

The rest of the paper is structured as follows. Section 2 briefly introduces
Dantzig-Wolfe decomposition and the column generation procedure. Section 3
presents approaches and issues to warmstart the column generation procedure. A
model based on a neural network is considered to generate an initial value for
the algorithm. The proposed approach is tested on large scale UC problems in
Section 4. Finally, in Section 5 conclusions and further extensions of this
work are discussed.

##  2 Dantzig-Wolfe decomposition and column generation

In this section, we briefly review Dantzig-Wolfe decomposition and the column
generation procedure. For further background, see [14].

Consider the following family of mixed-integer programmes parametrized by
$\omega$:

$\displaystyle z(\omega)=\min_{x}\ $ $\displaystyle\sum_{s\in
S}c_{s}^{T}x_{s}$ (2.1) s.t. $\displaystyle\sum_{s\in S}A_{s}x_{s}=a(\omega),$
$\displaystyle x_{s}\in X_{s}(\omega)\quad\forall s\in S,$

where $x=(x_{s})_{s\in S}$ is a vector of decision variables and
$X_{s}(\omega)$ is the set of feasible points given by

$X_{s}(\omega)=\\{x_{s}\in\mathbb{Z}^{n}\times\mathbb{R}^{m}\mid Dx_{s}\leq
d(\omega)\\}\quad s\in S.$

For simplicity, we assume that $X_{s}(\omega)$ is non-empty and bounded for
every $\omega$ and $s\in S$, and that the problem (2.1) has an optimal
solution for every $\omega$. To reduce clutter in what follows we drop the
dependence on $\omega$ unless it causes confusion.

In Dantzig-Wolfe decomposition, we consider a relaxation of (2.1), referred to
as the master problem (MP), by replacing $X_{s}$ with $\mathrm{conv}({X_{s}})$
for every $s\in S$. Let $\\{x_{si}:i\in I_{s}\\}$ denote the extreme points of
$X_{s}$. Given the boundedness assumption on $X_{s}$, the MP is written as

$\displaystyle\min_{p}\ $ $\displaystyle\sum_{s\in S,i\in
I_{s}}c_{s}^{T}x_{si}p_{si}$ (2.2) s.t. $\displaystyle\sum_{s\in S,i\in
I_{s}}A_{s}x_{si}p_{si}=a,$ $\displaystyle\sum_{i\in
I_{s}}p_{si}=1\quad\forall s\in S,$ $\displaystyle p_{si}\geq 0\quad\forall
s\in S,i\in I_{s}.$

Since the MP is too large to formulate and solve explicitly, a column
generation procedure is used. The restricted master problem (RMP) is defined
by replacing $I_{s}$ in the MP with a subset $\hat{I}_{s}\subset I_{s}$ for
every $s\in S$. Suppose that the RMP is feasible ($\hat{I}_{s}$’s must have
enough columns to satisfy the first constraint in (2.2)) and let
$y,\sigma_{s}$ for $s$ in $S$ be the optimal dual solution to the RMP
corresponding to the first and second constraint respectively. To find columns
to be added to the RMP, the following subproblems, known as the pricing
problems, are solved:

$r_{s}(y)=\min_{x_{s}}\\{(c_{s}-y^{T}A_{s})x_{s}:x_{s}\in X_{s}\\}\quad s\in
S.$ (2.3)

If $r_{s}(y)\geq\sigma_{s}$ for all $s\in S$, the RMP has found the optimal
solution to the MP. Otherwise, the solutions to the pricing problems are added
to the RMP and the above process is repeated. It follows from LP duality that
given a dual value $y$

$l(y)=a^{T}y+\sum_{s\in S}r_{s}(y)$ (2.4)

is a lower bound to the MP, which also bounds (2.1).

It is known that the naive column generation approach described above suffers
from instability. To mitigate the issue, quadratic regularization on the dual
variable may be added as described by [5]. Dualizing the RMP and adding
quadratic regularization on the dual variable we get

$\displaystyle\max_{y,\sigma}\ $ $\displaystyle a^{T}y+\sum_{s\in
S}\sigma_{s}-\frac{\mu}{2}\|y-\bar{y}\|_{2}^{2}$ (2.5) s.t.
$\displaystyle\sigma_{s}\leq(c_{s}-y^{T}A_{s})x_{si}\quad\forall s\in
S,i\in\hat{I}_{s},$ $\displaystyle y,\sigma:\mathrm{free},$

where $\bar{y}$ is a regularization centre and $\mu$ is a parameter to adjust
the strength of the regularization. In the regularized column generation
procedure, the regularized RMP is used in place of the RMP. The optimal
solution to (2.5) is computed and the pricing problems are solved based on
this solution. The regularization centre $\bar{y}$ may be updated to the
current dual value $y$ when the lower bound (2.4) gets improved.

##  3 Warmstarting

In the above algorithm, the dual variable $y$ is updated iteratively to
provide tighter lower bounds. Consider the case where a near-optimal dual
value is used as an initial point. It can be shown that the lower bound $l(y)$
is continuous (in fact concave) in the dual variable $y$. Thus the algorithm
would provide a near-optimal lower bound in the first iteration. Assuming the
lower bound is sufficiently tight, the algorithm would terminate as soon as a
primal feasible solution with sufficiently small suboptimality was found by
primal heuristics. It is expected that feeding a near-optimal dual value as an
initial dual value helps the algorithm to terminate in a shorter execution
time.

It is worth noting that regularization on the dual variable is crucial in this
context. Without regularization, even if a near-optimal dual value is used as
an initial point, the algorithm is likely to be quite unstable, yielding dual
values with poor lower bounds in the following iterations. See [5] for a
further discussion.

###  3.1 Approximation

One option to find a good dual value is to solve some approximation of (2.1).
For example, [11] used the LPR, which is a relaxation of (2.1) where the
binary variables are replaced with continuous variables. First the LPR was
solved and the corresponding optimal dual value $y$ was obtained. Then, this
was used as an initial dual value to warmstart the column generation
procedure. [12] constructed an approximation by dropping some constraints,
such as integrality constraints, minimum up/down time constraints and minimum
power output constraints. The optimal dual value to this problem was used as
the initial point for a subgradient method approach to Lagrangian relaxation.
These approaches do not require any offline computation or pre-trained models.
On the other hand solving the approximations requires a certain amount of
time.

###  3.2 Initialization methods

An alternative approach, which we explore below, is to train a model which
takes $\omega$ as input and outputs a dual value $y$ that provides a tight
lower bound. The training approach is based on the work by [10] where dual
decomposition was applied to solve parametrized stochastic two-stage
programming.

Consider a neural network model $f(\omega,\theta)$ which maps the pair of
problem parameter $\omega$ and model parameter $\theta$ to $y$. We shall learn
a value of the parameter $\theta$ so that given $\omega$ the model outputs a
dual value $y=f(\omega,\theta)$ that provides a tight lower bound
$\eqref{eq_lower_bound}$. In this section, we explicitly state the dependency
of the lower bound on $\omega$ as $l(\omega,y)$, which was suppressed in
(2.4). Assume that the distribution of $\omega$ is given (e.g. the empirical
distribution based on historical data). Our goal is to maximize the expected
lower bound

$\displaystyle p(\theta)=\mathrm{E}_{\omega}[l(\omega,f(\omega,\theta))].$

From (2.4) it follows that for each $\omega$

$\displaystyle l(\omega,y)=a^{T}y+\sum_{s\in
S}r_{s}(y)=\frac{1}{|S|}\sum_{s\in
S}(a^{T}y+|S|r_{s}(y))=\mathrm{E}_{t}[\widetilde{l}(\omega,y,t)],$

where

$\displaystyle\widetilde{l}(\omega,y,t)=a^{T}y+|S|r_{t}(y)$

and subproblem $t$ is sampled uniformly from $S$. Therefore, we get

$\displaystyle
p(\theta)=\mathrm{E}_{\omega,t}[\widetilde{l}(\omega,f(\omega,\theta),t)].$

We follow the standard approach used when a neural network is trained with the
stochastic gradient descent method. That is, we sample $\omega$ and $t$,
compute the gradient of $\widetilde{l}$ with respect to $\theta$ and make a
single gradient ascent step. We then resample $\omega$ and $t$ and repeat the
process again.

The gradient of $\widetilde{l}$ with respect to $\theta$ can be computed as
follows. Fix problem parameter $\omega$ and subproblem index $t$ and compute
the dual value $y=f(\omega,\theta)$ and the component
$\widetilde{l}(\omega,y,t)$ of the lower bound corresponding to subproblem
$t$. Suppose that the model output $f(\omega,\theta)$ is differentiable with
respect to $\theta$ and that the component $\widetilde{l}(\omega,y,t)$ of the
lower bound and the optimal value $r_{t}(y)$ of the pricing problem (2.3) are
differentiable with respect to $y$. Then the gradient of
$\widetilde{l}(\omega,y,t)$ with respect to $y$ is given by

$\frac{\partial\widetilde{l}}{\partial y}=a-|S|A_{t}x_{t}^{*}$

where $x_{t}^{*}$ is the solution to the pricing problem $t$ (2.3). Using the
chain rule, we obtain

$\frac{\partial\widetilde{l}}{\partial\theta}=J\frac{\partial\widetilde{l}}{\partial
y},$ (3.1)

where $J$ is the Jacobian matrix of the neural network model output
$y=f(\omega,\theta)$ with respect to $\theta$, which is typically given by
automatic differentiation.

It is worth noting that evaluating the above quantity (3.1) is much more
faster than solving the original problem or the MP to optimality. This is
because it requires the solution of only a single pricing problem, which is
typically substantially smaller than the original problem.

Once a model is trained, it can be used to compute an initial dual value in
the regularized column generation procedure. It is expected that unlike
solving an approximation such as the LPR this model will run quickly.
Nevertheless it produces a near-optimal initial dual value with which to
warmstart the column generation procedure.

##  4 Numerical experiments

In this section, the proposed method is compared with standard approaches. All
methods are implemented with Python. IBM ILOG
CPLEX111https://www.ibm.com/products/ilog-cplex-optimization-studio is used as
the optimization solver (to solve the pricing problems and the RMP) and
PyTorch to implement the neural networks.

###  4.1 Problem

In the experiments, we consider a setup in which UC problems are solved
repeatedly with a fixed set of generators but different demand forecasts. To
assess the scalability, we consider 3 different problem sizes; problems with
200, 600 and 1,000 generators. In all cases, the length of the planning
horizon is 48 hours with a time resolution of 1 hour. The generator data is
based on [4]. Since their sets of generators contain 200 generators at most,
we combine multiple sets to create larger ones. Each generator is distinct and
combining the sets of generators does not introduce symmetry. The demand data
is based on the historical demand data in the UK published by National Grid
ESO.222https://www.nationalgrideso.com/ For each set of generators, the demand
is scaled so that the median of the daily peaks becomes 50% of the total
generation capacity. A detailed description of the problem formulation is
found in the appendix.

Dantzig-Wolfe decomposition with the regularized column generation procedure
is used to solve the UC problems. Various initialization methods are compared,
and these are described in the next section. In each iteration of the column
generation procedure, after the pricing problems are solved, a primal
heuristic is run to find primal feasible solutions. Implementation details
including the primal heuristics are given in the appendix.

###  4.2 Initialization methods

We consider the following 5 methods to initialize the dual variable.

#### Neural Network

An initialization method based on a neural network is implemented as described
in the previous section. The model consists of 4 hidden layers of 1000 units
per layer, with skip connections between hidden layers. For each set of
generators, a single neural network is trained for 24 hours with 8 CPU cores.
The parameters of the neural network are kept in shared memory and updated
asynchronously in each worker. This follows the procedure used by [10], except
that only a single neural network is trained in each setup to reduce the
computational load. The hyperparameters on the model architecture and the
training are found in the appendix.

#### Random forest

As a baseline of a trainable model, a method based on a random forest ([5]) is
also evaluated. On the same training budget (24 hours on 8 CPU cores), as many
training instances as possible are solved to 0.25% optimality. Dantzig-Wolfe
decomposition with the LPR is used as a solver. The parameter values and the
corresponding optimal dual values of the training instances are saved. After
the timelimit, a random forest is trained by supervised learning to predict
the optimal dual value given a parameter value.

#### Nearest Neighbour

Nearest neighbour is tested to provide an alternative trainable baseline. We
use the dataset created to train the random forest model. At runtime, given a
test instance, its parameter is compared against those of the training
instances. The training instance with the closest parameter value in terms of
the Euclidean distance is selected, and the corresponding optimal dual value
is used as an initial dual value.

#### LPR

This is the method based on the LPR described in Section 3. Given a test
instance, we first solve the LPR by CPLEX. Then, the optimal dual value is
used as an initial value of the column generation procedure. This method does
not require any training. This is the method used by [11].

#### Coldstart

As a naive base line, column generation is run from $y=0$. This method does
not require any training.

###  4.3 Evaluation

To test the performances of the above methods, 40 test instances of each size
are created and solved by each method. The demand data to construct the test
instances are sampled from a different year than those of the training
instances. Each method is run sequentially on a single CPU core until a
solution within 1%, 0.5% and 0.25% suboptimality is found or the time limit of
5 minutes is reached.

Table 1 shows the qualities of the lower and upper bounds obtained after one
iteration of the column generation procedure and the required time (all times
in this section are wall clock times) to run the initialization methods. Out
of 40 test instances, 8 instances are sampled and the lower bounds of the
costs are computed by running CPLEX for 4 hours per instance. Then the lower
and upper bounds obtained from the column generation procedure are scaled so
that the lower bounds from CPLEX becomes 1.0 (the upper bound obtained from
CPLEX was very close to the lower bound from CPLEX, typically the gap was
0.04% in 1000-generator cases). The table reports the average of the scaled
bounds on the 8 instances.

Clearly, coldstart outputs dual values with poor lower bounds. Although the
lower bounds of the remaining 4 methods are comparable, the lower bounds
obtained from the method based on a neural network seems to be tighter than
the lower bounds from the others. In all cases, the upper bounds are looser
than the lower bounds. Also, methods with tighter lower bounds tend to have
tighter upper bounds. However, the upper bounds are likely to vary if we use a
different primal heuristic in the column generation procedure. The solution
time of the LPR grows rapidly and it takes significant time especially on
larger problems. On the other hand, the computational time of the other
methods remains small.

Table 1: Average gap and computational time (seconds). | size: 200 | size: 600 | size: 1000
---|---|---|---
method | lb | ub | time | lb | ub | time | lb | ub | time
network | 0.99945 | 1.00497 | 0.0 | 0.99948 | 1.00394 | 0.0 | 0.99959 | 1.00336 | 0.0
forest | 0.99911 | 1.00605 | 0.0 | 0.99911 | 1.00519 | 0.0 | 0.99922 | 1.00521 | 0.0
nearest | 0.99897 | 1.00646 | 0.0 | 0.99915 | 1.00507 | 0.0 | 0.99905 | 1.00512 | 0.0
LPR | 0.99874 | 1.00512 | 4.6 | 0.99904 | 1.00387 | 17.8 | 0.99918 | 1.00316 | 29.6
coldstart | 0.00153 | 1.02122 | 0.0 | 0.00105 | 1.01981 | 0.0 | 0.00141 | 1.01961 | 0.0

Table 2 shows the number of problems solved by the column generation procedure
within the time limit of 5 minutes, the average computational time and the
average number of iterations to close the optimality gap or to reach the time
limit. The computational time reported in this table includes the time to run
the initialization methods such as solving LPR (this does not contain time
spent to train the models). For the instances that are not solved within the
time limit, no penalty is added on the computation of the average and the time
is set to be the 5 minutes, and the number of iterations reached by the 5
minute timelimit is used. For comparison, we also solve the extensive form
(2.1) to the same tolerances without decomposition (i.e. by branch and bound)
using CPLEX.

Table 2: Number of problems solved, average computational time (seconds) and iterations. |  | 1% optimality | 0.5% optimality | 0.25% optimality
---|---|---|---|---
size | method | solved | time | iters | solved | time | iters | solved | time | iters
200 | network | 40/40 | 4.6 | 1.0 | 40/40 | 8.9 | 1.8 | 40/40 | 34.0 | 6.6
| forest | 40/40 | 4.5 | 1.0 | 40/40 | 11.0 | 2.4 | 40/40 | 39.0 | 7.9
| nearest | 40/40 | 6.5 | 1.3 | 40/40 | 20.0 | 3.8 | 40/40 | 52.2 | 9.4
| LPR | 40/40 | 7.9 | 1.0 | 40/40 | 16.0 | 3.0 | 40/40 | 45.4 | 8.8
| coldstart | 37/40 | 105.3 | 12.1 | 34/40 | 173.4 | 14.0 | 27/40 | 229.6 | 16.8
| cplex | 36/40 | 195.9 | 0.0 | 35/40 | 210.0 | 0.0 | 24/40 | 264.8 | 0.0
600 | network | 40/40 | 13.4 | 1.0 | 40/40 | 20.5 | 1.5 | 40/40 | 47.7 | 3.4
| forest | 40/40 | 14.4 | 1.2 | 40/40 | 30.1 | 2.3 | 40/40 | 63.2 | 4.8
| nearest | 40/40 | 17.6 | 1.2 | 40/40 | 40.5 | 2.8 | 40/40 | 81.9 | 5.7
| LPR | 40/40 | 28.6 | 1.0 | 40/40 | 43.3 | 2.1 | 40/40 | 87.7 | 5.4
| coldstart | 38/40 | 221.5 | 11.4 | 28/40 | 258.5 | 12.6 | 13/40 | 290.2 | 13.8
| cplex | 0/40 | 300.0 | 0.0 | 0/40 | 300.0 | 0.0 | 0/40 | 300.0 | 0.0
1000 | network | 40/40 | 21.9 | 1.0 | 40/40 | 34.4 | 1.5 | 40/40 | 72.0 | 3.1
| forest | 40/40 | 21.6 | 1.1 | 40/40 | 54.9 | 2.6 | 40/40 | 106.9 | 4.9
| nearest | 40/40 | 24.2 | 1.2 | 40/40 | 66.3 | 3.0 | 40/40 | 124.5 | 5.7
| LPR | 40/40 | 45.9 | 1.0 | 40/40 | 67.4 | 1.8 | 40/40 | 119.4 | 4.2
| coldstart | 23/40 | 280.3 | 10.9 | 16/40 | 295.5 | 11.5 | 6/40 | 298.7 | 11.9
| cplex | 0/40 | 300.0 | 0.0 | 0/40 | 300.0 | 0.0 | 0/40 | 300.0 | 0.0

These results show that with appropriate regularization and a fast heuristic
to generate suitable dual values, a column generation procedure can be
warmstarted to find solutions in fewer iterations and in shorter time. Solving
the problems without decomposition is the slowest. Among the remaining 5
methods, which are based on the column generation procedure, solving the
problem with coldstart is clearly the slowest in every case.

The method to generate initial dual values must balance its quality of the
dual values and the overall computational time. One extreme example is the
LPR, which produces dual values with a tight lower bound. However the solution
time grows rapidly as the problem size increases. The method based on a neural
network needs off-line training, but it runs very quickly when solving a new
instance, resulting in significant speed up of the overall algorithm.

A typical issue of learning-based heuristics is scalability. For example,
creating a dataset of optimal solutions can be prohibitive when large-scale
problems are to be solved. In such cases, methods based on supervised learning
that uses optimal solutions become impractical. The method used in our neural
network approach, on the other hand, can be trained efficiently based on the
decomposition, which allows the method to scale. This may account for the
performance difference between the neural network approach and the other
learning-based approaches, which use supervised learning.

These observations indicate that the warmstarting method with a neural network
tends to be competitive or outperform the other methods in terms of the
computational time in almost many cases. This trend is clearer with tighter
tolerances.

Table 3 shows a breakdown of the average computational time of the column
generation procedures based on the neural network and on the LPR for the
1000-generator instances for the 40 test cases. Initialization shows the time
required to run the initialization methods. In the LPR-based method, this is
the time to solve the LPR, while the method based on a neural network this is
the time to evaluate the neural network model (this does not include the time
to train the neural network). In cases of $1\%$ tolerance, the time required
to solve the LPR is longer than the sum of the time spent on the other
routines. However as the tolerance becomes tighter, the number of iterations
and the time spent in the other routines increases, and by $0.25\%$ tolerance
the LPR initialization time is only 25% of the total time.

Table 3: Breakdown of the average computational time (seconds). tolerance | method | initialization | RMP | subproblem | primal heuristic
---|---|---|---|---|---
1% | network | 0.0 | 0.0 | 17.2 | 4.7
| LPR | 29.3 | 0.0 | 11.8 | 4.8
0.5% | network | 0.0 | 0.1 | 25.2 | 9.1
| LPR | 29.3 | 0.1 | 25.4 | 12.6
0.25% | network | 0.0 | 0.3 | 52.6 | 19.0
| LPR | 29.3 | 0.7 | 62.8 | 26.7

##  5 Conclusion

We considered methods to accelerate Dantzig-Wolfe decomposition with a column
generation procedure to solve a parametrized UC problem. We trained a neural
network to generate values for the dual variables which yield tight lower
bounds. The training was done efficiently using dual decomposition. The model
was then used to warmstart the column generation procedure. When solving a new
instance, the neural network provided a dual value quickly, which was used as
the initial dual value of the column generation procedure. The initial dual
value gave a very tight lower bound and it successfully reduced the overall
computational time to find a solution of prescribed suboptimality. The
numerical experiments show that this approach scales well and can be effective
for solving large-scale unit commitment problems.

A potential extension of this work is a UC problem with a stochastic demand
forecast. In many solution methods, UC problems with deterministic demand
forecasts appear as subproblems and are solved repeatedly. See [12] and [11]
for instance. Typically the subproblems share the same problem structure and
only differ in the cost coefficients and the demand forecasts. In such cases,
it would be of interest to use data gathered from previous instances and
warmstart the algorithm when solving a new instance.

##  Appendix A Problem formulation

We follow one of the standard formulations in literature and formulate the
following constraints:

  * •

Load balance: Generators have to meet the all demand in each time period.

  * •

Reserve: To deal with contingencies, it is required to keep a sufficient
amount of back up in each time period, which can be activated quickly.

  * •

Power output bounds: Each generator’s power output has to be within its limit.

  * •

Ramp rate: Generators can only change their outputs within the ramp rates.

  * •

Minimum up/downtime: If switched on (off), each generator has to stay on (off)
for a given minimum period. This is to avoid thermal stress in the generators
which may cause wear and tear of the turbines.

To formulate the model, we use the following notation.

  * •

Sets

    * –

$G=\\{1,2,\ldots,n_{G}\\}$: set of generators

    * –

$T=\\{1,2,\ldots,n_{T}\\}$: set of time indices where decisions are taken

  * •

Parameters

    * –

$C^{\mathrm{nl}}_{g}$: no-load cost of generator $g$

    * –

$C^{\mathrm{mr}}_{g}$: marginal cost of generator $g$

    * –

$C^{\mathrm{up}}_{g}$: startup cost of generator $g$

    * –

$P^{\max/\min}_{g}$: maximum/minimum generation limit of generator $g$

    * –

$P^{\mathrm{ru}/\mathrm{rd}}_{g}$: operating ramp up/down limits of generator
$g$

    * –

$P^{\mathrm{su}/\mathrm{sd}}_{g}$: startup/shutdown ramp limits of generator
$g$

    * –

$T^{\mathrm{u}/\mathrm{d}}_{g}$: minimum uptime/downtime of generator $g$

    * –

$P^{\mathrm{d}}_{t}$: power demand at time $t$

    * –

$P^{\mathrm{r}}_{t}$: reserve requirement at time $t$

  * •

Variables

    * –

$\alpha_{gt}\in\\{0,1\\}$: 1 if generator $g$ is on in period $t$, and 0
otherwise

    * –

$\gamma_{gt}\in\\{0,1\\}$: 1 if generator $g$ starts up in period $t$, and 0
otherwise

    * –

$\eta_{gt}\in\\{0,1\\}$: 1 if generator $g$ shut down in period $t$, and 0
otherwise

    * –

$p_{gt}\geq 0$: power output of generator $g$ in period $t$

The objective is the total cost

$\min\sum_{t\in T}\sum_{g\in
G}\left(C^{\mathrm{nl}}_{g}\alpha_{gt}+C^{\mathrm{mr}}_{g}p_{gt}+C^{\mathrm{up}}_{g}\gamma_{gt}\right).$

This is to be minimized subject to the following constraints.

  * •

Load balance

$\sum_{g\in G}p_{gt}\geq P^{\mathrm{d}}_{t}\qquad t\in T$

  * •

Reserve

$\sum_{g\in G}(P^{\max}_{g}\alpha_{gt}-p_{gt})\geq P^{\mathrm{r}}_{t}\qquad
t\in T$

  * •

Power output bounds

$P^{\min}_{g}\alpha_{gt}\leq p_{gt}\leq P^{\max}_{g}\alpha_{gt}\qquad g\in
G,t\in T$

  * •

Ramp rate

$p_{gt}-p_{g\,t-1}\leq
P^{\mathrm{ru}}_{g}\alpha_{g\,t-1}+P^{\mathrm{su}}_{g}\gamma_{gt}\qquad g\in
G,t\in T\backslash\\{1\\}$ $p_{g\,t-1}-p_{gt}\leq
P^{\mathrm{rd}}_{g}\alpha_{gt}+P^{\mathrm{sd}}_{g}\eta_{gt}\qquad g\in G,t\in
T\backslash\\{1\\}$

  * •

Minimum up/downtime

$\sum_{i=\max\\{t-T^{\mathrm{u}}_{g}+1,1\\}}^{t}\gamma_{gi}\leq\alpha_{gt}\qquad
g\in G,t\in T$ $\sum_{i=\max\\{t-T^{\mathrm{u}}_{g}+1,1\\}}^{t}\eta_{gi}\leq
1-\alpha_{gt}\qquad g\in G,t\in T$

  * •

Polyhedral/Switching constraints (to enforce binaries to work as we expect)

$\alpha_{gt}-\alpha_{g\,t-1}=\gamma_{gt}-\eta_{gt}\qquad g\in G,t\in T$
$1\geq\gamma_{gt}+\eta_{gt}\qquad g\in G,t\in T$

##  Appendix B Implementation details

In the numerical experiment, a regularized column generation procedure is
used. Initially, the regularization centre is set to the dual value given by
the initialization method. In each iteration, a lower bound is evaluated using
the current dual value and the regularization centre is updated to the current
dual value if the lower bound gets improved. Furthermore, if the lower bound
gets improved the regularization parameter is divided by two, and otherwise
multiplied by two. This closely follows the implementation of a column
generation procedure described by [11].

In every iteration, after the pricing problems are solved, a primal heuristic
based on the local search is run. Given the solutions to the pricing problems,
this primal heuristic checks the feasibility and switches the cheapest
available generators on to make the solution feasible if necessary. Note that
the solutions to the pricing problems are feasible generator schedules and
infeasibility only arises from an insufficient generation capacity to meet the
demand or insufficient spinning reserve. This is similar to the methods used
by [15] and [9] in the sense that a feasible solution is constructed from
solutions to the pricing problems by committing additional generators. In our
experiments, we observe that this primal heuristic typically finds near-
optimal primal solutions if this is used with a warmstarted column generation
procedure. If the column generation procedure using the local search primal
heuristic fails to find a primal feasible solution within a given tolerance
within 30 iterations, we use in addition a column combination based primal
heuristic. The column combination heuristic is a popular primal heuristic in
the column generation procedure in general. The idea of this heuristics is to
solve (2.1) with restricted patterns of solutions, referred to as restricted
master IP by [13]. In the $k$th iteration, we replace $X_{s}$ in (2.1) with
$\\{x_{s}^{(k-l)}\\}_{l=0,1,2}$ where $x_{s}^{(i)}$ is a solution to the $s$th
pricing problem found in the $i$th iteration. The resulting problem is still
mixed-integer programming but the solution space is much smaller than the
original problem. In our numerical experiments, we observed that the problem
is typically feasible and there is no need to add artificial variables
(columns).

##  Appendix C Hyperparameters on neural networks

The model consists of 4 hidden layers of 1000 units per layer, with skip
connections between hidden layers. Tanh is used as an activation function.
This is because in earlier experiments we observed that Tanh performed better
than ReLU, which is a more standard activation function. As [6] proposed, the
weights on the residual connections are initialized to zero. The other weights
are initialized based on the methods by [7]. All biases are initialized to
zero.

Adam ([8]) is used to learn the parameters of the neural network. The learning
rate is initialized to $10^{-4}$. The model performance is evaluated every 5
minutes, and if it fails to improve for successive 15 minutes, the learning
rate is divided by 1.5.

## References

  * [1] W. Ackooij et al.  “Large-scale unit commitment under uncertainty: an updated literature survey”  In _Annals of Operations Research_ 271.1 New York: Springer US, 2018, pp. 11–85
  * [2] J.. Bard  “Short-Term Scheduling Of Thermal-Electric Generators Using”  In _Operations Research_ 36.5, 1988, pp. 756
  * [3] D. Bertsekas, G. Lauer, N. Sandell and T. Posbergh  “Optimal short-term scheduling of large-scale power systems”  In _IEEE Transactions on Automatic Control_ 28.1 IEEE, 1983, pp. 1–11
  * [4] A. Borghetti, A. Frangioni, F. Lacalandra and C.. Nucci  “Lagrangian Heuristics Based on Disaggregated Bundle Methods for Hydrothermal Unit Commitment”  In _IEEE Power Engineering Review_ 22.12 IEEE, 2002, pp. 60–60
  * [5] O Briant et al.  “Comparison of bundle and classical column generation”  In _Mathematical programming_ 113.2 Berlin/Heidelberg: Springer-Verlag, 2008, pp. 299–344
  * [6] S. De and S. Smith  “Batch Normalization Biases Residual Blocks Towards the Identity Function in Deep Networks”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 19964–19975
  * [7] X. Glorot and Y. Bengio  “Understanding the difficulty of training deep feedforward neural networks”  In _Journal of Machine Learning Research - Proceedings Track_ 9, 2010, pp. 249–256
  * [8] D.. Kingma and J. Ba  “Adam: A Method for Stochastic Optimization” ULR https://arxiv.org/abs/1606.01885 In _The International Conference on Learning Representations_ , 2015
  * [9] A. Merlin and P. Sandrin  “A New Method for Unit Commitment at Electricite De France”  In _IEEE Transactions on Power Apparatus and Systems_ PAS-102.5, 1983, pp. 1218–1225  DOI: 10.1109/TPAS.1983.318063
  * [10] V. Nair, K. Dvijotham, I. Dunning and O. Vinyals  “Learning Fast Optimizers for Contextual Stochastic Integer Programs”  In _Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence, UAI 2018, Monterey, California, USA, August 6-10, 2018_ , 2018, pp. 591–600
  * [11] T. Schulze, A. Grothey and K. McKinnon  “A Stabilised Scenario Decomposition Algorithm Applied to Stochastic Unit Commitment Problems”  In _European Journal of Operational Research_ 261.1 Elsevier B.V., 2017, pp. 247–259  DOI: 10.1016/j.ejor.2017.02.005
  * [12] S. Takriti, J.. Birge and E. Long  “A stochastic model for the unit commitment problem”  In _IEEE Transactions on Power Systems_ 11.3 IEEE, 1996, pp. 1497–1508
  * [13] F. Vanderbeck  “Implementing mixed integer column generation”  In _Column Generation_ Springer US, 2005, pp. 331–358
  * [14] F. Vanderbeck and M… Savelsbergh  “A generic view of Dantzig–Wolfe decomposition in mixed integer programming”  In _Operations Research Letters_ 34.3 Amsterdam: Elsevier B.V, 2006, pp. 296–306
  * [15] F. Zhuang and F.. Galiana  “Towards a more rigorous and practical unit commitment by Lagrangian relaxation”  In _IEEE Transactions on Power Systems_ 3.2, 1988, pp. 763–773  DOI: 10.1109/59.192933