# Time-Optimal Quantum Driving by Variational Circuit Learning

Tangyou Huang  Department of Physical Chemistry, University of the Basque
Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  International Center of
Quantum Artificial Intelligence for Science and Technology (QuArtist)
and Department of Physics, Shanghai University, 200444 Shanghai, China
Yongcheng Ding  Department of Physical Chemistry, University of the Basque
Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  Léonce Dupays  Department
of Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg,
Luxembourg  Yue Ban  TECNALIA, Basque Research and Technology Alliance (BRTA),
48160 Derio, Spain  Man-Hong Yung  Central Research Institute, 2012 Labs,
Huawei Technologies, Shenzhen 518129, China  Department of Physics, Southern
University of Science and Technology, Shenzhen 518055, China  Shenzhen
Institute for Quantum Science and Engineering, Southern University of Science
and Technology, Shenzhen 518055, China  Adolfo del Campo  Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg,
Luxembourg  Donostia International Physics Center, E-20018 San Sebastián,
Spain  Xi Chen  Department of Physical Chemistry, University of the Basque
Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center,
University of the Basque Country UPV/EHU, Barrio Sarriena, s/n, 48940 Leioa,
Spain

###### Abstract

The simulation of quantum dynamics on a digital quantum computer with
parameterized circuits has widespread applications in fundamental and applied
physics and chemistry. In this context, using the hybrid quantum-classical
algorithm, combining classical optimizers and quantum computers, is a
competitive strategy for solving specific problems. We put forward its use for
optimal quantum control. We simulate the wave-packet expansion of a trapped
quantum particle on a quantum device with a finite number of qubits. We then
use circuit learning based on gradient descent to work out the intrinsic
connection between the control phase transition and the quantum speed limit
imposed by unitary dynamics. We further discuss the robustness of our method
against errors and demonstrate the absence of barren plateaus in the circuit.
The combination of digital quantum simulation and hybrid circuit learning
opens up new prospects for quantum optimal control.

##  I Introduction

Following the vision by Feynman [1], quantum simulation has acquired a
potentially-disruptive role in the development of contemporary science and
technology, given the prospects of harnessing the advantage of using a quantum
computer (QC) for specific applications. In recent decades, quantum simulation
has been used to probe the dynamics of condensed matter systems [2, 3], for
quantum chemistry [4, 5], and as a test-bed for nonequilibrium statistical
mechanics, e.g., in studying thermalization and nonequilibrium behavior of
many-body systems [6, 7]. Quantum simulation is also expected to impact high-
energy physics, given the potential to facilitate the study of lattice gauge
theories [8] and gauge-gravity duality [9, 10], among other examples.

The use of a digital quantum simulator (DQS) based on the gate model offers a
prominent approach in current Noisy Intermediate-Scale Quantum (NISQ) devices
[11] and has gained relevance with theoretical and experimental progress [12,
13, 14]. In particular, DQS can be used to implement variational quantum
algorithms (VQAs), under development for quantum optimization [15], quantum
machine learning [16], and quantum control [17]. Their formulation generally
approximates the continuous time evolution by discrete, finite Trotter steps
[18, 19, 20, 21] implemented by a sequence of quantum gates, with controlled
accuracy, in principle. However, balancing the number of Trotter steps and
imperfections of quantum circuits in experiments is still a fundamental
challenge. In this sense, various optimization scenarios aim at quantum error
mitigation [22, 23, 24] for achieving a good precision of the quantum
simulation with limited quantum resources. Among those, the machine-learning-
enhanced optimization protocol [25, 26, 27, 28] utilizes a feedback loop
between the quantum device and a classical optimizer. This approach is
particularly useful in the field of hybrid quantum algorithms [29, 30] with
current quantum hardware. Nonetheless, solving the quantum optimal control
problem by VQAs in a NISQ device is still an open challenge [17].

In this work, we propose a circuit learning scheme based on gradient-descent
(GD) for time-optimal quantum control. As a concrete example, we consider a
quantum particle trapped in time-varying parabolic potential. We use a qubit
register and encode the spatial wave function using the basis of $n$-qubit
states. We then reproduce the exact state evolution on a designed quantum
circuit using a digital algorithm [31, 32]. We optimize the control function
to achieve maximum-fidelity control by using the GD-based circuit learning. We
further unveil the connection between a control phase transition and the
quantum speed limit, i.e., the minimum time for a quantum state to evolve into
a distinguishable state under a given dynamics. We demonstrate that the
fidelity-based GD method avoids a large number of measurements by comparison
to the reinforcement learning protocol [33] and show how it can be accelerated
by choosing different quantum quantities as the cost function.

In the following two sections, we introduce the quantum algorithm for the
circuit realization of a quadratic Hamiltonian and discuss the time-dependent
harmonic oscillator as an example. We then explore a fidelity-based GD method
for maximum-fidelity control in a nonequilibrium expansion process and
characterize the efficiency through various cost functions. The relation of
the quantum speed limit to the control phase transition is then discussed.
Finally, we establish the fault tolerance of our method against the quantum
errors in experiments and also address the problem of barren plateaus in the
parameterized circuit.

##  II PRELIMINARIES AND NOTATION

###  II.1 Time-dependent quantum harmonic oscillator

We exemplify our approach by considering the time-dependent harmonic
oscillator (TDHO), described by the Hamiltonian

$H(t)=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}(t)[\hat{x}-x_{0}(t)]^{2},$
(1)

where $\omega(t)$ and $x_{0}(t)$ are tunable and represent the trap frequency
and the location of the trap center, respectively. The TDHO is an ideal model
for benchmarking quantum control algorithms since its dynamics admits exact
closed-form solutions. In particular, we focus on the case with $x_{0}(t)=0$
and look for the expansion of the wave packet induced by a modulation of the
trap frequency $\omega(t)$ from an initial value $\omega_{0}$ to a final one
$\omega_{f}$. This model has many applications, including the cooling of a
particle in an optical trap [34], mechanical resonators [35], and tunable
transmon superconducting qubits [36]. The ground state of $H(0)$ evolves into
the time-dependent Gaussian state [34]

$\displaystyle\Psi(t,x)=\left(\frac{m\omega_{0}}{\pi\hbar
b^{2}}\right)^{1/4}\exp\left[-\frac{i}{2}\int_{0}^{t}\frac{\hbar\omega_{0}}{b^{2}}dt^{\prime}\right]$
$\displaystyle\times\exp\left[\frac{im}{2\hbar}\left(\frac{\dot{b}}{b}+i\frac{\omega_{0}}{b^{2}}\right)x^{2}\right],$
(2)

where the time-dependent scaling factor $b(t)>0$ characterizes the width of
the wave packet and satisfies the auxiliary equation

$\displaystyle\ddot{b}+\omega^{2}(t)b=\frac{\omega_{0}^{2}}{b^{3}}.$ (3)

A primary numerical solver of quantum dynamics is the so-called split-operator
method (SOM), also known as the split time propagation scheme [37]. For the
sake of convenience, one usually sets dimensionless variables based on
physical units of energy $\epsilon=\hbar\omega_{0}$, length
$b_{\text{HO}}=\sqrt{\hbar/m\omega_{0}}$, and time $\tau=1/\omega_{0}$. In a
classical computer, one defines a $N$-dimensional vector as $\Psi(\textbf{r})$
for encoding the amplitude of the wave function on the space grid
$\textbf{r}=[x_{0},~{}x_{1},\cdots,x_{N-1}]$. Note that the kinetic energy
operator $\hat{T}=\hat{p}^{2}/2$ and the potential operator
$\hat{V}=\omega^{2}(t)\hat{x}^{2}/2$ do not commute. Thus, the following
approximation stands for small $dt$ with an error $\mathcal{O}(dt^{3})$

$e^{-iHdt}\approx
e^{-\frac{i}{2}\hat{V}dt}e^{-i\hat{T}dt}e^{-\frac{i}{2}\hat{V}dt},$ (4)

evolving the wave function in a Trotter step. A common trick in implementing
this method uses forward and inverse Fourier transforms to change the
representation of the quantum state between the real space r and momentum
space k, in which the kinetic energy operator becomes diagonal in k,
simplifying the numerical calculation.

###  II.2 Time-optimal control

The frictionless expansion of quantum particles trapped in a time-varying
harmonic trap can be formulated as a time-optimal control problem by
minimizing the time of the process $t_{f}$ [34, 38, 39, 40, 41]. It follows
from Pontryagin’s maximum principle that the control Hamiltonian for all
$t\in[0,t_{f}]$ takes the form [38]

$H_{c}[x_{1},x_{2},p_{1},p_{2}]=p_{1}x_{2}+\frac{p_{2}}{x_{1}^{3}}-p_{2}x_{1}u(t),$
(5)

where the state $x_{1}=b$, $x_{2}=\dot{b}/\omega_{0}$, and the controller
$u(t)=\omega^{2}(t)/\omega_{0}^{2}$ are governed by the Ermakov equation (3).
Here, $p_{1}$ and $p_{2}$ are the conjugate momentum of $x_{1}$ and $x_{2}$,
respectively. Substituting the control Hamiltonian into the canonical equation
leads to the cost equations

$\displaystyle\dot{p_{1}}$ $\displaystyle=$
$\displaystyle\left(u+\frac{3}{x_{1}^{4}}\right)p_{2},$ (6)
$\displaystyle\dot{p_{2}}$ $\displaystyle=$ $\displaystyle-p_{1}.$ (7)

If the controller is bounded as $\delta_{1}\leq u(t)\leq\delta_{2}$, the time-
optimal control has a bang-bang form, i.e., it is a piece-wise function and
constant in each interval. For a specific problem with $b(0)=1$ and
$b(t_{f})=\sqrt{\omega_{0}/\omega_{f}}=\gamma$, consider the feasible three-
jump protocol

$u(t)=\begin{cases}1~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t=0\\\
\delta_{1}~{}~{}~{}~{}~{}~{}~{}0<t\leq t_{1}\\\
\delta_{2}~{}~{}~{}~{}~{}~{}~{}~{}t_{1}<t<t_{1}+t_{2}\\\
1/\gamma^{4}~{}~{}~{}~{}~{}t=t_{f}^{\text{opt}}=t_{1}+t_{2}\end{cases},$ (8)

where the switching time $t_{1}$ and the optimal operation time
$t_{f}^{\text{opt}}=t_{1}+t_{2}$ can be calculated by integrating the Ermakov
equation (3) by using boundary conditions. This yields the closed-form exact
time-optimal driving protocol $u(t)$ with

$\displaystyle t_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\delta_{1}}}\sinh^{-1}\sqrt{\frac{\delta_{1}(\gamma^{2}-1)(\gamma^{2}\delta_{2}-1)}{(\delta_{1}-\delta_{2})\gamma^{2}(1-\delta_{1})}},$
$\displaystyle t_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\delta_{2}}}\sin^{-1}\sqrt{\frac{\delta_{2}(\gamma^{2}-1)(1-\gamma^{2}\delta_{1})}{(\delta_{1}-\delta_{2})(1-\gamma^{4}\delta_{2})}}.$
(9)

Figure 1:  The circuit realization of time-evolution operator for a quadratic
Hamiltonian as defined by Eq. (22). The single phase gate and Pauli-X gate are
schematically represented by the symbols $R(\cdot)$, $X$, and the
corresponding control gates are shown by vertical lines with black circles.

###  II.3 Quantum speed limit

Quantum speed limits (QSLs) provide fundamental upper bounds on the speed of
quantum evolution [42]. They have wide-spread applications ranging from
quantum metrology to optimal control [43, 44, 45]. QSLs are formulated by
choosing a notion of distance between quantum states and identifying a maximum
speed of evolution. For isolated systems described by a time-independent
Hamiltonian, two seminal results are known. The Mandelstam-Tamm QSL determines
the maximum speed of evolution in terms of the energy dispersion [46], while
the Margolus-Levitin bound uses the mean energy above the ground state instead
[47]. The interplay of these bounds has recently been demonstrated in a
trapped system made of ultracold atoms that are suddenly quenched [48]. For a
generic driven system, only an analog of the Mandelstam-Tamm bound is known
[49, 50, 51].

Consider the quantum unitary dynamics generated by a time-dependent
Hamiltonian according to the Schrödinger equation. The distance between the
initial state and the time-dependent state in projective Hilbert space can be
quantified by the Bures angle

$\mathcal{L}(\psi_{0},\psi_{t})=\arccos(|\langle\psi_{0}|\psi_{t}\rangle|)\in[0,\pi/2].$
(10)

The minimum time scale required to sweep a given Bures angle is lower bounded
by

$\tau_{\rm QSL}=\frac{1}{\overline{\Delta E}}\mathcal{L}(\psi_{0},\psi_{t}),$
(11)

where the speed of evolution is set by the time-averaged energy dispersion

$\overline{\Delta
E}=\frac{1}{t}\int_{0}^{t}ds\sqrt{\langle\psi_{s}|\hat{H}(s)^{2}|\psi_{s}\rangle-\langle\psi_{s}|\hat{H}(s)|\psi_{s}\rangle^{2}}.$
(12)

The QSL $\tau_{\rm QSL}$ is thus approached by maximizing the energy
dispersion at all times. In control protocols for wave-packet expansion, we
consider the evolution of the ground state of the trap with initial trapping
frequency $\omega_{0}$ to the ground state of the trap with final frequency
$\omega_{f}$. Provided that the control protocol, specified by $u(t)$, has
unit efficiency in preparing the target state, the Bures angle is fixed, and
the corresponding QSL reads

$\tau_{\rm QSL}=\frac{\hbar}{\overline{\Delta
E}}\sqrt{\frac{2\gamma}{1+\gamma^{2}}}.$ (13)

We note that $\tau_{\rm QSL}$ should be distinguished from the minimum time
$t^{\rm opt}_{f}$ in the preceding time-optimal control with the trap
frequency bounded.

###  II.4  Fidelity susceptibility

Generally, the final state $\rho_{f}=|\Psi_{f}(x)\rangle\langle\Psi_{f}(x)|$
upon completion of a control protocol at $t=t_{f}$ differs from the target
state $\rho_{\rm tar}=|\Psi_{\rm tar}(x)\rangle\langle\Psi_{\rm tar}(x)|$ one
wishes to prepare. Let us consider the fidelity $F$ between these two states

$F(\rho_{\rm tar},\rho_{f})=\left[{\rm
Tr}(\sqrt{\sqrt{\rho_{\text{tar}}}\rho_{f}\sqrt{\rho_{\text{tar}}}})\right]^{2},$
(14)

where ${\rm Tr}(\cdot)$ denotes the trace operation.

The fidelity susceptibility $\chi_{f}$ quantifies the fidelity response to a
slight change of driving parameter [52, 53, 54]. For a functional Hamiltonian
$H(f_{t})$ parameterized by $f_{t}$, let $|\Psi_{0}(f_{t})\rangle$ be the
ground state. We assume $f_{t}$ to be a function of time and consider a
variation on the control function $f_{t}\to f_{t}+\delta f$, where $\delta
f\to 0$ is small enough to apply perturbation theory. As a result, the
perturbed ground state is $|\Psi_{0}(f+\delta f)\rangle$. The fidelity
susceptibility, without loss of generality, is defined as

$\chi_{f}\equiv\frac{-2\ln(F_{\delta f})}{\delta f},$ (15)

where the fidelity $F_{\delta f}=F[\rho_{0}(f),\rho_{0}(f+\delta f)]$. The
fidelity susceptibility quantifies the sensitivity of the fidelity to
variations of the control functions. In other words, the fidelity
susceptibility can be used as a cost function to accelerate the convergence of
the optimization process. For the sake of simplicity, we assume that $\delta
f\to 0$ is a time-independent real value.

##  III Quantum circuit realization of quadratic Hamiltonian

Next, we present the algorithm for the circuit realization of quadratic
Hamiltonians using a finite set of elementary quantum gates. We focus on DQS
of the continuous-variables system and encode a wave packet onto a $n$-qubit
register. Quantum states of this register can be described in binary notation

$|\Phi\rangle=\sum^{2^{n}-1}_{i=0}c_{i}|i\rangle,$ (16)

using the computational basis
$|i\rangle=|q_{n-1}\rangle\otimes\cdot\cdot\cdot\otimes|q_{1}\rangle\otimes|q_{0}\rangle$
with $q_{0},q_{1},\dots,q_{n-1}\in\\{0,1\\}$ , and the corresponding
amplitudes $c_{i}$ normalized as $\sum^{2^{n}-1}_{i=0}|c_{i}|^{2}\equiv 1$. To
solve the time-dependent Schrödinger equation on a quantum computer with a
$n$-qubit register, we discretize the continuous variables associated with the
spatial coordinate $x$ and time $t$, and subsequently map the coordinate space
$x$ into the Hilbert space of $n$ qubits. Specifically, the compact continuous
spatial domain $x\in[-L,L]$ is approximated by a lattice of $2^{n}$ points
spaced by a constant interval $dx=L/(2^{n-1}-1)$. A wave packet can be encoded
in the state of the $n$-qubit register as

$\displaystyle|\Phi\rangle$ $\displaystyle=$
$\displaystyle\sum^{2^{n}-1}_{i=0}\Psi(x_{i})|i\rangle$ (17) $\displaystyle=$
$\displaystyle\Psi(x_{0})|0\cdots 0\rangle+\dots+\Psi(x_{2^{n}-1})|1\cdots
1\rangle,$ (18)

which reproduces the vectorized wave function for the following quantum analog
of the SOM in the Hilbert space. As in numerical discretization methods, this
encoding of $\Psi(\textbf{r})$ provides, in principle, satisfactory accuracy
when the lattice length $dx$ is much smaller than any characteristic length
scale of the wave packet. The preparation of arbitrary initial qubit states
$|\Phi\rangle$ based on the initialized wave-packet $\Psi(\textbf{r})$ can be
approached by the variational quantum eigensolver (VQE)

$|\tilde{\Phi}\rangle=\prod_{i=1}^{p}\left[\prod_{q=0}^{n-1}\left(U^{q,i}\right)U_{\text{ENT}}\right]|+\rangle^{\otimes
n},$ (19)

where $|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ is a single-qubit
state, the unitary
$U^{q,i}(\theta)=R_{z}^{q,i}(\theta_{1}^{q,i})R_{x}^{q,i}(\theta_{2}^{q,i})R_{z}^{q,i}(\theta_{3}^{q,i})$
is a universal single-qubit gate, and $U_{\text{ENT}}$ represents CNOT gates
that entangle the neighboring qubits with periodic boundary conditions. In
this way, an approximated initial state $|\tilde{\Phi}\rangle$ is prepared by
optimizing $3pn$ parameters to minimize the cost function.

Next, consider the digital quantum simulation aimed at reproducing an
equivalent SOM. The dynamics of the wave packet is described by

$\Phi(t+dt)\approx e^{-iH(t)dt}\Phi(t),$ (20)

where $e^{-iHdt}$ is the time-evolution operator for the time step $dt$.

We use the quantum Fourier transform (QFT) as the quantum analog of the
inverse discrete Fourier transform, which is key to the efficient
implementation of the SOM. Hence, the wave function $|\Phi\rangle$ is evolved
as

$|\tilde{\Phi}(t+dt)\rangle=\mathcal{V}(t)_{dt/2}\text{QFT}\mathcal{T}(t)_{dt}\text{QFT}^{\dagger}\mathcal{V}(t)_{dt/2}|\tilde{\Phi}(t)\rangle,$
(21)

where $\mathcal{V}(t)$ and $\mathcal{T}(t)$ are the potential operator and the
kinetic-energy operator in the real space and momentum space, respectively. In
other words, they are both quadratic, and their diagonal elements can be
written as

$\mathcal{A}_{jj}=\exp\left\\{-idt[h(jdx+x_{0})^{2}+\sigma]\right\\},$ (22)

where other off-diagonal elements are zero. The preliminary result in Ref.
[31] demonstrates that the quadratic Hamiltonian can be exactly decomposed
into a quantum circuit. In Fig. 1, we plot the quantum circuit for
implementing a quadratic Hamiltonian in the computational basis for DQS. The
decomposition is verified by the DQS of nonadiabatic processes in molecular
systems [32].

##  IV Fidelity-based Gradient Descent

###  IV.1 Initial state preparation

The first step is to encode the information of the wave packet into the state
of the qubit register. The accuracy of the preparation of a target state of
qubits by VQE depends on the qubit number $n$ and the parameter depth $p$, see
Eq. (19). As we know, the depth $p$ and qubit number $n$ exponentially
increase the complexity of VQE. In Fig. 2, we compare the fidelity of state
preparation in the coefficient grid ${n,p}$ in (a). In (b), without loss of
accuracy, we choose $n=6$, $p=4$, and present the resulting $n$-qubits states
for the corresponding density (II.1) with $\omega_{0}=1$. In what follows, the
numerical results are produced by the quantum simulator statevector simulator
on the qiskit platform, which admits no errors, decoherence, and imperfections
at all. We will consider the noise of an actual quantum device in the
discussion.

Figure 2:  The fidelity of states preparation by VQE as a function of qubit
number $n$ and parameter depth $p$ in (a), and corresponding probability
distribution of qubits $|\Phi|^{2}$ compared with density of wave function
$|\Psi(x)|^{2}$ in (b), with fidelity $F=0.996$ for $n=6$ and $p=4$.

###  IV.2 The maximum-fidelity control

A parametric optimization problem is usually mapped into the minimization of a
given cost function to find its local minimum with the Gradient Descent (GD)
algorithm. The optimal solution $M=\\{m_{0},m_{1},...,m_{i}\\}$ is obtained by
minimizing the cost-value $c=J(M)$, and can be expressed as

$M^{\text{opt}}=\min_{c}J(M).$ (23)

In our case, the control function is the trap frequency $f(t)=\omega^{2}(t)$
that is piece-wise on the discrete time $t\in[0,t_{f}]$ involving
$N_{t}$-intervals. Accordingly, the control tuple
$f(t)=\\{f(0),f(dt),...,f(t_{f})\\}$ is constrained by
$\delta_{1}\leq|f(t)|\leq\delta_{2}$, $|f(t+dt)-f(t)|\leq\Delta f$ and the
boundary conditions $f(0)=1$ and $f(t_{f})=0.01$ are considered for trap
expansion with $\omega_{f}/\omega_{0}=1/10$. Then, we optimize
$\omega$-dependent parameters $M=\\{\theta_{3}(\omega)\\}_{t}$ in the circuit
shown in Fig. 1 by minimizing the loss value which can be generated by the
measurements on the qubits at $t=t_{f}$. Here, we exploit three different cost
functions: the infidelity (IF) $(1-F)$, the fidelity susceptibility (FS)
$\chi_{f}$, and the Bures angle (BA)
$\mathcal{L}(\psi_{t},\psi_{\text{tar}})$. For simplicity, we initialize the
controller $f(t)$ with a linear dependence of the form
$f(t)=(\omega_{f}^{2}-\omega_{0}^{2})(t/t_{f})+\omega_{0}^{2}$. This yields a
parametric constrained minimization problem in this case. We exemplify the
optimization process of finding maximum-fidelity policy for different cost
functions in Fig. 3 by using the optimizer $\tt{SLSQP}$ [55] based on the
$\tt{scipy}$ [56] platform.

Figure 3: (a) The infidelity $1-F$ as a function of the training iteration
using three different loss functions: fidelity susceptibility (FS) with
$\delta f=0.01$, infidelity (IF), and Bures angle (BA). (b) The infidelity as
a function of the iteration step for different $\delta f$ is compared when the
fidelity susceptibility (FS) is taken as a loss function. Parameters:
$N_{t}=50$, $\omega_{f}=0.1$, $n=6$, $p=4$, $\delta_{1}=10^{-6}$,
$\delta_{2}=1$, and $t_{f}=3.152$.

To this end, we assign the total time $t_{f}=t_{f}^{\text{opt}}$ calculated in
(9), and use the GD method to minimize the cost function for obtaining the
maximum-fidelity control function. In Fig. 3, we present infidelities versus
the GD iteration when using each cost function. The control with maximum
fidelity is obtained with their convergence. One can see that the learning
rate of FS outperforms the others with the same optimizer in (a), and we
compare the FS-based GD with various coefficients: $\delta
f=[10^{-4},10^{-3},10^{-2},10^{-1},10^{0}]$ in (b). Because the greatest
convergent rate of the optimization process arises when $\delta f\leq
10^{-3}$, we employ the FS $\chi_{f}$ as the cost function with the
coefficient $\delta f=10^{-3}$ to design the maximum fidelity policy in the GD
algorithm. It is worth emphasizing that reducing the training iteration is
equivalent to decreasing the accumulation of operation errors. As a result, we
can improve the accuracy of DQS with the same quantum volume.

By considering the trade-off between the complexity and accuracy of digitized
circuits, we analyze the fidelity achieved by the maximum-fidelity policy for
different numbers of Trotter steps $N_{t}$ and the constraints on control step
$\Delta f$. Trotter steps $N_{t}$ and the constraints on control step $\Delta
f$ determine the depth of the circuit and the continuity of the control
function, respectively. In Fig. 4, we present the fidelity density as a
function of $N_{t}$ and $\Delta f$ in (a) and show the maximum-fidelity
control protocols for various $N_{t}$ and $\Delta f$, where the counter curves
corresponding to the fidelity $F=0.999,0.99$ are marked. In addition, in Fig.
4(b), the results from circuit learning are compared with the the optimal-time
protocol $t_{f}=t_{f}^{\text{opt}}$ produced by bang-bang control (8). The
higher accuracy of DQS requires the larger Trotter step $N_{t}$ and higher
computation complexity. On the other hand, the control function becomes
smoother when $\Delta f\to 0$, with increasing the Trotter step and circuit
complexity. In this context, we choose $N_{t}=50$ and $\Delta f=1$ in the
following calculations.

Figure 4: (a) Fidelity as function of $\Delta f$ and $N_{t}$ with total time
$t_{f}=3.152$. The two dashed contour curves refer to $F=0.99,0.999$ and are
labeled by $-2,-3$, respectively. In (b), we present the trained fidelity-
optimal controls compared to the bang-bang control (solid gray line). The
fidelity takes values $F=0.84,0.98,0.998,0.9998$ in the four cases illustrated
in the legend from top to bottom. Other parameters are chosen as in Fig. 3

.

###  IV.3 Control phase transition at quantum speed limit

In the previous section, we described an efficient GD-based hybrid algorithm
to find the maximum-fidelity control in a quantum device by considering the
loss function $J$, Trotter step $N_{t}$, and step length $\Delta f$. Next, we
shall prove the control phase transition (CPT) appears at the critical point
due to the QSL. The use of optimal control to reach the QSL in quantum state
manipulation has been discussed in [43]. However, the authors in [43] intended
to find maximum fidelity control by reducing infidelity, which is a highly
time-consuming task. This can be improved by introducing the fidelity
susceptibility $\chi$, as we have shown in the previous section. In the space
of protocols, the control phase transition is associated with abrupt changes
in an optimal control function, satisfying given constraints as the duration
of the process is varied [33]. In particular, the maximum-fidelity control
function is unique when $t\leq t^{\rm opt}$. It is expected that the QSL can
be reached at the point of CPT.

Figure 5:  (a) The control phase diagram, with the fidelity-optimal control
sequence $f(t)=\omega^{2}(t)$ as a function of $t/t_{f}$ for different
$t_{f}\in[2,5]$. (b) Logarithm of the infidelity $\log_{10}(1-F)$ as a
function of ${t_{f},N_{t}}$. The solid line is included as a reference in (a)
and (b) and corresponds to the optimal bang-bang control with $t_{f}^{\rm
opt}=3.152$. (c) The trained fidelity-optimal controls for different values of
the total time $t_{f}$ are compared with the optimal bang-bang control, (thick
gray line). (d) The corresponding time-averaged energy dispersion
$\overline{\Delta E}$ for maximum-fidelity control. Parameters agree with
those in Fig. 3.

In practice, we consider the expansion process: from initial state $\psi_{0}$
($\omega_{0}=1$) to target states $\psi_{\text{tar}}$ ($\omega_{f}=0.1$) with
constrain $\delta_{1}=10^{-6}$, $\delta_{2}=1$ by assuming $\omega^{2}(t)>0$,
and subsequently generate the maximum-fidelity control sequence $f(t)$ for
$t_{f}\in[2,5]$, where other parameters are same as those in Fig. 3.
Consequently, we present the control phase diagram in Fig. 5(a) and illustrate
several control functions of selected $t_{f}$ compared with the optimal bang-
bang control (8). The control function is suddenly converted to a non-bang-
bang-typed phase at the transition point $t_{f}\approx t^{\text{opt}}_{f}$,
i.e., the control phase transition point. In addition, we note that there
exists only one solution of the maximum-fidelity control function when
$t_{f}\leq t^{\text{opt}}_{f}$ for fulfilling the requirement of time-energy
bound (11). Here, we confirm the minimum time $t_{f}^{\min}$ when the maximum-
fidelity is larger than 0.999, and compare the result for different $N_{t}$ in
(b) of Fig. 5. It is evident in Fig. 5(c) that the accuracy of optimal time
produced by circuit learning essentially depends on the $N_{t}$, and we set
$N_{t}=50$ for the criteria of $\log_{10}(1-F)\sim-4$ in the following
calculations.

We emphasize that the time-optimal driving obtained here differs from the QSL
but is closely related to it. Specifically, the time-optimal driving is
bounded in terms of the trap frequency by contrast to the QSL, which is
bounded in terms of the time-averaged standard deviation of the energy, see
Eq. (13). The former is a weaker and more conservative bound, as the energy
fluctuations can be upper-bounded in terms of the frequency. In Fig. 5(d), we
display the energy dispersion of maximum-fidelity control and compare it with
the corresponding time-optimal driving. The energy dispersion becomes
consistently unique when $t_{f}\leq t^{\rm opt}_{f}$ before the point of CPT.
Moreover, the energy dispersion for the maximum fidelity control is slightly
smaller than in the time-optimal bang-bang control with bounded trap
frequency. In this sense, one can approach the QSL in time-optimal driving
when an additional energy cost is allowed by relaxing the trap frequency
bound.

##  V Discussion

A significant source of discussion is the robustness of VQAs in an environment
with stochastic perturbations. In a real quantum computer with NISQ hardware,
imperfections are unavoidably induced as a result of a finite number of
measurements and a noisy environment. As emphasized, the previous results are
produced by the quantum simulator statevector simulator in the qiskit
platform, with no errors, decoherence, and imperfections at all. In this
section, we implement our method in a noise-associated quantum device
simulated by qasm simulator with $N_{m}$ measurement shots. The performance of
the GD algorithm depends on the tolerance of the optimizer to errors. Thus, we
shall balance the GD induced by noise and the parameter variance in a training
landscape. Let us recall the definition of noise in the framework of quantum
information processing. In general, the $n$-qubits register is coupled with an
environment $\varepsilon$, leading to the nonunitary evolution of the system.
Initially, we assume the density operators of the register
$\rho(t_{0})=\rho_{0}$ and the environment to be decoupled so that the
composite state is given by the tensor product $\rho\otimes\varepsilon$. For
any global unitary operator $U$ describing the dynamics of the composite
state, the reduced evolution of the register reads

$\rho(t)={\rm
Tr}[U(\rho(t_{0})\otimes\varepsilon)U^{\dagger}]\equiv\xi(\rho_{0}).$ (24)

This superoperator $\xi(\cdot)$ can be implemented for simulating a noise
model in a quantum circuit. The noisy quantum channel describes the nonunitary
evolution of the time-varying density state in the Kraus representation

$\rho(t)=\sum_{k}E_{k}\rho(t_{0})E_{k}^{\dagger},$ (25)

where $E_{k}$ satisfy the trace-preserving condition
$\sum_{k}E_{k}E_{k}^{\dagger}=\textbf{1}$. Since we perform the measurement on
the qubits register only at $t=t_{f}$, imperfections induced by any kind of
noise result in fluctuations of the measurement accuracy. In this sense, we
shall primarily consider the bit-flip error of measurements in a real quantum
computer. Assume the system’s noise flips $|0\rangle$ and $|1\rangle$ with
probability $\beta$. The superoperator for this bit flip noise can be
expressed as

$\mathcal{\xi}_{BF}(\rho)=(1-\beta)\rho+\beta X\rho X,$ (26)

where the corresponding Kraus operators are
$\\{\sqrt{1-\beta}\mathbb{I},\sqrt{\beta}X\\}$, in terms of the identity
$\mathbb{I}$ and the Pauli operator $X$.

Figure 6:  (a) The average gradient $\overline{|\partial_{\theta_{k}}J|}$ as a
function of the qubit number $n$ and number of layers given by the polynomial
function $N_{t}=5\times n$. (b) Infidelities as a function of the training
iteration for values of the noise strength $\beta=0,0.02,0.04,0.06$ based on
the $\tt SPSA$ optimizer. Other parameters are equal to those in Fig. 3.

Next, we discuss whether the barren plateau [57] phenomenon occurs here. The
barren plateau refers to the fact that the gradient of an observable vanishes
exponentially as a function of qubits number in a training landscape of VQAs.
It has been widely studied in various ansatze of deep circuits [58]. In
general, the gradient of an objective function is calculated by mean of the
parameter-shift rule [59, 60], expressed as
$\partial_{\theta_{k}}J=\frac{1}{2}[J({\theta_{k}}+\frac{\pi}{2})-J({\theta_{k}}-\frac{\pi}{2})]$
for an arbitrary trainable parameter $\theta_{k}$ in the circuit. In this
sense, we define the average of the absolute gradient over $N_{r}$ random
initializations

$\overline{|\partial_{\theta_{k}}J|}=\sum_{i=1}^{N_{r}}\frac{1}{2N_{r}}\left|J_{i}({\theta_{k}}+\frac{\pi}{2})-J_{i}({\theta_{k}}-\frac{\pi}{2})\right|.$
(27)

Since these three objective functions we proposed in the previous section are
all involve the fidelity, we here provide numerical analysis of the average
gradients of the fidelity $F$. Essentially, the probability distribution of
qubit states obeys a statistical precision of order $1/\sqrt{N_{m}}$, and the
fidelity defined in Eq. (14) meets the same criteria. The derivative of an
observable with respect to an arbitrary trainable $\theta_{k}$ in the circuit
is a linear function of the gradient with respect to the corresponding control
parameter $f_{k}=f(k)$ at $k$-th Trotter step:
$\partial_{f_{k}}J=2c\cdot\partial_{\theta_{k}}J$, with $c$ being a real
number. Thus, the gradient of the objective function with respect to the
control parameter $f(t)$ shares the analytic expression in Eq. (27), see also
the detail in Appendix A. Also, we calculate the average gradient Eq. (27)
over $N_{r}=50$ random initialization of $f(t)$ for various qubit number $n$,
while the Trotter step is taken as a polynomial function as $N_{t}=ploy(n)$.
In Fig. 6(a), we demonstrate that the barren plateau is avoided for the
average gradient of a correlating parameter $\theta_{k}$. The main reason for
the absence of barren plateau is the reduction of ansatz’s expressibility, due
to the strong correlation of parameter $\theta_{k}$ depending on the
controller $f(t)$ in our method, a common feature with the recent work in Ref.
[61]. In this regard, we set measurement shots $N_{m}=8192$ for statistical
accuracy and energy saving by considering the gradient magnitudes, as shown in
Fig. 6(a). Moreover, we apply the optimizer of simultaneous perturbation
stochastic approximation (SPSA), which is widely used for solving an
optimization problem with statistical noise [62, 63]. In Fig. 6(b), we present
the infidelity as a function of the training iteration for
$\beta=0,0.02,0.04,0.06$, where the infidelity for noise-free case ($\beta=0$)
convergent to $\sim 10^{-2}$ which obeys the criteria of $\sim
1/\sqrt{N_{m}}$. Moreover, the performance of SPSA is compared with various
optimizers in Appendix B.

Let us discuss the circuit complexity in terms of the qubits number $n$ and
the number of Trotter steps $N_{t}$. The whole circuit consists of $N_{t}$
circuit units that simulate each unitary operation $\hat{U}(dt)=e^{-iHdt}$, as
depicted in Fig. 1. In the absence of $\theta_{1}$ and $\theta_{2}$, the gate
number of each circuit unit, including the potential operator
$\mathcal{V}(dt)$ and the kinetic-energy operator $\mathcal{T}(dt)$ in real
space and momentum space, is $N_{unit}\sim 2n^{2}$. In addition, the operation
of QFT and iQFT requires $\sim n^{2}/2$ control-phase gates. Consequently, the
total number of gates for our ansatz is proportional to a quadratic function
of $n$, yielding $5N_{t}n^{2}/2$. To find the minimal-time control on the QC
with reasonable precision, one can increase the qubit number $n$ and Trotter
step $N_{t}$, with exponentially enlarged Hilbert space. But this leads to the
quadratic size increase in circuit complexity. Recently, alternative methods
inspired by the Grover-Rudolph algorithm [64] are worked out for the quantum
state preparation, which are expected to reduce its complexity in this
direction.

Finally, we discuss the control problem beyond the quadratic Hamiltonian on
which we have focused. For our case study, one can introduce perturbations of
the trap, e.g., a time-independent anharmonicity involving an operator
$x^{4}$, which is no longer quadratic. Although its exact decomposition into
quantum circuits does not exist, one can still approximate the evolution block
with arbitrary precision by the Solovay-Kitaev algorithm [65], placing it in
the block of $\mathcal{V}_{dt/2}$ before or after the evolution of the
quadratic Hamiltonian since it commutes with the harmonic potential operator.

##  VI Conclusion

To sum up, we propose the GD-based circuit learning to find the time-optimal
control problem, the driving of a quantum particle trapped in a time-varying
harmonic potential, and figure out its quantum speed limit in relation to the
control phase transition. First, we have constructed the digitized quantum
circuit of a time-dependent harmonic oscillator using a finite $n$-qubit
register. Second, we have demonstrated that the learning rate of circuit
optimization can be accelerated by considering various physical quantities,
such as the infidelity, Bures angle, and fidelity susceptibility, as cost
functions, thus reducing training iteration. Third, we have established the
relation between control phase transition and quantum speed limit. Finally, we
have established the error tolerance of our method by considering the presence
of measurement errors in a quantum computer. The absence of a barren plateau
is further justified in our ansatz, enabling the application of VQAs for a
class of tasks that is not affected by the fundamental limitations of NISQ
devices. As a heuristic example, we have demonstrated that quantum control can
be efficiently simulated and optimized using a NISQ device by combining
digital quantum simulation and hybrid circuit learning. Numerical experiments
prove that barren plateaus are avoided in the framework.

## ACKNOWLEDGMENTS

This work was financially supported by NSFC (12075145), STCSM (Grants No.
2019SHZDZX01-ZX04), EU FET Open Grant EPIQUS (899368), the Basque Government
through Grant No. IT1470-22, the project grant PID2021-126273NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” and “ERDF
Invest in your Future”, QUANTEK project (KK-2021/00070), and the BRTAQ project
(expendient no. KK-2022/00041). HTY acknowledges CSC fellowship
(202006890071). X.C. acknowledges ayudas para contratos Ramón y
Cajal–2015-2020 (RYC-2017-22482).

##  Appendix A  The parameter-shift rule

Figure 7: Schematic illustration of a quantum circuit with Trotter step
$N_{t}$. Each dashed-block denoted in Fig. 1) presents the implementation of
time-evolution $e^{-iH(t)dt}$ with quadratic Hamiltonian $H(t)$ (22). Here, an
arbitrary operator $U^{k}(f_{k})$ includes three $f_{k}$-correlated variables:
$\theta_{1}(f_{k})$, $\theta_{2}(f_{k})$ and $\theta_{3}(f_{k})$. The single-
qubit gate with a yellow shadow is used to calculate the gradient of the
objective function.

One recipe to find the partial derivative of an objective function $J(\Theta)$
in parametric quantum circuits (PQCs) is known as the parameter-shift rule
[59, 60]. In general, the expectation value of an observable $\hat{B}$ as a
function of a single parameter $\theta_{k}$ in a circuit is expressed as
$J(\theta_{k})=\langle\hat{B}(\theta_{k})\rangle$. We assume a sequence of
unitary operations represented as $U(\theta_{k})=U_{L}U_{k}(\theta_{k})U_{R}$,
for which we have

$\displaystyle J(\theta_{k})$ $\displaystyle=$ $\displaystyle\langle
0|U_{R}^{\dagger}U_{k}^{\dagger}(\theta_{k})U_{L}^{\dagger}\hat{B}U_{L}U_{k}(\theta_{k})U_{R}|0\rangle$
(28) $\displaystyle=$ $\displaystyle\langle
z|\mathcal{M}(\hat{B},\theta_{k})|z\rangle,$

where
$\mathcal{M}(\hat{B},\theta_{k})=U_{k}^{\dagger}(\theta_{k})U_{L}^{\dagger}\hat{B}U_{L}U_{k}(\theta_{k})$
and the basis reads $|z\rangle=U_{R}|0\rangle$. Consider a unitary operator
$U_{k}(\theta_{k})$ generated by a Pauli matrix $\sigma_{k}$ as
$U_{k}(\theta_{k})=\exp(-i\theta_{k}\sigma_{k}/2)$. The gradient of the
objective function is defined as

$\displaystyle\partial_{\theta_{k}}J(z;\theta_{k})$ $\displaystyle=$
$\displaystyle\langle
z|\partial_{\theta_{k}}\mathcal{M}(\hat{B},\theta_{k})|z\rangle$ (29)
$\displaystyle=$ $\displaystyle c\left[\langle
z|\mathcal{M}(\theta_{k}+s)|z\rangle-\langle
z|\mathcal{M}(\theta_{k}-s)|z\rangle\right],~{}~{}~{}$

where coefficient $c$ and shift $s$ are independent of $\theta_{k}$. The
gradient
$\partial_{\theta_{k}}U_{k}(\theta_{k})=-\frac{i}{2}U_{k}(\theta_{k})\sigma_{k}$,
and inserting it to (28), we have [57]
$\partial_{\theta_{k}}J=-\frac{i}{2}\langle
z|U_{k}^{\dagger}(\theta_{k})[\sigma_{k},\hat{B}]U_{k}(\theta_{k})|z\rangle$.
The commutation relation,
$[\sigma_{k},\hat{B}]=i\left(U_{k}^{\dagger}(\frac{\pi}{2})\hat{B}U_{k}(\frac{\pi}{2})-U_{k}^{\dagger}(-\frac{\pi}{2})\hat{B}U_{k}(-\frac{\pi}{2})\right)$,
enables us derive the analytical gradient as [59, 60]:

$\displaystyle\partial_{\theta_{k}}J$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle
z|\left[U_{k}(\theta_{k}^{+})\hat{B}U_{k}(\theta_{k}^{+})-U_{k}^{\dagger}(\theta_{k}^{-})\hat{B}U_{k}(\theta_{k}^{-})\right]|z\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left[\langle
z|\mathcal{M}(\theta_{k}^{+})|z\rangle-\langle
z|\mathcal{M}(\theta_{k}^{-})|z\rangle\right],$

with $\theta_{k}^{\pm}=\theta_{k}\pm\pi/2$. The above expression provides an
analytical evaluation of the gradient of an objective function involving Pauli
operators.

In Fig. 7, we schematically illustrate the deep circuit of our method, which
is composed of $N_{t}$ dashed blocks (refers to $N_{t}$ Trotter steps). In
each unitary operator $U^{k}(f_{k})$, all rotating parameters
$(\theta_{1},\theta_{2},\theta_{3})$ are $f_{k}$-correlated, as detailed in
Fig. 1. Let us now select an arbitrary single-qubit gate $R(\theta_{k})$, with
a gradient obeying the parameter shift rule in Eq. (A). According to the
algorithm in Fig. 1, the rotation $\theta_{k}$ is correlated with $f_{k}$ as a
liner form: $\theta_{k}=c_{1}\cdot f_{k}$, where $c_{1}$ is a real number.
Starting with Eq. (29), we have
$\partial_{f_{k}}J=c[J(f_{k}+s)-J(f_{k}-s)]=c[J(\theta^{\prime}_{k}+s)-J(\theta^{\prime}_{k}-s)]$,
where the gradient is independent of the initial angle
$\theta^{\prime}_{k}=\theta_{k}/c_{1}$. Consequently, we find
$\partial_{f_{k}}J=2c\cdot\partial_{\theta_{k}}J$, while the shift $f_{k}$ is
$s=\pi/2$. Furthermore, we introduce the notation

$\overline{|c|}=\sum_{i=1}^{N_{r}}\frac{1}{2N_{r}}\left|\frac{\partial_{f_{k}}J}{\partial_{\theta_{k}}J}\right|,$
(30)

where the absolute average value $\overline{|c|}$ over $N_{r}$ random
initialization. In Fig. 8(a), we calculate the average value $\overline{|c|}$
with $N_{r}=50$, which is irrelevant to $n$.

##  Appendix B Comparisons of optimizers

In Appendix (B), we compare the performance of several classical optimizers
for the same optimization task with statistical errors from a finite number of
measurements. We choose widely used optimizers, namely, $\tt SLSQP$, $\tt
COBYLA$ ,$\tt SPSA$ , $\tt BFQS$ based on the library of $\tt qiskit$. In Fig.
8(b), we present the infidelity as a function of the training iteration by
using various optimizers. The SPSA stands out for its performance in an
optimization task in the presence of bit-flip noise.

Figure 8: (a) The coefficient $\overline{|c|}$ as a function of qubit number
$n$. (b) Various training processes for using different optimizers are
illustrated. All calculations are proceeded with the parameters $\beta=0$,
$N_{m}=8192$, and others in Fig. 3.

## References

  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos, Science 334, 57 (2011).
  * Salathé  _et al._ [2015a] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff, Phys. Rev. X 5, 021027 (2015a).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Phys. Rev. X 8, 031022 (2018).
  * O’Malley  _et al._ [2016] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Gogolin  _et al._ [2011] C. Gogolin, M. P. Müller, and J. Eisert, Phys. Rev. Lett. 106, 040401 (2011).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Bañuls  _et al._ [2020] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, The European Physical Journal D 74, 165 (2020).
  * García-Álvarez  _et al._ [2017] L. García-Álvarez, I. L. Egusquiza, L. Lamata, A. del Campo, J. Sonner, and E. Solano, Phys. Rev. Lett. 119, 040501 (2017).
  * Luo  _et al._ [2019] Z. Luo, Y.-Z. You, J. Li, C.-M. Jian, D. Lu, C. Xu, B. Zeng, and R. Laflamme, npj Quantum Information 5, 53 (2019).
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Salathé  _et al._ [2015b] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff, Phys. Rev. X 5, 021027 (2015b).
  * Martinez  _et al._ [2016] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, Nature 534, 516 (2016).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and P. Zoller, Nature 569, 355 (2019).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1411.4028 (2014).
  * Pesah  _et al._ [2021] A. Pesah, M. Cerezo, S. Wang, T. Volkoff, A. T. Sornborger, and P. J. Coles, Phys. Rev. X 11, 041011 (2021).
  * Magann  _et al._ [2021] A. B. Magann, C. Arenz, M. D. Grace, T.-S. Ho, R. L. Kosut, J. R. McClean, H. A. Rabitz, and M. Sarovar, PRX Quantum 2, 010101 (2021).
  * Suzuki [1976] M. Suzuki, Communications in Mathematical Physics 51, 183 (1976).
  * Suzuki [1993] M. Suzuki, Proceedings of the Japan Academy, Series B 69, 161 (1993).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Suzuki [1990] M. Suzuki, Physics Letters A 146, 319 (1990).
  * Las Heras  _et al._ [2016] U. Las Heras, U. Alvarez-Rodriguez, E. Solano, and M. Sanz, Phys. Rev. Lett. 116, 230504 (2016).
  * McArdle  _et al._ [2019] S. McArdle, X. Yuan, and S. Benjamin, Phys. Rev. Lett. 122, 180501 (2019).
  * Liu  _et al._ [2020] Y.-X. Liu, J. Hines, Z. Li, A. Ajoy, and P. Cappellaro, Phys. Rev. A 102, 010601 (2020).
  * Bolens and Heyl [2021] A. Bolens and M. Heyl, Phys. Rev. Lett. 127, 110502 (2021).
  * Fösel  _et al._ [2018] T. Fösel, P. Tighineanu, T. Weiss, and F. Marquardt, Phys. Rev. X 8, 031084 (2018).
  * Bukov [2018] M. Bukov, Phys. Rev. B 98, 224305 (2018).
  * Albarrán-Arriagada  _et al._ [2018] F. Albarrán-Arriagada, J. C. Retamal, E. Solano, and L. Lamata, Phys. Rev. A 98, 042315 (2018).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * Giuliano Benenti [2007] G. S. Giuliano Benenti, arXiv , 0709.1704 (2007).
  * Ollitrault  _et al._ [2020] P. J. Ollitrault, G. Mazzola, and I. Tavernelli, Phys. Rev. Lett. 125, 260511 (2020).
  * Bukov  _et al._ [2018] M. Bukov, A. G. R. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Phys. Rev. X 8, 031086 (2018).
  * Chen  _et al._ [2010] X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. Lett. 104, 063002 (2010).
  * Li  _et al._ [2011] Y. Li, L.-A. Wu, and Z. D. Wang, Phys. Rev. A 83, 043804 (2011).
  * García-Ripoll  _et al._ [2020] J. García-Ripoll, A. Ruiz-Chamorro, and E. Torrontegui, Phys. Rev. Applied 14, 044035 (2020).
  * Kosloff [1988] R. Kosloff, The Journal of Physical Chemistry 92, 2087 (1988).
  * Stefanatos  _et al._ [2010] D. Stefanatos, J. Ruths, and J.-S. Li, Phys. Rev. A 82, 063422 (2010).
  * Hoffmann  _et al._ [2011] K. Hoffmann, P. Salamon, Y. Rezek, and R. Kosloff, EPL (Europhysics Letters) 96, 60015 (2011).
  * Huang  _et al._ [2020] T.-Y. Huang, B. A. Malomed, and X. Chen, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 053131 (2020).
  * Dupays  _et al._ [2021] L. Dupays, D. C. Spierings, A. M. Steinberg, and A. del Campo, Phys. Rev. Research 3, 033261 (2021).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Journal of Physics A: Mathematical and Theoretical 50, 453001 (2017).
  * Caneva  _et al._ [2009] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Phys. Rev. Lett. 103, 240501 (2009).
  * Funo  _et al._ [2017] K. Funo, J.-N. Zhang, C. Chatou, K. Kim, M. Ueda, and A. del Campo, Phys. Rev. Lett. 118, 100602 (2017).
  * Campbell and Deffner [2017] S. Campbell and S. Deffner, Phys. Rev. Lett. 118, 100601 (2017).
  * Mandelstam and Tamm [1945] L. Mandelstam and I. Tamm, J. Phys. USSR 9, 249 (1945).
  * Margolus and Levitin [1998] N. Margolus and L. B. Levitin, Physica D: Nonlinear Phenomena 120, 188 (1998).
  * Ness  _et al._ [2021] G. Ness, M. R. Lam, W. Alt, D. Meschede, Y. Sagi, and A. Alberti, Science Advances 7, eabj9119 (2021).
  * Anandan and Aharonov [1990] J. Anandan and Y. Aharonov, Phys. Rev. Lett. 65, 1697 (1990).
  * Uhlmann [1992] A. Uhlmann, Physics Letters A 161, 329 (1992).
  * Deffner and Lutz [2013] S. Deffner and E. Lutz, Journal of Physics A: Mathematical and Theoretical 46, 335302 (2013).
  * You  _et al._ [2007] W.-L. You, Y.-W. Li, and S.-J. Gu, Phys. Rev. E 76, 022101 (2007).
  * Gu  _et al._ [2008] S.-J. Gu, H.-M. Kwok, W.-Q. Ning, and H.-Q. Lin, Phys. Rev. B 77, 245109 (2008).
  * GU [2010] S.-J. GU, International Journal of Modern Physics B 24, 4371 (2010), https://doi.org/10.1142/S0217979210056335 .
  * Boggs and Tolle [1995] P. T. Boggs and J. W. Tolle,  _Sequential Quadratic Programming_, Vol. 4 (Cambridge University Press, 1995) p. 1–51.
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, and etc., Nature Methods 17, 261 (2020).
  * McClean  _et al._ [2018] J. McClean, S. Boixo, V. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Nature Reviews Physics 3, 625 (2021).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Phys. Rev. A 98, 032309 (2018).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Phys. Rev. A 99, 032331 (2019).
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, PRX Quantum 3, 010313 (2022).
  * Spall [1998] J. C. Spall, Johns Hopkins Apl Technical Digest 19, 482 (1998).
  * Bottou [2010] L. Bottou,  _Proceedings of COMPSTAT’2010_ , edited by Y. Lechevallier and G. Saporta (Physica-Verlag HD, Heidelberg, 2010) pp. 177–186.
  * Marin-Sanchez  _et al._ [2021] G. Marin-Sanchez, J. Gonzalez-Conde, and M. Sanz, arXiv preprint arXiv:2111.07933 (2021).
  * Kitaev [1997] A. Y. Kitaev, Russian Mathematical Surveys 52, 1191 (1997).


# Sliding-Mode Nash Equilibrium Seeking for a Quadratic Duopoly Game

Victor Hugo Pereira Rodrigues1, Tiago Roux Oliveira1, Miroslav Krstić2, Tamer
Başar3 1V. H. P. Rodrigues and T. R. Oliveira are with the Dept. of Eletronics
and Telecommunication Engineering (DETEL), State University of Rio de Janeiro
(UERJ), Rio de Janeiro – RJ, Brazil. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>Krstić is with the Dept. of Mechanical and Aerospace
Engineering, University of California – San Diego (UCSD), La Jolla – CA, USA.
Email<EMAIL_ADDRESS>Başar is with the Dept. of Electrical and Computer
Engineering, University of Illinois Urbana-Champaign, Urbana – IL, USA. Email:
<EMAIL_ADDRESS>
###### Abstract

This paper introduces a new method to achieve stable convergence to Nash
equilibrium in duopoly noncooperative games. Inspired by the recent fixed-time
Nash Equilibrium seeking (NES) [25] as well as prescribed-time extremum
seeking (ES) [29] and source seeking [28] schemes, our approach employs a
distributed sliding mode control (SMC) scheme, integrating extremum seeking
with sinusoidal perturbation signals to estimate the pseudogradients of
quadratic payoff functions. Notably, this is the first attempt to address
noncooperative games without relying on models, combining classical extremum
seeking with relay components instead of proportional control laws. We prove
finite-time convergence of the closed-loop average system to Nash equilibrium
using stability analysis techniques such as time-scaling, Lyapunov’s direct
method, and averaging theory for discontinuous systems. Additionally, we
quantify the size of residual sets around the Nash equilibrium and validate
our theoretical results through simulations.

##  I INTRODUCTION

Game theory provides a conceptual framework for examining social interactions
among competitive participants, utilizing mathematical models to grasp
strategic decision-making [10, 7]. Its utility extends across multiple domains
such as engineering systems, biological behaviors, and financial markets,
rendering it indispensable across diverse disciplines [12, 1, 8]. Considerable
research has been conducted on differential games, encompassing both
theoretical exploration and practical applications.

Games are commonly categorized as cooperative and noncooperative [6]. In
cooperative games, players establish enforceable agreements, whereas
noncooperative games center on individual player actions and Nash equilibria
[20]. Nash equilibrium, a pivotal concept in noncooperative game theory,
characterizes a situation where no player can improve its payoff independently
of other players in the game [6].

Efforts towards achieving convergence to Nash equilibrium have persisted for
decades [17, 5], with research delving into learning-based update schemes
[30]. Recent investigations have delved into extremum seeking (ES) methods for
real-time computation of Nash equilibria in noncooperative games [9],
facilitating stable convergence in the absence of any model information [16].

From the literature of ES based on sliding modes, we can mention the
publication from Özgüner and coauthors [23], where a periodic switching
function was introduced to provide a state-feedback control approach.
Following that, Oliveira et al. [21, 22] have expanded the application of this
previous methodology by considering only output feedback. Notably, the same
team has also introduced a different sliding mode-based ES controller using
monitoring functions [2, 3]. However, none of such ES control strategies
followed the original and simplest structure from the scheme with sinusoidal
periodic perturbations [16].

On the other hand, recently some interesting fixed-time results for Nash
Equilibrium seeking [25] (time-invariant, non-smooth) and prescribed-time
schemes (time-varying, smooth) for ES [29] and source seeking [28] have been
provided. In particular, the paper [25] proposes control laws which seem to be
similar to those in higher-order SMC. For instance, when considering second-
order SMC, the most known is the super-twisting algorithm (STA) [19]. Indeed,
there is an entire family of finite-time (or fixed-time) SMC schemes [18]
which could be explored for extremum seeking. In this context, bringing the
idea of a variable structure/sliding mode NES in finite (not fixed) time to
the attention of the large SMC community is still missing and we believe has
an important value. We have decided to start with the simplest (and more
popular) case: first-order SMC (FOSMC), simply a relay function. FOSMC is
discontinuous but can assure finite-time convergence.

To bridge the existing gaps, our paper extends the classical ES algorithm [16]
by incorporating relay components in place of the usual proportional control
laws with respect to the pseudogradient111The pseudogradient is essentially a
generalized notion of gradient that provides directional information allowing
the analysis of equilibrium points in complex strategic interactions of
concave $N$-players game. Hence, given a vector-valued function
$J(\theta):\mathbb{R}^{N}\to\mathbb{R}^{N}$, we define its pseudogradient $G$
as the vector-valued function $G(\theta):=\left[\frac{\partial
J_{1}}{\partial\theta_{1}}\,,\frac{\partial
J_{2}}{\partial\theta_{2}}\,,\ldots\,,\frac{\partial
J_{N}}{\partial\theta_{N}}\right]^{T}$. For more details, see [26, Eq. (3.9)].
estimates in order to derive the first sliding mode version of the classical
ES algorithm based on sinusoidal periodic perturbations.

We address the NES problem in two-player noncooperative games, employing the
proposed distributed variable structure ES controller. Despite the inherent
challenges of decentralized control in games, our approach ensures convergence
to Nash equilibrium when both players adopt the proposed ES algorithms via
sliding modes. We provide theoretical analysis for duopoly games, leveraging
time-scaling techniques, constructing a Lyapunov function, and employing
averaging methods to ensure closed-loop stability. Numerical simulations
demonstrate the efficacy of our approach, highlighting its superiority over
classical ES with linear proportional control laws [16] since now the average
closed-loop system is finite-time stable rather than only exponentially
stable.

##  II Duopoly Game with Quadratic Payoffs: General Formulation

In a duopoly game, the optimal outcomes for players P1 and P2, denoted by
$y_{1}(t)\in\mathbb{R}$ and $y_{2}(t)\in\mathbb{R}$, respectively, are not
solely determined by their individual actions or decision strategies,
represented by $\theta_{1}(t)\in\mathbb{R}$ and $\theta_{2}(t)\in\mathbb{R}$,
but also by their opponents’. Hence, the pair of decisions,
$\theta=[\theta_{1}\,,\theta_{2}]^{T}\in\mathbb{R}^{2}$, influences the payoff
functions $J_{1}$ and $J_{2}$ of the two players. An optimal choice
($\theta(t)\equiv\theta^{*}$), called Nash equilibrium [10], captures the
situation where no player benefits in terms of improvement of its own payoff
by deviating unilaterally from its own policy.

In particular, we examine duopoly games where each player’s payoff function
takes a quadratic form. This can be expressed for each player as a strictly
concave function of its individual action:

$\displaystyle J_{1}(\theta(t))$
$\displaystyle=\frac{H_{11}^{1}}{2}\theta_{1}^{2}(t)+\frac{H_{22}^{1}}{2}\theta_{2}^{2}(t)+H_{12}^{1}\theta_{1}(t)\theta_{2}(t)$
$\displaystyle\quad+h_{1}^{1}\theta_{1}(t)+h_{2}^{1}\theta_{2}(t)+c_{1},$ (1)
$\displaystyle J_{2}(\theta(t))$
$\displaystyle=\frac{H_{11}^{2}}{2}\theta_{1}^{2}(t)+\frac{H_{22}^{2}}{2}\theta_{2}^{2}(t)+H_{21}^{2}\theta_{1}(t)\theta_{2}(t)$
$\displaystyle\quad+h_{1}^{2}\theta_{1}(t)+h_{2}^{2}\theta_{2}(t)+c_{2},$ (2)

Here, $J_{1}(\theta),J_{2}(\theta):\mathbb{R}^{2}\to\mathbb{R}$ represent the
payoff functions for players 1 and 2, respectively, while
$\theta_{1}(t),\theta_{2}(t)\in\mathbb{R}$ denote the players’ decision
variables, and $H_{jk}^{i}$, $h_{j}^{i}$, $c_{i}\in\mathbb{R}$ are constants,
with $H_{ii}^{i}<0$, for all $i,j,k\in\\{1,2\\}$. In the analysis of the
duopoly game, the individual payoff functions exhibit strict concavity only in
self-action, deviating from concavity in both actions. This implies the
absence of a quadratic term in the opponent’s action, leading to
$H^{i}_{jj}=0$ for all $j\neq i$, alongside the general non-zero nature of all
other coefficients except for $c_{i}$, which are typically absent.

For the sake completeness, let us define the mathematical representation of
the Nash equilibrium
$\theta^{\ast}=[\theta^{\ast}_{1}\,,\theta^{\ast}_{2}]^{T}$ in a two-player
game as below, for all $\theta_{1}\,,\theta_{2}$:

$\displaystyle J_{1}(\theta_{1}^{\ast}\,,\theta_{2}^{\ast})\geq
J_{1}(\theta_{1}\,,\theta_{2}^{\ast})~{}~{}\mbox{and}~{}~{}~{}J_{2}(\theta_{1}^{\ast}\,,\theta_{2}^{\ast})\geq
J_{2}(\theta_{1}^{\ast}\,,\theta_{2}).\\!\\!$ (3)

Thus, no player has an incentive to independently deviate from $\theta^{*}$.
In the duopoly context, $\theta_{1}$ and $\theta_{2}$ are elements of
$\mathbb{R}$, representing the set of real numbers. In order to determine the
Nash equilibrium solution in strictly concave quadratic games involving two
players, where each action set is the entire real line, one would
differentiate $J_{i}$ with respect to $\theta_{i}(t)\,,\forall i=1,2$, setting
the resulting expressions equal to zero, and solving the set of equations thus
obtained. This set of equations, which also provides a sufficient condition
due to the strict concavity [26], is

$\displaystyle\begin{cases}H_{11}^{1}\theta_{1}^{*}+H_{12}^{1}\theta_{2}^{*}+h_{1}^{1}=0\\\
H_{21}^{2}\theta_{1}^{*}+H_{22}^{2}\theta_{2}^{*}+h_{2}^{2}=0\\\
\end{cases}\,,$ (4)

which can be written in the matrix form as

$\displaystyle\begin{bmatrix}H_{11}^{1}&H_{12}^{1}\\\
H_{21}^{2}&H_{22}^{2}\end{bmatrix}\begin{bmatrix}\theta_{1}^{*}\\\
\theta_{2}^{*}\end{bmatrix}=-\begin{bmatrix}h_{1}^{1}\\\
h_{2}^{2}\end{bmatrix}\,.$ (5)

Defining the Hessian matrix $H$, and vectors $\theta^{*}$ and $h$ by

$\displaystyle H:=\begin{bmatrix}H_{11}^{1}&H_{12}^{1}\\\
H_{21}^{2}&H_{22}^{2}\end{bmatrix}\,,\quad\theta^{*}:=\begin{bmatrix}\theta_{1}^{*}\\\
\theta_{2}^{*}\end{bmatrix}\,,\quad h:=\begin{bmatrix}h_{1}^{1}\\\
h_{2}^{2}\end{bmatrix}\,,$ (6)

there exists a unique Nash Equilibrium at $\theta^{*}\\!=\\!-H^{-1}h$, if $H$
is invertible (which holds in any realistic model of the duopoly game [6,
Chapter 4]):

$\displaystyle\begin{bmatrix}\theta_{1}^{*}\\\
\theta_{2}^{*}\end{bmatrix}=-\begin{bmatrix}H_{11}^{1}&H_{12}^{1}\\\
H_{21}^{2}&H_{22}^{2}\end{bmatrix}^{-1}\begin{bmatrix}h_{1}^{1}\\\
h_{2}^{2}\end{bmatrix}\,.$ (7)

###  II-A Nash Equilibrium Seeking

The primary goal of this paper is to develop distributed variable structure
policies via sliding modes to obtain NES in noncooperative duopoly games.
Throughout our analysis, it is assumed that $H_{ii}^{i}<0$ for both $i$ (a
natural consequence of strict concavity), the matrix $H$ in (6) has full rank,
and all the parameters are unknown.

Figure 1: Block diagram illustrating the NES strategy through distributed
event-triggered control policies used by each player.

The schematic diagram depicted in Fig. 1 provides an overview of the proposed
NES policy for each player, detailing their respective outputs as follows:

$\displaystyle\begin{cases}y_{1}(t)&=J_{1}(\theta(t))\,,\\\
y_{2}(t)&=J_{2}(\theta(t))\,.\end{cases}$ (8)

The probing signals are

$\displaystyle\begin{cases}S_{1}(t)&=a_{1}\sin(\omega_{1}t)\\\
S_{2}(t)&=a_{2}\sin(\omega_{2}t)\end{cases}\,,$ (9)

and the demodulating signals are represented by

$\displaystyle\begin{cases}M_{1}(t)&=\frac{2}{a_{1}}\sin(\omega_{1}t)\\\
M_{2}(t)&=\frac{2}{a_{2}}\sin(\omega_{2}t)\end{cases}\,,$ (10)

having non-zero constant amplitudes $a_{1}$ and $a_{2}$, both greater than
zero, and frequencies $\omega_{1}$ and $\omega_{2}$, satisfying
$\omega_{1}\neq\omega_{2}$. These probing frequencies $\omega_{i}$ can be
chosen as

$\omega_{i}=\omega_{i}^{\prime}\omega=\mathcal{O}(\omega)\,,\quad i=1\mbox{ or
}2\,,$ (11)

where $\omega$ represents a positive constant and $\omega_{i}^{\prime}$ is a
rational number. One potential selection is detailed in [11].

If we denote by $\hat{\theta}_{1}(t)$ and $\hat{\theta}_{2}(t)$ the estimates
of $\theta^{\ast}_{1}$ and $\theta^{\ast}_{2}$, respectively, we can define
the “estimation error” as follows:

$\displaystyle\begin{cases}\tilde{\theta}_{1}(t)=\hat{\theta}_{1}(t)-\theta_{1}^{*}\\\
\tilde{\theta}_{2}(t)=\hat{\theta}_{2}(t)-\theta_{2}^{*}\end{cases}\,.$ (12)

The estimation of the pseudogradient of the unknown payoff functions

$\displaystyle\begin{cases}\hat{G}_{1}(t)&=M_{1}(t)y_{1}(t)\\\
\hat{G}_{2}(t)&=M_{2}(t)y_{2}(t)\end{cases}\,,$ (13)

can be rewritten as

$\displaystyle\hat{G}_{1}(t)$
$\displaystyle=H^{1}_{11}\theta^{\ast}_{1}+H^{1}_{12}\theta^{\ast}_{2}+h^{1}_{1}$
$\displaystyle\quad+\mathcal{H}_{11}(t)\tilde{\theta}_{1}(t)+\mathcal{H}_{12}(t)\tilde{\theta}_{2}(t)$
$\displaystyle\quad+\frac{H^{1}_{11}}{a_{1}}\sin(\omega_{1}t)\tilde{\theta}_{1}^{2}(t)+\frac{H^{1}_{22}}{a_{1}}\sin(\omega_{1}t)\tilde{\theta}_{2}^{2}(t)$
$\displaystyle\quad+\frac{2H^{1}_{12}}{a_{1}}\sin(\omega_{1}t)\tilde{\theta}_{1}(t)\tilde{\theta}_{2}(t)+\delta_{1}(t)\,,$
(14) $\displaystyle\hat{G}_{2}(t)$
$\displaystyle=H^{2}_{21}\theta^{\ast}_{1}+H^{2}_{22}\theta^{\ast}_{2}+h^{2}_{2}$
$\displaystyle\quad+\mathcal{H}_{21}(t)\tilde{\theta}_{1}(t)+\mathcal{H}_{22}(t)\tilde{\theta}_{2}(t)$
$\displaystyle\quad+\frac{H^{2}_{11}}{a_{2}}\sin(\omega_{2}t)\tilde{\theta}_{1}^{2}(t)+\frac{H^{2}_{22}}{a_{2}}\sin(\omega_{2}t)\tilde{\theta}_{2}^{2}(t)$
$\displaystyle\quad+\frac{2H^{2}_{21}}{a_{2}}\sin(\omega_{2}t)\tilde{\theta}_{1}(t)\tilde{\theta}_{2}(t)+\delta_{2}(t)\,,$
(15)

with the time-varying parameters given by

$\displaystyle\mathcal{H}_{11}(t)$
$\displaystyle=H^{1}_{11}-H^{1}_{11}\cos(2\omega_{1}t)$
$\displaystyle\quad+\frac{2h^{1}_{1}}{a_{1}}\sin(\omega_{1}t)+\frac{2H^{1}_{11}\theta^{\ast}_{1}}{a_{1}}\sin(\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{12}a_{2}}{a_{1}}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad+\frac{2H^{1}_{12}\theta^{\ast}_{2}}{a_{1}}\sin(\omega_{1}t)+\frac{H^{1}_{12}a_{2}}{a_{1}}\cos[(\omega_{1}-\omega_{2})t]\,,$
(16) $\displaystyle\mathcal{H}_{12}(t)$
$\displaystyle=H^{1}_{12}-H^{1}_{12}\cos(2\omega_{1}t)$
$\displaystyle\quad+\frac{2h^{1}_{2}}{a_{1}}\sin(\omega_{1}t)+\frac{2H^{1}_{12}\theta^{\ast}_{1}}{a_{1}}\sin(\omega_{1}t)$
$\displaystyle\quad+\frac{2H^{1}_{22}\theta^{\ast}_{2}}{a_{1}}\sin(\omega_{1}t)+\frac{H^{1}_{22}a_{2}}{a_{1}}\cos[(\omega_{1}-\omega_{2})t]$
$\displaystyle\quad+\frac{H^{1}_{22}a_{2}}{a_{1}}\cos[(\omega_{1}+\omega_{2})t]\,,$
(17) $\displaystyle\mathcal{H}_{21}(t)$
$\displaystyle=H^{2}_{21}-H^{2}_{21}\cos(2\omega_{2}t)$
$\displaystyle\quad+\frac{2h^{2}_{1}}{a_{2}}\sin(\omega_{2}t)+\frac{2H^{2}_{21}\theta^{\ast}_{1}}{a_{2}}\sin(\omega_{2}t)$
$\displaystyle\quad+\frac{2H^{2}_{22}\theta^{\ast}_{2}}{a_{2}}\sin(\omega_{2}t)+\frac{H^{2}_{22}a_{1}}{a_{2}}\cos[(\omega_{2}-\omega_{1})t]$
$\displaystyle\quad+\frac{H^{2}_{22}a_{1}}{a_{2}}\cos[(\omega_{2}+\omega_{1})t]\,,$
(18) $\displaystyle\mathcal{H}_{22}(t)$
$\displaystyle=H^{2}_{22}-H^{2}_{22}\cos(2\omega_{2}t)$
$\displaystyle\quad+\frac{2h^{2}_{2}}{a_{2}}\sin(\omega_{2}t)+\frac{2H^{2}_{22}\theta^{\ast}_{2}}{a_{2}}\sin(\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{21}a_{1}}{a_{2}}\cos[(\omega_{2}+\omega_{1})t]$
$\displaystyle\quad+\frac{2H^{2}_{21}\theta^{\ast}_{1}}{a_{2}}\sin(\omega_{2}t)+\frac{H^{2}_{21}a_{1}}{a_{2}}\cos[(\omega_{2}-\omega_{1})t]\,,$
(19) $\displaystyle\delta_{1}(t)$
$\displaystyle=\frac{H^{1}_{11}a_{1}}{2}\sin(\omega_{1}t)+\frac{H^{1}_{11}a_{1}}{4}\cos(\omega_{1}t)+\frac{H^{1}_{11}a_{1}}{4}\cos(3\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{22}a_{2}^{2}}{2a_{1}}\sin(\omega_{1}t)-\frac{H^{1}_{22}a_{2}^{2}}{2a_{1}}\cos[(\omega_{1}-2\omega_{2})t]$
$\displaystyle\quad+\frac{H^{1}_{22}a_{2}^{2}}{2a_{1}}\cos[(\omega_{1}+2\omega_{2})t]+\frac{2H^{1}_{12}\theta^{\ast}_{1}\theta^{\ast}_{2}}{a_{1}}\sin(\omega_{1}t)$
$\displaystyle\quad-h^{1}_{1}\cos(2\omega_{1}t)$
$\displaystyle\quad+\frac{a_{2}h^{1}_{2}}{a_{1}}\cos[(\omega_{1}-\omega_{2})t]+\frac{a_{2}h^{1}_{2}}{a_{1}}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad-H^{1}_{11}\theta^{\ast}_{1}\cos(2\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{12}a_{2}\theta^{\ast}_{1}}{a_{1}}\cos[(\omega_{1}-\omega_{2})t]+\frac{H^{1}_{12}a_{2}\theta^{\ast}_{1}}{a_{1}}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad-H^{1}_{12}\theta^{\ast}_{2}\cos(2\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{22}a_{2}\theta^{\ast}_{2}}{a_{1}}\cos[(\omega_{1}-\omega_{2})t]-\frac{H^{1}_{22}a_{2}\theta^{\ast}_{2}}{a_{1}}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad+\frac{2c_{1}}{a_{1}}\sin(\omega_{1}t)+\frac{2h^{1}_{1}\theta^{\ast}_{1}}{a_{1}}\sin(\omega_{1}t)+\frac{2h^{1}_{2}\theta^{\ast}_{2}}{a_{1}}\sin(\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{11}\theta^{\ast
2}_{1}}{a_{1}}\sin(\omega_{1}t)+\frac{H^{1}_{22}\theta^{\ast
2}_{2}}{a_{1}}\sin(\omega_{1}t)$
$\displaystyle\quad+\frac{H^{1}_{12}a_{2}}{2}\cos(\omega_{2}t)-\frac{H^{1}_{12}a_{2}}{2}\cos[(2\omega_{1}-\omega_{2})t]$
$\displaystyle\quad+\frac{H^{1}_{12}a_{2}}{2}\cos(\omega_{2}t)+\frac{H^{1}_{12}a_{2}}{2}\cos[(2\omega_{1}+\omega_{2})t]\,,$
(20)

$\displaystyle\delta_{2}(t)$
$\displaystyle=\frac{H^{2}_{11}a_{2}}{2}\sin(\omega_{2}t)+\frac{H^{2}_{11}a_{2}}{4}\cos(\omega_{2}t)+\frac{H^{2}_{11}a_{2}}{4}\cos(3\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{22}a_{2}}{2}\sin(\omega_{2}t)-\frac{H^{2}_{22}a_{2}}{2}\cos[(\omega_{2}-2\omega_{1})t]$
$\displaystyle\quad+\frac{H^{2}_{22}a_{2}}{2}\cos[(\omega_{2}+2\omega_{1})t]+\frac{2H^{2}_{21}\theta^{\ast}_{1}\theta^{\ast}_{2}}{a_{2}}\sin(\omega_{2}t)$
$\displaystyle\quad-h^{2}_{2}\cos(2\omega_{2}t)$
$\displaystyle\quad+h^{1}_{2}\cos[(\omega_{2}-\omega_{1})t]+\frac{h^{2}_{2}}{a_{2}}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad-H^{2}_{11}\theta^{\ast}_{1}\cos(2\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{21}a_{2}\theta^{\ast}_{1}}{a_{1}}\cos[(\omega_{2}-\omega_{1})t]+H^{2}_{21}\theta^{\ast}_{1}\cos[(\omega_{1}+\omega_{2})t]$
$\displaystyle\quad-H^{2}_{21}\theta^{\ast}_{2}\cos(2\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{22}a_{1}\theta^{\ast}_{2}}{a_{2}}\cos[(\omega_{2}-\omega_{1})t]-\frac{H^{2}_{22}a_{1}\theta^{\ast}_{2}}{a_{2}}\cos[(\omega_{2}+\omega_{1})t]$
$\displaystyle\quad+\frac{2c_{2}}{a_{2}}\sin(\omega_{2}t)+\frac{2h^{2}_{1}\theta^{\ast}_{1}}{a_{2}}\sin(\omega_{2}t)+\frac{2h^{2}_{2}\theta^{\ast}_{2}}{a_{2}}\sin(\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{11}\theta^{\ast
2}_{1}}{a_{2}}\sin(\omega_{2}t)+\frac{H^{2}_{22}\theta^{\ast
2}_{2}}{a_{2}}\sin(\omega_{2}t)$
$\displaystyle\quad+\frac{H^{2}_{21}a_{1}}{2}\cos(\omega_{1}t)-\frac{H^{2}_{21}a_{1}}{2}\cos[(2\omega_{2}-\omega_{1})t]$
$\displaystyle\quad+\frac{H^{2}_{21}a_{1}}{2}\cos(\omega_{1}t)+\frac{H^{2}_{21}a_{1}}{2}\cos[(2\omega_{2}+\omega_{1})t]\,.$
(21)

It is worth noting that the quadratic terms in (14) and (15) can be
disregarded in a local analysis [15]. Moreover, from (4),
$H_{11}^{1}\theta_{1}^{*}+H_{12}^{1}\theta_{2}^{*}+h_{1}^{1}=0$ and
$H_{21}^{2}\theta_{1}^{*}+H_{22}^{2}\theta_{2}^{*}+h_{2}^{2}=0$. Therefore,
the pseudogradient estimates are locally expressed as

$\displaystyle\hat{G}_{1}(t)$
$\displaystyle=\mathcal{H}_{11}(t)\tilde{\theta}_{1}(t)+\mathcal{H}_{12}(t)\tilde{\theta}_{2}(t)+\delta_{1}(t)\,,$
(22) $\displaystyle\hat{G}_{2}(t)$
$\displaystyle=\mathcal{H}_{21}(t)\tilde{\theta}_{1}(t)+\mathcal{H}_{22}(t)\tilde{\theta}_{2}(t)+\delta_{2}(t)\,.$
(23)

Now, considering the time-derivative of (12) and referring to the NES scheme
illustrated in Fig. 1, we can derive the dynamics governing $\hat{\theta}(t)$
and $\tilde{\theta}(t)$ as follows:

$\displaystyle\frac{d\tilde{\theta}_{1}(t)}{dt}$
$\displaystyle=\frac{d\hat{\theta}_{1}(t)}{dt}=u_{1}(t)\,,$ (24)
$\displaystyle\frac{d\tilde{\theta}_{2}(t)}{dt}$
$\displaystyle=\frac{d\hat{\theta}_{2}(t)}{dt}=u_{2}(t)\,,$ (25)

where $u_{1}(t)$ and $u_{2}(t)$ represent the NES control laws to be
formulated.

###  II-B Distributed Variable Structure Actions for Nash Equilibrium Seeking
via Sliding Modes

In the linear approach, the continuous-time feedback law structures remain
fixed as

$\displaystyle\begin{cases}u_{1}(t)&=K_{1}\hat{G}_{1}(t)\\\
u_{2}(t)&=K_{2}\hat{G}_{2}(t)\end{cases}\,,\quad\text{for all }t\geq 0\,.$
(26)

with the constants $K_{1}$ and $K_{2}$ chosen such that

$\displaystyle K=\begin{bmatrix}K_{1}&0\\\ 0&K_{2}\end{bmatrix}\,,$ (27)

ensuring that $KH$ is Hurwitz.

In variable structure systems (VSS), the control is allowed to dynamically
change its structure, switching between various continuous functions of the
state. The challenge in variable structure design lies in selecting the
parameters for each structure and defining the switching logic. This concept
of structural flexibility enables the creation of novel system behaviors by
combining segments of trajectories from various structures. Additionally, a
fundamental aspect of variable structure systems is the ability to generate
trajectories that are not naturally present in any individual structure. These
trajectories define a unique form of motion known as the sliding mode.

In this paper, we design variable structure actions for Nash Equilibrium
Seeking through the following sliding mode control laws:

$\displaystyle\begin{cases}u_{1}&=K_{1}\text{sgn}(\hat{G}_{1}(t))\\\
u_{2}&=K_{2}\text{sgn}(\hat{G}_{2}(t))\end{cases}\,.$ (28)

###  II-C Closed-loop System

Plugging the proposed sliding mode based control laws in (28) into the time-
derivatives of (22), (23), (24) and equation (25) yields, for all $t\geq 0$,
the following dynamics governing $\hat{G}(t)$ and $\tilde{\theta}(t)$:

$\displaystyle\frac{d\hat{G}_{1}(t)}{dt}$
$\displaystyle=\mathcal{H}_{11}(t)K_{1}\text{sgn}(\hat{G}_{1}(t))+\mathcal{H}_{12}(t)K_{2}\text{sgn}(\hat{G}_{2}(t))+$
$\displaystyle\quad+\frac{d\mathcal{H}_{11}(t)}{dt}\tilde{\theta}_{1}(t)+\frac{d\mathcal{H}_{12}(t)}{dt}\tilde{\theta}_{2}(t)+\frac{d\delta_{1}(t)}{dt}\,,$
(29) $\displaystyle\frac{d\hat{G}_{2}(t)}{dt}$
$\displaystyle=\mathcal{H}_{21}(t)K_{1}\text{sgn}(\hat{G}_{1}(t))+\mathcal{H}_{22}(t)K_{2}\text{sgn}(\hat{G}_{2}(t))+$
$\displaystyle\quad+\frac{d\mathcal{H}_{21}(t)}{dt}\tilde{\theta}_{1}(t)+\frac{d\mathcal{H}_{22}(t)}{dt}\tilde{\theta}_{2}(t)+\frac{d\delta_{2}(t)}{dt}\,,$
(30) $\displaystyle\frac{d\tilde{\theta}_{1}(t)}{dt}$
$\displaystyle=K_{1}\text{sgn}(\hat{G}_{1}(t))\,,$ (31)
$\displaystyle\frac{d\tilde{\theta}_{2}(t)}{dt}$
$\displaystyle=K_{2}\text{sgn}(\hat{G}_{2}(t))\,.$ (32)

Hence, in a concise vector-form representation, we have

$\displaystyle\frac{d\hat{G}(t)}{dt}$
$\displaystyle=\mathcal{H}(t)K\text{sgn}(\hat{G}(t))+\frac{d\mathcal{H}(t)}{dt}\tilde{\theta}(t)+\frac{d\delta(t)}{dt}\,,$
(33) $\displaystyle\frac{d\tilde{\theta}(t)}{dt}$
$\displaystyle=K\text{sgn}(\hat{G}(t))\,,$ (34)

where

$\displaystyle\text{sgn}(\hat{G}(t))$
$\displaystyle:=\begin{bmatrix}\text{sgn}(\hat{G}_{1}(t))\\\
\text{sgn}(\hat{G}_{2}(t))\end{bmatrix}\,,\quad K:=\begin{bmatrix}K_{1}&0\\\
0&K_{2}\end{bmatrix}\,,$ $\displaystyle\mathcal{H}(t)$
$\displaystyle:=\begin{bmatrix}\mathcal{H}_{11}(t)&\mathcal{H}_{12}(t)\\\
\mathcal{H}_{21}(t)&\mathcal{H}_{22}(t)\end{bmatrix}\,,\mbox{ and
}\delta(t):=\begin{bmatrix}\delta_{1}(t)\\\ \delta_{2}(t)\end{bmatrix}\,.$
(35)

The closed-loop system described by (33) and (34) highlights a crucial point:
even if the product $\mathcal{H}(t)K$ on averaging sense would form a Hurwitz
matrix, the convergence to the equilibrium $\hat{G}\equiv 0$ and
$\tilde{\theta}\equiv 0$ would not be directly guaranteed through the
$\tilde{\theta}$–dynamics. On the other hand, it is important to note that the
disturbances $\delta(t)$ and $\frac{d\delta(t)}{dt}$ as well as the time-
varying matrix $\frac{d\mathcal{H}(t)}{dt}$ possess zero mean values, which
will be helpful in the following averaging analysis.

###  II-D Time-Scaling Procedure

To facilitate the stability analysis of the closed-loop system, we introduce a
convenient time scale. By examining (11), it is evident that the dither
frequencies (9) and (10), as well as their combinations, are rational.
Additionally, there exists a time period $T$ such that

$\displaystyle T$
$\displaystyle=2\pi\times\text{LCM}\left\\{\frac{1}{\omega_{1}},\frac{1}{\omega_{2}}\right\\}\,,$
(36)

where LCM denotes the least common multiple. Hence, we define a new time scale
for the dynamics (33) and (34) using the transformation $\bar{t}=\omega t$,
where

$\displaystyle\omega$ $\displaystyle:=\frac{2\pi}{T}\,.$ (37)

Consequently, the system (33) and (34) can be expressed as

$\displaystyle\frac{d\hat{G}(\bar{t})}{d\bar{t}}$
$\displaystyle=\frac{1}{\omega}\mathcal{H}(\bar{t})K\text{sgn}(\hat{G}(\bar{t}))+\frac{1}{\omega}\frac{d\mathcal{H}(\bar{t})}{d\bar{t}}\tilde{\theta}(\bar{t})+\frac{1}{\omega}\frac{d\delta(\bar{t})}{d\bar{t}}\,,$
(38) $\displaystyle\frac{d\tilde{\theta}(\bar{t})}{d\bar{t}}$
$\displaystyle=\frac{1}{\omega}K\text{sgn}(\hat{G}(\bar{t}))\,.$ (39)

Now, we can implement a suitable averaging mechanism within the transformed
time scale $\bar{t}$ based on the dynamics (38) and (39). Despite the
discontinuity on the right-hand side of equations (38) and (39), the closed-
loop system maintains its periodicity over time due to the periodic probing
and demodulation signals. This unique characteristic allows for the
application of the averaging results established by Plotnikov [24] to this
particular setup.

###  II-E Average Closed-Loop System

Let us define the augmented state as:

$\displaystyle
X^{T}(\bar{t}):=\begin{bmatrix}\hat{G}(\bar{t})\,,\tilde{\theta}(\bar{t})\end{bmatrix}\,,$
(40)

which reduces the system (38) and (39) to:

$\displaystyle\dfrac{dX(\bar{t})}{d\bar{t}}$
$\displaystyle=\dfrac{1}{\omega}\mathcal{F}\left(\bar{t},X,\dfrac{1}{\omega}\right)\,.$
(41)

The system (41) features a small parameter $1/\omega$ and a $T$-periodic
function $\mathcal{F}\left(\bar{t},X,\dfrac{1}{\omega}\right)$ in $\bar{t}$.
Therefore, the averaging theorem [24] can be applied to
$\mathcal{F}\left(\bar{t},X,\dfrac{1}{\omega}\right)$ at
$\displaystyle\lim_{\omega\to\infty}\dfrac{1}{\omega}=0$, following the
principles established by Plotnikov [24].

By employing the averaging of (41), we derive the following average system:

$\displaystyle\dfrac{dX^{\rm{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle=\dfrac{1}{\omega}\mathcal{F}^{\rm{av}}\left(X^{\rm{av}}\right)\,,$
(42) $\displaystyle\mathcal{F}^{\rm{av}}\left(X^{\rm{av}}\right)$
$\displaystyle=\dfrac{1}{T}\int_{0}^{T}\mathcal{F}\left(\xi,X^{\rm{av}},0\right)d\xi\,,$
(43)

where we “freeze” the average states of $\hat{G}(\bar{t})$ and
$\tilde{\theta}(\bar{t})$, resulting in the following equations, for all
$t\geq 0$:

$\displaystyle\frac{d\hat{G}^{\rm{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle=\frac{1}{\omega}HK\text{sgn}(\hat{G}^{\rm{av}}(\bar{t}))\,,$
(44) $\displaystyle\frac{d\tilde{\theta}^{\rm{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle=\frac{1}{\omega}K\text{sgn}(\hat{G}^{\rm{av}}(\bar{t}))\,,$
(45) $\displaystyle\hat{G}^{\rm{av}}(\bar{t})$
$\displaystyle=H\tilde{\theta}^{\rm{av}}(\bar{t})\,.$ (46)

Now, we introduce a necessary and sufficient condition regarding the matrix
$HK$.

###### Assumption 1

The matrix $HK=\begin{bmatrix}H_{11}^{1}K_{1}&H_{12}^{1}K_{2}\\\
H_{21}^{2}K_{1}&H_{22}^{2}K_{2}\end{bmatrix}$ is strictly diagonally dominant,
i.e., for $j\neq i$,

$\displaystyle|H_{ij}^{i}K_{j}|<|H_{ii}^{i}K_{i}|\,,\quad i\in\\{1,2\\}\,.$
(47)

Additionally,

$\displaystyle\text{sgn}(K_{i})=-\text{sgn}(H_{ii}^{i})\,.\quad
i\in\\{1,2\\}\,.$ (48)

According to Assumption 1, we can establish the existence and uniqueness of
the Nash equilibrium. This conclusion stems from the fact that strictly
diagonally dominant matrices, as proven by the Levy–Desplanques theorem [27,
13], are nonsingular.

##  III Stability Analysis

The next theorem guarantees finite-time stability of the proposed NES control
system via sliding modes, as depicted in Fig. 1.

###### Theorem 1

Consider the closed-loop average dynamics of the pseudogradient estimate (44),
and that Assumption 1 holds. For a sufficiently large $\omega>0$, defined in
(37), there exists a finite time
$t_{s}\in\left(0\,,\dfrac{\|\hat{G}^{\text{av}}(0)\|}{\min\\{|H_{11}^{1}K_{1}|-|H_{12}^{1}K_{2}|,|H_{22}^{2}K_{2}|-|H_{21}^{2}K_{1}|\\}}\right]$
such that the sliding motion on $\hat{G}_{\rm{av}}(t)\equiv 0$ occurs.
Furthermore, for the non-average system (34), there exist constants $m$ and
$M_{\theta}$ such that

$\displaystyle\|\theta(t)\mathbb{-}\theta^{\ast}\|\mathbb{\leq}\mathcal{O}\left(a\mathbb{+}\dfrac{1}{\omega}\right)+\begin{cases}M_{\theta}\|\theta(0)\mathbb{-}\theta^{\ast}\|\mathbb{-}mt,\forall
t\in[0,t_{s}),\\\ 0,\forall
t\in[t_{s},+\infty),\end{cases}\\!\\!\\!\\!\\!\\!\\!\\!$ (49)

where $a$ is defined as $a=\sqrt{a_{1}^{2}+a_{2}^{2}}$, with $a_{1}$ and
$a_{2}$ being specified in (9) and (10). The constant $m$ depends on both the
control gain $K$ and the Hessian $H$, while $M_{\theta}$ solely relies on the
Hessian $H$.

Proof: Let us examine the following candidate Lyapunov function for the
average system (44):

$\displaystyle V_{\text{av}}(\bar{t})$
$\displaystyle=[\hat{G}_{1}^{\text{av}}(\bar{t})]^{2}+[\hat{G}_{2}^{\text{av}}(\bar{t})]^{2}\,.$
(50)

The time derivative of (50) is given by

$\displaystyle\frac{dV_{\text{av}}(\bar{t})}{d\bar{t}}=2\hat{G}_{1}^{\text{av}}(\bar{t})\frac{d\hat{G}_{1}^{\text{av}}(\bar{t})}{d\bar{t}}+2\hat{G}_{2}^{\text{av}}(\bar{t})\frac{d\hat{G}_{2}^{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle=+\frac{2}{\omega}\hat{G}_{1}^{\text{av}}(\bar{t})(H_{11}^{1}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))+H_{12}^{1}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t})))+$
$\displaystyle\quad+\frac{2}{\omega}\hat{G}_{2}^{\text{av}}(\bar{t})(H_{21}^{2}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))+H_{22}^{2}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t})))$
$\displaystyle=\frac{2}{\omega}H_{11}^{1}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))\hat{G}_{1}^{\text{av}}(\bar{t})+\frac{2}{\omega}H_{12}^{1}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))\hat{G}_{1}^{\text{av}}(\bar{t})+$
$\displaystyle\quad+\frac{2}{\omega}H_{21}^{2}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))\hat{G}_{2}^{\text{av}}(\bar{t})+\frac{2}{\omega}H_{22}^{2}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))\hat{G}_{2}^{\text{av}}(\bar{t})$
$\displaystyle=\frac{2}{\omega}H_{11}^{1}K_{1}|\hat{G}_{1}^{\text{av}}(\bar{t})|+\frac{2}{\omega}H_{12}^{1}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))\hat{G}_{1}^{\text{av}}(\bar{t})+$
$\displaystyle\quad+\frac{2}{\omega}H_{21}^{2}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))\hat{G}_{2}^{\text{av}}(\bar{t})+\frac{2}{\omega}H_{22}^{2}K_{2}|\hat{G}_{2}^{\text{av}}(\bar{t})|\,.$
(51)

By using (48), equation (51) becomes

$\displaystyle\frac{dV_{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle=-\frac{2}{\omega}|H_{11}^{1}K_{1}||\hat{G}_{1}^{\text{av}}(\bar{t})|+\frac{2}{\omega}H_{12}^{1}K_{2}\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))\hat{G}_{1}^{\text{av}}(\bar{t})+$
$\displaystyle\quad+\frac{2}{\omega}H_{21}^{2}K_{1}\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))\hat{G}_{2}^{\text{av}}(\bar{t})-\frac{2}{\omega}|H_{22}^{2}K_{2}||\hat{G}_{2}^{\text{av}}(\bar{t})|\,,$
(52)

whose upper bound can be expressed as

$\displaystyle\frac{dV_{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle\leq-\frac{2}{\omega}|H_{11}^{1}K_{1}||\hat{G}_{1}^{\text{av}}(\bar{t})|+\frac{2}{\omega}|H_{12}^{1}K_{2}||\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))||\hat{G}_{1}^{\text{av}}(\bar{t})|+$
$\displaystyle\quad+\frac{2}{\omega}|H_{21}^{2}K_{1}||\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))||\hat{G}_{2}^{\text{av}}(\bar{t})|-\frac{2}{\omega}|H_{22}^{2}K_{2}||\hat{G}_{2}^{\text{av}}(\bar{t})|$
$\displaystyle=-\frac{2}{\omega}|H_{11}^{1}K_{1}||\hat{G}_{1}^{\text{av}}(\bar{t})|+\frac{2}{\omega}|H_{12}^{1}K_{2}||\text{sgn}(\hat{G}_{2}^{\text{av}}(\bar{t}))||\hat{G}_{1}^{\text{av}}(\bar{t})|+$
$\displaystyle\quad+\frac{2}{\omega}|H_{21}^{2}K_{1}||\text{sgn}(\hat{G}_{1}^{\text{av}}(\bar{t}))||\hat{G}_{2}^{\text{av}}(\bar{t})|-\frac{2}{\omega}|H_{22}^{2}K_{2}||\hat{G}_{2}^{\text{av}}(\bar{t})|$
$\displaystyle=-\frac{2}{\omega}(|H_{11}^{1}K_{1}|-|H_{12}^{1}K_{2}|)|\hat{G}_{1}^{\text{av}}(\bar{t})|+$
$\displaystyle\quad-\frac{2}{\omega}(|H_{22}^{2}K_{2}|-|H_{21}^{2}K_{1}|)|\hat{G}_{2}^{\text{av}}(\bar{t})|.$
(53)

From the strict diagonal dominance assumption, equation (47), one can
introduce positive constants $\alpha_{1}=|H_{11}^{1}K_{1}|-|H_{12}^{1}K_{2}|$
and $\alpha_{2}=|H_{22}^{2}K_{2}|-|H_{21}^{2}K_{1}|$ such that inequality (53)
is rewritten as

$\displaystyle\frac{dV_{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle\leq-\frac{2}{\omega}\alpha_{1}|\hat{G}_{1}^{\text{av}}(\bar{t})|-\frac{2}{\omega}\alpha_{2}|\hat{G}_{2}^{\text{av}}(\bar{t})|$
$\displaystyle\leq-\frac{2}{\omega}\alpha_{0}(|\hat{G}_{1}^{\text{av}}(\bar{t})|+|\hat{G}_{2}^{\text{av}}(\bar{t})|)\,,$
(54)

where $\alpha_{0}=\min\\{\alpha_{1}\,,\alpha_{2}\\}$. From equation (50) it is
easy to verify that
$\sqrt{V_{\text{av}}(\bar{t})}=\sqrt{[\hat{G}_{1}^{\text{av}}(\bar{t})]^{2}+[\hat{G}_{2}^{\text{av}}(\bar{t})]^{2}}\leq\sqrt{[\hat{G}_{1}^{\text{av}}(\bar{t})]^{2}}+\sqrt{[\hat{G}_{2}^{\text{av}}(\bar{t})]^{2}}=|\hat{G}_{1}^{\text{av}}(\bar{t})|+|\hat{G}_{2}^{\text{av}}(\bar{t})|$.

Thus,
$\sqrt{V_{\text{av}}(\bar{t})}\leq|\hat{G}_{1}^{\text{av}}(\bar{t})|+|\hat{G}_{2}^{\text{av}}(\bar{t})|$
and inequality (54) is upper-bounded by

$\displaystyle\frac{dV_{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle\leq-\frac{2}{\omega}\alpha_{0}\sqrt{V_{\text{av}}(\bar{t})}\,.$
(55)

By defining
$\tilde{V}_{\rm{av}}(\bar{t}):=\sqrt{V_{\text{av}}(\bar{t})}=\sqrt{[\hat{G}_{1}^{\text{av}}(\bar{t})]^{2}+[\hat{G}_{2}^{\text{av}}(\bar{t})]^{2}}=\|\hat{G}^{\text{av}}(\bar{t})\|$,
one can write its time-derivative as
$\dfrac{d\tilde{V}_{\text{av}}(\bar{t})}{d\bar{t}}=\dfrac{1}{2\sqrt{V_{\text{av}}(\bar{t})}}\dfrac{dV_{\text{av}}(\bar{t})}{d\bar{t}}$,
whose upper bound is

$\displaystyle\frac{d\tilde{V}_{\text{av}}(\bar{t})}{d\bar{t}}$
$\displaystyle\leq-\frac{\alpha_{0}}{\omega}\,,\quad\bar{t}\geq 0,$ (56)

by using inequality (55). From the Comparison Theorem [14, Theorem 7], there
exists an upper bound
$\bar{V}_{\text{av}}(\bar{t})\geq\tilde{V}_{\text{av}}(\bar{t})$ satisfying
the differential equation

$\displaystyle\frac{d\bar{V}_{\text{av}}(\bar{t})}{d\bar{t}}=-\frac{\alpha_{0}}{\omega}\,,\quad\bar{V}_{\text{av}}(0)=\tilde{V}_{\text{av}}(0)\geq
0\,,\quad\bar{t}\geq 0.$ (57)

The integration of both sides of the aforementioned differential equation
leads to the solution

$\displaystyle\bar{V}_{\text{av}}(\bar{t})-\bar{V}_{\text{av}}(0)=-\frac{\alpha_{0}}{\omega}\bar{t}\,,\quad\bar{t}\geq
0\,,$ (58)

or, equivalently,

$\displaystyle\bar{V}_{\text{av}}(\bar{t})=\|\hat{G}^{\text{av}}(0)\|-\frac{\alpha_{0}}{\omega}\bar{t}\,,\quad\bar{t}\geq
0\,.$ (59)

Hence, there exists
$\bar{t}_{1}=\frac{\omega}{\alpha_{0}}\|\hat{G}^{\text{av}}(0)\|$ such that
$\bar{V}_{\text{av}}(\bar{t})=0$, for all $t\geq t_{1}$. Now, it can be
concluded that $\exists\bar{t}_{s}\in(0\,,\bar{t}_{1}]$ when the ideal sliding
mode occurs, such that $\bar{V}_{\text{av}}(\bar{t})=0$,
$\forall\bar{t}\geq\bar{t}_{s}$. Consequently,

$\displaystyle\|\hat{G}^{\text{av}}(\bar{t})\|$
$\displaystyle\leq\|\hat{G}^{\text{av}}(0)\|-\dfrac{\alpha_{0}}{\omega}\bar{t}\,,\quad\forall\bar{t}\in[0,\bar{t}_{s})\,,$
(60) $\displaystyle\|\hat{G}^{\text{av}}(\bar{t})\|$
$\displaystyle=0\,,\quad\forall\bar{t}\geq\bar{t}_{s}\,.$ (61)

Since $H$ is invertible and, from (46),
$\tilde{\theta}_{\rm{av}}(\bar{t})=H^{-1}\hat{G}_{\rm{av}}(\bar{t})$, the
following inequality is established
$\|H^{-1}\|^{-1}\|\tilde{\theta}_{\rm{av}}(\bar{t})\|\leq\|H\tilde{\theta}_{\rm{av}}(\bar{t})\|\leq\|H\|\|\tilde{\theta}_{\rm{av}}(\bar{t})\|$
and, therefore, we can state that

$\displaystyle\|\tilde{\theta}^{\text{av}}(\bar{t})\|$
$\displaystyle\leq\|H^{-1}\|\|H\|\|\tilde{\theta}^{\text{av}}(0)\|-\dfrac{\|H^{-1}\|\alpha_{0}}{\omega}\bar{t}\,,\quad\forall\bar{t}\in[0,\bar{t}_{s})\,,$
(62) $\displaystyle\|\tilde{\theta}^{\text{av}}(\bar{t})\|$
$\displaystyle=0\,,\quad\forall\bar{t}\geq\bar{t}_{s}\,.$ (63)

Since (34) has a discontinuous right-hand side, but is also $T$-periodic in
$t$, and noting that the average system with state
$\tilde{\theta}_{\text{av}}(\bar{t})$ is finite-time stable according to (62)
and (63), we can invoke the averaging theorem in [24, Theorem 2] to conclude
that

$\displaystyle\|\tilde{\theta}(t)-\tilde{\theta}_{\text{av}}(t)\|\leq\mathcal{O}\left(\frac{1}{\omega}\right)\,.$
(64)

By applying the triangle inequality [4], we also obtain:

$\displaystyle\|\tilde{\theta}(t)\|\leq\|\tilde{\theta}_{\text{av}}(t)\|+\mathcal{O}\left(\frac{1}{\omega}\right)$
$\displaystyle\leq\begin{cases}\|H^{-1}\|\|H\|\|\tilde{\theta}^{\text{av}}(0)\|\\!-\\!\|H^{-1}\|\alpha_{0}t\\!+\\!\mathcal{O}\left(\dfrac{1}{\omega}\right)\,,\forall
t\in[0,t_{s})\,,\\\ \mathcal{O}\left(\dfrac{1}{\omega}\right)\,,\forall
t\in[t_{s}\,,+\infty)\,.\end{cases}$ (65)

Now, from (9) and Fig. 1, we can verify that

$\displaystyle\theta(t)-\theta^{\ast}=\tilde{\theta}(t)+S(t)\,,$ (66)

where $S(t):=[S_{1}(t)\,,S_{2}(t)]^{T}$, whose the Euclidean norm satisfies

$\displaystyle\|\theta(t)-\theta^{\ast}\|=\|\tilde{\theta}(t)+S(t)\|\leq\|\tilde{\theta}(t)\|+\|S(t)\|$
$\displaystyle\mathbb{\leq}\begin{cases}\|H^{-1}\|\|H\|\|\theta(0)\mathbb{-}\theta^{\ast}\|\mathbb{-}\|H^{-1}\|\alpha_{0}t\mathbb{+}\mathcal{O}\left(a\mathbb{+}\dfrac{1}{\omega}\right)\,,\forall
t\in[0,t_{s})\,,\\\
\mathcal{O}\left(a\mathbb{+}\dfrac{1}{\omega}\right)\,,\forall
t\in[t_{s}\,,+\infty)\,.\end{cases}$ (67)

leading to inequality (49), for appropriate $m$ and $M_{\theta}$. $\square$

##  IV SIMULATION RESULTS

This section presents simulation results in order to illustrate the variable
structure Nash equilibrium seeking based on sliding mode control. The
investigated system captures a noncooperative game involving two firms
operating in a duopoly market framework. These firms engage in competition
aimed at maximizing their profits $J_{i}(t)$ through the pricing strategy
$u_{i}(t)$ of their respective products, without sharing any information
between the players.

The plant parameters are consistent with those given in [9]: initial
conditions are set as $\hat{\theta}_{1}(0)=50$ and
$\hat{\theta}_{2}(0)=110/3$. Additionally, $S_{d}=100$, $p=0.20$, $m_{1}=30$
and $m_{2}=30$. According to (7), the parameter values in [9] for the duopoly
game lead to a unique Nash equilibrium given by:

$\displaystyle\theta^{\ast}$ $\displaystyle=[43.3333\,,~{}36.6667]^{T}\,,$
(68) $\displaystyle J^{\ast}$ $\displaystyle=[888.8889\,,~{}222.2222]^{T}\,.$
(69)

The controller parameters are chosen as follows: for the NES strategy, we set
$a_{1}=0.075$, $a_{2}=0.050$, $K_{1}=0.6$, $K_{2}=0.6$, $\omega_{1}=27$, and
$\omega_{2}=22$.

(a) Discontinuous control inputs.  (b) Pseudogradient estimates.

(c) Input signals of payoff functions, $\theta(t)$.  (d) Payoff functions.

Figure 2: Numerical simulations for NES in the noncooperative duopoly game
through distributed sliding mode control policies.

To achieve the Nash equilibrium (68) with the corresponding payoffs (69),
players P1 and P2 implement the proposed decentralized sliding mode NES
strategies employing sinusoidal perturbations to execute their optimal
actions.

The time evolution of the proposed approach is depicted in Fig. 2. Fig. 2(a)
shows the discontinuous behavior of the players’ actions, consistent with the
fact that each player estimates a distinct pseudogradient component and
decides independently to switch its action. In a simulation spanning 150
seconds, the pseudogradient estimate is quickly brought to a neighborhood near
the origin since the sliding surfaces are specified by
$\hat{G}_{1}^{\text{av}}=\hat{G}_{2}^{\text{av}}=0$; see Fig. 2(b).

Although the players employ distributed sliding mode strategies independently,
their collective effect drives the system toward the Nash equilibrium, as
depicted in Fig. 2(c). Additionally, Fig. 2(d) illustrates how, within the
described duopoly market structure and without information sharing, each
player can optimize its corresponding payoff $J_{i}(t)$ by adopting the
proposed variable structure Nash equilibrium seeking via sliding mode to
determine the product prices $\theta_{1}(t)$ and $\theta_{2}(t)$ in the non-
cooperative game.

##  V CONCLUSIONS

This paper has presented a novel approach to attain locally stable convergence
to Nash equilibrium in noncooperative games via a distributed sliding mode
control scheme employing extremum seeking based on sinusoidal periodic
perturbations. This methodology empowers players in a duopoly game to enhance
their payoffs without necessitating information sharing, as action updates
occur autonomously. Each player estimates and updates its sequence of actions
based on distinct pseudogradient components by using relay functions. The
proposed approach not only advances game theory but also provides a practical
means for decentralized decision-making in intricate systems. Finite-time
stability for the average closed-loop system is proved rather than exponential
stability, generally arrived at classical extremum seeking with proportional
control laws.

The results presented should be regarded as an SMC alternative to the
continuous but non-smooth Nash seeking algorithms with fixed-time convergence
properties in [25]. In the current paper, we have used
discontinuous/signum/relay functions of the estimates of the pseudogradient,
conducting averaging analysis of the discontinuous system with an averaging
theorem of Plotnikov for differential inclusions [24], and achieve finite-time
convergence (with a settling time dependent on the initial condition).
Specifically, the fixed-time result of [25] may look stronger since it is of
fixed-time nature (the settling time is independent of the initial condition),
but this conference paper provides an alternative approach, that of sliding
mode, to Nash equilibrium seeking.

## References

  * [1] S. Amina, G. A. Schwartz, and S. S. Sastry, “Security of interdependent and identical networked control systems,” _Automatica_ , vol. 49, pp. 186–192, 2013.
  * [2] N. O. Aminde, T. R. Oliveira, and L. Hsu, “Global output-feedback extremum seeking control via monitoring functions,” in _52nd IEEE Conference on Decision and Control (CDC)_ , 2013, pp. 1031–1036.
  * [3] ——, “Multivariable extremum seeking control via cyclic search and monitoring function,” _International Journal of Adaptive Control and Signal Processing_ , vol. 35, pp. 1217–1232, 2021.
  * [4] T. Apostol, _Mathematical Analysis - A Modern Approach to Advanced Calculus_. Addison-Wesley Publishing Company, 1957.
  * [5] T. Başar, “Relaxation techniques and asynchronous algorithms for on-line computation of noncooperative equilibria,” _J. Economic Dynamics and Control_ , vol. 11, pp. 531–549, 1987.
  * [6] T. Başar and G. J. Olsder, _Dynamic Noncooperative Game Theory_. SIAM Series in Classics in Applied Mathematics, 1999.
  * [7] T. Başar and G. Zaccour, Eds., _Handbook of Dynamic Game Theory, Volume I (Theory of Dynamic Games)_. Springer International Publishing, 2018.
  * [8] ——, _Handbook of Dynamic Game Theory, Volume II (Applications of Dynamic Games)_. Springer International Publishing, 2018.
  * [9] P. Frihauf, M. Krstić, and T. Başar, “Nash equilibrium seeking in noncooperative games,” _IEEE Trans. Automat. Contr._ , vol. 57, pp. 1192–1207, 2012.
  * [10] D. Fudenberg and J. Tirole, _Game Theory_. The MIT Press, 1991.
  * [11] A. Ghaffari, M. Krstić, and D. Nes̆ic, “Multivariable Newton-based extremum seeking,” _Automatica_ , vol. 48, pp. 1759–1767, 2012.
  * [12] Z. Han, D. Niyato, W. Saad, and T. Başar, _Game Theory for Next Generation Wireless and Communication Networks: Modeling, Analysis, and Design_. Cambridge University Press, 2019\.
  * [13] R. A. Horn and C. R. Johnson, _Matrix Analysis_. Cambridge, U.K.: Cambridge Univ. Press, 1985.
  * [14] H. K. Khalil, _Nonlinear Systems_. Prentice Hall, 2002.
  * [15] M. Krstić, “Extremum seeking control,” in _Encyclopedia of Systems and Control_. London: Springer, 2014, vol. 1, pp. 413–416.
  * [16] M. Krstić and H.-H. Wang, “Stability of extremum seeking feedback for general nonlinear dynamic systems,” _Automatica_ , vol. 36, pp. 595–601, 2000.
  * [17] S. Li and T. Başar, “Distributed learning algorithms for the computation of noncooperative equilibria,” _Automatica_ , vol. 23, pp. 523–533, 1987\.
  * [18] J. A. Moreno, “Lyapunov approach for analysis and design of second order sliding mode algorithms,” in _Sliding Modes after the First Decade of the 21st Century_ , ser. Lecture Notes in Control and Information Sciences, L. Fridman, J. Moreno, and R. Iriarte, Eds. Berlin, Heidelberg: Springer, 2011, vol. 412, pp. 113–149. [Online]. Available: https://doi.org/10.1007/978-3-642-22164-4_4
  * [19] J. A. Moreno and M. Osorio, “Strict lyapunov functions for the super-twisting algorithm,” _IEEE Transactions on Automatic Control_ , vol. 57, no. 4, pp. 1035–1040, 2012.
  * [20] J. F. Nash, “Noncooperative games,” _Annals of Mathematics_ , vol. 54, pp. 286–295, 1951.
  * [21] T. R. Oliveira, L. Hsu, and A. J. Peixoto, “Output-feedback global tracking for unknown control direction plants with application to extremum-seeking control,” _Automatica_ , vol. 47, pp. 2029–2038, 2011.
  * [22] T. R. Oliveira, A. J. Peixoto, and L. Hsu, “Global real-time optimization by output-feedback extremum-seeking control with sliding modes,” _Journal of the Franklin Institute_ , vol. 349, pp. 1397–1415, 2012.
  * [23] Y. Pan, Ü. Özgüner, and T. Acarman, “Stability and performance improvement of extremum seeking control with sliding mode,” _International Journal of Control_ , vol. 76, pp. 968–985, 2003.
  * [24] V. A. Plotnikov, “Averaging of differential inclusions,” _Ukrainian Mathematical Journal_ , vol. 31, pp. 454–457, 1980.
  * [25] J. I. Poveda, M. Krstić, and T. Başar, “Fixed-time Nash equilibrium seeking in time-varying networks,” _IEEE Trans. Automat. Contr._ , vol. 68, pp. 1954–1969, 2023.
  * [26] J. B. Rosen, “Existence and Uniqueness of Equilibrium Points for Concave N-Person Games,” _Econometrica_ , vol. 33, no. 3, pp. 520–534, 1965.
  * [27] O. Taussky, “A recurring theorem on determinants,” _Amer. Math. Monthly_ , vol. 56, no. 10, pp. 672–676, 1949.
  * [28] V. Todorovski and M. Krstić, “Practical prescribed-time seeking of a repulsive source by unicycle angular velocity tuning,” _Automatica_ , vol. 68, p. 111069, 2023.
  * [29] C. T. Yilmaz and M. Krstić, “Prescribed-time extremum seeking with chirpy probing for PDEs–Part I: Delay,” _IEEE American Control Conference_ , pp. 1000–1005, 2022.
  * [30] Q. Zhu, H. Tembine, and T. Başar, “Hybrid learning in stochastic games and its applications in network security,” in _Series on Computational Intelligence_ , F. L. Lewis and D. Liu, Eds. IEEE Press/Wiley, 2013, ch. 14, pp. 305–329.