On periodic solutions and attractors

for the Maxwell–Bloch equations

A. I. Komech

Faculty of Mathematics, University of Vienna

<EMAIL_ADDRESS>
###### Abstract

We consider the Maxwell–Bloch system which is a finite-dimensional
approximation of the coupled nonlinear Maxwell–Schrödinger equations. The
approximation consists of one-mode Maxwell field coupled to two-level
molecule. We construct time-periodic solutions to the factordynamics which is
due to the symmetry gauge group. For the corresponding solutions to the
Maxwell–Bloch system, the Maxwell field, current and the population inversion
are time-periodic, while the wave function acquires a unit factor in the
period. The proofs rely on high-amplitude asymptotics of the Maxwell field and
a development of suitable methods of differential topology: the transversality
and orientation arguments. We also prove the existence of the global compact
attractor.

MSC classification: 78A60, 37N20, 81V80, 70K28, 35P30.

Keywords: Maxwell–Bloch equations; laser; current; population inversion;
periodic solutions; Rabi solution; gauge group; Hopf fibration; Duhamel
representation; differential topology; homotopy; intersection theory;
transversality; orientation; attractor.

###### Contents

  1. 1 Introduction
  2. 2 The Maxwell–Bloch equations with one molecule
  3. 3 The high-amplitude asymptotics of the Maxwell field
  4. 4 Gauge symmetry and factordynamics in the Hopf fibration
  5. 5 Periodic solutions and fixed points of the Poincaré map
  6. 6 A priori bound for fixed points
  7. 7 The transversality and orientation

##  1 Introduction

The paper concerns an old problem of the quantum optics [28] on existence of
time-periodic solutions to the nonlinear coupled Maxwell–Bloch equations in
the case of time-periodic pumping. The Maxwell–Bloch equations are
traditionally used for the semiclassical description of the laser action [10,
13, 24, 26, 27]. The equations are finite-dimensional approximation of the
semiclassical Maxwell–Schrödinger system in a bounded cavity, see [1, 9, 14,
16].

For the two-level Schrödinger equation with a given harmonic Maxwell field,
time-periodic solutions were constructed by Rabi in the rotating wave
approximation [10, 21, 22, 28]. For various versions of the nonlinear
Maxwell–Bloch equations, the problems of existence of time-periodic solutions
and their stability are discussed since 1960’. For the model [20] with a time-
periodic source, the approximate Nth-order time-periodic solutions were
recently constructed in [30] by perturbation technique. For the
phenomenological model [2], time-periodic solutions were recently constructed
in [6] in the absence of time-periodic source [6, (1)]. The solutions are
constructed for small interaction constants relying on homotopy invariance of
the degree [4] and averaging arguments [5].

We obtain the results for the Maxwell–Bloch equations which are finite-
dimensional Galerkin approximation of the semiclassical Maxwell–Schrödinger
system in a bounded cavity, see [1, 9, 14, 16]. We prove the existence of
time-periodic solutions for the equations reduced by the symmetry gauge group.
The existence is proved for the corresponding factordynamics in the Hopf
fibration. For the corresponding solutions to the Maxwell–Bloch equations, the
Maxwell field, the current and the population inversion are time-periodic,
while the wave function acquires a unit factor in the period, see (5.6).

The time-periodic solutions correspond to fixed points of the corresponding
Poincaré map. Our construction of the fixed points relies on a detailed high-
amplitude asymptotics of the Maxwell field and suitable development of methods
of differential topology: the theory of transversality and homotopy invariance
of the orientation [8, 11].

The main issue in our development of the homotopy and transversality methods
of differential topology is that the homotopy invariance of the intersection
number is known for compact spaces, while in our case the space is not
compact. The intersection number can be changed in the homotopy if some of
intersection points run to infinity in a finite time. The key role in our
proofs is played by “a priori bound” for fixed points which prevents such
runaway and reduces the problem to the compact space. The bound follows from
the high-amplitude asymptotics of the Maxwell field which is established by a
novel careful analysis of the integral Duhamel representation. As a byproduct,
we obtain for the Maxwell–Bloch equations the global well-posedness and the
existence of a compact global attractor [7, 29].

Further we streamline the ideas of the theory of transversality and
orientation which lie behind the homotopy theory of intersection of smooth
manifolds [8, 11].
Acknowledgements. The author thanks V.A. Vassiliev who suggested the relation
to the intersection theory.

##  2 The Maxwell–Bloch equations with one molecule

We use the Heaviside–Lorentz units [12] (unrationalized Gaussian units), in
which the main physical constants (electron charge and mass, Planck’s
constant, and the speed of light in vacuum) read [31, p. 221]

$e\\!=\\!-4.8\\!\times\\!10^{-10}{\rm esu},\,{\rm
m}\\!=\\!9.1\\!\times\\!10^{-28}{\rm
g},\,\hbar\\!=\\!1.055\\!\times\\!10^{-27}\mbox{\rm erg$\cdot$
s},\,c\\!=\\!3.0\\!\times\\!10^{10}{\rm cm/s}.$

The Maxwell–Bloch equations describe the coupling of one-mode Maxwell field
with a two-level molecule in a bounded cavity $V\subset{\mathbb R}^{3}$:

${\bf A}(x,t)=A(t){\bf
X}(x),\quad\psi(x,t)=C_{1}(t)\varphi_{1}(x)+C_{2}(t)\varphi_{2}(x),\qquad x\in
V.$ (2.1)

Here ${\bf X}(x)$ and $\varphi_{l}$ are suitable normalized eigenfunctions of
the Laplace and Schrödinger operators:

$\left\\{\\!\\!\\!\\!\begin{array}[]{rl}\Delta{\bf
X}(x)\\!=\\!-\\!\frac{\Omega^{2}}{c^{2}}{\bf X}(x),{\rm div{\hskip
1.42262pt}}{\bf
X}(x)\\!=\\!0,\,x\\!\in\\!V;\\!\\!\\!\\!\\!&\\!\\!{\bf\nabla}{\bf
X}(x)\\!\wedge\\!{\bf X}(x)\\!=\\!0,\,{\bf n}(x)\\!\cdot\\!{\rm curl{\hskip
1.42262pt}}{\hskip 1.42262pt}{\bf X}(x)\\!=\\!0,\,x\\!\in\\!\partial V\\\ \\\
~{}\displaystyle{\mathbb H}\varphi_{l}(x)=\omega_{l}\varphi_{l}(x),\,\,\,x\in
V;\\!\\!&\\!\\!\varphi_{l}(x)=0,\,\,x\in\partial V\end{array}\right|.$ (2.2)

Here $\Omega$ is the frequency of the Maxwell mode, and ${\mathbb
H}:=-\frac{\hbar^{2}}{2{\rm m}}\Delta+e\Phi(x)$, where $\Phi(x)$ is the
molecular (ion’s) potential: we suppose that only one molecular electron is
involved. The semiclassical Maxwell–Bloch equations must be defined as the
Galerkin approximation of the damped-driven Maxwell–Schrödinger equations (see
[3, 9, 23] and also [15, 17, 18, 19] for the equations without damping and
pumping) :

$\frac{1}{c^{2}}\dot{A}(t)=H_{B},\quad\frac{1}{c^{2}}\dot{B}(t)=-H_{A}-\frac{\sigma}{c^{2}}B;\qquad
i\hbar\dot{C}_{l}(t)=H_{\overline{C}_{l}},\quad l=1,2.$ (2.3)

Here $\sigma>0$ is the electrical conductivity of the cavity medium, and the
Hamiltonian $H$ is

$H=H(A,B,C,t)={\cal H}(A{\bf X},B{\bf
X},C_{1}\varphi_{1}+C_{2}\varphi_{2},t),\qquad C=(C_{1},C_{2}),$ (2.4)

where ${\cal H}$ is the Hamiltonian of the Maxwell–Schrödinger equations with
a pumping. Neglecting the spin and scalar potential which can be easily added,
the Hamiltonian ${\cal H}$ reads as

${\cal H}({\bf A},{\bf B},\psi,t)=\frac{1}{2}[\|\frac{1}{c}{\bf
B}\|^{2}+\|{\rm curl{\hskip 1.42262pt}}{\bf A}\|^{2}]+\langle\psi,H({\bf
A},t)\psi\rangle,$ (2.5)

where $\|\cdot\|$ stands for the norm in the phase Hilbert space
$L^{2}({\mathbb R}^{3})\otimes{\mathbb R}^{3}$, and the brackets
$\langle\cdot,\cdot\rangle$ stand for the Hermitian inner product in
$L^{2}({\mathbb R}^{3})\otimes{\mathbb C}$. The Schrödinger operator reads as

$\displaystyle H({\bf A},t)$ $\displaystyle:=$
$\displaystyle\frac{1}{2m}[-i\hbar\nabla-\displaystyle\frac{e}{c}({\bf
A}(x)+{\bf A}_{p}(x,t))]^{2}+e\Phi(x)$ (2.6) $\displaystyle=$
$\displaystyle\displaystyle{\mathbb H}+i\frac{e\hbar}{mc}({\bf A}(x)+{\bf
A}_{p}(x,t))\nabla+\frac{e^{2}}{2mc^{2}}({\bf A}(x)+{\bf A}_{p}(x,t))^{2}$

where ${\bf A}_{p}(x,t)={\bf X}(x)A_{p}(t)$ is the pumping. Substituting (2.1)
into (2.5), we find

$H(A,B,C,t)=\frac{1}{2c^{2}}[B^{2}+\Omega^{2}A^{2}]+\langle\psi,H({\bf
A},t)\psi\rangle.$ (2.7)

Using (2.6), we obtain that

$\displaystyle\langle\psi,H({\bf
A},t)\psi\rangle=\hbar\omega_{1}|C_{1}|^{2}+\hbar\omega_{2}|C_{2}|^{2}$
$\displaystyle+i\frac{e\hbar}{mc}\sum_{l,l^{\prime}}\overline{C}_{l}C_{l^{\prime}}(A+A_{p}(t))\langle\varphi_{l}(x){\bf
X}(x),\nabla\varphi_{l^{\prime}}(x)\rangle+\frac{e^{2}}{2mc^{2}}(A(x)+A_{p}(x,t))^{2}{\bf
X}^{2}(x).\qquad\qquad\qquad\qquad\qquad$

Substituting into (2.7), we get

$\displaystyle
H(A,B,C,t)=\frac{1}{2c^{2}}[B^{2}+\Omega^{2}A^{2}]+\hbar\omega_{1}|C_{1}|^{2}+\hbar\omega_{2}|C_{2}|^{2}$
$\displaystyle+i\frac{\hbar}{mc}\sum_{l,l^{\prime}}\overline{C}_{l}C_{l^{\prime}}Q_{l,l^{\prime}}(A+A_{p}(t))+\frac{e^{2}}{2mc^{2}}(A(x)+A_{p}(x,t))^{2}{\bf
X}^{2}(x).\qquad\qquad$ (2.8)

where

$Q_{l,l^{\prime}}=e\langle{\bf
X}(x)\varphi_{l}(x),\nabla\varphi_{l^{\prime}}(x)\rangle.$ (2.9)

The last term on the right hand side of (2) is negligible with respect to the
first term because usually $\Omega^{2}\gg\frac{e^{2}}{m}$. For example,
$\Omega\approx 3\times 10^{15}s^{-1}$ for the Ruby laser [13, 27], while
$\frac{e^{2}}{m}\approx 3\times 10^{8}$. This is why the last term is
traditionally neglected [25, (44.13)], and we also will neglect this term in
the Hamiltonian. Moreover, we will use the standard dipol approximation which
physically means that the wavelength $\lambda=2\pi c/\Omega$ is negligible
with respect to the molecular size. In this case,

$Q_{l,l^{\prime}}\approx Q_{l,l^{\prime}}^{d}=e{\bf
X}(x_{*})\langle\varphi_{l}(x),\nabla\varphi_{l^{\prime}}(x)\rangle,$ (2.10)

where $x_{*}\in V$ is the location of the molecule. As the result, we accept
the Hamiltonian

$H(A,B,C,t)\\!=\\!\frac{1}{2c^{2}}[B^{2}\\!+\\!\Omega^{2}A^{2}]\\!+\\!\hbar\omega_{1}|C_{1}|^{2}\\!+\\!\hbar\omega_{2}|C_{2}|^{2}\\!+\\!i\frac{e\hbar}{mc}\sum_{l,l^{\prime}}Q_{l,l^{\prime}}^{d}\overline{C}_{l}C_{l^{\prime}}(A\\!+\\!A_{p}(t)).$
(2.11)

######  Remark 2.1.

We make the standard approximations to keep the succession of our results with
this commonly recognized practice.

The commutation $[H_{\Phi},x]=-\frac{\hbar^{2}}{m}\nabla$ implies well known
identity [25, (44.20)]:

$e\langle\varphi_{l},\nabla\varphi_{l^{\prime}}\rangle\\!=\\!-\frac{em}{\hbar^{2}}\langle\varphi_{l},[H,x]\varphi_{l^{\prime}}\rangle\\!=\\!-\\!\frac{em}{\hbar}[\omega_{l}-\omega_{l^{\prime}}]\langle\varphi_{l},x\varphi_{l^{\prime}}\rangle\\!=\\!-\\!\frac{m}{\hbar}\Delta_{ll^{\prime}}{\bf
P}^{ll^{\prime}}.$ (2.12)

Here $\Delta_{ll^{\prime}}\\!:=\\!\omega_{l}\\!-\\!\omega_{l^{\prime}}$ and
${\bf P}^{12}={\bf P}^{21}={\bf
P}=e\langle\varphi_{l},x\varphi_{l^{\prime}}\rangle\in{\mathbb R}^{3}$ since
we can assume both wave functions $\varphi_{l}$ to be real, where ${\bf P}$ is
the dipole moment (or polarization) of the molecule (this explains the term
“dipole approximation”). Hence, (2.11) and (2.10) give the approximation

$H(A,B,C,t)=\frac{1}{2c^{2}}[B^{2}+\Omega^{2}A^{2}]+\hbar\omega_{1}|C_{1}|^{2}+\hbar\omega_{2}|C_{2}|^{2}-\frac{2\Delta}{c}Q[A+A_{p}(t)]\,{\rm
Im{\hskip 1.42262pt}}[\overline{C}_{1}C_{2}],$ (2.13)

where $Q={\bf X}(x_{*}){\bf P}$ and $\Delta=\Delta_{21}$. Now the Hamilton
equations (2.3) read as

$\\!\\!\\!\left\\{\\!\\!\begin{array}[]{rcl}\dot{A}(t)\\!\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!\\!B(t),\quad\dot{B}(t)=-\Omega^{2}A(t)\\!-\\!\sigma
B(t)+cj(t)\\\ \\\
i\hbar\dot{C}_{1}(t)\\!\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!\\!\hbar\omega_{1}C_{1}(t)+iW(t)\,C_{2}(t),\,\,i\hbar\dot{C}_{2}(t)=\hbar\omega_{2}C_{2}(t)-iW(t)\,C_{1}(t)\end{array}\right|,$
(2.14)

where the current $j(t)$ and the function $W(t)$ are given by

$j(t)={2}q\,{\rm Im{\hskip 1.42262pt}}[\overline{C}_{1}(t)C_{2}(t)],\qquad
q={\Delta Q};\qquad W(t):=\frac{q}{c}[A(t)+A_{p}(t)].$ (2.15)

The charge conservation holds:

$|C_{1}(\tau)|^{2}+|C_{2}(\tau)|^{2}=1,\qquad t\in{\mathbb R}.$ (2.16)

It follows by differentiation from the last two equations of (2.14) since the
function $W(t)$ is real-valued.

To simplify further calculations, we replace the variable $A$ by ${\bf
A}=\Omega A$. Writing again $A$ instead of ${\bf A}$ and $A_{p}/\Omega$
instead of $A_{p}$, we obtain the system (2.14) in the form

$\\!\\!\left\\{\\!\\!\\!\begin{array}[]{rcl}\dot{A}(t)\\!\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!\\!\Omega
B(t),\quad\dot{B}(t)=-\Omega A(t)\\!-\\!\sigma B(t)+cj(t)\\\ \\\
i\hbar\dot{C}_{1}(t)\\!\\!\\!&\\!\\!\\!=\\!\\!\\!&\\!\\!\\!\hbar\omega_{1}C_{1}(t)\\!+\\!iW(t)\,C_{2}(t),\,\,i\hbar\dot{C}_{2}(t)\\!=\\!\hbar\omega_{2}C_{2}(t)\\!-\\!iW(t)\,C_{1}(t)\end{array}\right|,$
(2.17)

where the current $j(t)$ and the function $W(t)$ are given by

$j(t)=\varkappa\,{\rm Im{\hskip
1.42262pt}}[\overline{C}_{1}(t)C_{2}(t)],\qquad\varkappa=2q;\qquad
W(t):=\frac{q}{c\Omega}[A(t)+A_{p}(t)].$ (2.18)

##  3 The high-amplitude asymptotics of the Maxwell field

We prove the existence and uniqueness of global solutions to the system (2.17)
for every initial date, and the existence of bounded global attractor. It
suffices to prove suitable a priori bound for solutions
$X(t)=(A(t),B(t),C(t))$. Denote the phase space ${\mathbb X}={\mathbb
R}^{2}\times{\mathbb C}^{2}$.

######  Lemma 3.1.

The system (2.17) admits a unique solution $X(t)\in C^{1}([0,\infty),{\mathbb
X})$ for any initial state $X(0)\in{\mathbb X}$. The apriori bound holds,

$|X(t)|\leq C|X(0)|e^{-\gamma t}+D,\qquad t>0,$ (3.1)

where $\gamma>0$ and $C,D>0$ do not depend on initial state.

###### Proof.

The bound for $C_{l}(t)$ follows from the charge conservation (2.16). It
remains to prove the bound for the Maxwell amplitudes:

$|A(t)|+|B(t)|\leq C(X(0))e^{-\gamma t}+D,\qquad t>0.$ (3.2)

We rewrite the first line of (2.17) as

$\dot{X}^{1}(t)=MX^{1}(t)+J(t),\qquad t>0,\qquad
X^{1}(t)=\left(\\!\\!\\!\begin{array}[]{c}A(t)\\\
B(t)\end{array}\\!\\!\\!\right),$ (3.3)

where

$M=\left(\\!\\!\begin{array}[]{cc}0&\Omega\\\
-\Omega&-\sigma\end{array}\\!\\!\right),\qquad
J(t)=\left(\\!\\!\\!\begin{array}[]{c}0\\\
c^{2}j(t)\end{array}\\!\\!\\!\right).$ (3.4)

The Duhamel representation gives

$X^{1}(t)=e^{Mt}X^{1}(0)+I(t),\qquad I(t)=\int_{0}^{t}e^{M(t-s)}J(s)ds,\qquad
t>0.$ (3.5)

We can assume that $\sigma<2\Omega$, so the eigenvalues of the matrix $M$ are
$\lambda_{\pm}=-\gamma\pm i\tilde{\Omega}$, where

$\gamma=\sigma/2,\qquad\tilde{\Omega}=\sqrt{\Omega^{2}-\gamma^{2}}.$ (3.6)

It is easy to calculate the exponential of the matrix ${Mt}$ using the basis
of its complex eigenvectors:

$e^{Mt}=e^{\lambda_{+}t}e_{+}\otimes e_{+}^{*}+e^{\lambda_{-}t}e_{-}\otimes
e_{-}^{*},$ (3.7)

where $e_{\pm}\in{\mathbb C}^{2}$ are complex eigenvectors of the matrix $M$,
and $e_{\pm}^{*}$ is the dual basis to $e_{\pm}$. Hence,

$\displaystyle e^{\gamma t}e^{Mt}=e^{i\tilde{\Omega}t}e_{+}\otimes
e_{+}^{*}+e^{-i\tilde{\Omega}t}e_{-}\otimes e_{-}^{*}$ (3.8)

The eigenvalues $\lambda_{\pm}$ and eigenvectors $e_{\pm}$ are complex
conjugate: $e_{\pm}=e^{+}\pm ie^{-}$ where $e^{+}=(e_{+}+e_{-})/2$ and
$e^{-}=(e_{+}-e_{-})/(2i)$. Hence,

$\displaystyle e^{\gamma t}e^{Mt}e^{+}$ $\displaystyle=$ $\displaystyle
e^{i\tilde{\Omega}t}e_{+}/2+e^{-i\tilde{\Omega}t}e_{-}/2=e^{i\tilde{\Omega}t}(e^{+}+ie^{-})/2+e^{-i\tilde{\Omega}t}(e^{+}-ie^{-})/2$
(3.9) $\displaystyle=$
$\displaystyle\cos\tilde{\Omega}te^{+}-\sin\tilde{\Omega}te^{-}$

and similarly,

$\displaystyle e^{\gamma t}e^{Mt}e^{-}$ $\displaystyle=$ $\displaystyle
e^{i\tilde{\Omega}t}e_{+}/2i-e^{-i\tilde{\Omega}t}e_{-}/2i=e^{i\tilde{\Omega}t}(e^{+}+ie^{-})/2i-e^{-i\tilde{\Omega}t}(e^{+}-ie^{-})/2i$
(3.10) $\displaystyle=$
$\displaystyle\sin\tilde{\Omega}te^{+}+\cos\tilde{\Omega}te^{-}$

Define new coordinates $Z^{\pm}$ on the plane ${\mathbb R}^{2}$ corresponding
to the real basis $e^{\pm}$:

$X^{1}=Z^{+}e^{+}+Z^{-}e^{-},\qquad Y^{1}\in{\mathbb R}^{2}.$ (3.11)

Then (3.9) and (3.10) imply that

$e^{Mt}Z=e^{-\gamma t}R(t)Z,\qquad
R(t)=\left(\\!\\!\\!\begin{array}[]{rl}\cos\tilde{\Omega}t&\sin\tilde{\Omega}t\\\
-\sin\tilde{\Omega}t&\cos\tilde{\Omega}t\end{array}\\!\\!\\!\right),\qquad
Z=\left(\\!\\!\begin{array}[]{c}Z^{+}\\\ Z^{-}\end{array}\\!\\!\right).$
(3.12)

So, in the “$Z$-representation” (3.11) the solution (3.5) reads

$Z(t)\\!=\\!\left(\\!\\!\begin{array}[]{c}Z^{+}(t)\\\
Z^{-}(t)\end{array}\\!\\!\right)=e^{-\gamma
t}R(t)\left(\\!\\!\begin{array}[]{c}Z^{+}(0)\\\
Z^{-}(0)\end{array}\\!\\!\right)+\tilde{I}(t)$ (3.13)

Here $\tilde{I}(t)$ is the vector $I(t)$ in the $Z$-representation (3.11), and

$|\tilde{I}(t)|\leq\alpha c^{2}\varkappa\int_{0}^{t}e^{-\gamma t}=\frac{\alpha
c^{2}\varkappa}{\gamma}(1-e^{-\gamma t}),$ (3.14)

where $\alpha$ is the distortion coefficient:
$|\tilde{I}(t)|\leq\alpha|I(t)|$. Now (3.2) is proved with

$C\sim\varkappa,\qquad D\sim\frac{\alpha c^{2}\varkappa}{\gamma}.$ (3.15)

since $|X^{1}(t)|\sim|Z(t)|$. ∎

The methods [7, 29] allows us to deduce from the bound (3.1) the following
corollary.

######  Corollary 3.2.

There exists the compact global attractor ${\cal A}\subset{\mathbb X}$ for the
Maxwell–Bloch equations which is the minimal set such that

$X(t)\to{\cal A},\qquad t\to\infty$

for any solution $X(t)$ to (2.17).

The following high-amplitude asymptotics plays the crucial role in the sequel.

######  Corollary 3.3.

(3.13) implies that for bounded $t\in[0,T]$, we have in the polar coordinates
$Z^{+}(t)=r(t)\cos\phi(t)$, $Z^{-}(t)=r(t)\sin\phi(t)$,

$r(t)\\!=\\!e^{-\gamma
t}r(0)\\!+\\!O(1),\,\phi(t)\\!=\\!\phi(0)\\!+\\!\tilde{\Omega}t\\!+\\!O(t/r(0))\\!\\!\\!\\!\\!\\!\\!\mod
2\pi,\,\,r(0)\\!\to\\!\infty.$ (3.16)

##  4 Gauge symmetry and factordynamics in the Hopf fibration

The equations (2.17) are invariant with respect to the group $G$ of gauge
transformations

$g(\theta)(A,B,C)=(A,B,e^{i\theta}C),\qquad\theta\in{\mathbb R}.$ (4.1)

The charge conservation (2.16) means that $C(t)\in S^{3}$, where $S^{3}$ is
the 3D unit sphere in ${\mathbb C}^{2}$. The action (4.1) commutes with the
dynamics (2.17), hence the latter induces the corresponding factordynamics on
the factorspace ${\mathbb Y}={\mathbb X}/G$ which is diffeomorphic to
${\mathbb R}^{2}\times S^{2}$, since $S^{2}$ is the base of the Hopf fibration
$h:S^{3}\to S^{2}$, see [14]. The action (4.1) does not change the Maxwell
field. Let us write the factordynamics as

$\dot{Y}(t)=F(Y(t),t),\qquad t\in{\mathbb R};\qquad
Y(t)=(A(t),B(t),C_{*}(t))\in{\mathbb Y},$ (4.2)

where $C_{*}(t)=hC(t)\in S^{2}$. For $X=(A,B,C)\in{\mathbb X}$, the current
$j(X)=q{\rm Im{\hskip 1.42262pt}}[\overline{C}_{1}C_{2}]$ and the population
inversion $I:=|C_{2}|^{2}-|C_{1}|^{2}$ are invariant with respect to the gauge
transformations (4.1). So, $j_{*}(hC):=j(C)$ is a bounded smooth function on
$S^{2}$:

$|j_{*}(C_{*})|\leq\varkappa,\qquad C_{*}\in S^{2}.$ (4.3)

In particular, the first line of (2.17) can be written as

$\dot{A}(t)=\Omega B(t),\quad\dot{B}(t)=-\Omega A(t)\\!-\\!\sigma
B(t)+c^{2}j_{*}(C_{*}(t)).$ (4.4)

##  5 Periodic solutions and fixed points of the Poincaré map

In this section, we consider the factordynamics (4.2) with a time-periodic
pumping:

$A_{p}(x,t+T)=A_{p}(x,t),\quad t\in{\mathbb R};\qquad T=2\pi/\Omega_{p}.$
(5.1)

In this case, (4.2) is a time-periodic system:

$F(Y,t+T)=F(Y,t),\qquad Y\in{\mathbb Y},\qquad t\in{\mathbb R}.$ (5.2)

The main result of present paper is the following theorem. We will assume that

$0<\Omega^{2}-\gamma^{2}<\Omega_{p}.$ (5.3)

In particular, this condition holds for the pumping frequency close to
resonance, $\Omega_{p}\approx\Omega$.

######  Theorem 5.1.

Let (5.1) and (5.3) hold. Then the factordynamics (4.2) admits $T$-periodic
solutions $Y(t)$:

$Y(t+T)=Y(t),\qquad t\in{\mathbb R}.$ (5.4)

For $Y=(A,B,C_{*})\in{\mathbb Y}$, let us denote

$Y^{1}=(A,B)\in{\mathbb R}^{2},\qquad Y^{2}=C_{*}\in S^{2}.$ (5.5)

######  Remark 5.2.

For the time-periodic solutions (5.4), the corresponding solutions to the
Maxwell–Bloch equations (2.17) with initial data $X(0)=(X^{1}(0),X^{2}(0))$
where $X^{1}(0)=Y^{1}(0)$ and $X^{2}(0)\in h^{-1}Y^{2}(0)$, the Maxwell field,
the current and the population inversion are time-periodic, while the wave
function acquires a unit factor:

$X^{1}(t+T)=X^{1}(t),\qquad X^{2}(t+T)=e^{i\theta(t)}X^{2}(t),\quad
t,\,\theta(t)\in{\mathbb R}.$ (5.6)

Solutions to (4.2) admit the representation

$Y(t)=U(t)Y(0),\qquad t\in{\mathbb R},$ (5.7)

where $U(t):{\mathbb Y}\to{\mathbb Y}$ is the diffeomorphism, and $U(0)={\rm
Id}$ is the identity map. The map $U(T)$ is homotopic to the identity map
since $U(t)$ depends continuously on $t\in{\mathbb R}$. We will denote also by
$U$ the homotopy map $U:{\mathbb Y}\times[0,T]\to{\mathbb Y}$ defined by

$U:(Y,t)\mapsto U(t)Y.$ (5.8)

Theorem 5.1 is equivalent to the fact that the Poincaré map $U(T)$ admits at
least one fixed point:

$U(T)Y=Y.$ (5.9)

To construct such fixed point, we debelop the transversality theory of
differential topology [8, 11]. Let us present our strategy. We will consider
the manifold $S^{2}$ as the unit sphere in ${\mathbb R}^{3}$. Define the
“rotations” $R_{\varepsilon}$ of the phase space ${\mathbb Y}$ by

$R_{\varepsilon}Y=(R_{\varepsilon}^{F}Y^{1},R_{\varepsilon}^{S}Y^{2}),\qquad
Y=(Y^{1},Y^{2})\in{\mathbb Y}.$ (5.10)

Here $R_{\varepsilon}^{S}$ is the rotation of the space $R^{3}$ by the angle
$\varepsilon$ about an arbitrary vector $\xi\in S^{2}$, while
$R_{\varepsilon}^{F}$ stands for the rotation of the plane ${\mathbb R}^{2}$
by the angle $\varepsilon$ in the coordinates (3.11): for
$Y^{1}=Z^{+}e^{+}+Z^{-}e^{-}$,

$R_{\varepsilon}^{F}Y^{1}=Z^{+}_{\varepsilon}e^{+}+Z^{-}_{\varepsilon}e^{-},\qquad\left(\\!\\!\begin{array}[]{c}Z_{\varepsilon}^{+}\\\
Z_{\varepsilon}^{-}\end{array}\\!\\!\right)=R(\varepsilon)\left(\\!\\!\begin{array}[]{c}Z^{+}\\\
Z^{-}\end{array}\\!\\!\right),$ (5.11)

where $R(\varepsilon)$ is the matrix (3.12). For small $\varepsilon>0$, the
map $R_{\varepsilon}$ has exactly two fixed points
$p^{\pm}=(0,\pm\xi)\in{\mathbb Y}$. The points are nondegenerate, i.e., the
Jacobian $\partial(R_{\varepsilon}Y-Y)/\partial Y\neq 0$ at the points:

$\det[DR_{\varepsilon}(p^{\pm})-E]\neq 0.$ (5.12)

This follows from the fact that in suitable coordinates on $S^{2}$, the
differential
$DR_{\varepsilon}=\left(\\!\\!\\!\begin{array}[]{cc}R(\varepsilon)\\!\\!&\\!\\!0\\\
0\\!\\!&\\!\\!R(\varepsilon)\end{array}\\!\\!\\!\right)$ has only pure
imaginary eigenvalues. Further, define the maps

$U_{\varepsilon}(t)Y=U(t)R_{\varepsilon}Y.$ (5.13)

Note that $U_{\varepsilon}(T)Y=U(T)R_{\varepsilon}Y$, while
$U_{\varepsilon}(0)=R_{\varepsilon}$ has two nondegenerate fixed points:

$U_{\varepsilon}(0)p_{\pm}=p_{\pm},\qquad\det[DU_{\varepsilon}(0)(p^{\pm})-E]\neq
0.$ (5.14)

Further analysis relies on the study of the following sets of fixed points of
all maps $U_{\varepsilon}(t)$

$\Phi_{\varepsilon}=\\{(Y,t)\in{\mathbb
Y}_{T}:U_{\varepsilon}(t)Y=Y\\},\qquad{\mathbb Y}_{T}:={\mathbb
Y}\times[0,T].$ (5.15)

In the next sections we will prove the following proposition.

######  Proposition 5.3.

Let (5.1) and (5.3) hold. Then for $\varepsilon\in(0,\varepsilon_{*})$ with
sufficiently small $\varepsilon_{*}>0$,
i) The set $\Phi_{\varepsilon}$ contains some points $(Y_{\varepsilon},T)$.
ii) The set of all the points $Y_{\varepsilon}$ is bounded in ${\mathbb Y}$.

This proposition immediately implies Theorem 5.1. Indeed,
$U_{\varepsilon}(T)Y_{\varepsilon}=U(T)R_{\varepsilon}Y_{\varepsilon}=Y_{\varepsilon}$,
so any limit point $Y=\lim_{\varepsilon\to 0}Y_{\varepsilon}$ is the fixed
point of $U(T)$.

The proposition can be reformulated in terms of the intersection theory as
follows. Denote the direct product ${\mathbb W}={\mathbb Y}\times{\mathbb Y}$,
its diagonal $\Delta=\\{(Y,Y):Y\in{\mathbb Y}\\}$ and the graph
$\Gamma_{\varepsilon}(t)=\\{(Y,U_{\varepsilon}(t)Y:Y\in{\mathbb Y}\\}$ of the
map $U_{\varepsilon}(t):Y\mapsto U_{\varepsilon}(t)Y$. Every point of the
intersection $\Delta\cap\Gamma_{\varepsilon}(t)$ is the fixed point of
$U_{\varepsilon}(t)$ and vice versa. So, Proposition 5.3 can be reformulated
as follows: for $\varepsilon\in(0,\varepsilon_{*})$ with sufficiently small
$\varepsilon_{*}>0$,
i) $\Delta\cap\Gamma_{\varepsilon}(T)\neq\emptyset$.
ii) The sets $\Delta\cap\Gamma_{\varepsilon}(T)$ are uniformly bounded in
${\mathbb W}$.

##  6 A priori bound for fixed points

The boundedness in Proposition 5.3 ii) is provided by the following lemma. Let
us denote the subsets of ${\mathbb Y}$ which are “neighborhoods of infinity”:

${\mathbb Y}_{B}=\\{Y\in{\mathbb Y}:|Y^{1}|>B\\},\qquad B>0.$ (6.1)

Recall that we consider the manifold $S^{2}$ as the unit sphere in ${\mathbb
R}^{3}$, so the phase space ${\mathbb Y}$ is the subset of ${\mathbb R}^{5}$.

######  Lemma 6.1.

Let (5.1) and (5.3) hold. Then for large $B>0$ and sufficiently small
$\varepsilon_{*}>0$, the lower bound holds,

$\inf_{\begin{array}[]{c}Y\in{\mathbb Y}_{B}\\\
t\in[0,T]\end{array}}|U_{\varepsilon}(t)Y-Y|>0,\qquad\varepsilon\in(0,\varepsilon_{*}).$
(6.2)

###### Proof.

We deduce the bound from the asymptotics (3.16). In the representation (3.11),

$Y^{1}_{\varepsilon}(t):=U_{\varepsilon}(t)Y=Z^{+}_{\varepsilon}(t)e^{+}+Z^{-}_{\varepsilon}(t)e^{-},\qquad
Y^{1}_{\varepsilon}=Z^{+}_{\varepsilon}e^{+}+Z^{-}_{\varepsilon}e^{-}.$ (6.3)

The initial state $Y^{1}_{\varepsilon}(0)=R(\varepsilon)Y^{1}_{\varepsilon}$
is defined by (5.11):

$\left(\\!\\!\begin{array}[]{c}Z^{+}_{\varepsilon}(0)\\\
Z^{-}_{\varepsilon}(0)\end{array}\\!\\!\right)=R(\varepsilon)\left(\\!\\!\begin{array}[]{c}Z^{+}_{\varepsilon}\\\
Z^{-}_{\varepsilon}\end{array}\\!\\!\right).$ (6.4)

The second asymptotics in (3.16) implies that in the polar coordinates
$Z^{+}_{\varepsilon}=r_{\varepsilon}\cos\phi_{\varepsilon}$,
$Z^{-}_{\varepsilon}=r_{\varepsilon}\sin\phi_{\varepsilon}$ and
$Z^{+}_{\varepsilon}(t)=r_{\varepsilon}(t)\cos\phi_{\varepsilon}(t)$,
$Z^{-}_{\varepsilon}(t)=r_{\varepsilon}(t)\sin\phi_{\varepsilon}(t)$, we have
for $t\in[0,T]$

$\phi_{\varepsilon}(t)=\phi_{\varepsilon}(0)+\tilde{\Omega}t+O(t/r_{\varepsilon}(0))=\phi_{\varepsilon}+\tilde{\Omega}\varepsilon+\tilde{\Omega}t+O(t/r_{\varepsilon}(0))\mod
2\pi$ (6.5)

as $r_{\varepsilon}(0)\to\infty$ since
$\phi_{\varepsilon}(0)=\phi_{\varepsilon}+\tilde{\Omega}\varepsilon\mod 2\pi$
by (6.4). Rewrite as

$\phi_{\varepsilon}(t)-\phi_{\varepsilon}=\tilde{\Omega}\varepsilon+\tilde{\Omega}t+O(t/r_{\varepsilon}(0))\mod
2\pi,\qquad r_{\varepsilon}(0)\to\infty.$ (6.6)

By our assumption (5.3), we have
$\Theta:=\tilde{\Omega}\varepsilon+\tilde{\Omega}T<\Omega_{p}T=2\pi$ for small
$\varepsilon>0$, Hence, for $\delta:=2\pi-\Theta>0$, sufficiently small
$\varepsilon_{*}>0$ and large $B_{1}>0$,

$\phi_{\varepsilon}(t)-\phi_{\varepsilon}\leq 2\pi-\delta/2,\qquad
t\in[0,T],\qquad\varepsilon\in(0,\varepsilon_{*}],\quad|r_{\varepsilon}(0)|>B_{1}.$
(6.7)

On the other hand, for sufficiently large $B_{2}>0$.

$\phi_{\varepsilon}(t)-\phi_{\varepsilon}\geq\tilde{\Omega}\varepsilon>0,\qquad
t\in[0,T],\qquad\varepsilon\in(0,\varepsilon_{*}],\quad|r_{\varepsilon}(0)|>B_{2}.$
(6.8)

Finally, the first asymptotics in (3.16) implies that

$r_{\varepsilon}(t)\geq e^{-\gamma T}r_{\varepsilon}(0),\qquad
t\in[0,T],\quad|r_{\varepsilon}(0)|>B_{3}.$ (6.9)

Hence, (6.7) and (6.8) imply that

$|Z_{\varepsilon}(t)-Z|\geq d_{\varepsilon}>0,\qquad
t\in[0,T],\qquad\varepsilon\in(0,\varepsilon_{*}],\quad|r_{\varepsilon}(0)|>B=\max(B_{1},B_{2},B_{3}),$
(6.10)

where $d_{\varepsilon}\sim\tilde{\Omega}\varepsilon B$. It remains to note
that the constant $B$ in (6.10) does not depend on
$\varepsilon\in(0,\varepsilon_{*}]$, though $d_{\varepsilon}$ does. ∎

##  7 The transversality and orientation

The set $\Phi_{\varepsilon}$ coincides with the preimage of the diagonal
$\Delta$ by the map ${\mathbb Y}_{T}\to{\mathbb W}$ defined by
$\Gamma_{\varepsilon}:(Y,t)\to(Y,U_{\varepsilon}(t)Y)$:

$\Phi_{\varepsilon}=\Gamma_{\varepsilon}^{-1}\Delta.$ (7.1)

By Transversality Theorem [11, Section 3.2], there exist approximations
$\Gamma_{\varepsilon,\delta}:{\mathbb Y}_{T}\to{\mathbb W}$ with $\delta>0$,
which are transversal to $\Delta$, and

$\Gamma_{\varepsilon,\delta}(Y,t)\to\Gamma_{\varepsilon}(Y,t),\qquad\delta\to
0,$ (7.2)

where the convergence holds in $C_{b}({\mathbb Y}\times[0,T],{\mathbb W})$.
Then the corresponding preimage

$\Phi_{\varepsilon,\delta}:=\Gamma_{\varepsilon,\delta}^{-1}\Delta$ (7.3)

is a smooth 1D submanifold in ${\mathbb Y}_{T}$. The submanifold is the set of
all fixed points of the maps

$U_{\varepsilon,\delta}(t)Y=\Gamma_{\varepsilon,\delta}(Y,t),\qquad
t\in[0,T].$ (7.4)

The lower bound (6.2) and the convergence (7.2) imply similar lower bound for
small $\delta>0$:

$\inf_{\begin{array}[]{c}Y\in{\mathbb Y}_{B}\\\
t\in[0,T]\end{array}}|U_{\varepsilon,\delta}(t)Y-Y|>0,\qquad\varepsilon\in(0,\varepsilon_{*}),\quad\delta\in(0,\delta(\varepsilon)).$
(7.5)

Hence, the maps $U_{\varepsilon,\delta}(t)$ do not have fixed points in
${\mathbb Y}_{B}$ for all $t\in[0,T]$. Therefore,

$\sup_{(Y,t)\in\Phi_{\varepsilon,\delta}}|Y^{1}|\leq
B,\qquad\varepsilon\in(0,\varepsilon_{*}),\quad\delta\in(0,\delta(\varepsilon)),$
(7.6)

so the manifolds $\Phi_{\varepsilon,\delta}$ are compact and uniformly
bounded. Let us denote
$M_{\varepsilon,\delta}^{\pm}\subset\Phi_{\varepsilon,\delta}$ the smooth 1D
submanifolds containing the points $(p^{\pm},0)$, respectively.

######  Lemma 7.1.

The submanifolds $M_{\varepsilon,\delta}^{\pm}$ are diffeomorphic to $[0,1]$:

$M_{\varepsilon,\delta}^{\pm}=\\{q_{\varepsilon,\delta}^{\pm}(s):s\in[0,1]\\},\qquad
q^{\pm}(0)=(p^{\pm},0).$ (7.7)

###### Proof.

The submanifolds are compact by (7.6). The circles are forbidden since the
equation $Y=U_{\varepsilon,\delta}(t)Y$ with small $\varepsilon,\delta,t>0$
admits exactly one solution close to $(p^{+},0)$ or to $(p^{-},0)$. Indeed,
(5.14) and (7.2) imply that for sufficiently small $\varepsilon_{\\#}>0$, also

$\det[DU_{\varepsilon,\delta}(t)(p^{\pm})-E]\neq
0,\qquad\varepsilon\in(0,\varepsilon_{\\#}),\quad\delta\in(0,\delta_{\\#}(\varepsilon)),\quad
t\in(0,t_{\\#}(\varepsilon,\delta)).$ (7.8)

Here the matrix $DU_{\varepsilon,\delta}(t)(p^{\pm})-E$ is the Jacobian of the
map $Y\to U_{\varepsilon,\delta}(t)(p^{\pm})Y-Y$. Hence, by the implicit
function theorem, for small $\varepsilon_{\\#}>0$, the set
$\Phi_{\varepsilon,\delta}$ with
$\varepsilon\in(0,\varepsilon_{\\#}),\delta\in(0,\delta_{\\#}(\varepsilon))$,
and $t\in(0,t_{\\#}(\varepsilon,\delta))$, contains exactly two points
$(p^{\pm}(t),t)$ which are close to $(p^{\pm},0)$ respectively. ∎

######  Lemma 7.2.

$q^{\pm}_{\varepsilon,\delta}(1)\in{\mathbb Y}\times T$.

###### Proof.

We streamline the ideas behind the proof of Lemma 1.2 from [11] based on the
orientation arguments. We will omit below the signs $\pm$ not to overburden
the notation. We have the following alternative:

${\rm either}\quad q_{\varepsilon,\delta}(1)\in{\mathbb Y}\times T\qquad{\rm
or}\qquad q_{\varepsilon,\delta}(1)\in{\mathbb Y}\times 0.$ (7.9)

It remains to show that the latter is impossible. The transversality of
approximations $\Gamma_{\varepsilon,\delta}$ to $\Delta$ means, in particular,
that for each $q_{\varepsilon,\delta}(s)\in M_{\varepsilon,\delta}$, the
differential $D\Gamma_{\varepsilon,\delta}(q_{\varepsilon,\delta}(s))$ maps
the tangent spaces ${\cal T}_{q_{\varepsilon,\delta}(s)}{\mathbb Y}_{T}$ into
the space, transversal to the diagonal $\Delta$ at the points

$Q_{\varepsilon,\delta}(s):=\Gamma_{\varepsilon,\delta}(q_{\varepsilon,\delta}(s)).$
(7.10)

By the definition of the transversality,

$D\Gamma_{\varepsilon,\delta}(q_{\varepsilon,\delta}(s))[{\cal
T}_{q_{\varepsilon,\delta}(s)}{\mathbb Y}_{T}]+{\cal
T}_{Q_{\varepsilon,\delta}(t)}\Delta={\cal
T}_{Q_{\varepsilon,\delta}(s)}{\mathbb W},\qquad s\in[0,1].$ (7.11)

The key remark is that the normal bundle ${\cal N}M_{\varepsilon,\delta}$ is
mapped isomorphically into ${\cal N}\Delta$: the map

$D\Gamma_{\varepsilon,\delta}(q_{\varepsilon,\delta}(s)):{\cal
N}_{q_{\varepsilon,\delta}(s)}M_{\varepsilon,\delta}\to{\cal
N}_{Q_{\varepsilon,\delta}(s)}\Delta,\qquad s\in[0,1]$ (7.12)

is the isomorphism. This follows from (7.11) since the tangent to $\Delta$
vector $\partial_{s}Q_{\varepsilon,\delta}(s)$ is mapped into the tangent to
$M_{\varepsilon,\delta}$ vector $\partial_{s}q_{\varepsilon,\delta}(s)$: by
(7.10),

$\partial_{s}Q_{\varepsilon,\delta}(s):=D\Gamma_{\varepsilon,\delta}(q_{\varepsilon,\delta}(s))\partial_{s}q_{\varepsilon,\delta}(s).$
(7.13)

Now we can eliminate the second alternative in (7.9). First, note that the
space ${\mathbb Y}$ is orientable since both $D$ and $S^{2}$ are the ones.
Hence, the space ${\mathbb Y}_{T}$, the submanifold
$M_{\varepsilon,\delta}^{\pm}$ and its normal bundle ${\cal
N}M_{\varepsilon,\delta}^{\pm}$ also are orientable. In the latter case of
(7.9), the orientations of the “normal spaces” ${\cal
N}_{q_{\varepsilon,\delta}(s)}M_{\varepsilon,\delta}$ with $s=0$ and $s=1$ are
different. On the other hand, covering the curve $s\mapsto
Q_{\varepsilon,\delta}(s)$ in $\Delta$ by coordinate charts, we can construct
a continuous family of bases in ${\cal N}\Delta$. Then the isomorphism (7.12)
gives us the corresponding continuous family of the bases in ${\cal
N}_{q_{\varepsilon,\delta}(s)}M_{\varepsilon,\delta}$. This contradiction
proves the lemma. ∎

Lemma 7.2 implies that
$q^{\pm}_{\varepsilon,\delta}(1)=(Y^{\pm}_{\varepsilon,\delta},T)\in\Gamma_{\varepsilon,\delta}^{-1}\Delta$,
i.e.,

$\Gamma_{\varepsilon,\delta}Y^{\pm}_{\varepsilon,\delta}=(Y^{\pm}_{\varepsilon,\delta},U_{\varepsilon,\delta}(T)Y^{\pm}_{\varepsilon,\delta})\in\Delta.$
(7.14)

Equivalently,

$U_{\varepsilon,\delta}(T)Y^{\pm}_{\varepsilon,\delta}=Y^{\pm}_{\varepsilon,\delta}.$
(7.15)

The set of all $Y^{\pm}_{\varepsilon,\delta}$ is bounded by (7.6). Hence, for
any limit point as $\delta\to 0$, we obtain by (7.2) that

$U_{\varepsilon}(T)Y^{\pm}_{\varepsilon}=Y^{\pm}_{\varepsilon}.$ (7.16)

Thus, Proposition 5.3 is proved and Theorem 5.1 as well.

######  Remark 7.3.

i) The arguments above provide the existence of al least two different fixed
points and stationary solutions.
ii) The fixed points (7.15) are constructed as the ends of the curves
$s\mapsto q^{\pm}_{\varepsilon,\delta}(s)$ where each point
$q^{\pm}_{\varepsilon,\delta}(s)\in
M_{\varepsilon,\delta}^{\pm}\subset\Phi_{\varepsilon,\delta}$ is the fixed
point of $U_{\varepsilon,\delta}(t)$ with a suitable $t\in[0,T]$ by (7.1).

## References

  * [1] L. Allen, J.H. Eberly, Optical Resonance and Two-Level Atoms, Dover, New York, 1987.
  * [2] F.T. Arecchi, Chaos and generalized multistability in quantum optics, Phys. Scr. 9 (1985), 85–92.
  * [3] I. Bejenaru, D. Tataru, Global wellposedness in the energy space for the Maxwell–Schrödinger system, Commun. Math. Phys. 288 (2009), 145–198.
  * [4] M.R. Candido, J. Llibre, D.D. Novaes, Persistence of periodic solutions for higher order perturbed differential systems via Lyapunov–Schmidt reduction, Nonlinearity 30 (2017), 35–60.
  * [5] M.R. Candido, J. Llibre, Stability of periodic orbits in the averaging theory: Applications to Lorenz and Thomas’ differential systems, Int. J. Bifurcat. Chaos 28 (2018), 1830007–14.
  * [6] M. R. Candido, J. Llibre, C. Valls, New symmetric periodic solutions for the Maxwell–Bloch differential system, Math. Phys. Anal. Geom. 20 (2019), no. 16.
  * [7] V.V. Chepyzhov, M.I. Vishik, Attractors for Equations of Mathematical Physics, AMS, Providence, RI, 2000.
  * [8] V. Guillemin, A. Pollack, Differential Topology, Prentice Hall, New Jersey, 1974.
  * [9] Y. Guo, K. Nakamitsu, W. Strauss, Global finite-energy solutions of the Maxwell–Schrödinger system, Comm. Math. Phys. 170 (1995), no. 1, 181–196.
  * [10] H. Haken, Laser Theory, Springer, Berlin, 1984.
  * [11] M.W. Hirsch, Differential Topology, Springer, New York, 1976.
  * [12] R. D. Jackson, Classical Electrodynamics, Wiley, New York, 1999.
  * [13] W. Koechner, M. Bass, Solid-State Lasers, Springer, New York, 2003.
  * [14] A. Komech, On parametric resonance in the laser action, 2022. arXiv:2208.10179 [quant-ph].
  * [15] A.I. Komech, On quantum jumps and attractors of the Maxwell–Schrödinger equations, Annales mathematiques du Quebec 46 (2022), 139–159.
  * [16] A. Komech, On absorbing set for 3D Maxwell–Schrödinger damped driven equations in bounded region, preprint, 2021. arXiv: 2104.10723.
  * [17] A. Komech, Quantum Mechanics: Genesis and Achievements, Springer, Dordrecht, 2013.
  * [18] A.I. Komech, Lectures on Quantum Mechanics and Attractors, World Scientific, Singapore, 2022.
  * [19] A. Komech, E. Kopylova, Attractors of Hamiltonian Nonlinear Partial Differential Equations, Cambridge University Press, Cambridge, 2021.
  * [20] G.L. Lamb, Analytical descriptions of ultrashort optical pulse propagation in a resonant medium, Rev. Mod. Phys. 43 (1971), 99–124.
  * [21] H. Metcalf, Laser Cooling and Trapping, Springer, 1999.
  * [22] I.I. Rabi, Space quantization in a gyrating magnetic field, Physical Review 51 (1937), no. 8, 652–654.
  * [23] J.J. Sakurai, Advanced Quantum Mechanics, Pearson Education, Incorporated, 2006.
  * [24] M. Sargent III, M.O. Scully, W.E. Lamb Jr, Laser Physics, Addison Wesley, Reading, 1978.
  * [25] L.I. Schiff, Quantum Mechanics, McGraw-Hill, New York, 1955.
  * [26] A.E. Siegman, Lasers, University Science Books, 1986.
  * [27] W.T. Silfvast, Laser Fundamentals, Cambridge University Press, Cambridge, 2012.
  * [28] M.O. Scully, M.S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, 1997.
  * [29] R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, Springer, New York, 1997.
  * [30] J. Wei, X. Wang, X. Geng, Periodic and rational solutions of the reduced Maxwell–Bloch equations, Comm. Nonl. Sci. & Numer. Simulat 59 (2018), 1–14.
  * [31] G.K. Woodgate, Elementary Atomic Structure, Clarendon Press, Oxford, 2002.


# Demand Management for Peak to Average Ratio Minimization via Intraday Block
Pricing

Carolina Cortez, Andreas Kasis, and Stelios Timotheou This work was funded by
the European Union’s Horizon 2020 research and innovation program under grant
agreements No. 891101 (SmarTher Grid) and No. 739551 (KIOS CoE),  the project
EnergyXchange - Enterprises/0618/0164 and from the Republic of Cyprus through
the Directorate General for European Programs, Coordination, and
Development.C. Cortez, A. Kasis and S. Timotheou are with the KIOS Research
and Innovation Center of Excellence and the Department of Electrical and
Computer Engineering, University of Cyprus, CY-1678 Nicosia, Cyprus (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>timotheou.stelios@ucy.ac.cy).

###### Abstract

Price based demand response schemes may significantly improve power system
efficiency. Additionally, it is desired that such schemes yield improved power
operation, by reducing the peak consumption. This paper proposes the Intraday
Block Pricing (IBP) scheme, aiming to promote effective demand response among
consumers by charging their electricity usage based on intraday time-slots. To
design the prices associated with the proposed scheme, we formulate a bilevel
optimization problem that aims to minimize the Peak-to-Average Ratio (PAR) and
simultaneously benefit the consumers and the utility company. The bilevel
problem is converted into a single-level Mathematical Program with Equilibrium
Constraints (MPEC). The resulting MPEC is non-convex and includes nonlinear
constraints. Hence, to obtain a solution, it is relaxed into a Mixed Integer
Linear Program by dealing with all nonlinearities. To evaluate the
conservativeness of the proposed approach, a lower bound to the cost of the
original bilevel problem is obtained. The applicability of the proposed scheme
is demonstrated with simulations on various case studies, which exhibit a
significant reduction in PAR and economic gains for the utility company and
consumers. Furthermore, simulation results show that the solutions of the
original and relaxed problems are equivalent, demonstrating the effectiveness
of the proposed solution approach.

###### Index Terms:

Regulated electricity retail markets, optimal pricing scheme design, demand
response programs, peak-to-average ratio reduction, demand shifting
flexibility.

##  I Introduction

Motivation: Demand response programs are expected to become widespread over
the next years, particularly in the residential sector [1]. Although
individual household loads are typically small, their large numbers result in
a substantial aggregate consumption. In particular, the residential sector
accounts for almost 30% of the electricity market share of which nearly 80%
are currently being operated by regulated utility companies under concession
agreements [2]. In such agreements, the regulator allows a company to provide
power in a certain area. However, to avoid market power abuse, the regulator
determines the retail prices. Hence, the design of effective market mechanisms
that enable efficient demand response among households is highly important.
Such schemes may also take into account additional regulatory rules that aim
to improve the operation of the power network.

Price-based demand response programs use price signals to alter the behaviour
of electricity consumers, enabling the time shifting of electricity demand.
Such schemes may be used to achieve a more balanced daily load shape,
effectively reducing the peak-to-average ratio (PAR). The latter may yield
multiple benefits, since a reduced PAR is associated with rarer load-shedding
and power outage events, improved power quality, and more efficient use of the
power system assets [3, 1, 4, 5]. The need for PAR reduction is expected to
grow as more plug-in hybrid electric vehicles are incorporated in the power
grid, since those are expected to significantly add to the peak consumption
[5, 6]. The appropriate utilization of demand side resources through the
curtailment of PAR may also reduce the need for building expensive backup
generators and help to accommodate further penetration of intermittent
renewable generation.

Shifting flexible loads may also reduce household individual energy costs. In
particular, consumers may reduce their electricity bill by changing their
consumption patterns in accordance with the power  price signals. Hence, in
addition to improving the power system performance, price-based demand
response programs can also benefit consumers [1, 7]. Hence, the design of
suitable price schemes that will simultaneously minimize the PAR and reduce
the consumers costs is highly important.

Literature Review: Smart grid technologies are expected to enable the
widespread adoption of price-based programs through enhanced monitoring,
communication , and control. Among the proposed types of pricing schemes,
those based on time-varying prices are largely discussed as potential
solutions to promote demand response among consumers. For example, real-time
pricing (RTP) can reflect the continuous variation in the cost of electricity
supply, which is a desirable attribute that is attracting growing popularity
[8]. However, the volatility of prices might compromise the transparency of
the scheme, making it hard for consumers that are not generally equipped with
the required technology and knowledge (e.g. small households) to effectively
engage with dynamic markets. Previous studies showed that continuous price
changes might result in consumer fatigue, making them responsive to only
extreme price changes [9, 10]. Furthermore, RTP schemes may raise privacy
concerns due to the disclosure of confidential information from the
subscribers.

An alternative scheme, that yields more predictable price patterns, is called
Time-of-Use (ToU). The ToU pricing scheme includes a set of different prices
for different hours of the day, or even different seasons of the year, defined
at the beginning of a tariff cycle (e.g., a year), that  remains constant
throughout the cycle. ToU schemes have demonstrated a shift in demand from
periods of peak to periods of valley consumption [11]. A concern in such
schemesis that they may result in the concentration of flexible loads at
lower-priced periods, potentially causing larger load profile imbalances [12].

Another price-based program is the Inclining Block Rate (IBR) scheme, where
the electricity use is charged in consumption-blocks. In such schemes,
consumers that have higher electricity consumption throughout a charging
period, typically one or two months, are subjected to higher prices. This
strategy has been used in different markets to address supply shortage issues
[13] or promote social equality [14]. However, IBR schemes do not necessarily
reduce the PAR in power systems [15].

The use of price schemes as a means to reduce the PAR has been considered in
the literature. In particular, [5, 4, 16, 17] propose demand-side management
schemes that schedule the load profile of consumers through price signals  to
minimize the PAR. However, the proposed schemes raise issues of scalability
and transparency since consumers are required to engage with dynamic prices
and communicate their scheduled load profiles. Hence, an aim of this work  is
to design predictable prices for households that reduce individual electricity
bills and minimize the PAR of the power grid.

Contribution: This paper proposes a novel pricing scheme, called the Intraday
Block Pricing (IBP)  scheme that aims to enable  improved demand response. The
IBP scheme mechanism resembles that of the IBR pricing scheme; however,
instead of considering consumption-blocks on a monthly basis, it considers
consumption-blocks in intraday time-slots (e.g. one or two hours). The IBP
scheme aims to yield higher prices at times of high consumption, motivating a
shift of the demand to times of lower consumption. The selection of prices for
the IBP scheme aims to minimize the PAR and simultaneously provide a benefit
to the consumers participating, compared to a flat price scheme. The proposed
pricing scheme shares benefits with the RTP, ToU and , IBR schemes. Table I
provides a comparison of the IBP scheme with the existing schemes proposed in
the literature.

Table I: Comparison of Pricing Schemes. Features | Pricing Schemes
---|---
| Flat | RTP | ToU | IBR | IBP
Time |  |  |  |  | Intraday
Horizon | Constant | Hour | Day | Month | pricing
(Typically) |  |  |  |  | period
Energy |  |  |  | ✓ | ✓
conservation |  |  |
PAR |  |  |  |  |
reduction |  | ✓ | ✓ |  | ✓
Non concentration |  |  |  |  | ✓
effect of load shifting |  |  |  |
Predictability | ✓ |  | ✓ | ✓ | ✓
Promotion of |  |  |  | ✓ | ✓
social equity |  |  |
Non disclosure of |  |  |  |  |
private information | ✓ |  | ✓ | ✓ | ✓
Response to |  | ✓ |  |  |
real-time grid state |  |  |  |
Non scalability |  |  |  |  |
issues | ✓ |  | ✓ | ✓ | ✓
Simplicity | ✓ |  | ✓ | ✓ | ✓

The IBP scheme encourages consumers to manage their demand within each time-
slot, aiming to distribute the demand throughout the day. The price design
enables the minimization of the PAR of the power grid and simultaneously
ensures economic benefit for  the consumers and revenue adequacy for the
utility company. Since the prices directly depend on the demand at that time,
the IBP scheme prevents the concentration of loads at some particular time
period. An important advantage of the IBP scheme is its simplicity, since its
design enables price predictability, making it easy for consumers and
policymakers to implement.

For the implementation of the IBP scheme, we propose an optimization-based
approach that investigates, from the point of view of a regulated utility
company, how the parameters of the IBP scheme should be efficiently designed
in a residential market. To accommodate  the consumers behaviour in our
approach, an estimate of the demand response is integrated  into the design of
the IBP prices. The hierarchical relationship between the design of the IBP
prices (leader) and the demand response of residential consumers (followers)
is modeled as a bilevel optimization problem.  To solve this problem, it is
converted into an equivalent single-level Mathematical Program with
Equilibrium Constraints (MPEC). The resulting MPEC is non-convex and includes
nonlinear constraints. Hence, to obtain a solution, it is relaxed into a Mixed
Integer Linear Program (MILP). Finally, a lower bound to the bilevel
optimization problem is obtained  to evaluate the  conservativeness of the
solution resulting from the MILP problem.

The practicality of the proposed approach is demonstrated with realistic
simulations, which show that the presented scheme yields significantly reduced
PAR and economic gains for the utility company and the consumers. In addition,
the solution obtained from the relaxed MILP problem matches the estimated
lower bound to the original non-convex bilevel problem. Hence, the simulation
results demonstrate both the applicability of the proposed scheme and the
effectiveness of the solution approach.

In summary, the main contributions of this paper are:

  * •

The development of the IBP scheme which aims to enable improved price-based
demand response.

  * •

The formulation of a suitable bilevel optimization problem for the efficient
design of the IBP scheme parameters in a regulated residential market. The
problem aims to achieve PAR minimization and simultaneously ensure revenue
adequacy for the utility company and economic benefit for the consumers.

  * •

The development of a solution approach for the considered non-convex bilevel
problem, by relaxing it to a MILP problem.  To assess the conservativeness of
the solution approach, an iterative procedure to obtain a lower bound to the
original problem is also developed. Numerical results demonstrate that the
relaxed MILP problem  yields globally optimal solutions to the original
problem.

  * •

The provision of quantitative evidence  for the beneficial impact of the IBP
scheme on the residential consumers and the power system using suitable case
studies with realistic data.

Paper structure: The remainder of this paper is structured as follows. The
proposed IBP scheme is presented in Section II. Section III describes and
formulates the problem for the efficient selection of the IBP prices in a
regulated residential market, while Section IV outlines the solution approach.
Case studies that demonstrate the beneficial impact of the IBP scheme are
discussed in Section V. Finally, conclusions are drawn in Section VI.

##  II Intraday Block Pricing Scheme

The IBP scheme, illustrated in Fig. 1, consists of a price structure comprised
of the following parameters:

  * •

Time-slot - the day is divided into $T$ time-slots of equal duration. The set
of time slots is denoted by $\mathcal{T}$.

  * •

Clusters - consumers with similar preferences are classified into clusters.
Within each cluster $c\in\mathcal{C}$, where $\mathcal{C}$ is the set of
clusters, there are $n_{c}$ consumers, such that
$N\\!=\\!\sum_{c\in\mathcal{C}}n_{c}$ is the total number of consumers.

  * •

Consumption-blocks - the demand response to the IBP scheme for time-slot $t$
and cluster $c$, denoted by $d^{s}_{tc}$, is split among consumption-blocks
such that $d^{s}_{tc}=\sum_{f\in\mathcal{F}}d^{s}_{tcf}$, where $d^{s}_{tcf}$
is the demand response to the IBP scheme for time-slot $t$, cluster $c$, and
consumption-block $f$, and $\mathcal{F}$ the set of consumption-blocks. We
also define the set $\mathcal{F}^{*}=\\{1,\cdots,F-1\\}$, where
$F=|\mathcal{F}|$ 111Note that for a set $\mathcal{A}$, we use $|\mathcal{A}|$
to denote its cardinality..

  * •

Electricity prices - consumption-blocks are associated with positive
increasing price values satisfying
$\lambda^{s}_{f+1}=\lambda^{s}_{f}+\xi_{f}$,
$\forall(f\\!\in\\!\mathcal{F}^{*})$, where $\xi_{f}>0$ is the price increment
from consumption-block $f$ to $f+1$, and $\lambda^{s}_{f}$ the price
associated with consumption-block $f$.

  * •

Energy consumption breakpoints - each consumption-block $f$ is limited by an
amount of electricity demand $q_{f}>0$ such that $0\leq d^{s}_{tcf}\leq
q_{f}$, $\forall(f\in\mathcal{F}^{*})$ and $d^{s}_{tcF}\geq 0$. The last
consumption-block has no upper limit as the pricing scheme should not limit
the total demand. Note that since $\lambda_{f+1}>\lambda_{f}$, then
$d^{s}_{tc(f+1)}>0$ implies that $d^{s}_{tcf}=q_{f}$, i.e. for a non-zero
consumption at block $f+1$, block $f$ must have been filled.

Fig. 1: Schematic representation of the IBP scheme structure. The horizontal
and vertical axes indicate the consumer electricity demand per pricing period
and the IBP prices  respectively.

##  III IBP price structure design problem

In this section,  describe the considered problem, which aims to yield
appropriate design parameters for the proposed IBP scheme. The latter is
facilitated by formulating and solving a suitable optimization problem.

###  III-A General Problem Description

We consider a smart power system with multiple residential consumers and a
single utility company  that, following the directives of a regulatory entity,
designs the IBP prices. The utility company aims to replace an existing flat
pricing scheme with an (opt-in) IBP scheme. To increase the acceptance among
consumers and to promote a smooth transition towards the IBP scheme, an
incentive mechanism that implicitly protects consumer bills is integrated
into the design of the IBP prices. Revenue adequacy requirements are also
considered to ensure the long-term sustainability of the utility company, i.e.
the total revenue of the utility company must cover the cost of supply and
allow a fair return. The desired minimum rate of return, denoted by $r$, is
determined by the regulatory entity and it is typically based on the total
value of the utility company’s assets (e.g., plant, equipment, working
capital, and deductions for accumulated depreciation).

To estimate the supply cost for the residential consumers, we assume that its
value is proportional to the supplied energy. Thereby we consider exogenous
parameters $\lambda^{w}_{t}$ that reflect the energy cost rate per time-slot
$t$. This method is widely adopted in the literaturesince it captures the
operation of an existing wholesale market from where the utility company buys
electricity [18, 19, 20]. Typically, the energy cost rates (also known as
wholesale prices) result from a strictly increasing, convex cost function.
That is, the cost rates are high during peak demand due to the commitment of
expensive generators, while during low demand periods more efficient
generators are operating thus the cost rates are lower.

To ensure that consumers are motivated to act in a manner consistent with the
rules established by the IBP scheme, the utility company integrates  the
consumers demand response in the process of designing the IBP price structure.
Towards this, we consider a communication channel from consumers to the
utility company for billing purposes, which conforms with widely applied smart
metering practices [21, 22].  It should be noted that, beyond billing
purposes, household privacy is preserved, since no disclosure of information
about the existence or assignment of local energy devices is required. Based
on this communication channel, the utility company can use the periodically
collected data, associated with the flat prices, to estimate the load profile
of consumers. In addition, consumers are assumed to react rationally to
prices, aiming to maximize their welfare [23].

Although a large number of consumers is expected to implement the proposed
approach, we consider a relatively small number of consumer types to represent
them, both for simplicity and computational efficiency. Using the demand data
from flat pricing schemes, consumers are classified into clusters. The
clustering method that could be adopted by the utility company is itself a
widely discussed problem that is beyond the scope of this paper [21, 24].

To estimate the consumers response to the IBP scheme, we introduce an
optimization problem that aims to minimize the electricity bill and load
shifting inconvenience of consumers. The problem considers an aggregated,
technology-agnostic demand response model, similar to [18, 25, 26], and
assumes rational consumer behaviour. It is assumed that the overall energy
demand from consumers is inelastic, based on empirical evidence [27]. Instead,
our approach focuses on encouraging a time shift in consumption. Although
individual demand is in general stochastic, we assume that consumer clustering
results in a predictable aggregate demand response.

The utility company integrates the estimated demand response of consumers on
its price-design problem. The hierarchical relationship between the design of
the IBP price structure and the demand response of consumers is modeled as a
bilevel optimization problem, as illustrated in Fig. 2.

###  III-B Problem statement

In this section , we state the main problem considered in the paper associated
with selecting suitable parameters for the IBP scheme by the utility company.

######  Problem 1.

Select the parameters of the IBP price scheme to:

  1. (R1)

Minimize the PAR of the power grid.

  2. (R2)

Achieve revenue adequacy for the utility company.

  3. (R3)

Protect consumer bills, ensuring no increase in individual bills due to the
adoption of the IBP scheme.

  4. (R4)

Reflect the demand shifting flexibility of consumers.

Problem 1 aims to design the parameters of the IBP scheme such that
Specifications (R1)-(R4) are satisfied. In particular, it  intends to minimize
the PAR of the power system by promoting load shifting among households. In
addition, it aims to achieve revenue adequacy by ensuring that the monetary
amount collected from the consumers covers retail costs and allows a fair
return. A further objective is to promote the adoption of the IBP scheme by
the consumers, through an incentive mechanism that implicitly protects their
bills. Finally, an effective scheme should integrate the estimated demand
shifting flexibility of consumers in its price structure design.

Fig. 2: Bilevel optimization problem for the design of the IBP price
structure.

###  III-C Utility company problem

The optimization problem considered by the utility company aims to satisfy
Specifications (R1)-(R4), by suitable design of the IBP scheme parameters.
The problem can be stated as

$\displaystyle\min_{\Xi_{U}}$
$\displaystyle\frac{\displaystyle\max_{t\in\mathcal{T}}\Bigg{\\{}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}\\!n_{c}d^{s}_{tcf}\Bigg{\\}}}{\displaystyle\frac{1}{T}\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}n_{c}d^{s}_{tcf}}$
(1a) s.t.
$\displaystyle\lambda^{s}_{f}=\lambda^{s}_{1}+(f-1)\xi,$$\displaystyle\forall(f\\!\in\\!\mathcal{F}),$
(1b) $\displaystyle\check{q}\leq
q_{f}\leq\hat{q},$$\displaystyle\forall(f\\!\in\\!\mathcal{F}^{*}),$ (1c)
$\displaystyle\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}\\!\lambda^{s}_{f}n_{c}d^{s}_{tcf}\geq
r\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}\lambda^{w}_{t}n_{c}d^{s}_{tcf},$
(1d)
$\displaystyle\sum_{t\in\mathcal{T}}\sum_{f\in\mathcal{F}}\lambda^{s}_{f}d^{o}_{tcf}\leq\lambda^{o}\sum_{t\in\mathcal{T}}D_{tc},$
$\displaystyle\forall(c\\!\in\\!\mathcal{C}),$ (1e)

where $d^{0}_{tcf}$ represents the demand at time-slot $t$ for cluster $c$
associated with consumption-block $f\in\mathcal{F}$, given by

$\displaystyle
d^{0}_{tcf}\\!=\\!\\!\min\\!\\!\left(\\!\\!\max\\!\left(\\!D_{tc}\\!-\\!\sum_{j=1}^{f-1}\\!q_{j},\\!0\\!\right)\\!\\!,q_{f}\\!\\!\right)\\!\\!,\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F^{*}}),$
$\displaystyle
d^{0}_{tcF}=\max\left(D_{tc}-\sum_{f\in\mathcal{F}^{*}}q_{f},0\right)\\!\\!,\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$

which implies that
$\sum_{f\in\mathcal{F}}d^{o}_{tcf}=D_{tc},\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C})$.

In problem (1), the set of decision variables is denoted by
$\Xi_{U}\\!=\\!\\{\xi;\ q_{f},\ \forall f\in\mathcal{F}^{*};\
\lambda^{s}_{1}\\}$. In addition, $d^{s}_{tcf}$ denotes the estimated demand
response of block $f$, cluster $c$ at time $t$, under the IBP scheme. These
values are estimated by solving a lower level optimization problem, as
explained in Section III-D.

Ojective function (1a) aims to minimize the PAR,  following Specification
(R1). This is achieved by minimizing the peak (nominator of (1a)) over the
average (denominator of (1a)) system demand over the considered time horizon
$\mathcal{T}$.

The price structure of the IBP scheme is captured through constraints (1b) and
(1c). To make the price structure more intuitive and increase the chances of
its acceptance by the consumers, we consider fixed  price increments $\xi$,
such that $\xi_{f}=\xi$. We also consider that the energy consumption
breakpoints have known lower and upper bounds $\check{q}$, $\hat{q}$,
respectively. These bounds can be selected by the utility company or imposed
by the regulatory entity. It should be noted that the upper bound $\hat{q}$
should be selected to be sufficiently large to enable a broad set of solutions
for (1). A suitable option is to let
$\hat{q}=\max_{t\in\mathcal{T},c\in\mathcal{C}}{D_{tc}}$.

To satisfy  Specification (R2), constraint (1d) ensures that the total utility
revenue (left term) is at least as high as the total utility cost (right term)
incurred within the time horizon $T$. As mentioned earlier, the supply cost is
considered proportional to the energy supplied to the consumers,  with cost
rate $\lambda^{w}_{t}$. We also include a minimum rate of return $r\geq 1$,
determined by the regulator to incur a fair return to the utility company.

An incentive mechanism for the adoption of the IBP scheme ( Specification
(R3)) is presented in (1e). The key idea is to ensure that the consumers
electricity bill will not increase under the proposed scheme, compared to the
existing flat price scheme, providing an incentive to opt-in the IBP
scheme222Note also that problem (3) ensures that the resulting demand from
using the IBP scheme satisfies
$\sum_{t\in\mathcal{T}}\sum_{f\in\mathcal{F}}\\!\lambda^{s}_{f}d^{s}_{tcf}\leq\sum_{t\in\mathcal{T}}\sum_{f\in\mathcal{F}}\lambda^{s}_{f}d^{o}_{tcf},\forall
c\in\mathcal{C}$ since
$d^{s}_{tcf}=d^{0}_{tcf},d^{sh}_{tc}=0,\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F})$
is a feasible solution to it. The latter guarantees that the implementation of
the proposed scheme will not result in increased bills for consumers,
regardless  of whether they opt-in or not. . This is ensured by (1e), where
$d^{o}_{tcf}$ describes the distribution of $D_{tc}$ at different consumption-
blocks of the IBP scheme.

######  Remark 1.

The flat price $\lambda^{o}$ is selected to ensure that the market is budget
balanced under  a flat price and is given by

$\lambda^{o}=r\sum\limits_{t\in\mathcal{T}}\sum\limits_{c\in\mathcal{C}}\lambda^{w}_{t}{n_{c}D_{tc}}\Big{/}\sum\limits_{t\in\mathcal{T}}\sum\limits_{c\in\mathcal{C}}n_{c}D_{tc}.$
(2)

Equation (2) ensures that the revenue from selling energy at  a flat price
$\lambda^{o}$ equals the cost of the utility company. For consistency,
parameter $r$ assumes the same value as in (1d). It should be noted that (2)
ensures the existence of a solution to (1), by selecting the parameters of the
IBP scheme to enable the original flat price.

###  III-D Demand response model

The problem considered by the utility company to estimate the aggregated
demand response of consumers,  following Specification (R4) , is presented
below.

$\displaystyle\min_{\Xi_{L}}$
$\displaystyle\sum_{c\in\mathcal{C}}\\!n_{c}\left(\sum_{t\in\mathcal{T}}\sum_{f\in\mathcal{F}}\\!\lambda^{s}_{f}d^{s}_{tcf}\\!+\\!\sum_{t\in\mathcal{T}}\\!\frac{\tau_{c}}{2}{d^{sh}_{tc}}^{2}\right)\\!$
(3a) s.t.
$\displaystyle\sum_{f\in\mathcal{F}}d^{s}_{tcf}\\!-\\!d^{sh}_{tc}\\!=\\!D_{tc}:\rho_{tc},$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (3b)
$\displaystyle 0\leq d^{s}_{tcf}:\mu^{-}_{tcf},$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}),$
(3c) $\displaystyle d^{s}_{tcf}\\!\leq
q_{f}:\mu^{+}_{tcf},\quad$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*}),$
(3d) $\displaystyle-\sigma_{c}D_{tc}\leq
d^{sh}_{tc}:\phi^{-}_{tc},$$\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$
(3e) $\displaystyle
d^{sh}_{tc}\leq\sigma_{c}D_{tc}:\phi^{+}_{tc},$$\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$
(3f) $\displaystyle\sum_{t\in\mathcal{T}}d^{sh}_{tc}=0:\eta_{c},$
$\displaystyle\forall(c\\!\in\\!\mathcal{C}).$ (3g)

In problem (3), the set of decision variables is denoted by
$\Xi_{L}\\!=\\!\\{d^{sh}_{tc},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C});\ d^{s}_{tcf},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F})\\}$.
Variable $d^{sh}_{tc}$ denotes the energy shifted in cluster $c\in\mathcal{C}$
at time-slot $t\in\mathcal{T}$,  while variable $d^{s}_{tcf}$ denotes the
estimated demand response of block $f$, cluster $c$ at time $t$, under the IBP
scheme.

The objective function in (3a) represents the collective cost perceived by
consumers which includes the total consumer bill paid to the utility company
(first term) and the consumer shifting cost due to load shifting in time
(second term). Parameter $\tau_{c}$ indicates the cost coefficient of load
shifting at cluster $c$, describing the level of discomfort that a consumer
experiences when shifting his/her demand. The cost associated with load
shifting is described by a quadratic function since the rate of discomfort
increases when a greater shift in demand is experienced.  The objective
function (3a) inflicts a trade-off between the economic benefit provided by
the IBP scheme prices $\lambda^{s}_{f}$ and the discomfort associated with the
shifting cost coefficient $\tau_{c}$.

Constraints (3b)-(3g) express the time-shifting flexibility of different
clusters using the IBP scheme. Equality (3b) ensures power balance at every
time-slot $t\\!\in\\!\mathcal{T}$ and cluster $c\\!\in\\!\mathcal{C}$ such
that the demand supplied by the utility company to the consumers,
$\sum_{f\in\mathcal{F}}d^{s}_{tcf}$, minus the load shifting $d^{sh}_{tc}$
equals the baseline demand $D_{tc}$. The energy consumption partitioning among
the consumption-blocks of the IBP scheme is expressed through constraints (3c)
and (3d). Constraints (3e) and (3f) provide lower and upper bounds on load
shifting as a percentage $\sigma_{c}$ of the baseline demand. Finally,
constraint (3g) ensures that the demand shifting is energy neutral within the
considered time horizon. Parameters $\sigma_{c}$ in constraints (3e) and (3f)
quantify the amount of flexible demand that can be shifted at each time-slot
$t$.

Consumers can achieve economic benefit by shifting load from consumption-
blocks with higher prices to consumption-blocks with lower prices. Thus, the
economic benefit results from the price increment $\xi$ between consumption-
blocks, as in (1b). Finally, the value of $q_{f}$ limits the amount of demand
in each consumption-block $f$, thus it provides a guideline for the scheduling
of load throughout the day.

It should be noted that (3) is a strictly convex problem. In addition, it is
intrinsically decoupled among consumption clusters, since consumer responses
are independent between clusters. Note also that the dual variables of problem
(3) are indicated after the colon of each constraint. For convenience we
define the set of all dual variables $\Xi_{D}=\\{\mu^{-}_{tcf},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F});\
\mu^{+}_{tcf},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*});\
\phi^{-}_{tc},\ \phi^{+}_{tc},\ \rho_{tc},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C});\ \eta_{c},\
\forall(c\\!\in\\!\mathcal{C})\\}$.

###  III-E Bilevel Optimization Problem for the IBP scheme

The resulting bilevel optimization problem is comprised of the upper level
problem (1) and the lower level problem (3). The upper level problem,
describes the decision-making process of the utility company. For given
estimated demand response of consumers, the parameters of the IBP price
structure are selected to satisfy  Specifications (R1)-(R4). The lower level
problem estimates the demand response of the consumers to a given set of
parameters for the IBP scheme.

The problem defined in (1) and (3) is challenging for two reasons. First, the
problem is bilevel and hence cannot be handled by standard mathematical
programming solvers. Second, the objective is fractional and also there are
two non-convex constraints with bilinear terms ((1d) and (1e)). In the next
section, we develop a solution approach to deal with these issues.

##  IV Solution approach

To solve the bilevel optimization problem introduced in the previous section
we develop a two phase solution approach. First, we replace the lower level
problem (3), by its equivalent Karush-Kuhn-Tucker (KKT)  optimality conditions
transforming the bilevel problem into a single-level Mathematical Program with
Equilibrium Constraints (MPEC). Second, we relax the MPEC into a Mixed Integer
Linear Program (MILP) by dealing with all nonlinearities.

###  IV-A Single level MPEC

To convert the bilevel problem (1),(3) into a single level MPEC, the lower
level problem (3) is replaced by its equivalent KKT conditions, presented in
(4). This transformation is allowed since (3) includes a convex cost function
and linear constraints [28].

$\displaystyle
n_{c}\tau_{c}d^{sh}_{tc}-\\!\phi^{-}_{tc}+\phi^{+}_{tc}+\eta_{c}-\rho_{tc}=0,$
$\displaystyle\ \ \forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (4a)
$\displaystyle
n_{c}\lambda^{s}_{f}-\mu^{-}_{tcf}+\mu^{+}_{tcf}+\rho_{tc}=0,$$\displaystyle\
\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*}\\!),$
(4b) $\displaystyle
n_{c}\lambda^{s}_{f}-\mu^{-}_{tcf}+\rho_{tc}=0,$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!=\\!F),$
(4c)
$\displaystyle\mu^{-}_{tcf}d^{s}_{tcf}=0,$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}),$
(4d)
$\displaystyle\\!\mu^{+}_{tcf}\left(q_{f}\\!-\\!d^{s}_{tcf}\right)=0,$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*}),$
(4e) $\displaystyle\phi^{-}_{tc}\left(d^{sh}_{tc}+\sigma_{c}D_{tc}\right)=0,$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (4f)
$\displaystyle\phi^{+}_{tc}\left(-d^{sh}_{tc}+\sigma_{c}D_{tc}\right)=0,$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (4g)
$\displaystyle\mu^{-}_{tcf}\geq
0,$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}),$
(4h) $\displaystyle\mu^{+}_{tcf}\geq 0,\
$$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*}),$
(4i) $\displaystyle\phi^{-}_{tc}\geq 0,$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (4j)
$\displaystyle\phi^{+}_{tc}\geq 0,$
$\displaystyle\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (4k)
$\displaystyle\text{Constraints }\eqref{eq:LLs2}\text{-}\eqref{eq:LLs7}.$ (4l)

Equalities (4a) and (4b)-(4c) are the optimality KKT conditions resulting from
variables $d^{sh}_{tc}$ and $d^{s}_{tcf}$, respectively. Constraints (4d)-(4g)
are the complementary slackness conditions. Constraints (4h)-(4k) are the dual
feasibility conditions, while (4l) are the primal feasibility conditions [28].
Integrating the KKT conditions (4) into the bilevel problem yields the
following MPEC:

  * Objective function: (1a) (MPEC)

  * Constraints: (1b)-(1e), (3b)-(3g), (4a)-(4k).

  * Variables: $\Xi_{U}$, $\Xi_{L}$ and $\Xi_{D}$.

###  IV-B MILP Reformulation

The MPEC defined in Section IV-A has three types of nonlinearities. The first
nonlinearity type is associated with the presence of decision variables in the
denominator of the objective function (1a), yielding a fractional objective.
This fractional formulation can be transformed into a linear expression. In
particular the lower level problem (3) does not change the total demand of the
households, since summing (3b) over $t\in\mathcal{T}$ and considering (3g)
yields

$\displaystyle
D=\sum_{c\in\mathcal{C}}n_{c}D_{tc}=\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}n_{c}d^{s}_{tcf}$
(5)

where $D$ is a constant, allowing the denominator of the objective function
(1a) to be removed. Hence, objective (1a) can be equivalently expressed as

$\displaystyle\min d^{peak}$ (6a) $\displaystyle
d^{peak}\geq\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}d^{s}_{tcf},\qquad\forall(t\\!\in\\!\mathcal{T}),$
(6b)

where $d^{peak}$ is a new variable.

The second type of nonlinearity is found in constraints (1d) and (1e) due to
the products of variables $\lambda^{s}_{f}d^{s}_{tcf}$ and
$\lambda^{s}_{f}d^{o}_{tcf}$, respectively. To simplify the nonlinear term in
(1d), we substitute (1b) into (1d) and consider (5) yielding

$\displaystyle\sum_{t\in\mathcal{T}}$
$\displaystyle\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}n_{c}\lambda^{s}_{f}d^{s}_{tcf}=$
$\displaystyle=\lambda^{s}_{1}D+\xi\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}n_{c}(f-1)d^{s}_{tcf}.$
(7)

Replacing (IV-B) into (1d) and rearranging the terms yields

$\displaystyle\lambda^{s}_{1}D+\sum_{t\in\mathcal{T}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}(\xi
f-\xi-r\lambda^{w}_{t})n_{c}d^{s}_{tcf}\geq 0.$ (8)

Following a similar approach for constraint (1e) yields

$\displaystyle\left(\lambda^{s}_{1}\\!-\\!\lambda^{o}\right)\sum_{t\in\mathcal{T}}D_{tc}\\!+\\!\xi\\!\sum_{t\in\mathcal{T}}\sum_{f\in\mathcal{F}}(f-1)d^{o}_{tcf}\\!\leq\\!0,\
$ $\displaystyle\forall(c\in\mathcal{C}).$ (9)

All bilinear terms that appear in both (8) and (9) are functions of the
variable $\xi$. Hence, our approach to deal with these nonlinearities is to
fix $\xi$ and solve the problem over different values of
$\xi\in[\check{\xi},\hat{\xi}]$, where $\check{\xi}$ and $\hat{\xi}$ are
sufficiently broad lower and upper bounds of the price increment,
respectively. Consequently , $\xi$ is removed from the set of decision
variables originating from the upper level problem.

The last category of nonlinearities involves the complementary slackness
conditions (4d)-(4g) which are replaced with the respective linear constraints
(10). Note that  $M_{1}$,  $M_{2}$,  $M_{3}$,  $M_{4}$ should be positive
constants that are sufficiently large [20].

$\displaystyle
d^{s}_{tcf}\leq{\color[rgb]{0,0,0}\omega^{1}_{tcf}M_{1}},\quad\mu^{-}_{tcf}\leq(1-{\color[rgb]{0,0,0}\omega^{1}_{tcf}}){\color[rgb]{0,0,0}M_{1}},\
$ $\displaystyle\qquad{\color[rgb]{0,0,0}\omega^{1}_{tcf}}\in\\{0,1\\},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}),$
(10a) $\displaystyle
q_{f}\\!-\\!d^{s}_{tcf}\leq{\color[rgb]{0,0,0}\omega^{2}_{tcf}M_{2}},\quad\mu^{+}_{tcf}\leq(1-{\color[rgb]{0,0,0}\omega^{2}_{tcf}}){\color[rgb]{0,0,0}M_{2}},\
$ $\displaystyle{\color[rgb]{0,0,0}\omega^{2}_{tcf}}\in\\{0,1\\},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*}),$
(10b) $\displaystyle
d^{sh}_{tc}+\sigma_{c}D_{tc}\leq{\color[rgb]{0,0,0}\omega^{3}_{tc}M_{3}},\quad\phi^{-}_{tc}\leq(1-{\color[rgb]{0,0,0}\omega^{3}_{tc}}){\color[rgb]{0,0,0}M_{3}},$
$\displaystyle{\color[rgb]{0,0,0}\omega^{3}_{tc}}\in\\{0,1\\},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (10c)
$\displaystyle-d^{sh}_{tc}+\sigma_{c}D_{tc}\leq{\color[rgb]{0,0,0}\omega^{4}_{tc}M_{4}},\quad\phi^{+}_{tc}\leq(1-{\color[rgb]{0,0,0}\omega^{4}_{tc}}){\color[rgb]{0,0,0}M_{4}},$
$\displaystyle{\color[rgb]{0,0,0}\omega^{4}_{tc}}\in\\{0,1\\},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C}),$ (10d)

For convenience we define the set of all binary variables
$\Xi_{B}=\\{{\color[rgb]{0,0,0}\omega^{1}_{tcf}},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F});\
{\color[rgb]{0,0,0}\omega^{2}_{tcf}},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C},f\\!\in\\!\mathcal{F}^{*});\
{\color[rgb]{0,0,0}\omega^{3}_{tc}},\ {\color[rgb]{0,0,0}\omega^{4}_{tc}},\
\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C})\\}$. Hence, the MPEC
problem is reformulated as a MILP problem as follows:

  * Objective function: (6a) (MILP)

  * Constraints: (1b), (1c), (3b)-(3g), (4a)-(4b), (4h)-(4k), (6b), (8), (9) and (10).

  * Variables: $\Xi_{U}\backslash\\{\xi\\}$, $\Xi_{L}$, $\Xi_{D}$, $\Xi_{B}$ and $d^{peak}$.

In the above formulation , the complementary slackness conditions are treated
using (10), while variable $\xi$ is fixed. Hence, the resulting problem is a
Mixed Integer Linear Program (MILP) which can be solved using standard
optimization tools. The proposed MILP problem is a relaxed version of the MPEC
problem, thus there is no guarantee that we can obtain a value for $\xi$ that
globally minimizes the PAR. To evaluate the conservativeness of our solution
approach, we aim to compare the obtained solution with a theoretical lower
bound of the global minimum of the original problem. The approach to obtain
this lower bound is described in the following section.

###  IV-C Theoretical lower bound for the PAR minimization problem

In this section, we describe our approach to obtain a theoretical lower bound
for the PAR minimization problem, described by (1), (3). In particular, we
consider the load shifting flexibility associated with the lower level problem
(3), and study load profiles that minimize the PAR but possibly violate
(1d)-(1e).

Our approach to obtain the lower bound, is by setting $|\mathcal{F}|=2$ and
letting $\xi$ be sufficiently large, such that consumers have a strong price
incentive to limit their demand at $q_{1}$. Note that, given a fixed value of
$\xi$, the price values do not affect the resulting demand response, although
these are important to satisfy constraints (1d)-(1e). In addition, the case
$|\mathcal{F}|>2$ does not enable a reduced PAR compared to $|\mathcal{F}|=2$.
The latter follows since when $\xi$ is selected to be sufficiently large, then
consumers aim to keep their demand to no more than $q_{1}$ per time slot,
provided this is feasible, i.e. does not violate (3b)-(3g). A potential
addition of price blocks, when $|\mathcal{F}|>2$, would not hence provide
additional incentive for shifting the peak demand, since all shifting ability
is exhausted due to the large value of $\xi$. Considering a larger set of
price blocks could potentially result in improved demand distribution, and is
associated with the satisfaction of (1d)-(1e).

To obtain the value of $q_{1}$, we then consider an iterative  procedure
where, starting with $q_{1}^{(0)}=\check{q}$, the value of $q_{1}^{(i)}$ is
increased by a small positive constant $\varepsilon$ at each iteration $i$,
such that $q_{1}^{(i)}=q_{1}^{(i-1)}+\varepsilon$. The process stops when
$q_{1}^{(i)}\geq\hat{q}$ at some iteration $i$. This approach enables two
prices, with the second being prohibitively expensive, encouraging consumers
to avoid requiring more than $q_{1}$ energy units at some given time. By
considering a large number of values for $q_{1}\in[\check{q},\hat{q}]$, we
consider the variability of this scheme to its changes.

The lower bound for the PAR  is given by
$\underline{PAR}=\min_{{i\in\mathcal{I}}}\max_{t\in\mathcal{T}}{\color[rgb]{0,0,0}\underline{d}^{s,(i)}_{t}}$,
where
${\color[rgb]{0,0,0}\underline{d}^{s,(i)}_{t}}=\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}{{\color[rgb]{0,0,0}d^{s,(i)}_{tcf}}}$,
$\mathcal{I}$ is the set of iterations and ${d^{s,(i)}_{tcf}}$ the demand
associated with time $t$, cluster $c$ and block $f$ at iteration $i$. Note the
the lower bound to the PAR is associated with an iteration $i$ that minimizes
the peak demand when $q_{1}=q_{1}^{(i)}$.

##  V Case studies

This section considers realistic case studies to evaluate the potential impact
of the IBP scheme in minimizing the PAR and benefit the consumers and the
utility company.

###  V-A Simulation Setup

The examined case studies consider the IBP scheme in the context of a single
day time horizon with hourly time-slots (i.e., $T\\!=\\!24$). The case study
is built upon the power system of the United Kingdom (UK). The retail market
consists of $N\\!=\\!1,000$ residential consumers changing from flat pricing
and opting into the proposed pricing scheme. The utility company categorizes
its consumers into four clusters ($C\\!=\\!4$). The selected profiles and
their respective contribution follow from [29] which categorizes UK consumers
into four archetypes named afternoon actives ($c\\!=\\!1$), double risers
($c\\!=\\!2$), winter spinners ($c\\!=\\!3$) and home lunchers ($c\\!=\\!4$).
The hourly baseline demand $D_{tc}$,
$\forall(t\\!\in\\!\mathcal{T},c\\!\in\\!\mathcal{C})$ and the number of
consumers per cluster $n_{c}$, $\forall(c\\!\in\\!\mathcal{C})$ are shown in
Fig. 3(i). The stacked area graphs indicate the aggregated load (summation of
the individual demand of each consumer) with the contribution of each cluster
represented by a different color.

The assumed wholesale market prices follow the pattern of a typical winter day
in the UK as illustrated in Fig. 3(ii) [30]. The remaining operational costs
associated with the utility company’s activity are neglected; hence, it is
assumed that $r\\!=\\!1$.

Fig. 3: (i) Hourly baseline demand of each cluster. (ii) Wholesale prices.

The results from the optimal IBP price structure design approach proposed in
Section IV are compared against a reference case with flat price
$\lambda^{o}\\!=\\!0.083$ £/kWh, calculated using (2) for the considered
scenario. For simplicity, it is assumed that for each scenario the load
shifting flexibility parameters and shifting costs are equal among clusters,
that is $\sigma_{c}=\sigma$ and $\tau_{c}=\tau$ for all clusters
$c\in\mathcal{C}$, respectively. The lower and upper bounds of the energy
consumption breakpoints, such that $q_{f}\in[\check{q},\hat{q}]$, were
selected as follows

$\check{q}=\min_{t\in\mathcal{T},c\in\mathcal{C}}D_{tc},\ \
\hat{q}=\max_{t\in\mathcal{T},c\in\mathcal{C}}D_{tc}.\\\ $ (11)

Equation (11) ensures that all consumption-blocks span through the baseline
demand of different consumer clusters. Note that all optimization problems
have been solved using the mathematical programming software Gurobi [31].

###  V-B Impact of consumer shifting flexibility - $\sigma$

This subsection examines the impact of the consumer shifting flexibility on
the power grid PAR and the economic benefit perceived by the utility company
and consumers. A total of 104 scenarios are considered for varying price
increment values $\xi=\\{0,0.005,\cdots,0.055,0.06\\}$ £/kWh and different
penetration of flexible loads $\sigma=\\{10,20,30,40\\}$ %, for two and three
consumption-blocks ($|\mathcal{F}|\in\\{2,3\\}$). The implemented shifting
cost is $\tau=0.03$ £/kWh2. Note that the impact of $\tau$ is explored in
Section V-C.

Figure 4 depicts the PAR reduction as a percentage of the reference case for
increasing values of the price increment $\xi$. Each plotted line stands for a
different load shifting flexibility $\sigma$. Interestingly, all plots in Fig.
4 have similar behaviour. Initially, the PAR does not present  a significant
reduction for small values of $\xi$, since the price increment is not large
enough to offer an economic incentive that outgrows the consumer shifting
cost. Upon reaching a maximum reduction, the PAR remains constant and then
decreases as $\xi$ increases. This occurs since the revenue adequacy and bill
protection mechanisms, associated with constraints (1d) and (1e) respectively,
limit the strategic choices of the utility company.

The theoretical lower bound for the PAR minimization problem, obtained
following the approach presented in Section IV-C, is depicted in Fig. 4 with a
dashed line. It can be observed that the solution from the relaxed MILP
problem presented in Section IV-B is equivalent to the optimal solution of the
original bilevel problem. Note that, when the price increment $\xi$ becomes
large, then the set of values that minimize the PAR  is no longer feasible
since the constraints (1d)-(1e) are violated.

The presented approach yields the set of values for $\xi$ for each given price
structure and load shifting flexibility that enable s the largest PAR
reduction. For example, if the market described in Section V-A has a load
shifting flexibility of $\sigma=30\%$, then the utility company can
potentially achieve a 25% PAR reduction by selecting the parameters of the IBP
price structure such that (i) $F=2$, $\lambda_{1}=0.079$ £/kWh,
$\xi=\\{0.025,0.03,0.035,0.04\\}$ £/kWh, and $q=1.119$ kWh or (ii) $F=3$,
$\lambda_{1}=0.066$, $\xi=\\{0.025,0.03\\}$ £/kWh, and $q=0.560$ kWh.

Fig. 4: PAR reduction for increasing price increments $\xi$ and different load
shifting flexibility values $\sigma$ when $F=2$ (top) and $F=3$ (bottom).
Table II: Impact of the price structure parameters on variations in economic and physical quantities for selected scenarios Scenarios | S1 | S2 | S3
---|---|---|---
|  |  |  |  |  |
$\sigma$ [%] | 20 | 30 | 30
$\tau$ [£/kWh2] | 0.03 | 0.03 | 0.06
|  |  |  |  |  |
$\xi$ [£/kWh] | 0.03 | 0.03 | 0.03
$F$ [un.] | 2 | 3 | 2 | 3 | 2 | 3
$q_{f}$ [kWh] | 1.145 | 0.501 | 1.119 | 0.689 | 1.118 | 0.533
|  | 0.645 |  | 0.430 |  | 0.585
$\lambda^{s}_{f}$ [£/kWh] | 0.080 | 0.065 | 0.080 | 0.069 | 0.079 | 0.065
| 0.110 | 0.095 | 0.110 | 0.099 | 0.109 | 0.095
|  | 0.125 |  | 0.129 |  | 0.125
|  |  |  |  |  |
PAR [%]∗ | 18.72 | 18.72 | 24.84 | 24.84 | 20.09 | 20.09
Utility cost [%]∗ | 4.14 | 4.14 | 5.88 | 5.96 | 5.07 | 5.07
Consumer bill† [%]∗ | 3.61 | 4.14 | 4.27 | 5.96 | 4.78 | 4.97
Consumer total cost [%]∗ | 2.77 | 3.3 | 2.74 | 4.43 | 2.63 | 2.83
∗: Reduction in percentage compared to the reference case
†: Equivalent to utility revenue

Table II examines the interaction between the price structure design and the
demand response of consumers for three different scenarios (S1-S3) for two
($F=2$) and three ($F=3$) consumption-blocks. Table rows 1-2, 3-6 and 7-10
depict the consumer parameters, the optimal IBP price structure, and the
performance metrics of each scenario respectively. Note that Scenario S2 is
similar to S1 but has increased load shifting flexibility ($\sigma=30\%$
compared to $20\%$), while Scenario S3 is similar to S2 but has increased
shifting cost ($\tau=0.06$ compared to $\tau=0.03$). Comparing Scenarios S1
and S2 shows that a 50% increase in load shifting flexibility results in
approximately 40% reduction in terms of utility cost and consumer bills. The
same comparison for Scenarios S2 and S3 yields a 16% performance reduction in
both metrics. It should be noted that cases $F=2$ and $F=3$ yield the same
maximum PAR reduction, which is in line with the discussion in Section IV-C.
However, the additional flexibility offered by having the extra pricing block
may result in improved utility and consumer costs, as demonstrated in Table
II.

Considering the different metrics across all scenarios yields three important
observations. First, the PAR is reduced by more than 15% compared to the
reference case. Second, the utility operating cost, as reflected by the
consumer bills, is improved by 3-6%. Finally, the results indicate that the
reduction of the consumer bills outgrows the perceived increase on shifting
cost; hence, the consumer experience is improved with the implementation of
the IBP scheme.

To illustrate how consumers respond to the IBP scheme, Fig. 5 provides the
hourly demand of each cluster $c\in\mathcal{C}=\\{1,2,3,4\\}$, and the
aggregated load profile experienced by the grid under Scenario S1 for $F=2$
and $F=3$. The load shifting flexibility is $\sigma=20\%$.  The attainable
load profile for each cluster $c\in\mathcal{C}$, indicated by the red shaded
areas, are bounded from below and above by
$\underline{d}^{s}_{tc}=\min_{i\in\mathcal{I}}\sum_{f\in\mathcal{F}}{d^{s,(i)}_{tcf}}$
and
$\overline{d}^{s}_{tc}=\max_{i\in\mathcal{I}}\sum_{f\in\mathcal{F}}{d^{s,(i)}_{tcf}}$,
respectively. The aggregated attainable load profile are bounded by
$\underline{d}^{s}_{t}=\min_{i\in\mathcal{I}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}d^{s,(i)}_{tcf}$
and
$\overline{d}^{s}_{t}=\max_{i\in\mathcal{I}}\sum_{c\in\mathcal{C}}\sum_{f\in\mathcal{F}}d^{s,(i)}_{tcf}$.
Note that the iterations associated with the derivations of these bounds, and
similarly the definition of variables $d^{s,(i)}_{tcf}$ follow from the
approach presented in Section IV-C.

Fig. 5: Electricity demand profile under Scenario S1. The red dashed and blue
continuous lines depict the reference case, i.e. $D_{tc}$, and the demand
response to the IBP scheme after load shifting, i.e.
$\sum_{f\in\mathcal{F}}d^{s}_{tcf}$.

From Fig. 5 it can be observed that the utility company strategically selects
the breakpoint values, highlighted in the black dotted lines, to encourage
consumers to monitor their consumption within each pricing period and shift
their flexible load from periods of high to periods of low consumption. These
strategic pricing signals reduce the PAR and lead to decreased utility
operating costs and consumer bills.

###  V-C Impact of shifting cost - $\tau$

This section analyzes the impact of the consumer shifting cost. A total of 156
scenarios are considered for varying price increments
$\xi=\\{0,0.005,\cdots,0.055,0.06\\}$ £/kWh and different values of shifting
costs $\tau=\\{0,0.02,\cdots,0.08,0.1\\}$ £/kWh2, for two and three
consumption-blocks ($|\mathcal{F}|\in\\{2,3\\}$). The load shifting
flexibility is set to $\sigma=30\%$.

Figure 6 depicts the PAR reduction as a percentage of the reference case for
increasing values of the price increment value $\xi$. Each plotted line
illustrates the PAR reduction for different shifting cost coefficients $\tau$.
Similar to Fig. 4, the observed PAR reduction is a unimodal function of $\xi$.
In this case, however, the PAR reduction decreases from 30% to 16% as the
shifting cost coefficient increases from $\tau=0$ to $\tau=0.1$. The latter
follows since the willingness of consumers to shift their demand reduces as
the shifting cost increases, although their shifting flexibility remains the
same.

Fig. 6: PAR reduction as a function of the price increment $\xi$ for different
values of the load shifting cost $\tau$ when $F=2$ (top) and $F=3$ (bottom).
Fig. 7: Electricity demand profile when consumers have zero shifting cost
($\tau=0$ £/kWh2) for cluster 1 (left) and all consumer clusters (right).

An interesting case depicted in Fig. 6 is when $\tau=0$, i.e. when there is no
cost associated with demand shifting. In this context , the same PAR reduction
$(30\%)$ can be obtained for all price increments. This follows since the lack
of shifting cost maximises the amount of shifted demand. To illustrate this
scenario in more detail, Fig. 7 provides the hourly demand of cluster 1 and
the aggregated load profile experienced by the grid. The observed load
profiles can be perceived as unrealistic because the behavior of the consumers
is not reflected by the economic signals embedded in the price structure. This
result highlights the importance of the incorporation of the shifting cost in
(3) to realistically model the demand response of consumers.

##  VI Conclusions

We considered the problem of minimizing the PAR of the power system with
suitable price incentive schemes for consumers. In accordance to this, the
Intraday Block Pricing scheme has been proposed. Under the IBP scheme,
electricity charging is based on intraday pricing periods of equal duration,
where the consumption in each pricing period spans a number of pre-defined
blocks of increasing price. A bilevel optimization problem has been formulated
for the parameter design of the IBP scheme, aiming to minimize the PAR and
simultaneously protect the utility revenue and the consumer bills. A solution
approach has been developed to solve the considered bilevel optimization
problem,  by converting it into an equivalent single-level Mathematical
Program with Equilibrium Constraints, which is then relaxed into a Mixed
Integer Linear Program. In addition, a lower bound to its cost was obtained to
evaluate the conservativeness of the solution approach. The applicability and
implementability of the IBP scheme are demonstrated with numerical simulations
on a realistic setting, that demonstrate s a significant reduction in the PAR,
compared to flat pricing,  as well as simultaneous  benefits for the utility
company and the consumer. Furthermore, the simulation results yielded
equivalent solutions for the original and relaxed problems, demonstrating the
effectiveness of the proposed solution approach.

Future research aims to explore optimal pricing strategies in the context of
social welfare maximization and the implications on the pricing strategies and
profit of the retailer. In addition, we aim to investigate the optimal IBP
price structure design problem to accommodate further penetration of
intermittent renewable generation and plug-in hybrid electric vehicles.

## References

  * [1] P. Pinson, H. Madsen _et al._ , “Benefits and challenges of electrical demand response: A critical review,” _Renew. Sust. Energ. Rev._ , vol. 39, pp. 686–699, 2014.
  * [2] IEA, “Electricity Market Report - December 2020,” Paris, 2020. [Online]. Available: https://www.iea.org/reports/electricity-market-report-december-2020
  * [3] G. Strbac, “Demand side management: benefits and challenges,” _Energy Policy_ , vol. 36, no. 12, pp. 4419–4426, 2008.
  * [4] H. M. Soliman and A. Leon-Garcia, “Game-theoretic demand-side management with storage devices for the future smart grid,” _IEEE Transactions on Smart Grid_ , vol. 5, no. 3, pp. 1475–1485, 2014.
  * [5] A.-H. Mohsenian-Rad, V. W. Wong, J. Jatskevich, R. Schober, and A. Leon-Garcia, “Autonomous demand-side management based on game-theoretic energy consumption scheduling for the future smart grid,” _IEEE Trans. Smart Grid_ , vol. 1, no. 3, pp. 320–331, 2010.
  * [6] J. C. Gómez and M. M. Morcos, “Impact of ev battery chargers on the power quality of distribution systems,” _IEEE transactions on power delivery_ , vol. 18, no. 3, pp. 975–981, 2003.
  * [7] H. Allcott, “Rethinking real-time electricity pricing,” _Resour. Energy Econ._ , vol. 33, no. 4, pp. 820–842, 2011.
  * [8] T. Namerikawa, N. Okubo, R. Sato, Y. Okawa, and M. Ono, “Real-time pricing mechanism for electricity market with built-in incentive for participation,” _IEEE Trans. Smart Grid_ , vol. 6, no. 6, pp. 2714–2724, Nov. 2015.
  * [9] R. Deng, Z. Yang, M. Chow, and J. Chen, “A survey on demand response in smart grids: Mathematical models and approaches,” _IEEE Trans. Ind. Informat._ , vol. 11, no. 3, pp. 570–582, Jun. 2015.
  * [10] M. Behrangrad, “A review of demand side management business models in the electricity market,” _Renew. Sust. Energ. Rev._ , vol. 47, pp. 270–283, Jul. 2015.
  * [11] M. Fahrioglu and F. L. Alvarado, “Designing incentive compatible contracts for effective demand management,” _IEEE Trans. Power Syst._ , vol. 15, no. 4, pp. 1255–1260, Nov. 2000.
  * [12] D. Papadaskalopoulos and G. Strbac, “Nonlinear and randomized pricing for distributed management of flexible loads,” _IEEE Trans. Smart Grid_ , vol. 7, no. 2, pp. 1137–1146, Mar. 2016.
  * [13] B. Lin and X. Liu, “Electricity tariff reform and rebound effect of residential electricity consumption in China,” _Energy_ , vol. 59, pp. 240–247, Sep. 2013.
  * [14] M. Hung and B. Chie, “The long-run performance of increasing-block pricing in Taiwan’s residential electricity sector,” _Energy Policy_ , vol. 109, pp. 782–793, Oct. 2017.
  * [15] B. Zhou, R. Yang, C. Li, Y. Cao, Q. Wang, and J. Liu, “Multiobjective model of time-of-use and stepwise power tariff for residential consumers in regulated power markets,” _IEEE Syst. J._ , vol. 12, no. 3, pp. 2676–2687, Sep. 2018.
  * [16] Z. M. Fadlullah, D. M. Quan, N. Kato, and I. Stojmenovic, “Gtes: An optimized game-theoretic demand-side management scheme for smart grid,” _IEEE Systems journal_ , vol. 8, no. 2, pp. 588–597, 2013.
  * [17] P. Samadi, H. Mohsenian-Rad, V. W. Wong, and R. Schober, “Real-time pricing for demand response based on stochastic approximation,” _IEEE Transactions on Smart Grid_ , vol. 5, no. 2, pp. 789–798, 2014.
  * [18] R. de Sá Ferreira, L. A. Barroso, P. R. Lino, M. M. Carvalho, and P. Valenzuela, “Time-of-use tariff design under uncertainty in price-elasticities of electricity demand: A stochastic optimization approach,” _IEEE Trans. Smart Grid_ , vol. 4, no. 4, pp. 2285–2295, 2013\.
  * [19] A. Mohsenian-Rad and A. Leon-Garcia, “Optimal residential load control with price prediction in real-time electricity pricing environments,” _IEEE Trans. Smart Grid_ , vol. 1, no. 2, pp. 120–133, Sep. 2010.
  * [20] D. Qiu, D. Papadaskalopoulos, Y. Ye, and G. Strbac, “Investigating the effects of demand flexibility on electricity retailers’ business through a tri-level optimisation model,” _IET Gener. Transm. Distrib._ , vol. 14, no. 9, pp. 1739–1750, 2020.
  * [21] C. Feng, Y. Wang, K. Zheng, and Q. Chen, “Smart meter data-driven customizing price design for retailers,” _IEEE Trans. Smart Grid_ , vol. 11, no. 3, pp. 2043–2054, 2019.
  * [22] H. Akhavan-Hejazi and H. Mohsenian-Rad, “Power systems big data analytics: An assessment of paradigm shift barriers and prospects,” _Energy Reports_ , vol. 4, pp. 91–100, 2018.
  * [23] E. Bompard, Y. Ma, R. Napoli, and G. Abrate, “The demand elasticity impacts on the strategic bidding behavior of the electricity producers,” _IEEE Trans. Power Syst._ , vol. 22, no. 1, pp. 188–197, Feb. 2007.
  * [24] N. Al Khafaf, M. Jalili, and P. Sokolowski, “A novel clustering index to find optimal clusters size with application to segmentation of energy consumers,” _IEEE Trans. Industr. Inform._ , vol. 17, no. 1, pp. 346–355, 2021.
  * [25] S. Nojavan, K. Zare, and B. Mohammadi-Ivatloo, “Optimal stochastic energy management of retailer based on selling price determination under smart grid environment in the presence of demand response program,” _Applied energy_ , vol. 187, pp. 449–464, 2017.
  * [26] S. Maharjan, Q. Zhu, Y. Zhang, S. Gjessing, and T. Basar, “Dependable demand response management in the smart grid: A stackelberg game approach,” _IEEE Trans. Smart Grid_ , vol. 4, no. 1, pp. 120–132, 2013.
  * [27] P. Thorsnes, J. Williams, and R. Lawson, “Consumer responses to time varying prices for electricity,” _Energy Policy_ , vol. 49, pp. 552 – 561, 2012, special Section: Fuel Poverty Comes of Age: Commemorating 21 Years of Research and Policy.
  * [28] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge University Press, 2004.
  * [29] Z. Kmetty, “Load profile classification,” European Commission, Tech. Rep., 2016\.
  * [30] D. Qiu, Y. Ye, and D. Papadaskalopoulos, “Exploring the effects of local energy markets on electricity retailers and customers,” _Electr. Power Syst. Res._ , vol. 189, p. 106761, 2020.
  * [31] L. Gurobi Optimization, “Gurobi Optimizer Reference Manual,” 2019. [Online]. Available: http://www.gurobi.com